I’m using ‘special_bonds lj/coul’ to exclude the 1-2,1-3 and 1-4 interactions, but just find that no matter what three values I put after special_bonds, e.g. special_bonds lj/coul 0.0 0.0 0.0 or special_bonds lj/coul 1.0 0.0 0.0, the long-range electrostatic energy calculated in K-space does not change. Does that mean the contributions from 1-2, 1-3 and 1-4 interactions are not excluded from long-range electrostatic energy? What command should I use if I want to exclude the corresponding 1-2, 1-3, 1-4 contributions from long-range energy?

Dear all,
I'm using 'special_bonds lj/coul' to exclude the 1-2,1-3 and 1-4
interactions, but just find that no matter what three values I put after
special_bonds, e.g. special_bonds lj/coul 0.0 0.0 0.0 or special_bonds
lj/coul 1.0 0.0 0.0, the long-range electrostatic energy calculated in
K-space does not change. Does that mean the contributions from 1-2, 1-3 and
1-4 interactions are not excluded from long-range electrostatic energy? What

exactly.

command should I use if I want to exclude the corresponding 1-2, 1-3, 1-4
contributions from long-range energy?

why should that be desirable? and how should it be possible?

Dear Steve and Axel,
Thanks for your reply. I use both GROMACS and LAMMPS to simulate CH3OCH2CH2OCH3 in which 1-2, 1-3, and 1-4 non-bonded interactions are excluded. In GROMACS, if any of these non-bonded interactions are excluded, the coulombic energy in reciprocal space changes correspondingly. This is reasonable because if a pair of interaction is excluded, its contributions in both short-range and long-range should be excluded too. But LAMMPS just excludes the short-range energy if the corresponding value in special_bonds is set to zero. Axel replies that ‘why the exclusion from long-range energy is desirable’. My interpretation is that if not excluded, we are having a pair of interaction which equals to 0 at short range but have effects in long range. Isn’t it very weird? That’s why I asked if LAMMPS could exclude interactions from long-range energy just now. I’m not sure whether my interpretation is right.

The pair of particles that don't interact at short range,
do in fact interact at long range, with their periodic images.
Hence there is no change in the long-range contribution.
But other codes may partition the work differently. What
I'm confident of is that if you sum the short + long range
Coulombic energies LAMMPS generates you will
get the correct total Coulombic energy, which will
agree with Gromacs or other codes.

I’m using ‘fix rigid/nvt’ to simulate CO2 with fixed bonds and angles. But it outputs an error for my input script and initial configuration: “Fix rigid: Bad principal moments”.

I then decompose the initial configuraion into systems with only one molecule and simulate these one-molecule systems. Most of them do not generate this error any more. But there are still several one-molecule systems continuing to have this problem. I don’t know what causes this error. The structure of the only molecule involved in such systems are correct through visualization or calculation by hand.

I attached my test scripts (for one-CO2 system) with this message. Could you help check with it?

error. The structure of the only molecule involved in such systems are
correct through visualization or calculation by hand.

correct is a matter of opinion in this case.

I attached my test scripts (for one-CO2 system) with this message.
Could you help check with it?

the problem is that your molecules are not _exactly_
linear. the code in fix_rigid that tests for singularities
only recognizes molecules as linear, when one of the
principal moments is _exactly_ zero. i.e. by testing
for identity with 0.0. this is probably not a good idea
to begin with, since even some rounding errors can make
it become slightly different from 0.0.

in your test case the O-C-O angle is not 180 degrees,
that will have it go past the check for linear molecules,
but not the checks for degenerate principal moments.

the problem is that your molecules are not _exactly_
linear. the code in fix_rigid that tests for singularities
only recognizes molecules as linear, when one of the
principal moments is _exactly_ zero. i.e. by testing
for identity with 0.0. this is probably not a good idea
to begin with, since even some rounding errors can make
it become slightly different from 0.0.

I don't think this is the case. In fix_rigid.cppm there are
these lines:

// if any principal moment < scaled EPSILON, set to 0.0

double max;
max = MAX(inertia[ibody][0],inertia[ibody][1]);
max = MAX(max,inertia[ibody][2]);

