TABLE OF CONTENTS
- ABINIT/getdc1
- ABINIT/getgh1c
- ABINIT/getgh1c_mGGA
- ABINIT/getgh1dqc
- ABINIT/getgh1ndc
- ABINIT/m_getgh1c
- m_hamiltonian/getgh1c_setup
- m_hamiltonian/getgh1dqc_setup
- m_hamiltonian/rf_transgrid_and_pack
ABINIT/getdc1 [ Functions ]
NAME
getdc1
FUNCTION
Compute |delta_C^(1)> from one wave function C - PAW ONLY
Compute <G|delta_C^(1)> (dcwavef) and eventually <P_i| delta_C^(1)> (dcwaveprj) where P_i= non-local projector
delta_C^(1) is the variation of wavefunction only due to variation of overlap operator S.
delta_C^(1)=-1/2.Sum_j [ <C_j|S^(1)|C>.C_j
see PRB 78, 035105 (2008) [[cite:Audouze2008]], Eq. (42 and 40, term 2)
INPUTS
cgq(2,mcgq)=wavefunction coefficients for all bands j on present processor, at k+Q: cgq=< G |Cnk+q> cprjq(natom,mcprjq)= wave functions j at k+q projected with non-local projectors: cprjq=<P_i|Cnk+q> ibgq=shift to be applied on the location of data in the array cprjq icgq=shift to be applied on the location of data in the array cgq istwfk=option parameter that describes the storage of wfs mcgq=second dimension of the cgq array mcprjq=second dimension of the cprjq array mpi_enreg=information about MPI parallelization natom= number of atoms in cell nband=number of bands npw1=number of planewaves in basis sphere at k+Q nspinor=number of spinorial components of the wavefunctions opt_cprj=flag governing the computation of <P_i|delta_C^(1)> (P_i= non-local projector) s1cwave0(2,npw1*nspinor)=<G|S^(1)|C> where S^(1) is the first-order overlap operator
OUTPUT
dcwavef(2,npw1*nspinor)=change of wavefunction due to change of overlap PROJECTED ON PLANE-WAVES:
dcwavef is delta_C(1)=-1/2.Sum_{j}[<C0_k+q_j|S(1)|C0_k_i>.|C0_k+q_j>]
=== if optcprj=1 ===
dcwaveprj(natom,nspinor*optcprj)=change of wavefunction due to change of overlap PROJECTED ON NL-PROJECTORS:
SOURCE
1274 subroutine getdc1(band,band_procs,bands_treated_now,cgq,cprjq,dcwavef,dcwaveprj,& 1275 & ibgq,icgq,istwfk,mcgq,mcprjq,& 1276 & mpi_enreg,natom,nband,nband_me,npw1,nspinor,optcprj,s1cwave0) 1277 1278 !Arguments ------------------------------------ 1279 !scalars 1280 integer,intent(in) :: ibgq,icgq,istwfk,mcgq,mcprjq,natom,nband,npw1,nspinor,optcprj 1281 integer,intent(in) :: band, nband_me 1282 type(MPI_type),intent(in) :: mpi_enreg 1283 !arrays 1284 integer,intent(in) :: band_procs(nband),bands_treated_now(nband) 1285 real(dp),intent(in) :: cgq(2,mcgq),s1cwave0(2,npw1*nspinor) 1286 real(dp),intent(out) :: dcwavef(2,npw1*nspinor) 1287 type(pawcprj_type),intent(in) :: cprjq(natom,mcprjq) 1288 type(pawcprj_type),intent(inout) :: dcwaveprj(natom,nspinor*optcprj) 1289 1290 !Local variables------------------------------- 1291 !scalars 1292 integer, parameter :: tim_projbd=0 1293 integer :: ipw 1294 integer :: band_, ierr 1295 real(dp),parameter :: scal=-half 1296 !arrays 1297 real(dp), allocatable :: dummy(:,:),scprod(:,:) 1298 real(dp), allocatable :: dcwavef_tmp(:,:) 1299 type(pawcprj_type),allocatable :: dcwaveprj_tmp(:,:) 1300 1301 ! ********************************************************************* 1302 1303 DBG_ENTER("COLL") 1304 1305 ABI_MALLOC(dummy,(0,0)) 1306 ABI_MALLOC(scprod,(2,nband_me)) 1307 ABI_MALLOC(dcwavef_tmp,(2,npw1*nspinor)) 1308 if (optcprj == 1) then 1309 ABI_MALLOC(dcwaveprj_tmp,(natom,nspinor*optcprj)) 1310 call pawcprj_alloc(dcwaveprj_tmp, 0, dcwaveprj(:,1)%nlmn) 1311 end if 1312 1313 !=== 1- COMPUTE: <G|S^(1)|C_k> - Sum_j [<C_k+q,j|S^(1)|C_k>.<G|C_k+q,j>] 1314 !! using the projb routine 1315 !Note the subtlety: projbd is called with useoverlap=0 and s1cwave0 1316 !in order to get Sum[<cgq|s1|c>|cgq>]=Sum[<cgq|gs1>|cgq>] 1317 1318 ! run over procs in my pool which have a dcwavef to projbd 1319 do band_ = 1, nband 1320 if (bands_treated_now(band_) == 0) cycle 1321 dcwavef_tmp = zero 1322 1323 ! distribute dcwavef_tmp to my band pool 1324 ! everyone works on a single band s1cwave0 = <G|S^(1)|C_k> 1325 if (band_ == band) then 1326 !$OMP PARALLEL DO 1327 do ipw=1,npw1*nspinor 1328 dcwavef_tmp(1:2,ipw)=s1cwave0(1:2,ipw) 1329 end do 1330 end if 1331 call xmpi_bcast(dcwavef_tmp,band_procs(band_),mpi_enreg%comm_band,ierr) 1332 1333 ! get the projbd onto my processor's bands dcwavef = dcwavef - <cgq|dcwavef>|cgq> 1334 ! dcwavef = <G|S^(1)|C_k> - Sum_{MYj} [<C_k+q,j|S^(1)|C_k>.<G|C_k+q,j>] 1335 ! scprod = <C_k+q,j|S^(1)|C_k> for {MYj} 1336 call projbd(cgq,dcwavef_tmp,-1,icgq,0,istwfk,mcgq,0,nband_me,npw1,nspinor,& 1337 & dummy,scprod,0,tim_projbd,0,mpi_enreg%me_g0,mpi_enreg%comm_fft) 1338 1339 1340 ! sum all of the corrections 1341 ! dcwavef = Nprocband * <G|S^(1)|C_k> - Sum_{ALLj} [<C_k+q,j|S^(1)|C_k>.<G|C_k+q,j>] 1342 call xmpi_sum(dcwavef_tmp,mpi_enreg%comm_band,ierr) 1343 1344 ! save to my proc if it is my turn, and subtract Ntuple counted dcwavef 1345 if (band_ == band) then 1346 !=== 2- COMPUTE: <G|delta_C^(1)> = -1/2.Sum_j [<C_k+q,j|S^(1)|C_k>.<G|C_k+q,j>] by substraction 1347 ! tested this is equivalent to previous coding to within 1.e-18 accumulated error (probably in favor of this coding) 1348 !$OMP PARALLEL DO PRIVATE(ipw) SHARED(dcwavef,s1cwave0,dcwavef_tmp,npw1,nspinor) 1349 do ipw=1,npw1*nspinor 1350 dcwavef(1:2,ipw)= scal*(mpi_enreg%nproc_band*s1cwave0(1:2,ipw)-dcwavef_tmp(1:2,ipw)) 1351 end do 1352 end if 1353 1354 !=== 3- COMPUTE: <P_i|delta_C^(1)> = -1/2.Sum_j [<C_k+q,j|S^(1)|C_k>.<P_i|C_k+q,j>] 1355 ! as above everyone has to operate on each band band_ 1356 if (optcprj==1.and.mcprjq>0) then 1357 ! cprjq = <P_i|C_k+q,j> for MYj 1358 ! dcwaveprj_tmp = Sum_MYj [<C_k+q,j|S^(1)|C_k>.<P_i|C_k+q,j>] 1359 call pawcprj_lincom(scprod,cprjq(:,ibgq+1:ibgq+nspinor*nband_me),dcwaveprj_tmp,nband_me) 1360 1361 ! still need to mpisum the dcwaveprj to get linear combination of all bands, not just mine 1362 ! dcwaveprj = Sum_ALLj [<C_k+q,j|S^(1)|C_k,i>.<P_i|C_k+q,j>] 1363 call pawcprj_mpi_sum(dcwaveprj_tmp,mpi_enreg%comm_band,ierr) 1364 1365 if (band_ == band) then 1366 ! dcwaveprj = -1/2 dcwaveprj_tmp 1367 !TODO: check the correct order of scal and zero (alpha / beta) coefficients. 1368 ! Here dcwaveprj is squashed by the _tmp variable which is used in parallel 1369 call pawcprj_axpby(scal,zero,dcwaveprj_tmp,dcwaveprj) 1370 end if 1371 end if 1372 1373 end do ! procs in my band pool 1374 1375 ABI_FREE(dummy) 1376 ABI_FREE(scprod) 1377 ABI_FREE(dcwavef_tmp) 1378 if (optcprj == 1) then 1379 call pawcprj_free(dcwaveprj_tmp) 1380 ABI_FREE(dcwaveprj_tmp) 1381 end if 1382 1383 DBG_EXIT("COLL") 1384 1385 end subroutine getdc1
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 gvnlx1c. If required, <G|S^(1)|C> is returned in gs1c (S=overlap - PAW only)
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 sorted) 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_gvnlx1=option controlling the use of gvnlx1 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 gvnlx1=(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) 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) if (sij_opt=1) gs1c(2,npw1*nspinor)=<G|S^(1)|C> (S=overlap) on the k+q sphere.
SOURCE
119 subroutine getgh1c(berryopt,cwave,cwaveprj,gh1c,grad_berry,gs1c,gs_hamkq,& 120 gvnlx1,idir,ipert,lambda,mpi_enreg,optlocal,optnl,opt_gvnlx1,& 121 rf_hamkq,sij_opt,tim_getgh1c,usevnl,conj) 122 123 !Arguments ------------------------------------ 124 !scalars 125 logical,intent(in),optional :: conj 126 integer,intent(in) :: berryopt,idir,ipert,optlocal,optnl,opt_gvnlx1,sij_opt,tim_getgh1c,usevnl 127 real(dp),intent(in) :: lambda 128 type(MPI_type),intent(in) :: mpi_enreg 129 type(gs_hamiltonian_type),intent(inout),target :: gs_hamkq 130 type(rf_hamiltonian_type),intent(inout),target :: rf_hamkq 131 !arrays 132 real(dp),intent(in) :: grad_berry(:,:) 133 real(dp),intent(inout) :: cwave(2,gs_hamkq%npw_k*gs_hamkq%nspinor) 134 real(dp),intent(out) :: gh1c(2,gs_hamkq%npw_kp*gs_hamkq%nspinor) 135 real(dp),intent(out) :: gs1c(2,gs_hamkq%npw_kp*gs_hamkq%nspinor) 136 real(dp),intent(inout),target :: gvnlx1(2,gs_hamkq%npw_kp*gs_hamkq%nspinor) 137 type(pawcprj_type),intent(inout),target :: cwaveprj(:,:) 138 139 !Local variables------------------------------- 140 !scalars 141 integer,parameter :: level=16 142 integer :: choice,cplex1,cpopt,ipw,ipws,ispinor,istr,i1,i2,i3 143 integer :: my_nspinor,natom,ncpgr,nnlout=1,npw,npw1,paw_opt,signs 144 integer :: tim_fourwf,tim_nonlop,usecprj 145 logical :: compute_conjugate,has_kin,has_mGGA1,has_nd1,usevnl2 146 real(dp) :: weight !, cpu, wall, gflops 147 !character(len=500) :: msg 148 !arrays 149 real(dp) :: enlout(1),tsec(2),svectout_dum(1,1),vectout_dum(1,1) 150 real(dp),allocatable :: cwave_sp(:,:),cwavef1(:,:),cwavef2(:,:) 151 real(dp),allocatable :: gh1c_sp(:,:),gh1c1(:,:),gh1c2(:,:),gh1c3(:,:),gh1c4(:,:) 152 real(dp),allocatable :: gh1c_mGGA(:,:),gh1ndc(:,:),gvnl2(:,:) 153 real(dp),allocatable :: nonlop_out(:,:),vlocal1_tmp(:,:,:),work(:,:,:,:) 154 real(dp),ABI_CONTIGUOUS pointer :: gvnlx1_(:,:) 155 real(dp),pointer :: dkinpw(:),kinpw1(:) 156 type(pawcprj_type),allocatable,target :: cwaveprj_tmp(:,:) 157 type(pawcprj_type),pointer :: cwaveprj_ptr(:,:) 158 159 ! ********************************************************************* 160 161 DBG_ENTER("COLL") 162 163 ! Keep track of total time spent in getgh1c 164 call timab(196+tim_getgh1c,1,tsec) 165 166 !====================================================================== 167 !== Initialisations and compatibility tests 168 !====================================================================== 169 170 npw =gs_hamkq%npw_k 171 npw1 =gs_hamkq%npw_kp 172 natom=gs_hamkq%natom 173 174 ! Compatibility tests 175 if(gs_hamkq%usepaw==1.and.(ipert>=0.and.(ipert<=natom.or.ipert==natom+3.or.ipert==natom+4))) then 176 if ((optnl>=1.and.(.not.associated(rf_hamkq%e1kbfr))) .or. & 177 (optnl==2.and.(.not.associated(rf_hamkq%e1kbsc)))) then 178 ABI_BUG('ekb derivatives must be allocated for ipert<=natom or natom+3/4 !') 179 end if 180 end if 181 if(gs_hamkq%usepaw==1.and.(ipert==natom+2)) then 182 if ((optnl>=1.and.(.not.associated(rf_hamkq%e1kbfr))) .or. & 183 (optnl==2.and.(.not.associated(rf_hamkq%e1kbsc)))) then 184 ABI_BUG('ekb derivatives must be allocated for ipert=natom+2 !') 185 end if 186 if (usevnl==0) then 187 ABI_BUG('gvnlx1 must be allocated for ipert=natom+2 !') 188 end if 189 end if 190 if(ipert==natom+2.and.opt_gvnlx1==0) then 191 ABI_BUG('opt_gvnlx1=0 not compatible with ipert=natom+2 !') 192 end if 193 if (mpi_enreg%paral_spinor==1) then 194 ABI_BUG('Not compatible with parallelization over spinorial components !') 195 end if 196 197 ! Check sizes 198 my_nspinor=max(1,gs_hamkq%nspinor/mpi_enreg%nproc_spinor) 199 if (size(cwave)<2*npw*my_nspinor) then 200 ABI_BUG('wrong size for cwave!') 201 end if 202 if (size(gh1c)<2*npw1*my_nspinor) then 203 ABI_BUG('wrong size for gh1c!') 204 end if 205 if (usevnl/=0) then 206 if (size(gvnlx1)<2*npw1*my_nspinor) then 207 ABI_BUG('wrong size for gvnlx1!') 208 end if 209 end if 210 if (sij_opt==1) then 211 if (size(gs1c)<2*npw1*my_nspinor) then 212 ABI_BUG('wrong size for gs1c!') 213 end if 214 end if 215 if (berryopt>=4) then 216 if (size(grad_berry)<2*npw1*my_nspinor) then 217 ABI_BUG('wrong size for grad_berry!') 218 end if 219 end if 220 221 ! PAW: specific treatment for usecprj input arg 222 ! force it to zero if cwaveprj is not allocated 223 usecprj=gs_hamkq%usecprj ; ncpgr=0 224 if(gs_hamkq%usepaw==1) then 225 if (size(cwaveprj)==0) usecprj=0 226 if (usecprj/=0) then 227 ncpgr=cwaveprj(1,1)%ncpgr 228 if (size(cwaveprj)<gs_hamkq%natom*my_nspinor) then 229 ABI_BUG('wrong size for cwaveprj!') 230 end if 231 if(gs_hamkq%usepaw==1.and.(ipert>=0.and.(ipert<=natom.or.ipert==natom+3.or.ipert==natom+4))) then 232 if (ncpgr/=1)then 233 ABI_BUG('Projected WFs (cprj) derivatives are not correctly stored !') 234 end if 235 end if 236 end if 237 else 238 if(usecprj==1)then 239 ABI_BUG('usecprj==1 not allowed for NC psps !') 240 end if 241 end if 242 243 tim_nonlop=8 244 if (tim_getgh1c==1.and.ipert<=natom) tim_nonlop=7 245 if (tim_getgh1c==2.and.ipert<=natom) tim_nonlop=5 246 if (tim_getgh1c==1.and.ipert> natom) tim_nonlop=8 247 if (tim_getgh1c==2.and.ipert> natom) tim_nonlop=5 248 if (tim_getgh1c==3 ) tim_nonlop=0 249 250 compute_conjugate = .false. 251 if(present(conj)) compute_conjugate = conj 252 253 !====================================================================== 254 !== Apply the 1st-order local potential to the wavefunction 255 !====================================================================== 256 257 !Phonon perturbation 258 !or Electric field perturbation 259 !or Strain perturbation 260 !------------------------------------------- 261 if (ipert<=natom+5.and.ipert/=natom+1.and.optlocal>0) then !SPr deb 262 263 ABI_MALLOC(work,(2,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6)) 264 265 if (gs_hamkq%nvloc==1) then 266 267 weight=one ; tim_fourwf=4 268 call fourwf(rf_hamkq%cplex,rf_hamkq%vlocal1,cwave,gh1c,work,gs_hamkq%gbound_k,gs_hamkq%gbound_kp,& 269 gs_hamkq%istwf_k,gs_hamkq%kg_k,gs_hamkq%kg_kp,gs_hamkq%mgfft,mpi_enreg,1,gs_hamkq%ngfft,& 270 npw,npw1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,2,tim_fourwf,weight,weight,& 271 gpu_option=gs_hamkq%gpu_option) 272 if(gs_hamkq%nspinor==2)then 273 ABI_MALLOC(cwave_sp,(2,npw)) 274 ABI_MALLOC(gh1c_sp,(2,npw1)) 275 !$OMP PARALLEL DO 276 do ipw=1,npw 277 cwave_sp(1,ipw)=cwave(1,ipw+npw) 278 cwave_sp(2,ipw)=cwave(2,ipw+npw) 279 end do 280 call fourwf(rf_hamkq%cplex,rf_hamkq%vlocal1,cwave_sp,gh1c_sp,work,gs_hamkq%gbound_k,gs_hamkq%gbound_kp,& 281 gs_hamkq%istwf_k,gs_hamkq%kg_k,gs_hamkq%kg_kp,gs_hamkq%mgfft,mpi_enreg,1,gs_hamkq%ngfft,& 282 npw,npw1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,2,tim_fourwf,weight,weight,& 283 gpu_option=gs_hamkq%gpu_option) 284 !$OMP PARALLEL DO 285 do ipw=1,npw1 286 gh1c(1,ipw+npw1)=gh1c_sp(1,ipw) 287 gh1c(2,ipw+npw1)=gh1c_sp(2,ipw) 288 end do 289 ABI_FREE(cwave_sp) 290 ABI_FREE(gh1c_sp) 291 end if 292 else 293 ! Non-Collinear magnetism for nvloc=4 294 if (gs_hamkq%nspinor==2) then 295 weight=one ; tim_fourwf=4 296 ABI_MALLOC(gh1c1,(2,npw1)) 297 ABI_MALLOC(gh1c2,(2,npw1)) 298 ABI_MALLOC(gh1c3,(2,npw1)) 299 ABI_MALLOC(gh1c4,(2,npw1)) 300 gh1c1(:,:)=zero; gh1c2(:,:)=zero; gh1c3(:,:)=zero ; gh1c4(:,:)=zero 301 ABI_MALLOC(vlocal1_tmp,(rf_hamkq%cplex*gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6)) 302 !SPr: notation/dimension corrected vlocal_tmp -> vlocal1_tmp 303 ABI_MALLOC(cwavef1,(2,npw)) 304 ABI_MALLOC(cwavef2,(2,npw)) 305 do ipw=1,npw 306 cwavef1(1:2,ipw)=cwave(1:2,ipw) 307 cwavef2(1:2,ipw)=cwave(1:2,ipw+npw) 308 end do 309 ! gh1c1=v11*phi1 310 vlocal1_tmp(:,:,:)=rf_hamkq%vlocal1(:,:,:,1) 311 call fourwf(rf_hamkq%cplex,vlocal1_tmp,cwavef1,gh1c1,work,gs_hamkq%gbound_k,gs_hamkq%gbound_kp,& 312 gs_hamkq%istwf_k,gs_hamkq%kg_k,gs_hamkq%kg_kp,gs_hamkq%mgfft,mpi_enreg,1,gs_hamkq%ngfft,& 313 npw,npw1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,2,tim_fourwf,weight,weight,& 314 gpu_option=gs_hamkq%gpu_option) 315 ! gh1c2=v22*phi2 316 vlocal1_tmp(:,:,:)=rf_hamkq%vlocal1(:,:,:,2) 317 call fourwf(rf_hamkq%cplex,vlocal1_tmp,cwavef2,gh1c2,work,gs_hamkq%gbound_k,gs_hamkq%gbound_kp,& 318 gs_hamkq%istwf_k,gs_hamkq%kg_k,gs_hamkq%kg_kp,gs_hamkq%mgfft,mpi_enreg,1,gs_hamkq%ngfft,& 319 npw,npw1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,2,tim_fourwf,weight,weight,& 320 gpu_option=gs_hamkq%gpu_option) 321 ABI_FREE(vlocal1_tmp) 322 cplex1=2 323 ABI_MALLOC(vlocal1_tmp,(cplex1*gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6)) 324 ! gh1c3=(re(v12)-im(v12))*phi1 => v^21*phi1 325 if(rf_hamkq%cplex==1) then 326 do i3=1,gs_hamkq%n6 327 do i2=1,gs_hamkq%n5 328 do i1=1,gs_hamkq%n4 329 vlocal1_tmp(2*i1-1,i2,i3)= rf_hamkq%vlocal1(i1,i2,i3,3) 330 vlocal1_tmp(2*i1 ,i2,i3)=-rf_hamkq%vlocal1(i1,i2,i3,4) 331 end do 332 end do 333 end do 334 else 335 !SPr: modified definition of local potential components for cplex=2 (see dotprod_vn) 336 !also, v21==v12* not always holds (e.g. magnetic field perturbation) 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-1,i2,i3)= rf_hamkq%vlocal1(2*i1 ,i2,i3,4) 341 vlocal1_tmp(2*i1 ,i2,i3)=-rf_hamkq%vlocal1(2*i1-1,i2,i3,4) 342 end do 343 end do 344 end do 345 end if 346 call fourwf(cplex1,vlocal1_tmp,cwavef1,gh1c3,work,gs_hamkq%gbound_k,gs_hamkq%gbound_kp,& 347 gs_hamkq%istwf_k,gs_hamkq%kg_k,gs_hamkq%kg_kp,gs_hamkq%mgfft,mpi_enreg,1,gs_hamkq%ngfft,& 348 npw,npw1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,2,tim_fourwf,weight,weight,& 349 gpu_option=gs_hamkq%gpu_option) 350 ! gh1c4=(re(v12)+im(v12))*phi2 => v^12*phi2 351 if(rf_hamkq%cplex==1) then 352 do i3=1,gs_hamkq%n6 353 do i2=1,gs_hamkq%n5 354 do i1=1,gs_hamkq%n4 355 vlocal1_tmp(2*i1,i2,i3)=-vlocal1_tmp(2*i1,i2,i3) 356 end do 357 end do 358 end do 359 else 360 !for cplex=2 and time-reversal breaking perturbations,v21/=v12* 361 do i3=1,gs_hamkq%n6 362 do i2=1,gs_hamkq%n5 363 do i1=1,gs_hamkq%n4 364 vlocal1_tmp(2*i1-1,i2,i3)= rf_hamkq%vlocal1(2*i1-1,i2,i3,3) 365 vlocal1_tmp(2*i1 ,i2,i3)= rf_hamkq%vlocal1(2*i1 ,i2,i3,3) 366 end do 367 end do 368 end do 369 end if 370 call fourwf(cplex1,vlocal1_tmp,cwavef2,gh1c4,work,gs_hamkq%gbound_k,gs_hamkq%gbound_kp,& 371 gs_hamkq%istwf_k,gs_hamkq%kg_k,gs_hamkq%kg_kp,gs_hamkq%mgfft,mpi_enreg,1,gs_hamkq%ngfft,& 372 npw,npw1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,2,tim_fourwf,weight,weight,& 373 gpu_option=gs_hamkq%gpu_option) 374 ABI_FREE(vlocal1_tmp) 375 ! Build gh1c from pieces 376 ! gh1c_1 = (v11, v12) (psi1) matrix vector product 377 ! gh1c_2 = (v12*,v22) (psi2) 378 do ipw=1,npw1 379 gh1c(1:2,ipw) =gh1c1(1:2,ipw)+gh1c4(1:2,ipw) 380 gh1c(1:2,ipw+npw1)=gh1c3(1:2,ipw)+gh1c2(1:2,ipw) 381 end do 382 ABI_FREE(gh1c1) 383 ABI_FREE(gh1c2) 384 ABI_FREE(gh1c3) 385 ABI_FREE(gh1c4) 386 ABI_FREE(cwavef1) 387 ABI_FREE(cwavef2) 388 else 389 ABI_BUG('nspinor/=1 for Non-collinear calculations!') 390 end if 391 end if ! nvloc 392 393 ABI_FREE(work) 394 395 ! k-point perturbation (or no local part, i.e. optlocal=0) 396 ! ------------------------------------------- 397 else if (ipert==natom+1.or.optlocal==0) then 398 399 ! In the case of ddk operator, no local contribution (also because no self-consistency) 400 !$OMP PARALLEL DO 401 do ipw=1,npw1*my_nspinor 402 gh1c(:,ipw)=zero 403 end do 404 405 end if 406 407 !====================================================================== 408 !== Apply the 1st-order non-local potential to the wavefunction 409 !====================================================================== 410 411 !Use of gvnlx1 depends on usevnl 412 if (usevnl==1) then 413 gvnlx1_ => gvnlx1 414 else 415 ABI_MALLOC(gvnlx1_,(2,npw1*my_nspinor)) 416 end if 417 418 !Phonon perturbation 419 !------------------------------------------- 420 if (ipert<=natom.and.(optnl>0.or.sij_opt/=0)) then 421 422 ! PAW: 423 if (gs_hamkq%usepaw==1) then 424 425 if (usecprj==1) then 426 cwaveprj_ptr => cwaveprj 427 else 428 ABI_MALLOC(cwaveprj_tmp,(natom,my_nspinor)) 429 call pawcprj_alloc(cwaveprj_tmp,1,gs_hamkq%dimcprj) 430 cwaveprj_ptr => cwaveprj_tmp 431 end if 432 433 ! 1- Compute derivatives due to projectors |p_i>^(1) 434 ! Only displaced atom contributes 435 cpopt=-1+5*usecprj ; choice=2 ; signs=2 436 paw_opt=1;if (sij_opt/=0) paw_opt=sij_opt+3 437 call nonlop(choice,cpopt,cwaveprj_ptr,enlout,gs_hamkq,idir,(/lambda/),mpi_enreg,1,nnlout,& 438 & paw_opt,signs,gs1c,tim_nonlop,cwave,gvnlx1_,iatom_only=ipert) 439 440 ! 2- Compute derivatives due to frozen part of D_ij^(1) (independent of VHxc^(1)) 441 ! All atoms contribute 442 if (optnl>=1) then 443 ABI_MALLOC(nonlop_out,(2,npw1*my_nspinor)) 444 cpopt=1+3*usecprj ; choice=1 ; signs=2 ; paw_opt=1 445 call nonlop(choice,cpopt,cwaveprj_ptr,enlout,gs_hamkq,idir,(/lambda/),mpi_enreg,1,nnlout,& 446 & paw_opt,signs,svectout_dum,tim_nonlop,cwave,nonlop_out,enl=rf_hamkq%e1kbfr) 447 !$OMP PARALLEL DO 448 do ipw=1,npw1*my_nspinor 449 gvnlx1_(:,ipw)=gvnlx1_(:,ipw)+nonlop_out(:,ipw) 450 end do 451 ABI_FREE(nonlop_out) 452 end if 453 454 ! 3- Compute derivatives due to self-consistent part of D_ij^(1) (depending on VHxc^(1)) 455 ! All atoms contribute 456 if (optnl==2) then 457 ABI_MALLOC(gvnl2,(2,npw1*my_nspinor)) 458 cpopt=4 ; choice=1 ; signs=2 ; paw_opt=1 459 call nonlop(choice,cpopt,cwaveprj_ptr,enlout,gs_hamkq,idir,(/lambda/),mpi_enreg,1,nnlout,& 460 & paw_opt,signs,svectout_dum,tim_nonlop,cwave,gvnl2,enl=rf_hamkq%e1kbsc) 461 end if 462 463 if (usecprj==0) then 464 call pawcprj_free(cwaveprj_tmp) 465 ABI_FREE(cwaveprj_tmp) 466 end if 467 nullify(cwaveprj_ptr) 468 469 ! Norm-conserving psps: 470 else 471 ! Compute only derivatives due to projectors |p_i>^(1) 472 cpopt=-1 ; choice=2 ; signs=2 ; paw_opt=0 473 call nonlop(choice,cpopt,cwaveprj,enlout,gs_hamkq,idir,(/lambda/),mpi_enreg,1,nnlout,& 474 & paw_opt,signs,svectout_dum,tim_nonlop,cwave,gvnlx1_,iatom_only=ipert) 475 if (sij_opt==1) then 476 !$OMP PARALLEL DO 477 do ipw=1,npw1*my_nspinor 478 gs1c(:,ipw)=zero 479 end do 480 end if 481 end if 482 483 ! k-point perturbation 484 ! ------------------------------------------- 485 else if (ipert==natom+1.and.(optnl>0.or.sij_opt/=0)) then 486 487 tim_nonlop=8 ; signs=2 ; choice=5 488 if (gs_hamkq%usepaw==1) then 489 if (usecprj==1) then 490 cwaveprj_ptr => cwaveprj 491 else 492 ABI_MALLOC(cwaveprj_tmp,(natom,my_nspinor)) 493 call pawcprj_alloc(cwaveprj_tmp,1,gs_hamkq%dimcprj) 494 cwaveprj_ptr => cwaveprj_tmp 495 end if 496 cpopt=-1+5*usecprj; paw_opt=1; if (sij_opt/=0) paw_opt=sij_opt+3 497 ! JLJ: BUG (wrong result) of H^(1) if stored cprj are used in PAW DDKs with nspinor==2 (==1 works fine). 498 ! To be debugged, if someone has time... 499 if(gs_hamkq%nspinor==2) cpopt=-1 500 if(associated(gs_hamkq%vectornd)) cpopt=-1 501 call nonlop(choice,cpopt,cwaveprj_ptr,enlout,gs_hamkq,idir,(/lambda/),mpi_enreg,1,nnlout,& 502 & paw_opt,signs,gs1c,tim_nonlop,cwave,gvnlx1_) 503 if (usecprj==0) then 504 call pawcprj_free(cwaveprj_tmp) 505 ABI_FREE(cwaveprj_tmp) 506 end if 507 nullify(cwaveprj_ptr) 508 else 509 cpopt=-1 ; paw_opt=0 510 call nonlop(choice,cpopt,cwaveprj,enlout,gs_hamkq,idir,(/lambda/),mpi_enreg,1,nnlout,& 511 & paw_opt,signs,svectout_dum,tim_nonlop,cwave,gvnlx1_) 512 end if 513 !DEBUG 514 ! gvnlx1_(:,:)=zero 515 !ENDDEBUG 516 517 ! Electric field perturbation without Berry phase 518 ! ------------------------------------------- 519 else if(ipert==natom+2 .and. & 520 & (berryopt/=4 .and. berryopt/=6 .and. berryopt/=7 .and. & 521 & berryopt/=14 .and. berryopt/=16 .and. berryopt/=17) .and.(optnl>0.or.sij_opt/=0))then 522 ! gvnlx1 was already initialized in the calling routine, by reading a ddk file 523 ! It contains |i du^(0)/dk_band> 524 525 if (gs_hamkq%usepaw==1) then 526 if (usecprj==1) then 527 cwaveprj_ptr => cwaveprj 528 else 529 ABI_MALLOC(cwaveprj_tmp,(natom,my_nspinor)) 530 call pawcprj_alloc(cwaveprj_tmp,1,gs_hamkq%dimcprj) 531 cwaveprj_ptr => cwaveprj_tmp 532 end if 533 if (opt_gvnlx1==2.and.optnl>=1) then 534 535 ! PAW: Compute application of S^(0) to ddk WF 536 cpopt=-1 ; choice=1 ; paw_opt=3 ; signs=2 537 ABI_MALLOC(nonlop_out,(2,npw1*my_nspinor)) 538 call nonlop(choice,cpopt,cwaveprj_ptr,enlout,gs_hamkq,0,(/lambda/),mpi_enreg,1,nnlout,& 539 & paw_opt,signs,nonlop_out,tim_nonlop,gvnlx1_,vectout_dum) 540 !$OMP PARALLEL DO 541 do ipw=1,npw1*my_nspinor 542 gvnlx1_(:,ipw)=nonlop_out(:,ipw) 543 end do 544 545 ! PAW: Compute part of H^(1) due to derivative of S 546 cpopt=4*usecprj ; choice=51 ; paw_opt=3 ; signs=2 547 call nonlop(choice,cpopt,cwaveprj_ptr,enlout,gs_hamkq,idir,(/lambda/),mpi_enreg,1,nnlout,& 548 & paw_opt,signs,nonlop_out,tim_nonlop,cwave,vectout_dum) 549 if(compute_conjugate) then 550 !$OMP PARALLEL DO 551 do ipw=1,npw1*my_nspinor ! Note the multiplication by -i 552 gvnlx1_(1,ipw)=gvnlx1_(1,ipw)+nonlop_out(2,ipw) 553 gvnlx1_(2,ipw)=gvnlx1_(2,ipw)-nonlop_out(1,ipw) 554 end do 555 else 556 !$OMP PARALLEL DO 557 do ipw=1,npw1*my_nspinor ! Note the multiplication by i 558 gvnlx1_(1,ipw)=gvnlx1_(1,ipw)-nonlop_out(2,ipw) 559 gvnlx1_(2,ipw)=gvnlx1_(2,ipw)+nonlop_out(1,ipw) 560 end do 561 end if 562 563 ! PAW: Compute part of H^(1) due to derivative of electric field part of Dij 564 cpopt=2 ; choice=1 ; paw_opt=1 ; signs=2 565 call nonlop(choice,cpopt,cwaveprj_ptr,enlout,gs_hamkq,0,(/lambda/),mpi_enreg,1,nnlout,& 566 & paw_opt,signs,svectout_dum,tim_nonlop,cwave,nonlop_out,enl=rf_hamkq%e1kbfr) 567 !$OMP PARALLEL DO 568 do ipw=1,npw1*my_nspinor 569 gvnlx1_(:,ipw)=gvnlx1_(:,ipw)+nonlop_out(:,ipw) 570 end do 571 ABI_FREE(nonlop_out) 572 573 end if ! opt_gvnlx1==2 574 575 ! PAW: Compute derivatives due to part of D_ij^(1) depending on VHxc^(1) 576 if (optnl>=2) then 577 ABI_MALLOC(gvnl2,(2,npw1*my_nspinor)) 578 cpopt=-1+3*usecprj;if (opt_gvnlx1==2) cpopt=2 579 choice=1 ; paw_opt=1 ; signs=2 580 call nonlop(choice,cpopt,cwaveprj_ptr,enlout,gs_hamkq,0,(/lambda/),mpi_enreg,1,nnlout,& 581 & paw_opt,signs,svectout_dum,tim_nonlop,cwave,gvnl2,enl=rf_hamkq%e1kbsc) 582 end if 583 584 if (sij_opt==1) then 585 !$OMP PARALLEL DO 586 do ipw=1,npw1*my_nspinor 587 gs1c(:,ipw)=zero 588 end do 589 end if 590 if (usecprj==0) then 591 call pawcprj_free(cwaveprj_tmp) 592 ABI_FREE(cwaveprj_tmp) 593 end if 594 nullify(cwaveprj_ptr) 595 end if ! PAW 596 597 ! Electric field perturbation with Berry phase 598 ! ------------------------------------------- 599 else if(ipert==natom+2 .and. & 600 & (berryopt==4 .or. berryopt==6 .or. berryopt==7 .or. & 601 & berryopt==14 .or. berryopt==16 .or. berryopt==17 ) .and.(optnl>0.or.sij_opt/=0))then 602 603 if (optnl>=1) then 604 do ipw=1,npw1*my_nspinor 605 gvnlx1_(1,ipw)=-grad_berry(2,ipw) 606 gvnlx1_(2,ipw)= grad_berry(1,ipw) 607 end do 608 end if 609 if (sij_opt==1) then 610 !$OMP PARALLEL DO 611 do ipw=1,npw1*my_nspinor 612 gs1c(:,ipw)=zero 613 end do 614 end if 615 616 ! Strain perturbation 617 ! ------------------------------------------- 618 else if ((ipert==natom+3.or.ipert==natom+4).and.(optnl>0.or.sij_opt/=0)) then 619 620 istr=idir;if(ipert==natom+4) istr=istr+3 621 622 ! PAW: 623 if (gs_hamkq%usepaw==1) then 624 625 if (usecprj==1) then 626 cwaveprj_ptr => cwaveprj 627 else 628 ABI_MALLOC(cwaveprj_tmp,(natom,my_nspinor)) 629 call pawcprj_alloc(cwaveprj_tmp,1,gs_hamkq%dimcprj) 630 cwaveprj_ptr => cwaveprj_tmp 631 end if 632 633 ! 1- Compute derivatives due to projectors |p_i>^(1) 634 ! All atoms contribute 635 cpopt=-1+5*usecprj ; choice=3 ; signs=2 636 paw_opt=1;if (sij_opt/=0) paw_opt=sij_opt+3 637 call nonlop(choice,cpopt,cwaveprj_ptr,enlout,gs_hamkq,istr,(/lambda/),mpi_enreg,1,nnlout,& 638 & paw_opt,signs,gs1c,tim_nonlop,cwave,gvnlx1_) 639 640 ! 2- Compute derivatives due to frozen part of D_ij^(1) (independent of VHxc^(1)) 641 ! All atoms contribute 642 if (optnl>=1) then 643 ABI_MALLOC(nonlop_out,(2,npw1*my_nspinor)) 644 cpopt=1+3*usecprj ; choice=1 ; signs=2 ; paw_opt=1 645 call nonlop(choice,cpopt,cwaveprj_ptr,enlout,gs_hamkq,istr,(/lambda/),mpi_enreg,1,nnlout,& 646 & paw_opt,signs,svectout_dum,tim_nonlop,cwave,nonlop_out,enl=rf_hamkq%e1kbfr) 647 !$OMP PARALLEL DO 648 do ipw=1,npw1*my_nspinor 649 gvnlx1_(:,ipw)=gvnlx1_(:,ipw)+nonlop_out(:,ipw) 650 end do 651 ABI_FREE(nonlop_out) 652 end if 653 654 ! 3- Compute derivatives due to part of D_ij^(1) depending on VHxc^(1) 655 ! All atoms contribute 656 if (optnl>=2) then 657 ABI_MALLOC(gvnl2,(2,npw1*my_nspinor)) 658 cpopt=4 ; choice=1 ; signs=2 ; paw_opt=1 659 call nonlop(choice,cpopt,cwaveprj_ptr,enlout,gs_hamkq,istr,(/lambda/),mpi_enreg,1,nnlout,& 660 & paw_opt,signs,svectout_dum,tim_nonlop,cwave,gvnl2,enl=rf_hamkq%e1kbsc) 661 end if 662 663 if (usecprj==0) then 664 call pawcprj_free(cwaveprj_tmp) 665 ABI_FREE(cwaveprj_tmp) 666 end if 667 nullify(cwaveprj_ptr) 668 669 ! Norm-conserving psps: 670 else 671 ! Compute only derivatives due to projectors |p_i>^(1) 672 choice=3 ; cpopt=-1 ; signs=2 ; paw_opt=0 673 call nonlop(choice,cpopt,cwaveprj,enlout,gs_hamkq,istr,(/lambda/),mpi_enreg,1,nnlout,& 674 & paw_opt,signs,svectout_dum,tim_nonlop,cwave,gvnlx1_) 675 if (sij_opt==1) then 676 !$OMP PARALLEL DO 677 do ipw=1,npw1*my_nspinor 678 gs1c(:,ipw)=zero 679 end do 680 end if 681 end if 682 683 ! No non-local part 684 ! ------------------------------------------- 685 else if (usevnl>0.or.(sij_opt/=0)) then 686 687 if (optnl>=1) then 688 !$OMP PARALLEL DO 689 do ipw=1,npw1*my_nspinor 690 gvnlx1_(:,ipw)=zero 691 end do 692 end if 693 if (sij_opt/=0) then 694 !$OMP PARALLEL DO 695 do ipw=1,npw1*my_nspinor 696 gs1c(:,ipw)=zero 697 end do 698 end if 699 700 end if 701 702 !====================================================================== 703 !== Apply the 1st-order kinetic operator to the wavefunction 704 !== (add it to nl contribution) 705 !====================================================================== 706 707 !Phonon perturbation or Electric field perturbation 708 !------------------------------------------- 709 !No kinetic contribution 710 711 !k-point perturbation or Strain perturbation 712 !------------------------------------------- 713 714 usevnl2=allocated(gvnl2) 715 has_kin=(ipert==natom+1.or.ipert==natom+3.or.ipert==natom+4) 716 if (associated(gs_hamkq%kinpw_kp)) then 717 kinpw1 => gs_hamkq%kinpw_kp 718 else if (optnl>=1.or.usevnl2.or.has_kin) then 719 ABI_BUG('need kinpw1 allocated!') 720 end if 721 if (associated(rf_hamkq%dkinpw_k)) then 722 dkinpw => rf_hamkq%dkinpw_k 723 else if (has_kin) then 724 ABI_BUG('need dkinpw allocated!') 725 end if 726 727 if (has_kin) then ! This is the correct line 728 !DEBUG 729 !if (.false.)then 730 !ENDDEBUG 731 ! Remember that npw=npw1 for ddk perturbation 732 do ispinor=1,my_nspinor 733 !$OMP PARALLEL DO PRIVATE(ipw,ipws) SHARED(cwave,ispinor,gvnlx1_,dkinpw,kinpw1,npw,my_nspinor) 734 do ipw=1,npw 735 ipws=ipw+npw*(ispinor-1) 736 if(kinpw1(ipw)<huge(zero)*1.d-11)then 737 gvnlx1_(1,ipws)=gvnlx1_(1,ipws)+dkinpw(ipw)*cwave(1,ipws) 738 gvnlx1_(2,ipws)=gvnlx1_(2,ipws)+dkinpw(ipw)*cwave(2,ipws) 739 else 740 gvnlx1_(1,ipws)=zero 741 gvnlx1_(2,ipws)=zero 742 end if 743 end do 744 end do 745 end if 746 747 !====================================================================== 748 !== Apply the 1st-order nuclear dipole operator to the wavefunction 749 !== Only coded for DDK 750 !== (add it to nl contribution) 751 !====================================================================== 752 753 has_nd1=( (ipert .EQ. natom+1) .AND. ASSOCIATED(rf_hamkq%vectornd) ) 754 755 if (has_nd1) then 756 ABI_MALLOC(gh1ndc,(2,npw*my_nspinor)) 757 ! this is hard coded for ndat = 1 758 call getgh1ndc(cwave,gh1ndc,gs_hamkq%gbound_k,gs_hamkq%istwf_k,gs_hamkq%kg_k,& 759 & gs_hamkq%mgfft,mpi_enreg,1,gs_hamkq%ngfft,npw,gs_hamkq%nvloc,& 760 & gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,my_nspinor,rf_hamkq%vectornd,& 761 & gs_hamkq%gpu_option) 762 do ispinor=1,my_nspinor 763 do ipw=1,npw 764 ipws=ipw+npw*(ispinor-1) 765 gvnlx1_(1,ipws)=gvnlx1_(1,ipws)+gh1ndc(1,ipws) 766 gvnlx1_(2,ipws)=gvnlx1_(2,ipws)+gh1ndc(2,ipws) 767 end do 768 end do 769 ABI_FREE(gh1ndc) 770 end if 771 772 !====================================================================== 773 !== Apply the 1st-order mGGA operator to the wavefunction 774 !== Only coded for DDK 775 !== (add it to nl contribution) 776 !====================================================================== 777 778 has_mGGA1=( (ipert .EQ. natom+1) .AND. ASSOCIATED(rf_hamkq%vxctaulocal) ) 779 780 if (has_mGGA1) then 781 ABI_MALLOC(gh1c_mGGA,(2,npw*my_nspinor)) 782 ! this is hard coded for ndat = 1 783 call getgh1c_mGGA(cwave,dkinpw,gs_hamkq%gbound_k,gh1c_mGGA,gs_hamkq%gprimd,idir,gs_hamkq%istwf_k,& 784 & gs_hamkq%kg_k,kinpw1,gs_hamkq%mgfft,mpi_enreg,my_nspinor,gs_hamkq%n4,gs_hamkq%n5,& 785 & gs_hamkq%n6,1,gs_hamkq%ngfft,npw,gs_hamkq%nvloc,rf_hamkq%vxctaulocal,& 786 & gs_hamkq%gpu_option) 787 do ispinor=1,my_nspinor 788 do ipw=1,npw 789 ipws=ipw+npw*(ispinor-1) 790 gvnlx1_(1,ipws)=gvnlx1_(1,ipws)+gh1c_mGGA(1,ipws) 791 gvnlx1_(2,ipws)=gvnlx1_(2,ipws)+gh1c_mGGA(2,ipws) 792 end do 793 end do 794 ABI_FREE(gh1c_mGGA) 795 end if 796 797 !====================================================================== 798 !== Sum contributions to get the application of H^(1) to the wf 799 !====================================================================== 800 !Also filter the wavefunctions for large modified kinetic energy 801 802 !Add non-local+kinetic to local part 803 if (optnl>=1.or.has_kin) then 804 do ispinor=1,my_nspinor 805 ipws=(ispinor-1)*npw1 806 !$OMP PARALLEL DO PRIVATE(ipw) SHARED(gh1c,gvnlx1_,kinpw1,ipws,npw1) 807 do ipw=1+ipws,npw1+ipws 808 if(kinpw1(ipw-ipws)<huge(zero)*1.d-11)then 809 gh1c(1,ipw)=gh1c(1,ipw)+gvnlx1_(1,ipw) 810 gh1c(2,ipw)=gh1c(2,ipw)+gvnlx1_(2,ipw) 811 else 812 gh1c(1,ipw)=zero 813 gh1c(2,ipw)=zero 814 end if 815 end do 816 end do 817 end if 818 819 !PAW: add non-local part due to first order change of VHxc 820 if (usevnl2) then 821 do ispinor=1,my_nspinor 822 ipws=(ispinor-1)*npw1 823 !$OMP PARALLEL DO PRIVATE(ipw) SHARED(gh1c,gvnl2,kinpw1,ipws,npw1) 824 do ipw=1+ipws,npw1+ipws 825 if(kinpw1(ipw-ipws)<huge(zero)*1.d-11)then 826 gh1c(1,ipw)=gh1c(1,ipw)+gvnl2(1,ipw) 827 gh1c(2,ipw)=gh1c(2,ipw)+gvnl2(2,ipw) 828 end if 829 end do 830 end do 831 ABI_FREE(gvnl2) 832 end if 833 834 if (usevnl==1) then 835 nullify(gvnlx1_) 836 else 837 ABI_FREE(gvnlx1_) 838 end if 839 840 call timab(196+tim_getgh1c,2,tsec) 841 842 DBG_EXIT("COLL") 843 844 end subroutine getgh1c
ABINIT/getgh1c_mGGA [ Functions ]
NAME
getgh1c_mGGA
FUNCTION
Compute first order metaGGA contribution to <G|H|C> for input vector |C> expressed in reciprocal space. ONLY FOR DDK PERTURBATION
INPUTS
OUTPUT
gh1c_mGGA(2,npw_k*my_nspinor*ndat)=metaGGA contribution to <G|H1|C>
SIDE EFFECTS
SOURCE
2008 subroutine getgh1c_mGGA(cwavein,dkinpw,gbound_k,gh1c_mGGA,gprimd,idir,istwf_k,kg_k,& 2009 & kinpw1,mgfft,mpi_enreg,my_nspinor,n4,n5,n6,ndat,ngfft,npw_k,nvloc,vxctaulocal,gpu_option) 2010 2011 !Arguments ------------------------------------ 2012 !scalars 2013 integer,intent(in) :: idir,istwf_k,mgfft,my_nspinor,n4,n5,n6,ndat,npw_k,nvloc 2014 integer,intent(in),optional :: gpu_option 2015 type(MPI_type),intent(in) :: mpi_enreg 2016 !arrays 2017 integer,intent(in) :: gbound_k(2*mgfft+4),kg_k(3,npw_k),ngfft(18) 2018 real(dp),intent(in) :: gprimd(3,3) 2019 real(dp),intent(inout) :: cwavein(2,npw_k*my_nspinor*ndat) 2020 real(dp),intent(inout) :: gh1c_mGGA(2,npw_k*my_nspinor*ndat) 2021 real(dp),intent(inout) :: vxctaulocal(n4,n5,n6,nvloc,4) 2022 real(dp),pointer,intent(in) :: dkinpw(:),kinpw1(:) 2023 2024 !Local variables------------------------------- 2025 !scalars 2026 integer :: idat,ii,ipw,nspinortot,shift,gpu_option_ 2027 integer,parameter :: tim_fourwf=1 2028 real(dp) :: weight=one 2029 logical :: nspinor1TreatedByThisProc,nspinor2TreatedByThisProc 2030 !arrays 2031 real(dp),allocatable :: cwavein1(:,:),cwavein2(:,:),dgcwavef(:,:) 2032 real(dp),allocatable :: ghc1(:,:),ghc2(:,:),work(:,:,:,:) 2033 2034 if(present(gpu_option)) then 2035 gpu_option_=gpu_option 2036 else 2037 gpu_option_=0 2038 end if 2039 2040 gh1c_mGGA(:,:)=zero 2041 if (nvloc/=1) return 2042 2043 nspinortot=min(2,(1+mpi_enreg%paral_spinor)*my_nspinor) 2044 if (mpi_enreg%paral_spinor==0) then 2045 shift=npw_k 2046 nspinor1TreatedByThisProc=.true. 2047 nspinor2TreatedByThisProc=(nspinortot==2) 2048 else 2049 shift=0 2050 nspinor1TreatedByThisProc=(mpi_enreg%me_spinor==0) 2051 nspinor2TreatedByThisProc=(mpi_enreg%me_spinor==1) 2052 end if 2053 2054 ABI_MALLOC(work,(2,n4,n5,n6*ndat)) 2055 2056 if (nspinortot==1) then 2057 2058 ABI_MALLOC(ghc1,(2,npw_k*ndat)) 2059 ABI_MALLOC(ghc2,(2,npw_k*ndat)) 2060 ABI_MALLOC(dgcwavef,(2,npw_k*ndat)) 2061 2062 ! From -1/2 (grad vxctau) \cdot (grad \psi) = 2063 ! -1/2 (grad vxctau)\cdot(2\pi i (k + G)\psi), the 2064 do ii=1,3 2065 dgcwavef = zero 2066 do idat=1,ndat 2067 do ipw=1,npw_k 2068 dgcwavef(1,ipw+(idat-1)*npw_k)= cwavein(2,ipw+(idat-1)*npw_k)*two_pi*gprimd(ii,idir) 2069 dgcwavef(2,ipw+(idat-1)*npw_k)=-cwavein(1,ipw+(idat-1)*npw_k)*two_pi*gprimd(ii,idir) 2070 end do 2071 end do 2072 call fourwf(1,vxctaulocal(:,:,:,:,1+ii),dgcwavef,ghc1,work,gbound_k,gbound_k,& 2073 & istwf_k,kg_k,kg_k,mgfft,mpi_enreg,ndat,ngfft,npw_k,npw_k,n4,n5,n6,2,& 2074 & tim_fourwf,weight,weight,gpu_option=gpu_option_) 2075 gh1c_mGGA(:,:) = gh1c_mGGA(:,:) - half*ghc1 2076 end do ! ii 2077 2078 ! From -1/2 vxctau (grad . grad \psi), the k derivative is 2079 ! vxctau \times dkinpw_dir * \psi 2080 do ipw=1,npw_k 2081 if(kinpw1(ipw)<huge(zero)*1.d-11)then 2082 ghc1(:,ipw)=dkinpw(ipw)*cwavein(:,ipw) 2083 else 2084 ghc1(:,ipw) = zero 2085 end if 2086 end do 2087 call fourwf(1,vxctaulocal(:,:,:,:,1),ghc1,ghc2,work,gbound_k,gbound_k,& 2088 & istwf_k,kg_k,kg_k,mgfft,mpi_enreg,ndat,ngfft,npw_k,npw_k,n4,n5,n6,2,& 2089 & tim_fourwf,weight,weight,gpu_option=gpu_option_) 2090 2091 gh1c_mGGA = gh1c_mGGA + ghc2 2092 2093 ! JWZ debug 2094 ! gh1c_mGGA = zero 2095 2096 ABI_FREE(ghc1) 2097 ABI_FREE(ghc2) 2098 ABI_FREE(dgcwavef) 2099 2100 else ! nspinortot==2 2101 2102 ABI_MALLOC(cwavein1,(2,npw_k*ndat)) 2103 ABI_MALLOC(cwavein2,(2,npw_k*ndat)) 2104 do idat=1,ndat 2105 do ipw=1,npw_k 2106 cwavein1(1:2,ipw+(idat-1)*npw_k)=cwavein(1:2,ipw+(idat-1)*my_nspinor*npw_k) 2107 cwavein2(1:2,ipw+(idat-1)*npw_k)=cwavein(1:2,ipw+(idat-1)*my_nspinor*npw_k+shift) 2108 end do 2109 end do 2110 2111 ABI_MALLOC(ghc1,(2,npw_k*ndat)) 2112 ABI_MALLOC(ghc2,(2,npw_k*ndat)) 2113 2114 if (nspinor1TreatedByThisProc) then 2115 2116 ! call fourwf(1,vxctaulocal(:,:,:,:,1+idir),cwavein1,ghc1,work,gbound_k,gbound_k,& 2117 ! & istwf_k,kg_k,kg_k,mgfft,mpi_enreg,ndat,ngfft,npw_k,npw_k,n4,n5,n6,2,& 2118 ! & tim_fourwf,weight,weight,use_gpu_cuda=use_gpu_cuda_) 2119 2120 ! ! scale by -1/2 * 2\pi i = -i \pi 2121 ! do idat = 1, ndat 2122 ! do ipw = 1, npw_k 2123 ! gh1c_mGGA(1,ipw+(idat-1)*npw_k) = pi*ghc1(2,ipw+(idat-1)*npw_k) 2124 ! gh1c_mGGA(2,ipw+(idat-1)*npw_k) = -pi*ghc1(1,ipw+(idat-1)*npw_k) 2125 ! end do 2126 ! end do 2127 2128 ! From -1/2 vxctau (grad . grad \psi), the k derivative is 2129 ! vxctau \times dkinpw_dir * \psi 2130 do idat = 1, ndat 2131 do ipw=1,npw_k 2132 if(kinpw1(ipw)<huge(zero)*1.d-11)then 2133 ghc1(1,ipw)=dkinpw(ipw)*cwavein1(1,ipw+(idat-1)*npw_k) 2134 ghc1(2,ipw)=dkinpw(ipw)*cwavein1(2,ipw+(idat-1)*npw_k) 2135 else 2136 ghc1(:,ipw+(idat-1)*npw_k) = zero 2137 end if 2138 end do 2139 end do 2140 2141 call fourwf(1,vxctaulocal(:,:,:,:,1),ghc1,ghc2,work,gbound_k,gbound_k,& 2142 & istwf_k,kg_k,kg_k,mgfft,mpi_enreg,ndat,ngfft,npw_k,npw_k,n4,n5,n6,2,& 2143 & tim_fourwf,weight,weight,gpu_option=gpu_option_) 2144 2145 do idat = 1, ndat 2146 do ipw = 1, npw_k 2147 gh1c_mGGA(:,ipw+(idat-1)*npw_k) = gh1c_mGGA(:,ipw+(idat-1)*npw_k) + & 2148 & ghc2(:,ipw+(idat-1)*npw_k) 2149 end do 2150 end do 2151 2152 end if ! end spinor 1 2153 2154 if (nspinor2TreatedByThisProc) then 2155 2156 ABI_MALLOC(ghc2,(2,npw_k*ndat)) 2157 2158 ! call fourwf(1,vxctaulocal(:,:,:,:,1+idir),cwavein2,ghc2,work,gbound_k,gbound_k,& 2159 ! & istwf_k,kg_k,kg_k,mgfft,mpi_enreg,ndat,ngfft,npw_k,npw_k,n4,n5,n6,2,& 2160 ! & tim_fourwf,weight,weight,use_gpu_cuda=use_gpu_cuda_) 2161 2162 ! ! scale by -1/2 * 2\pi i = -i \pi 2163 ! do idat = 1, ndat 2164 ! do ipw = 1, npw_k 2165 ! gh1c_mGGA(1,ipw+(idat-1)*npw_k) = pi*ghc2(2,ipw+(idat-1)*npw_k) 2166 ! gh1c_mGGA(2,ipw+(idat-1)*npw_k) = -pi*ghc2(1,ipw+(idat-1)*npw_k) 2167 ! end do 2168 ! end do 2169 2170 ! From -1/2 vxctau (grad . grad \psi), the k derivative is 2171 ! vxctau \times dkinpw_dir * \psi 2172 do idat = 1, ndat 2173 do ipw=1,npw_k 2174 if(kinpw1(ipw)<huge(zero)*1.d-11)then 2175 ghc1(1,ipw)=dkinpw(ipw)*cwavein2(1,ipw+(idat-1)*npw_k) 2176 ghc1(2,ipw)=dkinpw(ipw)*cwavein2(2,ipw+(idat-1)*npw_k) 2177 else 2178 ghc1(:,ipw+(idat-1)*npw_k) = zero 2179 end if 2180 end do 2181 end do 2182 2183 call fourwf(1,vxctaulocal(:,:,:,:,1),ghc1,ghc2,work,gbound_k,gbound_k,& 2184 & istwf_k,kg_k,kg_k,mgfft,mpi_enreg,ndat,ngfft,npw_k,npw_k,n4,n5,n6,2,& 2185 & tim_fourwf,weight,weight,gpu_option=gpu_option_) 2186 2187 do idat = 1, ndat 2188 do ipw = 1, npw_k 2189 gh1c_mGGA(:,ipw+(idat-1)*npw_k) = gh1c_mGGA(:,ipw+(idat-1)*npw_k) + & 2190 & ghc2(:,ipw+(idat-1)*npw_k) 2191 end do 2192 end do 2193 2194 end if ! end spinor 2 2195 2196 ABI_FREE(ghc1) 2197 ABI_FREE(ghc2) 2198 2199 ABI_FREE(cwavein1) 2200 ABI_FREE(cwavein2) 2201 2202 end if ! nspinortot 2203 2204 ABI_FREE(work) 2205 2206 end subroutine getgh1c_mGGA
ABINIT/getgh1dqc [ Functions ]
NAME
getgh1dqc
FUNCTION
Computes <G|dH^(1)/dq_{gamma}|C> or <G|d^2H^(1)/dq_{gamma}dq_{delta}|C>
for input vector |C> expressed in reciprocal space.
dH^(1)/dq_{gamma} and d^2H^(1)/dq_{gamma}dq_{delta} are the first
and second q-gradient (at q=0) of the 1st-order perturbed Hamiltonian.
The first (second) derivative direction is inferred from idir (qdir1).
INPUTS
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) gs_hamkq <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k+q idir=first index of the perturbation ipert=type of the perturbation 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 1: local part of H^(1) is computed in gvloc1dqc optnl=0: non-local part of H^(1) is not computed 1: non-local part of H^(1) depending on VHxc^(1) is not computed in gvloc1dqc 2: non-local part of H^(1) is totally computed in gvloc1dqc qdir1= direction of the 1st q-gradient rf_hamkq <type(rf_hamiltonian_type)>=all data for the 1st-order Hamiltonian at k,k+q qdir2= (optional) direction of the 2nd q-gradient
OUTPUT
gh1dqc(2,npw1*nspinor)= <G|dH^(1)/dq_{\gamma}|C> on the k+q sphere
gvloc1dqc(2,npw1*nspinor)= local potential part of gh1dqc
gvnl1dqc(2,npw1*nspinor)= non local potential part of gh1dqc
SIDE EFFECTS
NOTES
Currently two Hamiltonian gradients at (q=0) are implemented:
ipert<=natom -> first q-derivative along reduced coordinates directions
of the atomic displacement perturbation hamiltonian
ipert==natom+3 or natom+4 -> second q-derivative along cartesian coordinates
of the metric perturbation hamiltonian.
Which is equivalent (except for an i factor) to the first
q-derivative along cartesian coordinates of the strain
perturbation hamiltonian.
SOURCE
1438 subroutine getgh1dqc(cwave,cwaveprj,gh1dqc,gvloc1dqc,gvnl1dqc,gs_hamkq,& 1439 & idir,ipert,mpi_enreg,optlocal,optnl,qdir1,rf_hamkq,& 1440 & qdir2) !optional 1441 1442 !Arguments ------------------------------------ 1443 !scalars 1444 integer,intent(in) :: idir,ipert,optlocal,optnl,qdir1 1445 integer,intent(in),optional :: qdir2 1446 type(MPI_type),intent(in) :: mpi_enreg 1447 type(gs_hamiltonian_type),intent(inout),target :: gs_hamkq 1448 type(rf_hamiltonian_type),intent(inout),target :: rf_hamkq 1449 1450 !arrays 1451 real(dp),intent(inout) :: cwave(2,gs_hamkq%npw_k*gs_hamkq%nspinor) 1452 real(dp),intent(out) :: gh1dqc(2,gs_hamkq%npw_kp*gs_hamkq%nspinor) 1453 real(dp),intent(out) :: gvloc1dqc(2,gs_hamkq%npw_kp*gs_hamkq%nspinor) 1454 real(dp),intent(out) :: gvnl1dqc(2,gs_hamkq%npw_kp*gs_hamkq%nspinor) 1455 real(dp),pointer :: dqdqkinpw(:),kinpw1(:) 1456 type(pawcprj_type),intent(inout),target :: cwaveprj(:,:) 1457 1458 !Local variables------------------------------- 1459 !scalars 1460 integer :: choice,cpopt,iidir,ipw,ipws,ispinor,my_nspinor,natom,nnlout 1461 integer :: npw,npw1,paw_opt,signs,tim_fourwf,tim_nonlop 1462 logical :: has_kin 1463 character(len=500) :: msg 1464 real(dp) :: lambda,weight 1465 1466 !arrays 1467 integer,parameter :: ngamma(3,3)=reshape((/1,6,5,9,2,4,8,7,3/),(/3,3/)) 1468 real(dp) :: enlout(1),svectout_dum(1,1) 1469 real(dp),ABI_CONTIGUOUS pointer :: gvnl1dqc_(:,:) 1470 real(dp), allocatable :: work(:,:,:,:) 1471 1472 ! ************************************************************************* 1473 1474 DBG_ENTER("COLL") 1475 1476 !====================================================================== 1477 !== Initialisations and compatibility tests 1478 !====================================================================== 1479 1480 npw =gs_hamkq%npw_k 1481 npw1 =gs_hamkq%npw_kp 1482 natom=gs_hamkq%natom 1483 1484 !Compatibility tests 1485 if (mpi_enreg%paral_spinor==1) then 1486 msg='Not compatible with parallelization over spinorial components !' 1487 ABI_BUG(msg) 1488 end if 1489 1490 !Check sizes 1491 my_nspinor=max(1,gs_hamkq%nspinor/mpi_enreg%nproc_spinor) 1492 if (size(cwave)<2*npw*my_nspinor) then 1493 msg='wrong size for cwave!' 1494 ABI_BUG(msg) 1495 end if 1496 if (size(gh1dqc)<2*npw1*my_nspinor) then 1497 msg='wrong size for gh1dqc!' 1498 ABI_BUG(msg) 1499 end if 1500 1501 !============================================================================= 1502 !== Apply the q-gradients of the 1st-order local potential to the wavefunction 1503 !============================================================================= 1504 1505 !Phonon and metric (strain) perturbation 1506 if (ipert<=natom+5.and.ipert/=natom+1.and.ipert/=natom+2.and.optlocal>0) then 1507 1508 ABI_MALLOC(work,(2,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6)) 1509 1510 weight=one ; tim_fourwf=4 1511 call fourwf(rf_hamkq%cplex,rf_hamkq%vlocal1,cwave,gvloc1dqc,work,gs_hamkq%gbound_k,gs_hamkq%gbound_kp,& 1512 & gs_hamkq%istwf_k,gs_hamkq%kg_k,gs_hamkq%kg_kp,gs_hamkq%mgfft,mpi_enreg,1,gs_hamkq%ngfft,& 1513 & npw,npw1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,2,tim_fourwf,weight,weight,& 1514 & gpu_option=gs_hamkq%gpu_option) 1515 1516 ABI_FREE(work) 1517 1518 else 1519 1520 !$OMP PARALLEL DO 1521 do ipw=1,npw1*my_nspinor 1522 gvloc1dqc(:,ipw)=zero 1523 end do 1524 1525 end if 1526 1527 !================================================================================ 1528 !== Apply the q-gradients of the 1st-order non-local potential to the wavefunction 1529 !================================================================================ 1530 1531 !Initializations 1532 lambda=zero 1533 nnlout=1 1534 tim_nonlop=0 1535 1536 !Allocations 1537 ABI_MALLOC(gvnl1dqc_,(2,npw1*my_nspinor)) 1538 1539 !Phonon perturbation 1540 !------------------------------------------- 1541 !1st q-gradient 1542 if (ipert<=natom.and..not.present(qdir2).and.optnl>0) then 1543 cpopt=-1 ; choice=22 ; signs=2 ; paw_opt=0 1544 call nonlop(choice,cpopt,cwaveprj,enlout,gs_hamkq,idir,(/lambda/),mpi_enreg,1,nnlout,& 1545 & paw_opt,signs,svectout_dum,tim_nonlop,cwave,gvnl1dqc_,iatom_only=ipert,qdir=qdir1) 1546 1547 !$OMP PARALLEL DO 1548 do ipw=1,npw1*my_nspinor 1549 gvnl1dqc(1,ipw)=gvnl1dqc_(1,ipw) 1550 gvnl1dqc(2,ipw)=gvnl1dqc_(2,ipw) 1551 end do 1552 1553 !2nd q-gradient 1554 else if (ipert<=natom.and.present(qdir2).and.optnl>0) then 1555 iidir=ngamma(idir,qdir2) 1556 cpopt=-1 ; choice=25 ; signs=2 ; paw_opt=0 1557 call nonlop(choice,cpopt,cwaveprj,enlout,gs_hamkq,iidir,(/lambda/),mpi_enreg,1,nnlout,& 1558 & paw_opt,signs,svectout_dum,tim_nonlop,cwave,gvnl1dqc_,iatom_only=ipert,qdir=qdir1) 1559 1560 !$OMP PARALLEL DO 1561 do ipw=1,npw1*my_nspinor 1562 gvnl1dqc(1,ipw)=gvnl1dqc_(1,ipw) 1563 gvnl1dqc(2,ipw)=gvnl1dqc_(2,ipw) 1564 end do 1565 1566 !Metric (strain) perturbation 1567 !------------------------------------------- 1568 else if ((ipert==natom+3.or.ipert==natom+4).and.optnl>0) then 1569 cpopt=-1 ; choice=33 ; signs=2 ; paw_opt=0 1570 call nonlop(choice,cpopt,cwaveprj,enlout,gs_hamkq,idir,(/lambda/),mpi_enreg,1,nnlout,& 1571 & paw_opt,signs,svectout_dum,tim_nonlop,cwave,gvnl1dqc_,qdir=qdir1) 1572 1573 !$OMP PARALLEL DO 1574 do ipw=1,npw1*my_nspinor 1575 gvnl1dqc(1,ipw)=gvnl1dqc_(1,ipw) 1576 gvnl1dqc(2,ipw)=gvnl1dqc_(2,ipw) 1577 end do 1578 1579 else 1580 1581 !$OMP PARALLEL DO 1582 do ipw=1,npw1*my_nspinor 1583 gvnl1dqc(:,ipw)=zero 1584 end do 1585 1586 end if 1587 1588 !============================================================================== 1589 !== Apply the q-gradients of the 1st-order kinetic operator to the wavefunction 1590 !== (add it to nl contribution) 1591 !============================================================================== 1592 1593 !Strain (metric) perturbation 1594 !------------------------------------------- 1595 has_kin=(ipert==natom+3.or.ipert==natom+4) 1596 if (associated(gs_hamkq%kinpw_kp)) then 1597 kinpw1 => gs_hamkq%kinpw_kp 1598 else if (has_kin) then 1599 msg='need kinpw1 allocated!' 1600 ABI_BUG(msg) 1601 end if 1602 if (associated(rf_hamkq%dkinpw_k)) then 1603 dqdqkinpw => rf_hamkq%dkinpw_k 1604 else if (has_kin) then 1605 msg='need dqdqkinpw allocated!' 1606 ABI_BUG(msg) 1607 end if 1608 1609 if (has_kin) then 1610 ! Remember that npw=npw1 1611 do ispinor=1,my_nspinor 1612 !$OMP PARALLEL DO PRIVATE(ipw,ipws) SHARED(cwave,ispinor,gvnl1dqc,dqdqkinpw,kinpw1,npw,my_nspinor) 1613 do ipw=1,npw 1614 ipws=ipw+npw*(ispinor-1) 1615 if(kinpw1(ipw)<huge(zero)*1.d-11)then 1616 gvnl1dqc(1,ipws)=gvnl1dqc(1,ipws)+dqdqkinpw(ipw)*cwave(1,ipws) 1617 gvnl1dqc(2,ipws)=gvnl1dqc(2,ipws)+dqdqkinpw(ipw)*cwave(2,ipws) 1618 else 1619 gvnl1dqc(1,ipws)=zero 1620 gvnl1dqc(2,ipws)=zero 1621 end if 1622 end do 1623 end do 1624 end if 1625 1626 !=================================================================================== 1627 !== Sum contributions to get the application of dH^(1)/dq or d^2H^(1)/dqdq to the wf 1628 !=================================================================================== 1629 1630 do ispinor=1,my_nspinor 1631 ipws=(ispinor-1)*npw1 1632 !$OMP PARALLEL DO PRIVATE(ipw) SHARED(gh1dqc,gvnl1dqc,kinpw1,ipws,npw1) 1633 do ipw=1+ipws,npw1+ipws 1634 if(kinpw1(ipw-ipws)<huge(zero)*1.d-11)then 1635 gh1dqc(1,ipw)=gvloc1dqc(1,ipw)+gvnl1dqc(1,ipw) 1636 gh1dqc(2,ipw)=gvloc1dqc(2,ipw)+gvnl1dqc(2,ipw) 1637 else 1638 gh1dqc(1,ipw)=zero 1639 gh1dqc(2,ipw)=zero 1640 end if 1641 end do 1642 end do 1643 1644 ABI_FREE(gvnl1dqc_) 1645 DBG_EXIT("COLL") 1646 1647 end subroutine getgh1dqc
ABINIT/getgh1ndc [ Functions ]
NAME
getgh1ndc
FUNCTION
Compute 1st order magnetic nuclear dipole moment contribution to <G|H|C> for input vector |C> expressed in reciprocal space. Only for DDK perturbation
INPUTS
OUTPUT
gh1ndc(2,npw_k*my_nspinor*ndat)=1st order A.p contribution to <G|H|C> for array of nuclear dipoles
SIDE EFFECTS
NOTES
This codes only the DDK response for A.p, so effectively A_ipert|C>. The nuclear dipole Hamiltonian (to first order in the nuclear dipole strength) is A.p where in atomic units A.p=\alpha^2 m x (r-R)/(r-R)^3 . p. Here the components of A have been precomputed in real space by make_vectornd. The first-order DDK contribution is d A.p/dk = A_idir where idir is the direction of the DDK perturbation, or 2\pi A_idir when k is given in reduced coords as is usual
SOURCE
1877 subroutine getgh1ndc(cwavein,gh1ndc,gbound_k,istwf_k,kg_k,mgfft,mpi_enreg,& 1878 & ndat,ngfft,npw_k,nvloc,n4,n5,n6,my_nspinor,vectornd,gpu_option) 1879 1880 !Arguments ------------------------------------ 1881 !scalars 1882 integer,intent(in) :: istwf_k,mgfft,my_nspinor,ndat,npw_k,nvloc,n4,n5,n6,gpu_option 1883 type(MPI_type),intent(in) :: mpi_enreg 1884 !arrays 1885 integer,intent(in) :: gbound_k(2*mgfft+4),kg_k(3,npw_k),ngfft(18) 1886 real(dp),intent(inout) :: cwavein(2,npw_k*my_nspinor*ndat) 1887 real(dp),intent(inout) :: gh1ndc(2,npw_k*my_nspinor*ndat) 1888 real(dp),intent(inout) :: vectornd(n4,n5,n6,nvloc) 1889 1890 !Local variables------------------------------- 1891 !scalars 1892 integer,parameter :: tim_fourwf=1 1893 integer :: idat,ipw,nspinortot,shift 1894 logical :: nspinor1TreatedByThisProc,nspinor2TreatedByThisProc 1895 real(dp) :: weight=one 1896 !arrays 1897 real(dp),allocatable :: cwavein1(:,:),cwavein2(:,:) 1898 real(dp),allocatable :: ghc1(:,:),ghc2(:,:) 1899 real(dp),allocatable :: work(:,:,:,:) 1900 1901 ! ********************************************************************* 1902 1903 gh1ndc(:,:)=zero 1904 if (nvloc/=1) return 1905 1906 nspinortot=min(2,(1+mpi_enreg%paral_spinor)*my_nspinor) 1907 if (mpi_enreg%paral_spinor==0) then 1908 shift=npw_k 1909 nspinor1TreatedByThisProc=.true. 1910 nspinor2TreatedByThisProc=(nspinortot==2) 1911 else 1912 shift=0 1913 nspinor1TreatedByThisProc=(mpi_enreg%me_spinor==0) 1914 nspinor2TreatedByThisProc=(mpi_enreg%me_spinor==1) 1915 end if 1916 1917 ABI_MALLOC(work,(2,n4,n5,n6*ndat)) 1918 1919 if (nspinortot==1) then 1920 1921 ABI_MALLOC(ghc1,(2,npw_k*ndat)) 1922 1923 ! apply vector potential in direction ipert to input wavefunction 1924 call fourwf(1,vectornd,cwavein,ghc1,work,gbound_k,gbound_k,& 1925 & istwf_k,kg_k,kg_k,mgfft,mpi_enreg,ndat,ngfft,npw_k,npw_k,n4,n5,n6,2,& 1926 & tim_fourwf,weight,weight,gpu_option=gpu_option) 1927 1928 ! scale by 2\pi\alpha^2 1929 gh1ndc=two_pi*FineStructureConstant2*ghc1 1930 1931 ABI_FREE(ghc1) 1932 1933 else ! nspinortot==2 1934 1935 ABI_MALLOC(cwavein1,(2,npw_k*ndat)) 1936 ABI_MALLOC(cwavein2,(2,npw_k*ndat)) 1937 do idat=1,ndat 1938 do ipw=1,npw_k 1939 cwavein1(1:2,ipw+(idat-1)*npw_k)=cwavein(1:2,ipw+(idat-1)*my_nspinor*npw_k) 1940 cwavein2(1:2,ipw+(idat-1)*npw_k)=cwavein(1:2,ipw+(idat-1)*my_nspinor*npw_k+shift) 1941 end do 1942 end do 1943 1944 if (nspinor1TreatedByThisProc) then 1945 1946 ABI_MALLOC(ghc1,(2,npw_k*ndat)) 1947 1948 call fourwf(1,vectornd,cwavein1,ghc1,work,gbound_k,gbound_k,& 1949 & istwf_k,kg_k,kg_k,mgfft,mpi_enreg,ndat,ngfft,npw_k,npw_k,n4,n5,n6,2,& 1950 & tim_fourwf,weight,weight,gpu_option=gpu_option) 1951 1952 do idat=1,ndat 1953 do ipw=1,npw_k 1954 gh1ndc(1:2,ipw+(idat-1)*npw_k)=two_pi*FineStructureConstant2*ghc1(1:2,ipw+(idat-1)*npw_k) 1955 end do 1956 end do 1957 1958 ABI_FREE(ghc1) 1959 1960 end if ! end spinor 1 1961 1962 if (nspinor2TreatedByThisProc) then 1963 1964 ABI_MALLOC(ghc2,(2,npw_k*ndat)) 1965 1966 call fourwf(1,vectornd,cwavein2,ghc2,work,gbound_k,gbound_k,& 1967 & istwf_k,kg_k,kg_k,mgfft,mpi_enreg,ndat,ngfft,npw_k,npw_k,n4,n5,n6,2,& 1968 & tim_fourwf,weight,weight,gpu_option=gpu_option) 1969 1970 do idat=1,ndat 1971 do ipw=1,npw_k 1972 gh1ndc(1:2,ipw+(idat-1)*npw_k+shift)=two_pi*FineStructureConstant2*ghc2(1:2,ipw+(idat-1)*npw_k) 1973 end do 1974 end do 1975 1976 ABI_FREE(ghc2) 1977 1978 end if ! end spinor 2 1979 1980 ABI_FREE(cwavein1) 1981 ABI_FREE(cwavein2) 1982 1983 end if ! nspinortot 1984 1985 ABI_FREE(work) 1986 1987 end subroutine getgh1ndc
ABINIT/m_getgh1c [ Modules ]
NAME
m_getgh1c
FUNCTION
COPYRIGHT
Copyright (C) 1998-2024 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 .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_getgh1c 23 24 use defs_basis 25 use m_abicore 26 use m_errors 27 use m_dtset 28 use m_xmpi 29 30 use defs_abitypes, only : MPI_type 31 use defs_datatypes, only : pseudopotential_type 32 use m_time, only : timab 33 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free, & 34 & pawcprj_copy, pawcprj_lincom, pawcprj_axpby, pawcprj_mpi_sum 35 use m_kg, only : kpgstr, mkkin, mkkpg, mkkin_metdqdq 36 use m_mkffnl, only : mkffnl 37 use m_pawfgr, only : pawfgr_type 38 use m_fft, only : fftpac, fourwf 39 use m_hamiltonian, only : gs_hamiltonian_type, rf_hamiltonian_type 40 use m_cgtools, only : projbd 41 use m_nonlop, only : nonlop 42 use m_fourier_interpol, only : transgrid 43 44 implicit none 45 46 private
m_hamiltonian/getgh1c_setup [ Functions ]
[ Top ] [ m_hamiltonian ] [ Functions ]
NAME
getgh1c_setup
FUNCTION
INPUTS
OUTPUT
SOURCE
968 subroutine getgh1c_setup(gs_hamkq, rf_hamkq, dtset, psps, kpoint, kpq, idir, ipert,& ! In 969 natom, rmet, gprimd, gmet, istwf_k, npw_k, npw1_k, & ! In 970 useylmgr1, kg_k, ylm_k, kg1_k, ylm1_k, ylmgr1_k, & ! In 971 dkinpw, nkpg, nkpg1, kpg_k, kpg1_k, kinpw1, ffnlk, ffnl1, ph3d, ph3d1,& ! Out 972 ddkinpw, dkinpw2, rf_hamk_dir2, ffnl1_test, & ! Optional 973 reuse_kpg_k, reuse_kpg1_k, reuse_ffnlk, reuse_ffnl1) ! Optional 974 975 !Arguments ------------------------------------ 976 !scalars 977 integer,intent(in) :: idir,ipert,istwf_k,npw_k,npw1_k,natom,useylmgr1 978 integer,intent(out) :: nkpg,nkpg1 979 type(gs_hamiltonian_type),intent(inout) :: gs_hamkq 980 type(rf_hamiltonian_type),intent(inout) :: rf_hamkq 981 type(rf_hamiltonian_type),intent(inout),optional :: rf_hamk_dir2 982 type(dataset_type),intent(in) :: dtset 983 type(pseudopotential_type),intent(in) :: psps 984 integer,optional,intent(in) :: reuse_kpg_k, reuse_kpg1_k, reuse_ffnlk, reuse_ffnl1 985 !arrays 986 integer,intent(in) :: kg_k(3,npw_k),kg1_k(3,npw1_k) 987 real(dp),intent(in) :: kpoint(3),kpq(3),gmet(3,3),gprimd(3,3),rmet(3,3) 988 real(dp),intent(in) :: ylm_k(npw_k,psps%mpsang*psps%mpsang*psps%useylm) 989 real(dp),intent(in) :: ylmgr1_k(npw1_k,3+6*((ipert-natom)/10),psps%mpsang*psps%mpsang*psps%useylm*useylmgr1) 990 real(dp),intent(in) :: ylm1_k(npw1_k,psps%mpsang*psps%mpsang*psps%useylm) 991 real(dp),allocatable,intent(out) :: dkinpw(:),kinpw1(:) 992 real(dp),allocatable,intent(inout) :: ffnlk(:,:,:,:),ffnl1(:,:,:,:), kpg_k(:,:), kpg1_k(:,:) 993 real(dp),allocatable,intent(out),optional :: dkinpw2(:),ddkinpw(:),ffnl1_test(:,:,:,:) 994 real(dp),allocatable,intent(out) :: ph3d(:,:,:),ph3d1(:,:,:) 995 996 !Local variables------------------------------- 997 !scalars 998 integer :: dimffnl1,dimffnlk,ider,idir0,idir1,idir2,istr,ntypat,print_info 999 integer :: reuse_ffnlk_, reuse_ffnl1_, reuse_kpg_k_, reuse_kpg1_k_ 1000 logical :: qne0 1001 !real(dp) :: cpu, wall, gflops 1002 !arrays 1003 real(dp) :: ylmgr_dum(1,1,1), tsec(2) 1004 1005 ! ************************************************************************* 1006 1007 ! MG: This routine is called **many times** in the EPH code for phonon and DDK perturbations 1008 ! Please, be extremely careful when adding extra stuff that may affect performance. 1009 1010 ! Keep track of total time spent in getgh1c_setup (use 195 slot) 1011 call timab(195, 1, tsec) 1012 !call cwtime(cpu, wall, gflops, "start") 1013 1014 reuse_ffnlk_ = 0; if (present(reuse_ffnlk)) reuse_ffnlk_ = reuse_ffnlk 1015 reuse_ffnl1_ = 0; if (present(reuse_ffnl1)) reuse_ffnl1_ = reuse_ffnl1 1016 reuse_kpg_k_ = 0; if (present(reuse_kpg_k)) reuse_kpg_k_ = reuse_kpg_k 1017 reuse_kpg1_k_ = 0; if (present(reuse_kpg1_k)) reuse_kpg1_k_ = reuse_kpg1_k 1018 1019 if(.not.present(ddkinpw) .and. ipert==natom+10) then 1020 ABI_BUG("ddkinpw is not optional for ipert=natom+10.") 1021 end if 1022 if(.not.present(dkinpw2) .and. ipert==natom+10 .and. idir>3) then 1023 ABI_BUG("dkinpw2 is not optional for ipert=natom+10 and idir>3.") 1024 end if 1025 if(.not.present(rf_hamk_dir2) .and. ((ipert==natom+10 .and. idir>3) .or. ipert==natom+11)) then 1026 ABI_BUG("rf_hamk_dir2 is not optional for ipert=natom+10 (with idir>3) or ipert=natom+11.") 1027 end if 1028 1029 ntypat = psps%ntypat 1030 qne0 = ((kpq(1)-kpoint(1))**2+(kpq(2)-kpoint(2))**2+(kpq(3)-kpoint(3))**2>=tol14) 1031 1032 ! Compute k+G vectors 1033 nkpg = 0; if (ipert >= 1 .and. ipert <= natom) nkpg = 3*dtset%nloalg(3) 1034 if (reuse_kpg_k_ == 0) then 1035 ABI_MALLOC(kpg_k, (npw_k, nkpg)) 1036 if (nkpg > 0) call mkkpg(kg_k, kpg_k, kpoint, nkpg, npw_k) 1037 else 1038 ABI_CHECK(all(shape(kpg_k) == [npw_k, nkpg]), "Wrong shape in input kpg_k") 1039 endif 1040 1041 ! Compute k+q+G vectors 1042 nkpg1 = 0; if (ipert >= 1 .and. ipert <= natom) nkpg1 = 3*dtset%nloalg(3) 1043 if (reuse_kpg1_k_ == 0) then 1044 ABI_MALLOC(kpg1_k, (npw1_k, nkpg1)) 1045 if (nkpg1 > 0) call mkkpg(kg1_k, kpg1_k, kpq(:), nkpg1, npw1_k) 1046 else 1047 ABI_CHECK(all(shape(kpg1_k) == [npw1_k, nkpg1]), "Wrong shape in input kpg1_k") 1048 endif 1049 1050 ! ===== Preparation of the non-local contributions 1051 dimffnlk =0; if (ipert<=natom) dimffnlk=1 1052 1053 ! Compute nonlocal form factors ffnlk at (k+G) 1054 ! (only for atomic displacement perturbation) 1055 if (reuse_ffnlk_ == 0) then 1056 ABI_MALLOC(ffnlk, (npw_k, dimffnlk, psps%lmnmax, ntypat)) 1057 if (ipert<=natom) then 1058 ider=0;idir0=0 1059 call mkffnl(psps%dimekb,dimffnlk,psps%ekb,ffnlk,psps%ffspl,& 1060 gmet,gprimd,ider,idir0,psps%indlmn,kg_k,kpg_k,kpoint,psps%lmnmax,& 1061 psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg,npw_k,ntypat,& 1062 psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm_k,ylmgr_dum) 1063 end if 1064 else 1065 ABI_CHECK(all(shape(ffnlk) == [npw_k, dimffnlk, psps%lmnmax, ntypat]), "Wrong shape in input ffnlk") 1066 end if 1067 1068 ! Compute nonlocal form factors ffnl1 at (k+q+G) 1069 !-- Atomic displacement perturbation 1070 if (ipert<=natom) then 1071 ider=0;idir0=0 1072 !-- k-point perturbation (1st-derivative) 1073 else if (ipert==natom+1) then 1074 ider=1;idir0=idir 1075 !-- k-point perturbation (2nd-derivative) 1076 else if (ipert==natom+10.or.ipert==natom+11) then 1077 ider=2;idir0=4 1078 !-- Electric field perturbation 1079 else if (ipert==natom+2) then 1080 if (psps%usepaw==1) then 1081 ider=1;idir0=idir 1082 else 1083 ider=0;idir0=0 1084 end if 1085 !-- Strain perturbation 1086 else if (ipert==natom+3.or.ipert==natom+4) then 1087 if (ipert==natom+3) istr=idir 1088 if (ipert==natom+4) istr=idir+3 1089 ider=1;idir0=-istr 1090 !-- Magnetic field perturbation ( SPr, Zeeman ) 1091 else if(ipert==natom+5)then 1092 ider=0;idir0=0 1093 end if 1094 1095 ! Compute nonlocal form factors ffnl1 at (k+q+G), for all atoms 1096 dimffnl1=1+ider 1097 if (ider==1.and.idir0==0) dimffnl1=2+2*psps%useylm 1098 if (ider==2.and.idir0==4) dimffnl1=3+7*psps%useylm 1099 1100 if (reuse_ffnl1_ == 0) then 1101 ABI_MALLOC(ffnl1, (npw1_k, dimffnl1, psps%lmnmax, ntypat)) 1102 1103 call mkffnl(psps%dimekb,dimffnl1,psps%ekb,ffnl1,psps%ffspl,gmet,gprimd,ider,idir0,& 1104 psps%indlmn,kg1_k,kpg1_k,kpq,psps%lmnmax,psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg1,& 1105 npw1_k,ntypat,psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm1_k,ylmgr1_k) 1106 else 1107 ABI_CHECK(all(shape(ffnl1) == [npw1_k, dimffnl1, psps%lmnmax, ntypat]), "Wrong shape in input ffnl1") 1108 end if 1109 1110 ! Compute ffnl for nonlop with signs = 1 1111 print_info = 0 1112 if (dtset%prtvol==-19.or.dtset%prtvol==-20.or.dtset%prtvol==-21.or.dtset%nonlinear_info>=3) then 1113 print_info = 1 1114 end if 1115 if (present(ffnl1_test).and.print_info/=0.and.(ipert==natom+10.or.ipert==natom+11)) then 1116 ABI_MALLOC(ffnl1_test,(npw1_k,dimffnl1,psps%lmnmax,psps%ntypat)) 1117 idir0 = 0 ! for nonlop with signs = 1 1118 call mkffnl(psps%dimekb,dimffnl1,psps%ekb,ffnl1_test,psps%ffspl,gs_hamkq%gmet,gs_hamkq%gprimd,ider,idir0,& 1119 psps%indlmn,kg1_k,kpg1_k,kpq,psps%lmnmax,psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg1,& 1120 npw1_k,psps%ntypat,psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm1_k,ylmgr1_k) 1121 end if 1122 1123 !===== Preparation of the kinetic contributions 1124 ! Note that not all these arrays should be allocated in the general case when wtk_k vanishes 1125 1126 ! Compute (1/2) (2 Pi)**2 (k+q+G)**2: 1127 ABI_MALLOC(kinpw1, (npw1_k)) 1128 kinpw1(:)=zero 1129 call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg1_k,kinpw1,kpq,npw1_k,0,0) 1130 1131 ABI_MALLOC(dkinpw,(npw_k)) ! 1st derivative (1st direction) 1132 dkinpw(:)=zero 1133 if(ipert==natom+10 .and. idir>3) then 1134 ABI_MALLOC(dkinpw2,(npw_k)) ! 1st derivative (2nd directions) 1135 dkinpw2(:)=zero 1136 end if 1137 if(ipert==natom+10) then 1138 ABI_MALLOC(ddkinpw,(npw_k)) ! 2nd derivative 1139 ddkinpw(:)=zero 1140 end if 1141 1142 ! -- k-point perturbation (1st-derivative) 1143 if (ipert==natom+1) then 1144 ! Compute the derivative of the kinetic operator vs k 1145 call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg_k,dkinpw,kpoint,npw_k,idir,0) ! 1st derivative 1146 end if 1147 1148 !-- k-point perturbation (2nd-derivative) 1149 if (ipert==natom+10.or.ipert==natom+11) then 1150 ! Compute the derivative of the kinetic operator vs k in kinpw, second and first orders 1151 if(ipert==natom+10 .and. idir<=3) then 1152 call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg_k,dkinpw,kpoint,npw_k,idir,0) ! 1st derivative 1153 call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg_k,ddkinpw,kpoint,npw_k,idir,idir) ! 2nd derivative 1154 else 1155 select case(idir) 1156 ! Diagonal terms: 1157 case(1) 1158 idir1 = 1 1159 idir2 = 1 1160 case(2) 1161 idir1 = 2 1162 idir2 = 2 1163 case(3) 1164 idir1 = 3 1165 idir2 = 3 1166 ! Upper triangular terms: 1167 case(4) 1168 idir1 = 2 1169 idir2 = 3 1170 case(5) 1171 idir1 = 1 1172 idir2 = 3 1173 case(6) 1174 idir1 = 1 1175 idir2 = 2 1176 ! Lower triangular terms: 1177 case(7) 1178 idir1 = 3 1179 idir2 = 2 1180 case(8) 1181 idir1 = 3 1182 idir2 = 1 1183 case(9) 1184 idir1 = 2 1185 idir2 = 1 1186 end select 1187 call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg_k,dkinpw,kpoint,npw_k,idir1,0) ! 1st derivative, idir1 1188 if(ipert==natom+10) then 1189 call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg_k,dkinpw2,kpoint,npw_k,idir2,0) ! 1st derivative, idir2 1190 call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg_k,ddkinpw,kpoint,npw_k,idir1,idir2) ! 2nd derivative 1191 end if 1192 end if 1193 end if 1194 1195 !-- Strain perturbation 1196 if (ipert==natom+3.or.ipert==natom+4) then 1197 if (ipert==natom+3) istr=idir 1198 if (ipert==natom+4) istr=idir+3 1199 ! Compute the derivative of the kinetic operator vs strain 1200 call kpgstr(dkinpw,dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,gprimd,istr,kg_k,kpoint,npw_k) 1201 end if 1202 1203 !===== Load the k/k+q dependent parts of the Hamiltonian 1204 ! Load k-dependent part in the Hamiltonian datastructure 1205 ABI_MALLOC(ph3d,(2,npw_k,gs_hamkq%matblk)) 1206 call gs_hamkq%load_k(kpt_k=kpoint,npw_k=npw_k,istwf_k=istwf_k,kg_k=kg_k,kpg_k=kpg_k,& 1207 ph3d_k=ph3d,compute_ph3d=.true.,compute_gbound=.true.) 1208 1209 if (size(ffnlk)>0) then 1210 call gs_hamkq%load_k(ffnl_k=ffnlk) 1211 else 1212 call gs_hamkq%load_k(ffnl_k=ffnl1) 1213 end if 1214 1215 ! Load k+q-dependent part in the Hamiltonian datastructure 1216 ! Note: istwf_k is imposed to 1 for RF calculations (should use istwf_kq instead) 1217 call gs_hamkq%load_kprime(kpt_kp=kpq,npw_kp=npw1_k,istwf_kp=istwf_k,& 1218 kinpw_kp=kinpw1,kg_kp=kg1_k,kpg_kp=kpg1_k,ffnl_kp=ffnl1,compute_gbound=.true.) 1219 1220 if (qne0) then 1221 ABI_MALLOC(ph3d1,(2,npw1_k,gs_hamkq%matblk)) 1222 call gs_hamkq%load_kprime(ph3d_kp=ph3d1,compute_ph3d=.true.) 1223 end if 1224 1225 ! Load k-dependent part in the 1st-order Hamiltonian datastructure 1226 call rf_hamkq%load_k(npw_k=npw_k,dkinpw_k=dkinpw) 1227 1228 if (ipert==natom+10) then 1229 call rf_hamkq%load_k(ddkinpw_k=ddkinpw) 1230 if (idir>3) call rf_hamk_dir2%load_k(dkinpw_k=dkinpw2,ddkinpw_k=ddkinpw) 1231 end if 1232 1233 call timab(195, 2, tsec) 1234 1235 end subroutine getgh1c_setup
m_hamiltonian/getgh1dqc_setup [ Functions ]
[ Top ] [ m_hamiltonian ] [ Functions ]
NAME
getgh1dqc_setup
FUNCTION
INPUTS
OUTPUT
SOURCE
1663 subroutine getgh1dqc_setup(gs_hamkq,rf_hamkq,dtset,psps,kpoint,kpq,idir,ipert,qdir1,& ! In 1664 & natom,rmet,rprimd,gprimd,gmet,istwf_k,npw_k,npw1_k,nylmgr,& ! In 1665 & useylmgr1,kg_k,ylm_k,kg1_k,ylm1_k,ylmgr1_k,& ! In 1666 & nkpg,nkpg1,kpg_k,kpg1_k,dqdqkinpw,kinpw1,ffnlk,ffnl1,ph3d,ph3d1,& ! Out 1667 & reuse_ffnlk,reuse_ffnl1,qdir2) ! Optional 1668 1669 !Arguments ------------------------------------ 1670 !scalars 1671 integer,intent(in) :: idir,ipert,istwf_k,natom,npw_k,npw1_k,nylmgr,qdir1,useylmgr1 1672 integer,intent(in),optional :: reuse_ffnlk,reuse_ffnl1,qdir2 1673 integer,intent(out) :: nkpg,nkpg1 1674 type(gs_hamiltonian_type),intent(inout) :: gs_hamkq 1675 type(rf_hamiltonian_type),intent(inout) :: rf_hamkq 1676 type(dataset_type),intent(in) :: dtset 1677 type(pseudopotential_type),intent(in) :: psps 1678 !arrays 1679 integer,intent(in) :: kg_k(3,npw_k),kg1_k(3,npw1_k) 1680 real(dp),intent(in) :: kpoint(3),kpq(3),gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3) 1681 real(dp),intent(in) :: ylm_k(npw_k,psps%mpsang*psps%mpsang*psps%useylm) 1682 ! real(dp),intent(in) :: ylmgr1_k(npw1_k,3+6*((ipert-natom)/10),psps%mpsang*psps%mpsang*psps%useylm*useylmgr1) 1683 real(dp),intent(in) :: ylmgr1_k(npw1_k,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr1) 1684 real(dp),intent(in) :: ylm1_k(npw1_k,psps%mpsang*psps%mpsang*psps%useylm) 1685 real(dp),allocatable,intent(out) :: dqdqkinpw(:),kinpw1(:) 1686 real(dp),allocatable,intent(inout) :: ffnlk(:,:,:,:),ffnl1(:,:,:,:) 1687 real(dp),allocatable,intent(out) :: kpg_k(:,:),kpg1_k(:,:),ph3d(:,:,:),ph3d1(:,:,:) 1688 1689 !Local variables------------------------------- 1690 !scalars 1691 integer :: dimffnl1,dimffnlk,ider,idir0,ig,mu,mua,mub,ntypat 1692 integer :: nu,nua,nub 1693 integer :: reuse_ffnlk_,reuse_ffnl1_ 1694 logical :: qne0 1695 !arrays 1696 integer,parameter :: alpha(6)=(/1,2,3,3,3,2/),beta(6)=(/1,2,3,2,1,1/) 1697 integer,parameter :: gamma(3,3)=reshape((/1,6,5,6,2,4,5,4,3/),(/3,3/)) 1698 real(dp) :: ylmgr_dum(1,1,1) 1699 real(dp),allocatable :: ffnl1_tmp(:,:,:,:) 1700 1701 1702 ! ************************************************************************* 1703 1704 reuse_ffnlk_ = 0; if (present(reuse_ffnlk)) reuse_ffnlk_ = reuse_ffnlk 1705 reuse_ffnl1_ = 0; if (present(reuse_ffnl1)) reuse_ffnl1_ = reuse_ffnl1 1706 1707 ntypat = psps%ntypat 1708 qne0=((kpq(1)-kpoint(1))**2+(kpq(2)-kpoint(2))**2+(kpq(3)-kpoint(3))**2>=tol14) 1709 1710 !Compute (k+G) vectors 1711 nkpg=0;if(ipert>=1.and.ipert<=natom) nkpg=3*dtset%nloalg(3) 1712 ABI_MALLOC(kpg_k,(npw_k,nkpg)) 1713 if (nkpg>0) then 1714 call mkkpg(kg_k,kpg_k,kpoint,nkpg,npw_k) 1715 end if 1716 1717 !Compute (k+q+G) vectors 1718 nkpg1=0;if(ipert>=1.and.ipert<=natom) nkpg1=3*dtset%nloalg(3) 1719 ABI_MALLOC(kpg1_k,(npw1_k,nkpg1)) 1720 if (nkpg1>0) then 1721 call mkkpg(kg1_k,kpg1_k,kpq(:),nkpg1,npw1_k) 1722 end if 1723 1724 !===== Preparation of the non-local contributions 1725 1726 dimffnlk=0;if (ipert<=natom) dimffnlk=1 1727 1728 !Compute nonlocal form factors ffnlk at (k+G) 1729 if (reuse_ffnlk_ == 0) then 1730 ABI_MALLOC(ffnlk,(npw_k,dimffnlk,psps%lmnmax,ntypat)) 1731 if (ipert<=natom) then 1732 ider=0;idir0=0 1733 call mkffnl(psps%dimekb,dimffnlk,psps%ekb,ffnlk,psps%ffspl,& 1734 & gmet,gprimd,ider,idir0,psps%indlmn,kg_k,kpg_k,kpoint,psps%lmnmax,& 1735 & psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg,npw_k,ntypat,& 1736 & psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm_k,ylmgr_dum) 1737 end if 1738 else 1739 ABI_CHECK(all(shape(ffnlk) == [npw_k, dimffnlk, psps%lmnmax, ntypat]), "Wrong shape in input ffnlk") 1740 end if 1741 1742 !Compute nonlocal form factors ffnl1 at (k+q+G) 1743 !TODO: For the second order gradients, this routine is called for each 3 directions of the 1744 !derivative and every time it calculates all the form factors derivatives. This could be 1745 !done just once. 1746 !-- 1st q-grad of atomic displacement perturbation 1747 if (ipert<=natom.and..not.present(qdir2)) then 1748 ider=1;idir0=qdir1 1749 !-- 2nd q-grad of atomic displacement perturbation 1750 else if (ipert<=natom.and.present(qdir2)) then 1751 ider=2;idir0=4 1752 !-- 2nd q-grad of metric (1st q-grad of strain) perturbation 1753 else if (ipert==natom+3.or.ipert==natom+4) then 1754 ider=2;idir0=4 1755 end if 1756 1757 !Compute nonlocal form factors ffnl1 at (k+q+G), for all atoms 1758 dimffnl1=1+ider 1759 if (ider==2.and.(idir0==0.or.idir0==4)) dimffnl1=3+7*psps%useylm 1760 1761 if (reuse_ffnl1_ == 0) then 1762 ABI_MALLOC(ffnl1,(npw1_k,dimffnl1,psps%lmnmax,ntypat)) 1763 call mkffnl(psps%dimekb,dimffnl1,psps%ekb,ffnl1,psps%ffspl,gmet,gprimd,ider,idir0,& 1764 & psps%indlmn,kg1_k,kpg1_k,kpq,psps%lmnmax,psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg1,& 1765 & npw1_k,ntypat,psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm1_k,ylmgr1_k) 1766 else 1767 ABI_CHECK(all(shape(ffnl1) == [npw1_k, dimffnl1, psps%lmnmax, ntypat]), "Wrong shape in input ffnl1") 1768 end if 1769 1770 1771 !Convert nonlocal form factors to cartesian coordinates. 1772 !For metric (strain) perturbation only. 1773 if (ipert==natom+3.or.ipert==natom+4) then 1774 ABI_MALLOC(ffnl1_tmp,(npw1_k,dimffnl1,psps%lmnmax,ntypat)) 1775 ffnl1_tmp=ffnl1 1776 1777 !First q-derivative 1778 ffnl1(:,2:4,:,:)=zero 1779 do mu=1,3 1780 do ig=1,npw1_k 1781 do nu=1,3 1782 ffnl1(ig,1+mu,:,:)=ffnl1(ig,1+mu,:,:)+ffnl1_tmp(ig,1+nu,:,:)*rprimd(mu,nu) 1783 end do 1784 end do 1785 end do 1786 1787 !Second q-derivative 1788 ffnl1(:,5:10,:,:)=zero 1789 do mu=1,6 1790 mua=alpha(mu);mub=beta(mu) 1791 do ig=1,npw1_k 1792 do nua=1,3 1793 do nub=1,3 1794 nu=gamma(nua,nub) 1795 ffnl1(ig,4+mu,:,:)=ffnl1(ig,4+mu,:,:)+ & 1796 & ffnl1_tmp(ig,4+nu,:,:)*rprimd(mua,nua)*rprimd(mub,nub) 1797 end do 1798 end do 1799 end do 1800 end do 1801 1802 ABI_FREE(ffnl1_tmp) 1803 end if 1804 1805 !===== Preparation of the kinetic contributions 1806 ! Compute (1/2) (2 Pi)**2 (k+q+G)**2: 1807 ABI_MALLOC(kinpw1,(npw1_k)) 1808 kinpw1(:)=zero 1809 call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg1_k,kinpw1,kpq,npw1_k,0,0) 1810 1811 ABI_MALLOC(dqdqkinpw,(npw_k)) 1812 !-- Metric (strain) perturbation 1813 if (ipert==natom+3.or.ipert==natom+4) then 1814 call mkkin_metdqdq(dqdqkinpw,dtset%effmass_free,gprimd,idir,kg_k,kpoint,npw_k,qdir1) 1815 else 1816 dqdqkinpw(:)=zero 1817 end if 1818 1819 !===== Load the k/k+q dependent parts of the Hamiltonian 1820 1821 !Load k-dependent part in the Hamiltonian datastructure 1822 ABI_MALLOC(ph3d,(2,npw_k,gs_hamkq%matblk)) 1823 call gs_hamkq%load_k(kpt_k=kpoint,npw_k=npw_k,istwf_k=istwf_k,kg_k=kg_k,kpg_k=kpg_k,& 1824 & ph3d_k=ph3d,compute_ph3d=.true.,compute_gbound=.true.) 1825 1826 if (size(ffnlk)>0) then 1827 call gs_hamkq%load_k(ffnl_k=ffnlk) 1828 else 1829 call gs_hamkq%load_k(ffnl_k=ffnl1) 1830 end if 1831 1832 !Load k+q-dependent part in the Hamiltonian datastructure 1833 ! Note: istwf_k is imposed to 1 for RF calculations (should use istwf_kq instead) 1834 call gs_hamkq%load_kprime(kpt_kp=kpq,npw_kp=npw1_k,istwf_kp=istwf_k,& 1835 & kinpw_kp=kinpw1,kg_kp=kg1_k,kpg_kp=kpg1_k,ffnl_kp=ffnl1,& 1836 & compute_gbound=.true.) 1837 1838 if (qne0) then 1839 ABI_MALLOC(ph3d1,(2,npw1_k,gs_hamkq%matblk)) 1840 call gs_hamkq%load_kprime(ph3d_kp=ph3d1,compute_ph3d=.true.) 1841 end if 1842 1843 !Load k-dependent part in the 1st-order Hamiltonian datastructure 1844 call rf_hamkq%load_k(npw_k=npw_k,dkinpw_k=dqdqkinpw) 1845 1846 end subroutine getgh1dqc_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
SOURCE
876 subroutine rf_transgrid_and_pack(isppol,nspden,usepaw,cplex,nfftf,nfft,ngfft,nvloc,& 877 & pawfgr,mpi_enreg,vtrial,vtrial1,vlocal,vlocal1) 878 879 !Arguments ------------------------------------ 880 !scalars 881 integer,intent(in) :: isppol,nspden,usepaw,cplex,nfftf,nfft,nvloc 882 type(pawfgr_type),intent(in) :: pawfgr 883 type(MPI_type),intent(in) :: mpi_enreg 884 !arrays 885 integer,intent(in) :: ngfft(18) 886 real(dp),intent(in),target :: vtrial(nfftf,nspden) 887 real(dp),intent(inout),target :: vtrial1(cplex*nfftf,nspden) 888 real(dp),intent(out) :: vlocal(ngfft(4),ngfft(5),ngfft(6),nvloc) 889 real(dp),intent(out) :: vlocal1(cplex*ngfft(4),ngfft(5),ngfft(6),nvloc) 890 891 !Local variables------------------------------- 892 !scalars 893 integer :: n1,n2,n3,n4,n5,n6,paral_kgb,ispden 894 !arrays 895 real(dp) :: rhodum(1) !, tsec(2) 896 real(dp), ABI_CONTIGUOUS pointer :: vtrial_ptr(:,:),vtrial1_ptr(:,:) 897 real(dp),allocatable :: cgrvtrial(:,:),cgrvtrial1(:,:),vlocal_tmp(:,:,:),vlocal1_tmp(:,:,:) 898 899 ! ************************************************************************* 900 901 n1=ngfft(1); n2=ngfft(2); n3=ngfft(3) 902 n4=ngfft(4); n5=ngfft(5); n6=ngfft(6) 903 paral_kgb = mpi_enreg%paral_kgb 904 905 if (nspden/=4) then 906 vtrial_ptr => vtrial 907 if (usepaw==0.or.pawfgr%usefinegrid==0) then 908 call fftpac(isppol,mpi_enreg,nspden,n1,n2,n3,n4,n5,n6,ngfft,vtrial_ptr,vlocal(:,:,:,1),2) 909 call fftpac(isppol,mpi_enreg,nspden,cplex*n1,n2,n3,cplex*n4,n5,n6,ngfft,vtrial1,vlocal1(:,:,:,1),2) 910 else 911 ABI_MALLOC(cgrvtrial,(nfft,nspden)) 912 call transgrid(1,mpi_enreg,nspden,-1,0,0,paral_kgb,pawfgr,rhodum,rhodum,cgrvtrial,vtrial_ptr) 913 call fftpac(isppol,mpi_enreg,nspden,n1,n2,n3,n4,n5,n6,ngfft,cgrvtrial,vlocal(:,:,:,1),2) 914 ABI_REMALLOC(cgrvtrial, (cplex*nfft, nspden)) 915 call transgrid(cplex,mpi_enreg,nspden,-1,0,0,paral_kgb,pawfgr,rhodum,rhodum,cgrvtrial,vtrial1) 916 call fftpac(isppol,mpi_enreg,nspden,cplex*n1,n2,n3,cplex*n4,n5,n6,ngfft,cgrvtrial,vlocal1(:,:,:,1),2) 917 ABI_FREE(cgrvtrial) 918 end if 919 nullify(vtrial_ptr) 920 else 921 ! nspden==4 non-collinear magnetism 922 vtrial_ptr => vtrial 923 vtrial1_ptr => vtrial1 924 ABI_MALLOC(vlocal_tmp,(n4,n5,n6)) 925 ABI_MALLOC(vlocal1_tmp,(cplex*n4,n5,n6)) 926 if (usepaw==0.or.pawfgr%usefinegrid==0) then 927 do ispden=1,nspden 928 call fftpac(ispden,mpi_enreg,nspden,n1,n2,n3,n4,n5,n6,ngfft,vtrial_ptr,vlocal_tmp,2) 929 vlocal(:,:,:,ispden)=vlocal_tmp(:,:,:) 930 call fftpac(ispden,mpi_enreg,nspden,cplex*n1,n2,n3,cplex*n4,n5,n6,ngfft,vtrial1_ptr,vlocal1_tmp,2) 931 vlocal1(:,:,:,ispden)=vlocal1_tmp(:,:,:) 932 end do 933 else 934 ! TODO FR EB check the correctness of the following lines for PAW calculations 935 ABI_MALLOC(cgrvtrial,(nfft,nspden)) 936 ABI_MALLOC(cgrvtrial1,(nfft,nspden)) 937 call transgrid(cplex,mpi_enreg,nspden,-1,0,0,paral_kgb,pawfgr,rhodum,rhodum,cgrvtrial,vtrial_ptr) 938 call transgrid(cplex,mpi_enreg,nspden,-1,0,0,paral_kgb,pawfgr,rhodum,rhodum,cgrvtrial1,vtrial1_ptr) 939 do ispden=1,nspden 940 call fftpac(ispden,mpi_enreg,nspden,n1,n2,n3,n4,n5,n6,ngfft,vtrial_ptr,vlocal_tmp,2) 941 vlocal(:,:,:,ispden)=vlocal_tmp(:,:,:) 942 call fftpac(ispden,mpi_enreg,nspden,n1,n2,n3,n4,n5,n6,ngfft,vtrial1_ptr,vlocal1_tmp,2) 943 vlocal1(:,:,:,ispden)=vlocal1_tmp(:,:,:) 944 end do 945 ABI_FREE(cgrvtrial) 946 end if 947 ABI_FREE(vlocal_tmp) 948 ABI_FREE(vlocal1_tmp) 949 end if !nspden 950 951 end subroutine rf_transgrid_and_pack