[lammps-users] NPT-stress/atom-

Hi,

I have some queries regarding the pressure values specified in fix NPT command and the units. Could anybody please comment on these wrt simulation of a solid (~5000 atoms)?

  1. The specified pressure in the command is the internal pressure (calculated from virial + kinetic terms) and should be equal to the external (or ambient) pressure under equilibrium. True or False?

  2. If I have to calculate the lattice constant of Aluminum under normal ambient conditions, can I use the following?

fix 1 all npt 300.0 300.0 0.1 xyz 1.0 1.0 0.1 drag 2.0 ! for 1 bar (or atm) ambient pressure, 300 K (metal units)

run 50000 !for 50 ps

  1. Is the following a proper representation of the physical condition where ambient pressure is 1 bar (or atm) and a tensile stress of 1001 bar is applied in the z-direction?

fix 2 all npt 300.0 300.0 0.1 aniso 1.0 1.0 1.0 1.0 -1000.0 -1000.0 0.1 drag 2.0 !corresponding to applied axial tensile stress of ~1001 atm or 100.1 MPa, while allowing lateral dimensions to change.

run 50000 ! for 50 ps

  1. When I do as in 2 or 3, fluctuations in pressure (which are expected) are orders of magnitude higher than the atmospheric pressure. Even the values averaged for 5000 time steps (5 ps) and 5000 atoms don’t agree with the specified stresses in the 3 orthogonal directions.

fix avg1 all ave/time 1 100 100 v_pr v_xstrs v_ystrs v_zstrs file Average1 ave running

run 5000

  1. Also, as the axial strain (z direction) is increased, the strain in one of the two transverse orthogonal directions (x or y) is much different than another even the model is symmetrically constructed (lattice fcc 3.986 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1).

  2. For calculating stress, I have averaged stress/atom for the z direction over the whole simulation domain and divided by the volume? Should the stress have only the virial components (as perhaps suggested in Cao et al. 2008, PRB, vol 77) or virial + kinetic energy components?

I’d really appreciate advise on these issues.

Thanks

Rohit Rai

Hi Rohit,

I'll take a stab at answering these questions and let Steve and the others
correct me if I write something in error.

1. The specified pressure in the command is the internal pressure (calculated
from virial + kinetic terms) and should be equal to the external (or ambient)
pressure under equilibrium. True or False?

True.

2. If I have to calculate the lattice constant of Aluminum under normal
ambient conditions, can I use the following?
fix 1 all npt 300.0 300.0 0.1 xyz 1.0 1.0 0.1 drag 2.0
! for 1 bar (or atm) ambient pressure, 300 K (metal units)
run 50000
!for 50 ps

That looks correct, although in metal units the pressure will be 1 bar which
is not exactly equal to 1 atm.

3. Is the following a proper representation of the physical condition where
ambient pressure is 1 bar (or atm) and a tensile stress of 1001 bar is applied
in the z-direction?
            fix 2 all npt 300.0 300.0 0.1 aniso 1.0 1.0
1.0 1.0 -1000.0 -1000.0 0.1 drag 2.0 !corresponding to applied
axial tensile stress of ~1001 atm or 100.1 MPa, while allowing lateral
dimensions to change.
run 50000
! for 50 ps

Again, this looks right except that bar != atm

4. When I do as in 2 or 3, fluctuations in pressure (which are expected) are
orders of magnitude higher than the atmospheric pressure. Even the values
averaged for 5000 time steps (5 ps) and 5000 atoms don¹t agree with the
specified stresses in the 3 orthogonal directions.
fix avg1 all ave/time 1 100 100 v_pr v_xstrs v_ystrs
v_zstrs file Average1 ave running
run 5000

Steve should be able to comment more authoritatively about this issue.
Certainly, there should be a way to estimate how large the fluctuations
should get for a given simulation. Whether the magnitude of fluctuation
exceeds the magnitude of your loading will depend on those conditions.

5. Also, as the axial strain (z direction) is increased, the strain in one of
the two transverse orthogonal directions (x or y) is much different than
another even the model is symmetrically constructed (lattice fcc 3.986 orient
x 1 0 0 orient y 0 1 0 orient z 0 0 1).

First, how different is different? Can you quantify how much one varies
compared to the other? Also, if the loading is sufficiently high to generate
material defects, you will see asymmetric stress fields.

6. For calculating stress, I have averaged stress/atom for the z direction
over the whole simulation domain and divided by the volume? Should the stress
have only the virial components (as perhaps suggested in Cao et al. 2008, PRB,
vol 77) or virial + kinetic energy components?

