Bug in units

Hello, LAMMPS users and developers!

I want to report a small but pretty nasty bug: the man page for “nano” set of units
says:

  • charge = multiple of electron charge (+1.0 is a proton)

while in fact in code it is

force->qelectron = 1.6021765e-19;

making it Coulombs!

Also, as I have just noticed, there’s a similar matter with style micro: coulombs instead of listed picocoulombs.

Best of luck,
Pavel Zun

force->qelectron = 1.6021765e-19;

The only thing this affects is the Kspace error estimation

(also the pair zbl potential), not even the KSpace calculation.

It also does not affect Coulombic pairwise calculations. That is

in the conversion factor qqr2e, which appears to be correct.

So the doc page units for charge in nano/micro are correct.

Steve

fixed the qelecton settings for micro/nano

thanks,
Steve

Hello Steve,

I’m not sure if qelectron is the offender, but the attached sctipt behaves… weird, like the particles are actually charged for about 1 e (in the file the charge is set to 0.00000000000000000008).

Also, slightly larger charge results in energies on the order of 10e+40 and crashes everything.

Same simulation behaved correctly in lj units.

Pavel

in_nano.dpd (2.01 KB)

I don’t know - maybe you are getting overlapping particles
and it blows up. If you are saying you have run the
same simulation in nano and LJ units, do you get identical
answers for some number of timesteps, then a blow-up?

I doubt that, b/c it should be possible to run identical

simulations (except for round-off issues) in any unit
system.

Steve

There shouldn’t be any overlapping particles, since I run minimization first. The weird thing is, it seems like the Coulomb forces are many orders of magnitude stronger than they’re supposed to be — even a very small charge results in ridiculous energies.

Also, I’m sorry that I can’t provide much feedback about this issue — I have since switched to real units, and everything works smoothly there.

Pavel

I repeat my previous Q - when you run the same problem
in 2 different units, are you getting identical answers, from step 0
onwards, even for the minimization? To make this comparison, you have to convert your
output (e.g. thermo output) from one set of units to the other,

e.g. ergs to Kcal/mole.

You can do the same for dump output, e.g. dumping the force
on each atom. You are saying it crashes eventually. I am
asking if it is identical from the very first timestep or min iteration.

If it isn’t then it’s likely your 2 input scripts are in fact not
identical, and trying to figure out why it crashes later on

is not the place to start looking.

Steve

Yeah, I’m sorry. By “identical” I meant “same script with different constants”. I didn’t keep all propotions from LJ when transitioning to nano.
On the other hand, I kept most stuff the same when I switched from nano to real — and it behaves very differently, i.e. it doesn’t behave “charged” when particles are charged for 10^-19 e. In order for that to happen, they need to be charged for about 1 e.
Also, nano doesn’t crash eventually — with more significant charges it crashes right away.

Sorry, at this point, I’m not sure what you are asking anymore.

a) I don’t see any issues with the doc page or internal implementation
of “nano” units.

b) I believe it should be possible with any two unit choices

in LAMMPS, e.g. nano vs lj, to have 2 input scripts, one

for each units, that produce “identical” simulations, for some
number of timesteps, until they diverge due to round-off differences.

If you are not trying to do (b), but simply have a new model

in nano units, and are asking why it doesn’t behave as you
expect, then that’s a different kind of Q.

Steve

Okay.

When I run the attached script, it crashes. Note that the particles are charged for 0.000000000000000008 electron charge each. The particles weigh about 1000 atomic mass units each.

Step Temp E_pair E_mol TotEng Press
78 300 1.1934027e+10 0 1.1934058e+10 76500.503
ERROR on proc 4: Particle on or inside fix wall surface (…/fix_wall_lj126.cpp:81)

This happens right after minimization.

When I set charge to zero, it runs normally.

Step Temp E_pair E_mol TotEng Press
123 300 549.10713 0 31607.526 1.1250615
200 328.98064 205.67379 0 34264.402 0.10671868

Note that pairwise energies and pressure drop drastically.

The most logical explanation of this behaviour is that the particles are not, in fact, charged for 0.000000000000000008 electron charge each. They seem to be charged for a significantly higher amount. Perhaps, for the same amount of coulombs.

I mentioned LJ units to show that the ridiculous energies for relatively small charges are not an intrinsic element of this script.

P. S.: It might become more clear if you run this script yourself as is, and then set the charge for type 1 particles to 0 and run it again. Alternatively, you can set it to 0.00000000000000000000000000000000008. That works too.

P. P. S.: This issue persists in the version dated 31Mar14.

Pavel

in_nano_test (1 KB)

I’ve never run a simulation in nano units. There could be
any number of problems with the model as you defined it

in the script you sent.

Your original email said:

I’m not sure if qelectron is the offender, but the attached sctipt behaves… weird, like the particles are actually charged for >about 1 e (in the file the charge is set to 0.00000000000000000008).
Also, slightly larger charge results in energies on the order of 10e+40 and crashes everything.
Same simulation behaved correctly in lj units.

Can you post two scripts, one for nano, one for LJ units.

Preferably simpler than the one you sent. E.g. no

minimzation, smaller number of particles. Just initialize
the system and run a couple timesteps of dynamics.
Show me that the temp, energy, pressure
is identical (or not) at each step for the 2 cases, once you

do the unit conversion for the output.

In fact, for the issue you are focusing on with
charge, you should be able to run a system
with 2 particles, separated by some distance.

I am guessing your two scripts will not be
for identical systems. Attached are 2 scripts
that do this test for an Au lattice. The
output is not exactly identical b/c they
were written for an old version of LAMMPS and
the conversion factors in the current code
are now more precise. But it’s very close. They illustrate

how to do the comparison I’m talking about.

You’ll also note that none of the numeric
values in the real-units scripts are nice round
LJ-like values. Many of the numbers in

your nano-units script are, as if they were copied
from a LJ-units script.

That makes me guess you are not in fact
running the same system with your 2 scripts.

Steve

in.au.lj (446 Bytes)

in.au.real (585 Bytes)

Cancel what I wrote - I was wrong. There is an invalid conversion
factor for charge in the nano unit setup in src/update.cpp. Related
to the confusion on whether a proton is input as 1.6e19 Coulomb or 1.

The latter is now what the doc page says.

This mixup probably happened when the person who contributed the

nano units code chunk assumed the units for charge would be Coulombs

and it got switched around at some point.

Try this line in update.cpp in the nano section:
force->qqr2e = 230.7078669;

and see if your problems go away.

Steve