(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)