4.21. DYNAMO: Molecular Dynamics structure factor calculation
Authors: Doug du Boulay
Contact: D. du Boulay, Materials and Structures Laboratory, Tokyo Institute of Technology, Nagatsuta, Midori Ku.
DYNAMO calculates structure factors from molecular dynamics data such as that from the program MXDORTO (K. Kawamura). In MD simulations atomic nuclei positions are recorded at pico second intervals, typically for the order of nanosecond duration simulation experiments. These data sets can be used to calculate the average structure factors that one might expect to see in equivalent diffraction experiment averaged over the timescales involved and also over the spatial extent of an entire crystal. The "experimental" MD structure factors can then be compared with structure factors calculated in the standard crystallographic manner, incorporating temperature factors to model the time and spatially averaged electron density.
4.21.1. Process
Structural crystallography is ordinarily focussed on the most minimal set of independent atom sites arranged in the assymmetric unit or, the arrangement of a symmetry extended subset of sites contained within one unit cell. Consequently crystallographic software like Xtal has generally been optimised to work with the most primitive set of unique atomic information. In contrast, molecular dynamics (MD)simulation studies dictate an approach wherein all atoms in the entire crystal are assumed to evolve fully independently, over time scales of the order of femto to pico seconds. Practical limitations involving storage space and evaluation times place severe restrictions on this ideal so generally it is sufficient to adopt "very large" unit cells with no internal symmetry (i.e. P1) but subject to periodic boundary conditions.
To compare MD simulation data with charge density experiments we need to convert the MD atomic site information into structure factors and thereby into charge densities. DYNAMO can assist with the first step and CRYLSQ in conjunction with FOURR enables the second.
The results of MD simulations are generally large data files consisting of the instantaneous positions of atomic nuclei for each step. At each instant an MD atom effective scattering contribution should thus be very well accounted for by the elemental form factors of Cromer and Mann (19??) determined for isotropic noninteracting and non-vibrating atoms. By summing the form factors of all the independent MD atoms across the entire MD cell, and by further averaging of the MD cells across repeated MD time steps very well averaged, effective scattering factors can be determined. This is a convenient way of converting very extensive MD data sets into a relatively simple representation in the form of time and spatially averaged structure factors, exactly the same kind of data measured experimentally in single crystal X-ray scattering experiments.
If the MD effective structure factors Fmd are treated as experimentally determined structure factors, then an effective scattering model can be refined to model the MD data using standard least squares programs, such as CRYLSQ. In the refinement however there is absolutely no need to refine the scale factor because the MD structure factor data have all been calculated on an absolute scale. There is also no extinction and no dispersion because the MD data do not incorporate such effects, by default. It should therefore, be sufficient to refine only the first and second moments of the density, i.e. the mean atomic positions and anisotropic temperature factors. Any residual features in the charge density after Fourier transforming the structure factor differences therefore reflect real deviations from the averaged, harmonic model and, not chemical or bonding redistributions. An anharmonic structure factor model may or may not account for such features, but in any case can not be undertaken within Xtal.
4.21.2. Calculation Performed
For extensive MD data sets the calculations performed by DYNAM can be rather time consuming. As a small measure of optimization, all reflections on the bdf are loaded into memory simultaneously. Unfortnately this puts limits on the number of reflections for which MD structure factors can be evaluated. As each atom site in each requested frame is read, its contribution to the effective MD cell scattering power is accumulated for every reflection stored in memory. After each frame a running average is performed on the MD structure factors. Finally the MD structure factors are normalised by the determinant of the cell transformation matrix to put them on the same scale as the archive BDF sites. Subsequently, if atom sites exist on the archive, then normal Fcal structure factors are calculated, and relative scales compared (though not applied to archive). In principle, if your archive model sites and molecular dyanmics sites are consistent then te scale factor should be 1.0.
The Fmd values can be stored as iether Frel, or Fcal, in case it is desired to compare the Fmd values with some experimental structure factors already existant as Frel on the archive. In that case it may be necessary to apply dispersion and extinction corrections to the Fmd values to match the experimental Frel terms.
4.21.3. File Assignments
Reads atom and reflection data from the input archive bdf
Reads MD data file
Writes calculated structure factors to the output archive bdf
4.21.4. Example
STARTX : spinel cell 8.2228 8.2228 8.2228 90 90 90 cellsd 0.0001 0.0001 0.0001 0 0 0 spaceg F_d_-3_m:2 celmol 8 celcon Li 8 celcon Mn 16 celcon O 32 exper *2 0.60470 ADDATM UOV .03500 PARENT ATOM Li1 .12500 .12500 .12500 $1 1.000 .00000 .00000 .00000 U Li1 .01295 .00095 ATOM Mn .50000 .50000 .50000 $1 1.000 .00000 .00000 .00000 UIJ Mn .00521 .00521 .00521 -.00083 -.00083 -.00083 SUIJ Mn .00005 .00005 .00005 .00004 .00004 .00004 ATOM O .26526 .26526 .26526 $1 1.000 .00008 .00008 .00008 UIJ O .01921 .01921 .01921 -.00461 -.00461 -.00461 SUIJ O .00021 .00021 .00021 .00021 .00021 .00021 ADDREF ffac : use interpolated form factors in DYNAM/FC limits *4 1.0 : *9 yes hklgen hkl frel sigf : generate reflection hkl values for DYNAM end SORTRF frel aver 1 : remove symmetry redundant reflections end : Be wary of rounding errors in accumulating averages over many : frames. <F(n-1)>(n-1)/n + F(n)/n Doesn't scale well. DYNAMO frel list file fort.81 : save Fmd as Frel, list all refs , read file frames 1 100 : start from frame 1 read 100 frames transf 4 0 0 0 4 0 0 0 4 : change cell from 8x8x8 -> 32x32x32 typinf O 2048 Mn 512 Mn 512 Li 512 : order matches fort.81 order end : improve the fcal model by refining aniso Uij CRYLSQ cy 20 wu fr p1 tl 0. l1 : no dispersion/extinction needed noref skf FOURR fdif grid 96 96 96 |
In this example reflections in an asymetric unit of reciprocal space are generated for the 8x8x8 cell entered in STARTX. In DYNAM those hkl values are transformed to the 32x32x32 unit cell setting and MD structure factors are calculated before being transformed back to the original 8x8x8 cell. This is followed by a normal structure factor calculation within DYNAM, for magnitude comparison purposes.
:STARTX and ADDATM as for example 1 :followed by ... EXPAND cell 3 : generate 5x5x5 cells worth of redundant sites STARTX upd : throw away all symmetry sgname P1 NEWCELL : transfrom archive cell from 8x8x8 to 32x32x32 transf 4 0 0 0 4 0 0 0 4 STARTX upd : override wierd centred symmetry sgname P1 EXPAND not : contract to 32x32x32 cell P1 independents ADDATM upd : update atomic constraints ADDREF ffac : use interpolated form factors in DYNAM/FC limits *4 0.7 : *9 yes : too many reflections for DYNAM to handle hklgen hkl frel sigf : generate empty hkl packets end SORTRF frel aver 1 : remove redundant reflection indices DYNAMO frel list file fort.81 : save Fmd as Frel, list all refs , read file frames 1 10 : start from frame 1 read 10 frames transf 1 0 0 0 1 0 0 0 1 : No change from MD cell to archive bdf typinf O 2048 Mn 512 Mn 512 Li 512 : order matches fort.81 order FOURR fdif pset 1 : use MD phases via pset 1 flag grid 192 192 192 |
Here structure factors are computed for the 32x32x32 P1 cell. There is no need for Friedel pairs as there isn't any dispersion. To correctly do the Fourier transformation, the MD phases are stored on the archive (in addition to the normal structure factor phases - which are nearly all zero owing to the imposed symmetry of the original 8x8x8 model). The MD phases are accessed in FOURR via the pset 1 flag. Alternatively, the MD structure factors (now Frel values) and the whole unit cell could then be transformed back to the 8x8x8 F_d_-3_m:2 cell for structure factor averaging and model comparisons.