TABLE OF CONTENTS


ABINIT/m_cgwf [ Modules ]

[ Top ] [ Modules ]

NAME

  m_cgwf

FUNCTION

  Conjugate-gradient eigensolver.

COPYRIGHT

  Copyright (C) 2008-2018 ABINIT group (DCA, XG, GMR, MT, MVeithen, ISouza, JIniguez)
  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_cgwf
28 
29  use defs_basis
30  use defs_abitypes
31  use m_errors
32  use m_xmpi
33  use m_abicore
34  use m_cgtools
35  use m_efield
36 
37  use m_time,          only : timab
38  use m_numeric_tools, only : rhophi
39  use m_pawcprj,       only : pawcprj_type, pawcprj_alloc, pawcprj_put, pawcprj_copy, &
40                              pawcprj_get, pawcprj_mpi_allgather, pawcprj_free, pawcprj_symkn
41  use m_hamiltonian,   only : gs_hamiltonian_type
42  use m_fock,          only : fock_set_ieigen,fock_set_getghc_call
43  use m_getghc,        only : getghc
44  use m_berrytk,       only : smatrix
45  use m_nonlop,      only : nonlop
46  use m_paw_overlap,   only : smatrix_k_paw
47  use m_cgprj,         only : getcprj
48 
49  implicit none
50 
51  private

ABINIT/make_grad_berry [ Functions ]

[ Top ] [ Functions ]

NAME

 make_grad_berry

FUNCTION

 compute gradient contribution from berry phase in finite
 electric field case

COPYRIGHT

 Copyright (C) 1998-2018 ABINIT group
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

