problems with Ewald sums and rigid body dynamics

I can receive attachments on the mail list. Please keep the
list in the loop.

If you remove the neigh_modify exclude molecule command
from your scripts then you get nearly identical energies
independent of cutoff or accuracy as you should. The dynamics
should also be correct, since forces within the rigid bodies will
cancel.

I'm not clear on whether your cutoff values are consistent with
excluding the short-range terms within the rigid bodies. How
large are the rigid bodies?

Paul can probably comment on whether it is legitimate to
exclude intra-molecule terms in this manner with the long-range
solver. Paul - if you look at these scripts as-is, the energy
changes dramatically with cutoff and precision. Should this
happen if the excluded interactions are wholly within the cutoff?

Steve

From the manual: "the exclude options only affect pairwise interactions". This means that you can't use the exclude options with kspace and expect to get correct answers, as you've re-discovered. The reason is that the kspace computes do not exclude charged particle interactions, therefore requiring the complete neighbor list to turn off coulombic interactions between bonded particles. Probably can't do the correct thing with your rigid body, long-range interaction, charged system in LAMMPS as it currently is. May be a clever way around this limitation.

Paul

I will emphasize however, that you can choose not to exclude the
Coulombic interactions
within the rigid bodies. Since they sum to 0.0 for each body, that
should not affect
the dynamics. It will simply give you a large energy value for each
rigid body. That
value should be constant so you could just subtract it off to get the
Kcal/mole value
you are expecting. I.e. you could pre-compute it. Or use the rerun
command on your
snapshots to compute that term.

Steve

Dear Steve
thanks a lot for your quick response,

with a complex workaround in the input file we have been able to exclude the intra-molecular lennard-jones so that we can avoid the card “neigh_modify exclude”, as you said everything is correct in this case and cutoff independent.

if anyhow we do not exclude lj, we sum up so many huge positive numbers that we end-up in a numerical error bigger than the
quantities that we need to compute.

Thus it is fundamental for us to find a way to

-keep the intramolecular coulomb (we can the easily subtract it afterwards as you said and it is small)
-BUT we have to avoid the intramolecular LJ

is there an easy way to do this?

Thanks in advance for any further suggestions

Andrea

up to now what we did is:

double each molecule, defining two additional types
and excluding interactions among the types that are lj in the same molecule:

units real

atom_style full

boundary p p p

read_data doppie.dat

mass 1 12.011
mass 2 0.0001
mass 3 0.0001
mass 4 0.0001
group lelj type 3 4

velocity all create 50.0 87287

pair_style lj/cut/coul/long 30.0 14.0
pair_coeff 1 1 0.0 1.0
pair_coeff 1 2 0.0 1.0
pair_coeff 1 3 0.0 1.0
pair_coeff 1 4 0.0 1.0
pair_coeff 2 2 0.0 1.0
pair_coeff 2 3 0.0 1.0
pair_coeff 2 4 0.0 1.0
pair_coeff 3 3 0.0298 3.4
pair_coeff 3 4 0.0298 3.5
pair_coeff 4 4 0.0298 3.6
kspace_style ewald 1.0e-10

neighbor 0.3 bin
neigh_modify every 20 delay 0 check no one 100000 page 1000000
neigh_modify exclude molecule lelj

fix 1 all rigid molecule langevin 50.0 50.0 1000.0 1948

dump 1 all xyz 500 pos_50_28.04.xyz

thermo 100
thermo_style custom step temp ke pe etotal
thermo_modify format 5 %22.14g

run 100000

I'm not sure why you defined several extra atom types (3 and 4), or
why you doubled the number of molecules.
(Perhaps I did not read the description of your system carefully enough.)
Have you tried something simple like this? ...

    group lelj molecule <> 50 250
    neigh_modify exclude molecule lelj

