Implementing dynamic charge redistribution and three-body 'pair-wise' potentials in Lammps

Hello,

I'm looking to implement in lammps the potential found on page two of
the attached paper, and I want to make sure I have the best strategy for
this. This first two terms present in the potential are standard
Coulombic and Born-Huggins-Mayer potentials that I can just used out of
the box.

However, for the Coulombic portion, dynamically redistributed charges
are used based on the positions of nearby atoms in every step of a
minimization. How would I implement this sort of calculation inside
lammps? It appears to be similar to a pair-wise calc, but it is not
calculating an energy.

Secondly, the third term is a three body term with a cut-off, so similar
to a pair-wise potential but with an additional atom. Is there any way
you would recommend I implement this sort of potential in lammps that
leverages existing code as much as possible?

Thanks,
Michal Plucinski
Master's Student at Dalhousie University

Step 2 sounds like charge equilibration. There are a couple

potentials in LAMMPS that do it (COMB, ReaxFF). And
Ray (CCd) is working on a more general stand-alone
capability for it. For step 3, look at other 3-body potentails
in the code, like SW or Tersoff.

Steve

Yes, it is a variation of charge equilibration. Since the charge
dependent terms are not quadratic, you can not use the linear
conjugate gradient (CG) method ReaxFF (fix qeq/reax) is using.
Instead you will have to use the damped dynamics (DD) method COMB (fix
qeq/comb) uses.

But you still can not just use fix qeq/comb right out of the box. You
will have to write a pair style to supply the derivatives of energy
with respect to charge and the neutrality constraint to fix qeq/comb
(also you will have to make slight modifications to fix qeq/comb).

A generalized variable charge fix capability will choose either the CG
or the DD solver depending upon the pair style used, but your pair
style still have to supply necessary information.

Ray

Hello,

I have implemented a new three body potential based off the stillinger webber potential, and I complies successfully, but when I attempt to use it with the ELASTIC tutorial, line 82 of in.elastic
variable d1 equal -(v_pxx1-{pxx0})/(v_delta/v_len0)*{cfac}
throws an error: Invalid thermo keyword in variable formula. This error only occurs when I include my potential in the pair_style. What could be the reason for this error?

I have attached the 3 files I developed for this potential:

BO.hk (1.71 KB)

ATT00001.htm (217 Bytes)

pair_hk.cpp (19.6 KB)

ATT00001.htm (217 Bytes)

pair_hk.h (3.82 KB)

ATT00003.htm (4.5 KB)

It’s probably something in the syntax of another variable that
you are referencing in the line you included in your email.

Steve

Hi Steve,

The exact same lammps code, when I include a different pair style (other than the one I developed) works correctly. When I call my code, that’s when the error appears. My confusion arises from the fact that setting a variable doesn’t seem to have anything to do with the potential of the system.

Michal

Hi Steve,

The exact same lammps code, when I include a different pair style (other
than the one I developed) works correctly. When I call my code, that's when
the error appears. My confusion arises from the fact that setting a variable
doesn't seem to have anything to do with the potential of the system.

have you checked, that your code gives meaningful answers with a simple system?
how did you modify the init.mod file?

you would also get this kind of message, if the output for a thermo
keyword is invalid, e.g. NaN or Inf.

axel.

I don’t know. I’m telling you why LAMMPS generates
that error. It has something to do with

the way you defined your variables. That’s where

I would start looking.

Steve

I had not been running optimizations, so the only change in my init.mod file is to set maxiter at zero. Additionally I am using atom_style charge and reading the input using read data from a input file I created from published experimental data.

At what point would thermo call a variable set by a potential? I want to try to track down the part of my code that is generating incorrect output, and I an unfortunately still not very familiar with the structure of lammps code.

When I try to run a minimization on the crystal form of the compound this potential is written for (B2O3), the first step appears to take forever (but does not crash). The initial values set for pxx, pxy, etc are all -nan, unlike with a different potential where they have values. What would this signify?

Thanks,
Michal

I had not been running optimizations, so the only change in my init.mod file
is to set maxiter at zero. Additionally I am using atom_style charge and
reading the input using read data from a input file I created from published
experimental data.

At what point would thermo call a variable set by a potential? I want to
try to track down the part of my code that is generating incorrect output,
and I an unfortunately still not very familiar with the structure of lammps
code.

When I try to run a minimization on the crystal form of the compound this
potential is written for (B2O3), the first step appears to take forever (but
does not crash). The initial values set for pxx, pxy, etc are all -nan,
unlike with a different potential where they have values. What would this
signify?

either that there is a bug in your computation of the virial component
and/or forces, or that there is a units mixup or something else that

a pxx value of -nan will certainly trigger that error message you saw.

axel.

Thank you for your help, what would the virial component be?

Thank you for your help, what would the virial component be?

i don't understand your question. you wrote the pair style, didn't you?
so you should know what you feed into virial accumulator array.

one way to find out where things go bad, is to enable floating point
exceptions. see e.g.

axel.

That’s a very good idea, thank you!

I take it that array is what ev_tally reads into? I hadn’t been able to find a good description of the methods and structure of the code online, so I’ve been attempting to figure out how to develop it by reading the source directly of various pair-styles.

Additionally, I am getting a value of -inf for my potential energy, and I am not totally clear on how that is different from -nan for the force values. Is it infinity vs undefined?

Thanks again!
Michal

That's a very good idea, thank you!

I take it that array is what ev_tally reads into? I hadn't been able to

yes.

find a good description of the methods and structure of the code online, so
I've been attempting to figure out how to develop it by reading the source
directly of various pair-styles.

mainly, you need to read through the header and implementation of the
base class: pair.h and pair.cpp

Additionally, I am getting a value of -inf for my potential energy, and I am
not totally clear on how that is different from -nan for the force values.
Is it infinity vs undefined?

+/-Inf is the result of an overflow or a division of a finite number by zero.
NaN is the result of an invalid math operation, e.g. the log() or
sqrt() of a negative number. this is all defined in the IEEE-754