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

update programming style to latest conventions, enable and apply clang-format

Browse Source
This commit is contained in:
Axel Kohlmeyer
2022-03-28 16:26:55 -04:00
parent 63caa8bb44
commit 438cba3604
1 changed files with 118 additions and 160 deletions
Download Patch File Download Diff File Expand all files Collapse all files

144
src/FEP/compute_fep.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories https://www.lammps.org/, Sandia National Laboratories
@ -34,8 +33,8 @@
#include "update.h" #include "update.h"
#include "variable.h" #include "variable.h"
#include <cstring>
#include <cmath> #include <cmath>
#include <cstring>
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
@ -44,8 +43,7 @@ enum{CHARGE};
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
ComputeFEP::ComputeFEP(LAMMPS *lmp, int narg, char **arg) : ComputeFEP::ComputeFEP(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, 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");
@ -54,8 +52,8 @@ ComputeFEP::ComputeFEP(LAMMPS *lmp, int narg, char **arg) :
size_vector = 3; size_vector = 3;
extvector = 0; extvector = 0;
const int ntypes = atom->ntypes;
vector = new double[size_vector]; vector = new double[size_vector];
fepinitflag = 0; // avoid init to run entirely when called by write_data 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);
@ -66,16 +64,15 @@ ComputeFEP::ComputeFEP(LAMMPS *lmp, int narg, char **arg) :
int iarg = 4; int iarg = 4;
while (iarg < narg) { while (iarg < narg) {
if (strcmp(arg[iarg], "pair") == 0) { if (strcmp(arg[iarg], "pair") == 0) {
if (iarg+6 > narg) error->all(FLERR, if (iarg + 6 > narg) error->all(FLERR, "Illegal pair attribute in compute fep");
"Illegal pair attribute in compute fep");
npert++; npert++;
iarg += 6; iarg += 6;
} else if (strcmp(arg[iarg], "atom") == 0) { } else if (strcmp(arg[iarg], "atom") == 0) {
if (iarg+4 > narg) error->all(FLERR, if (iarg + 4 > narg) error->all(FLERR, "Illegal atom attribute in compute fep");
"Illegal atom attribute in compute fep");
npert++; npert++;
iarg += 4; 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");
@ -92,13 +89,12 @@ ComputeFEP::ComputeFEP(LAMMPS *lmp, int narg, char **arg) :
perturb[npert].which = PAIR; perturb[npert].which = PAIR;
perturb[npert].pstyle = utils::strdup(arg[iarg + 1]); perturb[npert].pstyle = utils::strdup(arg[iarg + 1]);
perturb[npert].pparam = utils::strdup(arg[iarg + 2]); perturb[npert].pparam = utils::strdup(arg[iarg + 2]);
utils::bounds(FLERR,arg[iarg+3],1,atom->ntypes, utils::bounds(FLERR, arg[iarg + 3], 1, ntypes, perturb[npert].ilo, perturb[npert].ihi, error);
perturb[npert].ilo,perturb[npert].ihi,error); utils::bounds(FLERR, arg[iarg + 4], 1, ntypes, perturb[npert].jlo, perturb[npert].jhi, error);
utils::bounds(FLERR,arg[iarg+4],1,atom->ntypes,
perturb[npert].jlo,perturb[npert].jhi,error);
if (utils::strmatch(arg[iarg + 5], "^v_")) { if (utils::strmatch(arg[iarg + 5], "^v_")) {
perturb[npert].var = utils::strdup(arg[iarg + 5] + 2); perturb[npert].var = utils::strdup(arg[iarg + 5] + 2);
} else error->all(FLERR,"Illegal variable in compute fep"); } else
error->all(FLERR, "Illegal variable in compute fep");
npert++; npert++;
iarg += 6; iarg += 6;
} else if (strcmp(arg[iarg], "atom") == 0) { } else if (strcmp(arg[iarg], "atom") == 0) {
@ -106,15 +102,17 @@ ComputeFEP::ComputeFEP(LAMMPS *lmp, int narg, char **arg) :
if (strcmp(arg[iarg + 1], "charge") == 0) { if (strcmp(arg[iarg + 1], "charge") == 0) {
perturb[npert].aparam = CHARGE; perturb[npert].aparam = CHARGE;
chgflag = 1; chgflag = 1;
} else error->all(FLERR,"Illegal atom argument in compute fep"); } else
utils::bounds(FLERR,arg[iarg+2],1,atom->ntypes, error->all(FLERR, "Illegal atom argument in compute fep");
perturb[npert].ilo,perturb[npert].ihi,error); utils::bounds(FLERR, arg[iarg + 2], 1, ntypes, perturb[npert].ilo, perturb[npert].ihi, error);
if (utils::strmatch(arg[iarg + 3], "^v_")) { if (utils::strmatch(arg[iarg + 3], "^v_")) {
perturb[npert].var = utils::strdup(arg[iarg + 3] + 2); perturb[npert].var = utils::strdup(arg[iarg + 3] + 2);
} else error->all(FLERR,"Illegal variable in compute fep"); } else
error->all(FLERR, "Illegal variable in compute fep");
npert++; npert++;
iarg += 4; iarg += 4;
} else break; } else
break;
} }
// optional keywords // optional keywords
@ -137,10 +135,9 @@ ComputeFEP::ComputeFEP(LAMMPS *lmp, int narg, char **arg) :
// allocate pair style arrays // allocate pair style arrays
int ntype = atom->ntypes;
for (int m = 0; m < npert; m++) { for (int m = 0; m < npert; m++) {
if (perturb[m].which == PAIR) 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 // allocate space for charge, force, energy, virial arrays
@ -182,7 +179,8 @@ void ComputeFEP::init()
if (!fepinitflag) // avoid init to run entirely when called by write_data if (!fepinitflag) // avoid init to run entirely when called by write_data
fepinitflag = 1; fepinitflag = 1;
else return; else
return;
// setup and error checks // setup and error checks
@ -192,23 +190,22 @@ void ComputeFEP::init()
Perturb *pert = &perturb[m]; Perturb *pert = &perturb[m];
pert->ivar = input->variable->find(pert->var); pert->ivar = input->variable->find(pert->var);
if (pert->ivar < 0) if (pert->ivar < 0) error->all(FLERR, "Variable name for compute fep does not exist");
error->all(FLERR,"Variable name for compute fep does not exist");
if (!input->variable->equalstyle(pert->ivar)) 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) if (force->pair == nullptr) error->all(FLERR, "compute fep pair requires pair interactions");
error->all(FLERR,"compute fep pair requires pair interactions");
if (pert->which == PAIR) { if (pert->which == PAIR) {
pairflag = 1; pairflag = 1;
Pair *pair = force->pair_match(pert->pstyle, 1); Pair *pair = force->pair_match(pert->pstyle, 1);
if (pair == nullptr) error->all(FLERR,"compute fep pair style " if (pair == nullptr)
error->all(FLERR,
"compute fep pair style "
"does not exist"); "does not exist");
void *ptr = pair->extract(pert->pparam, pert->pdim); void *ptr = pair->extract(pert->pparam, pert->pdim);
if (ptr == nullptr) if (ptr == nullptr) error->all(FLERR, "compute fep pair style param not supported");
error->all(FLERR,"compute fep pair style param not supported");
pert->array = (double **) ptr; pert->array = (double **) ptr;
@ -220,21 +217,22 @@ void ComputeFEP::init()
for (i = pert->ilo; i <= pert->ihi; i++) 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++)
if (!pair->check_ijtype(i, j, pert->pstyle)) if (!pair->check_ijtype(i, j, pert->pstyle))
error->all(FLERR,"compute fep type pair range is not valid for " error->all(FLERR,
"compute fep type pair range is not valid for "
"pair hybrid sub-style"); "pair hybrid sub-style");
} }
} else if (pert->which == ATOM) { } else if (pert->which == ATOM) {
if (pert->aparam == CHARGE) { if (pert->aparam == CHARGE) {
if (!atom->q_flag) if (!atom->q_flag) error->all(FLERR, "compute fep requires atom attribute charge");
error->all(FLERR,"compute fep requires atom attribute charge");
} }
} }
} }
if (tailflag) { if (tailflag) {
if (force->pair->tail_flag == 0) 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"); "compute tail corrections");
} }
@ -244,41 +242,22 @@ void ComputeFEP::init()
if (ifixgpu >= 0) fixgpu = modify->fix[ifixgpu]; if (ifixgpu >= 0) fixgpu = modify->fix[ifixgpu];
if (comm->me == 0) { if (comm->me == 0) {
if (screen) { auto mesg = fmt::format("FEP settings ...\n temperature = {:f}\n", temp_fep);
fprintf(screen, "FEP settings ...\n"); mesg += fmt::format(" tail {}\n", (tailflag ? "yes" : "no"));
fprintf(screen, " temperature = %f\n", temp_fep);
fprintf(screen, " tail %s\n", (tailflag ? "yes":"no"));
for (int m = 0; m < npert; m++) { for (int m = 0; m < npert; m++) {
Perturb *pert = &perturb[m]; Perturb *pert = &perturb[m];
if (pert->which == PAIR) if (pert->which == PAIR)
fprintf(screen, " pair %s %s %d-%d %d-%d\n", pert->pstyle, mesg += fmt::format(" pair {} {} {}-{} {}-{}\n", pert->pstyle, pert->pparam, pert->ilo,
pert->pparam, pert->ihi, pert->jlo, pert->jhi);
pert->ilo, pert->ihi, pert->jlo, pert->jhi);
else if (pert->which == ATOM) else if (pert->which == ATOM)
fprintf(screen, " atom charge %d-%d\n", pert->ilo, pert->ihi); mesg += fmt::format(" atom charge {}-{}\n", pert->ilo, pert->ihi);
} }
utils::logmesg(lmp, mesg);
} }
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);
}
}
}
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void ComputeFEP::compute_vector() void ComputeFEP::compute_vector()
{ {
double pe0, pe1; double pe0, pe1;
@ -336,11 +315,9 @@ void ComputeFEP::compute_vector()
vector[0] = pe1 - pe0; vector[0] = pe1 - pe0;
vector[1] = exp(-(pe1 - pe0) / (force->boltz * temp_fep)); vector[1] = exp(-(pe1 - pe0) / (force->boltz * temp_fep));
vector[2] = domain->xprd * domain->yprd * domain->zprd; vector[2] = domain->xprd * domain->yprd * domain->zprd;
if (volumeflag) if (volumeflag) vector[1] *= vector[2];
vector[1] *= vector[2];
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
obtain pair energy from lammps accumulators obtain pair energy from lammps accumulators
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
@ -350,8 +327,7 @@ double ComputeFEP::compute_epair()
double eng, eng_pair; double eng, eng_pair;
eng = 0.0; eng = 0.0;
if (force->pair) if (force->pair) eng = force->pair->eng_vdwl + force->pair->eng_coul;
eng = force->pair->eng_vdwl + force->pair->eng_coul;
MPI_Allreduce(&eng, &eng_pair, 1, MPI_DOUBLE, MPI_SUM, world); MPI_Allreduce(&eng, &eng_pair, 1, MPI_DOUBLE, MPI_SUM, world);
if (tailflag) { if (tailflag) {
@ -364,7 +340,6 @@ double ComputeFEP::compute_epair()
return eng_pair; return eng_pair;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
apply perturbation to pair, atom parameters based on variable evaluation apply perturbation to pair, atom parameters based on variable evaluation
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
@ -393,9 +368,7 @@ void ComputeFEP::perturb_params()
for (i = 0; i < natom; i++) for (i = 0; i < natom; i++)
if (atype[i] >= pert->ilo && atype[i] <= pert->ihi) if (atype[i] >= pert->ilo && atype[i] <= pert->ihi)
if (mask[i] & groupbit) if (mask[i] & groupbit) q[i] += delta;
q[i] += delta;
} }
} }
} }
@ -411,7 +384,6 @@ void ComputeFEP::perturb_params()
if (chgflag && force->kspace) force->kspace->qsum_qsq(); if (chgflag && force->kspace) force->kspace->qsum_qsq();
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
backup pair parameters backup pair parameters
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
@ -424,13 +396,11 @@ void ComputeFEP::backup_params()
Perturb *pert = &perturb[m]; Perturb *pert = &perturb[m];
if (pert->which == PAIR) { if (pert->which == PAIR) {
for (i = pert->ilo; i <= pert->ihi; i++) 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_orig[i][j] = pert->array[i][j];
pert->array_orig[i][j] = pert->array[i][j];
} }
} }
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
restore pair parameters to original values restore pair parameters to original values
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
@ -443,8 +413,7 @@ void ComputeFEP::restore_params()
Perturb *pert = &perturb[m]; Perturb *pert = &perturb[m];
if (pert->which == PAIR) { if (pert->which == PAIR) {
for (i = pert->ilo; i <= pert->ihi; i++) 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];
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(); if (chgflag && force->kspace) force->kspace->qsum_qsq();
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
manage storage for charge, force, energy, virial arrays manage storage for charge, force, energy, virial arrays
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
@ -492,7 +460,6 @@ void ComputeFEP::deallocate_storage()
pvatom_orig = kvatom_orig = nullptr; pvatom_orig = kvatom_orig = nullptr;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
backup and restore arrays with charge, force, energy, virial backup and restore arrays with charge, force, energy, virial
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
@ -503,8 +470,7 @@ void ComputeFEP::backup_qfev()
int nall = atom->nlocal + atom->nghost; int nall = atom->nlocal + atom->nghost;
int natom = atom->nlocal; int natom = atom->nlocal;
if (force->newton || force->kspace->tip4pflag) if (force->newton || force->kspace->tip4pflag) natom += atom->nghost;
natom += atom->nghost;
double **f = atom->f; double **f = atom->f;
for (i = 0; i < natom; i++) { for (i = 0; i < natom; i++) {
@ -525,8 +491,7 @@ void ComputeFEP::backup_qfev()
if (update->eflag_atom) { if (update->eflag_atom) {
double *peatom = force->pair->eatom; double *peatom = force->pair->eatom;
for (i = 0; i < natom; i++) for (i = 0; i < natom; i++) peatom_orig[i] = peatom[i];
peatom_orig[i] = peatom[i];
} }
if (update->vflag_atom) { if (update->vflag_atom) {
double **pvatom = force->pair->vatom; double **pvatom = force->pair->vatom;
@ -542,8 +507,7 @@ void ComputeFEP::backup_qfev()
if (chgflag) { if (chgflag) {
double *q = atom->q; double *q = atom->q;
for (i = 0; i < nall; i++) for (i = 0; i < nall; i++) q_orig[i] = q[i];
q_orig[i] = q[i];
if (force->kspace) { if (force->kspace) {
energy_orig = force->kspace->energy; energy_orig = force->kspace->energy;
@ -556,8 +520,7 @@ void ComputeFEP::backup_qfev()
if (update->eflag_atom) { if (update->eflag_atom) {
double *keatom = force->kspace->eatom; double *keatom = force->kspace->eatom;
for (i = 0; i < natom; i++) for (i = 0; i < natom; i++) keatom_orig[i] = keatom[i];
keatom_orig[i] = keatom[i];
} }
if (update->vflag_atom) { if (update->vflag_atom) {
double **kvatom = force->kspace->vatom; double **kvatom = force->kspace->vatom;
@ -582,8 +545,7 @@ void ComputeFEP::restore_qfev()
int nall = atom->nlocal + atom->nghost; int nall = atom->nlocal + atom->nghost;
int natom = atom->nlocal; int natom = atom->nlocal;
if (force->newton || force->kspace->tip4pflag) if (force->newton || force->kspace->tip4pflag) natom += atom->nghost;
natom += atom->nghost;
double **f = atom->f; double **f = atom->f;
for (i = 0; i < natom; i++) { for (i = 0; i < natom; i++) {
@ -604,8 +566,7 @@ void ComputeFEP::restore_qfev()
if (update->eflag_atom) { if (update->eflag_atom) {
double *peatom = force->pair->eatom; double *peatom = force->pair->eatom;
for (i = 0; i < natom; i++) for (i = 0; i < natom; i++) peatom[i] = peatom_orig[i];
peatom[i] = peatom_orig[i];
} }
if (update->vflag_atom) { if (update->vflag_atom) {
double **pvatom = force->pair->vatom; double **pvatom = force->pair->vatom;
@ -621,8 +582,7 @@ void ComputeFEP::restore_qfev()
if (chgflag) { if (chgflag) {
double *q = atom->q; double *q = atom->q;
for (i = 0; i < nall; i++) for (i = 0; i < nall; i++) q[i] = q_orig[i];
q[i] = q_orig[i];
if (force->kspace) { if (force->kspace) {
force->kspace->energy = energy_orig; force->kspace->energy = energy_orig;
@ -635,8 +595,7 @@ void ComputeFEP::restore_qfev()
if (update->eflag_atom) { if (update->eflag_atom) {
double *keatom = force->kspace->eatom; double *keatom = force->kspace->eatom;
for (i = 0; i < natom; i++) for (i = 0; i < natom; i++) keatom[i] = keatom_orig[i];
keatom[i] = keatom_orig[i];
} }
if (update->vflag_atom) { if (update->vflag_atom) {
double **kvatom = force->kspace->vatom; double **kvatom = force->kspace->vatom;
@ -652,4 +611,3 @@ void ComputeFEP::restore_qfev()
} }
} }
} }