Fix efield/tip4p

LAMMPS LAMMPS Development
, ,
1

Dear all,

I am trying to simulate the water behavior under electric field, and I want to use TIP4P/ice model. According to the suggestions mentioned at TIP4P , kspace_modify slab , efield !
. I modified the code of fix efield in lammps, and re-compiled LAMMPS to conduct simulation.

The main modification is in FixEfieldTIP4P::post_force, line 353-429.

      if (itype != typeO) {
        f[i][0] += fx;
        f[i][1] += fy;
        f[i][2] += fz;

        domain->unmap(x[i],image[i],unwrap);
        fsum[0] -= fx*unwrap[0]+fy*unwrap[1]+fz*unwrap[2];
        fsum[1] += fx;
        fsum[2] += fy;
        fsum[3] += fz;

        if (evflag) {
          v[0] = fx*unwrap[0];
          v[1] = fy*unwrap[1];
          v[2] = fz*unwrap[2];
          v[3] = fx*unwrap[1];
          v[4] = fx*unwrap[2];
          v[5] = fy*unwrap[2];
          v_tally(i, v);
        }

      } else {
        iH1 = atom->map(tag[i] + 1);
        iH2 = atom->map(tag[i] + 2);
        if (iH1 == -1 || iH2 == -1)
         error->one(FLERR,"TIP4P hydrogen is missing");
        if (atom->type[iH1] != typeH || atom->type[iH2] != typeH)
         error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
        //iH1 = domain->closest_image(i,iH1);
        //iH2 = domain->closest_image(i,iH2);

        fO[0] = fx*(1 - alpha);
        fO[1] = fy*(1 - alpha);
        fO[2] = fz*(1 - alpha);

        fH[0] = 0.5 * alpha * fx;
        fH[1] = 0.5 * alpha * fy;
        fH[2] = 0.5 * alpha * fz;

        f[i][0] += fO[0];
        f[i][1] += fO[1];
        f[i][2] += fO[2];


        f[iH1][0] += fH[0];
        f[iH1][1] += fH[1];
        f[iH1][2] += fH[2];

        f[iH2][0] += fH[0];
        f[iH2][1] += fH[1];
        f[iH2][2] += fH[2];

        domain->unmap(x[iH1],image[iH1],unwrapH1);
        domain->unmap(x[iH2],image[iH2],unwrapH2);


        fsum[0] -= fO[0]*unwrap[0]+fO[1]*unwrap[1]+fO[2]*unwrap[2]+fH[0]*unwrapH1[0]+fH[1]*unwrapH1[1]+fH[2]*unwrapH1[2]+fH[0]*unwrapH2[0]+fH[1]*unwrapH2[1]+fH[2]*unwrapH2[2];
        fsum[1] += fx;
        fsum[2] += fy;
        fsum[3] += fz;

        if (evflag) {

          v[0] = fO[0]*unwrap[0]+fH[0]*unwrapH1[0]+fH[0]*unwrapH2[0];
          v[1] = fO[1]*unwrap[1]+fH[1]*unwrapH1[1]+fH[1]*unwrapH2[1];
          v[2] = fO[2]*unwrap[2]+fH[2]*unwrapH1[2]+fH[2]*unwrapH2[2];
          v[3] = fO[0]*unwrap[1]+fH[0]*unwrapH1[1]+fH[0]*unwrapH2[1];
          v[4] = fO[0]*unwrap[2]+fH[0]*unwrapH1[2]+fH[0]*unwrapH2[2];
          v[5] = fO[1]*unwrap[2]+fH[1]*unwrapH1[2]+fH[1]*unwrapH2[2];
          v_tally(i, v);
        }

      }

However, there seems something wrong in the results, because the water behaviors very differently comparing to those simulated with the original fix efield.

I have attached some necessary files.

I wonder what’s the problem and how can I fix it.

Any sugesstions are welcome.

Best regard,
Meng

dump.lammpstrj (2.1 MB)
fix_efield_tip4p.cpp (17.1 KB)
fix_efield_tip4p.h (2.9 KB)
AgEdown-0.05-rep.in (2.2 KB)
restart3.equil (763.1 KB)

2

If you write code on your own, you also have to debug it. That is part of the deal and your job.

For this kind of issue, it should be extremely straightforward to see if your computation is correct.

  • set up an input where you have only a small number of TIP4P water molecules (2-5) that you read from a data file. do not include any time integration or other settings. use pair style tip4p/cut and pair_modify compute no. force constants for both bond and angle should be 0. the only fix should be your custom fix. set up a custom dump that outputs atom ID, type, positions and forces and include a “run 0 post no” command as the last command in the input. this should contain only the forces from your fix in the dump output.
  • now set up a second input where you represent the positions for the O and M atom explicitly and put the charge on the M point and 0 on the O position. everything else would be the same except you would use fix efield and pair style coul/cut

now you can compare the results and the center of mass force for each molecule as well as the torque should be the same in both setups.

3

Thanks for the reply.

I am a littile confused about the modelling method of the real four site model. I have obtained the data file for the four site water model, and what should I do to maintain the relative position of each atom in the water molecule unchanged, should I use some fix rigid command? just like the ways mentioned in about water model TIP5P and six-site

Best regards,
Meng

4

For testing your fix you don’t need to move the atoms. So it doesn’t matter. Please make sure you read my suggestions carefully.

5

Thanks for your reply.

I will try that.

Sorry for bother you.

6

Dear axel,

I still have a small question :), and I am very sorry to bother you again.
In the second stage, I have generate a data file that containing the M point. Should I add bonds for the connection of M point and O point? If I should, how could I give the relative parameters, such as bond_coeff? And If there is a bond, I think there would also be impropers. I am a little confused about those.

