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 Methods used by pao_main.F
10 : !> \author Ole Schuett
11 : ! **************************************************************************************************
12 : MODULE pao_methods
13 : USE ao_util, ONLY: exp_radius
14 : USE atomic_kind_types, ONLY: atomic_kind_type,&
15 : get_atomic_kind
16 : USE basis_set_types, ONLY: gto_basis_set_type
17 : USE bibliography, ONLY: Kolafa2004,&
18 : Kuhne2007,&
19 : cite_reference
20 : USE cp_control_types, ONLY: dft_control_type
21 : USE cp_log_handling, ONLY: cp_get_default_logger,&
22 : cp_logger_type,&
23 : cp_to_string
24 : USE dbcsr_api, ONLY: &
25 : dbcsr_add, dbcsr_binary_read, dbcsr_checksum, dbcsr_complete_redistribute, &
26 : dbcsr_copy_into_existing, dbcsr_create, dbcsr_desymmetrize, dbcsr_distribution_get, &
27 : dbcsr_distribution_new, dbcsr_distribution_type, dbcsr_dot, dbcsr_filter, &
28 : dbcsr_get_block_p, dbcsr_get_info, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
29 : dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_p_type, &
30 : dbcsr_release, dbcsr_reserve_diag_blocks, dbcsr_scale, dbcsr_set, dbcsr_type
31 : USE dm_ls_scf_methods, ONLY: density_matrix_trs4,&
32 : ls_scf_init_matrix_S
33 : USE dm_ls_scf_qs, ONLY: ls_scf_dm_to_ks,&
34 : ls_scf_qs_atomic_guess,&
35 : matrix_ls_to_qs,&
36 : matrix_qs_to_ls
37 : USE dm_ls_scf_types, ONLY: ls_mstruct_type,&
38 : ls_scf_env_type
39 : USE iterate_matrix, ONLY: purify_mcweeny
40 : USE kinds, ONLY: default_path_length,&
41 : dp
42 : USE machine, ONLY: m_walltime
43 : USE mathlib, ONLY: binomial,&
44 : diamat_all
45 : USE message_passing, ONLY: mp_para_env_type
46 : USE pao_ml, ONLY: pao_ml_forces
47 : USE pao_param, ONLY: pao_calc_AB,&
48 : pao_param_count
49 : USE pao_types, ONLY: pao_env_type
50 : USE particle_types, ONLY: particle_type
51 : USE qs_energy_types, ONLY: qs_energy_type
52 : USE qs_environment_types, ONLY: get_qs_env,&
53 : qs_environment_type
54 : USE qs_initial_guess, ONLY: calculate_atomic_fock_matrix
55 : USE qs_kind_types, ONLY: get_qs_kind,&
56 : pao_descriptor_type,&
57 : pao_potential_type,&
58 : qs_kind_type,&
59 : set_qs_kind
60 : USE qs_ks_methods, ONLY: qs_ks_update_qs_env
61 : USE qs_ks_types, ONLY: qs_ks_did_change
62 : USE qs_rho_methods, ONLY: qs_rho_update_rho
63 : USE qs_rho_types, ONLY: qs_rho_get,&
64 : qs_rho_type
65 :
66 : !$ USE OMP_LIB, ONLY: omp_get_level
67 : #include "./base/base_uses.f90"
68 :
69 : IMPLICIT NONE
70 :
71 : PRIVATE
72 :
73 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pao_methods'
74 :
75 : PUBLIC :: pao_print_atom_info, pao_init_kinds
76 : PUBLIC :: pao_build_orthogonalizer, pao_build_selector
77 : PUBLIC :: pao_build_diag_distribution
78 : PUBLIC :: pao_build_matrix_X, pao_build_core_hamiltonian
79 : PUBLIC :: pao_test_convergence
80 : PUBLIC :: pao_calc_energy, pao_check_trace_ps
81 : PUBLIC :: pao_store_P, pao_add_forces, pao_guess_initial_P
82 : PUBLIC :: pao_check_grad
83 :
84 : CONTAINS
85 :
86 : ! **************************************************************************************************
87 : !> \brief Initialize qs kinds
88 : !> \param pao ...
89 : !> \param qs_env ...
90 : ! **************************************************************************************************
91 94 : SUBROUTINE pao_init_kinds(pao, qs_env)
92 : TYPE(pao_env_type), POINTER :: pao
93 : TYPE(qs_environment_type), POINTER :: qs_env
94 :
95 : CHARACTER(len=*), PARAMETER :: routineN = 'pao_init_kinds'
96 :
97 : INTEGER :: handle, i, ikind, pao_basis_size
98 : TYPE(gto_basis_set_type), POINTER :: basis_set
99 94 : TYPE(pao_descriptor_type), DIMENSION(:), POINTER :: pao_descriptors
100 94 : TYPE(pao_potential_type), DIMENSION(:), POINTER :: pao_potentials
101 94 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
102 :
103 94 : CALL timeset(routineN, handle)
104 94 : CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set)
105 :
106 224 : DO ikind = 1, SIZE(qs_kind_set)
107 : CALL get_qs_kind(qs_kind_set(ikind), &
108 : basis_set=basis_set, &
109 : pao_basis_size=pao_basis_size, &
110 : pao_potentials=pao_potentials, &
111 130 : pao_descriptors=pao_descriptors)
112 :
113 130 : IF (pao_basis_size < 1) THEN
114 : ! pao disabled for ikind, set pao_basis_size to size of primary basis
115 12 : CALL set_qs_kind(qs_kind_set(ikind), pao_basis_size=basis_set%nsgf)
116 : END IF
117 :
118 : ! initialize radii of Gaussians to speedup screeing later on
119 192 : DO i = 1, SIZE(pao_potentials)
120 192 : pao_potentials(i)%beta_radius = exp_radius(0, pao_potentials(i)%beta, pao%eps_pgf, 1.0_dp)
121 : END DO
122 372 : DO i = 1, SIZE(pao_descriptors)
123 18 : pao_descriptors(i)%beta_radius = exp_radius(0, pao_descriptors(i)%beta, pao%eps_pgf, 1.0_dp)
124 148 : pao_descriptors(i)%screening_radius = exp_radius(0, pao_descriptors(i)%screening, pao%eps_pgf, 1.0_dp)
125 : END DO
126 : END DO
127 94 : CALL timestop(handle)
128 94 : END SUBROUTINE pao_init_kinds
129 :
130 : ! **************************************************************************************************
131 : !> \brief Prints a one line summary for each atom.
132 : !> \param pao ...
133 : ! **************************************************************************************************
134 94 : SUBROUTINE pao_print_atom_info(pao)
135 : TYPE(pao_env_type), POINTER :: pao
136 :
137 : INTEGER :: iatom, natoms
138 94 : INTEGER, DIMENSION(:), POINTER :: pao_basis, param_cols, param_rows, &
139 94 : pri_basis
140 :
141 94 : CALL dbcsr_get_info(pao%matrix_Y, row_blk_size=pri_basis, col_blk_size=pao_basis)
142 94 : CPASSERT(SIZE(pao_basis) == SIZE(pri_basis))
143 94 : natoms = SIZE(pao_basis)
144 :
145 94 : CALL dbcsr_get_info(pao%matrix_X, row_blk_size=param_rows, col_blk_size=param_cols)
146 94 : CPASSERT(SIZE(param_rows) == natoms .AND. SIZE(param_cols) == natoms)
147 :
148 94 : IF (pao%iw_atoms > 0) THEN
149 12 : DO iatom = 1, natoms
150 : WRITE (pao%iw_atoms, "(A,I7,T20,A,I3,T45,A,I3,T65,A,I3)") &
151 9 : " PAO| atom: ", iatom, &
152 9 : " prim_basis: ", pri_basis(iatom), &
153 9 : " pao_basis: ", pao_basis(iatom), &
154 21 : " pao_params: ", (param_cols(iatom)*param_rows(iatom))
155 : END DO
156 : END IF
157 94 : END SUBROUTINE pao_print_atom_info
158 :
159 : ! **************************************************************************************************
160 : !> \brief Constructs matrix_N and its inverse.
161 : !> \param pao ...
162 : !> \param qs_env ...
163 : ! **************************************************************************************************
164 94 : SUBROUTINE pao_build_orthogonalizer(pao, qs_env)
165 : TYPE(pao_env_type), POINTER :: pao
166 : TYPE(qs_environment_type), POINTER :: qs_env
167 :
168 : CHARACTER(len=*), PARAMETER :: routineN = 'pao_build_orthogonalizer'
169 :
170 : INTEGER :: acol, arow, handle, i, iatom, j, k, N
171 : LOGICAL :: found
172 : REAL(dp) :: v, w
173 94 : REAL(dp), DIMENSION(:), POINTER :: evals
174 94 : REAL(dp), DIMENSION(:, :), POINTER :: A, block_N, block_N_inv, block_S
175 : TYPE(dbcsr_iterator_type) :: iter
176 94 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
177 :
178 94 : CALL timeset(routineN, handle)
179 :
180 94 : CALL get_qs_env(qs_env, matrix_s=matrix_s)
181 :
182 94 : CALL dbcsr_create(pao%matrix_N, template=matrix_s(1)%matrix, name="PAO matrix_N")
183 94 : CALL dbcsr_reserve_diag_blocks(pao%matrix_N)
184 :
185 94 : CALL dbcsr_create(pao%matrix_N_inv, template=matrix_s(1)%matrix, name="PAO matrix_N_inv")
186 94 : CALL dbcsr_reserve_diag_blocks(pao%matrix_N_inv)
187 :
188 : !$OMP PARALLEL DEFAULT(NONE) SHARED(pao,matrix_s) &
189 94 : !$OMP PRIVATE(iter,arow,acol,iatom,block_N,block_N_inv,block_S,found,N,A,evals,k,i,j,w,v)
190 : CALL dbcsr_iterator_start(iter, pao%matrix_N)
191 : DO WHILE (dbcsr_iterator_blocks_left(iter))
192 : CALL dbcsr_iterator_next_block(iter, arow, acol, block_N)
193 : iatom = arow; CPASSERT(arow == acol)
194 :
195 : CALL dbcsr_get_block_p(matrix=pao%matrix_N_inv, row=iatom, col=iatom, block=block_N_inv, found=found)
196 : CPASSERT(ASSOCIATED(block_N_inv))
197 :
198 : CALL dbcsr_get_block_p(matrix=matrix_s(1)%matrix, row=iatom, col=iatom, block=block_S, found=found)
199 : CPASSERT(ASSOCIATED(block_S))
200 :
201 : N = SIZE(block_S, 1); CPASSERT(SIZE(block_S, 1) == SIZE(block_S, 2)) ! primary basis size
202 : ALLOCATE (A(N, N), evals(N))
203 :
204 : ! take square root of atomic overlap matrix
205 : A = block_S
206 : CALL diamat_all(A, evals) !afterwards A contains the eigenvectors
207 : block_N = 0.0_dp
208 : block_N_inv = 0.0_dp
209 : DO k = 1, N
210 : ! NOTE: To maintain a consistent notation with the Berghold paper,
211 : ! the "_inv" is swapped: N^{-1}=sqrt(S); N=sqrt(S)^{-1}
212 : w = 1.0_dp/SQRT(evals(k))
213 : v = SQRT(evals(k))
214 : DO i = 1, N
215 : DO j = 1, N
216 : block_N(i, j) = block_N(i, j) + w*A(i, k)*A(j, k)
217 : block_N_inv(i, j) = block_N_inv(i, j) + v*A(i, k)*A(j, k)
218 : END DO
219 : END DO
220 : END DO
221 : DEALLOCATE (A, evals)
222 : END DO
223 : CALL dbcsr_iterator_stop(iter)
224 : !$OMP END PARALLEL
225 :
226 : ! store a copies of N and N_inv that are distributed according to pao%diag_distribution
227 : CALL dbcsr_create(pao%matrix_N_diag, &
228 : name="PAO matrix_N_diag", &
229 : dist=pao%diag_distribution, &
230 94 : template=matrix_s(1)%matrix)
231 94 : CALL dbcsr_reserve_diag_blocks(pao%matrix_N_diag)
232 94 : CALL dbcsr_complete_redistribute(pao%matrix_N, pao%matrix_N_diag)
233 : CALL dbcsr_create(pao%matrix_N_inv_diag, &
234 : name="PAO matrix_N_inv_diag", &
235 : dist=pao%diag_distribution, &
236 94 : template=matrix_s(1)%matrix)
237 94 : CALL dbcsr_reserve_diag_blocks(pao%matrix_N_inv_diag)
238 94 : CALL dbcsr_complete_redistribute(pao%matrix_N_inv, pao%matrix_N_inv_diag)
239 :
240 94 : CALL timestop(handle)
241 94 : END SUBROUTINE pao_build_orthogonalizer
242 :
243 : ! **************************************************************************************************
244 : !> \brief Build rectangular matrix to converert between primary and PAO basis.
245 : !> \param pao ...
246 : !> \param qs_env ...
247 : ! **************************************************************************************************
248 94 : SUBROUTINE pao_build_selector(pao, qs_env)
249 : TYPE(pao_env_type), POINTER :: pao
250 : TYPE(qs_environment_type), POINTER :: qs_env
251 :
252 : CHARACTER(len=*), PARAMETER :: routineN = 'pao_build_selector'
253 :
254 : INTEGER :: acol, arow, handle, i, iatom, ikind, M, &
255 : natoms
256 94 : INTEGER, DIMENSION(:), POINTER :: blk_sizes_aux, blk_sizes_pri
257 94 : REAL(dp), DIMENSION(:, :), POINTER :: block_Y
258 : TYPE(dbcsr_iterator_type) :: iter
259 94 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
260 94 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
261 94 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
262 :
263 94 : CALL timeset(routineN, handle)
264 :
265 : CALL get_qs_env(qs_env, &
266 : natom=natoms, &
267 : matrix_s=matrix_s, &
268 : qs_kind_set=qs_kind_set, &
269 94 : particle_set=particle_set)
270 :
271 94 : CALL dbcsr_get_info(matrix_s(1)%matrix, col_blk_size=blk_sizes_pri)
272 :
273 282 : ALLOCATE (blk_sizes_aux(natoms))
274 308 : DO iatom = 1, natoms
275 214 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
276 214 : CALL get_qs_kind(qs_kind_set(ikind), pao_basis_size=M)
277 214 : CPASSERT(M > 0)
278 214 : IF (blk_sizes_pri(iatom) < M) &
279 0 : CPABORT("PAO basis size exceeds primary basis size.")
280 522 : blk_sizes_aux(iatom) = M
281 : END DO
282 :
283 : CALL dbcsr_create(pao%matrix_Y, &
284 : template=matrix_s(1)%matrix, &
285 : matrix_type="N", &
286 : row_blk_size=blk_sizes_pri, &
287 : col_blk_size=blk_sizes_aux, &
288 94 : name="PAO matrix_Y")
289 94 : DEALLOCATE (blk_sizes_aux)
290 :
291 94 : CALL dbcsr_reserve_diag_blocks(pao%matrix_Y)
292 :
293 : !$OMP PARALLEL DEFAULT(NONE) SHARED(pao) &
294 94 : !$OMP PRIVATE(iter,arow,acol,block_Y,i,M)
295 : CALL dbcsr_iterator_start(iter, pao%matrix_Y)
296 : DO WHILE (dbcsr_iterator_blocks_left(iter))
297 : CALL dbcsr_iterator_next_block(iter, arow, acol, block_Y)
298 : M = SIZE(block_Y, 2) ! size of pao basis
299 : block_Y = 0.0_dp
300 : DO i = 1, M
301 : block_Y(i, i) = 1.0_dp
302 : END DO
303 : END DO
304 : CALL dbcsr_iterator_stop(iter)
305 : !$OMP END PARALLEL
306 :
307 94 : CALL timestop(handle)
308 94 : END SUBROUTINE pao_build_selector
309 :
310 : ! **************************************************************************************************
311 : !> \brief Creates new DBCSR distribution which spreads diagonal blocks evenly across ranks
312 : !> \param pao ...
313 : !> \param qs_env ...
314 : ! **************************************************************************************************
315 94 : SUBROUTINE pao_build_diag_distribution(pao, qs_env)
316 : TYPE(pao_env_type), POINTER :: pao
317 : TYPE(qs_environment_type), POINTER :: qs_env
318 :
319 : CHARACTER(len=*), PARAMETER :: routineN = 'pao_build_diag_distribution'
320 :
321 : INTEGER :: handle, iatom, natoms, pgrid_cols, &
322 : pgrid_rows
323 94 : INTEGER, DIMENSION(:), POINTER :: diag_col_dist, diag_row_dist
324 : TYPE(dbcsr_distribution_type) :: main_dist
325 94 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
326 :
327 94 : CALL timeset(routineN, handle)
328 :
329 94 : CALL get_qs_env(qs_env, natom=natoms, matrix_s=matrix_s)
330 :
331 : ! get processor grid from matrix_s
332 94 : CALL dbcsr_get_info(matrix=matrix_s(1)%matrix, distribution=main_dist)
333 94 : CALL dbcsr_distribution_get(main_dist, nprows=pgrid_rows, npcols=pgrid_cols)
334 :
335 : ! create new mapping of matrix-grid to processor-grid
336 470 : ALLOCATE (diag_row_dist(natoms), diag_col_dist(natoms))
337 308 : DO iatom = 1, natoms
338 214 : diag_row_dist(iatom) = MOD(iatom - 1, pgrid_rows)
339 308 : diag_col_dist(iatom) = MOD((iatom - 1)/pgrid_rows, pgrid_cols)
340 : END DO
341 :
342 : ! instanciate distribution object
343 : CALL dbcsr_distribution_new(pao%diag_distribution, template=main_dist, &
344 94 : row_dist=diag_row_dist, col_dist=diag_col_dist)
345 :
346 94 : DEALLOCATE (diag_row_dist, diag_col_dist)
347 :
348 94 : CALL timestop(handle)
349 188 : END SUBROUTINE pao_build_diag_distribution
350 :
351 : ! **************************************************************************************************
352 : !> \brief Creates the matrix_X
353 : !> \param pao ...
354 : !> \param qs_env ...
355 : ! **************************************************************************************************
356 94 : SUBROUTINE pao_build_matrix_X(pao, qs_env)
357 : TYPE(pao_env_type), POINTER :: pao
358 : TYPE(qs_environment_type), POINTER :: qs_env
359 :
360 : CHARACTER(len=*), PARAMETER :: routineN = 'pao_build_matrix_X'
361 :
362 : INTEGER :: handle, iatom, ikind, natoms
363 94 : INTEGER, DIMENSION(:), POINTER :: col_blk_size, row_blk_size
364 94 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
365 :
366 94 : CALL timeset(routineN, handle)
367 :
368 : CALL get_qs_env(qs_env, &
369 : natom=natoms, &
370 94 : particle_set=particle_set)
371 :
372 : ! determine block-sizes of matrix_X
373 470 : ALLOCATE (row_blk_size(natoms), col_blk_size(natoms))
374 308 : col_blk_size = 1
375 308 : DO iatom = 1, natoms
376 214 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
377 308 : CALL pao_param_count(pao, qs_env, ikind, nparams=row_blk_size(iatom))
378 : END DO
379 :
380 : ! build actual matrix_X
381 : CALL dbcsr_create(pao%matrix_X, &
382 : name="PAO matrix_X", &
383 : dist=pao%diag_distribution, &
384 : matrix_type="N", &
385 : row_blk_size=row_blk_size, &
386 94 : col_blk_size=col_blk_size)
387 94 : DEALLOCATE (row_blk_size, col_blk_size)
388 :
389 94 : CALL dbcsr_reserve_diag_blocks(pao%matrix_X)
390 94 : CALL dbcsr_set(pao%matrix_X, 0.0_dp)
391 :
392 94 : CALL timestop(handle)
393 94 : END SUBROUTINE pao_build_matrix_X
394 :
395 : ! **************************************************************************************************
396 : !> \brief Creates the matrix_H0 which contains the core hamiltonian
397 : !> \param pao ...
398 : !> \param qs_env ...
399 : ! **************************************************************************************************
400 94 : SUBROUTINE pao_build_core_hamiltonian(pao, qs_env)
401 : TYPE(pao_env_type), POINTER :: pao
402 : TYPE(qs_environment_type), POINTER :: qs_env
403 :
404 94 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
405 94 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
406 94 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
407 :
408 : CALL get_qs_env(qs_env, &
409 : matrix_s=matrix_s, &
410 : atomic_kind_set=atomic_kind_set, &
411 94 : qs_kind_set=qs_kind_set)
412 :
413 : ! allocate matrix_H0
414 : CALL dbcsr_create(pao%matrix_H0, &
415 : name="PAO matrix_H0", &
416 : dist=pao%diag_distribution, &
417 94 : template=matrix_s(1)%matrix)
418 94 : CALL dbcsr_reserve_diag_blocks(pao%matrix_H0)
419 :
420 : ! calculate initial atomic fock matrix H0
421 : ! Can't use matrix_ks from ls_scf_qs_atomic_guess(), because it's not rotationally invariant.
422 : ! getting H0 directly from the atomic code
423 : CALL calculate_atomic_fock_matrix(pao%matrix_H0, &
424 : atomic_kind_set, &
425 : qs_kind_set, &
426 94 : ounit=pao%iw)
427 :
428 94 : END SUBROUTINE pao_build_core_hamiltonian
429 :
430 : ! **************************************************************************************************
431 : !> \brief Test whether the PAO optimization has reached convergence
432 : !> \param pao ...
433 : !> \param ls_scf_env ...
434 : !> \param new_energy ...
435 : !> \param is_converged ...
436 : ! **************************************************************************************************
437 2622 : SUBROUTINE pao_test_convergence(pao, ls_scf_env, new_energy, is_converged)
438 : TYPE(pao_env_type), POINTER :: pao
439 : TYPE(ls_scf_env_type) :: ls_scf_env
440 : REAL(KIND=dp), INTENT(IN) :: new_energy
441 : LOGICAL, INTENT(OUT) :: is_converged
442 :
443 : REAL(KIND=dp) :: energy_diff, loop_eps, now, time_diff
444 :
445 : ! calculate progress
446 2622 : energy_diff = new_energy - pao%energy_prev
447 2622 : pao%energy_prev = new_energy
448 2622 : now = m_walltime()
449 2622 : time_diff = now - pao%step_start_time
450 2622 : pao%step_start_time = now
451 :
452 : ! convergence criterion
453 2622 : loop_eps = pao%norm_G/ls_scf_env%nelectron_total
454 2622 : is_converged = loop_eps < pao%eps_pao
455 :
456 2622 : IF (pao%istep > 1) THEN
457 2546 : IF (pao%iw > 0) WRITE (pao%iw, *) "PAO| energy improvement:", energy_diff
458 : ! IF(energy_diff>0.0_dp) CPWARN("PAO| energy increased")
459 :
460 : ! print one-liner
461 2546 : IF (pao%iw > 0) WRITE (pao%iw, '(A,I6,11X,F20.9,1X,E10.3,1X,E10.3,1X,F9.3)') &
462 1273 : " PAO| step ", &
463 1273 : pao%istep, &
464 1273 : new_energy, &
465 1273 : loop_eps, &
466 1273 : pao%linesearch%step_size, & !prev step, which let to the current energy
467 2546 : time_diff
468 : END IF
469 2622 : END SUBROUTINE pao_test_convergence
470 :
471 : ! **************************************************************************************************
472 : !> \brief Calculate the pao energy
473 : !> \param pao ...
474 : !> \param qs_env ...
475 : !> \param ls_scf_env ...
476 : !> \param energy ...
477 : ! **************************************************************************************************
478 11818 : SUBROUTINE pao_calc_energy(pao, qs_env, ls_scf_env, energy)
479 : TYPE(pao_env_type), POINTER :: pao
480 : TYPE(qs_environment_type), POINTER :: qs_env
481 : TYPE(ls_scf_env_type), TARGET :: ls_scf_env
482 : REAL(KIND=dp), INTENT(OUT) :: energy
483 :
484 : CHARACTER(len=*), PARAMETER :: routineN = 'pao_calc_energy'
485 :
486 : INTEGER :: handle, ispin
487 : REAL(KIND=dp) :: penalty, trace_PH
488 :
489 11818 : CALL timeset(routineN, handle)
490 :
491 : ! calculate matrix U, which determines the pao basis
492 11818 : CALL pao_calc_AB(pao, qs_env, ls_scf_env, gradient=.FALSE., penalty=penalty)
493 :
494 : ! calculat S, S_inv, S_sqrt, and S_sqrt_inv in the new pao basis
495 11818 : CALL pao_rebuild_S(qs_env, ls_scf_env)
496 :
497 : ! calculate the density matrix P in the pao basis
498 11818 : CALL pao_dm_trs4(qs_env, ls_scf_env)
499 :
500 : ! calculate the energy from the trace(PH) in the pao basis
501 11818 : energy = 0.0_dp
502 23636 : DO ispin = 1, ls_scf_env%nspins
503 11818 : CALL dbcsr_dot(ls_scf_env%matrix_p(ispin), ls_scf_env%matrix_ks(ispin), trace_PH)
504 23636 : energy = energy + trace_PH
505 : END DO
506 :
507 : ! add penalty term
508 11818 : energy = energy + penalty
509 :
510 11818 : IF (pao%iw > 0) THEN
511 5909 : WRITE (pao%iw, *) ""
512 5909 : WRITE (pao%iw, *) "PAO| energy:", energy, "penalty:", penalty
513 : END IF
514 11818 : CALL timestop(handle)
515 11818 : END SUBROUTINE pao_calc_energy
516 :
517 : ! **************************************************************************************************
518 : !> \brief Ensure that the number of electrons is correct.
519 : !> \param ls_scf_env ...
520 : ! **************************************************************************************************
521 10326 : SUBROUTINE pao_check_trace_PS(ls_scf_env)
522 : TYPE(ls_scf_env_type) :: ls_scf_env
523 :
524 : CHARACTER(len=*), PARAMETER :: routineN = 'pao_check_trace_PS'
525 :
526 : INTEGER :: handle, ispin
527 : REAL(KIND=dp) :: tmp, trace_PS
528 : TYPE(dbcsr_type) :: matrix_S_desym
529 :
530 10326 : CALL timeset(routineN, handle)
531 10326 : CALL dbcsr_create(matrix_S_desym, template=ls_scf_env%matrix_s, matrix_type="N")
532 10326 : CALL dbcsr_desymmetrize(ls_scf_env%matrix_s, matrix_S_desym)
533 :
534 10326 : trace_PS = 0.0_dp
535 20652 : DO ispin = 1, ls_scf_env%nspins
536 10326 : CALL dbcsr_dot(ls_scf_env%matrix_p(ispin), matrix_S_desym, tmp)
537 20652 : trace_PS = trace_PS + tmp
538 : END DO
539 :
540 10326 : CALL dbcsr_release(matrix_S_desym)
541 :
542 10326 : IF (ABS(ls_scf_env%nelectron_total - trace_PS) > 0.5) &
543 0 : CPABORT("Number of electrons wrong. Trace(PS) ="//cp_to_string(trace_PS))
544 :
545 10326 : CALL timestop(handle)
546 10326 : END SUBROUTINE pao_check_trace_PS
547 :
548 : ! **************************************************************************************************
549 : !> \brief Read primary density matrix from file.
550 : !> \param pao ...
551 : !> \param qs_env ...
552 : ! **************************************************************************************************
553 28 : SUBROUTINE pao_read_preopt_dm(pao, qs_env)
554 : TYPE(pao_env_type), POINTER :: pao
555 : TYPE(qs_environment_type), POINTER :: qs_env
556 :
557 : CHARACTER(len=*), PARAMETER :: routineN = 'pao_read_preopt_dm'
558 :
559 : INTEGER :: handle, ispin
560 : REAL(KIND=dp) :: cs_pos
561 : TYPE(dbcsr_distribution_type) :: dist
562 28 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, rho_ao
563 : TYPE(dbcsr_type) :: matrix_tmp
564 : TYPE(dft_control_type), POINTER :: dft_control
565 : TYPE(qs_energy_type), POINTER :: energy
566 : TYPE(qs_rho_type), POINTER :: rho
567 :
568 28 : CALL timeset(routineN, handle)
569 :
570 : CALL get_qs_env(qs_env, &
571 : dft_control=dft_control, &
572 : matrix_s=matrix_s, &
573 : rho=rho, &
574 28 : energy=energy)
575 :
576 28 : CALL qs_rho_get(rho, rho_ao=rho_ao)
577 :
578 28 : IF (dft_control%nspins /= 1) CPABORT("open shell not yet implemented")
579 :
580 28 : CALL dbcsr_get_info(matrix_s(1)%matrix, distribution=dist)
581 :
582 56 : DO ispin = 1, dft_control%nspins
583 28 : CALL dbcsr_binary_read(pao%preopt_dm_file, matrix_new=matrix_tmp, distribution=dist)
584 28 : cs_pos = dbcsr_checksum(matrix_tmp, pos=.TRUE.)
585 28 : IF (pao%iw > 0) WRITE (pao%iw, *) "PAO| Read restart DM "// &
586 14 : TRIM(pao%preopt_dm_file)//" with checksum: ", cs_pos
587 28 : CALL dbcsr_copy_into_existing(rho_ao(ispin)%matrix, matrix_tmp)
588 56 : CALL dbcsr_release(matrix_tmp)
589 : END DO
590 :
591 : ! calculate corresponding ks matrix
592 28 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
593 28 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
594 : CALL qs_ks_update_qs_env(qs_env, calculate_forces=.FALSE., &
595 28 : just_energy=.FALSE., print_active=.TRUE.)
596 28 : IF (pao%iw > 0) WRITE (pao%iw, *) "PAO| Quickstep energy from restart density:", energy%total
597 :
598 28 : CALL timestop(handle)
599 :
600 28 : END SUBROUTINE pao_read_preopt_dm
601 :
602 : ! **************************************************************************************************
603 : !> \brief Rebuilds S, S_inv, S_sqrt, and S_sqrt_inv in the pao basis
604 : !> \param qs_env ...
605 : !> \param ls_scf_env ...
606 : ! **************************************************************************************************
607 11818 : SUBROUTINE pao_rebuild_S(qs_env, ls_scf_env)
608 : TYPE(qs_environment_type), POINTER :: qs_env
609 : TYPE(ls_scf_env_type), TARGET :: ls_scf_env
610 :
611 : CHARACTER(len=*), PARAMETER :: routineN = 'pao_rebuild_S'
612 :
613 : INTEGER :: handle
614 11818 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
615 :
616 11818 : CALL timeset(routineN, handle)
617 :
618 11818 : CALL dbcsr_release(ls_scf_env%matrix_s_inv)
619 11818 : CALL dbcsr_release(ls_scf_env%matrix_s_sqrt)
620 11818 : CALL dbcsr_release(ls_scf_env%matrix_s_sqrt_inv)
621 :
622 11818 : CALL get_qs_env(qs_env, matrix_s=matrix_s)
623 11818 : CALL ls_scf_init_matrix_s(matrix_s(1)%matrix, ls_scf_env)
624 :
625 11818 : CALL timestop(handle)
626 11818 : END SUBROUTINE pao_rebuild_S
627 :
628 : ! **************************************************************************************************
629 : !> \brief Calculate density matrix using TRS4 purification
630 : !> \param qs_env ...
631 : !> \param ls_scf_env ...
632 : ! **************************************************************************************************
633 11818 : SUBROUTINE pao_dm_trs4(qs_env, ls_scf_env)
634 : TYPE(qs_environment_type), POINTER :: qs_env
635 : TYPE(ls_scf_env_type), TARGET :: ls_scf_env
636 :
637 : CHARACTER(len=*), PARAMETER :: routineN = 'pao_dm_trs4'
638 :
639 : CHARACTER(LEN=default_path_length) :: project_name
640 : INTEGER :: handle, ispin, nelectron_spin_real, nspin
641 : LOGICAL :: converged
642 : REAL(KIND=dp) :: homo_spin, lumo_spin, mu_spin
643 : TYPE(cp_logger_type), POINTER :: logger
644 11818 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
645 :
646 11818 : CALL timeset(routineN, handle)
647 11818 : logger => cp_get_default_logger()
648 11818 : project_name = logger%iter_info%project_name
649 11818 : nspin = ls_scf_env%nspins
650 :
651 11818 : CALL get_qs_env(qs_env, matrix_ks=matrix_ks)
652 23636 : DO ispin = 1, nspin
653 : CALL matrix_qs_to_ls(ls_scf_env%matrix_ks(ispin), matrix_ks(ispin)%matrix, &
654 11818 : ls_scf_env%ls_mstruct, covariant=.TRUE.)
655 :
656 11818 : nelectron_spin_real = ls_scf_env%nelectron_spin(ispin)
657 11818 : IF (ls_scf_env%nspins == 1) nelectron_spin_real = nelectron_spin_real/2
658 : CALL density_matrix_trs4(ls_scf_env%matrix_p(ispin), ls_scf_env%matrix_ks(ispin), &
659 : ls_scf_env%matrix_s_sqrt_inv, &
660 : nelectron_spin_real, ls_scf_env%eps_filter, homo_spin, lumo_spin, mu_spin, &
661 : dynamic_threshold=.FALSE., converged=converged, &
662 : max_iter_lanczos=ls_scf_env%max_iter_lanczos, &
663 11818 : eps_lanczos=ls_scf_env%eps_lanczos)
664 23636 : IF (.NOT. converged) CPABORT("TRS4 did not converge")
665 : END DO
666 :
667 11818 : IF (nspin == 1) CALL dbcsr_scale(ls_scf_env%matrix_p(1), 2.0_dp)
668 :
669 11818 : CALL timestop(handle)
670 11818 : END SUBROUTINE pao_dm_trs4
671 :
672 : ! **************************************************************************************************
673 : !> \brief Debugging routine for checking the analytic gradient.
674 : !> \param pao ...
675 : !> \param qs_env ...
676 : !> \param ls_scf_env ...
677 : ! **************************************************************************************************
678 2634 : SUBROUTINE pao_check_grad(pao, qs_env, ls_scf_env)
679 : TYPE(pao_env_type), POINTER :: pao
680 : TYPE(qs_environment_type), POINTER :: qs_env
681 : TYPE(ls_scf_env_type), TARGET :: ls_scf_env
682 :
683 : CHARACTER(len=*), PARAMETER :: routineN = 'pao_check_grad'
684 :
685 : INTEGER :: handle, i, iatom, j, natoms
686 2622 : INTEGER, DIMENSION(:), POINTER :: blk_sizes_col, blk_sizes_row
687 : LOGICAL :: found
688 : REAL(dp) :: delta, delta_max, eps, Gij_num
689 2622 : REAL(dp), DIMENSION(:, :), POINTER :: block_G, block_X
690 : TYPE(ls_mstruct_type), POINTER :: ls_mstruct
691 : TYPE(mp_para_env_type), POINTER :: para_env
692 :
693 2610 : IF (pao%check_grad_tol < 0.0_dp) RETURN ! no checking
694 :
695 12 : CALL timeset(routineN, handle)
696 :
697 12 : ls_mstruct => ls_scf_env%ls_mstruct
698 :
699 12 : CALL get_qs_env(qs_env, para_env=para_env, natom=natoms)
700 :
701 12 : eps = pao%num_grad_eps
702 12 : delta_max = 0.0_dp
703 :
704 12 : CALL dbcsr_get_info(pao%matrix_X, col_blk_size=blk_sizes_col, row_blk_size=blk_sizes_row)
705 :
706 : ! can not use an iterator here, because other DBCSR routines are called within loop.
707 38 : DO iatom = 1, natoms
708 26 : IF (pao%iw > 0) WRITE (pao%iw, *) 'PAO| checking gradient of atom ', iatom
709 26 : CALL dbcsr_get_block_p(matrix=pao%matrix_X, row=iatom, col=iatom, block=block_X, found=found)
710 :
711 26 : IF (ASSOCIATED(block_X)) THEN !only one node actually has the block
712 13 : CALL dbcsr_get_block_p(matrix=pao%matrix_G, row=iatom, col=iatom, block=block_G, found=found)
713 13 : CPASSERT(ASSOCIATED(block_G))
714 : END IF
715 :
716 586 : DO i = 1, blk_sizes_row(iatom)
717 1070 : DO j = 1, blk_sizes_col(iatom)
718 828 : SELECT CASE (pao%num_grad_order)
719 : CASE (2) ! calculate derivative to 2th order
720 306 : Gij_num = -eval_point(block_X, i, j, -eps, pao, ls_scf_env, qs_env)
721 306 : Gij_num = Gij_num + eval_point(block_X, i, j, +eps, pao, ls_scf_env, qs_env)
722 306 : Gij_num = Gij_num/(2.0_dp*eps)
723 :
724 : CASE (4) ! calculate derivative to 4th order
725 180 : Gij_num = eval_point(block_X, i, j, -2_dp*eps, pao, ls_scf_env, qs_env)
726 180 : Gij_num = Gij_num - 8_dp*eval_point(block_X, i, j, -1_dp*eps, pao, ls_scf_env, qs_env)
727 180 : Gij_num = Gij_num + 8_dp*eval_point(block_X, i, j, +1_dp*eps, pao, ls_scf_env, qs_env)
728 180 : Gij_num = Gij_num - eval_point(block_X, i, j, +2_dp*eps, pao, ls_scf_env, qs_env)
729 180 : Gij_num = Gij_num/(12.0_dp*eps)
730 :
731 : CASE (6) ! calculate derivative to 6th order
732 36 : Gij_num = -1_dp*eval_point(block_X, i, j, -3_dp*eps, pao, ls_scf_env, qs_env)
733 36 : Gij_num = Gij_num + 9_dp*eval_point(block_X, i, j, -2_dp*eps, pao, ls_scf_env, qs_env)
734 36 : Gij_num = Gij_num - 45_dp*eval_point(block_X, i, j, -1_dp*eps, pao, ls_scf_env, qs_env)
735 36 : Gij_num = Gij_num + 45_dp*eval_point(block_X, i, j, +1_dp*eps, pao, ls_scf_env, qs_env)
736 36 : Gij_num = Gij_num - 9_dp*eval_point(block_X, i, j, +2_dp*eps, pao, ls_scf_env, qs_env)
737 36 : Gij_num = Gij_num + 1_dp*eval_point(block_X, i, j, +3_dp*eps, pao, ls_scf_env, qs_env)
738 36 : Gij_num = Gij_num/(60.0_dp*eps)
739 :
740 : CASE DEFAULT
741 522 : CPABORT("Unsupported numerical derivative order: "//cp_to_string(pao%num_grad_order))
742 : END SELECT
743 :
744 1044 : IF (ASSOCIATED(block_X)) THEN
745 261 : delta = ABS(Gij_num - block_G(i, j))
746 261 : delta_max = MAX(delta_max, delta)
747 : !WRITE (*,*) "gradient check", iatom, i, j, Gij_num, block_G(i,j), delta
748 : END IF
749 : END DO
750 : END DO
751 : END DO
752 :
753 12 : CALL para_env%max(delta_max)
754 12 : IF (pao%iw > 0) WRITE (pao%iw, *) 'PAO| checked gradient, max delta:', delta_max
755 12 : IF (delta_max > pao%check_grad_tol) CALL cp_abort(__LOCATION__, &
756 0 : "Analytic and numeric gradients differ too much:"//cp_to_string(delta_max))
757 :
758 12 : CALL timestop(handle)
759 2622 : END SUBROUTINE pao_check_grad
760 :
761 : ! **************************************************************************************************
762 : !> \brief Helper routine for pao_check_grad()
763 : !> \param block_X ...
764 : !> \param i ...
765 : !> \param j ...
766 : !> \param eps ...
767 : !> \param pao ...
768 : !> \param ls_scf_env ...
769 : !> \param qs_env ...
770 : !> \return ...
771 : ! **************************************************************************************************
772 3096 : FUNCTION eval_point(block_X, i, j, eps, pao, ls_scf_env, qs_env) RESULT(energy)
773 : REAL(dp), DIMENSION(:, :), POINTER :: block_X
774 : INTEGER, INTENT(IN) :: i, j
775 : REAL(dp), INTENT(IN) :: eps
776 : TYPE(pao_env_type), POINTER :: pao
777 : TYPE(ls_scf_env_type), TARGET :: ls_scf_env
778 : TYPE(qs_environment_type), POINTER :: qs_env
779 : REAL(dp) :: energy
780 :
781 : REAL(dp) :: old_Xij
782 :
783 1548 : IF (ASSOCIATED(block_X)) THEN
784 774 : old_Xij = block_X(i, j) ! backup old block_X
785 774 : block_X(i, j) = block_X(i, j) + eps ! add perturbation
786 : END IF
787 :
788 : ! calculate energy
789 1548 : CALL pao_calc_energy(pao, qs_env, ls_scf_env, energy)
790 :
791 : ! restore old block_X
792 1548 : IF (ASSOCIATED(block_X)) THEN
793 774 : block_X(i, j) = old_Xij
794 : END IF
795 :
796 1548 : END FUNCTION eval_point
797 :
798 : ! **************************************************************************************************
799 : !> \brief Stores density matrix as initial guess for next SCF optimization.
800 : !> \param qs_env ...
801 : !> \param ls_scf_env ...
802 : ! **************************************************************************************************
803 556 : SUBROUTINE pao_store_P(qs_env, ls_scf_env)
804 : TYPE(qs_environment_type), POINTER :: qs_env
805 : TYPE(ls_scf_env_type), TARGET :: ls_scf_env
806 :
807 : CHARACTER(len=*), PARAMETER :: routineN = 'pao_store_P'
808 :
809 : INTEGER :: handle, ispin, istore
810 278 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
811 : TYPE(dft_control_type), POINTER :: dft_control
812 : TYPE(ls_mstruct_type), POINTER :: ls_mstruct
813 : TYPE(pao_env_type), POINTER :: pao
814 :
815 0 : IF (ls_scf_env%scf_history%nstore == 0) RETURN
816 278 : CALL timeset(routineN, handle)
817 278 : ls_mstruct => ls_scf_env%ls_mstruct
818 278 : pao => ls_scf_env%pao_env
819 278 : CALL get_qs_env(qs_env, dft_control=dft_control, matrix_s=matrix_s)
820 :
821 278 : ls_scf_env%scf_history%istore = ls_scf_env%scf_history%istore + 1
822 278 : istore = MOD(ls_scf_env%scf_history%istore - 1, ls_scf_env%scf_history%nstore) + 1
823 278 : IF (pao%iw > 0) WRITE (pao%iw, *) "PAO| Storing density matrix for ASPC guess in slot:", istore
824 :
825 : ! initialize storage
826 278 : IF (ls_scf_env%scf_history%istore <= ls_scf_env%scf_history%nstore) THEN
827 208 : DO ispin = 1, dft_control%nspins
828 208 : CALL dbcsr_create(ls_scf_env%scf_history%matrix(ispin, istore), template=matrix_s(1)%matrix)
829 : END DO
830 : END IF
831 :
832 : ! We are storing the density matrix in the non-orthonormal primary basis.
833 : ! While the orthonormal basis would yield better extrapolations,
834 : ! we simply can not afford to calculat S_sqrt in the primary basis.
835 556 : DO ispin = 1, dft_control%nspins
836 : ! transform into primary basis
837 : CALL matrix_ls_to_qs(ls_scf_env%scf_history%matrix(ispin, istore), ls_scf_env%matrix_p(ispin), &
838 556 : ls_scf_env%ls_mstruct, covariant=.FALSE., keep_sparsity=.FALSE.)
839 : END DO
840 :
841 278 : CALL timestop(handle)
842 278 : END SUBROUTINE pao_store_P
843 :
844 : ! **************************************************************************************************
845 : !> \brief Provide an initial guess for the density matrix
846 : !> \param pao ...
847 : !> \param qs_env ...
848 : !> \param ls_scf_env ...
849 : ! **************************************************************************************************
850 278 : SUBROUTINE pao_guess_initial_P(pao, qs_env, ls_scf_env)
851 : TYPE(pao_env_type), POINTER :: pao
852 : TYPE(qs_environment_type), POINTER :: qs_env
853 : TYPE(ls_scf_env_type), TARGET :: ls_scf_env
854 :
855 : CHARACTER(len=*), PARAMETER :: routineN = 'pao_guess_initial_P'
856 :
857 : INTEGER :: handle
858 :
859 278 : CALL timeset(routineN, handle)
860 :
861 278 : IF (ls_scf_env%scf_history%istore > 0) THEN
862 184 : CALL pao_aspc_guess_P(pao, qs_env, ls_scf_env)
863 184 : pao%need_initial_scf = .TRUE.
864 : ELSE
865 94 : IF (LEN_TRIM(pao%preopt_dm_file) > 0) THEN
866 28 : CALL pao_read_preopt_dm(pao, qs_env)
867 28 : pao%need_initial_scf = .FALSE.
868 28 : pao%preopt_dm_file = "" ! load only for first MD step
869 : ELSE
870 66 : CALL ls_scf_qs_atomic_guess(qs_env, ls_scf_env%energy_init)
871 66 : IF (pao%iw > 0) WRITE (pao%iw, '(A,F20.9)') &
872 33 : " PAO| Energy from initial atomic guess:", ls_scf_env%energy_init
873 66 : pao%need_initial_scf = .TRUE.
874 : END IF
875 : END IF
876 :
877 278 : CALL timestop(handle)
878 :
879 278 : END SUBROUTINE pao_guess_initial_P
880 :
881 : ! **************************************************************************************************
882 : !> \brief Run the Always Stable Predictor-Corrector to guess an initial density matrix
883 : !> \param pao ...
884 : !> \param qs_env ...
885 : !> \param ls_scf_env ...
886 : ! **************************************************************************************************
887 184 : SUBROUTINE pao_aspc_guess_P(pao, qs_env, ls_scf_env)
888 : TYPE(pao_env_type), POINTER :: pao
889 : TYPE(qs_environment_type), POINTER :: qs_env
890 : TYPE(ls_scf_env_type), TARGET :: ls_scf_env
891 :
892 : CHARACTER(len=*), PARAMETER :: routineN = 'pao_aspc_guess_P'
893 :
894 : INTEGER :: handle, iaspc, ispin, istore, naspc
895 : REAL(dp) :: alpha
896 184 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
897 : TYPE(dbcsr_type) :: matrix_P
898 : TYPE(dft_control_type), POINTER :: dft_control
899 : TYPE(ls_mstruct_type), POINTER :: ls_mstruct
900 :
901 184 : CALL timeset(routineN, handle)
902 184 : ls_mstruct => ls_scf_env%ls_mstruct
903 184 : CPASSERT(ls_scf_env%scf_history%istore > 0)
904 184 : CALL cite_reference(Kolafa2004)
905 184 : CALL cite_reference(Kuhne2007)
906 184 : CALL get_qs_env(qs_env, dft_control=dft_control, matrix_s=matrix_s)
907 :
908 184 : IF (pao%iw > 0) WRITE (pao%iw, *) "PAO| Calculating initial guess with ASPC"
909 :
910 184 : CALL dbcsr_create(matrix_P, template=matrix_s(1)%matrix)
911 :
912 184 : naspc = MIN(ls_scf_env%scf_history%istore, ls_scf_env%scf_history%nstore)
913 368 : DO ispin = 1, dft_control%nspins
914 : ! actual extrapolation
915 184 : CALL dbcsr_set(matrix_P, 0.0_dp)
916 392 : DO iaspc = 1, naspc
917 : alpha = (-1.0_dp)**(iaspc + 1)*REAL(iaspc, KIND=dp)* &
918 208 : binomial(2*naspc, naspc - iaspc)/binomial(2*naspc - 2, naspc - 1)
919 208 : istore = MOD(ls_scf_env%scf_history%istore - iaspc, ls_scf_env%scf_history%nstore) + 1
920 392 : CALL dbcsr_add(matrix_P, ls_scf_env%scf_history%matrix(ispin, istore), 1.0_dp, alpha)
921 : END DO
922 :
923 : ! transform back from primary basis into pao basis
924 368 : CALL matrix_qs_to_ls(ls_scf_env%matrix_p(ispin), matrix_P, ls_scf_env%ls_mstruct, covariant=.FALSE.)
925 : END DO
926 :
927 184 : CALL dbcsr_release(matrix_P)
928 :
929 : ! linear combination of P's is not idempotent. A bit of McWeeny is needed to ensure it is again
930 368 : DO ispin = 1, dft_control%nspins
931 184 : IF (dft_control%nspins == 1) CALL dbcsr_scale(ls_scf_env%matrix_p(ispin), 0.5_dp)
932 : ! to ensure that noisy blocks do not build up during MD (in particular with curvy) filter that guess a bit more
933 184 : CALL dbcsr_filter(ls_scf_env%matrix_p(ispin), ls_scf_env%eps_filter**(2.0_dp/3.0_dp))
934 : ! we could go to the orthonomal basis, but it seems not worth the trouble
935 : ! TODO : 10 iterations is a conservative upper bound, figure out when it fails
936 184 : CALL purify_mcweeny(ls_scf_env%matrix_p(ispin:ispin), ls_scf_env%matrix_s, ls_scf_env%eps_filter, 10)
937 368 : IF (dft_control%nspins == 1) CALL dbcsr_scale(ls_scf_env%matrix_p(ispin), 2.0_dp)
938 : END DO
939 :
940 184 : CALL pao_check_trace_PS(ls_scf_env) ! sanity check
941 :
942 : ! compute corresponding energy and ks matrix
943 184 : CALL ls_scf_dm_to_ks(qs_env, ls_scf_env, ls_scf_env%energy_init, iscf=0)
944 :
945 184 : CALL timestop(handle)
946 184 : END SUBROUTINE pao_aspc_guess_P
947 :
948 : ! **************************************************************************************************
949 : !> \brief Calculate the forces contributed by PAO
950 : !> \param qs_env ...
951 : !> \param ls_scf_env ...
952 : ! **************************************************************************************************
953 42 : SUBROUTINE pao_add_forces(qs_env, ls_scf_env)
954 : TYPE(qs_environment_type), POINTER :: qs_env
955 : TYPE(ls_scf_env_type), TARGET :: ls_scf_env
956 :
957 : CHARACTER(len=*), PARAMETER :: routineN = 'pao_add_forces'
958 :
959 : INTEGER :: handle, iatom, natoms
960 42 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: forces
961 : TYPE(mp_para_env_type), POINTER :: para_env
962 : TYPE(pao_env_type), POINTER :: pao
963 42 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
964 :
965 42 : CALL timeset(routineN, handle)
966 42 : pao => ls_scf_env%pao_env
967 :
968 42 : IF (pao%iw > 0) WRITE (pao%iw, *) "PAO| Adding forces."
969 :
970 42 : IF (pao%max_pao /= 0) THEN
971 20 : IF (pao%penalty_strength /= 0.0_dp) &
972 0 : CPABORT("PAO forces require PENALTY_STRENGTH or MAX_PAO set to zero")
973 20 : IF (pao%linpot_regu_strength /= 0.0_dp) &
974 0 : CPABORT("PAO forces require LINPOT_REGULARIZATION_STRENGTH or MAX_PAO set to zero")
975 20 : IF (pao%regularization /= 0.0_dp) &
976 0 : CPABORT("PAO forces require REGULARIZATION or MAX_PAO set to zero")
977 : END IF
978 :
979 : CALL get_qs_env(qs_env, &
980 : para_env=para_env, &
981 : particle_set=particle_set, &
982 42 : natom=natoms)
983 :
984 126 : ALLOCATE (forces(natoms, 3))
985 42 : CALL pao_calc_AB(pao, qs_env, ls_scf_env, gradient=.TRUE., forces=forces) ! without penalty terms
986 :
987 42 : IF (SIZE(pao%ml_training_set) > 0) &
988 18 : CALL pao_ml_forces(pao, qs_env, pao%matrix_G, forces)
989 :
990 42 : CALL para_env%sum(forces)
991 136 : DO iatom = 1, natoms
992 418 : particle_set(iatom)%f = particle_set(iatom)%f + forces(iatom, :)
993 : END DO
994 :
995 42 : DEALLOCATE (forces)
996 :
997 42 : CALL timestop(handle)
998 :
999 42 : END SUBROUTINE pao_add_forces
1000 :
1001 : END MODULE pao_methods
|
