comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-28 14:15:19 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Support ROKS active-space FCI solver (#5203)

Browse source 
Co-authored-by: Thomas D. Kuehne <tkuehne@cp2k.org>
This commit is contained in:
Dynamics of Condensed Matter 2026-05-15 14:15:44 +02:00 • committed by GitHub
parent 0f818b0565
commit 2256cec575
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
7 changed files with 102 additions and 15 deletions
Download patch file Download diff file Expand all files Collapse all files

12
src/qs_active_space_fci.F
View file

@ -73,13 +73,14 @@ CONTAINS
CHARACTER(LEN=512) :: error_text
INTEGER :: max_iter, ms2, n2, n4, nmo_active, &
nspins, status
LOGICAL :: ionode
LOGICAL :: ionode, restricted_orbitals
REAL(KIND=dp) :: energy_active, threshold
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eri_aa, eri_ab, eri_bb, fock_a, fock_b, &
p_beta
nmo_active = active_space_env%nmo_active
nspins = active_space_env%nspins
restricted_orbitals = active_space_env%restricted_orbitals
n2 = nmo_active*nmo_active
n4 = n2*n2
ionode = para_env%is_source()
@ -104,8 +105,13 @@ CONTAINS
ASSOCIATE (act_indices => active_space_env%active_orbitals(:, 2))
CALL subspace_matrix_to_array(active_space_env%fock_sub(2), fock_b(1:n2), act_indices, act_indices)
END ASSOCIATE
CALL eri_to_array(active_space_env%eri, eri_ab(1:n4), active_space_env%active_orbitals, 1, 2)
CALL eri_to_array(active_space_env%eri, eri_bb(1:n4), active_space_env%active_orbitals, 2, 2)
IF (restricted_orbitals) THEN
eri_ab(1:n4) = eri_aa
eri_bb(1:n4) = eri_aa
ELSE
CALL eri_to_array(active_space_env%eri, eri_ab(1:n4), active_space_env%active_orbitals, 1, 2)
CALL eri_to_array(active_space_env%eri, eri_bb(1:n4), active_space_env%active_orbitals, 2, 2)
END IF
END IF
status = 0

31
src/qs_active_space_methods.F
View file

@ -378,6 +378,7 @@ CONTAINS
active_space_env%nelec_total = nelec_total
active_space_env%nspins = nspins
active_space_env%multiplicity = dft_control%multiplicity
active_space_env%restricted_orbitals = dft_control%roks
! define the active/inactive space orbitals
CALL section_vals_val_get(as_input, "ACTIVE_ORBITALS", explicit=explicit, i_val=nmo_active)
@ -505,7 +506,9 @@ CONTAINS
! create canonical orbitals
CALL get_qs_env(qs_env, scf_control=scf_control)
IF (dft_control%roks .AND. scf_control%roks_scheme /= high_spin_roks) THEN
CPABORT("Unclear how we define MOs in the general restricted case ... stopping")
CALL cp_abort(__LOCATION__, &
"Only high-spin ROKS is supported for ACTIVE_SPACE FCI; "// &
"general ROKS MO definitions are not implemented.")
ELSE
IF (dft_control%do_admm) THEN
IF (dft_control%do_admm_mo) THEN
@ -641,7 +644,9 @@ CONTAINS
CASE (manual_selection)
! create canonical orbitals
IF (dft_control%roks) THEN
CPABORT("Unclear how we define MOs in the restricted case ... stopping")
CALL cp_abort(__LOCATION__, &
"Manual ACTIVE_SPACE orbital selection is not supported for ROKS; "// &
"use canonical high-spin ROKS.")
ELSE
IF (dft_control%do_admm) THEN
! For admm_mo, the auxiliary density is computed from the MOs, which never change
@ -2293,7 +2298,7 @@ CONTAINS
LOGICAL, INTENT(IN) :: restricted
INTEGER :: i, i1, i2, i3, i4, isym, iw, m1, m2, &
nmo, norb, nspins
ms2, nmo, norb, nspins
REAL(KIND=dp) :: checksum, esub
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: fmat
TYPE(cp_logger_type), POINTER :: logger
@ -2307,10 +2312,10 @@ CONTAINS
!
nspins = active_space_env%nspins
norb = SIZE(active_space_env%active_orbitals, 1)
ms2 = active_space_env%multiplicity - 1
IF (nspins == 1 .OR. restricted) THEN
! Closed shell or restricted open-shell
ASSOCIATE (ms2 => active_space_env%multiplicity, &
nelec => active_space_env%nelec_active)
ASSOCIATE (nelec => active_space_env%nelec_active)
IF (iw > 0) THEN
WRITE (iw, "(A,A,I4,A,I4,A,I2,A)") "&FCI", " NORB=", norb, ",NELEC=", nelec, ",MS2=", ms2, ","
@ -2353,8 +2358,7 @@ CONTAINS
END ASSOCIATE
ELSE
ASSOCIATE (ms2 => active_space_env%multiplicity, &
nelec => active_space_env%nelec_active)
ASSOCIATE (nelec => active_space_env%nelec_active)
IF (iw > 0) THEN
WRITE (iw, "(A,A,I4,A,I4,A,I2,A)") "&FCI", " NORB=", norb, ",NELEC=", nelec, ",MS2=", ms2, ","
@ -3055,7 +3059,6 @@ CONTAINS
iw = cp_logger_get_default_io_unit(logger)
CALL get_qs_env(qs_env, para_env=para_env, dft_control=dft_control)
IF (dft_control%roks) CPABORT("AS_SOLVER FCI is not supported for ROKS calculations.")
CALL section_vals_val_get(as_input, "SCF_EMBEDDING", l_val=do_scf_embedding)
active_space_env%do_scf_embedding = do_scf_embedding
@ -3466,7 +3469,7 @@ CONTAINS
CHARACTER(len=default_string_length) :: header
INTEGER :: handle, iw
LOGICAL :: ionode
LOGICAL :: ionode, restricted_orbitals
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eri_aa, eri_ab, eri_bb, s_ab
TYPE(cp_logger_type), POINTER :: logger
@ -3475,14 +3478,20 @@ CONTAINS
logger => cp_get_default_logger()
iw = cp_logger_get_default_io_unit(logger)
ionode = para_env%is_source()
restricted_orbitals = active_space_env%restricted_orbitals
ALLOCATE (eri_aa(active_space_env%nmo_active**4))
CALL eri_to_array(active_space_env%eri, eri_aa, active_space_env%active_orbitals, 1, 1)
IF (active_space_env%nspins == 2) THEN
ALLOCATE (eri_ab(active_space_env%nmo_active**4))
CALL eri_to_array(active_space_env%eri, eri_ab, active_space_env%active_orbitals, 1, 2)
ALLOCATE (eri_bb(active_space_env%nmo_active**4))
CALL eri_to_array(active_space_env%eri, eri_bb, active_space_env%active_orbitals, 2, 2)
IF (restricted_orbitals) THEN
eri_ab(:) = eri_aa
eri_bb(:) = eri_aa
ELSE
CALL eri_to_array(active_space_env%eri, eri_ab, active_space_env%active_orbitals, 1, 2)
CALL eri_to_array(active_space_env%eri, eri_bb, active_space_env%active_orbitals, 2, 2)
END IF
! get the overlap_ab matrix into Fortran array
ALLOCATE (s_ab(active_space_env%nmo_active**2))
ASSOCIATE (act_indices_a => active_space_env%active_orbitals(:, 1), &

1
src/qs_active_space_types.F
View file

@ -94,6 +94,7 @@ MODULE qs_active_space_types
INTEGER :: nmo_inactive = 0
INTEGER :: multiplicity = 0
INTEGER :: nspins = 0
LOGICAL :: restricted_orbitals = .FALSE.
LOGICAL :: molecule = .FALSE.
INTEGER :: model = 0
REAL(KIND=dp) :: energy_total = 0.0_dp

3
tests/QS/regtest-as-2/TEST_FILES.toml
View file

@ -23,5 +23,6 @@
{matcher="M092", tol=1e-8, ref=4.08599253}]
"h2_gpw_ht_roks.inp" = [{matcher="E_total", tol=1e-12, ref=-0.54761910555768},
{matcher="M092", tol=1e-8, ref=4.386293}]
{matcher="M092", tol=1e-8, ref=4.386293},
{matcher="FCIDUMP_MS2", tol=0.0, ref=1}]
#EOF

