Merge pull request #226 from akohlmey/pair-gauss-mixing

add mixing support for pair styles gauss and gauss/cut

doc/src/pair_gauss.txt

@@ -106,8 +106,31 @@ more instructions on how to use the accelerated styles effectively.
 
 [Mixing, shift, table, tail correction, restart, rRESPA info]:
 
-These pair styles do not support mixing. Thus, coefficients for all
-I,J pairs must be specified explicitly.
+For atom type pairs I,J and I != J, the A, B, H, sigma_h, r_mh
+parameters, and the cutoff distance for these pair styles can be mixed:
+A (energy units)
+sqrt(1/B) (distance units, see below)
+H (energy units)
+sigma_h (distance units)
+r_mh (distance units)
+cutoff (distance units):ul
+
+The default mix value is {geometric}.
+Only {arithmetic} and {geometric} mix values are supported.
+See the "pair_modify" command for details.
+
+The A and H parameters are mixed using the same rules normally
+used to mix the "epsilon" parameter in a Lennard Jones interaction.
+The sigma_h, r_mh, and the cutoff distance are mixed using the same
+rules used to mix the "sigma" parameter in a Lennard Jones interaction.
+The B parameter is converted to a distance (sigma), before mixing
+(using sigma=B^-0.5), and converted back to a coefficient
+afterwards (using B=sigma^2).
+Negative A values are converted to positive A values (using abs(A))
+before mixing, and converted back after mixing
+(by multiplying by sign(Ai)*sign(Aj)).
+This way, if either particle is repulsive (if Ai<0 or Aj<0),
+then the default interaction between both particles will be repulsive.
 
 The {gauss} style does not support the "pair_modify"_pair_modify.html
 shift option. There is no effect due to the Gaussian well beyond the

src/USER-MISC/pair_gauss_cut.cpp

