TABLE OF CONTENTS
ABINIT/getgh2c [ Functions ]
NAME
getgh2c
FUNCTION
Compute <G|H^(2)|C> (or <G|H^(2)-Eps.S^(2)|C>) for input vector |C> expressed in reciprocal space. (H^(2) is the 2nd-order pertubed Hamiltonian, S^(2) is the 2nd-order perturbed overlap operator). Result is put in array gh2c. If required, part of <G|K(2)+Vnonlocal^(2)|C> not depending on VHxc^(2) is also returned in gvnl2. If required, <G|S^(2)|C> is returned in gs2c (S=overlap - PAW only)
COPYRIGHT
Copyright (C) 2015-2018 ABINIT group (MT,JLJ) 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
cwavef(2,npw*nspinor)=input wavefunction, in reciprocal space cwaveprj(natom,nspinor*usecprj)=<p_lmn|C> coefficients for wavefunction |C> gs_hamkq <type(gs_hamiltonian_type)>=all data for the Hamiltonian idir=direction of the perturbation ipert=type of the perturbation lambda=real use to apply H^(2)-lambda.S^(2) mpi_enreg=information about MPI parallelization optlocal=0: local part of H^(2) is not computed in gh2c=<G|H^(2)|C> 1: local part of H^(2) is computed in gh2c=<G|H^(2)|C> optnl=0: non-local part of H^(2) is not computed in gh2c=<G|H^(2)|C> 1: non-local part of H^(2) depending on VHxc^(2) is not computed in gh2c=<G|H^(2)|C> 2: non-local part of H^(2) is totally computed in gh2c=<G|H^(2)|C> opt_gvnl2=option controlling the use of gvnl2 array: 0: not used 1: used as input: - used only for PAW and ipert=natom+11 At input: contains the ddk 1-st order WF (times i) rf_hamkq <type(rf_hamiltonian_type)>=all data for the 2nd-order Hamiltonian at k,k+q sij_opt= -PAW ONLY- if 0, only matrix elements <G|H^(2)|C> have to be computed (S=overlap) if 1, matrix elements <G|S^(2)|C> have to be computed in gs2c in addition to gh2c if -1, matrix elements <G|H^(2)-lambda.S^(2)|C> have to be computed in gh2c (gs2c not used) tim_getgh2c=timing code of the calling subroutine (can be set to 0 if not attributed) usevnl=1 if gvnl2=(part of <G|K^(2)+Vnl^(2)-lambda.S^(2)|C> not depending on VHxc^(2)) has to be input/output
OUTPUT
gh2c(2,npw1*nspinor)= <G|H^(2)|C> or <G|H^(2)-lambda.S^(2)|C>
(only kinetic+non-local parts if optlocal=0)
if (usevnl==1)
gvnl2(2,npw1*nspinor)= part of <G|K^(2)+Vnl^(2)|C> not depending on VHxc^(2) (sij_opt/=-1)
or part of <G|K^(2)+Vnl^(2)-lambda.S^(2)|C> not depending on VHxc^(2) (sij_opt==-1)
if (sij_opt=1)
gs2c(2,npw1*nspinor)=<G|S^(2)|C> (S=overlap).
PARENTS
m_rf2
CHILDREN
nonlop,pawcprj_alloc,pawcprj_free
SOURCE
62 #if defined HAVE_CONFIG_H 63 #include "config.h" 64 #endif 65 66 #include "abi_common.h" 67 68 subroutine getgh2c(cwavef,cwaveprj,gh2c,gs2c,gs_hamkq,gvnl2,idir,ipert,lambda,& 69 & mpi_enreg,optlocal,optnl,opt_gvnl2,rf_hamkq,sij_opt,tim_getgh2c,usevnl) 70 71 use defs_basis 72 use defs_abitypes 73 use m_profiling_abi 74 use m_errors 75 76 use m_pawcprj, only : pawcprj_type,pawcprj_alloc,pawcprj_free 77 use m_hamiltonian, only : gs_hamiltonian_type,rf_hamiltonian_type 78 79 !This section has been created automatically by the script Abilint (TD). 80 !Do not modify the following lines by hand. 81 #undef ABI_FUNC 82 #define ABI_FUNC 'getgh2c' 83 use interfaces_66_nonlocal 84 !End of the abilint section 85 86 implicit none 87 88 !Arguments ------------------------------------ 89 !scalars 90 integer,intent(in) :: idir,ipert,optlocal,optnl,opt_gvnl2,sij_opt,tim_getgh2c,usevnl 91 real(dp),intent(in) :: lambda 92 type(MPI_type),intent(in) :: mpi_enreg 93 type(gs_hamiltonian_type),intent(inout),target :: gs_hamkq 94 type(rf_hamiltonian_type),intent(inout),target :: rf_hamkq 95 !arrays 96 real(dp),intent(inout) :: cwavef(:,:) 97 real(dp),intent(inout),target :: gvnl2(:,:) 98 real(dp),intent(out) :: gh2c(:,:),gs2c(:,:) 99 type(pawcprj_type),intent(inout),target :: cwaveprj(:,:) 100 101 !Local variables------------------------------- 102 !scalars 103 integer :: choice,cpopt,idir1,idir2,idirc,ipw,ipws,ispinor,my_nspinor 104 integer :: natom,ncpgr,nnlout=1,npw,npw1,paw_opt,signs,tim_nonlop,usecprj 105 logical :: has_kin,has_vnl 106 real(dp) :: enlout_dum(1) 107 character(len=500) :: msg 108 !arrays 109 ! integer,parameter :: alpha(6)=(/1,2,3,3,3,2/),beta(6)=(/1,2,3,2,1,1/) 110 integer,parameter :: alpha(9)=(/1,2,3,2,1,1,3,3,2/),beta(9)=(/1,2,3,3,3,2,2,1,1/) 111 !real(dp) :: tsec(2) 112 real(dp) :: svectout_dum(1,1),vectout_dum(1,1) 113 real(dp),allocatable :: nonlop_out(:,:) 114 real(dp), pointer :: gvnl2_(:,:) 115 real(dp), pointer :: ddkinpw(:),kinpw1(:) 116 type(pawcprj_type),allocatable,target :: cwaveprj_tmp(:,:) 117 type(pawcprj_type),pointer :: cwaveprj_ptr(:,:) 118 119 ! ********************************************************************* 120 121 !Keep track of total time spent in getgh2c 122 !call timab(196+tim_getgh2c,1,tsec) 123 124 !====================================================================== 125 !== Initialisations and compatibility tests 126 !====================================================================== 127 128 npw =gs_hamkq%npw_k 129 npw1 =gs_hamkq%npw_kp 130 natom=gs_hamkq%natom 131 132 !Compatibility tests 133 if(ipert/=natom+10.and.ipert/=natom+11)then 134 msg='only ipert<>natom+10/natom+11 implemented!' 135 MSG_BUG(msg) 136 end if 137 if (mpi_enreg%paral_spinor==1) then 138 msg='Not compatible with parallelization over spinorial components!' 139 MSG_BUG(msg) 140 end if 141 if (gs_hamkq%nvloc>1) then 142 msg='Not compatible with nvloc=4 (non-coll. magnetism)!' 143 MSG_BUG(msg) 144 end if 145 if(ipert==natom+11.and.gs_hamkq%usepaw==1) then 146 if (optnl>=1.and.((.not.associated(rf_hamkq%e1kbfr)).or.(.not.associated(rf_hamkq%e1kbsc)))) then 147 msg='ekb derivatives must be allocated for ipert=natom+11 !' 148 MSG_BUG(msg) 149 end if 150 if (usevnl==0) then 151 msg='gvnl2 must be allocated for ipert=natom+11 !' 152 MSG_BUG(msg) 153 end if 154 if(opt_gvnl2==0) then 155 msg='opt_gvnl2=0 not compatible with ipert=natom+11 !' 156 MSG_BUG(msg) 157 end if 158 end if 159 160 !Check sizes 161 my_nspinor=max(1,gs_hamkq%nspinor/mpi_enreg%nproc_spinor) 162 if (size(cwavef)<2*npw*my_nspinor) then 163 msg='wrong size for cwavef!' 164 MSG_BUG(msg) 165 end if 166 if (size(gh2c)<2*npw1*my_nspinor) then 167 msg='wrong size for gh2c!' 168 MSG_BUG(msg) 169 end if 170 if (usevnl/=0) then 171 if (size(gvnl2)<2*npw1*my_nspinor) then 172 msg='wrong size for gvnl2!' 173 MSG_BUG(msg) 174 end if 175 end if 176 if (sij_opt==1) then 177 if (size(gs2c)<2*npw1*my_nspinor) then 178 msg='wrong size for gs2c!' 179 MSG_BUG(msg) 180 end if 181 end if 182 183 !PAW: specific treatment for usecprj input arg 184 ! force it to zero if cwaveprj is not allocated 185 usecprj=gs_hamkq%usecprj ; ncpgr=0 186 if(gs_hamkq%usepaw==1) then 187 if (size(cwaveprj)==0) usecprj=0 188 if (usecprj/=0) then 189 ncpgr=cwaveprj(1,1)%ncpgr 190 if (size(cwaveprj)<gs_hamkq%natom*my_nspinor) then 191 msg='wrong size for cwaveprj!' 192 MSG_BUG(msg) 193 end if 194 end if 195 else 196 if(usecprj==1)then 197 msg='usecprj==1 not allowed for NC psps !' 198 MSG_BUG(msg) 199 end if 200 end if 201 202 tim_nonlop=8 203 if (tim_getgh2c==1.and.ipert<=natom) tim_nonlop=7 204 if (tim_getgh2c==2.and.ipert<=natom) tim_nonlop=5 205 if (tim_getgh2c==1.and.ipert> natom) tim_nonlop=8 206 if (tim_getgh2c==2.and.ipert> natom) tim_nonlop=5 207 if (tim_getgh2c==3 ) tim_nonlop=0 208 209 idir1=alpha(idir);idir2=beta(idir) 210 211 !====================================================================== 212 !== Apply the 2nd-order local potential to the wavefunction 213 !====================================================================== 214 215 if (ipert/=natom+10.and.ipert/=natom+11.and.optlocal>0) then 216 msg='local part not implemented' 217 MSG_BUG(msg) 218 else 219 ! In the case of ddk operator, no local contribution (also because no self-consistency) 220 !$OMP PARALLEL DO 221 do ipw=1,npw1*my_nspinor 222 gh2c(:,ipw)=zero 223 end do 224 225 end if 226 227 !====================================================================== 228 !== Apply the 2st-order non-local potential to the wavefunction 229 !====================================================================== 230 231 has_vnl=(ipert==natom+10.or.ipert==natom+11) 232 233 !Use of gvnl2 depends on usevnl 234 if (usevnl==1) then 235 gvnl2_ => gvnl2 236 else 237 ABI_ALLOCATE(gvnl2_,(2,npw1*my_nspinor)) 238 end if 239 240 if (has_vnl.and.(optnl>0.or.sij_opt/=0)) then 241 242 idirc=3*(idir1-1)+idir2 !xx=1, xy=2, xz=3, yx=4, yy=5, yz=6, zx=7, zy=8, zz=9, (xyz,xyz)=(idir1,idir2) 243 244 ! d^2[H_nl]/dk1dk2 245 ! ------------------------------------------- 246 if (ipert==natom+10) then 247 if (gs_hamkq%usepaw==1) then 248 if (usecprj==1) then 249 cwaveprj_ptr => cwaveprj 250 else 251 ABI_DATATYPE_ALLOCATE(cwaveprj_tmp,(natom,my_nspinor)) 252 call pawcprj_alloc(cwaveprj_tmp,0,gs_hamkq%dimcprj) 253 cwaveprj_ptr => cwaveprj_tmp 254 end if 255 cpopt=-1+5*usecprj 256 choice=8; signs=2; paw_opt=1; if (sij_opt/=0) paw_opt=sij_opt+3 257 call nonlop(choice,cpopt,cwaveprj_ptr,enlout_dum,gs_hamkq,idirc,(/lambda/),mpi_enreg,1,nnlout,& 258 & paw_opt,signs,gs2c,tim_nonlop,cwavef,gvnl2_) 259 if (usecprj==0) then 260 call pawcprj_free(cwaveprj_tmp) 261 ABI_DATATYPE_DEALLOCATE(cwaveprj_tmp) 262 end if 263 nullify(cwaveprj_ptr) 264 else 265 choice=8; signs=2; cpopt=-1 ; paw_opt=0 266 call nonlop(choice,cpopt,cwaveprj,enlout_dum,gs_hamkq,idirc,(/zero/),mpi_enreg,1,nnlout,& 267 & paw_opt,signs,svectout_dum,tim_nonlop,cwavef,gvnl2_) 268 end if 269 270 ! d^2[H_nl]/dk1dE2 : Non-zero only in PAW 271 ! ------------------------------------------- 272 else if (ipert==natom+11.and.gs_hamkq%usepaw==1) then 273 274 ABI_ALLOCATE(nonlop_out,(2,npw1*my_nspinor)) 275 276 if (usecprj==1) then 277 cwaveprj_ptr => cwaveprj 278 else 279 ABI_DATATYPE_ALLOCATE(cwaveprj_tmp,(natom,my_nspinor)) 280 call pawcprj_alloc(cwaveprj_tmp,2,gs_hamkq%dimcprj) 281 cwaveprj_ptr => cwaveprj_tmp 282 end if 283 284 if (opt_gvnl2==1.and.optnl>=1) then 285 286 ! Compute application of dS/dk1 to i*d[cwavef]/dk2 287 cpopt=-1 ; choice=5 ; paw_opt=3 ; signs=2 288 call nonlop(choice,cpopt,cwaveprj_ptr,enlout_dum,gs_hamkq,idir1,(/zero/),mpi_enreg,1,nnlout,& 289 & paw_opt,signs,nonlop_out,tim_nonlop,gvnl2_,vectout_dum) 290 !$OMP PARALLEL DO 291 do ipw=1,npw1*my_nspinor 292 gvnl2_(:,ipw)=nonlop_out(:,ipw) 293 end do 294 295 ! Compute part of H^(2) due to derivative of projectors (idir1) and derivative of Dij (idir2) 296 cpopt=4*usecprj 297 choice=5 ; paw_opt=1 ; signs=2 298 call nonlop(choice,cpopt,cwaveprj_ptr,enlout_dum,gs_hamkq,idir1,(/zero/),mpi_enreg,1,nnlout,& 299 & paw_opt,signs,svectout_dum,tim_nonlop,cwavef,nonlop_out,& 300 & enl=rf_hamkq%e1kbfr+rf_hamkq%e1kbsc) 301 !$OMP PARALLEL DO 302 do ipw=1,npw1*my_nspinor 303 gvnl2_(:,ipw)=gvnl2_(:,ipw)+nonlop_out(:,ipw) 304 end do 305 306 end if ! opt_gvnl2==1 307 308 ! Compute derivatives due to projectors |d^2[p_i]/dk1dk2>,|d[p_i]/dk1>,|d[p_i]/dk2> 309 cpopt=-1+5*usecprj 310 choice=81 ; paw_opt=3 ; signs=2 311 call nonlop(choice,cpopt,cwaveprj_ptr,enlout_dum,gs_hamkq,idirc,(/lambda/),mpi_enreg,1,nnlout,& 312 & paw_opt,signs,nonlop_out,tim_nonlop,cwavef,vectout_dum) 313 !$OMP PARALLEL DO 314 do ipw=1,npw1*my_nspinor ! Note the multiplication by i 315 gvnl2_(1,ipw)=gvnl2_(1,ipw)-nonlop_out(2,ipw) 316 gvnl2_(2,ipw)=gvnl2_(2,ipw)+nonlop_out(1,ipw) 317 end do 318 ABI_DEALLOCATE(nonlop_out) 319 if (sij_opt==1) gs2c=zero 320 if (usecprj==0) then 321 call pawcprj_free(cwaveprj_tmp) 322 ABI_DATATYPE_DEALLOCATE(cwaveprj_tmp) 323 end if 324 nullify(cwaveprj_ptr) 325 end if ! ipert==natom+11 and gs_hamkq%usepaw==1 326 327 !No non-local part 328 !------------------------------------------- 329 else 330 331 if (optnl>=1) then 332 !$OMP PARALLEL DO 333 do ipw=1,npw1*my_nspinor 334 gvnl2_(:,ipw)=zero 335 end do 336 end if 337 if (sij_opt/=0) then 338 !$OMP PARALLEL DO 339 do ipw=1,npw1*my_nspinor 340 gs2c(:,ipw)=zero 341 end do 342 end if 343 344 end if 345 346 !====================================================================== 347 !== Apply the 2nd-order kinetic operator to the wavefunction 348 !====================================================================== 349 350 has_kin=(ipert==natom+10) 351 352 !k-point perturbation 353 !------------------------------------------- 354 if (associated(gs_hamkq%kinpw_kp)) then 355 kinpw1 => gs_hamkq%kinpw_kp 356 else if (optnl>=1.or.has_kin) then 357 msg='need kinpw1 allocated!' 358 MSG_BUG(msg) 359 end if 360 if (associated(rf_hamkq%ddkinpw_k)) then 361 ddkinpw => rf_hamkq%ddkinpw_k 362 else if (has_kin) then 363 msg='need ddkinpw allocated!' 364 MSG_BUG(msg) 365 end if 366 367 if (has_kin) then 368 do ispinor=1,my_nspinor 369 !$OMP PARALLEL DO PRIVATE(ipw,ipws) SHARED(cwavef,ispinor,gvnl2_,ddkinpw,kinpw1,npw,my_nspinor) 370 do ipw=1,npw 371 ipws=ipw+npw*(ispinor-1) 372 if(kinpw1(ipw)<huge(zero)*1.d-11)then 373 gvnl2_(1,ipws)=gvnl2_(1,ipws)+ddkinpw(ipw)*cwavef(1,ipws) 374 gvnl2_(2,ipws)=gvnl2_(2,ipws)+ddkinpw(ipw)*cwavef(2,ipws) 375 else 376 gvnl2_(1,ipws)=zero 377 gvnl2_(2,ipws)=zero 378 end if 379 end do 380 end do 381 end if 382 383 !====================================================================== 384 !== Sum contributions to get the application of H^(2) to the wf 385 !====================================================================== 386 !Also filter the wavefunctions for large modified kinetic energy 387 388 !Add non-local+kinetic to local part 389 if (optnl>=1.or.has_kin) then 390 do ispinor=1,my_nspinor 391 ipws=(ispinor-1)*npw1 392 !$OMP PARALLEL DO PRIVATE(ipw) SHARED(gh2c,gvnl2_,kinpw1,ipws,npw1) 393 do ipw=1+ipws,npw1+ipws 394 if(kinpw1(ipw-ipws)<huge(zero)*1.d-11)then 395 gh2c(1,ipw)=gh2c(1,ipw)+gvnl2_(1,ipw) 396 gh2c(2,ipw)=gh2c(2,ipw)+gvnl2_(2,ipw) 397 else 398 gh2c(1,ipw)=zero 399 gh2c(2,ipw)=zero 400 end if 401 end do 402 end do 403 end if 404 405 if (usevnl==1) then 406 nullify(gvnl2_) 407 else 408 ABI_DEALLOCATE(gvnl2_) 409 end if 410 411 !call timab(196+tim_getgh2c,2,tsec) 412 413 end subroutine getgh2c