1
tests/QS/regtest-as-fci/TEST_FILES.toml
View file

@ -1,3 +1,4 @@
"h2_as_fci.inp" = [{matcher="M011", tol=1e-12, ref=-1.157903329414417}]
"h2_as_fci_uks.inp" = [{matcher="M011", tol=1e-10, ref=-1.157902959723317}]
"h2_as_fci_roks.inp" = [{matcher="M011", tol=1e-12, ref=-0.547619105557678}]
#EOF

66
tests/QS/regtest-as-fci/h2_as_fci_roks.inp Normal file
View file

@ -0,0 +1,66 @@
&GLOBAL
PRINT_LEVEL LOW
PROJECT h2_as_fci_roks
RUN_TYPE ENERGY
&END GLOBAL
&FORCE_EVAL
METHOD QS
&DFT
CHARGE 1
MULTIPLICITY 2
ROKS T
&ACTIVE_SPACE
ACTIVE_ELECTRONS 1
ACTIVE_ORBITALS 2
AS_SOLVER FCI
MAX_ITER 1
SCF_EMBEDDING F
&ERI
METHOD GPW_HALF_TRANSFORM
PERIODICITY 0 0 0
&END ERI
&ERI_GPW
CUTOFF 500
&END ERI_GPW
&END ACTIVE_SPACE
&MGRID
CUTOFF 500
&END MGRID
&POISSON
PERIODIC NONE
POISSON_SOLVER ANALYTIC
&END POISSON
&QS
METHOD GPW
&END QS
&SCF
ADDED_MOS 3
MAX_SCF 10
&PRINT
&RESTART OFF
&END RESTART
&END PRINT
&END SCF
&XC
&HF 1.0
&END HF
&XC_FUNCTIONAL NONE
&END XC_FUNCTIONAL
&END XC
&END DFT
&SUBSYS
&CELL
ABC 6.0 6.0 6.0
PERIODIC NONE
&END CELL
&COORD
H 0.000 0.000 0.356
H 0.000 0.000 -0.356
&END COORD
&KIND H
BASIS_SET 6-31G*
POTENTIAL GTH-HF
&END KIND
&END SUBSYS
&END FORCE_EVAL

3
tests/matchers.py
View file

@ -223,6 +223,9 @@ registry["M089"] = GenericMatcher(r"Electronic density on regular grids:", col=7
registry["M090"] = GenericMatcher(r"Final localization:", col=3)
registry["M091"] = GenericMatcher(r"Ionization potentials for XPS", col=8)
registry["M092"] = GenericMatcher(r"FCIDUMP| Checksum:", col=3)
registry["FCIDUMP_MS2"] = GenericMatcher(
r"&FCI .*MS2=\s*([-+0-9]+),", col=1, regex=True
)
registry["M093"] = GenericMatcher(r"SPGR| SPACE GROUP NUMBER:", col=5)
registry["M094"] = GenericMatcher(r"KS CSR write|", col=4)
registry["M095"] = GenericMatcher(r"Fermi energy:", col=3)