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

Merge pull request #537 from lammps/neb 

minor changes to NEB doc pages and examples
This commit is contained in:
sjplimp 2017-06-20 09:38:32 -06:00 committed by GitHub
parent 326a8a1289 50c9167913
commit 84b530cca1
17 changed files with 59 additions and 64 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

52
doc/src/fix_neb.txt
View File

@ -14,10 +14,10 @@ fix ID group-ID neb Kspring keyword value :pre
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, see nudge keyword) :l
Kspring = spring constant for parallel nudging force (force/distance units or force units, see parallel keyword) :l
zero or more keyword/value pairs may be appended :l
keyword = {nudge} or {perp} or {ends} :l
{nudge} value = {neigh} or {ideal}
keyword = {parallel} or {perp} or {end} :l
{parallel} value = {neigh} or {ideal}
{neigh} = parallel nudging force based on distance to neighbor replicas (Kspring = force/distance units)
{ideal} = parallel nudging force based on interpolated ideal position (Kspring = force units)
{perp} value = {Kspring2}
@ -28,7 +28,7 @@ keyword = {nudge} or {perp} or {ends} :l
{last} = apply force to last replica
{last/efirst} = apply force to last replica and set its target energy to that of first replica
{last/efirst/middle} = same as {last/efirst} plus prevent middle replicas having lower energy than first replica
{Kspring3} = spring constant for target energy term (1/distance units) :pre
{Kspring3} = spring constant for target energy term (1/distance units) :pre,ule
[Examples:]
@ -59,22 +59,22 @@ interatomic force Fi = -Grad(V) for each replica I is altered. For
all intermediate replicas (i.e. for 1 < I < N, except the climbing
replica) the force vector becomes:
Fi = -Grad(V) + (Grad(V) dot T') T' + Fnudge_parallel + Fspring_perp :pre
Fi = -Grad(V) + (Grad(V) dot T') T' + Fnudge_parallel + Fnudge_perp :pre
T' is the unit "tangent" vector for replica I and is a function of Ri,
Ri-1, Ri+1, and the potential energy of the 3 replicas; it points
roughly in the direction of (Ri+i - Ri-1); see the
"(Henkelman1)"_#Henkelman1 paper for details. Ri gives the atomic
"(Henkelman1)"_#Henkelman1 paper for details. Ri are the atomic
coordinates of replica I; Ri-1 and Ri+1 are the coordinates of its
neighbor replicas. The term (Grad(V) dot T') is used to remove the
component of the gradient parallel to the path which would tend to
distribute the replica unevenly along the path. Fnudge_parallel is an
artificial nudging force which is applied only in the tangent
direction and which maintains the equal spacing between replicas (see
below for more information). Fspring_perp is an optional artificial
spring which is applied only perpendicular to the tangent and which
prevent the paths from forming acute kinks (see below for more
information).
below for more information). Fnudge_perp is an optional artificial
spring which is applied in a direction perpendicular to the tangent
direction and which prevent the paths from forming acute kinks (see
below for more information).
In the second stage of the NEB calculation, the interatomic force Fi
for the climbing replica (the replica of highest energy after the
@ -86,7 +86,7 @@ and the relaxation procedure is continued to a new converged MEP.
:line
The keyword {nudge} specifies how the parallel nudging force is
The keyword {parallel} specifies how the parallel nudging force is
computed. With a value of {neigh}, the parallel nudging force is
computed as in "(Henkelman1)"_#Henkelman1 by connecting each
intermediate replica with the previous and the next image:
@ -98,7 +98,7 @@ units.
With a value of {ideal}, the spring force is computed as suggested in
"(WeinenE)"_#WeinenE :
Fnudge_parallel = -{Kspring} * (RD-RDideal) / (2 * meanDist) :pre
where RD is the "reaction coordinate" see "neb"_neb.html section, and
@ -113,15 +113,16 @@ keeping the replicas equally spaced.
:line
The keyword {perp} adds a spring force perpendicular to the path in
order to prevent the path from becoming too kinky. It
can significantly improve the convergence of the NEB calculation when
the resolution is poor. I.e. when too few replicas are used; see
The keyword {perp} specifies if and how a perpendicual nudging force
is computed. It adds a spring force perpendicular to the path in
order to prevent the path from becoming too kinky. It can
significantly improve the convergence of the NEB calculation when the
resolution is poor. I.e. when few replicas are used; see
"(Maras)"_#Maras1 for details.
The perpendicular spring force is given by
Fspring_perp = {Kspring2} * F(Ri-1,Ri,Ri+1) (Ri+1 + Ri-1 - 2 Ri) :pre
Fnudge_perp = {Kspring2} * F(Ri-1,Ri,Ri+1) (Ri+1 + Ri-1 - 2 Ri) :pre
where {Kspring2} is the specified value. F(Ri-1 Ri R+1) is a smooth
scalar function of the angle Ri-1 Ri Ri+1. It is equal to 0.0 when
@ -133,11 +134,11 @@ force is added.
:line
By default, no nudging 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 or last replica, to enable them to relax
toward a MEP while constraining their energy.
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.
The interatomic force Fi for the specified replica becomes:
@ -177,9 +178,10 @@ only be done if a particular intermediate replica has a lower energy
than the first replica. This should effectively prevent the
intermediate replicas from over-relaxing.
After converging a NEB calculation using an {estyle} of {last/efirst/middle},
you should check that all intermediate replicas have a larger energy than the
first replica. If this is not the case, the path is probably not a MEP.
After converging a NEB calculation using an {estyle} of
{last/efirst/middle}, you should check that all intermediate replicas
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

4
examples/neb/README
View File

@ -8,9 +8,7 @@ mpirun -np 3 lmp_g++ -partition 3x1 -in in.neb.sivac
mpirun -np 8 lmp_g++ -partition 4x2 -in in.neb.hop1
mpirun -np 8 lmp_g++ -partition 4x2 -in in.neb.hop2
mpirun -np 8 lmp_g++ -partition 4x2 -in in.neb.hop1.end
mpirun -np 6 lmp_g++ -partition 3x2 -in in.neb.sivac
mpirun -np 9 lmp_g++ -partition 3x3 -in in.neb.sivac
mpirun -np 8 lmp_g++ -partition 4x2 -in in.neb.sivac
Note that more than 4 replicas should be used for a precise estimate
of the activation energy corresponding to a transition.

2
examples/neb/in.neb.hop1
View File

@ -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 #nudge ideal
fix 2 nebatoms neb 1.0 parallel ideal
fix 3 all enforce2d
thermo 100

2
examples/neb/in.neb.hop1.end
View File

@ -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 nudge ideal end first 1.0
fix 2 nebatoms neb 1.0 parallel ideal end first 1.0
fix 3 all enforce2d
thermo 100

0
examples/neb/log.19June17.neb.hop1.end.g++.4 → examples/neb/log.19Jun17.neb.hop1.end.g++.4
View File

0
examples/neb/log.19June17.neb.hop1.end.g++.8 → examples/neb/log.19Jun17.neb.hop1.end.g++.8
View File

9
examples/neb/log.19Jun17.neb.hop1.g++.4 Normal file
View File

@ -0,0 +1,9 @@
LAMMPS (19 May 2017)
Running on 4 partitions of processors
Step MaxReplicaForce MaxAtomForce GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 ... RDN PEN
0 4327.2753 2746.3378 0.082169072 4.9967651 4514.5424 0.42933428 0.42323635 1.8941131 0 -3.0535948 0.33333333 -2.6242605 0.66666667 -2.7623811 1 -3.0474969
87 0.095951502 0.052720903 0.005588927 0.065110105 0.12467831 0.0071014928 0.0022798007 2.3003372 0 -3.0535967 0.32435271 -3.0473127 0.62805027 -3.0464952 1 -3.048775
Climbing replica = 3
Step MaxReplicaForce MaxAtomForce GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 ... RDN PEN
87 0.14137277 0.11108954 0.005588927 0.065110105 0.12467831 0.0071014928 0.0022798007 2.3003372 0 -3.0535967 0.32435271 -3.0473127 0.62805027 -3.0464952 1 -3.048775
124 0.099583263 0.085936899 0.0044220372 0.023873795 0.091308308 0.0071061754 0.0022863931 2.308121 0 -3.0535968 0.32223905 -3.0473329 0.61673898 -3.0464906 1 -3.048777

9
examples/neb/log.19Jun17.neb.hop1.g++.8 Normal file
View File

@ -0,0 +1,9 @@
LAMMPS (19 May 2017)
Running on 4 partitions of processors
Step MaxReplicaForce MaxAtomForce GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 ... RDN PEN
0 4327.2753 2746.3378 0.082169072 4.9967651 4514.5424 0.42933428 0.42323635 1.8941131 0 -3.0535948 0.33333333 -2.6242605 0.66666667 -2.7623811 1 -3.0474969
87 0.095951792 0.052720902 0.0055889267 0.065110091 0.12467831 0.0071014928 0.0022798007 2.3003372 0 -3.0535967 0.32435271 -3.0473127 0.62805027 -3.0464952 1 -3.048775
Climbing replica = 3
Step MaxReplicaForce MaxAtomForce GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 ... RDN PEN
87 0.14137297 0.11108954 0.0055889267 0.065110091 0.12467831 0.0071014928 0.0022798007 2.3003372 0 -3.0535967 0.32435271 -3.0473127 0.62805027 -3.0464952 1 -3.048775
124 0.099582186 0.08593683 0.0044220345 0.023873731 0.091308197 0.0071061754 0.0022863931 2.3081211 0 -3.0535968 0.32223904 -3.0473329 0.61673896 -3.0464906 1 -3.048777

0
examples/neb/log.19June17.neb.hop2.g++.4 → examples/neb/log.19Jun17.neb.hop2.g++.4
View File

0
examples/neb/log.19June17.neb.hop2.g++.8 → examples/neb/log.19Jun17.neb.hop2.g++.8
View File

0
examples/neb/log.19June17.neb.sivac.g++.4 → examples/neb/log.19Jun17.neb.sivac.g++.4
View File

0
examples/neb/log.19June17.neb.sivac.g++.8 → examples/neb/log.19Jun17.neb.sivac.g++.8
View File

10
examples/neb/log.19June17.neb.hop1.g++.4
View File

@ -1,10 +0,0 @@
LAMMPS (19 May 2017)
Running on 4 partitions of processors
Step MaxReplicaForce MaxAtomForce GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 ... RDN PEN
0 4327.2753 2746.3378 0.082169072 4.9967651 4514.5424 0.42933428 0.42323635 1.8941131 0 -3.0535948 0.33333333 -2.6242605 0.66666667 -2.7623811 1 -3.0474969
100 0.10482184 0.085218486 0.0051952047 0.04785954 0.19041553 0.0070900402 0.0022691875 2.3031875 0 -3.0535967 0.31839181 -3.0473647 0.63987598 -3.0465067 1 -3.0487759
111 0.096708467 0.07803707 0.0048656875 0.03613038 0.19671332 0.0070871172 0.0022668002 2.3052945 0 -3.0535968 0.31853431 -3.0473633 0.64178871 -3.0465096 1 -3.0487764
Climbing replica = 3
Step MaxReplicaForce MaxAtomForce GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 ... RDN PEN
111 0.2023467 0.1777038 0.0048656875 0.03613038 0.19671332 0.0070871172 0.0022668002 2.3052945 0 -3.0535968 0.31853431 -3.0473633 0.64178871 -3.0465096 1 -3.0487764
179 0.096874474 0.090676856 0.0034851031 0.0094134782 0.093630619 0.0071047642 0.0022856172 2.3122768 0 -3.0535969 0.31577311 -3.0473955 0.61798541 -3.0464922 1 -3.0487778

10
examples/neb/log.19June17.neb.hop1.g++.8
View File

@ -1,10 +0,0 @@
LAMMPS (19 May 2017)
Running on 4 partitions of processors
Step MaxReplicaForce MaxAtomForce GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 ... RDN PEN
0 4327.2753 2746.3378 0.082169072 4.9967651 4514.5424 0.42933428 0.42323635 1.8941131 0 -3.0535948 0.33333333 -2.6242605 0.66666667 -2.7623811 1 -3.0474969
100 0.10482171 0.085218406 0.0051952008 0.047859379 0.1904156 0.0070900401 0.0022691875 2.3031875 0 -3.0535967 0.31839181 -3.0473647 0.639876 -3.0465067 1 -3.0487759
111 0.096708718 0.078036984 0.0048656841 0.036130268 0.1967134 0.0070871172 0.0022668002 2.3052946 0 -3.0535968 0.31853431 -3.0473633 0.64178873 -3.0465096 1 -3.0487764
Climbing replica = 3
Step MaxReplicaForce MaxAtomForce GradV0 GradV1 GradVc EBF EBR RDT RD1 PE1 RD2 PE2 ... RDN PEN
111 0.20234693 0.17770387 0.0048656841 0.036130268 0.1967134 0.0070871172 0.0022668002 2.3052946 0 -3.0535968 0.31853431 -3.0473633 0.64178873 -3.0465096 1 -3.0487764
178 0.09975409 0.093814031 0.0035463662 0.010006594 0.096949208 0.0071042931 0.0022851195 2.312004 0 -3.0535969 0.31607934 -3.0473923 0.618931 -3.0464926 1 -3.0487777

4
src/REPLICA/fix_neb.cpp
View File

@ -49,8 +49,6 @@ FixNEB::FixNEB(LAMMPS *lmp, int narg, char **arg) :
displacements(NULL)
{
if (narg < 4) error->all(FLERR,"Illegal fix neb command");
kspring = force->numeric(FLERR,arg[3]);
@ -68,7 +66,7 @@ FixNEB::FixNEB(LAMMPS *lmp, int narg, char **arg) :
int iarg = 4;
while (iarg < narg) {
if (strcmp(arg[iarg],"nudge") == 0) {
if (strcmp(arg[iarg],"parallel") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix neb command");
if (strcmp(arg[iarg+1],"ideal") == 0) {
NEBLongRange = true;

8
src/dump.cpp
View File

@ -30,7 +30,7 @@
using namespace LAMMPS_NS;
#if defined(LMP_USE_LIBC_QSORT)
#if defined(LMP_QSORT)
// allocate space for static class variable
Dump *Dump::dumpptr;
#else
@ -692,7 +692,7 @@ void Dump::sort()
index[idsort[i]-idlo] = i;
}
#if defined(LMP_USE_LIBC_QSORT)
#if defined(LMP_QSORT)
if (!reorderflag) {
dumpptr = this;
for (i = 0; i < nme; i++) index[i] = i;
@ -728,7 +728,8 @@ void Dump::sort()
memcpy(&buf[i*size_one],&bufsort[index[i]*size_one],nbytes);
}
#if defined(LMP_USE_LIBC_QSORT)
#if defined(LMP_QSORT)
/* ----------------------------------------------------------------------
compare two atom IDs
called via qsort() in sort() method
@ -795,6 +796,7 @@ int Dump::bufcompare_reverse(const void *pi, const void *pj)
compare two atom IDs
called via merge_sort() in sort() method
------------------------------------------------------------------------- */
int Dump::idcompare(const int i, const int j, void *ptr)
{
tagint *idsort = ((Dump *)ptr)->idsort;

13
src/irregular.cpp
View File

@ -23,17 +23,13 @@
using namespace LAMMPS_NS;
#if defined(LMP_USE_LIBC_QSORT)
#if defined(LMP_QSORT)
// allocate space for static class variable
// prototype for non-class function
int *Irregular::proc_recv_copy;
static int compare_standalone(const void *, const void *);
#else
#include "mergesort.h"
// prototype for non-class function
static int compare_standalone(const int, const int, void *);
#endif
@ -433,7 +429,7 @@ int Irregular::create_atom(int n, int *sizes, int *proclist, int sortflag)
for (i = 0; i < nrecv_proc; i++) order[i] = i;
#if defined(LMP_USE_LIBC_QSORT)
#if defined(LMP_QSORT)
proc_recv_copy = proc_recv;
qsort(order,nrecv_proc,sizeof(int),compare_standalone);
#else
@ -464,7 +460,7 @@ int Irregular::create_atom(int n, int *sizes, int *proclist, int sortflag)
return nrecvsize;
}
#if defined(LMP_USE_LIBC_QSORT)
#if defined(LMP_QSORT)
/* ----------------------------------------------------------------------
comparison function invoked by qsort()
@ -705,12 +701,13 @@ int Irregular::create_data(int n, int *proclist, int sortflag)
for (i = 0; i < nrecv_proc; i++) order[i] = i;
#if defined(LMP_USE_LIBC_QSORT)
#if defined(LMP_QSORT)
proc_recv_copy = proc_recv;
qsort(order,nrecv_proc,sizeof(int),compare_standalone);
#else
merge_sort(order,nrecv_proc,(void *)proc_recv,compare_standalone);
#endif
int j;
for (i = 0; i < nrecv_proc; i++) {
j = order[i];