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

src/KSPACE/ewald.cpp

@@ -439,7 +439,7 @@ void Ewald::compute(int eflag, int vflag)
       energy += ug[k] * (sfacrl_all[k]*sfacrl_all[k] +
                          sfacim_all[k]*sfacim_all[k]);
 
-    if (atom->natoms != natoms_original) {
+    if (qsum_update_flag || (atom->natoms != natoms_original)) {
       qsum_qsq(0);
       natoms_original = atom->natoms;
     }

src/KSPACE/msm.cpp

@@ -576,7 +576,7 @@ void MSM::compute(int eflag, int vflag)
     MPI_Allreduce(&energy,&energy_all,1,MPI_DOUBLE,MPI_SUM,world);
     energy = energy_all;
 
-    if (atom->natoms != natoms_original) {
+    if (qsum_update_flag || (atom->natoms != natoms_original)) {
       qsum_qsq(0);
       natoms_original = atom->natoms;
     }

src/KSPACE/pppm.cpp

@@ -677,7 +677,7 @@ void PPPM::compute(int eflag, int vflag)
     MPI_Allreduce(&energy,&energy_all,1,MPI_DOUBLE,MPI_SUM,world);
     energy = energy_all;
 
-    if (atom->natoms != natoms_original) {
+    if (qsum_update_flag || (atom->natoms != natoms_original)) {
       qsum_qsq(0);
       natoms_original = atom->natoms;
     }

src/KSPACE/pppm_cg.cpp

@@ -218,7 +218,7 @@ void PPPMCG::compute(int eflag, int vflag)
     MPI_Allreduce(&energy,&energy_all,1,MPI_DOUBLE,MPI_SUM,world);
     energy = energy_all;
 
-    if (atom->natoms != natoms_original) {
+    if (qsum_update_flag || (atom->natoms != natoms_original)) {
       qsum_qsq(0);
       natoms_original = atom->natoms;
     }

src/KSPACE/pppm_disp.cpp

@@ -1135,7 +1135,7 @@ void PPPMDisp::compute(int eflag, int vflag)
     energy_1 *= 0.5*volume;
     energy_6 *= 0.5*volume;
     
-    if (atom->natoms != natoms_original) {
+    if (qsum_update_flag || (atom->natoms != natoms_original)) {
       qsum_qsq(0);
       natoms_original = atom->natoms;
     }

src/KSPACE/pppm_stagger.cpp

@@ -212,7 +212,7 @@ void PPPMStagger::compute(int eflag, int vflag)
     MPI_Allreduce(&energy,&energy_all,1,MPI_DOUBLE,MPI_SUM,world);
     energy = energy_all;
 
-    if (atom->natoms != natoms_original) {
+    if (qsum_update_flag || (atom->natoms != natoms_original)) {
       qsum_qsq(0);
       natoms_original = atom->natoms;
     }

src/QEQ/fix_qeq.cpp

@@ -262,6 +262,7 @@ void FixQEq::reallocate_matrix()
 void FixQEq::init_list(int id, NeighList *ptr)
 {
   list = ptr;
+  if (force->kspace) force->kspace->qsum_update_flag = 1;
 }
 
 /* ---------------------------------------------------------------------- */

src/QEQ/fix_qeq_dynamic.cpp

@@ -30,6 +30,7 @@
 #include "force.h"
 #include "group.h"
 #include "pair.h"
+#include "kspace.h"
 #include "respa.h"
 #include "memory.h"
 #include "error.h"
@@ -167,6 +168,7 @@ void FixQEqDynamic::pre_force(int vflag)
     }
   }
 
+  if (force->kspace) force->kspace->setup();
 }
 
 /* ---------------------------------------------------------------------- */

src/QEQ/fix_qeq_point.cpp

@@ -28,6 +28,7 @@
 #include "neigh_request.h"
 #include "update.h"
 #include "force.h"
+#include "kspace.h"
 #include "group.h"
 #include "respa.h"
 #include "memory.h"
@@ -82,6 +83,8 @@ void FixQEqPoint::pre_force(int vflag)
   matvecs += CG(b_t, t); 	// CG on t - parallel
   calculate_Q();
 
+  if (force->kspace) force->kspace->setup();
+
 }
 
 /* ---------------------------------------------------------------------- */

src/QEQ/fix_qeq_shielded.cpp

@@ -28,6 +28,7 @@
 #include "neigh_request.h"
 #include "update.h"
 #include "force.h"
+#include "kspace.h"
 #include "group.h"
 #include "respa.h"
 #include "memory.h"
@@ -125,6 +126,8 @@ void FixQEqShielded::pre_force(int vflag)
   matvecs += CG(b_t, t); 	// CG on t - parallel
   calculate_Q();
 
+  if (force->kspace) force->kspace->setup();
+
 }
 
 /* ---------------------------------------------------------------------- */

src/QEQ/fix_qeq_slater.cpp

@@ -28,9 +28,9 @@
 #include "neigh_request.h"
 #include "update.h"
 #include "force.h"
-#include "kspace.h"
 #include "group.h"
 #include "pair.h"
+#include "kspace.h"
 #include "respa.h"
 #include "math_const.h"
 #include "memory.h"
@@ -101,6 +101,8 @@ void FixQEqSlater::pre_force(int vflag)
   matvecs += CG(b_t, t); 	// CG on t - parallel
   calculate_Q();
 
+  if (force->kspace) force->kspace->setup();
+
 }
 
 /* ---------------------------------------------------------------------- */

