Please reply to the LAMMPS mailing list, not to me.
Hello,
I apologize for sending my reply to wrong the email.
I still have a few questions regarding using Qeq in LAMMPS.
The first charge model I wanted to try out was qeq/point. Basically I wanted to implement a meam potential plus an electrostatic potential and also allow for charge equlibration. I was under the impression that the qeq fix did nothing in term of calculating energies and only varied charges on atoms. I therefore decided I needed to use the overlay potential to first define my MEAM + electrostatics potential. I chose to use coul/long as my electrostatic part. I defined the meam + electrostatic potential in the following way while also and defining a k_space style for the long range electrostatics
pair_style hybrid/overlay meam coul/long 4.8
pair_coeff * * meam /home/vella/Potentials/libraryLi_H_D_v2.meam Li2NN D/home/vella/Potentials/Li_H_v2.meam Li2NN D
pair_coeff * * coul/long
kspace_style pppm 1.0e-4
I then implemented qeq/point to allow for charge equilibration during the simulation.
fix 2 all qeq/point 1 4.8 1.0e-6 200 param1.qeq
My question now is this the correct way of implementing a MEAM + Electrostatics (with variable charges using a point charge model) in LAMMPS?
Next question: Equation 1, in the Streitz-Mintmire paper shows the form of the electrostatic part of the potential as the self part (E_i(q_i)) + the coulomb pair interaction Vij. This form exists before they choose a charge model. Therefore, shouldn’t the self part be part of the potential even if you choose a point charge model or am I missed something here?
Thank you,
Joseph Vella
Hello,
I apologize for sending my reply to wrong the email.
No worries. 
I still have a few questions regarding using Qeq in LAMMPS.
The first charge model I wanted to try out was qeq/point. Basically I
wanted to implement a meam potential plus an electrostatic potential and
also allow for charge equlibration. I was under the impression that the
qeq fix did nothing in term of calculating energies and only varied charges
on atoms. I therefore decided I needed to use the overlay potential to
first define my MEAM + electrostatics potential. I chose to use coul/long
as my electrostatic part. I defined the meam + electrostatic potential in
the following way while also and defining a k_space style for the long
range electrostaticspair_style hybrid/overlay meam coul/long 4.8
pair_coeff * * meam /home/vella/Potentials/libraryLi_H_D_v2.meam
Li2NN D/home/vella/Potentials/Li_H_v2.meam Li2NN D
pair_coeff * * coul/longkspace_style pppm 1.0e-4
I then implemented qeq/point to allow for charge equilibration during
the simulation.fix 2 all qeq/point 1 4.8 1.0e-6 200 param1.qeq
My question now is this the correct way of implementing a MEAM +
Electrostatics (with variable charges using a point charge model) in LAMMPS?
Looks okay, but your cutoff may be too short for both coul/long and
qeq/point.
Next question: Equation 1, in the Streitz-Mintmire paper shows the form
of the electrostatic part of the potential as the self part (E_i(q_i)) +
the coulomb pair interaction Vij. This form exists before they choose a
charge model. Therefore, shouldn't the self part be part of the potential
even if you choose a point charge model or am I missed something here?
Sorry I confused myself and you in my previous email. I was thinking about
another self term. Anyhow, if you follow the paper Eq 1 in the paper is
replaced by Eq 7, and this equation is the master equation that is
implemented. So you don't need E_i(q_i), or to put it more correctly, this
term is transformed to Sum_i(q_i*Chi_i).
Ray
Thank you for all your help Ray.
I’ll definitely extend the cutoff distance for coul/long and qeq/point.
A few more questions. (Sorry!)
Maybe I’m a little confused about the Sum_i(q_i*Chi_i) term in the energy.
If I understand correctly the energy should be E = E_es + E_eam where E_es is electrostatics and E_eam is from eam (or meam) in my case.
The E_es term can be written as (Eq 1) E_es =Sum_i(E(q_i)) + 1/2sum(V_ij).
If I assume a point charge model the second term in this equation is simply the coulombic interaction (which is what I use coul/long for), however don’t I still need to somehow account for the Sum_i(E(q_i)) term or am I just overthinking this? It is unclear to me what happens to this term when if I just want to use the point charge model.
Finally if I eventually decide to just use the full Strietz-Mintmire potential, I assume the pair_style coul/streitz takes care of the Sum_i(E(q_i)) term (which as you stated is eventually transformed to the Sum_i(q_i*Chi_i) term). Is this correct? (I would of course have to use qeq/slater for charge equilibration as the documentation indicates.)
Again, thank you for your help and patience!
Best,
Joseph Vella
Thank you for all your help Ray.
I'll definitely extend the cutoff distance for coul/long and qeq/point.
A few more questions. (Sorry!)
Maybe I'm a little confused about the Sum_i(q_i*Chi_i) term in the
energy.
If I understand correctly the energy should be E = E_es + E_eam where E_es
is electrostatics and E_eam is from eam (or meam) in my case.
The E_es term can be written as (Eq 1) E_es =Sum_i(E(q_i)) +
1/2sum(V_ij).
E_0 in Eq. 7 accounts for charge-independent terms (albeit eam, meam,
tersoff, etc).
If I assume a point charge model the second term in this equation is
simply the coulombic interaction (which is what I use coul/long for),
however don't I still need to somehow account for the Sum_i(E(q_i)) term
or am I just overthinking this? It is unclear to me what happens to this
term when if I just want to use the point charge model.
Eq 1 is replaced by Eq 7 and Sum_i(E(q_i)) is accounted for by
Sum_i(q_i*Chi_i).
It does not matter which charge model you use, Sum_i(q_i*Chi_i) is always
there.
Finally if I eventually decide to just use the full Strietz-Mintmire
potential, I assume the pair_style coul/streitz takes care of the Sum_i(E(q_i))
term (which as you stated is eventually transformed to the Sum_i(q_i*Chi_i)
term). Is this correct?
No. Sum_i(q_i*Chi_i) exists for all charge models. See 1991 Rappe &
Goddard and 1994 Rick, Stuart, and Berne.
Ray
Thank you Ray but I think there is still some sort of miscommunication, (due to me not formulating my question correctly).
I’ll focus on just using a point charge model for the time being.
I realize the sum_i(q_i*chi_i) exists for all charge models, my question is related to whether or not LAMMPS calculates this term and other charge dependent one-body terms during the energy calculation when using a point charge model.
Let’s say, for example I want to augment an existing MEAM potential with an electrostatic term.
I would therefore use the following command.
pair_style hybrid/overlay meam coul/long 12.0
pair_coeff * * meam /home/vella/Potentials/libraryLi_H_D_v2.meam Li2NN D/home/vella/Potentials/Li_H_v2.meam Li2NN D
pair_coeff * * coul/long
This amounts to the potential energy being:
E = E_meam + sum (Cqiqj/rij)
Now I would also like to allow for charge equilibration. Since I am using coul/long for my electrostatics term (which assumes point charges), I figured qeq/point would be the appropriate fix to use.
However, when looking at the form of the potential energy from the hybrid/overlay of the meam and coul/long, the only charge dependent term is the two-body term sum (Cqiqj/rij). It seems I am therefore missing the charge dependent term one-body terms (such as sum_i(q_i*chi_i)) when calculating the energy of the system. As indicated in Streitz-Mintmire paper these charge dependent one-body terms appear in the energy calculation regardless of the charge model assumed (although their forms will depend on the charge model).
Is there a simple way to add these charge dependent one-body terms into the calculation of the energy?
Thank you,
Joseph Vella
Okay, I see your argument now. You are raising a good question here: should the energy terms used to derive equilibrium charges and that used to calculate energy and forces be the same?
The answer is yes for most of the variable charge potentials. For example, COMB and ReaxFF both include the “self energy” term, i.e., chi_iq_i + 0.5 eta_iq_i^2, to describe the ionization energy and electron affinity - aka energy require to form a charge. However, when qeq/variant is used with a simple coul/variant that does not have the self energy term, should the qeq/variant additionally include it? I think yes, otherwise there would be a discrepancy in the terms used for energy/force calculations and charge equilibration.
A clean and easy modification would be to add a “fix_modify energy yes” to fix qeq/variants to add this self energy term to the pair_style that does not have this term already, e.g., coul/cut, coul/long, etc. Coul/streitz already includes this so it wouldn’t be necessary.
I will work on this if I find some time next week.
Thanks,
Ray
Great! Thank you so much for your help.
I suppose for now I will focus on using qeq/slater and coul/streitz since the self terms are included.
I suppose I can check the announcement page to see when the “fix_modify energy” yes is implemented for the other cases?
Thank you again,
Joseph Vella
A word or two to share…
In the gas adsorption community (more specifically when modeling CO2 adsorption by MOFs) people have been using a hybrid combination where sometimes QEq type methods are parametrized to simulate the MOF framework charges but the energy calculation within the GCMC part of the adsorption process is done using purely Coulombic terms). Part of the rationale there sometimes is to keep the computational cost of the GCMC part low (or so it is claimed). The other one could be the lack of available software with charge equilibration capabilities and/or the lack of expertise to deal with the potential numerical instabilities from atomic overlaps. I know of other cases as well where the same hybrid approaches have bee tried in different fields by some of the QEq pioneers (references excluded for professional reasons). Thus I believe full consistency is a bit blurry when it comes to this story. After all one should remember that partial atomic charges are not quantum observables therefore their mere definition is a bit of an ad-hoc construction.
Carlos
Hi Carlos,
You raise a good point. It seems to me that running GCMC without
including the QEq energy is formally incorrect and will generate
adsorption isotherms that differ from those obtained from a pure MD
method, such as thermodynamic integration or an explicit gas
reservoir. However, in practice the difference may well be negligible.
Also, I can imagine lots of practical difficulties with incorporating
QEq energy change for particle insertion. That said, it should be
possible for someone to try this out using the LAMMPS GCMC fix, using
the full_energy keyword, and a QEq variant that adds its energy to the
total potential energy.
Aidan