[lammps-users] Point Charge/Massless Atom

[_Joshua_D_Moore]

Dear all,

Is it possible to have a point charge in LAMMPS?

I've read (see below) that it might be possible to add massless atoms,
but making the mass zero leads to "ERROR: Invalid mass value", as the
mass must be > 0.0, which is clearly stated in the manual.

http://sourceforge.net/mailarchive/message.php?msg_id=33ea1b0a0903200736i1f522ec5x4c876001c9dcc442%40mail.gmail.com

I know that there is an option for TIP4P water, but what about a generic site?

Is there a way to add a massless atom?

In particular I want to use a model of N2, which has 2 atomic sites,
and 3 point charges, 2 of which are on the atomic sites and the 3rd is
at the COM of the N2 molecule. I want to use fix rigid, etc. if
possible.

Thanks in advance.

Josh

[_Ahmed_E_Ismail1]

Hi, Josh:

LAMMPS doesn’t have support for generic point charges; if it did, it wouldn’t need a special TIP4P potential!

As for how to add that support, it’s a bit of a challenge. For the N_2 molecule, you’d need to have the capability of adding multiple point charges for a single atom, because two atoms requiring three point charges is a bit of a challenge—unless you can find someway of making the center charge a “mass” site, and redistributing some of the molecular mass accordingly. (Probably wouldn’t satisfy the dynamics if you did, though.)

You’d also need to do a lot of recoding so that the potentials you want to use can adjust the locations of the charges according to their new positions.

–AEI

[sjplimp]

You could try removing the error message check on mass 0
and see if it would work with fix rigid. It clearly would blow up
if you tried integrating w/out fix rigid. But maybe it would just
work. Unlikely, but I haven't thought thru all the issues that
might arise.

Steve

[_Joshua_D_Moore]

Works!

Just commented out the three lines in atom.cpp.

The trick seems to be not to assign a velocity to the massless "atom".
I'm also not using the massless atom in the temperature calculation,
but the kinetic energy seems to not change if I do (because the
velocity is zero), so that's a good sign.

It gives me the total configurational energy that I was expecting so I
think it's working.

fix rigid, fix rigid/nvt seem to be working well too.

Thank you Steve and Ahmed for your help.

Josh

[sjplimp]

I'll have to think whether to allow this in the general version.
I don't really want to check for 0 mass in all the places where
it would be bad (integration, setting velocities, computing temp, etc).
But it is apparently OK to do it within rigid bodies ...

Steve

[_Joshua_D_Moore]

Hi,

This is sort of a follow-up to my post from two years ago, so I'm
putting it in the same thread.
http://lammps.sandia.gov/threads/msg15491.html

Previously we had worked out that you could use a massless site with
fix rigid, by modifying atom.cpp to comment out the lines which check
for the zero mass.

I am working with a 3 site N2 model, with the third site at the center
of mass which contains a partial charge.

The simulation runs and appears to be correct, so long as fix rigid is
used on the entire rigid body and the initial velocities are only set
for the non-massless sites.

My questions today are about the temperature compute.

1) Should the temperature compute be only on the non-massless sites,
or on the entire rigid body? I had presumed because doing it either
way gives the same kinetic energy it is ok either way. However the
temperature changes of course because of the change in degrees of
freedom. Checking the velocities, it appears that fix rigid assigns
velocities to all three sites in the rigid body, even though the
initial temperature was set to zero for the massless site. Should
this be the case? If initial velocites are set for the massless site
initially the simulation crashes, but runs with the massless site
initially at v=0. Since fix rigid is assigning velocites to the
massless site, I believe then that the temperature compute needs to be
on the entire rigid body, but might be wrong.

2) When using thermo_style one, setting the compute only on the
non-massless site gives a different temperature compared to using it
on the entire rigid body (although the kinetic energies are the same
as I mentioned). However, when using thermo_style multi, the
temperature is always the same no matter what I choose for the
temperature compute. Is this a bug? I wouldn't think changing the
thermo_style only would change the result.

I have placed example input and output here. Two examples are given.
Both are for the thermo_style one case (as I mentioned thermo_style
multi gives the same result). The thermo_style multi case agrees with
the case which includes the entire rigid body for thermo_style one.

http://gubbins.ncsu.edu/users/joshua/LAMMPS_N2_Question/nitrogen_question_lammps.tar.gz

One case uses only the N atoms

compute new nitrogen temp
thermo_modify temp new

while the other case uses both the N atoms and the massless site

compute new nitrogen_plus temp
thermo_modify temp new

Thanks in advance for any advice about the compute temperature you can give.

Josh

[akohlmey]

Hi,

This is sort of a follow-up to my post from two years ago, so I'm
putting it in the same thread.
LAMMPS Molecular Dynamics Simulator

Previously we had worked out that you could use a massless site with
fix rigid, by modifying atom.cpp to comment out the lines which check
for the zero mass.