src/USER-FEP/fix_adapt_fep.cpp

@@ -20,13 +20,14 @@
 #include "stdlib.h"
 #include "fix_adapt_fep.h"
 #include "atom.h"
-#include "comm.h"
 #include "update.h"
+#include "group.h"
 #include "modify.h"
 #include "force.h"
 #include "pair.h"
 #include "pair_hybrid.h"
 #include "kspace.h"
+#include "fix_store.h"
 #include "input.h"
 #include "variable.h"
 #include "math_const.h"
@@ -40,8 +41,6 @@ using namespace MathConst;
 enum{PAIR,KSPACE,ATOM};
 enum{DIAMETER,CHARGE};
 
-#undef ADAPT_DEBUG
-
 /* ---------------------------------------------------------------------- */
 
 FixAdaptFEP::FixAdaptFEP(LAMMPS *lmp, int narg, char **arg) :
@@ -51,6 +50,8 @@ FixAdaptFEP::FixAdaptFEP(LAMMPS *lmp, int narg, char **arg) :
   nevery = force->inumeric(FLERR,arg[3]);
   if (nevery < 0) error->all(FLERR,"Illegal fix adapt/fep command");
 
+  dynamic_group_allow = 1;
+
   // count # of adaptations
 
   nadapt = 0;
@@ -169,6 +170,8 @@ FixAdaptFEP::FixAdaptFEP(LAMMPS *lmp, int narg, char **arg) :
   for (int m = 0; m < nadapt; m++)
     if (adapt[m].which == PAIR)
       memory->create(adapt[m].array_orig,n+1,n+1,"adapt:array_orig");
+
+  id_fix_diam = id_fix_chg = NULL;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -184,6 +187,13 @@ FixAdaptFEP::~FixAdaptFEP()
     }
   }
   delete [] adapt;
+
+  // check nfix in case all fixes have already been deleted
+
+  if (id_fix_diam && modify->nfix) modify->delete_fix(id_fix_diam);
+  if (id_fix_chg && modify->nfix) modify->delete_fix(id_fix_chg);
+  delete [] id_fix_diam;
+  delete [] id_fix_chg;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -196,12 +206,92 @@ int FixAdaptFEP::setmask()
   return mask;
 }
 
