(git:71c3ab0)
Loading...
Searching...
No Matches
qs_localization_methods.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Localization methods such as 2x2 Jacobi rotations
10!> Steepest Decents
11!> Conjugate Gradient
12!> \par History
13!> Initial parallellization of jacobi (JVDV 07.2003)
14!> direct minimization using exponential parametrization (JVDV 09.2003)
15!> crazy rotations go fast (JVDV 10.2003)
16!> \author CJM (04.2003)
17! **************************************************************************************************
18MODULE qs_localization_methods
19 USE bibliography, ONLY: schreder2024_2,&
20 cite_reference
21 USE cell_types, ONLY: cell_type
22 USE cp_blacs_env, ONLY: cp_blacs_env_type
23 USE cp_cfm_basic_linalg, ONLY: cp_cfm_column_scale,&
24 cp_cfm_gemm,&
25 cp_cfm_rot_cols,&
26 cp_cfm_rot_rows,&
27 cp_cfm_scale,&
28 cp_cfm_scale_and_add,&
29 cp_cfm_schur_product,&
30 cp_cfm_trace
31 USE cp_cfm_diag, ONLY: cp_cfm_heevd
32 USE cp_cfm_types, ONLY: &
33 cp_cfm_create, cp_cfm_get_element, cp_cfm_get_info, cp_cfm_get_submatrix, cp_cfm_release, &
34 cp_cfm_set_all, cp_cfm_set_element, cp_cfm_set_submatrix, cp_cfm_to_cfm, cp_cfm_to_fm, &
35 cp_cfm_type, cp_fm_to_cfm
36 USE cp_dbcsr_api, ONLY: dbcsr_p_type
37 USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply
38 USE cp_external_control, ONLY: external_control
39 USE cp_fm_basic_linalg, ONLY: cp_fm_frobenius_norm,&
40 cp_fm_pdgeqpf,&
41 cp_fm_pdorgqr,&
42 cp_fm_scale,&
43 cp_fm_scale_and_add,&
44 cp_fm_trace,&
45 cp_fm_transpose,&
46 cp_fm_triangular_multiply
47 USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose
48 USE cp_fm_diag, ONLY: choose_eigv_solver,&
49 cp_fm_syevd
50 USE cp_fm_struct, ONLY: cp_fm_struct_create,&
51 cp_fm_struct_get,&
52 cp_fm_struct_release,&
53 cp_fm_struct_type
54 USE cp_fm_types, ONLY: &
55 cp_fm_create, cp_fm_get_element, cp_fm_get_info, cp_fm_get_submatrix, cp_fm_init_random, &
56 cp_fm_maxabsrownorm, cp_fm_maxabsval, cp_fm_release, cp_fm_set_all, cp_fm_set_submatrix, &
57 cp_fm_to_fm, cp_fm_to_fm_submat, cp_fm_type
58 USE cp_log_handling, ONLY: cp_logger_get_default_io_unit,&
59 cp_logger_get_default_unit_nr
60 USE kahan_sum, ONLY: accurate_sum
61 USE kinds, ONLY: dp
62 USE machine, ONLY: m_flush,&
63 m_walltime
64 USE mathconstants, ONLY: gaussi,&
65 pi,&
66 twopi,&
67 z_one,&
68 z_zero
69 USE matrix_exp, ONLY: exp_pade_real,&
70 get_nsquare_norder
71 USE message_passing, ONLY: mp_para_env_type
72 USE parallel_gemm_api, ONLY: parallel_gemm
73#include "./base/base_uses.f90"
74
75 IMPLICIT NONE
76 PUBLIC :: initialize_weights, crazy_rotations, &
77 direct_mini, rotate_orbitals, approx_l1_norm_sd, jacobi_rotations, scdm_qrfact, zij_matrix, &
78 jacobi_cg_edf_ls, cardoso_souloumiac, cardoso_souloumiac_pipek
79
80 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_localization_methods'
81
82 PRIVATE
83
84 TYPE set_c_1d_type
85 COMPLEX(KIND=dp), POINTER, DIMENSION(:) :: c_array => null()
86 END TYPE set_c_1d_type
87
88 TYPE set_c_2d_type
89 COMPLEX(KIND=dp), POINTER, DIMENSION(:, :) :: c_array => null()
90 END TYPE set_c_2d_type
91
92CONTAINS
93! **************************************************************************************************
94!> \brief ...
95!> \param C ...
96!> \param iterations ...
97!> \param eps ...
98!> \param converged ...
99!> \param sweeps ...
100! **************************************************************************************************
101 SUBROUTINE approx_l1_norm_sd(C, iterations, eps, converged, sweeps)
102 TYPE(cp_fm_type), INTENT(IN) :: c
103 INTEGER, INTENT(IN) :: iterations
104 REAL(kind=dp), INTENT(IN) :: eps
105 LOGICAL, INTENT(INOUT) :: converged
106 INTEGER, INTENT(INOUT) :: sweeps
107
108 CHARACTER(len=*), PARAMETER :: routinen = 'approx_l1_norm_sd'
109 INTEGER, PARAMETER :: taylor_order = 100
110 REAL(kind=dp), PARAMETER :: alpha = 0.1_dp, f2_eps = 0.01_dp
111
112 INTEGER :: handle, i, istep, k, n, ncol_local, &
113 nrow_local, output_unit, p
114 REAL(kind=dp) :: expfactor, f2, f2old, gnorm, tnorm
115 TYPE(cp_blacs_env_type), POINTER :: context
116 TYPE(cp_fm_struct_type), POINTER :: fm_struct_k_k
117 TYPE(cp_fm_type) :: ctmp, g, gp1, gp2, u
118 TYPE(mp_para_env_type), POINTER :: para_env
119
120 CALL timeset(routinen, handle)
121
122 NULLIFY (context, para_env, fm_struct_k_k)
123
124 output_unit = cp_logger_get_default_io_unit()
125
126 CALL cp_fm_struct_get(c%matrix_struct, nrow_global=n, ncol_global=k, &
127 nrow_local=nrow_local, ncol_local=ncol_local, &
128 para_env=para_env, context=context)
129 CALL cp_fm_struct_create(fm_struct_k_k, para_env=para_env, context=context, &
130 nrow_global=k, ncol_global=k)
131 CALL cp_fm_create(ctmp, c%matrix_struct)
132 CALL cp_fm_create(u, fm_struct_k_k)
133 CALL cp_fm_create(g, fm_struct_k_k)
134 CALL cp_fm_create(gp1, fm_struct_k_k)
135 CALL cp_fm_create(gp2, fm_struct_k_k)
136 !
137 ! printing
138 IF (output_unit > 0) THEN
139 WRITE (output_unit, '(1X)')
140 WRITE (output_unit, '(2X,A)') '-----------------------------------------------------------------------------'
141 WRITE (output_unit, '(A,I5)') ' Nbr iterations =', iterations
142 WRITE (output_unit, '(A,E10.2)') ' eps convergence =', eps
143 WRITE (output_unit, '(A,I5)') ' Max Taylor order =', taylor_order
144 WRITE (output_unit, '(A,E10.2)') ' f2 eps =', f2_eps
145 WRITE (output_unit, '(A,E10.2)') ' alpha =', alpha
146 WRITE (output_unit, '(A)') ' iteration approx_l1_norm g_norm rel_err'
147 END IF
148 !
149 f2old = 0.0_dp
150 converged = .false.
151 !
152 ! Start the steepest descent
153 DO istep = 1, iterations
154 !
155 !-------------------------------------------------------------------
156 ! compute f_2
157 ! f_2(x)=(x^2+eps)^1/2
158 f2 = 0.0_dp
159 DO p = 1, ncol_local ! p
160 DO i = 1, nrow_local ! i
161 f2 = f2 + sqrt(c%local_data(i, p)**2 + f2_eps)
162 END DO
163 END DO
164 CALL c%matrix_struct%para_env%sum(f2)
165 !-------------------------------------------------------------------
166 ! compute the derivative of f_2
167 ! f_2(x)=(x^2+eps)^1/2
168 DO p = 1, ncol_local ! p
169 DO i = 1, nrow_local ! i
170 ctmp%local_data(i, p) = c%local_data(i, p)/sqrt(c%local_data(i, p)**2 + f2_eps)
171 END DO
172 END DO
173 CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, ctmp, c, 0.0_dp, g)
174 ! antisymmetrize
175 CALL cp_fm_transpose(g, u)
176 CALL cp_fm_scale_and_add(-0.5_dp, g, 0.5_dp, u)
177 !
178 !-------------------------------------------------------------------
179 !
180 gnorm = cp_fm_frobenius_norm(g)
181 !
182 ! rescale for steepest descent
183 CALL cp_fm_scale(-alpha, g)
184 !
185 ! compute unitary transform
186 ! zeroth order
187 CALL cp_fm_set_all(u, 0.0_dp, 1.0_dp)
188 ! first order
189 expfactor = 1.0_dp
190 CALL cp_fm_scale_and_add(1.0_dp, u, expfactor, g)
191 tnorm = cp_fm_frobenius_norm(g)
192 IF (tnorm > 1.0e-10_dp) THEN
193 ! other orders
194 CALL cp_fm_to_fm(g, gp1)
195 DO i = 2, taylor_order
196 ! new power of G
197 CALL parallel_gemm('N', 'N', k, k, k, 1.0_dp, g, gp1, 0.0_dp, gp2)
198 CALL cp_fm_to_fm(gp2, gp1)
199 ! add to the taylor expansion so far
200 expfactor = expfactor/real(i, kind=dp)
201 CALL cp_fm_scale_and_add(1.0_dp, u, expfactor, gp1)
202 tnorm = cp_fm_frobenius_norm(gp1)
203 IF (tnorm*expfactor < 1.0e-10_dp) EXIT
204 END DO
205 END IF
206 !
207 ! incrementaly rotate the MOs
208 CALL parallel_gemm('N', 'N', n, k, k, 1.0_dp, c, u, 0.0_dp, ctmp)
209 CALL cp_fm_to_fm(ctmp, c)
210 !
211 ! printing
212 IF (output_unit > 0) THEN
213 WRITE (output_unit, '(10X,I4,E18.10,2E10.2)') istep, f2, gnorm, abs((f2 - f2old)/f2)
214 END IF
215 !
216 ! Are we done?
217 sweeps = istep
218 IF (abs((f2 - f2old)/f2) <= eps .AND. istep > 1) THEN
219 converged = .true.
220 EXIT
221 END IF
222 f2old = f2
223 END DO
224 !
225 ! here we should do one refine step to enforce C'*S*C=1 for any case
226 !
227 ! Print the final result
228 IF (output_unit > 0) WRITE (output_unit, '(A,E16.10)') ' sparseness function f2 = ', f2
229 ! deallocate
230 CALL cp_fm_struct_release(fm_struct_k_k)
231 CALL cp_fm_release(ctmp)
232 CALL cp_fm_release(u)
233 CALL cp_fm_release(g)
234 CALL cp_fm_release(gp1)
235 CALL cp_fm_release(gp2)
236
237 CALL timestop(handle)
238
239 END SUBROUTINE approx_l1_norm_sd
240! **************************************************************************************************
241!> \brief ...
242!> \param cell ...
243!> \param weights ...
244! **************************************************************************************************
245 SUBROUTINE initialize_weights(cell, weights)
246
247 TYPE(cell_type), POINTER :: cell
248 REAL(kind=dp), DIMENSION(:) :: weights
249
250 REAL(kind=dp), DIMENSION(3, 3) :: metric
251
252 cpassert(ASSOCIATED(cell))
253
254 metric = 0.0_dp
255 CALL dgemm('T', 'N', 3, 3, 3, 1._dp, cell%hmat(:, :), 3, cell%hmat(:, :), 3, 0.0_dp, metric(:, :), 3)
256
257 weights(1) = metric(1, 1) - metric(1, 2) - metric(1, 3)
258 weights(2) = metric(2, 2) - metric(1, 2) - metric(2, 3)
259 weights(3) = metric(3, 3) - metric(1, 3) - metric(2, 3)
260 weights(4) = metric(1, 2)
261 weights(5) = metric(1, 3)
262 weights(6) = metric(2, 3)
263
264 END SUBROUTINE initialize_weights
265
266! **************************************************************************************************
267!> \brief wrapper for the jacobi routines, should be removed if jacobi_rot_para
268!> can deal with serial para_envs.
269!> \param weights ...
270!> \param zij ...
271!> \param vectors ...
272!> \param para_env ...
273!> \param max_iter ...
274!> \param eps_localization ...
275!> \param sweeps ...
276!> \param out_each ...
277!> \param target_time ...
278!> \param start_time ...
279!> \param restricted ...
280!> \par History
281!> \author Joost VandeVondele (02.2010)
282! **************************************************************************************************
283 SUBROUTINE jacobi_rotations(weights, zij, vectors, para_env, max_iter, &
284 eps_localization, sweeps, out_each, target_time, start_time, restricted)
285
286 REAL(kind=dp), INTENT(IN) :: weights(:)
287 TYPE(cp_fm_type), INTENT(IN) :: zij(:, :), vectors
288 TYPE(mp_para_env_type), POINTER :: para_env
289 INTEGER, INTENT(IN) :: max_iter
290 REAL(kind=dp), INTENT(IN) :: eps_localization
291 INTEGER :: sweeps
292 INTEGER, INTENT(IN) :: out_each
293 REAL(dp) :: target_time, start_time
294 INTEGER :: restricted
295
296 IF (para_env%num_pe == 1) THEN
297 CALL jacobi_rotations_serial(weights, zij, vectors, max_iter, eps_localization, &
298 sweeps, out_each, restricted=restricted)
299 ELSE
300 CALL jacobi_rot_para(weights, zij, vectors, para_env, max_iter, eps_localization, &
301 sweeps, out_each, target_time, start_time, restricted=restricted)
302 END IF
303
304 END SUBROUTINE jacobi_rotations
305
306! **************************************************************************************************
307!> \brief this routine, private to the module is a serial backup, till we have jacobi_rot_para to work in serial
308!> while the routine below works in parallel, it is too slow to be useful
309!> \param weights ...
310!> \param zij ...
311!> \param vectors ...
312!> \param max_iter ...
313!> \param eps_localization ...
314!> \param sweeps ...
315!> \param out_each ...
316!> \param restricted ...
317! **************************************************************************************************
318 SUBROUTINE jacobi_rotations_serial(weights, zij, vectors, max_iter, eps_localization, sweeps, &
319 out_each, restricted)
320 REAL(kind=dp), INTENT(IN) :: weights(:)
321 TYPE(cp_fm_type), INTENT(IN) :: zij(:, :), vectors
322 INTEGER, INTENT(IN) :: max_iter
323 REAL(kind=dp), INTENT(IN) :: eps_localization
324 INTEGER :: sweeps
325 INTEGER, INTENT(IN) :: out_each
326 INTEGER :: restricted
327
328 CHARACTER(len=*), PARAMETER :: routinen = 'jacobi_rotations_serial'
329
330 COMPLEX(KIND=dp), POINTER :: mii(:), mij(:), mjj(:)
331 INTEGER :: dim2, handle, idim, istate, jstate, &
332 nstate, unit_nr
333 REAL(kind=dp) :: ct, st, t1, t2, theta, tolerance
334 TYPE(cp_cfm_type) :: c_rmat
335 TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: c_zij
336 TYPE(cp_fm_type) :: rmat
337
338 CALL timeset(routinen, handle)
339
340 dim2 = SIZE(zij, 2)
341 ALLOCATE (c_zij(dim2))
342 NULLIFY (mii, mij, mjj)
343 ALLOCATE (mii(dim2), mij(dim2), mjj(dim2))
344
345 CALL cp_fm_create(rmat, zij(1, 1)%matrix_struct)
346 CALL cp_fm_set_all(rmat, 0._dp, 1._dp)
347
348 CALL cp_cfm_create(c_rmat, zij(1, 1)%matrix_struct)
349 CALL cp_cfm_set_all(c_rmat, (0._dp, 0._dp), (1._dp, 0._dp))
350 DO idim = 1, dim2
351 CALL cp_cfm_create(c_zij(idim), zij(1, 1)%matrix_struct)
352 c_zij(idim)%local_data = cmplx(zij(1, idim)%local_data, &
353 zij(2, idim)%local_data, dp)
354 END DO
355
356 CALL cp_fm_get_info(rmat, nrow_global=nstate)
357 tolerance = 1.0e10_dp
358
359 sweeps = 0
360 unit_nr = -1
361 IF (rmat%matrix_struct%para_env%is_source()) THEN
362 unit_nr = cp_logger_get_default_unit_nr()
363 WRITE (unit_nr, '(T4,A )') " Localization by iterative Jacobi rotation"
364 END IF
365
366 IF (restricted > 0) THEN
367 unit_nr = cp_logger_get_default_unit_nr()
368 WRITE (unit_nr, '(T4,A,I2,A )') "JACOBI: for the ROKS method, the last ", restricted, " orbitals DO NOT ROTATE"
369 nstate = nstate - restricted
370 END IF
371
372 ! do jacobi sweeps until converged
373 DO WHILE (tolerance >= eps_localization .AND. sweeps < max_iter)
374 sweeps = sweeps + 1
375 t1 = m_walltime()
376
377 DO istate = 1, nstate
378 DO jstate = istate + 1, nstate
379 DO idim = 1, dim2
380 CALL cp_cfm_get_element(c_zij(idim), istate, istate, mii(idim))
381 CALL cp_cfm_get_element(c_zij(idim), istate, jstate, mij(idim))
382 CALL cp_cfm_get_element(c_zij(idim), jstate, jstate, mjj(idim))
383 END DO
384 CALL get_angle(mii, mjj, mij, weights, theta)
385 st = sin(theta)
386 ct = cos(theta)
387 CALL rotate_zij(istate, jstate, st, ct, c_zij)
388
389 CALL rotate_rmat(istate, jstate, st, ct, c_rmat)
390 END DO
391 END DO
392
393 CALL check_tolerance(c_zij, weights, tolerance)
394
395 t2 = m_walltime()
396 IF (unit_nr > 0 .AND. modulo(sweeps, out_each) == 0) THEN
397 WRITE (unit_nr, '(T4,A,I7,T30,A,E12.4,T60,A,F8.3)') &
398 "Iteration:", sweeps, "Tolerance:", tolerance, "Time:", t2 - t1
399 CALL m_flush(unit_nr)
400 END IF
401
402 END DO
403
404 DO idim = 1, dim2
405 zij(1, idim)%local_data = real(c_zij(idim)%local_data, dp)
406 zij(2, idim)%local_data = aimag(c_zij(idim)%local_data)
407 CALL cp_cfm_release(c_zij(idim))
408 END DO
409 DEALLOCATE (c_zij)
410 DEALLOCATE (mii, mij, mjj)
411 rmat%local_data = real(c_rmat%local_data, dp)
412
413 CALL rotate_orbitals(rmat, vectors)
414
415 CALL cp_cfm_release(c_rmat)
416 CALL cp_fm_release(rmat)
417
418 CALL timestop(handle)
419
420 END SUBROUTINE jacobi_rotations_serial
421! **************************************************************************************************
422!> \brief very similar to jacobi_rotations_serial with some extra output options
423!> \param weights ...
424!> \param c_zij ...
425!> \param max_iter ...
426!> \param c_rmat ...
427!> \param eps_localization ...
428!> \param tol_out ...
429!> \param jsweeps ...
430!> \param out_each ...
431!> \param c_zij_out ...
432!> \param grad_final ...
433! **************************************************************************************************
434 SUBROUTINE jacobi_rotations_serial_1(weights, c_zij, max_iter, c_rmat, eps_localization, &
435 tol_out, jsweeps, out_each, c_zij_out, grad_final)
436 REAL(kind=dp), INTENT(IN) :: weights(:)
437 TYPE(cp_cfm_type), INTENT(IN) :: c_zij(:)
438 INTEGER, INTENT(IN) :: max_iter
439 TYPE(cp_cfm_type), INTENT(IN) :: c_rmat
440 REAL(kind=dp), INTENT(IN), OPTIONAL :: eps_localization
441 REAL(kind=dp), INTENT(OUT), OPTIONAL :: tol_out
442 INTEGER, INTENT(OUT), OPTIONAL :: jsweeps
443 INTEGER, INTENT(IN), OPTIONAL :: out_each
444 TYPE(cp_cfm_type), INTENT(IN), OPTIONAL :: c_zij_out(:)
445 TYPE(cp_fm_type), INTENT(OUT), OPTIONAL, POINTER :: grad_final
446
447 CHARACTER(len=*), PARAMETER :: routinen = 'jacobi_rotations_serial_1'
448
449 COMPLEX(KIND=dp) :: mzii
450 COMPLEX(KIND=dp), POINTER :: mii(:), mij(:), mjj(:)
451 INTEGER :: dim2, handle, idim, istate, jstate, &
452 nstate, sweeps, unit_nr
453 REAL(kind=dp) :: alpha, avg_spread_ii, ct, spread_ii, st, &
454 sum_spread_ii, t1, t2, theta, tolerance
455 TYPE(cp_cfm_type) :: c_rmat_local
456 TYPE(cp_cfm_type), ALLOCATABLE :: c_zij_local(:)
457
458 CALL timeset(routinen, handle)
459
460 dim2 = SIZE(c_zij)
461 NULLIFY (mii, mij, mjj)
462 ALLOCATE (mii(dim2), mij(dim2), mjj(dim2))
463
464 ALLOCATE (c_zij_local(dim2))
465 CALL cp_cfm_create(c_rmat_local, c_rmat%matrix_struct)
466 CALL cp_cfm_set_all(c_rmat_local, (0.0_dp, 0.0_dp), (1.0_dp, 0.0_dp))
467 DO idim = 1, dim2
468 CALL cp_cfm_create(c_zij_local(idim), c_zij(idim)%matrix_struct)
469 c_zij_local(idim)%local_data = c_zij(idim)%local_data
470 END DO
471
472 CALL cp_cfm_get_info(c_rmat_local, nrow_global=nstate)
473 tolerance = 1.0e10_dp
474
475 IF (PRESENT(grad_final)) CALL cp_fm_set_all(grad_final, 0.0_dp)
476
477 sweeps = 0
478 IF (PRESENT(out_each)) THEN
479 unit_nr = -1
480 IF (c_rmat_local%matrix_struct%para_env%is_source()) THEN
481 unit_nr = cp_logger_get_default_unit_nr()
482 END IF
483 alpha = 0.0_dp
484 DO idim = 1, dim2
485 alpha = alpha + weights(idim)
486 END DO
487 END IF
488
489 ! do jacobi sweeps until converged
490 DO WHILE (sweeps < max_iter)
491 sweeps = sweeps + 1
492 IF (PRESENT(eps_localization)) THEN
493 IF (tolerance < eps_localization) EXIT
494 END IF
495 IF (PRESENT(out_each)) t1 = m_walltime()
496
497 DO istate = 1, nstate
498 DO jstate = istate + 1, nstate
499 DO idim = 1, dim2
500 CALL cp_cfm_get_element(c_zij_local(idim), istate, istate, mii(idim))
501 CALL cp_cfm_get_element(c_zij_local(idim), istate, jstate, mij(idim))
502 CALL cp_cfm_get_element(c_zij_local(idim), jstate, jstate, mjj(idim))
503 END DO
504 CALL get_angle(mii, mjj, mij, weights, theta)
505 st = sin(theta)
506 ct = cos(theta)
507 CALL rotate_zij(istate, jstate, st, ct, c_zij_local)
508
509 CALL rotate_rmat(istate, jstate, st, ct, c_rmat_local)
510 END DO
511 END DO
512
513 IF (PRESENT(grad_final)) THEN
514 CALL check_tolerance(c_zij_local, weights, tolerance, grad=grad_final)
515 ELSE
516 CALL check_tolerance(c_zij_local, weights, tolerance)
517 END IF
518 IF (PRESENT(tol_out)) tol_out = tolerance
519
520 IF (PRESENT(out_each)) THEN
521 t2 = m_walltime()
522 IF (unit_nr > 0 .AND. modulo(sweeps, out_each) == 0) THEN
523 sum_spread_ii = 0.0_dp
524 DO istate = 1, nstate
525 spread_ii = 0.0_dp
526 DO idim = 1, dim2
527 CALL cp_cfm_get_element(c_zij_local(idim), istate, istate, mzii)
528 spread_ii = spread_ii + weights(idim)* &
529 abs(mzii)**2/twopi/twopi
530 END DO
531 sum_spread_ii = sum_spread_ii + spread_ii
532 END DO
533 sum_spread_ii = alpha*nstate/twopi/twopi - sum_spread_ii
534 avg_spread_ii = sum_spread_ii/nstate
535 WRITE (unit_nr, '(T4,A,T26,A,T48,A,T64,A)') &
536 "Iteration", "Avg. Spread_ii", "Tolerance", "Time"
537 WRITE (unit_nr, '(T4,I7,T20,F20.10,T45,E12.4,T60,F8.3)') &
538 sweeps, avg_spread_ii, tolerance, t2 - t1
539 CALL m_flush(unit_nr)
540 END IF
541 IF (PRESENT(jsweeps)) jsweeps = sweeps
542 END IF
543
544 END DO
545
546 IF (PRESENT(c_zij_out)) THEN
547 DO idim = 1, dim2
548 CALL cp_cfm_to_cfm(c_zij_local(idim), c_zij_out(idim))
549 END DO
550 END IF
551 CALL cp_cfm_to_cfm(c_rmat_local, c_rmat)
552
553 DEALLOCATE (mii, mij, mjj)
554 DO idim = 1, dim2
555 CALL cp_cfm_release(c_zij_local(idim))
556 END DO
557 DEALLOCATE (c_zij_local)
558 CALL cp_cfm_release(c_rmat_local)
559
560 CALL timestop(handle)
561
562 END SUBROUTINE jacobi_rotations_serial_1
563! **************************************************************************************************
564!> \brief combine jacobi rotations (serial) and conjugate gradient with golden section line search
565!> for partially occupied wannier functions
566!> \param para_env ...
567!> \param weights ...
568!> \param zij ...
569!> \param vectors ...
570!> \param max_iter ...
571!> \param eps_localization ...
572!> \param iter ...
573!> \param out_each ...
574!> \param nextra ...
575!> \param do_cg ...
576!> \param nmo ...
577!> \param vectors_2 ...
578!> \param mos_guess ...
579! **************************************************************************************************
580 SUBROUTINE jacobi_cg_edf_ls(para_env, weights, zij, vectors, max_iter, eps_localization, &
581 iter, out_each, nextra, do_cg, nmo, vectors_2, mos_guess)
582 TYPE(mp_para_env_type), POINTER :: para_env
583 REAL(kind=dp), INTENT(IN) :: weights(:)
584 TYPE(cp_fm_type), INTENT(IN) :: zij(:, :), vectors
585 INTEGER, INTENT(IN) :: max_iter
586 REAL(kind=dp), INTENT(IN) :: eps_localization
587 INTEGER :: iter
588 INTEGER, INTENT(IN) :: out_each, nextra
589 LOGICAL, INTENT(IN) :: do_cg
590 INTEGER, INTENT(IN), OPTIONAL :: nmo
591 TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: vectors_2, mos_guess
592
593 CHARACTER(len=*), PARAMETER :: routinen = 'jacobi_cg_edf_ls'
594 COMPLEX(KIND=dp), PARAMETER :: cone = (1.0_dp, 0.0_dp), &
595 czero = (0.0_dp, 0.0_dp)
596 REAL(kind=dp), PARAMETER :: gold_sec = 0.3819_dp
597
598 COMPLEX(KIND=dp) :: cnorm2_gct, cnorm2_gct_cross, mzii
599 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: tmp_cmat
600 COMPLEX(KIND=dp), DIMENSION(:), POINTER :: arr_zii
601 COMPLEX(KIND=dp), DIMENSION(:, :), POINTER :: matrix_zii
602 INTEGER :: dim2, handle, icinit, idim, istate, line_search_count, line_searches, lsl, lsm, &
603 lsr, miniter, nao, ndummy, nocc, norextra, northo, nstate, unit_nr
604 INTEGER, DIMENSION(1) :: iloc
605 LOGICAL :: do_cinit_mo, do_cinit_random, &
606 do_u_guess_mo, new_direction
607 REAL(kind=dp) :: alpha, avg_spread_ii, beta, beta_pr, ds, ds_min, mintol, norm, norm2_gct, &
608 norm2_gct_cross, norm2_old, spread_ii, spread_sum, sum_spread_ii, t1, tol, tolc, weight
609 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: sum_spread
610 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: tmp_mat, tmp_mat_1
611 REAL(kind=dp), DIMENSION(50) :: energy, pos
612 REAL(kind=dp), DIMENSION(:), POINTER :: tmp_arr
613 TYPE(cp_blacs_env_type), POINTER :: context
614 TYPE(cp_cfm_type) :: c_tilde, ctrans_lambda, gct_old, &
615 grad_ctilde, skc, tmp_cfm, tmp_cfm_1, &
616 tmp_cfm_2, u, ul, v, vl, zdiag
617 TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: c_zij, zij_0
618 TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct
619 TYPE(cp_fm_type) :: id_nextra, matrix_u, matrix_v, &
620 matrix_v_all, rmat, tmp_fm, vectors_all
621
622 CALL timeset(routinen, handle)
623
624 dim2 = SIZE(zij, 2)
625 NULLIFY (context)
626 NULLIFY (matrix_zii, arr_zii)
627 NULLIFY (tmp_fm_struct)
628 NULLIFY (tmp_arr)
629
630 ALLOCATE (c_zij(dim2))
631
632 CALL cp_fm_get_info(zij(1, 1), nrow_global=nstate)
633
634 ALLOCATE (sum_spread(nstate))
635 ALLOCATE (matrix_zii(nstate, dim2))
636 matrix_zii = czero
637 sum_spread = 0.0_dp
638
639 alpha = 0.0_dp
640 DO idim = 1, dim2
641 alpha = alpha + weights(idim)
642 CALL cp_cfm_create(c_zij(idim), zij(1, 1)%matrix_struct)
643 c_zij(idim)%local_data = cmplx(zij(1, idim)%local_data, &
644 zij(2, idim)%local_data, dp)
645 END DO
646
647 ALLOCATE (zij_0(dim2))
648
649 CALL cp_cfm_create(u, zij(1, 1)%matrix_struct)
650 CALL cp_fm_create(matrix_u, zij(1, 1)%matrix_struct)
651
652 CALL cp_cfm_set_all(u, czero, cone)
653 CALL cp_fm_set_all(matrix_u, 0.0_dp, 1.0_dp)
654
655 CALL cp_fm_get_info(vectors, nrow_global=nao)
656 IF (nextra > 0) THEN
657 IF (PRESENT(mos_guess)) THEN
658 do_cinit_random = .false.
659 do_cinit_mo = .true.
660 CALL cp_fm_get_info(mos_guess, ncol_global=ndummy)
661 ELSE
662 do_cinit_random = .true.
663 do_cinit_mo = .false.
664 ndummy = nstate
665 END IF
666
667 IF (do_cinit_random) THEN
668 icinit = 1
669 do_u_guess_mo = .false.
670 ELSE IF (do_cinit_mo) THEN
671 icinit = 2
672 do_u_guess_mo = .true.
673 END IF
674
675 nocc = nstate - nextra
676 northo = nmo - nocc
677 norextra = nmo - nstate
678 CALL cp_fm_struct_get(zij(1, 1)%matrix_struct, context=context)
679
680 ALLOCATE (tmp_cmat(nstate, nstate))
681 CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nmo, ncol_global=nmo, &
682 para_env=para_env, context=context)
683 DO idim = 1, dim2
684 CALL cp_cfm_create(zij_0(idim), tmp_fm_struct)
685 CALL cp_cfm_set_all(zij_0(idim), czero, cone)
686 CALL cp_cfm_get_submatrix(c_zij(idim), tmp_cmat)
687 CALL cp_cfm_set_submatrix(zij_0(idim), tmp_cmat)
688 END DO
689 CALL cp_fm_struct_release(tmp_fm_struct)
690 DEALLOCATE (tmp_cmat)
691
692 CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nmo, ncol_global=nstate, &
693 para_env=para_env, context=context)
694 CALL cp_cfm_create(v, tmp_fm_struct)
695 CALL cp_fm_create(matrix_v, tmp_fm_struct)
696 CALL cp_cfm_create(zdiag, tmp_fm_struct)
697 CALL cp_fm_create(rmat, tmp_fm_struct)
698 CALL cp_fm_struct_release(tmp_fm_struct)
699 CALL cp_cfm_set_all(v, czero, cone)
700 CALL cp_fm_set_all(matrix_v, 0.0_dp, 1.0_dp)
701
702 CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nmo, ncol_global=ndummy, &
703 para_env=para_env, context=context)
704 CALL cp_fm_create(matrix_v_all, tmp_fm_struct)
705 CALL cp_fm_struct_release(tmp_fm_struct)
706 CALL cp_fm_set_all(matrix_v_all, 0._dp, 1._dp)
707
708 ALLOCATE (arr_zii(nstate))
709
710 CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=northo, ncol_global=nextra, &
711 para_env=para_env, context=context)
712 CALL cp_cfm_create(c_tilde, tmp_fm_struct)
713 CALL cp_cfm_create(grad_ctilde, tmp_fm_struct)
714 CALL cp_cfm_create(gct_old, tmp_fm_struct)
715 CALL cp_cfm_create(skc, tmp_fm_struct)
716 CALL cp_fm_struct_release(tmp_fm_struct)
717 CALL cp_cfm_set_all(c_tilde, czero)
718 CALL cp_cfm_set_all(gct_old, czero)
719 CALL cp_cfm_set_all(skc, czero)
720
721 CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=northo, ncol_global=nstate, &
722 para_env=para_env, context=context)
723 CALL cp_cfm_create(vl, tmp_fm_struct)
724 CALL cp_cfm_set_all(vl, czero)
725 CALL cp_fm_struct_release(tmp_fm_struct)
726
727 CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nextra, ncol_global=nextra, &
728 para_env=para_env, context=context)
729 CALL cp_fm_create(id_nextra, tmp_fm_struct)
730 CALL cp_cfm_create(ctrans_lambda, tmp_fm_struct)
731 CALL cp_fm_struct_release(tmp_fm_struct)
732 CALL cp_cfm_set_all(ctrans_lambda, czero)
733 CALL cp_fm_set_all(id_nextra, 0.0_dp, 1.0_dp)
734
735 CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nextra, ncol_global=nstate, &
736 para_env=para_env, context=context)
737 CALL cp_cfm_create(ul, tmp_fm_struct)
738 CALL cp_fm_struct_release(tmp_fm_struct)
739 CALL cp_cfm_set_all(ul, czero)
740
741 CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, ncol_global=nmo, &
742 para_env=para_env, context=context)
743 CALL cp_fm_create(vectors_all, tmp_fm_struct)
744 CALL cp_fm_struct_release(tmp_fm_struct)
745 ALLOCATE (tmp_mat(nao, nstate))
746 CALL cp_fm_get_submatrix(vectors, tmp_mat)
747 CALL cp_fm_set_submatrix(vectors_all, tmp_mat, 1, 1, nao, nstate)
748 DEALLOCATE (tmp_mat)
749 ALLOCATE (tmp_mat(nao, norextra))
750 CALL cp_fm_get_submatrix(vectors_2, tmp_mat)
751 CALL cp_fm_set_submatrix(vectors_all, tmp_mat, 1, nstate + 1, nao, norextra)
752 DEALLOCATE (tmp_mat)
753
754 ! initialize c_tilde
755 SELECT CASE (icinit)
756 CASE (1) ! random coefficients
757 !WRITE (*, *) "RANDOM INITIAL GUESS FOR C"
758 CALL cp_fm_create(tmp_fm, c_tilde%matrix_struct)
759 CALL cp_fm_init_random(tmp_fm, nextra)
760 CALL ortho_vectors(tmp_fm)
761 c_tilde%local_data = tmp_fm%local_data
762 CALL cp_fm_release(tmp_fm)
763 ALLOCATE (tmp_cmat(northo, nextra))
764 CALL cp_cfm_get_submatrix(c_tilde, tmp_cmat)
765 CALL cp_cfm_set_submatrix(v, tmp_cmat, nocc + 1, nocc + 1, northo, nextra)
766 DEALLOCATE (tmp_cmat)
767 CASE (2) ! MO based coeffs
768 CALL parallel_gemm("T", "N", nmo, ndummy, nao, 1.0_dp, vectors_all, mos_guess, 0.0_dp, matrix_v_all)
769 ALLOCATE (tmp_arr(nmo))
770 ALLOCATE (tmp_mat(nmo, ndummy))
771 ALLOCATE (tmp_mat_1(nmo, nstate))
772 ! normalize matrix_V_all
773 CALL cp_fm_get_submatrix(matrix_v_all, tmp_mat)
774 DO istate = 1, ndummy
775 tmp_arr(:) = tmp_mat(:, istate)
776 norm = norm2(tmp_arr)
777 tmp_arr(:) = tmp_arr(:)/norm
778 tmp_mat(:, istate) = tmp_arr(:)
779 END DO
780 CALL cp_fm_set_submatrix(matrix_v_all, tmp_mat)
781 CALL cp_fm_get_submatrix(matrix_v_all, tmp_mat_1, 1, 1, nmo, nstate)
782 CALL cp_fm_set_submatrix(matrix_v, tmp_mat_1)
783 DEALLOCATE (tmp_arr, tmp_mat, tmp_mat_1)
784 CALL cp_fm_to_cfm(msourcer=matrix_v, mtarget=v)
785 ALLOCATE (tmp_mat(northo, ndummy))
786 ALLOCATE (tmp_mat_1(northo, nextra))
787 CALL cp_fm_get_submatrix(matrix_v_all, tmp_mat, nocc + 1, 1, northo, ndummy)
788 ALLOCATE (tmp_arr(ndummy))
789 tmp_arr = 0.0_dp
790 DO istate = 1, ndummy
791 tmp_arr(istate) = norm2(tmp_mat(:, istate))
792 END DO
793 ! find edfs
794 DO istate = 1, nextra
795 iloc = maxloc(tmp_arr)
796 tmp_mat_1(:, istate) = tmp_mat(:, iloc(1))
797 tmp_arr(iloc(1)) = 0.0_dp
798 END DO
799
800 DEALLOCATE (tmp_arr, tmp_mat)
801
802 CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=northo, ncol_global=nextra, &
803 para_env=para_env, context=context)
804 CALL cp_fm_create(tmp_fm, tmp_fm_struct)
805 CALL cp_fm_struct_release(tmp_fm_struct)
806 CALL cp_fm_set_submatrix(tmp_fm, tmp_mat_1)
807 DEALLOCATE (tmp_mat_1)
808 CALL ortho_vectors(tmp_fm)
809 CALL cp_fm_to_cfm(msourcer=tmp_fm, mtarget=c_tilde)
810 CALL cp_fm_release(tmp_fm)
811 ! initialize U
812 IF (do_u_guess_mo) THEN
813 ALLOCATE (tmp_cmat(nocc, nstate))
814 CALL cp_cfm_get_submatrix(v, tmp_cmat, 1, 1, nocc, nstate)
815 CALL cp_cfm_set_submatrix(u, tmp_cmat, 1, 1, nocc, nstate)
816 DEALLOCATE (tmp_cmat)
817 ALLOCATE (tmp_cmat(northo, nstate))
818 CALL cp_cfm_get_submatrix(v, tmp_cmat, nocc + 1, 1, northo, nstate)
819 CALL cp_cfm_set_submatrix(vl, tmp_cmat, 1, 1, northo, nstate)
820 DEALLOCATE (tmp_cmat)
821 CALL parallel_gemm("C", "N", nextra, nstate, northo, cone, c_tilde, vl, czero, ul)
822 ALLOCATE (tmp_cmat(nextra, nstate))
823 CALL cp_cfm_get_submatrix(ul, tmp_cmat, 1, 1, nextra, nstate)
824 CALL cp_cfm_set_submatrix(u, tmp_cmat, nocc + 1, 1, nextra, nstate)
825 DEALLOCATE (tmp_cmat)
826 CALL cp_fm_create(tmp_fm, u%matrix_struct)
827 tmp_fm%local_data = real(u%local_data, kind=dp)
828 CALL ortho_vectors(tmp_fm)
829 CALL cp_fm_to_cfm(msourcer=tmp_fm, mtarget=u)
830 CALL cp_fm_release(tmp_fm)
831 CALL cp_cfm_to_fm(u, matrix_u)
832 END IF
833 ! reevaluate V
834 ALLOCATE (tmp_cmat(nocc, nstate))
835 CALL cp_cfm_get_submatrix(u, tmp_cmat, 1, 1, nocc, nstate)
836 CALL cp_cfm_set_submatrix(v, tmp_cmat, 1, 1, nocc, nstate)
837 DEALLOCATE (tmp_cmat)
838 ALLOCATE (tmp_cmat(nextra, nstate))
839 CALL cp_cfm_get_submatrix(u, tmp_cmat, nocc + 1, 1, nextra, nstate)
840 CALL cp_cfm_set_submatrix(ul, tmp_cmat, 1, 1, nextra, nstate)
841 DEALLOCATE (tmp_cmat)
842 CALL parallel_gemm("N", "N", northo, nstate, nextra, cone, c_tilde, ul, czero, vl)
843 ALLOCATE (tmp_cmat(northo, nstate))
844 CALL cp_cfm_get_submatrix(vl, tmp_cmat)
845 CALL cp_cfm_set_submatrix(v, tmp_cmat, nocc + 1, 1, northo, nstate)
846 DEALLOCATE (tmp_cmat)
847 END SELECT
848 ELSE
849 DO idim = 1, dim2
850 CALL cp_cfm_create(zij_0(idim), zij(1, 1)%matrix_struct)
851 CALL cp_cfm_to_cfm(c_zij(idim), zij_0(idim))
852 END DO
853 CALL cp_fm_create(rmat, zij(1, 1)%matrix_struct)
854 CALL cp_fm_set_all(rmat, 0._dp, 1._dp)
855 END IF
856
857 unit_nr = -1
858 IF (rmat%matrix_struct%para_env%is_source()) THEN
859 unit_nr = cp_logger_get_default_unit_nr()
860 WRITE (unit_nr, '(T4,A )') " Localization by combined Jacobi rotations and Non-Linear Conjugate Gradient"
861 END IF
862
863 norm2_old = 1.0e30_dp
864 ds_min = 1.0_dp
865 new_direction = .true.
866 iter = 0
867 line_searches = 0
868 line_search_count = 0
869 tol = 1.0e+20_dp
870 mintol = 1.0e+10_dp
871 miniter = 0
872
873 !IF (nextra > 0) WRITE(*,*) 'random_guess, MO_guess, U_guess, conjugate_gradient: ', &
874 ! do_cinit_random, do_cinit_mo, do_U_guess_mo, do_cg
875
876 ! do conjugate gradient until converged
877 DO WHILE (iter < max_iter)
878 iter = iter + 1
879 !WRITE(*,*) 'iter = ', iter
880 t1 = m_walltime()
881
882 IF (iter > 1) THEN
883 ! comput U
884 CALL cp_cfm_create(tmp_cfm, zij(1, 1)%matrix_struct)
885 CALL cp_cfm_create(tmp_cfm_2, zij(1, 1)%matrix_struct)
886 IF (para_env%num_pe == 1) THEN
887 CALL jacobi_rotations_serial_1(weights, c_zij, 1, tmp_cfm_2, tol_out=tol)
888 ELSE
889 CALL jacobi_rot_para_1(weights, c_zij, para_env, 1, tmp_cfm_2, tol_out=tol)
890 END IF
891 CALL parallel_gemm('N', 'N', nstate, nstate, nstate, cone, u, tmp_cfm_2, czero, tmp_cfm)
892 CALL cp_cfm_to_cfm(tmp_cfm, u)
893 CALL cp_cfm_release(tmp_cfm)
894 CALL cp_cfm_release(tmp_cfm_2)
895 END IF
896
897 IF (nextra > 0) THEN
898 ALLOCATE (tmp_cmat(nextra, nstate))
899 CALL cp_cfm_get_submatrix(u, tmp_cmat, nocc + 1, 1, nextra, nstate)
900 CALL cp_cfm_set_submatrix(ul, tmp_cmat)
901 DEALLOCATE (tmp_cmat)
902 IF (iter > 1) THEN
903 ! orthonormalize c_tilde
904 CALL cp_fm_create(tmp_fm, c_tilde%matrix_struct)
905 tmp_fm%local_data = real(c_tilde%local_data, kind=dp)
906 CALL ortho_vectors(tmp_fm)
907 CALL cp_fm_to_cfm(msourcer=tmp_fm, mtarget=c_tilde)
908 CALL cp_fm_release(tmp_fm)
909
910 ALLOCATE (tmp_cmat(nocc, nstate))
911 CALL cp_cfm_get_submatrix(u, tmp_cmat, 1, 1, nocc, nstate)
912 CALL cp_cfm_set_submatrix(v, tmp_cmat, 1, 1, nocc, nstate)
913 DEALLOCATE (tmp_cmat)
914 CALL parallel_gemm("N", "N", northo, nstate, nextra, cone, c_tilde, ul, czero, vl)
915 ALLOCATE (tmp_cmat(northo, nstate))
916 CALL cp_cfm_get_submatrix(vl, tmp_cmat)
917 CALL cp_cfm_set_submatrix(v, tmp_cmat, nocc + 1, 1, northo, nstate)
918 DEALLOCATE (tmp_cmat)
919 END IF
920
921 ! reset if new_direction
922 IF (new_direction .AND. mod(line_searches, 20) == 5) THEN
923 CALL cp_cfm_set_all(skc, czero)
924 CALL cp_cfm_set_all(gct_old, czero)
925 norm2_old = 1.0e30_dp
926 END IF
927
928 CALL cp_cfm_create(tmp_cfm, v%matrix_struct)
929 CALL cp_cfm_to_cfm(v, tmp_cfm)
930 CALL cp_cfm_create(tmp_cfm_1, v%matrix_struct)
931 ndummy = nmo
932 ELSE
933 CALL cp_cfm_create(tmp_cfm, zij(1, 1)%matrix_struct)
934 CALL cp_cfm_to_cfm(u, tmp_cfm)
935 CALL cp_cfm_create(tmp_cfm_1, zij(1, 1)%matrix_struct)
936 ndummy = nstate
937 END IF
938 ! update z_ij
939 DO idim = 1, dim2
940 ! 'tmp_cfm_1 = zij_0*tmp_cfm'
941 CALL parallel_gemm("N", "N", ndummy, nstate, ndummy, cone, zij_0(idim), &
942 tmp_cfm, czero, tmp_cfm_1)
943 ! 'c_zij = tmp_cfm_dagg*tmp_cfm_1'
944 CALL parallel_gemm("C", "N", nstate, nstate, ndummy, cone, tmp_cfm, tmp_cfm_1, &
945 czero, c_zij(idim))
946 END DO
947 CALL cp_cfm_release(tmp_cfm)
948 CALL cp_cfm_release(tmp_cfm_1)
949 ! compute spread
950 DO istate = 1, nstate
951 spread_ii = 0.0_dp
952 DO idim = 1, dim2
953 CALL cp_cfm_get_element(c_zij(idim), istate, istate, mzii)
954 spread_ii = spread_ii + weights(idim)* &
955 abs(mzii)**2/twopi/twopi
956 matrix_zii(istate, idim) = mzii
957 END DO
958 !WRITE(*,*) 'spread_ii', spread_ii
959 sum_spread(istate) = spread_ii
960 END DO
961 CALL c_zij(1)%matrix_struct%para_env%sum(spread_ii)
962 spread_sum = accurate_sum(sum_spread)
963
964 IF (nextra > 0) THEN
965 ! update c_tilde
966 CALL cp_cfm_set_all(zdiag, czero)
967 CALL cp_cfm_set_all(grad_ctilde, czero)
968 CALL cp_cfm_create(tmp_cfm, v%matrix_struct)
969 CALL cp_cfm_set_all(tmp_cfm, czero)
970 CALL cp_cfm_create(tmp_cfm_1, v%matrix_struct)
971 CALL cp_cfm_set_all(tmp_cfm_1, czero)
972 ALLOCATE (tmp_cmat(northo, nstate))
973 DO idim = 1, dim2
974 weight = weights(idim)
975 arr_zii = matrix_zii(:, idim)
976 ! tmp_cfm = zij_0*V
977 CALL parallel_gemm("N", "N", nmo, nstate, nmo, cone, &
978 zij_0(idim), v, czero, tmp_cfm)
979 ! tmp_cfm = tmp_cfm*diag_zij_dagg
980 CALL cp_cfm_column_scale(tmp_cfm, conjg(arr_zii))
981 ! tmp_cfm_1 = tmp_cfm*U_dagg
982 CALL parallel_gemm("N", "C", nmo, nstate, nstate, cone, tmp_cfm, &
983 u, czero, tmp_cfm_1)
984 CALL cp_cfm_scale(weight, tmp_cfm_1)
985 ! zdiag = zdiag + tmp_cfm_1'
986 CALL cp_cfm_scale_and_add(cone, zdiag, cone, tmp_cfm_1)
987
988 ! tmp_cfm = zij_0_dagg*V
989 CALL parallel_gemm("C", "N", nmo, nstate, nmo, cone, &
990 zij_0(idim), v, czero, tmp_cfm)
991
992 ! tmp_cfm = tmp_cfm*diag_zij
993 CALL cp_cfm_column_scale(tmp_cfm, arr_zii)
994 ! tmp_cfm_1 = tmp_cfm*U_dagg
995 CALL parallel_gemm("N", "C", nmo, nstate, nstate, cone, tmp_cfm, &
996 u, czero, tmp_cfm_1)
997 CALL cp_cfm_scale(weight, tmp_cfm_1)
998 ! zdiag = zdiag + tmp_cfm_1'
999 CALL cp_cfm_scale_and_add(cone, zdiag, cone, tmp_cfm_1)
1000 END DO ! idim
1001 CALL cp_cfm_release(tmp_cfm)
1002 CALL cp_cfm_release(tmp_cfm_1)
1003 DEALLOCATE (tmp_cmat)
1004 ALLOCATE (tmp_cmat(northo, nextra))
1005 CALL cp_cfm_get_submatrix(zdiag, tmp_cmat, nocc + 1, nocc + 1, &
1006 northo, nextra, .false.)
1007 ! 'grad_ctilde'
1008 CALL cp_cfm_set_submatrix(grad_ctilde, tmp_cmat)
1009 DEALLOCATE (tmp_cmat)
1010 ! ctrans_lambda = c_tilde_dagg*grad_ctilde
1011 CALL parallel_gemm("C", "N", nextra, nextra, northo, cone, c_tilde, grad_ctilde, czero, ctrans_lambda)
1012 !WRITE(*,*) "norm(ctrans_lambda) = ", cp_cfm_norm(ctrans_lambda, "F")
1013 ! 'grad_ctilde = - c_tilde*ctrans_lambda + grad_ctilde'
1014 CALL parallel_gemm("N", "N", northo, nextra, nextra, -cone, c_tilde, ctrans_lambda, cone, grad_ctilde)
1015 END IF ! nextra > 0
1016
1017 ! tolerance
1018 IF (nextra > 0) THEN
1019 tolc = 0.0_dp
1020 CALL cp_fm_create(tmp_fm, grad_ctilde%matrix_struct)
1021 CALL cp_cfm_to_fm(grad_ctilde, tmp_fm)
1022 CALL cp_fm_maxabsval(tmp_fm, tolc)
1023 CALL cp_fm_release(tmp_fm)
1024 !WRITE(*,*) 'tolc = ', tolc
1025 tol = tol + tolc
1026 END IF
1027 !WRITE(*,*) 'tol = ', tol
1028
1029 IF (nextra > 0) THEN
1030 !WRITE(*,*) 'new_direction: ', new_direction
1031 IF (new_direction) THEN
1032 line_searches = line_searches + 1
1033 IF (mintol > tol) THEN
1034 mintol = tol
1035 miniter = iter
1036 END IF
1037
1038 IF (unit_nr > 0 .AND. modulo(iter, out_each) == 0) THEN
1039 sum_spread_ii = alpha*nstate/twopi/twopi - spread_sum
1040 avg_spread_ii = sum_spread_ii/nstate
1041 WRITE (unit_nr, '(T4,A,T26,A,T48,A)') &
1042 "Iteration", "Avg. Spread_ii", "Tolerance"
1043 WRITE (unit_nr, '(T4,I7,T20,F20.10,T45,E12.4)') &
1044 iter, avg_spread_ii, tol
1045 CALL m_flush(unit_nr)
1046 END IF
1047 IF (tol < eps_localization) EXIT
1048
1049 IF (do_cg) THEN
1050 cnorm2_gct = czero
1051 cnorm2_gct_cross = czero
1052 CALL cp_cfm_trace(grad_ctilde, gct_old, cnorm2_gct_cross)
1053 norm2_gct_cross = real(cnorm2_gct_cross, kind=dp)
1054 gct_old%local_data = grad_ctilde%local_data
1055 CALL cp_cfm_trace(grad_ctilde, gct_old, cnorm2_gct)
1056 norm2_gct = real(cnorm2_gct, kind=dp)
1057 ! compute beta_pr
1058 beta_pr = (norm2_gct - norm2_gct_cross)/norm2_old
1059 norm2_old = norm2_gct
1060 beta = max(0.0_dp, beta_pr)
1061 !WRITE(*,*) 'beta = ', beta
1062 ! compute skc / ska = beta * skc / ska + grad_ctilde / G
1063 CALL cp_cfm_scale(beta, skc)
1064 CALL cp_cfm_scale_and_add(cone, skc, cone, gct_old)
1065 CALL cp_cfm_trace(skc, gct_old, cnorm2_gct_cross)
1066 norm2_gct_cross = real(cnorm2_gct_cross, kind=dp)
1067 IF (norm2_gct_cross <= 0.0_dp) THEN ! back to steepest ascent
1068 CALL cp_cfm_scale_and_add(czero, skc, cone, gct_old)
1069 END IF
1070 ELSE
1071 CALL cp_cfm_scale_and_add(czero, skc, cone, grad_ctilde)
1072 END IF
1073 line_search_count = 0
1074 END IF
1075
1076 line_search_count = line_search_count + 1
1077 !WRITE(*,*) 'line_search_count = ', line_search_count
1078 energy(line_search_count) = spread_sum
1079
1080 ! gold line search
1081 new_direction = .false.
1082 IF (line_search_count == 1) THEN
1083 lsl = 1
1084 lsr = 0
1085 lsm = 1
1086 pos(1) = 0.0_dp
1087 pos(2) = ds_min/gold_sec
1088 ds = pos(2)
1089 ELSE
1090 IF (line_search_count == 50) THEN
1091 IF (abs(energy(line_search_count) - energy(line_search_count - 1)) < 1.0e-4_dp) THEN
1092 cpwarn("Line search failed to converge properly")
1093 ds_min = 0.1_dp
1094 new_direction = .true.
1095 ds = pos(line_search_count)
1096 line_search_count = 0
1097 ELSE
1098 cpabort("No. of line searches exceeds 50")
1099 END IF
1100 ELSE
1101 IF (lsr == 0) THEN
1102 IF (energy(line_search_count - 1) > energy(line_search_count)) THEN
1103 lsr = line_search_count
1104 pos(line_search_count + 1) = pos(lsm) + (pos(lsr) - pos(lsm))*gold_sec
1105 ELSE
1106 lsl = lsm
1107 lsm = line_search_count
1108 pos(line_search_count + 1) = pos(line_search_count)/gold_sec
1109 END IF
1110 ELSE
1111 IF (pos(line_search_count) < pos(lsm)) THEN
1112 IF (energy(line_search_count) > energy(lsm)) THEN
1113 lsr = lsm
1114 lsm = line_search_count
1115 ELSE
1116 lsl = line_search_count
1117 END IF
1118 ELSE
1119 IF (energy(line_search_count) > energy(lsm)) THEN
1120 lsl = lsm
1121 lsm = line_search_count
1122 ELSE
1123 lsr = line_search_count
1124 END IF
1125 END IF
1126 IF (pos(lsr) - pos(lsm) > pos(lsm) - pos(lsl)) THEN
1127 pos(line_search_count + 1) = pos(lsm) + gold_sec*(pos(lsr) - pos(lsm))
1128 ELSE
1129 pos(line_search_count + 1) = pos(lsl) + gold_sec*(pos(lsm) - pos(lsl))
1130 END IF
1131 IF ((pos(lsr) - pos(lsl)) < 1.0e-3_dp*pos(lsr)) THEN
1132 new_direction = .true.
1133 END IF
1134 END IF ! lsr .eq. 0
1135 END IF ! line_search_count .eq. 50
1136 ! now go to the suggested point
1137 ds = pos(line_search_count + 1) - pos(line_search_count)
1138 !WRITE(*,*) 'lsl, lsr, lsm, ds = ', lsl, lsr, lsm, ds
1139 IF ((abs(ds) < 1.0e-10_dp) .AND. (lsl == 1)) THEN
1140 new_direction = .true.
1141 ds_min = 0.5_dp/alpha
1142 ELSE IF (abs(ds) > 10.0_dp) THEN
1143 new_direction = .true.
1144 ds_min = 0.5_dp/alpha
1145 ELSE
1146 ds_min = pos(line_search_count + 1)
1147 END IF
1148 END IF ! first step
1149 ! 'c_tilde = c_tilde + d*skc'
1150 CALL cp_cfm_scale(ds, skc)
1151 CALL cp_cfm_scale_and_add(cone, c_tilde, cone, skc)
1152 ELSE
1153 IF (mintol > tol) THEN
1154 mintol = tol
1155 miniter = iter
1156 END IF
1157 IF (unit_nr > 0 .AND. modulo(iter, out_each) == 0) THEN
1158 sum_spread_ii = alpha*nstate/twopi/twopi - spread_sum
1159 avg_spread_ii = sum_spread_ii/nstate
1160 WRITE (unit_nr, '(T4,A,T26,A,T48,A)') &
1161 "Iteration", "Avg. Spread_ii", "Tolerance"
1162 WRITE (unit_nr, '(T4,I7,T20,F20.10,T45,E12.4)') &
1163 iter, avg_spread_ii, tol
1164 CALL m_flush(unit_nr)
1165 END IF
1166 IF (tol < eps_localization) EXIT
1167 END IF ! nextra > 0
1168
1169 END DO ! iteration
1170
1171 IF ((unit_nr > 0) .AND. (iter == max_iter)) THEN
1172 WRITE (unit_nr, '(T4,A,T4,A)') "Min. Itr.", "Min. Tol."
1173 WRITE (unit_nr, '(T4,I7,T4,E12.4)') miniter, mintol
1174 CALL m_flush(unit_nr)
1175 END IF
1176
1177 CALL cp_cfm_to_fm(u, matrix_u)
1178
1179 IF (nextra > 0) THEN
1180 rmat%local_data = real(v%local_data, kind=dp)
1181 CALL rotate_orbitals_edf(rmat, vectors_all, vectors)
1182
1183 CALL cp_cfm_release(c_tilde)
1184 CALL cp_cfm_release(grad_ctilde)
1185 CALL cp_cfm_release(gct_old)
1186 CALL cp_cfm_release(skc)
1187 CALL cp_cfm_release(ul)
1188 CALL cp_cfm_release(zdiag)
1189 CALL cp_cfm_release(ctrans_lambda)
1190 CALL cp_fm_release(id_nextra)
1191 CALL cp_fm_release(vectors_all)
1192 CALL cp_cfm_release(v)
1193 CALL cp_fm_release(matrix_v)
1194 CALL cp_fm_release(matrix_v_all)
1195 CALL cp_cfm_release(vl)
1196 DEALLOCATE (arr_zii)
1197 ELSE
1198 rmat%local_data = matrix_u%local_data
1199 CALL rotate_orbitals(rmat, vectors)
1200 END IF
1201 DO idim = 1, dim2
1202 CALL cp_cfm_release(zij_0(idim))
1203 END DO
1204 DEALLOCATE (zij_0)
1205
1206 DO idim = 1, dim2
1207 zij(1, idim)%local_data = real(c_zij(idim)%local_data, dp)
1208 zij(2, idim)%local_data = aimag(c_zij(idim)%local_data)
1209 CALL cp_cfm_release(c_zij(idim))
1210 END DO
1211 DEALLOCATE (c_zij)
1212 CALL cp_fm_release(rmat)
1213 CALL cp_cfm_release(u)
1214 CALL cp_fm_release(matrix_u)
1215 DEALLOCATE (matrix_zii, sum_spread)
1216
1217 CALL timestop(handle)
1218
1219 END SUBROUTINE jacobi_cg_edf_ls
1220
1221! **************************************************************************************************
1222!> \brief ...
1223!> \param vmatrix ...
1224! **************************************************************************************************
1225 SUBROUTINE ortho_vectors(vmatrix)
1226
1227 TYPE(cp_fm_type), INTENT(IN) :: vmatrix
1228
1229 CHARACTER(LEN=*), PARAMETER :: routinen = 'ortho_vectors'
1230
1231 INTEGER :: handle, n, ncol
1232 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
1233 TYPE(cp_fm_type) :: overlap_vv
1234
1235 CALL timeset(routinen, handle)
1236
1237 NULLIFY (fm_struct_tmp)
1238
1239 CALL cp_fm_get_info(matrix=vmatrix, nrow_global=n, ncol_global=ncol)
1240
1241 CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ncol, ncol_global=ncol, &
1242 para_env=vmatrix%matrix_struct%para_env, &
1243 context=vmatrix%matrix_struct%context)
1244 CALL cp_fm_create(overlap_vv, fm_struct_tmp, "overlap_vv")
1245 CALL cp_fm_struct_release(fm_struct_tmp)
1246
1247 CALL parallel_gemm('T', 'N', ncol, ncol, n, 1.0_dp, vmatrix, vmatrix, 0.0_dp, overlap_vv)
1248 CALL cp_fm_cholesky_decompose(overlap_vv)
1249 CALL cp_fm_triangular_multiply(overlap_vv, vmatrix, n_cols=ncol, side='R', invert_tr=.true.)
1250
1251 CALL cp_fm_release(overlap_vv)
1252
1253 CALL timestop(handle)
1254
1255 END SUBROUTINE ortho_vectors
1256
1257! **************************************************************************************************
1258!> \brief ...
1259!> \param istate ...
1260!> \param jstate ...
1261!> \param st ...
1262!> \param ct ...
1263!> \param zij ...
1264! **************************************************************************************************
1265 SUBROUTINE rotate_zij(istate, jstate, st, ct, zij)
1266 INTEGER, INTENT(IN) :: istate, jstate
1267 REAL(kind=dp), INTENT(IN) :: st, ct
1268 TYPE(cp_cfm_type) :: zij(:)
1269
1270 INTEGER :: id
1271
1272! Locals
1273
1274 DO id = 1, SIZE(zij, 1)
1275 CALL cp_cfm_rot_cols(zij(id), istate, jstate, ct, st)
1276 CALL cp_cfm_rot_rows(zij(id), istate, jstate, ct, st)
1277 END DO
1278
1279 END SUBROUTINE rotate_zij
1280! **************************************************************************************************
1281!> \brief ...
1282!> \param istate ...
1283!> \param jstate ...
1284!> \param st ...
1285!> \param ct ...
1286!> \param rmat ...
1287! **************************************************************************************************
1288 SUBROUTINE rotate_rmat(istate, jstate, st, ct, rmat)
1289 INTEGER, INTENT(IN) :: istate, jstate
1290 REAL(kind=dp), INTENT(IN) :: st, ct
1291 TYPE(cp_cfm_type), INTENT(IN) :: rmat
1292
1293 CALL cp_cfm_rot_cols(rmat, istate, jstate, ct, st)
1294
1295 END SUBROUTINE rotate_rmat
1296! **************************************************************************************************
1297!> \brief ...
1298!> \param mii ...
1299!> \param mjj ...
1300!> \param mij ...
1301!> \param weights ...
1302!> \param theta ...
1303!> \param grad_ij ...
1304!> \param step ...
1305! **************************************************************************************************
1306 SUBROUTINE get_angle(mii, mjj, mij, weights, theta, grad_ij, step)
1307 COMPLEX(KIND=dp), POINTER :: mii(:), mjj(:), mij(:)
1308 REAL(kind=dp), INTENT(IN) :: weights(:)
1309 REAL(kind=dp), INTENT(OUT) :: theta
1310 REAL(kind=dp), INTENT(IN), OPTIONAL :: grad_ij, step
1311
1312 COMPLEX(KIND=dp) :: z11, z12, z22
1313 INTEGER :: dim_m, idim
1314 REAL(kind=dp) :: a12, b12, d2, ratio
1315
1316 a12 = 0.0_dp
1317 b12 = 0.0_dp
1318 dim_m = SIZE(mii)
1319 DO idim = 1, dim_m
1320 z11 = mii(idim)
1321 z22 = mjj(idim)
1322 z12 = mij(idim)
1323 a12 = a12 + weights(idim)*real(conjg(z12)*(z11 - z22), kind=dp)
1324 b12 = b12 + weights(idim)*real((z12*conjg(z12) - &
1325 0.25_dp*(z11 - z22)*(conjg(z11) - conjg(z22))), kind=dp)
1326 END DO
1327 IF (abs(b12) > 1.e-10_dp) THEN
1328 ratio = -a12/b12
1329 theta = 0.25_dp*atan(ratio)
1330 ELSE IF (abs(b12) < 1.e-10_dp) THEN
1331 b12 = 0.0_dp
1332 theta = 0.0_dp
1333 ELSE
1334 theta = 0.25_dp*pi
1335 END IF
1336 IF (PRESENT(grad_ij)) theta = theta + step*grad_ij
1337! Check second derivative info
1338 d2 = a12*sin(4._dp*theta) - b12*cos(4._dp*theta)
1339 IF (d2 <= 0._dp) THEN ! go to the maximum, not the minimum
1340 IF (theta > 0.0_dp) THEN ! make theta as small as possible
1341 theta = theta - 0.25_dp*pi
1342 ELSE
1343 theta = theta + 0.25_dp*pi
1344 END IF
1345 END IF
1346 END SUBROUTINE get_angle
1347! **************************************************************************************************
1348!> \brief ...
1349!> \param zij ...
1350!> \param weights ...
1351!> \param tolerance ...
1352!> \param grad ...
1353! **************************************************************************************************
1354 SUBROUTINE check_tolerance(zij, weights, tolerance, grad)
1355 TYPE(cp_cfm_type) :: zij(:)
1356 REAL(kind=dp), INTENT(IN) :: weights(:)
1357 REAL(kind=dp), INTENT(OUT) :: tolerance
1358 TYPE(cp_fm_type), INTENT(OUT), OPTIONAL :: grad
1359
1360 CHARACTER(len=*), PARAMETER :: routinen = 'check_tolerance'
1361
1362 INTEGER :: handle
1363 TYPE(cp_fm_type) :: force
1364
1365 CALL timeset(routinen, handle)
1366
1367! compute gradient at t=0
1368
1369 CALL cp_fm_create(force, zij(1)%matrix_struct)
1370 CALL cp_fm_set_all(force, 0._dp)
1371 CALL grad_at_0(zij, weights, force)
1372 CALL cp_fm_maxabsval(force, tolerance)
1373 IF (PRESENT(grad)) CALL cp_fm_to_fm(force, grad)
1374 CALL cp_fm_release(force)
1375
1376 CALL timestop(handle)
1377
1378 END SUBROUTINE check_tolerance
1379
1380! **************************************************************************************************
1381!> \brief ...
1382!> \param rmat ...
1383!> \param vectors ...
1384! **************************************************************************************************
1385 SUBROUTINE rotate_orbitals(rmat, vectors)
1386 TYPE(cp_fm_type), INTENT(IN) :: rmat, vectors
1387
1388 INTEGER :: k, n
1389 TYPE(cp_fm_type) :: wf
1390
1391 CALL cp_fm_create(wf, vectors%matrix_struct)
1392 CALL cp_fm_get_info(vectors, nrow_global=n, ncol_global=k)
1393 CALL parallel_gemm("N", "N", n, k, k, 1.0_dp, vectors, rmat, 0.0_dp, wf)
1394 CALL cp_fm_to_fm(wf, vectors)
1395 CALL cp_fm_release(wf)
1396 END SUBROUTINE rotate_orbitals
1397
1398! **************************************************************************************************
1399!> \brief ...
1400!> \param rmat ...
1401!> \param vectors ...
1402! **************************************************************************************************
1403 SUBROUTINE rotate_orbitals_cfm(rmat, vectors)
1404 TYPE(cp_cfm_type), INTENT(IN) :: rmat, vectors
1405
1406 INTEGER :: k, n
1407 TYPE(cp_cfm_type) :: wf
1408
1409 CALL cp_cfm_create(wf, vectors%matrix_struct)
1410 CALL cp_cfm_get_info(vectors, nrow_global=n, ncol_global=k)
1411 CALL parallel_gemm("N", "N", n, k, k, z_one, vectors, rmat, z_zero, wf)
1412 CALL cp_cfm_to_cfm(wf, vectors)
1413 CALL cp_cfm_release(wf)
1414 END SUBROUTINE rotate_orbitals_cfm
1415
1416! **************************************************************************************************
1417!> \brief ...
1418!> \param rmat ...
1419!> \param vec_all ...
1420!> \param vectors ...
1421! **************************************************************************************************
1422 SUBROUTINE rotate_orbitals_edf(rmat, vec_all, vectors)
1423 TYPE(cp_fm_type), INTENT(IN) :: rmat, vec_all, vectors
1424
1425 INTEGER :: k, l, n
1426 TYPE(cp_fm_type) :: wf
1427
1428 CALL cp_fm_create(wf, vectors%matrix_struct)
1429 CALL cp_fm_get_info(vec_all, nrow_global=n, ncol_global=k)
1430 CALL cp_fm_get_info(rmat, ncol_global=l)
1431
1432 CALL parallel_gemm("N", "N", n, l, k, 1.0_dp, vec_all, rmat, 0.0_dp, wf)
1433 CALL cp_fm_to_fm(wf, vectors)
1434 CALL cp_fm_release(wf)
1435 END SUBROUTINE rotate_orbitals_edf
1436! **************************************************************************************************
1437!> \brief ...
1438!> \param diag ...
1439!> \param weights ...
1440!> \param matrix ...
1441!> \param ndim ...
1442! **************************************************************************************************
1443 SUBROUTINE gradsq_at_0(diag, weights, matrix, ndim)
1444 COMPLEX(KIND=dp), DIMENSION(:, :), POINTER :: diag
1445 REAL(kind=dp), INTENT(IN) :: weights(:)
1446 TYPE(cp_fm_type), INTENT(IN) :: matrix
1447 INTEGER, INTENT(IN) :: ndim
1448
1449 COMPLEX(KIND=dp) :: zii, zjj
1450 INTEGER :: idim, istate, jstate, ncol_local, &
1451 nrow_global, nrow_local
1452 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1453 REAL(kind=dp) :: gradsq_ij
1454
1455 CALL cp_fm_get_info(matrix, nrow_local=nrow_local, &
1456 ncol_local=ncol_local, nrow_global=nrow_global, &
1457 row_indices=row_indices, col_indices=col_indices)
1458
1459 DO istate = 1, nrow_local
1460 DO jstate = 1, ncol_local
1461! get real and imaginary parts
1462 gradsq_ij = 0.0_dp
1463 DO idim = 1, ndim
1464 zii = diag(row_indices(istate), idim)
1465 zjj = diag(col_indices(jstate), idim)
1466 gradsq_ij = gradsq_ij + weights(idim)* &
1467 4.0_dp*real((conjg(zii)*zii + conjg(zjj)*zjj), kind=dp)
1468 END DO
1469 matrix%local_data(istate, jstate) = gradsq_ij
1470 END DO
1471 END DO
1472 END SUBROUTINE gradsq_at_0
1473! **************************************************************************************************
1474!> \brief ...
1475!> \param matrix_p ...
1476!> \param weights ...
1477!> \param matrix ...
1478! **************************************************************************************************
1479 SUBROUTINE grad_at_0(matrix_p, weights, matrix)
1480 TYPE(cp_cfm_type) :: matrix_p(:)
1481 REAL(kind=dp), INTENT(IN) :: weights(:)
1482 TYPE(cp_fm_type), INTENT(IN) :: matrix
1483
1484 COMPLEX(KIND=dp) :: zii, zij, zjj
1485 COMPLEX(KIND=dp), DIMENSION(:, :), POINTER :: diag
1486 INTEGER :: dim_m, idim, istate, jstate, ncol_local, &
1487 nrow_global, nrow_local
1488 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1489 REAL(kind=dp) :: grad_ij
1490
1491 NULLIFY (diag)
1492 CALL cp_fm_get_info(matrix, nrow_local=nrow_local, &
1493 ncol_local=ncol_local, nrow_global=nrow_global, &
1494 row_indices=row_indices, col_indices=col_indices)
1495 dim_m = SIZE(matrix_p, 1)
1496 ALLOCATE (diag(nrow_global, dim_m))
1497
1498 DO idim = 1, dim_m
1499 DO istate = 1, nrow_global
1500 CALL cp_cfm_get_element(matrix_p(idim), istate, istate, diag(istate, idim))
1501 END DO
1502 END DO
1503
1504 DO istate = 1, nrow_local
1505 DO jstate = 1, ncol_local
1506! get real and imaginary parts
1507 grad_ij = 0.0_dp
1508 DO idim = 1, dim_m
1509 zii = diag(row_indices(istate), idim)
1510 zjj = diag(col_indices(jstate), idim)
1511 zij = matrix_p(idim)%local_data(istate, jstate)
1512 grad_ij = grad_ij + weights(idim)* &
1513 REAL(4.0_dp*conjg(zij)*(zjj - zii), dp)
1514 END DO
1515 matrix%local_data(istate, jstate) = grad_ij
1516 END DO
1517 END DO
1518 DEALLOCATE (diag)
1519 END SUBROUTINE grad_at_0
1520
1521! return energy and maximum gradient in the current point
1522! **************************************************************************************************
1523!> \brief ...
1524!> \param weights ...
1525!> \param zij ...
1526!> \param tolerance ...
1527!> \param value ...
1528! **************************************************************************************************
1529 SUBROUTINE check_tolerance_new(weights, zij, tolerance, value)
1530 REAL(kind=dp), INTENT(IN) :: weights(:)
1531 TYPE(cp_fm_type), INTENT(IN) :: zij(:, :)
1532 REAL(kind=dp) :: tolerance, value
1533
1534 COMPLEX(KIND=dp) :: kii, kij, kjj
1535 COMPLEX(KIND=dp), DIMENSION(:, :), POINTER :: diag
1536 INTEGER :: idim, istate, jstate, ncol_local, &
1537 nrow_global, nrow_local
1538 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1539 REAL(kind=dp) :: grad_ij, ra, rb
1540
1541 NULLIFY (diag)
1542 CALL cp_fm_get_info(zij(1, 1), nrow_local=nrow_local, &
1543 ncol_local=ncol_local, nrow_global=nrow_global, &
1544 row_indices=row_indices, col_indices=col_indices)
1545 ALLOCATE (diag(nrow_global, SIZE(zij, 2)))
1546 value = 0.0_dp
1547 DO idim = 1, SIZE(zij, 2)
1548 DO istate = 1, nrow_global
1549 CALL cp_fm_get_element(zij(1, idim), istate, istate, ra)
1550 CALL cp_fm_get_element(zij(2, idim), istate, istate, rb)
1551 diag(istate, idim) = cmplx(ra, rb, dp)
1552 value = value + weights(idim) - weights(idim)*abs(diag(istate, idim))**2
1553 END DO
1554 END DO
1555 tolerance = 0.0_dp
1556 DO istate = 1, nrow_local
1557 DO jstate = 1, ncol_local
1558 grad_ij = 0.0_dp
1559 DO idim = 1, SIZE(zij, 2)
1560 kii = diag(row_indices(istate), idim)
1561 kjj = diag(col_indices(jstate), idim)
1562 ra = zij(1, idim)%local_data(istate, jstate)
1563 rb = zij(2, idim)%local_data(istate, jstate)
1564 kij = cmplx(ra, rb, dp)
1565 grad_ij = grad_ij + weights(idim)* &
1566 REAL(4.0_dp*conjg(kij)*(kjj - kii), dp)
1567 END DO
1568 tolerance = max(abs(grad_ij), tolerance)
1569 END DO
1570 END DO
1571 CALL zij(1, 1)%matrix_struct%para_env%max(tolerance)
1572
1573 DEALLOCATE (diag)
1574
1575 END SUBROUTINE check_tolerance_new
1576
1577! **************************************************************************************************
1578!> \brief yet another crazy try, computes the angles needed to rotate the orbitals first
1579!> and rotates them all at the same time (hoping for the best of course)
1580!> \param weights ...
1581!> \param zij ...
1582!> \param vectors ...
1583!> \param max_iter ...
1584!> \param max_crazy_angle ...
1585!> \param crazy_scale ...
1586!> \param crazy_use_diag ...
1587!> \param eps_localization ...
1588!> \param iterations ...
1589!> \param converged ...
1590! **************************************************************************************************
1591 SUBROUTINE crazy_rotations(weights, zij, vectors, max_iter, max_crazy_angle, crazy_scale, crazy_use_diag, &
1592 eps_localization, iterations, converged)
1593 REAL(kind=dp), INTENT(IN) :: weights(:)
1594 TYPE(cp_fm_type), INTENT(IN) :: zij(:, :), vectors
1595 INTEGER, INTENT(IN) :: max_iter
1596 REAL(kind=dp), INTENT(IN) :: max_crazy_angle
1597 REAL(kind=dp) :: crazy_scale
1598 LOGICAL :: crazy_use_diag
1599 REAL(kind=dp), INTENT(IN) :: eps_localization
1600 INTEGER :: iterations
1601 LOGICAL, INTENT(out), OPTIONAL :: converged
1602
1603 CHARACTER(len=*), PARAMETER :: routinen = 'crazy_rotations'
1604 COMPLEX(KIND=dp), PARAMETER :: cone = (1.0_dp, 0.0_dp), &
1605 czero = (0.0_dp, 0.0_dp)
1606
1607 COMPLEX(KIND=dp), DIMENSION(:), POINTER :: evals_exp
1608 COMPLEX(KIND=dp), DIMENSION(:, :), POINTER :: diag_z
1609 COMPLEX(KIND=dp), POINTER :: mii(:), mij(:), mjj(:)
1610 INTEGER :: dim2, handle, i, icol, idim, irow, &
1611 method, ncol_global, ncol_local, &
1612 norder, nrow_global, nrow_local, &
1613 nsquare, unit_nr
1614 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1615 LOGICAL :: do_emd
1616 REAL(kind=dp) :: eps_exp, limit_crazy_angle, maxeval, &
1617 norm, ra, rb, theta, tolerance, value
1618 REAL(kind=dp), DIMENSION(:), POINTER :: evals
1619 TYPE(cp_cfm_type) :: cmat_a, cmat_r, cmat_t1
1620 TYPE(cp_fm_type) :: mat_r, mat_t, mat_theta, mat_u
1621
1622 CALL timeset(routinen, handle)
1623 NULLIFY (row_indices, col_indices)
1624 CALL cp_fm_get_info(zij(1, 1), nrow_global=nrow_global, &
1625 ncol_global=ncol_global, &
1626 row_indices=row_indices, col_indices=col_indices, &
1627 nrow_local=nrow_local, ncol_local=ncol_local)
1628
1629 limit_crazy_angle = max_crazy_angle
1630
1631 NULLIFY (diag_z, evals, evals_exp, mii, mij, mjj)
1632 dim2 = SIZE(zij, 2)
1633 ALLOCATE (diag_z(nrow_global, dim2))
1634 ALLOCATE (evals(nrow_global))
1635 ALLOCATE (evals_exp(nrow_global))
1636
1637 CALL cp_cfm_create(cmat_a, zij(1, 1)%matrix_struct)
1638 CALL cp_cfm_create(cmat_r, zij(1, 1)%matrix_struct)
1639 CALL cp_cfm_create(cmat_t1, zij(1, 1)%matrix_struct)
1640
1641 CALL cp_fm_create(mat_u, zij(1, 1)%matrix_struct)
1642 CALL cp_fm_create(mat_t, zij(1, 1)%matrix_struct)
1643 CALL cp_fm_create(mat_r, zij(1, 1)%matrix_struct)
1644
1645 CALL cp_fm_create(mat_theta, zij(1, 1)%matrix_struct)
1646
1647 CALL cp_fm_set_all(mat_r, 0.0_dp, 1.0_dp)
1648 CALL cp_fm_set_all(mat_t, 0.0_dp)
1649 ALLOCATE (mii(dim2), mij(dim2), mjj(dim2))
1650 DO idim = 1, dim2
1651 CALL cp_fm_scale_and_add(1.0_dp, mat_t, weights(idim), zij(1, idim))
1652 CALL cp_fm_scale_and_add(1.0_dp, mat_t, weights(idim), zij(2, idim))
1653 END DO
1654 CALL cp_fm_syevd(mat_t, mat_u, evals)
1655 DO idim = 1, dim2
1656 ! rotate z's
1657 CALL parallel_gemm('N', 'N', nrow_global, nrow_global, nrow_global, 1.0_dp, zij(1, idim), mat_u, 0.0_dp, mat_t)
1658 CALL parallel_gemm('T', 'N', nrow_global, nrow_global, nrow_global, 1.0_dp, mat_u, mat_t, 0.0_dp, zij(1, idim))
1659 CALL parallel_gemm('N', 'N', nrow_global, nrow_global, nrow_global, 1.0_dp, zij(2, idim), mat_u, 0.0_dp, mat_t)
1660 CALL parallel_gemm('T', 'N', nrow_global, nrow_global, nrow_global, 1.0_dp, mat_u, mat_t, 0.0_dp, zij(2, idim))
1661 END DO
1662 ! collect rotation matrix
1663 CALL parallel_gemm('N', 'N', nrow_global, nrow_global, nrow_global, 1.0_dp, mat_r, mat_u, 0.0_dp, mat_t)
1664 CALL cp_fm_to_fm(mat_t, mat_r)
1665
1666 unit_nr = -1
1667 IF (cmat_a%matrix_struct%para_env%is_source()) THEN
1668 unit_nr = cp_logger_get_default_unit_nr()
1669 WRITE (unit_nr, '(T2,A7,A6,1X,A20,A12,A12,A12)') &
1670 "CRAZY| ", "Iter", "value ", "gradient", "Max. eval", "limit"
1671 END IF
1672
1673 iterations = 0
1674 tolerance = 1.0_dp
1675
1676 DO
1677 iterations = iterations + 1
1678 DO idim = 1, dim2
1679 DO i = 1, nrow_global
1680 CALL cp_fm_get_element(zij(1, idim), i, i, ra)
1681 CALL cp_fm_get_element(zij(2, idim), i, i, rb)
1682 diag_z(i, idim) = cmplx(ra, rb, dp)
1683 END DO
1684 END DO
1685 DO irow = 1, nrow_local
1686 DO icol = 1, ncol_local
1687 DO idim = 1, dim2
1688 ra = zij(1, idim)%local_data(irow, icol)
1689 rb = zij(2, idim)%local_data(irow, icol)
1690 mij(idim) = cmplx(ra, rb, dp)
1691 mii(idim) = diag_z(row_indices(irow), idim)
1692 mjj(idim) = diag_z(col_indices(icol), idim)
1693 END DO
1694 IF (row_indices(irow) /= col_indices(icol)) THEN
1695 CALL get_angle(mii, mjj, mij, weights, theta)
1696 theta = crazy_scale*theta
1697 IF (theta > limit_crazy_angle) theta = limit_crazy_angle
1698 IF (theta < -limit_crazy_angle) theta = -limit_crazy_angle
1699 IF (crazy_use_diag) THEN
1700 cmat_a%local_data(irow, icol) = -cmplx(0.0_dp, theta, dp)
1701 ELSE
1702 mat_theta%local_data(irow, icol) = -theta
1703 END IF
1704 ELSE
1705 IF (crazy_use_diag) THEN
1706 cmat_a%local_data(irow, icol) = czero
1707 ELSE
1708 mat_theta%local_data(irow, icol) = 0.0_dp
1709 END IF
1710 END IF
1711 END DO
1712 END DO
1713
1714 ! construct rotation matrix U based on A using diagonalization
1715 ! alternatively, exp based on repeated squaring could be faster
1716 IF (crazy_use_diag) THEN
1717 CALL cp_cfm_heevd(cmat_a, cmat_r, evals)
1718 maxeval = maxval(abs(evals))
1719 evals_exp(:) = exp((0.0_dp, -1.0_dp)*evals(:))
1720 CALL cp_cfm_to_cfm(cmat_r, cmat_t1)
1721 CALL cp_cfm_column_scale(cmat_t1, evals_exp)
1722 CALL parallel_gemm('N', 'C', nrow_global, nrow_global, nrow_global, cone, &
1723 cmat_t1, cmat_r, czero, cmat_a)
1724 mat_u%local_data = real(cmat_a%local_data, kind=dp) ! U is a real matrix
1725 ELSE
1726 do_emd = .false.
1727 method = 2
1728 eps_exp = 1.0_dp*epsilon(eps_exp)
1729 CALL cp_fm_maxabsrownorm(mat_theta, norm)
1730 maxeval = norm ! an upper bound
1731 CALL get_nsquare_norder(norm, nsquare, norder, eps_exp, method, do_emd)
1732 CALL exp_pade_real(mat_u, mat_theta, nsquare, norder)
1733 END IF
1734
1735 DO idim = 1, dim2
1736 ! rotate z's
1737 CALL parallel_gemm('N', 'N', nrow_global, nrow_global, nrow_global, 1.0_dp, zij(1, idim), mat_u, 0.0_dp, mat_t)
1738 CALL parallel_gemm('T', 'N', nrow_global, nrow_global, nrow_global, 1.0_dp, mat_u, mat_t, 0.0_dp, zij(1, idim))
1739 CALL parallel_gemm('N', 'N', nrow_global, nrow_global, nrow_global, 1.0_dp, zij(2, idim), mat_u, 0.0_dp, mat_t)
1740 CALL parallel_gemm('T', 'N', nrow_global, nrow_global, nrow_global, 1.0_dp, mat_u, mat_t, 0.0_dp, zij(2, idim))
1741 END DO
1742 ! collect rotation matrix
1743 CALL parallel_gemm('N', 'N', nrow_global, nrow_global, nrow_global, 1.0_dp, mat_r, mat_u, 0.0_dp, mat_t)
1744 CALL cp_fm_to_fm(mat_t, mat_r)
1745
1746 CALL check_tolerance_new(weights, zij, tolerance, value)
1747
1748 IF (unit_nr > 0) THEN
1749 WRITE (unit_nr, '(T2,A7,I6,1X,G20.15,E12.4,E12.4,E12.4)') &
1750 "CRAZY| ", iterations, value, tolerance, maxeval, limit_crazy_angle
1751 CALL m_flush(unit_nr)
1752 END IF
1753 IF (tolerance < eps_localization .OR. iterations >= max_iter) EXIT
1754 END DO
1755
1756 IF (PRESENT(converged)) converged = (tolerance < eps_localization)
1757
1758 CALL cp_cfm_release(cmat_a)
1759 CALL cp_cfm_release(cmat_r)
1760 CALL cp_cfm_release(cmat_t1)
1761
1762 CALL cp_fm_release(mat_u)
1763 CALL cp_fm_release(mat_t)
1764 CALL cp_fm_release(mat_theta)
1765
1766 CALL rotate_orbitals(mat_r, vectors)
1767
1768 CALL cp_fm_release(mat_r)
1769 DEALLOCATE (evals_exp, evals, diag_z)
1770 DEALLOCATE (mii, mij, mjj)
1771
1772 CALL timestop(handle)
1773
1774 END SUBROUTINE crazy_rotations
1775
1776! **************************************************************************************************
1777!> \brief use the exponential parametrization as described in to perform a direct mini
1778!> Gerd Berghold et al. PRB 61 (15), pag. 10040 (2000)
1779!> none of the input is modified for the time being, just finds the rotations
1780!> that minimizes, and throws it away afterwards :-)
1781!> apart from being expensive and not cleaned, this works fine
1782!> useful to try different spread functionals
1783!> \param weights ...
1784!> \param zij ...
1785!> \param vectors ...
1786!> \param max_iter ...
1787!> \param eps_localization ...
1788!> \param iterations ...
1789! **************************************************************************************************
1790 SUBROUTINE direct_mini(weights, zij, vectors, max_iter, eps_localization, iterations)
1791 REAL(kind=dp), INTENT(IN) :: weights(:)
1792 TYPE(cp_fm_type), INTENT(IN) :: zij(:, :), vectors
1793 INTEGER, INTENT(IN) :: max_iter
1794 REAL(kind=dp), INTENT(IN) :: eps_localization
1795 INTEGER :: iterations
1796
1797 CHARACTER(len=*), PARAMETER :: routinen = 'direct_mini'
1798 COMPLEX(KIND=dp), PARAMETER :: cone = (1.0_dp, 0.0_dp), &
1799 czero = (0.0_dp, 0.0_dp)
1800 REAL(kind=dp), PARAMETER :: gold_sec = 0.3819_dp
1801
1802 COMPLEX(KIND=dp) :: lk, ll, tmp
1803 COMPLEX(KIND=dp), DIMENSION(:), POINTER :: evals_exp
1804 COMPLEX(KIND=dp), DIMENSION(:, :), POINTER :: diag_z
1805 INTEGER :: handle, i, icol, idim, irow, &
1806 line_search_count, line_searches, lsl, &
1807 lsm, lsr, n, ncol_local, ndim, &
1808 nrow_local, output_unit
1809 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1810 LOGICAL :: new_direction
1811 REAL(kind=dp) :: a, b, beta_pr, c, denom, ds, ds_min, fa, &
1812 fb, fc, nom, normg, normg_cross, &
1813 normg_old, npos, omega, tol, val, x0, &
1814 x1, xa, xb, xc
1815 REAL(kind=dp), DIMENSION(150) :: energy, grad, pos
1816 REAL(kind=dp), DIMENSION(:), POINTER :: evals, fval, fvald
1817 TYPE(cp_cfm_type) :: cmat_a, cmat_b, cmat_m, cmat_r, cmat_t1, &
1818 cmat_t2, cmat_u
1819 TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: c_zij
1820 TYPE(cp_fm_type) :: matrix_a, matrix_g, matrix_g_old, &
1821 matrix_g_search, matrix_h, matrix_r, &
1822 matrix_t
1823
1824 NULLIFY (evals, evals_exp, diag_z, fval, fvald)
1825
1826 CALL timeset(routinen, handle)
1827 output_unit = cp_logger_get_default_io_unit()
1828
1829 n = zij(1, 1)%matrix_struct%nrow_global
1830 ndim = (SIZE(zij, 2))
1831
1832 IF (output_unit > 0) THEN
1833 WRITE (output_unit, '(T4,A )') "Localization by direct minimization of the functional; "
1834 WRITE (output_unit, '(T5,2A13,A20,A20,A10 )') " Line search ", " Iteration ", " Functional ", " Tolerance ", " ds Min "
1835 END IF
1836
1837 ALLOCATE (evals(n), evals_exp(n), diag_z(n, ndim), fval(n), fvald(n))
1838 ALLOCATE (c_zij(ndim))
1839
1840 ! create the three complex matrices Z
1841 DO idim = 1, ndim
1842 CALL cp_cfm_create(c_zij(idim), zij(1, 1)%matrix_struct)
1843 c_zij(idim)%local_data = cmplx(zij(1, idim)%local_data, &
1844 zij(2, idim)%local_data, dp)
1845 END DO
1846
1847 CALL cp_fm_create(matrix_a, zij(1, 1)%matrix_struct)
1848 CALL cp_fm_create(matrix_g, zij(1, 1)%matrix_struct)
1849 CALL cp_fm_create(matrix_t, zij(1, 1)%matrix_struct)
1850 CALL cp_fm_create(matrix_h, zij(1, 1)%matrix_struct)
1851 CALL cp_fm_create(matrix_g_search, zij(1, 1)%matrix_struct)
1852 CALL cp_fm_create(matrix_g_old, zij(1, 1)%matrix_struct)
1853 CALL cp_fm_create(matrix_r, zij(1, 1)%matrix_struct)
1854 CALL cp_fm_set_all(matrix_r, 0.0_dp, 1.0_dp)
1855
1856 CALL cp_fm_set_all(matrix_a, 0.0_dp)
1857! CALL cp_fm_init_random ( matrix_A )
1858
1859 CALL cp_cfm_create(cmat_a, zij(1, 1)%matrix_struct)
1860 CALL cp_cfm_create(cmat_u, zij(1, 1)%matrix_struct)
1861 CALL cp_cfm_create(cmat_r, zij(1, 1)%matrix_struct)
1862 CALL cp_cfm_create(cmat_t1, zij(1, 1)%matrix_struct)
1863 CALL cp_cfm_create(cmat_t2, zij(1, 1)%matrix_struct)
1864 CALL cp_cfm_create(cmat_b, zij(1, 1)%matrix_struct)
1865 CALL cp_cfm_create(cmat_m, zij(1, 1)%matrix_struct)
1866
1867 CALL cp_cfm_get_info(cmat_b, nrow_local=nrow_local, ncol_local=ncol_local, &
1868 row_indices=row_indices, col_indices=col_indices)
1869
1870 CALL cp_fm_set_all(matrix_g_old, 0.0_dp)
1871 CALL cp_fm_set_all(matrix_g_search, 0.0_dp)
1872 normg_old = 1.0e30_dp
1873 ds_min = 1.0_dp
1874 new_direction = .true.
1875 iterations = 0
1876 line_searches = 0
1877 line_search_count = 0
1878 DO
1879 iterations = iterations + 1
1880 ! compute U,R,evals given A
1881 cmat_a%local_data = cmplx(0.0_dp, matrix_a%local_data, dp) ! cmat_A is hermitian, evals are reals
1882 CALL cp_cfm_heevd(cmat_a, cmat_r, evals)
1883 evals_exp(:) = exp((0.0_dp, -1.0_dp)*evals(:))
1884 CALL cp_cfm_to_cfm(cmat_r, cmat_t1)
1885 CALL cp_cfm_column_scale(cmat_t1, evals_exp)
1886 CALL parallel_gemm('N', 'C', n, n, n, cone, cmat_t1, cmat_r, czero, cmat_u)
1887 cmat_u%local_data = real(cmat_u%local_data, kind=dp) ! enforce numerics, U is a real matrix
1888
1889 IF (new_direction .AND. mod(line_searches, 20) == 5) THEN ! reset with A .eq. 0
1890 DO idim = 1, ndim
1891 CALL parallel_gemm('N', 'N', n, n, n, cone, c_zij(idim), cmat_u, czero, cmat_t1)
1892 CALL parallel_gemm('C', 'N', n, n, n, cone, cmat_u, cmat_t1, czero, c_zij(idim))
1893 END DO
1894 ! collect rotation matrix
1895 matrix_h%local_data = real(cmat_u%local_data, kind=dp)
1896 CALL parallel_gemm('N', 'N', n, n, n, 1.0_dp, matrix_r, matrix_h, 0.0_dp, matrix_t)
1897 CALL cp_fm_to_fm(matrix_t, matrix_r)
1898
1899 CALL cp_cfm_set_all(cmat_u, czero, cone)
1900 CALL cp_cfm_set_all(cmat_r, czero, cone)
1901 CALL cp_cfm_set_all(cmat_a, czero)
1902 CALL cp_fm_set_all(matrix_a, 0.0_dp)
1903 evals(:) = 0.0_dp
1904 evals_exp(:) = exp((0.0_dp, -1.0_dp)*evals(:))
1905 CALL cp_fm_set_all(matrix_g_old, 0.0_dp)
1906 CALL cp_fm_set_all(matrix_g_search, 0.0_dp)
1907 normg_old = 1.0e30_dp
1908 END IF
1909
1910 ! compute Omega and M
1911 CALL cp_cfm_set_all(cmat_m, czero)
1912 omega = 0.0_dp
1913 DO idim = 1, ndim
1914 CALL parallel_gemm('N', 'N', n, n, n, cone, c_zij(idim), cmat_u, czero, cmat_t1) ! t1=ZU
1915 CALL parallel_gemm('C', 'N', n, n, n, cone, cmat_u, cmat_t1, czero, cmat_t2) ! t2=(U^T)ZU
1916 DO i = 1, n
1917 CALL cp_cfm_get_element(cmat_t2, i, i, diag_z(i, idim))
1918 SELECT CASE (2) ! allows for selection of different spread functionals
1919 CASE (1)
1920 fval(i) = -weights(idim)*log(abs(diag_z(i, idim))**2)
1921 fvald(i) = -weights(idim)/(abs(diag_z(i, idim))**2)
1922 CASE (2) ! corresponds to the jacobi setup
1923 fval(i) = weights(idim) - weights(idim)*abs(diag_z(i, idim))**2
1924 fvald(i) = -weights(idim)
1925 END SELECT
1926 omega = omega + fval(i)
1927 END DO
1928 DO icol = 1, ncol_local
1929 DO irow = 1, nrow_local
1930 tmp = cmat_t1%local_data(irow, icol)*conjg(diag_z(col_indices(icol), idim))
1931 cmat_m%local_data(irow, icol) = cmat_m%local_data(irow, icol) &
1932 + 4.0_dp*fvald(col_indices(icol))*real(tmp, kind=dp)
1933 END DO
1934 END DO
1935 END DO
1936
1937 ! compute Hessian diagonal approximation for the preconditioner
1938 IF (.true.) THEN
1939 CALL gradsq_at_0(diag_z, weights, matrix_h, ndim)
1940 ELSE
1941 CALL cp_fm_set_all(matrix_h, 1.0_dp)
1942 END IF
1943
1944 ! compute B
1945 DO icol = 1, ncol_local
1946 DO irow = 1, nrow_local
1947 ll = (0.0_dp, -1.0_dp)*evals(row_indices(irow))
1948 lk = (0.0_dp, -1.0_dp)*evals(col_indices(icol))
1949 IF (abs(ll - lk) < 0.5_dp) THEN ! use a series expansion to avoid loss of precision
1950 tmp = 1.0_dp
1951 cmat_b%local_data(irow, icol) = 0.0_dp
1952 DO i = 1, 16
1953 cmat_b%local_data(irow, icol) = cmat_b%local_data(irow, icol) + tmp
1954 tmp = tmp*(ll - lk)/(i + 1)
1955 END DO
1956 cmat_b%local_data(irow, icol) = cmat_b%local_data(irow, icol)*exp(lk)
1957 ELSE
1958 cmat_b%local_data(irow, icol) = (exp(lk) - exp(ll))/(lk - ll)
1959 END IF
1960 END DO
1961 END DO
1962 ! compute gradient matrix_G
1963
1964 CALL parallel_gemm('C', 'N', n, n, n, cone, cmat_m, cmat_r, czero, cmat_t1) ! t1=(M^T)(R^T)
1965 CALL parallel_gemm('C', 'N', n, n, n, cone, cmat_r, cmat_t1, czero, cmat_t2) ! t2=(R)t1
1966 CALL cp_cfm_schur_product(cmat_t2, cmat_b, cmat_t1)
1967 CALL parallel_gemm('N', 'C', n, n, n, cone, cmat_t1, cmat_r, czero, cmat_t2)
1968 CALL parallel_gemm('N', 'N', n, n, n, cone, cmat_r, cmat_t2, czero, cmat_t1)
1969 matrix_g%local_data = real(cmat_t1%local_data, kind=dp)
1970 CALL cp_fm_transpose(matrix_g, matrix_t)
1971 CALL cp_fm_scale_and_add(-1.0_dp, matrix_g, 1.0_dp, matrix_t)
1972 CALL cp_fm_maxabsval(matrix_g, tol)
1973
1974 ! from here on, minimizing technology
1975 IF (new_direction) THEN
1976 ! energy converged up to machine precision ?
1977 line_searches = line_searches + 1
1978 IF (output_unit > 0) THEN
1979 WRITE (output_unit, '(T5,I10,T18,I10,T31,2F20.6,F10.3)') line_searches, iterations, omega, tol, ds_min
1980 CALL m_flush(output_unit)
1981 END IF
1982 IF (tol < eps_localization .OR. iterations > max_iter) EXIT
1983
1984 IF (.true.) THEN ! do conjugate gradient CG
1985 CALL cp_fm_trace(matrix_g, matrix_g_old, normg_cross)
1986 normg_cross = normg_cross*0.5_dp ! takes into account the fact that A is antisymmetric
1987 ! apply the preconditioner
1988 DO icol = 1, ncol_local
1989 DO irow = 1, nrow_local
1990 matrix_g_old%local_data(irow, icol) = matrix_g%local_data(irow, icol)/matrix_h%local_data(irow, icol)
1991 END DO
1992 END DO
1993 CALL cp_fm_trace(matrix_g, matrix_g_old, normg)
1994 normg = normg*0.5_dp
1995 beta_pr = (normg - normg_cross)/normg_old
1996 normg_old = normg
1997 beta_pr = max(beta_pr, 0.0_dp)
1998 CALL cp_fm_scale_and_add(beta_pr, matrix_g_search, -1.0_dp, matrix_g_old)
1999 CALL cp_fm_trace(matrix_g_search, matrix_g_old, normg_cross)
2000 IF (normg_cross >= 0) THEN ! back to SD
2001 IF (matrix_a%matrix_struct%para_env%is_source()) THEN
2002 WRITE (cp_logger_get_default_unit_nr(), *) "!"
2003 END IF
2004 beta_pr = 0.0_dp
2005 CALL cp_fm_scale_and_add(beta_pr, matrix_g_search, -1.0_dp, matrix_g_old)
2006 END IF
2007 ELSE ! SD
2008 CALL cp_fm_scale_and_add(0.0_dp, matrix_g_search, -1.0_dp, matrix_g)
2009 END IF
2010 ! ds_min=1.0E-4_dp
2011 line_search_count = 0
2012 END IF
2013 line_search_count = line_search_count + 1
2014 energy(line_search_count) = omega
2015
2016 ! line search section
2017 SELECT CASE (3)
2018 CASE (1) ! two point line search
2019 SELECT CASE (line_search_count)
2020 CASE (1)
2021 pos(1) = 0.0_dp
2022 pos(2) = ds_min
2023 CALL cp_fm_trace(matrix_g, matrix_g_search, grad(1))
2024 grad(1) = grad(1)/2.0_dp
2025 new_direction = .false.
2026 CASE (2)
2027 new_direction = .true.
2028 x0 = pos(1) ! 0.0_dp
2029 c = energy(1)
2030 b = grad(1)
2031 x1 = pos(2)
2032 a = (energy(2) - b*x1 - c)/(x1**2)
2033 IF (a <= 0.0_dp) a = 1.0e-15_dp
2034 npos = -b/(2.0_dp*a)
2035 val = a*npos**2 + b*npos + c
2036 IF (val < energy(1) .AND. val <= energy(2)) THEN
2037 ! we go to a minimum, but ...
2038 ! we take a guard against too large steps
2039 pos(3) = min(npos, maxval(pos(1:2))*4.0_dp)
2040 ELSE ! just take an extended step
2041 pos(3) = maxval(pos(1:2))*2.0_dp
2042 END IF
2043 END SELECT
2044 CASE (2) ! 3 point line search
2045 SELECT CASE (line_search_count)
2046 CASE (1)
2047 new_direction = .false.
2048 pos(1) = 0.0_dp
2049 pos(2) = ds_min*0.8_dp
2050 CASE (2)
2051 new_direction = .false.
2052 IF (energy(2) > energy(1)) THEN
2053 pos(3) = ds_min*0.7_dp
2054 ELSE
2055 pos(3) = ds_min*1.4_dp
2056 END IF
2057 CASE (3)
2058 new_direction = .true.
2059 xa = pos(1)
2060 xb = pos(2)
2061 xc = pos(3)
2062 fa = energy(1)
2063 fb = energy(2)
2064 fc = energy(3)
2065 nom = (xb - xa)**2*(fb - fc) - (xb - xc)**2*(fb - fa)
2066 denom = (xb - xa)*(fb - fc) - (xb - xc)*(fb - fa)
2067 IF (abs(denom) <= 1.0e-18_dp*max(abs(fb - fc), abs(fb - fa))) THEN
2068 npos = xb
2069 ELSE
2070 npos = xb - 0.5_dp*nom/denom ! position of the stationary point
2071 END IF
2072 val = (npos - xa)*(npos - xb)*fc/((xc - xa)*(xc - xb)) + &
2073 (npos - xb)*(npos - xc)*fa/((xa - xb)*(xa - xc)) + &
2074 (npos - xc)*(npos - xa)*fb/((xb - xc)*(xb - xa))
2075 IF (val < fa .AND. val <= fb .AND. val <= fc) THEN ! OK, we go to a minimum
2076 ! we take a guard against too large steps
2077 pos(4) = max(maxval(pos(1:3))*0.01_dp, &
2078 min(npos, maxval(pos(1:3))*4.0_dp))
2079 ELSE ! just take an extended step
2080 pos(4) = maxval(pos(1:3))*2.0_dp
2081 END IF
2082 END SELECT
2083 CASE (3) ! golden section hunt
2084 new_direction = .false.
2085 IF (line_search_count == 1) THEN
2086 lsl = 1
2087 lsr = 0
2088 lsm = 1
2089 pos(1) = 0.0_dp
2090 pos(2) = ds_min/gold_sec
2091 ELSE
2092 IF (line_search_count == 150) cpabort("Too many")
2093 IF (lsr == 0) THEN
2094 IF (energy(line_search_count - 1) < energy(line_search_count)) THEN
2095 lsr = line_search_count
2096 pos(line_search_count + 1) = pos(lsm) + (pos(lsr) - pos(lsm))*gold_sec
2097 ELSE
2098 lsl = lsm
2099 lsm = line_search_count
2100 pos(line_search_count + 1) = pos(line_search_count)/gold_sec
2101 END IF
2102 ELSE
2103 IF (pos(line_search_count) < pos(lsm)) THEN
2104 IF (energy(line_search_count) < energy(lsm)) THEN
2105 lsr = lsm
2106 lsm = line_search_count
2107 ELSE
2108 lsl = line_search_count
2109 END IF
2110 ELSE
2111 IF (energy(line_search_count) < energy(lsm)) THEN
2112 lsl = lsm
2113 lsm = line_search_count
2114 ELSE
2115 lsr = line_search_count
2116 END IF
2117 END IF
2118 IF (pos(lsr) - pos(lsm) > pos(lsm) - pos(lsl)) THEN
2119 pos(line_search_count + 1) = pos(lsm) + gold_sec*(pos(lsr) - pos(lsm))
2120 ELSE
2121 pos(line_search_count + 1) = pos(lsl) + gold_sec*(pos(lsm) - pos(lsl))
2122 END IF
2123 IF ((pos(lsr) - pos(lsl)) < 1.0e-3_dp*pos(lsr)) THEN
2124 new_direction = .true.
2125 END IF
2126 END IF ! lsr .eq. 0
2127 END IF ! first step
2128 END SELECT
2129 ! now go to the suggested point
2130 ds_min = pos(line_search_count + 1)
2131 ds = pos(line_search_count + 1) - pos(line_search_count)
2132 CALL cp_fm_scale_and_add(1.0_dp, matrix_a, ds, matrix_g_search)
2133 END DO
2134
2135 ! collect rotation matrix
2136 matrix_h%local_data = real(cmat_u%local_data, kind=dp)
2137 CALL parallel_gemm('N', 'N', n, n, n, 1.0_dp, matrix_r, matrix_h, 0.0_dp, matrix_t)
2138 CALL cp_fm_to_fm(matrix_t, matrix_r)
2139 CALL rotate_orbitals(matrix_r, vectors)
2140 CALL cp_fm_release(matrix_r)
2141
2142 CALL cp_fm_release(matrix_a)
2143 CALL cp_fm_release(matrix_g)
2144 CALL cp_fm_release(matrix_h)
2145 CALL cp_fm_release(matrix_t)
2146 CALL cp_fm_release(matrix_g_search)
2147 CALL cp_fm_release(matrix_g_old)
2148 CALL cp_cfm_release(cmat_a)
2149 CALL cp_cfm_release(cmat_u)
2150 CALL cp_cfm_release(cmat_r)
2151 CALL cp_cfm_release(cmat_t1)
2152 CALL cp_cfm_release(cmat_t2)
2153 CALL cp_cfm_release(cmat_b)
2154 CALL cp_cfm_release(cmat_m)
2155
2156 DEALLOCATE (evals, evals_exp, fval, fvald)
2157
2158 DO idim = 1, SIZE(c_zij)
2159 zij(1, idim)%local_data = real(c_zij(idim)%local_data, dp)
2160 zij(2, idim)%local_data = aimag(c_zij(idim)%local_data)
2161 CALL cp_cfm_release(c_zij(idim))
2162 END DO
2163 DEALLOCATE (c_zij)
2164 DEALLOCATE (diag_z)
2165
2166 CALL timestop(handle)
2167
2168 END SUBROUTINE direct_mini
2169
2170! **************************************************************************************************
2171!> \brief Parallel algorithm for jacobi rotations
2172!> \param weights ...
2173!> \param zij ...
2174!> \param vectors ...
2175!> \param para_env ...
2176!> \param max_iter ...
2177!> \param eps_localization ...
2178!> \param sweeps ...
2179!> \param out_each ...
2180!> \param target_time ...
2181!> \param start_time ...
2182!> \param restricted ...
2183!> \par History
2184!> use allgather for improved performance
2185!> \author MI (11.2009)
2186! **************************************************************************************************
2187 SUBROUTINE jacobi_rot_para(weights, zij, vectors, para_env, max_iter, eps_localization, &
2188 sweeps, out_each, target_time, start_time, restricted)
2189
2190 REAL(kind=dp), INTENT(IN) :: weights(:)
2191 TYPE(cp_fm_type), INTENT(IN) :: zij(:, :), vectors
2192 TYPE(mp_para_env_type), POINTER :: para_env
2193 INTEGER, INTENT(IN) :: max_iter
2194 REAL(kind=dp), INTENT(IN) :: eps_localization
2195 INTEGER :: sweeps
2196 INTEGER, INTENT(IN) :: out_each
2197 REAL(dp) :: target_time, start_time
2198 INTEGER :: restricted
2199
2200 CHARACTER(len=*), PARAMETER :: routinen = 'jacobi_rot_para'
2201
2202 INTEGER :: dim2, handle, i, idim, ii, ilow1, ip, j, &
2203 nblock, nblock_max, ns_me, nstate, &
2204 output_unit
2205 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: ns_bound
2206 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: rotmat, z_ij_loc_im, z_ij_loc_re
2207 REAL(kind=dp) :: xstate
2208 TYPE(cp_fm_type) :: rmat
2209 TYPE(set_c_2d_type), DIMENSION(:), POINTER :: cz_ij_loc
2210
2211 CALL timeset(routinen, handle)
2212
2213 output_unit = cp_logger_get_default_io_unit()
2214
2215 NULLIFY (cz_ij_loc)
2216
2217 dim2 = SIZE(zij, 2)
2218
2219 CALL cp_fm_create(rmat, zij(1, 1)%matrix_struct)
2220 CALL cp_fm_set_all(rmat, 0._dp, 1._dp)
2221
2222 CALL cp_fm_get_info(rmat, nrow_global=nstate)
2223
2224 IF (restricted > 0) THEN
2225 IF (output_unit > 0) THEN
2226 WRITE (output_unit, '(T4,A,I2,A )') "JACOBI: for the ROKS method, the last ", restricted, " orbitals DO NOT ROTATE"
2227 END IF
2228 nstate = nstate - restricted
2229 END IF
2230
2231 ! Distribution of the states (XXXXX safe against more pe than states ??? XXXXX)
2232 xstate = real(nstate, dp)/real(para_env%num_pe, dp)
2233 ALLOCATE (ns_bound(0:para_env%num_pe - 1, 2))
2234 DO ip = 1, para_env%num_pe
2235 ns_bound(ip - 1, 1) = min(nstate, nint(xstate*(ip - 1))) + 1
2236 ns_bound(ip - 1, 2) = min(nstate, nint(xstate*ip))
2237 END DO
2238 nblock_max = 0
2239 DO ip = 0, para_env%num_pe - 1
2240 nblock = ns_bound(ip, 2) - ns_bound(ip, 1) + 1
2241 nblock_max = max(nblock_max, nblock)
2242 END DO
2243
2244 ! otbtain local part of the matrix (could be made faster, but is likely irrelevant).
2245 ALLOCATE (z_ij_loc_re(nstate, nblock_max))
2246 ALLOCATE (z_ij_loc_im(nstate, nblock_max))
2247 ALLOCATE (cz_ij_loc(dim2))
2248 DO idim = 1, dim2
2249 DO ip = 0, para_env%num_pe - 1
2250 nblock = ns_bound(ip, 2) - ns_bound(ip, 1) + 1
2251 CALL cp_fm_get_submatrix(zij(1, idim), z_ij_loc_re, 1, ns_bound(ip, 1), nstate, nblock)
2252 CALL cp_fm_get_submatrix(zij(2, idim), z_ij_loc_im, 1, ns_bound(ip, 1), nstate, nblock)
2253 IF (para_env%mepos == ip) THEN
2254 ALLOCATE (cz_ij_loc(idim)%c_array(nstate, nblock))
2255 DO i = 1, nblock
2256 DO j = 1, nstate
2257 cz_ij_loc(idim)%c_array(j, i) = cmplx(z_ij_loc_re(j, i), z_ij_loc_im(j, i), dp)
2258 END DO
2259 END DO
2260 END IF
2261 END DO ! ip
2262 END DO
2263 DEALLOCATE (z_ij_loc_re)
2264 DEALLOCATE (z_ij_loc_im)
2265
2266 ALLOCATE (rotmat(nstate, 2*nblock_max))
2267
2268 CALL jacobi_rot_para_core(weights, para_env, max_iter, sweeps, out_each, dim2, nstate, nblock_max, ns_bound, &
2269 cz_ij_loc, rotmat, output_unit, eps_localization=eps_localization, &
2270 target_time=target_time, start_time=start_time)
2271
2272 ilow1 = ns_bound(para_env%mepos, 1)
2273 ns_me = ns_bound(para_env%mepos, 2) - ns_bound(para_env%mepos, 1) + 1
2274 ALLOCATE (z_ij_loc_re(nstate, nblock_max))
2275 ALLOCATE (z_ij_loc_im(nstate, nblock_max))
2276 DO idim = 1, dim2
2277 DO ip = 0, para_env%num_pe - 1
2278 z_ij_loc_re = 0.0_dp
2279 z_ij_loc_im = 0.0_dp
2280 nblock = ns_bound(ip, 2) - ns_bound(ip, 1) + 1
2281 IF (ip == para_env%mepos) THEN
2282 ns_me = nblock
2283 DO i = 1, ns_me
2284 ii = ilow1 + i - 1
2285 DO j = 1, nstate
2286 z_ij_loc_re(j, i) = real(cz_ij_loc(idim)%c_array(j, i), dp)
2287 z_ij_loc_im(j, i) = aimag(cz_ij_loc(idim)%c_array(j, i))
2288 END DO
2289 END DO
2290 END IF
2291 CALL para_env%bcast(z_ij_loc_re, ip)
2292 CALL para_env%bcast(z_ij_loc_im, ip)
2293 CALL cp_fm_set_submatrix(zij(1, idim), z_ij_loc_re, 1, ns_bound(ip, 1), nstate, nblock)
2294 CALL cp_fm_set_submatrix(zij(2, idim), z_ij_loc_im, 1, ns_bound(ip, 1), nstate, nblock)
2295 END DO ! ip
2296 END DO
2297
2298 DO ip = 0, para_env%num_pe - 1
2299 z_ij_loc_re = 0.0_dp
2300 nblock = ns_bound(ip, 2) - ns_bound(ip, 1) + 1
2301 IF (ip == para_env%mepos) THEN
2302 ns_me = nblock
2303 DO i = 1, ns_me
2304 ii = ilow1 + i - 1
2305 DO j = 1, nstate
2306 z_ij_loc_re(j, i) = rotmat(j, i)
2307 END DO
2308 END DO
2309 END IF
2310 CALL para_env%bcast(z_ij_loc_re, ip)
2311 CALL cp_fm_set_submatrix(rmat, z_ij_loc_re, 1, ns_bound(ip, 1), nstate, nblock)
2312 END DO
2313
2314 DEALLOCATE (z_ij_loc_re)
2315 DEALLOCATE (z_ij_loc_im)
2316 DO idim = 1, dim2
2317 DEALLOCATE (cz_ij_loc(idim)%c_array)
2318 END DO
2319 DEALLOCATE (cz_ij_loc)
2320
2321 CALL para_env%sync()
2322 CALL rotate_orbitals(rmat, vectors)
2323 CALL cp_fm_release(rmat)
2324
2325 DEALLOCATE (rotmat)
2326 DEALLOCATE (ns_bound)
2327
2328 CALL timestop(handle)
2329
2330 END SUBROUTINE jacobi_rot_para
2331
2332! **************************************************************************************************
2333!> \brief almost identical to 'jacobi_rot_para' but with different inout variables
2334!> \param weights ...
2335!> \param czij ...
2336!> \param para_env ...
2337!> \param max_iter ...
2338!> \param rmat ...
2339!> \param tol_out ...
2340!> \author Soumya Ghosh (08/21)
2341! **************************************************************************************************
2342 SUBROUTINE jacobi_rot_para_1(weights, czij, para_env, max_iter, rmat, tol_out)
2343
2344 REAL(kind=dp), INTENT(IN) :: weights(:)
2345 TYPE(cp_cfm_type), INTENT(IN) :: czij(:)
2346 TYPE(mp_para_env_type), POINTER :: para_env
2347 INTEGER, INTENT(IN) :: max_iter
2348 TYPE(cp_cfm_type), INTENT(IN) :: rmat
2349 REAL(dp), INTENT(OUT), OPTIONAL :: tol_out
2350
2351 CHARACTER(len=*), PARAMETER :: routinen = 'jacobi_rot_para_1'
2352
2353 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: czij_array
2354 INTEGER :: dim2, handle, i, idim, ii, ilow1, ip, j, &
2355 nblock, nblock_max, ns_me, nstate, &
2356 sweeps
2357 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: ns_bound
2358 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: rotmat, z_ij_loc_re
2359 REAL(kind=dp) :: xstate
2360 TYPE(set_c_2d_type), DIMENSION(:), POINTER :: cz_ij_loc
2361
2362 CALL timeset(routinen, handle)
2363
2364 dim2 = SIZE(czij)
2365
2366 CALL cp_cfm_set_all(rmat, cmplx(0._dp, 0._dp, dp), cmplx(1._dp, 0._dp, dp))
2367
2368 CALL cp_cfm_get_info(rmat, nrow_global=nstate)
2369
2370 ! Distribution of the states (XXXXX safe against more pe than states ??? XXXXX)
2371 xstate = real(nstate, dp)/real(para_env%num_pe, dp)
2372 ALLOCATE (ns_bound(0:para_env%num_pe - 1, 2))
2373 DO ip = 1, para_env%num_pe
2374 ns_bound(ip - 1, 1) = min(nstate, nint(xstate*(ip - 1))) + 1
2375 ns_bound(ip - 1, 2) = min(nstate, nint(xstate*ip))
2376 END DO
2377 nblock_max = 0
2378 DO ip = 0, para_env%num_pe - 1
2379 nblock = ns_bound(ip, 2) - ns_bound(ip, 1) + 1
2380 nblock_max = max(nblock_max, nblock)
2381 END DO
2382
2383 ! otbtain local part of the matrix (could be made faster, but is likely irrelevant).
2384 ALLOCATE (czij_array(nstate, nblock_max))
2385 ALLOCATE (cz_ij_loc(dim2))
2386 DO idim = 1, dim2
2387 DO ip = 0, para_env%num_pe - 1
2388 nblock = ns_bound(ip, 2) - ns_bound(ip, 1) + 1
2389 ! cfm --> allocatable
2390 CALL cp_cfm_get_submatrix(czij(idim), czij_array, 1, ns_bound(ip, 1), nstate, nblock)
2391 IF (para_env%mepos == ip) THEN
2392 ns_me = nblock
2393 ALLOCATE (cz_ij_loc(idim)%c_array(nstate, ns_me))
2394 DO i = 1, ns_me
2395 DO j = 1, nstate
2396 cz_ij_loc(idim)%c_array(j, i) = czij_array(j, i)
2397 END DO
2398 END DO
2399 END IF
2400 END DO ! ip
2401 END DO
2402 DEALLOCATE (czij_array)
2403
2404 ALLOCATE (rotmat(nstate, 2*nblock_max))
2405
2406 CALL jacobi_rot_para_core(weights, para_env, max_iter, sweeps, 1, dim2, nstate, nblock_max, ns_bound, &
2407 cz_ij_loc, rotmat, 0, tol_out=tol_out)
2408
2409 ilow1 = ns_bound(para_env%mepos, 1)
2410 ns_me = ns_bound(para_env%mepos, 2) - ns_bound(para_env%mepos, 1) + 1
2411 ALLOCATE (z_ij_loc_re(nstate, nblock_max))
2412
2413 DO ip = 0, para_env%num_pe - 1
2414 z_ij_loc_re = 0.0_dp
2415 nblock = ns_bound(ip, 2) - ns_bound(ip, 1) + 1
2416 IF (ip == para_env%mepos) THEN
2417 ns_me = nblock
2418 DO i = 1, ns_me
2419 ii = ilow1 + i - 1
2420 DO j = 1, nstate
2421 z_ij_loc_re(j, i) = rotmat(j, i)
2422 END DO
2423 END DO
2424 END IF
2425 CALL para_env%bcast(z_ij_loc_re, ip)
2426 CALL cp_cfm_set_submatrix(rmat, cmplx(z_ij_loc_re, 0.0_dp, dp), 1, ns_bound(ip, 1), nstate, nblock)
2427 END DO
2428
2429 DEALLOCATE (z_ij_loc_re)
2430 DO idim = 1, dim2
2431 DEALLOCATE (cz_ij_loc(idim)%c_array)
2432 END DO
2433 DEALLOCATE (cz_ij_loc)
2434
2435 CALL para_env%sync()
2436
2437 DEALLOCATE (rotmat)
2438 DEALLOCATE (ns_bound)
2439
2440 CALL timestop(handle)
2441
2442 END SUBROUTINE jacobi_rot_para_1
2443
2444! **************************************************************************************************
2445!> \brief Parallel algorithm for jacobi rotations
2446!> \param weights ...
2447!> \param para_env ...
2448!> \param max_iter ...
2449!> \param sweeps ...
2450!> \param out_each ...
2451!> \param dim2 ...
2452!> \param nstate ...
2453!> \param nblock_max ...
2454!> \param ns_bound ...
2455!> \param cz_ij_loc ...
2456!> \param rotmat ...
2457!> \param output_unit ...
2458!> \param tol_out ...
2459!> \param eps_localization ...
2460!> \param target_time ...
2461!> \param start_time ...
2462!> \par History
2463!> split out to reuse with different input types
2464!> \author HF (05.2022)
2465! **************************************************************************************************
2466 SUBROUTINE jacobi_rot_para_core(weights, para_env, max_iter, sweeps, out_each, dim2, nstate, nblock_max, &
2467 ns_bound, cz_ij_loc, rotmat, output_unit, tol_out, eps_localization, target_time, start_time)
2468
2469 REAL(kind=dp), INTENT(IN) :: weights(:)
2470 TYPE(mp_para_env_type), POINTER :: para_env
2471 INTEGER, INTENT(IN) :: max_iter
2472 INTEGER, INTENT(OUT) :: sweeps
2473 INTEGER, INTENT(IN) :: out_each, dim2, nstate, nblock_max
2474 INTEGER, DIMENSION(0:, :), INTENT(IN) :: ns_bound
2475 TYPE(set_c_2d_type), DIMENSION(:), POINTER :: cz_ij_loc
2476 REAL(dp), DIMENSION(:, :), INTENT(OUT) :: rotmat
2477 INTEGER, INTENT(IN) :: output_unit
2478 REAL(dp), INTENT(OUT), OPTIONAL :: tol_out
2479 REAL(kind=dp), INTENT(IN), OPTIONAL :: eps_localization
2480 REAL(dp), OPTIONAL :: target_time, start_time
2481
2482 COMPLEX(KIND=dp) :: zi, zj
2483 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: c_array_me, c_array_partner
2484 COMPLEX(KIND=dp), POINTER :: mii(:), mij(:), mjj(:)
2485 INTEGER :: i, idim, ii, ik, il1, il2, il_recv, il_recv_partner, ilow1, ilow2, ip, ip_has_i, &
2486 ip_partner, ip_recv_from, ip_recv_partner, ipair, iperm, istate, iu1, iu2, iup1, iup2, j, &
2487 jj, jstate, k, kk, lsweep, n1, n2, npair, nperm, ns_me, ns_partner, ns_recv_from, &
2488 ns_recv_partner
2489 INTEGER, ALLOCATABLE, DIMENSION(:) :: rcount, rdispl
2490 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: list_pair
2491 LOGICAL :: should_stop
2492 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: gmat, rmat_loc, rmat_recv, rmat_send
2493 REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: rmat_recv_all
2494 REAL(kind=dp) :: ct, func, gmax, grad, ri, rj, st, t1, &
2495 t2, theta, tolerance, zc, zr
2496 TYPE(set_c_1d_type), DIMENSION(:), POINTER :: zdiag_all, zdiag_me
2497 TYPE(set_c_2d_type), DIMENSION(:), POINTER :: xyz_mix, xyz_mix_ns
2498
2499 NULLIFY (zdiag_all, zdiag_me)
2500 NULLIFY (xyz_mix, xyz_mix_ns)
2501 NULLIFY (mii, mij, mjj)
2502
2503 ALLOCATE (mii(dim2), mij(dim2), mjj(dim2))
2504
2505 ALLOCATE (rcount(para_env%num_pe))
2506 ALLOCATE (rdispl(para_env%num_pe))
2507
2508 tolerance = 1.0e10_dp
2509 sweeps = 0
2510
2511 ! number of processor pairs and number of permutations
2512 npair = (para_env%num_pe + 1)/2
2513 nperm = para_env%num_pe - mod(para_env%num_pe + 1, 2)
2514 ALLOCATE (list_pair(2, npair))
2515
2516 ! initialize rotation matrix
2517 rotmat = 0.0_dp
2518 DO i = ns_bound(para_env%mepos, 1), ns_bound(para_env%mepos, 2)
2519 ii = i - ns_bound(para_env%mepos, 1) + 1
2520 rotmat(i, ii) = 1.0_dp
2521 END DO
2522
2523 ALLOCATE (xyz_mix(dim2))
2524 ALLOCATE (xyz_mix_ns(dim2))
2525 ALLOCATE (zdiag_me(dim2))
2526 ALLOCATE (zdiag_all(dim2))
2527
2528 ns_me = ns_bound(para_env%mepos, 2) - ns_bound(para_env%mepos, 1) + 1
2529 IF (ns_me /= 0) THEN
2530 ALLOCATE (c_array_me(nstate, ns_me, dim2))
2531 DO idim = 1, dim2
2532 ALLOCATE (xyz_mix_ns(idim)%c_array(nstate, ns_me))
2533 END DO
2534 ALLOCATE (gmat(nstate, ns_me))
2535 END IF
2536
2537 DO idim = 1, dim2
2538 ALLOCATE (zdiag_me(idim)%c_array(nblock_max))
2539 zdiag_me(idim)%c_array = z_zero
2540 ALLOCATE (zdiag_all(idim)%c_array(para_env%num_pe*nblock_max))
2541 zdiag_all(idim)%c_array = z_zero
2542 END DO
2543 ALLOCATE (rmat_recv(nblock_max*2, nblock_max))
2544 ALLOCATE (rmat_send(nblock_max*2, nblock_max))
2545
2546 ! buffer for message passing
2547 ALLOCATE (rmat_recv_all(nblock_max*2, nblock_max, 0:para_env%num_pe - 1))
2548
2549 IF (output_unit > 0) THEN
2550 WRITE (output_unit, '(T4,A )') " Localization by iterative distributed Jacobi rotation"
2551 WRITE (output_unit, '(T20,A12,T32, A22,T60, A12,A8 )') "Iteration", "Functional", "Tolerance", " Time "
2552 END IF
2553
2554 DO lsweep = 1, max_iter + 1
2555 sweeps = lsweep
2556 IF (sweeps == max_iter + 1) THEN
2557 IF (output_unit > 0) THEN
2558 WRITE (output_unit, *) ' LOCALIZATION! loop did not converge within the maximum number of iterations.'
2559 WRITE (output_unit, *) ' Present Max. gradient = ', tolerance
2560 END IF
2561 EXIT
2562 END IF
2563 t1 = m_walltime()
2564
2565 DO iperm = 1, nperm
2566
2567 ! fix partners for this permutation, and get the number of states
2568 CALL eberlein(iperm, para_env, list_pair)
2569 ip_partner = -1
2570 ns_partner = 0
2571 DO ipair = 1, npair
2572 IF (list_pair(1, ipair) == para_env%mepos) THEN
2573 ip_partner = list_pair(2, ipair)
2574 EXIT
2575 ELSE IF (list_pair(2, ipair) == para_env%mepos) THEN
2576 ip_partner = list_pair(1, ipair)
2577 EXIT
2578 END IF
2579 END DO
2580 IF (ip_partner >= 0) THEN
2581 ns_partner = ns_bound(ip_partner, 2) - ns_bound(ip_partner, 1) + 1
2582 ELSE
2583 ns_partner = 0
2584 END IF
2585
2586 ! if there is a non-zero block connecting two partners, jacobi-sweep it.
2587 IF (ns_partner*ns_me /= 0) THEN
2588
2589 ALLOCATE (rmat_loc(ns_me + ns_partner, ns_me + ns_partner))
2590 rmat_loc = 0.0_dp
2591 DO i = 1, ns_me + ns_partner
2592 rmat_loc(i, i) = 1.0_dp
2593 END DO
2594
2595 ALLOCATE (c_array_partner(nstate, ns_partner, dim2))
2596
2597 DO idim = 1, dim2
2598 ALLOCATE (xyz_mix(idim)%c_array(ns_me + ns_partner, ns_me + ns_partner))
2599 DO i = 1, ns_me
2600 c_array_me(1:nstate, i, idim) = cz_ij_loc(idim)%c_array(1:nstate, i)
2601 END DO
2602 END DO
2603
2604 CALL para_env%sendrecv(msgin=c_array_me, dest=ip_partner, &
2605 msgout=c_array_partner, source=ip_partner)
2606
2607 n1 = ns_me
2608 n2 = ns_partner
2609 ilow1 = ns_bound(para_env%mepos, 1)
2610 iup1 = ns_bound(para_env%mepos, 1) + n1 - 1
2611 ilow2 = ns_bound(ip_partner, 1)
2612 iup2 = ns_bound(ip_partner, 1) + n2 - 1
2613 IF (ns_bound(para_env%mepos, 1) < ns_bound(ip_partner, 1)) THEN
2614 il1 = 1
2615 iu1 = n1
2616 iu1 = n1
2617 il2 = 1 + n1
2618 iu2 = n1 + n2
2619 ELSE
2620 il1 = 1 + n2
2621 iu1 = n1 + n2
2622 iu1 = n1 + n2
2623 il2 = 1
2624 iu2 = n2
2625 END IF
2626
2627 DO idim = 1, dim2
2628 DO i = 1, n1
2629 xyz_mix(idim)%c_array(il1:iu1, il1 + i - 1) = c_array_me(ilow1:iup1, i, idim)
2630 xyz_mix(idim)%c_array(il2:iu2, il1 + i - 1) = c_array_me(ilow2:iup2, i, idim)
2631 END DO
2632 DO i = 1, n2
2633 xyz_mix(idim)%c_array(il2:iu2, il2 + i - 1) = c_array_partner(ilow2:iup2, i, idim)
2634 xyz_mix(idim)%c_array(il1:iu1, il2 + i - 1) = c_array_partner(ilow1:iup1, i, idim)
2635 END DO
2636 END DO
2637
2638 DO istate = 1, n1 + n2
2639 DO jstate = istate + 1, n1 + n2
2640 DO idim = 1, dim2
2641 mii(idim) = xyz_mix(idim)%c_array(istate, istate)
2642 mij(idim) = xyz_mix(idim)%c_array(istate, jstate)
2643 mjj(idim) = xyz_mix(idim)%c_array(jstate, jstate)
2644 END DO
2645 CALL get_angle(mii, mjj, mij, weights, theta)
2646 st = sin(theta)
2647 ct = cos(theta)
2648 DO idim = 1, dim2
2649 DO i = 1, n1 + n2
2650 zi = ct*xyz_mix(idim)%c_array(i, istate) + st*xyz_mix(idim)%c_array(i, jstate)
2651 zj = -st*xyz_mix(idim)%c_array(i, istate) + ct*xyz_mix(idim)%c_array(i, jstate)
2652 xyz_mix(idim)%c_array(i, istate) = zi
2653 xyz_mix(idim)%c_array(i, jstate) = zj
2654 END DO
2655 DO i = 1, n1 + n2
2656 zi = ct*xyz_mix(idim)%c_array(istate, i) + st*xyz_mix(idim)%c_array(jstate, i)
2657 zj = -st*xyz_mix(idim)%c_array(istate, i) + ct*xyz_mix(idim)%c_array(jstate, i)
2658 xyz_mix(idim)%c_array(istate, i) = zi
2659 xyz_mix(idim)%c_array(jstate, i) = zj
2660 END DO
2661 END DO
2662
2663 DO i = 1, n1 + n2
2664 ri = ct*rmat_loc(i, istate) + st*rmat_loc(i, jstate)
2665 rj = ct*rmat_loc(i, jstate) - st*rmat_loc(i, istate)
2666 rmat_loc(i, istate) = ri
2667 rmat_loc(i, jstate) = rj
2668 END DO
2669 END DO
2670 END DO
2671
2672 k = nblock_max + 1
2673 CALL para_env%sendrecv(rotmat(1:nstate, 1:ns_me), ip_partner, &
2674 rotmat(1:nstate, k:k + n2 - 1), ip_partner)
2675
2676 IF (ilow1 < ilow2) THEN
2677 ! no longer compiles in official sdgb:
2678 ! probably inefficient:
2679 CALL dgemm("N", "N", nstate, n1, n2, 1.0_dp, rotmat(1:, k:), nstate, rmat_loc(1 + n1:, 1:n1), &
2680 n2, 0.0_dp, gmat(:, :), nstate)
2681 CALL dgemm("N", "N", nstate, n1, n1, 1.0_dp, rotmat(1:, 1:), nstate, rmat_loc(1:, 1:), &
2682 n1 + n2, 1.0_dp, gmat(:, :), nstate)
2683 ELSE
2684 CALL dgemm("N", "N", nstate, n1, n2, 1.0_dp, rotmat(1:, k:), nstate, &
2685 rmat_loc(1:, n2 + 1:), n1 + n2, 0.0_dp, gmat(:, :), nstate)
2686 ! no longer compiles in official sdgb:
2687 ! probably inefficient:
2688 CALL dgemm("N", "N", nstate, n1, n1, 1.0_dp, rotmat(1:, 1:), nstate, rmat_loc(n2 + 1:, n2 + 1:), &
2689 n1, 1.0_dp, gmat(:, :), nstate)
2690 END IF
2691
2692 CALL dcopy(nstate*n1, gmat(1, 1), 1, rotmat(1, 1), 1)
2693
2694 DO idim = 1, dim2
2695 DO i = 1, n1
2696 xyz_mix_ns(idim)%c_array(1:nstate, i) = z_zero
2697 END DO
2698
2699 DO istate = 1, n1
2700 DO jstate = 1, nstate
2701 DO i = 1, n2
2702 xyz_mix_ns(idim)%c_array(jstate, istate) = &
2703 xyz_mix_ns(idim)%c_array(jstate, istate) + &
2704 c_array_partner(jstate, i, idim)*rmat_loc(il2 + i - 1, il1 + istate - 1)
2705 END DO
2706 END DO
2707 END DO
2708 DO istate = 1, n1
2709 DO jstate = 1, nstate
2710 DO i = 1, n1
2711 xyz_mix_ns(idim)%c_array(jstate, istate) = xyz_mix_ns(idim)%c_array(jstate, istate) + &
2712 c_array_me(jstate, i, idim)*rmat_loc(il1 + i - 1, il1 + istate - 1)
2713 END DO
2714 END DO
2715 END DO
2716 END DO ! idim
2717
2718 DEALLOCATE (c_array_partner)
2719
2720 ELSE ! save my data
2721 DO idim = 1, dim2
2722 DO i = 1, ns_me
2723 xyz_mix_ns(idim)%c_array(1:nstate, i) = cz_ij_loc(idim)%c_array(1:nstate, i)
2724 END DO
2725 END DO
2726 END IF
2727
2728 DO idim = 1, dim2
2729 DO i = 1, ns_me
2730 cz_ij_loc(idim)%c_array(1:nstate, i) = z_zero
2731 END DO
2732 END DO
2733
2734 IF (ns_partner*ns_me /= 0) THEN
2735 ! transpose rotation matrix rmat_loc
2736 DO i = 1, ns_me + ns_partner
2737 DO j = i + 1, ns_me + ns_partner
2738 ri = rmat_loc(i, j)
2739 rmat_loc(i, j) = rmat_loc(j, i)
2740 rmat_loc(j, i) = ri
2741 END DO
2742 END DO
2743
2744 ! prepare for distribution
2745 DO i = 1, n1
2746 rmat_send(1:n1, i) = rmat_loc(il1:iu1, il1 + i - 1)
2747 END DO
2748 ik = nblock_max
2749 DO i = 1, n2
2750 rmat_send(ik + 1:ik + n1, i) = rmat_loc(il1:iu1, il2 + i - 1)
2751 END DO
2752 ELSE
2753 rmat_send = 0.0_dp
2754 END IF
2755
2756 ! collect data from all tasks (this takes some significant time)
2757 CALL para_env%allgather(rmat_send, rmat_recv_all)
2758
2759 ! update blocks everywhere
2760 DO ip = 0, para_env%num_pe - 1
2761
2762 ip_recv_from = mod(para_env%mepos - ip + para_env%num_pe, para_env%num_pe)
2763 rmat_recv(:, :) = rmat_recv_all(:, :, ip_recv_from)
2764
2765 ns_recv_from = ns_bound(ip_recv_from, 2) - ns_bound(ip_recv_from, 1) + 1
2766
2767 IF (ns_me /= 0) THEN
2768 IF (ns_recv_from /= 0) THEN
2769 !look for the partner of ip_recv_from
2770 ip_recv_partner = -1
2771 ns_recv_partner = 0
2772 DO ipair = 1, npair
2773 IF (list_pair(1, ipair) == ip_recv_from) THEN
2774 ip_recv_partner = list_pair(2, ipair)
2775 EXIT
2776 ELSE IF (list_pair(2, ipair) == ip_recv_from) THEN
2777 ip_recv_partner = list_pair(1, ipair)
2778 EXIT
2779 END IF
2780 END DO
2781
2782 IF (ip_recv_partner >= 0) THEN
2783 ns_recv_partner = ns_bound(ip_recv_partner, 2) - ns_bound(ip_recv_partner, 1) + 1
2784 END IF
2785 IF (ns_recv_partner > 0) THEN
2786 il1 = ns_bound(para_env%mepos, 1)
2787 il_recv = ns_bound(ip_recv_from, 1)
2788 il_recv_partner = ns_bound(ip_recv_partner, 1)
2789 ik = nblock_max
2790
2791 DO idim = 1, dim2
2792 DO i = 1, ns_recv_from
2793 ii = il_recv + i - 1
2794 DO j = 1, ns_me
2795 jj = j
2796 DO k = 1, ns_recv_from
2797 kk = il_recv + k - 1
2798 cz_ij_loc(idim)%c_array(ii, jj) = cz_ij_loc(idim)%c_array(ii, jj) + &
2799 rmat_recv(i, k)*xyz_mix_ns(idim)%c_array(kk, j)
2800 END DO
2801 END DO
2802 END DO
2803 DO i = 1, ns_recv_from
2804 ii = il_recv + i - 1
2805 DO j = 1, ns_me
2806 jj = j
2807 DO k = 1, ns_recv_partner
2808 kk = il_recv_partner + k - 1
2809 cz_ij_loc(idim)%c_array(ii, jj) = cz_ij_loc(idim)%c_array(ii, jj) + &
2810 rmat_recv(ik + i, k)*xyz_mix_ns(idim)%c_array(kk, j)
2811 END DO
2812 END DO
2813 END DO
2814 END DO ! idim
2815 ELSE
2816 il1 = ns_bound(para_env%mepos, 1)
2817 il_recv = ns_bound(ip_recv_from, 1)
2818 DO idim = 1, dim2
2819 DO j = 1, ns_me
2820 jj = j
2821 DO i = 1, ns_recv_from
2822 ii = il_recv + i - 1
2823 cz_ij_loc(idim)%c_array(ii, jj) = xyz_mix_ns(idim)%c_array(ii, j)
2824 END DO
2825 END DO
2826 END DO ! idim
2827 END IF
2828 END IF
2829 END IF ! ns_me
2830 END DO ! ip
2831
2832 IF (ns_partner*ns_me /= 0) THEN
2833 DEALLOCATE (rmat_loc)
2834 DO idim = 1, dim2
2835 DEALLOCATE (xyz_mix(idim)%c_array)
2836 END DO
2837 END IF
2838
2839 END DO ! iperm
2840
2841 ! calculate the max gradient
2842 DO idim = 1, dim2
2843 DO i = ns_bound(para_env%mepos, 1), ns_bound(para_env%mepos, 2)
2844 ii = i - ns_bound(para_env%mepos, 1) + 1
2845 zdiag_me(idim)%c_array(ii) = cz_ij_loc(idim)%c_array(i, ii)
2846 zdiag_me(idim)%c_array(ii) = cz_ij_loc(idim)%c_array(i, ii)
2847 END DO
2848 rcount(:) = SIZE(zdiag_me(idim)%c_array)
2849 rdispl(1) = 0
2850 DO ip = 2, para_env%num_pe
2851 rdispl(ip) = rdispl(ip - 1) + rcount(ip - 1)
2852 END DO
2853 ! collect all the diagonal elements in a replicated 1d array
2854 CALL para_env%allgatherv(zdiag_me(idim)%c_array, zdiag_all(idim)%c_array, rcount, rdispl)
2855 END DO
2856
2857 gmax = 0.0_dp
2858 DO j = ns_bound(para_env%mepos, 1), ns_bound(para_env%mepos, 2)
2859 k = j - ns_bound(para_env%mepos, 1) + 1
2860 DO i = 1, j - 1
2861 ! find the location of the diagonal element (i,i)
2862 DO ip = 0, para_env%num_pe - 1
2863 IF (i >= ns_bound(ip, 1) .AND. i <= ns_bound(ip, 2)) THEN
2864 ip_has_i = ip
2865 EXIT
2866 END IF
2867 END DO
2868 ii = nblock_max*ip_has_i + i - ns_bound(ip_has_i, 1) + 1
2869 ! mepos has the diagonal element (j,j), as well as the off diagonal (i,j)
2870 jj = nblock_max*para_env%mepos + j - ns_bound(para_env%mepos, 1) + 1
2871 grad = 0.0_dp
2872 DO idim = 1, dim2
2873 zi = zdiag_all(idim)%c_array(ii)
2874 zj = zdiag_all(idim)%c_array(jj)
2875 grad = grad + weights(idim)*real(4.0_dp*conjg(cz_ij_loc(idim)%c_array(i, k))*(zj - zi), dp)
2876 END DO
2877 gmax = max(gmax, abs(grad))
2878 END DO
2879 END DO
2880
2881 CALL para_env%max(gmax)
2882 tolerance = gmax
2883 IF (PRESENT(tol_out)) tol_out = tolerance
2884
2885 func = 0.0_dp
2886 DO i = ns_bound(para_env%mepos, 1), ns_bound(para_env%mepos, 2)
2887 k = i - ns_bound(para_env%mepos, 1) + 1
2888 DO idim = 1, dim2
2889 zr = real(cz_ij_loc(idim)%c_array(i, k), dp)
2890 zc = aimag(cz_ij_loc(idim)%c_array(i, k))
2891 func = func + weights(idim)*(1.0_dp - (zr*zr + zc*zc))/twopi/twopi
2892 END DO
2893 END DO
2894 CALL para_env%sum(func)
2895 t2 = m_walltime()
2896
2897 IF (output_unit > 0 .AND. modulo(sweeps, out_each) == 0) THEN
2898 WRITE (output_unit, '(T20,I12,T35,F20.10,T60,E12.4,F8.3)') sweeps, func, tolerance, t2 - t1
2899 CALL m_flush(output_unit)
2900 END IF
2901 IF (PRESENT(eps_localization)) THEN
2902 IF (tolerance < eps_localization) EXIT
2903 END IF
2904 IF (PRESENT(target_time) .AND. PRESENT(start_time)) THEN
2905 CALL external_control(should_stop, "LOC", target_time=target_time, start_time=start_time)
2906 IF (should_stop) EXIT
2907 END IF
2908
2909 END DO ! lsweep
2910
2911 ! buffer for message passing
2912 DEALLOCATE (rmat_recv_all)
2913
2914 DEALLOCATE (rmat_recv)
2915 DEALLOCATE (rmat_send)
2916 IF (ns_me > 0) THEN
2917 DEALLOCATE (c_array_me)
2918 END IF
2919 DO idim = 1, dim2
2920 DEALLOCATE (zdiag_me(idim)%c_array)
2921 DEALLOCATE (zdiag_all(idim)%c_array)
2922 END DO
2923 DEALLOCATE (zdiag_me)
2924 DEALLOCATE (zdiag_all)
2925 DEALLOCATE (xyz_mix)
2926 DO idim = 1, dim2
2927 IF (ns_me /= 0) THEN
2928 DEALLOCATE (xyz_mix_ns(idim)%c_array)
2929 END IF
2930 END DO
2931 DEALLOCATE (xyz_mix_ns)
2932 IF (ns_me /= 0) THEN
2933 DEALLOCATE (gmat)
2934 END IF
2935 DEALLOCATE (mii)
2936 DEALLOCATE (mij)
2937 DEALLOCATE (mjj)
2938 DEALLOCATE (list_pair)
2939
2940 END SUBROUTINE jacobi_rot_para_core
2941
2942! **************************************************************************************************
2943!> \brief ...
2944!> \param iperm ...
2945!> \param para_env ...
2946!> \param list_pair ...
2947! **************************************************************************************************
2948 SUBROUTINE eberlein(iperm, para_env, list_pair)
2949 INTEGER, INTENT(IN) :: iperm
2950 TYPE(mp_para_env_type), POINTER :: para_env
2951 INTEGER, DIMENSION(:, :) :: list_pair
2952
2953 INTEGER :: i, ii, jj, npair
2954
2955 npair = (para_env%num_pe + 1)/2
2956 IF (iperm == 1) THEN
2957!..set up initial ordering
2958 DO i = 0, para_env%num_pe - 1
2959 ii = ((i + 1) + 1)/2
2960 jj = mod((i + 1) + 1, 2) + 1
2961 list_pair(jj, ii) = i
2962 END DO
2963 IF (mod(para_env%num_pe, 2) == 1) list_pair(2, npair) = -1
2964 ELSE IF (mod(iperm, 2) == 0) THEN
2965!..a type shift
2966 jj = list_pair(1, npair)
2967 DO i = npair, 3, -1
2968 list_pair(1, i) = list_pair(1, i - 1)
2969 END DO
2970 list_pair(1, 2) = list_pair(2, 1)
2971 list_pair(2, 1) = jj
2972 ELSE
2973!..b type shift
2974 jj = list_pair(2, 1)
2975 DO i = 1, npair - 1
2976 list_pair(2, i) = list_pair(2, i + 1)
2977 END DO
2978 list_pair(2, npair) = jj
2979 END IF
2980
2981 END SUBROUTINE eberlein
2982
2983! **************************************************************************************************
2984!> \brief ...
2985!> \param vectors ...
2986!> \param op_sm_set ...
2987!> \param zij_fm_set ...
2988! **************************************************************************************************
2989 SUBROUTINE zij_matrix(vectors, op_sm_set, zij_fm_set)
2990
2991 TYPE(cp_fm_type), INTENT(IN) :: vectors
2992 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: op_sm_set
2993 TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: zij_fm_set
2994
2995 CHARACTER(len=*), PARAMETER :: routinen = 'zij_matrix'
2996
2997 INTEGER :: handle, i, j, nao, nmoloc
2998 TYPE(cp_fm_type) :: opvec
2999
3000 CALL timeset(routinen, handle)
3001
3002 ! get rows and cols of the input
3003 CALL cp_fm_get_info(vectors, nrow_global=nao, ncol_global=nmoloc)
3004 ! replicate the input kind of matrix
3005 CALL cp_fm_create(opvec, vectors%matrix_struct)
3006
3007 ! Compute zij here
3008 DO i = 1, SIZE(zij_fm_set, 2)
3009 DO j = 1, SIZE(zij_fm_set, 1)
3010 CALL cp_fm_set_all(zij_fm_set(j, i), 0.0_dp)
3011 CALL cp_dbcsr_sm_fm_multiply(op_sm_set(j, i)%matrix, vectors, opvec, ncol=nmoloc)
3012 CALL parallel_gemm("T", "N", nmoloc, nmoloc, nao, 1.0_dp, vectors, opvec, 0.0_dp, &
3013 zij_fm_set(j, i))
3014 END DO
3015 END DO
3016
3017 CALL cp_fm_release(opvec)
3018 CALL timestop(handle)
3019
3020 END SUBROUTINE zij_matrix
3021
3022! **************************************************************************************************
3023!> \brief ...
3024!> \param vectors ...
3025! **************************************************************************************************
3026 SUBROUTINE scdm_qrfact(vectors)
3027
3028 TYPE(cp_fm_type), INTENT(IN) :: vectors
3029
3030 CHARACTER(len=*), PARAMETER :: routinen = 'scdm_qrfact'
3031
3032 INTEGER :: handle, ncolt, nrowt
3033 REAL(kind=dp), DIMENSION(:), POINTER :: tau
3034 TYPE(cp_fm_struct_type), POINTER :: cstruct
3035 TYPE(cp_fm_type) :: ctp, qf, tmp
3036
3037 CALL timeset(routinen, handle)
3038
3039 ! Create Transpose of Coefficient Matrix vectors
3040 nrowt = vectors%matrix_struct%ncol_global
3041 ncolt = vectors%matrix_struct%nrow_global
3042
3043 CALL cp_fm_struct_create(cstruct, template_fmstruct=vectors%matrix_struct, &
3044 nrow_global=nrowt, ncol_global=ncolt)
3045 CALL cp_fm_create(ctp, cstruct)
3046 CALL cp_fm_struct_release(cstruct)
3047
3048 ALLOCATE (tau(nrowt))
3049
3050 CALL cp_fm_transpose(vectors, ctp)
3051
3052 ! Get QR decomposition of CTs
3053 CALL cp_fm_pdgeqpf(ctp, tau, nrowt, ncolt, 1, 1)
3054
3055 ! Construction of Q from the scalapack output
3056 CALL cp_fm_struct_create(cstruct, para_env=ctp%matrix_struct%para_env, &
3057 context=ctp%matrix_struct%context, nrow_global=ctp%matrix_struct%nrow_global, &
3058 ncol_global=ctp%matrix_struct%nrow_global)
3059 CALL cp_fm_create(qf, cstruct)
3060 CALL cp_fm_struct_release(cstruct)
3061 CALL cp_fm_to_fm_submat(ctp, qf, nrowt, nrowt, 1, 1, 1, 1)
3062
3063 ! Get Q
3064 CALL cp_fm_pdorgqr(qf, tau, nrowt, 1, 1)
3065
3066 ! Transform original coefficient matrix vectors
3067 CALL cp_fm_create(tmp, vectors%matrix_struct)
3068 CALL cp_fm_set_all(tmp, 0.0_dp, 1.0_dp)
3069 CALL cp_fm_to_fm(vectors, tmp)
3070 CALL parallel_gemm('N', 'N', ncolt, nrowt, nrowt, 1.0_dp, tmp, qf, 0.0_dp, vectors)
3071
3072 ! Cleanup
3073 CALL cp_fm_release(ctp)
3074 CALL cp_fm_release(tmp)
3075 CALL cp_fm_release(qf)
3076 DEALLOCATE (tau)
3077
3078 CALL timestop(handle)
3079
3080 END SUBROUTINE scdm_qrfact
3081
3082! **************************************************************************************************
3083!> \brief Achieves minimisation of the spread functional by simultaneous diagonalisation with Jacobi
3084!> rotations as presented in Cardoso & Souloumiac, SIAM J. Matrix Anal. Appl., 17(1), 161.
3085!> Generalizes the Jacobi algorithm to complex matrices.
3086!> \param weights array of weights for calculating the total spread
3087!> \param zij spread operator matrices
3088!> \param max_iter maximum number iterations
3089!> \param eps_localization numerical tolerance
3090!> \param sweeps counts number of sweeps required
3091!> \param out_each how often to print info
3092!> \param vectors complex vectors to be localized
3093!> \par History
3094!> 2020-04 created [LS]
3095!> \author Lukas Schreder
3096! **************************************************************************************************
3097 SUBROUTINE cardoso_souloumiac(weights, zij, max_iter, eps_localization, sweeps, &
3098 out_each, vectors)
3099
3100 REAL(kind=dp), INTENT(IN) :: weights(:)
3101 TYPE(cp_cfm_type), INTENT(INOUT) :: zij(:, :)
3102 INTEGER, INTENT(IN) :: max_iter
3103 REAL(kind=dp), INTENT(IN) :: eps_localization
3104 INTEGER :: sweeps
3105 INTEGER, INTENT(IN) :: out_each
3106 TYPE(cp_cfm_type), POINTER :: vectors
3107
3108 CHARACTER(len=*), PARAMETER :: routinen = 'cardoso_souloumiac'
3109
3110 COMPLEX(KIND=dp) :: s
3111 COMPLEX(KIND=dp), ALLOCATABLE :: mii(:), mij(:), mji(:), mjj(:)
3112 INTEGER :: dim1, dim2, handle, idim, istate, jdim, &
3113 jstate, nstate, unit_nr
3114 REAL(kind=dp) :: c, old_spread, spread, t1, t2, tolerance
3115 TYPE(cp_cfm_type), POINTER :: c_rmat, c_zij(:)
3116
3117 CALL timeset(routinen, handle)
3118
3119 dim1 = SIZE(zij, 1)
3120 dim2 = SIZE(zij, 2)
3121
3122 NULLIFY (c_rmat, c_zij)
3123 ALLOCATE (c_rmat, c_zij(dim1*dim2), mii(dim1*dim2), mij(dim1*dim2), mji(dim1*dim2), mjj(dim1*dim2))
3124 CALL cp_cfm_create(c_rmat, zij(1, 1)%matrix_struct)
3125 CALL cp_cfm_set_all(c_rmat, (0.0_dp, 0.0_dp), (1.0_dp, 0.0_dp)) ! start with the identity transformation
3126 DO idim = 1, dim2
3127 DO jdim = 1, dim1
3128 CALL cp_cfm_create(c_zij((idim - 1)*dim1 + jdim), zij(jdim, idim)%matrix_struct)
3129 CALL cp_cfm_to_cfm(zij(jdim, idim), c_zij((idim - 1)*dim1 + jdim))
3130 END DO
3131 END DO
3132
3133 CALL cp_cfm_get_info(c_rmat, nrow_global=nstate)
3134 tolerance = 1.0e10_dp
3135 old_spread = 1.0e10_dp
3136 sweeps = 0
3137 unit_nr = -1
3138 IF (c_rmat%matrix_struct%para_env%is_source()) THEN
3139 unit_nr = cp_logger_get_default_unit_nr()
3140 WRITE (unit_nr, "(T4,A )") " Localization by iterative Jacobi rotation using "// &
3141 "Cardoso-Souloumiac angles"
3142
3143 END IF
3144
3145 ! Jacobi sweeps until converged
3146 DO WHILE (tolerance >= eps_localization .AND. sweeps < max_iter)
3147 sweeps = sweeps + 1
3148 t1 = m_walltime()
3149 DO jstate = 1, nstate
3150 DO istate = jstate + 1, nstate
3151 DO idim = 1, dim1*dim2
3152 CALL cp_cfm_get_element(c_zij(idim), istate, istate, mii(idim))
3153 CALL cp_cfm_get_element(c_zij(idim), istate, jstate, mij(idim))
3154 CALL cp_cfm_get_element(c_zij(idim), jstate, istate, mji(idim))
3155 CALL cp_cfm_get_element(c_zij(idim), jstate, jstate, mjj(idim))
3156 END DO
3157 CALL get_cardoso_angles(mii, mij, mji, mjj, c, s, vectors%matrix_struct)
3158 DO idim = 1, dim1*dim2
3159 CALL cp_cfm_rot_cols(c_zij(idim), istate, jstate, c, real(s))
3160 CALL cp_cfm_rot_rows(c_zij(idim), istate, jstate, c, real(s))
3161 END DO
3162 CALL cp_cfm_rot_cols(c_rmat, istate, jstate, c, real(s))
3163 END DO
3164 END DO
3165 CALL check_tolerance(c_zij, weights, spread)
3166 CALL check_tolerance(c_zij, weights, tolerance)
3167 tolerance = abs(spread - old_spread)
3168 old_spread = spread
3169
3170 t2 = m_walltime()
3171 IF (unit_nr > 0 .AND. modulo(sweeps, out_each) == 0) THEN
3172 WRITE (unit_nr, "(T4,A,I7,A,E12.4,A,E12.4,A,F8.3)") &
3173 "Iteration:", sweeps, "Functional", spread, "Tolerance:", tolerance, "Time:", t2 - t1
3174 CALL m_flush(unit_nr)
3175 END IF
3176 END DO
3177
3178 DO idim = 1, dim2
3179 ! back to an interlaced matrix
3180 DO jdim = 1, dim1
3181 CALL cp_cfm_to_cfm(c_zij((idim - 1)*dim1 + jdim), zij(jdim, idim))
3182 CALL cp_cfm_release(c_zij((idim - 1)*dim1 + jdim))
3183 END DO
3184 END DO
3185
3186 CALL rotate_orbitals_cfm(c_rmat, vectors)
3187
3188 DEALLOCATE (c_zij, mii, mij, mji, mjj)
3189 CALL cp_cfm_release(c_rmat)
3190 DEALLOCATE (c_rmat)
3191
3192 CALL timestop(handle)
3193
3194 END SUBROUTINE cardoso_souloumiac
3195
3196! **************************************************************************************************
3197!> \brief Pipek-Mezey version of the Cardoso-Souloumiac PADE algorithm for complex-valued matrices.
3198!> \param zij spread operator matrices
3199!> \param vec complex vectors to be localised
3200!> \param sweeps counts number of sweeps required
3201!> \param max_iter maximum number iterations
3202!> \param eps numerical tolerance
3203!> \param out_each how often to print info
3204!> \par History
3205!> 2020-04 created [LS]
3206!> \author Lukas Schreder
3207! **************************************************************************************************
3208 SUBROUTINE cardoso_souloumiac_pipek(zij, vec, sweeps, max_iter, eps, out_each)
3209
3210 TYPE(cp_cfm_type), POINTER :: zij(:, :), vec
3211 INTEGER :: sweeps, max_iter
3212 REAL(dp) :: eps
3213 INTEGER :: out_each
3214
3215 CHARACTER(*), PARAMETER :: routinen = 'cardoso_souloumiac_pipek'
3216
3217 COMPLEX(dp) :: c_spread, s
3218 COMPLEX(dp), POINTER :: qii(:), qij(:), qji(:), qjj(:)
3219 INTEGER :: handle, i, j, k, n_dim, n_states, &
3220 output_unit
3221 REAL(dp) :: c, old_spread, spread, t1, t2, tol
3222 TYPE(cp_cfm_type), POINTER :: c_zij(:), rmat
3223
3224 CALL timeset(routinen, handle)
3225
3226 CALL cite_reference(schreder2024_2)
3227
3228 n_dim = SIZE(zij)
3229
3230 ALLOCATE (rmat)
3231 CALL cp_cfm_create(rmat, zij(1, 1)%matrix_struct)
3232
3233 CALL cp_cfm_set_all(rmat, (0.0_dp, 0.0_dp), (1.0_dp, 0.0_dp))
3234 ALLOCATE (c_zij(n_dim), qii(n_dim), qij(n_dim), qji(n_dim), qjj(n_dim))
3235
3236 ! build Qij matrix
3237 DO k = 1, n_dim
3238 c_zij(k) = zij(k, 1)
3239 END DO
3240
3241 CALL cp_cfm_get_info(c_zij(1), ncol_global=n_states)
3242
3243 tol = 1.0e10_dp
3244 c_spread = (0.0_dp, 0.0_dp)
3245 DO k = 1, n_dim
3246 DO i = 1, n_states
3247 CALL cp_cfm_get_element(c_zij(k), i, i, s)
3248 c_spread = c_spread + s*s
3249 END DO
3250 END DO
3251 spread = real(c_spread)
3252 old_spread = spread
3253
3254 sweeps = 0
3255 output_unit = cp_logger_get_default_unit_nr()
3256 WRITE (output_unit, "(T4,A )") " Localization by iterative Jacobi rotation using "// &
3257 "Cardoso-Souloumiac angles"
3258 WRITE (output_unit, "(T4,A )") " and Pipek-Mezey spread functional"
3259
3260 IF (output_unit > 0 .AND. modulo(sweeps, out_each) == 0) THEN
3261 WRITE (output_unit, "(T4,A,I7,A,E12.4,A,E12.4,A,F8.3)") &
3262 "Iteration:", sweeps, " Functional", spread, " Tolerance:", tol, " Time:", 0.0_dp
3263 END IF
3264
3265 DO WHILE (tol >= eps .AND. sweeps < max_iter)
3266 t1 = m_walltime()
3267 sweeps = sweeps + 1
3268
3269 DO i = 1, n_states
3270 DO j = i + 1, n_states
3271 DO k = 1, n_dim
3272 CALL cp_cfm_get_element(c_zij(k), i, i, qii(k))
3273 CALL cp_cfm_get_element(c_zij(k), i, j, qij(k))
3274 CALL cp_cfm_get_element(c_zij(k), j, i, qji(k))
3275 CALL cp_cfm_get_element(c_zij(k), j, j, qjj(k))
3276 END DO
3277 CALL get_cardoso_angles(qii, qij, qji, qjj, c, s, vec%matrix_struct)
3278 DO k = 1, n_dim
3279 CALL cp_cfm_rot_cols(c_zij(k), i, j, c, real(s))
3280 CALL cp_cfm_rot_rows(c_zij(k), i, j, c, real(s))
3281 END DO
3282 CALL cp_cfm_rot_cols(rmat, i, j, c, real(s))
3283 END DO
3284 END DO
3285
3286 c_spread = (0.0_dp, 0.0_dp)
3287 DO i = 1, n_states
3288 DO k = 1, n_dim
3289 CALL cp_cfm_get_element(c_zij(k), i, i, s)
3290 c_spread = c_spread + s*s
3291 END DO
3292 END DO
3293 spread = real(c_spread)
3294
3295 tol = abs(spread - old_spread)
3296 old_spread = spread
3297 t2 = m_walltime()
3298 IF (output_unit > 0 .AND. modulo(sweeps, out_each) == 0) THEN
3299 WRITE (output_unit, "(T4,A,I7,A,E12.4,A,E12.4,A,F8.3)") &
3300 "Iteration:", sweeps, " Functional", spread, " Tolerance:", tol, " Time:", t2 - t1
3301 END IF
3302 END DO
3303
3304 CALL rotate_orbitals_cfm(rmat, vec)
3305
3306 ! c_zij(k) was assigned (shallow copy of zij(k,1)); do not release contents
3307 DEALLOCATE (c_zij, qii, qij, qji, qjj)
3308 CALL cp_cfm_release(rmat)
3309 DEALLOCATE (rmat)
3310
3311 CALL timestop(handle)
3312
3313 END SUBROUTINE cardoso_souloumiac_pipek
3314
3315! **************************************************************************************************
3316!> \brief calculates the Jacobi angles needed in serial Cardoso-Souloumiac diagonalisation
3317!> \param mii ...
3318!> \param mij ...
3319!> \param mji ...
3320!> \param mjj ...
3321!> \param c ...
3322!> \param s ...
3323!> \param tmp_fm_struct ...
3324!> \par History
3325!> 2020-04 created [LS]
3326!> \author Lukas Schreder
3327! **************************************************************************************************
3328 SUBROUTINE get_cardoso_angles(mii, mij, mji, mjj, c, s, tmp_fm_struct)
3329
3330 COMPLEX(KIND=dp), DIMENSION(:) :: mii, mij, mji, mjj
3331 REAL(kind=dp), INTENT(out) :: c
3332 COMPLEX(KIND=dp), INTENT(out) :: s
3333 TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct
3334
3335 INTEGER :: dim_m, i, i_max
3336 REAL(kind=dp) :: r, x, y, z
3337 REAL(kind=dp), DIMENSION(3) :: evals
3338 TYPE(cp_cfm_type), POINTER :: c_gmat, hmat
3339 TYPE(cp_fm_struct_type), POINTER :: g_fm_struct, h_fm_struct
3340 TYPE(cp_fm_type), POINTER :: evects, gmat
3341
3342 dim_m = SIZE(mii)
3343
3344 CALL cp_fm_struct_create(h_fm_struct, nrow_global=1, ncol_global=3, &
3345 template_fmstruct=tmp_fm_struct)
3346 CALL cp_fm_struct_create(g_fm_struct, nrow_global=3, ncol_global=3, &
3347 template_fmstruct=tmp_fm_struct)
3348 ALLOCATE (hmat, c_gmat, gmat, evects)
3349 CALL cp_cfm_create(hmat, h_fm_struct)
3350 CALL cp_cfm_create(c_gmat, g_fm_struct)
3351 CALL cp_cfm_set_all(c_gmat, (0.0_dp, 0.0_dp), (0.0_dp, 0.0_dp))
3352 CALL cp_fm_create(gmat, g_fm_struct)
3353 CALL cp_fm_create(evects, g_fm_struct)
3354
3355 DO i = 1, dim_m
3356 CALL cp_cfm_set_element(hmat, 1, 1, (mii(i) - mjj(i)))
3357 CALL cp_cfm_set_element(hmat, 1, 2, (mij(i) + mji(i)))
3358 CALL cp_cfm_set_element(hmat, 1, 3, (0.0_dp, 1.0_dp)*(mji(i) - mij(i)))
3359 CALL cp_cfm_gemm("C", "N", 3, 3, 1, (1.0_dp, 0.0_dp), hmat, hmat, (1.0_dp, 0.0_dp), c_gmat)
3360 END DO
3361 CALL cp_cfm_to_fm(c_gmat, gmat)
3362
3363 ! find eigenvector with highest eigenvalue
3364 CALL choose_eigv_solver(gmat, evects, evals)
3365 i_max = maxloc(evals, 1)
3366 CALL cp_fm_get_element(evects, 1, i_max, x)
3367 CALL cp_fm_get_element(evects, 2, i_max, y)
3368 CALL cp_fm_get_element(evects, 3, i_max, z)
3369 IF (x < 0) THEN
3370 x = -x
3371 y = -y
3372 z = -z
3373 END IF
3374
3375 ! calculate the angles
3376 r = sqrt(x**2 + y**2 + z**2) ! always 1
3377 c = sqrt((x + r)/(2*r)) ! always real
3378 s = (y - gaussi*z)/sqrt(2*r*(x + r)) ! always complex
3379
3380 CALL cp_fm_struct_release(h_fm_struct)
3381 CALL cp_fm_struct_release(g_fm_struct)
3382 CALL cp_cfm_release(hmat)
3383 CALL cp_cfm_release(c_gmat)
3384 CALL cp_fm_release(gmat)
3385 CALL cp_fm_release(evects)
3386 DEALLOCATE (hmat, c_gmat, gmat, evects)
3387
3388 END SUBROUTINE get_cardoso_angles
3389
3390END MODULE qs_localization_methods
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition grid_common.h:125
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:225
cp_cfm_basic_linalg::cp_cfm_scale
Definition cp_cfm_basic_linalg.F:60
cp_cfm_types::cp_cfm_to_cfm
Definition cp_cfm_types.F:56
cp_fm_basic_linalg::cp_fm_trace
Definition cp_fm_basic_linalg.F:78
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:88
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:83
kahan_sum::accurate_sum
Definition kahan_sum.F:41
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:39
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition bibliography.F:21
bibliography::schreder2024_2
integer, save, public schreder2024_2
Definition bibliography.F:36
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cp_blacs_env
methods related to the blacs parallel environment
Definition cp_blacs_env.F:15
cp_cfm_basic_linalg
Basic linear algebra operations for complex full matrices.
Definition cp_cfm_basic_linalg.F:16
cp_cfm_basic_linalg::cp_cfm_scale_and_add
subroutine, public cp_cfm_scale_and_add(alpha, matrix_a, beta, matrix_b)
Scale and add two BLACS matrices (a = alpha*a + beta*b).
Definition cp_cfm_basic_linalg.F:240
cp_cfm_basic_linalg::cp_cfm_gemm
subroutine, public cp_cfm_gemm(transa, transb, m, n, k, alpha, matrix_a, matrix_b, beta, matrix_c, a_first_col, a_first_row, b_first_col, b_first_row, c_first_col, c_first_row)
Performs one of the matrix-matrix operations: matrix_c = alpha * op1( matrix_a ) * op2( matrix_b ) + ...
Definition cp_cfm_basic_linalg.F:559
cp_cfm_basic_linalg::cp_cfm_rot_rows
subroutine, public cp_cfm_rot_rows(matrix, irow, jrow, cs, sn)
Applies a planar rotation defined by cs and sn to the i'th and j'th rows.
Definition cp_cfm_basic_linalg.F:1254
cp_cfm_basic_linalg::cp_cfm_schur_product
subroutine, public cp_cfm_schur_product(matrix_a, matrix_b, matrix_c)
Computes the element-wise (Schur) product of two matrices: C = A \circ B .
Definition cp_cfm_basic_linalg.F:153
cp_cfm_basic_linalg::cp_cfm_column_scale
subroutine, public cp_cfm_column_scale(matrix_a, scaling)
Scales columns of the full matrix by corresponding factors.
Definition cp_cfm_basic_linalg.F:629
cp_cfm_basic_linalg::cp_cfm_rot_cols
subroutine, public cp_cfm_rot_cols(matrix, icol, jcol, cs, sn)
Applies a planar rotation defined by cs and sn to the i'th and j'th columnns.
Definition cp_cfm_basic_linalg.F:1301
cp_cfm_basic_linalg::cp_cfm_trace
subroutine, public cp_cfm_trace(matrix_a, matrix_b, trace)
Returns the trace of matrix_a^T matrix_b, i.e sum_{i,j}(matrix_a(i,j)*matrix_b(i,j)) .
Definition cp_cfm_basic_linalg.F:906
cp_cfm_diag
used for collecting diagonalization schemes available for cp_cfm_type
Definition cp_cfm_diag.F:14
cp_cfm_diag::cp_cfm_heevd
subroutine, public cp_cfm_heevd(matrix, eigenvectors, eigenvalues)
Perform a diagonalisation of a complex matrix.
Definition cp_cfm_diag.F:74
cp_cfm_types
Represents a complex full matrix distributed on many processors.
Definition cp_cfm_types.F:12
cp_cfm_types::cp_cfm_get_submatrix
subroutine, public cp_cfm_get_submatrix(fm, target_m, start_row, start_col, n_rows, n_cols, transpose)
Extract a sub-matrix from the full matrix: op(target_m)(1:n_rows,1:n_cols) = fm(start_row:start_row+n...
Definition cp_cfm_types.F:356
cp_cfm_types::cp_cfm_get_element
subroutine, public cp_cfm_get_element(matrix, irow_global, icol_global, alpha)
Get the matrix element by its global index.
Definition cp_cfm_types.F:255
cp_cfm_types::cp_cfm_release
subroutine, public cp_cfm_release(matrix)
Releases a full matrix.
Definition cp_cfm_types.F:182
cp_cfm_types::cp_cfm_set_element
subroutine, public cp_cfm_set_element(matrix, irow_global, icol_global, alpha)
Set the matrix element (irow_global,icol_global) of the full matrix to alpha.
Definition cp_cfm_types.F:302
cp_cfm_types::cp_fm_to_cfm
subroutine, public cp_fm_to_cfm(msourcer, msourcei, mtarget)
Construct a complex full matrix by taking its real and imaginary parts from two separate real-value f...
Definition cp_cfm_types.F:844
cp_cfm_types::cp_cfm_create
subroutine, public cp_cfm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
Creates a new full matrix with the given structure.
Definition cp_cfm_types.F:125
cp_cfm_types::cp_cfm_get_info
subroutine, public cp_cfm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, matrix_struct, para_env)
Returns information about a full matrix.
Definition cp_cfm_types.F:631
cp_cfm_types::cp_cfm_set_submatrix
subroutine, public cp_cfm_set_submatrix(matrix, new_values, start_row, start_col, n_rows, n_cols, alpha, beta, transpose)
Set a sub-matrix of the full matrix: matrix(start_row:start_row+n_rows,start_col:start_col+n_cols) = ...
Definition cp_cfm_types.F:494
cp_cfm_types::cp_cfm_set_all
subroutine, public cp_cfm_set_all(matrix, alpha, beta)
Set all elements of the full matrix to alpha. Besides, set all diagonal matrix elements to beta (if g...
Definition cp_cfm_types.F:202
cp_cfm_types::cp_cfm_to_fm
subroutine, public cp_cfm_to_fm(msource, mtargetr, mtargeti)
Copy real and imaginary parts of a complex full matrix into separate real-value full matrices.
Definition cp_cfm_types.F:792
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::cp_dbcsr_sm_fm_multiply
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
Definition cp_dbcsr_operations.F:548
cp_external_control
Routines to handle the external control of CP2K.
Definition cp_external_control.F:15
cp_external_control::external_control
subroutine, public external_control(should_stop, flag, globenv, target_time, start_time, force_check)
External manipulations during a run : when the <PROJECT_NAME>.EXIT_$runtype command is sent the progr...
Definition cp_external_control.F:91
cp_fm_basic_linalg
Basic linear algebra operations for full matrices.
Definition cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_pdgeqpf
subroutine, public cp_fm_pdgeqpf(matrix, tau, nrow, ncol, first_row, first_col)
compute a QR factorization with column pivoting of a M-by-N distributed matrix sub( A ) = A(IA:IA+M-1...
Definition cp_fm_basic_linalg.F:2438
cp_fm_basic_linalg::cp_fm_frobenius_norm
real(kind=dp) function, public cp_fm_frobenius_norm(matrix_a)
computes the Frobenius norm of matrix_a
Definition cp_fm_basic_linalg.F:621
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:1532
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:170
cp_fm_basic_linalg::cp_fm_scale
subroutine, public cp_fm_scale(alpha, matrix_a)
scales a matrix matrix_a = alpha * matrix_b
Definition cp_fm_basic_linalg.F:1502
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:1413
cp_fm_basic_linalg::cp_fm_pdorgqr
subroutine, public cp_fm_pdorgqr(matrix, tau, nrow, first_row, first_col)
generates an M-by-N real distributed matrix Q denoting A(IA:IA+M-1,JA:JA+N-1) with orthonormal column...
Definition cp_fm_basic_linalg.F:2510
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:61
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:573
cp_fm_diag::choose_eigv_solver
subroutine, public choose_eigv_solver(matrix, eigenvectors, eigenvalues, info)
Choose the Eigensolver depending on which library is available ELPA seems to be unstable for small sy...
Definition cp_fm_diag.F:245
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:132
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:444
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:355
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:1038
cp_fm_types::cp_fm_get_element
subroutine, public cp_fm_get_element(matrix, irow_global, icol_global, alpha, local)
returns an element of a fm this value is valid on every cpu using this call is expensive
Definition cp_fm_types.F:627
cp_fm_types::cp_fm_set_submatrix
subroutine, public cp_fm_set_submatrix(fm, new_values, start_row, start_col, n_rows, n_cols, alpha, beta, transpose)
sets a submatrix of a full matrix fm(start_row:start_row+n_rows,start_col:start_col+n_cols) = alpha*o...
Definition cp_fm_types.F:750
cp_fm_types::cp_fm_maxabsrownorm
subroutine, public cp_fm_maxabsrownorm(matrix, a_max)
find the maximum over the rows of the sum of the absolute values of the elements of a given row = || ...
Definition cp_fm_types.F:1174
cp_fm_types::cp_fm_to_fm_submat
subroutine, public cp_fm_to_fm_submat(msource, mtarget, nrow, ncol, s_firstrow, s_firstcol, t_firstrow, t_firstcol)
copy just a part ot the matrix
Definition cp_fm_types.F:1455
cp_fm_types::cp_fm_maxabsval
subroutine, public cp_fm_maxabsval(matrix, a_max, ir_max, ic_max)
find the maximum absolute value of the matrix element maxval(abs(matrix))
Definition cp_fm_types.F:1086
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition cp_fm_types.F:534
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:165
cp_fm_types::cp_fm_get_submatrix
subroutine, public cp_fm_get_submatrix(fm, target_m, start_row, start_col, n_rows, n_cols, transpose)
gets a submatrix of a full matrix op(target_m)(1:n_rows,1:n_cols) =fm(start_row:start_row+n_rows,...
Definition cp_fm_types.F:925
cp_fm_types::cp_fm_init_random
subroutine, public cp_fm_init_random(matrix, ncol, start_col)
fills a matrix with random numbers
Definition cp_fm_types.F:452
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_unit_nr
recursive integer function, public cp_logger_get_default_unit_nr(logger, local, skip_not_ionode)
asks the default unit number of the given logger. try to use cp_logger_get_unit_nr
Definition cp_log_handling.F:576
cp_log_handling::cp_logger_get_default_io_unit
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Definition cp_log_handling.F:523
kahan_sum
sums arrays of real/complex numbers with much reduced round-off as compared to a naive implementation...
Definition kahan_sum.F:29
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
machine
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
machine::m_flush
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition machine.F:124
machine::m_walltime
real(kind=dp) function, public m_walltime()
returns time from a real-time clock, protected against rolling early/easily
Definition machine.F:141
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
mathconstants::z_one
complex(kind=dp), parameter, public z_one
Definition mathconstants.F:143
mathconstants::gaussi
complex(kind=dp), parameter, public gaussi
Definition mathconstants.F:142
mathconstants::twopi
real(kind=dp), parameter, public twopi
Definition mathconstants.F:112
mathconstants::z_zero
complex(kind=dp), parameter, public z_zero
Definition mathconstants.F:143
matrix_exp
Routines for calculating a complex matrix exponential.
Definition matrix_exp.F:13
matrix_exp::exp_pade_real
subroutine, public exp_pade_real(exp_h, matrix, nsquare, npade)
exponential of a real matrix, calculated using pade approximation together with scaling and squaring
Definition matrix_exp.F:570
matrix_exp::get_nsquare_norder
subroutine, public get_nsquare_norder(norm, nsquare, norder, eps_exp, method, do_emd)
optimization function for pade/taylor order and number of squaring steps
Definition matrix_exp.F:244
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition parallel_gemm_api.F:14
qs_localization_methods
Localization methods such as 2x2 Jacobi rotations Steepest Decents Conjugate Gradient.
Definition qs_localization_methods.F:18
qs_localization_methods::approx_l1_norm_sd
subroutine, public approx_l1_norm_sd(c, iterations, eps, converged, sweeps)
...
Definition qs_localization_methods.F:102
qs_localization_methods::rotate_orbitals
subroutine, public rotate_orbitals(rmat, vectors)
...
Definition qs_localization_methods.F:1386
qs_localization_methods::cardoso_souloumiac
subroutine, public cardoso_souloumiac(weights, zij, max_iter, eps_localization, sweeps, out_each, vectors)
Achieves minimisation of the spread functional by simultaneous diagonalisation with Jacobi rotations ...
Definition qs_localization_methods.F:3099
qs_localization_methods::jacobi_rotations
subroutine, public jacobi_rotations(weights, zij, vectors, para_env, max_iter, eps_localization, sweeps, out_each, target_time, start_time, restricted)
wrapper for the jacobi routines, should be removed if jacobi_rot_para can deal with serial para_envs.
Definition qs_localization_methods.F:285
qs_localization_methods::zij_matrix
subroutine, public zij_matrix(vectors, op_sm_set, zij_fm_set)
...
Definition qs_localization_methods.F:2990
qs_localization_methods::direct_mini
subroutine, public direct_mini(weights, zij, vectors, max_iter, eps_localization, iterations)
use the exponential parametrization as described in to perform a direct mini Gerd Berghold et al....
Definition qs_localization_methods.F:1791
qs_localization_methods::initialize_weights
subroutine, public initialize_weights(cell, weights)
...
Definition qs_localization_methods.F:246
qs_localization_methods::crazy_rotations
subroutine, public crazy_rotations(weights, zij, vectors, max_iter, max_crazy_angle, crazy_scale, crazy_use_diag, eps_localization, iterations, converged)
yet another crazy try, computes the angles needed to rotate the orbitals first and rotates them all a...
Definition qs_localization_methods.F:1593
qs_localization_methods::scdm_qrfact
subroutine, public scdm_qrfact(vectors)
...
Definition qs_localization_methods.F:3027
qs_localization_methods::cardoso_souloumiac_pipek
subroutine, public cardoso_souloumiac_pipek(zij, vec, sweeps, max_iter, eps, out_each)
Pipek-Mezey version of the Cardoso-Souloumiac PADE algorithm for complex-valued matrices.
Definition qs_localization_methods.F:3209
qs_localization_methods::jacobi_cg_edf_ls
subroutine, public jacobi_cg_edf_ls(para_env, weights, zij, vectors, max_iter, eps_localization, iter, out_each, nextra, do_cg, nmo, vectors_2, mos_guess)
combine jacobi rotations (serial) and conjugate gradient with golden section line search for partiall...
Definition qs_localization_methods.F:582
xc
Exchange and Correlation functional calculations.
Definition xc.F:17
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
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_cfm_types::cp_cfm_type
Represent a complex full matrix.
Definition cp_cfm_types.F:69
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
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:114
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:746