From e03cc994673e6b43ef41368f772d46f9ba3b7086 Mon Sep 17 00:00:00 2001 From: Emile Maras Date: Fri, 2 Jun 2017 23:42:16 +0200 Subject: [PATCH] made the command options more lammps standard style --- doc/src/fix_neb.txt | 54 +++++++++++++++++++++------------ examples/neb/in.neb.hop1 | 2 +- examples/neb/in.neb.hop1freeend | 2 +- examples/neb/in.neb.hop2 | 2 +- examples/neb/in.neb.sivac | 2 +- src/REPLICA/fix_neb.cpp | 49 +++++++++++++++++------------- src/REPLICA/fix_neb.h | 2 +- 7 files changed, 67 insertions(+), 46 deletions(-) diff --git a/doc/src/fix_neb.txt b/doc/src/fix_neb.txt index a5c4bf4396..cd0fa23e85 100644 --- a/doc/src/fix_neb.txt +++ b/doc/src/fix_neb.txt @@ -12,22 +12,33 @@ fix neb command :h3 fix ID group-ID neb Kspring keyword value :pre -ID, group-ID are documented in "fix"_fix.html command -neb = style name of this fix command -Kspring = parallel spring constant (force/distance units) -keyword = {idealpos} or {neigh} or {perp} or {freeend} :ul - {idealpos} = each replica is attached with a spring to its interpolated ideal position (default) - {neigh} = each replica is connected with spring to the previous and next replica. - {perp} value = set spring constant for the perpendicular spring to {value} - {freeend} flag = set behavior for the end points - flag = {ini} or {final} or {finaleini} or {final2eini} +ID, group-ID are documented in "fix"_fix.html command :ulb,l +neb = style name of this fix command :l +Kspring = parallel spring constant (force/distance units or force units) :l +zero or more keyword/value pairs may be appended :l +keyword = {nudg_style} or {perp} or {freend} or {freend_k_spring} :l + {nudg_style} value = {neigh} or {idealpos} + {neigh} = the parallel nudging force is calculated from the distance to neighbouring replicas (in this case, Kspring is in force/distance units) + {idealpos} = the parallel nudging force is proportional to the distance between the replica and its interpolated ideal position (in this case Kspring is in force units) + {perp} value {none} or kspring2 + {none} = no perpendicular spring force is applied + {freeend} value = {none} or {ini} or {final} or {finaleini} or {final2eini} + {none} = no nudging force is apply to the first and last replicas + {ini} = set the first replica to be a free end + {final} = set the last replica to be a free end + {finaleini} = set the last replica to be a free end and set its target energy as that of the first replica + {final2eini} = same as {finaleini} plus prevent intermediate replicas to have a lower energy than the first replica + {freeend_kspring} value = kspring2 + kspring2 = spring constant of the perpendicular spring force (per distance units) +flag = set behavior for the end points + flag = :pre [Examples:] fix 1 active neb 10.0 fix 2 all neb 1.0 perp 1.0 freeend final -fix 1 all neb 1.0 neigh freeend final2eini :pre +fix 1 all neb 1.0 nudg_style idealpos freeend final2eini freend_kspring 1:pre [Description:] @@ -109,11 +120,11 @@ replica, a free end neb calculation with the value {finaleini} or :line -The keywords {idealpos} and {neigh} allow to specify how to parallel -spring force is computed. If the keyword {idealpos} is used or by -default, the spring force is computed as suggested in "(E)"_#E : +The keyword {nudg_style} allow to specify how to parallel +nudging force is computed. With a value of idealpos, the spring +force is computed as suggested in "(E)"_#E : -Fspringparallel=-{Kspring}* (RD-RDideal)/(2 meanDist) :pre +Fnudgparallel=-{Kspring}* (RD-RDideal)/(2 meanDist) :pre where RD is the "reaction coordinate" see "neb"_neb.html section, and RDideal is the ideal RD for which all the images are equally spaced @@ -121,13 +132,14 @@ RDideal is the ideal RD for which all the images are equally spaced is the replica number). The meanDist is the average distance between replicas. -If the keyword {neigh} is used, the parallel spring force is computed -as in "(Henkelman1)"_#Henkelman1 by connecting each intermediate -replica with the previous and the next image: +When {nudg_style} has a value of neigh (or by default), the parallel +nudging force is computed as in "(Henkelman1)"_#Henkelman1 by +connecting each intermediate replica with the previous and the next +image: -Fspringparallel= {Kspring}* (|Ri+1 - Ri| - |Ri - Ri-1|) :pre +Fnudgparallel= {Kspring}* (|Ri+1 - Ri| - |Ri - Ri-1|) :pre -The parallel spring force associated with the key word idealpos should +The parallel nudging force associated with the key word idealpos should usually be more efficient at keeping the images equally spaced. :line @@ -172,7 +184,9 @@ for more info on packages. [Default:] -none +The option defaults are nudg_style = neigh, perp = none, freeend = none and freend_kspring = 1. + +:line :link(Henkelman1) [(Henkelman1)] Henkelman and Jonsson, J Chem Phys, 113, 9978-9985 (2000). diff --git a/examples/neb/in.neb.hop1 b/examples/neb/in.neb.hop1 index 495b2faa93..b874d1ba32 100644 --- a/examples/neb/in.neb.hop1 +++ b/examples/neb/in.neb.hop1 @@ -51,7 +51,7 @@ set group nebatoms type 3 group nonneb subtract all nebatoms fix 1 lower setforce 0.0 0.0 0.0 -fix 2 nebatoms neb 1.0 neigh +fix 2 nebatoms neb 1.0 nudg_style idealpos fix 3 all enforce2d thermo 100 diff --git a/examples/neb/in.neb.hop1freeend b/examples/neb/in.neb.hop1freeend index fc0677f30a..fa90e9a98c 100644 --- a/examples/neb/in.neb.hop1freeend +++ b/examples/neb/in.neb.hop1freeend @@ -41,7 +41,7 @@ set group nebatoms type 3 group nonneb subtract all nebatoms fix 1 lower setforce 0.0 0.0 0.0 -fix 2 nebatoms neb 1.0 freeend ini +fix 2 nebatoms neb 1.0 nudg_style idealpos freeend ini fix 3 all enforce2d thermo 100 diff --git a/examples/neb/in.neb.hop2 b/examples/neb/in.neb.hop2 index 9a8986f454..242de759fa 100644 --- a/examples/neb/in.neb.hop2 +++ b/examples/neb/in.neb.hop2 @@ -53,7 +53,7 @@ set group nebatoms type 3 group nonneb subtract all nebatoms fix 1 lower setforce 0.0 0.0 0.0 -fix 2 nebatoms neb 1.0 neigh +fix 2 nebatoms neb 1.0 fix 3 all enforce2d thermo 100 diff --git a/examples/neb/in.neb.sivac b/examples/neb/in.neb.sivac index 2fdf46f9d5..941d063b90 100644 --- a/examples/neb/in.neb.sivac +++ b/examples/neb/in.neb.sivac @@ -66,7 +66,7 @@ minimize 1.0e-6 1.0e-4 1000 10000 reset_timestep 0 -fix 1 all neb 1.0 neigh +fix 1 all neb 1.0 thermo 100 diff --git a/src/REPLICA/fix_neb.cpp b/src/REPLICA/fix_neb.cpp index f66167f2a2..b17315ca0d 100644 --- a/src/REPLICA/fix_neb.cpp +++ b/src/REPLICA/fix_neb.cpp @@ -46,12 +46,13 @@ FixNEB::FixNEB(LAMMPS *lmp, int narg, char **arg) : displacements(NULL) { - NEBLongRange=true; - StandardNEB=PerpSpring=FreeEndIni=FreeEndFinal=false; + NEBLongRange=false; + StandardNEB=true; + PerpSpring=FreeEndIni=FreeEndFinal=false; FreeEndFinalWithRespToEIni=FinalAndInterWithRespToEIni=false; kspringPerp=0.0; - + kspring2=1.0; if (narg < 4) error->all(FLERR,"Illegal fix neb command, argument missing"); @@ -62,21 +63,23 @@ FixNEB::FixNEB(LAMMPS *lmp, int narg, char **arg) : int iarg =4; while (iarg < narg) { - if (strcmp (arg[iarg],"idealpos")==0) { - NEBLongRange = true; - iarg+=1; - } else if (strcmp (arg[iarg],"neigh")==0) { - NEBLongRange = false; - StandardNEB = true; - iarg+=1; - } else if (strcmp (arg[iarg],"perp")==0) { + if (strcmp (arg[iarg],"nudg_style")==0) { + if (strcmp (arg[iarg+1],"idealpos")==0) { + NEBLongRange = true; + iarg+=2;} + else if (strcmp (arg[iarg+1],"neigh")==0) { + NEBLongRange = false; + StandardNEB = true; + iarg+=2;} + else error->all(FLERR,"Illegal fix neb command. Unknown keyword");} + else if (strcmp (arg[iarg],"perp")==0) { PerpSpring=true; kspringPerp = force->numeric(FLERR,arg[iarg+1]); if (kspringPerp < 0.0) error->all(FLERR,"Illegal fix neb command. " "The perpendicular spring force was not provided properly"); - iarg+=2; - } else if (strcmp (arg[iarg],"freeend")==0) { + iarg+=2;} + else if (strcmp (arg[iarg],"freeend")==0) { if (strcmp (arg[iarg+1],"ini")==0) FreeEndIni=true; else if (strcmp (arg[iarg+1],"final")==0) @@ -86,8 +89,12 @@ FixNEB::FixNEB(LAMMPS *lmp, int narg, char **arg) : else if (strcmp (arg[iarg+1],"final2eini")==0) { FinalAndInterWithRespToEIni=true; FreeEndFinalWithRespToEIni=true;} - iarg+=2; - } else error->all(FLERR,"Illegal fix neb command. Unknown keyword"); + else if (strcmp (arg[iarg+1],"none")!=0) error->all(FLERR,"Illegal fix neb command. Unknown keyword"); + iarg+=2;} + else if (strcmp (arg[iarg],"freeend_kspring")==0) { + kspring2=force->numeric(FLERR,arg[iarg+1]); + iarg+=2; } + else error->all(FLERR,"Illegal fix neb command. Unknown keyword"); } // nreplica = number of partitions @@ -468,8 +475,8 @@ 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 - (veng-EIniIni); - else prefactor = -dot + (veng-EIniIni); + if (dot<0) prefactor = -dot - kspring2*(veng-EIniIni); + else prefactor = -dot + kspring2*(veng-EIniIni); for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { @@ -486,8 +493,8 @@ 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 - (veng-EFinalIni); - else prefactor = -dot + (veng-EFinalIni); + if (dot<0) prefactor = -dot - kspring2*(veng-EFinalIni); + else prefactor = -dot + kspring2*(veng-EFinalIni); for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { @@ -504,8 +511,8 @@ 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 - (veng-vIni); - else prefactor = -dot + (veng-vIni); + if (dot<0) prefactor = -dot - kspring2*(veng-vIni); + else prefactor = -dot + kspring2*(veng-vIni); for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { diff --git a/src/REPLICA/fix_neb.h b/src/REPLICA/fix_neb.h index 1582912dac..7e9e6db865 100644 --- a/src/REPLICA/fix_neb.h +++ b/src/REPLICA/fix_neb.h @@ -38,7 +38,7 @@ class FixNEB : public Fix { private: int me,nprocs,nprocs_universe; - double kspring,kspringPerp,EIniIni,EFinalIni; + double kspring,kspring2,kspringPerp,EIniIni,EFinalIni; bool StandardNEB,NEBLongRange,PerpSpring,FreeEndIni,FreeEndFinal; bool FreeEndFinalWithRespToEIni,FinalAndInterWithRespToEIni; int ireplica,nreplica;