Fix deform with NPT to model Uniaxial Tension - appears to give unstable result

I have used a fix-deform erate option to model uniaxial tension in an anisotropic pbc sample. I got what looked like correct results, except as I increment to ever increasing strains. the results are slightly off for the target incremental strains. (e.g 10% is actually 9.5%). I therefore tried running the same simulation with fix deform scale, and calculated the incremental scale factor based on an engineering strain target. using an index variable running from 0.01 to 0.1 (one to 10% strain). The runs takes in a file with atom positions, equilbrates the example using an NPT aniso fix, then runs a fix deform for 20,000 steps. the settles for 10,000 steps (in the test run), then samples output for 10,000 steps. Step size is 0.001 ps, and damping factors are those recommended by fix-npt documentation. First odd thing, is initial fix NPT outputs a change in all dimensions from about 144 a in the sample to 188 in the sample. It holds the correct Lx in each subsequent fix, but keeps growing in the other dimensions, gowing beyong 1000 angstroms. It should shrinking in the Ly and Lz directions if the simulation of behaving correctly.

I tried increasing the settling and sampling steps to 100,000 for each fix, but the problem only got worse, and blew up in the end with core dumps. The only difference I made in this run versus erate based run was changing fix deform to scale and adding in the scale calculation.

Here is the input file text and some initial log file results.

#------------------------ INITIALIZATION ----------------------------
# ------------------------ FORCE FIELDS ------------------------------
read_data Pd_in_poslammps.xyz
pair_style eam/alloy
pair_coeff * * PdH_Zhou.eam.alloy Pd H
mass 1 106.42
mass 2 1.00794
# ------------------------ VARIABLES -----------------------------

variable i index 1. 2. 3.
# this does strains from 0 to 3 variable ippress equal "v\_i\*1000\.0" variable ipppress equal "v\_ippress\-2000\.0" variable jppress equal "v\_j\*10000\.0" variable jpppress equal "v\_jppress\-10000\.0" variable ipp1 equal "v\_ipppress/10000\.0" variable jpp1 equal "v\_jpppress/10000\.0" variable time\_step equal 0\.001 variable time\_npt1 equal 20000 variable time\_npt2 equal 100000 variable time\_npt3 equal 100000 variable tdamp equal "v\_time\_step\*100" variable pdamp equal "v\_time\_step\*1000" variable strain equal 0\.0 variable e\_ramp equal "0\.0005" \#this should give 1 strain every 20000 ramp steps at .001 time steps (1 ps = 1000 timesteps)
#variable strain_final equal "v_e_ramp*v_time_npt1*v_time_step*v_i"

timestep ${time_step} # DO NOT CHANGE

# ------------------------ 0 Gpa ---------------------------------
# change_box all triclinic # add this in
# if original equilibrium where aniso

thermo 10000
thermo_style custom step lx ly lz vol press pxx pyy pzz pe etotal temp &
pxy pxz pyz xy xz yz
variable tmp equal "lx"
variable L0 equal \{tmp\} print " Initial Unit Length, L0 = {L0}"

# ------------ Start at 0 GPA
thermo 10000
thermo_style custom step lx ly lz vol press pxx pyy pzz pe etotal temp pxy &
pxz pyz xy xz yz cella cellb cellc cellalpha cellbeta cellgamma
reset_timestep 0
velocity all create 1 12345 mom yes rot no
fix 1 all npt temp \{temp1\} {temp1} \{tdamp\} x 0\. 0\. {pdamp} y 0. 0. \{pdamp\} & z 0\. 0\. {pdamp} xy 0.0 0.0 \{pdamp\} xz 0\. 0\. {pdamp} yz 0. 0. ${pdamp}

run \{time\_npt3\} unfix 1 \#variable eq2 equal "lx\*\(\(4/v\_natom\)^\(1/3\)\)" variable tmp2 equal "lx" variable Ls equal {tmp2}
variable L0 equal \{tmp2\} print "L0 = {L0}"
print "Ls = ${Ls}"

label start_of_loop_1 # starts at 0 GPA and works up
# ------------------------ strain loop ---------------------------------
variable strain_2 equal "(v_eq2-v_LO)/v_L0" #corrected was v_eq1
variable strain equal "0.01*v_i" #this is on that is used
fix 1 all npt temp \{temp1\} {temp1} \{tdamp\} y 0\. 0\. {pdamp} z 0. 0. \{pdamp\} & xy 0\. 0\. {pdamp} xz 0. 0. \{pdamp\} yz 0\. 0\. {pdamp}
print " eramp = \{e\_ramp\}" print "L0 = {L0}"
print "i = \{i\}" print "Ls = {Ls}"
variable g1 equal "v_i/100."
print "g1 = \{g1\}" \#variable g2 equal "1\.\+v\_g1" \#variable tmpLf equal "v\_L0\*v\_g2" \#variable Lf equal {tmpLf}
#print "Lf = \{Lf\}" \#variable scale\_new equal "v\_Lf\-v\_Ls/v\_Ls\+1" \#variable g1 equal "1\.\+\(1\.\+v\_i/100\.\)" variable g1a equal "\(1\.\+v\_g1\)" print "g1a = {g1a}"
variable g1b equal "(1.+v_g1a)"
variable g2 equal "v_L0*v_g1a"
print "g2 = \{g2\}" variable g3 equal "\(v\_g2\-v\_Ls\)" print "g3 = {g3}"

