Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Local FCI solver interface for active-space embedding.
10 : ! **************************************************************************************************
11 : MODULE qs_active_space_fci
12 : #if defined(__LIBFCI)
13 : USE ISO_C_BINDING, ONLY: C_CHAR, &
14 : C_DOUBLE, &
15 : C_INT, &
16 : C_NULL_CHAR
17 : #endif
18 : USE kinds, ONLY: dp
19 : USE message_passing, ONLY: mp_para_env_type
20 : USE qs_active_space_types, ONLY: active_space_type
21 : #if defined(__LIBFCI)
22 : USE qs_active_space_utils, ONLY: eri_to_array, &
23 : subspace_matrix_to_array
24 : #endif
25 : #include "./base/base_uses.f90"
26 :
27 : IMPLICIT NONE
28 : PRIVATE
29 :
30 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_active_space_fci'
31 :
32 : PUBLIC :: solve_active_space_fci
33 :
34 : #if defined(__LIBFCI)
35 : INTERFACE
36 : FUNCTION libfci_solve(nspins, norb, nelec, ms2, h1_a, h1_b, eri_aa, eri_ab, eri_bb, &
37 : max_iter, threshold, pspace, project_spin, verbose, energy, &
38 : rdm1_a, rdm1_b, err_msg, err_msg_len) RESULT(status) &
39 : BIND(C, name="libfci_solve")
40 : IMPORT :: C_CHAR, C_DOUBLE, C_INT
41 : INTEGER(C_INT), VALUE :: nspins, norb, nelec, ms2
42 : REAL(C_DOUBLE), DIMENSION(*), INTENT(IN) :: h1_a, h1_b, eri_aa, eri_ab, eri_bb
43 : INTEGER(C_INT), VALUE :: max_iter, pspace, project_spin, verbose
44 : REAL(C_DOUBLE), VALUE :: threshold
45 : REAL(C_DOUBLE), INTENT(OUT) :: energy
46 : REAL(C_DOUBLE), DIMENSION(*), INTENT(OUT) :: rdm1_a, rdm1_b
47 : CHARACTER(KIND=C_CHAR), DIMENSION(*), INTENT(OUT) :: err_msg
48 : INTEGER(C_INT), VALUE :: err_msg_len
49 : INTEGER(C_INT) :: status
50 : END FUNCTION libfci_solve
51 : END INTERFACE
52 : #endif
53 :
54 : CONTAINS
55 :
56 : ! **************************************************************************************************
57 : !> \brief Solve the current active-space Hamiltonian with the local FCI kernel.
58 : !> \param active_space_env active-space environment
59 : !> \param para_env parallel environment
60 : !> \param p_act_mo_a alpha or spin-summed active-space 1-RDM
61 : !> \param p_act_mo_b beta active-space 1-RDM for unrestricted calculations
62 : ! **************************************************************************************************
63 8 : SUBROUTINE solve_active_space_fci(active_space_env, para_env, p_act_mo_a, p_act_mo_b)
64 : TYPE(active_space_type), INTENT(INOUT), POINTER :: active_space_env
65 : TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env
66 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
67 : INTENT(OUT) :: p_act_mo_a
68 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
69 : INTENT(OUT), OPTIONAL :: p_act_mo_b
70 :
71 : #if defined(__LIBFCI)
72 : CHARACTER(KIND=C_CHAR), DIMENSION(512) :: c_error
73 : CHARACTER(LEN=512) :: error_text
74 : INTEGER :: max_iter, ms2, n2, n4, nmo_active, &
75 : nspins, status
76 : LOGICAL :: ionode, restricted_orbitals
77 : REAL(KIND=dp) :: energy_active, threshold
78 8 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eri_aa, eri_ab, eri_bb, fock_a, fock_b, &
79 8 : p_beta
80 :
81 8 : nmo_active = active_space_env%nmo_active
82 8 : nspins = active_space_env%nspins
83 8 : restricted_orbitals = active_space_env%restricted_orbitals
84 8 : n2 = nmo_active*nmo_active
85 8 : n4 = n2*n2
86 8 : ionode = para_env%is_source()
87 :
88 48 : ALLOCATE (fock_a(n2), eri_aa(n4), p_act_mo_a(n2))
89 56 : ALLOCATE (fock_b(MAX(1, n2)), eri_ab(MAX(1, n4)), eri_bb(MAX(1, n4)), p_beta(MAX(1, n2)))
90 8 : fock_b(:) = 0.0_dp
91 8 : eri_ab(:) = 0.0_dp
92 8 : eri_bb(:) = 0.0_dp
93 8 : p_beta(:) = 0.0_dp
94 :
95 : ASSOCIATE (act_indices => active_space_env%active_orbitals(:, 1))
96 8 : CALL subspace_matrix_to_array(active_space_env%fock_sub(1), fock_a, act_indices, act_indices)
97 : END ASSOCIATE
98 8 : CALL eri_to_array(active_space_env%eri, eri_aa, active_space_env%active_orbitals, 1, 1)
99 :
100 8 : IF (nspins == 2) THEN
101 6 : IF (.NOT. PRESENT(p_act_mo_b)) THEN
102 0 : CPABORT("Missing beta output buffer for local active-space FCI solver.")
103 : END IF
104 12 : ALLOCATE (p_act_mo_b(n2))
105 : ASSOCIATE (act_indices => active_space_env%active_orbitals(:, 2))
106 6 : CALL subspace_matrix_to_array(active_space_env%fock_sub(2), fock_b(1:n2), act_indices, act_indices)
107 : END ASSOCIATE
108 6 : IF (restricted_orbitals) THEN
109 34 : eri_ab(1:n4) = eri_aa
110 34 : eri_bb(1:n4) = eri_aa
111 : ELSE
112 4 : CALL eri_to_array(active_space_env%eri, eri_ab(1:n4), active_space_env%active_orbitals, 1, 2)
113 4 : CALL eri_to_array(active_space_env%eri, eri_bb(1:n4), active_space_env%active_orbitals, 2, 2)
114 : END IF
115 : END IF
116 :
117 8 : status = 0
118 8 : energy_active = 0.0_dp
119 8 : error_text = ""
120 4104 : c_error(:) = C_NULL_CHAR
121 8 : ms2 = active_space_env%multiplicity - 1
122 8 : max_iter = 2048
123 8 : threshold = 1.0E-8_dp
124 :
125 8 : IF (ionode) THEN
126 : status = libfci_solve(INT(nspins, C_INT), INT(nmo_active, C_INT), &
127 : INT(active_space_env%nelec_active, C_INT), INT(ms2, C_INT), &
128 : fock_a, fock_b, eri_aa, eri_ab, eri_bb, INT(max_iter, C_INT), &
129 : REAL(threshold, C_DOUBLE), INT(200, C_INT), INT(1, C_INT), INT(0, C_INT), &
130 4 : energy_active, p_act_mo_a, p_beta, c_error, INT(SIZE(c_error), C_INT))
131 4 : error_text = c_string(c_error)
132 : END IF
133 :
134 8 : CALL para_env%bcast(status, para_env%source)
135 8 : CALL para_env%bcast(error_text, para_env%source)
136 8 : IF (status /= 0) THEN
137 0 : IF (LEN_TRIM(error_text) > 0) THEN
138 0 : CPABORT("Local active-space FCI solver failed: "//TRIM(error_text))
139 : ELSE
140 0 : CPABORT("Local active-space FCI solver failed.")
141 : END IF
142 : END IF
143 :
144 8 : CALL para_env%bcast(energy_active, para_env%source)
145 8 : CALL para_env%bcast(p_act_mo_a, para_env%source)
146 8 : active_space_env%energy_active = energy_active
147 8 : IF (nspins == 2) THEN
148 30 : p_act_mo_b(:) = p_beta(1:n2)
149 6 : CALL para_env%bcast(p_act_mo_b, para_env%source)
150 : END IF
151 :
152 8 : DEALLOCATE (fock_a, fock_b, eri_aa, eri_ab, eri_bb, p_beta)
153 : #else
154 : MARK_USED(active_space_env)
155 : MARK_USED(para_env)
156 : MARK_USED(p_act_mo_a)
157 : IF (PRESENT(p_act_mo_b)) THEN
158 : MARK_USED(p_act_mo_b)
159 : END IF
160 : CPABORT("AS_SOLVER FCI requires CP2K to be built with LibFCI support.")
161 : #endif
162 :
163 8 : END SUBROUTINE solve_active_space_fci
164 :
165 : #if defined(__LIBFCI)
166 : ! **************************************************************************************************
167 : !> \brief Convert a C null-terminated string buffer to a Fortran string.
168 : !> \param buffer ...
169 : !> \return ...
170 : ! **************************************************************************************************
171 4 : FUNCTION c_string(buffer) RESULT(text)
172 : CHARACTER(KIND=C_CHAR), DIMENSION(:), INTENT(IN) :: buffer
173 : CHARACTER(LEN=512) :: text
174 :
175 : INTEGER :: i
176 :
177 4 : text = ""
178 4 : DO i = 1, MIN(SIZE(buffer), LEN(text))
179 4 : IF (buffer(i) == C_NULL_CHAR) EXIT
180 4 : text(i:i) = buffer(i)
181 : END DO
182 4 : END FUNCTION c_string
183 : #endif
184 :
185 : END MODULE qs_active_space_fci
|
