ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

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

Browse Source
This commit is contained in:
sjplimp
2012-06-20 14:59:19 +00:00
parent c1e0bec7f6
commit bf54b6b9ff
2 changed files with 97 additions and 41 deletions
Download Patch File Download Diff File Expand all files Collapse all files

136
src/MANYBODY/fix_qeq_comb.cpp
View File

@ -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++];
}
/* ---------------------------------------------------------------------- */

2
src/MANYBODY/fix_qeq_comb.h
View File

@ -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;