Equilibration

George_Smith · February 13, 2013, 6:25pm

Hello,

I am using reaxFF. I have a few questions:

(i) I figured that if I do NTP with 0.001 fs timestep, I have small fluctuation in pressure (+/- 500 atmosphere) and temperature… even with 0.1 fs timestep the pressure fluctuation is enormously large (100K atmosphere). The average pressure is always zero. Wonder if it is normal ?

(ii) After using NPT and heating the system to 500K, I run it a million step with the same above parameters (timestep = 0.001 fs). I see that temperature, energy, enthalpy are stable. However, fluctuations in volume and pressure is as before.

(iii) Finally, I run NVE to see if the system is equilibrated ? This time I run with timestep = 0.1 fs for 100K steps (10 ps). The temperature started rising and rose from 500 K to 700 K. Pressure fluctuation is now +/- 80K atmosphere. Being NVE, the volume and energy is almost constant.

From (ii) I would believe, the system is equilibrated. However, from (iii) it looks like my system is not equilibrated. If I choose the same small step size, I would need to do 10 million steps to equilibrate the system (because I typically see 10-25 ps for the equilibration in literature).

(iv) I have tried “velocity all create 300 12345 mom yes rot no” after minimization to quickly heat the system. But it never worked for my system. The program stops with error message. Basically system breaks down with huge movement of atoms. I guess SHAKE is not allowed for reaxFF so that is not an option.

I would like some speedup in the dynamics for my system and I am looking for clues if I could be doing something wrong or inefficient. Is there something about reaxFF as implemented within LAMMPS that could be the issue ?

Thanks,

George

Carlos_Campana · February 13, 2013, 6:36pm

George,
One very important thing to keep in mind with reaxff in lammps is to
make sure that you are using the same params as the original
implementation of the ff parametrization you are employing (assuming
you are not the developer of such a parametrization). The fortran
reaxff library in lammps is hard-coded with the param set of the 2008
hydrocarbon combustion paper of Van Duin et al. If your reaxff params
are not the same, you should use the USER-REAXC package combined with
a control file where the new params can be specified. Not employing
the correct params could affect not only your energy values but also
energy conservation within NVE runs, etc.
Just to keep in mind.
Carlos

Ray_Shan · February 13, 2013, 7:18pm

Hi George,

To add to Carlos' comments, the control settings for pair_style reax/c
are exactly the same as those in the stand-alone reaxff code, so using
the NULL keyword "cfile" command is the best starting point. More
comments below, thanks.

Ray

Hello,

I am using reaxFF. I have a few questions:

Which one? reax or reax/c?

(i) I figured that if I do NTP with 0.001 fs timestep, I have small
fluctuation in pressure (+/- 500 atmosphere) and temperature.... even with
0.1 fs timestep the pressure fluctuation is enormously large (100K
atmosphere). The average pressure is always zero. Wonder if it is normal ?

I think 0.001 fs is an unnecessary small time step for reaxff.
Pressure fluctuation is normal, and smaller the system larger the
fluctuation.

(ii) After using NPT and heating the system to 500K, I run it a million step
with the same above parameters (timestep = 0.001 fs). I see that
temperature, energy, enthalpy are stable. However, fluctuations in volume
and pressure is as before.

This adds up to only 1 ps. Hard to say if it is really equilibrated,
and by equilibrated I mean minimal fluctuation is observed.

(iii) Finally, I run NVE to see if the system is equilibrated ? This time I
run with timestep = 0.1 fs for 100K steps (10 ps). The temperature started
rising and rose from 500 K to 700 K. Pressure fluctuation is now +/- 80K
atmosphere. Being NVE, the volume and energy is almost constant.

"volume" is almost constant with NVE? Then you are definitely doing
something wrong otherwise volume should not fluctuate at all. If
total energy remains more or less constant, then it should be fine.

From (ii) I would believe, the system is equilibrated. However, from (iii)
it looks like my system is not equilibrated. If I choose the same small
step size, I would need to do 10 million steps to equilibrate the system
(because I typically see 10-25 ps for the equilibration in literature).

But they don't usually use a timestep as small as 0.001 fs.

(iv) I have tried "velocity all create 300 12345 mom yes rot no"
after minimization to quickly heat the system. But it never worked for my
system. The program stops with error message. Basically system breaks down
with huge movement of atoms. I guess SHAKE is not allowed for reaxFF so that
is not an option.

