Geometry optimization input parameters
New input file: input.yaml
To trigger a geometry optimization you have to add a geopt section to the input.yaml file. For example:
geopt: method: SQNM ncount_cluster_x: 200 frac_fluct: 2.0 forcemax: 1.E-4 betax: 1
The above example uses the stabilized quasi newton minimizer (SQNM).^{[1]} Due to efficiency and stability, SQNM is the recommended minimizer within the BigDFT suite.
The geometry optimizations can be monitored with the data/geopt.mon file.
Available minimization methods
- SDCG: A combination of Steepest Descent and Conjugate Gradient
- VSSD: Variable Stepsize Steepest Descent method
- LBFGS: Limited-memory BFGS
- BFGS: Broyden-Fletcher-Goldfarb-Shanno
- PBFGS: Same as BFGS with an initial Hessian obtained from a force field
- DIIS: Direct inversion of iterative subspace
- FIRE: Fast Inertial Relaxation Engine as described by Bitzek et al.
- SQNM: Stabilized quasi-Newton minimzer
SQNM
The SQNM minimizer supports additional keywords. The following example shows all the available keywords for SQNM.
geopt: method: SQNM ncount_cluster_x: 200 frac_fluct: 2.0 forcemax: 1.E-4 betax: 1 biomode: no beta_stretchx: 5e-1 nhistx: 10 cutoffratio: 1e-4 steepthresh: .5 trustr: 0.5 maxrise: 1e-6
Keywords
- ncount_cluster_x: Maximum number of force evaluations to be used for the geometry optimization.
- frac_fluct, forcemax: Convergence criteria for the geometry optimization. The geometry optimization stops either if the norm of the individual forces acting on any atom in the system is smaller than forcemax or if the forces get noisy. Noise is present because of the underlying integration grid. The parameter frac_fluct specifies how small the forces should become compared to the noise level to stop the geometry optimization. Values in between 1 and 10 is are reasonable for this parameter. A value of 2 means that the geometry optimization will stop when the largest atomic force is comparable to the 2 times the average noise in the forces themselves. For values of frac_fluct smaller than 1 one can under certain circumstances obtain better relaxed geometries but one risks that the geometry optimization will not converge since the forces are too noisy. In such a case one should closely monitor the progress of the geometry optimization by looking at the posout.* files which are written at each step of the geometry optimization and at the data/geopt.mon file.
- betax: This is the step size for the geometry optimization. For SQNM a value of 1 works fine for many systems. For good efficiency, this parameter must be tuned!
- biomode: yes or no. yes activates additional preconditioning for bio-molecules, like, for example, alanine dipeptide.
- beta_stretchx: only of any relevance if biomode=yes. beta_stretchx defines the step size parallel to bond directions
- nhistx: length of the history of gradients and positions used to build an approximate hessian. 10 works fine for most systems. Typically 5 <=nhistx<=20.
- cutoffratio: Determines which displacement vectors form the significant subspace. 10e-4 works well for most systems and typical production accuracy of forces.
- steepthresh If the maximum force component on any atom is larger than steepthresh, steepest descent instead of SQNM is used.
- trustr: Trustradius. Atoms are not allowed to be displaced by more than trustr within one iteration.
- maxrise: Energy is not allowed to rise more than by maxrise within one iteration.
For more information on the meaning of these keywords it is referred to the SQNM article.^{[1]}
Example geopt.mon file
#------------------- Geopt file opened, name: ./data/geopt.mon, timestamp: 2016-10-14 02:36:27.724 # COUNT IT GEOPT_METHOD ENERGY DIFF FMAX FNRM FRAC*FLUC FLUC ADD. INFO 0 0 GEOPT_SQNM 0.00000000000000E+00 0.00E+00 1.510E-02 2.30E-02 3.91E-04 1.95E-04 beta=8.00E-01 dim=000 maxd=0.0E+00 dsplr=0.00000E+00 dsplp=0.00000E+00 1 1 GEOPT_SQNM -5.56488295298244E+02 -5.56E+02 1.199E-02 1.82E-02 3.86E-04 1.93E-04 beta=8.00E-01 dim=000 maxd=1.2E-02 dsplr=1.83971E-02 dsplp=1.83971E-02 2 2 GEOPT_SQNM -5.56488436457742E+02 -1.41E-04 1.019E-02 1.49E-02 3.90E-04 1.95E-04 beta=8.00E-01 dim=001 maxd=9.7E-03 dsplr=3.29422E-02 dsplp=3.29422E-02 3 3 GEOPT_SQNM -5.56488685698921E+02 -2.49E-04 2.867E-03 5.74E-03 3.96E-04 1.98E-04 beta=8.00E-01 dim=002 maxd=1.9E-02 dsplr=7.16202E-02 dsplp=7.16202E-02 4 4 GEOPT_SQNM -5.56488744266590E+02 -5.86E-05 3.171E-03 5.60E-03 4.14E-04 2.07E-04 beta=8.00E-01 dim=003 maxd=1.5E-02 dsplr=1.02664E-01 dsplp=1.02664E-01 5 5 GEOPT_SQNM -5.56488777609738E+02 -3.33E-05 1.638E-03 2.87E-03 4.28E-04 2.14E-04 beta=8.00E-01 dim=003 maxd=7.4E-03 dsplr=1.16061E-01 dsplp=1.16061E-01 6 6 GEOPT_SQNM -5.56488802430900E+02 -2.48E-05 2.307E-03 2.86E-03 4.44E-04 2.22E-04 beta=8.80E-01 dim=004 maxd=8.2E-03 dsplr=1.31036E-01 dsplp=1.31036E-01 7 7 GEOPT_SQNM -5.56488812886668E+02 -1.05E-05 1.556E-03 2.16E-03 4.58E-04 2.29E-04 beta=9.68E-01 dim=004 maxd=4.3E-03 dsplr=1.37242E-01 dsplp=1.37242E-01 8 8 GEOPT_SQNM -5.56488824888933E+02 -1.20E-05 1.231E-03 1.95E-03 4.71E-04 2.36E-04 beta=1.06E+00 dim=005 maxd=9.2E-03 dsplr=1.50336E-01 dsplp=1.50336E-01 9 9 GEOPT_SQNM -5.56488843393915E+02 -1.85E-05 3.905E-04 8.05E-04 4.79E-04 2.39E-04 beta=1.17E+00 dim=006 maxd=1.3E-02 dsplr=1.69873E-01 dsplp=1.69873E-01 10 10 GEOPT_SQNM -5.56488845740426E+02 -2.35E-06 3.063E-04 6.66E-04 4.83E-04 2.42E-04 beta=1.29E+00 dim=006 maxd=4.5E-03 dsplr=1.76899E-01 dsplp=1.76899E-01 11 11 GEOPT_SQNM -5.56488848055483E+02 -2.32E-06 3.757E-04 6.05E-04 4.86E-04 2.43E-04 beta=1.10E+00 dim=007 maxd=2.7E-03 dsplr=1.81566E-01 dsplp=1.81566E-01 12 12 GEOPT_SQNM -5.56488848484636E+02 -4.29E-07 3.515E-04 5.53E-04 4.87E-04 2.44E-04 beta=1.20E+00 dim=007 maxd=8.3E-04 dsplr=1.83171E-01 dsplp=1.83171E-01 13 13 GEOPT_SQNM -5.56488850308366E+02 -1.82E-06 2.993E-04 5.13E-04 4.88E-04 2.44E-04 beta=1.33E+00 dim=008 maxd=3.8E-03 dsplr=1.89199E-01 dsplp=1.89199E-01 14 14 GEOPT_SQNM -5.56488850755728E+02 -4.47E-07 3.185E-04 5.05E-04 4.89E-04 2.45E-04 beta=1.13E+00 dim=007 maxd=1.1E-03 dsplr=1.90874E-01 dsplp=1.90874E-01 SQNM converged at iteration 14. Needed bigdft calls: 14
- beta: the current step size. Cannot get smaller than betax.
- dim: the current dimension of the significant subspace. Cannot get larger than nhistx.
- maxd: The largest displacement of any atom in this iteration.
- displr: Integrated displacement of the entire system including rejected steps.
- displp: Integrated physical displacement of the entire system (that is, excluding rejected steps)
Old input files
This section describes the old input.geopt file format. This format is deprecated and should not be used any more. This file contains all the parameters required for a geometry optimization. It is named input.geopt.
Example
BFGS # Geometry optimisation method 100 # Maximum number of force evaluations 1e-1 5e-4 # fract_fluct,forcemax 0.0 # random displacement amplitude 4.0 # Stepsize for the geometry optimisation
Entry descriptions
All lines are mandatory. Follows a description line by line:
The optimization method
- geopt_approach: This character string specifies the method used for the geometry optimization.
- VSSD: Variable Stepsize Steepest Descent method;
- SDCG: A combination of Steepest Descent and Conjugate Gradient;
- LBFGS: Limited Memory BFGS;
- BFGS: Preconditioned steepest descent with energy feedback. A preconditioning matrix is build up according to the BFGS algorithm. The initial Hessian is a diagonal matrix where the diagonal elements are the inverse of the step size. This method is usually the most efficient one.
- PBFGS: Same as BFGS except that an initial Hessian is obtained from a force field. Force field parameters are only available for systems consisting of H,C,N,O. For such systems the method is the most efficient one in general
- FIRE: The Fast Inertial Relaxation Engine as described by Bitzek et al.^{[2]} with optimized parameters. An additional line is required in the input.geopt file with two parameters: initial time step dtinit and the maximum time step dtmax. A good starting point to estimate dtmax is by using dtmax=Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): 0.25\pi \sqrt{ \mathrm{betax_{VSSD}}}
and dtinit=Failed to parse (PNG conversion failed; check for correct installation of latex and dvipng (or dvips + gs + convert)): 0.5
dtmax, and gradually increasing dtmax. betax is not used
- AB6MD: The molecular dynamic routines from ABINIT 6.
- DIIS: a Direct Inversion of Iterative Subspace. In that case, an additional parameter is required with the history length (default is 4) on the line of betax.
The basic parameters
- ncount_cluster_x: Maximum number of force evaluations to be used for the geometry optimization.
- frac_fluct, forcemax: Convergence criteria for the geometry optimization. The geometry optimization stops either if the norm of the individual forces acting on any atom in the system is smaller than forcemax or if the forces get noisy. Noise is present because of the underlying integration grid. The parameter frac_fluct specifies how small the forces should become compared to the noise level to stop the geometry optimization. Values in between 1 and 10 is are reasonable for this parameter. A value of 2 means that the geometry optimization will stop when the largest atomic force is comparable to the 2 times the average noise in the forces themselves. For values of frac_fluct smaller than 1 one can under certain circumstances obtain better relaxed geometries but one risks that the geometry optimization will not converge since the forces are too noisy. In such a case one should closely monitor the progress of the geometry optimization by looking at the posout.* files which are written at each step of the geometry optimization and at the data/geopt.mon file.
- randdis: This parameter allows to add random displacements with amplitude randdis to the atomic positions in the input file posinp. This can for instance be useful to break degeneracies (which would lead to convergence problems in the wavefunction optimization) in highly symmetric structures.
- betax: This is the stepsize for the geometry optimization. This stepsize is system dependent and it has therefore to be determined for each system. If the VSSD method is used one can start with a small stepsize of around 1 and VSSD will suggest than a better value for {betax} in the last line of the data/geopt.mon file. Whether the stepsize is correct can also be seen from the data/geopt.mon output of the SDCG method. In this case the average stepsize in terms of betax should be around 4 (after a brief initial period where it is around 8). In contrast to the SDCG method, the BFGS method is not very sensitive to the correct stepsize but nevertheless one should try to find reasonable values also in this case.
The specific parameters, when geopt_approach = AB6MD
In case of geopt_approach = AB6MD, the betax line should be replaced ionmov. It is an integer value as described in the ABINIT manual. Possible values are:
- 6: simple velocity-Verlet molecular dynamic.
- 7: quenched molecular dynamic, when the scalar product force / velocity becomes negative, the velocity is set to zero. The force criterion is tested at each step.
- 8: Nose-Hoover thermostat.
- 9: Langevin dynamic (adding a friction force and a Gaussian random force on atoms).
- 12: Isokinetic ensemble molecular dynamics. The equation of motion of the ions in contact with a thermostat are solved with the algorithm proposed by Zhang^{[3]}, as worked out by Minary et al.^{[4]}. The conservation of the kinetic energy is obtained within machine precision, at each step.
- 13: Iosthermal/isenthalpic ensemble. The equation of motion of the ions in contact with a thermostat and a barostat are solved with the algorithm proposed by Martyna, Tuckermann Tobias and Klein^{[5]}.
Additional lines are also mandatory:
- dtion: the time step for molecular dynamic, in atomic time units. One atomic time unit is 2.418884e-17 seconds, which is the value of Planck's constant in Hartree×sec. The following lines depend on the choosen value for ionmov.
- if ionmov = 6: one line containing the initial temperature in kelvin. If negative, the initial velocities are all zero. If positive, random speeds are chosen to match the given temperature (NOT IMPLEMENTED YET!).
- if ionmov > 7: one line containing two temperatures in kelvin. When different, the temperature is linearly change at each geometry step to go from the first value to the second.
- if ionmov = 8: one additional line containing the thermostat inertia coefficient for Nose-Hoover dynamic.
- if ionmov = 9: two additional lines containing first the friction coefficient and then a value in Bohr corresponding at a distance where the atoms can bounce on (see the ABINIT documentation for further details).
- if ionmov = 13: several additional lines containing first the number of thermostats in the chain of thermostats. Then a line with the mass of each thermostat in the chain. And finally a line with two values for the barostat mass, depending on optcell value (NOT IMPLEMENTED YET!).
- ↑ ^{1.0} ^{1.1} J. Chem. Phys. 142, 034112 (2015)
- ↑ Phys. Rev. Let. 97, 170201 (2006)
- ↑ J. Chem. Phys. 106, 6102 (1997)
- ↑ J. Chem. Phys. 188, 2510 (2003)
- ↑ Mol. Phys., 1996, p. 1117