(git:58e3e09)
pao_methods.F
Go to the documentation of this file.
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 ! **************************************************************************************************
13  USE ao_util, ONLY: exp_radius
14  USE atomic_kind_types, ONLY: atomic_kind_type,&
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
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
33  USE dm_ls_scf_qs, ONLY: ls_scf_dm_to_ks,&
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,&
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,&
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
55  USE qs_kind_types, ONLY: get_qs_kind,&
56  pao_descriptor_type,&
57  pao_potential_type,&
58  qs_kind_type,&
61  USE qs_ks_types, ONLY: qs_ks_did_change
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 
79  PUBLIC :: pao_test_convergence
82  PUBLIC :: pao_check_grad
83 
84 CONTAINS
85 
86 ! **************************************************************************************************
87 !> \brief Initialize qs kinds
88 !> \param pao ...
89 !> \param qs_env ...
90 ! **************************************************************************************************
91  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  TYPE(pao_descriptor_type), DIMENSION(:), POINTER :: pao_descriptors
100  TYPE(pao_potential_type), DIMENSION(:), POINTER :: pao_potentials
101  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
102 
103  CALL timeset(routinen, handle)
104  CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set)
105 
106  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, &
111  pao_descriptors=pao_descriptors)
112 
113  IF (pao_basis_size < 1) THEN
114  ! pao disabled for ikind, set pao_basis_size to size of primary basis
115  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  DO i = 1, SIZE(pao_potentials)
120  pao_potentials(i)%beta_radius = exp_radius(0, pao_potentials(i)%beta, pao%eps_pgf, 1.0_dp)
121  END DO
122  DO i = 1, SIZE(pao_descriptors)
123  pao_descriptors(i)%beta_radius = exp_radius(0, pao_descriptors(i)%beta, pao%eps_pgf, 1.0_dp)
124  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  CALL timestop(handle)
128  END SUBROUTINE pao_init_kinds
129 
130 ! **************************************************************************************************
131 !> \brief Prints a one line summary for each atom.
132 !> \param pao ...
133 ! **************************************************************************************************
134  SUBROUTINE pao_print_atom_info(pao)
135  TYPE(pao_env_type), POINTER :: pao
136 
137  INTEGER :: iatom, natoms
138  INTEGER, DIMENSION(:), POINTER :: pao_basis, param_cols, param_rows, &
139  pri_basis
140 
141  CALL dbcsr_get_info(pao%matrix_Y, row_blk_size=pri_basis, col_blk_size=pao_basis)
142  cpassert(SIZE(pao_basis) == SIZE(pri_basis))
143  natoms = SIZE(pao_basis)
144 
145  CALL dbcsr_get_info(pao%matrix_X, row_blk_size=param_rows, col_blk_size=param_cols)
146  cpassert(SIZE(param_rows) == natoms .AND. SIZE(param_cols) == natoms)
147 
148  IF (pao%iw_atoms > 0) THEN
149  DO iatom = 1, natoms
150  WRITE (pao%iw_atoms, "(A,I7,T20,A,I3,T45,A,I3,T65,A,I3)") &
151  " PAO| atom: ", iatom, &
152  " prim_basis: ", pri_basis(iatom), &
153  " pao_basis: ", pao_basis(iatom), &
154  " pao_params: ", (param_cols(iatom)*param_rows(iatom))
155  END DO
156  END IF
157  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  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  REAL(dp), DIMENSION(:), POINTER :: evals
174  REAL(dp), DIMENSION(:, :), POINTER :: a, block_n, block_n_inv, block_s
175  TYPE(dbcsr_iterator_type) :: iter
176  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
177 
178  CALL timeset(routinen, handle)
179 
180  CALL get_qs_env(qs_env, matrix_s=matrix_s)
181 
182  CALL dbcsr_create(pao%matrix_N, template=matrix_s(1)%matrix, name="PAO matrix_N")
183  CALL dbcsr_reserve_diag_blocks(pao%matrix_N)
184 
185  CALL dbcsr_create(pao%matrix_N_inv, template=matrix_s(1)%matrix, name="PAO matrix_N_inv")
186  CALL dbcsr_reserve_diag_blocks(pao%matrix_N_inv)
187 
188 !$OMP PARALLEL DEFAULT(NONE) SHARED(pao,matrix_s) &
189 !$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  template=matrix_s(1)%matrix)
231  CALL dbcsr_reserve_diag_blocks(pao%matrix_N_diag)
232  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  template=matrix_s(1)%matrix)
237  CALL dbcsr_reserve_diag_blocks(pao%matrix_N_inv_diag)
238  CALL dbcsr_complete_redistribute(pao%matrix_N_inv, pao%matrix_N_inv_diag)
239 
240  CALL timestop(handle)
241  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  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  INTEGER, DIMENSION(:), POINTER :: blk_sizes_aux, blk_sizes_pri
257  REAL(dp), DIMENSION(:, :), POINTER :: block_y
258  TYPE(dbcsr_iterator_type) :: iter
259  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
260  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
261  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
262 
263  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  particle_set=particle_set)
270 
271  CALL dbcsr_get_info(matrix_s(1)%matrix, col_blk_size=blk_sizes_pri)
272 
273  ALLOCATE (blk_sizes_aux(natoms))
274  DO iatom = 1, natoms
275  CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
276  CALL get_qs_kind(qs_kind_set(ikind), pao_basis_size=m)
277  cpassert(m > 0)
278  IF (blk_sizes_pri(iatom) < m) &
279  cpabort("PAO basis size exceeds primary basis size.")
280  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  name="PAO matrix_Y")
289  DEALLOCATE (blk_sizes_aux)
290 
291  CALL dbcsr_reserve_diag_blocks(pao%matrix_Y)
292 
293 !$OMP PARALLEL DEFAULT(NONE) SHARED(pao) &
294 !$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  CALL timestop(handle)
308  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  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  INTEGER, DIMENSION(:), POINTER :: diag_col_dist, diag_row_dist
324  TYPE(dbcsr_distribution_type) :: main_dist
325  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
326 
327  CALL timeset(routinen, handle)
328 
329  CALL get_qs_env(qs_env, natom=natoms, matrix_s=matrix_s)
330 
331  ! get processor grid from matrix_s
332  CALL dbcsr_get_info(matrix=matrix_s(1)%matrix, distribution=main_dist)
333  CALL dbcsr_distribution_get(main_dist, nprows=pgrid_rows, npcols=pgrid_cols)
334 
335  ! create new mapping of matrix-grid to processor-grid
336  ALLOCATE (diag_row_dist(natoms), diag_col_dist(natoms))
337  DO iatom = 1, natoms
338  diag_row_dist(iatom) = mod(iatom - 1, pgrid_rows)
339  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  row_dist=diag_row_dist, col_dist=diag_col_dist)
345 
346  DEALLOCATE (diag_row_dist, diag_col_dist)
347 
348  CALL timestop(handle)
349  END SUBROUTINE pao_build_diag_distribution
350 
351 ! **************************************************************************************************
352 !> \brief Creates the matrix_X
353 !> \param pao ...
354 !> \param qs_env ...
355 ! **************************************************************************************************
356  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  INTEGER, DIMENSION(:), POINTER :: col_blk_size, row_blk_size
364  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
365 
366  CALL timeset(routinen, handle)
367 
368  CALL get_qs_env(qs_env, &
369  natom=natoms, &
370  particle_set=particle_set)
371 
372  ! determine block-sizes of matrix_X
373  ALLOCATE (row_blk_size(natoms), col_blk_size(natoms))
374  col_blk_size = 1
375  DO iatom = 1, natoms
376  CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
377  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  col_blk_size=col_blk_size)
387  DEALLOCATE (row_blk_size, col_blk_size)
388 
389  CALL dbcsr_reserve_diag_blocks(pao%matrix_X)
390  CALL dbcsr_set(pao%matrix_X, 0.0_dp)
391 
392  CALL timestop(handle)
393  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  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  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
405  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
406  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  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  template=matrix_s(1)%matrix)
418  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  ounit=pao%iw)
427 
428  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  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  energy_diff = new_energy - pao%energy_prev
447  pao%energy_prev = new_energy
448  now = m_walltime()
449  time_diff = now - pao%step_start_time
450  pao%step_start_time = now
451 
452  ! convergence criterion
453  loop_eps = pao%norm_G/ls_scf_env%nelectron_total
454  is_converged = loop_eps < pao%eps_pao
455 
456  IF (pao%istep > 1) THEN
457  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  IF (pao%iw > 0) WRITE (pao%iw, '(A,I6,11X,F20.9,1X,E10.3,1X,E10.3,1X,F9.3)') &
462  " PAO| step ", &
463  pao%istep, &
464  new_energy, &
465  loop_eps, &
466  pao%linesearch%step_size, & !prev step, which let to the current energy
467  time_diff
468  END IF
469  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  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  CALL timeset(routinen, handle)
490 
491  ! calculate matrix U, which determines the pao basis
492  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  CALL pao_rebuild_s(qs_env, ls_scf_env)
496 
497  ! calculate the density matrix P in the pao basis
498  CALL pao_dm_trs4(qs_env, ls_scf_env)
499 
500  ! calculate the energy from the trace(PH) in the pao basis
501  energy = 0.0_dp
502  DO ispin = 1, ls_scf_env%nspins
503  CALL dbcsr_dot(ls_scf_env%matrix_p(ispin), ls_scf_env%matrix_ks(ispin), trace_ph)
504  energy = energy + trace_ph
505  END DO
506 
507  ! add penalty term
508  energy = energy + penalty
509 
510  IF (pao%iw > 0) THEN
511  WRITE (pao%iw, *) ""
512  WRITE (pao%iw, *) "PAO| energy:", energy, "penalty:", penalty
513  END IF
514  CALL timestop(handle)
515  END SUBROUTINE pao_calc_energy
516 
517 ! **************************************************************************************************
518 !> \brief Ensure that the number of electrons is correct.
519 !> \param ls_scf_env ...
520 ! **************************************************************************************************
521  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  CALL timeset(routinen, handle)
531  CALL dbcsr_create(matrix_s_desym, template=ls_scf_env%matrix_s, matrix_type="N")
532  CALL dbcsr_desymmetrize(ls_scf_env%matrix_s, matrix_s_desym)
533 
534  trace_ps = 0.0_dp
535  DO ispin = 1, ls_scf_env%nspins
536  CALL dbcsr_dot(ls_scf_env%matrix_p(ispin), matrix_s_desym, tmp)
537  trace_ps = trace_ps + tmp
538  END DO
539 
540  CALL dbcsr_release(matrix_s_desym)
541 
542  IF (abs(ls_scf_env%nelectron_total - trace_ps) > 0.5) &
543  cpabort("Number of electrons wrong. Trace(PS) ="//cp_to_string(trace_ps))
544 
545  CALL timestop(handle)
546  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  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  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  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  energy=energy)
575 
576  CALL qs_rho_get(rho, rho_ao=rho_ao)
577 
578  IF (dft_control%nspins /= 1) cpabort("open shell not yet implemented")
579 
580  CALL dbcsr_get_info(matrix_s(1)%matrix, distribution=dist)
581 
582  DO ispin = 1, dft_control%nspins
583  CALL dbcsr_binary_read(pao%preopt_dm_file, matrix_new=matrix_tmp, distribution=dist)
584  cs_pos = dbcsr_checksum(matrix_tmp, pos=.true.)
585  IF (pao%iw > 0) WRITE (pao%iw, *) "PAO| Read restart DM "// &
586  trim(pao%preopt_dm_file)//" with checksum: ", cs_pos
587  CALL dbcsr_copy_into_existing(rho_ao(ispin)%matrix, matrix_tmp)
588  CALL dbcsr_release(matrix_tmp)
589  END DO
590 
591  ! calculate corresponding ks matrix
592  CALL qs_rho_update_rho(rho, qs_env=qs_env)
593  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  just_energy=.false., print_active=.true.)
596  IF (pao%iw > 0) WRITE (pao%iw, *) "PAO| Quickstep energy from restart density:", energy%total
597 
598  CALL timestop(handle)
599 
600  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  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  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
615 
616  CALL timeset(routinen, handle)
617 
618  CALL dbcsr_release(ls_scf_env%matrix_s_inv)
619  CALL dbcsr_release(ls_scf_env%matrix_s_sqrt)
620  CALL dbcsr_release(ls_scf_env%matrix_s_sqrt_inv)
621 
622  CALL get_qs_env(qs_env, matrix_s=matrix_s)
623  CALL ls_scf_init_matrix_s(matrix_s(1)%matrix, ls_scf_env)
624 
625  CALL timestop(handle)
626  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  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  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
645 
646  CALL timeset(routinen, handle)
647  logger => cp_get_default_logger()
648  project_name = logger%iter_info%project_name
649  nspin = ls_scf_env%nspins
650 
651  CALL get_qs_env(qs_env, matrix_ks=matrix_ks)
652  DO ispin = 1, nspin
653  CALL matrix_qs_to_ls(ls_scf_env%matrix_ks(ispin), matrix_ks(ispin)%matrix, &
654  ls_scf_env%ls_mstruct, covariant=.true.)
655 
656  nelectron_spin_real = ls_scf_env%nelectron_spin(ispin)
657  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  eps_lanczos=ls_scf_env%eps_lanczos)
664  IF (.NOT. converged) cpabort("TRS4 did not converge")
665  END DO
666 
667  IF (nspin == 1) CALL dbcsr_scale(ls_scf_env%matrix_p(1), 2.0_dp)
668 
669  CALL timestop(handle)
670  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  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  INTEGER, DIMENSION(:), POINTER :: blk_sizes_col, blk_sizes_row
687  LOGICAL :: found
688  REAL(dp) :: delta, delta_max, eps, gij_num
689  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  IF (pao%check_grad_tol < 0.0_dp) RETURN ! no checking
694 
695  CALL timeset(routinen, handle)
696 
697  ls_mstruct => ls_scf_env%ls_mstruct
698 
699  CALL get_qs_env(qs_env, para_env=para_env, natom=natoms)
700 
701  eps = pao%num_grad_eps
702  delta_max = 0.0_dp
703 
704  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  DO iatom = 1, natoms
708  IF (pao%iw > 0) WRITE (pao%iw, *) 'PAO| checking gradient of atom ', iatom
709  CALL dbcsr_get_block_p(matrix=pao%matrix_X, row=iatom, col=iatom, block=block_x, found=found)
710 
711  IF (ASSOCIATED(block_x)) THEN !only one node actually has the block
712  CALL dbcsr_get_block_p(matrix=pao%matrix_G, row=iatom, col=iatom, block=block_g, found=found)
713  cpassert(ASSOCIATED(block_g))
714  END IF
715 
716  DO i = 1, blk_sizes_row(iatom)
717  DO j = 1, blk_sizes_col(iatom)
718  SELECT CASE (pao%num_grad_order)
719  CASE (2) ! calculate derivative to 2th order
720  gij_num = -eval_point(block_x, i, j, -eps, pao, ls_scf_env, qs_env)
721  gij_num = gij_num + eval_point(block_x, i, j, +eps, pao, ls_scf_env, qs_env)
722  gij_num = gij_num/(2.0_dp*eps)
723 
724  CASE (4) ! calculate derivative to 4th order
725  gij_num = eval_point(block_x, i, j, -2_dp*eps, pao, ls_scf_env, qs_env)
726  gij_num = gij_num - 8_dp*eval_point(block_x, i, j, -1_dp*eps, pao, ls_scf_env, qs_env)
727  gij_num = gij_num + 8_dp*eval_point(block_x, i, j, +1_dp*eps, pao, ls_scf_env, qs_env)
728  gij_num = gij_num - eval_point(block_x, i, j, +2_dp*eps, pao, ls_scf_env, qs_env)
729  gij_num = gij_num/(12.0_dp*eps)
730 
731  CASE (6) ! calculate derivative to 6th order
732  gij_num = -1_dp*eval_point(block_x, i, j, -3_dp*eps, pao, ls_scf_env, qs_env)
733  gij_num = gij_num + 9_dp*eval_point(block_x, i, j, -2_dp*eps, pao, ls_scf_env, qs_env)
734  gij_num = gij_num - 45_dp*eval_point(block_x, i, j, -1_dp*eps, pao, ls_scf_env, qs_env)
735  gij_num = gij_num + 45_dp*eval_point(block_x, i, j, +1_dp*eps, pao, ls_scf_env, qs_env)
736  gij_num = gij_num - 9_dp*eval_point(block_x, i, j, +2_dp*eps, pao, ls_scf_env, qs_env)
737  gij_num = gij_num + 1_dp*eval_point(block_x, i, j, +3_dp*eps, pao, ls_scf_env, qs_env)
738  gij_num = gij_num/(60.0_dp*eps)
739 
740  CASE DEFAULT
741  cpabort("Unsupported numerical derivative order: "//cp_to_string(pao%num_grad_order))
742  END SELECT
743 
744  IF (ASSOCIATED(block_x)) THEN
745  delta = abs(gij_num - block_g(i, j))
746  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  CALL para_env%max(delta_max)
754  IF (pao%iw > 0) WRITE (pao%iw, *) 'PAO| checked gradient, max delta:', delta_max
755  IF (delta_max > pao%check_grad_tol) CALL cp_abort(__location__, &
756  "Analytic and numeric gradients differ too much:"//cp_to_string(delta_max))
757 
758  CALL timestop(handle)
759  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  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  IF (ASSOCIATED(block_x)) THEN
784  old_xij = block_x(i, j) ! backup old block_X
785  block_x(i, j) = block_x(i, j) + eps ! add perturbation
786  END IF
787 
788  ! calculate energy
789  CALL pao_calc_energy(pao, qs_env, ls_scf_env, energy)
790 
791  ! restore old block_X
792  IF (ASSOCIATED(block_x)) THEN
793  block_x(i, j) = old_xij
794  END IF
795 
796  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  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  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  IF (ls_scf_env%scf_history%nstore == 0) RETURN
816  CALL timeset(routinen, handle)
817  ls_mstruct => ls_scf_env%ls_mstruct
818  pao => ls_scf_env%pao_env
819  CALL get_qs_env(qs_env, dft_control=dft_control, matrix_s=matrix_s)
820 
821  ls_scf_env%scf_history%istore = ls_scf_env%scf_history%istore + 1
822  istore = mod(ls_scf_env%scf_history%istore - 1, ls_scf_env%scf_history%nstore) + 1
823  IF (pao%iw > 0) WRITE (pao%iw, *) "PAO| Storing density matrix for ASPC guess in slot:", istore
824 
825  ! initialize storage
826  IF (ls_scf_env%scf_history%istore <= ls_scf_env%scf_history%nstore) THEN
827  DO ispin = 1, dft_control%nspins
828  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  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  ls_scf_env%ls_mstruct, covariant=.false., keep_sparsity=.false.)
839  END DO
840 
841  CALL timestop(handle)
842  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  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  CALL timeset(routinen, handle)
860 
861  IF (ls_scf_env%scf_history%istore > 0) THEN
862  CALL pao_aspc_guess_p(pao, qs_env, ls_scf_env)
863  pao%need_initial_scf = .true.
864  ELSE
865  IF (len_trim(pao%preopt_dm_file) > 0) THEN
866  CALL pao_read_preopt_dm(pao, qs_env)
867  pao%need_initial_scf = .false.
868  pao%preopt_dm_file = "" ! load only for first MD step
869  ELSE
870  CALL ls_scf_qs_atomic_guess(qs_env, ls_scf_env%energy_init)
871  IF (pao%iw > 0) WRITE (pao%iw, '(A,F20.9)') &
872  " PAO| Energy from initial atomic guess:", ls_scf_env%energy_init
873  pao%need_initial_scf = .true.
874  END IF
875  END IF
876 
877  CALL timestop(handle)
878 
879  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  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  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  CALL timeset(routinen, handle)
902  ls_mstruct => ls_scf_env%ls_mstruct
903  cpassert(ls_scf_env%scf_history%istore > 0)
904  CALL cite_reference(kolafa2004)
905  CALL cite_reference(kuhne2007)
906  CALL get_qs_env(qs_env, dft_control=dft_control, matrix_s=matrix_s)
907 
908  IF (pao%iw > 0) WRITE (pao%iw, *) "PAO| Calculating initial guess with ASPC"
909 
910  CALL dbcsr_create(matrix_p, template=matrix_s(1)%matrix)
911 
912  naspc = min(ls_scf_env%scf_history%istore, ls_scf_env%scf_history%nstore)
913  DO ispin = 1, dft_control%nspins
914  ! actual extrapolation
915  CALL dbcsr_set(matrix_p, 0.0_dp)
916  DO iaspc = 1, naspc
917  alpha = (-1.0_dp)**(iaspc + 1)*real(iaspc, kind=dp)* &
918  binomial(2*naspc, naspc - iaspc)/binomial(2*naspc - 2, naspc - 1)
919  istore = mod(ls_scf_env%scf_history%istore - iaspc, ls_scf_env%scf_history%nstore) + 1
920  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  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  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  DO ispin = 1, dft_control%nspins
931  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  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  CALL purify_mcweeny(ls_scf_env%matrix_p(ispin:ispin), ls_scf_env%matrix_s, ls_scf_env%eps_filter, 10)
937  IF (dft_control%nspins == 1) CALL dbcsr_scale(ls_scf_env%matrix_p(ispin), 2.0_dp)
938  END DO
939 
940  CALL pao_check_trace_ps(ls_scf_env) ! sanity check
941 
942  ! compute corresponding energy and ks matrix
943  CALL ls_scf_dm_to_ks(qs_env, ls_scf_env, ls_scf_env%energy_init, iscf=0)
944 
945  CALL timestop(handle)
946  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  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  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: forces
961  TYPE(mp_para_env_type), POINTER :: para_env
962  TYPE(pao_env_type), POINTER :: pao
963  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
964 
965  CALL timeset(routinen, handle)
966  pao => ls_scf_env%pao_env
967 
968  IF (pao%iw > 0) WRITE (pao%iw, *) "PAO| Adding forces."
969 
970  IF (pao%max_pao /= 0) THEN
971  IF (pao%penalty_strength /= 0.0_dp) &
972  cpabort("PAO forces require PENALTY_STRENGTH or MAX_PAO set to zero")
973  IF (pao%linpot_regu_strength /= 0.0_dp) &
974  cpabort("PAO forces require LINPOT_REGULARIZATION_STRENGTH or MAX_PAO set to zero")
975  IF (pao%regularization /= 0.0_dp) &
976  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  natom=natoms)
983 
984  ALLOCATE (forces(natoms, 3))
985  CALL pao_calc_ab(pao, qs_env, ls_scf_env, gradient=.true., forces=forces) ! without penalty terms
986 
987  IF (SIZE(pao%ml_training_set) > 0) &
988  CALL pao_ml_forces(pao, qs_env, pao%matrix_G, forces)
989 
990  CALL para_env%sum(forces)
991  DO iatom = 1, natoms
992  particle_set(iatom)%f = particle_set(iatom)%f + forces(iatom, :)
993  END DO
994 
995  DEALLOCATE (forces)
996 
997  CALL timestop(handle)
998 
999  END SUBROUTINE pao_add_forces
1000 
1001 END MODULE pao_methods
All kind of helpful little routines.
Definition: ao_util.F:14
real(kind=dp) function, public exp_radius(l, alpha, threshold, prefactor, epsabs, epsrel, rlow)
The radius of a primitive Gaussian function for a given threshold is calculated. g(r) = prefactor*r**...
Definition: ao_util.F:96
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition: bibliography.F:28
integer, save, public kuhne2007
Definition: bibliography.F:43
integer, save, public kolafa2004
Definition: bibliography.F:43
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
various routines to log and control the output. The idea is that decisions about where to log should ...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
lower level routines for linear scaling SCF
subroutine, public ls_scf_init_matrix_s(matrix_s, ls_scf_env)
initialize S matrix related properties (sqrt, inverse...) Might be factored-out since this seems comm...
subroutine, public density_matrix_trs4(matrix_p, matrix_ks, matrix_s_sqrt_inv, nelectron, threshold, e_homo, e_lumo, e_mu, dynamic_threshold, matrix_ks_deviation, max_iter_lanczos, eps_lanczos, converged)
compute the density matrix using a trace-resetting algorithm
Routines for a linear scaling quickstep SCF run based on the density matrix, with a focus on the inte...
Definition: dm_ls_scf_qs.F:15
subroutine, public matrix_ls_to_qs(matrix_qs, matrix_ls, ls_mstruct, covariant, keep_sparsity)
second link to QS, copy a LS matrix to QS matrix used to isolate QS style matrices from LS style will...
Definition: dm_ls_scf_qs.F:307
subroutine, public ls_scf_dm_to_ks(qs_env, ls_scf_env, energy_new, iscf)
use the density matrix in ls_scf_env to compute the new energy and KS matrix
Definition: dm_ls_scf_qs.F:546
subroutine, public matrix_qs_to_ls(matrix_ls, matrix_qs, ls_mstruct, covariant)
first link to QS, copy a QS matrix to LS matrix used to isolate QS style matrices from LS style will ...
Definition: dm_ls_scf_qs.F:213
subroutine, public ls_scf_qs_atomic_guess(qs_env, energy)
get an atomic initial guess
Definition: dm_ls_scf_qs.F:465
Types needed for a linear scaling quickstep SCF run based on the density matrix.
Routines useful for iterative matrix calculations.
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
integer, parameter, public default_path_length
Definition: kinds.F:58
Machine interface based on Fortran 2003 and POSIX.
Definition: machine.F:17
real(kind=dp) function, public m_walltime()
returns time from a real-time clock, protected against rolling early/easily
Definition: machine.F:123
Collection of simple mathematical functions and subroutines.
Definition: mathlib.F:15
elemental real(kind=dp) function, public binomial(n, k)
The binomial coefficient n over k for 0 <= k <= n is calculated, otherwise zero is returned.
Definition: mathlib.F:206
subroutine, public diamat_all(a, eigval, dac)
Diagonalize the symmetric n by n matrix a using the LAPACK library. Only the upper triangle of matrix...
Definition: mathlib.F:376
Interface to the message passing library MPI.
Methods used by pao_main.F.
Definition: pao_methods.F:12
subroutine, public pao_build_orthogonalizer(pao, qs_env)
Constructs matrix_N and its inverse.
Definition: pao_methods.F:165
subroutine, public pao_build_core_hamiltonian(pao, qs_env)
Creates the matrix_H0 which contains the core hamiltonian.
Definition: pao_methods.F:401
subroutine, public pao_guess_initial_p(pao, qs_env, ls_scf_env)
Provide an initial guess for the density matrix.
Definition: pao_methods.F:851
subroutine, public pao_test_convergence(pao, ls_scf_env, new_energy, is_converged)
Test whether the PAO optimization has reached convergence.
Definition: pao_methods.F:438
subroutine, public pao_add_forces(qs_env, ls_scf_env)
Calculate the forces contributed by PAO.
Definition: pao_methods.F:954
subroutine, public pao_check_trace_ps(ls_scf_env)
Ensure that the number of electrons is correct.
Definition: pao_methods.F:522
subroutine, public pao_init_kinds(pao, qs_env)
Initialize qs kinds.
Definition: pao_methods.F:92
subroutine, public pao_build_matrix_x(pao, qs_env)
Creates the matrix_X.
Definition: pao_methods.F:357
subroutine, public pao_check_grad(pao, qs_env, ls_scf_env)
Debugging routine for checking the analytic gradient.
Definition: pao_methods.F:679
subroutine, public pao_print_atom_info(pao)
Prints a one line summary for each atom.
Definition: pao_methods.F:135
subroutine, public pao_store_p(qs_env, ls_scf_env)
Stores density matrix as initial guess for next SCF optimization.
Definition: pao_methods.F:804
subroutine, public pao_build_selector(pao, qs_env)
Build rectangular matrix to converert between primary and PAO basis.
Definition: pao_methods.F:249
subroutine, public pao_calc_energy(pao, qs_env, ls_scf_env, energy)
Calculate the pao energy.
Definition: pao_methods.F:479
subroutine, public pao_build_diag_distribution(pao, qs_env)
Creates new DBCSR distribution which spreads diagonal blocks evenly across ranks.
Definition: pao_methods.F:316
Main module for PAO Machine Learning.
Definition: pao_ml.F:12
subroutine, public pao_ml_forces(pao, qs_env, matrix_G, forces)
Calculate forces contributed by machine learning.
Definition: pao_ml.F:622
Front-End for any PAO parametrization.
Definition: pao_param.F:12
subroutine, public pao_calc_ab(pao, qs_env, ls_scf_env, gradient, penalty, forces)
Takes current matrix_X and calculates the matrices A and B.
Definition: pao_param.F:70
subroutine, public pao_param_count(pao, qs_env, ikind, nparams)
Returns the number of parameters for given atomic kind.
Definition: pao_param.F:173
Factory routines for potentials used e.g. by pao_param_exp and pao_ml.
Types used by the PAO machinery.
Definition: pao_types.F:12
Define the data structure for the particle information.
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Routines to somehow generate an initial guess.
subroutine, public calculate_atomic_fock_matrix(matrix_f, atomic_kind_set, qs_kind_set, ounit)
returns a block diagonal fock matrix.
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_r3d_rs_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, U_of_dft_plus_u, J_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, J0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
subroutine, public set_qs_kind(qs_kind, paw_atom, ghost, floating, hard_radius, hard0_radius, covalent_radius, vdw_radius, lmax_rho0, zeff, no_optimize, dispersion, u_minus_j, reltmat, dftb_parameter, xtb_parameter, elec_conf, pao_basis_size)
Set the components of an atomic kind data set.
routines that build the Kohn-Sham matrix (i.e calculate the coulomb and xc parts
Definition: qs_ks_methods.F:22
subroutine, public qs_ks_update_qs_env(qs_env, calculate_forces, just_energy, print_active)
updates the Kohn Sham matrix of the given qs_env (facility method)
subroutine, public qs_ks_did_change(ks_env, s_mstruct_changed, rho_changed, potential_changed, full_reset)
tells that some of the things relevant to the ks calculation did change. has to be called when change...
Definition: qs_ks_types.F:919
methods of the rho structure (defined in qs_rho_types)
subroutine, public qs_rho_update_rho(rho_struct, qs_env, rho_xc_external, local_rho_set, task_list_external, task_list_external_soft, pw_env_external, para_env_external)
updates rho_r and rho_g to the rhorho_ao. if use_kinetic_energy_density also computes tau_r and tau_g...
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition: qs_rho_types.F:18
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition: qs_rho_types.F:229