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