Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief density matrix optimization using exponential transformations
10 : !> \par History
11 : !> 2012.05 created [Florian Schiffmann]
12 : !> \author Florian Schiffmann
13 : ! **************************************************************************************************
14 :
15 : MODULE 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 :
45 : CONTAINS
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 90 : 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 90 : CALL timeset(routineN, handle)
67 :
68 90 : 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 90 : IF (.NOT. ALLOCATED(ls_scf_env%curvy_data%matrix_dp)) THEN
76 18 : CALL init_curvy(ls_scf_env%curvy_data, ls_scf_env%matrix_s, ls_scf_env%nspins)
77 18 : ls_scf_env%curvy_data%line_search_step = 1
78 :
79 18 : IF (ls_scf_env%curvy_data%line_search_type == ls_scf_line_search_3point_2d) THEN
80 6 : DO i = 1, ls_scf_env%nspins
81 : CALL dbcsr_copy(ls_scf_env%curvy_data%matrix_psave(i, 1), &
82 6 : ls_scf_env%matrix_p(i))
83 : END DO
84 : END IF
85 18 : 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 18 : ls_scf_env%eps_filter)
88 18 : CALL purify_mcweeny(ls_scf_env%matrix_p, ls_scf_env%eps_filter, 3)
89 38 : DO i = 1, ls_scf_env%nspins
90 38 : 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 90 : 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 90 : 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 28 : 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 90 : ls_scf_env%curvy_data%energies(lsstep) = energy
105 90 : 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 90 : IF (lsstep <= 2) THEN
109 56 : CALL optimization_step(ls_scf_env%curvy_data, ls_scf_env)
110 34 : 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 28 : 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 6 : ls_scf_env%curvy_data%double_step_size)
116 6 : ls_scf_env%curvy_data%line_search_step = ls_scf_env%curvy_data%line_search_step + 1
117 6 : CALL timestop(handle)
118 6 : RETURN
119 : END IF
120 84 : 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 84 : ls_scf_env%eps_filter)
126 84 : 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 84 : IF (lsstep <= 3 .AND. ls_scf_env%curvy_data%line_search_type == ls_scf_line_search_3point_2d) THEN
130 18 : DO i = 1, ls_scf_env%nspins
131 : CALL dbcsr_copy(ls_scf_env%curvy_data%matrix_psave(i, lsstep), &
132 18 : ls_scf_env%matrix_p(i))
133 : END DO
134 : END IF
135 84 : check_conv = lsstep == 1
136 :
137 84 : 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 84 : 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 84 : CALL timeset(routineN, handle)
164 :
165 84 : IF (curvy_data%line_search_step == 1) THEN
166 168 : curvy_data%step_size = MAXVAL(curvy_data%step_size)
167 84 : 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 28 : 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 28 : curvy_data%cg_numer, curvy_data%cg_denom, curvy_data%min_shift)
174 28 : curvy_data%filter_factor = curvy_data%scale_filter*curvy_data%filter_factor
175 84 : step_size = curvy_data%step_size
176 84 : curvy_data%BCH_saved = 0
177 56 : ELSE IF (curvy_data%line_search_step == 2) THEN
178 84 : step_size = curvy_data%step_size
179 28 : IF (curvy_data%energies(1) - curvy_data%energies(2) > 0.0_dp) THEN
180 72 : curvy_data%step_size = curvy_data%step_size*2.0_dp
181 24 : curvy_data%double_step_size = .TRUE.
182 : ELSE
183 12 : curvy_data%step_size = curvy_data%step_size*0.5_dp
184 4 : curvy_data%double_step_size = .FALSE.
185 : END IF
186 84 : step_size = curvy_data%step_size
187 28 : ELSE IF (curvy_data%line_search_step == ls_scf_line_search_3point_2d) THEN
188 2 : CALL line_search_2d(curvy_data%energies, curvy_data%step_size)
189 6 : step_size = curvy_data%step_size
190 26 : ELSE IF (curvy_data%line_search_step == ls_scf_line_search_3point) THEN
191 26 : CALL line_search_3pnt(curvy_data%energies, curvy_data%step_size)
192 78 : 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 84 : curvy_data%n_bch_hist)
198 :
199 : ! line_search type has the value appropriate to the numeber of energy calculations needed
200 84 : curvy_data%line_search_step = MOD(curvy_data%line_search_step, curvy_data%line_search_type) + 1
201 84 : IF (curvy_data%line_search_step == 1) THEN
202 58 : DO ispin = 1, SIZE(curvy_data%matrix_p)
203 58 : CALL dbcsr_copy(curvy_data%matrix_p(ispin), ls_scf_env%matrix_p(ispin))
204 : END DO
205 : END IF
206 84 : CALL timestop(handle)
207 :
208 84 : 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 2 : 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 2 : logger => cp_get_default_logger()
229 2 : IF (energies(1) - energies(2) < 0._dp) THEN
230 0 : tmp_e = energies(2); energies(2) = energies(3); energies(3) = tmp_e
231 0 : step_size = step_size*2.0_dp
232 : END IF
233 2 : IF (logger%para_env%is_source()) THEN
234 1 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
235 : ELSE
236 : unit_nr = -1
237 : END IF
238 2 : s1 = 0.5_dp*step_size(1); s2 = step_size(1); s1sq = s1**2; s2sq = s2**2
239 14 : sys_lin_eq = 0.0_dp; sys_lin_eq(:, 6) = 1.0_dp
240 2 : 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 2 : 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 2 : sys_lin_eq(4, 3) = s1sq; sys_lin_eq(4, 5) = s1
243 2 : sys_lin_eq(5, 1) = s1sq; sys_lin_eq(5, 4) = s1
244 2 : sys_lin_eq(6, 3) = s2sq; sys_lin_eq(6, 5) = s2
245 :
246 2 : CALL invmat(sys_lin_eq, info)
247 86 : param = MATMUL(sys_lin_eq, energies)
248 2 : v1 = (param(2)*param(4))/(2.0_dp*param(1)) - param(5)
249 2 : v2 = -(param(2)**2)/(2.0_dp*param(1)) + 2.0_dp*param(3)
250 2 : step_size(2) = v1/v2
251 2 : step_size(1) = (-param(2)*step_size(2) - param(4))/(2.0_dp*param(1))
252 2 : IF (step_size(1) < 0.0_dp) step_size(1) = 1.0_dp
253 2 : 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 2 : param(3)*step_size(2)**2 + param(4)*step_size(1) + param(5)*step_size(2) + param(6)
258 2 : IF (unit_nr > 0) WRITE (unit_nr, "(t3,a,F10.5,F10.5,A,F20.9)") &
259 1 : " 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 2 : 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 26 : 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 26 : logger => cp_get_default_logger()
283 26 : IF (energies(1) - energies(2) < 0._dp) THEN
284 2 : tmp_e = energies(2); energies(2) = energies(3); energies(3) = tmp_e
285 6 : step_size = step_size*2.0_dp
286 : END IF
287 26 : IF (logger%para_env%is_source()) THEN
288 13 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
289 : ELSE
290 13 : unit_nr = -1
291 : END IF
292 26 : step1 = 0.5_dp*step_size(1)
293 26 : c = energies(1)
294 26 : a = (energies(3) + c - 2.0_dp*energies(2))/(2.0_dp*step1**2)
295 26 : b = (energies(2) - c - a*step1**2)/step1
296 26 : IF (a < 1.0E-12_dp) a = -1.0E-12_dp
297 26 : min_val = -b/(2.0_dp*a)
298 26 : e_pred = a*min_val**2 + b*min_val + c
299 26 : tmp = step_size(1)
300 26 : IF (e_pred < energies(1) .AND. e_pred < energies(2)) THEN
301 : step_size = MAX(-1.0_dp, &
302 54 : MIN(min_val, 10_dp*step_size))
303 : ELSE
304 24 : step_size = 1.0_dp
305 : END IF
306 26 : e_pred = a*(step_size(1))**2 + b*(step_size(1)) + c
307 26 : IF (unit_nr > 0) THEN
308 13 : WRITE (unit_nr, "(t3,a,f16.8,a,F20.9)") "Line Search: Step Size", step_size(1), " Predicted energy", e_pred
309 13 : CALL m_flush(unit_nr)
310 : END IF
311 26 : 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 28 : 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 56 : logger => cp_get_default_logger()
349 :
350 28 : IF (logger%para_env%is_source()) THEN
351 14 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
352 : ELSE
353 14 : unit_nr = -1
354 : END IF
355 28 : CALL timeset(routineN, handle)
356 28 : nspin = SIZE(matrix_p)
357 :
358 28 : CALL dbcsr_create(matrix_PKs, template=matrix_dp(1), matrix_type=dbcsr_type_no_symmetry)
359 28 : CALL dbcsr_create(matrix_Ax, template=matrix_dp(1), matrix_type=dbcsr_type_no_symmetry)
360 28 : CALL dbcsr_create(matrix_tmp, template=matrix_dp(1), matrix_type=dbcsr_type_no_symmetry)
361 28 : CALL dbcsr_create(matrix_tmp1, template=matrix_dp(1), matrix_type=dbcsr_type_no_symmetry)
362 28 : CALL dbcsr_create(matrix_res, template=matrix_dp(1), matrix_type=dbcsr_type_no_symmetry)
363 28 : CALL dbcsr_create(matrix_cg, template=matrix_dp(1), matrix_type=dbcsr_type_no_symmetry)
364 28 : CALL dbcsr_create(matrix_b, template=matrix_dp(1), matrix_type=dbcsr_type_no_symmetry)
365 28 : CALL dbcsr_create(matrix_dp_old, template=matrix_dp(1), matrix_type=dbcsr_type_no_symmetry)
366 :
367 58 : DO ispin = 1, nspin
368 30 : 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 30 : 0.0_dp, matrix_PKs, filter_eps=eps_filter)
373 30 : CALL dbcsr_transposed(matrix_b, matrix_PKs)
374 30 : 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 30 : 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 30 : CALL dbcsr_copy(matrix_res, matrix_cg)
381 :
382 : ! Precompute -FP-[FP]T which will be used throughout the CG iterations
383 30 : 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 30 : old_conv = dbcsr_frobenius_norm(matrix_res)
387 30 : shift = MIN(10.0_dp, MAX(min_shift, 0.05_dp*old_conv))
388 30 : conv_val = MAX(0.010_dp*old_conv, 100.0_dp*eps_filter)
389 30 : old_conv = 100.0_dp
390 30 : IF (fix_shift(ispin)) THEN
391 0 : shift = MAX(min_shift, MIN(10.0_dp, MAX(shift, curvy_shift(ispin) - 0.5_dp*curvy_shift(ispin))))
392 0 : curvy_shift(ispin) = shift
393 : END IF
394 :
395 : ! Begin the real optimization loop
396 30 : CALL dbcsr_set(matrix_dp(ispin), 0.0_dp)
397 30 : ncyc = 10
398 104 : DO i = 1, ncyc
399 :
400 : ! One step to compute: -FPD-DPF-DFP-PFD (not obvious but symmetry allows for some tricks)
401 104 : 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 104 : 0.0_dp, matrix_tmp, filter_eps=eps_filter)
406 104 : CALL commutator_symm(matrix_ks(ispin), matrix_tmp, matrix_tmp1, eps_filter, 2.0_dp)
407 104 : 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 104 : 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 104 : matrix_tmp, eps_filter, at_limit)
414 104 : CALL dbcsr_filter(matrix_cg, eps_filter)
415 :
416 : ! check for convergence of the newton step
417 104 : maxel = dbcsr_frobenius_norm(matrix_res)
418 104 : IF (unit_nr > 0) THEN
419 52 : WRITE (unit_nr, "(T3,A,F12.6)") "Convergence of Newton iteration ", maxel
420 52 : CALL m_flush(unit_nr)
421 : END IF
422 104 : at_limit = at_limit .OR. (old_conv/maxel < 1.01_dp)
423 104 : old_conv = maxel
424 104 : IF (i == ncyc .AND. maxel/conv_val > 5.0_dp) THEN
425 0 : fix_shift(ispin) = .TRUE.
426 0 : curvy_shift(ispin) = 4.0_dp*shift
427 : END IF
428 104 : 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 30 : CALL dbcsr_transposed(matrix_b, matrix_PKs)
433 : !compute b= -2*KsP+2*PKs=-(2*gradient)
434 30 : CALL dbcsr_copy(matrix_cg, matrix_b)
435 30 : CALL dbcsr_add(matrix_cg, matrix_PKs, 1.0_dp, -1.0_dp)
436 30 : cg_denom(ispin) = cg_numer(ispin)
437 30 : CALL dbcsr_dot(matrix_cg, matrix_dp(ispin), cg_numer(ispin))
438 30 : beta = cg_numer(ispin)/MAX(cg_denom(ispin), 1.0E-6_dp)
439 30 : IF (beta < 1.0_dp) THEN
440 28 : beta = MAX(0.0_dp, beta)
441 28 : CALL dbcsr_add(matrix_dp(ispin), matrix_dp_old, 1.0_dp, beta)
442 : END IF
443 58 : IF (unit_nr > 0) WRITE (unit_nr, "(A)") " "
444 : END DO
445 :
446 28 : CALL dbcsr_release(matrix_PKs)
447 28 : CALL dbcsr_release(matrix_dp_old)
448 28 : CALL dbcsr_release(matrix_b)
449 28 : CALL dbcsr_release(matrix_Ax)
450 28 : CALL dbcsr_release(matrix_tmp)
451 28 : CALL dbcsr_release(matrix_tmp1)
452 28 : CALL dbcsr_release(matrix_b)
453 28 : CALL dbcsr_release(matrix_res)
454 28 : CALL dbcsr_release(matrix_cg)
455 :
456 28 : IF (unit_nr > 0) CALL m_flush(unit_nr)
457 28 : CALL timestop(handle)
458 28 : 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 104 : 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 104 : at_limit = .FALSE.
486 104 : CALL dbcsr_dot(res, res, norm_rr)
487 104 : CALL dbcsr_dot(cg, Ax, norm_cA)
488 104 : lin_eq = 0.0_dp
489 104 : fac = norm_rr/norm_cA
490 104 : fac1 = fac
491 : ! Use a 3point line search and a fit to a quadratic function to determine optimal step size
492 416 : DO i = 1, 3
493 312 : CALL dbcsr_copy(tmp, res)
494 312 : CALL dbcsr_add(tmp, Ax, 1.0_dp, -fac)
495 312 : devi(i) = dbcsr_frobenius_norm(tmp)
496 1248 : lin_eq(i, :) = [fac**2, fac, 1.0_dp]
497 416 : fac = fac1 + fac1*((-1)**i)*0.5_dp
498 : END DO
499 104 : CALL invmat(lin_eq, info)
500 1352 : vec = MATMUL(lin_eq, devi)
501 104 : alpha = -vec(2)/(2.0_dp*vec(1))
502 104 : fac = SQRT(norm_rr/(norm_cA*alpha))
503 : !scale the previous matrices to match the step size
504 104 : CALL dbcsr_scale(Ax, fac)
505 104 : CALL dbcsr_scale(cg, fac)
506 104 : norm_cA = norm_cA*fac**2
507 :
508 : ! USe CG to get the new matrices
509 104 : alpha = norm_rr/norm_cA
510 104 : CALL dbcsr_add(res, Ax, 1.0_dp, -alpha)
511 104 : CALL dbcsr_dot(res, res, new_norm)
512 104 : IF (norm_rr < eps_filter*0.001_dp .OR. new_norm < eps_filter*0.001_dp) THEN
513 : beta = 0.0_dp
514 22 : at_limit = .TRUE.
515 : ELSE
516 : beta = new_norm/norm_rr
517 82 : CALL dbcsr_add(deltp, cg, 1.0_dp, alpha)
518 : END IF
519 104 : beta = new_norm/norm_rr
520 104 : CALL dbcsr_add(cg, res, beta, 1.0_dp)
521 :
522 208 : 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 6 : 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 8 : SELECT CASE (lsstep)
544 : CASE (3)
545 2 : CALL dbcsr_copy(matrix_p(1), matrix_psave(1, 1))
546 2 : IF (DOUBLE) THEN
547 2 : CALL dbcsr_copy(matrix_p(2), matrix_psave(2, 2))
548 : ELSE
549 0 : CALL dbcsr_copy(matrix_p(2), matrix_psave(2, 3))
550 : END IF
551 : CASE (4)
552 2 : IF (DOUBLE) THEN
553 2 : CALL dbcsr_copy(matrix_p(1), matrix_psave(1, 2))
554 : ELSE
555 0 : CALL dbcsr_copy(matrix_p(1), matrix_psave(1, 3))
556 : END IF
557 2 : CALL dbcsr_copy(matrix_p(2), matrix_psave(2, 1))
558 : CASE (5)
559 2 : CALL dbcsr_copy(matrix_p(1), matrix_psave(1, 1))
560 8 : IF (DOUBLE) THEN
561 2 : CALL dbcsr_copy(matrix_p(2), matrix_psave(2, 3))
562 : ELSE
563 0 : CALL dbcsr_copy(matrix_p(2), matrix_psave(2, 2))
564 : END IF
565 : END SELECT
566 :
567 6 : 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 208 : 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 208 : CALL timeset(routineN, handle)
591 :
592 208 : CALL dbcsr_create(work, template=a, matrix_type=dbcsr_type_no_symmetry)
593 :
594 208 : CALL dbcsr_multiply("N", "N", prefac, a, b, 0.0_dp, res, filter_eps=eps_filter)
595 208 : CALL dbcsr_transposed(work, res)
596 208 : CALL dbcsr_add(res, work, 1.0_dp, -1.0_dp)
597 :
598 208 : CALL dbcsr_release(work)
599 :
600 208 : CALL timestop(handle)
601 208 : 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 84 : 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 84 : CALL timeset(routineN, handle)
636 :
637 84 : logger => cp_get_default_logger()
638 84 : IF (logger%para_env%is_source()) THEN
639 42 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
640 : ELSE
641 42 : unit_nr = -1
642 : END IF
643 :
644 84 : CALL dbcsr_create(matrix, template=matrix_p_in(1), matrix_type=dbcsr_type_no_symmetry)
645 84 : CALL dbcsr_create(matrix_tmp, template=matrix_p_in(1), matrix_type=dbcsr_type_no_symmetry)
646 84 : nspin = SIZE(matrix_p_in)
647 :
648 174 : DO ispin = 1, nspin
649 90 : step_fac = 1.0_dp
650 90 : frob_norm = 1.0_dp
651 90 : nsave = 0
652 :
653 90 : CALL dbcsr_copy(matrix_tmp, matrix_p_in(ispin))
654 90 : 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 130 : DO i = 1, BCH_saved(ispin)
658 86 : step_fac = step_fac*step_size(ispin)
659 86 : CALL dbcsr_copy(matrix_tmp, matrix_p_out(ispin))
660 86 : CALL dbcsr_add(matrix_p_out(ispin), matrix_BCH(ispin, i), 1.0_dp, ifac(i)*step_fac)
661 86 : CALL dbcsr_add(matrix_tmp, matrix_p_out(ispin), 1.0_dp, -1.0_dp)
662 86 : frob_norm = dbcsr_frobenius_norm(matrix_tmp)
663 86 : IF (unit_nr > 0) WRITE (unit_nr, "(t3,a,i3,a,f16.8)") "BCH: step", i, " Norm of P_old-Pnew:", frob_norm
664 130 : IF (frob_norm < threshold) EXIT
665 : END DO
666 90 : 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 44 : save_BCH = BCH_saved(ispin) == 0 .AND. n_bch_hist > 0
670 86 : DO i = BCH_saved(ispin) + 1, 20
671 86 : 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 86 : 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 86 : CALL dbcsr_transposed(matrix_tmp, matrix)
680 86 : 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 86 : CALL dbcsr_copy(matrix_tmp, matrix_p_out(ispin))
684 86 : CALL dbcsr_add(matrix_p_out(ispin), matrix, 1.0_dp, ifac(i)*step_fac)
685 86 : IF (save_BCH .AND. i <= n_bch_hist) THEN
686 78 : CALL dbcsr_copy(matrix_BCH(ispin, i), matrix)
687 78 : nsave = i
688 : END IF
689 :
690 86 : 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 86 : frob_norm = dbcsr_frobenius_norm(matrix_tmp)
694 86 : IF (unit_nr > 0) WRITE (unit_nr, "(t3,a,i3,a,f16.8)") "BCH: step", i, " Norm of P_old-Pnew:", frob_norm
695 86 : 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 42 : CALL dbcsr_copy(matrix_tmp, matrix)
699 86 : CALL dbcsr_filter(matrix_tmp, threshold)
700 : END DO
701 44 : BCH_saved(ispin) = nsave
702 128 : IF (unit_nr > 0) WRITE (unit_nr, "(A)") " "
703 : END DO
704 :
705 84 : CALL purify_mcweeny(matrix_p_out, threshold, 1)
706 84 : IF (unit_nr > 0) CALL m_flush(unit_nr)
707 84 : CALL dbcsr_release(matrix_tmp)
708 84 : CALL dbcsr_release(matrix)
709 84 : CALL timestop(handle)
710 84 : 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 130 : 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 130 : CALL timeset(routineN, handle)
734 :
735 130 : CALL dbcsr_create(matrix_work, template=matrix(1), matrix_type=dbcsr_type_no_symmetry)
736 130 : CALL dbcsr_create(matrix_tmp, template=matrix(1), matrix_type=dbcsr_type_no_symmetry)
737 :
738 270 : DO ispin = 1, SIZE(matrix)
739 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix(ispin), matrix_trafo, &
740 140 : 0.0_dp, matrix_work, filter_eps=eps_filter)
741 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_trafo, matrix_work, &
742 140 : 0.0_dp, matrix_tmp, filter_eps=eps_filter)
743 : ! symmetrize results (this is again needed to make sure everything is stable)
744 140 : CALL dbcsr_transposed(matrix_work, matrix_tmp)
745 140 : CALL dbcsr_add(matrix_tmp, matrix_work, 0.5_dp, 0.5_dp)
746 270 : CALL dbcsr_copy(matrix(ispin), matrix_tmp)
747 : END DO
748 :
749 130 : CALL dbcsr_release(matrix_tmp)
750 130 : CALL dbcsr_release(matrix_work)
751 130 : CALL timestop(handle)
752 :
753 130 : END SUBROUTINE transform_matrix_orth
754 :
755 : ! **************************************************************************************************
756 : !> \brief ...
757 : !> \param curvy_data ...
758 : ! **************************************************************************************************
759 982 : SUBROUTINE deallocate_curvy_data(curvy_data)
760 : TYPE(ls_scf_curvy_type) :: curvy_data
761 :
762 : INTEGER :: i, j
763 :
764 982 : CALL release_dbcsr_array(curvy_data%matrix_dp)
765 982 : CALL release_dbcsr_array(curvy_data%matrix_p)
766 :
767 982 : IF (ALLOCATED(curvy_data%matrix_psave)) THEN
768 6 : DO i = 1, SIZE(curvy_data%matrix_psave, 1)
769 18 : DO j = 1, 3
770 16 : CALL dbcsr_release(curvy_data%matrix_psave(i, j))
771 : END DO
772 : END DO
773 2 : DEALLOCATE (curvy_data%matrix_psave)
774 : END IF
775 982 : IF (ALLOCATED(curvy_data%matrix_BCH)) THEN
776 38 : DO i = 1, SIZE(curvy_data%matrix_BCH, 1)
777 178 : DO j = 1, 7
778 160 : CALL dbcsr_release(curvy_data%matrix_BCH(i, j))
779 : END DO
780 : END DO
781 18 : DEALLOCATE (curvy_data%matrix_BCH)
782 : END IF
783 982 : END SUBROUTINE deallocate_curvy_data
784 :
785 : ! **************************************************************************************************
786 : !> \brief ...
787 : !> \param matrix ...
788 : ! **************************************************************************************************
789 1964 : SUBROUTINE release_dbcsr_array(matrix)
790 : TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:) :: matrix
791 :
792 : INTEGER :: i
793 :
794 1964 : IF (ALLOCATED(matrix)) THEN
795 76 : DO i = 1, SIZE(matrix)
796 76 : CALL dbcsr_release(matrix(i))
797 : END DO
798 36 : DEALLOCATE (matrix)
799 : END IF
800 1964 : END SUBROUTINE release_dbcsr_array
801 :
802 : ! **************************************************************************************************
803 : !> \brief ...
804 : !> \param curvy_data ...
805 : !> \param matrix_s ...
806 : !> \param nspins ...
807 : ! **************************************************************************************************
808 18 : 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 74 : ALLOCATE (curvy_data%matrix_dp(nspins))
816 56 : ALLOCATE (curvy_data%matrix_p(nspins))
817 38 : DO ispin = 1, nspins
818 : CALL dbcsr_create(curvy_data%matrix_dp(ispin), template=matrix_s, &
819 20 : matrix_type=dbcsr_type_no_symmetry)
820 20 : CALL dbcsr_set(curvy_data%matrix_dp(ispin), 0.0_dp)
821 : CALL dbcsr_create(curvy_data%matrix_p(ispin), template=matrix_s, &
822 20 : matrix_type=dbcsr_type_no_symmetry)
823 60 : curvy_data%fix_shift = .FALSE.
824 20 : curvy_data%double_step_size = .TRUE.
825 60 : curvy_data%shift = 1.0_dp
826 60 : curvy_data%BCH_saved = 0
827 60 : curvy_data%step_size = 0.60_dp
828 60 : curvy_data%cg_numer = 0.00_dp
829 78 : curvy_data%cg_denom = 0.00_dp
830 : END DO
831 18 : IF (curvy_data%line_search_type == ls_scf_line_search_3point_2d) THEN
832 24 : ALLOCATE (curvy_data%matrix_psave(nspins, 3))
833 6 : DO ispin = 1, nspins
834 18 : DO j = 1, 3
835 : CALL dbcsr_create(curvy_data%matrix_psave(ispin, j), template=matrix_s, &
836 16 : matrix_type=dbcsr_type_no_symmetry)
837 : END DO
838 : END DO
839 : END IF
840 18 : IF (curvy_data%n_bch_hist > 0) THEN
841 338 : ALLOCATE (curvy_data%matrix_BCH(nspins, curvy_data%n_bch_hist))
842 38 : DO ispin = 1, nspins
843 178 : DO j = 1, curvy_data%n_bch_hist
844 : CALL dbcsr_create(curvy_data%matrix_BCH(ispin, j), template=matrix_s, &
845 160 : matrix_type=dbcsr_type_no_symmetry)
846 : END DO
847 : END DO
848 : END IF
849 :
850 18 : END SUBROUTINE init_curvy
851 :
852 : END MODULE dm_ls_scf_curvy
|
