Specifying individual charge pair-interactions

Dear all,

I am trying to simulate a chain molecule (6 atoms long) and need to be able to scale the 1-4, 1-5 and 1-6 interactions.

Since special_bonds only permits up to 1-4 interactions, I can work around this limitation for the LJ part of the 1-5 and 1-6 interactions by specifying every atom in the simulation as an individual atom type, and then listing all the pair interactions via pair_coeff.

I’ve not been able to figure out how to deal with the coulomb component of the interactions however. Does anybody know if a similar thing can somehow be done for the coulomb interactions? I have not been able to figure this out.

Thank you in advance!

With best regards,

Tom

Dear all,

I am trying to simulate a chain molecule (6 atoms long) and need to be able
to scale the 1-4, 1-5 and 1-6 interactions.

Since special_bonds only permits up to 1-4 interactions, I can work around
this limitation for the LJ part of the 1-5 and 1-6 interactions by
specifying every atom in the simulation as an individual atom type, and then
listing all the pair interactions via pair_coeff.

Perhaps this is something which may not have occurred to you, but if
you have multiple chains, then keep in mind that these interactions
will apply to all pairs of atoms of the two types, not just to atoms
in the same chain.

If you need to customize the forces between atoms in the same chain,
without effecting atoms in different chains, I can post a new lammps
pair style that does this. It allows the user to specify different
set of Lennard Jones parameters for atoms in the same molecules and
different molecules, respectively. You could use this to scale down
interactions between atoms in the same chain.

I've not been able to figure out how to deal with the coulomb component of
the interactions however. Does anybody know if a similar thing can somehow
be done for the coulomb interactions? I have not been able to figure this
out.

If you are using long-range coulombics (kspace) I can't help you with
this part, I'm afraid. (I'm not sure it can be done.)

However if you are using short-range cutoffs or, implicit or debye
screened coulombics, then these interactions are calculated directly
by the code which handles that respective pair style (which handles
both electrostatic and VdW interactions in that case). This means we
can modify the coulombic interaction the same way. (That feature is
not included in my code, but you could add it.)

For implementation details, see:
http://lammps.sandia.gov/threads/msg27869.html

Let me know if you need this.
Cheers
Andrew

I'm sorry for the digression in m previous email.
Perhaps you were not interested in special inter-molecular (vs
intra-molecular) forces.

It sounds like you just want to weaken the coulombic force between a
given pair of atom types (regardless of which molecule they belong
to). To do that, I'm guessing that right now you must edit the code
or use pair_style "table", and include the coulombic forces in the
table

--Disregarding whether or not this would be a good idea--
...below is a hack which might accomplish what it sounds like you want
(without using pair_style table).

