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