+/* ----------------------------------------------------------------------
+   if need to restore per-atom quantities, create new fix STORE styles
+------------------------------------------------------------------------- */
+
+void FixAdaptFEP::post_constructor()
+{
+  if (!resetflag) return;
+  if (!diamflag && !chgflag) return;
+
+  // new id = fix-ID + FIX_STORE_ATTRIBUTE
+  // new fix group = group for this fix
+
+  id_fix_diam = NULL;
+  id_fix_chg = NULL;
+
+  char **newarg = new char*[5];
+  newarg[1] = group->names[igroup];
+  newarg[2] = (char *) "STORE";
+  newarg[3] = (char *) "1";
+  newarg[4] = (char *) "1";
+
+  if (diamflag) {
+    int n = strlen(id) + strlen("_FIX_STORE_DIAM") + 1;
+    id_fix_diam = new char[n];
+    strcpy(id_fix_diam,id);
+    strcat(id_fix_diam,"_FIX_STORE_DIAM");
+    newarg[0] = id_fix_diam;
+    modify->add_fix(5,newarg);
+    fix_diam = (FixStore *) modify->fix[modify->nfix-1];
+
+    if (fix_diam->restart_reset) fix_diam->restart_reset = 0;
+    else {
+      double *vec = fix_diam->vstore;
+      double *radius = atom->radius;
+      int *mask = atom->mask;
+      int nlocal = atom->nlocal;
+
+      for (int i = 0; i < nlocal; i++) {
+        if (mask[i] & groupbit) vec[i] = radius[i];
+        else vec[i] = 0.0;
+      }
+    }
+  }
+
+  if (chgflag) {
+    int n = strlen(id) + strlen("_FIX_STORE_CHG") + 1;
+    id_fix_chg = new char[n];
+    strcpy(id_fix_chg,id);
+    strcat(id_fix_chg,"_FIX_STORE_CHG");
+    newarg[0] = id_fix_chg;
+    modify->add_fix(5,newarg);
+    fix_chg = (FixStore *) modify->fix[modify->nfix-1];
+
+    if (fix_chg->restart_reset) fix_chg->restart_reset = 0;
+    else {
+      double *vec = fix_chg->vstore;
+      double *q = atom->q;
+      int *mask = atom->mask;
+      int nlocal = atom->nlocal;
+
+      for (int i = 0; i < nlocal; i++) {
+        if (mask[i] & groupbit) vec[i] = q[i];
+        else vec[i] = 0.0;
+      }
+    }
+  }
+
+  delete [] newarg;
+}
+
 /* ---------------------------------------------------------------------- */
 
 void FixAdaptFEP::init()
 {
   int i,j;
 
+  // allow a dynamic group only if ATOM attribute not used
+
+  if (group->dynamic[igroup])
+    for (int i = 0; i < nadapt; i++)
+      if (adapt[i].which == ATOM) 
+        error->all(FLERR,"Cannot use dynamic group with fix adapt/fep atom");
+
+  // when using kspace, we need to recompute some additional parameters in kspace->setup()
+  if (force->kspace) force->kspace->qsum_update_flag = 1;
+  
   // setup and error checks
 
   anypair = 0;
@@ -217,9 +307,18 @@ void FixAdaptFEP::init()
 
     if (ad->which == PAIR) {
       anypair = 1;
+      Pair *pair = NULL;
 
-      Pair *pair = force->pair_match(ad->pstyle,1);
+      if (lmp->suffix_enable) {
-      if (pair == NULL) error->all(FLERR,"Fix adapt/fep pair style does not exist");
+        char psuffix[128];
+        strcpy(psuffix,ad->pstyle);
+        strcat(psuffix,"/");
+        strcat(psuffix,lmp->suffix);
+        pair = force->pair_match(psuffix,1);
+      }
+      if (pair == NULL) pair = force->pair_match(ad->pstyle,1);
+      if (pair == NULL)
+        error->all(FLERR, "Fix adapt/fep pair style does not exist");
       void *ptr = pair->extract(ad->pparam,ad->pdim);
       if (ptr == NULL) 
         error->all(FLERR,"Fix adapt/fep pair style param not supported");
@@ -268,22 +367,18 @@ void FixAdaptFEP::init()
     }
   }
 
