Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief density matrix optimization using exponential transformations
10 : !> \par History
11 : !> 2012.05 created [Florian Schiffmann]
12 : !> \author Florian Schiffmann
13 : ! **************************************************************************************************
14 :
15 : MODULE dm_ls_scf_curvy
16 : USE bibliography, ONLY: Shao2003,&
17 : cite_reference
18 : USE cp_log_handling, ONLY: cp_get_default_logger,&
19 : cp_logger_get_default_unit_nr,&
20 : cp_logger_type
21 : USE dbcsr_api, ONLY: &
22 : dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_dot, dbcsr_filter, dbcsr_frobenius_norm, &
23 : dbcsr_multiply, dbcsr_norm, dbcsr_release, dbcsr_scale, dbcsr_set, dbcsr_transposed, &
24 : dbcsr_type, dbcsr_type_no_symmetry
25 : USE dm_ls_scf_types, ONLY: ls_scf_curvy_type,&
26 : ls_scf_env_type
27 : USE input_constants, ONLY: ls_scf_line_search_3point,&
28 : ls_scf_line_search_3point_2d
29 : USE iterate_matrix, ONLY: purify_mcweeny
30 : USE kinds, ONLY: dp
31 : USE machine, ONLY: m_flush
32 : USE mathconstants, ONLY: ifac
33 : USE mathlib, ONLY: invmat
34 : #include "./base/base_uses.f90"
35 :
36 : IMPLICIT NONE
37 :
38 : PRIVATE
39 :
40 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'dm_ls_scf_curvy'
41 :
42 : PUBLIC :: dm_ls_curvy_optimization, deallocate_curvy_data
43 :
44 : CONTAINS
45 :
46 : ! **************************************************************************************************
47 : !> \brief driver routine for Head-Gordon curvy step approach
48 : !> \param ls_scf_env ...
49 : !> \param energy ...
50 : !> \param check_conv ...
51 : !> \par History
52 : !> 2012.05 created [Florian Schiffmann]
53 : !> \author Florian Schiffmann
54 : ! **************************************************************************************************
55 :
56 90 : SUBROUTINE dm_ls_curvy_optimization(ls_scf_env, energy, check_conv)
57 : TYPE(ls_scf_env_type) :: ls_scf_env
58 : REAL(KIND=dp) :: energy
59 : LOGICAL :: check_conv
60 :
61 : CHARACTER(LEN=*), PARAMETER :: routineN = 'dm_ls_curvy_optimization'
62 :
63 : INTEGER :: handle, i, lsstep
64 :
65 90 : CALL timeset(routineN, handle)
66 :
67 90 : CALL cite_reference(Shao2003)
68 :
69 : ! Upon first call initialize all matrices needed curing optimization
70 : ! In addition transform P into orthonormal basis. Will be scaled by 0.5 in closed shell case
71 : ! Only to be done once as it will be stored and reused afterwards
72 : ! TRS4 might yield a non-idempotent P therefore McWeeny purification is applied on initial P
73 :
74 90 : IF (.NOT. ALLOCATED(ls_scf_env%curvy_data%matrix_dp)) THEN
75 18 : CALL init_curvy(ls_scf_env%curvy_data, ls_scf_env%matrix_s, ls_scf_env%nspins)
76 18 : ls_scf_env%curvy_data%line_search_step = 1
77 :
78 18 : IF (ls_scf_env%curvy_data%line_search_type == ls_scf_line_search_3point_2d) THEN
79 6 : DO i = 1, ls_scf_env%nspins
80 : CALL dbcsr_copy(ls_scf_env%curvy_data%matrix_psave(i, 1), &
81 6 : ls_scf_env%matrix_p(i))
82 : END DO
83 : END IF
84 18 : IF (ls_scf_env%nspins == 1) CALL dbcsr_scale(ls_scf_env%matrix_p(1), 0.5_dp)
85 : CALL transform_matrix_orth(ls_scf_env%matrix_p, ls_scf_env%matrix_s_sqrt, &
86 18 : ls_scf_env%eps_filter)
87 18 : CALL purify_mcweeny(ls_scf_env%matrix_p, ls_scf_env%eps_filter, 3)
88 38 : DO i = 1, ls_scf_env%nspins
89 38 : CALL dbcsr_copy(ls_scf_env%curvy_data%matrix_p(i), ls_scf_env%matrix_p(i))
90 : END DO
91 : END IF
92 :
93 90 : lsstep = ls_scf_env%curvy_data%line_search_step
94 :
95 : ! If new search direction has to be computed transform H into the orthnormal basis
96 :
97 90 : IF (ls_scf_env%curvy_data%line_search_step == 1) &
98 : CALL transform_matrix_orth(ls_scf_env%matrix_ks, ls_scf_env%matrix_s_sqrt_inv, &
99 28 : ls_scf_env%eps_filter)
100 :
101 : ! Set the energies for the line search and make sure to give the correct energy back to scf_main
102 90 : ls_scf_env%curvy_data%energies(lsstep) = energy
103 90 : IF (lsstep .NE. 1) energy = ls_scf_env%curvy_data%energies(1)
104 :
105 : ! start the optimization by calling the driver routine or simply combine saved P(2D line search)
106 90 : IF (lsstep .LE. 2) THEN
107 56 : CALL optimization_step(ls_scf_env%curvy_data, ls_scf_env)
108 34 : ELSE IF (lsstep == ls_scf_env%curvy_data%line_search_type) THEN
109 : ! line_search type has the value appropriate to the number of energy calculations needed
110 28 : CALL optimization_step(ls_scf_env%curvy_data, ls_scf_env)
111 : ELSE
112 : CALL new_p_from_save(ls_scf_env%matrix_p, ls_scf_env%curvy_data%matrix_psave, lsstep, &
113 6 : ls_scf_env%curvy_data%double_step_size)
114 6 : ls_scf_env%curvy_data%line_search_step = ls_scf_env%curvy_data%line_search_step + 1
115 6 : CALL timestop(handle)
116 6 : RETURN
117 : END IF
118 84 : lsstep = ls_scf_env%curvy_data%line_search_step
119 :
120 : ! transform new density matrix back into nonorthonormal basis (again scaling might apply)
121 :
122 : CALL transform_matrix_orth(ls_scf_env%matrix_p, ls_scf_env%matrix_s_sqrt_inv, &
123 84 : ls_scf_env%eps_filter)
124 84 : IF (ls_scf_env%nspins == 1) CALL dbcsr_scale(ls_scf_env%matrix_p(1), 2.0_dp)
125 :
126 : ! P-matrices only need to be stored in case of 2D line search
127 84 : IF (lsstep .LE. 3 .AND. ls_scf_env%curvy_data%line_search_type == ls_scf_line_search_3point_2d) THEN
128 18 : DO i = 1, ls_scf_env%nspins
129 : CALL dbcsr_copy(ls_scf_env%curvy_data%matrix_psave(i, lsstep), &
130 18 : ls_scf_env%matrix_p(i))
131 : END DO
132 : END IF
133 84 : check_conv = lsstep == 1
134 :
135 84 : CALL timestop(handle)
136 :
137 : END SUBROUTINE dm_ls_curvy_optimization
138 :
139 : ! **************************************************************************************************
140 : !> \brief low level routine for Head-Gordons curvy step approach
141 : !> computes gradients, performs a cg and line search,
142 : !> and evaluates the BCH series to obtain the new P matrix
143 : !> \param curvy_data ...
144 : !> \param ls_scf_env ...
145 : !> \par History
146 : !> 2012.05 created [Florian Schiffmann]
147 : !> \author Florian Schiffmann
148 : ! **************************************************************************************************
149 :
150 84 : SUBROUTINE optimization_step(curvy_data, ls_scf_env)
151 : TYPE(ls_scf_curvy_type) :: curvy_data
152 : TYPE(ls_scf_env_type) :: ls_scf_env
153 :
154 : CHARACTER(LEN=*), PARAMETER :: routineN = 'optimization_step'
155 :
156 : INTEGER :: handle, ispin
157 : REAL(KIND=dp) :: filter, step_size(2)
158 :
159 : ! Upon first line search step compute new search direction and apply CG if required
160 :
161 84 : CALL timeset(routineN, handle)
162 :
163 84 : IF (curvy_data%line_search_step == 1) THEN
164 196 : curvy_data%step_size = MAXVAL(curvy_data%step_size)
165 84 : curvy_data%step_size = MIN(MAX(0.10_dp, 0.5_dp*ABS(curvy_data%step_size(1))), 0.5_dp)
166 : ! Dynamic eps_filter for newton steps
167 : filter = MAX(ls_scf_env%eps_filter*curvy_data%min_filter, &
168 28 : ls_scf_env%eps_filter*curvy_data%filter_factor)
169 : CALL compute_direction_newton(curvy_data%matrix_p, ls_scf_env%matrix_ks, &
170 : curvy_data%matrix_dp, filter, curvy_data%fix_shift, curvy_data%shift, &
171 28 : curvy_data%cg_numer, curvy_data%cg_denom, curvy_data%min_shift)
172 28 : curvy_data%filter_factor = curvy_data%scale_filter*curvy_data%filter_factor
173 84 : step_size = curvy_data%step_size
174 84 : curvy_data%BCH_saved = 0
175 56 : ELSE IF (curvy_data%line_search_step == 2) THEN
176 84 : step_size = curvy_data%step_size
177 28 : IF (curvy_data%energies(1) - curvy_data%energies(2) .GT. 0.0_dp) THEN
178 72 : curvy_data%step_size = curvy_data%step_size*2.0_dp
179 24 : curvy_data%double_step_size = .TRUE.
180 : ELSE
181 12 : curvy_data%step_size = curvy_data%step_size*0.5_dp
182 4 : curvy_data%double_step_size = .FALSE.
183 : END IF
184 84 : step_size = curvy_data%step_size
185 28 : ELSE IF (curvy_data%line_search_step == ls_scf_line_search_3point_2d) THEN
186 2 : CALL line_search_2d(curvy_data%energies, curvy_data%step_size)
187 6 : step_size = curvy_data%step_size
188 26 : ELSE IF (curvy_data%line_search_step == ls_scf_line_search_3point) THEN
189 26 : CALL line_search_3pnt(curvy_data%energies, curvy_data%step_size)
190 78 : step_size = curvy_data%step_size
191 : END IF
192 :
193 : CALL update_p_exp(curvy_data%matrix_p, ls_scf_env%matrix_p, curvy_data%matrix_dp, &
194 : curvy_data%matrix_BCH, ls_scf_env%eps_filter, step_size, curvy_data%BCH_saved, &
195 84 : curvy_data%n_bch_hist)
196 :
197 : ! line_search type has the value appropriate to the numeber of energy calculations needed
198 84 : curvy_data%line_search_step = MOD(curvy_data%line_search_step, curvy_data%line_search_type) + 1
199 84 : IF (curvy_data%line_search_step == 1) THEN
200 58 : DO ispin = 1, SIZE(curvy_data%matrix_p)
201 58 : CALL dbcsr_copy(curvy_data%matrix_p(ispin), ls_scf_env%matrix_p(ispin))
202 : END DO
203 : END IF
204 84 : CALL timestop(handle)
205 :
206 84 : END SUBROUTINE optimization_step
207 :
208 : ! **************************************************************************************************
209 : !> \brief Perform a 6pnt-2D line search for spin polarized calculations.
210 : !> Fit a 2D parabolic function to 6 points
211 : !> \param energies ...
212 : !> \param step_size ...
213 : !> \par History
214 : !> 2012.05 created [Florian Schiffmann]
215 : !> \author Florian Schiffmann
216 : ! **************************************************************************************************
217 :
218 2 : SUBROUTINE line_search_2d(energies, step_size)
219 : REAL(KIND=dp) :: energies(6), step_size(2)
220 :
221 : INTEGER :: info, unit_nr
222 : REAL(KIND=dp) :: e_pred, param(6), s1, s1sq, s2, s2sq, &
223 : sys_lin_eq(6, 6), tmp_e, v1, v2
224 : TYPE(cp_logger_type), POINTER :: logger
225 :
226 2 : logger => cp_get_default_logger()
227 2 : IF (energies(1) - energies(2) .LT. 0._dp) THEN
228 0 : tmp_e = energies(2); energies(2) = energies(3); energies(3) = tmp_e
229 0 : step_size = step_size*2.0_dp
230 : END IF
231 2 : IF (logger%para_env%is_source()) THEN
232 1 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
233 : ELSE
234 : unit_nr = -1
235 : END IF
236 2 : s1 = 0.5_dp*step_size(1); s2 = step_size(1); s1sq = s1**2; s2sq = s2**2
237 14 : sys_lin_eq = 0.0_dp; sys_lin_eq(:, 6) = 1.0_dp
238 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
239 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
240 2 : sys_lin_eq(4, 3) = s1sq; sys_lin_eq(4, 5) = s1
241 2 : sys_lin_eq(5, 1) = s1sq; sys_lin_eq(5, 4) = s1
242 2 : sys_lin_eq(6, 3) = s2sq; sys_lin_eq(6, 5) = s2
243 :
244 2 : CALL invmat(sys_lin_eq, info)
245 86 : param = MATMUL(sys_lin_eq, energies)
246 2 : v1 = (param(2)*param(4))/(2.0_dp*param(1)) - param(5)
247 2 : v2 = -(param(2)**2)/(2.0_dp*param(1)) + 2.0_dp*param(3)
248 2 : step_size(2) = v1/v2
249 2 : step_size(1) = (-param(2)*step_size(2) - param(4))/(2.0_dp*param(1))
250 2 : IF (step_size(1) .LT. 0.0_dp) step_size(1) = 1.0_dp
251 2 : IF (step_size(2) .LT. 0.0_dp) step_size(2) = 1.0_dp
252 : ! step_size(1)=MIN(step_size(1),2.0_dp)
253 : ! step_size(2)=MIN(step_size(2),2.0_dp)
254 : e_pred = param(1)*step_size(1)**2 + param(2)*step_size(1)*step_size(2) + &
255 2 : param(3)*step_size(2)**2 + param(4)*step_size(1) + param(5)*step_size(2) + param(6)
256 2 : IF (unit_nr .GT. 0) WRITE (unit_nr, "(t3,a,F10.5,F10.5,A,F20.9)") &
257 1 : " Line Search: Step Size", step_size, " Predicted energy", e_pred
258 : e_pred = param(1)*s1**2 + param(2)*s2*s1*0.0_dp + &
259 : param(3)*s1**2*0.0_dp + param(4)*s1 + param(5)*s1*0.0_dp + param(6)
260 :
261 2 : END SUBROUTINE line_search_2d
262 :
263 : ! **************************************************************************************************
264 : !> \brief Perform a 3pnt line search
265 : !> \param energies ...
266 : !> \param step_size ...
267 : !> \par History
268 : !> 2012.05 created [Florian Schiffmann]
269 : !> \author Florian Schiffmann
270 : ! **************************************************************************************************
271 :
272 26 : SUBROUTINE line_search_3pnt(energies, step_size)
273 : REAL(KIND=dp) :: energies(3), step_size(2)
274 :
275 : INTEGER :: unit_nr
276 : REAL(KIND=dp) :: a, b, c, e_pred, min_val, step1, tmp, &
277 : tmp_e
278 : TYPE(cp_logger_type), POINTER :: logger
279 :
280 26 : logger => cp_get_default_logger()
281 26 : IF (energies(1) - energies(2) .LT. 0._dp) THEN
282 4 : tmp_e = energies(2); energies(2) = energies(3); energies(3) = tmp_e
283 12 : step_size = step_size*2.0_dp
284 : END IF
285 26 : IF (logger%para_env%is_source()) THEN
286 13 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
287 : ELSE
288 13 : unit_nr = -1
289 : END IF
290 26 : step1 = 0.5_dp*step_size(1)
291 26 : c = energies(1)
292 26 : a = (energies(3) + c - 2.0_dp*energies(2))/(2.0_dp*step1**2)
293 26 : b = (energies(2) - c - a*step1**2)/step1
294 26 : IF (a .LT. 1.0E-12_dp) a = -1.0E-12_dp
295 26 : min_val = -b/(2.0_dp*a)
296 26 : e_pred = a*min_val**2 + b*min_val + c
297 26 : tmp = step_size(1)
298 26 : IF (e_pred .LT. energies(1) .AND. e_pred .LT. energies(2)) THEN
299 : step_size = MAX(-1.0_dp, &
300 54 : MIN(min_val, 10_dp*step_size))
301 : ELSE
302 24 : step_size = 1.0_dp
303 : END IF
304 26 : e_pred = a*(step_size(1))**2 + b*(step_size(1)) + c
305 26 : IF (unit_nr .GT. 0) THEN
306 13 : WRITE (unit_nr, "(t3,a,f16.8,a,F20.9)") "Line Search: Step Size", step_size(1), " Predicted energy", e_pred
307 13 : CALL m_flush(unit_nr)
308 : END IF
309 26 : END SUBROUTINE line_search_3pnt
310 :
311 : ! **************************************************************************************************
312 : !> \brief Get a new search direction. Iterate to obtain a Newton like step
313 : !> Refine with a CG update of the search direction
314 : !> \param matrix_p ...
315 : !> \param matrix_ks ...
316 : !> \param matrix_dp ...
317 : !> \param eps_filter ...
318 : !> \param fix_shift ...
319 : !> \param curvy_shift ...
320 : !> \param cg_numer ...
321 : !> \param cg_denom ...
322 : !> \param min_shift ...
323 : !> \par History
324 : !> 2012.05 created [Florian Schiffmann]
325 : !> \author Florian Schiffmann
326 : ! **************************************************************************************************
327 :
328 28 : SUBROUTINE compute_direction_newton(matrix_p, matrix_ks, matrix_dp, eps_filter, fix_shift, &
329 : curvy_shift, cg_numer, cg_denom, min_shift)
330 : TYPE(dbcsr_type), DIMENSION(:) :: matrix_p, matrix_ks, matrix_dp
331 : REAL(KIND=dp) :: eps_filter
332 : LOGICAL :: fix_shift(2)
333 : REAL(KIND=dp) :: curvy_shift(2), cg_numer(2), &
334 : cg_denom(2), min_shift
335 :
336 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_direction_newton'
337 :
338 : INTEGER :: handle, i, ispin, ncyc, nspin, unit_nr
339 : LOGICAL :: at_limit
340 : REAL(KIND=dp) :: beta, conv_val, maxel, old_conv, shift
341 : TYPE(cp_logger_type), POINTER :: logger
342 : TYPE(dbcsr_type) :: matrix_Ax, matrix_b, matrix_cg, &
343 : matrix_dp_old, matrix_PKs, matrix_res, &
344 : matrix_tmp, matrix_tmp1
345 :
346 56 : logger => cp_get_default_logger()
347 :
348 28 : IF (logger%para_env%is_source()) THEN
349 14 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
350 : ELSE
351 14 : unit_nr = -1
352 : END IF
353 28 : CALL timeset(routineN, handle)
354 28 : nspin = SIZE(matrix_p)
355 :
356 28 : CALL dbcsr_create(matrix_PKs, template=matrix_dp(1), matrix_type=dbcsr_type_no_symmetry)
357 28 : CALL dbcsr_create(matrix_Ax, template=matrix_dp(1), matrix_type=dbcsr_type_no_symmetry)
358 28 : CALL dbcsr_create(matrix_tmp, template=matrix_dp(1), matrix_type=dbcsr_type_no_symmetry)
359 28 : CALL dbcsr_create(matrix_tmp1, template=matrix_dp(1), matrix_type=dbcsr_type_no_symmetry)
360 28 : CALL dbcsr_create(matrix_res, template=matrix_dp(1), matrix_type=dbcsr_type_no_symmetry)
361 28 : CALL dbcsr_create(matrix_cg, template=matrix_dp(1), matrix_type=dbcsr_type_no_symmetry)
362 28 : CALL dbcsr_create(matrix_b, template=matrix_dp(1), matrix_type=dbcsr_type_no_symmetry)
363 28 : CALL dbcsr_create(matrix_dp_old, template=matrix_dp(1), matrix_type=dbcsr_type_no_symmetry)
364 :
365 58 : DO ispin = 1, nspin
366 30 : CALL dbcsr_copy(matrix_dp_old, matrix_dp(ispin))
367 :
368 : ! Precompute some matrices to save work during iterations
369 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_p(ispin), matrix_ks(ispin), &
370 30 : 0.0_dp, matrix_PKs, filter_eps=eps_filter)
371 30 : CALL dbcsr_transposed(matrix_b, matrix_PKs)
372 30 : CALL dbcsr_copy(matrix_cg, matrix_b)
373 :
374 : ! Starting CG with guess 0-matrix gives -2*gradient=[Ks*P-(Ks*P)T] for cg_matrix in second step
375 30 : CALL dbcsr_add(matrix_cg, matrix_PKs, 2.0_dp, -2.0_dp)
376 :
377 : ! Residual matrix in first step=cg matrix. Keep Pks for later use in CG!
378 30 : CALL dbcsr_copy(matrix_res, matrix_cg)
379 :
380 : ! Precompute -FP-[FP]T which will be used throughout the CG iterations
381 30 : CALL dbcsr_add(matrix_b, matrix_PKs, -1.0_dp, -1.0_dp)
382 :
383 : ! Setup some values to check convergence and safety checks for eigenvalue shifting
384 30 : CALL dbcsr_norm(matrix_res, which_norm=2, norm_scalar=old_conv)
385 30 : old_conv = dbcsr_frobenius_norm(matrix_res)
386 30 : shift = MIN(10.0_dp, MAX(min_shift, 0.05_dp*old_conv))
387 30 : conv_val = MAX(0.010_dp*old_conv, 100.0_dp*eps_filter)
388 30 : old_conv = 100.0_dp
389 30 : IF (fix_shift(ispin)) THEN
390 0 : shift = MAX(min_shift, MIN(10.0_dp, MAX(shift, curvy_shift(ispin) - 0.5_dp*curvy_shift(ispin))))
391 0 : curvy_shift(ispin) = shift
392 : END IF
393 :
394 : ! Begin the real optimization loop
395 30 : CALL dbcsr_set(matrix_dp(ispin), 0.0_dp)
396 30 : ncyc = 10
397 104 : DO i = 1, ncyc
398 :
399 : ! One step to compute: -FPD-DPF-DFP-PFD (not obvious but symmetry allows for some tricks)
400 104 : 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 104 : 0.0_dp, matrix_tmp, filter_eps=eps_filter)
405 104 : CALL commutator_symm(matrix_ks(ispin), matrix_tmp, matrix_tmp1, eps_filter, 2.0_dp)
406 104 : 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 104 : 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 104 : matrix_tmp, eps_filter, at_limit)
413 104 : CALL dbcsr_filter(matrix_cg, eps_filter)
414 :
415 : ! check for convergence of the newton step
416 104 : maxel = dbcsr_frobenius_norm(matrix_res)
417 104 : IF (unit_nr .GT. 0) THEN
418 52 : WRITE (unit_nr, "(T3,A,F12.6)") "Convergence of Newton iteration ", maxel
419 52 : CALL m_flush(unit_nr)
420 : END IF
421 104 : at_limit = at_limit .OR. (old_conv/maxel .LT. 1.01_dp)
422 104 : old_conv = maxel
423 104 : IF (i == ncyc .AND. maxel/conv_val .GT. 5.0_dp) THEN
424 0 : fix_shift(ispin) = .TRUE.
425 0 : curvy_shift(ispin) = 4.0_dp*shift
426 : END IF
427 104 : 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 30 : CALL dbcsr_transposed(matrix_b, matrix_PKs)
432 : !compute b= -2*KsP+2*PKs=-(2*gradient)
433 30 : CALL dbcsr_copy(matrix_cg, matrix_b)
434 30 : CALL dbcsr_add(matrix_cg, matrix_PKs, 1.0_dp, -1.0_dp)
435 30 : cg_denom(ispin) = cg_numer(ispin)
436 30 : CALL dbcsr_dot(matrix_cg, matrix_dp(ispin), cg_numer(ispin))
437 30 : beta = cg_numer(ispin)/MAX(cg_denom(ispin), 1.0E-6_dp)
438 30 : IF (beta .LT. 1.0_dp) THEN
439 28 : beta = MAX(0.0_dp, beta)
440 28 : CALL dbcsr_add(matrix_dp(ispin), matrix_dp_old, 1.0_dp, beta)
441 : END IF
442 58 : IF (unit_nr .GT. 0) WRITE (unit_nr, "(A)") " "
443 : END DO
444 :
445 28 : CALL dbcsr_release(matrix_PKs)
446 28 : CALL dbcsr_release(matrix_dp_old)
447 28 : CALL dbcsr_release(matrix_b)
448 28 : CALL dbcsr_release(matrix_Ax)
449 28 : CALL dbcsr_release(matrix_tmp)
450 28 : CALL dbcsr_release(matrix_tmp1)
451 28 : CALL dbcsr_release(matrix_b)
452 28 : CALL dbcsr_release(matrix_res)
453 28 : CALL dbcsr_release(matrix_cg)
454 :
455 28 : IF (unit_nr .GT. 0) CALL m_flush(unit_nr)
456 28 : CALL timestop(handle)
457 28 : 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 104 : 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 104 : at_limit = .FALSE.
485 104 : CALL dbcsr_dot(res, res, norm_rr)
486 104 : CALL dbcsr_dot(cg, Ax, norm_cA)
487 104 : lin_eq = 0.0_dp
488 104 : fac = norm_rr/norm_cA
489 104 : fac1 = fac
490 : ! Use a 3point line search and a fit to a quadratic function to determine optimal step size
491 416 : DO i = 1, 3
492 312 : CALL dbcsr_copy(tmp, res)
493 312 : CALL dbcsr_add(tmp, Ax, 1.0_dp, -fac)
494 312 : devi(i) = dbcsr_frobenius_norm(tmp)
495 1248 : lin_eq(i, :) = (/fac**2, fac, 1.0_dp/)
496 416 : fac = fac1 + fac1*((-1)**i)*0.5_dp
497 : END DO
498 104 : CALL invmat(lin_eq, info)
499 1352 : vec = MATMUL(lin_eq, devi)
500 104 : alpha = -vec(2)/(2.0_dp*vec(1))
501 104 : fac = SQRT(norm_rr/(norm_cA*alpha))
502 : !scale the previous matrices to match the step size
503 104 : CALL dbcsr_scale(Ax, fac)
504 104 : CALL dbcsr_scale(cg, fac)
505 104 : norm_cA = norm_cA*fac**2
506 :
507 : ! USe CG to get the new matrices
508 104 : alpha = norm_rr/norm_cA
509 104 : CALL dbcsr_add(res, Ax, 1.0_dp, -alpha)
510 104 : CALL dbcsr_dot(res, res, new_norm)
511 104 : 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 22 : at_limit = .TRUE.
514 : ELSE
515 : beta = new_norm/norm_rr
516 82 : CALL dbcsr_add(deltp, cg, 1.0_dp, alpha)
517 : END IF
518 104 : beta = new_norm/norm_rr
519 104 : CALL dbcsr_add(cg, res, beta, 1.0_dp)
520 :
521 104 : 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 6 : 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 8 : SELECT CASE (lsstep)
543 : CASE (3)
544 2 : CALL dbcsr_copy(matrix_p(1), matrix_psave(1, 1))
545 2 : IF (DOUBLE) THEN
546 2 : CALL dbcsr_copy(matrix_p(2), matrix_psave(2, 2))
547 : ELSE
548 0 : CALL dbcsr_copy(matrix_p(2), matrix_psave(2, 3))
549 : END IF
550 : CASE (4)
551 2 : IF (DOUBLE) THEN
552 2 : CALL dbcsr_copy(matrix_p(1), matrix_psave(1, 2))
553 : ELSE
554 0 : CALL dbcsr_copy(matrix_p(1), matrix_psave(1, 3))
555 : END IF
556 2 : CALL dbcsr_copy(matrix_p(2), matrix_psave(2, 1))
557 : CASE (5)
558 2 : CALL dbcsr_copy(matrix_p(1), matrix_psave(1, 1))
559 8 : IF (DOUBLE) THEN
560 2 : CALL dbcsr_copy(matrix_p(2), matrix_psave(2, 3))
561 : ELSE
562 0 : CALL dbcsr_copy(matrix_p(2), matrix_psave(2, 2))
563 : END IF
564 : END SELECT
565 :
566 6 : 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 208 : 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 208 : CALL timeset(routineN, handle)
590 :
591 208 : CALL dbcsr_create(work, template=a, matrix_type=dbcsr_type_no_symmetry)
592 :
593 208 : CALL dbcsr_multiply("N", "N", prefac, a, b, 0.0_dp, res, filter_eps=eps_filter)
594 208 : CALL dbcsr_transposed(work, res)
595 208 : CALL dbcsr_add(res, work, 1.0_dp, -1.0_dp)
596 :
597 208 : CALL dbcsr_release(work)
598 :
599 208 : CALL timestop(handle)
600 208 : 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 84 : 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 84 : CALL timeset(routineN, handle)
635 :
636 84 : logger => cp_get_default_logger()
637 84 : IF (logger%para_env%is_source()) THEN
638 42 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
639 : ELSE
640 42 : unit_nr = -1
641 : END IF
642 :
643 84 : CALL dbcsr_create(matrix, template=matrix_p_in(1), matrix_type=dbcsr_type_no_symmetry)
644 84 : CALL dbcsr_create(matrix_tmp, template=matrix_p_in(1), matrix_type=dbcsr_type_no_symmetry)
645 84 : nspin = SIZE(matrix_p_in)
646 :
647 174 : DO ispin = 1, nspin
648 90 : step_fac = 1.0_dp
649 90 : frob_norm = 1.0_dp
650 90 : nsave = 0
651 :
652 90 : CALL dbcsr_copy(matrix_tmp, matrix_p_in(ispin))
653 90 : 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 130 : DO i = 1, BCH_saved(ispin)
657 86 : step_fac = step_fac*step_size(ispin)
658 86 : CALL dbcsr_copy(matrix_tmp, matrix_p_out(ispin))
659 86 : CALL dbcsr_add(matrix_p_out(ispin), matrix_BCH(ispin, i), 1.0_dp, ifac(i)*step_fac)
660 86 : CALL dbcsr_add(matrix_tmp, matrix_p_out(ispin), 1.0_dp, -1.0_dp)
661 86 : frob_norm = dbcsr_frobenius_norm(matrix_tmp)
662 86 : 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 130 : IF (frob_norm .LT. threshold) EXIT
664 : END DO
665 90 : 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 44 : save_BCH = BCH_saved(ispin) == 0 .AND. n_bch_hist .GT. 0
669 86 : DO i = BCH_saved(ispin) + 1, 20
670 86 : 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 86 : 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 86 : CALL dbcsr_transposed(matrix_tmp, matrix)
679 86 : 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 86 : CALL dbcsr_copy(matrix_tmp, matrix_p_out(ispin))
683 86 : CALL dbcsr_add(matrix_p_out(ispin), matrix, 1.0_dp, ifac(i)*step_fac)
684 86 : IF (save_BCH .AND. i .LE. n_bch_hist) THEN
685 78 : CALL dbcsr_copy(matrix_BCH(ispin, i), matrix)
686 78 : nsave = i
687 : END IF
688 :
689 86 : 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 86 : frob_norm = dbcsr_frobenius_norm(matrix_tmp)
693 86 : 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 86 : 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 42 : CALL dbcsr_copy(matrix_tmp, matrix)
698 86 : CALL dbcsr_filter(matrix_tmp, threshold)
699 : END DO
700 44 : BCH_saved(ispin) = nsave
701 128 : IF (unit_nr .GT. 0) WRITE (unit_nr, "(A)") " "
702 : END DO
703 :
704 84 : CALL purify_mcweeny(matrix_p_out, threshold, 1)
705 84 : IF (unit_nr .GT. 0) CALL m_flush(unit_nr)
706 84 : CALL dbcsr_release(matrix_tmp)
707 84 : CALL dbcsr_release(matrix)
708 84 : CALL timestop(handle)
709 84 : 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 130 : 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 130 : CALL timeset(routineN, handle)
733 :
734 130 : CALL dbcsr_create(matrix_work, template=matrix(1), matrix_type=dbcsr_type_no_symmetry)
735 130 : CALL dbcsr_create(matrix_tmp, template=matrix(1), matrix_type=dbcsr_type_no_symmetry)
736 :
737 270 : DO ispin = 1, SIZE(matrix)
738 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix(ispin), matrix_trafo, &
739 140 : 0.0_dp, matrix_work, filter_eps=eps_filter)
740 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_trafo, matrix_work, &
741 140 : 0.0_dp, matrix_tmp, filter_eps=eps_filter)
742 : ! symmetrize results (this is again needed to make sure everything is stable)
743 140 : CALL dbcsr_transposed(matrix_work, matrix_tmp)
744 140 : CALL dbcsr_add(matrix_tmp, matrix_work, 0.5_dp, 0.5_dp)
745 270 : CALL dbcsr_copy(matrix(ispin), matrix_tmp)
746 : END DO
747 :
748 130 : CALL dbcsr_release(matrix_tmp)
749 130 : CALL dbcsr_release(matrix_work)
750 130 : CALL timestop(handle)
751 :
752 130 : END SUBROUTINE
753 :
754 : ! **************************************************************************************************
755 : !> \brief ...
756 : !> \param curvy_data ...
757 : ! **************************************************************************************************
758 814 : SUBROUTINE deallocate_curvy_data(curvy_data)
759 : TYPE(ls_scf_curvy_type) :: curvy_data
760 :
761 : INTEGER :: i, j
762 :
763 814 : CALL release_dbcsr_array(curvy_data%matrix_dp)
764 814 : CALL release_dbcsr_array(curvy_data%matrix_p)
765 :
766 814 : IF (ALLOCATED(curvy_data%matrix_psave)) THEN
767 6 : DO i = 1, SIZE(curvy_data%matrix_psave, 1)
768 18 : DO j = 1, 3
769 16 : CALL dbcsr_release(curvy_data%matrix_psave(i, j))
770 : END DO
771 : END DO
772 2 : DEALLOCATE (curvy_data%matrix_psave)
773 : END IF
774 814 : IF (ALLOCATED(curvy_data%matrix_BCH)) THEN
775 38 : DO i = 1, SIZE(curvy_data%matrix_BCH, 1)
776 178 : DO j = 1, 7
777 160 : CALL dbcsr_release(curvy_data%matrix_BCH(i, j))
778 : END DO
779 : END DO
780 18 : DEALLOCATE (curvy_data%matrix_BCH)
781 : END IF
782 814 : END SUBROUTINE deallocate_curvy_data
783 :
784 : ! **************************************************************************************************
785 : !> \brief ...
786 : !> \param matrix ...
787 : ! **************************************************************************************************
788 1628 : SUBROUTINE release_dbcsr_array(matrix)
789 : TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:) :: matrix
790 :
791 : INTEGER :: i
792 :
793 1628 : IF (ALLOCATED(matrix)) THEN
794 76 : DO i = 1, SIZE(matrix)
795 76 : CALL dbcsr_release(matrix(i))
796 : END DO
797 36 : DEALLOCATE (matrix)
798 : END IF
799 1628 : END SUBROUTINE release_dbcsr_array
800 :
801 : ! **************************************************************************************************
802 : !> \brief ...
803 : !> \param curvy_data ...
804 : !> \param matrix_s ...
805 : !> \param nspins ...
806 : ! **************************************************************************************************
807 18 : 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 74 : ALLOCATE (curvy_data%matrix_dp(nspins))
815 74 : ALLOCATE (curvy_data%matrix_p(nspins))
816 38 : DO ispin = 1, nspins
817 : CALL dbcsr_create(curvy_data%matrix_dp(ispin), template=matrix_s, &
818 20 : matrix_type=dbcsr_type_no_symmetry)
819 20 : CALL dbcsr_set(curvy_data%matrix_dp(ispin), 0.0_dp)
820 : CALL dbcsr_create(curvy_data%matrix_p(ispin), template=matrix_s, &
821 20 : matrix_type=dbcsr_type_no_symmetry)
822 60 : curvy_data%fix_shift = .FALSE.
823 20 : curvy_data%double_step_size = .TRUE.
824 60 : curvy_data%shift = 1.0_dp
825 60 : curvy_data%BCH_saved = 0
826 60 : curvy_data%step_size = 0.60_dp
827 60 : curvy_data%cg_numer = 0.00_dp
828 78 : curvy_data%cg_denom = 0.00_dp
829 : END DO
830 18 : IF (curvy_data%line_search_type == ls_scf_line_search_3point_2d) THEN
831 24 : ALLOCATE (curvy_data%matrix_psave(nspins, 3))
832 6 : DO ispin = 1, nspins
833 18 : DO j = 1, 3
834 : CALL dbcsr_create(curvy_data%matrix_psave(ispin, j), template=matrix_s, &
835 16 : matrix_type=dbcsr_type_no_symmetry)
836 : END DO
837 : END DO
838 : END IF
839 18 : IF (curvy_data%n_bch_hist .GT. 0) THEN
840 338 : ALLOCATE (curvy_data%matrix_BCH(nspins, curvy_data%n_bch_hist))
841 38 : DO ispin = 1, nspins
842 178 : DO j = 1, curvy_data%n_bch_hist
843 : CALL dbcsr_create(curvy_data%matrix_BCH(ispin, j), template=matrix_s, &
844 160 : matrix_type=dbcsr_type_no_symmetry)
845 : END DO
846 : END DO
847 : END IF
848 :
849 18 : END SUBROUTINE init_curvy
850 :
851 : END MODULE dm_ls_scf_curvy
|
