1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Calculates integral matrices for RIGPW method 8!> \par History 9!> created JGH [08.2012] 10!> Dorothea Golze [02.2014] (1) extended, re-structured, cleaned 11!> (2) heavily debugged 12!> updated for RI JGH [08.2017] 13!> \authors JGH 14!> Dorothea Golze 15! ************************************************************************************************** 16MODULE ri_environment_methods 17 USE arnoldi_api, ONLY: arnoldi_conjugate_gradient 18 USE atomic_kind_types, ONLY: atomic_kind_type,& 19 get_atomic_kind_set 20 USE cp_blacs_env, ONLY: cp_blacs_env_type 21 USE cp_control_types, ONLY: dft_control_type 22 USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl 23 USE cp_dbcsr_diag, ONLY: cp_dbcsr_syevd 24 USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set 25 USE cp_para_types, ONLY: cp_para_env_type 26 USE dbcsr_api, ONLY: & 27 dbcsr_add, dbcsr_create, dbcsr_dot, dbcsr_get_info, dbcsr_iterator_blocks_left, & 28 dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, & 29 dbcsr_multiply, dbcsr_p_type, dbcsr_release, dbcsr_scale_by_vector, dbcsr_set, dbcsr_type, & 30 dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, dbcsr_type_symmetric 31 USE iterate_matrix, ONLY: invert_Hotelling 32 USE kinds, ONLY: dp 33 USE lri_environment_types, ONLY: allocate_lri_coefs,& 34 lri_density_create,& 35 lri_density_release,& 36 lri_density_type,& 37 lri_environment_type,& 38 lri_kind_type 39 USE message_passing, ONLY: mp_sum 40 USE pw_types, ONLY: pw_p_type 41 USE qs_collocate_density, ONLY: calculate_lri_rho_elec 42 USE qs_environment_types, ONLY: get_qs_env,& 43 qs_environment_type,& 44 set_qs_env 45 USE qs_ks_types, ONLY: qs_ks_env_type 46 USE qs_o3c_methods, ONLY: calculate_o3c_integrals,& 47 contract12_o3c 48 USE qs_o3c_types, ONLY: get_o3c_iterator_info,& 49 init_o3c_container,& 50 o3c_container_type,& 51 o3c_iterate,& 52 o3c_iterator_create,& 53 o3c_iterator_release,& 54 o3c_iterator_type,& 55 release_o3c_container 56 USE qs_overlap, ONLY: build_overlap_matrix 57 USE qs_rho_types, ONLY: qs_rho_get,& 58 qs_rho_type 59 USE util, ONLY: get_limit 60 61!$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num 62#include "./base/base_uses.f90" 63 64 IMPLICIT NONE 65 66 PRIVATE 67 68! ************************************************************************************************** 69 70 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ri_environment_methods' 71 72 PUBLIC :: build_ri_matrices, calculate_ri_densities, ri_metric_solver 73 74! ************************************************************************************************** 75 76CONTAINS 77 78! ************************************************************************************************** 79!> \brief creates and initializes an lri_env 80!> \param lri_env the lri_environment you want to create 81!> \param qs_env ... 82!> \param calculate_forces ... 83! ************************************************************************************************** 84 SUBROUTINE build_ri_matrices(lri_env, qs_env, calculate_forces) 85 86 TYPE(lri_environment_type), POINTER :: lri_env 87 TYPE(qs_environment_type), POINTER :: qs_env 88 LOGICAL, INTENT(IN) :: calculate_forces 89 90 CHARACTER(LEN=*), PARAMETER :: routineN = 'build_ri_matrices', & 91 routineP = moduleN//':'//routineN 92 93 CALL calculate_ri_integrals(lri_env, qs_env, calculate_forces) 94 95 END SUBROUTINE build_ri_matrices 96 97! ************************************************************************************************** 98!> \brief calculates integrals needed for the RI density fitting, 99!> integrals are calculated once, before the SCF starts 100!> \param lri_env ... 101!> \param qs_env ... 102!> \param calculate_forces ... 103! ************************************************************************************************** 104 SUBROUTINE calculate_ri_integrals(lri_env, qs_env, calculate_forces) 105 106 TYPE(lri_environment_type), POINTER :: lri_env 107 TYPE(qs_environment_type), POINTER :: qs_env 108 LOGICAL, INTENT(IN) :: calculate_forces 109 110 CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_ri_integrals', & 111 routineP = moduleN//':'//routineN 112 REAL(KIND=dp), DIMENSION(2), PARAMETER :: fx = (/0.0_dp, 1.0_dp/) 113 114 INTEGER :: handle, i, i1, i2, iatom, ispin, izero, & 115 j, j1, j2, jatom, m, n, nbas 116 INTEGER, DIMENSION(:, :), POINTER :: bas_ptr 117 REAL(KIND=dp) :: fpre 118 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eval 119 REAL(KIND=dp), DIMENSION(:, :), POINTER :: avec, fblk, fout 120 TYPE(cp_blacs_env_type), POINTER :: blacs_env 121 TYPE(cp_para_env_type), POINTER :: para_env 122 TYPE(dbcsr_iterator_type) :: iter 123 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p 124 TYPE(dbcsr_type) :: emat 125 TYPE(dbcsr_type), POINTER :: fmat 126 TYPE(dft_control_type), POINTER :: dft_control 127 TYPE(qs_ks_env_type), POINTER :: ks_env 128 TYPE(qs_rho_type), POINTER :: rho 129 130 CALL timeset(routineN, handle) 131 132 CPASSERT(ASSOCIATED(lri_env)) 133 CPASSERT(ASSOCIATED(qs_env)) 134 135 ! overlap matrices 136 CALL get_qs_env(qs_env=qs_env, ks_env=ks_env) 137 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control) 138 NULLIFY (rho, matrix_p) 139 ! orbital overlap matrix (needed for N constraints and forces) 140 ! we replicate this in order to not directly interact qith qs_env 141 IF (calculate_forces) THEN 142 ! charge constraint forces are calculated here 143 CALL get_qs_env(qs_env=qs_env, rho=rho) 144 CALL qs_rho_get(rho, rho_ao=matrix_p) 145 ALLOCATE (fmat) 146 CALL dbcsr_create(fmat, template=matrix_p(1)%matrix) 147 DO ispin = 1, dft_control%nspins 148 fpre = lri_env%ri_fit%ftrm1n(ispin)/lri_env%ri_fit%ntrm1n 149 CALL dbcsr_add(fmat, matrix_p(ispin)%matrix, fx(ispin), -fpre) 150 END DO 151 CALL build_overlap_matrix(ks_env=ks_env, matrix_s=lri_env%ob_smat, nderivative=0, & 152 basis_type_a="ORB", basis_type_b="ORB", sab_nl=lri_env%soo_list, & 153 calculate_forces=.TRUE., matrix_p=fmat) 154 CALL dbcsr_release(fmat) 155 DEALLOCATE (fmat) 156 ELSE 157 CALL build_overlap_matrix(ks_env=ks_env, matrix_s=lri_env%ob_smat, nderivative=0, & 158 basis_type_a="ORB", basis_type_b="ORB", sab_nl=lri_env%soo_list) 159 END IF 160 ! RI overlap matrix 161 CALL build_overlap_matrix(ks_env=ks_env, matrix_s=lri_env%ri_smat, nderivative=0, & 162 basis_type_a="RI_HXC", basis_type_b="RI_HXC", sab_nl=lri_env%saa_list) 163 IF (calculate_forces) THEN 164 ! calculate f*a(T) pseudo density matrix for forces 165 bas_ptr => lri_env%ri_fit%bas_ptr 166 avec => lri_env%ri_fit%avec 167 fout => lri_env%ri_fit%fout 168 ALLOCATE (fmat) 169 CALL dbcsr_create(fmat, template=lri_env%ri_smat(1)%matrix) 170 CALL cp_dbcsr_alloc_block_from_nbl(fmat, lri_env%saa_list) 171 CALL dbcsr_set(fmat, 0.0_dp) 172 CALL dbcsr_iterator_start(iter, fmat) 173 DO WHILE (dbcsr_iterator_blocks_left(iter)) 174 CALL dbcsr_iterator_next_block(iter, iatom, jatom, fblk) 175 i1 = bas_ptr(1, iatom) 176 i2 = bas_ptr(2, iatom) 177 j1 = bas_ptr(1, jatom) 178 j2 = bas_ptr(2, jatom) 179 IF (iatom <= jatom) THEN 180 DO ispin = 1, dft_control%nspins 181 DO j = j1, j2 182 m = j - j1 + 1 183 DO i = i1, i2 184 n = i - i1 + 1 185 fblk(n, m) = fblk(n, m) + fout(i, ispin)*avec(j, ispin) + avec(i, ispin)*fout(j, ispin) 186 END DO 187 END DO 188 END DO 189 ELSE 190 DO ispin = 1, dft_control%nspins 191 DO i = i1, i2 192 n = i - i1 + 1 193 DO j = j1, j2 194 m = j - j1 + 1 195 fblk(m, n) = fblk(m, n) + fout(i, ispin)*avec(j, ispin) + avec(i, ispin)*fout(j, ispin) 196 END DO 197 END DO 198 END DO 199 END IF 200 IF (iatom == jatom) THEN 201 fblk(:, :) = 0.25_dp*fblk(:, :) 202 ELSE 203 fblk(:, :) = 0.5_dp*fblk(:, :) 204 END IF 205 END DO 206 CALL dbcsr_iterator_stop(iter) 207 ! 208 CALL build_overlap_matrix(ks_env=ks_env, matrix_s=lri_env%ri_smat, nderivative=0, & 209 basis_type_a="RI_HXC", basis_type_b="RI_HXC", sab_nl=lri_env%saa_list, & 210 calculate_forces=.TRUE., matrix_p=fmat) 211 CALL dbcsr_release(fmat) 212 DEALLOCATE (fmat) 213 END IF 214 ! approximation (preconditioner) or exact inverse of RI overlap 215 CALL dbcsr_allocate_matrix_set(lri_env%ri_sinv, 1) 216 ALLOCATE (lri_env%ri_sinv(1)%matrix) 217 SELECT CASE (lri_env%ri_sinv_app) 218 CASE ("NONE") 219 ! 220 CASE ("INVS") 221 CALL dbcsr_create(lri_env%ri_sinv(1)%matrix, template=lri_env%ri_smat(1)%matrix) 222 CALL invert_Hotelling(lri_env%ri_sinv(1)%matrix, lri_env%ri_smat(1)%matrix, & 223 threshold=1.e-10_dp, use_inv_as_guess=.FALSE., & 224 norm_convergence=1.e-10_dp, filter_eps=1.e-12_dp, silent=.FALSE.) 225 CASE ("INVF") 226 CALL dbcsr_create(emat, matrix_type=dbcsr_type_no_symmetry, template=lri_env%ri_smat(1)%matrix) 227 CALL get_qs_env(qs_env=qs_env, para_env=para_env, blacs_env=blacs_env) 228 CALL dbcsr_get_info(lri_env%ri_smat(1)%matrix, nfullrows_total=nbas) 229 ALLOCATE (eval(nbas)) 230 CALL cp_dbcsr_syevd(lri_env%ri_smat(1)%matrix, emat, eval, para_env, blacs_env) 231 izero = 0 232 DO i = 1, nbas 233 IF (eval(i) < 1.0e-10_dp) THEN 234 eval(i) = 0.0_dp 235 izero = izero + 1 236 ELSE 237 eval(i) = SQRT(1.0_dp/eval(i)) 238 ENDIF 239 END DO 240 CALL dbcsr_scale_by_vector(emat, eval, side='right') 241 CALL dbcsr_create(lri_env%ri_sinv(1)%matrix, template=lri_env%ri_smat(1)%matrix) 242 CALL dbcsr_multiply("N", "T", 1.0_dp, emat, emat, 0.0_dp, lri_env%ri_sinv(1)%matrix) 243 DEALLOCATE (eval) 244 CALL dbcsr_release(emat) 245 CASE ("AINV") 246 CALL dbcsr_create(lri_env%ri_sinv(1)%matrix, template=lri_env%ri_smat(1)%matrix) 247 CALL invert_Hotelling(lri_env%ri_sinv(1)%matrix, lri_env%ri_smat(1)%matrix, & 248 threshold=1.e-5_dp, use_inv_as_guess=.FALSE., & 249 norm_convergence=1.e-4_dp, filter_eps=1.e-4_dp, silent=.FALSE.) 250 CASE DEFAULT 251 CPABORT("Unknown RI_SINV type") 252 END SELECT 253 254 ! solve Rx=n 255 CALL ri_metric_solver(mat=lri_env%ri_smat(1)%matrix, & 256 vecr=lri_env%ri_fit%nvec, & 257 vecx=lri_env%ri_fit%rm1n, & 258 matp=lri_env%ri_sinv(1)%matrix, & 259 solver=lri_env%ri_sinv_app, & 260 ptr=lri_env%ri_fit%bas_ptr) 261 262 ! calculate n(t)x 263 lri_env%ri_fit%ntrm1n = SUM(lri_env%ri_fit%nvec(:)*lri_env%ri_fit%rm1n(:)) 264 265 ! calculate 3c-overlap integrals 266 IF (ASSOCIATED(lri_env%o3c)) THEN 267 CALL release_o3c_container(lri_env%o3c) 268 ELSE 269 ALLOCATE (lri_env%o3c) 270 END IF 271 CALL init_o3c_container(lri_env%o3c, dft_control%nspins, & 272 lri_env%orb_basis, lri_env%orb_basis, lri_env%ri_basis, & 273 lri_env%soo_list, lri_env%soa_list) 274 275 CALL calculate_o3c_integrals(lri_env%o3c, calculate_forces, matrix_p) 276 277 CALL timestop(handle) 278 279 END SUBROUTINE calculate_ri_integrals 280! ************************************************************************************************** 281!> \brief solver for RI systems (R*x=n) 282!> \param mat ... 283!> \param vecr ... 284!> \param vecx ... 285!> \param matp ... 286!> \param solver ... 287!> \param ptr ... 288! ************************************************************************************************** 289 SUBROUTINE ri_metric_solver(mat, vecr, vecx, matp, solver, ptr) 290 291 TYPE(dbcsr_type) :: mat 292 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: vecr 293 REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: vecx 294 TYPE(dbcsr_type) :: matp 295 CHARACTER(LEN=*), INTENT(IN) :: solver 296 INTEGER, DIMENSION(:, :), INTENT(IN) :: ptr 297 298 CHARACTER(LEN=*), PARAMETER :: routineN = 'ri_metric_solver', & 299 routineP = moduleN//':'//routineN 300 301 INTEGER :: handle, max_iter, n 302 LOGICAL :: converged 303 REAL(KIND=dp) :: rerror, threshold 304 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: vect 305 306 CALL timeset(routineN, handle) 307 308 threshold = 1.e-10_dp 309 max_iter = 100 310 311 vecx(:) = vecr(:) 312 313 SELECT CASE (solver) 314 CASE ("NONE") 315 CALL arnoldi_conjugate_gradient(mat, vecx, & 316 converged=converged, threshold=threshold, max_iter=max_iter) 317 CASE ("INVS", "INVF") 318 converged = .FALSE. 319 CALL ri_matvec(matp, vecr, vecx, ptr) 320 CASE ("AINV") 321 CALL arnoldi_conjugate_gradient(mat, vecx, matrix_p=matp, & 322 converged=converged, threshold=threshold, max_iter=max_iter) 323 CASE DEFAULT 324 CPABORT("Unknown RI solver") 325 END SELECT 326 327 IF (.NOT. converged) THEN 328 ! get error 329 rerror = 0.0_dp 330 n = SIZE(vecr) 331 ALLOCATE (vect(n)) 332 CALL ri_matvec(mat, vecx, vect, ptr) 333 vect(:) = vect(:) - vecr(:) 334 rerror = MAXVAL(ABS(vect(:))) 335 DEALLOCATE (vect) 336 IF (rerror > threshold) THEN 337 CPWARN("RI solver: CG did not converge properly") 338 END IF 339 END IF 340 341 CALL timestop(handle) 342 343 END SUBROUTINE ri_metric_solver 344 345! ************************************************************************************************** 346!> \brief performs the fitting of the density and distributes the fitted 347!> density on the grid 348!> \param lri_env the lri environment 349!> lri_density the environment for the fitting 350!> pmatrix density matrix 351!> lri_rho_struct where the fitted density is stored 352!> \param qs_env ... 353!> \param pmatrix ... 354!> \param lri_rho_struct ... 355!> \param atomic_kind_set ... 356!> \param para_env ... 357! ************************************************************************************************** 358 SUBROUTINE calculate_ri_densities(lri_env, qs_env, pmatrix, & 359 lri_rho_struct, atomic_kind_set, para_env) 360 361 TYPE(lri_environment_type), POINTER :: lri_env 362 TYPE(qs_environment_type), POINTER :: qs_env 363 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: pmatrix 364 TYPE(qs_rho_type), POINTER :: lri_rho_struct 365 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 366 TYPE(cp_para_env_type), POINTER :: para_env 367 368 CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_ri_densities', & 369 routineP = moduleN//':'//routineN 370 371 INTEGER :: atom_a, handle, i1, i2, iatom, ikind, & 372 ispin, n, natom, nspin 373 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of 374 INTEGER, DIMENSION(:, :), POINTER :: bas_ptr 375 REAL(KIND=dp), DIMENSION(:), POINTER :: tot_rho_r 376 REAL(KIND=dp), DIMENSION(:, :), POINTER :: avec 377 TYPE(lri_density_type), POINTER :: lri_density 378 TYPE(lri_kind_type), DIMENSION(:), POINTER :: lri_coef 379 TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_g, rho_r 380 381 CALL timeset(routineN, handle) 382 383 nspin = SIZE(pmatrix, 1) 384 CALL contract12_o3c(lri_env%o3c, pmatrix) 385 CALL calculate_tvec_ri(lri_env, lri_env%o3c, para_env) 386 CALL calculate_avec_ri(lri_env, pmatrix) 387 ! 388 CALL get_qs_env(qs_env, lri_density=lri_density, natom=natom) 389 CALL lri_density_release(lri_density) 390 CALL lri_density_create(lri_density) 391 lri_density%nspin = nspin 392 ! allocate the arrays to hold RI expansion coefficients lri_coefs 393 CALL allocate_lri_coefs(lri_env, lri_density, atomic_kind_set) 394 ! reassign avec 395 avec => lri_env%ri_fit%avec 396 bas_ptr => lri_env%ri_fit%bas_ptr 397 ALLOCATE (atom_of_kind(natom), kind_of(natom)) 398 CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of) 399 DO ispin = 1, nspin 400 lri_coef => lri_density%lri_coefs(ispin)%lri_kinds 401 DO iatom = 1, natom 402 ikind = kind_of(iatom) 403 atom_a = atom_of_kind(iatom) 404 i1 = bas_ptr(1, iatom) 405 i2 = bas_ptr(2, iatom) 406 n = i2 - i1 + 1 407 lri_coef(ikind)%acoef(atom_a, 1:n) = avec(i1:i2, ispin) 408 END DO 409 END DO 410 CALL set_qs_env(qs_env, lri_density=lri_density) 411 DEALLOCATE (atom_of_kind, kind_of) 412 ! 413 CALL qs_rho_get(lri_rho_struct, rho_r=rho_r, rho_g=rho_g, tot_rho_r=tot_rho_r) 414 DO ispin = 1, nspin 415 lri_coef => lri_density%lri_coefs(ispin)%lri_kinds 416 CALL calculate_lri_rho_elec(rho_g(ispin), rho_r(ispin), qs_env, & 417 lri_coef, tot_rho_r(ispin), "RI_HXC", .FALSE.) 418 ENDDO 419 420 CALL timestop(handle) 421 422 END SUBROUTINE calculate_ri_densities 423 424! ************************************************************************************************** 425!> \brief assembles the global vector t=<P.T> 426!> \param lri_env the lri environment 427!> \param o3c ... 428!> \param para_env ... 429! ************************************************************************************************** 430 SUBROUTINE calculate_tvec_ri(lri_env, o3c, para_env) 431 432 TYPE(lri_environment_type), POINTER :: lri_env 433 TYPE(o3c_container_type), POINTER :: o3c 434 TYPE(cp_para_env_type), POINTER :: para_env 435 436 CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_tvec_ri', & 437 routineP = moduleN//':'//routineN 438 INTEGER, PARAMETER :: msweep = 32 439 440 INTEGER :: handle, i1, i2, ibl, ibu, il, ispin, & 441 isweep, it, iu, katom, m, ma, mba, & 442 mepos, natom, nspin, nsweep, nthread 443 INTEGER, DIMENSION(2, msweep) :: nlimit 444 INTEGER, DIMENSION(:, :), POINTER :: bas_ptr 445 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: ta 446 REAL(KIND=dp), DIMENSION(:, :), POINTER :: rm1t, tvec, tvl 447 TYPE(o3c_iterator_type) :: o3c_iterator 448 449 CALL timeset(routineN, handle) 450 451 nspin = SIZE(lri_env%ri_fit%tvec, 2) 452 bas_ptr => lri_env%ri_fit%bas_ptr 453 tvec => lri_env%ri_fit%tvec 454 tvec = 0.0_dp 455 natom = SIZE(bas_ptr, 2) 456 nthread = 1 457!$ nthread = omp_get_max_threads() 458 459 IF (natom < 1000) THEN 460 nsweep = 1 461 ELSE 462 nsweep = MIN(nthread, msweep) 463 END IF 464 465 nlimit = 0 466 DO isweep = 1, nsweep 467 nlimit(1:2, isweep) = get_limit(natom, nsweep, isweep - 1) 468 END DO 469 470 DO ispin = 1, nspin 471 DO isweep = 1, nsweep 472 il = nlimit(1, isweep) 473 iu = nlimit(2, isweep) 474 ma = iu - il + 1 475 IF (ma < 1) CYCLE 476 ibl = bas_ptr(1, il) 477 ibu = bas_ptr(2, iu) 478 mba = ibu - ibl + 1 479 ALLOCATE (ta(mba, nthread)) 480 ta = 0.0_dp 481 482 CALL o3c_iterator_create(o3c, o3c_iterator, nthread=nthread) 483 484!$OMP PARALLEL DEFAULT(NONE)& 485!$OMP SHARED (nthread,o3c_iterator,ispin,il,iu,ibl,ta,bas_ptr)& 486!$OMP PRIVATE (mepos,katom,tvl,i1,i2,m) 487 488 mepos = 0 489!$ mepos = omp_get_thread_num() 490 491 DO WHILE (o3c_iterate(o3c_iterator, mepos=mepos) == 0) 492 CALL get_o3c_iterator_info(o3c_iterator, mepos=mepos, katom=katom, tvec=tvl) 493 IF (katom < il .OR. katom > iu) CYCLE 494 i1 = bas_ptr(1, katom) - ibl + 1 495 i2 = bas_ptr(2, katom) - ibl + 1 496 m = i2 - i1 + 1 497 ta(i1:i2, mepos + 1) = ta(i1:i2, mepos + 1) + tvl(1:m, ispin) 498 END DO 499!$OMP END PARALLEL 500 CALL o3c_iterator_release(o3c_iterator) 501 502 ! sum over threads 503 DO it = 1, nthread 504 tvec(ibl:ibu, ispin) = tvec(ibl:ibu, ispin) + ta(1:mba, it) 505 END DO 506 DEALLOCATE (ta) 507 END DO 508 END DO 509 510 ! global summation 511 CALL mp_sum(tvec, para_env%group) 512 513 rm1t => lri_env%ri_fit%rm1t 514 DO ispin = 1, nspin 515 ! solve Rx=t 516 CALL ri_metric_solver(mat=lri_env%ri_smat(1)%matrix, & 517 vecr=tvec(:, ispin), & 518 vecx=rm1t(:, ispin), & 519 matp=lri_env%ri_sinv(1)%matrix, & 520 solver=lri_env%ri_sinv_app, & 521 ptr=bas_ptr) 522 END DO 523 524 CALL timestop(handle) 525 526 END SUBROUTINE calculate_tvec_ri 527! ************************************************************************************************** 528!> \brief performs the fitting of the density in the RI method 529!> \param lri_env the lri environment 530!> lri_density the environment for the fitting 531!> pmatrix density matrix 532!> \param pmatrix ... 533! ************************************************************************************************** 534 SUBROUTINE calculate_avec_ri(lri_env, pmatrix) 535 536 TYPE(lri_environment_type), POINTER :: lri_env 537 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: pmatrix 538 539 CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_avec_ri', & 540 routineP = moduleN//':'//routineN 541 542 INTEGER :: handle, ispin, nspin 543 REAL(KIND=dp) :: etr, nelec, nrm1t 544 REAL(KIND=dp), DIMENSION(:), POINTER :: nvec, rm1n 545 REAL(KIND=dp), DIMENSION(:, :), POINTER :: avec, rm1t, tvec 546 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: smatrix 547 548 CALL timeset(routineN, handle) 549 550 nspin = SIZE(pmatrix) 551 ! number of electrons 552 smatrix => lri_env%ob_smat 553 DO ispin = 1, nspin 554 etr = 0.0_dp 555 CALL dbcsr_dot(smatrix(1)%matrix, pmatrix(ispin)%matrix, etr) 556 lri_env%ri_fit%echarge(ispin) = etr 557 END DO 558 559 tvec => lri_env%ri_fit%tvec 560 rm1t => lri_env%ri_fit%rm1t 561 nvec => lri_env%ri_fit%nvec 562 rm1n => lri_env%ri_fit%rm1n 563 564 ! calculate lambda 565 DO ispin = 1, nspin 566 nelec = lri_env%ri_fit%echarge(ispin) 567 nrm1t = SUM(nvec(:)*rm1t(:, ispin)) 568 lri_env%ri_fit%lambda(ispin) = 2.0_dp*(nrm1t - nelec)/lri_env%ri_fit%ntrm1n 569 END DO 570 571 ! calculate avec = rm1t - lambda/2 * rm1n 572 avec => lri_env%ri_fit%avec 573 DO ispin = 1, nspin 574 avec(:, ispin) = rm1t(:, ispin) - 0.5_dp*lri_env%ri_fit%lambda(ispin)*rm1n(:) 575 END DO 576 577 CALL timestop(handle) 578 579 END SUBROUTINE calculate_avec_ri 580 581! ************************************************************************************************** 582!> \brief Mutiplies a replicated vector with a DBCSR matrix: vo = mat*vi 583!> \param mat ... 584!> \param vi ... 585!> \param vo ... 586!> \param ptr ... 587! ************************************************************************************************** 588 SUBROUTINE ri_matvec(mat, vi, vo, ptr) 589 590 TYPE(dbcsr_type) :: mat 591 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: vi 592 REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: vo 593 INTEGER, DIMENSION(:, :), INTENT(IN) :: ptr 594 595 CHARACTER(LEN=*), PARAMETER :: routineN = 'ri_matvec', routineP = moduleN//':'//routineN 596 597 CHARACTER :: matrix_type 598 INTEGER :: group, handle, iatom, jatom, m1, m2, mb, & 599 n1, n2, nb 600 LOGICAL :: symm 601 REAL(KIND=dp), DIMENSION(:, :), POINTER :: block 602 TYPE(dbcsr_iterator_type) :: iter 603 604 CALL timeset(routineN, handle) 605 606 CALL dbcsr_get_info(mat, matrix_type=matrix_type, group=group) 607 608 SELECT CASE (matrix_type) 609 CASE (dbcsr_type_no_symmetry) 610 symm = .FALSE. 611 CASE (dbcsr_type_symmetric) 612 symm = .TRUE. 613 CASE (dbcsr_type_antisymmetric) 614 CPABORT("NYI, antisymmetric matrix not permitted") 615 CASE DEFAULT 616 CPABORT("Unknown matrix type, ...") 617 END SELECT 618 619 vo(:) = 0.0_dp 620 CALL dbcsr_iterator_start(iter, mat) 621 DO WHILE (dbcsr_iterator_blocks_left(iter)) 622 CALL dbcsr_iterator_next_block(iter, iatom, jatom, block) 623 n1 = ptr(1, iatom) 624 n2 = ptr(2, iatom) 625 nb = n2 - n1 + 1 626 m1 = ptr(1, jatom) 627 m2 = ptr(2, jatom) 628 mb = m2 - m1 + 1 629 CPASSERT(nb == SIZE(block, 1)) 630 CPASSERT(mb == SIZE(block, 2)) 631 vo(n1:n2) = vo(n1:n2) + MATMUL(block, vi(m1:m2)) 632 IF (symm .AND. (iatom /= jatom)) THEN 633 vo(m1:m2) = vo(m1:m2) + MATMUL(TRANSPOSE(block), vi(n1:n2)) 634 END IF 635 END DO 636 CALL dbcsr_iterator_stop(iter) 637 638 CALL mp_sum(vo, group) 639 640 CALL timestop(handle) 641 642 END SUBROUTINE ri_matvec 643 644END MODULE ri_environment_methods 645