neighbor list problem

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

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.

hmm... a more natural way to me would be to
implement this as a multi-replica method with
multiple topologies right away. no need to mess
with changing between them.

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 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;
     }
  }

ouch! this looks ugly.

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.

well, that is only happening if the factor is != 0.0.
for 0.0, the pair is completely excluded from the neighbor list.

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?

yes. coul/long is applied to all pairs, the tables are not.
the neighbor list code will exclude particles with no
interactions right away. no point to compute what isn't there.

axel.

Hi Axel, thanks very much for your answer. In the method I described before the oxygen O2 is the closer to H at T=0.
However when the reaction is happening the bond change from O1-H to O2-H. Maybe in this case a new oxygen O3 is closer
to H rather than O1. In other word the hydrogen can diffuse during the simulation. Then we'd like to study systems with more than
one hydrogen.
With the multi-replica approach how can I define dynamic topologies in order to allow hydrogen diffusivity?

Regards,

Nicola

Hi Axel, thanks very much for your answer. In the method I described before
the oxygen O2 is the closer to H at T=0.
However when the reaction is happening the bond change from O1-H to O2-H.
Maybe in this case a new oxygen O3 is closer
to H rather than O1. In other word the hydrogen can diffuse during the
simulation. Then we'd like to study systems with more than
one hydrogen.
With the multi-replica approach how can I define dynamic topologies in order
to allow hydrogen diffusivity?

you don't really need a "dynamic" topology,
just swap them occasionally. there are two
parts to that:

1) you have to decide how many topologies you want to superimpose,
    say for example 4 per hydrogen. in my (>10 year old) recollection
    of EVB you don't simply swap topologies, but compute a superposition
    from multiple possible topologies (mimicking the VB approach to QM).
    each of these topologies becomes an individual replica and the actual
    force is a sum over the forces from the individual topologies multiplied
    with the superposition coefficients.

2) depending on some suitable criteria you need to update the set of
   topologies. usually after some rearrangement has happened, either
   after a "reaction" or just regular diffusion. these events would be rather
   rather infrequent and you would handle this like interrupting and
   restarting a simulation and write a little piece of code that. for testing
   purposes, it might be possible to simply set up a set of template input
   files and writing an external code that modifies data files.

conceptually, this would be similar to the parallel replica MD code,
where you run multiple decorrelated MDs and wait for a confirmed
event and then pick the trajectory where the event happened, copy
it to all replica, decorrelate, and continue in parallel. in your case
you would have to replace topologies, which is a bit more complicated.

please note, that there may be other ways to do the same thing, e.g.
by using the approach for the bond order potentials in the manybody
package, where the bond topology is part of the non-bonded calculation
and completely dynamical.

cheers,
     axel.

Dear Axel, all. Thanks for your suggestions, I am going to implement your suggestion and I am doing a preliminary test.
Suppose in my system I have one bond(O1-H) that can change during the simulation. In the mesvb approach the first step is to
identify a second atom, O2 which is the possible candidate for a reaction. So, I've defined a routine find_bond for this purpose.
After this I change the topologies

      swap_topologies(atom1,atom2,j);

and rebuild the neighbor list

    neighbor->build();

The routine swap_topologies will swap the properties of the two atoms, i.e. q, type, bond_type, num_bond, bond_atom, bondlist,molecule,
tag, sametage. Thes is something missing?

In order to test if this approach works I did two test. The first test consist in running two lammps simulation, one with the bond O1-H and the other
with the bond O2-H. The second test consist in running the following: mpirun -np 2 lmp -partition 2x1 -in lmp.in. I ran for 1 timestep with timestep=0.0.
So the energies should be the same.

I printed out bonds, kspace, eng_vdwl, eng_coul energies. Bonds and kspace are computed correctly while eng_vdwl and eng_coul are not after the swapping.
However if I loop over the neighbor list and I change the entries of jlist the energies are correct, e.g.

    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(i==oxygens_id_local[nh][k][0] && t==hydrogens_id[nh])
        {
          jlist[jj] = j_atom2[k];
        }
     if(ii==oxygens_id_local[nh][k][1] && t==hydrogens_id[nh]){
          jlist[jj] = j_atom1[k];
        }
}

This approach looks ugly and dangerous to me but it seems to work. This suggest me that the operation neighbor->build does not rebuild the topology correctly
when the configurations are swapped. Any suggestion on how rebuild correctly the neighbor list in the swapped configuration?

Thanks very much

Nicola