Feature : a new implementation of (variable cell) relaxation (#1464)

* feature : implemented a new CG relaxation method

* feature : move some files around

* fix : include parallel_common header

* fix cmakelist.txt

* fix test cases (use old method)

* update cell_factor for vc calculations

* fix a bug in line search

* add constrained optimization

* set atom displacement to be 0 for constrained relaxation

* add fixed_ibrav (needs further testing!)

* fix : fixed_axes = shape

* add "fixed_atoms" keyword, and some documentation

* fix : include input.h header in read_atoms.cpp

* use global variable for fixed_atoms

* doc : if relax_new, then only cg is allowed

* move old method to a separated class

Co-authored-by: wenfei-li <liwenfei@gmail.com>

docs/advanced/input_files/input-main.md

@@ -20,7 +20,7 @@
   [method_sto](#method_sto) | [nbands_sto](#nbands_sto) | [nche_sto](#nche_sto) | [emin_sto](#emin_sto) | [emax_sto](#emax_sto) | [seed_sto](#seed_sto) | [initsto_freq](#initsto_freq) | [npart_sto](#npart_sto)
 - [Geometry relaxation](#geometry-relaxation)
 
-  [relax_nmax](#relax_nmax) | [relax_method](#relax_method) | [relax_cg_thr](#relax_cg_thr) | [relax_bfgs_w1](#relax_bfgs_w1) | [relax_bfgs_w2](#relax_bfgs_w2) | [relax_bfgs_rmax](#relax_bfgs_rmax) | [relax_bfgs_rmin](#relax_bfgs_rmin) | [relax_bfgs_init](#relax_bfgs_init) | [cal_force](#cal_force) | [force_thr](#force_thr) | [force_thr_ev](#force_thr_ev) | [force_thr_ev2](#force_thr_ev2) | [cal_stress](#cal_stress) | [stress_thr](#stress_thr) | [press1, press2, press3](#press1-press2-press3) | [fixed_axes](#fixed_axes) | [cell_factor](#cell_factor)
+  [relax_nmax](#relax_nmax) | [relax_method](#relax_method) | [relax_cg_thr](#relax_cg_thr) | [relax_bfgs_w1](#relax_bfgs_w1) | [relax_bfgs_w2](#relax_bfgs_w2) | [relax_bfgs_rmax](#relax_bfgs_rmax) | [relax_bfgs_rmin](#relax_bfgs_rmin) | [relax_bfgs_init](#relax_bfgs_init) | [cal_force](#cal_force) | [force_thr](#force_thr) | [force_thr_ev](#force_thr_ev) | [force_thr_ev2](#force_thr_ev2) | [cal_stress](#cal_stress) | [stress_thr](#stress_thr) | [press1, press2, press3](#press1-press2-press3) | [fixed_axes](#fixed_axes) | [cell_factor](#cell_factor) | [fixed_ibrav](#fixed_ibrav) | [relax_new](#relax_new)
 - [Variables related to output information](#variables-related-to-output-information)
 
   [out_force](#out_force) | [out_mul](#out_mul) | [out_freq_elec](#out_freq_elec) | [out_freq_ion](#out_freq_ion) | [out_chg](#out_chg) | [out_pot](#out_pot) | [out_dm](#out_dm) | [out_wfc_pw](#out_wfc_pw) | [out_wfc_r](#out_wfc_r) | [out_wfc_lcao](#out_wfc_lcao) | [out_dos](#out_dos) | [out_band](#out_band) | [out_proj_band](#out_proj_band) | [out_stru](#out_stru) | [out_level](#out_level) | [out_alllog](#out_alllog) | [out_mat_hs](#out_mat_hs) | [out_mat_r](#out_mat_r) | [out_mat_hs2](#out_mat_hs2) | [out_element_info](#out_element_info) | [restart_save](#restart_save) | [restart_load](#restart_load) | [dft_plus_dmft](#dft_plus_dmft) | [rpa](#rpa)
@@ -141,24 +141,24 @@ This part of variables are used to control general system parameters.
 ### latname
 
 - **Type**: String
-- **Description**: Specifies the type of Bravias lattice. When set to `test`, the three lattice vectors are supplied explicitly in STRU file. When set to certain Bravais lattice type, there is no need to provide lattice vector, but a few lattice parameters might be required. For more information regarding this parameter, consult the [page on STRU file](stru.md).
-  Available options are:
-  - `test`: free strcture.
-  - `sc`: simple cubie.
-  - `fcc`: face-centered cubic.
-  - `bcc`: body-centered cubic.
-  - `hexagonal`: hexagonal.
-  - `trigonal`: trigonal.
-  - `st`: simple tetragonal.
-  - `bct`: body-centered tetragonal.
-  - `so`: orthorhombic.
-  - `baco`: base-centered orthorhombic.
-  - `fco`: face-centered orthorhombic.
-  - `bco`: body-centered orthorhombic.
-  - `sm`: simple monoclinic.
-  - `bacm`: base-centered monoclinic.
-  - `triclinic`: triclinic.
-- **Default**: `test`
+- **Description**: Specifies the type of Bravias lattice. When set to `none`, the three lattice vectors are supplied explicitly in STRU file. When set to certain Bravais lattice type, there is no need to provide lattice vector, but a few lattice parameters might be required. For more information regarding this parameter, consult the [page on STRU file](stru.md).
+  Available options are (correspondence with ibrav in QE is given in parenthesis):
+  - `none`: free strcture.
+  - `sc`: simple cubic. (1)
+  - `fcc`: face-centered cubic. (2)
+  - `bcc`: body-centered cubic. (3)
+  - `hexagonal`: hexagonal. (4)
+  - `trigonal`: trigonal. (5)
+  - `st`: simple tetragonal. (6)
+  - `bct`: body-centered tetragonal. (7)
+  - `so`: orthorhombic. (8)
+  - `baco`: base-centered orthorhombic. (9)
+  - `fco`: face-centered orthorhombic. (10)
+  - `bco`: body-centered orthorhombic. (11)
+  - `sm`: simple monoclinic. (12)
+  - `bacm`: base-centered monoclinic. (13)
+  - `triclinic`: triclinic. (14)
+- **Default**: `none`
 
 ### init_wfc
 
@@ -733,19 +733,36 @@ This part of variables are used to control the geometry relaxation.
 - **Description**:which axes are fixed when do cell relaxation. Possible choices are:
   - None : default; all can relax
   - volume : relaxation with fixed volume
+  - shape : fix shape but change volume (i.e. only lattice constant changes)
   - a : fix a axis during relaxation
   - b : fix b axis during relaxation
   - c : fix c axis during relaxation
   - ab : fix both a and b axes during relaxation
   - ac : fix both a and c axes during relaxation
   - bc : fix both b and c axes during relaxation
-  - abc : fix all three axes during relaxation
+
+> Note : fixed_axes = "shape" and "volume" are only available for [relax_new](#relax_new) = 1
 - **Default**: None
 
+### fixed_ibrav
+
+- **Type**: Boolean
+- **Description**: when set to true, the lattice type will be preserved during relaxation. Must be used along with [relax_new](#relax_new) set to true, and a specific [latname](#latname) must be provided
+
+> Note: it is possible to use fixed_ibrav with fixed_axes, but please make sure you know what you are doing. For example, if we are doing relaxation of a simple cubic lattic (latname = "sc"), and we use fixed_ibrav along with fixed_axes = "volume", then the cell is never allowed to move and as a result the relaxation never converges.
+
+- **Default**: False
+
+### fixed_atoms
+
+- **Type**: Boolean
+- **Description**: when set to true, the direct coordinates of atoms will be preserved during variable-cell relaxation. If set to false, users can still fix certain components of certain atoms by using the `m` keyword in `STRU` file. For the latter option, check the end of this [instruction](stru.md).
+- **Default**: False
+
 ### relax_method
 
 - **Type**: String
-- **Description**: The method to do geometry optimizations:
+- **Description**: The method to do geometry optimizations, note that if relax_new is set to 1, then only cg is available:
   - bfgs: using BFGS algorithm.
   - sd: using steepest-descent algorithm.
   - cg: using cg algorithm.
@@ -758,6 +775,13 @@ This part of variables are used to control the geometry relaxation.
 - **Description**: When move-method is set to 'cg-bfgs', a mixed cg-bfgs algorithm is used. The ions first move according to cg method, then switched to bfgs when maximum of force on atoms is reduced below cg-threshold. Unit is eV/Angstrom.
 - **Default**: 0.5
 
+### relax_new
+
+- **Type**: Boolean
+- **Description**: At around end of 2022 we mad a new implemention of the CG method for relax and cell-relax calculations. But the old implementation was also kept. To use the new method, set relax_new to be true. To use the old one, set it to be false.
+
+- **Default**: True
+
 ### cell_factor
 
 - **Type**: Real

docs/advanced/opt.md

@@ -69,4 +69,12 @@ then the first Al atom will not be allowed to move in z direction.
 Fixing atomic position is sometimes helpful during relaxation of isolated molecule/cluster, to prevent the system from drifting in space.
 
 ### Fixing Cell Parameters
-Sometimes we want to do variable-cell relaxation with some of the cell degrees of freedom fixed. This is achieved by the [fixed_axes](./input_files/input-main.md#fixed_axes) keyword. We offer the options of fixing certain axis(axes), or fixing the volume.
+Sometimes we want to do variable-cell relaxation with some of the cell degrees of freedom fixed. This is achieved by keywords such as [fixed_axes](./input_files/input-main.md#fixed_axes), [fixed_ibrav](./input_files/input-main.md#fixed_ibrav) and [fixed_atoms](./input_files/input-main.md#fixed_atoms). Specifically, if users are familiar with the `ISIF` option from VASP, then we offer the following correspondence:
+
+- ISIF = 0 : calculation = "relax"
+- ISIF = 1, 2 : calculation = "relax", cal_stress = 1
+- ISIF = 3 : calculation = "cell-relax"
+- ISIF = 4 : calculation = "cell-relax", fixed_axes = "volume"
+- ISIF = 5 : calculation = "cell-relax", fixed_axes = "volume", fixed_atoms = True
+- ISIF = 6 : calculation = "cell-relax", fixed_atoms = True
+- ISIF = 7 : calculation = "cell-realx", fixed_axes = "shape", fixed_atoms = True

source/CMakeLists.txt

@@ -15,7 +15,7 @@ add_subdirectory(module_xc)
 add_subdirectory(module_esolver)
 add_subdirectory(module_gint)
 add_subdirectory(src_io)
-add_subdirectory(module_relaxation)
+add_subdirectory(module_relax)
 add_subdirectory(src_lcao)
 add_subdirectory(src_parallel)
 add_subdirectory(src_pdiag)

source/Makefile.Objects

@@ -36,7 +36,9 @@ VPATH=./src_global:\
 ./module_gint:\
 ./src_pw:\
 ./src_lcao:\
-./module_relaxation:\
+./module_relax:\
+./module_relax/relax_old:\
+./module_relax/relax_new:\
 ./module_rpa:\
 ./module_vdw:\
 ./src_io:\
@@ -280,7 +282,7 @@ OBJS_PW=fft.o\
     pw_transform_k.o\
 
 OBJS_RELAXATION=bfgs_basic.o\
-    ions.o\
+    relax_driver.o\
     ions_move_basic.o\
     ions_move_bfgs.o\
     ions_move_cg.o\
@@ -290,7 +292,9 @@ OBJS_RELAXATION=bfgs_basic.o\
     lattice_change_cg.o\
     lattice_change_methods.o\
     variable_cell.o\
-    relaxation.o\
+    relax_old.o\
+    relax.o\
+    line_search.o\
 
 OBJS_RPA=rpa.o\
 

source/driver_run.cpp

@@ -73,8 +73,8 @@ void Driver::driver_run()
     }
     else // scf; cell relaxation; nscf; etc
     {
-        Ions ions;
-        ions.opt_ions(p_esolver);
+        Relax_Driver rl_driver;
+        rl_driver.relax_driver(p_esolver);
     }
     //---------------------------MD/Relax------------------
 

source/input.cpp

@@ -130,7 +130,7 @@ void Input::Default(void)
     read_file_dir = "auto";
     // pseudo_type = "auto"; // mohan add 2013-05-20 (xiaohui add 2013-06-23)
     wannier_card = "";
-    latname = "test";
+    latname = "none";
     // xiaohui modify 2015-09-15, relax -> scf
     // calculation = "relax";
     calculation = "scf";
@@ -195,15 +195,19 @@ void Input::Default(void)
     press3 = 0.0;
     cal_stress = false;
     fixed_axes = "None"; // pengfei 2018-11-9
+    fixed_ibrav = false;
+    fixed_atoms = false;
     relax_method = "cg"; // pengfei  2014-10-13
     relax_cg_thr = 0.5; // pengfei add 2013-08-15
     out_level = "ie";
     out_md_control = false;
+    relax_new = true;
     relax_bfgs_w1 = 0.01; // mohan add 2011-03-13
     relax_bfgs_w2 = 0.5;
     relax_bfgs_rmax = 0.8; // bohr
     relax_bfgs_rmin = 1e-5;
     relax_bfgs_init = 0.5; // bohr
+    relax_scale_force = 0.5;
     nbspline = -1;
     //----------------------------------------------------------
     // ecutwfc
@@ -804,6 +808,14 @@ bool Input::Read(const std::string &fn)
         {
             read_value(ifs, fixed_axes);
         }
+        else if (strcmp("fixed_ibrav", word) == 0)
+        {
+            read_value(ifs, fixed_ibrav);
+        }
+        else if (strcmp("fixed_atoms", word) == 0)
+        {
+            read_value(ifs, fixed_atoms);
+        }
         else if (strcmp("relax_method", word) == 0)
         {
             read_value(ifs, relax_method);
@@ -837,12 +849,15 @@ bool Input::Read(const std::string &fn)
         {
             read_value(ifs, relax_bfgs_init);
         }
-        //		else if (strcmp("gauss_pao_flag", word) == 0)
-        //		else if (strcmp("gauss_pao_flag", word) == 0)
-        //		else if (strcmp("gauss_pao_flag", word) == 0)
-        //		{
-        //			read_value(ifs, gauss_PAO_flag);
-        //		}
+        else if (strcmp("relax_scale_force", word) == 0)
+        {
+            read_value(ifs, relax_scale_force);
+        }
+        else if (strcmp("relax_new", word) == 0)
+        {
+            read_value(ifs, relax_new);
+        }
+
         //----------------------------------------------------------
         // plane waves
         //----------------------------------------------------------
@@ -2168,6 +2183,8 @@ void Input::Bcast()
     Parallel_Common::bcast_double(press3);
     Parallel_Common::bcast_bool(cal_stress);
     Parallel_Common::bcast_string(fixed_axes);
+    Parallel_Common::bcast_bool(fixed_ibrav);
+    Parallel_Common::bcast_bool(fixed_atoms);
     Parallel_Common::bcast_string(relax_method);
     Parallel_Common::bcast_double(relax_cg_thr); // pengfei add 2013-08-15
     Parallel_Common::bcast_string(out_level);
@@ -2177,6 +2194,8 @@ void Input::Bcast()
     Parallel_Common::bcast_double(relax_bfgs_rmax);
     Parallel_Common::bcast_double(relax_bfgs_rmin);
     Parallel_Common::bcast_double(relax_bfgs_init);
+    Parallel_Common::bcast_double(relax_scale_force);
+    Parallel_Common::bcast_bool(relax_new);
 
     Parallel_Common::bcast_bool(gamma_only);
     Parallel_Common::bcast_bool(gamma_only_local);

source/input.h

@@ -120,8 +120,17 @@ class Input
     bool cal_stress; // calculate the stress
 
     std::string fixed_axes; // which axes are fixed
+    bool fixed_ibrav; //whether to keep type of lattice; must be used along with latname
+    bool fixed_atoms; //whether to fix atoms during vc-relax
     std::string relax_method; // methods to move_ion: sd, bfgs, cg...
 
+    //For now, this is only relevant if we choose to use
+    //CG relaxation method. If set to true, then the new
+    //implementation will be used; if set to false, then
+    //the original implementation will be used
+    //Default is true
+    bool relax_new;
+
     double relax_cg_thr; // threshold when cg to bfgs, pengfei add 2011-08-15
 
     double relax_bfgs_w1; // wolfe condition 1
@@ -131,6 +140,8 @@ class Input
     double relax_bfgs_rmin; // trust radius min
     double relax_bfgs_init; // initial move
 
+    double relax_scale_force;
+
     //==========================================================
     // Planewave
     //==========================================================

source/input_conv.cpp

@@ -11,7 +11,7 @@
 #include "src_io/epsilon0_pwscf.h"
 #include "src_io/epsilon0_vasp.h"
 #include "src_io/optical.h"
-#include "module_relaxation/ions_move_basic.h"
+#include "module_relax/relax_old/ions_move_basic.h"
 #include "src_pw/global.h"
 #include "src_pw/occupy.h"
 #ifdef __EXX
@@ -50,6 +50,23 @@ void Input_Conv::Convert(void)
         GlobalV::global_orbital_dir = INPUT.orbital_dir + "/";
     // GlobalV::global_pseudo_type = INPUT.pseudo_type;
     GlobalC::ucell.setup(INPUT.latname, INPUT.ntype, INPUT.lmaxmax, INPUT.init_vel, INPUT.fixed_axes);
+    if(INPUT.fixed_ibrav && !INPUT.relax_new)
+    {
+        ModuleBase::WARNING_QUIT("Input_Conv","fixed_ibrav only available for relax_new = 1");
+    }
+    if(INPUT.latname=="none" && INPUT.fixed_ibrav)
+    {
+        ModuleBase::WARNING_QUIT("Input_Conv","to use fixed_ibrav, latname must be provided");
+    }
+    if(INPUT.calculation == "relax" && INPUT.fixed_atoms)
+    {
+        ModuleBase::WARNING_QUIT("Input_Conv","fixed_atoms is not meant to be used for calculation = relax");
+    }
+    if(INPUT.relax_new && INPUT.relax_method!="cg")
+    {
+        ModuleBase::WARNING_QUIT("Input_Conv","only CG has been implemented for relax_new");
+    }
+    GlobalV::fixed_atoms = INPUT.fixed_atoms;
 
     GlobalV::KSPACING = INPUT.kspacing;
     GlobalV::MIN_DIST_COEF = INPUT.min_dist_coef;
@@ -109,6 +126,9 @@ void Input_Conv::Convert(void)
     GlobalV::CAL_STRESS = INPUT.cal_stress;
 
     GlobalV::RELAX_METHOD = INPUT.relax_method;
+    GlobalV::relax_scale_force = INPUT.relax_scale_force;
+    GlobalV::relax_new = INPUT.relax_new;
+
     GlobalV::OUT_LEVEL = INPUT.out_level;
     Ions_Move_CG::RELAX_CG_THR = INPUT.relax_cg_thr; // pengfei add 2013-09-09
 
@@ -354,6 +374,11 @@ void Input_Conv::Convert(void)
         }
     }
 
+    if(GlobalV::CALCULATION=="cell-relax" && INPUT.cell_factor < 2.0)
+    {
+        INPUT.cell_factor = 2.0; //follows QE
+    }
+
 //----------------------------------------------------------
 // about exx, Peize Lin add 2018-06-20
 //----------------------------------------------------------

source/input_update.cpp

@@ -7,7 +7,7 @@
 #include "module_base/global_function.h"
 #include "module_base/global_variable.h"
 #include "input.h"
-#include "module_relaxation/ions_move_basic.h"
+#include "module_relax/relax_old/ions_move_basic.h"
 #include "src_io/optical.h"
 #ifdef __LCAO
 #include "src_lcao/FORCE_STRESS.h"

source/module_base/global_variable.cpp

@@ -50,6 +50,9 @@ double PRESS3 = 0.0;
 double PRESSURE = 0.0;
 std::string RELAX_METHOD = "bfgs";
 std::string OUT_LEVEL = "ie";
+double relax_scale_force = 0.5;
+bool relax_new = true;
+bool fixed_atoms = false;
 int OUT_FREQ_ELEC = 0;
 int OUT_FREQ_ION = 0;
 int RELAX_NMAX = 20;

source/module_base/global_variable.h

@@ -54,6 +54,11 @@ extern std::string OUT_LEVEL;
 extern int OUT_FREQ_ELEC;
 extern int OUT_FREQ_ION;
 
+extern double relax_scale_force;
+extern bool relax_new;
+
+extern bool fixed_atoms;
+
 extern int RELAX_NMAX; // 8.3
 extern int SCF_NMAX; // 8.4
 extern int MD_NSTEP;

source/module_base/matrix.cpp

@@ -285,6 +285,12 @@ void matrix::zero_out(void)
 		c[i] = 0.0;
 }
 
+void matrix::fill_out(const double x)
+{
+	const int size = nr*nc;
+	for(int i = 0; i < size; i++)
+		c[i] = x;
+}
 
 matrix transpose(const matrix &m)
 {
@@ -337,11 +343,13 @@ void matrix::get_extreme_eigen_values(double &ev_lower, double &ev_upper)const
 }
 
 // Peize Lin add 2017-05-27
-void matrix::reshape( const double nr_new, const double nc_new )
+void matrix::reshape( const double nr_new, const double nc_new, const bool flag_zero )
 {
 	assert( nr*nc == nr_new*nc_new );
 	nr=nr_new;
 	nc=nc_new;
+
+	if(flag_zero) this-> zero_out();
 }
 
 double trace_on(const matrix &A, const matrix &B)