TABLE OF CONTENTS
- ABINIT/m_rf2
- m_rf2/rf2_accumulate_bands
- m_rf2/rf2_apply_hamiltonian
- m_rf2/rf2_destroy
- m_rf2/rf2_getidir
- m_rf2/rf2_getidirs
- m_rf2/rf2_t
ABINIT/m_rf2 [ Modules ]
NAME
m_rf2
FUNCTION
This module defines structures and provides procedures used to compute the 2nd order Sternheimer equation.
COPYRIGHT
Copyright (C) 2015-2024 ABINIT group (LB,MT) 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
17 #if defined HAVE_CONFIG_H 18 #include "config.h" 19 #endif 20 21 #include "abi_common.h" 22 23 MODULE m_rf2 24 25 use defs_basis 26 use m_abicore 27 use m_errors 28 use m_hamiltonian 29 use m_cgtools 30 31 use defs_abitypes, only : MPI_type 32 use m_getgh1c, only : getgh1c 33 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_copy, pawcprj_free, pawcprj_output 34 use m_getghc, only : getghc 35 use m_nonlop, only : nonlop 36 use m_getgh2c, only : getgh2c 37 38 implicit none 39 40 private
m_rf2/rf2_accumulate_bands [ Functions ]
[ Top ] [ m_rf2 ] [ Functions ]
NAME
rf2_init
FUNCTION
Compute the scalar product < vi | v1j > and add it to rf2%lambda_mn If necessary, compute < vi | v2j > and add it to rf2%amn.
INPUTS
rf2 : rf2_t object containing all rf2 data choice : 4 possible values, see "select case (choice)" below gs_hamkq <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k+q mpi_enreg=information about MPI parallelization iband : band index for vi idir1 (used only if debug_mode/=0) : direction of the 1st perturbation idir2 (used only if debug_mode/=0) : direction of the 2nd perturbation ipert1 (used only if debug_mode/=0) : 1st perturbation ipert2 (used only if debug_mode/=0) : 2nd perturbation jband : band index for v1j and v2j debug_mode : if /=0 : all < vi | v1j > and < vi | v2j > are printed in std_out vi,v1j,v2j : input vectors
OUTPUT
NOTES
SOURCE
258 subroutine rf2_accumulate_bands(rf2,choice,gs_hamkq,mpi_enreg,iband,idir1,idir2,ipert1,ipert2,& 259 jband,debug_mode,vi,v1j,v2j,factor_in) 260 261 !Arguments --------------------------------------------- 262 !scalars 263 integer,intent(in) :: choice,iband,idir1,idir2,ipert1,ipert2,jband,debug_mode 264 type(rf2_t),intent(inout) :: rf2 265 type(gs_hamiltonian_type),intent(in) :: gs_hamkq 266 type(MPI_type),intent(in) :: mpi_enreg 267 real(dp), intent(in) :: factor_in 268 269 !arrays 270 real(dp),intent(in) :: vi(2,rf2%size_wf),v1j(2,rf2%size_wf),v2j(2,rf2%size_wf) 271 272 !Local variables --------------------------------------- 273 !scalars 274 integer :: nband_k,size_wf 275 real(dp) :: dotr,dot2r,doti,dot2i 276 character(len=500) :: msg 277 character(len=15) :: bra_i,ket_j,op1,op2 278 character(len=2) :: pert1,pert2 279 280 ! ************************************************************************* 281 282 DBG_ENTER("COLL") 283 284 nband_k = rf2%nband_k 285 size_wf = rf2%size_wf 286 287 call dotprod_g(dotr,doti,gs_hamkq%istwf_k,size_wf,2,vi,v1j,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 288 289 if(debug_mode/=0) then 290 if (ipert1 == gs_hamkq%natom+1) then 291 pert1 = "dk" 292 else 293 pert1 = "dE" 294 end if 295 if (ipert2 == gs_hamkq%natom+1) then 296 pert2 = "dk" 297 else 298 pert2 = "dE" 299 end if 300 select case (choice) 301 case(1) ! < u^(pert2) | H^(0) | u^(pert1) > 302 write(bra_i,'(2a,2(i1,a))') ' < du/',pert2,idir2,'(',iband,') | ' 303 write(ket_j,'(2a,2(i1,a))') ' | du/',pert1,idir1,'(',jband,') > ' 304 write(op1,'(a)') ' H^(0) ' 305 write(op2,'(a)') ' S^(0) ' 306 case(2) ! < u^(pert2) | H^(pert1) | u^(0) > 307 write(bra_i,'(2a,2(i1,a))') ' < du/',pert2,idir2,'(',iband,') | ' 308 write(ket_j,'(a,i1,a)') ' | u^(0) (',jband,') > ' 309 write(op1,'(2a,i1,a)') ' dH/',pert1,idir1,' ' 310 write(op2,'(2a,i1,a)') ' dS/',pert1,idir1,' ' 311 case(3) ! < u^(0) | H^(pert1pert2) | u^(0) > 312 write(bra_i,'(a,i1,a)') ' < u^(0) (',iband,') | ' 313 write(ket_j,'(a,i1,a)') ' | u^(0) (',jband,') > ' 314 write(op1,'(2(2a,i1),a)') 'd^2H/(',pert1,idir1,' ',pert2,idir2,')' 315 write(op2,'(2(2a,i1),a)') 'd^2S/(',pert1,idir1,' ',pert2,idir2,')' 316 case(4) ! < u^(0) | H^(pert2) | u^(pert1) > 317 write(bra_i,'(a,i1,a)') ' < u^(0) (',iband,') | ' 318 write(ket_j,'(2a,2(i1,a))') ' | du/',pert1,idir1,'(',jband,') > ' 319 write(op1,'(2a,i1,a)') ' dH/',pert2,idir2,' ' 320 write(op2,'(2a,i1,a)') ' dS/',pert2,idir2,' ' 321 end select 322 dot2r = dotr ; if (abs(dot2r)<tol9) dot2r = zero ! in order to hide the numerical noise 323 dot2i = doti ; if (abs(dot2i)<tol9) dot2i = zero ! in order to hide the numerical noise 324 write(msg,'(3a,2(a,es17.8E3))') bra_i,op1,ket_j,' = ',dot2r,',',dot2i 325 call wrtout(std_out,msg) 326 end if 327 328 rf2%lambda_mn(:,iband+(jband-1)*nband_k) = factor_in*(/dotr,doti/) + rf2%lambda_mn(:,iband+(jband-1)*nband_k) 329 330 if (choice == 1 .or. gs_hamkq%usepaw==1) then 331 call dotprod_g(dotr,doti,gs_hamkq%istwf_k,size_wf,2,vi,v2j,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 332 333 if(debug_mode/=0) then 334 dot2r = dotr ; if (abs(dot2r)<tol9) dot2r = zero ! in order to hide the numerical noise 335 dot2i = doti ; if (abs(dot2i)<tol9) dot2i = zero ! in order to hide the numerical noise 336 write(msg,'(3a,2(a,es17.8E3))') bra_i,op2,ket_j,' = ',dot2r,',',dot2i 337 call wrtout(std_out,msg) 338 end if 339 340 rf2%amn(:,iband+(jband-1)*nband_k) = factor_in*(/dotr,doti/) + rf2%amn(:,iband+(jband-1)*nband_k) 341 342 end if ! end choice 343 344 DBG_EXIT("COLL") 345 346 end subroutine rf2_accumulate_bands
m_rf2/rf2_apply_hamiltonian [ Functions ]
[ Top ] [ m_rf2 ] [ Functions ]
NAME
rf2_apply_hamiltonian
FUNCTION
Apply the KS Hamiltonian (or derivative) to an input wave function. If asked, it also does some checks.
INPUTS
rf2 : rf2_t object containing all rf2 data cg_jband (used only if debug_mode/=0) : array containing |u^(0)(jband)> for all bands cprj_jband(natom,nspinor*usecprj)= u^(0) wave functions for all bands projected with non-local projectors: cprj_jband=<p_i|u^(0)(jband)> cwave(2,size_wf) : input wave function |u> cwaveprj(natom,nspinor) : input wave function |u> projected with non-local projectors <p_i|u> eig0 : 0-order eigenvalue for the present wavefunction at k [enl]=optional (if not present, use hamk%ekb); non-local coeffs connecting projectors eig1_k_jband : first-order lagrange multipliers for the band j (Lambda^(1)_ji for all i) [ffnl1]=nonlocal form factors (needed for tests of ipert>=natom+10) [ffnl1_test]=nonlocal form factors used for tests (needed for tests of ipert>=natom+10) jband : band index of the input wavefunction gs_hamkq <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k+q gvnlx1(2,npw1*nspinor)= part of <G|K^(1)+Vnl^(1)|C> not depending on VHxc^(1) (sij_opt/=-1) or part of <G|K^(1)+Vnl^(1)-lambda.S^(1)|C> not depending on VHxc^(1) (sij_opt==-1) idir=direction of the perturbation ipert=type of the perturbation of the Hamiltonian : ipert = 0 : GS calculation, call of getghc ipert = natom+1,2 : 1st order calculation, call of getgh1c ipert = natom+10,11 : 2nd order calculation, call of getgh2c ikpt=number of the k-point isppol=1 index of current spin component mkmem =number of k points trated by this node (GS data). mpi_enreg=information about MPI parallelization nband_k=number of bands at this k point for that spin polarization nsppol=1 for unpolarized, 2 for spin-polarized debug_mode : if /=0 : some tests are done (see NOTES below). Wrong results are printed in std_out prtvol=control print volume and debugging output (for getghc) rf_hamk_idir <type(rf_hamiltonian_type)>=all data for the 1st-order Hamiltonian at k,q (here q=0) size_cprj=size of a cprj array (=gs_hamkq%nspinor) size_wf=size of a wavefunction array (=gs_hamkq%npw_k*gs_hamkq%nspinor)
OUTPUT
h_cwave(2,size_wf) : array containing H^(ipert)|cwave> s_cwave(2,size_wf) : array containing S^(ipert)|cwave>
NOTES
* Tests are done if debug_mode/=0. In that case : cg_jband(:,jband,1) = |u^(0)(jband)>
For ipert==natom+2 (electric field) : cg_jband(:,jband,2) = i * |du/dk_idir(jband)>
According to ipert, we check that :
-- ipert = 0 :
< u^(0)(jband)| H^(0) | u^(0)(jband)> = eig0(jband)
-- ipert = natom+1 or 2 :
< u^(0)(jband)| ( H^(1) - (eig0(jband)+eig0(iband))/2 * S^(1) ) | u^(0)(iband)> = Lambda^(1)(jband,iband)
--ipert >= natom+10 :
< u^(0) | H^(2) | u^(0) > from getgh2c is equal to nonlop with signs=1
* Use of cprj :
The cprj array is computed in dfpt_looppert.F90. In the context of rf2 calculation, it contains
<Proj_i^(0)|u^(0)> and <Proj_i^(1)|u^(0)> (dir=1,2 and 3) for all wavefunctions u^(0) or
<Proj_i^(0)|u^(1)> and <Proj_i^(1)|u^(1)> (dir=1,2 and 3) for all 1st-order wavefunctions u^(1).
Note that <Proj_i^(2)|u^(0)> is always computed on the fly.
SOURCE
415 subroutine rf2_apply_hamiltonian(cg_jband,cprj_jband,cwave,cwaveprj,h_cwave,s_cwave,eig0,eig1_k_jband,& 416 & jband,gs_hamkq,gvnlx1,idir,ipert,ikpt,isppol,mkmem,mpi_enreg,nband_k,nsppol,& 417 & debug_mode,prtvol,rf_hamk_idir,size_cprj,size_wf,& 418 & conj,enl,ffnl1,ffnl1_test) ! optional 419 420 !Arguments --------------------------------------------- 421 !scalars 422 logical,intent(in),optional :: conj 423 integer,intent(in) :: idir,ipert,ikpt,isppol,jband,mkmem,nband_k,nsppol,debug_mode,prtvol,size_wf,size_cprj 424 type(gs_hamiltonian_type),intent(inout) :: gs_hamkq 425 ! type(rf2_t),intent(in) :: rf2 426 type(rf_hamiltonian_type),intent(inout),target :: rf_hamk_idir 427 type(MPI_type),intent(in) :: mpi_enreg 428 429 !arrays 430 real(dp),intent(in),target :: cg_jband(2,size_wf*debug_mode*nband_k,2) 431 real(dp),intent(in),optional,target :: enl(gs_hamkq%dimekb1,gs_hamkq%dimekb2,gs_hamkq%nspinor**2,gs_hamkq%dimekbq) 432 real(dp),intent(in),optional :: ffnl1(:,:,:,:),ffnl1_test(:,:,:,:) 433 real(dp),intent(in) :: eig0(nband_k),eig1_k_jband(2*nband_k) 434 real(dp),intent(inout) :: gvnlx1(2,size_wf) 435 real(dp),intent(inout) :: cwave(2,size_wf),h_cwave(2,size_wf),s_cwave(2,size_wf) 436 type(pawcprj_type),intent(inout),target :: cprj_jband(:,:),cwaveprj(:,:) 437 438 !Local variables --------------------------------------- 439 !scalars 440 integer,parameter :: berryopt=0,tim_getghc=0,tim_getgh1c=0,tim_getgh2c=0,tim_nonlop=0 441 integer,parameter :: alpha(9)=(/1,2,3,2,1,1,3,3,2/),beta(9)=(/1,2,3,3,3,2,2,1,1/) 442 integer :: choice,cpopt,iatom,iband,idirc,idir1,idir2,idir_dum,natom,nnlout,paw_opt,signs,sij_opt 443 integer :: opt_gvnlx1,opt_gvnl2,optlocal,optnl,usevnl 444 logical :: compute_conjugate,has_cprj_jband,has_cwaveprj,pert_phon_elfd 445 real(dp) :: dotr,doti,dotr2,doti2,enlout(18*gs_hamkq%natom),tol_test 446 real(dp), pointer :: enl_ptr(:,:,:,:) 447 real(dp),allocatable,target :: enl_temp(:,:,:,:) 448 character(len=500) :: msg 449 450 !arrays 451 real(dp),target :: cwave_empty(0,0),svectout_dum(0,0) 452 real(dp),allocatable :: gvnlxc(:,:),iddk(:,:) 453 real(dp), ABI_CONTIGUOUS pointer :: cwave_i(:,:),cwave_j(:,:) 454 type(pawcprj_type),target :: cprj_empty(0,0) 455 type(pawcprj_type),pointer :: cprj_j(:,:) 456 ! ********************************************************************* 457 458 DBG_ENTER("COLL") 459 460 compute_conjugate = .false. 461 if(present(conj)) compute_conjugate = conj 462 463 ABI_UNUSED(ikpt) 464 ABI_UNUSED(isppol) 465 ABI_UNUSED(mkmem) 466 ABI_UNUSED(nsppol) 467 468 !Check sizes 469 if (size(cprj_jband)/=0) then 470 if (size(cprj_jband)/=nband_k*gs_hamkq%natom*gs_hamkq%nspinor*gs_hamkq%usecprj) then 471 write(msg,'(2(a,i10))') & 472 'Wrong cprj size! actual size = ',size(cprj_jband),& 473 ' good size = ',nband_k*gs_hamkq%natom*gs_hamkq%nspinor*gs_hamkq%usecprj 474 ABI_BUG(msg) 475 end if 476 end if 477 if (size(cwaveprj)/=0) then 478 if (size(cwaveprj)/=gs_hamkq%natom*gs_hamkq%nspinor*gs_hamkq%usecprj) then 479 write(msg,'(2(a,i10))') & 480 'Wrong cwaveprj size! actual size = ',size(cwaveprj),& 481 ' good size = ',gs_hamkq%natom*gs_hamkq%nspinor*gs_hamkq%usecprj 482 ABI_BUG(msg) 483 end if 484 end if 485 486 natom = gs_hamkq%natom 487 488 pert_phon_elfd = .false. 489 if (ipert>natom+11.and.ipert<=2*natom+11) pert_phon_elfd = .true. 490 usevnl = 0 491 opt_gvnlx1 = 0 492 opt_gvnl2 = 0 493 optnl = 2 494 sij_opt=1;if (gs_hamkq%usepaw==0) sij_opt=0 495 if (ipert==natom+2.or.(ipert==natom+11.and.gs_hamkq%usepaw==1).or.pert_phon_elfd) then 496 usevnl = 1 497 opt_gvnlx1 = 2 498 opt_gvnl2 = 1 499 end if 500 optlocal = 1 501 if (pert_phon_elfd) then 502 optlocal = 0 503 sij_opt = 0 504 optnl = 1 505 end if 506 tol_test=tol8 507 508 ! In the PAW case : manage cprj_jband, cwaveprj 509 has_cprj_jband=(gs_hamkq%usepaw==1.and.gs_hamkq%usecprj==1.and.size(cprj_jband)/=0) 510 has_cwaveprj=(gs_hamkq%usepaw==1.and.gs_hamkq%usecprj==1.and.size(cwaveprj)/=0) 511 cprj_j => cprj_empty 512 513 if (debug_mode/=0) then 514 if (ipert/=0) then 515 write(msg,'(5(a,i4))') 'RF2 TEST rf2_apply_hamiltonian ipert = ',ipert,' idir =',idir,' isppol =',isppol,& 516 & ' ikpt = ',ikpt,' jband = ',jband 517 else 518 write(msg,'(3(a,i4))') 'RF2 TEST rf2_apply_hamiltonian ipert = 0 idir = 0 isppol =',isppol,& 519 & ' ikpt = ',ikpt,' jband = ',jband 520 end if 521 call wrtout(std_out,msg) 522 end if 523 524 ! ******************************************************************************************* 525 ! apply H^(0) 526 ! ******************************************************************************************* 527 if (ipert == 0) then 528 529 ABI_MALLOC(gvnlxc,(2,size_wf)) 530 gvnlxc(:,:) = zero 531 532 ! Test if < u^(0) | H^(0) | u^(0) > = eig0(jband) 533 if(debug_mode/=0) then 534 cwave_j => cg_jband(:,1+(jband-1)*size_wf:jband*size_wf,1) 535 if (has_cprj_jband) cprj_j => cprj_jband(:,1+(jband-1)*size_cprj:jband*size_cprj) 536 cpopt = -1+3*gs_hamkq%usecprj*gs_hamkq%usepaw 537 call getghc(cpopt,cwave_j,cprj_j,h_cwave,s_cwave,gs_hamkq,gvnlxc,zero,mpi_enreg,1,prtvol,& 538 sij_opt,tim_getghc,0,select_k=KPRIME_H_KPRIME) 539 540 call dotprod_g(dotr,doti,gs_hamkq%istwf_k,size_wf,2,cwave_j,h_cwave,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 541 dotr = dotr - eig0(jband) 542 dotr = sqrt(dotr**2+doti**2) 543 if (dotr > tol_test) then 544 write(msg,'(a,es17.8E3)') 'RF2 TEST GETGHC : NOT PASSED dotr = ',dotr 545 call wrtout(ab_out,msg) 546 call wrtout(std_out,msg) 547 end if 548 end if ! end tests 549 550 cpopt=-1;if (has_cwaveprj) cpopt=2 551 call getghc(cpopt,cwave,cwaveprj,h_cwave,s_cwave,gs_hamkq,gvnlxc,zero,mpi_enreg,1,prtvol,& 552 sij_opt,tim_getghc,0,select_k=KPRIME_H_KPRIME) 553 ABI_FREE(gvnlxc) 554 555 ! ******************************************************************************************* 556 ! apply H^(1) 557 ! ******************************************************************************************* 558 else if (ipert<=natom+2) then 559 560 ! Test if < u^(0) | ( H^(1) - eps^(0) S^(1) ) | u^(0) > = eig^(1) 561 if(debug_mode/=0) then 562 ABI_MALLOC(iddk,(2,size_wf)) 563 cwave_j => cg_jband(:,1+(jband-1)*size_wf:jband*size_wf,1) 564 if (has_cprj_jband) cprj_j => cprj_jband(:,1+(jband-1)*size_cprj:jband*size_cprj) 565 iddk(:,:) = zero;if (ipert==natom+2) iddk(:,:)=cg_jband(:,1+(jband-1)*size_wf:jband*size_wf,2) 566 call getgh1c(berryopt,cwave_j,cprj_j,h_cwave,cwave_empty,s_cwave,gs_hamkq,iddk,idir,ipert,zero,& 567 mpi_enreg,optlocal,optnl,opt_gvnlx1,rf_hamk_idir,sij_opt,tim_getgh1c,usevnl,conj=compute_conjugate) 568 do iband=1,nband_k 569 cwave_i => cg_jband(:,1+(iband-1)*size_wf:iband*size_wf,1) 570 call dotprod_g(dotr,doti,gs_hamkq%istwf_k,size_wf,2,cwave_i,h_cwave,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 571 if (gs_hamkq%usepaw==1.and.ipert/=natom+2) then ! S^(1) is zero for ipert=natom+2 572 call dotprod_g(dotr2,doti2,gs_hamkq%istwf_k,size_wf,2,cwave_i,s_cwave,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 573 dotr = dotr - (eig0(iband)+eig0(jband))*dotr2/two 574 doti = doti - (eig0(iband)+eig0(jband))*doti2/two 575 end if 576 dotr = dotr - eig1_k_jband(1+2*(iband-1)) 577 doti = doti - eig1_k_jband(2+2*(iband-1)) 578 dotr = sqrt(dotr**2+doti**2) 579 if (dotr > tol_test) then 580 write(msg,'(4(a,i2),a,es17.8E3)') 'RF2 TEST GETGH1 : ipert=',ipert,' idir=',idir,& 581 ' jband=',jband,' iband=',iband,' NOT PASSED dotr = ',dotr 582 call wrtout(ab_out,msg) 583 call wrtout(std_out,msg) 584 end if 585 end do ! end iband 586 ABI_FREE(iddk) 587 end if ! end tests 588 589 call getgh1c(berryopt,cwave,cwaveprj,h_cwave,cwave_empty,s_cwave,gs_hamkq,gvnlx1,idir,ipert,zero,& 590 mpi_enreg,optlocal,optnl,opt_gvnlx1,rf_hamk_idir,sij_opt,tim_getgh1c,usevnl,conj=compute_conjugate) 591 592 ! ******************************************************************************************* 593 ! apply H^(2) 594 ! ******************************************************************************************* 595 else if (ipert==natom+10.or.ipert==natom+11.or.pert_phon_elfd) then 596 597 ! ******************************************************************************************* 598 ! Test if < u^(0) | H^(2) | u^(0) > from getgh2c is equal to nonlop with signs=1 599 if(debug_mode/=0.and.present(ffnl1).and.present(ffnl1_test)) then 600 601 cwave_j => cg_jband(:,1+(jband-1)*size_wf:jband*size_wf,1) 602 603 ABI_MALLOC(iddk,(2,size_wf)) 604 iddk(:,:) = cwave_j(:,:) 605 606 call getgh2c(cwave_j,cprj_empty,h_cwave,s_cwave,gs_hamkq,iddk,idir,ipert,zero,& 607 mpi_enreg,optlocal,optnl,opt_gvnl2,rf_hamk_idir,sij_opt,tim_getgh2c,usevnl,& 608 conj=compute_conjugate,optkin=0,enl=enl) 609 ABI_FREE(iddk) 610 611 call dotprod_g(dotr,doti,gs_hamkq%istwf_k,size_wf,2,cwave_j,h_cwave,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 612 613 idir1=alpha(idir);idir2=beta(idir) 614 idirc=3*(idir1-1)+idir2 615 616 enlout = zero 617 618 ! Change the pointer ffnl_k to ffnl1_test (idir_ffnl=0, for nonlop with signs=1) 619 call gs_hamkq%load_k(ffnl_k=ffnl1_test) 620 621 signs=1 622 dotr2 = 0 623 doti2 = 0 624 625 ! ************************ 626 ! IPERT == NATOM+10 627 ! ************************ 628 if (ipert==natom+10) then 629 630 cpopt=-1; choice=8; paw_opt=1; if (gs_hamkq%usepaw==0) paw_opt=0 ; nnlout = 6 631 632 call nonlop(choice,cpopt,cprj_empty,enlout,gs_hamkq,idirc,(/zero/),mpi_enreg,1,nnlout,& 633 & paw_opt,signs,svectout_dum,tim_nonlop,cwave_j,svectout_dum) 634 dotr2 = enlout(idir) 635 636 end if ! IPERT == NATOM+10 637 638 ! ************************ 639 ! IPERT == NATOM+11 640 ! ************************ 641 if (ipert==natom+11) then 642 643 if (present(enl)) then 644 enl_ptr => enl 645 else if (associated(rf_hamk_idir%e1kbfr).and.associated(rf_hamk_idir%e1kbsc).and.optnl==2) then 646 ABI_MALLOC(enl_temp,(gs_hamkq%dimekb1,gs_hamkq%dimekb2,gs_hamkq%nspinor**2,rf_hamk_idir%cplex)) 647 enl_temp(:,:,:,:) = rf_hamk_idir%e1kbfr(:,:,:,:) + rf_hamk_idir%e1kbsc(:,:,:,:) 648 enl_ptr => enl_temp 649 else if (associated(rf_hamk_idir%e1kbfr)) then 650 enl_ptr => rf_hamk_idir%e1kbfr 651 end if 652 653 ! Compute application of dS/dk1 to cwavef 654 ! sum_{i,j} s_ij < psi^(0) | d(|p_i><p_j|)/dk(idir1) | psi^(0) > 655 cpopt=-1 ; choice=5 ; paw_opt=3 ; nnlout=3 656 call nonlop(choice,cpopt,cprj_empty,enlout,gs_hamkq,idir_dum,(/zero/),mpi_enreg,1,nnlout,& 657 & paw_opt,signs,svectout_dum,tim_nonlop,cwave_j,svectout_dum) 658 dotr2 = enlout(idir1) 659 660 ! Compute part of H^(2) due to derivative of projectors (idir1) and derivative of Dij (idir2) 661 ! sum_{i,j} chi_ij(idir2) < psi^(0) | d(|p_i><p_j|)/dk(idir1) | psi^(0) > 662 cpopt=-1 ; choice=5 ; paw_opt=1 ; nnlout=3 663 call nonlop(choice,cpopt,cprj_empty,enlout,gs_hamkq,idir_dum,(/zero/),mpi_enreg,1,nnlout,& 664 & paw_opt,signs,svectout_dum,tim_nonlop,cwave_j,svectout_dum,enl=enl_ptr) 665 dotr2 = dotr2 + enlout(idir1) 666 667 ! Compute derivatives due to projectors |d^2[p_i]/dk1dk2>,|d[p_i]/dk1>,|d[p_i]/dk2> 668 ! i * sum_{i,j} < psi^(0) | (d(|p_i><dp_j/dk(idir2)|)/dk(idir1) | psi^(0) > 669 cpopt=-1 ; choice=81 ; paw_opt=3 ; nnlout=18 670 call nonlop(choice,cpopt,cprj_empty,enlout,gs_hamkq,idir_dum,(/zero/),mpi_enreg,1,nnlout,& 671 & paw_opt,signs,svectout_dum,tim_nonlop,cwave_j,svectout_dum) 672 dotr2 = dotr2 - enlout(2*idirc ) 673 doti2 = doti2 + enlout(2*idirc-1) 674 675 end if ! IPERT == NATOM+11 676 677 ! ************************ 678 ! PERT_PHON_ELFD 679 ! ************************ 680 if (pert_phon_elfd) then 681 682 if (present(enl)) then 683 enl_ptr => enl 684 else if (associated(rf_hamk_idir%e1kbfr).and.associated(rf_hamk_idir%e1kbsc).and.optnl==2) then 685 ABI_MALLOC(enl_temp,(gs_hamkq%dimekb1,gs_hamkq%dimekb2,gs_hamkq%nspinor**2,rf_hamk_idir%cplex)) 686 enl_temp(:,:,:,:) = rf_hamk_idir%e1kbfr(:,:,:,:) + rf_hamk_idir%e1kbsc(:,:,:,:) 687 enl_ptr => enl_temp 688 else if (associated(rf_hamk_idir%e1kbfr)) then 689 enl_ptr => rf_hamk_idir%e1kbfr 690 end if 691 692 iatom = gs_hamkq%atindx(ipert - natom - 11) ! Atoms are type-sorted in enlout() 693 694 ! Compute application of dS/dtau1 to i*d[cwavef]/dk2 695 ! sum_{i,j} s_ij < psi^(0) | d(|p_i><p_j|)/dtau(idir1) | i*psi^(k(idir2)) > 696 cpopt=-1 ; choice=2 ; paw_opt=3 ; nnlout=3*natom 697 call nonlop(choice,cpopt,cprj_empty,enlout,gs_hamkq,idir_dum,(/zero/),mpi_enreg,1,nnlout,& 698 & paw_opt,signs,svectout_dum,tim_nonlop,cwave_j,svectout_dum) 699 dotr2 = enlout(3*(iatom-1)+idir1) 700 701 ! Compute part of H^(2) due to derivative of projectors (idir1) and derivative of Dij (idir2) 702 ! sum_{i,j} chi_ij(idir2) < psi^(0) | d(|p_i><p_j|)/dtau(idir1) | psi^(0) > 703 cpopt=-1 ; choice=2 ; paw_opt=1 ; nnlout=3*natom 704 call nonlop(choice,cpopt,cprj_empty,enlout,gs_hamkq,idir_dum,(/zero/),mpi_enreg,1,nnlout,& 705 & paw_opt,signs,svectout_dum,tim_nonlop,cwave_j,svectout_dum,enl=enl_ptr) 706 dotr2 = dotr2 + enlout(3*(iatom-1)+idir1) 707 708 ! Compute derivatives due to projectors |d^2[p_i]/dtau1dk2>,|d[p_i]/dtau1>,|d[p_i]/dk2> 709 ! i * sum_{i,j} < psi^(0) | (d(|p_i><dp_j/dk(idir2)|)/dtau(idir1) | psi^(0) > 710 cpopt=-1 ; choice=54 ; paw_opt=3 ; nnlout=18*natom 711 call nonlop(choice,cpopt,cprj_empty,enlout,gs_hamkq,idir_dum,(/zero/),mpi_enreg,1,nnlout,& 712 & paw_opt,signs,svectout_dum,tim_nonlop,cwave_j,svectout_dum) 713 dotr2 = dotr2 - enlout(18*(iatom-1)+2*idirc ) 714 doti2 = doti2 + enlout(18*(iatom-1)+2*idirc-1) 715 716 end if ! PERT_PHON_ELFD 717 718 dotr = dotr - dotr2 719 doti = doti - doti2 720 dotr = sqrt(dotr**2+doti**2) 721 if (dotr > 10*tol_test) then 722 write(msg,'(a,es17.8E3)') 'RF2 TEST GETGH2C : NOT PASSED dotr = ',dotr 723 call wrtout(ab_out,msg) 724 call wrtout(std_out,msg) 725 end if 726 727 ! Change the pointer ffnl_k back to ffnl1 (idir_ffnl=4, for nonlop with signs=2) 728 call gs_hamkq%load_k(ffnl_k=ffnl1) 729 730 if (associated(enl_ptr)) then 731 nullify(enl_ptr) 732 end if 733 if (allocated(enl_temp)) then 734 ABI_FREE(enl_temp) 735 end if 736 737 end if ! end tests 738 ! ******************************************************************************************* 739 740 call getgh2c(cwave,cwaveprj,h_cwave,s_cwave,gs_hamkq,gvnlx1,idir,ipert,zero,& 741 mpi_enreg,optlocal,optnl,opt_gvnl2,rf_hamk_idir,sij_opt,tim_getgh2c,usevnl,& 742 conj=compute_conjugate,enl=enl) 743 744 else 745 746 h_cwave = zero 747 ABI_ERROR(" rf2_apply_hamiltonian can be used only for: 0<=ipert<=natom+2 and natom+10<=ipert<=2*natom+11.") 748 return 749 750 end if 751 752 DBG_EXIT("COLL") 753 754 end subroutine rf2_apply_hamiltonian
m_rf2/rf2_destroy [ Functions ]
[ Top ] [ m_rf2 ] [ Functions ]
NAME
rf2_init
FUNCTION
Free all allocated arrays in a rf2_t object
INPUTS
rf2 : the rf2_t object to destroy
OUTPUT
NOTES
SOURCE
775 subroutine rf2_destroy(rf2) 776 777 !Arguments --------------------------------------------- 778 !scalars 779 type(rf2_t),intent(inout) :: rf2 780 !arrays 781 782 !Local variables --------------------------------------- 783 !scalars 784 785 ! ************************************************************************* 786 787 if (associated(rf2%RHS_Stern)) then 788 ABI_FREE(rf2%RHS_Stern) 789 end if 790 if (allocated(rf2%dcwavef)) then 791 ABI_FREE(rf2%dcwavef) 792 end if 793 if (allocated(rf2%amn)) then 794 ABI_FREE(rf2%amn) 795 end if 796 if (allocated(rf2%lambda_mn)) then 797 ABI_FREE(rf2%lambda_mn) 798 end if 799 800 end subroutine rf2_destroy
m_rf2/rf2_getidir [ Functions ]
[ Top ] [ m_rf2 ] [ Functions ]
NAME
rf2_getidir
FUNCTION
Get the direction of the 2nd order perturbation according to inputs idir1 and idir2
INPUTS
idir1 : index of the 1st direction (1<=idir1<=3) idir2 : index of the 2nd direction (1<=idir2<=3)
OUTPUT
idir : integer between 1 and 9
NOTES
SOURCE
148 subroutine rf2_getidir(idir1,idir2,idir) 149 150 !Arguments --------------------------------------------- 151 !scalars 152 integer,intent(in) :: idir1,idir2 153 integer,intent(out) :: idir 154 integer :: dir_from_dirs(9) = (/1,6,5,9,2,4,8,7,3/) 155 156 ! ************************************************************************* 157 158 idir = dir_from_dirs(3*(idir1-1)+idir2) 159 160 end subroutine rf2_getidir
m_rf2/rf2_getidirs [ Functions ]
[ Top ] [ m_rf2 ] [ Functions ]
NAME
rf2_getidirs
FUNCTION
Get directions of the 1st and 2nd perturbations according to the input idir
INPUTS
idir : integer between 1 and 9
OUTPUT
idir1 : index of the 1st direction (1<=idir1<=3) idir2 : index of the 2nd direction (1<=idir2<=3)
NOTES
SOURCE
183 subroutine rf2_getidirs(idir,idir1,idir2) 184 185 !Arguments --------------------------------------------- 186 !scalars 187 integer,intent(in) :: idir 188 integer,intent(out) :: idir1,idir2 189 190 ! ************************************************************************* 191 192 select case(idir) 193 ! Diagonal terms : 194 case(1) 195 idir1 = 1 196 idir2 = 1 197 case(2) 198 idir1 = 2 199 idir2 = 2 200 case(3) 201 idir1 = 3 202 idir2 = 3 203 ! Upper triangular terms : 204 case(4) 205 idir1 = 2 206 idir2 = 3 207 case(5) 208 idir1 = 1 209 idir2 = 3 210 case(6) 211 idir1 = 1 212 idir2 = 2 213 ! Lower triangular terms : 214 case(7) 215 idir1 = 3 216 idir2 = 2 217 case(8) 218 idir1 = 3 219 idir2 = 1 220 case(9) 221 idir1 = 2 222 idir2 = 1 223 end select 224 225 end subroutine rf2_getidirs
m_rf2/rf2_t [ Types ]
NAME
rf2_t
FUNCTION
Datatype gathering all needed data for the computation of the 2nd order Sternheimer equation.
SOURCE
54 type,public :: rf2_t 55 56 !scalars 57 integer :: ndir ! number of directions to consider 58 integer :: nband_k ! number of bands 59 integer :: size_wf ! number of components in a wavefunction 60 integer :: size_cprj ! number of components of a cprj variable (i.e. <p_i|WF>) 61 62 !arrays 63 integer :: iperts(2) ! perturbations to compute 64 ! ipert = natom + 10 : 65 ! iperts(1) = iperts(2) = natom+1 (wavevector k) 66 ! 67 ! ipert = natom + 11 : 68 ! iperts(1) = natom+1 (wavevector k) 69 ! iperts(2) = natom+2 (electric field) 70 71 integer :: idirs(2) ! directions of the perturbations (ndir=1 : idirs(1)=idirs(2) , ndir=2 : idirs(1)/=idirs(2)) 72 73 real(dp), ABI_CONTIGUOUS pointer :: RHS_Stern(:,:) => null() ! need a pointer because is a ptr target 74 ! Right-hand side of the 2nd order Sternheimer equation, for every bands. 75 ! Namely, for a band "n" : 76 ! |(RHS_Stern)_n> = (H^(2)-epsilon_n S^(2)) |u^(0)_n> + 2(H^(1)-epsilon_n S^(1))|u^(1)_n> 77 ! - 2 sum_m ( lambda_mn^(1) ( S^(1) |u^(0)_m> + S^(0) |u^(1)_m> ) 78 ! ( /!\ : in the sum, the index m runs over occupied bands ) 79 ! where : 80 ! - epsilon_n = eigenvalue of the GS Hamiltonian 81 ! - lambda_mn^(1) = <u^(0)_m| (H^(1)-(epsilon_n+epsilon_m)/2 S^(1) |u^(0)_n> (1st order Lagrange multiplier) 82 ! - X^(2) = d^2X/(dk_dir1 dlambda_dir2) 83 ! - 2X^(1)Y^(1) = dX/dk_dir1 dY/dlambda_dir2 + dX/dlambda_dir2 dY/dk_dir1 84 ! lambda_dir2 = k_dir2 (ipert==natom+10) , E_dir2 (ipert==natom+11) 85 ! Note : P_c^* will be apply to |(RHS_Stern)_n> in dfpt_cgwf.F90 86 ! ** 87 ! Computed in "rf2_init" 88 89 real(dp),allocatable :: amn(:,:) 90 ! Scalar needed for the "orhtonormalization condition", see "dcwavef(:,:)" below 91 ! Namely : 92 ! A_mn = <u^(0)_m| S^(2) |u^(0)_n> + 2 <u^(1)_m| S^(0) |u^(1)_n> 93 ! + 2 <u^(1)_m| S^(1) |u^(0)_n> + 2 <u^(0)_m| S^(1) |u^(1)_n> 94 ! ** 95 ! Computed in "rf2_init", stored only for testing purposes 96 97 real(dp),allocatable :: dcwavef(:,:) 98 ! Vector needed to enforce the "orhtonormalization condition" on 2nd order wave functions. 99 ! Namely : 100 ! |dcwavef_n> = -1/2 sum_m A_mn |u^(0)_m> 101 ! ** 102 ! Computed in "rf2_init" 103 104 real(dp),allocatable :: lambda_mn(:,:) 105 ! 2nd order Lagrange multiplier. 106 ! Namely : 107 ! lambda_mn = <u^(0)_m| H^(2) |u^(0)_n> + 2 <u^(1)_m| H^(0) |u^(1)_n> 108 ! + 2 <u^(1)_m| H^(1) |u^(0)_n> + 2 <u^(0)_m| H^(1) |u^(1)_n> 109 ! - A_mn * (epsilon_m + epsilon_n) / 2 110 ! ** 111 ! Computed in "rf2_init" 112 113 end type rf2_t 114 115 public :: rf2_getidir 116 public :: rf2_getidirs 117 public :: rf2_accumulate_bands 118 public :: rf2_apply_hamiltonian 119 public :: rf2_destroy