1!! Copyright (C) 2002-2011 M. Marques, A. Castro, A. Rubio, G. Bertsch, M. Oliveira 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_fft_oct_m 22 use cube_function_oct_m 23 use cube_oct_m 24 use fft_oct_m 25 use fourier_space_oct_m 26 use global_oct_m 27 use loct_math_oct_m 28 use math_oct_m 29 use mesh_cube_parallel_map_oct_m 30 use mesh_oct_m 31 use messages_oct_m 32 use namespace_oct_m 33 use parser_oct_m 34 use poisson_cutoff_oct_m 35 use profiling_oct_m 36 use simul_box_oct_m 37 use splines_oct_m 38 use unit_oct_m 39 use unit_system_oct_m 40 41 implicit none 42 private 43 public :: & 44 poisson_fft_t, & 45 poisson_fft_init, & 46 poisson_fft_get_kernel, & 47 poisson_fft_end, & 48 dpoisson_fft_solve, & 49 zpoisson_fft_solve 50 51 integer, public, parameter :: & 52 POISSON_FFT_KERNEL_NONE = -1, & 53 POISSON_FFT_KERNEL_SPH = 0, & 54 POISSON_FFT_KERNEL_CYL = 1, & 55 POISSON_FFT_KERNEL_PLA = 2, & 56 POISSON_FFT_KERNEL_NOCUT = 3, & 57 POISSON_FFT_KERNEL_CORRECTED = 4, & 58 POISSON_FFT_KERNEL_HOCKNEY = 5 59 60 type poisson_fft_t 61 ! Components are public by default 62 type(fourier_space_op_t) :: coulb !< object for Fourier space operations 63 integer :: kernel !< choice of kernel, one of options above 64 FLOAT :: soft_coulb_param !< Soft-Coulomb parameter 65 end type poisson_fft_t 66contains 67 68 subroutine poisson_fft_init(this, namespace, mesh, cube, kernel, soft_coulb_param, fullcube) 69 type(poisson_fft_t), intent(out) :: this 70 type(namespace_t), intent(in) :: namespace 71 type(mesh_t), intent(in) :: mesh 72 type(cube_t), intent(inout) :: cube 73 integer, intent(in) :: kernel 74 FLOAT, optional, intent(in) :: soft_coulb_param 75 type(cube_t), optional, intent(in) :: fullcube !< needed for Hockney kernel 76 77 PUSH_SUB(poisson_fft_init) 78 79 this%kernel = kernel 80 this%soft_coulb_param = optional_default(soft_coulb_param, M_ZERO) 81 82 this%coulb%qq(1:mesh%sb%periodic_dim) = M_ZERO 83 this%coulb%singularity = M_ZERO 84 this%coulb%mu = M_ZERO 85 86 call poisson_fft_get_kernel(namespace, mesh, cube, this%coulb, kernel, soft_coulb_param, fullcube) 87 88 POP_SUB(poisson_fft_init) 89 end subroutine poisson_fft_init 90 91 subroutine poisson_fft_get_kernel(namespace, mesh, cube, coulb, kernel, soft_coulb_param, fullcube) 92 type(namespace_t), intent(in) :: namespace 93 type(mesh_t), intent(in) :: mesh 94 type(cube_t), intent(in) :: cube 95 type(fourier_space_op_t), intent(inout) :: coulb 96 integer, intent(in) :: kernel 97 FLOAT, optional, intent(in) :: soft_coulb_param 98 type(cube_t), optional, intent(in) :: fullcube !< needed for Hockney kernel 99 100 PUSH_SUB(poisson_fft_get_kernel) 101 102 if(coulb%mu > M_EPSILON) then 103 if(mesh%sb%dim /= 3 .or. kernel /= POISSON_FFT_KERNEL_NOCUT) then 104 message(1) = "The screened Coulomb potential is only implemented in 3D for PoissonFFTKernel=fft_nocut." 105 call messages_fatal(1) 106 end if 107 end if 108 109 110 if(kernel == POISSON_FFT_KERNEL_HOCKNEY) then 111 if (.not. present(fullcube)) then 112 message(1) = "Hockney's FFT-kernel needs cube of full unit cell " 113 call messages_fatal(1) 114 else 115 if (.not.associated(fullcube%fft)) then 116 message(1) = "Hockney's FFT-kernel needs PoissonSolver=fft" 117 call messages_fatal(1) 118 end if 119 end if 120 end if 121 122 123 select case(mesh%sb%dim) 124 case(1) 125 ASSERT(present(soft_coulb_param)) 126 select case(kernel) 127 case(POISSON_FFT_KERNEL_SPH) 128 call poisson_fft_build_1d_0d(namespace, mesh, cube, coulb, soft_coulb_param) 129 case(POISSON_FFT_KERNEL_NOCUT) 130 call poisson_fft_build_1d_1d(mesh, cube, coulb, soft_coulb_param) 131 case default 132 message(1) = "Invalid Poisson FFT kernel for 1D." 133 call messages_fatal(1) 134 end select 135 136 case(2) 137 select case(kernel) 138 case(POISSON_FFT_KERNEL_SPH) 139 call poisson_fft_build_2d_0d(namespace, mesh, cube, coulb) 140 case(POISSON_FFT_KERNEL_CYL) 141 call poisson_fft_build_2d_1d(namespace, mesh, cube, coulb) 142 case(POISSON_FFT_KERNEL_NOCUT) 143 call poisson_fft_build_2d_2d(mesh, cube, coulb) 144 case default 145 message(1) = "Invalid Poisson FFT kernel for 2D." 146 call messages_fatal(1) 147 end select 148 149 case(3) 150 select case(kernel) 151 case(POISSON_FFT_KERNEL_SPH, POISSON_FFT_KERNEL_CORRECTED) 152 call poisson_fft_build_3d_0d(namespace, mesh, cube, kernel, coulb) 153 154 case(POISSON_FFT_KERNEL_CYL) 155 call poisson_fft_build_3d_1d(namespace, mesh, cube, coulb) 156 157 case(POISSON_FFT_KERNEL_PLA) 158 call poisson_fft_build_3d_2d(namespace, mesh, cube, coulb) 159 160 case(POISSON_FFT_KERNEL_NOCUT) 161 call poisson_fft_build_3d_3d(mesh, cube, coulb) 162 163 case(POISSON_FFT_KERNEL_HOCKNEY) 164 call poisson_fft_build_3d_3d_hockney(mesh, cube, coulb, fullcube) 165 166 case default 167 message(1) = "Invalid Poisson FFT kernel for 3D." 168 call messages_fatal(1) 169 end select 170 end select 171 172 POP_SUB(poisson_fft_get_kernel) 173 end subroutine poisson_fft_get_kernel 174 175 !----------------------------------------------------------------- 176 177 subroutine get_cutoff(namespace, default_r_c, r_c) 178 type(namespace_t), intent(in) :: namespace 179 FLOAT, intent(in) :: default_r_c 180 FLOAT, intent(out) :: r_c 181 182 PUSH_SUB(get_cutoff) 183 184 call parse_variable(namespace, 'PoissonCutoffRadius', default_r_c, r_c, units_inp%length) 185 186 call messages_write('Info: Poisson Cutoff Radius =') 187 call messages_write(r_c, units = units_out%length, fmt = '(f6.1)') 188 call messages_info() 189 190 if ( r_c > default_r_c + M_EPSILON) then 191 call messages_write('Poisson cutoff radius is larger than cell size.', new_line = .true.) 192 call messages_write('You can see electrons in neighboring cell(s).') 193 call messages_warning() 194 end if 195 196 POP_SUB(get_cutoff) 197 end subroutine get_cutoff 198 199 !----------------------------------------------------------------- 200 subroutine poisson_fft_gg_transform(gg_in, temp, sb, qq, gg, modg2) 201 integer, intent(in) :: gg_in(:) 202 FLOAT, intent(in) :: temp(:) 203 type(simul_box_t), intent(in) :: sb 204 FLOAT, intent(in) :: qq(:) 205 FLOAT, intent(inout) :: gg(:) 206 FLOAT, intent(out) :: modg2 207 208! integer :: idir 209 210 ! no PUSH_SUB, called too frequently 211 212 gg(1:3) = gg_in(1:3) 213 gg(1:sb%periodic_dim) = gg(1:sb%periodic_dim) + qq(1:sb%periodic_dim) 214 gg(1:3) = gg(1:3) * temp(1:3) 215 gg(1:3) = matmul(sb%klattice_primitive(1:3,1:3),gg(1:3)) 216! MJV 27 jan 2015 this should not be necessary 217! do idir = 1, 3 218! gg(idir) = gg(idir) / lalg_nrm2(3, sb%klattice_primitive(1:3, idir)) 219! end do 220 221 modg2 = sum(gg(1:3)**2) 222 223 end subroutine poisson_fft_gg_transform 224 225 !----------------------------------------------------------------- 226 subroutine poisson_fft_build_3d_3d(mesh, cube, coulb) 227 type(mesh_t), intent(in) :: mesh 228 type(cube_t), intent(in) :: cube 229 type(fourier_space_op_t), intent(inout) :: coulb 230 231 integer :: ix, iy, iz, ixx(3), db(3) 232 FLOAT :: temp(3), modg2, inv_four_mu2 233 FLOAT :: gg(3) 234 FLOAT, allocatable :: fft_Coulb_FS(:,:,:) 235 236 PUSH_SUB(poisson_fft_build_3d_3d) 237 238 db(1:3) = cube%rs_n_global(1:3) 239 240 if(coulb%mu > M_EPSILON) then 241 inv_four_mu2 = M_ONE/((M_TWO*coulb%mu)**2) 242 end if 243 244 ! store the Fourier transform of the Coulomb interaction 245 SAFE_ALLOCATE(fft_Coulb_FS(1:cube%fs_n_global(1), 1:cube%fs_n_global(2), 1:cube%fs_n_global(3))) 246 fft_Coulb_FS = M_ZERO 247 248 temp(1:3) = M_TWO*M_PI/(db(1:3)*mesh%spacing(1:3)) 249 250 do ix = 1, cube%fs_n_global(1) 251 ixx(1) = pad_feq(ix, db(1), .true.) 252 do iy = 1, cube%fs_n_global(2) 253 ixx(2) = pad_feq(iy, db(2), .true.) 254 do iz = 1, cube%fs_n_global(3) 255 ixx(3) = pad_feq(iz, db(3), .true.) 256 257 call poisson_fft_gg_transform(ixx, temp, mesh%sb, coulb%qq, gg, modg2) 258 259 !HH not very elegant 260 if(cube%fft%library.eq.FFTLIB_NFFT) modg2=cube%Lfs(ix,1)**2+cube%Lfs(iy,2)**2+cube%Lfs(iz,3)**2 261 262 if(abs(modg2) > CNST(1e-6)) then 263 !Screened coulomb potential (erfc function) 264 if(coulb%mu > M_EPSILON) then 265 fft_Coulb_FS(ix, iy, iz) = M_FOUR*M_PI/modg2*(M_ONE-exp(-modg2*inv_four_mu2)) 266 else 267 fft_Coulb_FS(ix, iy, iz) = M_FOUR*M_PI/modg2 268 end if 269 else 270 !Screened coulomb potential (erfc function) 271 if(coulb%mu > M_EPSILON) then 272 !Analytical limit of 1/|q|^2*(1-exp(-|q|^2/4mu^2)) 273 fft_Coulb_FS(ix, iy, iz) = M_FOUR*M_PI*inv_four_mu2 274 else 275 !We use the user-defined value of the singularity 276 fft_Coulb_FS(ix, iy, iz) = coulb%singularity 277 end if 278 end if 279 280 end do 281 end do 282 283 end do 284 285 call dfourier_space_op_init(coulb, cube, fft_Coulb_FS) 286 287 SAFE_DEALLOCATE_A(fft_Coulb_FS) 288 POP_SUB(poisson_fft_build_3d_3d) 289 end subroutine poisson_fft_build_3d_3d 290 !----------------------------------------------------------------- 291 292 !----------------------------------------------------------------- 293 !> Kernel for Hockneys algorithm that solves the poisson equation 294 !! in a small box while respecting the periodicity of a larger box 295 !! A. Damle, L. Lin, L. Ying, JCTC, 2015 296 !! DOI: 10.1021/ct500985f, supplementary info 297 subroutine poisson_fft_build_3d_3d_hockney(mesh, cube, coulb, fullcube) 298 type(mesh_t), intent(in) :: mesh 299 type(cube_t), intent(in) :: cube 300 type(fourier_space_op_t), intent(inout) :: coulb 301 type(cube_t), intent(in) :: fullcube 302 303 integer :: ix, iy, iz, ixx(3), db(3), nfs(3), nrs(3), nfs_s(3), nrs_s(3) 304 FLOAT :: temp(3), modg2 305 FLOAT :: gg(3) 306 FLOAT, allocatable :: fft_Coulb_small_RS(:,:,:) 307 FLOAT, allocatable :: fft_Coulb_RS(:,:,:) 308 CMPLX, allocatable :: fft_Coulb_small_FS(:,:,:), fft_Coulb_FS(:,:,:) 309 310 PUSH_SUB(poisson_fft_build_3d_3d_hockney) 311 312 ! dimensions of large boxes 313 nfs(1:3) = fullcube%fs_n_global(1:3) 314 nrs(1:3) = fullcube%rs_n_global(1:3) 315 316 SAFE_ALLOCATE(fft_Coulb_FS(1:nfs(1),1:nfs(2),1:nfs(3))) 317 SAFE_ALLOCATE(fft_Coulb_RS(1:nrs(1),1:nrs(2),1:nrs(3))) 318 319 ! dimensions of small boxes x_s 320 nfs_s(1:3) = cube%fs_n_global(1:3) 321 nrs_s(1:3) = cube%rs_n_global(1:3) 322 323 SAFE_ALLOCATE(fft_Coulb_small_FS(1:nfs_s(1),1:nfs_s(2),1:nfs_s(3))) 324 SAFE_ALLOCATE(fft_Coulb_small_RS(1:nrs_s(1),1:nrs_s(2),1:nrs_s(3))) 325 326 ! build full periodic Coulomb potenital in Fourier space 327 fft_Coulb_FS = M_ZERO 328 329 db(1:3) = fullcube%rs_n_global(1:3) 330 temp(1:3) = M_TWO*M_PI/(db(1:3)*mesh%spacing(1:3)) 331 332 do ix = 1, nfs(1) 333 ixx(1) = pad_feq(ix, db(1), .true.) 334 do iy = 1, nfs(2) 335 ixx(2) = pad_feq(iy, db(2), .true.) 336 do iz = 1, nfs(3) 337 ixx(3) = pad_feq(iz, db(3), .true.) 338 339 call poisson_fft_gg_transform(ixx, temp, mesh%sb, coulb%qq, gg, modg2) 340 341 if(abs(modg2) > M_EPSILON) then 342 fft_Coulb_FS(ix, iy, iz) = M_ONE/modg2 343 else 344 fft_Coulb_FS(ix, iy, iz) = M_ZERO 345 end if 346 end do 347 end do 348 end do 349 350 do iz = 1, nfs(3) 351 do iy = 1, nfs(2) 352 do ix = 1, nfs(1) 353 fft_Coulb_FS(ix, iy, iz) = M_FOUR*M_PI*fft_Coulb_FS(ix, iy, iz) 354 end do 355 end do 356 end do 357 358 ! get periodic Coulomb potential in real space 359 call dfft_backward(fullcube%fft, fft_Coulb_FS, fft_Coulb_RS) 360 361 ! copy to small box by respecting this pattern 362 ! full periodic coulomb: |abc--------------------------xyz| 363 ! Hockney: |abcxyz| 364 do ix = 1, nrs_s(1) 365 if(ix.le.nrs_s(1)/2+1) then 366 ixx(1)=ix 367 else 368 ixx(1) = nrs(1) - nrs_s(1)+ix 369 end if 370 do iy = 1, nrs_s(2) 371 if(iy.le.nrs_s(2)/2+1) then 372 ixx(2)=iy 373 else 374 ixx(2) = nrs(2) - nrs_s(2)+iy 375 end if 376 do iz = 1, nrs_s(3) 377 if(iz.le.nrs_s(3)/2+1) then 378 ixx(3)=iz 379 else 380 ixx(3) = nrs(3) - nrs_s(3)+iz 381 end if 382 fft_Coulb_small_RS(ix,iy,iz) = fft_Coulb_RS(ixx(1),ixx(2),ixx(3)) 383 end do 384 end do 385 end do 386 ! make Hockney kernel in Fourier space 387 call dfft_forward(cube%fft, fft_Coulb_small_RS, fft_Coulb_small_FS) 388 !dummy copy for type conversion 389 fft_Coulb_small_RS(1:nfs_s(1),1:nfs_s(2),1:nfs_s(3)) = & 390 TOFLOAT( fft_Coulb_small_FS(1:nfs_s(1),1:nfs_s(2),1:nfs_s(3))) 391 392 393 call dfourier_space_op_init(coulb, cube, fft_Coulb_small_RS(1:nfs_s(1),1:nfs_s(2),1:nfs_s(3))) 394 395 SAFE_DEALLOCATE_A(fft_Coulb_FS) 396 SAFE_DEALLOCATE_A(fft_Coulb_RS) 397 SAFE_DEALLOCATE_A(fft_Coulb_small_FS) 398 SAFE_DEALLOCATE_A(fft_Coulb_small_RS) 399 400 POP_SUB(poisson_fft_build_3d_3d_hockney) 401 402 end subroutine poisson_fft_build_3d_3d_hockney 403 404 !----------------------------------------------------------------- 405 !> C. A. Rozzi et al., Phys. Rev. B 73, 205119 (2006), Table I 406 subroutine poisson_fft_build_3d_2d(namespace, mesh, cube, coulb) 407 type(namespace_t), intent(in) :: namespace 408 type(mesh_t), intent(in) :: mesh 409 type(cube_t), intent(in) :: cube 410 type(fourier_space_op_t), intent(inout) :: coulb 411 412 integer :: ix, iy, iz, ixx(3), db(3) 413 FLOAT :: temp(3), modg2 414 FLOAT :: gpar, gz, r_c, gg(3), default_r_c 415 FLOAT, allocatable :: fft_coulb_FS(:,:,:) 416 417 PUSH_SUB(poisson_fft_build_3d_2d) 418 419 db(1:3) = cube%rs_n_global(1:3) 420 421 !%Variable PoissonCutoffRadius 422 !%Type float 423 !%Section Hamiltonian::Poisson 424 !%Description 425 !% When <tt>PoissonSolver = fft</tt> and <tt>PoissonFFTKernel</tt> is neither <tt>multipole_corrections</tt> 426 !% nor <tt>fft_nocut</tt>, 427 !% this variable controls the distance after which the electron-electron interaction goes to zero. 428 !% A warning will be written if the value is too large and will cause spurious interactions between images. 429 !% The default is half of the FFT box max dimension in a finite direction. 430 !%End 431 432 default_r_c = db(3)*mesh%spacing(3)/M_TWO 433 call get_cutoff(namespace, default_r_c, r_c) 434 435 ! store the fourier transform of the Coulomb interaction 436 SAFE_ALLOCATE(fft_Coulb_FS(1:cube%fs_n_global(1), 1:cube%fs_n_global(2), 1:cube%fs_n_global(3))) 437 fft_Coulb_FS = M_ZERO 438 439 temp(1:3) = M_TWO*M_PI/(db(1:3)*mesh%spacing(1:3)) 440 441 do ix = 1, cube%fs_n_global(1) 442 ixx(1) = pad_feq(ix, db(1), .true.) 443 do iy = 1, cube%fs_n_global(2) 444 ixx(2) = pad_feq(iy, db(2), .true.) 445 do iz = 1, cube%fs_n_global(3) 446 ixx(3) = pad_feq(iz, db(3), .true.) 447 448 call poisson_fft_gg_transform(ixx, temp, mesh%sb, coulb%qq, gg, modg2) 449 450 if(abs(modg2) > M_EPSILON) then 451 gz = abs(gg(3)) 452 gpar = hypot(gg(1), gg(2)) 453 ! note: if gpar = 0, then modg2 = gz**2 454 fft_Coulb_FS(ix, iy, iz) = poisson_cutoff_3D_2D(gpar,gz,r_c)/modg2 455 else 456 fft_Coulb_FS(ix, iy, iz) = -M_HALF*r_c**2 457 end if 458 end do 459 end do 460 461 end do 462 463 do iz = 1, cube%fs_n_global(3) 464 do iy = 1, cube%fs_n_global(2) 465 do ix = 1, cube%fs_n_global(1) 466 fft_Coulb_FS(ix, iy, iz) = M_FOUR*M_PI*fft_Coulb_FS(ix, iy, iz) 467 end do 468 end do 469 end do 470 471 call dfourier_space_op_init(coulb, cube, fft_Coulb_FS) 472 473 SAFE_DEALLOCATE_A(fft_Coulb_FS) 474 POP_SUB(poisson_fft_build_3d_2d) 475 end subroutine poisson_fft_build_3d_2d 476 !----------------------------------------------------------------- 477 478 479 !----------------------------------------------------------------- 480 !> C. A. Rozzi et al., Phys. Rev. B 73, 205119 (2006), Table I 481 subroutine poisson_fft_build_3d_1d(namespace, mesh, cube, coulb) 482 type(namespace_t), intent(in) :: namespace 483 type(mesh_t), intent(in) :: mesh 484 type(cube_t), intent(in) :: cube 485 type(fourier_space_op_t), intent(inout) :: coulb 486 487 type(spline_t) :: cylinder_cutoff_f 488 FLOAT, allocatable :: x(:), y(:) 489 integer :: ix, iy, iz, ixx(3), db(3), k, ngp 490 FLOAT :: temp(3), modg2, xmax 491 FLOAT :: gperp, gx, gy, gz, r_c, gg(3), default_r_c 492 FLOAT, allocatable :: fft_coulb_FS(:,:,:) 493 494 PUSH_SUB(poisson_fft_build_3d_1d) 495 496 db(1:3) = cube%rs_n_global(1:3) 497 498 default_r_c = maxval(db(2:3)*mesh%spacing(2:3)/M_TWO) 499 call get_cutoff(namespace, default_r_c, r_c) 500 501 ! store the fourier transform of the Coulomb interaction 502 SAFE_ALLOCATE(fft_Coulb_FS(1:cube%fs_n_global(1), 1:cube%fs_n_global(2), 1:cube%fs_n_global(3))) 503 fft_Coulb_FS = M_ZERO 504 505 temp(1:3) = M_TWO*M_PI/(db(1:3)*mesh%spacing(1:3)) 506 507 if( mesh%sb%periodic_dim == 0 ) then 508 ngp = 8*db(2) 509 SAFE_ALLOCATE(x(1:ngp)) 510 SAFE_ALLOCATE(y(1:ngp)) 511 end if 512 513 514 do ix = 1, cube%fs_n_global(1) 515 ixx(1) = pad_feq(ix, db(1), .true.) 516 gx = temp(1)*ixx(1) 517 518 if( mesh%sb%periodic_dim == 0 ) then 519 call spline_init(cylinder_cutoff_f) 520 xmax = sqrt((temp(2)*db(2)/2)**2 + (temp(3)*db(3)/2)**2) 521 do k = 1, ngp 522 x(k) = (k-1)*(xmax/(ngp-1)) 523 y(k) = poisson_cutoff_3D_1D_finite(gx, x(k), M_TWO*mesh%sb%xsize, M_TWO*mesh%sb%rsize) 524 end do 525 call spline_fit(ngp, x, y, cylinder_cutoff_f) 526 end if 527 528 do iy = 1, cube%fs_n_global(2) 529 ixx(2) = pad_feq(iy, db(2), .true.) 530 do iz = 1, db(3) 531 ixx(3) = pad_feq(iz, db(3), .true.) 532 533 call poisson_fft_gg_transform(ixx, temp, mesh%sb, coulb%qq, gg, modg2) 534 535 if(abs(modg2) > M_EPSILON) then 536 gperp = hypot(gg(2), gg(3)) 537 if (mesh%sb%periodic_dim==1) then 538 fft_Coulb_FS(ix, iy, iz) = poisson_cutoff_3D_1D(abs(gx), gperp, r_c)/modg2 539 else if (mesh%sb%periodic_dim==0) then 540 gy = gg(2) 541 gz = gg(3) 542 if ((gz >= M_ZERO) .and. (gy >= M_ZERO)) then 543 fft_Coulb_FS(ix, iy, iz) = spline_eval(cylinder_cutoff_f, gperp) 544 end if 545 if ((gz >= M_ZERO) .and. (gy < M_ZERO)) then 546 fft_Coulb_FS(ix, iy, iz) = fft_Coulb_FS(ix, -ixx(2) + 1, iz) 547 end if 548 if ((gz < M_ZERO) .and. (gy >= M_ZERO)) then 549 fft_Coulb_FS(ix, iy, iz) = fft_Coulb_FS(ix, iy, -ixx(3) + 1) 550 end if 551 if ((gz < M_ZERO) .and. (gy < M_ZERO) ) then 552 fft_Coulb_FS(ix, iy, iz) = fft_Coulb_FS(ix, -ixx(2) + 1, -ixx(3) + 1) 553 end if 554 end if 555 556 else 557 if (mesh%sb%periodic_dim == 1) then 558 fft_Coulb_FS(ix, iy, iz) = -(M_HALF*log(r_c) - M_FOURTH)*r_c**2 559 else if (mesh%sb%periodic_dim == 0) then 560 fft_Coulb_FS(ix, iy, iz) = poisson_cutoff_3D_1D_finite(M_ZERO, M_ZERO, & 561 M_TWO*mesh%sb%xsize, M_TWO*mesh%sb%rsize) 562 end if 563 564 end if 565 end do 566 end do 567 568 if( mesh%sb%periodic_dim == 0 ) call spline_end(cylinder_cutoff_f) 569 end do 570 571 do iz = 1, cube%fs_n_global(3) 572 do iy = 1, cube%fs_n_global(2) 573 do ix = 1, cube%fs_n_global(1) 574 fft_Coulb_FS(ix, iy, iz) = M_FOUR*M_PI*fft_Coulb_FS(ix, iy, iz) 575 end do 576 end do 577 end do 578 579 call dfourier_space_op_init(coulb, cube, fft_Coulb_FS) 580 581 SAFE_DEALLOCATE_A(fft_Coulb_FS) 582 SAFE_DEALLOCATE_A(x) 583 SAFE_DEALLOCATE_A(y) 584 POP_SUB(poisson_fft_build_3d_1d) 585 end subroutine poisson_fft_build_3d_1d 586 !----------------------------------------------------------------- 587 588 589 !----------------------------------------------------------------- 590 !> C. A. Rozzi et al., Phys. Rev. B 73, 205119 (2006), Table I 591 subroutine poisson_fft_build_3d_0d(namespace, mesh, cube, kernel, coulb) 592 type(namespace_t), intent(in) :: namespace 593 type(mesh_t), intent(in) :: mesh 594 type(cube_t), intent(in) :: cube 595 integer, intent(in) :: kernel 596 type(fourier_space_op_t), intent(inout) :: coulb 597 598 integer :: ix, iy, iz, ixx(3), db(3), lx, ly, lz, n1, n2, n3 599 FLOAT :: temp(3), modg2 600 FLOAT :: r_c, gg(3), default_r_c 601 FLOAT, allocatable :: fft_coulb_FS(:,:,:) 602 603 PUSH_SUB(poisson_fft_build_3d_0d) 604 605 db(1:3) = cube%rs_n_global(1:3) 606 607 if (kernel /= POISSON_FFT_KERNEL_CORRECTED) then 608 default_r_c = maxval(db(1:3)*mesh%spacing(1:3)/M_TWO) 609 call get_cutoff(namespace, default_r_c, r_c) 610 end if 611 612 n1 = max(1, cube%fs_n(1)) 613 n2 = max(1, cube%fs_n(2)) 614 n3 = max(1, cube%fs_n(3)) 615 616 ! store the fourier transform of the Coulomb interaction 617 ! store only the relevant part if PFFT is used 618 SAFE_ALLOCATE(fft_Coulb_FS(1:n1,1:n2,1:n3)) 619 fft_Coulb_FS = M_ZERO 620 621 temp(1:3) = M_TWO*M_PI/(db(1:3)*mesh%spacing(1:3)) 622 do lx = 1, n1 623 ix = cube%fs_istart(1) + lx - 1 624 ixx(1) = pad_feq(ix, db(1), .true.) 625 do ly = 1, n2 626 iy = cube%fs_istart(2) + ly - 1 627 ixx(2) = pad_feq(iy, db(2), .true.) 628 do lz = 1, n3 629 iz = cube%fs_istart(3) + lz - 1 630 ixx(3) = pad_feq(iz, db(3), .true.) 631 632 call poisson_fft_gg_transform(ixx, temp, mesh%sb, coulb%qq, gg, modg2) 633 634 !HH 635 if(cube%fft%library.eq.FFTLIB_NFFT) then 636 modg2=cube%Lfs(ix,1)**2+cube%Lfs(iy,2)**2+cube%Lfs(iz,3)**2 637 r_c = default_r_c*M_TWO 638 end if 639 640 if(abs(modg2) > M_EPSILON) then 641 select case(kernel) 642 case(POISSON_FFT_KERNEL_SPH) 643 fft_Coulb_FS(lx, ly, lz) = poisson_cutoff_3D_0D(sqrt(modg2),r_c)/modg2 644 case(POISSON_FFT_KERNEL_CORRECTED) 645 fft_Coulb_FS(lx, ly, lz) = M_ONE/modg2 646 end select 647 else 648 select case(kernel) 649 case(POISSON_FFT_KERNEL_SPH) 650 fft_Coulb_FS(lx, ly, lz) = r_c**2/M_TWO 651 case (POISSON_FFT_KERNEL_CORRECTED) 652 fft_Coulb_FS(lx, ly, lz) = M_ZERO 653 end select 654 end if 655 end do 656 end do 657 end do 658 659 do iz = 1, cube%fs_n(3) 660 do iy = 1, cube%fs_n(2) 661 do ix = 1, cube%fs_n(1) 662 fft_Coulb_FS(ix, iy, iz) = M_FOUR*M_PI*fft_Coulb_FS(ix, iy, iz) 663 end do 664 end do 665 end do 666 667 call dfourier_space_op_init(coulb, cube, fft_coulb_fs, in_device = (kernel /= POISSON_FFT_KERNEL_CORRECTED)) 668 669 SAFE_DEALLOCATE_A(fft_Coulb_FS) 670 POP_SUB(poisson_fft_build_3d_0d) 671 end subroutine poisson_fft_build_3d_0d 672 !----------------------------------------------------------------- 673 674 675 !----------------------------------------------------------------- 676 !> A. Castro et al., Phys. Rev. B 80, 033102 (2009) 677 subroutine poisson_fft_build_2d_0d(namespace, mesh, cube, coulb) 678 type(namespace_t), intent(in) :: namespace 679 type(mesh_t), intent(in) :: mesh 680 type(cube_t), intent(in) :: cube 681 type(fourier_space_op_t), intent(inout) :: coulb 682 683 type(spline_t) :: besselintf 684 integer :: i, ix, iy, ixx(2), db(2), npoints 685 FLOAT :: temp(2), vec, r_c, maxf, dk, default_r_c 686 FLOAT, allocatable :: x(:), y(:) 687 FLOAT, allocatable :: fft_coulb_FS(:,:,:) 688 689 PUSH_SUB(poisson_fft_build_2d_0d) 690 691 db(1:2) = cube%rs_n_global(1:2) 692 693 default_r_c = maxval(db(1:2)*mesh%spacing(1:2)/M_TWO) 694 call get_cutoff(namespace, default_r_c, r_c) 695 696 call spline_init(besselintf) 697 698 ! store the fourier transform of the Coulomb interaction 699 SAFE_ALLOCATE(fft_Coulb_FS(1:cube%fs_n_global(1), 1:cube%fs_n_global(2), 1:cube%fs_n_global(3))) 700 fft_Coulb_FS = M_ZERO 701 temp(1:2) = M_TWO*M_PI/(db(1:2)*mesh%spacing(1:2)) 702 703 maxf = r_c * sqrt((temp(1)*db(1)/2)**2 + (temp(2)*db(2)/2)**2) 704 dk = CNST(0.25) ! This seems to be reasonable. 705 npoints = nint(maxf/dk) 706 SAFE_ALLOCATE(x(1:npoints)) 707 SAFE_ALLOCATE(y(1:npoints)) 708 x(1) = M_ZERO 709 y(1) = M_ZERO 710 do i = 2, npoints 711 x(i) = (i-1) * maxf / (npoints-1) 712 y(i) = y(i-1) + poisson_cutoff_2D_0D(x(i-1), x(i)) 713 end do 714 call spline_fit(npoints, x, y, besselintf) 715 716 do iy = 1, cube%fs_n_global(2) 717 ixx(2) = pad_feq(iy, db(2), .true.) 718 do ix = 1, cube%fs_n_global(1) 719 ixx(1) = pad_feq(ix, db(1), .true.) 720 vec = sqrt( (temp(1)*ixx(1))**2 + (temp(2)*ixx(2))**2) 721 if (vec > M_ZERO) then 722 fft_coulb_fs(ix, iy, 1) = (M_TWO * M_PI / vec) * spline_eval(besselintf, vec*r_c) 723 else 724 fft_coulb_fs(ix, iy, 1) = M_TWO * M_PI * r_c 725 end if 726 end do 727 end do 728 729 call dfourier_space_op_init(coulb, cube, fft_Coulb_FS) 730 731 SAFE_DEALLOCATE_A(fft_Coulb_FS) 732 SAFE_DEALLOCATE_A(x) 733 SAFE_DEALLOCATE_A(y) 734 call spline_end(besselintf) 735 POP_SUB(poisson_fft_build_2d_0d) 736 end subroutine poisson_fft_build_2d_0d 737 !----------------------------------------------------------------- 738 739 740 !----------------------------------------------------------------- 741 !> A. Castro et al., Phys. Rev. B 80, 033102 (2009) 742 subroutine poisson_fft_build_2d_1d(namespace, mesh, cube, coulb) 743 type(namespace_t), intent(in) :: namespace 744 type(mesh_t), intent(in) :: mesh 745 type(cube_t), intent(in) :: cube 746 type(fourier_space_op_t), intent(inout) :: coulb 747 748 integer :: ix, iy, ixx(2), db(2) 749 FLOAT :: temp(2), r_c, gx, gy, default_r_c 750 FLOAT, allocatable :: fft_coulb_FS(:,:,:) 751 752 PUSH_SUB(poisson_fft_build_2d_1d) 753 754 db(1:2) = cube%rs_n_global(1:2) 755 756 default_r_c = db(2)*mesh%spacing(2)/M_TWO 757 call get_cutoff(namespace, default_r_c, r_c) 758 759 ! store the fourier transform of the Coulomb interaction 760 SAFE_ALLOCATE(fft_Coulb_FS(1:cube%fs_n_global(1), 1:cube%fs_n_global(2), 1:cube%fs_n_global(3))) 761 fft_Coulb_FS = M_ZERO 762 temp(1:2) = M_TWO*M_PI/(db(1:2)*mesh%spacing(1:2)) 763 764 ! First, the term ix = 0 => gx = 0. 765 fft_coulb_fs(1, 1, 1) = -M_FOUR * r_c * (log(r_c)-M_ONE) 766 do iy = 2, cube%fs_n_global(2) 767 ixx(2) = pad_feq(iy, db(2), .true.) 768 gy = temp(2)*ixx(2) 769 fft_coulb_fs(1, iy, 1) = -M_FOUR * poisson_cutoff_intcoslog(r_c, gy, M_ONE ) 770 end do 771 772 do ix = 2, cube%fs_n_global(1) 773 ixx(1) = pad_feq(ix, db(1), .true.) 774 gx = temp(1)*ixx(1) 775 do iy = 1, db(2) 776 ixx(2) = pad_feq(iy, db(2), .true.) 777 gy = temp(2)*ixx(2) 778 fft_coulb_fs(ix, iy, 1) = poisson_cutoff_2d_1d(gy, gx, r_c) 779 end do 780 end do 781 782 call dfourier_space_op_init(coulb, cube, fft_Coulb_FS) 783 784 SAFE_DEALLOCATE_A(fft_Coulb_FS) 785 786 POP_SUB(poisson_fft_build_2d_1d) 787 end subroutine poisson_fft_build_2d_1d 788 !----------------------------------------------------------------- 789 790 791 !----------------------------------------------------------------- 792 !> A. Castro et al., Phys. Rev. B 80, 033102 (2009) 793 subroutine poisson_fft_build_2d_2d(mesh, cube, coulb) 794 type(mesh_t), intent(in) :: mesh 795 type(cube_t), intent(in) :: cube 796 type(fourier_space_op_t), intent(inout) :: coulb 797 798 integer :: ix, iy, ixx(2), db(2) 799 FLOAT :: temp(2), vec 800 FLOAT, allocatable :: fft_coulb_FS(:,:,:) 801 802 PUSH_SUB(poisson_fft_build_2d_2d) 803 804 db(1:2) = cube%rs_n_global(1:2) 805 806 ! store the fourier transform of the Coulomb interaction 807 SAFE_ALLOCATE(fft_Coulb_FS(1:cube%fs_n_global(1), 1:cube%fs_n_global(2), 1:cube%fs_n_global(3))) 808 fft_Coulb_FS = M_ZERO 809 temp(1:2) = M_TWO*M_PI/(db(1:2)*mesh%spacing(1:2)) 810 811 do iy = 1, cube%fs_n_global(2) 812 ixx(2) = pad_feq(iy, db(2), .true.) 813 do ix = 1, cube%fs_n_global(1) 814 ixx(1) = pad_feq(ix, db(1), .true.) 815 vec = sqrt( (temp(1)*ixx(1))**2 + (temp(2)*ixx(2))**2) 816 if (vec > M_ZERO) fft_coulb_fs(ix, iy, 1) = M_TWO * M_PI / vec 817 end do 818 end do 819 820 call dfourier_space_op_init(coulb, cube, fft_Coulb_FS) 821 822 SAFE_DEALLOCATE_A(fft_Coulb_FS) 823 POP_SUB(poisson_fft_build_2d_2d) 824 end subroutine poisson_fft_build_2d_2d 825 !----------------------------------------------------------------- 826 827 828 !----------------------------------------------------------------- 829 subroutine poisson_fft_build_1d_1d(mesh, cube, coulb, poisson_soft_coulomb_param) 830 type(mesh_t), intent(in) :: mesh 831 type(cube_t), intent(in) :: cube 832 type(fourier_space_op_t), intent(inout) :: coulb 833 FLOAT, intent(in) :: poisson_soft_coulomb_param 834 835 integer :: ix, ixx 836 FLOAT :: g 837 FLOAT, allocatable :: fft_coulb_fs(:, :, :) 838 839 PUSH_SUB(poisson_fft_build_1d_1d) 840 841 SAFE_ALLOCATE(fft_coulb_fs(1:cube%fs_n_global(1), 1:cube%fs_n_global(2), 1:cube%fs_n_global(3))) 842 fft_coulb_fs = M_ZERO 843 844 ! Fourier transform of Soft Coulomb interaction. 845 do ix = 1, cube%fs_n_global(1) 846 ixx = pad_feq(ix, cube%rs_n_global(1), .true.) 847 g = (ixx + coulb%qq(1))*M_PI/mesh%sb%lsize(1) 848 if(abs(g) > CNST(1e-6)) then 849 fft_coulb_fs(ix, 1, 1) = M_TWO * loct_bessel_k0(poisson_soft_coulomb_param*abs(g)) 850 end if 851 end do 852 853 call dfourier_space_op_init(coulb, cube, fft_coulb_fs) 854 SAFE_DEALLOCATE_A(fft_coulb_fs) 855 856 POP_SUB(poisson_fft_build_1d_1d) 857 end subroutine poisson_fft_build_1d_1d 858 !----------------------------------------------------------------- 859 860 861 !----------------------------------------------------------------- 862 subroutine poisson_fft_build_1d_0d(namespace, mesh, cube, coulb, poisson_soft_coulomb_param) 863 type(namespace_t), intent(in) :: namespace 864 type(mesh_t), intent(in) :: mesh 865 type(cube_t), intent(in) :: cube 866 type(fourier_space_op_t), intent(inout) :: coulb 867 FLOAT, intent(in) :: poisson_soft_coulomb_param 868 869 integer :: box(1), ixx(1), ix 870 FLOAT :: temp(1), g, r_c, default_r_c 871 FLOAT, allocatable :: fft_coulb_fs(:, :, :) 872 873 PUSH_SUB(poisson_fft_build_1d_0d) 874 875 box(1:1) = cube%rs_n_global(1:1) 876 877 default_r_c = box(1)*mesh%spacing(1)/M_TWO 878 call get_cutoff(namespace, default_r_c, r_c) 879 880 SAFE_ALLOCATE(fft_coulb_fs(1:cube%fs_n_global(1), 1:cube%fs_n_global(2), 1:cube%fs_n_global(3))) 881 fft_coulb_fs = M_ZERO 882 temp(1:1) = M_TWO*M_PI/(box(1:1)*mesh%spacing(1:1)) 883 884 ! Fourier transform of Soft Coulomb interaction. 885 do ix = 1, cube%fs_n_global(1) 886 ixx(1) = pad_feq(ix, box(1), .true.) 887 g = temp(1)*ixx(1) 888 fft_coulb_fs(ix, 1, 1) = poisson_cutoff_1D_0D(g, poisson_soft_coulomb_param, r_c) 889 end do 890 891 call dfourier_space_op_init(coulb, cube, fft_coulb_fs) 892 SAFE_DEALLOCATE_A(fft_coulb_fs) 893 894 POP_SUB(poisson_fft_build_1d_0d) 895 end subroutine poisson_fft_build_1d_0d 896 !----------------------------------------------------------------- 897 898 899 !----------------------------------------------------------------- 900 subroutine poisson_fft_end(this) 901 type(poisson_fft_t), intent(inout) :: this 902 903 PUSH_SUB(poisson_fft_end) 904 905 call fourier_space_op_end(this%coulb) 906 907 POP_SUB(poisson_fft_end) 908 end subroutine poisson_fft_end 909 910#include "undef.F90" 911#include "real.F90" 912#include "poisson_fft_inc.F90" 913#include "undef.F90" 914#include "complex.F90" 915#include "poisson_fft_inc.F90" 916 917 918end module poisson_fft_oct_m 919 920!! Local Variables: 921!! mode: f90 922!! coding: utf-8 923!! End: 924