INPUTS

  cg(2,mcg)=input wavefunctions
  cgq(2,mcgq) = wavefunctions at neighboring k points
  cprj_k(natom,nband_k*usepaw)=cprj at this k point
  dimlmn(natom)=lmn_size for each atom in input order
  dimlmn_srt(natom)=lmn_size for each atom sorted by type
  direc(2,npw*nspinor)=gradient vector
  gs_hamk <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k
  iband=index of band currently being treated
  icg=shift to be applied on the location of data in the array cg
  ikpt=number of the k-point currently being treated
  isppol=spin polarization currently treated
  natom=number of atoms in cell.
  mband =maximum number of bands
  mpw=maximum dimensioned size of npw
  mcg=second dimension of the cg array
  mcgq=second dimension of the cgq array
  mkgq = second dimension of pwnsfacq
  nkpt=number of k points
  mpi_enreg=information about MPI parallelization
  npw=number of planewaves in basis sphere at given k.
  nspinor=number of spinorial components of the wavefunctions (on current proc)
  nsppol=number of spin polarizations
  pwind(pwind_alloc,2,3) = array used to compute
           the overlap matrix smat between k-points (see initberry.f)
  pwind_alloc = first dimension of pwind
  pwnsfac(2,pwind_alloc) = phase factors for non-symmorphic translations
                           (see initberry.f)
  pwnsfacq(2,mkgq) = phase factors for the nearest neighbours of the
                     current k-point (electric field, MPI //)

OUTPUT

 grad_berry(2,npw*nspinor) :: contribution to gradient in finite electric field case

SIDE EFFECTS

  dtefield <type(efield_type)> = variables related to Berry phase calculations (see initberry.f)

NOTES

PARENTS

      cgwf

CHILDREN

      nonlop,pawcprj_alloc,pawcprj_copy,pawcprj_free,pawcprj_get
      pawcprj_symkn,smatrix,smatrix_k_paw

SOURCE

2110 subroutine make_grad_berry(cg,cgq,cprj_k,detovc,dimlmn,dimlmn_srt,direc,dtefield,grad_berry,&
2111 &                          gs_hamk,iband,icg,ikpt,isppol,mband,mcg,mcgq,mkgq,mpi_enreg,mpw,natom,nkpt,&
2112 &                          npw,npwarr,nspinor,nsppol,pwind,pwind_alloc,pwnsfac,pwnsfacq)
2113 
2114 
2115 !This section has been created automatically by the script Abilint (TD).
2116 !Do not modify the following lines by hand.
2117 #undef ABI_FUNC
2118 #define ABI_FUNC 'make_grad_berry'
2119 !End of the abilint section
2120 
2121   implicit none
2122 
2123   !Arguments ------------------------------------
2124 
2125   !scalars
2126   integer,intent(in) :: iband,icg,ikpt,isppol,mband,mcg,mcgq
2127   integer,intent(in) :: mkgq,mpw,natom,nkpt,npw,nspinor,nsppol,pwind_alloc
2128   type(gs_hamiltonian_type),intent(in) :: gs_hamk
2129   type(efield_type),intent(inout) :: dtefield
2130   type(MPI_type),intent(in) :: mpi_enreg
2131 
2132   !arrays
2133   integer,intent(in) :: dimlmn(natom),dimlmn_srt(natom)
2134   integer,intent(in) :: npwarr(nkpt),pwind(pwind_alloc,2,3)
2135   real(dp),intent(in) :: cg(2,mcg),cgq(2,mcgq)
2136   real(dp),intent(inout) :: direc(2,npw*nspinor)
2137   real(dp),intent(in) :: pwnsfac(2,pwind_alloc),pwnsfacq(2,mkgq)
2138   real(dp),intent(out) :: detovc(2,2,3),grad_berry(2,npw*nspinor)
2139   type(pawcprj_type),intent(in) :: cprj_k(natom,dtefield%mband_occ*gs_hamk%usepaw*dtefield%nspinor)
2140 
2141   !Local variables-------------------------------
2142   !scalars
2143   integer :: choice,cpopt,ddkflag,dimenlr1,iatom,icg1,icp2,idum1
2144   integer :: idir,ifor,ikgf,ikptf,ikpt2,ikpt2f,ipw,i_paw_band,ispinor,itrs,itypat,job
2145   integer :: klmn,mcg1_k,mcg_q,nbo,npw_k2,nspinortot,paw_opt,shiftbd,signs
2146   real(dp) :: fac
2147   character(len=500) :: message
2148   !arrays
2149   integer :: pwind_k(npw),sflag_k(dtefield%mband_occ)
2150   real(dp) :: cg1_k(2,npw*nspinor),dtm_k(2),pwnsfac_k(4,mpw)
2151   real(dp) :: smat_k(2,dtefield%mband_occ,dtefield%mband_occ)
2152   real(dp) :: smat_inv(2,dtefield%mband_occ,dtefield%mband_occ),svectout_dum(2,0)
2153   real(dp) :: dummy_enlout(0)
2154   real(dp),allocatable :: cgq_k(:,:),enl_rij(:,:,:,:),grad_berry_ev(:,:)
2155   real(dp),allocatable :: qijbkk(:,:,:,:),smat_k_paw(:,:,:)
2156   ! type(pawcprj_type) :: cprj_dum(1,1) ! was used in on-site dipole, now suppressed
2157   ! 15 June 2012 J Zwanziger
2158   type(pawcprj_type),allocatable :: cprj_kb(:,:),cprj_band_srt(:,:)
2159   type(pawcprj_type),allocatable :: cprj_fkn(:,:),cprj_ikn(:,:)
2160 
2161 
2162   ! *********************************************************************
2163 
2164   !DBG_ENTER("COLL")
2165 
2166   nbo = dtefield%mband_occ
2167 
2168   !allocations
2169 
2170   !Electric field: compute the gradient of the Berry phase part of the energy functional.
2171   !See PRL 89, 117602 (2002) [[cite:Souza2002]], grad_berry(:,:) is the second term of Eq. (4)
2172   grad_berry(:,:) = zero
2173   job = 11 ; shiftbd = 1
2174   mcg_q = mpw*mband*nspinor
2175   mcg1_k = npw*nspinor
2176 
2177   if (gs_hamk%usepaw /= 0) then
2178      dimenlr1 = gs_hamk%lmnmax*(gs_hamk%lmnmax+1)/2
2179      ABI_ALLOCATE(qijbkk,(dimenlr1,natom,nspinor**2,2))
2180      ABI_ALLOCATE(enl_rij,(nspinor*dimenlr1,natom,nspinor**2,1))
2181      ABI_ALLOCATE(smat_k_paw,(2,nbo,nbo))
2182      ABI_ALLOCATE(grad_berry_ev,(2,npw*nspinor))
2183      enl_rij = zero
2184      qijbkk = zero
2185      smat_k_paw = zero
2186      ABI_DATATYPE_ALLOCATE(cprj_kb,(natom,nbo*nspinor))
2187      call pawcprj_alloc(cprj_kb,0,dimlmn)
2188      ABI_DATATYPE_ALLOCATE(cprj_band_srt,(natom,nspinor))
2189      call pawcprj_alloc(cprj_band_srt,0,dimlmn_srt)
2190      if (nkpt /= dtefield%fnkpt) then
2191         ABI_DATATYPE_ALLOCATE(cprj_fkn,(natom,nbo*nspinor))
2192         ABI_DATATYPE_ALLOCATE(cprj_ikn,(natom,nbo*nspinor))
2193         call pawcprj_alloc(cprj_fkn,0,dimlmn)
2194         call pawcprj_alloc(cprj_ikn,0,dimlmn)
2195      else
2196         ABI_DATATYPE_ALLOCATE(cprj_fkn,(0,0))
2197         ABI_DATATYPE_ALLOCATE(cprj_ikn,(0,0))
2198      end if
2199   else
2200      ABI_ALLOCATE(qijbkk,(0,0,0,0))
2201      ABI_ALLOCATE(enl_rij,(0,0,0,0))
2202      ABI_ALLOCATE(smat_k_paw,(0,0,0))
2203      ABI_ALLOCATE(grad_berry_ev,(0,0))
2204      ABI_DATATYPE_ALLOCATE(cprj_kb,(0,0))
2205      ABI_DATATYPE_ALLOCATE(cprj_band_srt,(0,0))
2206      ABI_DATATYPE_ALLOCATE(cprj_fkn,(0,0))
2207      ABI_DATATYPE_ALLOCATE(cprj_ikn,(0,0))
2208   end if
2209 
2210   ikptf = dtefield%i2fbz(ikpt)
2211   ikgf = dtefield%fkgindex(ikptf)  ! this is the shift for pwind
2212 
2213   do idir = 1, 3
2214      !  skip idir values for which efield_dot(idir)=0
2215      if (abs(dtefield%efield_dot(idir)) < tol12) cycle
2216      !  Implicitly, we use the gradient multiplied by the number of k points in the FBZ
2217      fac = dtefield%efield_dot(idir)*dble(dtefield%fnkpt)/&
2218           &   (dble(dtefield%nstr(idir))*four_pi)
2219      do ifor = 1, 2
2220         !    Handle dtefield%i2fbz properly and ask whether t.r.s. is used
2221         ikpt2f = dtefield%ikpt_dk(ikptf,ifor,idir)
2222         if (dtefield%indkk_f2ibz(ikpt2f,6) == 1) then
2223            itrs = 10
2224         else
2225            itrs = 0
2226         end if
2227         ikpt2 = dtefield%indkk_f2ibz(ikpt2f,1)
2228         npw_k2 = npwarr(ikpt2)
2229         ABI_ALLOCATE(cgq_k,(2,nbo*nspinor*npw_k2))
2230         pwind_k(1:npw) = pwind(ikgf+1:ikgf+npw,ifor,idir)
2231         pwnsfac_k(1:2,1:npw) = pwnsfac(1:2,ikgf+1:ikgf+npw)
2232         sflag_k(:) = dtefield%sflag(:,ikpt+(isppol-1)*nkpt,ifor,idir)
2233         smat_k(:,:,:) = dtefield%smat(:,:,:,ikpt+(isppol-1)*nkpt,ifor,idir)
2234         if (mpi_enreg%nproc_cell > 1) then
2235            icg1 = dtefield%cgqindex(2,ifor+2*(idir-1),ikpt+(isppol-1)*nkpt)
2236            cgq_k(:,1:nbo*nspinor*npw_k2) = &
2237                 &       cgq(:,icg1+1:icg1+nbo*nspinor*npw_k2)
2238            idum1 = dtefield%cgqindex(3,ifor+2*(idir-1),ikpt+(isppol-1)*nkpt)
2239            pwnsfac_k(3:4,1:npw_k2) = pwnsfacq(1:2,idum1+1:idum1+npw_k2)
2240         else
2241            icg1 = dtefield%cgindex(ikpt2,isppol)
2242            cgq_k(:,1:nbo*nspinor*npw_k2) = &
2243                 &       cg(:,icg1+1:icg1+nbo*nspinor*npw_k2)
2244            idum1 = dtefield%fkgindex(ikpt2f)
2245            pwnsfac_k(3:4,1:npw_k2) = pwnsfac(1:2,idum1+1:idum1+npw_k2)
2246         end if
2247         if (gs_hamk%usepaw == 1) then
2248            icp2=nbo*(ikpt2-1)*nspinor
2249            call pawcprj_get(gs_hamk%atindx1,cprj_kb,dtefield%cprj,natom,1,icp2,ikpt,0,isppol,&
2250                 &       nbo,dtefield%fnkpt,natom,nbo,nbo,nspinor,nsppol,0,&
2251                 &       mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
2252            if (ikpt2 /= ikpt2f) then ! construct cprj_kb by symmetry
2253               call pawcprj_copy(cprj_kb,cprj_ikn)
2254               call pawcprj_symkn(cprj_fkn,cprj_ikn,dtefield%atom_indsym,dimlmn,-1,gs_hamk%indlmn,&
2255                    &         dtefield%indkk_f2ibz(ikpt2f,2),dtefield%indkk_f2ibz(ikpt2f,6),&
2256                    &         dtefield%fkptns(:,dtefield%i2fbz(ikpt2)),&
2257                    &         dtefield%lmax,dtefield%lmnmax,mband,natom,nbo,nspinor,&
2258                    &         dtefield%nsym,gs_hamk%ntypat,gs_hamk%typat,dtefield%zarot)
2259               call pawcprj_copy(cprj_fkn,cprj_kb)
2260            end if
2261            call smatrix_k_paw(cprj_k,cprj_kb,dtefield,idir,ifor,mband,natom,smat_k_paw,gs_hamk%typat)
2262         end if
2263 
2264         icg1 = 0 ; ddkflag = 1
2265         call smatrix(cg,cgq_k,cg1_k,ddkflag,dtm_k,icg,icg1,itrs,&
2266              &     job,iband,mcg,mcg_q,mcg1_k,iband,mpw,nbo,dtefield%nband_occ(isppol),&
2267              &     npw,npw_k2,nspinor,pwind_k,pwnsfac_k,sflag_k,&
2268              &     shiftbd,smat_inv,smat_k,smat_k_paw,gs_hamk%usepaw)
2269         ABI_DEALLOCATE(cgq_k)
2270         detovc(:,ifor,idir) = dtm_k(:) !store the determinant of the overlap
2271         if (sqrt(dtm_k(1)*dtm_k(1) + dtm_k(2)*dtm_k(2)) < tol12) then
2272            write(message,'(3a,i5,a,i3,a,a,a)') &
2273                 &       '  (electric field)',ch10,&
2274                 &       '  For k-point #',ikpt,' and band # ',iband,',',ch10,&
2275                 &       '  the determinant of the overlap matrix is found to be 0. Fixing...'
2276            !      REC try this:
2277            write(std_out,*)message,dtm_k(1:2)
2278            if(abs(dtm_k(1))<=1d-12)dtm_k(1)=1d-12
2279            if(abs(dtm_k(2))<=1d-12)dtm_k(2)=1d-12
2280            write(std_out,*)' Changing to:',dtm_k(1:2)
2281            !      REC       MSG_BUG(message)
2282         end if
2283 
2284         if (gs_hamk%usepaw == 1) then
2285            !      this loop applies discretized derivative of projectors
2286            !      note that qijb_kk is sorted by input atom order, but nonlop wants it sorted by type
2287            do iatom = 1, natom
2288               itypat = gs_hamk%typat(gs_hamk%atindx1(iatom))
2289               do klmn = 1, dtefield%lmn2_size(itypat)
2290                  !          note: D_ij-like terms have 4 spinor components: 11, 22, 12, and 21. Here the qijb is diagonal
2291                  !          in spin space so only the first two are nonzero and they are equal
2292                  do ispinor = 1, nspinor
2293                     qijbkk(klmn,iatom,ispinor,1) = dtefield%qijb_kk(1,klmn,gs_hamk%atindx1(iatom),idir)
2294                     qijbkk(klmn,  iatom,ispinor,2) = dtefield%qijb_kk(2,klmn,gs_hamk%atindx1(iatom),idir)
2295                     if (ifor > 1) qijbkk(klmn,iatom,ispinor,2) = -qijbkk(klmn,iatom,ispinor,2)
2296                  end do
2297               end do ! end loop over lmn2_size
2298            end do ! end loop over natom
2299 
2300            choice = 1
2301            signs = 2
2302            paw_opt = 1
2303            cpopt = 2 ! use cprj_kb in memory
2304            nspinortot=min(2,nspinor*(1+mpi_enreg%paral_spinor))
2305            do i_paw_band = 1, nbo
2306 
2307               call pawcprj_get(gs_hamk%atindx,cprj_band_srt,cprj_kb,natom,i_paw_band,0,ikpt,1,&
2308                    &         isppol,nbo,1,natom,1,nbo,nspinor,nsppol,0,&
2309                    &         mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
2310 
2311               ! Pass dummy_enlout to avoid aliasing (enl, enlout)
2312               call nonlop(choice,cpopt,cprj_band_srt,dummy_enlout,gs_hamk,idir,(/zero/),mpi_enreg,1,0,&
2313                    &         paw_opt,signs,svectout_dum,0,direc,grad_berry_ev,enl=qijbkk)
2314 
2315               !        Add i*fac*smat_inv(i_paw_band,iband)*grad_berry_ev to the gradient
2316               do ipw = 1, npw*nspinor
2317 
2318                  grad_berry(1,ipw) = grad_berry(1,ipw) - &
2319                       &           fac*(smat_inv(2,i_paw_band,iband)*grad_berry_ev(1,ipw) + &
2320                       &           smat_inv(1,i_paw_band,iband)*grad_berry_ev(2,ipw))
2321 
2322                  grad_berry(2,ipw) = grad_berry(2,ipw) + &
2323                       &           fac*(smat_inv(1,i_paw_band,iband)*grad_berry_ev(1,ipw) - &
2324                       &           smat_inv(2,i_paw_band,iband)*grad_berry_ev(2,ipw))
2325 
2326               end do
2327            end do
2328         end if ! end if PAW
2329 
2330         !    Add i*fac*cg1_k to the gradient
2331         do ipw = 1, npw*nspinor
2332            grad_berry(1,ipw) = grad_berry(1,ipw) - fac*cg1_k(2,ipw)
2333            grad_berry(2,ipw) = grad_berry(2,ipw) + fac*cg1_k(1,ipw)
2334         end do
2335         fac = -1._dp*fac
2336         dtefield%sflag(:,ikpt+(isppol-1)*nkpt,ifor,idir) = sflag_k(:)
2337         dtefield%sflag(iband,ikpt+(isppol-1)*nkpt,ifor,idir) = 0
2338         dtefield%smat(:,:,:,ikpt+(isppol-1)*nkpt,ifor,idir) = smat_k(:,:,:)
2339      end do  ! ifor
2340 
2341      !  if (gs_hamk%usepaw == 1) then
2342      !  !    call nonlop to apply on-site dipole <EV> part to direc
2343      !  !    note that rij is sorted by input atom order, but nonlop wants it sorted by type
2344      !  do iatom = 1, natom
2345      !  itypat = gs_hamk%typat(gs_hamk%atindx1(iatom))
2346      !  do klmn = 1, dtefield%lmn2_size(itypat)
2347      !  !        note: D_ij-like terms have 4 spinor components: 11, 22, 12, and 21. Here the enl_rij is diagonal
2348      !  !        in spin space so only the first two are nonzero and they are equal
2349      !  do ispinor = 1, nspinor
2350      !  if (nspinor == 1) then
2351      !  enl_rij(klmn,iatom,ispinor) = dtefield%rij(klmn,itypat,idir)
2352      !  else
2353      !  enl_rij(2*klmn-1,iatom,ispinor) = dtefield%rij(klmn,itypat,idir)
2354      !  end if
2355      !  end do
2356      !  end do ! end loop over lmn2_size
2357      !  end do ! end loop over natom
2358      !  cpopt = -1 ! compute cprj inside nonlop because we do not have them for direc
2359      !  call nonlop(choice,cpopt,cprj_dum,dummy_enlout,gs_hamk,idir,zero,mpi_enreg,1,0,&
2360      !  &           paw_opt,signs,svectout_dum,0,direc,grad_berry_ev,enl=enl_rij)
2361      !  grad_berry(:,:) = grad_berry(:,:) - dtefield%efield_dot(idir)*grad_berry_ev(:,:)/two_pi
2362      !  end if
2363 
2364   end do    ! idir
2365 
2366   !deallocations
2367   if(gs_hamk%usepaw /= 0) then
2368      call pawcprj_free(cprj_kb)
2369      call pawcprj_free(cprj_band_srt)
2370      if (nkpt /= dtefield%fnkpt) then
2371         call pawcprj_free(cprj_fkn)
2372         call pawcprj_free(cprj_ikn)
2373      end if
2374   end if
2375   ABI_DEALLOCATE(grad_berry_ev)
2376   ABI_DEALLOCATE(qijbkk)
2377   ABI_DEALLOCATE(enl_rij)
2378   ABI_DEALLOCATE(smat_k_paw)
2379   ABI_DATATYPE_DEALLOCATE(cprj_kb)
2380   ABI_DATATYPE_DEALLOCATE(cprj_band_srt)
2381   ABI_DATATYPE_DEALLOCATE(cprj_fkn)
2382   ABI_DATATYPE_DEALLOCATE(cprj_ikn)
2383 
2384   !DBG_EXIT("COLL")
2385 
2386 end subroutine make_grad_berry

m_cgwf/cgwf [ Functions ]

[ Top ] [ m_cgwf ] [ Functions ]

NAME

 cgwf

FUNCTION

 Update all wavefunction |C>, non self-consistently.
 also compute the corresponding H|C> and Vnl|C> (and S|C> if paw).
 Uses a conjugate-gradient algorithm.
 In case of paw, resolves a generalized eigenproblem using an
  overlap matrix (not used for norm conserving psps).

INPUTS

  berryopt == 4/14: electric field is on;
              all other values, no field is present
  chkexit= if non-zero, check whether the user wishes to exit
  cpus = CPU time limit
  filnam_ds1=name of input file (used for exit checking)
  gs_hamk <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k
  icg=shift to be applied on the location of data in the array cg
  igsc=shift to be applied on the location of data in the array gsc
  ikpt=number of the k-point
  inonsc=index of non self-consistent loop
  isppol=spin polarization currently treated
  mband =maximum number of bands
  mcg=second dimension of the cg array
  mcgq=second dimension of the cgq array
  mgsc=second dimension of the gsc array
  mkgq = second dimension of pwnsfacq
  mpi_enreg=information about MPI parallelization
  mpw=maximum dimensioned size of npw
  nband=number of bands.
  nbdblock=number of bands in a block
  nkpt=number of k points
  nline=number of line minimizations per band.
  npw=number of planewaves in basis sphere at given k.
  npwarr(nkpt)=number of planewaves in basis at this k point
  nspinor=number of spinorial components of the wavefunctions (on current proc)
  nsppol=number of spin polarizations
  ortalg=governs the choice of the algorithm for orthogonalisation.
  prtvol=control print volume and debugging output
  pwind(pwind_alloc,2,3) = array used to compute
           the overlap matrix smat between k-points (see initberry.f)
  pwind_alloc = first dimension of pwind
  pwnsfac(2,pwind_alloc) = phase factors for non-symmorphic translations
                           (see initberry.f)
  pwnsfacq(2,mkgq) = phase factors for the nearest neighbours of the
                     current k-point (electric field, MPI //)
  tolrde=tolerance on the ratio of differences of energies (for the line minimisation)
  tolwfr=tolerance on largest wf residual
  use_subovl=1 if the overlap matrix is not identity in WFs subspace
  wfoptalg=govern the choice of algorithm for wf optimisation
   (0, 1, 10 and 11 : in the present routine, usual CG algorithm ;
   (2 and 3 : use shifted square Hamiltonian)
  zshift(nband)=in case wfoptalg is 2 or 3, shift of the Hamiltonian

OUTPUT

  dphase_k(3) = change in Zak phase for the current k-point in case berryopt = 4/14,6/16,7/17 (electric (displacement) field)
  resid(nband)=wf residual for new states=|(H-e)|C>|^2 (hartree^2)
  subham(nband*(nband+1))=Hamiltonian expressed in sthe WFs subspace
  subovl(nband*(nband+1)*use_subovl)=overlap matrix expressed in sthe WFs subspace
  subvnl(nband*(nband+1)*(1-gs_hamk%usepaw))=non-local Hamiltonian expressed in sthe WFs subspace

SIDE EFFECTS

  cg(2,mcg)
    at input =wavefunction <G|C band,k> coefficients for ALL bands
    at output same as input except that
      the current band, with number 'band' has been updated
  dtefield <type(efield_type)> = variables related to Berry phase
      calculations (see initberry.f)
  quit= if 1, proceeds to smooth ending of the job.
  if(gs_hamk%usepaw==1)
   gsc(2,mgsc)=<G|S|C band,k> coefficients for ALL bands
               where S is the overlap matrix (used only for paw)

NOTES

  1) cg should not be filtered and normalized : it should already be OK at input !
  2) Not sure that that the generalized eigenproblem (when gs_hamk%usepaw=1)
     is compatible with wfoptalg=2 or 3 (use of shifted square  Hamiltonian) - to be verified

PARENTS

      vtowfk

CHILDREN

      cg_precon,cg_zaxpy,cg_zcopy,cg_zscal,dotprod_g,etheta,fock_set_ieigen
      getcprj,getghc,linemin,make_grad_berry,mksubham,pawcprj_alloc
      pawcprj_copy,pawcprj_free,pawcprj_get,pawcprj_mpi_allgather,pawcprj_put
      pawcprj_symkn,projbd,smatrix,smatrix_k_paw,sqnorm_g,timab,wrtout
      xmpi_allgather

SOURCE

 153 subroutine cgwf(berryopt,cg,cgq,chkexit,cpus,dphase_k,dtefield,&
 154 &                filnam_ds1,gsc,gs_hamk,icg,igsc,ikpt,inonsc,&
 155 &                isppol,mband,mcg,mcgq,mgsc,mkgq,mpi_enreg,&
 156 &                mpw,nband,nbdblock,nkpt,nline,npw,npwarr,&
 157 &                nspinor,nsppol,ortalg,prtvol,pwind,&
 158 &                pwind_alloc,pwnsfac,pwnsfacq,quit,resid,subham,subovl,&
 159 &                subvnl,tolrde,tolwfr,use_subovl,wfoptalg,zshift)
 160 
 161 
 162 !This section has been created automatically by the script Abilint (TD).
 163 !Do not modify the following lines by hand.
 164 #undef ABI_FUNC
 165 #define ABI_FUNC 'cgwf'
 166 !End of the abilint section
 167 
 168  implicit none
 169 
 170 !Arguments ------------------------------------
 171  integer,intent(in) :: berryopt,chkexit,icg,igsc,ikpt,inonsc,isppol
 172  integer,intent(in) :: mband,mcg,mcgq,mgsc,mkgq,mpw,nband,nbdblock,nkpt,nline
 173  integer,intent(in) :: npw,nspinor,nsppol,ortalg,prtvol,pwind_alloc
 174  integer,intent(in) :: use_subovl,wfoptalg
 175  integer,intent(in) :: quit
 176  real(dp),intent(in) :: cpus,tolrde,tolwfr
 177  character(len=*),intent(in) :: filnam_ds1
 178  type(MPI_type),intent(in) :: mpi_enreg
 179  type(efield_type),intent(inout) :: dtefield
 180  type(gs_hamiltonian_type),intent(inout) :: gs_hamk
 181 !arrays
 182  integer,intent(in) :: npwarr(nkpt),pwind(pwind_alloc,2,3)
 183  real(dp),intent(in) :: cgq(2,mcgq)
 184  real(dp),intent(in) :: pwnsfac(2,pwind_alloc),pwnsfacq(2,mkgq),zshift(nband)
 185  real(dp),intent(inout) :: cg(2,mcg),gsc(2,mgsc)
 186  real(dp),intent(inout) :: dphase_k(3)
 187  real(dp),intent(out) :: subham(nband*(nband+1)),subovl(nband*(nband+1)*use_subovl)
 188  real(dp),intent(out) :: subvnl(nband*(nband+1)*(1-gs_hamk%usepaw))
 189  real(dp),intent(out) :: resid(nband)
 190 
 191 !Local variables-------------------------------
 192  integer,parameter :: level=113,tim_getghc=1,tim_projbd=1
 193  integer,save :: nskip=0
 194  integer :: choice,counter,cpopt,ddkflag,dimenlc1,dimenlr1,dimenl2,iat,iatom,itypat
 195  integer :: iband,ibandmin,ibandmax,me_g0
 196  integer :: ibdblock,iblock,icg1,icg_shift,icp1,icp2,idir,idum1,ierr,ifor,igs,igsc_shift,ii,ikgf
 197  integer :: ikpt2,ikpt2f,ikptf,iline,iproc,ipw,ispinor,istwf_k,isubh,isubo,itrs
 198  integer :: job,mcg_q,me_distrb,natom,ncpgr,nblock,nproc_distrb,npw_k2
 199  integer :: optekin,paw_opt,signs,shiftbd,sij_opt,spaceComm_distrb
 200  integer :: use_vnl,useoverlap,wfopta10
 201  real(dp) :: chc,costh,deltae,deold,dhc,dhd,diff,dotgg,dotgp,doti,dotr
 202  real(dp) :: dphase_aux2,e0,e0_old,e1,e1_old,eval,gamma
 203  real(dp) :: lam0,lamold,root,sinth,sintn,swap,tan2th,theta,thetam
 204  real(dp) :: xnorm
 205  logical :: gen_eigenpb, finite_field
 206  character(len=500) :: message
 207  integer :: hel(2,3)
 208  integer,allocatable :: dimlmn(:),dimlmn_srt(:),ikptf_recv(:),pwind_k(:),sflag_k(:)
 209  real(dp) :: bcut(2,3),dphase_aux1(3),dtm_k(2),phase_end(3)
 210  real(dp) :: phase_init(3),tsec(2)
 211  real(dp),allocatable :: cg1_k(:,:),cgq_k(:,:),conjgr(:,:),cwavef(:,:)
 212  real(dp),allocatable :: detovc(:,:,:),detovd(:,:,:),direc(:,:),direc_tmp(:,:)
 213  real(dp),allocatable :: gh_direc(:,:),gh_direcws(:,:),ghc(:,:),ghc_all(:,:),ghcws(:,:)
 214  real(dp),allocatable :: grad_berry(:,:),grad_total(:,:),gs_direc(:,:)
 215  real(dp) :: gsc_dummy(0,0)
 216  real(dp),allocatable :: gvnlc(:,:),gvnl_direc(:,:),gvnl_dummy(:,:)
 217  real(dp),allocatable :: pcon(:),pwnsfac_k(:,:),scprod(:,:),scwavef(:,:)
 218  real(dp),allocatable :: smat_inv(:,:,:),smat_k(:,:,:),smat_k_paw(:,:,:),swork(:,:),vresid(:,:),work(:,:)
 219  real(dp),pointer :: kinpw(:)
 220  type(pawcprj_type) :: cprj_dum(1,1)
 221  type(pawcprj_type),allocatable :: cprj_k(:,:),cprj_kb(:,:)
 222  type(pawcprj_type),allocatable :: cprj_direc(:,:),cprj_band_srt(:,:),cprj_gat(:,:)
 223  type(pawcprj_type),allocatable :: cprj_fkn(:,:),cprj_ikn(:,:)
 224 
 225 ! *********************************************************************
 226 
 227  DBG_ENTER("COLL")
 228 
 229 !Starting the routine
 230  call timab(22,1,tsec)
 231 
 232 !Touching chkexit, cpus,filnam_ds to avoid warning for abirules. This is dirty...
 233  if(chkexit<0)then
 234    MSG_BUG('chkexit should be positive!')
 235  end if
 236 
 237  if(cpus<0 .and. filnam_ds1=='a')then
 238    MSG_BUG('cpus should be positive!')
 239  end if
 240 
 241 !======================================================================
 242 !========= LOCAL VARIABLES DEFINITIONS AND ALLOCATIONS ================
 243 !======================================================================
 244 
 245 !MPI data
 246  spaceComm_distrb=mpi_enreg%comm_cell
 247  nproc_distrb=xmpi_comm_size(spaceComm_distrb)
 248  me_distrb=mpi_enreg%me_kpt
 249  me_g0 = mpi_enreg%me_g0
 250 
 251 !if PAW, one has to solve a generalized eigenproblem (H|Psi>=Lambda.S|Psi>)
 252 !else,   one has to solve a classical eigenproblem   (H|Psi>=Lambda.|Psi>)
 253  gen_eigenpb=(gs_hamk%usepaw==1)
 254  useoverlap=0;if (gen_eigenpb) useoverlap=1
 255 
 256 !if PAW, no need to compute Vnl contributions
 257  use_vnl=0; if (gs_hamk%usepaw==0) use_vnl=1
 258  if (gen_eigenpb.and.(use_vnl==1)) then
 259    MSG_BUG("Error in cgwf: contact Abinit group")
 260  end if
 261 
 262 !Initializations and allocations
 263  isubh=1;isubo=1
 264  nblock=(nband-1)/nbdblock+1
 265  istwf_k=gs_hamk%istwf_k
 266  wfopta10=mod(wfoptalg,10)
 267  optekin=0;if (wfoptalg>=10) optekin=1
 268  natom=gs_hamk%natom
 269  cpopt=-1
 270  kinpw => gs_hamk%kinpw_k
 271 
 272  ABI_ALLOCATE(pcon,(npw))
 273  ABI_ALLOCATE(ghc,(2,npw*nspinor))
 274  ABI_ALLOCATE(gvnlc,(2,npw*nspinor))
 275  ABI_ALLOCATE(conjgr,(2,npw*nspinor))
 276  ABI_ALLOCATE(cwavef,(2,npw*nspinor))
 277  ABI_ALLOCATE(direc,(2,npw*nspinor))
 278  ABI_ALLOCATE(scprod,(2,nband))
 279 
 280  ABI_ALLOCATE(gh_direc,(2,npw*nspinor))
 281  ABI_ALLOCATE(gvnl_direc,(2,npw*nspinor))
 282  ABI_ALLOCATE(vresid,(2,npw*nspinor))
 283 
 284  if (gen_eigenpb)  then
 285    ABI_ALLOCATE(gs_direc,(2,npw*nspinor))
 286  else
 287    ABI_ALLOCATE(gs_direc,(0,0))
 288  end if
 289 
 290  if (gen_eigenpb) then
 291    ABI_ALLOCATE(scwavef,(2,npw*nspinor))
 292    ABI_ALLOCATE(direc_tmp,(2,npw*nspinor))
 293  end if
 294 
 295  if (gen_eigenpb.and.(inonsc==1))  then
 296    ABI_STAT_ALLOCATE(ghc_all,(2,nband*npw*nspinor), ierr)
 297    ABI_CHECK(ierr==0, "out-of-memory in ghg_all")
 298  end if
 299 
 300  if (wfopta10==2.or.wfopta10==3)  then
 301    ABI_ALLOCATE(work,(2,npw*nspinor))
 302  end if
 303 
 304  if (gen_eigenpb.and.(wfopta10==2.or.wfopta10==3))  then
 305    ABI_ALLOCATE(swork,(2,npw*nspinor))
 306  else
 307    ABI_ALLOCATE(swork,(0,0))
 308  end if
 309 
 310  if (wfopta10==2 .or. wfopta10==3) then
 311    ABI_ALLOCATE(ghcws,(2,npw*nspinor))
 312    ABI_ALLOCATE(gh_direcws,(2,npw*nspinor))
 313    ABI_ALLOCATE(gvnl_dummy,(2,npw*nspinor))
 314  else
 315    ABI_ALLOCATE(gvnl_dummy,(0,0))
 316  end if
 317 
 318 !if "generalized eigenproblem", not sure of wfoptalg=2,3 algorithms
 319  if ((gen_eigenpb).and.(wfopta10==2.or.wfopta10==3)) then
 320    write(message, '(a,a,a,a,a)' ) &
 321 &   'Conjugate gradient algorithm not tested with',ch10,&
 322 &   'wfoptalg=2 or 3 and usepaw==1 !',ch10,&
 323 &   'Program will continue at your own risk...'
 324    MSG_WARNING(message)
 325  end if
 326 
 327 !Electric field: definition of local variables:
 328 !detovc(1:2,ifor,idir) determinant of the overlap matrix
 329 !S_{nm}(k,k+dk)=<u_{n,k}|u_{m,k+dk}>, with the states at
 330 !k as bras (neighbor is specified by ifor and idir)
 331 !detovd                same as detovc but with <u_{n,k}| replaced by
 332 !<D| (search direction) in the band-th line
 333 !grad_berry(1:2,ipw)   Berry phase term contribution to the gradient vector
 334 !hel(ifor,idir)        helicity of the ellipse associated w/ (ifor,idir)
 335 !bcut(ifor,idir)       branch cut of the ellipse associated w/ (ifor,idir)
 336 !theta_min             optimal angle theta in line_minimization when electric
 337 !field is on
 338 !grad_total(1:2,ipw)   total gradient (zero field term plus Berry phase term)
 339 
 340  finite_field = ( (berryopt == 4) .or. (berryopt == 6) .or. (berryopt == 7) .or.   &
 341 & (berryopt == 14) .or. (berryopt == 16) .or. (berryopt == 17) )
 342  ncpgr = 0 ! do not think the cprj's here need gradients (no force computation in cgwf)
 343 
 344  if (finite_field) then
 345 
 346 !  ji : These could be a couple of new input variables (but it is OK to define them here)
 347    ikptf = dtefield%i2fbz(ikpt)
 348    ikgf = dtefield%fkgindex(ikptf)  ! this is the shift for pwind
 349    mcg_q = mpw*mband*nspinor
 350    ABI_ALLOCATE(detovc,(2,2,3))
 351    ABI_ALLOCATE(detovd,(2,2,3))
 352    ABI_ALLOCATE(grad_berry,(2,npw*nspinor))
 353    ABI_ALLOCATE(cg1_k,(2,mpw))
 354    ABI_ALLOCATE(cgq_k,(2,mcg_q))
 355    ABI_ALLOCATE(grad_total,(2,npw*nspinor))
 356    ABI_ALLOCATE(sflag_k,(mband))
 357    ABI_ALLOCATE(pwind_k,(mpw))
 358    ABI_ALLOCATE(pwnsfac_k,(4,mpw))
 359    ABI_ALLOCATE(smat_k,(2,mband,mband))
 360    ABI_ALLOCATE(smat_inv,(2,mband,mband))
 361 !  now the special features if using PAW
 362    if (gs_hamk%usepaw /= 0) then
 363      ABI_ALLOCATE(smat_k_paw,(2,gs_hamk%usepaw*mband,gs_hamk%usepaw*mband))
 364      smat_k_paw(:,:,:) = zero
 365 !    the following are arguments to nonlop used to apply the on-site dipole to direc vector
 366      choice = 1 ! only apply projectors
 367      paw_opt = 1 ! only apply Dij
 368      signs = 2 ! apply nonlop to vector in k space
 369 !    following two items are the nonlocal potential strength dij due to the on-site dipoles
 370      dimenlc1 = 2*gs_hamk%lmnmax*(gs_hamk%lmnmax+1)/2
 371      dimenlr1 = gs_hamk%lmnmax*(gs_hamk%lmnmax+1)/2
 372      dimenl2 = natom
 373 !    cprj structures for finite_field case
 374      ABI_ALLOCATE(dimlmn,(natom))
 375      do iatom = 1, natom
 376        itypat = gs_hamk%typat(iatom)
 377        dimlmn(iatom)=dtefield%lmn_size(itypat)
 378      end do
 379      ABI_ALLOCATE(dimlmn_srt,(natom))
 380      iatom = 0
 381      do itypat = 1, gs_hamk%ntypat
 382        do iat = 1, gs_hamk%nattyp(itypat)
 383          iatom = iatom + 1
 384          dimlmn_srt(iatom) = dtefield%lmn_size(itypat)
 385        end do
 386      end do
 387      ABI_ALLOCATE(ikptf_recv,(nproc_distrb))
 388      ABI_DATATYPE_ALLOCATE(cprj_k,(natom,mband*nspinor))
 389      ABI_DATATYPE_ALLOCATE(cprj_kb,(natom,mband*nspinor))
 390      ABI_DATATYPE_ALLOCATE(cprj_direc,(natom,mband*nspinor))
 391      ABI_DATATYPE_ALLOCATE(cprj_band_srt,(natom,nspinor))
 392      ABI_DATATYPE_ALLOCATE(cprj_gat,(natom,mband*nspinor*nproc_distrb))
 393      call pawcprj_alloc(cprj_k,ncpgr,dimlmn)
 394      call pawcprj_alloc(cprj_kb,ncpgr,dimlmn)
 395      call pawcprj_alloc(cprj_direc,ncpgr,dimlmn)
 396      call pawcprj_alloc(cprj_band_srt,ncpgr,dimlmn_srt)
 397      call pawcprj_alloc(cprj_gat,ncpgr,dimlmn)
 398      if (nkpt /= dtefield%fnkpt) then
 399        ABI_DATATYPE_ALLOCATE(cprj_fkn,(natom,mband*nspinor))
 400        ABI_DATATYPE_ALLOCATE(cprj_ikn,(natom,mband*nspinor))
 401        call pawcprj_alloc(cprj_fkn,ncpgr,dimlmn)
 402        call pawcprj_alloc(cprj_ikn,ncpgr,dimlmn)
 403      end if
 404    else
 405      ABI_ALLOCATE(smat_k_paw,(0,0,0))
 406      ABI_ALLOCATE(dimlmn,(0))
 407      ABI_ALLOCATE(dimlmn_srt,(0))
 408      ABI_ALLOCATE(ikptf_recv,(0))
 409      ABI_DATATYPE_ALLOCATE(cprj_k,(0,0))
 410      ABI_DATATYPE_ALLOCATE(cprj_kb,(0,0))
 411      ABI_DATATYPE_ALLOCATE(cprj_direc,(0,0))
 412      ABI_DATATYPE_ALLOCATE(cprj_band_srt,(0,0))
 413      ABI_DATATYPE_ALLOCATE(cprj_gat,(0,0))
 414    end if
 415  end if ! finite_field
 416 
 417  ! ======================================================================
 418  ! If generalized eigenproblem: has to know <g|S|c> for all
 419  ! bands (for orthogonalization purpose); take benefit of this
 420  ! calculation to compute <g|H|c> at the same time, and cprj_k if finite_field
 421  ! ======================================================================
 422  if (gen_eigenpb.and.inonsc==1) then
 423    do iblock=1,nblock
 424      ibandmin=1+(iblock-1)*nbdblock
 425      ibandmax=min(iblock*nbdblock,nband)
 426      do iband=ibandmin,ibandmax
 427        ibdblock=iband-(iblock-1)*nbdblock
 428        icg_shift=npw*nspinor*(iband-1)+icg
 429        igsc_shift=npw*nspinor*(iband-1)+igsc
 430 
 431        call cg_zcopy(npw*nspinor,cg(1,1+icg_shift),cwavef)
 432 
 433 !      Compute <g|H|c>
 434 !      By setting ieigen to iband, Fock contrib. of this iband to the energy will be calculated
 435        call fock_set_ieigen(gs_hamk%fockcommon,iband)
 436        sij_opt=1
 437        if (finite_field .and. gs_hamk%usepaw == 1) then
 438          call getghc(0,cwavef,cprj_band_srt,ghc,scwavef,gs_hamk,gvnlc,&
 439 &         eval,mpi_enreg,1,prtvol,sij_opt,tim_getghc,0)
 440          call pawcprj_put(gs_hamk%atindx,cprj_band_srt,cprj_k,natom,iband,0,ikpt,&
 441 &         1,isppol,mband,1,natom,1,mband,dimlmn,nspinor,nsppol,0,&
 442 &         mpicomm=spaceComm_distrb,proc_distrb=mpi_enreg%proc_distrb)
 443        else
 444          call getghc(cpopt,cwavef,cprj_dum,ghc,scwavef,gs_hamk,gvnlc,&
 445 &         eval,mpi_enreg,1,prtvol,sij_opt,tim_getghc,0)
 446        end if
 447 
 448        call cg_zcopy(npw*nspinor,ghc,ghc_all(1,1+icg_shift-icg))
 449        call cg_zcopy(npw*nspinor,scwavef,gsc(1,1+igsc_shift))
 450 
 451      end do
 452    end do
 453  end if
 454 
 455  ! Loop over blocks of bands. In the standard band-sequential algorithm, nblock=nband.
 456  do iblock=1,nblock
 457    counter=100*iblock*nbdblock+inonsc
 458 
 459    ! Loop over bands in a block
 460    ! This loop can be MPI-parallelized, over processors attached to the same k point
 461    ibandmin=1+(iblock-1)*nbdblock
 462    ibandmax=min(iblock*nbdblock,nband)
 463 
 464    ! Big iband loop
 465    do iband=ibandmin,ibandmax
 466      ibdblock=iband-(iblock-1)*nbdblock
 467      counter=100*iband+inonsc
 468      icg_shift=npw*nspinor*(iband-1)+icg
 469      igsc_shift=npw*nspinor*(iband-1)+igsc
 470 
 471      ! ======================================================================
 472      ! ========== INITIALISATION OF MINIMIZATION ITERATIONS =================
 473      ! ======================================================================
 474 
 475      if (prtvol>=10) then ! Tell us what is going on:
 476        write(message, '(a,i6,2x,a,i3,a)' )' --- cgwf is called for band',iband,'for',nline,' lines'
 477        call wrtout(std_out,message,'PERS')
 478      end if
 479 
 480      dotgp=one
 481      if (finite_field) then
 482        detovc(:,:,:) = zero ; detovd(:,:,:) = zero
 483        phase_init(:) = zero
 484        dphase_aux1(:) = zero
 485        phase_end(:) = zero
 486        bcut(:,:) = zero
 487        hel(:,:) = 0
 488      end if
 489 
 490      ! Extraction of the vector that is iteratively updated
 491      call cg_zcopy(npw*nspinor,cg(1,1+icg_shift),cwavef)
 492 
 493      ! If generalized eigenproblem: extraction of the overlap information
 494      if (gen_eigenpb) then
 495        call cg_zcopy(npw*nspinor,gsc(1,1+igsc_shift),scwavef)
 496      end if
 497 
 498      ! Normalize incoming wf (and S.wf, if generalized eigenproblem):
 499      ! WARNING : It might be interesting to skip the following operation.
 500      ! The associated routines should be reexamined to see whether cwavef is not already normalized.
 501      if (gen_eigenpb) then
 502        call dotprod_g(dotr,doti,istwf_k,npw*nspinor,2,cwavef,scwavef,me_g0,mpi_enreg%comm_spinorfft)
 503        dotr=sqrt(dotr**2+doti**2); xnorm=one/sqrt(dotr)
 504        call cg_zscal(npw*nspinor,(/xnorm,zero/),cwavef)
 505        call cg_zscal(npw*nspinor,(/xnorm,zero/),scwavef)
 506      else
 507        call sqnorm_g(dotr,istwf_k,npw*nspinor,cwavef,me_g0,mpi_enreg%comm_fft)
 508        xnorm=one/sqrt(abs(dotr))
 509        call cg_zscal(npw*nspinor,(/xnorm,zero/),cwavef)
 510      end if
 511 
 512      if (prtvol==-level) then
 513        write(message,'(a,f14.6)')' cgwf : xnorm = ',xnorm
 514        call wrtout(std_out,message,'PERS')
 515      end if
 516 
 517      ! Compute (or extract) <g|H|c>
 518      if (gen_eigenpb.and.(inonsc==1)) then
 519 
 520 !$OMP PARALLEL DO PRIVATE(ipw)
 521        do ipw=1,npw*nspinor
 522          ghc(1,ipw)=xnorm*ghc_all(1,ipw+icg_shift-icg)
 523          ghc(2,ipw)=xnorm*ghc_all(2,ipw+icg_shift-icg)
 524        end do
 525 
 526      else
 527 !      By setting ieigen to iband, Fock contrib. of this iband to the energy will be calculated
 528        call fock_set_ieigen(gs_hamk%fockcommon,iband)
 529        sij_opt=0
 530        call getghc(cpopt,cwavef,cprj_dum,ghc,gsc_dummy,gs_hamk,gvnlc,&
 531 &       eval,mpi_enreg,1,prtvol,sij_opt,tim_getghc,0)
 532      end if
 533 
 534 
 535      ! Minimisation of the residual: compute <G|(H-zshift)^2|C iband,k>
 536      if(wfopta10==2 .or. wfopta10==3) then
 537        ghcws(:,:)=ghc(:,:)
 538        if (gen_eigenpb) then
 539          sij_opt=1
 540          work(:,:)=ghc(:,:)-zshift(iband)*scwavef(:,:)
 541        else
 542          sij_opt=0
 543          work(:,:)=ghc(:,:)-zshift(iband)*cwavef(:,:)
 544        end if
 545        call getghc(cpopt,work,cprj_dum,ghc,swork,gs_hamk,gvnl_dummy,&
 546 &       eval,mpi_enreg,1,prtvol,sij_opt,tim_getghc,0)
 547        if (gen_eigenpb) then
 548          ghc(:,:)=ghc(:,:)-zshift(iband)*swork(:,:)
 549        else
 550          ghc(:,:)=ghc(:,:)-zshift(iband)*work(:,:)
 551        end if
 552      end if
 553 
 554      ! ======================================================================
 555      ! ====== BEGIN LOOP FOR A GIVEN BAND: MINIMIZATION ITERATIONS ==========
 556      ! ======================================================================
 557      if(nline/=0)then
 558        do iline=1,nline
 559 
 560          ! === COMPUTE THE RESIDUAL ===
 561 
 562          ! Compute lambda = <C|H|C> or <C|(H-zshift)**2|C>
 563 
 564          call dotprod_g(chc,doti,istwf_k,npw*nspinor,1,cwavef,ghc,me_g0,mpi_enreg%comm_spinorfft)
 565          lam0=chc
 566 
 567          ! Check that lam0 is decreasing on succeeding lines:
 568          if (.not.finite_field) then
 569            if (iline==1) then
 570              lamold=lam0
 571            else
 572              if (lam0 > lamold+tol12) then
 573                write(message, '(a,i8,a,1p,e14.6,a1,3x,a,1p,e14.6,a1)')&
 574 &               'New trial energy at line ',iline,' = ',lam0,ch10,&
 575 &               'is higher than former =',lamold,ch10
 576                MSG_WARNING(message)
 577              end if
 578              lamold=lam0
 579            end if
 580          end if
 581 
 582          ! Compute residual vector:
 583          ! Note that vresid is precomputed to garantee cancellation of errors
 584          ! and allow residuals to reach values as small as 1.0d-24 or better.
 585 
 586          if (wfopta10<=1) then
 587            eval=chc
 588            if (gen_eigenpb) then
 589 
 590 !$OMP PARALLEL DO
 591              do ipw=1,npw*nspinor
 592                vresid(1,ipw)=ghc(1,ipw)-chc*scwavef(1,ipw)
 593                vresid(2,ipw)=ghc(2,ipw)-chc*scwavef(2,ipw)
 594              end do
 595            else
 596 !$OMP PARALLEL DO
 597              do ipw=1,npw*nspinor
 598                vresid(1,ipw)=ghc(1,ipw)-chc*cwavef(1,ipw)
 599                vresid(2,ipw)=ghc(2,ipw)-chc*cwavef(2,ipw)
 600              end do
 601            end if
 602          else
 603            call dotprod_g(eval,doti,istwf_k,npw*nspinor,1,cwavef,ghcws,me_g0,mpi_enreg%comm_spinorfft)
 604            if (gen_eigenpb) then
 605 !$OMP PARALLEL DO
 606              do ipw=1,npw*nspinor
 607                vresid(1,ipw)=ghcws(1,ipw)-eval*scwavef(1,ipw)
 608                vresid(2,ipw)=ghcws(2,ipw)-eval*scwavef(2,ipw)
 609              end do
 610            else
 611 !$OMP PARALLEL DO
 612              do ipw=1,npw*nspinor
 613                vresid(1,ipw)=ghcws(1,ipw)-eval*cwavef(1,ipw)
 614                vresid(2,ipw)=ghcws(2,ipw)-eval*cwavef(2,ipw)
 615              end do
 616            end if
 617          end if
 618 
 619          ! Compute residual (squared) norm
 620          call sqnorm_g(resid(iband),istwf_k,npw*nspinor,vresid,me_g0,mpi_enreg%comm_fft)
 621 
 622          if (prtvol==-level) then
 623            write(message,'(a,i0,2f14.6)')' cgwf : iline,eval,resid = ',iline,eval,resid(iband)
 624            call wrtout(std_out,message,'PERS')
 625          end if
 626 
 627          ! ======================================================================
 628          ! ============== CHECK FOR CONVERGENCE CRITERIA ========================
 629          ! ======================================================================
 630 
 631          ! If residual sufficiently small stop line minimizations
 632          if (resid(iband)<tolwfr) then
 633            if (prtvol>=10) then
 634              write(message, '(a,i4,a,i2,a,es12.4)' ) &
 635 &             ' cgwf: band ',iband,' converged after ',iline,' line minimizations : resid =',resid(iband)
 636              call wrtout(std_out,message,'PERS')
 637            end if
 638            nskip=nskip+(nline-iline+1)  ! Number of two-way 3D ffts skipped
 639            exit                         ! Exit from the loop on iline
 640          end if
 641 
 642          ! If user require exiting the job, stop line minimisations
 643          if (quit==1) then
 644            write(message, '(a,i0)' )' cgwf : user require exiting => skip update of band ',iband
 645            call wrtout(std_out,message,'PERS')
 646 
 647            nskip=nskip+(nline-iline+1)  ! Number of two-way 3D ffts skipped
 648            exit                         ! Exit from the loop on iline
 649          end if
 650 
 651          ! ======================================================================
 652          ! =========== COMPUTE THE STEEPEST DESCENT DIRECTION ===================
 653          ! ======================================================================
 654 
 655          ! Compute the steepest descent direction
 656          if (gen_eigenpb) then
 657            call cg_zcopy(npw*nspinor,vresid,direc)  ! Store <G|H-lambda.S|C> in direc
 658          else
 659            call cg_zcopy(npw*nspinor,ghc,direc)     ! Store <G|H|C> in direc
 660          end if
 661 
 662          ! Electric field: compute the gradient of the Berry phase part of the energy functional.
 663          ! See PRL 89, 117602 (2002) [[cite:Souza2002]], grad_berry(:,:) is the second term of Eq. (4)
 664          if (finite_field) then
 665 
 666            call make_grad_berry(cg,cgq,cprj_k,detovc,dimlmn,dimlmn_srt,direc,dtefield,grad_berry,&
 667 &           gs_hamk,iband,icg,ikpt,isppol,mband,mcg,mcgq,mkgq,mpi_enreg,mpw,natom,nkpt,npw,npwarr,&
 668 &           nspinor,nsppol,pwind,pwind_alloc,pwnsfac,pwnsfacq)
 669 
 670 !          Add grad_berry to direc and store original gradient
 671            direc(:,:) = direc(:,:) + grad_berry(:,:)
 672            grad_total(:,:) = direc(:,:)
 673 !          DEBUG: check that grad_berry is orthogonal to the occupied manifold at k
 674 !          do jband = 1, dtefield%mband_occ
 675 !          dotr = zero  ;  doti = zero
 676 !          do ipw = 1, npw*nspinor
 677 !          if(.not.gen_eigenpb) then
 678 !          dotr = dotr + cg(1,icg + (jband-1)*npw*nspinor + ipw)*grad_berry(1,ipw) + &
 679 !          &                 cg(2,icg + (jband-1)*npw*nspinor + ipw)*grad_berry(2,ipw)
 680 !          doti = doti + cg(1,icg + (jband-1)*npw*nspinor + ipw)*grad_berry(2,ipw) - &
 681 !          &                 cg(2,icg + (jband-1)*npw*nspinor + ipw)*grad_berry(1,ipw)
 682 !          end if
 683 !          end do
 684 !          if ((abs(dotr) > tol12).or.(abs(doti) > tol12)) then
 685 !          write(std_out,'(a)')'cgwf-berry : ERROR (orthogonality)'
 686 !          write(std_out,'(3(2x,i3),2(5x,e16.9))')ikpt,iband,jband,dotr,doti
 687 !          stop
 688 !          end if
 689 !          end do
 690 !          ENDDEBUG
 691          end if   ! finite_field
 692 
 693          ! =========== PROJECT THE STEEPEST DESCENT DIRECTION ===================
 694          ! ========= OVER THE SUBSPACE ORTHOGONAL TO OTHER BANDS ================
 695 
 696          ! The following projection over the subspace orthogonal to occupied bands
 697          ! is optional. It is a bit more accurate, but doubles the number of N^3 ops.
 698          ! It is done only if ortalg>=0.
 699 
 700          ! Project the steepest descent direction:
 701          ! direc(2,npw)=<G|H|Cnk> - \sum_{(i<=n)} <G|H|Cik> , normalized.
 702 
 703          if(ortalg>=0)then
 704            if (gen_eigenpb) then
 705              call projbd(cg,direc,iband,icg,igsc,istwf_k,mcg,mgsc,nband,npw,nspinor,&
 706 &             gsc,scprod,0,tim_projbd,useoverlap,me_g0,mpi_enreg%comm_fft)
 707            else
 708              call projbd(cg,direc,-1   ,icg,igsc,istwf_k,mcg,mgsc,nband,npw,nspinor,&
 709 &             gsc,scprod,0,tim_projbd,useoverlap,me_g0,mpi_enreg%comm_fft)
 710            end if
 711          else
 712 !          For negative ortalg must still project current band out of conjugate vector (unneeded if gen_eigenpb)
 713            if (.not.gen_eigenpb) then
 714              call dotprod_g(dotr,doti,istwf_k,npw*nspinor,3,cwavef,direc,me_g0,mpi_enreg%comm_spinorfft)
 715              if(istwf_k==1)then
 716                call cg_zaxpy(npw*nspinor,-(/dotr,doti/),cwavef,direc)
 717              else
 718                call cg_zaxpy(npw*nspinor,(/-dotr,zero/),cwavef,direc)
 719              end if
 720            end if
 721          end if
 722 
 723          ! For a generalized eigenpb, store the steepest descent direction
 724          if (gen_eigenpb) direc_tmp=direc
 725 
 726          ! ======================================================================
 727          ! ======== PRECONDITION THE STEEPEST DESCENT DIRECTION =================
 728          ! ======================================================================
 729 
 730          ! If wfoptalg>=10, the precondition matrix is kept constant during iteration ; otherwise it is recomputed
 731          if (wfoptalg<10.or.iline==1) then
 732            call cg_precon(cwavef,zero,istwf_k,kinpw,npw,nspinor,me_g0,optekin,pcon,direc,mpi_enreg%comm_fft)
 733 
 734            ! Minimisation of the residual: must precondition twice
 735            ! (might make only one call, with modified precon routine - might also make a shift !!!)
 736            if(wfopta10==2 .or. wfopta10==3)then
 737              call cg_precon(cwavef,zero,istwf_k,kinpw,npw,nspinor,me_g0,optekin,pcon,direc,mpi_enreg%comm_fft)
 738              if(iline==1)then
 739 !$OMP PARALLEL DO
 740                do ipw=1,npw
 741                  pcon(ipw)=pcon(ipw)**2
 742                  pcon(ipw)=pcon(ipw)**2
 743                end do
 744              end if
 745            end if
 746          else
 747            do ispinor=1,nspinor
 748              igs=(ispinor-1)*npw
 749 !$OMP PARALLEL DO
 750              do ipw=1+igs,npw+igs
 751                direc(1,ipw)=direc(1,ipw)*pcon(ipw-igs)
 752                direc(2,ipw)=direc(2,ipw)*pcon(ipw-igs)
 753              end do
 754            end do
 755          end if
 756 
 757          ! ======= PROJECT THE PRECOND. STEEPEST DESCENT DIRECTION ==============
 758          ! ========= OVER THE SUBSPACE ORTHOGONAL TO OTHER BANDS ================
 759          ! Projecting again out all bands (not normalized).
 760          call projbd(cg,direc,-1,icg,igsc,istwf_k,mcg,mgsc,nband,npw,nspinor,&
 761 &         gsc,scprod,0,tim_projbd,useoverlap,me_g0,mpi_enreg%comm_fft)
 762 
 763          ! ======================================================================
 764          ! ================= COMPUTE THE CONJUGATE-GRADIENT =====================
 765          ! ======================================================================
 766 
 767          if (finite_field) then
 768            call dotprod_g(dotgg,doti,istwf_k,npw*nspinor,1,direc,grad_total,me_g0,mpi_enreg%comm_spinorfft)
 769            !DEBUG (electric field)
 770            !check that the dotproduct is real
 771            !if (abs(doti) > tol8) then
 772            !write(std_out,*) ' cgwf-berry: ERROR'
 773            !write(std_out,*) ' doti = ',doti
 774            !stop
 775            !end if
 776            !ENDDEBUG
 777          else
 778            if (gen_eigenpb) then
 779              call dotprod_g(dotgg,doti,istwf_k,npw*nspinor,1,direc,direc_tmp,me_g0,mpi_enreg%comm_spinorfft)
 780            else
 781              call dotprod_g(dotgg,doti,istwf_k,npw*nspinor,1,direc,ghc,me_g0,mpi_enreg%comm_spinorfft)
 782            end if
 783          end if
 784 
 785          ! MJV: added 5 Feb 2012 - causes divide by 0 on next iteration of iline
 786          if (abs(dotgg) < TINY(0.0_dp)*1.e50_dp) dotgg = TINY(0.0_dp)*1.e50_dp
 787 
 788          ! At first iteration, gamma is set to zero
 789          if (iline==1) then
 790            gamma=zero
 791            dotgp=dotgg
 792            call cg_zcopy(npw*nspinor,direc,conjgr)
 793 
 794          else
 795            gamma=dotgg/dotgp
 796            dotgp=dotgg
 797 
 798            if(prtvol==-level)then
 799              write(message,'(a,2es16.6)')' cgwf : dotgg,gamma = ',dotgg,gamma
 800              call wrtout(std_out,message,'PERS')
 801            end if
 802 
 803            ! Note: another way to compute gamma: Polak, Ribiere no real improvement ; to be more carrefully tested
 804            ! call dotprod_g(dotgg,doti,istwf_k,mpi_enreg,npw*nspinor,1,direc,direc_tmp)
 805            ! !direcp must be set to zero at the beginning
 806            ! direcp=direc-direcp
 807            ! call dotprod_g(dotgmg,doti,istwf_k,mpi_enreg,npw*nspinor,1,direcp,direc_tmp)
 808            ! direcp=direc;gamma=dotgmg/dotgp;dotgp=dotgmg
 809 
 810 !$OMP PARALLEL DO
 811            do ipw=1,npw*nspinor
 812              conjgr(1,ipw)=direc(1,ipw)+gamma*conjgr(1,ipw)
 813              conjgr(2,ipw)=direc(2,ipw)+gamma*conjgr(2,ipw)
 814            end do
 815            !call cg_zaxpby(npw*nspinor,cg_one,direc,(/gamma,zero/),conjgr)
 816          end if
 817 
 818          ! ======================================================================
 819          ! ============ PROJECTION OF THE CONJUGATED GRADIENT ===================
 820          ! ======================================================================
 821 
 822          if (gen_eigenpb) then
 823            call dotprod_g(dotr,doti,istwf_k,npw*nspinor,3,scwavef,conjgr,me_g0,mpi_enreg%comm_spinorfft)
 824          else
 825            call dotprod_g(dotr,doti,istwf_k,npw*nspinor,3,cwavef,conjgr,me_g0,mpi_enreg%comm_spinorfft)
 826          end if
 827 
 828          ! Project the conjugated gradient onto the current band
 829          ! MG: TODO: this is an hot spot that could be rewritten with BLAS! provided
 830          ! that direc --> conjgr
 831          if(istwf_k==1)then
 832 
 833 !$OMP PARALLEL DO
 834            do ipw=1,npw*nspinor
 835              direc(1,ipw)=conjgr(1,ipw)-(dotr*cwavef(1,ipw)-doti*cwavef(2,ipw))
 836              direc(2,ipw)=conjgr(2,ipw)-(dotr*cwavef(2,ipw)+doti*cwavef(1,ipw))
 837            end do
 838          else
 839 !$OMP PARALLEL DO
 840            do ipw=1,npw*nspinor
 841              direc(1,ipw)=conjgr(1,ipw)-dotr*cwavef(1,ipw)
 842              direc(2,ipw)=conjgr(2,ipw)-dotr*cwavef(2,ipw)
 843            end do
 844          end if
 845 
 846          ! In case of generalized eigenproblem, normalization of direction vector
 847          ! cannot be done here (because S|D> is not known here).
 848          if (.not.gen_eigenpb) then
 849            call sqnorm_g(dotr,istwf_k,npw*nspinor,direc,me_g0,mpi_enreg%comm_fft)
 850            xnorm=one/sqrt(abs(dotr))
 851            call cg_zscal(npw*nspinor,(/xnorm,zero/),direc)
 852            xnorm=one
 853          end if
 854 
 855          ! ======================================================================
 856          ! ===== COMPUTE CONTRIBUTIONS TO 1ST AND 2ND DERIVATIVES OF ENERGY =====
 857          ! ======================================================================
 858 
 859          ! Compute gh_direc = <G|H|D> and eventually gs_direc = <G|S|D>
 860          sij_opt=0;if (gen_eigenpb) sij_opt=1
 861 
 862          call getghc(cpopt,direc,cprj_dum,gh_direc,gs_direc,gs_hamk,gvnl_direc,&
 863 &         eval,mpi_enreg,1,prtvol,sij_opt,tim_getghc,0)
 864 
 865          if(wfopta10==2 .or. wfopta10==3)then
 866            ! Minimisation of the residual, so compute <G|(H-zshift)^2|D>
 867            gh_direcws(:,:)=gh_direc(:,:)
 868            if (gen_eigenpb) then
 869              sij_opt=1
 870              work(:,:)=gh_direc(:,:)-zshift(iband)*gs_direc(:,:)
 871            else
 872              sij_opt=0
 873              work(:,:)=gh_direc(:,:)-zshift(iband)*direc(:,:)
 874            end if
 875 
 876            call getghc(cpopt,work,cprj_dum,gh_direc,swork,gs_hamk,gvnl_dummy,&
 877 &           eval,mpi_enreg,1,prtvol,0,tim_getghc,0)
 878 
 879            if (gen_eigenpb) then
 880              gh_direc(:,:)=gh_direc(:,:)-zshift(iband)*swork(:,:)
 881            else
 882              gh_direc(:,:)=gh_direc(:,:)-zshift(iband)*work(:,:)
 883            end if
 884          end if
 885 
 886          ! In case of generalized eigenproblem, compute now the norm of the conjugated gradient
 887          if (gen_eigenpb) then
 888            call dotprod_g(dotr,doti,istwf_k,npw*nspinor,1,direc,gs_direc,me_g0,mpi_enreg%comm_spinorfft)
 889            xnorm=one/sqrt(abs(dotr))
 890          end if
 891 
 892          ! Compute dhc = Re{<D|H|C>}
 893          call dotprod_g(dhc,doti,istwf_k,npw*nspinor,1,direc,ghc,me_g0,mpi_enreg%comm_spinorfft)
 894          dhc=dhc*xnorm
 895 
 896          ! Compute <D|H|D> or <D|(H-zshift)^2|D>
 897          call dotprod_g(dhd,doti,istwf_k,npw*nspinor,1,direc,gh_direc,me_g0,mpi_enreg%comm_spinorfft)
 898          dhd=dhd*xnorm**2
 899 
 900          if(prtvol==-level)then
 901            write(message,'(a,3f14.6)') 'cgwf : chc,dhc,dhd=',chc,dhc,dhd
 902            call wrtout(std_out,message,'PERS')
 903          end if
 904 
 905          ! ======================================================================
 906          ! ======= COMPUTE MIXING FACTORS - CHECK FOR CONVERGENCE ===============
 907          ! ======================================================================
 908 
 909          if (.not.finite_field) then
 910            ! Compute tan(2 theta),sin(theta) and cos(theta)
 911            tan2th=2.0_dp*dhc/(chc-dhd)
 912 
 913            if (abs(tan2th)<1.d-05) then
 914              costh=1.0_dp-0.125_dp*tan2th**2
 915              sinth=0.5_dp*tan2th*(1.0_dp-0.375_dp*tan2th**2)
 916 
 917              ! Check that result is above machine precision
 918              if (abs(sinth)<epsilon(0._dp)) then
 919                write(message, '(a,es16.4)' ) ' cgwf: converged with tan2th=',tan2th
 920                call wrtout(std_out,message,'PERS')
 921                ! Number of one-way 3D ffts skipped
 922                nskip=nskip+2*(nline-iline)
 923                exit ! Exit from the loop on iline
 924              end if
 925 
 926            else
 927              root=sqrt(1.0_dp+tan2th**2)
 928              costh=sqrt(0.5_dp+0.5_dp/root)
 929              sinth=sign(sqrt(0.5_dp-0.5_dp/root),tan2th)
 930            end if
 931 
 932            ! Check for lower of two possible roots (same sign as curvature at theta where slope is zero)
 933            diff=(chc-dhd)
 934            ! Swap c and d if value of diff is positive
 935            if (diff>zero) then
 936              swap=costh
 937              costh=-sinth
 938              sinth=swap
 939              if(prtvol<0 .or. prtvol>=10)then
 940                write(message,*)'   Note: swap roots, iline,diff=',iline,diff
 941                call wrtout(std_out,message,'PERS')
 942              end if
 943            end if
 944 
 945          else
 946            ! In case the eletric field is on, the line minimization has to be done numerically
 947 
 948            ! Compute determinant of the overlap matrix where in the band-th line
 949            ! the wavefunction is replaced by the search direction
 950            job = 10 ; shiftbd = 0
 951            do idir = 1, 3
 952              ! do not do this for efield_dot(idir)=0
 953              if (abs(dtefield%efield_dot(idir)) < tol12) cycle
 954              do ifor = 1, 2
 955                ikpt2f = dtefield%ikpt_dk(ikptf,ifor,idir)
 956                if (dtefield%indkk_f2ibz(ikpt2f,6) == 1) then
 957                  itrs = 10
 958                else
 959                  itrs = 0
 960                end if
 961                ikpt2 = dtefield%indkk_f2ibz(ikpt2f,1)
 962                npw_k2 = npwarr(ikpt2)
 963                pwind_k(1:npw) = pwind(ikgf+1:ikgf+npw,ifor,idir)
 964                pwnsfac_k(1:2,1:npw) = pwnsfac(1:2,ikgf+1:ikgf+npw)
 965                sflag_k(:) = dtefield%sflag(:,ikpt+(isppol-1)*nkpt,ifor,idir)
 966                smat_k(:,:,:) = dtefield%smat(:,:,:,ikpt+(isppol-1)*nkpt,ifor,idir)
 967 
 968                if (mpi_enreg%nproc_cell > 1) then
 969                  icg1 = dtefield%cgqindex(2,ifor+2*(idir-1),ikpt+(isppol-1)*nkpt)
 970                  cgq_k(:,1:dtefield%mband_occ*nspinor*npw_k2) = &
 971 &                 cgq(:,icg1+1:icg1+dtefield%mband_occ*nspinor*npw_k2)
 972                  idum1 = dtefield%cgqindex(3,ifor+2*(idir-1),ikpt+(isppol-1)*nkpt)
 973                  pwnsfac_k(3:4,1:npw_k2) = pwnsfacq(1:2,idum1+1:idum1+npw_k2)
 974                else
 975                  icg1 = dtefield%cgindex(ikpt2,isppol)
 976                  cgq_k(:,1:dtefield%mband_occ*nspinor*npw_k2) = &
 977 &                 cg(:,icg1+1:icg1+dtefield%mband_occ*nspinor*npw_k2)
 978                  idum1=dtefield%fkgindex(ikpt2f)
 979                  pwnsfac_k(3:4,1:npw_k2) = pwnsfac(1:2,idum1+1:idum1+npw_k2)
 980                end if
 981 
 982                icg1 = 0 ; ddkflag = 0
 983                if (gen_eigenpb) then
 984 !$OMP PARALLEL DO
 985                  do ipw=1,npw*nspinor
 986                    direc_tmp(1,ipw)=direc(1,ipw)*xnorm
 987                    direc_tmp(2,ipw)=direc(2,ipw)*xnorm
 988                  end do
 989                  ! need cprj corresponding to direc_tmp in order to make smat_k_paw properly
 990                  call getcprj(1,0,direc_tmp,cprj_band_srt,&
 991 &                 gs_hamk%ffnl_k,0,gs_hamk%indlmn,gs_hamk%istwf_k,gs_hamk%kg_k,&
 992 &                 gs_hamk%kpg_k,gs_hamk%kpt_k,gs_hamk%lmnmax,gs_hamk%mgfft,&
 993 &                 mpi_enreg,gs_hamk%natom,gs_hamk%nattyp,gs_hamk%ngfft,gs_hamk%nloalg,&
 994 &                 gs_hamk%npw_k,gs_hamk%nspinor,gs_hamk%ntypat,gs_hamk%phkxred,gs_hamk%ph1d,&
 995 &                 gs_hamk%ph3d_k,gs_hamk%ucvol,gs_hamk%useylm)
 996 
 997                  call pawcprj_copy(cprj_k,cprj_direc)
 998                  call pawcprj_put(gs_hamk%atindx,cprj_band_srt,cprj_direc,gs_hamk%natom,iband,0,&
 999 &                 ikpt,1,isppol,mband,1,gs_hamk%natom,1,mband,dimlmn,gs_hamk%nspinor,nsppol,0,&
1000 &                 mpicomm=spaceComm_distrb,proc_distrb=mpi_enreg%proc_distrb)
1001 
1002 !                icp1=dtefield%mband_occ*(ikptf-1)
1003                  icp2=mband*nspinor*(ikpt2-1)
1004                  call pawcprj_get(gs_hamk%atindx,cprj_kb,dtefield%cprj,natom,1,icp2,ikpt,0,isppol,&
1005 &                 mband,dtefield%fnkpt,natom,mband,mband,nspinor,nsppol,0,&
1006 &                 mpicomm=spaceComm_distrb,proc_distrb=mpi_enreg%proc_distrb)
1007 
1008                  if (ikpt2 /= ikpt2f) then ! construct cprj_kb by symmetry
1009                    call pawcprj_copy(cprj_kb,cprj_ikn)
1010                    call pawcprj_symkn(cprj_fkn,cprj_ikn,dtefield%atom_indsym,dimlmn,-1,gs_hamk%indlmn,&
1011 &                   dtefield%indkk_f2ibz(ikpt2f,2),dtefield%indkk_f2ibz(ikpt2f,6),&
1012 &                   dtefield%fkptns(:,dtefield%i2fbz(ikpt2)),&
1013 &                   dtefield%lmax,dtefield%lmnmax,mband,natom,dtefield%mband_occ,nspinor,&
1014 &                   dtefield%nsym,gs_hamk%ntypat,gs_hamk%typat,dtefield%zarot)
1015                    call pawcprj_copy(cprj_fkn,cprj_kb)
1016                  end if
1017 
1018                  call smatrix_k_paw(cprj_direc,cprj_kb,dtefield,idir,ifor,mband,natom,smat_k_paw,gs_hamk%typat)
1019 
1020                  call smatrix(direc_tmp,cgq_k,cg1_k,ddkflag,dtm_k,icg1,icg1,&
1021 &                 itrs,job,iband,npw*nspinor,mcg_q,mpw,iband,&
1022 &                 mpw,dtefield%mband_occ,dtefield%nband_occ(isppol),&
1023 &                 npw,npw_k2,nspinor,pwind_k,pwnsfac_k,sflag_k,&
1024 &                 shiftbd,smat_inv,smat_k,smat_k_paw,gs_hamk%usepaw)
1025                else
1026                  call smatrix(direc,cgq_k,cg1_k,ddkflag,dtm_k,icg1,icg1,&
1027 &                 itrs,job,iband,npw*nspinor,mcg_q,mpw,iband,&
1028 &                 mpw,dtefield%mband_occ,dtefield%nband_occ(isppol),&
1029 &                 npw,npw_k2,nspinor,pwind_k,pwnsfac_k,sflag_k,&
1030 &                 shiftbd,smat_inv,smat_k,smat_k_paw,gs_hamk%usepaw)
1031                end if
1032                detovd(:,ifor,idir) = dtm_k(:) ! Store the determinant of the overlap
1033 !              matrix (required to compute theta_min)
1034 !              DEBUG
1035 !              write(std_out,*)'cgwf-berry: detovc and detovd'
1036 !              write(std_out,*)detovc(:,ifor,idir)
1037 !              write(std_out,*)detovd(:,ifor,idir)
1038 !              write(std_out,*)'smat_k'
1039 !              do jband = 1, 4
1040 !              write(std_out,'(4(2x,e14.6))')smat_k(1,jband,:)
1041 !              write(std_out,'(4(2x,e14.6))')smat_k(2,jband,:)
1042 !              write(std_out,*)
1043 !              end do
1044 !              ENDDEBUG
1045              end do  ! ifor
1046            end do    ! idir
1047 
1048            call linemin(bcut,chc,costh,detovc,detovd,dhc,dhd,&
1049 &           dphase_aux1,dtefield%efield_dot,iline,&
1050 &           dtefield%fnkpt,dtefield%nstr,hel,phase_end,&
1051 &           phase_init,dtefield%sdeg,sinth,thetam)
1052 !          DEBUG
1053 !          if (mpi_enreg%me == 1) then
1054 !          write(std_out,*)'after linemin '
1055 !          write(std_out,'(a,3(2x,f16.9))')'phase_init  = ',phase_init(:)
1056 !          write(std_out,'(a,3(2x,f16.9))')'phase_end   = ',phase_end(:)
1057 !          write(std_out,'(a,3(2x,f16.9))')'dphase_aux1 = ',dphase_aux1(:)
1058 !          write(std_out,*) 'thetam',thetam
1059 !          end if
1060 !          ENDDEBUG
1061          end if  ! finite_field
1062 
1063          ! ======================================================================
1064          ! =========== GENERATE NEW |wf>, H|wf>, Vnl|Wf>, S|Wf> ... =============
1065          ! ======================================================================
1066 
1067          sintn=sinth*xnorm
1068 
1069 !$OMP PARALLEL DO
1070          do ipw=1,npw*nspinor
1071            cwavef(1,ipw)=cwavef(1,ipw)*costh+direc(1,ipw)*sintn
1072            cwavef(2,ipw)=cwavef(2,ipw)*costh+direc(2,ipw)*sintn
1073          end do
1074 
1075 !        call cg_zaxpby(npw*nspinor,(/sintn,zero/),direc,(/costh,zero/),cwavef)
1076          call cg_zcopy(npw*nspinor,cwavef,cg(1,1+icg_shift))
1077 
1078          do ipw=1,npw*nspinor
1079            ghc(1,ipw)  =ghc(1,ipw)*costh + gh_direc(1,ipw)*sintn
1080            ghc(2,ipw)  =ghc(2,ipw)*costh + gh_direc(2,ipw)*sintn
1081          end do
1082 
1083 
1084          if (use_vnl==1) then
1085 !$OMP PARALLEL DO
1086            do ipw=1,npw*nspinor
1087              gvnlc(1,ipw)=gvnlc(1,ipw)*costh + gvnl_direc(1,ipw)*sintn
1088              gvnlc(2,ipw)=gvnlc(2,ipw)*costh + gvnl_direc(2,ipw)*sintn
1089            end do
1090 !          call cg_zaxpby(npw*nspinor,(/sintn,zero/),gvnl_direc,(/costh,zero/),gvnlc)
1091          end if
1092 
1093          if (gen_eigenpb) then
1094 !$OMP PARALLEL DO
1095            do ipw=1,npw*nspinor
1096              scwavef(1,ipw)=scwavef(1,ipw)*costh+gs_direc(1,ipw)*sintn
1097              scwavef(2,ipw)=scwavef(2,ipw)*costh+gs_direc(2,ipw)*sintn
1098 !            gsc(1,ipw+igsc_shift)=scwavef(1,ipw)
1099 !            gsc(2,ipw+igsc_shift)=scwavef(2,ipw)
1100            end do
1101 !          call cg_zaxpby(npw*nspinor,(/sintn,zero/),gs_direc,(/costh,zero/),scwavef)
1102            call cg_zcopy(npw*nspinor,scwavef,gsc(1,1+igsc_shift))
1103 
1104            if (finite_field) then  ! must update cprj for the new wavefunction
1105              call getcprj(1,0,cwavef,cprj_band_srt,&
1106 &             gs_hamk%ffnl_k,0,gs_hamk%indlmn,istwf_k,gs_hamk%kg_k,gs_hamk%kpg_k,gs_hamk%kpt_k,&
1107 &             gs_hamk%lmnmax,gs_hamk%mgfft,mpi_enreg,natom,gs_hamk%nattyp,&
1108 &             gs_hamk%ngfft,gs_hamk%nloalg,gs_hamk%npw_k,gs_hamk%nspinor,gs_hamk%ntypat,&
1109 &             gs_hamk%phkxred,gs_hamk%ph1d,gs_hamk%ph3d_k,gs_hamk%ucvol,gs_hamk%useylm)
1110              call pawcprj_put(gs_hamk%atindx,cprj_band_srt,cprj_k,gs_hamk%natom,iband,0,ikpt,&
1111 &             1,isppol,mband,1,gs_hamk%natom,1,mband,dimlmn,gs_hamk%nspinor,nsppol,0,&
1112 &             mpicomm=spaceComm_distrb,proc_distrb=mpi_enreg%proc_distrb)
1113            end if
1114          end if
1115 
1116          if(wfopta10==2 .or. wfopta10==3)then
1117            ! Need to keep track of ghcws, in order to avoid recomputing it
1118 !$OMP PARALLEL DO
1119            do ipw=1,npw*nspinor
1120              ghcws(1,ipw)=ghcws(1,ipw)*costh + gh_direcws(1,ipw)*sintn
1121              ghcws(2,ipw)=ghcws(2,ipw)*costh + gh_direcws(2,ipw)*sintn
1122            end do
1123 !          call cg_zaxpby(npw*nspinor,(/sintn,zero/),gh_direcws,(/costh,zero/),ghcws)
1124          end if
1125 
1126          ! ======================================================================
1127          ! =========== CHECK CONVERGENCE AGAINST TRIAL ENERGY ===================
1128          ! ======================================================================
1129 
1130          ! Compute delta(E)
1131          if (.not.finite_field) then
1132            deltae=chc*(costh**2-1._dp)+dhd*sinth**2+2._dp*costh*sinth*dhc
1133          else
1134            ! Compute deltae
1135            call etheta(bcut,chc,detovc,detovd,dhc,dhd,dtefield%efield_dot,e0,e1,&
1136 &           hel,dtefield%fnkpt,dtefield%nstr,dtefield%sdeg,thetam)
1137            theta = zero
1138 
1139            call etheta(bcut,chc,detovc,detovd,dhc,dhd,&
1140 &           dtefield%efield_dot,e0_old,e1_old,&
1141 &           hel,dtefield%fnkpt,dtefield%nstr,dtefield%sdeg,theta)
1142            deltae = e0 - e0_old
1143 !          DEBUG
1144 !          write(std_out,*) 'e0, e0_old, deltae', e0, e0_old, deltae
1145 !          ENDDEBUG
1146 !          Check that e0 is decreasing on succeeding lines:
1147 !          if (deltae > zero) then
1148            if (deltae > tol12) then ! exploring different checks for finit_field
1149              write(message, '(3a,i8,a,1p,e14.6,a1,3x,a,1p,e14.6,a1)')&
1150 &             '  (electric field)',ch10,&
1151 &             '  New trial energy at line',iline,' = ',e0,ch10,&
1152 &             '  is higher than former:',e0_old,ch10
1153              MSG_WARNING(message)
1154            end if
1155          end if         ! finite_field
1156 
1157 !        Check convergence and eventually exit
1158          if (iline==1) then
1159            deold=deltae
1160          else if (abs(deltae)<tolrde*abs(deold) .and. iline/=nline .and. wfopta10<2)then
1161            if(prtvol>=10)then
1162              write(message, '(a,i4,1x,a,1p,e12.4,a,e12.4,a)' ) &
1163 &             ' cgwf: line',iline,&
1164 &             ' deltae=',deltae,' < tolrde*',deold,' =>skip lines'
1165              call wrtout(std_out,message,'PERS')
1166            end if
1167            nskip=nskip+2*(nline-iline)  ! Number of one-way 3D ffts skipped
1168            exit                         ! Exit from the loop on iline
1169          end if
1170 
1171        end do ! END LOOP FOR A GIVEN BAND Note that there are three "exit" instructions inside
1172 
1173        ! Additional computations in case of electric field
1174        if (finite_field) then
1175          ! Bring present contribution to dphasek(idir) into [-pi,pi]
1176          do idir = 1, 3
1177            dphase_aux2 = mod(phase_end(idir) - phase_init(idir) + 100*two_pi,two_pi)
1178            if (dphase_aux2 > pi) dphase_aux2 = dphase_aux2 - two_pi
1179            ! DEBUG
1180            ! dphase_aux1(idir)=mod(dphase_aux1(idir)+100*two_pi,two_pi)
1181            ! if(dphase_aux1(idir)>pi) dphase_aux1(idir)=dphase_aux1(idir)-two_pi
1182            ! diff = dphase_aux2 - dphase_aux1(idir)
1183            ! if (abs(diff) > tol10) then
1184            ! write(std_out,*)'cgwf-berry: ERROR'
1185            ! write(std_out,'(a,3(2x,i3),f16.9)')'ikpt,iband,idir,diff',ikpt,iband,idir,diff
1186            ! stop
1187            ! end if
1188            ! write(100,*) idir, dphase_aux2
1189            ! ENDDEBUG
1190            dphase_k(idir) = dphase_k(idir) + dphase_aux2
1191            ! DEBUG
1192            ! write(std_out,*) 'idir,phase_init,phase_end,dphase_k'
1193            ! write(std_out,*) idir,phase_init(idir),phase_end(idir),dphase_k(idir)
1194            ! ENDDEBUG
1195          end do
1196        end if   ! finite_field
1197 
1198      else ! nline==0 , needs to provide a residual
1199        resid(iband)=-one
1200      end if ! End nline==0 case
1201 
1202      ! ======================================================================
1203      ! =============== END OF CURRENT BAND: CLEANING ========================
1204      ! ======================================================================
1205 
1206      ! It was checked that getghc is NOT needed here : equivalent results with the copy below.
1207      if(wfopta10==2 .or. wfopta10==3) ghc(:,:)=ghcws(:,:)
1208 
1209      if (finite_field) dtefield%sflag(:,ikpt + (isppol-1)*nkpt,:,:) = 0
1210 
1211      ! At the end of the treatment of a set of bands, write the number of one-way 3D ffts skipped
1212      if (xmpi_paral==0 .and. mpi_enreg%paral_kgb==0 .and. iband==nband .and. prtvol/=0) then
1213        write(message,'(a,i0)')' cgwf : number of one-way 3D ffts skipped in cgwf until now =',nskip
1214        call wrtout(std_out,message,'PERS')
1215      end if
1216 
1217    end do !  End big iband loop. iband in a block
1218 
1219    !  ======================================================================
1220    !  ============= COMPUTE HAMILTONIAN IN WFs SUBSPACE ====================
1221    !  ======================================================================
1222    call mksubham(cg,ghc,gsc,gvnlc,iblock,icg,igsc,istwf_k,&
1223 &   isubh,isubo,mcg,mgsc,nband,nbdblock,npw,&
1224 &   nspinor,subham,subovl,subvnl,use_subovl,use_vnl,me_g0)
1225 
1226  end do ! iblock End loop over block of bands
1227 
1228  if (allocated(dimlmn_srt)) then
1229    ABI_DEALLOCATE(dimlmn_srt)
1230  end if
1231 
1232  if (finite_field .and. gs_hamk%usepaw == 1) then ! store updated cprjs for this kpt
1233    ! switch from ikptf to ikpt
1234    ikptf = ikpt
1235    call xmpi_allgather(ikptf,ikptf_recv,spaceComm_distrb,ierr)
1236    call pawcprj_mpi_allgather(cprj_k,cprj_gat,natom,nspinor*mband,dimlmn,ncpgr,nproc_distrb,&
1237 &   spaceComm_distrb,ierr,rank_ordered=.true.)
1238    do iproc = 1, nproc_distrb
1239      icp2=nspinor*mband*(iproc-1)
1240      call pawcprj_get(gs_hamk%atindx1,cprj_k,cprj_gat,natom,1,icp2,ikpt,0,isppol,&
1241 &     mband,nproc_distrb,natom,mband,mband,nspinor,nsppol,0,&
1242 &     mpicomm=spaceComm_distrb,proc_distrb=mpi_enreg%proc_distrb)
1243 !    ikptf = ikptf_recv(iproc)
1244      icp1 = nspinor*mband*(ikptf_recv(iproc)-1)
1245      call pawcprj_put(gs_hamk%atindx1,cprj_k,dtefield%cprj,natom,1,icp1,ikpt,0,isppol,&
1246 &     mband,nkpt,natom,mband,mband,dimlmn,nspinor,nsppol,0,&
1247 &     mpicomm=spaceComm_distrb,proc_distrb=mpi_enreg%proc_distrb)
1248    end do
1249  end if
1250 
1251  ! Debugging ouputs
1252  if(prtvol==-level)then
1253    isubh=1
1254    if (use_vnl==1) write(message,'(a)') ' cgwf : isubh  subham(isubh:isubh+1)  subvnl(isubh:isubh+1)'
1255    if (use_vnl==0) write(message,'(a)') ' cgwf : isubh  subham(isubh:isubh+1)'
1256    do iband=1,nband
1257      do ii=1,iband
1258        if (use_vnl==1) then
1259          write(message,'(i5,4es16.6)')isubh,subham(isubh:isubh+1),subvnl(isubh:isubh+1)
1260        else
1261          write(message,'(i5,2es16.6)')isubh,subham(isubh:isubh+1)
1262        end if
1263        call wrtout(std_out,message,'PERS')
1264        isubh=isubh+2
1265      end do
1266    end do
1267  end if
1268 
1269  ! ===================
1270  ! FINAL DEALLOCATIONS
1271  ! ===================
1272  ABI_DEALLOCATE(conjgr)
1273  ABI_DEALLOCATE(cwavef)
1274  ABI_DEALLOCATE(direc)
1275  ABI_DEALLOCATE(pcon)
1276  ABI_DEALLOCATE(scprod)
1277  ABI_DEALLOCATE(ghc)
1278  ABI_DEALLOCATE(gvnlc)
1279  ABI_DEALLOCATE(gh_direc)
1280  ABI_DEALLOCATE(gvnl_direc)
1281  ABI_DEALLOCATE(vresid)
1282  ABI_DEALLOCATE(gs_direc)
1283 
1284  if (gen_eigenpb)  then
1285    ABI_DEALLOCATE(scwavef)
1286    ABI_DEALLOCATE(direc_tmp)
1287  end if
1288 
1289  if (gen_eigenpb.and.(inonsc==1))  then
1290    ABI_DEALLOCATE(ghc_all)
1291  end if
1292 
1293  if(wfopta10==2.or.wfopta10==3)  then
1294    ABI_DEALLOCATE(ghcws)
1295    ABI_DEALLOCATE(gh_direcws)
1296  end if
1297  ABI_DEALLOCATE(gvnl_dummy)
1298 
1299  if(wfopta10==2.or.wfopta10==3)  then
1300    ABI_DEALLOCATE(work)
1301  end if
1302 
1303  ABI_DEALLOCATE(swork)
1304 
1305  if (finite_field) then
1306    ABI_DEALLOCATE(cg1_k)
1307    ABI_DEALLOCATE(cgq_k)
1308    ABI_DEALLOCATE(detovc)
1309    ABI_DEALLOCATE(detovd)
1310    ABI_DEALLOCATE(grad_berry)
1311    ABI_DEALLOCATE(sflag_k)
1312    ABI_DEALLOCATE(smat_inv)
1313    ABI_DEALLOCATE(smat_k)
1314    ABI_DEALLOCATE(pwind_k)
1315    ABI_DEALLOCATE(pwnsfac_k)
1316    ABI_DEALLOCATE(grad_total)
1317    if (gs_hamk%usepaw /= 0) then
1318      call pawcprj_free(cprj_k)
1319      call pawcprj_free(cprj_kb)
1320      call pawcprj_free(cprj_direc)
1321      call pawcprj_free(cprj_band_srt)
1322      call pawcprj_free(cprj_gat)
1323      if (nkpt /= dtefield%fnkpt) then
1324        call pawcprj_free(cprj_fkn)
1325        call pawcprj_free(cprj_ikn)
1326        ABI_DATATYPE_DEALLOCATE(cprj_fkn)
1327        ABI_DATATYPE_DEALLOCATE(cprj_ikn)
1328      end if
1329    end if
1330    ABI_DEALLOCATE(smat_k_paw)
1331    ABI_DEALLOCATE(dimlmn)
1332    ABI_DATATYPE_DEALLOCATE(cprj_k)
1333    ABI_DATATYPE_DEALLOCATE(cprj_kb)
1334    ABI_DATATYPE_DEALLOCATE(cprj_direc)
1335    ABI_DATATYPE_DEALLOCATE(cprj_gat)
1336    ABI_DEALLOCATE(ikptf_recv)
1337    ABI_DATATYPE_DEALLOCATE(cprj_band_srt)
1338  end if
1339 
1340 ! Do not delete this line, needed to run with open MP
1341  write(unit=message,fmt=*) resid(1)
1342 
1343  call timab(22,2,tsec)
1344 
1345  DBG_EXIT("COLL")
1346 
1347 end subroutine cgwf

m_cgwf/etheta [ Functions ]

[ Top ] [ m_cgwf ] [ Functions ]

NAME

 etheta

FUNCTION

 Computes the energy per unit cell and its first derivative
 for a given angle theta. More precisely, computes only the part of
 the energy that changes with theta.

INPUTS

 bcut(ifor,idir) = branch cut of the ellipse associated with (ifor,idir)
 chc = <C|H_0|C> where |C> is the wavefunction of the current band
 detovc = determinant of the overlap matrix S
 detovd = determinant of the overlap matrix where for the band
          that is being updated <C| is replaced by <D| (search direction)
 dhc = Re[<D|H_0|C>]
 dhd = <D|H_0|D>
 efield_dot = reciprocal lattice coordinates of the electric field
 hel(ifor,idir) = helicity of the ellipse associated with (ifor,idir)
 nkpt = number of k-points
 nsppol = 1 for unpolarized, 2 for spin-polarized
 nstr(idir) = number of strings along the idir-th direction
 sdeg = spin degeneracy
 theta = value of the angle for which the energy (e0) and its
         derivative (e1) are computed

OUTPUT

 e0 = energy for the given value of theta
 e1 = derivative of the energy with respect to theta

PARENTS

      cgwf,linemin

CHILDREN

      rhophi

SOURCE

1730 subroutine etheta(bcut,chc,detovc,detovd,dhc,dhd,efield_dot,e0,e1,&
1731 &    hel,nkpt,nstr,sdeg,theta)
1732 
1733 
1734 !This section has been created automatically by the script Abilint (TD).
1735 !Do not modify the following lines by hand.
1736 #undef ABI_FUNC
1737 #define ABI_FUNC 'etheta'
1738 !End of the abilint section
1739 
1740  implicit none
1741 
1742 !Arguments ------------------------------------
1743 !scalars
1744  integer,intent(in) :: nkpt
1745  real(dp),intent(in) :: chc,dhc,dhd,sdeg,theta
1746  real(dp),intent(out) :: e0,e1
1747 !arrays
1748  integer,intent(in) :: hel(2,3),nstr(3)
1749  real(dp),intent(in) :: bcut(2,3),detovc(2,2,3),detovd(2,2,3),efield_dot(3)
1750 
1751 !Local variables -------------------------
1752 !scalars
1753  integer :: idir,ifor
1754  real(dp) :: c2theta,ctheta,dphase,gnorm,phase,rho,s2theta,sgn,stheta
1755 !arrays
1756  real(dp) :: dg_theta(2),g_theta(2)
1757 
1758 ! ***********************************************************************
1759 
1760  e0 = zero ; e1 = zero
1761 
1762  ctheta = cos(theta)
1763  stheta = sin(theta)
1764  c2theta = ctheta*ctheta - stheta*stheta   ! cos(2*theta)
1765  s2theta = two*ctheta*stheta               ! sin(2*theta)
1766 
1767  e0 = chc*ctheta*ctheta + dhd*stheta*stheta + dhc*s2theta
1768  e0 = e0*sdeg/nkpt
1769 
1770 !DEBUG
1771 !e0 = zero
1772 !ENDDEBUG
1773 
1774  e1 = (dhd - chc)*s2theta + two*dhc*c2theta
1775  e1 = e1*sdeg/nkpt
1776 
1777  sgn = -1_dp
1778  do idir = 1, 3
1779 
1780    if (abs(efield_dot(idir)) < tol12) cycle
1781 
1782    do ifor = 1, 2
1783 
1784      g_theta(:)  = ctheta*detovc(:,ifor,idir) + &
1785 &     stheta*detovd(:,ifor,idir)
1786      dg_theta(:) = -1_dp*stheta*detovc(:,ifor,idir) + &
1787 &     ctheta*detovd(:,ifor,idir)
1788 
1789 !    Compute E(theta)
1790 
1791      call rhophi(g_theta,phase,rho)
1792      if (theta >= bcut(ifor,idir)) phase = phase + hel(ifor,idir)*two_pi
1793 
1794 !    DEBUG
1795 !    unit = 100 + 10*idir + ifor
1796 !    write(unit,'(4(f16.9))')theta,g_theta(:),phase
1797 !    ENDDEBUG
1798 
1799      e0 = e0 + sgn*sdeg*efield_dot(idir)*phase/(two_pi*nstr(idir))
1800 
1801 
1802 !    Compute dE/dtheta
1803 
1804 !    imaginary part of the derivative of ln(g_theta)
1805      gnorm = g_theta(1)*g_theta(1) + g_theta(2)*g_theta(2)
1806      dphase = (dg_theta(2)*g_theta(1) - dg_theta(1)*g_theta(2))/gnorm
1807 
1808      e1 = e1 + sgn*sdeg*efield_dot(idir)*dphase/(two_pi*nstr(idir))
1809 
1810      sgn = -1_dp*sgn
1811 
1812    end do
1813  end do
1814 
1815 end subroutine etheta

m_cgwf/linemin [ Functions ]

[ Top ] [ m_cgwf ] [ Functions ]

NAME

 linemin

FUNCTION

 Performs the "line minimization" w.r.t. the angle theta on a unit circle
 to update the wavefunction associated with the current k-point and
 band label.
 This routine is used only when the electric field is on (otherwise it could
 in principle also be used, but there is a simpler procedure, as originally
 coded in abinit).

INPUTS

 chc = <C|H_0|C> where |C> is the wavefunction of the current band
 detovc = determinant of the overlap matrix S
 detovd = determinant of the overlap matrix where for the band
          that is being updated <C| is replaced by <D| (search direction)
 dhc = Re[<D|H_0|C>]
 dhd = <D|H_0|D>
 efield_dot = reciprocal lattice coordinates of the electric field
 iline = index of the current line minimization
 nkpt = number of k-points
 nstr(idir) = number of strings along the idir-th direction
 sdeg = spin degeneracy

OUTPUT

 bcut(ifor,idir) = branch cut of the ellipse associated with (ifor,idir)
 costh = cos(thetam)
 hel(ifor,idir) = helicity of the ellipse associated with (ifor,idir)
 phase_end = total change in Zak phase, must be equal to
             dphase_aux1 + n*two_pi
 sinth = sin(thetam)
 thetam = optimal angle theta in line_minimization

SIDE EFFECTS

 Input/Output
 dphase_aux1 = change in Zak phase accumulated during the loop over iline
               (can be used for debugging in cgwf.f)
 phase_init = initial Zak phase (before doing the first line minimization)

NOTES

 We are making the "frozen Hamiltonian approximation", i.e., the
 Hamiltonian does not change with theta (we are neglecting the dependence
 of the Hartree and exchange-correlation terms on theta; the original
 abinit routine does the same)

PARENTS

      cgwf

CHILDREN

      etheta,rhophi,wrtout

SOURCE

1405 subroutine linemin(bcut,chc,costh,detovc,detovd,dhc,dhd,dphase_aux1,&
1406 &  efield_dot,iline,nkpt,nstr,hel,phase_end,phase_init,sdeg,sinth,thetam)
1407 
1408 
1409 !This section has been created automatically by the script Abilint (TD).
1410 !Do not modify the following lines by hand.
1411 #undef ABI_FUNC
1412 #define ABI_FUNC 'linemin'
1413 !End of the abilint section
1414 
1415  implicit none
1416 
1417 !Arguments ------------------------------------
1418 !scalars
1419  integer,intent(in) :: iline,nkpt
1420  real(dp),intent(in) :: chc,dhc,dhd,sdeg
1421  real(dp),intent(out) :: costh,sinth,thetam
1422 !arrays
1423  integer,intent(in) :: nstr(3)
1424  integer,intent(out) :: hel(2,3)
1425  real(dp),intent(in) :: detovc(2,2,3),detovd(2,2,3),efield_dot(3)
1426  real(dp),intent(inout) :: dphase_aux1(3),phase_init(3)
1427  real(dp),intent(out) :: bcut(2,3),phase_end(3)
1428 
1429 !Local variables -------------------------
1430 !scalars
1431  integer :: idir,ifor,igrid,iter,maxiter,ngrid
1432  real(dp) :: aa,angle,bb,big_axis,cc,cphi_0,delta_theta,e0,e1
1433  real(dp) :: excentr,iab,phase0,phase_min,phi_0,rdum,sgn,small_axis,sphi_0
1434  real(dp) :: theta,theta_0,val
1435  logical :: flag_neg
1436  character(len=500) :: message
1437 !arrays
1438  real(dp) :: g_theta(2),theta_min(2),theta_try(2)
1439  real(dp) :: esave(251),e1save(251)   !!REC
1440 
1441 ! ***********************************************************************
1442 
1443 !Compute the helicity and the branch cut of the ellipse in the complex
1444 !plane associated with the overlap between a k-point and one of its neighbours
1445 
1446  do idir = 1, 3
1447 
1448    if (abs(efield_dot(idir)) < tol12) cycle
1449 
1450    do ifor = 1, 2
1451 
1452      aa = half*(detovc(1,ifor,idir)*detovc(1,ifor,idir) + &
1453 &     detovc(2,ifor,idir)*detovc(2,ifor,idir) + &
1454 &     detovd(1,ifor,idir)*detovd(1,ifor,idir) + &
1455 &     detovd(2,ifor,idir)*detovd(2,ifor,idir))
1456 
1457      bb = half*(detovc(1,ifor,idir)*detovc(1,ifor,idir) + &
1458 &     detovc(2,ifor,idir)*detovc(2,ifor,idir) - &
1459 &     detovd(1,ifor,idir)*detovd(1,ifor,idir) - &
1460 &     detovd(2,ifor,idir)*detovd(2,ifor,idir))
1461 
1462      cc = detovc(1,ifor,idir)*detovd(1,ifor,idir) + &
1463 &     detovc(2,ifor,idir)*detovd(2,ifor,idir)
1464 
1465      iab = detovc(1,ifor,idir)*detovd(2,ifor,idir) - &
1466 &     detovc(2,ifor,idir)*detovd(1,ifor,idir)
1467 
1468      if (iab >= zero) then
1469        hel(ifor,idir) = 1
1470      else
1471        hel(ifor,idir) = -1
1472      end if
1473 
1474      if (abs(bb) > tol8) then
1475        theta_0 = half*atan(cc/bb)
1476      else
1477        theta_0 = quarter*pi
1478      end if
1479 
1480      if (bb < zero) theta_0 = theta_0 + pi*half
1481 
1482      g_theta(:) = cos(theta_0)*detovc(:,ifor,idir) + &
1483 &     sin(theta_0)*detovd(:,ifor,idir)
1484 !    DEBUG
1485 !    write(std_out,*)'before rhophi, g_theta =',g_theta
1486 !    ENDDEBUG
1487      call rhophi(g_theta,phi_0,rdum)
1488 !    DEBUG
1489 !    write(std_out,*)'after rhophi, phi_0 = ',phi_0
1490 !    ENDDEBUG
1491 
1492      cphi_0 = cos(phi_0)
1493      sphi_0 = sin(phi_0)
1494 
1495      rdum = aa - sqrt(bb*bb + cc*cc)
1496      if (rdum < zero) rdum = zero
1497      small_axis = sqrt(rdum)
1498      big_axis = sqrt(aa + sqrt(bb*bb + cc*cc))
1499      excentr = hel(ifor,idir)*small_axis/big_axis
1500 
1501 !    Find angle for which phi = pi
1502      if (abs(excentr) > tol8) then
1503        angle = atan(tan(pi-phi_0)/excentr)
1504      else
1505        if (tan(pi-phi_0)*hel(ifor,idir) > zero) then
1506          angle = half*pi
1507        else
1508          angle = -0.5_dp*pi
1509        end if
1510      end if
1511      bcut(ifor,idir) = angle + theta_0
1512 
1513 
1514 !    Compute the branch-cut angle
1515      if (hel(ifor,idir) == 1) then
1516        if ((sphi_0 > 0).and.(cphi_0 > 0)) bcut(ifor,idir) = bcut(ifor,idir) + pi
1517        if ((sphi_0 < 0).and.(cphi_0 > 0)) bcut(ifor,idir) = bcut(ifor,idir) - pi
1518      else
1519        if ((sphi_0 > 0).and.(cphi_0 > 0)) bcut(ifor,idir) = bcut(ifor,idir) - pi
1520        if ((sphi_0 < 0).and.(cphi_0 > 0)) bcut(ifor,idir) = bcut(ifor,idir) + pi
1521      end if
1522 
1523      if (bcut(ifor,idir) > pi) bcut(ifor,idir) = bcut(ifor,idir) - two_pi
1524      if (bcut(ifor,idir) < -1_dp*pi) bcut(ifor,idir) = bcut(ifor,idir) + two_pi
1525 
1526 !    DEBUG
1527 !    write(std_out,'(a,2x,i3,2x,i3,5x,f16.9,5x,i2)')'linemin: ifor,idir,bcut,hel',&
1528 !    &   ifor,idir,bcut(ifor,idir),hel(ifor,idir)
1529 !    write(std_out,'(a,5x,f16.9,5x,f16.9)')'linemin: big_axis,small_axis ',&
1530 !    &     big_axis,small_axis
1531 !    ENDDEBUG
1532 
1533    end do   ! ifor
1534  end do   ! idir
1535 
1536 !---------------------------------------------------------------------------
1537 
1538 !Perform the "line minimization" w.r.t. the angle theta on a unit circle
1539 !to update the wavefunction associated with the current k-point and band label.
1540 
1541  ngrid = 250   ! initial number of subdivisions in [-pi/2,pi/2]
1542 !for finding extrema
1543  maxiter = 100
1544  delta_theta = pi/ngrid
1545 
1546 !DEBUG
1547 !write(std_out,*)'linemin: theta, e0, e1, e1fdiff'
1548 !ENDDEBUG
1549 
1550 
1551 !Get the interval where the absolute minimum of E(theta) is located
1552 
1553  val = huge(one)             ! large number
1554  flag_neg=.false.
1555  theta_min(:) = ten
1556  do igrid = 1, ngrid+1
1557 
1558    theta = (igrid - 1)*delta_theta - pi*half
1559    call etheta(bcut,chc,detovc,detovd,dhc,dhd,efield_dot,e0,e1,&
1560 &   hel,nkpt,nstr,sdeg,theta)
1561 
1562    esave(igrid)=e0      !!REC
1563    e1save(igrid)=e1     !!REC
1564 
1565 !  It is important to detect when the slope changes from negative to positive
1566 !  Moreover, a slope being extremely close to zero must be ignored
1567 
1568 !  DEBUG
1569 !  write(std_out,*)' igrid,e0,e1,val,theta_min(:)=',igrid,theta,e0,e1,val,theta_min(:)
1570 !  ENDDEBUG
1571 
1572 !  Store e1 and theta if negative ...
1573    if(e1 < -tol10)then
1574      theta_try(1)=theta
1575      flag_neg=.true.
1576    end if
1577 !  A change of sign is just happening
1578    if(e1 > tol10 .and. flag_neg)then
1579      theta_try(2)=theta
1580      flag_neg=.false.
1581 !    Still, must be better than the previous minimum in order to succeed
1582      if (e0 < val-tol10) then
1583        val=e0
1584        theta_min(:)=theta_try(:)
1585      end if
1586    end if
1587  end do
1588 
1589 !In case the minimum was not found
1590 
1591  if (abs(theta_min(1) - ten) < tol10) then
1592 !  REC start
1593    write(message,'(a,a)')ch10,&
1594 &   ' linemin: ERROR- cannot find theta_min.'
1595    call wrtout(std_out,message,'COLL')
1596    write(message,'(a,a)')ch10,&
1597 &   ' igrid      theta          esave(igrid)    e1save(igrid) '
1598    call wrtout(std_out,message,'COLL')
1599    do igrid = 1, ngrid+1
1600      theta = (igrid - 1)*delta_theta - pi*half
1601      write(std_out,'(i6,3f16.9)')igrid,theta,esave(igrid),e1save(igrid)
1602      !write(101,'(i6,3f16.9)')igrid,theta,esave(igrid),e1save(igrid)
1603    end do
1604    write(message,'(6a)')ch10,&
1605 &   ' linemin : ERROR - ',ch10,&
1606 &   '  Cannot find theta_min. No minimum exists : the field is too strong ! ',ch10,&
1607 &   '  Try decreasing difference between D and 4 Pi P by changing structure or D (only for fixed D calculation)'
1608    call wrtout(std_out,message,'COLL')
1609    message = ' linemin cannot find theta_min'
1610    MSG_ERROR(message)
1611  end if
1612 
1613 !Compute the mimum of E(theta)
1614 
1615 
1616  iter = 0
1617  do while ((delta_theta > tol8).and.(iter < maxiter))
1618    delta_theta = half*(theta_min(2) - theta_min(1))
1619    theta = theta_min(1) + delta_theta
1620    call etheta(bcut,chc,detovc,detovd,dhc,dhd,efield_dot,e0,e1,&
1621 &   hel,nkpt,nstr,sdeg,theta)
1622    if (e1 > zero) then
1623      theta_min(2) = theta
1624    else
1625      theta_min(1) = theta
1626    end if
1627    iter = iter + 1
1628 
1629 !  DEBUG
1630 !  write(std_out,'(a,2x,i3,2(2x,f16.9))')'iter,e0,e1 = ',iter,e0,e1
1631 !  ENDDEBUG
1632 
1633  end do
1634 
1635  costh = cos(theta)
1636  sinth = sin(theta)
1637 
1638  thetam = theta
1639 
1640 !DEBUG
1641 !write(std_out,*)'linemin : thetam = ',thetam
1642 !ENDDEBUG
1643 
1644 
1645 !---------------------------------------------------------------------------
1646 
1647 !Compute and store the change in electronic polarization
1648 
1649  sgn = one
1650  do idir = 1, 3
1651 
1652    if (abs(efield_dot(idir)) < tol12) cycle
1653 
1654    phase_end(idir) = zero
1655    do ifor = 1, 2
1656 
1657      g_theta(:) = detovc(:,ifor,idir)
1658 !    DEBUG
1659 !    write(std_out,*)'before rhophi (2nd call), g_theta =',g_theta
1660 !    ENDDEBUG
1661      call rhophi(g_theta,phase0,rdum)
1662 !    DEBUG
1663 !    write(std_out,*)'after rhophi, phase0 = ',phase0
1664 !    ENDDEBUG
1665 
1666      if(iline == 1) phase_init(idir) = phase_init(idir) + sgn*phase0
1667 
1668      g_theta(:) = costh*detovc(:,ifor,idir) + sinth*detovd(:,ifor,idir)
1669      call rhophi(g_theta,phase_min,rdum)
1670 
1671      phase_end(idir) = phase_end(idir) + sgn*phase_min
1672 
1673 !    Correct for branch cuts (remove jumps)
1674      if (bcut(ifor,idir) <= zero) phase0 = phase0 + hel(ifor,idir)*two_pi
1675      if(thetam >= bcut(ifor,idir)) phase_min = phase_min + hel(ifor,idir)*two_pi
1676 
1677      dphase_aux1(idir) = dphase_aux1(idir) + sgn*(phase_min - phase0)
1678 
1679      sgn = -1_dp*sgn
1680 
1681    end do   ! idir
1682  end do    ! ifor
1683 
1684 !DEBUG
1685 !write(std_out,'(a,3(2x,f16.9))')'dphase_aux1 = ',(dphase_aux1(idir),idir = 1, 3)
1686 !write(std_out,*)' linemin: debug, exit.'
1687 !ENDDEBUG
1688 
1689 end subroutine linemin

m_cgwf/mksubham [ Functions ]

[ Top ] [ m_cgwf ] [ Functions ]

NAME

 mksubham

FUNCTION

 Build the Hamiltonian matrix in the eigenfunctions subspace,
 for one given band (or for one given block of bands)

INPUTS

  cg(2,mcg)=wavefunctions
  gsc(2,mgsc)=<g|S|c> matrix elements (S=overlap)
  iblock=index of block of bands
  icg=shift to be applied on the location of data in the array cg
  igsc=shift to be applied on the location of data in the array cg
  istwf_k=input parameter that describes the storage of wfs
  mcg=second dimension of the cg array
  mgsc=second dimension of the gsc array
  nband_k=number of bands at this k point for that spin polarization
  nbdblock=number of bands in a block
  npw_k=number of plane waves at this k point
  nspinor=number of spinorial components of the wavefunctions
  use_subovl=1 if the overlap matrix is not identity in WFs subspace
  use_vnl= 1 if <C band,k|H|C band_prime,k> has to be computed
  me_g0=1 if this processors has G=0, 0 otherwise

OUTPUT

SIDE EFFECTS

  ghc(2,npw_k*nspinor)=<G|H|C band,k> for the current state
                       This is an input in non-blocked algorithm
                               an output in blocked algorithm
  gvnlc(2,npw_k*nspinor)=<G|Vnl|C band,k> for the current state
                       This is an input in non-blocked algorithm
                               an output in blocked algorithm
  isubh=index of current state in array subham
  isubo=index of current state in array subovl
  subham(nband_k*(nband_k+1))=Hamiltonian expressed in the WFs subspace
  subovl(nband_k*(nband_k+1)*use_subovl)=overlap matrix expressed in the WFs subspace
  subvnl(nband_k*(nband_k+1)*use_vnl)=non-local Hamiltonian expressed in the WFs subspace

PARENTS

      cgwf

CHILDREN

SOURCE

1865 subroutine mksubham(cg,ghc,gsc,gvnlc,iblock,icg,igsc,istwf_k,&
1866 &                    isubh,isubo,mcg,mgsc,nband_k,nbdblock,npw_k,&
1867 &                    nspinor,subham,subovl,subvnl,use_subovl,use_vnl,me_g0)
1868 
1869 
1870 !This section has been created automatically by the script Abilint (TD).
1871 !Do not modify the following lines by hand.
1872 #undef ABI_FUNC
1873 #define ABI_FUNC 'mksubham'
1874 !End of the abilint section
1875 
1876  implicit none
1877 
1878 !Arguments ------------------------------------
1879 !scalars
1880  integer,intent(in) :: iblock,icg,igsc,istwf_k,mcg,mgsc,nband_k
1881  integer,intent(in) :: nbdblock,npw_k,nspinor,use_subovl,use_vnl,me_g0
1882  integer,intent(inout) :: isubh,isubo
1883 !arrays
1884  real(dp),intent(in) :: cg(2,mcg)
1885  real(dp),intent(in) :: gsc(2,mgsc)
1886  real(dp),intent(inout) :: ghc(2,npw_k*nspinor),gvnlc(2,npw_k*nspinor)
1887  real(dp),intent(inout) :: subham(nband_k*(nband_k+1))
1888  real(dp),intent(inout) :: subovl(nband_k*(nband_k+1)*use_subovl)
1889  real(dp),intent(inout) :: subvnl(nband_k*(nband_k+1)*use_vnl)
1890 
1891 !Local variables-------------------------------
1892 !scalars
1893  integer :: iband,ibdblock,ii,ipw,ipw1,isp,iwavef,jwavef
1894  real(dp) :: cgimipw,cgreipw,chcim,chcre,cscim,cscre,cvcim,cvcre
1895 !real(dp) :: chc(2),cvc(2),csc(2)
1896 
1897 ! *********************************************************************
1898 
1899 !Loop over bands in a block This loop can be parallelized
1900  do iband=1+(iblock-1)*nbdblock,min(iblock*nbdblock,nband_k)
1901    ibdblock=iband-(iblock-1)*nbdblock
1902 
1903 !  Compute elements of subspace Hamiltonian <C(i)|H|C(n)> and <C(i)|Vnl|C(n)>
1904    if(istwf_k==1)then
1905 
1906      do ii=1,iband
1907        iwavef=(ii-1)*npw_k*nspinor+icg
1908        chcre=zero ; chcim=zero
1909        if (use_vnl==0) then
1910          do ipw=1,npw_k*nspinor
1911            cgreipw=cg(1,ipw+iwavef)
1912            cgimipw=cg(2,ipw+iwavef)
1913            chcre=chcre+cgreipw*ghc(1,ipw)+cgimipw*ghc(2,ipw)
1914            chcim=chcim+cgreipw*ghc(2,ipw)-cgimipw*ghc(1,ipw)
1915          end do
1916 !        chc = cg_zdotc(npw_k*nspinor,cg(1,1+iwavef),ghc)
1917        else
1918 #if 1
1919          do ipw=1,npw_k*nspinor
1920            cgreipw=cg(1,ipw+iwavef)
1921            cgimipw=cg(2,ipw+iwavef)
1922            chcre=chcre+cgreipw*ghc(1,ipw)+cgimipw*ghc(2,ipw)
1923            chcim=chcim+cgreipw*ghc(2,ipw)-cgimipw*ghc(1,ipw)
1924          end do
1925          cvcre=zero ; cvcim=zero
1926          do ipw=1,npw_k*nspinor
1927            cgreipw=cg(1,ipw+iwavef)
1928            cgimipw=cg(2,ipw+iwavef)
1929            cvcre=cvcre+cgreipw*gvnlc(1,ipw)+cgimipw*gvnlc(2,ipw)
1930            cvcim=cvcim+cgreipw*gvnlc(2,ipw)-cgimipw*gvnlc(1,ipw)
1931          end do
1932          subvnl(isubh  )=cvcre
1933          subvnl(isubh+1)=cvcim
1934 #else
1935 !        New version with BLAS1, will require some update of the refs.
1936          cvc = cg_zdotc(npw_k*nspinor,cg(1,1+iwavef),gvnlc)
1937          subvnl(isubh  )=cvc(1)
1938          subvnl(isubh+1)=cvc(2)
1939          chc = cg_zdotc(npw_k*nspinor,cg(1,1+iwavef),ghc)
1940          chcre = chc(1)
1941          chcim = chc(2)
1942 #endif
1943 !        Store real and imag parts in Hermitian storage mode:
1944        end if
1945        subham(isubh  )=chcre
1946        subham(isubh+1)=chcim
1947 !      subham(isubh  )=chc(1)
1948 !      subham(isubh+1)=chc(2)
1949        isubh=isubh+2
1950      end do
1951 
1952    else if(istwf_k>=2)then
1953      do ii=1,iband
1954        iwavef=(ii-1)*npw_k+icg
1955 !      Use the time-reversal symmetry, but should not double-count G=0
1956        if(istwf_k==2 .and. me_g0==1) then
1957          chcre = half*cg(1,1+iwavef)*ghc(1,1)
1958          if (use_vnl==1) cvcre=half*cg(1,1+iwavef)*gvnlc(1,1)
1959          ipw1=2
1960        else
1961          chcre=zero; ipw1=1
1962          if (use_vnl==1) cvcre=zero
1963        end if
1964        if (use_vnl==0) then
1965          do isp=1,nspinor
1966            do ipw=ipw1+(isp-1)*npw_k,npw_k*isp
1967              cgreipw=cg(1,ipw+iwavef)
1968              cgimipw=cg(2,ipw+iwavef)
1969              chcre=chcre+cgreipw*ghc(1,ipw)+cgimipw*ghc(2,ipw)
1970            end do
1971          end do
1972          chcre=two*chcre
1973        else
1974          do isp=1,nspinor
1975            do ipw=ipw1+(isp-1)*npw_k,npw_k*isp
1976              cgreipw=cg(1,ipw+iwavef)
1977              cgimipw=cg(2,ipw+iwavef)
1978              chcre=chcre+cgreipw*ghc(1,ipw)+cgimipw*ghc(2,ipw)
1979              cvcre=cvcre+cgreipw*gvnlc(1,ipw)+cgimipw*gvnlc(2,ipw)
1980            end do
1981          end do
1982          chcre=two*chcre
1983          cvcre=two*cvcre
1984 !        Store real and imag parts in Hermitian storage mode:
1985          subvnl(isubh  )=cvcre
1986          subvnl(isubh+1)=zero
1987        end if
1988        subham(isubh  )=chcre
1989        subham(isubh+1)=zero
1990        isubh=isubh+2
1991      end do
1992    end if
1993 
1994 !  Compute elements of subspace <C(i)|S|C(n)> (S=overlap matrix)
1995 !  <C(i)|S|C(n)> should be closed to Identity.
1996    if (use_subovl==1) then
1997      jwavef=(iband-1)*npw_k*nspinor+igsc
1998      if(istwf_k==1)then
1999        do ii=1,iband
2000          iwavef=(ii-1)*npw_k*nspinor+icg
2001          cscre=zero ; cscim=zero
2002          do ipw=1,npw_k*nspinor
2003            cgreipw=cg(1,ipw+iwavef)
2004            cgimipw=cg(2,ipw+iwavef)
2005            cscre=cscre+cgreipw*gsc(1,ipw+jwavef)+cgimipw*gsc(2,ipw+jwavef)
2006            cscim=cscim+cgreipw*gsc(2,ipw+jwavef)-cgimipw*gsc(1,ipw+jwavef)
2007          end do
2008 !        csc = cg_zdotc(npw_k*nspinor,cg(1,1+iwavef),gsc)
2009 !        subovl(isubo  )=csc(1)
2010 !        subovl(isubo+1)=csc(2)
2011 !        Store real and imag parts in Hermitian storage mode:
2012          subovl(isubo  )=cscre
2013          subovl(isubo+1)=cscim
2014          isubo=isubo+2
2015        end do
2016      else if(istwf_k>=2)then
2017        do ii=1,iband
2018          iwavef=(ii-1)*npw_k*nspinor+icg
2019          if(istwf_k==2 .and. me_g0==1)then
2020            cscre=half*cg(1,1+iwavef)*gsc(1,1+jwavef)
2021            ipw1=2
2022          else
2023            cscre=zero; ipw1=1
2024          end if
2025          do isp=1,nspinor
2026            do ipw=ipw1+(isp-1)*npw_k,npw_k*isp
2027              cgreipw=cg(1,ipw+iwavef)
2028              cgimipw=cg(2,ipw+iwavef)
2029              cscre=cscre+cg(1,ipw+iwavef)*gsc(1,ipw+jwavef)+cg(2,ipw+iwavef)*gsc(2,ipw+jwavef)
2030            end do
2031          end do
2032          cscre=two*cscre
2033 !        Store real and imag parts in Hermitian storage mode:
2034          subovl(isubo  )=cscre
2035          subovl(isubo+1)=zero
2036          isubo=isubo+2
2037        end do
2038      end if
2039    end if
2040 
2041  end do ! iband in a block
2042 
2043 end subroutine mksubham