From e11a5df1ad19f795bdff06386f91a679bfbe9ef9 Mon Sep 17 00:00:00 2001 From: dyzheng Date: Sun, 2 Aug 2026 11:06:13 +0800 Subject: [PATCH 1/3] Fix(xc): fix libxc gga_grad branching for nspin=4 GGA The gga_grad=3 path (default for nspin=4) incorrectly used the collinear approximation (gga_grad=1) within libxc instead of the Scalmani-Frisch gradient method (gga_grad=2). This caused inconsistent results for noncollinear GGA calculations. Fix: extend libxc_pot.cpp and libxc_tools.cpp gga_grad==2 checks to also match gga_grad==3, routing both through the SF gradient path. Also fix input_parameter.h: restore out_hsr/out_hsk/out_hsr_npz_compat members removed during cherry-pick that were present in newer develop. --- docs/advanced/input_files/input-main.md | 145 +++------ docs/advanced/scf/spin.md | 265 +++++++++++++++ docs/parameters.yaml | 142 ++++---- source/source_hamilt/module_xc/CMakeLists.txt | 1 + source/source_hamilt/module_xc/libxc_abacus.h | 80 +++-- source/source_hamilt/module_xc/libxc_pot.cpp | 25 +- .../source_hamilt/module_xc/libxc_setup.cpp | 77 ++--- .../source_hamilt/module_xc/libxc_tools.cpp | 225 ++++++++++++- .../module_xc/test/CMakeLists.txt | 1 + .../source_hamilt/module_xc/test/test_xc3.cpp | 3 +- .../source_hamilt/module_xc/test/test_xc5.cpp | 7 +- .../source_hamilt/module_xc/test/xc3_mock.h | 2 +- .../source_hamilt/module_xc/xc_functional.h | 16 +- .../xc_functional_gga_noncol_sf_builtin.cpp | 308 ++++++++++++++++++ .../xc_functional_gga_noncol_sf_builtin.h | 30 ++ source/source_hamilt/module_xc/xc_pot.cpp | 9 + .../module_parameter/input_parameter.h | 27 +- .../read_input_item_elec_stru.cpp | 76 +++-- 18 files changed, 1134 insertions(+), 305 deletions(-) create mode 100644 source/source_hamilt/module_xc/xc_functional_gga_noncol_sf_builtin.cpp create mode 100644 source/source_hamilt/module_xc/xc_functional_gga_noncol_sf_builtin.h diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index f04f903448..94aa78b935 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -168,9 +168,7 @@ - [out\_proj\_band](#out_proj_band) - [out\_stru](#out_stru) - [out\_level](#out_level) - - [out\_hsk](#out_hsk) - [out\_mat\_hs](#out_mat_hs) - - [out\_hsr](#out_hsr) - [out\_mat\_hs2](#out_mat_hs2) - [out\_mat\_tk](#out_mat_tk) - [out\_mat\_r](#out_mat_r) @@ -410,9 +408,6 @@ - [sc\_scf\_thr](#sc_scf_thr) - [sc\_direction\_only](#sc_direction_only) - [sc\_lambda\_strategy](#sc_lambda_strategy) - - [sc\_scan\_lambda\_start](#sc_scan_lambda_start) - - [sc\_scan\_lambda\_end](#sc_scan_lambda_end) - - [sc\_scan\_steps](#sc_scan_steps) - [vdW correction](#vdw-correction) - [vdw\_method](#vdw_method) - [vdw\_d4\_xc](#vdw_d4_xc) @@ -587,7 +582,7 @@ - md: perform molecular dynamics simulations - get_pchg: obtain partial (band-decomposed) charge densities (for LCAO basis only). See out_pchg for more information - get_wf: obtain real space wave functions (for LCAO basis only). See out_wfc_norm and out_wfc_re_im for more information - - get_s: obtain the overlap matrix formed by localized orbitals (for LCAO basis with multiple k points). The file name is OUT.${suffix}/sr_nao.csr, with the same file format as generated by out_hsr 1 + - get_s: obtain the overlap matrix formed by localized orbitals (for LCAO basis with multiple k points). the file name is SR.csr with file format being the same as that generated by out_mat_hs2 - gen_bessel: generates projectors, i.e., a series of Bessel functions, for the DeePKS method (for LCAO basis only) - gen_opt_abfs: generate opt-ABFs as discussed in this article - test_memory: obtain a rough estimation of memory consumption for the calculation @@ -616,7 +611,7 @@ - **Description**: Takes value 1, 0 or -1. - -1: No symmetry will be considered. It is recommended to set -1 for non-colinear + soc calculations, where time reversal symmetry is broken sometimes. - 0: Only time reversal symmetry would be considered in symmetry operations, which implied k point and -k point would be treated as a single k point with twice the weight. - - 1: Symmetry analysis will be performed to determine the type of Bravais lattice and associated symmetry operations (point groups, space groups, primitive cells, and irreducible k-points). For a magnetic system, the symmetry of the initial magnetic structure will be analyzed and preserved. + - 1: Symmetry analysis will be performed to determine the type of Bravais lattice and associated symmetry operations. (point groups, space groups, primitive cells, and irreducible k-points) > Note: When symmetry is enabled (value 1), k-points are reduced to the irreducible Brillouin zone (IBZ). For explicit k-point lists with custom weights (see KPT file), the custom weights are preserved during symmetry reduction. For Monkhorst-Pack grids, uniform weights are used. - **Default**: default @@ -1154,9 +1149,9 @@ For plane-wave basis, - cg: The conjugate-gradient (CG) method. + - bpcg: The BPCG method, which is a block-parallel Conjugate Gradient (CG) method, typically exhibits higher acceleration in a GPU environment. - dav: The Davidson algorithm. - dav_subspace: The Davidson algorithm without orthogonalization operation, this method is the most recommended for efficiency. `pw_diag_ndim` can be set to 2 for this method. - - bpcg: The BPCG method, which is a block-parallel Conjugate Gradient (CG) method, typically exhibits higher acceleration in a GPU environment. The BPCG method is currently under testing and is not recommended for use. For numerical atomic orbitals basis, @@ -1234,16 +1229,18 @@ ### xc_exch_ext - **Type**: Integer followed by Real values -- **Description**: Customized parameterization of the exchange part of an XC functional. The first value should be the Libxc ID of the original functional, followed by the complete list of external parameters required by the linked Libxc version. If unset, Libxc's own default parameters are used. For functional IDs and parameter definitions, refer to the Libxc documentation and source code. +- **Description**: Customized parameterization on the exchange part of XC functional. The first value should be the LibXC ID of the original functional, and latter values are external parameters. Default values are those of Perdew-Burke-Ernzerhof (PBE) functional. For more information on LibXC ID of functionals, please refer to LibXC. For parameters of functionals of interest, please refer to the source code of LibXC, such as PBE functional interface in LibXC: gga_x_pbe.c. > Note: Solely setting this keyword will take no effect on XC functionals. One should also set dft_functional to the corresponding functional to apply the customized parameterization. Presently this feature can only support parameterization on one exchange functional. +- **Default**: 101 0.8040 0.2195149727645171 ### xc_corr_ext - **Type**: Integer followed by Real values -- **Description**: Customized parameterization of the correlation part of an XC functional. The first value should be the Libxc ID of the original functional, followed by the complete list of external parameters required by the linked Libxc version. If unset, Libxc's own default parameters are used. For functional IDs and parameter definitions, refer to the Libxc documentation and source code. +- **Description**: Customized parameterization on the correlation part of XC functional. The first value should be the LibXC ID of the original functional, and latter values are external parameters. Default values are those of Perdew-Burke-Ernzerhof (PBE) functional. For more information on LibXC ID of functionals, please refer to LibXC. For parameters of functionals of interest, please refer to the source code of LibXC, such as PBE functional interface in LibXC: gga_c_pbe.c. > Note: Solely setting this keyword will take no effect on XC functionals. One should also set dft_functional to the corresponding functional to apply the customized parameterization. Presently this feature can only support parameterization on one correlation functional. +- **Default**: 130 0.06672455060314922 0.031090690869654895034 1.0 ### pseudo_rcut @@ -1909,7 +1906,7 @@ - **Type**: Boolean \[Integer\](optional) - **Availability**: *Numerical atomic orbital basis (multi-k points)* -- **Description**: Whether to output the density matrix with Bravias lattice vector R index into files in the folder OUT.${suffix}. The files are named as dmr{s}{spin index}{g}{geometry index}{_nao} + {".csr"}. Here, 's' refers to spin, where s1 means spin up channel while s2 means spin down channel, and the sparse matrix format 'csr' is mentioned in out_hsr. Finally, if out_app_flag is set to false, the file name contains the optional 'g' index for each ionic step that may have different geometries, and if out_app_flag is set to true, the density matrix with respect to Bravias lattice vector R accumulates during ionic steps: +- **Description**: Whether to output the density matrix with Bravias lattice vector R index into files in the folder OUT.${suffix}. The files are named as dmr{s}{spin index}{g}{geometry index}{_nao} + {".csr"}. Here, 's' refers to spin, where s1 means spin up channel while s2 means spin down channel, and the sparse matrix format 'csr' is mentioned in out_mat_hs2. Finally, if out_app_flag is set to false, the file name contains the optional 'g' index for each ionic step that may have different geometries, and if out_app_flag is set to true, the density matrix with respect to Bravias lattice vector R accumulates during ionic steps: - nspin = 1: dmrs1_nao.csr; - nspin = 2: dmrs1_nao.csr and dmrs2_nao.csr for the two spin channels. @@ -2006,60 +2003,32 @@ - m: molecular dynamics level, which does not print some information for simplicity. - **Default**: ie -### out_hsk - -- **Type**: Integer \[Integer\](optional) -- **Availability**: *Numerical atomic orbital basis* -- **Description**: Output the upper triangular part of the Hamiltonian and overlap matrices in reciprocal space for each k-point into files in the directory OUT.${suffix}. The first integer selects the format: - - 0: disabled; - - 1: text output; the optional second integer controls precision and defaults to 8; - - 2: reserved for binary output, which is not implemented yet; - - 3: NPZ output, which is not implemented for H(k)/S(k). - - The output is also controlled by out_freq_ion and out_app_flag. For more information, refer to hs_matrix.md. - - - Gamma-only, nspin = 1: hk_nao.txt for the Hamiltonian matrix and sk_nao.txt for the overlap matrix. - - Gamma-only, nspin = 2: hks1_nao.txt and hks2_nao.txt for the two spin channels of the Hamiltonian matrix, and sk_nao.txt for the overlap matrix. Only one overlap matrix is written because it is identical for both spin channels. - - Gamma-only, nspin = 4: not available with the gamma-only algorithm. - - Multi-k, nspin = 1: hk1_nao.txt for the Hamiltonian matrix and sk1_nao.txt for the overlap matrix at the first k-point. - - Multi-k, nspin = 2: hk1s1_nao.txt and hk1s2_nao.txt for the two spin channels of the Hamiltonian matrix, and sk1_nao.txt for the overlap matrix at the first k-point. Only one overlap matrix is written because it is identical for both spin channels. - - Multi-k, nspin = 4: hk1s4_nao.txt for the spinor Hamiltonian matrix and sk1_nao.txt for the spinor overlap matrix at the first k-point. - When out_app_flag is false, g followed by the one-based ionic-step index is inserted before _nao, for example hk1s1g1_nao.txt. - - > Note: In the 3.10-LTS version, the file names are data-0-H and data-0-S, etc. -- **Default**: 0 8 -- **Unit**: Ry - ### out_mat_hs - **Type**: Boolean \[Integer\](optional) - **Availability**: *Numerical atomic orbital basis* -- **Description**: Legacy alias for out_hsk 1, which outputs Hamiltonian and overlap matrices in reciprocal space for each k-point. The optional second integer controls text precision. If both out_hsk and out_mat_hs are present, out_hsk takes precedence. -- **Default**: False 8 -- **Unit**: Ry - -### out_hsr - -- **Type**: Integer \[Integer\](optional) -- **Availability**: *Numerical atomic orbital basis* -- **Description**: Output Hamiltonian and overlap matrices in real space, indexed by the Bravais lattice vector R, in the directory OUT.${suffix}. The first integer selects the format: - - 0: disabled; - - 1: text CSR output; the optional second integer controls precision and defaults to 8; - - 2: reserved for binary output, which is not implemented yet; - - 3: NPZ output using hrs1_nao.npz, hrs2_nao.npz when needed, and sr_nao.npz. - - For multi-k calculations, the output contains the individual real-space blocks stored for the Bravais lattice vectors R. For gamma-only calculations, the internal real-space contributions are folded into a single R = (0, 0, 0) block. This folded result cannot recover the original R-resolved contributions or interpolate arbitrary k points. Terms added only while constructing H(k) are not guaranteed to be present. +- **Description**: Whether to print the upper triangular part of the Hamiltonian matrices and overlap matrices for each k-point into files in the directory OUT.${suffix}. The second number controls precision. For more information, please refer to hs_matrix.md. Also controled by out_freq_ion and out_app_flag. + - For gamma only case: + - nspin = 1: hks1_nao.txt for the Hamiltonian matrix and sks1_nao.txt for the overlap matrix; + - nspin = 2: hks1_nao.txt and hks2_nao.txt for the Hamiltonian matrix and sks1_nao.txt for the overlap matrix. Note that the code will not output sks2_nao.txt because it is the same as sks1_nao.txt; + - nspin = 4: hks12_nao.txt for the Hamiltonian matrix and sks12_nao.txt for the overlap matrix. + - For multi-k points case: + - nspin = 1: hks1k1_nao.txt for the Hamiltonian matrix at the 1st k-point, and sks1k1_nao.txt for the overlap matrix for the 1st k-point, ...; + - nspin = 2: hks1k1_nao.txt and hks2k1_nao.txt for the two spin channels of the Hamiltonian matrix at the 1st k-point, and sks1k1_nao.txt for the overlap matrix for the 1st k-point. Note that the code will not output sks2k1_nao.txt because it is the same as sks1k1_nao.txt, ...; + - nspin = 4: hks12k1_nao.txt for the Hamiltonian matrix at the 1st k-point, and sks12k1_nao.txt for the overlap matrix for the 1st k-point, ...; - > Note: In the 3.10-LTS version, the file names are data-HR-sparse_SPIN0.csr and data-SR-sparse_SPIN0.csr, etc. -- **Default**: 0 8 + > Note: In the 3.10-LTS version, the file names are data-0-H and data-0-S, etc. +- **Default**: False 8 - **Unit**: Ry ### out_mat_hs2 - **Type**: Boolean \[Integer\](optional) -- **Availability**: *Numerical atomic orbital basis* -- **Description**: Legacy alias for out_hsr 1, which outputs Hamiltonian and overlap matrices in real space indexed by the Bravais lattice vector R. The optional second integer controls text precision. If both out_hsr and out_mat_hs2 are present, out_hsr takes precedence. -- **Default**: False 8 +- **Availability**: *Numerical atomic orbital basis (not gamma-only algorithm)* +- **Description**: Whether to print files containing the Hamiltonian matrix and overlap matrix into files in the directory OUT.${suffix}. For more information, please refer to hs_matrix.md. + + > Note: In the 3.10-LTS version, the file names are data-HR-sparse_SPIN0.csr and data-SR-sparse_SPIN0.csr, etc. +- **Default**: False [8] - **Unit**: Ry ### out_mat_tk @@ -2086,7 +2055,7 @@ - **Type**: Boolean \[Integer\](optional) - **Availability**: *Numerical atomic orbital basis (not gamma-only algorithm)* -- **Description**: Generate files containing the kinetic energy matrix. The optional second parameter controls text output precision. The format will be the same as the Hamiltonian matrix and overlap matrix as mentioned in out_hsr. The name of the files will be trs1_nao.csr and so on. Also controled by out_freq_ion and out_app_flag. +- **Description**: Generate files containing the kinetic energy matrix. The optional second parameter controls text output precision. The format will be the same as the Hamiltonian matrix and overlap matrix as mentioned in out_mat_hs2. The name of the files will be trs1_nao.csr and so on. Also controled by out_freq_ion and out_app_flag. > Note: In the 3.10-LTS version, the file name is data-TR-sparse_SPIN0.csr. - **Default**: False 8 @@ -2096,7 +2065,7 @@ - **Type**: Integer - **Availability**: *Numerical atomic orbital basis (not gamma-only algorithm)* -- **Description**: Whether to print files containing the derivatives of the Hamiltonian matrix. The format will be the same as the Hamiltonian matrix and overlap matrix as mentioned in out_hsr. The name of the files will be dhrxs1_nao.csr, dhrys1_nao.csr, dhrzs1_nao.csr and so on. Also controled by out_freq_ion and out_app_flag. +- **Description**: Whether to print files containing the derivatives of the Hamiltonian matrix. The format will be the same as the Hamiltonian matrix and overlap matrix as mentioned in out_mat_hs2. The name of the files will be dhrxs1_nao.csr, dhrys1_nao.csr, dhrzs1_nao.csr and so on. Also controled by out_freq_ion and out_app_flag. Format: <enable> [precision] [iat1 iat2 ...]. The first value (0/1) enables/disables output. The second optional value sets the output precision (default: 8). Starting from the third value, 1-based atom indices can be listed to restrict output to derivatives with respect to those specific atoms only; if no atom indices are given, all atoms are written. @@ -2163,7 +2132,7 @@ - **Type**: Integer - **Description**: Whether to print files containing the kinetic energy matrix T(R) in CSR format. - See out_hsr for format details. + See out_mat_hs2 for format details. - **Default**: 0 8 - **Unit**: Ry @@ -2172,7 +2141,7 @@ - **Type**: Integer - **Description**: Whether to print files containing the nonlocal pseudopotential matrix Vnl(R) in CSR format. - See out_hsr for format details. + See out_mat_hs2 for format details. - **Default**: 0 8 - **Unit**: Ry @@ -2181,7 +2150,7 @@ - **Type**: Integer - **Description**: Whether to print files containing the local pseudopotential matrix Vl(R) in CSR format. - See out_hsr for format details. + See out_mat_hs2 for format details. - **Default**: 0 8 - **Unit**: Ry @@ -2190,7 +2159,7 @@ - **Type**: Integer - **Description**: Whether to print files containing the Hartree matrix Vh(R) in CSR format. - See out_hsr for format details. + See out_mat_hs2 for format details. - **Default**: 0 8 - **Unit**: Ry @@ -2199,7 +2168,7 @@ - **Type**: Integer - **Description**: Whether to print files containing the XC matrix Vxc(R) in CSR format. - See out_hsr for format details. + See out_mat_hs2 for format details. - **Default**: 0 8 - **Unit**: Ry @@ -2208,7 +2177,7 @@ - **Type**: Integer - **Description**: Whether to print files containing the exact-exchange matrix Vexx(R) in CSR format. - See out_hsr for format details. + See out_mat_hs2 for format details. - **Default**: 0 8 - **Unit**: Ry @@ -2275,15 +2244,15 @@ - **Type**: Boolean - **Availability**: *Numerical atomic orbital basis (not gamma-only algorithm)* -- **Description**: Whether to print Hamiltonian matrices H(R) in NPZ format as hrs1_nao.npz and, for nspin = 2, hrs2_nao.npz. This feature does not work for gamma-only calculations. +- **Description**: Whether to print Hamiltonian matrices H(R) in npz format. This feature does not work for gamma-only calculations. - **Default**: False - **Unit**: Ry ### out_hsr_npz - **Type**: Boolean -- **Availability**: *Numerical atomic orbital basis* -- **Description**: Legacy alias for out_hsr 3, writing hrs1_nao.npz, hrs2_nao.npz when needed, and sr_nao.npz. If both out_hsr and out_hsr_npz are present, out_hsr takes precedence. Gamma-only calculations write the folded R = (0, 0, 0) representation. +- **Availability**: *Numerical atomic orbital basis (not gamma-only algorithm)* +- **Description**: Whether to print Hamiltonian matrices H(R) and overlap matrix S(R) in npz format. This feature does not work for gamma-only calculations. - **Default**: False - **Unit**: Ry @@ -2305,13 +2274,13 @@ - **Type**: Boolean - **Availability**: *Numerical atomic orbital basis (not gamma-only algorithm)* -- **Description**: Whether to output r(R), H(R), S(R), T(R), dH(R), dS(R), and wfc matrices in an append manner during molecular dynamics calculations. Check input parameters out_mat_r, out_hsr, out_mat_t, out_mat_dh, out_hsk and out_wfc_lcao for more information. +- **Description**: Whether to output r(R), H(R), S(R), T(R), dH(R), dS(R), and wfc matrices in an append manner during molecular dynamics calculations. Check input parameters out_mat_r, out_mat_hs2, out_mat_t, out_mat_dh, out_mat_hs and out_wfc_lcao for more information. - **Default**: true ### out_ndigits - **Type**: Integer -- **Availability**: *out_hsk 1 case presently.* +- **Availability**: *out_mat_hs 1 case presently.* - **Description**: Controls the length of decimal part of output data, such as charge density, Hamiltonian matrix, Overlap matrix and so on. - **Default**: 8 @@ -3707,8 +3676,8 @@ - **Type**: Integer - **Description**: Determines whether to calculate the plus U correction, which is especially important for correlated electrons. - - 1: Calculate plus U correction with radius-adjustable localized projections (with parameter onsite_radius). - - 2: Calculate plus U correction using first zeta of NAOs as projections (this is old method for testing). + - 1: Calculate plus U correction with radius-adjustable localized projections (with parameter onsite_radius). Supported for both PW and LCAO basis sets. + - 2: Calculate plus U correction using first zeta of NAOs as projections (this is old method for testing). Only available for LCAO basis. - 0: Do not calculate plus U correction. - **Default**: 0 @@ -3852,46 +3821,20 @@ - **Type**: Boolean - **Availability**: *sc_mag_switch is true* -- **Description**: When true, only the direction of the magnetic moment is constrained to the target direction, while the magnitude is allowed to vary freely. This is useful for studying magnetic anisotropy or when the magnitude of the moment is determined by the electronic structure rather than an external constraint. - - When false (default), both the direction and magnitude of the magnetic moment are constrained to the target values. +- **Description**: When true, only the direction of the magnetic moment is constrained to the target direction, while the magnitude is allowed to vary freely. This is useful for studying magnetic anisotropy or when the magnitude of the moment is determined by the electronic structure rather than an external constraint. When false (default), both the direction and magnitude of the magnetic moment are constrained to the target values. - **Default**: False ### sc_lambda_strategy - **Type**: String - **Availability**: *sc_mag_switch is true* -- **Description**: Lambda update strategy for spin-constrained DFT: - - bfgs: BFGS quasi-Newton method - - linear_response: linear response (Scheme B) - - augmented_lagrangian: augmented Lagrangian (Scheme C) +- **Description**: Lambda update strategy for spin-constrained DFT. Available options are: + - bfgs: BFGS quasi-Newton method (default, robust and well-tested) + - linear_response: linear response method (Scheme B) + - augmented_lagrangian: augmented Lagrangian method (Scheme C) - hybrid_delayed: hybrid delayed update (Scheme D) - - linear_scan: linear sweep of lambda for testing magnetic moment response - **Default**: bfgs -### sc_scan_lambda_start - -- **Type**: Float -- **Availability**: *sc_lambda_strategy is linear_scan* -- **Description**: Starting lambda value for linear_scan strategy. Only used when sc_lambda_strategy=linear_scan. -- **Default**: 0.0 -- **Unit**: eV/uB - -### sc_scan_lambda_end - -- **Type**: Float -- **Availability**: *sc_lambda_strategy is linear_scan* -- **Description**: Ending lambda value for linear_scan strategy. Only used when sc_lambda_strategy=linear_scan. -- **Default**: 1.0 -- **Unit**: eV/uB - -### sc_scan_steps - -- **Type**: Integer -- **Availability**: *sc_lambda_strategy is linear_scan* -- **Description**: Number of lambda values to scan. Only used when sc_lambda_strategy=linear_scan. -- **Default**: 20 - [back to top](#full-list-of-input-keywords) ## vdW correction diff --git a/docs/advanced/scf/spin.md b/docs/advanced/scf/spin.md index 1749db156d..10c783d868 100644 --- a/docs/advanced/scf/spin.md +++ b/docs/advanced/scf/spin.md @@ -28,6 +28,224 @@ If **"ocp=1"** and **"ocp_set"** is set in INPUT file, the occupations of states 2. **"nupdown"** If **"nupdown"** is set to non-zero, number of spin-up and spin-down electrons will be fixed, and Fermi energy level will split to E_Fermi_up and E_Fermi_down. By the way, total magnetization will also be fixed, and will be the value of **"nupdown"**. +## DeltaSpin (Spin-Constrained DFT) + +DeltaSpin is a spin-constrained DFT method that allows users to constrain the magnetic moments on individual atoms to target values during self-consistent field (SCF) calculations. This is useful for studying magnetic excitations, non-collinear magnetic structures, and systems where the magnetic ground state is not known a priori. + +The theoretical foundation and implementation details can be found in: + +- Cai Z, Wang K, Xu Y, et al., "A self-adaptive first-principles approach for magnetic excited states," *Quantum Frontiers* 2.1 (2023): 21. [DOI: 10.1007/s44214-023-00050-z](https://doi.org/10.1007/s44214-023-00050-z) +- Zheng D, Peng X, Huang Y, et al., "Integrating deep-learning-based magnetic model and non-collinear spin-constrained method: methodology, implementation and application," *npj Computational Materials* (2026). + +### Enabling DeltaSpin + +Set `sc_mag_switch 1` in the INPUT file. DeltaSpin is supported for both PW (`basis_type = pw`) and LCAO (`basis_type = lcao`) basis sets, with `nspin = 2` (collinear) or `nspin = 4` (non-collinear). + +### Specifying Target Magnetic Moments in STRU + +Target magnetic moments and constraint flags are specified per atom in the `ATOMIC_POSITIONS` section of the STRU file, using the `mag` (or `magmom`), `sc`, `lambda`, `angle1`, and `angle2` keywords after the atomic coordinates. + +#### Collinear (nspin=2) + +For collinear spin, only the z-component of the magnetic moment is constrained: + +``` +ATOMIC_POSITIONS +Direct + +Fe +0.0 +2 +0.00 0.00 0.00 mag 2.0 sc 1 +0.51 0.51 0.51 mag -2.0 sc 1 +``` + +- `mag 2.0`: target magnetic moment of 2.0 $\mu_B$ along z-axis +- `sc 1`: constrain the z-component (1 = constrained, 0 = unconstrained) + +#### Non-collinear (nspin=4), vector form + +For non-collinear spin, specify the magnetic moment as a vector (mx, my, mz): + +``` +ATOMIC_POSITIONS +Direct + +Fe +0.0 +2 +0.00 0.00 0.00 mag 2.0 0.0 0.0 sc 1 1 1 +0.51 0.51 0.51 mag 0.0 0.0 -2.0 sc 1 1 1 +``` + +- `mag 2.0 0.0 0.0`: target moment vector in Cartesian coordinates ($\mu_B$) +- `sc 1 1 1`: constrain x, y, z components respectively + +#### Non-collinear (nspin=4), angle form + +Alternatively, use `angle1` (polar angle $\theta$) and `angle2` (azimuthal angle $\phi$) in degrees to specify the direction: + +``` +0.00 0.00 0.00 mag 2.0 angle1 0 angle2 0 sc 1 1 1 +0.51 0.51 0.51 mag 2.0 angle1 180 angle2 0 sc 1 1 1 +``` + +The Cartesian components are computed as: +- $m_z = |\mathbf{m}| \cos\theta$ +- $m_x = |\mathbf{m}| \sin\theta \cos\phi$ +- $m_y = |\mathbf{m}| \sin\theta \sin\phi$ + +#### Providing initial Lagrange multipliers + +Initial lambda values (in eV/$\mu_B$) can be provided via the `lambda` keyword to accelerate convergence: + +``` +0.00 0.00 0.00 mag 2.0 lambda 0.01 0.0 0.0 sc 1 1 1 +``` + +A single value sets $\lambda_z$; three values set $\lambda_x$, $\lambda_y$, $\lambda_z$. + +#### Partial constraints + +Set `sc 0` for unconstrained components. For example, to constrain only the direction but not the magnitude (use with `sc_direction_only`): + +``` +0.00 0.00 0.00 mag 2.0 0.0 0.0 sc 1 1 0 +``` + +### DeltaSpin INPUT Parameters + +| Parameter | Type | Default | Description | +|-----------|------|---------|-------------| +| `sc_mag_switch` | Boolean | False | Enable DeltaSpin | +| `sc_thr` | Real | 1.0e-6 | Convergence criterion for lambda loop (RMS, in $\mu_B$) | +| `nsc` | Integer | 100 | Maximum number of lambda iterations | +| `nsc_min` | Integer | 2 | Minimum number of lambda iterations | +| `sc_scf_nmin` | Integer | 2 | Minimum outer SCF iterations before starting lambda loop | +| `alpha_trial` | Real | 0.01 | Initial trial step size for lambda (eV/$\mu_B^2$) | +| `sccut` | Real | 3.0 | Maximum step size for lambda (eV/$\mu_B$) | +| `sc_drop_thr` | Real | 1.0e-2 | Convergence ratio threshold for adaptive lambda loop | +| `sc_scf_thr` | Real | 1.0e-4 | Density error threshold for entering lambda loop | +| `sc_direction_only` | Boolean | False | Constrain only the direction, not the magnitude | +| `sc_lambda_strategy` | String | bfgs | Lambda update strategy (see below) | +| `decay_grad_switch` | Boolean | False | Enable gradient-based early exit | + +For full parameter details, see the [Spin-Constrained DFT](../input_files/input-main.md#spin-constrained-dft) section of the input keyword list. + +### Lambda Update Strategies + +The `sc_lambda_strategy` parameter controls how the Lagrange multipliers $\lambda$ are updated during the lambda loop: + +- **`bfgs`** (default): BFGS quasi-Newton method with line search. Robust and well-tested for both PW and LCAO. Uses `alpha_trial` and `sccut` to control step size. + +- **`linear_response`**: Linear response method (Scheme B). Estimates the magnetic susceptibility $\chi$ from the history of $(\lambda, M)$ pairs and performs a one-step Newton-like update: $\Delta\lambda = \beta (M_{\text{target}} - M) / \chi$, where $\beta$ is a mixing parameter. + +- **`augmented_lagrangian`**: Augmented Lagrangian method (Scheme C). Uses a penalty parameter $\mu$ that grows over iterations: $\lambda_{\text{new}} = \lambda + \mu (M - M_{\text{target}})$. The penalty increases until convergence is achieved. + +- **`hybrid_delayed`**: Hybrid delayed update (Scheme D). Two-phase approach: in the early phase (SCF not yet converged), lambda updates are gentle; in the late phase (SCF nearly converged), augmented Lagrangian updates are applied. + +### Direction-Only Mode + +When `sc_direction_only 1` is set, only the **direction** of the magnetic moment is constrained to match the target, while the magnitude is allowed to vary freely. This is useful for: + +- Studying magnetic anisotropy energy surfaces +- Cases where the moment magnitude is determined by the electronic structure +- Converging to the easy-axis direction without fixing the moment size + +In this mode, the lambda vector is projected to be perpendicular to the target moment direction at each iteration, ensuring it can only rotate the magnetization, not stretch it. + +### Combining DeltaSpin with DFT+U + +DeltaSpin can be combined with DFT+U for strongly correlated systems. When both `sc_mag_switch` and `dft_plus_u` are enabled: + +1. DFT+U occupation update runs first in each SCF iteration +2. DeltaSpin lambda loop runs after, constraining the magnetic moments +3. The DFT+U-corrected Hamiltonian is used by the lambda loop + +Example INPUT for PW DFT+U + DeltaSpin: + +``` +INPUT_PARAMETERS +calculation scf +basis_type pw +ecutwfc 50 +nspin 2 +dft_plus_u 1 +orbital_corr -1 2 +hubbard_u 0.0 4.0 +sc_mag_switch 1 +sc_thr 1.0e-6 +sc_scf_thr 1.0e-4 +sc_lambda_strategy bfgs +``` + +### Example: Collinear antiferromagnetic Fe + +INPUT file: + +``` +INPUT_PARAMETERS +calculation scf +basis_type pw +ecutwfc 50 +nspin 2 +sc_mag_switch 1 +sc_thr 1.0e-6 +``` + +STRU file: + +``` +ATOMIC_SPECIES +Fe 55.845 Fe.upf + +LATTICE_CONSTANT +8.190 + +LATTICE_VECTORS + 1.00 0.50 0.50 + 0.50 1.00 0.50 + 0.50 0.50 1.00 + +ATOMIC_POSITIONS +Direct + +Fe +0.0 +2 +0.00 0.00 0.00 mag 2.0 sc 1 +0.51 0.51 0.51 mag -2.0 sc 1 +``` + +### Example: Non-collinear constrained moments + +INPUT file: + +``` +INPUT_PARAMETERS +calculation scf +basis_type pw +ecutwfc 50 +nspin 4 +noncolin 1 +sc_mag_switch 1 +sc_direction_only 1 +sc_lambda_strategy bfgs +``` + +STRU file: + +``` +ATOMIC_POSITIONS +Direct + +Fe +0.0 +2 +0.00 0.00 0.00 mag 2.0 0.0 0.0 sc 1 1 0 +0.51 0.51 0.51 mag 0.0 0.0 2.0 sc 1 1 0 +``` + ## Noncollinear Spin Polarized Calculations The spin non-collinear polarization calculation corresponds to setting **"noncolin 1"**, in which case the coupling between spin up and spin down will be taken into account. In this case, nspin is automatically set to 4, which is usually not required to be specified manually. @@ -48,6 +266,53 @@ Note: different settings for "noncolin" and "lspinorb" correspond to different c **Special case**: `noncolin=0, lspinorb=1` is commonly used for non-magnetic materials with SOC effects (e.g., topological insulators, semiconductors with spin-orbit splitting). In this case, the magnetization is NOT automatically set, implying no magnetic moments in the system. +### Choosing gga_grad for Noncollinear Calculations + +When performing noncollinear spin calculations (`nspin=4`), the `gga_grad` parameter controls how the gradient of the magnetization is computed in GGA exchange-correlation functionals. This is critical for obtaining correct magnetic anisotropy energies and spin-orbit coupling effects. + +#### Available Options + +| gga_grad | Gradient Method | Description | When to Use | +|----------|----------------|-------------|-------------| +| 0 | Standard | Conventional GGA gradient (collinear approximation) | Collinear calculations only (`nspin<4`) | +| 1 | $\nabla\|\mathbf{m}\|$ | Gradient of magnetization magnitude | Noncollinear with SOC, when only magnitude variation matters | +| 2 | $\nabla\hat{\mathbf{m}}$ | Gradient of magnetization direction (default for `nspin=4`) | **Recommended** for most noncollinear calculations with SOC | +| 3 | Built-in SF | Special functional for noncollinear SOC (requires `v_xc_ncgga_sf_builtin`) | Advanced noncollinear SOC calculations | + +#### Recommendations + +**For noncollinear calculations with SOC (`lspinorb=1`)**: +- **Use `gga_grad=2`** (default when `nspin=4`): This computes the gradient of the magnetization direction unit vector $\hat{\mathbf{m}}$, which is physically correct for spin-orbit coupling effects. +- This is automatically set when you specify `nspin=4` or `lspinorb=1`, but you can explicitly set it to be clear. + +**For noncollinear calculations without SOC (`noncolin=1`, `lspinorb=0`)**: +- **Use `gga_grad=2`** for spin spiral calculations or frustrated magnetic systems +- The gradient of the magnetization direction is still important even without SOC + +**For collinear calculations (`nspin=1` or `nspin=2`)**: +- **Use `gga_grad=0`** (default): Standard collinear GGA gradient +- Do NOT use `gga_grad=1`, `2`, or `3` for collinear calculations + +#### Example INPUT for Noncollinear SOC + +``` +INPUT_PARAMETERS +calculation scf +basis_type pw +ecutwfc 50 +nspin 4 +noncolin 1 +lspinorb 1 +gga_grad 2 +``` + +#### Important Notes + +1. **Automatic setting**: When `nspin=4`, ABACUS automatically sets `gga_grad=2` if not explicitly specified +2. **Physical correctness**: Using the wrong `gga_grad` value can lead to incorrect magnetic anisotropy energies and spin textures +3. **Compatibility**: `gga_grad=3` requires the built-in special functional and may not be available for all XC functionals +4. **Performance**: `gga_grad=2` has similar computational cost to `gga_grad=0` for most systems + ## For the continuation job - Continuation job for "nspin 1" need file "SPIN1_CHG.cube" which is generated by setting "out_chg=1" in task before. By setting "init_chg file" in new job's INPUT file, charge density will start from file but not atomic. - Continuation job for "nspin 2" need files "SPIN1_CHG.cube" and "SPIN2_CHG.cube" which are generated by "out_chg 1" with "nspin 2", and refer to spin-up and spin-down charge densities respectively. It should be note that reading "SPIN1_CHG.cube" only for the continuation target magnetic moment job is not supported now. diff --git a/docs/parameters.yaml b/docs/parameters.yaml index 005fad3076..31f216e9dc 100644 --- a/docs/parameters.yaml +++ b/docs/parameters.yaml @@ -30,7 +30,7 @@ parameters: * md: perform molecular dynamics simulations * get_pchg: obtain partial (band-decomposed) charge densities (for LCAO basis only). See out_pchg for more information * get_wf: obtain real space wave functions (for LCAO basis only). See out_wfc_norm and out_wfc_re_im for more information - * get_s: obtain the overlap matrix formed by localized orbitals (for LCAO basis with multiple k points). The file name is OUT.${suffix}/sr_nao.csr, with the same file format as generated by out_hsr 1 + * get_s: obtain the overlap matrix formed by localized orbitals (for LCAO basis with multiple k points). the file name is SR.csr with file format being the same as that generated by out_mat_hs2 * gen_bessel: generates projectors, i.e., a series of Bessel functions, for the DeePKS method (for LCAO basis only) * gen_opt_abfs: generate opt-ABFs as discussed in this article * test_memory: obtain a rough estimation of memory consumption for the calculation @@ -63,7 +63,7 @@ parameters: Takes value 1, 0 or -1. * -1: No symmetry will be considered. It is recommended to set -1 for non-colinear + soc calculations, where time reversal symmetry is broken sometimes. * 0: Only time reversal symmetry would be considered in symmetry operations, which implied k point and -k point would be treated as a single k point with twice the weight. - * 1: Symmetry analysis will be performed to determine the type of Bravais lattice and associated symmetry operations (point groups, space groups, primitive cells, and irreducible k-points). For a magnetic system, the symmetry of the initial magnetic structure will be analyzed and preserved. + * 1: Symmetry analysis will be performed to determine the type of Bravais lattice and associated symmetry operations. (point groups, space groups, primitive cells, and irreducible k-points) [NOTE] When symmetry is enabled (value 1), k-points are reduced to the irreducible Brillouin zone (IBZ). For explicit k-point lists with custom weights (see KPT file), the custom weights are preserved during symmetry reduction. For Monkhorst-Pack grids, uniform weights are used. default_value: default @@ -533,9 +533,9 @@ parameters: For plane-wave basis, * cg: The conjugate-gradient (CG) method. + * bpcg: The BPCG method, which is a block-parallel Conjugate Gradient (CG) method, typically exhibits higher acceleration in a GPU environment. * dav: The Davidson algorithm. * dav_subspace: The Davidson algorithm without orthogonalization operation, this method is the most recommended for efficiency. `pw_diag_ndim` can be set to 2 for this method. - * bpcg: The BPCG method, which is a block-parallel Conjugate Gradient (CG) method, typically exhibits higher acceleration in a GPU environment. The BPCG method is currently under testing and is not recommended for use. For numerical atomic orbitals basis, @@ -621,20 +621,20 @@ parameters: category: Electronic structure type: Integer followed by Real values description: | - Customized parameterization of the exchange part of an XC functional. The first value should be the Libxc ID of the original functional, followed by the complete list of external parameters required by the linked Libxc version. If unset, Libxc's own default parameters are used. For functional IDs and parameter definitions, refer to the Libxc documentation and source code. + Customized parameterization on the exchange part of XC functional. The first value should be the LibXC ID of the original functional, and latter values are external parameters. Default values are those of Perdew-Burke-Ernzerhof (PBE) functional. For more information on LibXC ID of functionals, please refer to LibXC. For parameters of functionals of interest, please refer to the source code of LibXC, such as PBE functional interface in LibXC: gga_x_pbe.c. [NOTE] Solely setting this keyword will take no effect on XC functionals. One should also set dft_functional to the corresponding functional to apply the customized parameterization. Presently this feature can only support parameterization on one exchange functional. - default_value: "" + default_value: 101 0.8040 0.2195149727645171 unit: "" availability: "" - name: xc_corr_ext category: Electronic structure type: Integer followed by Real values description: | - Customized parameterization of the correlation part of an XC functional. The first value should be the Libxc ID of the original functional, followed by the complete list of external parameters required by the linked Libxc version. If unset, Libxc's own default parameters are used. For functional IDs and parameter definitions, refer to the Libxc documentation and source code. + Customized parameterization on the correlation part of XC functional. The first value should be the LibXC ID of the original functional, and latter values are external parameters. Default values are those of Perdew-Burke-Ernzerhof (PBE) functional. For more information on LibXC ID of functionals, please refer to LibXC. For parameters of functionals of interest, please refer to the source code of LibXC, such as PBE functional interface in LibXC: gga_c_pbe.c. [NOTE] Solely setting this keyword will take no effect on XC functionals. One should also set dft_functional to the corresponding functional to apply the customized parameterization. Presently this feature can only support parameterization on one correlation functional. - default_value: "" + default_value: 130 0.06672455060314922 0.031090690869654895034 1.0 unit: "" availability: "" - name: pseudo_rcut @@ -2953,7 +2953,7 @@ parameters: category: Output information type: "Boolean \\[Integer\\](optional)" description: | - Whether to output the density matrix with Bravias lattice vector R index into files in the folder OUT.${suffix}. The files are named as dmr{s}{spin index}{g}{geometry index}{_nao} + {".csr"}. Here, 's' refers to spin, where s1 means spin up channel while s2 means spin down channel, and the sparse matrix format 'csr' is mentioned in out_hsr. Finally, if out_app_flag is set to false, the file name contains the optional 'g' index for each ionic step that may have different geometries, and if out_app_flag is set to true, the density matrix with respect to Bravias lattice vector R accumulates during ionic steps: + Whether to output the density matrix with Bravias lattice vector R index into files in the folder OUT.${suffix}. The files are named as dmr{s}{spin index}{g}{geometry index}{_nao} + {".csr"}. Here, 's' refers to spin, where s1 means spin up channel while s2 means spin down channel, and the sparse matrix format 'csr' is mentioned in out_mat_hs2. Finally, if out_app_flag is set to false, the file name contains the optional 'g' index for each ionic step that may have different geometries, and if out_app_flag is set to true, the density matrix with respect to Bravias lattice vector R accumulates during ionic steps: * nspin = 1: dmrs1_nao.csr; * nspin = 2: dmrs1_nao.csr and dmrs2_nao.csr for the two spin channels. @@ -3065,61 +3065,34 @@ parameters: default_value: ie unit: "" availability: "" - - name: out_hsk - category: Output information - type: "Integer \\[Integer\\](optional)" - description: | - Output the upper triangular part of the Hamiltonian and overlap matrices in reciprocal space for each k-point into files in the directory OUT.${suffix}. The first integer selects the format: - * 0: disabled; - * 1: text output; the optional second integer controls precision and defaults to 8; - * 2: reserved for binary output, which is not implemented yet; - * 3: NPZ output, which is not implemented for H(k)/S(k). - - The output is also controlled by out_freq_ion and out_app_flag. For more information, refer to hs_matrix.md. - * Gamma-only, nspin = 1: hk_nao.txt for the Hamiltonian matrix and sk_nao.txt for the overlap matrix. - * Gamma-only, nspin = 2: hks1_nao.txt and hks2_nao.txt for the two spin channels of the Hamiltonian matrix, and sk_nao.txt for the overlap matrix. Only one overlap matrix is written because it is identical for both spin channels. - * Gamma-only, nspin = 4: not available with the gamma-only algorithm. - * Multi-k, nspin = 1: hk1_nao.txt for the Hamiltonian matrix and sk1_nao.txt for the overlap matrix at the first k-point. - * Multi-k, nspin = 2: hk1s1_nao.txt and hk1s2_nao.txt for the two spin channels of the Hamiltonian matrix, and sk1_nao.txt for the overlap matrix at the first k-point. Only one overlap matrix is written because it is identical for both spin channels. - * Multi-k, nspin = 4: hk1s4_nao.txt for the spinor Hamiltonian matrix and sk1_nao.txt for the spinor overlap matrix at the first k-point. - When out_app_flag is false, g followed by the one-based ionic-step index is inserted before _nao, for example hk1s1g1_nao.txt. - - [NOTE] In the 3.10-LTS version, the file names are data-0-H and data-0-S, etc. - default_value: 0 8 - unit: Ry - availability: Numerical atomic orbital basis - name: out_mat_hs category: Output information type: "Boolean \\[Integer\\](optional)" description: | - Legacy alias for out_hsk 1, which outputs Hamiltonian and overlap matrices in reciprocal space for each k-point. The optional second integer controls text precision. If both out_hsk and out_mat_hs are present, out_hsk takes precedence. - default_value: False 8 - unit: Ry - availability: Numerical atomic orbital basis - - name: out_hsr - category: Output information - type: "Integer \\[Integer\\](optional)" - description: | - Output Hamiltonian and overlap matrices in real space, indexed by the Bravais lattice vector R, in the directory OUT.${suffix}. The first integer selects the format: - * 0: disabled; - * 1: text CSR output; the optional second integer controls precision and defaults to 8; - * 2: reserved for binary output, which is not implemented yet; - * 3: NPZ output using hrs1_nao.npz, hrs2_nao.npz when needed, and sr_nao.npz. - - For multi-k calculations, the output contains the individual real-space blocks stored for the Bravais lattice vectors R. For gamma-only calculations, the internal real-space contributions are folded into a single R = (0, 0, 0) block. This folded result cannot recover the original R-resolved contributions or interpolate arbitrary k points. Terms added only while constructing H(k) are not guaranteed to be present. + Whether to print the upper triangular part of the Hamiltonian matrices and overlap matrices for each k-point into files in the directory OUT.${suffix}. The second number controls precision. For more information, please refer to hs_matrix.md. Also controled by out_freq_ion and out_app_flag. + * For gamma only case: + * nspin = 1: hks1_nao.txt for the Hamiltonian matrix and sks1_nao.txt for the overlap matrix; + * nspin = 2: hks1_nao.txt and hks2_nao.txt for the Hamiltonian matrix and sks1_nao.txt for the overlap matrix. Note that the code will not output sks2_nao.txt because it is the same as sks1_nao.txt; + * nspin = 4: hks12_nao.txt for the Hamiltonian matrix and sks12_nao.txt for the overlap matrix. + * For multi-k points case: + * nspin = 1: hks1k1_nao.txt for the Hamiltonian matrix at the 1st k-point, and sks1k1_nao.txt for the overlap matrix for the 1st k-point, ...; + * nspin = 2: hks1k1_nao.txt and hks2k1_nao.txt for the two spin channels of the Hamiltonian matrix at the 1st k-point, and sks1k1_nao.txt for the overlap matrix for the 1st k-point. Note that the code will not output sks2k1_nao.txt because it is the same as sks1k1_nao.txt, ...; + * nspin = 4: hks12k1_nao.txt for the Hamiltonian matrix at the 1st k-point, and sks12k1_nao.txt for the overlap matrix for the 1st k-point, ...; - [NOTE] In the 3.10-LTS version, the file names are data-HR-sparse_SPIN0.csr and data-SR-sparse_SPIN0.csr, etc. - default_value: 0 8 + [NOTE] In the 3.10-LTS version, the file names are data-0-H and data-0-S, etc. + default_value: False 8 unit: Ry availability: Numerical atomic orbital basis - name: out_mat_hs2 category: Output information type: "Boolean \\[Integer\\](optional)" description: | - Legacy alias for out_hsr 1, which outputs Hamiltonian and overlap matrices in real space indexed by the Bravais lattice vector R. The optional second integer controls text precision. If both out_hsr and out_mat_hs2 are present, out_hsr takes precedence. - default_value: False 8 + Whether to print files containing the Hamiltonian matrix and overlap matrix into files in the directory OUT.${suffix}. For more information, please refer to hs_matrix.md. + + [NOTE] In the 3.10-LTS version, the file names are data-HR-sparse_SPIN0.csr and data-SR-sparse_SPIN0.csr, etc. + default_value: "False [8]" unit: Ry - availability: Numerical atomic orbital basis + availability: Numerical atomic orbital basis (not gamma-only algorithm) - name: out_mat_tk category: Output information type: "Boolean \\[Integer\\](optional)" @@ -3144,7 +3117,7 @@ parameters: category: Output information type: "Boolean \\[Integer\\](optional)" description: | - Generate files containing the kinetic energy matrix. The optional second parameter controls text output precision. The format will be the same as the Hamiltonian matrix and overlap matrix as mentioned in out_hsr. The name of the files will be trs1_nao.csr and so on. Also controled by out_freq_ion and out_app_flag. + Generate files containing the kinetic energy matrix. The optional second parameter controls text output precision. The format will be the same as the Hamiltonian matrix and overlap matrix as mentioned in out_mat_hs2. The name of the files will be trs1_nao.csr and so on. Also controled by out_freq_ion and out_app_flag. [NOTE] In the 3.10-LTS version, the file name is data-TR-sparse_SPIN0.csr. default_value: False 8 @@ -3154,7 +3127,7 @@ parameters: category: Output information type: Integer description: | - Whether to print files containing the derivatives of the Hamiltonian matrix. The format will be the same as the Hamiltonian matrix and overlap matrix as mentioned in out_hsr. The name of the files will be dhrxs1_nao.csr, dhrys1_nao.csr, dhrzs1_nao.csr and so on. Also controled by out_freq_ion and out_app_flag. + Whether to print files containing the derivatives of the Hamiltonian matrix. The format will be the same as the Hamiltonian matrix and overlap matrix as mentioned in out_mat_hs2. The name of the files will be dhrxs1_nao.csr, dhrys1_nao.csr, dhrzs1_nao.csr and so on. Also controled by out_freq_ion and out_app_flag. Format: [precision] [iat1 iat2 ...]. The first value (0/1) enables/disables output. The second optional value sets the output precision (default: 8). Starting from the third value, 1-based atom indices can be listed to restrict output to derivatives with respect to those specific atoms only; if no atom indices are given, all atoms are written. @@ -3228,7 +3201,7 @@ parameters: description: | Whether to print files containing the kinetic energy matrix T(R) in CSR format. - See out_hsr for format details. + See out_mat_hs2 for format details. default_value: 0 8 unit: Ry availability: "" @@ -3238,7 +3211,7 @@ parameters: description: | Whether to print files containing the nonlocal pseudopotential matrix Vnl(R) in CSR format. - See out_hsr for format details. + See out_mat_hs2 for format details. default_value: 0 8 unit: Ry availability: "" @@ -3248,7 +3221,7 @@ parameters: description: | Whether to print files containing the local pseudopotential matrix Vl(R) in CSR format. - See out_hsr for format details. + See out_mat_hs2 for format details. default_value: 0 8 unit: Ry availability: "" @@ -3258,7 +3231,7 @@ parameters: description: | Whether to print files containing the Hartree matrix Vh(R) in CSR format. - See out_hsr for format details. + See out_mat_hs2 for format details. default_value: 0 8 unit: Ry availability: "" @@ -3268,7 +3241,7 @@ parameters: description: | Whether to print files containing the XC matrix Vxc(R) in CSR format. - See out_hsr for format details. + See out_mat_hs2 for format details. default_value: 0 8 unit: Ry availability: "" @@ -3278,7 +3251,7 @@ parameters: description: | Whether to print files containing the exact-exchange matrix Vexx(R) in CSR format. - See out_hsr for format details. + See out_mat_hs2 for format details. default_value: 0 8 unit: Ry availability: "" @@ -3349,7 +3322,7 @@ parameters: category: Output information type: Boolean description: | - Whether to print Hamiltonian matrices H(R) in NPZ format as hrs1_nao.npz and, for nspin = 2, hrs2_nao.npz. This feature does not work for gamma-only calculations. + Whether to print Hamiltonian matrices H(R) in npz format. This feature does not work for gamma-only calculations. default_value: "False" unit: Ry availability: Numerical atomic orbital basis (not gamma-only algorithm) @@ -3357,10 +3330,10 @@ parameters: category: Output information type: Boolean description: | - Legacy alias for out_hsr 3, writing hrs1_nao.npz, hrs2_nao.npz when needed, and sr_nao.npz. If both out_hsr and out_hsr_npz are present, out_hsr takes precedence. Gamma-only calculations write the folded R = (0, 0, 0) representation. + Whether to print Hamiltonian matrices H(R) and overlap matrix S(R) in npz format. This feature does not work for gamma-only calculations. default_value: "False" unit: Ry - availability: Numerical atomic orbital basis + availability: Numerical atomic orbital basis (not gamma-only algorithm) - name: out_dm_npz category: Output information type: Boolean @@ -3381,7 +3354,7 @@ parameters: category: Output information type: Boolean description: | - Whether to output r(R), H(R), S(R), T(R), dH(R), dS(R), and wfc matrices in an append manner during molecular dynamics calculations. Check input parameters out_mat_r, out_hsr, out_mat_t, out_mat_dh, out_hsk and out_wfc_lcao for more information. + Whether to output r(R), H(R), S(R), T(R), dH(R), dS(R), and wfc matrices in an append manner during molecular dynamics calculations. Check input parameters out_mat_r, out_mat_hs2, out_mat_t, out_mat_dh, out_mat_hs and out_wfc_lcao for more information. default_value: "true" unit: "" availability: Numerical atomic orbital basis (not gamma-only algorithm) @@ -3392,7 +3365,7 @@ parameters: Controls the length of decimal part of output data, such as charge density, Hamiltonian matrix, Overlap matrix and so on. default_value: "8" unit: "" - availability: out_hsk 1 case presently. + availability: out_mat_hs 1 case presently. - name: out_element_info category: Output information type: Boolean @@ -4500,8 +4473,45 @@ parameters: category: Spin-Constrained DFT type: Real description: | - Density error threshold for inner loop of spin-constrained SCF - default_value: "1.0e-4" + When the charge density error drho falls below sc_scf_thr, the DeltaSpin lambda optimization loop is activated. Should be 10-100x larger than scf_thr. Only used when sc_scf_thr_mode="threshold". For other activation modes, use sc_scf_thr_mode="immediate" or sc_scf_thr_mode="off". + default_value: "1.0e-3" + unit: "" + availability: sc_mag_switch is true + - name: sc_scf_thr_mode + category: Spin-Constrained DFT + type: String + description: | + Controls when the DeltaSpin lambda loop is activated. + * threshold (default): activate when drho < sc_scf_thr. The lambda loop starts once the charge density is reasonably stable. + * immediate: activate from the first iteration with valid wavefunctions (iter>=2). Used for PW basis where the first iteration cannot compute initial magnetic moments. Replaces the old convention of setting sc_scf_thr=10.0. + * off: never activate the lambda loop. Lambda values are loaded from the STRU file and used as constant constraints without optimization. Replaces the old convention of setting sc_scf_thr=1e-10. + default_value: "threshold" + unit: "" + availability: sc_mag_switch is true + - name: sc_direction_only + category: Spin-Constrained DFT + type: Boolean + description: | + When true, only the direction of the magnetic moment is constrained to the target direction, while the magnitude is allowed to vary freely. This is useful for studying magnetic anisotropy or when the magnitude of the moment is determined by the electronic structure rather than an external constraint. When false (default), both the direction and magnitude of the magnetic moment are constrained to the target values. + default_value: "False" + unit: "" + availability: sc_mag_switch is true + - name: sc_dir_phase1_steps + category: Spin-Constrained DFT + type: Integer + description: | + Number of SCF iterations for Phase 1 (magnitude constraint) in the direction_only two-phase strategy for collinear (nspin=2) calculations. During Phase 1, direction_only projection is temporarily disabled so BFGS can constrain the magnetic moment magnitude. After Phase 1, lambda decays and the system relaxes naturally. Minimum: 2. + default_value: "5" + unit: "" + availability: sc_mag_switch is true and sc_direction_only is true and nspin=2 + - name: sc_lambda_strategy + category: Spin-Constrained DFT + type: String + description: | + Lambda update strategy for spin-constrained DFT. Available options are: + * bfgs: BFGS quasi-Newton method (default, robust and well-tested) + * linear_scan: linear sweep of lambda for testing magnetic moment response + default_value: "bfgs" unit: "" availability: sc_mag_switch is true - name: sc_direction_only diff --git a/source/source_hamilt/module_xc/CMakeLists.txt b/source/source_hamilt/module_xc/CMakeLists.txt index 5f310d09d5..4a016ed04d 100644 --- a/source/source_hamilt/module_xc/CMakeLists.txt +++ b/source/source_hamilt/module_xc/CMakeLists.txt @@ -2,6 +2,7 @@ add_library( xc_ OBJECT xc_functional.cpp + xc_functional_gga_noncol_sf_builtin.cpp xc_pot.cpp xc_grad.cpp xc_lda_wrap.cpp diff --git a/source/source_hamilt/module_xc/libxc_abacus.h b/source/source_hamilt/module_xc/libxc_abacus.h index bc584b7518..23198d8cc5 100644 --- a/source/source_hamilt/module_xc/libxc_abacus.h +++ b/source/source_hamilt/module_xc/libxc_abacus.h @@ -24,30 +24,14 @@ namespace XC_Functional_Libxc //------------------- // sets functional type, which allows combination of LIBXC keyword connected by "+" - // for example, "XC_LDA_X+XC_LDA_C_PZ" + // for example: "XC_LDA_X+XC_LDA_C_PZ" extern std::pair> set_xc_type_libxc(const std::string& xc_func_in); - /** - * @brief instantiate the XC functional by its ID, and set the external parameters if provided. - * - * @param func_id libxc ID of functional, see https://libxc.gitlab.io/functionals/ for details - * @param xc_polarized 0: unpolarized, 1: spin-polarized - * @return std::vector - * - * @note the functionality of this method is extended by supporting the user-defined - * external parameters of xc. However, there are several functionals' external - * parameters are pre-defined in the code, which herein we call those are - * "in-built" parameters. If the same functional ID is found in both in-built - * and external parameters, the external parameters will overwrite the in-built ones. - * The external parameters can be passed here by keywords xc_exch_ext and - * xc_corr_ext in the input file. The expected format would be an XC ID - * followed by a list of parameters. - */ extern std::vector init_func( const std::vector &func_id, const int xc_polarized, - const double hybrid_alpha, - const double hse_omega); + const double hybrid_alpha = 0.0, + const double hse_omega = 0.0); extern void finish_func(std::vector &funcs); @@ -65,9 +49,9 @@ namespace XC_Functional_Libxc const int nspin, const bool domag, const bool domag_z, - const std::map* scaling_factor, - const double hybrid_alpha, - const double hse_omega); + const std::map* scaling_factor = nullptr, + const double hybrid_alpha = 0.0, + const double hse_omega = 0.0); // for mGGA functional extern std::tuple v_xc_meta( @@ -77,8 +61,8 @@ namespace XC_Functional_Libxc const double tpiba, const Charge* const chr, const int nspin, - const double hybrid_alpha, - const double hse_omega); + const double hybrid_alpha = 0.0, + const double hse_omega = 0.0); //------------------- @@ -157,6 +141,34 @@ namespace XC_Functional_Libxc const std::vector &amag, const ModuleBase::matrix &v); + extern std::vector compute_mag_part_nspin4( + const std::size_t nrxx, + const Charge* const chr); + + extern std::vector>> cal_gdr_sf( + const int nspin, + const std::size_t nrxx, + const std::vector &rho, + const std::vector &mag_part, + const double tpiba, + const Charge* const chr); + + extern std::vector> cal_dh_sf( + const int nspin, + const std::size_t nrxx, + const std::vector &sgn, + const std::vector>> &gdr, + const std::vector &vsigma, + const std::vector &mag_part, + const double tpiba, + const Charge* const chr); + + extern ModuleBase::matrix convert_v_nspin4_sf( + const std::size_t nrxx, + const Charge* const chr, + const std::vector &mag_part, + const ModuleBase::matrix &v); + //------------------- // libxc_lda_wrap.cpp @@ -169,8 +181,8 @@ namespace XC_Functional_Libxc double &exc, double &vxcup, double &vxcdw, - const double hybrid_alpha, - const double hse_omega); + const double hybrid_alpha = 0.0, + const double hse_omega = 0.0); //------------------- @@ -185,8 +197,8 @@ namespace XC_Functional_Libxc double &sxc, double &v1xc, double &v2xc, - const double hybrid_alpha, - const double hse_omega); + const double hybrid_alpha = 0.0, + const double hse_omega = 0.0); // the entire GGA functional, for nspin=2 case extern void gcxc_spin_libxc( @@ -201,8 +213,8 @@ namespace XC_Functional_Libxc double &v2xcup, double &v2xcdw, double &v2xcud, - const double hybrid_alpha, - const double hse_omega); + const double hybrid_alpha = 0.0, + const double hse_omega = 0.0); //------------------- @@ -219,8 +231,8 @@ namespace XC_Functional_Libxc double &v1xc, double &v2xc, double &v3xc, - const double &hybrid_alpha, - const double &hse_omega); + const double &hybrid_alpha = 0.0, + const double &hse_omega = 0.0); extern void tau_xc_spin( const std::vector &func_id, @@ -238,8 +250,8 @@ namespace XC_Functional_Libxc double &v2xcud, double &v3xcup, double &v3xcdw, - const double &hybrid_alpha, - const double &hse_omega); + const double &hybrid_alpha = 0.0, + const double &hse_omega = 0.0); } // namespace XC_Functional_Libxc diff --git a/source/source_hamilt/module_xc/libxc_pot.cpp b/source/source_hamilt/module_xc/libxc_pot.cpp index 6f605bbd48..a6fc73d9ed 100644 --- a/source/source_hamilt/module_xc/libxc_pot.cpp +++ b/source/source_hamilt/module_xc/libxc_pot.cpp @@ -66,7 +66,8 @@ std::tuple XC_Functional_Libxc::v_xc_libxc( / // converting rho std::vector rho; std::vector amag; - if(1==nspin || 2==nspin_in) + std::vector mag_part; + if(1==nspin || 2==PARAM.inp.nspin) { rho = XC_Functional_Libxc::convert_rho(nspin, nrxx, chr); } @@ -75,13 +76,24 @@ std::tuple XC_Functional_Libxc::v_xc_libxc( / std::tuple,std::vector> rho_amag = XC_Functional_Libxc::convert_rho_amag_nspin4(nspin, nrxx, chr); rho = std::get<0>(std::move(rho_amag)); amag = std::get<1>(std::move(rho_amag)); + if((PARAM.inp.gga_grad == 2 || PARAM.inp.gga_grad == 3) && (PARAM.globalv.domag || PARAM.globalv.domag_z)) + { + mag_part = XC_Functional_Libxc::compute_mag_part_nspin4(nrxx, chr); + } } std::vector>> gdr; std::vector sigma; if(is_gga) { - gdr = XC_Functional_Libxc::cal_gdr(nspin, nrxx, rho, tpiba, chr); + if(PARAM.inp.nspin==4 && (PARAM.inp.gga_grad == 2 || PARAM.inp.gga_grad == 3) && (PARAM.globalv.domag || PARAM.globalv.domag_z)) + { + gdr = XC_Functional_Libxc::cal_gdr_sf(nspin, nrxx, rho, mag_part, tpiba, chr); + } + else + { + gdr = XC_Functional_Libxc::cal_gdr(nspin, nrxx, rho, tpiba, chr); + } sigma = XC_Functional_Libxc::convert_sigma(gdr); } @@ -179,7 +191,14 @@ std::tuple XC_Functional_Libxc::v_xc_libxc( / if(4==nspin_in) { - v = XC_Functional_Libxc::convert_v_nspin4(nrxx, chr, amag, v); + if((PARAM.inp.gga_grad == 2 || PARAM.inp.gga_grad == 3) && (PARAM.globalv.domag || PARAM.globalv.domag_z)) + { + v = XC_Functional_Libxc::convert_v_nspin4_sf(nrxx, chr, mag_part, v); + } + else + { + v = XC_Functional_Libxc::convert_v_nspin4(nrxx, chr, amag, v); + } } //------------------------------------------------- diff --git a/source/source_hamilt/module_xc/libxc_setup.cpp b/source/source_hamilt/module_xc/libxc_setup.cpp index 9a2a4c0e17..3b29edf727 100644 --- a/source/source_hamilt/module_xc/libxc_setup.cpp +++ b/source/source_hamilt/module_xc/libxc_setup.cpp @@ -14,6 +14,7 @@ #include #include #include +#include bool not_supported_xc_with_laplacian(const std::string& xc_func_in) { @@ -199,7 +200,7 @@ const std::vector in_built_xc_func_ext_params(const int id, #ifdef __EXX // hybrid functionals case XC_HYB_GGA_XC_PBEH: - return {hybrid_alpha}; + return {hybrid_alpha, hse_omega, hse_omega}; case XC_HYB_GGA_XC_HSE06: return {hybrid_alpha, hse_omega, hse_omega}; // short-range of B88_X @@ -261,19 +262,20 @@ const std::vector in_built_xc_func_ext_params(const int id, const std::vector external_xc_func_ext_params(const int id) { - const auto& exch_ext = PARAM.inp.xc_exch_ext; - if (!exch_ext.empty() && static_cast(exch_ext.front()) == id) - { - return {exch_ext.begin() + 1, exch_ext.end()}; - } - - const auto& corr_ext = PARAM.inp.xc_corr_ext; - if (!corr_ext.empty() && static_cast(corr_ext.front()) == id) - { - return {corr_ext.begin() + 1, corr_ext.end()}; - } - - return {}; + const std::map> mymap = { + { + PARAM.inp.xc_exch_ext[0], + std::vector(PARAM.inp.xc_exch_ext.begin()+1, + PARAM.inp.xc_exch_ext.end()) + }, + { + PARAM.inp.xc_corr_ext[0], + std::vector(PARAM.inp.xc_corr_ext.begin()+1, + PARAM.inp.xc_corr_ext.end()) + } + }; + auto it = mymap.find(id); + return (it != mymap.end()) ? it->second : std::vector{}; } std::vector @@ -288,32 +290,33 @@ XC_Functional_Libxc::init_func(const std::vector &func_id, funcs.push_back({}); // create placeholder xc_func_init(&funcs.back(), id, xc_polarized); // instantiate the XC term - // Search for external parameters. User-supplied parameters take precedence - // over ABACUS built-in overrides. - const std::vector in_built_ext_params - = in_built_xc_func_ext_params(id, hybrid_alpha, hse_omega); - const std::vector external_ext_params - = external_xc_func_ext_params(id); - const std::vector& requested_ext_params - = external_ext_params.empty() ? in_built_ext_params : external_ext_params; + // search for external parameters + const std::vector in_built_ext_params = in_built_xc_func_ext_params(id, hybrid_alpha, hse_omega); + const std::vector external_ext_params = external_xc_func_ext_params(id); + // for temporary use, I name their size as n1 and n2 + const int n1 = in_built_ext_params.size(); + const int n2 = external_ext_params.size(); - if (!requested_ext_params.empty()) - { - const xc_func_info_type* info = xc_func_get_info(&funcs.back()); - const int nref = xc_func_info_get_n_ext_params(info); +// #ifdef __DEBUG // will the following assertion cause performance issue? + // assert the number of parameters should be either zero or the value from + // libxc function xc_func_info_get_n_ext_params, this is to avoid the undefined + // behavior of illegal memory access + const xc_func_info_type* info = xc_func_get_info(&funcs.back()); + const int nref = xc_func_info_get_n_ext_params(info); + assert ((n1 == 0) || (n1 == nref) || (n2 == 0) || (n2 == nref)); +// #endif - // xc_func_set_ext_params() reads exactly nref entries - if (requested_ext_params.size() != static_cast(nref)) - { - ModuleBase::WARNING_QUIT( - "XC_Functional_Libxc::init_func", - "Invalid number of external parameters for Libxc functional id " - + std::to_string(id) + ": got " - + std::to_string(requested_ext_params.size()) - + ", expected " + std::to_string(nref) + "."); - } + // external overwrites in-built if the same functional id is found in both maps + const double* xc_func_ext_params = + (n2 > 0) ? external_ext_params.data() : + (n1 > 0) ? in_built_ext_params.data() : + nullptr; // nullptr if no external parameters are found - xc_func_set_ext_params(&funcs.back(), requested_ext_params.data()); + // if there are no external parameters, do nothing, otherwise we set + if(xc_func_ext_params != nullptr) + { + // set the external parameters + xc_func_set_ext_params(&funcs.back(), const_cast(xc_func_ext_params)); } } return funcs; diff --git a/source/source_hamilt/module_xc/libxc_tools.cpp b/source/source_hamilt/module_xc/libxc_tools.cpp index 916732fa8d..e3694044f0 100644 --- a/source/source_hamilt/module_xc/libxc_tools.cpp +++ b/source/source_hamilt/module_xc/libxc_tools.cpp @@ -51,6 +51,28 @@ XC_Functional_Libxc::convert_rho_amag_nspin4( return std::make_tuple(std::move(rho), std::move(amag)); } +std::vector XC_Functional_Libxc::compute_mag_part_nspin4( + const std::size_t nrxx, + const Charge* const chr) +{ + std::vector mag_part(3 * nrxx, 0.0); + #ifdef _OPENMP + #pragma omp parallel for schedule(static, 1024) + #endif + for (std::size_t ir = 0; ir < nrxx; ++ir) + { + double mx = chr->rho[1][ir], my = chr->rho[2][ir], mz = chr->rho[3][ir]; + double amag = std::sqrt(mx * mx + my * my + mz * mz); + if (amag > 1e-12) + { + mag_part[ir] = mx / amag; + mag_part[ir + nrxx] = my / amag; + mag_part[ir + 2 * nrxx] = mz / amag; + } + } + return mag_part; +} + // calculating grho std::vector>> XC_Functional_Libxc::cal_gdr( @@ -213,22 +235,48 @@ std::pair XC_Functional_Libxc::convert_vtxc_v( if(func.info->family == XC_FAMILY_GGA || func.info->family == XC_FAMILY_HYB_GGA) { - const std::vector> dh = XC_Functional_Libxc::cal_dh(nspin, nrxx, sgn, gdr, vsigma, tpiba, chr); - - double rvtxc = 0.0; - #ifdef _OPENMP - #pragma omp parallel for collapse(2) reduction(+:rvtxc) schedule(static, 256) - #endif - for( int is=0; is mag_part_tmp = XC_Functional_Libxc::compute_mag_part_nspin4(nrxx, chr); + const std::vector> dh = XC_Functional_Libxc::cal_dh_sf(nspin, nrxx, sgn, gdr, vsigma, mag_part_tmp, tpiba, chr); + + constexpr int nspin4 = 4; + double rvtxc = 0.0; + #ifdef _OPENMP + #pragma omp parallel for collapse(2) reduction(+:rvtxc) schedule(static, 256) + #endif + for (int is = 0; is < nspin4; ++is) { - rvtxc += dh[is][ir] * rho[ir*nspin+is]; - v(is,ir) -= dh[is][ir]; + for (std::size_t ir = 0; ir < nrxx; ++ir) + { + double rho_ir = 0.0; + if (is == 0) { rho_ir = rho[ir * nspin + 0]; } + else { rho_ir = rho[ir * nspin + 0] * mag_part_tmp[ir + (is - 1) * nrxx]; } + rvtxc += dh[is][ir] * rho_ir; + v(is, ir) -= dh[is][ir]; + } } + vtxc -= rvtxc; } + else + { + const std::vector> dh = XC_Functional_Libxc::cal_dh(nspin, nrxx, sgn, gdr, vsigma, tpiba, chr); + + double rvtxc = 0.0; + #ifdef _OPENMP + #pragma omp parallel for collapse(2) reduction(+:rvtxc) schedule(static, 256) + #endif + for( int is=0; isfamily == XC_FAMILY_GGA || func.info->family == XC_FAMILY_HYB_GGA)) return std::make_pair(vtxc, std::move(v)); @@ -317,4 +365,159 @@ ModuleBase::matrix XC_Functional_Libxc::convert_v_nspin4( return v_nspin4; } +std::vector>> XC_Functional_Libxc::cal_gdr_sf( + const int nspin, + const std::size_t nrxx, + const std::vector &rho, + const std::vector &mag_part, + const double tpiba, + const Charge* const chr) +{ + std::vector>> gdr(nspin); + std::vector rhor(nrxx); + std::vector> rhog(chr->rhopw->npw); + std::vector> gdr_tmp(nrxx); + + #ifdef _OPENMP + #pragma omp parallel for schedule(static, 1024) + #endif + for (std::size_t ir = 0; ir < nrxx; ++ir) + { + rhor[ir] = rho[ir * nspin + 0]; + } + chr->rhopw->real2recip(rhor.data(), rhog.data()); + gdr[0].resize(nrxx); + XC_Functional::grad_rho(rhog.data(), gdr[0].data(), chr->rhopw, tpiba); + + gdr[1].resize(nrxx); + #ifdef _OPENMP + #pragma omp parallel for schedule(static, 1024) + #endif + for (std::size_t ir = 0; ir < nrxx; ++ir) + { + gdr_tmp[ir] = gdr[0][ir]; + gdr[0][ir] = 0.5 * gdr_tmp[ir]; + gdr[1][ir] = 0.5 * gdr_tmp[ir]; + } + + for (int is = 1; is <= 3; ++is) + { + chr->rhopw->real2recip(chr->rho[is], rhog.data()); + XC_Functional::grad_rho(rhog.data(), gdr_tmp.data(), chr->rhopw, tpiba); + const double* mp = mag_part.data() + (is - 1) * nrxx; + #ifdef _OPENMP + #pragma omp parallel for schedule(static, 1024) + #endif + for (std::size_t ir = 0; ir < nrxx; ++ir) + { + const ModuleBase::Vector3 g = 0.5 * gdr_tmp[ir] * mp[ir]; + gdr[0][ir] += g; + gdr[1][ir] -= g; + } + } + + return gdr; +} + +std::vector> XC_Functional_Libxc::cal_dh_sf( + const int nspin, + const std::size_t nrxx, + const std::vector &sgn, + const std::vector>> &gdr, + const std::vector &vsigma, + const std::vector &mag_part, + const double tpiba, + const Charge* const chr) +{ + std::vector> h1(nrxx), h2(nrxx); + #ifdef _OPENMP + #pragma omp parallel for schedule(static, 1024) + #endif + for (std::size_t ir = 0; ir < nrxx; ++ir) + { + h1[ir] = 2.0 * (gdr[0][ir] * vsigma[ir*3 ] * sgn[ir*2 ] * 2.0 + + gdr[1][ir] * vsigma[ir*3+1] * sgn[ir*2] * sgn[ir*2+1]); + h2[ir] = 2.0 * (gdr[1][ir] * vsigma[ir*3+2] * sgn[ir*2+1] * 2.0 + + gdr[0][ir] * vsigma[ir*3+1] * sgn[ir*2] * sgn[ir*2+1]); + } + + std::vector> tmp_h(nrxx); + std::vector dh0(nrxx), dh_mu(nrxx); + + #ifdef _OPENMP + #pragma omp parallel for schedule(static, 1024) + #endif + for (std::size_t ir = 0; ir < nrxx; ++ir) + { + tmp_h[ir] = 0.5 * (h1[ir] + h2[ir]); + } + XC_Functional::grad_dot(tmp_h.data(), dh0.data(), chr->rhopw, tpiba); + + std::vector> dh_total(4, std::vector(nrxx, 0.0)); +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 1024) +#endif + for (std::size_t ir = 0; ir < nrxx; ++ir) + { + dh_total[0][ir] = dh0[ir]; + } + + for (int mu = 1; mu < 4; ++mu) + { + const double* mp_mu = mag_part.data() + (mu - 1) * nrxx; + #ifdef _OPENMP + #pragma omp parallel for schedule(static, 1024) + #endif + for (std::size_t ir = 0; ir < nrxx; ++ir) + { + tmp_h[ir] = 0.5 * (h1[ir] - h2[ir]) * mp_mu[ir]; + } + XC_Functional::grad_dot(tmp_h.data(), dh_mu.data(), chr->rhopw, tpiba); + #ifdef _OPENMP + #pragma omp parallel for schedule(static, 1024) + #endif + for (std::size_t ir = 0; ir < nrxx; ++ir) + { + dh_total[mu][ir] = dh_mu[ir]; + } + } + + return dh_total; +} + +ModuleBase::matrix XC_Functional_Libxc::convert_v_nspin4_sf( + const std::size_t nrxx, + const Charge* const chr, + const std::vector &mag_part, + const ModuleBase::matrix &v) +{ + assert(PARAM.inp.nspin==4); + constexpr double vanishing_charge = 1.0e-10; + ModuleBase::matrix v_nspin4(PARAM.inp.nspin, nrxx); + + for (std::size_t ir = 0; ir < nrxx; ++ir) + { + v_nspin4(0, ir) = 0.5 * (v(0, ir) + v(1, ir)); + } + + if (PARAM.globalv.domag || PARAM.globalv.domag_z) + { + for (std::size_t ir = 0; ir < nrxx; ++ir) + { + double amag = std::sqrt(std::pow(chr->rho[1][ir], 2) + + std::pow(chr->rho[2][ir], 2) + + std::pow(chr->rho[3][ir], 2)); + if (amag > vanishing_charge) + { + const double vs = 0.5 * (v(0, ir) - v(1, ir)); + for (int ipol = 1; ipol < PARAM.inp.nspin; ++ipol) + { + v_nspin4(ipol, ir) = vs * mag_part[ir + (ipol - 1) * nrxx]; + } + } + } + } + return v_nspin4; +} + #endif diff --git a/source/source_hamilt/module_xc/test/CMakeLists.txt b/source/source_hamilt/module_xc/test/CMakeLists.txt index 34634eae02..e014d905cc 100644 --- a/source/source_hamilt/module_xc/test/CMakeLists.txt +++ b/source/source_hamilt/module_xc/test/CMakeLists.txt @@ -71,6 +71,7 @@ AddTest( ../xc_gga_corr.cpp ../xc_lda_corr.cpp ../xc_gga_exch.cpp ../xc_lda_exch.cpp ../xc_hcth.cpp ../xc_pot.cpp + ../xc_functional_gga_noncol_sf_builtin.cpp ../libxc_pot.cpp ../libxc_tools.cpp ../../../source_base/module_external/blas_connector_base.cpp ../../../source_base/module_external/blas_connector_vector.cpp ../../../source_base/module_external/blas_connector_matrix.cpp diff --git a/source/source_hamilt/module_xc/test/test_xc3.cpp b/source/source_hamilt/module_xc/test/test_xc3.cpp index 3d141accbc..93e9dfc3c0 100644 --- a/source/source_hamilt/module_xc/test/test_xc3.cpp +++ b/source/source_hamilt/module_xc/test/test_xc3.cpp @@ -4,7 +4,6 @@ #include "../exx_info.h" #include "xc3_mock.h" #include "source_base/matrix.h" -#include "source_cell/cal_ux.h" /************************************************ * unit test of functionals @@ -54,7 +53,7 @@ class XCTest_GRADCORR : public XCTest ucell.tpiba = 1; ucell.magnet.lsign_ = true; - unitcell::cal_ux(ucell, 4); + elecstate::cal_ux(ucell, 4); chr.rho = new double*[4]; chr.rho[0] = new double[5]; diff --git a/source/source_hamilt/module_xc/test/test_xc5.cpp b/source/source_hamilt/module_xc/test/test_xc5.cpp index 01aa214fea..e6e0013d4c 100644 --- a/source/source_hamilt/module_xc/test/test_xc5.cpp +++ b/source/source_hamilt/module_xc/test/test_xc5.cpp @@ -5,7 +5,6 @@ #include "../exx_info.h" #include "xc3_mock.h" #include "source_base/matrix.h" -#include "source_cell/cal_ux.h" #include "../../../source_base/parallel_reduce.h" /************************************************ @@ -48,7 +47,7 @@ class XCTest_VXC : public XCTest ucell.tpiba = 1; ucell.magnet.lsign_ = true; - unitcell::cal_ux(ucell, 4); + elecstate::cal_ux(ucell, 4); ucell.omega = 1; chr.rhopw = &(rhopw); @@ -152,7 +151,7 @@ class XCTest_VXC_Libxc : public XCTest ucell.tpiba = 1; ucell.magnet.lsign_ = true; - unitcell::cal_ux(ucell, 4); + elecstate::cal_ux(ucell, 4); ucell.omega = 1; chr.rhopw = &(rhopw); @@ -254,7 +253,7 @@ class XCTest_VXC_meta : public XCTest ucell.tpiba = 1; ucell.magnet.lsign_ = true; - unitcell::cal_ux(ucell, 4); + elecstate::cal_ux(ucell, 4); ucell.omega = 1; chr.rhopw = &(rhopw); diff --git a/source/source_hamilt/module_xc/test/xc3_mock.h b/source/source_hamilt/module_xc/test/xc3_mock.h index afa1017ee5..57625c8c52 100644 --- a/source/source_hamilt/module_xc/test/xc3_mock.h +++ b/source/source_hamilt/module_xc/test/xc3_mock.h @@ -192,7 +192,7 @@ SepPot::~SepPot(){} Sep_Cell::Sep_Cell() noexcept {} Sep_Cell::~Sep_Cell() noexcept {} -namespace unitcell +namespace elecstate { void cal_ux(UnitCell& ucell, const int nspin) { diff --git a/source/source_hamilt/module_xc/xc_functional.h b/source/source_hamilt/module_xc/xc_functional.h index 540469b81e..5808ea2646 100644 --- a/source/source_hamilt/module_xc/xc_functional.h +++ b/source/source_hamilt/module_xc/xc_functional.h @@ -249,15 +249,13 @@ class XC_Functional double* dh, const ModulePW::PW_Basis* rho_basis, const double tpiba); - - static void noncolin_rho( - double* rhoout1, - double* rhoout2, - double* seg, - const double* const* const rho, - const int nrxx, - const double* ux_, - const bool lsign_); + static void noncolin_rho(double* rhoout1, + double* rhoout2, + double* neg, + const double* const* const rho, + const int nrxx, + const double* ux_, + const bool lsign_); //------------------- // xc_lda_exch.cpp diff --git a/source/source_hamilt/module_xc/xc_functional_gga_noncol_sf_builtin.cpp b/source/source_hamilt/module_xc/xc_functional_gga_noncol_sf_builtin.cpp new file mode 100644 index 0000000000..eeaa792c77 --- /dev/null +++ b/source/source_hamilt/module_xc/xc_functional_gga_noncol_sf_builtin.cpp @@ -0,0 +1,308 @@ +#include "xc_functional_gga_noncol_sf_builtin.h" + +#include "source_base/parallel_reduce.h" +#include "source_base/timer.h" +#include "source_base/vector3.h" +#include "source_basis/module_pw/pw_basis.h" +#include "source_cell/unitcell.h" +#include "source_io/module_parameter/parameter.h" +#include "xc_functional.h" + +#include +#include + +namespace ModuleXC +{ +namespace NCGGA_SF_Builtin +{ + +std::tuple v_xc_ncgga_sf_builtin( + const int& nrxx, const double& omega, const double tpiba, const Charge* const chr) +{ + ModuleBase::TITLE("XC_Functional", "v_xc_ncgga_sf_builtin"); + ModuleBase::timer::start("XC_Functional", "v_xc_ncgga_sf_builtin"); + + if (PARAM.inp.nspin != 4 || (!PARAM.globalv.domag && !PARAM.globalv.domag_z)) + throw std::domain_error("v_xc_ncgga_sf_builtin requires NSPIN==4."); + + ModulePW::PW_Basis* rhopw = chr->rhopw; + const int npw = rhopw->npw; + const double e2 = ModuleBase::e2; + constexpr double vanishing = 1e-10; + constexpr double epsr = 1e-6; + const double fac = 0.5; + const bool is_gga = (XC_Functional::get_func_type() == 2 || XC_Functional::get_func_type() == 4); + + std::vector rhotmp1(nrxx), rhotmp2(nrxx), amag(nrxx); + std::vector mag_part(3 * nrxx, 0.0); + + for (int ir = 0; ir < nrxx; ++ir) + { + const double mx = chr->rho[1][ir], my = chr->rho[2][ir], mz = chr->rho[3][ir]; + amag[ir] = std::sqrt(mx * mx + my * my + mz * mz); + rhotmp1[ir] = 0.5 * (chr->rho[0][ir] + amag[ir]); + rhotmp2[ir] = 0.5 * (chr->rho[0][ir] - amag[ir]); + if (amag[ir] > 1e-12) + { + mag_part[ir] = mx / amag[ir]; + mag_part[ir + nrxx] = my / amag[ir]; + mag_part[ir + 2 * nrxx] = mz / amag[ir]; + } + } + for (int ir = 0; ir < nrxx; ++ir) + { + rhotmp1[ir] += fac * chr->rho_core[ir]; + rhotmp2[ir] += fac * chr->rho_core[ir]; + } + + std::vector> rhogsum1(npw), tmp_recip(npw); + rhopw->real2recip(chr->rho[0], rhogsum1.data()); + for (int ig = 0; ig < npw; ++ig) + rhogsum1[ig] += chr->rhog_core[ig]; + + std::vector> gdr1(nrxx), gdr2(nrxx); + std::vector> gdr_mag(nrxx); + XC_Functional::grad_rho(rhogsum1.data(), gdr1.data(), rhopw, tpiba); + + for (int ir = 0; ir < nrxx; ++ir) + { + gdr_mag[ir] = gdr1[ir]; + gdr1[ir] = 0.5 * gdr_mag[ir]; + gdr2[ir] = 0.5 * gdr_mag[ir]; + } + for (int is = 1; is <= 3; ++is) + { + rhopw->real2recip(chr->rho[is], tmp_recip.data()); + XC_Functional::grad_rho(tmp_recip.data(), gdr_mag.data(), rhopw, tpiba); + const double* mp = mag_part.data() + (is - 1) * nrxx; + for (int ir = 0; ir < nrxx; ++ir) + { + const ModuleBase::Vector3 g = 0.5 * gdr_mag[ir] * mp[ir]; + gdr1[ir] += g; + gdr2[ir] -= g; + } + } + + double etxc = 0, vtxc = 0; + ModuleBase::matrix v(PARAM.inp.nspin, nrxx); + + for (int ir = 0; ir < nrxx; ++ir) + { + const double arho = std::abs(chr->rho[0][ir] + chr->rho_core[ir]); + if (arho <= vanishing) + continue; + + double zeta = amag[ir] / arho; + if (std::abs(zeta) > 1.0) + zeta = (zeta > 0) ? 1.0 : -1.0; + double exc = 0, vxc[2] = {0, 0}; + XC_Functional::xc_spin(arho, zeta, exc, vxc[0], vxc[1]); + + v(0, ir) = e2 * 0.5 * (vxc[0] + vxc[1]); + vtxc += v(0, ir) * chr->rho[0][ir]; + + if (amag[ir] > vanishing) + { + const double vs = e2 * 0.5 * (vxc[0] - vxc[1]); + const double inv_a = 1.0 / amag[ir]; + for (int mu = 1; mu < 4; ++mu) + { + v(mu, ir) = vs * chr->rho[mu][ir] * inv_a; + vtxc += v(mu, ir) * chr->rho[mu][ir]; + } + } + etxc += e2 * exc * arho; + } + + if (is_gga) + { + double etxcgc = 0, vtxcgc = 0; + std::vector vup_gga(nrxx, 0), vdw_gga(nrxx, 0); + std::vector> h1(nrxx), h2(nrxx); + + for (int ir = 0; ir < nrxx; ++ir) + { + double sx = 0, v1xup = 0, v1xdw = 0, v2xup = 0, v2xdw = 0; + double sc = 0, v1cup = 0, v1cdw = 0, v2c = 0; + const double grho2a = gdr1[ir] * gdr1[ir]; + const double grho2b = gdr2[ir] * gdr2[ir]; + const double rh = rhotmp1[ir] + rhotmp2[ir]; + + XC_Functional::gcx_spin(rhotmp1[ir], rhotmp2[ir], grho2a, grho2b, + sx, v1xup, v1xdw, v2xup, v2xdw); + + if (rh > epsr) + { + double zeta = (rhotmp1[ir] - rhotmp2[ir]) / rh; + zeta = std::fabs(zeta); + const double grh2 = (gdr1[ir] + gdr2[ir]) * (gdr1[ir] + gdr2[ir]); + XC_Functional::gcc_spin(rh, zeta, grh2, sc, v1cup, v1cdw, v2c); + } + + vup_gga[ir] = e2 * (v1xup + v1cup); + vdw_gga[ir] = e2 * (v1xdw + v1cdw); + + const double v2cup = v2c, v2cdw = v2c, v2cud = v2c; + h1[ir] = e2 * ((v2xup + v2cup) * gdr1[ir] + v2cud * gdr2[ir]); + h2[ir] = e2 * ((v2xdw + v2cdw) * gdr2[ir] + v2cud * gdr1[ir]); + + vtxcgc += vup_gga[ir] * (rhotmp1[ir] - chr->rho_core[ir] * fac); + vtxcgc += vdw_gga[ir] * (rhotmp2[ir] - chr->rho_core[ir] * fac); + etxcgc += e2 * (sx + sc); + } + + for (int ir = 0; ir < nrxx; ++ir) + { + v(0, ir) += 0.5 * (vup_gga[ir] + vdw_gga[ir]); + const double vdiff = 0.5 * (vup_gga[ir] - vdw_gga[ir]); + for (int mu = 1; mu < 4; ++mu) + { + v(mu, ir) += vdiff * mag_part[ir + (mu - 1) * nrxx]; + } + } + + std::vector dh(nrxx); + std::vector> tmp_h(nrxx); + + for (int ir = 0; ir < nrxx; ++ir) + tmp_h[ir] = 0.5 * (h1[ir] + h2[ir]); + XC_Functional::grad_dot(tmp_h.data(), dh.data(), rhopw, tpiba); + for (int ir = 0; ir < nrxx; ++ir) + v(0, ir) -= dh[ir]; + double sum = 0; + for (int ir = 0; ir < nrxx; ++ir) + sum += dh[ir] * chr->rho[0][ir]; + vtxcgc -= sum; + + for (int mu = 1; mu < 4; ++mu) + { + const double* mp = mag_part.data() + (mu - 1) * nrxx; + for (int ir = 0; ir < nrxx; ++ir) + tmp_h[ir] = 0.5 * (h1[ir] - h2[ir]) * mp[ir]; + XC_Functional::grad_dot(tmp_h.data(), dh.data(), rhopw, tpiba); + for (int ir = 0; ir < nrxx; ++ir) + v(mu, ir) -= dh[ir]; + double sum_mu = 0; + for (int ir = 0; ir < nrxx; ++ir) + sum_mu += dh[ir] * chr->rho[mu][ir]; + vtxcgc -= sum_mu; + } + + etxc += etxcgc; + vtxc += vtxcgc; + } + +#ifdef __MPI + Parallel_Reduce::reduce_pool(etxc); + Parallel_Reduce::reduce_pool(vtxc); +#endif + etxc *= omega / rhopw->nxyz; + vtxc *= omega / rhopw->nxyz; + + ModuleBase::timer::end("XC_Functional", "v_xc_ncgga_sf_builtin"); + return std::make_tuple(etxc, vtxc, std::move(v)); +} + +void gradcorr_ncgga_sf_builtin(const Charge* const chr, ModulePW::PW_Basis* rhopw, + const UnitCell* ucell, std::vector& stress_gga) +{ + stress_gga.assign(9, 0.0); + + const int nrxx = rhopw->nrxx; + const int npw = rhopw->npw; + const double e2 = ModuleBase::e2; + constexpr double epsr = 1.0e-6; + constexpr double small = 1.0e-10; + const double fac = 0.5; + + std::vector rhotmp1(nrxx), rhotmp2(nrxx), amag_arr(nrxx); + std::vector mag_part(3 * nrxx, 0.0); + + for (int ir = 0; ir < nrxx; ++ir) + { + const double mx = chr->rho[1][ir], my = chr->rho[2][ir], mz = chr->rho[3][ir]; + double amag = std::sqrt(mx * mx + my * my + mz * mz); + amag_arr[ir] = amag; + rhotmp1[ir] = 0.5 * (chr->rho[0][ir] + amag); + rhotmp2[ir] = 0.5 * (chr->rho[0][ir] - amag); + if (amag > 1e-12) + { + mag_part[ir] = mx / amag; + mag_part[ir + nrxx] = my / amag; + mag_part[ir + 2 * nrxx] = mz / amag; + } + } + for (int ir = 0; ir < nrxx; ++ir) + { + rhotmp1[ir] += fac * chr->rho_core[ir]; + rhotmp2[ir] += fac * chr->rho_core[ir]; + } + + std::vector> rhogsum1(npw), tmp_recip(npw); + rhopw->real2recip(chr->rho[0], rhogsum1.data()); + for (int ig = 0; ig < npw; ++ig) + rhogsum1[ig] += chr->rhog_core[ig]; + + std::vector> gdr1(nrxx), gdr2(nrxx); + std::vector> gdr_mag(nrxx); + XC_Functional::grad_rho(rhogsum1.data(), gdr1.data(), rhopw, ucell->tpiba); + + for (int ir = 0; ir < nrxx; ++ir) + { + gdr_mag[ir] = gdr1[ir]; + gdr1[ir] = 0.5 * gdr_mag[ir]; + gdr2[ir] = 0.5 * gdr_mag[ir]; + } + for (int is = 1; is <= 3; ++is) + { + rhopw->real2recip(chr->rho[is], tmp_recip.data()); + XC_Functional::grad_rho(tmp_recip.data(), gdr_mag.data(), rhopw, ucell->tpiba); + const double* mp = mag_part.data() + (is - 1) * nrxx; + for (int ir = 0; ir < nrxx; ++ir) + { + const ModuleBase::Vector3 g = 0.5 * gdr_mag[ir] * mp[ir]; + gdr1[ir] += g; + gdr2[ir] -= g; + } + } + + for (int ir = 0; ir < nrxx; ++ir) + { + double sx = 0, v1xup = 0, v1xdw = 0, v2xup = 0, v2xdw = 0; + double sc = 0, v1cup = 0, v1cdw = 0, v2c = 0; + double v2cup = 0, v2cdw = 0, v2cud = 0; + const double grho2a = gdr1[ir] * gdr1[ir]; + const double grho2b = gdr2[ir] * gdr2[ir]; + const double rh = rhotmp1[ir] + rhotmp2[ir]; + + XC_Functional::gcx_spin(rhotmp1[ir], rhotmp2[ir], grho2a, grho2b, + sx, v1xup, v1xdw, v2xup, v2xdw); + + if (rh > epsr) + { + double zeta = (rhotmp1[ir] - rhotmp2[ir]) / rh; + zeta = std::fabs(zeta); + const double grh2 = (gdr1[ir] + gdr2[ir]) * (gdr1[ir] + gdr2[ir]); + XC_Functional::gcc_spin(rh, zeta, grh2, sc, v1cup, v1cdw, v2c); + v2cup = v2c; + v2cdw = v2c; + v2cud = v2c; + } + + double tt1[3] = {gdr1[ir].x, gdr1[ir].y, gdr1[ir].z}; + double tt2[3] = {gdr2[ir].x, gdr2[ir].y, gdr2[ir].z}; + for (int l = 0; l < 3; l++) + { + for (int m = 0; m < l + 1; m++) + { + int ind = l * 3 + m; + stress_gga[ind] += e2 * (tt1[l] * tt1[m] * v2xup + tt2[l] * tt2[m] * v2xdw); + stress_gga[ind] += e2 * (tt1[l] * tt1[m] * v2cup + tt2[l] * tt2[m] * v2cdw + + (tt1[l] * tt2[m] + tt2[l] * tt1[m]) * v2cud); + } + } + } +} + +} // namespace NCGGA_SF_Builtin +} // namespace ModuleXC diff --git a/source/source_hamilt/module_xc/xc_functional_gga_noncol_sf_builtin.h b/source/source_hamilt/module_xc/xc_functional_gga_noncol_sf_builtin.h new file mode 100644 index 0000000000..fa615224ee --- /dev/null +++ b/source/source_hamilt/module_xc/xc_functional_gga_noncol_sf_builtin.h @@ -0,0 +1,30 @@ +#ifndef XC_FUNCTIONAL_GGA_NONCOL_SF_BUILTIN_H +#define XC_FUNCTIONAL_GGA_NONCOL_SF_BUILTIN_H + +#include "source_base/matrix.h" +#include "source_estate/module_charge/charge.h" + +#include +#include + +namespace ModulePW +{ +class PW_Basis; +} +struct UnitCell; + +namespace ModuleXC +{ +namespace NCGGA_SF_Builtin +{ + +std::tuple v_xc_ncgga_sf_builtin( + const int& nrxx, const double& omega, const double tpiba, const Charge* const chr); + +void gradcorr_ncgga_sf_builtin(const Charge* const chr, ModulePW::PW_Basis* rhopw, + const UnitCell* ucell, std::vector& stress_gga); + +} // namespace NCGGA_SF_Builtin +} // namespace ModuleXC + +#endif diff --git a/source/source_hamilt/module_xc/xc_pot.cpp b/source/source_hamilt/module_xc/xc_pot.cpp index 1f8a3dcd32..6627b22698 100644 --- a/source/source_hamilt/module_xc/xc_pot.cpp +++ b/source/source_hamilt/module_xc/xc_pot.cpp @@ -9,6 +9,8 @@ #include "source_io/module_parameter/parameter.h" #include "xc_functional.h" +#include "xc_functional_gga_noncol_sf_builtin.h" + #ifdef __LIBXC #include "libxc_abacus.h" #ifdef __EXX @@ -16,6 +18,8 @@ #endif #endif + + // [etxc, vtxc, v] = XC_Functional::v_xc(...) std::tuple XC_Functional::v_xc( const int& nrxx, @@ -48,6 +52,11 @@ std::tuple XC_Functional::v_xc( #endif } + if (PARAM.inp.nspin == 4 && (PARAM.globalv.domag || PARAM.globalv.domag_z) && PARAM.inp.gga_grad == 3) + { + return ModuleXC::NCGGA_SF_Builtin::v_xc_ncgga_sf_builtin(nrxx, ucell->omega, ucell->tpiba, chr); + } + ModuleBase::timer::start("XC_Functional", "v_xc"); //Exchange-Correlation potential Vxc(r) from n(r) diff --git a/source/source_io/module_parameter/input_parameter.h b/source/source_io/module_parameter/input_parameter.h index 127e459548..44cc83a15e 100644 --- a/source/source_io/module_parameter/input_parameter.h +++ b/source/source_io/module_parameter/input_parameter.h @@ -69,7 +69,7 @@ struct Input_para std::string kmesh_type = "gamma"; ///< k-point mesh type for kspacing-generated k-point mesh: gamma or mp double min_dist_coef = 0.2; ///< allowed minimum distance between two atoms - std::string device = "cpu"; + std::string device = "auto"; std::string precision = "double"; std::string gint_precision = "double"; bool timer_enable_nvtx = false; @@ -87,6 +87,8 @@ struct Input_para bool pseudo_mesh = false; ///< 0: use msh to normalize radial wave functions; 1: ///< use mesh, which is used in QE. int nspin = 1; ///< LDA ; LSDA ; non-linear spin + int gga_grad = 3; ///< GGA gradient method for noncollinear spin (nspin=4): 1=collinear approx, 2=projected div(h), 3=Scalmani-Frisch transform (default, most accurate) + std::string device_memory_mode = ""; ///< GPU memory mode: "" (auto), "full_gpu" (all on GPU), "paged" (CPU storage + k-point paging) int pw_diag_nmax = 50; double pw_diag_thr = 0.01; ///< used in cg method bool diago_smooth_ethr = false; ///< smooth ethr for iter methods @@ -617,18 +619,23 @@ struct Input_para bool decay_grad_switch = false; ///< the switch to use the local approximation of gradient ///< decay, 0: no local approximation; 1: apply the method double sc_thr = 1e-06; ///< threshold for spin-constrained DFT in uB - int nsc = 100; ///< maximum number of inner lambda loop + int nsc = 5; ///< maximum number of inner lambda loop int nsc_min = 2; ///< minimum number of inner lambda loop double alpha_trial = 0.01; ///< initial trial step size for lambda in eV/uB^2 double sccut = 3.0; ///< restriction of step size in eV/uB - double sc_scf_thr = 1e-3; ///< minimum number of outer scf loop before initial lambda loop + double sc_scf_thr = 10; ///< density error threshold for activating the lambda loop in spin-constrained DFT + std::string sc_scf_thr_mode = "immediate"; ///< controls when the lambda loop activates: "threshold" (drho=2), or "off" (never activate) + int sc_dir_phase1_steps = 5; ///< number of Phase 1 iterations in direction_only two-phase strategy for collinear (nspin=2) double sc_drop_thr = 1e-3; ///< threshold for lambda-loop threshold cutoff in spin-constrained DFT - std::string sc_lambda_strategy = "bfgs"; ///< lambda update strategy: bfgs, bfgs2, linear_response, augmented_lagrangian, hybrid_delayed, linear_scan + std::string sc_lambda_strategy = "bfgs"; ///< lambda update strategy: bfgs or linear_scan bool sc_direction_only = false; ///< only optimize the direction of magnetization // linear_scan parameters double sc_scan_lambda_start = 0.0; ///< start value for lambda scan (eV/uB) double sc_scan_lambda_end = 1.0; ///< end value for lambda scan (eV/uB) int sc_scan_steps = 20; ///< number of steps in lambda scan + std::string sc_strategy = "normal"; ///< DeltaSpin execution strategy: "fast" (subspace), "accuracy" (full), "normal" (threshold-triggered) + std::string sc_acceleration_mode = "off"; ///< acceleration mode: "off", "first_order", "subspace" + double sc_acceleration_rms_thr = -1.0; ///< RMS threshold (uB) to activate acceleration, <0 disables // ============== #Parameters (18.Quasiatomic Orbital analysis) ========= ///<========================================================== @@ -725,12 +732,14 @@ struct Input_para * * Likewise, the correlation part can be found in corresponding files. * - * These vectors are empty unless the user explicitly requests an - * override. This leaves the version-specific default parameters under - * Libxc's control. + * PBE functional is used as the default functional for XCPNet. */ - std::vector xc_exch_ext = {}; - std::vector xc_corr_ext = {}; + // src/gga_x_pbe.c + std::vector xc_exch_ext = { + 101, 0.8040, 0.2195149727645171}; + // src/gga_c_pbe.c + std::vector xc_corr_ext = { + 130, 0.06672455060314922, 0.031090690869654895034, 1.00000}; // ============== #Parameters (24.td-ofdft) =========================== bool of_cd = false; ///< add CD potential or not diff --git a/source/source_io/module_parameter/read_input_item_elec_stru.cpp b/source/source_io/module_parameter/read_input_item_elec_stru.cpp index bc5ff4f4de..87f290f21f 100644 --- a/source/source_io/module_parameter/read_input_item_elec_stru.cpp +++ b/source/source_io/module_parameter/read_input_item_elec_stru.cpp @@ -53,9 +53,9 @@ void ReadInput::item_elec_stru() For plane-wave basis, * cg: The conjugate-gradient (CG) method. +* bpcg: The BPCG method, which is a block-parallel Conjugate Gradient (CG) method, typically exhibits higher acceleration in a GPU environment. * dav: The Davidson algorithm. * dav_subspace: The Davidson algorithm without orthogonalization operation, this method is the most recommended for efficiency. `pw_diag_ndim` can be set to 2 for this method. -* bpcg: The BPCG method, which is a block-parallel Conjugate Gradient (CG) method, typically exhibits higher acceleration in a GPU environment. The BPCG method is currently under testing and is not recommended for use. For numerical atomic orbitals basis, @@ -340,6 +340,28 @@ The other way is only available when compiling with LIBXC, and it allows for sup read_sync_string(input.dft_functional); this->add_item(item); } + { + Input_Item item("gga_grad"); + item.annotation = "1: collinear approximation (grad|m|); 2: projected div(h); 3: Scalmani-Frisch transform (default, most accurate)"; + item.category = "Electronic structure"; + item.type = "Integer"; + item.description = R"(Method for computing GGA gradients in noncollinear spin (nspin=4) calculations. +* 1: Collinear approximation. Treats nspin=4 as nspin=2 with |m|, computes grad(rho +/- |m|). Fastest but least accurate. +* 2: Improved gradient. Computes grad(m_mu)*m_hat_mu and projects div(h) via m_hat. Intermediate accuracy. +* 3: Scalmani-Frisch transform. Applies SF rotation to decompose gradients per spin channel and retains all cross-terms in div(h). Most accurate. +Only relevant when nspin=4 and a GGA or hybrid functional is used.)"; + item.default_value = "3"; + item.unit = ""; + item.availability = "Only relevant for nspin=4 with GGA/hybrid functionals"; + read_sync_int(input.gga_grad); + item.check_value = [](const Input_Item&, const Parameter& para) { + if (para.input.gga_grad < 1 || para.input.gga_grad > 3) + { + ModuleBase::WARNING_QUIT("ReadInput", "gga_grad must be 1, 2, or 3."); + } + }; + this->add_item(item); + } { Input_Item item("xc_temperature"); item.annotation = "temperature for finite temperature functionals"; @@ -354,14 +376,14 @@ The other way is only available when compiling with LIBXC, and it allows for sup } { Input_Item item("xc_exch_ext"); - item.annotation = "customize Libxc exchange functional parameters"; + item.annotation = "placeholder for xcpnet exchange functional"; item.category = "Electronic structure"; item.type = "Integer followed by Real values"; - item.description = "Customized parameterization of the exchange part of an XC functional. The first value should be the Libxc ID of the original functional, followed by the complete list of external parameters required by the linked Libxc version. If unset, Libxc's own default parameters are used. For functional IDs and parameter definitions, refer to the Libxc documentation and source code." + item.description = "Customized parameterization on the exchange part of XC functional. The first value should be the LibXC ID of the original functional, and latter values are external parameters. Default values are those of Perdew-Burke-Ernzerhof (PBE) functional. For more information on LibXC ID of functionals, please refer to LibXC. For parameters of functionals of interest, please refer to the source code of LibXC, such as PBE functional interface in LibXC: gga_x_pbe.c." "\n\n[NOTE] Solely setting this keyword will take no effect on XC functionals. One should also set " "dft_functional to the corresponding functional to apply the customized parameterization. " "Presently this feature can only support parameterization on one exchange functional."; - item.default_value = ""; + item.default_value = "101 0.8040 0.2195149727645171"; item.unit = ""; item.availability = ""; item.read_value = [](const Input_Item& item, Parameter& para) { @@ -371,15 +393,10 @@ The other way is only available when compiling with LIBXC, and it allows for sup [](const std::string& str) { return std::stod(str); }); }; item.check_value = [](const Input_Item& item, const Parameter& para) { - if (!item.is_read()) + // at least one value should be set + if (para.input.xc_exch_ext.empty()) { - return; - } - if (para.input.xc_exch_ext.size() < 2) - { - ModuleBase::WARNING_QUIT( - "ReadInput", - "xc_exch_ext requires a Libxc ID followed by external parameters."); + ModuleBase::WARNING_QUIT("ReadInput", "xc_exch_ext should not be empty."); } // the first value is actually an integer, not a double const double libxc_id_dbl = para.input.xc_exch_ext[0]; @@ -402,14 +419,14 @@ The other way is only available when compiling with LIBXC, and it allows for sup } { Input_Item item("xc_corr_ext"); - item.annotation = "customize Libxc correlation functional parameters"; + item.annotation = "placeholder for xcpnet exchange functional"; item.category = "Electronic structure"; item.type = "Integer followed by Real values"; - item.description = "Customized parameterization of the correlation part of an XC functional. The first value should be the Libxc ID of the original functional, followed by the complete list of external parameters required by the linked Libxc version. If unset, Libxc's own default parameters are used. For functional IDs and parameter definitions, refer to the Libxc documentation and source code." + item.description = "Customized parameterization on the correlation part of XC functional. The first value should be the LibXC ID of the original functional, and latter values are external parameters. Default values are those of Perdew-Burke-Ernzerhof (PBE) functional. For more information on LibXC ID of functionals, please refer to LibXC. For parameters of functionals of interest, please refer to the source code of LibXC, such as PBE functional interface in LibXC: gga_c_pbe.c." "\n\n[NOTE] Solely setting this keyword will take no effect on XC functionals. One should also set " "dft_functional to the corresponding functional to apply the customized parameterization. " "Presently this feature can only support parameterization on one correlation functional."; - item.default_value = ""; + item.default_value = "130 0.06672455060314922 0.031090690869654895034 1.0"; item.unit = ""; item.availability = ""; item.read_value = [](const Input_Item& item, Parameter& para) { @@ -419,15 +436,10 @@ The other way is only available when compiling with LIBXC, and it allows for sup [](const std::string& str) { return std::stod(str); }); }; item.check_value = [](const Input_Item& item, const Parameter& para) { - if (!item.is_read()) + // at least one value should be set + if (para.input.xc_corr_ext.empty()) { - return; - } - if (para.input.xc_corr_ext.size() < 2) - { - ModuleBase::WARNING_QUIT( - "ReadInput", - "xc_corr_ext requires a Libxc ID followed by external parameters."); + ModuleBase::WARNING_QUIT("ReadInput", "xc_corr_ext should not be empty."); } // the first value is actually an integer, not a double const double libxc_id_dbl = para.input.xc_corr_ext[0]; @@ -675,16 +687,24 @@ For systems that are difficult to converge, one could try increasing the value o item.availability = ""; read_sync_double(input.mixing_restart); item.reset_value = [](const Input_Item& item, Parameter& para) { - if (para.input.sc_mag_switch == 1) - {// for DeltaSpin calculation, the mixing_restart should be same as sc_scf_thr - if(para.input.sc_scf_thr != 10.0) + if (para.input.sc_mag_switch) + { + if (para.input.sc_direction_only) + { + para.input.mixing_restart = 0.0; + } + else if (para.input.sc_scf_thr_mode == "threshold") { para.input.mixing_restart = para.input.sc_scf_thr; } - else - {// no mixing_restart until oscillation happen in PW base + else if (para.input.sc_scf_thr_mode == "immediate") + { para.input.mixing_restart = para.input.scf_thr / 10.0; } + else // "off" + { + para.input.mixing_restart = 0.0; + } } }; this->add_item(item); From 1c200ab49ab379201e90c09939b0c2d7e569cad5 Mon Sep 17 00:00:00 2001 From: dyzheng Date: Sun, 2 Aug 2026 12:19:49 +0800 Subject: [PATCH 2/3] Fix(gga_grad): correct vtxc contribution, revert unrelated changes - Fix libxc_tools.cpp vtxc: use chr->rho[is][ir] (nspin=4 density) instead of libxc-format rho (up/dn spin density) - Revert libxc_setup.cpp to develop (the accel branch changes had map key collision, PBEH param count error, and weakened param validation) - Revert input_parameter.h: restore develop defaults for device, sc_* params, xc_exch_ext/xc_corr_ext; keep only gga_grad=3 - Revert docs to develop version --- docs/advanced/input_files/input-main.md | 145 +++++++--- docs/advanced/scf/spin.md | 265 ------------------ docs/parameters.yaml | 142 +++++----- .../source_hamilt/module_xc/libxc_setup.cpp | 77 +++-- .../source_hamilt/module_xc/libxc_tools.cpp | 4 +- .../module_parameter/input_parameter.h | 26 +- .../read_input_item_elec_stru.cpp | 91 +++--- 7 files changed, 258 insertions(+), 492 deletions(-) diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index 94aa78b935..f04f903448 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -168,7 +168,9 @@ - [out\_proj\_band](#out_proj_band) - [out\_stru](#out_stru) - [out\_level](#out_level) + - [out\_hsk](#out_hsk) - [out\_mat\_hs](#out_mat_hs) + - [out\_hsr](#out_hsr) - [out\_mat\_hs2](#out_mat_hs2) - [out\_mat\_tk](#out_mat_tk) - [out\_mat\_r](#out_mat_r) @@ -408,6 +410,9 @@ - [sc\_scf\_thr](#sc_scf_thr) - [sc\_direction\_only](#sc_direction_only) - [sc\_lambda\_strategy](#sc_lambda_strategy) + - [sc\_scan\_lambda\_start](#sc_scan_lambda_start) + - [sc\_scan\_lambda\_end](#sc_scan_lambda_end) + - [sc\_scan\_steps](#sc_scan_steps) - [vdW correction](#vdw-correction) - [vdw\_method](#vdw_method) - [vdw\_d4\_xc](#vdw_d4_xc) @@ -582,7 +587,7 @@ - md: perform molecular dynamics simulations - get_pchg: obtain partial (band-decomposed) charge densities (for LCAO basis only). See out_pchg for more information - get_wf: obtain real space wave functions (for LCAO basis only). See out_wfc_norm and out_wfc_re_im for more information - - get_s: obtain the overlap matrix formed by localized orbitals (for LCAO basis with multiple k points). the file name is SR.csr with file format being the same as that generated by out_mat_hs2 + - get_s: obtain the overlap matrix formed by localized orbitals (for LCAO basis with multiple k points). The file name is OUT.${suffix}/sr_nao.csr, with the same file format as generated by out_hsr 1 - gen_bessel: generates projectors, i.e., a series of Bessel functions, for the DeePKS method (for LCAO basis only) - gen_opt_abfs: generate opt-ABFs as discussed in this article - test_memory: obtain a rough estimation of memory consumption for the calculation @@ -611,7 +616,7 @@ - **Description**: Takes value 1, 0 or -1. - -1: No symmetry will be considered. It is recommended to set -1 for non-colinear + soc calculations, where time reversal symmetry is broken sometimes. - 0: Only time reversal symmetry would be considered in symmetry operations, which implied k point and -k point would be treated as a single k point with twice the weight. - - 1: Symmetry analysis will be performed to determine the type of Bravais lattice and associated symmetry operations. (point groups, space groups, primitive cells, and irreducible k-points) + - 1: Symmetry analysis will be performed to determine the type of Bravais lattice and associated symmetry operations (point groups, space groups, primitive cells, and irreducible k-points). For a magnetic system, the symmetry of the initial magnetic structure will be analyzed and preserved. > Note: When symmetry is enabled (value 1), k-points are reduced to the irreducible Brillouin zone (IBZ). For explicit k-point lists with custom weights (see KPT file), the custom weights are preserved during symmetry reduction. For Monkhorst-Pack grids, uniform weights are used. - **Default**: default @@ -1149,9 +1154,9 @@ For plane-wave basis, - cg: The conjugate-gradient (CG) method. - - bpcg: The BPCG method, which is a block-parallel Conjugate Gradient (CG) method, typically exhibits higher acceleration in a GPU environment. - dav: The Davidson algorithm. - dav_subspace: The Davidson algorithm without orthogonalization operation, this method is the most recommended for efficiency. `pw_diag_ndim` can be set to 2 for this method. + - bpcg: The BPCG method, which is a block-parallel Conjugate Gradient (CG) method, typically exhibits higher acceleration in a GPU environment. The BPCG method is currently under testing and is not recommended for use. For numerical atomic orbitals basis, @@ -1229,18 +1234,16 @@ ### xc_exch_ext - **Type**: Integer followed by Real values -- **Description**: Customized parameterization on the exchange part of XC functional. The first value should be the LibXC ID of the original functional, and latter values are external parameters. Default values are those of Perdew-Burke-Ernzerhof (PBE) functional. For more information on LibXC ID of functionals, please refer to LibXC. For parameters of functionals of interest, please refer to the source code of LibXC, such as PBE functional interface in LibXC: gga_x_pbe.c. +- **Description**: Customized parameterization of the exchange part of an XC functional. The first value should be the Libxc ID of the original functional, followed by the complete list of external parameters required by the linked Libxc version. If unset, Libxc's own default parameters are used. For functional IDs and parameter definitions, refer to the Libxc documentation and source code. > Note: Solely setting this keyword will take no effect on XC functionals. One should also set dft_functional to the corresponding functional to apply the customized parameterization. Presently this feature can only support parameterization on one exchange functional. -- **Default**: 101 0.8040 0.2195149727645171 ### xc_corr_ext - **Type**: Integer followed by Real values -- **Description**: Customized parameterization on the correlation part of XC functional. The first value should be the LibXC ID of the original functional, and latter values are external parameters. Default values are those of Perdew-Burke-Ernzerhof (PBE) functional. For more information on LibXC ID of functionals, please refer to LibXC. For parameters of functionals of interest, please refer to the source code of LibXC, such as PBE functional interface in LibXC: gga_c_pbe.c. +- **Description**: Customized parameterization of the correlation part of an XC functional. The first value should be the Libxc ID of the original functional, followed by the complete list of external parameters required by the linked Libxc version. If unset, Libxc's own default parameters are used. For functional IDs and parameter definitions, refer to the Libxc documentation and source code. > Note: Solely setting this keyword will take no effect on XC functionals. One should also set dft_functional to the corresponding functional to apply the customized parameterization. Presently this feature can only support parameterization on one correlation functional. -- **Default**: 130 0.06672455060314922 0.031090690869654895034 1.0 ### pseudo_rcut @@ -1906,7 +1909,7 @@ - **Type**: Boolean \[Integer\](optional) - **Availability**: *Numerical atomic orbital basis (multi-k points)* -- **Description**: Whether to output the density matrix with Bravias lattice vector R index into files in the folder OUT.${suffix}. The files are named as dmr{s}{spin index}{g}{geometry index}{_nao} + {".csr"}. Here, 's' refers to spin, where s1 means spin up channel while s2 means spin down channel, and the sparse matrix format 'csr' is mentioned in out_mat_hs2. Finally, if out_app_flag is set to false, the file name contains the optional 'g' index for each ionic step that may have different geometries, and if out_app_flag is set to true, the density matrix with respect to Bravias lattice vector R accumulates during ionic steps: +- **Description**: Whether to output the density matrix with Bravias lattice vector R index into files in the folder OUT.${suffix}. The files are named as dmr{s}{spin index}{g}{geometry index}{_nao} + {".csr"}. Here, 's' refers to spin, where s1 means spin up channel while s2 means spin down channel, and the sparse matrix format 'csr' is mentioned in out_hsr. Finally, if out_app_flag is set to false, the file name contains the optional 'g' index for each ionic step that may have different geometries, and if out_app_flag is set to true, the density matrix with respect to Bravias lattice vector R accumulates during ionic steps: - nspin = 1: dmrs1_nao.csr; - nspin = 2: dmrs1_nao.csr and dmrs2_nao.csr for the two spin channels. @@ -2003,32 +2006,60 @@ - m: molecular dynamics level, which does not print some information for simplicity. - **Default**: ie -### out_mat_hs +### out_hsk -- **Type**: Boolean \[Integer\](optional) +- **Type**: Integer \[Integer\](optional) - **Availability**: *Numerical atomic orbital basis* -- **Description**: Whether to print the upper triangular part of the Hamiltonian matrices and overlap matrices for each k-point into files in the directory OUT.${suffix}. The second number controls precision. For more information, please refer to hs_matrix.md. Also controled by out_freq_ion and out_app_flag. - - For gamma only case: - - nspin = 1: hks1_nao.txt for the Hamiltonian matrix and sks1_nao.txt for the overlap matrix; - - nspin = 2: hks1_nao.txt and hks2_nao.txt for the Hamiltonian matrix and sks1_nao.txt for the overlap matrix. Note that the code will not output sks2_nao.txt because it is the same as sks1_nao.txt; - - nspin = 4: hks12_nao.txt for the Hamiltonian matrix and sks12_nao.txt for the overlap matrix. - - For multi-k points case: - - nspin = 1: hks1k1_nao.txt for the Hamiltonian matrix at the 1st k-point, and sks1k1_nao.txt for the overlap matrix for the 1st k-point, ...; - - nspin = 2: hks1k1_nao.txt and hks2k1_nao.txt for the two spin channels of the Hamiltonian matrix at the 1st k-point, and sks1k1_nao.txt for the overlap matrix for the 1st k-point. Note that the code will not output sks2k1_nao.txt because it is the same as sks1k1_nao.txt, ...; - - nspin = 4: hks12k1_nao.txt for the Hamiltonian matrix at the 1st k-point, and sks12k1_nao.txt for the overlap matrix for the 1st k-point, ...; +- **Description**: Output the upper triangular part of the Hamiltonian and overlap matrices in reciprocal space for each k-point into files in the directory OUT.${suffix}. The first integer selects the format: + - 0: disabled; + - 1: text output; the optional second integer controls precision and defaults to 8; + - 2: reserved for binary output, which is not implemented yet; + - 3: NPZ output, which is not implemented for H(k)/S(k). + + The output is also controlled by out_freq_ion and out_app_flag. For more information, refer to hs_matrix.md. + + - Gamma-only, nspin = 1: hk_nao.txt for the Hamiltonian matrix and sk_nao.txt for the overlap matrix. + - Gamma-only, nspin = 2: hks1_nao.txt and hks2_nao.txt for the two spin channels of the Hamiltonian matrix, and sk_nao.txt for the overlap matrix. Only one overlap matrix is written because it is identical for both spin channels. + - Gamma-only, nspin = 4: not available with the gamma-only algorithm. + - Multi-k, nspin = 1: hk1_nao.txt for the Hamiltonian matrix and sk1_nao.txt for the overlap matrix at the first k-point. + - Multi-k, nspin = 2: hk1s1_nao.txt and hk1s2_nao.txt for the two spin channels of the Hamiltonian matrix, and sk1_nao.txt for the overlap matrix at the first k-point. Only one overlap matrix is written because it is identical for both spin channels. + - Multi-k, nspin = 4: hk1s4_nao.txt for the spinor Hamiltonian matrix and sk1_nao.txt for the spinor overlap matrix at the first k-point. + When out_app_flag is false, g followed by the one-based ionic-step index is inserted before _nao, for example hk1s1g1_nao.txt. > Note: In the 3.10-LTS version, the file names are data-0-H and data-0-S, etc. -- **Default**: False 8 +- **Default**: 0 8 - **Unit**: Ry -### out_mat_hs2 +### out_mat_hs - **Type**: Boolean \[Integer\](optional) -- **Availability**: *Numerical atomic orbital basis (not gamma-only algorithm)* -- **Description**: Whether to print files containing the Hamiltonian matrix and overlap matrix into files in the directory OUT.${suffix}. For more information, please refer to hs_matrix.md. +- **Availability**: *Numerical atomic orbital basis* +- **Description**: Legacy alias for out_hsk 1, which outputs Hamiltonian and overlap matrices in reciprocal space for each k-point. The optional second integer controls text precision. If both out_hsk and out_mat_hs are present, out_hsk takes precedence. +- **Default**: False 8 +- **Unit**: Ry + +### out_hsr + +- **Type**: Integer \[Integer\](optional) +- **Availability**: *Numerical atomic orbital basis* +- **Description**: Output Hamiltonian and overlap matrices in real space, indexed by the Bravais lattice vector R, in the directory OUT.${suffix}. The first integer selects the format: + - 0: disabled; + - 1: text CSR output; the optional second integer controls precision and defaults to 8; + - 2: reserved for binary output, which is not implemented yet; + - 3: NPZ output using hrs1_nao.npz, hrs2_nao.npz when needed, and sr_nao.npz. + + For multi-k calculations, the output contains the individual real-space blocks stored for the Bravais lattice vectors R. For gamma-only calculations, the internal real-space contributions are folded into a single R = (0, 0, 0) block. This folded result cannot recover the original R-resolved contributions or interpolate arbitrary k points. Terms added only while constructing H(k) are not guaranteed to be present. > Note: In the 3.10-LTS version, the file names are data-HR-sparse_SPIN0.csr and data-SR-sparse_SPIN0.csr, etc. -- **Default**: False [8] +- **Default**: 0 8 +- **Unit**: Ry + +### out_mat_hs2 + +- **Type**: Boolean \[Integer\](optional) +- **Availability**: *Numerical atomic orbital basis* +- **Description**: Legacy alias for out_hsr 1, which outputs Hamiltonian and overlap matrices in real space indexed by the Bravais lattice vector R. The optional second integer controls text precision. If both out_hsr and out_mat_hs2 are present, out_hsr takes precedence. +- **Default**: False 8 - **Unit**: Ry ### out_mat_tk @@ -2055,7 +2086,7 @@ - **Type**: Boolean \[Integer\](optional) - **Availability**: *Numerical atomic orbital basis (not gamma-only algorithm)* -- **Description**: Generate files containing the kinetic energy matrix. The optional second parameter controls text output precision. The format will be the same as the Hamiltonian matrix and overlap matrix as mentioned in out_mat_hs2. The name of the files will be trs1_nao.csr and so on. Also controled by out_freq_ion and out_app_flag. +- **Description**: Generate files containing the kinetic energy matrix. The optional second parameter controls text output precision. The format will be the same as the Hamiltonian matrix and overlap matrix as mentioned in out_hsr. The name of the files will be trs1_nao.csr and so on. Also controled by out_freq_ion and out_app_flag. > Note: In the 3.10-LTS version, the file name is data-TR-sparse_SPIN0.csr. - **Default**: False 8 @@ -2065,7 +2096,7 @@ - **Type**: Integer - **Availability**: *Numerical atomic orbital basis (not gamma-only algorithm)* -- **Description**: Whether to print files containing the derivatives of the Hamiltonian matrix. The format will be the same as the Hamiltonian matrix and overlap matrix as mentioned in out_mat_hs2. The name of the files will be dhrxs1_nao.csr, dhrys1_nao.csr, dhrzs1_nao.csr and so on. Also controled by out_freq_ion and out_app_flag. +- **Description**: Whether to print files containing the derivatives of the Hamiltonian matrix. The format will be the same as the Hamiltonian matrix and overlap matrix as mentioned in out_hsr. The name of the files will be dhrxs1_nao.csr, dhrys1_nao.csr, dhrzs1_nao.csr and so on. Also controled by out_freq_ion and out_app_flag. Format: <enable> [precision] [iat1 iat2 ...]. The first value (0/1) enables/disables output. The second optional value sets the output precision (default: 8). Starting from the third value, 1-based atom indices can be listed to restrict output to derivatives with respect to those specific atoms only; if no atom indices are given, all atoms are written. @@ -2132,7 +2163,7 @@ - **Type**: Integer - **Description**: Whether to print files containing the kinetic energy matrix T(R) in CSR format. - See out_mat_hs2 for format details. + See out_hsr for format details. - **Default**: 0 8 - **Unit**: Ry @@ -2141,7 +2172,7 @@ - **Type**: Integer - **Description**: Whether to print files containing the nonlocal pseudopotential matrix Vnl(R) in CSR format. - See out_mat_hs2 for format details. + See out_hsr for format details. - **Default**: 0 8 - **Unit**: Ry @@ -2150,7 +2181,7 @@ - **Type**: Integer - **Description**: Whether to print files containing the local pseudopotential matrix Vl(R) in CSR format. - See out_mat_hs2 for format details. + See out_hsr for format details. - **Default**: 0 8 - **Unit**: Ry @@ -2159,7 +2190,7 @@ - **Type**: Integer - **Description**: Whether to print files containing the Hartree matrix Vh(R) in CSR format. - See out_mat_hs2 for format details. + See out_hsr for format details. - **Default**: 0 8 - **Unit**: Ry @@ -2168,7 +2199,7 @@ - **Type**: Integer - **Description**: Whether to print files containing the XC matrix Vxc(R) in CSR format. - See out_mat_hs2 for format details. + See out_hsr for format details. - **Default**: 0 8 - **Unit**: Ry @@ -2177,7 +2208,7 @@ - **Type**: Integer - **Description**: Whether to print files containing the exact-exchange matrix Vexx(R) in CSR format. - See out_mat_hs2 for format details. + See out_hsr for format details. - **Default**: 0 8 - **Unit**: Ry @@ -2244,15 +2275,15 @@ - **Type**: Boolean - **Availability**: *Numerical atomic orbital basis (not gamma-only algorithm)* -- **Description**: Whether to print Hamiltonian matrices H(R) in npz format. This feature does not work for gamma-only calculations. +- **Description**: Whether to print Hamiltonian matrices H(R) in NPZ format as hrs1_nao.npz and, for nspin = 2, hrs2_nao.npz. This feature does not work for gamma-only calculations. - **Default**: False - **Unit**: Ry ### out_hsr_npz - **Type**: Boolean -- **Availability**: *Numerical atomic orbital basis (not gamma-only algorithm)* -- **Description**: Whether to print Hamiltonian matrices H(R) and overlap matrix S(R) in npz format. This feature does not work for gamma-only calculations. +- **Availability**: *Numerical atomic orbital basis* +- **Description**: Legacy alias for out_hsr 3, writing hrs1_nao.npz, hrs2_nao.npz when needed, and sr_nao.npz. If both out_hsr and out_hsr_npz are present, out_hsr takes precedence. Gamma-only calculations write the folded R = (0, 0, 0) representation. - **Default**: False - **Unit**: Ry @@ -2274,13 +2305,13 @@ - **Type**: Boolean - **Availability**: *Numerical atomic orbital basis (not gamma-only algorithm)* -- **Description**: Whether to output r(R), H(R), S(R), T(R), dH(R), dS(R), and wfc matrices in an append manner during molecular dynamics calculations. Check input parameters out_mat_r, out_mat_hs2, out_mat_t, out_mat_dh, out_mat_hs and out_wfc_lcao for more information. +- **Description**: Whether to output r(R), H(R), S(R), T(R), dH(R), dS(R), and wfc matrices in an append manner during molecular dynamics calculations. Check input parameters out_mat_r, out_hsr, out_mat_t, out_mat_dh, out_hsk and out_wfc_lcao for more information. - **Default**: true ### out_ndigits - **Type**: Integer -- **Availability**: *out_mat_hs 1 case presently.* +- **Availability**: *out_hsk 1 case presently.* - **Description**: Controls the length of decimal part of output data, such as charge density, Hamiltonian matrix, Overlap matrix and so on. - **Default**: 8 @@ -3676,8 +3707,8 @@ - **Type**: Integer - **Description**: Determines whether to calculate the plus U correction, which is especially important for correlated electrons. - - 1: Calculate plus U correction with radius-adjustable localized projections (with parameter onsite_radius). Supported for both PW and LCAO basis sets. - - 2: Calculate plus U correction using first zeta of NAOs as projections (this is old method for testing). Only available for LCAO basis. + - 1: Calculate plus U correction with radius-adjustable localized projections (with parameter onsite_radius). + - 2: Calculate plus U correction using first zeta of NAOs as projections (this is old method for testing). - 0: Do not calculate plus U correction. - **Default**: 0 @@ -3821,20 +3852,46 @@ - **Type**: Boolean - **Availability**: *sc_mag_switch is true* -- **Description**: When true, only the direction of the magnetic moment is constrained to the target direction, while the magnitude is allowed to vary freely. This is useful for studying magnetic anisotropy or when the magnitude of the moment is determined by the electronic structure rather than an external constraint. When false (default), both the direction and magnitude of the magnetic moment are constrained to the target values. +- **Description**: When true, only the direction of the magnetic moment is constrained to the target direction, while the magnitude is allowed to vary freely. This is useful for studying magnetic anisotropy or when the magnitude of the moment is determined by the electronic structure rather than an external constraint. + + When false (default), both the direction and magnitude of the magnetic moment are constrained to the target values. - **Default**: False ### sc_lambda_strategy - **Type**: String - **Availability**: *sc_mag_switch is true* -- **Description**: Lambda update strategy for spin-constrained DFT. Available options are: - - bfgs: BFGS quasi-Newton method (default, robust and well-tested) - - linear_response: linear response method (Scheme B) - - augmented_lagrangian: augmented Lagrangian method (Scheme C) +- **Description**: Lambda update strategy for spin-constrained DFT: + - bfgs: BFGS quasi-Newton method + - linear_response: linear response (Scheme B) + - augmented_lagrangian: augmented Lagrangian (Scheme C) - hybrid_delayed: hybrid delayed update (Scheme D) + - linear_scan: linear sweep of lambda for testing magnetic moment response - **Default**: bfgs +### sc_scan_lambda_start + +- **Type**: Float +- **Availability**: *sc_lambda_strategy is linear_scan* +- **Description**: Starting lambda value for linear_scan strategy. Only used when sc_lambda_strategy=linear_scan. +- **Default**: 0.0 +- **Unit**: eV/uB + +### sc_scan_lambda_end + +- **Type**: Float +- **Availability**: *sc_lambda_strategy is linear_scan* +- **Description**: Ending lambda value for linear_scan strategy. Only used when sc_lambda_strategy=linear_scan. +- **Default**: 1.0 +- **Unit**: eV/uB + +### sc_scan_steps + +- **Type**: Integer +- **Availability**: *sc_lambda_strategy is linear_scan* +- **Description**: Number of lambda values to scan. Only used when sc_lambda_strategy=linear_scan. +- **Default**: 20 + [back to top](#full-list-of-input-keywords) ## vdW correction diff --git a/docs/advanced/scf/spin.md b/docs/advanced/scf/spin.md index 10c783d868..1749db156d 100644 --- a/docs/advanced/scf/spin.md +++ b/docs/advanced/scf/spin.md @@ -28,224 +28,6 @@ If **"ocp=1"** and **"ocp_set"** is set in INPUT file, the occupations of states 2. **"nupdown"** If **"nupdown"** is set to non-zero, number of spin-up and spin-down electrons will be fixed, and Fermi energy level will split to E_Fermi_up and E_Fermi_down. By the way, total magnetization will also be fixed, and will be the value of **"nupdown"**. -## DeltaSpin (Spin-Constrained DFT) - -DeltaSpin is a spin-constrained DFT method that allows users to constrain the magnetic moments on individual atoms to target values during self-consistent field (SCF) calculations. This is useful for studying magnetic excitations, non-collinear magnetic structures, and systems where the magnetic ground state is not known a priori. - -The theoretical foundation and implementation details can be found in: - -- Cai Z, Wang K, Xu Y, et al., "A self-adaptive first-principles approach for magnetic excited states," *Quantum Frontiers* 2.1 (2023): 21. [DOI: 10.1007/s44214-023-00050-z](https://doi.org/10.1007/s44214-023-00050-z) -- Zheng D, Peng X, Huang Y, et al., "Integrating deep-learning-based magnetic model and non-collinear spin-constrained method: methodology, implementation and application," *npj Computational Materials* (2026). - -### Enabling DeltaSpin - -Set `sc_mag_switch 1` in the INPUT file. DeltaSpin is supported for both PW (`basis_type = pw`) and LCAO (`basis_type = lcao`) basis sets, with `nspin = 2` (collinear) or `nspin = 4` (non-collinear). - -### Specifying Target Magnetic Moments in STRU - -Target magnetic moments and constraint flags are specified per atom in the `ATOMIC_POSITIONS` section of the STRU file, using the `mag` (or `magmom`), `sc`, `lambda`, `angle1`, and `angle2` keywords after the atomic coordinates. - -#### Collinear (nspin=2) - -For collinear spin, only the z-component of the magnetic moment is constrained: - -``` -ATOMIC_POSITIONS -Direct - -Fe -0.0 -2 -0.00 0.00 0.00 mag 2.0 sc 1 -0.51 0.51 0.51 mag -2.0 sc 1 -``` - -- `mag 2.0`: target magnetic moment of 2.0 $\mu_B$ along z-axis -- `sc 1`: constrain the z-component (1 = constrained, 0 = unconstrained) - -#### Non-collinear (nspin=4), vector form - -For non-collinear spin, specify the magnetic moment as a vector (mx, my, mz): - -``` -ATOMIC_POSITIONS -Direct - -Fe -0.0 -2 -0.00 0.00 0.00 mag 2.0 0.0 0.0 sc 1 1 1 -0.51 0.51 0.51 mag 0.0 0.0 -2.0 sc 1 1 1 -``` - -- `mag 2.0 0.0 0.0`: target moment vector in Cartesian coordinates ($\mu_B$) -- `sc 1 1 1`: constrain x, y, z components respectively - -#### Non-collinear (nspin=4), angle form - -Alternatively, use `angle1` (polar angle $\theta$) and `angle2` (azimuthal angle $\phi$) in degrees to specify the direction: - -``` -0.00 0.00 0.00 mag 2.0 angle1 0 angle2 0 sc 1 1 1 -0.51 0.51 0.51 mag 2.0 angle1 180 angle2 0 sc 1 1 1 -``` - -The Cartesian components are computed as: -- $m_z = |\mathbf{m}| \cos\theta$ -- $m_x = |\mathbf{m}| \sin\theta \cos\phi$ -- $m_y = |\mathbf{m}| \sin\theta \sin\phi$ - -#### Providing initial Lagrange multipliers - -Initial lambda values (in eV/$\mu_B$) can be provided via the `lambda` keyword to accelerate convergence: - -``` -0.00 0.00 0.00 mag 2.0 lambda 0.01 0.0 0.0 sc 1 1 1 -``` - -A single value sets $\lambda_z$; three values set $\lambda_x$, $\lambda_y$, $\lambda_z$. - -#### Partial constraints - -Set `sc 0` for unconstrained components. For example, to constrain only the direction but not the magnitude (use with `sc_direction_only`): - -``` -0.00 0.00 0.00 mag 2.0 0.0 0.0 sc 1 1 0 -``` - -### DeltaSpin INPUT Parameters - -| Parameter | Type | Default | Description | -|-----------|------|---------|-------------| -| `sc_mag_switch` | Boolean | False | Enable DeltaSpin | -| `sc_thr` | Real | 1.0e-6 | Convergence criterion for lambda loop (RMS, in $\mu_B$) | -| `nsc` | Integer | 100 | Maximum number of lambda iterations | -| `nsc_min` | Integer | 2 | Minimum number of lambda iterations | -| `sc_scf_nmin` | Integer | 2 | Minimum outer SCF iterations before starting lambda loop | -| `alpha_trial` | Real | 0.01 | Initial trial step size for lambda (eV/$\mu_B^2$) | -| `sccut` | Real | 3.0 | Maximum step size for lambda (eV/$\mu_B$) | -| `sc_drop_thr` | Real | 1.0e-2 | Convergence ratio threshold for adaptive lambda loop | -| `sc_scf_thr` | Real | 1.0e-4 | Density error threshold for entering lambda loop | -| `sc_direction_only` | Boolean | False | Constrain only the direction, not the magnitude | -| `sc_lambda_strategy` | String | bfgs | Lambda update strategy (see below) | -| `decay_grad_switch` | Boolean | False | Enable gradient-based early exit | - -For full parameter details, see the [Spin-Constrained DFT](../input_files/input-main.md#spin-constrained-dft) section of the input keyword list. - -### Lambda Update Strategies - -The `sc_lambda_strategy` parameter controls how the Lagrange multipliers $\lambda$ are updated during the lambda loop: - -- **`bfgs`** (default): BFGS quasi-Newton method with line search. Robust and well-tested for both PW and LCAO. Uses `alpha_trial` and `sccut` to control step size. - -- **`linear_response`**: Linear response method (Scheme B). Estimates the magnetic susceptibility $\chi$ from the history of $(\lambda, M)$ pairs and performs a one-step Newton-like update: $\Delta\lambda = \beta (M_{\text{target}} - M) / \chi$, where $\beta$ is a mixing parameter. - -- **`augmented_lagrangian`**: Augmented Lagrangian method (Scheme C). Uses a penalty parameter $\mu$ that grows over iterations: $\lambda_{\text{new}} = \lambda + \mu (M - M_{\text{target}})$. The penalty increases until convergence is achieved. - -- **`hybrid_delayed`**: Hybrid delayed update (Scheme D). Two-phase approach: in the early phase (SCF not yet converged), lambda updates are gentle; in the late phase (SCF nearly converged), augmented Lagrangian updates are applied. - -### Direction-Only Mode - -When `sc_direction_only 1` is set, only the **direction** of the magnetic moment is constrained to match the target, while the magnitude is allowed to vary freely. This is useful for: - -- Studying magnetic anisotropy energy surfaces -- Cases where the moment magnitude is determined by the electronic structure -- Converging to the easy-axis direction without fixing the moment size - -In this mode, the lambda vector is projected to be perpendicular to the target moment direction at each iteration, ensuring it can only rotate the magnetization, not stretch it. - -### Combining DeltaSpin with DFT+U - -DeltaSpin can be combined with DFT+U for strongly correlated systems. When both `sc_mag_switch` and `dft_plus_u` are enabled: - -1. DFT+U occupation update runs first in each SCF iteration -2. DeltaSpin lambda loop runs after, constraining the magnetic moments -3. The DFT+U-corrected Hamiltonian is used by the lambda loop - -Example INPUT for PW DFT+U + DeltaSpin: - -``` -INPUT_PARAMETERS -calculation scf -basis_type pw -ecutwfc 50 -nspin 2 -dft_plus_u 1 -orbital_corr -1 2 -hubbard_u 0.0 4.0 -sc_mag_switch 1 -sc_thr 1.0e-6 -sc_scf_thr 1.0e-4 -sc_lambda_strategy bfgs -``` - -### Example: Collinear antiferromagnetic Fe - -INPUT file: - -``` -INPUT_PARAMETERS -calculation scf -basis_type pw -ecutwfc 50 -nspin 2 -sc_mag_switch 1 -sc_thr 1.0e-6 -``` - -STRU file: - -``` -ATOMIC_SPECIES -Fe 55.845 Fe.upf - -LATTICE_CONSTANT -8.190 - -LATTICE_VECTORS - 1.00 0.50 0.50 - 0.50 1.00 0.50 - 0.50 0.50 1.00 - -ATOMIC_POSITIONS -Direct - -Fe -0.0 -2 -0.00 0.00 0.00 mag 2.0 sc 1 -0.51 0.51 0.51 mag -2.0 sc 1 -``` - -### Example: Non-collinear constrained moments - -INPUT file: - -``` -INPUT_PARAMETERS -calculation scf -basis_type pw -ecutwfc 50 -nspin 4 -noncolin 1 -sc_mag_switch 1 -sc_direction_only 1 -sc_lambda_strategy bfgs -``` - -STRU file: - -``` -ATOMIC_POSITIONS -Direct - -Fe -0.0 -2 -0.00 0.00 0.00 mag 2.0 0.0 0.0 sc 1 1 0 -0.51 0.51 0.51 mag 0.0 0.0 2.0 sc 1 1 0 -``` - ## Noncollinear Spin Polarized Calculations The spin non-collinear polarization calculation corresponds to setting **"noncolin 1"**, in which case the coupling between spin up and spin down will be taken into account. In this case, nspin is automatically set to 4, which is usually not required to be specified manually. @@ -266,53 +48,6 @@ Note: different settings for "noncolin" and "lspinorb" correspond to different c **Special case**: `noncolin=0, lspinorb=1` is commonly used for non-magnetic materials with SOC effects (e.g., topological insulators, semiconductors with spin-orbit splitting). In this case, the magnetization is NOT automatically set, implying no magnetic moments in the system. -### Choosing gga_grad for Noncollinear Calculations - -When performing noncollinear spin calculations (`nspin=4`), the `gga_grad` parameter controls how the gradient of the magnetization is computed in GGA exchange-correlation functionals. This is critical for obtaining correct magnetic anisotropy energies and spin-orbit coupling effects. - -#### Available Options - -| gga_grad | Gradient Method | Description | When to Use | -|----------|----------------|-------------|-------------| -| 0 | Standard | Conventional GGA gradient (collinear approximation) | Collinear calculations only (`nspin<4`) | -| 1 | $\nabla\|\mathbf{m}\|$ | Gradient of magnetization magnitude | Noncollinear with SOC, when only magnitude variation matters | -| 2 | $\nabla\hat{\mathbf{m}}$ | Gradient of magnetization direction (default for `nspin=4`) | **Recommended** for most noncollinear calculations with SOC | -| 3 | Built-in SF | Special functional for noncollinear SOC (requires `v_xc_ncgga_sf_builtin`) | Advanced noncollinear SOC calculations | - -#### Recommendations - -**For noncollinear calculations with SOC (`lspinorb=1`)**: -- **Use `gga_grad=2`** (default when `nspin=4`): This computes the gradient of the magnetization direction unit vector $\hat{\mathbf{m}}$, which is physically correct for spin-orbit coupling effects. -- This is automatically set when you specify `nspin=4` or `lspinorb=1`, but you can explicitly set it to be clear. - -**For noncollinear calculations without SOC (`noncolin=1`, `lspinorb=0`)**: -- **Use `gga_grad=2`** for spin spiral calculations or frustrated magnetic systems -- The gradient of the magnetization direction is still important even without SOC - -**For collinear calculations (`nspin=1` or `nspin=2`)**: -- **Use `gga_grad=0`** (default): Standard collinear GGA gradient -- Do NOT use `gga_grad=1`, `2`, or `3` for collinear calculations - -#### Example INPUT for Noncollinear SOC - -``` -INPUT_PARAMETERS -calculation scf -basis_type pw -ecutwfc 50 -nspin 4 -noncolin 1 -lspinorb 1 -gga_grad 2 -``` - -#### Important Notes - -1. **Automatic setting**: When `nspin=4`, ABACUS automatically sets `gga_grad=2` if not explicitly specified -2. **Physical correctness**: Using the wrong `gga_grad` value can lead to incorrect magnetic anisotropy energies and spin textures -3. **Compatibility**: `gga_grad=3` requires the built-in special functional and may not be available for all XC functionals -4. **Performance**: `gga_grad=2` has similar computational cost to `gga_grad=0` for most systems - ## For the continuation job - Continuation job for "nspin 1" need file "SPIN1_CHG.cube" which is generated by setting "out_chg=1" in task before. By setting "init_chg file" in new job's INPUT file, charge density will start from file but not atomic. - Continuation job for "nspin 2" need files "SPIN1_CHG.cube" and "SPIN2_CHG.cube" which are generated by "out_chg 1" with "nspin 2", and refer to spin-up and spin-down charge densities respectively. It should be note that reading "SPIN1_CHG.cube" only for the continuation target magnetic moment job is not supported now. diff --git a/docs/parameters.yaml b/docs/parameters.yaml index 31f216e9dc..005fad3076 100644 --- a/docs/parameters.yaml +++ b/docs/parameters.yaml @@ -30,7 +30,7 @@ parameters: * md: perform molecular dynamics simulations * get_pchg: obtain partial (band-decomposed) charge densities (for LCAO basis only). See out_pchg for more information * get_wf: obtain real space wave functions (for LCAO basis only). See out_wfc_norm and out_wfc_re_im for more information - * get_s: obtain the overlap matrix formed by localized orbitals (for LCAO basis with multiple k points). the file name is SR.csr with file format being the same as that generated by out_mat_hs2 + * get_s: obtain the overlap matrix formed by localized orbitals (for LCAO basis with multiple k points). The file name is OUT.${suffix}/sr_nao.csr, with the same file format as generated by out_hsr 1 * gen_bessel: generates projectors, i.e., a series of Bessel functions, for the DeePKS method (for LCAO basis only) * gen_opt_abfs: generate opt-ABFs as discussed in this article * test_memory: obtain a rough estimation of memory consumption for the calculation @@ -63,7 +63,7 @@ parameters: Takes value 1, 0 or -1. * -1: No symmetry will be considered. It is recommended to set -1 for non-colinear + soc calculations, where time reversal symmetry is broken sometimes. * 0: Only time reversal symmetry would be considered in symmetry operations, which implied k point and -k point would be treated as a single k point with twice the weight. - * 1: Symmetry analysis will be performed to determine the type of Bravais lattice and associated symmetry operations. (point groups, space groups, primitive cells, and irreducible k-points) + * 1: Symmetry analysis will be performed to determine the type of Bravais lattice and associated symmetry operations (point groups, space groups, primitive cells, and irreducible k-points). For a magnetic system, the symmetry of the initial magnetic structure will be analyzed and preserved. [NOTE] When symmetry is enabled (value 1), k-points are reduced to the irreducible Brillouin zone (IBZ). For explicit k-point lists with custom weights (see KPT file), the custom weights are preserved during symmetry reduction. For Monkhorst-Pack grids, uniform weights are used. default_value: default @@ -533,9 +533,9 @@ parameters: For plane-wave basis, * cg: The conjugate-gradient (CG) method. - * bpcg: The BPCG method, which is a block-parallel Conjugate Gradient (CG) method, typically exhibits higher acceleration in a GPU environment. * dav: The Davidson algorithm. * dav_subspace: The Davidson algorithm without orthogonalization operation, this method is the most recommended for efficiency. `pw_diag_ndim` can be set to 2 for this method. + * bpcg: The BPCG method, which is a block-parallel Conjugate Gradient (CG) method, typically exhibits higher acceleration in a GPU environment. The BPCG method is currently under testing and is not recommended for use. For numerical atomic orbitals basis, @@ -621,20 +621,20 @@ parameters: category: Electronic structure type: Integer followed by Real values description: | - Customized parameterization on the exchange part of XC functional. The first value should be the LibXC ID of the original functional, and latter values are external parameters. Default values are those of Perdew-Burke-Ernzerhof (PBE) functional. For more information on LibXC ID of functionals, please refer to LibXC. For parameters of functionals of interest, please refer to the source code of LibXC, such as PBE functional interface in LibXC: gga_x_pbe.c. + Customized parameterization of the exchange part of an XC functional. The first value should be the Libxc ID of the original functional, followed by the complete list of external parameters required by the linked Libxc version. If unset, Libxc's own default parameters are used. For functional IDs and parameter definitions, refer to the Libxc documentation and source code. [NOTE] Solely setting this keyword will take no effect on XC functionals. One should also set dft_functional to the corresponding functional to apply the customized parameterization. Presently this feature can only support parameterization on one exchange functional. - default_value: 101 0.8040 0.2195149727645171 + default_value: "" unit: "" availability: "" - name: xc_corr_ext category: Electronic structure type: Integer followed by Real values description: | - Customized parameterization on the correlation part of XC functional. The first value should be the LibXC ID of the original functional, and latter values are external parameters. Default values are those of Perdew-Burke-Ernzerhof (PBE) functional. For more information on LibXC ID of functionals, please refer to LibXC. For parameters of functionals of interest, please refer to the source code of LibXC, such as PBE functional interface in LibXC: gga_c_pbe.c. + Customized parameterization of the correlation part of an XC functional. The first value should be the Libxc ID of the original functional, followed by the complete list of external parameters required by the linked Libxc version. If unset, Libxc's own default parameters are used. For functional IDs and parameter definitions, refer to the Libxc documentation and source code. [NOTE] Solely setting this keyword will take no effect on XC functionals. One should also set dft_functional to the corresponding functional to apply the customized parameterization. Presently this feature can only support parameterization on one correlation functional. - default_value: 130 0.06672455060314922 0.031090690869654895034 1.0 + default_value: "" unit: "" availability: "" - name: pseudo_rcut @@ -2953,7 +2953,7 @@ parameters: category: Output information type: "Boolean \\[Integer\\](optional)" description: | - Whether to output the density matrix with Bravias lattice vector R index into files in the folder OUT.${suffix}. The files are named as dmr{s}{spin index}{g}{geometry index}{_nao} + {".csr"}. Here, 's' refers to spin, where s1 means spin up channel while s2 means spin down channel, and the sparse matrix format 'csr' is mentioned in out_mat_hs2. Finally, if out_app_flag is set to false, the file name contains the optional 'g' index for each ionic step that may have different geometries, and if out_app_flag is set to true, the density matrix with respect to Bravias lattice vector R accumulates during ionic steps: + Whether to output the density matrix with Bravias lattice vector R index into files in the folder OUT.${suffix}. The files are named as dmr{s}{spin index}{g}{geometry index}{_nao} + {".csr"}. Here, 's' refers to spin, where s1 means spin up channel while s2 means spin down channel, and the sparse matrix format 'csr' is mentioned in out_hsr. Finally, if out_app_flag is set to false, the file name contains the optional 'g' index for each ionic step that may have different geometries, and if out_app_flag is set to true, the density matrix with respect to Bravias lattice vector R accumulates during ionic steps: * nspin = 1: dmrs1_nao.csr; * nspin = 2: dmrs1_nao.csr and dmrs2_nao.csr for the two spin channels. @@ -3065,34 +3065,61 @@ parameters: default_value: ie unit: "" availability: "" - - name: out_mat_hs + - name: out_hsk category: Output information - type: "Boolean \\[Integer\\](optional)" + type: "Integer \\[Integer\\](optional)" description: | - Whether to print the upper triangular part of the Hamiltonian matrices and overlap matrices for each k-point into files in the directory OUT.${suffix}. The second number controls precision. For more information, please refer to hs_matrix.md. Also controled by out_freq_ion and out_app_flag. - * For gamma only case: - * nspin = 1: hks1_nao.txt for the Hamiltonian matrix and sks1_nao.txt for the overlap matrix; - * nspin = 2: hks1_nao.txt and hks2_nao.txt for the Hamiltonian matrix and sks1_nao.txt for the overlap matrix. Note that the code will not output sks2_nao.txt because it is the same as sks1_nao.txt; - * nspin = 4: hks12_nao.txt for the Hamiltonian matrix and sks12_nao.txt for the overlap matrix. - * For multi-k points case: - * nspin = 1: hks1k1_nao.txt for the Hamiltonian matrix at the 1st k-point, and sks1k1_nao.txt for the overlap matrix for the 1st k-point, ...; - * nspin = 2: hks1k1_nao.txt and hks2k1_nao.txt for the two spin channels of the Hamiltonian matrix at the 1st k-point, and sks1k1_nao.txt for the overlap matrix for the 1st k-point. Note that the code will not output sks2k1_nao.txt because it is the same as sks1k1_nao.txt, ...; - * nspin = 4: hks12k1_nao.txt for the Hamiltonian matrix at the 1st k-point, and sks12k1_nao.txt for the overlap matrix for the 1st k-point, ...; + Output the upper triangular part of the Hamiltonian and overlap matrices in reciprocal space for each k-point into files in the directory OUT.${suffix}. The first integer selects the format: + * 0: disabled; + * 1: text output; the optional second integer controls precision and defaults to 8; + * 2: reserved for binary output, which is not implemented yet; + * 3: NPZ output, which is not implemented for H(k)/S(k). + + The output is also controlled by out_freq_ion and out_app_flag. For more information, refer to hs_matrix.md. + * Gamma-only, nspin = 1: hk_nao.txt for the Hamiltonian matrix and sk_nao.txt for the overlap matrix. + * Gamma-only, nspin = 2: hks1_nao.txt and hks2_nao.txt for the two spin channels of the Hamiltonian matrix, and sk_nao.txt for the overlap matrix. Only one overlap matrix is written because it is identical for both spin channels. + * Gamma-only, nspin = 4: not available with the gamma-only algorithm. + * Multi-k, nspin = 1: hk1_nao.txt for the Hamiltonian matrix and sk1_nao.txt for the overlap matrix at the first k-point. + * Multi-k, nspin = 2: hk1s1_nao.txt and hk1s2_nao.txt for the two spin channels of the Hamiltonian matrix, and sk1_nao.txt for the overlap matrix at the first k-point. Only one overlap matrix is written because it is identical for both spin channels. + * Multi-k, nspin = 4: hk1s4_nao.txt for the spinor Hamiltonian matrix and sk1_nao.txt for the spinor overlap matrix at the first k-point. + When out_app_flag is false, g followed by the one-based ionic-step index is inserted before _nao, for example hk1s1g1_nao.txt. [NOTE] In the 3.10-LTS version, the file names are data-0-H and data-0-S, etc. - default_value: False 8 + default_value: 0 8 unit: Ry availability: Numerical atomic orbital basis - - name: out_mat_hs2 + - name: out_mat_hs category: Output information type: "Boolean \\[Integer\\](optional)" description: | - Whether to print files containing the Hamiltonian matrix and overlap matrix into files in the directory OUT.${suffix}. For more information, please refer to hs_matrix.md. + Legacy alias for out_hsk 1, which outputs Hamiltonian and overlap matrices in reciprocal space for each k-point. The optional second integer controls text precision. If both out_hsk and out_mat_hs are present, out_hsk takes precedence. + default_value: False 8 + unit: Ry + availability: Numerical atomic orbital basis + - name: out_hsr + category: Output information + type: "Integer \\[Integer\\](optional)" + description: | + Output Hamiltonian and overlap matrices in real space, indexed by the Bravais lattice vector R, in the directory OUT.${suffix}. The first integer selects the format: + * 0: disabled; + * 1: text CSR output; the optional second integer controls precision and defaults to 8; + * 2: reserved for binary output, which is not implemented yet; + * 3: NPZ output using hrs1_nao.npz, hrs2_nao.npz when needed, and sr_nao.npz. + + For multi-k calculations, the output contains the individual real-space blocks stored for the Bravais lattice vectors R. For gamma-only calculations, the internal real-space contributions are folded into a single R = (0, 0, 0) block. This folded result cannot recover the original R-resolved contributions or interpolate arbitrary k points. Terms added only while constructing H(k) are not guaranteed to be present. [NOTE] In the 3.10-LTS version, the file names are data-HR-sparse_SPIN0.csr and data-SR-sparse_SPIN0.csr, etc. - default_value: "False [8]" + default_value: 0 8 unit: Ry - availability: Numerical atomic orbital basis (not gamma-only algorithm) + availability: Numerical atomic orbital basis + - name: out_mat_hs2 + category: Output information + type: "Boolean \\[Integer\\](optional)" + description: | + Legacy alias for out_hsr 1, which outputs Hamiltonian and overlap matrices in real space indexed by the Bravais lattice vector R. The optional second integer controls text precision. If both out_hsr and out_mat_hs2 are present, out_hsr takes precedence. + default_value: False 8 + unit: Ry + availability: Numerical atomic orbital basis - name: out_mat_tk category: Output information type: "Boolean \\[Integer\\](optional)" @@ -3117,7 +3144,7 @@ parameters: category: Output information type: "Boolean \\[Integer\\](optional)" description: | - Generate files containing the kinetic energy matrix. The optional second parameter controls text output precision. The format will be the same as the Hamiltonian matrix and overlap matrix as mentioned in out_mat_hs2. The name of the files will be trs1_nao.csr and so on. Also controled by out_freq_ion and out_app_flag. + Generate files containing the kinetic energy matrix. The optional second parameter controls text output precision. The format will be the same as the Hamiltonian matrix and overlap matrix as mentioned in out_hsr. The name of the files will be trs1_nao.csr and so on. Also controled by out_freq_ion and out_app_flag. [NOTE] In the 3.10-LTS version, the file name is data-TR-sparse_SPIN0.csr. default_value: False 8 @@ -3127,7 +3154,7 @@ parameters: category: Output information type: Integer description: | - Whether to print files containing the derivatives of the Hamiltonian matrix. The format will be the same as the Hamiltonian matrix and overlap matrix as mentioned in out_mat_hs2. The name of the files will be dhrxs1_nao.csr, dhrys1_nao.csr, dhrzs1_nao.csr and so on. Also controled by out_freq_ion and out_app_flag. + Whether to print files containing the derivatives of the Hamiltonian matrix. The format will be the same as the Hamiltonian matrix and overlap matrix as mentioned in out_hsr. The name of the files will be dhrxs1_nao.csr, dhrys1_nao.csr, dhrzs1_nao.csr and so on. Also controled by out_freq_ion and out_app_flag. Format: [precision] [iat1 iat2 ...]. The first value (0/1) enables/disables output. The second optional value sets the output precision (default: 8). Starting from the third value, 1-based atom indices can be listed to restrict output to derivatives with respect to those specific atoms only; if no atom indices are given, all atoms are written. @@ -3201,7 +3228,7 @@ parameters: description: | Whether to print files containing the kinetic energy matrix T(R) in CSR format. - See out_mat_hs2 for format details. + See out_hsr for format details. default_value: 0 8 unit: Ry availability: "" @@ -3211,7 +3238,7 @@ parameters: description: | Whether to print files containing the nonlocal pseudopotential matrix Vnl(R) in CSR format. - See out_mat_hs2 for format details. + See out_hsr for format details. default_value: 0 8 unit: Ry availability: "" @@ -3221,7 +3248,7 @@ parameters: description: | Whether to print files containing the local pseudopotential matrix Vl(R) in CSR format. - See out_mat_hs2 for format details. + See out_hsr for format details. default_value: 0 8 unit: Ry availability: "" @@ -3231,7 +3258,7 @@ parameters: description: | Whether to print files containing the Hartree matrix Vh(R) in CSR format. - See out_mat_hs2 for format details. + See out_hsr for format details. default_value: 0 8 unit: Ry availability: "" @@ -3241,7 +3268,7 @@ parameters: description: | Whether to print files containing the XC matrix Vxc(R) in CSR format. - See out_mat_hs2 for format details. + See out_hsr for format details. default_value: 0 8 unit: Ry availability: "" @@ -3251,7 +3278,7 @@ parameters: description: | Whether to print files containing the exact-exchange matrix Vexx(R) in CSR format. - See out_mat_hs2 for format details. + See out_hsr for format details. default_value: 0 8 unit: Ry availability: "" @@ -3322,7 +3349,7 @@ parameters: category: Output information type: Boolean description: | - Whether to print Hamiltonian matrices H(R) in npz format. This feature does not work for gamma-only calculations. + Whether to print Hamiltonian matrices H(R) in NPZ format as hrs1_nao.npz and, for nspin = 2, hrs2_nao.npz. This feature does not work for gamma-only calculations. default_value: "False" unit: Ry availability: Numerical atomic orbital basis (not gamma-only algorithm) @@ -3330,10 +3357,10 @@ parameters: category: Output information type: Boolean description: | - Whether to print Hamiltonian matrices H(R) and overlap matrix S(R) in npz format. This feature does not work for gamma-only calculations. + Legacy alias for out_hsr 3, writing hrs1_nao.npz, hrs2_nao.npz when needed, and sr_nao.npz. If both out_hsr and out_hsr_npz are present, out_hsr takes precedence. Gamma-only calculations write the folded R = (0, 0, 0) representation. default_value: "False" unit: Ry - availability: Numerical atomic orbital basis (not gamma-only algorithm) + availability: Numerical atomic orbital basis - name: out_dm_npz category: Output information type: Boolean @@ -3354,7 +3381,7 @@ parameters: category: Output information type: Boolean description: | - Whether to output r(R), H(R), S(R), T(R), dH(R), dS(R), and wfc matrices in an append manner during molecular dynamics calculations. Check input parameters out_mat_r, out_mat_hs2, out_mat_t, out_mat_dh, out_mat_hs and out_wfc_lcao for more information. + Whether to output r(R), H(R), S(R), T(R), dH(R), dS(R), and wfc matrices in an append manner during molecular dynamics calculations. Check input parameters out_mat_r, out_hsr, out_mat_t, out_mat_dh, out_hsk and out_wfc_lcao for more information. default_value: "true" unit: "" availability: Numerical atomic orbital basis (not gamma-only algorithm) @@ -3365,7 +3392,7 @@ parameters: Controls the length of decimal part of output data, such as charge density, Hamiltonian matrix, Overlap matrix and so on. default_value: "8" unit: "" - availability: out_mat_hs 1 case presently. + availability: out_hsk 1 case presently. - name: out_element_info category: Output information type: Boolean @@ -4473,45 +4500,8 @@ parameters: category: Spin-Constrained DFT type: Real description: | - When the charge density error drho falls below sc_scf_thr, the DeltaSpin lambda optimization loop is activated. Should be 10-100x larger than scf_thr. Only used when sc_scf_thr_mode="threshold". For other activation modes, use sc_scf_thr_mode="immediate" or sc_scf_thr_mode="off". - default_value: "1.0e-3" - unit: "" - availability: sc_mag_switch is true - - name: sc_scf_thr_mode - category: Spin-Constrained DFT - type: String - description: | - Controls when the DeltaSpin lambda loop is activated. - * threshold (default): activate when drho < sc_scf_thr. The lambda loop starts once the charge density is reasonably stable. - * immediate: activate from the first iteration with valid wavefunctions (iter>=2). Used for PW basis where the first iteration cannot compute initial magnetic moments. Replaces the old convention of setting sc_scf_thr=10.0. - * off: never activate the lambda loop. Lambda values are loaded from the STRU file and used as constant constraints without optimization. Replaces the old convention of setting sc_scf_thr=1e-10. - default_value: "threshold" - unit: "" - availability: sc_mag_switch is true - - name: sc_direction_only - category: Spin-Constrained DFT - type: Boolean - description: | - When true, only the direction of the magnetic moment is constrained to the target direction, while the magnitude is allowed to vary freely. This is useful for studying magnetic anisotropy or when the magnitude of the moment is determined by the electronic structure rather than an external constraint. When false (default), both the direction and magnitude of the magnetic moment are constrained to the target values. - default_value: "False" - unit: "" - availability: sc_mag_switch is true - - name: sc_dir_phase1_steps - category: Spin-Constrained DFT - type: Integer - description: | - Number of SCF iterations for Phase 1 (magnitude constraint) in the direction_only two-phase strategy for collinear (nspin=2) calculations. During Phase 1, direction_only projection is temporarily disabled so BFGS can constrain the magnetic moment magnitude. After Phase 1, lambda decays and the system relaxes naturally. Minimum: 2. - default_value: "5" - unit: "" - availability: sc_mag_switch is true and sc_direction_only is true and nspin=2 - - name: sc_lambda_strategy - category: Spin-Constrained DFT - type: String - description: | - Lambda update strategy for spin-constrained DFT. Available options are: - * bfgs: BFGS quasi-Newton method (default, robust and well-tested) - * linear_scan: linear sweep of lambda for testing magnetic moment response - default_value: "bfgs" + Density error threshold for inner loop of spin-constrained SCF + default_value: "1.0e-4" unit: "" availability: sc_mag_switch is true - name: sc_direction_only diff --git a/source/source_hamilt/module_xc/libxc_setup.cpp b/source/source_hamilt/module_xc/libxc_setup.cpp index 3b29edf727..9a2a4c0e17 100644 --- a/source/source_hamilt/module_xc/libxc_setup.cpp +++ b/source/source_hamilt/module_xc/libxc_setup.cpp @@ -14,7 +14,6 @@ #include #include #include -#include bool not_supported_xc_with_laplacian(const std::string& xc_func_in) { @@ -200,7 +199,7 @@ const std::vector in_built_xc_func_ext_params(const int id, #ifdef __EXX // hybrid functionals case XC_HYB_GGA_XC_PBEH: - return {hybrid_alpha, hse_omega, hse_omega}; + return {hybrid_alpha}; case XC_HYB_GGA_XC_HSE06: return {hybrid_alpha, hse_omega, hse_omega}; // short-range of B88_X @@ -262,20 +261,19 @@ const std::vector in_built_xc_func_ext_params(const int id, const std::vector external_xc_func_ext_params(const int id) { - const std::map> mymap = { - { - PARAM.inp.xc_exch_ext[0], - std::vector(PARAM.inp.xc_exch_ext.begin()+1, - PARAM.inp.xc_exch_ext.end()) - }, - { - PARAM.inp.xc_corr_ext[0], - std::vector(PARAM.inp.xc_corr_ext.begin()+1, - PARAM.inp.xc_corr_ext.end()) - } - }; - auto it = mymap.find(id); - return (it != mymap.end()) ? it->second : std::vector{}; + const auto& exch_ext = PARAM.inp.xc_exch_ext; + if (!exch_ext.empty() && static_cast(exch_ext.front()) == id) + { + return {exch_ext.begin() + 1, exch_ext.end()}; + } + + const auto& corr_ext = PARAM.inp.xc_corr_ext; + if (!corr_ext.empty() && static_cast(corr_ext.front()) == id) + { + return {corr_ext.begin() + 1, corr_ext.end()}; + } + + return {}; } std::vector @@ -290,33 +288,32 @@ XC_Functional_Libxc::init_func(const std::vector &func_id, funcs.push_back({}); // create placeholder xc_func_init(&funcs.back(), id, xc_polarized); // instantiate the XC term - // search for external parameters - const std::vector in_built_ext_params = in_built_xc_func_ext_params(id, hybrid_alpha, hse_omega); - const std::vector external_ext_params = external_xc_func_ext_params(id); - // for temporary use, I name their size as n1 and n2 - const int n1 = in_built_ext_params.size(); - const int n2 = external_ext_params.size(); + // Search for external parameters. User-supplied parameters take precedence + // over ABACUS built-in overrides. + const std::vector in_built_ext_params + = in_built_xc_func_ext_params(id, hybrid_alpha, hse_omega); + const std::vector external_ext_params + = external_xc_func_ext_params(id); + const std::vector& requested_ext_params + = external_ext_params.empty() ? in_built_ext_params : external_ext_params; -// #ifdef __DEBUG // will the following assertion cause performance issue? - // assert the number of parameters should be either zero or the value from - // libxc function xc_func_info_get_n_ext_params, this is to avoid the undefined - // behavior of illegal memory access - const xc_func_info_type* info = xc_func_get_info(&funcs.back()); - const int nref = xc_func_info_get_n_ext_params(info); - assert ((n1 == 0) || (n1 == nref) || (n2 == 0) || (n2 == nref)); -// #endif + if (!requested_ext_params.empty()) + { + const xc_func_info_type* info = xc_func_get_info(&funcs.back()); + const int nref = xc_func_info_get_n_ext_params(info); - // external overwrites in-built if the same functional id is found in both maps - const double* xc_func_ext_params = - (n2 > 0) ? external_ext_params.data() : - (n1 > 0) ? in_built_ext_params.data() : - nullptr; // nullptr if no external parameters are found + // xc_func_set_ext_params() reads exactly nref entries + if (requested_ext_params.size() != static_cast(nref)) + { + ModuleBase::WARNING_QUIT( + "XC_Functional_Libxc::init_func", + "Invalid number of external parameters for Libxc functional id " + + std::to_string(id) + ": got " + + std::to_string(requested_ext_params.size()) + + ", expected " + std::to_string(nref) + "."); + } - // if there are no external parameters, do nothing, otherwise we set - if(xc_func_ext_params != nullptr) - { - // set the external parameters - xc_func_set_ext_params(&funcs.back(), const_cast(xc_func_ext_params)); + xc_func_set_ext_params(&funcs.back(), requested_ext_params.data()); } } return funcs; diff --git a/source/source_hamilt/module_xc/libxc_tools.cpp b/source/source_hamilt/module_xc/libxc_tools.cpp index e3694044f0..1af2d61bcc 100644 --- a/source/source_hamilt/module_xc/libxc_tools.cpp +++ b/source/source_hamilt/module_xc/libxc_tools.cpp @@ -250,8 +250,8 @@ std::pair XC_Functional_Libxc::convert_vtxc_v( for (std::size_t ir = 0; ir < nrxx; ++ir) { double rho_ir = 0.0; - if (is == 0) { rho_ir = rho[ir * nspin + 0]; } - else { rho_ir = rho[ir * nspin + 0] * mag_part_tmp[ir + (is - 1) * nrxx]; } + if (is == 0) { rho_ir = chr->rho[0][ir]; } + else { rho_ir = chr->rho[is][ir]; } rvtxc += dh[is][ir] * rho_ir; v(is, ir) -= dh[is][ir]; } diff --git a/source/source_io/module_parameter/input_parameter.h b/source/source_io/module_parameter/input_parameter.h index 44cc83a15e..2d9b8d42a4 100644 --- a/source/source_io/module_parameter/input_parameter.h +++ b/source/source_io/module_parameter/input_parameter.h @@ -69,7 +69,7 @@ struct Input_para std::string kmesh_type = "gamma"; ///< k-point mesh type for kspacing-generated k-point mesh: gamma or mp double min_dist_coef = 0.2; ///< allowed minimum distance between two atoms - std::string device = "auto"; + std::string device = "cpu"; std::string precision = "double"; std::string gint_precision = "double"; bool timer_enable_nvtx = false; @@ -88,7 +88,6 @@ struct Input_para ///< use mesh, which is used in QE. int nspin = 1; ///< LDA ; LSDA ; non-linear spin int gga_grad = 3; ///< GGA gradient method for noncollinear spin (nspin=4): 1=collinear approx, 2=projected div(h), 3=Scalmani-Frisch transform (default, most accurate) - std::string device_memory_mode = ""; ///< GPU memory mode: "" (auto), "full_gpu" (all on GPU), "paged" (CPU storage + k-point paging) int pw_diag_nmax = 50; double pw_diag_thr = 0.01; ///< used in cg method bool diago_smooth_ethr = false; ///< smooth ethr for iter methods @@ -619,23 +618,18 @@ struct Input_para bool decay_grad_switch = false; ///< the switch to use the local approximation of gradient ///< decay, 0: no local approximation; 1: apply the method double sc_thr = 1e-06; ///< threshold for spin-constrained DFT in uB - int nsc = 5; ///< maximum number of inner lambda loop + int nsc = 100; ///< maximum number of inner lambda loop int nsc_min = 2; ///< minimum number of inner lambda loop double alpha_trial = 0.01; ///< initial trial step size for lambda in eV/uB^2 double sccut = 3.0; ///< restriction of step size in eV/uB - double sc_scf_thr = 10; ///< density error threshold for activating the lambda loop in spin-constrained DFT - std::string sc_scf_thr_mode = "immediate"; ///< controls when the lambda loop activates: "threshold" (drho=2), or "off" (never activate) - int sc_dir_phase1_steps = 5; ///< number of Phase 1 iterations in direction_only two-phase strategy for collinear (nspin=2) + double sc_scf_thr = 1e-3; ///< minimum number of outer scf loop before initial lambda loop double sc_drop_thr = 1e-3; ///< threshold for lambda-loop threshold cutoff in spin-constrained DFT - std::string sc_lambda_strategy = "bfgs"; ///< lambda update strategy: bfgs or linear_scan + std::string sc_lambda_strategy = "bfgs"; ///< lambda update strategy: bfgs, bfgs2, linear_response, augmented_lagrangian, hybrid_delayed, linear_scan bool sc_direction_only = false; ///< only optimize the direction of magnetization // linear_scan parameters double sc_scan_lambda_start = 0.0; ///< start value for lambda scan (eV/uB) double sc_scan_lambda_end = 1.0; ///< end value for lambda scan (eV/uB) int sc_scan_steps = 20; ///< number of steps in lambda scan - std::string sc_strategy = "normal"; ///< DeltaSpin execution strategy: "fast" (subspace), "accuracy" (full), "normal" (threshold-triggered) - std::string sc_acceleration_mode = "off"; ///< acceleration mode: "off", "first_order", "subspace" - double sc_acceleration_rms_thr = -1.0; ///< RMS threshold (uB) to activate acceleration, <0 disables // ============== #Parameters (18.Quasiatomic Orbital analysis) ========= ///<========================================================== @@ -732,14 +726,12 @@ struct Input_para * * Likewise, the correlation part can be found in corresponding files. * - * PBE functional is used as the default functional for XCPNet. + * These vectors are empty unless the user explicitly requests an + * override. This leaves the version-specific default parameters under + * Libxc's control. */ - // src/gga_x_pbe.c - std::vector xc_exch_ext = { - 101, 0.8040, 0.2195149727645171}; - // src/gga_c_pbe.c - std::vector xc_corr_ext = { - 130, 0.06672455060314922, 0.031090690869654895034, 1.00000}; + std::vector xc_exch_ext = {}; + std::vector xc_corr_ext = {}; // ============== #Parameters (24.td-ofdft) =========================== bool of_cd = false; ///< add CD potential or not diff --git a/source/source_io/module_parameter/read_input_item_elec_stru.cpp b/source/source_io/module_parameter/read_input_item_elec_stru.cpp index 87f290f21f..5ddaeb8bd1 100644 --- a/source/source_io/module_parameter/read_input_item_elec_stru.cpp +++ b/source/source_io/module_parameter/read_input_item_elec_stru.cpp @@ -53,9 +53,9 @@ void ReadInput::item_elec_stru() For plane-wave basis, * cg: The conjugate-gradient (CG) method. -* bpcg: The BPCG method, which is a block-parallel Conjugate Gradient (CG) method, typically exhibits higher acceleration in a GPU environment. * dav: The Davidson algorithm. * dav_subspace: The Davidson algorithm without orthogonalization operation, this method is the most recommended for efficiency. `pw_diag_ndim` can be set to 2 for this method. +* bpcg: The BPCG method, which is a block-parallel Conjugate Gradient (CG) method, typically exhibits higher acceleration in a GPU environment. The BPCG method is currently under testing and is not recommended for use. For numerical atomic orbitals basis, @@ -340,28 +340,6 @@ The other way is only available when compiling with LIBXC, and it allows for sup read_sync_string(input.dft_functional); this->add_item(item); } - { - Input_Item item("gga_grad"); - item.annotation = "1: collinear approximation (grad|m|); 2: projected div(h); 3: Scalmani-Frisch transform (default, most accurate)"; - item.category = "Electronic structure"; - item.type = "Integer"; - item.description = R"(Method for computing GGA gradients in noncollinear spin (nspin=4) calculations. -* 1: Collinear approximation. Treats nspin=4 as nspin=2 with |m|, computes grad(rho +/- |m|). Fastest but least accurate. -* 2: Improved gradient. Computes grad(m_mu)*m_hat_mu and projects div(h) via m_hat. Intermediate accuracy. -* 3: Scalmani-Frisch transform. Applies SF rotation to decompose gradients per spin channel and retains all cross-terms in div(h). Most accurate. -Only relevant when nspin=4 and a GGA or hybrid functional is used.)"; - item.default_value = "3"; - item.unit = ""; - item.availability = "Only relevant for nspin=4 with GGA/hybrid functionals"; - read_sync_int(input.gga_grad); - item.check_value = [](const Input_Item&, const Parameter& para) { - if (para.input.gga_grad < 1 || para.input.gga_grad > 3) - { - ModuleBase::WARNING_QUIT("ReadInput", "gga_grad must be 1, 2, or 3."); - } - }; - this->add_item(item); - } { Input_Item item("xc_temperature"); item.annotation = "temperature for finite temperature functionals"; @@ -376,14 +354,14 @@ Only relevant when nspin=4 and a GGA or hybrid functional is used.)"; } { Input_Item item("xc_exch_ext"); - item.annotation = "placeholder for xcpnet exchange functional"; + item.annotation = "customize Libxc exchange functional parameters"; item.category = "Electronic structure"; item.type = "Integer followed by Real values"; - item.description = "Customized parameterization on the exchange part of XC functional. The first value should be the LibXC ID of the original functional, and latter values are external parameters. Default values are those of Perdew-Burke-Ernzerhof (PBE) functional. For more information on LibXC ID of functionals, please refer to LibXC. For parameters of functionals of interest, please refer to the source code of LibXC, such as PBE functional interface in LibXC: gga_x_pbe.c." + item.description = "Customized parameterization of the exchange part of an XC functional. The first value should be the Libxc ID of the original functional, followed by the complete list of external parameters required by the linked Libxc version. If unset, Libxc's own default parameters are used. For functional IDs and parameter definitions, refer to the Libxc documentation and source code." "\n\n[NOTE] Solely setting this keyword will take no effect on XC functionals. One should also set " "dft_functional to the corresponding functional to apply the customized parameterization. " "Presently this feature can only support parameterization on one exchange functional."; - item.default_value = "101 0.8040 0.2195149727645171"; + item.default_value = ""; item.unit = ""; item.availability = ""; item.read_value = [](const Input_Item& item, Parameter& para) { @@ -393,10 +371,15 @@ Only relevant when nspin=4 and a GGA or hybrid functional is used.)"; [](const std::string& str) { return std::stod(str); }); }; item.check_value = [](const Input_Item& item, const Parameter& para) { - // at least one value should be set - if (para.input.xc_exch_ext.empty()) + if (!item.is_read()) { - ModuleBase::WARNING_QUIT("ReadInput", "xc_exch_ext should not be empty."); + return; + } + if (para.input.xc_exch_ext.size() < 2) + { + ModuleBase::WARNING_QUIT( + "ReadInput", + "xc_exch_ext requires a Libxc ID followed by external parameters."); } // the first value is actually an integer, not a double const double libxc_id_dbl = para.input.xc_exch_ext[0]; @@ -419,14 +402,14 @@ Only relevant when nspin=4 and a GGA or hybrid functional is used.)"; } { Input_Item item("xc_corr_ext"); - item.annotation = "placeholder for xcpnet exchange functional"; + item.annotation = "customize Libxc correlation functional parameters"; item.category = "Electronic structure"; item.type = "Integer followed by Real values"; - item.description = "Customized parameterization on the correlation part of XC functional. The first value should be the LibXC ID of the original functional, and latter values are external parameters. Default values are those of Perdew-Burke-Ernzerhof (PBE) functional. For more information on LibXC ID of functionals, please refer to LibXC. For parameters of functionals of interest, please refer to the source code of LibXC, such as PBE functional interface in LibXC: gga_c_pbe.c." + item.description = "Customized parameterization of the correlation part of an XC functional. The first value should be the Libxc ID of the original functional, followed by the complete list of external parameters required by the linked Libxc version. If unset, Libxc's own default parameters are used. For functional IDs and parameter definitions, refer to the Libxc documentation and source code." "\n\n[NOTE] Solely setting this keyword will take no effect on XC functionals. One should also set " "dft_functional to the corresponding functional to apply the customized parameterization. " "Presently this feature can only support parameterization on one correlation functional."; - item.default_value = "130 0.06672455060314922 0.031090690869654895034 1.0"; + item.default_value = ""; item.unit = ""; item.availability = ""; item.read_value = [](const Input_Item& item, Parameter& para) { @@ -436,10 +419,15 @@ Only relevant when nspin=4 and a GGA or hybrid functional is used.)"; [](const std::string& str) { return std::stod(str); }); }; item.check_value = [](const Input_Item& item, const Parameter& para) { - // at least one value should be set - if (para.input.xc_corr_ext.empty()) + if (!item.is_read()) + { + return; + } + if (para.input.xc_corr_ext.size() < 2) { - ModuleBase::WARNING_QUIT("ReadInput", "xc_corr_ext should not be empty."); + ModuleBase::WARNING_QUIT( + "ReadInput", + "xc_corr_ext requires a Libxc ID followed by external parameters."); } // the first value is actually an integer, not a double const double libxc_id_dbl = para.input.xc_corr_ext[0]; @@ -513,6 +501,21 @@ Only relevant when nspin=4 and a GGA or hybrid functional is used.)"; }; this->add_item(item); } + { + Input_Item item("gga_grad"); + item.annotation = "GGA gradient method for nspin=4: 1 collinear approx, 2 projected div(h), 3 Scalmani-Frisch"; + item.category = "Electronic structure"; + item.type = "Integer"; + item.default_value = "3"; + read_sync_int(input.gga_grad); + item.check_value = [](const Input_Item& item, const Parameter& para) { + if (para.input.gga_grad < 1 || para.input.gga_grad > 3) + { + ModuleBase::WARNING_QUIT("ReadInput", "gga_grad must be 1, 2, or 3."); + } + }; + this->add_item(item); + } { Input_Item item("smearing_method"); item.annotation = "type of smearing_method: gauss; fd; fixed; mp; mp2; mv"; @@ -687,24 +690,16 @@ For systems that are difficult to converge, one could try increasing the value o item.availability = ""; read_sync_double(input.mixing_restart); item.reset_value = [](const Input_Item& item, Parameter& para) { - if (para.input.sc_mag_switch) - { - if (para.input.sc_direction_only) - { - para.input.mixing_restart = 0.0; - } - else if (para.input.sc_scf_thr_mode == "threshold") + if (para.input.sc_mag_switch == 1) + {// for DeltaSpin calculation, the mixing_restart should be same as sc_scf_thr + if(para.input.sc_scf_thr != 10.0) { para.input.mixing_restart = para.input.sc_scf_thr; } - else if (para.input.sc_scf_thr_mode == "immediate") - { + else + {// no mixing_restart until oscillation happen in PW base para.input.mixing_restart = para.input.scf_thr / 10.0; } - else // "off" - { - para.input.mixing_restart = 0.0; - } } }; this->add_item(item); From 2c234bf7016c1d0470ee857b81ab40574d0d0297 Mon Sep 17 00:00:00 2001 From: dyzheng Date: Sun, 2 Aug 2026 12:30:46 +0800 Subject: [PATCH 3/3] docs(gga_grad): add comments explaining SF formulas and branching logic - libxc_pot.cpp: comments for SF gradient decomposition, gga_grad branching - libxc_tools.cpp: comment for nspin=4 SF dh path with 4-component potential - xc_pot.cpp: comment for non-libxc SF builtin dispatch - xc_functional_gga_noncol_sf_builtin.cpp: header comment block with Scalmani-Frisch reference, chain-rule formulas, and step-by-step labels --- source/source_hamilt/module_xc/libxc_pot.cpp | 15 +++++++++++ .../source_hamilt/module_xc/libxc_tools.cpp | 5 ++++ .../xc_functional_gga_noncol_sf_builtin.cpp | 27 +++++++++++++++++++ source/source_hamilt/module_xc/xc_pot.cpp | 4 +++ 4 files changed, 51 insertions(+) diff --git a/source/source_hamilt/module_xc/libxc_pot.cpp b/source/source_hamilt/module_xc/libxc_pot.cpp index a6fc73d9ed..a035edb1ec 100644 --- a/source/source_hamilt/module_xc/libxc_pot.cpp +++ b/source/source_hamilt/module_xc/libxc_pot.cpp @@ -64,6 +64,10 @@ std::tuple XC_Functional_Libxc::v_xc_libxc( / }(); // converting rho + // For nspin=4, the charge density has 4 components: + // rho[0] = total charge, rho[1..3] = magnetization (mx, my, mz) + // libxc works with spin-up/spin-down densities: + // rho_up = 0.5*(rho[0] + |m|), rho_dn = 0.5*(rho[0] - |m|) std::vector rho; std::vector amag; std::vector mag_part; @@ -76,6 +80,11 @@ std::tuple XC_Functional_Libxc::v_xc_libxc( / std::tuple,std::vector> rho_amag = XC_Functional_Libxc::convert_rho_amag_nspin4(nspin, nrxx, chr); rho = std::get<0>(std::move(rho_amag)); amag = std::get<1>(std::move(rho_amag)); + // gga_grad=2/3: compute magnetization unit vector m_hat = m/|m| + // needed for the Scalmani-Frisch (SF) gradient decomposition: + // grad(rho_up) = 0.5 grad(rho[0]) + 0.5 m_hat_mu * grad(m_mu) + // grad(rho_dn) = 0.5 grad(rho[0]) - 0.5 m_hat_mu * grad(m_mu) + // gga_grad=1 (else): collinear approximation, uses grad(|m|) only if((PARAM.inp.gga_grad == 2 || PARAM.inp.gga_grad == 3) && (PARAM.globalv.domag || PARAM.globalv.domag_z)) { mag_part = XC_Functional_Libxc::compute_mag_part_nspin4(nrxx, chr); @@ -86,6 +95,9 @@ std::tuple XC_Functional_Libxc::v_xc_libxc( / std::vector sigma; if(is_gga) { + // gga_grad=2/3: use SF method to compute spin-up/spin-down gradients + // via the chain-rule decomposition. This is more accurate than the + // collinear approximation (gga_grad=1) for noncollinear magnetism. if(PARAM.inp.nspin==4 && (PARAM.inp.gga_grad == 2 || PARAM.inp.gga_grad == 3) && (PARAM.globalv.domag || PARAM.globalv.domag_z)) { gdr = XC_Functional_Libxc::cal_gdr_sf(nspin, nrxx, rho, mag_part, tpiba, chr); @@ -191,6 +203,9 @@ std::tuple XC_Functional_Libxc::v_xc_libxc( / if(4==nspin_in) { + // gga_grad=2/3: convert libxc spin-up/spin-down potential back to + // nspin=4 representation via SF method (v_total, v_mag_hat) + // gga_grad=1: standard conversion using |m| decomposition only if((PARAM.inp.gga_grad == 2 || PARAM.inp.gga_grad == 3) && (PARAM.globalv.domag || PARAM.globalv.domag_z)) { v = XC_Functional_Libxc::convert_v_nspin4_sf(nrxx, chr, mag_part, v); diff --git a/source/source_hamilt/module_xc/libxc_tools.cpp b/source/source_hamilt/module_xc/libxc_tools.cpp index 1af2d61bcc..6b4395e244 100644 --- a/source/source_hamilt/module_xc/libxc_tools.cpp +++ b/source/source_hamilt/module_xc/libxc_tools.cpp @@ -235,6 +235,11 @@ std::pair XC_Functional_Libxc::convert_vtxc_v( if(func.info->family == XC_FAMILY_GGA || func.info->family == XC_FAMILY_HYB_GGA) { + // For nspin=4 with domag and gga_grad=2/3, use the SF (Scalmani-Frisch) + // gradient correction. The potential v has 4 components: + // v(0,:) = total density channel, v(1..3,:) = magnetization channels + // cal_dh_sf returns dh[4][nrxx] in the same 4-component representation. + // The vtxc contribution sums dh[is] * rho[is] for each channel. if(PARAM.inp.nspin==4 && (PARAM.inp.gga_grad == 2 || PARAM.inp.gga_grad == 3) && (PARAM.globalv.domag || PARAM.globalv.domag_z)) { std::vector mag_part_tmp = XC_Functional_Libxc::compute_mag_part_nspin4(nrxx, chr); diff --git a/source/source_hamilt/module_xc/xc_functional_gga_noncol_sf_builtin.cpp b/source/source_hamilt/module_xc/xc_functional_gga_noncol_sf_builtin.cpp index eeaa792c77..adec74973e 100644 --- a/source/source_hamilt/module_xc/xc_functional_gga_noncol_sf_builtin.cpp +++ b/source/source_hamilt/module_xc/xc_functional_gga_noncol_sf_builtin.cpp @@ -25,6 +25,26 @@ std::tuple v_xc_ncgga_sf_builtin( if (PARAM.inp.nspin != 4 || (!PARAM.globalv.domag && !PARAM.globalv.domag_z)) throw std::domain_error("v_xc_ncgga_sf_builtin requires NSPIN==4."); + // ====================================================================== + // Scalmani-Frisch (SF) builtin for noncollinear GGA (gga_grad=3) + // + // Reference: Scalmani & Frisch, JCTC 8, 1069 (2012) + // + // Key formulas: + // rho_up = 0.5*(rho + |m|), rho_dn = 0.5*(rho - |m|) + // m_hat = m / |m| (magnetization unit vector) + // + // Gradient decomposition (chain rule via m_hat): + // grad(rho_up) = 0.5 grad(rho) + 0.5 m_hat_mu * grad(m_mu) + // grad(rho_dn) = 0.5 grad(rho) - 0.5 m_hat_mu * grad(m_mu) + // + // For GGA, the sigma = |grad(rho_up)|^2 etc. are computed from gdr1,gdr2, + // then passed to spin-polarized XC functionals (xc_spin, gcx_spin, gcc_spin). + // The resulting potential is converted back to nspin=4 representation: + // v_tot = 0.5*(v_up + v_dn) + // v_mag = 0.5*(v_up - v_dn) * m_hat + // ====================================================================== + ModulePW::PW_Basis* rhopw = chr->rhopw; const int npw = rhopw->npw; const double e2 = ModuleBase::e2; @@ -33,6 +53,7 @@ std::tuple v_xc_ncgga_sf_builtin( const double fac = 0.5; const bool is_gga = (XC_Functional::get_func_type() == 2 || XC_Functional::get_func_type() == 4); + // Step 1: decompose charge density into spin-up/spin-down and compute m_hat std::vector rhotmp1(nrxx), rhotmp2(nrxx), amag(nrxx); std::vector mag_part(3 * nrxx, 0.0); @@ -83,6 +104,9 @@ std::tuple v_xc_ncgga_sf_builtin( } } + // Step 3: LDA contribution (xc_spin) and convert to nspin=4 potential + // v(0,:) = 0.5*(vxc[0] + vxc[1]) -- total density channel + // v(1..3,:) = 0.5*(vxc[0] - vxc[1]) * m_mu / |m| -- magnetization channels double etxc = 0, vtxc = 0; ModuleBase::matrix v(PARAM.inp.nspin, nrxx); @@ -114,6 +138,9 @@ std::tuple v_xc_ncgga_sf_builtin( etxc += e2 * exc * arho; } + // Step 4: GGA contribution (gcx_spin + gcc_spin) and divergence correction (grad_dot) + // The gradient correction uses the same SF-decomposed gdr1, gdr2 and converts + // back to nspin=4 via m_hat, matching the LDA conversion above. if (is_gga) { double etxcgc = 0, vtxcgc = 0; diff --git a/source/source_hamilt/module_xc/xc_pot.cpp b/source/source_hamilt/module_xc/xc_pot.cpp index 6627b22698..89568f71f8 100644 --- a/source/source_hamilt/module_xc/xc_pot.cpp +++ b/source/source_hamilt/module_xc/xc_pot.cpp @@ -52,6 +52,10 @@ std::tuple XC_Functional::v_xc( #endif } + // For non-libxc builds: gga_grad=3 uses the SF builtin that computes the + // full nspin=4 GGA potential via the Scalmani-Frisch transformation. + // This path handles spin-up/spin-down decomposition and gradient + // corrections internally, returning (etxc, vtxc, v) directly. if (PARAM.inp.nspin == 4 && (PARAM.globalv.domag || PARAM.globalv.domag_z) && PARAM.inp.gga_grad == 3) { return ModuleXC::NCGGA_SF_Builtin::v_xc_ncgga_sf_builtin(nrxx, ucell->omega, ucell->tpiba, chr);