Dependence of MSD on super cell shape

LAMMPS LAMMPS General Discussion
1

Dear LAMMPS developers and users

I have been bugged for some time with an issue, which I have not found a satisfactory solution to thus far. I am basically using LAMMPS to generate representative configurations of atomic displacements of a crystal for the purpose of including the effect of atomic motion into calculations of electron scattering for the purpose of transmission electron microscopy and electron energy loss spectroscopy. The main quantity we are interested in from the MD perspective is the MSD of each atom in the super cell, which governs the strength of incoherent scattering.

What we usually do is that we set up an orthogonal super cell of the shape of about 3-5 nm in lateral dimensions (x and y) and then match the experimental thickness of the specimen in the last dimension (z), typically several tens of nm. We apply periodic boundary conditions in all directions, since surface effects play a minor role. We apply a NVT or langevin thermostat and sample atomic displacements at a temperature of 300 K.

This set up of the model system has worked out nicely for most of our calculations, but I have recently started to go to more elongated super cells (about a 1:10 to 1:20 ratio of lateral to z dimensions) in a calculation for diamond. Here I have noted that the MSD shows a behavior, which is not expected for a crystal at thermal equilibrium. The MSD does not fluctuate around its mean (like in the attached example for a 10x10x10 super cell), but it seems that the mean itself is time-dependent to some extent along all three directions (see the other attached plot for 10x10x170). I have done a more detailed analysis and the MSD of the atoms also shows a strong dependence on where the atom is located along the z coordinate in the super cell. Furthermore the deviations are mostly connected with lower frequencies below about 8 THz. So I suspect that it is primarily acoustic modes in our model system, which are not representative of the acoustic vibrational modes in crystalline diamond.

Through testing different variables (initial conditions, length of the trajectory, thermostat, force field, boundary conditions) I have developed the understanding that the origin of the problems is simply the unequal side lengths of the super cell, since the ratio of lateral size to z length seems to be the only parameter, which meaningfully affects the results. For a ratio of 1:1 everything looks as expected with equal MSD along all directions, while I observe gradually larger and larger deviations from the expected MSD for periods of time for smaller and smaller ratios. I think this observation has to do with the unequal implicit grid of allowed wave vectors \Delta \mathbf{k} = (\Delta k_x, \Delta k_y, \Delta k_z), \Delta k_x = \Delta k_y \neq \Delta k_z, for vibrational modes in the super cell, but I don’t have full intuition how this comes about. I have tried to find information on this topic in the literature, but to no avail. It seems that most MD simulations of solids work with equal side lengths of the super cell, so I have not found any mentioning of the issue.

I would appreciate to hear your thoughts and comments on this issue and if you have encountered something similar. I am also grateful for any hints on how one can deal with this issue and/or tips for additional reading. The most obvious solution is to me that we need to rethink our calculation strategy and always use super cells which have a ratio of lateral to z dimensions of about 1:1, but this brings other problems which I would rather avoid, so I think a good understanding of the origin of our issue is necessary to start with.

Thanks for reading this huge wall of text and I look forward to any replies. Have a nice day!

Kind Regards,
Paul

direction projected MSD as a function of time step for super cell size 10x10x10
image

direction projected MSD as a function of time step for super cell size 10x10x170
image

LAMMPS input file lmp.in

units            metal
boundary         p p p
atom_style       atomic

timestep         0.0005


lattice          diamond 3.57329 
region           simbox block 0 10 0 10 0 170
create_box       1 simbox
create_atoms     1 box

mass             1 12.011

group            c type 1

pair_style       tersoff
pair_coeff       * * C.tersoff C

thermo           1
thermo_style     custom step temp ke pe econserve press pxx pyy pzz

# Set equilibration time, measurement time and other variables
variable         tequi equal 10000
variable         tmeas equal 400000
variable         ttot  equal ${tequi}+${tmeas}
variable         taurelax equal "dt*100"
variable         temp equal 300
variable         inittemp equal 2*${temp}
variable         random_seed equal 983632

velocity         all create ${inittemp} ${random_seed} loop local rot yes dist gaussian mom yes
fix              mynve all nve
fix              mylang all langevin ${temp} ${temp} ${taurelax} 4643963 zero yes
fix              mypos all recenter INIT INIT INIT
fix              mymom all momentum 1 linear 1 1 1 angular rescale

compute          mymsd1 c msd com yes average yes
fix              myave1 c ave/time 1 1 1 c_mymsd1 file msd1_c.dat mode vector

run              ${tequi}

compute          mymsd2 c msd com yes average yes
fix              myave2 c ave/time 1 1 1 c_mymsd2 file msd2_c.dat mode vector

run              ${tmeas}

Potential file C.tersoff created from the Tersoff potential for carbon included in LAMMPS