(See http://lammps.sandia.gov/doc/group.html)

Other than neighmodify exclude, there is a new pair style
"lj/charmm/coul/charmm/inter" which allows you to specify different
lennard-jones parameters for inter and intra-molecular LJ
interactions. You could use this to switch off intramolecular LJ.
However I think this would have the same effect as using neighmodify
exclude.

In any case, the documentation for that new pair style is located here:
http://www.moltemplate.org/lammps_code/pair_lj_charmm_coul_charmm_inter.html

You can download the code for this pair style here:
http://www.moltemplate.org/lammps_code/index.html

Hasan Babaei and I have independently tested this code recently and it works.
(This pair style is used frequently in the "moltemplate" examples
online and I will continue to support and update it. It's my fault it
hasn't yet been added to LAMMPS. I've been too busy convert the
documentation to lammps format. I will do this after grant season.)

Cheers
Andrew

is there an easy way to do this?

no - I don't think there currently is.
There would have to be logic added to the
neighbor list to encode the excluded intra-molecular
pairs so that their Coulomb interaction could
be handled correctly in the pairwise interactions
when used with a KSpace solver. Similar to
what is done for the special bonds 1-2,1-3,1-4
interactions.

Steve

I think that pair_style I posted in the previous email
(lj/charmm/coul/charmm/inter) does what you asked for. Unfortunately
this pair-style calculates short-range electrostatics incorrectly (not
long-range). So you have to tell it NOT to compute electrostatic
interactions (by setting the coulombic-cut-offs to zero) and then
combine it with pair_style coul/long using the "hybrid/overlay"
keyword. Download the code at
http://www.moltemplate.org/lammps_code/
recompile LAMMPS, and try this:

pair_style hybrid/overlay &
lj/charmm/coul/charmm/inter 29.5 30 0 0 &
coul/long 14.0

pair_coeff lj/charmm/coul/charmm/inter 1 1 0 0 0 0 0.0298 3.4
pair_coeff lj/charmm/coul/charmm/inter 1 2 0 0 0 0 0.0298 3.5
pair_coeff lj/charmm/coul/charmm/inter 2 2 0 0 0 0 0.0298 3.6

pair_coeff coul/long * *

# (The "pair_style" command was supposed to fit on a single line,
# but I was worried it might get cut into pieces during email.
# Fortunately LAMMPS allows you to use the "&" character
# to split one command into multiple lines.)
# Note that lj/charmm/coul/charmm/inter ALSO computes the
# coulombic interaction, so you have to disable that by setting
# the coulombic cutoffs to zero. (That's what the "0 0" does
# in the middle line.) This should not effect the coulombic
# interactions calculated by pair_style coul/long.

Hi everybody,
I tried to use the “lj/charmm/coul/charmm/inter” as suggested by Andrew, this is the input

units real
atom_style full

boundary p p p

read_data initpos.dat

mass 1 12.011
mass 2 0.0001

velocity all create 0.00000 87287

pair_style hybrid/overlay lj/charmm/coul/charmm/inter 29.5 30 0 0 coul/long 14.0
pair_coeff lj/charmm/coul/charmm/inter 1 1 0 0 0 0 0.0298 3.4
pair_coeff lj/charmm/coul/charmm/inter 1 2 0 0 0 0 0.0298 3.5
pair_coeff lj/charmm/coul/charmm/inter 2 2 0 0 0 0 0.0298 3.6
pair_coeff coul/long * *

kspace_style ewald 1.0e-10

neighbor 0.3 bin
neigh_modify every 20 delay 0 check no one 20000 page 200000

fix 1 all rigid molecule

thermo 1
thermo_style custom step temp ke pe etotal
thermo_modify format 5 %22.14g
run 1

and this is the output:

LAMMPS (13 Oct 2012)
Scanning data file …
Reading data file …
orthogonal box = (0 0 0) to (28.0816 28.0816 28.0816)
2 by 2 by 2 MPI processor grid
2880 atoms
Finding 1-2 1-3 1-4 neighbors …
0 = max # of 1-2 neighbors
0 = max # of 1-3 neighbors
0 = max # of 1-4 neighbors
1 = max # of special neighbors
ERROR: Numeric index is out of bounds (force.cpp:666)

In a normal run, after calculating the number of neighbors, the codes usually displays the number of rigid bodies and then it initializes the Ewald mesh. So it seems that this error occurs prior to the Ewald setup.

Actually I can avoid to make NPT dynamics but to study the rotational melting transition I will need to make NVT calculations.

Thanks for your help.
Andrea

Some atom type index, probably in your input script, is invalid
I would guess.

Steve

I left the same input coordinates and the same input file, just replacing the pair_coeff and the pair_style

pair_style hybrid/overlay &
lj/charmm/coul/charmm/inter 29.5 30 0 0 &
coul/long 14.0
pair_coeff lj/charmm/coul/charmm/inter 1 1 0 0 0 0 0.0298 3.4
pair_coeff lj/charmm/coul/charmm/inter 1 2 0 0 0 0 0.0298 3.5
pair_coeff lj/charmm/coul/charmm/inter 2 2 0 0 0 0 0.0298 3.6
pair_coeff coul/long * *

with the old ones:

pair_style lj/cut/coul/long 14.0 30.0
pair_coeff 1 1 0.0298 3.4
pair_coeff 1 2 0.0298 3.5
pair_coeff 2 2 0.0298 3.6

and the code works, maybe is a problem with the "pair_coeff coul/long * * " instruction?

I left the same input coordinates and the same input file, just replacing
the pair_coeff and the pair_style

pair_style hybrid/overlay &
lj/charmm/coul/charmm/inter 29.5 30 0 0 &

coul/long 14.0

Sorry I was away from email the past couple of days.
1) First of all, this is supposed to be a single command. (I was using
"&" characters to split the line up into multiple lines, but they have
to be contiguous.) Try getting rid of this blank line and see if it
works.

