TABLE OF CONTENTS


ABINIT/dfpt_ciflexowf [ Functions ]

[ Top ] [ 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 ]

[ Top ] [ 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 ]

[ Top ] [ 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 ]

[ Top ] [ 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 ]

[ Top ] [ 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 ]

[ Top ] [ 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