diff --git a/src/MANYBODY/Install.sh b/src/MANYBODY/Install.sh
index 2946e0e3e..efa85e3aa 100644
--- a/src/MANYBODY/Install.sh
+++ b/src/MANYBODY/Install.sh
@@ -1,58 +1,58 @@
# Install/unInstall package files in LAMMPS
# for unInstall, also unInstall/Install OPT package if installed
# so it will remove OPT files that depend on MANYBODY files,
# then replace others
if (test $1 = 1) then
- cp fix_qeq.cpp ..
+ cp fix_qeq_comb.cpp ..
cp pair_airebo.cpp ..
cp pair_comb.cpp ..
cp pair_eam.cpp ..
cp pair_eam_alloy.cpp ..
cp pair_eam_fs.cpp ..
cp pair_sw.cpp ..
cp pair_tersoff.cpp ..
cp pair_tersoff_zbl.cpp ..
- cp fix_qeq.h ..
+ cp fix_qeq_comb.h ..
cp pair_airebo.h ..
cp pair_comb.h ..
cp pair_eam.h ..
cp pair_eam_alloy.h ..
cp pair_eam_fs.h ..
cp pair_sw.h ..
cp pair_tersoff.h ..
cp pair_tersoff_zbl.h ..
if (test -e ../pair_lj_cut_opt.h) then
cd ../OPT; sh Install.sh 1
fi
elif (test $1 = 0) then
- rm ../fix_qeq.cpp
+ rm ../fix_qeq_comb.cpp
rm ../pair_airebo.cpp
rm ../pair_comb.cpp
rm ../pair_eam.cpp
rm ../pair_eam_alloy.cpp
rm ../pair_eam_fs.cpp
rm ../pair_sw.cpp
rm ../pair_tersoff.cpp
rm ../pair_tersoff_zbl.cpp
- rm ../pair_qeq.h
+ rm ../pair_qeq_comb.h
rm ../pair_airebo.h
rm ../pair_comb.h
rm ../pair_eam.h
rm ../pair_eam_alloy.h
rm ../pair_eam_fs.h
rm ../pair_sw.h
rm ../pair_tersoff.h
rm ../pair_tersoff_zbl.h
if (test -e ../pair_eam_opt.h) then
cd ../OPT; sh Install.sh 0; sh Install.sh 1
fi
fi
diff --git a/src/MANYBODY/fix_qeq_comb.cpp b/src/MANYBODY/fix_qeq_comb.cpp
new file mode 100644
index 000000000..b7c87e340
--- /dev/null
+++ b/src/MANYBODY/fix_qeq_comb.cpp
@@ -0,0 +1,253 @@
+/* ----------------------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ Copyright (2003) Sandia Corporation. Under the terms of Contract
+ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+ certain rights in this software. This software is distributed under
+ the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+ Contributing author: Tzu-Ray Shan (U Florida, rayshan@ufl.edu)
+------------------------------------------------------------------------- */
+
+#include "mpi.h"
+#include "math.h"
+#include "stdlib.h"
+#include "string.h"
+#include "fix_qeq_comb.h"
+#include "atom.h"
+#include "force.h"
+#include "group.h"
+#include "respa.h"
+#include "pair_comb.h"
+#include "update.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+#define MIN(A,B) ((A) < (B)) ? (A) : (B)
+#define MAX(A,B) ((A) > (B)) ? (A) : (B)
+
+/* ---------------------------------------------------------------------- */
+
+FixQEQComb::FixQEQComb(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
+{
+ if (narg < 5) error->all("Illegal fix qeq/comb command");
+
+ peratom_flag = 1;
+ size_peratom_cols = 0;
+ peratom_freq = 1;
+
+ nevery = force->inumeric(arg[3]);
+ precision = force->numeric(arg[4]);
+
+ if (nevery <= 0 || precision <= 0.0)
+ error->all("Illegal fix qeq/comb command");
+
+ MPI_Comm_rank(world,&me);
+
+ // optional args
+
+ fp = NULL;
+
+ int iarg = 5;
+ while (iarg < narg) {
+ if (strcmp(arg[iarg],"file") == 0) {
+ if (iarg+2 > narg) error->all("Illegal fix qeq/comb command");
+ if (me == 0) {
+ fp = fopen(arg[iarg+1],"w");
+ if (fp == NULL) {
+ char str[128];
+ sprintf(str,"Cannot open fix qeq/comb file %s",arg[iarg+1]);
+ error->one(str);
+ }
+ }
+ iarg += 2;
+ } else error->all("Illegal fix qeq/comb command");
+ }
+
+ nmax = atom->nmax;
+ qf = (double *) memory->smalloc(nmax*sizeof(double),"qeq:qf");
+ q1 = (double *) memory->smalloc(nmax*sizeof(double),"qeq:q1");
+ q2 = (double *) memory->smalloc(nmax*sizeof(double),"qeq:q2");
+ vector_atom = qf;
+
+ // zero the vector since dump may access it on timestep 0
+ // zero the vector since a variable may access it before first run
+
+ int nlocal = atom->nlocal;
+ for (int i = 0; i < nlocal; i++) qf[i] = 0.0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+FixQEQComb::~FixQEQComb()
+{
+ if (me == 0 && fp) fclose(fp);
+ memory->sfree(qf);
+ memory->sfree(q1);
+ memory->sfree(q2);
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixQEQComb::setmask()
+{
+ int mask = 0;
+ mask |= POST_FORCE;
+ mask |= POST_FORCE_RESPA;
+ return mask;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixQEQComb::init()
+{
+ if (!atom->q_flag)
+ error->all("Atoms must have charge to use fix qeq/comb");
+
+ comb = (PairComb *) force->pair_match("comb",1);
+ if (comb == NULL) error->all("Must use pair_style comb with fix qeq/comb");
+
+ if (strcmp(update->integrate_style,"respa") == 0)
+ nlevels_respa = ((Respa *) update->integrate)->nlevels;
+
+ ngroup = group->count(igroup);
+ if (ngroup == 0.0) error->all("Fix qeq/comb group has no atoms");
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixQEQComb::setup(int vflag)
+{
+ firstflag = 1;
+ if (strcmp(update->integrate_style,"verlet") == 0)
+ post_force(vflag);
+ else {
+ ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
+ post_force_respa(vflag,nlevels_respa-1,0);
+ ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
+ }
+ firstflag = 0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixQEQComb::post_force(int vflag)
+{
+ int i,iloop,loopmax;
+ double heatpq,qmass,dtq,dtq2;
+ double enegchkall,enegmaxall;
+
+ if (update->ntimestep % nevery) return;
+
+ // reallocate work arrays if necessary
+ // qf = charge force
+ // q1 = charge displacement
+ // q2 = tmp storage of charge force for next iteration
+
+ if (atom->nmax > nmax) {
+ memory->sfree(qf);
+ memory->sfree(q1);
+ memory->sfree(q2);
+ nmax = atom->nmax;
+ qf = (double *) memory->smalloc(nmax*sizeof(double),"qeq:qf");
+ q1 = (double *) memory->smalloc(nmax*sizeof(double),"qeq:q1");
+ q2 = (double *) memory->smalloc(nmax*sizeof(double),"qeq:q2");
+ vector_atom = qf;
+ }
+
+ // more loops for first-time charge equilibrium
+
+ iloop = 0;
+ if (firstflag) loopmax = 1000;
+ else loopmax = 500;
+
+ // charge-equilibration loop
+
+ if (me == 0 && fp)
+ fprintf(fp,"Charge equilibration on step %d\n",update->ntimestep);
+
+ heatpq = 0.01;
+ qmass = 0.06;
+ dtq = 0.040;
+ dtq2 = 0.5*dtq*dtq/qmass;
+
+ double enegchk = 0.0;
+ double enegtot = 0.0;
+ double enegmax = 0.0;
+
+ double *q = atom->q;
+ int *mask = atom->mask;
+ int nlocal = atom->nlocal;
+
+ for (i = 0; i < nlocal; i++)
+ q1[i] = q2[i] = qf[i] = 0.0;
+
+ for (iloop = 0; iloop < loopmax; iloop ++ ) {
+ for (i = 0; i < nlocal; i++)
+ if (mask[i] & groupbit) {
+ q1[i] = qf[i]*dtq2 - heatpq*q1[i];
+ q[i] += q1[i];
+ }
+
+ enegtot = comb->yasu_char(qf,igroup);
+ enegtot /= ngroup;
+ enegchk = enegmax = 0.0;
+
+ for (i = 0; i < nlocal ; i++)
+ if (mask[i] & groupbit) {
+ q2[i] = enegtot-qf[i];
+ enegmax = MAX(enegmax,fabs(q2[i]));
+ enegchk += fabs(q2[i]);
+ qf[i] = q2[i];
+ }
+
+ MPI_Allreduce(&enegchk,&enegchkall,1,MPI_DOUBLE,MPI_SUM,world);
+ enegchk = enegchkall/ngroup;
+ MPI_Allreduce(&enegmax,&enegmaxall,1,MPI_DOUBLE,MPI_MAX,world);
+ enegmax = enegmaxall;
+
+ if (enegchk <= precision && enegmax <= 10.0*precision) break;
+
+ if (me == 0 && fp)
+ fprintf(fp," iteration: %d, enegtot %.6g, "
+ "enegmax %.6g, fq deviation: %.6g\n",
+ iloop,enegtot,enegmax,enegchk);
+
+ for (i = 0; i < nlocal; i++)
+ if (mask[i] & groupbit)
+ q1[i] += qf[i]*dtq2 - heatpq*q1[i];
+ }
+
+ if (me == 0 && fp) {
+ if (iloop == loopmax)
+ fprintf(fp,"Charges did not converge in %d iterations\n",iloop);
+ else
+ fprintf(fp,"Charges converged in %d iterations to %.10f tolerance\n",
+ iloop,enegchk);
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixQEQComb::post_force_respa(int vflag, int ilevel, int iloop)
+{
+ if (ilevel == nlevels_respa-1) post_force(vflag);
+}
+
+/* ----------------------------------------------------------------------
+ memory usage of local atom-based arrays
+------------------------------------------------------------------------- */
+
+double FixQEQComb::memory_usage()
+{
+ double bytes = atom->nmax*3 * sizeof(double);
+ return bytes;
+}
diff --git a/src/MANYBODY/fix_qeq_comb.h b/src/MANYBODY/fix_qeq_comb.h
new file mode 100644
index 000000000..76c694ef4
--- /dev/null
+++ b/src/MANYBODY/fix_qeq_comb.h
@@ -0,0 +1,54 @@
+/* ----------------------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ Copyright (2003) Sandia Corporation. Under the terms of Contract
+ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+ certain rights in this software. This software is distributed under
+ the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef FIX_CLASS
+
+FixStyle(qeq/comb,FixQEQComb)
+
+#else
+
+#ifndef LMP_FIX_QEQ_COMB_H
+#define LMP_FIX_QEQ_COMB_H
+
+#include "stdio.h"
+#include "fix.h"
+
+namespace LAMMPS_NS {
+
+class FixQEQComb : public Fix {
+ public:
+ FixQEQComb(class LAMMPS *, int, char **);
+ ~FixQEQComb();
+ int setmask();
+ void init();
+ void setup(int);
+ void post_force(int);
+ void post_force_respa(int,int,int);
+ double memory_usage();
+
+ private:
+ int me,firstflag;
+ double precision;
+ int nlevels_respa;
+ double ngroup;
+ FILE *fp;
+
+ class PairComb *comb;
+ int nmax;
+ double *qf,*q1,*q2;
+};
+
+}
+
+#endif
+#endif