Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief density matrix optimization using exponential transformations
10 : !> \par History
11 : !> 2012.05 created [Florian Schiffmann]
12 : !> \author Florian Schiffmann
13 : ! **************************************************************************************************
14 :
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_dot, dbcsr_filter, dbcsr_multiply, &
20 : dbcsr_release, dbcsr_scale, dbcsr_set, dbcsr_transposed, dbcsr_type, dbcsr_type_no_symmetry
21 : USE cp_dbcsr_contrib, ONLY: dbcsr_frobenius_norm
22 : USE cp_log_handling, ONLY: cp_get_default_logger,&
23 : cp_logger_get_default_unit_nr,&
24 : cp_logger_type
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 66 : curvy_data%step_size = curvy_data%step_size*2.0_dp
179 22 : curvy_data%double_step_size = .TRUE.
180 : ELSE
181 18 : curvy_data%step_size = curvy_data%step_size*0.5_dp
182 6 : 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 : old_conv = dbcsr_frobenius_norm(matrix_res)
385 30 : shift = MIN(10.0_dp, MAX(min_shift, 0.05_dp*old_conv))
386 30 : conv_val = MAX(0.010_dp*old_conv, 100.0_dp*eps_filter)
387 30 : old_conv = 100.0_dp
388 30 : IF (fix_shift(ispin)) THEN
389 0 : shift = MAX(min_shift, MIN(10.0_dp, MAX(shift, curvy_shift(ispin) - 0.5_dp*curvy_shift(ispin))))
390 0 : curvy_shift(ispin) = shift
391 : END IF
392 :
393 : ! Begin the real optimization loop
394 30 : CALL dbcsr_set(matrix_dp(ispin), 0.0_dp)
395 30 : ncyc = 10
396 104 : DO i = 1, ncyc
397 :
398 : ! One step to compute: -FPD-DPF-DFP-PFD (not obvious but symmetry allows for some tricks)
399 104 : CALL commutator_symm(matrix_b, matrix_cg, matrix_Ax, eps_filter, 1.0_dp)
400 :
401 : ! Compute the missing bits 2*(FDP+PDF) (again use symmetry to compute as a commutator)
402 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_cg, matrix_p(ispin), &
403 104 : 0.0_dp, matrix_tmp, filter_eps=eps_filter)
404 104 : CALL commutator_symm(matrix_ks(ispin), matrix_tmp, matrix_tmp1, eps_filter, 2.0_dp)
405 104 : CALL dbcsr_add(matrix_Ax, matrix_tmp1, 1.0_dp, 1.0_dp)
406 :
407 : ! Apply the shift and hope it's enough to stabilize the CG iterations
408 104 : CALL dbcsr_add(matrix_Ax, matrix_cg, 1.0_dp, shift)
409 :
410 : CALL compute_cg_matrices(matrix_Ax, matrix_res, matrix_cg, matrix_dp(ispin), &
411 104 : matrix_tmp, eps_filter, at_limit)
412 104 : CALL dbcsr_filter(matrix_cg, eps_filter)
413 :
414 : ! check for convergence of the newton step
415 104 : maxel = dbcsr_frobenius_norm(matrix_res)
416 104 : IF (unit_nr .GT. 0) THEN
417 52 : WRITE (unit_nr, "(T3,A,F12.6)") "Convergence of Newton iteration ", maxel
418 52 : CALL m_flush(unit_nr)
419 : END IF
420 104 : at_limit = at_limit .OR. (old_conv/maxel .LT. 1.01_dp)
421 104 : old_conv = maxel
422 104 : IF (i == ncyc .AND. maxel/conv_val .GT. 5.0_dp) THEN
423 0 : fix_shift(ispin) = .TRUE.
424 0 : curvy_shift(ispin) = 4.0_dp*shift
425 : END IF
426 104 : IF (maxel .LT. conv_val .OR. at_limit) EXIT
427 : END DO
428 :
429 : ! Refine the Newton like search direction with a preconditioned cg update
430 30 : CALL dbcsr_transposed(matrix_b, matrix_PKs)
431 : !compute b= -2*KsP+2*PKs=-(2*gradient)
432 30 : CALL dbcsr_copy(matrix_cg, matrix_b)
433 30 : CALL dbcsr_add(matrix_cg, matrix_PKs, 1.0_dp, -1.0_dp)
434 30 : cg_denom(ispin) = cg_numer(ispin)
435 30 : CALL dbcsr_dot(matrix_cg, matrix_dp(ispin), cg_numer(ispin))
436 30 : beta = cg_numer(ispin)/MAX(cg_denom(ispin), 1.0E-6_dp)
437 30 : IF (beta .LT. 1.0_dp) THEN
438 28 : beta = MAX(0.0_dp, beta)
439 28 : CALL dbcsr_add(matrix_dp(ispin), matrix_dp_old, 1.0_dp, beta)
440 : END IF
441 58 : IF (unit_nr .GT. 0) WRITE (unit_nr, "(A)") " "
442 : END DO
443 :
444 28 : CALL dbcsr_release(matrix_PKs)
445 28 : CALL dbcsr_release(matrix_dp_old)
446 28 : CALL dbcsr_release(matrix_b)
447 28 : CALL dbcsr_release(matrix_Ax)
448 28 : CALL dbcsr_release(matrix_tmp)
449 28 : CALL dbcsr_release(matrix_tmp1)
450 28 : CALL dbcsr_release(matrix_b)
451 28 : CALL dbcsr_release(matrix_res)
452 28 : CALL dbcsr_release(matrix_cg)
453 :
454 28 : IF (unit_nr .GT. 0) CALL m_flush(unit_nr)
455 28 : CALL timestop(handle)
456 28 : END SUBROUTINE compute_direction_newton
457 :
458 : ! **************************************************************************************************
459 : !> \brief compute the optimal step size of the current cycle and update the
460 : !> matrices needed to solve the system of linear equations
461 : !> \param Ax ...
462 : !> \param res ...
463 : !> \param cg ...
464 : !> \param deltp ...
465 : !> \param tmp ...
466 : !> \param eps_filter ...
467 : !> \param at_limit ...
468 : !> \par History
469 : !> 2012.05 created [Florian Schiffmann]
470 : !> \author Florian Schiffmann
471 : ! **************************************************************************************************
472 :
473 104 : SUBROUTINE compute_cg_matrices(Ax, res, cg, deltp, tmp, eps_filter, at_limit)
474 : TYPE(dbcsr_type) :: Ax, res, cg, deltp, tmp
475 : REAL(KIND=dp) :: eps_filter
476 : LOGICAL :: at_limit
477 :
478 : INTEGER :: i, info
479 : REAL(KIND=dp) :: alpha, beta, devi(3), fac, fac1, &
480 : lin_eq(3, 3), new_norm, norm_cA, &
481 : norm_rr, vec(3)
482 :
483 104 : at_limit = .FALSE.
484 104 : CALL dbcsr_dot(res, res, norm_rr)
485 104 : CALL dbcsr_dot(cg, Ax, norm_cA)
486 104 : lin_eq = 0.0_dp
487 104 : fac = norm_rr/norm_cA
488 104 : fac1 = fac
489 : ! Use a 3point line search and a fit to a quadratic function to determine optimal step size
490 416 : DO i = 1, 3
491 312 : CALL dbcsr_copy(tmp, res)
492 312 : CALL dbcsr_add(tmp, Ax, 1.0_dp, -fac)
493 312 : devi(i) = dbcsr_frobenius_norm(tmp)
494 1248 : lin_eq(i, :) = (/fac**2, fac, 1.0_dp/)
495 416 : fac = fac1 + fac1*((-1)**i)*0.5_dp
496 : END DO
497 104 : CALL invmat(lin_eq, info)
498 1352 : vec = MATMUL(lin_eq, devi)
499 104 : alpha = -vec(2)/(2.0_dp*vec(1))
500 104 : fac = SQRT(norm_rr/(norm_cA*alpha))
501 : !scale the previous matrices to match the step size
502 104 : CALL dbcsr_scale(Ax, fac)
503 104 : CALL dbcsr_scale(cg, fac)
504 104 : norm_cA = norm_cA*fac**2
505 :
506 : ! USe CG to get the new matrices
507 104 : alpha = norm_rr/norm_cA
508 104 : CALL dbcsr_add(res, Ax, 1.0_dp, -alpha)
509 104 : CALL dbcsr_dot(res, res, new_norm)
510 104 : IF (norm_rr .LT. eps_filter*0.001_dp .OR. new_norm .LT. eps_filter*0.001_dp) THEN
511 : beta = 0.0_dp
512 22 : at_limit = .TRUE.
513 : ELSE
514 : beta = new_norm/norm_rr
515 82 : CALL dbcsr_add(deltp, cg, 1.0_dp, alpha)
516 : END IF
517 104 : beta = new_norm/norm_rr
518 104 : CALL dbcsr_add(cg, res, beta, 1.0_dp)
519 :
520 104 : END SUBROUTINE compute_cg_matrices
521 :
522 : ! **************************************************************************************************
523 : !> \brief Only for 2D line search. Use saved P-components to construct new
524 : !> test density matrix. Takes care as well, whether step_size
525 : !> increased or decreased during 2nd step and combines matrices accordingly
526 : !> \param matrix_p ...
527 : !> \param matrix_psave ...
528 : !> \param lsstep ...
529 : !> \param DOUBLE ...
530 : !> \par History
531 : !> 2012.05 created [Florian Schiffmann]
532 : !> \author Florian Schiffmann
533 : ! **************************************************************************************************
534 :
535 6 : SUBROUTINE new_p_from_save(matrix_p, matrix_psave, lsstep, DOUBLE)
536 : TYPE(dbcsr_type), DIMENSION(:) :: matrix_p
537 : TYPE(dbcsr_type), DIMENSION(:, :) :: matrix_psave
538 : INTEGER :: lsstep
539 : LOGICAL :: DOUBLE
540 :
541 8 : SELECT CASE (lsstep)
542 : CASE (3)
543 2 : CALL dbcsr_copy(matrix_p(1), matrix_psave(1, 1))
544 2 : IF (DOUBLE) THEN
545 2 : CALL dbcsr_copy(matrix_p(2), matrix_psave(2, 2))
546 : ELSE
547 0 : CALL dbcsr_copy(matrix_p(2), matrix_psave(2, 3))
548 : END IF
549 : CASE (4)
550 2 : IF (DOUBLE) THEN
551 2 : CALL dbcsr_copy(matrix_p(1), matrix_psave(1, 2))
552 : ELSE
553 0 : CALL dbcsr_copy(matrix_p(1), matrix_psave(1, 3))
554 : END IF
555 2 : CALL dbcsr_copy(matrix_p(2), matrix_psave(2, 1))
556 : CASE (5)
557 2 : CALL dbcsr_copy(matrix_p(1), matrix_psave(1, 1))
558 8 : IF (DOUBLE) THEN
559 2 : CALL dbcsr_copy(matrix_p(2), matrix_psave(2, 3))
560 : ELSE
561 0 : CALL dbcsr_copy(matrix_p(2), matrix_psave(2, 2))
562 : END IF
563 : END SELECT
564 :
565 6 : END SUBROUTINE new_p_from_save
566 :
567 : ! **************************************************************************************************
568 : !> \brief computes a commutator exploiting symmetry RES=k*[A,B]=k*[AB-(AB)T]
569 : !> \param a ...
570 : !> \param b ...
571 : !> \param res ...
572 : !> \param eps_filter filtering threshold for sparse matrices
573 : !> \param prefac prefactor k in above equation
574 : !> \par History
575 : !> 2012.05 created [Florian Schiffmann]
576 : !> \author Florian Schiffmann
577 : ! **************************************************************************************************
578 :
579 208 : SUBROUTINE commutator_symm(a, b, res, eps_filter, prefac)
580 : TYPE(dbcsr_type) :: a, b, res
581 : REAL(KIND=dp) :: eps_filter, prefac
582 :
583 : CHARACTER(LEN=*), PARAMETER :: routineN = 'commutator_symm'
584 :
585 : INTEGER :: handle
586 : TYPE(dbcsr_type) :: work
587 :
588 208 : CALL timeset(routineN, handle)
589 :
590 208 : CALL dbcsr_create(work, template=a, matrix_type=dbcsr_type_no_symmetry)
591 :
592 208 : CALL dbcsr_multiply("N", "N", prefac, a, b, 0.0_dp, res, filter_eps=eps_filter)
593 208 : CALL dbcsr_transposed(work, res)
594 208 : CALL dbcsr_add(res, work, 1.0_dp, -1.0_dp)
595 :
596 208 : CALL dbcsr_release(work)
597 :
598 208 : CALL timestop(handle)
599 208 : END SUBROUTINE commutator_symm
600 :
601 : ! **************************************************************************************************
602 : !> \brief Use the BCH update to get the new idempotent P
603 : !> Numerics don't allow for perfect idempotency, therefore a mc weeny
604 : !> step is used to make sure we stay close to the idempotent surface
605 : !> \param matrix_p_in ...
606 : !> \param matrix_p_out ...
607 : !> \param matrix_dp ...
608 : !> \param matrix_BCH ...
609 : !> \param threshold ...
610 : !> \param step_size ...
611 : !> \param BCH_saved ...
612 : !> \param n_bch_hist ...
613 : !> \par History
614 : !> 2012.05 created [Florian Schiffmann]
615 : !> \author Florian Schiffmann
616 : ! **************************************************************************************************
617 :
618 84 : SUBROUTINE update_p_exp(matrix_p_in, matrix_p_out, matrix_dp, matrix_BCH, threshold, step_size, &
619 : BCH_saved, n_bch_hist)
620 : TYPE(dbcsr_type), DIMENSION(:) :: matrix_p_in, matrix_p_out, matrix_dp
621 : TYPE(dbcsr_type), DIMENSION(:, :) :: matrix_BCH
622 : REAL(KIND=dp) :: threshold, step_size(2)
623 : INTEGER :: BCH_saved(2), n_bch_hist
624 :
625 : CHARACTER(LEN=*), PARAMETER :: routineN = 'update_p_exp'
626 :
627 : INTEGER :: handle, i, ispin, nsave, nspin, unit_nr
628 : LOGICAL :: save_BCH
629 : REAL(KIND=dp) :: frob_norm, step_fac
630 : TYPE(cp_logger_type), POINTER :: logger
631 : TYPE(dbcsr_type) :: matrix, matrix_tmp
632 :
633 84 : CALL timeset(routineN, handle)
634 :
635 84 : logger => cp_get_default_logger()
636 84 : IF (logger%para_env%is_source()) THEN
637 42 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
638 : ELSE
639 42 : unit_nr = -1
640 : END IF
641 :
642 84 : CALL dbcsr_create(matrix, template=matrix_p_in(1), matrix_type=dbcsr_type_no_symmetry)
643 84 : CALL dbcsr_create(matrix_tmp, template=matrix_p_in(1), matrix_type=dbcsr_type_no_symmetry)
644 84 : nspin = SIZE(matrix_p_in)
645 :
646 174 : DO ispin = 1, nspin
647 90 : step_fac = 1.0_dp
648 90 : frob_norm = 1.0_dp
649 90 : nsave = 0
650 :
651 90 : CALL dbcsr_copy(matrix_tmp, matrix_p_in(ispin))
652 90 : CALL dbcsr_copy(matrix_p_out(ispin), matrix_p_in(ispin))
653 : ! If a BCH history is used make good use of it and do a few steps as a copy and scale update of P
654 : ! else BCH_saved will be 0 and loop is skipped
655 130 : DO i = 1, BCH_saved(ispin)
656 86 : step_fac = step_fac*step_size(ispin)
657 86 : CALL dbcsr_copy(matrix_tmp, matrix_p_out(ispin))
658 86 : CALL dbcsr_add(matrix_p_out(ispin), matrix_BCH(ispin, i), 1.0_dp, ifac(i)*step_fac)
659 86 : CALL dbcsr_add(matrix_tmp, matrix_p_out(ispin), 1.0_dp, -1.0_dp)
660 86 : frob_norm = dbcsr_frobenius_norm(matrix_tmp)
661 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
662 130 : IF (frob_norm .LT. threshold) EXIT
663 : END DO
664 90 : IF (frob_norm .LT. threshold) CYCLE
665 :
666 : ! 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
667 44 : save_BCH = BCH_saved(ispin) == 0 .AND. n_bch_hist .GT. 0
668 86 : DO i = BCH_saved(ispin) + 1, 20
669 86 : step_fac = step_fac*step_size(ispin)
670 : !allow for a bit of matrix magic here by exploiting matrix and matrix_tmp
671 : !matrix_tmp is alway the previous order of the BCH series
672 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_tmp, matrix_dp(ispin), &
673 86 : 0.0_dp, matrix, filter_eps=threshold)
674 :
675 : !(anti)symmetry allows to sum the transposed instead of the full commutator, matrix becomes the latest result
676 :
677 86 : CALL dbcsr_transposed(matrix_tmp, matrix)
678 86 : CALL dbcsr_add(matrix, matrix_tmp, 1.0_dp, 1.0_dp)
679 :
680 : !Finally, add the new BCH order to P, but store the previous one for a convergence check
681 86 : CALL dbcsr_copy(matrix_tmp, matrix_p_out(ispin))
682 86 : CALL dbcsr_add(matrix_p_out(ispin), matrix, 1.0_dp, ifac(i)*step_fac)
683 86 : IF (save_BCH .AND. i .LE. n_bch_hist) THEN
684 78 : CALL dbcsr_copy(matrix_BCH(ispin, i), matrix)
685 78 : nsave = i
686 : END IF
687 :
688 86 : CALL dbcsr_add(matrix_tmp, matrix_p_out(ispin), 1.0_dp, -1.0_dp)
689 :
690 : !Stop the BCH-series if two successive P's differ by less the threshold
691 86 : frob_norm = dbcsr_frobenius_norm(matrix_tmp)
692 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
693 86 : IF (frob_norm .LT. threshold) EXIT
694 :
695 : !Copy the latest BCH-matrix on matrix tmp, so we can cycle with all matrices in place
696 42 : CALL dbcsr_copy(matrix_tmp, matrix)
697 86 : CALL dbcsr_filter(matrix_tmp, threshold)
698 : END DO
699 44 : BCH_saved(ispin) = nsave
700 128 : IF (unit_nr .GT. 0) WRITE (unit_nr, "(A)") " "
701 : END DO
702 :
703 84 : CALL purify_mcweeny(matrix_p_out, threshold, 1)
704 84 : IF (unit_nr .GT. 0) CALL m_flush(unit_nr)
705 84 : CALL dbcsr_release(matrix_tmp)
706 84 : CALL dbcsr_release(matrix)
707 84 : CALL timestop(handle)
708 84 : END SUBROUTINE update_p_exp
709 :
710 : ! **************************************************************************************************
711 : !> \brief performs a transformation of a matrix back to/into orthonormal basis
712 : !> in case of P a scaling of 0.5 has to be applied for closed shell case
713 : !> \param matrix matrix to be transformed
714 : !> \param matrix_trafo transformation matrix
715 : !> \param eps_filter filtering threshold for sparse matrices
716 : !> \par History
717 : !> 2012.05 created [Florian Schiffmann]
718 : !> \author Florian Schiffmann
719 : ! **************************************************************************************************
720 :
721 130 : SUBROUTINE transform_matrix_orth(matrix, matrix_trafo, eps_filter)
722 : TYPE(dbcsr_type), DIMENSION(:) :: matrix
723 : TYPE(dbcsr_type) :: matrix_trafo
724 : REAL(KIND=dp) :: eps_filter
725 :
726 : CHARACTER(LEN=*), PARAMETER :: routineN = 'transform_matrix_orth'
727 :
728 : INTEGER :: handle, ispin
729 : TYPE(dbcsr_type) :: matrix_tmp, matrix_work
730 :
731 130 : CALL timeset(routineN, handle)
732 :
733 130 : CALL dbcsr_create(matrix_work, template=matrix(1), matrix_type=dbcsr_type_no_symmetry)
734 130 : CALL dbcsr_create(matrix_tmp, template=matrix(1), matrix_type=dbcsr_type_no_symmetry)
735 :
736 270 : DO ispin = 1, SIZE(matrix)
737 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix(ispin), matrix_trafo, &
738 140 : 0.0_dp, matrix_work, filter_eps=eps_filter)
739 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_trafo, matrix_work, &
740 140 : 0.0_dp, matrix_tmp, filter_eps=eps_filter)
741 : ! symmetrize results (this is again needed to make sure everything is stable)
742 140 : CALL dbcsr_transposed(matrix_work, matrix_tmp)
743 140 : CALL dbcsr_add(matrix_tmp, matrix_work, 0.5_dp, 0.5_dp)
744 270 : CALL dbcsr_copy(matrix(ispin), matrix_tmp)
745 : END DO
746 :
747 130 : CALL dbcsr_release(matrix_tmp)
748 130 : CALL dbcsr_release(matrix_work)
749 130 : CALL timestop(handle)
750 :
751 130 : END SUBROUTINE
752 :
753 : ! **************************************************************************************************
754 : !> \brief ...
755 : !> \param curvy_data ...
756 : ! **************************************************************************************************
757 882 : SUBROUTINE deallocate_curvy_data(curvy_data)
758 : TYPE(ls_scf_curvy_type) :: curvy_data
759 :
760 : INTEGER :: i, j
761 :
762 882 : CALL release_dbcsr_array(curvy_data%matrix_dp)
763 882 : CALL release_dbcsr_array(curvy_data%matrix_p)
764 :
765 882 : IF (ALLOCATED(curvy_data%matrix_psave)) THEN
766 6 : DO i = 1, SIZE(curvy_data%matrix_psave, 1)
767 18 : DO j = 1, 3
768 16 : CALL dbcsr_release(curvy_data%matrix_psave(i, j))
769 : END DO
770 : END DO
771 2 : DEALLOCATE (curvy_data%matrix_psave)
772 : END IF
773 882 : IF (ALLOCATED(curvy_data%matrix_BCH)) THEN
774 38 : DO i = 1, SIZE(curvy_data%matrix_BCH, 1)
775 178 : DO j = 1, 7
776 160 : CALL dbcsr_release(curvy_data%matrix_BCH(i, j))
777 : END DO
778 : END DO
779 18 : DEALLOCATE (curvy_data%matrix_BCH)
780 : END IF
781 882 : END SUBROUTINE deallocate_curvy_data
782 :
783 : ! **************************************************************************************************
784 : !> \brief ...
785 : !> \param matrix ...
786 : ! **************************************************************************************************
787 1764 : SUBROUTINE release_dbcsr_array(matrix)
788 : TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:) :: matrix
789 :
790 : INTEGER :: i
791 :
792 1764 : IF (ALLOCATED(matrix)) THEN
793 76 : DO i = 1, SIZE(matrix)
794 76 : CALL dbcsr_release(matrix(i))
795 : END DO
796 36 : DEALLOCATE (matrix)
797 : END IF
798 1764 : END SUBROUTINE release_dbcsr_array
799 :
800 : ! **************************************************************************************************
801 : !> \brief ...
802 : !> \param curvy_data ...
803 : !> \param matrix_s ...
804 : !> \param nspins ...
805 : ! **************************************************************************************************
806 18 : SUBROUTINE init_curvy(curvy_data, matrix_s, nspins)
807 : TYPE(ls_scf_curvy_type) :: curvy_data
808 : TYPE(dbcsr_type) :: matrix_s
809 : INTEGER :: nspins
810 :
811 : INTEGER :: ispin, j
812 :
813 74 : ALLOCATE (curvy_data%matrix_dp(nspins))
814 74 : ALLOCATE (curvy_data%matrix_p(nspins))
815 38 : DO ispin = 1, nspins
816 : CALL dbcsr_create(curvy_data%matrix_dp(ispin), template=matrix_s, &
817 20 : matrix_type=dbcsr_type_no_symmetry)
818 20 : CALL dbcsr_set(curvy_data%matrix_dp(ispin), 0.0_dp)
819 : CALL dbcsr_create(curvy_data%matrix_p(ispin), template=matrix_s, &
820 20 : matrix_type=dbcsr_type_no_symmetry)
821 60 : curvy_data%fix_shift = .FALSE.
822 20 : curvy_data%double_step_size = .TRUE.
823 60 : curvy_data%shift = 1.0_dp
824 60 : curvy_data%BCH_saved = 0
825 60 : curvy_data%step_size = 0.60_dp
826 60 : curvy_data%cg_numer = 0.00_dp
827 78 : curvy_data%cg_denom = 0.00_dp
828 : END DO
829 18 : IF (curvy_data%line_search_type == ls_scf_line_search_3point_2d) THEN
830 24 : ALLOCATE (curvy_data%matrix_psave(nspins, 3))
831 6 : DO ispin = 1, nspins
832 18 : DO j = 1, 3
833 : CALL dbcsr_create(curvy_data%matrix_psave(ispin, j), template=matrix_s, &
834 16 : matrix_type=dbcsr_type_no_symmetry)
835 : END DO
836 : END DO
837 : END IF
838 18 : IF (curvy_data%n_bch_hist .GT. 0) THEN
839 338 : ALLOCATE (curvy_data%matrix_BCH(nspins, curvy_data%n_bch_hist))
840 38 : DO ispin = 1, nspins
841 178 : DO j = 1, curvy_data%n_bch_hist
842 : CALL dbcsr_create(curvy_data%matrix_BCH(ispin, j), template=matrix_s, &
843 160 : matrix_type=dbcsr_type_no_symmetry)
844 : END DO
845 : END DO
846 : END IF
847 :
848 18 : END SUBROUTINE init_curvy
849 :
850 : END MODULE dm_ls_scf_curvy
|
