possible bug in lj/charmm/coul/charmm ?

Hi All,

I’m simulating water using lj/charmm/coul/charmm. The script is attached, which is just the same as in.spce in lammps/bench/POTENTIALS except the pair-potential is changed from lj/cut/coul/long to lj/charmm/coul/charmm.

The problem is water molecules are just vibrating around their lattice positions with out moving anywhere. Changing the potential to any other is working fine.

In the mailing list also a couple of users also reported the freezing problems with using lj/charmm/coul/charmm potential.

Lammps script and data file are attached with this email. water.xyz file can be used to visualize the problem.

Thanks.

water.data (65.4 KB)

water.lmp (563 Bytes)

Hi All,

I'm simulating water using lj/charmm/coul/charmm. The script is attached,
which is just the same as in.spce in lammps/bench/POTENTIALS except the
pair-potential is changed from lj/cut/coul/long to lj/charmm/coul/charmm.

The problem is water molecules are just vibrating around their lattice
positions with out moving anywhere. Changing the potential to any other is
working fine.

In the mailing list also a couple of users also reported the freezing
problems with using lj/charmm/coul/charmm potential.

through making this switch you also change your water potential.
different water potentials have different melting points.
so, if you want to determine whether there is a bug or not, you have
to prove that the forces are wrong, i.e. not according the potential
function in the documentation or that that function is not correct.

axel.

Hi Axel,

Thanks for the reply.

Even at 1000K temp, the water molecules are vibrating only around their lattice cites with out moving anywhere. I have observed similar behaviour with ionic liquids also. I noticed this problem few months ago, but did not investigate further.

So I’m sure, its not a melting point problem, but a possible bug …

Best,

Hi Axel,

Thanks for the reply.
Even at 1000K temp, the water molecules are vibrating only around their
lattice cites with out moving anywhere. I have observed similar behaviour
with ionic liquids also. I noticed this problem few months ago, but did not
investigate further.

so what do you expect to happen now? especially, after you publicly
state that you don't consider it worth spending time on it yourself.

more importantly, you don't note which version and which platform this
is on. for all we know, you could be running a massively outdated
version of LAMMPS and the possible bug is possibly fixed already.

So I'm sure, its not a melting point problem, but a possible bug .....

LAMMPS is open source and depends on its users to help making it
better. your speculation is not very helpful and not exactly enticing
to even look into the issue, the way you are presenting it.

axel.

Hi all

so what do you expect to happen now? especially, after you publicly
state that you don’t consider it worth spending time on it yourself.

Axel, please show me where did I say that ?

more importantly, you don’t note which version and which platform this
is on. for all we know, you could be running a massively outdated
version of LAMMPS and the possible bug is possibly fixed already.

I’m using 6 Apr 2015-ICMS. The same problem was reported in Sep-2009 - http://lammps.sandia.gov/threads/msg08733.html

As I’m sure, this is not a lammps-version problem I didn’t mention the version I’m using.

LAMMPS is open source and depends on its users to help making it
better. your speculation is not very helpful and not exactly enticing
to even look into the issue, the way you are presenting it.

Only after a thorough investigation (temperature as high as 1000K and different liquids) I speculated that it could be a bug. I have posted the question only after spending enough time on it, to see if others have also experienced any similar problems. If you don’t have time, I request you not to be offensive.

Thank you very much.

Hi all

so what do you expect to happen now? especially, after you publicly

state that you don't consider it worth spending time on it yourself.

Axel, please show me where did I say that ?

"I noticed this problem few months ago, but did not investigate further."

more importantly, you don't note which version and which platform this

is on. for all we know, you could be running a massively outdated
version of LAMMPS and the possible bug is possibly fixed already.

I'm using 6 Apr 2015-ICMS. The same problem was reported in Sep-2009 -
http://lammps.sandia.gov/threads/msg08733.html
As I'm sure, this is not a lammps-version problem I didn't mention the
version I'm using.

if you don't know what causes the problem, how can you be certain,
that this is the case?
how do you know for certain, that it has not been corrected since?

don't you see that this is what provokes an aggressive reaction? on
one hand you say, that you suspect there is a problem, but don't know
what it is, and can only describe some symptoms, but on the other
hand, you are sure to know that this hasn't been addressed since and
that you are competent to judge what information matters and what not.

LAMMPS is open source and depends on its users to help making it

better. your speculation is not very helpful and not exactly enticing
to even look into the issue, the way you are presenting it.

Only after a thorough investigation (temperature as high as 1000K and
different liquids) I speculated that it could be a bug. I have posted the
question only after spending enough time on it, to see if others have also
experienced any similar problems. If you don't have time, I request you not
to be offensive.