The stress/atom compute includes the kinetic contribution to the stress
tensor. This value you calculate should match LAMMPS' estimate for pressure
in the z-direction, pzz. As to whether the kinetic contribution should be
included, that's a much longer conversation. Check out articles by myself,
Min Zhou, Delph, and others.

Jon Zimmerman

Hi Jon,

Thanks for your answers. I used MD simulations to calculate Young's modulus for Al and SiC (meam potential). In both cases, I am getting values (with 5000 atoms, pbc) that are about half of the values reported elsewhere: ~38 GPa for Al single crystal (~72 GPa reported), and ~300 GPa for SiC single crystal (~700 GPa reported). I am guessing it may have something to do with the large pressure fluctuations. As you can see below, pressure fluctuates from -721 to +1117 bar which is orders of magnitude higher than the applied pressure (1 bar).
Steve, could you please comment on the order of fluctuations and what can I do to bring it down? Drag parameter is already at the higher limit of the recommended values in the documentation. The averages for 5000 iterations are given in the attached file (Average1). Does a different pressure control (or a different choice of units) work better? Are there other probable reasons for getting low values of Young's modulus? I am using virial contributions to the stress/atom compute for calculating Young's modulus (although in these runs kinetic components are included).

Jon: Could you please give me references to the articles on virial and kinetic contributions that you mentioned in your previous email?

For SiC:
Input:
atom_style atomic
units metal
lattice diamond 4.3 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1
region reg1 block 0 5 0 5 0 25
create_box 2 reg1
create_atoms 1 box basis 1 1 basis 2 1 basis 3 1 basis 4 1 basis 5 2 basis 6 2 basis 7 2 basis 8 2

neighbor 2.0 bin
neigh_modify every 20 delay 0 check no

read_restart restart.50000

mass 1 28 #SiC
mass 2 12 #C

pair_style meam
pair_coeff * * library.meam Si C sic.meam Si C

compute peratom all stress/atom #virial
compute p all reduce sum c_peratom[1] c_peratom[2] c_peratom[3]
variable xstrs equal -c_p[1]/vol
variable ystrs equal -c_p[2]/vol
variable zstrs equal -c_p[3]/vol
variable pr equal press

velocity all create 300.0 5812775
fix 1 all npt 300.0 300.0 0.1 xyz 1.0 1.0 0.1 drag 2.0
fix avg1 all ave/time 1 100 100 v_pr v_xstrs v_ystrs v_zstrs file Average1 ave running

Output:

Lx~21.866576 Ly~21.866576 Lz~109.33288

     Step Temp Press xstrs ystrs zstrs
   50000 300 -33.652298 -189.42975 4.6243777 83.848477
   51000 310.04981 -588.01579 -536.28733 -617.02856 -610.73149
   52000 288.06588 841.65598 754.03135 767.35958 1003.577
   53000 287.9016 1117.5542 1191.5007 1015.2914 1145.8706
   54000 303.03016 372.62781 417.2084 342.05846 358.61657
   55000 305.9671 -592.8846 -716.11467 -582.06652 -480.4726
   56000 301.774 -721.82782 -779.00698 -653.12894 -733.34755
   57000 293.94093 253.28756 210.67405 316.63017 232.55846
   58000 298.73428 -82.846809 -101.55199 55.999367 -202.9878
   59000 299.43185 -311.61448 -247.28776 -286.65983 -400.89586
   60000 293.16826 640.25115 748.88023 668.1746 503.69861

Thanks
Rohit

Senior R&D Engineer
Hawthorne &York Intermational
4645 S 35th Street
Phoenix, AZ 85040
Ph. +1 602 275 4633, Ext 210
Fax: +1 602 273 4503

Average1 (2.05 KB)

Steve, could you please comment on the order of fluctuations and what can I do to bring it down? Drag parameter is already at the higher limit of the recommended values in the documentation. The averages for 5000 iterations are given in the attached file (Average1). Does a different pressure control (or a different choice of units) work better? Are there other probable reasons for getting low values of Young's modulus? I am using virial contributions to the stress/atom compute for calculating Young's modulus (although in these runs kinetic components are included).

Pressures fluctuate. Especially for solids. You might try using a bigger
system. If your pressure is equilibrated (not osillating wildly, box volume
is staying constant), then I doubt another barostat will help. Units
make no difference.

Steve