(git:ccc2433)
dkh_main.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 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_upper_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  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  CALL timeset(routinen, handle)
116  NULLIFY (blacs_env)
117 
118  !-----------------------------------------------------------------------
119  ! Construct the matrix structure
120  !-----------------------------------------------------------------------
121 
122  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  ncol_global=n)
127 
128  !-----------------------------------------------------------------------
129  ! Allocate some matrices
130  !-----------------------------------------------------------------------
131 
132  ALLOCATE (e(n))
133  ALLOCATE (aa(n))
134  ALLOCATE (rr(n))
135  ALLOCATE (tt(n))
136  ALLOCATE (ev0t(n))
137 
138  CALL cp_fm_create(matrix_eig, matrix_full)
139  CALL cp_fm_create(matrix_aux, matrix_full)
140  CALL cp_fm_create(matrix_aux2, matrix_full)
141  CALL cp_fm_create(matrix_rev, matrix_full)
142  CALL cp_fm_create(matrix_se, matrix_full)
143  CALL cp_fm_create(matrix_ev1, matrix_full)
144  CALL cp_fm_create(matrix_ev2, matrix_full)
145  CALL cp_fm_create(matrix_sinv, matrix_full)
146  CALL cp_fm_create(matrix_ev3, matrix_full)
147  CALL cp_fm_create(matrix_ev4, matrix_full)
148  CALL cp_fm_create(matrix_pe1p, matrix_full)
149 
150  !-----------------------------------------------------------------------
151  ! Now with Cholesky decomposition
152  !-----------------------------------------------------------------------
153 
154  CALL cp_fm_to_fm(matrix_s, matrix_sinv)
155  CALL cp_fm_cholesky_decompose(matrix_sinv, n)
156 
157  !-----------------------------------------------------------------------
158  ! Calculate matrix representation from nonrelativistic T matrix
159  !-----------------------------------------------------------------------
160 
161  CALL cp_fm_cholesky_reduce(matrix_t, matrix_sinv)
162  CALL cp_fm_syevd(matrix_t, matrix_eig, tt)
163 
164  !-----------------------------------------------------------------------
165  ! Calculate kinetic part of Hamiltonian in T-basis
166  !-----------------------------------------------------------------------
167 
168  CALL kintegral(n, ev0t, tt, e)
169 
170  !-----------------------------------------------------------------------
171  ! Calculate reverse transformation matrix revt
172  !-----------------------------------------------------------------------
173 
174  CALL cp_fm_to_fm(matrix_eig, matrix_rev)
175  CALL cp_fm_triangular_multiply(matrix_sinv, matrix_rev, transpose_tr=.true.)
176 
177  !-----------------------------------------------------------------------
178  ! Calculate kinetic part of the Hamiltonian
179  !-----------------------------------------------------------------------
180 
181  CALL cp_fm_to_fm(matrix_rev, matrix_aux)
182  CALL cp_fm_column_scale(matrix_aux, ev0t)
183  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  DO i = 1, n
191  aa(i) = sqrt((c_light_au**2 + e(i))/(2.0_dp*e(i)))
192  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  CALL cp_fm_cholesky_reduce(matrix_v, matrix_sinv)
200  CALL cp_fm_upper_to_full(matrix_v, matrix_aux)
201  CALL parallel_gemm("T", "N", n, n, n, 1.0_dp, matrix_eig, matrix_v, 0.0_dp, matrix_aux)
202  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  CALL cp_fm_cholesky_reduce(matrix_pvp, matrix_sinv)
209  CALL cp_fm_upper_to_full(matrix_pvp, matrix_aux)
210  CALL parallel_gemm("T", "N", n, n, n, 1.0_dp, matrix_eig, matrix_pvp, 0.0_dp, matrix_aux)
211  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  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  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  IF (dkh_order .GE. 3) THEN
230  CALL peven1p(n, matrix_pe1p, matrix_v, matrix_pvp, matrix_aux, matrix_aux2, aa, rr, tt)
231  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  CALL parallel_gemm("N", "N", n, n, n, 1.0_dp, matrix_rev, matrix_ev3, 0.0_dp, matrix_aux)
238  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  IF (dkh_order .GE. 4) THEN
245  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  CALL parallel_gemm("N", "N", n, n, n, 1.0_dp, matrix_rev, matrix_ev1, 0.0_dp, matrix_aux)
264  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  CALL parallel_gemm("N", "N", n, n, n, 1.0_dp, matrix_rev, matrix_ev2, 0.0_dp, matrix_aux)
271  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  CALL cp_fm_scale_and_add(1.0_dp, matrix_ev1, 1.0_dp, matrix_ev2)
278  CALL cp_fm_upper_to_full(matrix_ev1, matrix_aux)
279  CALL cp_fm_to_fm(matrix_ev1, matrix_v)
280  IF (dkh_order .GE. 3) THEN
281  CALL cp_fm_scale_and_add(1.0_dp, matrix_v, 1.0_dp, matrix_ev3)
282  IF (dkh_order .GE. 4) THEN
283  CALL cp_fm_scale_and_add(1.0_dp, matrix_v, 1.0_dp, matrix_ev4)
284  END IF
285  END IF
286 
287  CALL cp_fm_release(matrix_eig)
288  CALL cp_fm_release(matrix_aux)
289  CALL cp_fm_release(matrix_aux2)
290  CALL cp_fm_release(matrix_rev)
291  CALL cp_fm_release(matrix_se)
292  CALL cp_fm_release(matrix_ev1)
293  CALL cp_fm_release(matrix_ev2)
294  CALL cp_fm_release(matrix_sinv)
295  CALL cp_fm_release(matrix_ev3)
296  CALL cp_fm_release(matrix_ev4)
297  CALL cp_fm_release(matrix_pe1p)
298 
299  CALL cp_fm_struct_release(matrix_full)
300 
301  DEALLOCATE (ev0t, e, aa, rr, tt)
302 
303  CALL timestop(handle)
304 
305  END SUBROUTINE dkh_full_transformation
306 
307 ! **************************************************************************************************
308 !> \brief ...
309 !> \param n ...
310 !> \param ev0t ...
311 !> \param tt ...
312 !> \param e ...
313 ! **************************************************************************************************
314  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  prea = 1/(c_light_au**2)
325  con2 = prea + prea
326  con = 1.0_dp/prea
327 
328  DO i = 1, n
329  IF (tt(i) .LT. 0.0_dp) THEN
330  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  ev0t(i) = tt(i)
337  ratio = tt(i)/c_light_au
338  IF (ratio .LE. 0.02_dp) THEN
339  tv1 = tt(i)
340  tv2 = -tv1*tt(i)*prea*0.5_dp
341  tv3 = -tv2*tt(i)*prea
342  tv4 = -tv3*tt(i)*prea*1.25_dp
343  ev0t(i) = tv1 + tv2 + tv3 + tv4
344  ELSE
345  ev0t(i) = con*(sqrt(1.0_dp + con2*tt(i)) - 1.0_dp)
346  END IF
347  e(i) = ev0t(i) + con
348  END DO
349 
350  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  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  CALL cp_fm_to_fm(matrix_v, matrix_aux)
368  CALL cp_fm_column_scale(matrix_aux, aa)
369  CALL cp_fm_transpose(matrix_aux, matrix_ev1)
370  CALL cp_fm_column_scale(matrix_ev1, aa)
371 
372  CALL cp_fm_to_fm(matrix_pvp, matrix_aux)
373  CALL cp_fm_column_scale(matrix_aux, aa)
374  CALL cp_fm_column_scale(matrix_aux, rr)
375  CALL cp_fm_transpose(matrix_aux, matrix_aux2)
376  CALL cp_fm_column_scale(matrix_aux2, aa)
377  CALL cp_fm_column_scale(matrix_aux2, rr)
378 
379  CALL cp_fm_scale_and_add(1.0_dp, matrix_ev1, 1.0_dp, matrix_aux2)
380 
381  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  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  INTEGER, DIMENSION(:), POINTER :: row_indices
404  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  DO i = 1, n
410  vec_ar(i) = aa(i)*rr(i)
411  vec_arrt(i) = vec_ar(i)*rr(i)*tt(i)
412  END DO
413 
414  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  ncol_global=1)
419 
420  CALL cp_fm_create(vec_a, vec_full)
421 
422  CALL cp_fm_get_info(matrix_v, nrow_local=nrow_local, &
423  row_indices=row_indices)
424 
425  DO i = 1, nrow_local
426  vec_a%local_data(i, 1) = vec_arrt(row_indices(i))
427  END DO
428 
429  CALL cp_fm_syrk('U', 'N', 1, 1.0_dp, vec_a, 1, 1, 0.0_dp, matrix_aux)
430  CALL cp_fm_upper_to_full(matrix_aux, matrix_aux2)
431  CALL cp_fm_schur_product(matrix_v, matrix_aux, matrix_pe1p)
432 
433  DO i = 1, nrow_local
434  vec_a%local_data(i, 1) = vec_ar(row_indices(i))
435  END DO
436 
437  CALL cp_fm_syrk('U', 'N', 1, 1.0_dp, vec_a, 1, 1, 0.0_dp, matrix_aux)
438  CALL cp_fm_upper_to_full(matrix_aux, matrix_aux2)
439  CALL cp_fm_schur_product(matrix_pvp, matrix_aux, matrix_aux2)
440 
441  CALL cp_fm_scale_and_add(4.0_dp, matrix_pe1p, 1.0_dp, matrix_aux2)
442 
443  CALL cp_fm_release(vec_a)
444  CALL cp_fm_struct_release(vec_full)
445 
446  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  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  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  ncol_global=n)
497 
498  CALL cp_fm_create(matrix_aux2, matrix_full)
499  CALL cp_fm_create(matrix_ava, matrix_full)
500  CALL cp_fm_create(matrix_avva, matrix_full)
501  CALL cp_fm_create(matrix_apvpa, matrix_full)
502  CALL cp_fm_create(matrix_apvvpa, matrix_full)
503 
504  CALL cp_fm_to_fm(matrix_v, matrix_ava)
505  CALL cp_fm_to_fm(matrix_v, matrix_avva)
506  CALL cp_fm_to_fm(matrix_pvp, matrix_apvpa)
507  CALL cp_fm_to_fm(matrix_pvp, matrix_apvvpa)
508 
509  ! Calculate v = A V A:
510 
511  CALL mat_axa(matrix_v, matrix_ava, n, aa, matrix_aux)
512 
513  ! Calculate pvp = A P V P A:
514 
515  CALL mat_arxra(matrix_pvp, matrix_apvpa, n, aa, rr, matrix_aux)
516 
517  ! Calculate vh = A V~ A:
518 
519  CALL mat_1_over_h(matrix_v, matrix_avva, e, matrix_aux)
520  CALL cp_fm_to_fm(matrix_avva, matrix_aux2)
521  CALL mat_axa(matrix_aux2, matrix_avva, n, aa, matrix_aux)
522 
523  ! Calculate pvph = A P V~ P A:
524 
525  CALL mat_1_over_h(matrix_pvp, matrix_apvvpa, e, matrix_aux)
526  CALL cp_fm_to_fm(matrix_apvvpa, matrix_aux2)
527  CALL mat_arxra(matrix_aux2, matrix_apvvpa, n, aa, rr, matrix_aux)
528 
529  ! Calculate w1o1:
530 
531  CALL parallel_gemm("N", "N", n, n, n, -1.0_dp, matrix_apvvpa, matrix_ava, 0.0_dp, matrix_aux2)
532  CALL mat_muld(matrix_aux2, matrix_apvvpa, matrix_apvpa, n, 1.0_dp, 1.0_dp, tt, rr, matrix_aux)
533  CALL mat_mulm(matrix_aux2, matrix_avva, matrix_ava, n, 1.0_dp, 1.0_dp, tt, rr, matrix_aux)
534  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  CALL parallel_gemm("N", "N", n, n, n, 1.0_dp, matrix_apvpa, matrix_avva, 0.0_dp, matrix_ev2)
539  CALL mat_muld(matrix_ev2, matrix_apvpa, matrix_apvvpa, n, -1.0_dp, 1.0_dp, tt, rr, matrix_aux)
540  CALL mat_mulm(matrix_ev2, matrix_ava, matrix_avva, n, -1.0_dp, 1.0_dp, tt, rr, matrix_aux)
541  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  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  CALL cp_fm_release(matrix_aux2)
554  CALL cp_fm_release(matrix_ava)
555  CALL cp_fm_release(matrix_avva)
556  CALL cp_fm_release(matrix_apvpa)
557  CALL cp_fm_release(matrix_apvvpa)
558 
559  CALL cp_fm_struct_release(matrix_full)
560 
561 ! WRITE (*,*) "CAW: DKH2 with even2c (Alex)"
562 ! WRITE (*,*) "JT: Now available in cp2k"
563 
564  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  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  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  ncol_global=n)
639 
640  CALL cp_fm_create(matrix_aux2, matrix_full)
641  CALL cp_fm_create(matrix_w1w1, matrix_full)
642  CALL cp_fm_create(matrix_w1e1w1, matrix_full)
643  CALL cp_fm_create(matrix_avva, matrix_full)
644  CALL cp_fm_create(matrix_apvvpa, matrix_full)
645 
646  CALL cp_fm_to_fm(matrix_v, matrix_avva)
647  CALL cp_fm_to_fm(matrix_pvp, matrix_apvvpa)
648 
649  ! Calculate vh = A V~ A:
650 
651  CALL mat_1_over_h(matrix_v, matrix_avva, e, matrix_aux)
652  CALL cp_fm_to_fm(matrix_avva, matrix_aux2)
653  CALL mat_axa(matrix_aux2, matrix_avva, n, aa, matrix_aux)
654 
655  ! Calculate pvph = A P V~ P A:
656 
657  CALL mat_1_over_h(matrix_pvp, matrix_apvvpa, e, matrix_aux)
658  CALL cp_fm_to_fm(matrix_apvvpa, matrix_aux2)
659  CALL mat_arxra(matrix_aux2, matrix_apvvpa, n, aa, rr, matrix_aux)
660 
661  ! Calculate w1w1:
662 
663  CALL parallel_gemm("N", "N", n, n, n, 1.0_dp, matrix_apvvpa, matrix_avva, 0.0_dp, matrix_w1w1)
664  CALL mat_muld(matrix_w1w1, matrix_apvvpa, matrix_apvvpa, n, -1.0_dp, 1.0_dp, tt, rr, matrix_aux2)
665  CALL mat_mulm(matrix_w1w1, matrix_avva, matrix_avva, n, -1.0_dp, 1.0_dp, tt, rr, matrix_aux2)
666  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  CALL mat_muld(matrix_aux, matrix_apvvpa, matrix_pe1p, n, 1.0_dp, 0.0_dp, tt, rr, matrix_aux2)
671  CALL parallel_gemm("N", "N", n, n, n, 1.0_dp, matrix_avva, matrix_pe1p, 0.0_dp, matrix_aux2)
672  CALL parallel_gemm("N", "N", n, n, n, 1.0_dp, matrix_aux, matrix_avva, 0.0_dp, matrix_w1e1w1)
673  CALL mat_muld(matrix_w1e1w1, matrix_aux, matrix_apvvpa, n, -1.0_dp, 1.0_dp, tt, rr, matrix_ev3)
674  CALL parallel_gemm("N", "N", n, n, n, -1.0_dp, matrix_aux2, matrix_avva, 1.0_dp, matrix_w1e1w1)
675  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  CALL parallel_gemm("N", "N", n, n, n, 0.5_dp, matrix_w1w1, matrix_ev1, 0.0_dp, matrix_ev3)
682  CALL parallel_gemm("N", "N", n, n, n, 0.5_dp, matrix_ev1, matrix_w1w1, 1.0_dp, matrix_ev3)
683  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  CALL cp_fm_release(matrix_aux2)
690  CALL cp_fm_release(matrix_avva)
691  CALL cp_fm_release(matrix_apvvpa)
692  CALL cp_fm_release(matrix_w1w1)
693  CALL cp_fm_release(matrix_w1e1w1)
694 
695  CALL cp_fm_struct_release(matrix_full)
696 
697 ! WRITE (*,*) "CAW: DKH3 with even3b (Alex)"
698 ! WRITE (*,*) "JT: Now available in cp2k"
699 
700  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  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  REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: o1w1, pvp, pvph, scr_1, scr_2, scr_3, &
767  scr_4, scrh_1, scrh_2, scrh_3, scrh_4, &
768  sum_1, sum_2, v, vh, w1o1, w1w1
769 
770 !C-----------------------------------------------------------------------
771 !C 1. General Structures and Patterns for DKH4
772 !C-----------------------------------------------------------------------
773 
774  ALLOCATE (v(n, n))
775  ALLOCATE (pvp(n, n))
776  ALLOCATE (vh(n, n))
777  ALLOCATE (pvph(n, n))
778  v = 0.0_dp
779  pvp = 0.0_dp
780  vh = 0.0_dp
781  pvph = 0.0_dp
782  v(1:n, 1:n) = vv(1:n, 1:n)
783  vh(1:n, 1:n) = vv(1:n, 1:n)
784  pvp(1:n, 1:n) = gg(1:n, 1:n)
785  pvph(1:n, 1:n) = gg(1:n, 1:n)
786 
787  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  ALLOCATE (w1w1(n, n))
808  w1w1 = 0.0_dp
809  ALLOCATE (w1o1(n, n))
810  w1o1 = 0.0_dp
811  ALLOCATE (o1w1(n, n))
812  o1w1 = 0.0_dp
813  ALLOCATE (sum_1(n, n))
814  sum_1 = 0.0_dp
815  ALLOCATE (sum_2(n, n))
816  sum_2 = 0.0_dp
817  ALLOCATE (scr_1(n, n))
818  scr_1 = 0.0_dp
819  ALLOCATE (scr_2(n, n))
820  scr_2 = 0.0_dp
821  ALLOCATE (scr_3(n, n))
822  scr_3 = 0.0_dp
823  ALLOCATE (scr_4(n, n))
824  scr_4 = 0.0_dp
825  ALLOCATE (scrh_1(n, n))
826  scrh_1 = 0.0_dp
827  ALLOCATE (scrh_2(n, n))
828  scrh_2 = 0.0_dp
829  ALLOCATE (scrh_3(n, n))
830  scrh_3 = 0.0_dp
831  ALLOCATE (scrh_4(n, n))
832  scrh_4 = 0.0_dp
833 
834  ! Calculate w1w1:
835  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  CALL dgemm("N", "N", n, n, n, 1.0_dp, vh, n, pvph, n, 1.0_dp, w1w1, n)
839 
840  ! Calculate w1o1:
841  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  CALL dgemm("N", "N", n, n, n, -1.0_dp, vh, n, pvp, n, 1.0_dp, w1o1, n)
845  ! Calculate o1w1:
846  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  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  CALL dgemm("N", "N", n, n, n, 1.0_dp, vh, n, e1, n, 0.0_dp, scr_1, n)
858  CALL dgemm("N", "N", n, n, n, 1.0_dp, pvph, n, e1, n, 0.0_dp, scr_2, n)
859  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  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  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  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  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  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  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  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  CALL dgemm("N", "N", n, n, n, -0.5_dp, scrh_4, n, scr_3, n, 1.0_dp, sum_1, n)
891  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  CALL dgemm("N", "N", n, n, n, 1.0_dp, vh, n, pe1p, n, 0.0_dp, scr_2, n)
897  CALL dgemm("N", "N", n, n, n, 1.0_dp, e1, n, pvph, n, 0.0_dp, scr_3, n)
898  CALL dgemm("N", "N", n, n, n, 1.0_dp, e1, n, vh, n, 0.0_dp, scr_4, n)
899 
900  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  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  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  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  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  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  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  CALL dgemm("N", "N", n, n, n, -0.5_dp, scr_4, n, scrh_3, n, 1.0_dp, sum_1, n)
929  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  CALL dgemm("N", "N", n, n, n, 0.125_dp, w1w1, n, w1o1, n, 0.0_dp, sum_2, n)
938  CALL dgemm("N", "N", n, n, n, -0.375_dp, w1w1, n, o1w1, n, 1.0_dp, sum_2, n)
939  CALL dgemm("N", "N", n, n, n, 0.375_dp, w1o1, n, w1w1, n, 1.0_dp, sum_2, n)
940  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  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  DEALLOCATE (v, pvp, vh, pvph, w1w1, w1o1, o1w1, sum_1, sum_2)
953  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  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  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  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  row_indices=row_indices, col_indices=col_indices)
1003 
1004  DO j = 1, ncol_local
1005  DO i = 1, nrow_local
1006  matrix_aux%local_data(i, j) = 1/(e(row_indices(i)) + e(col_indices(j)))
1007  END DO
1008  END DO
1009 
1010  CALL cp_fm_schur_product(matrix_p, matrix_aux, matrix_pp)
1011 
1012  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  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  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  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  ncol_global=1)
1050 
1051  CALL cp_fm_create(vec_a, vec_full)
1052 
1053  CALL cp_fm_get_info(matrix_x, nrow_local=nrow_local, &
1054  row_indices=row_indices)
1055 
1056  DO i = 1, nrow_local
1057  vec_a%local_data(i, 1) = a(row_indices(i))
1058  END DO
1059 
1060  CALL cp_fm_syrk('U', 'N', 1, 1.0_dp, vec_a, 1, 1, 0.0_dp, matrix_aux)
1061  CALL cp_fm_upper_to_full(matrix_aux, matrix_axa)
1062  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  CALL cp_fm_release(vec_a)
1071  CALL cp_fm_struct_release(vec_full)
1072 
1073  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  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  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  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  ncol_global=1)
1113 
1114  CALL cp_fm_get_info(matrix_aux, nrow_local=nrow_local, &
1115  row_indices=row_indices)
1116 
1117  CALL cp_fm_create(vec_a, vec_full)
1118 
1119  DO i = 1, nrow_local
1120  vec_a%local_data(i, 1) = a(row_indices(i))*r(row_indices(i))
1121  END DO
1122 
1123  CALL cp_fm_syrk('U', 'N', 1, 1.0_dp, vec_a, 1, 1, 0.0_dp, matrix_aux)
1124  CALL cp_fm_upper_to_full(matrix_aux, matrix_axa)
1125  CALL cp_fm_schur_product(matrix_x, matrix_aux, matrix_axa)
1126 
1127  CALL cp_fm_release(vec_a)
1128  CALL cp_fm_struct_release(vec_full)
1129 
1130  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  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  REAL(KIND=dp), DIMENSION(n) :: vec
1171 
1172  CALL cp_fm_to_fm(matrix_q, matrix_aux)
1173 
1174  DO i = 1, n
1175  vec(i) = 2.0_dp*t(i)*rr(i)*rr(i)
1176  END DO
1177  CALL cp_fm_column_scale(matrix_aux, vec)
1178 
1179  CALL parallel_gemm("N", "N", n, n, n, alpha, matrix_aux, matrix_r, beta, matrix_p)
1180 
1181  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  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  REAL(KIND=dp), DIMENSION(n) :: vec
1222 
1223  CALL cp_fm_to_fm(matrix_q, matrix_aux)
1224 
1225  DO i = 1, n
1226  vec(i) = 0.5_dp/(t(i)*rr(i)*rr(i))
1227  END DO
1228 
1229  CALL cp_fm_column_scale(matrix_aux, vec)
1230 
1231  CALL parallel_gemm("N", "N", n, n, n, alpha, matrix_aux, matrix_r, beta, matrix_p)
1232 
1233  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  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  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: aa, e, ev0t, rr, tt
1287  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: aux, eig, ev1, ev1t, ev2, ev2t, ev3, &
1288  ev3t, ev4, ev4t, ove, pev1tp, pvpt, &
1289  revt, sinv, vt
1290 
1291  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  ALLOCATE (eig(n, n))
1300  ALLOCATE (sinv(n, n))
1301  ALLOCATE (revt(n, n))
1302  ALLOCATE (aux(n, n))
1303  ALLOCATE (ove(n, n))
1304  ALLOCATE (ev0t(n))
1305  ALLOCATE (e(n))
1306  ALLOCATE (aa(n))
1307  ALLOCATE (rr(n))
1308  ALLOCATE (tt(n))
1309  ALLOCATE (ev1t(n, n))
1310  ALLOCATE (ev2t(n, n))
1311  ALLOCATE (ev3t(n, n))
1312  ALLOCATE (ev4t(n, n))
1313  ALLOCATE (vt(n, n))
1314  ALLOCATE (pvpt(n, n))
1315  ALLOCATE (pev1tp(n, n))
1316  ALLOCATE (ev1(n, n))
1317  ALLOCATE (ev2(n, n))
1318  ALLOCATE (ev3(n, n))
1319  ALLOCATE (ev4(n, n))
1320 
1321  !-----------------------------------------------------------------------
1322  ! Schmidt-orthogonalize overlap matrix
1323  !-----------------------------------------------------------------------
1324 
1325  CALL sog(n, s, sinv)
1326 
1327  !-----------------------------------------------------------------------
1328  ! Calculate matrix representation from nonrelativistic T matrix
1329  !-----------------------------------------------------------------------
1330 
1331  CALL dkh_diag(h, n, eig, tt, sinv, aux, 0)
1332 
1333  !-----------------------------------------------------------------------
1334  ! Calculate kinetic part of Hamiltonian in T-basis
1335  !-----------------------------------------------------------------------
1336 
1337  CALL kintegral_a(n, ev0t, tt, e)
1338 
1339  !-----------------------------------------------------------------------
1340  ! Calculate reverse transformation matrix revt
1341  !-----------------------------------------------------------------------
1342 
1343  CALL dgemm("N", "N", n, n, n, 1.0_dp, sinv, n, eig, n, 0.0_dp, aux, n)
1344  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  h = 0.0_dp
1351  DO i = 1, n
1352  DO j = 1, i
1353  DO k = 1, n
1354  h(i, j) = h(i, j) + revt(i, k)*revt(j, k)*ev0t(k)
1355  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  DO i = 1, n
1365  aa(i) = sqrt((c_light_au**2 + e(i))/(2.0_dp*e(i)))
1366  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  CALL trsm(v, sinv, ove, n, aux)
1374  CALL trsm(ove, eig, vt, n, aux)
1375 
1376  !-----------------------------------------------------------------------
1377  ! Transform pVp integrals to T-basis (pVp -> pVpt)
1378  !-----------------------------------------------------------------------
1379 
1380  CALL trsm(pvp, sinv, ove, n, aux)
1381  CALL trsm(ove, eig, pvpt, n, aux)
1382 
1383  !-----------------------------------------------------------------------
1384  ! Calculate even1 in T-basis
1385  !-----------------------------------------------------------------------
1386 
1387  IF (dkh_order .GE. 1) THEN
1388  CALL even1_a(n, ev1t, vt, pvpt, aa, rr)
1389 
1390  !----------------------------------------------------------------------
1391  ! Transform even1 back to position space
1392  !----------------------------------------------------------------------
1393 
1394  CALL dgemm("N", "N", n, n, n, 1.0_dp, revt, n, ev1t, n, 0.0_dp, aux, n)
1395  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  IF (dkh_order .GE. 2) THEN
1403  CALL even2c_a(n, ev2t, vt, pvpt, aa, rr, tt, e)
1404 
1405  !-----------------------------------------------------------------------
1406  ! Transform even2 back to position space
1407  !-----------------------------------------------------------------------
1408 
1409  aux = 0.0_dp
1410  CALL dgemm("N", "N", n, n, n, 1.0_dp, revt, n, ev2t, n, 0.0_dp, aux, n)
1411  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  IF (dkh_order .GE. 3) THEN
1419  CALL peven1p_a(n, pev1tp, vt, pvpt, aa, rr, tt)
1420  CALL even3b_a(n, ev3t, ev1t, pev1tp, vt, pvpt, aa, rr, tt, e)
1421 
1422  !-----------------------------------------------------------------------
1423  ! Transform even3 back to position space
1424  !-----------------------------------------------------------------------
1425  aux = 0.0_dp
1426  CALL dgemm("N", "N", n, n, n, 1.0_dp, revt, n, ev3t, n, 0.0_dp, aux, n)
1427  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  IF (dkh_order .GE. 4) THEN
1434  CALL even4a_a(n, ev4t, ev1t, pev1tp, vt, pvpt, aa, rr, tt, e)
1435 
1436  !-----------------------------------------------------------------------
1437  ! Transform even4 back to position space
1438  !-----------------------------------------------------------------------
1439  aux = 0.0_dp
1440  CALL dgemm("N", "N", n, n, n, 1.0_dp, revt, n, ev4t, n, 0.0_dp, aux, n)
1441  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  IF (dkh_order .GE. 4) THEN
1446  cpabort("DKH 4")
1447  END IF
1448  !-----------------------------------------------------------------------
1449  ! Calculate v in position space
1450  !-----------------------------------------------------------------------
1451 
1452  IF (dkh_order .GE. 1) THEN
1453  CALL mat_add2(v, 0.0_dp, 1.0_dp, ev1, n)
1454  END IF
1455  IF (dkh_order .GE. 2) THEN
1456  CALL mat_add2(v, 1.0_dp, 1.0_dp, ev2, n)
1457  END IF
1458  IF (dkh_order .GE. 3) THEN
1459  CALL mat_add2(v, 1.0_dp, 1.0_dp, ev3, n)
1460  END IF
1461  IF (dkh_order .GE. 4) THEN
1462  CALL mat_add2(v, 1.0_dp, 1.0_dp, ev4, n)
1463  END IF
1464 
1465  !-----------------------------------------------------------------------
1466 
1467  DEALLOCATE (eig, sinv, revt, ove, aux, vt, pvpt, ev1, ev2, ev3, ev4, ev1t, ev2t, ev3t, ev4t, pev1tp)
1468  DEALLOCATE (ev0t, e, aa, rr, tt)
1469 
1470  END SUBROUTINE dkh_atom_transformation
1471 
1472 ! **************************************************************************************************
1473 !> \brief ...
1474 !> \param n ...
1475 !> \param ev0t ...
1476 !> \param tt ...
1477 !> \param e ...
1478 ! **************************************************************************************************
1479  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  DO i = 1, n
1491  IF (tt(i) .LT. 0.0_dp) THEN
1492  WRITE (*, *) ' dkh_main.F | tt(', i, ') = ', tt(i)
1493  END IF
1494 
1495  ! Calculate some constants
1496 
1497  prea = 1/(c_light_au**2)
1498  con2 = prea + prea
1499  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  ev0t(i) = tt(i)
1505  ratio = tt(i)/c_light_au
1506  IF (ratio .LE. 0.02_dp) THEN
1507  tv1 = tt(i)
1508  tv2 = -tv1*tt(i)*prea/2.0_dp
1509  tv3 = -tv2*tt(i)*prea
1510  tv4 = -tv3*tt(i)*prea*1.25_dp
1511  ev0t(i) = tv1 + tv2 + tv3 + tv4
1512  ELSE
1513  ev0t(i) = con*(sqrt(1.0_dp + con2*tt(i)) - 1.0_dp)
1514  END IF
1515  e(i) = ev0t(i) + con
1516  END DO
1517 
1518  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  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  DO i = 1, n
1554  DO j = 1, i
1555  ev1t(i, j) = vt(i, j)*aa(i)*aa(j) + pvpt(i, j)*aa(i)*rr(i)*aa(j)*rr(j)
1556  ev1t(j, i) = ev1t(i, j)
1557  END DO
1558  END DO
1559 
1560  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  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  DO i = 1, n
1597  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  pvpt(i, j)*aa(i)*rr(i)*aa(j)*rr(j)
1600  pev1tp(j, i) = pev1tp(i, j)
1601  END DO
1602  END DO
1603 
1604  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  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  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: o1w1, pvp, pvph, v, vh, w1o1
1659 
1660 !-----------------------------------------------------------------------
1661 ! 1. General Structures and Patterns for DKH2
1662 !-----------------------------------------------------------------------
1663 
1664  ALLOCATE (v(n, n))
1665  ALLOCATE (pvp(n, n))
1666  ALLOCATE (vh(n, n))
1667  ALLOCATE (pvph(n, n))
1668  v = 0.0_dp
1669  pvp = 0.0_dp
1670  vh = 0.0_dp
1671  pvph = 0.0_dp
1672  v(1:n, 1:n) = vv(1:n, 1:n)
1673  vh(1:n, 1:n) = vv(1:n, 1:n)
1674  pvp(1:n, 1:n) = gg(1:n, 1:n)
1675  pvph(1:n, 1:n) = gg(1:n, 1:n)
1676 
1677  ev2 = 0.0_dp
1678  ! Calculate v = A V A:
1679 
1680  CALL mat_axa_a(v, n, aa)
1681 
1682  ! Calculate pvp = A P V P A:
1683 
1684  CALL mat_arxra_a(pvp, n, aa, rr)
1685 
1686  ! Calculate vh = A V~ A:
1687 
1688  CALL mat_1_over_h_a(vh, n, e)
1689  CALL mat_axa_a(vh, n, aa)
1690 
1691  ! Calculate pvph = A P V~ P A:
1692 
1693  CALL mat_1_over_h_a(pvph, n, e)
1694  CALL mat_arxra_a(pvph, n, aa, rr)
1695 
1696  ! Create/Initialize necessary matrices:
1697  ALLOCATE (w1o1(n, n))
1698  ALLOCATE (o1w1(n, n))
1699  w1o1 = 0.0_dp
1700  o1w1 = 0.0_dp
1701 
1702  ! Calculate w1o1:
1703  CALL dgemm("N", "N", n, n, n, -1.0_dp, pvph, n, v, n, 0.0_dp, w1o1, n)
1704  CALL mat_muld_a(w1o1, pvph, pvp, n, 1.0_dp, 1.0_dp, tt, rr)
1705  CALL mat_mulm_a(w1o1, vh, v, n, 1.0_dp, 1.0_dp, tt, rr)
1706  CALL dgemm("N", "N", n, n, n, -1.0_dp, vh, n, pvp, n, 1.0_dp, w1o1, n)
1707  ! Calculate o1w1:
1708  CALL dgemm("N", "N", n, n, n, 1.0_dp, pvp, n, vh, n, 0.0_dp, o1w1, n)
1709  CALL mat_muld_a(o1w1, pvp, pvph, n, -1.0_dp, 1.0_dp, tt, rr)
1710  CALL mat_mulm_a(o1w1, v, vh, n, -1.0_dp, 1.0_dp, tt, rr)
1711  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  CALL mat_add(ev2, 0.5_dp, w1o1, -0.5_dp, o1w1, n)
1719 
1720  !-----------------------------------------------------------------------
1721  ! 3. Finish up the stuff!!
1722  !-----------------------------------------------------------------------
1723 
1724  DEALLOCATE (v, vh, pvp, pvph, w1o1, o1w1)
1725 
1726 ! WRITE (*,*) "CAW: DKH2 with even2c (Alex)"
1727 ! WRITE (*,*) "!JT: Now available in cp2k"
1728 
1729  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  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  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  ALLOCATE (vh(n, n))
1795  ALLOCATE (pvph(n, n))
1796  vh = 0.0_dp
1797  pvph = 0.0_dp
1798  vh(1:n, 1:n) = vv(1:n, 1:n)
1799  pvph(1:n, 1:n) = gg(1:n, 1:n)
1800 
1801  ev3 = 0.0_dp
1802 
1803  ! Calculate vh = A V~ A:
1804 
1805  CALL mat_1_over_h_a(vh, n, e)
1806  CALL mat_axa_a(vh, n, aa)
1807 
1808  ! Calculate pvph = A P V~ P A:
1809 
1810  CALL mat_1_over_h_a(pvph, n, e)
1811  CALL mat_arxra_a(pvph, n, aa, rr)
1812 
1813  ! Create/Initialize necessary matrices:
1814  ALLOCATE (w1w1(n, n))
1815  ALLOCATE (w1e1w1(n, n))
1816  ALLOCATE (scr_1(n, n))
1817  ALLOCATE (scr_2(n, n))
1818  w1w1 = 0.0_dp
1819  w1e1w1 = 0.0_dp
1820  scr_1 = 0.0_dp
1821  scr_2 = 0.0_dp
1822 
1823  ! Calculate w1w1:
1824  CALL dgemm("N", "N", n, n, n, 1.0_dp, pvph, n, vh, n, 0.0_dp, w1w1, n)
1825  CALL mat_muld_a(w1w1, pvph, pvph, n, -1.0_dp, 1.0_dp, tt, rr)
1826  CALL mat_mulm_a(w1w1, vh, vh, n, -1.0_dp, 1.0_dp, tt, rr)
1827  CALL dgemm("N", "N", n, n, n, 1.0_dp, vh, n, pvph, n, 1.0_dp, w1w1, n)
1828 
1829  ! Calculate w1e1w1:
1830  CALL mat_muld_a(scr_1, pvph, pe1p, n, 1.0_dp, 0.0_dp, tt, rr)
1831  CALL dgemm("N", "N", n, n, n, 1.0_dp, vh, n, pe1p, n, 0.0_dp, scr_2, n)
1832  CALL dgemm("N", "N", n, n, n, 1.0_dp, scr_1, n, vh, n, 0.0_dp, w1e1w1, n)
1833  CALL mat_muld_a(w1e1w1, scr_1, pvph, n, -1.0_dp, 1.0_dp, tt, rr)
1834  CALL dgemm("N", "N", n, n, n, -1.0_dp, scr_2, n, vh, n, 1.0_dp, w1e1w1, n)
1835  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  CALL dgemm("N", "N", n, n, n, 0.5_dp, w1w1, n, e1, n, 0.0_dp, ev3, n)
1842  CALL dgemm("N", "N", n, n, n, 0.5_dp, e1, n, w1w1, n, 1.0_dp, ev3, n)
1843  CALL mat_add2(ev3, 1.0_dp, -1.0_dp, w1e1w1, n)
1844 
1845  !-----------------------------------------------------------------------
1846  ! 3. Finish up the stuff!!
1847  !-----------------------------------------------------------------------
1848 
1849  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  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  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  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: o1w1, pvp, pvph, scr_1, scr_2, scr_3, &
1923  scr_4, scrh_1, scrh_2, scrh_3, scrh_4, &
1924  sum_1, sum_2, v, vh, w1o1, w1w1
1925 
1926 !C-----------------------------------------------------------------------
1927 !C 1. General Structures and Patterns for DKH4
1928 !C-----------------------------------------------------------------------
1929 
1930  ALLOCATE (v(n, n))
1931  ALLOCATE (pvp(n, n))
1932  ALLOCATE (vh(n, n))
1933  ALLOCATE (pvph(n, n))
1934  v = 0.0_dp
1935  pvp = 0.0_dp
1936  vh = 0.0_dp
1937  pvph = 0.0_dp
1938  v(1:n, 1:n) = vv(1:n, 1:n)
1939  vh(1:n, 1:n) = vv(1:n, 1:n)
1940  pvp(1:n, 1:n) = gg(1:n, 1:n)
1941  pvph(1:n, 1:n) = gg(1:n, 1:n)
1942 
1943  ev4 = 0.0_dp
1944  ! Calculate v = A V A:
1945 
1946  CALL mat_axa_a(v, n, aa)
1947 
1948  ! Calculate pvp = A P V P A:
1949 
1950  CALL mat_arxra_a(pvp, n, aa, rr)
1951 
1952  ! Calculate vh = A V~ A:
1953 
1954  CALL mat_1_over_h_a(vh, n, e)
1955  CALL mat_axa_a(vh, n, aa)
1956 
1957  ! Calculate pvph = A P V~ P A:
1958 
1959  CALL mat_1_over_h_a(pvph, n, e)
1960  CALL mat_arxra_a(pvph, n, aa, rr)
1961 
1962  ! Create/Initialize necessary matrices:
1963  ALLOCATE (w1w1(n, n))
1964  w1w1 = 0.0_dp
1965  ALLOCATE (w1o1(n, n))
1966  w1o1 = 0.0_dp
1967  ALLOCATE (o1w1(n, n))
1968  o1w1 = 0.0_dp
1969  ALLOCATE (sum_1(n, n))
1970  sum_1 = 0.0_dp
1971  ALLOCATE (sum_2(n, n))
1972  sum_2 = 0.0_dp
1973  ALLOCATE (scr_1(n, n))
1974  scr_1 = 0.0_dp
1975  ALLOCATE (scr_2(n, n))
1976  scr_2 = 0.0_dp
1977  ALLOCATE (scr_3(n, n))
1978  scr_3 = 0.0_dp
1979  ALLOCATE (scr_4(n, n))
1980  scr_4 = 0.0_dp
1981  ALLOCATE (scrh_1(n, n))
1982  scrh_1 = 0.0_dp
1983  ALLOCATE (scrh_2(n, n))
1984  scrh_2 = 0.0_dp
1985  ALLOCATE (scrh_3(n, n))
1986  scrh_3 = 0.0_dp
1987  ALLOCATE (scrh_4(n, n))
1988  scrh_4 = 0.0_dp
1989 
1990  ! Calculate w1w1:
1991  CALL dgemm("N", "N", n, n, n, 1.0_dp, pvph, n, vh, n, 0.0_dp, w1w1, n)
1992  CALL mat_muld_a(w1w1, pvph, pvph, n, -1.0_dp, 1.0_dp, tt, rr)
1993  CALL mat_mulm_a(w1w1, vh, vh, n, -1.0_dp, 1.0_dp, tt, rr)
1994  CALL dgemm("N", "N", n, n, n, 1.0_dp, vh, n, pvph, n, 1.0_dp, w1w1, n)
1995 
1996  ! Calculate w1o1:
1997  CALL dgemm("N", "N", n, n, n, -1.0_dp, pvph, n, v, n, 0.0_dp, w1o1, n)
1998  CALL mat_muld_a(w1o1, pvph, pvp, n, 1.0_dp, 1.0_dp, tt, rr)
1999  CALL mat_mulm_a(w1o1, vh, v, n, 1.0_dp, 1.0_dp, tt, rr)
2000  CALL dgemm("N", "N", n, n, n, -1.0_dp, vh, n, pvp, n, 1.0_dp, w1o1, n)
2001  ! Calculate o1w1:
2002  CALL dgemm("N", "N", n, n, n, 1.0_dp, pvp, n, vh, n, 0.0_dp, o1w1, n)
2003  CALL mat_muld_a(o1w1, pvp, pvph, n, -1.0_dp, 1.0_dp, tt, rr)
2004  CALL mat_mulm_a(o1w1, v, vh, n, -1.0_dp, 1.0_dp, tt, rr)
2005  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  CALL dgemm("N", "N", n, n, n, 1.0_dp, vh, n, e1, n, 0.0_dp, scr_1, n)
2014  CALL dgemm("N", "N", n, n, n, 1.0_dp, pvph, n, e1, n, 0.0_dp, scr_2, n)
2015  CALL dgemm("N", "N", n, n, n, 1.0_dp, pe1p, n, vh, n, 0.0_dp, scr_3, n)
2016  CALL mat_muld_a(scr_4, pe1p, pvph, n, 1.0_dp, 0.0_dp, tt, rr)
2017 
2018  CALL mat_muld_a(scrh_1, pvph, pe1p, n, 1.0_dp, 0.0_dp, tt, rr)
2019  CALL mat_1_over_h_a(scrh_1, n, e)
2020  CALL dgemm("N", "N", n, n, n, 1.0_dp, vh, n, pe1p, n, 0.0_dp, scrh_2, n)
2021  CALL mat_1_over_h_a(scrh_2, n, e)
2022  CALL dgemm("N", "N", n, n, n, 1.0_dp, e1, n, pvph, n, 0.0_dp, scrh_3, n)
2023  CALL mat_1_over_h_a(scrh_3, n, e)
2024  CALL dgemm("N", "N", n, n, n, 1.0_dp, e1, n, vh, n, 0.0_dp, scrh_4, n)
2025  CALL mat_1_over_h_a(scrh_4, n, e)
2026 
2027  ! 2a) sum_1 = 1/2 W2W1E1 ( [1]-[8] )
2028 
2029  CALL dgemm("N", "N", n, n, n, 0.5_dp, scrh_1, n, scr_1, n, 0.0_dp, sum_1, n)
2030  CALL mat_muld_a(sum_1, scrh_1, scr_2, n, -0.5_dp, 1.0_dp, tt, rr)
2031  CALL dgemm("N", "N", n, n, n, -0.5_dp, scrh_2, n, scr_1, n, 1.0_dp, sum_1, n)
2032  CALL mat_muld_a(sum_1, scrh_2, scr_2, n, 0.5_dp, 1.0_dp, tt, rr)
2033  CALL dgemm("N", "N", n, n, n, -0.5_dp, scrh_3, n, scr_1, n, 1.0_dp, sum_1, n)
2034  CALL mat_muld_a(sum_1, scrh_3, scr_2, n, 0.5_dp, 1.0_dp, tt, rr)
2035  CALL mat_mulm_a(sum_1, scrh_4, scr_1, n, 0.5_dp, 1.0_dp, tt, rr)
2036  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  CALL mat_muld_a(sum_1, scrh_1, scr_3, n, -0.5_dp, 1.0_dp, tt, rr)
2041  CALL mat_muld_a(sum_1, scrh_1, scr_4, n, 0.5_dp, 1.0_dp, tt, rr)
2042  CALL mat_muld_a(sum_1, scrh_2, scr_3, n, 0.5_dp, 1.0_dp, tt, rr)
2043  CALL mat_muld_a(sum_1, scrh_2, scr_4, n, -0.5_dp, 1.0_dp, tt, rr)
2044  CALL mat_muld_a(sum_1, scrh_3, scr_3, n, 0.5_dp, 1.0_dp, tt, rr)
2045  CALL mat_muld_a(sum_1, scrh_3, scr_4, n, -0.5_dp, 1.0_dp, tt, rr)
2046  CALL dgemm("N", "N", n, n, n, -0.5_dp, scrh_4, n, scr_3, n, 1.0_dp, sum_1, n)
2047  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  CALL mat_muld_a(scr_1, pvph, pe1p, n, 1.0_dp, 0.0_dp, tt, rr)
2052  CALL dgemm("N", "N", n, n, n, 1.0_dp, vh, n, pe1p, n, 0.0_dp, scr_2, n)
2053  CALL dgemm("N", "N", n, n, n, 1.0_dp, e1, n, pvph, n, 0.0_dp, scr_3, n)
2054  CALL dgemm("N", "N", n, n, n, 1.0_dp, e1, n, vh, n, 0.0_dp, scr_4, n)
2055 
2056  CALL dgemm("N", "N", n, n, n, 1.0_dp, vh, n, e1, n, 0.0_dp, scrh_1, n)
2057  CALL mat_1_over_h_a(scrh_1, n, e)
2058  CALL dgemm("N", "N", n, n, n, 1.0_dp, pvph, n, e1, n, 0.0_dp, scrh_2, n)
2059  CALL mat_1_over_h_a(scrh_2, n, e)
2060  CALL dgemm("N", "N", n, n, n, 1.0_dp, pe1p, n, vh, n, 0.0_dp, scr_3, n)
2061  CALL mat_1_over_h_a(scrh_3, n, e)
2062  CALL mat_muld_a(scrh_4, pe1p, pvph, n, 1.0_dp, 0.0_dp, tt, rr)
2063  CALL mat_1_over_h_a(scrh_4, n, e)
2064 
2065  ! 2c) sum_1 = - 1/2 W1E1W2 (+ sum_1) ( [17]-[24] )
2066 
2067  CALL dgemm("N", "N", n, n, n, 0.5_dp, scr_1, n, scrh_1, n, 0.0_dp, sum_1, n)
2068  CALL mat_muld_a(sum_1, scr_1, scrh_2, n, -0.5_dp, 1.0_dp, tt, rr)
2069  CALL dgemm("N", "N", n, n, n, -0.5_dp, scr_2, n, scrh_1, n, 1.0_dp, sum_1, n)
2070  CALL mat_muld_a(sum_1, scr_2, scrh_2, n, 0.5_dp, 1.0_dp, tt, rr)
2071  CALL mat_muld_a(sum_1, scr_1, scrh_3, n, -0.5_dp, 1.0_dp, tt, rr)
2072  CALL mat_muld_a(sum_1, scr_1, scrh_4, n, 0.5_dp, 1.0_dp, tt, rr)
2073  CALL mat_muld_a(sum_1, scr_2, scrh_3, n, 0.5_dp, 1.0_dp, tt, rr)
2074  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  CALL dgemm("N", "N", n, n, n, -0.5_dp, scr_3, n, scrh_1, n, 0.0_dp, sum_1, n)
2079  CALL mat_muld_a(sum_1, scr_3, scrh_2, n, 0.5_dp, 1.0_dp, tt, rr)
2080  CALL mat_mulm_a(sum_1, scr_4, scrh_1, n, 0.5_dp, 1.0_dp, tt, rr)
2081  CALL dgemm("N", "N", n, n, n, -0.5_dp, scr_4, n, scrh_2, n, 1.0_dp, sum_1, n)
2082  CALL mat_muld_a(sum_1, scr_3, scrh_3, n, 0.5_dp, 1.0_dp, tt, rr)
2083  CALL mat_muld_a(sum_1, scr_3, scrh_4, n, -0.5_dp, 1.0_dp, tt, rr)
2084  CALL dgemm("N", "N", n, n, n, -0.5_dp, scr_4, n, scrh_3, n, 1.0_dp, sum_1, n)
2085  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  CALL dgemm("N", "N", n, n, n, 0.125_dp, w1w1, n, w1o1, n, 0.0_dp, sum_2, n)
2094  CALL dgemm("N", "N", n, n, n, -0.375_dp, w1w1, n, o1w1, n, 1.0_dp, sum_2, n)
2095  CALL dgemm("N", "N", n, n, n, 0.375_dp, w1o1, n, w1w1, n, 1.0_dp, sum_2, n)
2096  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  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  DEALLOCATE (v, pvp, vh, pvph, w1w1, w1o1, o1w1, sum_1, sum_2)
2109  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  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  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  DO i = 1, n
2157  DO j = 1, n
2158  p(i, j) = p(i, j)/(e(i) + e(j))
2159  END DO
2160  END DO
2161 
2162  END SUBROUTINE mat_1_over_h_a
2163 
2164 ! **************************************************************************************************
2165 !> \brief ...
2166 !> \param p ...
2167 !> \param n ...
2168 !> \param a ...
2169 ! **************************************************************************************************
2170  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  DO i = 1, n
2189  DO j = 1, n
2190  p(i, j) = p(i, j)*a(i)*a(j)
2191  END DO
2192  END DO
2193 
2194  END SUBROUTINE mat_axa_a
2195 
2196 ! **************************************************************************************************
2197 !> \brief ...
2198 !> \param p ...
2199 !> \param n ...
2200 !> \param a ...
2201 !> \param r ...
2202 ! **************************************************************************************************
2203  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  DO i = 1, n
2223  DO j = 1, n
2224  p(i, j) = p(i, j)*a(i)*a(j)*r(i)*r(j)
2225  END DO
2226  END DO
2227 
2228  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  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  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: qtemp
2267 
2268  ALLOCATE (qtemp(n, n))
2269 
2270  DO i = 1, n
2271  DO j = 1, n
2272  qtemp(i, j) = q(i, j)*2.0_dp*t(j)*rr(j)*rr(j)
2273  END DO
2274  END DO
2275 
2276  CALL dgemm("N", "N", n, n, n, alpha, qtemp, n, r, n, beta, p, n)
2277 
2278  DEALLOCATE (qtemp)
2279 
2280  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  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  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: qtemp
2319 
2320  ALLOCATE (qtemp(n, n))
2321 
2322  DO i = 1, n
2323  DO j = 1, n
2324  qtemp(i, j) = q(i, j)*0.5_dp/(t(j)*rr(j)*rr(j))
2325  END DO
2326  END DO
2327 
2328  CALL dgemm("N", "N", n, n, n, alpha, qtemp, n, r, n, beta, p, n)
2329 
2330  DEALLOCATE (qtemp)
2331 
2332  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  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  DO i = 1, n
2372  DO j = 1, n
2373  p(i, j) = alpha*p(i, j) + beta*r(i, j)
2374  END DO
2375  END DO
2376 
2377  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  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  DO i = 1, n
2420  DO j = 1, n
2421  p(i, j) = alpha*q(i, j) + beta*r(i, j)
2422  END DO
2423  END DO
2424 
2425  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  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  ij = 0
2450  DO i = 1, n
2451  DO j = 1, i
2452  ij = ij + 1
2453  c(i, j) = 0.0_dp
2454  c(j, i) = 0.0_dp
2455  h(i, j) = 0.0_dp
2456  h(j, i) = 0.0_dp
2457  END DO
2458  END DO
2459  DO i = 1, n
2460  DO l = 1, n
2461  DO k = 1, n
2462  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  DO i = 1, n
2469  DO j = 1, i
2470  ij = ij + 1
2471  DO l = 1, n
2472  c(i, j) = h(i, l)*b(l, j) + c(i, j)
2473  c(j, i) = c(i, j)
2474  END DO
2475  END DO
2476  END DO
2477 
2478  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  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  eig = 0.0_dp
2503  aux = 0.0_dp
2504 
2505  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  aux = 0.0_dp
2508 
2509  CALL dgemm("T", "N", n, n, n, 1.0_dp, matrix_sinv_pgf, n, eig, n, 0.0_dp, aux, n)
2510 
2511  n2 = 3*n - 1
2512 
2513  CALL jacob2(aux, eig, ew, n, ic)
2514 
2515  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  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  tol = 1.0e-15
2538  ext_norm = 0.0_dp
2539  u1 = real(n, kind=dp)
2540  DO i = 1, n
2541  eigv(i, i) = 1.0_dp
2542  eigw(i) = sogt(i, i)
2543  DO j = 1, i
2544  IF (i .NE. j) THEN
2545  eigv(i, j) = 0.0_dp
2546  eigv(j, i) = 0.0_dp
2547  ext_norm = ext_norm + sogt(i, j)*sogt(i, j)
2548  END IF
2549  END DO
2550  END DO
2551 
2552  IF (ext_norm .GT. 0.0_dp) THEN
2553  ext_norm = sqrt(2.0_dp*ext_norm)
2554  thr_min = ext_norm*tol/u1
2555  ind = 0
2556  thr = ext_norm
2557 
2558  DO
2559  thr = thr/u1
2560  DO
2561  l = 1
2562  DO
2563  m = l + 1
2564  DO
2565  IF ((abs(sogt(m, l)) - thr) .GE. 0.0_dp) THEN
2566  ind = 1
2567  x = 0.5_dp*(eigw(l) - eigw(m))
2568  y = -sogt(m, l)/sqrt(sogt(m, l)*sogt(m, l) + x*x)
2569  IF (x .LT. 0.0_dp) y = -y
2570 
2571  IF (y .GT. 1.0_dp) y = 1.0_dp
2572  IF (y .LT. -1.0_dp) y = -1.0_dp
2573  xy = 1.0_dp - y*y
2574  sint = y/sqrt(2.0_dp*(1.0_dp + sqrt(xy)))
2575  sint2 = sint*sint
2576  cost2 = 1.0_dp - sint2
2577  cost = sqrt(cost2)
2578  sincs = sint*cost
2579 
2580  DO i = 1, n
2581  IF ((i - m) .NE. 0) THEN
2582  IF ((i - m) .LT. 0) THEN
2583  im = m
2584  mm = i
2585  ELSE
2586  im = i
2587  mm = m
2588  END IF
2589  IF ((i - l) .NE. 0) THEN
2590  IF ((i - l) .LT. 0) THEN
2591  il = l
2592  ll = i
2593  ELSE
2594  il = i
2595  ll = l
2596  END IF
2597  x = sogt(il, ll)*cost - sogt(im, mm)*sint
2598  sogt(im, mm) = sogt(il, ll)*sint + sogt(im, mm)*cost
2599  sogt(il, ll) = x
2600  END IF
2601  END IF
2602 
2603  x = eigv(i, l)*cost - eigv(i, m)*sint
2604  eigv(i, m) = eigv(i, l)*sint + eigv(i, m)*cost
2605  eigv(i, l) = x
2606  END DO
2607 
2608  x = 2.0_dp*sogt(m, l)*sincs
2609  y = eigw(l)*cost2 + eigw(m)*sint2 - x
2610  x = eigw(l)*sint2 + eigw(m)*cost2 + x
2611  sogt(m, l) = (eigw(l) - eigw(m))*sincs + sogt(m, l)*(cost2 - sint2)
2612  eigw(l) = y
2613  eigw(m) = x
2614  END IF
2615  IF ((m - n) .EQ. 0) EXIT
2616  m = m + 1
2617  END DO
2618  IF ((l - m + 1) .EQ. 0) EXIT
2619  l = l + 1
2620  END DO
2621  IF ((ind - 1) .NE. 0.0_dp) EXIT
2622  ind = 0
2623  END DO
2624  IF ((thr - thr_min) .LE. 0.0_dp) EXIT
2625  END DO
2626  END IF
2627 
2628  IF (ic .NE. 0) THEN
2629  DO i = 1, n
2630  DO j = 1, n
2631  IF ((eigw(i) - eigw(j)) .GT. 0.0_dp) THEN
2632  x = eigw(i)
2633  eigw(i) = eigw(j)
2634  eigw(j) = x
2635  DO k = 1, n
2636  y = eigv(k, i)
2637  eigv(k, i) = eigv(k, j)
2638  eigv(k, j) = y
2639  END DO
2640  END IF
2641  END DO
2642  END DO
2643 
2644  END IF
2645 
2646  END SUBROUTINE jacob2
2647 
2648 ! **************************************************************************************************
2649 !> \brief ...
2650 !> \param n ...
2651 !> \param matrix_s_pgf ...
2652 !> \param matrix_sinv_pgf ...
2653 ! **************************************************************************************************
2654  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  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: a
2663  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  ALLOCATE (a(n))
2675  ALLOCATE (g(n, n))
2676 
2677  DO jn = 1, n
2678  diag_s = matrix_s_pgf(jn, jn)
2679  g(jn, jn) = 1.0_dp
2680 
2681  IF (jn .NE. 1) THEN
2682  DO j = 1, jn - 1
2683  scalar = 0.0_dp
2684  DO i = 1, j
2685  scalar = scalar + matrix_s_pgf(i, jn)*g(i, j)
2686  END DO
2687  diag_s = diag_s - scalar*scalar
2688  a(j) = scalar
2689  END DO
2690 
2691  DO j = 1, jn - 1
2692  row_sum = 0.0_dp
2693  DO k = j, jn - 1
2694  row_sum = row_sum + a(k)*g(j, k)
2695  END DO
2696  g(j, jn) = -row_sum
2697  END DO
2698  END IF
2699 
2700  diag_s = 1.0_dp/sqrt(diag_s)
2701  DO i = 1, jn
2702  g(i, jn) = g(i, jn)*diag_s
2703  END DO
2704  END DO
2705 
2706  DO j = 1, n
2707  DO i = 1, j
2708  matrix_sinv_pgf(j, i) = 0.0_dp
2709  matrix_sinv_pgf(i, j) = g(i, j)
2710  END DO
2711  END DO
2712 
2713  DEALLOCATE (a, g)
2714 
2715  END SUBROUTINE sog
2716 
2717 END MODULE dkh_main
dgemm
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Definition: grid_cpu_task_list.c:195
cp_blacs_env
methods related to the blacs parallel environment
Definition: cp_blacs_env.F:15
cp_fm_basic_linalg
basic linear algebra operations for full matrices
Definition: cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_column_scale
subroutine, public cp_fm_column_scale(matrixa, scaling)
scales column i of matrix a with scaling(i)
Definition: cp_fm_basic_linalg.F:1882
cp_fm_basic_linalg::cp_fm_transpose
subroutine, public cp_fm_transpose(matrix, matrixt)
transposes a matrix matrixt = matrix ^ T
Definition: cp_fm_basic_linalg.F:1710
cp_fm_basic_linalg::cp_fm_syrk
subroutine, public cp_fm_syrk(uplo, trans, k, alpha, matrix_a, ia, ja, beta, matrix_c)
performs a rank-k update of a symmetric matrix_c matrix_c = beta * matrix_c + alpha * matrix_a * tran...
Definition: cp_fm_basic_linalg.F:675
cp_fm_basic_linalg::cp_fm_upper_to_full
subroutine, public cp_fm_upper_to_full(matrix, work)
given an upper triangular matrix computes the corresponding full matrix
Definition: cp_fm_basic_linalg.F:1759
cp_fm_basic_linalg::cp_fm_schur_product
subroutine, public cp_fm_schur_product(matrix_a, matrix_b, matrix_c)
computes the schur product of two matrices c_ij = a_ij * b_ij
Definition: cp_fm_basic_linalg.F:728
cp_fm_basic_linalg::cp_fm_scale_and_add
subroutine, public cp_fm_scale_and_add(alpha, matrix_a, beta, matrix_b)
calc A <- alpha*A + beta*B optimized for alpha == 1.0 (just add beta*B) and beta == 0....
Definition: cp_fm_basic_linalg.F:167
cp_fm_basic_linalg::cp_fm_triangular_multiply
subroutine, public cp_fm_triangular_multiply(triangular_matrix, matrix_b, side, transpose_tr, invert_tr, uplo_tr, unit_diag_tr, n_rows, n_cols, alpha)
multiplies in place by a triangular matrix: matrix_b = alpha op(triangular_matrix) matrix_b or (if si...
Definition: cp_fm_basic_linalg.F:1592
cp_fm_cholesky
various cholesky decomposition related routines
Definition: cp_fm_cholesky.F:14
cp_fm_cholesky::cp_fm_cholesky_decompose
subroutine, public cp_fm_cholesky_decompose(matrix, n, info_out)
used to replace a symmetric positive def. matrix M with its cholesky decomposition U: M = U^T * U,...
Definition: cp_fm_cholesky.F:50
cp_fm_cholesky::cp_fm_cholesky_reduce
subroutine, public cp_fm_cholesky_reduce(matrix, matrixb, itype)
reduce a matrix pencil A,B to normal form B has to be cholesky decomposed with cp_fm_cholesky_decompo...
Definition: cp_fm_cholesky.F:192
cp_fm_diag
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
Definition: cp_fm_diag.F:17
cp_fm_diag::cp_fm_syevd
subroutine, public cp_fm_syevd(matrix, eigenvectors, eigenvalues, info)
Computes all eigenvalues and vectors of a real symmetric matrix significantly faster than syevx,...
Definition: cp_fm_diag.F:413
cp_fm_struct
represent the structure of a full matrix
Definition: cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition: cp_fm_struct.F:125
cp_fm_struct::cp_fm_struct_get
subroutine, public cp_fm_struct_get(fmstruct, para_env, context, descriptor, ncol_block, nrow_block, nrow_global, ncol_global, first_p_pos, row_indices, col_indices, nrow_local, ncol_local, nrow_locals, ncol_locals, local_leading_dimension)
returns the values of various attributes of the matrix structure
Definition: cp_fm_struct.F:409
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition: cp_fm_struct.F:320
cp_fm_types
represent a full matrix distributed on many processors
Definition: cp_fm_types.F:15
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition: cp_fm_types.F:1016
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition: cp_fm_types.F:167
dkh_main
Definition: dkh_main.F:7
dkh_main::dkh_atom_transformation
subroutine, public dkh_atom_transformation(s, v, h, pVp, n, dkh_order)
...
Definition: dkh_main.F:1245
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition: parallel_gemm_api.F:14
physcon
Definition of physical constants:
Definition: physcon.F:68
physcon::c_light_au
real(kind=dp), parameter, public c_light_au
Definition: physcon.F:90
qs_environment_types
Definition: qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Definition: qs_environment_types.F:502