In other words, get rid of the blank line between:

lj/charmm/coul/charmm/inter 29.5 30 0 0 &

and

coul/long 14.0

pair_coeff lj/charmm/coul/charmm/inter 1 1 0 0 0 0 0.0298 3.4
pair_coeff lj/charmm/coul/charmm/inter 1 2 0 0 0 0 0.0298 3.5
pair_coeff lj/charmm/coul/charmm/inter 2 2 0 0 0 0 0.0298 3.6
pair_coeff coul/long * *

with the old ones:

pair_style lj/cut/coul/long 14.0 30.0
pair_coeff 1 1 0.0298 3.4
pair_coeff 1 2 0.0298 3.5
pair_coeff 2 2 0.0298 3.6
and the code works, maybe is a problem with the "pair_coeff coul/long * * "
instruction?

It could be. I've actually never tried using pair_coeff coul/long before.

What do the last few lines of your log.lammps file say?

You are bug testing my code for me (thanks). I got the idea for
doing it that way from an example at this web page.

2) Also, if the suggestions above do not work, then try getting rid of
the wildcard characters * *, and see if it helps. In other words, try
this:

Hi Andrew,
I tried your suggestions but I have the same problem:

  1. the blank line was just a result of the cut-and-paste of the input file into the mail (see input1 and output1)

  2. I also tried to substitute the * * with the explicit three combinations but the result is the same.

  3. this combination works only if I switch off the kspace calculation (see input3_bis and output3_bis) naturally the energy is not the correct one. If I ask LAMMPS to calculate also long range coulombic contribution an error message appear saying that the Kspace command is incompatible with the pair_style (see input3 and output3).

I also attached the coordinate file inipos.dat.
thanks for your help.
Andrea

initpos.dat (220 KB)

input1 (697 Bytes)

input3 (563 Bytes)

input3 (563 Bytes)

output1 (376 Bytes)

output_3 (448 Bytes)

output3_bis (1.21 KB)

Hi Andrea,

please see the example of pair_coeff in the documentation: the atom
types should go first, then the sub-style, and finally, the arguments:

http://lammps.sandia.gov/doc/pair_hybrid.html

The syntax for pair_coeff with pair style hybrid/overlay is supposed to be

pair_coeff 1 1 lj/charm/coul/charmm/inter arg1 arg2 ...
pair_coeff 1 2 lj/charm/coul/charmm/inter arg1 arg2 ...
pair_coeff 2 2 lj/charm/coul/charmm/inter arg1 arg2 ...
pair_coeff * * coul/long

Hope that helps,
-Trung

I'm sorry, that was definitely my fault. She was following my
instructions (which were wrong).
Thanks Trung!

In addition, there was also one other problem. (I was using "0 0" as
the inner and outer cutoff for the coulombic interactions.
Apparently, the coulombic cutoffs must be different.) So I had to
replace "0 0" with something like "0 0.000000001". This seems to work
fine. (See attached file "input1_000000001" Excerpt below:)

pair_lj_charmm_coul_charmm_inter.cpp (36.4 KB)

input1_corrected (705 Bytes)

input1_000000001 (715 Bytes)

Dear all,
following the last Andrew and Trung instructions the numerical results for the Fullerite cohesive energy seem to be correct.
Andrew, soon or later I will also test the NPT dynamics with the payr_style “lj/charmm/coul/charmm/inter” and I will refer to you for any possible bug report.
Thanks a lot for your quick help and suggestions.
Cheers. Andrea