Dear users, I'm trying to implement the empirical valence bond method in LAMMPS.

This should be useful to study phenomena like proton diffusivity. Let say we have 3 atoms, O1, O2, H.

Initially there is a bond O1-H, this will be our 1 configuration. The second configuration will be with the bond O2-H.

With this approach I'd like to compute, at each step, forces and energy for the two configuration and to choose which is the most energetically favorable.

So, the idea is to define a new fix and to use the methods pre_force and post_force to manage those computations.

Particularly, the method post_force will compute forces and energy in the second configuration and will decide if the reaction is happening.

In my tests I used a BaZrO3 system with 1 hydrogen atom. My pair interactions are defined as

pair_style hybrid/overlay coul/long 10.0 table linear 2001

pair_coeff * * coul/long

pair_coeff 1 3 table bazro3-ff.pot tablepot_O2_Ba 10.0

pair_coeff 1 4 table bazro3-ff.pot tablepot_O1_Ba 10.0

pair_coeff 2 3 table bazro3-ff.pot tablepot_O2_Zr 10.0

pair_coeff 2 4 table bazro3-ff.pot tablepot_O1_Zr 10.0

pair_coeff 3 3 table bazro3-ff.pot tablepot_O2_O2 10.0

pair_coeff 3 4 table bazro3-ff.pot tablepot_O2_O1 10.0

pair_coeff 3 5 table bazro3-ff.pot tablepot_H1_O2 10.0

pair_coeff 4 4 table bazro3-ff.pot tablepot_O1_O1 10.0

pair_coeff 4 5 table bazro3-ff.pot tablepot_H1_O1 10.0

In order to implement the computation of forces and energy in the second configuration I had to swap bond type, atom type, charges.

However, in order to compute the interactions in the new configuration it is also necessary to modify the neighbor list. The following

code is the implementation

for(int ii=0; ii<inum; ii++){

i=ilist[ii];

jlist = firstneigh[i];

jnum = numneigh[i];

for (int jj = 0; jj < jnum; jj++) {

j = jlist[jj];

t = j&NEIGHMASK;

if((t==oxygens_id[0][0] && ii==hydrogens_id[0]) || (ii==oxygens_id[0][0] && t==hydrogens_id[0]))jlist[jj]=t;

if((t==oxygens_id[0][1] && ii==hydrogens_id[0]) || (ii==oxygens_id[0][1] && t==hydrogens_id[0]))jlist[jj]=j_atom1;

}

}

I need to do this kind of trick because in the pair_style there is this line

factor_lj = special_lj[sbmask(j)];

which tell to exclude the O1-H from the pair interaction computation. In the second configuration I want to include O1-H and exclude O2-H which is the new bond.

This change produce the desired results in the coul/long case but not in the tablepot case.

So, those two potential are using the same neighbor list? From what I'm seeing the answer is no.

I did some debugging and the variable jnum seems to have different values during the execution. Is this what is expected?

Thanks a lot for you help

Nicola