@@ -39,6 +39,7 @@ using namespace MathConst;
 PairGaussCut::PairGaussCut(LAMMPS *lmp) : Pair(lmp)
 {
   respa_enable = 0;
+  writedata = 1;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -96,7 +97,6 @@ void PairGaussCut::compute(int eflag, int vflag)
 
     for (jj = 0; jj < jnum; jj++) {
       j = jlist[jj];
-
       factor_lj = special_lj[sbmask(j)];
       j &= NEIGHMASK;
 
@@ -221,7 +221,21 @@ void PairGaussCut::coeff(int narg, char **arg)
 
 double PairGaussCut::init_one(int i, int j)
 {
-  if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
+  if (setflag[i][j] == 0) {
+    hgauss[i][j] = mix_energy(fabs(hgauss[i][i]), fabs(hgauss[j][j]),
+			      fabs(sigmah[i][i]), fabs(sigmah[j][j]));
+
+    // If either of the particles is repulsive (ie, if hgauss > 0),
+    // then the interaction between both is repulsive.
+    double sign_hi = (hgauss[i][i] >= 0.0) ? 1.0 : -1.0;
+    double sign_hj = (hgauss[j][j] >= 0.0) ? 1.0 : -1.0;
+    hgauss[i][j] *= MAX(sign_hi, sign_hj);
+
+    sigmah[i][j] = mix_distance(sigmah[i][i], sigmah[j][j]);
+    rmh[i][j] = mix_distance(rmh[i][i], rmh[j][j]);
+    cut[i][j] = mix_distance(cut[i][i], cut[j][j]);
+  }
+
   pgauss[i][j] = hgauss[i][j] / sqrt(MY_2PI) / sigmah[i][j];
 
   if (offset_flag) {
@@ -334,6 +348,27 @@ void PairGaussCut::read_restart_settings(FILE *fp)
   MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
 }
 
+/* ----------------------------------------------------------------------
+   proc 0 writes to data file
+------------------------------------------------------------------------- */
+
+void PairGaussCut::write_data(FILE *fp)
+{
+  for (int i = 1; i <= atom->ntypes; i++)
+    fprintf(fp,"%d %g %g %g\n",i,hgauss[i][i],rmh[i][i],sigmah[i][i]);
+}
+
+/* ----------------------------------------------------------------------
+   proc 0 writes all pairs to data file
+------------------------------------------------------------------------- */
+
+void PairGaussCut::write_data_all(FILE *fp)
+{
+  for (int i = 1; i <= atom->ntypes; i++)
+    for (int j = i; j <= atom->ntypes; j++)
+      fprintf(fp,"%d %d %g %g %g %g\n",i,j,hgauss[i][j],rmh[i][j],sigmah[i][j],cut[i][j]);
+}
+
 /* ---------------------------------------------------------------------- */
 
 double PairGaussCut::single(int i, int j, int itype, int jtype, double rsq,

src/USER-MISC/pair_gauss_cut.h

@@ -42,6 +42,8 @@ class PairGaussCut : public Pair {
   virtual void read_restart(FILE *);
   virtual void write_restart_settings(FILE *);
   virtual void read_restart_settings(FILE *);
+  virtual void write_data(FILE *fp);
+  virtual void write_data_all(FILE *fp);
 
   virtual double memory_usage();
 

src/pair_gauss.cpp

@@ -34,10 +34,11 @@ using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */
 
-PairGauss::PairGauss(LAMMPS *lmp) :Pair(lmp)
+PairGauss::PairGauss(LAMMPS *lmp) : Pair(lmp)
 {
   nextra = 1;
   pvector = new double[1];
+  writedata = 1;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -187,7 +188,7 @@ void PairGauss::coeff(int narg, char **arg)
     error->all(FLERR,"Incorrect args for pair coefficients");
   if (!allocated) allocate();
 
-  int ilo, ihi, jlo, jhi;
+  int ilo,ihi,jlo,jhi;
   force->bounds(arg[0],atom->ntypes,ilo,ihi);
   force->bounds(arg[1],atom->ntypes,jlo,jhi);
 
@@ -204,7 +205,7 @@ void PairGauss::coeff(int narg, char **arg)
       b[i][j] = b_one;
       cut[i][j] = cut_one;
       setflag[i][j] = 1;
-      count++ ;
+      count++;
     }
   }
 
@@ -217,12 +218,27 @@ void PairGauss::coeff(int narg, char **arg)
 
 double PairGauss::init_one(int i, int j)
 {
+  if (setflag[i][j] == 0) {
+    double sign_bi = (b[i][i] >= 0.0) ? 1.0 : -1.0;
+    double sign_bj = (b[j][j] >= 0.0) ? 1.0 : -1.0;
+    double si = sqrt(0.5/fabs(b[i][i]));
+    double sj = sqrt(0.5/fabs(b[j][j]));
+    double sij = mix_distance(si, sj);
+    b[i][j] = 0.5 / (sij*sij);
+    b[i][j] *= MAX(sign_bi, sign_bj);
 
-  // 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");
+    // Negative "a" values are useful for simulating repulsive particles.
+    // If either of the particles is repulsive (a<0), then the
+    // interaction between both is repulsive.
+    double sign_ai = (a[i][i] >= 0.0) ? 1.0 : -1.0;
+    double sign_aj = (a[j][j] >= 0.0) ? 1.0 : -1.0;
+    a[i][j] = mix_energy(fabs(a[i][i]), fabs(a[j][j]), si, sj);
+    a[i][j] *= MIN(sign_ai, sign_aj);
 
+    cut[i][j] = mix_distance(cut[i][i],cut[j][j]);
+  }
+
+  // cutoff correction to energy
   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;
 
@@ -260,7 +276,6 @@ void PairGauss::write_restart(FILE *fp)
 void PairGauss::read_restart(FILE *fp)
 {
   read_restart_settings(fp);
-
   allocate();
 
   int i,j;
@@ -309,6 +324,27 @@ void PairGauss::read_restart_settings(FILE *fp)
   MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
 }
 
+/* ----------------------------------------------------------------------
+   proc 0 writes to data file
+------------------------------------------------------------------------- */
+
+void PairGauss::write_data(FILE *fp)
+{
+  for (int i = 1; i <= atom->ntypes; i++)
+    fprintf(fp,"%d %g %g\n",i,a[i][i],b[i][i]);
+}
+
+/* ----------------------------------------------------------------------
+   proc 0 writes all pairs to data file
+------------------------------------------------------------------------- */
+
+void PairGauss::write_data_all(FILE *fp)
+{
+  for (int i = 1; i <= atom->ntypes; i++)
+    for (int j = i; j <= atom->ntypes; j++)
+      fprintf(fp,"%d %d %g %g %g\n",i,j,a[i][j],b[i][j],cut[i][j]);
+}
+
 /* ---------------------------------------------------------------------- */
 
 double PairGauss::single(int i, int j, int itype, int jtype, double rsq,

src/pair_gauss.h

@@ -36,6 +36,8 @@ class PairGauss : public Pair {
   void read_restart(FILE *);
   void write_restart_settings(FILE *);
   void read_restart_settings(FILE *);
+  void write_data(FILE *fp);
+  void write_data_all(FILE *fp);
   double single(int, int, int, int, double, double, double, double &);
   void *extract(const char *, int &);