diff --git a/src/MANYBODY/fix_qeq_comb.cpp b/src/MANYBODY/fix_qeq_comb.cpp index 26d91e83aa..e8a603ce3d 100644 --- a/src/MANYBODY/fix_qeq_comb.cpp +++ b/src/MANYBODY/fix_qeq_comb.cpp @@ -5,7 +5,7 @@ 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 + 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. @@ -20,16 +20,20 @@ #include "math.h" #include "stdlib.h" #include "string.h" +#include "pair_comb.h" #include "fix_qeq_comb.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "neigh_request.h" #include "atom.h" +#include "comm.h" #include "force.h" #include "group.h" #include "respa.h" -#include "pair_comb.h" #include "update.h" #include "memory.h" #include "error.h" - +#include "iostream" using namespace LAMMPS_NS; using namespace FixConst; @@ -60,17 +64,17 @@ FixQEQComb::FixQEQComb(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) if (strcmp(arg[iarg],"file") == 0) { if (iarg+2 > narg) error->all(FLERR,"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(FLERR,str); - } + 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(FLERR,str); + } } iarg += 2; } else error->all(FLERR,"Illegal fix qeq/comb command"); } - + nmax = atom->nmax; memory->create(qf,nmax,"qeq:qf"); memory->create(q1,nmax,"qeq:q1"); @@ -82,6 +86,10 @@ FixQEQComb::FixQEQComb(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) qf[i] = 0.0; + + comb = NULL; + + comm_forward = 1; } /* ---------------------------------------------------------------------- */ @@ -112,13 +120,19 @@ void FixQEQComb::init() error->all(FLERR,"Fix qeq/comb requires atom attribute q"); comb = (PairComb *) force->pair_match("comb",1); - if (comb == NULL) error->all(FLERR,"Must use pair_style comb with fix qeq/comb"); + if (comb == NULL) + error->all(FLERR,"Must use pair_style comb with fix qeq/comb"); if (strstr(update->integrate_style,"respa")) nlevels_respa = ((Respa *) update->integrate)->nlevels; ngroup = group->count(igroup); if (ngroup == 0) error->all(FLERR,"Fix qeq/comb group has no atoms"); + + int irequest = neighbor->request(this); + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->full = 1; + neighbor->requests[irequest]->ghost = 1; } /* ---------------------------------------------------------------------- */ @@ -140,7 +154,7 @@ void FixQEQComb::setup(int vflag) void FixQEQComb::post_force(int vflag) { - int i,iloop,loopmax; + int i,ii,iloop,loopmax,inum,*ilist; double heatpq,qmass,dtq,dtq2; double enegchkall,enegmaxall; @@ -161,77 +175,90 @@ void FixQEQComb::post_force(int vflag) memory->create(q2,nmax,"qeq:q2"); vector_atom = qf; } - + // more loops for first-time charge equilibrium - iloop = 0; - if (firstflag) loopmax = 5000; - else loopmax = 2000; + iloop = 0; + if (firstflag) loopmax = 500; + else loopmax = 200; // charge-equilibration loop if (me == 0 && fp) fprintf(fp,"Charge equilibration on step " BIGINT_FORMAT "\n", - update->ntimestep); - + update->ntimestep); + heatpq = 0.05; - qmass = 0.000548580; - dtq = 0.0006; + qmass = 0.016; + dtq = 0.01; dtq2 = 0.5*dtq*dtq/qmass; double enegchk = 0.0; - double enegtot = 0.0; + double enegtot = 0.0; double enegmax = 0.0; double *q = atom->q; int *mask = atom->mask; + int *tag = atom->tag; int nlocal = atom->nlocal; - for (i = 0; i < nlocal; i++) + inum = comb->list->inum; + ilist = comb->list->ilist; + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; q1[i] = q2[i] = qf[i] = 0.0; + } for (iloop = 0; iloop < loopmax; iloop ++ ) { - for (i = 0; i < nlocal; i++) + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; if (mask[i] & groupbit) { - q1[i] += qf[i]*dtq2 - heatpq*q1[i]; - q[i] += q1[i]; + q1[i] += qf[i]*dtq2 - heatpq*q1[i]; + q[i] += q1[i]; } + } + comm->forward_comm_fix(this); - enegtot = comb->yasu_char(qf,igroup); + if(comb) enegtot = comb->yasu_char(qf,igroup); enegtot /= ngroup; enegchk = enegmax = 0.0; - - for (i = 0; i < nlocal ; i++) + + for (ii = 0; ii < inum ; ii++) { + i = ilist[ii]; if (mask[i] & groupbit) { - q2[i] = enegtot-qf[i]; - enegmax = MAX(enegmax,fabs(q2[i])); - enegchk += fabs(q2[i]); - qf[i] = q2[i]; + 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 <= 100.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++) + "enegmax %.6g, fq deviation: %.6g\n", + iloop,enegtot,enegmax,enegchk); + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; if (mask[i] & groupbit) - q1[i] += qf[i]*dtq2 - heatpq*q1[i]; - } - + 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); + iloop,enegchk); } } @@ -251,3 +278,30 @@ double FixQEQComb::memory_usage() double bytes = atom->nmax*3 * sizeof(double); return bytes; } +/* ---------------------------------------------------------------------- */ + +int FixQEQComb::pack_comm(int n, int *list, double *buf, int pbc_flag, int *pbc) +{ + int i,j,m; + + m = 0; + for (i = 0; i < n; i ++) { + j = list[i]; + buf[m++] = atom->q[j]; + } + return 1; +} + +/* ---------------------------------------------------------------------- */ + +void FixQEQComb::unpack_comm(int n, int first, double *buf) +{ + int i,m,last; + + m = 0; + last = first + n ; + for (i = first; i < last; i++) atom->q[i] = buf[m++]; +} + +/* ---------------------------------------------------------------------- */ + diff --git a/src/MANYBODY/fix_qeq_comb.h b/src/MANYBODY/fix_qeq_comb.h index 7452fcfa0e..62d040e64f 100644 --- a/src/MANYBODY/fix_qeq_comb.h +++ b/src/MANYBODY/fix_qeq_comb.h @@ -35,6 +35,8 @@ class FixQEQComb : public Fix { virtual void post_force(int); void post_force_respa(int,int,int); double memory_usage(); + int pack_comm(int , int *, double *, int, int *); + void unpack_comm(int , int , double *); protected: int me,firstflag;