Strange phenomena during increasing of temperature

Dear all,

I am reproducing thermal properties of strontium titanate from a literature. My input script arranged the system (3645 atoms) according to the experimental lattice parameters and did thermalization/equilibration in NPT for 20000 steps. After that, the temperature was increased from 298K to 2000K with NPT thermostat (fix 1 all npt temp 298 2000 100.0 iso 1.01325 1.01325 1000.0). During the increasing of temperature and at temperature around 1800K, the temperature suddenly shoot up drastically and then immediately shrink to zero (and kinetic energy as well) and the volume collapsed to the volume at room temperature. I used VMD to visualize the atoms for the whole process and found out the system expanded due to temperature increased, which I expected, but it suddenly collapsed and all atoms were not moving. The melting point is 2080K, but the author manage to simulate up to 2200K with other (Japanese) MD code. I wonder why would this strange phenomena occured?

Thank you for your time.

I attached my input script for your perusal:

strontium titanate system

atom_style charge
units real
dimension 3
boundary p p p

lattice

lattice custom 3.9051 a1 1.0 0.0 0.0 a2 0.0 1.0 0.0 a3 0.0 0.0 1.0 basis 0.0 0.0 0.0 basis 0.5 0.5 0.5 &
basis 0.5 0.5 0.0 basis 0.5 0.0 0.5 basis 0.0 0.5 0.5
region simbox block 0 9 0 9 0 9
create_box 3 simbox
create_atoms 1 box basis 1 1 basis 2 2 basis 3 3 basis 4 3 basis 5 3
group sr type 1
group ti type 2
group o type 3
mass 1 87.62
mass 2 47.867
mass 3 15.999
set group sr charge 1.2
set group ti charge 2.4
set group o charge -1.2

force-field

pair_style hybrid/overlay born 11.0 coul/long 11.0 morse 11.0
pair_coeff * * coul/long
pair_coeff 1 1 born 0.32 0.32 2.396 100.0 0.0
pair_coeff 1 2 born 0.34 0.34 2.253 250.0 0.0
pair_coeff 1 3 born 0.32 0.32 3.124 200.0 0.0
pair_coeff 2 2 born 0.36 0.36 2.110 625.0 0.0
pair_coeff 2 3 born 0.34 0.34 2.981 500.0 0.0
pair_coeff 3 3 born 0.32 0.32 3.852 400.0 0.0
pair_coeff 1 3 morse 2.41 1.18 2.7615
pair_coeff 2 3 morse 4.23 3.82 2.1923

kspace_style ewald 0.0001
pair_modify tail yes shift no

initialize

timestep 1.0
velocity all create 298.0 1234567 mom yes rot yes dist gaussian
neighbor 4.0 bin
neigh_modify every 1 delay 0 check yes

output

variable N equal step
variable pote equal pe
variable Etotal equal etotal
variable T equal temp
variable Press equal press
variable V equal vol
variable kine equal ke
variable Lx equal lx/9
variable Ly equal ly/9
variable Lz equal lz/9
fix extra all print 1 “{N} {T} {V} {Press} {kine} {pote} {Etotal}" file strontium_titanate.out fix cell all print 1 "{N} {T} {Lx} {Ly} {Lz}” file strontium_titanate.cell
dump snapshot all custom 1000 strontium_titanate.lammpstrj id q type x y z

equilibration & thermalization

fix NPT all npt temp 298.0 298.0 100.0 iso 1.01325 1.01325 1000.0
run 20000

fix NPT all npt temp 298.0 2000.0 100.0 iso 1.01325 1.01325 1000.0
run 340000

fix NPT all npt temp 2000.0 2000.0 100.0 iso 1.01325 1.01325 1000.0
run 220000

Regards,
Christopher

I don't know. I would monitor thermo output at a fine time scale near
when the transition happens and see if anything goes haywire. Ditto with
dumps.

Steve

Dear Dr. Steve and all,