variable scale_new equal "1+v_g3/v_Ls" #added 1 here in front
#variable scale_new equal "(v_L0*(1+(1+v_i))-v_Ls)/v_Ls
print " scalenew = \{scale\_new\}" fix 2 all deform 1 x scale {scale_new} remap v flip yes # scale 1% longer in x

run ${time_npt1}

unfix 1
unfix 2
variable tmp3 equal "lx"
variable Ls equal \{tmp3\} print "Ls : {Ls}"

fix 1 all npt temp \{temp1\} {temp1} \{tdamp\} y 0\. 0\. {pdamp} z 0. 0. \{pdamp\} & xy 0\. 0\. {pdamp} xz 0. 0. \{pdamp\} yz 0\. 0\. {pdamp}

run \{time\_npt2\} unfix 1 fix 4 all npt temp {temp1} \{temp1\} {tdamp} y 0. 0. \{pdamp\} z 0\. 0\. {pdamp} &
xy 0. 0. \{pdamp\} xz 0\. 0\. {pdamp} yz 0. 0. ${pdamp}

thermo 5000
thermo_style custom step lx ly lz vol press pxx pyy pzz pe etotal temp pxy &
pxz pyz xy xz yz
run \{time\_npt3\} write\_restart {temp1}K_+x+${strain}_strain.equil
unfix 4

next i
jump SELF start_of_loop_1

# SIMULATION DONE
print "All done"

******************************************8

Part of the log file is below:

thermo 10000
thermo_style custom step lx ly lz vol press pxx pyy pzz pe etotal temp pxy pxz pyz xy xz yz
variable tmp equal "lx"
variable L0 equal \{tmp\} variable L0 equal 148\.997671378212 print " Initial Unit Length, L0 = {L0}"
Initial Unit Length, L0 = 148.997671378212

# ------------ Start at 0 GPA
thermo 10000
thermo_style custom step lx ly lz vol press pxx pyy pzz pe etotal temp pxy pxz pyz xy xz yz cella cellb cellc cellalpha cellbeta cellgamma
reset_timestep 0
velocity all create 1 12345 mom yes rot no
fix 1 all npt temp \{temp1\} {temp1} \{tdamp\} x 0\. 0\. {pdamp} y 0. 0. \{pdamp\} z 0\. 0\. {pdamp} xy 0.0 0.0 \{pdamp\} xz 0\. 0\. {pdamp} yz 0. 0. \{pdamp\} fix 1 all npt temp 400 {temp1} \{tdamp\} x 0\. 0\. {pdamp} y 0. 0. \{pdamp\} z 0\. 0\. {pdamp} xy 0.0 0.0 \{pdamp\} xz 0\. 0\. {pdamp} yz 0. 0. \{pdamp\} fix 1 all npt temp 400 400 {tdamp} x 0. 0. \{pdamp\} y 0\. 0\. {pdamp} z 0. 0. \{pdamp\} xy 0\.0 0\.0 {pdamp} xz 0. 0. \{pdamp\} yz 0\. 0\. {pdamp}
fix 1 all npt temp 400 400 0.1 x 0. 0. \{pdamp\} y 0\. 0\. {pdamp} z 0. 0. \{pdamp\} xy 0\.0 0\.0 {pdamp} xz 0. 0. \{pdamp\} yz 0\. 0\. {pdamp}
fix 1 all npt temp 400 400 0.1 x 0. 0. 1 y 0. 0. \{pdamp\} z 0\. 0\. {pdamp} xy 0.0 0.0 \{pdamp\} xz 0\. 0\. {pdamp} yz 0. 0. \{pdamp\} fix 1 all npt temp 400 400 0\.1 x 0\. 0\. 1 y 0\. 0\. 1 z 0\. 0\. {pdamp} xy 0.0 0.0 \{pdamp\} xz 0\. 0\. {pdamp} yz 0. 0. \{pdamp\} fix 1 all npt temp 400 400 0\.1 x 0\. 0\. 1 y 0\. 0\. 1 z 0\. 0\. 1 xy 0\.0 0\.0 {pdamp} xz 0. 0. \{pdamp\} yz 0\. 0\. {pdamp}
fix 1 all npt temp 400 400 0.1 x 0. 0. 1 y 0. 0. 1 z 0. 0. 1 xy 0.0 0.0 1 xz 0. 0. \{pdamp\} yz 0\. 0\. {pdamp}
fix 1 all npt temp 400 400 0.1 x 0. 0. 1 y 0. 0. 1 z 0. 0. 1 xy 0.0 0.0 1 xz 0. 0. 1 yz 0. 0. ${pdamp}
fix 1 all npt temp 400 400 0.1 x 0. 0. 1 y 0. 0. 1 z 0. 0. 1 xy 0.0 0.0 1 xz 0. 0. 1 yz 0. 0. 1

