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