From 166a9c5bc791c523bdd3dbc351ab128e7617f741 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Tue, 7 Oct 2014 23:03:23 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@12621 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/KSPACE/ewald.cpp | 2 +- src/KSPACE/msm.cpp | 2 +- src/KSPACE/pppm.cpp | 2 +- src/KSPACE/pppm_cg.cpp | 2 +- src/KSPACE/pppm_disp.cpp | 2 +- src/KSPACE/pppm_stagger.cpp | 2 +- src/QEQ/fix_qeq.cpp | 1 + src/QEQ/fix_qeq_dynamic.cpp | 2 + src/QEQ/fix_qeq_point.cpp | 3 + src/QEQ/fix_qeq_shielded.cpp | 3 + src/QEQ/fix_qeq_slater.cpp | 4 +- src/USER-FEP/fix_adapt_fep.cpp | 202 +++++++++++++++++++++++++-------- src/USER-FEP/fix_adapt_fep.h | 3 + src/kspace.cpp | 9 +- src/kspace.h | 2 + 15 files changed, 184 insertions(+), 57 deletions(-) diff --git a/src/KSPACE/ewald.cpp b/src/KSPACE/ewald.cpp index 1acf04e220..054bf343e4 100644 --- a/src/KSPACE/ewald.cpp +++ b/src/KSPACE/ewald.cpp @@ -439,7 +439,7 @@ void Ewald::compute(int eflag, int vflag) energy += ug[k] * (sfacrl_all[k]*sfacrl_all[k] + sfacim_all[k]*sfacim_all[k]); - if (atom->natoms != natoms_original) { + if (qsum_update_flag || (atom->natoms != natoms_original)) { qsum_qsq(0); natoms_original = atom->natoms; } diff --git a/src/KSPACE/msm.cpp b/src/KSPACE/msm.cpp index 9ba0c58191..7a33c2a765 100644 --- a/src/KSPACE/msm.cpp +++ b/src/KSPACE/msm.cpp @@ -576,7 +576,7 @@ void MSM::compute(int eflag, int vflag) MPI_Allreduce(&energy,&energy_all,1,MPI_DOUBLE,MPI_SUM,world); energy = energy_all; - if (atom->natoms != natoms_original) { + if (qsum_update_flag || (atom->natoms != natoms_original)) { qsum_qsq(0); natoms_original = atom->natoms; } diff --git a/src/KSPACE/pppm.cpp b/src/KSPACE/pppm.cpp index a2672ca7b7..f0ec0b9db5 100644 --- a/src/KSPACE/pppm.cpp +++ b/src/KSPACE/pppm.cpp @@ -677,7 +677,7 @@ void PPPM::compute(int eflag, int vflag) MPI_Allreduce(&energy,&energy_all,1,MPI_DOUBLE,MPI_SUM,world); energy = energy_all; - if (atom->natoms != natoms_original) { + if (qsum_update_flag || (atom->natoms != natoms_original)) { qsum_qsq(0); natoms_original = atom->natoms; } diff --git a/src/KSPACE/pppm_cg.cpp b/src/KSPACE/pppm_cg.cpp index be84a5fc0b..3ebcb5d3a4 100644 --- a/src/KSPACE/pppm_cg.cpp +++ b/src/KSPACE/pppm_cg.cpp @@ -218,7 +218,7 @@ void PPPMCG::compute(int eflag, int vflag) MPI_Allreduce(&energy,&energy_all,1,MPI_DOUBLE,MPI_SUM,world); energy = energy_all; - if (atom->natoms != natoms_original) { + if (qsum_update_flag || (atom->natoms != natoms_original)) { qsum_qsq(0); natoms_original = atom->natoms; } diff --git a/src/KSPACE/pppm_disp.cpp b/src/KSPACE/pppm_disp.cpp index 37fa0b46f6..78611360bd 100755 --- a/src/KSPACE/pppm_disp.cpp +++ b/src/KSPACE/pppm_disp.cpp @@ -1135,7 +1135,7 @@ void PPPMDisp::compute(int eflag, int vflag) energy_1 *= 0.5*volume; energy_6 *= 0.5*volume; - if (atom->natoms != natoms_original) { + if (qsum_update_flag || (atom->natoms != natoms_original)) { qsum_qsq(0); natoms_original = atom->natoms; } diff --git a/src/KSPACE/pppm_stagger.cpp b/src/KSPACE/pppm_stagger.cpp index f0ee7e10dc..bc1a1ffe24 100755 --- a/src/KSPACE/pppm_stagger.cpp +++ b/src/KSPACE/pppm_stagger.cpp @@ -212,7 +212,7 @@ void PPPMStagger::compute(int eflag, int vflag) MPI_Allreduce(&energy,&energy_all,1,MPI_DOUBLE,MPI_SUM,world); energy = energy_all; - if (atom->natoms != natoms_original) { + if (qsum_update_flag || (atom->natoms != natoms_original)) { qsum_qsq(0); natoms_original = atom->natoms; } diff --git a/src/QEQ/fix_qeq.cpp b/src/QEQ/fix_qeq.cpp index 9e97182d99..c2aab73c35 100644 --- a/src/QEQ/fix_qeq.cpp +++ b/src/QEQ/fix_qeq.cpp @@ -262,6 +262,7 @@ void FixQEq::reallocate_matrix() void FixQEq::init_list(int id, NeighList *ptr) { list = ptr; + if (force->kspace) force->kspace->qsum_update_flag = 1; } /* ---------------------------------------------------------------------- */ diff --git a/src/QEQ/fix_qeq_dynamic.cpp b/src/QEQ/fix_qeq_dynamic.cpp index f483d3d32f..848ca501c3 100644 --- a/src/QEQ/fix_qeq_dynamic.cpp +++ b/src/QEQ/fix_qeq_dynamic.cpp @@ -30,6 +30,7 @@ #include "force.h" #include "group.h" #include "pair.h" +#include "kspace.h" #include "respa.h" #include "memory.h" #include "error.h" @@ -167,6 +168,7 @@ void FixQEqDynamic::pre_force(int vflag) } } + if (force->kspace) force->kspace->setup(); } /* ---------------------------------------------------------------------- */ diff --git a/src/QEQ/fix_qeq_point.cpp b/src/QEQ/fix_qeq_point.cpp index 1e4ce1d41f..5f25332c44 100644 --- a/src/QEQ/fix_qeq_point.cpp +++ b/src/QEQ/fix_qeq_point.cpp @@ -28,6 +28,7 @@ #include "neigh_request.h" #include "update.h" #include "force.h" +#include "kspace.h" #include "group.h" #include "respa.h" #include "memory.h" @@ -82,6 +83,8 @@ void FixQEqPoint::pre_force(int vflag) matvecs += CG(b_t, t); // CG on t - parallel calculate_Q(); + if (force->kspace) force->kspace->setup(); + } /* ---------------------------------------------------------------------- */ diff --git a/src/QEQ/fix_qeq_shielded.cpp b/src/QEQ/fix_qeq_shielded.cpp index dfc110f75e..a2c8611af7 100644 --- a/src/QEQ/fix_qeq_shielded.cpp +++ b/src/QEQ/fix_qeq_shielded.cpp @@ -28,6 +28,7 @@ #include "neigh_request.h" #include "update.h" #include "force.h" +#include "kspace.h" #include "group.h" #include "respa.h" #include "memory.h" @@ -125,6 +126,8 @@ void FixQEqShielded::pre_force(int vflag) matvecs += CG(b_t, t); // CG on t - parallel calculate_Q(); + if (force->kspace) force->kspace->setup(); + } /* ---------------------------------------------------------------------- */ diff --git a/src/QEQ/fix_qeq_slater.cpp b/src/QEQ/fix_qeq_slater.cpp index 88e4b00175..d3ece69836 100644 --- a/src/QEQ/fix_qeq_slater.cpp +++ b/src/QEQ/fix_qeq_slater.cpp @@ -28,9 +28,9 @@ #include "neigh_request.h" #include "update.h" #include "force.h" -#include "kspace.h" #include "group.h" #include "pair.h" +#include "kspace.h" #include "respa.h" #include "math_const.h" #include "memory.h" @@ -101,6 +101,8 @@ void FixQEqSlater::pre_force(int vflag) matvecs += CG(b_t, t); // CG on t - parallel calculate_Q(); + if (force->kspace) force->kspace->setup(); + } /* ---------------------------------------------------------------------- */ diff --git a/src/USER-FEP/fix_adapt_fep.cpp b/src/USER-FEP/fix_adapt_fep.cpp index c2364889eb..794278efdc 100644 --- a/src/USER-FEP/fix_adapt_fep.cpp +++ b/src/USER-FEP/fix_adapt_fep.cpp @@ -20,13 +20,14 @@ #include "stdlib.h" #include "fix_adapt_fep.h" #include "atom.h" -#include "comm.h" #include "update.h" +#include "group.h" #include "modify.h" #include "force.h" #include "pair.h" #include "pair_hybrid.h" #include "kspace.h" +#include "fix_store.h" #include "input.h" #include "variable.h" #include "math_const.h" @@ -40,8 +41,6 @@ using namespace MathConst; enum{PAIR,KSPACE,ATOM}; enum{DIAMETER,CHARGE}; -#undef ADAPT_DEBUG - /* ---------------------------------------------------------------------- */ FixAdaptFEP::FixAdaptFEP(LAMMPS *lmp, int narg, char **arg) : @@ -51,6 +50,8 @@ FixAdaptFEP::FixAdaptFEP(LAMMPS *lmp, int narg, char **arg) : nevery = force->inumeric(FLERR,arg[3]); if (nevery < 0) error->all(FLERR,"Illegal fix adapt/fep command"); + dynamic_group_allow = 1; + // count # of adaptations nadapt = 0; @@ -169,6 +170,8 @@ FixAdaptFEP::FixAdaptFEP(LAMMPS *lmp, int narg, char **arg) : for (int m = 0; m < nadapt; m++) if (adapt[m].which == PAIR) memory->create(adapt[m].array_orig,n+1,n+1,"adapt:array_orig"); + + id_fix_diam = id_fix_chg = NULL; } /* ---------------------------------------------------------------------- */ @@ -184,6 +187,13 @@ FixAdaptFEP::~FixAdaptFEP() } } delete [] adapt; + + // check nfix in case all fixes have already been deleted + + if (id_fix_diam && modify->nfix) modify->delete_fix(id_fix_diam); + if (id_fix_chg && modify->nfix) modify->delete_fix(id_fix_chg); + delete [] id_fix_diam; + delete [] id_fix_chg; } /* ---------------------------------------------------------------------- */ @@ -196,12 +206,92 @@ int FixAdaptFEP::setmask() return mask; } +/* ---------------------------------------------------------------------- + if need to restore per-atom quantities, create new fix STORE styles +------------------------------------------------------------------------- */ + +void FixAdaptFEP::post_constructor() +{ + if (!resetflag) return; + if (!diamflag && !chgflag) return; + + // new id = fix-ID + FIX_STORE_ATTRIBUTE + // new fix group = group for this fix + + id_fix_diam = NULL; + id_fix_chg = NULL; + + char **newarg = new char*[5]; + newarg[1] = group->names[igroup]; + newarg[2] = (char *) "STORE"; + newarg[3] = (char *) "1"; + newarg[4] = (char *) "1"; + + if (diamflag) { + int n = strlen(id) + strlen("_FIX_STORE_DIAM") + 1; + id_fix_diam = new char[n]; + strcpy(id_fix_diam,id); + strcat(id_fix_diam,"_FIX_STORE_DIAM"); + newarg[0] = id_fix_diam; + modify->add_fix(5,newarg); + fix_diam = (FixStore *) modify->fix[modify->nfix-1]; + + if (fix_diam->restart_reset) fix_diam->restart_reset = 0; + else { + double *vec = fix_diam->vstore; + double *radius = atom->radius; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) vec[i] = radius[i]; + else vec[i] = 0.0; + } + } + } + + if (chgflag) { + int n = strlen(id) + strlen("_FIX_STORE_CHG") + 1; + id_fix_chg = new char[n]; + strcpy(id_fix_chg,id); + strcat(id_fix_chg,"_FIX_STORE_CHG"); + newarg[0] = id_fix_chg; + modify->add_fix(5,newarg); + fix_chg = (FixStore *) modify->fix[modify->nfix-1]; + + if (fix_chg->restart_reset) fix_chg->restart_reset = 0; + else { + double *vec = fix_chg->vstore; + double *q = atom->q; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) vec[i] = q[i]; + else vec[i] = 0.0; + } + } + } + + delete [] newarg; +} + /* ---------------------------------------------------------------------- */ void FixAdaptFEP::init() { int i,j; + // allow a dynamic group only if ATOM attribute not used + + if (group->dynamic[igroup]) + for (int i = 0; i < nadapt; i++) + if (adapt[i].which == ATOM) + error->all(FLERR,"Cannot use dynamic group with fix adapt/fep atom"); + + // when using kspace, we need to recompute some additional parameters in kspace->setup() + if (force->kspace) force->kspace->qsum_update_flag = 1; + // setup and error checks anypair = 0; @@ -217,9 +307,18 @@ void FixAdaptFEP::init() if (ad->which == PAIR) { anypair = 1; + Pair *pair = NULL; - Pair *pair = force->pair_match(ad->pstyle,1); - if (pair == NULL) error->all(FLERR,"Fix adapt/fep pair style does not exist"); + if (lmp->suffix_enable) { + char psuffix[128]; + strcpy(psuffix,ad->pstyle); + strcat(psuffix,"/"); + strcat(psuffix,lmp->suffix); + pair = force->pair_match(psuffix,1); + } + if (pair == NULL) pair = force->pair_match(ad->pstyle,1); + if (pair == NULL) + error->all(FLERR, "Fix adapt/fep pair style does not exist"); void *ptr = pair->extract(ad->pparam,ad->pdim); if (ptr == NULL) error->all(FLERR,"Fix adapt/fep pair style param not supported"); @@ -268,22 +367,18 @@ void FixAdaptFEP::init() } } -#ifdef ADAPT_DEBUG - if (comm->me == 0 && screen) { - for (int m = 0; m < nadapt; m++) { - Adapt *ad = &adapt[m]; - if (ad->which == PAIR && ad->pdim == 2) { - fprintf(screen, "###ADAPT original %s %s\n", ad->pstyle, ad->pparam); - fprintf(screen, "###ADAPT I J old_param\n"); - for (i = ad->ilo; i <= ad->ihi; i++) - for (j = MAX(ad->jlo,i); j <= ad->jhi; j++) - fprintf(screen, "###ADAPT %2d %2d %9.5f\n", i, j, - ad->array_orig[i][j]); - } - } + // fixes that store initial per-atom values + + if (id_fix_diam) { + int ifix = modify->find_fix(id_fix_diam); + if (ifix < 0) error->all(FLERR,"Could not find fix adapt storage fix ID"); + fix_diam = (FixStore *) modify->fix[ifix]; + } + if (id_fix_chg) { + int ifix = modify->find_fix(id_fix_chg); + if (ifix < 0) error->all(FLERR,"Could not find fix adapt storage fix ID"); + fix_chg = (FixStore *) modify->fix[ifix]; } -#endif - } /* ---------------------------------------------------------------------- */ @@ -349,17 +444,6 @@ void FixAdaptFEP::change_settings() for (i = ad->ilo; i <= ad->ihi; i++) for (j = MAX(ad->jlo,i); j <= ad->jhi; j++) ad->array[i][j] = value; - -#ifdef ADAPT_DEBUG - if (comm->me == 0 && screen) { - fprintf(screen, "###ADAPT change %s %s\n", ad->pstyle, ad->pparam); - fprintf(screen, "###ADAPT I J param\n"); - for (i = ad->ilo; i <= ad->ihi; i++) - for (j = MAX(ad->jlo,i); j <= ad->jhi; j++) - fprintf(screen, "###ADAPT %2d %2d %9.5f\n", i, j, ad->array[i][j]); - } -#endif - } // set kspace scale factor @@ -369,7 +453,7 @@ void FixAdaptFEP::change_settings() } else if (ad->which == ATOM) { - // set radius from diameter + // reset radius from diameter // also scale rmass to new value if (ad->aparam == DIAMETER) { @@ -381,15 +465,16 @@ void FixAdaptFEP::change_settings() double *radius = atom->radius; double *rmass = atom->rmass; int *mask = atom->mask; - int natom = atom->nlocal + atom->nghost; + int nlocal = atom->nlocal; + int nall = nlocal + atom->nghost; if (mflag == 0) { - for (i = 0; i < natom; i++) + for (i = 0; i < nall; i++) if (atype[i] >= ad->ilo && atype[i] <= ad->ihi) if (mask[i] & groupbit) radius[i] = 0.5*value; } else { - for (i = 0; i < natom; i++) + for (i = 0; i < nall; i++) if (atype[i] >= ad->ilo && atype[i] <= ad->ihi) if (mask[i] & groupbit) { density = rmass[i] / (4.0*MY_PI/3.0 * @@ -403,22 +488,12 @@ void FixAdaptFEP::change_settings() int *atype = atom->type; double *q = atom->q; int *mask = atom->mask; - int natom = atom->nlocal + atom->nghost; + int nlocal = atom->nlocal; + int nall = nlocal + atom->nghost; - for (i = 0; i < natom; i++) + for (i = 0; i < nall; i++) if (atype[i] >= ad->ilo && atype[i] <= ad->ihi) if (mask[i] & groupbit) q[i] = value; - -#ifdef ADAPT_DEBUG - if (comm->me == 0 && screen) { - fprintf(screen, "###ADAPT change charge\n"); - fprintf(screen, "###ADAPT atom I q\n"); - for (i = 0; i < atom->nlocal; i++) - if (atype[i] >= ad->ilo && atype[i] <= ad->ihi) - if (mask[i] & groupbit) - fprintf(screen, "###ADAPT %5d %2d %9.5f\n", i, atype[i], q[i]); - } -#endif } } } @@ -430,6 +505,10 @@ void FixAdaptFEP::change_settings() // and also offset and tail corrections if (anypair) force->pair->reinit(); + + if (force->kspace) + force->kspace->setup(); + } /* ---------------------------------------------------------------------- @@ -452,9 +531,36 @@ void FixAdaptFEP::restore_settings() *kspace_scale = 1.0; } else if (ad->which == ATOM) { + if (diamflag) { + double density; + double *vec = fix_diam->vstore; + double *radius = atom->radius; + double *rmass = atom->rmass; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + density = rmass[i] / (4.0*MY_PI/3.0 * + radius[i]*radius[i]*radius[i]); + radius[i] = vec[i]; + rmass[i] = 4.0*MY_PI/3.0 * radius[i]*radius[i]*radius[i] * density; + } + } + if (chgflag) { + double *vec = fix_chg->vstore; + double *q = atom->q; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) q[i] = vec[i]; + } } } if (anypair) force->pair->reinit(); + + if (force->kspace) force->kspace->setup(); } diff --git a/src/USER-FEP/fix_adapt_fep.h b/src/USER-FEP/fix_adapt_fep.h index 0fca92fcaa..cb159566cc 100644 --- a/src/USER-FEP/fix_adapt_fep.h +++ b/src/USER-FEP/fix_adapt_fep.h @@ -32,6 +32,7 @@ class FixAdaptFEP : public Fix { FixAdaptFEP(class LAMMPS *, int, char **); ~FixAdaptFEP(); int setmask(); + void post_constructor(); void init(); void setup_pre_force(int); void pre_force(int); @@ -40,6 +41,8 @@ class FixAdaptFEP : public Fix { private: int nadapt,resetflag,scaleflag,afterflag; int anypair; + char *id_fix_diam,*id_fix_chg; + class FixStore *fix_diam,*fix_chg; struct Adapt { int which,ivar; diff --git a/src/kspace.cpp b/src/kspace.cpp index e25a53b759..7727c24c7a 100644 --- a/src/kspace.cpp +++ b/src/kspace.cpp @@ -68,6 +68,8 @@ KSpace::KSpace(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) suffix_flag = Suffix::NONE; adjust_cutoff_flag = 1; scalar_pressure_flag = 0; + qsum_update_flag = 0; + warn_neutral = 1; accuracy_absolute = -1.0; accuracy_real_6 = -1.0; @@ -284,8 +286,11 @@ void KSpace::qsum_qsq(int flag) if (fabs(qsum) > SMALL) { char str[128]; sprintf(str,"System is not charge neutral, net charge = %g",qsum); - if (flag) error->all(FLERR,str); - else if (comm->me == 0) error->warning(FLERR,str); + if (warn_neutral && (comm->me == 0)) { + if (flag) error->all(FLERR,str); + else error->warning(FLERR,str); + } + warn_neutral = 0; } } diff --git a/src/kspace.h b/src/kspace.h index 6d8ae294f5..466a1de239 100644 --- a/src/kspace.h +++ b/src/kspace.h @@ -51,6 +51,7 @@ class KSpace : protected Pointers { // for LJ coefficients int slabflag; int scalar_pressure_flag; // 1 if using MSM fast scalar pressure + int qsum_update_flag; // 1 if setup() needs to call qsum_qsq() double slab_volfactor; @@ -168,6 +169,7 @@ class KSpace : protected Pointers { int minorder,overlap_allowed; int adjust_cutoff_flag; int suffix_flag; // suffix compatibility flag + int warn_neutral; // warn about non-neutral system if 1 bigint natoms_original; double scale,qqrd2e; double qsum,qsqsum,q2;