excluded intramolecular interactions in kspace

LAMMPS LAMMPS Mailing List Mirror
1

Dear LAMMPS community,

I am interested in how LAMMPS handles the excluded intramolecular interactions (e.g. atoms in a dihedral) when using kspace (e.g. Ewald). The efficient Ewald implementation includes spurious self-self and excluded interaction terms, which should be corrected with terms of the kind qi*qj*alpha/sqrt(pi) for self-self and qi*qj*erf(alpha*rij)/rij for excluded intramolecular interactions.

Could someone please point me to the appropriate source code file where these corrections are implemented? I cannot seem to find them in the pair, kspace or bonded (e.g. dihedral) files.

Also, on theoretical grounds, I don't understand how the 1-4 excluded terms in dihedral_charmm.cpp (e.g. line ~259) may be included with a simple coulombic pair-wise term if one is utilizing kspace (e.g. Ewald). If kspace is used, this 1-4 term should be the real-space erfc() term instead. Utilizing the pair-wise Coulombic form leads me to believe that the excluded terms were handled somehow differently than the traditional method that I understand (e.g. as well explained in the DL_POLY2 Manual). If the coulomb form is used, then one would have to correct for the fact that one assumed the system is embedded in a medium of infinite dielectric constant (leading to the so called 'polarization' term of the square of the dipole moment of the system). But there is no indication that the charges on the 1-4 interactions in kspace routines are excluded (which would be nontrivial, if not impossible).

Best Regards,
Harold Hatch

2

Hello, Harold:

The self-self terms are included in ewald.cpp under the compute routine, while the issues with 1-4 (and 1-2 and 1-3) pairs are dealt with using the special_bonds command.

Those should answer most of your questions.

—AEI

3

The simplest bit of code that does this is
in pair_lj_cut_coul_long.cpp, Similar code is
in any pair style with a "long". Here it is:

      prefactor = qqrd2e * qtmp*q[j]/r;
      forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
      if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;

Factor_coul is a value from 0 to 1 (inclusive), which is setup by the
special_bonds command and is peculiar to this pair of atoms. So it is
being applied only to the closest pair of all the images of these 2 atoms.

The long-range solver (Ewald, PPPM) thus needs to do nothing
special. There is also no self-self interaction in LAMMPS, since
there is no I-I entry in the neighbor list.

Steve

4

Thank you Steve and Ahmed for your helpful comments!

I am interested in how excluded intramolecular interactions (e.g. 1-2,1-3 dihedrals) are handled in kspace by LAMMPS, and would like to know where exactly this is implemented in the source code. I understand how they are handled in real space by weights and/or removal from neighbor lists. But in kspace, it is either impossible or terribly inefficient to remove these interactions. Thus, all charge-charge interactions are calculated in kspace, and the spurious intramolecular interactions are removed in real space with a term like the following

-qi * qj * erf(g_ewald * |rij|) / (|rij|)

for each excluded intramolecular i,j particle pair. Note this is the error function, erf, and NOT the complementary erfc which is used in the non-bonded real space term. This term is not as well-known as the self term, and is not explained in Frenkel and Smit or Allen and Tildesley for simplicity. DL_POLY has it in their manual: http://wanglab.bu.edu/DLPOLY2/node65.html , eq. 2.131, last term.

Regarding Ahmed's reply, I have indeed found the self terms in ewald.cpp, about line 285: ( energy -= g_ewald*qsqsum/1.772453851 +). This is the well known correction due to the spurious interaction of a point charge with the Gaussian cloud centered about the same particle, and is different from the excluded intramolecular terms.

Best Regards,
Harold Hatch

Ahmed wrote "The self-self terms are included in ewald.cpp under the compute routine, while the issues with 1-4 (and 1-2 and 1-3) pairs are dealt with using the special_bonds command."

5

I'm not clear if you are clarifying how LAMMPS does this,
or saying there is something missing in the LAMMPS formulation.

Steve

6

Steve:

The argument is that there is a term missing, comparable to the first term in Eq. 2.5 of Essmann et al. (1995).

—AEI

7

There is a >99% chance the terms I expect are either in the LAMMPS code somewhere that I can’t find, or someone really smart figured out a trick to include these terms in a different way. So really I am asking for an explanation of the source code, which is asking for a lot. I need to collect these terms to include them in a custom calculation I am doing with LAMMPS, and I am hoping someone familiar with the intricacies of Ewald methods in LAMMPS would point me in the right direction.

