TABLE OF CONTENTS
ABINIT/pspnl_hgh_rec [ Functions ]
NAME
pspnl_hgh_rec
FUNCTION
This routine computes the matrices S_kk'^{l,A} of the projectors
(it is the exponential of the overlap matrix).
it coorresponds to the matrix:
$$\left[(U_l)^{-1}*Exp(-temperature D_l )*U_l* (g_l)^{-1} -Identity\right]_kk'
where (U_l)^-1* D_l* U_l = h^lg_l.
$g_l = <f^l_k|f^l_{k'}>$ is the overlap matrix between projectors
and $h^l_{kk'}$ is the strength matrix of the projectors.
It calulates also the strength eigenvalues and eigenvectors of $h$,
used in the calculus of non-local energy
COPYRIGHT
Copyright (C) 2008-2018 ABINIT group ( ). This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
temperature=4*rtrotter/beta=4*rtrotter*tsmear: the effective temp. in recursion psps <type(pseudopotential_type)>=variables related to pseudo-potentials debug=debug variable
OUTPUT
nlrec%mat_exp_psp_nl=the matrix of the exponential of the projectors:
for any psp, for any angular moment:
h_ij=strength; g_ij=ovelap => exp(-h.g/temp/4p).g^(-1)
nlrec%radii=Local radii of nl psp
nlrec%pspinfo(:,:) for any typat: (momang,typat)=number of projectors
nlrec%eival(:,:,:) for any psp, any momang, the eigenvalues of the
strength matrix H: eival(:,mang,psp)
nlrec%eivec(:,:,:,:)for any psp, any momang, the eigenvectors of the
strength matrix H: eivec(:,:,mang,psp)
SIDE EFFECTS
PARENTS
m_rec
CHILDREN
dgetrf,dgetri,dsyev,exp_mat,gamma_function,set2unit,wrtout
NOTES
at this time :
SOURCE
54 #if defined HAVE_CONFIG_H 55 #include "config.h" 56 #endif 57 58 #include "abi_common.h" 59 60 subroutine pspnl_hgh_rec(psps,temperature,nlrec,debug) 61 62 use m_profiling_abi 63 64 use defs_basis 65 use defs_datatypes 66 use defs_rectypes 67 use m_exp_mat, only : exp_mat 68 use m_numeric_tools, only : set2unit 69 use m_linalg_interfaces 70 71 !This section has been created automatically by the script Abilint (TD). 72 !Do not modify the following lines by hand. 73 #undef ABI_FUNC 74 #define ABI_FUNC 'pspnl_hgh_rec' 75 use interfaces_14_hidewrite 76 use interfaces_28_numeric_noabirule 77 !End of the abilint section 78 79 implicit none 80 81 !Arguments ----------------------------------- 82 !scalars 83 real(dp),intent(in) :: temperature 84 logical,intent(in) :: debug 85 type(pseudopotential_type),intent(in) :: psps 86 type(nlpsprec_type),intent(inout) :: nlrec 87 !arrays 88 !Local variables------------------------------- 89 !scalars 90 integer,parameter :: maxsize=3 91 integer,parameter :: lwork=(1+32)*maxsize 92 integer :: iangol,ipseudo,info,nproj 93 integer :: g_mat_size,ii,nproj2 94 real(dp) :: denom_1,denom_2,tot_proj 95 character(len=500) :: msg 96 !arrays 97 integer :: ipvt(1) 98 !real(dp) :: rwork(2*maxsize) 99 real(dp) :: h_mat_init(3,3), rework(lwork) 100 real(dp), allocatable :: g_mat(:,:),h_mat(:,:),eig_val_h(:) 101 real(dp), allocatable :: identity(:,:),inv_g_mat(:,:),u_mat(:,:) 102 103 complex(dpc),allocatable :: hg_mat(:,:) 104 ! ************************************************************************* 105 106 if(debug)then 107 write(msg,'(80a,a,a)') ('=',ii=1,80),ch10,' pspnl_hgh_rec : enter ' 108 call wrtout(std_out,msg,'COLL') 109 end if 110 111 !--For any pseudo potential: 112 do ipseudo = 1, psps%npsp !--Loop on the pseudos 113 write(msg,'(a,80a)')' pseudo file',('-',ii=1,10) 114 call wrtout(std_out,msg,'COLL') 115 write(msg,'(a)') psps%filpsp(ipseudo) 116 call wrtout(std_out,msg,'COLL') 117 118 ! --For any angular moment: 119 do iangol = 0,psps%mpsang-1 !--Loop on the angular moment 120 121 ! --Local radius 122 nlrec%radii(iangol+1,ipseudo) = psps%gth_params%psppar(iangol+1,0,ipseudo) 123 124 ! --Strenghts of non-local projectors (matrix h) 125 ! --Diagonal part: 126 h_mat_init = zero 127 h_mat_init(1,1) = psps%gth_params%psppar(iangol+1,1,ipseudo) 128 h_mat_init(2,2) = psps%gth_params%psppar(iangol+1,2,ipseudo) 129 h_mat_init(3,3) = psps%gth_params%psppar(iangol+1,3,ipseudo) 130 ! --Out-diagonal part 131 ! --Depending on angular moment the projectors 132 ! strength is calculated differently 133 select case(iangol) 134 case(0) 135 h_mat_init(1,2) = -half*sqrt(3.d0/5.d0)*h_mat_init(2,2) 136 h_mat_init(1,3) = half*sqrt(5.d0/21.d0)*h_mat_init(3,3) 137 h_mat_init(2,3) = -half*sqrt(100.d0/63.d0)*h_mat_init(3,3) 138 case(1) 139 h_mat_init(1,2) = -half*sqrt(5.d0/7.d0)*h_mat_init(2,2) 140 h_mat_init(1,3) = sixth*sqrt(35.d0/11.d0)*h_mat_init(3,3) 141 h_mat_init(2,3) = -14.d0/six/sqrt(11.d0) *h_mat_init(3,3) 142 case(2) 143 h_mat_init(1,2) = -half*sqrt(7.d0/9.d0)*h_mat_init(2,2) 144 h_mat_init(1,3) = half*sqrt(63.d0/143.d0)*h_mat_init(3,3) 145 h_mat_init(2,3) = -nine/sqrt(143.d0)*h_mat_init(3,3) 146 case(3) 147 h_mat_init(1,2) = zero; h_mat_init(1,3) = zero; h_mat_init(2,3) = zero; 148 case default 149 write(msg,'(a)')' error angular: moment component' 150 call wrtout(std_out,msg,'COLL') 151 end select 152 153 154 155 ! --Real dimensions of projectors. 156 g_mat_size = count(abs((/ (h_mat_init(ii,ii),ii=1,3) /))>1.d-8) 157 nlrec%pspinfo(iangol+1,ipseudo) = g_mat_size 158 write(msg,'(a,i2,a,i2)')' ang. moment=',iangol,', N projectors=',g_mat_size 159 call wrtout(std_out,msg,'COLL') 160 if (g_mat_size>0) then 161 ! --Identity matrix 162 ABI_ALLOCATE(identity,(g_mat_size,g_mat_size)) 163 call set2unit(identity) 164 ! identity = zero 165 ! identity(:,1) = one 166 ! identity(:,:) = cshift(array=identity,shift=(/ (-ii,ii=0,g_mat_size) /), dim=2 ) 167 168 169 ! ############## CALCULOUS OF THE EIGEN_SPACE OF THE PROJECTORS STRENGTHS ################## 170 ! --Inverse of the matrix h 171 ABI_ALLOCATE(eig_val_h,(g_mat_size)) 172 ABI_ALLOCATE(u_mat,(g_mat_size,g_mat_size)) 173 ! --u-mat will contain the eigenvectors of h_mat_init 174 u_mat = h_mat_init(:g_mat_size,:g_mat_size) 175 176 ! write(std_out,*)'hmat_init' 177 ! do ii=1,g_mat_size 178 ! write(std_out,*)h_mat_init(ii,:) 179 ! end do 180 call DSYEV('v','u',g_mat_size,u_mat,g_mat_size,eig_val_h,rework,lwork,info) 181 182 ! --THE DIAGONAL MATRIX IS GIVEN BY D=U^t.H.U 183 ! (eival=transpose(eivec).h_mat_init.eivec) 184 write(msg,'(a,3d10.3)')' eigenvalues=',eig_val_h 185 call wrtout(std_out,msg,'COLL') 186 ! write(std_out,*)'autovec';write(std_out,*)u_mat 187 188 nlrec%eival(:g_mat_size,1+iangol,ipseudo) = eig_val_h 189 nlrec%eivec(:g_mat_size,:g_mat_size,1+iangol,ipseudo) = u_mat 190 ABI_DEALLOCATE(eig_val_h) 191 ABI_DEALLOCATE(u_mat) 192 193 ! ##########END CALCULOUS OF THE EIGEN_SPACE OF THE PROJECTORS STRENGTHS ################## 194 195 ABI_ALLOCATE(g_mat,(g_mat_size,g_mat_size)) 196 ABI_ALLOCATE(inv_g_mat,(g_mat_size,g_mat_size)) 197 ABI_ALLOCATE(h_mat,(g_mat_size,g_mat_size)) 198 ABI_ALLOCATE(hg_mat,(g_mat_size,g_mat_size)) 199 200 g_mat(:,:) = one 201 h_mat(:,:) = zero 202 h_mat(:,:) = h_mat_init(:g_mat_size,:g_mat_size) 203 204 ! ------------------------------------------------------- 205 ! --Matrix of the overlap between projetors (matrix g) 206 ! and the h matrix of strength 207 do nproj = 1,g_mat_size-1 208 do nproj2 = 1+nproj,g_mat_size 209 tot_proj = zero 210 ! --Analytic value of overlap 211 ! g_ij=Gamma[-1/2+i+j+l]/Sqrt(Gamma[-1/2+i+iangol]*Gamma[-1/2+j+iangol]) 212 call gamma_function(-half+real(nproj+nproj2+iangol,dp),tot_proj) 213 call gamma_function(-half+real(iangol+2*nproj,dp),denom_1) 214 call gamma_function(-half+real(iangol+2*nproj2,dp),denom_2) 215 216 g_mat(nproj,nproj2) = tot_proj/sqrt(denom_1*denom_2) 217 g_mat(nproj2,nproj) = g_mat(nproj,nproj2) 218 219 h_mat(nproj,nproj2) = h_mat_init(nproj,nproj2) 220 h_mat(nproj2,nproj) = h_mat_init(nproj,nproj2) 221 end do 222 end do 223 224 ! --Inverse of the overlap matrix g 225 inv_g_mat = g_mat 226 call DGETRF(g_mat_size,g_mat_size,inv_g_mat,g_mat_size,ipvt,info) 227 call DGETRI(g_mat_size,inv_g_mat,g_mat_size,ipvt,rework,lwork,info) 228 229 230 ! ----------------------------------------------------------- 231 ! --Now it calculates the exponential of the matrix h.g 232 hg_mat = matmul(h_mat,g_mat) 233 234 call exp_mat(hg_mat,g_mat_size,-one/temperature) 235 236 ! --(exp(h.g)-Identity).(g^-1) 237 hg_mat = hg_mat-identity(:,:) 238 239 240 ! --results on output 241 nlrec%mat_exp_psp_nl(:g_mat_size,:g_mat_size,1+iangol,ipseudo) = matmul(real(hg_mat,dp),inv_g_mat) 242 243 ! write(std_out,*) nlrec%mat_exp_psp_nl(:g_mat_size,:g_mat_size,1+iangol,ipseudo) 244 245 ABI_DEALLOCATE(g_mat) 246 ABI_DEALLOCATE(hg_mat) 247 ABI_DEALLOCATE(h_mat) 248 ABI_DEALLOCATE(inv_g_mat) 249 ABI_DEALLOCATE(identity) 250 end if 251 252 end do !enddo on angular moment 253 end do !enddo on pseudos 254 255 256 if(debug)then 257 write(msg,'(80a,a,a)') ('=',ii=1,80),ch10,' pspnl_hgh_rec : exit ' 258 call wrtout(std_out,msg,'COLL') 259 end if 260 261 262 end subroutine pspnl_hgh_rec