... or setting the mass to a ridiculously low value like 1e-10.

I am working with a 3 site N2 model, with the third site at the center
of mass which contains a partial charge.

The simulation runs and appears to be correct, so long as fix rigid is
used on the entire rigid body and the initial velocities are only set
for the non-massless sites.

My questions today are about the temperature compute.

1) Should the temperature compute be only on the non-massless sites,
or on the entire rigid body? I had presumed because doing it either
way gives the same kinetic energy it is ok either way. However the
temperature changes of course because of the change in degrees of
freedom. Checking the velocities, it appears that fix rigid assigns
velocities to all three sites in the rigid body, even though the
initial temperature was set to zero for the massless site. Should
this be the case? If initial velocites are set for the massless site
initially the simulation crashes, but runs with the massless site
initially at v=0. Since fix rigid is assigning velocites to the
massless site, I believe then that the temperature compute needs to be
on the entire rigid body, but might be wrong.

you get the crash when computing a kinetic energy from a temperature,
since that requires a division by the mass. no problem the other way
around, since that results in a multiplication by zero.

the temperature compute *has* to cover the entire rigid body. fix
rigid will print a warning, if you don't. you have to cover the entire
rigid body, since otherwise fix rigid cannot compute the number of
degrees of freedom it removes from the system. by adding your massless
atoms, you add degrees of freedom, but by including them in a rigid
body they are removed from it again. if they are "hidden" from fix
rigid in the temperature compute, fix rigid will not remove enough
degrees of freedom and you get too low a temperature.

2) When using thermo_style one, setting the compute only on the
non-massless site gives a different temperature compared to using it
on the entire rigid body (although the kinetic energies are the same
as I mentioned). However, when using thermo_style multi, the
temperature is always the same no matter what I choose for the
temperature compute. Is this a bug? I wouldn't think changing the
thermo_style only would change the result.

this seems to be a bug that has been fixed. i get the desired behavior
with the current version.

I have placed example input and output here. Two examples are given.
Both are for the thermo_style one case (as I mentioned thermo_style
multi gives the same result). The thermo_style multi case agrees with
the case which includes the entire rigid body for thermo_style one.

http://gubbins.ncsu.edu/users/joshua/LAMMPS_N2_Question/nitrogen_question_lammps.tar.gz

there is one "anomaly" in this input. you specify the pppm accuracy as
1.0d-8 and not 1.0e-8. is that intentional? if not, you should be
aware that 1.0d-8 will be parsed as 1.0 and not 10**-8 and hence your
kspace calculation is much less accurate than what you might expect.
1.0d-8 is only a valid representation for 10**-8 in fortran. in C/C++
all exponential numbers are assumed to be double precision by default.

axel.

[sjplimp]

Once your massless atom is in the rigid body, there
is no reason I can think of to worry about it being massless.
You can use fix rigid/nvt if you want to thermostat the bodies.
You can also monitor the temperature of the rigid bodies
from output of the fix. The fix will not care that one particle
in the rigid body is massless. Nothing you do with other
commands in LAMMPS should include the rigid body atoms,
e.g. time integration, thermostatting on the solvent.

A accuracy of 1.0e-8 for PPPM is way beyond what is needed
for typical problems. 1.0e-4 is sufficient.

Steve

[akohlmey]

A accuracy of 1.0e-8 for PPPM is way beyond what is needed
for typical problems. 1.0e-4 is sufficient.

we had plenty of cases here on the list that have shown that 1.0e-4 is
on the sloppy side and that 1.0e-5 would be better.

i'd also like to point out that i've found the accuracy estimates are
good for systems where you have a typical and evenly distributed
density of charged particles. for the coarse grain systems that i
occasionally deal with, where only few particles are charged and they
are not evenly distributed in space, the estimator is not working so
well and using suitable parameters manually leads to an estimated
accuracy in the neighborhood of 1.e-8 quite frequently.

[_Joshua_D_Moore]

Dear Axel and Steve,

Thank you for your quick comments. Sorry for my delay in responding.

Axel, you are right, the 1.0d-8 was me thinking in Fortran

We are doing some adsorption with the other part in GCMC, so probably
will work out what exact precision we need. The 10^-8 was kind of
just an extreme case for testing, but that is a good point of checking
this. No need to use such an extreme case if not needed.

Just to confirm, it is proper to include all particles in any compute
for the rigid N2 with 3 sites, even if they have zero mass? This
makes sense.

Thanks for checking with the latest version of LAMMPS. I should have
checked this before sending my message. I was using the Jan 20
version.

[akohlmey]

Dear Axel and Steve,

Thank you for your quick comments. Sorry for my delay in responding.

Axel, you are right, the 1.0d-8 was me thinking in Fortran

We are doing some adsorption with the other part in GCMC, so probably
will work out what exact precision we need. The 10^-8 was kind of
just an extreme case for testing, but that is a good point of checking
this. No need to use such an extreme case if not needed.