I am really a newbie in LAMMPS, and lots of thing are not conversant. Sorry for the inconvenience.

Best regards,
Meng

7

Again, you didn’t pay attention to all details that I mentioned in my suggestion. For debugging you want to turn off any interactions but the forces (and stress) from the external field. Since the atoms are not supposed to be moving, it won’t matter, anyway.

8

Dear axel,

Thanks for your help, and I have found the actual reason that results the wrong compuation of the total force acting on each H and O.

After comparing the force results of fix efield, fix efild/tip4p, and without any efield. I found that the force acting on M point should be decomposed to two H atoms and one O atom in the water molecule; however, some of the H atoms’ force were not changed after the decomposion. I think that this is contributed to the parallel operation of LAMMPS. These H atoms might be on the different cores of its O atom. But after checked some documents, I didn’t find the way to resolve that. Could you please give me some suggestions?

Thanks for your help and sorry to bother you again.

Best regards,
Meng

9

As a general tip, try running your tests with differing numbers of processes (i.e. mpirun -np X with different X) and with different settings of newton on/off to detect if you have a problem in communication. This will help give some evidence for whether your code is treating “local” and “ghost” atoms correctly.

10

Ok, the problem is that some atoms, when the water molecule straddles a sub-domain boundary (or a periodic boundary, which is internally the same), are so-called ghost atoms, i.e. its local index is >= atom->nlocal.

To deal with that you have to do the following:

  • you must collect the forces in a separate array that needs to be dimensioned to [atom->nmax][3]
  • this array must be grown or reallocated, if needed before it is used, since the number or local and ghost atoms may change
  • this array must be zeroed in every step for local and ghost atoms before computing the efield force
  • the forces are then added to the atoms in the array, even if they are not local
  • you then must do a so-called reverse communication, where the forces from the ghost atoms are summed to their corresponding local atoms.
  • for that you need to set the comm_reverse flag in the fix constructor, implement pack/unpack_reverse_comm() and call comm->reverse_comm(this) in the fix after all efield forces have been computed.
  • after the reverse communication you can add the forces in the custom array to the atom->f array while looping only over local atoms (that are member of the fix group).

If you need some more inspiration, you can look at the fixes in the QEQ package. Since you don’t need to keep the per-atom data around for multiple time steps (you clear the array), you don’t have to do the more complex memory management with a callback done by the qeq fixes.

Overall, this is quite some complication, though that makes defining TIP4P water with an explicit point M and fix rigid rather attractive (where you can use the existing fix efield) unless you have to do very large systems or very long trajectories where you cannot afford to spend the extra computational effort.

For more info on the communication, please also see:
https://docs.lammps.org/Developer_par_comm.html
https://docs.lammps.org/Developer_comm_ops.html

11

Dear srtee,

Thanks for your suggestion, I have run the tests with different number of cores. Results show that the number of H with a wrong force is different, which means the problem is probably due to the unccorect consieration of ghost atoms.

Best regards,
Meng

12

Dear axel,

Thanks for your suggestion.

These procedures are really complex and hard for me, haha. Because I only considered constant electric field, in other words, the electric forces for all of the H atoms are similar. Therefore, I think I can use a simpler way. That is, neglecting which O atom that the specific H belongs to. Specifying the decomposed electric force at M point according to the atom type. I would try this method, and if you think this simplification might result in some errors, please feel free to point it out for me. Thanks again.

Best regards,
Meng

13

Before you go too far, you may want to compare your simulation speeds with the explicit four-point water model and with the tip4p style to see if all of this work is worthwhile.

If you spend (for example) a week writing out the efield modifications for tip4p to save (for example) two days of simulation time, then that doesn’t pay off – especially since there’s a chance you will have a subtle bug writing your efield modifications. There are other ways to speed up your simulation – for example Intel acceleration can be used with an explicit four-point water model, but not with tip4p styles (if I recall correctly).

14

Dear srtee,

Thanks for your reminding, I have’t compared the computing cost of the tip4p style and four-point water model yet. There are 1560 water molecules in my system. And I would check the speed before I go to the next step work.

Best regards,
Meng

15

I do think so. If you are not projecting out the forces, there is no meaning in writing a new fix, but you could just use the existing fix efield. The error will be about the same.

This is how LAMMPS is designed and why it scales to a very large number of processors because the effort remains nearly constant.

There is a simpler way, though, that doesn’t have such good parallel scaling and which will not work well with very large systems (it uses more memory and gets slower the more processors and atoms you use). The reason I didn’t mention it before is, that with this approach your code is not likely to be considered for inclusion into LAMMPS once it is done:

  • allocate a fieldforce array in the fix that has the dimensions [atom->natoms][3]
  • zero the array
  • compute the forces on each atom and store them in fieldforce[tag[i]][*]
  • perform an MPI_Allreduce() with MPI_SUM on &fieldforce[0][0] with length 3*atom->natoms
  • now loop over the local atoms and add the combined forces to the local atoms with
f[i][0] += fieldforce[tag[i]][0];
f[i][1] += fieldforce[tag[i]][1];
f[i][2] += fieldforce[tag[i]][2];

This will work because in case of forces on atoms that are ghost atoms, there will be non-zero force contributions in the fieldforce array on multiple processors and with the MPI_Allreduce() those arrays are summed by item across all MPI processes and then distributed again.

Going the path of the reverse communication looks far more complex than it actually is. The complicated stuff is all hidden from you in the functions calling the pack/unpack functions and the comm->reverse_comm() call.