implement external efield in fix/dplr

Author: Han Wang

source/lmp/fix_dplr.cpp

@@ -34,7 +34,7 @@ is_key (const string& input)
 
 
 FixDPLR::FixDPLR(LAMMPS *lmp, int narg, char **arg) 
-    :Fix(lmp, narg, arg)
+    :Fix(lmp, narg, arg), efield(3, 0.0)
 {
   virial_flag = 1;
 
@@ -54,6 +54,13 @@ FixDPLR::FixDPLR(LAMMPS *lmp, int narg, char **arg)
       model = string(arg[iarg+1]);
       iarg += 2;
     }
+    if (string(arg[iarg]) == string("efield")) {
+      if (iarg+3 > narg) error->all(FLERR,"Illegal fix adapt command, efield should be provided 3 float numbers");
+      efield[0] = atof(arg[iarg+1]);
+      efield[1] = atof(arg[iarg+2]);
+      efield[2] = atof(arg[iarg+3]);
+      iarg += 4;
+    }
     if (string(arg[iarg]) == string("type_associate")) {
       int iend = iarg+1;
       while (iend < narg && (! is_key(arg[iend]) )) {
@@ -347,6 +354,7 @@ void FixDPLR::post_force(int vflag)
   int nall = nlocal + nghost;
   vector<FLOAT_PREC> dcoord(nall*3, 0.0), dbox(9, 0.0), dfele(nlocal*3, 0.0);
   vector<int> dtype(nall, 0);
+  // set values for dcoord, dbox, dfele
   {
     int *type = atom->type;
     for (int ii = 0; ii < nall; ++ii){
@@ -366,9 +374,30 @@ void FixDPLR::post_force(int vflag)
       }
     }
     assert(dfele_.size() == nlocal * 3);
+    // revise force according to efield
     for (int ii = 0; ii < nlocal*3; ++ii){
       dfele[ii] = dfele_[ii];
     }
+    // revise force and virial according to efield
+    double * q = atom->q;
+    imageint *image = atom->image;
+    double unwrap[3];
+    double v[6];
+    for (int ii = 0; ii < nlocal; ++ii){
+      for (int dd = 0; dd < 3; ++dd){
+	dfele[ii*3+dd] += q[ii] * efield[dd];
+      }
+      domain->unmap(x[ii],image[ii],unwrap);
+      if (evflag) {
+	v[0] = q[ii] * efield[0] *unwrap[0];
+	v[1] = q[ii] * efield[1] *unwrap[1];
+	v[2] = q[ii] * efield[2] *unwrap[2];
+	v[3] = q[ii] * efield[0] *unwrap[1];
+	v[4] = q[ii] * efield[0] *unwrap[2];
+	v[5] = q[ii] * efield[1] *unwrap[2];
+	v_tally(ii, v);
+      }
+    }    
   }
   // lmp nlist
   NeighList * list = pair_nnp->list;

source/lmp/fix_dplr.h

@@ -45,6 +45,7 @@ namespace LAMMPS_NS {
     map<int,int > bk_type_asso;
     vector<FLOAT_PREC> dipole_recd;
     vector<double> dfcorr_buff;
+    vector<double> efield;
     void get_valid_pairs(vector<pair<int,int> >& pairs);
   };
 }