Thanks and I’ll monitor the atoms and thermo output carefully near the weird transition. I suspect the strange phenomena may be due to Nose-Hoover thermostat’s problem. Therefore I rerun the simulation with Langevin and Berendsen thermostats. Both Langevin and Berendsen thermostats were used along with nve for time-integrations. With Langevin thermostat, same thing happened, during the increasing of temperature (at temperature around 1800K), the temperature suddenly shoot up drastically and immediately shrink to zero, and simultaneously the volume collapsed to the volume at room temperature. However, with Berendsen thermostat, the simulation managed to drive the temperature to 2000K without any strange phenomena, but after 100000 steps of equilibration in NPT strange phenomena happened, the temperature suddenly shoot up drastically and immediately shrink to zero, and at the same time the volume also followed the same trend, i.e. it expended drastically and immediately collapsed to the volume at room temperature. I really wonder why the system cannot be heated to the desired temperature(2000K)?

Thank you again for your time.

Regards,
Christopher

From your description, it is not thermostat dependent. Rather your system

is undergoing some dramatic transition near that temperature. If it is not
physical (e.g. a melting transition), then it could be b/c your timestep is
too big for that high a temperature (atoms get too close together), or
you are not reneighboring often enough (try delay 0 check yes in neigh_modify),
etc.

Steve

Dear Dr. Steve and all,

I used 1ps timestep, and did re-neighboring at every timestep (every 1 delay 0 check yes) as you can see from my attached input script, which can reproduce the problem. I used VMD to visualize the atoms and they moved very fast at high temperature, hence the probability of 2 atoms overlapping at high temperature is high. So do you think I should use a smaller timestep? The author even used 2ps timestep and able to simulate up to 2200K without phase transition. May I know what possible else may cause this strange phenomena, i.e. during the increasing of temperature and at temperature around 1800K, the temperature suddenly shoot up drastically and then immediately shrink to zero (and kinetic energy as well) and simultaneously the volume collapsed to the volume at room temperature?

Thank you for your time.

Regards,
Christopher

Dear Dr. Steve and all,

I used 1ps timestep, and did re-neighboring at every timestep (every 1 delay 0 check yes) as you can see from my attached input script, which can reproduce the problem. I used VMD to visualize the atoms and they moved very fast at high temperature, hence the probability of 2 atoms overlapping at high temperature is high. So do you think I should use a smaller timestep? The author even used 2ps timestep and able to simulate up to 2200K without phase transition. May I know what possible else may cause this strange phenomena? (i.e. during the increasing of temperature and at temperature around 1800K, the temperature suddenly shoot up drastically and then immediately shrink to zero (and kinetic energy as well) and simultaneously the volume collapsed to the volume at room temperature) ?

Thank you for your time.

Regards,
Christopher

Dear Dr. Steve and all,

I used 1ps timestep, and did re-neighboring at every timestep (every 1 delay
0 check yes) as you can see from my attached input script, which can
reproduce the problem. I used VMD to visualize the atoms and they moved very
fast at high temperature, hence the probability of 2 atoms overlapping at
high temperature is high. So do you think I should use a smaller timestep?
The author even used 2ps timestep and able to simulate up to 2200K without
phase transition. May I know what possible else may cause this strange

that means that there have to be some "safeguards" in the MD program
that prohibit excessive forces. i know that several MD package have an
option to place a limit or cap on forces to avoid this behavior. given the
mass of oxygens and the high temperature, i would consider a 2ps timestep
impossible without some modifications to prohibit excessive forces.

does the paper say what MD code was used?

phenomena, i.e. during the increasing of temperature and at temperature
around 1800K, the temperature suddenly shoot up drastically and then
immediately shrink to zero (and kinetic energy as well) and simultaneously
the volume collapsed to the volume at room temperature?

i think the best way to proceed would be to make an estimate of the box
size at 2200K and then try to simulate at that temperature with nvt or
nve+langevin without going through the npt expansion phase, and see
if you can run a stable MD. there could be something happening that
ticks the npt box rescaling off, and this way you avoid it. if you still get
the crashes, that one has to investigate the potential, or implement a
fix limit that can limit forces/velocities etc. _and_ identify the particles
that cause the problem without crashing the MD every time and making
debugging so difficult.

