NPT with rigid and non rigid bodies

LAMMPS LAMMPS Mailing List Mirror
1

Dear list,
I am working with CO2 + CH4 mixtures using TraPPE force fields. Since the CO2 is a 3-particle rigid model and the CH4 is a single particle model, I need to simulate a mixture of rigid and non rigid bodies. I’ve tried to carry out my NPT simulations following the three options reported in fix rigid documentation (1 - rigid/npt with dilate all for CO2 and nvt for CH4; 2 - npt with dilate all for CH4 and rigid/nvt for CO2 and; 3 - fix press berendsen for all, rigid/nvt for CO2 and nvt for CH4). All three options gives similar results, with issues in the temperature control (high temperature peaks appears even after long equilibration runs - 10 000 000 fs). I’ve tried to identify the causes by:

  • Using hybrid pair style (lj/cut/coul/long for CO2 and lj/cut for CH4);
  • Applying smaller time steps (0,1 fs);
  • Removing the thermostat and barostat (nve for CH4 and rigid/nve for CO2);
  • Adjusting the Tdamp and Pdamp parameters (in a range from 1 to 10 000).

None of those trials succeed and the temperature still reached very high peaks, moving the average to values higher than desired. It is worth mentioning that the pure compounds simulations (with npt for CH4 and rigid/npt CO2) gives reasonable results, comparable to the NIST standard values. Moreover, the mixture simulation goes ok by adding small flexibility to the CO2 molecule and using a single fix npt for all particles. This works even with very large energy constants for bonds and angles, which means that the CO2 molecule stays almost rigid.

Any help would be very appreciated. My input files are attached.

Best regards,

Cassiano Aimoli

in.lammps (1.74 KB)

lammps.data (92.8 KB)

1 Like
2

- Removing the thermostat and barostat (nve for CH4 and rigid/nve for CO2);

If you can't even run NVE for the system w/out large temperature
excursions, or possibly the timestep is too large, which is easy to
do with small rigid bodies. Putting a thermostat or barostat on it
won't fix the underlying problem. Maybe Trung has another idea. Can
you run a pure CO2 or pure CH4 system w/out issue?

Steve

3

Hi Cassiano,
Which property are you looking to compute in the end? The rigid bonds
in the TraPPE FF are only an advantage when combined with
Configurational Bias Monte Carlo techniques. For response properties
like viscosities, it has been shown by the original authors of the
TraPPE FF that its better to switch to harmonic bonds taken from the
Amber FF for example. Harmonic bonds are a better representation when
it comes to calculating properties such as the one I have mentioned.
Keep that in mind.
Cheers,
Carlos

4

Cassiano,
In stead of using rigid try to use a high bondLength (~5000 for Kb) value for the flexible. See what happens.
Amir

5

Simply increasing the rigidity of the bond is not a good idea Amir.
Doing this implies reducing accordingly the time step to get proper
integration of the EOM. One faces simulations taking a long time.
Carlos

6

Hello Carlos,
Thanks for your comments. At the moment I am trying to perform a comparison of CO2 models results in density calculation. Besides TraPPE and other rigid models, I am also simulating flexible and coarse grain models, for the sake of comparison.

Regards,

Cassiano

7

Well, you'll probably find that using a united atom representation of
CH4 will work just fine for density purposes. In the end your system
is nothing but a bunch of balls and rods moving around.
Carlos

8

Dear Amir,
Thank you for your response.
I actually did it (using large bond and angle constants - 10 000), and it runs ok with fix npt. Nevertheless, I still want to run this with rigid CO2.

Regards,

Cassiano

9

Dear list,
I am working with CO2 + CH4 mixtures using TraPPE force fields. Since the
CO2 is a 3-particle rigid model and the CH4 is a single particle model, I
need to simulate a mixture of rigid and non rigid bodies. I've tried to
carry out my NPT simulations following the three options reported in fix
rigid documentation (1 - rigid/npt with dilate all for CO2 and nvt for CH4;
2 - npt with dilate all for CH4 and rigid/nvt for CO2 and; 3 - fix press
berendsen for all, rigid/nvt for CO2 and nvt for CH4). All three options
gives similar results, with issues in the temperature control (high
temperature peaks appears even after long equilibration runs - 10 000 000
fs). I've tried to identify the causes by:

- Using hybrid pair style (lj/cut/coul/long for CO2 and lj/cut for CH4);

there is no reason why this should make a difference.
the forces and energies are the same whether you use
lj/cut/coul/long for all or a hybrid pair style.

- Applying smaller time steps (0,1 fs);

...and does an even smaller one (0.01fs) help?

- Removing the thermostat and barostat (nve for CH4 and rigid/nve for CO2);
- Adjusting the Tdamp and Pdamp parameters (in a range from 1 to 10 000).

None of those trials succeed and the temperature still reached very high
peaks, moving the average to values higher than desired. It is worth

ok. that means, that you have instabilities in the integration and/or atoms
getting too close. how did you derive the interaction cross-terms between
the CO2 atoms and the CH4 sites? simply by applying a mixing rule?
or did you derive specific parameters?
it is quite conceivable, that with applying a mixing rule, you get too stiff
mixed interactions. in the case of pure CO2 you also have the self-self
repulsion due to coulomb and the C vs O attraction, which would prefer
a slipped parallel configuration. the CH4 site isn't set up for that.

in order to investigate this, you should try to analyse the trajectory around
such a "spike" and check which conformation exactly is problematic.
the fix may be simply to make the cross interactions between CO2 and
CH4 softer (on the repulsive part while keeping the attractive part similar).

mentioning that the pure compounds simulations (with npt for CH4 and
rigid/npt CO2) gives reasonable results, comparable to the NIST standard
values. Moreover, the mixture simulation goes ok by adding small flexibility
to the CO2 molecule and using a single fix npt for all particles. This works
even with very large energy constants for bonds and angles, which means that
the CO2 molecule stays almost rigid.

