Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Utility functions for CNEO-DFT
10 : !> (see J. Chem. Theory Comput. 2025, 21, 16, 7865–7877)
11 : !> \par History
12 : !> 08.2025 created [zc62]
13 : !> \author Zehua Chen
14 : ! **************************************************************************************************
15 : MODULE qs_cneo_utils
16 : USE ao_util, ONLY: trace_r_AxB
17 : USE basis_set_types, ONLY: get_gto_basis_set,&
18 : gto_basis_set_type
19 : USE kinds, ONLY: dp
20 : USE memory_utilities, ONLY: reallocate
21 : USE orbital_pointers, ONLY: indso,&
22 : nsoset
23 : USE qs_harmonics_atom, ONLY: get_none0_cg_list,&
24 : harmonics_atom_type
25 : USE spherical_harmonics, ONLY: clebsch_gordon,&
26 : clebsch_gordon_deallocate,&
27 : clebsch_gordon_init
28 : #include "./base/base_uses.f90"
29 :
30 : IMPLICIT NONE
31 :
32 : PRIVATE
33 :
34 : ! *** Global parameters ***
35 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_cneo_utils'
36 :
37 : ! *** Public subroutines ***
38 : PUBLIC :: atom_solve_cneo, cneo_gather, cneo_scatter, create_harmonics_atom_cneo, &
39 : create_my_CG_cneo, get_maxl_CG_cneo
40 :
41 : CONTAINS
42 :
43 : ! **************************************************************************************************
44 : !> \brief Mostly copied from qs_rho_atom_methods::init_rho_atom
45 : !> \param my_CG ...
46 : !> \param lcleb ...
47 : !> \param maxl ...
48 : !> \param llmax ...
49 : ! **************************************************************************************************
50 8 : SUBROUTINE create_my_CG_cneo(my_CG, lcleb, maxl, llmax)
51 :
52 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: my_CG
53 : INTEGER, INTENT(IN) :: lcleb, maxl, llmax
54 :
55 : INTEGER :: il, iso, iso1, iso2, l1, l1l2, l2, lc1, &
56 : lc2, lp, m1, m2, mm, mp
57 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: rga
58 :
59 : ! *** allocate calculate the CG coefficients up to the maxl ***
60 8 : CALL clebsch_gordon_init(lcleb)
61 :
62 24 : ALLOCATE (rga(lcleb, 2))
63 32 : DO lc1 = 0, maxl
64 104 : DO iso1 = nsoset(lc1 - 1) + 1, nsoset(lc1)
65 72 : l1 = indso(1, iso1)
66 72 : m1 = indso(2, iso1)
67 312 : DO lc2 = 0, maxl
68 936 : DO iso2 = nsoset(lc2 - 1) + 1, nsoset(lc2)
69 648 : l2 = indso(1, iso2)
70 648 : m2 = indso(2, iso2)
71 648 : CALL clebsch_gordon(l1, m1, l2, m2, rga)
72 648 : IF (l1 + l2 > llmax) THEN
73 : l1l2 = llmax
74 : ELSE
75 : l1l2 = l1 + l2
76 : END IF
77 648 : mp = m1 + m2
78 648 : mm = m1 - m2
79 648 : IF (m1*m2 < 0 .OR. (m1*m2 == 0 .AND. (m1 < 0 .OR. m2 < 0))) THEN
80 288 : mp = -ABS(mp)
81 288 : mm = -ABS(mm)
82 : ELSE
83 360 : mp = ABS(mp)
84 360 : mm = ABS(mm)
85 : END IF
86 2304 : DO lp = MOD(l1 + l2, 2), l1l2, 2
87 1440 : il = lp/2 + 1
88 1440 : IF (ABS(mp) <= lp) THEN
89 1024 : IF (mp >= 0) THEN
90 688 : iso = nsoset(lp - 1) + lp + 1 + mp
91 : ELSE
92 336 : iso = nsoset(lp - 1) + lp + 1 - ABS(mp)
93 : END IF
94 1024 : my_CG(iso1, iso2, iso) = rga(il, 1)
95 : END IF
96 2088 : IF (mp /= mm .AND. ABS(mm) <= lp) THEN
97 480 : IF (mm >= 0) THEN
98 320 : iso = nsoset(lp - 1) + lp + 1 + mm
99 : ELSE
100 160 : iso = nsoset(lp - 1) + lp + 1 - ABS(mm)
101 : END IF
102 480 : my_CG(iso1, iso2, iso) = rga(il, 2)
103 : END IF
104 : END DO
105 : END DO ! iso2
106 : END DO ! lc2
107 : END DO ! iso1
108 : END DO ! lc1
109 8 : DEALLOCATE (rga)
110 8 : CALL clebsch_gordon_deallocate()
111 :
112 8 : END SUBROUTINE create_my_CG_cneo
113 :
114 : ! **************************************************************************************************
115 : !> \brief Mostly copied from qs_harmonics_atom::create_harmonics_atom
116 : !> \param harmonics ...
117 : !> \param my_CG ...
118 : !> \param llmax ...
119 : !> \param maxs ...
120 : !> \param max_s_harm ...
121 : ! **************************************************************************************************
122 8 : SUBROUTINE create_harmonics_atom_cneo(harmonics, my_CG, llmax, maxs, max_s_harm)
123 :
124 : TYPE(harmonics_atom_type), POINTER :: harmonics
125 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: my_CG
126 : INTEGER, INTENT(IN) :: llmax, maxs, max_s_harm
127 :
128 : INTEGER :: i, is
129 :
130 8 : CPASSERT(ASSOCIATED(harmonics))
131 :
132 8 : harmonics%max_s_harm = max_s_harm
133 8 : harmonics%llmax = llmax
134 :
135 8 : NULLIFY (harmonics%my_CG, harmonics%my_CG_dxyz, harmonics%my_CG_dxyz_asym)
136 8 : CALL reallocate(harmonics%my_CG, 1, maxs, 1, maxs, 1, max_s_harm)
137 :
138 208 : DO i = 1, max_s_harm
139 2008 : DO is = 1, maxs
140 36200 : harmonics%my_CG(1:maxs, is, i) = my_CG(1:maxs, is, i)
141 : END DO
142 : END DO
143 :
144 8 : END SUBROUTINE create_harmonics_atom_cneo
145 :
146 : ! **************************************************************************************************
147 : !> \brief Mostly copied from qs_harmonics_atom::get_maxl_CG
148 : !> \param harmonics ...
149 : !> \param orb_basis ...
150 : !> \param llmax ...
151 : !> \param max_s_harm ...
152 : ! **************************************************************************************************
153 16 : SUBROUTINE get_maxl_CG_cneo(harmonics, orb_basis, llmax, max_s_harm)
154 :
155 : TYPE(harmonics_atom_type), POINTER :: harmonics
156 : TYPE(gto_basis_set_type), POINTER :: orb_basis
157 : INTEGER, INTENT(IN) :: llmax, max_s_harm
158 :
159 : INTEGER :: is1, is2, itmp, max_iso_not0, nset
160 8 : INTEGER, DIMENSION(:), POINTER :: lmax, lmin
161 :
162 0 : CPASSERT(ASSOCIATED(harmonics))
163 :
164 8 : CALL get_gto_basis_set(gto_basis_set=orb_basis, lmax=lmax, lmin=lmin, nset=nset)
165 :
166 : ! *** Assign indices for the non null CG coefficients ***
167 8 : max_iso_not0 = 0
168 80 : DO is1 = 1, nset
169 728 : DO is2 = 1, nset
170 : CALL get_none0_cg_list(harmonics%my_CG, &
171 : lmin(is1), lmax(is1), lmin(is2), lmax(is2), &
172 648 : max_s_harm, llmax, max_iso_not0=itmp)
173 720 : max_iso_not0 = MAX(max_iso_not0, itmp)
174 : END DO ! is2
175 : END DO ! is1
176 8 : harmonics%max_iso_not0 = max_iso_not0
177 :
178 8 : END SUBROUTINE get_maxl_CG_cneo
179 :
180 : ! **************************************************************************************************
181 : !> \brief Mostly copied from atom_utils::atom_solve
182 : !> \param hmat ...
183 : !> \param f ...
184 : !> \param umat ...
185 : !> \param orb ...
186 : !> \param ener ...
187 : !> \param pmat ...
188 : !> \param r ...
189 : !> \param dist ...
190 : !> \param nb ...
191 : !> \param nv ...
192 : ! **************************************************************************************************
193 645 : SUBROUTINE atom_solve_cneo(hmat, f, umat, orb, ener, pmat, r, dist, nb, nv)
194 :
195 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: hmat
196 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: f
197 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: umat
198 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: orb
199 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: ener
200 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: pmat
201 : REAL(KIND=dp), DIMENSION(3), INTENT(INOUT) :: r
202 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: dist
203 : INTEGER, INTENT(IN) :: nb, nv
204 :
205 : INTEGER :: info, lwork, m, n
206 645 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: w, work
207 645 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: a, b, h_fx
208 :
209 645 : CPASSERT(nb >= nv)
210 :
211 356685 : orb = 0._dp
212 645 : n = nb
213 645 : m = nv
214 645 : IF (n > 0 .AND. m > 0) THEN
215 645 : lwork = MAX(n*n, n + 100)
216 7095 : ALLOCATE (a(m, m), b(n, m), w(m), work(lwork))
217 2580 : IF (DOT_PRODUCT(f, f) /= 0.0_dp) THEN
218 2396 : ALLOCATE (h_fx(n, n))
219 : h_fx(1:n, 1:n) = hmat(1:n, 1:n) + f(1)*dist(1:n, 1:n, 1) + &
220 331247 : f(2)*dist(1:n, 1:n, 2) + f(3)*dist(1:n, 1:n, 3)
221 599 : CALL dgemm("N", "N", n, m, n, 1.0_dp, h_fx, n, umat, n, 0.0_dp, b, n)
222 599 : DEALLOCATE (h_fx)
223 : ELSE
224 46 : CALL dgemm("N", "N", n, m, n, 1.0_dp, hmat, n, umat, n, 0.0_dp, b, n)
225 : END IF
226 645 : CALL dgemm("T", "N", m, m, n, 1.0_dp, umat, n, b, n, 0.0_dp, a, m)
227 645 : CALL dsyev("V", "U", m, a, m, w, work, lwork, info)
228 645 : CALL dgemm("N", "N", n, m, m, 1.0_dp, umat, n, a, m, 0.0_dp, b, n)
229 :
230 645 : m = MIN(m, SIZE(orb, 2))
231 356685 : orb(1:n, 1:m) = b(1:n, 1:m)
232 15480 : ener(1:m) = w(1:m)
233 :
234 645 : DEALLOCATE (a, b, w, work)
235 :
236 : ! calculate the density matrix using the orbital with the lowest orbital energy
237 356685 : pmat = 0.0_dp
238 645 : CALL dger(n, n, 1.0_dp, orb(:, 1), 1, orb(:, 1), 1, pmat, n)
239 : ! calculate the expectation position (basis center as the origin)
240 : r = [trace_r_AxB(dist(1:n, 1:n, 1), n, pmat, n, n, n), &
241 : trace_r_AxB(dist(1:n, 1:n, 2), n, pmat, n, n, n), &
242 3225 : trace_r_AxB(dist(1:n, 1:n, 3), n, pmat, n, n, n)]
243 : END IF
244 :
245 645 : END SUBROUTINE atom_solve_cneo
246 :
247 : ! **************************************************************************************************
248 : !> \brief Mostly copied from qs_oce_methods::prj_gather
249 : !> \param ain ...
250 : !> \param aout ...
251 : !> \param nbas ...
252 : !> \param n2oindex ...
253 : ! **************************************************************************************************
254 92 : SUBROUTINE cneo_gather(ain, aout, nbas, n2oindex)
255 :
256 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: ain
257 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: aout
258 : INTEGER, INTENT(IN) :: nbas
259 : INTEGER, DIMENSION(:), POINTER :: n2oindex
260 :
261 : INTEGER :: i, ip, j, jp
262 :
263 2208 : DO i = 1, nbas
264 2116 : ip = n2oindex(i)
265 50876 : DO j = 1, nbas
266 48668 : jp = n2oindex(j)
267 50784 : aout(j, i) = ain(jp, ip)
268 : END DO
269 : END DO
270 :
271 92 : END SUBROUTINE cneo_gather
272 :
273 : ! **************************************************************************************************
274 : !> \brief Mostly copied from qs_oce_methods::prj_scatter
275 : !> \param ain ...
276 : !> \param aout ...
277 : !> \param nbas ...
278 : !> \param n2oindex ...
279 : ! **************************************************************************************************
280 156 : SUBROUTINE cneo_scatter(ain, aout, nbas, n2oindex)
281 :
282 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: ain
283 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: aout
284 : INTEGER, INTENT(IN) :: nbas
285 : INTEGER, DIMENSION(:), POINTER :: n2oindex
286 :
287 : INTEGER :: i, ip, j, jp
288 :
289 3744 : DO i = 1, nbas
290 3588 : ip = n2oindex(i)
291 86268 : DO j = 1, nbas
292 82524 : jp = n2oindex(j)
293 86112 : aout(jp, ip) = aout(jp, ip) + ain(j, i)
294 : END DO
295 : END DO
296 :
297 156 : END SUBROUTINE cneo_scatter
298 :
299 : END MODULE qs_cneo_utils
|