velocity create should work, and if your run stops with error message,
then something terribly wrong must have occurred. What kind of error
message?

George_Smith · February 13, 2013, 9:30pm

Hi George,

To add to Carlos' comments, the control settings for pair_style reax/c
are exactly the same as those in the stand-alone reaxff code, so using
the NULL keyword "cfile" command is the best starting point. More
comments below, thanks.

Ray

Hi Ray,

Thanks ! I have not built lammps with user-reaxc. I shall do that soon and
see how I can run it. Just curious if it is preferable to use pair_style
reaxc to pair_style reax ? or they are exactly the same... from my reading
I thought there is no difference between them except for the implementation.

> Hello,
>
> I am using reaxFF. I have a few questions:

Which one? reax or reax/c?

reax

> (i) I figured that if I do NTP with 0.001 fs timestep, I have small
> fluctuation in pressure (+/- 500 atmosphere) and temperature.... even
with
> 0.1 fs timestep the pressure fluctuation is enormously large (100K
> atmosphere). The average pressure is always zero. Wonder if it is normal
?

I think 0.001 fs is an unnecessary small time step for reaxff.
Pressure fluctuation is normal, and smaller the system larger the
fluctuation.

So, you believe a pressure fluctuation of +/- 100,000 atmosphere is ok. I
see some waviness in structure which I interpreted as an effect of high
pressure fluctuations.

>
> (ii) After using NPT and heating the system to 500K, I run it a million
step
> with the same above parameters (timestep = 0.001 fs). I see that
> temperature, energy, enthalpy are stable. However, fluctuations in volume
> and pressure is as before.

This adds up to only 1 ps. Hard to say if it is really equilibrated,
and by equilibrated I mean minimal fluctuation is observed.

>
> (iii) Finally, I run NVE to see if the system is equilibrated ? This
time I
> run with timestep = 0.1 fs for 100K steps (10 ps). The temperature
started
> rising and rose from 500 K to 700 K. Pressure fluctuation is now +/- 80K
> atmosphere. Being NVE, the volume and energy is almost constant.
>

"volume" is almost constant with NVE? Then you are definitely doing
something wrong otherwise volume should not fluctuate at all. If
total energy remains more or less constant, then it should be fine.

I was inaccurate in my earlier mail. The volume is constant. Energy does

fluctuate very little. But the pressure fluctuation is very large ~ +/-
80,000 atmosphere and temperature increased from 500 Kelvin to 700 Kelvin.

From (ii) I would believe, the system is equilibrated. However, from (iii)
> it looks like my system is not equilibrated. If I choose the same small
> step size, I would need to do 10 million steps to equilibrate the system
> (because I typically see 10-25 ps for the equilibration in literature).

But they don't usually use a timestep as small as 0.001 fs.

>
> (iv) I have tried "velocity all create 300 12345 mom yes rot no"
> after minimization to quickly heat the system. But it never worked for my
> system. The program stops with error message. Basically system breaks
down
> with huge movement of atoms. I guess SHAKE is not allowed for reaxFF so
that
> is not an option.

velocity create should work, and if your run stops with error message,
then something terribly wrong must have occurred. What kind of error
message?

ia(i1,2) = 21 reax_defs.h::MBONDDEF = 20 after
  705 of 2232958 pairs completed.
STOP Too many bonds on atom. Increase MBONDDEF
ia(i1,2) = 21 reax_defs.h::MBONDDEF = 20 after
  540 of 2232856 pairs completed.
STOP Too many bonds on atom. Increase MBONDDEF
ia(i1,2) = 21 reax_defs.h::MBONDDEF = 20 after
  729 of 2307413 pairs completed.
STOP Too many bonds on atom. Increase MBONDDEF

Ray_Shan · February 14, 2013, 2:07am

George,

Thanks ! I have not built lammps with user-reaxc. I shall do that soon and
see how I can run it. Just curious if it is preferable to use pair_style
reaxc to pair_style reax ? or they are exactly the same... from my reading I
thought there is no difference between them except for the implementation.

reax/c compiles easier, allocates memory in a more efficient way, has
more freedom in control settings, and runs faster.

I think 0.001 fs is an unnecessary small time step for reaxff.
Pressure fluctuation is normal, and smaller the system larger the
fluctuation.

So, you believe a pressure fluctuation of +/- 100,000 atmosphere is ok. I
see some waviness in structure which I interpreted as an effect of high
pressure fluctuations.

