TABLE OF CONTENTS
- ABINIT/dfpt_ciflexowf
- ABINIT/dfpt_ddmdqwf
- ABINIT/dfpt_isdqfr
- ABINIT/dfpt_isdqwf
- ABINIT/dfpt_qdrpwf.f90
- ABINIT/m_dfpt_lwwf
ABINIT/dfpt_ciflexowf [ Functions ]
NAME
dfpt_ciflexowf
FUNCTION
This routine computes the band and kpt resolved contributions to the flexoelectric tensor.
INPUTS
cg(2,mpw*nspinor*mband*mkmem*nsppol)=planewave coefficients of wavefunctions at k cplex: if 1, several magnitudes are REAL, if 2, COMPLEX dtset <type(dataset_type)>=all input variables for this dataset gs_hamkq <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k gsqcut=large sphere cut-off icg=shift to be applied on the location of data in the array cg ikpt=number of the k-point indkpt1(nkpt_rbz)=non-symmetrized indices of the k-points isppol=1 for unpolarized, 2 for spin-polarized istwf_k=parameter that describes the storage of wfs kg_k(3,npw_k)=reduced planewave coordinates. kpt(3)=reduced coordinates of k point kxc(nfft,nkxc)=exchange and correlation kernel mkmem =number of k points treated by this node mpi_enreg=information about MPI parallelization mpw=maximum dimensioned size of npw or wfs at k nattyp(ntypat)= # atoms of each type. nband_k=number of bands at this k point for that spin polarization nefipert=number of electric field perturbations nfft=(effective) number of FFT grid points (for this proc) ngfft(1:18)=integer array with FFT box dimensions and other nkpt_rbz= number of k-points in the RBZ nkxc=second dimension of the kxc array. If /=0, the XC kernel must be computed. npw_k=number of plane waves at this k point nq1grad=number of q1 (q_{\gamma}) gradients nq1q2grad=number of q1q2 2nd order gradients nspden=number of spin-density components nsppol=1 for unpolarized, 2 for spin-polarized nstrpert=number of strain perturbations nylmgr=second dimension of ylmgr_k occ_k(nband_k)=occupation number for each band (usually 2) for each k. pert_efield(3,nefipert)=array with the info for the electric field perturbations pert_strain(6,nstrpert)=array with the info for the strain perturbations ph1d(2,3*(2*dtset%mgfft+1)*dtset%natom)=1-dimensional phases psps <type(pseudopotential_type)>=variables related to pseudopotentials q1grad(3,nq1grad)=array with the info for the q1 (q_{\gamma}) gradients q1q2grad(4,nq1q2grad)=array with the info for the q1q2 2nd order gradients rhog(2,nfftf)=array for Fourier transform of GS electron density rhor(nfftf,nspden)=array for GS electron density in electrons/bohr**3. rmet(3,3)=real space metric (bohr**2) ucvol=unit cell volume in bohr**3. useylmgr= if 1 use the derivative of spherical harmonics vhxc1_efield(nefipert,cplex*nfft)= electrostatic potential generated by a first order electric field density vhxc1_strain(nstrpert,cplex*nfft)= electrostatic potential generated by a first order strain density wfk_t_ddk(nq1grad)=unit numbers for the ddk wf1 files wfk_t_dkdk(nq1q2grad)=unit numbers for the dkdk wf1 files wfk_t_efield(nefipert)=unit numbers for the electric field wf1 files wfk_t_strain(3,3)=unit numbers for the strain displacement wf1 files wtk_k=weight assigned to the k point. ylm_k(npw_k,psps%mpsang*psps%mpsang*psps%useylm)=real spherical harmonics for the k point ylmgr_k(npw_k,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr)= k-gradients of real spherical harmonics for the k point
OUTPUT
elflexowf_k(2,3,3,3,3)=wave function dependent part of the electronic flexoelectric tensor for the k-point kpt elflexowf_t1_k(2,3,3,3,3)=t1 term (see notes) of elflexowf_k elflexowf_t2_k(2,3,3,3,3)=t2 term (see notes) of elflexowf_k elflexowf_t3_k(2,3,3,3,3)=t3 term (see notes) of elflexowf_k elflexowf_t4_k(2,3,3,3,3)=t5 term (see notes) of elflexowf_k elflexowf_t5_k(2,3,3,3,3)=t5 term (see notes) of elflexowf_k
SIDE EFFECTS
NOTES
SOURCE
1015 subroutine dfpt_ciflexowf(cg,cplex,dtset,elflexowf_k,elflexowf_t1_k,elflexowf_t2_k,& 1016 & elflexowf_t3_k,elflexowf_t4_k,elflexowf_t5_k, & 1017 & gs_hamkq,gsqcut,icg,ikpt,indkpt1,isppol,istwf_k, & 1018 & kg_k,kpt,kxc,mkmem, & 1019 & mpi_enreg,mpw,nattyp,nband_k,nefipert,nfft,ngfft,nkpt_rbz,nkxc, & 1020 & npw_k,nq1grad, & 1021 & nq1q2grad,nspden,nsppol,nstrpert,nylmgr,occ_k, & 1022 & pert_efield,pert_strain,ph1d,psps,q1grad,q1q2grad,rhog,rhor,rmet,ucvol,useylmgr, & 1023 & vhxc1_efield,vhxc1_strain,wfk_t_efield,wfk_t_ddk, & 1024 & wfk_t_dkdk,wfk_t_strain,wtk_k,ylm_k,ylmgr_k) 1025 1026 !Arguments ------------------------------------ 1027 !scalars 1028 integer,intent(in) :: cplex,icg,ikpt,isppol,istwf_k 1029 integer,intent(in) :: mkmem,mpw,nband_k,nefipert,nfft 1030 integer,intent(in) :: nkpt_rbz,nkxc,npw_k,nq1grad,nq1q2grad,nspden,nsppol,nstrpert,nylmgr 1031 integer,intent(in) :: useylmgr 1032 real(dp),intent(in) :: gsqcut,ucvol,wtk_k 1033 type(dataset_type),intent(in) :: dtset 1034 type(gs_hamiltonian_type),intent(inout) :: gs_hamkq 1035 type(MPI_type),intent(in) :: mpi_enreg 1036 type(pseudopotential_type),intent(in) :: psps 1037 1038 !arrays 1039 integer,intent(in) :: indkpt1(nkpt_rbz),kg_k(3,npw_k),nattyp(dtset%ntypat),ngfft(18) 1040 integer,intent(in) :: pert_efield(3,nefipert),pert_strain(6,nstrpert) 1041 integer,intent(in) :: q1grad(3,nq1grad),q1q2grad(4,nq1q2grad) 1042 real(dp),intent(in) :: cg(2,mpw*dtset%nspinor*dtset%mband*mkmem*nsppol) 1043 real(dp),intent(out) :: elflexowf_k(2,3,3,3,3) 1044 real(dp),intent(out) :: elflexowf_t1_k(2,3,3,3,3),elflexowf_t2_k(2,3,3,3,3) 1045 real(dp),intent(out) :: elflexowf_t3_k(2,3,3,3,3),elflexowf_t4_k(2,3,3,3,3) 1046 real(dp),intent(out) :: elflexowf_t5_k(2,3,3,3,3) 1047 real(dp),intent(in) :: kpt(3),occ_k(nband_k),kxc(nfft,nkxc) 1048 real(dp),intent(in) :: ph1d(2,3*(2*dtset%mgfft+1)*dtset%natom) 1049 real(dp),intent(in) :: rhog(2,nfft),rhor(nfft,nspden),rmet(3,3) 1050 real(dp),intent(in) :: vhxc1_strain(nstrpert,cplex*nfft) 1051 real(dp),intent(in) :: vhxc1_efield(nefipert,cplex*nfft) 1052 real(dp),intent(in) :: ylm_k(npw_k,psps%mpsang*psps%mpsang*psps%useylm) 1053 real(dp),intent(in) :: ylmgr_k(npw_k,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr) 1054 type(wfk_t),intent(inout) :: wfk_t_ddk(nq1grad),wfk_t_dkdk(nq1q2grad) 1055 type(wfk_t),intent(inout) :: wfk_t_efield(nefipert),wfk_t_strain(3,3) 1056 1057 !Local variables------------------------------- 1058 !scalars 1059 integer :: berryopt,g0term 1060 integer :: iband,idir,idirq1,idirq2,iefipert,ii,iq1grad,iq1q2grad,ipert,istr,istrpert,jband,ka,kb 1061 integer :: nkpg,nkpg1,npw_disk 1062 integer :: opt_gvnl1,opthartdqdq,optlocal,optnl 1063 integer ::sij_opt,tim_getgh1c,usevnl,useylmgr1 1064 real(dp) :: cprodi,cprodr,doti,dotr,dum_lambda 1065 character(len=500) :: msg 1066 type(rf_hamiltonian_type) :: rf_hamkq 1067 type(pawfgr_type) :: pawfgr 1068 1069 !arrays 1070 real(dp) :: dum_grad_berry(1,1),dum_gs1(1,1),dum_gvnl1(1,1) 1071 real(dp),allocatable :: cg1_ddk(:,:),cg1_dkdk(:,:),cg1_dkdk_ar(:,:,:),cg1_efield(:,:) 1072 real(dp),allocatable :: cg1_strain(:,:),cg1_strain_ar(:,:,:,:) 1073 real(dp),allocatable :: ci_vefielddag_cj(:,:,:,:),cj_h1vstrain_ci(:,:,:,:) 1074 real(dp),allocatable :: cwave0i(:,:),cwave0j(:,:) 1075 real(dp),allocatable :: c0_VefielddQ_c1strain_bks(:,:,:,:,:) 1076 real(dp),allocatable :: c1dkdk_c1strain_bks(:,:,:,:,:) 1077 real(dp),allocatable :: c1efield_q1gradH0_c1strain_bks(:,:,:,:,:) 1078 real(dp),allocatable :: c1efield_dQcalHstrain_c0_bks(:,:,:,:,:) 1079 real(dp),allocatable :: c1efield_Hmetricdqdq_c0_bks(:,:,:,:,:) 1080 real(dp),allocatable :: dkinpw(:) 1081 real(dp),allocatable :: ffnlk(:,:,:,:),ffnl1(:,:,:,:) 1082 real(dp),allocatable :: gh1dqdqc(:,:),gvloc1dqdqc(:,:),gvnl1dqdqc(:,:),gv1c(:,:) 1083 real(dp),allocatable :: kinpw1(:),kpg_k(:,:),kpg1_k(:,:),ph3d(:,:,:),ph3d1(:,:,:) 1084 real(dp),allocatable :: vhart1dqdq(:) 1085 real(dp),allocatable :: dum_vlocal(:,:,:,:),vlocal1(:,:,:,:),vlocal1dqdq(:,:,:,:),dum_vpsp(:) 1086 real(dp),allocatable :: vpsp1(:),vpsp1dqdq(:),vxc1dqdq(:) 1087 real(dp),allocatable :: dum_ylmgr1_k(:,:,:),part_ylmgr_k(:,:,:) 1088 type(pawcprj_type),allocatable :: dum_cwaveprj(:,:) 1089 1090 ! ************************************************************************* 1091 1092 DBG_ENTER("COLL") 1093 1094 !Not valid for PAW 1095 if (psps%usepaw==1) then 1096 msg=' This routine cannot be used for PAW (use pawnst3 instead) !' 1097 ABI_BUG(msg) 1098 end if 1099 1100 if(dtset%prtvol>2)then 1101 write(msg,'(2a,i5,2x,a,3f9.5)')ch10,' Electronic FxE tensor calculation; k pt #',ikpt,'k=',kpt(:) 1102 call wrtout(std_out,msg) 1103 end if 1104 1105 !Additional definitions 1106 tim_getgh1c=0 1107 1108 !Additional allocations 1109 ABI_MALLOC(cwave0i,(2,npw_k*dtset%nspinor)) 1110 ABI_MALLOC(cwave0j,(2,npw_k*dtset%nspinor)) 1111 ABI_MALLOC(gv1c,(2,npw_k*dtset%nspinor)) 1112 ABI_MALLOC(vlocal1,(cplex*ngfft(4),ngfft(5),ngfft(6),gs_hamkq%nvloc)) 1113 ABI_MALLOC(dum_vpsp,(nfft)) 1114 ABI_MALLOC(dum_vlocal,(ngfft(4),ngfft(5),ngfft(6),gs_hamkq%nvloc)) 1115 ABI_MALLOC(vpsp1,(cplex*nfft)) 1116 ABI_MALLOC(dum_cwaveprj,(0,0)) 1117 1118 !-------------------------------------------------------------------------------------- 1119 !Calculate first terms involving only ground state wavefunctions 1120 !-------------------------------------------------------------------------------------- 1121 1122 !-------------------------------------------------------------------------------------- 1123 !Electric field 1st order potential matrix element: 1124 ! < u_{i,k}^{(0)} | (V^{E_{\alpha}})^{\dagger} | u_{j,k}^{(0)} > = 1125 ! (< u_{j,k}^{(0)} | V^{E_{\alpha}} | u_{i,k}^{(0)} >)^* 1126 !-------------------------------------------------------------------------------------- 1127 1128 !Specific definitions 1129 ipert=dtset%natom+2 1130 useylmgr1=0 1131 dum_lambda=zero 1132 berryopt=0;optlocal=1;optnl=0;usevnl=0;opt_gvnl1=1;sij_opt=0 1133 1134 !Specific allocations 1135 ABI_MALLOC(ci_vefielddag_cj,(2,nefipert,nband_k,nband_k)) 1136 ABI_MALLOC(dum_ylmgr1_k,(npw_k,3+6*((ipert-dtset%natom)/10),psps%mpsang*psps%mpsang*psps%useylm*useylmgr1)) 1137 1138 do iefipert=1,nefipert 1139 ipert=pert_efield(1,iefipert) 1140 idir=pert_efield(2,iefipert) 1141 1142 !Set up local potential vlocal1 with proper dimensioning, from vhxc1_efield 1143 vpsp1=vhxc1_efield(iefipert,:) 1144 call rf_transgrid_and_pack(isppol,nspden,psps%usepaw,cplex,nfft,dtset%nfft,dtset%ngfft,& 1145 & gs_hamkq%nvloc,pawfgr,mpi_enreg,dum_vpsp,vpsp1,dum_vlocal,vlocal1) 1146 1147 !Initializes rf_hamiltonian (the k-dependent part is prepared in getgh1c_setup) 1148 call init_rf_hamiltonian(cplex,gs_hamkq,ipert,rf_hamkq,& 1149 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab) 1150 call rf_hamkq%load_spin(isppol,vlocal1=vlocal1,with_nonlocal=.false.) 1151 1152 !Set up the ground-state Hamiltonian, and some parts of the 1st-order Hamiltonian 1153 call getgh1c_setup(gs_hamkq,rf_hamkq,dtset,psps,& ! In 1154 kpt,kpt,idir,ipert,dtset%natom,rmet,gs_hamkq%gprimd,gs_hamkq%gmet,istwf_k,& ! In 1155 npw_k,npw_k,useylmgr1,kg_k,ylm_k,kg_k,ylm_k,dum_ylmgr1_k,& ! In 1156 dkinpw,nkpg,nkpg1,kpg_k,kpg1_k,kinpw1,ffnlk,ffnl1,ph3d,ph3d1) ! Out 1157 1158 !LOOP OVER KET BANDS 1159 do iband=1,nband_k 1160 1161 if(mpi_enreg%proc_distrb(ikpt,iband,isppol) /= mpi_enreg%me_kpt) cycle 1162 1163 !Read ket ground-state wavefunctions 1164 cwave0i(:,:)=cg(:,1+(iband-1)*npw_k*dtset%nspinor+icg:iband*npw_k*dtset%nspinor+icg) 1165 1166 !Compute <g | V^{E_{\alpha}} | u_{i,k}^{(0)} > 1167 call getgh1c(berryopt,cwave0i,dum_cwaveprj,gv1c,dum_grad_berry,& 1168 & dum_gs1,gs_hamkq,dum_gvnl1,idir,ipert,dum_lambda,mpi_enreg,optlocal,& 1169 & optnl,opt_gvnl1,rf_hamkq,sij_opt,tim_getgh1c,usevnl) 1170 1171 !LOOP OVER BRA BANDS 1172 do jband=1,nband_k 1173 1174 !Read bra ground-state wavefunctions 1175 cwave0j(:,:)=cg(:,1+(jband-1)*npw_k*dtset%nspinor+icg:jband*npw_k*dtset%nspinor+icg) 1176 1177 !Compute ci_vefielddag_cj=(< u_{j,k}^{(0)} | V^{E_{\alpha}} | u_{i,k}^{(0)} >)^* 1178 call dotprod_g(dotr,doti,istwf_k,npw_k*dtset%nspinor,2,cwave0j,gv1c,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 1179 1180 ci_vefielddag_cj(1,iefipert,jband,iband)=dotr 1181 ci_vefielddag_cj(2,iefipert,jband,iband)=-doti 1182 1183 end do !jband 1184 1185 end do !iband 1186 1187 !Clean the electric field rf_hamiltonian 1188 call rf_hamkq%free() 1189 1190 !Deallocations 1191 ABI_FREE(kpg_k) 1192 ABI_FREE(kpg1_k) 1193 ABI_FREE(dkinpw) 1194 ABI_FREE(kinpw1) 1195 ABI_FREE(ffnlk) 1196 ABI_FREE(ffnl1) 1197 ABI_FREE(ph3d) 1198 1199 end do 1200 !End loop over electric field perturbations 1201 1202 !Deallocations 1203 ABI_FREE(dum_ylmgr1_k) 1204 1205 !----------------------------------------------------------------------------------------------- 1206 ! Strain 1st order hamiltonian + 1st order potential matrix element: 1207 ! < u_{j,k}^{(0)} | H^{n_{\beta\delta}}+V^{n_{\beta\delta}} | u_{i,k}^{(0)} > 1208 !----------------------------------------------------------------------------------------------- 1209 1210 !Specific definitions 1211 ipert=dtset%natom+3 1212 useylmgr1=1 1213 dum_lambda=zero 1214 berryopt=0;optlocal=1;optnl=1;usevnl=0;opt_gvnl1=0;sij_opt=0 1215 g0term=1 1216 1217 !Specific allocations 1218 ABI_MALLOC(cj_h1vstrain_ci,(2,nstrpert,nband_k,nband_k)) 1219 ABI_MALLOC(part_ylmgr_k,(npw_k,3+6*((ipert-dtset%natom)/10), psps%mpsang*psps%mpsang*psps%useylm*useylmgr1)) 1220 part_ylmgr_k(:,:,:)=ylmgr_k(:,1:3+6*((ipert-dtset%natom)/10),:) 1221 1222 !LOOP OVER STRAIN PERTURBATIONS 1223 do istrpert=1,nstrpert 1224 ipert=pert_strain(1,istrpert) 1225 idir=pert_strain(2,istrpert) 1226 istr=idir 1227 if(ipert==dtset%natom+4) istr=idir+3 1228 1229 !Get first-order local part of the pseudopotential 1230 call vlocalstr(gs_hamkq%gmet,gs_hamkq%gprimd,gsqcut,istr,dtset%mgfft,mpi_enreg,& 1231 & psps%mqgrid_vl,dtset%natom,nattyp,nfft,ngfft,dtset%ntypat,ph1d,psps%qgrid_vl,& 1232 & ucvol,psps%vlspl,vpsp1,g0term=g0term) 1233 1234 !Set up local potential vlocal1 with proper dimensioning, from vpsp1 + vhxc1_strain 1235 vpsp1=vpsp1+vhxc1_strain(istrpert,:) 1236 call rf_transgrid_and_pack(isppol,nspden,psps%usepaw,cplex,nfft,dtset%nfft,dtset%ngfft,& 1237 & gs_hamkq%nvloc,pawfgr,mpi_enreg,dum_vpsp,vpsp1,dum_vlocal,vlocal1) 1238 1239 !Initialize rf_hamiltonian (the k-dependent part is prepared in getgh1c_setup) 1240 call init_rf_hamiltonian(cplex,gs_hamkq,ipert,rf_hamkq,& 1241 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab) 1242 call rf_hamkq%load_spin(isppol,vlocal1=vlocal1,with_nonlocal=.true.) 1243 1244 !Set up the ground-state Hamiltonian, and some parts of the 1st-order Hamiltonian 1245 call getgh1c_setup(gs_hamkq,rf_hamkq,dtset,psps,& ! In 1246 kpt,kpt,idir,ipert,dtset%natom,rmet,gs_hamkq%gprimd,gs_hamkq%gmet,istwf_k,& ! In 1247 npw_k,npw_k,useylmgr1,kg_k,ylm_k,kg_k,ylm_k,part_ylmgr_k,& ! In 1248 dkinpw,nkpg,nkpg1,kpg_k,kpg1_k,kinpw1,ffnlk,ffnl1,ph3d,ph3d1) ! Out 1249 1250 !LOOP OVER KET BANDS 1251 do iband=1,nband_k 1252 1253 if(mpi_enreg%proc_distrb(ikpt,iband,isppol) /= mpi_enreg%me_kpt) cycle 1254 1255 !Read ket ground-state wavefunctions 1256 cwave0i(:,:)=cg(:,1+(iband-1)*npw_k*dtset%nspinor+icg:iband*npw_k*dtset%nspinor+icg) 1257 1258 !Compute < g | H^{n_{\beta\delta}}+V^{n_{\beta\delta}} | u_{i,k}^{(0)} > 1259 call getgh1c(berryopt,cwave0i,dum_cwaveprj,gv1c,dum_grad_berry,& 1260 & dum_gs1,gs_hamkq,dum_gvnl1,idir,ipert,dum_lambda,mpi_enreg,optlocal,& 1261 & optnl,opt_gvnl1,rf_hamkq,sij_opt,tim_getgh1c,usevnl) 1262 1263 !LOOP OVER BRA BANDS 1264 do jband=1,nband_k 1265 1266 !Read bra ground-state wavefunctions 1267 cwave0j(:,:)=cg(:,1+(jband-1)*npw_k*dtset%nspinor+icg:jband*npw_k*dtset%nspinor+icg) 1268 1269 !Compute = cj_h1vstrain_ci 1270 !< u_{j,k}^{(0)} | H^{n_{\beta\delta}}+V^{n_{\beta\delta}} | u_{i,k}^{(0)} > 1271 call dotprod_g(dotr,doti,istwf_k,npw_k*dtset%nspinor,2,cwave0j,gv1c,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 1272 cj_h1vstrain_ci(1,istrpert,jband,iband)=dotr 1273 cj_h1vstrain_ci(2,istrpert,jband,iband)=doti 1274 1275 end do !jband 1276 1277 end do !iband 1278 1279 !Clean the atomic displacement rf_hamiltonian 1280 call rf_hamkq%free() 1281 1282 !Deallocations 1283 ABI_FREE(kpg_k) 1284 ABI_FREE(kpg1_k) 1285 ABI_FREE(dkinpw) 1286 ABI_FREE(kinpw1) 1287 ABI_FREE(ffnlk) 1288 ABI_FREE(ffnl1) 1289 ABI_FREE(ph3d) 1290 1291 end do !istrpert 1292 !End loop over strain perturbations 1293 1294 !Deallocations 1295 ABI_FREE(cwave0j) 1296 ABI_FREE(vpsp1) 1297 ABI_FREE(vlocal1) 1298 1299 1300 !---------------------------------------------------------------------------------------- 1301 ! Terms that involve first order response functions 1302 !---------------------------------------------------------------------------------------- 1303 1304 !Allocation of bks (band, k-point and spin) dependent terms 1305 ABI_MALLOC(c1efield_q1gradH0_c1strain_bks,(2,nband_k,nefipert,nq1grad,nstrpert)) 1306 ABI_MALLOC(c1efield_dQcalHstrain_c0_bks,(2,nband_k,nefipert,nq1grad,nstrpert)) 1307 ABI_MALLOC(c0_VefielddQ_c1strain_bks,(2,nband_k,nefipert,nq1grad,nstrpert)) 1308 ABI_MALLOC(c1dkdk_c1strain_bks,(2,nband_k,nefipert,nq1grad,nstrpert)) 1309 ABI_MALLOC(c1efield_Hmetricdqdq_c0_bks,(2,nband_k,nefipert,nq1grad,nstrpert)) 1310 c1efield_dQcalHstrain_c0_bks=zero 1311 c0_VefielddQ_c1strain_bks=zero 1312 1313 !Allocation of wf1s 1314 ABI_MALLOC(cg1_strain,(2,npw_k*dtset%nspinor)) 1315 ABI_MALLOC(cg1_strain_ar,(3,3,2,npw_k*dtset%nspinor)) 1316 ABI_MALLOC(cg1_efield,(2,npw_k*dtset%nspinor)) 1317 ABI_MALLOC(cg1_ddk,(2,npw_k*dtset%nspinor)) 1318 ABI_MALLOC(cg1_dkdk,(2,npw_k*dtset%nspinor)) 1319 ABI_MALLOC(cg1_dkdk_ar,(nq1q2grad,2,npw_k*dtset%nspinor)) 1320 1321 1322 !Check correspondance with the data in wf1 files 1323 !Electric field 1324 do iefipert=1,nefipert 1325 !k-point index check 1326 ii = wfk_t_efield(iefipert)%findk(kpt(:)) 1327 ABI_CHECK(ii == indkpt1(ikpt), "ii != indkpt1 in electric field wf1 file") 1328 !npw check 1329 npw_disk = wfk_t_efield(iefipert)%hdr%npwarr(ii) 1330 if (npw_k /= npw_disk) then 1331 write(unit=msg,fmt='(a,i5,a,i5,a,a,i5,a,a,i5)')& 1332 & 'For ikpt = ',ikpt,', and pertcase= ',pert_efield(3,iefipert),ch10,& 1333 & 'the number of plane waves in the wf1 file is equal to', npw_disk,ch10,& 1334 & 'while it should be ',npw_k 1335 ABI_BUG(msg) 1336 end if 1337 end do 1338 1339 !ddk 1340 do iq1grad=1,nq1grad 1341 !k-point index check 1342 ii = wfk_t_ddk(iq1grad)%findk( kpt(:)) 1343 ABI_CHECK(ii == indkpt1(ikpt), "ii != indkpt1 in ddk wf1 file") 1344 !npw check 1345 npw_disk = wfk_t_ddk(iq1grad)%hdr%npwarr(ii) 1346 if (npw_k /= npw_disk) then 1347 write(unit=msg,fmt='(a,i5,a,i5,a,a,i5,a,a,i5)')& 1348 & 'For ikpt = ',ikpt,', and pertcase= ',q1grad(3,iq1grad),ch10,& 1349 & 'the number of plane waves in the wf1 file is equal to', npw_disk,ch10,& 1350 & 'while it should be ',npw_k 1351 ABI_BUG(msg) 1352 end if 1353 end do 1354 1355 !strain 1356 do istrpert=1,nstrpert 1357 ka=pert_strain(3,istrpert) 1358 kb=pert_strain(4,istrpert) 1359 !k-point index check 1360 ii = wfk_t_strain(ka,kb)%findk( kpt(:)) 1361 ABI_CHECK(ii == indkpt1(ikpt), "ii != indkpt1 in strain wf1 file") 1362 !npw check 1363 npw_disk = wfk_t_strain(ka,kb)%hdr%npwarr(ii) 1364 if (npw_k /= npw_disk) then 1365 write(unit=msg,fmt='(a,i5,a,i5,a,a,i5,a,a,i5)')& 1366 & 'For ikpt = ',ikpt,', and pertcase= ',pert_strain(5,istrpert),ch10,& 1367 & 'the number of plane waves in the wf1 file is equal to', npw_disk,ch10,& 1368 & 'while it should be ',npw_k 1369 ABI_BUG(msg) 1370 end if 1371 end do 1372 1373 !dkdk 1374 do iq1q2grad=1,nq1q2grad 1375 !k-point index check 1376 ii = wfk_t_dkdk(iq1q2grad)%findk(kpt(:)) 1377 ABI_CHECK(ii == indkpt1(ikpt), "ii != indkpt1 in dkdk wf1 file") 1378 !npw check 1379 npw_disk = wfk_t_dkdk(iq1q2grad)%hdr%npwarr(ii) 1380 if (npw_k /= npw_disk) then 1381 write(unit=msg,fmt='(a,i5,a,i5,a,a,i5,a,a,i5)')& 1382 & 'For ikpt = ',ikpt,', and pertcase= ',q1q2grad(4,iq1q2grad),ch10,& 1383 & 'the number of plane waves in the wf1 file is equal to', npw_disk,ch10,& 1384 & 'while it should be ',npw_k 1385 ABI_BUG(msg) 1386 end if 1387 end do 1388 1389 !-------------------------------------------------------------------------------------- 1390 !q1-gradient of gs Hamiltonian: 1391 ! < u_{i,k}^{E_{\alpha}} | \partial_{gamma} H^{(0)} | u_{i,k}^{n_{\beta\delta}} > 1392 !-------------------------------------------------------------------------------------- 1393 1394 !Specific allocations 1395 ABI_MALLOC(vlocal1,(cplex*ngfft(4),ngfft(5),ngfft(6),gs_hamkq%nvloc)) 1396 1397 !Specific definitions 1398 vlocal1=zero 1399 useylmgr1=1 1400 dum_lambda=zero 1401 berryopt=0;optlocal=0;optnl=1;usevnl=0;opt_gvnl1=0;sij_opt=0 1402 1403 !LOOP OVER Q1-GRADIENT 1404 do iq1grad=1,nq1grad 1405 ipert=q1grad(1,iq1grad) 1406 idir=q1grad(2,iq1grad) 1407 1408 !Initializes rf_hamiltonian (the k-dependent part is prepared in getgh1c_setup) 1409 call init_rf_hamiltonian(cplex,gs_hamkq,ipert,rf_hamkq,& 1410 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab) 1411 call rf_hamkq%load_spin(isppol,vlocal1=vlocal1,with_nonlocal=.true.) 1412 1413 !Set up the ground-state Hamiltonian, and some parts of the 1st-order Hamiltonian 1414 call getgh1c_setup(gs_hamkq,rf_hamkq,dtset,psps,& ! In 1415 kpt,kpt,idir,ipert,dtset%natom,rmet,gs_hamkq%gprimd,gs_hamkq%gmet,istwf_k,& ! In 1416 npw_k,npw_k,useylmgr1,kg_k,ylm_k,kg_k,ylm_k,part_ylmgr_k,& ! In 1417 dkinpw,nkpg,nkpg1,kpg_k,kpg1_k,kinpw1,ffnlk,ffnl1,ph3d,ph3d1) ! Out 1418 1419 !LOOP OVER BANDS 1420 do iband=1,nband_k 1421 1422 if(mpi_enreg%proc_distrb(ikpt,iband,isppol) /= mpi_enreg%me_kpt) cycle 1423 1424 !LOOP OVER STRAIN PERTURBATION 1425 do istrpert=1,nstrpert 1426 ka=pert_strain(3,istrpert) 1427 kb=pert_strain(4,istrpert) 1428 1429 !Read strain field wf1 1430 if (ka>=kb) then 1431 call wfk_t_strain(ka,kb)%read_bks( iband, indkpt1(ikpt), & 1432 & isppol, xmpio_single, cg_bks=cg1_strain) 1433 cg1_strain_ar(ka,kb,:,:)=cg1_strain 1434 else 1435 cg1_strain=cg1_strain_ar(kb,ka,:,:) 1436 end if 1437 1438 !Compute < g |\partial_{gamma} H^{(0)} | u_{i,k}^{n_{\beta\delta}} > 1439 call getgh1c(berryopt,cg1_strain,dum_cwaveprj,gv1c,dum_grad_berry,& 1440 & dum_gs1,gs_hamkq,dum_gvnl1,idir,ipert,dum_lambda,mpi_enreg,optlocal,& 1441 & optnl,opt_gvnl1,rf_hamkq,sij_opt,tim_getgh1c,usevnl) 1442 1443 !LOOP OVER ELECTRIC FIELD PERTURBATION 1444 do iefipert=1,nefipert 1445 1446 !Read electric field wf1 1447 call wfk_t_efield(iefipert)%read_bks( iband, indkpt1(ikpt), & 1448 & isppol, xmpio_single, cg_bks=cg1_efield) 1449 1450 !calculate: 1451 ! < u_{i,k}^{E_{\alpha}} | \partial_{gamma} H^{(0)} | u_{i,k}^{n_{\beta\delta}} > 1452 call dotprod_g(dotr,doti,istwf_k,npw_k*dtset%nspinor,2,cg1_efield,gv1c, & 1453 & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 1454 1455 c1efield_q1gradH0_c1strain_bks(1,iband,iefipert,iq1grad,istrpert)=dotr 1456 c1efield_q1gradH0_c1strain_bks(2,iband,iefipert,iq1grad,istrpert)=doti 1457 1458 end do !iefipert 1459 1460 end do !istrpert 1461 1462 end do !iband 1463 1464 !Clean the ddk rf_hamiltonian 1465 call rf_hamkq%free() 1466 1467 !Deallocations 1468 ABI_FREE(kpg_k) 1469 ABI_FREE(kpg1_k) 1470 ABI_FREE(dkinpw) 1471 ABI_FREE(kinpw1) 1472 ABI_FREE(ffnlk) 1473 ABI_FREE(ffnl1) 1474 ABI_FREE(ph3d) 1475 1476 end do !iq1grad 1477 1478 !Deallocations 1479 ABI_FREE(gv1c) 1480 ABI_FREE(vlocal1) 1481 !ABI_FREE(ph3d1) !it is only allocated if kpt and kpq are different. Not the case. 1482 ABI_FREE(part_ylmgr_k) 1483 1484 !-------------------------------------------------------------------------------------- 1485 !Other three terms 1486 !-------------------------------------------------------------------------------------- 1487 !LOOP OVER BANDS 1488 do iband=1,nband_k 1489 1490 if(mpi_enreg%proc_distrb(ikpt,iband,isppol) /= mpi_enreg%me_kpt) cycle 1491 1492 !LOOP OVER ELECTRIC FIELD PERTURBATION 1493 do iefipert=1,nefipert 1494 1495 !Read electric field wf1 1496 call wfk_t_efield(iefipert)%read_bks( iband, indkpt1(ikpt), & 1497 & isppol, xmpio_single, cg_bks=cg1_efield) 1498 1499 !LOOP OVER q1-GRADIENT 1500 do iq1grad=1,nq1grad 1501 1502 !LOOP OVER BANDS 1503 do jband=1,nband_k 1504 1505 !Read ddk wf1 1506 call wfk_t_ddk(iq1grad)%read_bks( jband, indkpt1(ikpt), & 1507 & isppol, xmpio_single, cg_bks=cg1_ddk) 1508 1509 !Calculate: < u_{i,k}^{E_{\alpha}} | u_{j,k}^{k_{\gamma}} > 1510 call dotprod_g(dotr,doti,istwf_k,npw_k*dtset%nspinor,2,cg1_efield,cg1_ddk, & 1511 & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 1512 1513 !LOOP OVER STRAIN PERTURBATION 1514 do istrpert=1,nstrpert 1515 1516 !Calculate: -\sum_{j} < u_{i,k}^{E_{\alpha}} | u_{j,k}^{k_{\gamma}} > * 1517 ! < u_{j,k}^{(0)} | H^{n_{\beta\delta}}+V^{n_{\beta\delta}} | u_{i,k}^{(0)} > 1518 cprodr=dotr*cj_h1vstrain_ci(1,istrpert,jband,iband) - & 1519 & doti*cj_h1vstrain_ci(2,istrpert,jband,iband) 1520 cprodi=dotr*cj_h1vstrain_ci(2,istrpert,jband,iband) + & 1521 & doti*cj_h1vstrain_ci(1,istrpert,jband,iband) 1522 1523 c1efield_dQcalHstrain_c0_bks(1,iband,iefipert,iq1grad,istrpert)= & 1524 & c1efield_dQcalHstrain_c0_bks(1,iband,iefipert,iq1grad,istrpert)-cprodr 1525 c1efield_dQcalHstrain_c0_bks(2,iband,iefipert,iq1grad,istrpert)= & 1526 & c1efield_dQcalHstrain_c0_bks(2,iband,iefipert,iq1grad,istrpert)-cprodi 1527 1528 end do !istrpert 1529 1530 end do !jband 1531 1532 end do !iq1grad 1533 1534 end do !iefipert 1535 1536 !LOOP OVER STRAIN PERTURBATION 1537 do istrpert=1,nstrpert 1538 ka=pert_strain(3,istrpert) 1539 kb=pert_strain(4,istrpert) 1540 1541 !Read strain field wf1 1542 if (ka>=kb) then 1543 call wfk_t_strain(ka,kb)%read_bks( iband, indkpt1(ikpt), & 1544 & isppol, xmpio_single, cg_bks=cg1_strain) 1545 cg1_strain_ar(ka,kb,:,:)=cg1_strain 1546 else 1547 cg1_strain=cg1_strain_ar(kb,ka,:,:) 1548 end if 1549 1550 !LOOP OVER q1-GRADIENT 1551 do iq1grad=1,nq1grad 1552 1553 !LOOP OVER BANDS 1554 do jband=1,nband_k 1555 1556 !Read ddk wf1 1557 call wfk_t_ddk(iq1grad)%read_bks( jband, indkpt1(ikpt), & 1558 & isppol, xmpio_single, cg_bks=cg1_ddk) 1559 1560 !Calculate: < u_{j,k}^{k_{\gamma}} | u_{i,k}^{n_{\beta\delta}} > 1561 call dotprod_g(dotr,doti,istwf_k,npw_k*dtset%nspinor,2,cg1_ddk,cg1_strain, & 1562 & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 1563 1564 !LOOP OVER ELECTRIC FIELD PERTURBATION 1565 do iefipert=1,nefipert 1566 !Calculate: -\sum_{j} (< u_{j,k}^{(0)} | V^{E_{\alpha}} | u_{i,k}^{(0)} >)^* x 1567 ! < u_{j,k}^{k_{\gamma}} | u_{i,k}^{n_{\beta\delta}} > 1568 1569 cprodr=dotr*ci_vefielddag_cj(1,iefipert,jband,iband) - & 1570 & doti*ci_vefielddag_cj(2,iefipert,jband,iband) 1571 cprodi=dotr*ci_vefielddag_cj(2,iefipert,jband,iband) + & 1572 & doti*ci_vefielddag_cj(1,iefipert,jband,iband) 1573 1574 c0_VefielddQ_c1strain_bks(1,iband,iefipert,iq1grad,istrpert)= & 1575 & c0_VefielddQ_c1strain_bks(1,iband,iefipert,iq1grad,istrpert)-cprodr 1576 c0_VefielddQ_c1strain_bks(2,iband,iefipert,iq1grad,istrpert)= & 1577 & c0_VefielddQ_c1strain_bks(2,iband,iefipert,iq1grad,istrpert)-cprodi 1578 1579 end do !iefipert 1580 1581 end do !jband 1582 1583 end do !iq1grad 1584 1585 !LOOP OVER q1q2 SECOND ORDER GRADIENT 1586 do iq1q2grad=1,nq1q2grad 1587 1588 !Read dkdk wf1 1589 if (iq1q2grad<=6) then 1590 call wfk_t_dkdk(iq1q2grad)%read_bks( iband, indkpt1(ikpt), & 1591 & isppol, xmpio_single, cg_bks=cg1_dkdk) 1592 cg1_dkdk_ar(iq1q2grad,:,:)=cg1_dkdk 1593 else 1594 cg1_dkdk=cg1_dkdk_ar(iq1q2grad-3,:,:) 1595 end if 1596 1597 !Calculate: -i 1/2 < u_{i,k}^{k_{\alpha},k_{\gamma}} | u_{i,k}^{n_{\beta\delta}} > 1598 call dotprod_g(dotr,doti,istwf_k,npw_k*dtset%nspinor,2,cg1_dkdk,cg1_strain, & 1599 & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 1600 1601 idirq1=q1q2grad(2,iq1q2grad) 1602 idirq2=q1q2grad(3,iq1q2grad) 1603 c1dkdk_c1strain_bks(1,iband,idirq1,idirq2,istrpert)= half*doti 1604 c1dkdk_c1strain_bks(2,iband,idirq1,idirq2,istrpert)=-half*dotr 1605 1606 end do !iq1q2grad 1607 1608 end do !istrpert 1609 1610 end do !iband 1611 1612 !Deallocations 1613 ABI_FREE(ci_vefielddag_cj) 1614 ABI_FREE(cj_h1vstrain_ci) 1615 ABI_FREE(cg1_strain) 1616 ABI_FREE(cg1_strain_ar) 1617 ABI_FREE(cg1_ddk) 1618 ABI_FREE(cg1_dkdk) 1619 ABI_FREE(cg1_dkdk_ar) 1620 1621 !-------------------------------------------------------------------------------------- 1622 ! q1-gradient of the strain first order Hamiltonian: 1623 ! i/2 < u_{i,k}^{E_{\alpha} | H^{n_{\beta\delta}}_{\gamma} | u_{i,k}^{(0)} > 1624 ! Computed as the second order q-gradient of the first order metric Hamiltonian: 1625 ! 1/2 < u_{i,k}^{E_{\alpha} | H^{(\beta)}_{\gamma\delta} | u_{i,k}^{(0)} > 1626 !-------------------------------------------------------------------------------------- 1627 1628 !Specific allocations 1629 ABI_MALLOC(vhart1dqdq,(2*nfft)) 1630 ABI_MALLOC(vpsp1dqdq,(2*nfft)) 1631 ABI_MALLOC(vxc1dqdq,(2*nfft)) 1632 ABI_MALLOC(vlocal1dqdq,(2*ngfft(4),ngfft(5),ngfft(6),gs_hamkq%nvloc)) 1633 ABI_MALLOC(gh1dqdqc,(2,npw_k*dtset%nspinor)) 1634 ABI_MALLOC(gvloc1dqdqc,(2,npw_k*dtset%nspinor)) 1635 ABI_MALLOC(gvnl1dqdqc,(2,npw_k*dtset%nspinor)) 1636 1637 !Specific definitions 1638 useylmgr1=1;optlocal=1;optnl=1;opthartdqdq=1; 1639 1640 !LOOP OVER STRAIN PERTURBATIONS 1641 do istrpert=1,nstrpert 1642 ipert=pert_strain(1,istrpert) 1643 idir=pert_strain(6,istrpert) 1644 1645 !LOOP OVER Q1-GRADIENT 1646 do iq1grad=1,nq1grad 1647 1648 !Get 2nd q-gradient of first-order local part of the pseudopotential and of the Hartree 1649 !(and XC if GGA) contribution from ground state density 1650 call dfpt_vmetdqdq(2,gs_hamkq%gmet,gs_hamkq%gprimd,gsqcut,idir,ipert,kxc,mpi_enreg, & 1651 & psps%mqgrid_vl,dtset%natom, & 1652 & nattyp,nfft,ngfft,dtset%ntypat,ngfft(1),ngfft(2),ngfft(3),nkxc,nspden,opthartdqdq, & 1653 & ph1d,q1grad(2,iq1grad),psps%qgrid_vl,& 1654 & dtset%qptn,rhog,rhor,ucvol,psps%vlspl,vhart1dqdq,vpsp1dqdq,vxc1dqdq) 1655 1656 !Merge the local contributions 1657 vpsp1dqdq=vpsp1dqdq+vhart1dqdq+vxc1dqdq 1658 1659 !Set up q-gradient of strain potential vlocal1dqdq with proper dimensioning 1660 call rf_transgrid_and_pack(isppol,nspden,psps%usepaw,2,nfft,dtset%nfft,dtset%ngfft,& 1661 & gs_hamkq%nvloc,pawfgr,mpi_enreg,dum_vpsp,vpsp1dqdq,dum_vlocal,vlocal1dqdq) 1662 1663 !Initialize rf_hamiltonian (the k-dependent part is prepared in getgh1c_setup) 1664 call init_rf_hamiltonian(2,gs_hamkq,ipert,rf_hamkq,& 1665 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab) 1666 call rf_hamkq%load_spin(isppol,vlocal1=vlocal1dqdq,with_nonlocal=.true.) 1667 1668 !Set up the ground-state Hamiltonian, and some parts of the 1st-order Hamiltonian 1669 call getgh1dqc_setup(gs_hamkq,rf_hamkq,dtset,psps,kpt,kpt,idir,ipert,q1grad(2,iq1grad), & 1670 & dtset%natom,rmet,gs_hamkq%gprimd,gs_hamkq%gmet,istwf_k,npw_k,npw_k,nylmgr,useylmgr1,kg_k, & 1671 & ylm_k,kg_k,ylm_k,ylmgr_k,nkpg,nkpg1,kpg_k,kpg1_k,dkinpw,kinpw1,ffnlk,ffnl1,ph3d,ph3d1) 1672 1673 !LOOP OVER BANDS 1674 do iband=1,nband_k 1675 1676 if(mpi_enreg%proc_distrb(ikpt,iband,isppol) /= mpi_enreg%me_kpt) cycle 1677 1678 !Read ket ground-state wavefunctions 1679 cwave0i(:,:)=cg(:,1+(iband-1)*npw_k*dtset%nspinor+icg:iband*npw_k*dtset%nspinor+icg) 1680 1681 !Compute < g | H^{(\beta)}_{\gamma\delta} | u_{i,k}^{(0)} > 1682 call getgh1dqc(cwave0i,dum_cwaveprj,gh1dqdqc,gvloc1dqdqc,gvnl1dqdqc,gs_hamkq, & 1683 & idir,ipert,mpi_enreg,optlocal,optnl,q1grad(2,iq1grad),rf_hamkq) 1684 1685 !LOOP OVER ELECTRIC FIELD PERTURBATION 1686 do iefipert=1,nefipert 1687 1688 !Read electric field wf1 1689 call wfk_t_efield(iefipert)%read_bks( iband, indkpt1(ikpt), & 1690 & isppol, xmpio_single, cg_bks=cg1_efield) 1691 1692 !calculate: 1693 ! < u_{i,k}^{E_{\alpha}} | H^{(\beta)}_{\gamma\delta} | u_{i,k}^{(0)} > 1694 call dotprod_g(dotr,doti,istwf_k,npw_k*dtset%nspinor,2,cg1_efield,gh1dqdqc, & 1695 & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 1696 1697 c1efield_Hmetricdqdq_c0_bks(1,iband,iefipert,iq1grad,istrpert)=half*dotr 1698 c1efield_Hmetricdqdq_c0_bks(2,iband,iefipert,iq1grad,istrpert)=half*doti 1699 1700 end do !iefipert 1701 1702 end do !iband 1703 1704 !Clean the rf_hamiltonian 1705 call rf_hamkq%free() 1706 1707 !Deallocations 1708 ABI_FREE(kpg_k) 1709 ABI_FREE(kpg1_k) 1710 ABI_FREE(dkinpw) 1711 ABI_FREE(kinpw1) 1712 ABI_FREE(ffnlk) 1713 ABI_FREE(ffnl1) 1714 ABI_FREE(ph3d) 1715 1716 end do !iq1grad 1717 1718 end do !istrpert 1719 1720 !Deallocations 1721 ABI_FREE(dum_cwaveprj) 1722 ABI_FREE(gh1dqdqc) 1723 ABI_FREE(gvloc1dqdqc) 1724 ABI_FREE(gvnl1dqdqc) 1725 ABI_FREE(vpsp1dqdq) 1726 ABI_FREE(vlocal1dqdq) 1727 ABI_FREE(dum_vpsp) 1728 ABI_FREE(dum_vlocal) 1729 ABI_FREE(cwave0i) 1730 ABI_FREE(cg1_efield) 1731 ABI_FREE(vhart1dqdq) 1732 ABI_FREE(vxc1dqdq) 1733 1734 !-------------------------------------------------------------------------------------- 1735 ! Acumulates all the wf dependent terms of the flexoelectric tensor 1736 !-------------------------------------------------------------------------------------- 1737 elflexowf_k=zero 1738 elflexowf_t1_k=zero 1739 elflexowf_t2_k=zero 1740 elflexowf_t3_k=zero 1741 elflexowf_t4_k=zero 1742 elflexowf_t5_k=zero 1743 do istrpert=1,nstrpert 1744 ka=pert_strain(3,istrpert) 1745 kb=pert_strain(4,istrpert) 1746 do iq1grad=1,nq1grad 1747 do iefipert=1,nefipert 1748 do iband=1,nband_k 1749 1750 if(mpi_enreg%proc_distrb(ikpt,iband,isppol) /= mpi_enreg%me_kpt) cycle 1751 1752 !All terms toghether except T4 that needs further treatment 1753 elflexowf_k(1,iefipert,iq1grad,ka,kb)=elflexowf_k(1,iefipert,iq1grad,ka,kb) + & 1754 & occ_k(iband) * ( c1efield_q1gradH0_c1strain_bks(1,iband,iefipert,iq1grad,istrpert) + & !T1 1755 & c1efield_dQcalHstrain_c0_bks(1,iband,iefipert,iq1grad,istrpert) + & !T2 1756 & c0_VefielddQ_c1strain_bks(1,iband,iefipert,iq1grad,istrpert) + & !T3 1757 & c1dkdk_c1strain_bks(1,iband,iefipert,iq1grad,istrpert) ) !T5 1758 1759 elflexowf_k(2,iefipert,iq1grad,ka,kb)=elflexowf_k(2,iefipert,iq1grad,ka,kb) + & 1760 & occ_k(iband) * ( c1efield_q1gradH0_c1strain_bks(2,iband,iefipert,iq1grad,istrpert) + & !T1 1761 & c1efield_dQcalHstrain_c0_bks(2,iband,iefipert,iq1grad,istrpert) + & !T2 1762 & c0_VefielddQ_c1strain_bks(2,iband,iefipert,iq1grad,istrpert) + & !T3 1763 & c1dkdk_c1strain_bks(2,iband,iefipert,iq1grad,istrpert) ) !T5 1764 1765 !Separate them 1766 !T1 1767 elflexowf_t1_k(1,iefipert,iq1grad,ka,kb)=elflexowf_t1_k(1,iefipert,iq1grad,ka,kb) + & 1768 & occ_k(iband) * c1efield_q1gradH0_c1strain_bks(1,iband,iefipert,iq1grad,istrpert) 1769 1770 elflexowf_t1_k(2,iefipert,iq1grad,ka,kb)=elflexowf_t1_k(2,iefipert,iq1grad,ka,kb) + & 1771 & occ_k(iband) * c1efield_q1gradH0_c1strain_bks(2,iband,iefipert,iq1grad,istrpert) 1772 1773 !T2 1774 elflexowf_t2_k(1,iefipert,iq1grad,ka,kb)=elflexowf_t2_k(1,iefipert,iq1grad,ka,kb) + & 1775 & occ_k(iband) * c1efield_dQcalHstrain_c0_bks(1,iband,iefipert,iq1grad,istrpert) 1776 1777 elflexowf_t2_k(2,iefipert,iq1grad,ka,kb)=elflexowf_t2_k(2,iefipert,iq1grad,ka,kb) + & 1778 & occ_k(iband) * c1efield_dQcalHstrain_c0_bks(2,iband,iefipert,iq1grad,istrpert) 1779 1780 !T3 1781 elflexowf_t3_k(1,iefipert,iq1grad,ka,kb)=elflexowf_t3_k(1,iefipert,iq1grad,ka,kb) + & 1782 & occ_k(iband) * c0_VefielddQ_c1strain_bks(1,iband,iefipert,iq1grad,istrpert) 1783 1784 elflexowf_t3_k(2,iefipert,iq1grad,ka,kb)=elflexowf_t3_k(2,iefipert,iq1grad,ka,kb) + & 1785 & occ_k(iband) * c0_VefielddQ_c1strain_bks(2,iband,iefipert,iq1grad,istrpert) 1786 1787 !T4 has type-I ordering for the array indexes 1788 elflexowf_t4_k(1,iefipert,ka,kb,iq1grad)=elflexowf_t4_k(1,iefipert,ka,kb,iq1grad) + & 1789 & occ_k(iband) * c1efield_Hmetricdqdq_c0_bks(1,iband,iefipert,iq1grad,istrpert) 1790 1791 elflexowf_t4_k(2,iefipert,ka,kb,iq1grad)=elflexowf_t4_k(2,iefipert,ka,kb,iq1grad) + & 1792 & occ_k(iband) * c1efield_Hmetricdqdq_c0_bks(2,iband,iefipert,iq1grad,istrpert) 1793 1794 !T5 1795 elflexowf_t5_k(1,iefipert,iq1grad,ka,kb)=elflexowf_t5_k(1,iefipert,iq1grad,ka,kb) + & 1796 & occ_k(iband) * c1dkdk_c1strain_bks(1,iband,iefipert,iq1grad,istrpert) 1797 1798 elflexowf_t5_k(2,iefipert,iq1grad,ka,kb)=elflexowf_t5_k(2,iefipert,iq1grad,ka,kb) + & 1799 & occ_k(iband) * c1dkdk_c1strain_bks(2,iband,iefipert,iq1grad,istrpert) 1800 1801 end do 1802 end do 1803 end do 1804 end do 1805 1806 !scale by the k-point weight 1807 elflexowf_t1_k=elflexowf_t1_k * wtk_k 1808 elflexowf_t2_k=elflexowf_t2_k * wtk_k 1809 elflexowf_t3_k=elflexowf_t3_k * wtk_k 1810 elflexowf_t4_k=elflexowf_t4_k * wtk_k 1811 elflexowf_t5_k=elflexowf_t5_k * wtk_k 1812 elflexowf_k=elflexowf_k * wtk_k 1813 1814 !Deallocations 1815 ABI_FREE(c1efield_q1gradH0_c1strain_bks) 1816 ABI_FREE(c1efield_dQcalHstrain_c0_bks) 1817 ABI_FREE(c0_VefielddQ_c1strain_bks) 1818 ABI_FREE(c1efield_Hmetricdqdq_c0_bks) 1819 ABI_FREE(c1dkdk_c1strain_bks) 1820 1821 1822 DBG_EXIT("COLL") 1823 1824 end subroutine dfpt_ciflexowf
ABINIT/dfpt_ddmdqwf [ Functions ]
NAME
dfpt_ddmdqwf
FUNCTION
This routine computes the band and kpt resolved contributions to the first q-gradient of the dynamical matrix.
INPUTS
atindx(natom)=index table for atoms (see gstate.f) cg(2,mpw*nspinor*mband*mkmem*nsppol)=planewave coefficients of wavefunctions at k cplex: if 1, several magnitudes are REAL, if 2, COMPLEX dtset <type(dataset_type)>=all input variables for this dataset gs_hamkq <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k gsqcut=large sphere cut-off icg=shift to be applied on the location of data in the array cg ikpt=number of the k-point indkpt1(nkpt_rbz)=non-symmetrized indices of the k-points isppol=1 for unpolarized, 2 for spin-polarized istwf_k=parameter that describes the storage of wfs kg_k(3,npw_k)=reduced planewave coordinates. kpt(3)=reduced coordinates of k point mkmem =number of k points treated by this node mpi_enreg=information about MPI parallelization mpw=maximum dimensioned size of npw or wfs at k natpert=number of atomic displacement perturbations nattyp(ntypat)= # atoms of each type. nband_k=number of bands at this k point for that spin polarization nfft=(effective) number of FFT grid points (for this proc) ngfft(1:18)=integer array with FFT box dimensions and other nkpt_rbz= number of k-points in the RBZ npw_k=number of plane waves at this k point nq1grad=number of q1 (q_{\gamma}) gradients nspden=number of spin-density components nsppol=1 for unpolarized, 2 for spin-polarized nylmgr=second dimension of ylmgr_k occ_k(nband_k)=occupation number for each band (usually 2) for each k. pert_atdis(3,natpert)=array with the info for the atomic displacement perturbations ph1d(2,3*(2*dtset%mgfft+1)*dtset%natom)=1-dimensional phases psps <type(pseudopotential_type)>=variables related to pseudopotentials q1grad(3,nq1grad)=array with the info for the q1 (q_{\gamma}) gradients rmet(3,3)=real space metric (bohr**2) ucvol=unit cell volume in bohr**3. useylmgr= if 1 use the derivative of spherical harmonics vhxc1_atdis(natdis,cplex*nfft)= electrostatic potential generated by a first order atomic displacement density wfk_t_ddk(nq1grad)=unit numbers for the ddk wf1 files wfk_t_atdis(natpert)=unit numbers for the atomic displacement wf1 files wtk_k=weight assigned to the k point. xred(3,natom)=reduced dimensionless atomic coordinates ylm_k(npw_k,psps%mpsang*psps%mpsang*psps%useylm)=real spherical harmonics for the k point ylmgr_k(npw_k,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr)= k-gradients of real spherical harmonics for the k point
OUTPUT
ddmdqwf_k(2,natpert,natpert,nq1grad)=wave function dependent part of the first q-gradient of the dynamical matrix for the k-point kpt. ddmdqwf_t1_k(2,natpert,natpert,nq1grad)=t1 term (see notes) of ddmdqwf_k ddmdqwf_t2_k(2,natpert,natpert,nq1grad)=t2 (cc of t4) term (see notes) of ddmdqwf_k ddmdqwf_t3_k(2,natpert,natpert,nq1grad)=t3 (cc of t5) term (see notes) of ddmdqwf_k
SIDE EFFECTS
NOTES
SOURCE
1895 subroutine dfpt_ddmdqwf(atindx,cg,cplex,ddmdqwf_k,ddmdqwf_t1_k,ddmdqwf_t2_k,& 1896 & ddmdqwf_t3_k,dtset, & 1897 & gs_hamkq,gsqcut,icg,ikpt,indkpt1,isppol,istwf_k, & 1898 & kg_k,kpt,mkmem, & 1899 & mpi_enreg,mpw,natpert,nattyp,nband_k,nfft,ngfft,nkpt_rbz, & 1900 & npw_k,nq1grad,nspden,nsppol,nylmgr,occ_k, & 1901 & pert_atdis,ph1d,psps,q1grad,rmet,ucvol,useylmgr, & 1902 & vhxc1_atdis,wfk_t_atdis,wfk_t_ddk, & 1903 & wtk_k,xred,ylm_k,ylmgr_k) 1904 1905 !Arguments ------------------------------------ 1906 !scalars 1907 integer,intent(in) :: cplex,icg,ikpt,isppol,istwf_k 1908 integer,intent(in) :: mkmem,mpw,natpert,nband_k,nfft 1909 integer,intent(in) :: nkpt_rbz,npw_k,nq1grad,nspden,nsppol,nylmgr 1910 integer,intent(in) :: useylmgr 1911 real(dp),intent(in) :: gsqcut,ucvol,wtk_k 1912 type(dataset_type),intent(in) :: dtset 1913 type(gs_hamiltonian_type),intent(inout) :: gs_hamkq 1914 type(MPI_type),intent(in) :: mpi_enreg 1915 type(pseudopotential_type),intent(in) :: psps 1916 1917 !arrays 1918 integer,intent(in) :: atindx(dtset%natom) 1919 integer,intent(in) :: indkpt1(nkpt_rbz),kg_k(3,npw_k),nattyp(dtset%ntypat),ngfft(18) 1920 integer,intent(in) :: pert_atdis(3,natpert) 1921 integer,intent(in) :: q1grad(3,nq1grad) 1922 real(dp),intent(in) :: cg(2,mpw*dtset%nspinor*dtset%mband*mkmem*nsppol) 1923 real(dp),intent(out) :: ddmdqwf_k(2,natpert,natpert,nq1grad) 1924 real(dp),intent(out) :: ddmdqwf_t1_k(2,natpert,natpert,nq1grad) 1925 real(dp),intent(out) :: ddmdqwf_t2_k(2,natpert,natpert,nq1grad) 1926 real(dp),intent(out) :: ddmdqwf_t3_k(2,natpert,natpert,nq1grad) 1927 real(dp),intent(in) :: kpt(3),occ_k(nband_k) 1928 real(dp),intent(in) :: ph1d(2,3*(2*dtset%mgfft+1)*dtset%natom) 1929 real(dp),intent(in) :: rmet(3,3) 1930 real(dp),intent(in) :: vhxc1_atdis(natpert,cplex*nfft) 1931 real(dp),intent(in) :: xred(3,dtset%natom) 1932 real(dp),intent(in) :: ylm_k(npw_k,psps%mpsang*psps%mpsang*psps%useylm) 1933 real(dp),intent(in) :: ylmgr_k(npw_k,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr) 1934 type(wfk_t),intent(inout) :: wfk_t_ddk(nq1grad),wfk_t_atdis(natpert) 1935 1936 !Local variables------------------------------- 1937 !scalars 1938 integer :: berryopt,iatpert,iband,idir,ii,ipert,iq1grad 1939 integer :: jatpert,jband,jdir,jpert,nkpg,nkpg1,npw_disk,nylmgreff 1940 integer :: opt_gvnl1,optlocal,optnl,sij_opt,tim_getgh1c,usevnl,useylmgr1 1941 real(dp) :: cprodi,cprodr,doti,dotr,dum_lambda 1942 character(len=500) :: msg 1943 type(pawfgr_type) :: pawfgr 1944 type(rf_hamiltonian_type) :: rf_hamkq 1945 !arrays 1946 real(dp) :: dum_grad_berry(1,1),dum_gs1(1,1),dum_gvnl1(1,1) 1947 real(dp),allocatable :: c1atdis_dQHatdis_c0_bks(:,:,:,:,:) 1948 real(dp),allocatable :: c1atdis_Hatdisdq_c0_bks(:,:,:,:,:) 1949 real(dp),allocatable :: c1atdis_q1gradH0_c1atdis_bks(:,:,:,:,:) 1950 real(dp),allocatable :: cg1_iatdis(:,:),cg1_jatdis(:,:),cg1_ddk(:,:) 1951 real(dp),allocatable :: ci_h1vatdis_cj(:,:,:,:) 1952 real(dp),allocatable :: cwave0i(:,:),cwave0j(:,:) 1953 real(dp),allocatable :: dkinpw(:) 1954 real(dp),allocatable :: ffnlk(:,:,:,:),ffnl1(:,:,:,:) 1955 real(dp),allocatable :: gh1dqc(:,:),gvloc1dqc(:,:),gvnl1dqc(:,:),gv1c(:,:) 1956 real(dp),allocatable :: kinpw1(:),kpg_k(:,:),kpg1_k(:,:),ph3d(:,:,:),ph3d1(:,:,:) 1957 real(dp),allocatable :: dum_vlocal(:,:,:,:),vlocal1(:,:,:,:),vlocal1dq(:,:,:,:), dum_vpsp(:) 1958 real(dp),allocatable :: vpsp1(:),vpsp1dq(:), dum_ylmgr1_k(:,:,:),part_ylmgr_k(:,:,:) 1959 type(pawcprj_type),allocatable :: dum_cwaveprj(:,:) 1960 1961 1962 ! ************************************************************************* 1963 1964 DBG_ENTER("COLL") 1965 1966 if(dtset%prtvol>2)then 1967 write(msg,'(2a,i5,2x,a,3f9.5)')ch10,' First q-moment IFCs calculation; k pt #',ikpt,'k=',kpt(:) 1968 call wrtout(std_out,msg) 1969 end if 1970 1971 !Additional definitions 1972 tim_getgh1c=0 1973 1974 !Additional allocations 1975 ABI_MALLOC(cwave0i,(2,npw_k*dtset%nspinor)) 1976 ABI_MALLOC(cwave0j,(2,npw_k*dtset%nspinor)) 1977 ABI_MALLOC(gv1c,(2,npw_k*dtset%nspinor)) 1978 ABI_MALLOC(vlocal1,(cplex*ngfft(4),ngfft(5),ngfft(6),gs_hamkq%nvloc)) 1979 ABI_MALLOC(dum_vpsp,(nfft)) 1980 ABI_MALLOC(dum_vlocal,(ngfft(4),ngfft(5),ngfft(6),gs_hamkq%nvloc)) 1981 ABI_MALLOC(vpsp1,(cplex*nfft)) 1982 ABI_MALLOC(dum_cwaveprj,(0,0)) 1983 1984 !----------------------------------------------------------------------------------------------- 1985 ! Atomic displacement 1st order hamiltonian + 1st order potential matrix element: 1986 ! < u_{i,k}^{(0)} | H^{\tau_{\kappa\beta}}+V^{\tau_{\kappa\beta}} | u_{j,k}^{(0)} > 1987 !----------------------------------------------------------------------------------------------- 1988 1989 !Specific definitions 1990 ipert=1 1991 useylmgr1=0 1992 dum_lambda=zero 1993 berryopt=0;optlocal=1;optnl=1;usevnl=0;opt_gvnl1=0;sij_opt=0 1994 1995 !Specific allocations 1996 ABI_MALLOC(ci_h1vatdis_cj,(2,natpert,nband_k,nband_k)) 1997 ABI_MALLOC(dum_ylmgr1_k,(npw_k,3+6*((ipert-dtset%natom)/10), psps%mpsang*psps%mpsang*psps%useylm*useylmgr1)) 1998 1999 !LOOP OVER ATOMIC DISPLACEMENT PERTURBATIONS 2000 do iatpert= 1, natpert 2001 ipert=pert_atdis(1,iatpert) 2002 idir=pert_atdis(2,iatpert) 2003 2004 !Get first-order local part of the pseudopotential 2005 call dfpt_vlocal(atindx,cplex,gs_hamkq%gmet,gsqcut,idir,ipert,mpi_enreg, & 2006 & psps%mqgrid_vl,dtset%natom,& 2007 & nattyp,nfft,ngfft,dtset%ntypat,ngfft(1),ngfft(2),ngfft(3), & 2008 & ph1d,psps%qgrid_vl,& 2009 & dtset%qptn,ucvol,psps%vlspl,vpsp1,xred) 2010 2011 !Set up local potential vlocal1 with proper dimensioning, from vpsp1 + vhxc1_atdis 2012 vpsp1=vpsp1+vhxc1_atdis(iatpert,:) 2013 call rf_transgrid_and_pack(isppol,nspden,psps%usepaw,cplex,nfft,dtset%nfft,dtset%ngfft,& 2014 & gs_hamkq%nvloc,pawfgr,mpi_enreg,dum_vpsp,vpsp1,dum_vlocal,vlocal1) 2015 2016 !Initialize rf_hamiltonian (the k-dependent part is prepared in getgh1c_setup) 2017 call init_rf_hamiltonian(cplex,gs_hamkq,ipert,rf_hamkq,& 2018 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab) 2019 call rf_hamkq%load_spin(isppol,vlocal1=vlocal1,with_nonlocal=.true.) 2020 2021 !Set up the ground-state Hamiltonian, and some parts of the 1st-order Hamiltonian 2022 call getgh1c_setup(gs_hamkq,rf_hamkq,dtset,psps,& ! In 2023 kpt,kpt,idir,ipert,dtset%natom,rmet,gs_hamkq%gprimd,gs_hamkq%gmet,istwf_k,& ! In 2024 npw_k,npw_k,useylmgr1,kg_k,ylm_k,kg_k,ylm_k,dum_ylmgr1_k,& ! In 2025 dkinpw,nkpg,nkpg1,kpg_k,kpg1_k,kinpw1,ffnlk,ffnl1,ph3d,ph3d1) ! Out 2026 2027 !LOOP OVER KET BANDS 2028 do iband=1,nband_k 2029 2030 if(mpi_enreg%proc_distrb(ikpt,iband,isppol) /= mpi_enreg%me_kpt) cycle 2031 2032 !Read ket ground-state wavefunctions 2033 cwave0i(:,:)=cg(:,1+(iband-1)*npw_k*dtset%nspinor+icg:iband*npw_k*dtset%nspinor+icg) 2034 2035 !Compute < g | H^{\tau_{\kappa\beta}}+V^{\tau_{\kappa\beta}} | u_{i,k}^{(0)} > 2036 call getgh1c(berryopt,cwave0i,dum_cwaveprj,gv1c,dum_grad_berry,& 2037 & dum_gs1,gs_hamkq,dum_gvnl1,idir,ipert,dum_lambda,mpi_enreg,optlocal,& 2038 & optnl,opt_gvnl1,rf_hamkq,sij_opt,tim_getgh1c,usevnl) 2039 2040 !LOOP OVER BRA BANDS 2041 do jband=1,nband_k 2042 2043 !Read bra ground-state wavefunctions 2044 cwave0j(:,:)=cg(:,1+(jband-1)*npw_k*dtset%nspinor+icg:jband*npw_k*dtset%nspinor+icg) 2045 2046 !Compute = ci_h1vatdisdag_cj 2047 !(< u_{j,k}^{(0)} | H^{\tau_{\kappa\beta}}+V^{\tau_{\kappa\beta}} | u_{i,k}^{(0)} >)^* 2048 call dotprod_g(dotr,doti,istwf_k,npw_k*dtset%nspinor,2,cwave0j,gv1c,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 2049 ci_h1vatdis_cj(1,iatpert,jband,iband)=dotr 2050 ci_h1vatdis_cj(2,iatpert,jband,iband)=doti 2051 2052 end do !iband 2053 2054 end do !jband 2055 2056 !Clean the atomic displacement rf_hamiltonian 2057 call rf_hamkq%free() 2058 2059 !Deallocations 2060 ABI_FREE(kpg_k) 2061 ABI_FREE(kpg1_k) 2062 ABI_FREE(dkinpw) 2063 ABI_FREE(kinpw1) 2064 ABI_FREE(ffnlk) 2065 ABI_FREE(ffnl1) 2066 ABI_FREE(ph3d) 2067 2068 end do !iatpert 2069 2070 !Deallocations 2071 ABI_FREE(dum_ylmgr1_k) 2072 ABI_FREE(cwave0j) 2073 ABI_FREE(vpsp1) 2074 ABI_FREE(vlocal1) 2075 2076 !---------------------------------------------------------------------------------------- 2077 ! Terms that involve first order response functions 2078 !---------------------------------------------------------------------------------------- 2079 !Allocation of bks (band, k-point and spin) dependent terms 2080 ABI_MALLOC(c1atdis_q1gradH0_c1atdis_bks,(2,nband_k,natpert,natpert,nq1grad)) 2081 ABI_MALLOC(c1atdis_dQHatdis_c0_bks,(2,nband_k,natpert,natpert,nq1grad)) 2082 ABI_MALLOC(c1atdis_Hatdisdq_c0_bks,(2,nband_k,natpert,natpert,nq1grad)) 2083 c1atdis_dQHatdis_c0_bks=zero 2084 2085 !Allocation of wf1s 2086 ABI_MALLOC(cg1_iatdis,(2,npw_k*dtset%nspinor)) 2087 ABI_MALLOC(cg1_jatdis,(2,npw_k*dtset%nspinor)) 2088 ABI_MALLOC(cg1_ddk,(2,npw_k*dtset%nspinor)) 2089 2090 !Check correspondance with the data in wf1 files 2091 !Atomic displacements 2092 do iatpert=1,natpert 2093 !k-point index check 2094 ii = wfk_t_atdis(iatpert)%findk( kpt(:)) 2095 ABI_CHECK(ii == indkpt1(ikpt), "ii != indkpt in atomic displacement wf1 file") 2096 !npw check 2097 npw_disk = wfk_t_atdis(iatpert)%hdr%npwarr(ii) 2098 if (npw_k /= npw_disk) then 2099 write(unit=msg,fmt='(a,i5,a,i5,a,a,i5,a,a,i5)')& 2100 & 'For ikpt = ',ikpt,', and pertcase= ',pert_atdis(3,iatpert),ch10,& 2101 & 'the number of plane waves in the wf1 file is equal to', npw_disk,ch10,& 2102 & 'while it should be ',npw_k 2103 ABI_BUG(msg) 2104 end if 2105 end do 2106 2107 !ddk 2108 do iq1grad=1,nq1grad 2109 !k-point index check 2110 ii = wfk_t_ddk(iq1grad)%findk( kpt(:)) 2111 ABI_CHECK(ii == indkpt1(ikpt), "ii != indkpt1 in ddk wf1 file") 2112 !npw check 2113 npw_disk = wfk_t_ddk(iq1grad)%hdr%npwarr(ii) 2114 if (npw_k /= npw_disk) then 2115 write(unit=msg,fmt='(a,i5,a,i5,a,a,i5,a,a,i5)')& 2116 & 'For ikpt = ',ikpt,', and pertcase= ',q1grad(3,iq1grad),ch10,& 2117 & 'the number of plane waves in the wf1 file is equal to', npw_disk,ch10,& 2118 & 'while it should be ',npw_k 2119 ABI_BUG(msg) 2120 end if 2121 end do 2122 2123 2124 !-------------------------------------------------------------------------------------- 2125 !q1-gradient of gs Hamiltonian: 2126 ! < u_{i,k}^{\tau_{\kappa\alpha}} | \partial_{gamma} H^{(0)} | u_{i,k}^{\tau_{\kappa'\beta}} > 2127 !-------------------------------------------------------------------------------------- 2128 2129 !Specific allocations 2130 ipert=dtset%natom+1 2131 ABI_MALLOC(vlocal1,(cplex*ngfft(4),ngfft(5),ngfft(6),gs_hamkq%nvloc)) 2132 ABI_MALLOC(part_ylmgr_k,(npw_k,3+6*((ipert-dtset%natom)/10), psps%mpsang*psps%mpsang*psps%useylm*useylmgr)) 2133 part_ylmgr_k(:,:,:)=ylmgr_k(:,1:3+6*((ipert-dtset%natom)/10),:) 2134 2135 !Specific definitions 2136 vlocal1=zero 2137 useylmgr1=1 2138 dum_lambda=zero 2139 berryopt=0;optlocal=0;optnl=1;usevnl=0;opt_gvnl1=0;sij_opt=0 2140 2141 !LOOP OVER Q1-GRADIENT 2142 do iq1grad=1,nq1grad 2143 ipert=q1grad(1,iq1grad) 2144 idir=q1grad(2,iq1grad) 2145 2146 !Initializes rf_hamiltonian (the k-dependent part is prepared in getgh1c_setup) 2147 call init_rf_hamiltonian(cplex,gs_hamkq,ipert,rf_hamkq,& 2148 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab) 2149 call rf_hamkq%load_spin(isppol,vlocal1=vlocal1,with_nonlocal=.true.) 2150 2151 !Set up the ground-state Hamiltonian, and some parts of the 1st-order Hamiltonian 2152 call getgh1c_setup(gs_hamkq,rf_hamkq,dtset,psps,& ! In 2153 kpt,kpt,idir,ipert,dtset%natom,rmet,gs_hamkq%gprimd,gs_hamkq%gmet,istwf_k,& ! In 2154 npw_k,npw_k,useylmgr1,kg_k,ylm_k,kg_k,ylm_k,part_ylmgr_k,& ! In 2155 dkinpw,nkpg,nkpg1,kpg_k,kpg1_k,kinpw1,ffnlk,ffnl1,ph3d,ph3d1) ! Out 2156 2157 !LOOP OVER BANDS 2158 do iband=1,nband_k 2159 2160 if(mpi_enreg%proc_distrb(ikpt,iband,isppol) /= mpi_enreg%me_kpt) cycle 2161 2162 !LOOP OVER KET ATOMIC DISPERSION PERTURBATION 2163 do jatpert=1,natpert 2164 2165 !Read atomic displacement wf1 2166 call wfk_t_atdis(jatpert)%read_bks( iband, indkpt1(ikpt), & 2167 & isppol, xmpio_single, cg_bks=cg1_jatdis) 2168 2169 !Compute < g |\partial_{gamma} H^{(0)} | u_{i,k}^{\tau_{\kappa'\beta}} > 2170 call getgh1c(berryopt,cg1_jatdis,dum_cwaveprj,gv1c,dum_grad_berry,& 2171 & dum_gs1,gs_hamkq,dum_gvnl1,idir,ipert,dum_lambda,mpi_enreg,optlocal,& 2172 & optnl,opt_gvnl1,rf_hamkq,sij_opt,tim_getgh1c,usevnl) 2173 2174 !LOOP OVER BRA ATOMIC DISPLACEMENT PERTURBATION 2175 do iatpert=1,natpert 2176 2177 !Read atomic displacement wf1 2178 call wfk_t_atdis(iatpert)%read_bks( iband, indkpt1(ikpt), & 2179 & isppol, xmpio_single, cg_bks=cg1_iatdis) 2180 2181 !calculate: 2182 !< u_{i,k}^{\tau_{\kappa\beta}} | \partial_{gamma} H^{(0)} | | u_{i,k}^{\tau_{\kappa'\beta}} > 2183 call dotprod_g(dotr,doti,istwf_k,npw_k*dtset%nspinor,2,cg1_iatdis,gv1c, & 2184 & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 2185 2186 c1atdis_q1gradH0_c1atdis_bks(1,iband,iatpert,jatpert,iq1grad)= dotr 2187 c1atdis_q1gradH0_c1atdis_bks(2,iband,iatpert,jatpert,iq1grad)= doti 2188 2189 end do !iatpert 2190 2191 end do !jatpert 2192 2193 end do !iband 2194 2195 !Clean the ddk rf_hamiltonian 2196 call rf_hamkq%free() 2197 2198 !Deallocations 2199 ABI_FREE(kpg_k) 2200 ABI_FREE(kpg1_k) 2201 ABI_FREE(dkinpw) 2202 ABI_FREE(kinpw1) 2203 ABI_FREE(ffnlk) 2204 ABI_FREE(ffnl1) 2205 ABI_FREE(ph3d) 2206 2207 end do !iq1grad 2208 2209 !Deallocations 2210 ABI_FREE(cg1_jatdis) 2211 ABI_FREE(gv1c) 2212 ABI_FREE(vlocal1) 2213 2214 !--------------------------------------------------------------------------------------------------- 2215 ! < u_{i,k}^{\tau_{\kappa\alpha}} | \partial_{gamma} \hat(Q)H^{\tau_{\kappa'\beta}} | u_{i,k}^{(0)} > 2216 !--------------------------------------------------------------------------------------------------- 2217 2218 !LOOP OVER BANDS 2219 do iband=1,nband_k 2220 2221 if(mpi_enreg%proc_distrb(ikpt,iband,isppol) /= mpi_enreg%me_kpt) cycle 2222 2223 !LOOP OVER BRA ATOMIC DISPLACEMENT PERTURBATION 2224 do iatpert=1,natpert 2225 2226 !Read atomic displacement wf1 2227 call wfk_t_atdis(iatpert)%read_bks( iband, indkpt1(ikpt), & 2228 & isppol, xmpio_single, cg_bks=cg1_iatdis) 2229 2230 !LOOP OVER q1-GRADIENT 2231 do iq1grad=1,nq1grad 2232 2233 !LOOP OVER BANDS 2234 do jband=1,nband_k 2235 2236 !Read ddk wf1 2237 call wfk_t_ddk(iq1grad)%read_bks( jband, indkpt1(ikpt), & 2238 & isppol, xmpio_single, cg_bks=cg1_ddk) 2239 2240 !Calculate: < u_{i,k}^{\tau_{\kappa\alpha}} | u_{j,k}^{k_{\gamma}} > 2241 call dotprod_g(dotr,doti,istwf_k,npw_k*dtset%nspinor,2,cg1_iatdis,cg1_ddk, & 2242 & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 2243 2244 !LOOP OVER HAMILTONIAN ATOMIC DISPLACEMENT PERTURBATION 2245 do jatpert=1,natpert 2246 2247 !Calculate: -\sum_{j} < u_{i,k}^{\tau_{\kappa\alpha}} | u_{j,k}^{k_{\gamma}} > * 2248 ! < u_{j,k}^{(0)} |H^{\tau_{\kappa'\beta}}+V^{\tau_{\kappa'\beta}} | u_{i,k}^{(0)} > 2249 cprodr=dotr*ci_h1vatdis_cj(1,jatpert,jband,iband) - & 2250 & doti*ci_h1vatdis_cj(2,jatpert,jband,iband) 2251 cprodi=dotr*ci_h1vatdis_cj(2,jatpert,jband,iband) + & 2252 & doti*ci_h1vatdis_cj(1,jatpert,jband,iband) 2253 2254 c1atdis_dQHatdis_c0_bks(1,iband,iatpert,jatpert,iq1grad)= & 2255 & c1atdis_dQHatdis_c0_bks(1,iband,iatpert,jatpert,iq1grad)-cprodr 2256 c1atdis_dQHatdis_c0_bks(2,iband,iatpert,jatpert,iq1grad)= & 2257 & c1atdis_dQHatdis_c0_bks(2,iband,iatpert,jatpert,iq1grad)-cprodi 2258 2259 end do !jatpert 2260 2261 end do !jband 2262 2263 end do !iq1grad 2264 2265 end do !iatpert 2266 2267 end do !iband 2268 2269 !-------------------------------------------------------------------------------------- 2270 ! q1-gradient of atomic displacement 1st-order Hamiltonian: 2271 ! <u_{i,k}^{\tau_{\kappa\alpha}} | H^{\tau_{\kappa'\beta}}_{\gamma} | u_{i,k}^{(0)} > 2272 !-------------------------------------------------------------------------------------- 2273 2274 !Specific allocations 2275 ABI_MALLOC(vpsp1dq,(2*nfft)) 2276 ABI_MALLOC(vlocal1dq,(2*ngfft(4),ngfft(5),ngfft(6),gs_hamkq%nvloc)) 2277 ABI_MALLOC(gh1dqc,(2,npw_k*dtset%nspinor)) 2278 ABI_MALLOC(gvloc1dqc,(2,npw_k*dtset%nspinor)) 2279 ABI_MALLOC(gvnl1dqc,(2,npw_k*dtset%nspinor)) 2280 2281 !Specific definitions 2282 useylmgr1=1;optlocal=1;optnl=1 2283 nylmgreff=3+6*((ipert-dtset%natom)/10) 2284 2285 !LOOP OVER HAMILTONIAN ATOMIC DISPLACEMENT PERTURBATIONS 2286 do jatpert= 1, natpert 2287 jpert=pert_atdis(1,jatpert) 2288 jdir=pert_atdis(2,jatpert) 2289 2290 !LOOP OVER Q1-GRADIENT 2291 do iq1grad=1,nq1grad 2292 2293 !Get q-gradient of first-order local part of the pseudopotential 2294 call dfpt_vlocaldq(atindx,2,gs_hamkq%gmet,gsqcut,jdir,jpert,mpi_enreg, & 2295 & psps%mqgrid_vl,dtset%natom,& 2296 & nattyp,nfft,ngfft,dtset%ntypat,ngfft(1),ngfft(2),ngfft(3), & 2297 & ph1d,q1grad(2,iq1grad),psps%qgrid_vl,& 2298 & dtset%qptn,ucvol,psps%vlspl,vpsp1dq) 2299 2300 !Set up q-gradient of local potential vlocal1dq with proper dimensioning 2301 call rf_transgrid_and_pack(isppol,nspden,psps%usepaw,2,nfft,dtset%nfft,dtset%ngfft,& 2302 & gs_hamkq%nvloc,pawfgr,mpi_enreg,dum_vpsp,vpsp1dq,dum_vlocal,vlocal1dq) 2303 2304 !Initialize rf_hamiltonian (the k-dependent part is prepared in getgh1c_setup) 2305 call init_rf_hamiltonian(2,gs_hamkq,jpert,rf_hamkq,& 2306 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab) 2307 call rf_hamkq%load_spin(isppol,vlocal1=vlocal1dq,with_nonlocal=.true.) 2308 2309 !Set up the ground-state Hamiltonian, and some parts of the 1st-order Hamiltonian 2310 call getgh1dqc_setup(gs_hamkq,rf_hamkq,dtset,psps,kpt,kpt,jdir,jpert,q1grad(2,iq1grad), & 2311 & dtset%natom,rmet,gs_hamkq%gprimd,gs_hamkq%gmet,istwf_k,npw_k,npw_k,nylmgreff,useylmgr1,kg_k, & 2312 & ylm_k,kg_k,ylm_k,part_ylmgr_k,nkpg,nkpg1,kpg_k,kpg1_k,dkinpw,kinpw1,ffnlk,ffnl1,ph3d,ph3d1) 2313 2314 !LOOP OVER BANDS 2315 do iband=1,nband_k 2316 2317 if(mpi_enreg%proc_distrb(ikpt,iband,isppol) /= mpi_enreg%me_kpt) cycle 2318 2319 !Read ket ground-state wavefunctions 2320 cwave0i(:,:)=cg(:,1+(iband-1)*npw_k*dtset%nspinor+icg:iband*npw_k*dtset%nspinor+icg) 2321 2322 !Compute < g |H^{\tau_{\kappa'\beta}}_{\gamma} | u_{i,k}^{(0)} > 2323 call getgh1dqc(cwave0i,dum_cwaveprj,gh1dqc,gvloc1dqc,gvnl1dqc,gs_hamkq, & 2324 & jdir,jpert,mpi_enreg,optlocal,optnl,q1grad(2,iq1grad),rf_hamkq) 2325 2326 !LOOP OVER BRA ATOMIC DISPLACEMENT PERTURBATION 2327 do iatpert=1,natpert 2328 2329 !Read atomic displacement wf1 2330 call wfk_t_atdis(iatpert)%read_bks( iband, indkpt1(ikpt), & 2331 & isppol, xmpio_single, cg_bks=cg1_iatdis) 2332 2333 !Calculate: 2334 !<u_{i,k}^{\tau_{\kappa\alpha}} | H^{\tau_{\kappa'\beta}}_{\gamma} | u_{i,k}^{(0)} > 2335 call dotprod_g(dotr,doti,istwf_k,npw_k*dtset%nspinor,2,cg1_iatdis,gh1dqc, & 2336 & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 2337 2338 c1atdis_Hatdisdq_c0_bks(1,iband,iatpert,jatpert,iq1grad)= doti 2339 c1atdis_Hatdisdq_c0_bks(2,iband,iatpert,jatpert,iq1grad)= -dotr 2340 2341 end do !iatpert 2342 2343 end do !iband 2344 2345 !Clean the rf_hamiltonian 2346 call rf_hamkq%free() 2347 2348 !Deallocations 2349 ABI_FREE(kpg_k) 2350 ABI_FREE(kpg1_k) 2351 ABI_FREE(dkinpw) 2352 ABI_FREE(kinpw1) 2353 ABI_FREE(ffnlk) 2354 ABI_FREE(ffnl1) 2355 ABI_FREE(ph3d) 2356 2357 end do !iq1grad 2358 2359 end do !jatpert 2360 2361 !Deallocations 2362 ABI_FREE(dum_cwaveprj) 2363 ABI_FREE(gh1dqc) 2364 ABI_FREE(gvloc1dqc) 2365 ABI_FREE(gvnl1dqc) 2366 ABI_FREE(vpsp1dq) 2367 ABI_FREE(vlocal1dq) 2368 ABI_FREE(dum_vpsp) 2369 ABI_FREE(dum_vlocal) 2370 ABI_FREE(part_ylmgr_k) 2371 2372 !-------------------------------------------------------------------------------------- 2373 ! Acumulates all the wf dependent terms of the quadrupole tensor 2374 !-------------------------------------------------------------------------------------- 2375 ddmdqwf_k=zero 2376 ddmdqwf_t1_k=zero 2377 ddmdqwf_t2_k=zero 2378 ddmdqwf_t3_k=zero 2379 2380 do iq1grad=1,nq1grad 2381 do jatpert=1,natpert 2382 do iatpert=1,natpert 2383 do iband=1,nband_k 2384 2385 if(mpi_enreg%proc_distrb(ikpt,iband,isppol) /= mpi_enreg%me_kpt) cycle 2386 2387 !All terms toghether (Terms 4 and 5 are computed as hermitian of 2 and 3) 2388 ddmdqwf_k(1,iatpert,jatpert,iq1grad)=ddmdqwf_k(1,iatpert,jatpert,iq1grad) + & 2389 & wtk_k * occ_k(iband) * & 2390 & ( c1atdis_q1gradH0_c1atdis_bks(1,iband,iatpert,jatpert,iq1grad) + & 2391 & c1atdis_dQHatdis_c0_bks(1,iband,iatpert,jatpert,iq1grad) + & 2392 & c1atdis_dQHatdis_c0_bks(1,iband,jatpert,iatpert,iq1grad) + & 2393 & c1atdis_Hatdisdq_c0_bks(1,iband,iatpert,jatpert,iq1grad) + & 2394 & c1atdis_Hatdisdq_c0_bks(1,iband,jatpert,iatpert,iq1grad) ) 2395 2396 ddmdqwf_k(2,iatpert,jatpert,iq1grad)=ddmdqwf_k(2,iatpert,jatpert,iq1grad) + & 2397 & wtk_k * occ_k(iband) * & 2398 & ( c1atdis_q1gradH0_c1atdis_bks(2,iband,iatpert,jatpert,iq1grad) + & 2399 & c1atdis_dQHatdis_c0_bks(2,iband,iatpert,jatpert,iq1grad) - & 2400 & c1atdis_dQHatdis_c0_bks(2,iband,jatpert,iatpert,iq1grad) + & 2401 & c1atdis_Hatdisdq_c0_bks(2,iband,iatpert,jatpert,iq1grad) - & 2402 & c1atdis_Hatdisdq_c0_bks(2,iband,jatpert,iatpert,iq1grad) ) 2403 2404 !Separate them 2405 !T1 2406 ddmdqwf_t1_k(1,iatpert,jatpert,iq1grad)=ddmdqwf_t1_k(1,iatpert,jatpert,iq1grad) + & 2407 & wtk_k * occ_k(iband) * & 2408 & c1atdis_q1gradH0_c1atdis_bks(1,iband,iatpert,jatpert,iq1grad) 2409 2410 ddmdqwf_t1_k(2,iatpert,jatpert,iq1grad)=ddmdqwf_t1_k(2,iatpert,jatpert,iq1grad) + & 2411 & wtk_k * occ_k(iband) * & 2412 & c1atdis_q1gradH0_c1atdis_bks(2,iband,iatpert,jatpert,iq1grad) 2413 2414 !T2+T4 2415 ddmdqwf_t2_k(1,iatpert,jatpert,iq1grad)=ddmdqwf_t2_k(1,iatpert,jatpert,iq1grad) + & 2416 & wtk_k * occ_k(iband) * & 2417 & (c1atdis_dQHatdis_c0_bks(1,iband,iatpert,jatpert,iq1grad) + & 2418 & c1atdis_dQHatdis_c0_bks(1,iband,jatpert,iatpert,iq1grad)) 2419 2420 ddmdqwf_t2_k(2,iatpert,jatpert,iq1grad)=ddmdqwf_t2_k(2,iatpert,jatpert,iq1grad) + & 2421 & wtk_k * occ_k(iband) * & 2422 & (c1atdis_dQHatdis_c0_bks(2,iband,iatpert,jatpert,iq1grad) - & 2423 & c1atdis_dQHatdis_c0_bks(2,iband,jatpert,iatpert,iq1grad)) 2424 2425 !T3+T5 2426 ddmdqwf_t3_k(1,iatpert,jatpert,iq1grad)=ddmdqwf_t3_k(1,iatpert,jatpert,iq1grad) + & 2427 & wtk_k * occ_k(iband) * & 2428 & (c1atdis_Hatdisdq_c0_bks(1,iband,iatpert,jatpert,iq1grad) + & 2429 & c1atdis_Hatdisdq_c0_bks(1,iband,jatpert,iatpert,iq1grad)) 2430 2431 ddmdqwf_t3_k(2,iatpert,jatpert,iq1grad)=ddmdqwf_t3_k(2,iatpert,jatpert,iq1grad) + & 2432 & wtk_k * occ_k(iband) * & 2433 & (c1atdis_Hatdisdq_c0_bks(2,iband,iatpert,jatpert,iq1grad) - & 2434 & c1atdis_Hatdisdq_c0_bks(2,iband,jatpert,iatpert,iq1grad)) 2435 2436 end do 2437 end do 2438 end do 2439 end do 2440 2441 !Deallocations 2442 ABI_FREE(ci_h1vatdis_cj) 2443 ABI_FREE(c1atdis_q1gradH0_c1atdis_bks) 2444 ABI_FREE(c1atdis_dQHatdis_c0_bks) 2445 ABI_FREE(c1atdis_Hatdisdq_c0_bks) 2446 ABI_FREE(cwave0i) 2447 ABI_FREE(cg1_iatdis) 2448 ABI_FREE(cg1_ddk) 2449 2450 DBG_EXIT("COLL") 2451 2452 end subroutine dfpt_ddmdqwf
ABINIT/dfpt_isdqfr [ Functions ]
NAME
dfpt_isdqfr
FUNCTION
This routine computes the frozen wf contribution to the q-gradient of the internal strain tensor
INPUTS
atindx(natom)=index table for atoms (see gstate.f) cg(2,mpw*nspinor*mband*mkmem*nsppol)=planewave coefficients of wavefunctions at k cplex: if 1, several magnitudes are REAL, if 2, COMPLEX dtset <type(dataset_type)>=all input variables for this dataset gs_hamkq <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k gsqcut=large sphere cut-off icg=shift to be applied on the location of data in the array cg ikpt=number of the k-point isppol=1 for unpolarized, 2 for spin-polarized istwf_k=parameter that describes the storage of wfs kg_k(3,npw_k)=reduced planewave coordinates. kpt(3)=reduced coordinates of k point matom= number of atoms in the cell mkmem =number of k points treated by this node mpi_enreg=information about MPI parallelization mpw=maximum dimensioned size of npw or wfs at k natpert=number of atomic displacement perturbations nattyp(ntypat)= # atoms of each type. nband_k=number of bands at this k point for that spin polarization nfft=(effective) number of FFT grid points (for this proc) ngfft(1:18)=integer array with FFT box dimensions and other npw_k=number of plane waves at this k point nq1grad=number of q1 (q_{\gamma}) gradients nspden=number of spin-density components nsppol=1 for unpolarized, 2 for spin-polarized nstrpert=number of strain perturbations nylmgr=second dimension of ylmgr_k occ_k(nband_k)=occupation number for each band (usually 2) for each k. pert_atdis(3,natpert)=array with the info for the atomic displacement perturbations pert_strain(6,nstrpert)=array with the info for the strain perturbations ph1d(2,3*(2*dtset%mgfft+1)*dtset%natom)=1-dimensional phases psps <type(pseudopotential_type)>=variables related to pseudopotentials rmet(3,3)=real space metric (bohr**2) ucvol=unit cell volume in bohr**3. useylmgr= if 1 use the derivative of spherical harmonics wtk_k=weight assigned to the k point. ylm_k(npw_k,psps%mpsang*psps%mpsang*psps%useylm)=real spherical harmonics for the k point ylmgr_k(npw_k,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr)= k-gradients of real spherical harmonics for the k point
OUTPUT
frwfdq_k(2,matom,3,3,3,nq1grad)=frozen wave function dependent part of the 1st q-gradient of internal strain tensor
SIDE EFFECTS
NOTES
SOURCE
3479 subroutine dfpt_isdqfr(atindx,cg,dtset,frwfdq_k,gs_hamkq,gsqcut,icg,ikpt,& 3480 & isppol,istwf_k,kg_k,kpt,mkmem,mpi_enreg,matom,mpw,natpert,nattyp,nband_k,nfft,& 3481 & ngfft,npw_k,nq1grad,nspden,nsppol,nstrpert,nylmgr,occ_k,pert_atdis, & 3482 & pert_strain,ph1d,psps,rmet,ucvol,useylmgr,wtk_k,ylm_k,ylmgr_k) 3483 3484 !Arguments ------------------------------------ 3485 !scalars 3486 integer,intent(in) :: icg,ikpt,isppol,istwf_k 3487 integer,intent(in) :: matom,mkmem,mpw,natpert,nband_k,nfft 3488 integer,intent(in) :: npw_k,nq1grad,nspden,nsppol,nstrpert,nylmgr 3489 integer,intent(in) :: useylmgr 3490 real(dp),intent(in) :: gsqcut,ucvol,wtk_k 3491 type(dataset_type),intent(in) :: dtset 3492 type(gs_hamiltonian_type),intent(inout) :: gs_hamkq 3493 type(MPI_type),intent(in) :: mpi_enreg 3494 type(pseudopotential_type),intent(in) :: psps 3495 3496 !arrays 3497 integer,intent(in) :: atindx(dtset%natom) 3498 integer,intent(in) :: kg_k(3,npw_k),nattyp(dtset%ntypat),ngfft(18) 3499 integer,intent(in) :: pert_atdis(3,natpert),pert_strain(6,nstrpert) 3500 real(dp),intent(in) :: cg(2,mpw*dtset%nspinor*dtset%mband*mkmem*nsppol) 3501 real(dp),intent(out) :: frwfdq_k(2,matom,3,3,3,nq1grad) 3502 real(dp),intent(in) :: kpt(3),occ_k(nband_k) 3503 real(dp),intent(in) :: ph1d(2,3*(2*dtset%mgfft+1)*dtset%natom) 3504 real(dp),intent(in) :: rmet(3,3) 3505 real(dp),intent(in) :: ylm_k(npw_k,psps%mpsang*psps%mpsang*psps%useylm) 3506 real(dp),intent(in) :: ylmgr_k(npw_k,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr) 3507 3508 !Local variables------------------------------- 3509 !scalars 3510 integer :: iatdir,iatom,iatpert,iband,idir,ipert,ipw,iq1grad,iq2grad,istrpert,ka,kb 3511 integer :: nkpg,nkpg1,nylmgrpart,optlocal,optnl,tim_getgh1c,useylmgr1 3512 real(dp) :: doti,dotr 3513 type(pawfgr_type) :: pawfgr 3514 type(rf_hamiltonian_type) :: rf_hamkq 3515 3516 !arrays 3517 real(dp),allocatable :: c0_hatdisdq_c0_bks(:,:,:,:) 3518 real(dp),allocatable :: cwave0i(:,:) 3519 real(dp),allocatable :: dkinpw(:) 3520 real(dp),allocatable :: ffnlk(:,:,:,:),ffnl1(:,:,:,:) 3521 real(dp),allocatable :: frwfdq_bks(:,:,:,:,:,:,:) 3522 real(dp),allocatable :: gh1dqc(:,:),gh1dqpkc(:,:),gvloc1dqc(:,:),gvnl1dqc(:,:) 3523 real(dp),allocatable :: c0_ghatdisdqdq_pk_c0(:,:,:,:,:,:) 3524 real(dp),allocatable :: kinpw1(:),kpg_k(:,:),kpg1_k(:,:),kpg_pk(:,:),ph3d(:,:,:),ph3d1(:,:,:) 3525 real(dp),allocatable :: dum_vlocal(:,:,:,:),vlocal1dq(:,:,:,:), dum_vpsp(:) 3526 real(dp),allocatable :: vpsp1dq(:),part_ylmgr_k(:,:,:) 3527 type(pawcprj_type),allocatable :: dum_cwaveprj(:,:) 3528 3529 3530 DBG_ENTER("COLL") 3531 3532 !Additional definitions 3533 tim_getgh1c=0 3534 3535 !Additional allocations 3536 ABI_MALLOC(cwave0i,(2,npw_k*dtset%nspinor)) 3537 ABI_MALLOC(dum_vpsp,(nfft)) 3538 ABI_MALLOC(dum_vlocal,(ngfft(4),ngfft(5),ngfft(6),gs_hamkq%nvloc)) 3539 ABI_MALLOC(dum_cwaveprj,(0,0)) 3540 ABI_MALLOC(vpsp1dq,(2*nfft)) 3541 ABI_MALLOC(vlocal1dq,(2*ngfft(4),ngfft(5),ngfft(6),gs_hamkq%nvloc)) 3542 ABI_MALLOC(gh1dqc,(2,npw_k*dtset%nspinor)) 3543 ABI_MALLOC(gvloc1dqc,(2,npw_k*dtset%nspinor)) 3544 ABI_MALLOC(gvnl1dqc,(2,npw_k*dtset%nspinor)) 3545 3546 !Additional definitions 3547 useylmgr1=1;optlocal=1;optnl=1 3548 3549 !----------------------------------------------------------------------------------------------- 3550 ! q1-gradient of atomic displacement 1st order hamiltonian: 3551 ! < u_{i,k}^{(0)} | H^{\tau_{\kappa\alpha}_{\gamma} | u_{i,k}^{(0)} > 3552 !----------------------------------------------------------------------------------------------- 3553 3554 !Specific allocations 3555 ABI_MALLOC(c0_hatdisdq_c0_bks,(2,nband_k,3,natpert)) 3556 ABI_MALLOC(part_ylmgr_k,(npw_k,3, psps%mpsang*psps%mpsang*psps%useylm*useylmgr1)) 3557 part_ylmgr_k(:,:,:)=ylmgr_k(:,1:3,:) 3558 3559 !Specific definitions 3560 nylmgrpart=3 3561 3562 !LOOP OVER HAMILTONIAN ATOMIC DISPLACEMENT PERTURBATIONS 3563 do iatpert=1,natpert 3564 ipert=pert_atdis(1,iatpert) 3565 idir=pert_atdis(2,iatpert) 3566 3567 !LOOP OVER Q1-GRADIENT 3568 do iq1grad=1,3 3569 3570 !Get q-gradient of first-order local part of the pseudopotential 3571 call dfpt_vlocaldq(atindx,2,gs_hamkq%gmet,gsqcut,idir,ipert,mpi_enreg, & 3572 & psps%mqgrid_vl,dtset%natom,& 3573 & nattyp,nfft,ngfft,dtset%ntypat,ngfft(1),ngfft(2),ngfft(3), & 3574 & ph1d,iq1grad,psps%qgrid_vl,& 3575 & dtset%qptn,ucvol,psps%vlspl,vpsp1dq) 3576 3577 !Set up q-gradient of local potential vlocal1dq with proper dimensioning 3578 call rf_transgrid_and_pack(isppol,nspden,psps%usepaw,2,nfft,dtset%nfft,dtset%ngfft,& 3579 & gs_hamkq%nvloc,pawfgr,mpi_enreg,dum_vpsp,vpsp1dq,dum_vlocal,vlocal1dq) 3580 3581 !Initialize rf_hamiltonian (the k-dependent part is prepared in getgh1c_setup) 3582 call init_rf_hamiltonian(2,gs_hamkq,ipert,rf_hamkq,& 3583 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab) 3584 call rf_hamkq%load_spin(isppol,vlocal1=vlocal1dq,with_nonlocal=.true.) 3585 3586 !Set up the ground-state Hamiltonian, and some parts of the 1st-order Hamiltonian 3587 call getgh1dqc_setup(gs_hamkq,rf_hamkq,dtset,psps,kpt,kpt,idir,ipert,iq1grad, & 3588 & dtset%natom,rmet,gs_hamkq%gprimd,gs_hamkq%gmet,istwf_k,npw_k,npw_k,nylmgrpart,useylmgr1,kg_k, & 3589 & ylm_k,kg_k,ylm_k,part_ylmgr_k,nkpg,nkpg1,kpg_k,kpg1_k,dkinpw,kinpw1,ffnlk,ffnl1,ph3d,ph3d1) 3590 3591 !LOOP OVER BANDS 3592 do iband=1,nband_k 3593 3594 if(mpi_enreg%proc_distrb(ikpt,iband,isppol) /= mpi_enreg%me_kpt) cycle 3595 3596 !Read ket ground-state wavefunctions 3597 cwave0i(:,:)=cg(:,1+(iband-1)*npw_k*dtset%nspinor+icg:iband*npw_k*dtset%nspinor+icg) 3598 3599 !Compute < g |H^{\tau_{\kappa\alpha}}_{\gamma} | u_{i,k}^{(0)} > 3600 call getgh1dqc(cwave0i,dum_cwaveprj,gh1dqc,gvloc1dqc,gvnl1dqc,gs_hamkq, & 3601 & idir,ipert,mpi_enreg,optlocal,optnl,iq1grad,rf_hamkq) 3602 3603 !Calculate: 3604 !<u_{i,k}^{(0)} | H^{\tau_{\kappa\alpha}}_{\gamma} | u_{i,k}^{(0)} > 3605 call dotprod_g(dotr,doti,istwf_k,npw_k*dtset%nspinor,2,cwave0i,gh1dqc, & 3606 & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 3607 3608 c0_hatdisdq_c0_bks(1,iband,iq1grad,iatpert)=dotr 3609 c0_hatdisdq_c0_bks(2,iband,iq1grad,iatpert)=doti 3610 3611 end do !iband 3612 3613 !Clean the rf_hamiltonian 3614 call rf_hamkq%free() 3615 3616 !Deallocations 3617 ABI_FREE(kpg_k) 3618 ABI_FREE(kpg1_k) 3619 ABI_FREE(dkinpw) 3620 ABI_FREE(kinpw1) 3621 ABI_FREE(ffnlk) 3622 ABI_FREE(ffnl1) 3623 ABI_FREE(ph3d) 3624 3625 end do !iq1grad 3626 3627 end do !iatpert 3628 3629 ABI_FREE(part_ylmgr_k) 3630 3631 !----------------------------------------------------------------------------------------------- 3632 ! 2nd q-gradient of atomic displacement 1st order hamiltonian * momentum operator : 3633 ! <u_{i,k}^{(0)} | H^{\tau_{\kappa\alpha}}_{\gamma\delta} (k+G)_{\beta} | u_{i,k}^{(0)} > 3634 !----------------------------------------------------------------------------------------------- 3635 3636 !Specific allocations 3637 ABI_MALLOC(gh1dqpkc,(2,npw_k*dtset%nspinor)) 3638 ABI_MALLOC(c0_ghatdisdqdq_pk_c0,(2,nband_k,3,3,3,natpert)) 3639 3640 !Generate k+G vectors 3641 nkpg=3 3642 ABI_MALLOC(kpg_pk,(npw_k,nkpg)) 3643 call mkkpg(kg_k,kpg_pk,kpt,nkpg,npw_k) 3644 3645 !LOOP OVER ATOMIC DISPLACEMENT PERTURBATIONS 3646 do iatpert=1,natpert 3647 ipert=pert_atdis(1,iatpert) 3648 idir=pert_atdis(2,iatpert) 3649 3650 !LOOP OVER Q1-GRADIENT 3651 do iq1grad=1,3 3652 3653 !LOOP OVER Q2-GRADIENT 3654 do iq2grad=1,3 3655 3656 !Get q-gradient of first-order local part of the pseudopotential 3657 call dfpt_vlocaldqdq(atindx,2,gs_hamkq%gmet,gsqcut,idir,ipert,mpi_enreg, & 3658 & psps%mqgrid_vl,dtset%natom,& 3659 & nattyp,nfft,ngfft,dtset%ntypat,ngfft(1),ngfft(2),ngfft(3), & 3660 & ph1d,iq1grad,iq2grad,psps%qgrid_vl,& 3661 & dtset%qptn,ucvol,psps%vlspl,vpsp1dq) 3662 3663 !Set up q-gradient of local potential vlocal1dq with proper dimensioning 3664 call rf_transgrid_and_pack(isppol,nspden,psps%usepaw,2,nfft,dtset%nfft,dtset%ngfft,& 3665 & gs_hamkq%nvloc,pawfgr,mpi_enreg,dum_vpsp,vpsp1dq,dum_vlocal,vlocal1dq) 3666 3667 !Initialize rf_hamiltonian (the k-dependent part is prepared in getgh1c_setup) 3668 call init_rf_hamiltonian(2,gs_hamkq,ipert,rf_hamkq,& 3669 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab) 3670 call rf_hamkq%load_spin(isppol,vlocal1=vlocal1dq,with_nonlocal=.true.) 3671 3672 !Set up the ground-state Hamiltonian, and some parts of the 1st-order Hamiltonian 3673 call getgh1dqc_setup(gs_hamkq,rf_hamkq,dtset,psps,kpt,kpt,idir,ipert,iq1grad, & 3674 & dtset%natom,rmet,gs_hamkq%gprimd,gs_hamkq%gmet,istwf_k,npw_k,npw_k,nylmgr,useylmgr1,kg_k, & 3675 & ylm_k,kg_k,ylm_k,ylmgr_k,nkpg,nkpg1,kpg_k,kpg1_k,dkinpw,kinpw1,ffnlk,ffnl1,ph3d,ph3d1, & 3676 & qdir2=iq2grad) 3677 3678 !LOOP OVER BANDS 3679 do iband=1,nband_k 3680 3681 if(mpi_enreg%proc_distrb(ikpt,iband,isppol) /= mpi_enreg%me_kpt) cycle 3682 3683 !Read ket ground-state wavefunctions 3684 cwave0i(:,:)=cg(:,1+(iband-1)*npw_k*dtset%nspinor+icg:iband*npw_k*dtset%nspinor+icg) 3685 3686 !Compute < g |H^{\tau_{\kappa\alpha}}_{\gamma\delta} | u_{i,k}^{(0)} > 3687 call getgh1dqc(cwave0i,dum_cwaveprj,gh1dqc,gvloc1dqc,gvnl1dqc,gs_hamkq, & 3688 & idir,ipert,mpi_enreg,optlocal,optnl,iq1grad,rf_hamkq,qdir2=iq2grad) 3689 3690 !LOOP OVER ONE OF THE STRAIN DIRECTIONS 3691 do ka=1,3 3692 do ipw=1,npw_k 3693 gh1dqpkc(:,ipw)=gh1dqc(:,ipw)*kpg_pk(ipw,ka) 3694 end do 3695 3696 call dotprod_g(dotr,doti,istwf_k,npw_k*dtset%nspinor,2,cwave0i,gh1dqpkc, & 3697 & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 3698 3699 c0_ghatdisdqdq_pk_c0(1,iband,ka,iq1grad,iq2grad,iatpert)=dotr 3700 c0_ghatdisdqdq_pk_c0(2,iband,ka,iq1grad,iq2grad,iatpert)=doti 3701 end do 3702 3703 end do !iband 3704 3705 !Clean the rf_hamiltonian 3706 call rf_hamkq%free() 3707 3708 !Deallocations 3709 ABI_FREE(kpg_k) 3710 ABI_FREE(kpg1_k) 3711 ABI_FREE(dkinpw) 3712 ABI_FREE(kinpw1) 3713 ABI_FREE(ffnlk) 3714 ABI_FREE(ffnl1) 3715 ABI_FREE(ph3d) 3716 3717 end do !iq2grad 3718 3719 end do !iq1grad 3720 3721 end do !iatpert 3722 3723 !Deallocations 3724 ABI_FREE(dum_cwaveprj) 3725 ABI_FREE(gh1dqc) 3726 ABI_FREE(gh1dqpkc) 3727 ABI_FREE(gvloc1dqc) 3728 ABI_FREE(gvnl1dqc) 3729 ABI_FREE(vpsp1dq) 3730 ABI_FREE(vlocal1dq) 3731 ABI_FREE(dum_vpsp) 3732 ABI_FREE(dum_vlocal) 3733 ABI_FREE(kpg_pk) 3734 ABI_FREE(cwave0i) 3735 3736 !-------------------------------------------------------------------------------------- 3737 ! Acumulates the three frozen wf terms of the q-gradient of the internal strain 3738 !-------------------------------------------------------------------------------------- 3739 3740 !Specific allocations and definitions 3741 ABI_MALLOC(frwfdq_bks,(2,nband_k,matom,3,3,3,nq1grad)) 3742 ! fac=pi*pi/ucvol 3743 3744 3745 !LOOP OVER ATOMIC DISPLACEMENT PERTURBATIONS 3746 do iatpert= 1, natpert 3747 iatom=pert_atdis(1,iatpert) 3748 iatdir=pert_atdis(2,iatpert) 3749 3750 !LOOP OVER STRAIN PERTURBATIONS 3751 do istrpert= 1, nstrpert 3752 ka=pert_strain(3,istrpert) 3753 kb=pert_strain(4,istrpert) 3754 3755 !LOOP OVER Q1-GRADIENT 3756 do iq1grad=1,3 3757 3758 !LOOP OVER BANDS 3759 do iband=1,nband_k 3760 3761 if(mpi_enreg%proc_distrb(ikpt,iband,isppol) /= mpi_enreg%me_kpt) cycle 3762 3763 !First accumulate the term involving the 2nd q-gradient of atdis Hamiltonian: 3764 !Take here into account the -i*(-i)^{\dagger} prefactors. 3765 frwfdq_bks(1,iband,iatom,iatdir,ka,kb,iq1grad)=c0_ghatdisdqdq_pk_c0(1,iband,ka,iq1grad,kb,iatpert) 3766 frwfdq_bks(2,iband,iatom,iatdir,ka,kb,iq1grad)=-c0_ghatdisdqdq_pk_c0(2,iband,ka,iq1grad,kb,iatpert) 3767 3768 !Next complete the other two terms involving the 1st q-gradient of atdis Hamiltonian: 3769 !<u_{i,k}^{(0)} | H^{\tau_{\kappa\alpha}}_{\gamma} \frac{\delta_{\beta\delta}}{2} | u_{i,k}^{(0)} > 3770 !<u_{i,k}^{(0)} | H^{\tau_{\kappa\alpha}}_{\delta} \frac{\delta_{\beta\gamma}}{2} | u_{i,k}^{(0)} > 3771 !-------------------------------------------------------------------------- 3772 if (ka==kb) then 3773 frwfdq_bks(1,iband,iatom,iatdir,ka,kb,iq1grad)=frwfdq_bks(1,iband,iatom,iatdir,ka,kb,iq1grad)+ & 3774 & half*c0_hatdisdq_c0_bks(1,iband,iq1grad,iatpert) 3775 frwfdq_bks(2,iband,iatom,iatdir,ka,kb,iq1grad)=frwfdq_bks(2,iband,iatom,iatdir,ka,kb,iq1grad)- & 3776 & half*c0_hatdisdq_c0_bks(2,iband,iq1grad,iatpert) 3777 end if 3778 if (ka==iq1grad) then 3779 frwfdq_bks(1,iband,iatom,iatdir,ka,kb,iq1grad)=frwfdq_bks(1,iband,iatom,iatdir,ka,kb,iq1grad)+ & 3780 & half*c0_hatdisdq_c0_bks(1,iband,kb,iatpert) 3781 frwfdq_bks(2,iband,iatom,iatdir,ka,kb,iq1grad)=frwfdq_bks(2,iband,iatom,iatdir,ka,kb,iq1grad)- & 3782 & half*c0_hatdisdq_c0_bks(2,iband,kb,iatpert) 3783 end if 3784 3785 !Take into account the two pi factor from the term 3786 !(\hat{p}_{k\beta + \frac{q_{\beta}}{2}}) appearing before the double q-derivation 3787 frwfdq_bks(:,iband,iatom,iatdir,ka,kb,iq1grad)=frwfdq_bks(:,iband,iatom,iatdir,ka,kb,iq1grad)* & 3788 & two_pi 3789 3790 end do !iband 3791 3792 end do !iq1grad 3793 3794 end do !istrpert 3795 3796 end do !iatpert 3797 3798 !-------------------------------------------------------------------------------------- 3799 ! Acumulate the frwf terms of this ikpt 3800 !-------------------------------------------------------------------------------------- 3801 frwfdq_k=zero 3802 do iatpert= 1, natpert 3803 iatom=pert_atdis(1,iatpert) 3804 iatdir=pert_atdis(2,iatpert) 3805 do istrpert= 1, nstrpert 3806 ka=pert_strain(3,istrpert) 3807 kb=pert_strain(4,istrpert) 3808 do iq1grad=1,3 3809 do iband=1,nband_k 3810 3811 if(mpi_enreg%proc_distrb(ikpt,iband,isppol) /= mpi_enreg%me_kpt) cycle 3812 3813 frwfdq_k(:,iatom,iatdir,ka,kb,iq1grad)=frwfdq_k(:,iatom,iatdir,ka,kb,iq1grad) + & 3814 & occ_k(iband) * frwfdq_bks(:,iband,iatom,iatdir,ka,kb,iq1grad) 3815 3816 end do 3817 end do 3818 end do 3819 end do 3820 3821 !scale by the k-point weight 3822 frwfdq_k=frwfdq_k * wtk_k 3823 3824 !Deallocations 3825 ABI_FREE(c0_hatdisdq_c0_bks) 3826 ABI_FREE(c0_ghatdisdqdq_pk_c0) 3827 ABI_FREE(frwfdq_bks) 3828 3829 3830 DBG_EXIT("COLL") 3831 3832 end subroutine dfpt_isdqfr
ABINIT/dfpt_isdqwf [ Functions ]
NAME
dfpt_isdqwf
FUNCTION
This routine computes the band and kpt resolved contributions to the flexoelectric tensor.
INPUTS
atindx(natom)=index table for atoms (see gstate.f) cg(2,mpw*nspinor*mband*mkmem*nsppol)=planewave coefficients of wavefunctions at k cplex: if 1, several magnitudes are REAL, if 2, COMPLEX dtset <type(dataset_type)>=all input variables for this dataset gs_hamkq <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k gsqcut=large sphere cut-off icg=shift to be applied on the location of data in the array cg ikpt=number of the k-point indkpt1(nkpt_rbz)=non-symmetrized indices of the k-points isppol=1 for unpolarized, 2 for spin-polarized istwf_k=parameter that describes the storage of wfs kg_k(3,npw_k)=reduced planewave coordinates. kpt(3)=reduced coordinates of k point kxc(nfft,nkxc)=exchange and correlation kernel matom= number of atoms in the unit cell mkmem =number of k points treated by this node mpi_enreg=information about MPI parallelization mpw=maximum dimensioned size of npw or wfs at k natpert=number of atomic displacement perturbations nattyp(ntypat)= # atoms of each type. nband_k=number of bands at this k point for that spin polarization nfft=(effective) number of FFT grid points (for this proc) ngfft(1:18)=integer array with FFT box dimensions and other nkpt_rbz= number of k-points in the RBZ nkxc=second dimension of the kxc array. If /=0, the XC kernel must be computed. npw_k=number of plane waves at this k point nq1grad=number of q1 (q_{\gamma}) gradients nspden=number of spin-density components nsppol=1 for unpolarized, 2 for spin-polarized nstrpert=number of strain perturbations nylmgr=second dimension of ylmgr_k occ_k(nband_k)=occupation number for each band (usually 2) for each k. pert_atdis(3,natpert)=array with the info for the atomic displacement perturbations pert_strain(6,nstrpert)=array with the info for the strain perturbations ph1d(2,3*(2*dtset%mgfft+1)*dtset%natom)=1-dimensional phases psps <type(pseudopotential_type)>=variables related to pseudopotentials q1grad(3,nq1grad)=array with the info for the q1 (q_{\gamma}) gradients rhog(2,nfftf)=array for Fourier transform of GS electron density rhor(nfftf,nspden)=array for GS electron density in electrons/bohr**3. rmet(3,3)=real space metric (bohr**2) ucvol=unit cell volume in bohr**3. useylmgr= if 1 use the derivative of spherical harmonics vhxc1_atdis(natpert,cplex*nfft)= electrostatic potential generated by a first order atomic displacement density vhxc1_strain(nstrpert,cplex*nfft)= electrostatic potential generated by a first order strain density wfk_t_ddk(nq1grad)=unit numbers for the ddk wf1 files wfk_t_atdis(natpert)=unit numbers for the atomic displacement wf1 files wfk_t_strain(3,3)=unit numbers for the strain wf1 files wtk_k=weight assigned to the k point. xred(3,natom)=reduced dimensionless atomic coordinates ylm_k(npw_k,psps%mpsang*psps%mpsang*psps%useylm)=real spherical harmonics for the k point ylmgr_k(npw_k,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr)= k-gradients of real spherical harmonics for the k point
OUTPUT
isdqwf_k(2,matom,3,3,3,nq1grad)=wave function dependent part of the electronic flexoelectric tensor for the k-point kpt isdqwf_t1_k(2,matom,3,nq1grad,3,3)=t1 term (see notes) of isdqwf_k isdqwf_t2_k(2,matom,3,nq1grad,3,3)=t2 term (see notes) of isdqwf_k isdqwf_t3_k(2,matom,3,nq1grad,3,3)=t3 term (see notes) of isdqwf_k isdqwf_t4_k(2,matom,3,3,3,nq1grad)=t4 term (see notes) of isdqwf_k isdqwf_t5_k(2,matom,3,nq1grad,3,3)=t5 term (see notes) of isdqwf_k
SIDE EFFECTS
NOTES
SOURCE
2536 subroutine dfpt_isdqwf(atindx,cg,cplex,dtset,gs_hamkq,gsqcut,icg,ikpt,indkpt1,isdqwf_k, & 2537 & isdqwf_t1_k,isdqwf_t2_k,isdqwf_t3_k,isdqwf_t4_k,isdqwf_t5_k,isppol,istwf_k, & 2538 & kg_k,kpt,kxc,matom,mkmem,mpi_enreg,mpw,natpert,nattyp,nband_k,nfft,ngfft,nkpt_rbz,nkxc, & 2539 & npw_k,nq1grad,nspden,nsppol,nstrpert,nylmgr,occ_k, & 2540 & pert_atdis,pert_strain,ph1d,psps,q1grad,rhog,rhor,rmet,ucvol,useylmgr, & 2541 & vhxc1_atdis,vhxc1_strain,wfk_t_atdis,wfk_t_ddk, & 2542 & wfk_t_strain,wtk_k,xred,ylm_k,ylmgr_k) 2543 2544 !Arguments ------------------------------------ 2545 !scalars 2546 integer,intent(in) :: cplex,icg,ikpt,isppol,istwf_k 2547 integer,intent(in) :: matom,mkmem,mpw,natpert,nband_k,nfft 2548 integer,intent(in) :: nkpt_rbz,nkxc,npw_k,nq1grad,nspden,nsppol,nstrpert,nylmgr 2549 integer,intent(in) :: useylmgr 2550 real(dp),intent(in) :: gsqcut,ucvol,wtk_k 2551 type(dataset_type),intent(in) :: dtset 2552 type(gs_hamiltonian_type),intent(inout) :: gs_hamkq 2553 type(MPI_type),intent(in) :: mpi_enreg 2554 type(pseudopotential_type),intent(in) :: psps 2555 2556 !arrays 2557 integer,intent(in) :: atindx(dtset%natom) 2558 integer,intent(in) :: indkpt1(nkpt_rbz),kg_k(3,npw_k),nattyp(dtset%ntypat),ngfft(18) 2559 integer,intent(in) :: pert_atdis(3,natpert),pert_strain(6,nstrpert) 2560 integer,intent(in) :: q1grad(3,nq1grad) 2561 real(dp),intent(in) :: cg(2,mpw*dtset%nspinor*dtset%mband*mkmem*nsppol) 2562 real(dp),intent(out) :: isdqwf_k(2,matom,3,nq1grad,3,3) 2563 real(dp),intent(out) :: isdqwf_t1_k(2,matom,3,nq1grad,3,3),isdqwf_t2_k(2,matom,3,nq1grad,3,3) 2564 real(dp),intent(out) :: isdqwf_t3_k(2,matom,3,nq1grad,3,3),isdqwf_t4_k(2,matom,3,3,3,nq1grad) 2565 real(dp),intent(out) :: isdqwf_t5_k(2,matom,3,nq1grad,3,3) 2566 real(dp),intent(in) :: kpt(3),occ_k(nband_k),kxc(nfft,nkxc) 2567 real(dp),intent(in) :: ph1d(2,3*(2*dtset%mgfft+1)*dtset%natom) 2568 real(dp),intent(in) :: rhog(2,nfft),rhor(nfft,nspden),rmet(3,3) 2569 real(dp),intent(in) :: vhxc1_atdis(natpert,cplex*nfft) 2570 real(dp),intent(in) :: vhxc1_strain(nstrpert,cplex*nfft) 2571 real(dp),intent(in) :: xred(3,dtset%natom) 2572 real(dp),intent(in) :: ylm_k(npw_k,psps%mpsang*psps%mpsang*psps%useylm) 2573 real(dp),intent(in) :: ylmgr_k(npw_k,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr) 2574 type(wfk_t),intent(inout) :: wfk_t_ddk(nq1grad) 2575 type(wfk_t),intent(inout) :: wfk_t_atdis(natpert),wfk_t_strain(3,3) 2576 2577 !Local variables------------------------------- 2578 !scalars 2579 integer :: iatdir,iatom,iatpert,iq1grad,berryopt,g0term,iband,idir,ii,ipert,istr,istrpert,jband 2580 integer :: ka,kb,nkpg,nkpg1,npw_disk,nylmgrpart 2581 integer :: opt_gvnl1,opthartdqdq,optlocal,optnl,sij_opt,tim_getgh1c,usevnl,useylmgr1 2582 real(dp) :: cprodr,cprodi,doti,dotr,dum_lambda 2583 character(len=500) :: msg 2584 type(rf_hamiltonian_type) :: rf_hamkq 2585 type(pawfgr_type) :: pawfgr 2586 2587 !arrays 2588 real(dp) :: dum_grad_berry(1,1),dum_gs1(1,1),dum_gvnl1(1,1) 2589 real(dp),allocatable :: cg1_atdis(:,:) 2590 real(dp),allocatable :: c1atdis_dQcalHstrain_c0_bks(:,:,:,:,:) 2591 real(dp),allocatable :: c1atdis_Hmetricdqdq_c0_bks(:,:,:,:,:) 2592 real(dp),allocatable :: c1atdis_q1gradH0_c1strain_bks(:,:,:,:,:) 2593 real(dp),allocatable :: cg1_ddk(:,:) 2594 real(dp),allocatable :: cg1_strain(:,:),cg1_strain_ar(:,:,:,:) 2595 real(dp),allocatable :: c0_HatdisdagdQ_c1strain_bks(:,:,:,:,:) 2596 real(dp),allocatable :: c0_Hatdisdqdag_c1strain_bks(:,:,:,:,:) 2597 real(dp),allocatable :: ci_h1vatdisdag_cj(:,:,:,:),cj_h1vstrain_ci(:,:,:,:) 2598 real(dp),allocatable :: cwave0i(:,:),cwave0j(:,:) 2599 real(dp),allocatable :: dkinpw(:) 2600 real(dp),allocatable :: ffnlk(:,:,:,:),ffnl1(:,:,:,:) 2601 real(dp),allocatable :: gh1dqc(:,:),gh1dqdqc(:,:) 2602 real(dp),allocatable :: gvloc1dqc(:,:),gvloc1dqdqc(:,:) 2603 real(dp),allocatable :: gvnl1dqc(:,:),gvnl1dqdqc(:,:),gv1c(:,:) 2604 real(dp),allocatable :: kinpw1(:),kpg_k(:,:),kpg1_k(:,:),ph3d(:,:,:),ph3d1(:,:,:) 2605 real(dp),allocatable :: vhart1dqdq(:) 2606 real(dp),allocatable :: dum_vlocal(:,:,:,:),vlocal1(:,:,:,:) 2607 real(dp),allocatable :: vlocal1dq(:,:,:,:),vlocal1dqdq(:,:,:,:),dum_vpsp(:) 2608 real(dp),allocatable :: vpsp1(:),vpsp1dq(:),vpsp1dqdq(:),vxc1dqdq(:) 2609 real(dp),allocatable :: dum_ylmgr1_k(:,:,:),part_ylmgr_k(:,:,:) 2610 type(pawcprj_type),allocatable :: dum_cwaveprj(:,:) 2611 2612 ! ************************************************************************* 2613 2614 DBG_ENTER("COLL") 2615 2616 if(dtset%prtvol>2)then 2617 write(msg,'(2a,i5,2x,a,3f9.5)')ch10,' First q-moment internal strain calculation; k pt #',ikpt,'k=',kpt(:) 2618 call wrtout(std_out,msg) 2619 end if 2620 2621 !Additional definitions 2622 tim_getgh1c=0 2623 2624 !Additional allocations 2625 ABI_MALLOC(cwave0i,(2,npw_k*dtset%nspinor)) 2626 ABI_MALLOC(cwave0j,(2,npw_k*dtset%nspinor)) 2627 ABI_MALLOC(gv1c,(2,npw_k*dtset%nspinor)) 2628 ABI_MALLOC(vlocal1,(cplex*ngfft(4),ngfft(5),ngfft(6),gs_hamkq%nvloc)) 2629 ABI_MALLOC(dum_vpsp,(nfft)) 2630 ABI_MALLOC(dum_vlocal,(ngfft(4),ngfft(5),ngfft(6),gs_hamkq%nvloc)) 2631 ABI_MALLOC(vpsp1,(cplex*nfft)) 2632 ABI_MALLOC(dum_cwaveprj,(0,0)) 2633 2634 !-------------------------------------------------------------------------------------- 2635 !Calculate first terms involving only ground state wavefunctions 2636 !-------------------------------------------------------------------------------------- 2637 2638 !----------------------------------------------------------------------------------------------- 2639 ! Strain 1st order hamiltonian + 1st order potential matrix element: 2640 ! < u_{j,k}^{(0)} | H^{n_{\beta\delta}}+V^{n_{\beta\delta}} | u_{i,k}^{(0)} > 2641 !----------------------------------------------------------------------------------------------- 2642 2643 !Specific definitions 2644 ipert=dtset%natom+3 2645 useylmgr1=1 2646 dum_lambda=zero 2647 berryopt=0;optlocal=1;optnl=1;usevnl=0;opt_gvnl1=0;sij_opt=0 2648 g0term=1 2649 2650 !Specific allocations 2651 ABI_MALLOC(cj_h1vstrain_ci,(2,nstrpert,nband_k,nband_k)) 2652 ABI_MALLOC(part_ylmgr_k,(npw_k,3+6*((ipert-dtset%natom)/10), psps%mpsang*psps%mpsang*psps%useylm*useylmgr1)) 2653 part_ylmgr_k(:,:,:)=ylmgr_k(:,1:3+6*((ipert-dtset%natom)/10),:) 2654 2655 !LOOP OVER STRAIN PERTURBATIONS 2656 do istrpert=1,nstrpert 2657 ipert=pert_strain(1,istrpert) 2658 idir=pert_strain(2,istrpert) 2659 istr=idir 2660 if(ipert==dtset%natom+4) istr=idir+3 2661 2662 !Get first-order local part of the pseudopotential 2663 call vlocalstr(gs_hamkq%gmet,gs_hamkq%gprimd,gsqcut,istr,dtset%mgfft,mpi_enreg,& 2664 & psps%mqgrid_vl,dtset%natom,nattyp,nfft,ngfft,dtset%ntypat,ph1d,psps%qgrid_vl,& 2665 & ucvol,psps%vlspl,vpsp1,g0term=g0term) 2666 2667 !Set up local potential vlocal1 with proper dimensioning, from vpsp1 + vhxc1_strain 2668 vpsp1=vpsp1+vhxc1_strain(istrpert,:) 2669 call rf_transgrid_and_pack(isppol,nspden,psps%usepaw,cplex,nfft,dtset%nfft,dtset%ngfft,& 2670 & gs_hamkq%nvloc,pawfgr,mpi_enreg,dum_vpsp,vpsp1,dum_vlocal,vlocal1) 2671 2672 !Initialize rf_hamiltonian (the k-dependent part is prepared in getgh1c_setup) 2673 call init_rf_hamiltonian(cplex,gs_hamkq,ipert,rf_hamkq,& 2674 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab) 2675 call rf_hamkq%load_spin(isppol,vlocal1=vlocal1,with_nonlocal=.true.) 2676 2677 !Set up the ground-state Hamiltonian, and some parts of the 1st-order Hamiltonian 2678 call getgh1c_setup(gs_hamkq,rf_hamkq,dtset,psps,& ! In 2679 kpt,kpt,idir,ipert,dtset%natom,rmet,gs_hamkq%gprimd,gs_hamkq%gmet,istwf_k,& ! In 2680 npw_k,npw_k,useylmgr1,kg_k,ylm_k,kg_k,ylm_k,part_ylmgr_k,& ! In 2681 dkinpw,nkpg,nkpg1,kpg_k,kpg1_k,kinpw1,ffnlk,ffnl1,ph3d,ph3d1) ! Out 2682 2683 !LOOP OVER KET BANDS 2684 do iband=1,nband_k 2685 2686 if(mpi_enreg%proc_distrb(ikpt,iband,isppol) /= mpi_enreg%me_kpt) cycle 2687 2688 !Read ket ground-state wavefunctions 2689 cwave0i(:,:)=cg(:,1+(iband-1)*npw_k*dtset%nspinor+icg:iband*npw_k*dtset%nspinor+icg) 2690 2691 !Compute < g | H^{n_{\beta\delta}}+V^{n_{\beta\delta}} | u_{i,k}^{(0)} > 2692 call getgh1c(berryopt,cwave0i,dum_cwaveprj,gv1c,dum_grad_berry,& 2693 & dum_gs1,gs_hamkq,dum_gvnl1,idir,ipert,dum_lambda,mpi_enreg,optlocal,& 2694 & optnl,opt_gvnl1,rf_hamkq,sij_opt,tim_getgh1c,usevnl) 2695 2696 !LOOP OVER BRA BANDS 2697 do jband=1,nband_k 2698 2699 !Read bra ground-state wavefunctions 2700 cwave0j(:,:)=cg(:,1+(jband-1)*npw_k*dtset%nspinor+icg:jband*npw_k*dtset%nspinor+icg) 2701 2702 !Compute = cj_h1vstrain_ci 2703 !< u_{j,k}^{(0)} | H^{n_{\beta\delta}}+V^{n_{\beta\delta}} | u_{i,k}^{(0)} > 2704 call dotprod_g(dotr,doti,istwf_k,npw_k*dtset%nspinor,2,cwave0j,gv1c,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 2705 cj_h1vstrain_ci(1,istrpert,jband,iband)=dotr 2706 cj_h1vstrain_ci(2,istrpert,jband,iband)=doti 2707 2708 end do !jband 2709 2710 end do !iband 2711 2712 !Clean the atomic displacement rf_hamiltonian 2713 call rf_hamkq%free() 2714 2715 !Deallocations 2716 ABI_FREE(kpg_k) 2717 ABI_FREE(kpg1_k) 2718 ABI_FREE(dkinpw) 2719 ABI_FREE(kinpw1) 2720 ABI_FREE(ffnlk) 2721 ABI_FREE(ffnl1) 2722 ABI_FREE(ph3d) 2723 2724 end do !istrpert 2725 2726 !----------------------------------------------------------------------------------------------- 2727 ! Atomic displacement 1st order hamiltonian + 1st order potential matrix element: 2728 ! < u_{i,k}^{(0)} | (H^{\tau_{\kappa\alpha}}+V^{\tau_{\kappa\alpha}})^{\dagger} | u_{j,k}^{(0)} > 2729 ! calculated as 2730 ! (< u_{j,k}^{(0)} | H^{\tau_{\kappa\alpha}}+V^{\tau_{\kappa\alpha}} | u_{i,k}^{(0)} >)^* 2731 !----------------------------------------------------------------------------------------------- 2732 2733 !Specific definitions 2734 ipert=1 2735 useylmgr1=0 2736 dum_lambda=zero 2737 berryopt=0;optlocal=1;optnl=1;usevnl=0;opt_gvnl1=0;sij_opt=0 2738 2739 !Specific allocations 2740 ABI_MALLOC(ci_h1vatdisdag_cj,(2,natpert,nband_k,nband_k)) 2741 ABI_MALLOC(dum_ylmgr1_k,(npw_k,3+6*((ipert-dtset%natom)/10), psps%mpsang*psps%mpsang*psps%useylm*useylmgr1)) 2742 2743 !LOOP OVER ATOMIC DISPLACEMENT PERTURBATIONS 2744 do iatpert= 1, natpert 2745 ipert=pert_atdis(1,iatpert) 2746 idir=pert_atdis(2,iatpert) 2747 2748 !Get first-order local part of the pseudopotential 2749 call dfpt_vlocal(atindx,cplex,gs_hamkq%gmet,gsqcut,idir,ipert,mpi_enreg, & 2750 & psps%mqgrid_vl,dtset%natom,& 2751 & nattyp,nfft,ngfft,dtset%ntypat,ngfft(1),ngfft(2),ngfft(3), & 2752 & ph1d,psps%qgrid_vl,& 2753 & dtset%qptn,ucvol,psps%vlspl,vpsp1,xred) 2754 2755 !Set up local potential vlocal1 with proper dimensioning, from vpsp1 + vhxc1_atdis 2756 vpsp1=vpsp1+vhxc1_atdis(iatpert,:) 2757 call rf_transgrid_and_pack(isppol,nspden,psps%usepaw,cplex,nfft,dtset%nfft,dtset%ngfft,& 2758 & gs_hamkq%nvloc,pawfgr,mpi_enreg,dum_vpsp,vpsp1,dum_vlocal,vlocal1) 2759 2760 !Initialize rf_hamiltonian (the k-dependent part is prepared in getgh1c_setup) 2761 call init_rf_hamiltonian(cplex,gs_hamkq,ipert,rf_hamkq,& 2762 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab) 2763 call rf_hamkq%load_spin(isppol,vlocal1=vlocal1,with_nonlocal=.true.) 2764 2765 !Set up the ground-state Hamiltonian, and some parts of the 1st-order Hamiltonian 2766 call getgh1c_setup(gs_hamkq,rf_hamkq,dtset,psps,& ! In 2767 kpt,kpt,idir,ipert,dtset%natom,rmet,gs_hamkq%gprimd,gs_hamkq%gmet,istwf_k,& ! In 2768 npw_k,npw_k,useylmgr1,kg_k,ylm_k,kg_k,ylm_k,dum_ylmgr1_k,& ! In 2769 dkinpw,nkpg,nkpg1,kpg_k,kpg1_k,kinpw1,ffnlk,ffnl1,ph3d,ph3d1) ! Out 2770 2771 !LOOP OVER KET BANDS 2772 do iband=1,nband_k 2773 2774 if(mpi_enreg%proc_distrb(ikpt,iband,isppol) /= mpi_enreg%me_kpt) cycle 2775 2776 !Read ket ground-state wavefunctions 2777 cwave0i(:,:)=cg(:,1+(iband-1)*npw_k*dtset%nspinor+icg:iband*npw_k*dtset%nspinor+icg) 2778 2779 !Compute < g | H^{\tau_{\kappa\beta}}+V^{\tau_{\kappa\beta}} | u_{i,k}^{(0)} > 2780 call getgh1c(berryopt,cwave0i,dum_cwaveprj,gv1c,dum_grad_berry,& 2781 & dum_gs1,gs_hamkq,dum_gvnl1,idir,ipert,dum_lambda,mpi_enreg,optlocal,& 2782 & optnl,opt_gvnl1,rf_hamkq,sij_opt,tim_getgh1c,usevnl) 2783 2784 !LOOP OVER BRA BANDS 2785 do jband=1,nband_k 2786 2787 !Read bra ground-state wavefunctions 2788 cwave0j(:,:)=cg(:,1+(jband-1)*npw_k*dtset%nspinor+icg:jband*npw_k*dtset%nspinor+icg) 2789 2790 !Compute = ci_h1vatdisdag_cj 2791 !(< u_{j,k}^{(0)} | H^{\tau_{\kappa\beta}}+V^{\tau_{\kappa\beta}} | u_{i,k}^{(0)} >)^* 2792 call dotprod_g(dotr,doti,istwf_k,npw_k*dtset%nspinor,2,cwave0j,gv1c,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 2793 ci_h1vatdisdag_cj(1,iatpert,jband,iband)=dotr 2794 ci_h1vatdisdag_cj(2,iatpert,jband,iband)=-doti 2795 2796 end do !jband 2797 2798 end do !iband 2799 2800 !Clean the atomic displacement rf_hamiltonian 2801 call rf_hamkq%free() 2802 2803 !Deallocations 2804 ABI_FREE(kpg_k) 2805 ABI_FREE(kpg1_k) 2806 ABI_FREE(dkinpw) 2807 ABI_FREE(kinpw1) 2808 ABI_FREE(ffnlk) 2809 ABI_FREE(ffnl1) 2810 ABI_FREE(ph3d) 2811 2812 end do !iatpert 2813 2814 !Deallocations 2815 ABI_FREE(dum_ylmgr1_k) 2816 ABI_FREE(cwave0j) 2817 ABI_FREE(vpsp1) 2818 ABI_FREE(vlocal1) 2819 2820 2821 !---------------------------------------------------------------------------------------- 2822 ! Terms that involve first order response functions 2823 !---------------------------------------------------------------------------------------- 2824 2825 !Allocation of bks (band, k-point and spin) dependent terms 2826 ABI_MALLOC(c1atdis_q1gradH0_c1strain_bks,(2,nband_k,natpert,nq1grad,nstrpert)) 2827 ABI_MALLOC(c1atdis_dQcalHstrain_c0_bks,(2,nband_k,natpert,nq1grad,nstrpert)) 2828 ABI_MALLOC(c0_HatdisdagdQ_c1strain_bks,(2,nband_k,natpert,nq1grad,nstrpert)) 2829 ABI_MALLOC(c0_Hatdisdqdag_c1strain_bks,(2,nband_k,natpert,nq1grad,nstrpert)) 2830 ABI_MALLOC(c1atdis_Hmetricdqdq_c0_bks,(2,nband_k,natpert,nq1grad,nstrpert)) 2831 c1atdis_dQcalHstrain_c0_bks=zero 2832 c0_HatdisdagdQ_c1strain_bks=zero 2833 2834 !Allocation of wf1s 2835 ABI_MALLOC(cg1_strain,(2,npw_k*dtset%nspinor)) 2836 ABI_MALLOC(cg1_strain_ar,(3,3,2,npw_k*dtset%nspinor)) 2837 ABI_MALLOC(cg1_atdis,(2,npw_k*dtset%nspinor)) 2838 ABI_MALLOC(cg1_ddk,(2,npw_k*dtset%nspinor)) 2839 2840 !Check correspondance with the data in wf1 files 2841 !Atomic displacements 2842 do iatpert=1,natpert 2843 !k-point index check 2844 ii = wfk_t_atdis(iatpert)%findk( kpt(:)) 2845 ABI_CHECK(ii == indkpt1(ikpt), "ii != indkpt in atomic displacement wf1 file") 2846 !npw check 2847 npw_disk = wfk_t_atdis(iatpert)%hdr%npwarr(ii) 2848 if (npw_k /= npw_disk) then 2849 write(unit=msg,fmt='(a,i5,a,i5,a,a,i5,a,a,i5)')& 2850 & 'For ikpt = ',ikpt,', and pertcase= ',pert_atdis(3,iatpert),ch10,& 2851 & 'the number of plane waves in the wf1 file is equal to', npw_disk,ch10,& 2852 & 'while it should be ',npw_k 2853 ABI_BUG(msg) 2854 end if 2855 end do 2856 2857 !ddk 2858 do iq1grad=1,nq1grad 2859 !k-point index check 2860 ii = wfk_t_ddk(iq1grad)%findk( kpt(:)) 2861 ABI_CHECK(ii == indkpt1(ikpt), "ii != indkpt1 in ddk wf1 file") 2862 !npw check 2863 npw_disk = wfk_t_ddk(iq1grad)%hdr%npwarr(ii) 2864 if (npw_k /= npw_disk) then 2865 write(unit=msg,fmt='(a,i5,a,i5,a,a,i5,a,a,i5)')& 2866 & 'For ikpt = ',ikpt,', and pertcase= ',q1grad(3,iq1grad),ch10,& 2867 & 'the number of plane waves in the wf1 file is equal to', npw_disk,ch10,& 2868 & 'while it should be ',npw_k 2869 ABI_BUG(msg) 2870 end if 2871 end do 2872 2873 !strain 2874 do istrpert=1,nstrpert 2875 ka=pert_strain(3,istrpert) 2876 kb=pert_strain(4,istrpert) 2877 !k-point index check 2878 ii = wfk_t_strain(ka,kb)%findk( kpt(:)) 2879 ABI_CHECK(ii == indkpt1(ikpt), "ii != indkpt1 in strain wf1 file") 2880 !npw check 2881 npw_disk = wfk_t_strain(ka,kb)%hdr%npwarr(ii) 2882 if (npw_k /= npw_disk) then 2883 write(unit=msg,fmt='(a,i5,a,i5,a,a,i5,a,a,i5)')& 2884 & 'For ikpt = ',ikpt,', and pertcase= ',pert_strain(5,istrpert),ch10,& 2885 & 'the number of plane waves in the wf1 file is equal to', npw_disk,ch10,& 2886 & 'while it should be ',npw_k 2887 ABI_BUG(msg) 2888 end if 2889 end do 2890 2891 !-------------------------------------------------------------------------------------- 2892 !q1-gradient of gs Hamiltonian: 2893 ! < u_{i,k}^{\tau_{\kappa\alpha}} | \partial_{gamma} H^{(0)} | u_{i,k}^{n_{\beta\delta}} > 2894 !-------------------------------------------------------------------------------------- 2895 2896 !Specific allocations 2897 ABI_MALLOC(vlocal1,(cplex*ngfft(4),ngfft(5),ngfft(6),gs_hamkq%nvloc)) 2898 2899 !Specific definitions 2900 vlocal1=zero 2901 useylmgr1=1 2902 dum_lambda=zero 2903 berryopt=0;optlocal=0;optnl=1;usevnl=0;opt_gvnl1=0;sij_opt=0 2904 2905 !LOOP OVER Q1-GRADIENT 2906 do iq1grad=1,nq1grad 2907 ipert=q1grad(1,iq1grad) 2908 idir=q1grad(2,iq1grad) 2909 2910 !Initializes rf_hamiltonian (the k-dependent part is prepared in getgh1c_setup) 2911 call init_rf_hamiltonian(cplex,gs_hamkq,ipert,rf_hamkq,& 2912 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab) 2913 call rf_hamkq%load_spin(isppol,vlocal1=vlocal1,with_nonlocal=.true.) 2914 2915 !Set up the ground-state Hamiltonian, and some parts of the 1st-order Hamiltonian 2916 call getgh1c_setup(gs_hamkq,rf_hamkq,dtset,psps,& ! In 2917 kpt,kpt,idir,ipert,dtset%natom,rmet,gs_hamkq%gprimd,gs_hamkq%gmet,istwf_k,& ! In 2918 npw_k,npw_k,useylmgr1,kg_k,ylm_k,kg_k,ylm_k,part_ylmgr_k,& ! In 2919 dkinpw,nkpg,nkpg1,kpg_k,kpg1_k,kinpw1,ffnlk,ffnl1,ph3d,ph3d1) ! Out 2920 2921 !LOOP OVER BANDS 2922 do iband=1,nband_k 2923 2924 if(mpi_enreg%proc_distrb(ikpt,iband,isppol) /= mpi_enreg%me_kpt) cycle 2925 2926 !LOOP OVER STRAIN PERTURBATION 2927 do istrpert=1,nstrpert 2928 ka=pert_strain(3,istrpert) 2929 kb=pert_strain(4,istrpert) 2930 2931 !Read strain field wf1 2932 if (ka>=kb) then 2933 call wfk_t_strain(ka,kb)%read_bks( iband, indkpt1(ikpt), & 2934 & isppol, xmpio_single, cg_bks=cg1_strain) 2935 cg1_strain_ar(ka,kb,:,:)=cg1_strain 2936 else 2937 cg1_strain=cg1_strain_ar(kb,ka,:,:) 2938 end if 2939 2940 !Compute < g |\partial_{gamma} H^{(0)} | u_{i,k}^{n_{\beta\delta}} > 2941 call getgh1c(berryopt,cg1_strain,dum_cwaveprj,gv1c,dum_grad_berry,& 2942 & dum_gs1,gs_hamkq,dum_gvnl1,idir,ipert,dum_lambda,mpi_enreg,optlocal,& 2943 & optnl,opt_gvnl1,rf_hamkq,sij_opt,tim_getgh1c,usevnl) 2944 2945 !LOOP OVER ATOMIC DISPLACEMENT PERTURBATION 2946 do iatpert=1,natpert 2947 2948 !Read atomic displacement wf1 2949 call wfk_t_atdis(iatpert)%read_bks( iband, indkpt1(ikpt), & 2950 & isppol, xmpio_single, cg_bks=cg1_atdis) 2951 2952 !calculate: 2953 ! < u_{i,k}^{\tau_{\kappa\alpha}} | \partial_{gamma} H^{(0)} | u_{i,k}^{n_{\beta\delta}} > 2954 call dotprod_g(dotr,doti,istwf_k,npw_k*dtset%nspinor,2,cg1_atdis,gv1c, & 2955 & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 2956 2957 c1atdis_q1gradH0_c1strain_bks(1,iband,iatpert,iq1grad,istrpert)=dotr 2958 c1atdis_q1gradH0_c1strain_bks(2,iband,iatpert,iq1grad,istrpert)=doti 2959 2960 end do !iefipert 2961 2962 end do !istrpert 2963 2964 end do !iband 2965 2966 !Clean the ddk rf_hamiltonian 2967 call rf_hamkq%free() 2968 2969 !Deallocations 2970 ABI_FREE(kpg_k) 2971 ABI_FREE(kpg1_k) 2972 ABI_FREE(dkinpw) 2973 ABI_FREE(kinpw1) 2974 ABI_FREE(ffnlk) 2975 ABI_FREE(ffnl1) 2976 ABI_FREE(ph3d) 2977 2978 end do !iq1grad 2979 2980 !Deallocations 2981 ABI_FREE(gv1c) 2982 ABI_FREE(vlocal1) 2983 !ABI_FREE(ph3d1) !it is only allocated if kpt and kpq are different. Not the case. 2984 2985 !-------------------------------------------------------------------------------------- 2986 !Two terms involving a q-gradient of projector operators 2987 !-------------------------------------------------------------------------------------- 2988 !LOOP OVER BANDS 2989 do iband=1,nband_k 2990 2991 if(mpi_enreg%proc_distrb(ikpt,iband,isppol) /= mpi_enreg%me_kpt) cycle 2992 2993 !LOOP OVER ATOMIC DISPLACEMENT PERTURBATION 2994 do iatpert=1,natpert 2995 2996 !Read atomic displacement wf1 2997 call wfk_t_atdis(iatpert)%read_bks( iband, indkpt1(ikpt), & 2998 & isppol, xmpio_single, cg_bks=cg1_atdis) 2999 3000 !LOOP OVER q1-GRADIENT 3001 do iq1grad=1,nq1grad 3002 3003 !LOOP OVER BANDS 3004 do jband=1,nband_k 3005 3006 !Read ddk wf1 3007 call wfk_t_ddk(iq1grad)%read_bks( jband, indkpt1(ikpt), & 3008 & isppol, xmpio_single, cg_bks=cg1_ddk) 3009 3010 !Calculate: < u_{i,k}^{\tau_{\kappa\alpha}} | u_{j,k}^{k_{\gamma}} > 3011 call dotprod_g(dotr,doti,istwf_k,npw_k*dtset%nspinor,2,cg1_atdis,cg1_ddk, & 3012 & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 3013 3014 !LOOP OVER STRAIN PERTURBATION 3015 do istrpert=1,nstrpert 3016 3017 !Calculate: -\sum_{j} < u_{i,k}^{\tau_{\kappa\alpha}} | u_{j,k}^{k_{\gamma}} > * 3018 ! < u_{j,k}^{(0)} | H^{n_{\beta\delta}}+V^{n_{\beta\delta}} | u_{i,k}^{(0)} > 3019 cprodr=dotr*cj_h1vstrain_ci(1,istrpert,jband,iband) - & 3020 & doti*cj_h1vstrain_ci(2,istrpert,jband,iband) 3021 cprodi=dotr*cj_h1vstrain_ci(2,istrpert,jband,iband) + & 3022 & doti*cj_h1vstrain_ci(1,istrpert,jband,iband) 3023 3024 c1atdis_dQcalHstrain_c0_bks(1,iband,iatpert,iq1grad,istrpert)= & 3025 & c1atdis_dQcalHstrain_c0_bks(1,iband,iatpert,iq1grad,istrpert)-cprodr 3026 c1atdis_dQcalHstrain_c0_bks(2,iband,iatpert,iq1grad,istrpert)= & 3027 & c1atdis_dQcalHstrain_c0_bks(2,iband,iatpert,iq1grad,istrpert)-cprodi 3028 3029 end do !istrpert 3030 3031 end do !jband 3032 3033 end do !iq1grad 3034 3035 end do !iatpert 3036 3037 !LOOP OVER STRAIN PERTURBATION 3038 do istrpert=1,nstrpert 3039 ka=pert_strain(3,istrpert) 3040 kb=pert_strain(4,istrpert) 3041 3042 !Read strain field wf1 3043 if (ka>=kb) then 3044 call wfk_t_strain(ka,kb)%read_bks( iband, indkpt1(ikpt), & 3045 & isppol, xmpio_single, cg_bks=cg1_strain) 3046 cg1_strain_ar(ka,kb,:,:)=cg1_strain 3047 else 3048 cg1_strain=cg1_strain_ar(kb,ka,:,:) 3049 end if 3050 3051 !LOOP OVER q1-GRADIENT 3052 do iq1grad=1,nq1grad 3053 3054 !LOOP OVER BANDS 3055 do jband=1,nband_k 3056 3057 !Read ddk wf1 3058 call wfk_t_ddk(iq1grad)%read_bks( jband, indkpt1(ikpt), & 3059 & isppol, xmpio_single, cg_bks=cg1_ddk) 3060 3061 !Calculate: < u_{j,k}^{k_{\gamma}} | u_{i,k}^{n_{\beta\delta}} > 3062 call dotprod_g(dotr,doti,istwf_k,npw_k*dtset%nspinor,2,cg1_ddk,cg1_strain, & 3063 & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 3064 3065 !LOOP OVER ATOMIC DISPLACEMENT PERTURBATION 3066 do iatpert=1,natpert 3067 3068 !Calculate: -\sum_{j} 3069 ! < u_{i,k}^{(0)} | (H^{\tau_{\kappa\alpha}}+V^{\tau_{\kappa\alpha}})^{\dagger} | u_{j,k}^{(0)} > 3070 ! < u_{j,k}^{k_{\gamma}} | u_{i,k}^{n_{\beta\delta}} > 3071 cprodr=dotr*ci_h1vatdisdag_cj(1,iatpert,jband,iband) - & 3072 & doti*ci_h1vatdisdag_cj(2,iatpert,jband,iband) 3073 cprodi=dotr*ci_h1vatdisdag_cj(2,iatpert,jband,iband) + & 3074 & doti*ci_h1vatdisdag_cj(1,iatpert,jband,iband) 3075 3076 c0_HatdisdagdQ_c1strain_bks(1,iband,iatpert,iq1grad,istrpert)= & 3077 & c0_HatdisdagdQ_c1strain_bks(1,iband,iatpert,iq1grad,istrpert)-cprodr 3078 c0_HatdisdagdQ_c1strain_bks(2,iband,iatpert,iq1grad,istrpert)= & 3079 & c0_HatdisdagdQ_c1strain_bks(2,iband,iatpert,iq1grad,istrpert)-cprodi 3080 3081 end do !iatpert 3082 3083 end do !jband 3084 3085 end do !iq1grad 3086 3087 end do !istrpert 3088 3089 end do !iband 3090 3091 !Deallocations 3092 ABI_FREE(ci_h1vatdisdag_cj) 3093 ABI_FREE(cj_h1vstrain_ci) 3094 ABI_FREE(cg1_ddk) 3095 3096 !-------------------------------------------------------------------------------------- 3097 !Two terms involving a q-gradient of first order Hamiltonians 3098 !-------------------------------------------------------------------------------------- 3099 3100 !-------------------------------------------------------------------------------------- 3101 ! q1-gradient of the strain first order Hamiltonian: 3102 ! i/2 < u_{i,k}^{\tau_{\kappa\alpha} | H^{n_{\beta\delta}}_{\gamma} | u_{i,k}^{(0)} > 3103 ! Computed as the second order q-gradient of the first order metric Hamiltonian: 3104 ! 1/2 < u_{i,k}^{\tau_{\kappa\alpha} | H^{(\beta)}_{\gamma\delta} | u_{i,k}^{(0)} > 3105 !-------------------------------------------------------------------------------------- 3106 3107 !Specific allocations 3108 ABI_MALLOC(vhart1dqdq,(2*nfft)) 3109 ABI_MALLOC(vpsp1dqdq,(2*nfft)) 3110 ABI_MALLOC(vxc1dqdq,(2*nfft)) 3111 ABI_MALLOC(vlocal1dqdq,(2*ngfft(4),ngfft(5),ngfft(6),gs_hamkq%nvloc)) 3112 ABI_MALLOC(gh1dqdqc,(2,npw_k*dtset%nspinor)) 3113 ABI_MALLOC(gvloc1dqdqc,(2,npw_k*dtset%nspinor)) 3114 ABI_MALLOC(gvnl1dqdqc,(2,npw_k*dtset%nspinor)) 3115 3116 !Specific definitions 3117 useylmgr1=1;optlocal=1;optnl=1;opthartdqdq=1; 3118 3119 !LOOP OVER STRAIN PERTURBATIONS 3120 do istrpert=1,nstrpert 3121 ipert=pert_strain(1,istrpert) 3122 idir=pert_strain(6,istrpert) 3123 3124 !LOOP OVER Q1-GRADIENT 3125 do iq1grad=1,nq1grad 3126 3127 !Get 2nd q-gradient of first-order local part of the pseudopotential and of the Hartree 3128 !(and XC if GGA) contribution from ground state density 3129 call dfpt_vmetdqdq(2,gs_hamkq%gmet,gs_hamkq%gprimd,gsqcut,idir,ipert,kxc,mpi_enreg, & 3130 & psps%mqgrid_vl,dtset%natom, & 3131 & nattyp,nfft,ngfft,dtset%ntypat,ngfft(1),ngfft(2),ngfft(3),nkxc,nspden,opthartdqdq, & 3132 & ph1d,q1grad(2,iq1grad),psps%qgrid_vl,& 3133 & dtset%qptn,rhog,rhor,ucvol,psps%vlspl,vhart1dqdq,vpsp1dqdq,vxc1dqdq) 3134 3135 !Merge the local contributions 3136 vpsp1dqdq=vpsp1dqdq+vhart1dqdq+vxc1dqdq 3137 3138 !Set up q-gradient of strain potential vlocal1dqdq with proper dimensioning 3139 call rf_transgrid_and_pack(isppol,nspden,psps%usepaw,2,nfft,dtset%nfft,dtset%ngfft,& 3140 & gs_hamkq%nvloc,pawfgr,mpi_enreg,dum_vpsp,vpsp1dqdq,dum_vlocal,vlocal1dqdq) 3141 3142 !Initialize rf_hamiltonian (the k-dependent part is prepared in getgh1c_setup) 3143 call init_rf_hamiltonian(2,gs_hamkq,ipert,rf_hamkq,& 3144 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab) 3145 call rf_hamkq%load_spin(isppol,vlocal1=vlocal1dqdq,with_nonlocal=.true.) 3146 3147 !Set up the ground-state Hamiltonian, and some parts of the 1st-order Hamiltonian 3148 call getgh1dqc_setup(gs_hamkq,rf_hamkq,dtset,psps,kpt,kpt,idir,ipert,q1grad(2,iq1grad), & 3149 & dtset%natom,rmet,gs_hamkq%gprimd,gs_hamkq%gmet,istwf_k,npw_k,npw_k,nylmgr,useylmgr1,kg_k, & 3150 & ylm_k,kg_k,ylm_k,ylmgr_k,nkpg,nkpg1,kpg_k,kpg1_k,dkinpw,kinpw1,ffnlk,ffnl1,ph3d,ph3d1) 3151 3152 !LOOP OVER BANDS 3153 do iband=1,nband_k 3154 3155 if(mpi_enreg%proc_distrb(ikpt,iband,isppol) /= mpi_enreg%me_kpt) cycle 3156 3157 !Read ket ground-state wavefunctions 3158 cwave0i(:,:)=cg(:,1+(iband-1)*npw_k*dtset%nspinor+icg:iband*npw_k*dtset%nspinor+icg) 3159 3160 !Compute < g | H^{(\beta)}_{\gamma\delta} | u_{i,k}^{(0)} > 3161 call getgh1dqc(cwave0i,dum_cwaveprj,gh1dqdqc,gvloc1dqdqc,gvnl1dqdqc,gs_hamkq, & 3162 & idir,ipert,mpi_enreg,optlocal,optnl,q1grad(2,iq1grad),rf_hamkq) 3163 3164 !LOOP OVER ATOMIC DISPLACEMENT PERTURBATION 3165 do iatpert=1,natpert 3166 3167 !Read atomic displacement wf1 3168 call wfk_t_atdis(iatpert)%read_bks( iband, indkpt1(ikpt), & 3169 & isppol, xmpio_single, cg_bks=cg1_atdis) 3170 3171 !calculate: 3172 ! < u_{i,k}^{\tau_{\kappa\alpha}} | H^{(\beta)}_{\gamma\delta} | u_{i,k}^{(0)} > 3173 call dotprod_g(dotr,doti,istwf_k,npw_k*dtset%nspinor,2,cg1_atdis,gh1dqdqc, & 3174 & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 3175 3176 c1atdis_Hmetricdqdq_c0_bks(1,iband,iatpert,iq1grad,istrpert)=half*dotr 3177 c1atdis_Hmetricdqdq_c0_bks(2,iband,iatpert,iq1grad,istrpert)=half*doti 3178 3179 end do !iatpert 3180 3181 end do !iband 3182 3183 !Clean the rf_hamiltonian 3184 call rf_hamkq%free() 3185 3186 !Deallocations 3187 ABI_FREE(kpg_k) 3188 ABI_FREE(kpg1_k) 3189 ABI_FREE(dkinpw) 3190 ABI_FREE(kinpw1) 3191 ABI_FREE(ffnlk) 3192 ABI_FREE(ffnl1) 3193 ABI_FREE(ph3d) 3194 3195 end do !iq1grad 3196 3197 end do !istrpert 3198 3199 !Deallocations 3200 ABI_FREE(cg1_atdis) 3201 ABI_FREE(vhart1dqdq) 3202 ABI_FREE(gh1dqdqc) 3203 ABI_FREE(gvloc1dqdqc) 3204 ABI_FREE(gvnl1dqdqc) 3205 ABI_FREE(vpsp1dqdq) 3206 ABI_FREE(vlocal1dqdq) 3207 ABI_FREE(vxc1dqdq) 3208 3209 !-------------------------------------------------------------------------------------- 3210 ! q1-gradient of atomic displacement 1st-order Hamiltonian: 3211 ! < u_{i,k}^{(0)} | (H^{\tau_{\kappa\alpha}}_{\gamma})^{dagger} | u_{i,k}^{n_{\beta\delta}} > 3212 ! calculated as: 3213 ! (<u_{i,k}^{n_{\beta\delta} | H^{\tau_{\kappa\alpha}}_{\gamma} | u_{i,k}^{(0)} >)* 3214 !-------------------------------------------------------------------------------------- 3215 3216 !Specific allocations 3217 ABI_MALLOC(vpsp1dq,(2*nfft)) 3218 ABI_MALLOC(vlocal1dq,(2*ngfft(4),ngfft(5),ngfft(6),gs_hamkq%nvloc)) 3219 ABI_MALLOC(gh1dqc,(2,npw_k*dtset%nspinor)) 3220 ABI_MALLOC(gvloc1dqc,(2,npw_k*dtset%nspinor)) 3221 ABI_MALLOC(gvnl1dqc,(2,npw_k*dtset%nspinor)) 3222 3223 !Specific definitions 3224 useylmgr1=1;nylmgrpart=3;optlocal=1;optnl=1 3225 3226 !LOOP OVER ATOMIC DISPLACEMENT PERTURBATIONS 3227 do iatpert= 1, natpert 3228 ipert=pert_atdis(1,iatpert) 3229 idir=pert_atdis(2,iatpert) 3230 3231 !LOOP OVER Q1-GRADIENT 3232 do iq1grad=1,nq1grad 3233 3234 !Get q-gradient of first-order local part of the pseudopotential 3235 call dfpt_vlocaldq(atindx,2,gs_hamkq%gmet,gsqcut,idir,ipert,mpi_enreg, & 3236 & psps%mqgrid_vl,dtset%natom,& 3237 & nattyp,nfft,ngfft,dtset%ntypat,ngfft(1),ngfft(2),ngfft(3), & 3238 & ph1d,q1grad(2,iq1grad),psps%qgrid_vl,& 3239 & dtset%qptn,ucvol,psps%vlspl,vpsp1dq) 3240 3241 !Set up q-gradient of local potential vlocal1dq with proper dimensioning 3242 call rf_transgrid_and_pack(isppol,nspden,psps%usepaw,2,nfft,dtset%nfft,dtset%ngfft,& 3243 & gs_hamkq%nvloc,pawfgr,mpi_enreg,dum_vpsp,vpsp1dq,dum_vlocal,vlocal1dq) 3244 3245 !Initialize rf_hamiltonian (the k-dependent part is prepared in getgh1c_setup) 3246 call init_rf_hamiltonian(2,gs_hamkq,ipert,rf_hamkq,& 3247 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab) 3248 call rf_hamkq%load_spin(isppol,vlocal1=vlocal1dq,with_nonlocal=.true.) 3249 3250 !Set up the ground-state Hamiltonian, and some parts of the 1st-order Hamiltonian 3251 call getgh1dqc_setup(gs_hamkq,rf_hamkq,dtset,psps,kpt,kpt,idir,ipert,q1grad(2,iq1grad), & 3252 & dtset%natom,rmet,gs_hamkq%gprimd,gs_hamkq%gmet,istwf_k,npw_k,npw_k,nylmgrpart,useylmgr1,kg_k, & 3253 & ylm_k,kg_k,ylm_k,part_ylmgr_k,nkpg,nkpg1,kpg_k,kpg1_k,dkinpw,kinpw1,ffnlk,ffnl1,ph3d,ph3d1) 3254 3255 !LOOP OVER BANDS 3256 do iband=1,nband_k 3257 3258 if(mpi_enreg%proc_distrb(ikpt,iband,isppol) /= mpi_enreg%me_kpt) cycle 3259 3260 !Read ket ground-state wavefunctions 3261 cwave0i(:,:)=cg(:,1+(iband-1)*npw_k*dtset%nspinor+icg:iband*npw_k*dtset%nspinor+icg) 3262 3263 !Compute < g |H^{\tau_{\kappa\beta}}_{\gamma} | u_{i,k}^{(0)} > 3264 call getgh1dqc(cwave0i,dum_cwaveprj,gh1dqc,gvloc1dqc,gvnl1dqc,gs_hamkq, & 3265 & idir,ipert,mpi_enreg,optlocal,optnl,q1grad(2,iq1grad),rf_hamkq) 3266 3267 !LOOP OVER STRAIN PERTURBATION 3268 do istrpert=1,nstrpert 3269 ka=pert_strain(3,istrpert) 3270 kb=pert_strain(4,istrpert) 3271 3272 !Read strain field wf1 3273 if (ka>=kb) then 3274 call wfk_t_strain(ka,kb)%read_bks( iband, indkpt1(ikpt), & 3275 & isppol, xmpio_single, cg_bks=cg1_strain) 3276 cg1_strain_ar(ka,kb,:,:)=cg1_strain 3277 else 3278 cg1_strain=cg1_strain_ar(kb,ka,:,:) 3279 end if 3280 3281 !Calculate: 3282 !(<u_{i,k}^{n_{\beta\delta} | H^{\tau_{\kappa\alpha}}_{\gamma} | u_{i,k}^{(0)} >)* 3283 call dotprod_g(dotr,doti,istwf_k,npw_k*dtset%nspinor,2,cg1_strain,gh1dqc, & 3284 & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 3285 3286 !Apply the -i factor that has been factorized in the 3287 !H^{\tau_{\kappa\beta}}_{\gamma} terms. (And consider the complex 3288 !conjugate too) 3289 c0_Hatdisdqdag_c1strain_bks(1,iband,iatpert,iq1grad,istrpert)= doti 3290 c0_Hatdisdqdag_c1strain_bks(2,iband,iatpert,iq1grad,istrpert)= dotr 3291 3292 end do !istrpert 3293 3294 end do !iband 3295 3296 !Clean the rf_hamiltonian 3297 call rf_hamkq%free() 3298 3299 !Deallocations 3300 ABI_FREE(kpg_k) 3301 ABI_FREE(kpg1_k) 3302 ABI_FREE(dkinpw) 3303 ABI_FREE(kinpw1) 3304 ABI_FREE(ffnlk) 3305 ABI_FREE(ffnl1) 3306 ABI_FREE(ph3d) 3307 3308 end do !iq1grad 3309 3310 end do !iatpert 3311 3312 !Deallocations 3313 ABI_FREE(cg1_strain) 3314 ABI_FREE(cg1_strain_ar) 3315 ABI_FREE(part_ylmgr_k) 3316 ABI_FREE(dum_cwaveprj) 3317 ABI_FREE(gh1dqc) 3318 ABI_FREE(gvloc1dqc) 3319 ABI_FREE(gvnl1dqc) 3320 ABI_FREE(vpsp1dq) 3321 ABI_FREE(vlocal1dq) 3322 ABI_FREE(dum_vpsp) 3323 ABI_FREE(dum_vlocal) 3324 3325 !-------------------------------------------------------------------------------------- 3326 ! Acumulate all the wf dependent terms of the flexoelectric tensor 3327 !-------------------------------------------------------------------------------------- 3328 isdqwf_k=zero 3329 isdqwf_t1_k=zero 3330 isdqwf_t2_k=zero 3331 isdqwf_t3_k=zero 3332 isdqwf_t4_k=zero 3333 isdqwf_t5_k=zero 3334 do istrpert=1,nstrpert 3335 ka=pert_strain(3,istrpert) 3336 kb=pert_strain(4,istrpert) 3337 do iq1grad=1,nq1grad 3338 do iatpert=1,natpert 3339 iatom=pert_atdis(1,iatpert) 3340 iatdir=pert_atdis(2,iatpert) 3341 do iband=1,nband_k 3342 3343 if(mpi_enreg%proc_distrb(ikpt,iband,isppol) /= mpi_enreg%me_kpt) cycle 3344 3345 !All terms toghether except T4 that needs further treatment 3346 isdqwf_k(1,iatom,iatdir,iq1grad,ka,kb)=isdqwf_k(1,iatom,iatdir,iq1grad,ka,kb) + & 3347 & occ_k(iband) * ( c1atdis_q1gradH0_c1strain_bks(1,iband,iatpert,iq1grad,istrpert) + & !T1 3348 & c1atdis_dQcalHstrain_c0_bks(1,iband,iatpert,iq1grad,istrpert) + & !T2 3349 & c0_HatdisdagdQ_c1strain_bks(1,iband,iatpert,iq1grad,istrpert) + & !T3 3350 & c0_Hatdisdqdag_c1strain_bks(1,iband,iatpert,iq1grad,istrpert) ) !T5 3351 isdqwf_k(2,iatom,iatdir,iq1grad,ka,kb)=isdqwf_k(2,iatom,iatdir,iq1grad,ka,kb) + & 3352 & occ_k(iband) * ( c1atdis_q1gradH0_c1strain_bks(2,iband,iatpert,iq1grad,istrpert) + & !T1 3353 & c1atdis_dQcalHstrain_c0_bks(2,iband,iatpert,iq1grad,istrpert) + & !T2 3354 & c0_HatdisdagdQ_c1strain_bks(2,iband,iatpert,iq1grad,istrpert) + & !T3 3355 & c0_Hatdisdqdag_c1strain_bks(2,iband,iatpert,iq1grad,istrpert) ) !T5 3356 3357 !Separate them 3358 !T1 3359 isdqwf_t1_k(1,iatom,iatdir,iq1grad,ka,kb)=isdqwf_t1_k(1,iatom,iatdir,iq1grad,ka,kb) + & 3360 & occ_k(iband) * c1atdis_q1gradH0_c1strain_bks(1,iband,iatpert,iq1grad,istrpert) 3361 3362 isdqwf_t1_k(2,iatom,iatdir,iq1grad,ka,kb)=isdqwf_t1_k(2,iatom,iatdir,iq1grad,ka,kb) + & 3363 & occ_k(iband) * c1atdis_q1gradH0_c1strain_bks(2,iband,iatpert,iq1grad,istrpert) 3364 3365 !T2 3366 isdqwf_t2_k(1,iatom,iatdir,iq1grad,ka,kb)=isdqwf_t2_k(1,iatom,iatdir,iq1grad,ka,kb) + & 3367 & occ_k(iband) * c1atdis_dQcalHstrain_c0_bks(1,iband,iatpert,iq1grad,istrpert) 3368 3369 isdqwf_t2_k(2,iatom,iatdir,iq1grad,ka,kb)=isdqwf_t2_k(2,iatom,iatdir,iq1grad,ka,kb) + & 3370 & occ_k(iband) * c1atdis_dQcalHstrain_c0_bks(2,iband,iatpert,iq1grad,istrpert) 3371 3372 !T3 3373 isdqwf_t3_k(1,iatom,iatdir,iq1grad,ka,kb)=isdqwf_t3_k(1,iatom,iatdir,iq1grad,ka,kb) + & 3374 & occ_k(iband) * c0_HatdisdagdQ_c1strain_bks(1,iband,iatpert,iq1grad,istrpert) 3375 3376 isdqwf_t3_k(2,iatom,iatdir,iq1grad,ka,kb)=isdqwf_t3_k(2,iatom,iatdir,iq1grad,ka,kb) + & 3377 & occ_k(iband) * c0_HatdisdagdQ_c1strain_bks(2,iband,iatpert,iq1grad,istrpert) 3378 3379 !T4 has type-I ordering for the array indexes 3380 isdqwf_t4_k(1,iatom,iatdir,ka,kb,iq1grad)=isdqwf_t4_k(1,iatom,iatdir,ka,kb,iq1grad) + & 3381 & occ_k(iband) * c1atdis_Hmetricdqdq_c0_bks(1,iband,iatpert,iq1grad,istrpert) 3382 3383 isdqwf_t4_k(2,iatom,iatdir,ka,kb,iq1grad)=isdqwf_t4_k(2,iatom,iatdir,ka,kb,iq1grad) + & 3384 & occ_k(iband) * c1atdis_Hmetricdqdq_c0_bks(2,iband,iatpert,iq1grad,istrpert) 3385 3386 !T5 3387 isdqwf_t5_k(1,iatom,iatdir,iq1grad,ka,kb)=isdqwf_t5_k(1,iatom,iatdir,iq1grad,ka,kb) + & 3388 & occ_k(iband) * c0_Hatdisdqdag_c1strain_bks(1,iband,iatpert,iq1grad,istrpert) 3389 3390 isdqwf_t5_k(2,iatom,iatdir,iq1grad,ka,kb)=isdqwf_t5_k(2,iatom,iatdir,iq1grad,ka,kb) + & 3391 & occ_k(iband) * c0_Hatdisdqdag_c1strain_bks(2,iband,iatpert,iq1grad,istrpert) 3392 3393 end do 3394 end do 3395 end do 3396 end do 3397 3398 !scale by the k-point weight 3399 isdqwf_t1_k=isdqwf_t1_k * wtk_k 3400 isdqwf_t2_k=isdqwf_t2_k * wtk_k 3401 isdqwf_t3_k=isdqwf_t3_k * wtk_k 3402 isdqwf_t4_k=isdqwf_t4_k * wtk_k 3403 isdqwf_t5_k=isdqwf_t5_k * wtk_k 3404 isdqwf_k=isdqwf_k * wtk_k 3405 3406 !Deallocations 3407 ABI_FREE(c1atdis_q1gradH0_c1strain_bks) 3408 ABI_FREE(c1atdis_dQcalHstrain_c0_bks) 3409 ABI_FREE(c0_HatdisdagdQ_c1strain_bks) 3410 ABI_FREE(c1atdis_Hmetricdqdq_c0_bks) 3411 ABI_FREE(c0_Hatdisdqdag_c1strain_bks) 3412 ABI_FREE(cwave0i) 3413 3414 DBG_EXIT("COLL") 3415 3416 end subroutine dfpt_isdqwf
ABINIT/dfpt_qdrpwf.f90 [ Functions ]
NAME
dfpt_qdrpwf
FUNCTION
This routine computes the band and kpt resolved contributions to the quadrupole tensor.
INPUTS
atindx(natom)=index table for atoms (see gstate.f) cg(2,mpw*nspinor*mband*mkmem*nsppol)=planewave coefficients of wavefunctions at k cplex: if 1, several magnitudes are REAL, if 2, COMPLEX dtset <type(dataset_type)>=all input variables for this dataset gs_hamkq <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k gsqcut=large sphere cut-off icg=shift to be applied on the location of data in the array cg ikpt=number of the k-point indkpt1(nkpt_rbz)=non-symmetrized indices of the k-points isppol=1 for unpolarized, 2 for spin-polarized istwf_k=parameter that describes the storage of wfs kg_k(3,npw_k)=reduced planewave coordinates. kpt(3)=reduced coordinates of k point mkmem =number of k points treated by this node mpi_enreg=information about MPI parallelization mpw=maximum dimensioned size of npw or wfs at k natpert=number of atomic displacement perturbations nattyp(ntypat)= # atoms of each type. nband_k=number of bands at this k point for that spin polarization nfft=(effective) number of FFT grid points (for this proc) ngfft(1:18)=integer array with FFT box dimensions and other nkpt_rbz= number of k-points in the RBZ npw_k=number of plane waves at this k point nq1grad=number of q1 (q_{\gamma}) gradients nq2grad=number of q2 (q_{\delta}) gradients nq1q2grad=number of q1q2 2nd order gradients nspden=number of spin-density components nsppol=1 for unpolarized, 2 for spin-polarized nylmgr=second dimension of ylmgr_k occ_k(nband_k)=occupation number for each band (usually 2) for each k. pert_atdis(3,natpert)=array with the info for the atomic displacement perturbations ph1d(2,3*(2*dtset%mgfft+1)*dtset%natom)=1-dimensional phases psps <type(pseudopotential_type)>=variables related to pseudopotentials q1grad(3,nq1grad)=array with the info for the q1 (q_{\gamma}) gradients q2grad(3,nq2grad)=array with the info for the q2 (q_{\gamma}) gradients q1q2grad(4,nq1q2grad)=array with the info for the q1q2 2nd order gradients rmet(3,3)=real space metric (bohr**2) ucvol=unit cell volume in bohr**3. useylmgr= if 1 use the derivative of spherical harmonics vhxc1_atdis(natpert,cplex*nfft)= electrostatic potential generated by a first order atomic displacement density vhxc1_efield(nq2grad,cplex*nfft)= electrostatic potential generated by a first order electric field density wfk_t_atdis(natpert)=unit numbers for the atomic displacement wf1 files wfk_t_efield(nq2grad)=unit numbers for the electric field wf1 files wfk_t_ddk(nq1grad)=unit numbers for the ddk wf1 files wfk_t_dkdk(nq1q2grad)=unit numbers for the dkdk wf1 files wtk_k=weight assigned to the k point. xred(3,natom)=reduced dimensionless atomic coordinates ylm_k(npw_k,psps%mpsang*psps%mpsang*psps%useylm)=real spherical harmonics for the k point ylmgr_k((npw_k,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr)= k-gradients of real spherical harmonics for the k point
OUTPUT
qdrpwf_k(2,natpert,nq2grad,nq1grad)= wave function dependent part of the quadrupole tensor for the k-point kpt qdrpwf_t1_k(2,natpert,nq2grad,nq1grad)= t1 term (see notes) of qdrpwf_k qdrpwf_t2_k(2,natpert,nq2grad,nq1grad)= t2 term (see notes) of qdrpwf_k qdrpwf_t3_k(2,natpert,nq2grad,nq1grad)= t3 term (see notes) of qdrpwf_k qdrpwf_t4_k(2,natpert,nq2grad,nq1grad)= t4 term (see notes) of qdrpwf_k qdrpwf_t5_k(2,natpert,nq2grad,nq1grad)= t5 term (see notes) of qdrpwf_k
SIDE EFFECTS
NOTES
SOURCE
140 subroutine dfpt_qdrpwf(atindx,cg,cplex,dtset,gs_hamkq,gsqcut,icg,ikpt,indkpt1,isppol,istwf_k, & 141 & kg_k,kpt,mkmem,mpi_enreg,mpw,natpert,nattyp,nband_k,nfft,ngfft,nkpt_rbz, & 142 & npw_k,nq1grad,nq2grad,nq1q2grad,nspden,nsppol,nylmgr,occ_k,pert_atdis,ph1d,psps,qdrpwf_k, & 143 & qdrpwf_t1_k,qdrpwf_t2_k,qdrpwf_t3_k,qdrpwf_t4_k,qdrpwf_t5_k,q1grad,q2grad,q1q2grad,& 144 & rmet,ucvol,useylmgr,vhxc1_atdis,vhxc1_efield,wfk_t_atdis,wfk_t_efield,wfk_t_ddk, & 145 & wfk_t_dkdk,wtk_k,xred,ylm_k,ylmgr_k) 146 147 !Arguments ------------------------------------ 148 !scalars 149 integer,intent(in) :: cplex,icg,ikpt,isppol,istwf_k 150 integer,intent(in) :: mkmem,mpw,natpert,nband_k,nfft 151 integer,intent(in) :: nkpt_rbz,npw_k,nq1grad,nq2grad,nq1q2grad,nspden,nsppol,nylmgr 152 integer,intent(in) :: useylmgr 153 real(dp),intent(in) :: gsqcut,ucvol,wtk_k 154 type(dataset_type),intent(in) :: dtset 155 type(gs_hamiltonian_type),intent(inout) :: gs_hamkq 156 type(MPI_type),intent(in) :: mpi_enreg 157 type(pseudopotential_type),intent(in) :: psps 158 159 !arrays 160 integer,intent(in) :: atindx(dtset%natom),indkpt1(nkpt_rbz) 161 integer,intent(in) :: kg_k(3,npw_k),nattyp(dtset%ntypat),ngfft(18) 162 integer,intent(in) :: pert_atdis(3,natpert),q1grad(3,nq1grad) 163 integer,intent(in) :: q2grad(3,nq2grad),q1q2grad(4,nq1q2grad) 164 real(dp),intent(in) :: cg(2,mpw*dtset%nspinor*dtset%mband*mkmem*nsppol) 165 real(dp),intent(in) :: kpt(3),occ_k(nband_k) 166 real(dp),intent(in) :: ph1d(2,3*(2*dtset%mgfft+1)*dtset%natom) 167 real(dp),intent(out) :: qdrpwf_k(2,natpert,nq2grad,nq1grad) 168 real(dp),intent(out) :: qdrpwf_t1_k(2,natpert,nq2grad,nq1grad) 169 real(dp),intent(out) :: qdrpwf_t2_k(2,natpert,nq2grad,nq1grad) 170 real(dp),intent(out) :: qdrpwf_t3_k(2,natpert,nq2grad,nq1grad) 171 real(dp),intent(out) :: qdrpwf_t4_k(2,natpert,nq2grad,nq1grad) 172 real(dp),intent(out) :: qdrpwf_t5_k(2,natpert,nq2grad,nq1grad) 173 real(dp),intent(in) :: rmet(3,3) 174 real(dp),intent(in) :: vhxc1_atdis(natpert,cplex*nfft) 175 real(dp),intent(in) :: vhxc1_efield(nq2grad,cplex*nfft) 176 real(dp),intent(in) :: xred(3,dtset%natom) 177 real(dp),intent(in) :: ylm_k(npw_k,psps%mpsang*psps%mpsang*psps%useylm) 178 real(dp),intent(in) :: ylmgr_k(npw_k,nylmgr,psps%mpsang*psps%mpsang*psps%useylm*useylmgr) 179 type(wfk_t),intent(inout) :: wfk_t_atdis(natpert),wfk_t_efield(nq2grad) 180 type(wfk_t),intent(inout) :: wfk_t_ddk(nq1grad),wfk_t_dkdk(nq1q2grad) 181 182 !Local variables------------------------------- 183 !scalars 184 integer :: berryopt,iatpert,iband,idir,idirq1,idirq2,ii,ipert 185 integer :: iq1grad,iq2grad,iq1q2grad 186 integer :: jband,nkpg,nkpg1,npw_disk,nq2grad_3d,sij_opt,opt_gvnl1,optlocal,optnl 187 integer :: tim_getgh1c,usevnl 188 integer :: useylmgr1 189 real(dp) :: cprodi,cprodr,doti,dotr,dum_lambda 190 character(len=500) :: msg 191 type(rf_hamiltonian_type) :: rf_hamkq 192 type(pawfgr_type) :: pawfgr 193 194 !arrays 195 real(dp) :: dum_grad_berry(1,1),dum_gs1(1,1),dum_gvnl1(1,1) 196 real(dp),allocatable :: cg1_atdis(:,:),cg1_efield(:,:),cg1_ddk(:,:),cg1_dkdk(:,:) 197 real(dp),allocatable :: cg1_dkdk_ar(:,:,:) 198 real(dp),allocatable :: ci_h1vatdisdag_cj(:,:,:,:),cj_vefield_ci(:,:,:,:) 199 real(dp),allocatable :: cwave0i(:,:),cwave0j(:,:) 200 real(dp),allocatable :: c0_calHatdisdagdQ_c1efield_bks(:,:,:,:,:) 201 real(dp),allocatable :: c0_Hatdisdq_c1efield_bks(:,:,:,:,:) 202 real(dp),allocatable :: c1atdis_c1dkdk_bks(:,:,:,:,:) 203 real(dp),allocatable :: c1atdis_dQVefield_c0_bks(:,:,:,:,:) 204 real(dp),allocatable :: c1atdis_q1gradH0_c1efield_bks(:,:,:,:,:) 205 real(dp),allocatable :: dkinpw(:) 206 real(dp),allocatable :: ffnlk(:,:,:,:),ffnl1(:,:,:,:) 207 real(dp),allocatable :: gh1dqc(:,:),gvloc1dqc(:,:),gvnl1dqc(:,:),gv1c(:,:) 208 real(dp),allocatable :: kinpw1(:),kpg_k(:,:),kpg1_k(:,:),ph3d(:,:,:),ph3d1(:,:,:) 209 real(dp),allocatable :: dum_vlocal(:,:,:,:),vlocal1(:,:,:,:),vlocal1dq(:,:,:,:), dum_vpsp(:) 210 real(dp),allocatable :: vpsp1(:),vpsp1dq(:) 211 real(dp),allocatable :: dum_ylmgr1_k(:,:,:) 212 type(pawcprj_type),allocatable :: dum_cwaveprj(:,:) 213 214 ! ************************************************************************* 215 216 DBG_ENTER("COLL") 217 218 !Not valid for PAW 219 if (psps%usepaw==1) then 220 msg=' This routine cannot be used for PAW (use pawnst3 instead) !' 221 ABI_BUG(msg) 222 end if 223 224 if(dtset%prtvol>2)then 225 write(msg,'(2a,i5,2x,a,3f9.5)')ch10,' Quadrupoles calculation; k pt #',ikpt,'k=',kpt(:) 226 call wrtout(std_out,msg) 227 end if 228 229 !Additional definitions 230 tim_getgh1c=0 231 232 !Additional allocations 233 ABI_MALLOC(cwave0i,(2,npw_k*dtset%nspinor)) 234 ABI_MALLOC(cwave0j,(2,npw_k*dtset%nspinor)) 235 ABI_MALLOC(gv1c,(2,npw_k*dtset%nspinor)) 236 ABI_MALLOC(vlocal1,(cplex*ngfft(4),ngfft(5),ngfft(6),gs_hamkq%nvloc)) 237 ABI_MALLOC(dum_vpsp,(nfft)) 238 ABI_MALLOC(dum_vlocal,(ngfft(4),ngfft(5),ngfft(6),gs_hamkq%nvloc)) 239 ABI_MALLOC(vpsp1,(cplex*nfft)) 240 ABI_MALLOC(dum_cwaveprj,(0,0)) 241 242 243 !-------------------------------------------------------------------------------------- 244 !Electric field 1st order potential matrix element: 245 ! < u_{j,k}^{(0)} | V^{E_{\delta}} | u_{i,k}^{(0)} > 246 !-------------------------------------------------------------------------------------- 247 248 !Specific definitions 249 ipert=dtset%natom+2 250 useylmgr1=0 251 dum_lambda=zero 252 berryopt=0;optlocal=1;optnl=0;usevnl=0;opt_gvnl1=1;sij_opt=0 253 254 !Specific allocations 255 ABI_MALLOC(cj_vefield_ci,(2,nq2grad,nband_k,nband_k)) 256 ABI_MALLOC(dum_ylmgr1_k,(npw_k,3+6*((ipert-dtset%natom)/10),psps%mpsang*psps%mpsang*psps%useylm*useylmgr1)) 257 258 259 !LOOP OVER ELECTRIC FIELD PERTURBATIONS 260 do iq2grad=1,nq2grad 261 ipert=q2grad(1,iq2grad) 262 idir=q2grad(2,iq2grad) 263 264 !Set up local potential vlocal1 with proper dimensioning, from vhxc1_efield 265 vpsp1=vhxc1_efield(iq2grad,:) 266 call rf_transgrid_and_pack(isppol,nspden,psps%usepaw,cplex,nfft,dtset%nfft,dtset%ngfft,& 267 & gs_hamkq%nvloc,pawfgr,mpi_enreg,dum_vpsp,vpsp1,dum_vlocal,vlocal1) 268 269 270 !Initializes rf_hamiltonian (the k-dependent part is prepared in getgh1c_setup) 271 call init_rf_hamiltonian(cplex,gs_hamkq,ipert,rf_hamkq,& 272 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab) 273 call rf_hamkq%load_spin(isppol,vlocal1=vlocal1,with_nonlocal=.false.) 274 275 !Set up the ground-state Hamiltonian, and some parts of the 1st-order Hamiltonian 276 call getgh1c_setup(gs_hamkq,rf_hamkq,dtset,psps,& ! In 277 kpt,kpt,idir,ipert,dtset%natom,rmet,gs_hamkq%gprimd,gs_hamkq%gmet,istwf_k,& ! In 278 npw_k,npw_k,useylmgr1,kg_k,ylm_k,kg_k,ylm_k,dum_ylmgr1_k,& ! In 279 dkinpw,nkpg,nkpg1,kpg_k,kpg1_k,kinpw1,ffnlk,ffnl1,ph3d,ph3d1) ! Out 280 281 !LOOP OVER KET BANDS 282 do iband=1,nband_k 283 284 if(mpi_enreg%proc_distrb(ikpt,iband,isppol) /= mpi_enreg%me_kpt) cycle 285 286 !Read ket ground-state wavefunctions 287 cwave0i(:,:)=cg(:,1+(iband-1)*npw_k*dtset%nspinor+icg:iband*npw_k*dtset%nspinor+icg) 288 289 !Compute <g | V^{E_{\delta}} | u_{i,k}^{(0)} > 290 call getgh1c(berryopt,cwave0i,dum_cwaveprj,gv1c,dum_grad_berry,& 291 & dum_gs1,gs_hamkq,dum_gvnl1,idir,ipert,dum_lambda,mpi_enreg,optlocal,& 292 & optnl,opt_gvnl1,rf_hamkq,sij_opt,tim_getgh1c,usevnl) 293 294 !LOOP OVER BRA BANDS 295 do jband=1,nband_k 296 297 !Read bra ground-state wavefunctions 298 cwave0j(:,:)=cg(:,1+(jband-1)*npw_k*dtset%nspinor+icg:jband*npw_k*dtset%nspinor+icg) 299 300 !Compute cj_vefield_ci= < u_{j,k}^{(0)} | V^{E_{\delta}} | u_{i,k}^{(0)} > 301 call dotprod_g(dotr,doti,istwf_k,npw_k*dtset%nspinor,2,cwave0j,gv1c,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 302 cj_vefield_ci(1,iq2grad,jband,iband)=dotr 303 cj_vefield_ci(2,iq2grad,jband,iband)=doti 304 305 end do !jband 306 307 end do !iband 308 309 !Clean the electric field rf_hamiltonian 310 call rf_hamkq%free() 311 312 !Deallocations 313 ABI_FREE(kpg_k) 314 ABI_FREE(kpg1_k) 315 ABI_FREE(dkinpw) 316 ABI_FREE(kinpw1) 317 ABI_FREE(ffnlk) 318 ABI_FREE(ffnl1) 319 ABI_FREE(ph3d) 320 321 end do !iq2grad 322 !End loop over electric field perturbations 323 324 !Deallocations 325 ABI_FREE(dum_ylmgr1_k) 326 327 !----------------------------------------------------------------------------------------------- 328 ! Atomic displacement 1st order hamiltonian + 1st order potential matrix element: 329 ! < u_{i,k}^{(0)} | (H^{\tau_{\kappa\beta}}+V^{\tau_{\kappa\beta}})^{\dagger} | u_{j,k}^{(0)} > 330 ! calculated as 331 ! (< u_{j,k}^{(0)} | H^{\tau_{\kappa\beta}}+V^{\tau_{\kappa\beta}} | u_{i,k}^{(0)} >)^* 332 !----------------------------------------------------------------------------------------------- 333 334 !Specific definitions 335 ipert=1 336 useylmgr1=0 337 dum_lambda=zero 338 berryopt=0;optlocal=1;optnl=1;usevnl=0;opt_gvnl1=0;sij_opt=0 339 340 !Specific allocations 341 ABI_MALLOC(ci_h1vatdisdag_cj,(2,natpert,nband_k,nband_k)) 342 ABI_MALLOC(dum_ylmgr1_k,(npw_k,3+6*((ipert-dtset%natom)/10), psps%mpsang*psps%mpsang*psps%useylm*useylmgr1)) 343 344 !LOOP OVER ATOMIC DISPLACEMENT PERTURBATIONS 345 do iatpert= 1, natpert 346 ipert=pert_atdis(1,iatpert) 347 idir=pert_atdis(2,iatpert) 348 349 !Get first-order local part of the pseudopotential 350 call dfpt_vlocal(atindx,cplex,gs_hamkq%gmet,gsqcut,idir,ipert,mpi_enreg, & 351 & psps%mqgrid_vl,dtset%natom,& 352 & nattyp,nfft,ngfft,dtset%ntypat,ngfft(1),ngfft(2),ngfft(3), & 353 & ph1d,psps%qgrid_vl,& 354 & dtset%qptn,ucvol,psps%vlspl,vpsp1,xred) 355 356 !Set up local potential vlocal1 with proper dimensioning, from vpsp1 + vhxc1_atdis 357 vpsp1=vpsp1+vhxc1_atdis(iatpert,:) 358 call rf_transgrid_and_pack(isppol,nspden,psps%usepaw,cplex,nfft,dtset%nfft,dtset%ngfft,& 359 & gs_hamkq%nvloc,pawfgr,mpi_enreg,dum_vpsp,vpsp1,dum_vlocal,vlocal1) 360 361 !Initialize rf_hamiltonian (the k-dependent part is prepared in getgh1c_setup) 362 call init_rf_hamiltonian(cplex,gs_hamkq,ipert,rf_hamkq,& 363 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab) 364 call rf_hamkq%load_spin(isppol,vlocal1=vlocal1,with_nonlocal=.true.) 365 366 !Set up the ground-state Hamiltonian, and some parts of the 1st-order Hamiltonian 367 call getgh1c_setup(gs_hamkq,rf_hamkq,dtset,psps,& ! In 368 kpt,kpt,idir,ipert,dtset%natom,rmet,gs_hamkq%gprimd,gs_hamkq%gmet,istwf_k,& ! In 369 npw_k,npw_k,useylmgr1,kg_k,ylm_k,kg_k,ylm_k,dum_ylmgr1_k,& ! In 370 dkinpw,nkpg,nkpg1,kpg_k,kpg1_k,kinpw1,ffnlk,ffnl1,ph3d,ph3d1) ! Out 371 372 !LOOP OVER KET BANDS 373 do iband=1,nband_k 374 375 if(mpi_enreg%proc_distrb(ikpt,iband,isppol) /= mpi_enreg%me_kpt) cycle 376 377 !Read ket ground-state wavefunctions 378 cwave0i(:,:)=cg(:,1+(iband-1)*npw_k*dtset%nspinor+icg:iband*npw_k*dtset%nspinor+icg) 379 380 !Compute < g | H^{\tau_{\kappa\beta}}+V^{\tau_{\kappa\beta}} | u_{i,k}^{(0)} > 381 call getgh1c(berryopt,cwave0i,dum_cwaveprj,gv1c,dum_grad_berry,& 382 & dum_gs1,gs_hamkq,dum_gvnl1,idir,ipert,dum_lambda,mpi_enreg,optlocal,& 383 & optnl,opt_gvnl1,rf_hamkq,sij_opt,tim_getgh1c,usevnl) 384 385 !LOOP OVER BRA BANDS 386 do jband=1,nband_k 387 388 !Read bra ground-state wavefunctions 389 cwave0j(:,:)=cg(:,1+(jband-1)*npw_k*dtset%nspinor+icg:jband*npw_k*dtset%nspinor+icg) 390 391 !Compute = ci_h1vatdisdag_cj 392 !(< u_{j,k}^{(0)} | H^{\tau_{\kappa\beta}}+V^{\tau_{\kappa\beta}} | u_{i,k}^{(0)} >)^* 393 call dotprod_g(dotr,doti,istwf_k,npw_k*dtset%nspinor,2,cwave0j,gv1c,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 394 ci_h1vatdisdag_cj(1,iatpert,jband,iband)=dotr 395 ci_h1vatdisdag_cj(2,iatpert,jband,iband)=-doti 396 397 end do !jband 398 399 end do !iband 400 401 !Clean the atomic displacement rf_hamiltonian 402 call rf_hamkq%free() 403 404 !Deallocations 405 ABI_FREE(kpg_k) 406 ABI_FREE(kpg1_k) 407 ABI_FREE(dkinpw) 408 ABI_FREE(kinpw1) 409 ABI_FREE(ffnlk) 410 ABI_FREE(ffnl1) 411 ABI_FREE(ph3d) 412 413 end do !iatpert 414 !End loop over atomic displacement perturbations 415 416 !Deallocations 417 ABI_FREE(dum_ylmgr1_k) 418 ABI_FREE(cwave0j) 419 ABI_FREE(vpsp1) 420 ABI_FREE(vlocal1) 421 422 !---------------------------------------------------------------------------------------- 423 ! Terms that involve first order response functions 424 !---------------------------------------------------------------------------------------- 425 426 !Allocation of bks (band, k-point and spin) dependent terms 427 ABI_MALLOC(c1atdis_dQVefield_c0_bks,(2,nband_k,natpert,nq2grad,nq1grad)) 428 !TODO:For the moment c1atdis_c1dkdk_bks is allocated for all three directions, i.e., 429 !nq2grad=3. This will have to be modified in the future when ABINIT enables to calculate specific 430 !components of the d2_dkdk 431 nq2grad_3d=3 432 ABI_MALLOC(c1atdis_c1dkdk_bks,(2,nband_k,natpert,nq2grad_3d,nq1grad)) 433 ABI_MALLOC(c0_calHatdisdagdQ_c1efield_bks,(2,nband_k,natpert,nq2grad,nq1grad)) 434 ABI_MALLOC(c1atdis_q1gradH0_c1efield_bks,(2,nband_k,natpert,nq2grad,nq1grad)) 435 ABI_MALLOC(c0_Hatdisdq_c1efield_bks,(2,nband_k,natpert,nq2grad,nq1grad)) 436 c1atdis_dQVefield_c0_bks=zero 437 c0_calHatdisdagdQ_c1efield_bks=zero 438 439 !Allocation of wf1s 440 ABI_MALLOC(cg1_atdis,(2,npw_k*dtset%nspinor)) 441 ABI_MALLOC(cg1_efield,(2,npw_k*dtset%nspinor)) 442 ABI_MALLOC(cg1_ddk,(2,npw_k*dtset%nspinor)) 443 ABI_MALLOC(cg1_dkdk,(2,npw_k*dtset%nspinor)) 444 ABI_MALLOC(cg1_dkdk_ar,(nq1q2grad,2,npw_k*dtset%nspinor)) 445 446 !Check correspondance with the data in wf1 files 447 !Atomic displacements 448 do iatpert=1,natpert 449 !k-point index check 450 ii = wfk_t_atdis(iatpert)%findk(kpt(:)) 451 ABI_CHECK(ii == indkpt1(ikpt), sjoin("ii != indkpt in atomic displacement wf1 file for iatpert:", itoa(iatpert))) 452 !npw check 453 npw_disk = wfk_t_atdis(iatpert)%hdr%npwarr(ii) 454 if (npw_k /= npw_disk) then 455 write(unit=msg,fmt='(a,i5,a,i5,a,a,i5,a,a,i5)')& 456 & 'For ikpt = ',ikpt,', and pertcase= ',pert_atdis(3,iatpert),ch10,& 457 & 'the number of plane waves in the wf1 file is equal to', npw_disk,ch10,& 458 & 'while it should be ',npw_k 459 ABI_BUG(msg) 460 end if 461 end do 462 463 !Electric field 464 do iq2grad=1,nq2grad 465 !k-point index check 466 ii = wfk_t_efield(iq2grad)%findk(kpt(:)) 467 ABI_CHECK(ii == indkpt1(ikpt), sjoin("ii != indkpt1 in electric field wf1 file for iq2grad:", itoa(iq2grad))) 468 !npw check 469 npw_disk = wfk_t_efield(iq2grad)%hdr%npwarr(ii) 470 if (npw_k /= npw_disk) then 471 write(unit=msg,fmt='(a,i5,a,i5,a,a,i5,a,a,i5)')& 472 & 'For ikpt = ',ikpt,', and pertcase= ',q2grad(3,iq2grad),ch10,& 473 & 'the number of plane waves in the wf1 file is equal to', npw_disk,ch10,& 474 & 'while it should be ',npw_k 475 ABI_BUG(msg) 476 end if 477 end do 478 479 !ddk 480 do iq1grad=1,nq1grad 481 !k-point index check 482 ii = wfk_t_ddk(iq1grad)%findk(kpt(:)) 483 ABI_CHECK(ii == indkpt1(ikpt), sjoin("ii != indkpt1 in ddk wf1 file for iq1grad:", itoa(iq1grad))) 484 !npw check 485 npw_disk = wfk_t_ddk(iq1grad)%hdr%npwarr(ii) 486 if (npw_k /= npw_disk) then 487 write(unit=msg,fmt='(a,i5,a,i5,a,a,i5,a,a,i5)')& 488 & 'For ikpt = ',ikpt,', and pertcase= ',q1grad(3,iq1grad),ch10,& 489 & 'the number of plane waves in the wf1 file is equal to', npw_disk,ch10,& 490 & 'while it should be ',npw_k 491 ABI_BUG(msg) 492 end if 493 end do 494 495 !dkdk 496 do iq1q2grad=1,nq1q2grad 497 !k-point index check 498 ii = wfk_t_dkdk(iq1q2grad)%findk(kpt(:)) 499 ABI_CHECK(ii == indkpt1(ikpt), sjoin("ii != indkpt1 in dkdk wf1 file for iq1q2grad:", itoa(iq1q2grad))) 500 !npw check 501 npw_disk = wfk_t_dkdk(iq1q2grad)%hdr%npwarr(ii) 502 if (npw_k /= npw_disk) then 503 write(unit=msg,fmt='(a,i5,a,i5,a,a,i5,a,a,i5)')& 504 & 'For ikpt = ',ikpt,', and pertcase= ',q1q2grad(4,iq1q2grad),ch10,& 505 & 'the number of plane waves in the wf1 file is equal to', npw_disk,ch10,& 506 & 'while it should be ',npw_k 507 ABI_BUG(msg) 508 end if 509 end do 510 511 !-------------------------------------------------------------------------------------- 512 !q1-gradient of gs Hamiltonian: 513 ! < u_{i,k}^{\tau_{\kappa\beta}} | \partial_{gamma} H^{(0)} | u_{i,k}^{E_{\delta}} > 514 !-------------------------------------------------------------------------------------- 515 516 !Specific allocations 517 ABI_MALLOC(vlocal1,(cplex*ngfft(4),ngfft(5),ngfft(6),gs_hamkq%nvloc)) 518 519 !Specific definitions 520 vlocal1=zero 521 useylmgr1=1 522 dum_lambda=zero 523 berryopt=0;optlocal=0;optnl=1;usevnl=0;opt_gvnl1=0;sij_opt=0 524 525 !LOOP OVER Q1-GRADIENT 526 do iq1grad=1,nq1grad 527 ipert=q1grad(1,iq1grad) 528 idir=q1grad(2,iq1grad) 529 530 !Initializes rf_hamiltonian (the k-dependent part is prepared in getgh1c_setup) 531 call init_rf_hamiltonian(cplex,gs_hamkq,ipert,rf_hamkq,& 532 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab) 533 call rf_hamkq%load_spin(isppol,vlocal1=vlocal1,with_nonlocal=.true.) 534 535 !Set up the ground-state Hamiltonian, and some parts of the 1st-order Hamiltonian 536 call getgh1c_setup(gs_hamkq,rf_hamkq,dtset,psps,& ! In 537 kpt,kpt,idir,ipert,dtset%natom,rmet,gs_hamkq%gprimd,gs_hamkq%gmet,istwf_k,& ! In 538 npw_k,npw_k,useylmgr1,kg_k,ylm_k,kg_k,ylm_k,ylmgr_k,& ! In 539 dkinpw,nkpg,nkpg1,kpg_k,kpg1_k,kinpw1,ffnlk,ffnl1,ph3d,ph3d1) ! Out 540 541 !LOOP OVER BANDS 542 do iband=1,nband_k 543 544 if(mpi_enreg%proc_distrb(ikpt,iband,isppol) /= mpi_enreg%me_kpt) cycle 545 546 !LOOP OVER ELECTRIC FIELD PERTURBATION 547 do iq2grad=1,nq2grad 548 549 !Read electric field wf1 550 call wfk_t_efield(iq2grad)%read_bks( iband, indkpt1(ikpt), & 551 & isppol, xmpio_single, cg_bks=cg1_efield) 552 553 !Compute < g |\partial_{gamma} H^{(0)} | u_{i,k}^{E_{\delta}} > 554 call getgh1c(berryopt,cg1_efield,dum_cwaveprj,gv1c,dum_grad_berry,& 555 & dum_gs1,gs_hamkq,dum_gvnl1,idir,ipert,dum_lambda,mpi_enreg,optlocal,& 556 & optnl,opt_gvnl1,rf_hamkq,sij_opt,tim_getgh1c,usevnl) 557 558 !LOOP OVER ATOMIC DISPLACEMENT PERTURBATION 559 do iatpert=1,natpert 560 561 !Read atomic displacement wf1 562 call wfk_t_atdis(iatpert)%read_bks( iband, indkpt1(ikpt), & 563 & isppol, xmpio_single, cg_bks=cg1_atdis) 564 565 !calculate: 566 !< u_{i,k}^{\tau_{\kappa\beta}} | \partial_{gamma} H^{(0)} | u_{i,k}^{E_{\delta}} > 567 call dotprod_g(dotr,doti,istwf_k,npw_k*dtset%nspinor,2,cg1_atdis,gv1c, & 568 & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 569 570 c1atdis_q1gradH0_c1efield_bks(1,iband,iatpert,iq2grad,iq1grad)= dotr 571 c1atdis_q1gradH0_c1efield_bks(2,iband,iatpert,iq2grad,iq1grad)= doti 572 573 end do !iatpert 574 575 end do !iq2grad 576 577 end do !iband 578 579 !Clean the ddk rf_hamiltonian 580 call rf_hamkq%free() 581 582 !Deallocations 583 ABI_FREE(kpg_k) 584 ABI_FREE(kpg1_k) 585 ABI_FREE(dkinpw) 586 ABI_FREE(kinpw1) 587 ABI_FREE(ffnlk) 588 ABI_FREE(ffnl1) 589 ABI_FREE(ph3d) 590 591 end do !iq1grad 592 593 !Deallocations 594 ABI_FREE(gv1c) 595 ABI_FREE(vlocal1) 596 !ABI_FREE(ph3d1) !it is only allocated if kpt and kpq are different. Not the case. 597 598 !-------------------------------------------------------------------------------------- 599 !Other three therms 600 !-------------------------------------------------------------------------------------- 601 602 !LOOP OVER BANDS 603 do iband=1,nband_k 604 605 if(mpi_enreg%proc_distrb(ikpt,iband,isppol) /= mpi_enreg%me_kpt) cycle 606 607 !LOOP OVER ATOMIC DISPLACEMENT PERTURBATION 608 do iatpert=1,natpert 609 610 !Read atomic displacement wf1 611 call wfk_t_atdis(iatpert)%read_bks( iband, indkpt1(ikpt), & 612 & isppol, xmpio_single, cg_bks=cg1_atdis) 613 614 !LOOP OVER q1-GRADIENT 615 do iq1grad=1,nq1grad 616 617 !LOOP OVER BANDS 618 do jband=1,nband_k 619 620 !Read ddk wf1 621 call wfk_t_ddk(iq1grad)%read_bks( jband, indkpt1(ikpt), & 622 & isppol, xmpio_single, cg_bks=cg1_ddk) 623 624 !Calculate: < u_{i,k}^{\tau_{\kappa\beta}} | u_{j,k}^{k_{\gamma}} > 625 call dotprod_g(dotr,doti,istwf_k,npw_k*dtset%nspinor,2,cg1_atdis,cg1_ddk, & 626 & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 627 628 !LOOP OVER ELECTRIC FIELD PERTURBATION 629 do iq2grad=1,nq2grad 630 631 !Calculate: -\sum_{j} < u_{i,k}^{\tau_{\kappa\beta}} | u_{j,k}^{k_{\gamma}} > * 632 ! < u_{j,k}^{(0)} | V^{E_{\delta}} | u_{i,k}^{(0)} > 633 cprodr= dotr*cj_vefield_ci(1,iq2grad,jband,iband) - & 634 & doti*cj_vefield_ci(2,iq2grad,jband,iband) 635 cprodi= dotr*cj_vefield_ci(2,iq2grad,jband,iband) + & 636 & doti*cj_vefield_ci(1,iq2grad,jband,iband) 637 638 c1atdis_dQVefield_c0_bks(1,iband,iatpert,iq2grad,iq1grad)= & 639 & c1atdis_dQVefield_c0_bks(1,iband,iatpert,iq2grad,iq1grad)-cprodr 640 c1atdis_dQVefield_c0_bks(2,iband,iatpert,iq2grad,iq1grad)= & 641 & c1atdis_dQVefield_c0_bks(2,iband,iatpert,iq2grad,iq1grad)-cprodi 642 643 644 end do !iq2grad 645 646 end do !jband 647 648 end do !iq1grad 649 650 !LOOP OVER q1q2 SECOND ORDER GRADIENT 651 do iq1q2grad=1,nq1q2grad 652 653 !Read dkdk wf1 654 if (iq1q2grad<=6) then 655 call wfk_t_dkdk(iq1q2grad)%read_bks( iband, indkpt1(ikpt), & 656 & isppol, xmpio_single, cg_bks=cg1_dkdk) 657 cg1_dkdk_ar(iq1q2grad,:,:)=cg1_dkdk 658 else 659 cg1_dkdk=cg1_dkdk_ar(iq1q2grad-3,:,:) 660 end if 661 662 !Calculate: i/2 < u_{i,k}^{\tau_{\kappa\beta}} | u_{i,k}^{k_{\gamma},k_{\delta}} > 663 call dotprod_g(dotr,doti,istwf_k,npw_k*dtset%nspinor,2,cg1_atdis,cg1_dkdk, & 664 & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 665 666 idirq1=q1q2grad(2,iq1q2grad) 667 idirq2=q1q2grad(3,iq1q2grad) 668 c1atdis_c1dkdk_bks(1,iband,iatpert,idirq2,idirq1)= -half*doti 669 c1atdis_c1dkdk_bks(2,iband,iatpert,idirq2,idirq1)= half*dotr 670 671 end do !iq1q2grad 672 673 end do !iatpert 674 675 !LOOP OVER ELECTRIC FIELD PERTURBATION 676 do iq2grad=1,nq2grad 677 678 !Read electric field wf1 679 call wfk_t_efield(iq2grad)%read_bks( iband, indkpt1(ikpt), & 680 & isppol, xmpio_single, cg_bks=cg1_efield) 681 682 !LOOP OVER q1-GRADIENT 683 do iq1grad=1,nq1grad 684 685 !LOOP OVER BANDS 686 do jband=1,nband_k 687 688 !Read ddk wf1 689 call wfk_t_ddk(iq1grad)%read_bks( jband, indkpt1(ikpt), & 690 & isppol, xmpio_single, cg_bks=cg1_ddk) 691 692 !Calculate: < u_{j,k}^{k_{\gamma}} | u_{i,k}^{E_{\delta}} > 693 call dotprod_g(dotr,doti,istwf_k,npw_k*dtset%nspinor,2,cg1_ddk,cg1_efield, & 694 & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 695 696 !LOOP OVER ATOMIC DISPLACEMENT PERTURBATION 697 do iatpert=1,natpert 698 699 !Calculate: -\sum_{j} 700 ! < u_{i,k}^{(0)} | (H^{\tau_{\kappa\beta}}+V^{\tau_{\kappa\beta}})^{\dagger} | u_{j,k}^{(0)} > 701 ! < u_{j,k}^{k_{\gamma}} | u_{i,k}^{E_{\delta}} > 702 cprodr=dotr*ci_h1vatdisdag_cj(1,iatpert,jband,iband) - & 703 & doti*ci_h1vatdisdag_cj(2,iatpert,jband,iband) 704 cprodi=dotr*ci_h1vatdisdag_cj(2,iatpert,jband,iband) + & 705 & doti*ci_h1vatdisdag_cj(1,iatpert,jband,iband) 706 707 c0_calHatdisdagdQ_c1efield_bks(1,iband,iatpert,iq2grad,iq1grad)= & 708 & c0_calHatdisdagdQ_c1efield_bks(1,iband,iatpert,iq2grad,iq1grad)-cprodr 709 c0_calHatdisdagdQ_c1efield_bks(2,iband,iatpert,iq2grad,iq1grad)= & 710 & c0_calHatdisdagdQ_c1efield_bks(2,iband,iatpert,iq2grad,iq1grad)-cprodi 711 712 end do !iatpert 713 714 end do !jband 715 716 end do !iq1grad 717 718 end do !iq2grad 719 720 end do !iband 721 722 !-------------------------------------------------------------------------------------- 723 ! q1-gradient of atomic displacement 1st-order Hamiltonian: 724 ! < u_{i,k}^{(0)} | (H^{\tau_{\kappa\beta}}_{\gamma})^{dagger} | u_{i,k}^{E_{\delta}} > 725 ! calculated as: 726 ! (<u_{i,k}^{E_{\delta} | H^{\tau_{\kappa\beta}}_{\gamma} | u_{i,k}^{(0)} >)* 727 !-------------------------------------------------------------------------------------- 728 729 !Specific allocations 730 ABI_MALLOC(vpsp1dq,(2*nfft)) 731 ABI_MALLOC(vlocal1dq,(2*ngfft(4),ngfft(5),ngfft(6),gs_hamkq%nvloc)) 732 ABI_MALLOC(gh1dqc,(2,npw_k*dtset%nspinor)) 733 ABI_MALLOC(gvloc1dqc,(2,npw_k*dtset%nspinor)) 734 ABI_MALLOC(gvnl1dqc,(2,npw_k*dtset%nspinor)) 735 736 !Specific definitions 737 useylmgr1=1;optlocal=1;optnl=1 738 739 !LOOP OVER ATOMIC DISPLACEMENT PERTURBATIONS 740 do iatpert= 1, natpert 741 ipert=pert_atdis(1,iatpert) 742 idir=pert_atdis(2,iatpert) 743 744 !LOOP OVER Q1-GRADIENT 745 do iq1grad=1,nq1grad 746 747 !Get q-gradient of first-order local part of the pseudopotential 748 call dfpt_vlocaldq(atindx,2,gs_hamkq%gmet,gsqcut,idir,ipert,mpi_enreg, & 749 & psps%mqgrid_vl,dtset%natom,& 750 & nattyp,nfft,ngfft,dtset%ntypat,ngfft(1),ngfft(2),ngfft(3), & 751 & ph1d,q1grad(2,iq1grad),psps%qgrid_vl,& 752 & dtset%qptn,ucvol,psps%vlspl,vpsp1dq) 753 754 755 !Set up q-gradient of local potential vlocal1dq with proper dimensioning 756 call rf_transgrid_and_pack(isppol,nspden,psps%usepaw,2,nfft,dtset%nfft,dtset%ngfft,& 757 & gs_hamkq%nvloc,pawfgr,mpi_enreg,dum_vpsp,vpsp1dq,dum_vlocal,vlocal1dq) 758 759 !Initialize rf_hamiltonian (the k-dependent part is prepared in getgh1c_setup) 760 call init_rf_hamiltonian(2,gs_hamkq,ipert,rf_hamkq,& 761 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab) 762 call rf_hamkq%load_spin(isppol,vlocal1=vlocal1dq,with_nonlocal=.true.) 763 764 !Set up the ground-state Hamiltonian, and some parts of the 1st-order Hamiltonian 765 call getgh1dqc_setup(gs_hamkq,rf_hamkq,dtset,psps,kpt,kpt,idir,ipert,q1grad(2,iq1grad), & 766 & dtset%natom,rmet,gs_hamkq%gprimd,gs_hamkq%gmet,istwf_k,npw_k,npw_k,nylmgr,useylmgr1,kg_k, & 767 & ylm_k,kg_k,ylm_k,ylmgr_k,nkpg,nkpg1,kpg_k,kpg1_k,dkinpw,kinpw1,ffnlk,ffnl1,ph3d,ph3d1) 768 769 !LOOP OVER BANDS 770 do iband=1,nband_k 771 772 if(mpi_enreg%proc_distrb(ikpt,iband,isppol) /= mpi_enreg%me_kpt) cycle 773 774 !Read ket ground-state wavefunctions 775 cwave0i(:,:)=cg(:,1+(iband-1)*npw_k*dtset%nspinor+icg:iband*npw_k*dtset%nspinor+icg) 776 777 !Compute < g |H^{\tau_{\kappa\beta}}_{\gamma} | u_{i,k}^{(0)} > 778 call getgh1dqc(cwave0i,dum_cwaveprj,gh1dqc,gvloc1dqc,gvnl1dqc,gs_hamkq, & 779 & idir,ipert,mpi_enreg,optlocal,optnl,q1grad(2,iq1grad),rf_hamkq) 780 781 !LOOP OVER ELECTRIC FIELD PERTURBATION 782 do iq2grad=1,nq2grad 783 784 !Read electric field wf1 785 call wfk_t_efield(iq2grad)%read_bks( iband, indkpt1(ikpt), & 786 & isppol, xmpio_single, cg_bks=cg1_efield) 787 788 !Calculate: 789 !(<u_{i,k}^{E_{\delta} | H^{\tau_{\kappa\beta}}_{\gamma} | u_{i,k}^{(0)} >)* 790 call dotprod_g(dotr,doti,istwf_k,npw_k*dtset%nspinor,2,cg1_efield,gh1dqc, & 791 & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 792 793 !Apply the -i factor that has been factorized in the 794 !H^{\tau_{\kappa\beta}}_{\gamma} terms. (And consider the complex 795 !conjugate too) 796 c0_Hatdisdq_c1efield_bks(1,iband,iatpert,iq2grad,iq1grad)= doti 797 c0_Hatdisdq_c1efield_bks(2,iband,iatpert,iq2grad,iq1grad)= dotr 798 799 end do !iq2grad 800 801 end do !iband 802 803 !Clean the rf_hamiltonian 804 call rf_hamkq%free() 805 806 !Deallocations 807 ABI_FREE(kpg_k) 808 ABI_FREE(kpg1_k) 809 ABI_FREE(dkinpw) 810 ABI_FREE(kinpw1) 811 ABI_FREE(ffnlk) 812 ABI_FREE(ffnl1) 813 ABI_FREE(ph3d) 814 815 end do !iq1grad 816 817 end do !iatpert 818 819 !Deallocations 820 ABI_FREE(dum_cwaveprj) 821 ABI_FREE(gh1dqc) 822 ABI_FREE(gvloc1dqc) 823 ABI_FREE(gvnl1dqc) 824 ABI_FREE(vpsp1dq) 825 ABI_FREE(vlocal1dq) 826 ABI_FREE(dum_vpsp) 827 ABI_FREE(dum_vlocal) 828 829 830 !-------------------------------------------------------------------------------------- 831 ! Acumulates all the wf dependent terms of the quadrupole tensor 832 !-------------------------------------------------------------------------------------- 833 qdrpwf_k=zero 834 qdrpwf_t1_k=zero 835 qdrpwf_t2_k=zero 836 qdrpwf_t3_k=zero 837 qdrpwf_t4_k=zero 838 qdrpwf_t5_k=zero 839 do iq1grad=1,nq1grad 840 idirq1=q1grad(2,iq1grad) 841 do iq2grad=1,nq2grad 842 idirq2=q2grad(2,iq2grad) 843 do iatpert=1,natpert 844 do iband=1,nband_k 845 846 if(mpi_enreg%proc_distrb(ikpt,iband,isppol) /= mpi_enreg%me_kpt) cycle 847 848 !All terms toghether 849 qdrpwf_k(1,iatpert,iq2grad,iq1grad)=qdrpwf_k(1,iatpert,iq2grad,iq1grad) + & 850 & wtk_k * occ_k(iband) * & 851 & ( c1atdis_q1gradH0_c1efield_bks(1,iband,iatpert,iq2grad,iq1grad) + & !T1 852 & c1atdis_dQVefield_c0_bks(1,iband,iatpert,iq2grad,iq1grad) + & !T2 853 & c0_calHatdisdagdQ_c1efield_bks(1,iband,iatpert,iq2grad,iq1grad) + & !T3 854 & c0_Hatdisdq_c1efield_bks(1,iband,iatpert,iq2grad,iq1grad) + & !T4 855 & c1atdis_c1dkdk_bks(1,iband,iatpert,idirq2,idirq1) ) !T5 856 857 qdrpwf_k(2,iatpert,iq2grad,iq1grad)=qdrpwf_k(2,iatpert,iq2grad,iq1grad) + & 858 & wtk_k * occ_k(iband) * & 859 & ( c1atdis_q1gradH0_c1efield_bks(2,iband,iatpert,iq2grad,iq1grad) + & !T1 860 & c1atdis_dQVefield_c0_bks(2,iband,iatpert,iq2grad,iq1grad) + & !T2 861 & c0_calHatdisdagdQ_c1efield_bks(2,iband,iatpert,iq2grad,iq1grad) + & !T3 862 & c0_Hatdisdq_c1efield_bks(2,iband,iatpert,iq2grad,iq1grad) + & !T4 863 & c1atdis_c1dkdk_bks(2,iband,iatpert,idirq2,idirq1) ) !T5 864 865 !Separate them 866 !T1 867 qdrpwf_t1_k(1,iatpert,iq2grad,iq1grad)=qdrpwf_t1_k(1,iatpert,iq2grad,iq1grad) + & 868 & wtk_k * occ_k(iband) * & 869 & c1atdis_q1gradH0_c1efield_bks(1,iband,iatpert,iq2grad,iq1grad) 870 871 qdrpwf_t1_k(2,iatpert,iq2grad,iq1grad)=qdrpwf_t1_k(2,iatpert,iq2grad,iq1grad) + & 872 & wtk_k * occ_k(iband) * & 873 & c1atdis_q1gradH0_c1efield_bks(2,iband,iatpert,iq2grad,iq1grad) 874 875 !T2 876 qdrpwf_t2_k(1,iatpert,iq2grad,iq1grad)=qdrpwf_t2_k(1,iatpert,iq2grad,iq1grad) + & 877 & wtk_k * occ_k(iband) * & 878 & c1atdis_dQVefield_c0_bks(1,iband,iatpert,iq2grad,iq1grad) 879 880 qdrpwf_t2_k(2,iatpert,iq2grad,iq1grad)=qdrpwf_t2_k(2,iatpert,iq2grad,iq1grad) + & 881 & wtk_k * occ_k(iband) * & 882 & c1atdis_dQVefield_c0_bks(2,iband,iatpert,iq2grad,iq1grad) 883 884 !T3 885 qdrpwf_t3_k(1,iatpert,iq2grad,iq1grad)=qdrpwf_t3_k(1,iatpert,iq2grad,iq1grad) + & 886 & wtk_k * occ_k(iband) * & 887 & c0_calHatdisdagdQ_c1efield_bks(1,iband,iatpert,iq2grad,iq1grad) 888 889 qdrpwf_t3_k(2,iatpert,iq2grad,iq1grad)=qdrpwf_t3_k(2,iatpert,iq2grad,iq1grad) + & 890 & wtk_k * occ_k(iband) * & 891 & c0_calHatdisdagdQ_c1efield_bks(2,iband,iatpert,iq2grad,iq1grad) 892 893 !T4 894 qdrpwf_t4_k(1,iatpert,iq2grad,iq1grad)=qdrpwf_t4_k(1,iatpert,iq2grad,iq1grad) + & 895 & wtk_k * occ_k(iband) * & 896 & c0_Hatdisdq_c1efield_bks(1,iband,iatpert,iq2grad,iq1grad) 897 898 qdrpwf_t4_k(2,iatpert,iq2grad,iq1grad)=qdrpwf_t4_k(2,iatpert,iq2grad,iq1grad) + & 899 & wtk_k * occ_k(iband) * & 900 & c0_Hatdisdq_c1efield_bks(2,iband,iatpert,iq2grad,iq1grad) 901 902 !T5 903 qdrpwf_t5_k(1,iatpert,iq2grad,iq1grad)=qdrpwf_t5_k(1,iatpert,iq2grad,iq1grad) + & 904 & wtk_k * occ_k(iband) * c1atdis_c1dkdk_bks(1,iband,iatpert,idirq2,idirq1) 905 906 qdrpwf_t5_k(2,iatpert,iq2grad,iq1grad)=qdrpwf_t5_k(2,iatpert,iq2grad,iq1grad) + & 907 & wtk_k * occ_k(iband) * c1atdis_c1dkdk_bks(2,iband,iatpert,idirq2,idirq1) 908 909 end do 910 end do 911 end do 912 end do 913 914 915 !Deallocations 916 ABI_FREE(cj_vefield_ci) 917 ABI_FREE(ci_h1vatdisdag_cj) 918 ABI_FREE(c1atdis_dQVefield_c0_bks) 919 ABI_FREE(c1atdis_c1dkdk_bks) 920 ABI_FREE(c0_calHatdisdagdQ_c1efield_bks) 921 ABI_FREE(c1atdis_q1gradH0_c1efield_bks) 922 ABI_FREE(c0_Hatdisdq_c1efield_bks) 923 ABI_FREE(cwave0i) 924 ABI_FREE(cg1_atdis) 925 ABI_FREE(cg1_efield) 926 ABI_FREE(cg1_ddk) 927 ABI_FREE(cg1_dkdk) 928 ABI_FREE(cg1_dkdk_ar) 929 930 DBG_EXIT("COLL") 931 932 end subroutine dfpt_qdrpwf
ABINIT/m_dfpt_lwwf [ Modules ]
NAME
m_dfpt_lwwf
FUNCTION
First-order response function contributions to the quadrupole and flexoelectric tensors
COPYRIGHT
Copyright (C) 2018-2022 ABINIT group (MR,MS) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
NOTES
SOURCE
18 #if defined HAVE_CONFIG_H 19 #include "config.h" 20 #endif 21 22 #include "abi_common.h" 23 24 module m_dfpt_lwwf 25 26 use defs_basis 27 use defs_abitypes 28 use defs_datatypes 29 use m_dtset 30 use m_errors 31 use m_profiling_abi 32 use m_hamiltonian 33 use m_cgtools 34 use m_pawcprj 35 use m_pawfgr 36 use m_wfk 37 use m_xmpi 38 use m_getgh1c 39 use m_mklocl 40 41 use m_fstrings, only : itoa, sjoin 42 use m_io_tools, only : file_exists 43 use m_time, only : cwtime 44 use m_kg, only : mkkpg 45 46 implicit none 47 48 private