i am not trying to be offensive, but you are not making it easy to
take your speculation seriously. i am trying to challenge you to
provide better report. i understand your frustration, but you have to
understand that your way of reporting a problem, is equally
frustrating if not more. regardless of the amount of time you spent on
it, your "investigations" are simply not very convincing, and it would
be *so* easy to provide a more consistent and convincing problem
report:

like i mentioned initially, try to map out the energy function and
forces, and compare to the provided formula,
e.g. take a system with only two or three atoms at different distances
and with different configurations, where you can compute the energy
and forces by hand and compare what should be with what LAMMPS gives
you. that will take you one day (or two if you are very careful) and
will help to confirm if there is something wrong with energies, forces
or both. then you can look at LJ and Coulomb separately as well and
narrow down whether either one or both are off.

this way, you don't spend a lot of time doing useless simulations, and
people like me don't have to spend a massive amount of time doing
trivial things that you could do as easily. in other words, if you
want to get help, you have to help others to help you.

axel.

Hi All,

Here is the bug - the coulomb forces are calculated wrongly when the distance between two atoms ® is r_in < r < r_out

Lets consider only two atoms. We set lj_epsilon equal to zero so the atoms are interacting only via coulombic interactions.

The pair potential used is - pair_style lj/charmm/coul/charmm 14 18.

As the distance between atoms systematically increased from 8 to 20 Angstroms, this is how the energy and the force acting on them are varying.

distance - energy - force (3 columns)

8 -41.507964 5.1884955
9 -36.895968 4.099552
10 -33.206371 3.3206371
11 -30.18761 2.7443282
12 -27.671976 2.305998
13 -25.543362 1.964874
14 -23.718836 1.6942026

15 -19.243479 6.7380391 Notice these three forces

16 -11.348567 8.4616273
17 -3.5826825 6.3955188
18 0 0
19 0 0
20 0 0

  1. As the distance is increased from 8 to 14 (r_in=14) both the energy and forces are decreasing as they should.

  2. For the distances 15, 16, 17 (r_in < r < r_out) the energy is decreasing as it should, but the force which should be smaller than the force when the distance is below 14, is actually increasing.

So with the lj/charmm/coul/charmm the energies are calculated correctly but the forces are wrong when the distance ® is r_in < r < r_out

The bug must be the force calculation of S®. Notice S® is also used in LJ, so the bug could be in LJ forces also for r_in < r < r_out

The script to reproduce the above data is:

units real
atom_style full
atom_modify map array

lattice sc 1.0
region box block 0 100 0 100 0 100
create_box 2 box

variable dist equal 8.0 # this dist value is increased from 8 to 20 in steps of 1

create_atoms 1 single 0 0 0
create_atoms 2 single 0 0 ${dist}
mass * 1.0

set type 1 charge 1.0
set type 2 charge -1.0

pair_style lj/charmm/coul/charmm 14 18
pair_coeff * * 0.0 1.0

variable f1 equal fz[1]
variable f2 equal fz[2]

thermo_style custom v_dist ecoul v_f1

run 0

Best,

Hi All,

Here is the bug - the coulomb forces are calculated wrongly when the
distance between two atoms (r) is r_in < r < r_out

Lets consider only two atoms. We set lj_epsilon equal to zero so the atoms
are interacting only via coulombic interactions.

The pair potential used is - pair_style lj/charmm/coul/charmm 14 18.

As the distance between atoms systematically increased from 8 to 20
Angstroms, this is how the energy and the force acting on them are varying.

distance - energy - force (3 columns)
8 -41.507964 5.1884955
9 -36.895968 4.099552
10 -33.206371 3.3206371
11 -30.18761 2.7443282
12 -27.671976 2.305998
13 -25.543362 1.964874
14 -23.718836 1.6942026
15 -19.243479 6.7380391 Notice these three forces
16 -11.348567 8.4616273
17 -3.5826825 6.3955188
18 0 0
19 0 0
20 0 0

1. As the distance is increased from 8 to 14 (r_in=14) both the energy and
forces are decreasing as they should.
2. For the distances 15, 16, 17 (r_in < r < r_out) the energy is decreasing
as it should, but the force which should be smaller than the force when the
distance is below 14, is actually increasing.

So with the lj/charmm/coul/charmm the energies are calculated correctly
but the forces are wrong when the distance (r) is r_in < r < r_out

great! now *this* is a bug report that is convincing and providing
specific information that one can work with.

...and it only took you an hour! doesn't it feel good to do a good job? :wink:

axel.

The problem is in this part of the code

