TABLE OF CONTENTS


ABINIT/dfpt_cgwf [ Functions ]

[ Top ] [ Functions ]

NAME

 dfpt_cgwf

FUNCTION

 Update one single wavefunction (cwavef), non self-consistently.
 Uses a conjugate-gradient algorithm.
 Try to keep close to the formulas in PRB55, 10337 (1997) [[cite:Gonze1997]], for the
 non-self-consistent case, except that we are computing here
 the second-derivative of the total energy, and not E(2). There
 is a factor of 2 between the two quantities ...
 The wavefunction that is generated is always orthogonal to cgq .
 It is orthogonal to the active Hilbert space, and will be complemented
 by contributions from the active space in the calling routine, if needed.

INPUTS

  band=which particular band we are converging.
  berryopt=option for Berry phase
  cgq(2,mcgq)=wavefunction coefficients for ALL bands at k+Q
  cwave0(2,npw*nspinor)=GS wavefunction at k, in reciprocal space
  cwaveprj0(natom,nspinor*usecprj)=GS wave function at k projected with nl projectors
  eig0nk=0-order eigenvalue for the present wavefunction at k
  eig0_kq(nband)=GS eigenvalues at k+Q (hartree)
  grad_berry(2,mpw1,dtefield%mband_occ) = the gradient of the Berry phase term
  gscq(2,mgscq)=<g|S|Cnk+q> coefficients for ALL bands (PAW) at k+Q
  gs_hamkq <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k+Q
  icgq=shift to be applied on the location of data in the array cgq
  igscq=shift to be applied on the location of data in the array gscq
  idir=direction of the perturbation
  ipert=type of the perturbation
  mcgq=second dimension of the cgq array
  mgscq=second dimension of gscq
  mpi_enreg=information about MPI parallelization
  mpw1=maximum number of planewave for first-order wavefunctions
  natom=number of atoms in cell.
  nband=number of bands.
  nbdbuf=number of buffer bands for the minimisation
  nline=number of line minimizations per band.
  npw=number of planewaves in basis sphere at given k.
  npw1=number of planewaves in basis sphere at k+Q
  nspinor=number of spinorial components of the wavefunctions
  opt_gvnl1=option controlling the use of gvnl1 array:
            0: used as an output
            1: used as an input:    - used only for ipert=natom+2
                 NCPP: contains the ddk 1-st order WF
                 PAW: contains frozen part of 1st-order hamiltonian
            2: used as input/ouput:    - used only for PAW and ipert=natom+2
                 At input: contains the ddk 1-st order WF (times i)
                 At output: contains frozen part of 1st-order hamiltonian
  prtvol=control print volume and debugging output
  quit= if 1, proceeds to smooth ending of the job.
  dfpt_sciss=scissor shift (Ha)
  tolrde=tolerance on the ratio of differences of energies (for the line minimisation)
  tolwfr=tolerance on largest wf residual
  usedcwavef=flag controlling the use of dcwavef array (PAW only):
             0: not used (not allocated)
             1: used as input
             2: used as output
  wfoptalg=govern the choice of algorithm for wf optimisation (0 or 10, at present)

OUTPUT

  eig1_k(2*nband**2)=matrix of first-order eigenvalues (hartree)
                     eig1(:,ii,jj)=<C0 ii|H1|C0 jj> for norm-conserving psps
                     eig1(:,ii,jj)=<C0 ii|H1-(eig0_k+eig0_k+q)/2.S(1)|C0 jj> for PAW
  ghc(2,npw1*nspinor)=<G|H0-eig0_k.I|C1 band,k> (NCPP) or <G|H0-eig0_k.S0|C1 band,k> (PAW)
  gvnlc(2,npw1*nspinor)=<G|Vnl|C1 band,k>
  gvnl1(2,npw1*nspinor)=  part of <G|K1+Vnl1|C0 band,k> not depending on VHxc1           (NCPP)
                       or part of <G|K1+Vnl1-eig0k.S1|C0 band,k> not depending on VHxc1 (PAW)
  resid=wf residual for current band
  gh1c_n= <G|H1|C0 band,k> (NCPP) or <G|H1-eig0k.S1|C0 band,k> (PAW).
          This vector is not projected on the subspace orthogonal to the WF.
  === if gs_hamkq%usepaw==1 ===
  gsc(2,npw1*nspinor*usepaw)=<G|S0|C1 band,k>

SIDE EFFECTS

  Input/Output:
  cwavef(2,npw1*nspinor)=first-order  wavefunction at k,q, in reciprocal space (updated)
  === if gs_hamkq%usepaw==1 ===
  cwaveprj(natom,nspinor)= wave functions at k projected with nl projectors
  === if usedcwavef>0 ===
  dcwavef(2,npw1*nspinor)=change of wavefunction due to change of overlap:
         dcwavef is delta_Psi(1)=-1/2.Sum_{j}[<C0_k+q_j|S(1)|C0_k_i>.|C0_k+q_j>]
         see PRB 78, 035105 (2008) [[cite:Audouze2008]], Eq. (42)
         input if usedcwavef=1, output if usedcwavef=2

PARENTS

      dfpt_vtowfk

CHILDREN

      cg_precon,cg_zaxpy,cg_zcopy,dotprod_g,getdc1,getgh1c,getghc
      pawcprj_alloc,pawcprj_axpby,pawcprj_free,pawcprj_set_zero,projbd
      sqnorm_g,timab,wrtout

