TABLE OF CONTENTS
- ABINIT/m_gwls_polarisability
- m_hamiltonian/epsilon_k
- m_hamiltonian/matrix_function_epsilon_k
- m_hamiltonian/Pk
- m_hamiltonian/set_dielectric_function_frequency
ABINIT/m_gwls_polarisability [ Modules ]
NAME
m_gwls_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_polarisability 24 ! local modules 25 use m_gwls_utility 26 use m_gwls_wf 27 use m_gwls_valenceWavefunctions 28 use m_gwls_hamiltonian 29 use m_gwls_lineqsolver 30 31 ! abinit modules 32 use defs_basis 33 use m_errors 34 use m_abicore 35 use m_bandfft_kpt 36 37 use m_time, only : timab 38 39 implicit none 40 save 41 private
m_hamiltonian/epsilon_k [ Functions ]
[ Top ] [ m_hamiltonian ] [ Functions ]
NAME
epsilon_k
FUNCTION
.
INPUTS
OUTPUT
SOURCE
656 subroutine epsilon_k(psi_out,psi_in,omega) 657 658 implicit none 659 660 real(dp), intent(out) :: psi_out(2,npw_k) 661 real(dp), intent(in) :: psi_in(2,npw_k), omega(2) 662 663 ! ************************************************************************* 664 665 psi_out = psi_in 666 call sqrt_vc_k(psi_out) 667 call Pk(psi_out,omega) 668 call sqrt_vc_k(psi_out) 669 670 psi_out = psi_in - psi_out 671 end subroutine epsilon_k
m_hamiltonian/matrix_function_epsilon_k [ Functions ]
[ Top ] [ m_hamiltonian ] [ Functions ]
NAME
matrix_function_epsilon_k
FUNCTION
.
INPUTS
OUTPUT
SOURCE
714 subroutine matrix_function_epsilon_k(vector_out,vector_in,Hsize) 715 !---------------------------------------------------------------------------------------------------- 716 ! This function is a simple wrapper around epsilon_k to be fed to the Lanczos 717 ! algorithm. 718 !---------------------------------------------------------------------------------------------------- 719 implicit none 720 integer, intent(in) :: Hsize 721 complex(dpc), intent(out) :: vector_out(Hsize) 722 complex(dpc), intent(in) :: vector_in(Hsize) 723 724 725 ! local variables 726 real(dp) :: psik (2,npw_k) 727 real(dp) :: psik2(2,npw_k) 728 729 ! ************************************************************************* 730 731 ! convert from one format to the other 732 psik(1,:) = dble (vector_in(:)) 733 psik(2,:) = dimag(vector_in(:)) 734 735 ! act on vector 736 737 738 call epsilon_k(psik2 ,psik, matrix_function_omega) 739 740 ! convert back 741 vector_out = cmplx_1*psik2(1,:)+cmplx_i*psik2(2,:) 742 743 end subroutine matrix_function_epsilon_k
m_hamiltonian/Pk [ Functions ]
[ Top ] [ m_hamiltonian ] [ Functions ]
NAME
Pk
FUNCTION
.
INPUTS
OUTPUT
SOURCE
81 subroutine Pk(psi_inout,omega) 82 !=============================================================================== 83 ! 84 ! This routine applies the polarizability operator to an arbitrary state psi_inout. 85 ! 86 ! 87 ! A note about parallelism: 88 ! ------------------------- 89 ! The input/output is in "linear algebra" configuration, which is to say 90 ! that ALL processors have a fraction of the G-vectors for this state. Internally, 91 ! this routine will parallelise over bands and FFT, thus using the "FFT" 92 ! configuration of the data. For more information, see the lobpcgwf.F90, and 93 ! the article by F. Bottin et al. 94 !=============================================================================== 95 implicit none 96 97 real(dp), intent(inout) :: psi_inout(2,npw_k) 98 real(dp), intent(in) :: omega(2) 99 100 logical :: omega_imaginary 101 102 integer :: v, mb, iblk 103 104 real(dp):: norm_omega 105 real(dp), allocatable :: psik(:,:) 106 real(dp), allocatable :: psik_alltoall(:,:), psik_wrk_alltoall(:,:) 107 real(dp), allocatable :: psik_in_alltoall(:,:), psik_tmp_alltoall(:,:) 108 109 real(dp), allocatable :: psik_ext(:,:), psik_ext_alltoall(:,:) 110 111 112 real(dp), allocatable :: psir(:,:,:,:), psir_ext(:,:,:,:) 113 114 integer :: cplex 115 integer :: recy_i 116 117 118 integer :: mpi_band_rank 119 120 real(dp) :: list_SQMR_frequencies(2) 121 real(dp) :: list_QMR_frequencies(2,2) 122 123 124 integer, save :: icounter = 0 125 real(dp), save :: total_time1 = zero, total_time2 = zero, total_time = zero 126 127 128 integer :: num_op_v, i_op_v, case_op_v 129 real(dp):: factor_op_v 130 real(dp):: lbda 131 real(dp):: zz(2) 132 133 character(len=500) :: message 134 135 real(dp) :: tsec(2) 136 integer :: GWLS_TIMAB, OPTION_TIMAB 137 138 ! ************************************************************************* 139 140 GWLS_TIMAB = 1524 141 OPTION_TIMAB = 1 142 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 143 144 145 146 icounter = icounter + 1 147 call cpu_time(total_time1) 148 149 150 !======================================== 151 ! Allocate work arrays and define 152 ! important parameters 153 !======================================== 154 GWLS_TIMAB = 1525 155 OPTION_TIMAB = 1 156 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 157 158 159 cplex = 2 ! complex potential 160 161 mpi_band_rank = mpi_enreg%me_band 162 163 ABI_MALLOC(psik, (2,npw_kb)) 164 ABI_MALLOC(psik_alltoall, (2,npw_g)) 165 ABI_MALLOC(psik_wrk_alltoall, (2,npw_g)) 166 ABI_MALLOC(psik_tmp_alltoall, (2,npw_g)) 167 ABI_MALLOC(psik_in_alltoall, (2,npw_g)) 168 169 ABI_MALLOC(psir, (2,n4,n5,n6)) 170 ABI_MALLOC(psir_ext,(2,n4,n5,n6)) 171 172 173 OPTION_TIMAB = 2 174 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 175 176 177 !-------------------------------------------------------------------------------- 178 ! 179 ! The polarizability acting on a state | PSI > is given by 180 ! 181 ! Pk | PSI > = 2 sum_{v} | h_v^* phi_v > ( factor of 2 comes from sum on spin) 182 ! 183 ! | h_v > = Pc . OPERATOR_v . Pc | PSI^* phi_v > 184 ! 185 ! There are multiple cases to consider: 186 ! 187 ! CASE: 188 ! 1) |omega| = 0 189 ! OPERATOR_v = 1/[H-ev] 190 ! prefactor = -4 191 ! 192 ! 2) omega = i lbda (imaginary) 193 ! OPERATOR_v = [H-ev]/(lbda^2+[H-ev]^2) 194 ! prefactor = -4 195 ! 196 ! 3) omega = lbda (real) (USE SQMR) 197 ! OPERATOR_v = { 1/[H-ev+lbda]+ 1/[H-ev-lbda] } 198 ! prefactor = -2 199 ! 200 ! 4) omega = lbda (real) (USE QMR) 201 ! OPERATOR_v = { 1/[H-ev+lbda]+ 1/[H-ev-lbda] } 202 ! prefactor = -2 203 ! 204 ! It simplifies the code below to systematize the algorithm. 205 ! 206 !-------------------------------------------------------------------------------- 207 208 ! Check which part of omega is non-zero. Default is that omega is real. 209 if (abs(omega(2)) < 1.0d-12) then 210 omega_imaginary=.false. 211 norm_omega = omega(1) 212 213 elseif (abs(omega(1)) < 1.0d-12 .and. abs(omega(2)) > 1.0d-12) then 214 omega_imaginary = .true. 215 norm_omega = omega(2) 216 217 else 218 write(message,"(a,es16.8,3a)")& 219 "omega=",omega,",",ch10,& 220 "but either it's real or imaginary part need to be 0 for the polarisability routine to work." 221 ABI_ERROR(message) 222 end if 223 224 225 !----------------------------------------------------------------- 226 ! I) Prepare global values depending on the CASE being considered 227 !----------------------------------------------------------------- 228 229 if (norm_omega < 1.0D-12 ) then 230 case_op_v = 1 231 num_op_v = 1 232 factor_op_v = -4.0_dp 233 234 else if (omega_imaginary) then 235 case_op_v = 2 236 num_op_v = 1 237 factor_op_v = -4.0_dp 238 239 else if( .not. activate_inf_shift_poles) then 240 case_op_v = 3 241 num_op_v = 2 242 factor_op_v = -2.0_dp 243 244 else 245 case_op_v = 4 246 num_op_v = 2 247 factor_op_v = -2.0_dp 248 249 inf_shift_poles = dtset%zcut 250 251 252 write(message,*) " inf_shift_poles = ",inf_shift_poles 253 call wrtout(std_out,message,'COLL') 254 end if 255 256 257 write(message,10)" " 258 call wrtout(std_out,message,'COLL') 259 write(message,10) " Pk: applying the polarizability on states" 260 call wrtout(std_out,message,'COLL') 261 write(message,10) " ==============================================================" 262 call wrtout(std_out,message,'COLL') 263 write(message,12) " CASE : ",case_op_v 264 call wrtout(std_out,message,'COLL') 265 266 267 268 !----------------------------------------------------------------- 269 ! II) copy conjugate of initial wavefunction in local array, 270 ! and set inout array to zero. Each FFT row of processors 271 ! must have a copy of the initial wavefunction in FFT 272 ! configuration! 273 !----------------------------------------------------------------- 274 call cpu_time(time1) 275 276 GWLS_TIMAB = 1526 277 OPTION_TIMAB = 1 278 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 279 280 281 ABI_MALLOC(psik_ext,(2,npw_kb)) 282 ABI_MALLOC(psik_ext_alltoall,(2,npw_g)) 283 284 ! fill the array psik_ext with copies of the external state 285 do mb = 1, blocksize 286 psik_ext(:,(mb-1)*npw_k+1:mb*npw_k) = psi_inout(:,:) 287 end do 288 289 ! change configuration of the data, from LA to FFT 290 call wf_block_distribute(psik_ext, psik_ext_alltoall,1) ! LA -> FFT 291 ! Now every row of FFT processors has a copy of the external state. 292 293 ! Copy the external state to the real space format, appropriate for real space products to be 294 ! used later. 295 296 call g_to_r(psir_ext,psik_ext_alltoall) 297 psir_ext(2,:,:,:) = -psir_ext(2,:,:,:) 298 299 300 ! Don't need these arrays anymore... 301 ABI_FREE(psik_ext) 302 ABI_FREE(psik_ext_alltoall) 303 304 OPTION_TIMAB = 2 305 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 306 307 call cpu_time(time2) 308 time_fft = time_fft + time2-time1 309 counter_fft = counter_fft+1 310 311 ! set external state to zero, ready to start cumulating the answer 312 psi_inout = zero 313 314 !----------------------------------------------------------------- 315 ! III) Iterate on all valence bands 316 !----------------------------------------------------------------- 317 318 ! Loop on all blocks of eigenstates 319 do iblk = 1, nbdblock 320 321 322 ! What is the valence band index for this block and this row of FFT processors? It is not clear from the 323 ! code in lobpcgwf.F90; I'm going to *guess*. 324 325 v = (iblk-1)*blocksize + mpi_band_rank + 1 ! CAREFUL! This is a guess. Revisit this if code doesn't work as intended. 326 327 328 !Solving of Sternheiner equation 329 write(message,12) " band :", v 330 call wrtout(std_out,message,'COLL') 331 write(message,14) " eigenvalue (Ha) : ",eig(v) 332 call wrtout(std_out,message,'COLL') 333 write(message,14) " Re[omega] (Ha) : ",omega(1) 334 call wrtout(std_out,message,'COLL') 335 write(message,14) " Im[omega] (Ha) : ",omega(2) 336 call wrtout(std_out,message,'COLL') 337 338 !----------------------------------------------------------------- 339 ! IV) prepare some arrays, if they are needed 340 !----------------------------------------------------------------- 341 if (case_op_v == 3) then 342 343 list_SQMR_frequencies(1) = eig(v) - norm_omega 344 list_SQMR_frequencies(2) = eig(v) + norm_omega 345 346 else if (case_op_v == 4) then 347 list_QMR_frequencies(:,1) = (/eig(v)-norm_omega,-inf_shift_poles/) 348 list_QMR_frequencies(:,2) = (/eig(v)+norm_omega, inf_shift_poles/) 349 end if 350 351 352 !----------------------------------------------------------------- 353 ! V) Compute the real-space product of the input wavefunction 354 ! with the valence wavefunction. 355 !----------------------------------------------------------------- 356 ! Unfortunately, the input wavefunctions will be FFT-transformed k-> r nbandv times 357 ! because fourwf cannot multiply a real-space potential with a real-space wavefunction! 358 call cpu_time(time1) 359 360 GWLS_TIMAB = 1527 361 OPTION_TIMAB = 1 362 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 363 364 call gr_to_g(psik_in_alltoall,psir_ext,valence_wavefunctions_FFT(:,:,iblk)) 365 366 OPTION_TIMAB = 2 367 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 368 369 370 call cpu_time(time2) 371 time_fft = time_fft + time2-time1 372 counter_fft = counter_fft+1 373 374 !----------------------------------------------------------------- 375 ! VI) Project out to conduction space 376 !----------------------------------------------------------------- 377 call cpu_time(time1) 378 379 GWLS_TIMAB = 1528 380 OPTION_TIMAB = 1 381 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 382 383 !call pc_k(psik_in) 384 call pc_k_valence_kernel(psik_in_alltoall) 385 386 OPTION_TIMAB = 2 387 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 388 389 call cpu_time(time2) 390 time_proj = time_proj + time2-time1 391 counter_proj= counter_proj+1 392 393 394 !----------------------------------------------------------------- 395 ! VII) Loop on potential valence operators, 396 !----------------------------------------------------------------- 397 do i_op_v = 1, num_op_v 398 399 if (case_op_v == 1) then 400 401 ! frequency is zero, operator to apply is 1/[H-ev] 402 call cpu_time(time1) 403 GWLS_TIMAB = 1529 404 OPTION_TIMAB = 1 405 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 406 407 call sqmr(psik_in_alltoall,psik_tmp_alltoall,eig(v),1) 408 409 OPTION_TIMAB = 2 410 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 411 412 413 call cpu_time(time2) 414 time_sqmr = time_sqmr+ time2-time1 415 counter_sqmr= counter_sqmr+1 416 417 ! If we are constructing the $\hat \epsilon(i\omega = 0)$ matrix (and the Lanczos basis at the same time), 418 ! keep the Sternheimer solutions for use in the projected Sternheimer section (in LA configuration). 419 if(write_solution .and. ((iblk-1)*blocksize < nbandv)) then 420 call wf_block_distribute(psik, psik_tmp_alltoall, 2) ! FFT -> LA 421 do mb=1,blocksize 422 v = (iblk-1)*blocksize+mb 423 if(v <= nbandv) then 424 if(dtset%gwls_recycle == 1) then 425 Sternheimer_solutions_zero(:,:,index_solution,v) = psik(:,(mb-1)*npw_k+1:mb*npw_k) 426 end if 427 if(dtset%gwls_recycle == 2) then 428 recy_i = (index_solution-1)*nbandv + v 429 !BUG : On petrus, NAG 5.3.1 + OpenMPI 1.6.2 cause read(...,rec=i) to read the data written by write(...,rec=i+1). 430 !Workaround compatible only with nag : write(recy_unit,rec=recy_i+1). 431 write(recy_unit,rec=recy_i) psik(:,(mb-1)*npw_k+1:mb*npw_k) 432 end if 433 end if 434 end do 435 end if 436 437 else if (case_op_v == 2) then 438 439 ! frequency purely imaginary, operator to apply is 440 ! [H-ev]/(lbda^2+[H-ev]^2) 441 call cpu_time(time1) 442 443 GWLS_TIMAB = 1533 444 OPTION_TIMAB = 1 445 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 446 447 psik_wrk_alltoall = psik_in_alltoall 448 call Hpsik(psik_wrk_alltoall,eig(v)) 449 450 OPTION_TIMAB = 2 451 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 452 453 454 call cpu_time(time2) 455 time_H = time_H + time2-time1 456 counter_H = counter_H+1 457 458 call cpu_time(time1) 459 460 GWLS_TIMAB = 1530 461 OPTION_TIMAB = 1 462 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 463 464 call sqmr(psik_wrk_alltoall,psik_tmp_alltoall,eig(v),0,norm_omega,omega_imaginary) 465 466 OPTION_TIMAB = 2 467 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 468 469 470 call cpu_time(time2) 471 time_sqmr = time_sqmr+ time2-time1 472 counter_sqmr= counter_sqmr+1 473 474 else if (case_op_v == 3) then 475 476 ! frequency purely real, operator to apply is 477 ! 1/[H-ev +/- lbda] 478 ! TREATED WITH SQMR 479 480 481 lbda = list_SQMR_frequencies(i_op_v) 482 call cpu_time(time1) 483 484 GWLS_TIMAB = 1531 485 OPTION_TIMAB = 1 486 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 487 488 call sqmr(psik_in_alltoall,psik_tmp_alltoall,lbda,1) 489 490 OPTION_TIMAB = 2 491 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 492 493 494 495 call cpu_time(time2) 496 time_sqmr = time_sqmr+ time2-time1 497 counter_sqmr= counter_sqmr+1 498 499 else if (case_op_v == 4) then 500 ! frequency purely real, operator to apply is 501 ! 1/[H-ev +/- lbda] 502 ! TREATED WITH QMR 503 504 zz(:) = list_QMR_frequencies(:,i_op_v) 505 506 call cpu_time(time1) 507 GWLS_TIMAB = 1532 508 OPTION_TIMAB = 1 509 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 510 511 call qmr(psik_in_alltoall,psik_tmp_alltoall,zz) !,0) 512 513 OPTION_TIMAB = 2 514 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 515 516 517 call cpu_time(time2) 518 time_sqmr = time_sqmr+ time2-time1 519 counter_sqmr= counter_sqmr+1 520 521 end if 522 !----------------------------------------------------------------- 523 ! VIII) Project on conduction states 524 !----------------------------------------------------------------- 525 call cpu_time(time1) 526 GWLS_TIMAB = 1528 527 OPTION_TIMAB = 1 528 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 529 530 call pc_k_valence_kernel(psik_tmp_alltoall) 531 532 OPTION_TIMAB = 2 533 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 534 535 536 call cpu_time(time2) 537 time_proj = time_proj + time2-time1 538 counter_proj= counter_proj+1 539 540 !----------------------------------------------------------------- 541 ! IX) Conjugate result, and express in denpot format 542 !----------------------------------------------------------------- 543 544 call cpu_time(time1) 545 546 GWLS_TIMAB = 1526 547 OPTION_TIMAB = 1 548 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 549 550 551 call g_to_r(psir,psik_tmp_alltoall) 552 psir(2,:,:,:) = -psir(2,:,:,:) 553 554 555 OPTION_TIMAB = 2 556 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 557 558 call cpu_time(time2) 559 time_fft = time_fft + time2-time1 560 counter_fft = counter_fft+1 561 562 !----------------------------------------------------------------- 563 ! X) Multiply by valence state in real space 564 !----------------------------------------------------------------- 565 call cpu_time(time1) 566 GWLS_TIMAB = 1527 567 OPTION_TIMAB = 1 568 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 569 570 call gr_to_g(psik_alltoall,psir,valence_wavefunctions_FFT(:,:,iblk)) 571 572 573 OPTION_TIMAB = 2 574 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 575 576 577 call cpu_time(time2) 578 time_fft = time_fft + time2-time1 579 counter_fft = counter_fft+1 580 581 !----------------------------------------------------------------- 582 ! XI) Return to LA configuration, and cumulate the sum 583 !----------------------------------------------------------------- 584 call wf_block_distribute(psik, psik_alltoall,2) ! FFT -> LA 585 586 do mb = 1, blocksize 587 588 v = (iblk-1)*blocksize + mb 589 590 if (v <= nbandv) then 591 ! only add contributions from valence 592 psi_inout(:,:) = psi_inout(:,:) + psik(:,(mb-1)*npw_k+1:mb*npw_k) 593 end if 594 595 end do 596 597 598 end do ! i_op_v 599 600 end do ! iblk 601 602 603 !----------------------------------------------------------------- 604 ! XII) account for prefactor 605 !----------------------------------------------------------------- 606 psi_inout = factor_op_v*psi_inout 607 608 609 GWLS_TIMAB = 1525 610 OPTION_TIMAB = 1 611 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 612 613 614 ABI_FREE(psik) 615 ABI_FREE(psik_alltoall) 616 ABI_FREE(psik_wrk_alltoall) 617 ABI_FREE(psik_tmp_alltoall) 618 ABI_FREE(psik_in_alltoall) 619 620 ABI_FREE(psir) 621 ABI_FREE(psir_ext) 622 623 OPTION_TIMAB = 2 624 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 625 626 627 call cpu_time(total_time2) 628 629 total_time = total_time + total_time2 - total_time1 630 631 632 GWLS_TIMAB = 1524 633 OPTION_TIMAB = 2 634 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 635 636 10 format(A) 637 12 format(A,I5) 638 14 format(A,F24.16) 639 640 end subroutine Pk
m_hamiltonian/set_dielectric_function_frequency [ Functions ]
[ Top ] [ m_hamiltonian ] [ Functions ]
NAME
set_dielectric_function_frequency
FUNCTION
.
INPUTS
OUTPUT
SOURCE
687 subroutine set_dielectric_function_frequency(omega) 688 !---------------------------------------------------------------------------------------------------- 689 ! This routine sets the value of the module's frequency. 690 !---------------------------------------------------------------------------------------------------- 691 implicit none 692 real(dp), intent(in) :: omega(2) 693 694 ! ************************************************************************* 695 696 matrix_function_omega(:) = omega(:) 697 698 end subroutine set_dielectric_function_frequency