TABLE OF CONTENTS


ABINIT/corrmetalwf1 [ Functions ]

[ Top ] [ Functions ]

NAME

 corrmetalwf1

FUNCTION

 Response function calculation only:
 Correct 1st-order wave-function, taking into account "metallic" occupations.
 1st-order WF orthogonal to C_n,k+q, restore the "active space" content of the first-order WF.

INPUTS

  cg(2,mcgq)=planewave coefficients of wavefunctions at k+q
  cprjq(natom,mcprjq)= wave functions at k+q projected with non-local projectors
  cwavef(2,npw1*nspinor)= 1st-order wave-function before correction
  cwaveprj(natom,nspinor)= 1st-order wave-function before correction
                           projected on NL projectors (PAW)
  eig1(2*nband**2)=first-order eigenvalues (hartree)
  fermie1=derivative of fermi energy wrt (strain) perturbation
  ghc(2,npw1*nspinor)=<G|H0-eig0_k.I|C1 band,k> (NCPP) or <G|H0-eig0_k.S0|C1 band,k> (PAW)
                      (C1 before correction)
  iband=index of current band
  ibgq=shift to be applied on the location of data in the array cprjq
  icgq=shift to be applied on the location of data in the array cgq
  istwf_k=option parameter that describes the storage of wfs
  mcgq=second dimension of the cgq array
  mcprjq=second dimension of the cprjq array
  mpi_enreg=information about MPI parallelization
  natom=number of atoms in cell
  nband=number of bands
  npw1=number of plane waves at this k+q point
  nspinor=number of spinorial components of the wavefunctions
  occ(nband)=occupation number for each band for each k.
  rocceig(nband,nband)= (occ_kq(m)-occ_k(n))/(eig0_kq(m)-eig0_k(n)),
    if this ratio has been attributed to the band n (second argument), zero otherwise
  timcount=index used to accumulate timing (0 from dfpt_vtowfk, 1 from dfpt_nstwf)
  usepaw=flag for PAW

OUTPUT

  cwave1(2,npw1*nspinor)= 1st-order wave-function after correction
  cwaveprj1(natom,nspinor)= 1st-order wave-function after correction
                            projected on NL projectors (PAW)
  edocc(nband)=correction to 2nd-order total energy coming from changes of occupations
  wf_corrected=flag put to 1 if input cwave1 is effectively different from output cwavef

NOTES

  Was part of dfpt_vtowfk before.

PARENTS

      dfpt_vtowfk

CHILDREN

      cg_zcopy,dotprod_g,pawcprj_copy,pawcprj_zaxpby,timab

SOURCE

