Small changes to fix adapt for solid free energies

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@9150 f3b2605a-c512-4ea7-a41b-209d697bcdaa

src/compute_ti.cpp

@@ -16,6 +16,7 @@
 ------------------------------------------------------------------------- */
 
 #include "mpi.h"
+#include "atom.h"
 #include "string.h"
 #include "compute_ti.h"
 #include "update.h"
@@ -40,18 +41,23 @@ ComputeTI::ComputeTI(LAMMPS *lmp, int narg, char **arg) :
   if (narg < 4) error->all(FLERR,"Illegal compute ti command");
 
   peflag = 1;
+  peratom_flag = 1;
+  peatomflag = 1; 
   scalar_flag = 1;
   extscalar = 1;
   timeflag = 1;
 
   // terms come in triplets
+  // changed to quadruplets to include atom type
 
-  nterms = (narg-3) / 3;
+  nterms = (narg-3) / 4;
-  if (narg != 3*nterms + 3) error->all(FLERR,"Illegal compute ti command");
+  if (narg != 4*nterms + 3) error->all(FLERR,"Illegal compute ti command");
 
   which = new int[nterms];
   ivar1 = new int[nterms];
   ivar2 = new int[nterms];
+  ilo = new int[nterms];
+  ihi = new int[nterms];
   var1 = new char*[nterms];
   var2 = new char*[nterms];
   pptr = new Pair*[nterms];
@@ -65,15 +71,17 @@ ComputeTI::ComputeTI(LAMMPS *lmp, int narg, char **arg) :
 
   int iarg = 3;
   while (iarg < narg) {
-    if (iarg+3 > narg) error->all(FLERR,"Illegal compute ti command");
+    if (iarg+4 > narg) error->all(FLERR,"Illegal compute ti command");
     if (strcmp(arg[iarg],"kspace") == 0) which[nterms] = KSPACE;
     else if (strcmp(arg[iarg],"tail") == 0) which[nterms] = TAIL;
     else {
-      which[nterms] = PAIR;
+      which[nterms] = PAIR;} 
-      int n = strlen(arg[iarg]) + 1;
+    int n = strlen(arg[iarg]) + 1;
-      pstyle[nterms] = new char[n];
+    pstyle[nterms] = new char[n];
-      strcpy(pstyle[nterms],arg[iarg]);
+    strcpy(pstyle[nterms],arg[iarg]);
-    }
+    force->bounds(arg[iarg+1],atom->ntypes,ilo[nterms],ihi[nterms]);
+    
+    iarg += 1;
 
     if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) {
       int n = strlen(&arg[iarg+1][2]) + 1;
@@ -105,6 +113,8 @@ ComputeTI::~ComputeTI()
   delete [] ivar2;
   delete [] var1;
   delete [] var2;
+  delete [] ilo;
+  delete [] ihi;
   delete [] pptr;
   delete [] pstyle;
 }
@@ -153,22 +163,79 @@ double ComputeTI::compute_scalar()
   double dUdl = 0.0;
 
   for (int m = 0; m < nterms; m++) {
+    int total_flag = 0;
+    if ((ihi[m]-ilo[m])==atom->ntypes) total_flag = 1; 
+    eng = 0.0;
     value1 = input->variable->compute_equal(ivar1[m]);
     value2 = input->variable->compute_equal(ivar2[m]);
     if (value1 == 0.0) continue;
 
     if (which[m] == PAIR) {
-      eng = pptr[m]->eng_vdwl + pptr[m]->eng_coul;
+      int ntypes = atom->ntypes;
-      MPI_Allreduce(&eng,&engall,1,MPI_DOUBLE,MPI_SUM,world);
+      int *mask = atom->mask; 
+      if (total_flag) {
+        eng = pptr[m]->eng_vdwl + pptr[m]->eng_coul;
+        MPI_Allreduce(&eng,&engall,1,MPI_DOUBLE,MPI_SUM,world);
+      } 
+      else { 
+        int nlocal = atom->nlocal;
+        int npair = nlocal;
+        int *mask = atom->mask;
+        int *type = atom->type;
+
+        double *eatom = pptr[m]->eatom;
+
+        if (force->newton) npair += atom->nghost;
+        for(int i = 0; i < npair; i++)
+        {     
+          if ((ilo[m]<=type[i])&(ihi[m]>=type[i]))  
+            eng += eatom[i];  
+        }  
+        MPI_Allreduce(&eng,&engall,1,MPI_DOUBLE,MPI_SUM,world);
+      }
       dUdl += engall/value1 * value2;
 
     } else if (which[m] == TAIL) {
       double vol = domain->xprd*domain->yprd*domain->zprd;
-      eng = force->pair->etail / vol;
+      if (total_flag) 
+        eng = force->pair->etail / vol;
+      else {
+        eng = 0; 
+        for (int it = 1; it <= atom->ntypes; it++) {
+          int jt; 
+          if ((it >= ilo[m])&&(it <=ihi[m])) jt = it;
+          else jt = ilo[m];
+          for (; jt <=ihi[m];jt++) {
+            if ((force->pair->tail_flag)&&(force->pair->setflag[it][jt])) {
+              double cut = force->pair->init_one(it,jt);
+              eng += force->pair->etail_ij;
+            }
+            if (it !=jt) eng += force->pair->etail_ij; 
+          }
+        }
+        eng /= vol; 
+      }
       dUdl += eng/value1 * value2;
 
     } else if (which[m] == KSPACE) {
-      eng = force->kspace->energy;
+
+      int ntypes = atom->ntypes;
+      int *mask = atom->mask; 
+      if (total_flag) 
+        eng = force->kspace->energy;
+      else { 
+        int nlocal = atom->nlocal;
+        int npair = nlocal;
+        int *mask = atom->mask;
+        int *type = atom->type;
+        double *eatom = force->kspace->eatom;
+        eng = 0;
+        for(int i = 0; i < nlocal; i++)
+          if ((ilo[m]<=type[i])&(ihi[m]>=type[i])) 
+	    eng += eatom[i];  
+        MPI_Allreduce(&eng,&engall,1,MPI_DOUBLE,MPI_SUM,world);
+        eng = engall;
+      }
       dUdl += eng/value1 * value2;
     }
   }