There are two reasons I’m confident LAMMPS isn’t missing these terms, aside from LAMMPS being extremely well tested by others

  1. I think water with shake requires these excluded intramolecular terms, and LAMMPS has been well tested for this case. I compared my own water codes to LAMMPS and got the same answer. Although if you use rigid body formalism you don’t need to worry about these exclusion terms.
  2. I compared a protein simulation between LAMMPS and NAMD, so both LAMMPS and NAMD would have to be wrong if LAMMPS didn’t have these terms.

By the way, LAMMPS is a beautiful code and is my favorite MD package! I apologize for insinuating there is something wrong with it, just because I can’t find out how it was done. I’ve had people use my codes and tell me there are bugs, but they simply didn’t understand the code :slight_smile:

Thank you,
Harold Hatch

8

thanks, but we could easily have a bug or missing term somewhere,
especially if it's a small contribution.

Maybe Paul can comment on this thread ...

Steve

9

LAMMPS is not missing a term. We just haven't shown where the term is in the source code to Harold's satisfaction. I'll give it a shot.

First, we need to be clear about what term we're talking about. We're talking about the term that Harold refers to in this paragraph:

> all charge-charge interactions are
> calculated in kspace, and the spurious intramolecular interactions
> are removed in real space with a term like the following
>
> -qi * qj * erf(g_ewald * |rij|) / (|rij|)
>
> for each excluded intramolecular i,j particle pair. Note this is
> the error function, erf, and NOT the complementary erfc which is
> used in the non-bonded real space term. This term is not as
> well-known as the self term, and is not explained in Frenkel and
> Smit or Allen and Tildesley for simplicity. DL_POLY has it in
> their manual:
> http://wanglab.bu.edu/DLPOLY2/node65.html , eq. 2.131, last term.

Let's call it the "intramolecular self term," which is not to be confused with the regular "self term": qi*qj*alpha/sqrt(pi). LAMMPS correctly includes both of these.

Next, note that we've so far been talking only about energies. These terms take a different form when converted to forces. The regular "self term" is independent of particle positions, so it does not contribute to the force. Ahmed and Harold have already identified the "self term" energy contribution:

> Regarding Ahmed's reply, I have indeed found the self terms in
> ewald.cpp, about line 285: ( energy -= g_ewald*qsqsum/1.772453851
> +). This is the well known correction due to the spurious
> interaction of a point charge with the Gaussian cloud centered
> about the same particle, and is different from the excluded
> intramolecular terms.

But the "intramolecular self term" contributes to both the energy and the force. We've already agreed upon the correct expression for the energy (-qi * qj * erf(g_ewald * |rij|) / (|rij|)), but no one in this thread has given the correct expression for the force yet. So, we just need to come up with the correct expression for the force (left as an exercise for the reader), point to the code that includes the energy and force versions of the term, and convince everyone that the energy and force versions of the term in LAMMPS are the same as what we've identified the energy and force versions of the term should be (left as exercises for the reader). Let's start with the LAMMMPS code for the contribution to the force. Steve has already pointed us to the code that does this:

>> The simplest bit of code that does this is in
>> pair_lj_cut_coul_long.cpp, Similar code is in any pair style with
>> a "long". Here it is:
>>
>> prefactor = qqrd2e * qtmp*q[j]/r;
>> forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
>> if (factor_coul< 1.0) forcecoul -=
>> (1.0-factor_coul)*prefactor;
>>
>> Factor_coul is a value from 0 to 1 (inclusive), which is setup by
>> the special_bonds command and is peculiar to this pair of atoms.
>> So it is being applied only to the closest pair of all the images
>> of these 2 atoms.

The LAMMPS contribution to the energy can be found in the same function about 40 lines further down (line 186 in pari_lj_cut_coul_long.cpp in the version I'm looking at):

      if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;

But this doesn't look the same as "intramolecular self term" contribution to the energy that we've agreed upon above. Let's assume that factor_coul is zero, meaning that we want to turn off the interaction. Stepping through some algebra, we get:

prefactor = qqrd2e * qtmp*q[j]/r;
ecoul -= (1.0-factor_coul)*prefactor;
ecoul -= qqrd2e * qtmp*q[j]/r;

which tells us we've subtracted off the entire interaction, including the kspace and real space contributions. But this is OK since LAMMPS has already added both the real and kspace contributions to this term. Both the real space and kspace contributions were added in the same way as all of the other non-cancelled terms. So the form of the equation that we identified for the "intramolecular self term" energy never shows up, but in the end, the math is done correctly in LAMMPS.

In the interest of time, I'll leave the corresponding analysis for the force term as an exercise for the reader.

Paul

10

Paul, thank you for the explanation of this elegant implementation!