(git:ed6f26b)
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-2025 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_dbcsr_api, ONLY: &
19 dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_filter, dbcsr_multiply, dbcsr_release, &
20 dbcsr_scale, dbcsr_set, dbcsr_transposed, dbcsr_type, dbcsr_type_no_symmetry
21 USE cp_dbcsr_contrib, ONLY: dbcsr_dot,&
22 dbcsr_frobenius_norm
23 USE cp_log_handling, ONLY: cp_get_default_logger,&
24 cp_logger_get_default_unit_nr,&
25 cp_logger_type
26 USE dm_ls_scf_types, ONLY: ls_scf_curvy_type,&
27 ls_scf_env_type
28 USE input_constants, ONLY: ls_scf_line_search_3point,&
29 ls_scf_line_search_3point_2d
30 USE iterate_matrix, ONLY: purify_mcweeny
31 USE kinds, ONLY: dp
32 USE machine, ONLY: m_flush
33 USE mathconstants, ONLY: ifac
34 USE mathlib, ONLY: invmat
35#include "./base/base_uses.f90"
36
37 IMPLICIT NONE
38
39 PRIVATE
40
41 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'dm_ls_scf_curvy'
42
43 PUBLIC :: dm_ls_curvy_optimization, deallocate_curvy_data
44
45CONTAINS
46
47! **************************************************************************************************
48!> \brief driver routine for Head-Gordon curvy step approach
49!> \param ls_scf_env ...
50!> \param energy ...
51!> \param check_conv ...
52!> \par History
53!> 2012.05 created [Florian Schiffmann]
54!> \author Florian Schiffmann
55! **************************************************************************************************
56
57 SUBROUTINE dm_ls_curvy_optimization(ls_scf_env, energy, check_conv)
58 TYPE(ls_scf_env_type) :: ls_scf_env
59 REAL(kind=dp) :: energy
60 LOGICAL :: check_conv
61
62 CHARACTER(LEN=*), PARAMETER :: routinen = 'dm_ls_curvy_optimization'
63
64 INTEGER :: handle, i, lsstep
65
66 CALL timeset(routinen, handle)
67
68 CALL cite_reference(shao2003)
69
70! Upon first call initialize all matrices needed curing optimization
71! In addition transform P into orthonormal basis. Will be scaled by 0.5 in closed shell case
72! Only to be done once as it will be stored and reused afterwards
73! TRS4 might yield a non-idempotent P therefore McWeeny purification is applied on initial P
74
75 IF (.NOT. ALLOCATED(ls_scf_env%curvy_data%matrix_dp)) THEN
76 CALL init_curvy(ls_scf_env%curvy_data, ls_scf_env%matrix_s, ls_scf_env%nspins)
77 ls_scf_env%curvy_data%line_search_step = 1
78
79 IF (ls_scf_env%curvy_data%line_search_type == ls_scf_line_search_3point_2d) THEN
80 DO i = 1, ls_scf_env%nspins
81 CALL dbcsr_copy(ls_scf_env%curvy_data%matrix_psave(i, 1), &
82 ls_scf_env%matrix_p(i))
83 END DO
84 END IF
85 IF (ls_scf_env%nspins == 1) CALL dbcsr_scale(ls_scf_env%matrix_p(1), 0.5_dp)
86 CALL transform_matrix_orth(ls_scf_env%matrix_p, ls_scf_env%matrix_s_sqrt, &
87 ls_scf_env%eps_filter)
88 CALL purify_mcweeny(ls_scf_env%matrix_p, ls_scf_env%eps_filter, 3)
89 DO i = 1, ls_scf_env%nspins
90 CALL dbcsr_copy(ls_scf_env%curvy_data%matrix_p(i), ls_scf_env%matrix_p(i))
91 END DO
92 END IF
93
94 lsstep = ls_scf_env%curvy_data%line_search_step
95
96! If new search direction has to be computed transform H into the orthnormal basis
97
98 IF (ls_scf_env%curvy_data%line_search_step == 1) &
99 CALL transform_matrix_orth(ls_scf_env%matrix_ks, ls_scf_env%matrix_s_sqrt_inv, &
100 ls_scf_env%eps_filter)
101
102! Set the energies for the line search and make sure to give the correct energy back to scf_main
103 ls_scf_env%curvy_data%energies(lsstep) = energy
104 IF (lsstep .NE. 1) energy = ls_scf_env%curvy_data%energies(1)
105
106! start the optimization by calling the driver routine or simply combine saved P(2D line search)
107 IF (lsstep .LE. 2) THEN
108 CALL optimization_step(ls_scf_env%curvy_data, ls_scf_env)
109 ELSE IF (lsstep == ls_scf_env%curvy_data%line_search_type) THEN
110! line_search type has the value appropriate to the number of energy calculations needed
111 CALL optimization_step(ls_scf_env%curvy_data, ls_scf_env)
112 ELSE
113 CALL new_p_from_save(ls_scf_env%matrix_p, ls_scf_env%curvy_data%matrix_psave, lsstep, &
114 ls_scf_env%curvy_data%double_step_size)
115 ls_scf_env%curvy_data%line_search_step = ls_scf_env%curvy_data%line_search_step + 1
116 CALL timestop(handle)
117 RETURN
118 END IF
119 lsstep = ls_scf_env%curvy_data%line_search_step
120
121! transform new density matrix back into nonorthonormal basis (again scaling might apply)
122
123 CALL transform_matrix_orth(ls_scf_env%matrix_p, ls_scf_env%matrix_s_sqrt_inv, &
124 ls_scf_env%eps_filter)
125 IF (ls_scf_env%nspins == 1) CALL dbcsr_scale(ls_scf_env%matrix_p(1), 2.0_dp)
126
127! P-matrices only need to be stored in case of 2D line search
128 IF (lsstep .LE. 3 .AND. ls_scf_env%curvy_data%line_search_type == ls_scf_line_search_3point_2d) THEN
129 DO i = 1, ls_scf_env%nspins
130 CALL dbcsr_copy(ls_scf_env%curvy_data%matrix_psave(i, lsstep), &
131 ls_scf_env%matrix_p(i))
132 END DO
133 END IF
134 check_conv = lsstep == 1
135
136 CALL timestop(handle)
137
138 END SUBROUTINE dm_ls_curvy_optimization
139
140! **************************************************************************************************
141!> \brief low level routine for Head-Gordons curvy step approach
142!> computes gradients, performs a cg and line search,
143!> and evaluates the BCH series to obtain the new P matrix
144!> \param curvy_data ...
145!> \param ls_scf_env ...
146!> \par History
147!> 2012.05 created [Florian Schiffmann]
148!> \author Florian Schiffmann
149! **************************************************************************************************
150
151 SUBROUTINE optimization_step(curvy_data, ls_scf_env)
152 TYPE(ls_scf_curvy_type) :: curvy_data
153 TYPE(ls_scf_env_type) :: ls_scf_env
154
155 CHARACTER(LEN=*), PARAMETER :: routinen = 'optimization_step'
156
157 INTEGER :: handle, ispin
158 REAL(kind=dp) :: filter, step_size(2)
159
160! Upon first line search step compute new search direction and apply CG if required
161
162 CALL timeset(routinen, handle)
163
164 IF (curvy_data%line_search_step == 1) THEN
165 curvy_data%step_size = maxval(curvy_data%step_size)
166 curvy_data%step_size = min(max(0.10_dp, 0.5_dp*abs(curvy_data%step_size(1))), 0.5_dp)
167! Dynamic eps_filter for newton steps
168 filter = max(ls_scf_env%eps_filter*curvy_data%min_filter, &
169 ls_scf_env%eps_filter*curvy_data%filter_factor)
170 CALL compute_direction_newton(curvy_data%matrix_p, ls_scf_env%matrix_ks, &
171 curvy_data%matrix_dp, filter, curvy_data%fix_shift, curvy_data%shift, &
172 curvy_data%cg_numer, curvy_data%cg_denom, curvy_data%min_shift)
173 curvy_data%filter_factor = curvy_data%scale_filter*curvy_data%filter_factor
174 step_size = curvy_data%step_size
175 curvy_data%BCH_saved = 0
176 ELSE IF (curvy_data%line_search_step == 2) THEN
177 step_size = curvy_data%step_size
178 IF (curvy_data%energies(1) - curvy_data%energies(2) .GT. 0.0_dp) THEN
179 curvy_data%step_size = curvy_data%step_size*2.0_dp
180 curvy_data%double_step_size = .true.
181 ELSE
182 curvy_data%step_size = curvy_data%step_size*0.5_dp
183 curvy_data%double_step_size = .false.
184 END IF
185 step_size = curvy_data%step_size
186 ELSE IF (curvy_data%line_search_step == ls_scf_line_search_3point_2d) THEN
187 CALL line_search_2d(curvy_data%energies, curvy_data%step_size)
188 step_size = curvy_data%step_size
189 ELSE IF (curvy_data%line_search_step == ls_scf_line_search_3point) THEN
190 CALL line_search_3pnt(curvy_data%energies, curvy_data%step_size)
191 step_size = curvy_data%step_size
192 END IF
193
194 CALL update_p_exp(curvy_data%matrix_p, ls_scf_env%matrix_p, curvy_data%matrix_dp, &
195 curvy_data%matrix_BCH, ls_scf_env%eps_filter, step_size, curvy_data%BCH_saved, &
196 curvy_data%n_bch_hist)
197
198! line_search type has the value appropriate to the numeber of energy calculations needed
199 curvy_data%line_search_step = mod(curvy_data%line_search_step, curvy_data%line_search_type) + 1
200 IF (curvy_data%line_search_step == 1) THEN
201 DO ispin = 1, SIZE(curvy_data%matrix_p)
202 CALL dbcsr_copy(curvy_data%matrix_p(ispin), ls_scf_env%matrix_p(ispin))
203 END DO
204 END IF
205 CALL timestop(handle)
206
207 END SUBROUTINE optimization_step
208
209! **************************************************************************************************
210!> \brief Perform a 6pnt-2D line search for spin polarized calculations.
211!> Fit a 2D parabolic function to 6 points
212!> \param energies ...
213!> \param step_size ...
214!> \par History
215!> 2012.05 created [Florian Schiffmann]
216!> \author Florian Schiffmann
217! **************************************************************************************************
218
219 SUBROUTINE line_search_2d(energies, step_size)
220 REAL(kind=dp) :: energies(6), step_size(2)
221
222 INTEGER :: info, unit_nr
223 REAL(kind=dp) :: e_pred, param(6), s1, s1sq, s2, s2sq, &
224 sys_lin_eq(6, 6), tmp_e, v1, v2
225 TYPE(cp_logger_type), POINTER :: logger
226
227 logger => cp_get_default_logger()
228 IF (energies(1) - energies(2) .LT. 0._dp) THEN
229 tmp_e = energies(2); energies(2) = energies(3); energies(3) = tmp_e
230 step_size = step_size*2.0_dp
231 END IF
232 IF (logger%para_env%is_source()) THEN
233 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
234 ELSE
235 unit_nr = -1
236 END IF
237 s1 = 0.5_dp*step_size(1); s2 = step_size(1); s1sq = s1**2; s2sq = s2**2
238 sys_lin_eq = 0.0_dp; sys_lin_eq(:, 6) = 1.0_dp
239 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
240 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
241 sys_lin_eq(4, 3) = s1sq; sys_lin_eq(4, 5) = s1
242 sys_lin_eq(5, 1) = s1sq; sys_lin_eq(5, 4) = s1
243 sys_lin_eq(6, 3) = s2sq; sys_lin_eq(6, 5) = s2
244
245 CALL invmat(sys_lin_eq, info)
246 param = matmul(sys_lin_eq, energies)
247 v1 = (param(2)*param(4))/(2.0_dp*param(1)) - param(5)
248 v2 = -(param(2)**2)/(2.0_dp*param(1)) + 2.0_dp*param(3)
249 step_size(2) = v1/v2
250 step_size(1) = (-param(2)*step_size(2) - param(4))/(2.0_dp*param(1))
251 IF (step_size(1) .LT. 0.0_dp) step_size(1) = 1.0_dp
252 IF (step_size(2) .LT. 0.0_dp) step_size(2) = 1.0_dp
253! step_size(1)=MIN(step_size(1),2.0_dp)
254! step_size(2)=MIN(step_size(2),2.0_dp)
255 e_pred = param(1)*step_size(1)**2 + param(2)*step_size(1)*step_size(2) + &
256 param(3)*step_size(2)**2 + param(4)*step_size(1) + param(5)*step_size(2) + param(6)
257 IF (unit_nr .GT. 0) WRITE (unit_nr, "(t3,a,F10.5,F10.5,A,F20.9)") &
258 " Line Search: Step Size", step_size, " Predicted energy", e_pred
259 e_pred = param(1)*s1**2 + param(2)*s2*s1*0.0_dp + &
260 param(3)*s1**2*0.0_dp + param(4)*s1 + param(5)*s1*0.0_dp + param(6)
261
262 END SUBROUTINE line_search_2d
263
264! **************************************************************************************************
265!> \brief Perform a 3pnt line search
266!> \param energies ...
267!> \param step_size ...
268!> \par History
269!> 2012.05 created [Florian Schiffmann]
270!> \author Florian Schiffmann
271! **************************************************************************************************
272
273 SUBROUTINE line_search_3pnt(energies, step_size)
274 REAL(kind=dp) :: energies(3), step_size(2)
275
276 INTEGER :: unit_nr
277 REAL(kind=dp) :: a, b, c, e_pred, min_val, step1, tmp, &
278 tmp_e
279 TYPE(cp_logger_type), POINTER :: logger
280
281 logger => cp_get_default_logger()
282 IF (energies(1) - energies(2) .LT. 0._dp) THEN
283 tmp_e = energies(2); energies(2) = energies(3); energies(3) = tmp_e
284 step_size = step_size*2.0_dp
285 END IF
286 IF (logger%para_env%is_source()) THEN
287 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
288 ELSE
289 unit_nr = -1
290 END IF
291 step1 = 0.5_dp*step_size(1)
292 c = energies(1)
293 a = (energies(3) + c - 2.0_dp*energies(2))/(2.0_dp*step1**2)
294 b = (energies(2) - c - a*step1**2)/step1
295 IF (a .LT. 1.0e-12_dp) a = -1.0e-12_dp
296 min_val = -b/(2.0_dp*a)
297 e_pred = a*min_val**2 + b*min_val + c
298 tmp = step_size(1)
299 IF (e_pred .LT. energies(1) .AND. e_pred .LT. energies(2)) THEN
300 step_size = max(-1.0_dp, &
301 min(min_val, 10_dp*step_size))
302 ELSE
303 step_size = 1.0_dp
304 END IF
305 e_pred = a*(step_size(1))**2 + b*(step_size(1)) + c
306 IF (unit_nr .GT. 0) THEN
307 WRITE (unit_nr, "(t3,a,f16.8,a,F20.9)") "Line Search: Step Size", step_size(1), " Predicted energy", e_pred
308 CALL m_flush(unit_nr)
309 END IF
310 END SUBROUTINE line_search_3pnt
311
312! **************************************************************************************************
313!> \brief Get a new search direction. Iterate to obtain a Newton like step
314!> Refine with a CG update of the search direction
315!> \param matrix_p ...
316!> \param matrix_ks ...
317!> \param matrix_dp ...
318!> \param eps_filter ...
319!> \param fix_shift ...
320!> \param curvy_shift ...
321!> \param cg_numer ...
322!> \param cg_denom ...
323!> \param min_shift ...
324!> \par History
325!> 2012.05 created [Florian Schiffmann]
326!> \author Florian Schiffmann
327! **************************************************************************************************
328
329 SUBROUTINE compute_direction_newton(matrix_p, matrix_ks, matrix_dp, eps_filter, fix_shift, &
330 curvy_shift, cg_numer, cg_denom, min_shift)
331 TYPE(dbcsr_type), DIMENSION(:) :: matrix_p, matrix_ks, matrix_dp
332 REAL(kind=dp) :: eps_filter
333 LOGICAL :: fix_shift(2)
334 REAL(kind=dp) :: curvy_shift(2), cg_numer(2), &
335 cg_denom(2), min_shift
336
337 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_direction_newton'
338
339 INTEGER :: handle, i, ispin, ncyc, nspin, unit_nr
340 LOGICAL :: at_limit
341 REAL(kind=dp) :: beta, conv_val, maxel, old_conv, shift
342 TYPE(cp_logger_type), POINTER :: logger
343 TYPE(dbcsr_type) :: matrix_ax, matrix_b, matrix_cg, &
344 matrix_dp_old, matrix_pks, matrix_res, &
345 matrix_tmp, matrix_tmp1
346
347 logger => cp_get_default_logger()
348
349 IF (logger%para_env%is_source()) THEN
350 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
351 ELSE
352 unit_nr = -1
353 END IF
354 CALL timeset(routinen, handle)
355 nspin = SIZE(matrix_p)
356
357 CALL dbcsr_create(matrix_pks, template=matrix_dp(1), matrix_type=dbcsr_type_no_symmetry)
358 CALL dbcsr_create(matrix_ax, template=matrix_dp(1), matrix_type=dbcsr_type_no_symmetry)
359 CALL dbcsr_create(matrix_tmp, template=matrix_dp(1), matrix_type=dbcsr_type_no_symmetry)
360 CALL dbcsr_create(matrix_tmp1, template=matrix_dp(1), matrix_type=dbcsr_type_no_symmetry)
361 CALL dbcsr_create(matrix_res, template=matrix_dp(1), matrix_type=dbcsr_type_no_symmetry)
362 CALL dbcsr_create(matrix_cg, template=matrix_dp(1), matrix_type=dbcsr_type_no_symmetry)
363 CALL dbcsr_create(matrix_b, template=matrix_dp(1), matrix_type=dbcsr_type_no_symmetry)
364 CALL dbcsr_create(matrix_dp_old, template=matrix_dp(1), matrix_type=dbcsr_type_no_symmetry)
365
366 DO ispin = 1, nspin
367 CALL dbcsr_copy(matrix_dp_old, matrix_dp(ispin))
368
369! Precompute some matrices to save work during iterations
370 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_p(ispin), matrix_ks(ispin), &
371 0.0_dp, matrix_pks, filter_eps=eps_filter)
372 CALL dbcsr_transposed(matrix_b, matrix_pks)
373 CALL dbcsr_copy(matrix_cg, matrix_b)
374
375! Starting CG with guess 0-matrix gives -2*gradient=[Ks*P-(Ks*P)T] for cg_matrix in second step
376 CALL dbcsr_add(matrix_cg, matrix_pks, 2.0_dp, -2.0_dp)
377
378! Residual matrix in first step=cg matrix. Keep Pks for later use in CG!
379 CALL dbcsr_copy(matrix_res, matrix_cg)
380
381! Precompute -FP-[FP]T which will be used throughout the CG iterations
382 CALL dbcsr_add(matrix_b, matrix_pks, -1.0_dp, -1.0_dp)
383
384! Setup some values to check convergence and safety checks for eigenvalue shifting
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
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:192
iterate_matrix::purify_mcweeny
Definition iterate_matrix.F:58
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition bibliography.F:21
bibliography::shao2003
integer, save, public shao2003
Definition bibliography.F:36
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_transposed
subroutine, public dbcsr_transposed(transposed, normal, shallow_data_copy, transpose_distribution, use_distribution)
...
Definition cp_dbcsr_api.F:1219
cp_dbcsr_api::dbcsr_scale
subroutine, public dbcsr_scale(matrix, alpha_scalar)
...
Definition cp_dbcsr_api.F:1165
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:368
cp_dbcsr_api::dbcsr_multiply
subroutine, public dbcsr_multiply(transa, transb, alpha, matrix_a, matrix_b, beta, matrix_c, first_row, last_row, first_column, last_column, first_k, last_k, retain_sparsity, filter_eps, flop)
...
Definition cp_dbcsr_api.F:1073
cp_dbcsr_api::dbcsr_filter
subroutine, public dbcsr_filter(matrix, eps)
...
Definition cp_dbcsr_api.F:644
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1181
cp_dbcsr_api::dbcsr_release
subroutine, public dbcsr_release(matrix)
...
Definition cp_dbcsr_api.F:1119
cp_dbcsr_api::dbcsr_add
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
Definition cp_dbcsr_api.F:251
cp_dbcsr_contrib
Definition cp_dbcsr_contrib.F:8
cp_dbcsr_contrib::dbcsr_frobenius_norm
real(dp) function, public dbcsr_frobenius_norm(matrix)
Compute the frobenius norm of a dbcsr matrix.
Definition cp_dbcsr_contrib.F:127
cp_dbcsr_contrib::dbcsr_dot
subroutine, public dbcsr_dot(matrix_a, matrix_b, trace)
Computes the dot product of two matrices, also known as the trace of their matrix product.
Definition cp_dbcsr_contrib.F:367
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:58
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:981
input_constants::ls_scf_line_search_3point_2d
integer, parameter, public ls_scf_line_search_3point_2d
Definition input_constants.F:981
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:130
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_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:174
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