diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index e10171a1d0..0d9df7f76f 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -1530,27 +1530,45 @@ These variables are used to control the geometry relaxation. ### relax_method - **Type**: Vector of string -- **Description**: The methods to do geometry optimization. - the first element: - - cg: using the conjugate gradient (CG) algorithm. Note that there are two implementations of the conjugate gradient (CG) method, see [relax_new](#relax_new). - - bfgs : using the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm. - - lbfgs: using the Limited-memory Broyden–Fletcher–Goldfarb–Shanno (LBFGS) algorithm. - - cg_bfgs: using the CG method for the initial steps, and switching to BFGS method when the force convergence is smaller than [relax_cg_thr](#relax_cg_thr). - - sd: using the steepest descent (SD) algorithm. - - fire: the Fast Inertial Relaxation Engine method (FIRE), a kind of molecular-dynamics-based relaxation algorithm, is implemented in the molecular dynamics (MD) module. The algorithm can be used by setting [calculation](#calculation) to `md` and [md_type](#md_type) to `fire`. Also ionic velocities should be set in this case. See [fire](../md.md#fire) for more details. - - the second element: - when the first element is bfgs, if the second parameter is 1, it indicates the use of the new BFGS algorithm; if the second parameter is not 1, it indicates the use of the old BFGS algorithm. -- **Default**: cg 1 -- **Note**:In the 3.10-LTS version, the type of this parameter is std::string. It can be set to "cg","bfgs","cg_bfgs","bfgs_trad","lbfgs","sd","fire". +- **Description**: The methods to do geometry optimization. The available algorithms depend on the [relax_new](#relax_new) setting. + + **First element** (algorithm selection): + - `cg`: Conjugate gradient (CG) algorithm. Available for both `relax_new = True` (default, simultaneous optimization) and `relax_new = False` (nested optimization). See [relax_new](#relax_new) for implementation details. + - `bfgs`: Broyden–Fletcher–Goldfarb–Shanno (BFGS) quasi-Newton algorithm. **Only available when `relax_new = False`**. + - `lbfgs`: Limited-memory BFGS algorithm, suitable for large systems. **Only available when `relax_new = False`**. + - `cg_bfgs`: Mixed method starting with CG and switching to BFGS when force convergence reaches [relax_cg_thr](#relax_cg_thr). **Only available when `relax_new = False`**. + - `sd`: Steepest descent algorithm. **Only available when `relax_new = False`**. Not recommended for production use. + - `fire`: Fast Inertial Relaxation Engine method, a molecular-dynamics-based relaxation algorithm. Use by setting [calculation](#calculation) to `md` and [md_type](#md_type) to `fire`. Ionic velocities must be set in STRU file. See [fire](../md.md#fire) for details. + + **Second element** (BFGS variant, only when first element is `bfgs`): + - `1`: Traditional BFGS that updates the Hessian matrix B and then inverts it. + - `2` or omitted: Default BFGS that directly updates the inverse Hessian (recommended). + +- **Default**: `cg 1` +- **Note**: In the 3.10-LTS version, the type of this parameter is std::string. It can be set to "cg", "bfgs", "cg_bfgs", "bfgs_trad", "lbfgs", "sd", "fire". ### relax_new - **Type**: Boolean -- **Description**: At around the end of 2022 we made a new implementation of the Conjugate Gradient (CG) method for `relax` and `cell-relax` calculations. But the old implementation was also kept. - - True: use the new implementation of CG method for `relax` and `cell-relax` calculations. - - False: use the old implementation of CG method for `relax` and `cell-relax` calculations. +- **Description**: Controls which implementation of geometry relaxation to use. At the end of 2022, a new implementation of the Conjugate Gradient (CG) method was introduced for `relax` and `cell-relax` calculations, while the old implementation was kept for backward compatibility. + + - **True** (default): Use the new CG implementation with the following features: + - Simultaneous optimization of ionic positions and cell parameters (for `cell-relax`) + - Line search algorithm for step size determination + - Only CG algorithm is available (`relax_method` must be `cg`) + - Supports advanced cell constraints: `fixed_axes = "shape"`, `"volume"`, `"a"`, `"b"`, `"c"`, etc. + - Supports `fixed_ibrav` to maintain lattice type + - More efficient for variable-cell relaxation + - Step size controlled by [relax_scale_force](#relax_scale_force) + + - **False**: Use the old implementation with the following features: + - Nested optimization procedure: ionic positions optimized first, then cell parameters (for `cell-relax`) + - Multiple algorithms available: `cg`, `bfgs`, `lbfgs`, `sd`, `cg_bfgs` + - Limited cell constraints: only `fixed_axes = "volume"` is supported + - Traditional approach with separate ionic and cell optimization steps + - **Default**: True +- **Recommendation**: Use `relax_new = True` (default) for most cases, especially for `cell-relax` calculations. Use `relax_new = False` only if you need BFGS/LBFGS algorithms or for reproducing old results. ### relax_scale_force @@ -1568,7 +1586,8 @@ These variables are used to control the geometry relaxation. ### relax_cg_thr - **Type**: Real -- **Description**: When move-method is set to `cg_bfgs`, a mixed algorithm of conjugate gradient (CG) method and Broyden–Fletcher–Goldfarb–Shanno (BFGS) method is used. The ions first move according to CG method, then switched to BFGS method when the maximum of force on atoms is reduced below the CG force threshold, which is set by this parameter. +- **Availability**: Only used when `relax_new = False` and `relax_method = cg_bfgs` +- **Description**: When `relax_method` is set to `cg_bfgs`, a mixed algorithm of conjugate gradient (CG) and Broyden–Fletcher–Goldfarb–Shanno (BFGS) is used. The ions first move according to the CG method, then switch to the BFGS method when the maximum force on atoms is reduced below this threshold. - **Default**: 0.5 - **Unit**: eV/Angstrom @@ -1604,33 +1623,38 @@ These variables are used to control the geometry relaxation. ### relax_bfgs_w1 - **Type**: Real -- **Description**: Controls the Wolfe condition for Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm used in geometry relaxation. You can look into the paper Phys.Chem.Chem.Phys.,2000,2,2177 for more information. +- **Availability**: Only used when `relax_new = False` and `relax_method` is `bfgs` or `cg_bfgs` +- **Description**: Controls the Wolfe condition for the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm used in geometry relaxation. This parameter sets the sufficient decrease condition (c1 in Wolfe conditions). For more information, see Phys. Chem. Chem. Phys., 2000, 2, 2177. - **Default**: 0.01 ### relax_bfgs_w2 - **Type**: Real -- **Description**: Controls the Wolfe condition for Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm used in geometry relaxation. You can look into the paper Phys.Chem.Chem.Phys.,2000,2,2177 for more information. +- **Availability**: Only used when `relax_new = False` and `relax_method` is `bfgs` or `cg_bfgs` +- **Description**: Controls the Wolfe condition for the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm used in geometry relaxation. This parameter sets the curvature condition (c2 in Wolfe conditions). For more information, see Phys. Chem. Chem. Phys., 2000, 2, 2177. - **Default**: 0.5 ### relax_bfgs_rmax - **Type**: Real -- **Description**: For geometry optimization. It stands for the maximal movement of all the atoms. The sum of the movements from all atoms can be increased during the optimization steps. However, it can not be larger than `relax_bfgs_rmax`. -- **Unit**: Bohr +- **Availability**: Only used when `relax_new = False` and `relax_method` is `bfgs` or `cg_bfgs` +- **Description**: Maximum allowed total displacement of all atoms during geometry optimization. The sum of atomic displacements can increase during optimization steps but cannot exceed this value. - **Default**: 0.8 +- **Unit**: Bohr ### relax_bfgs_rmin - **Type**: Real -- **Description**: In old bfgs algorithm, it indicates the minimal movement of all the atoms. When the movement of all the atoms is smaller than relax_bfgs_rmin Bohr, and the force convergence is still not achieved, the calculation will break down. In the current default bfgs algorithm, this parameter is not used. +- **Availability**: Only used when `relax_new = False` and `relax_method = bfgs 1` (traditional BFGS) +- **Description**: Minimum allowed total displacement of all atoms. When the total atomic displacement falls below this value and force convergence is not achieved, the calculation will terminate. **Note**: This parameter is not used in the default BFGS algorithm (`relax_method = bfgs 2` or `bfgs`). - **Default**: 1e-5 - **Unit**: Bohr ### relax_bfgs_init - **Type**: Real -- **Description**: For geometry optimization. It stands for the sum of initial movements of all of the atoms. +- **Availability**: Only used when `relax_new = False` and `relax_method` is `bfgs` or `cg_bfgs` +- **Description**: Initial total displacement of all atoms in the first BFGS step. This sets the scale for the initial movement. - **Default**: 0.5 - **Unit**: Bohr @@ -1659,31 +1683,38 @@ These variables are used to control the geometry relaxation. ### fixed_axes - **Type**: String -- **Availability**: Only used when `calculation` set to `cell-relax` -- **Description**: Axes that are fixed during cell relaxation. Possible choices are: - - None**: default; all of the axes 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 - -> Note : fixed_axes = "shape" and "volume" are only available for [relax_new](#relax_new) = True +- **Availability**: Only used when `calculation` is set to `cell-relax` +- **Description**: Specifies which cell degrees of freedom are fixed during variable-cell relaxation. The available options depend on the [relax_new](#relax_new) setting: + + **When `relax_new = True` (default)**, all options are available: + - `None`: Default; all cell parameters can relax freely + - `volume`: Relaxation with fixed volume (allows shape changes) + - `shape`: Fix shape but allow volume changes (hydrostatic pressure only) + - `a`: Fix the a-axis lattice vector during relaxation + - `b`: Fix the b-axis lattice vector during relaxation + - `c`: Fix the c-axis lattice vector 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 + + **When `relax_new = False`**, all options are now available: + - `None`: Default; all cell parameters can relax freely + - `volume`: Relaxation with fixed volume (allows shape changes). Volume is preserved by rescaling the lattice after each update. + - `shape`: Fix shape but allow volume changes (hydrostatic pressure only). Stress tensor is replaced with isotropic pressure. + - `a`, `b`, `c`, `ab`, `ac`, `bc`: Fix specific lattice vectors. Gradients for fixed vectors are set to zero. - **Default**: None +- **Note**: For VASP users, see the [ISIF correspondence table](../opt.md#fixing-cell-parameters) in the geometry optimization documentation. Both implementations now support all constraint types. ### fixed_ibrav - **Type**: Boolean -- **Availability**: Must be used along with [relax_new](#relax_new) set to True, and a specific [latname](#latname) must be provided +- **Availability**: Can be used with both `relax_new = True` and `relax_new = False`. A specific [latname](#latname) must be provided. - **Description**: - - True: the lattice type will be preserved during relaxation + - True: the lattice type will be preserved during relaxation. The lattice vectors are reconstructed to match the specified Bravais lattice type after each update. - False: No restrictions are exerted during relaxation in terms of lattice type -> 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 lattice (`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. +> 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 lattice (`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. When both are used, `fixed_ibrav` is applied first, then `fixed_axes = "volume"` rescaling is applied. - **Default**: False diff --git a/docs/advanced/opt.md b/docs/advanced/opt.md index 9c2023b655..4f0bdf32a0 100644 --- a/docs/advanced/opt.md +++ b/docs/advanced/opt.md @@ -2,41 +2,76 @@ By setting `calculation` to be `relax` or `cell-relax`, ABACUS supports structural relaxation and variable-cell relaxation. -Current implementation of variable-cell relaxation in ABACUS now follows a nested procedure: fixed cell structural relaxation will be performed, followed by an update of the cell parameters, and the process is repeated until convergence is achieved. +ABACUS provides two implementations for variable-cell relaxation, controlled by the [relax_new](./input_files/input-main.md#relax_new) parameter: -An example of the variable cell relaxation can be found in our [repository](https://github.com/deepmodeling/abacus-develop/tree/develop/examples/relax/pw_al), which is provided with the reference output file log.ref. Note that in log.ref, each ionic step is labelled in the following manner: +- **New implementation** (`relax_new = True`, default since v3.8): Uses a simultaneous conjugate gradient (CG) optimization for both ionic positions and cell parameters. Both degrees of freedom are optimized together in each step. + +- **Old implementation** (`relax_new = False`): Follows a nested procedure where fixed-cell structural relaxation is performed first, followed by an update of the cell parameters, and the process is repeated until convergence is achieved. + +An example of the variable cell relaxation can be found in our [repository](https://github.com/deepmodeling/abacus-develop/tree/develop/examples/relax/pw_al), which is provided with the reference output file log.ref. When using the old implementation (`relax_new = False`), each ionic step is labelled in the following manner: ``` ------------------------------------------- RELAX CELL : 3 RELAX IONS : 1 (in total: 15) ------------------------------------------- -``` +``` indicating that this is the first ionic step of the 3rd cell configuration, and it is the 15-th ionic step in total. ## Optimization Algorithms -In the nested procedure mentioned above, we used CG method to perform cell relaxation, while offering four different algorithms for doing fixed-cell structural relaxation: BFGS, SD(steepest descent), CG(conjugate gradient), as well as a mixed CG-BFGS method. The optimziation algorithm can be selected using keyword [relax_method](./input_files/input-main.md#relax_method). We also provide a [list of keywords](./input_files/input-main.md#geometry-relaxation) for controlling the relaxation process. + +ABACUS offers multiple optimization algorithms for structural relaxation, which can be selected using the [relax_method](./input_files/input-main.md#relax_method) keyword. The available algorithms and their behavior depend on the [relax_new](./input_files/input-main.md#relax_new) setting: + +### Algorithm Availability + +**New implementation** (`relax_new = True`, default): +- **CG (Conjugate Gradient)**: Simultaneous optimization of both ionic positions and cell parameters using CG with line search. This is the only algorithm available for the new implementation. + +**Old implementation** (`relax_new = False`): +- **CG (Conjugate Gradient)**: For ionic relaxation; CG is also used for cell parameter optimization in the nested procedure +- **BFGS**: Quasi-Newton method for ionic relaxation +- **LBFGS**: Limited-memory BFGS for ionic relaxation +- **SD (Steepest Descent)**: Simple gradient descent for ionic relaxation +- **CG-BFGS**: Mixed method that starts with CG and switches to BFGS when force convergence reaches the threshold set by [relax_cg_thr](./input_files/input-main.md#relax_cg_thr) + +We also provide a [list of keywords](./input_files/input-main.md#geometry-relaxation) for controlling the relaxation process. ### BFGS method -The [BFGS method](https://en.wikipedia.org/wiki/Broyden%E2%80%93Fletcher%E2%80%93Goldfarb%E2%80%93Shanno_algorithm) is a quasi-Newton method for solving nonlinear optimization problem. It belongs to the class of quasi-Newton method where the Hessian matrix is approximated during the optimization process. If the initial point is not far from the extrema, BFGS tends to work better than gradient-based methods. +The [BFGS method](https://en.wikipedia.org/wiki/Broyden%E2%80%93Fletcher%E2%80%93Goldfarb%E2%80%93Shanno_algorithm) is a quasi-Newton method for solving nonlinear optimization problems. It belongs to the class of quasi-Newton methods where the Hessian matrix is approximated during the optimization process. If the initial point is not far from the extrema, BFGS tends to work better than gradient-based methods. + +**Note**: BFGS is only available with the old implementation (`relax_new = False`). + +ABACUS provides two BFGS implementations, controlled by the second element of [relax_method](./input_files/input-main.md#relax_method): + +- **Default BFGS** (`relax_method = bfgs 2` or `relax_method = bfgs`): Updates the inverse of the approximate Hessian matrix B directly. This is the recommended implementation. -There is an alternative traditional BFGS method, which can be called by using the keyword 'bfgs_trad'. The bfgs_trad method is a quasi-Newton method that substitute an approximate matrix B for the Hessian matrix. The main difference between 'bfgs' and 'bfgs_trad' is that 'bfgs' updates the inverse of matrix B while 'bfgs_trad' updates matrix B and obtains the inverse of B by solving the matrix eigenvalues and taking the reciprocal of the eigenvalues. Both methods are mathematically equivalent, but in some cases, 'bfgs_trad' performs better. +- **Traditional BFGS** (`relax_method = bfgs 1`): Updates the approximate Hessian matrix B itself, then obtains the inverse by solving matrix eigenvalues and taking their reciprocals. Both methods are mathematically equivalent, but in some cases the traditional variant may perform better. -In ABACUS, we implemented the BFGS method for doing fixed-cell structural relaxation. Users can choose which implementation of BFGS to call by adding the 'bfgs_trad' or 'bfgs' parameter. +### LBFGS method + +The [L-BFGS (Limited-memory BFGS)](https://en.wikipedia.org/wiki/Limited-memory_BFGS) method is a memory-efficient variant of BFGS that stores only a few vectors representing the Hessian approximation instead of the full matrix. This makes it particularly suitable for large systems with many atoms. + +**Note**: LBFGS is only available with the old implementation (`relax_new = False`). Set `relax_method = lbfgs` to use this method. ### SD method -The [SD (steepest descent) method](https://en.wikipedia.org/wiki/Gradient_descent) is one of the simplest first-order optimization methods, where in each step the motion is along the direction of the gradient, where the function descents the fastest. +The [SD (steepest descent) method](https://en.wikipedia.org/wiki/Gradient_descent) is one of the simplest first-order optimization methods, where in each step the motion is along the direction of the gradient, where the function descends the fastest. + +**Note**: SD is only available with the old implementation (`relax_new = False`). -In practice, SD method may take many iterations to converge, and is generally not used. +In practice, the SD method may take many iterations to converge, and is generally not recommended for production calculations. ### CG method The [CG (conjugate gradient) method](https://en.wikipedia.org/wiki/Conjugate_gradient_method) is one of the most widely used methods for solving optimization problems. -In ABACUS, we implemented the CG method for doing fixed-cell structural relaxation as well as the optimization of cell parameters. +ABACUS provides two implementations of the CG method: + +- **New CG implementation** (`relax_new = True`, default): Performs simultaneous optimization of both ionic positions and cell parameters using a line search algorithm. This implementation is more efficient for `cell-relax` calculations as it optimizes all degrees of freedom together. The step size can be controlled by [relax_scale_force](./input_files/input-main.md#relax_scale_force). + +- **Old CG implementation** (`relax_new = False`): Uses a nested procedure where ionic positions are optimized first using CG, followed by cell parameter optimization (also using CG) in `cell-relax` calculations. This is the traditional approach where the two optimization steps are separated. ## Constrained Optimization @@ -71,7 +106,26 @@ 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 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: + +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). + +**Available constraints by implementation:** + +- **New implementation** (`relax_new = True`): + - `fixed_axes = "shape"`: Only allows volume changes (hydrostatic pressure), cell shape is fixed + - `fixed_axes = "volume"`: Allows shape changes but keeps volume constant + - `fixed_axes = "a"`, `"b"`, `"c"`, etc.: Fix specific lattice vectors or combinations + - `fixed_ibrav = True`: Maintain the Bravais lattice type during relaxation + +- **Old implementation** (`relax_new = False`): + - **All `fixed_axes` options now supported**: "shape", "volume", "a", "b", "c", "ab", "ac", "bc", "abc" + - **`fixed_ibrav` now supported**: Maintains Bravais lattice type during relaxation + - Can combine `fixed_axes` with `fixed_ibrav` for constrained relaxation + - **Implementation approach**: Uses post-update constraint enforcement (volume rescaling and lattice reconstruction after each CG step) + +**VASP ISIF correspondence:** + +If you are familiar with the `ISIF` option from VASP, here is the correspondence: - ISIF = 0 : calculation = "relax" - ISIF = 1, 2 : calculation = "relax", cal_stress = 1 @@ -79,7 +133,7 @@ Sometimes we want to do variable-cell relaxation with some of the cell degrees o - 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 +- ISIF = 7 : calculation = "cell-relax", fixed_axes = "shape", fixed_atoms = True ### Stop Geometry Optimization Manually diff --git a/source/source_cell/unitcell.cpp b/source/source_cell/unitcell.cpp index e79a4206e8..bae4b1b2dd 100644 --- a/source/source_cell/unitcell.cpp +++ b/source/source_cell/unitcell.cpp @@ -428,18 +428,7 @@ void UnitCell::setup(const std::string& latname_in, this->lc[0] = 1; this->lc[1] = 1; this->lc[2] = 1; - if (!PARAM.inp.relax_new) { - ModuleBase::WARNING_QUIT( - "Input", - "there are bugs in the old implementation; set relax_new to be " - "1 for fixed_volume relaxation"); - } } else if (fixed_axes_in == "shape") { - if (!PARAM.inp.relax_new) { - ModuleBase::WARNING_QUIT( - "Input", - "set relax_new to be 1 for fixed_shape relaxation"); - } this->lc[0] = 1; this->lc[1] = 1; this->lc[2] = 1; diff --git a/source/source_relax/lattice_change_basic.cpp b/source/source_relax/lattice_change_basic.cpp index d6cfaec9ef..47b0c283e1 100644 --- a/source/source_relax/lattice_change_basic.cpp +++ b/source/source_relax/lattice_change_basic.cpp @@ -4,6 +4,7 @@ #include "source_base/global_variable.h" #include "source_base/parallel_common.h" #include "source_io/module_parameter/parameter.h" +#include "source_cell/update_cell.h" int Lattice_Change_Basic::dim = 0; bool Lattice_Change_Basic::converged = true; @@ -24,13 +25,27 @@ void Lattice_Change_Basic::setup_gradient(const UnitCell &ucell, double *lat, do { ModuleBase::TITLE("Lattice_Change_Basic", "setup_gradient"); - if (Lattice_Change_Basic::fixed_axes == "volume") + // Apply fixed_axes constraints to stress tensor + if (Lattice_Change_Basic::fixed_axes == "shape") + { + // Shape constraint: only volume can change (isotropic expansion/contraction) + // Replace stress with pure hydrostatic pressure + double pressure = (stress(0, 0) + stress(1, 1) + stress(2, 2)) / 3.0; + stress.zero_out(); + stress(0, 0) = pressure; + stress(1, 1) = pressure; + stress(2, 2) = pressure; + } + else if (Lattice_Change_Basic::fixed_axes == "volume") { + // Volume constraint: only shape can change + // Remove hydrostatic pressure component from stress double stress_aver = (stress(0, 0) + stress(1, 1) + stress(2, 2)) / 3.0; stress(0, 0) = stress(0, 0) - stress_aver; stress(1, 1) = stress(1, 1) - stress_aver; stress(2, 2) = stress(2, 2) - stress_aver; } + // Note: Axis constraints ("a", "b", "c", etc.) are handled via ucell.lc[] flags below lat[0] = ucell.latvec.e11 * ucell.lat0; lat[1] = ucell.latvec.e12 * ucell.lat0; @@ -42,24 +57,48 @@ void Lattice_Change_Basic::setup_gradient(const UnitCell &ucell, double *lat, do lat[7] = ucell.latvec.e32 * ucell.lat0; lat[8] = ucell.latvec.e33 * ucell.lat0; + // Calculate gradients for each lattice vector, or zero them if fixed if (ucell.lc[0] == 1) { grad[0] = -(lat[0] * stress(0, 0) + lat[1] * stress(1, 0) + lat[2] * stress(2, 0)); grad[1] = -(lat[0] * stress(0, 1) + lat[1] * stress(1, 1) + lat[2] * stress(2, 1)); grad[2] = -(lat[0] * stress(0, 2) + lat[1] * stress(1, 2) + lat[2] * stress(2, 2)); } + else + { + // Zero out gradient for fixed lattice vector a + grad[0] = 0.0; + grad[1] = 0.0; + grad[2] = 0.0; + } + if (ucell.lc[1] == 1) { grad[3] = -(lat[3] * stress(0, 0) + lat[4] * stress(1, 0) + lat[5] * stress(2, 0)); grad[4] = -(lat[3] * stress(0, 1) + lat[4] * stress(1, 1) + lat[5] * stress(2, 1)); grad[5] = -(lat[3] * stress(0, 2) + lat[4] * stress(1, 2) + lat[5] * stress(2, 2)); } + else + { + // Zero out gradient for fixed lattice vector b + grad[3] = 0.0; + grad[4] = 0.0; + grad[5] = 0.0; + } + if (ucell.lc[2] == 1) { grad[6] = -(lat[6] * stress(0, 0) + lat[7] * stress(1, 0) + lat[8] * stress(2, 0)); grad[7] = -(lat[6] * stress(0, 1) + lat[7] * stress(1, 1) + lat[8] * stress(2, 1)); grad[8] = -(lat[6] * stress(0, 2) + lat[7] * stress(1, 2) + lat[8] * stress(2, 2)); } + else + { + // Zero out gradient for fixed lattice vector c + grad[6] = 0.0; + grad[7] = 0.0; + grad[8] = 0.0; + } // grad[0] = -stress(0,0); grad[1] = -stress(0,1); grad[2] = -stress(0,2); // grad[3] = -stress(1,0); grad[4] = -stress(1,1); grad[5] = -stress(1,2); @@ -117,17 +156,57 @@ void Lattice_Change_Basic::change_lattice(UnitCell &ucell, double *move, double ucell.latvec.e33 = (move[8] + lat[8]) / ucell.lat0; } - ucell.a1.x = ucell.latvec.e11; - ucell.a1.y = ucell.latvec.e12; - ucell.a1.z = ucell.latvec.e13; - ucell.a2.x = ucell.latvec.e21; - ucell.a2.y = ucell.latvec.e22; - ucell.a2.z = ucell.latvec.e23; - ucell.a3.x = ucell.latvec.e31; - ucell.a3.y = ucell.latvec.e32; - ucell.a3.z = ucell.latvec.e33; - - ucell.omega = std::abs(ucell.latvec.Det()) * ucell.lat0 * ucell.lat0 * ucell.lat0; + // Apply post-update constraints + // Order matters: fixed_ibrav first, then volume rescaling + + // 1. Enforce Bravais lattice symmetry if fixed_ibrav is set + if (PARAM.inp.fixed_ibrav) + { + unitcell::remake_cell(ucell.lat); + } + + // 2. Enforce volume constraint by rescaling lattice + if (Lattice_Change_Basic::fixed_axes == "volume") + { + double omega_old = ucell.omega; // Volume before this update + double omega_new = std::abs(ucell.latvec.Det()) * ucell.lat0 * ucell.lat0 * ucell.lat0; + + // Rescale lattice to maintain constant volume + if (omega_new > 1e-10) // Avoid division by zero + { + double scale = std::pow(omega_old / omega_new, 1.0 / 3.0); + ucell.latvec *= scale; + + // Update a1, a2, a3 vectors + ucell.a1.x = ucell.latvec.e11; + ucell.a1.y = ucell.latvec.e12; + ucell.a1.z = ucell.latvec.e13; + ucell.a2.x = ucell.latvec.e21; + ucell.a2.y = ucell.latvec.e22; + ucell.a2.z = ucell.latvec.e23; + ucell.a3.x = ucell.latvec.e31; + ucell.a3.y = ucell.latvec.e32; + ucell.a3.z = ucell.latvec.e33; + + // Recalculate omega (should be equal to omega_old now) + ucell.omega = std::abs(ucell.latvec.Det()) * ucell.lat0 * ucell.lat0 * ucell.lat0; + } + } + else + { + // Update a1, a2, a3 vectors for non-volume-constrained cases + ucell.a1.x = ucell.latvec.e11; + ucell.a1.y = ucell.latvec.e12; + ucell.a1.z = ucell.latvec.e13; + ucell.a2.x = ucell.latvec.e21; + ucell.a2.y = ucell.latvec.e22; + ucell.a2.z = ucell.latvec.e23; + ucell.a3.x = ucell.latvec.e31; + ucell.a3.y = ucell.latvec.e32; + ucell.a3.z = ucell.latvec.e33; + + ucell.omega = std::abs(ucell.latvec.Det()) * ucell.lat0 * ucell.lat0 * ucell.lat0; + } ucell.GT = ucell.latvec.Inverse(); ucell.G = ucell.GT.Transpose(); diff --git a/source/source_relax/test/CMakeLists.txt b/source/source_relax/test/CMakeLists.txt index 5f8a600cc6..f55360dd00 100644 --- a/source/source_relax/test/CMakeLists.txt +++ b/source/source_relax/test/CMakeLists.txt @@ -31,22 +31,23 @@ list(APPEND cell_source_files ) AddTest( TARGET lattice_change_methods_test - LIBS parameter ${math_libs} base device - SOURCES lattice_change_methods_test.cpp ../lattice_change_methods.cpp ../lattice_change_basic.cpp + LIBS parameter ${math_libs} base device + SOURCES lattice_change_methods_test.cpp ../lattice_change_methods.cpp ../lattice_change_basic.cpp mock_remake_cell.cpp ) AddTest( TARGET lattice_change_basic_test - LIBS parameter ${math_libs} base device - SOURCES lattice_change_basic_test.cpp ../lattice_change_basic.cpp + LIBS parameter ${math_libs} base device + SOURCES lattice_change_basic_test.cpp ../lattice_change_basic.cpp mock_remake_cell.cpp ) AddTest( TARGET lattice_change_cg_test - LIBS parameter ${math_libs} base device - SOURCES lattice_change_cg_test.cpp - ../lattice_change_cg.cpp + LIBS parameter ${math_libs} base device + SOURCES lattice_change_cg_test.cpp + ../lattice_change_cg.cpp ../lattice_change_basic.cpp + mock_remake_cell.cpp ../../source_io/orb_io.cpp ) diff --git a/source/source_relax/test/ions_move_methods_test.cpp b/source/source_relax/test/ions_move_methods_test.cpp index f0027be236..464377f643 100644 --- a/source/source_relax/test/ions_move_methods_test.cpp +++ b/source/source_relax/test/ions_move_methods_test.cpp @@ -19,6 +19,64 @@ * - Ions_Move_Methods::get_update_iter() */ + // Mock the remake_cell function from update_cell.h +namespace unitcell +{ + // Track if remake_cell was called and with what lattice + static bool remake_cell_called = false; + static std::string remake_cell_latName; + static ModuleBase::Matrix3 remake_cell_latvec; + + void remake_cell(Lattice& lat) + { + remake_cell_called = true; + remake_cell_latName = lat.latName; + remake_cell_latvec = lat.latvec; + + // Mock implementation: enforce simple cubic structure for "sc" + if (lat.latName == "sc") + { + double celldm = std::sqrt(lat.latvec.e11 * lat.latvec.e11 + + lat.latvec.e12 * lat.latvec.e12 + + lat.latvec.e13 * lat.latvec.e13); + lat.latvec.Zero(); + lat.latvec.e11 = celldm; + lat.latvec.e22 = celldm; + lat.latvec.e33 = celldm; + } + // Mock implementation: enforce FCC structure for "fcc" + else if (lat.latName == "fcc") + { + double celldm = std::sqrt(lat.latvec.e11 * lat.latvec.e11 + + lat.latvec.e12 * lat.latvec.e12 + + lat.latvec.e13 * lat.latvec.e13) / std::sqrt(2.0); + lat.latvec.e11 = -celldm; + lat.latvec.e12 = 0.0; + lat.latvec.e13 = celldm; + lat.latvec.e21 = 0.0; + lat.latvec.e22 = celldm; + lat.latvec.e23 = celldm; + lat.latvec.e31 = -celldm; + lat.latvec.e32 = celldm; + lat.latvec.e33 = 0.0; + } + } + + // Helper function to reset mock state + void reset_remake_cell_mock() + { + remake_cell_called = false; + remake_cell_latName = ""; + remake_cell_latvec.Zero(); + } + + // Helper function to check if remake_cell was called + bool was_remake_cell_called() + { + return remake_cell_called; + } +} + // Define a fixture for the tests class IonsMoveMethodsTest : public ::testing::Test { diff --git a/source/source_relax/test/lattice_change_basic_test.cpp b/source/source_relax/test/lattice_change_basic_test.cpp index 4a12fca3fb..b40af6d88a 100644 --- a/source/source_relax/test/lattice_change_basic_test.cpp +++ b/source/source_relax/test/lattice_change_basic_test.cpp @@ -1,10 +1,13 @@ #include "source_relax/lattice_change_basic.h" +#include "mock_remake_cell.h" #include "for_test.h" #include "gtest/gtest.h" +#include "gmock/gmock.h" #define private public #include "source_io/module_parameter/parameter.h" #undef private + /************************************************ * unit tests of namespace Lattice_Change_Basic ***********************************************/ @@ -30,11 +33,17 @@ class LatticeChangeBasicTest : public ::testing::Test { // Initialize variables before each test stress.create(3, 3); + // Reset mock state before each test + unitcell::reset_remake_cell_mock(); + // Reset fixed_ibrav to default + PARAM.input.fixed_ibrav = false; } virtual void TearDown() { // Clean up after each test + unitcell::reset_remake_cell_mock(); + PARAM.input.fixed_ibrav = false; } }; @@ -525,3 +534,584 @@ TEST_F(LatticeChangeBasicTest, SetupEtotJudgementFalse) EXPECT_DOUBLE_EQ(80.0, Lattice_Change_Basic::etot); EXPECT_DOUBLE_EQ(-20.0, Lattice_Change_Basic::ediff); } + +// ============================================================================ +// NEW TESTS FOR SHAPE CONSTRAINT, VOLUME RESCALING, AND FIXED_IBRAV +// ============================================================================ + +// Test the setup_gradient function with fixed_axes = "shape" +TEST_F(LatticeChangeBasicTest, SetupGradientShape) +{ + // Initialize variables + ucell.lc[0] = 1; + ucell.lc[1] = 1; + ucell.lc[2] = 1; + + // Non-isotropic stress tensor + stress(0, 0) = 1.0; + stress(0, 1) = 2.0; + stress(0, 2) = 3.0; + stress(1, 0) = 4.0; + stress(1, 1) = 5.0; + stress(1, 2) = 6.0; + stress(2, 0) = 7.0; + stress(2, 1) = 8.0; + stress(2, 2) = 9.0; + + Lattice_Change_Basic::fixed_axes = "shape"; + + // Call setup_gradient method + Lattice_Change_Basic::setup_gradient(ucell, lat, grad, stress); + + // Check that stress becomes isotropic (only diagonal, all equal) + // Average pressure = (1 + 5 + 9) / 3 = 5 + EXPECT_DOUBLE_EQ(stress(0, 0), 5.0); + EXPECT_DOUBLE_EQ(stress(1, 1), 5.0); + EXPECT_DOUBLE_EQ(stress(2, 2), 5.0); + + // Off-diagonal elements should be zero + EXPECT_DOUBLE_EQ(stress(0, 1), 0.0); + EXPECT_DOUBLE_EQ(stress(0, 2), 0.0); + EXPECT_DOUBLE_EQ(stress(1, 0), 0.0); + EXPECT_DOUBLE_EQ(stress(1, 2), 0.0); + EXPECT_DOUBLE_EQ(stress(2, 0), 0.0); + EXPECT_DOUBLE_EQ(stress(2, 1), 0.0); +} + +// Test volume constraint rescaling in change_lattice +TEST_F(LatticeChangeBasicTest, ChangeLatticeVolumeRescaling) +{ + // Initialize variables + ucell.lc[0] = 1; + ucell.lc[1] = 1; + ucell.lc[2] = 1; + ucell.lat0 = 10.0; + + // Set initial lattice (cubic, volume = 1000) + ucell.latvec.e11 = 1.0; + ucell.latvec.e12 = 0.0; + ucell.latvec.e13 = 0.0; + ucell.latvec.e21 = 0.0; + ucell.latvec.e22 = 1.0; + ucell.latvec.e23 = 0.0; + ucell.latvec.e31 = 0.0; + ucell.latvec.e32 = 0.0; + ucell.latvec.e33 = 1.0; + + ucell.omega = 1000.0; // Initial volume + + lat[0] = 10.0; + lat[1] = 0.0; + lat[2] = 0.0; + lat[3] = 0.0; + lat[4] = 10.0; + lat[5] = 0.0; + lat[6] = 0.0; + lat[7] = 0.0; + lat[8] = 10.0; + + // Apply a move that would change volume (expand by 10%) + move[0] = 1.0; + move[1] = 0.0; + move[2] = 0.0; + move[3] = 0.0; + move[4] = 1.0; + move[5] = 0.0; + move[6] = 0.0; + move[7] = 0.0; + move[8] = 1.0; + + Lattice_Change_Basic::fixed_axes = "volume"; + + // Call change_lattice method + Lattice_Change_Basic::change_lattice(ucell, move, lat); + + // Check that volume is preserved (should still be 1000) + EXPECT_NEAR(ucell.omega, 1000.0, 1e-8); + + // Check that lattice vectors were rescaled uniformly + double expected_scale = std::pow(1000.0 / 1331.0, 1.0/3.0); // (old_vol / new_vol)^(1/3) + EXPECT_NEAR(ucell.latvec.e11, 1.1 * expected_scale, 1e-10); + EXPECT_NEAR(ucell.latvec.e22, 1.1 * expected_scale, 1e-10); + EXPECT_NEAR(ucell.latvec.e33, 1.1 * expected_scale, 1e-10); +} + +// Test volume constraint with non-cubic cell +TEST_F(LatticeChangeBasicTest, ChangeLatticeVolumeRescalingNonCubic) +{ + // Initialize variables + ucell.lc[0] = 1; + ucell.lc[1] = 1; + ucell.lc[2] = 1; + ucell.lat0 = 10.0; + + // Set initial lattice (non-cubic, volume = 1200) + ucell.latvec.e11 = 1.0; + ucell.latvec.e12 = 0.0; + ucell.latvec.e13 = 0.0; + ucell.latvec.e21 = 0.0; + ucell.latvec.e22 = 1.2; + ucell.latvec.e23 = 0.0; + ucell.latvec.e31 = 0.0; + ucell.latvec.e32 = 0.0; + ucell.latvec.e33 = 1.0; + + ucell.omega = 1200.0; // Initial volume + + lat[0] = 10.0; + lat[1] = 0.0; + lat[2] = 0.0; + lat[3] = 0.0; + lat[4] = 12.0; + lat[5] = 0.0; + lat[6] = 0.0; + lat[7] = 0.0; + lat[8] = 10.0; + + // Apply a move that would change volume + move[0] = 0.5; + move[1] = 0.0; + move[2] = 0.0; + move[3] = 0.0; + move[4] = -0.5; + move[5] = 0.0; + move[6] = 0.0; + move[7] = 0.0; + move[8] = 0.3; + + Lattice_Change_Basic::fixed_axes = "volume"; + + // Call change_lattice method + Lattice_Change_Basic::change_lattice(ucell, move, lat); + + // Check that volume is preserved + EXPECT_NEAR(ucell.omega, 1200.0, 1e-8); +} + +// Test change_lattice without volume constraint (should change volume) +TEST_F(LatticeChangeBasicTest, ChangeLatticeNoVolumeConstraint) +{ + // Initialize variables + ucell.lc[0] = 1; + ucell.lc[1] = 1; + ucell.lc[2] = 1; + ucell.lat0 = 10.0; + + // Set initial lattice + ucell.latvec.e11 = 1.0; + ucell.latvec.e12 = 0.0; + ucell.latvec.e13 = 0.0; + ucell.latvec.e21 = 0.0; + ucell.latvec.e22 = 1.0; + ucell.latvec.e23 = 0.0; + ucell.latvec.e31 = 0.0; + ucell.latvec.e32 = 0.0; + ucell.latvec.e33 = 1.0; + + ucell.omega = 1000.0; // Initial volume + + lat[0] = 10.0; + lat[1] = 0.0; + lat[2] = 0.0; + lat[3] = 0.0; + lat[4] = 10.0; + lat[5] = 0.0; + lat[6] = 0.0; + lat[7] = 0.0; + lat[8] = 10.0; + + // Apply a move that changes volume + move[0] = 1.0; + move[1] = 0.0; + move[2] = 0.0; + move[3] = 0.0; + move[4] = 1.0; + move[5] = 0.0; + move[6] = 0.0; + move[7] = 0.0; + move[8] = 1.0; + + Lattice_Change_Basic::fixed_axes = "None"; + + // Call change_lattice method + Lattice_Change_Basic::change_lattice(ucell, move, lat); + + // Check that volume DID change (should be 1331) + EXPECT_NEAR(ucell.omega, 1331.0, 1e-8); + + // Check lattice vectors + EXPECT_DOUBLE_EQ(ucell.latvec.e11, 1.1); + EXPECT_DOUBLE_EQ(ucell.latvec.e22, 1.1); + EXPECT_DOUBLE_EQ(ucell.latvec.e33, 1.1); +} + +// Test fixed_ibrav with simple cubic lattice +TEST_F(LatticeChangeBasicTest, ChangeLatticeFixedIbravSimpleCubic) +{ + // Initialize variables + ucell.lc[0] = 1; + ucell.lc[1] = 1; + ucell.lc[2] = 1; + ucell.lat0 = 10.0; + ucell.latName = "sc"; + + // Set initial lattice (slightly distorted cubic) + ucell.latvec.e11 = 1.0; + ucell.latvec.e12 = 0.01; // Small distortion + ucell.latvec.e13 = 0.0; + ucell.latvec.e21 = 0.0; + ucell.latvec.e22 = 1.0; + ucell.latvec.e23 = 0.0; + ucell.latvec.e31 = 0.0; + ucell.latvec.e32 = 0.0; + ucell.latvec.e33 = 1.0; + + lat[0] = 10.0; + lat[1] = 0.1; + lat[2] = 0.0; + lat[3] = 0.0; + lat[4] = 10.0; + lat[5] = 0.0; + lat[6] = 0.0; + lat[7] = 0.0; + lat[8] = 10.0; + + // Apply a move + move[0] = 0.1; + move[1] = 0.0; + move[2] = 0.0; + move[3] = 0.0; + move[4] = 0.1; + move[5] = 0.0; + move[6] = 0.0; + move[7] = 0.0; + move[8] = 0.1; + + PARAM.input.fixed_ibrav = true; + Lattice_Change_Basic::fixed_axes = "None"; + + // Verify remake_cell was not called yet + EXPECT_FALSE(unitcell::was_remake_cell_called()); + + // Call change_lattice method + Lattice_Change_Basic::change_lattice(ucell, move, lat); + + // Verify remake_cell was called + EXPECT_TRUE(unitcell::was_remake_cell_called()); + + // Check that lattice is now perfect cubic (all diagonal, equal) + // This is enforced by the mock remake_cell function + EXPECT_NEAR(ucell.latvec.e11, ucell.latvec.e22, 1e-10); + EXPECT_NEAR(ucell.latvec.e22, ucell.latvec.e33, 1e-10); + EXPECT_NEAR(ucell.latvec.e12, 0.0, 1e-10); + EXPECT_NEAR(ucell.latvec.e13, 0.0, 1e-10); + EXPECT_NEAR(ucell.latvec.e21, 0.0, 1e-10); + EXPECT_NEAR(ucell.latvec.e23, 0.0, 1e-10); + EXPECT_NEAR(ucell.latvec.e31, 0.0, 1e-10); + EXPECT_NEAR(ucell.latvec.e32, 0.0, 1e-10); + + // Reset for other tests + PARAM.input.fixed_ibrav = false; +} + +// Test fixed_ibrav with FCC lattice +TEST_F(LatticeChangeBasicTest, ChangeLatticeFixedIbravFCC) +{ + // Initialize variables + ucell.lc[0] = 1; + ucell.lc[1] = 1; + ucell.lc[2] = 1; + ucell.lat0 = 10.0; + ucell.latName = "fcc"; + + // Set initial lattice (slightly distorted FCC) + double celldm = 1.0; + ucell.latvec.e11 = -celldm + 0.01; // Small distortion + ucell.latvec.e12 = 0.0; + ucell.latvec.e13 = celldm; + ucell.latvec.e21 = 0.0; + ucell.latvec.e22 = celldm; + ucell.latvec.e23 = celldm; + ucell.latvec.e31 = -celldm; + ucell.latvec.e32 = celldm; + ucell.latvec.e33 = 0.0; + + lat[0] = (-celldm + 0.01) * ucell.lat0; + lat[1] = 0.0; + lat[2] = celldm * ucell.lat0; + lat[3] = 0.0; + lat[4] = celldm * ucell.lat0; + lat[5] = celldm * ucell.lat0; + lat[6] = -celldm * ucell.lat0; + lat[7] = celldm * ucell.lat0; + lat[8] = 0.0; + + // Apply a small move + for (int i = 0; i < 9; i++) move[i] = 0.01 * ucell.lat0; + + PARAM.input.fixed_ibrav = true; + Lattice_Change_Basic::fixed_axes = "None"; + + // Verify remake_cell was not called yet + EXPECT_FALSE(unitcell::was_remake_cell_called()); + + // Call change_lattice method + Lattice_Change_Basic::change_lattice(ucell, move, lat); + + // Verify remake_cell was called + EXPECT_TRUE(unitcell::was_remake_cell_called()); + + // Check that lattice maintains FCC structure + // For FCC: a1 = (-a, 0, a), a2 = (0, a, a), a3 = (-a, a, 0) + // All should have same magnitude + double mag1 = std::sqrt(ucell.latvec.e11*ucell.latvec.e11 + + ucell.latvec.e12*ucell.latvec.e12 + + ucell.latvec.e13*ucell.latvec.e13); + double mag2 = std::sqrt(ucell.latvec.e21*ucell.latvec.e21 + + ucell.latvec.e22*ucell.latvec.e22 + + ucell.latvec.e23*ucell.latvec.e23); + double mag3 = std::sqrt(ucell.latvec.e31*ucell.latvec.e31 + + ucell.latvec.e32*ucell.latvec.e32 + + ucell.latvec.e33*ucell.latvec.e33); + + EXPECT_NEAR(mag1, mag2, 1e-10); + EXPECT_NEAR(mag2, mag3, 1e-10); + + // Check FCC structure: e11 should be negative, e13 positive, e12 = 0 + EXPECT_LT(ucell.latvec.e11, 0.0); + EXPECT_GT(ucell.latvec.e13, 0.0); + EXPECT_NEAR(ucell.latvec.e12, 0.0, 1e-10); + + // Reset for other tests + PARAM.input.fixed_ibrav = false; +} + +// Test combination of fixed_axes = "volume" and fixed_ibrav +TEST_F(LatticeChangeBasicTest, ChangeLatticeVolumeAndIbrav) +{ + // Initialize variables + ucell.lc[0] = 1; + ucell.lc[1] = 1; + ucell.lc[2] = 1; + ucell.lat0 = 10.0; + ucell.latName = "sc"; + + // Set initial lattice (cubic) + ucell.latvec.e11 = 1.0; + ucell.latvec.e12 = 0.0; + ucell.latvec.e13 = 0.0; + ucell.latvec.e21 = 0.0; + ucell.latvec.e22 = 1.0; + ucell.latvec.e23 = 0.0; + ucell.latvec.e31 = 0.0; + ucell.latvec.e32 = 0.0; + ucell.latvec.e33 = 1.0; + + ucell.omega = 1000.0; // Initial volume + + lat[0] = 10.0; + lat[1] = 0.0; + lat[2] = 0.0; + lat[3] = 0.0; + lat[4] = 10.0; + lat[5] = 0.0; + lat[6] = 0.0; + lat[7] = 0.0; + lat[8] = 10.0; + + // Apply a move that would change volume and break symmetry + move[0] = 1.0; + move[1] = 0.1; // Try to break symmetry + move[2] = 0.0; + move[3] = 0.0; + move[4] = 0.8; + move[5] = 0.0; + move[6] = 0.0; + move[7] = 0.0; + move[8] = 1.2; + + PARAM.input.fixed_ibrav = true; + Lattice_Change_Basic::fixed_axes = "volume"; + + // Verify remake_cell was not called yet + EXPECT_FALSE(unitcell::was_remake_cell_called()); + + // Call change_lattice method + Lattice_Change_Basic::change_lattice(ucell, move, lat); + + // Verify remake_cell was called (should be called before volume rescaling) + EXPECT_TRUE(unitcell::was_remake_cell_called()); + + // Check that volume is preserved + EXPECT_NEAR(ucell.omega, 1000.0, 1e-8); + + // Check that lattice is cubic (all diagonal, equal) + EXPECT_NEAR(ucell.latvec.e11, ucell.latvec.e22, 1e-10); + EXPECT_NEAR(ucell.latvec.e22, ucell.latvec.e33, 1e-10); + EXPECT_NEAR(ucell.latvec.e12, 0.0, 1e-10); + EXPECT_NEAR(ucell.latvec.e13, 0.0, 1e-10); + EXPECT_NEAR(ucell.latvec.e21, 0.0, 1e-10); + EXPECT_NEAR(ucell.latvec.e23, 0.0, 1e-10); + EXPECT_NEAR(ucell.latvec.e31, 0.0, 1e-10); + EXPECT_NEAR(ucell.latvec.e32, 0.0, 1e-10); + + // Reset for other tests + PARAM.input.fixed_ibrav = false; +} + +// Test axis constraint with fixed_axes = "a" +TEST_F(LatticeChangeBasicTest, SetupGradientAxisA) +{ + // Initialize variables + ucell.lc[0] = 0; // First axis fixed + ucell.lc[1] = 1; + ucell.lc[2] = 1; + + stress(0, 0) = 1.0; + stress(0, 1) = 2.0; + stress(0, 2) = 3.0; + stress(1, 0) = 4.0; + stress(1, 1) = 5.0; + stress(1, 2) = 6.0; + stress(2, 0) = 7.0; + stress(2, 1) = 8.0; + stress(2, 2) = 9.0; + + Lattice_Change_Basic::fixed_axes = "a"; + + // Call setup_gradient method + Lattice_Change_Basic::setup_gradient(ucell, lat, grad, stress); + + // Check that gradient for first lattice vector is zero + EXPECT_DOUBLE_EQ(grad[0], 0.0); + EXPECT_DOUBLE_EQ(grad[1], 0.0); + EXPECT_DOUBLE_EQ(grad[2], 0.0); + + // Check that gradients for other vectors are non-zero + EXPECT_NE(grad[3], 0.0); + EXPECT_NE(grad[4], 0.0); + EXPECT_NE(grad[5], 0.0); + EXPECT_NE(grad[6], 0.0); + EXPECT_NE(grad[7], 0.0); + EXPECT_NE(grad[8], 0.0); +} + +// Test that fixed axis doesn't move in change_lattice +TEST_F(LatticeChangeBasicTest, ChangeLatticeFixedAxisA) +{ + // Initialize variables + ucell.lc[0] = 0; // First axis fixed + ucell.lc[1] = 1; + ucell.lc[2] = 1; + ucell.lat0 = 10.0; + + // Set initial lattice + ucell.latvec.e11 = 1.0; + ucell.latvec.e12 = 0.1; + ucell.latvec.e13 = 0.2; + ucell.latvec.e21 = 0.0; + ucell.latvec.e22 = 1.0; + ucell.latvec.e23 = 0.0; + ucell.latvec.e31 = 0.0; + ucell.latvec.e32 = 0.0; + ucell.latvec.e33 = 1.0; + + // Save initial first lattice vector + double initial_e11 = ucell.latvec.e11; + double initial_e12 = ucell.latvec.e12; + double initial_e13 = ucell.latvec.e13; + + lat[0] = 10.0; + lat[1] = 1.0; + lat[2] = 2.0; + lat[3] = 0.0; + lat[4] = 10.0; + lat[5] = 0.0; + lat[6] = 0.0; + lat[7] = 0.0; + lat[8] = 10.0; + + // Try to move all lattice vectors + move[0] = 1.0; + move[1] = 0.5; + move[2] = 0.3; + move[3] = 0.5; + move[4] = 1.0; + move[5] = 0.0; + move[6] = 0.0; + move[7] = 0.0; + move[8] = 1.0; + + Lattice_Change_Basic::fixed_axes = "a"; + + // Call change_lattice method + Lattice_Change_Basic::change_lattice(ucell, move, lat); + + // Check that first lattice vector didn't change + EXPECT_DOUBLE_EQ(ucell.latvec.e11, initial_e11); + EXPECT_DOUBLE_EQ(ucell.latvec.e12, initial_e12); + EXPECT_DOUBLE_EQ(ucell.latvec.e13, initial_e13); + + // Check that other lattice vectors did change + EXPECT_NE(ucell.latvec.e21, 0.0); + EXPECT_NE(ucell.latvec.e22, 1.0); + EXPECT_NE(ucell.latvec.e33, 1.0); +} + +// Test that remake_cell is NOT called when fixed_ibrav = false +TEST_F(LatticeChangeBasicTest, ChangeLatticeNoFixedIbrav) +{ + // Initialize variables + ucell.lc[0] = 1; + ucell.lc[1] = 1; + ucell.lc[2] = 1; + ucell.lat0 = 10.0; + ucell.latName = "sc"; + + // Set initial lattice + ucell.latvec.e11 = 1.0; + ucell.latvec.e12 = 0.01; // Small distortion + ucell.latvec.e13 = 0.0; + ucell.latvec.e21 = 0.0; + ucell.latvec.e22 = 1.0; + ucell.latvec.e23 = 0.0; + ucell.latvec.e31 = 0.0; + ucell.latvec.e32 = 0.0; + ucell.latvec.e33 = 1.0; + + lat[0] = 10.0; + lat[1] = 0.1; + lat[2] = 0.0; + lat[3] = 0.0; + lat[4] = 10.0; + lat[5] = 0.0; + lat[6] = 0.0; + lat[7] = 0.0; + lat[8] = 10.0; + + // Apply a move + move[0] = 0.1; + move[1] = 0.0; + move[2] = 0.0; + move[3] = 0.0; + move[4] = 0.1; + move[5] = 0.0; + move[6] = 0.0; + move[7] = 0.0; + move[8] = 0.1; + + PARAM.input.fixed_ibrav = false; // Explicitly set to false + Lattice_Change_Basic::fixed_axes = "None"; + + // Verify remake_cell was not called yet + EXPECT_FALSE(unitcell::was_remake_cell_called()); + + // Call change_lattice method + Lattice_Change_Basic::change_lattice(ucell, move, lat); + + // Verify remake_cell was NOT called + EXPECT_FALSE(unitcell::was_remake_cell_called()); + + // Check that distortion remains (e12 should still be non-zero) + EXPECT_NE(ucell.latvec.e12, 0.0); +} diff --git a/source/source_relax/test/lattice_change_cg_test.cpp b/source/source_relax/test/lattice_change_cg_test.cpp index cdea39282b..01c863b8ee 100644 --- a/source/source_relax/test/lattice_change_cg_test.cpp +++ b/source/source_relax/test/lattice_change_cg_test.cpp @@ -1,5 +1,6 @@ #include "for_test.h" #include "gtest/gtest.h" +#include "mock_remake_cell.h" #define private public #include "source_relax/lattice_change_basic.h" #include "source_relax/lattice_change_cg.h" diff --git a/source/source_relax/test/lattice_change_methods_test.cpp b/source/source_relax/test/lattice_change_methods_test.cpp index 24d9b515ad..64173e110e 100644 --- a/source/source_relax/test/lattice_change_methods_test.cpp +++ b/source/source_relax/test/lattice_change_methods_test.cpp @@ -1,4 +1,5 @@ #include "source_relax/lattice_change_methods.h" +#include "mock_remake_cell.h" #include "for_test.h" #include "gtest/gtest.h" diff --git a/source/source_relax/test/mock_remake_cell.cpp b/source/source_relax/test/mock_remake_cell.cpp new file mode 100644 index 0000000000..ee51ba8653 --- /dev/null +++ b/source/source_relax/test/mock_remake_cell.cpp @@ -0,0 +1,58 @@ +#include "mock_remake_cell.h" +#include + +namespace unitcell +{ + // Mock state variables + bool remake_cell_called = false; + std::string remake_cell_latName = ""; + ModuleBase::Matrix3 remake_cell_latvec; + + void remake_cell(Lattice& lat) + { + remake_cell_called = true; + remake_cell_latName = lat.latName; + remake_cell_latvec = lat.latvec; + + // Mock implementation: enforce simple cubic structure for "sc" + if (lat.latName == "sc") + { + double celldm = std::sqrt(lat.latvec.e11 * lat.latvec.e11 + + lat.latvec.e12 * lat.latvec.e12 + + lat.latvec.e13 * lat.latvec.e13); + lat.latvec.Zero(); + lat.latvec.e11 = celldm; + lat.latvec.e22 = celldm; + lat.latvec.e33 = celldm; + } + // Mock implementation: enforce FCC structure for "fcc" + else if (lat.latName == "fcc") + { + double celldm = std::sqrt(lat.latvec.e11 * lat.latvec.e11 + + lat.latvec.e12 * lat.latvec.e12 + + lat.latvec.e13 * lat.latvec.e13) / std::sqrt(2.0); + lat.latvec.e11 = -celldm; + lat.latvec.e12 = 0.0; + lat.latvec.e13 = celldm; + lat.latvec.e21 = 0.0; + lat.latvec.e22 = celldm; + lat.latvec.e23 = celldm; + lat.latvec.e31 = -celldm; + lat.latvec.e32 = celldm; + lat.latvec.e33 = 0.0; + } + // For other lattice types, do nothing (just track the call) + } + + void reset_remake_cell_mock() + { + remake_cell_called = false; + remake_cell_latName = ""; + remake_cell_latvec.Zero(); + } + + bool was_remake_cell_called() + { + return remake_cell_called; + } +} diff --git a/source/source_relax/test/mock_remake_cell.h b/source/source_relax/test/mock_remake_cell.h new file mode 100644 index 0000000000..c2d4d63bd8 --- /dev/null +++ b/source/source_relax/test/mock_remake_cell.h @@ -0,0 +1,22 @@ +#ifndef MOCK_REMAKE_CELL_H +#define MOCK_REMAKE_CELL_H + +#include "source_cell/unitcell_data.h" +#include + +namespace unitcell +{ + // Mock state tracking + extern bool remake_cell_called; + extern std::string remake_cell_latName; + extern ModuleBase::Matrix3 remake_cell_latvec; + + // Mock implementation of remake_cell + void remake_cell(Lattice& lat); + + // Helper functions for testing + void reset_remake_cell_mock(); + bool was_remake_cell_called(); +} + +#endif // MOCK_REMAKE_CELL_H