LCOV - code coverage report
Current view: top level - src - pao_methods.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 98.0 % 348 341
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 19 19

            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        11802 :    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        11802 :       CALL timeset(routineN, handle)
     505              : 
     506              :       ! calculate matrix U, which determines the pao basis
     507        11802 :       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        11802 :       CALL pao_rebuild_S(qs_env, ls_scf_env)
     511              : 
     512              :       ! calculate the density matrix P in the pao basis
     513        11802 :       CALL pao_dm_trs4(qs_env, ls_scf_env)
     514              : 
     515              :       ! calculate the energy from the trace(PH) in the pao basis
     516        11802 :       energy = 0.0_dp
     517        23604 :       DO ispin = 1, ls_scf_env%nspins
     518        11802 :          CALL dbcsr_dot(ls_scf_env%matrix_p(ispin), ls_scf_env%matrix_ks(ispin), trace_PH)
     519        23604 :          energy = energy + trace_PH
     520              :       END DO
     521              : 
     522              :       ! add penalty term
     523        11802 :       energy = energy + penalty
     524              : 
     525        11802 :       IF (pao%iw > 0) THEN
     526         5901 :          WRITE (pao%iw, *) ""
     527         5901 :          WRITE (pao%iw, *) "PAO| energy:", energy, "penalty:", penalty
     528              :       END IF
     529        11802 :       CALL timestop(handle)
     530        11802 :    END SUBROUTINE pao_calc_energy
     531              : 
     532              : ! **************************************************************************************************
     533              : !> \brief Ensure that the number of electrons is correct.
     534              : !> \param ls_scf_env ...
     535              : ! **************************************************************************************************
     536        10322 :    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        10322 :       CALL timeset(routineN, handle)
     546        10322 :       CALL dbcsr_create(matrix_S_desym, template=ls_scf_env%matrix_s, matrix_type="N")
     547        10322 :       CALL dbcsr_desymmetrize(ls_scf_env%matrix_s, matrix_S_desym)
     548              : 
     549        10322 :       trace_PS = 0.0_dp
     550        20644 :       DO ispin = 1, ls_scf_env%nspins
     551        10322 :          CALL dbcsr_dot(ls_scf_env%matrix_p(ispin), matrix_S_desym, tmp)
     552        20644 :          trace_PS = trace_PS + tmp
     553              :       END DO
     554              : 
     555        10322 :       CALL dbcsr_release(matrix_S_desym)
     556              : 
     557        10322 :       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        10322 :       CALL timestop(handle)
     561        10322 :    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        11802 :    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        11802 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s
     630              : 
     631        11802 :       CALL timeset(routineN, handle)
     632              : 
     633        11802 :       CALL dbcsr_release(ls_scf_env%matrix_s_inv)
     634        11802 :       CALL dbcsr_release(ls_scf_env%matrix_s_sqrt)
     635        11802 :       CALL dbcsr_release(ls_scf_env%matrix_s_sqrt_inv)
     636              : 
     637        11802 :       CALL get_qs_env(qs_env, matrix_s=matrix_s)
     638        11802 :       CALL ls_scf_init_matrix_s(matrix_s(1)%matrix, ls_scf_env)
     639              : 
     640        11802 :       CALL timestop(handle)
     641        11802 :    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        11802 :    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        11802 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_ks
     660              : 
     661        11802 :       CALL timeset(routineN, handle)
     662        11802 :       logger => cp_get_default_logger()
     663        11802 :       project_name = logger%iter_info%project_name
     664        11802 :       nspin = ls_scf_env%nspins
     665              : 
     666        11802 :       CALL get_qs_env(qs_env, matrix_ks=matrix_ks)
     667        23604 :       DO ispin = 1, nspin
     668              :          CALL matrix_qs_to_ls(ls_scf_env%matrix_ks(ispin), matrix_ks(ispin)%matrix, &
     669        11802 :                               ls_scf_env%ls_mstruct, covariant=.TRUE.)
     670              : 
     671        11802 :          nelectron_spin_real = ls_scf_env%nelectron_spin(ispin)
     672        11802 :          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        11802 :                                   eps_lanczos=ls_scf_env%eps_lanczos)
     679        23604 :          IF (.NOT. converged) CPABORT("TRS4 did not converge")
     680              :       END DO
     681              : 
     682        11802 :       IF (nspin == 1) CALL dbcsr_scale(ls_scf_env%matrix_p(1), 2.0_dp)
     683              : 
     684        11802 :       CALL timestop(handle)
     685        11802 :    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
        

Generated by: LCOV version 2.0-1