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)