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.
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.
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:
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.
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.
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?
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
is bad happening.
a pxx value of -nan will certainly trigger that error message you saw.
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?
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
standard and its updates.