if (rsq > cut_coul_innersq) {
switch1 = (cut_coulsq-rsq) * (cut_coulsq-rsq) *
(cut_coulsq + 2.0rsq - 3.0cut_coul_innersq) / denom_coul;
switch2 = 12.0*rsq * (cut_coulsq-rsq) *
(rsq-cut_coul_innersq) / denom_coul;
forcecoul *= switch1 + switch2;
}

Derivative of S® is a very different expression from what is being coded above.
http://lammps.sandia.gov/doc/pair_charmm.html

All the potentials (lj and coul) which used charmm type has the same bug.

pair_lj_charmm_coul_charmm.cpp
pair_lj_charmm_coul_charmm_implicit.cpp
pair_lj_charmm_coul_long.cpp
pair_lj_charmm_coul_msm.cpp
pair_lj_charmm_coul_msm.cpp

Best,

The problem is in this part of the code

if (rsq > cut_coul_innersq) {
            switch1 = (cut_coulsq-rsq) * (cut_coulsq-rsq) *
              (cut_coulsq + 2.0*rsq - 3.0*cut_coul_innersq) / denom_coul;
            switch2 = 12.0*rsq * (cut_coulsq-rsq) *
              (rsq-cut_coul_innersq) / denom_coul;
            forcecoul *= switch1 + switch2;
          }

Derivative of S(r) is a very different expression from what is being coded
above.
http://lammps.sandia.gov/doc/pair_charmm.html

All the potentials (lj and coul) which used charmm type has the same bug.

no, they don't. you have to look at the derivative of C(r)*S(r) and
LJ(r)*S(r), respectively.

again, you are over-eagerly making unproven speculations and jumping
to incorrect conclusions.
you have to stop doing this, or you'll sooner or later are making such
a mistake where it really matters and then discredit yourself.

the problem *only* exists for pair_lj_charmm_coul_charmm.cpp and
*only* for the coulomb part.
for coul/long and coul/msm the switching function does not apply, and
coul/charmm/implicit is correct.

this and the fact that the LJ switching is done correctly, can be
easily seen from using the pair_write command and plotting the
corresponding potential functions and comparing them to
lj/cut/coul/cut

thus the only necessary change is the patch (and the corresponding
changes to the derived classes in KOKKOS and USER-OMP). i have already
included them into LAMMPS-ICMS and will forward them to steve for
inclusion into the upstream LAMMPS version ASAP.

axel.

diff --git a/src/MOLECULE/pair_lj_charmm_coul_charmm.cpp
b/src/MOLECULE/pair_lj_charmm_coul_charmm.cpp
index 0d08a67..f12bc8f 100644
--- a/src/MOLECULE/pair_lj_charmm_coul_charmm.cpp
+++ b/src/MOLECULE/pair_lj_charmm_coul_charmm.cpp
@@ -125,9 +125,7 @@ void PairLJCharmmCoulCharmm::compute(int eflag, int vflag)
           if (rsq > cut_coul_innersq) {
             switch1 = (cut_coulsq-rsq) * (cut_coulsq-rsq) *
               (cut_coulsq + 2.0*rsq - 3.0*cut_coul_innersq) / denom_coul;
- switch2 = 12.0*rsq * (cut_coulsq-rsq) *
- (rsq-cut_coul_innersq) / denom_coul;
- forcecoul *= switch1 + switch2;
+ forcecoul *= switch1;
           }
         } else forcecoul = 0.0;

@@ -470,9 +468,7 @@ double PairLJCharmmCoulCharmm::single(int i, int
j, int itype, int jtype,
     if (rsq > cut_coul_innersq) {
       switch1 = (cut_coulsq-rsq) * (cut_coulsq-rsq) *
         (cut_coulsq + 2.0*rsq - 3.0*cut_coul_innersq) / denom_coul;
- switch2 = 12.0*rsq * (cut_coulsq-rsq) *
- (rsq-cut_coul_innersq) / denom_coul;
- forcecoul *= switch1 + switch2;
+ forcecoul *= switch1;
     }
   } else forcecoul = 0.0;

Sridhar - thanks for finding the issue with lj/charmm/coul/charmm.

I think it must have been there for a long time

(most CHARMM calculations use the coul/long variant).

We’ll post a patch …

Steve

I guess we must also thank Axel, who goaded Sridhar to delve into it :-).

Thanks Axel, Steve and all,

The new version of LAMMPS-ICMS is working fine.

Axel, so the switch2 term shouldn’t be there in the coulomb force calculations.

Thanks again for a very quick response. Let me take this opportunity to also thank Alex for developing topotools, which is a great tool for lammps-vmd community.

Best,

patch for this is being posted today.

Steve