NPT ensemble dynamic relaxation process, temperature rise problem

LAMMPS LAMMPS General Discussion
1

Dear LAMMPS users,
Hello,
Excuse me: Why does the temperature keep rising during my dynamic relaxation of the NPT ensemble at 10 K? The temperature keeps increasing from 10 K with the number of running steps. What is causing this? How to solve this difficult problem?
The following is my script:

---------------------------------- Initialize Simulation -----------------------------

atom_style full

units metal

boundary s s p

atom_modify map array sort 0 0.0

------------------ -------------------- Create Atoms -----------------------------------

read_data data.lmp

------------------------------------ Create Regions -----------------------------------

region bto_fixed1 block 0.0 80.32 0 32.128 0 32.128 units box # 0 20 0 8 0 8 units lattice #

region bto_middle block 80.32 240.96 0 32.128 0 32.128 units box # 20 60 0 8 0 8 units lattice #

region bto_fixed2 block 240.96 257.024 0 32.128 0 32.128 units box # 60 64 0 8 0 8 units lattice #

region bto_mobile union 2 bto_middle bto_fixed2

region bto_total union 3 bto_middle bto_fixed1 bto_fixed2

group Gr_bto_fixed1 region bto_fixed1

group Gr_bto_middle region bto_middle

group Gr_bto_fixed2 region bto_fixed2

group Gr_bto_mobile region bto_mobile

group Gr_bto_total region bto_total

#------------------------------------ potential ------------------------------

pair_style born/coul/wolf 0.25 12.0 12.0

pair_coeff * * 0.0 1.0 0.0 0.0 0.0

pair_coeff 1 3 1061.3 0.374 0.0 0.0 0.0

pair_coeff 2 3 3769.93 0.2589 0.0 0.0 0.0

pair_coeff 3 3 4740.0 0.2686 0.0 160.0 0.0

#------------------------------------ neighbor -------------------------------

neighbor 2.0 bin

neigh_modify every 1 delay 0 check yes

#---------------------------------- Equilibration ----------------------------

velocity Gr_bto_total create 10 12345 mom yes rot no

#--------------------------Thermo information shows--------------------------------------

thermo_style custom step temp etotal ke pe cella cellb cellc cellalpha cellbeta cellgamma vol press

#---------------------------------------Display thermo-------------------------

timestep 1.0e-5

thermo 1

fix 1 Gr_bto_total npt temp 10 10 0.001 z 0 0 0.001

dump 1 all custom 10 dumps/stab1_*.txt id type x y z q

run 15000

unfix 1

undump 1

2

When running a simulation with a thermostat yet the temperature increases (significantly) then this is a sign of a problem with your model or simulation settings.
It is not really a difficult problem but rather a fundamental problem, i.e. you have to have understood certain fundamental aspects about MD simulations. This is not the place to give you an in-depth discussion. For that you need to discuss with your adviser/supervisor/tutor.

The short version is, usually this comes down to either a too large timestep or bad potential parameters. In your input you are using a timestep of 0.01 femtoseconds, which is extremely small. The maximum timestep possible depends on the mass of the particles and the steepness/stiffness of the potentials. For a typical atomic/molecular system, a timestep of up to 0.5 femtoseconds should always be possible (i.e. when you have hydrogen atoms). That number scales with sqrt(mass), so for a system with only oxygen atoms your possible timestep would be 2fs (=sqrt(16)*0.5).

So why do you choose such a miniscule timestep? and what happens if you use a larger one?

The fact that you still get a continuing increase in temperature despite using a thermostat under these conditions suggests that there is a problem with your force field parameters that lets atoms get too close resulting in too large forces and then to fast motions which will destabilize the time integration and thus lead to divergence and crashes. The increase in temperature is a “precursor” to that.

The matter of the fact is that you should be able to have a conserved total energy without a thermostat, i.e. using fix nve. This is a necessary requirement for a meaningful MD simulation. For more on that, please read up on the topic in your MD simulation text book.

3

Dear Prof.akohlmey,
Thank you for your reply. Your reply is very detailed. After reading your explanation, I checked my simulation script. As you said, my problem is that the time step is too small. With your help, I have Solved my problem.
Jackie

4

You are not making sense. If raising the time step relaxes the system into the desired temperature, then your original description was incorrect and the temperature would eventually relax as well; only your simulation was too short to determine it reliably.

5

After I changed the time step, I used the NPT ensemble to relax at 10 K. During the relaxation process, my temperature rose from 10K to 350K in 1000 steps, and then decreased from 350K to about 10K. After 1000 steps, the temperature was relaxed. The temperature in Henan is around 10K.

6

You should see the exact same behavior with a shorter timestep, only you need to run for (many) more steps. That is why I asked about why you chose the shorter timestep. Since you changed it from the default value, you must have changed it for a reason.

7

Yes, I observed. After 1000 steps-1500000 steps, the temperature has been around 10K.