(git:374b731)
Loading...
Searching...
No Matches
dm_ls_scf_curvy.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief density matrix optimization using exponential transformations
10!> \par History
11!> 2012.05 created [Florian Schiffmann]
12!> \author Florian Schiffmann
13! **************************************************************************************************
14
15MODULE dm_ls_scf_curvy
16 USE bibliography, ONLY: shao2003,&
17 cite_reference
18 USE cp_log_handling, ONLY: cp_get_default_logger,&
19 cp_logger_get_default_unit_nr,&
20 cp_logger_type
21 USE dbcsr_api, ONLY: &
22 dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_dot, dbcsr_filter, dbcsr_frobenius_norm, &
23 dbcsr_multiply, dbcsr_norm, dbcsr_release, dbcsr_scale, dbcsr_set, dbcsr_transposed, &
24 dbcsr_type, dbcsr_type_no_symmetry
25 USE dm_ls_scf_types, ONLY: ls_scf_curvy_type,&
26 ls_scf_env_type
27 USE input_constants, ONLY: ls_scf_line_search_3point,&
28 ls_scf_line_search_3point_2d
29 USE iterate_matrix, ONLY: purify_mcweeny
30 USE kinds, ONLY: dp
31 USE machine, ONLY: m_flush
32 USE mathconstants, ONLY: ifac
33 USE mathlib, ONLY: invmat
34#include "./base/base_uses.f90"
35
36 IMPLICIT NONE
37
38 PRIVATE
39
40 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'dm_ls_scf_curvy'
41
42 PUBLIC :: dm_ls_curvy_optimization, deallocate_curvy_data
43
44CONTAINS
45
46! **************************************************************************************************
47!> \brief driver routine for Head-Gordon curvy step approach
48!> \param ls_scf_env ...
49!> \param energy ...
50!> \param check_conv ...
51!> \par History
52!> 2012.05 created [Florian Schiffmann]
53!> \author Florian Schiffmann
54! **************************************************************************************************
55
56 SUBROUTINE dm_ls_curvy_optimization(ls_scf_env, energy, check_conv)
57 TYPE(ls_scf_env_type) :: ls_scf_env
58 REAL(kind=dp) :: energy
59 LOGICAL :: check_conv
60
61 CHARACTER(LEN=*), PARAMETER :: routinen = 'dm_ls_curvy_optimization'
62
63 INTEGER :: handle, i, lsstep
64
65 CALL timeset(routinen, handle)
66
67 CALL cite_reference(shao2003)
68
69! Upon first call initialize all matrices needed curing optimization
70! In addition transform P into orthonormal basis. Will be scaled by 0.5 in closed shell case
71! Only to be done once as it will be stored and reused afterwards
72! TRS4 might yield a non-idempotent P therefore McWeeny purification is applied on initial P
73
74 IF (.NOT. ALLOCATED(ls_scf_env%curvy_data%matrix_dp)) THEN
75 CALL init_curvy(ls_scf_env%curvy_data, ls_scf_env%matrix_s, ls_scf_env%nspins)
76 ls_scf_env%curvy_data%line_search_step = 1
77
78 IF (ls_scf_env%curvy_data%line_search_type == ls_scf_line_search_3point_2d) THEN
79 DO i = 1, ls_scf_env%nspins
80 CALL dbcsr_copy(ls_scf_env%curvy_data%matrix_psave(i, 1), &
81 ls_scf_env%matrix_p(i))
82 END DO
83 END IF
84 IF (ls_scf_env%nspins == 1) CALL dbcsr_scale(ls_scf_env%matrix_p(1), 0.5_dp)
85 CALL transform_matrix_orth(ls_scf_env%matrix_p, ls_scf_env%matrix_s_sqrt, &
86 ls_scf_env%eps_filter)
87 CALL purify_mcweeny(ls_scf_env%matrix_p, ls_scf_env%eps_filter, 3)
88 DO i = 1, ls_scf_env%nspins
89 CALL dbcsr_copy(ls_scf_env%curvy_data%matrix_p(i), ls_scf_env%matrix_p(i))
90 END DO
91 END IF
92
93 lsstep = ls_scf_env%curvy_data%line_search_step
94
95! If new search direction has to be computed transform H into the orthnormal basis
96
97 IF (ls_scf_env%curvy_data%line_search_step == 1) &
98 CALL transform_matrix_orth(ls_scf_env%matrix_ks, ls_scf_env%matrix_s_sqrt_inv, &
99 ls_scf_env%eps_filter)
100
101! Set the energies for the line search and make sure to give the correct energy back to scf_main
102 ls_scf_env%curvy_data%energies(lsstep) = energy
103 IF (lsstep .NE. 1) energy = ls_scf_env%curvy_data%energies(1)
104
105! start the optimization by calling the driver routine or simply combine saved P(2D line search)
106 IF (lsstep .LE. 2) THEN
107 CALL optimization_step(ls_scf_env%curvy_data, ls_scf_env)
108 ELSE IF (lsstep == ls_scf_env%curvy_data%line_search_type) THEN
109! line_search type has the value appropriate to the number of energy calculations needed
110 CALL optimization_step(ls_scf_env%curvy_data, ls_scf_env)
111 ELSE
112 CALL new_p_from_save(ls_scf_env%matrix_p, ls_scf_env%curvy_data%matrix_psave, lsstep, &
113 ls_scf_env%curvy_data%double_step_size)
114 ls_scf_env%curvy_data%line_search_step = ls_scf_env%curvy_data%line_search_step + 1
115 CALL timestop(handle)
116 RETURN
117 END IF
118 lsstep = ls_scf_env%curvy_data%line_search_step
119
120! transform new density matrix back into nonorthonormal basis (again scaling might apply)
121
122 CALL transform_matrix_orth(ls_scf_env%matrix_p, ls_scf_env%matrix_s_sqrt_inv, &
123 ls_scf_env%eps_filter)
124 IF (ls_scf_env%nspins == 1) CALL dbcsr_scale(ls_scf_env%matrix_p(1), 2.0_dp)
125
126! P-matrices only need to be stored in case of 2D line search
127 IF (lsstep .LE. 3 .AND. ls_scf_env%curvy_data%line_search_type == ls_scf_line_search_3point_2d) THEN
128 DO i = 1, ls_scf_env%nspins
129 CALL dbcsr_copy(ls_scf_env%curvy_data%matrix_psave(i, lsstep), &
130 ls_scf_env%matrix_p(i))
131 END DO
132 END IF
133 check_conv = lsstep == 1
134
135 CALL timestop(handle)
136
137 END SUBROUTINE dm_ls_curvy_optimization
138
139! **************************************************************************************************
140!> \brief low level routine for Head-Gordons curvy step approach
141!> computes gradients, performs a cg and line search,
142!> and evaluates the BCH series to obtain the new P matrix
143!> \param curvy_data ...
144!> \param ls_scf_env ...
145!> \par History
146!> 2012.05 created [Florian Schiffmann]
147!> \author Florian Schiffmann
148! **************************************************************************************************
149
150 SUBROUTINE optimization_step(curvy_data, ls_scf_env)
151 TYPE(ls_scf_curvy_type) :: curvy_data
152 TYPE(ls_scf_env_type) :: ls_scf_env
153
154 CHARACTER(LEN=*), PARAMETER :: routinen = 'optimization_step'
155
156 INTEGER :: handle, ispin
157 REAL(kind=dp) :: filter, step_size(2)
158
159! Upon first line search step compute new search direction and apply CG if required
160
161 CALL timeset(routinen, handle)
162
163 IF (curvy_data%line_search_step == 1) THEN
164 curvy_data%step_size = maxval(curvy_data%step_size)
165 curvy_data%step_size = min(max(0.10_dp, 0.5_dp*abs(curvy_data%step_size(1))), 0.5_dp)
166! Dynamic eps_filter for newton steps
167 filter = max(ls_scf_env%eps_filter*curvy_data%min_filter, &
168 ls_scf_env%eps_filter*curvy_data%filter_factor)
169 CALL compute_direction_newton(curvy_data%matrix_p, ls_scf_env%matrix_ks, &
170 curvy_data%matrix_dp, filter, curvy_data%fix_shift, curvy_data%shift, &
171 curvy_data%cg_numer, curvy_data%cg_denom, curvy_data%min_shift)
172 curvy_data%filter_factor = curvy_data%scale_filter*curvy_data%filter_factor
173 step_size = curvy_data%step_size
174 curvy_data%BCH_saved = 0
175 ELSE IF (curvy_data%line_search_step == 2) THEN
176 step_size = curvy_data%step_size
177 IF (curvy_data%energies(1) - curvy_data%energies(2) .GT. 0.0_dp) THEN
178 curvy_data%step_size = curvy_data%step_size*2.0_dp
179 curvy_data%double_step_size = .true.
180 ELSE
181 curvy_data%step_size = curvy_data%step_size*0.5_dp
182 curvy_data%double_step_size = .false.
183 END IF
184 step_size = curvy_data%step_size
185 ELSE IF (curvy_data%line_search_step == ls_scf_line_search_3point_2d) THEN
186 CALL line_search_2d(curvy_data%energies, curvy_data%step_size)
187 step_size = curvy_data%step_size
188 ELSE IF (curvy_data%line_search_step == ls_scf_line_search_3point) THEN
189 CALL line_search_3pnt(curvy_data%energies, curvy_data%step_size)
190 step_size = curvy_data%step_size
191 END IF
192
193 CALL update_p_exp(curvy_data%matrix_p, ls_scf_env%matrix_p, curvy_data%matrix_dp, &
194 curvy_data%matrix_BCH, ls_scf_env%eps_filter, step_size, curvy_data%BCH_saved, &
195 curvy_data%n_bch_hist)
196
197! line_search type has the value appropriate to the numeber of energy calculations needed
198 curvy_data%line_search_step = mod(curvy_data%line_search_step, curvy_data%line_search_type) + 1
199 IF (curvy_data%line_search_step == 1) THEN
200 DO ispin = 1, SIZE(curvy_data%matrix_p)
201 CALL dbcsr_copy(curvy_data%matrix_p(ispin), ls_scf_env%matrix_p(ispin))
202 END DO
203 END IF
204 CALL timestop(handle)
205
206 END SUBROUTINE optimization_step
207
208! **************************************************************************************************
209!> \brief Perform a 6pnt-2D line search for spin polarized calculations.
210!> Fit a 2D parabolic function to 6 points
211!> \param energies ...
212!> \param step_size ...
213!> \par History
214!> 2012.05 created [Florian Schiffmann]
215!> \author Florian Schiffmann
216! **************************************************************************************************
217
218 SUBROUTINE line_search_2d(energies, step_size)
219 REAL(kind=dp) :: energies(6), step_size(2)
220
221 INTEGER :: info, unit_nr
222 REAL(kind=dp) :: e_pred, param(6), s1, s1sq, s2, s2sq, &
223 sys_lin_eq(6, 6), tmp_e, v1, v2
224 TYPE(cp_logger_type), POINTER :: logger
225
226 logger => cp_get_default_logger()
227 IF (energies(1) - energies(2) .LT. 0._dp) THEN
228 tmp_e = energies(2); energies(2) = energies(3); energies(3) = tmp_e
229 step_size = step_size*2.0_dp
230 END IF
231 IF (logger%para_env%is_source()) THEN
232 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
233 ELSE
234 unit_nr = -1
235 END IF
236 s1 = 0.5_dp*step_size(1); s2 = step_size(1); s1sq = s1**2; s2sq = s2**2
237 sys_lin_eq = 0.0_dp; sys_lin_eq(:, 6) = 1.0_dp
238 sys_lin_eq(2, 1) = s1sq; sys_lin_eq(2, 2) = s1sq; sys_lin_eq(2, 3) = s1sq; sys_lin_eq(2, 4) = s1; sys_lin_eq(2, 5) = s1
239 sys_lin_eq(3, 1) = s2sq; sys_lin_eq(3, 2) = s2sq; sys_lin_eq(3, 3) = s2sq; sys_lin_eq(3, 4) = s2; sys_lin_eq(3, 5) = s2
240 sys_lin_eq(4, 3) = s1sq; sys_lin_eq(4, 5) = s1
241 sys_lin_eq(5, 1) = s1sq; sys_lin_eq(5, 4) = s1
242 sys_lin_eq(6, 3) = s2sq; sys_lin_eq(6, 5) = s2
243
244 CALL invmat(sys_lin_eq, info)
245 param = matmul(sys_lin_eq, energies)
246 v1 = (param(2)*param(4))/(2.0_dp*param(1)) - param(5)
247 v2 = -(param(2)**2)/(2.0_dp*param(1)) + 2.0_dp*param(3)
248 step_size(2) = v1/v2
249 step_size(1) = (-param(2)*step_size(2) - param(4))/(2.0_dp*param(1))
250 IF (step_size(1) .LT. 0.0_dp) step_size(1) = 1.0_dp
251 IF (step_size(2) .LT. 0.0_dp) step_size(2) = 1.0_dp
252! step_size(1)=MIN(step_size(1),2.0_dp)
253! step_size(2)=MIN(step_size(2),2.0_dp)
254 e_pred = param(1)*step_size(1)**2 + param(2)*step_size(1)*step_size(2) + &
255 param(3)*step_size(2)**2 + param(4)*step_size(1) + param(5)*step_size(2) + param(6)
256 IF (unit_nr .GT. 0) WRITE (unit_nr, "(t3,a,F10.5,F10.5,A,F20.9)") &
257 " Line Search: Step Size", step_size, " Predicted energy", e_pred
258 e_pred = param(1)*s1**2 + param(2)*s2*s1*0.0_dp + &
259 param(3)*s1**2*0.0_dp + param(4)*s1 + param(5)*s1*0.0_dp + param(6)
260
261 END SUBROUTINE line_search_2d
262
263! **************************************************************************************************
264!> \brief Perform a 3pnt line search
265!> \param energies ...
266!> \param step_size ...
267!> \par History
268!> 2012.05 created [Florian Schiffmann]
269!> \author Florian Schiffmann
270! **************************************************************************************************
271
272 SUBROUTINE line_search_3pnt(energies, step_size)
273 REAL(kind=dp) :: energies(3), step_size(2)
274
275 INTEGER :: unit_nr
276 REAL(kind=dp) :: a, b, c, e_pred, min_val, step1, tmp, &
277 tmp_e
278 TYPE(cp_logger_type), POINTER :: logger
279
280 logger => cp_get_default_logger()
281 IF (energies(1) - energies(2) .LT. 0._dp) THEN
282 tmp_e = energies(2); energies(2) = energies(3); energies(3) = tmp_e
283 step_size = step_size*2.0_dp
284 END IF
285 IF (logger%para_env%is_source()) THEN
286 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
287 ELSE
288 unit_nr = -1
289 END IF
290 step1 = 0.5_dp*step_size(1)
291 c = energies(1)
292 a = (energies(3) + c - 2.0_dp*energies(2))/(2.0_dp*step1**2)
293 b = (energies(2) - c - a*step1**2)/step1
294 IF (a .LT. 1.0e-12_dp) a = -1.0e-12_dp
295 min_val = -b/(2.0_dp*a)
296 e_pred = a*min_val**2 + b*min_val + c
297 tmp = step_size(1)
298 IF (e_pred .LT. energies(1) .AND. e_pred .LT. energies(2)) THEN
299 step_size = max(-1.0_dp, &
300 min(min_val, 10_dp*step_size))
301 ELSE
302 step_size = 1.0_dp
303 END IF
304 e_pred = a*(step_size(1))**2 + b*(step_size(1)) + c
305 IF (unit_nr .GT. 0) THEN
306 WRITE (unit_nr, "(t3,a,f16.8,a,F20.9)") "Line Search: Step Size", step_size(1), " Predicted energy", e_pred
307 CALL m_flush(unit_nr)
308 END IF
309 END SUBROUTINE line_search_3pnt
310
311! **************************************************************************************************
312!> \brief Get a new search direction. Iterate to obtain a Newton like step
313!> Refine with a CG update of the search direction
314!> \param matrix_p ...
315!> \param matrix_ks ...
316!> \param matrix_dp ...
317!> \param eps_filter ...
318!> \param fix_shift ...
319!> \param curvy_shift ...
320!> \param cg_numer ...
321!> \param cg_denom ...
322!> \param min_shift ...
323!> \par History
324!> 2012.05 created [Florian Schiffmann]
325!> \author Florian Schiffmann
326! **************************************************************************************************
327
328 SUBROUTINE compute_direction_newton(matrix_p, matrix_ks, matrix_dp, eps_filter, fix_shift, &
329 curvy_shift, cg_numer, cg_denom, min_shift)
330 TYPE(dbcsr_type), DIMENSION(:) :: matrix_p, matrix_ks, matrix_dp
331 REAL(kind=dp) :: eps_filter
332 LOGICAL :: fix_shift(2)
333 REAL(kind=dp) :: curvy_shift(2), cg_numer(2), &
334 cg_denom(2), min_shift
335
336 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_direction_newton'
337
338 INTEGER :: handle, i, ispin, ncyc, nspin, unit_nr
339 LOGICAL :: at_limit
340 REAL(kind=dp) :: beta, conv_val, maxel, old_conv, shift
341 TYPE(cp_logger_type), POINTER :: logger
342 TYPE(dbcsr_type) :: matrix_ax, matrix_b, matrix_cg, &
343 matrix_dp_old, matrix_pks, matrix_res, &
344 matrix_tmp, matrix_tmp1
345
346 logger => cp_get_default_logger()
347
348 IF (logger%para_env%is_source()) THEN
349 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
350 ELSE
351 unit_nr = -1
352 END IF
353 CALL timeset(routinen, handle)
354 nspin = SIZE(matrix_p)
355
356 CALL dbcsr_create(matrix_pks, template=matrix_dp(1), matrix_type=dbcsr_type_no_symmetry)
357 CALL dbcsr_create(matrix_ax, template=matrix_dp(1), matrix_type=dbcsr_type_no_symmetry)
358 CALL dbcsr_create(matrix_tmp, template=matrix_dp(1), matrix_type=dbcsr_type_no_symmetry)
359 CALL dbcsr_create(matrix_tmp1, template=matrix_dp(1), matrix_type=dbcsr_type_no_symmetry)
360 CALL dbcsr_create(matrix_res, template=matrix_dp(1), matrix_type=dbcsr_type_no_symmetry)
361 CALL dbcsr_create(matrix_cg, template=matrix_dp(1), matrix_type=dbcsr_type_no_symmetry)
362 CALL dbcsr_create(matrix_b, template=matrix_dp(1), matrix_type=dbcsr_type_no_symmetry)
363 CALL dbcsr_create(matrix_dp_old, template=matrix_dp(1), matrix_type=dbcsr_type_no_symmetry)
364
365 DO ispin = 1, nspin
366 CALL dbcsr_copy(matrix_dp_old, matrix_dp(ispin))
367
368! Precompute some matrices to save work during iterations
369 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_p(ispin), matrix_ks(ispin), &
370 0.0_dp, matrix_pks, filter_eps=eps_filter)
371 CALL dbcsr_transposed(matrix_b, matrix_pks)
372 CALL dbcsr_copy(matrix_cg, matrix_b)
373
374! Starting CG with guess 0-matrix gives -2*gradient=[Ks*P-(Ks*P)T] for cg_matrix in second step
375 CALL dbcsr_add(matrix_cg, matrix_pks, 2.0_dp, -2.0_dp)
376
377! Residual matrix in first step=cg matrix. Keep Pks for later use in CG!
378 CALL dbcsr_copy(matrix_res, matrix_cg)
379
380! Precompute -FP-[FP]T which will be used throughout the CG iterations
381 CALL dbcsr_add(matrix_b, matrix_pks, -1.0_dp, -1.0_dp)
382
383! Setup some values to check convergence and safety checks for eigenvalue shifting
384 CALL dbcsr_norm(matrix_res, which_norm=2, norm_scalar=old_conv)
385 old_conv = dbcsr_frobenius_norm(matrix_res)
386 shift = min(10.0_dp, max(min_shift, 0.05_dp*old_conv))
387 conv_val = max(0.010_dp*old_conv, 100.0_dp*eps_filter)
388 old_conv = 100.0_dp
389 IF (fix_shift(ispin)) THEN
390 shift = max(min_shift, min(10.0_dp, max(shift, curvy_shift(ispin) - 0.5_dp*curvy_shift(ispin))))
391 curvy_shift(ispin) = shift
392 END IF
393
394! Begin the real optimization loop
395 CALL dbcsr_set(matrix_dp(ispin), 0.0_dp)
396 ncyc = 10
397 DO i = 1, ncyc
398
399! One step to compute: -FPD-DPF-DFP-PFD (not obvious but symmetry allows for some tricks)
400 CALL commutator_symm(matrix_b, matrix_cg, matrix_ax, eps_filter, 1.0_dp)
401
402! Compute the missing bits 2*(FDP+PDF) (again use symmetry to compute as a commutator)
403 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_cg, matrix_p(ispin), &
404 0.0_dp, matrix_tmp, filter_eps=eps_filter)
405 CALL commutator_symm(matrix_ks(ispin), matrix_tmp, matrix_tmp1, eps_filter, 2.0_dp)
406 CALL dbcsr_add(matrix_ax, matrix_tmp1, 1.0_dp, 1.0_dp)
407
408! Apply the shift and hope it's enough to stabilize the CG iterations
409 CALL dbcsr_add(matrix_ax, matrix_cg, 1.0_dp, shift)
410
411 CALL compute_cg_matrices(matrix_ax, matrix_res, matrix_cg, matrix_dp(ispin), &
412 matrix_tmp, eps_filter, at_limit)
413 CALL dbcsr_filter(matrix_cg, eps_filter)
414
415! check for convergence of the newton step
416 maxel = dbcsr_frobenius_norm(matrix_res)
417 IF (unit_nr .GT. 0) THEN
418 WRITE (unit_nr, "(T3,A,F12.6)") "Convergence of Newton iteration ", maxel
419 CALL m_flush(unit_nr)
420 END IF
421 at_limit = at_limit .OR. (old_conv/maxel .LT. 1.01_dp)
422 old_conv = maxel
423 IF (i == ncyc .AND. maxel/conv_val .GT. 5.0_dp) THEN
424 fix_shift(ispin) = .true.
425 curvy_shift(ispin) = 4.0_dp*shift
426 END IF
427 IF (maxel .LT. conv_val .OR. at_limit) EXIT
428 END DO
429
430! Refine the Newton like search direction with a preconditioned cg update
431 CALL dbcsr_transposed(matrix_b, matrix_pks)
432 !compute b= -2*KsP+2*PKs=-(2*gradient)
433 CALL dbcsr_copy(matrix_cg, matrix_b)
434 CALL dbcsr_add(matrix_cg, matrix_pks, 1.0_dp, -1.0_dp)
435 cg_denom(ispin) = cg_numer(ispin)
436 CALL dbcsr_dot(matrix_cg, matrix_dp(ispin), cg_numer(ispin))
437 beta = cg_numer(ispin)/max(cg_denom(ispin), 1.0e-6_dp)
438 IF (beta .LT. 1.0_dp) THEN
439 beta = max(0.0_dp, beta)
440 CALL dbcsr_add(matrix_dp(ispin), matrix_dp_old, 1.0_dp, beta)
441 END IF
442 IF (unit_nr .GT. 0) WRITE (unit_nr, "(A)") " "
443 END DO
444
445 CALL dbcsr_release(matrix_pks)
446 CALL dbcsr_release(matrix_dp_old)
447 CALL dbcsr_release(matrix_b)
448 CALL dbcsr_release(matrix_ax)
449 CALL dbcsr_release(matrix_tmp)
450 CALL dbcsr_release(matrix_tmp1)
451 CALL dbcsr_release(matrix_b)
452 CALL dbcsr_release(matrix_res)
453 CALL dbcsr_release(matrix_cg)
454
455 IF (unit_nr .GT. 0) CALL m_flush(unit_nr)
456 CALL timestop(handle)
457 END SUBROUTINE compute_direction_newton
458
459! **************************************************************************************************
460!> \brief compute the optimal step size of the current cycle and update the
461!> matrices needed to solve the system of linear equations
462!> \param Ax ...
463!> \param res ...
464!> \param cg ...
465!> \param deltp ...
466!> \param tmp ...
467!> \param eps_filter ...
468!> \param at_limit ...
469!> \par History
470!> 2012.05 created [Florian Schiffmann]
471!> \author Florian Schiffmann
472! **************************************************************************************************
473
474 SUBROUTINE compute_cg_matrices(Ax, res, cg, deltp, tmp, eps_filter, at_limit)
475 TYPE(dbcsr_type) :: ax, res, cg, deltp, tmp
476 REAL(kind=dp) :: eps_filter
477 LOGICAL :: at_limit
478
479 INTEGER :: i, info
480 REAL(kind=dp) :: alpha, beta, devi(3), fac, fac1, &
481 lin_eq(3, 3), new_norm, norm_ca, &
482 norm_rr, vec(3)
483
484 at_limit = .false.
485 CALL dbcsr_dot(res, res, norm_rr)
486 CALL dbcsr_dot(cg, ax, norm_ca)
487 lin_eq = 0.0_dp
488 fac = norm_rr/norm_ca
489 fac1 = fac
490! Use a 3point line search and a fit to a quadratic function to determine optimal step size
491 DO i = 1, 3
492 CALL dbcsr_copy(tmp, res)
493 CALL dbcsr_add(tmp, ax, 1.0_dp, -fac)
494 devi(i) = dbcsr_frobenius_norm(tmp)
495 lin_eq(i, :) = (/fac**2, fac, 1.0_dp/)
496 fac = fac1 + fac1*((-1)**i)*0.5_dp
497 END DO
498 CALL invmat(lin_eq, info)
499 vec = matmul(lin_eq, devi)
500 alpha = -vec(2)/(2.0_dp*vec(1))
501 fac = sqrt(norm_rr/(norm_ca*alpha))
502!scale the previous matrices to match the step size
503 CALL dbcsr_scale(ax, fac)
504 CALL dbcsr_scale(cg, fac)
505 norm_ca = norm_ca*fac**2
506
507! USe CG to get the new matrices
508 alpha = norm_rr/norm_ca
509 CALL dbcsr_add(res, ax, 1.0_dp, -alpha)
510 CALL dbcsr_dot(res, res, new_norm)
511 IF (norm_rr .LT. eps_filter*0.001_dp .OR. new_norm .LT. eps_filter*0.001_dp) THEN
512 beta = 0.0_dp
513 at_limit = .true.
514 ELSE
515 beta = new_norm/norm_rr
516 CALL dbcsr_add(deltp, cg, 1.0_dp, alpha)
517 END IF
518 beta = new_norm/norm_rr
519 CALL dbcsr_add(cg, res, beta, 1.0_dp)
520
521 END SUBROUTINE compute_cg_matrices
522
523! **************************************************************************************************
524!> \brief Only for 2D line search. Use saved P-components to construct new
525!> test density matrix. Takes care as well, whether step_size
526!> increased or decreased during 2nd step and combines matrices accordingly
527!> \param matrix_p ...
528!> \param matrix_psave ...
529!> \param lsstep ...
530!> \param DOUBLE ...
531!> \par History
532!> 2012.05 created [Florian Schiffmann]
533!> \author Florian Schiffmann
534! **************************************************************************************************
535
536 SUBROUTINE new_p_from_save(matrix_p, matrix_psave, lsstep, DOUBLE)
537 TYPE(dbcsr_type), DIMENSION(:) :: matrix_p
538 TYPE(dbcsr_type), DIMENSION(:, :) :: matrix_psave
539 INTEGER :: lsstep
540 LOGICAL :: double
541
542 SELECT CASE (lsstep)
543 CASE (3)
544 CALL dbcsr_copy(matrix_p(1), matrix_psave(1, 1))
545 IF (double) THEN
546 CALL dbcsr_copy(matrix_p(2), matrix_psave(2, 2))
547 ELSE
548 CALL dbcsr_copy(matrix_p(2), matrix_psave(2, 3))
549 END IF
550 CASE (4)
551 IF (double) THEN
552 CALL dbcsr_copy(matrix_p(1), matrix_psave(1, 2))
553 ELSE
554 CALL dbcsr_copy(matrix_p(1), matrix_psave(1, 3))
555 END IF
556 CALL dbcsr_copy(matrix_p(2), matrix_psave(2, 1))
557 CASE (5)
558 CALL dbcsr_copy(matrix_p(1), matrix_psave(1, 1))
559 IF (double) THEN
560 CALL dbcsr_copy(matrix_p(2), matrix_psave(2, 3))
561 ELSE
562 CALL dbcsr_copy(matrix_p(2), matrix_psave(2, 2))
563 END IF
564 END SELECT
565
566 END SUBROUTINE new_p_from_save
567
568! **************************************************************************************************
569!> \brief computes a commutator exploiting symmetry RES=k*[A,B]=k*[AB-(AB)T]
570!> \param a ...
571!> \param b ...
572!> \param res ...
573!> \param eps_filter filtering threshold for sparse matrices
574!> \param prefac prefactor k in above equation
575!> \par History
576!> 2012.05 created [Florian Schiffmann]
577!> \author Florian Schiffmann
578! **************************************************************************************************
579
580 SUBROUTINE commutator_symm(a, b, res, eps_filter, prefac)
581 TYPE(dbcsr_type) :: a, b, res
582 REAL(kind=dp) :: eps_filter, prefac
583
584 CHARACTER(LEN=*), PARAMETER :: routinen = 'commutator_symm'
585
586 INTEGER :: handle
587 TYPE(dbcsr_type) :: work
588
589 CALL timeset(routinen, handle)
590
591 CALL dbcsr_create(work, template=a, matrix_type=dbcsr_type_no_symmetry)
592
593 CALL dbcsr_multiply("N", "N", prefac, a, b, 0.0_dp, res, filter_eps=eps_filter)
594 CALL dbcsr_transposed(work, res)
595 CALL dbcsr_add(res, work, 1.0_dp, -1.0_dp)
596
597 CALL dbcsr_release(work)
598
599 CALL timestop(handle)
600 END SUBROUTINE commutator_symm
601
602! **************************************************************************************************
603!> \brief Use the BCH update to get the new idempotent P
604!> Numerics don't allow for perfect idempotency, therefore a mc weeny
605!> step is used to make sure we stay close to the idempotent surface
606!> \param matrix_p_in ...
607!> \param matrix_p_out ...
608!> \param matrix_dp ...
609!> \param matrix_BCH ...
610!> \param threshold ...
611!> \param step_size ...
612!> \param BCH_saved ...
613!> \param n_bch_hist ...
614!> \par History
615!> 2012.05 created [Florian Schiffmann]
616!> \author Florian Schiffmann
617! **************************************************************************************************
618
619 SUBROUTINE update_p_exp(matrix_p_in, matrix_p_out, matrix_dp, matrix_BCH, threshold, step_size, &
620 BCH_saved, n_bch_hist)
621 TYPE(dbcsr_type), DIMENSION(:) :: matrix_p_in, matrix_p_out, matrix_dp
622 TYPE(dbcsr_type), DIMENSION(:, :) :: matrix_bch
623 REAL(kind=dp) :: threshold, step_size(2)
624 INTEGER :: bch_saved(2), n_bch_hist
625
626 CHARACTER(LEN=*), PARAMETER :: routinen = 'update_p_exp'
627
628 INTEGER :: handle, i, ispin, nsave, nspin, unit_nr
629 LOGICAL :: save_bch
630 REAL(kind=dp) :: frob_norm, step_fac
631 TYPE(cp_logger_type), POINTER :: logger
632 TYPE(dbcsr_type) :: matrix, matrix_tmp
633
634 CALL timeset(routinen, handle)
635
636 logger => cp_get_default_logger()
637 IF (logger%para_env%is_source()) THEN
638 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
639 ELSE
640 unit_nr = -1
641 END IF
642
643 CALL dbcsr_create(matrix, template=matrix_p_in(1), matrix_type=dbcsr_type_no_symmetry)
644 CALL dbcsr_create(matrix_tmp, template=matrix_p_in(1), matrix_type=dbcsr_type_no_symmetry)
645 nspin = SIZE(matrix_p_in)
646
647 DO ispin = 1, nspin
648 step_fac = 1.0_dp
649 frob_norm = 1.0_dp
650 nsave = 0
651
652 CALL dbcsr_copy(matrix_tmp, matrix_p_in(ispin))
653 CALL dbcsr_copy(matrix_p_out(ispin), matrix_p_in(ispin))
654! If a BCH history is used make good use of it and do a few steps as a copy and scale update of P
655! else BCH_saved will be 0 and loop is skipped
656 DO i = 1, bch_saved(ispin)
657 step_fac = step_fac*step_size(ispin)
658 CALL dbcsr_copy(matrix_tmp, matrix_p_out(ispin))
659 CALL dbcsr_add(matrix_p_out(ispin), matrix_bch(ispin, i), 1.0_dp, ifac(i)*step_fac)
660 CALL dbcsr_add(matrix_tmp, matrix_p_out(ispin), 1.0_dp, -1.0_dp)
661 frob_norm = dbcsr_frobenius_norm(matrix_tmp)
662 IF (unit_nr .GT. 0) WRITE (unit_nr, "(t3,a,i3,a,f16.8)") "BCH: step", i, " Norm of P_old-Pnew:", frob_norm
663 IF (frob_norm .LT. threshold) EXIT
664 END DO
665 IF (frob_norm .LT. threshold) cycle
666
667! If the copy and scale isn't enough compute a few more BCH steps. 20 seems high but except of the first step it will never be close
668 save_bch = bch_saved(ispin) == 0 .AND. n_bch_hist .GT. 0
669 DO i = bch_saved(ispin) + 1, 20
670 step_fac = step_fac*step_size(ispin)
671 !allow for a bit of matrix magic here by exploiting matrix and matrix_tmp
672 !matrix_tmp is alway the previous order of the BCH series
673 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_tmp, matrix_dp(ispin), &
674 0.0_dp, matrix, filter_eps=threshold)
675
676 !(anti)symmetry allows to sum the transposed instead of the full commutator, matrix becomes the latest result
677
678 CALL dbcsr_transposed(matrix_tmp, matrix)
679 CALL dbcsr_add(matrix, matrix_tmp, 1.0_dp, 1.0_dp)
680
681 !Finally, add the new BCH order to P, but store the previous one for a convergence check
682 CALL dbcsr_copy(matrix_tmp, matrix_p_out(ispin))
683 CALL dbcsr_add(matrix_p_out(ispin), matrix, 1.0_dp, ifac(i)*step_fac)
684 IF (save_bch .AND. i .LE. n_bch_hist) THEN
685 CALL dbcsr_copy(matrix_bch(ispin, i), matrix)
686 nsave = i
687 END IF
688
689 CALL dbcsr_add(matrix_tmp, matrix_p_out(ispin), 1.0_dp, -1.0_dp)
690
691 !Stop the BCH-series if two successive P's differ by less the threshold
692 frob_norm = dbcsr_frobenius_norm(matrix_tmp)
693 IF (unit_nr .GT. 0) WRITE (unit_nr, "(t3,a,i3,a,f16.8)") "BCH: step", i, " Norm of P_old-Pnew:", frob_norm
694 IF (frob_norm .LT. threshold) EXIT
695
696 !Copy the latest BCH-matrix on matrix tmp, so we can cycle with all matrices in place
697 CALL dbcsr_copy(matrix_tmp, matrix)
698 CALL dbcsr_filter(matrix_tmp, threshold)
699 END DO
700 bch_saved(ispin) = nsave
701 IF (unit_nr .GT. 0) WRITE (unit_nr, "(A)") " "
702 END DO
703
704 CALL purify_mcweeny(matrix_p_out, threshold, 1)
705 IF (unit_nr .GT. 0) CALL m_flush(unit_nr)
706 CALL dbcsr_release(matrix_tmp)
707 CALL dbcsr_release(matrix)
708 CALL timestop(handle)
709 END SUBROUTINE update_p_exp
710
711! **************************************************************************************************
712!> \brief performs a transformation of a matrix back to/into orthonormal basis
713!> in case of P a scaling of 0.5 has to be applied for closed shell case
714!> \param matrix matrix to be transformed
715!> \param matrix_trafo transformation matrix
716!> \param eps_filter filtering threshold for sparse matrices
717!> \par History
718!> 2012.05 created [Florian Schiffmann]
719!> \author Florian Schiffmann
720! **************************************************************************************************
721
722 SUBROUTINE transform_matrix_orth(matrix, matrix_trafo, eps_filter)
723 TYPE(dbcsr_type), DIMENSION(:) :: matrix
724 TYPE(dbcsr_type) :: matrix_trafo
725 REAL(kind=dp) :: eps_filter
726
727 CHARACTER(LEN=*), PARAMETER :: routinen = 'transform_matrix_orth'
728
729 INTEGER :: handle, ispin
730 TYPE(dbcsr_type) :: matrix_tmp, matrix_work
731
732 CALL timeset(routinen, handle)
733
734 CALL dbcsr_create(matrix_work, template=matrix(1), matrix_type=dbcsr_type_no_symmetry)
735 CALL dbcsr_create(matrix_tmp, template=matrix(1), matrix_type=dbcsr_type_no_symmetry)
736
737 DO ispin = 1, SIZE(matrix)
738 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix(ispin), matrix_trafo, &
739 0.0_dp, matrix_work, filter_eps=eps_filter)
740 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_trafo, matrix_work, &
741 0.0_dp, matrix_tmp, filter_eps=eps_filter)
742 ! symmetrize results (this is again needed to make sure everything is stable)
743 CALL dbcsr_transposed(matrix_work, matrix_tmp)
744 CALL dbcsr_add(matrix_tmp, matrix_work, 0.5_dp, 0.5_dp)
745 CALL dbcsr_copy(matrix(ispin), matrix_tmp)
746 END DO
747
748 CALL dbcsr_release(matrix_tmp)
749 CALL dbcsr_release(matrix_work)
750 CALL timestop(handle)
751
752 END SUBROUTINE
753
754! **************************************************************************************************
755!> \brief ...
756!> \param curvy_data ...
757! **************************************************************************************************
758 SUBROUTINE deallocate_curvy_data(curvy_data)
759 TYPE(ls_scf_curvy_type) :: curvy_data
760
761 INTEGER :: i, j
762
763 CALL release_dbcsr_array(curvy_data%matrix_dp)
764 CALL release_dbcsr_array(curvy_data%matrix_p)
765
766 IF (ALLOCATED(curvy_data%matrix_psave)) THEN
767 DO i = 1, SIZE(curvy_data%matrix_psave, 1)
768 DO j = 1, 3
769 CALL dbcsr_release(curvy_data%matrix_psave(i, j))
770 END DO
771 END DO
772 DEALLOCATE (curvy_data%matrix_psave)
773 END IF
774 IF (ALLOCATED(curvy_data%matrix_BCH)) THEN
775 DO i = 1, SIZE(curvy_data%matrix_BCH, 1)
776 DO j = 1, 7
777 CALL dbcsr_release(curvy_data%matrix_BCH(i, j))
778 END DO
779 END DO
780 DEALLOCATE (curvy_data%matrix_BCH)
781 END IF
782 END SUBROUTINE deallocate_curvy_data
783
784! **************************************************************************************************
785!> \brief ...
786!> \param matrix ...
787! **************************************************************************************************
788 SUBROUTINE release_dbcsr_array(matrix)
789 TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:) :: matrix
790
791 INTEGER :: i
792
793 IF (ALLOCATED(matrix)) THEN
794 DO i = 1, SIZE(matrix)
795 CALL dbcsr_release(matrix(i))
796 END DO
797 DEALLOCATE (matrix)
798 END IF
799 END SUBROUTINE release_dbcsr_array
800
801! **************************************************************************************************
802!> \brief ...
803!> \param curvy_data ...
804!> \param matrix_s ...
805!> \param nspins ...
806! **************************************************************************************************
807 SUBROUTINE init_curvy(curvy_data, matrix_s, nspins)
808 TYPE(ls_scf_curvy_type) :: curvy_data
809 TYPE(dbcsr_type) :: matrix_s
810 INTEGER :: nspins
811
812 INTEGER :: ispin, j
813
814 ALLOCATE (curvy_data%matrix_dp(nspins))
815 ALLOCATE (curvy_data%matrix_p(nspins))
816 DO ispin = 1, nspins
817 CALL dbcsr_create(curvy_data%matrix_dp(ispin), template=matrix_s, &
818 matrix_type=dbcsr_type_no_symmetry)
819 CALL dbcsr_set(curvy_data%matrix_dp(ispin), 0.0_dp)
820 CALL dbcsr_create(curvy_data%matrix_p(ispin), template=matrix_s, &
821 matrix_type=dbcsr_type_no_symmetry)
822 curvy_data%fix_shift = .false.
823 curvy_data%double_step_size = .true.
824 curvy_data%shift = 1.0_dp
825 curvy_data%BCH_saved = 0
826 curvy_data%step_size = 0.60_dp
827 curvy_data%cg_numer = 0.00_dp
828 curvy_data%cg_denom = 0.00_dp
829 END DO
830 IF (curvy_data%line_search_type == ls_scf_line_search_3point_2d) THEN
831 ALLOCATE (curvy_data%matrix_psave(nspins, 3))
832 DO ispin = 1, nspins
833 DO j = 1, 3
834 CALL dbcsr_create(curvy_data%matrix_psave(ispin, j), template=matrix_s, &
835 matrix_type=dbcsr_type_no_symmetry)
836 END DO
837 END DO
838 END IF
839 IF (curvy_data%n_bch_hist .GT. 0) THEN
840 ALLOCATE (curvy_data%matrix_BCH(nspins, curvy_data%n_bch_hist))
841 DO ispin = 1, nspins
842 DO j = 1, curvy_data%n_bch_hist
843 CALL dbcsr_create(curvy_data%matrix_BCH(ispin, j), template=matrix_s, &
844 matrix_type=dbcsr_type_no_symmetry)
845 END DO
846 END DO
847 END IF
848
849 END SUBROUTINE init_curvy
850
851END MODULE dm_ls_scf_curvy
iterate_matrix::purify_mcweeny
Definition iterate_matrix.F:53
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition bibliography.F:28
bibliography::shao2003
integer, save, public shao2003
Definition bibliography.F:43
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_unit_nr
recursive integer function, public cp_logger_get_default_unit_nr(logger, local, skip_not_ionode)
asks the default unit number of the given logger. try to use cp_logger_get_unit_nr
Definition cp_log_handling.F:567
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
dm_ls_scf_curvy
density matrix optimization using exponential transformations
Definition dm_ls_scf_curvy.F:15
dm_ls_scf_curvy::dm_ls_curvy_optimization
subroutine, public dm_ls_curvy_optimization(ls_scf_env, energy, check_conv)
driver routine for Head-Gordon curvy step approach
Definition dm_ls_scf_curvy.F:57
dm_ls_scf_curvy::deallocate_curvy_data
subroutine, public deallocate_curvy_data(curvy_data)
...
Definition dm_ls_scf_curvy.F:759
dm_ls_scf_types
Types needed for a linear scaling quickstep SCF run based on the density matrix.
Definition dm_ls_scf_types.F:15
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::ls_scf_line_search_3point
integer, parameter, public ls_scf_line_search_3point
Definition input_constants.F:968
input_constants::ls_scf_line_search_3point_2d
integer, parameter, public ls_scf_line_search_3point_2d
Definition input_constants.F:968
iterate_matrix
Routines useful for iterative matrix calculations.
Definition iterate_matrix.F:13
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
machine
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
machine::m_flush
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition machine.F:106
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::ifac
real(kind=dp), dimension(0:maxfac), parameter, public ifac
Definition mathconstants.F:49
mathconstants::fac
real(kind=dp), dimension(0:maxfac), parameter, public fac
Definition mathconstants.F:37
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::invmat
subroutine, public invmat(a, info)
returns inverse of matrix using the lapack routines DGETRF and DGETRI
Definition mathlib.F:543
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
dm_ls_scf_types::ls_scf_curvy_type
Definition dm_ls_scf_types.F:55
dm_ls_scf_types::ls_scf_env_type
Definition dm_ls_scf_types.F:88