which hints at having some close contact between CO2 and CH4
leading to an instability in the rotation part of the rigid integrator
and then pushing sites too close.

axel.

10

Hi Cassiano,

just want to add to Axel's comments regarding the possible instability
of the rotational part of the rigid integrator, I suggest you try
different (greater) values of Titer via the option tparam. The default
value of Titer is 1, you can increase it to 3, 4, 5, etc. and see if
the issue persists.

Regards,
-Trung

PS. Is there any special reason for using ewald instead of pppm for
kspace style?

11

Dear all,
Applying smaller timesteps (0.05fs) and adjusting the Titer to 3 appears to fix the problem. Thank you all for the help.

Trung,
There is no special reason for using Ewald. Actually, I am using PPPM and it is working fine. The Ewald appears on my input file because I was trying to track any possible cause of the temperature instability (changing kspace style, pair styles, initial conditions etc. I just forgot to change it back to PPPM before submitting this…)

12

First of all, I think there are probably ways to get the system to
work using "fix rigid". There are many things to try before you give
up.
--- however ---
Instead I would recommend using "fix shake".

I read somewhere that "fix rigid" requires much smaller timesteps than
"fix shake". (It was probably in a post by Axel.) It seems like
SHAKE does not work for larger cyclic molecules. However it works
well to constrain small (3-site rigid) water molecules. So you should
be able to use it to make CO2 molecules rigid as well.
http://lammps.sandia.gov/doc/fix_shake.html

---example---
There is an example how to use shake for SPCE water below:
group spce type 1 2
fix fSHAKE spce shake 0.0001 10 100 b 1 a 1
timestep 2.0
fix fxnpt all npt temp 300.0 300.0 100.0 iso 1.0 1.0 1000.0 drag 1.0
run 50000
(Note: Unlike "fix rigid", there is no need to use "neigh_modify
exclude" if you use "fix shake".)

These examples were taken from here:
http://www.moltemplate.org/examples/nanotube+water/spce.lt
http://www.moltemplate.org/examples/waterSPCE+Na+Cl/run.in.npt
(Just replace "@atom:O", "@atom:H", "@bond:OH" and "@angle:HOH" with
"1", "2", "1", and "1".)

I hope this helps
Andrew

13

Thanks Andrew,
But according to previous messages (http://lammps.sandia.gov/threads/msg14914.html) the fix shake does not work properly for linear molecules, like CO2.
Regards,

Cassiano

14

Thanks Andrew,
But according to previous messages
(http://lammps.sandia.gov/threads/msg14914.html) the fix shake does not work
properly for linear molecules, like CO2.

Interesting. I did not know. Thanks!

In my last email I mentioned there are other things you could try to
get "fix rigid" working. Here they are:

1) Use "fix npt" (instead of nvt) for CH4.
In your input script I noticed you are mixing "fix npt" with "fix nvt".

2) Try leaving out "dilate all" from the end of the "fix npt/rigid" command.

In other words, try using:

fix 1 CO2 rigid/npt molecule temp 373.15 373.15 100.0 iso 300.0 300.0 1000.0

fix 2 CH4 npt temp 373.15 373.15 100.0 iso 300.0 300.0 1000.0 dilate CH4

Earlier, you were using
"fix 1 CO2 rigid/npt molecule ... dilate all"

Wouldn't using this command change the internal distances in your CO2
molecules? (So would "dilate CO2") Is this really what you want
LAMMPS to do? (Does this even make sense? This seems wrong.)

The documentation for "fix rigid/npt" does not explain what happens if
you leave out the "dilate" keyword. (Is this equivalent to "dilate
all" or "dilate none"?) Anyway, try removing "dilate all", and see
what happens. (I would like to know what this does.)

3) Alternatively, if this fails, try using the old "fix rigid" on your
CO2 molecules:

fix 1 CO2 rigid molecule
fix 2 CH4 npt temp 373.15 373.15 100.0 iso 300.0 300.0 1000.0 dilate CH4

To get the temperature correct, you might need to use a combination of
"compute temp" and "fix_modify"

compute tempCH4 CH4 temp
fix_modify 2 temp tempCH4

(Credit to Trung for this suggestion.)

3) If you are using fix rigid (or rigid/npt), it's a good idea to
start the simulation with the CO2 molecules in the lowest-energy
conformation (straight, with the bond-lengths at the minimum-energy
value). This should turn off internal forces in CO2 molecules which
might contribute to the virial. (If you only have one rigid object
you can use "neigh_modify exclude" to get rid of these contributions.
But you have many rigid objects.) This is probably not necessary, but
it does not hurt.

4) Of coarse, after all this, you are still having trouble with "fix
rigid", I suppose I would try using "fix shake" to constrain the C-O
bonds and use harmonic angle forces to maintain the 180 angle. If
that is not numerically stable, I suppose you could use SHAKE to
constrain the O-O distance, and remove the C-O distance constraints,
and the harmonic angle constraint.

In any case, using a 0.05 fs timestep seems painful. Hopefully you
find a better way.

Andrew

15

As much as I understand the rigidity of the model will not affect the density significantly (especially at gaseous densities, which your input configuration seems to be at). What will be important is the LJ parameters for the model. As someone said before, TraPPE was developed for use with CBMC moves in an MC simulation. A comparison of a flexible version of TraPPE with other flexible CO2 models using MD is a valid comparison of the model parameters and how they predict density.

16

Using a different code (dlpoly) I was able to run TraPPE CO2 with shaked bond lengths and a harmonic potential for the angle. I know shake in the two programs is implemented differently though. In dlpoly it is possible to shake whole alkane molecules without a problem.