1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief The implicit (generalized) Poisson solver 8!> \par History 9!> 06.2014 created [Hossein Bani-Hashemian] 10!> 11.2015 - dealt with missing grid points of periodic grids while performing dct; 11!> - revised solver for Neumann and mixed boundary setups. 12!> \author Hossein Bani-Hashemian 13! ************************************************************************************************** 14MODULE ps_implicit_methods 15 USE bibliography, ONLY: BaniHashemian2016,& 16 cite_reference 17 USE cp_log_handling, ONLY: cp_get_default_logger,& 18 cp_logger_get_default_unit_nr,& 19 cp_logger_type 20 USE dct, ONLY: & 21 dct_type, dct_type_init, neumannX, neumannXY, neumannXYZ, neumannXZ, neumannY, neumannYZ, & 22 neumannZ, pw_expand, pw_shrink 23 USE dielectric_methods, ONLY: dielectric_create 24 USE dielectric_types, ONLY: dielectric_type 25 USE dirichlet_bc_methods, ONLY: dirichlet_boundary_region_setup 26 USE dirichlet_bc_types, ONLY: dbc_tile_release 27 USE kahan_sum, ONLY: accurate_sum 28 USE kinds, ONLY: dp,& 29 int_8 30 USE mathconstants, ONLY: fourpi,& 31 pi 32 USE message_passing, ONLY: mp_sum 33 USE ps_implicit_types, ONLY: MIXED_BC,& 34 MIXED_PERIODIC_BC,& 35 NEUMANN_BC,& 36 PERIODIC_BC,& 37 ps_implicit_type 38 USE pw_grid_types, ONLY: pw_grid_type 39 USE pw_methods, ONLY: pw_axpy,& 40 pw_copy,& 41 pw_derive,& 42 pw_integral_ab,& 43 pw_scale,& 44 pw_transfer,& 45 pw_zero 46 USE pw_poisson_types, ONLY: greens_fn_type,& 47 pw_poisson_parameter_type,& 48 pw_poisson_type 49 USE pw_pool_types, ONLY: pw_pool_create,& 50 pw_pool_create_pw,& 51 pw_pool_give_back_pw,& 52 pw_pool_release,& 53 pw_pool_type 54 USE pw_types, ONLY: COMPLEXDATA1D,& 55 REALDATA3D,& 56 REALSPACE,& 57 RECIPROCALSPACE,& 58 pw_p_type,& 59 pw_type 60#include "../base/base_uses.f90" 61 62 IMPLICIT NONE 63 PRIVATE 64 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ps_implicit_methods' 65 66 PUBLIC ps_implicit_create, & 67 implicit_poisson_solver_periodic, & 68 implicit_poisson_solver_neumann, & 69 implicit_poisson_solver_mixed_periodic, & 70 implicit_poisson_solver_mixed 71 72 INTERFACE ps_implicit_compute_ehartree 73 MODULE PROCEDURE compute_ehartree_periodic_bc, & 74 compute_ehartree_mixed_bc 75 END INTERFACE ps_implicit_compute_ehartree 76 77 REAL(dp), PRIVATE, PARAMETER :: large_error = 1.0E4_dp 78 79CONTAINS 80 81! ************************************************************************************************** 82!> \brief Creates implicit Poisson solver environment 83!> \param pw_pool pool of pw grid 84!> \param poisson_params poisson_env parameters 85!> \param dct_pw_grid discrete cosine transform (extended) grid 86!> \param green green function for FFT based inverse Laplacian 87!> \param ps_implicit_env implicit env to be created 88!> \par History 89!> 06.2014 created [Hossein Bani-Hashemian] 90!> \author Mohammad Hossein Bani-Hashemian 91! ************************************************************************************************** 92 SUBROUTINE ps_implicit_create(pw_pool, poisson_params, dct_pw_grid, green, ps_implicit_env) 93 94 TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool 95 TYPE(pw_poisson_parameter_type), INTENT(INOUT) :: poisson_params 96 TYPE(pw_grid_type), INTENT(IN), POINTER :: dct_pw_grid 97 TYPE(greens_fn_type), INTENT(IN), POINTER :: green 98 TYPE(ps_implicit_type), INTENT(INOUT), POINTER :: ps_implicit_env 99 100 CHARACTER(LEN=*), PARAMETER :: routineN = 'ps_implicit_create', & 101 routineP = moduleN//':'//routineN 102 103 INTEGER :: boundary_condition, handle, j, & 104 n_contacts, neumann_directions 105 TYPE(pw_pool_type), POINTER :: pw_pool_xpndd 106 107 CALL timeset(routineN, handle) 108 109 CALL cite_reference(BaniHashemian2016) 110 111 IF (.NOT. ASSOCIATED(ps_implicit_env)) THEN 112 ALLOCATE (ps_implicit_env) 113 114 ps_implicit_env%do_dbc_cube = poisson_params%dbc_params%do_dbc_cube 115 boundary_condition = poisson_params%ps_implicit_params%boundary_condition 116 neumann_directions = poisson_params%ps_implicit_params%neumann_directions 117 118! create dielectric 119 NULLIFY (ps_implicit_env%dielectric) 120 SELECT CASE (boundary_condition) 121 CASE (PERIODIC_BC, MIXED_PERIODIC_BC) 122 CALL dielectric_create(ps_implicit_env%dielectric, pw_pool, poisson_params%dielectric_params) 123 CASE (NEUMANN_BC, MIXED_BC) 124 CALL pw_pool_create(pw_pool_xpndd, pw_grid=dct_pw_grid) 125 CALL dielectric_create(ps_implicit_env%dielectric, pw_pool_xpndd, poisson_params%dielectric_params) 126 CALL pw_pool_release(pw_pool_xpndd) 127 END SELECT 128 129! initial guess 130 NULLIFY (ps_implicit_env%initial_guess) 131 132! v_eps 133 NULLIFY (ps_implicit_env%v_eps) 134 CALL pw_pool_create_pw(pw_pool, ps_implicit_env%v_eps, use_data=REALDATA3D, in_space=REALSPACE) 135 CALL pw_zero(ps_implicit_env%v_eps) 136 137! constraint charge 138 NULLIFY (ps_implicit_env%cstr_charge) 139 SELECT CASE (boundary_condition) 140 CASE (MIXED_PERIODIC_BC) 141 CALL pw_pool_create_pw(pw_pool, ps_implicit_env%cstr_charge, use_data=REALDATA3D, in_space=REALSPACE) 142 CALL pw_zero(ps_implicit_env%cstr_charge) 143 CASE (MIXED_BC) 144 CALL pw_pool_create(pw_pool_xpndd, pw_grid=dct_pw_grid) 145 CALL pw_pool_create_pw(pw_pool_xpndd, ps_implicit_env%cstr_charge, use_data=REALDATA3D, in_space=REALSPACE) 146 CALL pw_zero(ps_implicit_env%cstr_charge) 147 CALL pw_pool_release(pw_pool_xpndd) 148 END SELECT 149 150! initialize energies 151 ps_implicit_env%ehartree = 0.0_dp 152 ps_implicit_env%electric_enthalpy = 0.0_dp 153! times called 154 ps_implicit_env%times_called = 0 155 156! dct env 157 IF (boundary_condition .EQ. MIXED_BC .OR. boundary_condition .EQ. NEUMANN_BC) THEN 158 CALL dct_type_init(pw_pool%pw_grid, neumann_directions, ps_implicit_env%dct_env) 159 END IF 160 161! prepare dirichlet bc 162 CALL dirichlet_boundary_region_setup(pw_pool, poisson_params, ps_implicit_env%contacts) 163 CALL ps_implicit_prepare_blocks(pw_pool, dct_pw_grid, green, poisson_params, ps_implicit_env) 164 ! release tiles if they are not supposed to be written into cube files 165 IF ((boundary_condition .EQ. MIXED_PERIODIC_BC .OR. boundary_condition .EQ. MIXED_BC) .AND. & 166 (.NOT. poisson_params%dbc_params%do_dbc_cube)) THEN 167 n_contacts = SIZE(ps_implicit_env%contacts) 168 DO j = 1, n_contacts 169 CALL dbc_tile_release(ps_implicit_env%contacts(j)%dirichlet_bc, pw_pool) 170 END DO 171 END IF 172 173 END IF 174 175 CALL timestop(handle) 176 177 END SUBROUTINE ps_implicit_create 178 179! ************************************************************************************************** 180!> \brief implicit Poisson solver for periodic boundary conditions 181!> \param poisson_env poisson environment 182!> \param density electron density 183!> \param v_new electrostatic potential 184!> \param ehartree Hartree energy 185!> \par History 186!> 07.2014 created [Hossein Bani-Hashemian] 187!> \author Mohammad Hossein Bani-Hashemian 188! ************************************************************************************************** 189 SUBROUTINE implicit_poisson_solver_periodic(poisson_env, density, v_new, ehartree) 190 191 TYPE(pw_poisson_type), POINTER :: poisson_env 192 TYPE(pw_type), INTENT(IN), POINTER :: density 193 TYPE(pw_type), INTENT(INOUT), POINTER :: v_new 194 REAL(dp), INTENT(OUT), OPTIONAL :: ehartree 195 196 CHARACTER(LEN=*), PARAMETER :: routineN = 'implicit_poisson_solver_periodic', & 197 routineP = moduleN//':'//routineN 198 199 INTEGER :: handle, iter, max_iter, outp_unit, & 200 times_called 201 LOGICAL :: reached_max_iter, reached_tol, & 202 use_zero_initial_guess 203 REAL(dp) :: nabs_error, omega, pres_error, tol 204 TYPE(dielectric_type), POINTER :: dielectric 205 TYPE(greens_fn_type), POINTER :: green 206 TYPE(ps_implicit_type), POINTER :: ps_implicit_env 207 TYPE(pw_pool_type), POINTER :: pw_pool 208 TYPE(pw_type), POINTER :: g, PxQAinvxres, QAinvxres, res_new, & 209 res_old, v_old 210 211 CALL timeset(routineN, handle) 212 213 pw_pool => poisson_env%pw_pools(poisson_env%pw_level)%pool 214 dielectric => poisson_env%implicit_env%dielectric 215 green => poisson_env%green_fft 216 ps_implicit_env => poisson_env%implicit_env 217 218 tol = poisson_env%parameters%ps_implicit_params%tol 219 omega = poisson_env%parameters%ps_implicit_params%omega 220 max_iter = poisson_env%parameters%ps_implicit_params%max_iter 221 use_zero_initial_guess = poisson_env%parameters%ps_implicit_params%zero_initial_guess 222 times_called = ps_implicit_env%times_called 223 224! check if this is the first scf iteration 225 IF (times_called .EQ. 0) CALL ps_implicit_initial_guess_create(ps_implicit_env, pw_pool) 226 227 CALL pw_pool_create_pw(pw_pool, g, use_data=REALDATA3D, in_space=REALSPACE) 228 CALL pw_pool_create_pw(pw_pool, v_old, use_data=REALDATA3D, in_space=REALSPACE) 229 CALL pw_pool_create_pw(pw_pool, res_old, use_data=REALDATA3D, in_space=REALSPACE) 230 CALL pw_pool_create_pw(pw_pool, res_new, use_data=REALDATA3D, in_space=REALSPACE) 231 CALL pw_pool_create_pw(pw_pool, QAinvxres, use_data=REALDATA3D, in_space=REALSPACE) 232 CALL pw_pool_create_pw(pw_pool, PxQAinvxres, use_data=REALDATA3D, in_space=REALSPACE) 233 234 IF (use_zero_initial_guess) THEN 235 CALL pw_zero(v_old) 236 ELSE 237 CALL pw_copy(ps_implicit_env%initial_guess, v_old) 238 END IF 239 240 g%cr3d = fourpi*density%cr3d/dielectric%eps%cr3d 241 242! res_old = g - \Delta(v_old) - P(v_old) 243 CALL apply_poisson_operator_fft(pw_pool, green, dielectric, v_old, res_old) 244 CALL pw_scale(res_old, -1.0_dp) 245 CALL pw_axpy(g, res_old) 246 247! evaluate \Delta^-1(res_old) 248 CALL apply_inv_laplace_operator_fft(pw_pool, green, res_old, QAinvxres) 249 250 iter = 1 251 DO 252 253! v_new = v_old + \omega * QAinvxres_old 254 CALL pw_scale(QAinvxres, omega) 255 CALL pw_copy(QAinvxres, v_new) 256 CALL pw_axpy(v_old, v_new) 257 258! res_new = res_old - \omega * ( \Delta(QAinvxres_old) + P(QAinvxres_old) ) 259! = (1 - \omega) * res_old - \omega * PxQAinvxres 260 CALL apply_P_operator(pw_pool, dielectric, QAinvxres, PxQAinvxres) 261 CALL pw_copy(PxQAinvxres, res_new) 262 CALL pw_scale(res_new, -1.0_dp) 263 CALL pw_axpy(res_old, res_new, 1.0_dp - omega) 264 265! compute the error 266 CALL ps_implicit_compute_error_fft(pw_pool, green, res_new, v_old, v_new, QAinvxres, & 267 pres_error, nabs_error) 268! output 269 CALL ps_implicit_output(iter, pres_error, nabs_error, outp_unit) 270 IF (PRESENT(ehartree)) THEN 271 CALL ps_implicit_compute_ehartree(density, v_new, ehartree) 272 CALL ps_implicit_report_ehartree(ps_implicit_env, outp_unit, ehartree) 273 ps_implicit_env%ehartree = ehartree 274 ELSE 275 IF (outp_unit .GT. 0) WRITE (outp_unit, '(A1,/)') 276 END IF 277 278 iter = iter + 1 279 reached_max_iter = iter .GT. max_iter 280 reached_tol = pres_error .LE. tol 281 IF (pres_error .GT. large_error) & 282 CPABORT("Poisson solver did not converge.") 283 ps_implicit_env%times_called = ps_implicit_env%times_called + 1 284 IF (reached_max_iter .OR. reached_tol) EXIT 285 286! v_old = v_new, res_old = res_new 287 CALL pw_copy(v_new, v_old) 288 CALL pw_copy(res_new, res_old) 289 290 END DO 291 CALL ps_implicit_print_convergence_msg(iter, max_iter, outp_unit) 292 293 IF ((times_called .NE. 0) .AND. (.NOT. use_zero_initial_guess)) & 294 CALL pw_copy(v_new, ps_implicit_env%initial_guess) 295 296 IF (PRESENT(ehartree)) ehartree = ps_implicit_env%ehartree 297! compute the extra contribution to the Hamiltonian due to the presence of dielectric 298 CALL ps_implicit_compute_veps(pw_pool, dielectric, v_new, ps_implicit_env%v_eps) 299 300 CALL pw_pool_give_back_pw(pw_pool, g) 301 CALL pw_pool_give_back_pw(pw_pool, v_old) 302 CALL pw_pool_give_back_pw(pw_pool, res_old) 303 CALL pw_pool_give_back_pw(pw_pool, res_new) 304 CALL pw_pool_give_back_pw(pw_pool, QAinvxres) 305 CALL pw_pool_give_back_pw(pw_pool, PxQAinvxres) 306 307 CALL timestop(handle) 308 309 END SUBROUTINE implicit_poisson_solver_periodic 310 311! ************************************************************************************************** 312!> \brief implicit Poisson solver: zero-average solution of the Poisson equation 313!> subject to homogeneous Neumann boundary conditions 314!> \param poisson_env poisson environment 315!> \param density electron density 316!> \param v_new electrostatic potential 317!> \param ehartree Hartree energy 318!> \par History 319!> 02.2015 created [Hossein Bani-Hashemian] 320!> 11.2015 revised [Hossein Bani-Hashemian] 321!> \author Mohammad Hossein Bani-Hashemian 322! ************************************************************************************************** 323 SUBROUTINE implicit_poisson_solver_neumann(poisson_env, density, v_new, ehartree) 324 325 TYPE(pw_poisson_type), POINTER :: poisson_env 326 TYPE(pw_type), INTENT(IN), POINTER :: density 327 TYPE(pw_type), INTENT(INOUT), POINTER :: v_new 328 REAL(dp), INTENT(OUT), OPTIONAL :: ehartree 329 330 CHARACTER(LEN=*), PARAMETER :: routineN = 'implicit_poisson_solver_neumann', & 331 routineP = moduleN//':'//routineN 332 333 INTEGER :: handle, iter, max_iter, & 334 neumann_directions, outp_unit, & 335 times_called 336 LOGICAL :: reached_max_iter, reached_tol, & 337 use_zero_initial_guess 338 REAL(dp) :: nabs_error, omega, pres_error, tol, & 339 vol_scfac 340 TYPE(dct_type), POINTER :: dct_env 341 TYPE(dielectric_type), POINTER :: dielectric 342 TYPE(greens_fn_type), POINTER :: green 343 TYPE(ps_implicit_type), POINTER :: ps_implicit_env 344 TYPE(pw_pool_type), POINTER :: pw_pool, pw_pool_xpndd 345 TYPE(pw_type), POINTER :: density_xpndd, g, PxQAinvxres, & 346 QAinvxres, res_new, res_old, & 347 v_eps_xpndd, v_new_xpndd, v_old 348 349 CALL timeset(routineN, handle) 350 351 pw_pool => poisson_env%pw_pools(poisson_env%pw_level)%pool 352 dielectric => poisson_env%implicit_env%dielectric 353 green => poisson_env%green_fft 354 ps_implicit_env => poisson_env%implicit_env 355 dct_env => ps_implicit_env%dct_env 356 357 tol = poisson_env%parameters%ps_implicit_params%tol 358 omega = poisson_env%parameters%ps_implicit_params%omega 359 max_iter = poisson_env%parameters%ps_implicit_params%max_iter 360 use_zero_initial_guess = poisson_env%parameters%ps_implicit_params%zero_initial_guess 361 neumann_directions = poisson_env%parameters%ps_implicit_params%neumann_directions 362 times_called = ps_implicit_env%times_called 363 364 SELECT CASE (neumann_directions) 365 CASE (neumannXYZ) 366 vol_scfac = 8.0_dp 367 CASE (neumannXY, neumannXZ, neumannYZ) 368 vol_scfac = 4.0_dp 369 CASE (neumannX, neumannY, neumannZ) 370 vol_scfac = 2.0_dp 371 END SELECT 372 373 CALL pw_pool_create(pw_pool_xpndd, pw_grid=poisson_env%dct_pw_grid) 374 375! check if this is the first scf iteration 376 IF (times_called .EQ. 0) CALL ps_implicit_initial_guess_create(ps_implicit_env, pw_pool_xpndd) 377 378 CALL pw_pool_create_pw(pw_pool_xpndd, g, use_data=REALDATA3D, in_space=REALSPACE) 379 CALL pw_pool_create_pw(pw_pool_xpndd, v_old, use_data=REALDATA3D, in_space=REALSPACE) 380 CALL pw_pool_create_pw(pw_pool_xpndd, res_old, use_data=REALDATA3D, in_space=REALSPACE) 381 CALL pw_pool_create_pw(pw_pool_xpndd, res_new, use_data=REALDATA3D, in_space=REALSPACE) 382 CALL pw_pool_create_pw(pw_pool_xpndd, QAinvxres, use_data=REALDATA3D, in_space=REALSPACE) 383 CALL pw_pool_create_pw(pw_pool_xpndd, PxQAinvxres, use_data=REALDATA3D, in_space=REALSPACE) 384 CALL pw_pool_create_pw(pw_pool_xpndd, density_xpndd, use_data=REALDATA3D, in_space=REALSPACE) 385 CALL pw_pool_create_pw(pw_pool_xpndd, v_new_xpndd, use_data=REALDATA3D, in_space=REALSPACE) 386 CALL pw_pool_create_pw(pw_pool_xpndd, v_eps_xpndd, use_data=REALDATA3D, in_space=REALSPACE) 387 388 IF (use_zero_initial_guess) THEN 389 CALL pw_zero(v_old) 390 ELSE 391 CALL pw_copy(ps_implicit_env%initial_guess, v_old) 392 END IF 393 394 CALL pw_expand(neumann_directions, & 395 dct_env%recv_msgs_bnds, dct_env%dests_expand, dct_env%srcs_expand, & 396 dct_env%flipg_stat, dct_env%bounds_shftd, density, density_xpndd) 397 CALL pw_expand(neumann_directions, & 398 dct_env%recv_msgs_bnds, dct_env%dests_expand, dct_env%srcs_expand, & 399 dct_env%flipg_stat, dct_env%bounds_shftd, v_new, v_new_xpndd) 400 401 g%cr3d = fourpi*density_xpndd%cr3d/dielectric%eps%cr3d 402 403! res_old = g - \Delta(v_old) - P(v_old) 404 CALL apply_poisson_operator_dct(pw_pool_xpndd, green, dielectric, v_old, res_old) 405 CALL pw_scale(res_old, -1.0_dp) 406 CALL pw_axpy(g, res_old) 407 408! evaluate \Delta^-1(res_old) 409 CALL apply_inv_laplace_operator_dct(pw_pool_xpndd, green, res_old, QAinvxres) 410 411 iter = 1 412 DO 413 414! v_new = v_old + \omega * QAinvxres_old 415 CALL pw_scale(QAinvxres, omega) 416 CALL pw_copy(QAinvxres, v_new_xpndd) 417 CALL pw_axpy(v_old, v_new_xpndd) 418 419! res_new = res_old - \omega * ( \Delta(QAinvxres_old) + P(QAinvxres_old) ) 420! = (1 - \omega) * res_old - \omega * PxQAinvxres 421 CALL apply_P_operator(pw_pool_xpndd, dielectric, QAinvxres, PxQAinvxres) 422 CALL pw_copy(PxQAinvxres, res_new) 423 CALL pw_scale(res_new, -1.0_dp) 424 CALL pw_axpy(res_old, res_new, 1.0_dp - omega) 425 426! compute the error 427 CALL ps_implicit_compute_error_dct(pw_pool_xpndd, green, res_new, v_old, v_new_xpndd, QAinvxres, & 428 pres_error, nabs_error) 429! output 430 CALL ps_implicit_output(iter, pres_error, nabs_error, outp_unit) 431 IF (PRESENT(ehartree)) THEN 432 CALL ps_implicit_compute_ehartree(density_xpndd, v_new_xpndd, ehartree) 433 CALL ps_implicit_report_ehartree(ps_implicit_env, outp_unit, ehartree/vol_scfac) 434 ps_implicit_env%ehartree = ehartree/vol_scfac 435 ELSE 436 IF (outp_unit .GT. 0) WRITE (outp_unit, '(A1,/)') 437 END IF 438 439 iter = iter + 1 440 reached_max_iter = iter .GT. max_iter 441 reached_tol = pres_error .LE. tol 442 IF (pres_error .GT. large_error) & 443 CPABORT("Poisson solver did not converge.") 444 ps_implicit_env%times_called = ps_implicit_env%times_called + 1 445 IF (reached_max_iter .OR. reached_tol) EXIT 446 447! v_old = v_new, res_old = res_new 448 CALL pw_copy(v_new_xpndd, v_old) 449 CALL pw_copy(res_new, res_old) 450 451 END DO 452 CALL ps_implicit_print_convergence_msg(iter, max_iter, outp_unit) 453 454 CALL pw_shrink(neumann_directions, dct_env%dests_shrink, dct_env%srcs_shrink, & 455 dct_env%bounds_local_shftd, v_new_xpndd, v_new) 456 457 IF ((times_called .NE. 0) .AND. (.NOT. use_zero_initial_guess)) & 458 CALL pw_copy(v_new_xpndd, ps_implicit_env%initial_guess) 459 460 IF (PRESENT(ehartree)) ehartree = ps_implicit_env%ehartree 461! compute the extra contribution to the Hamiltonian due to the presence of dielectric 462! veps has to be computed for the expanded data and then shrinked otherwise we loose accuracy 463 CALL ps_implicit_compute_veps(pw_pool_xpndd, dielectric, v_new_xpndd, v_eps_xpndd) 464 CALL pw_shrink(neumann_directions, dct_env%dests_shrink, dct_env%srcs_shrink, & 465 dct_env%bounds_local_shftd, v_eps_xpndd, ps_implicit_env%v_eps) 466 467 CALL pw_pool_give_back_pw(pw_pool_xpndd, g) 468 CALL pw_pool_give_back_pw(pw_pool_xpndd, v_old) 469 CALL pw_pool_give_back_pw(pw_pool_xpndd, res_old) 470 CALL pw_pool_give_back_pw(pw_pool_xpndd, res_new) 471 CALL pw_pool_give_back_pw(pw_pool_xpndd, QAinvxres) 472 CALL pw_pool_give_back_pw(pw_pool_xpndd, PxQAinvxres) 473 CALL pw_pool_give_back_pw(pw_pool_xpndd, density_xpndd) 474 CALL pw_pool_give_back_pw(pw_pool_xpndd, v_new_xpndd) 475 CALL pw_pool_give_back_pw(pw_pool_xpndd, v_eps_xpndd) 476 CALL pw_pool_release(pw_pool_xpndd) 477 478 CALL timestop(handle) 479 480 END SUBROUTINE implicit_poisson_solver_neumann 481 482! ************************************************************************************************** 483!> \brief implicit Poisson solver for mixed-periodic boundary conditions (periodic + Dirichlet) 484!> \param poisson_env poisson environment 485!> \param density electron density 486!> \param v_new electrostatic potential 487!> \param electric_enthalpy electric enthalpy 488!> \par History 489!> 07.2014 created [Hossein Bani-Hashemian] 490!> \author Mohammad Hossein Bani-Hashemian 491! ************************************************************************************************** 492 SUBROUTINE implicit_poisson_solver_mixed_periodic(poisson_env, density, v_new, electric_enthalpy) 493 494 TYPE(pw_poisson_type), POINTER :: poisson_env 495 TYPE(pw_type), INTENT(IN), POINTER :: density 496 TYPE(pw_type), INTENT(INOUT), POINTER :: v_new 497 REAL(dp), INTENT(OUT), OPTIONAL :: electric_enthalpy 498 499 CHARACTER(LEN=*), PARAMETER :: routineN = 'implicit_poisson_solver_mixed_periodic', & 500 routineP = moduleN//':'//routineN 501 502 INTEGER :: data_size, handle, iter, j, lb1, lb2, lb3, max_iter, n_contacts, n_tiles_tot, ng, & 503 ngpts_local, nt, nt_tot, outp_unit, times_called, ub1, ub2, ub3 504 INTEGER(KIND=int_8) :: ngpts 505 INTEGER, DIMENSION(2, 3) :: bounds_local 506 INTEGER, DIMENSION(3) :: npts_local 507 LOGICAL :: reached_max_iter, reached_tol, & 508 use_zero_initial_guess 509 REAL(dp) :: Axvbar_avg, ehartree, eta, g_avg, & 510 gminusAxvbar_avg, nabs_error, omega, & 511 pres_error, tol 512 REAL(dp), ALLOCATABLE, DIMENSION(:) :: Btxlambda_new, Btxlambda_old, Bxv_bar, Bxv_new, & 513 lambda0, lambda_new, lambda_newNeta, lambda_old, QSxlambda, v_bar1D, v_D, v_new1D, w 514 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: B, Bt, QS, Rinv 515 REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: Btxlambda_new3D, Btxlambda_old3D 516 TYPE(dielectric_type), POINTER :: dielectric 517 TYPE(greens_fn_type), POINTER :: green 518 TYPE(ps_implicit_type), POINTER :: ps_implicit_env 519 TYPE(pw_grid_type), POINTER :: pw_grid 520 TYPE(pw_poisson_parameter_type), POINTER :: poisson_params 521 TYPE(pw_pool_type), POINTER :: pw_pool 522 TYPE(pw_type), POINTER :: Axvbar, g, PxQAinvxres, QAinvxres, & 523 res_new, res_old, v_old 524 525 CALL timeset(routineN, handle) 526 527 pw_pool => poisson_env%pw_pools(poisson_env%pw_level)%pool 528 pw_grid => pw_pool%pw_grid 529 poisson_params => poisson_env%parameters 530 dielectric => poisson_env%implicit_env%dielectric 531 green => poisson_env%green_fft 532 ps_implicit_env => poisson_env%implicit_env 533 534 ngpts_local = pw_grid%ngpts_local 535 ngpts = pw_grid%ngpts 536 npts_local = pw_grid%npts_local 537 bounds_local = pw_grid%bounds_local 538 tol = poisson_params%ps_implicit_params%tol 539 omega = poisson_params%ps_implicit_params%omega 540 max_iter = poisson_params%ps_implicit_params%max_iter 541 use_zero_initial_guess = poisson_params%ps_implicit_params%zero_initial_guess 542 times_called = ps_implicit_env%times_called 543 544 n_contacts = SIZE(ps_implicit_env%contacts) 545 n_tiles_tot = 0 546 DO j = 1, n_contacts 547 n_tiles_tot = n_tiles_tot + ps_implicit_env%contacts(j)%dirichlet_bc%n_tiles 548 END DO 549 550 IF (pw_grid%para%blocked) THEN 551 data_size = PRODUCT(npts_local) 552 ELSE IF (pw_grid%para%ray_distribution) THEN 553 data_size = ngpts_local 554 ELSE ! parallel run with np = 1 555 data_size = PRODUCT(npts_local) 556 END IF 557 558! check if this is the first scf iteration 559 IF (times_called .EQ. 0) CALL ps_implicit_initial_guess_create(ps_implicit_env, pw_pool) 560 561 ALLOCATE (B(n_tiles_tot, data_size)) 562 ALLOCATE (Bt(data_size, n_tiles_tot)) 563 ALLOCATE (QS(n_tiles_tot, n_tiles_tot)) 564 ALLOCATE (Rinv(n_tiles_tot + 1, n_tiles_tot + 1)) 565 566 B(:, :) = ps_implicit_env%B 567 Bt(:, :) = ps_implicit_env%Bt 568 QS(:, :) = ps_implicit_env%QS 569 Rinv(:, :) = ps_implicit_env%Rinv 570 CALL get_voltage(poisson_env%parameters%dbc_params%time, ps_implicit_env%v_D, ps_implicit_env%osc_frac, & 571 ps_implicit_env%frequency, ps_implicit_env%phase, v_D) 572 573 lb1 = bounds_local(1, 1); ub1 = bounds_local(2, 1) 574 lb2 = bounds_local(1, 2); ub2 = bounds_local(2, 2) 575 lb3 = bounds_local(1, 3); ub3 = bounds_local(2, 3) 576 577 ALLOCATE (lambda0(n_tiles_tot), lambda_old(n_tiles_tot), lambda_new(n_tiles_tot)) 578 ALLOCATE (Btxlambda_old(data_size), Btxlambda_new(data_size)) 579 ALLOCATE (Btxlambda_old3D(lb1:ub1, lb2:ub2, lb3:ub3), Btxlambda_new3D(lb1:ub1, lb2:ub2, lb3:ub3)) 580 ALLOCATE (QSxlambda(n_tiles_tot)) 581 ALLOCATE (w(n_tiles_tot + 1)) 582 ALLOCATE (lambda_newNeta(n_tiles_tot + 1)) 583 ALLOCATE (v_bar1D(data_size)) 584 ALLOCATE (Bxv_bar(n_tiles_tot)) 585 586 ALLOCATE (v_new1D(data_size)) 587 ALLOCATE (Bxv_new(n_tiles_tot)) 588 589 CALL pw_pool_create_pw(pw_pool, g, use_data=REALDATA3D, in_space=REALSPACE) 590 CALL pw_pool_create_pw(pw_pool, v_old, use_data=REALDATA3D, in_space=REALSPACE) 591 CALL pw_pool_create_pw(pw_pool, res_old, use_data=REALDATA3D, in_space=REALSPACE) 592 CALL pw_pool_create_pw(pw_pool, res_new, use_data=REALDATA3D, in_space=REALSPACE) 593 CALL pw_pool_create_pw(pw_pool, QAinvxres, use_data=REALDATA3D, in_space=REALSPACE) 594 CALL pw_pool_create_pw(pw_pool, PxQAinvxres, use_data=REALDATA3D, in_space=REALSPACE) 595 CALL pw_pool_create_pw(pw_pool, Axvbar, use_data=REALDATA3D, in_space=REALSPACE) 596 597 IF (use_zero_initial_guess) THEN 598 CALL pw_zero(v_old) 599 lambda0 = 0.0_dp 600 ELSE 601 CALL pw_copy(ps_implicit_env%initial_guess, v_old) 602 lambda0(:) = ps_implicit_env%initial_lambda 603 END IF 604 605 g%cr3d = fourpi*density%cr3d/dielectric%eps%cr3d 606 g_avg = accurate_sum(g%cr3d)/ngpts 607 608 lambda_old(:) = lambda0 609 610! res_old = g - \Delta(v_old) - P(v_old) - B^t * \lambda_old 611 CALL apply_poisson_operator_fft(pw_pool, green, dielectric, v_old, res_old) 612 CALL pw_scale(res_old, -1.0_dp) 613 CALL pw_axpy(g, res_old) 614 IF (data_size .NE. 0) THEN 615 CALL DGEMV('N', data_size, n_tiles_tot, 1.0_dp, Bt, data_size, lambda_old, 1, 0.0_dp, Btxlambda_old, 1) 616 END IF 617 CALL convert_1dto3d(ps_implicit_env%idx_1dto3d, Btxlambda_old, Btxlambda_old3D) 618 res_old%cr3d = res_old%cr3d-Btxlambda_old3D 619 620! evaluate \Delta^-1(res_old) 621 CALL apply_inv_laplace_operator_fft(pw_pool, green, res_old, QAinvxres) 622 623 iter = 1 624 DO 625 626! v_new (v_bar) = v_old + \omega * QAinvxres_old 627 CALL pw_scale(QAinvxres, omega) 628 CALL pw_copy(QAinvxres, v_new) 629 CALL pw_axpy(v_old, v_new) 630 631! evaluate 1^t * (g - \Delta(\bar{v}) - P(\bar{v})) 632! = 1^t * (g - P(\bar{v})) 633 CALL apply_P_operator(pw_pool, dielectric, v_new, Axvbar) 634 Axvbar_avg = accurate_sum(Axvbar%cr3d)/ngpts 635 gminusAxvbar_avg = g_avg - Axvbar_avg 636 CALL mp_sum(gminusAxvbar_avg, pw_grid%para%group) 637 638! evaluate Q_S * \lambda + v_D - B * \bar{v} 639 CALL DGEMV('N', n_tiles_tot, n_tiles_tot, 1.0_dp, QS, n_tiles_tot, lambda_old, 1, 0.0_dp, QSxlambda, 1) 640 v_bar1D(ps_implicit_env%idx_1dto3d) = RESHAPE(v_new%cr3d, (/data_size/)) 641 CALL DGEMV('N', n_tiles_tot, data_size, 1.0_dp, B, n_tiles_tot, v_bar1D, 1, 0.0_dp, Bxv_bar, 1) 642 CALL mp_sum(Bxv_bar, pw_grid%para%group) 643! solve R [\lambda; \eta] = [Q_S * \lambda + v_D - B * \bar{v}; 1^t * (g - \Delta(\bar{v}) - P(\bar{v}))] 644 w = 0.0_dp 645 w(:) = (/QSxlambda + v_D - Bxv_bar, gminusAxvbar_avg/) 646 CALL DGEMV('N', n_tiles_tot + 1, n_tiles_tot + 1, 1.0_dp, Rinv, n_tiles_tot + 1, w, 1, 0.0_dp, lambda_newNeta, 1) 647 lambda_new(:) = lambda_newNeta(1:n_tiles_tot) 648 eta = lambda_newNeta(n_tiles_tot + 1) 649 650! v_new = v_bar + 1 * \eta 651 v_new%cr3d = v_new%cr3d+eta/ngpts 652 653! evaluate B^t * \lambda_new 654 IF (data_size .NE. 0) THEN 655 CALL DGEMV('N', data_size, n_tiles_tot, 1.0_dp, Bt, data_size, lambda_new, 1, 0.0_dp, Btxlambda_new, 1) 656 END IF 657 CALL convert_1dto3d(ps_implicit_env%idx_1dto3d, Btxlambda_new, Btxlambda_new3D) 658 659! res_new = res_old - \omega * ( \Delta(QAinvxres_old) + P(QAinvxres_old) ) - B^t * ( \lambda_new - \lambda_old ) 660! = (1 - \omega) * res_old - \omega * P(QAinvxres_old) - B^t * ( \lambda_new - \lambda_old ) 661 CALL pw_zero(res_new) 662 CALL apply_P_operator(pw_pool, dielectric, QAinvxres, PxQAinvxres) 663 CALL pw_axpy(PxQAinvxres, res_new, -1.0_dp) 664 CALL pw_axpy(res_old, res_new, 1.0_dp - omega) 665 res_new%cr3d = res_new%cr3d+Btxlambda_old3D-Btxlambda_new3D 666 667! compute the error 668 CALL ps_implicit_compute_error_fft(pw_pool, green, res_new, v_old, v_new, QAinvxres, & 669 pres_error, nabs_error) 670! output 671 CALL ps_implicit_output(iter, pres_error, nabs_error, outp_unit) 672 IF (PRESENT(electric_enthalpy)) THEN 673 CALL ps_implicit_compute_ehartree(dielectric, density, Btxlambda_new3D, v_new, ehartree, electric_enthalpy) 674 CALL ps_implicit_report_ehartree(ps_implicit_env, outp_unit, ehartree) 675 ps_implicit_env%ehartree = ehartree 676 ps_implicit_env%electric_enthalpy = electric_enthalpy 677 ELSE 678 IF (outp_unit .GT. 0) WRITE (outp_unit, '(A1,/)') 679 END IF 680 681! verbose output 682 IF (poisson_params%dbc_params%verbose_output) THEN 683 v_new1D(ps_implicit_env%idx_1dto3d) = RESHAPE(v_new%cr3d, (/data_size/)) 684 CALL DGEMV('N', n_tiles_tot, data_size, 1.0_dp, B, n_tiles_tot, v_new1D, 1, 0.0_dp, Bxv_new, 1) 685 CALL mp_sum(Bxv_new, pw_grid%para%group) 686 IF (outp_unit .GT. 0) THEN 687 WRITE (outp_unit, '(T3,A,A)') "======== verbose ", REPEAT('=', 61) 688 WRITE (outp_unit, '(T20,A)') "Drgn tile vhartree lambda " 689 WRITE (outp_unit, '(T19,A)') REPEAT('-', 46) 690 nt_tot = 1 691 DO ng = 1, n_contacts 692 DO nt = 1, ps_implicit_env%contacts(ng)%dirichlet_bc%n_tiles 693 WRITE (outp_unit, '(T17,I6,5X,I6,3X,E13.4,E13.4)') ng, nt, Bxv_new(nt_tot), lambda_new(nt_tot) 694 nt_tot = nt_tot + 1 695 END DO 696 END DO 697 WRITE (outp_unit, '(T3,A)') REPEAT('=', 78) 698 END IF 699 END IF 700 701! check the convergence 702 iter = iter + 1 703 reached_max_iter = iter .GT. max_iter 704 reached_tol = pres_error .LE. tol 705 ps_implicit_env%times_called = ps_implicit_env%times_called + 1 706 IF (pres_error .GT. large_error) & 707 CPABORT("Poisson solver did not converge.") 708 IF (reached_max_iter .OR. reached_tol) EXIT 709 710! update 711 CALL pw_copy(v_new, v_old) 712 lambda_old(:) = lambda_new 713 CALL pw_copy(res_new, res_old) 714 Btxlambda_old3D(:, :, :) = Btxlambda_new3D 715 716 END DO 717 CALL ps_implicit_print_convergence_msg(iter, max_iter, outp_unit) 718 719 IF ((times_called .NE. 0) .AND. (.NOT. use_zero_initial_guess)) THEN 720 CALL pw_copy(v_new, ps_implicit_env%initial_guess) 721 ps_implicit_env%initial_lambda = lambda_new 722 END IF 723 724 ps_implicit_env%cstr_charge%cr3d = Btxlambda_new3D 725 IF (PRESENT(electric_enthalpy)) electric_enthalpy = ps_implicit_env%electric_enthalpy 726! compute the extra contribution to the Hamiltonian due to the presence of dielectric 727 CALL ps_implicit_compute_veps(pw_pool, dielectric, v_new, ps_implicit_env%v_eps) 728 729 CALL pw_pool_give_back_pw(pw_pool, g) 730 CALL pw_pool_give_back_pw(pw_pool, v_old) 731 CALL pw_pool_give_back_pw(pw_pool, res_old) 732 CALL pw_pool_give_back_pw(pw_pool, res_new) 733 CALL pw_pool_give_back_pw(pw_pool, QAinvxres) 734 CALL pw_pool_give_back_pw(pw_pool, PxQAinvxres) 735 CALL pw_pool_give_back_pw(pw_pool, Axvbar) 736 737 CALL timestop(handle) 738 739 END SUBROUTINE implicit_poisson_solver_mixed_periodic 740 741! ************************************************************************************************** 742!> \brief implicit Poisson solver for mixed boundary conditions (Neumann + Dirichlet) 743!> \param poisson_env poisson environment 744!> \param density electron density 745!> \param v_new electrostatic potential 746!> \param electric_enthalpy electric enthalpy 747!> \par History 748!> 10.2014 created [Hossein Bani-Hashemian] 749!> 11.2015 revised [Hossein Bani-Hashemian] 750!> \author Mohammad Hossein Bani-Hashemian 751! ************************************************************************************************** 752 SUBROUTINE implicit_poisson_solver_mixed(poisson_env, density, v_new, electric_enthalpy) 753 754 TYPE(pw_poisson_type), POINTER :: poisson_env 755 TYPE(pw_type), INTENT(IN), POINTER :: density 756 TYPE(pw_type), INTENT(INOUT), POINTER :: v_new 757 REAL(dp), INTENT(OUT), OPTIONAL :: electric_enthalpy 758 759 CHARACTER(LEN=*), PARAMETER :: routineN = 'implicit_poisson_solver_mixed', & 760 routineP = moduleN//':'//routineN 761 762 INTEGER :: data_size, handle, iter, j, lb1, lb2, lb3, max_iter, n_contacts, n_tiles_tot, & 763 neumann_directions, ng, ngpts_local, nt, nt_tot, outp_unit, times_called, ub1, ub2, ub3 764 INTEGER(KIND=int_8) :: ngpts 765 INTEGER, DIMENSION(2, 3) :: bounds_local 766 INTEGER, DIMENSION(3) :: npts_local 767 LOGICAL :: reached_max_iter, reached_tol, & 768 use_zero_initial_guess 769 REAL(dp) :: Axvbar_avg, ehartree, eta, g_avg, & 770 gminusAxvbar_avg, nabs_error, omega, & 771 pres_error, tol, vol_scfac 772 REAL(dp), ALLOCATABLE, DIMENSION(:) :: Btxlambda_new, Btxlambda_old, Bxv_bar, Bxv_new, & 773 lambda0, lambda_new, lambda_newNeta, lambda_old, QSxlambda, v_bar1D, v_D, v_new1D, w 774 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: B, Bt, QS, Rinv 775 REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: Btxlambda_new3D, Btxlambda_old3D 776 TYPE(dct_type), POINTER :: dct_env 777 TYPE(dielectric_type), POINTER :: dielectric 778 TYPE(greens_fn_type), POINTER :: green 779 TYPE(ps_implicit_type), POINTER :: ps_implicit_env 780 TYPE(pw_grid_type), POINTER :: dct_pw_grid, pw_grid 781 TYPE(pw_poisson_parameter_type), POINTER :: poisson_params 782 TYPE(pw_pool_type), POINTER :: pw_pool, pw_pool_xpndd 783 TYPE(pw_type), POINTER :: Axvbar, density_xpndd, g, PxQAinvxres, & 784 QAinvxres, res_new, res_old, & 785 v_eps_xpndd, v_new_xpndd, v_old 786 787 CALL timeset(routineN, handle) 788 789 pw_pool => poisson_env%pw_pools(poisson_env%pw_level)%pool 790 pw_grid => pw_pool%pw_grid 791 poisson_params => poisson_env%parameters 792 dielectric => poisson_env%implicit_env%dielectric 793 green => poisson_env%green_fft 794 ps_implicit_env => poisson_env%implicit_env 795 dct_env => ps_implicit_env%dct_env 796 797 dct_pw_grid => poisson_env%dct_pw_grid 798 ngpts_local = dct_pw_grid%ngpts_local 799 ngpts = dct_pw_grid%ngpts 800 npts_local = dct_pw_grid%npts_local 801 bounds_local = dct_pw_grid%bounds_local 802 tol = poisson_params%ps_implicit_params%tol 803 omega = poisson_params%ps_implicit_params%omega 804 max_iter = poisson_params%ps_implicit_params%max_iter 805 use_zero_initial_guess = poisson_params%ps_implicit_params%zero_initial_guess 806 neumann_directions = poisson_params%ps_implicit_params%neumann_directions 807 times_called = ps_implicit_env%times_called 808 809 SELECT CASE (neumann_directions) 810 CASE (neumannXYZ) 811 vol_scfac = 8.0_dp 812 CASE (neumannXY, neumannXZ, neumannYZ) 813 vol_scfac = 4.0_dp 814 CASE (neumannX, neumannY, neumannZ) 815 vol_scfac = 2.0_dp 816 END SELECT 817 818 n_contacts = SIZE(ps_implicit_env%contacts) 819 n_tiles_tot = 0 820 DO j = 1, n_contacts 821 n_tiles_tot = n_tiles_tot + ps_implicit_env%contacts(j)%dirichlet_bc%n_tiles 822 END DO 823 824 IF (dct_pw_grid%para%blocked) THEN 825 data_size = PRODUCT(npts_local) 826 ELSE IF (dct_pw_grid%para%ray_distribution) THEN 827 data_size = ngpts_local 828 ELSE ! parallel run with np = 1 829 data_size = PRODUCT(npts_local) 830 END IF 831 832 CALL pw_pool_create(pw_pool_xpndd, pw_grid=dct_pw_grid) 833 834! check if this is the first scf iteration 835 IF (times_called .EQ. 0) CALL ps_implicit_initial_guess_create(ps_implicit_env, pw_pool_xpndd) 836 837 ALLOCATE (B(n_tiles_tot, data_size)) 838 ALLOCATE (Bt(data_size, n_tiles_tot)) 839 ALLOCATE (QS(n_tiles_tot, n_tiles_tot)) 840 ALLOCATE (Rinv(n_tiles_tot + 1, n_tiles_tot + 1)) 841 842 B(:, :) = ps_implicit_env%B 843 Bt(:, :) = ps_implicit_env%Bt 844 QS(:, :) = ps_implicit_env%QS 845 Rinv(:, :) = ps_implicit_env%Rinv 846 CALL get_voltage(poisson_env%parameters%dbc_params%time, ps_implicit_env%v_D, ps_implicit_env%osc_frac, & 847 ps_implicit_env%frequency, ps_implicit_env%phase, v_D) 848 849 lb1 = bounds_local(1, 1); ub1 = bounds_local(2, 1) 850 lb2 = bounds_local(1, 2); ub2 = bounds_local(2, 2) 851 lb3 = bounds_local(1, 3); ub3 = bounds_local(2, 3) 852 853 ALLOCATE (lambda0(n_tiles_tot), lambda_old(n_tiles_tot), lambda_new(n_tiles_tot)) 854 ALLOCATE (Btxlambda_old(data_size), Btxlambda_new(data_size)) 855 ALLOCATE (Btxlambda_old3D(lb1:ub1, lb2:ub2, lb3:ub3), Btxlambda_new3D(lb1:ub1, lb2:ub2, lb3:ub3)) 856 ALLOCATE (QSxlambda(n_tiles_tot)) 857 ALLOCATE (w(n_tiles_tot + 1)) 858 ALLOCATE (lambda_newNeta(n_tiles_tot + 1)) 859 ALLOCATE (v_bar1D(data_size)) 860 ALLOCATE (Bxv_bar(n_tiles_tot)) 861 862 ALLOCATE (v_new1D(data_size)) 863 ALLOCATE (Bxv_new(n_tiles_tot)) 864 865 CALL pw_pool_create_pw(pw_pool_xpndd, g, use_data=REALDATA3D, in_space=REALSPACE) 866 CALL pw_pool_create_pw(pw_pool_xpndd, v_old, use_data=REALDATA3D, in_space=REALSPACE) 867 CALL pw_pool_create_pw(pw_pool_xpndd, res_old, use_data=REALDATA3D, in_space=REALSPACE) 868 CALL pw_pool_create_pw(pw_pool_xpndd, res_new, use_data=REALDATA3D, in_space=REALSPACE) 869 CALL pw_pool_create_pw(pw_pool_xpndd, QAinvxres, use_data=REALDATA3D, in_space=REALSPACE) 870 CALL pw_pool_create_pw(pw_pool_xpndd, PxQAinvxres, use_data=REALDATA3D, in_space=REALSPACE) 871 CALL pw_pool_create_pw(pw_pool_xpndd, Axvbar, use_data=REALDATA3D, in_space=REALSPACE) 872 CALL pw_pool_create_pw(pw_pool_xpndd, density_xpndd, use_data=REALDATA3D, in_space=REALSPACE) 873 CALL pw_pool_create_pw(pw_pool_xpndd, v_new_xpndd, use_data=REALDATA3D, in_space=REALSPACE) 874 CALL pw_pool_create_pw(pw_pool_xpndd, v_eps_xpndd, use_data=REALDATA3D, in_space=REALSPACE) 875 876 IF (use_zero_initial_guess) THEN 877 CALL pw_zero(v_old) 878 lambda0 = 0.0_dp 879 ELSE 880 CALL pw_copy(ps_implicit_env%initial_guess, v_old) 881 lambda0(:) = ps_implicit_env%initial_lambda 882 END IF 883 884 CALL pw_expand(neumann_directions, & 885 dct_env%recv_msgs_bnds, dct_env%dests_expand, dct_env%srcs_expand, & 886 dct_env%flipg_stat, dct_env%bounds_shftd, density, density_xpndd) 887 CALL pw_expand(neumann_directions, & 888 dct_env%recv_msgs_bnds, dct_env%dests_expand, dct_env%srcs_expand, & 889 dct_env%flipg_stat, dct_env%bounds_shftd, v_new, v_new_xpndd) 890 891 g%cr3d = fourpi*density_xpndd%cr3d/dielectric%eps%cr3d 892 g_avg = accurate_sum(g%cr3d)/ngpts 893 894 lambda_old(:) = lambda0 895 896! res_old = g - \Delta(v_old) - P(v_old) - B^t * \lambda_old 897 CALL apply_poisson_operator_dct(pw_pool_xpndd, green, dielectric, v_old, res_old) 898 CALL pw_scale(res_old, -1.0_dp) 899 CALL pw_axpy(g, res_old) 900 IF (data_size .NE. 0) THEN 901 CALL DGEMV('N', data_size, n_tiles_tot, 1.0_dp, Bt, data_size, lambda_old, 1, 0.0_dp, Btxlambda_old, 1) 902 END IF 903 CALL convert_1dto3d(ps_implicit_env%idx_1dto3d, Btxlambda_old, Btxlambda_old3D) 904 res_old%cr3d = res_old%cr3d-Btxlambda_old3D 905 906! evaluate \Delta^-1(res_old) 907 CALL apply_inv_laplace_operator_dct(pw_pool_xpndd, green, res_old, QAinvxres) 908 909 iter = 1 910 DO 911 912! v_new (v_bar) = v_old + \omega * QAinvxres_old 913 CALL pw_scale(QAinvxres, omega) 914 CALL pw_copy(QAinvxres, v_new_xpndd) 915 CALL pw_axpy(v_old, v_new_xpndd) 916 917! evaluate 1^t * (g - \Delta(\bar{v}) - P(\bar{v})) 918! = 1^t * (g - P(\bar{v})) 919 CALL apply_P_operator(pw_pool_xpndd, dielectric, v_new_xpndd, Axvbar) 920 Axvbar_avg = accurate_sum(Axvbar%cr3d)/ngpts 921 gminusAxvbar_avg = g_avg - Axvbar_avg 922 CALL mp_sum(gminusAxvbar_avg, dct_pw_grid%para%group) 923 924! evaluate Q_S * \lambda + v_D - B * \bar{v} 925 CALL DGEMV('N', n_tiles_tot, n_tiles_tot, 1.0_dp, QS, n_tiles_tot, lambda_old, 1, 0.0_dp, QSxlambda, 1) 926 v_bar1D(ps_implicit_env%idx_1dto3d) = RESHAPE(v_new_xpndd%cr3d, (/data_size/)) 927 CALL DGEMV('N', n_tiles_tot, data_size, 1.0_dp, B, n_tiles_tot, v_bar1D, 1, 0.0_dp, Bxv_bar, 1) 928 CALL mp_sum(Bxv_bar, dct_pw_grid%para%group) 929! solve R [\lambda; \eta] = [Q_S * \lambda + v_D - B * \bar{v}; 1^t * (g - \Delta(\bar{v}) - P(\bar{v}))] 930 w = 0.0_dp 931 w(:) = (/QSxlambda + v_D - Bxv_bar, gminusAxvbar_avg/) 932 CALL DGEMV('N', n_tiles_tot + 1, n_tiles_tot + 1, 1.0_dp, Rinv, n_tiles_tot + 1, w, 1, 0.0_dp, lambda_newNeta, 1) 933 lambda_new(:) = lambda_newNeta(1:n_tiles_tot) 934 eta = lambda_newNeta(n_tiles_tot + 1) 935 936! v_new = v_bar + 1 * \eta 937 v_new_xpndd%cr3d = v_new_xpndd%cr3d+eta/ngpts 938 939! evaluate B^t * \lambda_new 940 IF (data_size .NE. 0) THEN 941 CALL DGEMV('N', data_size, n_tiles_tot, 1.0_dp, Bt, data_size, lambda_new, 1, 0.0_dp, Btxlambda_new, 1) 942 END IF 943 CALL convert_1dto3d(ps_implicit_env%idx_1dto3d, Btxlambda_new, Btxlambda_new3D) 944 945! res_new = res_old - \omega * ( \Delta(QAinvxres_old) + P(QAinvxres_old) ) - B^t * ( \lambda_new - \lambda_old ) 946! = (1 - \omega) * res_old - \omega * P(QAinvxres_old) - B^t * ( \lambda_new - \lambda_old ) 947 CALL pw_zero(res_new) 948 CALL apply_P_operator(pw_pool_xpndd, dielectric, QAinvxres, PxQAinvxres) 949 CALL pw_axpy(PxQAinvxres, res_new, -1.0_dp) 950 CALL pw_axpy(res_old, res_new, 1.0_dp - omega) 951 res_new%cr3d = res_new%cr3d-Btxlambda_new3D+Btxlambda_old3D 952 953! compute the error 954 CALL ps_implicit_compute_error_dct(pw_pool_xpndd, green, res_new, v_old, v_new_xpndd, QAinvxres, & 955 pres_error, nabs_error) 956! output 957 CALL ps_implicit_output(iter, pres_error, nabs_error, outp_unit) 958 IF (PRESENT(electric_enthalpy)) THEN 959 CALL ps_implicit_compute_ehartree(dielectric, density_xpndd, Btxlambda_new3D, v_new_xpndd, & 960 ehartree, electric_enthalpy) 961 CALL ps_implicit_report_ehartree(ps_implicit_env, outp_unit, ehartree/vol_scfac) 962 ps_implicit_env%ehartree = ehartree/vol_scfac 963 ps_implicit_env%electric_enthalpy = electric_enthalpy/vol_scfac 964 ELSE 965 IF (outp_unit .GT. 0) WRITE (outp_unit, '(A1,/)') 966 END IF 967 968! verbose output 969 IF (poisson_params%dbc_params%verbose_output) THEN 970 v_new1D(ps_implicit_env%idx_1dto3d) = RESHAPE(v_new_xpndd%cr3d, (/data_size/)) 971 CALL DGEMV('N', n_tiles_tot, data_size, 1.0_dp, B, n_tiles_tot, v_new1D, 1, 0.0_dp, Bxv_new, 1) 972 CALL mp_sum(Bxv_new, pw_grid%para%group) 973 IF (outp_unit .GT. 0) THEN 974 WRITE (outp_unit, '(T3,A)') "======== verbose "//REPEAT('=', 61) 975 WRITE (outp_unit, '(T20,A)') "Drgn tile vhartree lambda " 976 WRITE (outp_unit, '(T19,A)') REPEAT('-', 46) 977 nt_tot = 1 978 DO ng = 1, n_contacts 979 DO nt = 1, ps_implicit_env%contacts(ng)%dirichlet_bc%n_tiles 980 WRITE (outp_unit, '(T17,I6,5X,I6,3X,E13.4,E13.4)') ng, nt, Bxv_new(nt_tot), lambda_new(nt_tot) 981 nt_tot = nt_tot + 1 982 END DO 983 END DO 984 WRITE (outp_unit, '(T3,A)') REPEAT('=', 78) 985 END IF 986 END IF 987 988! check the convergence 989 iter = iter + 1 990 reached_max_iter = iter .GT. max_iter 991 reached_tol = pres_error .LE. tol 992 ps_implicit_env%times_called = ps_implicit_env%times_called + 1 993 IF (pres_error .GT. large_error) & 994 CPABORT("Poisson solver did not converge.") 995 IF (reached_max_iter .OR. reached_tol) EXIT 996 997! update 998 CALL pw_copy(v_new_xpndd, v_old) 999 lambda_old(:) = lambda_new 1000 CALL pw_copy(res_new, res_old) 1001 Btxlambda_old3D(:, :, :) = Btxlambda_new3D 1002 1003 END DO 1004 CALL ps_implicit_print_convergence_msg(iter, max_iter, outp_unit) 1005 1006 CALL pw_shrink(neumann_directions, dct_env%dests_shrink, dct_env%srcs_shrink, & 1007 dct_env%bounds_local_shftd, v_new_xpndd, v_new) 1008 1009 IF ((times_called .NE. 0) .AND. (.NOT. use_zero_initial_guess)) THEN 1010 CALL pw_copy(v_new_xpndd, ps_implicit_env%initial_guess) 1011 ps_implicit_env%initial_lambda = lambda_new 1012 END IF 1013 1014 ps_implicit_env%cstr_charge%cr3d = Btxlambda_new3D 1015 IF (PRESENT(electric_enthalpy)) electric_enthalpy = ps_implicit_env%electric_enthalpy 1016! compute the extra contribution to the Hamiltonian due to the presence of dielectric 1017 CALL ps_implicit_compute_veps(pw_pool_xpndd, dielectric, v_new_xpndd, v_eps_xpndd) 1018 CALL pw_shrink(neumann_directions, dct_env%dests_shrink, dct_env%srcs_shrink, & 1019 dct_env%bounds_local_shftd, v_eps_xpndd, ps_implicit_env%v_eps) 1020 1021 CALL pw_pool_give_back_pw(pw_pool_xpndd, g) 1022 CALL pw_pool_give_back_pw(pw_pool_xpndd, v_old) 1023 CALL pw_pool_give_back_pw(pw_pool_xpndd, res_old) 1024 CALL pw_pool_give_back_pw(pw_pool_xpndd, res_new) 1025 CALL pw_pool_give_back_pw(pw_pool_xpndd, QAinvxres) 1026 CALL pw_pool_give_back_pw(pw_pool_xpndd, PxQAinvxres) 1027 CALL pw_pool_give_back_pw(pw_pool_xpndd, Axvbar) 1028 CALL pw_pool_give_back_pw(pw_pool_xpndd, density_xpndd) 1029 CALL pw_pool_give_back_pw(pw_pool_xpndd, v_new_xpndd) 1030 CALL pw_pool_give_back_pw(pw_pool_xpndd, v_eps_xpndd) 1031 CALL pw_pool_release(pw_pool_xpndd) 1032 1033 CALL timestop(handle) 1034 1035 END SUBROUTINE implicit_poisson_solver_mixed 1036 1037! ************************************************************************************************** 1038!> \brief allocates and zeroises initial guess for implicit (iterative) Poisson solver 1039!> \param ps_implicit_env the implicit env contaning the initial guess 1040!> \param pw_pool pool of pw grid 1041!> \par History 1042!> 06.2014 created [Hossein Bani-Hashemian] 1043!> \author Mohammad Hossein Bani-Hashemian 1044! ************************************************************************************************** 1045 SUBROUTINE ps_implicit_initial_guess_create(ps_implicit_env, pw_pool) 1046 1047 TYPE(ps_implicit_type), INTENT(INOUT), POINTER :: ps_implicit_env 1048 TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool 1049 1050 CHARACTER(LEN=*), PARAMETER :: routineN = 'ps_implicit_initial_guess_create', & 1051 routineP = moduleN//':'//routineN 1052 1053 INTEGER :: handle, n_tiles_tot 1054 1055 CALL timeset(routineN, handle) 1056 1057 n_tiles_tot = SIZE(ps_implicit_env%v_D) 1058 CALL pw_pool_create_pw(pw_pool, ps_implicit_env%initial_guess, & 1059 use_data=REALDATA3D, in_space=REALSPACE) 1060 CALL pw_zero(ps_implicit_env%initial_guess) 1061 ALLOCATE (ps_implicit_env%initial_lambda(n_tiles_tot)) 1062 ps_implicit_env%initial_lambda = 0.0_dp 1063 1064 CALL timestop(handle) 1065 1066 END SUBROUTINE ps_implicit_initial_guess_create 1067 1068! ************************************************************************************************** 1069!> \brief prepare blocks B, Bt, QS, R^-1, v_D 1070!> \param pw_pool_orig original pw grid 1071!> \param dct_pw_grid DCT (extended) grid 1072!> \param green green functions for FFT based inverse Laplacian 1073!> \param poisson_params paramaters of the poisson_env 1074!> \param ps_implicit_env the implicit_env that stores the blocks 1075!> \par History 1076!> 10.2014 created [Hossein Bani-Hashemian] 1077!> \author Mohammad Hossein Bani-Hashemian 1078! ************************************************************************************************** 1079 SUBROUTINE ps_implicit_prepare_blocks(pw_pool_orig, dct_pw_grid, green, & 1080 poisson_params, ps_implicit_env) 1081 1082 TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool_orig 1083 TYPE(pw_grid_type), INTENT(IN), POINTER :: dct_pw_grid 1084 TYPE(greens_fn_type), INTENT(IN), POINTER :: green 1085 TYPE(pw_poisson_parameter_type), INTENT(IN) :: poisson_params 1086 TYPE(ps_implicit_type), INTENT(INOUT), POINTER :: ps_implicit_env 1087 1088 CHARACTER(LEN=*), PARAMETER :: routineN = 'ps_implicit_prepare_blocks', & 1089 routineP = moduleN//':'//routineN 1090 1091 INTEGER :: data_size, handle, i, indx1, indx2, info, j, k, l, lb1, lb2, lb3, n_contacts, & 1092 n_tiles, n_tiles_tot, neumann_directions, ngpts_local, ub1, ub2, ub3, unit_nr 1093 INTEGER(KIND=int_8) :: ngpts 1094 INTEGER, ALLOCATABLE, DIMENSION(:) :: ipiv 1095 INTEGER, DIMENSION(2, 3) :: bounds, bounds_local 1096 INTEGER, DIMENSION(3) :: npts, npts_local 1097 LOGICAL :: done_preparing 1098 REAL(dp) :: tile_volume, vol_scfac 1099 REAL(dp), ALLOCATABLE, DIMENSION(:) :: Bxunit_vec, test_vec, work_arr 1100 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: QAinvxBt, R 1101 TYPE(cp_logger_type), POINTER :: logger 1102 TYPE(dct_type), POINTER :: dct_env 1103 TYPE(pw_grid_type), POINTER :: pw_grid_orig 1104 TYPE(pw_pool_type), POINTER :: pw_pool_xpndd 1105 TYPE(pw_type), POINTER :: pw_in, pw_out, tile_pw 1106 1107 CALL timeset(routineN, handle) 1108 1109 pw_grid_orig => pw_pool_orig%pw_grid 1110 1111 logger => cp_get_default_logger() 1112 IF (logger%para_env%mepos .EQ. logger%para_env%source) THEN 1113 unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.) 1114 ELSE 1115 unit_nr = -1 1116 ENDIF 1117 1118 SELECT CASE (poisson_params%ps_implicit_params%boundary_condition) 1119 CASE (MIXED_BC) 1120 1121 ngpts_local = dct_pw_grid%ngpts_local 1122 ngpts = dct_pw_grid%ngpts 1123 npts_local = dct_pw_grid%npts_local 1124 npts = dct_pw_grid%npts 1125 bounds_local = dct_pw_grid%bounds_local 1126 bounds = dct_pw_grid%bounds 1127 dct_env => ps_implicit_env%dct_env 1128 1129 neumann_directions = poisson_params%ps_implicit_params%neumann_directions 1130 1131 SELECT CASE (neumann_directions) 1132 CASE (neumannXYZ) 1133 vol_scfac = 8.0_dp 1134 CASE (neumannXY, neumannXZ, neumannYZ) 1135 vol_scfac = 4.0_dp 1136 CASE (neumannX, neumannY, neumannZ) 1137 vol_scfac = 2.0_dp 1138 END SELECT 1139 1140! evaluate indices for converting 3D arrays into 1D arrays 1141 lb1 = bounds_local(1, 1); ub1 = bounds_local(2, 1) 1142 lb2 = bounds_local(1, 2); ub2 = bounds_local(2, 2) 1143 lb3 = bounds_local(1, 3); ub3 = bounds_local(2, 3) 1144 1145 IF (dct_pw_grid%para%blocked) THEN 1146 data_size = PRODUCT(npts_local) 1147 ELSE IF (dct_pw_grid%para%ray_distribution) THEN 1148 data_size = ngpts_local 1149 ELSE ! parallel run with np = 1 1150 data_size = PRODUCT(npts_local) 1151 END IF 1152 1153 ALLOCATE (ps_implicit_env%idx_1dto3d(data_size)) 1154 l = 1 1155 DO k = lb3, ub3 1156 DO j = lb2, ub2 1157 DO i = lb1, ub1 1158 ps_implicit_env%idx_1dto3d(l) = (i - lb1 + 1) + & 1159 (j - lb2)*npts_local(1) + & 1160 (k - lb3)*npts_local(1)*npts_local(2) 1161 l = l + 1 1162 END DO 1163 END DO 1164 END DO 1165 1166 n_contacts = SIZE(ps_implicit_env%contacts) 1167 n_tiles_tot = 0 1168 DO j = 1, n_contacts 1169 n_tiles_tot = n_tiles_tot + ps_implicit_env%contacts(j)%dirichlet_bc%n_tiles 1170 END DO 1171 1172 ALLOCATE (ps_implicit_env%B(n_tiles_tot, data_size)) 1173 ALLOCATE (ps_implicit_env%Bt(data_size, n_tiles_tot)) 1174 ALLOCATE (ps_implicit_env%QS(n_tiles_tot, n_tiles_tot)) 1175 ALLOCATE (ps_implicit_env%Rinv(n_tiles_tot + 1, n_tiles_tot + 1)) 1176 ALLOCATE (ps_implicit_env%v_D(n_tiles_tot)) 1177 ALLOCATE (ps_implicit_env%osc_frac(n_tiles_tot)) 1178 ALLOCATE (ps_implicit_env%frequency(n_tiles_tot)) 1179 ALLOCATE (ps_implicit_env%phase(n_tiles_tot)) 1180 1181 ALLOCATE (QAinvxBt(data_size, n_tiles_tot)) 1182 ALLOCATE (Bxunit_vec(n_tiles_tot)) 1183 ALLOCATE (test_vec(n_tiles_tot)) 1184 ALLOCATE (R(n_tiles_tot + 1, n_tiles_tot + 1)) 1185 ALLOCATE (work_arr(n_tiles_tot + 1), ipiv(n_tiles_tot + 1)) ! LAPACK work and ipiv arrays 1186 1187! prepare pw_pool for evaluating inverse Laplacian of tile_pw's using DCT 1188 CALL pw_pool_create(pw_pool_xpndd, pw_grid=dct_pw_grid) 1189 1190! set up B, B^t, (\Delta^-1)*B^t 1191 indx1 = 1 1192 DO j = 1, n_contacts 1193 n_tiles = ps_implicit_env%contacts(j)%dirichlet_bc%n_tiles 1194 indx2 = indx1 + n_tiles - 1 1195 DO i = 1, n_tiles 1196 tile_pw => ps_implicit_env%contacts(j)%dirichlet_bc%tiles(i)%tile%tile_pw 1197 1198 CALL pw_pool_create_pw(pw_pool_xpndd, pw_in, use_data=REALDATA3D, in_space=REALSPACE) 1199 CALL pw_expand(neumann_directions, & 1200 dct_env%recv_msgs_bnds, dct_env%dests_expand, dct_env%srcs_expand, & 1201 dct_env%flipg_stat, dct_env%bounds_shftd, tile_pw, pw_in) 1202 1203 tile_volume = ps_implicit_env%contacts(j)%dirichlet_bc%tiles(i)%tile%volume 1204 CALL pw_scale(pw_in, 1.0_dp/(vol_scfac*tile_volume)) ! normalize tile_pw 1205 ps_implicit_env%Bt(ps_implicit_env%idx_1dto3d, indx1 + i - 1) = RESHAPE(pw_in%cr3d, (/data_size/)) 1206 1207 CALL pw_pool_create_pw(pw_pool_xpndd, pw_out, use_data=REALDATA3D, in_space=REALSPACE) 1208 CALL apply_inv_laplace_operator_dct(pw_pool_xpndd, green, pw_in, pw_out) 1209 QAinvxBt(ps_implicit_env%idx_1dto3d, indx1 + i - 1) = RESHAPE(pw_out%cr3d, (/data_size/)) 1210 ! the electrostatic potential has opposite sign by internal convention 1211 ps_implicit_env%v_D(indx1 + i - 1) = -1.0_dp*ps_implicit_env%contacts(j)%dirichlet_bc%v_D 1212 ps_implicit_env%osc_frac(indx1 + i - 1) = ps_implicit_env%contacts(j)%dirichlet_bc%osc_frac 1213 ps_implicit_env%frequency(indx1 + i - 1) = ps_implicit_env%contacts(j)%dirichlet_bc%frequency 1214 ps_implicit_env%phase(indx1 + i - 1) = ps_implicit_env%contacts(j)%dirichlet_bc%phase 1215 1216 CALL pw_pool_give_back_pw(pw_pool_xpndd, pw_in) 1217 CALL pw_pool_give_back_pw(pw_pool_xpndd, pw_out) 1218 END DO 1219 indx1 = indx2 + 1 1220 END DO 1221 ps_implicit_env%B = TRANSPOSE(ps_implicit_env%Bt) 1222 1223! evaluate QS = - B*(\Delta^-1)*B^t 1224 IF (data_size .NE. 0) THEN 1225 CALL DGEMM('N', 'N', n_tiles_tot, n_tiles_tot, data_size, & 1226 -1.0_dp, ps_implicit_env%B, n_tiles_tot, QAinvxBt, & 1227 data_size, 0.0_dp, ps_implicit_env%QS, n_tiles_tot) 1228 END IF 1229 CALL mp_sum(ps_implicit_env%QS, pw_grid_orig%para%group) 1230 1231! evaluate B*1 1232 Bxunit_vec(:) = SUM(ps_implicit_env%B, 2)/ngpts 1233 CALL mp_sum(Bxunit_vec, pw_grid_orig%para%group) 1234! set up R = [QS B*1; (B*1)^t 0] 1235 R = 0.0_dp 1236 R(1:n_tiles_tot, 1:n_tiles_tot) = ps_implicit_env%QS 1237 R(1:n_tiles_tot, n_tiles_tot + 1) = Bxunit_vec 1238 R(n_tiles_tot + 1, 1:n_tiles_tot) = Bxunit_vec 1239! evaluate R^(-1) 1240 ps_implicit_env%Rinv = R 1241 CALL DGETRF(n_tiles_tot + 1, n_tiles_tot + 1, ps_implicit_env%Rinv, n_tiles_tot + 1, ipiv, info) 1242 IF (info .NE. 0) & 1243 CALL cp_abort(__LOCATION__, & 1244 "R is (nearly) singular! Either two Dirichlet constraints are identical or "// & 1245 "you need to reduce the number of tiles.") 1246 CALL DGETRI(n_tiles_tot + 1, ps_implicit_env%Rinv, n_tiles_tot + 1, ipiv, work_arr, n_tiles_tot + 1, info) 1247 IF (info .NE. 0) & 1248 CPABORT("Inversion of R failed!") 1249 1250 DEALLOCATE (QAinvxBt, Bxunit_vec, R, work_arr, ipiv) 1251 CALL pw_pool_release(pw_pool_xpndd) 1252 1253 done_preparing = .TRUE. 1254 CALL mp_sum(done_preparing, pw_grid_orig%para%group) 1255 IF ((unit_nr .GT. 0) .AND. done_preparing) THEN 1256 WRITE (unit_nr, "(T3,A,/,T3,A,/,A)") "POISSON| ... Done. ", REPEAT('-', 78) 1257 END IF 1258 1259 CASE (MIXED_PERIODIC_BC) 1260 1261 ngpts_local = pw_grid_orig%ngpts_local 1262 ngpts = pw_grid_orig%ngpts 1263 npts_local = pw_grid_orig%npts_local 1264 npts = pw_grid_orig%npts 1265 bounds_local = pw_grid_orig%bounds_local 1266 bounds = pw_grid_orig%bounds 1267 dct_env => ps_implicit_env%dct_env 1268 1269! evaluate indices for converting 3D arrays into 1D arrays 1270 lb1 = bounds_local(1, 1); ub1 = bounds_local(2, 1) 1271 lb2 = bounds_local(1, 2); ub2 = bounds_local(2, 2) 1272 lb3 = bounds_local(1, 3); ub3 = bounds_local(2, 3) 1273 1274 IF (pw_grid_orig%para%blocked) THEN 1275 data_size = PRODUCT(npts_local) 1276 ELSE IF (pw_grid_orig%para%ray_distribution) THEN 1277 data_size = ngpts_local 1278 ELSE ! parallel run with np = 1 1279 data_size = PRODUCT(npts_local) 1280 END IF 1281 1282 ALLOCATE (ps_implicit_env%idx_1dto3d(data_size)) 1283 l = 1 1284 DO k = lb3, ub3 1285 DO j = lb2, ub2 1286 DO i = lb1, ub1 1287 ps_implicit_env%idx_1dto3d(l) = (i - lb1 + 1) + & 1288 (j - lb2)*npts_local(1) + & 1289 (k - lb3)*npts_local(1)*npts_local(2) 1290 l = l + 1 1291 END DO 1292 END DO 1293 END DO 1294 1295 n_contacts = SIZE(ps_implicit_env%contacts) 1296 n_tiles_tot = 0 1297 DO j = 1, n_contacts 1298 n_tiles_tot = n_tiles_tot + ps_implicit_env%contacts(j)%dirichlet_bc%n_tiles 1299 END DO 1300 1301 ALLOCATE (ps_implicit_env%B(n_tiles_tot, data_size)) 1302 ALLOCATE (ps_implicit_env%Bt(data_size, n_tiles_tot)) 1303 ALLOCATE (ps_implicit_env%QS(n_tiles_tot, n_tiles_tot)) 1304 ALLOCATE (ps_implicit_env%Rinv(n_tiles_tot + 1, n_tiles_tot + 1)) 1305 ALLOCATE (ps_implicit_env%v_D(n_tiles_tot)) 1306 ALLOCATE (ps_implicit_env%osc_frac(n_tiles_tot)) 1307 ALLOCATE (ps_implicit_env%frequency(n_tiles_tot)) 1308 ALLOCATE (ps_implicit_env%phase(n_tiles_tot)) 1309 1310 ALLOCATE (QAinvxBt(data_size, n_tiles_tot)) 1311 ALLOCATE (Bxunit_vec(n_tiles_tot)) 1312 ALLOCATE (test_vec(n_tiles_tot)) 1313 ALLOCATE (R(n_tiles_tot + 1, n_tiles_tot + 1)) 1314 ALLOCATE (work_arr(n_tiles_tot + 1), ipiv(n_tiles_tot + 1)) 1315 1316! set up B, B^t, (\Delta^-1)*B^t 1317 indx1 = 1 1318 DO j = 1, n_contacts 1319 n_tiles = ps_implicit_env%contacts(j)%dirichlet_bc%n_tiles 1320 indx2 = indx1 + n_tiles - 1 1321 DO i = 1, n_tiles 1322 tile_pw => ps_implicit_env%contacts(j)%dirichlet_bc%tiles(i)%tile%tile_pw 1323 1324 CALL pw_pool_create_pw(pw_pool_orig, pw_in, use_data=REALDATA3D, in_space=REALSPACE) 1325 CALL pw_copy(tile_pw, pw_in) 1326 1327 tile_volume = ps_implicit_env%contacts(j)%dirichlet_bc%tiles(i)%tile%volume 1328 CALL pw_scale(pw_in, 1.0_dp/tile_volume) ! normalize tile_pw 1329 ps_implicit_env%Bt(ps_implicit_env%idx_1dto3d, indx1 + i - 1) = RESHAPE(pw_in%cr3d, (/data_size/)) 1330 1331 CALL pw_pool_create_pw(pw_pool_orig, pw_out, use_data=REALDATA3D, in_space=REALSPACE) 1332 CALL apply_inv_laplace_operator_fft(pw_pool_orig, green, pw_in, pw_out) 1333 QAinvxBt(ps_implicit_env%idx_1dto3d, indx1 + i - 1) = RESHAPE(pw_out%cr3d, (/data_size/)) 1334 ! the electrostatic potential has opposite sign by internal convention 1335 ps_implicit_env%v_D(indx1 + i - 1) = -1.0_dp*ps_implicit_env%contacts(j)%dirichlet_bc%v_D 1336 ps_implicit_env%osc_frac(indx1 + i - 1) = ps_implicit_env%contacts(j)%dirichlet_bc%osc_frac 1337 ps_implicit_env%frequency(indx1 + i - 1) = ps_implicit_env%contacts(j)%dirichlet_bc%frequency 1338 ps_implicit_env%phase(indx1 + i - 1) = ps_implicit_env%contacts(j)%dirichlet_bc%phase 1339 1340 CALL pw_pool_give_back_pw(pw_pool_orig, pw_in) 1341 CALL pw_pool_give_back_pw(pw_pool_orig, pw_out) 1342 END DO 1343 indx1 = indx2 + 1 1344 END DO 1345 ps_implicit_env%B = TRANSPOSE(ps_implicit_env%Bt) 1346 1347! evaluate QS = - B*(\Delta^-1)*B^t 1348 IF (data_size .NE. 0) THEN 1349 CALL DGEMM('N', 'N', n_tiles_tot, n_tiles_tot, data_size, & 1350 -1.0_dp, ps_implicit_env%B, n_tiles_tot, QAinvxBt, & 1351 data_size, 0.0_dp, ps_implicit_env%QS, n_tiles_tot) 1352 END IF 1353 CALL mp_sum(ps_implicit_env%QS, pw_grid_orig%para%group) 1354 1355! evaluate B*1 1356 Bxunit_vec(:) = SUM(ps_implicit_env%B, 2)/ngpts 1357 CALL mp_sum(Bxunit_vec, pw_grid_orig%para%group) 1358! set up R = [QS B*1; (B*1)^t 0] 1359 R = 0.0_dp 1360 R(1:n_tiles_tot, 1:n_tiles_tot) = ps_implicit_env%QS 1361 R(1:n_tiles_tot, n_tiles_tot + 1) = Bxunit_vec 1362 R(n_tiles_tot + 1, 1:n_tiles_tot) = Bxunit_vec 1363! evaluate R^(-1) 1364 ps_implicit_env%Rinv = R 1365 CALL DGETRF(n_tiles_tot + 1, n_tiles_tot + 1, ps_implicit_env%Rinv, n_tiles_tot + 1, ipiv, info) 1366 IF (info .NE. 0) & 1367 CALL cp_abort(__LOCATION__, & 1368 "R is (nearly) singular! Either two Dirichlet constraints are identical or "// & 1369 "you need to reduce the number of tiles.") 1370 CALL DGETRI(n_tiles_tot + 1, ps_implicit_env%Rinv, n_tiles_tot + 1, ipiv, work_arr, n_tiles_tot + 1, info) 1371 IF (info .NE. 0) & 1372 CPABORT("Inversion of R failed!") 1373 1374 DEALLOCATE (QAinvxBt, Bxunit_vec, R, work_arr, ipiv) 1375 1376 done_preparing = .TRUE. 1377 CALL mp_sum(done_preparing, pw_grid_orig%para%group) 1378 IF ((unit_nr .GT. 0) .AND. done_preparing) THEN 1379 WRITE (unit_nr, "(T3,A,/,T3,A,/,A)") "POISSON| ... Done. ", REPEAT('-', 78) 1380 END IF 1381 1382 CASE (PERIODIC_BC, NEUMANN_BC) 1383 1384 ALLOCATE (ps_implicit_env%idx_1dto3d(1)) 1385 ALLOCATE (ps_implicit_env%B(1, 1)) 1386 ALLOCATE (ps_implicit_env%Bt(1, 1)) 1387 ALLOCATE (ps_implicit_env%QS(1, 1)) 1388 ALLOCATE (ps_implicit_env%Rinv(1, 1)) 1389 ALLOCATE (ps_implicit_env%v_D(1)) 1390 ALLOCATE (ps_implicit_env%osc_frac(1)) 1391 ALLOCATE (ps_implicit_env%frequency(1)) 1392 ALLOCATE (ps_implicit_env%phase(1)) 1393 1394 ps_implicit_env%idx_1dto3d = 1 1395 ps_implicit_env%B = 0.0_dp 1396 ps_implicit_env%Bt = 0.0_dp 1397 ps_implicit_env%QS = 0.0_dp 1398 ps_implicit_env%Rinv = 0.0_dp 1399 ps_implicit_env%v_D = 0.0_dp 1400 1401 CASE DEFAULT 1402 CALL cp_abort(__LOCATION__, & 1403 "Please specify the type of boundary conditions using the "// & 1404 "input file keyword BOUNDARY_CONDITIONS.") 1405 END SELECT 1406 1407 CALL timestop(handle) 1408 1409 END SUBROUTINE ps_implicit_prepare_blocks 1410 1411! ************************************************************************************************** 1412!> \brief Evaluates the action of the operator P on a given matrix v, defined 1413!> as: P(v) := - \nabla_r(\ln(\eps)) \cdot \nabla_r(v) 1414!> \param pw_pool pool of pw grid 1415!> \param dielectric dielectric_type contaning eps 1416!> \param v input matrix 1417!> \param Pxv action of the operator P on v 1418!> \par History 1419!> 07.2014 created [Hossein Bani-Hashemian] 1420!> \author Mohammad Hossein Bani-Hashemian 1421! ************************************************************************************************** 1422 SUBROUTINE apply_P_operator(pw_pool, dielectric, v, Pxv) 1423 1424 TYPE(pw_pool_type), POINTER :: pw_pool 1425 TYPE(dielectric_type), INTENT(IN), POINTER :: dielectric 1426 TYPE(pw_type), INTENT(IN), POINTER :: v 1427 TYPE(pw_type), INTENT(INOUT), POINTER :: Pxv 1428 1429 CHARACTER(LEN=*), PARAMETER :: routineN = 'apply_P_operator', & 1430 routineP = moduleN//':'//routineN 1431 1432 INTEGER :: handle, i 1433 TYPE(pw_p_type), DIMENSION(3) :: dln_eps, dv 1434 1435 CALL timeset(routineN, handle) 1436 1437 dln_eps = dielectric%dln_eps 1438 DO i = 1, 3 1439 CALL pw_pool_create_pw(pw_pool, dv(i)%pw, & 1440 use_data=REALDATA3D, in_space=REALSPACE) 1441 END DO 1442 1443 CALL derive_fft(v, dv, pw_pool) 1444 Pxv%cr3d = -(dv(1)%pw%cr3d*dln_eps(1)%pw%cr3d+ & 1445 dv(2)%pw%cr3d*dln_eps(2)%pw%cr3d+ & 1446 dv(3)%pw%cr3d*dln_eps(3)%pw%cr3d) 1447 1448 DO i = 1, 3 1449 CALL pw_pool_give_back_pw(pw_pool, dv(i)%pw) 1450 END DO 1451 1452 CALL timestop(handle) 1453 1454 END SUBROUTINE apply_P_operator 1455 1456! ************************************************************************************************** 1457!> \brief Evaluates the action of the inverse of the Laplace operator on a given 3d matrix 1458!> \param pw_pool pool of pw grid 1459!> \param green green functions for FFT based inverse Laplacian 1460!> \param pw_in pw_in (density) 1461!> \param pw_out pw_out (potential) 1462!> \par History 1463!> 07.2014 created [Hossein Bani-Hashemian] 1464!> \author Mohammad Hossein Bani-Hashemian 1465! ************************************************************************************************** 1466 SUBROUTINE apply_inv_laplace_operator_fft(pw_pool, green, pw_in, pw_out) 1467 1468 TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool 1469 TYPE(greens_fn_type), INTENT(IN), POINTER :: green 1470 TYPE(pw_type), INTENT(IN), POINTER :: pw_in 1471 TYPE(pw_type), INTENT(INOUT), POINTER :: pw_out 1472 1473 CHARACTER(LEN=*), PARAMETER :: routineN = 'apply_inv_laplace_operator_fft', & 1474 routineP = moduleN//':'//routineN 1475 1476 INTEGER :: handle, ig, ng 1477 REAL(dp) :: prefactor 1478 TYPE(pw_grid_type), POINTER :: pw_grid 1479 TYPE(pw_type), POINTER :: pw_in_gs 1480 1481 CALL timeset(routineN, handle) 1482 1483! here I devide by fourpi to cancel out the prefactor fourpi in influence_fn 1484 prefactor = 1.0_dp/fourpi 1485 1486 pw_grid => pw_pool%pw_grid 1487 ng = SIZE(pw_grid%gsq) 1488 1489 CALL pw_pool_create_pw(pw_pool, pw_in_gs, & 1490 use_data=COMPLEXDATA1D, in_space=RECIPROCALSPACE) 1491 1492 CALL pw_transfer(pw_in, pw_in_gs) 1493 DO ig = 1, ng 1494 pw_in_gs%cc(ig) = prefactor*pw_in_gs%cc(ig)*green%influence_fn%cc(ig) 1495 END DO 1496 CALL pw_transfer(pw_in_gs, pw_out) 1497 1498 CALL pw_pool_give_back_pw(pw_pool, pw_in_gs) 1499 1500 CALL timestop(handle) 1501 1502 END SUBROUTINE apply_inv_laplace_operator_fft 1503 1504! ************************************************************************************************** 1505!> \brief Evaluates the action of the inverse of the Laplace operator on a given 1506!> 3d matrix using DCT-I 1507!> \param pw_pool pool of pw grid 1508!> \param green the greens_fn_type data holding a valid dct_influence_fn 1509!> \param pw_in pw_in (density) 1510!> \param pw_out pw_out (potential) 1511!> \par History 1512!> 07.2014 created [Hossein Bani-Hashemian] 1513!> 11.2015 revised [Hossein Bani-Hashemian] 1514!> \author Mohammad Hossein Bani-Hashemian 1515! ************************************************************************************************** 1516 SUBROUTINE apply_inv_laplace_operator_dct(pw_pool, green, pw_in, pw_out) 1517 1518 TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool 1519 TYPE(greens_fn_type), INTENT(IN), POINTER :: green 1520 TYPE(pw_type), INTENT(IN), POINTER :: pw_in 1521 TYPE(pw_type), INTENT(INOUT), POINTER :: pw_out 1522 1523 CHARACTER(LEN=*), PARAMETER :: routineN = 'apply_inv_laplace_operator_dct', & 1524 routineP = moduleN//':'//routineN 1525 1526 INTEGER :: handle, ig, ng 1527 REAL(dp) :: prefactor 1528 TYPE(pw_grid_type), POINTER :: pw_grid 1529 TYPE(pw_type), POINTER :: pw_in_gs 1530 1531 CALL timeset(routineN, handle) 1532 1533! here I devide by fourpi to cancel out the prefactor fourpi in influence_fn 1534 prefactor = 1.0_dp/fourpi 1535 1536 pw_grid => pw_pool%pw_grid 1537 ng = SIZE(pw_grid%gsq) 1538 1539 CALL pw_pool_create_pw(pw_pool, pw_in_gs, & 1540 use_data=COMPLEXDATA1D, in_space=RECIPROCALSPACE) 1541 1542 CALL pw_transfer(pw_in, pw_in_gs) 1543 DO ig = 1, ng 1544 pw_in_gs%cc(ig) = prefactor*pw_in_gs%cc(ig)*green%dct_influence_fn%cc(ig) 1545 END DO 1546 CALL pw_transfer(pw_in_gs, pw_out) 1547 1548 CALL pw_pool_give_back_pw(pw_pool, pw_in_gs) 1549 1550 CALL timestop(handle) 1551 1552 END SUBROUTINE apply_inv_laplace_operator_dct 1553 1554! ************************************************************************************************** 1555!> \brief Evaluates the action of the Laplace operator on a given 3d matrix 1556!> \param pw_pool pool of pw grid 1557!> \param green green functions for FFT based inverse Laplacian 1558!> \param pw_in pw_in (potential) 1559!> \param pw_out pw_out (density) 1560!> \par History 1561!> 07.2014 created [Hossein Bani-Hashemian] 1562!> \author Mohammad Hossein Bani-Hashemian 1563! ************************************************************************************************** 1564 SUBROUTINE apply_laplace_operator_fft(pw_pool, green, pw_in, pw_out) 1565 1566 TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool 1567 TYPE(greens_fn_type), INTENT(IN), POINTER :: green 1568 TYPE(pw_type), INTENT(IN), POINTER :: pw_in 1569 TYPE(pw_type), INTENT(INOUT), POINTER :: pw_out 1570 1571 CHARACTER(LEN=*), PARAMETER :: routineN = 'apply_laplace_operator_fft', & 1572 routineP = moduleN//':'//routineN 1573 1574 INTEGER :: g0_index, handle, ig, ng 1575 LOGICAL :: have_g0 1576 REAL(dp) :: prefactor 1577 TYPE(pw_grid_type), POINTER :: pw_grid 1578 TYPE(pw_type), POINTER :: pw_in_gs 1579 1580 CALL timeset(routineN, handle) 1581 1582! here I multiply by fourpi to cancel out the prefactor fourpi in influence_fn 1583 prefactor = fourpi 1584 1585 pw_grid => pw_pool%pw_grid 1586 ng = SIZE(pw_in%pw_grid%gsq) 1587 have_g0 = green%influence_fn%pw_grid%have_g0 1588 1589 CALL pw_pool_create_pw(pw_pool, pw_in_gs, & 1590 use_data=COMPLEXDATA1D, in_space=RECIPROCALSPACE) 1591 1592 CALL pw_transfer(pw_in, pw_in_gs) 1593 1594 IF (have_g0) THEN 1595 g0_index = green%influence_fn%pw_grid%first_gne0 - 1 1596 pw_in_gs%cc(g0_index) = 0.0_dp 1597 END IF 1598 DO ig = green%influence_fn%pw_grid%first_gne0, ng 1599 pw_in_gs%cc(ig) = prefactor*(pw_in_gs%cc(ig)/green%influence_fn%cc(ig)) 1600 END DO 1601 1602 CALL pw_transfer(pw_in_gs, pw_out) 1603 1604 CALL pw_pool_give_back_pw(pw_pool, pw_in_gs) 1605 1606 CALL timestop(handle) 1607 1608 END SUBROUTINE apply_laplace_operator_fft 1609 1610! ************************************************************************************************** 1611!> \brief Evaluates the action of the Laplace operator on a given 3d matrix using DCT-I 1612!> \param pw_pool pool of pw grid 1613!> \param green the greens_fn_type data holding a valid dct_influence_fn 1614!> \param pw_in pw_in (potential) 1615!> \param pw_out pw_out (density) 1616!> \par History 1617!> 07.2014 created [Hossein Bani-Hashemian] 1618!> 11.2015 revised [Hossein Bani-Hashemian] 1619!> \author Mohammad Hossein Bani-Hashemian 1620! ************************************************************************************************** 1621 SUBROUTINE apply_laplace_operator_dct(pw_pool, green, pw_in, pw_out) 1622 1623 TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool 1624 TYPE(greens_fn_type), INTENT(IN), POINTER :: green 1625 TYPE(pw_type), INTENT(IN), POINTER :: pw_in 1626 TYPE(pw_type), INTENT(INOUT), POINTER :: pw_out 1627 1628 CHARACTER(LEN=*), PARAMETER :: routineN = 'apply_laplace_operator_dct', & 1629 routineP = moduleN//':'//routineN 1630 1631 INTEGER :: g0_index, handle, ig, ng 1632 LOGICAL :: have_g0 1633 REAL(dp) :: prefactor 1634 TYPE(pw_grid_type), POINTER :: pw_grid 1635 TYPE(pw_type), POINTER :: pw_in_gs 1636 1637 CALL timeset(routineN, handle) 1638 1639! here I multiply by fourpi to cancel out the prefactor fourpi in influence_fn 1640 prefactor = fourpi 1641 1642 pw_grid => pw_pool%pw_grid 1643 ng = SIZE(pw_in%pw_grid%gsq) 1644 have_g0 = green%dct_influence_fn%pw_grid%have_g0 1645 1646 CALL pw_pool_create_pw(pw_pool, pw_in_gs, & 1647 use_data=COMPLEXDATA1D, in_space=RECIPROCALSPACE) 1648 1649 CALL pw_transfer(pw_in, pw_in_gs) 1650 1651 IF (have_g0) THEN 1652 g0_index = green%dct_influence_fn%pw_grid%first_gne0 - 1 1653 pw_in_gs%cc(g0_index) = 0.0_dp 1654 END IF 1655 DO ig = green%dct_influence_fn%pw_grid%first_gne0, ng 1656 pw_in_gs%cc(ig) = prefactor*(pw_in_gs%cc(ig)/green%dct_influence_fn%cc(ig)) 1657 END DO 1658 1659 CALL pw_transfer(pw_in_gs, pw_out) 1660 1661 CALL pw_pool_give_back_pw(pw_pool, pw_in_gs) 1662 1663 CALL timestop(handle) 1664 1665 END SUBROUTINE apply_laplace_operator_dct 1666 1667! ************************************************************************************************** 1668!> \brief Evaluates the action of the generalized Poisson operator on a given 3d matrix. 1669!> \param pw_pool pool of pw grid 1670!> \param green green functions for FFT based inverse Laplacian 1671!> \param dielectric dielectric environment 1672!> \param v potential 1673!> \param density density 1674!> \par History 1675!> 07.2014 created [Hossein Bani-Hashemian] 1676!> \author Mohammad Hossein Bani-Hashemian 1677! ************************************************************************************************** 1678 SUBROUTINE apply_poisson_operator_fft(pw_pool, green, dielectric, v, density) 1679 1680 TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool 1681 TYPE(greens_fn_type), INTENT(IN), POINTER :: green 1682 TYPE(dielectric_type), INTENT(IN), POINTER :: dielectric 1683 TYPE(pw_type), INTENT(IN), POINTER :: v 1684 TYPE(pw_type), INTENT(INOUT), POINTER :: density 1685 1686 CHARACTER(LEN=*), PARAMETER :: routineN = 'apply_poisson_operator_fft', & 1687 routineP = moduleN//':'//routineN 1688 1689 INTEGER :: handle 1690 TYPE(pw_type), POINTER :: Pxv 1691 1692 CALL timeset(routineN, handle) 1693 1694 CALL pw_pool_create_pw(pw_pool, Pxv, & 1695 use_data=REALDATA3D, in_space=REALSPACE) 1696 1697 CALL apply_P_operator(pw_pool, dielectric, v, Pxv) 1698 CALL apply_laplace_operator_fft(pw_pool, green, v, density) 1699 CALL pw_axpy(Pxv, density) 1700 1701 CALL pw_pool_give_back_pw(pw_pool, Pxv) 1702 1703 CALL timestop(handle) 1704 1705 END SUBROUTINE apply_poisson_operator_fft 1706 1707! ************************************************************************************************** 1708!> \brief Evaluates the action of the generalized Poisson operator on a given 1709!> 3d matrix using DCT-I. 1710!> \param pw_pool pool of pw grid 1711!> \param green the greens_fn_type data holding a valid dct_influence_fn 1712!> \param dielectric dielectric environment 1713!> \param v potential 1714!> \param density density 1715!> \par History 1716!> 07.2014 created [Hossein Bani-Hashemian] 1717!> 11.2015 revised [Hossein Bani-Hashemian] 1718!> \author Mohammad Hossein Bani-Hashemian 1719! ************************************************************************************************** 1720 SUBROUTINE apply_poisson_operator_dct(pw_pool, green, dielectric, v, density) 1721 1722 TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool 1723 TYPE(greens_fn_type), INTENT(IN), POINTER :: green 1724 TYPE(dielectric_type), INTENT(IN), POINTER :: dielectric 1725 TYPE(pw_type), INTENT(IN), POINTER :: v 1726 TYPE(pw_type), INTENT(INOUT), POINTER :: density 1727 1728 CHARACTER(LEN=*), PARAMETER :: routineN = 'apply_poisson_operator_dct', & 1729 routineP = moduleN//':'//routineN 1730 1731 INTEGER :: handle 1732 TYPE(pw_type), POINTER :: Pxv 1733 1734 CALL timeset(routineN, handle) 1735 1736 CALL pw_pool_create_pw(pw_pool, Pxv, & 1737 use_data=REALDATA3D, in_space=REALSPACE) 1738 1739 CALL apply_P_operator(pw_pool, dielectric, v, Pxv) 1740 CALL apply_laplace_operator_dct(pw_pool, green, v, density) 1741 CALL pw_axpy(Pxv, density) 1742 1743 CALL pw_pool_give_back_pw(pw_pool, Pxv) 1744 1745 CALL timestop(handle) 1746 1747 END SUBROUTINE apply_poisson_operator_dct 1748 1749! ************************************************************************************************** 1750!> \brief Computes the extra contribution (v_eps) 1751!> v_eps = - \frac{1}{8*\pi} * |\nabla_r(v)|^2 * \frac{d \eps}{d \rho} 1752!> to the functional derivative of the Hartree energy wrt the density, being 1753!> attributed to the dependecy of the dielectric constant to the charge density. 1754!> [see V. M. Sanchez, M. Sued, and D. A. Scherlis, J. Chem. Phys. 131, 174108 (2009)] 1755!> \param pw_pool pool of the original plane-wave grid 1756!> \param dielectric dielectric environment 1757!> \param v Hartree potential 1758!> \param v_eps v_eps 1759!> \par History 1760!> 08.2014 created [Hossein Bani-Hashemian] 1761!> \author Mohammad Hossein Bani-Hashemian 1762! ************************************************************************************************** 1763 SUBROUTINE ps_implicit_compute_veps(pw_pool, dielectric, v, v_eps) 1764 1765 TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool 1766 TYPE(dielectric_type), INTENT(IN), POINTER :: dielectric 1767 TYPE(pw_type), INTENT(IN), POINTER :: v 1768 TYPE(pw_type), INTENT(INOUT), POINTER :: v_eps 1769 1770 CHARACTER(LEN=*), PARAMETER :: routineN = 'ps_implicit_compute_veps', & 1771 routineP = moduleN//':'//routineN 1772 1773 INTEGER :: handle, i 1774 REAL(dp) :: eightpi 1775 TYPE(pw_p_type), DIMENSION(3) :: dv 1776 TYPE(pw_type), POINTER :: deps_drho, dv2 1777 1778 CALL timeset(routineN, handle) 1779 1780 eightpi = 2*fourpi 1781 deps_drho => dielectric%deps_drho 1782 1783 CALL pw_pool_create_pw(pw_pool, dv2, & 1784 use_data=REALDATA3D, in_space=REALSPACE) 1785 DO i = 1, 3 1786 CALL pw_pool_create_pw(pw_pool, dv(i)%pw, & 1787 use_data=REALDATA3D, in_space=REALSPACE) 1788 END DO 1789 1790 CALL derive_fft(v, dv, pw_pool) 1791 1792! evaluate |\nabla_r(v)|^2 1793 dv2%cr3d = dv(1)%pw%cr3d**2 + dv(2)%pw%cr3d**2 + dv(3)%pw%cr3d**2 1794 1795 v_eps%cr3d = -(1.0_dp/eightpi)*(dv2%cr3d*deps_drho%cr3d) 1796 1797 CALL pw_pool_give_back_pw(pw_pool, dv2) 1798 DO i = 1, 3 1799 CALL pw_pool_give_back_pw(pw_pool, dv(i)%pw) 1800 END DO 1801 1802 CALL timestop(handle) 1803 1804 END SUBROUTINE ps_implicit_compute_veps 1805 1806! ************************************************************************************************** 1807!> \brief Computes the Hartree energy 1808!> \param density electronic density 1809!> \param v Hartree potential 1810!> \param ehartree Hartree energy 1811!> \par History 1812!> 06.2015 created [Hossein Bani-Hashemian] 1813!> \author Mohammad Hossein Bani-Hashemian 1814! ************************************************************************************************** 1815 SUBROUTINE compute_ehartree_periodic_bc(density, v, ehartree) 1816 1817 TYPE(pw_type), INTENT(IN), POINTER :: density, v 1818 REAL(dp), INTENT(OUT) :: ehartree 1819 1820 CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_ehartree_periodic_bc', & 1821 routineP = moduleN//':'//routineN 1822 1823 INTEGER :: handle 1824 1825 CALL timeset(routineN, handle) 1826 1827! E_H = \frac{1}{2} * \int \rho * v dr 1828 ehartree = 0.5_dp*pw_integral_ab(density, v) 1829 1830 CALL timestop(handle) 1831 1832 END SUBROUTINE compute_ehartree_periodic_bc 1833 1834! ************************************************************************************************** 1835!> \brief Computes the Hartree energy 1836!> \param dielectric dielectric environment 1837!> \param density electronic density 1838!> \param Btxlambda B^t * \lambda (\lambda is the vector of Lagrange multipliers 1839!> and B^t is the transpose of the boundary operator 1840!> \param v Hartree potential 1841!> \param ehartree Hartree energy 1842!> \param electric_enthalpy electric enthalpy 1843!> \par History 1844!> 06.2015 created [Hossein Bani-Hashemian] 1845!> \author Mohammad Hossein Bani-Hashemian 1846! ************************************************************************************************** 1847 SUBROUTINE compute_ehartree_mixed_bc(dielectric, density, Btxlambda, v, ehartree, electric_enthalpy) 1848 1849 TYPE(dielectric_type), INTENT(IN), POINTER :: dielectric 1850 TYPE(pw_type), INTENT(IN), POINTER :: density 1851 REAL(dp), ALLOCATABLE, DIMENSION(:, :, :), & 1852 INTENT(IN) :: Btxlambda 1853 TYPE(pw_type), INTENT(IN), POINTER :: v 1854 REAL(dp), INTENT(OUT) :: ehartree, electric_enthalpy 1855 1856 CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_ehartree_mixed_bc', & 1857 routineP = moduleN//':'//routineN 1858 1859 INTEGER :: handle 1860 REAL(dp) :: dvol, ehartree_rho, ehartree_rho_cstr 1861 TYPE(pw_grid_type), POINTER :: pw_grid 1862 TYPE(pw_type), POINTER :: eps 1863 1864 CALL timeset(routineN, handle) 1865 1866 pw_grid => v%pw_grid 1867 eps => dielectric%eps 1868 1869 dvol = pw_grid%dvol 1870 1871! E_H = \frac{1}{2} * \int \rho * v dr + \frac{1}{8 \pi} * \int Btxlambda * v dr 1872! the sign of the second term depends on the sign chosen for the Lagrange multiplier 1873! term in the variational form 1874 ehartree_rho = accurate_sum(density%cr3d*v%cr3d) 1875 ehartree_rho_cstr = accurate_sum(eps%cr3d*Btxlambda*v%cr3d/fourpi) 1876 ehartree_rho = 0.5_dp*ehartree_rho*dvol 1877 ehartree_rho_cstr = 0.5_dp*ehartree_rho_cstr*dvol 1878 CALL mp_sum(ehartree_rho, pw_grid%para%group) 1879 CALL mp_sum(ehartree_rho_cstr, pw_grid%para%group) 1880 electric_enthalpy = ehartree_rho + ehartree_rho_cstr 1881 ehartree = ehartree_rho - ehartree_rho_cstr 1882 1883 CALL timestop(handle) 1884 1885 END SUBROUTINE compute_ehartree_mixed_bc 1886 1887! ************************************************************************************************** 1888!> \brief Computes the (normalized) preconditioned residual norm error and the 1889!> normalized absolute error 1890!> \param pw_pool pool of the original plane-wave grid 1891!> \param green greens functions for FFT based inverse Laplacian 1892!> \param res_new residual 1893!> \param v_old old v 1894!> \param v_new new v 1895!> \param QAinvxres_new Delta^-1(res_new) 1896!> \param pres_error preconditioned residual norm error 1897!> \param nabs_error normalized absolute error 1898!> \par History 1899!> 07.2014 created [Hossein Bani-Hashemian] 1900!> \author Mohammad Hossein Bani-Hashemian 1901! ************************************************************************************************** 1902 SUBROUTINE ps_implicit_compute_error_fft(pw_pool, green, res_new, v_old, v_new, & 1903 QAinvxres_new, pres_error, nabs_error) 1904 1905 TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool 1906 TYPE(greens_fn_type), INTENT(IN), POINTER :: green 1907 TYPE(pw_type), INTENT(IN), POINTER :: res_new, v_old, v_new 1908 TYPE(pw_type), INTENT(INOUT), POINTER :: QAinvxres_new 1909 REAL(dp), INTENT(OUT) :: pres_error, nabs_error 1910 1911 CHARACTER(LEN=*), PARAMETER :: routineN = 'ps_implicit_compute_error_fft', & 1912 routineP = moduleN//':'//routineN 1913 1914 INTEGER :: handle 1915 REAL(dp) :: vol 1916 1917 CALL timeset(routineN, handle) 1918 1919 vol = pw_pool%pw_grid%vol 1920 1921! evaluate \Delta^-1(res) = \Delta^-1 (g - \Delta(v_new) - P(v_new) + Bt \lambda) 1922 CALL apply_inv_laplace_operator_fft(pw_pool, green, res_new, QAinvxres_new) 1923! (normalized) preconditioned residual norm error : 1924 pres_error = accurate_sum(QAinvxres_new%cr3d(:, :, :)**2) 1925 CALL mp_sum(pres_error, pw_pool%pw_grid%para%group) 1926 pres_error = SQRT(pres_error)/vol 1927 1928! normalized absolute error : 1929! nabs_error := \frac{\| v_old - v_new \|}{volume} 1930 nabs_error = accurate_sum(ABS(v_old%cr3d-v_new%cr3d)**2) 1931 CALL mp_sum(nabs_error, pw_pool%pw_grid%para%group) 1932 nabs_error = SQRT(nabs_error)/vol 1933 1934 CALL timestop(handle) 1935 1936 END SUBROUTINE ps_implicit_compute_error_fft 1937 1938! ************************************************************************************************** 1939!> \brief Computes the (normalized) preconditioned residual norm error and the 1940!> normalized absolute error 1941!> \param pw_pool pool of the original plane-wave grid 1942!> \param green the greens_fn_type data holding a valid dct_influence_fn 1943!> \param res_new residual 1944!> \param v_old old v 1945!> \param v_new new v 1946!> \param QAinvxres_new Delta^-1(res_new) 1947!> \param pres_error preconditioned residual norm error 1948!> \param nabs_error normalized absolute error 1949!> \par History 1950!> 07.2014 created [Hossein Bani-Hashemian] 1951!> 11.2015 revised [Hossein Bani-Hashemian] 1952!> \author Mohammad Hossein Bani-Hashemian 1953! ************************************************************************************************** 1954 SUBROUTINE ps_implicit_compute_error_dct(pw_pool, green, res_new, v_old, v_new, & 1955 QAinvxres_new, pres_error, nabs_error) 1956 1957 TYPE(pw_pool_type), INTENT(IN), POINTER :: pw_pool 1958 TYPE(greens_fn_type), INTENT(IN), POINTER :: green 1959 TYPE(pw_type), INTENT(IN), POINTER :: res_new, v_old, v_new 1960 TYPE(pw_type), INTENT(INOUT), POINTER :: QAinvxres_new 1961 REAL(dp), INTENT(OUT) :: pres_error, nabs_error 1962 1963 CHARACTER(LEN=*), PARAMETER :: routineN = 'ps_implicit_compute_error_dct', & 1964 routineP = moduleN//':'//routineN 1965 1966 INTEGER :: handle 1967 REAL(dp) :: vol 1968 1969 CALL timeset(routineN, handle) 1970 1971 vol = pw_pool%pw_grid%vol 1972 1973! evaluate \Delta^-1(res) = \Delta^-1 (g - \Delta(v_new) - P(v_new) + Bt \lambda) 1974 CALL apply_inv_laplace_operator_dct(pw_pool, green, res_new, QAinvxres_new) 1975! (normalized) preconditioned residual norm error : 1976 pres_error = accurate_sum(QAinvxres_new%cr3d(:, :, :)**2) 1977 CALL mp_sum(pres_error, pw_pool%pw_grid%para%group) 1978 pres_error = SQRT(pres_error)/vol 1979 1980! normalized absolute error : 1981! nabs_error := \frac{\| v_old - v_new \|}{volume} 1982 nabs_error = accurate_sum(ABS(v_old%cr3d-v_new%cr3d)**2) 1983 CALL mp_sum(nabs_error, pw_pool%pw_grid%para%group) 1984 nabs_error = SQRT(nabs_error)/vol 1985 1986 CALL timestop(handle) 1987 1988 END SUBROUTINE ps_implicit_compute_error_dct 1989 1990! ************************************************************************************************** 1991!> \brief output of the implicit (iterative) Poisson solver 1992!> \param iter current iteration 1993!> \param pres_error preconditioned residual norm error 1994!> \param nabs_error normalized absolute error 1995!> \param outp_unit output unit 1996!> \par History 1997!> 07.2014 created [Hossein Bani-Hashemian] 1998!> \author Mohammad Hossein Bani-Hashemian 1999! ************************************************************************************************** 2000 SUBROUTINE ps_implicit_output(iter, pres_error, nabs_error, outp_unit) 2001 2002 INTEGER, INTENT(IN) :: iter 2003 REAL(dp), INTENT(IN) :: pres_error, nabs_error 2004 INTEGER, INTENT(OUT) :: outp_unit 2005 2006 CHARACTER(LEN=*), PARAMETER :: routineN = 'ps_implicit_output', & 2007 routineP = moduleN//':'//routineN 2008 INTEGER, PARAMETER :: low_print_level = 1 2009 2010 INTEGER :: handle 2011 TYPE(cp_logger_type), POINTER :: logger 2012 2013 CALL timeset(routineN, handle) 2014 2015 logger => cp_get_default_logger() 2016 IF (logger%para_env%mepos .EQ. logger%para_env%source) THEN 2017 outp_unit = cp_logger_get_default_unit_nr(logger, local=.TRUE.) 2018 ELSE 2019 outp_unit = -1 2020 ENDIF 2021 2022 IF (logger%iter_info%print_level .GT. low_print_level) THEN 2023 IF ((outp_unit .GT. 0) .AND. (iter .EQ. 1)) THEN 2024 WRITE (outp_unit, '(T3,A)') & 2025 "POISSON| iter pres error nabs error E_hartree delta E" 2026 END IF 2027 2028 IF (outp_unit .GT. 0) THEN 2029 WRITE (outp_unit, '(T3,A,I6,5X,E13.4,3X,E13.4)', ADVANCE='NO') & 2030 "POISSON| ", iter, pres_error, nabs_error 2031 END IF 2032 END IF 2033 2034 CALL timestop(handle) 2035 2036 END SUBROUTINE ps_implicit_output 2037 2038! ************************************************************************************************** 2039!> \brief reports the Hartree energy in every iteration 2040!> \param ps_implicit_env the implicit poisson solver environment 2041!> \param outp_unit output unit 2042!> \param ehartree Hartree energy 2043!> \par History 2044!> 07.2014 created [Hossein Bani-Hashemian] 2045!> \author Mohammad Hossein Bani-Hashemian 2046! ************************************************************************************************** 2047 SUBROUTINE ps_implicit_report_ehartree(ps_implicit_env, outp_unit, ehartree) 2048 2049 TYPE(ps_implicit_type) :: ps_implicit_env 2050 INTEGER, INTENT(IN) :: outp_unit 2051 REAL(dp), INTENT(IN) :: ehartree 2052 2053 CHARACTER(LEN=*), PARAMETER :: routineN = 'ps_implicit_report_ehartree', & 2054 routineP = moduleN//':'//routineN 2055 INTEGER, PARAMETER :: low_print_level = 1 2056 2057 INTEGER :: handle 2058 TYPE(cp_logger_type), POINTER :: logger 2059 2060 logger => cp_get_default_logger() 2061 CALL timeset(routineN, handle) 2062 IF (logger%iter_info%print_level .GT. low_print_level) THEN 2063 IF (outp_unit .GT. 0) WRITE (outp_unit, '(F19.10,E10.2)') & 2064 ehartree, ehartree - ps_implicit_env%ehartree 2065 END IF 2066 CALL timestop(handle) 2067 2068 END SUBROUTINE ps_implicit_report_ehartree 2069 2070! ************************************************************************************************** 2071!> \brief reports the final number of iteration 2072!> \param iter the iteration number after exiting the main loop 2073!> \param max_iter maximum number of iterations 2074!> \param outp_unit output unit 2075!> \par History 2076!> 02.2016 created [Hossein Bani-Hashemian] 2077!> \author Mohammad Hossein Bani-Hashemian 2078! ************************************************************************************************** 2079 SUBROUTINE ps_implicit_print_convergence_msg(iter, max_iter, outp_unit) 2080 2081 INTEGER, INTENT(IN) :: iter, max_iter, outp_unit 2082 2083 CHARACTER(LEN=*), PARAMETER :: routineN = 'ps_implicit_print_convergence_msg', & 2084 routineP = moduleN//':'//routineN 2085 2086 CHARACTER(LEN=12) :: msg 2087 INTEGER :: handle, last_iter 2088 2089 CALL timeset(routineN, handle) 2090 2091 last_iter = iter - 1 2092 IF (last_iter .EQ. 1) THEN 2093 msg = " iteration. " 2094 ELSE 2095 msg = " iterations." 2096 END IF 2097 2098 IF (outp_unit .GT. 0) THEN 2099 IF (last_iter .EQ. max_iter) THEN 2100 WRITE (outp_unit, '(T3,A)') & 2101 "POISSON| No convergence achieved within the maximum number of iterations." 2102 END IF 2103 IF (last_iter .LT. max_iter) THEN 2104 WRITE (outp_unit, '(T3,A,I0,A)') & 2105 "POISSON| Poisson solver converged in ", last_iter, msg 2106 END IF 2107 END IF 2108 CALL timestop(handle) 2109 2110 END SUBROUTINE ps_implicit_print_convergence_msg 2111 2112! ************************************************************************************************** 2113!> \brief computes the derivative of a function using FFT 2114!> \param f input funcition 2115!> \param df derivative of f 2116!> \param pw_pool pool of plane-wave grid 2117! ************************************************************************************************** 2118 SUBROUTINE derive_fft(f, df, pw_pool) 2119 2120 TYPE(pw_type), POINTER :: f 2121 TYPE(pw_p_type), DIMENSION(3), INTENT(INOUT) :: df 2122 TYPE(pw_pool_type), POINTER :: pw_pool 2123 2124 CHARACTER(LEN=*), PARAMETER :: routineN = 'derive_fft', routineP = moduleN//':'//routineN 2125 2126 INTEGER :: handle, i 2127 INTEGER, DIMENSION(3) :: nd 2128 TYPE(pw_p_type), DIMENSION(2) :: work_gs 2129 2130 CALL timeset(routineN, handle) 2131 2132 DO i = 1, 2 2133 NULLIFY (work_gs(i)%pw) 2134 CALL pw_pool_create_pw(pw_pool, work_gs(i)%pw, & 2135 use_data=COMPLEXDATA1D, in_space=RECIPROCALSPACE) 2136 END DO 2137 2138 CALL pw_transfer(f, work_gs(1)%pw) 2139 DO i = 1, 3 2140 nd(:) = 0 2141 nd(i) = 1 2142 CALL pw_copy(work_gs(1)%pw, work_gs(2)%pw) 2143 CALL pw_derive(work_gs(2)%pw, nd(:)) 2144 CALL pw_transfer(work_gs(2)%pw, df(i)%pw) 2145 END DO 2146 2147 DO i = 1, 2 2148 CALL pw_pool_give_back_pw(pw_pool, work_gs(i)%pw) 2149 END DO 2150 2151 CALL timestop(handle) 2152 2153 END SUBROUTINE derive_fft 2154 2155! ************************************************************************************************** 2156!> \brief converts a 1D array to a 3D array (contiguous layout) 2157!> \param idx_1dto3d mapping of indices 2158!> \param arr1d input 1D array 2159!> \param arr3d input 3D array 2160! ************************************************************************************************** 2161 SUBROUTINE convert_1dto3d(idx_1dto3d, arr1d, arr3d) 2162 2163 INTEGER, DIMENSION(:), INTENT(IN), POINTER :: idx_1dto3d 2164 REAL(dp), ALLOCATABLE, DIMENSION(:), INTENT(IN) :: arr1d 2165 REAL(dp), ALLOCATABLE, DIMENSION(:, :, :), & 2166 INTENT(INOUT) :: arr3d 2167 2168 CHARACTER(LEN=*), PARAMETER :: routineN = 'convert_1dto3d', routineP = moduleN//':'//routineN 2169 2170 INTEGER :: handle, i, j, k, l, lb1, lb2, lb3, & 2171 npts1, npts2, npts3, ub1, ub2, ub3 2172 2173 CALL timeset(routineN, handle) 2174 2175 lb1 = LBOUND(arr3d, 1); ub1 = UBOUND(arr3d, 1) 2176 lb2 = LBOUND(arr3d, 2); ub2 = UBOUND(arr3d, 2) 2177 lb3 = LBOUND(arr3d, 3); ub3 = UBOUND(arr3d, 3) 2178 2179 npts1 = ub1 - lb1 + 1 2180 npts2 = ub2 - lb2 + 1 2181 npts3 = ub3 - lb3 + 1 2182 2183 DO l = 1, SIZE(idx_1dto3d) 2184 k = ((idx_1dto3d(l) - 1)/(npts1*npts2)) + lb3 2185 j = ((idx_1dto3d(l) - 1) - (k - lb3)*npts1*npts2)/npts1 + lb2 2186 i = idx_1dto3d(l) - ((j - lb2)*npts1 + (k - lb3)*npts1*npts2) + lb1 - 1 2187 arr3d(i, j, k) = arr1d(l) 2188 END DO 2189 2190 CALL timestop(handle) 2191 2192 END SUBROUTINE convert_1dto3d 2193 2194! ************************************************************************************************** 2195!> \brief Returns the voltage of a tile. In case an alternating field is used, the oltage is a funtion of time 2196!> \param time ... 2197!> \param v_D ... 2198!> \param osc_frac ... 2199!> \param frequency ... 2200!> \param phase ... 2201!> \param v_D_new ... 2202! ************************************************************************************************** 2203 SUBROUTINE get_voltage(time, v_D, osc_frac, frequency, phase, v_D_new) 2204 2205 REAL(dp), INTENT(IN) :: time 2206 REAL(dp), DIMENSION(:), INTENT(IN) :: v_D, osc_frac, frequency, phase 2207 REAL(dp), ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: v_D_new 2208 2209 CHARACTER(LEN=*), PARAMETER :: routineN = 'get_voltage', routineP = moduleN//':'//routineN 2210 2211 INTEGER :: handle, i 2212 2213 CALL timeset(routineN, handle) 2214 2215 ALLOCATE (v_D_new(SIZE(v_D))) 2216 2217 DO i = 1, SIZE(v_D) 2218 v_D_new(i) = v_D(i)*(1 - osc_frac(i)) + & 2219 v_D(i)*osc_frac(i)*COS(2*pi*time*frequency(i) + phase(i)) 2220 END DO 2221 2222 CALL timestop(handle) 2223 2224 END SUBROUTINE get_voltage 2225 2226END MODULE ps_implicit_methods 2227