lj/charmm/coul/charmm problem with tip3p

I ran your script with fix nvt and fix shake replaced by
fix nve. It seems to run fine and equilibrate to some
high temperature (2800K). Then I ran with fix nvt and
fix shake and it seemed to be equlibrating to 300K
as requested.

So when you say "frozen" in your original email, you
don't mean that the system loses energy over time
and the temperature decreases?

Do you simply mean that the water molecules aren't diffusing?
That could be a lot of things. Hi density. Using a longer
cutoff (you are using 12 with charmm and 10 with coul/cut
in your input script), etc. Are you using the SHAKE params
for water with a cutoff and not for long-range water? And with
the correct cutoff?

Steve

Thanks for helping me with my problem. With frozen I mean that the water molecules aren't diffusing at 300 K. Density is adjusted to 1000 kg/m^3 and this works also fine in an npt simulation, but molecules were still not diffusing.

For the cutoff values I took the parameters from Lopes et al. [2006] (because in the end I want to simulate a-quartz plus water as described in his paper). What do you mean with SHAKE params for water with a cutoff? Are there different ones? My bond angle and bond distance is fixed with SHAKE quite well with nearly no deviations from initial values.

In a nve simulation my energy is not conserving. I start from an equilibrated one at 270 K and after some steps (500000) my system explodes. I attached some graphs with energy, pressure and temperature vs. time.

Again, when I'm using a simple lj/cut/coul/cut pair style with a cutoff of 10.0 Ang everything is fine.

I also asked Brian if he solved the problem, but he still has no solution for it. Here is what he said:

>> Robert,

>> No, unfortunately I never got it to give correct behavior for water. It was not critical that I got it working, since I could >> use lj/charmm/coul/long with a slab geometry. I did not want to spend too much time trying to figure it out.

>> I just ran another quick NPT simulation with a newer version (11Feb11), and it still looks like it is freezing at 1500K. I did not try to run a NAMD or CHARMM simulation with the same cut off distances. That would prove that this behavior is not due to the switched coulomb potential, but likely due to a bug.

>> Good luck,
>> Brian

Greetings Robert

time_vs_temp.png

time_vs_epair.png

time_vs_etot.png

time_vs_press.png

See section 4.7 of the manual. TIP3P has different parameterizations
for rigid vs flexible water, and for cutoff vs long-range Coulombics.
You are also not free to set the cutoff to whatever you wish (even
for long-range, there is still a Vdwl cutoff), since the TIP3P model
is parameterized for a specific cutoff. I don't recall what that cutoff
is, but we should list it in that section of the manual.

Steve

Steve,

I don't know what I'm missing. After extensive studying of literature regarding TIP3P water models and simulation of them I still cannot figure out where the mistake is in my simulation. I ran a 1ns simulation of TIP3P water with lj/charmm/coul/long and the appropriate parameters from the manual, the result was a self-diffusivity of 4.22e-9 m^2/s which is close to other results in the literature (see Mark et al. 2001, J. Phys. Chem. A 2001) and the small difference could be explained through the discrepancy in the density and temperature. I used 1.0 g/cm^3 and 300 K and Mark et al. used 0.998 g/cm^3 and 299 K. Nonetheless the value I get is in a quite good agreement.

I looked in the Lopes et al. and MacKerell et al. 1998 publications to find some details about the modified TIP3P model and I choose my cutoff and neighbor list parameters according to their values (10-12 Ang force switching, rigid water molecules with SHAKE and 2 Ang neighbor skin). In contrast to the previous simulation with lj/charmm/coul/charmm I get a diffusivity with lj/charmm/coul/charmm of 6.377e-13 m^2/s which is far too small. I also talked with a college who had the same problems with TIP3P and lj/charmm/coul/charmm some time ago, he said when simulating near a surface this leads to unrealistic water layering and it seems it has something to do with the cutoff distance.

My intention is to simulate Quartz and TIP3P water as Lopes et al. did in their publication. Unfortunately they were using the CHARMM suite, this is not possible in our group and I thought it was a good idea to use LAMMPS because it has all the functional forms I need for a CHARMM style simulation.

Now I ask myself If I can make a good TIP3P simulation with lj/charmm/coul/long why this is not possible with lj/charmm/coul/charmm and the appropriate parameters for it?

Best regards
Robert

robert,

Steve,