if (inertia[ibody][0] < EPSILON*max) inertia[ibody][0] = 0.0;
if (inertia[ibody][1] < EPSILON*max) inertia[ibody][1] = 0.0;
if (inertia[ibody][2] < EPSILON*max) inertia[ibody][2] = 0.0;

so I think it is checking against EPSILON which is 1.0e-7.
If this is too strict for your geometry, you might be sloppy
as Axel says. You might also try relaxing the EPSILON
setting (top of file, recompile) and see if that helps.

the problem is that your molecules are not _exactly_
linear. the code in fix_rigid that tests for singularities
only recognizes molecules as linear, when one of the
principal moments is _exactly_ zero. i.e. by testing
for identity with 0.0. this is probably not a good idea
to begin with, since even some rounding errors can make
it become slightly different from 0.0.

I don't think this is the case. In fix_rigid.cppm there are
these lines:

but there are also these lines a little bit further down,
and they cause the problems.

double norm;
for (ibody = 0; ibody < nbody; ibody++) {
if (inertia[ibody][0] == 0.0) {
if (fabs(all[ibody][0]) > TOLERANCE)
error->all("Fix rigid: Bad principal moments");
} else {
if (fabs((all[ibody][0]-inertia[ibody][0])/inertia[ibody][0]) >
TOLERANCE) error->all("Fix rigid: Bad principal moments");
}
if (inertia[ibody][1] == 0.0) {
if (fabs(all[ibody][1]) > TOLERANCE)
error->all("Fix rigid: Bad principal moments");
} else {
if (fabs((all[ibody][1]-inertia[ibody][1])/inertia[ibody][1]) >
TOLERANCE) error->all("Fix rigid: Bad principal moments");
}
if (inertia[ibody][2] == 0.0) {
if (fabs(all[ibody][2]) > TOLERANCE)
error->all("Fix rigid: Bad principal moments");
} else {
if (fabs((all[ibody][2]-inertia[ibody][2])/inertia[ibody][2]) >
TOLERANCE) error->all("Fix rigid: Bad principal moments");
}
norm = (inertia[ibody][0] + inertia[ibody][1] + inertia[ibody][2]) / 3.0;
if (fabs(all[ibody][3]/norm) > TOLERANCE ||
fabs(all[ibody][4]/norm) > TOLERANCE ||

Thanks for your replies. After I modified the coordinates of all CO2 slightly and made them more strictly linear, this error does not show up any more.

However, there is coming a new problem. The system just crashes after several steps and outputs "ERROR: out of range atoms - cannot compute PPPM’. All the thermo outputs look nice before crash and then suddenly those velocity-related properties, e.g kinetic energy, temperature and pressure, become ‘NAN’.

I checked the LAMMPS mailing list archive and found someone met with the same problem around 2 months ago. There has been a long discussion on it. http://lammps.sandia.gov/threads/msg13985.html
And finally the cause of the problem was attributed to some conversion mistakes and the too small Tdamp in fix_rigid/nvt.

My test case is slightly different from Josh’s test system in that discussion. I use ‘real’ units and Josh’s is ‘metal’ and the vdw potentials are slightly different. ‘fix_rigid’ and ‘fix_rigid/nve’ seem to be OK for my system . But when fix_rigid/nvt is used, my system crashes at the 2nd step no matter whether I set Tdamp = 100fs or Tdamp = 100000fs. I then convert my system to the unit ‘metal’. The crash does not occur at Tdamp = 0.1ps, but does still happen at very small Tdamp (eg. 0.001ps) as pointed out in that discussion.

Does it mean there are still bugs in fix_rigid/nvt with ‘real’ units? I attach my test systems (just having one molecule) with this message.

when I decreased the time step from 2fs to 1fs in 'real' units, the crash seems to disappear for Tdamp = 100fs and Tdamp = 100000fs. Bigger time steps, e.g. 1.2fs and 1.5fs, however cause crashing after 2 or 3 steps.

Hi, Trung,
Thanks for your reply. In the ‘metal’ units, the time step I used was 0.002ps which was 2fs and the system did not crash. I thought it was equivalent to use 2fs in ‘real’ units and use 0.002ps in ‘metal’ units. Is my interpretation wrong? Is there any special thing in the code making the two, ‘real’ and ‘metal’, not equivalent simply by unit conversions?