-#ifdef ADAPT_DEBUG
+  // fixes that store initial per-atom values
-  if (comm->me == 0 && screen) {
-    for (int m = 0; m < nadapt; m++) {
-      Adapt *ad = &adapt[m];
-      if (ad->which == PAIR && ad->pdim == 2) {
-        fprintf(screen, "###ADAPT original %s %s\n", ad->pstyle, ad->pparam);
-        fprintf(screen, "###ADAPT  I  J   old_param\n");
-        for (i = ad->ilo; i <= ad->ihi; i++)
-          for (j = MAX(ad->jlo,i); j <= ad->jhi; j++)
-            fprintf(screen, "###ADAPT %2d %2d %9.5f\n", i, j,
-                    ad->array_orig[i][j]);
-      }
-    }
-  }
-#endif
   
+  if (id_fix_diam) {
+    int ifix = modify->find_fix(id_fix_diam);
+    if (ifix < 0) error->all(FLERR,"Could not find fix adapt storage fix ID");
+    fix_diam = (FixStore *) modify->fix[ifix];
+  }
+  if (id_fix_chg) {
+    int ifix = modify->find_fix(id_fix_chg);
+    if (ifix < 0) error->all(FLERR,"Could not find fix adapt storage fix ID");
+    fix_chg = (FixStore *) modify->fix[ifix];
+  }
 }
 
 /* ---------------------------------------------------------------------- */
@@ -349,17 +444,6 @@ void FixAdaptFEP::change_settings()
           for (i = ad->ilo; i <= ad->ihi; i++)
             for (j = MAX(ad->jlo,i); j <= ad->jhi; j++)
               ad->array[i][j] = value;
-
-#ifdef ADAPT_DEBUG
-        if (comm->me == 0 && screen) {
-          fprintf(screen, "###ADAPT change %s %s\n", ad->pstyle, ad->pparam);
-          fprintf(screen, "###ADAPT  I  J   param\n");
-          for (i = ad->ilo; i <= ad->ihi; i++)
-            for (j = MAX(ad->jlo,i); j <= ad->jhi; j++)
-              fprintf(screen, "###ADAPT %2d %2d %9.5f\n", i, j, ad->array[i][j]);
-        }
-#endif
-
       }
 
     // set kspace scale factor
