Dear LAMMPS developers,
We are developing a model that requires charge equilibratation, similar to QEQ, but where electronegativity(chi) and atomic hardness (eta) parameters are predicted for each atom individually and not depends only on species type as it is in fix qeq/point, for example.
Is it possible to somehow re-use classes from QEQ package for such charge equilibration? What would be your suggestion for code design ?
Best regards,
Yury
It should be possible with small changes.
From what I can see quickly using grep, it appears that you may be able to declare the function compute_eneg() in qeq/dynamic and qeq/fire as virtual and then override it in a custom derived class.
Since in those two cases there is only one location in exactly that function that accesses the per-type eta and chi properties. Similar for qeq/point, qeq/slater, and qeq/shielded the init_matvec() function.
Finally, for all variants, you would also need to override the sparse_matvec() function.
In all cases, the modification to accommodate a per-atom parameter instead of a per-type parameter should be obvious and straightforward.
You can also come join the ELECTRODE family! We have a conjugate gradient method and our electric potential evaluation includes the Kspace contributions. You’d just have to add routines to read in chi and eta as per-atom variables and then feed them into the CG iterations.
@akohlmey and @srtee Thank you for your reply.
Currently I’m leaning towards extending/overriding fix qeq/point class. I have already some working option, but I want to ask further advice about code design.
Since in our case chi and J (=eta) are also environment dependent, they will contribute to the force via q[i]*d_chi/d_r +0.5*q[i]**2 * d_J/d_r derivatives. That means we need to have equilibrated charges at the moment of this contribution to be computed. So, the current working scheme / method call sequence is following:
# in.lammps
...
pair_style hybrid/overlay pace/qeq coul/cut 7
pair_coeff * * pace/qeq potential.yaml Cl Na
pair_coeff * * coul/cut
fix 2 all qeq/point/pace 1 7.0 1e-6 100 pace/qeq
...
pace/qeq::compute (in pseudocode):
// stage 1: compute only chi and J
for i = 0 .. list->inum-1
compute chi and J for atom i, store in global per-atom array chi[i] and eta[i]
end for
// stage 2: charge equilibration
explicit call fix::pre_force to do charge equilibration, using chi[i] and eta[i]
// stage 3: full compute
for i = 0 .. list->inum-1
compute short range energy e_atom, pair neighbour forces `fij`
compute compute pair neighbours d_chi/d_r and d_J/d_r
// add local QEQ contribution to forces
fij -= q[i]*d_chi/d_r +0.5*q[i]**2 * d_J/d_r
accumulate `fij` into global `f` array
end for
pair coul/cut is responsible for Coulomb interation between different atoms.
fix qeq/point/pace has no direct pre- or post-force calls from LAMMPS, but being called from pair::compute.
This is the current draft implementation scheme. Would this match the LAMMPS philosophy and could it be merged into the main branch after a pull request?
And another quick question: is it possible to have non-zero total charge for non-periodic system in charge equilibration ?
You should still call the fix’s pre_force directly from the timestep loop, rather than from pair.
Firstly this lets the code parallelize better. The QEQ procedure requires at least some all-to-all communication and so every process needs to reach and leave it at the same time. Doing so in the pre_force part of the timestep is preferable, as compared to sticking it in pair where different processes might take different amounts of time to do the prior pair operations (or the subsequent bond, etc. operations).
Secondly, it makes for easier testing and coding. Calling a fix from a pair means you have to make them friend classes, and put in switches to turn the fix call on and off if you want to test them separately. LAMMPS already builds in encapsulation between fixes and pairs that make independent development and testing easier, so don’t throw that away.
Lastly, if you call your QEQ from inside pair pace/qeq, then anyone who does pair hybrid coul/cut ... pace/qeq will get wrong results as compared to pair hybrid pace/qeq ... coul/cut. You could certainly write check code to raise an error with the former command, but then you need to interact with pair hybrid and so on, whereas calling the fix’s pre_force within the usual timestep loop just needs your pace/qeq to check for fix qeq/..., which is much more straightforward.
I totally agree with you, but for the moment I do not understand how to deal with the requirement that I need to do fix::qeq right in the middle of pair::compute (i.e. stage 2), as I need chi and J computed…
Why not calculate (or retrieve) chi and J in the fix’s pre_force as well?
EDIT: and after that you can also add the “charge reaction forces” in via post_force. That way you can do it all in a fix, instead of having to code up a pair style, since I am guessing you will not use any of the “machinery” around pair styles (pair coeffs, mixing, tail and shift options, etc).
This is good point. So in pre_force I will have access to neighbour list as in pair::compute, right?
I would need for that pair derivatives, i.e. for each central atom i I have pair derivatives wrt. it’s each neighbors positions. So, in pair:: compute I still need to accumulate them and communication between pair style and fix is unavoidable
You have to create a neighbor list request within the fix. However, If the settings of the request are identical to the ones from the pair style, the neighbor list code will return a pointer to the same neighbor list. You can monitor this in the neighbor list summary output at the beginning of the first run (copy from …).
Fix qeq/reaxff supports electric fields, which uses a per-atom chi, see lammps/src/REAXFF/fix_qeq_reaxff.cpp at c5d9f901d9e8336602dee9558ffd6a5aab30e4d1 · lammps/lammps · GitHub. It wouldn’t be hard to extend that to values read from a file. Additionally qeq/reaxff is ported to Kokkos.
Dear LAMMPS developers,
we have another question regarding the Qeq model. How do you deal with the PBC case? According to our knowledge, different matrix elements should be considered in this case (i.e. reciprocal and real space terms, self-interaction, neutralization term, etc.) if Ewald summation is used. However, we did not find any mention of those terms in the function compute_H of fix_qeq_point.
Best regards,
Matteo
It helps other forum readers and posters if one thread is used to discuss one topic, so please start a new thread with your question. Thanks!