TABLE OF CONTENTS
ABINIT/calc_rpa_functional [ Functions ]
NAME
calc_rpa_functional
FUNCTION
Routine used to calculate the RPA approximation to the correlation energy from the irreducible polarizability.
COPYRIGHT
Copyright (C) 2008-2018 ABINIT group (FB) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
INPUTS
iq=index of the q-point in the array Qmesh%ibz where epsilon^-1 has to be calculated Ep<em1params_t>=Structure with parameters and dimensions related to the inverse dielectric matrix. Pvc<vcoul_t>=Structure gathering data on the Coulombian interaction Qmesh<kmesh_t>=Data type with information on the q-sampling Dtfil<Datafiles_type)>=variables related to files spaceComm=MPI communicator.
OUTPUT
SIDE EFFECTS
NOTES
PARENTS
screening
CHILDREN
coeffs_gausslegint,wrtout,xginv,xheev,xmpi_sum_master
SOURCE
38 #if defined HAVE_CONFIG_H 39 #include "config.h" 40 #endif 41 42 #include "abi_common.h" 43 44 subroutine calc_rpa_functional(gwrpacorr,iqcalc,iq,Ep,Pvc,Qmesh,Dtfil,gmet,chi0,spaceComm,ec_rpa) 45 46 use defs_basis 47 use defs_abitypes 48 use m_xmpi 49 use m_errors 50 use m_profiling_abi 51 52 use m_gwdefs, only : GW_TOLQ0, em1params_t 53 use m_io_tools, only : open_file 54 use m_numeric_tools, only : coeffs_gausslegint 55 use m_abilasi, only : xginv, xheev 56 use m_geometry, only : normv 57 use m_bz_mesh, only : kmesh_t 58 use m_vcoul, only : vcoul_t 59 60 !This section has been created automatically by the script Abilint (TD). 61 !Do not modify the following lines by hand. 62 #undef ABI_FUNC 63 #define ABI_FUNC 'calc_rpa_functional' 64 use interfaces_14_hidewrite 65 !End of the abilint section 66 67 implicit none 68 69 !Arguments ------------------------------------ 70 !scalars 71 integer,intent(in) :: iqcalc,iq,gwrpacorr,spaceComm 72 type(kmesh_t),intent(in) :: Qmesh 73 type(vcoul_t),intent(in) :: Pvc 74 type(Datafiles_type),intent(in) :: Dtfil 75 type(em1params_t),intent(in) :: Ep 76 !arrays 77 real(dp),intent(in) :: gmet(3,3) 78 real(dp),intent(inout) :: ec_rpa(gwrpacorr) 79 complex(gwpc),intent(inout) :: chi0(Ep%npwe,Ep%npwe,Ep%nomega) 80 81 !Local variables------------------------------- 82 !scalars 83 integer :: ig1,ig2,ilambda,io,master,rank,nprocs,unt,ierr 84 real(dp) :: ecorr 85 real(dp) :: lambda 86 logical :: qeq0 87 character(len=500) :: msg 88 !arrays 89 real(dp),allocatable :: z(:),zl(:),zlw(:),zw(:) 90 complex(gwpc),allocatable :: chi0_diag(:),chitmp(:,:) 91 real(gwpc),allocatable :: eig(:) 92 93 ! ************************************************************************* 94 95 DBG_ENTER("COLL") 96 97 ! initialize MPI data 98 master=0 99 rank = xmpi_comm_rank(spaceComm) 100 nprocs = xmpi_comm_size(spaceComm) 101 102 !if (rank==master) then ! presently only master has chi0 in screening 103 104 ! vc_sqrt contains vc^{1/2}(q,G), complex-valued to allow for a possible cutoff 105 qeq0=(normv(Qmesh%ibz(:,iq),gmet,'G')<GW_TOLQ0) 106 107 ! Calculate Gauss-Legendre quadrature knots and weights for the omega integration 108 ABI_ALLOCATE(zw,(Ep%nomegaei)) 109 ABI_ALLOCATE(z,(Ep%nomegaei)) 110 call coeffs_gausslegint(zero,one,z,zw,Ep%nomegaei) 111 112 ! Calculate Gauss-Legendre quadrature knots and weights for the lambda integration 113 ABI_ALLOCATE(zlw,(gwrpacorr)) 114 ABI_ALLOCATE(zl,(gwrpacorr)) 115 call coeffs_gausslegint(zero,one,zl,zlw,gwrpacorr) 116 117 118 ABI_ALLOCATE(chi0_diag,(Ep%npwe)) 119 ABI_STAT_ALLOCATE(chitmp,(Ep%npwe,Ep%npwe), ierr) 120 ABI_CHECK(ierr==0, "out-of-memory in chitmp") 121 122 do io=2,Ep%nomega 123 124 if(gwrpacorr==1) then ! exact integration over the coupling constant 125 126 if(modulo(io-2,nprocs)/=rank) cycle ! distributing the workload 127 128 do ig2=1,Ep%npwe 129 do ig1=1,Ep%npwe 130 chitmp(ig1,ig2) = Pvc%vc_sqrt(ig1,iq) * Pvc%vc_sqrt(ig2,iq) * chi0(ig1,ig2,io) 131 end do !ig1 132 end do !ig2 133 ABI_ALLOCATE(eig,(Ep%npwe)) 134 call xheev('V','U',Ep%npwe,chitmp,eig) 135 136 do ig1=1,Ep%npwe 137 ec_rpa(:) = ec_rpa(:) & 138 & - zw(io-1) / ( z(io-1) * z(io-1) ) & 139 & * Qmesh%wt(iq) * (-log( 1.0_dp-eig(ig1) ) - eig(ig1) ) / (2.0_dp * pi ) 140 end do 141 ABI_DEALLOCATE(eig) 142 143 else ! numerical integration over the coupling constant 144 145 if(modulo( (ilambda-1)+gwrpacorr*(io-2),nprocs)/=rank) cycle ! distributing the workload 146 147 do ilambda=1,gwrpacorr 148 lambda=zl(ilambda) 149 do ig1=1,Ep%npwe 150 chi0_diag(ig1) = Pvc%vc_sqrt(ig1,iq)**2 * chi0(ig1,ig1,io) 151 end do 152 153 do ig2=1,Ep%npwe 154 do ig1=1,Ep%npwe 155 chitmp(ig1,ig2) = - lambda * Pvc%vc_sqrt(ig1,iq) * Pvc%vc_sqrt(ig1,iq) * chi0(ig1,ig2,io) 156 end do !ig1 157 chitmp(ig2,ig2) = chitmp(ig2,ig2) + 1.0_dp 158 end do !ig2 159 call xginv(chitmp(:,:),Ep%npwe) 160 chitmp(:,:) = matmul( chi0(:,:,io) , chitmp(:,:) ) 161 162 do ig1=1,Ep%npwe 163 chi0_diag(ig1) = Pvc%vc_sqrt(ig1,iq) * Pvc%vc_sqrt(ig1,iq) * chitmp(ig1,ig1) - chi0_diag(ig1) 164 end do 165 166 do ig1=1,Ep%npwe 167 ec_rpa(ilambda) = ec_rpa(ilambda) & 168 & - zw(io-1) / ( z(io-1) * z(io-1) ) * Qmesh%wt(iq) * real( chi0_diag(ig1) ) / (2.0_dp * pi ) 169 end do 170 171 end do ! ilambda 172 173 end if ! exact or numerical integration over the coupling constant 174 175 end do ! io 176 177 178 ! Output the correlation energy when the last q-point to be calculated is reached 179 ! This would allow for a manual parallelization over q-points 180 if(iqcalc==Ep%nqcalc) then 181 182 call xmpi_sum_master(ec_rpa,master,spaceComm,ierr) 183 184 if(rank==master) then 185 ecorr = sum( zlw(:)*ec_rpa(:) ) 186 if (open_file(dtfil%fnameabo_rpa, msg, newunit=unt) /=0) then 187 MSG_ERROR(msg) 188 end if 189 write(unt,'(a,(2x,f14.8))') '#RPA',ecorr 190 write(msg,'(2a,(2x,f14.8))') ch10,' RPA energy [Ha] :',ecorr 191 call wrtout(std_out,msg,'COLL') 192 call wrtout(ab_out,msg,'COLL') 193 if(gwrpacorr>1) then 194 do ilambda=1,gwrpacorr 195 write(unt,'(i6,2x,f10.6,2x,e13.6)') ilambda,zl(ilambda),ec_rpa(ilambda) 196 write(msg,'(i6,2x,f10.6,2x,e13.6)') ilambda,zl(ilambda),ec_rpa(ilambda) 197 call wrtout(std_out,msg,'COLL') 198 call wrtout(ab_out,msg,'COLL') 199 end do 200 end if 201 close(unt) 202 end if 203 204 end if 205 206 ABI_DEALLOCATE(chi0_diag) 207 ABI_DEALLOCATE(chitmp) 208 ABI_DEALLOCATE(zl) 209 ABI_DEALLOCATE(zlw) 210 ABI_DEALLOCATE(z) 211 ABI_DEALLOCATE(zw) 212 213 DBG_EXIT("COLL") 214 215 end subroutine calc_rpa_functional