Unverified Commit f18d0507 authored by Axel Kohlmeyer's avatar Axel Kohlmeyer
Browse files

fix issues with lj/long pair styles when not using long-range for dispersion

parent c8f2634b
Loading
Loading
Loading
Loading
+12 −13
Original line number Diff line number Diff line
@@ -90,6 +90,8 @@ void PairBuckLongCoulLong::settings(int narg, char **arg)
    error->warning(FLERR,"Using largest cutoff for buck/long/coul/long");
  if (!*(++arg))
    error->all(FLERR,"Cutoffs missing in pair_style buck/long/coul/long");
  if (!((ewald_order^ewald_off) & (1<<6)))
    dispersionflag = 0;
  if (ewald_off & (1<<6))
    error->all(FLERR,"LJ6 off not supported in pair_style buck/long/coul/long");
  if (!((ewald_order^ewald_off) & (1<<1)))
@@ -387,6 +389,7 @@ void PairBuckLongCoulLong::write_restart_settings(FILE *fp)
  fwrite(&ncoultablebits,sizeof(int),1,fp);
  fwrite(&tabinner,sizeof(double),1,fp);
  fwrite(&ewald_order,sizeof(int),1,fp);
  fwrite(&dispersionflag,sizeof(int),1,fp);
}

/* ----------------------------------------------------------------------
@@ -403,6 +406,7 @@ void PairBuckLongCoulLong::read_restart_settings(FILE *fp)
    utils::sfread(FLERR,&ncoultablebits,sizeof(int),1,fp,NULL,error);
    utils::sfread(FLERR,&tabinner,sizeof(double),1,fp,NULL,error);
    utils::sfread(FLERR,&ewald_order,sizeof(int),1,fp,NULL,error);
    utils::sfread(FLERR,&dispersionflag,sizeof(int),1,fp,NULL,error);
  }
  MPI_Bcast(&cut_buck_global,1,MPI_DOUBLE,0,world);
  MPI_Bcast(&cut_coul,1,MPI_DOUBLE,0,world);
@@ -411,6 +415,7 @@ void PairBuckLongCoulLong::read_restart_settings(FILE *fp)
  MPI_Bcast(&ncoultablebits,1,MPI_INT,0,world);
  MPI_Bcast(&tabinner,1,MPI_DOUBLE,0,world);
  MPI_Bcast(&ewald_order,1,MPI_INT,0,world);
  MPI_Bcast(&dispersionflag,1,MPI_INT,0,world);
}

/* ----------------------------------------------------------------------
@@ -529,8 +534,7 @@ void PairBuckLongCoulLong::compute(int eflag, int vflag)
            if (eflag) ecoul = qiqj*(etable[k]+f*detable[k]-t.f);
          }
        }
      }
      else force_coul = ecoul = 0.0;
      } else force_coul = ecoul = 0.0;

      if (rsq < cut_bucksqi[typej]) {                        // buckingham
        double rn = r2inv*r2inv*r2inv,
@@ -543,16 +547,14 @@ void PairBuckLongCoulLong::compute(int eflag, int vflag)
              force_buck =
                r*expr*buck1i[typej]-g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq;
              if (eflag) evdwl = expr*buckai[typej]-g6*((a2+1.0)*a2+0.5)*x2;
            }
            else {                                        // special case
            } else {                                        // special case
              double f = special_lj[ni], t = rn*(1.0-f);
              force_buck = f*r*expr*buck1i[typej]-
                g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq+t*buck2i[typej];
              if (eflag) evdwl = f*expr*buckai[typej] -
                           g6*((a2+1.0)*a2+0.5)*x2+t*buckci[typej];
            }
          }
          else {                                              //table real space
          } else {                                              //table real space
            union_int_float_t disp_t;
            disp_t.f = rsq;
            const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
@@ -560,21 +562,18 @@ void PairBuckLongCoulLong::compute(int eflag, int vflag)
            if (ni == 0) {
              force_buck = r*expr*buck1i[typej]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*buckci[typej];
              if (eflag) evdwl = expr*buckai[typej]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*buckci[typej];
            }
            else {                                             //speial case
            } else {                                             //special case
              double f = special_lj[ni], t = rn*(1.0-f);
              force_buck = f*r*expr*buck1i[typej] -(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*buckci[typej] +t*buck2i[typej];
              if (eflag) evdwl = f*expr*buckai[typej] -(edisptable[disp_k]+f_disp*dedisptable[disp_k])*buckci[typej]+t*buckci[typej];
            }
          }
        }
        else {                                                // cut
        } else {                                                // cut
          if (ni == 0) {
            force_buck = r*expr*buck1i[typej]-rn*buck2i[typej];
            if (eflag) evdwl = expr*buckai[typej] -
                         rn*buckci[typej]-offseti[typej];
          }
          else {                                        // special case
          } else {                                        // special case
            double f = special_lj[ni];
            force_buck = f*(r*expr*buck1i[typej]-rn*buck2i[typej]);
            if (eflag)
@@ -1018,7 +1017,7 @@ double PairBuckLongCoulLong::single(int i, int j, int itype, int jtype,
        g6*((a2+1.0)*a2+0.5)*x2+t*buck_c[itype][jtype];
    } else {                                                // cut
      force_buck =
        buck1[itype][jtype]*r*expr-factor_buck*buck_c[itype][jtype]*r6inv;
        factor_buck*(buck1[itype][jtype]*r*expr-buck2[itype][jtype]*r6inv);
      eng += buck_a[itype][jtype]*expr-
        factor_buck*(buck_c[itype][jtype]*r6inv-offset[itype][jtype]);
    }
+5 −0
Original line number Diff line number Diff line
@@ -92,6 +92,8 @@ void PairLJLongCoulLong::settings(int narg, char **arg)
    error->warning(FLERR,"Using largest cutoff for lj/long/coul/long");
  if (!*(++arg))
    error->all(FLERR,"Cutoffs missing in pair_style lj/long/coul/long");
  if (!((ewald_order^ewald_off) & (1<<6)))
    dispersionflag = 0;
  if (!((ewald_order^ewald_off) & (1<<1)))
    error->all(FLERR,
               "Coulomb cut not supported in pair_style lj/long/coul/long");
@@ -386,6 +388,7 @@ void PairLJLongCoulLong::write_restart_settings(FILE *fp)
  fwrite(&ncoultablebits,sizeof(int),1,fp);
  fwrite(&tabinner,sizeof(double),1,fp);
  fwrite(&ewald_order,sizeof(int),1,fp);
  fwrite(&dispersionflag,sizeof(int),1,fp);
}

/* ----------------------------------------------------------------------
@@ -402,6 +405,7 @@ void PairLJLongCoulLong::read_restart_settings(FILE *fp)
    utils::sfread(FLERR,&ncoultablebits,sizeof(int),1,fp,NULL,error);
    utils::sfread(FLERR,&tabinner,sizeof(double),1,fp,NULL,error);
    utils::sfread(FLERR,&ewald_order,sizeof(int),1,fp,NULL,error);
    utils::sfread(FLERR,&dispersionflag,sizeof(int),1,fp,NULL,error);
  }
  MPI_Bcast(&cut_lj_global,1,MPI_DOUBLE,0,world);
  MPI_Bcast(&cut_coul,1,MPI_DOUBLE,0,world);
@@ -410,6 +414,7 @@ void PairLJLongCoulLong::read_restart_settings(FILE *fp)
  MPI_Bcast(&ncoultablebits,1,MPI_INT,0,world);
  MPI_Bcast(&tabinner,1,MPI_DOUBLE,0,world);
  MPI_Bcast(&ewald_order,1,MPI_INT,0,world);
  MPI_Bcast(&dispersionflag,1,MPI_INT,0,world);
}

/* ----------------------------------------------------------------------
+5 −0
Original line number Diff line number Diff line
@@ -1440,6 +1440,8 @@ void PairLJLongTIP4PLong::settings(int narg, char **arg)
  if (!comm->me && ewald_order==((1<<1)|(1<<6)))
    error->warning(FLERR,
                   "Using largest cutoff for pair_style lj/long/tip4p/long");
  if (!((ewald_order^ewald_off) & (1<<6)))
    dispersionflag = 0;
  if (!((ewald_order^ewald_off)&(1<<1)))
    error->all(FLERR,
               "Coulomb cut not supported in pair_style lj/long/tip4p/long");
@@ -1532,6 +1534,7 @@ void PairLJLongTIP4PLong::write_restart_settings(FILE *fp)
  fwrite(&ncoultablebits,sizeof(int),1,fp);
  fwrite(&tabinner,sizeof(double),1,fp);
  fwrite(&ewald_order,sizeof(int),1,fp);
  fwrite(&dispersionflag,sizeof(int),1,fp);
}

/* ----------------------------------------------------------------------
@@ -1554,6 +1557,7 @@ void PairLJLongTIP4PLong::read_restart_settings(FILE *fp)
    utils::sfread(FLERR,&ncoultablebits,sizeof(int),1,fp,NULL,error);
    utils::sfread(FLERR,&tabinner,sizeof(double),1,fp,NULL,error);
    utils::sfread(FLERR,&ewald_order,sizeof(int),1,fp,NULL,error);
    utils::sfread(FLERR,&dispersionflag,sizeof(int),1,fp,NULL,error);
  }

  MPI_Bcast(&typeO,1,MPI_INT,0,world);
@@ -1569,6 +1573,7 @@ void PairLJLongTIP4PLong::read_restart_settings(FILE *fp)
  MPI_Bcast(&ncoultablebits,1,MPI_INT,0,world);
  MPI_Bcast(&tabinner,1,MPI_DOUBLE,0,world);
  MPI_Bcast(&ewald_order,1,MPI_INT,0,world);
  MPI_Bcast(&dispersionflag,1,MPI_INT,0,world);
}

/* ----------------------------------------------------------------------