From 7b818ace888aa51735d449523cdbbf9bc7b86d89 Mon Sep 17 00:00:00 2001 From: Thomas Swinburne Date: Sun, 4 Dec 2022 17:20:41 +0100 Subject: [PATCH 01/25] NEARLY working E_neb --- src/REPLICA/fix_neb.cpp | 85 ++++++++++++++++++++++++++++++++++++----- src/REPLICA/fix_neb.h | 4 +- 2 files changed, 78 insertions(+), 11 deletions(-) diff --git a/src/REPLICA/fix_neb.cpp b/src/REPLICA/fix_neb.cpp index 3ebf64d609..e627778bd6 100644 --- a/src/REPLICA/fix_neb.cpp +++ b/src/REPLICA/fix_neb.cpp @@ -33,6 +33,7 @@ #include #include +#include using namespace LAMMPS_NS; using namespace FixConst; @@ -46,7 +47,7 @@ enum{SINGLE_PROC_DIRECT,SINGLE_PROC_MAP,MULTI_PROC}; FixNEB::FixNEB(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), - id_pe(nullptr), pe(nullptr), nlenall(nullptr), xprev(nullptr), xnext(nullptr), + id_pe(nullptr), pe(nullptr), nlenall(nullptr), vengall(nullptr), xprev(nullptr), xnext(nullptr), fnext(nullptr), springF(nullptr), tangent(nullptr), xsend(nullptr), xrecv(nullptr), fsend(nullptr), frecv(nullptr), tagsend(nullptr), tagrecv(nullptr), xsendall(nullptr), xrecvall(nullptr), fsendall(nullptr), frecvall(nullptr), @@ -62,6 +63,7 @@ FixNEB::FixNEB(LAMMPS *lmp, int narg, char **arg) : // optional params NEBLongRange = false; + EqualForceNEB = false; StandardNEB = true; PerpSpring = FreeEndIni = FreeEndFinal = false; FreeEndFinalWithRespToEIni = FinalAndInterWithRespToEIni = false; @@ -75,9 +77,15 @@ FixNEB::FixNEB(LAMMPS *lmp, int narg, char **arg) : if (iarg+2 > narg) error->all(FLERR,"Illegal fix neb command"); if (strcmp(arg[iarg+1],"ideal") == 0) { NEBLongRange = true; + EqualForceNEB = false; + StandardNEB = false; + } else if (strcmp(arg[iarg+1],"equal") == 0) { + NEBLongRange = false; + EqualForceNEB = true; StandardNEB = false; } else if (strcmp(arg[iarg+1],"neigh") == 0) { NEBLongRange = false; + EqualForceNEB = false; StandardNEB = true; } else error->all(FLERR,"Illegal fix neb command"); iarg += 2; @@ -135,7 +143,7 @@ FixNEB::FixNEB(LAMMPS *lmp, int narg, char **arg) : else procnext = -1; uworld = universe->uworld; - if (NEBLongRange) { + if (NEBLongRange or EqualForceNEB) { int *iroots = new int[nreplica]; MPI_Group uworldgroup,rootgroup; @@ -189,10 +197,11 @@ FixNEB::~FixNEB() memory->destroy(counts); memory->destroy(displacements); - if (NEBLongRange) { + if (NEBLongRange or EqualForceNEB) { if (rootworld != MPI_COMM_NULL) MPI_Comm_free(&rootworld); memory->destroy(nlenall); } + if(EqualForceNEB) memory->destroy(vengall); } /* ---------------------------------------------------------------------- */ @@ -546,8 +555,10 @@ void FixNEB::min_post_force(int /*vflag*/) double lentot = 0; double meanDist,idealPos,lenuntilIm,lenuntilClimber; + double meanDistBeforeClimber,meanDistAfterClimber; lenuntilClimber=0; - if (NEBLongRange) { + + if (NEBLongRange or EqualForceNEB) { if (cmode == SINGLE_PROC_DIRECT || cmode == SINGLE_PROC_MAP) { MPI_Allgather(&nlen,1,MPI_DOUBLE,&nlenall[0],1,MPI_DOUBLE,uworld); } else { @@ -568,9 +579,8 @@ void FixNEB::min_post_force(int /*vflag*/) if (rclimber>0) { for (int i = 0; i < rclimber; i++) lenuntilClimber += nlenall[i]; - double meanDistBeforeClimber = lenuntilClimber/rclimber; - double meanDistAfterClimber = - (lentot-lenuntilClimber)/(nreplica-rclimber-1); + meanDistBeforeClimber = lenuntilClimber/rclimber; + meanDistAfterClimber = (lentot-lenuntilClimber)/(nreplica-rclimber-1); if (ireplica0) { + for (int i = 0; i < rclimber-1; i++) + ElenuntilClimber += std::abs(vengall[i+1]-vengall[i]); + + meanEDistBeforeClimber = ElenuntilClimber/rclimber; + meanEDistAfterClimber = (Elentot-ElenuntilClimber)/(nreplica-rclimber-1); + if (ireplicacreate(tagrecv,maxlocal,"neb:tagrecv"); } - if (NEBLongRange) { + if (NEBLongRange or EqualForceNEB) { memory->destroy(nlenall); memory->create(nlenall,nreplica,"neb:nlenall"); } + if (EqualForceNEB) { + memory->destroy(vengall); + memory->create(vengall,nreplica,"neb:vengall"); + } } diff --git a/src/REPLICA/fix_neb.h b/src/REPLICA/fix_neb.h index a5aa831163..6cfc2ed110 100644 --- a/src/REPLICA/fix_neb.h +++ b/src/REPLICA/fix_neb.h @@ -39,7 +39,7 @@ class FixNEB : public Fix { private: int me, nprocs, nprocs_universe; double kspring, kspringIni, kspringFinal, kspringPerp, EIniIni, EFinalIni; - bool StandardNEB, NEBLongRange, PerpSpring, FreeEndIni, FreeEndFinal; + bool StandardNEB, NEBLongRange, EqualForceNEB, PerpSpring, FreeEndIni, FreeEndFinal; bool FreeEndFinalWithRespToEIni, FinalAndInterWithRespToEIni; int ireplica, nreplica; int procnext, procprev; @@ -53,7 +53,7 @@ class FixNEB : public Fix { int nebatoms; int ntotal; // total # of atoms, NEB or not int maxlocal; // size of xprev,xnext,tangent arrays - double *nlenall; + double *nlenall, *vengall; double **xprev, **xnext, **fnext, **springF; double **tangent; double **xsend, **xrecv; // coords to send/recv to/from other replica From 0d2e0f5b3687cf759aa37fb9f581132c386c0dee Mon Sep 17 00:00:00 2001 From: Thomas Swinburne Date: Mon, 5 Dec 2022 13:58:05 +0100 Subject: [PATCH 02/25] working ? --- src/REPLICA/fix_neb.cpp | 163 ++++++++++++++++------------------------ src/REPLICA/fix_neb.h | 5 +- 2 files changed, 69 insertions(+), 99 deletions(-) diff --git a/src/REPLICA/fix_neb.cpp b/src/REPLICA/fix_neb.cpp index e627778bd6..662be541fa 100644 --- a/src/REPLICA/fix_neb.cpp +++ b/src/REPLICA/fix_neb.cpp @@ -64,6 +64,7 @@ FixNEB::FixNEB(LAMMPS *lmp, int narg, char **arg) : NEBLongRange = false; EqualForceNEB = false; + EqualForceNEBDone = false; StandardNEB = true; PerpSpring = FreeEndIni = FreeEndFinal = false; FreeEndFinalWithRespToEIni = FinalAndInterWithRespToEIni = false; @@ -313,15 +314,7 @@ void FixNEB::min_post_force(int /*vflag*/) } if (FreeEndIni && ireplica == 0 && (update->ntimestep == 0)) EIniIni = veng; - /* if (FreeEndIni && ireplica == 0) { - // if (me == 0 ) - if (update->ntimestep == 0) { - EIniIni = veng; - // if (cmode == MULTI_PROC) - // MPI_Bcast(&EIniIni,1,MPI_DOUBLE,0,world); - } - }*/ - + // communicate atoms to/from adjacent replicas to fill xprev,xnext inter_replica_comm(); @@ -553,93 +546,7 @@ void FixNEB::min_post_force(int /*vflag*/) } } - double lentot = 0; - double meanDist,idealPos,lenuntilIm,lenuntilClimber; - double meanDistBeforeClimber,meanDistAfterClimber; - lenuntilClimber=0; - - if (NEBLongRange or EqualForceNEB) { - if (cmode == SINGLE_PROC_DIRECT || cmode == SINGLE_PROC_MAP) { - MPI_Allgather(&nlen,1,MPI_DOUBLE,&nlenall[0],1,MPI_DOUBLE,uworld); - } else { - if (me == 0) - MPI_Allgather(&nlen,1,MPI_DOUBLE,&nlenall[0],1,MPI_DOUBLE,rootworld); - MPI_Bcast(nlenall,nreplica,MPI_DOUBLE,0,world); - } - - lenuntilIm = 0; - for (int i = 0; i < ireplica; i++) - lenuntilIm += nlenall[i]; - - for (int i = 0; i < nreplica; i++) - lentot += nlenall[i]; - - meanDist = lentot/(nreplica -1); - - if (rclimber>0) { - for (int i = 0; i < rclimber; i++) - lenuntilClimber += nlenall[i]; - meanDistBeforeClimber = lenuntilClimber/rclimber; - meanDistAfterClimber = (lentot-lenuntilClimber)/(nreplica-rclimber-1); - if (ireplica0) { - for (int i = 0; i < rclimber-1; i++) - ElenuntilClimber += std::abs(vengall[i+1]-vengall[i]); - - meanEDistBeforeClimber = ElenuntilClimber/rclimber; - meanEDistAfterClimber = (Elentot-ElenuntilClimber)/(nreplica-rclimber-1); - if (ireplica0) { + for (int i = 0; i < rclimber; i++) + lenuntilClimber += nlenall[i]; + meanDistBeforeClimber = lenuntilClimber/rclimber; + meanDistAfterClimber = (lentot-lenuntilClimber)/(nreplica-rclimber-1); + if (ireplica0 and !EqualForceNEBDone) { + double lenEtot = 0; + + if (cmode == SINGLE_PROC_DIRECT || cmode == SINGLE_PROC_MAP) { + MPI_Allgather(&veng,1,MPI_DOUBLE,&vengall[0],1,MPI_DOUBLE,uworld); + } else { + if (me == 0) + MPI_Allgather(&veng,1,MPI_DOUBLE,&vengall[0],1,MPI_DOUBLE,rootworld); + MPI_Bcast(vengall,nreplica,MPI_DOUBLE,0,world); + } + + actualPos = 0; + for (int i = 0; i < ireplica; i++) + actualPos += std::abs(nlenall[i+1]-nlenall[i]); + for (int i = 0; i < nreplica -1; i++) + lenEtot += std::abs(nlenall[i+1]-nlenall[i]); + + actualPos *= lentot/lenEtot; + EqualForceNEBDone = true; + } + +} + /* ---------------------------------------------------------------------- reallocate xprev,xnext,tangent arrays if necessary reallocate communication arrays if necessary diff --git a/src/REPLICA/fix_neb.h b/src/REPLICA/fix_neb.h index 6cfc2ed110..962b053414 100644 --- a/src/REPLICA/fix_neb.h +++ b/src/REPLICA/fix_neb.h @@ -38,8 +38,8 @@ class FixNEB : public Fix { private: int me, nprocs, nprocs_universe; - double kspring, kspringIni, kspringFinal, kspringPerp, EIniIni, EFinalIni; - bool StandardNEB, NEBLongRange, EqualForceNEB, PerpSpring, FreeEndIni, FreeEndFinal; + double kspring, kspringIni, kspringFinal, kspringPerp, EIniIni, EFinalIni, idealPos, actualPos, meanDist; + bool StandardNEB, NEBLongRange, EqualForceNEB, PerpSpring, FreeEndIni, FreeEndFinal, EqualForceNEBDone; bool FreeEndFinalWithRespToEIni, FinalAndInterWithRespToEIni; int ireplica, nreplica; int procnext, procprev; @@ -68,6 +68,7 @@ class FixNEB : public Fix { int *counts, *displacements; // used for MPI_Gather void inter_replica_comm(); + void calculate_ideal_positions(); void reallocate(); }; From 1048c26900be7c4d99d5f459238ecb7178689275 Mon Sep 17 00:00:00 2001 From: Thomas Swinburne Date: Mon, 5 Dec 2022 14:07:03 +0100 Subject: [PATCH 03/25] working but not always optimal.. --- src/REPLICA/fix_neb.cpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/REPLICA/fix_neb.cpp b/src/REPLICA/fix_neb.cpp index 662be541fa..5b5de3086e 100644 --- a/src/REPLICA/fix_neb.cpp +++ b/src/REPLICA/fix_neb.cpp @@ -860,7 +860,7 @@ void FixNEB::calculate_ideal_positions() if (EqualForceNEB and rclimber>0 and !EqualForceNEBDone) { double lenEtot = 0; - + if (cmode == SINGLE_PROC_DIRECT || cmode == SINGLE_PROC_MAP) { MPI_Allgather(&veng,1,MPI_DOUBLE,&vengall[0],1,MPI_DOUBLE,uworld); } else { @@ -872,13 +872,12 @@ void FixNEB::calculate_ideal_positions() actualPos = 0; for (int i = 0; i < ireplica; i++) actualPos += std::abs(nlenall[i+1]-nlenall[i]); - for (int i = 0; i < nreplica -1; i++) + for (int i = 0; i < nreplica-1; i++) lenEtot += std::abs(nlenall[i+1]-nlenall[i]); actualPos *= lentot/lenEtot; EqualForceNEBDone = true; } - } /* ---------------------------------------------------------------------- From d69f22176f267d08e170981ea771d2328894c03f Mon Sep 17 00:00:00 2001 From: Thomas Swinburne Date: Mon, 5 Dec 2022 15:16:40 +0100 Subject: [PATCH 04/25] a version which works in initial tests --- src/REPLICA/fix_neb.cpp | 22 +++++++------ src/REPLICA/fix_neb.h | 4 ++- src/REPLICA/neb.cpp | 70 +++++++++++++++++++++++++++++++++++++++++ 3 files changed, 86 insertions(+), 10 deletions(-) diff --git a/src/REPLICA/fix_neb.cpp b/src/REPLICA/fix_neb.cpp index 5b5de3086e..938ed956de 100644 --- a/src/REPLICA/fix_neb.cpp +++ b/src/REPLICA/fix_neb.cpp @@ -64,7 +64,6 @@ FixNEB::FixNEB(LAMMPS *lmp, int narg, char **arg) : NEBLongRange = false; EqualForceNEB = false; - EqualForceNEBDone = false; StandardNEB = true; PerpSpring = FreeEndIni = FreeEndFinal = false; FreeEndFinalWithRespToEIni = FinalAndInterWithRespToEIni = false; @@ -227,6 +226,11 @@ void FixNEB::init() rclimber = -1; + // turn off equal force mode, NEB command turns it on if running "equal" mode + + equal_force = -1; + + // nebatoms = # of atoms in fix group = atoms with inter-replica forces bigint count = group->count(igroup); @@ -576,7 +580,7 @@ void FixNEB::min_post_force(int /*vflag*/) if (ireplica == rclimber) prefactor = -2.0*dot; else { if (NEBLongRange or EqualForceNEB) { - prefactor = -dot - kspring*(actualPos-idealPos)/(2*meanDist); + prefactor = -dot - kspring*(actualPos-idealPos)/2; } else if (StandardNEB) { prefactor = -dot + kspring*(nlen-plen); } @@ -856,11 +860,12 @@ void FixNEB::calculate_ideal_positions() else idealPos = lenuntilClimber+ (ireplica-rclimber)*meanDistAfterClimber; } else idealPos = ireplica * meanDist; + idealPos /= meanDist; + actualPos /= meanDist; } - if (EqualForceNEB and rclimber>0 and !EqualForceNEBDone) { + if (EqualForceNEB and rclimber>0 and equal_force>0) { double lenEtot = 0; - if (cmode == SINGLE_PROC_DIRECT || cmode == SINGLE_PROC_MAP) { MPI_Allgather(&veng,1,MPI_DOUBLE,&vengall[0],1,MPI_DOUBLE,uworld); } else { @@ -871,12 +876,11 @@ void FixNEB::calculate_ideal_positions() actualPos = 0; for (int i = 0; i < ireplica; i++) - actualPos += std::abs(nlenall[i+1]-nlenall[i]); + actualPos += std::abs(vengall[i+1]-vengall[i]); for (int i = 0; i < nreplica-1; i++) - lenEtot += std::abs(nlenall[i+1]-nlenall[i]); - - actualPos *= lentot/lenEtot; - EqualForceNEBDone = true; + lenEtot += std::abs(vengall[i+1]-vengall[i]); + actualPos *= (nreplica-1)/lenEtot; + idealPos = ireplica; } } diff --git a/src/REPLICA/fix_neb.h b/src/REPLICA/fix_neb.h index 962b053414..2eaaf8d32e 100644 --- a/src/REPLICA/fix_neb.h +++ b/src/REPLICA/fix_neb.h @@ -28,6 +28,8 @@ class FixNEB : public Fix { public: double veng, plen, nlen, dotpath, dottangrad, gradlen, dotgrad; int rclimber; + int equal_force; + bool EqualForceNEB; FixNEB(class LAMMPS *, int, char **); ~FixNEB() override; @@ -39,7 +41,7 @@ class FixNEB : public Fix { private: int me, nprocs, nprocs_universe; double kspring, kspringIni, kspringFinal, kspringPerp, EIniIni, EFinalIni, idealPos, actualPos, meanDist; - bool StandardNEB, NEBLongRange, EqualForceNEB, PerpSpring, FreeEndIni, FreeEndFinal, EqualForceNEBDone; + bool StandardNEB, NEBLongRange, PerpSpring, FreeEndIni, FreeEndFinal; bool FreeEndFinalWithRespToEIni, FinalAndInterWithRespToEIni; int ireplica, nreplica; int procnext, procprev; diff --git a/src/REPLICA/neb.cpp b/src/REPLICA/neb.cpp index 9c5befaf36..163375b7d5 100644 --- a/src/REPLICA/neb.cpp +++ b/src/REPLICA/neb.cpp @@ -317,6 +317,7 @@ void NEB::run() update->minimize->init(); fneb->rclimber = top; update->minimize->setup(); + if (me_universe == 0) { if (uscreen) { @@ -369,6 +370,75 @@ void NEB::run() if (update->minimize->stop_condition) break; } + // If running in equal spacing, perform additional NEB cycle of n2steps + + if (fneb->EqualForceNEB) { + + if (me_universe == 0) { + if (uscreen) + fprintf(uscreen,"Setting up equal force run....\n"); + if (ulogfile) + fprintf(ulogfile,"Setting up equal force run....\n"); + } + + update->beginstep = update->firststep = update->ntimestep; + update->endstep = update->laststep = update->firststep + n2steps; + update->nsteps = n2steps; + update->max_eval = n2steps; + if (update->laststep < 0) error->universe_all(FLERR,"Too many timesteps"); + + update->minimize->init(); + fneb->equal_force = 1; + fneb->rclimber = top; + update->minimize->setup(); + + if (me_universe == 0) { + if (uscreen) { + fmt::print(uscreen," Step {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} ", + "MaxReplicaForce", "MaxAtomForce", "GradV0","GradV1","GradVc","EBF", "EBR", "RDT"); + for (int i = 1; i <= nreplica; ++i) + fmt::print(uscreen, "{:^14} {:^14} ", "RD"+std::to_string(i), "PE"+std::to_string(i)); + + if (verbose) { + for (int i = 1; i <= nreplica; ++i) { + auto idx = std::to_string(i); + fmt::print(uscreen, "{:^12}{:^12}{:^12} {:^12} {:^12}{:^12} ", + "pathangle"+idx, "angletangrad"+idx, "anglegrad"+idx, "gradV"+idx, + "RepForce"+idx, "MaxAtomForce"+idx); + } + } + fprintf(uscreen,"\n"); + } + + if (ulogfile) { + fmt::print(ulogfile," Step {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} ", + "MaxReplicaForce", "MaxAtomForce", "GradV0","GradV1","GradVc","EBF", "EBR", "RDT"); + for (int i = 1; i <= nreplica; ++i) + fmt::print(ulogfile, "{:^14} {:^14} ", "RD"+std::to_string(i), "PE"+std::to_string(i)); + + if (verbose) { + for (int i = 1; i <= nreplica; ++i) { + auto idx = std::to_string(i); + fmt::print(ulogfile, "{:^12}{:^12}{:^12} {:^12} {:^12}{:^12} ", + "pathangle"+idx, "angletangrad"+idx, "anglegrad"+idx, "gradV"+idx, + "RepForce"+idx, "MaxAtomForce"+idx); + } + } + fprintf(ulogfile,"\n"); + } + } + print_status(); + + timer->init(); + timer->barrier_start(); + + while (update->minimize->niter < n2steps) { + update->minimize->run(nevery); + print_status(); + if (update->minimize->stop_condition) break; + } + } + timer->barrier_stop(); update->minimize->cleanup(); From 978db4b737147ebbc6168276b352dc9e346b5ad9 Mon Sep 17 00:00:00 2001 From: Thomas Swinburne Date: Mon, 5 Dec 2022 15:19:33 +0100 Subject: [PATCH 05/25] even simpler, still works --- src/REPLICA/fix_neb.cpp | 8 +---- src/REPLICA/fix_neb.h | 4 +-- src/REPLICA/neb.cpp | 69 ----------------------------------------- 3 files changed, 2 insertions(+), 79 deletions(-) diff --git a/src/REPLICA/fix_neb.cpp b/src/REPLICA/fix_neb.cpp index 938ed956de..97ca00eea4 100644 --- a/src/REPLICA/fix_neb.cpp +++ b/src/REPLICA/fix_neb.cpp @@ -226,11 +226,6 @@ void FixNEB::init() rclimber = -1; - // turn off equal force mode, NEB command turns it on if running "equal" mode - - equal_force = -1; - - // nebatoms = # of atoms in fix group = atoms with inter-replica forces bigint count = group->count(igroup); @@ -864,7 +859,7 @@ void FixNEB::calculate_ideal_positions() actualPos /= meanDist; } - if (EqualForceNEB and rclimber>0 and equal_force>0) { + if (EqualForceNEB and rclimber>0) { double lenEtot = 0; if (cmode == SINGLE_PROC_DIRECT || cmode == SINGLE_PROC_MAP) { MPI_Allgather(&veng,1,MPI_DOUBLE,&vengall[0],1,MPI_DOUBLE,uworld); @@ -873,7 +868,6 @@ void FixNEB::calculate_ideal_positions() MPI_Allgather(&veng,1,MPI_DOUBLE,&vengall[0],1,MPI_DOUBLE,rootworld); MPI_Bcast(vengall,nreplica,MPI_DOUBLE,0,world); } - actualPos = 0; for (int i = 0; i < ireplica; i++) actualPos += std::abs(vengall[i+1]-vengall[i]); diff --git a/src/REPLICA/fix_neb.h b/src/REPLICA/fix_neb.h index 2eaaf8d32e..442bf31199 100644 --- a/src/REPLICA/fix_neb.h +++ b/src/REPLICA/fix_neb.h @@ -28,8 +28,6 @@ class FixNEB : public Fix { public: double veng, plen, nlen, dotpath, dottangrad, gradlen, dotgrad; int rclimber; - int equal_force; - bool EqualForceNEB; FixNEB(class LAMMPS *, int, char **); ~FixNEB() override; @@ -41,7 +39,7 @@ class FixNEB : public Fix { private: int me, nprocs, nprocs_universe; double kspring, kspringIni, kspringFinal, kspringPerp, EIniIni, EFinalIni, idealPos, actualPos, meanDist; - bool StandardNEB, NEBLongRange, PerpSpring, FreeEndIni, FreeEndFinal; + bool StandardNEB, NEBLongRange, EqualForceNEB, PerpSpring, FreeEndIni, FreeEndFinal; bool FreeEndFinalWithRespToEIni, FinalAndInterWithRespToEIni; int ireplica, nreplica; int procnext, procprev; diff --git a/src/REPLICA/neb.cpp b/src/REPLICA/neb.cpp index 163375b7d5..5ce45475ba 100644 --- a/src/REPLICA/neb.cpp +++ b/src/REPLICA/neb.cpp @@ -370,75 +370,6 @@ void NEB::run() if (update->minimize->stop_condition) break; } - // If running in equal spacing, perform additional NEB cycle of n2steps - - if (fneb->EqualForceNEB) { - - if (me_universe == 0) { - if (uscreen) - fprintf(uscreen,"Setting up equal force run....\n"); - if (ulogfile) - fprintf(ulogfile,"Setting up equal force run....\n"); - } - - update->beginstep = update->firststep = update->ntimestep; - update->endstep = update->laststep = update->firststep + n2steps; - update->nsteps = n2steps; - update->max_eval = n2steps; - if (update->laststep < 0) error->universe_all(FLERR,"Too many timesteps"); - - update->minimize->init(); - fneb->equal_force = 1; - fneb->rclimber = top; - update->minimize->setup(); - - if (me_universe == 0) { - if (uscreen) { - fmt::print(uscreen," Step {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} ", - "MaxReplicaForce", "MaxAtomForce", "GradV0","GradV1","GradVc","EBF", "EBR", "RDT"); - for (int i = 1; i <= nreplica; ++i) - fmt::print(uscreen, "{:^14} {:^14} ", "RD"+std::to_string(i), "PE"+std::to_string(i)); - - if (verbose) { - for (int i = 1; i <= nreplica; ++i) { - auto idx = std::to_string(i); - fmt::print(uscreen, "{:^12}{:^12}{:^12} {:^12} {:^12}{:^12} ", - "pathangle"+idx, "angletangrad"+idx, "anglegrad"+idx, "gradV"+idx, - "RepForce"+idx, "MaxAtomForce"+idx); - } - } - fprintf(uscreen,"\n"); - } - - if (ulogfile) { - fmt::print(ulogfile," Step {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} ", - "MaxReplicaForce", "MaxAtomForce", "GradV0","GradV1","GradVc","EBF", "EBR", "RDT"); - for (int i = 1; i <= nreplica; ++i) - fmt::print(ulogfile, "{:^14} {:^14} ", "RD"+std::to_string(i), "PE"+std::to_string(i)); - - if (verbose) { - for (int i = 1; i <= nreplica; ++i) { - auto idx = std::to_string(i); - fmt::print(ulogfile, "{:^12}{:^12}{:^12} {:^12} {:^12}{:^12} ", - "pathangle"+idx, "angletangrad"+idx, "anglegrad"+idx, "gradV"+idx, - "RepForce"+idx, "MaxAtomForce"+idx); - } - } - fprintf(ulogfile,"\n"); - } - } - print_status(); - - timer->init(); - timer->barrier_start(); - - while (update->minimize->niter < n2steps) { - update->minimize->run(nevery); - print_status(); - if (update->minimize->stop_condition) break; - } - } - timer->barrier_stop(); update->minimize->cleanup(); From 8dd9682ce2383c8829e2007deb4cfa2714b8c29f Mon Sep 17 00:00:00 2001 From: Thomas Swinburne Date: Mon, 5 Dec 2022 15:28:35 +0100 Subject: [PATCH 06/25] remove iostream --- src/REPLICA/fix_neb.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/REPLICA/fix_neb.cpp b/src/REPLICA/fix_neb.cpp index 97ca00eea4..ed4eb9c135 100644 --- a/src/REPLICA/fix_neb.cpp +++ b/src/REPLICA/fix_neb.cpp @@ -33,7 +33,6 @@ #include #include -#include using namespace LAMMPS_NS; using namespace FixConst; From 5f735467fe39ef378d6d1c624a96ab7bd7b6d89b Mon Sep 17 00:00:00 2001 From: tomswinburne Date: Mon, 5 Dec 2022 16:53:50 +0100 Subject: [PATCH 07/25] added documentation --- doc/src/fix_neb.rst | 16 +++++++++++++++- 1 file changed, 15 insertions(+), 1 deletion(-) diff --git a/doc/src/fix_neb.rst b/doc/src/fix_neb.rst index f319269220..19b12604a4 100644 --- a/doc/src/fix_neb.rst +++ b/doc/src/fix_neb.rst @@ -18,9 +18,10 @@ Syntax .. parsed-literal:: - *parallel* value = *neigh* or *ideal* + *parallel* value = *neigh* or *ideal* or *equal* *neigh* = parallel nudging force based on distance to neighbor replicas (Kspring = force/distance units) *ideal* = parallel nudging force based on interpolated ideal position (Kspring = force units) + *equal* = parallel nudging force based on interpolated ideal position before climbing, then interpolated ideal energy whilst climbing (Kspring = force units) *perp* value = *Kspring2* *Kspring2* = spring constant for perpendicular nudging force (force/distance units) *end* values = estyle Kspring3 @@ -124,6 +125,19 @@ in force units. Note that the *ideal* form of nudging can often be more effective at keeping the replicas equally spaced. +With a value of *equal* the spring force is computed as for *ideal*, +before the climbing stage, then is computed to promote equidistant +spacing in energy rather than distance: + +.. parsed-literal:: + Fnudge_parallel = -\ *Kspring* \* (ED-EDideal) / (2 \* meanEDist) + +where ED is the sum of absolute (nonnegative) energy differences +between knots, EDideal = (I-1)\*meanEdist and meanEdist is the +average absolute energy difference. This form of nudging is +intended to aid schemes which integrate forces along NEB +pathways such as :doc:`fix_pafi `. + ---------- The keyword *perp* specifies if and how a perpendicular nudging force From 3c08bbb7902ba5f9393f483dc40826aa7b8cd5da Mon Sep 17 00:00:00 2001 From: tomswinburne Date: Mon, 5 Dec 2022 16:57:39 +0100 Subject: [PATCH 08/25] added documentation --- doc/src/fix_neb.rst | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/doc/src/fix_neb.rst b/doc/src/fix_neb.rst index 19b12604a4..4cf862fdb1 100644 --- a/doc/src/fix_neb.rst +++ b/doc/src/fix_neb.rst @@ -125,17 +125,17 @@ in force units. Note that the *ideal* form of nudging can often be more effective at keeping the replicas equally spaced. -With a value of *equal* the spring force is computed as for *ideal*, +With a value of *equal* the spring force is computed as for *ideal* before the climbing stage, then is computed to promote equidistant spacing in energy rather than distance: .. parsed-literal:: Fnudge_parallel = -\ *Kspring* \* (ED-EDideal) / (2 \* meanEDist) -where ED is the sum of absolute (nonnegative) energy differences -between knots, EDideal = (I-1)\*meanEdist and meanEdist is the -average absolute energy difference. This form of nudging is -intended to aid schemes which integrate forces along NEB +where ED is the cumulative sum of absolute energy differences +|E(Ri+1)-E(Ri)| for i`. ---------- From 18af945f4d919f2eaa10db8a5b216b1ec904b784 Mon Sep 17 00:00:00 2001 From: tomswinburne Date: Mon, 5 Dec 2022 16:59:26 +0100 Subject: [PATCH 09/25] added documentation --- doc/src/fix_neb.rst | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/doc/src/fix_neb.rst b/doc/src/fix_neb.rst index 4cf862fdb1..0c6d0323d7 100644 --- a/doc/src/fix_neb.rst +++ b/doc/src/fix_neb.rst @@ -133,10 +133,11 @@ spacing in energy rather than distance: Fnudge_parallel = -\ *Kspring* \* (ED-EDideal) / (2 \* meanEDist) where ED is the cumulative sum of absolute energy differences -|E(Ri+1)-E(Ri)| for i`. +ED=sum(i`, by providing optimal +quadrature points. ---------- From 0007788e0119d83cbfd671741a3709e9275dfc40 Mon Sep 17 00:00:00 2001 From: Thomas Swinburne Date: Wed, 7 Dec 2022 14:20:45 +0100 Subject: [PATCH 10/25] cleaner, minimal changes --- src/REPLICA/fix_neb.cpp | 78 ++++++++++++++++++----------------------- 1 file changed, 35 insertions(+), 43 deletions(-) diff --git a/src/REPLICA/fix_neb.cpp b/src/REPLICA/fix_neb.cpp index ed4eb9c135..6f5d027b85 100644 --- a/src/REPLICA/fix_neb.cpp +++ b/src/REPLICA/fix_neb.cpp @@ -822,44 +822,9 @@ calculate ideal positions */ void FixNEB::calculate_ideal_positions() { - double lentot = 0; - double lenuntilClimber = 0; - double meanDistBeforeClimber,meanDistAfterClimber; + double lentot,lenuntilClimber,meanDistBeforeClimber,meanDistAfterClimber; - if (NEBLongRange or EqualForceNEB) { - if (cmode == SINGLE_PROC_DIRECT || cmode == SINGLE_PROC_MAP) { - MPI_Allgather(&nlen,1,MPI_DOUBLE,&nlenall[0],1,MPI_DOUBLE,uworld); - } else { - if (me == 0) - MPI_Allgather(&nlen,1,MPI_DOUBLE,&nlenall[0],1,MPI_DOUBLE,rootworld); - MPI_Bcast(nlenall,nreplica,MPI_DOUBLE,0,world); - } - - actualPos = 0; - for (int i = 0; i < ireplica; i++) - actualPos += nlenall[i]; - - for (int i = 0; i < nreplica; i++) - lentot += nlenall[i]; - - meanDist = lentot/(nreplica-1); - - if (rclimber>0) { - for (int i = 0; i < rclimber; i++) - lenuntilClimber += nlenall[i]; - meanDistBeforeClimber = lenuntilClimber/rclimber; - meanDistAfterClimber = (lentot-lenuntilClimber)/(nreplica-rclimber-1); - if (ireplica0) { - double lenEtot = 0; + if(EqualForceNEB and rclimber>0.) { if (cmode == SINGLE_PROC_DIRECT || cmode == SINGLE_PROC_MAP) { MPI_Allgather(&veng,1,MPI_DOUBLE,&vengall[0],1,MPI_DOUBLE,uworld); } else { @@ -867,14 +832,41 @@ void FixNEB::calculate_ideal_positions() MPI_Allgather(&veng,1,MPI_DOUBLE,&vengall[0],1,MPI_DOUBLE,rootworld); MPI_Bcast(vengall,nreplica,MPI_DOUBLE,0,world); } - actualPos = 0; - for (int i = 0; i < ireplica; i++) - actualPos += std::abs(vengall[i+1]-vengall[i]); for (int i = 0; i < nreplica-1; i++) - lenEtot += std::abs(vengall[i+1]-vengall[i]); - actualPos *= (nreplica-1)/lenEtot; - idealPos = ireplica; + nlenall[i] = std::abs(vengall[i+1]-vengall[i]); + nlenall[nreplica-1] = 0.0; + } else if(NEBLongRange or EqualForceNEB) { + if (cmode == SINGLE_PROC_DIRECT || cmode == SINGLE_PROC_MAP) { + MPI_Allgather(&nlen,1,MPI_DOUBLE,&nlenall[0],1,MPI_DOUBLE,uworld); + } else { + if (me == 0) + MPI_Allgather(&nlen,1,MPI_DOUBLE,&nlenall[0],1,MPI_DOUBLE,rootworld); + MPI_Bcast(nlenall,nreplica,MPI_DOUBLE,0,world); + } } + actualPos = 0.; + for (int i = 0; i < ireplica; i++) + actualPos += nlenall[i]; + + lentot = 0.; + for (int i = 0; i < nreplica; i++) + lentot += nlenall[i]; + + meanDist = lentot/(nreplica-1); + + lenuntilClimber = 0.; + if (rclimber>0) { + for (int i = 0; i < rclimber; i++) + lenuntilClimber += nlenall[i]; + meanDistBeforeClimber = lenuntilClimber/rclimber; + meanDistAfterClimber = (lentot-lenuntilClimber)/(nreplica-rclimber-1); + if (ireplica Date: Fri, 9 Dec 2022 16:35:51 +0100 Subject: [PATCH 11/25] small addition --- src/REPLICA/fix_neb.cpp | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/src/REPLICA/fix_neb.cpp b/src/REPLICA/fix_neb.cpp index 6f5d027b85..e512c88423 100644 --- a/src/REPLICA/fix_neb.cpp +++ b/src/REPLICA/fix_neb.cpp @@ -64,6 +64,7 @@ FixNEB::FixNEB(LAMMPS *lmp, int narg, char **arg) : NEBLongRange = false; EqualForceNEB = false; StandardNEB = true; + PerpSpring = FreeEndIni = FreeEndFinal = false; FreeEndFinalWithRespToEIni = FinalAndInterWithRespToEIni = false; kspringPerp = 0.0; @@ -74,6 +75,7 @@ FixNEB::FixNEB(LAMMPS *lmp, int narg, char **arg) : while (iarg < narg) { if (strcmp(arg[iarg],"parallel") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix neb command"); + if (strcmp(arg[iarg+1],"ideal") == 0) { NEBLongRange = true; EqualForceNEB = false; @@ -142,6 +144,7 @@ FixNEB::FixNEB(LAMMPS *lmp, int narg, char **arg) : else procnext = -1; uworld = universe->uworld; + if (NEBLongRange or EqualForceNEB) { int *iroots = new int[nreplica]; MPI_Group uworldgroup,rootgroup; @@ -818,13 +821,16 @@ void FixNEB::inter_replica_comm() } /* -calculate ideal positions +Calculate ideal positions for parallel "ideal" or "equal" */ void FixNEB::calculate_ideal_positions() { + // Skip unless "ideal" or "equal" + if (not (EqualForceNEB or NEBLongRange)) return; + double lentot,lenuntilClimber,meanDistBeforeClimber,meanDistAfterClimber; - if(EqualForceNEB and rclimber>0.) { + if (EqualForceNEB and rclimber>0.) { if (cmode == SINGLE_PROC_DIRECT || cmode == SINGLE_PROC_MAP) { MPI_Allgather(&veng,1,MPI_DOUBLE,&vengall[0],1,MPI_DOUBLE,uworld); } else { @@ -844,6 +850,7 @@ void FixNEB::calculate_ideal_positions() MPI_Bcast(nlenall,nreplica,MPI_DOUBLE,0,world); } } + actualPos = 0.; for (int i = 0; i < ireplica; i++) actualPos += nlenall[i]; @@ -867,6 +874,7 @@ void FixNEB::calculate_ideal_positions() } else idealPos = ireplica * meanDist; idealPos /= meanDist; actualPos /= meanDist; + } /* ---------------------------------------------------------------------- From 1d601f23b102772d1a0d79fe2c131f104108cc10 Mon Sep 17 00:00:00 2001 From: Thomas Swinburne Date: Fri, 9 Dec 2022 16:49:30 +0100 Subject: [PATCH 12/25] terse screen output --- src/REPLICA/fix_neb.cpp | 4 ++-- src/REPLICA/neb.cpp | 41 +++++++++++++++++++++++++++++++---------- src/REPLICA/neb.h | 2 +- 3 files changed, 34 insertions(+), 13 deletions(-) diff --git a/src/REPLICA/fix_neb.cpp b/src/REPLICA/fix_neb.cpp index e512c88423..a3be7531cb 100644 --- a/src/REPLICA/fix_neb.cpp +++ b/src/REPLICA/fix_neb.cpp @@ -64,7 +64,7 @@ FixNEB::FixNEB(LAMMPS *lmp, int narg, char **arg) : NEBLongRange = false; EqualForceNEB = false; StandardNEB = true; - + PerpSpring = FreeEndIni = FreeEndFinal = false; FreeEndFinalWithRespToEIni = FinalAndInterWithRespToEIni = false; kspringPerp = 0.0; @@ -122,7 +122,7 @@ FixNEB::FixNEB(LAMMPS *lmp, int narg, char **arg) : } else error->all(FLERR,"Illegal fix neb command"); iarg += 3; - + } else error->all(FLERR,"Illegal fix neb command"); } diff --git a/src/REPLICA/neb.cpp b/src/REPLICA/neb.cpp index 5ce45475ba..9f67a07fda 100644 --- a/src/REPLICA/neb.cpp +++ b/src/REPLICA/neb.cpp @@ -61,6 +61,7 @@ NEB::NEB(LAMMPS *lmp, double etol_in, double ftol_in, int n1steps_in, n2steps = n2steps_in; nevery = nevery_in; verbose = false; + terse = false; // replica info @@ -147,6 +148,7 @@ void NEB::command(int narg, char **arg) int iarg = 5; int filecmd = 0; verbose=false; + terse=false; while (iarg < narg) { if (strcmp(arg[iarg],"final") == 0) { if (iarg + 2 > narg) error->universe_all(FLERR,"Illegal NEB command: missing arguments"); @@ -166,8 +168,12 @@ void NEB::command(int narg, char **arg) } else if (strcmp(arg[iarg],"verbose") == 0) { verbose=true; ++iarg; + } else if (strcmp(arg[iarg],"terse") == 0) { + terse=true; + ++iarg; } else error->universe_all(FLERR,fmt::format("Unknown NEB command keyword: {}", arg[iarg])); } + if(terse and verbose) terse = false; // verbose overrides terse.... if (!filecmd) error->universe_all(FLERR, "NEB is missing 'final', 'each', or 'none' keyword"); @@ -230,8 +236,11 @@ void NEB::run() if (uscreen) { fmt::print(uscreen," Step {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} ", "MaxReplicaForce", "MaxAtomForce", "GradV0","GradV1","GradVc","EBF", "EBR", "RDT"); - for (int i = 1; i <= nreplica; ++i) - fmt::print(uscreen, "{:^14} {:^14} ", "RD"+std::to_string(i), "PE"+std::to_string(i)); + + if (not terse) { + for (int i = 1; i <= nreplica; ++i) + fmt::print(uscreen, "{:^14} {:^14} ", "RD"+std::to_string(i), "PE"+std::to_string(i)); + } if (verbose) { for (int i = 1; i <= nreplica; ++i) { @@ -247,8 +256,11 @@ void NEB::run() if (ulogfile) { fmt::print(ulogfile," Step {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} ", "MaxReplicaForce", "MaxAtomForce", "GradV0","GradV1","GradVc","EBF", "EBR", "RDT"); - for (int i = 1; i <= nreplica; ++i) - fmt::print(ulogfile, "{:^14} {:^14} ", "RD"+std::to_string(i), "PE"+std::to_string(i)); + + if (not terse) { + for (int i = 1; i <= nreplica; ++i) + fmt::print(ulogfile, "{:^14} {:^14} ", "RD"+std::to_string(i), "PE"+std::to_string(i)); + } if (verbose) { for (int i = 1; i <= nreplica; ++i) { @@ -258,6 +270,7 @@ void NEB::run() "RepForce"+idx, "MaxAtomForce"+idx); } } + fprintf(ulogfile,"\n"); } } @@ -323,9 +336,12 @@ void NEB::run() if (uscreen) { fmt::print(uscreen," Step {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} ", "MaxReplicaForce", "MaxAtomForce", "GradV0","GradV1","GradVc","EBF", "EBR", "RDT"); - for (int i = 1; i <= nreplica; ++i) - fmt::print(uscreen, "{:^14} {:^14} ", "RD"+std::to_string(i), "PE"+std::to_string(i)); - + + if (not terse) { + for (int i = 1; i <= nreplica; ++i) + fmt::print(uscreen, "{:^14} {:^14} ", "RD"+std::to_string(i), "PE"+std::to_string(i)); + } + if (verbose) { for (int i = 1; i <= nreplica; ++i) { auto idx = std::to_string(i); @@ -340,8 +356,11 @@ void NEB::run() if (ulogfile) { fmt::print(ulogfile," Step {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} ", "MaxReplicaForce", "MaxAtomForce", "GradV0","GradV1","GradVc","EBF", "EBR", "RDT"); - for (int i = 1; i <= nreplica; ++i) - fmt::print(ulogfile, "{:^14} {:^14} ", "RD"+std::to_string(i), "PE"+std::to_string(i)); + + if (not terse) { + for (int i = 1; i <= nreplica; ++i) + fmt::print(ulogfile, "{:^14} {:^14} ", "RD"+std::to_string(i), "PE"+std::to_string(i)); + } if (verbose) { for (int i = 1; i <= nreplica; ++i) { @@ -652,7 +671,9 @@ void NEB::print_status() std::string mesg = fmt::format("{:10} {:<14.8g} {:<14.8g} ",update->ntimestep,fmaxreplica,fmaxatom); mesg += fmt::format("{:<14.8g} {:<14.8g} {:<14.8g} ",gradvnorm0,gradvnorm1,gradvnormc); mesg += fmt::format("{:<14.8g} {:<14.8g} {:<14.8g} ",ebf,ebr,endpt); - for (int i = 0; i < nreplica; i++) mesg += fmt::format("{:<14.8g} {:<14.8g} ",rdist[i],all[i][0]); + if (not terse) { + for (int i = 0; i < nreplica; i++) mesg += fmt::format("{:<14.8g} {:<14.8g} ",rdist[i],all[i][0]); + } if (verbose) { mesg += fmt::format("{:<12.5g} {:<12.5g} {:<12.5g} {:<12.5g} {:<12.5g} {:<12.5g}", NAN,180-acos(all[0][5])*todeg,180-acos(all[0][6])*todeg, diff --git a/src/REPLICA/neb.h b/src/REPLICA/neb.h index cad0ac7122..7faf744161 100644 --- a/src/REPLICA/neb.h +++ b/src/REPLICA/neb.h @@ -37,7 +37,7 @@ class NEB : public Command { private: int me, me_universe; // my proc ID in world and universe int ireplica, nreplica; - bool verbose; + bool verbose, terse; MPI_Comm uworld; MPI_Comm roots; // MPI comm with 1 root proc from each world FILE *fp; From 249ac0b34ec2faebb5b9e312c8d49c1997376fe0 Mon Sep 17 00:00:00 2001 From: Thomas Swinburne Date: Fri, 9 Dec 2022 16:52:01 +0100 Subject: [PATCH 13/25] ignore --- .gitignore | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.gitignore b/.gitignore index bd2d0ea705..a89f75cc9a 100644 --- a/.gitignore +++ b/.gitignore @@ -55,3 +55,5 @@ out/RelWithDebInfo out/Release out/x86 out/x64 +src/Makefile.package-e +src/Makefile.package.settings-e From 6113bd5aa49585d6d8af1d7beaaa4b34801d9494 Mon Sep 17 00:00:00 2001 From: Thomas Swinburne Date: Mon, 12 Dec 2022 11:30:21 +0100 Subject: [PATCH 14/25] small clean up --- src/REPLICA/fix_neb.cpp | 21 +++++++++++---------- src/REPLICA/fix_neb.h | 2 +- 2 files changed, 12 insertions(+), 11 deletions(-) diff --git a/src/REPLICA/fix_neb.cpp b/src/REPLICA/fix_neb.cpp index a3be7531cb..109f061901 100644 --- a/src/REPLICA/fix_neb.cpp +++ b/src/REPLICA/fix_neb.cpp @@ -70,12 +70,10 @@ FixNEB::FixNEB(LAMMPS *lmp, int narg, char **arg) : kspringPerp = 0.0; kspringIni = 1.0; kspringFinal = 1.0; - int iarg = 4; while (iarg < narg) { if (strcmp(arg[iarg],"parallel") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix neb command"); - if (strcmp(arg[iarg+1],"ideal") == 0) { NEBLongRange = true; EqualForceNEB = false; @@ -828,7 +826,8 @@ void FixNEB::calculate_ideal_positions() // Skip unless "ideal" or "equal" if (not (EqualForceNEB or NEBLongRange)) return; - double lentot,lenuntilClimber,meanDistBeforeClimber,meanDistAfterClimber; + double Elentot,lentot,lenuntilClimber; + double meanDist,meanDistBeforeClimber,meanDistAfterClimber; if (EqualForceNEB and rclimber>0.) { if (cmode == SINGLE_PROC_DIRECT || cmode == SINGLE_PROC_MAP) { @@ -838,9 +837,11 @@ void FixNEB::calculate_ideal_positions() MPI_Allgather(&veng,1,MPI_DOUBLE,&vengall[0],1,MPI_DOUBLE,rootworld); MPI_Bcast(vengall,nreplica,MPI_DOUBLE,0,world); } + for (int i = 0; i < nreplica-1; i++) nlenall[i] = std::abs(vengall[i+1]-vengall[i]); nlenall[nreplica-1] = 0.0; + } else if(NEBLongRange or EqualForceNEB) { if (cmode == SINGLE_PROC_DIRECT || cmode == SINGLE_PROC_MAP) { MPI_Allgather(&nlen,1,MPI_DOUBLE,&nlenall[0],1,MPI_DOUBLE,uworld); @@ -848,17 +849,17 @@ void FixNEB::calculate_ideal_positions() if (me == 0) MPI_Allgather(&nlen,1,MPI_DOUBLE,&nlenall[0],1,MPI_DOUBLE,rootworld); MPI_Bcast(nlenall,nreplica,MPI_DOUBLE,0,world); - } + } } - - actualPos = 0.; - for (int i = 0; i < ireplica; i++) - actualPos += nlenall[i]; - + lentot = 0.; for (int i = 0; i < nreplica; i++) lentot += nlenall[i]; - + + actualPos = 0.; + for (int i = 0; i < ireplica; i++) + actualPos += nlenall[i]; + meanDist = lentot/(nreplica-1); lenuntilClimber = 0.; diff --git a/src/REPLICA/fix_neb.h b/src/REPLICA/fix_neb.h index 442bf31199..ad79d62c2b 100644 --- a/src/REPLICA/fix_neb.h +++ b/src/REPLICA/fix_neb.h @@ -38,7 +38,7 @@ class FixNEB : public Fix { private: int me, nprocs, nprocs_universe; - double kspring, kspringIni, kspringFinal, kspringPerp, EIniIni, EFinalIni, idealPos, actualPos, meanDist; + double kspring, kspringIni, kspringFinal, kspringPerp, EIniIni, EFinalIni, idealPos, actualPos; bool StandardNEB, NEBLongRange, EqualForceNEB, PerpSpring, FreeEndIni, FreeEndFinal; bool FreeEndFinalWithRespToEIni, FinalAndInterWithRespToEIni; int ireplica, nreplica; From b7db402c2d84ff40661e32993c73a30babd27ebf Mon Sep 17 00:00:00 2001 From: Thomas Swinburne Date: Thu, 12 Jan 2023 18:30:15 +0100 Subject: [PATCH 15/25] post-axel updates --- doc/src/neb.rst | 7 ++++++- src/REPLICA/fix_neb.cpp | 37 ++++++++++++++----------------------- src/REPLICA/fix_neb.h | 5 +++-- src/REPLICA/neb.cpp | 38 +++++++++++++++++--------------------- src/REPLICA/neb.h | 3 ++- 5 files changed, 42 insertions(+), 48 deletions(-) diff --git a/doc/src/neb.rst b/doc/src/neb.rst index 69a0877f3a..1ba1bd1cc2 100644 --- a/doc/src/neb.rst +++ b/doc/src/neb.rst @@ -29,7 +29,7 @@ Syntax *none* arg = no argument all replicas assumed to already have their initial coords -keyword = *verbose* +keyword = *verbose* or *terse* Examples """""""" @@ -347,6 +347,11 @@ energy gradient of image i. ReplicaForce is the two-norm of the 3N-length force vector (including nudging forces) for replica i. MaxAtomForce is the maximum force component of any atom in replica i. +Alternatively, a restricted print out can be obtained by adding the +terse keyword, which omits per-replica information. This typically +fits on one line of a command terminal, aiding visual inspection of +an ongoing NEB calculation. + When a NEB calculation does not converge properly, the supplementary information can help understanding what is going wrong. For instance when the path angle becomes acute, the definition of tangent used in diff --git a/src/REPLICA/fix_neb.cpp b/src/REPLICA/fix_neb.cpp index 109f061901..ba49ee65b3 100644 --- a/src/REPLICA/fix_neb.cpp +++ b/src/REPLICA/fix_neb.cpp @@ -61,10 +61,7 @@ FixNEB::FixNEB(LAMMPS *lmp, int narg, char **arg) : // optional params - NEBLongRange = false; - EqualForceNEB = false; - StandardNEB = true; - + neb_mode = 0; PerpSpring = FreeEndIni = FreeEndFinal = false; FreeEndFinalWithRespToEIni = FinalAndInterWithRespToEIni = false; kspringPerp = 0.0; @@ -75,17 +72,11 @@ FixNEB::FixNEB(LAMMPS *lmp, int narg, char **arg) : if (strcmp(arg[iarg],"parallel") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix neb command"); if (strcmp(arg[iarg+1],"ideal") == 0) { - NEBLongRange = true; - EqualForceNEB = false; - StandardNEB = false; + neb_mode = IDEAL; } else if (strcmp(arg[iarg+1],"equal") == 0) { - NEBLongRange = false; - EqualForceNEB = true; - StandardNEB = false; + neb_mode = EQUAL; } else if (strcmp(arg[iarg+1],"neigh") == 0) { - NEBLongRange = false; - EqualForceNEB = false; - StandardNEB = true; + neb_mode = NEIGHBOR; } else error->all(FLERR,"Illegal fix neb command"); iarg += 2; @@ -143,7 +134,7 @@ FixNEB::FixNEB(LAMMPS *lmp, int narg, char **arg) : uworld = universe->uworld; - if (NEBLongRange or EqualForceNEB) { + if (neb_mode==IDEAL || neb_mode==EQUAL) { int *iroots = new int[nreplica]; MPI_Group uworldgroup,rootgroup; @@ -197,11 +188,11 @@ FixNEB::~FixNEB() memory->destroy(counts); memory->destroy(displacements); - if (NEBLongRange or EqualForceNEB) { + if (neb_mode==IDEAL or neb_mode==EQUAL) { if (rootworld != MPI_COMM_NULL) MPI_Comm_free(&rootworld); memory->destroy(nlenall); } - if(EqualForceNEB) memory->destroy(vengall); + if(neb_mode==EQUAL) memory->destroy(vengall); } /* ---------------------------------------------------------------------- */ @@ -574,9 +565,9 @@ void FixNEB::min_post_force(int /*vflag*/) if (ireplica == rclimber) prefactor = -2.0*dot; else { - if (NEBLongRange or EqualForceNEB) { + if (neb_mode==IDEAL||neb_mode==EQUAL) { prefactor = -dot - kspring*(actualPos-idealPos)/2; - } else if (StandardNEB) { + } else { prefactor = -dot + kspring*(nlen-plen); } @@ -824,12 +815,12 @@ Calculate ideal positions for parallel "ideal" or "equal" void FixNEB::calculate_ideal_positions() { // Skip unless "ideal" or "equal" - if (not (EqualForceNEB or NEBLongRange)) return; + if (!(neb_mode==IDEAL||neb_mode==EQUAL)) return; double Elentot,lentot,lenuntilClimber; double meanDist,meanDistBeforeClimber,meanDistAfterClimber; - if (EqualForceNEB and rclimber>0.) { + if (neb_mode==EQUAL and rclimber>0.) { if (cmode == SINGLE_PROC_DIRECT || cmode == SINGLE_PROC_MAP) { MPI_Allgather(&veng,1,MPI_DOUBLE,&vengall[0],1,MPI_DOUBLE,uworld); } else { @@ -842,7 +833,7 @@ void FixNEB::calculate_ideal_positions() nlenall[i] = std::abs(vengall[i+1]-vengall[i]); nlenall[nreplica-1] = 0.0; - } else if(NEBLongRange or EqualForceNEB) { + } else if(neb_mode==IDEAL || neb_mode==EQUAL) { if (cmode == SINGLE_PROC_DIRECT || cmode == SINGLE_PROC_MAP) { MPI_Allgather(&nlen,1,MPI_DOUBLE,&nlenall[0],1,MPI_DOUBLE,uworld); } else { @@ -914,11 +905,11 @@ void FixNEB::reallocate() memory->create(tagrecv,maxlocal,"neb:tagrecv"); } - if (NEBLongRange or EqualForceNEB) { + if (neb_mode==IDEAL or neb_mode==EQUAL) { memory->destroy(nlenall); memory->create(nlenall,nreplica,"neb:nlenall"); } - if (EqualForceNEB) { + if (neb_mode==EQUAL) { memory->destroy(vengall); memory->create(vengall,nreplica,"neb:vengall"); } diff --git a/src/REPLICA/fix_neb.h b/src/REPLICA/fix_neb.h index ad79d62c2b..b35e998c2c 100644 --- a/src/REPLICA/fix_neb.h +++ b/src/REPLICA/fix_neb.h @@ -37,9 +37,10 @@ class FixNEB : public Fix { void min_post_force(int) override; private: - int me, nprocs, nprocs_universe; + enum {NEIGHBOR=0, IDEAL=1, EQUAL=2}; + int me, nprocs, nprocs_universe, neb_mode; double kspring, kspringIni, kspringFinal, kspringPerp, EIniIni, EFinalIni, idealPos, actualPos; - bool StandardNEB, NEBLongRange, EqualForceNEB, PerpSpring, FreeEndIni, FreeEndFinal; + bool PerpSpring, FreeEndIni, FreeEndFinal; bool FreeEndFinalWithRespToEIni, FinalAndInterWithRespToEIni; int ireplica, nreplica; int procnext, procprev; diff --git a/src/REPLICA/neb.cpp b/src/REPLICA/neb.cpp index 9f67a07fda..c10d882de0 100644 --- a/src/REPLICA/neb.cpp +++ b/src/REPLICA/neb.cpp @@ -60,8 +60,6 @@ NEB::NEB(LAMMPS *lmp, double etol_in, double ftol_in, int n1steps_in, n1steps = n1steps_in; n2steps = n2steps_in; nevery = nevery_in; - verbose = false; - terse = false; // replica info @@ -147,8 +145,7 @@ void NEB::command(int narg, char **arg) // process file-style setting to setup initial configs for all replicas int iarg = 5; int filecmd = 0; - verbose=false; - terse=false; + print_mode = 0; // normal while (iarg < narg) { if (strcmp(arg[iarg],"final") == 0) { if (iarg + 2 > narg) error->universe_all(FLERR,"Illegal NEB command: missing arguments"); @@ -166,14 +163,13 @@ void NEB::command(int narg, char **arg) filecmd = 1; ++iarg; } else if (strcmp(arg[iarg],"verbose") == 0) { - verbose=true; + print_mode = VERBOSE; ++iarg; } else if (strcmp(arg[iarg],"terse") == 0) { - terse=true; + print_mode = TERSE; ++iarg; } else error->universe_all(FLERR,fmt::format("Unknown NEB command keyword: {}", arg[iarg])); } - if(terse and verbose) terse = false; // verbose overrides terse.... if (!filecmd) error->universe_all(FLERR, "NEB is missing 'final', 'each', or 'none' keyword"); @@ -200,7 +196,7 @@ void NEB::run() error->universe_all(FLERR,"NEB requires use of exactly one fix neb instance"); fneb = dynamic_cast(fixes[0]); - if (verbose) numall =7; + if (print_mode==VERBOSE) numall =7; else numall = 4; memory->create(all,nreplica,numall,"neb:all"); rdist = new double[nreplica]; @@ -237,12 +233,12 @@ void NEB::run() fmt::print(uscreen," Step {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} ", "MaxReplicaForce", "MaxAtomForce", "GradV0","GradV1","GradVc","EBF", "EBR", "RDT"); - if (not terse) { + if (print_mode!=TERSE) { for (int i = 1; i <= nreplica; ++i) fmt::print(uscreen, "{:^14} {:^14} ", "RD"+std::to_string(i), "PE"+std::to_string(i)); } - if (verbose) { + if (print_mode==VERBOSE) { for (int i = 1; i <= nreplica; ++i) { auto idx = std::to_string(i); fmt::print(uscreen, "{:^12}{:^12}{:^12} {:^12} {:^12}{:^12} ", @@ -257,12 +253,12 @@ void NEB::run() fmt::print(ulogfile," Step {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} ", "MaxReplicaForce", "MaxAtomForce", "GradV0","GradV1","GradVc","EBF", "EBR", "RDT"); - if (not terse) { + if (print_mode!=TERSE) { for (int i = 1; i <= nreplica; ++i) fmt::print(ulogfile, "{:^14} {:^14} ", "RD"+std::to_string(i), "PE"+std::to_string(i)); } - if (verbose) { + if (print_mode==VERBOSE) { for (int i = 1; i <= nreplica; ++i) { auto idx = std::to_string(i); fmt::print(ulogfile, "{:^12}{:^12}{:^12} {:^12} {:^12}{:^12} ", @@ -337,12 +333,12 @@ void NEB::run() fmt::print(uscreen," Step {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} ", "MaxReplicaForce", "MaxAtomForce", "GradV0","GradV1","GradVc","EBF", "EBR", "RDT"); - if (not terse) { + if (print_mode!=TERSE) { for (int i = 1; i <= nreplica; ++i) fmt::print(uscreen, "{:^14} {:^14} ", "RD"+std::to_string(i), "PE"+std::to_string(i)); } - if (verbose) { + if (print_mode==VERBOSE) { for (int i = 1; i <= nreplica; ++i) { auto idx = std::to_string(i); fmt::print(uscreen, "{:^12}{:^12}{:^12} {:^12} {:^12}{:^12} ", @@ -357,12 +353,12 @@ void NEB::run() fmt::print(ulogfile," Step {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} ", "MaxReplicaForce", "MaxAtomForce", "GradV0","GradV1","GradVc","EBF", "EBR", "RDT"); - if (not terse) { + if (print_mode!=TERSE) { for (int i = 1; i <= nreplica; ++i) fmt::print(ulogfile, "{:^14} {:^14} ", "RD"+std::to_string(i), "PE"+std::to_string(i)); } - if (verbose) { + if (print_mode==VERBOSE) { for (int i = 1; i <= nreplica; ++i) { auto idx = std::to_string(i); fmt::print(ulogfile, "{:^12}{:^12}{:^12} {:^12} {:^12}{:^12} ", @@ -605,7 +601,7 @@ void NEB::print_status() double fmaxatom; MPI_Allreduce(&fnorminf,&fmaxatom,1,MPI_DOUBLE,MPI_MAX,roots); - if (verbose) { + if (print_mode==VERBOSE) { freplica = new double[nreplica]; MPI_Allgather(&fnorm2,1,MPI_DOUBLE,&freplica[0],1,MPI_DOUBLE,roots); fmaxatomInRepl = new double[nreplica]; @@ -618,7 +614,7 @@ void NEB::print_status() one[2] = fneb->nlen; one[3] = fneb->gradlen; - if (verbose) { + if (print_mode==VERBOSE) { one[4] = fneb->dotpath; one[5] = fneb->dottangrad; one[6] = fneb->dotgrad; @@ -671,10 +667,10 @@ void NEB::print_status() std::string mesg = fmt::format("{:10} {:<14.8g} {:<14.8g} ",update->ntimestep,fmaxreplica,fmaxatom); mesg += fmt::format("{:<14.8g} {:<14.8g} {:<14.8g} ",gradvnorm0,gradvnorm1,gradvnormc); mesg += fmt::format("{:<14.8g} {:<14.8g} {:<14.8g} ",ebf,ebr,endpt); - if (not terse) { + if (print_mode!=TERSE) { for (int i = 0; i < nreplica; i++) mesg += fmt::format("{:<14.8g} {:<14.8g} ",rdist[i],all[i][0]); } - if (verbose) { + if (print_mode==VERBOSE) { mesg += fmt::format("{:<12.5g} {:<12.5g} {:<12.5g} {:<12.5g} {:<12.5g} {:<12.5g}", NAN,180-acos(all[0][5])*todeg,180-acos(all[0][6])*todeg, all[0][3],freplica[0],fmaxatomInRepl[0]); @@ -694,7 +690,7 @@ void NEB::print_status() fflush(universe->ulogfile); } } - if (verbose) { + if (print_mode==VERBOSE) { delete[] freplica; delete[] fmaxatomInRepl; } diff --git a/src/REPLICA/neb.h b/src/REPLICA/neb.h index 7faf744161..fe114684aa 100644 --- a/src/REPLICA/neb.h +++ b/src/REPLICA/neb.h @@ -35,9 +35,10 @@ class NEB : public Command { double ebf, ebr; // forward and reverse energy barriers private: + enum {NORMAL=0,TERSE=1,VERBOSE=2}; + int print_mode; // output verbosity int me, me_universe; // my proc ID in world and universe int ireplica, nreplica; - bool verbose, terse; MPI_Comm uworld; MPI_Comm roots; // MPI comm with 1 root proc from each world FILE *fp; From bdf6cdd32727a62f0ba0a84ce59ceb37c6b2b955 Mon Sep 17 00:00:00 2001 From: Thomas Swinburne Date: Thu, 12 Jan 2023 18:34:54 +0100 Subject: [PATCH 16/25] found two or -> || --- src/REPLICA/fix_neb.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/REPLICA/fix_neb.cpp b/src/REPLICA/fix_neb.cpp index ba49ee65b3..0ef4c2f6e2 100644 --- a/src/REPLICA/fix_neb.cpp +++ b/src/REPLICA/fix_neb.cpp @@ -188,7 +188,7 @@ FixNEB::~FixNEB() memory->destroy(counts); memory->destroy(displacements); - if (neb_mode==IDEAL or neb_mode==EQUAL) { + if (neb_mode==IDEAL || neb_mode==EQUAL) { if (rootworld != MPI_COMM_NULL) MPI_Comm_free(&rootworld); memory->destroy(nlenall); } @@ -905,7 +905,7 @@ void FixNEB::reallocate() memory->create(tagrecv,maxlocal,"neb:tagrecv"); } - if (neb_mode==IDEAL or neb_mode==EQUAL) { + if (neb_mode==IDEAL || neb_mode==EQUAL) { memory->destroy(nlenall); memory->create(nlenall,nreplica,"neb:nlenall"); } From c8f380ffbb61fee25efbbf74018dcdce252774a8 Mon Sep 17 00:00:00 2001 From: Thomas Swinburne Date: Thu, 12 Jan 2023 21:24:40 +0100 Subject: [PATCH 17/25] small changes- still not compiling on Windows... --- src/REPLICA/fix_neb.cpp | 2 +- src/REPLICA/neb.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/REPLICA/fix_neb.cpp b/src/REPLICA/fix_neb.cpp index 0ef4c2f6e2..c9a6d05442 100644 --- a/src/REPLICA/fix_neb.cpp +++ b/src/REPLICA/fix_neb.cpp @@ -61,7 +61,7 @@ FixNEB::FixNEB(LAMMPS *lmp, int narg, char **arg) : // optional params - neb_mode = 0; + neb_mode = NEIGHBOR; // default setting PerpSpring = FreeEndIni = FreeEndFinal = false; FreeEndFinalWithRespToEIni = FinalAndInterWithRespToEIni = false; kspringPerp = 0.0; diff --git a/src/REPLICA/neb.cpp b/src/REPLICA/neb.cpp index c10d882de0..dd7cc38247 100644 --- a/src/REPLICA/neb.cpp +++ b/src/REPLICA/neb.cpp @@ -145,7 +145,7 @@ void NEB::command(int narg, char **arg) // process file-style setting to setup initial configs for all replicas int iarg = 5; int filecmd = 0; - print_mode = 0; // normal + print_mode = NORMAL; // normal while (iarg < narg) { if (strcmp(arg[iarg],"final") == 0) { if (iarg + 2 > narg) error->universe_all(FLERR,"Illegal NEB command: missing arguments"); From a6234ab3be336873f9a25b9b26c3278b41d572c3 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 12 Jan 2023 18:50:09 -0500 Subject: [PATCH 18/25] move enum to .cpp file and away from header --- src/REPLICA/fix_neb.cpp | 3 ++- src/REPLICA/fix_neb.h | 1 - src/REPLICA/neb.cpp | 1 + src/REPLICA/neb.h | 3 +-- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/REPLICA/fix_neb.cpp b/src/REPLICA/fix_neb.cpp index c9a6d05442..68e60ca334 100644 --- a/src/REPLICA/fix_neb.cpp +++ b/src/REPLICA/fix_neb.cpp @@ -38,7 +38,8 @@ using namespace LAMMPS_NS; using namespace FixConst; using namespace MathConst; -enum{SINGLE_PROC_DIRECT,SINGLE_PROC_MAP,MULTI_PROC}; +enum { SINGLE_PROC_DIRECT, SINGLE_PROC_MAP, MULTI_PROC }; +enum { NEIGHBOR, IDEAL, EQUAL }; #define BUFSIZE 8 diff --git a/src/REPLICA/fix_neb.h b/src/REPLICA/fix_neb.h index b35e998c2c..036d90b828 100644 --- a/src/REPLICA/fix_neb.h +++ b/src/REPLICA/fix_neb.h @@ -37,7 +37,6 @@ class FixNEB : public Fix { void min_post_force(int) override; private: - enum {NEIGHBOR=0, IDEAL=1, EQUAL=2}; int me, nprocs, nprocs_universe, neb_mode; double kspring, kspringIni, kspringFinal, kspringPerp, EIniIni, EFinalIni, idealPos, actualPos; bool PerpSpring, FreeEndIni, FreeEndFinal; diff --git a/src/REPLICA/neb.cpp b/src/REPLICA/neb.cpp index dd7cc38247..b65a08d116 100644 --- a/src/REPLICA/neb.cpp +++ b/src/REPLICA/neb.cpp @@ -41,6 +41,7 @@ using namespace MathConst; #define CHUNK 1024 #define ATTRIBUTE_PERLINE 4 +enum { NORMAL, TERSE, VERBOSE }; /* ---------------------------------------------------------------------- */ NEB::NEB(LAMMPS *lmp) : Command(lmp), all(nullptr), rdist(nullptr) {} diff --git a/src/REPLICA/neb.h b/src/REPLICA/neb.h index fe114684aa..db5388ac20 100644 --- a/src/REPLICA/neb.h +++ b/src/REPLICA/neb.h @@ -35,8 +35,7 @@ class NEB : public Command { double ebf, ebr; // forward and reverse energy barriers private: - enum {NORMAL=0,TERSE=1,VERBOSE=2}; - int print_mode; // output verbosity + int print_mode; // output verbosity int me, me_universe; // my proc ID in world and universe int ireplica, nreplica; MPI_Comm uworld; From 7ac611b671ca85aa56e9d13d609d9c969aee8af8 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 12 Jan 2023 18:50:58 -0500 Subject: [PATCH 19/25] enable and apply clang-format --- src/REPLICA/fix_neb.cpp | 628 +++++++++++++++++++--------------------- src/REPLICA/neb.cpp | 381 ++++++++++++------------ 2 files changed, 500 insertions(+), 509 deletions(-) diff --git a/src/REPLICA/fix_neb.cpp b/src/REPLICA/fix_neb.cpp index 68e60ca334..1ef3019b92 100644 --- a/src/REPLICA/fix_neb.cpp +++ b/src/REPLICA/fix_neb.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -46,23 +45,21 @@ enum { NEIGHBOR, IDEAL, EQUAL }; /* ---------------------------------------------------------------------- */ FixNEB::FixNEB(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg), - id_pe(nullptr), pe(nullptr), nlenall(nullptr), vengall(nullptr), xprev(nullptr), xnext(nullptr), - fnext(nullptr), springF(nullptr), tangent(nullptr), xsend(nullptr), xrecv(nullptr), - fsend(nullptr), frecv(nullptr), tagsend(nullptr), tagrecv(nullptr), - xsendall(nullptr), xrecvall(nullptr), fsendall(nullptr), frecvall(nullptr), - tagsendall(nullptr), tagrecvall(nullptr), counts(nullptr), - displacements(nullptr) + Fix(lmp, narg, arg), id_pe(nullptr), pe(nullptr), nlenall(nullptr), vengall(nullptr), + xprev(nullptr), xnext(nullptr), fnext(nullptr), springF(nullptr), tangent(nullptr), + xsend(nullptr), xrecv(nullptr), fsend(nullptr), frecv(nullptr), tagsend(nullptr), + tagrecv(nullptr), xsendall(nullptr), xrecvall(nullptr), fsendall(nullptr), frecvall(nullptr), + tagsendall(nullptr), tagrecvall(nullptr), counts(nullptr), displacements(nullptr) { - if (narg < 4) error->all(FLERR,"Illegal fix neb command"); + if (narg < 4) error->all(FLERR, "Illegal fix neb command"); - kspring = utils::numeric(FLERR,arg[3],false,lmp); - if (kspring <= 0.0) error->all(FLERR,"Illegal fix neb command"); + kspring = utils::numeric(FLERR, arg[3], false, lmp); + if (kspring <= 0.0) error->all(FLERR, "Illegal fix neb command"); // optional params - neb_mode = NEIGHBOR; // default setting + neb_mode = NEIGHBOR; // default setting PerpSpring = FreeEndIni = FreeEndFinal = false; FreeEndFinalWithRespToEIni = FinalAndInterWithRespToEIni = false; kspringPerp = 0.0; @@ -70,50 +67,53 @@ FixNEB::FixNEB(LAMMPS *lmp, int narg, char **arg) : kspringFinal = 1.0; int iarg = 4; while (iarg < narg) { - if (strcmp(arg[iarg],"parallel") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix neb command"); - if (strcmp(arg[iarg+1],"ideal") == 0) { + if (strcmp(arg[iarg], "parallel") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal fix neb command"); + if (strcmp(arg[iarg + 1], "ideal") == 0) { neb_mode = IDEAL; - } else if (strcmp(arg[iarg+1],"equal") == 0) { + } else if (strcmp(arg[iarg + 1], "equal") == 0) { neb_mode = EQUAL; - } else if (strcmp(arg[iarg+1],"neigh") == 0) { + } else if (strcmp(arg[iarg + 1], "neigh") == 0) { neb_mode = NEIGHBOR; - } else error->all(FLERR,"Illegal fix neb command"); + } else + error->all(FLERR, "Illegal fix neb command"); iarg += 2; - } else if (strcmp(arg[iarg],"perp") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix neb command"); + } else if (strcmp(arg[iarg], "perp") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal fix neb command"); PerpSpring = true; - kspringPerp = utils::numeric(FLERR,arg[iarg+1],false,lmp); + kspringPerp = utils::numeric(FLERR, arg[iarg + 1], false, lmp); if (kspringPerp == 0.0) PerpSpring = false; - if (kspringPerp < 0.0) error->all(FLERR,"Illegal fix neb command"); + if (kspringPerp < 0.0) error->all(FLERR, "Illegal fix neb command"); iarg += 2; - } else if (strcmp (arg[iarg],"end") == 0) { - if (iarg+3 > narg) error->all(FLERR,"Illegal fix neb command"); - if (strcmp(arg[iarg+1],"first") == 0) { + } else if (strcmp(arg[iarg], "end") == 0) { + if (iarg + 3 > narg) error->all(FLERR, "Illegal fix neb command"); + if (strcmp(arg[iarg + 1], "first") == 0) { FreeEndIni = true; - kspringIni = utils::numeric(FLERR,arg[iarg+2],false,lmp); - } else if (strcmp(arg[iarg+1],"last") == 0) { + kspringIni = utils::numeric(FLERR, arg[iarg + 2], false, lmp); + } else if (strcmp(arg[iarg + 1], "last") == 0) { FreeEndFinal = true; FinalAndInterWithRespToEIni = false; FreeEndFinalWithRespToEIni = false; - kspringFinal = utils::numeric(FLERR,arg[iarg+2],false,lmp); - } else if (strcmp(arg[iarg+1],"last/efirst") == 0) { + kspringFinal = utils::numeric(FLERR, arg[iarg + 2], false, lmp); + } else if (strcmp(arg[iarg + 1], "last/efirst") == 0) { FreeEndFinal = false; FinalAndInterWithRespToEIni = false; FreeEndFinalWithRespToEIni = true; - kspringFinal = utils::numeric(FLERR,arg[iarg+2],false,lmp); - } else if (strcmp(arg[iarg+1],"last/efirst/middle") == 0) { + kspringFinal = utils::numeric(FLERR, arg[iarg + 2], false, lmp); + } else if (strcmp(arg[iarg + 1], "last/efirst/middle") == 0) { FreeEndFinal = false; FinalAndInterWithRespToEIni = true; FreeEndFinalWithRespToEIni = true; - kspringFinal = utils::numeric(FLERR,arg[iarg+2],false,lmp); - } else error->all(FLERR,"Illegal fix neb command"); + kspringFinal = utils::numeric(FLERR, arg[iarg + 2], false, lmp); + } else + error->all(FLERR, "Illegal fix neb command"); iarg += 3; - - } else error->all(FLERR,"Illegal fix neb command"); + + } else + error->all(FLERR, "Illegal fix neb command"); } // nreplica = number of partitions @@ -128,31 +128,35 @@ FixNEB::FixNEB(LAMMPS *lmp, int narg, char **arg) : nreplica = universe->nworlds; ireplica = universe->iworld; - if (ireplica > 0) procprev = universe->root_proc[ireplica-1]; - else procprev = -1; - if (ireplica < nreplica-1) procnext = universe->root_proc[ireplica+1]; - else procnext = -1; + if (ireplica > 0) + procprev = universe->root_proc[ireplica - 1]; + else + procprev = -1; + if (ireplica < nreplica - 1) + procnext = universe->root_proc[ireplica + 1]; + else + procnext = -1; uworld = universe->uworld; - if (neb_mode==IDEAL || neb_mode==EQUAL) { + if ((neb_mode == IDEAL) || (neb_mode == EQUAL)) { int *iroots = new int[nreplica]; - MPI_Group uworldgroup,rootgroup; + MPI_Group uworldgroup, rootgroup; - for (int i=0; iroot_proc[i]; + for (int i = 0; i < nreplica; i++) iroots[i] = universe->root_proc[i]; MPI_Comm_group(uworld, &uworldgroup); MPI_Group_incl(uworldgroup, nreplica, iroots, &rootgroup); MPI_Comm_create(uworld, rootgroup, &rootworld); if (rootgroup != MPI_GROUP_NULL) MPI_Group_free(&rootgroup); if (uworldgroup != MPI_GROUP_NULL) MPI_Group_free(&uworldgroup); - delete [] iroots; + delete[] iroots; } // create a new compute pe style // id = fix-ID + pe, compute group = all - id_pe = utils::strdup(std::string(id)+"_pe"); - modify->add_compute(std::string(id_pe)+" all pe"); + id_pe = utils::strdup(std::string(id) + "_pe"); + modify->add_compute(std::string(id_pe) + " all pe"); // initialize local storage @@ -165,7 +169,7 @@ FixNEB::FixNEB(LAMMPS *lmp, int narg, char **arg) : FixNEB::~FixNEB() { modify->delete_compute(id_pe); - delete [] id_pe; + delete[] id_pe; memory->destroy(xprev); memory->destroy(xnext); @@ -189,11 +193,11 @@ FixNEB::~FixNEB() memory->destroy(counts); memory->destroy(displacements); - if (neb_mode==IDEAL || neb_mode==EQUAL) { + if ((neb_mode == IDEAL) || (neb_mode == EQUAL)) { if (rootworld != MPI_COMM_NULL) MPI_Comm_free(&rootworld); memory->destroy(nlenall); } - if(neb_mode==EQUAL) memory->destroy(vengall); + if (neb_mode == EQUAL) memory->destroy(vengall); } /* ---------------------------------------------------------------------- */ @@ -210,8 +214,7 @@ int FixNEB::setmask() void FixNEB::init() { int icompute = modify->find_compute(id_pe); - if (icompute < 0) - error->all(FLERR,"Potential energy ID for fix neb does not exist"); + if (icompute < 0) error->all(FLERR, "Potential energy ID for fix neb does not exist"); pe = modify->compute[icompute]; // turn off climbing mode, NEB command turns it on after init() @@ -221,33 +224,34 @@ void FixNEB::init() // nebatoms = # of atoms in fix group = atoms with inter-replica forces bigint count = group->count(igroup); - if (count > MAXSMALLINT) error->all(FLERR,"Too many active NEB atoms"); + if (count > MAXSMALLINT) error->all(FLERR, "Too many active NEB atoms"); nebatoms = count; // comm mode for inter-replica exchange of coords - if (nreplica == nprocs_universe && - nebatoms == atom->natoms && atom->sortfreq == 0) + if (nreplica == nprocs_universe && nebatoms == atom->natoms && atom->sortfreq == 0) cmode = SINGLE_PROC_DIRECT; - else if (nreplica == nprocs_universe) cmode = SINGLE_PROC_MAP; - else cmode = MULTI_PROC; + else if (nreplica == nprocs_universe) + cmode = SINGLE_PROC_MAP; + else + cmode = MULTI_PROC; // ntotal = total # of atoms in system, NEB atoms or not - if (atom->natoms > MAXSMALLINT) error->all(FLERR,"Too many atoms for NEB"); + if (atom->natoms > MAXSMALLINT) error->all(FLERR, "Too many atoms for NEB"); ntotal = atom->natoms; if (atom->nmax > maxlocal) reallocate(); if ((cmode == MULTI_PROC) && (counts == nullptr)) { - memory->create(xsendall,ntotal,3,"neb:xsendall"); - memory->create(xrecvall,ntotal,3,"neb:xrecvall"); - memory->create(fsendall,ntotal,3,"neb:fsendall"); - memory->create(frecvall,ntotal,3,"neb:frecvall"); - memory->create(tagsendall,ntotal,"neb:tagsendall"); - memory->create(tagrecvall,ntotal,"neb:tagrecvall"); - memory->create(counts,nprocs,"neb:counts"); - memory->create(displacements,nprocs,"neb:displacements"); + memory->create(xsendall, ntotal, 3, "neb:xsendall"); + memory->create(xrecvall, ntotal, 3, "neb:xrecvall"); + memory->create(fsendall, ntotal, 3, "neb:fsendall"); + memory->create(frecvall, ntotal, 3, "neb:frecvall"); + memory->create(tagsendall, ntotal, "neb:tagsendall"); + memory->create(tagrecvall, ntotal, "neb:tagrecvall"); + memory->create(counts, nprocs, "neb:counts"); + memory->create(displacements, nprocs, "neb:displacements"); } } @@ -259,60 +263,57 @@ void FixNEB::min_setup(int vflag) // trigger potential energy computation on next timestep - pe->addstep(update->ntimestep+1); + pe->addstep(update->ntimestep + 1); } /* ---------------------------------------------------------------------- */ void FixNEB::min_post_force(int /*vflag*/) { - double vprev,vnext; - double delxp,delyp,delzp,delxn,delyn,delzn; - double vIni=0.0; + double vprev, vnext; + double delxp, delyp, delzp, delxn, delyn, delzn; + double vIni = 0.0; vprev = vnext = veng = pe->compute_scalar(); - if (ireplica < nreplica-1 && me == 0) - MPI_Send(&veng,1,MPI_DOUBLE,procnext,0,uworld); + if (ireplica < nreplica - 1 && me == 0) MPI_Send(&veng, 1, MPI_DOUBLE, procnext, 0, uworld); if (ireplica > 0 && me == 0) - MPI_Recv(&vprev,1,MPI_DOUBLE,procprev,0,uworld,MPI_STATUS_IGNORE); + MPI_Recv(&vprev, 1, MPI_DOUBLE, procprev, 0, uworld, MPI_STATUS_IGNORE); - if (ireplica > 0 && me == 0) - MPI_Send(&veng,1,MPI_DOUBLE,procprev,0,uworld); - if (ireplica < nreplica-1 && me == 0) - MPI_Recv(&vnext,1,MPI_DOUBLE,procnext,0,uworld,MPI_STATUS_IGNORE); + if (ireplica > 0 && me == 0) MPI_Send(&veng, 1, MPI_DOUBLE, procprev, 0, uworld); + if (ireplica < nreplica - 1 && me == 0) + MPI_Recv(&vnext, 1, MPI_DOUBLE, procnext, 0, uworld, MPI_STATUS_IGNORE); if (cmode == MULTI_PROC) { - MPI_Bcast(&vprev,1,MPI_DOUBLE,0,world); - MPI_Bcast(&vnext,1,MPI_DOUBLE,0,world); + MPI_Bcast(&vprev, 1, MPI_DOUBLE, 0, world); + MPI_Bcast(&vnext, 1, MPI_DOUBLE, 0, world); } - if (FreeEndFinal && ireplica == nreplica-1 && (update->ntimestep == 0)) EFinalIni = veng; + if (FreeEndFinal && (ireplica == nreplica - 1) && (update->ntimestep == 0)) EFinalIni = veng; - if (ireplica == 0) vIni=veng; + if (ireplica == 0) vIni = veng; if (FreeEndFinalWithRespToEIni) { - if (cmode == SINGLE_PROC_DIRECT || cmode == SINGLE_PROC_MAP) { + if ((cmode == SINGLE_PROC_DIRECT) || (cmode == SINGLE_PROC_MAP)) { int procFirst; - procFirst=universe->root_proc[0]; - MPI_Bcast(&vIni,1,MPI_DOUBLE,procFirst,uworld); + procFirst = universe->root_proc[0]; + MPI_Bcast(&vIni, 1, MPI_DOUBLE, procFirst, uworld); } else { - if (me == 0) - MPI_Bcast(&vIni,1,MPI_DOUBLE,0,rootworld); + if (me == 0) MPI_Bcast(&vIni, 1, MPI_DOUBLE, 0, rootworld); - MPI_Bcast(&vIni,1,MPI_DOUBLE,0,world); + MPI_Bcast(&vIni, 1, MPI_DOUBLE, 0, world); } } if (FreeEndIni && ireplica == 0 && (update->ntimestep == 0)) EIniIni = veng; - + // communicate atoms to/from adjacent replicas to fill xprev,xnext inter_replica_comm(); // trigger potential energy computation on next timestep - pe->addstep(update->ntimestep+1); + pe->addstep(update->ntimestep + 1); double **x = atom->x; int *mask = atom->mask; @@ -331,25 +332,24 @@ void FixNEB::min_post_force(int /*vflag*/) dotgrad = gradlen = dotpath = dottangrad = 0.0; - if (ireplica == nreplica-1) { + if (ireplica == nreplica - 1) { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { delxp = x[i][0] - xprev[i][0]; delyp = x[i][1] - xprev[i][1]; delzp = x[i][2] - xprev[i][2]; - domain->minimum_image(delxp,delyp,delzp); - plen += delxp*delxp + delyp*delyp + delzp*delzp; - dottangrad += delxp* f[i][0]+ delyp*f[i][1]+delzp*f[i][2]; - gradlen += f[i][0]*f[i][0] + f[i][1]*f[i][1] + f[i][2]*f[i][2]; - if (FreeEndFinal||FreeEndFinalWithRespToEIni) { - tangent[i][0]=delxp; - tangent[i][1]=delyp; - tangent[i][2]=delzp; - tlen += tangent[i][0]*tangent[i][0] + - tangent[i][1]*tangent[i][1] + tangent[i][2]*tangent[i][2]; - dot += f[i][0]*tangent[i][0] + f[i][1]*tangent[i][1] + - f[i][2]*tangent[i][2]; + domain->minimum_image(delxp, delyp, delzp); + plen += delxp * delxp + delyp * delyp + delzp * delzp; + dottangrad += delxp * f[i][0] + delyp * f[i][1] + delzp * f[i][2]; + gradlen += f[i][0] * f[i][0] + f[i][1] * f[i][1] + f[i][2] * f[i][2]; + if (FreeEndFinal || FreeEndFinalWithRespToEIni) { + tangent[i][0] = delxp; + tangent[i][1] = delyp; + tangent[i][2] = delzp; + tlen += tangent[i][0] * tangent[i][0] + tangent[i][1] * tangent[i][1] + + tangent[i][2] * tangent[i][2]; + dot += f[i][0] * tangent[i][0] + f[i][1] * tangent[i][1] + f[i][2] * tangent[i][2]; } } @@ -359,44 +359,41 @@ void FixNEB::min_post_force(int /*vflag*/) delxn = xnext[i][0] - x[i][0]; delyn = xnext[i][1] - x[i][1]; delzn = xnext[i][2] - x[i][2]; - domain->minimum_image(delxn,delyn,delzn); - nlen += delxn*delxn + delyn*delyn + delzn*delzn; - gradnextlen += fnext[i][0]*fnext[i][0] - + fnext[i][1]*fnext[i][1] +fnext[i][2] * fnext[i][2]; - dotgrad += f[i][0]*fnext[i][0] - + f[i][1]*fnext[i][1] + f[i][2]*fnext[i][2]; - dottangrad += delxn*f[i][0]+ delyn*f[i][1] + delzn*f[i][2]; - gradlen += f[i][0]*f[i][0] + f[i][1]*f[i][1] + f[i][2]*f[i][2]; + domain->minimum_image(delxn, delyn, delzn); + nlen += delxn * delxn + delyn * delyn + delzn * delzn; + gradnextlen += + fnext[i][0] * fnext[i][0] + fnext[i][1] * fnext[i][1] + fnext[i][2] * fnext[i][2]; + dotgrad += f[i][0] * fnext[i][0] + f[i][1] * fnext[i][1] + f[i][2] * fnext[i][2]; + dottangrad += delxn * f[i][0] + delyn * f[i][1] + delzn * f[i][2]; + gradlen += f[i][0] * f[i][0] + f[i][1] * f[i][1] + f[i][2] * f[i][2]; if (FreeEndIni) { - tangent[i][0]=delxn; - tangent[i][1]=delyn; - tangent[i][2]=delzn; - tlen += tangent[i][0]*tangent[i][0] + - tangent[i][1]*tangent[i][1] + tangent[i][2]*tangent[i][2]; - dot += f[i][0]*tangent[i][0] + f[i][1]*tangent[i][1] + - f[i][2]*tangent[i][2]; + tangent[i][0] = delxn; + tangent[i][1] = delyn; + tangent[i][2] = delzn; + tlen += tangent[i][0] * tangent[i][0] + tangent[i][1] * tangent[i][1] + + tangent[i][2] * tangent[i][2]; + dot += f[i][0] * tangent[i][0] + f[i][1] * tangent[i][1] + f[i][2] * tangent[i][2]; } } } else { // not the first or last replica - double vmax = MAX(fabs(vnext-veng),fabs(vprev-veng)); - double vmin = MIN(fabs(vnext-veng),fabs(vprev-veng)); - + double vmax = MAX(fabs(vnext - veng), fabs(vprev - veng)); + double vmin = MIN(fabs(vnext - veng), fabs(vprev - veng)); for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { delxp = x[i][0] - xprev[i][0]; delyp = x[i][1] - xprev[i][1]; delzp = x[i][2] - xprev[i][2]; - domain->minimum_image(delxp,delyp,delzp); - plen += delxp*delxp + delyp*delyp + delzp*delzp; + domain->minimum_image(delxp, delyp, delzp); + plen += delxp * delxp + delyp * delyp + delzp * delzp; delxn = xnext[i][0] - x[i][0]; delyn = xnext[i][1] - x[i][1]; delzn = xnext[i][2] - x[i][2]; - domain->minimum_image(delxn,delyn,delzn); + domain->minimum_image(delxn, delyn, delzn); if (vnext > veng && veng > vprev) { tangent[i][0] = delxn; @@ -408,35 +405,33 @@ void FixNEB::min_post_force(int /*vflag*/) tangent[i][2] = delzp; } else { if (vnext > vprev) { - tangent[i][0] = vmax*delxn + vmin*delxp; - tangent[i][1] = vmax*delyn + vmin*delyp; - tangent[i][2] = vmax*delzn + vmin*delzp; + tangent[i][0] = vmax * delxn + vmin * delxp; + tangent[i][1] = vmax * delyn + vmin * delyp; + tangent[i][2] = vmax * delzn + vmin * delzp; } else if (vnext < vprev) { - tangent[i][0] = vmin*delxn + vmax*delxp; - tangent[i][1] = vmin*delyn + vmax*delyp; - tangent[i][2] = vmin*delzn + vmax*delzp; - } else { // vnext == vprev, e.g. for potentials that do not compute an energy + tangent[i][0] = vmin * delxn + vmax * delxp; + tangent[i][1] = vmin * delyn + vmax * delyp; + tangent[i][2] = vmin * delzn + vmax * delzp; + } else { // vnext == vprev, e.g. for potentials that do not compute an energy tangent[i][0] = delxn + delxp; tangent[i][1] = delyn + delyp; tangent[i][2] = delzn + delzp; } } - nlen += delxn*delxn + delyn*delyn + delzn*delzn; - tlen += tangent[i][0]*tangent[i][0] + - tangent[i][1]*tangent[i][1] + tangent[i][2]*tangent[i][2]; - gradlen += f[i][0]*f[i][0] + f[i][1]*f[i][1] + f[i][2]*f[i][2]; - dotpath += delxp*delxn + delyp*delyn + delzp*delzn; - dottangrad += tangent[i][0]*f[i][0] + - tangent[i][1]*f[i][1] + tangent[i][2]*f[i][2]; - gradnextlen += fnext[i][0]*fnext[i][0] + - fnext[i][1]*fnext[i][1] +fnext[i][2] * fnext[i][2]; - dotgrad += f[i][0]*fnext[i][0] + f[i][1]*fnext[i][1] + - f[i][2]*fnext[i][2]; + nlen += delxn * delxn + delyn * delyn + delzn * delzn; + tlen += tangent[i][0] * tangent[i][0] + tangent[i][1] * tangent[i][1] + + tangent[i][2] * tangent[i][2]; + gradlen += f[i][0] * f[i][0] + f[i][1] * f[i][1] + f[i][2] * f[i][2]; + dotpath += delxp * delxn + delyp * delyn + delzp * delzn; + dottangrad += tangent[i][0] * f[i][0] + tangent[i][1] * f[i][1] + tangent[i][2] * f[i][2]; + gradnextlen += + fnext[i][0] * fnext[i][0] + fnext[i][1] * fnext[i][1] + fnext[i][2] * fnext[i][2]; + dotgrad += f[i][0] * fnext[i][0] + f[i][1] * fnext[i][1] + f[i][2] * fnext[i][2]; - springF[i][0] = kspringPerp*(delxn-delxp); - springF[i][1] = kspringPerp*(delyn-delyp); - springF[i][2] = kspringPerp*(delzn-delzp); + springF[i][0] = kspringPerp * (delxn - delxp); + springF[i][1] = kspringPerp * (delyn - delyp); + springF[i][2] = kspringPerp * (delzn - delzp); } } @@ -449,7 +444,7 @@ void FixNEB::min_post_force(int /*vflag*/) bufin[5] = dotpath; bufin[6] = dottangrad; bufin[7] = dotgrad; - MPI_Allreduce(bufin,bufout,BUFSIZE,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(bufin, bufout, BUFSIZE, MPI_DOUBLE, MPI_SUM, world); nlen = sqrt(bufout[0]); plen = sqrt(bufout[1]); tlen = sqrt(bufout[2]); @@ -462,7 +457,7 @@ void FixNEB::min_post_force(int /*vflag*/) // normalize tangent vector if (tlen > 0.0) { - double tleninv = 1.0/tlen; + double tleninv = 1.0 / tlen; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { tangent[i][0] *= tleninv; @@ -473,125 +468,127 @@ void FixNEB::min_post_force(int /*vflag*/) // first or last replica has no change to forces, just return - if (ireplica > 0 && ireplica < nreplica-1) - dottangrad = dottangrad/(tlen*gradlen); - if (ireplica == 0) - dottangrad = dottangrad/(nlen*gradlen); - if (ireplica == nreplica-1) - dottangrad = dottangrad/(plen*gradlen); - if (ireplica < nreplica-1) - dotgrad = dotgrad /(gradlen*gradnextlen); + if (ireplica > 0 && ireplica < nreplica - 1) dottangrad = dottangrad / (tlen * gradlen); + if (ireplica == 0) dottangrad = dottangrad / (nlen * gradlen); + if (ireplica == nreplica - 1) dottangrad = dottangrad / (plen * gradlen); + if (ireplica < nreplica - 1) dotgrad = dotgrad / (gradlen * gradnextlen); if (FreeEndIni && ireplica == 0) { if (tlen > 0.0) { double dotall; - MPI_Allreduce(&dot,&dotall,1,MPI_DOUBLE,MPI_SUM,world); - dot=dotall/tlen; + MPI_Allreduce(&dot, &dotall, 1, MPI_DOUBLE, MPI_SUM, world); + dot = dotall / tlen; - if (dot<0) prefactor = -dot - kspringIni*(veng-EIniIni); - else prefactor = -dot + kspringIni*(veng-EIniIni); + if (dot < 0) + prefactor = -dot - kspringIni * (veng - EIniIni); + else + prefactor = -dot + kspringIni * (veng - EIniIni); for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - f[i][0] += prefactor *tangent[i][0]; - f[i][1] += prefactor *tangent[i][1]; - f[i][2] += prefactor *tangent[i][2]; + f[i][0] += prefactor * tangent[i][0]; + f[i][1] += prefactor * tangent[i][1]; + f[i][2] += prefactor * tangent[i][2]; } } } - if (FreeEndFinal && ireplica == nreplica -1) { + if (FreeEndFinal && ireplica == nreplica - 1) { if (tlen > 0.0) { double dotall; - MPI_Allreduce(&dot,&dotall,1,MPI_DOUBLE,MPI_SUM,world); - dot=dotall/tlen; + MPI_Allreduce(&dot, &dotall, 1, MPI_DOUBLE, MPI_SUM, world); + dot = dotall / tlen; - if (veng 0.0) { double dotall; - MPI_Allreduce(&dot,&dotall,1,MPI_DOUBLE,MPI_SUM,world); - dot=dotall/tlen; - if (veng xprev and x -> xnext if (cmode == SINGLE_PROC_DIRECT) { - if (ireplica > 0) - MPI_Irecv(xprev[0],3*nlocal,MPI_DOUBLE,procprev,0,uworld,&request); - if (ireplica < nreplica-1) - MPI_Send(x[0],3*nlocal,MPI_DOUBLE,procnext,0,uworld); - if (ireplica > 0) MPI_Wait(&request,MPI_STATUS_IGNORE); - if (ireplica < nreplica-1) - MPI_Irecv(xnext[0],3*nlocal,MPI_DOUBLE,procnext,0,uworld,&request); - if (ireplica > 0) - MPI_Send(x[0],3*nlocal,MPI_DOUBLE,procprev,0,uworld); - if (ireplica < nreplica-1) MPI_Wait(&request,MPI_STATUS_IGNORE); + if (ireplica > 0) MPI_Irecv(xprev[0], 3 * nlocal, MPI_DOUBLE, procprev, 0, uworld, &request); + if (ireplica < nreplica - 1) MPI_Send(x[0], 3 * nlocal, MPI_DOUBLE, procnext, 0, uworld); + if (ireplica > 0) MPI_Wait(&request, MPI_STATUS_IGNORE); + if (ireplica < nreplica - 1) + MPI_Irecv(xnext[0], 3 * nlocal, MPI_DOUBLE, procnext, 0, uworld, &request); + if (ireplica > 0) MPI_Send(x[0], 3 * nlocal, MPI_DOUBLE, procprev, 0, uworld); + if (ireplica < nreplica - 1) MPI_Wait(&request, MPI_STATUS_IGNORE); - if (ireplica < nreplica-1) - MPI_Irecv(fnext[0],3*nlocal,MPI_DOUBLE,procnext,0,uworld,&request); - if (ireplica > 0) - MPI_Send(f[0],3*nlocal,MPI_DOUBLE,procprev,0,uworld); - if (ireplica < nreplica-1) MPI_Wait(&request,MPI_STATUS_IGNORE); + if (ireplica < nreplica - 1) + MPI_Irecv(fnext[0], 3 * nlocal, MPI_DOUBLE, procnext, 0, uworld, &request); + if (ireplica > 0) MPI_Send(f[0], 3 * nlocal, MPI_DOUBLE, procprev, 0, uworld); + if (ireplica < nreplica - 1) MPI_Wait(&request, MPI_STATUS_IGNORE); return; } @@ -669,16 +661,16 @@ void FixNEB::inter_replica_comm() } if (ireplica > 0) { - MPI_Irecv(xrecv[0],3*nebatoms,MPI_DOUBLE,procprev,0,uworld,&requests[0]); - MPI_Irecv(tagrecv,nebatoms,MPI_LMP_TAGINT,procprev,0,uworld,&requests[1]); + MPI_Irecv(xrecv[0], 3 * nebatoms, MPI_DOUBLE, procprev, 0, uworld, &requests[0]); + MPI_Irecv(tagrecv, nebatoms, MPI_LMP_TAGINT, procprev, 0, uworld, &requests[1]); } - if (ireplica < nreplica-1) { - MPI_Send(xsend[0],3*nebatoms,MPI_DOUBLE,procnext,0,uworld); - MPI_Send(tagsend,nebatoms,MPI_LMP_TAGINT,procnext,0,uworld); + if (ireplica < nreplica - 1) { + MPI_Send(xsend[0], 3 * nebatoms, MPI_DOUBLE, procnext, 0, uworld); + MPI_Send(tagsend, nebatoms, MPI_LMP_TAGINT, procnext, 0, uworld); } if (ireplica > 0) { - MPI_Waitall(2,requests,statuses); + MPI_Waitall(2, requests, statuses); for (i = 0; i < nebatoms; i++) { m = atom->map(tagrecv[i]); xprev[m][0] = xrecv[i][0]; @@ -686,19 +678,19 @@ void FixNEB::inter_replica_comm() xprev[m][2] = xrecv[i][2]; } } - if (ireplica < nreplica-1) { - MPI_Irecv(xrecv[0],3*nebatoms,MPI_DOUBLE,procnext,0,uworld,&requests[0]); - MPI_Irecv(frecv[0],3*nebatoms,MPI_DOUBLE,procnext,0,uworld,&requests[0]); - MPI_Irecv(tagrecv,nebatoms,MPI_LMP_TAGINT,procnext,0,uworld,&requests[1]); + if (ireplica < nreplica - 1) { + MPI_Irecv(xrecv[0], 3 * nebatoms, MPI_DOUBLE, procnext, 0, uworld, &requests[0]); + MPI_Irecv(frecv[0], 3 * nebatoms, MPI_DOUBLE, procnext, 0, uworld, &requests[0]); + MPI_Irecv(tagrecv, nebatoms, MPI_LMP_TAGINT, procnext, 0, uworld, &requests[1]); } if (ireplica > 0) { - MPI_Send(xsend[0],3*nebatoms,MPI_DOUBLE,procprev,0,uworld); - MPI_Send(fsend[0],3*nebatoms,MPI_DOUBLE,procprev,0,uworld); - MPI_Send(tagsend,nebatoms,MPI_LMP_TAGINT,procprev,0,uworld); + MPI_Send(xsend[0], 3 * nebatoms, MPI_DOUBLE, procprev, 0, uworld); + MPI_Send(fsend[0], 3 * nebatoms, MPI_DOUBLE, procprev, 0, uworld); + MPI_Send(tagsend, nebatoms, MPI_LMP_TAGINT, procprev, 0, uworld); } - if (ireplica < nreplica-1) { - MPI_Waitall(2,requests,statuses); + if (ireplica < nreplica - 1) { + MPI_Waitall(2, requests, statuses); for (i = 0; i < nebatoms; i++) { m = atom->map(tagrecv[i]); xnext[m][0] = xrecv[i][0]; @@ -732,42 +724,39 @@ void FixNEB::inter_replica_comm() m++; } - MPI_Gather(&m,1,MPI_INT,counts,1,MPI_INT,0,world); + MPI_Gather(&m, 1, MPI_INT, counts, 1, MPI_INT, 0, world); displacements[0] = 0; - for (i = 0; i < nprocs-1; i++) - displacements[i+1] = displacements[i] + counts[i]; - MPI_Gatherv(tagsend,m,MPI_LMP_TAGINT, - tagsendall,counts,displacements,MPI_LMP_TAGINT,0,world); + for (i = 0; i < nprocs - 1; i++) displacements[i + 1] = displacements[i] + counts[i]; + MPI_Gatherv(tagsend, m, MPI_LMP_TAGINT, tagsendall, counts, displacements, MPI_LMP_TAGINT, 0, + world); for (i = 0; i < nprocs; i++) counts[i] *= 3; - for (i = 0; i < nprocs-1; i++) - displacements[i+1] = displacements[i] + counts[i]; + for (i = 0; i < nprocs - 1; i++) displacements[i + 1] = displacements[i] + counts[i]; if (xsend) { - MPI_Gatherv(xsend[0],3*m,MPI_DOUBLE, - xsendall[0],counts,displacements,MPI_DOUBLE,0,world); - MPI_Gatherv(fsend[0],3*m,MPI_DOUBLE, - fsendall[0],counts,displacements,MPI_DOUBLE,0,world); + MPI_Gatherv(xsend[0], 3 * m, MPI_DOUBLE, xsendall[0], counts, displacements, MPI_DOUBLE, 0, + world); + MPI_Gatherv(fsend[0], 3 * m, MPI_DOUBLE, fsendall[0], counts, displacements, MPI_DOUBLE, 0, + world); } else { - MPI_Gatherv(nullptr,3*m,MPI_DOUBLE, - xsendall[0],counts,displacements,MPI_DOUBLE,0,world); - MPI_Gatherv(nullptr,3*m,MPI_DOUBLE, - fsendall[0],counts,displacements,MPI_DOUBLE,0,world); + MPI_Gatherv(nullptr, 3 * m, MPI_DOUBLE, xsendall[0], counts, displacements, MPI_DOUBLE, 0, + world); + MPI_Gatherv(nullptr, 3 * m, MPI_DOUBLE, fsendall[0], counts, displacements, MPI_DOUBLE, 0, + world); } if (ireplica > 0 && me == 0) { - MPI_Irecv(xrecvall[0],3*nebatoms,MPI_DOUBLE,procprev,0,uworld,&requests[0]); - MPI_Irecv(tagrecvall,nebatoms,MPI_LMP_TAGINT,procprev,0,uworld, - &requests[1]); + MPI_Irecv(xrecvall[0], 3 * nebatoms, MPI_DOUBLE, procprev, 0, uworld, &requests[0]); + MPI_Irecv(tagrecvall, nebatoms, MPI_LMP_TAGINT, procprev, 0, uworld, &requests[1]); } - if (ireplica < nreplica-1 && me == 0) { - MPI_Send(xsendall[0],3*nebatoms,MPI_DOUBLE,procnext,0,uworld); - MPI_Send(tagsendall,nebatoms,MPI_LMP_TAGINT,procnext,0,uworld); + if (ireplica < nreplica - 1 && me == 0) { + MPI_Send(xsendall[0], 3 * nebatoms, MPI_DOUBLE, procnext, 0, uworld); + MPI_Send(tagsendall, nebatoms, MPI_LMP_TAGINT, procnext, 0, uworld); } if (ireplica > 0) { - if (me == 0) MPI_Waitall(2,requests,statuses); + if (me == 0) MPI_Waitall(2, requests, statuses); - MPI_Bcast(tagrecvall,nebatoms,MPI_INT,0,world); - MPI_Bcast(xrecvall[0],3*nebatoms,MPI_DOUBLE,0,world); + MPI_Bcast(tagrecvall, nebatoms, MPI_INT, 0, world); + MPI_Bcast(xrecvall[0], 3 * nebatoms, MPI_DOUBLE, 0, world); for (i = 0; i < nebatoms; i++) { m = atom->map(tagrecvall[i]); @@ -778,24 +767,23 @@ void FixNEB::inter_replica_comm() } } - if (ireplica < nreplica-1 && me == 0) { - MPI_Irecv(xrecvall[0],3*nebatoms,MPI_DOUBLE,procnext,0,uworld,&requests[0]); - MPI_Irecv(frecvall[0],3*nebatoms,MPI_DOUBLE,procnext,0,uworld,&requests[0]); - MPI_Irecv(tagrecvall,nebatoms,MPI_LMP_TAGINT,procnext,0,uworld, - &requests[1]); + if (ireplica < nreplica - 1 && me == 0) { + MPI_Irecv(xrecvall[0], 3 * nebatoms, MPI_DOUBLE, procnext, 0, uworld, &requests[0]); + MPI_Irecv(frecvall[0], 3 * nebatoms, MPI_DOUBLE, procnext, 0, uworld, &requests[0]); + MPI_Irecv(tagrecvall, nebatoms, MPI_LMP_TAGINT, procnext, 0, uworld, &requests[1]); } if (ireplica > 0 && me == 0) { - MPI_Send(xsendall[0],3*nebatoms,MPI_DOUBLE,procprev,0,uworld); - MPI_Send(fsendall[0],3*nebatoms,MPI_DOUBLE,procprev,0,uworld); - MPI_Send(tagsendall,nebatoms,MPI_LMP_TAGINT,procprev,0,uworld); + MPI_Send(xsendall[0], 3 * nebatoms, MPI_DOUBLE, procprev, 0, uworld); + MPI_Send(fsendall[0], 3 * nebatoms, MPI_DOUBLE, procprev, 0, uworld); + MPI_Send(tagsendall, nebatoms, MPI_LMP_TAGINT, procprev, 0, uworld); } - if (ireplica < nreplica-1) { - if (me == 0) MPI_Waitall(2,requests,statuses); + if (ireplica < nreplica - 1) { + if (me == 0) MPI_Waitall(2, requests, statuses); - MPI_Bcast(tagrecvall,nebatoms,MPI_INT,0,world); - MPI_Bcast(xrecvall[0],3*nebatoms,MPI_DOUBLE,0,world); - MPI_Bcast(frecvall[0],3*nebatoms,MPI_DOUBLE,0,world); + MPI_Bcast(tagrecvall, nebatoms, MPI_INT, 0, world); + MPI_Bcast(xrecvall[0], 3 * nebatoms, MPI_DOUBLE, 0, world); + MPI_Bcast(frecvall[0], 3 * nebatoms, MPI_DOUBLE, 0, world); for (i = 0; i < nebatoms; i++) { m = atom->map(tagrecvall[i]); @@ -816,58 +804,52 @@ Calculate ideal positions for parallel "ideal" or "equal" void FixNEB::calculate_ideal_positions() { // Skip unless "ideal" or "equal" - if (!(neb_mode==IDEAL||neb_mode==EQUAL)) return; + if (!((neb_mode == IDEAL) || (neb_mode == EQUAL))) return; - double Elentot,lentot,lenuntilClimber; - double meanDist,meanDistBeforeClimber,meanDistAfterClimber; + double lentot, lenuntilClimber; + double meanDist, meanDistBeforeClimber, meanDistAfterClimber; - if (neb_mode==EQUAL and rclimber>0.) { - if (cmode == SINGLE_PROC_DIRECT || cmode == SINGLE_PROC_MAP) { - MPI_Allgather(&veng,1,MPI_DOUBLE,&vengall[0],1,MPI_DOUBLE,uworld); + if ((neb_mode == EQUAL) && (rclimber > 0.0)) { + if ((cmode == SINGLE_PROC_DIRECT) || (cmode == SINGLE_PROC_MAP)) { + MPI_Allgather(&veng, 1, MPI_DOUBLE, &vengall[0], 1, MPI_DOUBLE, uworld); } else { - if (me == 0) - MPI_Allgather(&veng,1,MPI_DOUBLE,&vengall[0],1,MPI_DOUBLE,rootworld); - MPI_Bcast(vengall,nreplica,MPI_DOUBLE,0,world); + if (me == 0) MPI_Allgather(&veng, 1, MPI_DOUBLE, &vengall[0], 1, MPI_DOUBLE, rootworld); + MPI_Bcast(vengall, nreplica, MPI_DOUBLE, 0, world); } - for (int i = 0; i < nreplica-1; i++) - nlenall[i] = std::abs(vengall[i+1]-vengall[i]); - nlenall[nreplica-1] = 0.0; + for (int i = 0; i < nreplica - 1; i++) nlenall[i] = std::abs(vengall[i + 1] - vengall[i]); + nlenall[nreplica - 1] = 0.0; - } else if(neb_mode==IDEAL || neb_mode==EQUAL) { - if (cmode == SINGLE_PROC_DIRECT || cmode == SINGLE_PROC_MAP) { - MPI_Allgather(&nlen,1,MPI_DOUBLE,&nlenall[0],1,MPI_DOUBLE,uworld); + } else if ((neb_mode == IDEAL) || (neb_mode == EQUAL)) { + if ((cmode == SINGLE_PROC_DIRECT) || (cmode == SINGLE_PROC_MAP)) { + MPI_Allgather(&nlen, 1, MPI_DOUBLE, &nlenall[0], 1, MPI_DOUBLE, uworld); } else { - if (me == 0) - MPI_Allgather(&nlen,1,MPI_DOUBLE,&nlenall[0],1,MPI_DOUBLE,rootworld); - MPI_Bcast(nlenall,nreplica,MPI_DOUBLE,0,world); - } + if (me == 0) MPI_Allgather(&nlen, 1, MPI_DOUBLE, &nlenall[0], 1, MPI_DOUBLE, rootworld); + MPI_Bcast(nlenall, nreplica, MPI_DOUBLE, 0, world); + } } - + lentot = 0.; - for (int i = 0; i < nreplica; i++) - lentot += nlenall[i]; - + for (int i = 0; i < nreplica; i++) lentot += nlenall[i]; + actualPos = 0.; - for (int i = 0; i < ireplica; i++) - actualPos += nlenall[i]; - - meanDist = lentot/(nreplica-1); + for (int i = 0; i < ireplica; i++) actualPos += nlenall[i]; + + meanDist = lentot / (nreplica - 1); lenuntilClimber = 0.; - if (rclimber>0) { - for (int i = 0; i < rclimber; i++) - lenuntilClimber += nlenall[i]; - meanDistBeforeClimber = lenuntilClimber/rclimber; - meanDistAfterClimber = (lentot-lenuntilClimber)/(nreplica-rclimber-1); - if (ireplica 0) { + for (int i = 0; i < rclimber; i++) lenuntilClimber += nlenall[i]; + meanDistBeforeClimber = lenuntilClimber / rclimber; + meanDistAfterClimber = (lentot - lenuntilClimber) / (nreplica - rclimber - 1); + if (ireplica < rclimber) idealPos = ireplica * meanDistBeforeClimber; else - idealPos = lenuntilClimber+ (ireplica-rclimber)*meanDistAfterClimber; - } else idealPos = ireplica * meanDist; + idealPos = lenuntilClimber + (ireplica - rclimber) * meanDistAfterClimber; + } else + idealPos = ireplica * meanDist; idealPos /= meanDist; actualPos /= meanDist; - } /* ---------------------------------------------------------------------- @@ -885,11 +867,11 @@ void FixNEB::reallocate() memory->destroy(fnext); memory->destroy(springF); - memory->create(xprev,maxlocal,3,"neb:xprev"); - memory->create(xnext,maxlocal,3,"neb:xnext"); - memory->create(tangent,maxlocal,3,"neb:tangent"); - memory->create(fnext,maxlocal,3,"neb:fnext"); - memory->create(springF,maxlocal,3,"neb:springF"); + memory->create(xprev, maxlocal, 3, "neb:xprev"); + memory->create(xnext, maxlocal, 3, "neb:xnext"); + memory->create(tangent, maxlocal, 3, "neb:tangent"); + memory->create(fnext, maxlocal, 3, "neb:fnext"); + memory->create(springF, maxlocal, 3, "neb:springF"); if (cmode != SINGLE_PROC_DIRECT) { memory->destroy(xsend); @@ -898,20 +880,20 @@ void FixNEB::reallocate() memory->destroy(frecv); memory->destroy(tagsend); memory->destroy(tagrecv); - memory->create(xsend,maxlocal,3,"neb:xsend"); - memory->create(fsend,maxlocal,3,"neb:fsend"); - memory->create(xrecv,maxlocal,3,"neb:xrecv"); - memory->create(frecv,maxlocal,3,"neb:frecv"); - memory->create(tagsend,maxlocal,"neb:tagsend"); - memory->create(tagrecv,maxlocal,"neb:tagrecv"); + memory->create(xsend, maxlocal, 3, "neb:xsend"); + memory->create(fsend, maxlocal, 3, "neb:fsend"); + memory->create(xrecv, maxlocal, 3, "neb:xrecv"); + memory->create(frecv, maxlocal, 3, "neb:frecv"); + memory->create(tagsend, maxlocal, "neb:tagsend"); + memory->create(tagrecv, maxlocal, "neb:tagrecv"); } - if (neb_mode==IDEAL || neb_mode==EQUAL) { + if ((neb_mode == IDEAL) || (neb_mode == EQUAL)) { memory->destroy(nlenall); - memory->create(nlenall,nreplica,"neb:nlenall"); + memory->create(nlenall, nreplica, "neb:nlenall"); } - if (neb_mode==EQUAL) { + if (neb_mode == EQUAL) { memory->destroy(vengall); - memory->create(vengall,nreplica,"neb:vengall"); + memory->create(vengall, nreplica, "neb:vengall"); } } diff --git a/src/REPLICA/neb.cpp b/src/REPLICA/neb.cpp index b65a08d116..2e89e0047e 100644 --- a/src/REPLICA/neb.cpp +++ b/src/REPLICA/neb.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -42,6 +41,7 @@ using namespace MathConst; #define ATTRIBUTE_PERLINE 4 enum { NORMAL, TERSE, VERBOSE }; + /* ---------------------------------------------------------------------- */ NEB::NEB(LAMMPS *lmp) : Command(lmp), all(nullptr), rdist(nullptr) {} @@ -50,11 +50,12 @@ NEB::NEB(LAMMPS *lmp) : Command(lmp), all(nullptr), rdist(nullptr) {} internal NEB constructor, called from TAD ------------------------------------------------------------------------- */ -NEB::NEB(LAMMPS *lmp, double etol_in, double ftol_in, int n1steps_in, - int n2steps_in, int nevery_in, double *buf_init, double *buf_final) - : Command(lmp), fp(nullptr), all(nullptr), rdist(nullptr) +NEB::NEB(LAMMPS *lmp, double etol_in, double ftol_in, int n1steps_in, int n2steps_in, int nevery_in, + double *buf_init, double *buf_final) : + Command(lmp), + fp(nullptr), all(nullptr), rdist(nullptr) { - double delx,dely,delz; + double delx, dely, delz; etol = etol_in; ftol = ftol_in; @@ -68,22 +69,22 @@ NEB::NEB(LAMMPS *lmp, double etol_in, double ftol_in, int n1steps_in, ireplica = universe->iworld; me_universe = universe->me; uworld = universe->uworld; - MPI_Comm_rank(world,&me); + MPI_Comm_rank(world, &me); // generate linear interpolate replica - double fraction = ireplica/(nreplica-1.0); + double fraction = ireplica / (nreplica - 1.0); double **x = atom->x; int nlocal = atom->nlocal; int ii = 0; for (int i = 0; i < nlocal; i++) { delx = buf_final[ii] - buf_init[ii]; - dely = buf_final[ii+1] - buf_init[ii+1]; - delz = buf_final[ii+2] - buf_init[ii+2]; - domain->minimum_image(delx,dely,delz); - x[i][0] = buf_init[ii] + fraction*delx; - x[i][1] = buf_init[ii+1] + fraction*dely; - x[i][2] = buf_init[ii+2] + fraction*delz; + dely = buf_final[ii + 1] - buf_init[ii + 1]; + delz = buf_final[ii + 2] - buf_init[ii + 2]; + domain->minimum_image(delx, dely, delz); + x[i][0] = buf_init[ii] + fraction * delx; + x[i][1] = buf_init[ii + 1] + fraction * dely; + x[i][2] = buf_init[ii + 2] + fraction * delz; ii += 3; } } @@ -96,8 +97,10 @@ NEB::~NEB() memory->destroy(all); delete[] rdist; if (fp) { - if (compressed) platform::pclose(fp); - else fclose(fp); + if (compressed) + platform::pclose(fp); + else + fclose(fp); } } @@ -108,22 +111,22 @@ NEB::~NEB() void NEB::command(int narg, char **arg) { if (domain->box_exist == 0) - error->universe_all(FLERR,"NEB command before simulation box is defined"); + error->universe_all(FLERR, "NEB command before simulation box is defined"); - if (narg < 6) error->universe_all(FLERR,"Illegal NEB command: missing argument(s)"); + if (narg < 6) error->universe_all(FLERR, "Illegal NEB command: missing argument(s)"); - etol = utils::numeric(FLERR,arg[0],false,lmp); - ftol = utils::numeric(FLERR,arg[1],false,lmp); - n1steps = utils::inumeric(FLERR,arg[2],false,lmp); - n2steps = utils::inumeric(FLERR,arg[3],false,lmp); - nevery = utils::inumeric(FLERR,arg[4],false,lmp); + etol = utils::numeric(FLERR, arg[0], false, lmp); + ftol = utils::numeric(FLERR, arg[1], false, lmp); + n1steps = utils::inumeric(FLERR, arg[2], false, lmp); + n2steps = utils::inumeric(FLERR, arg[3], false, lmp); + nevery = utils::inumeric(FLERR, arg[4], false, lmp); // error checks if (etol < 0.0) error->universe_all(FLERR, fmt::format("Illegal NEB energy tolerance: {}", etol)); if (ftol < 0.0) error->universe_all(FLERR, fmt::format("Illegal NEB force tolerance: {}", ftol)); if (nevery <= 0) - error->universe_all(FLERR,fmt::format("Illegal NEB command every parameter: {}", nevery)); + error->universe_all(FLERR, fmt::format("Illegal NEB command every parameter: {}", nevery)); if (n1steps % nevery) error->all(FLERR, fmt::format("NEB N1 value {} incompatible with every {}", n1steps, nevery)); if (n2steps % nevery) @@ -135,41 +138,42 @@ void NEB::command(int narg, char **arg) ireplica = universe->iworld; me_universe = universe->me; uworld = universe->uworld; - MPI_Comm_rank(world,&me); + MPI_Comm_rank(world, &me); // error checks - if (nreplica == 1) error->universe_all(FLERR,"Cannot use NEB with a single replica"); + if (nreplica == 1) error->universe_all(FLERR, "Cannot use NEB with a single replica"); if (atom->map_style == Atom::MAP_NONE) - error->universe_all(FLERR,"Cannot use NEB without an atom map"); + error->universe_all(FLERR, "Cannot use NEB without an atom map"); // process file-style setting to setup initial configs for all replicas int iarg = 5; int filecmd = 0; - print_mode = NORMAL; // normal + print_mode = NORMAL; // normal while (iarg < narg) { - if (strcmp(arg[iarg],"final") == 0) { - if (iarg + 2 > narg) error->universe_all(FLERR,"Illegal NEB command: missing arguments"); - inpfile = arg[iarg+1]; - readfile(inpfile,0); + if (strcmp(arg[iarg], "final") == 0) { + if (iarg + 2 > narg) error->universe_all(FLERR, "Illegal NEB command: missing arguments"); + inpfile = arg[iarg + 1]; + readfile(inpfile, 0); filecmd = 1; iarg += 2; - } else if (strcmp(arg[iarg],"each") == 0) { - if (iarg + 2 > narg) error->universe_all(FLERR,"Illegal NEB command: missing arguments"); - inpfile = arg[iarg+1]; - readfile(inpfile,1); + } else if (strcmp(arg[iarg], "each") == 0) { + if (iarg + 2 > narg) error->universe_all(FLERR, "Illegal NEB command: missing arguments"); + inpfile = arg[iarg + 1]; + readfile(inpfile, 1); filecmd = 1; iarg += 2; - } else if (strcmp(arg[iarg],"none") == 0) { + } else if (strcmp(arg[iarg], "none") == 0) { filecmd = 1; ++iarg; - } else if (strcmp(arg[iarg],"verbose") == 0) { + } else if (strcmp(arg[iarg], "verbose") == 0) { print_mode = VERBOSE; ++iarg; - } else if (strcmp(arg[iarg],"terse") == 0) { + } else if (strcmp(arg[iarg], "terse") == 0) { print_mode = TERSE; ++iarg; - } else error->universe_all(FLERR,fmt::format("Unknown NEB command keyword: {}", arg[iarg])); + } else + error->universe_all(FLERR, fmt::format("Unknown NEB command keyword: {}", arg[iarg])); } if (!filecmd) error->universe_all(FLERR, "NEB is missing 'final', 'each', or 'none' keyword"); @@ -188,18 +192,22 @@ void NEB::run() // create MPI communicator for root proc from each world int color; - if (me == 0) color = 0; - else color = 1; - MPI_Comm_split(uworld,color,0,&roots); + if (me == 0) + color = 0; + else + color = 1; + MPI_Comm_split(uworld, color, 0, &roots); auto fixes = modify->get_fix_by_style("^neb$"); if (fixes.size() != 1) - error->universe_all(FLERR,"NEB requires use of exactly one fix neb instance"); + error->universe_all(FLERR, "NEB requires use of exactly one fix neb instance"); fneb = dynamic_cast(fixes[0]); - if (print_mode==VERBOSE) numall =7; - else numall = 4; - memory->create(all,nreplica,numall,"neb:all"); + if (print_mode == VERBOSE) + numall = 7; + else + numall = 4; + memory->create(all, nreplica, numall, "neb:all"); rdist = new double[nreplica]; // initialize LAMMPS @@ -212,63 +220,65 @@ void NEB::run() lmp->init(); if (update->minimize->searchflag) - error->universe_all(FLERR,"NEB requires a damped dynamics minimizer"); + error->universe_all(FLERR, "NEB requires a damped dynamics minimizer"); // setup regular NEB minimization FILE *uscreen = universe->uscreen; FILE *ulogfile = universe->ulogfile; - if (me_universe == 0 && uscreen) - fprintf(uscreen,"Setting up regular NEB ...\n"); + if (me_universe == 0 && uscreen) fprintf(uscreen, "Setting up regular NEB ...\n"); update->beginstep = update->firststep = update->ntimestep; update->endstep = update->laststep = update->firststep + n1steps; update->nsteps = n1steps; update->max_eval = n1steps; - if (update->laststep < 0) error->universe_all(FLERR,"Too many timesteps for NEB"); + if (update->laststep < 0) error->universe_all(FLERR, "Too many timesteps for NEB"); update->minimize->setup(); if (me_universe == 0) { if (uscreen) { - fmt::print(uscreen," Step {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} ", - "MaxReplicaForce", "MaxAtomForce", "GradV0","GradV1","GradVc","EBF", "EBR", "RDT"); - - if (print_mode!=TERSE) { + fmt::print(uscreen, " Step {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} ", + "MaxReplicaForce", "MaxAtomForce", "GradV0", "GradV1", "GradVc", "EBF", "EBR", + "RDT"); + + if (print_mode != TERSE) { for (int i = 1; i <= nreplica; ++i) - fmt::print(uscreen, "{:^14} {:^14} ", "RD"+std::to_string(i), "PE"+std::to_string(i)); + fmt::print(uscreen, "{:^14} {:^14} ", "RD" + std::to_string(i), "PE" + std::to_string(i)); } - if (print_mode==VERBOSE) { + if (print_mode == VERBOSE) { for (int i = 1; i <= nreplica; ++i) { auto idx = std::to_string(i); - fmt::print(uscreen, "{:^12}{:^12}{:^12} {:^12} {:^12}{:^12} ", - "pathangle"+idx, "angletangrad"+idx, "anglegrad"+idx, "gradV"+idx, - "RepForce"+idx, "MaxAtomForce"+idx); + fmt::print(uscreen, "{:^12}{:^12}{:^12} {:^12} {:^12}{:^12} ", "pathangle" + idx, + "angletangrad" + idx, "anglegrad" + idx, "gradV" + idx, "RepForce" + idx, + "MaxAtomForce" + idx); } } - fprintf(uscreen,"\n"); + fprintf(uscreen, "\n"); } if (ulogfile) { - fmt::print(ulogfile," Step {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} ", - "MaxReplicaForce", "MaxAtomForce", "GradV0","GradV1","GradVc","EBF", "EBR", "RDT"); - - if (print_mode!=TERSE) { + fmt::print(ulogfile, " Step {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} ", + "MaxReplicaForce", "MaxAtomForce", "GradV0", "GradV1", "GradVc", "EBF", "EBR", + "RDT"); + + if (print_mode != TERSE) { for (int i = 1; i <= nreplica; ++i) - fmt::print(ulogfile, "{:^14} {:^14} ", "RD"+std::to_string(i), "PE"+std::to_string(i)); + fmt::print(ulogfile, "{:^14} {:^14} ", "RD" + std::to_string(i), + "PE" + std::to_string(i)); } - if (print_mode==VERBOSE) { + if (print_mode == VERBOSE) { for (int i = 1; i <= nreplica; ++i) { auto idx = std::to_string(i); - fmt::print(ulogfile, "{:^12}{:^12}{:^12} {:^12} {:^12}{:^12} ", - "pathangle"+idx, "angletangrad"+idx, "anglegrad"+idx, "gradV"+idx, - "RepForce"+idx, "MaxAtomForce"+idx); + fmt::print(ulogfile, "{:^12}{:^12}{:^12} {:^12} {:^12}{:^12} ", "pathangle" + idx, + "angletangrad" + idx, "anglegrad" + idx, "gradV" + idx, "RepForce" + idx, + "MaxAtomForce" + idx); } } - fprintf(ulogfile,"\n"); + fprintf(ulogfile, "\n"); } } print_status(); @@ -308,66 +318,65 @@ void NEB::run() // setup climbing NEB minimization // must reinitialize minimizer so it re-creates its fix MINIMIZE - if (me_universe == 0 && uscreen) - fprintf(uscreen,"Setting up climbing ...\n"); + if (me_universe == 0 && uscreen) fprintf(uscreen, "Setting up climbing ...\n"); if (me_universe == 0) { - if (uscreen) - fprintf(uscreen,"Climbing replica = %d\n",top+1); - if (ulogfile) - fprintf(ulogfile,"Climbing replica = %d\n",top+1); + if (uscreen) fprintf(uscreen, "Climbing replica = %d\n", top + 1); + if (ulogfile) fprintf(ulogfile, "Climbing replica = %d\n", top + 1); } update->beginstep = update->firststep = update->ntimestep; update->endstep = update->laststep = update->firststep + n2steps; update->nsteps = n2steps; update->max_eval = n2steps; - if (update->laststep < 0) error->universe_all(FLERR,"Too many timesteps"); + if (update->laststep < 0) error->universe_all(FLERR, "Too many timesteps"); update->minimize->init(); fneb->rclimber = top; update->minimize->setup(); - if (me_universe == 0) { if (uscreen) { - fmt::print(uscreen," Step {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} ", - "MaxReplicaForce", "MaxAtomForce", "GradV0","GradV1","GradVc","EBF", "EBR", "RDT"); - - if (print_mode!=TERSE) { + fmt::print(uscreen, " Step {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} ", + "MaxReplicaForce", "MaxAtomForce", "GradV0", "GradV1", "GradVc", "EBF", "EBR", + "RDT"); + + if (print_mode != TERSE) { for (int i = 1; i <= nreplica; ++i) - fmt::print(uscreen, "{:^14} {:^14} ", "RD"+std::to_string(i), "PE"+std::to_string(i)); + fmt::print(uscreen, "{:^14} {:^14} ", "RD" + std::to_string(i), "PE" + std::to_string(i)); } - - if (print_mode==VERBOSE) { + + if (print_mode == VERBOSE) { for (int i = 1; i <= nreplica; ++i) { auto idx = std::to_string(i); - fmt::print(uscreen, "{:^12}{:^12}{:^12} {:^12} {:^12}{:^12} ", - "pathangle"+idx, "angletangrad"+idx, "anglegrad"+idx, "gradV"+idx, - "RepForce"+idx, "MaxAtomForce"+idx); + fmt::print(uscreen, "{:^12}{:^12}{:^12} {:^12} {:^12}{:^12} ", "pathangle" + idx, + "angletangrad" + idx, "anglegrad" + idx, "gradV" + idx, "RepForce" + idx, + "MaxAtomForce" + idx); } } - fprintf(uscreen,"\n"); + fprintf(uscreen, "\n"); } if (ulogfile) { - fmt::print(ulogfile," Step {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} ", - "MaxReplicaForce", "MaxAtomForce", "GradV0","GradV1","GradVc","EBF", "EBR", "RDT"); - - if (print_mode!=TERSE) { + fmt::print(ulogfile, " Step {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} {:^14} ", + "MaxReplicaForce", "MaxAtomForce", "GradV0", "GradV1", "GradVc", "EBF", "EBR", + "RDT"); + + if (print_mode != TERSE) { for (int i = 1; i <= nreplica; ++i) - fmt::print(ulogfile, "{:^14} {:^14} ", "RD"+std::to_string(i), "PE"+std::to_string(i)); + fmt::print(ulogfile, "{:^14} {:^14} ", "RD" + std::to_string(i), + "PE" + std::to_string(i)); } - if (print_mode==VERBOSE) { + if (print_mode == VERBOSE) { for (int i = 1; i <= nreplica; ++i) { auto idx = std::to_string(i); - fmt::print(ulogfile, "{:^12}{:^12}{:^12} {:^12} {:^12}{:^12} ", - "pathangle"+idx, "angletangrad"+idx, "anglegrad"+idx, "gradV"+idx, - "RepForce"+idx, "MaxAtomForce"+idx); + fmt::print(ulogfile, "{:^12}{:^12}{:^12} {:^12} {:^12}{:^12} ", "pathangle" + idx, + "angletangrad" + idx, "anglegrad" + idx, "gradV" + idx, "RepForce" + idx, + "MaxAtomForce" + idx); } } - fprintf(ulogfile,"\n"); + fprintf(ulogfile, "\n"); } } print_status(); @@ -415,14 +424,14 @@ void NEB::run() void NEB::readfile(char *file, int flag) { - int i,nchunk,eofflag,nlines; + int i, nchunk, eofflag, nlines; tagint tag; - char *eof,*start,*next,*buf; + char *eof, *start, *next, *buf; char line[MAXLINE]; - double delx,dely,delz; + double delx, dely, delz; if (me_universe == 0 && universe->uscreen) - fprintf(universe->uscreen,"Reading NEB coordinate file(s) ...\n"); + fprintf(universe->uscreen, "Reading NEB coordinate file(s) ...\n"); // flag = 0, universe root reads header of file, bcast to universe // flag = 1, each replica's root reads header of file, bcast to world @@ -432,38 +441,37 @@ void NEB::readfile(char *file, int flag) if (me_universe == 0) { open(file); while (true) { - eof = fgets(line,MAXLINE,fp); - if (eof == nullptr) error->one(FLERR,"Unexpected end of NEB file"); - start = &line[strspn(line," \t\n\v\f\r")]; + eof = fgets(line, MAXLINE, fp); + if (eof == nullptr) error->one(FLERR, "Unexpected end of NEB file"); + start = &line[strspn(line, " \t\n\v\f\r")]; if (*start != '\0' && *start != '#') break; } - int rv = sscanf(line,"%d",&nlines); + int rv = sscanf(line, "%d", &nlines); if (rv != 1) nlines = -1; } - MPI_Bcast(&nlines,1,MPI_INT,0,uworld); - if (nlines < 0) - error->universe_all(FLERR,"Incorrectly formatted NEB file"); + MPI_Bcast(&nlines, 1, MPI_INT, 0, uworld); + if (nlines < 0) error->universe_all(FLERR, "Incorrectly formatted NEB file"); } else { if (me == 0) { if (ireplica) { open(file); while (true) { - eof = fgets(line,MAXLINE,fp); - if (eof == nullptr) error->one(FLERR,"Unexpected end of NEB file"); - start = &line[strspn(line," \t\n\v\f\r")]; + eof = fgets(line, MAXLINE, fp); + if (eof == nullptr) error->one(FLERR, "Unexpected end of NEB file"); + start = &line[strspn(line, " \t\n\v\f\r")]; if (*start != '\0' && *start != '#') break; } - int rv = sscanf(line,"%d",&nlines); + int rv = sscanf(line, "%d", &nlines); if (rv != 1) nlines = -1; - } else nlines = 0; + } else + nlines = 0; } - MPI_Bcast(&nlines,1,MPI_INT,0,world); - if (nlines < 0) - error->universe_all(FLERR,"Incorrectly formatted NEB file"); + MPI_Bcast(&nlines, 1, MPI_INT, 0, world); + if (nlines < 0) error->universe_all(FLERR, "Incorrectly formatted NEB file"); } - auto buffer = new char[CHUNK*MAXLINE]; - double fraction = ireplica/(nreplica-1.0); + auto buffer = new char[CHUNK * MAXLINE]; + double fraction = ireplica / (nreplica - 1.0); double **x = atom->x; int nlocal = atom->nlocal; @@ -473,32 +481,31 @@ void NEB::readfile(char *file, int flag) int ncount = 0, nread = 0; while (nread < nlines) { - nchunk = MIN(nlines-nread,CHUNK); + nchunk = MIN(nlines - nread, CHUNK); if (flag == 0) - eofflag = utils::read_lines_from_file(fp,nchunk,MAXLINE,buffer, - universe->me,universe->uworld); + eofflag = + utils::read_lines_from_file(fp, nchunk, MAXLINE, buffer, universe->me, universe->uworld); else - eofflag = utils::read_lines_from_file(fp,nchunk,MAXLINE,buffer,me,world); - if (eofflag) error->all(FLERR,"Unexpected end of NEB file"); + eofflag = utils::read_lines_from_file(fp, nchunk, MAXLINE, buffer, me, world); + if (eofflag) error->all(FLERR, "Unexpected end of NEB file"); buf = buffer; - next = strchr(buf,'\n'); + next = strchr(buf, '\n'); *next = '\0'; int nwords = utils::count_words(utils::trim_comment(buf)); *next = '\n'; - if (nwords != ATTRIBUTE_PERLINE) - error->all(FLERR,"Incorrect atom format in NEB file"); + if (nwords != ATTRIBUTE_PERLINE) error->all(FLERR, "Incorrect atom format in NEB file"); // loop over lines of atom coords // tokenize the line into values for (i = 0; i < nchunk; i++) { - next = strchr(buf,'\n'); + next = strchr(buf, '\n'); *next = '\0'; try { - ValueTokenizer values(buf," \t\n\r\f"); + ValueTokenizer values(buf, " \t\n\r\f"); // adjust atom coord based on replica fraction // for flag = 0, interpolate for intermediate and final replicas @@ -519,12 +526,12 @@ void NEB::readfile(char *file, int flag) dely = values.next_double() - x[m][1]; delz = values.next_double() - x[m][2]; - domain->minimum_image(delx,dely,delz); + domain->minimum_image(delx, dely, delz); if (flag == 0) { - x[m][0] += fraction*delx; - x[m][1] += fraction*dely; - x[m][2] += fraction*delz; + x[m][0] += fraction * delx; + x[m][1] += fraction * dely; + x[m][2] += fraction * delz; } else { x[m][0] += delx; x[m][1] += dely; @@ -532,7 +539,7 @@ void NEB::readfile(char *file, int flag) } } } catch (std::exception &e) { - error->universe_one(FLERR,"Incorrectly formatted NEB file: " + std::string(e.what())); + error->universe_one(FLERR, "Incorrectly formatted NEB file: " + std::string(e.what())); } buf = next + 1; } @@ -543,14 +550,12 @@ void NEB::readfile(char *file, int flag) if (flag == 0) { int ntotal; - MPI_Allreduce(&ncount,&ntotal,1,MPI_INT,MPI_SUM,uworld); - if (ntotal != nreplica*nlines) - error->universe_all(FLERR,"Invalid atom IDs in NEB file"); + MPI_Allreduce(&ncount, &ntotal, 1, MPI_INT, MPI_SUM, uworld); + if (ntotal != nreplica * nlines) error->universe_all(FLERR, "Invalid atom IDs in NEB file"); } else { int ntotal; - MPI_Allreduce(&ncount,&ntotal,1,MPI_INT,MPI_SUM,world); - if (ntotal != nlines) - error->all(FLERR,"Invalid atom IDs in NEB file"); + MPI_Allreduce(&ncount, &ntotal, 1, MPI_INT, MPI_SUM, world); + if (ntotal != nlines) error->all(FLERR, "Invalid atom IDs in NEB file"); } // clean up @@ -558,13 +563,17 @@ void NEB::readfile(char *file, int flag) if (flag == 0) { if (me_universe == 0) { - if (compressed) platform::pclose(fp); - else fclose(fp); + if (compressed) + platform::pclose(fp); + else + fclose(fp); } } else { if (me == 0 && ireplica) { - if (compressed) platform::pclose(fp); - else fclose(fp); + if (compressed) + platform::pclose(fp); + else + fclose(fp); } } fp = nullptr; @@ -581,11 +590,11 @@ void NEB::open(char *file) if (platform::has_compress_extension(file)) { compressed = 1; fp = platform::compressed_read(file); - if (!fp) error->one(FLERR,"Cannot open compressed file {}: {}", file, utils::getsyserror()); - } else fp = fopen(file,"r"); + if (!fp) error->one(FLERR, "Cannot open compressed file {}: {}", file, utils::getsyserror()); + } else + fp = fopen(file, "r"); - if (fp == nullptr) - error->one(FLERR,"Cannot open file {}: {}", file, utils::getsyserror()); + if (fp == nullptr) error->one(FLERR, "Cannot open file {}: {}", file, utils::getsyserror()); } /* ---------------------------------------------------------------------- @@ -597,16 +606,16 @@ void NEB::print_status() { double fnorm2 = sqrt(update->minimize->fnorm_sqr()); double fmaxreplica; - MPI_Allreduce(&fnorm2,&fmaxreplica,1,MPI_DOUBLE,MPI_MAX,roots); + MPI_Allreduce(&fnorm2, &fmaxreplica, 1, MPI_DOUBLE, MPI_MAX, roots); double fnorminf = update->minimize->fnorm_inf(); double fmaxatom; - MPI_Allreduce(&fnorminf,&fmaxatom,1,MPI_DOUBLE,MPI_MAX,roots); + MPI_Allreduce(&fnorminf, &fmaxatom, 1, MPI_DOUBLE, MPI_MAX, roots); - if (print_mode==VERBOSE) { + if (print_mode == VERBOSE) { freplica = new double[nreplica]; - MPI_Allgather(&fnorm2,1,MPI_DOUBLE,&freplica[0],1,MPI_DOUBLE,roots); + MPI_Allgather(&fnorm2, 1, MPI_DOUBLE, &freplica[0], 1, MPI_DOUBLE, roots); fmaxatomInRepl = new double[nreplica]; - MPI_Allgather(&fnorminf,1,MPI_DOUBLE,&fmaxatomInRepl[0],1,MPI_DOUBLE,roots); + MPI_Allgather(&fnorminf, 1, MPI_DOUBLE, &fmaxatomInRepl[0], 1, MPI_DOUBLE, roots); } double one[7]; @@ -615,23 +624,20 @@ void NEB::print_status() one[2] = fneb->nlen; one[3] = fneb->gradlen; - if (print_mode==VERBOSE) { + if (print_mode == VERBOSE) { one[4] = fneb->dotpath; one[5] = fneb->dottangrad; one[6] = fneb->dotgrad; } if (output->thermo->normflag) one[0] /= atom->natoms; - if (me == 0) - MPI_Allgather(one,numall,MPI_DOUBLE,&all[0][0],numall,MPI_DOUBLE,roots); - MPI_Bcast(&all[0][0],numall*nreplica,MPI_DOUBLE,0,world); + if (me == 0) MPI_Allgather(one, numall, MPI_DOUBLE, &all[0][0], numall, MPI_DOUBLE, roots); + MPI_Bcast(&all[0][0], numall * nreplica, MPI_DOUBLE, 0, world); rdist[0] = 0.0; - for (int i = 1; i < nreplica; i++) - rdist[i] = rdist[i-1] + all[i][1]; - double endpt = rdist[nreplica-1] = rdist[nreplica-2] + all[nreplica-2][2]; - for (int i = 1; i < nreplica; i++) - rdist[i] /= endpt; + for (int i = 1; i < nreplica; i++) rdist[i] = rdist[i - 1] + all[i][1]; + double endpt = rdist[nreplica - 1] = rdist[nreplica - 2] + all[nreplica - 2][2]; + for (int i = 1; i < nreplica; i++) rdist[i] /= endpt; // look up GradV for the initial, final, and climbing replicas // these are identical to fnorm2, but to be safe we @@ -642,13 +648,13 @@ void NEB::print_status() int irep; irep = 0; gradvnorm0 = all[irep][3]; - irep = nreplica-1; + irep = nreplica - 1; gradvnorm1 = all[irep][3]; irep = fneb->rclimber; if (irep > -1) { gradvnormc = all[irep][3]; - ebf = all[irep][0]-all[0][0]; - ebr = all[irep][0]-all[nreplica-1][0]; + ebf = all[irep][0] - all[0][0]; + ebr = all[irep][0] - all[nreplica - 1][0]; } else { double vmax = all[0][0]; int top = 0; @@ -659,29 +665,32 @@ void NEB::print_status() } irep = top; gradvnormc = all[irep][3]; - ebf = all[irep][0]-all[0][0]; - ebr = all[irep][0]-all[nreplica-1][0]; + ebf = all[irep][0] - all[0][0]; + ebr = all[irep][0] - all[nreplica - 1][0]; } if (me_universe == 0) { - constexpr double todeg=180.0/MY_PI; - std::string mesg = fmt::format("{:10} {:<14.8g} {:<14.8g} ",update->ntimestep,fmaxreplica,fmaxatom); - mesg += fmt::format("{:<14.8g} {:<14.8g} {:<14.8g} ",gradvnorm0,gradvnorm1,gradvnormc); - mesg += fmt::format("{:<14.8g} {:<14.8g} {:<14.8g} ",ebf,ebr,endpt); - if (print_mode!=TERSE) { - for (int i = 0; i < nreplica; i++) mesg += fmt::format("{:<14.8g} {:<14.8g} ",rdist[i],all[i][0]); + constexpr double todeg = 180.0 / MY_PI; + std::string mesg = + fmt::format("{:10} {:<14.8g} {:<14.8g} ", update->ntimestep, fmaxreplica, fmaxatom); + mesg += fmt::format("{:<14.8g} {:<14.8g} {:<14.8g} ", gradvnorm0, gradvnorm1, gradvnormc); + mesg += fmt::format("{:<14.8g} {:<14.8g} {:<14.8g} ", ebf, ebr, endpt); + if (print_mode != TERSE) { + for (int i = 0; i < nreplica; i++) + mesg += fmt::format("{:<14.8g} {:<14.8g} ", rdist[i], all[i][0]); } - if (print_mode==VERBOSE) { - mesg += fmt::format("{:<12.5g} {:<12.5g} {:<12.5g} {:<12.5g} {:<12.5g} {:<12.5g}", - NAN,180-acos(all[0][5])*todeg,180-acos(all[0][6])*todeg, - all[0][3],freplica[0],fmaxatomInRepl[0]); - for (int i = 1; i < nreplica-1; i++) - mesg += fmt::format("{:<12.5g} {:<12.5g} {:<12.5g} {:<12.5g} {:<12.5g} {:<12.5g}", - 180-acos(all[i][4])*todeg,180-acos(all[i][5])*todeg, - 180-acos(all[i][6])*todeg,all[i][3],freplica[i],fmaxatomInRepl[i]); - mesg += fmt::format("{:<12.5g} {:<12.5g} {:<12.5g} {:<12.5g} {:<12.5g} {:<12.5g}", - NAN,180-acos(all[nreplica-1][5])*todeg,NAN,all[nreplica-1][3], - freplica[nreplica-1],fmaxatomInRepl[nreplica-1]); + if (print_mode == VERBOSE) { + mesg += fmt::format("{:<12.5g} {:<12.5g} {:<12.5g} {:<12.5g} {:<12.5g} {:<12.5g}", NAN, + 180 - acos(all[0][5]) * todeg, 180 - acos(all[0][6]) * todeg, all[0][3], + freplica[0], fmaxatomInRepl[0]); + for (int i = 1; i < nreplica - 1; i++) + mesg += + fmt::format("{:<12.5g} {:<12.5g} {:<12.5g} {:<12.5g} {:<12.5g} {:<12.5g}", + 180 - acos(all[i][4]) * todeg, 180 - acos(all[i][5]) * todeg, + 180 - acos(all[i][6]) * todeg, all[i][3], freplica[i], fmaxatomInRepl[i]); + mesg += fmt::format("{:<12.5g} {:<12.5g} {:<12.5g} {:<12.5g} {:<12.5g} {:<12.5g}", NAN, + 180 - acos(all[nreplica - 1][5]) * todeg, NAN, all[nreplica - 1][3], + freplica[nreplica - 1], fmaxatomInRepl[nreplica - 1]); } mesg += "\n"; @@ -691,7 +700,7 @@ void NEB::print_status() fflush(universe->ulogfile); } } - if (print_mode==VERBOSE) { + if (print_mode == VERBOSE) { delete[] freplica; delete[] fmaxatomInRepl; } From f27c7a9135cb33abfb043a25c5af8cc2ca5b69e1 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 12 Jan 2023 22:58:16 -0500 Subject: [PATCH 20/25] rework neb docs to use .. math:: and :math: in sphinx --- doc/src/fix_neb.rst | 114 ++++++++++++++++++++++++-------------------- doc/src/neb.rst | 57 +++++++++++----------- 2 files changed, 91 insertions(+), 80 deletions(-) diff --git a/doc/src/fix_neb.rst b/doc/src/fix_neb.rst index 0c6d0323d7..8e9ba3939c 100644 --- a/doc/src/fix_neb.rst +++ b/doc/src/fix_neb.rst @@ -60,37 +60,37 @@ replica having the highest energy relaxes toward the saddle point relaxation is performed. A key purpose of the nudging forces is to keep the replicas equally -spaced. During the NEB calculation, the 3N-length vector of -interatomic force Fi = -Grad(V) for each replica I is altered. For -all intermediate replicas (i.e. for 1 < I < N, except the climbing -replica) the force vector becomes: +spaced. During the NEB calculation, the :math:`3N`-length vector of +interatomic force :math:`F_i = -\nabla V` for each replica *i* is +altered. For all intermediate replicas (i.e. for :math:`1 < i < N`, +except the climbing replica) the force vector becomes: -.. parsed-literal:: +.. math:: - Fi = -Grad(V) + (Grad(V) dot T') T' + Fnudge_parallel + Fnudge_perp + F_i = -\nabla V + (\nabla V \cdot T') T' + F_\parallel + F_\perp -T' is the unit "tangent" vector for replica I and is a function of Ri, -Ri-1, Ri+1, and the potential energy of the 3 replicas; it points -roughly in the direction of (Ri+i - Ri-1); see the -:ref:`(Henkelman1) ` paper for details. Ri are the atomic -coordinates of replica I; Ri-1 and Ri+1 are the coordinates of its -neighbor replicas. The term (Grad(V) dot T') is used to remove the -component of the gradient parallel to the path which would tend to -distribute the replica unevenly along the path. Fnudge_parallel is an -artificial nudging force which is applied only in the tangent -direction and which maintains the equal spacing between replicas (see -below for more information). Fnudge_perp is an optional artificial -spring which is applied in a direction perpendicular to the tangent -direction and which prevent the paths from forming acute kinks (see -below for more information). +T' is the unit "tangent" vector for replica *i* and is a function of +:math:`R_i, R_{i-1}, R_{i+1}`, and the potential energy of the 3 +replicas; it points roughly in the direction of :math:`R_{i+i} - +R_{i-1}`; see the :ref:`(Henkelman1) ` paper for details. +:math:`R_i` are the atomic coordinates of replica *i*; :math:`R_{i-1}` +and :math:`R_{i+1}` are the coordinates of its neighbor replicas. The +term :math:`\nabla V \cdot T'` is used to remove the component of the +gradient parallel to the path which would tend to distribute the replica +unevenly along the path. :math:`F_\parallel` is an artificial nudging +force which is applied only in the tangent direction and which maintains +the equal spacing between replicas (see below for more information). +:math:`F_\perp` is an optional artificial spring which is applied in a +direction perpendicular to the tangent direction and which prevent the +paths from forming acute kinks (see below for more information). -In the second stage of the NEB calculation, the interatomic force Fi +In the second stage of the NEB calculation, the interatomic force :math:`F_i` for the climbing replica (the replica of highest energy after the first stage) is changed to: -.. parsed-literal:: +.. math:: - Fi = -Grad(V) + 2 (Grad(V) dot T') T' + Fnudge_perp + F_i = -\nabla V + 2 (\nabla V \cdot T') T' + F_\perp and the relaxation procedure is continued to a new converged MEP. @@ -101,26 +101,26 @@ computed. With a value of *neigh*, the parallel nudging force is computed as in :ref:`(Henkelman1) ` by connecting each intermediate replica with the previous and the next image: -.. parsed-literal:: +.. math:: - Fnudge_parallel = *Kspring* \* (\|Ri+1 - Ri\| - \|Ri - Ri-1\|) + F_\parallel = Kspring \cdot \left(\left|R_{i+1} - R_i\right| - \left|R_i - R_{i-1}\right|\right) -Note that in this case the specified *Kspring* is in force/distance -units. +Note that in this case the specified *Kspring* is in +force/distance units. With a value of *ideal*, the spring force is computed as suggested in ref`(WeinanE) ` -.. parsed-literal:: +.. math:: - Fnudge_parallel = -\ *Kspring* \* (RD-RDideal) / (2 \* meanDist) + F_\parallel = -Kspring \cdot (RD - RD_{ideal}) / (2 \cdot meanDist) -where RD is the "reaction coordinate" see :doc:`neb ` section, and -RDideal is the ideal RD for which all the images are equally spaced. -I.e. RDideal = (I-1)\*meanDist when the climbing replica is off, where -I is the replica number). The meanDist is the average distance -between replicas. Note that in this case the specified *Kspring* is -in force units. +where *RD* is the "reaction coordinate" see :doc:`neb ` section, +and :math:`RD_{ideal}` is the ideal *RD* for which all the images are +equally spaced. I.e. :math:`RD_{ideal} = (i-1) \cdot meanDist` when the +climbing replica is off, where *i* is the replica number). The +*meanDist* is the average distance between replicas. Note that in this +case the specified *Kspring* is in force units. Note that the *ideal* form of nudging can often be more effective at keeping the replicas equally spaced. @@ -129,15 +129,20 @@ With a value of *equal* the spring force is computed as for *ideal* before the climbing stage, then is computed to promote equidistant spacing in energy rather than distance: -.. parsed-literal:: - Fnudge_parallel = -\ *Kspring* \* (ED-EDideal) / (2 \* meanEDist) +.. math:: -where ED is the cumulative sum of absolute energy differences -ED=sum(i`, by providing optimal -quadrature points. + F_\parallel = -Kspring \cdot (ED - ED_{ideal}) / (2 \cdot meanEDist) + +where *ED* is the cumulative sum of absolute energy differences: + +.. math:: + + ED = \sum_{i`, by providing +optimal quadrature points. ---------- @@ -150,14 +155,16 @@ resolution is poor. I.e. when few replicas are used; see The perpendicular spring force is given by -.. parsed-literal:: +.. math:: - Fnudge_perp = *Kspring2* \* F(Ri-1,Ri,Ri+1) (Ri+1 + Ri-1 - 2 Ri) + F_\perp = K_{spring2} \cdot F(R_{i-1},R_i,R_{i+1}) (R_{i+1} + R_{i-1} - 2 R_i) -where *Kspring2* is the specified value. F(Ri-1 Ri R+1) is a smooth -scalar function of the angle Ri-1 Ri Ri+1. It is equal to 0.0 when -the path is straight and is equal to 1 when the angle Ri-1 Ri Ri+1 is -acute. F(Ri-1 Ri R+1) is defined in :ref:`(Jonsson) `. +where *Kspring2* is the specified value. :math:`F(R_{i-1}, R_i, +R_{i+1})` is a smooth scalar function of the angle :math:`R_{i-1} R_i +R_{i+1}`. It is equal to 0.0 when the path is straight and is equal to +1 when the angle :math:`R_{i-1} R_i R_{i+1}` is acute. +:math:`F(R_{i-1}, R_i, R_{i+1})` is defined in :ref:`(Jonsson) +`. If *Kspring2* is set to 0.0 (the default) then no perpendicular spring force is added. @@ -171,12 +178,13 @@ forces can be applied to the first and/or last replicas, to enable them to relax toward a MEP while constraining their energy E to the target energy ETarget. -If ETarget>E, the interatomic force Fi for the specified replica becomes: +If :math:`E_{Target} > E`, the interatomic force :math:`F_i` for the +specified replica becomes: -.. parsed-literal:: +.. math:: - Fi = -Grad(V) + (Grad(V) dot T' + (E-ETarget)\*Kspring3) T', *when* Grad(V) dot T' < 0 - Fi = -Grad(V) + (Grad(V) dot T' + (ETarget- E)\*Kspring3) T', *when* Grad(V) dot T' > 0 + F_i & = -\nabla V + (\nabla V \cdot T' + (E - E_{Target}) \cdot K_{spring3}) T', \qquad \textrm{when} \quad \nabla V \cdot T' < 0 \\ + F_i & = -\nabla V + (\nabla V \cdot T' + (E_{Target} - E) \cdot K_{spring3}) T', \qquad \textrm{when} \quad \nabla V \cdot T' > 0 The "spring" constant on the difference in energies is the specified *Kspring3* value. diff --git a/doc/src/neb.rst b/doc/src/neb.rst index 1ba1bd1cc2..33c9ac8c27 100644 --- a/doc/src/neb.rst +++ b/doc/src/neb.rst @@ -29,7 +29,7 @@ Syntax *none* arg = no argument all replicas assumed to already have their initial coords -keyword = *verbose* or *terse* +* keyword = *verbose* or *terse* Examples """""""" @@ -47,19 +47,21 @@ Perform a nudged elastic band (NEB) calculation using multiple replicas of a system. Two or more replicas must be used; the first and last are the end points of the transition path. -NEB is a method for finding both the atomic configurations and height -of the energy barrier associated with a transition state, e.g. for an -atom to perform a diffusive hop from one energy basin to another in a +NEB is a method for finding both the atomic configurations and height of +the energy barrier associated with a transition state, e.g. for an atom +to perform a diffusive hop from one energy basin to another in a coordinated fashion with its neighbors. The implementation in LAMMPS -follows the discussion in these 4 papers: :ref:`(HenkelmanA) `, -:ref:`(HenkelmanB) `, :ref:`(Nakano) ` and :ref:`(Maras) `. +follows the discussion in these 4 papers: :ref:`(HenkelmanA) +`, :ref:`(HenkelmanB) `, :ref:`(Nakano) +` and :ref:`(Maras) `. Each replica runs on a partition of one or more processors. Processor -partitions are defined at run-time using the :doc:`-partition command-line switch `. Note that if you have MPI installed, you -can run a multi-replica simulation with more replicas (partitions) -than you have physical processors, e.g you can run a 10-replica -simulation on just one or two processors. You will simply not get the -performance speed-up you would see with one or more physical +partitions are defined at run-time using the :doc:`-partition +command-line switch `. Note that if you have MPI +installed, you can run a multi-replica simulation with more replicas +(partitions) than you have physical processors, e.g you can run a +10-replica simulation on just one or two processors. You will simply +not get the performance speed-up you would see with one or more physical processors per replica. See the :doc:`Howto replica ` doc page for further discussion. @@ -330,25 +332,26 @@ the fix neb command. The forward (reverse) energy barrier is the potential energy of the highest replica minus the energy of the first (last) replica. -Supplementary information for all replicas can be printed out to the -screen and master log.lammps file by adding the verbose keyword. This -information include the following. The "path angle" (pathangle) for -the replica i which is the angle between the 3N-length vectors (Ri-1 - -Ri) and (Ri+1 - Ri) (where Ri is the atomic coordinates of replica -i). A "path angle" of 180 indicates that replicas i-1, i and i+1 are +Supplementary information for all replicas can be printed out by adding +the *verbose* keyword. This information include the following. The +"path angle" (pathangle) for the replica i which is the angle between +the 3N-length vectors :math:`(R_{i-1} - R_i)` and :math:`(R_{i+1} - +R_i)` (where :math:`R_i` is the atomic coordinates of replica *i*). A +"path angle" of 180 indicates that replicas *i*-1, *i* and *i*+1 are aligned. "angletangrad" is the angle between the 3N-length tangent -vector and the 3N-length force vector at image i. The tangent vector -is calculated as in :ref:`(HenkelmanA) ` for all intermediate -replicas and at R2 - R1 and RM - RM-1 for the first and last replica, -respectively. "anglegrad" is the angle between the 3N-length energy -gradient vector of replica i and that of replica i+1. It is not -defined for the final replica and reads nan. gradV is the norm of the -energy gradient of image i. ReplicaForce is the two-norm of the -3N-length force vector (including nudging forces) for replica i. -MaxAtomForce is the maximum force component of any atom in replica i. +vector and the 3N-length force vector at image *i*. The tangent vector +is calculated as in :ref:`(HenkelmanA) ` for all +intermediate replicas and at R2 - R1 and RM - RM-1 for the first and +last replica, respectively. "anglegrad" is the angle between the +3N-length energy gradient vector of replica *i* and that of replica +*i*+1. It is not defined for the final replica and reads nan. gradV is +the norm of the energy gradient of image *i* (:math:`\nabla V`). +ReplicaForce is the two-norm of the 3N-length force vector (including +nudging forces) for replica *i*. MaxAtomForce is the maximum force +component of any atom in replica *i*. Alternatively, a restricted print out can be obtained by adding the -terse keyword, which omits per-replica information. This typically +*terse* keyword, which omits per-replica information. This typically fits on one line of a command terminal, aiding visual inspection of an ongoing NEB calculation. From b97a0c62e42b1cbdcf9a86b347ce950b04cffa44 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 12 Jan 2023 23:00:52 -0500 Subject: [PATCH 21/25] whitespace --- doc/src/fix_neb.rst | 2 +- doc/src/neb.rst | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/doc/src/fix_neb.rst b/doc/src/fix_neb.rst index 8e9ba3939c..4d0d301540 100644 --- a/doc/src/fix_neb.rst +++ b/doc/src/fix_neb.rst @@ -126,7 +126,7 @@ Note that the *ideal* form of nudging can often be more effective at keeping the replicas equally spaced. With a value of *equal* the spring force is computed as for *ideal* -before the climbing stage, then is computed to promote equidistant +before the climbing stage, then is computed to promote equidistant spacing in energy rather than distance: .. math:: diff --git a/doc/src/neb.rst b/doc/src/neb.rst index 33c9ac8c27..4703cc7f60 100644 --- a/doc/src/neb.rst +++ b/doc/src/neb.rst @@ -350,9 +350,9 @@ ReplicaForce is the two-norm of the 3N-length force vector (including nudging forces) for replica *i*. MaxAtomForce is the maximum force component of any atom in replica *i*. -Alternatively, a restricted print out can be obtained by adding the +Alternatively, a restricted print out can be obtained by adding the *terse* keyword, which omits per-replica information. This typically -fits on one line of a command terminal, aiding visual inspection of +fits on one line of a command terminal, aiding visual inspection of an ongoing NEB calculation. When a NEB calculation does not converge properly, the supplementary From e1a8a70a6cf9293eaf2dca7b21be561ddda5474c Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 14 Jan 2023 07:06:26 -0500 Subject: [PATCH 22/25] replace individual *verbose* / *terse* keywords with *verbosity* setting --- doc/src/Fortran.rst | 3 +- doc/src/neb.rst | 77 ++++++++++++--------- examples/neb/log.14Jan23.neb.hop1.end.g++.4 | 11 +++ examples/neb/log.14Jan23.neb.hop1.end.g++.8 | 10 +++ examples/neb/log.14Jan23.neb.hop1.g++.4 | 9 +++ examples/neb/log.14Jan23.neb.hop1.g++.8 | 10 +++ examples/neb/log.14Jan23.neb.hop2.g++.4 | 9 +++ examples/neb/log.14Jan23.neb.hop2.g++.8 | 9 +++ examples/neb/log.14Jan23.neb.sivac.g++.4 | 18 +++++ examples/neb/log.14Jan23.neb.sivac.g++.8 | 20 ++++++ examples/neb/log.19Jun17.neb.hop1.end.g++.4 | 11 --- examples/neb/log.19Jun17.neb.hop1.end.g++.8 | 11 --- examples/neb/log.19Jun17.neb.hop1.g++.4 | 9 --- examples/neb/log.19Jun17.neb.hop1.g++.8 | 9 --- examples/neb/log.19Jun17.neb.hop2.g++.4 | 12 ---- examples/neb/log.19Jun17.neb.hop2.g++.8 | 12 ---- examples/neb/log.19Jun17.neb.sivac.g++.4 | 17 ----- examples/neb/log.19Jun17.neb.sivac.g++.8 | 18 ----- src/REPLICA/neb.cpp | 28 +++++--- 19 files changed, 162 insertions(+), 141 deletions(-) create mode 100644 examples/neb/log.14Jan23.neb.hop1.end.g++.4 create mode 100644 examples/neb/log.14Jan23.neb.hop1.end.g++.8 create mode 100644 examples/neb/log.14Jan23.neb.hop1.g++.4 create mode 100644 examples/neb/log.14Jan23.neb.hop1.g++.8 create mode 100644 examples/neb/log.14Jan23.neb.hop2.g++.4 create mode 100644 examples/neb/log.14Jan23.neb.hop2.g++.8 create mode 100644 examples/neb/log.14Jan23.neb.sivac.g++.4 create mode 100644 examples/neb/log.14Jan23.neb.sivac.g++.8 delete mode 100644 examples/neb/log.19Jun17.neb.hop1.end.g++.4 delete mode 100644 examples/neb/log.19Jun17.neb.hop1.end.g++.8 delete mode 100644 examples/neb/log.19Jun17.neb.hop1.g++.4 delete mode 100644 examples/neb/log.19Jun17.neb.hop1.g++.8 delete mode 100644 examples/neb/log.19Jun17.neb.hop2.g++.4 delete mode 100644 examples/neb/log.19Jun17.neb.hop2.g++.8 delete mode 100644 examples/neb/log.19Jun17.neb.sivac.g++.4 delete mode 100644 examples/neb/log.19Jun17.neb.sivac.g++.8 diff --git a/doc/src/Fortran.rst b/doc/src/Fortran.rst index 10947bfef9..d20c9d1df5 100644 --- a/doc/src/Fortran.rst +++ b/doc/src/Fortran.rst @@ -47,7 +47,8 @@ Fortran code in order to uses the Fortran interface. A working example can be found together with equivalent examples in C and C++ in the ``examples/COUPLE/simple`` folder of the LAMMPS distribution. -.. admonitions:: Fortran compiler compatibility +.. admonition:: Fortran compiler compatibility + :class: note A fully Fortran 2003 compatible Fortran compiler is required. This means that currently only GNU Fortran 9 and later are diff --git a/doc/src/neb.rst b/doc/src/neb.rst index 4703cc7f60..a959eee01b 100644 --- a/doc/src/neb.rst +++ b/doc/src/neb.rst @@ -8,7 +8,7 @@ Syntax .. parsed-literal:: - neb etol ftol N1 N2 Nevery file-style arg keyword + neb etol ftol N1 N2 Nevery file-style arg keyword values * etol = stopping tolerance for energy (energy units) * ftol = stopping tolerance for force (force units) @@ -29,7 +29,15 @@ Syntax *none* arg = no argument all replicas assumed to already have their initial coords -* keyword = *verbose* or *terse* +* zero or more keyword/value pairs may be appended +* keyword = *verbosity* + + .. parsed-literal:: + + *verbosity* value = *verbose* or *default* or *terse* + *verbose* = very detailed per-replica output + *default* = some per-replica output + *terse* = only global state output Examples """""""" @@ -304,12 +312,26 @@ and restart files. When running with multiple partitions (each of which is a replica in this case), the print-out to the screen and master log.lammps file -contains a line of output, printed once every *Nevery* timesteps. It -contains the timestep, the maximum force per replica, the maximum -force per atom (in any replica), potential gradients in the initial, -final, and climbing replicas, the forward and backward energy -barriers, the total reaction coordinate (RDT), and the normalized -reaction coordinate and potential energy of each replica. +contains a line of output, printed once every *Nevery* timesteps. The +amount of information printed in this line can be selected with the +*verbosity* keyword. Available options are *terse*, *default*, and +*verbose*. + +With the *terse* setting, it contains the timestep, the maximum force of +a replica, the maximum force per atom (in any replica), potential +gradients in the initial, final, and climbing replicas, the forward and +backward energy barriers, the total reaction coordinate (RDT). + +With the *default* setting, additionally the normalized +reaction coordinate and potential energy of each replica are printed. + +With the *verbose* setting, additional per-replica properties are +printed: the "path angle" (pathangle), the angle between the 3N-length +tangent vector and the 3N-length force vector at image *i* +(angletangrad), the angle between the 3N-length energy gradient vector +of replica *i* and that of replica *i*\ +1 (anglegrad), the norm of the +energy gradient (gradV), the the two-norm of the 3N-length force vector +(RepForce), and the maximum force component of any atom (MaxAtomForce). The "maximum force per replica" is the two-norm of the 3N-length force vector for the atoms in each replica, maximized across replicas, which @@ -332,28 +354,21 @@ the fix neb command. The forward (reverse) energy barrier is the potential energy of the highest replica minus the energy of the first (last) replica. -Supplementary information for all replicas can be printed out by adding -the *verbose* keyword. This information include the following. The -"path angle" (pathangle) for the replica i which is the angle between -the 3N-length vectors :math:`(R_{i-1} - R_i)` and :math:`(R_{i+1} - -R_i)` (where :math:`R_i` is the atomic coordinates of replica *i*). A -"path angle" of 180 indicates that replicas *i*-1, *i* and *i*+1 are -aligned. "angletangrad" is the angle between the 3N-length tangent -vector and the 3N-length force vector at image *i*. The tangent vector -is calculated as in :ref:`(HenkelmanA) ` for all -intermediate replicas and at R2 - R1 and RM - RM-1 for the first and -last replica, respectively. "anglegrad" is the angle between the -3N-length energy gradient vector of replica *i* and that of replica -*i*+1. It is not defined for the final replica and reads nan. gradV is -the norm of the energy gradient of image *i* (:math:`\nabla V`). -ReplicaForce is the two-norm of the 3N-length force vector (including -nudging forces) for replica *i*. MaxAtomForce is the maximum force -component of any atom in replica *i*. - -Alternatively, a restricted print out can be obtained by adding the -*terse* keyword, which omits per-replica information. This typically -fits on one line of a command terminal, aiding visual inspection of -an ongoing NEB calculation. +The "path angle" (pathangle) for the replica i which is the angle +between the 3N-length vectors :math:`(R_{i-1} - R_i)` and +:math:`(R_{i+1} - R_i)` (where :math:`R_i` is the atomic coordinates of +replica *i*). A "path angle" of 180 indicates that replicas *i*\ -1, *i* +and *i*\ +1 are aligned. "angletangrad" is the angle between the +3N-length tangent vector and the 3N-length force vector at image +*i*. The tangent vector is calculated as in :ref:`(HenkelmanA) +` for all intermediate replicas and at R2 - R1 and RM - RM-1 +for the first and last replica, respectively. "anglegrad" is the angle +between the 3N-length energy gradient vector of replica *i* and that of +replica *i*\ +1. It is not defined for the final replica and reads nan. +gradV is the norm of the energy gradient of image *i* (:math:`\nabla +V`). ReplicaForce is the two-norm of the 3N-length force vector +(including nudging forces) for replica *i*. MaxAtomForce is the maximum +force component of any atom in replica *i*. When a NEB calculation does not converge properly, the supplementary information can help understanding what is going wrong. For instance @@ -435,7 +450,7 @@ Related commands Default """"""" -none +*verbosity* = *default* ---------- diff --git a/examples/neb/log.14Jan23.neb.hop1.end.g++.4 b/examples/neb/log.14Jan23.neb.hop1.end.g++.4 new file mode 100644 index 0000000000..b9dcc3cdd9 --- /dev/null +++ b/examples/neb/log.14Jan23.neb.hop1.end.g++.4 @@ -0,0 +1,11 @@ +LAMMPS (22 Dec 2022) +Running on 4 partitions of processors + Step MaxReplicaForce MaxAtomForce GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 RD3 PE3 RD4 PE4 + 0 229.26196 21515.76 2.9774577 4.4127369 233.11559 0.023301843 0.0224626 1.4763579 0 -3.048332 0.33333333 -3.0250302 0.66666667 -3.0291888 1 -3.0474928 + 100 0.11027532 0.0072949206 3.0967938 0.024201563 0.38551033 0.0017583261 0.0021866943 1.7710358 0 -3.0483469 0.31192818 -3.0465886 0.61093022 -3.0466143 1 -3.0487752 + 130 0.09954083 0.0056973977 3.0927626 0.015664388 0.37491833 0.0017573704 0.0021913201 1.7713726 0 -3.048342 0.31428487 -3.0465846 0.61762817 -3.0466296 1 -3.048776 +Climbing replica = 2 + Step MaxReplicaForce MaxAtomForce GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 RD3 PE3 RD4 PE4 + 130 0.37838747 0.12267051 3.0927626 0.015664388 0.37491833 0.0017573704 0.0021913201 1.7713726 0 -3.048342 0.31428487 -3.0465846 0.61762817 -3.0466296 1 -3.048776 + 230 0.14827532 0.015635771 3.1230352 0.0083310795 0.14070676 0.0018353858 0.0022776729 1.769345 0 -3.0483349 0.39530832 -3.0464995 0.64330942 -3.046694 1 -3.0487772 + 279 0.099842536 0.0074295981 3.1400424 0.0068603912 0.095180758 0.0018424455 0.0022860594 1.7685036 0 -3.0483338 0.41228692 -3.0464914 0.65522525 -3.0467282 1 -3.0487775 diff --git a/examples/neb/log.14Jan23.neb.hop1.end.g++.8 b/examples/neb/log.14Jan23.neb.hop1.end.g++.8 new file mode 100644 index 0000000000..bb31e10924 --- /dev/null +++ b/examples/neb/log.14Jan23.neb.hop1.end.g++.8 @@ -0,0 +1,10 @@ +LAMMPS (22 Dec 2022) +Running on 8 partitions of processors + Step MaxReplicaForce MaxAtomForce GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 RD3 PE3 RD4 PE4 RD5 PE5 RD6 PE6 RD7 PE7 RD8 PE8 + 0 353.31221 52709.758 2.9774577 4.4127369 354.16808 0.035499368 0.034660125 1.4763579 0 -3.048332 0.14285714 -3.0450092 0.28571429 -3.0318848 0.42857143 -3.0128327 0.57142857 -3.0157894 0.71428571 -3.0354827 0.85714286 -3.0458819 1 -3.0474928 + 100 0.18371802 0.016897952 3.0181603 0.026618267 0.11111596 0.0018567895 0.0022827846 1.8009397 0 -3.0483491 0.13062498 -3.0471842 0.26462893 -3.0466755 0.40763124 -3.0464923 0.53986455 -3.0465208 0.68351825 -3.0468184 0.81489539 -3.0476023 1 -3.0487751 + 182 0.099603091 0.0039050116 3.0068142 0.010832038 0.074630455 0.001848721 0.0022882994 1.8108868 0 -3.0483371 0.13698322 -3.0471393 0.27599082 -3.0466464 0.41841798 -3.0464884 0.55670019 -3.046546 0.69926219 -3.0469046 0.83686246 -3.0478727 1 -3.0487767 +Climbing replica = 4 + Step MaxReplicaForce MaxAtomForce GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 RD3 PE3 RD4 PE4 RD5 PE5 RD6 PE6 RD7 PE7 RD8 PE8 + 182 0.093479786 0.0046539741 3.0068142 0.010832038 0.074630455 0.001848721 0.0022882994 1.8108868 0 -3.0483371 0.13698322 -3.0471393 0.27599082 -3.0466464 0.41841798 -3.0464884 0.55670019 -3.046546 0.69926219 -3.0469046 0.83686246 -3.0478727 1 -3.0487767 + 183 0.093479786 0.0046539741 3.0068142 0.010832038 0.074630455 0.001848721 0.0022882994 1.8108868 0 -3.0483371 0.13698322 -3.0471393 0.27599082 -3.0466464 0.41841798 -3.0464884 0.55670019 -3.046546 0.69926219 -3.0469046 0.83686246 -3.0478727 1 -3.0487767 diff --git a/examples/neb/log.14Jan23.neb.hop1.g++.4 b/examples/neb/log.14Jan23.neb.hop1.g++.4 new file mode 100644 index 0000000000..8e48a7c9a4 --- /dev/null +++ b/examples/neb/log.14Jan23.neb.hop1.g++.4 @@ -0,0 +1,9 @@ +LAMMPS (22 Dec 2022) +Running on 4 partitions of processors + Step MaxReplicaForce MaxAtomForce GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 RD3 PE3 RD4 PE4 + 0 4327.2753 7542371.4 0.082169072 4.9967651 4514.5424 0.42933428 0.42323635 1.8941131 0 -3.0535948 0.33333333 -2.6242605 0.66666667 -2.7623811 1 -3.0474969 + 87 0.095951502 0.0027794936 0.005588927 0.065110105 0.12467831 0.0071014928 0.0022798007 2.3003372 0 -3.0535967 0.32435271 -3.0473127 0.62805027 -3.0464952 1 -3.048775 +Climbing replica = 3 + Step MaxReplicaForce MaxAtomForce GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 RD3 PE3 RD4 PE4 + 87 0.14137277 0.012340886 0.005588927 0.065110105 0.12467831 0.0071014928 0.0022798007 2.3003372 0 -3.0535967 0.32435271 -3.0473127 0.62805027 -3.0464952 1 -3.048775 + 124 0.099583263 0.0073851506 0.0044220372 0.023873795 0.091308308 0.0071061754 0.0022863931 2.308121 0 -3.0535968 0.32223905 -3.0473329 0.61673898 -3.0464906 1 -3.048777 diff --git a/examples/neb/log.14Jan23.neb.hop1.g++.8 b/examples/neb/log.14Jan23.neb.hop1.g++.8 new file mode 100644 index 0000000000..43b5325f34 --- /dev/null +++ b/examples/neb/log.14Jan23.neb.hop1.g++.8 @@ -0,0 +1,10 @@ +LAMMPS (22 Dec 2022) +Running on 8 partitions of processors + Step MaxReplicaForce MaxAtomForce GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 RD3 PE3 RD4 PE4 RD5 PE5 RD6 PE6 RD7 PE7 RD8 PE8 + 0 12172.47 48315818 0.082169072 4.9967651 12254.987 1.1034547 1.0973568 1.8941131 0 -3.0535948 0.14285714 -3.0398884 0.28571429 -2.8426874 0.42857143 -1.9501401 0.57142857 -2.1439119 0.71428571 -2.9187789 0.85714286 -3.0422538 1 -3.0474969 + 100 0.14433533 0.0056834028 0.0051714625 0.047051976 0.11774619 0.0071018066 0.002280994 2.3772214 0 -3.0535967 0.12651739 -3.0512933 0.26516133 -3.0480936 0.41814034 -3.0467867 0.55534109 -3.0464949 0.69779433 -3.0465932 0.83399724 -3.0472934 1 -3.0487759 + 147 0.098667494 0.0031400895 0.0040851637 0.017209593 0.096695061 0.007103925 0.002284434 2.3892288 0 -3.0535968 0.13266908 -3.051146 0.27206276 -3.0479976 0.42113052 -3.046774 0.55968697 -3.0464929 0.7017833 -3.0466074 0.8397451 -3.0473745 1 -3.0487774 +Climbing replica = 5 + Step MaxReplicaForce MaxAtomForce GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 RD3 PE3 RD4 PE4 RD5 PE5 RD6 PE6 RD7 PE7 RD8 PE8 + 147 0.10065405 0.0077673218 0.0040851637 0.017209593 0.096695061 0.007103925 0.002284434 2.3892288 0 -3.0535968 0.13266908 -3.051146 0.27206276 -3.0479976 0.42113052 -3.046774 0.55968697 -3.0464929 0.7017833 -3.0466074 0.8397451 -3.0473745 1 -3.0487774 + 150 0.099859382 0.0075603521 0.0040615957 0.016809825 0.096007238 0.0071040699 0.0022845963 2.3894645 0 -3.0535969 0.13277764 -3.0511437 0.27216411 -3.0479964 0.42112067 -3.0467741 0.55997348 -3.0464928 0.70181216 -3.0466076 0.83985415 -3.0473761 1 -3.0487774 diff --git a/examples/neb/log.14Jan23.neb.hop2.g++.4 b/examples/neb/log.14Jan23.neb.hop2.g++.4 new file mode 100644 index 0000000000..c608a35ca8 --- /dev/null +++ b/examples/neb/log.14Jan23.neb.hop2.g++.4 @@ -0,0 +1,9 @@ +LAMMPS (22 Dec 2022) +Running on 4 partitions of processors + Step MaxReplicaForce MaxAtomForce GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 RD3 PE3 RD4 PE4 + 0 14.104748 108.56876 0.1227071 4.999238 8.2087606 0.0018276223 0.00064050211 0.98401186 0 -3.0514921 0.33333333 -3.0496673 0.66666667 -3.0496645 1 -3.050305 + 46 0.049691778 0.001449144 0.0041836599 0.0070455894 0.94610709 0.0014279446 0.0014276971 1.1698782 0 -3.0514942 0.3026746 -3.0502925 0.63457326 -3.0500662 1 -3.0514939 +Climbing replica = 3 + Step MaxReplicaForce MaxAtomForce GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 RD3 PE3 RD4 PE4 + 46 0.94646532 0.8692782 0.0041836599 0.0070455894 0.94610709 0.0014279446 0.0014276971 1.1698782 0 -3.0514942 0.3026746 -3.0502925 0.63457326 -3.0500662 1 -3.0514939 + 88 0.042941386 0.00046017358 0.0021155242 0.0031463438 0.0095074767 0.0016016953 0.0016016258 1.1725045 0 -3.0514943 0.26292929 -3.0504854 0.50610293 -3.0498926 1 -3.0514943 diff --git a/examples/neb/log.14Jan23.neb.hop2.g++.8 b/examples/neb/log.14Jan23.neb.hop2.g++.8 new file mode 100644 index 0000000000..65f071d8b4 --- /dev/null +++ b/examples/neb/log.14Jan23.neb.hop2.g++.8 @@ -0,0 +1,9 @@ +LAMMPS (22 Dec 2022) +Running on 8 partitions of processors + Step MaxReplicaForce MaxAtomForce GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 RD3 PE3 RD4 PE4 RD5 PE5 RD6 PE6 RD7 PE7 RD8 PE8 + 0 15.805883 137.15721 0.1227071 4.999238 15.828539 0.0021830032 0.00099588309 0.98401186 0 -3.0514921 0.14285714 -3.0509383 0.28571429 -3.049949 0.42857143 -3.0493091 0.57142857 -3.0493602 0.71428571 -3.0498345 0.85714286 -3.0502192 1 -3.050305 + 73 0.049714315 0.0011853232 0.0036385263 0.0057874304 0.079285231 0.0016005788 0.001600399 1.2115403 0 -3.0514942 0.13043211 -3.0512538 0.25481496 -3.0505642 0.37535798 -3.0500676 0.49844851 -3.0498936 0.63571599 -3.0500698 0.80272589 -3.0508376 1 -3.051494 +Climbing replica = 5 + Step MaxReplicaForce MaxAtomForce GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 RD3 PE3 RD4 PE4 RD5 PE5 RD6 PE6 RD7 PE7 RD8 PE8 + 73 0.079737673 0.0060506501 0.0036385263 0.0057874304 0.079285231 0.0016005788 0.001600399 1.2115403 0 -3.0514942 0.13043211 -3.0512538 0.25481496 -3.0505642 0.37535798 -3.0500676 0.49844851 -3.0498936 0.63571599 -3.0500698 0.80272589 -3.0508376 1 -3.051494 + 82 0.046923812 0.00207513 0.0033938491 0.0054082765 0.046102157 0.0016015266 0.0016013694 1.2125153 0 -3.0514942 0.13012408 -3.0512565 0.25462869 -3.0505663 0.37576663 -3.0500667 0.50363184 -3.0498927 0.63838032 -3.0500776 0.80516707 -3.0508541 1 -3.0514941 diff --git a/examples/neb/log.14Jan23.neb.sivac.g++.4 b/examples/neb/log.14Jan23.neb.sivac.g++.4 new file mode 100644 index 0000000000..b2e5250229 --- /dev/null +++ b/examples/neb/log.14Jan23.neb.sivac.g++.4 @@ -0,0 +1,18 @@ +LAMMPS (22 Dec 2022) +Running on 4 partitions of processors + Step MaxReplicaForce MaxAtomForce GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 RD3 PE3 RD4 PE4 + 0 7.5525391 2.671788 0.16683659 7.5525391 7.5525391 1.5383951 0 1.6207355 0 -2213.3343 0.33333333 -2212.7428 0.66666667 -2212.2247 1 -2211.7959 + 10 0.24005275 0.0013324036 0.036483049 0.24005275 0.68351722 0.42916118 0.41794425 1.6989349 0 -2213.3365 0.32909183 -2212.9587 0.65386736 -2212.9073 1 -2213.3253 + 20 0.07940898 0.00026889621 0.024706844 0.07940898 0.71637784 0.41387872 0.41157886 1.7343662 0 -2213.3369 0.32478734 -2212.9621 0.65348766 -2212.923 1 -2213.3346 + 30 0.094973706 6.994258e-05 0.015145947 0.035267404 0.7535772 0.40072717 0.40024605 1.7504612 0 -2213.3372 0.32705584 -2212.9584 0.65894506 -2212.9365 1 -2213.3367 + 40 0.027727472 1.9827557e-05 0.011618173 0.022562656 0.76133752 0.39614635 0.39591731 1.7547519 0 -2213.3373 0.32873163 -2212.9562 0.66124255 -2212.9411 1 -2213.337 + 50 0.019429351 9.0662903e-06 0.0087135562 0.015391975 0.76952681 0.39274846 0.3926388 1.7578616 0 -2213.3373 0.33022595 -2212.9543 0.66307279 -2212.9446 1 -2213.3372 + 60 0.018992557 2.6317454e-06 0.005342608 0.0086165759 0.777596 0.38936859 0.38933361 1.7610433 0 -2213.3374 0.33187549 -2212.9523 0.66497619 -2212.948 1 -2213.3373 + 63 0.0097571737 1.6189866e-06 0.0047788465 0.0076143824 0.77864965 0.38888882 0.3888615 1.7615283 0 -2213.3374 0.33212054 -2212.952 0.66525325 -2212.9485 1 -2213.3373 +Climbing replica = 3 + Step MaxReplicaForce MaxAtomForce GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 RD3 PE3 RD4 PE4 + 63 0.77864965 0.096633568 0.0047788465 0.0076143824 0.77864965 0.38888882 0.3888615 1.7615283 0 -2213.3374 0.33212054 -2212.952 0.66525325 -2212.9485 1 -2213.3373 + 73 0.098863403 0.0011310488 0.0027890076 0.0042744972 0.036283915 0.5102483 0.51023975 1.76072 0 -2213.3374 0.2757098 -2213.0417 0.50430234 -2212.8271 1 -2213.3374 + 83 0.031748825 0.00015683838 0.0020809514 0.0031533485 0.0099572603 0.51014575 0.51014111 1.7602531 0 -2213.3374 0.26034606 -2213.0674 0.50354802 -2212.8272 1 -2213.3374 + 93 0.011560904 6.778749e-06 0.0013983005 0.0020949841 0.0052524482 0.51011026 0.51010823 1.7601164 0 -2213.3374 0.2537964 -2213.0786 0.50388122 -2212.8273 1 -2213.3374 + 96 0.0080547159 4.4226058e-06 0.001318412 0.0019726613 0.0048040798 0.510108 0.5101062 1.7601158 0 -2213.3374 0.25347618 -2213.0792 0.50393582 -2212.8273 1 -2213.3374 diff --git a/examples/neb/log.14Jan23.neb.sivac.g++.8 b/examples/neb/log.14Jan23.neb.sivac.g++.8 new file mode 100644 index 0000000000..b699c6f14c --- /dev/null +++ b/examples/neb/log.14Jan23.neb.sivac.g++.8 @@ -0,0 +1,20 @@ +LAMMPS (22 Dec 2022) +Running on 8 partitions of processors + Step MaxReplicaForce MaxAtomForce GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 RD3 PE3 RD4 PE4 RD5 PE5 RD6 PE6 RD7 PE7 RD8 PE8 + 0 7.5525391 2.671788 0.16683659 7.5525391 7.5525391 1.5383951 0 1.6207355 0 -2213.3343 0.14285714 -2213.1848 0.28571429 -2212.8577 0.42857143 -2212.5353 0.57142857 -2212.3131 0.71428571 -2212.192 0.85714286 -2212.0797 1 -2211.7959 + 10 0.45237221 0.0038899054 0.043713415 0.45237221 0.33984237 0.50434364 0.47974444 1.681593 0 -2213.3362 0.145972 -2213.2393 0.28963267 -2213.0224 0.43234899 -2212.8512 0.57436617 -2212.8319 0.71564291 -2212.9699 0.85294239 -2213.1766 1 -2213.3116 + 20 0.17466055 0.00038449492 0.027733064 0.16054263 0.29142661 0.49877454 0.49547267 1.7524163 0 -2213.3369 0.14387205 -2213.2456 0.28439667 -2213.0291 0.42233932 -2212.8573 0.5588361 -2212.8381 0.69513708 -2212.9781 0.83702043 -2213.2011 1 -2213.3336 + 30 0.045422642 0.00010721083 0.01818213 0.045422642 0.28112975 0.49759387 0.49675173 1.7746784 0 -2213.3371 0.14290957 -2213.2488 0.28281709 -2213.0326 0.41989281 -2212.859 0.55582211 -2212.8395 0.69368801 -2212.9853 0.83973606 -2213.2161 1 -2213.3363 + 40 0.038021578 6.0627173e-05 0.013142053 0.027303941 0.28657705 0.49667648 0.49635337 1.7848045 0 -2213.3372 0.1420496 -2213.2509 0.28169834 -2213.0351 0.41895631 -2212.8598 0.55576002 -2212.8406 0.69554601 -2212.9918 0.84293215 -2213.2247 1 -2213.3369 + 50 0.023478387 2.5440651e-05 0.0094370881 0.017117841 0.2948129 0.4956627 0.49552819 1.7913648 0 -2213.3373 0.14130551 -2213.2524 0.28080367 -2213.0371 0.41862494 -2212.8601 0.55666336 -2212.8417 0.69796419 -2212.9977 0.84561062 -2213.2308 1 -2213.3372 + 60 0.015568348 1.3385057e-05 0.0065356733 0.010880239 0.30476692 0.49454143 0.49448589 1.7961707 0 -2213.3374 0.14069048 -2213.2537 0.28018046 -2213.0385 0.41875119 -2212.86 0.55806463 -2212.8428 0.70047903 -2213.0032 0.84789117 -2213.2356 1 -2213.3373 + 70 0.011117036 8.589416e-06 0.0046665959 0.0074340441 0.31343723 0.4935777 0.49355166 1.7991441 0 -2213.3374 0.14033265 -2213.2544 0.27992696 -2213.0392 0.41915513 -2212.8596 0.55940642 -2212.8438 0.7024513 -2213.0073 0.8494878 -2213.2387 1 -2213.3373 + 77 0.0095032996 6.4872648e-06 0.003778822 0.0059131255 0.31843187 0.49301961 0.49300316 1.8005356 0 -2213.3374 0.14019726 -2213.2547 0.27989276 -2213.0394 0.41949002 -2212.8594 0.56021116 -2212.8444 0.70351827 -2213.0094 0.85030233 -2213.2402 1 -2213.3374 +Climbing replica = 5 + Step MaxReplicaForce MaxAtomForce GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 RD3 PE3 RD4 PE4 RD5 PE5 RD6 PE6 RD7 PE7 RD8 PE8 + 77 0.31843187 0.016533023 0.003778822 0.0059131255 0.31843187 0.49301961 0.49300316 1.8005356 0 -2213.3374 0.14019726 -2213.2547 0.27989276 -2213.0394 0.41949002 -2212.8594 0.56021116 -2212.8444 0.70351827 -2213.0094 0.85030233 -2213.2402 1 -2213.3374 + 87 0.068521294 0.00041696849 0.002590902 0.0039658354 0.024993269 0.51001428 0.51000692 1.8017659 0 -2213.3374 0.13988527 -2213.2554 0.27815879 -2213.0427 0.4060317 -2212.8709 0.50551576 -2212.8274 0.69112296 -2212.99 0.84991798 -2213.2398 1 -2213.3374 + 97 0.031487615 0.00015286198 0.0020335677 0.0030841188 0.0070128682 0.51009653 0.51009209 1.8020354 0 -2213.3374 0.13898618 -2213.2569 0.27382107 -2213.0504 0.39686059 -2212.8798 0.50256203 -2212.8273 0.68247072 -2212.9764 0.84621774 -2213.2346 1 -2213.3374 + 107 0.016814587 4.4986184e-05 0.0013659976 0.0020482331 0.0024996635 0.51009833 0.51009639 1.8028763 0 -2213.3374 0.13616152 -2213.2608 0.26706798 -2213.0624 0.38875801 -2212.8885 0.50194006 -2212.8273 0.67498881 -2212.9648 0.84072701 -2213.2269 1 -2213.3374 + 117 0.013420011 2.244148e-05 0.00091637825 0.0013624439 0.0016781664 0.51009641 0.51009556 1.8036672 0 -2213.3374 0.13322116 -2213.2647 0.26208966 -2213.0712 0.38454183 -2212.8933 0.50207828 -2212.8273 0.67150999 -2212.9596 0.83745606 -2213.2222 1 -2213.3374 + 127 0.0091080864 1.3208992e-05 0.00070311367 0.0010407514 0.0013146017 0.5100957 0.5100952 1.804094 0 -2213.3374 0.13160053 -2213.2667 0.25966338 -2213.0755 0.38274006 -2212.8954 0.50217662 -2212.8273 0.6702885 -2212.9577 0.83622287 -2213.2204 1 -2213.3374 diff --git a/examples/neb/log.19Jun17.neb.hop1.end.g++.4 b/examples/neb/log.19Jun17.neb.hop1.end.g++.4 deleted file mode 100644 index 4878b86566..0000000000 --- a/examples/neb/log.19Jun17.neb.hop1.end.g++.4 +++ /dev/null @@ -1,11 +0,0 @@ -LAMMPS (19 May 2017) -Running on 4 partitions of processors -Step MaxReplicaForce MaxAtomForce GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 ... RDN PEN -0 229.26196 146.68251 2.9774577 4.4127369 233.11559 0.023301843 0.0224626 1.4763579 0 -3.048332 0.33333333 -3.0250302 0.66666667 -3.0291888 1 -3.0474928 -100 0.11027532 0.085410308 3.0967938 0.024201563 0.38551033 0.0017583261 0.0021866943 1.7710358 0 -3.0483469 0.31192818 -3.0465886 0.61093022 -3.0466143 1 -3.0487752 -130 0.09954083 0.075481108 3.0927626 0.015664388 0.37491833 0.0017573704 0.0021913201 1.7713726 0 -3.048342 0.31428487 -3.0465846 0.61762817 -3.0466296 1 -3.048776 -Climbing replica = 2 -Step MaxReplicaForce MaxAtomForce GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 ... RDN PEN -130 0.37838747 0.3502435 3.0927626 0.015664388 0.37491833 0.0017573704 0.0021913201 1.7713726 0 -3.048342 0.31428487 -3.0465846 0.61762817 -3.0466296 1 -3.048776 -230 0.22757286 0.12027481 3.1250243 0.0081260569 0.14019507 0.0018364585 0.002278918 1.76926 0 -3.0483347 0.39730698 -3.0464983 0.64450769 -3.0466973 1 -3.0487772 -278 0.096184498 0.085088496 3.1405655 0.0068164307 0.093861113 0.0018426056 0.002286256 1.7684765 0 -3.0483338 0.41277997 -3.0464912 0.65562984 -3.0467294 1 -3.0487775 diff --git a/examples/neb/log.19Jun17.neb.hop1.end.g++.8 b/examples/neb/log.19Jun17.neb.hop1.end.g++.8 deleted file mode 100644 index 62344b3da5..0000000000 --- a/examples/neb/log.19Jun17.neb.hop1.end.g++.8 +++ /dev/null @@ -1,11 +0,0 @@ -LAMMPS (19 May 2017) -Running on 4 partitions of processors -Step MaxReplicaForce MaxAtomForce GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 ... RDN PEN -0 229.26196 146.68251 2.9774577 4.4127369 233.11559 0.023301843 0.0224626 1.4763579 0 -3.048332 0.33333333 -3.0250302 0.66666667 -3.0291888 1 -3.0474928 -100 0.11375359 0.085350745 3.0966418 0.0236765 0.38531777 0.0017582606 0.0021868783 1.7710738 0 -3.0483467 0.31201141 -3.0465884 0.61117406 -3.0466149 1 -3.0487753 -119 0.09996986 0.078639268 3.0937691 0.017444108 0.3780308 0.0017574935 0.0021899317 1.7713574 0 -3.0483433 0.31354192 -3.0465858 0.61555533 -3.0466249 1 -3.0487758 -Climbing replica = 2 -Step MaxReplicaForce MaxAtomForce GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 ... RDN PEN -119 0.3793192 0.35281863 3.0937691 0.017444108 0.3780308 0.0017574935 0.0021899317 1.7713574 0 -3.0483433 0.31354192 -3.0465858 0.61555533 -3.0466249 1 -3.0487758 -219 0.20159133 0.12247026 3.1244061 0.0085896057 0.13938632 0.0018362816 0.0022783681 1.7693295 0 -3.048335 0.39646633 -3.0464988 0.64277703 -3.0466925 1 -3.0487771 -266 0.099868725 0.086180598 3.1401661 0.0070922949 0.095128081 0.001842608 0.002286044 1.7685191 0 -3.048334 0.41231024 -3.0464914 0.65425179 -3.0467252 1 -3.0487774 diff --git a/examples/neb/log.19Jun17.neb.hop1.g++.4 b/examples/neb/log.19Jun17.neb.hop1.g++.4 deleted file mode 100644 index e2984c031c..0000000000 --- a/examples/neb/log.19Jun17.neb.hop1.g++.4 +++ /dev/null @@ -1,9 +0,0 @@ -LAMMPS (19 May 2017) -Running on 4 partitions of processors -Step MaxReplicaForce MaxAtomForce GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 ... RDN PEN -0 4327.2753 2746.3378 0.082169072 4.9967651 4514.5424 0.42933428 0.42323635 1.8941131 0 -3.0535948 0.33333333 -2.6242605 0.66666667 -2.7623811 1 -3.0474969 -87 0.095951502 0.052720903 0.005588927 0.065110105 0.12467831 0.0071014928 0.0022798007 2.3003372 0 -3.0535967 0.32435271 -3.0473127 0.62805027 -3.0464952 1 -3.048775 -Climbing replica = 3 -Step MaxReplicaForce MaxAtomForce GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 ... RDN PEN -87 0.14137277 0.11108954 0.005588927 0.065110105 0.12467831 0.0071014928 0.0022798007 2.3003372 0 -3.0535967 0.32435271 -3.0473127 0.62805027 -3.0464952 1 -3.048775 -124 0.099583263 0.085936899 0.0044220372 0.023873795 0.091308308 0.0071061754 0.0022863931 2.308121 0 -3.0535968 0.32223905 -3.0473329 0.61673898 -3.0464906 1 -3.048777 diff --git a/examples/neb/log.19Jun17.neb.hop1.g++.8 b/examples/neb/log.19Jun17.neb.hop1.g++.8 deleted file mode 100644 index d1be1284fa..0000000000 --- a/examples/neb/log.19Jun17.neb.hop1.g++.8 +++ /dev/null @@ -1,9 +0,0 @@ -LAMMPS (19 May 2017) -Running on 4 partitions of processors -Step MaxReplicaForce MaxAtomForce GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 ... RDN PEN -0 4327.2753 2746.3378 0.082169072 4.9967651 4514.5424 0.42933428 0.42323635 1.8941131 0 -3.0535948 0.33333333 -2.6242605 0.66666667 -2.7623811 1 -3.0474969 -87 0.095951792 0.052720902 0.0055889267 0.065110091 0.12467831 0.0071014928 0.0022798007 2.3003372 0 -3.0535967 0.32435271 -3.0473127 0.62805027 -3.0464952 1 -3.048775 -Climbing replica = 3 -Step MaxReplicaForce MaxAtomForce GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 ... RDN PEN -87 0.14137297 0.11108954 0.0055889267 0.065110091 0.12467831 0.0071014928 0.0022798007 2.3003372 0 -3.0535967 0.32435271 -3.0473127 0.62805027 -3.0464952 1 -3.048775 -124 0.099582186 0.08593683 0.0044220345 0.023873731 0.091308197 0.0071061754 0.0022863931 2.3081211 0 -3.0535968 0.32223904 -3.0473329 0.61673896 -3.0464906 1 -3.048777 diff --git a/examples/neb/log.19Jun17.neb.hop2.g++.4 b/examples/neb/log.19Jun17.neb.hop2.g++.4 deleted file mode 100644 index c6b6cbe2ce..0000000000 --- a/examples/neb/log.19Jun17.neb.hop2.g++.4 +++ /dev/null @@ -1,12 +0,0 @@ -LAMMPS (19 May 2017) -Running on 4 partitions of processors -Step MaxReplicaForce MaxAtomForce GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 ... RDN PEN -0 14.104748 10.419633 0.1227071 4.999238 8.2087606 0.0018276223 0.00064050211 0.98401186 0 -3.0514921 0.33333333 -3.0496673 0.66666667 -3.0496645 1 -3.050305 -100 0.24646695 0.10792196 0.0077146918 0.058733261 0.63504706 0.001516756 0.0015151635 1.165391 0 -3.0514939 0.2890334 -3.0503533 0.59718494 -3.0499771 1 -3.0514923 -200 0.061777741 0.050288749 0.0047486883 0.0095236035 0.88698597 0.0014465772 0.0014462528 1.1692938 0 -3.0514941 0.29975094 -3.0503052 0.62768286 -3.0500476 1 -3.0514938 -261 0.048699591 0.038138604 0.0040083594 0.0074854409 0.95722712 0.0014243579 0.0014241377 1.1696848 0 -3.0514942 0.30525481 -3.0502812 0.6357998 -3.0500698 1 -3.051494 -Climbing replica = 3 -Step MaxReplicaForce MaxAtomForce GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 ... RDN PEN -261 0.95753855 0.94297239 0.0040083594 0.0074854409 0.95722712 0.0014243579 0.0014241377 1.1696848 0 -3.0514942 0.30525481 -3.0502812 0.6357998 -3.0500698 1 -3.051494 -361 0.072509627 0.06580631 0.0027545765 0.0044749366 0.016746483 0.0016018879 0.0016017805 1.1704611 0 -3.0514943 0.28176307 -3.0503855 0.50355454 -3.0498924 1 -3.0514942 -381 0.04884836 0.040787876 0.0023445904 0.0035162935 0.017959209 0.0016017716 0.0016016898 1.1713862 0 -3.0514943 0.27120138 -3.0504399 0.50428218 -3.0498925 1 -3.0514942 diff --git a/examples/neb/log.19Jun17.neb.hop2.g++.8 b/examples/neb/log.19Jun17.neb.hop2.g++.8 deleted file mode 100644 index c6b6cbe2ce..0000000000 --- a/examples/neb/log.19Jun17.neb.hop2.g++.8 +++ /dev/null @@ -1,12 +0,0 @@ -LAMMPS (19 May 2017) -Running on 4 partitions of processors -Step MaxReplicaForce MaxAtomForce GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 ... RDN PEN -0 14.104748 10.419633 0.1227071 4.999238 8.2087606 0.0018276223 0.00064050211 0.98401186 0 -3.0514921 0.33333333 -3.0496673 0.66666667 -3.0496645 1 -3.050305 -100 0.24646695 0.10792196 0.0077146918 0.058733261 0.63504706 0.001516756 0.0015151635 1.165391 0 -3.0514939 0.2890334 -3.0503533 0.59718494 -3.0499771 1 -3.0514923 -200 0.061777741 0.050288749 0.0047486883 0.0095236035 0.88698597 0.0014465772 0.0014462528 1.1692938 0 -3.0514941 0.29975094 -3.0503052 0.62768286 -3.0500476 1 -3.0514938 -261 0.048699591 0.038138604 0.0040083594 0.0074854409 0.95722712 0.0014243579 0.0014241377 1.1696848 0 -3.0514942 0.30525481 -3.0502812 0.6357998 -3.0500698 1 -3.051494 -Climbing replica = 3 -Step MaxReplicaForce MaxAtomForce GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 ... RDN PEN -261 0.95753855 0.94297239 0.0040083594 0.0074854409 0.95722712 0.0014243579 0.0014241377 1.1696848 0 -3.0514942 0.30525481 -3.0502812 0.6357998 -3.0500698 1 -3.051494 -361 0.072509627 0.06580631 0.0027545765 0.0044749366 0.016746483 0.0016018879 0.0016017805 1.1704611 0 -3.0514943 0.28176307 -3.0503855 0.50355454 -3.0498924 1 -3.0514942 -381 0.04884836 0.040787876 0.0023445904 0.0035162935 0.017959209 0.0016017716 0.0016016898 1.1713862 0 -3.0514943 0.27120138 -3.0504399 0.50428218 -3.0498925 1 -3.0514942 diff --git a/examples/neb/log.19Jun17.neb.sivac.g++.4 b/examples/neb/log.19Jun17.neb.sivac.g++.4 deleted file mode 100644 index 0d9880ca81..0000000000 --- a/examples/neb/log.19Jun17.neb.sivac.g++.4 +++ /dev/null @@ -1,17 +0,0 @@ -LAMMPS (19 May 2017) -Running on 4 partitions of processors -Step MaxReplicaForce MaxAtomForce GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 ... RDN PEN -0 7.5525391 1.6345605 0.16683659 7.5525391 7.5525391 1.5383951 0 1.6207355 0 -2213.3343 0.33333333 -2212.7428 0.66666667 -2212.2247 1 -2211.7959 -10 0.24005275 0.036502104 0.036483049 0.24005275 0.68351722 0.42916118 0.41794425 1.6989349 0 -2213.3365 0.32909183 -2212.9587 0.65386736 -2212.9073 1 -2213.3253 -20 0.07940898 0.016398055 0.024706844 0.07940898 0.71637784 0.41387872 0.41157886 1.7343662 0 -2213.3369 0.32478734 -2212.9621 0.65348766 -2212.923 1 -2213.3346 -30 0.094973707 0.0083631681 0.015145947 0.035267404 0.7535772 0.40072717 0.40024605 1.7504612 0 -2213.3372 0.32705584 -2212.9584 0.65894506 -2212.9365 1 -2213.3367 -40 0.027727472 0.0044528145 0.011618173 0.022562656 0.76133752 0.39614635 0.39591731 1.7547519 0 -2213.3373 0.32873163 -2212.9562 0.66124255 -2212.9411 1 -2213.337 -50 0.019429348 0.0030110281 0.0087135563 0.015391975 0.76952681 0.39274846 0.3926388 1.7578616 0 -2213.3373 0.33022595 -2212.9543 0.66307279 -2212.9446 1 -2213.3372 -60 0.019009471 0.0016234562 0.0053426307 0.0086166186 0.77759617 0.38936861 0.38933364 1.7610433 0 -2213.3374 0.33187548 -2212.9523 0.66497617 -2212.948 1 -2213.3373 -63 0.0097365134 0.0012734598 0.004777604 0.0076121987 0.77865149 0.38888778 0.38886047 1.7615294 0 -2213.3374 0.33212107 -2212.952 0.66525385 -2212.9485 1 -2213.3373 -Climbing replica = 3 -Step MaxReplicaForce MaxAtomForce GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 ... RDN PEN -63 0.77865149 0.31085821 0.004777604 0.0076121987 0.77865149 0.38888778 0.38886047 1.7615294 0 -2213.3374 0.33212107 -2212.952 0.66525385 -2212.9485 1 -2213.3373 -73 0.098175496 0.033609035 0.0027886955 0.0042742148 0.036594003 0.51024838 0.51023983 1.7607181 0 -2213.3374 0.27574151 -2213.0416 0.50432348 -2212.8271 1 -2213.3374 -83 0.03341862 0.012760857 0.0020868177 0.0031625649 0.010189924 0.51014634 0.51014168 1.7602562 0 -2213.3374 0.26045338 -2213.0672 0.50355193 -2212.8272 1 -2213.3374 -93 0.0097374358 0.0028416114 0.0014003718 0.0020986584 0.0053485291 0.51011052 0.51010848 1.7601202 0 -2213.3374 0.25397887 -2213.0783 0.50388111 -2212.8273 1 -2213.3374 diff --git a/examples/neb/log.19Jun17.neb.sivac.g++.8 b/examples/neb/log.19Jun17.neb.sivac.g++.8 deleted file mode 100644 index 260eb9e18b..0000000000 --- a/examples/neb/log.19Jun17.neb.sivac.g++.8 +++ /dev/null @@ -1,18 +0,0 @@ -LAMMPS (19 May 2017) -Running on 4 partitions of processors -Step MaxReplicaForce MaxAtomForce GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 ... RDN PEN -0 7.5525391 1.6345605 0.16683659 7.5525391 7.5525391 1.5383951 0 1.6207355 0 -2213.3343 0.33333333 -2212.7428 0.66666667 -2212.2247 1 -2211.7959 -10 0.24005275 0.036502104 0.036483049 0.24005275 0.68351722 0.42916118 0.41794425 1.6989349 0 -2213.3365 0.32909183 -2212.9587 0.65386736 -2212.9073 1 -2213.3253 -20 0.07940898 0.016398055 0.024706844 0.07940898 0.71637784 0.41387872 0.41157886 1.7343662 0 -2213.3369 0.32478734 -2212.9621 0.65348766 -2212.923 1 -2213.3346 -30 0.094973708 0.0083631681 0.015145947 0.035267404 0.7535772 0.40072717 0.40024605 1.7504612 0 -2213.3372 0.32705584 -2212.9584 0.65894506 -2212.9365 1 -2213.3367 -40 0.027727472 0.0044528144 0.011618173 0.022562656 0.76133752 0.39614635 0.39591731 1.7547519 0 -2213.3373 0.32873163 -2212.9562 0.66124255 -2212.9411 1 -2213.337 -50 0.019429341 0.0030110281 0.0087135565 0.015391975 0.7695268 0.39274846 0.3926388 1.7578616 0 -2213.3373 0.33022595 -2212.9543 0.66307279 -2212.9446 1 -2213.3372 -60 0.019048963 0.0016262345 0.0053426844 0.0086167196 0.77759655 0.38936867 0.3893337 1.7610433 0 -2213.3374 0.33187545 -2212.9523 0.66497615 -2212.948 1 -2213.3373 -63 0.0097037048 0.0012761841 0.0047749367 0.0076075138 0.77865545 0.38888554 0.38885827 1.7615318 0 -2213.3374 0.33212221 -2212.952 0.66525512 -2212.9485 1 -2213.3373 -Climbing replica = 3 -Step MaxReplicaForce MaxAtomForce GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 ... RDN PEN -63 0.77865545 0.3108551 0.0047749367 0.0076075138 0.77865545 0.38888554 0.38885827 1.7615318 0 -2213.3374 0.33212221 -2212.952 0.66525512 -2212.9485 1 -2213.3373 -73 0.098595989 0.033659485 0.0027927196 0.0042813387 0.038224344 0.51024759 0.51023901 1.7607156 0 -2213.3374 0.27595612 -2213.0413 0.50453988 -2212.8271 1 -2213.3374 -83 0.033344977 0.012868685 0.0020880608 0.0031645847 0.010250413 0.51014677 0.5101421 1.7602601 0 -2213.3374 0.26053624 -2213.067 0.50358775 -2212.8272 1 -2213.3374 -93 0.013254873 0.0038176141 0.0014928226 0.0022407967 0.0058577818 0.51011371 0.51011138 1.7601272 0 -2213.3374 0.25452741 -2213.0774 0.50382161 -2212.8273 1 -2213.3374 -95 0.0099964951 0.0031053214 0.0014131665 0.0021184362 0.0053683638 0.51011105 0.51010897 1.7601232 0 -2213.3374 0.2540975 -2213.0781 0.50387313 -2212.8273 1 -2213.3374 diff --git a/src/REPLICA/neb.cpp b/src/REPLICA/neb.cpp index 2e89e0047e..26ceb3a5fd 100644 --- a/src/REPLICA/neb.cpp +++ b/src/REPLICA/neb.cpp @@ -40,7 +40,7 @@ using namespace MathConst; #define CHUNK 1024 #define ATTRIBUTE_PERLINE 4 -enum { NORMAL, TERSE, VERBOSE }; +enum { DEFAULT, TERSE, VERBOSE }; /* ---------------------------------------------------------------------- */ @@ -149,16 +149,18 @@ void NEB::command(int narg, char **arg) // process file-style setting to setup initial configs for all replicas int iarg = 5; int filecmd = 0; - print_mode = NORMAL; // normal + print_mode = DEFAULT; while (iarg < narg) { if (strcmp(arg[iarg], "final") == 0) { - if (iarg + 2 > narg) error->universe_all(FLERR, "Illegal NEB command: missing arguments"); + if (iarg + 2 > narg) + error->universe_all(FLERR, "Illegal NEB final command: missing arguments"); inpfile = arg[iarg + 1]; readfile(inpfile, 0); filecmd = 1; iarg += 2; } else if (strcmp(arg[iarg], "each") == 0) { - if (iarg + 2 > narg) error->universe_all(FLERR, "Illegal NEB command: missing arguments"); + if (iarg + 2 > narg) + error->universe_all(FLERR, "Illegal NEB each command: missing arguments"); inpfile = arg[iarg + 1]; readfile(inpfile, 1); filecmd = 1; @@ -166,12 +168,18 @@ void NEB::command(int narg, char **arg) } else if (strcmp(arg[iarg], "none") == 0) { filecmd = 1; ++iarg; - } else if (strcmp(arg[iarg], "verbose") == 0) { - print_mode = VERBOSE; - ++iarg; - } else if (strcmp(arg[iarg], "terse") == 0) { - print_mode = TERSE; - ++iarg; + } else if (strcmp(arg[iarg], "verbosity") == 0) { + if (iarg + 2 > narg) + error->universe_all(FLERR, "Illegal NEB verbosity command: missing arguments"); + if (strcmp(arg[iarg+1], "verbose") == 0) + print_mode = VERBOSE; + else if (strcmp(arg[iarg+1], "default") == 0) + print_mode = DEFAULT; + else if (strcmp(arg[iarg+1], "terse") == 0) + print_mode = TERSE; + else + error->universe_all(FLERR, fmt::format("Unknown NEB verbosity option {}", arg[iarg+1])); + iarg += 2; } else error->universe_all(FLERR, fmt::format("Unknown NEB command keyword: {}", arg[iarg])); } From 07da78dfe818af41ef35568378d8bb6b49d735ec Mon Sep 17 00:00:00 2001 From: Thomas Swinburne Date: Mon, 16 Jan 2023 14:12:52 +0100 Subject: [PATCH 23/25] Documentation update after suggestions of @athomps --- doc/src/fix_neb.rst | 29 ++++++++++++++++++----------- 1 file changed, 18 insertions(+), 11 deletions(-) diff --git a/doc/src/fix_neb.rst b/doc/src/fix_neb.rst index 4d0d301540..1d4c28f12b 100644 --- a/doc/src/fix_neb.rst +++ b/doc/src/fix_neb.rst @@ -120,14 +120,17 @@ and :math:`RD_{ideal}` is the ideal *RD* for which all the images are equally spaced. I.e. :math:`RD_{ideal} = (i-1) \cdot meanDist` when the climbing replica is off, where *i* is the replica number). The *meanDist* is the average distance between replicas. Note that in this -case the specified *Kspring* is in force units. - -Note that the *ideal* form of nudging can often be more effective at -keeping the replicas equally spaced. +case the specified *Kspring* is in force units. When the climbing replica +is on, :math:`RD_{ideal}` and :math:`meanDist` are calculated separately +each side of the climbing image. Note that the *ideal* form of nudging +can often be more effective at keeping the replicas equally spaced before +climbing, then equally spaced either side of the climbing image whilst +climbing. With a value of *equal* the spring force is computed as for *ideal* -before the climbing stage, then is computed to promote equidistant -spacing in energy rather than distance: +before the climbing stage, then during the climbing stage is modified +to promote equidistant absolute differnces in energy, rather than +distance, each side of the climbing image: .. math:: @@ -137,12 +140,16 @@ where *ED* is the cumulative sum of absolute energy differences: .. math:: - ED = \sum_{i`, by providing -optimal quadrature points. +*meanEdist* is the average absolute energy difference between +replicas up to the climbing image or from the climbing image +to the final image, for images before or after the climbing +image respectively. :math:`ED_{ideal}` is the corresponding +cumulative sum of average absolute energy differences in +each case, just as for *ideal*. This form of nudging is +to aid schemes which integrate forces along, or near to, +NEB pathways such as :doc:`fix_pafi `. ---------- From a8d0f94a5ae9b1221070dc88e548438f5813c5b0 Mon Sep 17 00:00:00 2001 From: Thomas Swinburne Date: Mon, 16 Jan 2023 14:21:31 +0100 Subject: [PATCH 24/25] small clean up --- doc/src/fix_neb.rst | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/doc/src/fix_neb.rst b/doc/src/fix_neb.rst index 1d4c28f12b..ec5d89a1b0 100644 --- a/doc/src/fix_neb.rst +++ b/doc/src/fix_neb.rst @@ -128,9 +128,10 @@ climbing, then equally spaced either side of the climbing image whilst climbing. With a value of *equal* the spring force is computed as for *ideal* -before the climbing stage, then during the climbing stage is modified -to promote equidistant absolute differnces in energy, rather than -distance, each side of the climbing image: +when the climbing replica is off, promoting equidistance. When the climbing +replica is on, the spring force is computed to promote equidistant +absolute differences in energy, rather than distance, each side of +the climbing image: .. math:: @@ -147,8 +148,8 @@ replicas up to the climbing image or from the climbing image to the final image, for images before or after the climbing image respectively. :math:`ED_{ideal}` is the corresponding cumulative sum of average absolute energy differences in -each case, just as for *ideal*. This form of nudging is -to aid schemes which integrate forces along, or near to, +each case, in close analogy to *ideal*. This form of nudging +is to aid schemes which integrate forces along, or near to, NEB pathways such as :doc:`fix_pafi `. ---------- From 503c51c070d5b3d10202c7dd2071e33f0367b396 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 16 Jan 2023 08:37:52 -0500 Subject: [PATCH 25/25] whitespace --- doc/src/fix_neb.rst | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/doc/src/fix_neb.rst b/doc/src/fix_neb.rst index ec5d89a1b0..5f22acfb2b 100644 --- a/doc/src/fix_neb.rst +++ b/doc/src/fix_neb.rst @@ -120,17 +120,17 @@ and :math:`RD_{ideal}` is the ideal *RD* for which all the images are equally spaced. I.e. :math:`RD_{ideal} = (i-1) \cdot meanDist` when the climbing replica is off, where *i* is the replica number). The *meanDist* is the average distance between replicas. Note that in this -case the specified *Kspring* is in force units. When the climbing replica +case the specified *Kspring* is in force units. When the climbing replica is on, :math:`RD_{ideal}` and :math:`meanDist` are calculated separately -each side of the climbing image. Note that the *ideal* form of nudging +each side of the climbing image. Note that the *ideal* form of nudging can often be more effective at keeping the replicas equally spaced before -climbing, then equally spaced either side of the climbing image whilst -climbing. +climbing, then equally spaced either side of the climbing image whilst +climbing. With a value of *equal* the spring force is computed as for *ideal* -when the climbing replica is off, promoting equidistance. When the climbing -replica is on, the spring force is computed to promote equidistant -absolute differences in energy, rather than distance, each side of +when the climbing replica is off, promoting equidistance. When the climbing +replica is on, the spring force is computed to promote equidistant +absolute differences in energy, rather than distance, each side of the climbing image: .. math:: @@ -141,15 +141,15 @@ where *ED* is the cumulative sum of absolute energy differences: .. math:: - ED = \sum_{i`. ----------