I have a molecule of the form AB_{5}. The intermolecular potential between B atoms, depends on the intermolecular A-A distance. Thus, I need to access the index numbers of the A atoms while calculating the force between B atoms. I’ve tried:
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
factor_lj = special_lj[sbmask(j)];
factor_coul = special_coul[sbmask(j)];
j &= NEIGHMASK;
jtype = type[j];
if (jtype != 1) {
jA = atom->map(atom->tag[j] - jtype + 1);
if ((atom->molecule[jA] != atom->molecule[j])) {
moljA=atom->molecule[jA];
molj=atom->molecule[j];
fprintf(screen,“mol jA is %d, mol j is %d \n”,moljA,molj);
}
}
}
where jtype=1 for atom A, and 2, 3, 4, 5 or 6 for the different B atoms. Everything seems to work when I run on a single processor, but fails when I run on multiple processors. How can I get the correct index number of the A atom bonded to the B atom of interest?
Any help or suggestions would be greatly appreciated.
I have a molecule of the form AB_{5}. The intermolecular potential between B
atoms, depends on the intermolecular A-A distance. Thus, I need to access
the index numbers of the A atoms while calculating the force between B
atoms. I've tried:
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
factor_lj = special_lj[sbmask(j)];
factor_coul = special_coul[sbmask(j)];
j &= NEIGHMASK;
jtype = type[j];
if (jtype != 1) {
jA = atom->map(atom->tag[j] - jtype + 1);
if ((atom->molecule[jA] != atom->molecule[j])) {
moljA=atom->molecule[jA];
molj=atom->molecule[j];
fprintf(screen,"mol jA is %d, mol j is %d \n",moljA,molj);
}
}
}
where jtype=1 for atom A, and 2, 3, 4, 5 or 6 for the different B atoms.
Everything seems to work when I run on a single processor, but fails when I
run on multiple processors. How can I get the correct index number of the A
atom bonded to the B atom of interest?
i don't think that it really works correctly even for a single
processor. if so, only by chance. there are several issues with the
code as it is.
- hardcoded atom types are a *very* bad idea. have a look at the tip4p
styles how they communicate the O-type and H-type so they can have any
type number.
- you have to take properties of the neighbor list into account,
especially the treatment of "ghost" atoms vs. local atoms and the
differences between ghost atoms when newton_pair is on or when it is
off
- please also note that bonded atoms (i.e. those for which an explicit
bond was defined) will not show up in the neighborlist for non-bonded
interactions.
- the best way to address your issue is to first construct a special
neighbor list only for A type atoms with a suitably chosen cutoff.
several manybody pair styles do something like that.
that should get you started. there is probably more.