From 9a71efc5d58fd83fecd6132a6426d2fd2b6592b8 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 13 Dec 2017 15:19:46 -0500 Subject: [PATCH] fix neb bugfix from Emile Maras NEB was not working fine when using multiple proc per replica and the keywords last/efirst or last/efirst/middle I have corrected this in the enclosed fix_neb.cpp I also slightly modified the nudging for this free end so that it would be applied only when the target energy is larger than the energy. Anyway if the target energy is lower than the energy, the replica should relax toward the target energy without adding any nudging. I also modified the documentation according to this change. --- doc/src/fix_neb.txt | 17 ++++++------ src/REPLICA/fix_neb.cpp | 59 ++++++++++++++++++++++------------------- 2 files changed, 39 insertions(+), 37 deletions(-) diff --git a/doc/src/fix_neb.txt b/doc/src/fix_neb.txt index 73b3e31266..be4766eae9 100644 --- a/doc/src/fix_neb.txt +++ b/doc/src/fix_neb.txt @@ -138,17 +138,17 @@ By default, no additional forces act on the first and last replicas during the NEB relaxation, so these replicas simply relax toward their respective local minima. By using the key word {end}, additional forces can be applied to the first and/or last replicas, to enable -them to relax toward a MEP while constraining their energy. +them to relax toward a MEP while constraining their energy E to the +target energy ETarget. -The interatomic force Fi for the specified replica becomes: +If ETarget>E, the interatomic force Fi for the specified replica becomes: 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 :pre -where E is the current energy of the replica and ETarget is the target -energy. The "spring" constant on the difference in energies is the -specified {Kspring3} value. +The "spring" constant on the difference in energies is the specified +{Kspring3} value. When {estyle} is specified as {first}, the force is applied to the first replica. When {estyle} is specified as {last}, the force is @@ -183,10 +183,9 @@ After converging a NEB calculation using an {estyle} of have a larger energy than the first replica. If this is not the case, the path is probably not a MEP. -Finally, note that if the last replica converges toward a local -minimum which has a larger energy than the energy of the first -replica, a NEB calculation using an {estyle} of {last/efirst} or -{last/efirst/middle} cannot reach final convergence. +Finally, note that the last replica may never reach the target energy +if it is stuck in a local minima which has a larger energy than the +target energy. [Restart, fix_modify, output, run start/stop, minimize info:] diff --git a/src/REPLICA/fix_neb.cpp b/src/REPLICA/fix_neb.cpp index 13a6ad1d62..98e5539872 100644 --- a/src/REPLICA/fix_neb.cpp +++ b/src/REPLICA/fix_neb.cpp @@ -94,26 +94,26 @@ FixNEB::FixNEB(LAMMPS *lmp, int narg, char **arg) : if (iarg+3 > narg) error->all(FLERR,"Illegal fix neb command"); if (strcmp(arg[iarg+1],"first") == 0) { FreeEndIni = true; - kspringIni = force->numeric(FLERR,arg[iarg+2]); + kspringIni = force->numeric(FLERR,arg[iarg+2]); } else if (strcmp(arg[iarg+1],"last") == 0) { FreeEndFinal = true; FinalAndInterWithRespToEIni = false; FreeEndFinalWithRespToEIni = false; - kspringFinal = force->numeric(FLERR,arg[iarg+2]); + kspringFinal = force->numeric(FLERR,arg[iarg+2]); } else if (strcmp(arg[iarg+1],"last/efirst") == 0) { FreeEndFinal = false; FinalAndInterWithRespToEIni = false; FreeEndFinalWithRespToEIni = true; - kspringFinal = force->numeric(FLERR,arg[iarg+2]); + kspringFinal = force->numeric(FLERR,arg[iarg+2]); } else if (strcmp(arg[iarg+1],"last/efirst/middle") == 0) { FreeEndFinal = false; FinalAndInterWithRespToEIni = true; FreeEndFinalWithRespToEIni = true; - kspringFinal = force->numeric(FLERR,arg[iarg+2]); + kspringFinal = force->numeric(FLERR,arg[iarg+2]); } else error->all(FLERR,"Illegal fix neb command"); iarg += 3; - + } else error->all(FLERR,"Illegal fix neb command"); } @@ -298,12 +298,14 @@ void FixNEB::min_post_force(int vflag) if (ireplica == 0) vIni=veng; if (FreeEndFinalWithRespToEIni) { - if (me == 0) { + if (cmode == SINGLE_PROC_DIRECT || cmode == SINGLE_PROC_MAP) { int procFirst; procFirst=universe->root_proc[0]; MPI_Bcast(&vIni,1,MPI_DOUBLE,procFirst,uworld); - } - if (cmode == MULTI_PROC) { + }else { + if (me == 0) + MPI_Bcast(&vIni,1,MPI_DOUBLE,0,rootworld); + MPI_Bcast(&vIni,1,MPI_DOUBLE,0,world); } } @@ -313,8 +315,8 @@ void FixNEB::min_post_force(int vflag) // if (me == 0 ) if (update->ntimestep == 0) { EIniIni = veng; - // if (cmode == MULTI_PROC) - // MPI_Bcast(&EIniIni,1,MPI_DOUBLE,0,world); + // if (cmode == MULTI_PROC) + // MPI_Bcast(&EIniIni,1,MPI_DOUBLE,0,world); } }*/ @@ -360,7 +362,7 @@ void FixNEB::min_post_force(int vflag) 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] + + dot += f[i][0]*tangent[i][0] + f[i][1]*tangent[i][1] + f[i][2]*tangent[i][2]; } } @@ -383,9 +385,9 @@ void FixNEB::min_post_force(int vflag) tangent[i][0]=delxn; tangent[i][1]=delyn; tangent[i][2]=delzn; - tlen += tangent[i][0]*tangent[i][0] + + 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] + + dot += f[i][0]*tangent[i][0] + f[i][1]*tangent[i][1] + f[i][2]*tangent[i][2]; } } @@ -431,15 +433,15 @@ void FixNEB::min_post_force(int vflag) } nlen += delxn*delxn + delyn*delyn + delzn*delzn; - tlen += tangent[i][0]*tangent[i][0] + + 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] + + 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] + + 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] + + 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); @@ -514,9 +516,10 @@ void FixNEB::min_post_force(int vflag) MPI_Allreduce(&dot,&dotall,1,MPI_DOUBLE,MPI_SUM,world); dot=dotall/tlen; - if (dot<0) prefactor = -dot - kspringFinal*(veng-EFinalIni); - else prefactor = -dot + kspringFinal*(veng-EFinalIni); - + if (veng