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