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 : MODULE dkh_main
8 : USE cp_blacs_env, ONLY: cp_blacs_env_type
9 : USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,&
10 : cp_fm_scale_and_add,&
11 : cp_fm_schur_product,&
12 : cp_fm_syrk,&
13 : cp_fm_transpose,&
14 : cp_fm_triangular_multiply,&
15 : cp_fm_uplo_to_full
16 : USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose,&
17 : cp_fm_cholesky_reduce
18 : USE cp_fm_diag, ONLY: cp_fm_syevd
19 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
20 : cp_fm_struct_get,&
21 : cp_fm_struct_release,&
22 : cp_fm_struct_type
23 : USE cp_fm_types, ONLY: cp_fm_create,&
24 : cp_fm_get_info,&
25 : cp_fm_release,&
26 : cp_fm_to_fm,&
27 : cp_fm_type
28 : USE kinds, ONLY: dp
29 : USE parallel_gemm_api, ONLY: parallel_gemm
30 : USE physcon, ONLY: c_light_au
31 : USE qs_environment_types, ONLY: get_qs_env,&
32 : qs_environment_type
33 : #include "./base/base_uses.f90"
34 :
35 : IMPLICIT NONE
36 : PRIVATE
37 :
38 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'dkh_main'
39 :
40 : PUBLIC :: dkh_atom_transformation
41 :
42 : CONTAINS
43 :
44 : ! **************************************************************************************************
45 : !> \brief 2th order DKH calculations
46 : !>
47 : !> \param qs_env ...
48 : !> \param matrix_s ...
49 : !> \param matrix_v ...
50 : !> \param matrix_t ...
51 : !> \param matrix_pVp ...
52 : !> \param n ...
53 : !> \param dkh_order ...
54 : !> \par Literature
55 : !> M. Reiher, A. Wolf, J. Chem. Phys. 121 (2004) 10944-10956
56 : !> A. Wolf, M. Reiher, B. A. Hess, J. Chem. Phys. 117 (2002) 9215-9226
57 : !>
58 : !>\par Note
59 : !> based on subroutines for DKH1 to DKH5 by
60 : !> A. Wolf, M. Reiher, B. A. Hess
61 : !>
62 : !> INPUT:
63 : !> qs_env (:) The quickstep environment
64 : !> n Number of primitive gaussians
65 : !> matrix_s (:,:) overlap matrix
66 : !> matrix_pVp (:,:) pVp matrix
67 : !>
68 : !> IN_OUT:
69 : !> matrix_v (:,:) input: nonrelativistic potential energy matrix
70 : !> output: (ev1+ev2)
71 : !> matrix_t (:,:) input: kinetic energy matrix
72 : !> output: kinetic part of hamiltonian
73 : !> in position space
74 : !>
75 : !> INTERNAL
76 : !> sinv (:,:) inverted, orthogonalized overlap matrix
77 : !> ev0t (:) DKH-even0 matrix in T-basis
78 : !> e (:) e=SQRT(p^2c^2+c^4)
79 : !> eig (:,:) eigenvectors of sinv' h sinv
80 : !> tt (:) eigenvalues of sinv' h sinv
81 : !> revt (:,:) reverse transformation matrix T-basis -> position space
82 : !> aa (:) kinematical factors f. DKH SQRT((c^2+e(i))/(2.0*e(i)))
83 : !> rr (:) kinematical factors f. DKH c/(c^2+e(i))
84 : !> vt (:,:) non relativistic potential matrix in T-basis
85 : !> pvpt (:,:) pvp integral matrix in T-basis
86 : !> ev1t (:,:) DKH-even1 matrix in T-basis
87 : !> evt2 (:,:) DKH-even2 matrix in T-basis
88 : !> ev1 (:,:) DKH-even1 matrix in position space
89 : !> ev2 (:,:) DKH-even2 matrix in position space
90 : !> ove (:,:) scratch
91 : !> aux (:,:) scratch
92 : !> c_light_au velocity of light 137 a.u.
93 : !> prea prefactor, 1/137^2
94 : !> con2 prefactor, 2/137^2
95 : !> con prefactor, 137^2
96 : !> \author
97 : !> Jens Thar, Barbara Kirchner: Uni Bonn (09/2006)
98 : !> Markus Reiher: ETH Zurich (09/2006)
99 : !>
100 : ! **************************************************************************************************
101 0 : SUBROUTINE DKH_full_transformation(qs_env, matrix_s, matrix_v, matrix_t, matrix_pVp, n, dkh_order)
102 : TYPE(qs_environment_type), POINTER :: qs_env
103 : TYPE(cp_fm_type), INTENT(IN) :: matrix_s, matrix_v, matrix_t, matrix_pVp
104 : INTEGER, INTENT(IN) :: n, dkh_order
105 :
106 : CHARACTER(LEN=*), PARAMETER :: routineN = 'DKH_full_transformation'
107 :
108 : INTEGER :: handle, i
109 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: aa, e, ev0t, rr, tt
110 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
111 : TYPE(cp_fm_struct_type), POINTER :: matrix_full
112 : TYPE(cp_fm_type) :: matrix_aux, matrix_aux2, matrix_eig, matrix_ev1, matrix_ev2, matrix_ev3, &
113 : matrix_ev4, matrix_pe1p, matrix_rev, matrix_se, matrix_sinv
114 :
115 0 : CALL timeset(routineN, handle)
116 0 : NULLIFY (blacs_env)
117 :
118 : !-----------------------------------------------------------------------
119 : ! Construct the matrix structure
120 : !-----------------------------------------------------------------------
121 :
122 0 : CALL get_qs_env(qs_env, blacs_env=blacs_env)
123 : CALL cp_fm_struct_create(fmstruct=matrix_full, &
124 : context=blacs_env, &
125 : nrow_global=n, &
126 0 : ncol_global=n)
127 :
128 : !-----------------------------------------------------------------------
129 : ! Allocate some matrices
130 : !-----------------------------------------------------------------------
131 :
132 0 : ALLOCATE (e(n))
133 0 : ALLOCATE (aa(n))
134 0 : ALLOCATE (rr(n))
135 0 : ALLOCATE (tt(n))
136 0 : ALLOCATE (ev0t(n))
137 :
138 0 : CALL cp_fm_create(matrix_eig, matrix_full)
139 0 : CALL cp_fm_create(matrix_aux, matrix_full)
140 0 : CALL cp_fm_create(matrix_aux2, matrix_full)
141 0 : CALL cp_fm_create(matrix_rev, matrix_full)
142 0 : CALL cp_fm_create(matrix_se, matrix_full)
143 0 : CALL cp_fm_create(matrix_ev1, matrix_full)
144 0 : CALL cp_fm_create(matrix_ev2, matrix_full)
145 0 : CALL cp_fm_create(matrix_sinv, matrix_full)
146 0 : CALL cp_fm_create(matrix_ev3, matrix_full)
147 0 : CALL cp_fm_create(matrix_ev4, matrix_full)
148 0 : CALL cp_fm_create(matrix_pe1p, matrix_full)
149 :
150 : !-----------------------------------------------------------------------
151 : ! Now with Cholesky decomposition
152 : !-----------------------------------------------------------------------
153 :
154 0 : CALL cp_fm_to_fm(matrix_s, matrix_sinv)
155 0 : CALL cp_fm_cholesky_decompose(matrix_sinv, n)
156 :
157 : !-----------------------------------------------------------------------
158 : ! Calculate matrix representation from nonrelativistic T matrix
159 : !-----------------------------------------------------------------------
160 :
161 0 : CALL cp_fm_cholesky_reduce(matrix_t, matrix_sinv)
162 0 : CALL cp_fm_syevd(matrix_t, matrix_eig, tt)
163 :
164 : !-----------------------------------------------------------------------
165 : ! Calculate kinetic part of Hamiltonian in T-basis
166 : !-----------------------------------------------------------------------
167 :
168 0 : CALL kintegral(n, ev0t, tt, e)
169 :
170 : !-----------------------------------------------------------------------
171 : ! Calculate reverse transformation matrix revt
172 : !-----------------------------------------------------------------------
173 :
174 0 : CALL cp_fm_to_fm(matrix_eig, matrix_rev)
175 0 : CALL cp_fm_triangular_multiply(matrix_sinv, matrix_rev, transpose_tr=.TRUE.)
176 :
177 : !-----------------------------------------------------------------------
178 : ! Calculate kinetic part of the Hamiltonian
179 : !-----------------------------------------------------------------------
180 :
181 0 : CALL cp_fm_to_fm(matrix_rev, matrix_aux)
182 0 : CALL cp_fm_column_scale(matrix_aux, ev0t)
183 0 : CALL parallel_gemm("N", "T", n, n, n, 1.0_dp, matrix_rev, matrix_aux, 0.0_dp, matrix_t)
184 :
185 : !-----------------------------------------------------------------------
186 : ! Calculate kinematical factors for DKH
187 : ! only vectors present - will be done on every CPU
188 : !-----------------------------------------------------------------------
189 :
190 0 : DO i = 1, n
191 0 : aa(i) = SQRT((c_light_au**2 + e(i))/(2.0_dp*e(i)))
192 0 : rr(i) = SQRT(c_light_au**2)/(c_light_au**2 + e(i))
193 : END DO
194 :
195 : !-----------------------------------------------------------------------
196 : ! Transform v integrals to T-basis (v -> v(t))
197 : !-----------------------------------------------------------------------
198 :
199 0 : CALL cp_fm_cholesky_reduce(matrix_v, matrix_sinv)
200 0 : CALL cp_fm_uplo_to_full(matrix_v, matrix_aux)
201 0 : CALL parallel_gemm("T", "N", n, n, n, 1.0_dp, matrix_eig, matrix_v, 0.0_dp, matrix_aux)
202 0 : CALL parallel_gemm("N", "N", n, n, n, 1.0_dp, matrix_aux, matrix_eig, 0.0_dp, matrix_v)
203 :
204 : !-----------------------------------------------------------------------
205 : ! Transform pVp integrals to T-basis (pVp -> pVp(t))
206 : !-----------------------------------------------------------------------
207 :
208 0 : CALL cp_fm_cholesky_reduce(matrix_pVp, matrix_sinv)
209 0 : CALL cp_fm_uplo_to_full(matrix_pVp, matrix_aux)
210 0 : CALL parallel_gemm("T", "N", n, n, n, 1.0_dp, matrix_eig, matrix_pVp, 0.0_dp, matrix_aux)
211 0 : CALL parallel_gemm("N", "N", n, n, n, 1.0_dp, matrix_aux, matrix_eig, 0.0_dp, matrix_pVp)
212 :
213 : !-----------------------------------------------------------------------
214 : ! Calculate even1 in T-basis
215 : !-----------------------------------------------------------------------
216 :
217 0 : CALL even1(matrix_ev1, matrix_v, matrix_pvp, aa, rr, matrix_aux, matrix_aux2)
218 :
219 : !-----------------------------------------------------------------------
220 : ! Calculate even2 in T-basis
221 : !-----------------------------------------------------------------------
222 :
223 0 : CALL even2c(n, matrix_ev2, matrix_v, matrix_pVp, aa, rr, tt, e, matrix_aux)
224 :
225 : !-----------------------------------------------------------------------
226 : ! Calculate even3 in T-basis, only if requested
227 : !-----------------------------------------------------------------------
228 :
229 0 : IF (dkh_order >= 3) THEN
230 0 : CALL peven1p(n, matrix_pe1p, matrix_v, matrix_pvp, matrix_aux, matrix_aux2, aa, rr, tt)
231 0 : CALL even3b(n, matrix_ev3, matrix_ev1, matrix_pe1p, matrix_v, matrix_pvp, aa, rr, tt, e, matrix_aux)
232 :
233 : !-----------------------------------------------------------------------
234 : ! Transform even3 back to position space
235 : !-----------------------------------------------------------------------
236 :
237 0 : CALL parallel_gemm("N", "N", n, n, n, 1.0_dp, matrix_rev, matrix_ev3, 0.0_dp, matrix_aux)
238 0 : CALL parallel_gemm("N", "T", n, n, n, 1.0_dp, matrix_aux, matrix_rev, 0.0_dp, matrix_ev3)
239 :
240 : !-----------------------------------------------------------------------
241 : ! Calculate even4 in T-basis, only if requested
242 : !-----------------------------------------------------------------------
243 :
244 0 : IF (dkh_order >= 4) THEN
245 0 : CPABORT("DKH order greater than 3 not yet available")
246 : ! CALL even4a(n,matrix_ev4%local_data,matrix_ev2%local_data,matrix_pe1p%local_data,matrix_v%local_data,&
247 : ! matrix_pvp%local_data,aa,rr,tt,e)
248 :
249 : !-----------------------------------------------------------------------
250 : ! Transform even4 back to position space
251 : !-----------------------------------------------------------------------
252 :
253 : ! CALL parallel_gemm("N","N",n,n,n,1.0_dp,matrix_rev,matrix_ev4,0.0_dp,matrix_aux)
254 : ! CALL parallel_gemm("N","T",n,n,n,1.0_dp,matrix_aux,matrix_rev,0.0_dp,matrix_ev4)
255 :
256 : END IF
257 : END IF
258 :
259 : !----------------------------------------------------------------------
260 : ! Transform even1 back to position space
261 : !----------------------------------------------------------------------
262 :
263 0 : CALL parallel_gemm("N", "N", n, n, n, 1.0_dp, matrix_rev, matrix_ev1, 0.0_dp, matrix_aux)
264 0 : CALL parallel_gemm("N", "T", n, n, n, 1.0_dp, matrix_aux, matrix_rev, 0.0_dp, matrix_ev1)
265 :
266 : !-----------------------------------------------------------------------
267 : ! Transform even2 back to position space
268 : !-----------------------------------------------------------------------
269 :
270 0 : CALL parallel_gemm("N", "N", n, n, n, 1.0_dp, matrix_rev, matrix_ev2, 0.0_dp, matrix_aux)
271 0 : CALL parallel_gemm("N", "T", n, n, n, 1.0_dp, matrix_aux, matrix_rev, 0.0_dp, matrix_ev2)
272 :
273 : !-----------------------------------------------------------------------
274 : ! Calculate v in position space
275 : !-----------------------------------------------------------------------
276 :
277 0 : CALL cp_fm_scale_and_add(1.0_dp, matrix_ev1, 1.0_dp, matrix_ev2)
278 0 : CALL cp_fm_uplo_to_full(matrix_ev1, matrix_aux)
279 0 : CALL cp_fm_to_fm(matrix_ev1, matrix_v)
280 0 : IF (dkh_order >= 3) THEN
281 0 : CALL cp_fm_scale_and_add(1.0_dp, matrix_v, 1.0_dp, matrix_ev3)
282 0 : IF (dkh_order >= 4) THEN
283 0 : CALL cp_fm_scale_and_add(1.0_dp, matrix_v, 1.0_dp, matrix_ev4)
284 : END IF
285 : END IF
286 :
287 0 : CALL cp_fm_release(matrix_eig)
288 0 : CALL cp_fm_release(matrix_aux)
289 0 : CALL cp_fm_release(matrix_aux2)
290 0 : CALL cp_fm_release(matrix_rev)
291 0 : CALL cp_fm_release(matrix_se)
292 0 : CALL cp_fm_release(matrix_ev1)
293 0 : CALL cp_fm_release(matrix_ev2)
294 0 : CALL cp_fm_release(matrix_sinv)
295 0 : CALL cp_fm_release(matrix_ev3)
296 0 : CALL cp_fm_release(matrix_ev4)
297 0 : CALL cp_fm_release(matrix_pe1p)
298 :
299 0 : CALL cp_fm_struct_release(matrix_full)
300 :
301 0 : DEALLOCATE (ev0t, e, aa, rr, tt)
302 :
303 0 : CALL timestop(handle)
304 :
305 0 : END SUBROUTINE DKH_full_transformation
306 :
307 : ! **************************************************************************************************
308 : !> \brief ...
309 : !> \param n ...
310 : !> \param ev0t ...
311 : !> \param tt ...
312 : !> \param e ...
313 : ! **************************************************************************************************
314 0 : SUBROUTINE kintegral(n, ev0t, tt, e)
315 : INTEGER, INTENT(IN) :: n
316 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: ev0t
317 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: tt
318 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: e
319 :
320 : INTEGER :: i
321 : REAL(KIND=dp) :: con, con2, prea, ratio, tv1, tv2, tv3, &
322 : tv4
323 :
324 0 : prea = 1/(c_light_au**2)
325 0 : con2 = prea + prea
326 0 : con = 1.0_dp/prea
327 :
328 0 : DO i = 1, n
329 0 : IF (tt(i) < 0.0_dp) THEN
330 0 : WRITE (*, *) ' dkh_main.F | tt(', i, ') = ', tt(i)
331 : END IF
332 :
333 : ! If T is sufficiently small, use series expansion to avoid
334 : ! cancellation, otherwise calculate SQRT directly
335 :
336 0 : ev0t(i) = tt(i)
337 0 : ratio = tt(i)/c_light_au
338 0 : IF (ratio <= 0.02_dp) THEN
339 0 : tv1 = tt(i)
340 0 : tv2 = -tv1*tt(i)*prea*0.5_dp
341 0 : tv3 = -tv2*tt(i)*prea
342 0 : tv4 = -tv3*tt(i)*prea*1.25_dp
343 0 : ev0t(i) = tv1 + tv2 + tv3 + tv4
344 : ELSE
345 0 : ev0t(i) = con*(SQRT(1.0_dp + con2*tt(i)) - 1.0_dp)
346 : END IF
347 0 : e(i) = ev0t(i) + con
348 : END DO
349 :
350 0 : END SUBROUTINE kintegral
351 :
352 : ! **************************************************************************************************
353 : !> \brief 1st order DKH-approximation
354 : !> \param matrix_ev1 even1 output matrix
355 : !> \param matrix_v potential matrix v in T-space
356 : !> \param matrix_pvp pvp matrix in T-space
357 : !> \param aa A-factors (diagonal)
358 : !> \param rr R-factors (diagonal)
359 : !> \param matrix_aux ...
360 : !> \param matrix_aux2 ...
361 : ! **************************************************************************************************
362 0 : SUBROUTINE even1(matrix_ev1, matrix_v, matrix_pvp, aa, rr, matrix_aux, matrix_aux2)
363 : TYPE(cp_fm_type), INTENT(IN) :: matrix_ev1, matrix_v, matrix_pVp
364 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: aa, rr
365 : TYPE(cp_fm_type), INTENT(IN) :: matrix_aux, matrix_aux2
366 :
367 0 : CALL cp_fm_to_fm(matrix_v, matrix_aux)
368 0 : CALL cp_fm_column_scale(matrix_aux, aa)
369 0 : CALL cp_fm_transpose(matrix_aux, matrix_ev1)
370 0 : CALL cp_fm_column_scale(matrix_ev1, aa)
371 :
372 0 : CALL cp_fm_to_fm(matrix_pVp, matrix_aux)
373 0 : CALL cp_fm_column_scale(matrix_aux, aa)
374 0 : CALL cp_fm_column_scale(matrix_aux, rr)
375 0 : CALL cp_fm_transpose(matrix_aux, matrix_aux2)
376 0 : CALL cp_fm_column_scale(matrix_aux2, aa)
377 0 : CALL cp_fm_column_scale(matrix_aux2, rr)
378 :
379 0 : CALL cp_fm_scale_and_add(1.0_dp, matrix_ev1, 1.0_dp, matrix_aux2)
380 :
381 0 : END SUBROUTINE even1
382 :
383 : ! **************************************************************************************************
384 : !> \brief 1st order DKH-approximation
385 : !> \param n dimension of matrices
386 : !> \param matrix_pe1p peven1p output matrix
387 : !> \param matrix_v potential matrix v in T-space
388 : !> \param matrix_pvp pvp matrix in T-space
389 : !> \param matrix_aux ...
390 : !> \param matrix_aux2 ...
391 : !> \param aa A-factors (diagonal)
392 : !> \param rr R-factors (diagonal)
393 : !> \param tt ...
394 : ! **************************************************************************************************
395 0 : SUBROUTINE peven1p(n, matrix_pe1p, matrix_v, matrix_pvp, matrix_aux, matrix_aux2, aa, rr, tt)
396 :
397 : INTEGER, INTENT(IN) :: n
398 : TYPE(cp_fm_type), INTENT(IN) :: matrix_pe1p, matrix_v, matrix_pvp, &
399 : matrix_aux, matrix_aux2
400 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: aa, rr, tt
401 :
402 : INTEGER :: i, nrow_local
403 0 : INTEGER, DIMENSION(:), POINTER :: row_indices
404 0 : REAL(KIND=dp), DIMENSION(n) :: vec_ar, vec_arrt
405 : TYPE(cp_blacs_env_type), POINTER :: context
406 : TYPE(cp_fm_struct_type), POINTER :: vec_full
407 : TYPE(cp_fm_type) :: vec_a
408 :
409 0 : DO i = 1, n
410 0 : vec_ar(i) = aa(i)*rr(i)
411 0 : vec_arrt(i) = vec_ar(i)*rr(i)*tt(i)
412 : END DO
413 :
414 0 : CALL cp_fm_struct_get(matrix_v%matrix_struct, context=context)
415 : CALL cp_fm_struct_create(fmstruct=vec_full, &
416 : context=context, &
417 : nrow_global=n, &
418 0 : ncol_global=1)
419 :
420 0 : CALL cp_fm_create(vec_a, vec_full)
421 :
422 : CALL cp_fm_get_info(matrix_v, nrow_local=nrow_local, &
423 0 : row_indices=row_indices)
424 :
425 0 : DO i = 1, nrow_local
426 0 : vec_a%local_data(i, 1) = vec_arrt(row_indices(i))
427 : END DO
428 :
429 0 : CALL cp_fm_syrk('U', 'N', 1, 1.0_dp, vec_a, 1, 1, 0.0_dp, matrix_aux)
430 0 : CALL cp_fm_uplo_to_full(matrix_aux, matrix_aux2)
431 0 : CALL cp_fm_schur_product(matrix_v, matrix_aux, matrix_pe1p)
432 :
433 0 : DO i = 1, nrow_local
434 0 : vec_a%local_data(i, 1) = vec_ar(row_indices(i))
435 : END DO
436 :
437 0 : CALL cp_fm_syrk('U', 'N', 1, 1.0_dp, vec_a, 1, 1, 0.0_dp, matrix_aux)
438 0 : CALL cp_fm_uplo_to_full(matrix_aux, matrix_aux2)
439 0 : CALL cp_fm_schur_product(matrix_pvp, matrix_aux, matrix_aux2)
440 :
441 0 : CALL cp_fm_scale_and_add(4.0_dp, matrix_pe1p, 1.0_dp, matrix_aux2)
442 :
443 0 : CALL cp_fm_release(vec_a)
444 0 : CALL cp_fm_struct_release(vec_full)
445 :
446 0 : END SUBROUTINE peven1p
447 :
448 : ! **************************************************************************************************
449 : !> \brief 2nd order DK-approximation (original DK-transformation with U = SQRT(1+W^2) + W)
450 : !> \param n Dimension of matrices
451 : !> \param matrix_ev2 even2 output matrix = final result
452 : !> \param matrix_v symmetric (n x n)-matrix containing (A V A)
453 : !> \param matrix_pVp symmetric (n x n)-matrix containing (A P V P A)
454 : !> \param aa A-Factors (DIAGONAL
455 : !> \param rr R-Factors (DIAGONAL)
456 : !> \param tt Nonrel. kinetic Energy (DIAGONAL)
457 : !> \param e Rel. Energy = SQRT(p^2*c^2 + c^4) (DIAGONAL)
458 : !> \param matrix_aux ...
459 : ! **************************************************************************************************
460 0 : SUBROUTINE even2c(n, matrix_ev2, matrix_v, matrix_pVp, aa, rr, tt, e, matrix_aux)
461 :
462 : !***********************************************************************
463 : ! *
464 : ! Alexander Wolf, last modified: 20.02.2002 - DKH2 *
465 : ! *
466 : ! *
467 : ! Version: 1.1 (20.2.2002) : Usage of SR mat_add included *
468 : ! 1.0 (6.2.2002) *
469 : ! Modification history: *
470 : ! 30.09.2006 Jens Thar: deleted obsolete F77 memory manager *
471 : ! 2008 Jens Thar: transfer to CP2K *
472 : ! *
473 : ! ev2 = 1/2 [W1,O1] *
474 : ! *
475 : !***********************************************************************
476 :
477 : INTEGER, INTENT(IN) :: n
478 : TYPE(cp_fm_type), INTENT(IN) :: matrix_ev2, matrix_v, matrix_pVp
479 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: aa, rr, tt, e
480 : TYPE(cp_fm_type), INTENT(IN) :: matrix_aux
481 :
482 : TYPE(cp_blacs_env_type), POINTER :: context
483 : TYPE(cp_fm_struct_type), POINTER :: matrix_full
484 : TYPE(cp_fm_type) :: matrix_apVpa, matrix_apVVpa, &
485 : matrix_aux2, matrix_ava, matrix_avva
486 :
487 : ! result intermediate result of even2-calculation
488 : !-----------------------------------------------------------------------
489 : ! 1. General Structures and Patterns for DKH2
490 : !-----------------------------------------------------------------------
491 :
492 0 : CALL cp_fm_struct_get(matrix_v%matrix_struct, context=context)
493 : CALL cp_fm_struct_create(fmstruct=matrix_full, &
494 : context=context, &
495 : nrow_global=n, &
496 0 : ncol_global=n)
497 :
498 0 : CALL cp_fm_create(matrix_aux2, matrix_full)
499 0 : CALL cp_fm_create(matrix_ava, matrix_full)
500 0 : CALL cp_fm_create(matrix_avva, matrix_full)
501 0 : CALL cp_fm_create(matrix_apVpa, matrix_full)
502 0 : CALL cp_fm_create(matrix_apVVpa, matrix_full)
503 :
504 0 : CALL cp_fm_to_fm(matrix_v, matrix_ava)
505 0 : CALL cp_fm_to_fm(matrix_v, matrix_avva)
506 0 : CALL cp_fm_to_fm(matrix_pVp, matrix_apVpa)
507 0 : CALL cp_fm_to_fm(matrix_pVp, matrix_apVVpa)
508 :
509 : ! Calculate v = A V A:
510 :
511 0 : CALL mat_axa(matrix_v, matrix_ava, n, aa, matrix_aux)
512 :
513 : ! Calculate pvp = A P V P A:
514 :
515 0 : CALL mat_arxra(matrix_pVp, matrix_apVpa, n, aa, rr, matrix_aux)
516 :
517 : ! Calculate vh = A V~ A:
518 :
519 0 : CALL mat_1_over_h(matrix_v, matrix_avva, e, matrix_aux)
520 0 : CALL cp_fm_to_fm(matrix_avva, matrix_aux2)
521 0 : CALL mat_axa(matrix_aux2, matrix_avva, n, aa, matrix_aux)
522 :
523 : ! Calculate pvph = A P V~ P A:
524 :
525 0 : CALL mat_1_over_h(matrix_pVp, matrix_apVVpa, e, matrix_aux)
526 0 : CALL cp_fm_to_fm(matrix_apVVpa, matrix_aux2)
527 0 : CALL mat_arxra(matrix_aux2, matrix_apVVpa, n, aa, rr, matrix_aux)
528 :
529 : ! Calculate w1o1:
530 :
531 0 : CALL parallel_gemm("N", "N", n, n, n, -1.0_dp, matrix_apVVpa, matrix_ava, 0.0_dp, matrix_aux2)
532 0 : CALL mat_muld(matrix_aux2, matrix_apVVpa, matrix_apVpa, n, 1.0_dp, 1.0_dp, tt, rr, matrix_aux)
533 0 : CALL mat_mulm(matrix_aux2, matrix_avva, matrix_ava, n, 1.0_dp, 1.0_dp, tt, rr, matrix_aux)
534 0 : CALL parallel_gemm("N", "N", n, n, n, -1.0_dp, matrix_avva, matrix_apVpa, 1.0_dp, matrix_aux2)
535 :
536 : ! Calculate o1w1 (already stored in ev2):
537 :
538 0 : CALL parallel_gemm("N", "N", n, n, n, 1.0_dp, matrix_apVpa, matrix_avva, 0.0_dp, matrix_ev2)
539 0 : CALL mat_muld(matrix_ev2, matrix_apVpa, matrix_apVVpa, n, -1.0_dp, 1.0_dp, tt, rr, matrix_aux)
540 0 : CALL mat_mulm(matrix_ev2, matrix_ava, matrix_avva, n, -1.0_dp, 1.0_dp, tt, rr, matrix_aux)
541 0 : CALL parallel_gemm("N", "N", n, n, n, 1.0_dp, matrix_ava, matrix_apVVpa, 1.0_dp, matrix_ev2)
542 :
543 : !-----------------------------------------------------------------------
544 : ! 2. 1/2 [W1,O1] = 1/2 W1O1 - 1/2 O1W1
545 : !-----------------------------------------------------------------------
546 :
547 0 : CALL cp_fm_scale_and_add(-0.5_dp, matrix_ev2, 0.5_dp, matrix_aux2)
548 :
549 : !-----------------------------------------------------------------------
550 : ! 3. Finish up the stuff!!
551 : !-----------------------------------------------------------------------
552 :
553 0 : CALL cp_fm_release(matrix_aux2)
554 0 : CALL cp_fm_release(matrix_ava)
555 0 : CALL cp_fm_release(matrix_avva)
556 0 : CALL cp_fm_release(matrix_apVpa)
557 0 : CALL cp_fm_release(matrix_apVVpa)
558 :
559 0 : CALL cp_fm_struct_release(matrix_full)
560 :
561 : ! WRITE (*,*) "CAW: DKH2 with even2c (Alex)"
562 : ! WRITE (*,*) "JT: Now available in cp2k"
563 :
564 0 : END SUBROUTINE even2c
565 :
566 : ! **************************************************************************************************
567 : !> \brief ...
568 : !> \param n ...
569 : !> \param matrix_ev3 ...
570 : !> \param matrix_ev1 ...
571 : !> \param matrix_pe1p ...
572 : !> \param matrix_v ...
573 : !> \param matrix_pVp ...
574 : !> \param aa ...
575 : !> \param rr ...
576 : !> \param tt ...
577 : !> \param e ...
578 : !> \param matrix_aux ...
579 : ! **************************************************************************************************
580 0 : SUBROUTINE even3b(n, matrix_ev3, matrix_ev1, matrix_pe1p, matrix_v, matrix_pVp, aa, rr, tt, e, matrix_aux)
581 :
582 : !***********************************************************************
583 : ! *
584 : ! Alexander Wolf, last modified: 20.2.2002 - DKH3 *
585 : ! *
586 : ! 3rd order DK-approximation (generalised DK-transformation) *
587 : ! *
588 : ! Version: 1.1 (20.2.2002) : Usage of SR mat_add included *
589 : ! 1.0 (7.2.2002) *
590 : ! *
591 : ! ev3 = 1/2 [W1,[W1,E1]] *
592 : ! *
593 : ! Modification history: *
594 : ! 30.09.2006 Jens Thar: deleted obsolete F77 memory manager *
595 : ! *
596 : ! ---- Meaning of Parameters ---- *
597 : ! *
598 : ! n in Dimension of matrices *
599 : ! ev3 out even3 output matrix = final result *
600 : ! e1 in E1 = even1-operator *
601 : ! pe1p in pE1p *
602 : ! vv in potential v *
603 : ! gg in pvp *
604 : ! aa in A-Factors (DIAGONAL) *
605 : ! rr in R-Factors (DIAGONAL) *
606 : ! tt in Nonrel. kinetic Energy (DIAGONAL) *
607 : ! e in Rel. Energy = SQRT(p^2*c^2 + c^4) (DIAGONAL) *
608 : ! result intermediate result of even2-calculation
609 : ! vh symmetric (n x n)-matrix containing (A V~ A)
610 : ! pvph symmetric (n x n)-matrix containing (A P V~ P A)
611 : ! e1 E1
612 : ! pe1p pE1p
613 : ! w1w1 (W1)^2
614 : ! w1e1w1 W1*E1*W1
615 : ! scr_i temporary (n x n)-scratch-matrices (i=1,2)
616 : ! *
617 : !***********************************************************************
618 :
619 : INTEGER, INTENT(IN) :: n
620 : TYPE(cp_fm_type), INTENT(IN) :: matrix_ev3, matrix_ev1, matrix_pe1p, &
621 : matrix_v, matrix_pVp
622 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: aa, rr, tt, e
623 : TYPE(cp_fm_type), INTENT(IN) :: matrix_aux
624 :
625 : TYPE(cp_blacs_env_type), POINTER :: context
626 : TYPE(cp_fm_struct_type), POINTER :: matrix_full
627 : TYPE(cp_fm_type) :: matrix_apVVpa, matrix_aux2, matrix_avva, &
628 : matrix_w1e1w1, matrix_w1w1
629 :
630 : !-----------------------------------------------------------------------
631 : ! 1. General Structures and Patterns for DKH3
632 : !-----------------------------------------------------------------------
633 :
634 0 : CALL cp_fm_struct_get(matrix_v%matrix_struct, context=context)
635 : CALL cp_fm_struct_create(fmstruct=matrix_full, &
636 : context=context, &
637 : nrow_global=n, &
638 0 : ncol_global=n)
639 :
640 0 : CALL cp_fm_create(matrix_aux2, matrix_full)
641 0 : CALL cp_fm_create(matrix_w1w1, matrix_full)
642 0 : CALL cp_fm_create(matrix_w1e1w1, matrix_full)
643 0 : CALL cp_fm_create(matrix_avva, matrix_full)
644 0 : CALL cp_fm_create(matrix_apVVpa, matrix_full)
645 :
646 0 : CALL cp_fm_to_fm(matrix_v, matrix_avva)
647 0 : CALL cp_fm_to_fm(matrix_pVp, matrix_apVVpa)
648 :
649 : ! Calculate vh = A V~ A:
650 :
651 0 : CALL mat_1_over_h(matrix_v, matrix_avva, e, matrix_aux)
652 0 : CALL cp_fm_to_fm(matrix_avva, matrix_aux2)
653 0 : CALL mat_axa(matrix_aux2, matrix_avva, n, aa, matrix_aux)
654 :
655 : ! Calculate pvph = A P V~ P A:
656 :
657 0 : CALL mat_1_over_h(matrix_pVp, matrix_apVVpa, e, matrix_aux)
658 0 : CALL cp_fm_to_fm(matrix_apVVpa, matrix_aux2)
659 0 : CALL mat_arxra(matrix_aux2, matrix_apVVpa, n, aa, rr, matrix_aux)
660 :
661 : ! Calculate w1w1:
662 :
663 0 : CALL parallel_gemm("N", "N", n, n, n, 1.0_dp, matrix_apVVpa, matrix_avva, 0.0_dp, matrix_w1w1)
664 0 : CALL mat_muld(matrix_w1w1, matrix_apVVpa, matrix_apVVpa, n, -1.0_dp, 1.0_dp, tt, rr, matrix_aux2)
665 0 : CALL mat_mulm(matrix_w1w1, matrix_avva, matrix_avva, n, -1.0_dp, 1.0_dp, tt, rr, matrix_aux2)
666 0 : CALL parallel_gemm("N", "N", n, n, n, 1.0_dp, matrix_avva, matrix_apVVpa, 1.0_dp, matrix_w1w1)
667 :
668 : ! Calculate w1e1w1: (warning: ev3 is scratch array)
669 :
670 0 : CALL mat_muld(matrix_aux, matrix_apVVpa, matrix_pe1p, n, 1.0_dp, 0.0_dp, tt, rr, matrix_aux2)
671 0 : CALL parallel_gemm("N", "N", n, n, n, 1.0_dp, matrix_avva, matrix_pe1p, 0.0_dp, matrix_aux2)
672 0 : CALL parallel_gemm("N", "N", n, n, n, 1.0_dp, matrix_aux, matrix_avva, 0.0_dp, matrix_w1e1w1)
673 0 : CALL mat_muld(matrix_w1e1w1, matrix_aux, matrix_apVVpa, n, -1.0_dp, 1.0_dp, tt, rr, matrix_ev3)
674 0 : CALL parallel_gemm("N", "N", n, n, n, -1.0_dp, matrix_aux2, matrix_avva, 1.0_dp, matrix_w1e1w1)
675 0 : CALL mat_muld(matrix_w1e1w1, matrix_aux2, matrix_apVVpa, n, 1.0_dp, 1.0_dp, tt, rr, matrix_ev3)
676 :
677 : !-----------------------------------------------------------------------
678 : ! 2. ev3 = 1/2 (W1^2)E1 + 1/2 E1(W1^2) - W1E1W1
679 : !-----------------------------------------------------------------------
680 :
681 0 : CALL parallel_gemm("N", "N", n, n, n, 0.5_dp, matrix_w1w1, matrix_ev1, 0.0_dp, matrix_ev3)
682 0 : CALL parallel_gemm("N", "N", n, n, n, 0.5_dp, matrix_ev1, matrix_w1w1, 1.0_dp, matrix_ev3)
683 0 : CALL cp_fm_scale_and_add(1.0_dp, matrix_ev3, -1.0_dp, matrix_w1e1w1)
684 :
685 : !-----------------------------------------------------------------------
686 : ! 3. Finish up the stuff!!
687 : !-----------------------------------------------------------------------
688 :
689 0 : CALL cp_fm_release(matrix_aux2)
690 0 : CALL cp_fm_release(matrix_avva)
691 0 : CALL cp_fm_release(matrix_apVVpa)
692 0 : CALL cp_fm_release(matrix_w1w1)
693 0 : CALL cp_fm_release(matrix_w1e1w1)
694 :
695 0 : CALL cp_fm_struct_release(matrix_full)
696 :
697 : ! WRITE (*,*) "CAW: DKH3 with even3b (Alex)"
698 : ! WRITE (*,*) "JT: Now available in cp2k"
699 :
700 0 : END SUBROUTINE even3b
701 :
702 : !-----------------------------------------------------------------------
703 :
704 : ! **************************************************************************************************
705 : !> \brief ...
706 : !> \param n ...
707 : !> \param ev4 ...
708 : !> \param e1 ...
709 : !> \param pe1p ...
710 : !> \param vv ...
711 : !> \param gg ...
712 : ! **************************************************************************************************
713 0 : SUBROUTINE even4a(n, ev4, e1, pe1p, vv, gg)
714 :
715 : !***********************************************************************
716 : ! *
717 : ! Alexander Wolf, last modified: 25.02.2002 -- DKH4 *
718 : ! *
719 : ! 4th order DK-approximation (scalar = spin-free) *
720 : ! *
721 : ! Version: 1.2 (25.2.2002) : Elegant (short) way of calculation *
722 : ! included *
723 : ! 1.1 (20.2.2002) : Usage of SR mat_add included *
724 : ! 1.0 (8.2.2002) *
725 : ! *
726 : ! ev4 = 1/2 [W2,[W1,E1]] + 1/8 [W1,[W1,[W1,O1]]] = *
727 : ! *
728 : ! = sum_1 + sum_2 *
729 : ! *
730 : ! *
731 : ! Modification history: *
732 : ! 30.09.2006 Jens Thar: deleted obsolete F77 memory manager *
733 : ! (not working yet) *
734 : ! *
735 : ! ---- Meaning of Parameters ---- *
736 : ! *
737 : ! n in Dimension of matrices *
738 : ! ev4 out even4 output matrix = final result *
739 : ! e1 in E1 *
740 : ! pe1p in p(E1)p *
741 : ! vv in potential v *
742 : ! gg in pvp *
743 : ! aa in A-Factors (DIAGONAL) *
744 : ! rr in R-Factors (DIAGONAL) *
745 : ! tt in Nonrel. kinetic Energy (DIAGONAL) *
746 : ! e in Rel. Energy = SQRT(p^2*c^2 + c^4) (DIAGONAL) *
747 : ! v symmetric (n x n)-matrix containing (A V A) *
748 : ! pvp symmetric (n x n)-matrix containing (A P V P A) *
749 : ! vh symmetric (n x n)-matrix containing (A V~ A) *
750 : ! pvph symmetric (n x n)-matrix containing (A P V~ P A) *
751 : ! w1w1 (W1)^2 *
752 : ! w1o1 W1*O1 (2-component formulation) *
753 : ! o1w1 O1*W1 (2-component formulation) *
754 : ! e1 symmetric (n x n)-matrix containing E1 *
755 : ! pe1p symmetric (n x n)-matrix containing p(E1)p *
756 : ! sum_i 2 addends defined above (i=1,2) *
757 : ! scr_i temporary (n x n)-scratch-matrices (i=1,..,4) *
758 : ! scrh_i temp. (n x n)-scr.-mat. with energy-denom. (i=1,..,4) *
759 : ! *
760 : !***********************************************************************
761 :
762 : INTEGER, INTENT(IN) :: n
763 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: ev4
764 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: e1, pe1p, vv, gg
765 :
766 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: o1w1, pvp, pvph, scr_1, scr_2, scr_3, &
767 0 : scr_4, scrh_1, scrh_2, scrh_3, scrh_4, &
768 0 : sum_1, sum_2, v, vh, w1o1, w1w1
769 :
770 : !C-----------------------------------------------------------------------
771 : !C 1. General Structures and Patterns for DKH4
772 : !C-----------------------------------------------------------------------
773 :
774 0 : ALLOCATE (v(n, n))
775 0 : ALLOCATE (pVp(n, n))
776 0 : ALLOCATE (vh(n, n))
777 0 : ALLOCATE (pVph(n, n))
778 0 : v = 0.0_dp
779 0 : pVp = 0.0_dp
780 0 : vh = 0.0_dp
781 0 : pVph = 0.0_dp
782 0 : v(1:n, 1:n) = vv(1:n, 1:n)
783 0 : vh(1:n, 1:n) = vv(1:n, 1:n)
784 0 : pvp(1:n, 1:n) = gg(1:n, 1:n)
785 0 : pvph(1:n, 1:n) = gg(1:n, 1:n)
786 :
787 0 : ev4 = 0.0_dp
788 : ! Calculate v = A V A:
789 :
790 : ! CALL mat_axa(v,n,aa)
791 :
792 : ! Calculate pvp = A P V P A:
793 :
794 : ! CALL mat_arxra(pvp,n,aa,rr)
795 :
796 : ! Calculate vh = A V~ A:
797 :
798 : ! CALL mat_1_over_h(vh,n,e)
799 : ! CALL mat_axa(vh,n,aa)
800 :
801 : ! Calculate pvph = A P V~ P A:
802 :
803 : ! CALL mat_1_over_h(pvph,n,e)
804 : ! CALL mat_arxra(pvph,n,aa,rr)
805 :
806 : ! Create/Initialize necessary matrices:
807 0 : ALLOCATE (w1w1(n, n))
808 0 : w1w1 = 0.0_dp
809 0 : ALLOCATE (w1o1(n, n))
810 0 : w1o1 = 0.0_dp
811 0 : ALLOCATE (o1w1(n, n))
812 0 : o1w1 = 0.0_dp
813 0 : ALLOCATE (sum_1(n, n))
814 0 : sum_1 = 0.0_dp
815 0 : ALLOCATE (sum_2(n, n))
816 0 : sum_2 = 0.0_dp
817 0 : ALLOCATE (scr_1(n, n))
818 0 : scr_1 = 0.0_dp
819 0 : ALLOCATE (scr_2(n, n))
820 0 : scr_2 = 0.0_dp
821 0 : ALLOCATE (scr_3(n, n))
822 0 : scr_3 = 0.0_dp
823 0 : ALLOCATE (scr_4(n, n))
824 0 : scr_4 = 0.0_dp
825 0 : ALLOCATE (scrh_1(n, n))
826 0 : scrh_1 = 0.0_dp
827 0 : ALLOCATE (scrh_2(n, n))
828 0 : scrh_2 = 0.0_dp
829 0 : ALLOCATE (scrh_3(n, n))
830 0 : scrh_3 = 0.0_dp
831 0 : ALLOCATE (scrh_4(n, n))
832 0 : scrh_4 = 0.0_dp
833 :
834 : ! Calculate w1w1:
835 0 : CALL dgemm("N", "N", n, n, n, 1.0_dp, pvph, n, vh, n, 0.0_dp, w1w1, n)
836 : ! CALL mat_muld(w1w1,pvph,pvph,n, -1.0_dp,1.0_dp,tt,rr)
837 : ! CALL mat_mulm(w1w1,vh, vh,n, -1.0_dp,1.0_dp,tt,rr)
838 0 : CALL dgemm("N", "N", n, n, n, 1.0_dp, vh, n, pvph, n, 1.0_dp, w1w1, n)
839 :
840 : ! Calculate w1o1:
841 0 : CALL dgemm("N", "N", n, n, n, -1.0_dp, pvph, n, v, n, 0.0_dp, w1o1, n)
842 : ! CALL mat_muld(w1o1,pvph,pvp,n, 1.0_dp,1.0_dp,tt,rr)
843 : ! CALL mat_mulm(w1o1,vh, v,n, 1.0_dp,1.0_dp,tt,rr)
844 0 : CALL dgemm("N", "N", n, n, n, -1.0_dp, vh, n, pvp, n, 1.0_dp, w1o1, n)
845 : ! Calculate o1w1:
846 0 : CALL dgemm("N", "N", n, n, n, 1.0_dp, pvp, n, vh, n, 0.0_dp, o1w1, n)
847 : ! CALL mat_muld(o1w1,pvp,pvph,n, -1.0_dp,1.0_dp,tt,rr)
848 : ! CALL mat_mulm(o1w1,v, vh,n, -1.0_dp,1.0_dp,tt,rr)
849 0 : CALL dgemm("N", "N", n, n, n, 1.0_dp, v, n, pvph, n, 1.0_dp, o1w1, n)
850 :
851 : !-----------------------------------------------------------------------
852 : ! 2. sum_1 = 1/2 [W2,[W1,E1]] = 1/2 (W2W1E1 - W2E1W1 - W1E1W2 + E1W1W2)
853 : !-----------------------------------------------------------------------
854 :
855 : ! scr_i & scrh_i for steps 2a (W2W1E1) and 2b (W2E1W1):
856 :
857 0 : CALL dgemm("N", "N", n, n, n, 1.0_dp, vh, n, e1, n, 0.0_dp, scr_1, n)
858 0 : CALL dgemm("N", "N", n, n, n, 1.0_dp, pvph, n, e1, n, 0.0_dp, scr_2, n)
859 0 : CALL dgemm("N", "N", n, n, n, 1.0_dp, pe1p, n, vh, n, 0.0_dp, scr_3, n)
860 : ! CALL mat_muld(scr_4, pe1p,pvph,n,1.0_dp,0.0_dp,tt,rr)
861 :
862 : ! CALL mat_muld(scrh_1,pvph,pe1p,n,1.0_dp,0.0_dp,tt,rr)
863 : ! CALL mat_1_over_h(scrh_1,n,e)
864 0 : CALL dgemm("N", "N", n, n, n, 1.0_dp, vh, n, pe1p, n, 0.0_dp, scrh_2, n)
865 : ! CALL mat_1_over_h(scrh_2,n,e)
866 0 : CALL dgemm("N", "N", n, n, n, 1.0_dp, e1, n, pvph, n, 0.0_dp, scrh_3, n)
867 : ! CALL mat_1_over_h(scrh_3,n,e)
868 0 : CALL dgemm("N", "N", n, n, n, 1.0_dp, e1, n, vh, n, 0.0_dp, scrh_4, n)
869 : ! CALL mat_1_over_h(scrh_4,n,e)
870 :
871 : ! 2a) sum_1 = 1/2 W2W1E1 ( [1]-[8] )
872 :
873 0 : CALL dgemm("N", "N", n, n, n, 0.5_dp, scrh_1, n, scr_1, n, 0.0_dp, sum_1, n)
874 : ! CALL mat_muld(sum_1,scrh_1,scr_2,n,-0.5_dp,1.0_dp,tt,rr)
875 0 : CALL dgemm("N", "N", n, n, n, -0.5_dp, scrh_2, n, scr_1, n, 1.0_dp, sum_1, n)
876 : ! CALL mat_muld(sum_1,scrh_2,scr_2,n, 0.5_dp,1.0_dp,tt,rr)
877 0 : CALL dgemm("N", "N", n, n, n, -0.5_dp, scrh_3, n, scr_1, n, 1.0_dp, sum_1, n)
878 : ! CALL mat_muld(sum_1,scrh_3,scr_2,n, 0.5_dp,1.0_dp,tt,rr)
879 : ! CALL mat_mulm(sum_1,scrh_4,scr_1,n, 0.5_dp,1.0_dp,tt,rr)
880 0 : CALL dgemm("N", "N", n, n, n, -0.5_dp, scrh_4, n, scr_2, n, 1.0_dp, sum_1, n)
881 :
882 : ! 2b) sum_1 = - 1/2 W2E1W1 (+ sum_1) ( [9]-[16] )
883 :
884 : ! CALL mat_muld(sum_1,scrh_1,scr_3,n,-0.5_dp,1.0_dp,tt,rr)
885 : ! CALL mat_muld(sum_1,scrh_1,scr_4,n, 0.5_dp,1.0_dp,tt,rr)
886 : ! CALL mat_muld(sum_1,scrh_2,scr_3,n, 0.5_dp,1.0_dp,tt,rr)
887 : ! CALL mat_muld(sum_1,scrh_2,scr_4,n,-0.5_dp,1.0_dp,tt,rr)
888 : ! CALL mat_muld(sum_1,scrh_3,scr_3,n, 0.5_dp,1.0_dp,tt,rr)
889 : ! CALL mat_muld(sum_1,scrh_3,scr_4,n,-0.5_dp,1.0_dp,tt,rr)
890 0 : CALL dgemm("N", "N", n, n, n, -0.5_dp, scrh_4, n, scr_3, n, 1.0_dp, sum_1, n)
891 0 : CALL dgemm("N", "N", n, n, n, 0.5_dp, scrh_4, n, scr_4, n, 1.0_dp, sum_1, n)
892 :
893 : ! scr_i & scrh_i for steps 2c (W1E1W2) and 2d (E1W1W2):
894 :
895 : ! CALL mat_muld(scr_1, pvph,pe1p,n,1.0_dp,0.0_dp,tt,rr)
896 0 : CALL dgemm("N", "N", n, n, n, 1.0_dp, vh, n, pe1p, n, 0.0_dp, scr_2, n)
897 0 : CALL dgemm("N", "N", n, n, n, 1.0_dp, e1, n, pvph, n, 0.0_dp, scr_3, n)
898 0 : CALL dgemm("N", "N", n, n, n, 1.0_dp, e1, n, vh, n, 0.0_dp, scr_4, n)
899 :
900 0 : CALL dgemm("N", "N", n, n, n, 1.0_dp, vh, n, e1, n, 0.0_dp, scrh_1, n)
901 : ! CALL mat_1_over_h(scrh_1,n,e)
902 0 : CALL dgemm("N", "N", n, n, n, 1.0_dp, pvph, n, e1, n, 0.0_dp, scrh_2, n)
903 : ! CALL mat_1_over_h(scrh_2,n,e)
904 0 : CALL dgemm("N", "N", n, n, n, 1.0_dp, pe1p, n, vh, n, 0.0_dp, scr_3, n)
905 : ! CALL mat_1_over_h(scrh_3,n,e)
906 : ! CALL mat_muld(scrh_4,pe1p,pvph,n,1.0_dp,0.0_dp,tt,rr)
907 : ! CALL mat_1_over_h(scrh_4,n,e)
908 :
909 : ! 2c) sum_1 = - 1/2 W1E1W2 (+ sum_1) ( [17]-[24] )
910 :
911 0 : CALL dgemm("N", "N", n, n, n, 0.5_dp, scr_1, n, scrh_1, n, 0.0_dp, sum_1, n)
912 : ! CALL mat_muld(sum_1,scr_1,scrh_2,n,-0.5_dp,1.0_dp,tt,rr)
913 0 : CALL dgemm("N", "N", n, n, n, -0.5_dp, scr_2, n, scrh_1, n, 1.0_dp, sum_1, n)
914 : ! CALL mat_muld(sum_1,scr_2,scrh_2,n, 0.5_dp,1.0_dp,tt,rr)
915 : ! CALL mat_muld(sum_1,scr_1,scrh_3,n,-0.5_dp,1.0_dp,tt,rr)
916 : ! CALL mat_muld(sum_1,scr_1,scrh_4,n, 0.5_dp,1.0_dp,tt,rr)
917 : ! CALL mat_muld(sum_1,scr_2,scrh_3,n, 0.5_dp,1.0_dp,tt,rr)
918 : ! CALL mat_muld(sum_1,scr_2,scrh_4,n,-0.5_dp,1.0_dp,tt,rr)
919 :
920 : ! 2d) sum_1 = 1/2 E1W1W2 (+ sum_1) ( [25]-[32] )
921 :
922 0 : CALL dgemm("N", "N", n, n, n, -0.5_dp, scr_3, n, scrh_1, n, 0.0_dp, sum_1, n)
923 : ! CALL mat_muld(sum_1,scr_3,scrh_2,n, 0.5_dp,1.0_dp,tt,rr)
924 : ! CALL mat_mulm(sum_1,scr_4,scrh_1,n, 0.5_dp,1.0_dp,tt,rr)
925 0 : CALL dgemm("N", "N", n, n, n, -0.5_dp, scr_4, n, scrh_2, n, 1.0_dp, sum_1, n)
926 : ! CALL mat_muld(sum_1,scr_3,scrh_3,n, 0.5_dp,1.0_dp,tt,rr)
927 : ! CALL mat_muld(sum_1,scr_3,scrh_4,n,-0.5_dp,1.0_dp,tt,rr)
928 0 : CALL dgemm("N", "N", n, n, n, -0.5_dp, scr_4, n, scrh_3, n, 1.0_dp, sum_1, n)
929 0 : CALL dgemm("N", "N", n, n, n, 0.5_dp, scr_4, n, scrh_4, n, 1.0_dp, sum_1, n)
930 :
931 : !-----------------------------------------------------------------------
932 : ! 3. sum_2 = 1/8 [W1,[W1,[W1,O1]]] =
933 : !
934 : ! = 1/8 ( (W1^3)O1 - 3(W1^2)O1W1 + 3 W1O1(W1^2) - O1(W1^3) )
935 : !-----------------------------------------------------------------------
936 :
937 0 : CALL dgemm("N", "N", n, n, n, 0.125_dp, w1w1, n, w1o1, n, 0.0_dp, sum_2, n)
938 0 : CALL dgemm("N", "N", n, n, n, -0.375_dp, w1w1, n, o1w1, n, 1.0_dp, sum_2, n)
939 0 : CALL dgemm("N", "N", n, n, n, 0.375_dp, w1o1, n, w1w1, n, 1.0_dp, sum_2, n)
940 0 : CALL dgemm("N", "N", n, n, n, -0.125_dp, o1w1, n, w1w1, n, 1.0_dp, sum_2, n)
941 :
942 : !-----------------------------------------------------------------------
943 : ! 4. result = sum_1 + sum_2
944 : !-----------------------------------------------------------------------
945 :
946 0 : CALL mat_add(ev4, 1.0_dp, sum_1, 1.0_dp, sum_2, n)
947 :
948 : !-----------------------------------------------------------------------
949 : ! 5. Finish up the stuff!!
950 : !-----------------------------------------------------------------------
951 :
952 0 : DEALLOCATE (v, pvp, vh, pvph, w1w1, w1o1, o1w1, sum_1, sum_2)
953 0 : DEALLOCATE (scr_1, scr_2, scr_3, scr_4, scrh_1, scrh_2, scrh_3, scrh_4)
954 :
955 : ! WRITE (*,*) "CAW: DKH4 with even4a (Alex)"
956 : ! WRITE (*,*) "JT: Now available in cp2k"
957 :
958 0 : END SUBROUTINE even4a
959 :
960 : !-----------------------------------------------------------------------
961 : ! -
962 : ! Matrix routines for DKH-procedure -
963 : ! Alexander Wolf -
964 : ! modified: Jens Thar: Mem manager deleted -
965 : ! This file contains the -
966 : ! following subroutines: -
967 : ! 1. mat_1_over_h -
968 : ! 2. mat_axa -
969 : ! 3. mat_arxra -
970 : ! 4. mat_mulm -
971 : ! 5. mat_muld -
972 : ! 6. mat_add -
973 : ! -
974 : !-----------------------------------------------------------------------
975 :
976 : ! **************************************************************************************************
977 : !> \brief ...
978 : !> \param matrix_p ...
979 : !> \param matrix_pp ...
980 : !> \param e ...
981 : !> \param matrix_aux ...
982 : ! **************************************************************************************************
983 0 : SUBROUTINE mat_1_over_h(matrix_p, matrix_pp, e, matrix_aux)
984 :
985 : !***********************************************************************
986 : ! *
987 : ! 2. SR mat_1_over_h: Transform matrix p into matrix p/(e(i)+e(j)) *
988 : ! *
989 : ! p in REAL(:,:) : matrix p *
990 : ! e in REAL(:) : rel. energy (diagonal) *
991 : ! *
992 : !***********************************************************************
993 :
994 : TYPE(cp_fm_type), INTENT(IN) :: matrix_p, matrix_pp
995 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: e
996 : TYPE(cp_fm_type), INTENT(IN) :: matrix_aux
997 :
998 : INTEGER :: i, j, ncol_local, nrow_local
999 0 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1000 :
1001 : CALL cp_fm_get_info(matrix_aux, nrow_local=nrow_local, ncol_local=ncol_local, &
1002 0 : row_indices=row_indices, col_indices=col_indices)
1003 :
1004 0 : DO j = 1, ncol_local
1005 0 : DO i = 1, nrow_local
1006 0 : matrix_aux%local_data(i, j) = 1/(e(row_indices(i)) + e(col_indices(j)))
1007 : END DO
1008 : END DO
1009 :
1010 0 : CALL cp_fm_schur_product(matrix_p, matrix_aux, matrix_pp)
1011 :
1012 0 : END SUBROUTINE mat_1_over_h
1013 :
1014 : ! **************************************************************************************************
1015 : !> \brief ...
1016 : !> \param matrix_x ...
1017 : !> \param matrix_axa ...
1018 : !> \param n ...
1019 : !> \param a ...
1020 : !> \param matrix_aux ...
1021 : ! **************************************************************************************************
1022 0 : SUBROUTINE mat_axa(matrix_x, matrix_axa, n, a, matrix_aux)
1023 :
1024 : !C***********************************************************************
1025 : !C *
1026 : !C 3. SR mat_axa: Transform matrix p into matrix a*p*a *
1027 : !C *
1028 : !C p in REAL(:,:): matrix p *
1029 : !C a in REAL(:) : A-factors (diagonal) *
1030 : !CJT n in INTEGER : dimension of matrix p *
1031 : !C *
1032 : !C***********************************************************************
1033 :
1034 : TYPE(cp_fm_type), INTENT(IN) :: matrix_x, matrix_axa
1035 : INTEGER, INTENT(IN) :: n
1036 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: a
1037 : TYPE(cp_fm_type), INTENT(IN) :: matrix_aux
1038 :
1039 : INTEGER :: i, nrow_local
1040 0 : INTEGER, DIMENSION(:), POINTER :: row_indices
1041 : TYPE(cp_blacs_env_type), POINTER :: context
1042 : TYPE(cp_fm_struct_type), POINTER :: vec_full
1043 : TYPE(cp_fm_type) :: vec_a
1044 :
1045 0 : CALL cp_fm_struct_get(matrix_x%matrix_struct, context=context)
1046 : CALL cp_fm_struct_create(fmstruct=vec_full, &
1047 : context=context, &
1048 : nrow_global=n, &
1049 0 : ncol_global=1)
1050 :
1051 0 : CALL cp_fm_create(vec_a, vec_full)
1052 :
1053 : CALL cp_fm_get_info(matrix_x, nrow_local=nrow_local, &
1054 0 : row_indices=row_indices)
1055 :
1056 0 : DO i = 1, nrow_local
1057 0 : vec_a%local_data(i, 1) = a(row_indices(i))
1058 : END DO
1059 :
1060 0 : CALL cp_fm_syrk('U', 'N', 1, 1.0_dp, vec_a, 1, 1, 0.0_dp, matrix_aux)
1061 0 : CALL cp_fm_uplo_to_full(matrix_aux, matrix_axa)
1062 0 : CALL cp_fm_schur_product(matrix_x, matrix_aux, matrix_axa)
1063 :
1064 : ! DO i=1,n
1065 : ! DO j=1,n
1066 : ! p(i,j)=p(i,j)*a(i)*a(j)
1067 : ! END DO
1068 : ! END DO
1069 :
1070 0 : CALL cp_fm_release(vec_a)
1071 0 : CALL cp_fm_struct_release(vec_full)
1072 :
1073 0 : END SUBROUTINE mat_axa
1074 :
1075 : ! **************************************************************************************************
1076 : !> \brief ...
1077 : !> \param matrix_x ...
1078 : !> \param matrix_axa ...
1079 : !> \param n ...
1080 : !> \param a ...
1081 : !> \param r ...
1082 : !> \param matrix_aux ...
1083 : ! **************************************************************************************************
1084 0 : SUBROUTINE mat_arxra(matrix_x, matrix_axa, n, a, r, matrix_aux)
1085 :
1086 : !C***********************************************************************
1087 : !C *
1088 : !C 4. SR mat_arxra: Transform matrix p into matrix a*r*p*r*a *
1089 : !C *
1090 : !C p in REAL(:,:) : matrix p *
1091 : !C a in REAL(:) : A-factors (diagonal) *
1092 : !C r in REAL(:) : R-factors (diagonal) *
1093 : !C n in INTEGER : dimension of matrix p *
1094 : !C *
1095 : !C***********************************************************************
1096 :
1097 : TYPE(cp_fm_type), INTENT(IN) :: matrix_x, matrix_axa
1098 : INTEGER, INTENT(IN) :: n
1099 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: a, r
1100 : TYPE(cp_fm_type), INTENT(IN) :: matrix_aux
1101 :
1102 : INTEGER :: i, nrow_local
1103 0 : INTEGER, DIMENSION(:), POINTER :: row_indices
1104 : TYPE(cp_blacs_env_type), POINTER :: context
1105 : TYPE(cp_fm_struct_type), POINTER :: vec_full
1106 : TYPE(cp_fm_type) :: vec_a
1107 :
1108 0 : CALL cp_fm_struct_get(matrix_x%matrix_struct, context=context)
1109 : CALL cp_fm_struct_create(fmstruct=vec_full, &
1110 : context=context, &
1111 : nrow_global=n, &
1112 0 : ncol_global=1)
1113 :
1114 : CALL cp_fm_get_info(matrix_aux, nrow_local=nrow_local, &
1115 0 : row_indices=row_indices)
1116 :
1117 0 : CALL cp_fm_create(vec_a, vec_full)
1118 :
1119 0 : DO i = 1, nrow_local
1120 0 : vec_a%local_data(i, 1) = a(row_indices(i))*r(row_indices(i))
1121 : END DO
1122 :
1123 0 : CALL cp_fm_syrk('U', 'N', 1, 1.0_dp, vec_a, 1, 1, 0.0_dp, matrix_aux)
1124 0 : CALL cp_fm_uplo_to_full(matrix_aux, matrix_axa)
1125 0 : CALL cp_fm_schur_product(matrix_x, matrix_aux, matrix_axa)
1126 :
1127 0 : CALL cp_fm_release(vec_a)
1128 0 : CALL cp_fm_struct_release(vec_full)
1129 :
1130 0 : END SUBROUTINE mat_arxra
1131 :
1132 : ! **************************************************************************************************
1133 : !> \brief ...
1134 : !> \param matrix_p ...
1135 : !> \param matrix_q ...
1136 : !> \param matrix_r ...
1137 : !> \param n ...
1138 : !> \param alpha ...
1139 : !> \param beta ...
1140 : !> \param t ...
1141 : !> \param rr ...
1142 : !> \param matrix_aux ...
1143 : ! **************************************************************************************************
1144 0 : SUBROUTINE mat_mulm(matrix_p, matrix_q, matrix_r, n, alpha, beta, t, rr, matrix_aux)
1145 :
1146 : !C***********************************************************************
1147 : !C *
1148 : !C 5. SR mat_mulm: Multiply matrices according to: *
1149 : !C *
1150 : !C p = alpha*q*(..P^2..)*r + beta*p *
1151 : !C *
1152 : !C p out REAL(:,:): matrix p *
1153 : !C q in REAL(:,:): matrix q *
1154 : !C r in REAL(:,.): matrix r *
1155 : !C n in INTEGER : dimension n of matrices *
1156 : !C alpha in REAL(dp) : *
1157 : !C beta in REAL(dp) : *
1158 : !C t in REAL(:) : non-rel. kinetic energy (diagonal) *
1159 : !C rr in REAL(:) : R-factors (diagonal) *
1160 : !C *
1161 : !C***********************************************************************
1162 :
1163 : TYPE(cp_fm_type), INTENT(IN) :: matrix_p, matrix_q, matrix_r
1164 : INTEGER, INTENT(IN) :: n
1165 : REAL(KIND=dp), INTENT(IN) :: alpha, beta
1166 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: t, rr
1167 : TYPE(cp_fm_type), INTENT(IN) :: matrix_aux
1168 :
1169 : INTEGER :: i
1170 0 : REAL(KIND=dp), DIMENSION(n) :: vec
1171 :
1172 0 : CALL cp_fm_to_fm(matrix_q, matrix_aux)
1173 :
1174 0 : DO i = 1, n
1175 0 : vec(i) = 2.0_dp*t(i)*rr(i)*rr(i)
1176 : END DO
1177 0 : CALL cp_fm_column_scale(matrix_aux, vec)
1178 :
1179 0 : CALL parallel_gemm("N", "N", n, n, n, alpha, matrix_aux, matrix_r, beta, matrix_p)
1180 :
1181 0 : END SUBROUTINE mat_mulm
1182 :
1183 : ! **************************************************************************************************
1184 : !> \brief ...
1185 : !> \param matrix_p ...
1186 : !> \param matrix_q ...
1187 : !> \param matrix_r ...
1188 : !> \param n ...
1189 : !> \param alpha ...
1190 : !> \param beta ...
1191 : !> \param t ...
1192 : !> \param rr ...
1193 : !> \param matrix_aux ...
1194 : ! **************************************************************************************************
1195 0 : SUBROUTINE mat_muld(matrix_p, matrix_q, matrix_r, n, alpha, beta, t, rr, matrix_aux)
1196 :
1197 : !C***********************************************************************
1198 : !C *
1199 : !C 16. SR mat_muld: Multiply matrices according to: *
1200 : !C *
1201 : !C p = alpha*q*(..1/P^2..)*r + beta*p *
1202 : !C *
1203 : !C p out REAL(:,:): matrix p *
1204 : !C q in REAL(:,:): matrix q *
1205 : !C r in REAL(:,:): matrix r *
1206 : !C n in INTEGER : Dimension of all matrices *
1207 : !C alpha in REAL(dp) : *
1208 : !C beta in REAL(dp) : *
1209 : !C t in REAL(:) : non-rel. kinetic energy (diagonal) *
1210 : !C rr in REAL(:) : R-factors (diagonal) *
1211 : !C *
1212 : !C***********************************************************************
1213 :
1214 : TYPE(cp_fm_type), INTENT(IN) :: matrix_p, matrix_q, matrix_r
1215 : INTEGER, INTENT(IN) :: n
1216 : REAL(KIND=dp), INTENT(IN) :: alpha, beta
1217 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: t, rr
1218 : TYPE(cp_fm_type), INTENT(IN) :: matrix_aux
1219 :
1220 : INTEGER :: i
1221 0 : REAL(KIND=dp), DIMENSION(n) :: vec
1222 :
1223 0 : CALL cp_fm_to_fm(matrix_q, matrix_aux)
1224 :
1225 0 : DO i = 1, n
1226 0 : vec(i) = 0.5_dp/(t(i)*rr(i)*rr(i))
1227 : END DO
1228 :
1229 0 : CALL cp_fm_column_scale(matrix_aux, vec)
1230 :
1231 0 : CALL parallel_gemm("N", "N", n, n, n, alpha, matrix_aux, matrix_r, beta, matrix_p)
1232 :
1233 0 : END SUBROUTINE mat_muld
1234 :
1235 : ! **************************************************************************************************
1236 : !> \brief ...
1237 : !> \param s ...
1238 : !> \param v ...
1239 : !> \param h ...
1240 : !> \param pVp ...
1241 : !> \param n ...
1242 : !> \param dkh_order ...
1243 : ! **************************************************************************************************
1244 272 : SUBROUTINE DKH_atom_transformation(s, v, h, pVp, n, dkh_order)
1245 :
1246 : !-----------------------------------------------------------------------
1247 : ! *
1248 : ! INPUT: *
1249 : ! n Number of primitive gaussians *
1250 : ! s (:,:) overlap matrix *
1251 : ! pVp (:,:) pVp matrix *
1252 : ! *
1253 : ! IN_OUT: *
1254 : ! v (:,:) input: nonrelativistic potential energy matrix *
1255 : ! output: (ev1+ev2) *
1256 : ! h (:,:) input: kinetic energy matrix *
1257 : ! output: kinetic part of hamiltonian in position space *
1258 : ! *
1259 : ! INTERNAL *
1260 : ! sinv (:,:) inverted, orthogonalized overlap matrix *
1261 : ! ev0t (:) DKH-even0 matrix in T-basis *
1262 : ! e (:) e=SQRT(p^2c^2+c^4) *
1263 : ! eig (:,:) eigenvectors of sinv' h sinv *
1264 : ! tt (:) eigenvalues of sinv' h sinv *
1265 : ! revt (:,:) reverse transformation matrix T-basis -> position space*
1266 : ! aa (:) kinematical factors f. DKH SQRT((c^2+e(i))/(2.0*e(i))) *
1267 : ! rr (:) kinematical factors f. DKH c/(c^2+e(i)) *
1268 : ! vt (:,:) non relativistic potential matrix in T-basis *
1269 : ! pvpt (:,:) pvp integral matrix in T-basis *
1270 : ! ev1t (:,:) DKH-even1 matrix in T-basis *
1271 : ! evt2 (:,:) DKH-even2 matrix in T-basis *
1272 : ! ev1 (:,:) DKH-even1 matrix in position space *
1273 : ! ev2 (:,:) DKH-even2 matrix in position space *
1274 : ! ove (:,:) scratch *
1275 : ! aux (:,:) scratch *
1276 : ! c_light_au velocity of light 137 a.u. *
1277 : ! prea prefactor, 1/137^2 *
1278 : ! con2 prefactor, 2/137^2 *
1279 : ! con prefactor, 137^2 *
1280 : !-----------------------------------------------------------------------
1281 :
1282 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: s, v, h, pVp
1283 : INTEGER, INTENT(IN) :: n, dkh_order
1284 :
1285 : INTEGER :: i, j, k
1286 272 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: aa, e, ev0t, rr, tt
1287 272 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: aux, eig, ev1, ev1t, ev2, ev2t, ev3, &
1288 272 : ev3t, ev4, ev4t, ove, pev1tp, pVpt, &
1289 272 : revt, sinv, vt
1290 :
1291 272 : IF (dkh_order < 0) RETURN
1292 :
1293 : !CAW pp: p^2-values (in momentum-space), stored as matrix!!
1294 :
1295 : !-----------------------------------------------------------------------
1296 : ! Allocate some matrices
1297 : !-----------------------------------------------------------------------
1298 :
1299 1088 : ALLOCATE (eig(n, n))
1300 816 : ALLOCATE (sinv(n, n))
1301 816 : ALLOCATE (revt(n, n))
1302 816 : ALLOCATE (aux(n, n))
1303 816 : ALLOCATE (ove(n, n))
1304 816 : ALLOCATE (ev0t(n))
1305 544 : ALLOCATE (e(n))
1306 544 : ALLOCATE (aa(n))
1307 544 : ALLOCATE (rr(n))
1308 544 : ALLOCATE (tt(n))
1309 816 : ALLOCATE (ev1t(n, n))
1310 816 : ALLOCATE (ev2t(n, n))
1311 816 : ALLOCATE (ev3t(n, n))
1312 816 : ALLOCATE (ev4t(n, n))
1313 816 : ALLOCATE (vt(n, n))
1314 816 : ALLOCATE (pVpt(n, n))
1315 816 : ALLOCATE (pev1tp(n, n))
1316 816 : ALLOCATE (ev1(n, n))
1317 816 : ALLOCATE (ev2(n, n))
1318 816 : ALLOCATE (ev3(n, n))
1319 816 : ALLOCATE (ev4(n, n))
1320 :
1321 : !-----------------------------------------------------------------------
1322 : ! Schmidt-orthogonalize overlap matrix
1323 : !-----------------------------------------------------------------------
1324 :
1325 272 : CALL sog(n, s, sinv)
1326 :
1327 : !-----------------------------------------------------------------------
1328 : ! Calculate matrix representation from nonrelativistic T matrix
1329 : !-----------------------------------------------------------------------
1330 :
1331 272 : CALL dkh_diag(h, n, eig, tt, sinv, aux, 0)
1332 :
1333 : !-----------------------------------------------------------------------
1334 : ! Calculate kinetic part of Hamiltonian in T-basis
1335 : !-----------------------------------------------------------------------
1336 :
1337 272 : CALL kintegral_a(n, ev0t, tt, e)
1338 :
1339 : !-----------------------------------------------------------------------
1340 : ! Calculate reverse transformation matrix revt
1341 : !-----------------------------------------------------------------------
1342 :
1343 272 : CALL dgemm("N", "N", n, n, n, 1.0_dp, sinv, n, eig, n, 0.0_dp, aux, n)
1344 70072 : CALL dgemm("N", "N", n, n, n, 1.0_dp, s, n, aux, n, 0.0_dp, revt, n)
1345 :
1346 : !-----------------------------------------------------------------------
1347 : ! Calculate kinetic part of the Hamiltonian
1348 : !-----------------------------------------------------------------------
1349 :
1350 223656 : h = 0.0_dp
1351 6950 : DO i = 1, n
1352 118642 : DO j = 1, i
1353 4116802 : DO k = 1, n
1354 3998432 : h(i, j) = h(i, j) + revt(i, k)*revt(j, k)*ev0t(k)
1355 4110124 : h(j, i) = h(i, j)
1356 : END DO
1357 : END DO
1358 : END DO
1359 :
1360 : !-----------------------------------------------------------------------
1361 : ! Calculate kinematical factors for DKH
1362 : !-----------------------------------------------------------------------
1363 :
1364 6950 : DO i = 1, n
1365 6678 : aa(i) = SQRT((c_light_au**2 + e(i))/(2.0_dp*e(i)))
1366 6950 : rr(i) = SQRT(c_light_au**2)/(c_light_au**2 + e(i))
1367 : END DO
1368 :
1369 : !-----------------------------------------------------------------------
1370 : ! Transform v integrals to T-basis (v -> vt)
1371 : !-----------------------------------------------------------------------
1372 :
1373 272 : CALL trsm(v, sinv, ove, n, aux)
1374 272 : CALL trsm(ove, eig, vt, n, aux)
1375 :
1376 : !-----------------------------------------------------------------------
1377 : ! Transform pVp integrals to T-basis (pVp -> pVpt)
1378 : !-----------------------------------------------------------------------
1379 :
1380 272 : CALL trsm(pVp, sinv, ove, n, aux)
1381 272 : CALL trsm(ove, eig, pVpt, n, aux)
1382 :
1383 : !-----------------------------------------------------------------------
1384 : ! Calculate even1 in T-basis
1385 : !-----------------------------------------------------------------------
1386 :
1387 272 : IF (dkh_order >= 1) THEN
1388 266 : CALL even1_a(n, ev1t, vt, pvpt, aa, rr)
1389 :
1390 : !----------------------------------------------------------------------
1391 : ! Transform even1 back to position space
1392 : !----------------------------------------------------------------------
1393 :
1394 266 : CALL dgemm("N", "N", n, n, n, 1.0_dp, revt, n, ev1t, n, 0.0_dp, aux, n)
1395 266 : CALL dgemm("N", "T", n, n, n, 1.0_dp, aux, n, revt, n, 0.0_dp, ev1, n)
1396 : END IF
1397 :
1398 : !-----------------------------------------------------------------------
1399 : ! Calculate even2 in T-basis
1400 : !-----------------------------------------------------------------------
1401 :
1402 272 : IF (dkh_order >= 2) THEN
1403 258 : CALL even2c_a(n, ev2t, vt, pvpt, aa, rr, tt, e)
1404 :
1405 : !-----------------------------------------------------------------------
1406 : ! Transform even2 back to position space
1407 : !-----------------------------------------------------------------------
1408 :
1409 214066 : aux = 0.0_dp
1410 258 : CALL dgemm("N", "N", n, n, n, 1.0_dp, revt, n, ev2t, n, 0.0_dp, aux, n)
1411 258 : CALL dgemm("N", "T", n, n, n, 1.0_dp, aux, n, revt, n, 0.0_dp, ev2, n)
1412 : END IF
1413 :
1414 : !-----------------------------------------------------------------------
1415 : ! Calculate even3 in T-basis, only if requested
1416 : !-----------------------------------------------------------------------
1417 :
1418 272 : IF (dkh_order >= 3) THEN
1419 164 : CALL peven1p_a(n, pev1tp, vt, pvpt, aa, rr, tt)
1420 164 : CALL even3b_a(n, ev3t, ev1t, pev1tp, vt, pvpt, aa, rr, tt, e)
1421 :
1422 : !-----------------------------------------------------------------------
1423 : ! Transform even3 back to position space
1424 : !-----------------------------------------------------------------------
1425 144476 : aux = 0.0_dp
1426 164 : CALL dgemm("N", "N", n, n, n, 1.0_dp, revt, n, ev3t, n, 0.0_dp, aux, n)
1427 164 : CALL dgemm("N", "T", n, n, n, 1.0_dp, aux, n, revt, n, 0.0_dp, ev3, n)
1428 :
1429 : !-----------------------------------------------------------------------
1430 : ! Calculate even4 in T-basis, only if requested
1431 : !-----------------------------------------------------------------------
1432 :
1433 164 : IF (dkh_order >= 4) THEN
1434 0 : CALL even4a_a(n, ev4t, ev1t, pev1tp, vt, pvpt, aa, rr, tt, e)
1435 :
1436 : !-----------------------------------------------------------------------
1437 : ! Transform even4 back to position space
1438 : !-----------------------------------------------------------------------
1439 0 : aux = 0.0_dp
1440 0 : CALL dgemm("N", "N", n, n, n, 1.0_dp, revt, n, ev4t, n, 0.0_dp, aux, n)
1441 0 : CALL dgemm("N", "T", n, n, n, 1.0_dp, aux, n, revt, n, 0.0_dp, ev4, n)
1442 : END IF
1443 : END IF
1444 :
1445 272 : IF (dkh_order >= 4) THEN
1446 0 : CPABORT("DKH 4")
1447 : END IF
1448 : !-----------------------------------------------------------------------
1449 : ! Calculate v in position space
1450 : !-----------------------------------------------------------------------
1451 :
1452 272 : IF (dkh_order >= 1) THEN
1453 266 : CALL mat_add2(v, 0.0_dp, 1.0_dp, ev1, n)
1454 : END IF
1455 272 : IF (dkh_order >= 2) THEN
1456 258 : CALL mat_add2(v, 1.0_dp, 1.0_dp, ev2, n)
1457 : END IF
1458 272 : IF (dkh_order >= 3) THEN
1459 164 : CALL mat_add2(v, 1.0_dp, 1.0_dp, ev3, n)
1460 : END IF
1461 272 : IF (dkh_order >= 4) THEN
1462 0 : CALL mat_add2(v, 1.0_dp, 1.0_dp, ev4, n)
1463 : END IF
1464 :
1465 : !-----------------------------------------------------------------------
1466 :
1467 272 : DEALLOCATE (eig, sinv, revt, ove, aux, vt, pVpt, ev1, ev2, ev3, ev4, ev1t, ev2t, ev3t, ev4t, pev1tp)
1468 272 : DEALLOCATE (ev0t, e, aa, rr, tt)
1469 :
1470 272 : END SUBROUTINE dkh_atom_transformation
1471 :
1472 : ! **************************************************************************************************
1473 : !> \brief ...
1474 : !> \param n ...
1475 : !> \param ev0t ...
1476 : !> \param tt ...
1477 : !> \param e ...
1478 : ! **************************************************************************************************
1479 272 : SUBROUTINE kintegral_a(n, ev0t, tt, e)
1480 :
1481 : INTEGER, INTENT(IN) :: n
1482 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: ev0t
1483 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: tt
1484 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: e
1485 :
1486 : INTEGER :: i
1487 : REAL(KIND=dp) :: con, con2, prea, ratio, tv1, tv2, tv3, &
1488 : tv4
1489 :
1490 6950 : DO i = 1, n
1491 6678 : IF (tt(i) < 0.0_dp) THEN
1492 0 : WRITE (*, *) ' dkh_main.F | tt(', i, ') = ', tt(i)
1493 : END IF
1494 :
1495 : ! Calculate some constants
1496 :
1497 6678 : prea = 1/(c_light_au**2)
1498 6678 : con2 = prea + prea
1499 6678 : con = 1.0_dp/prea
1500 :
1501 : ! If T is sufficiently small, use series expansion to avoid
1502 : ! cancellation, otherwise calculate SQRT directly
1503 :
1504 6678 : ev0t(i) = tt(i)
1505 6678 : ratio = tt(i)/c_light_au
1506 6678 : IF (ratio <= 0.02_dp) THEN
1507 1330 : tv1 = tt(i)
1508 1330 : tv2 = -tv1*tt(i)*prea/2.0_dp
1509 1330 : tv3 = -tv2*tt(i)*prea
1510 1330 : tv4 = -tv3*tt(i)*prea*1.25_dp
1511 1330 : ev0t(i) = tv1 + tv2 + tv3 + tv4
1512 : ELSE
1513 5348 : ev0t(i) = con*(SQRT(1.0_dp + con2*tt(i)) - 1.0_dp)
1514 : END IF
1515 6950 : e(i) = ev0t(i) + con
1516 : END DO
1517 :
1518 272 : END SUBROUTINE kintegral_a
1519 :
1520 : ! **************************************************************************************************
1521 : !> \brief ...
1522 : !> \param n ...
1523 : !> \param ev1t ...
1524 : !> \param vt ...
1525 : !> \param pvpt ...
1526 : !> \param aa ...
1527 : !> \param rr ...
1528 : ! **************************************************************************************************
1529 266 : SUBROUTINE even1_a(n, ev1t, vt, pvpt, aa, rr)
1530 :
1531 : !-----------------------------------------------------------------------
1532 : ! -
1533 : ! 1st order DKH-approximation -
1534 : ! -
1535 : ! n in dimension of matrices -
1536 : ! ev1t out even1 output matrix -
1537 : ! vt in potential matrix v in T-space -
1538 : ! pvpt in pvp matrix in T-space -
1539 : ! aa in A-factors (diagonal) -
1540 : ! rr in R-factors (diagonal) -
1541 : ! -
1542 : !-----------------------------------------------------------------------
1543 :
1544 : INTEGER, INTENT(IN) :: n
1545 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: ev1t
1546 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: vt, pvpt
1547 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: aa, rr
1548 :
1549 : INTEGER :: i, j
1550 :
1551 : !-----------------------------------------------------------------------
1552 :
1553 6790 : DO i = 1, n
1554 116356 : DO j = 1, i
1555 109566 : ev1t(i, j) = vt(i, j)*aa(i)*aa(j) + pVpt(i, j)*aa(i)*rr(i)*aa(j)*rr(j)
1556 116090 : ev1t(j, i) = ev1t(i, j)
1557 : END DO
1558 : END DO
1559 :
1560 266 : END SUBROUTINE even1_a
1561 :
1562 : ! **************************************************************************************************
1563 : !> \brief ...
1564 : !> \param n ...
1565 : !> \param pev1tp ...
1566 : !> \param vt ...
1567 : !> \param pvpt ...
1568 : !> \param aa ...
1569 : !> \param rr ...
1570 : !> \param tt ...
1571 : ! **************************************************************************************************
1572 164 : SUBROUTINE peven1p_a(n, pev1tp, vt, pvpt, aa, rr, tt)
1573 :
1574 : !-----------------------------------------------------------------------
1575 : ! -
1576 : ! 1st order DKH-approximation -
1577 : ! -
1578 : ! n in dimension of matrices -
1579 : ! pev1tp out peven1p output matrix -
1580 : ! vt in potential matrix v in T-space -
1581 : ! pvpt in pvp matrix in T-space -
1582 : ! aa in A-factors (diagonal) -
1583 : ! rr in R-factors (diagonal) -
1584 : ! -
1585 : !-----------------------------------------------------------------------
1586 :
1587 : INTEGER, INTENT(IN) :: n
1588 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: pev1tp
1589 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: vt, pvpt
1590 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: aa, rr, tt
1591 :
1592 : INTEGER :: i, j
1593 :
1594 : !-----------------------------------------------------------------------
1595 :
1596 4290 : DO i = 1, n
1597 76446 : DO j = 1, i
1598 : pev1tp(i, j) = 4.0_dp*vt(i, j)*aa(i)*aa(j)*rr(i)*rr(i)*rr(j)*rr(j)*tt(i)*tt(j) + &
1599 72156 : pVpt(i, j)*aa(i)*rr(i)*aa(j)*rr(j)
1600 76282 : pev1tp(j, i) = pev1tp(i, j)
1601 : END DO
1602 : END DO
1603 :
1604 164 : END SUBROUTINE peven1p_a
1605 :
1606 : ! **************************************************************************************************
1607 : !> \brief ...
1608 : !> \param n ...
1609 : !> \param ev2 ...
1610 : !> \param vv ...
1611 : !> \param gg ...
1612 : !> \param aa ...
1613 : !> \param rr ...
1614 : !> \param tt ...
1615 : !> \param e ...
1616 : ! **************************************************************************************************
1617 258 : SUBROUTINE even2c_a(n, ev2, vv, gg, aa, rr, tt, e)
1618 :
1619 : !***********************************************************************
1620 : ! *
1621 : ! Alexander Wolf, last modified: 20.02.2002 - DKH2 *
1622 : ! *
1623 : ! 2nd order DK-approximation ( original DK-transformation with *
1624 : ! U = SQRT(1+W^2) + W ) *
1625 : ! *
1626 : ! Version: 1.1 (20.2.2002) : Usage of SR mat_add included *
1627 : ! 1.0 (6.2.2002) *
1628 : ! Modification history: *
1629 : ! 30.09.2006 Jens Thar: deleted obsolete F77 memory manager *
1630 : ! *
1631 : ! ev2 = 1/2 [W1,O1] *
1632 : ! *
1633 : ! ---- Meaning of Parameters ---- *
1634 : ! *
1635 : ! n in Dimension of matrices *
1636 : ! ev2 out even2 output matrix = final result *
1637 : ! vv in potential v *
1638 : ! gg in pvp *
1639 : ! aa in A-Factors (DIAGONAL) *
1640 : ! rr in R-Factors (DIAGONAL) *
1641 : ! tt in Nonrel. kinetic Energy (DIAGONAL) *
1642 : ! e in Rel. Energy = SQRT(p^2*c^2 + c^4) (DIAGONAL) *
1643 : ! result intermediate result of even2-calculation *
1644 : ! v symmetric (n x n)-matrix containing (A V A) *
1645 : ! pvp symmetric (n x n)-matrix containing (A P V P A) *
1646 : ! vh symmetric (n x n)-matrix containing (A V~ A) *
1647 : ! pvph symmetric (n x n)-matrix containing (A P V~ P A) *
1648 : ! w1o1 W1*O1 (2-component form) *
1649 : ! o1w1 O1*W1 (2-component form) *
1650 : ! *
1651 : !***********************************************************************
1652 :
1653 : INTEGER, INTENT(IN) :: n
1654 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: ev2
1655 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: vv, gg
1656 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: aa, rr, tt, e
1657 :
1658 258 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: o1w1, pvp, pvph, v, vh, w1o1
1659 :
1660 : !-----------------------------------------------------------------------
1661 : ! 1. General Structures and Patterns for DKH2
1662 : !-----------------------------------------------------------------------
1663 :
1664 1032 : ALLOCATE (v(n, n))
1665 1032 : ALLOCATE (pVp(n, n))
1666 1032 : ALLOCATE (vh(n, n))
1667 1032 : ALLOCATE (pVph(n, n))
1668 214066 : v = 0.0_dp
1669 214066 : pVp = 0.0_dp
1670 214066 : vh = 0.0_dp
1671 214066 : pVph = 0.0_dp
1672 214066 : v(1:n, 1:n) = vv(1:n, 1:n)
1673 214066 : vh(1:n, 1:n) = vv(1:n, 1:n)
1674 214066 : pvp(1:n, 1:n) = gg(1:n, 1:n)
1675 214066 : pvph(1:n, 1:n) = gg(1:n, 1:n)
1676 :
1677 214066 : ev2 = 0.0_dp
1678 : ! Calculate v = A V A:
1679 :
1680 258 : CALL mat_axa_a(v, n, aa)
1681 :
1682 : ! Calculate pvp = A P V P A:
1683 :
1684 258 : CALL mat_arxra_a(pvp, n, aa, rr)
1685 :
1686 : ! Calculate vh = A V~ A:
1687 :
1688 258 : CALL mat_1_over_h_a(vh, n, e)
1689 258 : CALL mat_axa_a(vh, n, aa)
1690 :
1691 : ! Calculate pvph = A P V~ P A:
1692 :
1693 258 : CALL mat_1_over_h_a(pvph, n, e)
1694 258 : CALL mat_arxra_a(pvph, n, aa, rr)
1695 :
1696 : ! Create/Initialize necessary matrices:
1697 1032 : ALLOCATE (w1o1(n, n))
1698 1032 : ALLOCATE (o1w1(n, n))
1699 214066 : w1o1 = 0.0_dp
1700 214066 : o1w1 = 0.0_dp
1701 :
1702 : ! Calculate w1o1:
1703 258 : CALL dgemm("N", "N", n, n, n, -1.0_dp, pvph, n, v, n, 0.0_dp, w1o1, n)
1704 258 : CALL mat_muld_a(w1o1, pvph, pvp, n, 1.0_dp, 1.0_dp, tt, rr)
1705 258 : CALL mat_mulm_a(w1o1, vh, v, n, 1.0_dp, 1.0_dp, tt, rr)
1706 258 : CALL dgemm("N", "N", n, n, n, -1.0_dp, vh, n, pvp, n, 1.0_dp, w1o1, n)
1707 : ! Calculate o1w1:
1708 258 : CALL dgemm("N", "N", n, n, n, 1.0_dp, pvp, n, vh, n, 0.0_dp, o1w1, n)
1709 258 : CALL mat_muld_a(o1w1, pvp, pvph, n, -1.0_dp, 1.0_dp, tt, rr)
1710 258 : CALL mat_mulm_a(o1w1, v, vh, n, -1.0_dp, 1.0_dp, tt, rr)
1711 258 : CALL dgemm("N", "N", n, n, n, 1.0_dp, v, n, pvph, n, 1.0_dp, o1w1, n)
1712 : ! Calculate in symmetric pakets
1713 :
1714 : !-----------------------------------------------------------------------
1715 : ! 2. 1/2 [W1,O1] = 1/2 W1O1 - 1/2 O1W1
1716 : !-----------------------------------------------------------------------
1717 :
1718 258 : CALL mat_add(ev2, 0.5_dp, w1o1, -0.5_dp, o1w1, n)
1719 :
1720 : !-----------------------------------------------------------------------
1721 : ! 3. Finish up the stuff!!
1722 : !-----------------------------------------------------------------------
1723 :
1724 258 : DEALLOCATE (v, vh, pvp, pvph, w1o1, o1w1)
1725 :
1726 : ! WRITE (*,*) "CAW: DKH2 with even2c (Alex)"
1727 : ! WRITE (*,*) "!JT: Now available in cp2k"
1728 :
1729 258 : END SUBROUTINE even2c_a
1730 :
1731 : ! **************************************************************************************************
1732 : !> \brief ...
1733 : !> \param n ...
1734 : !> \param ev3 ...
1735 : !> \param e1 ...
1736 : !> \param pe1p ...
1737 : !> \param vv ...
1738 : !> \param gg ...
1739 : !> \param aa ...
1740 : !> \param rr ...
1741 : !> \param tt ...
1742 : !> \param e ...
1743 : ! **************************************************************************************************
1744 164 : SUBROUTINE even3b_a(n, ev3, e1, pe1p, vv, gg, aa, rr, tt, e)
1745 :
1746 : !***********************************************************************
1747 : ! *
1748 : ! Alexander Wolf, last modified: 20.2.2002 - DKH3 *
1749 : ! *
1750 : ! 3rd order DK-approximation (generalised DK-transformation) *
1751 : ! *
1752 : ! Version: 1.1 (20.2.2002) : Usage of SR mat_add included *
1753 : ! 1.0 (7.2.2002) *
1754 : ! *
1755 : ! ev3 = 1/2 [W1,[W1,E1]] *
1756 : ! *
1757 : ! Modification history: *
1758 : ! 30.09.2006 Jens Thar: deleted obsolete F77 memory manager *
1759 : ! *
1760 : ! ---- Meaning of Parameters ---- *
1761 : ! *
1762 : ! n in Dimension of matrices *
1763 : ! ev3 out even3 output matrix = final result *
1764 : ! e1 in E1 = even1-operator *
1765 : ! pe1p in pE1p *
1766 : ! vv in potential v *
1767 : ! gg in pvp *
1768 : ! aa in A-Factors (DIAGONAL) *
1769 : ! rr in R-Factors (DIAGONAL) *
1770 : ! tt in Nonrel. kinetic Energy (DIAGONAL) *
1771 : ! e in Rel. Energy = SQRT(p^2*c^2 + c^4) (DIAGONAL) *
1772 : ! result intermediate result of even2-calculation *
1773 : ! vh symmetric (n x n)-matrix containing (A V~ A) *
1774 : ! pvph symmetric (n x n)-matrix containing (A P V~ P A) *
1775 : ! e1 E1 *
1776 : ! pe1p pE1p *
1777 : ! w1w1 (W1)^2 *
1778 : ! w1e1w1 W1*E1*W1 *
1779 : ! scr_i temporary (n x n)-scratch-matrices (i=1,2) *
1780 : ! *
1781 : !***********************************************************************
1782 :
1783 : INTEGER, INTENT(IN) :: n
1784 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: ev3
1785 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: e1, pe1p, vv, gg
1786 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: aa, rr, tt, e
1787 :
1788 164 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: pvph, scr_1, scr_2, vh, w1e1w1, w1w1
1789 :
1790 : !-----------------------------------------------------------------------
1791 : ! 1. General Structures and Patterns for DKH3
1792 : !-----------------------------------------------------------------------
1793 :
1794 656 : ALLOCATE (vh(n, n))
1795 492 : ALLOCATE (pVph(n, n))
1796 144476 : vh = 0.0_dp
1797 144476 : pVph = 0.0_dp
1798 144476 : vh(1:n, 1:n) = vv(1:n, 1:n)
1799 144476 : pvph(1:n, 1:n) = gg(1:n, 1:n)
1800 :
1801 144476 : ev3 = 0.0_dp
1802 :
1803 : ! Calculate vh = A V~ A:
1804 :
1805 164 : CALL mat_1_over_h_a(vh, n, e)
1806 164 : CALL mat_axa_a(vh, n, aa)
1807 :
1808 : ! Calculate pvph = A P V~ P A:
1809 :
1810 164 : CALL mat_1_over_h_a(pvph, n, e)
1811 164 : CALL mat_arxra_a(pvph, n, aa, rr)
1812 :
1813 : ! Create/Initialize necessary matrices:
1814 492 : ALLOCATE (w1w1(n, n))
1815 492 : ALLOCATE (w1e1w1(n, n))
1816 492 : ALLOCATE (scr_1(n, n))
1817 492 : ALLOCATE (scr_2(n, n))
1818 144476 : w1w1 = 0.0_dp
1819 144476 : w1e1w1 = 0.0_dp
1820 144476 : scr_1 = 0.0_dp
1821 144476 : scr_2 = 0.0_dp
1822 :
1823 : ! Calculate w1w1:
1824 164 : CALL dgemm("N", "N", n, n, n, 1.0_dp, pvph, n, vh, n, 0.0_dp, w1w1, n)
1825 164 : CALL mat_muld_a(w1w1, pvph, pvph, n, -1.0_dp, 1.0_dp, tt, rr)
1826 164 : CALL mat_mulm_a(w1w1, vh, vh, n, -1.0_dp, 1.0_dp, tt, rr)
1827 164 : CALL dgemm("N", "N", n, n, n, 1.0_dp, vh, n, pvph, n, 1.0_dp, w1w1, n)
1828 :
1829 : ! Calculate w1e1w1:
1830 164 : CALL mat_muld_a(scr_1, pvph, pe1p, n, 1.0_dp, 0.0_dp, tt, rr)
1831 164 : CALL dgemm("N", "N", n, n, n, 1.0_dp, vh, n, pe1p, n, 0.0_dp, scr_2, n)
1832 164 : CALL dgemm("N", "N", n, n, n, 1.0_dp, scr_1, n, vh, n, 0.0_dp, w1e1w1, n)
1833 164 : CALL mat_muld_a(w1e1w1, scr_1, pvph, n, -1.0_dp, 1.0_dp, tt, rr)
1834 164 : CALL dgemm("N", "N", n, n, n, -1.0_dp, scr_2, n, vh, n, 1.0_dp, w1e1w1, n)
1835 164 : CALL mat_muld_a(w1e1w1, scr_2, pvph, n, 1.0_dp, 1.0_dp, tt, rr)
1836 :
1837 : !-----------------------------------------------------------------------
1838 : ! 2. ev3 = 1/2 (W1^2)E1 + 1/2 E1(W1^2) - W1E1W1
1839 : !-----------------------------------------------------------------------
1840 :
1841 164 : CALL dgemm("N", "N", n, n, n, 0.5_dp, w1w1, n, e1, n, 0.0_dp, ev3, n)
1842 164 : CALL dgemm("N", "N", n, n, n, 0.5_dp, e1, n, w1w1, n, 1.0_dp, ev3, n)
1843 164 : CALL mat_add2(ev3, 1.0_dp, -1.0_dp, w1e1w1, n)
1844 :
1845 : !-----------------------------------------------------------------------
1846 : ! 3. Finish up the stuff!!
1847 : !-----------------------------------------------------------------------
1848 :
1849 164 : DEALLOCATE (vh, pvph, w1w1, w1e1w1, scr_1, scr_2)
1850 :
1851 : ! WRITE (*,*) "CAW: DKH3 with even3b (Alex)"
1852 : ! WRITE (*,*) "JT: Now available in cp2k"
1853 :
1854 164 : END SUBROUTINE even3b_a
1855 :
1856 : ! **************************************************************************************************
1857 : !> \brief ...
1858 : !> \param n ...
1859 : !> \param ev4 ...
1860 : !> \param e1 ...
1861 : !> \param pe1p ...
1862 : !> \param vv ...
1863 : !> \param gg ...
1864 : !> \param aa ...
1865 : !> \param rr ...
1866 : !> \param tt ...
1867 : !> \param e ...
1868 : ! **************************************************************************************************
1869 0 : SUBROUTINE even4a_a(n, ev4, e1, pe1p, vv, gg, aa, rr, tt, e)
1870 :
1871 : !***********************************************************************
1872 : ! *
1873 : ! Alexander Wolf, last modified: 25.02.2002 -- DKH4 *
1874 : ! *
1875 : ! 4th order DK-approximation (scalar = spin-free) *
1876 : ! *
1877 : ! Version: 1.2 (25.2.2002) : Elegant (short) way of calculation *
1878 : ! included *
1879 : ! 1.1 (20.2.2002) : Usage of SR mat_add included *
1880 : ! 1.0 (8.2.2002) *
1881 : ! *
1882 : ! ev4 = 1/2 [W2,[W1,E1]] + 1/8 [W1,[W1,[W1,O1]]] = *
1883 : ! *
1884 : ! = sum_1 + sum_2 *
1885 : ! *
1886 : ! *
1887 : ! Modification history: *
1888 : ! 30.09.2006 Jens Thar: deleted obsolete F77 memory manager *
1889 : ! *
1890 : ! ---- Meaning of Parameters ---- *
1891 : ! *
1892 : ! n in Dimension of matrices *
1893 : ! ev4 out even4 output matrix = final result *
1894 : ! e1 in E1 *
1895 : ! pe1p in p(E1)p *
1896 : ! vv in potential v *
1897 : ! gg in pvp *
1898 : ! aa in A-Factors (DIAGONAL) *
1899 : ! rr in R-Factors (DIAGONAL) *
1900 : ! tt in Nonrel. kinetic Energy (DIAGONAL) *
1901 : ! e in Rel. Energy = SQRT(p^2*c^2 + c^4) (DIAGONAL) *
1902 : ! v symmetric (n x n)-matrix containing (A V A) *
1903 : ! pvp symmetric (n x n)-matrix containing (A P V P A) *
1904 : ! vh symmetric (n x n)-matrix containing (A V~ A) *
1905 : ! pvph symmetric (n x n)-matrix containing (A P V~ P A) *
1906 : ! w1w1 (W1)^2 *
1907 : ! w1o1 W1*O1 (2-component formulation) *
1908 : ! o1w1 O1*W1 (2-component formulation) *
1909 : ! e1 symmetric (n x n)-matrix containing E1 *
1910 : ! pe1p symmetric (n x n)-matrix containing p(E1)p *
1911 : ! sum_i 2 addends defined above (i=1,2) *
1912 : ! scr_i temporary (n x n)-scratch-matrices (i=1,..,4) *
1913 : ! scrh_i temp. (n x n)-scr.-mat. with energy-denom. (i=1,..,4) *
1914 : ! *
1915 : !***********************************************************************
1916 :
1917 : INTEGER, INTENT(IN) :: n
1918 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: ev4
1919 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: e1, pe1p, vv, gg
1920 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: aa, rr, tt, e
1921 :
1922 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: o1w1, pvp, pvph, scr_1, scr_2, scr_3, &
1923 0 : scr_4, scrh_1, scrh_2, scrh_3, scrh_4, &
1924 0 : sum_1, sum_2, v, vh, w1o1, w1w1
1925 :
1926 : !C-----------------------------------------------------------------------
1927 : !C 1. General Structures and Patterns for DKH4
1928 : !C-----------------------------------------------------------------------
1929 :
1930 0 : ALLOCATE (v(n, n))
1931 0 : ALLOCATE (pVp(n, n))
1932 0 : ALLOCATE (vh(n, n))
1933 0 : ALLOCATE (pVph(n, n))
1934 0 : v = 0.0_dp
1935 0 : pVp = 0.0_dp
1936 0 : vh = 0.0_dp
1937 0 : pVph = 0.0_dp
1938 0 : v(1:n, 1:n) = vv(1:n, 1:n)
1939 0 : vh(1:n, 1:n) = vv(1:n, 1:n)
1940 0 : pvp(1:n, 1:n) = gg(1:n, 1:n)
1941 0 : pvph(1:n, 1:n) = gg(1:n, 1:n)
1942 :
1943 0 : ev4 = 0.0_dp
1944 : ! Calculate v = A V A:
1945 :
1946 0 : CALL mat_axa_a(v, n, aa)
1947 :
1948 : ! Calculate pvp = A P V P A:
1949 :
1950 0 : CALL mat_arxra_a(pvp, n, aa, rr)
1951 :
1952 : ! Calculate vh = A V~ A:
1953 :
1954 0 : CALL mat_1_over_h_a(vh, n, e)
1955 0 : CALL mat_axa_a(vh, n, aa)
1956 :
1957 : ! Calculate pvph = A P V~ P A:
1958 :
1959 0 : CALL mat_1_over_h_a(pvph, n, e)
1960 0 : CALL mat_arxra_a(pvph, n, aa, rr)
1961 :
1962 : ! Create/Initialize necessary matrices:
1963 0 : ALLOCATE (w1w1(n, n))
1964 0 : w1w1 = 0.0_dp
1965 0 : ALLOCATE (w1o1(n, n))
1966 0 : w1o1 = 0.0_dp
1967 0 : ALLOCATE (o1w1(n, n))
1968 0 : o1w1 = 0.0_dp
1969 0 : ALLOCATE (sum_1(n, n))
1970 0 : sum_1 = 0.0_dp
1971 0 : ALLOCATE (sum_2(n, n))
1972 0 : sum_2 = 0.0_dp
1973 0 : ALLOCATE (scr_1(n, n))
1974 0 : scr_1 = 0.0_dp
1975 0 : ALLOCATE (scr_2(n, n))
1976 0 : scr_2 = 0.0_dp
1977 0 : ALLOCATE (scr_3(n, n))
1978 0 : scr_3 = 0.0_dp
1979 0 : ALLOCATE (scr_4(n, n))
1980 0 : scr_4 = 0.0_dp
1981 0 : ALLOCATE (scrh_1(n, n))
1982 0 : scrh_1 = 0.0_dp
1983 0 : ALLOCATE (scrh_2(n, n))
1984 0 : scrh_2 = 0.0_dp
1985 0 : ALLOCATE (scrh_3(n, n))
1986 0 : scrh_3 = 0.0_dp
1987 0 : ALLOCATE (scrh_4(n, n))
1988 0 : scrh_4 = 0.0_dp
1989 :
1990 : ! Calculate w1w1:
1991 0 : CALL dgemm("N", "N", n, n, n, 1.0_dp, pvph, n, vh, n, 0.0_dp, w1w1, n)
1992 0 : CALL mat_muld_a(w1w1, pvph, pvph, n, -1.0_dp, 1.0_dp, tt, rr)
1993 0 : CALL mat_mulm_a(w1w1, vh, vh, n, -1.0_dp, 1.0_dp, tt, rr)
1994 0 : CALL dgemm("N", "N", n, n, n, 1.0_dp, vh, n, pvph, n, 1.0_dp, w1w1, n)
1995 :
1996 : ! Calculate w1o1:
1997 0 : CALL dgemm("N", "N", n, n, n, -1.0_dp, pvph, n, v, n, 0.0_dp, w1o1, n)
1998 0 : CALL mat_muld_a(w1o1, pvph, pvp, n, 1.0_dp, 1.0_dp, tt, rr)
1999 0 : CALL mat_mulm_a(w1o1, vh, v, n, 1.0_dp, 1.0_dp, tt, rr)
2000 0 : CALL dgemm("N", "N", n, n, n, -1.0_dp, vh, n, pvp, n, 1.0_dp, w1o1, n)
2001 : ! Calculate o1w1:
2002 0 : CALL dgemm("N", "N", n, n, n, 1.0_dp, pvp, n, vh, n, 0.0_dp, o1w1, n)
2003 0 : CALL mat_muld_a(o1w1, pvp, pvph, n, -1.0_dp, 1.0_dp, tt, rr)
2004 0 : CALL mat_mulm_a(o1w1, v, vh, n, -1.0_dp, 1.0_dp, tt, rr)
2005 0 : CALL dgemm("N", "N", n, n, n, 1.0_dp, v, n, pvph, n, 1.0_dp, o1w1, n)
2006 :
2007 : !-----------------------------------------------------------------------
2008 : ! 2. sum_1 = 1/2 [W2,[W1,E1]] = 1/2 (W2W1E1 - W2E1W1 - W1E1W2 + E1W1W2)
2009 : !-----------------------------------------------------------------------
2010 :
2011 : ! scr_i & scrh_i for steps 2a (W2W1E1) and 2b (W2E1W1):
2012 :
2013 0 : CALL dgemm("N", "N", n, n, n, 1.0_dp, vh, n, e1, n, 0.0_dp, scr_1, n)
2014 0 : CALL dgemm("N", "N", n, n, n, 1.0_dp, pvph, n, e1, n, 0.0_dp, scr_2, n)
2015 0 : CALL dgemm("N", "N", n, n, n, 1.0_dp, pe1p, n, vh, n, 0.0_dp, scr_3, n)
2016 0 : CALL mat_muld_a(scr_4, pe1p, pvph, n, 1.0_dp, 0.0_dp, tt, rr)
2017 :
2018 0 : CALL mat_muld_a(scrh_1, pvph, pe1p, n, 1.0_dp, 0.0_dp, tt, rr)
2019 0 : CALL mat_1_over_h_a(scrh_1, n, e)
2020 0 : CALL dgemm("N", "N", n, n, n, 1.0_dp, vh, n, pe1p, n, 0.0_dp, scrh_2, n)
2021 0 : CALL mat_1_over_h_a(scrh_2, n, e)
2022 0 : CALL dgemm("N", "N", n, n, n, 1.0_dp, e1, n, pvph, n, 0.0_dp, scrh_3, n)
2023 0 : CALL mat_1_over_h_a(scrh_3, n, e)
2024 0 : CALL dgemm("N", "N", n, n, n, 1.0_dp, e1, n, vh, n, 0.0_dp, scrh_4, n)
2025 0 : CALL mat_1_over_h_a(scrh_4, n, e)
2026 :
2027 : ! 2a) sum_1 = 1/2 W2W1E1 ( [1]-[8] )
2028 :
2029 0 : CALL dgemm("N", "N", n, n, n, 0.5_dp, scrh_1, n, scr_1, n, 0.0_dp, sum_1, n)
2030 0 : CALL mat_muld_a(sum_1, scrh_1, scr_2, n, -0.5_dp, 1.0_dp, tt, rr)
2031 0 : CALL dgemm("N", "N", n, n, n, -0.5_dp, scrh_2, n, scr_1, n, 1.0_dp, sum_1, n)
2032 0 : CALL mat_muld_a(sum_1, scrh_2, scr_2, n, 0.5_dp, 1.0_dp, tt, rr)
2033 0 : CALL dgemm("N", "N", n, n, n, -0.5_dp, scrh_3, n, scr_1, n, 1.0_dp, sum_1, n)
2034 0 : CALL mat_muld_a(sum_1, scrh_3, scr_2, n, 0.5_dp, 1.0_dp, tt, rr)
2035 0 : CALL mat_mulm_a(sum_1, scrh_4, scr_1, n, 0.5_dp, 1.0_dp, tt, rr)
2036 0 : CALL dgemm("N", "N", n, n, n, -0.5_dp, scrh_4, n, scr_2, n, 1.0_dp, sum_1, n)
2037 :
2038 : ! 2b) sum_1 = - 1/2 W2E1W1 (+ sum_1) ( [9]-[16] )
2039 :
2040 0 : CALL mat_muld_a(sum_1, scrh_1, scr_3, n, -0.5_dp, 1.0_dp, tt, rr)
2041 0 : CALL mat_muld_a(sum_1, scrh_1, scr_4, n, 0.5_dp, 1.0_dp, tt, rr)
2042 0 : CALL mat_muld_a(sum_1, scrh_2, scr_3, n, 0.5_dp, 1.0_dp, tt, rr)
2043 0 : CALL mat_muld_a(sum_1, scrh_2, scr_4, n, -0.5_dp, 1.0_dp, tt, rr)
2044 0 : CALL mat_muld_a(sum_1, scrh_3, scr_3, n, 0.5_dp, 1.0_dp, tt, rr)
2045 0 : CALL mat_muld_a(sum_1, scrh_3, scr_4, n, -0.5_dp, 1.0_dp, tt, rr)
2046 0 : CALL dgemm("N", "N", n, n, n, -0.5_dp, scrh_4, n, scr_3, n, 1.0_dp, sum_1, n)
2047 0 : CALL dgemm("N", "N", n, n, n, 0.5_dp, scrh_4, n, scr_4, n, 1.0_dp, sum_1, n)
2048 :
2049 : ! scr_i & scrh_i for steps 2c (W1E1W2) and 2d (E1W1W2):
2050 :
2051 0 : CALL mat_muld_a(scr_1, pvph, pe1p, n, 1.0_dp, 0.0_dp, tt, rr)
2052 0 : CALL dgemm("N", "N", n, n, n, 1.0_dp, vh, n, pe1p, n, 0.0_dp, scr_2, n)
2053 0 : CALL dgemm("N", "N", n, n, n, 1.0_dp, e1, n, pvph, n, 0.0_dp, scr_3, n)
2054 0 : CALL dgemm("N", "N", n, n, n, 1.0_dp, e1, n, vh, n, 0.0_dp, scr_4, n)
2055 :
2056 0 : CALL dgemm("N", "N", n, n, n, 1.0_dp, vh, n, e1, n, 0.0_dp, scrh_1, n)
2057 0 : CALL mat_1_over_h_a(scrh_1, n, e)
2058 0 : CALL dgemm("N", "N", n, n, n, 1.0_dp, pvph, n, e1, n, 0.0_dp, scrh_2, n)
2059 0 : CALL mat_1_over_h_a(scrh_2, n, e)
2060 0 : CALL dgemm("N", "N", n, n, n, 1.0_dp, pe1p, n, vh, n, 0.0_dp, scr_3, n)
2061 0 : CALL mat_1_over_h_a(scrh_3, n, e)
2062 0 : CALL mat_muld_a(scrh_4, pe1p, pvph, n, 1.0_dp, 0.0_dp, tt, rr)
2063 0 : CALL mat_1_over_h_a(scrh_4, n, e)
2064 :
2065 : ! 2c) sum_1 = - 1/2 W1E1W2 (+ sum_1) ( [17]-[24] )
2066 :
2067 0 : CALL dgemm("N", "N", n, n, n, 0.5_dp, scr_1, n, scrh_1, n, 0.0_dp, sum_1, n)
2068 0 : CALL mat_muld_a(sum_1, scr_1, scrh_2, n, -0.5_dp, 1.0_dp, tt, rr)
2069 0 : CALL dgemm("N", "N", n, n, n, -0.5_dp, scr_2, n, scrh_1, n, 1.0_dp, sum_1, n)
2070 0 : CALL mat_muld_a(sum_1, scr_2, scrh_2, n, 0.5_dp, 1.0_dp, tt, rr)
2071 0 : CALL mat_muld_a(sum_1, scr_1, scrh_3, n, -0.5_dp, 1.0_dp, tt, rr)
2072 0 : CALL mat_muld_a(sum_1, scr_1, scrh_4, n, 0.5_dp, 1.0_dp, tt, rr)
2073 0 : CALL mat_muld_a(sum_1, scr_2, scrh_3, n, 0.5_dp, 1.0_dp, tt, rr)
2074 0 : CALL mat_muld_a(sum_1, scr_2, scrh_4, n, -0.5_dp, 1.0_dp, tt, rr)
2075 :
2076 : ! 2d) sum_1 = 1/2 E1W1W2 (+ sum_1) ( [25]-[32] )
2077 :
2078 0 : CALL dgemm("N", "N", n, n, n, -0.5_dp, scr_3, n, scrh_1, n, 0.0_dp, sum_1, n)
2079 0 : CALL mat_muld_a(sum_1, scr_3, scrh_2, n, 0.5_dp, 1.0_dp, tt, rr)
2080 0 : CALL mat_mulm_a(sum_1, scr_4, scrh_1, n, 0.5_dp, 1.0_dp, tt, rr)
2081 0 : CALL dgemm("N", "N", n, n, n, -0.5_dp, scr_4, n, scrh_2, n, 1.0_dp, sum_1, n)
2082 0 : CALL mat_muld_a(sum_1, scr_3, scrh_3, n, 0.5_dp, 1.0_dp, tt, rr)
2083 0 : CALL mat_muld_a(sum_1, scr_3, scrh_4, n, -0.5_dp, 1.0_dp, tt, rr)
2084 0 : CALL dgemm("N", "N", n, n, n, -0.5_dp, scr_4, n, scrh_3, n, 1.0_dp, sum_1, n)
2085 0 : CALL dgemm("N", "N", n, n, n, 0.5_dp, scr_4, n, scrh_4, n, 1.0_dp, sum_1, n)
2086 :
2087 : !-----------------------------------------------------------------------
2088 : ! 3. sum_2 = 1/8 [W1,[W1,[W1,O1]]] =
2089 : !
2090 : ! = 1/8 ( (W1^3)O1 - 3(W1^2)O1W1 + 3 W1O1(W1^2) - O1(W1^3) )
2091 : !-----------------------------------------------------------------------
2092 :
2093 0 : CALL dgemm("N", "N", n, n, n, 0.125_dp, w1w1, n, w1o1, n, 0.0_dp, sum_2, n)
2094 0 : CALL dgemm("N", "N", n, n, n, -0.375_dp, w1w1, n, o1w1, n, 1.0_dp, sum_2, n)
2095 0 : CALL dgemm("N", "N", n, n, n, 0.375_dp, w1o1, n, w1w1, n, 1.0_dp, sum_2, n)
2096 0 : CALL dgemm("N", "N", n, n, n, -0.125_dp, o1w1, n, w1w1, n, 1.0_dp, sum_2, n)
2097 :
2098 : !-----------------------------------------------------------------------
2099 : ! 4. result = sum_1 + sum_2
2100 : !-----------------------------------------------------------------------
2101 :
2102 0 : CALL mat_add(ev4, 1.0_dp, sum_1, 1.0_dp, sum_2, n)
2103 :
2104 : !-----------------------------------------------------------------------
2105 : ! 5. Finish up the stuff!!
2106 : !-----------------------------------------------------------------------
2107 :
2108 0 : DEALLOCATE (v, pvp, vh, pvph, w1w1, w1o1, o1w1, sum_1, sum_2)
2109 0 : DEALLOCATE (scr_1, scr_2, scr_3, scr_4, scrh_1, scrh_2, scrh_3, scrh_4)
2110 :
2111 : ! WRITE (*,*) "CAW: DKH4 with even4a (Alex)"
2112 : ! WRITE (*,*) "JT: Now available in cp2k"
2113 :
2114 0 : END SUBROUTINE even4a_a
2115 :
2116 : !-----------------------------------------------------------------------
2117 : ! -
2118 : ! Matrix routines for DKH-procedure -
2119 : ! Alexander Wolf -
2120 : ! modified: Jens Thar: Mem manager deleted -
2121 : ! This file contains the -
2122 : ! following subroutines: -
2123 : ! 1. mat_1_over_h -
2124 : ! 2. mat_axa -
2125 : ! 3. mat_arxra -
2126 : ! 4. mat_mulm -
2127 : ! 5. mat_muld -
2128 : ! 6. mat_add -
2129 : ! -
2130 : !-----------------------------------------------------------------------
2131 :
2132 : ! **************************************************************************************************
2133 : !> \brief ...
2134 : !> \param p ...
2135 : !> \param n ...
2136 : !> \param e ...
2137 : ! **************************************************************************************************
2138 844 : SUBROUTINE mat_1_over_h_a(p, n, e)
2139 :
2140 : !***********************************************************************
2141 : ! *
2142 : ! 2. SR mat_1_over_h: Transform matrix p into matrix p/(e(i)+e(j)) *
2143 : ! *
2144 : ! p in REAL(:,:) : matrix p *
2145 : ! e in REAL(:) : rel. energy (diagonal) *
2146 : ! n in INTEGER *
2147 : ! *
2148 : !***********************************************************************
2149 :
2150 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: p
2151 : INTEGER, INTENT(IN) :: n
2152 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: e
2153 :
2154 : INTEGER :: i, j
2155 :
2156 21756 : DO i = 1, n
2157 717084 : DO j = 1, n
2158 716240 : p(i, j) = p(i, j)/(e(i) + e(j))
2159 : END DO
2160 : END DO
2161 :
2162 844 : END SUBROUTINE mat_1_over_h_a
2163 :
2164 : ! **************************************************************************************************
2165 : !> \brief ...
2166 : !> \param p ...
2167 : !> \param n ...
2168 : !> \param a ...
2169 : ! **************************************************************************************************
2170 680 : SUBROUTINE mat_axa_a(p, n, a)
2171 :
2172 : !C***********************************************************************
2173 : !C *
2174 : !C 3. SR mat_axa: Transform matrix p into matrix a*p*a *
2175 : !C *
2176 : !C p in REAL(:,:): matrix p *
2177 : !C a in REAL(:) : A-factors (diagonal) *
2178 : !CJT n in INTEGER : dimension of matrix p *
2179 : !C *
2180 : !C***********************************************************************
2181 :
2182 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: p
2183 : INTEGER, INTENT(IN) :: n
2184 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: a
2185 :
2186 : INTEGER :: i, j
2187 :
2188 17466 : DO i = 1, n
2189 572608 : DO j = 1, n
2190 571928 : p(i, j) = p(i, j)*a(i)*a(j)
2191 : END DO
2192 : END DO
2193 :
2194 680 : END SUBROUTINE mat_axa_a
2195 :
2196 : ! **************************************************************************************************
2197 : !> \brief ...
2198 : !> \param p ...
2199 : !> \param n ...
2200 : !> \param a ...
2201 : !> \param r ...
2202 : ! **************************************************************************************************
2203 680 : SUBROUTINE mat_arxra_a(p, n, a, r)
2204 :
2205 : !C***********************************************************************
2206 : !C *
2207 : !C 4. SR mat_arxra: Transform matrix p into matrix a*r*p*r*a *
2208 : !C *
2209 : !C p in REAL(:,:) : matrix p *
2210 : !C a in REAL(:) : A-factors (diagonal) *
2211 : !C r in REAL(:) : R-factors (diagonal) *
2212 : !C n in INTEGER : dimension of matrix p *
2213 : !C *
2214 : !C***********************************************************************
2215 :
2216 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: p
2217 : INTEGER, INTENT(IN) :: n
2218 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: a, r
2219 :
2220 : INTEGER :: i, j
2221 :
2222 17466 : DO i = 1, n
2223 572608 : DO j = 1, n
2224 571928 : p(i, j) = p(i, j)*a(i)*a(j)*r(i)*r(j)
2225 : END DO
2226 : END DO
2227 :
2228 680 : END SUBROUTINE mat_arxra_a
2229 :
2230 : ! **************************************************************************************************
2231 : !> \brief ...
2232 : !> \param p ...
2233 : !> \param q ...
2234 : !> \param r ...
2235 : !> \param n ...
2236 : !> \param alpha ...
2237 : !> \param beta ...
2238 : !> \param t ...
2239 : !> \param rr ...
2240 : ! **************************************************************************************************
2241 680 : SUBROUTINE mat_mulm_a(p, q, r, n, alpha, beta, t, rr)
2242 :
2243 : !C***********************************************************************
2244 : !C *
2245 : !C 5. SR mat_mulm: Multiply matrices according to: *
2246 : !C *
2247 : !C p = alpha*q*(..P^2..)*r + beta*p *
2248 : !C *
2249 : !C p out REAL(:,:): matrix p *
2250 : !C q in REAL(:,:): matrix q *
2251 : !C r in REAL(:,.): matrix r *
2252 : !C n in INTEGER : dimension n of matrices *
2253 : !C alpha in REAL(dp) : *
2254 : !C beta in REAL(dp) : *
2255 : !C t in REAL(:) : non-rel. kinetic energy (diagonal) *
2256 : !C rr in REAL(:) : R-factors (diagonal) *
2257 : !C *
2258 : !C***********************************************************************
2259 :
2260 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: p, q, r
2261 : INTEGER, INTENT(IN) :: n
2262 : REAL(KIND=dp), INTENT(IN) :: alpha, beta
2263 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: t, rr
2264 :
2265 : INTEGER :: i, j
2266 680 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: qtemp
2267 :
2268 2720 : ALLOCATE (qtemp(n, n))
2269 :
2270 17466 : DO i = 1, n
2271 572608 : DO j = 1, n
2272 571928 : qtemp(i, j) = q(i, j)*2.0_dp*t(j)*rr(j)*rr(j)
2273 : END DO
2274 : END DO
2275 :
2276 680 : CALL dgemm("N", "N", n, n, n, alpha, qtemp, n, r, n, beta, p, n)
2277 :
2278 680 : DEALLOCATE (qtemp)
2279 :
2280 680 : END SUBROUTINE mat_mulm_a
2281 :
2282 : ! **************************************************************************************************
2283 : !> \brief ...
2284 : !> \param p ...
2285 : !> \param q ...
2286 : !> \param r ...
2287 : !> \param n ...
2288 : !> \param alpha ...
2289 : !> \param beta ...
2290 : !> \param t ...
2291 : !> \param rr ...
2292 : ! **************************************************************************************************
2293 1172 : SUBROUTINE mat_muld_a(p, q, r, n, alpha, beta, t, rr)
2294 :
2295 : !C***********************************************************************
2296 : !C *
2297 : !C 16. SR mat_muld: Multiply matrices according to: *
2298 : !C *
2299 : !C p = alpha*q*(..1/P^2..)*r + beta*p *
2300 : !C *
2301 : !C p out REAL(:,:): matrix p *
2302 : !C q in REAL(:,:): matrix q *
2303 : !C r in REAL(:,:): matrix r *
2304 : !C n in INTEGER : Dimension of all matrices *
2305 : !C alpha in REAL(dp) : *
2306 : !C beta in REAL(dp) : *
2307 : !C t in REAL(:) : non-rel. kinetic energy (diagonal) *
2308 : !C rr in REAL(:) : R-factors (diagonal) *
2309 : !C *
2310 : !C***********************************************************************
2311 :
2312 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: p, q, r
2313 : INTEGER, INTENT(IN) :: n
2314 : REAL(KIND=dp), INTENT(IN) :: alpha, beta
2315 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: t, rr
2316 :
2317 : INTEGER :: i, j
2318 1172 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: qtemp
2319 :
2320 4688 : ALLOCATE (qtemp(n, n))
2321 :
2322 30336 : DO i = 1, n
2323 1006036 : DO j = 1, n
2324 1004864 : qtemp(i, j) = q(i, j)*0.5_dp/(t(j)*rr(j)*rr(j))
2325 : END DO
2326 : END DO
2327 :
2328 1172 : CALL dgemm("N", "N", n, n, n, alpha, qtemp, n, r, n, beta, p, n)
2329 :
2330 1172 : DEALLOCATE (qtemp)
2331 :
2332 1172 : END SUBROUTINE mat_muld_a
2333 :
2334 : ! **************************************************************************************************
2335 : !> \brief ...
2336 : !> \param p ...
2337 : !> \param alpha ...
2338 : !> \param beta ...
2339 : !> \param r ...
2340 : !> \param n ...
2341 : ! **************************************************************************************************
2342 852 : SUBROUTINE mat_add2(p, alpha, beta, r, n)
2343 :
2344 : !C***********************************************************************
2345 : !C *
2346 : !C 19. SR mat_add: Add two matrices of the same size according to: *
2347 : !C *
2348 : !C p = alpha*p + beta*r *
2349 : !C *
2350 : !C and store them in the first *
2351 : !C p out REAL(:,:) : matrix p *
2352 : !C r in REAL(:,:) : matrix r *
2353 : !C alpha in REAL(dp) *
2354 : !C beta in REAL(dp) *
2355 : !C *
2356 : !C Matrix p must already exist before calling this SR!! *
2357 : !C *
2358 : !C [written by: Alexander Wolf, 20.2.2002, v1.0] *
2359 : !C *
2360 : !C***********************************************************************
2361 :
2362 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: p
2363 : REAL(KIND=dp), INTENT(IN) :: alpha, beta
2364 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: r
2365 : INTEGER, INTENT(IN) :: n
2366 :
2367 : INTEGER :: i, j
2368 :
2369 : !C Add matrices:
2370 :
2371 21958 : DO i = 1, n
2372 722416 : DO j = 1, n
2373 721564 : p(i, j) = alpha*p(i, j) + beta*r(i, j)
2374 : END DO
2375 : END DO
2376 :
2377 852 : END SUBROUTINE mat_add2
2378 :
2379 : ! **************************************************************************************************
2380 : !> \brief ...
2381 : !> \param p ...
2382 : !> \param alpha ...
2383 : !> \param q ...
2384 : !> \param beta ...
2385 : !> \param r ...
2386 : !> \param n ...
2387 : ! **************************************************************************************************
2388 258 : SUBROUTINE mat_add(p, alpha, q, beta, r, n)
2389 :
2390 : !C***********************************************************************
2391 : !C *
2392 : !C 19. SR mat_add: Add two matrices of the same size according to: *
2393 : !C *
2394 : !C p = alpha*q + beta*r *
2395 : !C *
2396 : !C p out REAL(:,:) : matrix p *
2397 : !C q in REAL(:,:) : matrix q *
2398 : !C r in REAL(:,:) : matrix r *
2399 : !C alpha in REAL(dp) *
2400 : !C beta in REAL(dp) *
2401 : !C *
2402 : !C Matrix p must already exist before calling this SR!! *
2403 : !C *
2404 : !C [written by: Alexander Wolf, 20.2.2002, v1.0] *
2405 : !C *
2406 : !C***********************************************************************
2407 :
2408 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: p
2409 : REAL(KIND=dp), INTENT(IN) :: alpha
2410 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: q
2411 : REAL(KIND=dp), INTENT(IN) :: beta
2412 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: r
2413 : INTEGER, INTENT(IN) :: n
2414 :
2415 : INTEGER :: i, j
2416 :
2417 : ! Add matrices:
2418 :
2419 6588 : DO i = 1, n
2420 214066 : DO j = 1, n
2421 213808 : p(i, j) = alpha*q(i, j) + beta*r(i, j)
2422 : END DO
2423 : END DO
2424 :
2425 258 : END SUBROUTINE mat_add
2426 :
2427 : ! **************************************************************************************************
2428 : !> \brief ...
2429 : !> \param W ...
2430 : !> \param B ...
2431 : !> \param C ...
2432 : !> \param N ...
2433 : !> \param H ...
2434 : ! **************************************************************************************************
2435 1088 : SUBROUTINE TRSM(W, B, C, N, H)
2436 :
2437 : REAL(KIND=dp), DIMENSION(:, :) :: W, B, C
2438 : INTEGER :: N
2439 : REAL(KIND=dp), DIMENSION(:, :) :: H
2440 :
2441 : INTEGER :: I, IJ, J, K, L
2442 :
2443 : !C
2444 : !C TRANSFORM SYMMETRIC matrix A by UNITARY TRANSFORMATION
2445 : !C IN B. RESULT IS IN C
2446 : !C
2447 : !CAW C = B^{dagger} * A * B
2448 :
2449 1088 : IJ = 0
2450 27800 : DO I = 1, N
2451 474568 : DO J = 1, I
2452 446768 : IJ = IJ + 1
2453 446768 : C(I, J) = 0.0_dp
2454 446768 : C(J, I) = 0.0_dp
2455 446768 : H(I, J) = 0.0_dp
2456 473480 : H(J, I) = 0.0_dp
2457 : END DO
2458 : END DO
2459 27800 : DO I = 1, N
2460 894624 : DO L = 1, N
2461 32014168 : DO K = 1, N
2462 31987456 : H(I, L) = B(K, I)*W(K, L) + H(I, L)
2463 : END DO
2464 : END DO
2465 : END DO
2466 :
2467 : IJ = 0
2468 27800 : DO I = 1, N
2469 474568 : DO J = 1, I
2470 16440496 : IJ = IJ + 1
2471 16467208 : DO L = 1, N
2472 15993728 : C(I, J) = H(I, L)*B(L, J) + C(I, J)
2473 16440496 : C(J, I) = C(I, J)
2474 : END DO
2475 : END DO
2476 : END DO
2477 :
2478 1088 : END SUBROUTINE TRSM
2479 :
2480 : ! **************************************************************************************************
2481 : !> \brief ...
2482 : !> \param matrix_t_pgf ...
2483 : !> \param n ...
2484 : !> \param eig ...
2485 : !> \param ew ...
2486 : !> \param matrix_sinv_pgf ...
2487 : !> \param aux ...
2488 : !> \param ic ...
2489 : ! **************************************************************************************************
2490 272 : SUBROUTINE dkh_diag(matrix_t_pgf, n, eig, ew, matrix_sinv_pgf, aux, ic)
2491 :
2492 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: matrix_t_pgf
2493 : INTEGER, INTENT(IN) :: n
2494 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: eig
2495 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: ew
2496 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: matrix_sinv_pgf
2497 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: aux
2498 : INTEGER :: ic
2499 :
2500 : INTEGER :: n2
2501 :
2502 223656 : eig = 0.0_dp
2503 223656 : aux = 0.0_dp
2504 :
2505 70184 : CALL dgemm("N", "N", n, n, n, 1.0_dp, matrix_t_pgf, n, matrix_sinv_pgf, n, 0.0_dp, eig, n)
2506 :
2507 223656 : aux = 0.0_dp
2508 :
2509 272 : CALL dgemm("T", "N", n, n, n, 1.0_dp, matrix_sinv_pgf, n, eig, n, 0.0_dp, aux, n)
2510 :
2511 272 : n2 = 3*n - 1
2512 :
2513 272 : CALL JACOB2(AUX, EIG, EW, N, IC)
2514 :
2515 272 : END SUBROUTINE dkh_diag
2516 :
2517 : ! **************************************************************************************************
2518 : !> \brief ...
2519 : !> \param sogt ...
2520 : !> \param eigv ...
2521 : !> \param eigw ...
2522 : !> \param n ...
2523 : !> \param ic ...
2524 : ! **************************************************************************************************
2525 272 : SUBROUTINE JACOB2(sogt, eigv, eigw, n, ic)
2526 :
2527 : INTEGER, INTENT(IN) :: n
2528 : REAL(KIND=dp), DIMENSION(n), INTENT(OUT) :: eigw
2529 : REAL(KIND=dp), DIMENSION(n, n), INTENT(OUT) :: eigv
2530 : REAL(KIND=dp), DIMENSION(n, n), INTENT(INOUT) :: sogt
2531 : INTEGER, INTENT(IN) :: ic
2532 :
2533 : INTEGER :: i, il, im, ind, j, k, l, ll, m, mm
2534 : REAL(KIND=dp) :: cost, cost2, ext_norm, sincs, sint, &
2535 : sint2, thr, thr_min, tol, u1, x, xy, y
2536 :
2537 272 : tol = 1.0E-15
2538 272 : ext_norm = 0.0_dp
2539 272 : u1 = REAL(n, KIND=dp)
2540 6950 : DO i = 1, n
2541 6678 : eigv(i, i) = 1.0_dp
2542 6678 : eigw(i) = sogt(i, i)
2543 118642 : DO j = 1, i
2544 118370 : IF (i /= j) THEN
2545 105014 : eigv(i, j) = 0.0_dp
2546 105014 : eigv(j, i) = 0.0_dp
2547 105014 : ext_norm = ext_norm + sogt(i, j)*sogt(i, j)
2548 : END IF
2549 : END DO
2550 : END DO
2551 :
2552 272 : IF (ext_norm > 0.0_dp) THEN
2553 268 : ext_norm = SQRT(2.0_dp*ext_norm)
2554 268 : thr_min = ext_norm*tol/u1
2555 268 : ind = 0
2556 268 : thr = ext_norm
2557 :
2558 : DO
2559 3866 : thr = thr/u1
2560 : DO
2561 8536 : l = 1
2562 : DO
2563 194558 : m = l + 1
2564 2960976 : DO
2565 3155534 : IF ((ABS(sogt(m, l)) - thr) >= 0.0_dp) THEN
2566 261490 : ind = 1
2567 261490 : x = 0.5_dp*(eigw(l) - eigw(m))
2568 261490 : y = -sogt(m, l)/SQRT(sogt(m, l)*sogt(m, l) + x*x)
2569 261490 : IF (x < 0.0_dp) y = -y
2570 :
2571 261490 : IF (y > 1.0_dp) y = 1.0_dp
2572 261490 : IF (y < -1.0_dp) y = -1.0_dp
2573 261490 : xy = 1.0_dp - y*y
2574 261490 : sint = y/SQRT(2.0_dp*(1.0_dp + SQRT(xy)))
2575 261490 : sint2 = sint*sint
2576 261490 : cost2 = 1.0_dp - sint2
2577 261490 : cost = SQRT(cost2)
2578 261490 : sincs = sint*cost
2579 :
2580 9463576 : DO i = 1, n
2581 9202086 : IF ((i - m) /= 0) THEN
2582 8940596 : IF ((i - m) < 0) THEN
2583 : im = m
2584 : mm = i
2585 : ELSE
2586 3037970 : im = i
2587 3037970 : mm = m
2588 : END IF
2589 8940596 : IF ((i - l) /= 0) THEN
2590 8679106 : IF ((i - l) < 0) THEN
2591 : il = l
2592 : ll = i
2593 : ELSE
2594 5310436 : il = i
2595 5310436 : ll = l
2596 : END IF
2597 8679106 : x = sogt(il, ll)*cost - sogt(im, mm)*sint
2598 8679106 : sogt(im, mm) = sogt(il, ll)*sint + sogt(im, mm)*cost
2599 8679106 : sogt(il, ll) = x
2600 : END IF
2601 : END IF
2602 :
2603 9202086 : x = eigv(i, l)*cost - eigv(i, m)*sint
2604 9202086 : eigv(i, m) = eigv(i, l)*sint + eigv(i, m)*cost
2605 9463576 : eigv(i, l) = x
2606 : END DO
2607 :
2608 261490 : x = 2.0_dp*sogt(m, l)*sincs
2609 261490 : y = eigw(l)*cost2 + eigw(m)*sint2 - x
2610 261490 : x = eigw(l)*sint2 + eigw(m)*cost2 + x
2611 261490 : sogt(m, l) = (eigw(l) - eigw(m))*sincs + sogt(m, l)*(cost2 - sint2)
2612 261490 : eigw(l) = y
2613 261490 : eigw(m) = x
2614 : END IF
2615 3155534 : IF ((m - n) == 0) EXIT
2616 2960976 : m = m + 1
2617 : END DO
2618 194558 : IF ((l - m + 1) == 0) EXIT
2619 8536 : l = l + 1
2620 : END DO
2621 8536 : IF ((ind - 1) /= 0.0_dp) EXIT
2622 3866 : ind = 0
2623 : END DO
2624 3866 : IF ((thr - thr_min) <= 0.0_dp) EXIT
2625 : END DO
2626 : END IF
2627 :
2628 272 : IF (ic /= 0) THEN
2629 0 : DO i = 1, n
2630 0 : DO j = 1, n
2631 0 : IF ((eigw(i) - eigw(j)) > 0.0_dp) THEN
2632 0 : x = eigw(i)
2633 0 : eigw(i) = eigw(j)
2634 0 : eigw(j) = x
2635 0 : DO k = 1, n
2636 0 : y = eigv(k, i)
2637 0 : eigv(k, i) = eigv(k, j)
2638 0 : eigv(k, j) = y
2639 : END DO
2640 : END IF
2641 : END DO
2642 : END DO
2643 :
2644 : END IF
2645 :
2646 272 : END SUBROUTINE JACOB2
2647 :
2648 : ! **************************************************************************************************
2649 : !> \brief ...
2650 : !> \param n ...
2651 : !> \param matrix_s_pgf ...
2652 : !> \param matrix_sinv_pgf ...
2653 : ! **************************************************************************************************
2654 272 : SUBROUTINE SOG(n, matrix_s_pgf, matrix_sinv_pgf)
2655 :
2656 : INTEGER :: n
2657 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: matrix_s_pgf
2658 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: matrix_sinv_pgf
2659 :
2660 : INTEGER :: i, j, jn, k
2661 : REAL(KIND=dp) :: diag_s, row_sum, scalar
2662 272 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: a
2663 272 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: g
2664 :
2665 : ! SUBROUTINE TO CALCULATE TRANSFORMATION TO SCHMIDT-
2666 : ! ORTHOGONALIZED BASIS.
2667 : ! sinv-1*matrix_s_pgf*sinv = "orthogonal matrix"
2668 : ! n dimension of matrices
2669 : ! matrix_s_pgf original overlap matrix
2670 : ! matrix_sinv_pgf new overlap matrix
2671 : ! g scratch
2672 : ! a scratch
2673 :
2674 816 : ALLOCATE (a(n))
2675 1088 : ALLOCATE (g(n, n))
2676 :
2677 6950 : DO jn = 1, n
2678 6678 : diag_s = matrix_s_pgf(jn, jn)
2679 6678 : g(jn, jn) = 1.0_dp
2680 :
2681 6678 : IF (jn /= 1) THEN
2682 111420 : DO j = 1, jn - 1
2683 : scalar = 0.0_dp
2684 1400594 : DO i = 1, j
2685 1400594 : scalar = scalar + matrix_s_pgf(i, jn)*g(i, j)
2686 : END DO
2687 105014 : diag_s = diag_s - scalar*scalar
2688 111420 : a(j) = scalar
2689 : END DO
2690 :
2691 111420 : DO j = 1, jn - 1
2692 : row_sum = 0.0_dp
2693 1400594 : DO k = j, jn - 1
2694 1400594 : row_sum = row_sum + a(k)*g(j, k)
2695 : END DO
2696 111420 : g(j, jn) = -row_sum
2697 : END DO
2698 : END IF
2699 :
2700 6678 : diag_s = 1.0_dp/SQRT(diag_s)
2701 118642 : DO i = 1, jn
2702 118370 : g(i, jn) = g(i, jn)*diag_s
2703 : END DO
2704 : END DO
2705 :
2706 6950 : DO j = 1, n
2707 118642 : DO i = 1, j
2708 111692 : matrix_sinv_pgf(j, i) = 0.0_dp
2709 118370 : matrix_sinv_pgf(i, j) = g(i, j)
2710 : END DO
2711 : END DO
2712 :
2713 272 : DEALLOCATE (a, g)
2714 :
2715 272 : END SUBROUTINE SOG
2716 :
2717 : END MODULE dkh_main
|
