From dba5a5ad3fb50dfda6b0f77dd9e27c622e8dc85b Mon Sep 17 00:00:00 2001 From: athomps Date: Wed, 12 Dec 2012 22:30:41 +0000 Subject: [PATCH] Small changes to fix adapt for solid free energies git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@9150 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/compute_ti.cpp | 91 ++++++++++++++++++++++++++++++++++++++++------ src/fix_adapt.cpp | 20 +++++++++- src/fix_adapt.h | 5 +++ src/pair_gauss.cpp | 6 ++- 4 files changed, 108 insertions(+), 14 deletions(-) diff --git a/src/compute_ti.cpp b/src/compute_ti.cpp index 749a169c45..c2177f3dd9 100644 --- a/src/compute_ti.cpp +++ b/src/compute_ti.cpp @@ -16,6 +16,7 @@ ------------------------------------------------------------------------- */ #include "mpi.h" +#include "atom.h" #include "string.h" #include "compute_ti.h" #include "update.h" @@ -40,18 +41,23 @@ ComputeTI::ComputeTI(LAMMPS *lmp, int narg, char **arg) : if (narg < 4) error->all(FLERR,"Illegal compute ti command"); peflag = 1; + peratom_flag = 1; + peatomflag = 1; scalar_flag = 1; extscalar = 1; timeflag = 1; // terms come in triplets + // changed to quadruplets to include atom type - nterms = (narg-3) / 3; - if (narg != 3*nterms + 3) error->all(FLERR,"Illegal compute ti command"); + nterms = (narg-3) / 4; + if (narg != 4*nterms + 3) error->all(FLERR,"Illegal compute ti command"); which = new int[nterms]; ivar1 = new int[nterms]; ivar2 = new int[nterms]; + ilo = new int[nterms]; + ihi = new int[nterms]; var1 = new char*[nterms]; var2 = new char*[nterms]; pptr = new Pair*[nterms]; @@ -65,15 +71,17 @@ ComputeTI::ComputeTI(LAMMPS *lmp, int narg, char **arg) : int iarg = 3; while (iarg < narg) { - if (iarg+3 > narg) error->all(FLERR,"Illegal compute ti command"); + if (iarg+4 > narg) error->all(FLERR,"Illegal compute ti command"); if (strcmp(arg[iarg],"kspace") == 0) which[nterms] = KSPACE; else if (strcmp(arg[iarg],"tail") == 0) which[nterms] = TAIL; else { - which[nterms] = PAIR; - int n = strlen(arg[iarg]) + 1; - pstyle[nterms] = new char[n]; - strcpy(pstyle[nterms],arg[iarg]); - } + which[nterms] = PAIR;} + int n = strlen(arg[iarg]) + 1; + pstyle[nterms] = new char[n]; + strcpy(pstyle[nterms],arg[iarg]); + force->bounds(arg[iarg+1],atom->ntypes,ilo[nterms],ihi[nterms]); + + iarg += 1; if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) { int n = strlen(&arg[iarg+1][2]) + 1; @@ -105,6 +113,8 @@ ComputeTI::~ComputeTI() delete [] ivar2; delete [] var1; delete [] var2; + delete [] ilo; + delete [] ihi; delete [] pptr; delete [] pstyle; } @@ -153,22 +163,79 @@ double ComputeTI::compute_scalar() double dUdl = 0.0; for (int m = 0; m < nterms; m++) { + int total_flag = 0; + if ((ihi[m]-ilo[m])==atom->ntypes) total_flag = 1; + eng = 0.0; value1 = input->variable->compute_equal(ivar1[m]); value2 = input->variable->compute_equal(ivar2[m]); if (value1 == 0.0) continue; if (which[m] == PAIR) { - eng = pptr[m]->eng_vdwl + pptr[m]->eng_coul; - MPI_Allreduce(&eng,&engall,1,MPI_DOUBLE,MPI_SUM,world); + int ntypes = atom->ntypes; + int *mask = atom->mask; + if (total_flag) { + eng = pptr[m]->eng_vdwl + pptr[m]->eng_coul; + MPI_Allreduce(&eng,&engall,1,MPI_DOUBLE,MPI_SUM,world); + } + else { + int nlocal = atom->nlocal; + int npair = nlocal; + int *mask = atom->mask; + int *type = atom->type; + + double *eatom = pptr[m]->eatom; + + if (force->newton) npair += atom->nghost; + for(int i = 0; i < npair; i++) + { + if ((ilo[m]<=type[i])&(ihi[m]>=type[i])) + eng += eatom[i]; + } + MPI_Allreduce(&eng,&engall,1,MPI_DOUBLE,MPI_SUM,world); + } dUdl += engall/value1 * value2; } else if (which[m] == TAIL) { double vol = domain->xprd*domain->yprd*domain->zprd; - eng = force->pair->etail / vol; + if (total_flag) + eng = force->pair->etail / vol; + else { + eng = 0; + for (int it = 1; it <= atom->ntypes; it++) { + int jt; + if ((it >= ilo[m])&&(it <=ihi[m])) jt = it; + else jt = ilo[m]; + for (; jt <=ihi[m];jt++) { + if ((force->pair->tail_flag)&&(force->pair->setflag[it][jt])) { + double cut = force->pair->init_one(it,jt); + eng += force->pair->etail_ij; + } + if (it !=jt) eng += force->pair->etail_ij; + } + } + eng /= vol; + } dUdl += eng/value1 * value2; } else if (which[m] == KSPACE) { - eng = force->kspace->energy; + + int ntypes = atom->ntypes; + int *mask = atom->mask; + if (total_flag) + eng = force->kspace->energy; + else { + int nlocal = atom->nlocal; + int npair = nlocal; + int *mask = atom->mask; + int *type = atom->type; + double *eatom = force->kspace->eatom; + eng = 0; + for(int i = 0; i < nlocal; i++) + if ((ilo[m]<=type[i])&(ihi[m]>=type[i])) + eng += eatom[i]; + MPI_Allreduce(&eng,&engall,1,MPI_DOUBLE,MPI_SUM,world); + eng = engall; + } dUdl += eng/value1 * value2; } } diff --git a/src/fix_adapt.cpp b/src/fix_adapt.cpp index 22bed86e6a..2f0c735efc 100644 --- a/src/fix_adapt.cpp +++ b/src/fix_adapt.cpp @@ -33,7 +33,7 @@ using namespace FixConst; using namespace MathConst; enum{PAIR,KSPACE,ATOM}; -enum{DIAMETER}; +enum{DIAMETER,CHARGE}; /* ---------------------------------------------------------------------- */ @@ -71,6 +71,7 @@ FixAdapt::FixAdapt(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) nadapt = 0; diamflag = 0; + chgflag = 0; iarg = 4; while (iarg < narg) { @@ -110,6 +111,9 @@ FixAdapt::FixAdapt(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) if (strcmp(arg[iarg+1],"diameter") == 0) { adapt[nadapt].aparam = DIAMETER; diamflag = 1; + } else if (strcmp(arg[iarg+1],"charge") == 0) { + adapt[nadapt].aparam = CHARGE; + chgflag = 1; } else error->all(FLERR,"Illegal fix adapt command"); if (strstr(arg[iarg+2],"v_") == arg[iarg+2]) { int n = strlen(&arg[iarg+2][2]) + 1; @@ -229,6 +233,10 @@ void FixAdapt::init() if (!atom->radius_flag) error->all(FLERR,"Fix adapt requires atom attribute diameter"); } + if (ad->aparam == CHARGE) { + if (!atom->q_flag) + error->all(FLERR,"Fix adapt requires atom attribute charge"); + } } } @@ -334,6 +342,16 @@ void FixAdapt::change_settings() radius[i]*radius[i]*radius[i] * density; } } + } else if (ad->aparam == CHARGE) { + double *q = atom->q; + int *mask = atom->mask; + int nlocal = atom->nlocal; + for (i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { +// cout<<"Fix adapt called now!\n"; +// cout<<" Charge of atom: "<all(FLERR,"All pair coeffs are not set"); + + // 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"); 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;