src/fix_adapt.cpp

@@ -33,7 +33,7 @@ using namespace FixConst;
 using namespace MathConst;
 
 enum{PAIR,KSPACE,ATOM};
-enum{DIAMETER};
+enum{DIAMETER,CHARGE};
 
 /* ---------------------------------------------------------------------- */
 
@@ -71,6 +71,7 @@ FixAdapt::FixAdapt(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
 
   nadapt = 0;
   diamflag = 0;
+  chgflag = 0;
 
   iarg = 4;
   while (iarg < narg) {
@@ -110,6 +111,9 @@ FixAdapt::FixAdapt(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
       if (strcmp(arg[iarg+1],"diameter") == 0) {
         adapt[nadapt].aparam = DIAMETER;
         diamflag = 1;
+      } else if (strcmp(arg[iarg+1],"charge") == 0) {
+        adapt[nadapt].aparam = CHARGE; 
+        chgflag = 1; 
       } else error->all(FLERR,"Illegal fix adapt command");
       if (strstr(arg[iarg+2],"v_") == arg[iarg+2]) {
         int n = strlen(&arg[iarg+2][2]) + 1;
@@ -229,6 +233,10 @@ void FixAdapt::init()
         if (!atom->radius_flag)
           error->all(FLERR,"Fix adapt requires atom attribute diameter");
       }
+      if (ad->aparam == CHARGE) {
+	if (!atom->q_flag)
+	  error->all(FLERR,"Fix adapt requires atom attribute charge");
+      }
     }
   }
 
@@ -334,6 +342,16 @@ void FixAdapt::change_settings()
                 radius[i]*radius[i]*radius[i] * density;
             }
         }
+      } else if (ad->aparam == CHARGE) {
+        double *q = atom->q; 
+	int *mask = atom->mask;
+	int nlocal = atom->nlocal;
+	for (i = 0; i < nlocal; i++)
+	  if (mask[i] & groupbit) {
+//              cout<<"Fix adapt called now!\n";
+//            cout<<" Charge of atom: "<<i<<" is: "<<value; 
+            q[i] = value; 
+	  }
       }
     }
   }

src/fix_adapt.h

@@ -27,6 +27,7 @@ namespace LAMMPS_NS {
 class FixAdapt : public Fix {
  public:
   int diamflag;        // 1 if atom diameters will vary, for AtomVecGranular
+  int chgflag; 
 
   FixAdapt(class LAMMPS *, int, char **);
   ~FixAdapt();
@@ -99,4 +100,8 @@ E: Fix adapt requires atom attribute diameter
 
 The atom style being used does not specify an atom diameter.
 
+E: Fix adapt requires atom attribute charge 
+
+The atom style being used does not specify an atom charge.
+
 */

src/pair_gauss.cpp

@@ -217,7 +217,11 @@ void PairGauss::coeff(int narg, char **arg)
 
 double PairGauss::init_one(int i, int j)
 {
-  if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
+ 
+  // This error is triggered when ti is performed on lj/cut tail
+  // in presence of extra atom type for tether sites
+  // "i = 2 j = 1 ERROR: All pair coeffs are not set (pair_gauss.cpp:223)"
+  //  if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
 
   if (offset_flag) offset[i][j] = a[i][j]*exp(-b[i][j]*cut[i][j]*cut[i][j]);
   else offset[i][j] = 0.0;