diff --git a/doc/minimize.html b/doc/minimize.html index 4b2dfbdbe..2c88f4df1 100644 --- a/doc/minimize.html +++ b/doc/minimize.html @@ -1,228 +1,228 @@ <HTML> <CENTER><A HREF = "http://lammps.sandia.gov">LAMMPS WWW Site</A> - <A HREF = "Manual.html">LAMMPS Documentation</A> - <A HREF = "Section_commands.html#comm">LAMMPS Commands</A> </CENTER> <HR> <H3>minimize command </H3> <P><B>Syntax:</B> </P> <PRE>minimize etol ftol maxiter maxeval </PRE> <UL><LI>etol = stopping tolerance for energy (unitless) <LI>ftol = stopping tolerance for force (force units) <LI>maxiter = max iterations of minimizer <LI>maxeval = max number of total force/energy evaluations </UL> <P><B>Examples:</B> </P> <PRE>minimize 1.0e-4 1.0e-6 100 1000 minimize 0.0 1.0e-8 1000 100000 </PRE> <P><B>Description:</B> </P> <P>Perform an energy minimization of the system, by iteratively adjusting atom coordinates. Iterations are terminated when one of the stopping criteria is satisfied. At that point the configuration will hopefully be in local potential energy minimum. Mathematically speaking, the configuration should approximate a critical point for the objective function (see below), which may or may not be a local minimum. </P> <P>The minimization algorithm used is set by the <A HREF = "min_style.html">min_style</A> command. Other options are set by the <A HREF = "min_modify.html">min_modify</A> command. Minimize commands can be interspersed with <A HREF = "run.html">run</A> 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. </P> <P>Alternate means of relaxing a system are to run dynamics with a small or <A HREF = "fix_nve_limit.html">limited timestep</A>. Or dynamics can be run using <A HREF = "fix_viscous.html">fix viscous</A> to impose a damping force that slowly drains all kinetic energy from the system. The <A HREF = "pair_soft.html">pair_style soft</A> potential can be used to un-overlap atoms while running dynamics. </P> <P>A minimization involves an outer iteration loop which sets the search direction along which atom coordinates are changed. An inner iteration is then performed using a line search algorithm. The line search typically evaluates forces and energies several times to set new coordinates. Currently, a backtracking algorithm is used which may not be optimal in terms of the number of force evaulations performed, but appears to be more robust than previous line searches we've tried. The backtracking method is described in Nocedal and Wright's Numerical Optimization (Procedure 3.1 on p 41). </P> <P>The objective function being minimized is the total potential energy of the system as a function of the N atom coordinates: </P> <CENTER><IMG SRC = "Eqs/min_energy.jpg"> </CENTER> <P>where the first term is the sum of all non-bonded <A HREF = "pair_style.html">pairwise interactions</A> including <A HREF = "kspace_style.html">long-range Coulombic interactions</A>, the 2nd thru 5th terms are <A HREF = "bond_style.html">bond</A>, <A HREF = "angle_style.html">angle</A>, <A HREF = "dihedral_style.html">dihedral</A>, and <A HREF = "improper_style.html">improper</A> interactions respectively, and the last term is energy due to <A HREF = "fix.html">fixes</A> which can act as constraints or apply force to atoms, such as thru interaction with a wall. See the discussion below about -which fix commands affect minimization. +how fix commands affect minimization. </P> <P>The starting point for the minimization is the current configuration of the atoms. </P> <HR> <P>The minimization procedure stops if any of several criteria are met: </P> <UL><LI>the change in energy between outer iterations is less than <I>etol</I> <LI>the 2-norm (length) of the global force vector is less than the <I>ftol</I> <LI>the line search fails because the step distance backtracks to 0.0 <LI>the number of outer iterations exceeds <I>maxiter</I> <LI>the number of total force evaluations exceeds <I>maxeval</I> </UL> <P>For the first criterion, the specified energy tolerance <I>etol</I> is unitless; it is met when the energy change between successive iterations divided by the energy magnitude is less than or equal to the tolerance. For example, a setting of 1.0e-4 for <I>etol</I> means an energy tolerance of one part in 10^4. </P> <P>For the second criterion, the specified force tolerance <I>ftol</I> is in force units, since it is the length of the global force vector for all atoms, e.g. a vector of size 3N for N atoms. Since many of the components will be near zero after minimization, you can think of <I>ftol</I> as an upper bound on the final force on any component of any atom. For example, a setting of 1.0e-4 for <I>ftol</I> means no x, y, or z component of force on any atom will be larger than 1.0e-4 (in force units) after minimization. </P> <P>Either or both of the <I>etol</I> and <I>ftol</I> values can be set to 0.0, in which case some other criterion will terminate the minimization. </P> <P>During a minimization, the outer iteration count is treated as a timestep. Output is triggered by this timestep, e.g. thermodynamic output or dump and restart files. </P> <P>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 and iteration counts. An example is as follows: </P> <PRE>Minimization stats: Stopping criterion = max iterations Energy initial, next-to-last, final = -0.626828169302 -2.82642039062 -2.82643549739 Force two-norm initial, final = 2052.1 91.9642 Force max component initial, final = 346.048 9.78056 Final line search alpha, max atom move = 2.23899e-06 2.18986e-05 Iterations, force evaluations = 2000 12724 </PRE> <P>The 3 energy values are for before and after the minimization and on the next-to-last iteration. This is what the <I>etol</I> parameter checks. </P> <P>The two-norm force values are the length of the global force vector before and after minimization. This is what the <I>ftol</I> parameter checks. </P> <P>The max-component force values are the absolute value of the largest component (x,y,z) in the global force vector. </P> <P>The alpha parameter for the line-search, when multiplied by the max force component (on the last iteration), gives the max distance any atom moved during the last iteration. Alpha will be 0.0 if the line search could not reduce the energy. Even if alpha is non-zero, if the "max atom move" distance is tiny compared to typical atom coordinates, then it is possible the last iteration effectively caused no atom movement and thus the evaluated energy did not change and the minimizer terminated. Said another way, even with non-zero forces, it's possible the effect of those forces is to move atoms a distance less than machine precision, so that the energy cannot be further reduced. </P> <P>The iterations and force evaluation values are what is checked by the <I>maxiter</I> and <I>maxeval</I> parameters. </P> <HR> <P>IMPORTANT NOTE: It is highly recommended that you use a <A HREF = "pair_style.html">pair style</A> that goes to 0.0 at the cutoff distance when performing minimization (even if you later change it when running dynamics). If this is not done, the total energy of the system will have discontinuities when the relative distance between any pair of atoms changes from cutoff+epsilon to cutoff-epsilon and the minimizer may behave poorly. </P> <P>Note that a cutoff Lennard-Jones potential (and others) can be shifted so that its energy is 0.0 at the cutoff via the <A HREF = "pair_modify">pair_modify</A> command. See the doc pages for inidividual <A HREF = "pair_style.html">pair styles</A> for details. Note that Coulombic potentials always have a cutoff, unless versions with a long-range component are used (e.g. <A HREF = "pair_lj.html">pair_style lj/cut/coul/long</A>). The CHARMM potentials go to 0.0 at the cutoff (e.g. <A HREF = "pair_charmm.html</I>">pair_style lj/charmm/coul/charmm</A>, as do the GROMACS potentials (e.g. <A HREF = "pair_gromacs.html</I>">pair_style lj/gromacs</A>. </P> <P>If a soft potential (<A HREF = "pair_soft.html">pair_style soft</A>) is used the Astop value is used for the prefactor (no time dependence). </P> <P>Only some fixes (typically those that apply force constraints) are invoked during minimization. See the doc pages for individual <A HREF = "fix.html">fix</A> commands to see which ones are relevant. </P> <P>IMPORTANT NOTE: Some fixes which are invoked during minimization have an associated potential energy. For that energy to be included in the total potential energy of the system (the quantity being minimized), you MUST enable the <A HREF = "fix_modify.html">fix_modify</A> <I>energy</I> option for that fix. The doc pages for individual <A HREF = "fix.html">fix</A> commands specify if this should be done. </P> <HR> <P><B>Restrictions:</B> </P> <P>Features that are not yet implemented are listed here, in case someone knows how they could be coded: </P> <P>It is an error to use <A HREF = "fix_shake.html">fix shake</A> with minimization because it turns off bonds that should be included in the potential energy of the system. The effect of a fix shake can be approximated during a minimization by using stiff spring constants for the bonds and/or angles that would normally be constrained by the SHAKE algorithm. </P> <P><A HREF = "fix_rigid.html">Fix rigid</A> is also not supported by minimization. It is not an error to have it defined, but the energy minimization will not keep the defined body(s) rigid during the minimization. Note that if bonds, angles, etc internal to a rigid body have been turned off (e.g. via <A HREF = "neigh_modify.html">neigh_modify exclude</A>), they will not contribute to the potential energy which is probably not what is desired. </P> <P>The volume of the simulation domain is not allowed to change during a minimization. Ideally we would allow a fix such as <I>npt</I> to impose an external pressure that would be included in the minimization (i.e. allow the box dimensions to change), but this has not yet been implemented. </P> <P>Pair potentials that produce torque on a particle (e.g. <A HREF = "pair_gran.html">granular potentials</A> or the <A HREF = "pair_gayberne.html">GayBerne potential</A> for ellipsoidal particles) are not relaxed by a minimization. More specifically, radial relaxations are induced, but no rotations are induced by a minimization, so such a system will not fully relax. </P> <P><B>Related commands:</B> </P> <P><A HREF = "min_modify.html">min_modify</A>, <A HREF = "min_style.html">min_style</A>, <A HREF = "run_style.html">run_style</A> </P> <P><B>Default:</B> none </P> </HTML> diff --git a/doc/minimize.txt b/doc/minimize.txt index bf9fcd0e3..7085e21d3 100644 --- a/doc/minimize.txt +++ b/doc/minimize.txt @@ -1,223 +1,223 @@ "LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c :link(lws,http://lammps.sandia.gov) :link(ld,Manual.html) :link(lc,Section_commands.html#comm) :line minimize command :h3 [Syntax:] minimize etol ftol maxiter maxeval :pre etol = stopping tolerance for energy (unitless) ftol = stopping tolerance for force (force units) maxiter = max iterations of minimizer maxeval = max number of total force/energy evaluations :ul [Examples:] minimize 1.0e-4 1.0e-6 100 1000 minimize 0.0 1.0e-8 1000 100000 :pre [Description:] Perform an energy minimization of the system, by iteratively adjusting atom coordinates. Iterations are terminated when one of the stopping criteria is satisfied. At that point the configuration will hopefully be in local potential energy minimum. Mathematically speaking, 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 "min_style"_min_style.html command. Other options are set by the "min_modify"_min_modify.html command. Minimize commands can be interspersed with "run"_run.html 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 "limited timestep"_fix_nve_limit.html. Or dynamics can be run using "fix viscous"_fix_viscous.html to impose a damping force that slowly drains all kinetic energy from the system. The "pair_style soft"_pair_soft.html potential can be used to un-overlap atoms while running dynamics. A minimization involves an outer iteration loop which sets the search direction along which atom coordinates are changed. An inner iteration is then performed using a line search algorithm. The line search typically evaluates forces and energies several times to set new coordinates. Currently, a backtracking algorithm is used which may not be optimal in terms of the number of force evaulations performed, but appears to be more robust than previous line searches we've tried. The backtracking method is described in Nocedal and Wright's Numerical Optimization (Procedure 3.1 on p 41). The objective function being minimized is the total potential energy of the system as a function of the N atom coordinates: :c,image(Eqs/min_energy.jpg) where the first term is the sum of all non-bonded "pairwise interactions"_pair_style.html including "long-range Coulombic interactions"_kspace_style.html, the 2nd thru 5th terms are "bond"_bond_style.html, "angle"_angle_style.html, "dihedral"_dihedral_style.html, and "improper"_improper_style.html interactions respectively, and the last term is energy due to "fixes"_fix.html which can act as constraints or apply force to atoms, such as thru interaction with a wall. See the discussion below about -which fix commands affect minimization. +how fix commands affect minimization. The starting point for the minimization is the current configuration of the atoms. :line The minimization procedure stops if any of several criteria are met: the change in energy between outer iterations is less than {etol} the 2-norm (length) of the global force vector is less than the {ftol} the line search fails because the step distance backtracks to 0.0 the number of outer iterations exceeds {maxiter} the number of total force evaluations exceeds {maxeval} :ul For the first criterion, the specified energy tolerance {etol} is unitless; it is met when the energy change between successive iterations divided by the energy magnitude is less than or equal to the tolerance. For example, a setting of 1.0e-4 for {etol} means an energy tolerance of one part in 10^4. For the second criterion, the specified force tolerance {ftol} is in force units, since it is the length of the global force vector for all atoms, e.g. a vector of size 3N for N atoms. Since many of the components will be near zero after minimization, you can think of {ftol} as an upper bound on the final force on any component of any atom. For example, a setting of 1.0e-4 for {ftol} means no x, y, or z component of force on any atom will be larger than 1.0e-4 (in force units) after minimization. Either or both of the {etol} and {ftol} values can be set to 0.0, in which case some other criterion will terminate the minimization. During a minimization, the outer iteration count is treated as a timestep. Output is triggered by this timestep, e.g. thermodynamic output or dump and restart files. 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 and iteration counts. An example is as follows: Minimization stats: Stopping criterion = max iterations Energy initial, next-to-last, final = -0.626828169302 -2.82642039062 -2.82643549739 Force two-norm initial, final = 2052.1 91.9642 Force max component initial, final = 346.048 9.78056 Final line search alpha, max atom move = 2.23899e-06 2.18986e-05 Iterations, force evaluations = 2000 12724 :pre The 3 energy values are for before and after the minimization and on the next-to-last iteration. This is what the {etol} parameter checks. The two-norm force values are the length of the global force vector before and after minimization. This is what the {ftol} parameter checks. The max-component force values are the absolute value of the largest component (x,y,z) in the global force vector. The alpha parameter for the line-search, when multiplied by the max force component (on the last iteration), gives the max distance any atom moved during the last iteration. Alpha will be 0.0 if the line search could not reduce the energy. Even if alpha is non-zero, if the "max atom move" distance is tiny compared to typical atom coordinates, then it is possible the last iteration effectively caused no atom movement and thus the evaluated energy did not change and the minimizer terminated. Said another way, even with non-zero forces, it's possible the effect of those forces is to move atoms a distance less than machine precision, so that the energy cannot be further reduced. The iterations and force evaluation values are what is checked by the {maxiter} and {maxeval} parameters. :line IMPORTANT NOTE: It is highly recommended that you use a "pair style"_pair_style.html that goes to 0.0 at the cutoff distance when performing minimization (even if you later change it when running dynamics). If this is not done, the total energy of the system will have discontinuities when the relative distance between any pair of atoms changes from cutoff+epsilon to cutoff-epsilon and the minimizer may behave poorly. Note that a cutoff Lennard-Jones potential (and others) can be shifted so that its energy is 0.0 at the cutoff via the "pair_modify"_pair_modify command. See the doc pages for inidividual "pair styles"_pair_style.html for details. Note that Coulombic potentials always have a cutoff, unless versions with a long-range component are used (e.g. "pair_style lj/cut/coul/long"_pair_lj.html). The CHARMM potentials go to 0.0 at the cutoff (e.g. "pair_style lj/charmm/coul/charmm"_pair_charmm.html}, as do the GROMACS potentials (e.g. "pair_style lj/gromacs"_pair_gromacs.html}. If a soft potential ("pair_style soft"_pair_soft.html) is used the Astop value is used for the prefactor (no time dependence). Only some fixes (typically those that apply force constraints) are invoked during minimization. See the doc pages for individual "fix"_fix.html commands to see which ones are relevant. IMPORTANT NOTE: Some fixes which are invoked during minimization have an associated potential energy. For that energy to be included in the total potential energy of the system (the quantity being minimized), you MUST enable the "fix_modify"_fix_modify.html {energy} option for that fix. The doc pages for individual "fix"_fix.html commands specify if this should be done. :line [Restrictions:] Features that are not yet implemented are listed here, in case someone knows how they could be coded: It is an error to use "fix shake"_fix_shake.html with minimization because it turns off bonds that should be included in the potential energy of the system. The effect of a fix shake can be approximated during a minimization by using stiff spring constants for the bonds and/or angles that would normally be constrained by the SHAKE algorithm. "Fix rigid"_fix_rigid.html is also not supported by minimization. It is not an error to have it defined, but the energy minimization will not keep the defined body(s) rigid during the minimization. Note that if bonds, angles, etc internal to a rigid body have been turned off (e.g. via "neigh_modify exclude"_neigh_modify.html), they will not contribute to the potential energy which is probably not what is desired. The volume of the simulation domain is not allowed to change during a minimization. Ideally we would allow a fix such as {npt} to impose an external pressure that would be included in the minimization (i.e. allow the box dimensions to change), but this has not yet been implemented. Pair potentials that produce torque on a particle (e.g. "granular potentials"_pair_gran.html or the "GayBerne potential"_pair_gayberne.html for ellipsoidal particles) are not relaxed by a minimization. More specifically, radial relaxations are induced, but no rotations are induced by a minimization, so such a system will not fully relax. [Related commands:] "min_modify"_min_modify.html, "min_style"_min_style.html, "run_style"_run_style.html [Default:] none