@@ -369,7 +453,7 @@ void FixAdaptFEP::change_settings()
 
     } else if (ad->which == ATOM) {
 
-      // set radius from diameter
+      // reset radius from diameter
       // also scale rmass to new value
 
       if (ad->aparam == DIAMETER) {
@@ -381,15 +465,16 @@ void FixAdaptFEP::change_settings()
         double *radius = atom->radius;
         double *rmass = atom->rmass;
         int *mask = atom->mask;
-        int natom = atom->nlocal + atom->nghost;
+        int nlocal = atom->nlocal;
+        int nall = nlocal + atom->nghost;
 
         if (mflag == 0) {
-          for (i = 0; i < natom; i++)
+          for (i = 0; i < nall; i++)
             if (atype[i] >= ad->ilo && atype[i] <= ad->ihi)
               if (mask[i] & groupbit)
                 radius[i] = 0.5*value;
         } else {
-          for (i = 0; i < natom; i++)
+          for (i = 0; i < nall; i++)
             if (atype[i] >= ad->ilo && atype[i] <= ad->ihi)
               if (mask[i] & groupbit) {
                 density = rmass[i] / (4.0*MY_PI/3.0 *
@@ -403,22 +488,12 @@ void FixAdaptFEP::change_settings()
         int *atype = atom->type;
         double *q = atom->q; 
         int *mask = atom->mask;
-        int natom = atom->nlocal + atom->nghost;
+        int nlocal = atom->nlocal;
+        int nall = nlocal + atom->nghost;
 
-        for (i = 0; i < natom; i++)
+        for (i = 0; i < nall; i++)
           if (atype[i] >= ad->ilo && atype[i] <= ad->ihi)
             if (mask[i] & groupbit) q[i] = value;
-
-#ifdef ADAPT_DEBUG
-        if (comm->me == 0 && screen) {
-          fprintf(screen, "###ADAPT change charge\n");
-          fprintf(screen, "###ADAPT  atom  I   q\n");
-          for (i = 0; i < atom->nlocal; i++)
-            if (atype[i] >= ad->ilo && atype[i] <= ad->ihi)
-              if (mask[i] & groupbit)
-                fprintf(screen, "###ADAPT %5d %2d %9.5f\n", i, atype[i], q[i]);
-        }
-#endif
       }
     }
   }
@@ -430,6 +505,10 @@ void FixAdaptFEP::change_settings()
   // and also offset and tail corrections
 
   if (anypair) force->pair->reinit();
+
+  if (force->kspace)
+    force->kspace->setup();
+
 }
 
 /* ----------------------------------------------------------------------
@@ -452,9 +531,36 @@ void FixAdaptFEP::restore_settings()
       *kspace_scale = 1.0;
 
     } else if (ad->which == ATOM) {
+      if (diamflag) {
+        double density;
 
+        double *vec = fix_diam->vstore;
+        double *radius = atom->radius;
+        double *rmass = atom->rmass;
+        int *mask = atom->mask;
+        int nlocal = atom->nlocal;
+
+        for (int i = 0; i < nlocal; i++)
+          if (mask[i] & groupbit) {
+            density = rmass[i] / (4.0*MY_PI/3.0 *
+                                  radius[i]*radius[i]*radius[i]);
+            radius[i] = vec[i];
+            rmass[i] = 4.0*MY_PI/3.0 * radius[i]*radius[i]*radius[i] * density;
+          }
+      }
+      if (chgflag) {
+        double *vec = fix_chg->vstore;
+        double *q = atom->q;
+        int *mask = atom->mask;
+        int nlocal = atom->nlocal;
+
+        for (int i = 0; i < nlocal; i++)
+          if (mask[i] & groupbit) q[i] = vec[i];
+      }
     }
   }
 
   if (anypair) force->pair->reinit();
+
+  if (force->kspace) force->kspace->setup();
 }

src/USER-FEP/fix_adapt_fep.h

@@ -32,6 +32,7 @@ class FixAdaptFEP : public Fix {
   FixAdaptFEP(class LAMMPS *, int, char **);
   ~FixAdaptFEP();
   int setmask();
+  void post_constructor();
   void init();
   void setup_pre_force(int);
   void pre_force(int);
@@ -40,6 +41,8 @@ class FixAdaptFEP : public Fix {
  private:
   int nadapt,resetflag,scaleflag,afterflag;
   int anypair;
+  char *id_fix_diam,*id_fix_chg;
+  class FixStore *fix_diam,*fix_chg;
 
   struct Adapt {
     int which,ivar;

src/kspace.cpp

@@ -68,6 +68,8 @@ KSpace::KSpace(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
   suffix_flag = Suffix::NONE;
   adjust_cutoff_flag = 1;
   scalar_pressure_flag = 0;
+  qsum_update_flag = 0;
+  warn_neutral = 1;
 
   accuracy_absolute = -1.0;
   accuracy_real_6 = -1.0;
@@ -284,8 +286,11 @@ void KSpace::qsum_qsq(int flag)
   if (fabs(qsum) > SMALL) {
     char str[128];
     sprintf(str,"System is not charge neutral, net charge = %g",qsum);
+    if (warn_neutral && (comm->me == 0)) {
       if (flag) error->all(FLERR,str);
-    else if (comm->me == 0) error->warning(FLERR,str);
+      else error->warning(FLERR,str);
+    }
+    warn_neutral = 0;
   }
 }
 

src/kspace.h

@@ -51,6 +51,7 @@ class KSpace : protected Pointers {
                                  // for LJ coefficients
   int slabflag;
   int scalar_pressure_flag;      // 1 if using MSM fast scalar pressure
+  int qsum_update_flag;          // 1 if setup() needs to call qsum_qsq()
   double slab_volfactor;
 
 
@@ -168,6 +169,7 @@ class KSpace : protected Pointers {
   int minorder,overlap_allowed;
   int adjust_cutoff_flag;
   int suffix_flag;                  // suffix compatibility flag
+  int warn_neutral;                 // warn about non-neutral system if 1
   bigint natoms_original;
   double scale,qqrd2e;
   double qsum,qsqsum,q2;