SOURCE

 150 subroutine dfpt_cgwf(band,berryopt,cgq,cwavef,cwave0,cwaveprj,cwaveprj0,rf2,dcwavef,&
 151 & eig0nk,eig0_kq,eig1_k,ghc,gh1c_n,grad_berry,gsc,gscq,&
 152 & gs_hamkq,gvnlc,gvnl1,icgq,idir,ipert,igscq,&
 153 & mcgq,mgscq,mpi_enreg,mpw1,natom,nband,nbdbuf,nline_in,npw,npw1,nspinor,&
 154 & opt_gvnl1,prtvol,quit,resid,rf_hamkq,dfpt_sciss,tolrde,tolwfr,&
 155 & usedcwavef,wfoptalg,nlines_done)
 156 
 157 
 158 !This section has been created automatically by the script Abilint (TD).
 159 !Do not modify the following lines by hand.
 160 #undef ABI_FUNC
 161 #define ABI_FUNC 'dfpt_cgwf'
 162 !End of the abilint section
 163 
 164  implicit none
 165 
 166 !Arguments ------------------------------------
 167 !scalars
 168  integer,intent(in) :: band,berryopt
 169  integer,intent(in) :: icgq,idir,igscq,ipert,mcgq,mgscq,mpw1,natom,nband
 170  integer,intent(in) :: nbdbuf,nline_in,npw,npw1,nspinor,opt_gvnl1
 171  integer,intent(in) :: prtvol,quit,usedcwavef,wfoptalg
 172  integer,intent(inout) :: nlines_done
 173  real(dp),intent(in) :: eig0nk,dfpt_sciss,tolrde,tolwfr
 174  real(dp),intent(out) :: resid
 175  type(MPI_type),intent(in) :: mpi_enreg
 176  type(rf2_t), intent(in) :: rf2
 177  type(gs_hamiltonian_type),intent(inout) :: gs_hamkq
 178  type(rf_hamiltonian_type),intent(inout) :: rf_hamkq
 179 !arrays
 180  real(dp),intent(in) :: cgq(2,mcgq),eig0_kq(nband)
 181  real(dp),intent(in) :: grad_berry(2,mpw1*nspinor,nband),gscq(2,mgscq)
 182  real(dp),intent(inout) :: cwave0(2,npw*nspinor),cwavef(2,npw1*nspinor)
 183  real(dp),intent(inout) :: dcwavef(2,npw1*nspinor*((usedcwavef+1)/2))
 184  real(dp),intent(inout) :: eig1_k(2*nband**2)
 185  real(dp),intent(out) :: gh1c_n(2,npw1*nspinor)
 186  real(dp),intent(out) :: ghc(2,npw1*nspinor)
 187  real(dp),intent(out) :: gsc(2,npw1*nspinor*gs_hamkq%usepaw)
 188  real(dp),intent(inout) :: gvnl1(2,npw1*nspinor),gvnlc(2,npw1*nspinor)
 189  type(pawcprj_type),intent(inout) :: cwaveprj(natom,nspinor)
 190  type(pawcprj_type),intent(inout) :: cwaveprj0(natom,nspinor*gs_hamkq%usecprj)
 191 
 192 !Local variables-------------------------------
 193 !scalars
 194  integer,parameter :: level=15,tim_getgh1c=1,tim_getghc=2,tim_projbd=2
 195  integer,save :: nskip=0
 196  integer :: cpopt,iband,igs,iline,indx_cgq,ipw,me_g0,comm_fft
 197  integer :: ipws,ispinor,istwf_k,jband,nline,optlocal,optnl,shift_band,sij_opt
 198  integer :: test_is_ok,useoverlap,usepaw,usevnl,usetolrde
 199  real(dp) :: d2edt2,d2te,d2teold,dedt,deltae,deold,dotgg
 200  real(dp) :: dotgp,doti,dotr,eshift,eshiftkq,gamma,optekin,prod1,prod2
 201  real(dp) :: theta,tol_restart,u1h0me0u1
 202  logical :: gen_eigenpb
 203  character(len=500) :: msg
 204 !arrays
 205  real(dp) :: dummy(0,0),tsec(2)
 206  real(dp),allocatable :: conjgr(:,:),cwaveq(:,:),cwwork(:,:),direc(:,:)
 207  real(dp),allocatable :: gberry(:,:),gh1c(:,:),gh_direc(:,:),gresid(:,:),gvnl1_saved(:,:)
 208  real(dp),allocatable :: gs1c(:,:),gvnl_direc(:,:),pcon(:),sconjgr(:,:)
 209  real(dp),allocatable :: scprod(:,:),work(:,:),work1(:,:),work2(:,:)
 210  real(dp),pointer :: kinpw1(:)
 211  type(pawcprj_type),allocatable :: conjgrprj(:,:)
 212  type(pawcprj_type) :: cprj_dummy(0,0)
 213 
 214 ! *********************************************************************
 215 
 216  DBG_ENTER("COLL")
 217 
 218  call timab(122,1,tsec)
 219 
 220 !======================================================================
 221 !========= LOCAL VARIABLES DEFINITIONS AND ALLOCATIONS ================
 222 !====================================================================
 223 
 224  nline = nline_in
 225  usetolrde = 1
 226 
 227 ! LB-23/04/17:
 228 ! For ipert=natom+10 or ipert=natom+11, the Sternheimer equation is non-self-consistent, so we have
 229 ! to solve a true linear problem (A.X = B) for each kpoint and band. In this case, the conjugate
 230 ! gradient algorithm can find the exact solution (within the numerical precision) with only ONE call
 231 ! of dfpt_cgwf (per kpoint and band...). The solution is found with at most N iterations, N being the dimension of X.
 232 ! In order to avoid useless scfcv loops (and many calls of rf2_init, which can be time consuming),
 233 ! we want to leave this routine only if 'tolwfr' is reached. Consequently, 'tolrde' is not used and 'nline' is set to 100.
 234 ! One could use nline=npw1*nspinor (>> 100 !) instead, but when the method cannot converge (i.e when tolwfr is lower than the
 235 ! numerical noise) the program could be stuck here for a very long time.
 236 ! NOTE : This is also true for ipert==natom+1, but a lot of references in the test suite have to be changed...
 237  if(ipert==natom+10.or.ipert==natom+11) then
 238    nline = 100 ! The default value is only 4... This should be sufficient to converge with nstep=1 or 2
 239    if (nline_in > 100) nline = nline_in ! Keep the possibility to increase nline
 240    usetolrde = 0  ! see below
 241  end if
 242 
 243 !Tell us what is going on:
 244  if (prtvol>=10) then
 245    write(msg,'(a,i6,2x,a,i3,a)')' --- dfpt_cgwf is called for band',band,'for',nline,' lines'
 246    call wrtout(std_out,msg,'PERS')
 247  end if
 248 
 249  me_g0 = mpi_enreg%me_g0
 250  comm_fft = mpi_enreg%comm_fft
 251 
 252 !if PAW, one has to solve a generalized eigenproblem
 253  usepaw=gs_hamkq%usepaw
 254  gen_eigenpb=(usepaw==1)
 255  useoverlap=0;if (gen_eigenpb) useoverlap=1
 256 
 257 !Use scissor shift on 0-order eigenvalue
 258  eshift=eig0nk-dfpt_sciss
 259 
 260 !Additional initializations
 261  istwf_k=gs_hamkq%istwf_k
 262  optekin=0;if (wfoptalg>=10) optekin=1
 263  tol_restart=tol12;if (gen_eigenpb) tol_restart=tol8
 264  if (ipert == natom+10 .or. ipert == natom+11) then
 265    tol_restart = tol7
 266  end if
 267  kinpw1 => gs_hamkq%kinpw_kp
 268 
 269 !Memory allocations
 270  ABI_ALLOCATE(gh1c,(2,npw1*nspinor))
 271  ABI_ALLOCATE(pcon,(npw1))
 272  ABI_ALLOCATE(scprod,(2,nband))
 273 
 274  if (berryopt== 4.or.berryopt== 6.or.berryopt== 7.or.&
 275 & berryopt==14.or.berryopt==16.or.berryopt==17) then
 276    ABI_ALLOCATE(gberry,(2,npw1*nspinor))
 277    gberry(:,1:npw1*nspinor)=grad_berry(:,1:npw1*nspinor,band)
 278  else
 279    ABI_ALLOCATE(gberry,(0,0))
 280  end if
 281 
 282  shift_band=(band-1)*npw1*nspinor
 283 
 284 !Several checking statements
 285  if (prtvol==-level.or.prtvol==-19.or.prtvol==-20) then
 286    write(msg,'(a)') " ** cgwf3 : debugging mode, tests will be done"
 287    ! Search CGWF3_WARNING in the log file to find errors (if any)
 288    call wrtout(std_out,msg,'PERS')
 289    ABI_ALLOCATE(work,(2,npw1*nspinor))
 290    ABI_ALLOCATE(work1,(2,npw1*nspinor))
 291 !  ===== Check <Psi_k+q^(0)|S(0)|Psi_k+q^(0)>=delta
 292    if (.not.gen_eigenpb) work1(:,:)=cgq(:,1+npw1*nspinor*(band-1)+icgq:npw1*nspinor*band+icgq)
 293    if (     gen_eigenpb) work1(:,:)=gscq(:,1+npw1*nspinor*(band-1)+igscq:npw1*nspinor*band+igscq)
 294    do iband=1,nband
 295      work(:,:)=cgq(:,1+npw1*nspinor*(iband-1)+icgq:npw1*nspinor*iband+icgq)
 296      call dotprod_g(dotr,doti,istwf_k,npw1*nspinor,2,work1,work,me_g0,mpi_enreg%comm_spinorfft)
 297      test_is_ok=1
 298      if(iband==band) then
 299        if(abs(dotr-1)>tol12) test_is_ok=0
 300      else
 301        if(abs(dotr)>tol12) test_is_ok=0
 302      end if
 303      if(abs(doti)>tol12) test_is_ok=0
 304      if(test_is_ok/=1) then
 305        write(msg,'(a,i3,a,2es22.15)') "CGWF3_WARNING : <Psi_k+q,i^(0)|S(0)|Psi_k+q,j^(0)> for band j=",iband," is ",dotr,doti
 306        call wrtout(std_out,msg,'PERS')
 307      end if
 308    end do
 309 !  ===== Check Pc.Psi_k+q^(0)=0
 310    do iband=1,nband
 311      work(:,:)=cgq(:,1+npw1*nspinor*(iband-1)+icgq:npw1*nspinor*iband+icgq)
 312      call projbd(cgq,work,-1,icgq,igscq,istwf_k,mcgq,mgscq,nband,npw1,nspinor,&
 313 &     gscq,scprod,0,tim_projbd,useoverlap,me_g0,comm_fft)
 314      call sqnorm_g(dotr,istwf_k,npw1*nspinor,work,me_g0,comm_fft)
 315      if(sqrt(dotr)>tol12) then
 316        write(msg,'(a,i3,a,es22.15)') "CGWF3_WARNING : Norm of Pc.Psi_k+q_j^(0) for band j=",iband," is ",sqrt(dotr)
 317        call wrtout(std_out,msg,'PERS')
 318      end if
 319    end do
 320 !  ===== Check Pc.Psi_k^(0)=0
 321    work(:,:)=cwave0(:,:)
 322    call projbd(cgq,work,-1,icgq,igscq,istwf_k,mcgq,mgscq,nband,npw1,nspinor,&
 323 &   gscq,scprod,0,tim_projbd,useoverlap,me_g0,comm_fft)
 324    call sqnorm_g(dotr,istwf_k,npw1*nspinor,work,me_g0,comm_fft)
 325    if(sqrt(dotr)>tol12) then
 326      write(msg,'(a,i3,a,es22.15)') "CGWF3_WARNING : Norm of Pc.Psi_k^(0) for band ",band," is ",sqrt(dotr)
 327      call wrtout(std_out,msg,'PERS')
 328    end if
 329 !  ===== Check Pc.dcwavef=0 (for 2nd order only)
 330    if(ipert==natom+10.or.ipert==natom+11) then
 331      work(:,:)=rf2%dcwavef(:,1+shift_band:npw1*nspinor+shift_band)
 332      call projbd(cgq,work,-1,icgq,igscq,istwf_k,mcgq,mgscq,nband,npw1,nspinor,&
 333 &     gscq,scprod,0,tim_projbd,useoverlap,me_g0,comm_fft)
 334      call sqnorm_g(dotr,istwf_k,npw1*nspinor,work,me_g0,comm_fft)
 335      if(sqrt(dotr)>tol10) then
 336        write(msg,'(a,i3,a,es22.15)') "CGWF3_WARNING : Norm of Pc.dcwavef for band ",band," is ",sqrt(dotr)
 337        call wrtout(std_out,msg,'PERS')
 338      end if
 339    end if
 340 !  ===== Check Pc^*.S(0).Psi_k+q^(0)=0
 341    if (gen_eigenpb) then
 342      do iband=1,nband
 343        work(:,:)=gscq(:,1+npw1*nspinor*(iband-1)+igscq:npw1*nspinor*iband+igscq)
 344        call projbd(gscq,work,-1,igscq,icgq,istwf_k,mgscq,mcgq,nband,npw1,nspinor,&
 345 &       cgq,scprod,0,tim_projbd,useoverlap,me_g0,comm_fft)
 346        call sqnorm_g(dotr,istwf_k,npw1*nspinor,work,me_g0,comm_fft)
 347        if(sqrt(dotr)>tol12) then
 348          write(msg,'(a,i3,a,es22.15)') "CGWF3_WARNING : Norm of Pc^*.S(0).Psi_k+q_j^(0) for band j=",iband," is ",sqrt(dotr)
 349          call wrtout(std_out,msg,'PERS')
 350        end if
 351      end do
 352    end if
 353    ABI_DEALLOCATE(work)
 354    ABI_DEALLOCATE(work1)
 355  end if
 356 
 357 !======================================================================
 358 !========== INITIALISATION OF MINIMIZATION ITERATIONS =================
 359 !======================================================================
 360 
 361  if(ipert/=natom+10.and.ipert/=natom+11) then
 362 !  The following is needed for first order perturbations only
 363 !  Otherwise, the work is already done in rf2_init (called in dfpt_vtowfk.F90)
 364 
 365 !  Compute H(1) applied to GS wavefunction Psi(0)
 366    if (gen_eigenpb) then
 367      sij_opt=1
 368      ABI_ALLOCATE(gs1c,(2,npw1*nspinor))
 369    else
 370      ABI_ALLOCATE(gs1c,(0,0))
 371      sij_opt=0
 372    end if
 373    usevnl=1; optlocal=1; optnl=2
 374    if (prtvol==-level.or.prtvol==-19) then
 375      ABI_ALLOCATE(gvnl1_saved,(2,npw1*nspinor))
 376      gvnl1_saved(:,:) = gvnl1(:,:)
 377    end if
 378    call getgh1c(berryopt,cwave0,cwaveprj0,gh1c,gberry,gs1c,gs_hamkq,gvnl1,idir,ipert,eshift,&
 379 &   mpi_enreg,optlocal,optnl,opt_gvnl1,rf_hamkq,sij_opt,tim_getgh1c,usevnl)
 380 
 381    if (gen_eigenpb) then
 382      if (ipert/=natom+2) then  ! S^(1) is zero for ipert=natom+2
 383 !$OMP PARALLEL
 384 !$OMP DO
 385        do ipw=1,npw1*nspinor
 386          gh1c (1:2,ipw)=gh1c (1:2,ipw)-eshift*gs1c(1:2,ipw)
 387        end do
 388 !$OMP END DO NOWAIT
 389        if (opt_gvnl1/=1) then
 390 !$OMP DO
 391          do ipw=1,npw1*nspinor
 392            gvnl1(1:2,ipw)=gvnl1(1:2,ipw)-eshift*gs1c(1:2,ipw)
 393          end do
 394 !$OMP END DO NOWAIT
 395        end if
 396 !$OMP END PARALLEL
 397      end if
 398 
 399 !    If generalized eigenPb and dcwavef requested, compute it:
 400 !    dcwavef is delta_Psi(1)=-1/2.Sum_{j}[<C0_k+q_j|S(1)|C0_k_i>.|C0_k+q_j>]
 401 !    see PRB 78, 035105 (2008) [[cite:Audouze2008]], Eq. (42)
 402      if (usedcwavef==2) then
 403        call getdc1(cgq,cprj_dummy,dcwavef,cprj_dummy,0,icgq,istwf_k,mcgq,0,&
 404 &       mpi_enreg,natom,nband,npw1,nspinor,0,gs1c)
 405      end if
 406    end if
 407 
 408  else ! 2nd order case (wrt k perturbation)
 409    ! Copy RHS_Stern(:,:) of the given band in gh1c
 410    gh1c(:,:)=rf2%RHS_Stern(:,1+shift_band:npw1*nspinor+shift_band)
 411  end if
 412 
 413 !Check that Pc^*.(H^(0)-E.S^(0)).delta_Psi^(1) is zero ! This is a consequence of P_c delta_Psi^(1) = 0
 414  if (prtvol==-level.and.usedcwavef==2) then
 415    ABI_ALLOCATE(cwwork,(2,npw1*nspinor))
 416    cwwork=dcwavef
 417 !  - Apply H^(0)-E.S^(0) to delta_Psi^(1)
 418    sij_opt=0;if (gen_eigenpb) sij_opt=-1
 419    cpopt=-1
 420    ABI_ALLOCATE(work,(2,npw1*nspinor))
 421    ABI_ALLOCATE(work1,(2,npw1*nspinor*((sij_opt+1)/2)))
 422    ABI_ALLOCATE(work2,(2,npw1*nspinor))
 423    call getghc(cpopt,cwwork,conjgrprj,work,work1,gs_hamkq,work2,eshift,mpi_enreg,&
 424 &   1,prtvol,sij_opt,tim_getghc,0,select_k=KPRIME_H_KPRIME)
 425    cwwork=work
 426    ABI_DEALLOCATE(work)
 427    ABI_DEALLOCATE(work1)
 428    ABI_DEALLOCATE(work2)
 429 !  -Apply Pc^*
 430    call projbd(gscq,cwwork,-1,igscq,icgq,istwf_k,mgscq,mcgq,nband,npw1,nspinor,&
 431 &   cgq,scprod,0,tim_projbd,useoverlap,me_g0,comm_fft)
 432    call sqnorm_g(dotr,istwf_k,npw1*nspinor,cwwork,me_g0,comm_fft)
 433    ABI_DEALLOCATE(cwwork)
 434    if(sqrt(dotr)>tol12) then
 435      write(msg,'(a,i3,a,es22.15)') 'CGWF3_WARNING : |Pc^*.(H^(0)-E.S^(0)).delta_Psi^(1)| (band ',band,')=',sqrt(dotr)
 436      call wrtout(std_out,msg,'PERS')
 437    end if
 438  end if
 439 
 440  call cg_zcopy(npw1*nspinor,gh1c,gh1c_n)
 441 
 442 !Projecting out all bands
 443 !While we could avoid calculating all the eig1_k to obtain the perturbed density,
 444 !we do need all of the matrix element when outputing the full 1st-order wfn.
 445 !Note the subtlety:
 446 !-For the generalized eigenPb, S|cgq> is used in place of |cgq>,
 447 !in order to apply P_c+ projector (see PRB 73, 235101 (2006) [[cite:Audouze2006]], Eq. (71), (72))
 448  if(gen_eigenpb)then
 449    call projbd(gscq,gh1c,-1,igscq,icgq,istwf_k,mgscq,mcgq,nband,npw1,nspinor,&
 450 &   cgq,scprod,0,tim_projbd,useoverlap,me_g0,comm_fft)
 451  else
 452    call projbd(cgq,gh1c,-1,icgq,0,istwf_k,mcgq,mgscq,nband,npw1,nspinor,&
 453 &   dummy,scprod,0,tim_projbd,useoverlap,me_g0,comm_fft)
 454  end if
 455 
 456  if(ipert/=natom+10.and.ipert/=natom+11) then
 457  ! For ipert=natom+10 or natom+11, this is done in rf2_init
 458 
 459    !The array eig1_k contains:
 460    !<u_(iband,k+q)^(0)|H_(k+q,k)^(1)|u_(jband,k)^(0)>                           (NC psps)
 461    !or <u_(iband,k+q)^(0)|H_(k+q,k)^(1)-(eig0_k+eig0_k+q)/2.S^(1)|u_(jband,k)^(0)> (PAW)
 462    jband=(band-1)*2*nband
 463    if (gen_eigenpb) then
 464      indx_cgq=icgq
 465      do iband=1,nband
 466        eshiftkq=half*(eig0_kq(iband)-eig0nk)
 467        call dotprod_g(dotr,doti,istwf_k,npw1*nspinor,2,cgq(:,indx_cgq+1:indx_cgq+npw1*nspinor),gs1c,&
 468 &       me_g0,mpi_enreg%comm_spinorfft)
 469        eig1_k(2*iband-1+jband)=scprod(1,iband)-eshiftkq*dotr
 470        eig1_k(2*iband  +jband)=scprod(2,iband)-eshiftkq*doti
 471        indx_cgq=indx_cgq+npw1*nspinor
 472      end do
 473    else
 474      do iband=1,nband
 475        eig1_k(2*iband-1+jband)=scprod(1,iband)
 476        eig1_k(2*iband  +jband)=scprod(2,iband)
 477      end do
 478    end if
 479 
 480   !No more need of gs1c
 481    ABI_DEALLOCATE(gs1c)
 482 
 483  end if
 484 
 485 !Filter the wavefunctions for large modified kinetic energy (see routine mkkin.f)
 486  do ispinor=1,nspinor
 487    ipws=(ispinor-1)*npw1
 488 !$OMP PARALLEL DO PRIVATE(ipw) SHARED(cwavef,kinpw1,ipws,npw1)
 489    do ipw=1+ipws,npw1+ipws
 490      if(kinpw1(ipw-ipws)>huge(zero)*1.d-11)then
 491        cwavef(1:2,ipw)=zero
 492      end if
 493    end do
 494  end do
 495 
 496 !Apply the orthogonality condition: <C1 k,q|C0 k+q>=0 (NCPP) or <C1 k,q|S0|C0 k+q>=0 (PAW)
 497 !Project out all bands from cwavef, i.e. apply P_c projector on cwavef
 498 !(this is needed when there are some partially or unoccupied states)
 499  call projbd(cgq,cwavef,-1,icgq,igscq,istwf_k,mcgq,mgscq,nband,npw1,nspinor,&
 500 & gscq,scprod,0,tim_projbd,useoverlap,me_g0,comm_fft)
 501 
 502  if(ipert/=natom+10.and.ipert/=natom+11) then
 503 !    If PAW, the orthogonality condition is
 504 !    <C1 k,q|S0|C0 k+q>+1/2<C0 k|S1|C0 k+q>=0
 505    if (usepaw==1.and.usedcwavef>0) then
 506 !$OMP PARALLEL DO
 507      do ipw=1,npw1*nspinor
 508        cwavef(1:2,ipw)=cwavef(1:2,ipw)+dcwavef(1:2,ipw)
 509      end do
 510    end if
 511  else ! In 2nd order case, dcwavef/=0 even in NC, and it is already computed in rf2_init (called in dfpt_vtowfk.F90)
 512    do ipw=1,npw1*nspinor
 513      cwavef(:,ipw)=cwavef(:,ipw)+rf2%dcwavef(:,ipw+shift_band)
 514    end do
 515  end if
 516 
 517 !Treat the case of buffer bands
 518  if(band>max(1,nband-nbdbuf))then
 519    cwavef=zero
 520    ghc   =zero
 521    gvnlc =zero
 522    if (gen_eigenpb) gsc=zero
 523    if (usedcwavef==2) dcwavef=zero
 524    if (usepaw==1) then
 525      call pawcprj_set_zero(cwaveprj)
 526    end if
 527    if (usedcwavef==2) dcwavef=zero
 528 !  A small negative residual will be associated with these
 529    resid=-0.1_dp
 530 !  Number of one-way 3D ffts skipped
 531    nskip=nskip+nline
 532 
 533  else
 534 !  If not a buffer band, perform the optimisation
 535 
 536    ABI_ALLOCATE(conjgr,(2,npw1*nspinor))
 537    ABI_ALLOCATE(direc,(2,npw1*nspinor))
 538    ABI_ALLOCATE(gresid,(2,npw1*nspinor))
 539    ABI_ALLOCATE(cwaveq,(2,npw1*nspinor))
 540    if (usepaw==1) then
 541      ABI_DATATYPE_ALLOCATE(conjgrprj,(natom,nspinor))
 542      call pawcprj_alloc(conjgrprj,0,gs_hamkq%dimcprj)
 543    else
 544      ABI_DATATYPE_ALLOCATE(conjgrprj,(0,0))
 545    end if
 546 
 547    cwaveq(:,:)=cgq(:,1+npw1*nspinor*(band-1)+icgq:npw1*nspinor*band+icgq)
 548    dotgp=one
 549 
 550 !  Here apply H(0) at k+q to input orthogonalized 1st-order wfs
 551    sij_opt=0;if (gen_eigenpb) sij_opt=1
 552    cpopt=-1+usepaw
 553    call getghc(cpopt,cwavef,cwaveprj,ghc,gsc,gs_hamkq,gvnlc,eshift,mpi_enreg,1,&
 554 &   prtvol,sij_opt,tim_getghc,0,select_k=KPRIME_H_KPRIME)
 555 
 556 !  ghc also includes the eigenvalue shift
 557    if (gen_eigenpb) then
 558      call cg_zaxpy(npw1*nspinor,(/-eshift,zero/),gsc,ghc)
 559    else
 560      call cg_zaxpy(npw1*nspinor,(/-eshift,zero/),cwavef,ghc)
 561    end if
 562 
 563 !  Initialize resid, in case of nline==0
 564    resid=zero
 565 
 566 !  ======================================================================
 567 !  ====== BEGIN LOOP FOR A GIVEN BAND: MINIMIZATION ITERATIONS ==========
 568 !  ======================================================================
 569    do iline=1,nline
 570 
 571 !    ======================================================================
 572 !    ================= COMPUTE THE RESIDUAL ===============================
 573 !    ======================================================================
 574 !    Note that gresid (=steepest-descent vector, Eq.(26) of PRB 55, 10337 (1996) [[cite:Gonze1997]])
 575 !    is precomputed to garantee cancellation of errors
 576 !    and allow residuals to reach values as small as 1.0d-24 or better.
 577      if (berryopt== 4.or.berryopt== 6.or.berryopt== 7.or.&
 578 &     berryopt==14.or.berryopt==16.or.berryopt==17) then
 579        if (ipert==natom+2) then
 580          if (opt_gvnl1/=1) gvnl1=zero
 581 !$OMP PARALLEL DO
 582          do ipw=1,npw1*nspinor
 583            gresid(1:2,ipw)=-ghc(1:2,ipw)-gh1c(1:2,ipw)
 584          end do
 585        else
 586 !$OMP PARALLEL DO
 587          do ipw=1,npw1*nspinor
 588            gresid(1,ipw)=-ghc(1,ipw)-gh1c(1,ipw)+gberry(2,ipw)
 589            gresid(2,ipw)=-ghc(2,ipw)-gh1c(2,ipw)-gberry(1,ipw)
 590          end do
 591        end if
 592      else
 593 !$OMP PARALLEL DO
 594        do ipw=1,npw1*nspinor
 595          gresid(1:2,ipw)=-ghc(1:2,ipw)-gh1c(1:2,ipw)
 596        end do
 597      end if
 598 
 599 !    ======================================================================
 600 !    =========== PROJECT THE STEEPEST DESCENT DIRECTION ===================
 601 !    ========= OVER THE SUBSPACE ORTHOGONAL TO OTHER BANDS ================
 602 !    ======================================================================
 603 !    Project all bands from gresid into direc:
 604 !    The following projection over the subspace orthogonal to occupied bands
 605 !    is not optional in the RF case, unlike the GS case.
 606 !    However, the order of operations could be changed, so that
 607 !    as to make it only applied at the beginning, to H(1) psi(0),
 608 !    so, THIS IS TO BE REEXAMINED
 609 !    Note the subtlety:
 610 !    -For the generalized eigenPb, S|cgq> is used in place of |cgq>,
 611 !    in order to apply P_c+ projector (see PRB 73, 235101 (2006) [[cite:Audouze2006]], Eq. (71), (72)
 612      if (gen_eigenpb) then
 613        call projbd(gscq,gresid,-1,igscq,icgq,istwf_k,mgscq,mcgq,nband,npw1,nspinor,&
 614 &       cgq,scprod,0,tim_projbd,useoverlap,me_g0,comm_fft)
 615      else
 616        call projbd(cgq,gresid,-1,icgq,0,istwf_k,mcgq,mgscq,nband,npw1,nspinor,&
 617 &       dummy,scprod,0,tim_projbd,useoverlap,me_g0,comm_fft)
 618      end if
 619 
 620      call cg_zcopy(npw1*nspinor,gresid,direc)
 621 
 622 !    ======================================================================
 623 !    ============== CHECK FOR CONVERGENCE CRITERIA ========================
 624 !    ======================================================================
 625 
 626 !    Compute second-order derivative of the energy using a variational expression
 627      call dotprod_g(prod1,doti,istwf_k,npw1*nspinor,1,cwavef,gresid,me_g0,mpi_enreg%comm_spinorfft)
 628      call dotprod_g(prod2,doti,istwf_k,npw1*nspinor,1,cwavef,gh1c,me_g0,mpi_enreg%comm_spinorfft)
 629      d2te=two*(-prod1+prod2)
 630 
 631 !    DEBUG Keep this debugging feature !
 632 !     write(msg,'(a,f14.6,a,f14.6)') 'prod1 = ',prod1,' prod2 = ',prod2
 633 !     call wrtout(std_out,msg,'PERS')
 634 !    ENDDEBUG
 635 
 636 !    Compute <u_m(1)|H(0)-e_m(0)|u_m(1)>
 637 !    (<u_m(1)|H(0)-e_m(0).S|u_m(1)> if gen. eigenPb),
 638 !    that should be positive,
 639 !    except when the eigenvalue eig_mk(0) is higher than
 640 !    the lowest non-treated eig_mk+q(0). For insulators, this
 641 !    has no influence, but for metallic occupations,
 642 !    the conjugate gradient algorithm breaks down. The solution adopted here
 643 !    is very crude, and rely upon the fact that occupancies of such
 644 !    levels should be smaller and smaller with increasing nband, so that
 645 !    a convergence study will give the right result.
 646 !    The same trick is also used later.
 647      u1h0me0u1=-prod1-prod2
 648 !    Some tolerance is allowed, to account for very small numerical inaccuracies and cancellations.
 649      if(u1h0me0u1<-tol_restart)then
 650        if (prtvol==-level.or.prtvol==-19) then
 651          write(msg,'(a,es22.13e3)') '  cgwf3: u1h0me0u1 = ',u1h0me0u1
 652          call wrtout(std_out,msg,'PERS')
 653        end if
 654        cwavef =zero
 655        ghc    =zero
 656        gvnlc  =zero
 657        if (gen_eigenpb) gsc(:,:)=zero
 658        if (usepaw==1) then
 659          call pawcprj_set_zero(cwaveprj)
 660        end if
 661 !      A negative residual will be the signal of this problem ...
 662        resid=-one
 663        call wrtout(std_out,' dfpt_cgwf: problem of minimisation (likely metallic), set resid to -1','PERS')
 664 !      Number of one-way 3D ffts skipped
 665        nskip=nskip+(nline-iline+1)
 666 !      Exit from the loop on iline
 667        exit
 668      end if
 669 
 670 !    Compute residual (squared) norm
 671      call sqnorm_g(resid,istwf_k,npw1*nspinor,gresid,me_g0,comm_fft)
 672      if (prtvol==-level.or.prtvol==-19)then
 673        write(msg,'(a,a,i3,f14.6,a,a,4es12.4)') ch10,&
 674 &       ' dfpt_cgwf : iline,eshift     =',iline,eshift,ch10,&
 675 &       '         resid,prod1,prod2,d2te=',resid,prod1,prod2,d2te
 676        call wrtout(std_out,msg,'PERS')
 677      end if
 678 
 679 !    If residual sufficiently small stop line minimizations
 680      if (resid<tolwfr) then
 681        if(prtvol>=10)then
 682          write(msg,'(a,i4,a,i2,a,es12.4)')' dfpt_cgwf: band',band,' converged after ',iline,' line minimizations: resid = ',resid
 683          call wrtout(std_out,msg,'PERS')
 684        end if
 685        nskip=nskip+(nline-iline+1)  ! Number of two-way 3D ffts skipped
 686        exit                         ! Exit from the loop on iline
 687      end if
 688 
 689 !    If user require exiting the job, stop line minimisations
 690      if (quit==1) then
 691        write(msg,'(a,i0)')' dfpt_cgwf: user require exiting => skip update of band ',band
 692        call wrtout(std_out,msg,'PERS')
 693        nskip=nskip+(nline-iline+1)  ! Number of two-way 3D ffts skipped
 694        exit                         ! Exit from the loop on iline
 695      end if
 696 
 697 !    Check that d2te is decreasing on succeeding lines:
 698      if (iline/=1) then
 699        !if (d2te>d2teold+tol12) then
 700        if (d2te>d2teold+tol6) then
 701          write(msg,'(a,i0,a,e14.6,a,e14.6)')'New trial energy at line ',iline,'=',d2te,'is higher than former:',d2teold
 702          MSG_WARNING(msg)
 703        end if
 704      end if
 705      d2teold=d2te
 706 
 707 !    DEBUG Keep this debugging feature !
 708 !    call sqnorm_g(dotr,istwf_k,npw1*nspinor,direc,me_g0,comm_fft)
 709 !    write(std_out,*)' dfpt_cgwf : before precon, direc**2=',dotr
 710 !    if (gen_eigenpb) then
 711 !    call dotprod_g(dotr,doti,istwf_k,npw1*nspinor,1,cwaveq,&
 712 !    &                 gscq(:,1+npw1*nspinor*(band-1)+igscq:npw1*nspinor*band+igscq),me_g0,mpi_enreg%comm_spinorfft)
 713 !    else
 714 !    call sqnorm_g(dotr,istwf_k,npw1*nspinor,cwaveq,me_g0,comm_fft)
 715 !    end if
 716 !    write(std_out,*)' dfpt_cgwf : before precon, cwaveq**2=',dotr
 717 !    ENDDEBUG
 718 
 719 !    ======================================================================
 720 !    ======== PRECONDITION THE STEEPEST DESCENT DIRECTION =================
 721 !    ======================================================================
 722 
 723 !    If wfoptalg>=10, the precondition matrix is kept constant
 724 !    during iteration ; otherwise it is recomputed
 725      if (wfoptalg<10.or.iline==1) then
 726        call cg_precon(cwaveq,zero,istwf_k,kinpw1,npw1,nspinor,me_g0,0,pcon,direc,mpi_enreg%comm_fft)
 727      else
 728        do ispinor=1,nspinor
 729          igs=(ispinor-1)*npw1
 730 !$OMP PARALLEL DO PRIVATE(ipw) SHARED(igs,npw,direc,pcon)
 731          do ipw=1+igs,npw1+igs
 732            direc(1:2,ipw)=direc(1:2,ipw)*pcon(ipw-igs)
 733          end do
 734        end do
 735      end if
 736 
 737 !    DEBUG Keep this debugging feature !
 738 !    call sqnorm_g(dotr,istwf_k,npw1*nspinor,direc,me_g0,comm_fft)
 739 !    write(std_out,*)' dfpt_cgwf : after precon, direc**2=',dotr
 740 !    ENDDEBUG
 741 
 742 !    ======================================================================
 743 !    ======= PROJECT THE PRECOND. STEEPEST DESCENT DIRECTION ==============
 744 !    ========= OVER THE SUBSPACE ORTHOGONAL TO OTHER BANDS ================
 745 !    ======================================================================
 746 
 747 !    Projecting again out all bands:
 748 !    -For the simple eigenPb, gscq is used as dummy argument
 749      call projbd(cgq,direc,-1,icgq,igscq,istwf_k,mcgq,mgscq,nband,npw1,nspinor,&
 750 &     gscq,scprod,0,tim_projbd,useoverlap,me_g0,comm_fft)
 751 
 752 !    DEBUG Keep this debugging feature !
 753 !    call sqnorm_g(dotr,istwf_k,npw1*nspinor,direc,me_g0,comm_fft)
 754 !    write(std_out,*)' dfpt_cgwf : after projbd, direc**2=',dotr
 755 !    ENDDEBUG
 756 
 757 !    ======================================================================
 758 !    ================= COMPUTE THE CONJUGATE-GRADIENT =====================
 759 !    ======================================================================
 760 
 761 !    get dot of direction vector with residual vector
 762      call dotprod_g(dotgg,doti,istwf_k,npw1*nspinor,1,direc,gresid,me_g0,mpi_enreg%comm_spinorfft)
 763 
 764 !    At first iteration, gamma is set to zero
 765      if (iline==1) then
 766        gamma=zero
 767        dotgp=dotgg
 768        call cg_zcopy(npw1*nspinor,direc,conjgr)
 769      else
 770 !      At next iterations, h = g + gamma * h
 771        gamma=dotgg/dotgp
 772        dotgp=dotgg
 773        if (prtvol==-level.or.prtvol==-19)then
 774          write(msg,'(a,2es16.6)') 'dfpt_cgwf: dotgg,gamma = ',dotgg,gamma
 775          call wrtout(std_out,msg,'PERS')
 776        end if
 777 !$OMP PARALLEL DO
 778        do ipw=1,npw1*nspinor
 779          conjgr(1:2,ipw)=direc(1:2,ipw)+gamma*conjgr(1:2,ipw)
 780        end do
 781        if (prtvol==-level.or.prtvol==-19)then
 782          call wrtout(std_out,'dfpt_cgwf: conjugate direction has been found','PERS')
 783        end if
 784      end if
 785 
 786 !    ======================================================================
 787 !    ===== COMPUTE CONTRIBUTIONS TO 1ST AND 2ND DERIVATIVES OF ENERGY =====
 788 !    ======================================================================
 789 !    ...along the search direction
 790 
 791 !    Compute dedt, Eq.(29) of of PRB55, 10337 (1997) [[cite:Gonze1997]],
 792 !    with an additional factor of 2 for the difference
 793 !    between E(2) and the 2DTE
 794      call dotprod_g(dedt,doti,istwf_k,npw1*nspinor,1,conjgr,gresid,me_g0,mpi_enreg%comm_spinorfft)
 795      dedt=-two*two*dedt
 796      if((prtvol==-level.or.prtvol==-19.or.prtvol==-20).and.dedt-tol14>0) then
 797        write(msg,'(a)') '  CGWF3_WARNING : dedt>0'
 798        call wrtout(std_out,msg,'PERS')
 799      end if
 800      ABI_ALLOCATE(gvnl_direc,(2,npw1*nspinor))
 801      ABI_ALLOCATE(gh_direc,(2,npw1*nspinor))
 802      if (gen_eigenpb)  then
 803        ABI_ALLOCATE(sconjgr,(2,npw1*nspinor))
 804      else
 805        ABI_ALLOCATE(sconjgr,(0,0))
 806      end if
 807      sij_opt=0;if (gen_eigenpb) sij_opt=1
 808      cpopt=-1+usepaw
 809      call getghc(cpopt,conjgr,conjgrprj,gh_direc,sconjgr,gs_hamkq,gvnl_direc,&
 810 &     eshift,mpi_enreg,1,prtvol,sij_opt,tim_getghc,0,select_k=KPRIME_H_KPRIME)
 811 
 812 !    ghc also includes the eigenvalue shift
 813      if (gen_eigenpb) then
 814 !$OMP PARALLEL DO
 815        do ipw=1,npw1*nspinor
 816          gh_direc(1:2,ipw)=gh_direc(1:2,ipw)-eshift*sconjgr(1:2,ipw)
 817        end do
 818      else
 819 !$OMP PARALLEL DO
 820        do ipw=1,npw1*nspinor
 821          gh_direc(1:2,ipw)=gh_direc(1:2,ipw)-eshift*conjgr(1:2,ipw)
 822        end do
 823      end if
 824 
 825 !    compute d2edt2, Eq.(30) of of PRB55, 10337 (1997) [[cite:Gonze1997]],
 826 !    with an additional factor of 2 for the difference
 827 !    between E(2) and the 2DTE, and neglect of local fields (SC terms)
 828      call dotprod_g(d2edt2,doti,istwf_k,npw1*nspinor,1,conjgr,gh_direc,me_g0,mpi_enreg%comm_spinorfft)
 829      d2edt2=two*two*d2edt2
 830      if(prtvol==-level.or.prtvol==-19)then
 831        write(msg,'(a,2es14.6)') 'dfpt_cgwf: dedt,d2edt2=',dedt,d2edt2
 832        call wrtout(std_out,msg,'PERS')
 833      end if
 834 
 835 !    ======================================================================
 836 !    ======= COMPUTE MIXING FACTOR - CHECK FOR CONVERGENCE ===============
 837 !    ======================================================================
 838 
 839 !    see Eq.(31) of PRB55, 10337 (1997) [[cite:Gonze1997]]
 840 !
 841      if(d2edt2<-tol_restart)then
 842 !      This may happen when the eigenvalue eig_mk(0) is higher than
 843 !      the lowest non-treated eig_mk+q(0). The solution adopted here
 844 !      is very crude, and rely upon the fact that occupancies of such
 845 !      levels should be smaller and smaller with increasing nband, so that
 846 !      a convergence study will give the right result.
 847        theta=zero
 848 
 849        cwavef=zero
 850        ghc   =zero
 851        gvnlc =zero
 852        if (gen_eigenpb) gsc=zero
 853        if (usepaw==1) then
 854          call pawcprj_set_zero(cwaveprj)
 855        end if
 856 !      A negative residual will be the signal of this problem ...
 857        resid=-two
 858        call wrtout(std_out,' dfpt_cgwf: problem of minimisation (likely metallic), set resid to -2',"PERS")
 859      else
 860 !      Here, the value of theta that gives the minimum
 861        theta=-dedt/d2edt2
 862 !      DEBUG
 863 !      write(std_out,*)' dfpt_cgwf: dedt,d2edt2=',dedt,d2edt2
 864 !      ENDDEBUG
 865      end if
 866 
 867 !    Check that result is above machine precision
 868      if (one+theta==one) then
 869        write(msg, '(a,es16.4)' ) ' dfpt_cgwf: converged with theta=',theta
 870        call wrtout(std_out,msg,'PERS')
 871        nskip=nskip+2*(nline-iline) ! Number of one-way 3D ffts skipped
 872        exit                        ! Exit from the loop on iline
 873      end if
 874 
 875 !    ======================================================================
 876 !    ================ GENERATE NEW |wf>, H|wf>, Vnl|Wf ... ================
 877 !    ======================================================================
 878 
 879      call cg_zaxpy(npw1*nspinor,(/theta,zero/),conjgr,cwavef)
 880      call cg_zaxpy(npw1*nspinor,(/theta,zero/),gh_direc,ghc)
 881      call cg_zaxpy(npw1*nspinor,(/theta,zero/),gvnl_direc,gvnlc)
 882 
 883      if (gen_eigenpb) then
 884        call cg_zaxpy(npw1*nspinor,(/theta,zero/),sconjgr,gsc)
 885      end if
 886      if (usepaw==1) then
 887        call pawcprj_axpby(theta,one,conjgrprj,cwaveprj)
 888      end if
 889 
 890      ABI_DEALLOCATE(gh_direc)
 891      ABI_DEALLOCATE(gvnl_direc)
 892      ABI_DEALLOCATE(sconjgr)
 893 
 894 !    ======================================================================
 895 !    =========== CHECK CONVERGENCE AGAINST TRIAL ENERGY ===================
 896 !    ======================================================================
 897 
 898      if(usetolrde/=0) then
 899 !      Check reduction in trial energy deltae, Eq.(28) of PRB55, 10337 (1997) [[cite:Gonze1997]]
 900        deltae=half*d2edt2*theta**2+theta*dedt
 901 
 902        if (iline==1) then
 903          deold=deltae
 904 !        The extra factor of two should be removed !
 905        else if (abs(deltae)<tolrde*two*abs(deold) .and. iline/=nline) then
 906          if(prtvol>=10.or.prtvol==-level.or.prtvol==-19)then
 907            write(msg, '(a,i4,1x,a,1p,e12.4,a,e12.4,a)' ) &
 908 &           ' dfpt_cgwf: line',iline,' deltae=',deltae,' < tolrde*',deold,' =>skip lines'
 909            call wrtout(std_out,msg,'PERS')
 910          end if
 911          nskip=nskip+2*(nline-iline) ! Number of one-way 3D ffts skipped
 912          exit                        ! Exit from the loop on iline
 913        end if
 914      end if
 915 
 916 !    ======================================================================
 917 !    ================== END LOOP FOR GIVEN BAND ===========================
 918 !    ======================================================================
 919 
 920 !    Note that there are three "exit" instruction inside the loop.
 921      nlines_done = nlines_done + 1
 922    end do ! iline
 923 
 924 !  Check that final cwavef (Psi^(1)) satisfies the orthogonality condition
 925    if (prtvol==-level.or.prtvol==-19) then
 926      sij_opt=0 ; usevnl=1 ; optlocal=1 ; optnl=2 ; if (gen_eigenpb)  sij_opt=1
 927      ABI_ALLOCATE(work,(2,npw1*nspinor))
 928      ABI_ALLOCATE(work1,(2,npw1*nspinor))
 929      ABI_ALLOCATE(work2,(2,npw1*nspinor*sij_opt))
 930      do iband=1,nband
 931        if (gen_eigenpb) then
 932          work(:,:)=gscq(:,1+npw1*nspinor*(iband-1)+igscq:npw1*nspinor*iband+igscq)
 933        else
 934          work(:,:)=cgq(:,1+npw1*nspinor*(iband-1)+icgq:npw1*nspinor*iband+icgq)
 935        end if
 936        ! Compute : <Psi^(0)_i,k+q|Psi^(1)_j,k,q>
 937        call dotprod_g(prod1,prod2,istwf_k,npw1*nspinor,2,work,cwavef,me_g0,mpi_enreg%comm_spinorfft)
 938 
 939        if (ipert/=natom+10.and.ipert/=natom+11) then
 940          if (gen_eigenpb) then
 941            call getgh1c(berryopt,cwave0,cwaveprj0,work1,gberry,work2,gs_hamkq,gvnl1_saved,idir,ipert,eshift,&
 942 &           mpi_enreg,optlocal,optnl,opt_gvnl1,rf_hamkq,sij_opt,tim_getgh1c,usevnl)
 943            work(:,:)=cgq(:,1+npw1*nspinor*(iband-1)+icgq:npw1*nspinor*iband+icgq)
 944            call dotprod_g(dotr,doti,istwf_k,npw1*nspinor,2,work,work2,me_g0,mpi_enreg%comm_spinorfft)
 945          else
 946            dotr=zero; doti=zero
 947          end if
 948          dotr=prod1+half*dotr
 949          doti=prod2+half*doti
 950        else if(prtvol==-19) then ! 2nd order case
 951          dotr=prod1+half*rf2%amn(1,iband+(band-1)*nband)
 952          doti=prod2+half*rf2%amn(2,iband+(band-1)*nband)
 953        else
 954          write(msg,'(a)') 'CGWF3_WARNING : Use prtvol=-19 to test orthogonality for ipert=natom+10 or +11'
 955          call wrtout(std_out,msg,'COLL')
 956        end if
 957        dotr=sqrt(dotr**2+doti**2)
 958        if(dotr>tol10) then
 959 !         if (gen_eigenpb) then
 960 !           write(msg,'(2a,i3,a,2es22.15)') 'CGWF3_WARNING : <Psi^(1)_i,k,q|S^(0)|Psi^(0)_j,k+q>',&
 961 !&           '+ 1/2<Psi^(0)_i,k|S^(1)|Psi^(0)_j,k+q>, for j= ',iband,' is ',dotr,doti
 962 !         else
 963          write(msg,'(a,i3,a,es22.15)') 'CGWF3_WARNING : |<Psi^(0)_i,k+q|Psi^(1)_j,k,q>+amn(i,j)/2|, for j= ',iband,' is ',dotr
 964 !         end if
 965          call wrtout(std_out,msg,'PERS')
 966        end if
 967      end do
 968      ABI_DEALLOCATE(work)
 969      ABI_DEALLOCATE(work1)
 970      ABI_DEALLOCATE(work2)
 971      if (ipert/=natom+10.and.ipert/=natom+11) then
 972        ABI_DEALLOCATE(gvnl1_saved)
 973      end if
 974    end if
 975 
 976 !  Check that final cwavef Psi^(1) is Pc.Psi^(1)+delta_Psi^(1)
 977    if (prtvol==-level.or.prtvol==-19)then
 978      ABI_ALLOCATE(cwwork,(2,npw1*nspinor))
 979      cwwork=cwavef
 980 !    -Apply Pc to Psi^(1)
 981      if (gen_eigenpb) then
 982        call projbd(cgq,cwwork,-1,icgq,igscq,istwf_k,mcgq,mgscq,nband,npw1,nspinor,&
 983 &       gscq,scprod,0,tim_projbd,useoverlap,me_g0,comm_fft)
 984      else
 985        call projbd(cgq,cwwork,-1,icgq,0,istwf_k,mcgq,mgscq,nband,npw1,nspinor,&
 986 &       dummy,scprod,0,tim_projbd,useoverlap,me_g0,comm_fft)
 987      end if
 988 !    -Add delta_Psi^(1)
 989      if (usedcwavef>0) cwwork=cwwork+dcwavef
 990      if(ipert==natom+10.or.ipert==natom+11) cwwork=cwwork+rf2%dcwavef(:,1+shift_band:npw1*nspinor+shift_band)
 991 !    -Compare to Psi^(1)
 992      cwwork=cwwork-cwavef
 993      call sqnorm_g(dotr,istwf_k,npw1*nspinor,cwwork,me_g0,comm_fft)
 994      ABI_DEALLOCATE(cwwork)
 995      if(sqrt(dotr)>tol10) then
 996 !       if (gen_eigenpb) then
 997 !         write(msg,'(a,i3,a,es22.15)') &
 998 !&         'CGWF3_WARNING : |(Pc.Psi^(1)_i,k,q + delta_Psi^(1)_i,k) - Psi^(1)_i,k,q|^2 (band ',band,')=',dotr
 999 !       else
1000        write(msg,'(a,es22.15)') 'CGWF3_WARNING : |(Pc.Psi^(1)_i,k,q + delta_Psi^(1)_i,k) - Psi^(1)_i,k,q| = ',sqrt(dotr)
1001 !       end if
1002        call wrtout(std_out,msg,'PERS')
1003      end if
1004    end if
1005 
1006 !  Check that final cwavef (Psi^(1)) solves the Sternheimer equation
1007    if(prtvol==-level.or.prtvol==-19.or.prtvol==-20)then
1008      ABI_ALLOCATE(cwwork,(2,npw1*nspinor))
1009      cwwork=cwavef
1010 !    -Apply Pc to Psi^(1)
1011      if (gen_eigenpb) then
1012        call projbd(cgq,cwwork,-1,icgq,igscq,istwf_k,mcgq,mgscq,nband,npw1,nspinor,&
1013 &       gscq,scprod,0,tim_projbd,useoverlap,me_g0,comm_fft)
1014      else
1015        call projbd(cgq,cwwork,-1,icgq,0,istwf_k,mcgq,mgscq,nband,npw1,nspinor,&
1016 &       dummy,scprod,0,tim_projbd,useoverlap,me_g0,comm_fft)
1017      end if
1018 !    - Apply H^(0)-E.S^(0)
1019      sij_opt=0;if (gen_eigenpb) sij_opt=1
1020      cpopt=-1
1021      ABI_ALLOCATE(work,(2,npw1*nspinor))
1022      ABI_ALLOCATE(work1,(2,npw1*nspinor*((sij_opt+1)/2)))
1023      ABI_ALLOCATE(work2,(2,npw1*nspinor))
1024      call getghc(cpopt,cwwork,conjgrprj,work,work1,gs_hamkq,work2,eshift,&
1025 &     mpi_enreg,1,prtvol,sij_opt,tim_getghc,0,select_k=KPRIME_H_KPRIME)
1026      if (gen_eigenpb) then
1027        cwwork=work-eshift*work1
1028      else
1029        cwwork=work-eshift*cwwork
1030      end if
1031      ABI_DEALLOCATE(work)
1032      ABI_DEALLOCATE(work1)
1033      ABI_DEALLOCATE(work2)
1034 !  The following is not mandatory, as Pc has been already applied to Psi^(1)
1035 !  and Pc^* H^(0) Pc = Pc^* H^(0) = H^(0) Pc (same for S^(0)).
1036 !  However, in PAW, to apply Pc^* here seems to reduce the numerical error
1037 !    -Apply Pc^*
1038      if (gen_eigenpb) then
1039        call projbd(gscq,cwwork,-1,igscq,icgq,istwf_k,mgscq,mcgq,nband,npw1,nspinor,&
1040 &       cgq,scprod,0,tim_projbd,useoverlap,me_g0,comm_fft)
1041      else
1042        call projbd(cgq,cwwork,-1,icgq,0,istwf_k,mcgq,mgscq,nband,npw1,nspinor,&
1043 &       dummy,scprod,0,tim_projbd,useoverlap,me_g0,comm_fft)
1044      end if
1045 !    - Add Pc^*(H^(1)-E.S^(1)).Psi^(0)
1046      cwwork=cwwork+gh1c
1047      call sqnorm_g(dotr,istwf_k,npw1*nspinor,cwwork,me_g0,comm_fft)
1048      ABI_DEALLOCATE(cwwork)
1049      write(msg,'(a,i3,a,es22.15,2a,i4)') &
1050 &     '*** CGWF3 Sternheimer equation test for band ',band,'=',sqrt(dotr),ch10,&
1051 &     'It should go to zero for large nline : nlines_done = ',nlines_done
1052      call wrtout(std_out,msg,'PERS')
1053    end if
1054 
1055 !  Check that < Psi^(0) | ( H^(0)-eps^(0) S^(0) ) | Psi^(1) > is in agreement with eig^(1)
1056    if(prtvol==-level.or.prtvol==-19.or.prtvol==-20)then
1057      ABI_ALLOCATE(cwwork,(2,npw1*nspinor))
1058      cwwork=cwavef
1059 !    - Apply H^(0)-E.S^(0)
1060      sij_opt=0;if (gen_eigenpb) sij_opt=1
1061      cpopt=-1
1062      ABI_ALLOCATE(work,(2,npw1*nspinor))
1063      ABI_ALLOCATE(work1,(2,npw1*nspinor*((sij_opt+1)/2)))
1064      ABI_ALLOCATE(work2,(2,npw1*nspinor))
1065      call getghc(cpopt,cwwork,conjgrprj,work,work1,gs_hamkq,work2,eshift,&
1066 &     mpi_enreg,1,prtvol,sij_opt,tim_getghc,0,select_k=KPRIME_H_KPRIME)
1067      if (gen_eigenpb) then
1068        cwwork=work-eshift*work1
1069      else
1070        cwwork=work-eshift*cwwork
1071      end if
1072      ABI_DEALLOCATE(work1)
1073      ABI_DEALLOCATE(work2)
1074      cwwork=cwwork+gh1c_n
1075      jband=(band-1)*2*nband
1076      do iband=1,nband
1077        work(:,:)=cgq(:,1+npw1*nspinor*(iband-1)+icgq:npw1*nspinor*iband+icgq)
1078        call dotprod_g(dotr,doti,istwf_k,npw1*nspinor,2,work,cwwork,me_g0,mpi_enreg%comm_spinorfft)
1079        dotr = dotr - eig1_k(2*iband-1+jband)
1080        doti = doti - eig1_k(2*iband  +jband)
1081        dotr = sqrt(dotr**2+doti**2)
1082        if (dotr > tol8) then
1083          write(msg,'(2(a,i3),a,es22.15)') &
1084 &         'CGWF3_WARNING < Psi^(0) | ( H^(0)-eps^(0) S^(0) ) | Psi^(1) > for i=',iband,' j=',band,&
1085          ' : ',sqrt(dotr**2+doti**2)
1086          call wrtout(std_out,msg,'PERS')
1087        end if
1088      end do
1089 !     write(msg,'(a)') '< Psi^(0) | ( H^(0)-eps^(0) S^(0) ) | Psi^(1) > is done.'
1090 !     call wrtout(std_out,msg,'PERS')
1091      ABI_DEALLOCATE(work)
1092      ABI_DEALLOCATE(cwwork)
1093    end if
1094 
1095    if (allocated(gh_direc))  then
1096      ABI_DEALLOCATE(gh_direc)
1097    end if
1098    if (allocated(gvnl_direc))  then
1099      ABI_DEALLOCATE(gvnl_direc)
1100    end if
1101    ABI_DEALLOCATE(conjgr)
1102    ABI_DEALLOCATE(cwaveq)
1103    ABI_DEALLOCATE(direc)
1104    ABI_DEALLOCATE(gresid)
1105    if (usepaw==1) then
1106      call pawcprj_free(conjgrprj)
1107    end if
1108    ABI_DATATYPE_DEALLOCATE(conjgrprj)
1109 
1110  end if  ! End condition of not being a buffer band
1111 
1112 !At the end of the treatment of a set of bands, write the number of one-way 3D ffts skipped
1113  if (xmpi_paral==1 .and. band==nband .and. prtvol>=10) then
1114    write(msg,'(a,i0)')' dfpt_cgwf: number of one-way 3D ffts skipped in cgwf3 until now =',nskip
1115    call wrtout(std_out,msg,'PERS')
1116  end if
1117 
1118  ABI_DEALLOCATE(gh1c)
1119  ABI_DEALLOCATE(pcon)
1120  ABI_DEALLOCATE(scprod)
1121  ABI_DEALLOCATE(gberry)
1122 
1123  call timab(122,2,tsec)
1124 
1125  DBG_EXIT("COLL")
1126 
1127 end subroutine dfpt_cgwf

ABINIT/m_dfpt_cgwf [ Modules ]

[ Top ] [ Modules ]

NAME

  m_dfpt_cgwf

FUNCTION

 Update one single wavefunction (cwavef), non self-consistently.
 Uses a conjugate-gradient algorithm.

COPYRIGHT

  Copyright (C) 1999-2018 ABINIT group (XG,DRH,XW,FJ,MT,LB)
  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

22 #if defined HAVE_CONFIG_H
23 #include "config.h"
24 #endif
25 
26 #include "abi_common.h"
27 
28 module m_dfpt_cgwf
29 
30  use defs_basis
31  use defs_abitypes
32  use m_abicore
33  use m_errors
34  use m_xmpi
35  use m_cgtools
36  use m_rf2
37 
38  use m_time,        only : timab
39  use m_pawcprj,     only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_set_zero, pawcprj_axpby
40  use m_hamiltonian, only : gs_hamiltonian_type, rf_hamiltonian_type, KPRIME_H_KPRIME
41  use m_getghc,      only : getghc
42  use m_getgh1c,     only : getgh1c, getdc1
43 
44  implicit none
45 
46  private