834 subroutine corrmetalwf1(cgq,cprjq,cwavef,cwave1,cwaveprj,cwaveprj1,edocc,eig1,fermie1,ghc,iband, &
835 &          ibgq,icgq,istwf_k,mcgq,mcprjq,mpi_enreg,natom,nband,npw1,nspinor,occ,rocceig,timcount,&
836 &          usepaw,wf_corrected)
837 
838 
839 !This section has been created automatically by the script Abilint (TD).
840 !Do not modify the following lines by hand.
841 #undef ABI_FUNC
842 #define ABI_FUNC 'corrmetalwf1'
843 !End of the abilint section
844 
845  implicit none
846 
847 !Arguments ------------------------------------
848 !scalars
849  integer,intent(in) :: iband,ibgq,icgq,istwf_k,mcgq,mcprjq,natom,nband,npw1,nspinor,timcount,usepaw
850  integer,intent(out) :: wf_corrected
851  real(dp),intent(in) :: fermie1
852  type(MPI_type),intent(in) :: mpi_enreg
853 !arrays
854  real(dp),intent(in) :: cgq(2,mcgq),cwavef(2,npw1*nspinor)
855  real(dp),intent(in) :: eig1(2*nband**2),ghc(2,npw1*nspinor),occ(nband),rocceig(nband,nband)
856  real(dp),intent(out) :: cwave1(2,npw1*nspinor),edocc(nband)
857  type(pawcprj_type),intent(in) :: cprjq(natom,mcprjq),cwaveprj(natom,nspinor*usepaw)
858  type(pawcprj_type),intent(inout) :: cwaveprj1(natom,nspinor*usepaw) !vz_i
859 
860 !Local variables-------------------------------
861 !scalars
862  integer :: ibandkq,index_cgq,index_cprjq,index_eig1,ii
863  real(dp) :: facti,factr,invocc
864 !arrays
865  real(dp) :: tsec(2)
866  real(dp),allocatable :: cwcorr(:,:)
867 
868 ! *********************************************************************
869 
870  DBG_ENTER("COLL")
871 
872  call timab(214+timcount,1,tsec)
873 
874 !At this stage, the 1st order function cwavef is orthogonal to cgq (unlike when it is input to dfpt_cgwf).
875 !Here, restore the "active space" content of the 1st-order wavefunction, to give cwave1 .
876 
877 !First copy input WF into output WF
878  wf_corrected=0
879  call cg_zcopy(npw1*nspinor,cwavef,cwave1)
880 
881  if (usepaw==1) then
882    call pawcprj_copy(cwaveprj,cwaveprj1)
883  end if
884 
885 !Correct WF only for occupied states
886  if (abs(occ(iband)) > tol8) then
887    invocc=one/occ(iband)
888 
889    edocc(iband)=zero
890 
891 !  Loop over WF at k+q subspace
892    do ibandkq=1,nband
893 
894 !    Select bands with variable occupation
895      if (abs(rocceig(ibandkq,iband))>tol8) then
896 
897        wf_corrected=1
898 
899        index_eig1=2*ibandkq-1+(iband-1)*2*nband
900        index_cgq=npw1*nspinor*(ibandkq-1)+icgq
901 
902        if(ibandkq==iband) then
903          factr=rocceig(ibandkq,iband)*invocc*(eig1(index_eig1)-fermie1)
904        else
905          factr=rocceig(ibandkq,iband)*invocc*eig1(index_eig1)
906        end if
907        facti= rocceig(ibandkq,iband)*invocc*eig1(index_eig1+1)
908 
909 !      Apply correction to 1st-order WF
910 !$OMP PARALLEL DO PRIVATE(ii) SHARED(cgq,cwave1,facti,factr,index_cgq,npw1,nspinor)
911        do ii=1,npw1*nspinor
912          cwave1(1,ii)=cwave1(1,ii)+(factr*cgq(1,ii+index_cgq)-facti*cgq(2,ii+index_cgq))
913          cwave1(2,ii)=cwave1(2,ii)+(facti*cgq(1,ii+index_cgq)+factr*cgq(2,ii+index_cgq))
914        end do
915 
916 !      In the PAW case, also apply correction to projected WF
917        if (usepaw==1) then
918          index_cprjq=nspinor*(ibandkq-1)+ibgq
919          call pawcprj_zaxpby((/factr,facti/),(/one,zero/),cprjq(:,index_cprjq+1:index_cprjq+nspinor),cwaveprj1)
920        end if
921 
922 !      The factor of two is needed because we compute the 2DTE, and not E(2)
923        edocc(iband)=edocc(iband)-two*(factr*eig1(index_eig1)+facti*eig1(index_eig1+1))
924 
925      end if ! Variable occupations
926    end do ! Loop over k+q subspace
927  end if ! occupied states
928 
929 !In the PAW case, compute <Psi^(1)_ortho|H-Eig0_k.S|Psi^(1)_parallel> contribution to 2DTE
930  if (usepaw==1.and.wf_corrected==1) then
931    ABI_ALLOCATE(cwcorr,(2,npw1*nspinor))
932 !$OMP WORKSHARE
933    cwcorr(:,:)=cwave1(:,:)-cwavef(:,:)
934 !$OMP END WORKSHARE
935    call dotprod_g(factr,facti,istwf_k,npw1*nspinor,1,cwcorr,ghc,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
936    edocc(iband)=edocc(iband)+four*factr
937    ABI_DEALLOCATE(cwcorr)
938  end if
939 
940  call timab(214+timcount,2,tsec)
941 
942  DBG_EXIT("COLL")
943 
944 end subroutine corrmetalwf1

ABINIT/dfpt_vtowfk [ Functions ]

[ Top ] [ Functions ]

NAME

 dfpt_vtowfk

FUNCTION

 This routine compute the partial density at a given k-point,
 for a given spin-polarization, from a fixed potential (vlocal1).

INPUTS

  cg(2,mpw*nspinor*mband*mkmem*nsppol)=planewave coefficients of wavefunctions
  cgq(2,mcgq)=array for planewave coefficients of wavefunctions.
  cg1(2,mpw1*nspinor*mband*mk1mem*nsppol)=pw coefficients of RF wavefunctions at k,q.
  cplex=1 if rhoaug1 is real, 2 if rhoaug1 is complex
  cprj(natom,nspinor*mband*mkmem*nsppol*usecprj)= wave functions at k
              projected with non-local projectors: cprj=<p_i|Cnk>
  cprjq(natom,mcprjq)= wave functions at k+q projected with non-local projectors: cprjq=<p_i|Cnk+q>
  dim_eig2rf = dimension for the second order eigenvalues
  dtfil <type(datafiles_type)>=variables related to files
  dtset <type(dataset_type)>=all input variables for this dataset
  eig0_k(nband_k)=GS eigenvalues at k (hartree)
  eig0_kq(nband_k)=GS eigenvalues at k+Q (hartree)
  fermie1=derivative of fermi energy wrt (strain) perturbation
  grad_berry(2,mpw1,dtefield%mband_occ) = the gradient of the Berry phase term
  gs_hamkq <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k+q
  ibg=shift to be applied on the location of data in the array cprj
  ibgq=shift to be applied on the location of data in the array cprjq
  ibg1=shift to be applied on the location of data in the array cprj1
  icg=shift to be applied on the location of data in the array cg
  icgq=shift to be applied on the location of data in the array cgq
  icg1=shift to be applied on the location of data in the array cg1
  idir=direction of the current perturbation
  ikpt=number of the k-point
  ipert=type of the perturbation
  isppol=1 index of current spin component
  mband=maximum number of bands
  mcgq=second dimension of the cgq array
  mcprjq=second dimension of the cprjq array
  mkmem =number of k points trated by this node (GS data).
  mk1mem =number of k points treated by this node (RF data)
  mpi_enreg=information about MPI parallelization
  mpw=maximum dimensioned size of npw or wfs at k
  mpw1=maximum dimensioned size of npw for wfs at k+q (also for 1-order wfs).
  natom=number of atoms in cell.
  nband_k=number of bands at this k point for that spin polarization
  ncpgr=number of gradients stored in cprj array (cprj=<p_i|Cnk>)
  nnsclo_now=number of non-self-consistent loops for the current vtrial
    (often 1 for SCF calculation, =nstep for non-SCF calculations)
  npw_k=number of plane waves at this k point
  npw1_k=number of plane waves at this k+q point
  nspinor=number of spinorial components of the wavefunctions
  nsppol=1 for unpolarized, 2 for spin-polarized
  n4,n5,n6 used for dimensioning real space arrays
  occ_k(nband_k)=occupation number for each band (usually 2) for each k.
  prtvol=control print volume and debugging output
  psps <type(pseudopotential_type)>=variables related to pseudopotentials
  rf_hamkq <type(rf_hamiltonian_type)>=all data for the 1st-order Hamiltonian at k,q
  rf_hamk_dir2 <type(rf_hamiltonian_type)>= (used only when ipert=natom+11, so q=0)
    same as rf_hamkq, but the direction of the perturbation is different
  rhoaug1(cplex*n4,n5,n6,nspden)= density in electrons/bohr**3,
   on the augmented fft grid. (cumulative, so input as well as output)
  rocceig(nband_k,nband_k)= (occ_kq(m)-occ_k(n))/(eig0_kq(m)-eig0_k(n)),
    if this ratio has been attributed to the band n (second argument), zero otherwise
  ddk<wfk_t>=struct info for DDK file.
  wtk_k=weight assigned to the k point.

OUTPUT

  cg1(2,mpw1*nspinor*mband*mk1mem*nsppol)=pw coefficients of RF
    wavefunctions at k,q. They are orthogonalized to the occupied states.
  cg1_active(2,mpw1*nspinor*mband*mk1mem*nsppol*dim_eig2rf)=pw coefficients of RF
    wavefunctions at k,q. They are orthogonalized to the active. Only needed for ieigrf/=0
  edocc_k(nband_k)=correction to 2nd-order total energy coming
      from changes of occupation
  eeig0_k(nband_k)=zero-order eigenvalues contribution to 2nd-order total
      energy from all bands at this k point.
  eig1_k(2*nband_k**2)=first-order eigenvalues (hartree)
  ek0_k(nband_k)=0-order kinetic energy contribution to 2nd-order total
      energy from all bands at this k point.
  ek1_k(nband_k)=1st-order kinetic energy contribution to 2nd-order total
      energy from all bands at this k point.
  eloc0_k(nband_k)=zero-order local contribution to 2nd-order total energy
      from all bands at this k point.
  enl0_k(nband_k)=zero-order non-local contribution to 2nd-order total energy
      from all bands at this k point.
  enl1_k(nband_k)=first-order non-local contribution to 2nd-order total energy
      from all bands at this k point.
  gh1c_set(2,mpw1*nspinor*mband*mk1mem*nsppol*dim_eig2rf)= set of <G|H^{(1)}|nK>
  gh0c1_set(2,mpw1*nspinor*mband*mk1mem*nsppol*dim_eig2rf)= set of <G|H^{(0)}k+q-eig^{(0)}nk|\Psi^{(1)}kq>
      The wavefunction is orthogonal to the active space (for metals). It is not coherent with cg1.
  resid_k(nband_k)=residuals for each band over all k points,
  rhoaug1(cplex*n4,n5,n6,nspden)= density in electrons/bohr**3,
   on the augmented fft grid. (cumulative, so input as well as output).
  ==== if (gs_hamkq%usepaw==1) ====
    cprj1(natom,nspinor*mband*mk1mem*nsppol*usecprj)=
              1st-order wave functions at k,q projected with non-local projectors:
                       cprj1=<p_i|C1nk,q> where p_i is a non-local projector
    pawrhoij1(natom) <type(pawrhoij_type)>= 1st-order paw rhoij occupancies and related data
                                            (cumulative, so input as well as output)

PARENTS

      dfpt_vtorho

CHILDREN

      cg_zcopy,corrmetalwf1,dfpt_accrho,dfpt_cgwf,dotprod_g,getgsc
      matrixelmt_g,meanvalue_g,pawcprj_alloc,pawcprj_copy,pawcprj_free
      pawcprj_get,pawcprj_put,rf2_destroy,rf2_init,sqnorm_g,status,timab
      wfk_read_bks,wrtout

SOURCE

180 subroutine dfpt_vtowfk(cg,cgq,cg1,cg1_active,cplex,cprj,cprjq,cprj1,&
181 & dim_eig2rf,dtfil,dtset,&
182 & edocc_k,eeig0_k,eig0_k,eig0_kq,eig1_k,&
183 & ek0_k,ek1_k,eloc0_k,enl0_k,enl1_k,&
184 & fermie1,ffnl1,ffnl1_test,gh0c1_set,gh1c_set,grad_berry,gs_hamkq,&
185 & ibg,ibgq,ibg1,icg,icgq,icg1,idir,ikpt,ipert,&
186 & isppol,mband,mcgq,mcprjq,mkmem,mk1mem,&
187 & mpi_enreg,mpw,mpw1,natom,nband_k,ncpgr,&
188 & nnsclo_now,npw_k,npw1_k,nspinor,nsppol,&
189 & n4,n5,n6,occ_k,pawrhoij1,prtvol,psps,resid_k,rf_hamkq,rf_hamk_dir2,rhoaug1,rocceig,&
190 & ddk_f,wtk_k,nlines_done,cg1_out)
191 
192 
193 !This section has been created automatically by the script Abilint (TD).
194 !Do not modify the following lines by hand.
195 #undef ABI_FUNC
196 #define ABI_FUNC 'dfpt_vtowfk'
197 !End of the abilint section
198 
199  implicit none
200 
201 !Arguments ------------------------------------
202 !scalars
203  integer,intent(in) :: cplex,dim_eig2rf,ibg
204  integer,intent(in) :: ibg1,ibgq,icg,icg1,icgq,idir,ikpt,ipert,isppol
205  integer,intent(in) :: mband,mcgq,mcprjq,mk1mem,mkmem
206  integer,intent(in) :: mpw,mpw1,n4,n5,n6,natom,ncpgr
207  integer,intent(in) :: nnsclo_now,nspinor,nsppol,prtvol
208  integer,optional,intent(in) :: cg1_out
209  integer,intent(in) :: nband_k,npw1_k,npw_k
210  integer,intent(inout) :: nlines_done
211  real(dp),intent(in) :: fermie1,wtk_k
212  type(MPI_type),intent(in) :: mpi_enreg
213  type(datafiles_type),intent(in) :: dtfil
214  type(dataset_type),intent(in) :: dtset
215  type(gs_hamiltonian_type),intent(inout) :: gs_hamkq
216  type(rf_hamiltonian_type),intent(inout) :: rf_hamkq,rf_hamk_dir2
217  type(pseudopotential_type),intent(in) :: psps
218 !arrays
219  real(dp),intent(in) :: cg(2,mpw*nspinor*mband*mkmem*nsppol),cgq(2,mcgq)
220  real(dp),intent(in) :: eig0_k(nband_k),eig0_kq(nband_k)
221  real(dp),intent(in) :: ffnl1(:,:,:,:),ffnl1_test(:,:,:,:)
222  real(dp),intent(in) :: grad_berry(2,mpw1*nspinor,nband_k)
223  real(dp),intent(in) :: occ_k(nband_k),rocceig(nband_k,nband_k)
224  real(dp),intent(inout) :: cg1(2,mpw1*nspinor*mband*mk1mem*nsppol)
225  real(dp),intent(inout) :: rhoaug1(cplex*n4,n5,n6,gs_hamkq%nvloc)
226  real(dp),intent(inout) :: cg1_active(2,mpw1*nspinor*mband*mk1mem*nsppol*dim_eig2rf)
227  real(dp),intent(inout) :: gh1c_set(2,mpw1*nspinor*mband*mk1mem*nsppol*dim_eig2rf)
228  real(dp),intent(inout) :: gh0c1_set(2,mpw1*nspinor*mband*mk1mem*nsppol*dim_eig2rf)
229  real(dp),intent(inout) :: edocc_k(nband_k),eeig0_k(nband_k),eig1_k(2*nband_k**2)
230  real(dp),intent(out) :: ek0_k(nband_k),eloc0_k(nband_k)
231  real(dp),intent(inout) :: ek1_k(nband_k)
232  real(dp),intent(out) :: enl0_k(nband_k),enl1_k(nband_k)
233  real(dp),intent(out) :: resid_k(nband_k)
234  type(pawcprj_type),intent(in) :: cprj(natom,nspinor*mband*mkmem*nsppol*gs_hamkq%usecprj)
235  type(pawcprj_type),intent(in) :: cprjq(natom,mcprjq)
236  type(pawcprj_type),intent(inout) :: cprj1(natom,nspinor*mband*mk1mem*nsppol*gs_hamkq%usecprj)
237  type(pawrhoij_type),intent(inout) :: pawrhoij1(natom*gs_hamkq%usepaw)
238  type(wfk_t),intent(inout) :: ddk_f(4)
239 
240 !Local variables-------------------------------
241 !scalars
242  integer,parameter :: level=14,tim_fourwf=5
243  integer,save :: nskip=0
244  integer :: counter,iband,idir0,ierr,iexit,igs,igscq,ii,dim_dcwf,inonsc
245  integer :: iorder_cprj,iorder_cprj1,ipw,iscf_mod,ispinor,me,mgscq,nkpt_max
246  integer :: option,opt_gvnl1,quit,test_ddk
247  integer :: tocceig,usedcwavef,ptr,shift_band
248  real(dp) :: aa,ai,ar,eig0nk,resid,residk,scprod,energy_factor
249  character(len=500) :: message
250  type(rf2_t) :: rf2
251 !arrays
252  real(dp) :: tsec(2)
253  real(dp),allocatable :: cwave0(:,:),cwave1(:,:),cwavef(:,:)
254  real(dp),allocatable :: dcwavef(:,:),gh1c_n(:,:),gh0c1(:,:)
255  real(dp),allocatable :: gsc(:,:),gscq(:,:),gvnl1(:,:),gvnlc(:,:)
256  real(dp),pointer :: kinpw1(:)
257  type(pawcprj_type),allocatable :: cwaveprj(:,:),cwaveprj0(:,:),cwaveprj1(:,:)
258 
259 ! *********************************************************************
260 
261  DBG_ENTER('COLL')
262 
263 !Keep track of total time spent in dfpt_vtowfk
264  call timab(128,1,tsec)
265 
266  nkpt_max=50; if (xmpi_paral==1) nkpt_max=-1
267 
268  if(prtvol>2 .or. ikpt<=nkpt_max)then
269    write(message,'(2a,i5,2x,a,3f9.5,2x,a)')ch10,' Non-SCF iterations; k pt #',ikpt,'k=',&
270 &   gs_hamkq%kpt_k(:),'band residuals:'
271    call wrtout(std_out,message,'PERS')
272  end if
273 
274 !Initializations and allocations
275  me=mpi_enreg%me_kpt
276  quit=0
277 
278 !The value of iscf must be modified if ddk perturbation
279  iscf_mod=dtset%iscf;if(ipert==natom+1.or.ipert==natom+10.or.ipert==natom+11) iscf_mod=-3
280 
281  kinpw1 => gs_hamkq%kinpw_kp
282  ABI_ALLOCATE(gh0c1,(2,npw1_k*nspinor))
283  ABI_ALLOCATE(gvnlc,(2,npw1_k*nspinor))
284  ABI_ALLOCATE(gvnl1,(2,npw1_k*nspinor))
285  ABI_ALLOCATE(cwave0,(2,npw_k*nspinor))
286  ABI_ALLOCATE(cwavef,(2,npw1_k*nspinor))
287  ABI_ALLOCATE(cwave1,(2,npw1_k*nspinor))
288  ABI_ALLOCATE(gh1c_n,(2,npw1_k*nspinor))
289  if (gs_hamkq%usepaw==1) then
290    ABI_ALLOCATE(gsc,(2,npw1_k*nspinor))
291  else
292    ABI_ALLOCATE(gsc,(0,0))
293  end if
294 
295 !Read the npw and kg records of wf files
296  call status(0,dtfil%filstat,iexit,level,'before WffRead')
297  test_ddk=0
298  if ((ipert==natom+2.and.sum((dtset%qptn(1:3))**2)<1.0d-7.and.&
299 & (dtset%berryopt/= 4.and.dtset%berryopt/= 6.and.dtset%berryopt/= 7.and.&
300 & dtset%berryopt/=14.and.dtset%berryopt/=16.and.dtset%berryopt/=17)).or.&
301 & ipert==natom+10.or.ipert==natom+11) then
302    test_ddk=1
303    if(ipert==natom+10.or.ipert==natom+11) test_ddk=0
304  end if
305 
306 !Additional stuff for PAW
307  ABI_DATATYPE_ALLOCATE(cwaveprj0,(0,0))
308  if (gs_hamkq%usepaw==1) then
309 !  1-Compute all <g|S|Cnk+q>
310    igscq=0
311    mgscq=mpw1*nspinor*mband
312    ABI_STAT_ALLOCATE(gscq,(2,mgscq), ierr)
313    ABI_CHECK(ierr==0, "out of memory in gscq")
314 
315    call getgsc(cgq,cprjq,gs_hamkq,gscq,ibgq,icgq,igscq,ikpt,isppol,mcgq,mcprjq,&
316 &   mgscq,mpi_enreg,natom,nband_k,npw1_k,dtset%nspinor,select_k=KPRIME_H_KPRIME)
317 !  2-Initialize additional scalars/arrays
318    iorder_cprj=0;iorder_cprj1=0
319    dim_dcwf=npw1_k*nspinor;if (ipert==natom+2.or.ipert==natom+10.or.ipert==natom+11) dim_dcwf=0
320    ABI_ALLOCATE(dcwavef,(2,dim_dcwf))
321    if (gs_hamkq%usecprj==1) then
322      ABI_DATATYPE_DEALLOCATE(cwaveprj0)
323      ABI_DATATYPE_ALLOCATE(cwaveprj0,(natom,nspinor))
324      call pawcprj_alloc(cwaveprj0,1,gs_hamkq%dimcprj)
325    end if
326    ABI_DATATYPE_ALLOCATE(cwaveprj,(natom,nspinor))
327    ABI_DATATYPE_ALLOCATE(cwaveprj1,(natom,nspinor))
328    call pawcprj_alloc(cwaveprj ,0,gs_hamkq%dimcprj)
329    call pawcprj_alloc(cwaveprj1,0,gs_hamkq%dimcprj)
330  else
331    igscq=0;mgscq=0;dim_dcwf=0
332    ABI_ALLOCATE(gscq,(0,0))
333    ABI_ALLOCATE(dcwavef,(0,0))
334    ABI_DATATYPE_ALLOCATE(cwaveprj,(0,0))
335    ABI_DATATYPE_ALLOCATE(cwaveprj1,(0,0))
336  end if
337 
338  energy_factor=two
339  if(ipert==natom+10.or.ipert==natom+11) energy_factor=six
340 
341 !For rf2 perturbation :
342  if(ipert==natom+10.or.ipert==natom+11) then
343    call rf2_init(cg,cprj,rf2,dtset,dtfil,eig0_k,eig1_k,ffnl1,ffnl1_test,gs_hamkq,ibg,icg,idir,ikpt,ipert,isppol,mkmem,&
344    mpi_enreg,mpw,nband_k,nsppol,rf_hamkq,rf_hamk_dir2,occ_k,rocceig,ddk_f)
345  end if
346 
347  call timab(139,1,tsec)
348 
349 !======================================================================
350 !==================  LOOP OVER BANDS ==================================
351 !======================================================================
352 
353  do iband=1,nband_k
354 
355 !  Skip bands not treated by current proc
356    if( (mpi_enreg%proc_distrb(ikpt, iband,isppol)/=me)) cycle
357 
358 !  Get ground-state wavefunctions
359    ptr = 1+(iband-1)*npw_k*nspinor+icg
360    call cg_zcopy(npw_k*nspinor,cg(1,ptr),cwave0)
361 
362 !  Get PAW ground state projected WF (cprj)
363    if (gs_hamkq%usepaw==1.and.gs_hamkq%usecprj==1.and.ipert/=natom+10.and.ipert/=natom+11) then
364      idir0 = idir
365      if(ipert==natom+3.or.ipert==natom+4) idir0 =1
366      call pawcprj_get(gs_hamkq%atindx1,cwaveprj0,cprj,natom,iband,ibg,ikpt,iorder_cprj,&
367 &     isppol,mband,mkmem,natom,1,nband_k,nspinor,nsppol,dtfil%unpaw,&
368 &     mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb,&
369 &     icpgr=idir0,ncpgr=ncpgr)
370    end if
371 
372 !  Get first-order wavefunctions
373    ptr = 1+(iband-1)*npw1_k*nspinor+icg1
374    call cg_zcopy(npw1_k*nspinor,cg1(1,ptr),cwavef)
375 
376 !  Read PAW projected 1st-order WF (cprj)
377 !  Unuseful for the time being (will be recomputed in dfpt_cgwf)
378 !  if (gs_hamkq%usepaw==1.and.gs_hamkq%usecprj==1) then
379 !  call pawcprj_get(gs_hamkq%atindx1,cwaveprj,cprj1,natom,iband,ibg1,ikpt,iorder_cprj1,&
380 !  &    isppol,mband,mk1mem,natom,1,nband_k,nspinor,nsppol,dtfil%unpaw1,
381 !  &    mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
382 !  end if
383 
384 !  Filter the wavefunctions for large modified kinetic energy
385 !  The GS wavefunctions should already be non-zero
386    do ispinor=1,nspinor
387      igs=(ispinor-1)*npw1_k
388      do ipw=1+igs,npw1_k+igs
389        if(kinpw1(ipw-igs)>huge(zero)*1.d-11)then
390          cwavef(1,ipw)=zero
391          cwavef(2,ipw)=zero
392        end if
393      end do
394    end do
395 
396 
397 !  If electric field, the derivative of the wf should be read, and multiplied by i.
398    if(test_ddk==1) then
399      ii = wfk_findk(ddk_f(1), gs_hamkq%kpt_k)
400      ABI_CHECK(ii == ikpt, "ii != ikpt")
401      call wfk_read_bks(ddk_f(1), iband, ikpt, isppol, xmpio_single, cg_bks=gvnl1)
402 
403 !    Multiplication by -i
404 !    MVeithen 021212 : use + i instead,
405 !    See X. Gonze, Phys. Rev. B 55, 10337 (1997) [[cite:Gonze1997]] Eq. (79)
406 !    the operator used to compute the first-order derivative
407 !    of the wavefunctions with respect to an electric field
408 !    is $+i \frac{d}{dk}$
409 !    This change will affect the computation of the 2dtes from non
410 !    stationary expressions, see dfpt_nstdy.f and dfpt_nstwf.f
411      do ipw=1,npw1_k*nspinor
412 !      aa=gvnl1(1,ipw)
413 !      gvnl1(1,ipw)=gvnl1(2,ipw)
414 !      gvnl1(2,ipw)=-aa
415        aa=gvnl1(1,ipw)
416        gvnl1(1,ipw)=-gvnl1(2,ipw)
417        gvnl1(2,ipw)=aa
418      end do
419    end if
420 
421 !  Unlike in GS calculations, the inonsc loop is inside the band loop
422 !  nnsclo_now=number of non-self-consistent loops for the current vtrial
423 !  (often 1 for SCF calculation, =nstep for non-SCF calculations)
424    do inonsc=1,nnsclo_now
425 
426      counter=100*iband+inonsc
427 
428 !    Note that the following translation occurs in the called routine :
429 !    iband->band, nband_k->nband, npw_k->npw, npw1_k->npw1
430      eig0nk=eig0_k(iband)
431      usedcwavef=gs_hamkq%usepaw;if (dim_dcwf==0) usedcwavef=0
432      if (inonsc==1) usedcwavef=2*usedcwavef
433      opt_gvnl1=0;if (ipert==natom+2) opt_gvnl1=1
434      if (ipert==natom+2.and.gs_hamkq%usepaw==1.and.inonsc==1) opt_gvnl1=2
435 
436      if ( (ipert/=natom+10 .and. ipert/=natom+11) .or. abs(occ_k(iband))>tol8 ) then
437        call dfpt_cgwf(iband,dtset%berryopt,cgq,cwavef,cwave0,cwaveprj,cwaveprj0,rf2,dcwavef,&
438 &       eig0nk,eig0_kq,eig1_k,gh0c1,gh1c_n,grad_berry,gsc,gscq,gs_hamkq,gvnlc,gvnl1,icgq,&
439 &       idir,ipert,igscq,mcgq,mgscq,mpi_enreg,mpw1,natom,nband_k,dtset%nbdbuf,dtset%nline,&
440 &       npw_k,npw1_k,nspinor,opt_gvnl1,prtvol,quit,resid,rf_hamkq,dtset%dfpt_sciss,dtset%tolrde,&
441 &       dtset%tolwfr,usedcwavef,dtset%wfoptalg,nlines_done)
442        resid_k(iband)=resid
443      else
444        resid_k(iband)=zero
445      end if
446 
447      if (ipert/=natom+10 .and. ipert/= natom+11) then
448 !    At this stage, the 1st order function cwavef is orthogonal to cgq (unlike
449 !    when it is input to dfpt_cgwf). Here, restore the "active space" content
450 !    of the first-order wavefunction, to give cwave1.
451 !    PAW: note that dcwavef (1st-order change of WF due to overlap change)
452 !         remains in the subspace orthogonal to cgq
453        if (dtset%prtfull1wf>0) then
454          call full_active_wf1(cgq,cprjq,cwavef,cwave1,cwaveprj,cwaveprj1,eig1_k,fermie1,&
455 &         eig0nk,eig0_kq,dtset%elph2_imagden,iband,ibgq,icgq,mcgq,mcprjq,mpi_enreg,natom,nband_k,npw1_k,nspinor,&
456 &         0,gs_hamkq%usepaw)
457          edocc_k=zero
458          tocceig=1
459        else
460          call corrmetalwf1(cgq,cprjq,cwavef,cwave1,cwaveprj,cwaveprj1,edocc_k,eig1_k,fermie1,gh0c1,&
461 &         iband,ibgq,icgq,gs_hamkq%istwf_k,mcgq,mcprjq,mpi_enreg,natom,nband_k,npw1_k,nspinor,&
462 &         occ_k,rocceig,0,gs_hamkq%usepaw,tocceig)
463        end if
464      else
465        tocceig=0
466        call cg_zcopy(npw1_k*nspinor,cwavef,cwave1)
467        if (gs_hamkq%usepaw==1) then
468          call pawcprj_copy(cwaveprj,cwaveprj1)
469        end if
470      end if
471 
472      if (abs(occ_k(iband))<= tol8) then
473        ek0_k(iband)=zero
474        ek1_k(iband)=zero
475        eeig0_k(iband)=zero
476        enl0_k(iband)=zero
477        enl1_k(iband)=zero
478        eloc0_k(iband)=zero
479        nskip=nskip+1
480      else
481 !      Compute the 0-order kinetic operator contribution (with cwavef)
482        call meanvalue_g(ar,kinpw1,0,gs_hamkq%istwf_k,mpi_enreg,npw1_k,nspinor,cwavef,cwavef,0)
483 !      There is an additional factor of 2 with respect to the bare matrix element
484        ek0_k(iband)=energy_factor*ar
485 !      Compute the 1-order kinetic operator contribution (with cwave1 and cwave0), if needed.
486 !      Note that this is called only for ddk or strain, so that npw1_k=npw_k
487        if(ipert==natom+1 .or. ipert==natom+3 .or. ipert==natom+4)then
488          call matrixelmt_g(ai,ar,rf_hamkq%dkinpw_k,gs_hamkq%istwf_k,0,npw_k,nspinor,cwave1,cwave0,&
489 &         mpi_enreg%me_g0, mpi_enreg%comm_fft)
490 !        There is an additional factor of 4 with respect to the bare matrix element
491          ek1_k(iband)=two*energy_factor*ar
492        end if
493 
494 !      Compute eigenvalue part of total energy (with cwavef)
495        if (gs_hamkq%usepaw==1) then
496          call dotprod_g(scprod,ai,gs_hamkq%istwf_k,npw1_k*nspinor,1,cwavef,gsc,mpi_enreg%me_g0,&
497 &         mpi_enreg%comm_spinorfft)
498        else
499          call sqnorm_g(scprod,gs_hamkq%istwf_k,npw1_k*nspinor,cwavef,mpi_enreg%me_g0,&
500 &         mpi_enreg%comm_fft)
501        end if
502        eeig0_k(iband)=-energy_factor*(eig0_k(iband)- (dtset%dfpt_sciss) )*scprod
503 
504 !      Compute nonlocal psp contributions to nonlocal energy:
505 !      <G|Vnl|C1nk(perp)> is contained in gvnlc (with cwavef)
506        call dotprod_g(scprod,ai,gs_hamkq%istwf_k,npw1_k*nspinor,1,cwavef,gvnlc,mpi_enreg%me_g0,&
507 &       mpi_enreg%comm_spinorfft)
508        enl0_k(iband)=energy_factor*scprod
509 
510        if(ipert/=natom+10.and.ipert/=natom+11) then
511 !        <G|Vnl1|Cnk> is contained in gvnl1 (with cwave1)
512          call dotprod_g(scprod,ai,gs_hamkq%istwf_k,npw1_k*nspinor,1,cwave1,gvnl1,mpi_enreg%me_g0,&
513 &         mpi_enreg%comm_spinorfft)
514          enl1_k(iband)=two*energy_factor*scprod
515        end if
516 
517 !      Removal of the 1st-order kinetic energy from the 1st-order non-local part.
518        if(ipert==natom+1 .or. ipert==natom+3 .or. ipert==natom+4) then
519          enl1_k(iband)=enl1_k(iband)-ek1_k(iband)
520        end if
521 
522 !      Accumulate 1st-order density (only at the last inonsc)
523 !      Accumulate zero-order potential part of the 2nd-order total energy
524 !   BUGFIX from Max Stengel: need to initialize eloc at each inonsc iteration, in case nnonsc > 1
525        eloc0_k(iband) = zero
526        option=2;if (iscf_mod>0.and.inonsc==nnsclo_now) option=3
527        call dfpt_accrho(counter,cplex,cwave0,cwave1,cwavef,cwaveprj0,cwaveprj1,eloc0_k(iband),&
528 &       dtfil%filstat,gs_hamkq,iband,idir,ipert,isppol,dtset%kptopt,mpi_enreg,natom,nband_k,ncpgr,&
529 &       npw_k,npw1_k,nspinor,occ_k,option,pawrhoij1,prtvol,rhoaug1,tim_fourwf,tocceig,wtk_k)
530        if(ipert==natom+10.or.ipert==natom+11) eloc0_k(iband)=energy_factor*eloc0_k(iband)/two
531 
532        if(ipert==natom+10.or.ipert==natom+11) then
533          shift_band=(iband-1)*npw1_k*nspinor
534          call dotprod_g(scprod,ai,gs_hamkq%istwf_k,npw1_k*nspinor,1,cwave1,&
535 &         rf2%RHS_Stern(:,1+shift_band:npw1_k*nspinor+shift_band),mpi_enreg%me_g0, mpi_enreg%comm_spinorfft)
536          ek1_k(iband)=two*energy_factor*scprod
537        end if
538 
539      end if ! End of non-zero occupation
540 
541 !    Exit loop over inonsc if converged and if non-self-consistent
542      if (iscf_mod<0 .and. resid<dtset%tolwfr) exit
543 
544    end do ! End loop over inonsc
545 
546 !  Get first-order eigenvalues and wavefunctions
547    ptr = 1+(iband-1)*npw1_k*nspinor+icg1
548    if (.not. present(cg1_out)) then
549      call cg_zcopy(npw1_k*nspinor,cwave1,cg1(1,ptr))
550    end if
551    if(dim_eig2rf > 0) then
552      if (.not. present(cg1_out)) then
553        cg1_active(:,1+(iband-1)*npw1_k*nspinor+icg1:iband*npw1_k*nspinor+icg1)=cwavef(:,:)
554      end if
555      gh1c_set(:,1+(iband-1)*npw1_k*nspinor+icg1:iband*npw1_k*nspinor+icg1)=gh1c_n(:,:)
556      gh0c1_set(:,1+(iband-1)*npw1_k*nspinor+icg1:iband*npw1_k*nspinor+icg1)=gh0c1(:,:)
557    end if
558 
559 !  PAW: write first-order projected wavefunctions
560    if (psps%usepaw==1.and.gs_hamkq%usecprj==1) then
561      call pawcprj_put(gs_hamkq%atindx,cwaveprj,cprj1,natom,iband,ibg1,ikpt,iorder_cprj1,isppol,&
562 &     mband,mk1mem,natom,1,nband_k,gs_hamkq%dimcprj,nspinor,nsppol,dtfil%unpaw1,&
563 &     mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb,to_be_gathered=.true.)
564    end if
565 
566  end do
567 
568 !======================================================================
569 !==================  END LOOP OVER BANDS ==============================
570 !======================================================================
571 
572 !For rf2 perturbation
573  if(ipert==natom+10.or.ipert==natom+11) call rf2_destroy(rf2)
574 
575 !Find largest resid over bands at this k point
576  residk=maxval(resid_k(:))
577  if (prtvol>2 .or. ikpt<=nkpt_max) then
578    do ii=0,(nband_k-1)/8
579      write(message,'(1p,8e10.2)')(resid_k(iband),iband=1+ii*8,min(nband_k,8+ii*8))
580      call wrtout(std_out,message,'PERS')
581    end do
582  end if
583 
584  call timab(139,2,tsec)
585  call timab(130,1,tsec)
586 
587  ABI_DEALLOCATE(cwave0)
588  ABI_DEALLOCATE(cwavef)
589  ABI_DEALLOCATE(cwave1)
590  ABI_DEALLOCATE(gh0c1)
591  ABI_DEALLOCATE(gvnlc)
592  ABI_DEALLOCATE(gvnl1)
593  ABI_DEALLOCATE(gh1c_n)
594 
595  if (gs_hamkq%usepaw==1) then
596    call pawcprj_free(cwaveprj)
597    call pawcprj_free(cwaveprj1)
598    if (gs_hamkq%usecprj==1) then
599      call pawcprj_free(cwaveprj0)
600    end if
601  end if
602  ABI_DEALLOCATE(dcwavef)
603  ABI_DEALLOCATE(gscq)
604  ABI_DEALLOCATE(gsc)
605  ABI_DATATYPE_DEALLOCATE(cwaveprj0)
606  ABI_DATATYPE_DEALLOCATE(cwaveprj)
607  ABI_DATATYPE_DEALLOCATE(cwaveprj1)
608 
609 !###################################################################
610 
611 !Write the number of one-way 3D ffts skipped until now (in case of fixed occupation numbers)
612  if(iscf_mod>0 .and. (prtvol>2 .or. ikpt<=nkpt_max))then
613    write(message,'(a,i0)')' dfpt_vtowfk : number of one-way 3D ffts skipped in vtowfk3 until now =',nskip
614    call wrtout(std_out,message,'PERS')
615  end if
616 
617  if(prtvol<=2 .and. ikpt==nkpt_max+1)then
618    write(message,'(3a)') ch10,' dfpt_vtowfk : prtvol=0, 1 or 2, do not print more k-points.',ch10
619    call wrtout(std_out,message,'PERS')
620  end if
621 
622  if (residk>dtset%tolwfr .and. iscf_mod<=0 .and. iscf_mod/=-3) then
623    write(message,'(a,2i0,a,es13.5)')'Wavefunctions not converged for nnsclo,ikpt=',nnsclo_now,&
624 &   ikpt,' max resid=',residk
625    MSG_WARNING(message)
626  end if
627 
628  call timab(130,2,tsec)
629  call timab(128,2,tsec)
630 
631  DBG_EXIT('COLL')
632 
633 end subroutine dfpt_vtowfk

ABINIT/full_active_wf1 [ Functions ]

[ Top ] [ Functions ]

NAME

 full_active_wf1

FUNCTION

 Response function calculation only:
 Restore the full "active space" contribution to the 1st-order wavefunctions.
 The 1st-order WF corrected in this way will no longer be othogonal to the other occupied states.
 This routine will be only used in a non self-consistent calculation of the
 1st-order WF for post-processing purposes. Therefore, it does not compute
 the contribution of the 2DTE coming from the change of occupations.

INPUTS

  cg(2,mcgq)=planewave coefficients of wavefunctions at k+q
  cprjq(natom,mcprjq)= wave functions at k+q projected with non-local projectors
  cwavef(2,npw1*nspinor)= 1st-order wave-function before correction
  cwaveprj(natom,nspinor)= 1st-order wave-function before correction
                           projected on NL projectors (PAW)
  eig1(2*nband**2)=first-order eigenvalues (hartree)
  fermie1=derivative of fermi energy wrt (strain) perturbation
  eig0nk=energy of the band at k being corrected
  eig0_kq(nband)=energies of the bands at k+q
  elph2_imagden=imaginary parameter to broaden the energy denominators
  iband=index of current band
  ibgq=shift to be applied on the location of data in the array cprjq
  icgq=shift to be applied on the location of data in the array cgq
  mcgq=second dimension of the cgq array
  mcprjq=second dimension of the cprjq array
  mpi_enreg=information about MPI parallelization
  natom=number of atoms in cell
  nband=number of bands
  npw1=number of plane waves at this k+q point
  nspinor=number of spinorial components of the wavefunctions
  timcount=index used to accumulate timing (0 from dfpt_vtowfk, 1 from dfpt_nstwf)
  usepaw=flag for PAW

OUTPUT

  cwave1(2,npw1*nspinor)= 1st-order wave-function after correction
  cwaveprj1(natom,nspinor)= 1st-order wave-function after correction
                            projected on NL projectors (PAW)

PARENTS

      dfpt_vtowfk

CHILDREN

      cg_zcopy,dotprod_g,pawcprj_copy,pawcprj_zaxpby,timab

SOURCE

686 subroutine full_active_wf1(cgq,cprjq,cwavef,cwave1,cwaveprj,cwaveprj1,eig1,&
687 &               fermie1,eig0nk,eig0_kq,elph2_imagden,&
688 &               iband,ibgq,icgq,mcgq,mcprjq,mpi_enreg,natom,nband,npw1,&
689 &               nspinor,timcount,usepaw)
690 
691 
692 !This section has been created automatically by the script Abilint (TD).
693 !Do not modify the following lines by hand.
694 #undef ABI_FUNC
695 #define ABI_FUNC 'full_active_wf1'
696 !End of the abilint section
697 
698  implicit none
699 
700 !Arguments ------------------------------------
701 !scalars
702  integer,intent(in) :: iband,ibgq,icgq,mcgq,mcprjq,natom,nband,npw1,nspinor,timcount,usepaw
703  real(dp),intent(in) :: fermie1, eig0nk
704  real(dp),intent(in) :: elph2_imagden
705  type(MPI_type),intent(in) :: mpi_enreg
706 !arrays
707  real(dp),intent(in) :: cgq(2,mcgq),cwavef(2,npw1*nspinor)
708  real(dp),intent(in) :: eig0_kq(nband)
709  real(dp),intent(in) :: eig1(2*nband**2)
710  real(dp),intent(out) :: cwave1(2,npw1*nspinor)
711  type(pawcprj_type),intent(in) :: cprjq(natom,mcprjq),cwaveprj(natom,nspinor*usepaw)
712  type(pawcprj_type),intent(inout) :: cwaveprj1(natom,nspinor*usepaw) !vz_i
713 
714 !Local variables-------------------------------
715 !scalars
716  integer :: ibandkq,index_cgq,index_cprjq,index_eig1,ii
717  real(dp) :: facti,factr,invocc,eta,delta_E,inv_delta_E,gkkr
718 !arrays
719  real(dp) :: tsec(2)
720  real(dp),allocatable :: cwcorr(:,:)
721 
722 ! *********************************************************************
723 
724  DBG_ENTER("COLL")
725 
726  call timab(214+timcount,1,tsec)
727 
728 !At this stage, the 1st order function cwavef is orthogonal to cgq (unlike when it is input to dfpt_cgwf).
729 !Here, restore the "active space" content of the 1st-order wavefunction, to give cwave1 .
730 
731 !First copy input WF into output WF
732  call cg_zcopy(npw1*nspinor,cwavef,cwave1)
733 
734  if (usepaw==1) then
735    call pawcprj_copy(cwaveprj,cwaveprj1)
736  end if
737 
738  eta = elph2_imagden
739 
740 !Loop over WF at k+q subspace
741  do ibandkq=1,nband
742 
743    delta_E = eig0nk - eig0_kq(ibandkq)
744    inv_delta_E = delta_E / ( delta_E ** 2 + eta ** 2)
745 
746    index_eig1=2*ibandkq-1+(iband-1)*2*nband
747    index_cgq=npw1*nspinor*(ibandkq-1)+icgq
748 
749    if(ibandkq==iband) then
750      gkkr = eig1(index_eig1) - fermie1
751    else
752      gkkr = eig1(index_eig1)
753    end if
754    factr = inv_delta_E * gkkr
755    facti = inv_delta_E * eig1(index_eig1+1)
756 
757 !  Apply correction to 1st-order WF
758 !$OMP PARALLEL DO PRIVATE(ii) SHARED(cgq,cwave1,facti,factr,index_cgq,npw1,nspinor)
759    do ii=1,npw1*nspinor
760      cwave1(1,ii)=cwave1(1,ii)+(factr*cgq(1,ii+index_cgq)-facti*cgq(2,ii+index_cgq))
761      cwave1(2,ii)=cwave1(2,ii)+(facti*cgq(1,ii+index_cgq)+factr*cgq(2,ii+index_cgq))
762    end do
763 
764 !  In the PAW case, also apply correction to projected WF
765    if (usepaw==1) then
766      index_cprjq=nspinor*(ibandkq-1)+ibgq
767      call pawcprj_zaxpby((/factr,facti/),(/one,zero/),cprjq(:,index_cprjq+1:index_cprjq+nspinor),cwaveprj1)
768    end if
769 
770  end do ! Loop over k+q subspace
771 
772  call timab(214+timcount,2,tsec)
773 
774  DBG_EXIT("COLL")
775 
776 end subroutine full_active_wf1

ABINIT/m_dfpt_vtowfk [ Modules ]

[ Top ] [ Modules ]

NAME

  m_dfpt_vtowfk

FUNCTION

COPYRIGHT

  Copyright (C) 1999-2018 ABINIT group (XG, AR, DRH, MB, MVer,XW, MT, GKA)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_dfpt_vtowfk
28 
29  use defs_basis
30  use defs_datatypes
31  use defs_abitypes
32  use m_abicore
33  use m_errors
34  use m_xmpi
35  use m_cgtools
36  use m_wfk
37  use m_rf2
38  use m_rf2_init,         only : rf2_init
39 
40  use m_dtfil,        only : status
41  use m_time,         only : timab
42  use m_pawrhoij,     only : pawrhoij_type
43  use m_pawcprj,      only : pawcprj_type, pawcprj_alloc, pawcprj_put, pawcprj_free, pawcprj_get, pawcprj_copy, pawcprj_zaxpby
44  use m_hamiltonian,  only : gs_hamiltonian_type, rf_hamiltonian_type, KPRIME_H_KPRIME
45  use m_spacepar,     only : meanvalue_g
46  use m_dfpt_mkrho,   only : dfpt_accrho
47  use m_dfpt_cgwf,    only : dfpt_cgwf
48  use m_getghc,       only : getgsc
49 
50  implicit none
51 
52  private