diff --git a/src/REPLICA/fix_neb.cpp b/src/REPLICA/fix_neb.cpp index d10d29b3ce..652ba83755 100644 --- a/src/REPLICA/fix_neb.cpp +++ b/src/REPLICA/fix_neb.cpp @@ -168,6 +168,17 @@ void FixNEB::min_post_force(int vflag) pe->addstep(update->ntimestep+1); + // Compute norm of GradV for log output + + double **f = atom->f; + double fsq = 0.0; + for (int i = 0; i < nlocal; i++) { + fsq += f[i][0]*f[i][0]+f[i][1]*f[i][1]+f[i][2]*f[i][2]; + } + + MPI_Allreduce(&fsq,&gradvnorm,1,MPI_DOUBLE,MPI_MAX,world); + gradvnorm = sqrt(gradvnorm); + // if this is first or last replica, no change to forces, just return if (ireplica == 0 || ireplica == nreplica-1) { @@ -274,13 +285,12 @@ void FixNEB::min_post_force(int vflag) // see Henkelman & Jonsson 2000 paper, eqs 3,4,12 // see Henkelman & Jonsson 2000a paper, eq 5 - double **f = atom->f; - double dot = 0.0; - for (int i = 0; i < nlocal; i++) + for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) dot += f[i][0]*tangent[i][0] + f[i][1]*tangent[i][1] + f[i][2]*tangent[i][2]; + } double prefactor; if (ireplica == rclimber) prefactor = -2.0*dot; diff --git a/src/REPLICA/fix_neb.h b/src/REPLICA/fix_neb.h index 7e7d97c635..066259b166 100644 --- a/src/REPLICA/fix_neb.h +++ b/src/REPLICA/fix_neb.h @@ -28,6 +28,7 @@ class FixNEB : public Fix { public: double veng,plen,nlen; int rclimber; + double gradvnorm; FixNEB(class LAMMPS *, int, char **); ~FixNEB(); diff --git a/src/REPLICA/neb.cpp b/src/REPLICA/neb.cpp index dc3453dd49..bed7a179ac 100644 --- a/src/REPLICA/neb.cpp +++ b/src/REPLICA/neb.cpp @@ -138,9 +138,11 @@ void NEB::command(int narg, char **arg) if (me_universe == 0) { if (universe->uscreen) fprintf(universe->uscreen,"Step MaxReplicaForce MaxAtomForce " + "GradV0 GradV1 GradVc " "RD1 PE1 RD2 PE2 ... RDN PEN\n"); if (universe->ulogfile) fprintf(universe->ulogfile,"Step MaxReplicaForce MaxAtomForce " + "GradV0 GradV1 GradVc " "RD1 PE1 RD2 PE2 ... RDN PEN\n"); } print_status(); @@ -203,9 +205,11 @@ void NEB::command(int narg, char **arg) if (me_universe == 0) { if (universe->uscreen) fprintf(universe->uscreen,"Step MaxReplicaForce MaxAtomForce " + "GradV0 GradV1 GradVc " "RD1 PE1 RD2 PE2 ... RDN PEN\n"); if (universe->ulogfile) fprintf(universe->ulogfile,"Step MaxReplicaForce MaxAtomForce " + "GradV0 GradV1 GradVc " "RD1 PE1 RD2 PE2 ... RDN PEN\n"); } print_status(); @@ -381,11 +385,43 @@ void NEB::print_status() 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 better to be safe we + // take them straight from fix_neb + + double gradvnorm0, gradvnorm1, gradvnormc; + + int irep; + irep = 0; + if (me_universe == irep) gradvnorm0 = fneb->gradvnorm; + MPI_Bcast(&gradvnorm0,1,MPI_DOUBLE,irep,uworld); + irep = nreplica-1; + if (me_universe == irep) gradvnorm1 = fneb->gradvnorm; + MPI_Bcast(&gradvnorm1,1,MPI_DOUBLE,irep,uworld); + irep = fneb->rclimber; + if (irep > -1) { + if (me_universe == irep) gradvnormc = fneb->gradvnorm; + MPI_Bcast(&gradvnormc,1,MPI_DOUBLE,irep,uworld); + } else { + double vmax = all[0][0]; + int top = 0; + for (int m = 1; m < nreplica; m++) + if (vmax < all[m][0]) { + vmax = all[m][0]; + top = m; + } + irep = top; + if (me_universe == irep) gradvnormc = fneb->gradvnorm; + MPI_Bcast(&gradvnormc,1,MPI_DOUBLE,irep,uworld); + } + if (me_universe == 0) { if (universe->uscreen) { fprintf(universe->uscreen,"%d %g %g ",update->ntimestep, fmaxreplica,fmaxatom); + fprintf(universe->uscreen,"%g %g %g ", + gradvnorm0,gradvnorm1,gradvnormc); for (int i = 0; i < nreplica; i++) fprintf(universe->uscreen,"%g %g ",rdist[i],all[i][0]); fprintf(universe->uscreen,"\n"); @@ -393,6 +429,8 @@ void NEB::print_status() if (universe->ulogfile) { fprintf(universe->ulogfile,"%d %g %g ",update->ntimestep, fmaxreplica,fmaxatom); + fprintf(universe->ulogfile,"%g %g %g ", + gradvnorm0,gradvnorm1,gradvnormc); for (int i = 0; i < nreplica; i++) fprintf(universe->ulogfile,"%g %g ",rdist[i],all[i][0]); fprintf(universe->ulogfile,"\n");