jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

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.
This commit is contained in:
Axel Kohlmeyer 2017-12-13 15:19:46 -05:00
parent 2f857c6eda
commit 9a71efc5d5
2 changed files with 39 additions and 37 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

17
doc/src/fix_neb.txt
View File

@ -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:]

59
src/REPLICA/fix_neb.cpp
View File

@ -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<EFinalIni) {
if (dot<0) prefactor = -dot - kspringFinal*(veng-EFinalIni);
else prefactor = -dot + kspringFinal*(veng-EFinalIni);
}
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
f[i][0] += prefactor *tangent[i][0];
@ -531,10 +534,10 @@ void FixNEB::min_post_force(int vflag)
double dotall;
MPI_Allreduce(&dot,&dotall,1,MPI_DOUBLE,MPI_SUM,world);
dot=dotall/tlen;
if (dot<0) prefactor = -dot - kspringFinal*(veng-vIni);
else prefactor = -dot + kspringFinal*(veng-vIni);
if (veng<vIni) {
if (dot<0) prefactor = -dot - kspringFinal*(veng-vIni);
else prefactor = -dot + kspringFinal*(veng-vIni);
}
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
f[i][0] += prefactor *tangent[i][0];
@ -589,7 +592,7 @@ void FixNEB::min_post_force(int vflag)
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
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];
dotSpringTangent += springF[i][0]*tangent[i][0] +
springF[i][1]*tangent[i][1] + springF[i][2]*tangent[i][2];}
@ -625,11 +628,11 @@ void FixNEB::min_post_force(int vflag)
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
f[i][0] += prefactor*tangent[i][0] +
f[i][0] += prefactor*tangent[i][0] +
AngularContr*(springF[i][0] - dotSpringTangent*tangent[i][0]);
f[i][1] += prefactor*tangent[i][1] +
f[i][1] += prefactor*tangent[i][1] +
AngularContr*(springF[i][1] - dotSpringTangent*tangent[i][1]);
f[i][2] += prefactor*tangent[i][2] +
f[i][2] += prefactor*tangent[i][2] +
AngularContr*(springF[i][2] - dotSpringTangent*tangent[i][2]);
}
}