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