Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Common routines for PAO parametrizations.
10 : !> \author Ole Schuett
11 : ! **************************************************************************************************
12 : MODULE pao_param_methods
13 : USE cp_control_types, ONLY: dft_control_type
14 : USE cp_dbcsr_api, ONLY: &
15 : dbcsr_add, dbcsr_complete_redistribute, dbcsr_create, dbcsr_get_block_p, dbcsr_get_info, &
16 : dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, dbcsr_iterator_start, &
17 : dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_multiply, dbcsr_p_type, dbcsr_release, &
18 : dbcsr_reserve_diag_blocks, dbcsr_scale, dbcsr_type
19 : USE cp_log_handling, ONLY: cp_to_string
20 : USE dm_ls_scf_qs, ONLY: matrix_decluster
21 : USE dm_ls_scf_types, ONLY: ls_mstruct_type,&
22 : ls_scf_env_type
23 : USE kinds, ONLY: dp
24 : USE message_passing, ONLY: mp_comm_type
25 : USE pao_types, ONLY: pao_env_type
26 : USE qs_environment_types, ONLY: get_qs_env,&
27 : qs_environment_type
28 : USE qs_rho_types, ONLY: qs_rho_get,&
29 : qs_rho_type
30 : #include "./base/base_uses.f90"
31 :
32 : IMPLICIT NONE
33 :
34 : PRIVATE
35 :
36 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pao_param_methods'
37 :
38 : PUBLIC :: pao_calc_grad_lnv_wrt_U, pao_calc_AB_from_U, pao_calc_grad_lnv_wrt_AB
39 :
40 : CONTAINS
41 :
42 : ! **************************************************************************************************
43 : !> \brief Helper routine, calculates partial derivative dE/dU
44 : !> \param qs_env ...
45 : !> \param ls_scf_env ...
46 : !> \param matrix_M_diag the derivate wrt U, matrix uses pao%diag_distribution
47 : ! **************************************************************************************************
48 2430 : SUBROUTINE pao_calc_grad_lnv_wrt_U(qs_env, ls_scf_env, matrix_M_diag)
49 : TYPE(qs_environment_type), POINTER :: qs_env
50 : TYPE(ls_scf_env_type), TARGET :: ls_scf_env
51 : TYPE(dbcsr_type) :: matrix_M_diag
52 :
53 : CHARACTER(len=*), PARAMETER :: routineN = 'pao_calc_grad_lnv_wrt_U'
54 :
55 : INTEGER :: handle
56 : REAL(KIND=dp) :: filter_eps
57 2430 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
58 : TYPE(dbcsr_type) :: matrix_M, matrix_Ma, matrix_Mb, matrix_NM
59 : TYPE(ls_mstruct_type), POINTER :: ls_mstruct
60 : TYPE(pao_env_type), POINTER :: pao
61 :
62 2430 : CALL timeset(routineN, handle)
63 :
64 2430 : ls_mstruct => ls_scf_env%ls_mstruct
65 2430 : pao => ls_scf_env%pao_env
66 2430 : filter_eps = ls_scf_env%eps_filter
67 2430 : CALL get_qs_env(qs_env, matrix_s=matrix_s)
68 :
69 2430 : CALL pao_calc_grad_lnv_wrt_AB(qs_env, ls_scf_env, matrix_Ma, matrix_Mb)
70 :
71 : ! Calculation uses distr. of matrix_s, afterwards we redistribute to pao%diag_distribution.
72 2430 : CALL dbcsr_create(matrix_M, template=matrix_s(1)%matrix, matrix_type="N")
73 2430 : CALL dbcsr_reserve_diag_blocks(matrix_M)
74 :
75 2430 : CALL dbcsr_create(matrix_NM, template=ls_mstruct%matrix_A, matrix_type="N")
76 :
77 : CALL dbcsr_multiply("N", "N", 1.0_dp, pao%matrix_N_inv, matrix_Ma, &
78 2430 : 1.0_dp, matrix_NM, filter_eps=filter_eps)
79 :
80 : CALL dbcsr_multiply("N", "N", 1.0_dp, pao%matrix_N, matrix_Mb, &
81 2430 : 1.0_dp, matrix_NM, filter_eps=filter_eps)
82 :
83 : CALL dbcsr_multiply("N", "T", 1.0_dp, matrix_NM, pao%matrix_Y, &
84 2430 : 1.0_dp, matrix_M, filter_eps=filter_eps)
85 :
86 : !---------------------------------------------------------------------------
87 : ! redistribute using pao%diag_distribution
88 : CALL dbcsr_create(matrix_M_diag, &
89 : name="PAO matrix_M", &
90 : matrix_type="N", &
91 : dist=pao%diag_distribution, &
92 2430 : template=matrix_s(1)%matrix)
93 2430 : CALL dbcsr_reserve_diag_blocks(matrix_M_diag)
94 2430 : CALL dbcsr_complete_redistribute(matrix_M, matrix_M_diag)
95 :
96 : !---------------------------------------------------------------------------
97 : ! cleanup:
98 2430 : CALL dbcsr_release(matrix_M)
99 2430 : CALL dbcsr_release(matrix_Ma)
100 2430 : CALL dbcsr_release(matrix_Mb)
101 2430 : CALL dbcsr_release(matrix_NM)
102 :
103 2430 : CALL timestop(handle)
104 2430 : END SUBROUTINE pao_calc_grad_lnv_wrt_U
105 :
106 : ! **************************************************************************************************
107 : !> \brief Takes current matrix_X and calculates the matrices A and B.
108 : !> \param pao ...
109 : !> \param qs_env ...
110 : !> \param ls_scf_env ...
111 : !> \param matrix_U_diag ...
112 : ! **************************************************************************************************
113 13070 : SUBROUTINE pao_calc_AB_from_U(pao, qs_env, ls_scf_env, matrix_U_diag)
114 : TYPE(pao_env_type), POINTER :: pao
115 : TYPE(qs_environment_type), POINTER :: qs_env
116 : TYPE(ls_scf_env_type), TARGET :: ls_scf_env
117 : TYPE(dbcsr_type) :: matrix_U_diag
118 :
119 : CHARACTER(len=*), PARAMETER :: routineN = 'pao_calc_AB_from_U'
120 :
121 : INTEGER :: acol, arow, handle, iatom
122 : LOGICAL :: found
123 13070 : REAL(dp), DIMENSION(:, :), POINTER :: block_A, block_B, block_N, block_N_inv, &
124 13070 : block_U, block_Y
125 : TYPE(dbcsr_iterator_type) :: iter
126 13070 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
127 : TYPE(dbcsr_type) :: matrix_U
128 : TYPE(ls_mstruct_type), POINTER :: ls_mstruct
129 :
130 13070 : CALL timeset(routineN, handle)
131 13070 : CALL get_qs_env(qs_env, matrix_s=matrix_s)
132 13070 : ls_mstruct => ls_scf_env%ls_mstruct
133 :
134 : ! --------------------------------------------------------------------------------------------
135 : ! sanity check matrix U
136 13070 : CALL pao_assert_unitary(pao, matrix_U_diag)
137 :
138 : ! --------------------------------------------------------------------------------------------
139 : ! redistribute matrix_U_diag from diag_distribution to distribution of matrix_s
140 13070 : CALL get_qs_env(qs_env, matrix_s=matrix_s)
141 13070 : CALL dbcsr_create(matrix_U, matrix_type="N", template=matrix_s(1)%matrix)
142 13070 : CALL dbcsr_reserve_diag_blocks(matrix_U)
143 13070 : CALL dbcsr_complete_redistribute(matrix_U_diag, matrix_U)
144 :
145 : ! --------------------------------------------------------------------------------------------
146 : ! calculate matrix A and B from matrix U
147 : ! Multiplying diagonal matrices is a local operation.
148 : ! To take advantage of this we're using an iterator instead of calling dbcsr_multiply().
149 : !$OMP PARALLEL DEFAULT(NONE) SHARED(pao,ls_mstruct,matrix_U) &
150 13070 : !$OMP PRIVATE(iter,arow,acol,iatom,block_U,block_Y,block_A,block_B,block_N,block_N_inv,found)
151 : CALL dbcsr_iterator_start(iter, matrix_U)
152 : DO WHILE (dbcsr_iterator_blocks_left(iter))
153 : CALL dbcsr_iterator_next_block(iter, arow, acol, block_U)
154 : iatom = arow; CPASSERT(arow == acol)
155 :
156 : CALL dbcsr_get_block_p(matrix=pao%matrix_Y, row=iatom, col=iatom, block=block_Y, found=found)
157 : CPASSERT(ASSOCIATED(block_Y))
158 :
159 : CALL dbcsr_get_block_p(matrix=ls_mstruct%matrix_A, row=iatom, col=iatom, block=block_A, found=found)
160 : CALL dbcsr_get_block_p(matrix=pao%matrix_N_inv, row=iatom, col=iatom, block=block_N_inv, found=found)
161 : CPASSERT(ASSOCIATED(block_A) .AND. ASSOCIATED(block_N_inv))
162 :
163 : CALL dbcsr_get_block_p(matrix=ls_mstruct%matrix_B, row=iatom, col=iatom, block=block_B, found=found)
164 : CALL dbcsr_get_block_p(matrix=pao%matrix_N, row=iatom, col=iatom, block=block_N, found=found)
165 : CPASSERT(ASSOCIATED(block_B) .AND. ASSOCIATED(block_N))
166 :
167 : block_A = MATMUL(MATMUL(block_N_inv, block_U), block_Y)
168 : block_B = MATMUL(MATMUL(block_N, block_U), block_Y)
169 :
170 : END DO
171 : CALL dbcsr_iterator_stop(iter)
172 : !$OMP END PARALLEL
173 :
174 13070 : CALL dbcsr_release(matrix_U)
175 :
176 13070 : CALL timestop(handle)
177 13070 : END SUBROUTINE pao_calc_AB_from_U
178 :
179 : ! **************************************************************************************************
180 : !> \brief Debugging routine, check unitaryness of U
181 : !> \param pao ...
182 : !> \param matrix_U ...
183 : ! **************************************************************************************************
184 18002 : SUBROUTINE pao_assert_unitary(pao, matrix_U)
185 : TYPE(pao_env_type), POINTER :: pao
186 : TYPE(dbcsr_type) :: matrix_U
187 :
188 : CHARACTER(len=*), PARAMETER :: routineN = 'pao_assert_unitary'
189 :
190 : INTEGER :: acol, arow, group_handle, handle, i, &
191 : iatom, M, N
192 13070 : INTEGER, DIMENSION(:), POINTER :: blk_sizes_pao, blk_sizes_pri
193 : REAL(dp) :: delta_max
194 13070 : REAL(dp), DIMENSION(:, :), POINTER :: block_test, tmp1, tmp2
195 : TYPE(dbcsr_iterator_type) :: iter
196 : TYPE(mp_comm_type) :: group
197 :
198 10604 : IF (pao%check_unitary_tol < 0.0_dp) RETURN ! no checking
199 :
200 2466 : CALL timeset(routineN, handle)
201 2466 : delta_max = 0.0_dp
202 :
203 2466 : CALL dbcsr_get_info(pao%matrix_Y, row_blk_size=blk_sizes_pri, col_blk_size=blk_sizes_pao)
204 :
205 : !$OMP PARALLEL DEFAULT(NONE) SHARED(pao,matrix_U,blk_sizes_pri,blk_sizes_pao,delta_max) &
206 2466 : !$OMP PRIVATE(iter,arow,acol,iatom,N,M,block_test,tmp1,tmp2)
207 : CALL dbcsr_iterator_start(iter, matrix_U)
208 : DO WHILE (dbcsr_iterator_blocks_left(iter))
209 : CALL dbcsr_iterator_next_block(iter, arow, acol, block_test)
210 : iatom = arow; CPASSERT(arow == acol)
211 : N = blk_sizes_pri(iatom) ! size of primary basis
212 : M = blk_sizes_pao(iatom) ! size of pao basis
213 :
214 : ! we only need the upper left "PAO-corner" to be unitary
215 : ALLOCATE (tmp1(N, M), tmp2(M, M))
216 : tmp1 = block_test(:, 1:M)
217 : tmp2 = MATMUL(TRANSPOSE(tmp1), tmp1)
218 : DO i = 1, M
219 : tmp2(i, i) = tmp2(i, i) - 1.0_dp
220 : END DO
221 :
222 : !$OMP ATOMIC
223 : delta_max = MAX(delta_max, MAXVAL(ABS(tmp2)))
224 :
225 : DEALLOCATE (tmp1, tmp2)
226 : END DO
227 : CALL dbcsr_iterator_stop(iter)
228 : !$OMP END PARALLEL
229 :
230 2466 : CALL dbcsr_get_info(matrix_U, group=group_handle)
231 2466 : CALL group%set_handle(group_handle)
232 :
233 2466 : CALL group%max(delta_max)
234 2466 : IF (pao%iw > 0) WRITE (pao%iw, *) 'PAO| checked unitaryness, max delta:', delta_max
235 2466 : IF (delta_max > pao%check_unitary_tol) &
236 0 : CPABORT("Found bad unitaryness:"//cp_to_string(delta_max))
237 :
238 2466 : CALL timestop(handle)
239 13070 : END SUBROUTINE pao_assert_unitary
240 :
241 : ! **************************************************************************************************
242 : !> \brief Helper routine, calculates partial derivative dE/dA and dE/dB.
243 : !> As energy functional serves the definition by LNV (Li, Nunes, Vanderbilt).
244 : !> \param qs_env ...
245 : !> \param ls_scf_env ...
246 : !> \param matrix_Ma the derivate wrt A, matrix uses s_matrix-distribution.
247 : !> \param matrix_Mb the derivate wrt B, matrix uses s_matrix-distribution.
248 : ! **************************************************************************************************
249 2662 : SUBROUTINE pao_calc_grad_lnv_wrt_AB(qs_env, ls_scf_env, matrix_Ma, matrix_Mb)
250 : TYPE(qs_environment_type), POINTER :: qs_env
251 : TYPE(ls_scf_env_type), TARGET :: ls_scf_env
252 : TYPE(dbcsr_type) :: matrix_Ma, matrix_Mb
253 :
254 : CHARACTER(len=*), PARAMETER :: routineN = 'pao_calc_grad_lnv_wrt_AB'
255 :
256 : INTEGER :: handle, nspin
257 2662 : INTEGER, DIMENSION(:), POINTER :: pao_blk_sizes
258 : REAL(KIND=dp) :: filter_eps
259 2662 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s, rho_ao
260 : TYPE(dbcsr_type) :: matrix_HB, matrix_HPS, matrix_M, matrix_M1, matrix_M1_dc, matrix_M2, &
261 : matrix_M2_dc, matrix_M3, matrix_M3_dc, matrix_PA, matrix_PH, matrix_PHP, matrix_PSP, &
262 : matrix_SB, matrix_SP
263 : TYPE(dft_control_type), POINTER :: dft_control
264 : TYPE(ls_mstruct_type), POINTER :: ls_mstruct
265 : TYPE(pao_env_type), POINTER :: pao
266 : TYPE(qs_rho_type), POINTER :: rho
267 :
268 2662 : CALL timeset(routineN, handle)
269 :
270 2662 : ls_mstruct => ls_scf_env%ls_mstruct
271 2662 : pao => ls_scf_env%pao_env
272 :
273 : CALL get_qs_env(qs_env, &
274 : rho=rho, &
275 : matrix_ks=matrix_ks, &
276 : matrix_s=matrix_s, &
277 2662 : dft_control=dft_control)
278 2662 : CALL qs_rho_get(rho, rho_ao=rho_ao)
279 2662 : nspin = dft_control%nspins
280 2662 : filter_eps = ls_scf_env%eps_filter
281 :
282 2662 : CALL dbcsr_get_info(ls_mstruct%matrix_A, col_blk_size=pao_blk_sizes)
283 :
284 2662 : IF (nspin /= 1) CPABORT("open shell not yet implemented")
285 : !TODO: handle openshell case properly
286 :
287 : ! Notation according to equation (4.6) on page 50 from:
288 : ! https://dx.doi.org/10.3929%2Fethz-a-010819495
289 :
290 : !---------------------------------------------------------------------------
291 : ! calculate need products in pao basis
292 2662 : CALL dbcsr_create(matrix_PH, template=ls_scf_env%matrix_s, matrix_type="N")
293 : CALL dbcsr_multiply("N", "N", 1.0_dp, ls_scf_env%matrix_p(1), ls_scf_env%matrix_ks(1), &
294 2662 : 0.0_dp, matrix_PH, filter_eps=filter_eps)
295 :
296 2662 : CALL dbcsr_create(matrix_PHP, template=ls_scf_env%matrix_s, matrix_type="N")
297 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_PH, ls_scf_env%matrix_p(1), &
298 2662 : 0.0_dp, matrix_PHP, filter_eps=filter_eps)
299 :
300 2662 : CALL dbcsr_create(matrix_SP, template=ls_scf_env%matrix_s, matrix_type="N")
301 : CALL dbcsr_multiply("N", "N", 1.0_dp, ls_scf_env%matrix_s, ls_scf_env%matrix_p(1), &
302 2662 : 0.0_dp, matrix_SP, filter_eps=filter_eps)
303 :
304 2662 : IF (nspin == 1) CALL dbcsr_scale(matrix_SP, 0.5_dp)
305 :
306 2662 : CALL dbcsr_create(matrix_HPS, template=ls_scf_env%matrix_s, matrix_type="N")
307 : CALL dbcsr_multiply("N", "T", 1.0_dp, ls_scf_env%matrix_ks(1), matrix_SP, &
308 2662 : 0.0_dp, matrix_HPS, filter_eps=filter_eps)
309 :
310 2662 : CALL dbcsr_create(matrix_PSP, template=ls_scf_env%matrix_s, matrix_type="N")
311 : CALL dbcsr_multiply("N", "N", 1.0_dp, ls_scf_env%matrix_p(1), matrix_SP, &
312 2662 : 0.0_dp, matrix_PSP, filter_eps=filter_eps)
313 :
314 : !---------------------------------------------------------------------------
315 : ! M1 = dE_lnv / dP_pao
316 2662 : CALL dbcsr_create(matrix_M1, template=ls_scf_env%matrix_s, matrix_type="N")
317 :
318 : CALL dbcsr_multiply("N", "T", 3.0_dp, ls_scf_env%matrix_ks(1), matrix_SP, &
319 2662 : 1.0_dp, matrix_M1, filter_eps=filter_eps)
320 :
321 : CALL dbcsr_multiply("N", "N", 3.0_dp, matrix_SP, ls_scf_env%matrix_ks(1), &
322 2662 : 1.0_dp, matrix_M1, filter_eps=filter_eps)
323 :
324 : CALL dbcsr_multiply("N", "T", -2.0_dp, matrix_HPS, matrix_SP, &
325 2662 : 1.0_dp, matrix_M1, filter_eps=filter_eps)
326 :
327 : CALL dbcsr_multiply("N", "N", -2.0_dp, matrix_SP, matrix_HPS, &
328 2662 : 1.0_dp, matrix_M1, filter_eps=filter_eps)
329 :
330 : CALL dbcsr_multiply("N", "T", -2.0_dp, matrix_SP, matrix_HPS, &
331 2662 : 1.0_dp, matrix_M1, filter_eps=filter_eps)
332 :
333 : ! reverse possible molecular clustering
334 : CALL dbcsr_create(matrix_M1_dc, &
335 : template=matrix_s(1)%matrix, &
336 : row_blk_size=pao_blk_sizes, &
337 2662 : col_blk_size=pao_blk_sizes)
338 2662 : CALL matrix_decluster(matrix_M1_dc, matrix_M1, ls_mstruct)
339 :
340 : !---------------------------------------------------------------------------
341 : ! M2 = dE_lnv / dH
342 2662 : CALL dbcsr_create(matrix_M2, template=ls_scf_env%matrix_s, matrix_type="N")
343 :
344 2662 : CALL dbcsr_add(matrix_M2, matrix_PSP, 1.0_dp, 3.0_dp)
345 :
346 : CALL dbcsr_multiply("N", "N", -2.0_dp, matrix_PSP, matrix_SP, &
347 2662 : 1.0_dp, matrix_M2, filter_eps=filter_eps)
348 :
349 : ! reverse possible molecular clustering
350 : CALL dbcsr_create(matrix_M2_dc, &
351 : template=matrix_s(1)%matrix, &
352 : row_blk_size=pao_blk_sizes, &
353 2662 : col_blk_size=pao_blk_sizes)
354 2662 : CALL matrix_decluster(matrix_M2_dc, matrix_M2, ls_mstruct)
355 :
356 : !---------------------------------------------------------------------------
357 : ! M3 = dE_lnv / dS
358 2662 : CALL dbcsr_create(matrix_M3, template=ls_scf_env%matrix_s, matrix_type="N")
359 :
360 2662 : CALL dbcsr_add(matrix_M3, matrix_PHP, 1.0_dp, 3.0_dp)
361 :
362 : CALL dbcsr_multiply("N", "N", -2.0_dp, matrix_PHP, matrix_SP, &
363 2662 : 1.0_dp, matrix_M3, filter_eps=filter_eps)
364 :
365 : CALL dbcsr_multiply("N", "T", -2.0_dp, matrix_PSP, matrix_PH, &
366 2662 : 1.0_dp, matrix_M3, filter_eps=filter_eps)
367 :
368 : ! reverse possible molecular clustering
369 : CALL dbcsr_create(matrix_M3_dc, &
370 : template=matrix_s(1)%matrix, &
371 : row_blk_size=pao_blk_sizes, &
372 2662 : col_blk_size=pao_blk_sizes)
373 2662 : CALL matrix_decluster(matrix_M3_dc, matrix_M3, ls_mstruct)
374 :
375 : !---------------------------------------------------------------------------
376 : ! assemble Ma and Mb
377 : ! matrix_Ma = dE_lnv / dA = P * A * M1
378 : ! matrix_Mb = dE_lnv / dB = H * B * M2 + S * B * M3
379 2662 : CALL dbcsr_create(matrix_Ma, template=ls_mstruct%matrix_A, matrix_type="N")
380 2662 : CALL dbcsr_reserve_diag_blocks(matrix_Ma)
381 2662 : CALL dbcsr_create(matrix_Mb, template=ls_mstruct%matrix_B, matrix_type="N")
382 2662 : CALL dbcsr_reserve_diag_blocks(matrix_Mb)
383 :
384 : !---------------------------------------------------------------------------
385 : ! combine M1 with matrices from primary basis
386 2662 : CALL dbcsr_create(matrix_PA, template=ls_mstruct%matrix_A, matrix_type="N")
387 : CALL dbcsr_multiply("N", "N", 1.0_dp, rho_ao(1)%matrix, ls_mstruct%matrix_A, &
388 2662 : 0.0_dp, matrix_PA, filter_eps=filter_eps)
389 :
390 : ! matrix_Ma = P * A * M1
391 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_PA, matrix_M1_dc, &
392 2662 : 0.0_dp, matrix_Ma, filter_eps=filter_eps)
393 :
394 : !---------------------------------------------------------------------------
395 : ! combine M2 with matrices from primary basis
396 2662 : CALL dbcsr_create(matrix_HB, template=ls_mstruct%matrix_B, matrix_type="N")
397 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_ks(1)%matrix, ls_mstruct%matrix_B, &
398 2662 : 0.0_dp, matrix_HB, filter_eps=filter_eps)
399 :
400 : ! matrix_Mb = H * B * M2
401 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_HB, matrix_M2_dc, &
402 2662 : 0.0_dp, matrix_Mb, filter_eps=filter_eps)
403 :
404 : !---------------------------------------------------------------------------
405 : ! combine M3 with matrices from primary basis
406 2662 : CALL dbcsr_create(matrix_SB, template=ls_mstruct%matrix_B, matrix_type="N")
407 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s(1)%matrix, ls_mstruct%matrix_B, &
408 2662 : 0.0_dp, matrix_SB, filter_eps=filter_eps)
409 :
410 2662 : IF (nspin == 1) CALL dbcsr_scale(matrix_SB, 0.5_dp)
411 :
412 : ! matrix_Mb += S * B * M3
413 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_SB, matrix_M3_dc, &
414 2662 : 1.0_dp, matrix_Mb, filter_eps=filter_eps)
415 :
416 2662 : IF (nspin == 1) CALL dbcsr_scale(matrix_Ma, 2.0_dp)
417 2662 : IF (nspin == 1) CALL dbcsr_scale(matrix_Mb, 2.0_dp)
418 :
419 : !---------------------------------------------------------------------------
420 : ! cleanup: TODO release matrices as early as possible
421 2662 : CALL dbcsr_release(matrix_PH)
422 2662 : CALL dbcsr_release(matrix_PHP)
423 2662 : CALL dbcsr_release(matrix_SP)
424 2662 : CALL dbcsr_release(matrix_HPS)
425 2662 : CALL dbcsr_release(matrix_PSP)
426 2662 : CALL dbcsr_release(matrix_M)
427 2662 : CALL dbcsr_release(matrix_M1)
428 2662 : CALL dbcsr_release(matrix_M2)
429 2662 : CALL dbcsr_release(matrix_M3)
430 2662 : CALL dbcsr_release(matrix_M1_dc)
431 2662 : CALL dbcsr_release(matrix_M2_dc)
432 2662 : CALL dbcsr_release(matrix_M3_dc)
433 2662 : CALL dbcsr_release(matrix_PA)
434 2662 : CALL dbcsr_release(matrix_HB)
435 2662 : CALL dbcsr_release(matrix_SB)
436 :
437 2662 : CALL timestop(handle)
438 2662 : END SUBROUTINE pao_calc_grad_lnv_wrt_AB
439 :
440 : END MODULE pao_param_methods
|