# DATE: 2011-04-26 UNITS: metal CONTRIBUTOR: Aidan Thompson, [email protected] CITATION: Tersoff, Phys Rev B, 39, 5566-5568 (1989)

# Si and C mixture, parameterized for Tersoff potential
# this file is from Rutuparna.Narulkar @ okstate.edu
# values are from Phys Rev B, 39, 5566-5568 (1989)
# and errata (PRB 41, 3248)

# Tersoff parameters for various elements and mixtures
# multiple entries can be added to this file, LAMMPS reads the ones it needs
# these entries are in LAMMPS "metal" units:
#   A,B = eV; lambda1,lambda2,lambda3 = 1/Angstroms; R,D = Angstroms
#   other quantities are unitless

# format of a single entry (one or more lines):
#   element 1, element 2, element 3,
#               m, gamma, lambda3, c, d, costheta0, n,
#               beta, lambda2, B, R, D, lambda1, A

# Tersoff's nomenclature -      -       -       c         d        h          n           beta            mu        B       R = (R+S)/2     D = (S-R)/2   Lambda      A

C       C       C       3       1       0       38049.0   4.3484  −0.7482     0.72751     1.5724E-07      2.2119  393.747   1.95            0.15          3.4879   1393.6
2

Here you are basically simulating an infinite crystal using a very anisotropic unit cell. For this task, I would rather use an isotropic unit cell. To me, the only explanation for the observed difference are numerical errors. Could you share more details about your calculations? These are essential to pin-point the possible origin of the observed behaviour of the MSD.
PS I believe the first plot refers to the 10x10x10 supercell.

3

This is an intuitive (and nice!) result if you consider the Debye model of a solid’s heat capacity. If you simulate a crystal in a periodic box of length L, its phonon spectrum is bounded to a maximum wavelength of L (and a matching minimum frequency). Thus the heat capacity and corresponding thermal diffusivity are also artificially lowered.

Increasing the box size increases the phonon spectrum available, and thus the heat capacity and thermal diffusivity. I suspect the particular effect of anisotropy is not large unless the crystal lattice itself is significantly anisotropic. Certainly your data shows the x- and y-diffusivity increasing with z-length, indicating that phonons thermalize well between those different dimensions.

This phenomenon is also known for liquids in MD (through a hydrodynamic model for displacements rather than a phonon model) as a “finite size correction to diffusivity”, and you will sometimes see the named “Yeh and Hummer correction”. If you search for the “finite size” keyword you may find a good description of the solid-state case.

If not – it seems like you have a very nice test case with which you can study this effect! I would try to establish a relationship first between the box size and heat capacity – for which the Debye model should give you sensible results – and then between the heat capacity and the MSD, which would be much harder but I think a basic empirical relationship would be something like the MSD scaling as the square root of the heat capacity. Then you can extrapolate between calculations at different box sizes using these relationships (see the liquid literature for similar calculations).

1 Like
4

Thank you very much for your answer, @hothello. Yes, I am trying to simulate an infinite crystal with an anisotropic super cell. Though I observe the same behavior when I remove the periodic boundary along z. I also think that using an isotropic super cell is the cleanest solution, but then we run into some troubles with our subsequent electron scattering calculations, since we thus far need to match the experimental thickness with the box size. So I am still curious if there is anything we can do to improve the simulation.

Thanks also for pointing out my copy-paste mistake regarding the super cell size of the first image. i have updated my answer and also included my LAMMPS input file and potential for reproducibility. Since I am a new user, I cannot upload files, so I just pasted the file contents.

5

Thank you very much for your extensive answer, @srtee! I see what you are hinting at. However I think I have not stated our goal with these simulations clear enough, since we are not directly interested in computing the MSD or the heat capacity. We are rather interested in generating some sensible “snapshots” of the atomic structure using MD. These structure snapshots are subsequently used as an input fur further electron scattering calculations. For us the MSD is therefore more of a “benchmark tool” to see if everything behaves as we expect it to in the MD simulation. And we essentially want the MSD to be equal in all three dimensions for all atoms without any larger (time-dependent) deviations, since in an infinite crystal, this would not be the case. The catch is that for now we have always needed to match the experimental thickness of the sample we want to model with the z length of the super cell.

So in short I am mainly wondering if this idea is fundamentally flawed due to the anisotropic super cell we are using (I start to think that it is, but it is not yet clear why exactly?) or if there is something to understand to make atoms move more isotropically for our purposes.

6

Why not directly examine the radial distribution function (RDF)? It would be much more relevant to a subsequent scattering study, I would have thought.

I don’t have an intuitive sense of how your variables affect the RDF but my hunch is you will find they are very similar for the different box aspect ratios you are trying.