No, I don't think a fluctuation of 100K atom is normal. What I meant
was pressure fluctuation is a normal phenomenon, and fluctuations are
larger for smaller systems. You probably have a relatively small
system.

velocity create should work, and if your run stops with error message,
then something terribly wrong must have occurred. What kind of error
message?

ia(i1,2) = 21 reax_defs.h::MBONDDEF = 20 after
705 of 2232958 pairs completed.
STOP Too many bonds on atom. Increase MBONDDEF
ia(i1,2) = 21 reax_defs.h::MBONDDEF = 20 after
540 of 2232856 pairs completed.
STOP Too many bonds on atom. Increase MBONDDEF
ia(i1,2) = 21 reax_defs.h::MBONDDEF = 20 after
729 of 2307413 pairs completed.
STOP Too many bonds on atom. Increase MBONDDEF

This error only shows up when maximum number of bonds around a atoms
has exceeded 20. You may increase the #defined variable in /lib/reax
and re-compile, but more often it indicates some bad dynamics resulted
from bad structure. I'd suggest you examine your structure if atoms
were too close to each other.

Cheers,
Ray

gilbert1 · February 14, 2013, 3:45am

George and Ray,

I work with organic polymers (CHNOS) and MBONDDEF always needs to be at least 30. When I build solvated polymers with lots of water I usually need to increase the hbond term also.

If you are working with more than C and H make certain that the order of atom type on the pair coeff line: * * ffield.reax 1 2 3 4 5 (this line is for C H O N and S)

matches the order of atom types (masses) in your data file.It is very easy to miss this and normally results in a lack of temperature control because you are applying the wrong parameters to the atoms in your structure. The methods I use for building my polymers often interchanges the order of oxygen and nitrogen.

I usually do NVE for equilibration and have not had to use a timestep less than 0.1 fs. But you should also look at the the neigh_modify command to make sure you are updating the neighbor list in a timely manner.

You did not say if you were doing gas phase or liquid phase simulations, but large pressure fluctuations are common in the liquid phase but decrease with the number of particles simulated.

Kevin

George_Smith · February 15, 2013, 5:17pm

Thanks,

George

George and Ray,

I work with organic polymers (CHNOS) and MBONDDEF always needs to be
at least 30. When I build solvated polymers with lots of water I usually
need to increase the hbond term also.

If you are working with more than C and H make certain that the order
of atom type on the pair coeff line: * * ffield.reax 1 2 3 4 5 (this
line is for C H O N and S)

matches the order of atom types (masses) in your data file.It is very
easy to miss this and normally results in a lack of temperature control
because you are applying the wrong parameters to the atoms in your
structure. The methods I use for building my polymers often interchanges
the order of oxygen and nitrogen.

I usually do NVE for equilibration and have not had to use a timestep
less than 0.1 fs. But you should also look at the the neigh_modify
command to make sure you are updating the neighbor list in a timely
manner.

I use the following, I guess I picked up from other example files:
pair_style reax 10.0 0 1 1.0e-6
neighbor 2.5 bin
neigh_modify every 10 delay 0 check no

You did not say if you were doing gas phase or liquid phase
simulations, but large pressure fluctuations are common in the liquid
phase but decrease with the number of particles simulated.

It is a periodic solid.

gilbert1 · February 15, 2013, 6:39pm

George,

It is not the pair_style line that is important - I use the same one. It is the pair_coeff line where you map the order of atom types in your file to the order of atom types in the reax parameter file. Each ffield.reax file has a specific order of atom types (atomic masses). For the 2010 hydrocarbon file is is C,H,O,N,S. For the Carbon/Ni file it is C,H,Ni. You can read the file (ffield.reax) into a text editor to determine the order of parameters.

On the pair_coeff line you tell lammps how to map the atoms in your input file to the atom types in the ffield.reax file. So if my input file has Carbon as the first element, followed by Nitrogen, then Oxygen and finally Hydrogen I would have a pair_coeff line of:

pair_coeff * * ffield.reax 1 4 3 2

I hate to admit how many times I forgot to make this change in the Lammps input file as I built new polymers because the only way I would find out was the lack of temperature control.

I also got in the habit of first minimizing my lammps input file with reax and then comparing the output structure with the input structure bond by bond. If the distances changed I would know there was a problem. The original structure had been minimized with an MMFF94 force field and for most bonds there was very little difference in the structures minimized with MMFF and ReaxFF.

Kevin