1!! Copyright (C) 2002-2011 M. Marques, A. Castro, A. Rubio, 2!! G. Bertsch, M. Oliveira 3!! 4!! This program is free software; you can redistribute it and/or modify 5!! it under the terms of the GNU General Public License as published by 6!! the Free Software Foundation; either version 2, or (at your option) 7!! any later version. 8!! 9!! This program is distributed in the hope that it will be useful, 10!! but WITHOUT ANY WARRANTY; without even the implied warranty of 11!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 12!! GNU General Public License for more details. 13!! 14!! You should have received a copy of the GNU General Public License 15!! along with this program; if not, write to the Free Software 16!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 17!! 02110-1301, USA. 18!! 19 20#include "global.h" 21 22module poisson_oct_m 23 use batch_oct_m 24 use comm_oct_m 25 use cube_oct_m 26 use cube_function_oct_m 27 use derivatives_oct_m 28 use fft_oct_m 29 use fourier_space_oct_m 30 use global_oct_m 31 use index_oct_m 32 use io_oct_m 33 use io_function_oct_m 34 use loct_math_oct_m 35 use mesh_oct_m 36 use mesh_cube_parallel_map_oct_m 37 use mesh_function_oct_m 38 use messages_oct_m 39 use mpi_oct_m 40 use multicomm_oct_m 41 use namespace_oct_m 42#ifdef HAVE_OPENMP 43 use omp_lib 44#endif 45 use par_vec_oct_m 46 use parser_oct_m 47 use partition_oct_m 48 use photon_mode_oct_m 49 use poisson_cg_oct_m 50 use poisson_corrections_oct_m 51 use poisson_isf_oct_m 52 use poisson_fft_oct_m 53 use poisson_fmm_oct_m 54 use poisson_psolver_oct_m 55 use poisson_multigrid_oct_m 56 use poisson_no_oct_m 57 use profiling_oct_m 58 use simul_box_oct_m 59 use submesh_oct_m 60 use types_oct_m 61 use unit_oct_m 62 use unit_system_oct_m 63 use varinfo_oct_m 64 65#ifdef HAVE_POKE 66 use poke 67#endif 68 69 implicit none 70 71 private 72 public :: & 73 poisson_t, & 74 poisson_fmm_t, & 75 poisson_get_solver, & 76 poisson_init, & 77 poisson_init_sm, & 78 dpoisson_solve, & 79 zpoisson_solve, & 80 dpoisson_solve_sm, & 81 zpoisson_solve_sm, & 82 poisson_solve_batch, & 83 poisson_solver_is_iterative, & 84 poisson_solver_has_free_bc, & 85 poisson_end, & 86 poisson_test, & 87 poisson_is_multigrid, & 88 poisson_slave_work, & 89 poisson_async_init, & 90 poisson_async_end, & 91 dpoisson_solve_start, & 92 dpoisson_solve_finish, & 93 zpoisson_solve_start, & 94 zpoisson_solve_finish, & 95 poisson_build_kernel, & 96 poisson_is_async 97 98 integer, public, parameter :: & 99 POISSON_DIRECT_SUM = -1, & 100 POISSON_FMM = -4, & 101 POISSON_FFT = 0, & 102 POISSON_CG = 5, & 103 POISSON_CG_CORRECTED = 6, & 104 POISSON_MULTIGRID = 7, & 105 POISSON_ISF = 8, & 106 POISSON_PSOLVER = 10, & 107 POISSON_POKE = 11, & 108 POISSON_NO = -99, & 109 POISSON_NULL = -999 110 111 type poisson_t 112 private 113 type(derivatives_t), pointer, public :: der 114 integer, public :: method = POISSON_NULL 115 integer, public :: kernel 116 type(cube_t), public :: cube 117 type(mesh_cube_parallel_map_t), public :: mesh_cube_map 118 type(mg_solver_t) :: mg 119 type(poisson_fft_t), public :: fft_solver 120 FLOAT, public :: poisson_soft_coulomb_param 121 logical :: all_nodes_default 122 type(poisson_corr_t) :: corrector 123 type(poisson_isf_t) :: isf_solver 124 type(poisson_psolver_t) :: psolver_solver 125 type(poisson_no_t) :: no_solver 126 integer :: nslaves 127 logical, public :: is_dressed = .false. 128 type(photon_mode_t), public :: photons 129 type(poisson_fmm_t) :: params_fmm 130#ifdef HAVE_MPI2 131 integer :: intercomm 132 type(mpi_grp_t) :: local_grp 133 logical :: root 134#endif 135#ifdef HAVE_POKE 136 type(PokeGrid) :: poke_grid 137 type(PokeSolver) :: poke_solver 138#endif 139 end type poisson_t 140 141 integer, parameter :: & 142 CMD_FINISH = 1, & 143 CMD_POISSON_SOLVE = 2 144 145contains 146 147 !----------------------------------------------------------------- 148 subroutine poisson_init(this, namespace, der, mc, qtot, label, solver, verbose, force_serial, force_cmplx) 149 type(poisson_t), intent(out) :: this 150 type(namespace_t), intent(in) :: namespace 151 type(derivatives_t), target, intent(in) :: der 152 type(multicomm_t), intent(in) :: mc 153 FLOAT, intent(in) :: qtot !< total charge 154 character(len=*), optional, intent(in) :: label 155 integer, optional, intent(in) :: solver 156 logical, optional, intent(in) :: verbose 157 logical, optional, intent(in) :: force_serial 158 logical, optional, intent(in) :: force_cmplx 159 160 logical :: need_cube, isf_data_is_parallel 161 integer :: default_solver, default_kernel, box(MAX_DIM), fft_type, fft_library 162 FLOAT :: fft_alpha 163 character(len=60) :: str 164 165 if(this%method /= POISSON_NULL) return ! already initialized 166 167 PUSH_SUB(poisson_init) 168 169 if(optional_default(verbose,.true.)) then 170 str = "Hartree" 171 if(present(label)) str = trim(str) // trim(label) 172 call messages_print_stress(stdout, trim(str)) 173 end if 174 175 this%nslaves = 0 176 this%der => der 177 178 !%Variable DressedOrbitals 179 !%Type logical 180 !%Default false 181 !%Section Hamiltonian::Poisson 182 !%Description 183 !% Allows for the calculation of coupled elecron-photon problems 184 !% by applying the dressed orbital approach. Details can be found in 185 !% https://arxiv.org/abs/1812.05562 186 !% At the moment, N electrons in d (<=3) spatial dimensions, coupled 187 !% to one photon mode can be described. The photon mode is included by 188 !% raising the orbital dimension to d+1 and changing the particle interaction 189 !% kernel and the local potential, where the former is included automatically, 190 !% but the latter needs to by added by hand as a user_defined_potential! 191 !% Coordinate 1-d: electron; coordinate d+1: photon. 192 !%End 193 call parse_variable(namespace, 'DressedOrbitals', .false., this%is_dressed) 194 call messages_print_var_value(stdout, 'DressedOrbitals', this%is_dressed) 195 if (this%is_dressed) then 196 call messages_experimental('Dressed Orbitals') 197 ASSERT(qtot > M_ZERO) 198 call photon_mode_init(this%photons, namespace, der%mesh, der%mesh%sb%dim-1, qtot) 199 if (this%photons%nmodes > 1) then 200 call messages_not_implemented('DressedOrbitals for more than one photon mode.') 201 end if 202 end if 203 204#ifdef HAVE_MPI 205 if(.not.optional_default(force_serial,.false.)) then 206 !%Variable ParallelizationPoissonAllNodes 207 !%Type logical 208 !%Default true 209 !%Section Execution::Parallelization 210 !%Description 211 !% When running in parallel, this variable selects whether the 212 !% Poisson solver should divide the work among all nodes or only 213 !% among the parallelization-in-domains groups. 214 !%End 215 216 call parse_variable(namespace, 'ParallelizationPoissonAllNodes', .true., this%all_nodes_default) 217 else 218 this%all_nodes_default = .false. 219 end if 220#endif 221 222 !%Variable PoissonSolver 223 !%Type integer 224 !%Section Hamiltonian::Poisson 225 !%Description 226 !% Defines which method to use to solve the Poisson equation. Some incompatibilities apply depending on 227 !% dimensionality, periodicity, etc. 228 !% For a comparison of the accuracy and performance of the methods in Octopus, see P Garcia-Risueño, 229 !% J Alberdi-Rodriguez <i>et al.</i>, <i>J. Comp. Chem.</i> <b>35</b>, 427-444 (2014) 230 !% or <a href=http://arxiv.org/abs/1211.2092>arXiV</a>. 231 !% Defaults: 232 !% <br> 1D and 2D: <tt>fft</tt>. 233 !% <br> 3D: <tt>cg_corrected</tt> if curvilinear, <tt>isf</tt> if not periodic, <tt>fft</tt> if periodic. 234 !% <br> Dressed orbitals: <tt>direct_sum</tt>. 235 !%Option NoPoisson -99 236 !% Do not use a Poisson solver at all. 237 !%Option FMM -4 238 !% (Experimental) Fast multipole method. Requires FMM library. 239 !%Option direct_sum -1 240 !% Direct evaluation of the Hartree potential (only for finite systems). 241 !%Option fft 0 242 !% The Poisson equation is solved using FFTs. A cutoff technique 243 !% for the Poisson kernel is selected so the proper boundary 244 !% conditions are imposed according to the periodicity of the 245 !% system. This can be overridden by the <tt>PoissonFFTKernel</tt> 246 !% variable. To choose the FFT library use <tt>FFTLibrary</tt> 247 !%Option cg 5 248 !% Conjugate gradients (only for finite systems). 249 !%Option cg_corrected 6 250 !% Conjugate gradients, corrected for boundary conditions (only for finite systems). 251 !%Option multigrid 7 252 !% Multigrid method (only for finite systems). 253 !%Option isf 8 254 !% Interpolating Scaling Functions Poisson solver (only for finite systems). 255 !%Option psolver 10 256 !% Solver based on Interpolating Scaling Functions as implemented in the PSolver library. 257 !% Parallelization in k-points requires <tt>PoissonSolverPSolverParallelData</tt> = no. 258 !% Requires the PSolver external library. 259 !%Option poke 11 260 !% (Experimental) Solver from the Poke library. 261 !%End 262 263 default_solver = POISSON_FFT 264 265 if(der%mesh%sb%dim == 3 .and. der%mesh%sb%periodic_dim == 0) default_solver = POISSON_ISF 266 267 if(der%mesh%sb%dim > 3) default_solver = POISSON_CG_CORRECTED 268 269#ifdef HAVE_CLFFT 270 ! this is disabled, since the difference between solvers are big 271 ! enough to cause problems with the tests. 272 ! if(accel_is_enabled()) default_solver = POISSON_FFT 273#endif 274 275 if(der%mesh%use_curvilinear) then 276 select case(der%mesh%sb%dim) 277 case(1) 278 default_solver = POISSON_DIRECT_SUM 279 case(2) 280 default_solver = POISSON_DIRECT_SUM 281 case(3) 282 default_solver = POISSON_CG_CORRECTED 283 end select 284 end if 285 286 if (this%is_dressed) default_solver = POISSON_DIRECT_SUM 287 288 if(.not.present(solver)) then 289 call parse_variable(namespace, 'PoissonSolver', default_solver, this%method) 290 else 291 this%method = solver 292 end if 293 if(.not.varinfo_valid_option('PoissonSolver', this%method)) call messages_input_error(namespace, 'PoissonSolver') 294 if(optional_default(verbose,.true.)) then 295 select case(this%method) 296 case (POISSON_DIRECT_SUM) 297 str = "direct sum" 298 case (POISSON_FMM) 299 str = "fast multipole method" 300 case (POISSON_FFT) 301 str = "fast Fourier transform" 302 case (POISSON_CG) 303 str = "conjugate gradients" 304 case (POISSON_CG_CORRECTED) 305 str = "conjugate gradients, corrected" 306 case (POISSON_MULTIGRID) 307 str = "multigrid" 308 case (POISSON_ISF) 309 str = "interpolating scaling functions" 310 case (POISSON_PSOLVER) 311 str = "interpolating scaling functions (from BigDFT)" 312 case (POISSON_NO) 313 str = "no Poisson solver - Hartree set to 0" 314 case (POISSON_POKE) 315 str = "Poke library" 316 end select 317 write(message(1),'(a,a,a)') "The chosen Poisson solver is '", trim(str), "'" 318 call messages_info(1) 319 end if 320 321 if(this%method /= POISSON_FFT) then 322 this%kernel = POISSON_FFT_KERNEL_NONE 323 else 324 325 ! Documentation in cube.F90 326 call parse_variable(namespace, 'FFTLibrary', FFTLIB_FFTW, fft_library) 327 328 !%Variable PoissonFFTKernel 329 !%Type integer 330 !%Section Hamiltonian::Poisson 331 !%Description 332 !% Defines which kernel is used to impose the correct boundary 333 !% conditions when using FFTs to solve the Poisson equation. The 334 !% default is selected depending on the dimensionality and 335 !% periodicity of the system: 336 !% <br>In 1D, <tt>spherical</tt> if finite, <tt>fft_nocut</tt> if periodic. 337 !% <br>In 2D, <tt>spherical</tt> if finite, <tt>cylindrical</tt> if 1D-periodic, <tt>fft_nocut</tt> if 2D-periodic. 338 !% <br>In 3D, <tt>spherical</tt> if finite, <tt>cylindrical</tt> if 1D-periodic, <tt>planar</tt> if 2D-periodic, 339 !% <tt>fft_nocut</tt> if 3D-periodic. 340 !% See C. A. Rozzi et al., <i>Phys. Rev. B</i> <b>73</b>, 205119 (2006) for 3D implementation and 341 !% A. Castro et al., <i>Phys. Rev. B</i> <b>80</b>, 033102 (2009) for 2D implementation. 342 !%Option spherical 0 343 !% FFTs using spherical cutoff (in 2D or 3D). 344 !%Option cylindrical 1 345 !% FFTs using cylindrical cutoff (in 2D or 3D). 346 !%Option planar 2 347 !% FFTs using planar cutoff (in 3D). 348 !%Option fft_nocut 3 349 !% FFTs without using a cutoff (for fully periodic systems). 350 !%Option multipole_correction 4 351 !% The boundary conditions are imposed by using a multipole expansion. Only appropriate for finite systems. 352 !% Further specification occurs with variables <tt>PoissonSolverBoundaries</tt> and <tt>PoissonSolverMaxMultipole</tt>. 353 !%End 354 355 select case(der%mesh%sb%dim) 356 case(1) 357 if(der%mesh%sb%periodic_dim == 0) then 358 default_kernel = POISSON_FFT_KERNEL_SPH 359 else 360 default_kernel = POISSON_FFT_KERNEL_NOCUT 361 end if 362 case(2) 363 if (der%mesh%sb%periodic_dim == 2) then 364 default_kernel = POISSON_FFT_KERNEL_NOCUT 365 else if (der%mesh%sb%periodic_dim > 0) then 366 default_kernel = der%mesh%sb%periodic_dim 367 else 368 default_kernel = POISSON_FFT_KERNEL_SPH 369 end if 370 case(3) 371 default_kernel = der%mesh%sb%periodic_dim 372 end select 373 374 call parse_variable(namespace, 'PoissonFFTKernel', default_kernel, this%kernel) 375 if(.not.varinfo_valid_option('PoissonFFTKernel', this%kernel)) call messages_input_error(namespace, 'PoissonFFTKernel') 376 377 if(optional_default(verbose,.true.)) & 378 call messages_print_var_option(stdout, "PoissonFFTKernel", this%kernel) 379 380 end if 381 382 !We assume the developer knows what he is doing by providing the solver option 383 if(.not. present(solver)) then 384 if(der%mesh%sb%periodic_dim > 0 .and. this%method == POISSON_DIRECT_SUM) then 385 message(1) = 'A periodic system may not use the direct_sum Poisson solver.' 386 call messages_fatal(1) 387 end if 388 389 if(der%mesh%sb%periodic_dim > 0 .and. this%method == POISSON_CG_CORRECTED) then 390 message(1) = 'A periodic system may not use the cg_corrected Poisson solver.' 391 call messages_fatal(1) 392 end if 393 394 if(der%mesh%sb%periodic_dim > 0 .and. this%method == POISSON_CG) then 395 message(1) = 'A periodic system may not use the cg Poisson solver.' 396 call messages_fatal(1) 397 end if 398 399 if(der%mesh%sb%periodic_dim > 0 .and. this%method == POISSON_MULTIGRID) then 400 message(1) = 'A periodic system may not use the multigrid Poisson solver.' 401 call messages_fatal(1) 402 end if 403 404 select case(der%mesh%sb%dim) 405 case(1) 406 407 select case(der%mesh%sb%periodic_dim) 408 case(0) 409 if( (this%method /= POISSON_FFT) .and. (this%method /= POISSON_DIRECT_SUM)) then 410 message(1) = 'A finite 1D system may only use fft or direct_sum Poisson solvers.' 411 call messages_fatal(1) 412 end if 413 case(1) 414 if(this%method /= POISSON_FFT) then 415 message(1) = 'A periodic 1D system may only use the fft Poisson solver.' 416 call messages_fatal(1) 417 end if 418 end select 419 420 if(der%mesh%use_curvilinear .and. this%method /= POISSON_DIRECT_SUM) then 421 message(1) = 'If curvilinear coordinates are used in 1D, then the only working' 422 message(2) = 'Poisson solver is direct_sum.' 423 call messages_fatal(2) 424 end if 425 426 case(2) 427 428 if ((this%method /= POISSON_FFT) .and. (this%method /= POISSON_DIRECT_SUM)) then 429 message(1) = 'A 2D system may only use fft or direct_sum solvers.' 430 call messages_fatal(1) 431 end if 432 433 if(der%mesh%use_curvilinear .and. (this%method /= POISSON_DIRECT_SUM) ) then 434 message(1) = 'If curvilinear coordinates are used in 2D, then the only working' 435 message(2) = 'Poisson solver is direct_sum.' 436 call messages_fatal(2) 437 end if 438 439 case(3) 440 441 if(der%mesh%sb%periodic_dim > 0 .and. this%method == POISSON_FMM) then 442 call messages_not_implemented('FMM for periodic systems') 443 end if 444 445 if(der%mesh%sb%periodic_dim > 0 .and. this%method == POISSON_ISF) then 446 call messages_write('The ISF solver can only be used for finite systems.') 447 call messages_fatal() 448 end if 449 450 if(der%mesh%sb%periodic_dim > 0 .and. this%method == POISSON_FFT .and. & 451 this%kernel /= der%mesh%sb%periodic_dim .and. this%kernel >=0 .and. this%kernel <=3) then 452 write(message(1), '(a,i1,a)')'The system is periodic in ', der%mesh%sb%periodic_dim ,' dimension(s),' 453 write(message(2), '(a,i1,a)')'but Poisson solver is set for ', this%kernel, ' dimensions.' 454 call messages_warning(2) 455 end if 456 457 if(der%mesh%sb%periodic_dim > 0 .and. this%method == POISSON_FFT .and. & 458 this%kernel == POISSON_FFT_KERNEL_CORRECTED) then 459 write(message(1), '(a,i1,a)')'PoissonFFTKernel = multipole_correction cannot be used for periodic systems.' 460 call messages_fatal(1) 461 end if 462 463 if(der%mesh%use_curvilinear .and. (this%method/=POISSON_CG_CORRECTED)) then 464 message(1) = 'If curvilinear coordinates are used, then the only working' 465 message(2) = 'Poisson solver is cg_corrected.' 466 call messages_fatal(2) 467 end if 468 469 if( (der%mesh%sb%box_shape == MINIMUM) .and. (this%method == POISSON_CG_CORRECTED) ) then 470 message(1) = 'When using the "minimum" box shape and the "cg_corrected"' 471 message(2) = 'Poisson solver, we have observed "sometimes" some non-' 472 message(3) = 'negligible error. You may want to check that the "fft" or "cg"' 473 message(4) = 'solver are providing, in your case, the same results.' 474 call messages_warning(4) 475 end if 476 477 if (this%method == POISSON_FMM) then 478 call messages_experimental('FMM Poisson solver') 479 end if 480 end select 481 end if 482 483 if (this%method == POISSON_PSOLVER) then 484#if !((defined HAVE_LIBISF) || (defined HAVE_PSOLVER)) 485 message(1) = "The PSolver Poisson solver cannot be used since the code was not compiled with the PSolver libary." 486 call messages_fatal(1) 487#endif 488#ifdef HAVE_LIBISF 489 message(1) = "The use of versions older than 1.8 of the PSolver library (previously known as LibISF)" 490 message(2) = "are deprecated and will be removed in the next major release." 491 call messages_warning(2) 492#endif 493 end if 494 495 if(optional_default(verbose,.true.)) & 496 call messages_print_stress(stdout) 497 498 ! Now that we know the method, we check if we need a cube and its dimentions 499 need_cube = .false. 500 fft_type = FFT_REAL 501 if(optional_default(force_cmplx, .false.)) fft_type = FFT_COMPLEX 502 503 if (this%method == POISSON_ISF .or. this%method == POISSON_PSOLVER) then 504 fft_type = FFT_NONE 505 box(:) = der%mesh%idx%ll(:) 506 need_cube = .true. 507 end if 508 509 if (this%method == POISSON_PSOLVER .and. multicomm_have_slaves(mc)) then 510 call messages_not_implemented('Task parallelization with LibISF Poisson solver') 511 end if 512 513 if ( multicomm_strategy_is_parallel(mc, P_STRATEGY_KPOINTS) ) then 514 ! Documentation in poisson_psolver.F90 515 call parse_variable(namespace, 'PoissonSolverPSolverParallelData', .true., isf_data_is_parallel) 516 if ( this%method == POISSON_PSOLVER .and. isf_data_is_parallel ) then 517 call messages_not_implemented("k-point parallelization with PSolver library and PoissonSolverPSolverParallelData = yes") 518 end if 519 if ( this%method == POISSON_FFT .and. fft_library == FFTLIB_PFFT ) then 520 call messages_not_implemented("k-point parallelization with PFFT library for Poisson solver") 521 end if 522 end if 523 524 if (this%method == POISSON_FFT) then 525 526 need_cube = .true. 527 528 !%Variable DoubleFFTParameter 529 !%Type float 530 !%Default 2.0 531 !%Section Mesh::FFTs 532 !%Description 533 !% For solving the Poisson equation in Fourier space, and for applying the local potential 534 !% in Fourier space, an auxiliary cubic mesh is built. This mesh will be larger than 535 !% the circumscribed cube of the usual mesh by a factor <tt>DoubleFFTParameter</tt>. See 536 !% the section that refers to Poisson equation, and to the local potential for details 537 !% [the default value of two is typically good]. 538 !%End 539 call parse_variable(namespace, 'DoubleFFTParameter', M_TWO, fft_alpha) 540 if (fft_alpha < M_ONE .or. fft_alpha > M_THREE ) then 541 write(message(1), '(a,f12.5,a)') "Input: '", fft_alpha, & 542 "' is not a valid DoubleFFTParameter" 543 message(2) = '1.0 <= DoubleFFTParameter <= 3.0' 544 call messages_fatal(2) 545 end if 546 547 if (der%mesh%sb%dim /= 3 .and. fft_library == FFTLIB_PFFT) then 548 call messages_not_implemented('PFFT support for dimensionality other than 3') 549 end if 550 if (der%mesh%sb%periodic_dim /= 0 .and. fft_library == FFTLIB_PFFT) then 551 call messages_not_implemented('PFFT support for periodic systems') 552 end if 553 554 select case (der%mesh%sb%dim) 555 556 case (1) 557 select case(this%kernel) 558 case(POISSON_FFT_KERNEL_SPH) 559 call mesh_double_box(der%mesh%sb, der%mesh, fft_alpha, box) 560 case(POISSON_FFT_KERNEL_NOCUT) 561 box = der%mesh%idx%ll 562 end select 563 564 case (2) 565 select case(this%kernel) 566 case(POISSON_FFT_KERNEL_SPH) 567 call mesh_double_box(der%mesh%sb, der%mesh, fft_alpha, box) 568 box(1:2) = maxval(box) 569 case(POISSON_FFT_KERNEL_CYL) 570 call mesh_double_box(der%mesh%sb, der%mesh, fft_alpha, box) 571 case(POISSON_FFT_KERNEL_NOCUT) 572 box(:) = der%mesh%idx%ll(:) 573 end select 574 575 case (3) 576 select case(this%kernel) 577 case(POISSON_FFT_KERNEL_SPH) 578 call mesh_double_box(der%mesh%sb, der%mesh, fft_alpha, box) 579 box(:) = maxval(box) 580 case(POISSON_FFT_KERNEL_CYL) 581 call mesh_double_box(der%mesh%sb, der%mesh, fft_alpha, box) 582 box(2) = maxval(box(2:3)) ! max of finite directions 583 box(3) = maxval(box(2:3)) ! max of finite directions 584 case(POISSON_FFT_KERNEL_CORRECTED) 585 box(:) = der%mesh%idx%ll(:) 586 case(POISSON_FFT_KERNEL_PLA, POISSON_FFT_KERNEL_NOCUT) 587 call mesh_double_box(der%mesh%sb, der%mesh, fft_alpha, box) 588 end select 589 590 end select 591 592 end if 593 594 if(this%method == POISSON_POKE) then 595#ifndef HAVE_POKE 596 call messages_write('Octopus was compiled without Poke support, you cannot use', new_line = .true.) 597 call messages_write(" 'PoissonSolver = poke'. ") 598 call messages_fatal() 599#endif 600 601 call messages_experimental('Poke library') 602 ASSERT(der%mesh%sb%dim == 3) 603 box(1:der%mesh%sb%dim) = der%mesh%idx%ll(1:der%mesh%sb%dim) 604 need_cube = .true. 605 fft_type = FFTLIB_NONE 606 end if 607 608 ! Create the cube 609 if (need_cube) then 610 call cube_init(this%cube, box, der%mesh%sb, namespace, fft_type = fft_type, & 611 need_partition=.not.der%mesh%parallel_in_domains) 612 if (this%cube%parallel_in_domains .and. this%method == POISSON_FFT) then 613 call mesh_cube_parallel_map_init(this%mesh_cube_map, der%mesh, this%cube) 614 end if 615 end if 616 617 if(this%method == POISSON_POKE) then 618 619#ifdef HAVE_POKE 620 this%poke_grid = PokeGrid(der%mesh%spacing, this%cube%rs_n) 621 if(der%mesh%sb%periodic_dim > 0) then 622 call this%poke_grid%set_boundaries(POKE_BOUNDARIES_PERIODIC) 623 else 624 call this%poke_grid%set_boundaries(POKE_BOUNDARIES_FREE) 625 end if 626 this%poke_solver = PokeSolver(this%poke_grid) 627 call this%poke_solver%build() 628#endif 629 end if 630 631 if (this%is_dressed .and. .not. this%method == POISSON_DIRECT_SUM) then 632 write(message(1), '(a)')'Dressed Orbital calculation currently only working with direct sum Poisson solver.' 633 call messages_fatal(1) 634 end if 635 636 call poisson_kernel_init(this, namespace, mc%master_comm) 637 638 POP_SUB(poisson_init) 639 end subroutine poisson_init 640 641 !----------------------------------------------------------------- 642 subroutine poisson_end(this) 643 type(poisson_t), intent(inout) :: this 644 645 logical :: has_cube 646 647 PUSH_SUB(poisson_end) 648 649 has_cube = .false. 650 651 select case(this%method) 652 case(POISSON_FFT) 653 call poisson_fft_end(this%fft_solver) 654 if(this%kernel == POISSON_FFT_KERNEL_CORRECTED) call poisson_corrections_end(this%corrector) 655 has_cube = .true. 656 657 case(POISSON_CG_CORRECTED, POISSON_CG) 658 call poisson_cg_end() 659 call poisson_corrections_end(this%corrector) 660 661 case(POISSON_MULTIGRID) 662 call poisson_multigrid_end(this%mg) 663 664 case(POISSON_ISF) 665 call poisson_isf_end(this%isf_solver) 666 has_cube = .true. 667 668 case(POISSON_PSOLVER) 669 call poisson_psolver_end(this%psolver_solver) 670 has_cube = .true. 671 672 case(POISSON_FMM) 673 call poisson_fmm_end(this%params_fmm) 674 675 case(POISSON_NO) 676 call poisson_no_end(this%no_solver) 677 678 case(POISSON_POKE) 679#ifdef HAVE_POKE 680 call this%poke_grid%end() 681 call this%poke_solver%end() 682#endif 683 684 end select 685 this%method = POISSON_NULL 686 687 if (has_cube) then 688 if (this%cube%parallel_in_domains) then 689 call mesh_cube_parallel_map_end(this%mesh_cube_map) 690 end if 691 call cube_end(this%cube) 692 end if 693 694 if (this%is_dressed) then 695 call photon_mode_end(this%photons) 696 end if 697 this%is_dressed = .false. 698 699 POP_SUB(poisson_end) 700 end subroutine poisson_end 701 702 !----------------------------------------------------------------- 703 704 subroutine zpoisson_solve_real_and_imag_separately(this, pot, rho, all_nodes, kernel) 705 type(poisson_t), intent(in) :: this 706 CMPLX, intent(inout) :: pot(:) !< pot(mesh%np) 707 CMPLX, intent(in) :: rho(:) !< rho(mesh%np) 708 logical, optional, intent(in) :: all_nodes 709 type(fourier_space_op_t), optional, intent(in) :: kernel 710 711 FLOAT, allocatable :: aux1(:), aux2(:) 712 type(derivatives_t), pointer :: der 713 logical :: all_nodes_value 714 715 type(profile_t), save :: prof 716 717 der => this%der 718 719 PUSH_SUB(zpoisson_solve_real_and_imag_separately) 720 721 call profiling_in(prof, 'POISSON_RE_IM_SOLVE') 722 723 if(present(kernel)) then 724 ASSERT(.not. any(abs(kernel%qq(:))>CNST(1e-8))) 725 end if 726 727 all_nodes_value = optional_default(all_nodes, this%all_nodes_default) 728 729 SAFE_ALLOCATE(aux1(1:der%mesh%np)) 730 SAFE_ALLOCATE(aux2(1:der%mesh%np)) 731 ! first the real part 732 aux1(1:der%mesh%np) = real(rho(1:der%mesh%np)) 733 aux2(1:der%mesh%np) = real(pot(1:der%mesh%np)) 734 call dpoisson_solve(this, aux2, aux1, all_nodes=all_nodes_value, kernel=kernel) 735 pot(1:der%mesh%np) = aux2(1:der%mesh%np) 736 737 ! now the imaginary part 738 aux1(1:der%mesh%np) = aimag(rho(1:der%mesh%np)) 739 aux2(1:der%mesh%np) = aimag(pot(1:der%mesh%np)) 740 call dpoisson_solve(this, aux2, aux1, all_nodes=all_nodes_value, kernel=kernel) 741 pot(1:der%mesh%np) = pot(1:der%mesh%np) + M_zI*aux2(1:der%mesh%np) 742 743 SAFE_DEALLOCATE_A(aux1) 744 SAFE_DEALLOCATE_A(aux2) 745 746 call profiling_out(prof) 747 748 POP_SUB(zpoisson_solve_real_and_imag_separately) 749 end subroutine zpoisson_solve_real_and_imag_separately 750 751 !----------------------------------------------------------------- 752 753 subroutine zpoisson_solve(this, pot, rho, all_nodes, kernel) 754 type(poisson_t), intent(in) :: this 755 CMPLX, intent(inout) :: pot(:) !< pot(mesh%np) 756 CMPLX, intent(in) :: rho(:) !< rho(mesh%np) 757 logical, optional, intent(in) :: all_nodes 758 type(fourier_space_op_t), optional, intent(in) :: kernel 759 760 logical :: all_nodes_value 761 type(profile_t), save :: prof 762 763 PUSH_SUB(zpoisson_solve) 764 765 all_nodes_value = optional_default(all_nodes, this%all_nodes_default) 766 767 ASSERT(ubound(pot, dim = 1) == this%der%mesh%np_part .or. ubound(pot, dim = 1) == this%der%mesh%np) 768 ASSERT(ubound(rho, dim = 1) == this%der%mesh%np_part .or. ubound(rho, dim = 1) == this%der%mesh%np) 769 770 ASSERT(this%method /= POISSON_NULL) 771 772 if(this%method == POISSON_FFT .and. this%kernel /= POISSON_FFT_KERNEL_CORRECTED & 773 .and. .not. this%is_dressed) then 774 !The default (real) Poisson solver is used for OEP and Sternheimer calls were we do not need 775 !a complex-to-xomplex FFT as these parts use the normal Coulomb potential 776 if(this%cube%fft%type == FFT_COMPLEX) then 777 !We add the profiling here, as the other path uses dpoisson_solve 778 call profiling_in(prof, 'ZPOISSON_SOLVE') 779 call zpoisson_fft_solve(this%fft_solver, this%der%mesh, this%cube, pot, rho, this%mesh_cube_map, kernel=kernel) 780 call profiling_out(prof) 781 else 782 call zpoisson_solve_real_and_imag_separately(this, pot, rho, all_nodes_value, kernel=kernel) 783 end if 784 else 785 call zpoisson_solve_real_and_imag_separately(this, pot, rho, all_nodes_value, kernel = kernel) 786 end if 787 788 POP_SUB(zpoisson_solve) 789 end subroutine zpoisson_solve 790 791 792 !----------------------------------------------------------------- 793 794 subroutine poisson_solve_batch(this, potb, rhob, all_nodes, kernel) 795 type(poisson_t), intent(inout) :: this 796 type(batch_t), intent(inout) :: potb 797 type(batch_t), intent(inout) :: rhob 798 logical, optional, intent(in) :: all_nodes 799 type(fourier_space_op_t), optional, intent(in) :: kernel 800 801 integer :: ii 802 803 PUSH_SUB(poisson_solve_batch) 804 805 ASSERT(potb%nst_linear == rhob%nst_linear) 806 ASSERT(potb%type() == rhob%type()) 807 808 if(potb%type() == TYPE_FLOAT) then 809 do ii = 1, potb%nst_linear 810 call dpoisson_solve(this, potb%dff_linear(:, ii), rhob%dff_linear(:, ii), all_nodes, kernel=kernel) 811 end do 812 else 813 do ii = 1, potb%nst_linear 814 call zpoisson_solve(this, potb%zff_linear(:, ii), rhob%zff_linear(:, ii), all_nodes, kernel=kernel) 815 end do 816 end if 817 818 POP_SUB(poisson_solve_batch) 819 end subroutine poisson_solve_batch 820 821 !----------------------------------------------------------------- 822 823 !> Calculates the Poisson equation. 824 !! Given the density returns the corresponding potential. 825 !! 826 !! Different solvers are available that can be chosen in the input file 827 !! with the "PoissonSolver" parameter 828 subroutine dpoisson_solve(this, pot, rho, all_nodes, kernel) 829 type(poisson_t), intent(in) :: this 830 FLOAT, intent(inout) :: pot(:) !< Local size of the \b potential vector. 831 FLOAT, intent(inout) :: rho(:) !< Local size of the \b density (rho) vector. 832 !> Is the Poisson solver allowed to utilise 833 !! all nodes or only the domain nodes for 834 !! its calculations? (Defaults to .true.) 835 logical, optional, intent(in) :: all_nodes 836 type(fourier_space_op_t), optional, intent(in) :: kernel 837 838 type(derivatives_t), pointer :: der 839 type(cube_function_t) :: crho, cpot 840 FLOAT, allocatable :: rho_corrected(:), vh_correction(:) 841 logical :: all_nodes_value 842 type(profile_t), save :: prof 843 844 call profiling_in(prof, 'POISSON_SOLVE') 845 PUSH_SUB(dpoisson_solve) 846 847 der => this%der 848 849 ASSERT(ubound(pot, dim = 1) == der%mesh%np_part .or. ubound(pot, dim = 1) == der%mesh%np) 850 ASSERT(ubound(rho, dim = 1) == der%mesh%np_part .or. ubound(rho, dim = 1) == der%mesh%np) 851 852 ! Check optional argument and set to default if necessary. 853 all_nodes_value = optional_default(all_nodes, this%all_nodes_default) 854 855 ASSERT(this%method /= POISSON_NULL) 856 857 if(present(kernel)) then 858 ASSERT(this%method == POISSON_FFT) 859 end if 860 861 select case(this%method) 862 case(POISSON_DIRECT_SUM) 863 if ( (this%is_dressed .and. this%der%mesh%sb%dim - 1 > 3) .or. this%der%mesh%sb%dim > 3) then 864 message(1) = "Direct sum Poisson solver only available for 1, 2, or 3 dimensions." 865 call messages_fatal(1) 866 end if 867 call poisson_solve_direct(this, pot, rho) 868 869 case(POISSON_FMM) 870 call poisson_fmm_solve(this%params_fmm, this%der, pot, rho) 871 872 case(POISSON_CG) 873 call poisson_cg1(der, this%corrector, pot, rho) 874 875 case(POISSON_CG_CORRECTED) 876 SAFE_ALLOCATE(rho_corrected(1:der%mesh%np)) 877 SAFE_ALLOCATE(vh_correction(1:der%mesh%np_part)) 878 879 call correct_rho(this%corrector, der, rho, rho_corrected, vh_correction) 880 881 pot(1:der%mesh%np) = pot(1:der%mesh%np) - vh_correction(1:der%mesh%np) 882 call poisson_cg2(der, pot, rho_corrected) 883 pot(1:der%mesh%np) = pot(1:der%mesh%np) + vh_correction(1:der%mesh%np) 884 885 SAFE_DEALLOCATE_A(rho_corrected) 886 SAFE_DEALLOCATE_A(vh_correction) 887 888 case(POISSON_MULTIGRID) 889 call poisson_multigrid_solver(this%mg, der, pot, rho) 890 891 case(POISSON_FFT) 892 if(this%kernel /= POISSON_FFT_KERNEL_CORRECTED) then 893 call dpoisson_fft_solve(this%fft_solver, der%mesh, this%cube, pot, rho, this%mesh_cube_map, kernel=kernel) 894 else 895 SAFE_ALLOCATE(rho_corrected(1:der%mesh%np)) 896 SAFE_ALLOCATE(vh_correction(1:der%mesh%np_part)) 897 898 call correct_rho(this%corrector, der, rho, rho_corrected, vh_correction) 899 call dpoisson_fft_solve(this%fft_solver, der%mesh, this%cube, pot, rho_corrected, this%mesh_cube_map, & 900 average_to_zero = .true., kernel=kernel) 901 902 pot(1:der%mesh%np) = pot(1:der%mesh%np) + vh_correction(1:der%mesh%np) 903 SAFE_DEALLOCATE_A(rho_corrected) 904 SAFE_DEALLOCATE_A(vh_correction) 905 end if 906 907 case(POISSON_ISF) 908 call poisson_isf_solve(this%isf_solver, der%mesh, this%cube, pot, rho, all_nodes_value) 909 910 911 case(POISSON_PSOLVER) 912 if (this%psolver_solver%datacode == "G") then 913 ! Global version 914 call poisson_psolver_global_solve(this%psolver_solver, der%mesh, this%cube, pot, rho) 915 else ! "D" Distributed version 916 call poisson_psolver_parallel_solve(this%psolver_solver, der%mesh, this%cube, pot, rho, this%mesh_cube_map) 917 end if 918 919 case(POISSON_POKE) 920 call cube_function_null(crho) 921 call cube_function_null(cpot) 922 call dcube_function_alloc_RS(this%cube, crho) 923 call dcube_function_alloc_RS(this%cube, cpot) 924 call dmesh_to_cube(der%mesh, rho, this%cube, crho) 925#if HAVE_POKE 926 call this%poke_solver%solve(crho%drs, cpot%drs) 927#endif 928 call dcube_to_mesh(this%cube, cpot, der%mesh, pot) 929 call dcube_function_free_RS(this%cube, crho) 930 call dcube_function_free_RS(this%cube, cpot) 931 932 case(POISSON_NO) 933 call poisson_no_solve(this%no_solver, der%mesh, this%cube, pot, rho) 934 end select 935 936 937 ! Add extra terms for dressed interaction 938 if (this%is_dressed .and. this%method /= POISSON_NO) then 939 call photon_mode_add_poisson_terms(this%photons, der%mesh, rho, pot) 940 end if 941 942 POP_SUB(dpoisson_solve) 943 call profiling_out(prof) 944 end subroutine dpoisson_solve 945 946 !----------------------------------------------------------------- 947 subroutine poisson_init_sm(this, namespace, main, der, sm, method) 948 type(poisson_t), intent(out) :: this 949 type(namespace_t), intent(in) :: namespace 950 type(poisson_t), intent(in) :: main 951 type(derivatives_t), target, intent(in) :: der 952 type(submesh_t), intent(inout) :: sm 953 integer, optional, intent(in) :: method 954 955 integer :: default_solver 956 integer :: box(MAX_DIM) 957 958 if(this%method /= POISSON_NULL) return ! already initialized 959 960 PUSH_SUB(poisson_init_sm) 961 962 this%is_dressed = .false. 963 964 this%nslaves = 0 965 this%der => der 966 967#ifdef HAVE_MPI 968 this%all_nodes_default = main%all_nodes_default 969#endif 970 971 default_solver = POISSON_DIRECT_SUM 972 this%method = default_solver 973 if(present(method)) this%method = method 974 975 if(der%mesh%use_curvilinear) then 976 call messages_not_implemented("Submesh Poisson solver with curvilinear mesh") 977 end if 978 979 this%kernel = POISSON_FFT_KERNEL_NONE 980 981 nullify(sm%cube_map%map) 982 983 select case(this%method) 984 case(POISSON_DIRECT_SUM) 985 !Nothing to be done 986 case(POISSON_ISF) 987 !TODO: Add support for domain parrallelization 988 ASSERT(.not. der%mesh%parallel_in_domains) 989 call submesh_get_cube_dim(sm, box, der%dim) 990 call submesh_init_cube_map(sm, der%dim) 991 call cube_init(this%cube, box, der%mesh%sb, namespace, fft_type = FFT_NONE, & 992 need_partition=.not.der%mesh%parallel_in_domains) 993 call poisson_isf_init(this%isf_solver, namespace, der%mesh, this%cube, mpi_world%comm, init_world = this%all_nodes_default) 994 end select 995 996 POP_SUB(poisson_init_sm) 997 end subroutine poisson_init_sm 998 999 !----------------------------------------------------------------- 1000 !> This routine checks the Hartree solver selected in the input 1001 !! file by calculating numerically and analytically the Hartree 1002 !! potential originated by a Gaussian distribution of charge. 1003 !! This only makes sense for finite systems. 1004 subroutine poisson_test(this, mesh, namespace, repetitions) 1005 type(poisson_t), intent(in) :: this 1006 type(mesh_t), intent(in) :: mesh 1007 type(namespace_t), intent(in) :: namespace 1008 integer, intent(in) :: repetitions 1009 1010 FLOAT, allocatable :: rho(:), vh(:), vh_exact(:), rhop(:), xx(:, :) 1011 FLOAT :: alpha, beta, rr, delta, ralpha, hartree_nrg_num, & 1012 hartree_nrg_analyt, lcl_hartree_nrg 1013 FLOAT :: total_charge 1014 integer :: ip, idir, ierr, iunit, nn, n_gaussians, itime 1015 1016 PUSH_SUB(poisson_test) 1017 1018 if(mesh%sb%dim == 1) then 1019 call messages_not_implemented('Poisson test for 1D case') 1020 end if 1021 1022 n_gaussians = 1 1023 1024 SAFE_ALLOCATE( rho(1:mesh%np)) 1025 SAFE_ALLOCATE( rhop(1:mesh%np)) 1026 SAFE_ALLOCATE( vh(1:mesh%np)) 1027 SAFE_ALLOCATE(vh_exact(1:mesh%np)) 1028 SAFE_ALLOCATE(xx(1:mesh%sb%dim, 1:n_gaussians)) 1029 1030 rho = M_ZERO; vh = M_ZERO; vh_exact = M_ZERO; rhop = M_ZERO 1031 1032 alpha = CNST(4.0)*mesh%spacing(1) 1033 write(message(1),'(a,f14.6)') "Info: The alpha value is ", alpha 1034 write(message(2),'(a)') " Higher values of alpha lead to more physical densities and more reliable results." 1035 call messages_info(2) 1036 beta = M_ONE / ( alpha**mesh%sb%dim * sqrt(M_PI)**mesh%sb%dim ) 1037 1038 write(message(1), '(a)') 'Building the Gaussian distribution of charge...' 1039 call messages_info(1) 1040 1041 rho = M_ZERO 1042 do nn = 1, n_gaussians 1043 do idir = 1, mesh%sb%dim 1044 xx(idir, nn) = M_ZERO 1045 end do 1046 1047 rr = sqrt(sum(xx(:, nn)*xx(:,nn))) 1048 do ip = 1, mesh%np 1049 call mesh_r(mesh, ip, rr, origin = xx(:, nn)) 1050 rhop(ip) = beta*exp(-(rr/alpha)**2) 1051 end do 1052 1053 rhop = (-1)**nn * rhop 1054 do ip = 1, mesh%np 1055 rho(ip) = rho(ip) + rhop(ip) 1056 end do 1057 end do 1058 1059 total_charge = dmf_integrate(mesh, rho) 1060 1061 write(message(1), '(a,f14.6)') 'Total charge of the Gaussian distribution', total_charge 1062 call messages_info(1) 1063 1064 ! This builds analytically its potential 1065 vh_exact = M_ZERO 1066 do nn = 1, n_gaussians 1067 do ip = 1, mesh%np 1068 call mesh_r(mesh, ip, rr, origin = xx(:, nn)) 1069 select case(mesh%sb%dim) 1070 case(3) 1071 if(rr > R_SMALL) then 1072 vh_exact(ip) = vh_exact(ip) + (-1)**nn * loct_erf(rr/alpha)/rr 1073 else 1074 vh_exact(ip) = vh_exact(ip) + (-1)**nn * (M_TWO/sqrt(M_PI))/alpha 1075 end if 1076 case(2) 1077 ralpha = rr**2/(M_TWO*alpha**2) 1078 if(ralpha < CNST(100.0)) then 1079 vh_exact(ip) = vh_exact(ip) + (-1)**nn * beta * (M_PI)**(M_THREE*M_HALF) * alpha * exp(-rr**2/(M_TWO*alpha**2)) * & 1080 loct_bessel_in(0, rr**2/(M_TWO*alpha**2)) 1081 else 1082 vh_exact(ip) = vh_exact(ip) + (-1)**nn * beta * (M_PI)**(M_THREE*M_HALF) * alpha * & 1083 (M_ONE/sqrt(M_TWO*M_PI*ralpha)) 1084 end if 1085 end select 1086 end do 1087 end do 1088 1089 ! This calculates the numerical potential 1090 do itime = 1, repetitions 1091 call dpoisson_solve(this, vh, rho) 1092 end do 1093 1094 ! Output results 1095 iunit = io_open("hartree_results", namespace, action='write') 1096 delta = dmf_nrm2(mesh, vh-vh_exact) 1097 write(iunit, '(a,f19.13)' ) 'Hartree test (abs.) = ', delta 1098 delta = delta/dmf_nrm2(mesh, vh_exact) 1099 write(iunit, '(a,f19.13)' ) 'Hartree test (rel.) = ', delta 1100 1101 ! Calculate the numerical Hartree energy (serially) 1102 lcl_hartree_nrg = M_ZERO 1103 do ip = 1, mesh%np 1104 lcl_hartree_nrg = lcl_hartree_nrg + rho(ip) * vh(ip) 1105 end do 1106 lcl_hartree_nrg = lcl_hartree_nrg * mesh%spacing(1) * mesh%spacing(2) * mesh%spacing(3)/M_TWO 1107#ifdef HAVE_MPI 1108 call MPI_Reduce(lcl_hartree_nrg, hartree_nrg_num, 1, & 1109 MPI_FLOAT, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) 1110 if(mpi_err /= 0) then 1111 write(message(1),'(a)') "MPI error in MPI_Reduce; subroutine poisson_test of file poisson.F90" 1112 call messages_warning(1) 1113 end if 1114#else 1115 hartree_nrg_num = lcl_hartree_nrg 1116#endif 1117 1118 ! Calculate the analytical Hartree energy (serially, discrete - not exactly exact) 1119 lcl_hartree_nrg = M_ZERO 1120 do ip = 1, mesh%np 1121 lcl_hartree_nrg = lcl_hartree_nrg + rho(ip) * vh_exact(ip) 1122 end do 1123 lcl_hartree_nrg = lcl_hartree_nrg * mesh%spacing(1) * mesh%spacing(2) * mesh%spacing(3)/M_TWO 1124#ifdef HAVE_MPI 1125 call MPI_Reduce(lcl_hartree_nrg, hartree_nrg_analyt, 1, & 1126 MPI_FLOAT, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) 1127 if(mpi_err /= 0) then 1128 write(message(1),'(a)') "MPI error in MPI_Reduce; subroutine poisson_test of file poisson.F90" 1129 call messages_warning(1) 1130 end if 1131#else 1132 hartree_nrg_analyt = lcl_hartree_nrg 1133#endif 1134 1135 write(iunit, '(a,f19.13)' ) 1136 1137 if (mpi_world%rank == 0) then 1138 write(iunit,'(a,f19.13)') 'Hartree Energy (numerical) =',hartree_nrg_num,'Hartree Energy (analytical) =',hartree_nrg_analyt 1139 end if 1140 1141 call io_close(iunit) 1142 1143 call dio_function_output (io_function_fill_how('AxisX'), ".", "poisson_test_rho", namespace, & 1144 mesh, rho, unit_one, ierr) 1145 call dio_function_output (io_function_fill_how('AxisX'), ".", "poisson_test_exact", namespace, & 1146 mesh, vh_exact, unit_one, ierr) 1147 call dio_function_output (io_function_fill_how('AxisX'), ".", "poisson_test_numerical", namespace, & 1148 mesh, vh, unit_one, ierr) 1149 call dio_function_output (io_function_fill_how('AxisY'), ".", "poisson_test_rho", namespace, & 1150 mesh, rho, unit_one, ierr) 1151 call dio_function_output (io_function_fill_how('AxisY'), ".", "poisson_test_exact", namespace, & 1152 mesh, vh_exact, unit_one, ierr) 1153 call dio_function_output (io_function_fill_how('AxisY'), ".", "poisson_test_numerical", namespace, & 1154 mesh, vh, unit_one, ierr) 1155 ! not dimensionless, but no need for unit conversion for a test routine 1156 1157 SAFE_DEALLOCATE_A(rho) 1158 SAFE_DEALLOCATE_A(rhop) 1159 SAFE_DEALLOCATE_A(vh) 1160 SAFE_DEALLOCATE_A(vh_exact) 1161 SAFE_DEALLOCATE_A(xx) 1162 1163 POP_SUB(poisson_test) 1164 end subroutine poisson_test 1165 1166 ! ----------------------------------------------------------------- 1167 1168 logical pure function poisson_solver_is_iterative(this) result(iterative) 1169 type(poisson_t), intent(in) :: this 1170 1171 iterative = this%method == POISSON_CG .or. this%method == POISSON_CG_CORRECTED .or. this%method == POISSON_MULTIGRID 1172 end function poisson_solver_is_iterative 1173 1174 ! ----------------------------------------------------------------- 1175 1176 logical pure function poisson_is_multigrid(this) result(is_multigrid) 1177 type(poisson_t), intent(in) :: this 1178 1179 is_multigrid = (this%method == POISSON_MULTIGRID) 1180 1181 end function poisson_is_multigrid 1182 1183 ! ----------------------------------------------------------------- 1184 1185 logical pure function poisson_solver_has_free_bc(this) result(free_bc) 1186 type(poisson_t), intent(in) :: this 1187 1188 free_bc = .true. 1189 1190 if (this%method == POISSON_FFT .and. & 1191 this%kernel /= POISSON_FFT_KERNEL_SPH .and. this%kernel /= POISSON_FFT_KERNEL_CORRECTED) then 1192 free_bc = .false. 1193 end if 1194 1195 end function poisson_solver_has_free_bc 1196 1197 !----------------------------------------------------------------- 1198 1199 integer pure function poisson_get_solver(this) result (solver) 1200 type(poisson_t), intent(in) :: this 1201 1202 solver = this%method 1203 end function poisson_get_solver 1204 1205 !----------------------------------------------------------------- 1206 1207 subroutine poisson_async_init(this, mc) 1208 type(poisson_t), intent(inout) :: this 1209 type(multicomm_t), intent(in) :: mc 1210 1211 PUSH_SUB(poisson_async_init) 1212 1213#ifdef HAVE_MPI2 1214 if(multicomm_have_slaves(mc)) then 1215 1216 call mpi_grp_init(this%local_grp, mc%group_comm(P_STRATEGY_STATES)) 1217 1218 this%root = (this%local_grp%rank == 0) 1219 1220 this%intercomm = mc%slave_intercomm 1221 call MPI_Comm_remote_size(this%intercomm, this%nslaves, mpi_err) 1222 1223 end if 1224#endif 1225 1226 POP_SUB(poisson_async_init) 1227 1228 end subroutine poisson_async_init 1229 1230 !----------------------------------------------------------------- 1231 1232 subroutine poisson_async_end(this, mc) 1233 type(poisson_t), intent(inout) :: this 1234 type(multicomm_t), intent(in) :: mc 1235 1236#ifdef HAVE_MPI2 1237 integer :: islave 1238#endif 1239 1240 PUSH_SUB(poisson_async_end) 1241 1242#ifdef HAVE_MPI2 1243 if(multicomm_have_slaves(mc)) then 1244 1245 ! send the finish signal 1246 do islave = this%local_grp%rank, this%nslaves - 1, this%local_grp%size 1247 call MPI_Send(M_ONE, 1, MPI_FLOAT, islave, CMD_FINISH, this%intercomm, mpi_err) 1248 end do 1249 1250 end if 1251#endif 1252 1253 POP_SUB(poisson_async_end) 1254 1255 end subroutine poisson_async_end 1256 1257 !----------------------------------------------------------------- 1258 1259 subroutine poisson_slave_work(this) 1260 type(poisson_t), intent(inout) :: this 1261 1262#ifdef HAVE_MPI2 1263 FLOAT, allocatable :: rho(:), pot(:) 1264 logical :: done 1265 integer :: status(MPI_STATUS_SIZE) 1266 type(profile_t), save :: prof, bcast_prof, wait_prof 1267 integer :: bcast_root 1268 1269 PUSH_SUB(poisson_slave_work) 1270 call profiling_in(prof, "SLAVE_WORK") 1271 1272 SAFE_ALLOCATE(rho(1:this%der%mesh%np)) 1273 SAFE_ALLOCATE(pot(1:this%der%mesh%np)) 1274 done = .false. 1275 1276 do while(.not. done) 1277 1278 call profiling_in(wait_prof, "SLAVE_WAIT") 1279 call MPI_Recv(rho(1), this%der%mesh%np, MPI_FLOAT, MPI_ANY_SOURCE, MPI_ANY_TAG, this%intercomm, status(1), mpi_err) 1280 call profiling_out(wait_prof) 1281 1282 ! The tag of the message tells us what we have to do. 1283 select case(status(MPI_TAG)) 1284 1285 case(CMD_FINISH) 1286 done = .true. 1287 1288 case(CMD_POISSON_SOLVE) 1289 call dpoisson_solve(this, pot, rho) 1290 1291 call profiling_in(bcast_prof, "SLAVE_BROADCAST") 1292 bcast_root = MPI_PROC_NULL 1293 if(this%root) bcast_root = MPI_ROOT 1294 call MPI_Bcast(pot(1), this%der%mesh%np, MPI_FLOAT, bcast_root, this%intercomm, mpi_err) 1295 call profiling_out(bcast_prof) 1296 1297 end select 1298 1299 end do 1300 1301 SAFE_DEALLOCATE_A(pot) 1302 SAFE_DEALLOCATE_A(rho) 1303 1304 call profiling_out(prof) 1305 POP_SUB(poisson_slave_work) 1306#endif 1307 end subroutine poisson_slave_work 1308 1309 !---------------------------------------------------------------- 1310 1311 logical pure function poisson_is_async(this) result(async) 1312 type(poisson_t), intent(in) :: this 1313 1314 async = (this%nslaves > 0) 1315 1316 end function poisson_is_async 1317 1318 !---------------------------------------------------------------- 1319 1320 subroutine poisson_build_kernel(this, namespace, sb, coulb, qq, mu, singul) 1321 type(poisson_t), intent(in) :: this 1322 type(namespace_t),intent(in) :: namespace 1323 type(simul_box_t),intent(in) :: sb 1324 type(fourier_space_op_t), intent(inout) :: coulb 1325 FLOAT, intent(in) :: qq(:) 1326 FLOAT, intent(in) :: mu 1327 FLOAT, optional, intent(in) :: singul 1328 1329 PUSH_SUB(poisson_build_kernel) 1330 1331 if(simul_box_is_periodic(sb)) then 1332 ASSERT(ubound(qq, 1) >= sb%periodic_dim) 1333 ASSERT(this%method == POISSON_FFT) 1334 end if 1335 1336 if(mu > M_EPSILON) then 1337 if(this%method /= POISSON_FFT) then 1338 write(message(1),'(a)') "Poisson solver with range separation is only implemented with FFT." 1339 call messages_fatal(1) 1340 end if 1341 coulb%mu = mu 1342 end if 1343 1344 !TODO: this should be a select case supporting other kernels. 1345 ! This means that we need an abstract object for kernels. 1346 select case(this%method) 1347 case(POISSON_FFT) 1348 !We only reinitialize the poisson sover if needed 1349 if(any(abs(coulb%qq(1:sb%periodic_dim) - qq(1:sb%periodic_dim)) > M_EPSILON)) then 1350 call fourier_space_op_end(coulb) 1351 coulb%qq(1:sb%periodic_dim) = qq(1:sb%periodic_dim) 1352 !We must define the singularity if we specify a q vector and we do not use the short-range Coulomb potential 1353 coulb%singularity = optional_default(singul, M_ZERO) 1354 call poisson_fft_get_kernel(namespace, this%der%mesh, this%cube, coulb, this%kernel, & 1355 this%poisson_soft_coulomb_param) 1356 end if 1357 case default 1358 call messages_not_implemented("poisson_build_kernel with other methods than FFT") 1359 end select 1360 1361 1362 POP_SUB(poisson_build_kernel) 1363 end subroutine poisson_build_kernel 1364 1365#include "poisson_init_inc.F90" 1366#include "poisson_direct_inc.F90" 1367#include "poisson_direct_sm_inc.F90" 1368 1369#include "undef.F90" 1370#include "real.F90" 1371#include "poisson_inc.F90" 1372#include "undef.F90" 1373#include "complex.F90" 1374#include "poisson_inc.F90" 1375 1376end module poisson_oct_m 1377 1378!! Local Variables: 1379!! mode: f90 1380!! coding: utf-8 1381!! End: 1382