(Again, I am assuming you don't use long-range/kspace coulombics.)

(I can't guarantee this will work, but try it.)

Consider the pair_lj_charmm_coul_charmm.cpp file, for example. (I'm
using the June 17th version, but hopefully it hasn't changed since
then.)

  line# 144 of that file contains:
  fpair = (factor_coul*forcecoul + factor_lj*forcelj) * r2inv;

  line# 466 contains a similar line.
  fforce = (factor_coul*forcecoul + factor_lj*forcelj) * r2inv;

I assume you have 6 atom types, (one for each position in the chain),
and you want to weaken the coulombic interaction by 50% between atoms
of seperated by a distance of more than 4 (pairs 1,5, or 2,6 or 1,6),
then insert the following text right before BOTH of these lines (144
and 466):

if (((itype-jtype) >= 4) || ((jtype-jtype) >= 4))
  forcecoul *= 0.5;

Now, modify the code that calculates the coulombic energies.
if (evflag) ev_tally(i,j,nlocal,newton_pair,
                     evdwl,ecoul,fpair,delx,dely,delz);
Insert the following text before this line:
if (((itype-jtype) >= 4) || ((jtype-jtype) >= 4))
  ecoul *= 0.5;

Line 477 of the original code contains this line:
eng += factor_coul*phicoul;

Right before it, insert this text:
if (((itype-jtype) >= 4) || ((jtype-jtype) >= 4))
  phicoul *= 0.5;

Similar customizations can be made to other other pair styles which
use short-range coulombics, such as
  pair_lj_cut_coul_cut.cpp
  pair_lj_cut_coul_debye.cpp
  pair_lj_charmm_coul_charmm_implicit.cpp

(I haven't thought about tail-corrections. Hopefully you can ignore them)

----- INTER-vs.-INTRA molecular forces (again) ----

Now, if you only want to weaken INTRA-molecular coulombic interactions
(ie between atoms in the same molecule), then do it this way:

if ((atoms->molecule[i] == atom->molecule[j]) &&
    (((itype-jtype) >= 4) || ((jtype-jtype) >= 4)))
  forcecoul *= 0.5;

if ((atoms->molecule[i] == atom->molecule[j]) &&
    (((itype-jtype) >= 4) || ((jtype-jtype) >= 4)))
  ecoul *= 0.5;

if ((atoms->molecule[i] == atom->molecule[j]) &&
    (((itype-jtype) >= 4) || ((jtype-jtype) >= 4)))
  phicoul *= 0.5;

(I hope I counted the parenthesis correctly and don't have any mistakes.)

Again, I can't guarantee this will work (look for things I might have
overlooked). If not, hopefully this points you in the correct
direction.
Cheers.

Andrew

And let me know if you want me to post my code for customizing forces
between different molecules.
(Unfortunately it only works for only non-coulombic interactions.)

Dear Andrew,

Thank you very much for your kind replies.

Based on my limited understanding of C/C++, I wonder if I can simply apply a…

if ((atoms->molecule[i] == atom->molecule[j]) &&
(((itype-jtype) >= 4) || ((jtype-jtype) >= 4)))
forcecoul *= 0.5;

…style adjustment to pair_lj_cut_coul_cut. Hopefully it’s possible to reference atoms->molecule[i] within the pertinent function.

To clarify, I am indeed mainly interested in modifying the intra-molecular interaction, vis-^¬˚∆˙©ƒ∂-vis

Dear Andrew,

Thank you very much for your kind replies.

Based on my limited understanding of C/C++, I wonder if I can simply apply a…

if ((atoms->molecule[i] == atom->molecule[j]) &&
(((itype-jtype) >= 4) || ((jtype-jtype) >= 4)))
forcecoul *= 0.5;

…style adjustment to pair_lj_cut_coul_cut. Hopefully it’s possible to reference atoms->molecule[i] within the pertinent function.

To clarify, I am indeed mainly interested in modifying the intra-molecular interaction, vis-a-vis scaling the intra-molecular interactions whilst not effecting inter-molecular interactions. A simple coulomb cut-off will suffice in this case.

Thank you very much for your very useful suggestions. I try to implement this and let you know how I get on.

Best regards.

Tom

Dear Andrew,

Many thanks for forwarding the code and your suggestions. I very much appreciate it.

I have been working on the problem with mixed results.

If I do it by atom type, then it doesn’t crash but I don’t get the characteristics that I’m looking for in terms of diffusion coefficient, etc. To rule out the possibility that it’s the code causing this (such as a change in the order of the atom types by lammps), I tried to do it in terms of atom number, but then it crashes, and I can’t figure out why.

In terms of atom type, I have modified the pair_lj_cut_coul_cut file in a simple way with:

(On a side note, I amwondering whether turning off the coulombic
interactions makes sense for what you are trying to do. I will let
you worry about this.)

Also, make sure you remember to give each molecule a unique ID number.
If you use the "full" atom_style, then the molecule ID number is in
the 2nd column of the "Atoms" section of your LAMMPS data file I
think.

Answers to specific questions below

So I re-wrote it in terms of atomic number:

--------------------------------------------------------------
    for (jj = 0; jj < jnum; jj++) {
      j = jlist[jj];
      factor_lj = special_lj[sbmask(j)];
      factor_coul = special_coul[sbmask(j)];
      j &= NEIGHMASK;

      delx = xtmp - x[j][0];
      dely = ytmp - x[j][1];
      delz = ztmp - x[j][2];
      rsq = delx*delx + dely*dely + delz*delz;
      jtype = type[j];

      if (rsq < cutsq[itype][jtype]) {
        r2inv = 1.0/rsq;

      if (rsq < cut_coulsq[itype][jtype]) {
        forcecoul = qqrd2e * qtmp*q[j]*sqrt(r2inv); // calculate the
interaction force
      } else forcecoul = 0.0;

      if (rsq < cut_ljsq[itype][jtype]) {
        r6inv = r2inv*r2inv*r2inv;
        forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
      } else forcelj = 0.0;

      // Adjust the internal parameters for the 6-site molecule
      if (atom->molecule[i] == atom->molecule[j]) {
        if ((i-j == 1) || (j-i == 1)) { // 1-2 interaction
            forcecoul *= 0.0;
            forcelj *= 0.0;
        }
        else if ((i-j == 2) || (j-i == 2)) { // 1-3 interaction
            forcecoul *= 0.0;
            forcelj *= 0.0;
        }
        else if ((i-j == 3) || (j-i == 3)) { // 1-4 interaction
            forcecoul *= 0.0;
            forcelj *= 0.0;
        }
        else if ((i-j == 4) || (j-i == 4)) { // 1-5 interaction
            forcecoul *= 0.0;
        }
        else if ((i-j == 5) || (j-i == 5)) { // 1-6 interaction
            forcecoul *= 0.0;
        }
      }

--------------------------------------------------------------
This compiles, but then crashes...

Do you (or anyone else who's read this far) know what I'm doing wrong for the
by-atomic-number case? Does the by-type method look ok?
Much appreciating your help.

As for why it crashes, it would help to find out which line it's
crashing on using a debugger (see gdb instructions below). But
perhaps it does not matter (see below).

  Unfortunately the variables "i" and "j" in your code above can not
be trusted...
What I mean is: the variables "i" and "j" in your code will not be
identical to the numbers in the first column of the "Atoms" section of
your LAMMPS DATA file.) LAMMPS has its own internal numbering system.
For this reason I would not try this approach (even if you can find a
way to keep it from crashing).

  Are the itype and jtype variables trustworthy?
I hope I did not lead you astray.
I am not sure if itype and jtype variables correspond to the atom type
numbers in your LAMMPS data file (usually the 3rd column of the
"Atoms" section). I THINK they do, but you will have to test this.
It's possible that these numbers are shifted by 1. (IE itype=0
corresponds to "1" in the DATA file.)

You have to check the identity of each atom and make sure you are
tuning off interactions between the correct pairs of atoms.

You can do this using printf statements (see suggestions below),
or GDB (instructions on gdb below).

One way to identify each aom is by looking at the coordinates.
Although they will have shifted, the relative position between each
atom should be the same, so you can tell which atoms are at the
endpoints by adding some printf statements to your code which print
out x[i][0],x[i][1],x[i][2] and x[j][0],x[j][1],x[j][2]. For details
see below.

Dear Andrew,

Many thanks for forwarding the code and your suggestions. I very much
appreciate it.

I have been working on the problem with mixed results.

If I do it by atom type, then it doesn't crash but I don't get the
characteristics that I'm looking for in terms of diffusion coefficient, etc.
To rule out the possibility that it's the code causing this (such as a
change in the order of the atom types by lammps), I tried to do it in terms
of atom number, but then it crashes, and I can't figure out why.

In terms of atom type, I have modified the pair_lj_cut_coul_cut file in a
simple way with:

--------------------------------------------------------------
    for (jj = 0; jj < jnum; jj++) {
      j = jlist[jj];
      factor_lj = special_lj[sbmask(j)];
      factor_coul = special_coul[sbmask(j)];
      j &= NEIGHMASK;

      delx = xtmp - x[j][0];
      dely = ytmp - x[j][1];
      delz = ztmp - x[j][2];
      rsq = delx*delx + dely*dely + delz*delz;
      jtype = type[j];

      if (rsq < cutsq[itype][jtype]) {
        r2inv = 1.0/rsq;

        if (rsq < cut_coulsq[itype][jtype]) {
           forcecoul = qqrd2e * qtmp*q[j]*sqrt(r2inv); // calculate the
interaction force
        } else forcecoul = 0.0;

        if (rsq < cut_ljsq[itype][jtype]) {
          r6inv = r2inv*r2inv*r2inv;
          forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
        } else forcelj = 0.0;

        // Adjust the internal parameters for the 6-site molecule
          if (atom->molecule[i] == atom->molecule[j]) {
              if (itype == 1) { // *** A
                  if (jtype == 1) { // A
                      forcecoul *= 0.0;
                  }
                  else if (jtype == 2) { // B
                      forcecoul *= 0.0;
                  }
                  else if (jtype == 3) { // C
                      forcecoul *= 0.0;
                  }
              }
              if (itype == 2) { // *** B
                  if (jtype == 1) { // A
                      forcecoul *= 0.0;
                  }
                  else if (jtype == 2) { // B
                      forcecoul *= 0.0;
                      forcelj *= 0.0;
                  }
                  else if (jtype == 3) { // C
                      forcecoul *= 0.0;
                      forcelj *= 0.0;
                  }
              }
              if (itype == 3) { // *** C
                  if (jtype == 1) { // A
                      forcecoul *= 0.0;
                  }
                  else if (jtype == 2) { // B
                      forcecoul *= 0.0;
                      forcelj *= 0.0;
                  }
                  else if (jtype == 3) { // C
                      forcecoul *= 0.0;
                      forcelj *= 0.0;
                  }

Try adding some printf (or fprintf) statements here in your code:

fprintf(stderr,"itypex[%d]=%d, (x,y,z)=(%g,%g,%g)\n", i,
itype,x[i][0],x[i][1],x[i][2]);
fprintf(stderr,"jtypex[%d]=%d, (x,y,z)=(%g,%g,%g)\n", j,
jtype,x[j][0],x[j][1],x[j][2]);
fprintf(stderr,"fourcelj=%g, fourcecoul=%g\n", fourcelj, fourcecoul);

              }

This way you can check to see that you are turning off the forces
between the correct pairs of atoms. If you put them in the correct
place then they are only printed between atoms in the same molecule.

Perhaps there is a better way to do this, but this should work.

--------------------------------------------------------------

(the testing 6-site molecules have the pattern A-B-C-C-B-A)

Okay. Keep in mind that each time you turn off interactions between a
pair of atom types, such as A and B, you are actually turning off
interactions between four pairs of atoms in each chain: (1st,2nd),
(1st,5th), (2nd, 6th), and (2nd,6th) atoms. I hope this is what you
want. The printf statements above should help figure out what is
hapenning.

Cheers
Andrew

Makefile.ubuntu_parallel.txt (649 Bytes)

Makefile.ubuntu_dbg (657 Bytes)

Regarding the modification of lj/cut/coul/cut interaction:

Andrew, thank you again for your kind email. It has been incredibly helpful.

I didn’t follow the debugging route for the current situation (although this knowledge in general is very useful), because I found that the by-molecule detection of intra-molecular interaction was suitable. Through testing, I found that only the 1-5 and 1-6 interactions (ie, A-B and A-A interactions in the A-B-C-C-B-A molecule) were being considered by that particular point in the program. So it’s possible to deal with these interactions here, and use special_bonds to deal with 1-4 and shorter interactions.

So for this aspect, I think it is working!!

Regarding the modification of lj/cut/coul/long interaction:

It turns out that I do need to include ewald interaction (I was mistaken earlier when I said it wasn’t necessary), so I have tried to implement this in a similar way into the lj/cut/coul/long pair interaction style. In order to check my understanding, If anyone could confirm the following things for me certainly I’d appreciate it:

  1. Looking at the code for lj/cut/coul/long, it seems that the coulomb force is zero whenever “rsq < cut_coulsq”, which means that a given molecule only interacts with molecules (and their pbc ewald images) that are within this coulomb cut-off distance.

  2. If I set the force to zero when dealing with atoms from the same molecule, as described in the previous emails, this will ensure there is no 1-5 or 1-6 interaction of this atom with any of the other atoms in the molecule, whether short-range or through ewald.

  3. I can take out the short-range coulomb force whilst preserving the ewald image interaction by adding the following lines within the “if (rsq < cut_coulsq)” function:

// Adjust the internal parameters for 1-5 and 1-6 interactions
if (atom->molecule[i] == atom->molecule[j]) {
if (itype == 1) { // *** A
if (jtype == 1) { // - A
forcecoul -= qqrd2e * qtmp*q[j]*sqrt(r2inv);
if (factor_coul < 1.0) forcecoul += (1.0-factor_coul)*forcecoul;

}
else if (jtype == 2) { // - B
forcecoul -= qqrd2e * qtmp*q[j]*sqrt(r2inv);
if (factor_coul < 1.0) forcecoul += (1.0-factor_coul)forcecoul;
}
}
if (itype == 2) { // *** B
if (jtype == 1) { // - A
forcecoul -= qqrd2e * qtmp
q[j]*sqrt(r2inv);
if (factor_coul < 1.0) forcecoul += (1.0-factor_coul)*forcecoul;
}
}
}

(note the reversal of signs “-=” and “+=” in attempt to cancel-out the short-range force)

If anyone is able to highlight a flaw in the above reasoning / methodology, I would be very grateful to learn about it.

With sincere thanks,

Tom