Quasi-static stretching at 0K

Dear all,

I want to simulate quasi-static stretch along y-direction at 0K for alpha-quartz. The structure is relaxed after cooling from room temperature. After that, I want to apply quasi-static stretching along y-direction only such that the y-side is increased by 0.00001 each step. I hope to loop this for 10 times and after each time, the structure should be relaxed. I found from online that we can put minimize command in a loop with stretching the structure incrementally. However, I am not very clear on how to fulfil this.

I attach my script in the end where I output intermediate stress&strain in a file. There are some problems in the file:

Fix print output for fix def1

1.00000000011046e-05 -0.0817324513709522 -0.36564318648842 -0.0960313824874549

1.00000000011046e-05 0.174520072522468 -0.1658028766301 0.226471384221925

1.00000000011046e-05 -0.564709020717339 -0.566876156108682 -0.4047267590484

1.00000000011046e-05 1.33305850224231 0.383342237194663 1.07880524464239

1.00000000011046e-05 -3.26166178333311 -2.10974464192211 -2.82491780546663

1.00000000011046e-05 7.11027748421183 3.47216717139645 5.94526335103269

1.00000000011046e-05 -14.7009575476715 -8.22056269333739 -12.3104819875476

1.00000000011046e-05 19.805079249192 10.1792292234112 16.9360698950953

1.00000000011046e-05 -38.8945644375914 -21.8577323625349 -32.7842146897005

1.00000000011046e-05 35.6893543000372 18.4028751919197 31.1743306861121

1st column is strain which should by 0.00001 each time; 2-4th columns are pressures along x, y, z direction which seems not be relaxed.

I have checked that, until unfix 2, the structure is indeed relaxed. The question is how to simulate quasi-static stretching then.

Many thanks

Bingyu Cui

#Initialization############################################################

units metal

dimension 3

boundary p p p

atom_style charge

Atom definition##########################################################

lattice custom 5.4047 &

a1 0.9092 0.0000 0.0000 &

a2 -0.4546 0.7873 0.0000 &

a3 0.0000 0.0000 1.0000 &

basis 0.4697 0.0000 0.0000 &

basis 0.0000 0.4697 0.6667 &

basis 0.5303 0.5303 0.3333 &

basis 0.4133 0.2672 0.1188 &

basis 0.2672 0.4133 0.5479 &

basis 0.7328 0.1461 0.7855 &

basis 0.5867 0.8539 0.2145 &

basis 0.8539 0.5867 0.4521 &

basis 0.1461 0.7328 0.8812

region simbox block 0.0 7.0 0.0 7.0 0.0 7.0 units lattice

create_box 2 simbox

create_atoms 1 box basis 1 1 basis 2 1 basis 3 1 basis 4 2 basis 5 2 basis 6 2 basis 7 2 basis 8 2 basis 9 2

mass 1 28.0855

mass 2 15.9994

group siliconatoms type 1

group oxygenatoms type 2

set group siliconatoms charge 2.4

set group oxygenatoms charge -1.2

Atoms interactions settings##################################

Si type 1, O type 2

pair_style table linear 64644 pppm

pair_coeff 1 1 potential.Ewald SiSi

pair_coeff 1 2 potential.Ewald SiO

pair_coeff 2 2 potential.Ewald OO

kspace_style pppm 1.0e-6

kspace_modify gewald 0.31465093412

kspace_modify kmax/ewald 23.19 13.39 17.01

neighbor 2.0 bin

neigh_modify every 1 delay 0 check yes

pair_write 1 2 10001 r 0.51 10.17 tableEW.txt Si_O 2.4 -1.2

pair_write 2 2 10001 r 0.51 10.17 tableEW.txt O_O -1.2 -1.2

pair_write 1 1 10001 r 0.51 10.17 tableEW.txt Si_Si 2.4 2.4

Equilibrate

velocity all create 298.0 635845 dist gaussian

fix 0 all npt temp 298.0 0.1 0.1 iso 1.01325 0.0 1.0

thermo_style custom step temp pe etotal press vol density lx ly lz pxx pyy pzz pxy pxz pyz fmax

thermo 10000

dump 1 all atom 10000 dump.qua

run 80000

fix 2 all box/relax x 0.0 y 0.0 z 0.0 vmax 0.005

min_style cg

minimize 1e-25 1e-25 5000 10000

unfix 0

unfix 2

reset_timestep 0

variable tmp equal “ly”

variable L0 equal ${tmp}

print “Initial Length, L0: ${L0}”

fix 3 all npt temp 0.1 0.1 0.1 x 0.0 0.0 1.0 z 0.0 0.0 1.0

min_style cg

minimize 1e-25 1e-25 5000 10000

change_box all y scale 1.00001 remap

variable strain equal “(ly-v_L0)/v_L0”

variable p1 equal “v_strain”

variable p2 equal “-pxx/10000”

variable p3 equal “-pyy/10000”

variable p4 equal “-pzz/10000”

fix def1 all print 1 “{p1} {p2} {p3} {p4}” file Al_SC_100.def1.txt screen no

thermo 1

thermo_style custom step v_strain temp v_p2 v_p3 v_p4 ke pe press

run 10

I don’t see any looping in your code, which is what you said you
were asking about. I suggest you don’t start with 100 line input
script. Start with a very simple script, and put a loop around
a “minimize, run, change_box” sequence of 3 commands.
Visualize the results, look at the thermo output, including box
size. When you know how to do it, then add back the
complexity of your 100 line script, one feature at a time.

Steve