Error encountered in parallel run of modified tip4p pair type

LAMMPS LAMMPS Mailing List Mirror
1

Dear LAMMPS developers,

We modified the pair style lj/cut/tip4p/long, such that the virtue site (instead of oxygen) in the TIP4P model will be used for computing Lennard-Jones (LJ) interaction.
Two files were modified, pair_lj_cut_tip4p_long.cpp and pair_lj_cut_tip4p_long.h, in the “settings” function so that it reads in additional parameters.

nwtype = force->inumeric(FLERR,arg[7]);

if (nwtype < 1) error->all(FLERR,"tip4p_yw requires at least one wall atom type\n”);

if (narg < 8+nwtype) error->all(FLERR,"not enough parameters for pair tip4p_yw, check number wall types specified\n”);

typeW = new int[nwtype];
for(int i=0; i<nwtype; i++)
{
typeW[i] = force->inumeric(FLERR,arg[8+i]);
}

To compute LJ interaction using TIP4P virtual site (instead of oxygen), we modified the “compute” function. The modified codes are attached.

The modified code works fine with serial runs. But it crashed with parallel runs with the error message:
“ERROR on proc 1: TIP4P hydrogen is missing (…/pair_lj_cut_tip4p_long.cpp:201)”.
And when we checked the source code, this error message will be printed if the two atoms following an oxygen atom are not hydrogen atoms.
This error does not happen when the code runs in serial version. This means that in parallel runs, the atom types were changed, or not communicated correctly among cores.

Could someone guide me the directions to debug this error?

More explanation about using this modified pair type is as follows:

pair_style lj/cut/tip4p/long otype htype btype atype qhist ljcut qcut nwtype wtype1 wtype2 … wtypenw

Syntax
otype, htype = atom types for TIP4P O and H
bytpe, atype = bond and angle types for TIP4P waters
qdist = distance from O atom to the massless charge M (distance units)
ljcut = global cutoff for the LJ interaction (distance units)
qcut = global cutoff for the Coulombic interaction (distance units)
nwtype = number of types of atoms that interact with M, not O, in computing the LJ interaction
wtype1, wtype2, …, wtypenw = types of atoms that interact with M, not O, in computing the LJ interaction

Example
Pair_style lj/cut/tip4pwwa/long 3 4 1 1 0.15 14.0 14.0 2 1 2

Thank you in advance.

Yanbin

pair_lj_cut_tip4p_long.cpp (20.8 KB)

pair_lj_cut_tip4p_long.h (3.35 KB)

2

Dear LAMMPS developers,

We modified the pair style lj/cut/tip4p/long, such that the virtue site (instead of oxygen) in the TIP4P model will be used for computing Lennard-Jones (LJ) interaction.

This doesn’t make sense. Why can’t you treat those atoms as regular and use an adjusted geometry?

Axel

3

Hi, Axel,

Thanks for your comments.

Regarding your question “Why can’t you treat those atoms as regular and use an adjusted geometry?”, it is because only LJ interaction between water and selected atoms uses the virtue site (instead of oxygen).

For example, we have a system composed of water, ion and hexagonal boron-nitride surface.
In computing water-boron and water-nitrogen LJ interactions, virtue site (instead of oxygen) is used.
In computing water-water and water-ion LJ interaction, oxygen position is used. We would like to keep the water-water interaction intact as provided in the TIP4P water model.

Please let us know if there is anything unclear.
We appreciate if you could guide us the directions to debug the error.

Thank you.

Yanbin

4

Re: debugging, if it works in serial but not parallel, then

I suggest you find the simplest 2-proc problem that

fails reliably. Fewest # of atoms, fewest # of timesteps.

Then print out the info that is triggering the error. If

some values are bogus, e.g. atom IDs or types leading

to this error, then you have a clue what to track during

the timesteps and communicatino leading up to

that step. I also highly recommend that you first run

valgrind (on each proc) to see if there are memory

errors in both the serial and parallel cases. That often

finds issues you never imagined were there.

Steve

5

Dear Dr. Plimpton,

Thank you for your guidance.

We figured that the error happens when we try to compute the virtual site for all the oxygen atoms in the neighbor list. Here is the demonstration:
Consider a system composed of boron, nitrogen, and water. In computing water-boron and water-nitrogen LJ interaction, virtual site (instead of oxygen) will be used.
In our previous code, we compute the virtual sites for all oxygen atoms in the neighbor list of boron/nitrogen. But for some oxygen atoms in the neighbor list, its hydrogen atoms may not be in the list. Error happens when we try to access those hydrogen atoms, with the error message “ERROR on proc 1: TIP4P hydrogen is missing (…/pair_lj_cut_tip4p_long.cpp:201)”.
This bug is fixed by only computing virtual site when the distance between boron/nitrogen atom and oxygen is less than cut_coulplus, where cut_coulplus=cut_coul + 2.0*q_dist and qdist=0.15 angstrom.

The bug was fixed a couple of months ago and posted at
https://github.com/wuyb02/ff-bn-w/tree/master/lammps
But I figured that I shall conclude this post with a solution.

Sincerely,
Yanbin Wu