cheers,
    axel.

Dear Christopher,

AFAICS, your potential function has a well at r=0. Is that possible that
at elevated temperature atoms fall in there?

Best regards,

Vasily Pisarev

Dear Christopher,

AFAICS, your potential function has a well at r=0. Is that possible that
at elevated temperature atoms fall in there?

well, that brings me to another suspicion:

are you sure, christopher, that you have entered the correct
parameters in the correct units?

if i plot the functional form of the potential with you parameters
if get a totally wrong shape. your first parameter in particular seems
to be completely off.

cheers,
    axel.

Dear Dr. Axel, Vasily Pisarev and all,

Thank you for your reply.
Initially, I also suspected that the potential parameters was wrong. However, when I reproduced the lattice constants with the potential parameters, I have no doubt, because I was able to reproduce the lattice constants at temperature range from room temperature to 1500K, and the results were close to the results from the literature. But even so, when I plotted the potential functions (Born and Morse) with the parameters, I saw something wrong with the Born potential. It has a minimum at r=0, as pointed out by Vasily Pisarev. However the Morse potential and overall (Born+Morse) potential are good. I attached the graphs of Born, Morse and Born+Morse potentials vs r and their derivatives(force) vs r for your perusal.
The first parameter in Born style is 4.184 KJ/mol/A times rho(A), so it is equals to rho in Kcal/mol unit. If I use other parameter values or other units for the parameters, I get wrong lattice constants.

The author used MXDORTO, a Japanese MD code developed by Kawamura. It has a manual in Japanese, but unfortunately I don’t understand Japanese language. Perhaps it has a function which could limit the forces/velocities.
I followed your instructions: I made an estimate of volume at 2000K (in point of fact I used the volume from literature), and use NVT thermostat to keep it there. After 20 000 steps, it crashed: the temperature shot up suddenly and immediately fell to zero, and reported atom lost. After that, I used a smaller timestep (0.5) and re-run the MD, this time it crashed after 110 000 steps. Since I want to extract the lattice constants by averaging over large steps at 2000K, so I made an estimation on volume and used NPT thermostat to keep the temperature at 2000K. It crashed too after 28 000 steps. I attached the input script below which can reproduce the problems. By right the thermostat should keep the temperature fluctuating around the desired temperature, I wonder why would this strange phenomena/crash occurred?

Thank you for your time. I really appreciate for your help. These experiences could not be learnt from books.

Regards,
Christopher

strontium titanate system

atom_style charge
units real
dimension 3
boundary p p p

lattice

lattice custom 3.9650 a1 1.0 0.0 0.0 a2 0.0 1.0 0.0 a3 0.0 0.0 1.0 basis 0.0 0.0 0.0 basis 0.5 0.5 0.5 &
basis 0.5 0.5 0.0 basis 0.5 0.0 0.5 basis 0.0 0.5 0.5
region simbox block 0 4 0 4 0 4
create_box 3 simbox
create_atoms 1 box basis 1 1 basis 2 2 basis 3 3 basis 4 3 basis 5 3
group sr type 1
group ti type 2
group o type 3
mass 1 87.62
mass 2 47.867
mass 3 15.999
set group sr charge 1.2
set group ti charge 2.4
set group o charge -1.2

force-field

pair_style hybrid/overlay born 11.0 coul/long 11.0 morse 11.0
pair_coeff * * coul/long
pair_coeff 1 1 born 0.32 0.32 2.396 100.0 0.0
pair_coeff 1 2 born 0.34 0.34 2.253 250.0 0.0
pair_coeff 1 3 born 0.32 0.32 3.124 200.0 0.0
pair_coeff 2 2 born 0.36 0.36 2.110 625.0 0.0
pair_coeff 2 3 born 0.34 0.34 2.981 500.0 0.0
pair_coeff 3 3 born 0.32 0.32 3.852 400.0 0.0
pair_coeff 1 3 morse 2.41 1.18 2.7615
pair_coeff 2 3 morse 4.23 3.82 2.1923

