Line data Source code
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 : ! **************************************************************************************************
18 : MODULE 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 :
92 : CONTAINS
93 : ! **************************************************************************************************
94 : !> \brief ...
95 : !> \param C ...
96 : !> \param iterations ...
97 : !> \param eps ...
98 : !> \param converged ...
99 : !> \param sweeps ...
100 : ! **************************************************************************************************
101 210 : 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 30 : CALL timeset(routineN, handle)
121 :
122 30 : NULLIFY (context, para_env, fm_struct_k_k)
123 :
124 30 : 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 30 : para_env=para_env, context=context)
129 : CALL cp_fm_struct_create(fm_struct_k_k, para_env=para_env, context=context, &
130 30 : nrow_global=k, ncol_global=k)
131 30 : CALL cp_fm_create(CTmp, C%matrix_struct)
132 30 : CALL cp_fm_create(U, fm_struct_k_k)
133 30 : CALL cp_fm_create(G, fm_struct_k_k)
134 30 : CALL cp_fm_create(Gp1, fm_struct_k_k)
135 30 : CALL cp_fm_create(Gp2, fm_struct_k_k)
136 : !
137 : ! printing
138 30 : IF (output_unit > 0) THEN
139 15 : WRITE (output_unit, '(1X)')
140 15 : WRITE (output_unit, '(2X,A)') '-----------------------------------------------------------------------------'
141 15 : WRITE (output_unit, '(A,I5)') ' Nbr iterations =', iterations
142 15 : WRITE (output_unit, '(A,E10.2)') ' eps convergence =', eps
143 15 : WRITE (output_unit, '(A,I5)') ' Max Taylor order =', taylor_order
144 15 : WRITE (output_unit, '(A,E10.2)') ' f2 eps =', f2_eps
145 15 : WRITE (output_unit, '(A,E10.2)') ' alpha =', alpha
146 15 : WRITE (output_unit, '(A)') ' iteration approx_l1_norm g_norm rel_err'
147 : END IF
148 : !
149 30 : f2old = 0.0_dp
150 30 : converged = .FALSE.
151 : !
152 : ! Start the steepest descent
153 1480 : DO istep = 1, iterations
154 : !
155 : !-------------------------------------------------------------------
156 : ! compute f_2
157 : ! f_2(x)=(x^2+eps)^1/2
158 1470 : f2 = 0.0_dp
159 22748 : DO p = 1, ncol_local ! p
160 1100956 : DO i = 1, nrow_local ! i
161 1099486 : f2 = f2 + SQRT(C%local_data(i, p)**2 + f2_eps)
162 : END DO
163 : END DO
164 1470 : 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 22748 : DO p = 1, ncol_local ! p
169 1100956 : DO i = 1, nrow_local ! i
170 1099486 : 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 1470 : CALL parallel_gemm('T', 'N', k, k, n, 1.0_dp, CTmp, C, 0.0_dp, G)
174 : ! antisymmetrize
175 1470 : CALL cp_fm_transpose(G, U)
176 1470 : CALL cp_fm_scale_and_add(-0.5_dp, G, 0.5_dp, U)
177 : !
178 : !-------------------------------------------------------------------
179 : !
180 1470 : gnorm = cp_fm_frobenius_norm(G)
181 : !
182 : ! rescale for steepest descent
183 1470 : CALL cp_fm_scale(-alpha, G)
184 : !
185 : ! compute unitary transform
186 : ! zeroth order
187 1470 : CALL cp_fm_set_all(U, 0.0_dp, 1.0_dp)
188 : ! first order
189 1470 : expfactor = 1.0_dp
190 1470 : CALL cp_fm_scale_and_add(1.0_dp, U, expfactor, G)
191 1470 : tnorm = cp_fm_frobenius_norm(G)
192 1470 : IF (tnorm > 1.0E-10_dp) THEN
193 : ! other orders
194 1470 : CALL cp_fm_to_fm(G, Gp1)
195 4866 : DO i = 2, taylor_order
196 : ! new power of G
197 4866 : CALL parallel_gemm('N', 'N', k, k, k, 1.0_dp, G, Gp1, 0.0_dp, Gp2)
198 4866 : CALL cp_fm_to_fm(Gp2, Gp1)
199 : ! add to the taylor expansion so far
200 4866 : expfactor = expfactor/REAL(i, KIND=dp)
201 4866 : CALL cp_fm_scale_and_add(1.0_dp, U, expfactor, Gp1)
202 4866 : tnorm = cp_fm_frobenius_norm(Gp1)
203 4866 : IF (tnorm*expfactor < 1.0E-10_dp) EXIT
204 : END DO
205 : END IF
206 : !
207 : ! incrementaly rotate the MOs
208 1470 : CALL parallel_gemm('N', 'N', n, k, k, 1.0_dp, C, U, 0.0_dp, CTmp)
209 1470 : CALL cp_fm_to_fm(CTmp, C)
210 : !
211 : ! printing
212 1470 : IF (output_unit > 0) THEN
213 735 : 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 1470 : sweeps = istep
218 1470 : IF (ABS((f2 - f2old)/f2) <= eps .AND. istep > 1) THEN
219 20 : converged = .TRUE.
220 20 : EXIT
221 : END IF
222 1460 : 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 30 : IF (output_unit > 0) WRITE (output_unit, '(A,E16.10)') ' sparseness function f2 = ', f2
229 : ! deallocate
230 30 : CALL cp_fm_struct_release(fm_struct_k_k)
231 30 : CALL cp_fm_release(CTmp)
232 30 : CALL cp_fm_release(U)
233 30 : CALL cp_fm_release(G)
234 30 : CALL cp_fm_release(Gp1)
235 30 : CALL cp_fm_release(Gp2)
236 :
237 30 : CALL timestop(handle)
238 :
239 30 : END SUBROUTINE approx_l1_norm_sd
240 : ! **************************************************************************************************
241 : !> \brief ...
242 : !> \param cell ...
243 : !> \param weights ...
244 : ! **************************************************************************************************
245 504 : 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 504 : CPASSERT(ASSOCIATED(cell))
253 :
254 504 : metric = 0.0_dp
255 504 : CALL dgemm('T', 'N', 3, 3, 3, 1._dp, cell%hmat(:, :), 3, cell%hmat(:, :), 3, 0.0_dp, metric(:, :), 3)
256 :
257 504 : weights(1) = METRIC(1, 1) - METRIC(1, 2) - METRIC(1, 3)
258 504 : weights(2) = METRIC(2, 2) - METRIC(1, 2) - METRIC(2, 3)
259 504 : weights(3) = METRIC(3, 3) - METRIC(1, 3) - METRIC(2, 3)
260 504 : weights(4) = METRIC(1, 2)
261 504 : weights(5) = METRIC(1, 3)
262 504 : weights(6) = METRIC(2, 3)
263 :
264 504 : 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 390 : 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 390 : IF (para_env%num_pe == 1) THEN
297 : CALL jacobi_rotations_serial(weights, zij, vectors, max_iter, eps_localization, &
298 0 : sweeps, out_each, restricted=restricted)
299 : ELSE
300 : CALL jacobi_rot_para(weights, zij, vectors, para_env, max_iter, eps_localization, &
301 390 : sweeps, out_each, target_time, start_time, restricted=restricted)
302 : END IF
303 :
304 390 : 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 0 : 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 0 : CALL timeset(routineN, handle)
339 :
340 0 : dim2 = SIZE(zij, 2)
341 0 : ALLOCATE (c_zij(dim2))
342 : NULLIFY (mii, mij, mjj)
343 0 : ALLOCATE (mii(dim2), mij(dim2), mjj(dim2))
344 :
345 0 : CALL cp_fm_create(rmat, zij(1, 1)%matrix_struct)
346 0 : CALL cp_fm_set_all(rmat, 0._dp, 1._dp)
347 :
348 0 : CALL cp_cfm_create(c_rmat, zij(1, 1)%matrix_struct)
349 0 : CALL cp_cfm_set_all(c_rmat, (0._dp, 0._dp), (1._dp, 0._dp))
350 0 : DO idim = 1, dim2
351 0 : 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 0 : zij(2, idim)%local_data, dp)
354 : END DO
355 :
356 0 : CALL cp_fm_get_info(rmat, nrow_global=nstate)
357 0 : tolerance = 1.0e10_dp
358 :
359 0 : sweeps = 0
360 0 : unit_nr = -1
361 0 : IF (rmat%matrix_struct%para_env%is_source()) THEN
362 0 : unit_nr = cp_logger_get_default_unit_nr()
363 0 : WRITE (unit_nr, '(T4,A )') " Localization by iterative Jacobi rotation"
364 : END IF
365 :
366 0 : IF (restricted > 0) THEN
367 0 : unit_nr = cp_logger_get_default_unit_nr()
368 0 : WRITE (unit_nr, '(T4,A,I2,A )') "JACOBI: for the ROKS method, the last ", restricted, " orbitals DO NOT ROTATE"
369 0 : nstate = nstate - restricted
370 : END IF
371 :
372 : ! do jacobi sweeps until converged
373 0 : DO WHILE (tolerance >= eps_localization .AND. sweeps < max_iter)
374 0 : sweeps = sweeps + 1
375 0 : t1 = m_walltime()
376 :
377 0 : DO istate = 1, nstate
378 0 : DO jstate = istate + 1, nstate
379 0 : DO idim = 1, dim2
380 0 : CALL cp_cfm_get_element(c_zij(idim), istate, istate, mii(idim))
381 0 : CALL cp_cfm_get_element(c_zij(idim), istate, jstate, mij(idim))
382 0 : CALL cp_cfm_get_element(c_zij(idim), jstate, jstate, mjj(idim))
383 : END DO
384 0 : CALL get_angle(mii, mjj, mij, weights, theta)
385 0 : st = SIN(theta)
386 0 : ct = COS(theta)
387 0 : CALL rotate_zij(istate, jstate, st, ct, c_zij)
388 :
389 0 : CALL rotate_rmat(istate, jstate, st, ct, c_rmat)
390 : END DO
391 : END DO
392 :
393 0 : CALL check_tolerance(c_zij, weights, tolerance)
394 :
395 0 : t2 = m_walltime()
396 0 : 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 0 : "Iteration:", sweeps, "Tolerance:", tolerance, "Time:", t2 - t1
399 0 : CALL m_flush(unit_nr)
400 : END IF
401 :
402 : END DO
403 :
404 0 : DO idim = 1, dim2
405 0 : zij(1, idim)%local_data = REAL(c_zij(idim)%local_data, dp)
406 0 : zij(2, idim)%local_data = AIMAG(c_zij(idim)%local_data)
407 0 : CALL cp_cfm_release(c_zij(idim))
408 : END DO
409 0 : DEALLOCATE (c_zij)
410 0 : DEALLOCATE (mii, mij, mjj)
411 0 : rmat%local_data = REAL(c_rmat%local_data, dp)
412 :
413 0 : CALL rotate_orbitals(rmat, vectors)
414 :
415 0 : CALL cp_cfm_release(c_rmat)
416 0 : CALL cp_fm_release(rmat)
417 :
418 0 : CALL timestop(handle)
419 :
420 0 : 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 0 : SUBROUTINE jacobi_rotations_serial_1(weights, c_zij, max_iter, c_rmat, eps_localization, &
435 0 : 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 0 : CALL timeset(routineN, handle)
459 :
460 0 : dim2 = SIZE(c_zij)
461 : NULLIFY (mii, mij, mjj)
462 0 : ALLOCATE (mii(dim2), mij(dim2), mjj(dim2))
463 :
464 0 : ALLOCATE (c_zij_local(dim2))
465 0 : CALL cp_cfm_create(c_rmat_local, c_rmat%matrix_struct)
466 0 : CALL cp_cfm_set_all(c_rmat_local, (0.0_dp, 0.0_dp), (1.0_dp, 0.0_dp))
467 0 : DO idim = 1, dim2
468 0 : CALL cp_cfm_create(c_zij_local(idim), c_zij(idim)%matrix_struct)
469 0 : c_zij_local(idim)%local_data = c_zij(idim)%local_data
470 : END DO
471 :
472 0 : CALL cp_cfm_get_info(c_rmat_local, nrow_global=nstate)
473 0 : tolerance = 1.0e10_dp
474 :
475 0 : IF (PRESENT(grad_final)) CALL cp_fm_set_all(grad_final, 0.0_dp)
476 :
477 0 : sweeps = 0
478 0 : IF (PRESENT(out_each)) THEN
479 0 : unit_nr = -1
480 0 : IF (c_rmat_local%matrix_struct%para_env%is_source()) THEN
481 0 : unit_nr = cp_logger_get_default_unit_nr()
482 : END IF
483 0 : alpha = 0.0_dp
484 0 : DO idim = 1, dim2
485 0 : alpha = alpha + weights(idim)
486 : END DO
487 : END IF
488 :
489 : ! do jacobi sweeps until converged
490 0 : DO WHILE (sweeps < max_iter)
491 0 : sweeps = sweeps + 1
492 0 : IF (PRESENT(eps_localization)) THEN
493 0 : IF (tolerance < eps_localization) EXIT
494 : END IF
495 0 : IF (PRESENT(out_each)) t1 = m_walltime()
496 :
497 0 : DO istate = 1, nstate
498 0 : DO jstate = istate + 1, nstate
499 0 : DO idim = 1, dim2
500 0 : CALL cp_cfm_get_element(c_zij_local(idim), istate, istate, mii(idim))
501 0 : CALL cp_cfm_get_element(c_zij_local(idim), istate, jstate, mij(idim))
502 0 : CALL cp_cfm_get_element(c_zij_local(idim), jstate, jstate, mjj(idim))
503 : END DO
504 0 : CALL get_angle(mii, mjj, mij, weights, theta)
505 0 : st = SIN(theta)
506 0 : ct = COS(theta)
507 0 : CALL rotate_zij(istate, jstate, st, ct, c_zij_local)
508 :
509 0 : CALL rotate_rmat(istate, jstate, st, ct, c_rmat_local)
510 : END DO
511 : END DO
512 :
513 0 : IF (PRESENT(grad_final)) THEN
514 0 : CALL check_tolerance(c_zij_local, weights, tolerance, grad=grad_final)
515 : ELSE
516 0 : CALL check_tolerance(c_zij_local, weights, tolerance)
517 : END IF
518 0 : IF (PRESENT(tol_out)) tol_out = tolerance
519 :
520 0 : IF (PRESENT(out_each)) THEN
521 0 : t2 = m_walltime()
522 0 : IF (unit_nr > 0 .AND. MODULO(sweeps, out_each) == 0) THEN
523 0 : sum_spread_ii = 0.0_dp
524 0 : DO istate = 1, nstate
525 : spread_ii = 0.0_dp
526 0 : DO idim = 1, dim2
527 0 : CALL cp_cfm_get_element(c_zij_local(idim), istate, istate, mzii)
528 : spread_ii = spread_ii + weights(idim)* &
529 0 : ABS(mzii)**2/twopi/twopi
530 : END DO
531 0 : sum_spread_ii = sum_spread_ii + spread_ii
532 : END DO
533 0 : sum_spread_ii = alpha*nstate/twopi/twopi - sum_spread_ii
534 0 : avg_spread_ii = sum_spread_ii/nstate
535 : WRITE (unit_nr, '(T4,A,T26,A,T48,A,T64,A)') &
536 0 : "Iteration", "Avg. Spread_ii", "Tolerance", "Time"
537 : WRITE (unit_nr, '(T4,I7,T20,F20.10,T45,E12.4,T60,F8.3)') &
538 0 : sweeps, avg_spread_ii, tolerance, t2 - t1
539 0 : CALL m_flush(unit_nr)
540 : END IF
541 0 : IF (PRESENT(jsweeps)) jsweeps = sweeps
542 : END IF
543 :
544 : END DO
545 :
546 0 : IF (PRESENT(c_zij_out)) THEN
547 0 : DO idim = 1, dim2
548 0 : CALL cp_cfm_to_cfm(c_zij_local(idim), c_zij_out(idim))
549 : END DO
550 : END IF
551 0 : CALL cp_cfm_to_cfm(c_rmat_local, c_rmat)
552 :
553 0 : DEALLOCATE (mii, mij, mjj)
554 0 : DO idim = 1, dim2
555 0 : CALL cp_cfm_release(c_zij_local(idim))
556 : END DO
557 0 : DEALLOCATE (c_zij_local)
558 0 : CALL cp_cfm_release(c_rmat_local)
559 :
560 0 : CALL timestop(handle)
561 :
562 0 : 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 2 : 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 2 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: tmp_cmat
600 2 : COMPLEX(KIND=dp), DIMENSION(:), POINTER :: arr_zii
601 2 : 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 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: tmp_mat, tmp_mat_1
611 : REAL(KIND=dp), DIMENSION(50) :: energy, pos
612 2 : 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 2 : 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 2 : CALL timeset(routineN, handle)
623 :
624 2 : dim2 = SIZE(zij, 2)
625 2 : NULLIFY (context)
626 2 : NULLIFY (matrix_zii, arr_zii)
627 2 : NULLIFY (tmp_fm_struct)
628 2 : NULLIFY (tmp_arr)
629 :
630 12 : ALLOCATE (c_zij(dim2))
631 :
632 2 : CALL cp_fm_get_info(zij(1, 1), nrow_global=nstate)
633 :
634 6 : ALLOCATE (sum_spread(nstate))
635 8 : ALLOCATE (matrix_zii(nstate, dim2))
636 266 : matrix_zii = czero
637 2 : sum_spread = 0.0_dp
638 :
639 2 : alpha = 0.0_dp
640 8 : DO idim = 1, dim2
641 6 : alpha = alpha + weights(idim)
642 6 : 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 5813 : zij(2, idim)%local_data, dp)
645 : END DO
646 :
647 10 : ALLOCATE (zij_0(dim2))
648 :
649 2 : CALL cp_cfm_create(U, zij(1, 1)%matrix_struct)
650 2 : CALL cp_fm_create(matrix_U, zij(1, 1)%matrix_struct)
651 :
652 2 : CALL cp_cfm_set_all(U, czero, cone)
653 2 : CALL cp_fm_set_all(matrix_U, 0.0_dp, 1.0_dp)
654 :
655 2 : CALL cp_fm_get_info(vectors, nrow_global=nao)
656 2 : IF (nextra > 0) THEN
657 2 : IF (PRESENT(mos_guess)) THEN
658 2 : do_cinit_random = .FALSE.
659 2 : do_cinit_mo = .TRUE.
660 2 : CALL cp_fm_get_info(mos_guess, ncol_global=ndummy)
661 : ELSE
662 0 : do_cinit_random = .TRUE.
663 0 : do_cinit_mo = .FALSE.
664 0 : 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 2 : icinit = 2
672 2 : do_U_guess_mo = .TRUE.
673 : END IF
674 :
675 2 : nocc = nstate - nextra
676 2 : northo = nmo - nocc
677 2 : norextra = nmo - nstate
678 2 : CALL cp_fm_struct_get(zij(1, 1)%matrix_struct, context=context)
679 :
680 8 : ALLOCATE (tmp_cmat(nstate, nstate))
681 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nmo, ncol_global=nmo, &
682 2 : para_env=para_env, context=context)
683 8 : DO idim = 1, dim2
684 6 : CALL cp_cfm_create(zij_0(idim), tmp_fm_struct)
685 6 : CALL cp_cfm_set_all(zij_0(idim), czero, cone)
686 6 : CALL cp_cfm_get_submatrix(c_zij(idim), tmp_cmat)
687 8 : CALL cp_cfm_set_submatrix(zij_0(idim), tmp_cmat)
688 : END DO
689 2 : CALL cp_fm_struct_release(tmp_fm_struct)
690 2 : DEALLOCATE (tmp_cmat)
691 :
692 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nmo, ncol_global=nstate, &
693 2 : para_env=para_env, context=context)
694 2 : CALL cp_cfm_create(V, tmp_fm_struct)
695 2 : CALL cp_fm_create(matrix_V, tmp_fm_struct)
696 2 : CALL cp_cfm_create(zdiag, tmp_fm_struct)
697 2 : CALL cp_fm_create(rmat, tmp_fm_struct)
698 2 : CALL cp_fm_struct_release(tmp_fm_struct)
699 2 : CALL cp_cfm_set_all(V, czero, cone)
700 2 : 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 2 : para_env=para_env, context=context)
704 2 : CALL cp_fm_create(matrix_V_all, tmp_fm_struct)
705 2 : CALL cp_fm_struct_release(tmp_fm_struct)
706 2 : CALL cp_fm_set_all(matrix_V_all, 0._dp, 1._dp)
707 :
708 6 : ALLOCATE (arr_zii(nstate))
709 :
710 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=northo, ncol_global=nextra, &
711 2 : para_env=para_env, context=context)
712 2 : CALL cp_cfm_create(c_tilde, tmp_fm_struct)
713 2 : CALL cp_cfm_create(grad_ctilde, tmp_fm_struct)
714 2 : CALL cp_cfm_create(Gct_old, tmp_fm_struct)
715 2 : CALL cp_cfm_create(skc, tmp_fm_struct)
716 2 : CALL cp_fm_struct_release(tmp_fm_struct)
717 2 : CALL cp_cfm_set_all(c_tilde, czero)
718 2 : CALL cp_cfm_set_all(Gct_old, czero)
719 2 : CALL cp_cfm_set_all(skc, czero)
720 :
721 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=northo, ncol_global=nstate, &
722 2 : para_env=para_env, context=context)
723 2 : CALL cp_cfm_create(VL, tmp_fm_struct)
724 2 : CALL cp_cfm_set_all(VL, czero)
725 2 : 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 2 : para_env=para_env, context=context)
729 2 : CALL cp_fm_create(id_nextra, tmp_fm_struct)
730 2 : CALL cp_cfm_create(ctrans_lambda, tmp_fm_struct)
731 2 : CALL cp_fm_struct_release(tmp_fm_struct)
732 2 : CALL cp_cfm_set_all(ctrans_lambda, czero)
733 2 : 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 2 : para_env=para_env, context=context)
737 2 : CALL cp_cfm_create(UL, tmp_fm_struct)
738 2 : CALL cp_fm_struct_release(tmp_fm_struct)
739 2 : CALL cp_cfm_set_all(UL, czero)
740 :
741 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, ncol_global=nmo, &
742 2 : para_env=para_env, context=context)
743 2 : CALL cp_fm_create(vectors_all, tmp_fm_struct)
744 2 : CALL cp_fm_struct_release(tmp_fm_struct)
745 8 : ALLOCATE (tmp_mat(nao, nstate))
746 2 : CALL cp_fm_get_submatrix(vectors, tmp_mat)
747 2 : CALL cp_fm_set_submatrix(vectors_all, tmp_mat, 1, 1, nao, nstate)
748 2 : DEALLOCATE (tmp_mat)
749 8 : ALLOCATE (tmp_mat(nao, norextra))
750 2 : CALL cp_fm_get_submatrix(vectors_2, tmp_mat)
751 2 : CALL cp_fm_set_submatrix(vectors_all, tmp_mat, 1, nstate + 1, nao, norextra)
752 2 : DEALLOCATE (tmp_mat)
753 :
754 : ! initialize c_tilde
755 14 : SELECT CASE (icinit)
756 : CASE (1) ! random coefficients
757 : !WRITE (*, *) "RANDOM INITIAL GUESS FOR C"
758 0 : CALL cp_fm_create(tmp_fm, c_tilde%matrix_struct)
759 0 : CALL cp_fm_init_random(tmp_fm, nextra)
760 0 : CALL ortho_vectors(tmp_fm)
761 0 : c_tilde%local_data = tmp_fm%local_data
762 0 : CALL cp_fm_release(tmp_fm)
763 0 : ALLOCATE (tmp_cmat(northo, nextra))
764 0 : CALL cp_cfm_get_submatrix(c_tilde, tmp_cmat)
765 0 : CALL cp_cfm_set_submatrix(V, tmp_cmat, nocc + 1, nocc + 1, northo, nextra)
766 0 : DEALLOCATE (tmp_cmat)
767 : CASE (2) ! MO based coeffs
768 2 : CALL parallel_gemm("T", "N", nmo, ndummy, nao, 1.0_dp, vectors_all, mos_guess, 0.0_dp, matrix_V_all)
769 6 : ALLOCATE (tmp_arr(nmo))
770 8 : ALLOCATE (tmp_mat(nmo, ndummy))
771 8 : ALLOCATE (tmp_mat_1(nmo, nstate))
772 : ! normalize matrix_V_all
773 2 : CALL cp_fm_get_submatrix(matrix_V_all, tmp_mat)
774 88 : DO istate = 1, ndummy
775 4558 : tmp_arr(:) = tmp_mat(:, istate)
776 4558 : norm = NORM2(tmp_arr)
777 4558 : tmp_arr(:) = tmp_arr(:)/norm
778 4560 : tmp_mat(:, istate) = tmp_arr(:)
779 : END DO
780 2 : CALL cp_fm_set_submatrix(matrix_V_all, tmp_mat)
781 2 : CALL cp_fm_get_submatrix(matrix_V_all, tmp_mat_1, 1, 1, nmo, nstate)
782 2 : CALL cp_fm_set_submatrix(matrix_V, tmp_mat_1)
783 2 : DEALLOCATE (tmp_arr, tmp_mat, tmp_mat_1)
784 2 : CALL cp_fm_to_cfm(msourcer=matrix_V, mtarget=V)
785 8 : ALLOCATE (tmp_mat(northo, ndummy))
786 8 : ALLOCATE (tmp_mat_1(northo, nextra))
787 2 : CALL cp_fm_get_submatrix(matrix_V_all, tmp_mat, nocc + 1, 1, northo, ndummy)
788 6 : ALLOCATE (tmp_arr(ndummy))
789 88 : tmp_arr = 0.0_dp
790 88 : DO istate = 1, ndummy
791 1034 : tmp_arr(istate) = NORM2(tmp_mat(:, istate))
792 : END DO
793 : ! find edfs
794 6 : DO istate = 1, nextra
795 180 : iloc = MAXLOC(tmp_arr)
796 48 : tmp_mat_1(:, istate) = tmp_mat(:, iloc(1))
797 6 : tmp_arr(iloc(1)) = 0.0_dp
798 : END DO
799 :
800 2 : DEALLOCATE (tmp_arr, tmp_mat)
801 :
802 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=northo, ncol_global=nextra, &
803 2 : para_env=para_env, context=context)
804 2 : CALL cp_fm_create(tmp_fm, tmp_fm_struct)
805 2 : CALL cp_fm_struct_release(tmp_fm_struct)
806 2 : CALL cp_fm_set_submatrix(tmp_fm, tmp_mat_1)
807 2 : DEALLOCATE (tmp_mat_1)
808 2 : CALL ortho_vectors(tmp_fm)
809 2 : CALL cp_fm_to_cfm(msourcer=tmp_fm, mtarget=c_tilde)
810 2 : CALL cp_fm_release(tmp_fm)
811 : ! initialize U
812 2 : IF (do_U_guess_mo) THEN
813 8 : ALLOCATE (tmp_cmat(nocc, nstate))
814 2 : CALL cp_cfm_get_submatrix(V, tmp_cmat, 1, 1, nocc, nstate)
815 2 : CALL cp_cfm_set_submatrix(U, tmp_cmat, 1, 1, nocc, nstate)
816 2 : DEALLOCATE (tmp_cmat)
817 8 : ALLOCATE (tmp_cmat(northo, nstate))
818 2 : CALL cp_cfm_get_submatrix(V, tmp_cmat, nocc + 1, 1, northo, nstate)
819 2 : CALL cp_cfm_set_submatrix(VL, tmp_cmat, 1, 1, northo, nstate)
820 2 : DEALLOCATE (tmp_cmat)
821 2 : CALL parallel_gemm("C", "N", nextra, nstate, northo, cone, c_tilde, VL, czero, UL)
822 8 : ALLOCATE (tmp_cmat(nextra, nstate))
823 2 : CALL cp_cfm_get_submatrix(UL, tmp_cmat, 1, 1, nextra, nstate)
824 2 : CALL cp_cfm_set_submatrix(U, tmp_cmat, nocc + 1, 1, nextra, nstate)
825 2 : DEALLOCATE (tmp_cmat)
826 2 : CALL cp_fm_create(tmp_fm, U%matrix_struct)
827 1937 : tmp_fm%local_data = REAL(U%local_data, KIND=dp)
828 2 : CALL ortho_vectors(tmp_fm)
829 2 : CALL cp_fm_to_cfm(msourcer=tmp_fm, mtarget=U)
830 2 : CALL cp_fm_release(tmp_fm)
831 4 : CALL cp_cfm_to_fm(U, matrix_U)
832 : END IF
833 : ! reevaluate V
834 8 : ALLOCATE (tmp_cmat(nocc, nstate))
835 2 : CALL cp_cfm_get_submatrix(U, tmp_cmat, 1, 1, nocc, nstate)
836 2 : CALL cp_cfm_set_submatrix(V, tmp_cmat, 1, 1, nocc, nstate)
837 2 : DEALLOCATE (tmp_cmat)
838 8 : ALLOCATE (tmp_cmat(nextra, nstate))
839 2 : CALL cp_cfm_get_submatrix(U, tmp_cmat, nocc + 1, 1, nextra, nstate)
840 2 : CALL cp_cfm_set_submatrix(UL, tmp_cmat, 1, 1, nextra, nstate)
841 2 : DEALLOCATE (tmp_cmat)
842 2 : CALL parallel_gemm("N", "N", northo, nstate, nextra, cone, c_tilde, UL, czero, VL)
843 8 : ALLOCATE (tmp_cmat(northo, nstate))
844 2 : CALL cp_cfm_get_submatrix(VL, tmp_cmat)
845 2 : CALL cp_cfm_set_submatrix(V, tmp_cmat, nocc + 1, 1, northo, nstate)
846 6 : DEALLOCATE (tmp_cmat)
847 : END SELECT
848 : ELSE
849 0 : DO idim = 1, dim2
850 0 : CALL cp_cfm_create(zij_0(idim), zij(1, 1)%matrix_struct)
851 0 : CALL cp_cfm_to_cfm(c_zij(idim), zij_0(idim))
852 : END DO
853 0 : CALL cp_fm_create(rmat, zij(1, 1)%matrix_struct)
854 0 : CALL cp_fm_set_all(rmat, 0._dp, 1._dp)
855 : END IF
856 :
857 2 : unit_nr = -1
858 2 : IF (rmat%matrix_struct%para_env%is_source()) THEN
859 1 : unit_nr = cp_logger_get_default_unit_nr()
860 1 : WRITE (unit_nr, '(T4,A )') " Localization by combined Jacobi rotations and Non-Linear Conjugate Gradient"
861 : END IF
862 :
863 2 : norm2_old = 1.0E30_dp
864 2 : ds_min = 1.0_dp
865 2 : new_direction = .TRUE.
866 2 : iter = 0
867 2 : line_searches = 0
868 2 : line_search_count = 0
869 2 : tol = 1.0E+20_dp
870 2 : mintol = 1.0E+10_dp
871 2 : 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 34 : DO WHILE (iter < max_iter)
878 34 : iter = iter + 1
879 : !WRITE(*,*) 'iter = ', iter
880 34 : t1 = m_walltime()
881 :
882 34 : IF (iter > 1) THEN
883 : ! comput U
884 32 : CALL cp_cfm_create(tmp_cfm, zij(1, 1)%matrix_struct)
885 32 : CALL cp_cfm_create(tmp_cfm_2, zij(1, 1)%matrix_struct)
886 32 : IF (para_env%num_pe == 1) THEN
887 0 : CALL jacobi_rotations_serial_1(weights, c_zij, 1, tmp_cfm_2, tol_out=tol)
888 : ELSE
889 32 : CALL jacobi_rot_para_1(weights, c_zij, para_env, 1, tmp_cfm_2, tol_out=tol)
890 : END IF
891 32 : CALL parallel_gemm('N', 'N', nstate, nstate, nstate, cone, U, tmp_cfm_2, czero, tmp_cfm)
892 32 : CALL cp_cfm_to_cfm(tmp_cfm, U)
893 32 : CALL cp_cfm_release(tmp_cfm)
894 32 : CALL cp_cfm_release(tmp_cfm_2)
895 : END IF
896 :
897 34 : IF (nextra > 0) THEN
898 136 : ALLOCATE (tmp_cmat(nextra, nstate))
899 34 : CALL cp_cfm_get_submatrix(U, tmp_cmat, nocc + 1, 1, nextra, nstate)
900 34 : CALL cp_cfm_set_submatrix(UL, tmp_cmat)
901 34 : DEALLOCATE (tmp_cmat)
902 34 : IF (iter > 1) THEN
903 : ! orthonormalize c_tilde
904 32 : CALL cp_fm_create(tmp_fm, c_tilde%matrix_struct)
905 448 : tmp_fm%local_data = REAL(c_tilde%local_data, KIND=dp)
906 32 : CALL ortho_vectors(tmp_fm)
907 32 : CALL cp_fm_to_cfm(msourcer=tmp_fm, mtarget=c_tilde)
908 32 : CALL cp_fm_release(tmp_fm)
909 :
910 128 : ALLOCATE (tmp_cmat(nocc, nstate))
911 32 : CALL cp_cfm_get_submatrix(U, tmp_cmat, 1, 1, nocc, nstate)
912 32 : CALL cp_cfm_set_submatrix(V, tmp_cmat, 1, 1, nocc, nstate)
913 32 : DEALLOCATE (tmp_cmat)
914 32 : CALL parallel_gemm("N", "N", northo, nstate, nextra, cone, c_tilde, UL, czero, VL)
915 128 : ALLOCATE (tmp_cmat(northo, nstate))
916 32 : CALL cp_cfm_get_submatrix(VL, tmp_cmat)
917 32 : CALL cp_cfm_set_submatrix(V, tmp_cmat, nocc + 1, 1, northo, nstate)
918 64 : DEALLOCATE (tmp_cmat)
919 : END IF
920 :
921 : ! reset if new_direction
922 34 : IF (new_direction .AND. MOD(line_searches, 20) == 5) THEN
923 0 : CALL cp_cfm_set_all(skc, czero)
924 0 : CALL cp_cfm_set_all(Gct_old, czero)
925 0 : norm2_old = 1.0E30_dp
926 : END IF
927 :
928 34 : CALL cp_cfm_create(tmp_cfm, V%matrix_struct)
929 34 : CALL cp_cfm_to_cfm(V, tmp_cfm)
930 34 : CALL cp_cfm_create(tmp_cfm_1, V%matrix_struct)
931 68 : ndummy = nmo
932 : ELSE
933 0 : CALL cp_cfm_create(tmp_cfm, zij(1, 1)%matrix_struct)
934 0 : CALL cp_cfm_to_cfm(U, tmp_cfm)
935 0 : CALL cp_cfm_create(tmp_cfm_1, zij(1, 1)%matrix_struct)
936 0 : ndummy = nstate
937 : END IF
938 : ! update z_ij
939 136 : 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 102 : 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 136 : czero, c_zij(idim))
946 : END DO
947 34 : CALL cp_cfm_release(tmp_cfm)
948 34 : CALL cp_cfm_release(tmp_cfm_1)
949 : ! compute spread
950 1496 : DO istate = 1, nstate
951 1462 : spread_ii = 0.0_dp
952 5848 : DO idim = 1, dim2
953 4386 : CALL cp_cfm_get_element(c_zij(idim), istate, istate, mzii)
954 : spread_ii = spread_ii + weights(idim)* &
955 4386 : ABS(mzii)**2/twopi/twopi
956 5848 : matrix_zii(istate, idim) = mzii
957 : END DO
958 : !WRITE(*,*) 'spread_ii', spread_ii
959 1496 : sum_spread(istate) = spread_ii
960 : END DO
961 34 : CALL c_zij(1)%matrix_struct%para_env%sum(spread_ii)
962 34 : spread_sum = accurate_sum(sum_spread)
963 :
964 34 : IF (nextra > 0) THEN
965 : ! update c_tilde
966 34 : CALL cp_cfm_set_all(zdiag, czero)
967 34 : CALL cp_cfm_set_all(grad_ctilde, czero)
968 34 : CALL cp_cfm_create(tmp_cfm, V%matrix_struct)
969 34 : CALL cp_cfm_set_all(tmp_cfm, czero)
970 34 : CALL cp_cfm_create(tmp_cfm_1, V%matrix_struct)
971 34 : CALL cp_cfm_set_all(tmp_cfm_1, czero)
972 136 : ALLOCATE (tmp_cmat(northo, nstate))
973 136 : DO idim = 1, dim2
974 102 : weight = weights(idim)
975 8874 : arr_zii = matrix_zii(:, idim)
976 : ! tmp_cfm = zij_0*V
977 : CALL parallel_gemm("N", "N", nmo, nstate, nmo, cone, &
978 102 : zij_0(idim), V, czero, tmp_cfm)
979 : ! tmp_cfm = tmp_cfm*diag_zij_dagg
980 4488 : 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 102 : U, czero, tmp_cfm_1)
984 102 : CALL cp_cfm_scale(weight, tmp_cfm_1)
985 : ! zdiag = zdiag + tmp_cfm_1'
986 102 : 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 102 : zij_0(idim), V, czero, tmp_cfm)
991 :
992 : ! tmp_cfm = tmp_cfm*diag_zij
993 102 : 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 102 : U, czero, tmp_cfm_1)
997 102 : CALL cp_cfm_scale(weight, tmp_cfm_1)
998 : ! zdiag = zdiag + tmp_cfm_1'
999 136 : CALL cp_cfm_scale_and_add(cone, zdiag, cone, tmp_cfm_1)
1000 : END DO ! idim
1001 34 : CALL cp_cfm_release(tmp_cfm)
1002 34 : CALL cp_cfm_release(tmp_cfm_1)
1003 34 : DEALLOCATE (tmp_cmat)
1004 136 : ALLOCATE (tmp_cmat(northo, nextra))
1005 : CALL cp_cfm_get_submatrix(zdiag, tmp_cmat, nocc + 1, nocc + 1, &
1006 34 : northo, nextra, .FALSE.)
1007 : ! 'grad_ctilde'
1008 34 : CALL cp_cfm_set_submatrix(grad_ctilde, tmp_cmat)
1009 34 : DEALLOCATE (tmp_cmat)
1010 : ! ctrans_lambda = c_tilde_dagg*grad_ctilde
1011 34 : 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 102 : 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 34 : IF (nextra > 0) THEN
1019 : tolc = 0.0_dp
1020 34 : CALL cp_fm_create(tmp_fm, grad_ctilde%matrix_struct)
1021 34 : CALL cp_cfm_to_fm(grad_ctilde, tmp_fm)
1022 34 : CALL cp_fm_maxabsval(tmp_fm, tolc)
1023 34 : CALL cp_fm_release(tmp_fm)
1024 : !WRITE(*,*) 'tolc = ', tolc
1025 34 : tol = tol + tolc
1026 : END IF
1027 : !WRITE(*,*) 'tol = ', tol
1028 :
1029 36 : IF (nextra > 0) THEN
1030 : !WRITE(*,*) 'new_direction: ', new_direction
1031 34 : IF (new_direction) THEN
1032 6 : line_searches = line_searches + 1
1033 6 : IF (mintol > tol) THEN
1034 4 : mintol = tol
1035 4 : miniter = iter
1036 : END IF
1037 :
1038 6 : IF (unit_nr > 0 .AND. MODULO(iter, out_each) == 0) THEN
1039 0 : sum_spread_ii = alpha*nstate/twopi/twopi - spread_sum
1040 0 : avg_spread_ii = sum_spread_ii/nstate
1041 : WRITE (unit_nr, '(T4,A,T26,A,T48,A)') &
1042 0 : "Iteration", "Avg. Spread_ii", "Tolerance"
1043 : WRITE (unit_nr, '(T4,I7,T20,F20.10,T45,E12.4)') &
1044 0 : iter, avg_spread_ii, tol
1045 0 : CALL m_flush(unit_nr)
1046 : END IF
1047 6 : IF (tol < eps_localization) EXIT
1048 :
1049 4 : IF (do_cg) THEN
1050 : cnorm2_Gct = czero
1051 : cnorm2_Gct_cross = czero
1052 4 : CALL cp_cfm_trace(grad_ctilde, Gct_old, cnorm2_Gct_cross)
1053 4 : norm2_Gct_cross = REAL(cnorm2_Gct_cross, KIND=dp)
1054 56 : Gct_old%local_data = grad_ctilde%local_data
1055 4 : CALL cp_cfm_trace(grad_ctilde, Gct_old, cnorm2_Gct)
1056 4 : norm2_Gct = REAL(cnorm2_Gct, KIND=dp)
1057 : ! compute beta_pr
1058 4 : beta_pr = (norm2_Gct - norm2_Gct_cross)/norm2_old
1059 4 : norm2_old = norm2_Gct
1060 4 : beta = MAX(0.0_dp, beta_pr)
1061 : !WRITE(*,*) 'beta = ', beta
1062 : ! compute skc / ska = beta * skc / ska + grad_ctilde / G
1063 4 : CALL cp_cfm_scale(beta, skc)
1064 4 : CALL cp_cfm_scale_and_add(cone, skc, cone, Gct_old)
1065 4 : CALL cp_cfm_trace(skc, Gct_old, cnorm2_Gct_cross)
1066 4 : norm2_Gct_cross = REAL(cnorm2_Gct_cross, KIND=dp)
1067 4 : IF (norm2_Gct_cross <= 0.0_dp) THEN ! back to steepest ascent
1068 0 : CALL cp_cfm_scale_and_add(czero, skc, cone, Gct_old)
1069 : END IF
1070 : ELSE
1071 0 : CALL cp_cfm_scale_and_add(czero, skc, cone, grad_ctilde)
1072 : END IF
1073 : line_search_count = 0
1074 : END IF
1075 :
1076 32 : line_search_count = line_search_count + 1
1077 : !WRITE(*,*) 'line_search_count = ', line_search_count
1078 32 : energy(line_search_count) = spread_sum
1079 :
1080 : ! gold line search
1081 32 : new_direction = .FALSE.
1082 32 : IF (line_search_count == 1) THEN
1083 4 : lsl = 1
1084 4 : lsr = 0
1085 4 : lsm = 1
1086 4 : pos(1) = 0.0_dp
1087 4 : pos(2) = ds_min/gold_sec
1088 4 : ds = pos(2)
1089 : ELSE
1090 28 : IF (line_search_count == 50) THEN
1091 0 : IF (ABS(energy(line_search_count) - energy(line_search_count - 1)) < 1.0E-4_dp) THEN
1092 0 : CPWARN("Line search failed to converge properly")
1093 0 : ds_min = 0.1_dp
1094 0 : new_direction = .TRUE.
1095 : ds = pos(line_search_count)
1096 0 : line_search_count = 0
1097 : ELSE
1098 0 : CPABORT("No. of line searches exceeds 50")
1099 : END IF
1100 : ELSE
1101 28 : IF (lsr == 0) THEN
1102 28 : IF (energy(line_search_count - 1) > energy(line_search_count)) THEN
1103 0 : lsr = line_search_count
1104 0 : pos(line_search_count + 1) = pos(lsm) + (pos(lsr) - pos(lsm))*gold_sec
1105 : ELSE
1106 28 : lsl = lsm
1107 28 : lsm = line_search_count
1108 28 : pos(line_search_count + 1) = pos(line_search_count)/gold_sec
1109 : END IF
1110 : ELSE
1111 0 : IF (pos(line_search_count) < pos(lsm)) THEN
1112 0 : IF (energy(line_search_count) > energy(lsm)) THEN
1113 : lsr = lsm
1114 : lsm = line_search_count
1115 : ELSE
1116 0 : lsl = line_search_count
1117 : END IF
1118 : ELSE
1119 0 : IF (energy(line_search_count) > energy(lsm)) THEN
1120 : lsl = lsm
1121 : lsm = line_search_count
1122 : ELSE
1123 0 : lsr = line_search_count
1124 : END IF
1125 : END IF
1126 0 : IF (pos(lsr) - pos(lsm) > pos(lsm) - pos(lsl)) THEN
1127 0 : pos(line_search_count + 1) = pos(lsm) + gold_sec*(pos(lsr) - pos(lsm))
1128 : ELSE
1129 0 : pos(line_search_count + 1) = pos(lsl) + gold_sec*(pos(lsm) - pos(lsl))
1130 : END IF
1131 0 : IF ((pos(lsr) - pos(lsl)) < 1.0E-3_dp*pos(lsr)) THEN
1132 0 : 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 28 : ds = pos(line_search_count + 1) - pos(line_search_count)
1138 : !WRITE(*,*) 'lsl, lsr, lsm, ds = ', lsl, lsr, lsm, ds
1139 28 : IF ((ABS(ds) < 1.0E-10_dp) .AND. (lsl == 1)) THEN
1140 0 : new_direction = .TRUE.
1141 0 : ds_min = 0.5_dp/alpha
1142 28 : ELSE IF (ABS(ds) > 10.0_dp) THEN
1143 4 : new_direction = .TRUE.
1144 4 : 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 32 : CALL cp_cfm_scale(ds, skc)
1151 32 : CALL cp_cfm_scale_and_add(cone, c_tilde, cone, skc)
1152 : ELSE
1153 0 : IF (mintol > tol) THEN
1154 0 : mintol = tol
1155 0 : miniter = iter
1156 : END IF
1157 0 : IF (unit_nr > 0 .AND. MODULO(iter, out_each) == 0) THEN
1158 0 : sum_spread_ii = alpha*nstate/twopi/twopi - spread_sum
1159 0 : avg_spread_ii = sum_spread_ii/nstate
1160 : WRITE (unit_nr, '(T4,A,T26,A,T48,A)') &
1161 0 : "Iteration", "Avg. Spread_ii", "Tolerance"
1162 : WRITE (unit_nr, '(T4,I7,T20,F20.10,T45,E12.4)') &
1163 0 : iter, avg_spread_ii, tol
1164 0 : CALL m_flush(unit_nr)
1165 : END IF
1166 0 : IF (tol < eps_localization) EXIT
1167 : END IF ! nextra > 0
1168 :
1169 : END DO ! iteration
1170 :
1171 2 : IF ((unit_nr > 0) .AND. (iter == max_iter)) THEN
1172 0 : WRITE (unit_nr, '(T4,A,T4,A)') "Min. Itr.", "Min. Tol."
1173 0 : WRITE (unit_nr, '(T4,I7,T4,E12.4)') miniter, mintol
1174 0 : CALL m_flush(unit_nr)
1175 : END IF
1176 :
1177 2 : CALL cp_cfm_to_fm(U, matrix_U)
1178 :
1179 2 : IF (nextra > 0) THEN
1180 2324 : rmat%local_data = REAL(V%local_data, KIND=dp)
1181 2 : CALL rotate_orbitals_edf(rmat, vectors_all, vectors)
1182 :
1183 2 : CALL cp_cfm_release(c_tilde)
1184 2 : CALL cp_cfm_release(grad_ctilde)
1185 2 : CALL cp_cfm_release(Gct_old)
1186 2 : CALL cp_cfm_release(skc)
1187 2 : CALL cp_cfm_release(UL)
1188 2 : CALL cp_cfm_release(zdiag)
1189 2 : CALL cp_cfm_release(ctrans_lambda)
1190 2 : CALL cp_fm_release(id_nextra)
1191 2 : CALL cp_fm_release(vectors_all)
1192 2 : CALL cp_cfm_release(V)
1193 2 : CALL cp_fm_release(matrix_V)
1194 2 : CALL cp_fm_release(matrix_V_all)
1195 2 : CALL cp_cfm_release(VL)
1196 2 : DEALLOCATE (arr_zii)
1197 : ELSE
1198 0 : rmat%local_data = matrix_U%local_data
1199 0 : CALL rotate_orbitals(rmat, vectors)
1200 : END IF
1201 8 : DO idim = 1, dim2
1202 8 : CALL cp_cfm_release(zij_0(idim))
1203 : END DO
1204 2 : DEALLOCATE (zij_0)
1205 :
1206 8 : DO idim = 1, dim2
1207 5811 : zij(1, idim)%local_data = REAL(c_zij(idim)%local_data, dp)
1208 5811 : zij(2, idim)%local_data = AIMAG(c_zij(idim)%local_data)
1209 8 : CALL cp_cfm_release(c_zij(idim))
1210 : END DO
1211 2 : DEALLOCATE (c_zij)
1212 2 : CALL cp_fm_release(rmat)
1213 2 : CALL cp_cfm_release(U)
1214 2 : CALL cp_fm_release(matrix_U)
1215 2 : DEALLOCATE (matrix_zii, sum_spread)
1216 :
1217 2 : CALL timestop(handle)
1218 :
1219 10 : END SUBROUTINE jacobi_cg_edf_ls
1220 :
1221 : ! **************************************************************************************************
1222 : !> \brief ...
1223 : !> \param vmatrix ...
1224 : ! **************************************************************************************************
1225 108 : 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 36 : CALL timeset(routineN, handle)
1236 :
1237 36 : NULLIFY (fm_struct_tmp)
1238 :
1239 36 : 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 36 : context=vmatrix%matrix_struct%context)
1244 36 : CALL cp_fm_create(overlap_vv, fm_struct_tmp, "overlap_vv")
1245 36 : CALL cp_fm_struct_release(fm_struct_tmp)
1246 :
1247 36 : CALL parallel_gemm('T', 'N', ncol, ncol, n, 1.0_dp, vmatrix, vmatrix, 0.0_dp, overlap_vv)
1248 36 : CALL cp_fm_cholesky_decompose(overlap_vv)
1249 36 : CALL cp_fm_triangular_multiply(overlap_vv, vmatrix, n_cols=ncol, side='R', invert_tr=.TRUE.)
1250 :
1251 36 : CALL cp_fm_release(overlap_vv)
1252 :
1253 36 : CALL timestop(handle)
1254 :
1255 36 : 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 0 : 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 0 : DO id = 1, SIZE(zij, 1)
1275 0 : CALL cp_cfm_rot_cols(zij(id), istate, jstate, ct, st)
1276 0 : CALL cp_cfm_rot_rows(zij(id), istate, jstate, ct, st)
1277 : END DO
1278 :
1279 0 : END SUBROUTINE rotate_zij
1280 : ! **************************************************************************************************
1281 : !> \brief ...
1282 : !> \param istate ...
1283 : !> \param jstate ...
1284 : !> \param st ...
1285 : !> \param ct ...
1286 : !> \param rmat ...
1287 : ! **************************************************************************************************
1288 0 : 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 0 : CALL cp_cfm_rot_cols(rmat, istate, jstate, ct, st)
1294 :
1295 0 : 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 2342024 : 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 2342024 : a12 = 0.0_dp
1317 2342024 : b12 = 0.0_dp
1318 2342024 : dim_m = SIZE(mii)
1319 9509066 : DO idim = 1, dim_m
1320 7167042 : z11 = mii(idim)
1321 7167042 : z22 = mjj(idim)
1322 7167042 : z12 = mij(idim)
1323 7167042 : a12 = a12 + weights(idim)*REAL(CONJG(z12)*(z11 - z22), KIND=dp)
1324 : b12 = b12 + weights(idim)*REAL((z12*CONJG(z12) - &
1325 9509066 : 0.25_dp*(z11 - z22)*(CONJG(z11) - CONJG(z22))), KIND=dp)
1326 : END DO
1327 2342024 : IF (ABS(b12) > 1.e-10_dp) THEN
1328 2342024 : ratio = -a12/b12
1329 2342024 : theta = 0.25_dp*ATAN(ratio)
1330 0 : ELSE IF (ABS(b12) < 1.e-10_dp) THEN
1331 0 : b12 = 0.0_dp
1332 0 : theta = 0.0_dp
1333 : ELSE
1334 0 : theta = 0.25_dp*pi
1335 : END IF
1336 2342024 : IF (PRESENT(grad_ij)) theta = theta + step*grad_ij
1337 : ! Check second derivative info
1338 2342024 : d2 = a12*SIN(4._dp*theta) - b12*COS(4._dp*theta)
1339 2342024 : IF (d2 <= 0._dp) THEN ! go to the maximum, not the minimum
1340 3322 : IF (theta > 0.0_dp) THEN ! make theta as small as possible
1341 1629 : theta = theta - 0.25_dp*pi
1342 : ELSE
1343 1693 : theta = theta + 0.25_dp*pi
1344 : END IF
1345 : END IF
1346 2342024 : END SUBROUTINE get_angle
1347 : ! **************************************************************************************************
1348 : !> \brief ...
1349 : !> \param zij ...
1350 : !> \param weights ...
1351 : !> \param tolerance ...
1352 : !> \param grad ...
1353 : ! **************************************************************************************************
1354 104 : 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 104 : CALL timeset(routineN, handle)
1366 :
1367 : ! compute gradient at t=0
1368 :
1369 104 : CALL cp_fm_create(force, zij(1)%matrix_struct)
1370 104 : CALL cp_fm_set_all(force, 0._dp)
1371 104 : CALL grad_at_0(zij, weights, force)
1372 104 : CALL cp_fm_maxabsval(force, tolerance)
1373 104 : IF (PRESENT(grad)) CALL cp_fm_to_fm(force, grad)
1374 104 : CALL cp_fm_release(force)
1375 :
1376 104 : CALL timestop(handle)
1377 :
1378 104 : END SUBROUTINE check_tolerance
1379 :
1380 : ! **************************************************************************************************
1381 : !> \brief ...
1382 : !> \param rmat ...
1383 : !> \param vectors ...
1384 : ! **************************************************************************************************
1385 1128 : 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 564 : CALL cp_fm_create(wf, vectors%matrix_struct)
1392 564 : CALL cp_fm_get_info(vectors, nrow_global=n, ncol_global=k)
1393 564 : CALL parallel_gemm("N", "N", n, k, k, 1.0_dp, vectors, rmat, 0.0_dp, wf)
1394 564 : CALL cp_fm_to_fm(wf, vectors)
1395 564 : CALL cp_fm_release(wf)
1396 564 : END SUBROUTINE rotate_orbitals
1397 :
1398 : ! **************************************************************************************************
1399 : !> \brief ...
1400 : !> \param rmat ...
1401 : !> \param vectors ...
1402 : ! **************************************************************************************************
1403 24 : 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 12 : CALL cp_cfm_create(wf, vectors%matrix_struct)
1410 12 : CALL cp_cfm_get_info(vectors, nrow_global=n, ncol_global=k)
1411 12 : CALL parallel_gemm("N", "N", n, k, k, z_one, vectors, rmat, z_zero, wf)
1412 12 : CALL cp_cfm_to_cfm(wf, vectors)
1413 12 : CALL cp_cfm_release(wf)
1414 12 : END SUBROUTINE rotate_orbitals_cfm
1415 :
1416 : ! **************************************************************************************************
1417 : !> \brief ...
1418 : !> \param rmat ...
1419 : !> \param vec_all ...
1420 : !> \param vectors ...
1421 : ! **************************************************************************************************
1422 6 : 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 2 : CALL cp_fm_create(wf, vectors%matrix_struct)
1429 2 : CALL cp_fm_get_info(vec_all, nrow_global=n, ncol_global=k)
1430 2 : CALL cp_fm_get_info(rmat, ncol_global=l)
1431 :
1432 2 : CALL parallel_gemm("N", "N", n, l, k, 1.0_dp, vec_all, rmat, 0.0_dp, wf)
1433 2 : CALL cp_fm_to_fm(wf, vectors)
1434 2 : CALL cp_fm_release(wf)
1435 2 : END SUBROUTINE rotate_orbitals_edf
1436 : ! **************************************************************************************************
1437 : !> \brief ...
1438 : !> \param diag ...
1439 : !> \param weights ...
1440 : !> \param matrix ...
1441 : !> \param ndim ...
1442 : ! **************************************************************************************************
1443 154 : 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 154 : 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 154 : row_indices=row_indices, col_indices=col_indices)
1458 :
1459 308 : DO istate = 1, nrow_local
1460 616 : DO jstate = 1, ncol_local
1461 : ! get real and imaginary parts
1462 308 : gradsq_ij = 0.0_dp
1463 2156 : DO idim = 1, ndim
1464 1848 : zii = diag(row_indices(istate), idim)
1465 1848 : zjj = diag(col_indices(jstate), idim)
1466 : gradsq_ij = gradsq_ij + weights(idim)* &
1467 2156 : 4.0_dp*REAL((CONJG(zii)*zii + CONJG(zjj)*zjj), KIND=dp)
1468 : END DO
1469 462 : matrix%local_data(istate, jstate) = gradsq_ij
1470 : END DO
1471 : END DO
1472 154 : END SUBROUTINE gradsq_at_0
1473 : ! **************************************************************************************************
1474 : !> \brief ...
1475 : !> \param matrix_p ...
1476 : !> \param weights ...
1477 : !> \param matrix ...
1478 : ! **************************************************************************************************
1479 104 : 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 104 : COMPLEX(KIND=dp), DIMENSION(:, :), POINTER :: diag
1486 : INTEGER :: dim_m, idim, istate, jstate, ncol_local, &
1487 : nrow_global, nrow_local
1488 104 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1489 : REAL(KIND=dp) :: grad_ij
1490 :
1491 104 : NULLIFY (diag)
1492 : CALL cp_fm_get_info(matrix, nrow_local=nrow_local, &
1493 : ncol_local=ncol_local, nrow_global=nrow_global, &
1494 104 : row_indices=row_indices, col_indices=col_indices)
1495 104 : dim_m = SIZE(matrix_p, 1)
1496 416 : ALLOCATE (diag(nrow_global, dim_m))
1497 :
1498 728 : DO idim = 1, dim_m
1499 23192 : DO istate = 1, nrow_global
1500 23088 : CALL cp_cfm_get_element(matrix_p(idim), istate, istate, diag(istate, idim))
1501 : END DO
1502 : END DO
1503 :
1504 1976 : DO istate = 1, nrow_local
1505 69368 : DO jstate = 1, ncol_local
1506 : ! get real and imaginary parts
1507 : grad_ij = 0.0_dp
1508 471744 : DO idim = 1, dim_m
1509 404352 : zii = diag(row_indices(istate), idim)
1510 404352 : zjj = diag(col_indices(jstate), idim)
1511 404352 : zij = matrix_p(idim)%local_data(istate, jstate)
1512 : grad_ij = grad_ij + weights(idim)* &
1513 471744 : REAL(4.0_dp*CONJG(zij)*(zjj - zii), dp)
1514 : END DO
1515 69264 : matrix%local_data(istate, jstate) = grad_ij
1516 : END DO
1517 : END DO
1518 104 : DEALLOCATE (diag)
1519 104 : 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 4986 : 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 4986 : COMPLEX(KIND=dp), DIMENSION(:, :), POINTER :: diag
1536 : INTEGER :: idim, istate, jstate, ncol_local, &
1537 : nrow_global, nrow_local
1538 4986 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1539 : REAL(KIND=dp) :: grad_ij, ra, rb
1540 :
1541 4986 : 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 4986 : row_indices=row_indices, col_indices=col_indices)
1545 19944 : ALLOCATE (diag(nrow_global, SIZE(zij, 2)))
1546 4986 : value = 0.0_dp
1547 20058 : DO idim = 1, SIZE(zij, 2)
1548 203172 : DO istate = 1, nrow_global
1549 183114 : CALL cp_fm_get_element(zij(1, idim), istate, istate, ra)
1550 183114 : CALL cp_fm_get_element(zij(2, idim), istate, istate, rb)
1551 183114 : diag(istate, idim) = CMPLX(ra, rb, dp)
1552 198186 : value = value + weights(idim) - weights(idim)*ABS(diag(istate, idim))**2
1553 : END DO
1554 : END DO
1555 4986 : tolerance = 0.0_dp
1556 35410 : DO istate = 1, nrow_local
1557 504826 : DO jstate = 1, ncol_local
1558 1879089 : grad_ij = 0.0_dp
1559 1879089 : DO idim = 1, SIZE(zij, 2)
1560 1409673 : kii = diag(row_indices(istate), idim)
1561 1409673 : kjj = diag(col_indices(jstate), idim)
1562 1409673 : ra = zij(1, idim)%local_data(istate, jstate)
1563 1409673 : rb = zij(2, idim)%local_data(istate, jstate)
1564 1409673 : kij = CMPLX(ra, rb, dp)
1565 : grad_ij = grad_ij + weights(idim)* &
1566 1879089 : REAL(4.0_dp*CONJG(kij)*(kjj - kii), dp)
1567 : END DO
1568 499840 : tolerance = MAX(ABS(grad_ij), tolerance)
1569 : END DO
1570 : END DO
1571 4986 : CALL zij(1, 1)%matrix_struct%para_env%max(tolerance)
1572 :
1573 4986 : DEALLOCATE (diag)
1574 :
1575 4986 : 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 136 : 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 136 : 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 136 : 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 136 : CALL timeset(routineN, handle)
1623 136 : 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 136 : nrow_local=nrow_local, ncol_local=ncol_local)
1628 :
1629 136 : limit_crazy_angle = max_crazy_angle
1630 :
1631 136 : NULLIFY (diag_z, evals, evals_exp, mii, mij, mjj)
1632 136 : dim2 = SIZE(zij, 2)
1633 544 : ALLOCATE (diag_z(nrow_global, dim2))
1634 408 : ALLOCATE (evals(nrow_global))
1635 408 : ALLOCATE (evals_exp(nrow_global))
1636 :
1637 136 : CALL cp_cfm_create(cmat_A, zij(1, 1)%matrix_struct)
1638 136 : CALL cp_cfm_create(cmat_R, zij(1, 1)%matrix_struct)
1639 136 : CALL cp_cfm_create(cmat_t1, zij(1, 1)%matrix_struct)
1640 :
1641 136 : CALL cp_fm_create(mat_U, zij(1, 1)%matrix_struct)
1642 136 : CALL cp_fm_create(mat_t, zij(1, 1)%matrix_struct)
1643 136 : CALL cp_fm_create(mat_R, zij(1, 1)%matrix_struct)
1644 :
1645 136 : CALL cp_fm_create(mat_theta, zij(1, 1)%matrix_struct)
1646 :
1647 136 : CALL cp_fm_set_all(mat_R, 0.0_dp, 1.0_dp)
1648 136 : CALL cp_fm_set_all(mat_t, 0.0_dp)
1649 680 : ALLOCATE (mii(dim2), mij(dim2), mjj(dim2))
1650 550 : DO idim = 1, dim2
1651 414 : CALL cp_fm_scale_and_add(1.0_dp, mat_t, weights(idim), zij(1, idim))
1652 550 : CALL cp_fm_scale_and_add(1.0_dp, mat_t, weights(idim), zij(2, idim))
1653 : END DO
1654 136 : CALL cp_fm_syevd(mat_t, mat_U, evals)
1655 550 : DO idim = 1, dim2
1656 : ! rotate z's
1657 414 : 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 414 : 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 414 : 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 550 : 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 136 : CALL parallel_gemm('N', 'N', nrow_global, nrow_global, nrow_global, 1.0_dp, mat_R, mat_U, 0.0_dp, mat_t)
1664 136 : CALL cp_fm_to_fm(mat_t, mat_R)
1665 :
1666 136 : unit_nr = -1
1667 136 : IF (cmat_A%matrix_struct%para_env%is_source()) THEN
1668 68 : unit_nr = cp_logger_get_default_unit_nr()
1669 : WRITE (unit_nr, '(T2,A7,A6,1X,A20,A12,A12,A12)') &
1670 68 : "CRAZY| ", "Iter", "value ", "gradient", "Max. eval", "limit"
1671 : END IF
1672 :
1673 136 : iterations = 0
1674 136 : tolerance = 1.0_dp
1675 :
1676 : DO
1677 4986 : iterations = iterations + 1
1678 20058 : DO idim = 1, dim2
1679 203172 : DO i = 1, nrow_global
1680 183114 : CALL cp_fm_get_element(zij(1, idim), i, i, ra)
1681 183114 : CALL cp_fm_get_element(zij(2, idim), i, i, rb)
1682 198186 : diag_z(i, idim) = CMPLX(ra, rb, dp)
1683 : END DO
1684 : END DO
1685 35410 : DO irow = 1, nrow_local
1686 504826 : DO icol = 1, ncol_local
1687 1879089 : DO idim = 1, dim2
1688 1409673 : ra = zij(1, idim)%local_data(irow, icol)
1689 1409673 : rb = zij(2, idim)%local_data(irow, icol)
1690 1409673 : mij(idim) = CMPLX(ra, rb, dp)
1691 1409673 : mii(idim) = diag_z(row_indices(irow), idim)
1692 1879089 : mjj(idim) = diag_z(col_indices(icol), idim)
1693 : END DO
1694 499840 : IF (row_indices(irow) /= col_indices(icol)) THEN
1695 438992 : CALL get_angle(mii, mjj, mij, weights, theta)
1696 438992 : theta = crazy_scale*theta
1697 438992 : IF (theta > limit_crazy_angle) theta = limit_crazy_angle
1698 438992 : IF (theta < -limit_crazy_angle) theta = -limit_crazy_angle
1699 438992 : IF (crazy_use_diag) THEN
1700 0 : cmat_A%local_data(irow, icol) = -CMPLX(0.0_dp, theta, dp)
1701 : ELSE
1702 438992 : mat_theta%local_data(irow, icol) = -theta
1703 : END IF
1704 : ELSE
1705 30424 : IF (crazy_use_diag) THEN
1706 0 : cmat_A%local_data(irow, icol) = czero
1707 : ELSE
1708 30424 : 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 4986 : IF (crazy_use_diag) THEN
1717 0 : CALL cp_cfm_heevd(cmat_A, cmat_R, evals)
1718 0 : maxeval = MAXVAL(ABS(evals))
1719 0 : evals_exp(:) = EXP((0.0_dp, -1.0_dp)*evals(:))
1720 0 : CALL cp_cfm_to_cfm(cmat_R, cmat_t1)
1721 0 : CALL cp_cfm_column_scale(cmat_t1, evals_exp)
1722 : CALL parallel_gemm('N', 'C', nrow_global, nrow_global, nrow_global, cone, &
1723 0 : cmat_t1, cmat_R, czero, cmat_A)
1724 0 : mat_U%local_data = REAL(cmat_A%local_data, KIND=dp) ! U is a real matrix
1725 : ELSE
1726 4986 : do_emd = .FALSE.
1727 4986 : method = 2
1728 4986 : eps_exp = 1.0_dp*EPSILON(eps_exp)
1729 4986 : CALL cp_fm_maxabsrownorm(mat_theta, norm)
1730 4986 : maxeval = norm ! an upper bound
1731 4986 : CALL get_nsquare_norder(norm, nsquare, norder, eps_exp, method, do_emd)
1732 4986 : CALL exp_pade_real(mat_U, mat_theta, nsquare, norder)
1733 : END IF
1734 :
1735 20058 : DO idim = 1, dim2
1736 : ! rotate z's
1737 15072 : 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 15072 : 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 15072 : 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 20058 : 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 4986 : CALL parallel_gemm('N', 'N', nrow_global, nrow_global, nrow_global, 1.0_dp, mat_R, mat_U, 0.0_dp, mat_t)
1744 4986 : CALL cp_fm_to_fm(mat_t, mat_R)
1745 :
1746 4986 : CALL check_tolerance_new(weights, zij, tolerance, value)
1747 :
1748 4986 : IF (unit_nr > 0) THEN
1749 : WRITE (unit_nr, '(T2,A7,I6,1X,G20.15,E12.4,E12.4,E12.4)') &
1750 2493 : "CRAZY| ", iterations, value, tolerance, maxeval, limit_crazy_angle
1751 2493 : CALL m_flush(unit_nr)
1752 : END IF
1753 4986 : IF (tolerance < eps_localization .OR. iterations >= max_iter) EXIT
1754 : END DO
1755 :
1756 136 : IF (PRESENT(converged)) converged = (tolerance < eps_localization)
1757 :
1758 136 : CALL cp_cfm_release(cmat_A)
1759 136 : CALL cp_cfm_release(cmat_R)
1760 136 : CALL cp_cfm_release(cmat_T1)
1761 :
1762 136 : CALL cp_fm_release(mat_U)
1763 136 : CALL cp_fm_release(mat_T)
1764 136 : CALL cp_fm_release(mat_theta)
1765 :
1766 136 : CALL rotate_orbitals(mat_R, vectors)
1767 :
1768 136 : CALL cp_fm_release(mat_R)
1769 136 : DEALLOCATE (evals_exp, evals, diag_z)
1770 136 : DEALLOCATE (mii, mij, mjj)
1771 :
1772 136 : CALL timestop(handle)
1773 :
1774 272 : 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 2 : 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 2 : COMPLEX(KIND=dp), DIMENSION(:), POINTER :: evals_exp
1804 2 : 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 2 : 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 2 : 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 2 : 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 2 : NULLIFY (evals, evals_exp, diag_z, fval, fvald)
1825 :
1826 2 : CALL timeset(routineN, handle)
1827 2 : output_unit = cp_logger_get_default_io_unit()
1828 :
1829 2 : n = zij(1, 1)%matrix_struct%nrow_global
1830 2 : ndim = (SIZE(zij, 2))
1831 :
1832 2 : IF (output_unit > 0) THEN
1833 1 : WRITE (output_unit, '(T4,A )') "Localization by direct minimization of the functional; "
1834 1 : WRITE (output_unit, '(T5,2A13,A20,A20,A10 )') " Line search ", " Iteration ", " Functional ", " Tolerance ", " ds Min "
1835 : END IF
1836 :
1837 20 : ALLOCATE (evals(n), evals_exp(n), diag_z(n, ndim), fval(n), fvald(n))
1838 18 : ALLOCATE (c_zij(ndim))
1839 :
1840 : ! create the three complex matrices Z
1841 14 : DO idim = 1, ndim
1842 12 : 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 62 : zij(2, idim)%local_data, dp)
1845 : END DO
1846 :
1847 2 : CALL cp_fm_create(matrix_A, zij(1, 1)%matrix_struct)
1848 2 : CALL cp_fm_create(matrix_G, zij(1, 1)%matrix_struct)
1849 2 : CALL cp_fm_create(matrix_T, zij(1, 1)%matrix_struct)
1850 2 : CALL cp_fm_create(matrix_H, zij(1, 1)%matrix_struct)
1851 2 : CALL cp_fm_create(matrix_G_search, zij(1, 1)%matrix_struct)
1852 2 : CALL cp_fm_create(matrix_G_old, zij(1, 1)%matrix_struct)
1853 2 : CALL cp_fm_create(matrix_R, zij(1, 1)%matrix_struct)
1854 2 : CALL cp_fm_set_all(matrix_R, 0.0_dp, 1.0_dp)
1855 :
1856 2 : CALL cp_fm_set_all(matrix_A, 0.0_dp)
1857 : ! CALL cp_fm_init_random ( matrix_A )
1858 :
1859 2 : CALL cp_cfm_create(cmat_A, zij(1, 1)%matrix_struct)
1860 2 : CALL cp_cfm_create(cmat_U, zij(1, 1)%matrix_struct)
1861 2 : CALL cp_cfm_create(cmat_R, zij(1, 1)%matrix_struct)
1862 2 : CALL cp_cfm_create(cmat_t1, zij(1, 1)%matrix_struct)
1863 2 : CALL cp_cfm_create(cmat_t2, zij(1, 1)%matrix_struct)
1864 2 : CALL cp_cfm_create(cmat_B, zij(1, 1)%matrix_struct)
1865 2 : 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 2 : row_indices=row_indices, col_indices=col_indices)
1869 :
1870 2 : CALL cp_fm_set_all(matrix_G_old, 0.0_dp)
1871 2 : CALL cp_fm_set_all(matrix_G_search, 0.0_dp)
1872 2 : normg_old = 1.0E30_dp
1873 2 : ds_min = 1.0_dp
1874 2 : new_direction = .TRUE.
1875 2 : Iterations = 0
1876 2 : line_searches = 0
1877 2 : line_search_count = 0
1878 154 : DO
1879 154 : iterations = iterations + 1
1880 : ! compute U,R,evals given A
1881 770 : cmat_A%local_data = CMPLX(0.0_dp, matrix_A%local_data, dp) ! cmat_A is hermitian, evals are reals
1882 154 : CALL cp_cfm_heevd(cmat_A, cmat_R, evals)
1883 770 : evals_exp(:) = EXP((0.0_dp, -1.0_dp)*evals(:))
1884 154 : CALL cp_cfm_to_cfm(cmat_R, cmat_t1)
1885 154 : CALL cp_cfm_column_scale(cmat_t1, evals_exp)
1886 154 : CALL parallel_gemm('N', 'C', n, n, n, cone, cmat_t1, cmat_R, czero, cmat_U)
1887 770 : cmat_U%local_data = REAL(cmat_U%local_data, KIND=dp) ! enforce numerics, U is a real matrix
1888 :
1889 154 : IF (new_direction .AND. MOD(line_searches, 20) == 5) THEN ! reset with A .eq. 0
1890 0 : DO idim = 1, ndim
1891 0 : CALL parallel_gemm('N', 'N', n, n, n, cone, c_zij(idim), cmat_U, czero, cmat_t1)
1892 0 : 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 0 : matrix_H%local_data = REAL(cmat_U%local_data, KIND=dp)
1896 0 : CALL parallel_gemm('N', 'N', n, n, n, 1.0_dp, matrix_R, matrix_H, 0.0_dp, matrix_T)
1897 0 : CALL cp_fm_to_fm(matrix_T, matrix_R)
1898 :
1899 0 : CALL cp_cfm_set_all(cmat_U, czero, cone)
1900 0 : CALL cp_cfm_set_all(cmat_R, czero, cone)
1901 0 : CALL cp_cfm_set_all(cmat_A, czero)
1902 0 : CALL cp_fm_set_all(matrix_A, 0.0_dp)
1903 0 : evals(:) = 0.0_dp
1904 0 : evals_exp(:) = EXP((0.0_dp, -1.0_dp)*evals(:))
1905 0 : CALL cp_fm_set_all(matrix_G_old, 0.0_dp)
1906 0 : CALL cp_fm_set_all(matrix_G_search, 0.0_dp)
1907 0 : normg_old = 1.0E30_dp
1908 : END IF
1909 :
1910 : ! compute Omega and M
1911 154 : CALL cp_cfm_set_all(cmat_M, czero)
1912 154 : omega = 0.0_dp
1913 1078 : DO idim = 1, ndim
1914 924 : CALL parallel_gemm('N', 'N', n, n, n, cone, c_zij(idim), cmat_U, czero, cmat_t1) ! t1=ZU
1915 924 : CALL parallel_gemm('C', 'N', n, n, n, cone, cmat_U, cmat_t1, czero, cmat_t2) ! t2=(U^T)ZU
1916 2772 : DO i = 1, n
1917 1848 : 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 1848 : fval(i) = weights(idim) - weights(idim)*ABS(diag_z(i, idim))**2
1924 1848 : fvald(i) = -weights(idim)
1925 : END SELECT
1926 2772 : omega = omega + fval(i)
1927 : END DO
1928 2926 : DO icol = 1, ncol_local
1929 4620 : DO irow = 1, nrow_local
1930 1848 : 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 3696 : + 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 154 : 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 462 : DO icol = 1, ncol_local
1946 770 : DO irow = 1, nrow_local
1947 308 : ll = (0.0_dp, -1.0_dp)*evals(row_indices(irow))
1948 308 : lk = (0.0_dp, -1.0_dp)*evals(col_indices(icol))
1949 616 : IF (ABS(ll - lk) < 0.5_dp) THEN ! use a series expansion to avoid loss of precision
1950 176 : tmp = 1.0_dp
1951 176 : cmat_B%local_data(irow, icol) = 0.0_dp
1952 2992 : DO i = 1, 16
1953 2816 : cmat_B%local_data(irow, icol) = cmat_B%local_data(irow, icol) + tmp
1954 2992 : tmp = tmp*(ll - lk)/(i + 1)
1955 : END DO
1956 176 : cmat_B%local_data(irow, icol) = cmat_B%local_data(irow, icol)*EXP(lk)
1957 : ELSE
1958 132 : 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 154 : CALL parallel_gemm('C', 'N', n, n, n, cone, cmat_M, cmat_R, czero, cmat_t1) ! t1=(M^T)(R^T)
1965 154 : CALL parallel_gemm('C', 'N', n, n, n, cone, cmat_R, cmat_t1, czero, cmat_t2) ! t2=(R)t1
1966 154 : CALL cp_cfm_schur_product(cmat_t2, cmat_B, cmat_t1)
1967 154 : CALL parallel_gemm('N', 'C', n, n, n, cone, cmat_t1, cmat_R, czero, cmat_t2)
1968 154 : CALL parallel_gemm('N', 'N', n, n, n, cone, cmat_R, cmat_t2, czero, cmat_t1)
1969 770 : matrix_G%local_data = REAL(cmat_t1%local_data, KIND=dp)
1970 154 : CALL cp_fm_transpose(matrix_G, matrix_T)
1971 154 : CALL cp_fm_scale_and_add(-1.0_dp, matrix_G, 1.0_dp, matrix_T)
1972 154 : CALL cp_fm_maxabsval(matrix_G, tol)
1973 :
1974 : ! from here on, minimizing technology
1975 154 : IF (new_direction) THEN
1976 : ! energy converged up to machine precision ?
1977 6 : line_searches = line_searches + 1
1978 6 : IF (output_unit > 0) THEN
1979 3 : WRITE (output_unit, '(T5,I10,T18,I10,T31,2F20.6,F10.3)') line_searches, Iterations, Omega, tol, ds_min
1980 3 : CALL m_flush(output_unit)
1981 : END IF
1982 6 : IF (tol < eps_localization .OR. iterations > max_iter) EXIT
1983 :
1984 : IF (.TRUE.) THEN ! do conjugate gradient CG
1985 4 : CALL cp_fm_trace(matrix_G, matrix_G_old, normg_cross)
1986 4 : normg_cross = normg_cross*0.5_dp ! takes into account the fact that A is antisymmetric
1987 : ! apply the preconditioner
1988 12 : DO icol = 1, ncol_local
1989 20 : DO irow = 1, nrow_local
1990 16 : 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 4 : CALL cp_fm_trace(matrix_G, matrix_G_old, normg)
1994 4 : normg = normg*0.5_dp
1995 4 : beta_pr = (normg - normg_cross)/normg_old
1996 4 : normg_old = normg
1997 4 : beta_pr = MAX(beta_pr, 0.0_dp)
1998 4 : CALL cp_fm_scale_and_add(beta_pr, matrix_G_search, -1.0_dp, matrix_G_old)
1999 4 : CALL cp_fm_trace(matrix_G_search, matrix_G_old, normg_cross)
2000 4 : IF (normg_cross >= 0) THEN ! back to SD
2001 0 : IF (matrix_A%matrix_struct%para_env%is_source()) THEN
2002 0 : WRITE (cp_logger_get_default_unit_nr(), *) "!"
2003 : END IF
2004 0 : beta_pr = 0.0_dp
2005 0 : 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 152 : line_search_count = line_search_count + 1
2014 152 : 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 152 : new_direction = .FALSE.
2085 152 : IF (line_search_count == 1) THEN
2086 4 : lsl = 1
2087 4 : lsr = 0
2088 4 : lsm = 1
2089 4 : pos(1) = 0.0_dp
2090 4 : pos(2) = ds_min/gold_sec
2091 : ELSE
2092 148 : IF (line_search_count == 150) CPABORT("Too many")
2093 148 : IF (lsr == 0) THEN
2094 28 : IF (energy(line_search_count - 1) < energy(line_search_count)) THEN
2095 4 : lsr = line_search_count
2096 4 : pos(line_search_count + 1) = pos(lsm) + (pos(lsr) - pos(lsm))*gold_sec
2097 : ELSE
2098 24 : lsl = lsm
2099 24 : lsm = line_search_count
2100 24 : pos(line_search_count + 1) = pos(line_search_count)/gold_sec
2101 : END IF
2102 : ELSE
2103 120 : IF (pos(line_search_count) < pos(lsm)) THEN
2104 38 : IF (energy(line_search_count) < energy(lsm)) THEN
2105 : lsr = lsm
2106 : lsm = line_search_count
2107 : ELSE
2108 28 : lsl = line_search_count
2109 : END IF
2110 : ELSE
2111 82 : IF (energy(line_search_count) < energy(lsm)) THEN
2112 : lsl = lsm
2113 : lsm = line_search_count
2114 : ELSE
2115 64 : lsr = line_search_count
2116 : END IF
2117 : END IF
2118 120 : IF (pos(lsr) - pos(lsm) > pos(lsm) - pos(lsl)) THEN
2119 78 : pos(line_search_count + 1) = pos(lsm) + gold_sec*(pos(lsr) - pos(lsm))
2120 : ELSE
2121 42 : pos(line_search_count + 1) = pos(lsl) + gold_sec*(pos(lsm) - pos(lsl))
2122 : END IF
2123 120 : IF ((pos(lsr) - pos(lsl)) < 1.0E-3_dp*pos(lsr)) THEN
2124 4 : 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 152 : ds_min = pos(line_search_count + 1)
2131 152 : ds = pos(line_search_count + 1) - pos(line_search_count)
2132 154 : CALL cp_fm_scale_and_add(1.0_dp, matrix_A, ds, matrix_G_search)
2133 : END DO
2134 :
2135 : ! collect rotation matrix
2136 10 : matrix_H%local_data = REAL(cmat_U%local_data, KIND=dp)
2137 2 : CALL parallel_gemm('N', 'N', n, n, n, 1.0_dp, matrix_R, matrix_H, 0.0_dp, matrix_T)
2138 2 : CALL cp_fm_to_fm(matrix_T, matrix_R)
2139 2 : CALL rotate_orbitals(matrix_R, vectors)
2140 2 : CALL cp_fm_release(matrix_R)
2141 :
2142 2 : CALL cp_fm_release(matrix_A)
2143 2 : CALL cp_fm_release(matrix_G)
2144 2 : CALL cp_fm_release(matrix_H)
2145 2 : CALL cp_fm_release(matrix_T)
2146 2 : CALL cp_fm_release(matrix_G_search)
2147 2 : CALL cp_fm_release(matrix_G_old)
2148 2 : CALL cp_cfm_release(cmat_A)
2149 2 : CALL cp_cfm_release(cmat_U)
2150 2 : CALL cp_cfm_release(cmat_R)
2151 2 : CALL cp_cfm_release(cmat_t1)
2152 2 : CALL cp_cfm_release(cmat_t2)
2153 2 : CALL cp_cfm_release(cmat_B)
2154 2 : CALL cp_cfm_release(cmat_M)
2155 :
2156 2 : DEALLOCATE (evals, evals_exp, fval, fvald)
2157 :
2158 14 : DO idim = 1, SIZE(c_zij)
2159 60 : zij(1, idim)%local_data = REAL(c_zij(idim)%local_data, dp)
2160 60 : zij(2, idim)%local_data = AIMAG(c_zij(idim)%local_data)
2161 14 : CALL cp_cfm_release(c_zij(idim))
2162 : END DO
2163 2 : DEALLOCATE (c_zij)
2164 2 : DEALLOCATE (diag_z)
2165 :
2166 2 : CALL timestop(handle)
2167 :
2168 6 : 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 390 : 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 390 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: ns_bound
2206 390 : 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 390 : CALL timeset(routineN, handle)
2212 :
2213 390 : output_unit = cp_logger_get_default_io_unit()
2214 :
2215 : NULLIFY (cz_ij_loc)
2216 :
2217 390 : dim2 = SIZE(zij, 2)
2218 :
2219 390 : CALL cp_fm_create(rmat, zij(1, 1)%matrix_struct)
2220 390 : CALL cp_fm_set_all(rmat, 0._dp, 1._dp)
2221 :
2222 390 : CALL cp_fm_get_info(rmat, nrow_global=nstate)
2223 :
2224 390 : IF (restricted > 0) THEN
2225 0 : IF (output_unit > 0) THEN
2226 0 : WRITE (output_unit, '(T4,A,I2,A )') "JACOBI: for the ROKS method, the last ", restricted, " orbitals DO NOT ROTATE"
2227 : END IF
2228 0 : nstate = nstate - restricted
2229 : END IF
2230 :
2231 : ! Distribution of the states (XXXXX safe against more pe than states ??? XXXXX)
2232 390 : xstate = REAL(nstate, dp)/REAL(para_env%num_pe, dp)
2233 1560 : ALLOCATE (ns_bound(0:para_env%num_pe - 1, 2))
2234 1170 : DO ip = 1, para_env%num_pe
2235 780 : ns_bound(ip - 1, 1) = MIN(nstate, NINT(xstate*(ip - 1))) + 1
2236 1170 : ns_bound(ip - 1, 2) = MIN(nstate, NINT(xstate*ip))
2237 : END DO
2238 390 : nblock_max = 0
2239 1170 : DO ip = 0, para_env%num_pe - 1
2240 780 : nblock = ns_bound(ip, 2) - ns_bound(ip, 1) + 1
2241 1170 : 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 1560 : ALLOCATE (z_ij_loc_re(nstate, nblock_max))
2246 1170 : ALLOCATE (z_ij_loc_im(nstate, nblock_max))
2247 2364 : ALLOCATE (cz_ij_loc(dim2))
2248 1584 : DO idim = 1, dim2
2249 3972 : DO ip = 0, para_env%num_pe - 1
2250 2388 : nblock = ns_bound(ip, 2) - ns_bound(ip, 1) + 1
2251 2388 : CALL cp_fm_get_submatrix(zij(1, idim), z_ij_loc_re, 1, ns_bound(ip, 1), nstate, nblock)
2252 2388 : CALL cp_fm_get_submatrix(zij(2, idim), z_ij_loc_im, 1, ns_bound(ip, 1), nstate, nblock)
2253 3582 : IF (para_env%mepos == ip) THEN
2254 4752 : ALLOCATE (cz_ij_loc(idim)%c_array(nstate, nblock))
2255 5034 : DO i = 1, nblock
2256 37866 : DO j = 1, nstate
2257 36672 : 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 390 : DEALLOCATE (z_ij_loc_re)
2264 390 : DEALLOCATE (z_ij_loc_im)
2265 :
2266 1560 : 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 390 : target_time=target_time, start_time=start_time)
2271 :
2272 390 : ilow1 = ns_bound(para_env%mepos, 1)
2273 390 : ns_me = ns_bound(para_env%mepos, 2) - ns_bound(para_env%mepos, 1) + 1
2274 1170 : ALLOCATE (z_ij_loc_re(nstate, nblock_max))
2275 1170 : ALLOCATE (z_ij_loc_im(nstate, nblock_max))
2276 1584 : DO idim = 1, dim2
2277 3972 : DO ip = 0, para_env%num_pe - 1
2278 2388 : z_ij_loc_re = 0.0_dp
2279 2388 : z_ij_loc_im = 0.0_dp
2280 2388 : nblock = ns_bound(ip, 2) - ns_bound(ip, 1) + 1
2281 2388 : IF (ip == para_env%mepos) THEN
2282 5034 : ns_me = nblock
2283 5034 : DO i = 1, ns_me
2284 36672 : ii = ilow1 + i - 1
2285 37866 : DO j = 1, nstate
2286 32832 : z_ij_loc_re(j, i) = REAL(cz_ij_loc(idim)%c_array(j, i), dp)
2287 36672 : 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 2388 : CALL para_env%bcast(z_ij_loc_re, ip)
2292 2388 : CALL para_env%bcast(z_ij_loc_im, ip)
2293 2388 : CALL cp_fm_set_submatrix(zij(1, idim), z_ij_loc_re, 1, ns_bound(ip, 1), nstate, nblock)
2294 3582 : 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 1170 : DO ip = 0, para_env%num_pe - 1
2299 780 : z_ij_loc_re = 0.0_dp
2300 780 : nblock = ns_bound(ip, 2) - ns_bound(ip, 1) + 1
2301 780 : IF (ip == para_env%mepos) THEN
2302 1633 : ns_me = nblock
2303 1633 : DO i = 1, ns_me
2304 11780 : ii = ilow1 + i - 1
2305 12170 : DO j = 1, nstate
2306 11780 : z_ij_loc_re(j, i) = rotmat(j, i)
2307 : END DO
2308 : END DO
2309 : END IF
2310 780 : CALL para_env%bcast(z_ij_loc_re, ip)
2311 1170 : CALL cp_fm_set_submatrix(rmat, z_ij_loc_re, 1, ns_bound(ip, 1), nstate, nblock)
2312 : END DO
2313 :
2314 390 : DEALLOCATE (z_ij_loc_re)
2315 390 : DEALLOCATE (z_ij_loc_im)
2316 1584 : DO idim = 1, dim2
2317 1584 : DEALLOCATE (cz_ij_loc(idim)%c_array)
2318 : END DO
2319 390 : DEALLOCATE (cz_ij_loc)
2320 :
2321 390 : CALL para_env%sync()
2322 390 : CALL rotate_orbitals(rmat, vectors)
2323 390 : CALL cp_fm_release(rmat)
2324 :
2325 390 : DEALLOCATE (rotmat)
2326 390 : DEALLOCATE (ns_bound)
2327 :
2328 390 : CALL timestop(handle)
2329 :
2330 1170 : 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 32 : 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 32 : 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 32 : 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 32 : CALL timeset(routineN, handle)
2363 :
2364 32 : dim2 = SIZE(czij)
2365 :
2366 32 : CALL cp_cfm_set_all(rmat, CMPLX(0._dp, 0._dp, dp), CMPLX(1._dp, 0._dp, dp))
2367 :
2368 32 : CALL cp_cfm_get_info(rmat, nrow_global=nstate)
2369 :
2370 : ! Distribution of the states (XXXXX safe against more pe than states ??? XXXXX)
2371 32 : xstate = REAL(nstate, dp)/REAL(para_env%num_pe, dp)
2372 128 : ALLOCATE (ns_bound(0:para_env%num_pe - 1, 2))
2373 96 : DO ip = 1, para_env%num_pe
2374 64 : ns_bound(ip - 1, 1) = MIN(nstate, NINT(xstate*(ip - 1))) + 1
2375 96 : ns_bound(ip - 1, 2) = MIN(nstate, NINT(xstate*ip))
2376 : END DO
2377 32 : nblock_max = 0
2378 96 : DO ip = 0, para_env%num_pe - 1
2379 64 : nblock = ns_bound(ip, 2) - ns_bound(ip, 1) + 1
2380 96 : 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 128 : ALLOCATE (czij_array(nstate, nblock_max))
2385 192 : ALLOCATE (cz_ij_loc(dim2))
2386 128 : DO idim = 1, dim2
2387 320 : DO ip = 0, para_env%num_pe - 1
2388 192 : nblock = ns_bound(ip, 2) - ns_bound(ip, 1) + 1
2389 : ! cfm --> allocatable
2390 192 : CALL cp_cfm_get_submatrix(czij(idim), czij_array, 1, ns_bound(ip, 1), nstate, nblock)
2391 288 : IF (para_env%mepos == ip) THEN
2392 96 : ns_me = nblock
2393 384 : ALLOCATE (cz_ij_loc(idim)%c_array(nstate, ns_me))
2394 2160 : DO i = 1, ns_me
2395 90912 : DO j = 1, nstate
2396 90816 : 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 32 : DEALLOCATE (czij_array)
2403 :
2404 128 : 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 32 : cz_ij_loc, rotmat, 0, tol_out=tol_out)
2408 :
2409 32 : ilow1 = ns_bound(para_env%mepos, 1)
2410 32 : ns_me = ns_bound(para_env%mepos, 2) - ns_bound(para_env%mepos, 1) + 1
2411 128 : ALLOCATE (z_ij_loc_re(nstate, nblock_max))
2412 :
2413 96 : DO ip = 0, para_env%num_pe - 1
2414 64 : z_ij_loc_re = 0.0_dp
2415 64 : nblock = ns_bound(ip, 2) - ns_bound(ip, 1) + 1
2416 64 : IF (ip == para_env%mepos) THEN
2417 720 : ns_me = nblock
2418 720 : DO i = 1, ns_me
2419 30272 : ii = ilow1 + i - 1
2420 30304 : DO j = 1, nstate
2421 30272 : z_ij_loc_re(j, i) = rotmat(j, i)
2422 : END DO
2423 : END DO
2424 : END IF
2425 64 : CALL para_env%bcast(z_ij_loc_re, ip)
2426 62048 : 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 32 : DEALLOCATE (z_ij_loc_re)
2430 128 : DO idim = 1, dim2
2431 128 : DEALLOCATE (cz_ij_loc(idim)%c_array)
2432 : END DO
2433 32 : DEALLOCATE (cz_ij_loc)
2434 :
2435 32 : CALL para_env%sync()
2436 :
2437 32 : DEALLOCATE (rotmat)
2438 32 : DEALLOCATE (ns_bound)
2439 :
2440 32 : CALL timestop(handle)
2441 :
2442 64 : 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 1266 : SUBROUTINE jacobi_rot_para_core(weights, para_env, max_iter, sweeps, out_each, dim2, nstate, nblock_max, &
2467 422 : 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 422 : 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 422 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: gmat, rmat_loc, rmat_recv, rmat_send
2493 422 : 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 422 : TYPE(set_c_1d_type), DIMENSION(:), POINTER :: zdiag_all, zdiag_me
2497 422 : TYPE(set_c_2d_type), DIMENSION(:), POINTER :: xyz_mix, xyz_mix_ns
2498 :
2499 422 : NULLIFY (zdiag_all, zdiag_me)
2500 422 : NULLIFY (xyz_mix, xyz_mix_ns)
2501 : NULLIFY (mii, mij, mjj)
2502 :
2503 2110 : ALLOCATE (mii(dim2), mij(dim2), mjj(dim2))
2504 :
2505 1266 : ALLOCATE (rcount(para_env%num_pe))
2506 844 : ALLOCATE (rdispl(para_env%num_pe))
2507 :
2508 422 : tolerance = 1.0e10_dp
2509 422 : sweeps = 0
2510 :
2511 : ! number of processor pairs and number of permutations
2512 422 : npair = (para_env%num_pe + 1)/2
2513 422 : nperm = para_env%num_pe - MOD(para_env%num_pe + 1, 2)
2514 1266 : ALLOCATE (list_pair(2, npair))
2515 :
2516 : ! initialize rotation matrix
2517 87366 : rotmat = 0.0_dp
2518 2353 : DO i = ns_bound(para_env%mepos, 1), ns_bound(para_env%mepos, 2)
2519 1931 : ii = i - ns_bound(para_env%mepos, 1) + 1
2520 2353 : rotmat(i, ii) = 1.0_dp
2521 : END DO
2522 :
2523 2556 : ALLOCATE (xyz_mix(dim2))
2524 2134 : ALLOCATE (xyz_mix_ns(dim2))
2525 2556 : ALLOCATE (zdiag_me(dim2))
2526 2134 : ALLOCATE (zdiag_all(dim2))
2527 :
2528 422 : ns_me = ns_bound(para_env%mepos, 2) - ns_bound(para_env%mepos, 1) + 1
2529 422 : IF (ns_me /= 0) THEN
2530 2070 : ALLOCATE (c_array_me(nstate, ns_me, dim2))
2531 1680 : DO idim = 1, dim2
2532 5478 : ALLOCATE (xyz_mix_ns(idim)%c_array(nstate, ns_me))
2533 : END DO
2534 1656 : ALLOCATE (gmat(nstate, ns_me))
2535 : END IF
2536 :
2537 1712 : DO idim = 1, dim2
2538 3870 : ALLOCATE (zdiag_me(idim)%c_array(nblock_max))
2539 7542 : zdiag_me(idim)%c_array = z_zero
2540 3870 : ALLOCATE (zdiag_all(idim)%c_array(para_env%num_pe*nblock_max))
2541 14216 : zdiag_all(idim)%c_array = z_zero
2542 : END DO
2543 1688 : ALLOCATE (rmat_recv(nblock_max*2, nblock_max))
2544 1266 : ALLOCATE (rmat_send(nblock_max*2, nblock_max))
2545 :
2546 : ! buffer for message passing
2547 2110 : ALLOCATE (rmat_recv_all(nblock_max*2, nblock_max, 0:para_env%num_pe - 1))
2548 :
2549 422 : IF (output_unit > 0) THEN
2550 195 : WRITE (output_unit, '(T4,A )') " Localization by iterative distributed Jacobi rotation"
2551 195 : WRITE (output_unit, '(T20,A12,T32, A22,T60, A12,A8 )') "Iteration", "Functional", "Tolerance", " Time "
2552 : END IF
2553 :
2554 87474 : DO lsweep = 1, max_iter + 1
2555 87474 : sweeps = lsweep
2556 87474 : IF (sweeps == max_iter + 1) THEN
2557 158 : IF (output_unit > 0) THEN
2558 63 : WRITE (output_unit, *) ' LOCALIZATION! loop did not converge within the maximum number of iterations.'
2559 63 : WRITE (output_unit, *) ' Present Max. gradient = ', tolerance
2560 : END IF
2561 : EXIT
2562 : END IF
2563 87316 : t1 = m_walltime()
2564 :
2565 174632 : DO iperm = 1, nperm
2566 :
2567 : ! fix partners for this permutation, and get the number of states
2568 87316 : CALL eberlein(iperm, para_env, list_pair)
2569 87316 : ip_partner = -1
2570 87316 : ns_partner = 0
2571 87316 : DO ipair = 1, npair
2572 87316 : IF (list_pair(1, ipair) == para_env%mepos) THEN
2573 43658 : ip_partner = list_pair(2, ipair)
2574 43658 : EXIT
2575 43658 : ELSE IF (list_pair(2, ipair) == para_env%mepos) THEN
2576 43658 : ip_partner = list_pair(1, ipair)
2577 43658 : EXIT
2578 : END IF
2579 : END DO
2580 87316 : IF (ip_partner >= 0) THEN
2581 87316 : 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 87316 : IF (ns_partner*ns_me /= 0) THEN
2588 :
2589 349200 : ALLOCATE (rmat_loc(ns_me + ns_partner, ns_me + ns_partner))
2590 87300 : rmat_loc = 0.0_dp
2591 701456 : DO i = 1, ns_me + ns_partner
2592 701456 : rmat_loc(i, i) = 1.0_dp
2593 : END DO
2594 :
2595 436500 : ALLOCATE (c_array_partner(nstate, ns_partner, dim2))
2596 :
2597 351300 : DO idim = 1, dim2
2598 1056000 : ALLOCATE (xyz_mix(idim)%c_array(ns_me + ns_partner, ns_me + ns_partner))
2599 1284915 : DO i = 1, ns_me
2600 7980156 : 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 87300 : msgout=c_array_partner, source=ip_partner)
2606 :
2607 87300 : n1 = ns_me
2608 87300 : n2 = ns_partner
2609 87300 : ilow1 = ns_bound(para_env%mepos, 1)
2610 87300 : iup1 = ns_bound(para_env%mepos, 1) + n1 - 1
2611 87300 : ilow2 = ns_bound(ip_partner, 1)
2612 87300 : iup2 = ns_bound(ip_partner, 1) + n2 - 1
2613 87300 : IF (ns_bound(para_env%mepos, 1) < ns_bound(ip_partner, 1)) THEN
2614 43650 : il1 = 1
2615 : iu1 = n1
2616 43650 : iu1 = n1
2617 43650 : il2 = 1 + n1
2618 43650 : iu2 = n1 + n2
2619 : ELSE
2620 43650 : il1 = 1 + n2
2621 : iu1 = n1 + n2
2622 43650 : iu1 = n1 + n2
2623 43650 : il2 = 1
2624 43650 : iu2 = n2
2625 : END IF
2626 :
2627 351300 : DO idim = 1, dim2
2628 1197615 : DO i = 1, n1
2629 4386354 : xyz_mix(idim)%c_array(il1:iu1, il1 + i - 1) = c_array_me(ilow1:iup1, i, idim)
2630 4527417 : xyz_mix(idim)%c_array(il2:iu2, il1 + i - 1) = c_array_me(ilow2:iup2, i, idim)
2631 : END DO
2632 1284915 : DO i = 1, n2
2633 4386354 : xyz_mix(idim)%c_array(il2:iu2, il2 + i - 1) = c_array_partner(ilow2:iup2, i, idim)
2634 4527417 : 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 701456 : DO istate = 1, n1 + n2
2639 2604488 : DO jstate = istate + 1, n1 + n2
2640 7751958 : DO idim = 1, dim2
2641 5848926 : mii(idim) = xyz_mix(idim)%c_array(istate, istate)
2642 5848926 : mij(idim) = xyz_mix(idim)%c_array(istate, jstate)
2643 7751958 : mjj(idim) = xyz_mix(idim)%c_array(jstate, jstate)
2644 : END DO
2645 1903032 : CALL get_angle(mii, mjj, mij, weights, theta)
2646 1903032 : st = SIN(theta)
2647 1903032 : ct = COS(theta)
2648 7751958 : DO idim = 1, dim2
2649 52347834 : DO i = 1, n1 + n2
2650 46498908 : zi = ct*xyz_mix(idim)%c_array(i, istate) + st*xyz_mix(idim)%c_array(i, jstate)
2651 46498908 : zj = -st*xyz_mix(idim)%c_array(i, istate) + ct*xyz_mix(idim)%c_array(i, jstate)
2652 46498908 : xyz_mix(idim)%c_array(i, istate) = zi
2653 52347834 : xyz_mix(idim)%c_array(i, jstate) = zj
2654 : END DO
2655 54250866 : DO i = 1, n1 + n2
2656 46498908 : zi = ct*xyz_mix(idim)%c_array(istate, i) + st*xyz_mix(idim)%c_array(jstate, i)
2657 46498908 : zj = -st*xyz_mix(idim)%c_array(istate, i) + ct*xyz_mix(idim)%c_array(jstate, i)
2658 46498908 : xyz_mix(idim)%c_array(istate, i) = zi
2659 52347834 : xyz_mix(idim)%c_array(jstate, i) = zj
2660 : END DO
2661 : END DO
2662 :
2663 17423024 : DO i = 1, n1 + n2
2664 14905836 : ri = ct*rmat_loc(i, istate) + st*rmat_loc(i, jstate)
2665 14905836 : rj = ct*rmat_loc(i, jstate) - st*rmat_loc(i, istate)
2666 14905836 : rmat_loc(i, istate) = ri
2667 16808868 : rmat_loc(i, jstate) = rj
2668 : END DO
2669 : END DO
2670 : END DO
2671 :
2672 87300 : k = nblock_max + 1
2673 : CALL para_env%sendrecv(rotmat(1:nstate, 1:ns_me), ip_partner, &
2674 5121676 : rotmat(1:nstate, k:k + n2 - 1), ip_partner)
2675 :
2676 87300 : 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 1476102 : 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 43650 : 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 43650 : 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 1151398 : n1, 1.0_dp, gmat(:, :), nstate)
2690 : END IF
2691 :
2692 87300 : CALL dcopy(nstate*n1, gmat(1, 1), 1, rotmat(1, 1), 1)
2693 :
2694 351300 : DO idim = 1, dim2
2695 1197615 : DO i = 1, n1
2696 7980156 : xyz_mix_ns(idim)%c_array(1:nstate, i) = z_zero
2697 : END DO
2698 :
2699 1197615 : DO istate = 1, n1
2700 7980156 : DO jstate = 1, nstate
2701 33924672 : DO i = 1, n2
2702 : xyz_mix_ns(idim)%c_array(jstate, istate) = &
2703 : xyz_mix_ns(idim)%c_array(jstate, istate) + &
2704 32991057 : c_array_partner(jstate, i, idim)*rmat_loc(il2 + i - 1, il1 + istate - 1)
2705 : END DO
2706 : END DO
2707 : END DO
2708 1284915 : DO istate = 1, n1
2709 7980156 : DO jstate = 1, nstate
2710 34789089 : DO i = 1, n1
2711 : xyz_mix_ns(idim)%c_array(jstate, istate) = xyz_mix_ns(idim)%c_array(jstate, istate) + &
2712 33855474 : 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 87300 : DEALLOCATE (c_array_partner)
2719 :
2720 : ELSE ! save my data
2721 64 : DO idim = 1, dim2
2722 88 : DO i = 1, ns_me
2723 120 : 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 351364 : DO idim = 1, dim2
2729 1285003 : DO i = 1, ns_me
2730 7980252 : cz_ij_loc(idim)%c_array(1:nstate, i) = z_zero
2731 : END DO
2732 : END DO
2733 :
2734 87316 : IF (ns_partner*ns_me /= 0) THEN
2735 : ! transpose rotation matrix rmat_loc
2736 701456 : DO i = 1, ns_me + ns_partner
2737 2604488 : DO j = i + 1, ns_me + ns_partner
2738 1903032 : ri = rmat_loc(i, j)
2739 1903032 : rmat_loc(i, j) = rmat_loc(j, i)
2740 2517188 : rmat_loc(j, i) = ri
2741 : END DO
2742 : END DO
2743 :
2744 : ! prepare for distribution
2745 394378 : DO i = 1, n1
2746 1519752 : rmat_send(1:n1, i) = rmat_loc(il1:iu1, il1 + i - 1)
2747 : END DO
2748 87300 : ik = nblock_max
2749 394378 : DO i = 1, n2
2750 1479114 : rmat_send(ik + 1:ik + n1, i) = rmat_loc(il1:iu1, il2 + i - 1)
2751 : END DO
2752 : ELSE
2753 16 : rmat_send = 0.0_dp
2754 : END IF
2755 :
2756 : ! collect data from all tasks (this takes some significant time)
2757 87316 : CALL para_env%allgather(rmat_send, rmat_recv_all)
2758 :
2759 : ! update blocks everywhere
2760 261948 : DO ip = 0, para_env%num_pe - 1
2761 :
2762 174632 : ip_recv_from = MOD(para_env%mepos - IP + para_env%num_pe, para_env%num_pe)
2763 6507920 : rmat_recv(:, :) = rmat_recv_all(:, :, ip_recv_from)
2764 :
2765 174632 : ns_recv_from = ns_bound(ip_recv_from, 2) - ns_bound(ip_recv_from, 1) + 1
2766 :
2767 261948 : IF (ns_me /= 0) THEN
2768 174616 : IF (ns_recv_from /= 0) THEN
2769 : !look for the partner of ip_recv_from
2770 174608 : ip_recv_partner = -1
2771 174608 : ns_recv_partner = 0
2772 174608 : DO ipair = 1, npair
2773 174608 : IF (list_pair(1, ipair) == ip_recv_from) THEN
2774 87308 : ip_recv_partner = list_pair(2, ipair)
2775 87308 : EXIT
2776 87300 : 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 174608 : IF (ip_recv_partner >= 0) THEN
2783 174608 : ns_recv_partner = ns_bound(ip_recv_partner, 2) - ns_bound(ip_recv_partner, 1) + 1
2784 : END IF
2785 174608 : IF (ns_recv_partner > 0) THEN
2786 174600 : il1 = ns_bound(para_env%mepos, 1)
2787 174600 : il_recv = ns_bound(ip_recv_from, 1)
2788 174600 : il_recv_partner = ns_bound(ip_recv_partner, 1)
2789 174600 : ik = nblock_max
2790 :
2791 702600 : DO idim = 1, dim2
2792 2395230 : DO i = 1, ns_recv_from
2793 1867230 : ii = il_recv + i - 1
2794 9177771 : DO j = 1, ns_me
2795 33855474 : jj = j
2796 35722704 : DO k = 1, ns_recv_from
2797 27072933 : kk = il_recv + k - 1
2798 : cz_ij_loc(idim)%c_array(ii, jj) = cz_ij_loc(idim)%c_array(ii, jj) + &
2799 33855474 : rmat_recv(i, k)*xyz_mix_ns(idim)%c_array(kk, j)
2800 : END DO
2801 : END DO
2802 : END DO
2803 2569830 : DO i = 1, ns_recv_from
2804 1867230 : ii = il_recv + i - 1
2805 9177771 : DO j = 1, ns_me
2806 32991057 : jj = j
2807 34858287 : DO k = 1, ns_recv_partner
2808 26208516 : kk = il_recv_partner + k - 1
2809 : cz_ij_loc(idim)%c_array(ii, jj) = cz_ij_loc(idim)%c_array(ii, jj) + &
2810 32991057 : 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 8 : il1 = ns_bound(para_env%mepos, 1)
2817 8 : il_recv = ns_bound(ip_recv_from, 1)
2818 32 : DO idim = 1, dim2
2819 56 : DO j = 1, ns_me
2820 48 : jj = j
2821 72 : DO i = 1, ns_recv_from
2822 24 : ii = il_recv + i - 1
2823 48 : 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 174632 : IF (ns_partner*ns_me /= 0) THEN
2833 87300 : DEALLOCATE (rmat_loc)
2834 351300 : DO idim = 1, dim2
2835 351300 : DEALLOCATE (xyz_mix(idim)%c_array)
2836 : END DO
2837 : END IF
2838 :
2839 : END DO ! iperm
2840 :
2841 : ! calculate the max gradient
2842 351364 : DO idim = 1, dim2
2843 1197687 : DO i = ns_bound(para_env%mepos, 1), ns_bound(para_env%mepos, 2)
2844 933639 : ii = i - ns_bound(para_env%mepos, 1) + 1
2845 933639 : zdiag_me(idim)%c_array(ii) = cz_ij_loc(idim)%c_array(i, ii)
2846 1197687 : zdiag_me(idim)%c_array(ii) = cz_ij_loc(idim)%c_array(i, ii)
2847 : END DO
2848 792144 : rcount(:) = SIZE(zdiag_me(idim)%c_array)
2849 264048 : rdispl(1) = 0
2850 528096 : DO ip = 2, para_env%num_pe
2851 528096 : rdispl(ip) = rdispl(ip - 1) + rcount(ip - 1)
2852 : END DO
2853 : ! collect all the diagonal elements in a replicated 1d array
2854 3521164 : CALL para_env%allgatherv(zdiag_me(idim)%c_array, zdiag_all(idim)%c_array, rcount, rdispl)
2855 : END DO
2856 :
2857 87316 : gmax = 0.0_dp
2858 394402 : DO j = ns_bound(para_env%mepos, 1), ns_bound(para_env%mepos, 2)
2859 307086 : k = j - ns_bound(para_env%mepos, 1) + 1
2860 1345918 : DO i = 1, j - 1
2861 : ! find the location of the diagonal element (i,i)
2862 1095233 : DO ip = 0, para_env%num_pe - 1
2863 1095233 : 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 951516 : 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 951516 : jj = nblock_max*para_env%mepos + j - ns_bound(para_env%mepos, 1) + 1
2871 951516 : grad = 0.0_dp
2872 3875979 : DO idim = 1, dim2
2873 2924463 : zi = zdiag_all(idim)%c_array(ii)
2874 2924463 : zj = zdiag_all(idim)%c_array(jj)
2875 3875979 : grad = grad + weights(idim)*REAL(4.0_dp*CONJG(cz_ij_loc(idim)%c_array(i, k))*(zj - zi), dp)
2876 : END DO
2877 1258602 : gmax = MAX(gmax, ABS(grad))
2878 : END DO
2879 : END DO
2880 :
2881 87316 : CALL para_env%max(gmax)
2882 87316 : tolerance = gmax
2883 87316 : IF (PRESENT(tol_out)) tol_out = tolerance
2884 :
2885 87316 : func = 0.0_dp
2886 394402 : DO i = ns_bound(para_env%mepos, 1), ns_bound(para_env%mepos, 2)
2887 307086 : k = i - ns_bound(para_env%mepos, 1) + 1
2888 1328041 : DO idim = 1, dim2
2889 933639 : zr = REAL(cz_ij_loc(idim)%c_array(i, k), dp)
2890 933639 : zc = AIMAG(cz_ij_loc(idim)%c_array(i, k))
2891 1240725 : func = func + weights(idim)*(1.0_dp - (zr*zr + zc*zc))/twopi/twopi
2892 : END DO
2893 : END DO
2894 87316 : CALL para_env%sum(func)
2895 87316 : t2 = m_walltime()
2896 :
2897 87316 : IF (output_unit > 0 .AND. MODULO(sweeps, out_each) == 0) THEN
2898 444 : WRITE (output_unit, '(T20,I12,T35,F20.10,T60,E12.4,F8.3)') sweeps, func, tolerance, t2 - t1
2899 444 : CALL m_flush(output_unit)
2900 : END IF
2901 87316 : IF (PRESENT(eps_localization)) THEN
2902 87284 : IF (tolerance < eps_localization) EXIT
2903 : END IF
2904 87316 : IF (PRESENT(target_time) .AND. PRESENT(start_time)) THEN
2905 87020 : CALL external_control(should_stop, "LOC", target_time=target_time, start_time=start_time)
2906 87020 : IF (should_stop) EXIT
2907 : END IF
2908 :
2909 : END DO ! lsweep
2910 :
2911 : ! buffer for message passing
2912 422 : DEALLOCATE (rmat_recv_all)
2913 :
2914 422 : DEALLOCATE (rmat_recv)
2915 422 : DEALLOCATE (rmat_send)
2916 422 : IF (ns_me > 0) THEN
2917 414 : DEALLOCATE (c_array_me)
2918 : END IF
2919 1712 : DO idim = 1, dim2
2920 1290 : DEALLOCATE (zdiag_me(idim)%c_array)
2921 1712 : DEALLOCATE (zdiag_all(idim)%c_array)
2922 : END DO
2923 422 : DEALLOCATE (zdiag_me)
2924 422 : DEALLOCATE (zdiag_all)
2925 422 : DEALLOCATE (xyz_mix)
2926 1712 : DO idim = 1, dim2
2927 1712 : IF (ns_me /= 0) THEN
2928 1266 : DEALLOCATE (xyz_mix_ns(idim)%c_array)
2929 : END IF
2930 : END DO
2931 422 : DEALLOCATE (xyz_mix_ns)
2932 422 : IF (ns_me /= 0) THEN
2933 414 : DEALLOCATE (gmat)
2934 : END IF
2935 422 : DEALLOCATE (mii)
2936 422 : DEALLOCATE (mij)
2937 422 : DEALLOCATE (mjj)
2938 422 : DEALLOCATE (list_pair)
2939 :
2940 844 : END SUBROUTINE jacobi_rot_para_core
2941 :
2942 : ! **************************************************************************************************
2943 : !> \brief ...
2944 : !> \param iperm ...
2945 : !> \param para_env ...
2946 : !> \param list_pair ...
2947 : ! **************************************************************************************************
2948 87316 : 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 87316 : npair = (para_env%num_pe + 1)/2
2956 87316 : IF (iperm == 1) THEN
2957 : !..set up initial ordering
2958 261948 : DO I = 0, para_env%num_pe - 1
2959 174632 : II = ((i + 1) + 1)/2
2960 174632 : JJ = MOD((i + 1) + 1, 2) + 1
2961 261948 : list_pair(JJ, II) = i
2962 : END DO
2963 87316 : IF (MOD(para_env%num_pe, 2) == 1) list_pair(2, npair) = -1
2964 0 : ELSE IF (MOD(iperm, 2) == 0) THEN
2965 : !..a type shift
2966 0 : jj = list_pair(1, npair)
2967 0 : DO I = npair, 3, -1
2968 0 : list_pair(1, I) = list_pair(1, I - 1)
2969 : END DO
2970 0 : list_pair(1, 2) = list_pair(2, 1)
2971 0 : list_pair(2, 1) = jj
2972 : ELSE
2973 : !..b type shift
2974 0 : jj = list_pair(2, 1)
2975 0 : DO I = 1, npair - 1
2976 0 : list_pair(2, I) = list_pair(2, I + 1)
2977 : END DO
2978 0 : list_pair(2, npair) = jj
2979 : END IF
2980 :
2981 87316 : END SUBROUTINE eberlein
2982 :
2983 : ! **************************************************************************************************
2984 : !> \brief ...
2985 : !> \param vectors ...
2986 : !> \param op_sm_set ...
2987 : !> \param zij_fm_set ...
2988 : ! **************************************************************************************************
2989 532 : 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 532 : CALL timeset(routineN, handle)
3001 :
3002 : ! get rows and cols of the input
3003 532 : CALL cp_fm_get_info(vectors, nrow_global=nao, ncol_global=nmoloc)
3004 : ! replicate the input kind of matrix
3005 532 : CALL cp_fm_create(opvec, vectors%matrix_struct)
3006 :
3007 : ! Compute zij here
3008 2164 : DO i = 1, SIZE(zij_fm_set, 2)
3009 5428 : DO j = 1, SIZE(zij_fm_set, 1)
3010 3264 : CALL cp_fm_set_all(zij_fm_set(j, i), 0.0_dp)
3011 3264 : 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 4896 : zij_fm_set(j, i))
3014 : END DO
3015 : END DO
3016 :
3017 532 : CALL cp_fm_release(opvec)
3018 532 : CALL timestop(handle)
3019 :
3020 532 : END SUBROUTINE zij_matrix
3021 :
3022 : ! **************************************************************************************************
3023 : !> \brief ...
3024 : !> \param vectors ...
3025 : ! **************************************************************************************************
3026 38 : 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 38 : REAL(KIND=dp), DIMENSION(:), POINTER :: tau
3034 : TYPE(cp_fm_struct_type), POINTER :: cstruct
3035 : TYPE(cp_fm_type) :: CTp, Qf, tmp
3036 :
3037 38 : CALL timeset(routineN, handle)
3038 :
3039 : ! Create Transpose of Coefficient Matrix vectors
3040 38 : nrowT = vectors%matrix_struct%ncol_global
3041 38 : ncolT = vectors%matrix_struct%nrow_global
3042 :
3043 : CALL cp_fm_struct_create(cstruct, template_fmstruct=vectors%matrix_struct, &
3044 38 : nrow_global=nrowT, ncol_global=ncolT)
3045 38 : CALL cp_fm_create(CTp, cstruct)
3046 38 : CALL cp_fm_struct_release(cstruct)
3047 :
3048 114 : ALLOCATE (tau(nrowT))
3049 :
3050 38 : CALL cp_fm_transpose(vectors, CTp)
3051 :
3052 : ! Get QR decomposition of CTs
3053 38 : 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 38 : ncol_global=CTp%matrix_struct%nrow_global)
3059 38 : CALL cp_fm_create(Qf, cstruct)
3060 38 : CALL cp_fm_struct_release(cstruct)
3061 38 : CALL cp_fm_to_fm_submat(CTp, Qf, nrowT, nrowT, 1, 1, 1, 1)
3062 :
3063 : ! Get Q
3064 38 : CALL cp_fm_pdorgqr(Qf, tau, nrowT, 1, 1)
3065 :
3066 : ! Transform original coefficient matrix vectors
3067 38 : CALL cp_fm_create(tmp, vectors%matrix_struct)
3068 38 : CALL cp_fm_set_all(tmp, 0.0_dp, 1.0_dp)
3069 38 : CALL cp_fm_to_fm(vectors, tmp)
3070 38 : CALL parallel_gemm('N', 'N', ncolT, nrowT, nrowT, 1.0_dp, tmp, Qf, 0.0_dp, vectors)
3071 :
3072 : ! Cleanup
3073 38 : CALL cp_fm_release(CTp)
3074 38 : CALL cp_fm_release(tmp)
3075 38 : CALL cp_fm_release(Qf)
3076 38 : DEALLOCATE (tau)
3077 :
3078 38 : CALL timestop(handle)
3079 :
3080 152 : 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 6 : 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 6 : 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 6 : CALL timeset(routineN, handle)
3118 :
3119 6 : dim1 = SIZE(zij, 1)
3120 6 : dim2 = SIZE(zij, 2)
3121 :
3122 6 : NULLIFY (c_rmat, c_zij)
3123 84 : ALLOCATE (c_rmat, c_zij(dim1*dim2), mii(dim1*dim2), mij(dim1*dim2), mji(dim1*dim2), mjj(dim1*dim2))
3124 6 : CALL cp_cfm_create(c_rmat, zij(1, 1)%matrix_struct)
3125 6 : 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 24 : DO idim = 1, dim2
3127 60 : DO jdim = 1, dim1
3128 36 : CALL cp_cfm_create(c_zij((idim - 1)*dim1 + jdim), zij(jdim, idim)%matrix_struct)
3129 54 : CALL cp_cfm_to_cfm(zij(jdim, idim), c_zij((idim - 1)*dim1 + jdim))
3130 : END DO
3131 : END DO
3132 :
3133 6 : CALL cp_cfm_get_info(c_rmat, nrow_global=nstate)
3134 6 : tolerance = 1.0e10_dp
3135 6 : old_spread = 1.0e10_dp
3136 6 : sweeps = 0
3137 6 : unit_nr = -1
3138 6 : IF (c_rmat%matrix_struct%para_env%is_source()) THEN
3139 3 : unit_nr = cp_logger_get_default_unit_nr()
3140 : WRITE (unit_nr, "(T4,A )") " Localization by iterative Jacobi rotation using "// &
3141 3 : "Cardoso-Souloumiac angles"
3142 :
3143 : END IF
3144 :
3145 : ! Jacobi sweeps until converged
3146 58 : DO WHILE (tolerance >= eps_localization .AND. sweeps < max_iter)
3147 52 : sweeps = sweeps + 1
3148 52 : t1 = m_walltime()
3149 1924 : DO jstate = 1, nstate
3150 34684 : DO istate = jstate + 1, nstate
3151 229320 : DO idim = 1, dim1*dim2
3152 196560 : CALL cp_cfm_get_element(c_zij(idim), istate, istate, mii(idim))
3153 196560 : CALL cp_cfm_get_element(c_zij(idim), istate, jstate, mij(idim))
3154 196560 : CALL cp_cfm_get_element(c_zij(idim), jstate, istate, mji(idim))
3155 229320 : CALL cp_cfm_get_element(c_zij(idim), jstate, jstate, mjj(idim))
3156 : END DO
3157 32760 : CALL get_cardoso_angles(mii, mij, mji, mjj, c, s, vectors%matrix_struct)
3158 229320 : DO idim = 1, dim1*dim2
3159 196560 : CALL cp_cfm_rot_cols(c_zij(idim), istate, jstate, c, REAL(s))
3160 229320 : CALL cp_cfm_rot_rows(c_zij(idim), istate, jstate, c, REAL(s))
3161 : END DO
3162 67392 : CALL cp_cfm_rot_cols(c_rmat, istate, jstate, c, REAL(s))
3163 : END DO
3164 : END DO
3165 52 : CALL check_tolerance(c_zij, weights, spread)
3166 52 : CALL check_tolerance(c_zij, weights, tolerance)
3167 52 : tolerance = ABS(spread - old_spread)
3168 52 : old_spread = spread
3169 :
3170 52 : t2 = m_walltime()
3171 58 : 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 26 : "Iteration:", sweeps, "Functional", spread, "Tolerance:", tolerance, "Time:", t2 - t1
3174 26 : CALL m_flush(unit_nr)
3175 : END IF
3176 : END DO
3177 :
3178 24 : DO idim = 1, dim2
3179 : ! back to an interlaced matrix
3180 60 : DO jdim = 1, dim1
3181 36 : CALL cp_cfm_to_cfm(c_zij((idim - 1)*dim1 + jdim), zij(jdim, idim))
3182 54 : CALL cp_cfm_release(c_zij((idim - 1)*dim1 + jdim))
3183 : END DO
3184 : END DO
3185 :
3186 6 : CALL rotate_orbitals_cfm(c_rmat, vectors)
3187 :
3188 6 : DEALLOCATE (c_zij, mii, mij, mji, mjj)
3189 6 : CALL cp_cfm_release(c_rmat)
3190 6 : DEALLOCATE (c_rmat)
3191 :
3192 6 : CALL timestop(handle)
3193 :
3194 18 : 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 6 : 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 6 : TYPE(cp_cfm_type), POINTER :: c_zij(:), rmat
3223 :
3224 6 : CALL timeset(routineN, handle)
3225 :
3226 6 : CALL cite_reference(Schreder2024_2)
3227 :
3228 18 : n_dim = SIZE(zij)
3229 :
3230 6 : ALLOCATE (rmat)
3231 6 : CALL cp_cfm_create(rmat, zij(1, 1)%matrix_struct)
3232 :
3233 6 : CALL cp_cfm_set_all(rmat, (0.0_dp, 0.0_dp), (1.0_dp, 0.0_dp))
3234 120 : ALLOCATE (c_zij(n_dim), qii(n_dim), qij(n_dim), qji(n_dim), qjj(n_dim))
3235 :
3236 : ! build Qij matrix
3237 78 : DO k = 1, n_dim
3238 78 : c_zij(k) = zij(k, 1)
3239 : END DO
3240 :
3241 6 : CALL cp_cfm_get_info(c_zij(1), ncol_global=n_states)
3242 :
3243 6 : tol = 1.0e10_dp
3244 6 : c_spread = (0.0_dp, 0.0_dp)
3245 78 : DO k = 1, n_dim
3246 2382 : DO i = 1, n_states
3247 2304 : CALL cp_cfm_get_element(c_zij(k), i, i, s)
3248 2376 : c_spread = c_spread + s*s
3249 : END DO
3250 : END DO
3251 6 : spread = REAL(c_spread)
3252 6 : old_spread = spread
3253 :
3254 6 : sweeps = 0
3255 6 : output_unit = cp_logger_get_default_unit_nr()
3256 : WRITE (output_unit, "(T4,A )") " Localization by iterative Jacobi rotation using "// &
3257 6 : "Cardoso-Souloumiac angles"
3258 6 : WRITE (output_unit, "(T4,A )") " and Pipek-Mezey spread functional"
3259 :
3260 6 : 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 6 : "Iteration:", sweeps, " Functional", spread, " Tolerance:", tol, " Time:", 0.0_dp
3263 : END IF
3264 :
3265 36 : DO WHILE (tol >= eps .AND. sweeps < max_iter)
3266 30 : t1 = m_walltime()
3267 30 : sweeps = sweeps + 1
3268 :
3269 990 : DO i = 1, n_states
3270 15870 : DO j = i + 1, n_states
3271 193440 : DO k = 1, n_dim
3272 178560 : CALL cp_cfm_get_element(c_zij(k), i, i, qii(k))
3273 178560 : CALL cp_cfm_get_element(c_zij(k), i, j, qij(k))
3274 178560 : CALL cp_cfm_get_element(c_zij(k), j, i, qji(k))
3275 193440 : CALL cp_cfm_get_element(c_zij(k), j, j, qjj(k))
3276 : END DO
3277 14880 : CALL get_cardoso_angles(qii, qij, qji, qjj, c, s, vec%matrix_struct)
3278 193440 : DO k = 1, n_dim
3279 178560 : CALL cp_cfm_rot_cols(c_zij(k), i, j, c, REAL(s))
3280 193440 : CALL cp_cfm_rot_rows(c_zij(k), i, j, c, REAL(s))
3281 : END DO
3282 30720 : CALL cp_cfm_rot_cols(rmat, i, j, c, REAL(s))
3283 : END DO
3284 : END DO
3285 :
3286 30 : c_spread = (0.0_dp, 0.0_dp)
3287 990 : DO i = 1, n_states
3288 12510 : DO k = 1, n_dim
3289 11520 : CALL cp_cfm_get_element(c_zij(k), i, i, s)
3290 12480 : c_spread = c_spread + s*s
3291 : END DO
3292 : END DO
3293 30 : spread = REAL(c_spread)
3294 :
3295 30 : tol = ABS(spread - old_spread)
3296 30 : old_spread = spread
3297 30 : t2 = m_walltime()
3298 36 : 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 30 : "Iteration:", sweeps, " Functional", spread, " Tolerance:", tol, " Time:", t2 - t1
3301 : END IF
3302 : END DO
3303 :
3304 6 : CALL rotate_orbitals_cfm(rmat, vec)
3305 :
3306 : ! c_zij(k) was assigned (shallow copy of zij(k,1)); do not release contents
3307 6 : DEALLOCATE (c_zij, qii, qij, qji, qjj)
3308 6 : CALL cp_cfm_release(rmat)
3309 6 : DEALLOCATE (rmat)
3310 :
3311 6 : CALL timestop(handle)
3312 :
3313 18 : 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 47640 : 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 47640 : dim_m = SIZE(mii)
3343 :
3344 : CALL cp_fm_struct_create(h_fm_struct, nrow_global=1, ncol_global=3, &
3345 47640 : template_fmstruct=tmp_fm_struct)
3346 : CALL cp_fm_struct_create(G_fm_struct, nrow_global=3, ncol_global=3, &
3347 47640 : template_fmstruct=tmp_fm_struct)
3348 47640 : ALLOCATE (hmat, c_Gmat, Gmat, evects)
3349 47640 : CALL cp_cfm_create(hmat, h_fm_struct)
3350 47640 : CALL cp_cfm_create(c_Gmat, G_fm_struct)
3351 47640 : CALL cp_cfm_set_all(c_Gmat, (0.0_dp, 0.0_dp), (0.0_dp, 0.0_dp))
3352 47640 : CALL cp_fm_create(Gmat, G_fm_struct)
3353 47640 : CALL cp_fm_create(evects, G_fm_struct)
3354 :
3355 422760 : DO i = 1, dim_m
3356 375120 : CALL cp_cfm_set_element(hmat, 1, 1, (mii(i) - mjj(i)))
3357 375120 : CALL cp_cfm_set_element(hmat, 1, 2, (mij(i) + mji(i)))
3358 375120 : CALL cp_cfm_set_element(hmat, 1, 3, (0.0_dp, 1.0_dp)*(mji(i) - mij(i)))
3359 422760 : 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 47640 : CALL cp_cfm_to_fm(c_Gmat, Gmat)
3362 :
3363 : ! find eigenvector with highest eigenvalue
3364 47640 : CALL choose_eigv_solver(Gmat, evects, evals)
3365 190560 : i_max = MAXLOC(evals, 1)
3366 47640 : CALL cp_fm_get_element(evects, 1, i_max, x)
3367 47640 : CALL cp_fm_get_element(evects, 2, i_max, y)
3368 47640 : CALL cp_fm_get_element(evects, 3, i_max, z)
3369 47640 : IF (x < 0) THEN
3370 47288 : x = -x
3371 47288 : y = -y
3372 47288 : z = -z
3373 : END IF
3374 :
3375 : ! calculate the angles
3376 47640 : r = SQRT(x**2 + y**2 + z**2) ! always 1
3377 47640 : c = SQRT((x + r)/(2*r)) ! always real
3378 47640 : s = (y - gaussi*z)/SQRT(2*r*(x + r)) ! always complex
3379 :
3380 47640 : CALL cp_fm_struct_release(h_fm_struct)
3381 47640 : CALL cp_fm_struct_release(G_fm_struct)
3382 47640 : CALL cp_cfm_release(hmat)
3383 47640 : CALL cp_cfm_release(c_Gmat)
3384 47640 : CALL cp_fm_release(Gmat)
3385 47640 : CALL cp_fm_release(evects)
3386 47640 : DEALLOCATE (hmat, c_Gmat, Gmat, evects)
3387 :
3388 47640 : END SUBROUTINE get_cardoso_angles
3389 :
3390 0 : END MODULE qs_localization_methods
|