I don't know what I'm missing. After extensive studying of literature
regarding TIP3P water models and simulation of them I still cannot
figure out where the mistake is in my simulation. I ran a 1ns simulation
of TIP3P water with lj/charmm/coul/long and the appropriate parameters
from the manual, the result was a self-diffusivity of 4.22e-9 m^2/s
which is close to other results in the literature (see Mark et al. 2001,
J. Phys. Chem. A 2001) and the small difference could be explained
through the discrepancy in the density and temperature. I used 1.0
g/cm^3 and 300 K and Mark et al. used 0.998 g/cm^3 and 299 K.
Nonetheless the value I get is in a quite good agreement.

I looked in the Lopes et al. and MacKerell et al. 1998 publications to
find some details about the modified TIP3P model and I choose my cutoff
and neighbor list parameters according to their values (10-12 Ang force
switching, rigid water molecules with SHAKE and 2 Ang neighbor skin). In
contrast to the previous simulation with lj/charmm/coul/charmm I get a
diffusivity with lj/charmm/coul/charmm of 6.377e-13 m^2/s which is far
too small. I also talked with a college who had the same problems with
TIP3P and lj/charmm/coul/charmm some time ago, he said when simulating
near a surface this leads to unrealistic water layering and it seems it
has something to do with the cutoff distance.

My intention is to simulate Quartz and TIP3P water as Lopes et al. did
in their publication. Unfortunately they were using the CHARMM suite,
this is not possible in our group and I thought it was a good idea to
use LAMMPS because it has all the functional forms I need for a CHARMM
style simulation.

please note that the TIP3P potential in CHARMM is actually _not_
the TIP3P potential, but a modified one that is also known as TIPS2.
i suspect that a lot of the confusion about TIP3P is resulting from this,
also people forget that some properties are very difficult to converge well.

self-diffusion is particularly nasty, as it can appear to be converged.
check out: http://klein-group.icms.temple.edu/akohlmey/files/talk-trieste2004-water.pdf

Now I ask myself If I can make a good TIP3P simulation with
lj/charmm/coul/long why this is not possible with lj/charmm/coul/charmm
and the appropriate parameters for it?

there is one thing that i don't get. why _do_ you want to simulate
with lj/charmm/coul/charmm when lj/charmm/coul/long is the better
solution in any case. i can confirm the report of your colleague that
any kind of coulomb force truncation - however smart or smooth -
will cause a _significant_ spurious potential perpendicular to the surface.

check out: http://link.aip.org/link/doi/10.1063/1.474295

cheers,
    axel.

I want to use lj/charmm/coul/charmm because Lopes et al. used this in his publication from 2006. It is written there that "Long-range interactions in CHARMM were truncated due to the Ewald summation algorithm being incompatible with the covalent linkages between the primary and image atoms. [...] Force shift and force switch smoothing of the electrostatic and Lennard-Jones (LJ) interactions [...] were performed with the force switch initiated at 10 Ang, nonbonded interactions truncated at 12 Ang and nonbond pair lists were maintained heuristically to 14 Ang."

On the one hand I don't understand why Ewald summation should be incompatible with covalent linkages between primary and image atoms and on the other hand I want to do the same simulation as Lopes did, which means that I have to use the same parameters.

Beside I plotted and annotated the energy curve of separating two water molecules. There is something strange happening in the smoothing area in both plots. One is using lj/charmm/coul/long and the other is using lj/charmm/coul/charmm.

Robert

tip3p-dimer-lj-charmm-coul-long-energy.eps (15.2 KB)

tip3p-dimer-lj-charmm-coul-charmm-energy.eps (15.4 KB)

I don't know the answer to your Q. If you suspect a bug in
the LAMMPS implementation of the potential (possible, this
one is likley not used very much - as Axel says, most people
do long-range water, not cutoff), then you could use pair_write
to plot the potential and see if anything looks suspect. That actually
evaluates the potential at every distance, so it is a good test.

Steve

I would guess the bumps in the plots are simply what happens
when some atoms in the water molecules move across the cutoff.
LAMMPS uses an atom-based cutoff, not a group (or molecule-based) cutoff,
meaning you can have some atoms in the two waters interacting
and not others.

The bad bumps in the cutoff case are why people use long-range
Coulombics (in almost all models, not just water).

Steve

I'm a little bit confused and disappointed. What is now the essence of our discussion? Is there a way to simulate modified TIP3P water with a spherical shifted-force cutoff method in LAMMPS as it is done in the literature (but with other MD programs) or is it not possible with LAMMPS?