run ${time_npt3}
run 100000
Neighbor list info ...
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 5.65
ghost atom cutoff = 5.65
binsize = 2.825, bins = 57 40 38
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair eam/alloy, perpetual
attributes: half, newton on
pair build: half/bin/newton/tri
stencil: half/bin/3d/newton/tri
bin: standard
Per MPI rank memory allocation (min/avg/max) = 6.99 | 7.089 | 7.378 Mbytes
Step Lx Ly Lz Volume Press Pxx Pyy Pzz PotEng TotEng Temp Pxy Pxz Pyz Xy Xz Yz Cella Cellb Cellc CellAlpha CellBeta CellGamma
0 148.99767 104.37835 105.43542 1639745.5 0.090241289 0.08289879 0.093591156 0.09423392 -3.6657508 -2.1657645 1 -9.4078381e-05 9.3704351e-06 0.00092587676 8.1796564 2.6852845 -5.8388412 148.99767 104.69836 105.63111 93.045028 88.543308 85.519148
10000 188.2189 131.62297 132.93218 3293252.8 13.492851 13.565163 13.461007 13.452382 -0.10131911 605.10372 403.4737 -0.017653675 0.0062981554 -0.01199816 10.288597 3.4120096 -7.3947481 188.2189 132.02448 133.18142 93.058671 88.531964 85.530439
20000 238.12935 166.27939 167.91134 6648617.8 6.5735322 6.5966533 6.5634868 6.5604566 -0.088896712 597.0055 398.06656 -0.0059747578 0.0045521216 -0.0041540252 12.973093 4.3442658 -9.3747831 238.12935 166.78471 168.22894 93.069592 88.520256 85.538829
30000 297.60249 207.6176 209.63887 12953065 3.4032543 3.410874 3.4000245 3.3988644 0.081940433 601.9706 401.26275 -0.0025279953 0.0017423412 -0.0016654449 16.176209 5.4594083 -11.729023 297.60249 208.24682 210.03769 93.075664 88.510571 85.544887
40000 369.30367 257.49532 259.9835 24722862 1.77421 1.7768516 1.7731552 1.7726233 0.098950312 598.79299 399.13299 -0.0010299758 0.00066704853 -0.00064199077 20.036718 6.8015035 -14.569706 369.30367 258.27371 260.48024 93.080544 88.503757 85.550557
50000 455.99773 317.82791 320.88073 46504862 0.94670619 0.94763268 0.94634484 0.94614105 0.071176373 600.95331 400.59174 -0.00047015247 0.00031058175 -0.00027849237 24.704475 8.4232432 -18.0037 455.99773 318.78659 321.49577 93.08406 88.498669 85.555388
60000 560.89087 390.85037 394.59057 86503881 0.50767892 0.5080037 0.50754546 0.5074876 0.046350665 599.43947 399.59905 -0.0002076571 0.00013741277 -0.00010013905 30.351305 10.387159 -22.160042 560.89087 392.02705 395.3488 93.086843 88.494472 85.559634
70000 687.79254 479.20871 483.78612 1.5945406e+08 0.27585924 0.27597613 0.27581409 0.27578751 0.026021716 600.40621 400.2571 -9.14631e-05 4.7789192e-05 -4.0887271e-05 37.181541 12.761189 -27.186582 687.79254 480.649 484.71741 93.088755 88.491396 85.563341
80000 841.25202 586.07134 591.66165 2.9170914e+08 0.15062882 0.15067209 0.15061171 0.15060265 0.016037892 599.76658 399.83733 -3.9830966e-05 2.3077042e-05 -1.9098258e-05 45.442869 15.629979 -33.265221 841.25202 587.83048 592.80215 93.090257 88.489149 85.566264
90000 1026.7071 715.22584 722.03669 5.3021138e+08 0.08292307 0.082939879 0.082917076 0.082912256 0.0081777663 600.1373 400.08972 -1.8779122e-05 9.7506989e-06 -7.7379304e-06 55.425899 19.097748 -40.612021 1026.7071 717.37022 723.43005 93.091501 88.487279 85.568762
100000 1250.6668 871.2043 879.49213 9.5828261e+08 0.045865108 0.045870742 0.045863223 0.04586136 0.0045833624 599.93386 399.95649 -8.03656e-06 4.4894031e-06 -3.7459076e-06 67.481862 23.283703 -49.484427 1250.6668 873.81391 881.19081 93.092496 88.485898 85.57082
Loop time of 565.737 on 10 procs for 100000 steps with 109891 atoms

Performance: 15272.117 tau/day, 176.761 timesteps/s
100.0% CPU use with 10 MPI tasks x no OpenMP threads

MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total

Hi! Have you used this approach successfully with other material systems? Was your system well-equilibrated at 0 bar (I’m assuming you’re using metal units because I don’t see a unit command anywhere) prior to this simulation? If it was equilibrated at any pressure higher than 0 bar (or simulated with NVT) prior to this uniaxial tension simulation, then the expansion you’re seeing is due to the zero pressure. Are you studying loading and unloading behavior?

Will