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 (C*qi*qj/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 (C*qi*qj/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_i*q_i + 0.5 eta_i*q_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