I also plotted the energy curves with the pair_write command and there is till the bump in the lj/charmm/coul/charmm style between 10 and 12 Ang, which is neither in lj/charmm/coul/long nor in lj/cut/coul/cut.

Best regards
Robert

tip3p-pair-energy-lj-cut-coul-cut.eps (28.7 KB)

tip3p-pair-energy-lj-charmm-coul-charmm.eps (28.7 KB)

tip3p-pair-energy-lj-charmm-coul-long.eps (29 KB)

I'm sorry, I forgot to include the charges. Here is another plot of pair_write but with correct charges. The bump in lj/charmm/coul/charmm isn't ok, or?

Regards
Robert

tip3p-pair-energy-lj-cut-coul-cut.eps (27.7 KB)

tip3p-pair-energy-lj-charmm-coul-charmm.eps (27.7 KB)

tip3p-pair-energy-lj-charmm-coul-long.eps (27.9 KB)

Again sorry, the filename is also erroneous. The title of the respective plot is however correct.

sorry it's late, so don't get confused by me because I am it already :slight_smile:

Regards
Robert

I'm sorry, I forgot to include the charges. Here is another plot of
pair_write but with correct charges. The bump in lj/charmm/coul/charmm isn't
ok, or?

no this does not look ok. :wink:

this looks a lot like some pre-factor is wrong.

axel.

well, I finally managed to run a simple TIP3P water simulation and also the lopes forcefield works now as expected. The reason was indeed the potential switching, which is very unsuitable in many CHARMM simulation and a force shifting or switching should be used. First I was a little bit confused about potential and force shifting or switching, but after reading the Steinbach & Brooks publication from 1994 a lot of misunderstandings was eliminated.

Because of the above mentioned reason I programmed a force switching method as implemented in recent CHARMM versions for LAMMPS (namely lj/fsw/coul/fsw) according to the Steinbach and Brooks [1994] publication and now everything works fine :slight_smile:

Best regards
Robert

glad you got to the bottom of this - do you have
a new (or modified) pair style that does this
that we could release in LAMMPS?

Steve

sure, what do you need for that?

It is based on the lj/charmm/coul/charmm pair style and I tried to code it in regular LAMMPS fashion. So every constant which is not distance dependent is pre-initialized. I also checked if the forces and energies are in agreement with theory, so I hope there is no bug in it :slight_smile:

regards
robert

I just need the *.cpp and *.h file and a doc page with
the altered formula. For the latter you can probably just
add it to doc/pair_charmm.txt.

I attached my *.cpp and *.h files, a modified doc/pair_charmm.txt with modified equations and a plot of the force switching energy and force curves. I hope thats enough. If there are some problems or misunderstandings in the doc/pair_charmm.txt please let me know.

by the way the force switching for vdW considered in the plots (namely between two TIP3P H atoms) is in comparison to the coulombic forces much smaller. That's why in the vdW energy plot both curves (theoretical and unswitched) are nearly identical, the difference between these two is in the order of 1e-18.

regards
robert

pair_lj_fsw_coul_fsw.tar.gz (7.67 KB)

coul_fsw_energy.eps (90.9 KB)

coul_fsw_force.eps (43.7 KB)

lj_fsw_force.eps (47 KB)

lj_fsw_energy.eps (91.6 KB)

Recently I red the MacKerell [1998] publication about CHARMM and there was written that "Nonbonded interaction terms are included for all atoms separated by three or more covalent bonds. No general scaling of the electrostatic or Lennard-Jones interactions for atoms separated by three bonds (the so-called 1-4 term) is used. In specific cases there is scaling of the 1-4 Lennard-Jones term;".

If I want to use the all-atom CHARMM forcefield from MacKerell in LAMMPS, is the special_bonds charmm command then still correct? As far is I understand this will only consider 1-4 interactions from dihedrals in the data file and the rest is ignored.

The Lennard-Jones scaling is already implemented in the lj/charmm pair_style but does this work independently from the special_bonds charmm command?

regards
robert

The special_bonds command has a dihedral yes/no option
which differentiates between 1-4 interactions that are part
of a specific listed dihedral or not. I think this is what
you are referring to below?

Steve

OK, thanks for the explanation. But there is one thing left. Let us assume I have a pair_style like this

pair_style lj/charmm/coul/charmm
pair_coeff 1 1 100.0 2.0 150.0 3.5
special_bonds charmm dihedral yes

is then the already scaled LJ energy scaled additionally with the weighting factor of the dihedral?

Robert