Changing the LAMMPS source code for applying electric field

LAMMPS LAMMPS Development
1

Hello,

I have a question and I would be thankful if you could help me.

I am trying to change a small part of the Lammps source code to add an Electric Field parameters (ex, ey, ez) to the reaxff_nonbonded.cpp and fix_qeq_reaxff files. These changes solve the problem of unphysical long-range charge transfer which have seen in qeq method. I could add efield to the fix_qeq_reaxff file successfully, but I got some errors for the reaxff_nonbonded.cpp file. What I did was add these command lines to the reaxff_nonbonded.cpp file:

double x, y , z;
double ex_mtn = system->efield->ex;
double ey_mtn = system->efield->ey;
double ez_mtn = system->efield->ez;

x = system->my_atoms[i].x[0];
y = system->my_atoms[i].x[1];
z = system->my_atoms[i].x[2];
en_tmp = KCALpMOL_to_EV * ((system->reax_param.sbp[type_i].chi * q +(system->reax_param.sbp[type_i].eta / 2.) * SQR(q))-q*(xex_mtn+yey_mtn+z*ez_mtn));

In addition, I added " class FixEfield *efield; " in “structure reax_system {}” section of the “reaxff_type.h” file.

I attached the changed files. In these files, I have shown added lines by "//…nina " mark.
I got the error: invalid use of incomplete type ‘class ReaxFF::FixEfield’.

I tried different things but none of them worked. I will be thankful if anyone can help me solve this error.
reaxff_types.h (8.2 KB)
reaxff_nonbonded.cpp (20.5 KB)
fix_qeq_reaxff.cpp (30.3 KB)

2

Unfortunately, it is not so simple to change C++ source code.

You can certainly keep trying and learning in your own spare time, but it looks in your code like you are simply trying to remove the electric field contribution altogether to the electrostatic potential that is put into QEq.

If that is actually what you are trying to do, you make that change in your LAMMPS input script instead of trying to recompile. Simply change your input script from

fix efield all efield v_ex v_ey v_ez // using generic variable names

to

variable e_force_x atom q*v_ex*23.06 // double check force conversion factor
variable e_force_y atom q*v_ey*23.06
variable e_force_z atom q*v_ez*23.06
fix efield all addforce v_e_force_x v_e_force_y v_e_force_z

This calculates the electric field forces (\mathbf{F} = q \mathbf{E}) as x, y and z components in per-atom variables, which are then imposed on the system’s dynamics via fix addforce. Since ReaxFF doesn’t “pick it up” as an electric field, ReaxFF QEq will simply neglect the field imposed in this way.

3

Hello Srtee,

Thank you for your response. I apologize I was not clear enough. Actually, I do not want to remove the electric field contribution.
Chen and Martinez proposed an analytical solution to account for atom polarization, charge conservation, and electronegativity equalization using the charge fluctuation model. They coupled electrostatic energy E(q) [picture 1] with an external electric field to calculate electrostatic properties such as multipole moments and polarizability. The modification made to the electrostatic energy is of the form [picture 2]

To take into account the new distribution of atomic charges and the effect of the coupling with the field, the electronegativity of the atoms has also been modified according to [picture 3].

picture 1:

Screenshot 2023-10-03 163222
Screenshot 2023-10-03 163222894×170 20.2 KB

picture 2:
2

picture 3:
3

What I wanna do is adding the electric field [picture 3] to the electrostatic potential which is defined in nonbonded file in LAMMPS. I have provided some documents in the previous post.

Thank you.

4

What are the \mathbf{R}_\nu and R_{i \nu}?

5

Rv is the spatial positions of atoms (Rx, Ry and Rz).
Riv is the spatial positions of atoms with a counter for number of atoms.

6

Well, this is already implemented in fix qeq/reaxff (see the fix’s documentation, “Restrictions” section – albeit maybe we should cite Chen and Martinez – and source code). If the electric field contribution is somehow not being applied, that is a bug, which you should file with a minimal working example on the LAMMPS GitHub.

Please note that field-driven overpolarization is a known problem with QEq itself. See for example https://pubs.acs.org/doi/abs/10.1021/acs.jctc.1c00975 .

7

Thank you for your explanation. Yes, as you know, in MD interactions between atoms are typically reduced beyond a certain cutoff distance to limit computational cost. But qeq abruptly truncates interactions, which shows unphysical behavior for the system. To mitigate this problem tapering function is induced, which gradually reduces the contribution of interactions (What exactly does ACKS2 method).

As I do not have access to the ACKS2 potential file for my system, I am trying to apply above mentioned changes to the source code of the LAMMPS.

8

Once again. LAMMPS’s fix qeq/reaxff already subtracts the field term \chi \mapsto \chi - \mathbf{E} \cdot \mathbf{x} in its current code. It has done so since Sept 2021 (git blame here) for fix efield with raw numbers and since 2023 (including the version you are trying to modify) with equal- and atom-style variables in fix efield. I have highlighted the relevant lines of source code lines in the expandable section below.

So in the “above mentioned changes”, are you trying to get LAMMPS to do something that it is already doing? That would double-count the electric field contribution. Are you seeing weird results from your current LAMMPS runs – which already include the electric field contribution? Then you can use the fix addforce method above to impose the electric field as an abstract force and not have it affect the QEq charges.

If you are not sure what to do, then you should find a collaborator who will help you (not me, I am a stranger on the Internet). QEq and ACKS2 are both “wrong” models, but QEq is more wrong than ACKS2. As George Box said, all models are wrong, but some models are useful. But as he said in greater detail, which many LAMMPS users should consider very carefully:

Since all models are wrong the scientist must be alert to what is importantly wrong. It is inappropriate to be concerned about safety from mice when there are tigers abroad.

Overpolarization is QEq’s “tiger”, and you will not fix it with a “mousetrap” of this or that electric field correction (indeed, the JCTC paper I cited above demonstrates precisely that). If you don’t understand why, you must find someone (whom you trust) to help you understand.

Source Code Walkthrough

In line 584

if (efield) get_chi_field();

the pre_force function calls get_chi_field if fix qeq/reaxff detects that an efield is being used. The source code of the get_chi_field function starts at line 1091, and line 1131 calculates the actual dot product:

        chi_field[i] = factor*(ex*unwrap[0] + ey*unwrap[1] + ez*unwrap[2]);

which is then subtracted from the initial solution vector guess in the init_matvec() function at line 625:

      if (efield) b_s[i] -= chi_field[i];

The init_storage function performs a similar get_chi_field followed by subtraction.

2 Likes
9

Tgank you dear Srtee for your explanation. You are right. If I want to use “fix acks2/reaxff”, should I use the reaxff potential file which are optimize for this method? Or I can use the reaxff potential file which I have been using for qeq method.

Thank you in advance for your help.

10

No. This has been explained to you before: Two ACKS2-specific parameters - #5 by stamoor