TABLE OF CONTENTS
ABINIT/getgh1c [ Functions ]
NAME
getgh1c
FUNCTION
Compute <G|H^(1)|C> (or <G|H^(1)-lambda.S^(1)|C>) for input vector |C> expressed in reciprocal space. (H^(1) is the 1st-order pertubed Hamiltonian, S^(1) is the 1st-order perturbed overlap operator). Result is put in array gh1c. If required, part of <G|K(1)+Vnonlocal^(1)|C> not depending on VHxc^(1) is also returned in gvnl1c. If required, <G|S^(1)|C> is returned in gs1c (S=overlap - PAW only)
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (XG, DRH, MT, SPr) 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
berryopt=option for Berry phase cwave(2,npw*nspinor)=input wavefunction, in reciprocal space cwaveprj(natom,nspinor*usecprj)=<p_lmn|C> coefficients for wavefunction |C> (and 1st derivatives) if not allocated or size=0, they are locally computed (and not sotred) dkinpw(npw)=derivative of the (modified) kinetic energy for each plane wave at k (Hartree) grad_berry(2,npw1*nspinor*(berryopt/4))= the gradient of the Berry phase term gs_hamkq <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k+q idir=direction of the perturbation ipert=type of the perturbation lambda=real use to apply H^(1)-lambda.S^(1) mpi_enreg=information about MPI parallelization npw=number of planewaves in basis sphere at given k. npw1=number of planewaves in basis sphere at k+q optlocal=0: local part of H^(1) is not computed in gh1c=<G|H^(1)|C> 1: local part of H^(1) is computed in gh1c=<G|H^(1)|C> optnl=0: non-local part of H^(1) is not computed in gh1c=<G|H^(1)|C> 1: non-local part of H^(1) depending on VHxc^(1) is not computed in gh1c=<G|H^(1)|C> 2: non-local part of H^(1) is totally computed in gh1c=<G|H^(1)|C> opt_gvnl1=option controlling the use of gvnl1 array: 0: used as an output 1: used as an input: (only for ipert=natom+2) NCPP: contains the ddk 1-st order WF PAW: contains frozen part of 1st-order hamiltonian 2: used as input/ouput: - used only for PAW and ipert=natom+2 At input: contains the ddk 1-st order WF (times i) At output: contains frozen part of 1st-order hamiltonian rf_hamkq <type(rf_hamiltonian_type)>=all data for the 1st-order Hamiltonian at k,k+q sij_opt= -PAW ONLY- if 0, only matrix elements <G|H^(1)|C> have to be computed (S=overlap) if 1, matrix elements <G|S^(1)|C> have to be computed in gs1c in addition to gh1c if -1, matrix elements <G|H^(1)-lambda.S^(1)|C> have to be computed in gh1c (gs1c not used) tim_getgh1c=timing code of the calling subroutine (can be set to 0 if not attributed) usevnl=1 if gvnl1=(part of <G|K^(1)+Vnl^(1)-lambda.S^(1)|C> not depending on VHxc^(1)) has to be input/output
OUTPUT
gh1c(2,npw1*nspinor)= <G|H^(1)|C> or <G|H^(1)-lambda.S^(1)|C> on the k+q sphere (only kinetic+non-local parts if optlocal=0) if (usevnl==1) gvnl1(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) if (sij_opt=1) gs1c(2,npw1*nspinor)=<G|S^(1)|C> (S=overlap) on the k+q sphere.
PARENTS
dfpt_cgwf,dfpt_nstpaw,dfpt_nstwf,dfpt_wfkfermi,m_gkk,m_phgamma,m_phpi m_rf2,m_sigmaph
CHILDREN
kpgstr,load_k_hamiltonian,load_k_rf_hamiltonian,load_kprime_hamiltonian mkffnl,mkkin,mkkpg
SOURCE
74 #if defined HAVE_CONFIG_H 75 #include "config.h" 76 #endif 77 78 #include "abi_common.h" 79 80 subroutine getgh1c(berryopt,cwave,cwaveprj,gh1c,grad_berry,gs1c,gs_hamkq,& 81 & gvnl1,idir,ipert,lambda,mpi_enreg,optlocal,optnl,opt_gvnl1,& 82 & rf_hamkq,sij_opt,tim_getgh1c,usevnl) 83 84 use defs_basis 85 use defs_abitypes 86 use m_profiling_abi 87 use m_errors 88 89 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_copy 90 use m_hamiltonian, only : gs_hamiltonian_type,rf_hamiltonian_type 91 use m_kg, only : kpgstr, mkkin 92 93 !This section has been created automatically by the script Abilint (TD). 94 !Do not modify the following lines by hand. 95 #undef ABI_FUNC 96 #define ABI_FUNC 'getgh1c' 97 use interfaces_18_timing 98 use interfaces_53_ffts 99 use interfaces_66_nonlocal 100 !End of the abilint section 101 102 implicit none 103 104 !Arguments ------------------------------------ 105 !scalars 106 integer,intent(in) :: berryopt,idir,ipert,optlocal,optnl,opt_gvnl1,sij_opt,tim_getgh1c,usevnl 107 real(dp),intent(in) :: lambda 108 type(MPI_type),intent(in) :: mpi_enreg 109 type(gs_hamiltonian_type),intent(inout),target :: gs_hamkq 110 type(rf_hamiltonian_type),intent(inout),target :: rf_hamkq 111 !arrays 112 real(dp),intent(in) :: grad_berry(:,:) 113 real(dp),intent(inout) :: cwave(2,gs_hamkq%npw_k*gs_hamkq%nspinor) 114 real(dp),intent(out) :: gh1c(2,gs_hamkq%npw_kp*gs_hamkq%nspinor) 115 real(dp),intent(out) :: gs1c(2,gs_hamkq%npw_kp*gs_hamkq%nspinor) 116 real(dp),intent(inout),target :: gvnl1(2,gs_hamkq%npw_kp*gs_hamkq%nspinor) 117 type(pawcprj_type),intent(inout),target :: cwaveprj(:,:) 118 119 !Local variables------------------------------- 120 !scalars 121 integer,parameter :: level=16 122 integer :: choice,cplex1,cpopt,ipw,ipws,ispinor,istr,i1,i2,i3 123 integer :: my_nspinor,natom,ncpgr,nnlout=1,npw,npw1,paw_opt,signs 124 integer :: tim_fourwf,tim_nonlop,usecprj 125 logical :: has_kin,usevnl2 126 real(dp) :: weight 127 character(len=500) :: msg 128 !arrays 129 real(dp) :: enlout(1),tsec(2),svectout_dum(1,1),vectout_dum(1,1) 130 real(dp),allocatable :: cwave_sp(:,:),cwavef1(:,:),cwavef2(:,:) 131 real(dp),allocatable :: gh1c_sp(:,:),gh1c1(:,:),gh1c2(:,:),gh1c3(:,:),gh1c4(:,:),gvnl2(:,:) 132 real(dp),allocatable :: nonlop_out(:,:),vlocal1_tmp(:,:,:),work(:,:,:,:) 133 real(dp),ABI_CONTIGUOUS pointer :: gvnl1_(:,:) 134 real(dp),pointer :: dkinpw(:),kinpw1(:) 135 type(pawcprj_type),allocatable,target :: cwaveprj_tmp(:,:) 136 type(pawcprj_type),pointer :: cwaveprj_ptr(:,:) 137 138 ! ********************************************************************* 139 140 !Keep track of total time spent in getgh1c 141 call timab(196+tim_getgh1c,1,tsec) 142 143 !====================================================================== 144 !== Initialisations and compatibility tests 145 !====================================================================== 146 147 npw =gs_hamkq%npw_k 148 npw1 =gs_hamkq%npw_kp 149 natom=gs_hamkq%natom 150 151 !Compatibility tests 152 if(gs_hamkq%usepaw==1.and.(ipert>=0.and.(ipert<=natom.or.ipert==natom+3.or.ipert==natom+4))) then 153 if ((optnl>=1.and.(.not.associated(rf_hamkq%e1kbfr))).or. & 154 & (optnl>=2.and.(.not.associated(rf_hamkq%e1kbsc))))then 155 msg='ekb derivatives must be allocated for ipert<=natom or natom+3/4 !' 156 MSG_BUG(msg) 157 end if 158 end if 159 if(gs_hamkq%usepaw==1.and.(ipert==natom+2)) then 160 if ((optnl>=1.and.(.not.associated(rf_hamkq%e1kbfr))).or. & 161 & (optnl>=2.and.(.not.associated(rf_hamkq%e1kbsc))))then 162 msg='ekb derivatives must be allocated for ipert=natom+2 !' 163 MSG_BUG(msg) 164 end if 165 if (usevnl==0) then 166 msg='gvnl1 must be allocated for ipert=natom+2 !' 167 MSG_BUG(msg) 168 end if 169 end if 170 if(ipert==natom+2.and.opt_gvnl1==0) then 171 msg='opt_gvnl1=0 not compatible with ipert=natom+2 !' 172 MSG_BUG(msg) 173 end if 174 if (mpi_enreg%paral_spinor==1) then 175 msg='Not compatible with parallelization over spinorial components !' 176 MSG_BUG(msg) 177 end if 178 179 !Check sizes 180 my_nspinor=max(1,gs_hamkq%nspinor/mpi_enreg%nproc_spinor) 181 if (size(cwave)<2*npw*my_nspinor) then 182 msg='wrong size for cwave!' 183 MSG_BUG(msg) 184 end if 185 if (size(gh1c)<2*npw1*my_nspinor) then 186 msg='wrong size for gh1c!' 187 MSG_BUG(msg) 188 end if 189 if (usevnl/=0) then 190 if (size(gvnl1)<2*npw1*my_nspinor) then 191 msg='wrong size for gvnl1!' 192 MSG_BUG(msg) 193 end if 194 end if 195 if (sij_opt==1) then 196 if (size(gs1c)<2*npw1*my_nspinor) then 197 msg='wrong size for gs1c!' 198 MSG_BUG(msg) 199 end if 200 end if 201 if (berryopt>=4) then 202 if (size(grad_berry)<2*npw1*my_nspinor) then 203 msg='wrong size for grad_berry!' 204 MSG_BUG(msg) 205 end if 206 end if 207 208 !PAW: specific treatment for usecprj input arg 209 ! force it to zero if cwaveprj is not allocated 210 usecprj=gs_hamkq%usecprj ; ncpgr=0 211 if(gs_hamkq%usepaw==1) then 212 if (size(cwaveprj)==0) usecprj=0 213 if (usecprj/=0) then 214 ncpgr=cwaveprj(1,1)%ncpgr 215 if (size(cwaveprj)<gs_hamkq%natom*my_nspinor) then 216 msg='wrong size for cwaveprj!' 217 MSG_BUG(msg) 218 end if 219 if(gs_hamkq%usepaw==1.and.(ipert>=0.and.(ipert<=natom.or.ipert==natom+3.or.ipert==natom+4))) then 220 if (ncpgr/=1)then 221 msg='Projected WFs (cprj) derivatives are not correctly stored !' 222 MSG_BUG(msg) 223 end if 224 end if 225 end if 226 else 227 if(usecprj==1)then 228 msg='usecprj==1 not allowed for NC psps !' 229 MSG_BUG(msg) 230 end if 231 end if 232 233 tim_nonlop=8 234 if (tim_getgh1c==1.and.ipert<=natom) tim_nonlop=7 235 if (tim_getgh1c==2.and.ipert<=natom) tim_nonlop=5 236 if (tim_getgh1c==1.and.ipert> natom) tim_nonlop=8 237 if (tim_getgh1c==2.and.ipert> natom) tim_nonlop=5 238 if (tim_getgh1c==3 ) tim_nonlop=0 239 240 !====================================================================== 241 !== Apply the 1st-order local potential to the wavefunction 242 !====================================================================== 243 244 !Phonon perturbation 245 !or Electric field perturbation 246 !or Strain perturbation 247 !------------------------------------------- 248 if (ipert<=natom+5.and.ipert/=natom+1.and.optlocal>0) then !SPr deb 249 250 ABI_ALLOCATE(work,(2,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6)) 251 252 if (gs_hamkq%nvloc==1) then 253 254 weight=one ; tim_fourwf=4 255 call fourwf(rf_hamkq%cplex,rf_hamkq%vlocal1,cwave,gh1c,work,gs_hamkq%gbound_k,gs_hamkq%gbound_kp,& 256 & gs_hamkq%istwf_k,gs_hamkq%kg_k,gs_hamkq%kg_kp,gs_hamkq%mgfft,mpi_enreg,1,gs_hamkq%ngfft,& 257 & npw,npw1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,2,mpi_enreg%paral_kgb,tim_fourwf,weight,weight,& 258 & use_gpu_cuda=gs_hamkq%use_gpu_cuda) 259 if(gs_hamkq%nspinor==2)then 260 ABI_ALLOCATE(cwave_sp,(2,npw)) 261 ABI_ALLOCATE(gh1c_sp,(2,npw1)) 262 !$OMP PARALLEL DO 263 do ipw=1,npw 264 cwave_sp(1,ipw)=cwave(1,ipw+npw) 265 cwave_sp(2,ipw)=cwave(2,ipw+npw) 266 end do 267 call fourwf(rf_hamkq%cplex,rf_hamkq%vlocal1,cwave_sp,gh1c_sp,work,gs_hamkq%gbound_k,gs_hamkq%gbound_kp,& 268 & gs_hamkq%istwf_k,gs_hamkq%kg_k,gs_hamkq%kg_kp,gs_hamkq%mgfft,mpi_enreg,1,gs_hamkq%ngfft,& 269 & npw,npw1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,2,mpi_enreg%paral_kgb,tim_fourwf,weight,weight,& 270 & use_gpu_cuda=gs_hamkq%use_gpu_cuda) 271 !$OMP PARALLEL DO 272 do ipw=1,npw1 273 gh1c(1,ipw+npw1)=gh1c_sp(1,ipw) 274 gh1c(2,ipw+npw1)=gh1c_sp(2,ipw) 275 end do 276 ABI_DEALLOCATE(cwave_sp) 277 ABI_DEALLOCATE(gh1c_sp) 278 end if 279 else ! Non-Collinear magnetism for nvloc=4 280 if (gs_hamkq%nspinor==2) then 281 weight=one ; tim_fourwf=4 282 ABI_ALLOCATE(gh1c1,(2,npw1)) 283 ABI_ALLOCATE(gh1c2,(2,npw1)) 284 ABI_ALLOCATE(gh1c3,(2,npw1)) 285 ABI_ALLOCATE(gh1c4,(2,npw1)) 286 gh1c1(:,:)=zero; gh1c2(:,:)=zero; gh1c3(:,:)=zero ; gh1c4(:,:)=zero 287 ABI_ALLOCATE(vlocal1_tmp,(rf_hamkq%cplex*gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6)) !SPr: notation/dimension corrected vlocal_tmp -> vlocal1_tmp 288 ABI_ALLOCATE(cwavef1,(2,npw)) 289 ABI_ALLOCATE(cwavef2,(2,npw)) 290 do ipw=1,npw 291 cwavef1(1:2,ipw)=cwave(1:2,ipw) 292 cwavef2(1:2,ipw)=cwave(1:2,ipw+npw) 293 end do 294 ! gh1c1=v11*phi1 295 vlocal1_tmp(:,:,:)=rf_hamkq%vlocal1(:,:,:,1) 296 call fourwf(rf_hamkq%cplex,vlocal1_tmp,cwavef1,gh1c1,work,gs_hamkq%gbound_k,gs_hamkq%gbound_kp,& 297 & gs_hamkq%istwf_k,gs_hamkq%kg_k,gs_hamkq%kg_kp,gs_hamkq%mgfft,mpi_enreg,1,gs_hamkq%ngfft,& 298 & npw,npw1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,2,mpi_enreg%paral_kgb,tim_fourwf,weight,weight,& 299 & use_gpu_cuda=gs_hamkq%use_gpu_cuda) 300 ! gh1c2=v22*phi2 301 vlocal1_tmp(:,:,:)=rf_hamkq%vlocal1(:,:,:,2) 302 call fourwf(rf_hamkq%cplex,vlocal1_tmp,cwavef2,gh1c2,work,gs_hamkq%gbound_k,gs_hamkq%gbound_kp,& 303 & gs_hamkq%istwf_k,gs_hamkq%kg_k,gs_hamkq%kg_kp,gs_hamkq%mgfft,mpi_enreg,1,gs_hamkq%ngfft,& 304 & npw,npw1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,2,mpi_enreg%paral_kgb,tim_fourwf,weight,weight,& 305 & use_gpu_cuda=gs_hamkq%use_gpu_cuda) 306 ABI_DEALLOCATE(vlocal1_tmp) 307 cplex1=2 308 ABI_ALLOCATE(vlocal1_tmp,(cplex1*gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6)) 309 ! gh1c3=(re(v12)-im(v12))*phi1 => v^21*phi1 310 if(rf_hamkq%cplex==1) then 311 do i3=1,gs_hamkq%n6 312 do i2=1,gs_hamkq%n5 313 do i1=1,gs_hamkq%n4 314 vlocal1_tmp(2*i1-1,i2,i3)= rf_hamkq%vlocal1(i1,i2,i3,3) 315 vlocal1_tmp(2*i1 ,i2,i3)=-rf_hamkq%vlocal1(i1,i2,i3,4) 316 end do 317 end do 318 end do 319 else 320 !SPr: modified definition of local potential components for cplex=2 (see dotprod_vn) 321 !also, v21==v12* not always holds (e.g. magnetic field perturbation) 322 do i3=1,gs_hamkq%n6 323 do i2=1,gs_hamkq%n5 324 do i1=1,gs_hamkq%n4 325 vlocal1_tmp(2*i1-1,i2,i3)= rf_hamkq%vlocal1(2*i1 ,i2,i3,4) 326 vlocal1_tmp(2*i1 ,i2,i3)=-rf_hamkq%vlocal1(2*i1-1,i2,i3,4) 327 end do 328 end do 329 end do 330 end if 331 call fourwf(cplex1,vlocal1_tmp,cwavef1,gh1c3,work,gs_hamkq%gbound_k,gs_hamkq%gbound_kp,& 332 & gs_hamkq%istwf_k,gs_hamkq%kg_k,gs_hamkq%kg_kp,gs_hamkq%mgfft,mpi_enreg,1,gs_hamkq%ngfft,& 333 & npw,npw1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,2,mpi_enreg%paral_kgb,tim_fourwf,weight,weight,& 334 & use_gpu_cuda=gs_hamkq%use_gpu_cuda) 335 ! gh1c4=(re(v12)+im(v12))*phi2 => v^12*phi2 336 if(rf_hamkq%cplex==1) then 337 do i3=1,gs_hamkq%n6 338 do i2=1,gs_hamkq%n5 339 do i1=1,gs_hamkq%n4 340 vlocal1_tmp(2*i1,i2,i3)=-vlocal1_tmp(2*i1,i2,i3) 341 end do 342 end do 343 end do 344 else 345 !for cplex=2 and time-reversal breaking perturbations,v21/=v12* 346 do i3=1,gs_hamkq%n6 347 do i2=1,gs_hamkq%n5 348 do i1=1,gs_hamkq%n4 349 vlocal1_tmp(2*i1-1,i2,i3)= rf_hamkq%vlocal1(2*i1-1,i2,i3,3) 350 vlocal1_tmp(2*i1 ,i2,i3)= rf_hamkq%vlocal1(2*i1 ,i2,i3,3) 351 end do 352 end do 353 end do 354 end if 355 call fourwf(cplex1,vlocal1_tmp,cwavef2,gh1c4,work,gs_hamkq%gbound_k,gs_hamkq%gbound_kp,& 356 & gs_hamkq%istwf_k,gs_hamkq%kg_k,gs_hamkq%kg_kp,gs_hamkq%mgfft,mpi_enreg,1,gs_hamkq%ngfft,& 357 & npw,npw1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,2,mpi_enreg%paral_kgb,tim_fourwf,weight,weight,& 358 & use_gpu_cuda=gs_hamkq%use_gpu_cuda) 359 ABI_DEALLOCATE(vlocal1_tmp) 360 ! Build gh1c from pieces 361 ! gh1c_1 = (v11, v12) (psi1) matrix vector product 362 ! gh1c_2 = (v12*,v22) (psi2) 363 do ipw=1,npw1 364 gh1c(1:2,ipw) =gh1c1(1:2,ipw)+gh1c4(1:2,ipw) 365 gh1c(1:2,ipw+npw1)=gh1c3(1:2,ipw)+gh1c2(1:2,ipw) 366 end do 367 ABI_DEALLOCATE(gh1c1) 368 ABI_DEALLOCATE(gh1c2) 369 ABI_DEALLOCATE(gh1c3) 370 ABI_DEALLOCATE(gh1c4) 371 ABI_DEALLOCATE(cwavef1) 372 ABI_DEALLOCATE(cwavef2) 373 else 374 msg='nspinor/=1 for Non-collinear calculations!' 375 MSG_BUG(msg) 376 end if 377 end if ! nvloc 378 379 ABI_DEALLOCATE(work) 380 381 ! k-point perturbation (or no local part, i.e. optlocal=0) 382 ! ------------------------------------------- 383 else if (ipert==natom+1.or.optlocal==0) then 384 385 ! In the case of ddk operator, no local contribution (also because no self-consistency) 386 !$OMP PARALLEL DO 387 do ipw=1,npw1*my_nspinor 388 gh1c(:,ipw)=zero 389 end do 390 391 end if 392 393 !====================================================================== 394 !== Apply the 1st-order non-local potential to the wavefunction 395 !====================================================================== 396 397 !Use of gvnl1 depends on usevnl 398 if (usevnl==1) then 399 gvnl1_ => gvnl1 400 else 401 ABI_ALLOCATE(gvnl1_,(2,npw1*my_nspinor)) 402 end if 403 usevnl2=(gs_hamkq%usepaw==1.and.optnl>=2.and.& 404 & ((ipert>0.and.(ipert<=natom.or.ipert==natom+3.or.ipert==natom+4)).or.(ipert==natom+2))) 405 406 !Phonon perturbation 407 !------------------------------------------- 408 if (ipert<=natom.and.(optnl>0.or.sij_opt/=0)) then 409 410 ! PAW: 411 if (gs_hamkq%usepaw==1) then 412 413 if (usecprj==1) then 414 cwaveprj_ptr => cwaveprj 415 else 416 ABI_DATATYPE_ALLOCATE(cwaveprj_tmp,(natom,my_nspinor)) 417 call pawcprj_alloc(cwaveprj_tmp,1,gs_hamkq%dimcprj) 418 cwaveprj_ptr => cwaveprj_tmp 419 end if 420 421 ! 1- Compute derivatives due to projectors |p_i>^(1) 422 ! Only displaced atom contributes 423 cpopt=-1+5*usecprj ; choice=2 ; signs=2 424 paw_opt=1;if (sij_opt/=0) paw_opt=sij_opt+3 425 call nonlop(choice,cpopt,cwaveprj_ptr,enlout,gs_hamkq,idir,(/lambda/),mpi_enreg,1,nnlout,& 426 & paw_opt,signs,gs1c,tim_nonlop,cwave,gvnl1_,iatom_only=ipert) 427 ! 2- Compute derivatives due to frozen part of D_ij^(1) (independent of VHxc^(1)) 428 ! All atoms contribute 429 if (optnl>=1) then 430 ABI_ALLOCATE(gvnl2,(2,npw1*my_nspinor)) 431 cpopt=1+3*usecprj ; choice=1 ; signs=2 ; paw_opt=1 432 call nonlop(choice,cpopt,cwaveprj_ptr,enlout,gs_hamkq,idir,(/lambda/),mpi_enreg,1,nnlout,& 433 & paw_opt,signs,svectout_dum,tim_nonlop,cwave,gvnl2,enl=rf_hamkq%e1kbfr) 434 !$OMP PARALLEL DO 435 do ipw=1,npw1*my_nspinor 436 gvnl1_(:,ipw)=gvnl1_(:,ipw)+gvnl2(:,ipw) 437 end do 438 end if 439 ! 3- Compute derivatives due to part of D_ij^(1) depending on VHxc^(1) 440 ! All atoms contribute 441 if (optnl>=2) then 442 cpopt=4 ; choice=1 ; signs=2 ; paw_opt=1 443 call nonlop(choice,cpopt,cwaveprj_ptr,enlout,gs_hamkq,idir,(/lambda/),mpi_enreg,1,nnlout,& 444 & paw_opt,signs,svectout_dum,tim_nonlop,cwave,gvnl2,enl=rf_hamkq%e1kbsc) 445 else if (optnl==1) then 446 ABI_DEALLOCATE(gvnl2) 447 end if 448 449 if (usecprj==0) then 450 call pawcprj_free(cwaveprj_tmp) 451 ABI_DATATYPE_DEALLOCATE(cwaveprj_tmp) 452 end if 453 nullify(cwaveprj_ptr) 454 455 ! Norm-conserving psps: 456 else 457 ! Compute only derivatives due to projectors |p_i>^(1) 458 cpopt=-1 ; choice=2 ; signs=2 ; paw_opt=0 459 call nonlop(choice,cpopt,cwaveprj,enlout,gs_hamkq,idir,(/lambda/),mpi_enreg,1,nnlout,& 460 & paw_opt,signs,svectout_dum,tim_nonlop,cwave,gvnl1_,iatom_only=ipert) 461 if (sij_opt==1) then 462 !$OMP PARALLEL DO 463 do ipw=1,npw1*my_nspinor 464 gs1c(:,ipw)=zero 465 end do 466 end if 467 end if 468 469 ! k-point perturbation 470 ! ------------------------------------------- 471 else if (ipert==natom+1.and.(optnl>0.or.sij_opt/=0)) then 472 473 ! Remember, q=0, so can take all RF data... 474 tim_nonlop=8 ; signs=2 ; choice=5 475 if (gs_hamkq%usepaw==1) then 476 if (usecprj==1) then 477 cwaveprj_ptr => cwaveprj 478 else 479 ABI_DATATYPE_ALLOCATE(cwaveprj_tmp,(natom,my_nspinor)) 480 call pawcprj_alloc(cwaveprj_tmp,1,gs_hamkq%dimcprj) 481 cwaveprj_ptr => cwaveprj_tmp 482 end if 483 cpopt=-1+5*usecprj; paw_opt=1; if (sij_opt/=0) paw_opt=sij_opt+3 484 ! JLJ: BUG (wrong result) of H^(1) if stored cprj are used in PAW DDKs with nspinor==2 (==1 works fine). 485 ! To be debugged, if someone has time... 486 if(gs_hamkq%nspinor==2) cpopt=-1 487 call nonlop(choice,cpopt,cwaveprj_ptr,enlout,gs_hamkq,idir,(/lambda/),mpi_enreg,1,nnlout,& 488 & paw_opt,signs,gs1c,tim_nonlop,cwave,gvnl1_) 489 if (usecprj==0) then 490 call pawcprj_free(cwaveprj_tmp) 491 ABI_DATATYPE_DEALLOCATE(cwaveprj_tmp) 492 end if 493 nullify(cwaveprj_ptr) 494 else 495 cpopt=-1 ; paw_opt=0 496 call nonlop(choice,cpopt,cwaveprj,enlout,gs_hamkq,idir,(/lambda/),mpi_enreg,1,nnlout,& 497 & paw_opt,signs,svectout_dum,tim_nonlop,cwave,gvnl1_) 498 end if 499 500 ! Electric field perturbation without Berry phase 501 ! ------------------------------------------- 502 else if(ipert==natom+2 .and. & 503 & (berryopt/=4 .and. berryopt/=6 .and. berryopt/=7 .and. & 504 & berryopt/=14 .and. berryopt/=16 .and. berryopt/=17) .and.(optnl>0.or.sij_opt/=0))then 505 ! gvnl1 was already initialized in the calling routine, by reading a ddk file 506 ! It contains |i du^(0)/dk_band> 507 508 if (gs_hamkq%usepaw==1) then 509 if (usecprj==1) then 510 cwaveprj_ptr => cwaveprj 511 else 512 ABI_DATATYPE_ALLOCATE(cwaveprj_tmp,(natom,my_nspinor)) 513 call pawcprj_alloc(cwaveprj_tmp,1,gs_hamkq%dimcprj) 514 cwaveprj_ptr => cwaveprj_tmp 515 end if 516 if (opt_gvnl1==2.and.optnl>=1) then 517 518 ! PAW: Compute application of S^(0) to ddk WF 519 cpopt=-1 ; choice=1 ; paw_opt=3 ; signs=2 520 ABI_ALLOCATE(nonlop_out,(2,npw1*my_nspinor)) 521 call nonlop(choice,cpopt,cwaveprj_ptr,enlout,gs_hamkq,0,(/lambda/),mpi_enreg,1,nnlout,& 522 & paw_opt,signs,nonlop_out,tim_nonlop,gvnl1_,vectout_dum) 523 !$OMP PARALLEL DO 524 do ipw=1,npw1*my_nspinor 525 gvnl1_(:,ipw)=nonlop_out(:,ipw) 526 end do 527 528 ! PAW: Compute part of H^(1) due to derivative of S 529 cpopt=4*usecprj ; choice=51 ; paw_opt=3 ; signs=2 530 call nonlop(choice,cpopt,cwaveprj_ptr,enlout,gs_hamkq,idir,(/lambda/),mpi_enreg,1,nnlout,& 531 & paw_opt,signs,nonlop_out,tim_nonlop,cwave,vectout_dum) 532 533 ! Note the multiplication by i 534 !$OMP PARALLEL DO 535 do ipw=1,npw1*my_nspinor 536 gvnl1_(1,ipw)=gvnl1_(1,ipw)-nonlop_out(2,ipw) 537 gvnl1_(2,ipw)=gvnl1_(2,ipw)+nonlop_out(1,ipw) 538 end do 539 540 ! PAW: Compute part of H^(1) due to derivative of electric field part of Dij 541 cpopt=2 ; choice=1 ; paw_opt=1 ; signs=2 542 call nonlop(choice,cpopt,cwaveprj_ptr,enlout,gs_hamkq,0,(/lambda/),mpi_enreg,1,nnlout,& 543 & paw_opt,signs,svectout_dum,tim_nonlop,cwave,nonlop_out,enl=rf_hamkq%e1kbfr) 544 !$OMP PARALLEL DO 545 do ipw=1,npw1*my_nspinor 546 gvnl1_(:,ipw)=gvnl1_(:,ipw)+nonlop_out(:,ipw) 547 end do 548 ABI_DEALLOCATE(nonlop_out) 549 550 end if ! opt_gvnl1==2 551 552 ! PAW: Compute derivatives due to part of D_ij^(1) depending on VHxc^(1) 553 if (optnl>=2) then 554 ABI_ALLOCATE(gvnl2,(2,npw1*my_nspinor)) 555 cpopt=-1+3*usecprj;if (opt_gvnl1==2) cpopt=2 556 choice=1 ; paw_opt=1 ; signs=2 557 call nonlop(choice,cpopt,cwaveprj_ptr,enlout,gs_hamkq,0,(/lambda/),mpi_enreg,1,nnlout,& 558 & paw_opt,signs,svectout_dum,tim_nonlop,cwave,gvnl2,enl=rf_hamkq%e1kbsc) 559 end if 560 if (sij_opt==1) then 561 !$OMP PARALLEL DO 562 do ipw=1,npw1*my_nspinor 563 gs1c(:,ipw)=zero 564 end do 565 end if 566 if (usecprj==0) then 567 call pawcprj_free(cwaveprj_tmp) 568 ABI_DATATYPE_DEALLOCATE(cwaveprj_tmp) 569 end if 570 nullify(cwaveprj_ptr) 571 end if ! PAW 572 573 ! Electric field perturbation with Berry phase 574 ! ------------------------------------------- 575 else if(ipert==natom+2 .and. & 576 & (berryopt==4 .or. berryopt==6 .or. berryopt==7 .or. & 577 & berryopt==14 .or. berryopt==16 .or. berryopt==17 ) .and.(optnl>0.or.sij_opt/=0))then 578 579 if (optnl>=1) then 580 do ipw=1,npw1*my_nspinor 581 gvnl1_(1,ipw)=-grad_berry(2,ipw) 582 gvnl1_(2,ipw)= grad_berry(1,ipw) 583 end do 584 end if 585 if (sij_opt==1) then 586 !$OMP PARALLEL DO 587 do ipw=1,npw1*my_nspinor 588 gs1c(:,ipw)=zero 589 end do 590 end if 591 592 ! Strain perturbation 593 ! ------------------------------------------- 594 else if ((ipert==natom+3.or.ipert==natom+4).and.(optnl>0.or.sij_opt/=0)) then 595 596 ! Remember, q=0, so can take all RF data 597 signs=2 ; istr=idir;if(ipert==natom+4) istr=istr+3 598 599 ! PAW: 600 if (gs_hamkq%usepaw==1) then 601 602 if (usecprj==1) then 603 cwaveprj_ptr => cwaveprj 604 else 605 ABI_DATATYPE_ALLOCATE(cwaveprj_tmp,(natom,my_nspinor)) 606 call pawcprj_alloc(cwaveprj_tmp,1,gs_hamkq%dimcprj) 607 cwaveprj_ptr => cwaveprj_tmp 608 end if 609 610 ! 1- Compute derivatives due to projectors |p_i>^(1) 611 ! All atoms contribute 612 cpopt=-1+5*usecprj ; choice=3 613 paw_opt=1;if (sij_opt/=0) paw_opt=sij_opt+3 614 call nonlop(choice,cpopt,cwaveprj_ptr,enlout,gs_hamkq,istr,(/lambda/),mpi_enreg,1,nnlout,& 615 & paw_opt,signs,gs1c,tim_nonlop,cwave,gvnl1_) 616 ! 2- Compute derivatives due to frozen part of D_ij^(1) (independent of VHxc^(1)) 617 ! All atoms contribute 618 if (optnl>=1) then 619 ABI_ALLOCATE(gvnl2,(2,npw1*my_nspinor)) 620 gvnl2 = zero 621 cpopt=1+3*usecprj ; choice=1 ; paw_opt=1 622 call nonlop(choice,cpopt,cwaveprj_ptr,enlout,gs_hamkq,istr,(/lambda/),mpi_enreg,1,nnlout,& 623 & paw_opt,signs,svectout_dum,tim_nonlop,cwave,gvnl2,& 624 & enl=rf_hamkq%e1kbfr) 625 !$OMP PARALLEL DO 626 do ipw=1,npw1*my_nspinor 627 gvnl1_(:,ipw)=gvnl1_(:,ipw)+gvnl2(:,ipw) 628 end do 629 end if 630 ! 3- Compute derivatives due to part of D_ij^(1) depending on VHxc^(1) 631 ! All atoms contribute 632 if (optnl>=2) then 633 cpopt=4 ; choice=1 ; paw_opt=1 634 call nonlop(choice,cpopt,cwaveprj_ptr,enlout,gs_hamkq,istr,(/lambda/),mpi_enreg,1,nnlout,& 635 & paw_opt,signs,svectout_dum,tim_nonlop,cwave,gvnl2,enl=rf_hamkq%e1kbsc) 636 637 else if (optnl==1) then 638 ABI_DEALLOCATE(gvnl2) 639 end if 640 if (usecprj==0) then 641 call pawcprj_free(cwaveprj_tmp) 642 ABI_DATATYPE_DEALLOCATE(cwaveprj_tmp) 643 end if 644 nullify(cwaveprj_ptr) 645 646 ! Norm-conserving psps: 647 else 648 ! Compute only derivatives due to projectors |p_i>^(1) 649 choice=3 ; cpopt=-1 ; paw_opt=0 650 call nonlop(choice,cpopt,cwaveprj,enlout,gs_hamkq,istr,(/lambda/),mpi_enreg,1,nnlout,& 651 & paw_opt,signs,svectout_dum,tim_nonlop,cwave,gvnl1_) 652 if (sij_opt==1) then 653 !$OMP PARALLEL DO 654 do ipw=1,npw1*my_nspinor 655 gs1c(:,ipw)=zero 656 end do 657 end if 658 end if 659 660 ! No non-local part 661 ! ------------------------------------------- 662 else if (usevnl>0.or.(sij_opt/=0)) then 663 664 if (optnl>=1) then 665 !$OMP PARALLEL DO 666 do ipw=1,npw1*my_nspinor 667 gvnl1_(:,ipw)=zero 668 end do 669 end if 670 if (sij_opt/=0) then 671 !$OMP PARALLEL DO 672 do ipw=1,npw1*my_nspinor 673 gs1c(:,ipw)=zero 674 end do 675 end if 676 677 end if 678 679 !====================================================================== 680 !== Apply the 1st-order kinetic operator to the wavefunction 681 !== (add it to nl contribution) 682 !====================================================================== 683 684 !Phonon perturbation or Electric field perturbation 685 !------------------------------------------- 686 !No kinetic contribution 687 688 !k-point perturbation or Strain perturbation 689 !------------------------------------------- 690 691 has_kin=(ipert==natom+1.or.ipert==natom+3.or.ipert==natom+4) 692 if (associated(gs_hamkq%kinpw_kp)) then 693 kinpw1 => gs_hamkq%kinpw_kp 694 else if (optnl>=1.or.usevnl2.or.has_kin) then 695 msg='need kinpw1 allocated!' 696 MSG_BUG(msg) 697 end if 698 if (associated(rf_hamkq%dkinpw_k)) then 699 dkinpw => rf_hamkq%dkinpw_k 700 else if (has_kin) then 701 msg='need dkinpw allocated!' 702 MSG_BUG(msg) 703 end if 704 705 if (has_kin) then 706 ! Remember that npw=npw1 for ddk perturbation 707 do ispinor=1,my_nspinor 708 !$OMP PARALLEL DO PRIVATE(ipw,ipws) SHARED(cwave,ispinor,gvnl1_,dkinpw,kinpw1,npw,my_nspinor) 709 do ipw=1,npw 710 ipws=ipw+npw*(ispinor-1) 711 if(kinpw1(ipw)<huge(zero)*1.d-11)then 712 gvnl1_(1,ipws)=gvnl1_(1,ipws)+dkinpw(ipw)*cwave(1,ipws) 713 gvnl1_(2,ipws)=gvnl1_(2,ipws)+dkinpw(ipw)*cwave(2,ipws) 714 else 715 gvnl1_(1,ipws)=zero 716 gvnl1_(2,ipws)=zero 717 end if 718 end do 719 end do 720 end if 721 722 !====================================================================== 723 !== Sum contributions to get the application of H^(1) to the wf 724 !====================================================================== 725 !Also filter the wavefunctions for large modified kinetic energy 726 727 !Add non-local+kinetic to local part 728 if (optnl>=1.or.has_kin) then 729 do ispinor=1,my_nspinor 730 ipws=(ispinor-1)*npw1 731 !$OMP PARALLEL DO PRIVATE(ipw) SHARED(gh1c,gvnl1_,kinpw1,ipws,npw1) 732 do ipw=1+ipws,npw1+ipws 733 if(kinpw1(ipw-ipws)<huge(zero)*1.d-11)then 734 gh1c(1,ipw)=gh1c(1,ipw)+gvnl1_(1,ipw) 735 gh1c(2,ipw)=gh1c(2,ipw)+gvnl1_(2,ipw) 736 else 737 gh1c(1,ipw)=zero 738 gh1c(2,ipw)=zero 739 end if 740 end do 741 end do 742 end if 743 744 !PAW: add non-local part due to first order change of VHxc 745 if (usevnl2) then 746 do ispinor=1,my_nspinor 747 ipws=(ispinor-1)*npw1 748 !$OMP PARALLEL DO PRIVATE(ipw) SHARED(gh1c,gvnl2,kinpw1,ipws,npw1) 749 do ipw=1+ipws,npw1+ipws 750 if(kinpw1(ipw-ipws)<huge(zero)*1.d-11)then 751 gh1c(1,ipw)=gh1c(1,ipw)+gvnl2(1,ipw) 752 gh1c(2,ipw)=gh1c(2,ipw)+gvnl2(2,ipw) 753 end if 754 end do 755 end do 756 ABI_DEALLOCATE(gvnl2) 757 end if 758 759 if (usevnl==1) then 760 nullify(gvnl1_) 761 else 762 ABI_DEALLOCATE(gvnl1_) 763 end if 764 765 call timab(196+tim_getgh1c,1,tsec) 766 767 end subroutine getgh1c
m_hamiltonian/getgh1c_setup [ Functions ]
[ Top ] [ m_hamiltonian ] [ Functions ]
NAME
getgh1c_setup
FUNCTION
INPUTS
OUTPUT
PARENTS
dfpt_vtorho,m_gkk,m_phgamma,m_phpi,m_sigmaph
CHILDREN
kpgstr,load_k_hamiltonian,load_k_rf_hamiltonian,load_kprime_hamiltonian mkffnl,mkkin,mkkpg
SOURCE
922 subroutine getgh1c_setup(gs_hamkq,rf_hamkq,dtset,psps,kpoint,kpq,idir,ipert,& ! In 923 & natom,rmet,gprimd,gmet,istwf_k,npw_k,npw1_k,& ! In 924 & useylmgr1,kg_k,ylm_k,kg1_k,ylm1_k,ylmgr1_k,& ! In 925 & dkinpw,nkpg,nkpg1,kpg_k,kpg1_k,kinpw1,ffnlk,ffnl1,ph3d,ph3d1,& ! Out 926 & ddkinpw,dkinpw2,rf_hamk_dir2) ! Optional 927 928 use defs_basis 929 use defs_datatypes 930 use defs_abitypes 931 use m_profiling_abi 932 use m_errors 933 934 use m_kg, only : mkkin, kpgstr 935 use m_hamiltonian, only : gs_hamiltonian_type, rf_hamiltonian_type,& 936 & load_k_hamiltonian,load_kprime_hamiltonian,& 937 & load_k_rf_hamiltonian 938 939 !This section has been created automatically by the script Abilint (TD). 940 !Do not modify the following lines by hand. 941 #undef ABI_FUNC 942 #define ABI_FUNC 'getgh1c_setup' 943 use interfaces_66_nonlocal 944 !End of the abilint section 945 946 implicit none 947 948 !Arguments ------------------------------------ 949 !scalars 950 integer,intent(in) :: idir,ipert,istwf_k,npw_k,npw1_k,natom,useylmgr1 951 integer,intent(out) :: nkpg,nkpg1 952 type(gs_hamiltonian_type),intent(inout) :: gs_hamkq 953 type(rf_hamiltonian_type),intent(inout) :: rf_hamkq 954 type(rf_hamiltonian_type),intent(inout),optional :: rf_hamk_dir2 955 type(dataset_type),intent(in) :: dtset 956 type(pseudopotential_type),intent(in) :: psps 957 !arrays 958 integer,intent(in) :: kg_k(3,npw_k),kg1_k(3,npw1_k) 959 real(dp),intent(in) :: kpoint(3),kpq(3),gmet(3,3),gprimd(3,3),rmet(3,3) 960 real(dp),intent(in) :: ylm_k(npw_k,psps%mpsang*psps%mpsang*psps%useylm) 961 real(dp),intent(in) :: ylmgr1_k(npw1_k,3+6*((ipert-natom)/10),psps%mpsang*psps%mpsang*psps%useylm*useylmgr1) 962 real(dp),intent(in) :: ylm1_k(npw1_k,psps%mpsang*psps%mpsang*psps%useylm) 963 real(dp),allocatable,intent(out) :: dkinpw(:),kinpw1(:),ffnlk(:,:,:,:),ffnl1(:,:,:,:) 964 real(dp),allocatable,intent(out),optional :: dkinpw2(:),ddkinpw(:) 965 real(dp),allocatable,intent(out) :: kpg_k(:,:),kpg1_k(:,:),ph3d(:,:,:),ph3d1(:,:,:) 966 967 !Local variables------------------------------- 968 !scalars 969 integer :: dimffnl1,dimffnlk,ider,idir0,idir1,idir2,istr,ntypat 970 logical :: qne0 971 !arrays 972 real(dp) :: ylmgr_dum(1,1,1) 973 974 ! ************************************************************************* 975 976 if(.not.present(ddkinpw) .and. ipert==natom+10) then 977 MSG_BUG("ddkinpw is not optional for ipert=natom+10.") 978 end if 979 if(.not.present(dkinpw2) .and. ipert==natom+10 .and. idir>3) then 980 MSG_BUG("dkinpw2 is not optional for ipert=natom+10 and idir>3.") 981 end if 982 if(.not.present(rf_hamk_dir2) .and. ((ipert==natom+10 .and. idir>3) .or. ipert==natom+11)) then 983 MSG_BUG("rf_hamk_dir2 is not optional for ipert=natom+10 (with idir>3) or ipert=natom+11.") 984 end if 985 986 ntypat = psps%ntypat 987 qne0=((kpq(1)-kpoint(1))**2+(kpq(2)-kpoint(2))**2+(kpq(3)-kpoint(3))**2>=tol14) 988 989 !Compute (k+G) vectors 990 nkpg=0;if(ipert>=1.and.ipert<=natom) nkpg=3*dtset%nloalg(3) 991 ABI_ALLOCATE(kpg_k,(npw_k,nkpg)) 992 if (nkpg>0) then 993 call mkkpg(kg_k,kpg_k,kpoint,nkpg,npw_k) 994 end if 995 996 !Compute (k+q+G) vectors 997 nkpg1=0;if(ipert>=1.and.ipert<=natom) nkpg1=3*dtset%nloalg(3) 998 ABI_ALLOCATE(kpg1_k,(npw1_k,nkpg1)) 999 if (nkpg1>0) then 1000 call mkkpg(kg1_k,kpg1_k,kpq(:),nkpg1,npw1_k) 1001 end if 1002 1003 !===== Preparation of the non-local contributions 1004 1005 dimffnlk=0;if (ipert<=natom) dimffnlk=1 1006 ABI_ALLOCATE(ffnlk,(npw_k,dimffnlk,psps%lmnmax,ntypat)) 1007 1008 !Compute nonlocal form factors ffnlk at (k+G) 1009 !(only for atomic displacement perturbation) 1010 if (ipert<=natom) then 1011 ider=0;idir0=0 1012 call mkffnl(psps%dimekb,dimffnlk,psps%ekb,ffnlk,psps%ffspl,& 1013 & gmet,gprimd,ider,idir0,psps%indlmn,kg_k,kpg_k,kpoint,psps%lmnmax,& 1014 & psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg,npw_k,ntypat,& 1015 & psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm_k,ylmgr_dum) 1016 end if 1017 1018 !Compute nonlocal form factors ffnl1 at (k+q+G) 1019 !-- Atomic displacement perturbation 1020 if (ipert<=natom) then 1021 ider=0;idir0=0 1022 !-- k-point perturbation (1st-derivative) 1023 else if (ipert==natom+1) then 1024 ider=1;idir0=idir 1025 !-- k-point perturbation (2nd-derivative) 1026 else if (ipert==natom+10.or.ipert==natom+11) then 1027 ider=2;idir0=4 1028 !-- Electric field perturbation 1029 else if (ipert==natom+2) then 1030 if (psps%usepaw==1) then 1031 ider=1;idir0=idir 1032 else 1033 ider=0;idir0=0 1034 end if 1035 !-- Strain perturbation 1036 else if (ipert==natom+3.or.ipert==natom+4) then 1037 if (ipert==natom+3) istr=idir 1038 if (ipert==natom+4) istr=idir+3 1039 ider=1;idir0=-istr 1040 !-- Magnetic field perturbation ( SPr, Zeeman ) 1041 else if(ipert==natom+5)then 1042 ider=0;idir0=0 1043 end if 1044 1045 !Compute nonlocal form factors ffnl1 at (k+q+G), for all atoms 1046 dimffnl1=1+ider 1047 if (ider==1.and.idir0==0) dimffnl1=2+2*psps%useylm 1048 if (ider==2.and.idir0==4) dimffnl1=3+7*psps%useylm 1049 ABI_ALLOCATE(ffnl1,(npw1_k,dimffnl1,psps%lmnmax,ntypat)) 1050 call mkffnl(psps%dimekb,dimffnl1,psps%ekb,ffnl1,psps%ffspl,gmet,gprimd,ider,idir0,& 1051 & psps%indlmn,kg1_k,kpg1_k,kpq,psps%lmnmax,psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg1,& 1052 & npw1_k,ntypat,psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm1_k,ylmgr1_k) 1053 1054 !===== Preparation of the kinetic contributions 1055 1056 !Note that not all these arrays should be allocated in the general case when wtk_k vanishes 1057 1058 !Compute (1/2) (2 Pi)**2 (k+q+G)**2: 1059 ABI_ALLOCATE(kinpw1,(npw1_k)) 1060 kinpw1(:)=zero 1061 call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg1_k,kinpw1,kpq,npw1_k,0,0) 1062 1063 ABI_ALLOCATE(dkinpw,(npw_k)) ! 1st derivative (1st direction) 1064 dkinpw(:)=zero 1065 if(ipert==natom+10 .and. idir>3) then 1066 ABI_ALLOCATE(dkinpw2,(npw_k)) ! 1st derivative (2nd directions) 1067 dkinpw2(:)=zero 1068 end if 1069 if(ipert==natom+10) then 1070 ABI_ALLOCATE(ddkinpw,(npw_k)) ! 2nd derivative 1071 ddkinpw(:)=zero 1072 end if 1073 1074 !-- k-point perturbation (1st-derivative) 1075 if (ipert==natom+1) then 1076 ! Compute the derivative of the kinetic operator vs k 1077 call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg_k,dkinpw,kpoint,npw_k,idir,0) ! 1st derivative 1078 end if 1079 1080 !-- k-point perturbation (2nd-derivative) 1081 if (ipert==natom+10.or.ipert==natom+11) then 1082 ! Compute the derivative of the kinetic operator vs k in kinpw, second and first orders 1083 if(ipert==natom+10 .and. idir<=3) then 1084 call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg_k,dkinpw,kpoint,npw_k,idir,0) ! 1st derivative 1085 call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg_k,ddkinpw,kpoint,npw_k,idir,idir) ! 2nd derivative 1086 else 1087 select case(idir) 1088 ! Diagonal terms : 1089 case(1) 1090 idir1 = 1 1091 idir2 = 1 1092 case(2) 1093 idir1 = 2 1094 idir2 = 2 1095 case(3) 1096 idir1 = 3 1097 idir2 = 3 1098 ! Upper triangular terms : 1099 case(4) 1100 idir1 = 2 1101 idir2 = 3 1102 case(5) 1103 idir1 = 1 1104 idir2 = 3 1105 case(6) 1106 idir1 = 1 1107 idir2 = 2 1108 ! Lower triangular terms : 1109 case(7) 1110 idir1 = 3 1111 idir2 = 2 1112 case(8) 1113 idir1 = 3 1114 idir2 = 1 1115 case(9) 1116 idir1 = 2 1117 idir2 = 1 1118 end select 1119 call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg_k,dkinpw,kpoint,npw_k,idir1,0) ! 1st derivative, idir1 1120 if(ipert==natom+10) then 1121 call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg_k,dkinpw2,kpoint,npw_k,idir2,0) ! 1st derivative, idir2 1122 call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg_k,ddkinpw,kpoint,npw_k,idir1,idir2) ! 2nd derivative 1123 end if 1124 end if 1125 end if 1126 1127 !-- Strain perturbation 1128 if (ipert==natom+3.or.ipert==natom+4) then 1129 if (ipert==natom+3) istr=idir 1130 if (ipert==natom+4) istr=idir+3 1131 ! Compute the derivative of the kinetic operator vs strain 1132 call kpgstr(dkinpw,dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,gprimd,istr,kg_k,kpoint,npw_k) 1133 end if 1134 1135 !===== Load the k/k+q dependent parts of the Hamiltonian 1136 1137 !Load k-dependent part in the Hamiltonian datastructure 1138 ABI_ALLOCATE(ph3d,(2,npw_k,gs_hamkq%matblk)) 1139 call load_k_hamiltonian(gs_hamkq,kpt_k=kpoint,npw_k=npw_k,istwf_k=istwf_k,kg_k=kg_k,kpg_k=kpg_k,& 1140 & ph3d_k=ph3d,compute_ph3d=.true.,compute_gbound=.true.) 1141 if (size(ffnlk)>0) then 1142 call load_k_hamiltonian(gs_hamkq,ffnl_k=ffnlk) 1143 else 1144 call load_k_hamiltonian(gs_hamkq,ffnl_k=ffnl1) 1145 end if 1146 1147 !Load k+q-dependent part in the Hamiltonian datastructure 1148 ! Note: istwf_k is imposed to 1 for RF calculations (should use istwf_kq instead) 1149 call load_kprime_hamiltonian(gs_hamkq,kpt_kp=kpq,npw_kp=npw1_k,istwf_kp=istwf_k,& 1150 & kinpw_kp=kinpw1,kg_kp=kg1_k,kpg_kp=kpg1_k,ffnl_kp=ffnl1,& 1151 & compute_gbound=.true.) 1152 if (qne0) then 1153 ABI_ALLOCATE(ph3d1,(2,npw1_k,gs_hamkq%matblk)) 1154 call load_kprime_hamiltonian(gs_hamkq,ph3d_kp=ph3d1,compute_ph3d=.true.) 1155 end if 1156 1157 !Load k-dependent part in the 1st-order Hamiltonian datastructure 1158 call load_k_rf_hamiltonian(rf_hamkq,npw_k=npw_k,dkinpw_k=dkinpw) 1159 if (ipert==natom+10) then 1160 call load_k_rf_hamiltonian(rf_hamkq,ddkinpw_k=ddkinpw) 1161 if (idir>3) then 1162 call load_k_rf_hamiltonian(rf_hamk_dir2,dkinpw_k=dkinpw2,ddkinpw_k=ddkinpw) 1163 end if 1164 end if 1165 1166 end subroutine getgh1c_setup
m_hamiltonian/rf_transgrid_and_pack [ Functions ]
[ Top ] [ m_hamiltonian ] [ Functions ]
NAME
rf_transgrid_and_pack
FUNCTION
Set up local potential vlocal1 with proper dimensioning, from vtrial1 taking into account the spin. Same thing for vlocal from vtrial.
INPUTS
isppol=Spin index. nspden=Number of density components usepaw=1 if PAW, 0 for NC. cplex=1 if DFPT potential is real, 2 for complex nfftf=Number of FFT points on the FINE grid treated by this processor nfft=Number of FFT points on the COARSE grid treated by this processor ngfft(18)=Info on the coarse grid. pawfgr <type(pawfgr_type)>=fine grid parameters and related data mpi_enreg=information about MPI parallelization vtrial(nfftf,nspden)=GS Vtrial(r) on the DENSE mesh vtrial1(cplex*nfftf,nspden)=INPUT RF Vtrial(r) on the DENSE mesh
OUTPUT
vlocal(n4,n5,n6,nvloc)= GS local potential in real space, on the augmented coarse fft grid vlocal1(cplex*n4,n5,n6,nvloc)= RF local potential in real space, on the augmented coarse fft grid
PARENTS
dfpt_vtorho,m_gkk,m_phgamma,m_phpi,m_sigmaph
CHILDREN
kpgstr,load_k_hamiltonian,load_k_rf_hamiltonian,load_kprime_hamiltonian mkffnl,mkkin,mkkpg
SOURCE
806 subroutine rf_transgrid_and_pack(isppol,nspden,usepaw,cplex,nfftf,nfft,ngfft,nvloc,& 807 & pawfgr,mpi_enreg,vtrial,vtrial1,vlocal,vlocal1) 808 809 use defs_basis 810 use defs_abitypes 811 use m_profiling_abi 812 use m_errors 813 814 use m_pawfgr, only : pawfgr_type 815 use m_kg, only : mkkin, kpgstr 816 817 !This section has been created automatically by the script Abilint (TD). 818 !Do not modify the following lines by hand. 819 #undef ABI_FUNC 820 #define ABI_FUNC 'rf_transgrid_and_pack' 821 use interfaces_53_ffts 822 use interfaces_65_paw 823 !End of the abilint section 824 825 implicit none 826 827 !Arguments ------------------------------------ 828 !scalars 829 integer,intent(in) :: isppol,nspden,usepaw,cplex,nfftf,nfft,nvloc 830 type(pawfgr_type),intent(in) :: pawfgr 831 type(MPI_type),intent(in) :: mpi_enreg 832 !arrays 833 integer,intent(in) :: ngfft(18) 834 real(dp),intent(in),target :: vtrial(nfftf,nspden) 835 real(dp),intent(inout),target :: vtrial1(cplex*nfftf,nspden) 836 real(dp),intent(out) :: vlocal(ngfft(4),ngfft(5),ngfft(6),nvloc) 837 real(dp),intent(out) :: vlocal1(cplex*ngfft(4),ngfft(5),ngfft(6),nvloc) 838 839 !Local variables------------------------------- 840 !scalars 841 integer :: n1,n2,n3,n4,n5,n6,paral_kgb,ispden 842 !arrays 843 real(dp) :: rhodum(1) 844 real(dp), ABI_CONTIGUOUS pointer :: vtrial_ptr(:,:),vtrial1_ptr(:,:) 845 real(dp),allocatable :: cgrvtrial(:,:),cgrvtrial1(:,:),vlocal_tmp(:,:,:),vlocal1_tmp(:,:,:) 846 847 ! ************************************************************************* 848 849 n1=ngfft(1); n2=ngfft(2); n3=ngfft(3) 850 n4=ngfft(4); n5=ngfft(5); n6=ngfft(6) 851 paral_kgb = mpi_enreg%paral_kgb 852 853 if (nspden/=4) then 854 vtrial_ptr=>vtrial 855 if (usepaw==0.or.pawfgr%usefinegrid==0) then 856 call fftpac(isppol,mpi_enreg,nspden,n1,n2,n3,n4,n5,n6,ngfft,vtrial_ptr,vlocal(:,:,:,1),2) 857 call fftpac(isppol,mpi_enreg,nspden,cplex*n1,n2,n3,cplex*n4,n5,n6,ngfft,vtrial1,vlocal1(:,:,:,1),2) 858 else 859 ABI_ALLOCATE(cgrvtrial,(nfft,nspden)) 860 call transgrid(1,mpi_enreg,nspden,-1,0,0,paral_kgb,pawfgr,rhodum,rhodum,cgrvtrial,vtrial_ptr) 861 call fftpac(isppol,mpi_enreg,nspden,n1,n2,n3,n4,n5,n6,ngfft,cgrvtrial,vlocal(:,:,:,1),2) 862 ABI_DEALLOCATE(cgrvtrial) 863 ABI_ALLOCATE(cgrvtrial,(cplex*nfft,nspden)) 864 call transgrid(cplex,mpi_enreg,nspden,-1,0,0,paral_kgb,pawfgr,rhodum,rhodum,cgrvtrial,vtrial1) 865 call fftpac(isppol,mpi_enreg,nspden,cplex*n1,n2,n3,cplex*n4,n5,n6,ngfft,cgrvtrial,vlocal1(:,:,:,1),2) 866 ABI_DEALLOCATE(cgrvtrial) 867 end if 868 nullify(vtrial_ptr) 869 else ! nspden==4 non-collinear magnetism 870 vtrial_ptr=>vtrial 871 vtrial1_ptr=>vtrial1 872 ABI_ALLOCATE(vlocal_tmp,(n4,n5,n6)) 873 ABI_ALLOCATE(vlocal1_tmp,(cplex*n4,n5,n6)) 874 if (usepaw==0.or.pawfgr%usefinegrid==0) then 875 do ispden=1,nspden 876 call fftpac(ispden,mpi_enreg,nspden,n1,n2,n3,n4,n5,n6,ngfft,vtrial_ptr,vlocal_tmp,2) 877 vlocal(:,:,:,ispden)=vlocal_tmp(:,:,:) 878 call fftpac(ispden,mpi_enreg,nspden,cplex*n1,n2,n3,cplex*n4,n5,n6,ngfft,vtrial1_ptr,vlocal1_tmp,2) 879 vlocal1(:,:,:,ispden)=vlocal1_tmp(:,:,:) 880 end do 881 else ! TODO FR EB check the correctness of the following lines for PAW calculations 882 ABI_ALLOCATE(cgrvtrial,(nfft,nspden)) 883 ABI_ALLOCATE(cgrvtrial1,(nfft,nspden)) 884 call transgrid(cplex,mpi_enreg,nspden,-1,0,0,paral_kgb,pawfgr,rhodum,rhodum,cgrvtrial,vtrial_ptr) 885 call transgrid(cplex,mpi_enreg,nspden,-1,0,0,paral_kgb,pawfgr,rhodum,rhodum,cgrvtrial1,vtrial1_ptr) 886 do ispden=1,nspden 887 call fftpac(ispden,mpi_enreg,nspden,n1,n2,n3,n4,n5,n6,ngfft,vtrial_ptr,vlocal_tmp,2) 888 vlocal(:,:,:,ispden)=vlocal_tmp(:,:,:) 889 call fftpac(ispden,mpi_enreg,nspden,n1,n2,n3,n4,n5,n6,ngfft,vtrial1_ptr,vlocal1_tmp,2) 890 vlocal1(:,:,:,ispden)=vlocal1_tmp(:,:,:) 891 end do 892 ABI_DEALLOCATE(cgrvtrial) 893 end if 894 ABI_DEALLOCATE(vlocal_tmp) 895 ABI_DEALLOCATE(vlocal1_tmp) 896 end if !nspden 897 898 end subroutine rf_transgrid_and_pack