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 Nabla V_KS (local part if PSP) on the different grids 8!> \par History 9!> created 06-2007 [RD] 10!> \author RD 11! ************************************************************************************************** 12MODULE qs_linres_epr_nablavks 13 USE atomic_kind_types, ONLY: atomic_kind_type,& 14 get_atomic_kind 15 USE cell_types, ONLY: cell_type 16 USE cp_control_types, ONLY: dft_control_type 17 USE cp_log_handling, ONLY: cp_get_default_logger,& 18 cp_logger_type 19 USE cp_output_handling, ONLY: cp_p_file,& 20 cp_print_key_finished_output,& 21 cp_print_key_should_output,& 22 cp_print_key_unit_nr 23 USE cp_para_types, ONLY: cp_para_env_type 24 USE cp_realspace_grid_cube, ONLY: cp_pw_to_cube 25 USE dbcsr_api, ONLY: dbcsr_p_type 26 USE external_potential_types, ONLY: all_potential_type,& 27 get_potential,& 28 gth_potential_type,& 29 sgp_potential_type 30 USE hartree_local_methods, ONLY: calculate_Vh_1center 31 USE input_section_types, ONLY: section_get_ivals,& 32 section_vals_get_subs_vals,& 33 section_vals_type,& 34 section_vals_val_get 35 USE kinds, ONLY: default_string_length,& 36 dp 37 USE mathconstants, ONLY: rootpi,& 38 twopi 39 USE particle_list_types, ONLY: particle_list_type 40 USE particle_types, ONLY: particle_type 41 USE pw_env_types, ONLY: pw_env_get,& 42 pw_env_type 43 USE pw_methods, ONLY: pw_axpy,& 44 pw_copy,& 45 pw_derive,& 46 pw_transfer,& 47 pw_zero 48 USE pw_poisson_methods, ONLY: pw_poisson_solve 49 USE pw_poisson_types, ONLY: pw_poisson_type 50 USE pw_pool_types, ONLY: pw_pool_create_pw,& 51 pw_pool_give_back_pw,& 52 pw_pool_type 53 USE pw_types, ONLY: COMPLEXDATA1D,& 54 REALDATA3D,& 55 REALSPACE,& 56 RECIPROCALSPACE,& 57 pw_p_type,& 58 pw_type 59 USE qs_environment_types, ONLY: get_qs_env,& 60 qs_environment_type 61 USE qs_gapw_densities, ONLY: prepare_gapw_den 62 USE qs_grid_atom, ONLY: grid_atom_type 63 USE qs_harmonics_atom, ONLY: harmonics_atom_type 64 USE qs_kind_types, ONLY: get_qs_kind,& 65 qs_kind_type 66 USE qs_ks_methods, ONLY: calc_rho_tot_gspace 67 USE qs_ks_types, ONLY: qs_ks_env_type 68 USE qs_linres_types, ONLY: epr_env_type,& 69 get_epr_env,& 70 nablavks_atom_type 71! R0 72 USE qs_local_rho_types, ONLY: rhoz_type 73 USE qs_rho0_types, ONLY: rho0_atom_type 74 USE qs_rho_atom_methods, ONLY: calculate_rho_atom_coeff 75 USE qs_rho_atom_types, ONLY: get_rho_atom,& 76 rho_atom_coeff,& 77 rho_atom_type 78! R0 79 USE qs_rho_types, ONLY: qs_rho_get,& 80 qs_rho_p_type,& 81 qs_rho_type 82 USE qs_subsys_types, ONLY: qs_subsys_get,& 83 qs_subsys_type 84 USE qs_vxc, ONLY: qs_vxc_create 85 USE qs_vxc_atom, ONLY: calculate_vxc_atom 86 USE util, ONLY: get_limit 87#include "./base/base_uses.f90" 88 89 IMPLICIT NONE 90 91 PRIVATE 92 PUBLIC :: epr_nablavks 93 94 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_linres_epr_nablavks' 95 96CONTAINS 97 98! ************************************************************************************************** 99!> \brief Evaluates Nabla V_KS on the grids 100!> \param epr_env ... 101!> \param qs_env ... 102!> \par History 103!> 06.2006 created [RD] 104!> \author RD 105! ************************************************************************************************** 106 SUBROUTINE epr_nablavks(epr_env, qs_env) 107 108 TYPE(epr_env_type) :: epr_env 109 TYPE(qs_environment_type), POINTER :: qs_env 110 111 CHARACTER(LEN=*), PARAMETER :: routineN = 'epr_nablavks', routineP = moduleN//':'//routineN 112 113 CHARACTER(LEN=default_string_length) :: ext, filename 114 COMPLEX(KIND=dp) :: gtemp 115 INTEGER :: bo_atom(2), ia, iat, iatom, idir, iexp, & 116 ig, ikind, ir, iso, ispin, ix, iy, iz, & 117 natom, nexp_ppl, nkind, nspins, & 118 output_unit, unit_nr 119 INTEGER, DIMENSION(2, 3) :: bo 120 INTEGER, DIMENSION(:), POINTER :: atom_list 121 LOGICAL :: gapw, gapw_xc, gth_gspace, ionode, & 122 make_soft, mpi_io, paw_atom 123 REAL(KIND=dp) :: alpha, alpha_core, arg, charge, ehartree, exc, exp_rap, gapw_max_alpha, & 124 hard_radius, hard_value, soft_value, sqrt_alpha, sqrt_rap 125 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: vh1_rad_h, vh1_rad_s 126 REAL(KIND=dp), DIMENSION(3) :: rap, ratom, roffset, rpoint 127 REAL(KIND=dp), DIMENSION(:), POINTER :: cexp_ppl, rho_rad_z 128 REAL(KIND=dp), DIMENSION(:, :), POINTER :: rho_rad_0 129 TYPE(all_potential_type), POINTER :: all_potential 130 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 131 TYPE(cell_type), POINTER :: cell 132 TYPE(cp_logger_type), POINTER :: logger 133 TYPE(cp_para_env_type), POINTER :: para_env 134 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao 135 TYPE(dft_control_type), POINTER :: dft_control 136 TYPE(grid_atom_type), POINTER :: grid_atom 137 TYPE(gth_potential_type), POINTER :: gth_potential 138 TYPE(harmonics_atom_type), POINTER :: harmonics 139 TYPE(nablavks_atom_type), DIMENSION(:), POINTER :: nablavks_atom_set 140 TYPE(particle_list_type), POINTER :: particles 141 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 142 TYPE(pw_env_type), POINTER :: pw_env 143 TYPE(pw_p_type), DIMENSION(:), POINTER :: rho1_r, rho2_r, rho_r, v_rspace_new, & 144 v_tau_rspace 145 TYPE(pw_p_type), POINTER :: rho_tot_gspace, v_coulomb_gspace, v_coulomb_gtemp, & 146 v_coulomb_rtemp, v_hartree_gspace, v_hartree_gtemp, v_hartree_rtemp, v_xc_gtemp, & 147 v_xc_rtemp, wf_r 148 TYPE(pw_poisson_type), POINTER :: poisson_env 149 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 150 TYPE(pw_type), POINTER :: pwx, pwy, pwz 151 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 152 TYPE(qs_ks_env_type), POINTER :: ks_env 153 TYPE(qs_rho_p_type), DIMENSION(:, :), POINTER :: nablavks_set 154 TYPE(qs_rho_type), POINTER :: rho, rho_xc 155 TYPE(qs_subsys_type), POINTER :: subsys 156 TYPE(rho0_atom_type), DIMENSION(:), POINTER :: rho0_atom_set 157 TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: rho_rad_h, rho_rad_s 158 TYPE(rho_atom_coeff), DIMENSION(:, :), POINTER :: nablavks_vec_rad_h, nablavks_vec_rad_s 159 TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set 160 TYPE(rho_atom_type), POINTER :: rho_atom 161 TYPE(rhoz_type), DIMENSION(:), POINTER :: rhoz_set 162 TYPE(section_vals_type), POINTER :: g_section, input, lr_section, xc_section 163 TYPE(sgp_potential_type), POINTER :: sgp_potential 164 165! R0 166! R0 167! R0 168! R0 169! R0 170! R0 171 172 NULLIFY (auxbas_pw_pool) 173 NULLIFY (cell) 174 NULLIFY (dft_control) 175 NULLIFY (g_section) 176 NULLIFY (logger) 177 NULLIFY (lr_section) 178 NULLIFY (nablavks_set) 179 NULLIFY (nablavks_atom_set) ! R0 180 NULLIFY (nablavks_vec_rad_h) ! R0 181 NULLIFY (nablavks_vec_rad_s) ! R0 182 NULLIFY (para_env) 183 NULLIFY (particle_set) 184 NULLIFY (particles) 185 NULLIFY (poisson_env) 186 NULLIFY (pw_env) 187 NULLIFY (pwx) ! R0 188 NULLIFY (pwy) ! R0 189 NULLIFY (pwz) ! R0 190 NULLIFY (rho) 191 NULLIFY (rho_xc) 192 NULLIFY (rho0_atom_set) 193 NULLIFY (rho_atom_set) 194 NULLIFY (rhoz_set) 195 NULLIFY (subsys) 196 NULLIFY (v_rspace_new) 197 NULLIFY (v_tau_rspace) 198 NULLIFY (xc_section) 199 NULLIFY (input) 200 NULLIFY (ks_env) 201 NULLIFY (rho_r, rho_ao, rho1_r, rho2_r) 202 203 logger => cp_get_default_logger() 204 lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES") 205 ionode = logger%para_env%ionode 206 207 output_unit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", & 208 extension=".linresLog") 209 210! ------------------------------------- 211! Read settings 212! ------------------------------------- 213 214 g_section => section_vals_get_subs_vals(lr_section, & 215 "EPR%PRINT%G_TENSOR") 216 217 CALL section_vals_val_get(g_section, "gapw_max_alpha", r_val=gapw_max_alpha) 218 219 gth_gspace = .TRUE. 220 221! ------------------------------------- 222! Get nablavks arrays 223! ------------------------------------- 224 225 CALL get_epr_env(epr_env, nablavks_set=nablavks_set, & ! R0 226 nablavks_atom_set=nablavks_atom_set) ! R0 227 ! R0 228 229 DO ispin = 1, SIZE(nablavks_set, 2) 230 DO idir = 1, SIZE(nablavks_set, 1) 231 CALL qs_rho_get(nablavks_set(idir, ispin)%rho, rho_r=rho_r) 232 CALL pw_zero(rho_r(1)%pw) 233 END DO 234 END DO 235 236 CALL qs_rho_get(nablavks_set(1, 1)%rho, rho_r=rho_r) 237 pwx => rho_r(1)%pw 238 CALL qs_rho_get(nablavks_set(2, 1)%rho, rho_r=rho_r) 239 pwy => rho_r(1)%pw 240 CALL qs_rho_get(nablavks_set(3, 1)%rho, rho_r=rho_r) 241 pwz => rho_r(1)%pw 242 243 roffset = -REAL(MODULO(pwx%pw_grid%npts, 2), dp)*pwx%pw_grid%dr/2.0_dp 244 245! ------------------------------------- 246! Get grids / atom info 247! ------------------------------------- 248 249 CALL get_qs_env(qs_env=qs_env, & 250 atomic_kind_set=atomic_kind_set, & 251 qs_kind_set=qs_kind_set, & 252 input=input, & 253 cell=cell, & 254 dft_control=dft_control, & 255 para_env=para_env, & 256 particle_set=particle_set, & 257 pw_env=pw_env, & 258 rho=rho, & 259 rho_xc=rho_xc, & 260 rho_atom_set=rho_atom_set, & 261 rho0_atom_set=rho0_atom_set, & 262 rhoz_set=rhoz_set, & 263 subsys=subsys, & 264 ks_env=ks_env) 265 266 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, & 267 poisson_env=poisson_env) 268 269 CALL qs_subsys_get(subsys, particles=particles) 270 271 gapw = dft_control%qs_control%gapw 272 gapw_xc = dft_control%qs_control%gapw_xc 273 nkind = SIZE(atomic_kind_set) 274 nspins = dft_control%nspins 275 276! ------------------------------------- 277! Add Hartree potential 278! ------------------------------------- 279 280 ALLOCATE (v_hartree_gspace) 281 ALLOCATE (v_hartree_gtemp) 282 ALLOCATE (v_hartree_rtemp) 283 ALLOCATE (rho_tot_gspace) 284 285 CALL pw_pool_create_pw(auxbas_pw_pool, v_hartree_gspace%pw, & 286 use_data=COMPLEXDATA1D, in_space=RECIPROCALSPACE) 287 CALL pw_pool_create_pw(auxbas_pw_pool, v_hartree_gtemp%pw, & 288 use_data=COMPLEXDATA1D, in_space=RECIPROCALSPACE) 289 CALL pw_pool_create_pw(auxbas_pw_pool, v_hartree_rtemp%pw, & 290 use_data=REALDATA3D, in_space=REALSPACE) 291 CALL pw_pool_create_pw(auxbas_pw_pool, rho_tot_gspace%pw, & 292 use_data=COMPLEXDATA1D, in_space=RECIPROCALSPACE) 293 294 IF (gapw) THEN 295 ! need to rebuild the coeff ! 296 CALL qs_rho_get(rho, rho_ao_kp=rho_ao) 297 CALL calculate_rho_atom_coeff(qs_env, rho_ao) 298 CALL prepare_gapw_den(qs_env, do_rho0=.TRUE.) 299 END IF 300 301 CALL pw_zero(rho_tot_gspace%pw) 302 303 CALL calc_rho_tot_gspace(rho_tot_gspace, qs_env, rho, & 304 skip_nuclear_density=.NOT. gapw) 305 306 CALL pw_poisson_solve(poisson_env, rho_tot_gspace%pw, ehartree, & 307 v_hartree_gspace%pw) 308 309 ! ------------------------------------- 310 ! Atomic grids part 311 ! ------------------------------------- 312 313 IF (gapw) THEN 314 315 DO ikind = 1, nkind ! loop over atom types 316 317 NULLIFY (atom_list) 318 NULLIFY (grid_atom) 319 NULLIFY (harmonics) 320 NULLIFY (rho_rad_z) 321 322 rho_rad_z => rhoz_set(ikind)%r_coef 323 324 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom) 325 CALL get_qs_kind(qs_kind_set(ikind), & 326 grid_atom=grid_atom, & 327 harmonics=harmonics, & 328 hard_radius=hard_radius, & 329 paw_atom=paw_atom, & 330 zeff=charge, & 331 alpha_core_charge=alpha_core) 332 333 IF (paw_atom) THEN 334 335 ALLOCATE (vh1_rad_h(grid_atom%nr, harmonics%max_iso_not0)) 336 ALLOCATE (vh1_rad_s(grid_atom%nr, harmonics%max_iso_not0)) 337 338 ! DO iat = 1, natom ! natom = # atoms for ikind 339 ! 340 ! iatom = atom_list(iat) 341 ! ratom = particle_set(iatom)%r 342 ! 343 ! DO ig = v_hartree_gspace%pw%pw_grid%first_gne0,v_hartree_gspace%pw%pw_grid%ngpts_cut_local 344 ! 345 ! gtemp = fourpi * charge / cell%deth * & 346 ! EXP ( - v_hartree_gspace%pw%pw_grid%gsq(ig) / (4.0_dp * alpha_core) ) & 347 ! / v_hartree_gspace%pw%pw_grid%gsq(ig) 348 ! 349 ! arg = DOT_PRODUCT(v_hartree_gspace%pw%pw_grid%g(:,ig),ratom) 350 ! 351 ! gtemp = gtemp * CMPLX(COS(arg),-SIN(arg),KIND=dp) 352 ! 353 ! v_hartree_gspace%pw%cc(ig) = v_hartree_gspace%pw%cc(ig) + gtemp 354 ! END DO 355 ! IF ( v_hartree_gspace%pw%pw_grid%have_g0 ) v_hartree_gspace%pw%cc(1) = 0.0_dp 356 ! 357 ! END DO 358 359 bo_atom = get_limit(natom, para_env%num_pe, para_env%mepos) 360 361 DO iat = bo_atom(1), bo_atom(2) ! natomkind = # atoms for ikind 362 363 iatom = atom_list(iat) 364 ratom = particle_set(iatom)%r 365 366 nablavks_vec_rad_h => nablavks_atom_set(iatom)%nablavks_vec_rad_h 367 nablavks_vec_rad_s => nablavks_atom_set(iatom)%nablavks_vec_rad_s 368 369 DO ispin = 1, nspins 370 DO idir = 1, 3 371 nablavks_vec_rad_h(idir, ispin)%r_coef(:, :) = 0.0_dp 372 nablavks_vec_rad_s(idir, ispin)%r_coef(:, :) = 0.0_dp 373 END DO ! idir 374 END DO ! ispin 375 376 rho_atom => rho_atom_set(iatom) 377 NULLIFY (rho_rad_h, rho_rad_s, rho_rad_0) 378 CALL get_rho_atom(rho_atom=rho_atom, rho_rad_h=rho_rad_h, & 379 rho_rad_s=rho_rad_s) 380 rho_rad_0 => rho0_atom_set(iatom)%rho0_rad_h%r_coef 381 vh1_rad_h = 0.0_dp 382 vh1_rad_s = 0.0_dp 383 384 CALL calculate_Vh_1center(vh1_rad_h, vh1_rad_s, rho_rad_h, rho_rad_s, rho_rad_0, rho_rad_z, grid_atom) 385 386 DO ir = 2, grid_atom%nr 387 388 IF (grid_atom%rad(ir) >= hard_radius) CYCLE 389 390 DO ia = 1, grid_atom%ng_sphere 391 392 ! hard part 393 394 DO idir = 1, 3 395 hard_value = 0.0_dp 396 DO iso = 1, harmonics%max_iso_not0 397 hard_value = hard_value + & 398 vh1_rad_h(ir, iso)*harmonics%dslm_dxyz(idir, ia, iso) + & 399 harmonics%slm(ia, iso)* & 400 (vh1_rad_h(ir - 1, iso) - vh1_rad_h(ir, iso))/ & 401 (grid_atom%rad(ir - 1) - grid_atom%rad(ir))* & 402 (harmonics%a(idir, ia)) 403 END DO 404 nablavks_vec_rad_h(idir, 1)%r_coef(ir, ia) = hard_value 405 END DO 406 407 ! soft part 408 409 DO idir = 1, 3 410 soft_value = 0.0_dp 411 DO iso = 1, harmonics%max_iso_not0 412 soft_value = soft_value + & 413 vh1_rad_s(ir, iso)*harmonics%dslm_dxyz(idir, ia, iso) + & 414 harmonics%slm(ia, iso)* & 415 (vh1_rad_s(ir - 1, iso) - vh1_rad_s(ir, iso))/ & 416 (grid_atom%rad(ir - 1) - grid_atom%rad(ir))* & 417 (harmonics%a(idir, ia)) 418 END DO 419 nablavks_vec_rad_s(idir, 1)%r_coef(ir, ia) = soft_value 420 END DO 421 422 END DO ! ia 423 424 END DO ! ir 425 426 DO idir = 1, 3 427 nablavks_vec_rad_h(idir, 2)%r_coef(:, :) = nablavks_vec_rad_h(idir, 1)%r_coef(:, :) 428 nablavks_vec_rad_s(idir, 2)%r_coef(:, :) = nablavks_vec_rad_s(idir, 1)%r_coef(:, :) 429 END DO 430 431 END DO ! iat 432 433 DEALLOCATE (vh1_rad_h) 434 DEALLOCATE (vh1_rad_s) 435 436 END IF ! paw_atom 437 438 END DO ! ikind 439 440 END IF ! gapw 441 442 CALL pw_copy(v_hartree_gspace%pw, v_hartree_gtemp%pw) 443 CALL pw_derive(v_hartree_gtemp%pw, (/1, 0, 0/)) 444 CALL pw_transfer(v_hartree_gtemp%pw, v_hartree_rtemp%pw) 445 CALL pw_copy(v_hartree_rtemp%pw, pwx) 446 447 CALL pw_copy(v_hartree_gspace%pw, v_hartree_gtemp%pw) 448 CALL pw_derive(v_hartree_gtemp%pw, (/0, 1, 0/)) 449 CALL pw_transfer(v_hartree_gtemp%pw, v_hartree_rtemp%pw) 450 CALL pw_copy(v_hartree_rtemp%pw, pwy) 451 452 CALL pw_copy(v_hartree_gspace%pw, v_hartree_gtemp%pw) 453 CALL pw_derive(v_hartree_gtemp%pw, (/0, 0, 1/)) 454 CALL pw_transfer(v_hartree_gtemp%pw, v_hartree_rtemp%pw) 455 CALL pw_copy(v_hartree_rtemp%pw, pwz) 456 457 CALL pw_pool_give_back_pw(auxbas_pw_pool, v_hartree_gspace%pw) 458 CALL pw_pool_give_back_pw(auxbas_pw_pool, v_hartree_gtemp%pw) 459 CALL pw_pool_give_back_pw(auxbas_pw_pool, v_hartree_rtemp%pw) 460 CALL pw_pool_give_back_pw(auxbas_pw_pool, rho_tot_gspace%pw) 461 462 DEALLOCATE (v_hartree_gspace) 463 DEALLOCATE (v_hartree_gtemp) 464 DEALLOCATE (v_hartree_rtemp) 465 DEALLOCATE (rho_tot_gspace) 466 467! ------------------------------------- 468! Add Coulomb potential 469! ------------------------------------- 470 471 DO ikind = 1, nkind ! loop over atom types 472 473 NULLIFY (atom_list) 474 NULLIFY (grid_atom) 475 NULLIFY (harmonics) 476 477 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom) 478 CALL get_qs_kind(qs_kind_set(ikind), & 479 grid_atom=grid_atom, & 480 harmonics=harmonics, & 481 hard_radius=hard_radius, & 482 gth_potential=gth_potential, & 483 sgp_potential=sgp_potential, & 484 all_potential=all_potential, & 485 paw_atom=paw_atom) 486 487 IF (ASSOCIATED(gth_potential)) THEN 488 489 NULLIFY (cexp_ppl) 490 491 CALL get_potential(potential=gth_potential, & 492 zeff=charge, & 493 alpha_ppl=alpha, & 494 nexp_ppl=nexp_ppl, & 495 cexp_ppl=cexp_ppl) 496 497 sqrt_alpha = SQRT(alpha) 498 499 IF (gapw .AND. paw_atom .AND. alpha > gapw_max_alpha) THEN 500 make_soft = .TRUE. 501 ELSE 502 make_soft = .FALSE. 503 END IF 504 505 ! ------------------------------------- 506 ! PW grid part 507 ! ------------------------------------- 508 509 IF (gth_gspace) THEN 510 511 ALLOCATE (v_coulomb_gspace) 512 ALLOCATE (v_coulomb_gtemp) 513 ALLOCATE (v_coulomb_rtemp) 514 515 CALL pw_pool_create_pw(auxbas_pw_pool, v_coulomb_gspace%pw, & 516 use_data=COMPLEXDATA1D, in_space=RECIPROCALSPACE) 517 CALL pw_pool_create_pw(auxbas_pw_pool, v_coulomb_gtemp%pw, & 518 use_data=COMPLEXDATA1D, in_space=RECIPROCALSPACE) 519 CALL pw_pool_create_pw(auxbas_pw_pool, v_coulomb_rtemp%pw, & 520 use_data=REALDATA3D, in_space=REALSPACE) 521 522 CALL pw_zero(v_coulomb_gspace%pw) 523 524 DO iat = 1, natom ! natom = # atoms for ikind 525 526 iatom = atom_list(iat) 527 ratom = particle_set(iatom)%r 528 529 DO ig = v_coulomb_gspace%pw%pw_grid%first_gne0, v_coulomb_gspace%pw%pw_grid%ngpts_cut_local 530 gtemp = 0.0_dp 531 ! gtemp = - fourpi * charge / cell%deth * & 532 ! EXP ( - v_coulomb_gspace%pw%pw_grid%gsq(ig) / (4.0_dp * alpha) ) & 533 ! / v_coulomb_gspace%pw%pw_grid%gsq(ig) 534 535 IF (.NOT. make_soft) THEN 536 537 SELECT CASE (nexp_ppl) 538 CASE (1) 539 gtemp = gtemp + & 540 (twopi)**(1.5_dp)/(cell%deth*(2.0_dp*alpha)**(1.5_dp))* & 541 EXP(-v_coulomb_gspace%pw%pw_grid%gsq(ig)/(4.0_dp*alpha))*( & 542 ! C1 543 +cexp_ppl(1) & 544 ) 545 CASE (2) 546 gtemp = gtemp + & 547 (twopi)**(1.5_dp)/(cell%deth*(2.0_dp*alpha)**(1.5_dp))* & 548 EXP(-v_coulomb_gspace%pw%pw_grid%gsq(ig)/(4.0_dp*alpha))*( & 549 ! C1 550 +cexp_ppl(1) & 551 ! C2 552 + cexp_ppl(2)/(2.0_dp*alpha)* & 553 (3.0_dp - v_coulomb_gspace%pw%pw_grid%gsq(ig)/(2.0_dp*alpha)) & 554 ) 555 CASE (3) 556 gtemp = gtemp + & 557 (twopi)**(1.5_dp)/(cell%deth*(2.0_dp*alpha)**(1.5_dp))* & 558 EXP(-v_coulomb_gspace%pw%pw_grid%gsq(ig)/(4.0_dp*alpha))*( & 559 ! C1 560 +cexp_ppl(1) & 561 ! C2 562 + cexp_ppl(2)/(2.0_dp*alpha)* & 563 (3.0_dp - v_coulomb_gspace%pw%pw_grid%gsq(ig)/(2.0_dp*alpha)) & 564 ! C3 565 + cexp_ppl(3)/(2.0_dp*alpha)**2* & 566 (15.0_dp - 10.0_dp*v_coulomb_gspace%pw%pw_grid%gsq(ig)/(2.0_dp*alpha) & 567 + (v_coulomb_gspace%pw%pw_grid%gsq(ig)/(2.0_dp*alpha))**2) & 568 ) 569 CASE (4) 570 gtemp = gtemp + & 571 (twopi)**(1.5_dp)/(cell%deth*(2.0_dp*alpha)**(1.5_dp))* & 572 EXP(-v_coulomb_gspace%pw%pw_grid%gsq(ig)/(4.0_dp*alpha))*( & 573 ! C1 574 +cexp_ppl(1) & 575 ! C2 576 + cexp_ppl(2)/(2.0_dp*alpha)* & 577 (3.0_dp - v_coulomb_gspace%pw%pw_grid%gsq(ig)/(2.0_dp*alpha)) & 578 ! C3 579 + cexp_ppl(3)/(2.0_dp*alpha)**2* & 580 (15.0_dp - 10.0_dp*v_coulomb_gspace%pw%pw_grid%gsq(ig)/(2.0_dp*alpha) & 581 + (v_coulomb_gspace%pw%pw_grid%gsq(ig)/(2.0_dp*alpha))**2) & 582 ! C4 583 + cexp_ppl(4)/(2.0_dp*alpha)**3* & 584 (105.0_dp - 105.0_dp*v_coulomb_gspace%pw%pw_grid%gsq(ig)/(2.0_dp*alpha) & 585 + 21.0_dp*(v_coulomb_gspace%pw%pw_grid%gsq(ig)/(2.0_dp*alpha))**2 & 586 - (v_coulomb_gspace%pw%pw_grid%gsq(ig)/(2.0_dp*alpha))**3) & 587 ) 588 END SELECT 589 590 END IF 591 592 arg = DOT_PRODUCT(v_coulomb_gspace%pw%pw_grid%g(:, ig), ratom) 593 594 gtemp = gtemp*CMPLX(COS(arg), -SIN(arg), KIND=dp) 595 v_coulomb_gspace%pw%cc(ig) = v_coulomb_gspace%pw%cc(ig) + gtemp 596 END DO 597 IF (v_coulomb_gspace%pw%pw_grid%have_g0) v_coulomb_gspace%pw%cc(1) = 0.0_dp 598 599 END DO 600 601 CALL pw_copy(v_coulomb_gspace%pw, v_coulomb_gtemp%pw) 602 CALL pw_derive(v_coulomb_gtemp%pw, (/1, 0, 0/)) 603 CALL pw_transfer(v_coulomb_gtemp%pw, v_coulomb_rtemp%pw) 604 CALL pw_axpy(v_coulomb_rtemp%pw, pwx) 605 606 CALL pw_copy(v_coulomb_gspace%pw, v_coulomb_gtemp%pw) 607 CALL pw_derive(v_coulomb_gtemp%pw, (/0, 1, 0/)) 608 CALL pw_transfer(v_coulomb_gtemp%pw, v_coulomb_rtemp%pw) 609 CALL pw_axpy(v_coulomb_rtemp%pw, pwy) 610 611 CALL pw_copy(v_coulomb_gspace%pw, v_coulomb_gtemp%pw) 612 CALL pw_derive(v_coulomb_gtemp%pw, (/0, 0, 1/)) 613 CALL pw_transfer(v_coulomb_gtemp%pw, v_coulomb_rtemp%pw) 614 CALL pw_axpy(v_coulomb_rtemp%pw, pwz) 615 616 CALL pw_pool_give_back_pw(auxbas_pw_pool, v_coulomb_gspace%pw) 617 CALL pw_pool_give_back_pw(auxbas_pw_pool, v_coulomb_gtemp%pw) 618 CALL pw_pool_give_back_pw(auxbas_pw_pool, v_coulomb_rtemp%pw) 619 620 DEALLOCATE (v_coulomb_gspace) 621 DEALLOCATE (v_coulomb_gtemp) 622 DEALLOCATE (v_coulomb_rtemp) 623 ELSE 624 625 ! Attic of the atomic parallellisation 626 ! 627 ! bo(2) 628 ! bo = get_limit(natom, para_env%num_pe, para_env%mepos) 629 ! DO iat = bo(1),bo(2) ! natom = # atoms for ikind 630 ! DO ix = lbound(pwx%cr3d,1), ubound(pwx%cr3d,1) 631 ! DO iy = lbound(pwx%cr3d,2), ubound(pwx%cr3d,2) 632 ! DO iz = lbound(pwx%cr3d,3), ubound(pwx%cr3d,3) 633 634 bo = pwx%pw_grid%bounds_local 635 636 DO iat = 1, natom ! natom = # atoms for ikind 637 638 iatom = atom_list(iat) 639 ratom = particle_set(iatom)%r 640 641 DO ix = bo(1, 1), bo(2, 1) 642 DO iy = bo(1, 2), bo(2, 2) 643 DO iz = bo(1, 3), bo(2, 3) 644 rpoint = (/REAL(ix, dp)*pwx%pw_grid%dr(1), & 645 REAL(iy, dp)*pwx%pw_grid%dr(2), & 646 REAL(iz, dp)*pwx%pw_grid%dr(3)/) 647 rpoint = rpoint + roffset 648 rap = rpoint - ratom 649 rap(1) = MODULO(rap(1), cell%hmat(1, 1)) - cell%hmat(1, 1)/2._dp 650 rap(2) = MODULO(rap(2), cell%hmat(2, 2)) - cell%hmat(2, 2)/2._dp 651 rap(3) = MODULO(rap(3), cell%hmat(3, 3)) - cell%hmat(3, 3)/2._dp 652 sqrt_rap = SQRT(DOT_PRODUCT(rap, rap)) 653 exp_rap = EXP(-alpha*sqrt_rap**2) 654 sqrt_rap = MAX(sqrt_rap, 1.e-10_dp) 655 ! d_x 656 657 pwx%cr3d(ix, iy, iz) = pwx%cr3d(ix, iy, iz) + charge*( & 658 -2.0_dp*sqrt_alpha*EXP(-sqrt_rap**2*sqrt_alpha**2)*rap(1) & 659 /(rootpi*sqrt_rap**2) & 660 + erf(sqrt_rap*sqrt_alpha)*rap(1) & 661 /sqrt_rap**3) 662 663 ! d_y 664 665 pwy%cr3d(ix, iy, iz) = pwy%cr3d(ix, iy, iz) + charge*( & 666 -2.0_dp*sqrt_alpha*EXP(-sqrt_rap**2*sqrt_alpha**2)*rap(2) & 667 /(rootpi*sqrt_rap**2) & 668 + erf(sqrt_rap*sqrt_alpha)*rap(2) & 669 /sqrt_rap**3) 670 671 ! d_z 672 673 pwz%cr3d(ix, iy, iz) = pwz%cr3d(ix, iy, iz) + charge*( & 674 -2.0_dp*sqrt_alpha*EXP(-sqrt_rap**2*sqrt_alpha**2)*rap(3) & 675 /(rootpi*sqrt_rap**2) & 676 + erf(sqrt_rap*sqrt_alpha)*rap(3) & 677 /sqrt_rap**3) 678 679 IF (make_soft) CYCLE 680 681 ! d_x 682 683 DO iexp = 1, nexp_ppl 684 pwx%cr3d(ix, iy, iz) = pwx%cr3d(ix, iy, iz) + ( & 685 -2.0_dp*alpha*rap(1)*exp_rap* & 686 cexp_ppl(iexp)*(sqrt_rap**2)**(iexp - 1)) 687 IF (iexp > 1) THEN 688 pwx%cr3d(ix, iy, iz) = pwx%cr3d(ix, iy, iz) + ( & 689 2.0_dp*exp_rap*cexp_ppl(iexp)* & 690 (sqrt_rap**2)**(iexp - 2)*REAL(iexp - 1, dp)*rap(1)) 691 END IF 692 END DO 693 694 ! d_y 695 696 DO iexp = 1, nexp_ppl 697 pwy%cr3d(ix, iy, iz) = pwy%cr3d(ix, iy, iz) + ( & 698 -2.0_dp*alpha*rap(2)*exp_rap* & 699 cexp_ppl(iexp)*(sqrt_rap**2)**(iexp - 1)) 700 IF (iexp > 1) THEN 701 pwy%cr3d(ix, iy, iz) = pwy%cr3d(ix, iy, iz) + ( & 702 2.0_dp*exp_rap*cexp_ppl(iexp)* & 703 (sqrt_rap**2)**(iexp - 2)*REAL(iexp - 1, dp)*rap(2)) 704 END IF 705 END DO 706 707 ! d_z 708 709 DO iexp = 1, nexp_ppl 710 pwz%cr3d(ix, iy, iz) = pwz%cr3d(ix, iy, iz) + ( & 711 -2.0_dp*alpha*rap(3)*exp_rap* & 712 cexp_ppl(iexp)*(sqrt_rap**2)**(iexp - 1)) 713 IF (iexp > 1) THEN 714 pwz%cr3d(ix, iy, iz) = pwz%cr3d(ix, iy, iz) + ( & 715 2.0_dp*exp_rap*cexp_ppl(iexp)* & 716 (sqrt_rap**2)**(iexp - 2)*REAL(iexp - 1, dp)*rap(3)) 717 END IF 718 END DO 719 720 END DO ! iz 721 END DO ! iy 722 END DO ! ix 723 END DO ! iat 724 END IF ! gth_gspace 725 726 ! ------------------------------------- 727 ! Atomic grids part 728 ! ------------------------------------- 729 730 IF (gapw .AND. paw_atom) THEN 731 732 bo_atom = get_limit(natom, para_env%num_pe, para_env%mepos) 733 734 DO iat = bo_atom(1), bo_atom(2) ! natom = # atoms for ikind 735 736 iatom = atom_list(iat) 737 738 nablavks_vec_rad_h => nablavks_atom_set(iatom)%nablavks_vec_rad_h 739 nablavks_vec_rad_s => nablavks_atom_set(iatom)%nablavks_vec_rad_s 740 741 DO ir = 1, grid_atom%nr 742 743 IF (grid_atom%rad(ir) >= hard_radius) CYCLE 744 745 exp_rap = EXP(-alpha*grid_atom%rad(ir)**2) 746 747 DO ia = 1, grid_atom%ng_sphere 748 749 DO idir = 1, 3 750 hard_value = 0.0_dp 751 hard_value = charge*( & 752 -2.0_dp*sqrt_alpha*EXP(-grid_atom%rad(ir)**2*sqrt_alpha**2) & 753 *grid_atom%rad(ir)*harmonics%a(idir, ia) & 754 /(rootpi*grid_atom%rad(ir)**2) & 755 + erf(grid_atom%rad(ir)*sqrt_alpha) & 756 *grid_atom%rad(ir)*harmonics%a(idir, ia) & 757 /grid_atom%rad(ir)**3) 758 soft_value = hard_value 759 DO iexp = 1, nexp_ppl 760 hard_value = hard_value + ( & 761 -2.0_dp*alpha*grid_atom%rad(ir)*harmonics%a(idir, ia) & 762 *exp_rap*cexp_ppl(iexp)*(grid_atom%rad(ir)**2)**(iexp - 1)) 763 IF (iexp > 1) THEN 764 hard_value = hard_value + ( & 765 2.0_dp*exp_rap*cexp_ppl(iexp) & 766 *(grid_atom%rad(ir)**2)**(iexp - 2)*REAL(iexp - 1, dp) & 767 *grid_atom%rad(ir)*harmonics%a(idir, ia)) 768 END IF 769 END DO 770 nablavks_vec_rad_h(idir, 1)%r_coef(ir, ia) = & 771 nablavks_vec_rad_h(idir, 1)%r_coef(ir, ia) + hard_value 772 IF (make_soft) THEN 773 nablavks_vec_rad_s(idir, 1)%r_coef(ir, ia) = & 774 nablavks_vec_rad_s(idir, 1)%r_coef(ir, ia) + soft_value 775 ELSE 776 nablavks_vec_rad_s(idir, 1)%r_coef(ir, ia) = & 777 nablavks_vec_rad_s(idir, 1)%r_coef(ir, ia) + hard_value 778 END IF 779 END DO 780 781 END DO ! ia 782 END DO ! ir 783 784 DO ispin = 2, nspins 785 DO idir = 1, 3 786 nablavks_vec_rad_h(idir, ispin)%r_coef(:, :) = nablavks_vec_rad_h(idir, 1)%r_coef(:, :) 787 nablavks_vec_rad_s(idir, ispin)%r_coef(:, :) = nablavks_vec_rad_s(idir, 1)%r_coef(:, :) 788 END DO 789 END DO 790 791 END DO 792 793 END IF 794 795 ELSE IF (ASSOCIATED(sgp_potential)) THEN 796 797 CPABORT("EPR with SGP potentials is not implemented") 798 799 ELSE IF (ASSOCIATED(all_potential)) THEN 800 801 CALL get_potential(potential=all_potential, & 802 alpha_core_charge=alpha, & 803 zeff=charge) 804 805 sqrt_alpha = SQRT(alpha) 806 807 ! ------------------------------------- 808 ! Atomic grids part 809 ! ------------------------------------- 810 811 bo_atom = get_limit(natom, para_env%num_pe, para_env%mepos) 812 813 DO iat = bo_atom(1), bo_atom(2) ! natom = # atoms for ikind 814 815 iatom = atom_list(iat) 816 817 nablavks_vec_rad_h => nablavks_atom_set(iatom)%nablavks_vec_rad_h 818 819 DO ir = 1, grid_atom%nr 820 821 IF (grid_atom%rad(ir) >= hard_radius) CYCLE 822 823 DO ia = 1, grid_atom%ng_sphere 824 825 DO idir = 1, 3 826 hard_value = 0.0_dp 827 hard_value = charge*( & 828 2.0_dp*sqrt_alpha*EXP(-grid_atom%rad(ir)**2*sqrt_alpha**2) & 829 *grid_atom%rad(ir)*harmonics%a(idir, ia) & 830 /(rootpi*grid_atom%rad(ir)**2) & 831 + erfc(grid_atom%rad(ir)*sqrt_alpha) & 832 *grid_atom%rad(ir)*harmonics%a(idir, ia) & 833 /grid_atom%rad(ir)**3) 834 nablavks_vec_rad_h(idir, 1)%r_coef(ir, ia) = & 835 nablavks_vec_rad_h(idir, 1)%r_coef(ir, ia) + hard_value 836 END DO 837 838 END DO ! ia 839 END DO ! ir 840 841 DO ispin = 2, nspins 842 DO idir = 1, 3 843 nablavks_vec_rad_h(idir, ispin)%r_coef(:, :) = nablavks_vec_rad_h(idir, 1)%r_coef(:, :) 844 END DO 845 END DO 846 847 END DO 848 849 ELSE 850 CYCLE 851 END IF 852 853 END DO 854 855 DO idir = 1, 3 856 CALL qs_rho_get(nablavks_set(idir, 1)%rho, rho_r=rho1_r) 857 CALL qs_rho_get(nablavks_set(idir, 2)%rho, rho_r=rho2_r) 858 CALL pw_copy(rho1_r(1)%pw, rho2_r(1)%pw) 859 END DO 860 861! ------------------------------------- 862! Add V_xc potential 863! ------------------------------------- 864 865 ALLOCATE (v_xc_gtemp) 866 ALLOCATE (v_xc_rtemp) 867 868 CALL pw_pool_create_pw(auxbas_pw_pool, v_xc_gtemp%pw, & 869 use_data=COMPLEXDATA1D, in_space=RECIPROCALSPACE) 870 CALL pw_pool_create_pw(auxbas_pw_pool, v_xc_rtemp%pw, & 871 use_data=REALDATA3D, in_space=REALSPACE) 872 873 xc_section => section_vals_get_subs_vals(input, "PROPERTIES%LINRES%EPR%PRINT%G_TENSOR%XC") 874 ! a possible vdW section in xc_section will be ignored 875 876 IF (gapw_xc) THEN 877 CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho_xc, xc_section=xc_section, & 878 vxc_rho=v_rspace_new, vxc_tau=v_tau_rspace, exc=exc, just_energy=.FALSE.) 879 ELSE 880 CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho, xc_section=xc_section, & 881 vxc_rho=v_rspace_new, vxc_tau=v_tau_rspace, exc=exc, just_energy=.FALSE.) 882 END IF 883 884 IF (ASSOCIATED(v_rspace_new)) THEN 885 886 DO ispin = 1, nspins 887 888 CALL pw_transfer(v_rspace_new(ispin)%pw, v_xc_gtemp%pw) 889 CALL pw_derive(v_xc_gtemp%pw, (/1, 0, 0/)) 890 CALL pw_transfer(v_xc_gtemp%pw, v_xc_rtemp%pw) 891 CALL qs_rho_get(nablavks_set(1, ispin)%rho, rho_r=rho_r) 892 CALL pw_axpy(v_xc_rtemp%pw, rho_r(1)%pw) 893 894 CALL pw_transfer(v_rspace_new(ispin)%pw, v_xc_gtemp%pw) 895 CALL pw_derive(v_xc_gtemp%pw, (/0, 1, 0/)) 896 CALL pw_transfer(v_xc_gtemp%pw, v_xc_rtemp%pw) 897 CALL qs_rho_get(nablavks_set(2, ispin)%rho, rho_r=rho_r) 898 CALL pw_axpy(v_xc_rtemp%pw, rho_r(1)%pw) 899 900 CALL pw_transfer(v_rspace_new(ispin)%pw, v_xc_gtemp%pw) 901 CALL pw_derive(v_xc_gtemp%pw, (/0, 0, 1/)) 902 CALL pw_transfer(v_xc_gtemp%pw, v_xc_rtemp%pw) 903 CALL qs_rho_get(nablavks_set(3, ispin)%rho, rho_r=rho_r) 904 CALL pw_axpy(v_xc_rtemp%pw, rho_r(1)%pw) 905 906 CALL pw_pool_give_back_pw(auxbas_pw_pool, v_rspace_new(ispin)%pw) 907 908 END DO 909 910 DEALLOCATE (v_rspace_new) 911 END IF 912 913 IF (ASSOCIATED(v_tau_rspace)) THEN 914 915 DO ispin = 1, nspins 916 917 CALL pw_transfer(v_tau_rspace(ispin)%pw, v_xc_gtemp%pw) 918 CALL pw_derive(v_xc_gtemp%pw, (/1, 0, 0/)) 919 CALL pw_transfer(v_xc_gtemp%pw, v_xc_rtemp%pw) 920 CALL qs_rho_get(nablavks_set(1, ispin)%rho, rho_r=rho_r) 921 CALL pw_axpy(v_xc_rtemp%pw, rho_r(1)%pw) 922 923 CALL pw_transfer(v_tau_rspace(ispin)%pw, v_xc_gtemp%pw) 924 CALL pw_derive(v_xc_gtemp%pw, (/0, 1, 0/)) 925 CALL pw_transfer(v_xc_gtemp%pw, v_xc_rtemp%pw) 926 CALL qs_rho_get(nablavks_set(2, ispin)%rho, rho_r=rho_r) 927 CALL pw_axpy(v_xc_rtemp%pw, rho_r(1)%pw) 928 929 CALL pw_transfer(v_tau_rspace(ispin)%pw, v_xc_gtemp%pw) 930 CALL pw_derive(v_xc_gtemp%pw, (/0, 0, 1/)) 931 CALL pw_transfer(v_xc_gtemp%pw, v_xc_rtemp%pw) 932 CALL qs_rho_get(nablavks_set(3, ispin)%rho, rho_r=rho_r) 933 CALL pw_axpy(v_xc_rtemp%pw, rho_r(1)%pw) 934 935 CALL pw_pool_give_back_pw(auxbas_pw_pool, v_tau_rspace(ispin)%pw) 936 937 END DO 938 939 DEALLOCATE (v_tau_rspace) 940 END IF 941 942 CALL pw_pool_give_back_pw(auxbas_pw_pool, v_xc_gtemp%pw) 943 CALL pw_pool_give_back_pw(auxbas_pw_pool, v_xc_rtemp%pw) 944 945 DEALLOCATE (v_xc_gtemp) 946 DEALLOCATE (v_xc_rtemp) 947 948 IF (gapw .OR. gapw_xc) THEN 949 CALL calculate_vxc_atom(qs_env=qs_env, energy_only=.FALSE., & 950 gradient_atom_set=nablavks_atom_set) 951 END IF 952 953! ------------------------------------- 954! Write Nabla V_KS (local) to cubes 955! ------------------------------------- 956 957 IF (BTEST(cp_print_key_should_output(logger%iter_info, lr_section, & 958 "EPR%PRINT%NABLAVKS_CUBES"), cp_p_file)) THEN 959 ALLOCATE (wf_r) 960 CALL pw_pool_create_pw(auxbas_pw_pool, wf_r%pw, & 961 use_data=REALDATA3D, & 962 in_space=REALSPACE) 963 DO idir = 1, 3 964 CALL pw_zero(wf_r%pw) 965 CALL qs_rho_get(nablavks_set(idir, 1)%rho, rho_r=rho_r) 966 CALL pw_copy(rho_r(1)%pw, wf_r%pw) ! RA 967 filename = "nablavks" 968 mpi_io = .TRUE. 969 WRITE (ext, '(a2,I1,a5)') "_d", idir, ".cube" 970 unit_nr = cp_print_key_unit_nr(logger, lr_section, "EPR%PRINT%NABLAVKS_CUBES", & 971 extension=TRIM(ext), middle_name=TRIM(filename), & 972 log_filename=.FALSE., file_position="REWIND", & 973 mpi_io=mpi_io) 974 CALL cp_pw_to_cube(wf_r%pw, unit_nr, "NABLA V_KS ", & 975 particles=particles, & 976 stride=section_get_ivals(lr_section, & 977 "EPR%PRINT%NABLAVKS_CUBES%STRIDE"), & 978 mpi_io=mpi_io) 979 CALL cp_print_key_finished_output(unit_nr, logger, lr_section, & 980 "EPR%PRINT%NABLAVKS_CUBES", mpi_io=mpi_io) 981 END DO 982 CALL pw_pool_give_back_pw(auxbas_pw_pool, wf_r%pw) 983 DEALLOCATE (wf_r) 984 END IF 985 986 CALL cp_print_key_finished_output(output_unit, logger, lr_section, & 987 "PRINT%PROGRAM_RUN_INFO") 988 989 END SUBROUTINE epr_nablavks 990 991END MODULE qs_linres_epr_nablavks 992 993