BTW: if you are looking for extreme accuracy, you also have to turn
off tabulation for coulomb.

Just to confirm, it is proper to include all particles in any compute
for the rigid N2 with 3 sites, even if they have zero mass?

yes. have a look at the comments in FixRigid::dof().

axel.

[sjplimp]

Just to confirm, it is proper to include all particles in any compute
for the rigid N2 with 3 sites, even if they have zero mass? This
makes sense.

What "compute" do you mean? If you mean a termperature compute,
then you will only get translational temperature of the rigid bodies.
If you want the full temperature, including rotational, then you should
access it through fix rigid, which does that computation.

What other compute uses the mass of atoms inside rigid bodies in
a way that tells you anything meaningful?

Steve

[_Joshua_D_Moore]

Dear Steve,

I am confusing things. Sorry.

I think I understand what you are getting at. fix rigid, rigid/nve
and rigid/small produce a global scalar temperature that includes the
translational rotational degrees of freedom. fix rigid/nvt does not.
It produces a global scalar but is the cumulative energy change due to
the thermostatting from the manual. Is this what you were referring
to?

So you can't compute the full temperature with fix rigid/nvt?

Or is (where nitrogen_plus includes all particles in the rigid body).

compute new nitrogen_plus temp
thermo_modify temp new

doing the trick?

I think you may be right and I am missing the rotational degrees of
freedom in the temperature because my first step is like 297 K and not
the 300 K which I set with the velocity command?

I want to be able to compute the temperature of only the rigid bodies
separately, because I will be putting these inside of an adsorbate
which is frozen (zero velocity) and is not included in any time
integration.

Thanks for your patience.

Josh

[akohlmey]

Dear Steve,

I am confusing things. Sorry.

I think I understand what you are getting at. fix rigid, rigid/nve
and rigid/small produce a global scalar temperature that includes the
translational rotational degrees of freedom. fix rigid/nvt does not.

how do you come to that conclusion?

It produces a global scalar but is the cumulative energy change due to
the thermostatting from the manual. Is this what you were referring
to?

you are not making any sense here.

So you can't compute the full temperature with fix rigid/nvt?

sure you can.

Or is (where nitrogen_plus includes all particles in the rigid body).

compute new nitrogen_plus temp
thermo_modify temp new

doing the trick?

that is just what you'll use as the thermo output.
fix rigid/nvt will create a temperature compute for its own purposes
anyway and it will cover all atoms in the rigid bodies.

I think you may be right and I am missing the rotational degrees of
freedom in the temperature because my first step is like 297 K and not
the 300 K which I set with the velocity command?

i wouldn't worry so much about that when using a thermostat.

I want to be able to compute the temperature of only the rigid bodies
separately, because I will be putting these inside of an adsorbate
which is frozen (zero velocity) and is not included in any time
integration.

so you need to change the thermo compute to cover all atoms that are
time integrated. the rest will work automatically.

axel.

[sjplimp]

Comments below.

Steve

Dear Steve,

I am confusing things. Sorry.

I think I understand what you are getting at. fix rigid, rigid/nve
and rigid/small produce a global scalar temperature that includes the
translational rotational degrees of freedom. fix rigid/nvt does not.
It produces a global scalar but is the cumulative energy change due to
the thermostatting from the manual. Is this what you were referring
to?

I was really referring to the different ways that fix nvt and fix nvt/rigid
will measure the temperature of the rigid bodies and do the thermostatting.

Fix nvt with a temperature compute like compute temp, will compute
the temp (via the compute) by summing by 3/2 kT = sum (1/2 m v^2) for
all the atoms. Then it will remove degrees of freedom (DOF) to
leave only 6 DOF per rigid body. That is a reasonable estimate
of the temp of the rigid bodies. However, the thermostatting will simply
rescale the v of all the atoms, which is not a good way to thermostat,
b/c the velocities of atoms in the bodies are reset by fix rigid.

You can use fix rigid with fix langevin (and no fix nve) and do
a reasonable job of thermostatting at least the translational
DOF of the bodies.

So you can't compute the full temperature with fix rigid/nvt?

It's not an output of the fix, as you note. But compute temp
will remove DOF and give a good estimate of the temp.

Or is (where nitrogen_plus includes all particles in the rigid body).

compute new nitrogen_plus temp
thermo_modify temp new

doing the trick?

Assuming the group n_plus has your rigid bodies,
then yes that is the estimated temp. The code for
computing the exact temp (trans + rot DOF) is
in fix rigid, so we need to think about how to
make that value accessible from fix rigid/nvt. Trung
chose to make the energy delta accessible which
is also a useful value. The issue is that there
is also a global array made visible, which means
we can't also have a global vector of 2 values,
like temp and energy.

[_Joshua_D_Moore2]

Thanks Steve and Axel. I think I understand now