diff --git a/doc/src/min_modify.rst b/doc/src/min_modify.rst index 092e655d7d..1f137abe6a 100644 --- a/doc/src/min_modify.rst +++ b/doc/src/min_modify.rst @@ -65,32 +65,30 @@ highly overlapped atoms from being moved long distances (e.g. through another atom) due to large forces. The choice of line search algorithm for the *cg* and *sd* minimization -styles can be selected via the *line* keyword. -The default *quadratic* line search algorithm starts out using -the robust backtracking method described below. However, once -the system gets close to a local -minimum and the linesearch steps get small, so that the energy -is approximately quadratic in the step length, it uses the -estimated location of zero gradient as the linesearch step, -provided the energy change is downhill. -This becomes more efficient than backtracking -for highly-converged relaxations. The *forcezero* -line search algorithm is similar to *quadratic*\ . -It may be more efficient than *quadratic* on some systems. +styles can be selected via the *line* keyword. The default +*quadratic* line search algorithm starts out using the robust +backtracking method described below. However, once the system gets +close to a local minimum and the linesearch steps get small, so that +the energy is approximately quadratic in the step length, it uses the +estimated location of zero gradient as the linesearch step, provided +the energy change is downhill. This becomes more efficient than +backtracking for highly-converged relaxations. The *forcezero* line +search algorithm is similar to *quadratic*\ . It may be more +efficient than *quadratic* on some systems. -The backtracking search is robust and should always find a local energy -minimum. However, it will "converge" when it can no longer reduce the -energy of the system. Individual atom forces may still be larger than -desired at this point, because the energy change is measured as the -difference of two large values (energy before and energy after) and -that difference may be smaller than machine epsilon even if atoms -could move in the gradient direction to reduce forces further. +The backtracking search is robust and should always find a local +energy minimum. However, it will "converge" when it can no longer +reduce the energy of the system. Individual atom forces may still be +larger than desired at this point, because the energy change is +measured as the difference of two large values (energy before and +energy after) and that difference may be smaller than machine epsilon +even if atoms could move in the gradient direction to reduce forces +further. -The choice of a norm can be modified for the min styles *cg*\ , *sd*\ , -*quickmin*\ , *fire*\ , *spin*\ , *spin/cg* and *spin/lbfgs* using -the *norm* keyword. -The default *two* norm computes the 2-norm (Euclidean length) of the -global force vector: +The choice of a norm can be modified for the min styles *cg*\ , *sd*\ +, *quickmin*\ , *fire*\ , *spin*\ , *spin/cg* and *spin/lbfgs* using +the *norm* keyword. The default *two* norm computes the 2-norm +(Euclidean length) of the global force vector: .. image:: Eqs/norm_two.jpg :align: center @@ -111,57 +109,57 @@ all atoms in the system: For the min styles *spin*\ , *spin/cg* and *spin/lbfgs*\ , the force norm is replaced by the spin-torque norm. -Keywords *alpha\_damp* and *discrete\_factor* only make sense when -a :doc:`min\_spin ` command is declared. -Keyword *alpha\_damp* defines an analog of a magnetic Gilbert -damping. It defines a relaxation rate toward an equilibrium for -a given magnetic system. -Keyword *discrete\_factor* defines a discretization factor for the -adaptive timestep used in the *spin* minimization. -See :doc:`min\_spin ` for more information about those +Keywords *alpha\_damp* and *discrete\_factor* should only be used when +a :doc:`min\_spin ` command is declared. Keyword +*alpha\_damp* defines an analog of a magnetic Gilbert damping. It +defines a relaxation rate toward an equilibrium for a given magnetic +system. Keyword *discrete\_factor* defines a discretization factor +for the adaptive timestep used in the *spin* minimization. See +:doc:`min\_spin ` for more information about those quantities. The choice of a line search algorithm for the *spin/cg* and -*spin/lbfgs* styles can be specified via the *line* keyword. -The *spin\_cubic* and *spin\_none* only make sense when one of those -two minimization styles is declared. -The *spin\_cubic* performs the line search based on a cubic interpolation -of the energy along the search direction. The *spin\_none* keyword -deactivates the line search procedure. -The *spin\_none* is a default value for *line* keyword for both *spin/lbfgs* -and *spin/cg*\ . Convergence of *spin/lbfgs* can be more robust if -*spin\_cubic* line search is used. +*spin/lbfgs* styles can be specified via the *line* keyword. The +*spin\_cubic* and *spin\_none* only make sense when one of those two +minimization styles is declared. The *spin\_cubic* performs the line +search based on a cubic interpolation of the energy along the search +direction. The *spin\_none* keyword deactivates the line search +procedure. The *spin\_none* is a default value for *line* keyword for +both *spin/lbfgs* and *spin/cg*\ . Convergence of *spin/lbfgs* can be +more robust if *spin\_cubic* line search is used. -The Newton *integrator* used for *fire2* minimization can be selected to be -either the symplectic Euler (\ *eulerimplicit*\ ) or velocity Verlet (\ *verlet*\ ). -*tmax* define the maximum value for the adaptive timestep -during a *fire2* minimization. It is multiplication factor applied -to the current :doc:`timestep ` (not in time unit). For example, -*tmax* = 4.0 in metal :doc:`units ` means that the maximum value -the timestep can reach during a minimization is 4fs (with the default -:doc:`timestep ` value). Note that parameters defaults has been -chosen to be reliable in most cases, but one should consider adjusting -:doc:`timestep ` and *tmax* to optimize the minimization for large -or complex systems. -Others parameters of the *fire2* minimization can be tuned (\ *tmin*\ , *delaystep*\ , -*dtgrow*\ , *dtshrink*\ , *alpha0*\ , and *alphashrink*\ ). Please refer to the article -describing the *fire2* :doc:`min\_style `. +The Newton *integrator* used for *fire2* minimization can be selected +to be either the symplectic Euler (\ *eulerimplicit*\ ) or velocity +Verlet (\ *verlet*\ ). *tmax* define the maximum value for the +adaptive timestep during a *fire2* minimization. It is multiplication +factor applied to the current :doc:`timestep ` (not in time +unit). For example, *tmax* = 4.0 in metal :doc:`units ` means +that the maximum value the timestep can reach during a minimization is +4fs (with the default :doc:`timestep ` value). Note that +parameter defaults has been chosen to be reliable in most cases, but +one should consider adjusting :doc:`timestep ` and *tmax* to +optimize the minimization for large or complex systems. Other +parameters of the *fire2* minimization can be tuned (\ *tmin*\ , +*delaystep*\ , *dtgrow*\ , *dtshrink*\ , *alpha0*\ , and +*alphashrink*\ ). Please refer to the article describing the *fire2* +:doc:`min\_style `. -An additional stopping criteria *vdfmax* is added in order to avoid unnecessary looping -when it is reasonable to think the system will not be relaxed further. -Note that in this case the system will NOT be relaxed. This could -happen when the system comes to be stuck in a local basin of the phase space. -*vdfmax* is the maximum number of consecutive iterations with P(t) < 0. -For debugging purposes, it is possible to switch off the inertia correction -(\ *halfstepback* = *no*\ ) and the initial delay (\ *initialdelay* = *no*\ ). +An additional stopping criteria *vdfmax* is added in order to avoid +unnecessary looping when it is reasonable to think the system will not +be relaxed further. Note that in this case the system will NOT be +relaxed. This could happen when the system comes to be stuck in a +local basin of the phase space. *vdfmax* is the maximum number of +consecutive iterations with P(t) < 0. For debugging purposes, it is +possible to switch off the inertia correction (\ *halfstepback* = +*no*\ ) and the initial delay (\ *initialdelay* = *no*\ ). Restrictions """""""""""" -For magnetic GNEB calculations, only *spin\_none* value for *line* keyword can be used -when styles *spin/cg* and *spin/lbfgs* are employed. -See :doc:`neb/spin ` for more explanation. +For magnetic GNEB calculations, only *spin\_none* value for *line* +keyword can be used when styles *spin/cg* and *spin/lbfgs* are +employed. See :doc:`neb/spin ` for more explanation. Related commands """""""""""""""" @@ -173,11 +171,11 @@ Default The option defaults are dmax = 0.1, line = quadratic and norm = two. -For the *spin*\ , *spin/cg* and *spin/lbfgs* styles, the -option defaults are alpha\_damp = 1.0, discrete\_factor = 10.0, -line = spin\_none, and norm = euclidean. +For the *spin*\ , *spin/cg* and *spin/lbfgs* styles, the option +defaults are alpha\_damp = 1.0, discrete\_factor = 10.0, line = +spin\_none, and norm = euclidean. -For the *fire2* style, the option defaults are -integrator = eulerimplicit, tmax = 10.0, tmin = 0.02, -delaystep = 20, dtgrow = 1.1, dtshrink = 0.5, alpha0 = 0.25, alphashrink = 0.99, -vdfmax = 2000, halfstepback = yes and initialdelay = yes. +For the *fire2* style, the option defaults are integrator = +eulerimplicit, tmax = 10.0, tmin = 0.02, delaystep = 20, dtgrow = 1.1, +dtshrink = 0.5, alpha0 = 0.25, alphashrink = 0.99, vdfmax = 2000, +halfstepback = yes and initialdelay = yes. diff --git a/doc/src/min_style.rst b/doc/src/min_style.rst index 9921369b56..af8aa95425 100644 --- a/doc/src/min_style.rst +++ b/doc/src/min_style.rst @@ -26,8 +26,8 @@ Examples Description """"""""""" -Choose a minimization algorithm to use when a :doc:`minimize ` -command is performed. +Choose a minimization algorithm to use when a :doc:`minimize +` command is performed. Style *cg* is the Polak-Ribiere version of the conjugate gradient (CG) algorithm. At each iteration the force gradient is combined with the @@ -55,55 +55,55 @@ descent will not converge as quickly as CG, but may be more robust in some situations. Style *quickmin* is a damped dynamics method described in -:ref:`(Sheppard) `, where the damping parameter is related to the -projection of the velocity vector along the current force vector for -each atom. The velocity of each atom is initialized to 0.0 by this -style, at the beginning of a minimization. +:ref:`(Sheppard) `, where the damping parameter is related +to the projection of the velocity vector along the current force +vector for each atom. The velocity of each atom is initialized to 0.0 +by this style, at the beginning of a minimization. -Style *fire* is a damped dynamics method described in -:ref:`(Bitzek) `, which is similar to *quickmin* but adds a variable -timestep and alters the projection operation to maintain components of -the velocity non-parallel to the current force vector. The velocity -of each atom is initialized to 0.0 by this style, at the beginning of -a minimization. +Style *fire* is a damped dynamics method described in :ref:`(Bitzek) +`, which is similar to *quickmin* but adds a variable timestep +and alters the projection operation to maintain components of the +velocity non-parallel to the current force vector. The velocity of +each atom is initialized to 0.0 by this style, at the beginning of a +minimization. -Style *fire2* is an optimization of the style *fire*\ , including different -time integration schemes, described in :ref:`(Guenole) `. +Style *fire2* is an optimization of the style *fire*\ , including +different time integration schemes, described in :ref:`(Guenole) +`. -Style *spin* is a damped spin dynamics with an adaptive -timestep. +Style *spin* is a damped spin dynamics with an adaptive timestep. -Style *spin/cg* uses an orthogonal spin optimization (OSO) -combined to a conjugate gradient (CG) approach to minimize spin -configurations. +Style *spin/cg* uses an orthogonal spin optimization (OSO) combined to +a conjugate gradient (CG) approach to minimize spin configurations. -Style *spin/lbfgs* uses an orthogonal spin optimization (OSO) -combined to a limited-memory Broyden-Fletcher-Goldfarb-Shanno -(LBFGS) approach to minimize spin configurations. +Style *spin/lbfgs* uses an orthogonal spin optimization (OSO) combined +to a limited-memory Broyden-Fletcher-Goldfarb-Shanno (LBFGS) approach +to minimize spin configurations. -See the :doc:`min/spin ` doc page for more information -about the *spin*\ , *spin/cg* and *spin/lbfgs* styles. +See the :doc:`min/spin ` doc page for more information about +the *spin*\ , *spin/cg* and *spin/lbfgs* styles. -Either the *quickmin*\ , *fire* and *fire2* styles are useful in the context of -nudged elastic band (NEB) calculations via the :doc:`neb ` command. +Either the *quickmin*\ , *fire* and *fire2* styles are useful in the +context of nudged elastic band (NEB) calculations via the :doc:`neb +` command. -Either the *spin*\ , *spin/cg* and *spin/lbfgs* styles are useful -in the context of magnetic geodesic nudged elastic band (GNEB) calculations -via the :doc:`neb/spin ` command. +Either the *spin*\ , *spin/cg* and *spin/lbfgs* styles are useful in +the context of magnetic geodesic nudged elastic band (GNEB) +calculations via the :doc:`neb/spin ` command. .. note:: The damped dynamic minimizers use whatever timestep you have - defined via the :doc:`timestep ` command. Often they will - converge more quickly if you use a timestep about 10x larger than you - would normally use for dynamics simulations. + defined via the :doc:`timestep ` command. Often they + will converge more quickly if you use a timestep about 10x larger + than you would normally use for dynamics simulations. .. note:: The *quickmin*\ , *fire*\ , *hftn*\ , and *cg/kk* styles do not yet - support the use of the :doc:`fix box/relax ` command or - minimizations involving the electron radius in :doc:`eFF ` - models. + support the use of the :doc:`fix box/relax ` command + or minimizations involving the electron radius in :doc:`eFF + ` models. ---------- @@ -112,17 +112,19 @@ via the :doc:`neb/spin ` command. Styles with a *gpu*\ , *intel*\ , *kk*\ , *omp*\ , or *opt* suffix are functionally the same as the corresponding style without the suffix. They have been optimized to run faster, depending on your available -hardware, as discussed on the :doc:`Speed packages ` doc -page. The accelerated styles take the same arguments and should +hardware, as discussed on the :doc:`Speed packages ` +doc page. The accelerated styles take the same arguments and should produce the same results, except for round-off and precision issues. These accelerated styles are part of the GPU, USER-INTEL, KOKKOS, USER-OMP and OPT packages, respectively. They are only enabled if -LAMMPS was built with those packages. See the :doc:`Build package ` doc page for more info. +LAMMPS was built with those packages. See the :doc:`Build package +` doc page for more info. You can specify the accelerated styles explicitly in your input script -by including their suffix, or you can use the :doc:`-suffix command-line switch ` when you invoke LAMMPS, or you can use the -:doc:`suffix ` command in your input script. +by including their suffix, or you can use the :doc:`-suffix +command-line switch ` when you invoke LAMMPS, or you can +use the :doc:`suffix ` command in your input script. See the :doc:`Speed packages ` doc page for more instructions on how to use the accelerated styles effectively. diff --git a/doc/src/minimize.rst b/doc/src/minimize.rst index 81a454ead3..7ed77994af 100644 --- a/doc/src/minimize.rst +++ b/doc/src/minimize.rst @@ -38,24 +38,25 @@ be in local potential energy minimum. More precisely, the configuration should approximate a critical point for the objective function (see below), which may or may not be a local minimum. -The minimization algorithm used is set by the -:doc:`min\_style ` command. Other options are set by the -:doc:`min\_modify ` command. Minimize commands can be -interspersed with :doc:`run ` commands to alternate between -relaxation and dynamics. The minimizers bound the distance atoms move -in one iteration, so that you can relax systems with highly overlapped -atoms (large energies and forces) by pushing the atoms off of each -other. +The minimization algorithm used is set by the :doc:`min\_style +` command. Other options are set by the :doc:`min\_modify +` command. Minimize commands can be interspersed with +:doc:`run ` commands to alternate between relaxation and +dynamics. The minimizers bound the distance atoms move in one +iteration, so that you can relax systems with highly overlapped atoms +(large energies and forces) by pushing the atoms off of each other. Alternate means of relaxing a system are to run dynamics with a small or :doc:`limited timestep `. Or dynamics can be run using :doc:`fix viscous ` to impose a damping force that -slowly drains all kinetic energy from the system. The :doc:`pair\_style soft ` potential can be used to un-overlap atoms while -running dynamics. +slowly drains all kinetic energy from the system. The +:doc:`pair\_style soft ` potential can be used to +un-overlap atoms while running dynamics. Note that you can minimize some atoms in the system while holding the -coordinates of other atoms fixed by applying :doc:`fix setforce ` to the other atoms. See a fuller -discussion of using fixes while minimizing below. +coordinates of other atoms fixed by applying :doc:`fix setforce +` to the other atoms. See a fuller discussion of using +fixes while minimizing below. The :doc:`minimization styles ` *cg*\ , *sd*\ , and *hftn* involves an outer iteration loop which sets the search direction along @@ -69,15 +70,15 @@ backtracking method is described in Nocedal and Wright's Numerical Optimization (Procedure 3.1 on p 41). The :doc:`minimization styles ` *quickmin*\ , *fire* and -*fire2* perform damped dynamics using an Euler integration step. -Thus they require a :doc:`timestep ` be defined. +*fire2* perform damped dynamics using an Euler integration step. Thus +they require a :doc:`timestep ` be defined. .. note:: The damped dynamic minimizers use whatever timestep you have - defined via the :doc:`timestep ` command. Often they will - converge more quickly if you use a timestep about 10x larger than you - would normally use for dynamics simulations. + defined via the :doc:`timestep ` command. Often they + will converge more quickly if you use a timestep about 10x larger + than you would normally use for dynamics simulations. ---------- @@ -90,13 +91,15 @@ coordinates: .. image:: Eqs/min_energy.jpg :align: center -where the first term is the sum of all non-bonded :doc:`pairwise interactions ` including :doc:`long-range Coulombic interactions `, the 2nd through 5th terms are -:doc:`bond `, :doc:`angle `, -:doc:`dihedral `, and :doc:`improper ` -interactions respectively, and the last term is energy due to -:doc:`fixes ` which can act as constraints or apply force to atoms, -such as through interaction with a wall. See the discussion below about -how fix commands affect minimization. +where the first term is the sum of all non-bonded :doc:`pairwise +interactions ` including :doc:`long-range Coulombic +interactions `, the 2nd through 5th terms are :doc:`bond +`, :doc:`angle `, :doc:`dihedral +`, and :doc:`improper ` interactions +respectively, and the last term is energy due to :doc:`fixes ` +which can act as constraints or apply force to atoms, such as through +interaction with a wall. See the discussion below about how fix +commands affect minimization. The starting point for the minimization is the current configuration of the atoms. @@ -126,9 +129,9 @@ The minimization procedure stops if any of several criteria are met: .. note:: You can also use the :doc:`fix halt ` command to specify - a general criterion for exiting a minimization, that is a calculation - performed on the state of the current system, as defined by an - :doc:`equal-style variable `. + a general criterion for exiting a minimization, that is a + calculation performed on the state of the current system, as + defined by an :doc:`equal-style variable `. For the first criterion, the specified energy tolerance *etol* is unitless; it is met when the energy change between successive @@ -163,8 +166,8 @@ freedom, such as from the :doc:`fix box/relax ` command. Following minimization, a statistical summary is printed that lists which convergence criterion caused the minimizer to stop, as well as -information about the energy, force, final line search, and -iteration counts. An example is as follows: +information about the energy, force, final line search, and iteration +counts. An example is as follows: .. parsed-literal:: diff --git a/src/min_fire2.cpp b/src/min_fire2.cpp index 6629e64c15..c34ccc57d2 100644 --- a/src/min_fire2.cpp +++ b/src/min_fire2.cpp @@ -11,6 +11,11 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Contributing authors: Julien Guénolé, RWTH Aachen University and + Erik Bitzek, FAU Erlangen-Nuernberg +------------------------------------------------------------------------- */ + #include #include "min_fire2.h" #include "universe.h" @@ -55,7 +60,6 @@ void MinFire2::init() alpha = alpha0; last_negative = ntimestep_start = update->ntimestep; vdotf_negatif = 0; - } /* ---------------------------------------------------------------------- */ @@ -260,6 +264,7 @@ int MinFire2::iterate(int maxiter) x[i][2] -= 0.5 * dtv * v[i][2]; } } + for (int i = 0; i < nlocal; i++) v[i][0] = v[i][1] = v[i][2] = 0.0; } @@ -273,6 +278,7 @@ int MinFire2::iterate(int maxiter) vmax = MAX(vmax,fabs(v[i][2])); if (dtvone*vmax > dmax) dtvone = dmax/vmax; } + MPI_Allreduce(&dtvone,&dtv,1,MPI_DOUBLE,MPI_MIN,world); // min dtv over replicas, if necessary @@ -333,7 +339,7 @@ int MinFire2::iterate(int maxiter) // Velocity Verlet integration - }else if (integrator == 1) { + } else if (integrator == 1) { dtf = 0.5 * dtv * force->ftm2v; @@ -391,7 +397,7 @@ int MinFire2::iterate(int maxiter) // Standard Euler integration - }else if (integrator == 3) { + } else if (integrator == 3) { dtf = dtv * force->ftm2v; @@ -430,7 +436,6 @@ int MinFire2::iterate(int maxiter) eprevious = ecurrent; ecurrent = energy_force(0); neval++; - } // energy tolerance criterion