From 3274ff38b0ae5be5c0eda99c7b728fd626b8a7d2 Mon Sep 17 00:00:00 2001 From: HE Zilong <159878975+he-zilong@users.noreply.github.com> Date: Tue, 26 May 2026 17:05:23 +0800 Subject: [PATCH] Expose LBFGS print control as input keyword, increase default verbosity, and expand document for optimization (#5274) --- .../optimization/geometry_and_cell_opt.md | 305 ++++++++++++++---- src/input_constants.F | 7 + src/motion/cp_lbfgs_geo.F | 11 +- src/motion/gopt_param_types.F | 19 +- src/start/input_cp2k_motion.F | 44 ++- 5 files changed, 303 insertions(+), 83 deletions(-) diff --git a/docs/methods/optimization/geometry_and_cell_opt.md b/docs/methods/optimization/geometry_and_cell_opt.md index ad126314bd..99f6f2370b 100644 --- a/docs/methods/optimization/geometry_and_cell_opt.md +++ b/docs/methods/optimization/geometry_and_cell_opt.md @@ -32,10 +32,20 @@ only internal atomic coordinates should relax. Use `CELL_OPT` when the equilibrium cell volume, shape, pressure, or residual stress is part of the question. It is most common for bulk crystals, strained solids, two-dimensional materials with -relaxed in-plane lattice constants, and structures prepared for fixed-cell molecular dynamics. A -cell optimization is not a replacement for finite-temperature pressure sampling; if the desired -quantity is a thermal average at finite temperature, use an appropriate NPT molecular-dynamics -workflow instead. +relaxed in-plane lattice constants, and structures prepared for fixed-cell molecular dynamics. + +```{note} +An optimization is not a replacement for finite-temperature pressure sampling; if the desired +quantity is a thermal average at finite temperature, use a molecular-dynamics workflow with +appropriate ensemble instead. + +For the usual electronic-structure methods based on the Born-Oppenheimer approximation +(which, by neglecting nuclear motion, provides the concept of "potential-energy surface" +in the first place), there is no involvement of the temperature (kinetic energy of the +nuclei) in an optimization, and the energy that is being optimized does not include +components like zero-point energy and thermal corrections (harmonic or anharmonic). +The notion "optimization at 0 K" is really a misnomer in this regard. +``` ## Basic setup @@ -58,7 +68,7 @@ block: &GEO_OPT TYPE MINIMIZATION OPTIMIZER BFGS - MAX_ITER 200 + MAX_ITER 200 MAX_DR 3.0E-03 RMS_DR 1.5E-03 MAX_FORCE 4.5E-04 @@ -87,7 +97,7 @@ sense as well: &MOTION &CELL_OPT OPTIMIZER BFGS - MAX_ITER 200 + MAX_ITER 200 EXTERNAL_PRESSURE 1.01325E+00 PRESSURE_TOLERANCE 100.0 MAX_DR 3.0E-03 @@ -98,6 +108,14 @@ sense as well: &END MOTION ``` +```{note} +The availability of analytical stress tensor depends on the method for evaluating energy and force. +Some electronic-structure methods may only support numerical stress tensor from finite differences, +which would take quite more time at each iteration of optimization. This can be confirmed by setting +[FORCE_EVAL/PRINT/STRESS_TENSOR](#CP2K_INPUT.FORCE_EVAL.PRINT.STRESS_TENSOR) and seeing if numerical +stress is mentioned in the output. +``` + For `GEO_OPT`, [TYPE](#CP2K_INPUT.MOTION.GEO_OPT.TYPE) selects the optimization target. `MINIMIZATION` searches for a local minimum and is the normal choice. `TRANSITION_STATE` activates transition-state search machinery; it requires additional method-specific settings in @@ -108,30 +126,72 @@ treated as an ordinary minimization. optimization iteration can require more than one force evaluation, depending on the optimizer, line search, and the selected `CELL_OPT` mode. +An optimization may terminate if one of the following events takes place: + +- The convergence criteria (see below) are satisfied before reaching `MAX_ITER`; +- The `MAX_ITER` is reached, regardless of whether the convergence criteria are satisfied or not; +- Some external control for the program is activated, such as hitting the global walltime as + specified by [WALLTIME](#CP2K_INPUT.GLOBAL.WALLTIME), or discovering a file in the working + directory with filename `EXIT` or `EXIT_GEO` that is created by user in case of emergency. + ## Starting structure and cell Start from a chemically and physically sensible structure. Severe close contacts, unrealistic coordination environments, or a badly prepared cell can lead to unstable forces, SCF failures, -unphysical atom displacements, or optimizer steps that are difficult to recover from. +unphysical atom displacements, or optimizer steps that are difficult to recover from. For instance, +"amorphous cell" structures generated by randomly packing molecules and/or atoms in a box can be +extremely challenging for optimization, especially if the chemical bonds between atoms are +significantly scrambled. Same for snapshots from a high-temperature molecular-dynamics simulation in +vaporized or molten conditions with lots of broken chemical bonds. For `CELL_OPT`, the initial cell matters as much as the initial coordinates. The starting volume and -shape should be close enough to the expected structure that the pressure and stress are not -dominated by preparation artefacts. For slabs, two-dimensional materials, or systems with vacuum, do -not relax the vacuum direction unless this is physically intended; constrain the appropriate cell +shape should be close enough to the expected structure and density, such that the pressure and +stress are not dominated by preparation artefacts. + +A visualization of the structure and cell in modelling programs, with the box and the neighboring +periodic images displayed, will be very helpful; neither large vacuous gaps nor crowded cluster of +atoms should occur near the boundary of the box on the directions consistent with the periodicity, +and conversely, sufficient vacuum space on the non-periodic directions is crucial for eliminating +unwanted interactions across the boundary. + +For surface slabs, two- or one-dimensional materials, or systems with vacuum, do not relax the +vacuum direction unless physically intended; it is better to constrain the appropriate cell components or use a fixed-cell optimization after choosing the desired cell. -For variable-cell optimizations, a conventional cell representation with A along the X-axis and B on -the XY plane is expected, while non-standard cell orientations can usually still work. - ```{warning} -An experimental `cif` should not be used as the initial structure without validation; refinement and -careful revision is necessary to make it "computation-ready" during pre-processing: +**Do not use experimental structure blindly.** -- Watch out for crystallographic disorder and fractional occupation of atoms: using the - superposition of all atoms as if every occupancy is 1.00 is a frequent pitfall as it introduces - contacting or even overlapping atoms. -- Beware of composition: the atomic structure may not match the intended charge-neutral chemical - formula, especially hydrogen atoms, which may be missing or duplicated due to disorder. +If available, **computational** materials databases are the most recommended avenue +to obtain structures that are "computation-ready", or even better, already optimized +with some electronic-structure methods. On the other hand, structures that are from +**experimental** characterization are frequently not "computation-ready", and thus +should not be subject to optimization without careful validation in pre-processing. +This can be prominent for `cif` and `pdb` structures determined by powder or single- +crystal XRD which can be affected by sample quality and thermal motion. + +- Watch out for crystallographic disorder and atoms with low resolution or fractional + occupation: using the superposition of all atoms as if every occupancy is 1.00 is + highly likely to introduce contacting or even overlapping atoms. +- Beware of composition: the atomic structure may not match the intended macroscopic, + charge-neutral chemical formula, owing to missing or duplicated hydrogen atoms, + small counter ions, solvent or ligand molecules. + +Possible resolutions vary from simple manual editing in the modelling stage, to +utilization of supercells and enumeration of special quasirandom structures (common +for materials with dopants), and to more rigorous XRD refinement and application of +quantum crystallography methods. +``` + +```{note} +For variable-cell optimizations, the [CELL_OPT](#CP2K_INPUT.MOTION.CELL_OPT) section +clarifies that a cell representation with the first vector A along the X-axis and +the second vector B on the XY plane is expected. This is to ensure that the cell +vectors are updated correctly in response to the external pressure in the current +implementation; even when an unconventional vector definition could sometimes feel +easy (for instance the primitive rhombohedral cell of the face-centered cubic lattice), +an appropriate linear transformation to the cell and the atomic coordinates should +still be applied to create the model for input. For more details, see online +[github discussion](https://github.com/cp2k/cp2k/issues/3384). ``` ## Choosing an optimizer @@ -184,19 +244,17 @@ preliminary optimization with looser settings before tightening the convergence The most important atomic convergence criteria are [MAX_FORCE](#CP2K_INPUT.MOTION.GEO_OPT.MAX_FORCE), [RMS_FORCE](#CP2K_INPUT.MOTION.GEO_OPT.RMS_FORCE), [MAX_DR](#CP2K_INPUT.MOTION.GEO_OPT.MAX_DR), and -[RMS_DR](#CP2K_INPUT.MOTION.GEO_OPT.RMS_DR). Force criteria are given in Hartree/Bohr, while -displacement criteria are given in Bohr. For `CELL_OPT`, the pressure target and pressure -convergence are controlled by [EXTERNAL_PRESSURE](#CP2K_INPUT.MOTION.CELL_OPT.EXTERNAL_PRESSURE) and -[PRESSURE_TOLERANCE](#CP2K_INPUT.MOTION.CELL_OPT.PRESSURE_TOLERANCE). The optimization is considered -converged only when all active criteria are satisfied. +[RMS_DR](#CP2K_INPUT.MOTION.GEO_OPT.RMS_DR). Force criteria are given in unit of Hartree/Bohr, while +displacement criteria are given in Bohr. In addition, for `CELL_OPT`, the pressure target and +pressure convergence in bar are specified by +[EXTERNAL_PRESSURE](#CP2K_INPUT.MOTION.CELL_OPT.EXTERNAL_PRESSURE) and +[PRESSURE_TOLERANCE](#CP2K_INPUT.MOTION.CELL_OPT.PRESSURE_TOLERANCE) respectively. -```{note} -The convergence criteria of the LBFGS optimizer are not the same with that of the BFGS and CG -optimizers, so you may see a converged report without all criteria listed above meeting the -expectation. -``` - -A typical convergence report for a variable-cell optimization looks like: +A typical convergence report for a variable-cell optimization, as controlled by +[PROGRAM_RUN_INFO](#CP2K_INPUT.MOTION.CELL_OPT.PRINT.PROGRAM_RUN_INFO), looks like the following +where forces are reported as gradients. For a fixed-cell optimization, the +[PROGRAM_RUN_INFO](#CP2K_INPUT.MOTION.GEO_OPT.PRINT.PROGRAM_RUN_INFO) does not mention the pressure +but is otherwise the same. ```none OPT| ************************************************************************** @@ -232,7 +290,44 @@ A typical convergence report for a variable-cell optimization looks like: OPT| ************************************************************************** ``` -In the output, forces are reported as gradients. +Satisfying all of these criteria within `MAX_ITER` steps is a sufficient condition for the +convergence of optimization. With the `BFGS` and `CG` optimizers, this is also the necessary +condition. With the `LBFGS` optimizer, however, this is not a necessary condition because `LBFGS` +has its own extra pair of criteria named +[WANTED_PROJ_GRADIENT](#CP2K_INPUT.MOTION.GEO_OPT.LBFGS.WANTED_PROJ_GRADIENT) and +[WANTED_REL_F_ERROR](#CP2K_INPUT.MOTION.GEO_OPT.LBFGS.WANTED_REL_F_ERROR) which may override the +aforementioned general ones. It is possible to see an optimization task with the `LBFGS` optimizer +finish with the message below, even when one or more of the general criteria have not been met yet; +try restarting the optimization task with tightened criteria or alternative optimizers until +satisfactory convergence. + +```none + *********************************************** + * Specific L-BFGS convergence criteria + * WANTED_PROJ_GRADIENT and WANTED_REL_F_ERROR + * satisfied .... run CONVERGED! + *********************************************** +``` + +The `LBFGS` optimizer has its own [PRINT_LEVEL](#CP2K_INPUT.MOTION.GEO_OPT.LBFGS.PRINT_LEVEL) which +controls the verbosity of details printed to the main output. In addition, a separate `iterate.dat` +file is also created by `LBFGS` optimizer where the columns `projg` and `f` can be checked against +the `WANTED_PROJ_GRADIENT` and `WANTED_REL_F_ERROR` criteria. Omitting some headers, the main body +of an `iterate.dat` excerpt is reproduced below. + +```none + it nf nseg nact sub itls stepl tstep projg f + 0 1 - - - - - - 4.632D-02 -1.082D+00 + 1 2 1 0 --- 0 1.0D+00 6.6D-02 5.016D-02 -1.087D+00 + 2 4 1 0 --- 1 1.9D+00 1.3D-01 5.707D-02 -1.097D+00 + 3 6 1 0 --- 1 1.7D+00 1.3D-01 6.212D-02 -1.108D+00 + 4 8 1 0 --- 1 1.5D+00 1.3D-01 6.413D-02 -1.120D+00 + 5 10 1 0 --- 1 1.5D+00 1.3D-01 6.112D-02 -1.132D+00 + 6 11 3 2 con 0 1.0D+00 1.3D-01 4.979D-02 -1.143D+00 + 7 12 3 2 con 0 1.0D+00 1.3D-01 2.440D-02 -1.150D+00 + 8 14 1 0 con 1 5.1D-01 6.6D-02 3.685D-03 -1.152D+00 + 9 15 1 0 con 0 1.0D+00 1.2D-02 7.569D-04 -1.152D+00 +``` Geometry and cell optimization normally try to lower the energy, but convergence is determined by the active force, displacement, and pressure criteria rather than by the total energy alone. @@ -241,6 +336,15 @@ is usually acceptable if all active convergence criteria are satisfied. In some may indicate that the structure is not a real minimum, which can be checked by vibrational analysis where imaginary frequencies reveal unstable modes. +Generally speaking, the usual convergence criteria may be satisfied within dozens or a few hundreds +of iterations. The default `MAX_ITER` value of 200 is mostly sufficient and may be raised to about +300 ~ 400 to deal with difficult cases like massive flexible structures. Running a single-point +calculation with `RUN_TYPE ENERGY_FORCE` beforehand, and then multiplying the elapsed time by +expected iterations, provides an estimate for the total time cost of an optimization task. To +monitor a long optimization, do not focus on these numeric criteria alone; it is recommended to +download the geometry (trajectory) file and observe the evolution of the structure and cell in a +visualizer program occasionally, and if anything goes wrong, halt the program in time. + ## Constraints, cell degrees of freedom, and symmetry Atomic constraints are defined in [MOTION/CONSTRAINT](#CP2K_INPUT.MOTION.CONSTRAINT). A common @@ -264,7 +368,9 @@ contains atom indices in the order used in [COORD](#CP2K_INPUT.FORCE_EVAL.SUBSYS are useful, for example, for fixing the bottom layers of a slab, holding part of a support structure, or imposing a chemically intended restriction. Check them carefully before production calculations: accidental constraints are a common reason for apparently converged but physically -wrong structures. +wrong structures. For example, in a model of surface or micropore adsorption (physisorption or +chemisorption), the substrate should be allowed to relax locally rather than fully fixed to account +for the realistic interaction between the substrate (adsorbent) and the adsorbate. For `CELL_OPT`, also decide which cell degrees of freedom should be optimized. Useful keywords include [CONSTRAINT](#CP2K_INPUT.MOTION.CELL_OPT.CONSTRAINT) for fixing selected cell components, @@ -280,25 +386,58 @@ symmetry; otherwise they can prevent the structure from relaxing to a lower-symm periodic fixed-cell `GEO_OPT`, see [KEEP_SPACE_GROUP](#CP2K_INPUT.MOTION.GEO_OPT.KEEP_SPACE_GROUP) and related symmetry keywords. +```{note} +The `KEEP_SYMMETRY` keyword should always be used in conjunction with the cell symmetry +specified by [FORCE_EVAL/SUBSYS/CELL/SYMMETRY](#CP2K_INPUT.FORCE_EVAL.SUBSYS.CELL.SYMMETRY). +These keywords act only on the cell vectors and their gradients in a variable-cell +optimization, and they both assume a conventional cell where A is along the X-axis and +B is on the XY plane as noted in "Starting structure and cell" above. + +To allow for the use of `KEEP_SPACE_GROUP`, the `spglib` library should be installed +and linked to the CP2K build in order to detect and preserve the space group. Usually +`KEEP_SPACE_GROUP` is used together with `KEEP_SYMMETRY`. +``` + ## Electronic-structure settings +Density functional theory (DFT) is an electronic-structure (wavefunction) method routinely used for +optimization. This section elaborate on the relevant aspects, assuming basic knowledge about the +method which can be found at [](../dft/index.md). + ### SCF quality -For ab-initio geometry or cell optimization, every optimization step requires an -electronic-structure calculation. The force quality depends on the SCF convergence, grid settings, -basis set, pseudopotentials, and method-specific thresholds. For `CELL_OPT`, the stress tensor and -pressure must also be accurate enough for reliable cell updates. Important Quickstep settings -include [EPS_DEFAULT](#CP2K_INPUT.FORCE_EVAL.DFT.QS.EPS_DEFAULT), -[CUTOFF](#CP2K_INPUT.FORCE_EVAL.DFT.MGRID.CUTOFF), -[REL_CUTOFF](#CP2K_INPUT.FORCE_EVAL.DFT.MGRID.REL_CUTOFF), and -[EPS_SCF](#CP2K_INPUT.FORCE_EVAL.DFT.SCF.EPS_SCF). +Every optimization step requires an electronic-structure calculation, and the force quality depends +on the SCF convergence, grid settings, basis set, pseudopotentials, and method-specific thresholds. +For `CELL_OPT`, the stress tensor and pressure must also be accurate enough for reliable cell +updates. (The SCF convergence for each iteration is not to be confused with the convergence of +geometry in terms of stepsize, gradient and pressure across iterations in the optimization process.) Optimization is driven by forces, and `CELL_OPT` is additionally driven by stress, rather than by -total energies alone. Therefore, simply choosing `EPS_SCF`, `CUTOFF`, `REL_CUTOFF`, or `EPS_DEFAULT` -from a final static-energy test may not be reliable. These settings usually need to be at least as -strict as, and often stricter than, those used for routine single-point energy calculations, because -noisy or poorly converged forces and stresses can slow down the optimization, cause oscillations, -produce line-search failures, or lead to an unreliable structure or cell. +total energies alone. Therefore, simply adjusting settings and tuning parameters from a final +static-energy test may not be sufficient. They usually need to be at least as strict as, and often +stricter than, those used for routine single-point energy calculations, because noisy or poorly +converged forces and stresses can slow down the optimization, cause oscillations, produce +line-search failures, or lead to an unreliable structure or cell. + +Important Quickstep settings include [CUTOFF](#CP2K_INPUT.FORCE_EVAL.DFT.MGRID.CUTOFF), +[REL_CUTOFF](#CP2K_INPUT.FORCE_EVAL.DFT.MGRID.REL_CUTOFF), +[EPS_DEFAULT](#CP2K_INPUT.FORCE_EVAL.DFT.QS.EPS_DEFAULT), and +[EPS_SCF](#CP2K_INPUT.FORCE_EVAL.DFT.SCF.EPS_SCF). + +```{note} +The mesh spacing of the plane-wave and real-space grid for electronic-structure +calculations is determined both by the `CUTOFF` setting and by the dimensions of +the cell. A variable-cell optimization may witness significant variations of the +cell, which will also affect the number of grid points and introduce artificial +discontinuities to the energy, force, stress and other properties even if the +change in structure is small. + +It is recommended for better stability to set up a reference cell in the section +[FORCE_EVAL/SUBSYS/CELL/CELL_REF](#CP2K_INPUT.FORCE_EVAL.SUBSYS.CELL.CELL_REF), +with the dimension of the reference cell larger than the expected upper bound within +reach during the whole variable-cell optimization process, which will be kept constant +so that the number of grid points is held fixed and the change in energy more smooth. +``` ### Extrapolation @@ -328,31 +467,71 @@ this warning appears together with abnormal changes in the geometry, cell, total stress, or SCF behaviour, stop using that density-matrix-based method and switch to a safer alternative. +```{note} +The applicability of extrapolation schemes depends on versions. To exemplify, extrapolation +was formerly available for gamma-only calculations, but the latest versions of CP2K has +started supporting extrapolation with full k-point sampling. However, when k-points are +requested to be reduced by symmetry, the dynamic refreshing of symmetry of k-points alongside +geometry updates could make extrapolation unavailable and automatically fall back to `USE_GUESS`. +``` + +```{danger} +**Be responsible, and do not ignore SCF convergence failure blindly.** + +By default, a failure in SCF convergence aborts the program with a message like: +`SCF run NOT converged. To continue the calculation regardless, please set the keyword +IGNORE_CONVERGENCE_FAILURE.` Setting +[IGNORE_CONVERGENCE_FAILURE](#CP2K_INPUT.FORCE_EVAL.DFT.SCF.IGNORE_CONVERGENCE_FAILURE) +to `.TRUE.` turns it into a warning `SCF run NOT converged` that allows for optimization +to continue. However, bad SCF convergence leads to unreliable energy, force, and stress +on the current step, introducing error to the updated structure and the extrapolated +wavefunction on the next step. An optimization process with most or all steps failing +to converge SCF cycles does not yield trustworthy results in the end. + +Therefore, `IGNORE_CONVERGENCE_FAILURE` should only be considered as a **last resort** +out of desperation, rather than a universal remedy used on a regular basis or even as +the default. There are dozens of measures available towards SCF convergence on top of +a realistic, chemically sensible model structure and appropriate net charge, spin +multiplicity, atomic magnetization, and k-point sampling; please try achieving SCF +convergence on the starting structure in a single-point calculation in the first place. +``` + ## Output files and restarts -The main output file contains the optimization progress and convergence checks. CP2K also writes -geometry and restart files whose names are based on [PROJECT_NAME](#CP2K_INPUT.GLOBAL.PROJECT_NAME) -(or `PROJECT`): +The main output file contains the optimization progress and convergence checks. Geometry and restart +files are also written with names based on the global +[PROJECT_NAME](#CP2K_INPUT.GLOBAL.PROJECT_NAME) (or `PROJECT`) as well as relevant printkeys under +the [MOTION/PRINT](#CP2K_INPUT.MOTION.PRINT) section. -- `project-pos-1.xyz` contains the sequence of geometries visited during the optimization. The last - frame is the latest optimized geometry. -- Cell information and stress output, when requested, are controlled by print keys such as - [MOTION/PRINT/CELL](#CP2K_INPUT.MOTION.PRINT.CELL) and - [MOTION/PRINT/STRESS](#CP2K_INPUT.MOTION.PRINT.STRESS). -- `project-FINAL-1_{cycles}.cif` is a CIF file containing the cell information of the final - structure. -- `project-FINAL-1_{cycles}.xyz` is an extented XYZ file containing the cell information of final - structure +- `project-pos-1.xyz` contains the sequence of geometries visited during the optimization. This is + controlled by the [TRAJECTORY](#CP2K_INPUT.MOTION.PRINT.TRAJECTORY) section with the default + format of XYZ (or XMOL). Note that the original XYZ specification does not convey periodicity, and + not all of the visualizers can parse the additional cell information in the comment line of XYZ + (e.g. in extended XYZ specification); switching to other formats using the `FORMAT` keyword may be + needed for visualization. +- `project-FINAL-1_{iter}.cif` and `project-FINAL-1_{iter}.xyz` is a CIF and an extended XYZ file + for the final structure, where `{iter}` is the number of iterations (steps) during optimization + until reaching convergence or hitting `MAX_ITER`. These are controlled by + [FINAL_STRUCTURE](#CP2K_INPUT.MOTION.PRINT.FINAL_STRUCTURE). +- `project-1.cell` contains the cell vectors during optimization, as controlled by + [CELL](#CP2K_INPUT.MOTION.PRINT.CELL). +- `project-1.stress` contains the stress tensors during optimization, as controlled by + [STRESS](#CP2K_INPUT.MOTION.PRINT.STRESS). - `project-1.restart` is a CP2K restart input containing the latest geometry, cell, and relevant - settings. -- `project-1.restart.bak-*` are backup restart files from previous optimization steps. + settings, as controlled by [RESTART](#CP2K_INPUT.MOTION.PRINT.RESTART). When `BACKUP_COPIES` under + this section is non-zero, backup restart files from previous optimization steps will also be + preserved as `project-1.restart.bak-*`. -If a job stops before convergence, continue from the latest restart file, for example: +The restart file can be used directly as a new input file like ```none cp2k -i project-1.restart -o project-restart.out ``` +Alternatively, the [EXT_RESTART](#CP2K_INPUT.EXT_RESTART) section can be used for finer restart +controls. These will be convenient in case an optimization does not converge and a new optimization +starting from the latest structure is needed. + For expensive electronic-structure calculations, wavefunction restart files can reduce the cost of continuation or follow-up calculations; see [FORCE_EVAL/DFT/SCF/PRINT/RESTART](#CP2K_INPUT.FORCE_EVAL.DFT.SCF.PRINT.RESTART). For BFGS @@ -361,7 +540,7 @@ optimization driver. ## Practical checklist -Before trusting an optimized structure or cell, it's suggested to check that: +Before trusting an optimized structure or cell, it is suggested to check that: - All active force and displacement convergence criteria are satisfied; - For `CELL_OPT`, the pressure criterion is satisfied and the final cell is physically sensible; diff --git a/src/input_constants.F b/src/input_constants.F index b90bcda043..bbd8e758e6 100644 --- a/src/input_constants.F +++ b/src/input_constants.F @@ -555,6 +555,13 @@ MODULE input_constants do_second_rotation_step = 1, & do_third_rotation_step = 2 + ! See variable "iprint" in src/motion/cp_lbfgs.F for explanation + INTEGER, PARAMETER, PUBLIC :: silent_lbfgs = 0, & + low_lbfgs = 1, & + medium_lbfgs = 99, & + high_lbfgs = 100, & + debug_lbfgs = 101 + INTEGER, PARAMETER, PUBLIC :: xc_funct_no_shortcut = 0, & xc_funct_blyp = 1, & xc_funct_pade = 2, & diff --git a/src/motion/cp_lbfgs_geo.F b/src/motion/cp_lbfgs_geo.F index 7395cf2228..ef4abb4441 100644 --- a/src/motion/cp_lbfgs_geo.F +++ b/src/motion/cp_lbfgs_geo.F @@ -80,7 +80,6 @@ CONTAINS INTEGER :: handle, iter_nr, its, output_unit LOGICAL :: converged, should_stop - REAL(KIND=dp) :: trust_radius TYPE(cell_type), POINTER :: cell TYPE(cp_lbfgs_opt_gopt_type), POINTER :: optimizer TYPE(cp_logger_type), POINTER :: logger @@ -118,12 +117,14 @@ CONTAINS IF (gopt_env%type_id == default_ts_method_id) & CPABORT("BFGS method not yet working with DIMER") - CALL section_vals_val_get(geo_section, "LBFGS%TRUST_RADIUS", r_val=trust_radius) ALLOCATE (optimizer) CALL cp_opt_gopt_create(optimizer, para_env=para_env, obj_funct=gopt_env, & - x0=x0, wanted_relative_f_delta=gopt_param%wanted_rel_f_error, & - wanted_projected_gradient=gopt_param%wanted_proj_gradient, m=gopt_param%max_h_rank, & - max_f_per_iter=gopt_param%max_f_per_iter, trust_radius=trust_radius) + x0=x0, m=gopt_param%max_h_rank, & + print_every=gopt_param%print_control, & + wanted_relative_f_delta=gopt_param%wanted_rel_f_error, & + wanted_projected_gradient=gopt_param%wanted_proj_gradient, & + max_f_per_iter=gopt_param%max_f_per_iter, & + trust_radius=gopt_param%trust_radius) CALL cp_iterate(logger%iter_info, increment=0, iter_nr_out=iter_nr) converged = .FALSE. diff --git a/src/motion/gopt_param_types.F b/src/motion/gopt_param_types.F index 8170db6546..fee4071cd6 100644 --- a/src/motion/gopt_param_types.F +++ b/src/motion/gopt_param_types.F @@ -57,11 +57,11 @@ MODULE gopt_param_types INTEGER :: method_id = default_cell_method_id, type_id = default_cell_method_id INTEGER :: ts_method_id = 0, shellcore_method_id = 0 INTEGER :: max_f_per_iter = 0, max_iter = 0, max_h_rank = 0 - INTEGER :: max_steep_steps = 0 - REAL(KIND=dp) :: restart_limit = 0.0_dp - REAL(KIND=dp) :: wanted_proj_gradient = 0.0_dp, wanted_rel_f_error = 0.0_dp - REAL(KIND=dp) :: max_dr = 0.0_dp, max_force = 0.0_dp, rms_dr = 0.0_dp, rms_force = 0.0_dp - REAL(KIND=dp) :: dimer_angle_tol = 0.0_dp + INTEGER :: max_steep_steps = 0, print_control = 0 + REAL(KIND=dp) :: restart_limit = 0.0_dp + REAL(KIND=dp) :: trust_radius = -1.0_dp, wanted_proj_gradient = 0.0_dp, wanted_rel_f_error = 0.0_dp + REAL(KIND=dp) :: max_dr = 0.0_dp, max_force = 0.0_dp, rms_dr = 0.0_dp, rms_force = 0.0_dp + REAL(KIND=dp) :: dimer_angle_tol = 0.0_dp TYPE(cg_ls_param_type) :: cg_ls = cg_ls_param_type() END TYPE gopt_param_type @@ -81,6 +81,8 @@ CONTAINS TYPE(section_vals_type), POINTER :: gopt_section INTEGER, INTENT(IN), OPTIONAL :: type_id + LOGICAL :: explicit + CPASSERT(ASSOCIATED(gopt_section)) IF (PRESENT(type_id)) THEN @@ -100,8 +102,15 @@ CONTAINS CASE (default_lbfgs_method_id) CALL section_vals_val_get(gopt_section, "LBFGS%MAX_H_RANK", i_val=gopt_param%max_h_rank) CALL section_vals_val_get(gopt_section, "LBFGS%MAX_F_PER_ITER", i_val=gopt_param%max_f_per_iter) + CALL section_vals_val_get(gopt_section, "LBFGS%TRUST_RADIUS", r_val=gopt_param%trust_radius) CALL section_vals_val_get(gopt_section, "LBFGS%WANTED_PROJ_GRADIENT", r_val=gopt_param%wanted_proj_gradient) CALL section_vals_val_get(gopt_section, "LBFGS%WANTED_REL_F_ERROR", r_val=gopt_param%wanted_rel_f_error) + CALL section_vals_val_get(gopt_section, "LBFGS%__CONTROL_VAL", explicit=explicit) + IF (explicit) THEN + CALL section_vals_val_get(gopt_section, "LBFGS%__CONTROL_VAL", i_val=gopt_param%print_control) + ELSE + CALL section_vals_val_get(gopt_section, "LBFGS%PRINT_LEVEL", i_val=gopt_param%print_control) + END IF CASE (default_bfgs_method_id) ! Do nothing CASE (default_cg_method_id) diff --git a/src/start/input_cp2k_motion.F b/src/start/input_cp2k_motion.F index 1f016fc69b..818e45beaa 100644 --- a/src/start/input_cp2k_motion.F +++ b/src/start/input_cp2k_motion.F @@ -29,7 +29,7 @@ MODULE input_cp2k_motion medium_print_level USE cp_units, ONLY: cp_unit_to_cp2k USE input_constants, ONLY: & - default_bfgs_method_id, default_cg_method_id, default_dimer_method_id, & + debug_lbfgs, default_bfgs_method_id, default_cg_method_id, default_dimer_method_id, & default_lbfgs_method_id, default_minimization_method_id, default_ts_method_id, & do_mc_gemc_npt, do_mc_gemc_nvt, do_mc_traditional, do_mc_virial, fix_none, fix_x, fix_xy, & fix_xz, fix_y, fix_yz, fix_z, fmt_id_pdb, fmt_id_xyz, gaussian, helium_cell_shape_cube, & @@ -37,10 +37,10 @@ MODULE input_cp2k_motion helium_mdist_exponential, helium_mdist_gaussian, helium_mdist_linear, & helium_mdist_quadratic, helium_mdist_singlev, helium_mdist_uniform, & helium_sampling_ceperley, helium_sampling_worm, helium_solute_intpot_mwater, helium_solute_intpot_nnp, & - helium_solute_intpot_none, integrate_exact, integrate_numeric, ls_2pnt, ls_3pnt, ls_fit, & - ls_gold, ls_none, matrix_init_cholesky, matrix_init_diagonal, numerical, perm_cycle, & - perm_plain, propagator_pimd, propagator_rpmd, propagator_cmd, propagator_bcmd, transformation_normal, & - transformation_stage + helium_solute_intpot_none, high_lbfgs, integrate_exact, integrate_numeric, low_lbfgs, ls_2pnt, ls_3pnt, ls_fit, & + ls_gold, ls_none, matrix_init_cholesky, matrix_init_diagonal, medium_lbfgs, numerical, perm_cycle, & + perm_plain, propagator_pimd, propagator_rpmd, propagator_cmd, propagator_bcmd, silent_lbfgs, & + transformation_normal, transformation_stage USE input_cp2k_constraints, ONLY: create_constraint_section USE input_cp2k_free_energy, ONLY: create_fe_section USE input_cp2k_md, ONLY: create_md_section @@ -65,7 +65,8 @@ MODULE input_cp2k_motion real_t, & char_t USE kinds, ONLY: dp - USE string_utilities, ONLY: s2a + USE string_utilities, ONLY: newline, & + s2a #include "../base/base_uses.f90" IMPLICIT NONE @@ -1018,10 +1019,10 @@ CONTAINS "of the simulation cell. As is noted in FORCE_EVAL/SUBSYS/CELL, the "// & "program convention is that the first cell vector A lies along the "// & "X-axis and the second cell vector B is in the XY plane, such that "// & - "the cell vector matrix is a lower triangle. The algorithm support "// & - "for updating the three upper triangular components during a cell "// & - "optimization is not complete or tested, so the input structure has "// & - "to be prepared accordingly with these three components precisely 0 "// & + "the cell vector matrix is a lower triangle. There is no complete, "// & + "official algorithm support and/or tests for updating the three "// & + "upper triangular components during a cell optimization; please "// & + "prepare input accordingly with these three components precisely 0 "// & "even for cases like the primitive rhombohedral cell of the FCC lattice.", & just_optimizers=.TRUE., use_model_hessian=.FALSE.) @@ -1482,6 +1483,21 @@ CONTAINS CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) + CALL keyword_create(keyword, __LOCATION__, name="PRINT_LEVEL", & + description="How much output is written out by the LBFGS algorithm. ", & + usage="PRINT_LEVEL MEDIUM", & + enum_c_vals=s2a("SILENT", "LOW", "MEDIUM", "HIGH", "DEBUG"), & + enum_desc=s2a("Almost no output", & + "Little output about f and |proj g| every iteration", & + "Quite some output about details every iteration", & + "Lots of output about changes of active set and final x", & + "Everything is written out, useful for debugging purposes only"), & + enum_i_vals=[silent_lbfgs, low_lbfgs, medium_lbfgs, & + high_lbfgs, debug_lbfgs], & + default_i_val=low_print_level) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + CALL keyword_create(keyword, __LOCATION__, name="WANTED_PROJ_GRADIENT", & description="Convergence criterion (overrides the general ones):"// & " Requested norm threshold of the gradient multiplied"// & @@ -1508,6 +1524,14 @@ CONTAINS CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) + CALL keyword_create(keyword, __LOCATION__, name="__CONTROL_VAL", & + description="Hidden parameter that controls the printing behavior "// & + "of the LBFGS optimizer for advanced debug purposes. This option "// & + "overrides PRINT_LEVEL setting if explicit.", & + default_i_val=-1) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + END SUBROUTINE create_lbfgs_section ! **************************************************************************************************