From 356091e1e0682da97e57468e120593aa28770eb2 Mon Sep 17 00:00:00 2001 From: Jinzhe Zeng Date: Sat, 26 Mar 2022 02:03:33 -0400 Subject: [PATCH 1/4] support kspace plugin I need such feature in the development of deepmd-kit. --- doc/src/Developer_plugins.rst | 2 +- doc/src/plugin.rst | 2 +- src/PLUGIN/plugin.cpp | 17 ++++++++++++++++- 3 files changed, 18 insertions(+), 3 deletions(-) diff --git a/doc/src/Developer_plugins.rst b/doc/src/Developer_plugins.rst index 264a77ed5b..c4378b6c85 100644 --- a/doc/src/Developer_plugins.rst +++ b/doc/src/Developer_plugins.rst @@ -68,7 +68,7 @@ Members of ``lammpsplugin_t`` * - author - String with the name and email of the author * - creator.v1 - - Pointer to factory function for pair, bond, angle, dihedral, improper or command styles + - Pointer to factory function for pair, bond, angle, dihedral, improper, command or kspace styles * - creator.v2 - Pointer to factory function for compute, fix, or region styles * - handle diff --git a/doc/src/plugin.rst b/doc/src/plugin.rst index 3ce9e41870..8f689e885a 100644 --- a/doc/src/plugin.rst +++ b/doc/src/plugin.rst @@ -17,7 +17,7 @@ Syntax *load* file = load plugin(s) from shared object in *file* *unload* style name = unload plugin *name* of style *style* - *style* = *pair* or *bond* or *angle* or *dihedral* or *improper* or *compute* or *fix* or *region* or *command* + *style* = *pair* or *bond* or *angle* or *dihedral* or *improper* or *compute* or *fix* or *region* or *command* or *kspace* *list* = print a list of currently loaded plugins *clear* = unload all currently loaded plugins diff --git a/src/PLUGIN/plugin.cpp b/src/PLUGIN/plugin.cpp index 3b28a32dbf..761b18760f 100644 --- a/src/PLUGIN/plugin.cpp +++ b/src/PLUGIN/plugin.cpp @@ -242,6 +242,15 @@ void plugin_register(lammpsplugin_t *plugin, void *ptr) } (*command_map)[plugin->name] = (Input::CommandCreator) plugin->creator.v1; + } else if (pstyle == "kspace") { + auto kspace_map = lmp->force->kspace_map; + if (kspace_map->find(plugin->name) != kspace_map->end()) { + if (lmp->comm->me == 0) + lmp->error->warning(FLERR, "Overriding built-in kspace style {} from plugin", + plugin->name); + } + (*kspace_map)[plugin->name] = (Force::KSpaceCreator) plugin->creator.v1; + } else { utils::logmesg(lmp, "Loading plugins for {} styles not yet implemented\n", pstyle); pluginlist.pop_back(); @@ -265,7 +274,7 @@ void plugin_unload(const char *style, const char *name, LAMMPS *lmp) (strcmp(style, "angle") != 0) && (strcmp(style, "dihedral") != 0) && (strcmp(style, "improper") != 0) && (strcmp(style, "compute") != 0) && (strcmp(style, "fix") != 0) && (strcmp(style, "region") != 0) && - (strcmp(style, "command") != 0)) { + (strcmp(style, "command") != 0) && (strcmp(style, "kspace") != 0)) { if (me == 0) utils::logmesg(lmp, "Ignoring unload: {} is not a supported plugin style\n", style); return; @@ -382,6 +391,12 @@ void plugin_unload(const char *style, const char *name, LAMMPS *lmp) auto command_map = lmp->input->command_map; auto found = command_map->find(name); if (found != command_map->end()) command_map->erase(name); + + } else if (pstyle == "kspace") { + + auto kspace_map = lmp->force->kspace_map; + auto found = kspace_map->find(name); + if (found != kspace_map->end()) kspace_map->erase(name); } // if reference count is down to zero, close DSO handle. From c52bffda9cbacddb59e9215a191a63be800861e6 Mon Sep 17 00:00:00 2001 From: Colin Denniston Date: Sun, 27 Mar 2022 11:26:23 -0400 Subject: [PATCH 2/4] remove defunct LBtype argument from lb/fluid doc --- doc/src/fix_lb_fluid.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/src/fix_lb_fluid.rst b/doc/src/fix_lb_fluid.rst index d0de406156..0191f14b1c 100644 --- a/doc/src/fix_lb_fluid.rst +++ b/doc/src/fix_lb_fluid.rst @@ -8,7 +8,7 @@ Syntax .. parsed-literal:: - fix ID group-ID lb/fluid nevery LBtype viscosity density keyword values ... + fix ID group-ID lb/fluid nevery viscosity density keyword values ... * ID, group-ID are documented in :doc:`fix ` command * lb/fluid = style name of this fix command From 63caa8bb441fd7fd6ba721ac868b15390efec314 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 28 Mar 2022 15:12:25 -0400 Subject: [PATCH 3/4] reorder so that kspace follows improper and comes before compute --- doc/src/Developer_plugins.rst | 2 +- doc/src/plugin.rst | 2 +- src/PLUGIN/plugin.cpp | 35 +++++++++++++++++------------------ 3 files changed, 19 insertions(+), 20 deletions(-) diff --git a/doc/src/Developer_plugins.rst b/doc/src/Developer_plugins.rst index c4378b6c85..96bb872929 100644 --- a/doc/src/Developer_plugins.rst +++ b/doc/src/Developer_plugins.rst @@ -68,7 +68,7 @@ Members of ``lammpsplugin_t`` * - author - String with the name and email of the author * - creator.v1 - - Pointer to factory function for pair, bond, angle, dihedral, improper, command or kspace styles + - Pointer to factory function for pair, bond, angle, dihedral, improper, kspace, or command styles * - creator.v2 - Pointer to factory function for compute, fix, or region styles * - handle diff --git a/doc/src/plugin.rst b/doc/src/plugin.rst index 8f689e885a..83eea789d1 100644 --- a/doc/src/plugin.rst +++ b/doc/src/plugin.rst @@ -17,7 +17,7 @@ Syntax *load* file = load plugin(s) from shared object in *file* *unload* style name = unload plugin *name* of style *style* - *style* = *pair* or *bond* or *angle* or *dihedral* or *improper* or *compute* or *fix* or *region* or *command* or *kspace* + *style* = *pair* or *bond* or *angle* or *dihedral* or *improper* or *kspace* or *compute* or *fix* or *region* or *command* *list* = print a list of currently loaded plugins *clear* = unload all currently loaded plugins diff --git a/src/PLUGIN/plugin.cpp b/src/PLUGIN/plugin.cpp index 761b18760f..d651a16760 100644 --- a/src/PLUGIN/plugin.cpp +++ b/src/PLUGIN/plugin.cpp @@ -208,6 +208,14 @@ void plugin_register(lammpsplugin_t *plugin, void *ptr) } (*improper_map)[plugin->name] = (Force::ImproperCreator) plugin->creator.v1; + } else if (pstyle == "kspace") { + auto kspace_map = lmp->force->kspace_map; + if (kspace_map->find(plugin->name) != kspace_map->end()) { + if (lmp->comm->me == 0) + lmp->error->warning(FLERR, "Overriding built-in kspace style {} from plugin", plugin->name); + } + (*kspace_map)[plugin->name] = (Force::KSpaceCreator) plugin->creator.v1; + } else if (pstyle == "compute") { auto compute_map = lmp->modify->compute_map; if (compute_map->find(plugin->name) != compute_map->end()) { @@ -242,15 +250,6 @@ void plugin_register(lammpsplugin_t *plugin, void *ptr) } (*command_map)[plugin->name] = (Input::CommandCreator) plugin->creator.v1; - } else if (pstyle == "kspace") { - auto kspace_map = lmp->force->kspace_map; - if (kspace_map->find(plugin->name) != kspace_map->end()) { - if (lmp->comm->me == 0) - lmp->error->warning(FLERR, "Overriding built-in kspace style {} from plugin", - plugin->name); - } - (*kspace_map)[plugin->name] = (Force::KSpaceCreator) plugin->creator.v1; - } else { utils::logmesg(lmp, "Loading plugins for {} styles not yet implemented\n", pstyle); pluginlist.pop_back(); @@ -272,9 +271,9 @@ void plugin_unload(const char *style, const char *name, LAMMPS *lmp) // ignore unload request from unsupported style categories if ((strcmp(style, "pair") != 0) && (strcmp(style, "bond") != 0) && (strcmp(style, "angle") != 0) && (strcmp(style, "dihedral") != 0) && - (strcmp(style, "improper") != 0) && (strcmp(style, "compute") != 0) && - (strcmp(style, "fix") != 0) && (strcmp(style, "region") != 0) && - (strcmp(style, "command") != 0) && (strcmp(style, "kspace") != 0)) { + (strcmp(style, "improper") != 0) && (strcmp(style, "kspace") != 0) && + (strcmp(style, "compute") != 0) && (strcmp(style, "fix") != 0) && + (strcmp(style, "region") != 0) && (strcmp(style, "command") != 0)) { if (me == 0) utils::logmesg(lmp, "Ignoring unload: {} is not a supported plugin style\n", style); return; @@ -356,6 +355,12 @@ void plugin_unload(const char *style, const char *name, LAMMPS *lmp) if ((lmp->force->improper_style != nullptr) && (lmp->force->improper_match(name) != nullptr)) lmp->force->create_improper("none", 0); + } else if (pstyle == "kspace") { + + auto kspace_map = lmp->force->kspace_map; + auto found = kspace_map->find(name); + if (found != kspace_map->end()) kspace_map->erase(name); + } else if (pstyle == "compute") { auto compute_map = lmp->modify->compute_map; @@ -391,12 +396,6 @@ void plugin_unload(const char *style, const char *name, LAMMPS *lmp) auto command_map = lmp->input->command_map; auto found = command_map->find(name); if (found != command_map->end()) command_map->erase(name); - - } else if (pstyle == "kspace") { - - auto kspace_map = lmp->force->kspace_map; - auto found = kspace_map->find(name); - if (found != kspace_map->end()) kspace_map->erase(name); } // if reference count is down to zero, close DSO handle. From 438cba3604d9257513f151bebaac41f99e2a81cf Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 28 Mar 2022 16:26:55 -0400 Subject: [PATCH 4/4] update programming style to latest conventions, enable and apply clang-format --- src/FEP/compute_fep.cpp | 278 +++++++++++++++++----------------------- 1 file changed, 118 insertions(+), 160 deletions(-) diff --git a/src/FEP/compute_fep.cpp b/src/FEP/compute_fep.cpp index 787c8a29e6..70779d1cb2 100644 --- a/src/FEP/compute_fep.cpp +++ b/src/FEP/compute_fep.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -34,51 +33,49 @@ #include "update.h" #include "variable.h" -#include #include +#include using namespace LAMMPS_NS; -enum{PAIR,ATOM}; -enum{CHARGE}; +enum { PAIR, ATOM }; +enum { CHARGE }; /* ---------------------------------------------------------------------- */ -ComputeFEP::ComputeFEP(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg) +ComputeFEP::ComputeFEP(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) { - if (narg < 5) error->all(FLERR,"Illegal number of arguments in compute fep"); + if (narg < 5) error->all(FLERR, "Illegal number of arguments in compute fep"); scalar_flag = 0; vector_flag = 1; size_vector = 3; extvector = 0; + const int ntypes = atom->ntypes; vector = new double[size_vector]; - fepinitflag = 0; // avoid init to run entirely when called by write_data - temp_fep = utils::numeric(FLERR,arg[3],false,lmp); + temp_fep = utils::numeric(FLERR, arg[3], false, lmp); // count # of perturbations npert = 0; int iarg = 4; while (iarg < narg) { - if (strcmp(arg[iarg],"pair") == 0) { - if (iarg+6 > narg) error->all(FLERR, - "Illegal pair attribute in compute fep"); + if (strcmp(arg[iarg], "pair") == 0) { + if (iarg + 6 > narg) error->all(FLERR, "Illegal pair attribute in compute fep"); npert++; iarg += 6; - } else if (strcmp(arg[iarg],"atom") == 0) { - if (iarg+4 > narg) error->all(FLERR, - "Illegal atom attribute in compute fep"); + } else if (strcmp(arg[iarg], "atom") == 0) { + if (iarg + 4 > narg) error->all(FLERR, "Illegal atom attribute in compute fep"); npert++; iarg += 4; - } else break; + } else + break; } - if (npert == 0) error->all(FLERR,"Illegal syntax in compute fep"); + if (npert == 0) error->all(FLERR, "Illegal syntax in compute fep"); perturb = new Perturb[npert]; // parse keywords @@ -88,33 +85,34 @@ ComputeFEP::ComputeFEP(LAMMPS *lmp, int narg, char **arg) : iarg = 4; while (iarg < narg) { - if (strcmp(arg[iarg],"pair") == 0) { + if (strcmp(arg[iarg], "pair") == 0) { perturb[npert].which = PAIR; - perturb[npert].pstyle = utils::strdup(arg[iarg+1]); - perturb[npert].pparam = utils::strdup(arg[iarg+2]); - utils::bounds(FLERR,arg[iarg+3],1,atom->ntypes, - perturb[npert].ilo,perturb[npert].ihi,error); - utils::bounds(FLERR,arg[iarg+4],1,atom->ntypes, - perturb[npert].jlo,perturb[npert].jhi,error); - if (utils::strmatch(arg[iarg+5],"^v_")) { - perturb[npert].var = utils::strdup(arg[iarg+5]+2); - } else error->all(FLERR,"Illegal variable in compute fep"); + perturb[npert].pstyle = utils::strdup(arg[iarg + 1]); + perturb[npert].pparam = utils::strdup(arg[iarg + 2]); + utils::bounds(FLERR, arg[iarg + 3], 1, ntypes, perturb[npert].ilo, perturb[npert].ihi, error); + utils::bounds(FLERR, arg[iarg + 4], 1, ntypes, perturb[npert].jlo, perturb[npert].jhi, error); + if (utils::strmatch(arg[iarg + 5], "^v_")) { + perturb[npert].var = utils::strdup(arg[iarg + 5] + 2); + } else + error->all(FLERR, "Illegal variable in compute fep"); npert++; iarg += 6; - } else if (strcmp(arg[iarg],"atom") == 0) { + } else if (strcmp(arg[iarg], "atom") == 0) { perturb[npert].which = ATOM; - if (strcmp(arg[iarg+1],"charge") == 0) { + if (strcmp(arg[iarg + 1], "charge") == 0) { perturb[npert].aparam = CHARGE; chgflag = 1; - } else error->all(FLERR,"Illegal atom argument in compute fep"); - utils::bounds(FLERR,arg[iarg+2],1,atom->ntypes, - perturb[npert].ilo,perturb[npert].ihi,error); - if (utils::strmatch(arg[iarg+3],"^v_")) { - perturb[npert].var = utils::strdup(arg[iarg+3]+2); - } else error->all(FLERR,"Illegal variable in compute fep"); + } else + error->all(FLERR, "Illegal atom argument in compute fep"); + utils::bounds(FLERR, arg[iarg + 2], 1, ntypes, perturb[npert].ilo, perturb[npert].ihi, error); + if (utils::strmatch(arg[iarg + 3], "^v_")) { + perturb[npert].var = utils::strdup(arg[iarg + 3] + 2); + } else + error->all(FLERR, "Illegal variable in compute fep"); npert++; iarg += 4; - } else break; + } else + break; } // optional keywords @@ -123,24 +121,23 @@ ComputeFEP::ComputeFEP(LAMMPS *lmp, int narg, char **arg) : volumeflag = 0; while (iarg < narg) { - if (strcmp(arg[iarg],"tail") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal optional keyword in compute fep"); - tailflag = utils::logical(FLERR,arg[iarg+1],false,lmp); + if (strcmp(arg[iarg], "tail") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal optional keyword in compute fep"); + tailflag = utils::logical(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"volume") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal optional keyword in compute fep"); - volumeflag = utils::logical(FLERR,arg[iarg+1],false,lmp); + } else if (strcmp(arg[iarg], "volume") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal optional keyword in compute fep"); + volumeflag = utils::logical(FLERR, arg[iarg + 1], false, lmp); iarg += 2; } else - error->all(FLERR,"Illegal optional keyword in compute fep"); + error->all(FLERR, "Illegal optional keyword in compute fep"); } // allocate pair style arrays - int ntype = atom->ntypes; for (int m = 0; m < npert; m++) { if (perturb[m].which == PAIR) - memory->create(perturb[m].array_orig,ntype+1,ntype+1,"fep:array_orig"); + memory->create(perturb[m].array_orig, ntypes + 1, ntypes + 1, "fep:array_orig"); } // allocate space for charge, force, energy, virial arrays @@ -159,17 +156,17 @@ ComputeFEP::ComputeFEP(LAMMPS *lmp, int narg, char **arg) : ComputeFEP::~ComputeFEP() { - delete [] vector; + delete[] vector; for (int m = 0; m < npert; m++) { - delete [] perturb[m].var; + delete[] perturb[m].var; if (perturb[m].which == PAIR) { - delete [] perturb[m].pstyle; - delete [] perturb[m].pparam; + delete[] perturb[m].pstyle; + delete[] perturb[m].pparam; memory->destroy(perturb[m].array_orig); } } - delete [] perturb; + delete[] perturb; deallocate_storage(); } @@ -178,11 +175,12 @@ ComputeFEP::~ComputeFEP() void ComputeFEP::init() { - int i,j; + int i, j; if (!fepinitflag) // avoid init to run entirely when called by write_data - fepinitflag = 1; - else return; + fepinitflag = 1; + else + return; // setup and error checks @@ -192,49 +190,49 @@ void ComputeFEP::init() Perturb *pert = &perturb[m]; pert->ivar = input->variable->find(pert->var); - if (pert->ivar < 0) - error->all(FLERR,"Variable name for compute fep does not exist"); + if (pert->ivar < 0) error->all(FLERR, "Variable name for compute fep does not exist"); if (!input->variable->equalstyle(pert->ivar)) - error->all(FLERR,"Variable for compute fep is of invalid style"); + error->all(FLERR, "Variable for compute fep is of invalid style"); - if (force->pair == nullptr) - error->all(FLERR,"compute fep pair requires pair interactions"); + if (force->pair == nullptr) error->all(FLERR, "compute fep pair requires pair interactions"); if (pert->which == PAIR) { pairflag = 1; - Pair *pair = force->pair_match(pert->pstyle,1); - if (pair == nullptr) error->all(FLERR,"compute fep pair style " - "does not exist"); - void *ptr = pair->extract(pert->pparam,pert->pdim); - if (ptr == nullptr) - error->all(FLERR,"compute fep pair style param not supported"); + Pair *pair = force->pair_match(pert->pstyle, 1); + if (pair == nullptr) + error->all(FLERR, + "compute fep pair style " + "does not exist"); + void *ptr = pair->extract(pert->pparam, pert->pdim); + if (ptr == nullptr) error->all(FLERR, "compute fep pair style param not supported"); pert->array = (double **) ptr; // if pair hybrid, test that ilo,ihi,jlo,jhi are valid for sub-style - if ((strcmp(force->pair_style,"hybrid") == 0 || - strcmp(force->pair_style,"hybrid/overlay") == 0)) { + if ((strcmp(force->pair_style, "hybrid") == 0 || + strcmp(force->pair_style, "hybrid/overlay") == 0)) { PairHybrid *pair = (PairHybrid *) force->pair; for (i = pert->ilo; i <= pert->ihi; i++) - for (j = MAX(pert->jlo,i); j <= pert->jhi; j++) - if (!pair->check_ijtype(i,j,pert->pstyle)) - error->all(FLERR,"compute fep type pair range is not valid for " + for (j = MAX(pert->jlo, i); j <= pert->jhi; j++) + if (!pair->check_ijtype(i, j, pert->pstyle)) + error->all(FLERR, + "compute fep type pair range is not valid for " "pair hybrid sub-style"); } } else if (pert->which == ATOM) { if (pert->aparam == CHARGE) { - if (!atom->q_flag) - error->all(FLERR,"compute fep requires atom attribute charge"); + if (!atom->q_flag) error->all(FLERR, "compute fep requires atom attribute charge"); } } } if (tailflag) { if (force->pair->tail_flag == 0) - error->all(FLERR,"Compute fep tail when pair style does not " + error->all(FLERR, + "Compute fep tail when pair style does not " "compute tail corrections"); } @@ -244,65 +242,46 @@ void ComputeFEP::init() if (ifixgpu >= 0) fixgpu = modify->fix[ifixgpu]; if (comm->me == 0) { - if (screen) { - fprintf(screen, "FEP settings ...\n"); - fprintf(screen, " temperature = %f\n", temp_fep); - fprintf(screen, " tail %s\n", (tailflag ? "yes":"no")); - for (int m = 0; m < npert; m++) { - Perturb *pert = &perturb[m]; - if (pert->which == PAIR) - fprintf(screen, " pair %s %s %d-%d %d-%d\n", pert->pstyle, - pert->pparam, - pert->ilo, pert->ihi, pert->jlo, pert->jhi); - else if (pert->which == ATOM) - fprintf(screen, " atom charge %d-%d\n", pert->ilo, pert->ihi); - } - } - if (logfile) { - fprintf(logfile, "FEP settings ...\n"); - fprintf(logfile, " temperature = %f\n", temp_fep); - fprintf(logfile, " tail %s\n", (tailflag ? "yes":"no")); - for (int m = 0; m < npert; m++) { - Perturb *pert = &perturb[m]; - if (pert->which == PAIR) - fprintf(logfile, " pair %s %s %d-%d %d-%d\n", pert->pstyle, - pert->pparam, - pert->ilo, pert->ihi, pert->jlo, pert->jhi); - else if (pert->which == ATOM) - fprintf(logfile, " atom charge %d-%d\n", pert->ilo, pert->ihi); - } + auto mesg = fmt::format("FEP settings ...\n temperature = {:f}\n", temp_fep); + mesg += fmt::format(" tail {}\n", (tailflag ? "yes" : "no")); + for (int m = 0; m < npert; m++) { + Perturb *pert = &perturb[m]; + if (pert->which == PAIR) + mesg += fmt::format(" pair {} {} {}-{} {}-{}\n", pert->pstyle, pert->pparam, pert->ilo, + pert->ihi, pert->jlo, pert->jhi); + else if (pert->which == ATOM) + mesg += fmt::format(" atom charge {}-{}\n", pert->ilo, pert->ihi); } + utils::logmesg(lmp, mesg); } - } /* ---------------------------------------------------------------------- */ - void ComputeFEP::compute_vector() { - double pe0,pe1; + double pe0, pe1; eflag = 1; vflag = 0; invoked_vector = update->ntimestep; - if (atom->nmax > nmax) { // reallocate working arrays if necessary + if (atom->nmax > nmax) { // reallocate working arrays if necessary deallocate_storage(); allocate_storage(); } - backup_qfev(); // backup charge, force, energy, virial array values - backup_params(); // backup pair parameters + backup_qfev(); // backup charge, force, energy, virial array values + backup_params(); // backup pair parameters timer->stamp(); if (force->pair && force->pair->compute_flag) { - force->pair->compute(eflag,vflag); + force->pair->compute(eflag, vflag); timer->stamp(Timer::PAIR); } if (chgflag && force->kspace && force->kspace->compute_flag) { - force->kspace->compute(eflag,vflag); + force->kspace->compute(eflag, vflag); timer->stamp(Timer::KSPACE); } @@ -315,11 +294,11 @@ void ComputeFEP::compute_vector() timer->stamp(); if (force->pair && force->pair->compute_flag) { - force->pair->compute(eflag,vflag); + force->pair->compute(eflag, vflag); timer->stamp(Timer::PAIR); } if (chgflag && force->kspace && force->kspace->compute_flag) { - force->kspace->compute(eflag,vflag); + force->kspace->compute(eflag, vflag); timer->stamp(Timer::KSPACE); } @@ -330,17 +309,15 @@ void ComputeFEP::compute_vector() pe1 = compute_epair(); - restore_qfev(); // restore charge, force, energy, virial array values - restore_params(); // restore pair parameters + restore_qfev(); // restore charge, force, energy, virial array values + restore_params(); // restore pair parameters - vector[0] = pe1-pe0; - vector[1] = exp(-(pe1-pe0)/(force->boltz*temp_fep)); + vector[0] = pe1 - pe0; + vector[1] = exp(-(pe1 - pe0) / (force->boltz * temp_fep)); vector[2] = domain->xprd * domain->yprd * domain->zprd; - if (volumeflag) - vector[1] *= vector[2]; + if (volumeflag) vector[1] *= vector[2]; } - /* ---------------------------------------------------------------------- obtain pair energy from lammps accumulators ------------------------------------------------------------------------- */ @@ -350,9 +327,8 @@ double ComputeFEP::compute_epair() double eng, eng_pair; eng = 0.0; - if (force->pair) - eng = force->pair->eng_vdwl + force->pair->eng_coul; - MPI_Allreduce(&eng,&eng_pair,1,MPI_DOUBLE,MPI_SUM,world); + if (force->pair) eng = force->pair->eng_vdwl + force->pair->eng_coul; + MPI_Allreduce(&eng, &eng_pair, 1, MPI_DOUBLE, MPI_SUM, world); if (tailflag) { double volume = domain->xprd * domain->yprd * domain->zprd; @@ -364,28 +340,27 @@ double ComputeFEP::compute_epair() return eng_pair; } - /* ---------------------------------------------------------------------- apply perturbation to pair, atom parameters based on variable evaluation ------------------------------------------------------------------------- */ void ComputeFEP::perturb_params() { - int i,j; + int i, j; for (int m = 0; m < npert; m++) { Perturb *pert = &perturb[m]; double delta = input->variable->compute_equal(pert->ivar); - if (pert->which == PAIR) { // modify pair parameters + if (pert->which == PAIR) { // modify pair parameters for (i = pert->ilo; i <= pert->ihi; i++) - for (j = MAX(pert->jlo,i); j <= pert->jhi; j++) + for (j = MAX(pert->jlo, i); j <= pert->jhi; j++) pert->array[i][j] = pert->array_orig[i][j] + delta; } else if (pert->which == ATOM) { - if (pert->aparam == CHARGE) { // modify charges + if (pert->aparam == CHARGE) { // modify charges int *atype = atom->type; double *q = atom->q; int *mask = atom->mask; @@ -393,9 +368,7 @@ void ComputeFEP::perturb_params() for (i = 0; i < natom; i++) if (atype[i] >= pert->ilo && atype[i] <= pert->ihi) - if (mask[i] & groupbit) - q[i] += delta; - + if (mask[i] & groupbit) q[i] += delta; } } } @@ -411,40 +384,36 @@ void ComputeFEP::perturb_params() if (chgflag && force->kspace) force->kspace->qsum_qsq(); } - /* ---------------------------------------------------------------------- backup pair parameters ------------------------------------------------------------------------- */ void ComputeFEP::backup_params() { - int i,j; + int i, j; for (int m = 0; m < npert; m++) { Perturb *pert = &perturb[m]; if (pert->which == PAIR) { for (i = pert->ilo; i <= pert->ihi; i++) - for (j = MAX(pert->jlo,i); j <= pert->jhi; j++) - pert->array_orig[i][j] = pert->array[i][j]; + for (j = MAX(pert->jlo, i); j <= pert->jhi; j++) pert->array_orig[i][j] = pert->array[i][j]; } } } - /* ---------------------------------------------------------------------- restore pair parameters to original values ------------------------------------------------------------------------- */ void ComputeFEP::restore_params() { - int i,j; + int i, j; for (int m = 0; m < npert; m++) { Perturb *pert = &perturb[m]; if (pert->which == PAIR) { for (i = pert->ilo; i <= pert->ihi; i++) - for (j = MAX(pert->jlo,i); j <= pert->jhi; j++) - pert->array[i][j] = pert->array_orig[i][j]; + for (j = MAX(pert->jlo, i); j <= pert->jhi; j++) pert->array[i][j] = pert->array_orig[i][j]; } } @@ -455,7 +424,6 @@ void ComputeFEP::restore_params() if (chgflag && force->kspace) force->kspace->qsum_qsq(); } - /* ---------------------------------------------------------------------- manage storage for charge, force, energy, virial arrays ------------------------------------------------------------------------- */ @@ -463,14 +431,14 @@ void ComputeFEP::restore_params() void ComputeFEP::allocate_storage() { nmax = atom->nmax; - memory->create(f_orig,nmax,3,"fep:f_orig"); - memory->create(peatom_orig,nmax,"fep:peatom_orig"); - memory->create(pvatom_orig,nmax,6,"fep:pvatom_orig"); + memory->create(f_orig, nmax, 3, "fep:f_orig"); + memory->create(peatom_orig, nmax, "fep:peatom_orig"); + memory->create(pvatom_orig, nmax, 6, "fep:pvatom_orig"); if (chgflag) { - memory->create(q_orig,nmax,"fep:q_orig"); + memory->create(q_orig, nmax, "fep:q_orig"); if (force->kspace) { - memory->create(keatom_orig,nmax,"fep:keatom_orig"); - memory->create(kvatom_orig,nmax,6,"fep:kvatom_orig"); + memory->create(keatom_orig, nmax, "fep:keatom_orig"); + memory->create(kvatom_orig, nmax, 6, "fep:kvatom_orig"); } } } @@ -492,7 +460,6 @@ void ComputeFEP::deallocate_storage() pvatom_orig = kvatom_orig = nullptr; } - /* ---------------------------------------------------------------------- backup and restore arrays with charge, force, energy, virial ------------------------------------------------------------------------- */ @@ -503,8 +470,7 @@ void ComputeFEP::backup_qfev() int nall = atom->nlocal + atom->nghost; int natom = atom->nlocal; - if (force->newton || force->kspace->tip4pflag) - natom += atom->nghost; + if (force->newton || force->kspace->tip4pflag) natom += atom->nghost; double **f = atom->f; for (i = 0; i < natom; i++) { @@ -525,8 +491,7 @@ void ComputeFEP::backup_qfev() if (update->eflag_atom) { double *peatom = force->pair->eatom; - for (i = 0; i < natom; i++) - peatom_orig[i] = peatom[i]; + for (i = 0; i < natom; i++) peatom_orig[i] = peatom[i]; } if (update->vflag_atom) { double **pvatom = force->pair->vatom; @@ -542,8 +507,7 @@ void ComputeFEP::backup_qfev() if (chgflag) { double *q = atom->q; - for (i = 0; i < nall; i++) - q_orig[i] = q[i]; + for (i = 0; i < nall; i++) q_orig[i] = q[i]; if (force->kspace) { energy_orig = force->kspace->energy; @@ -556,8 +520,7 @@ void ComputeFEP::backup_qfev() if (update->eflag_atom) { double *keatom = force->kspace->eatom; - for (i = 0; i < natom; i++) - keatom_orig[i] = keatom[i]; + for (i = 0; i < natom; i++) keatom_orig[i] = keatom[i]; } if (update->vflag_atom) { double **kvatom = force->kspace->vatom; @@ -582,8 +545,7 @@ void ComputeFEP::restore_qfev() int nall = atom->nlocal + atom->nghost; int natom = atom->nlocal; - if (force->newton || force->kspace->tip4pflag) - natom += atom->nghost; + if (force->newton || force->kspace->tip4pflag) natom += atom->nghost; double **f = atom->f; for (i = 0; i < natom; i++) { @@ -604,8 +566,7 @@ void ComputeFEP::restore_qfev() if (update->eflag_atom) { double *peatom = force->pair->eatom; - for (i = 0; i < natom; i++) - peatom[i] = peatom_orig[i]; + for (i = 0; i < natom; i++) peatom[i] = peatom_orig[i]; } if (update->vflag_atom) { double **pvatom = force->pair->vatom; @@ -621,8 +582,7 @@ void ComputeFEP::restore_qfev() if (chgflag) { double *q = atom->q; - for (i = 0; i < nall; i++) - q[i] = q_orig[i]; + for (i = 0; i < nall; i++) q[i] = q_orig[i]; if (force->kspace) { force->kspace->energy = energy_orig; @@ -635,8 +595,7 @@ void ComputeFEP::restore_qfev() if (update->eflag_atom) { double *keatom = force->kspace->eatom; - for (i = 0; i < natom; i++) - keatom[i] = keatom_orig[i]; + for (i = 0; i < natom; i++) keatom[i] = keatom_orig[i]; } if (update->vflag_atom) { double **kvatom = force->kspace->vatom; @@ -652,4 +611,3 @@ void ComputeFEP::restore_qfev() } } } -