All,
I am trying to equilibrate an Ice 1h crystalline structure to 250K in LAMMPS. When I tried to do this at 250K directly, the structure lost it's crystalline shape and became disordered post equilibration.
Based on the advise given by previous researchers who have worked with ice structures in LAMMPS, this should be done in small incremental steps from 10, 100, 150, 200, to finally 250K using a 0.5fs timestep. A minimum of 60ps relaxation is recommended at each temp.
Can anyone kindly help me understand what a script for this part of the simulation will look like? Below I post, the first few steps of the equilibration using fix npt (chosen after careful consideration of other ensembles and trials with fix nvt) as per my understanding. Please let me know if this looks correct. Any insight is appreciated.
#First step equilibration at 50K
timestep 0.5
velocity all create 50 34455 dist gaussian mom yes rot yes
fix 1 all shake 1e-6 10 1000 b 1 a 1
fix 2 all npt temp 50 50 $(100.0*dt) iso 0 0 $(1000.0*dt)
thermo_style custom step temp pe ke etotal press vol density
thermo 1000
run 120000
unfix 2
#Second step equilibration to 100K
fix 2 all npt temp 50 100 $(100.0*dt) iso 0 0 $(1000.0*dt)
thermo_style custom step temp pe ke etotal press vol density
thermo 1000
run 5000
unfix 2
#Third step equilibration at 100K
fix 2 all npt temp 100 100 $(100.0*dt) iso 0 0 $(1000.0*dt)
thermo_style custom step temp pe ke etotal press vol density
thermo 1000
run 120000
unfix 2
J. Javier Burgos-Mármol
Thank you so much for the reply. I will check out the threads.Could you advise on how to determine the required number of steps or how much dT/dt is acceptable?
As I understand it, the rates at which am doing the simulation are the main reason the structure is not holding it's crystalline structure.
Good afternoon. Will you maybe use loop run for every step equilibration?
What if you add minimization energy after each run?
Good luck and best wishes from me.
Example script(party).
variable myTemp equal temp
label loop
variable a loop 21
run 1000
print "Temperature is ${myTemp}"
if "${myTemp} > 100.0" then "jump SELF break"
next a
jump SELF loop
label break
print "ALL DONE, Final temperature is ${myTemp}"
- Badges
- Science topic
- Similar topics
- Physical Sciences
- Materials Science