TABLE OF CONTENTS
- ABINIT/m_gwls_model_polarisability
- m_hamiltonian/cleanup_Pk_model
- m_hamiltonian/epsilon_k_model
- m_hamiltonian/matrix_function_epsilon_model_operator
- m_hamiltonian/Pk_model
- m_hamiltonian/Pk_model_implementation_1
- m_hamiltonian/setup_Pk_model
ABINIT/m_gwls_model_polarisability [ Modules ]
NAME
m_gwls_model_polarisability
FUNCTION
.
COPYRIGHT
Copyright (C) 2009-2024 ABINIT group (JLJ, BR, MC) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 23 module m_gwls_model_polarisability 24 25 ! local modules 26 use m_gwls_utility 27 use m_gwls_wf 28 use m_gwls_valenceWavefunctions 29 use m_gwls_hamiltonian 30 use m_gwls_lineqsolver 31 32 ! abinit modules 33 use defs_basis 34 use m_abicore 35 use m_bandfft_kpt 36 use m_errors 37 38 use m_time, only : timab 39 40 implicit none 41 save 42 private
m_hamiltonian/cleanup_Pk_model [ Functions ]
[ Top ] [ m_hamiltonian ] [ Functions ]
NAME
cleanup_Pk_model
FUNCTION
.
INPUTS
OUTPUT
SOURCE
261 subroutine cleanup_Pk_model() 262 263 implicit none 264 265 ! ************************************************************************* 266 267 if (allocated(model_Y)) then 268 ABI_FREE(model_Y) 269 end if 270 271 if (allocated(model_Y_LA)) then 272 ABI_FREE(model_Y_LA) 273 end if 274 275 276 277 if (allocated(sqrt_density)) then 278 ABI_FREE(sqrt_density) 279 end if 280 281 if ( allocated(psir_model)) then 282 ABI_FREE(psir_model) 283 end if 284 285 if (allocated(psir_ext_model)) then 286 ABI_FREE(psir_ext_model) 287 end if 288 289 end subroutine cleanup_Pk_model
m_hamiltonian/epsilon_k_model [ Functions ]
[ Top ] [ m_hamiltonian ] [ Functions ]
NAME
epsilon_k_model
FUNCTION
.
INPUTS
OUTPUT
SOURCE
85 subroutine epsilon_k_model(psi_out,psi_in) 86 87 implicit none 88 89 real(dp), intent(out) :: psi_out(2,npw_k) 90 real(dp), intent(in) :: psi_in(2,npw_k) 91 92 real(dp) :: psik(2,npw_k) 93 94 ! ************************************************************************* 95 96 psik = psi_in 97 call sqrt_vc_k(psik) 98 call Pk_model(psi_out ,psik) 99 call sqrt_vc_k(psi_out) 100 psi_out = psi_in - psi_out 101 102 end subroutine epsilon_k_model
m_hamiltonian/matrix_function_epsilon_model_operator [ Functions ]
[ Top ] [ m_hamiltonian ] [ Functions ]
NAME
matrix_function_epsilon_model_operator
FUNCTION
.
INPUTS
OUTPUT
SOURCE
593 subroutine matrix_function_epsilon_model_operator(vector_out,vector_in,Hsize) 594 !---------------------------------------------------------------------------------------------------- 595 ! This function returns the action of the operator epsilon_model on a given vector. 596 ! It is assumed that the frequency has been set in the module using setup_Pk_model. 597 ! 598 ! 599 !---------------------------------------------------------------------------------------------------- 600 implicit none 601 integer, intent(in) :: Hsize 602 complex(dpc), intent(out) :: vector_out(Hsize) 603 complex(dpc), intent(in) :: vector_in(Hsize) 604 605 ! local variables 606 real(dp) :: psik (2,Hsize) 607 real(dp) :: psik2(2,Hsize) 608 609 ! ************************************************************************* 610 611 ! convert from one format to the other 612 psik(1,:) = dble (vector_in(:)) 613 psik(2,:) = dimag(vector_in(:)) 614 615 call epsilon_k_model(psik2 ,psik) 616 617 ! Act with epsilon_model 618 vector_out = cmplx_1*psik2(1,:)+cmplx_i*psik2(2,:) 619 620 end subroutine matrix_function_epsilon_model_operator
m_hamiltonian/Pk_model [ Functions ]
[ Top ] [ m_hamiltonian ] [ Functions ]
NAME
Pk_model
FUNCTION
.
INPUTS
OUTPUT
SOURCE
305 subroutine Pk_model(psi_out,psi_in) 306 !------------------------------------------------------------------------------------------------------------------------ 307 ! Returns the action of a frequency-dependent model susceptibility 308 !------------------------------------------------------------------------------------------------------------------------ 309 implicit none 310 311 real(dp), intent(out) :: psi_out(2,npw_k) 312 real(dp), intent(in) :: psi_in(2,npw_k) 313 314 ! ************************************************************************* 315 316 if ( dielectric_model_type == 1) then 317 call Pk_model_implementation_1(psi_out ,psi_in) 318 else if ( dielectric_model_type == 2) then 319 ! call Pk_model_implementation_2(psi_out ,psi_in) 320 end if 321 322 end subroutine Pk_model
m_hamiltonian/Pk_model_implementation_1 [ Functions ]
[ Top ] [ m_hamiltonian ] [ Functions ]
NAME
Pk_model_implementation_1
FUNCTION
.
INPUTS
OUTPUT
SOURCE
338 subroutine Pk_model_implementation_1(psi_out,psi_in) 339 !------------------------------------------------------------------------------------------------------------------------ 340 ! Returns the action of a frequency-dependent model susceptibility 341 !------------------------------------------------------------------------------------------------------------------------ 342 implicit none 343 344 real(dp), intent(out) :: psi_out(2,npw_k) 345 real(dp), intent(in) :: psi_in(2,npw_k) 346 347 integer :: v 348 349 real(dp), allocatable :: psik(:,:), psik_g(:,:) 350 351 integer, save :: icounter = 0 352 353 integer :: mb, iblk 354 355 integer :: mpi_band_rank 356 357 real(dp) :: time1, time2 358 real(dp) :: total_time1, total_time2 359 360 real(dp), save :: fft_time = zero 361 real(dp), save :: projection_time = zero 362 real(dp), save :: Y_time = zero 363 real(dp), save :: total_time = zero 364 365 366 real(dp) :: tsec(2) 367 integer :: GWLS_TIMAB, OPTION_TIMAB 368 369 ! ************************************************************************* 370 371 GWLS_TIMAB = 1534 372 OPTION_TIMAB = 1 373 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 374 375 376 377 call cpu_time(total_time1) 378 icounter = icounter + 1 379 380 GWLS_TIMAB = 1535 381 OPTION_TIMAB = 1 382 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 383 384 385 ABI_MALLOC(psik, (2,npw_kb)) 386 ABI_MALLOC(psik_g, (2,npw_g)) 387 388 OPTION_TIMAB = 2 389 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 390 391 392 ! initialize the output to zero 393 psi_out = zero 394 395 ! MPI information 396 mpi_band_rank = mpi_enreg%me_band 397 398 399 !----------------------------------------------------------------- 400 ! Put a copy of the external state on every row of FFT processors. 401 ! 402 ! The, copy conjugate of initial wavefunction in local array, 403 ! and set inout array to zero. 404 !----------------------------------------------------------------- 405 call cpu_time(time1) 406 GWLS_TIMAB = 1536 407 OPTION_TIMAB = 1 408 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 409 410 ! fill the array psik_ext with copies of the external state 411 do mb = 1, blocksize 412 psik(:,(mb-1)*npw_k+1:mb*npw_k) = psi_in(:,:) 413 end do 414 415 ! change configuration of the data, from LA to FFT 416 call wf_block_distribute(psik, psik_g,1) ! LA -> FFT 417 ! Now every row of FFT processors has a copy of the external state. 418 419 420 ! Store external state in real-space format, inside module-defined work array psir_ext_model 421 ! Each row of FFT processors will have a copy! 422 call g_to_r(psir_ext_model,psik_g) 423 424 ! Conjugate the external wavefunction; result will be conjugated again later, 425 ! insuring we are in fact acting on psi, and not psi^*. This conjugation 426 ! is only done for algorithmic convenience. 427 psir_ext_model(2,:,:,:) = -psir_ext_model(2,:,:,:) 428 429 430 OPTION_TIMAB = 2 431 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 432 433 434 call cpu_time(time2) 435 fft_time = fft_time+time2-time1 436 437 438 ! Loop on all blocks of eigenstates 439 do iblk = 1, nbdblock 440 441 v = (iblk-1)*blocksize + mpi_band_rank + 1 ! CAREFUL! This is a guess. Revisit this if code doesn't work as intended. 442 443 call cpu_time(time1) 444 GWLS_TIMAB = 1537 445 OPTION_TIMAB = 1 446 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 447 448 449 ! Multiply valence state by external state, yielding |psik_g> = | phi_v x psi_in^* > 450 call gr_to_g(psik_g, psir_ext_model, valence_wavefunctions_FFT(:,:,iblk)) 451 452 OPTION_TIMAB = 2 453 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 454 455 456 call cpu_time(time2) 457 fft_time = fft_time+time2-time1 458 459 ! Project out to conduction space 460 call cpu_time(time1) 461 GWLS_TIMAB = 1538 462 OPTION_TIMAB = 1 463 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 464 465 call pc_k_valence_kernel(psik_g) 466 467 OPTION_TIMAB = 2 468 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 469 470 471 call cpu_time(time2) 472 projection_time = projection_time+time2-time1 473 474 475 ! act with model susceptibility 476 call cpu_time(time1) 477 GWLS_TIMAB = 1539 478 OPTION_TIMAB = 1 479 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 480 481 482 psik_g(1,:) = psik_g(1,:)*model_Y(:) 483 psik_g(2,:) = psik_g(2,:)*model_Y(:) 484 485 OPTION_TIMAB = 2 486 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 487 488 489 call cpu_time(time2) 490 Y_time = Y_time+time2-time1 491 492 ! Project out to conduction space, again! 493 call cpu_time(time1) 494 GWLS_TIMAB = 1538 495 OPTION_TIMAB = 1 496 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 497 498 call pc_k_valence_kernel(psik_g) 499 500 OPTION_TIMAB = 2 501 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 502 call cpu_time(time2) 503 504 projection_time= projection_time+time2-time1 505 506 ! Express result in real space, in module-defined work array psir_model 507 call cpu_time(time1) 508 GWLS_TIMAB = 1536 509 OPTION_TIMAB = 1 510 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 511 512 call g_to_r(psir_model,psik_g) 513 ! conjugate the result, cancelling the initial conjugation described earlier. 514 psir_model(2,:,:,:) = -psir_model(2,:,:,:) 515 516 OPTION_TIMAB = 2 517 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 518 519 call cpu_time(time2) 520 fft_time = fft_time+time2-time1 521 522 ! Multiply by valence state in real space 523 call cpu_time(time1) 524 GWLS_TIMAB = 1537 525 OPTION_TIMAB = 1 526 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 527 528 call gr_to_g(psik_g,psir_model, valence_wavefunctions_FFT(:,:,iblk)) 529 530 OPTION_TIMAB = 2 531 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 532 call cpu_time(time2) 533 fft_time = fft_time+time2-time1 534 535 536 537 ! Return to linear algebra format, and add condtribution 538 GWLS_TIMAB = 1540 539 OPTION_TIMAB = 1 540 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 541 542 call wf_block_distribute(psik, psik_g,2) ! FFT -> LA 543 544 do mb = 1, blocksize 545 546 v = (iblk-1)*blocksize + mb 547 if ( v <= nbandv) then 548 psi_out(:,:) = psi_out(:,:) + psik(:,(mb-1)*npw_k+1:mb*npw_k) 549 end if 550 end do 551 552 OPTION_TIMAB = 2 553 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 554 555 556 end do ! iblk 557 558 call cpu_time(total_time2) 559 total_time = total_time + total_time2-total_time1 560 561 GWLS_TIMAB = 1535 562 OPTION_TIMAB = 1 563 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 564 565 ABI_FREE(psik) 566 ABI_FREE(psik_g) 567 568 OPTION_TIMAB = 2 569 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 570 571 572 GWLS_TIMAB = 1534 573 OPTION_TIMAB = 2 574 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 575 576 577 end subroutine Pk_model_implementation_1
m_hamiltonian/setup_Pk_model [ Functions ]
[ Top ] [ m_hamiltonian ] [ Functions ]
NAME
setup_Pk_model
FUNCTION
.
INPUTS
OUTPUT
SOURCE
118 subroutine setup_Pk_model(omega,epsilon_0) 119 !--------------------------------------------------------------- 120 ! 121 ! This subroutine prepares a global array in order to 122 ! act with the model susceptibility, given by 123 ! 124 ! Pk_model(r,r',i omega) = sum_v phi_v(r) Y(r-r',i omega) phi^*_v(r') 125 ! 126 ! 127 ! This subroutine computes the Fourier transform of Y(r,i omega), 128 ! Y(G,omega), for a given omega and epsilon_0. 129 ! 130 ! It is assumed that omega and epsilon_0 >= 0. Also, the model 131 ! describes an IMAGINARY frequency i omega. 132 !--------------------------------------------------------------- 133 real(dp), intent(in) :: epsilon_0, omega 134 135 real(dp) :: theta, R_omega 136 real(dp) :: x, y 137 138 139 integer :: ig 140 141 real(dp),allocatable :: G_array(:) 142 143 ! ************************************************************************* 144 145 if (.not. allocated(psir_model)) then 146 ABI_MALLOC(psir_model, (2,n4,n5,n6)) 147 end if 148 149 if (.not. allocated(psir_ext_model)) then 150 ABI_MALLOC(psir_ext_model, (2,n4,n5,n6)) 151 end if 152 153 R_omega = 2.0_dp*sqrt(epsilon_0**2+omega**2) 154 if(abs(omega) > tol16 .or. abs(epsilon_0) > tol16) then 155 theta = atan2(omega,epsilon_0) 156 else 157 theta = zero 158 end if 159 160 x = sqrt(R_omega)*cos(0.5_dp*theta) 161 y = sqrt(R_omega)*sin(0.5_dp*theta) 162 163 164 if (.not. allocated(model_Y)) then 165 ABI_MALLOC(model_Y, (npw_g)) 166 end if 167 168 if (.not. allocated(model_Y_LA)) then 169 ABI_MALLOC(model_Y_LA, (npw_k)) 170 end if 171 172 !================================================================================ 173 ! Compute model_Y, in FFT configuration 174 !================================================================================ 175 176 ABI_MALLOC(G_array,(npw_g)) 177 G_array(:) = sqrt(2.0_dp*kinpw_gather(:)) 178 179 model_Y(:) = zero 180 do ig = 1, npw_g 181 182 183 if (G_array(ig) > tol12) then 184 ! G != 0. 185 model_Y(ig) = -4.0_dp/G_array(ig)* & 186 ((G_array(ig)+y)/((G_array(ig)+y)**2+x**2) & 187 + (G_array(ig)-y)/((G_array(ig)-y)**2+x**2)) 188 189 190 else 191 if ( abs(epsilon_0) < tol12 ) then 192 model_Y(ig) = zero 193 else 194 model_Y(ig) = -4.0_dp*epsilon_0/(epsilon_0**2+omega**2) 195 end if 196 end if 197 end do ! ig 198 199 ABI_FREE(G_array) 200 201 !================================================================================ 202 ! Compute model_Y_LA, in LA configuration 203 !================================================================================ 204 205 ABI_MALLOC(G_array,(npw_k)) 206 G_array(:) = sqrt(2.0_dp*kinpw(:)) 207 208 model_Y_LA(:) = zero 209 do ig = 1, npw_k 210 211 if (G_array(ig) > tol12) then 212 ! G != 0. 213 model_Y_LA(ig) = -4.0_dp/G_array(ig)* & 214 ((G_array(ig)+y)/((G_array(ig)+y)**2+x**2) & 215 + (G_array(ig)-y)/((G_array(ig)-y)**2+x**2)) 216 217 218 else 219 if ( abs(epsilon_0) < tol12 ) then 220 model_Y_LA(ig) = zero 221 else 222 model_Y_LA(ig) = -4.0_dp*epsilon_0/(epsilon_0**2+omega**2) 223 end if 224 end if 225 end do ! ig 226 227 ABI_FREE(G_array) 228 229 230 231 if (dielectric_model_type == 2) then 232 233 ABI_BUG('dielectric_model_type == 2 not properly implemented. Review code or input!') 234 235 !ABI_MALLOC(sqrt_density,(2,n4,n5,n6)) 236 !sqrt_density(:,:,:,:) = zero 237 !do v= 1, nbandv 238 ! sqrt_density(1,:,:,:) = sqrt_density(1,:,:,:) + valence_wfr(1,:,:,:,v)**2+valence_wfr(2,:,:,:,v)**2 239 !end do 240 !sqrt_density(1,:,:,:) = sqrt(sqrt_density(1,:,:,:)) 241 242 end if 243 244 245 end subroutine setup_Pk_model