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; }