1!! Copyright (C) 2013 J. Alberdi-Rodriguez 2!! 3!! This program is free software; you can redistribute it and/or modify 4!! it under the terms of the GNU General Public License as published by 5!! the Free Software Foundation; either version 2, or (at your option) 6!! any later version. 7!! 8!! This program is distributed in the hope that it will be useful, 9!! but WITHOUT ANY WARRANTY; without even the implied warranty of 10!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 11!! GNU General Public License for more details. 12!! 13!! You should have received a copy of the GNU General Public License 14!! along with this program; if not, write to the Free Software 15!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 16!! 02110-1301, USA. 17!! 18 19#include "global.h" 20 21module poisson_psolver_oct_m 22 use cube_function_oct_m 23 use cube_oct_m 24 use fourier_space_oct_m 25 use global_oct_m 26 use kpoints_oct_m 27 use mesh_cube_parallel_map_oct_m 28 use mesh_oct_m 29#ifdef HAVE_PSOLVER 30 use mesh_function_oct_m 31#endif 32 use messages_oct_m 33 use mpi_oct_m 34 use namespace_oct_m 35 use parser_oct_m 36 use profiling_oct_m 37 38 !! Support for PSolver from BigDFT 39 40#ifdef HAVE_LIBISF 41 use poisson_solver 42 use dynamic_memory 43#endif 44 45#ifdef HAVE_PSOLVER 46 use poisson_solver 47 use dictionaries, dict_set => set 48 use yaml_output, only: yaml_map 49 use wrapper_MPI 50#endif 51 52 53 implicit none 54 55 private 56 public :: & 57 poisson_psolver_t, & 58 poisson_psolver_init, & 59 poisson_psolver_end, & 60 poisson_psolver_reinit, & 61 poisson_psolver_global_solve, & 62 poisson_psolver_parallel_solve, & 63 poisson_psolver_get_dims 64 65 type poisson_psolver_t 66 private 67 type(fourier_space_op_t) :: coulb !< object for Fourier space operations 68#if (defined HAVE_LIBISF) || (defined HAVE_PSOLVER) 69 type(coulomb_operator) :: kernel !< choice of kernel, one of options above 70#endif 71 !> Indicates the boundary conditions (BC) of the problem: 72 !! 'F' free BC, isolated systems. 73 !! The program calculates the solution as if the given density is 74 !! "alone" in R^3 space. 75 !! 'S' surface BC, isolated in y direction, periodic in xz plane 76 !! The given density is supposed to be periodic in the xz plane, 77 !! so the dimensions in these direction mus be compatible with the FFT 78 !! Beware of the fact that the isolated direction is y! 79 !! 'P' periodic BC. 80 !! The density is supposed to be periodic in all the three directions, 81 !! then all the dimensions must be compatible with the FFT. 82 !! No need for setting up the kernel. 83 character(len = 1) :: geocode = "F" !< 'F' free boundary contition 84 !> Indicates the distribution of the data of the input/output array: 85 !! 'G' global data. Each process has the whole array of the density 86 !! which will be overwritten with the whole array of the potential 87 !! 'D' distributed data. Each process has only the needed part of the density 88 !! and of the potential. The data distribution is such that each processor 89 !! has the xy planes needed for the calculation AND for the evaluation of the 90 !! gradient, needed for XC part, and for the White-Bird correction, which 91 !! may lead up to 8 planes more on each side. Due to this fact, the information 92 !! between the processors may overlap. 93 character(len = 1), public :: datacode = "G" 94 95 integer :: isf_order !< order of the interpolating scaling functions used in the decomposition 96 integer :: localnscatterarr(5) 97 integer :: rs_n_global(3) !< total size of the fft in each direction in real space 98 integer :: rs_istart(1:3) !< where does the local portion of the function start in real space 99#ifdef HAVE_PSOLVER 100 type(dictionary), pointer :: inputs !input parameters 101#endif 102 double precision :: offset 103 end type poisson_psolver_t 104 105#ifdef HAVE_PSOLVER 106 logical, save :: flib_initialized = .false. 107#endif 108 109contains 110 111 !----------------------------------------------------------------- 112 subroutine poisson_psolver_init(this, namespace, mesh, cube, mu, qq) 113 type(poisson_psolver_t), intent(out) :: this 114 type(namespace_t), intent(in) :: namespace 115 type(mesh_t), intent(inout) :: mesh 116 type(cube_t), intent(inout) :: cube 117 FLOAT, intent(in) :: mu 118 FLOAT, intent(in) :: qq(1:MAX_DIM) 119 120 logical data_is_parallel 121#ifdef HAVE_PSOLVER 122 FLOAT :: alpha, beta, gamma 123 FLOAT :: modq2 124 type(mpi_environment) mpi_env 125#endif 126 127 PUSH_SUB(poisson_psolver_init) 128 129#ifdef HAVE_PSOLVER 130 if(.not.flib_initialized) then 131 call f_lib_initialize() 132 flib_initialized = .true. 133 end if 134 135 call dict_init(this%inputs) 136#elif HAVE_LIBISF 137 call f_lib_initialize() 138#endif 139 140 select case(mesh%sb%periodic_dim) 141 case(0) 142 ! Free BC 143 this%geocode = "F" 144 case(1) 145 ! Wire BC 146 this%geocode = "W" 147 call messages_not_implemented("PSolver support for 1D periodic boundary conditions.") 148 case(2) 149 ! Surface BC 150 this%geocode = "S" 151 call messages_not_implemented("PSolver support for 2D periodic boundary conditions.") 152 case(3) 153 ! Periodic BC 154 this%geocode = "P" 155 call messages_experimental("PSolver support for 3D periodic boundary conditions.") 156 end select 157 158 159#ifdef HAVE_PSOLVER 160 !Verbosity switch 161 call dict_set(this%inputs//'setup'//'verbose', .false.) 162 !Order of the Interpolating Scaling Function family 163 call dict_set(this%inputs//'kernel'//'isf_order', 16) 164 !Mu screening parameter 165 call dict_set(this%inputs//'kernel'//'screening', mu) 166 !Calculation of the stress tensor 167 call dict_set(this%inputs//'kernel'//'stress_tensor', .false.) 168#else 169 this%isf_order = 16 170#endif 171 172 !%Variable PoissonSolverPSolverParallelData 173 !%Type logical 174 !%Section Hamiltonian::Poisson::PSolver 175 !%Default yes 176 !%Description 177 !% Indicates whether data is partitioned within the PSolver library. 178 !% If data is distributed among processes, Octopus uses parallel data-structures 179 !% and, thus, less memory. 180 !% If "yes", data is parallelized. The <i>z</i>-axis of the input vector 181 !% is split among the MPI processes. 182 !% If "no", entire input and output vector is saved in all the MPI processes. 183 !% If k-points parallelization is used, "no" must be selected. 184 !%End 185 call parse_variable(namespace, 'PoissonSolverPSolverParallelData', .true., data_is_parallel) 186 187 call messages_obsolete_variable(namespace, 'PoissonSolverISFParallelData', 'PoissonSolverPSolverParallelData') 188 189 if (data_is_parallel) then 190#ifdef HAVE_PSOLVER 191 call dict_set(this%inputs//'setup'//'global_data', .false.) 192#endif 193 this%datacode = "D" 194 else 195#ifdef HAVE_PSOLVER 196 call dict_set(this%inputs//'setup'//'global_data', .true.) 197#endif 198 this%datacode = "G" 199 end if 200 201#ifdef HAVE_PSOLVER 202 call dict_set(this%inputs//'setup'//'verbose', debug%info) 203 204 alpha = mesh%sb%alpha*M_PI/(CNST(180.0)) 205 beta = mesh%sb%beta*M_PI/(CNST(180.0)) 206 gamma = mesh%sb%gamma*M_PI/(CNST(180.0)) 207 208 ! Previously, pkernel_init set the communicator used within PSolver to comm_world. 209 ! This can be overwritten by passing an optional argument of type(mpi_environment) 210 ! to pkernel_init(). This data type is defined within the wrapper_MPI module of 211 ! the Futile library. Futile is a prerequisit for PSolver. 212 213 ! TODO: check that cube%mpi_grp corresponds to the correct mpi group, when parallelizing, e.g., over systems !! 214 215 call mpi_environment_set(mpi_env, cube%mpi_grp%rank, cube%mpi_grp%size, cube%mpi_grp%comm, cube%mpi_grp%size ) 216 217 this%kernel = pkernel_init(cube%mpi_grp%rank, cube%mpi_grp%size, this%inputs, this%geocode, cube%rs_n_global, & 218 mesh%spacing, alpha_bc = alpha, beta_ac = beta, gamma_ab = gamma, mpi_env = mpi_env) 219 call pkernel_set(this%kernel, verbose=debug%info) 220 221 !G=0 component 222 modq2 = sum(qq(1:mesh%sb%periodic_dim)**2) 223 if (modq2 > M_EPSILON) then 224 this%offset = M_ONE/modq2 225 else 226 this%offset = M_ZERO 227 end if 228 229 !Screened coulomb potential (erfc function) 230 if (mu > M_EPSILON) then 231 if(modq2 > M_EPSILON) then 232 this%offset = this%offset*(M_ONE - exp(-modq2/((M_TWO*mu)**2))) 233 else 234 !Analytical limit of 1/|q|^2*(1-exp(-|q|^2/4mu^2)) 235 this%offset = M_ONE/((M_TWO*mu)**2) 236 end if 237 end if 238 this%offset = this%offset*M_FOUR*M_PI 239 240#elif HAVE_LIBISF 241 this%kernel = pkernel_init(.false., cube%mpi_grp%rank, cube%mpi_grp%size, 0, this%geocode, cube%rs_n_global, & 242 mesh%spacing, this%isf_order) 243 call pkernel_set(this%kernel, .false.) 244#endif 245 246 POP_SUB(poisson_psolver_init) 247 end subroutine poisson_psolver_init 248 249 !----------------------------------------------------------------- 250 subroutine poisson_psolver_end(this) 251 type(poisson_psolver_t), intent(inout) :: this 252 253 PUSH_SUB(poisson_psolver_end) 254 255#if (defined HAVE_LIBISF) || (defined HAVE_PSOLVER) 256 call pkernel_free(this%kernel) 257 call f_lib_finalize() 258#endif 259#ifdef HAVE_PSOLVER 260 call dict_free(this%inputs) 261#endif 262 263 POP_SUB(poisson_psolver_end) 264 end subroutine poisson_psolver_end 265 266 !----------------------------------------------------------------- 267 subroutine poisson_psolver_reinit(this, mesh, cube, mu, qq_in) 268 type(poisson_psolver_t), intent(inout) :: this 269 type(cube_t), intent(inout) :: cube 270 type(mesh_t), intent(inout) :: mesh 271 FLOAT, intent(in) :: mu 272 FLOAT, intent(in) :: qq_in(1:MAX_DIM) 273 274 FLOAT :: alpha, beta, gamma 275 FLOAT :: qq_abs(1:MAX_DIM), qq(1:MAX_DIM) 276 FLOAT :: modq2 277 integer :: idim 278 279 PUSH_SUB(poisson_psolver_reinit) 280 281#ifdef HAVE_PSOLVER 282 call pkernel_free(this%kernel) 283#endif 284 285 !We might change the cell angles 286 alpha = mesh%sb%alpha*M_PI/(CNST(180.0)) 287 beta = mesh%sb%beta*M_PI/(CNST(180.0)) 288 gamma = mesh%sb%gamma*M_PI/(CNST(180.0)) 289 290#ifdef HAVE_PSOLVER 291 call dict_set(this%inputs//'kernel'//'screening',mu) 292 293 this%kernel = pkernel_init(cube%mpi_grp%rank, cube%mpi_grp%size, this%inputs,& 294 this%geocode,cube%rs_n_global,mesh%spacing, & 295 alpha_bc = alpha, beta_ac = beta, gamma_ab = gamma) 296 call pkernel_set(this%kernel,verbose=debug%info) 297#endif 298 299 !G=0 component 300 !We remove potential umklapp 301 do idim = 1, mesh%sb%periodic_dim 302 qq(idim) = qq_in(idim) - anint(qq_in(idim) + M_HALF*CNST(1e-8)) 303 end do 304 call kpoints_to_absolute(mesh%sb%klattice, qq, qq_abs, mesh%sb%periodic_dim) 305 modq2 = sum(qq_abs(1:mesh%sb%periodic_dim)**2) 306 if (modq2 > M_EPSILON) then 307 this%offset = M_ONE/modq2 308 else 309 this%offset = M_ZERO 310 end if 311 312 !Screened coulomb potential (erfc function) 313 if (mu > M_EPSILON) then 314 if (modq2 > M_EPSILON) then 315 this%offset = this%offset*(M_ONE - exp(-modq2/((M_TWO*mu)**2))) 316 else 317 !Analytical limit of 1/|q|^2*(1-exp(-|q|^2/4mu^2)) 318 this%offset = M_ONE/((M_TWO*mu)**2) 319 end if 320 end if 321 this%offset = this%offset*M_FOUR*M_PI 322 323 POP_SUB(poisson_psolver_reinit) 324 end subroutine poisson_psolver_reinit 325 326 !----------------------------------------------------------------- 327 subroutine poisson_psolver_parallel_solve(this, mesh, cube, pot, rho, mesh_cube_map) 328 type(poisson_psolver_t), intent(in), target :: this 329 type(mesh_t), intent(in) :: mesh 330 type(cube_t), intent(in) :: cube 331 FLOAT, intent(out) :: pot(:) 332 FLOAT, intent(in) :: rho(:) 333 type(mesh_cube_parallel_map_t), intent(in) :: mesh_cube_map 334 335 type(profile_t), save :: prof 336 type(cube_function_t) :: cf 337#if (defined HAVE_LIBISF) || (defined HAVE_PSOLVER) 338 double precision :: hartree_energy !< Hartree energy 339#endif 340 !> offset: Total integral on the supercell of the final potential on output. 341 !! To be used only in the periodic case, ignored for other boundary conditions. 342#ifdef HAVE_PSOLVER 343 FLOAT :: offset 344#elif HAVE_LIBISF 345 double precision :: offset 346#endif 347 348 !> pot_ion: additional external potential that is added to the output 349 !! when the XC parameter ixc/=0 and sumpion=.true. 350 !! When sumpion=.true., it is always provided in the distributed form, 351 !! clearly without the overlapping terms which are needed only for the XC part 352 double precision, allocatable :: pot_ion(:,:,:) 353 354 character(len=3) :: quiet 355 356#ifdef HAVE_PSOLVER 357 type(coulomb_operator), pointer :: kernel_pointer 358#elif HAVE_LIBISF 359 double precision :: strten(6) 360#endif 361 362 PUSH_SUB(poisson_psolver_parallel_solve) 363 364 call cube_function_null(cf) 365 call dcube_function_alloc_RS(cube, cf) 366 367 call dmesh_to_cube_parallel(mesh, rho, cube, cf, mesh_cube_map) 368 369 SAFE_ALLOCATE(pot_ion(1:cube%rs_n(1),1:cube%rs_n(2),1:cube%rs_n(3))) 370 371 if(.not.debug%info) then 372 quiet = "YES" 373 else 374 quiet = "NO " 375 end if 376 377#ifdef HAVE_PSOLVER 378 !The offset is the integral over space of the potential 379 !this%offset is the G=0 component of the (screened) Coulomb potential 380 !The G=0 component of the Hartree therefore needs to be 381 ! multiplied by the G=0 component of the density 382 if(this%offset > M_EPSILON) then 383 offset = this%offset*dmf_integrate(mesh,rho) 384 end if 385#endif 386 387 call profiling_in(prof,"PSOLVER_LIBRARY") 388#ifdef HAVE_PSOLVER 389 kernel_pointer => this%kernel 390 call H_potential(this%datacode, kernel_pointer, cf%dRS, pot_ion, hartree_energy, offset, .false., quiet = quiet) 391#elif HAVE_LIBISF 392 call H_potential(this%datacode, this%kernel, cf%dRS, pot_ion, hartree_energy, offset, .false., quiet = quiet, & 393 stress_tensor = strten) 394#endif 395 call profiling_out(prof) 396 SAFE_DEALLOCATE_A(pot_ion) 397 398 call dcube_to_mesh_parallel(cube, cf, mesh, pot, mesh_cube_map) 399 400 call dcube_function_free_RS(cube, cf) 401 402 POP_SUB(poisson_psolver_parallel_solve) 403 end subroutine poisson_psolver_parallel_solve 404 405 !----------------------------------------------------------------- 406 subroutine poisson_psolver_global_solve(this, mesh, cube, pot, rho) 407 type(poisson_psolver_t), intent(in), target :: this 408 type(mesh_t), intent(in) :: mesh 409 type(cube_t), intent(in) :: cube 410 FLOAT, intent(out) :: pot(:) 411 FLOAT, intent(in) :: rho(:) 412 413 character(len=3) :: quiet 414 type(profile_t), save :: prof 415 type(cube_function_t) :: cf 416#if (defined HAVE_LIBISF) || (defined HAVE_PSOLVER) 417 double precision :: hartree_energy !< Hartree energy 418#endif 419 !> offset: Total integral on the supercell of the final potential on output 420 !! To be used only in the periodic case, ignored for other boundary conditions. 421#ifdef HAVE_PSOLVER 422 FLOAT :: offset 423#elif HAVE_LIBISF 424 double precision :: offset 425#endif 426 427 !> pot_ion: additional external potential 428 !! that is added to the output when the XC parameter ixc/=0 and sumpion=.true. 429 !! When sumpion=.true., it is always provided in the distributed form, 430 !! clearly without the overlapping terms which are needed only for the XC part 431 double precision, allocatable :: pot_ion(:,:,:) 432 433#ifdef HAVE_PSOLVER 434 type(coulomb_operator), pointer :: kernel_pointer 435#elif HAVE_LIBISF 436 double precision :: strten(6) 437#endif 438 439 PUSH_SUB(poisson_psolver_global_solve) 440 441 call cube_function_null(cf) 442 call dcube_function_alloc_RS(cube, cf) 443 444 if(mesh%parallel_in_domains) then 445 call dmesh_to_cube(mesh, rho, cube, cf, local=.true.) 446 else 447 call dmesh_to_cube(mesh, rho, cube, cf) 448 end if 449 450 SAFE_ALLOCATE(pot_ion(1:cube%rs_n(1),1:cube%rs_n(2),1:cube%rs_n(3))) 451 452 if (.not. debug%info) then 453 quiet = "YES" 454 else 455 quiet = "NO " 456 end if 457 458#ifdef HAVE_PSOLVER 459 !The offset is the integral over space of the potential 460 !this%offset is the G=0 component of the (screened) Coulomb potential 461 !The G=0 component of the Hartree therefore needs to be 462 ! multiplied by the G=0 component of the density 463 if(this%offset > M_EPSILON) then 464 offset = this%offset*dmf_integrate(mesh, rho) 465 end if 466#endif 467 468 call profiling_in(prof,"PSOLVER_LIBRARY") 469 470#ifdef HAVE_PSOLVER 471 kernel_pointer => this%kernel 472 call H_potential(this%datacode, kernel_pointer, & 473 cf%dRS, pot_ion, hartree_energy, offset, .false., & 474 quiet = quiet) !optional argument 475#elif HAVE_LIBISF 476 call H_potential(this%datacode, this%kernel, & 477 cf%dRS, pot_ion, hartree_energy, offset, .false., & 478 quiet = quiet, stress_tensor = strten) !optional argument 479#endif 480 SAFE_DEALLOCATE_A(pot_ion) 481 482 if(mesh%parallel_in_domains) then 483 call dcube_to_mesh(cube, cf, mesh, pot, local=.true.) 484 else 485 call dcube_to_mesh(cube, cf, mesh, pot) 486 end if 487 488 call dcube_function_free_RS(cube, cf) 489 490 POP_SUB(poisson_psolver_global_solve) 491 end subroutine poisson_psolver_global_solve 492 493 subroutine poisson_psolver_get_dims(this, cube) 494 type(poisson_psolver_t), intent(inout) :: this 495 type(cube_t), intent(inout) :: cube 496 497#if (defined HAVE_LIBISF) || (defined HAVE_PSOLVER) 498 !> ixc eXchange-Correlation code. Indicates the XC functional to be used 499 !! for calculating XC energies and potential. 500 !! ixc=0 indicates that no XC terms are computed. The XC functional codes follow 501 !! the ABINIT convention. 502 !> n3d third dimension of the density. For distributed data, it takes into account 503 !! the enlarging needed for calculating the XC functionals. 504 !! For global data it is simply equal to n03. 505 !! When there are too many processes and there is no room for the density n3d=0 506 integer :: n3d 507 !> n3p third dimension for the potential. The same as n3d, but without 508 !! taking into account the enlargment for the XC part. For non-GGA XC, n3p=n3d. 509 integer :: n3p 510 !> n3pi Dimension of the pot_ion array, always with distributed data. 511 !! For distributed data n3pi=n3p 512 integer :: n3pi 513 !> i3xcsh Shift of the density that must be performed to enter in the 514 !! non-overlapping region. Useful for recovering the values of the potential 515 !! when using GGA XC functionals. If the density starts from rhopot(1,1,1), 516 !! the potential starts from rhopot(1,1,i3xcsh+1). 517 !! For non-GGA XCs and for global distribution data i3xcsh=0 518 integer :: i3xcsh 519 !> i3s Starting point of the density effectively treated by each processor 520 !! in the third direction. 521 !! It takes into account also the XC enlarging. The array rhopot will correspond 522 !! To the planes of third coordinate from i3s to i3s+n3d-1. 523 !! The potential to the planes from i3s+i3xcsh to i3s+i3xcsh+n3p-1 524 !! The array pot_ion to the planes from i3s+i3xcsh to i3s+i3xcsh+n3pi-1 525 !! For global disposition i3s is equal to distributed case with i3xcsh=0. 526 integer :: i3s 527 528 !> use_gradient: .true. if functional is using the gradient. 529 logical :: use_gradient = .false. 530 !> use_wb_corr: .true. if functional is using WB corrections. 531 logical :: use_wb_corr = .false. 532#endif 533 534 PUSH_SUB(poisson_psolver_get_dims) 535 536 !! Get the dimensions of the cube 537 538#ifdef HAVE_PSOLVER 539 call PS_dim4allocation(this%geocode, this%datacode, cube%mpi_grp%rank, cube%mpi_grp%size, & 540 cube%rs_n_global(1), cube%rs_n_global(2), cube%rs_n_global(3), & 541 use_gradient, use_wb_corr, & 542 0, n3d, n3p, n3pi, i3xcsh, i3s) 543 this%localnscatterarr(:) = (/ n3d, n3p, n3pi, i3xcsh, i3s /) 544#elif HAVE_LIBISF 545 call PS_dim4allocation(this%geocode, this%datacode, cube%mpi_grp%rank, cube%mpi_grp%size, & 546 cube%rs_n_global(1), cube%rs_n_global(2), cube%rs_n_global(3), & 547 use_gradient, use_wb_corr, & 548 n3d, n3p, n3pi, i3xcsh, i3s) 549 this%localnscatterarr(:) = (/ n3d, n3p, n3pi, i3xcsh, i3s /) 550#endif 551 cube%rs_n(1:2) = cube%rs_n_global(1:2) 552 cube%rs_n(3) = this%localnscatterarr(1) 553 cube%rs_istart(1:2) = 1 554 cube%rs_istart(3) = this%localnscatterarr(5) 555 556 !! With PSolver we don`t care about the Fourier space and its dimensions 557 !! We`ll put as in RS 558 cube%fs_n_global(1) = cube%rs_n_global(1) 559 cube%fs_n_global(2) = cube%rs_n_global(2) 560 cube%fs_n_global(3) = cube%rs_n_global(3) 561 cube%fs_n(1:2) = cube%rs_n_global(1:2) 562 cube%fs_n(3) = this%localnscatterarr(1) 563 cube%fs_istart(1:2) = 1 564 cube%fs_istart(3) = this%localnscatterarr(5) 565 566 POP_SUB(poisson_psolver_get_dims) 567 end subroutine poisson_psolver_get_dims 568 569end module poisson_psolver_oct_m 570 571!! Local Variables: 572!! mode: f90 573!! coding: utf-8 574!! End: 575