kspace_style ewald 0.0001
pair_modify tail yes shift no

initialize

timestep 1.0
velocity all create 2000.0 1234567 mom yes rot yes dist gaussian
neighbor 4.0 bin
neigh_modify every 1 delay 0 check yes

output

variable N equal step
variable pote equal pe
variable Etotal equal etotal
variable T equal temp
variable Press equal press
variable V equal vol
variable kine equal ke
variable Lx equal lx/4
variable Ly equal ly/4
variable Lz equal lz/4

fix extra all print 1 “{N} {T} {V} {Press} {kine} {pote} {Etotal}" file strontium_titanate.out fix cell all print 1 "{N} {T} {Lx} {Ly} {Lz}” file strontium_titanate.cell

dump snapshot all custom 1000 strontium_titanate.lammpstrj id q type x y z

equilibration & thermalization

#fix NVT all nvt temp 2000.0 2000.0 100.0
#run 200000
#unfix NVT

fix NPT all npt temp 2000.0 2000.0 100.0 iso 1.01325 1.01325 1000.0
run 220000

Date: Sun, 10 Apr 2011 14:52:02 -0400
Subject: Re: [lammps-users] Strange phenomena during increasing of temperature
From: akohlmey@…12…24…
To: pisarevvv@…24…
CC: christopher_gohwf@…8…; [email protected]

Dear Christopher,

AFAICS, your potential function has a well at r=0. Is that possible that
at elevated temperature atoms fall in there?

well, that brings me to another suspicion:

are you sure, christopher, that you have entered the correct
parameters in the correct units?

if i plot the functional form of the potential with you parameters
if get a totally wrong shape. your first parameter in particular seems
to be completely off.

cheers,
axel.

<[email protected]...> wrote:
> Dear Dr. Steve and all,
>
> I used 1ps timestep, and did re-neighboring at every timestep (every 1 delay
> 0 check yes) as you can see from my attached input script, which can
> reproduce the problem. I used VMD to visualize the atoms and they moved very
> fast at high temperature, hence the probability of 2 atoms overlapping at
> high temperature is high. So do you think I should use a smaller timestep?
> The author even used 2ps timestep and able to simulate up to 2200K without
> phase transition. May I know what possible else may cause this strange
 
that means that there have to be some "safeguards" in the MD program
that prohibit excessive forces. i know that several MD package have an
option to place a limit or cap on forces to avoid this behavior. given the
mass of oxygens and the high temperature, i would consider a 2ps timestep
impossible without some modifications to prohibit excessive forces.
 
does the paper say what MD code was used?
 
> phenomena, i.e. during the increasing of temperature and at temperature
> around 1800K, the temperature suddenly shoot up drastically and then
> immediately shrink to zero (and kinetic energy as well) and simultaneously
> the volume collapsed to the volume at room temperature?
 
i think the best way to proceed would be to make an estimate of the box
size at 2200K and then try to simulate at that temperature with nvt or
nve+langevin without going through the npt expansion phase, and see
if you can run a stable MD. there could be something happening that
ticks the npt box rescaling off, and this way you avoid it. if you still get
the crashes, that one has to investigate the potential, or implement a
fix limit that can limit forces/velocities etc. _and_ identify the particles
that cause the problem without crashing the MD every time and making
debugging so difficult.
 
cheers,
    axel.
  

U_born.pdf (61.5 KB)

U_morse.pdf (32.4 KB)

U_born+morse.pdf (62.5 KB)

F_born.pdf (61.6 KB)

F_morse.pdf (32.2 KB)

F_born+morse.pdf (61.6 KB)

I don't know your system or what potential you are using, but a 1 psec
timestep is enormous, by say a factor of 1000, at least for atomic
systems. No atomistic potential in LAMMPS will run properly
with that timestep.

Steve

Dear Dr. Axel, Vasily Pisarev and all,

