From 2256cec57518a66cedf6d7165fb4dd3f977fff17 Mon Sep 17 00:00:00 2001 From: Dynamics of Condensed Matter <30792324+DCM-Uni-Paderborn@users.noreply.github.com> Date: Fri, 15 May 2026 14:15:44 +0200 Subject: [PATCH] Support ROKS active-space FCI solver (#5203) Co-authored-by: Thomas D. Kuehne --- src/qs_active_space_fci.F | 12 +++- src/qs_active_space_methods.F | 31 ++++++---- src/qs_active_space_types.F | 1 + tests/QS/regtest-as-2/TEST_FILES.toml | 3 +- tests/QS/regtest-as-fci/TEST_FILES.toml | 1 + tests/QS/regtest-as-fci/h2_as_fci_roks.inp | 66 ++++++++++++++++++++++ tests/matchers.py | 3 + 7 files changed, 102 insertions(+), 15 deletions(-) create mode 100644 tests/QS/regtest-as-fci/h2_as_fci_roks.inp diff --git a/src/qs_active_space_fci.F b/src/qs_active_space_fci.F index 998ad6014c..cea8dab2fb 100644 --- a/src/qs_active_space_fci.F +++ b/src/qs_active_space_fci.F @@ -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 diff --git a/src/qs_active_space_methods.F b/src/qs_active_space_methods.F index 2b26aba609..a4bcfd1e69 100644 --- a/src/qs_active_space_methods.F +++ b/src/qs_active_space_methods.F @@ -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), & diff --git a/src/qs_active_space_types.F b/src/qs_active_space_types.F index c227c46905..8042e71158 100644 --- a/src/qs_active_space_types.F +++ b/src/qs_active_space_types.F @@ -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 diff --git a/tests/QS/regtest-as-2/TEST_FILES.toml b/tests/QS/regtest-as-2/TEST_FILES.toml index 7259b9e20f..bf35104f27 100644 --- a/tests/QS/regtest-as-2/TEST_FILES.toml +++ b/tests/QS/regtest-as-2/TEST_FILES.toml @@ -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 diff --git a/tests/QS/regtest-as-fci/TEST_FILES.toml b/tests/QS/regtest-as-fci/TEST_FILES.toml index 7524b40076..380c5b114c 100644 --- a/tests/QS/regtest-as-fci/TEST_FILES.toml +++ b/tests/QS/regtest-as-fci/TEST_FILES.toml @@ -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 diff --git a/tests/QS/regtest-as-fci/h2_as_fci_roks.inp b/tests/QS/regtest-as-fci/h2_as_fci_roks.inp new file mode 100644 index 0000000000..3f33d66a42 --- /dev/null +++ b/tests/QS/regtest-as-fci/h2_as_fci_roks.inp @@ -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 diff --git a/tests/matchers.py b/tests/matchers.py index f9f5e27e17..9d6785cd69 100644 --- a/tests/matchers.py +++ b/tests/matchers.py @@ -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)