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.

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.

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.

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

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.

References

  • Henry, Norman and Lonsdale, Kathleen, 1965. Eds. International Tables for X-Ray Crystallography Vol. I. Birmingham, England: Kynoch Press.