Thank you for your reply.
Initially, I also suspected that the potential parameters was wrong.
However, when I reproduced the lattice constants with the potential
parameters, I have no doubt, because I was able to reproduce the
lattice constants at temperature range from room temperature to 1500K,
and the results were close to the results from the literature. But
even so, when I plotted the potential functions (Born and Morse) with
the parameters, I saw something wrong with the Born potential. It has
a minimum at r=0, as pointed out by Vasily Pisarev. However the Morse
potential and overall (Born+Morse) potential are good. I attached the
graphs of Born, Morse and Born+Morse potentials vs r and their
derivatives(force) vs r for your perusal.

50 kcal/mole (the barrier height) is about 0.5 eV/atom and your
temperature is ca. 0.2 eV/atom. I think it's quite possible for particle
to cross it. I would replace potential with a tabbed one and add a
repulsive branch at small distances.

Best regards,
Vasily Pisarev

christopher,

Dear Dr. Axel, Vasily Pisarev and all,

Thank you for your reply.
Initially, I also suspected that the potential parameters was wrong.
However, when I reproduced the lattice constants with the potential
parameters, I have no doubt, because I was able to reproduce the lattice
constants at temperature range from room temperature to 1500K, and the
results were close to the results from the literature. But even so, when I

have you considered the fact, that you may be getting the
"right" results for the wrong reasons, e.g. that other components
of the potential parameters keep the atoms so far apart that the
unphysically attractive part of your born potentials does not have
any effect. if you look at the example potential parameters in the
pair_style born documentation, you see that the first parameter
is significantly larger and that give a reasonable shape of the potential.
whereas your born potentials don't look very reasonable.

plotted the potential functions (Born and Morse) with the parameters, I saw
something wrong with the Born potential. It has a minimum at r=0, as pointed
out by Vasily Pisarev. However the Morse potential and overall (Born+Morse)
potential are good. I attached the graphs of Born, Morse and Born+Morse
potentials vs r and their derivatives(force) vs r for your perusal.

The first parameter in Born style is 4.184 KJ/mol/A times rho(A), so it is

are you sure about the "times rho(A)"?

equals to rho in Kcal/mol unit. If I use other parameter values or other
units for the parameters, I get wrong lattice constants.

wrong with respect to experiment, or with respect to other simulations
using the same potentials?

The author used MXDORTO, a Japanese MD code developed by Kawamura. It has a
manual in Japanese, but unfortunately I don't understand Japanese language.
Perhaps it has a function which could limit the forces/velocities.
I followed your instructions: I made an estimate of volume at 2000K (in
point of fact I used the volume from literature), and use NVT thermostat to
keep it there. After 20 000 steps, it crashed: the temperature shot up
suddenly and immediately fell to zero, and reported atom lost. After that, I
used a smaller timestep (0.5) and re-run the MD, this time it crashed after
110 000 steps. Since I want to extract the lattice constants by averaging
over large steps at 2000K, so I made an estimation on volume and used NPT
thermostat to keep the temperature at 2000K. It crashed too after 28 000
steps. I attached the input script below which can reproduce the problems.
By right the thermostat should keep the temperature fluctuating around the
desired temperature, I wonder why would this strange phenomena/crash
occurred?

i maintain, that your born potentials "look wrong" to me
and the symptoms you describe, do match what i would
expect from using potentials with that shape.

i would look up some other systems using similar potentials. the documentation
lists the famous Tosi-Fumi papers. i would try some from there and make some
empirical comparisons.

cheers,
    axel.

Dear Vasily Pisarev,

Thank you for your reply. Ya, the problem come from the Born potential. May I know you do you mean by a tabbed potential? In order to add a repulsive branch at small distance, one have to modify the potential to have a negative gradient at small distance, may I know how would you achieve that?

Thank you.

Regards,
Christopher

Dear Christopher,

Sorry, I mean tabulated potential. And since the attractive branch at
r->0 is anyway unphysical, you can tabulate some repulsive potential
instead. Or you can add LJ potential with a small cutoff radius just to
compensate for the close-range attraction.

Best regards,
Vasily Pisarev