(git:374b731)
Loading...
Searching...
No Matches
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!--------------------------------------------------------------------------------------------------!
7MODULE 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
42CONTAINS
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
2717END 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_fm_types::cp_fm_release
Definition cp_fm_types.F:90
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:85
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:38
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
cp_blacs_env::cp_blacs_env_type
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
Definition cp_blacs_env.F:53
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:116
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:215