Logo Search packages:      
Sourcecode: abinit version File versions  Download package

scfcv3.F90

!{\src2tex{textfont=tt}}
!!****f* ABINIT/scfcv3
!! NAME
!! scfcv3
!!
!! FUNCTION
!! Conducts set of passes or overall iterations of preconditioned
!! conjugate gradient algorithm to converge wavefunctions to
!! optimum and optionally to compute mixed derivatives of energy.
!!
!! COPYRIGHT
!! Copyright (C) 1999-2007 ABINIT group (XG, DRH, MB, XW)
!! 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
!!  atindx(natom)=index table for atoms (see scfcv.f)
!!  atindx1(natom)=index table for atoms, inverse of atindx (see scfcv.f)
!!  cg(2,mpw*nspinor*mband*mkmem*nsppol)=pw coefficients of GS wavefunctions at k.
!!  cgq(2,mpw1*nspinor*mband*mkqmem*nsppol)=pw coefficients of GS wavefunctions at k+q.
!!  cplex: if 1, real space 1-order functions on FFT grid are REAL, if 2, COMPLEX
!!  cpus= cpu time limit in seconds
!!  densymop_rf <type(dens_sym_operator_type)>=the density symmetrization
!!   operator (response-function)
!!  doccde_rbz(mband*nkpt_rbz*nsppol)=derivative of occ_rbz wrt the energy
!!  docckqde(mband*nkpt_rbz*nsppol)=derivative of occkq wrt the energy
!!  dtfil <type(datafiles_type)>=variables related to files
!!  dtset <type(dataset_type)>=all input variables for this dataset
!!  eew=2nd derivative of Ewald energy (hartree)
!!  efrhar=Contribution from frozen-wavefunction, hartree energy,
!!           to the second-derivative of total energy.
!!  efrkin=Contribution from frozen-wavefunction, kinetic energy,
!!           to the second-derivative of total energy.
!!  efrloc=Contribution from frozen-wavefunction, local potential,
!!           to the second-derivative of total energy.
!!  efrnl=Contribution from frozen-wavefunction, non-local potential,
!!           to the second-derivative of total energy.
!!  efrx1=Contribution from frozen-wavefunction, xc core correction(1),
!!           to the second-derivative of total energy.
!!  efrx2=Contribution from frozen-wavefunction, xc core correction(2),
!!           to the second-derivative of total energy.
!!  eigenq(mband*nkpt_rbz*nsppol)=GS eigenvalues at k+q (hartree)
!!  eigen0(mband*nkpt_rbz*nsppol)=GS eigenvalues at k (hartree)
!!  eii=2nd derivative of pseudopotential core energy (hartree)
!!  fermie=fermi energy (Hartree)
!!  fform=integer specifying data structure of wavefunction files
!!  hdr <type(hdr_type)>=the header of wf, den and pot files
!!  idir=direction of the current perturbation
!!  indkpt1(nkpt_rbz)=non-symmetrized indices of the k-points
!!  indsy1(4,nsym1,natom)=indirect indexing array for atom labels
!!  ipert=type of the perturbation
!!  irrzon1(nfft**(1-1/nsym1),2,nspden/nsppol)=irreducible zone data for RF symmetries
!!  istwfk_rbz(nkpt_rbz)=input option parameter that describes the storage of wfs
!!  kg(3,mpw*mkmem)=reduced planewave coordinates at k
!!  kg1(3,mpw1*mk1mem)=reduced planewave coordinates at k+q, with RF k points
!!  kpt_rbz(3,nkpt_rbz)=reduced coordinates of k points.
!!  kxc(nfft,nkxc)=exchange and correlation kernel (see rhohxc.f)
!!  mgfft=maximum size of 1D FFTs
!!  mkmem =number of k points which can fit in memory (GS data); 0 if use disk
!!  mkqmem =number of k+q points which can fit in memory (GS data); 0 if use disk
!!  mk1mem =number of k points which can fit in memory (RF data); 0 if use disk
!!  mpert=maximum number of ipert
!!  mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
!!  mpw=maximum dimensioned size of npw for wfs at k.
!!  mpw1=maximum dimensioned size of npw for wfs at k+q (also for 1-order wfs).
!!  nattyp(ntypat)= # atoms of each type.
!!  nband_rbz(nkpt_rbz*nsppol)=number of bands at each RF k point, for each polarization
!!  nfft=(effective) number of FFT grid points (for this processor)
!!  nkpt=number of k points in the full BZ
!!  nkpt_rbz=number of k points in the reduced BZ for this perturbation
!!  nkxc=second dimension of the kxc array.
!!  mpi_enreg=informations about MPI parallelization
!!  npwarr(nkpt_rbz)=number of planewaves in basis at this GS k point
!!  npwar1(nkpt_rbz)=number of planewaves in basis at this RF k+q point
!!  nspden=number of spin-density components
!!  nspinor=number of spinorial components of the wavefunctions
!!  nsym1=number of symmetry elements in space group consistent with perturbation
!!  n3xccc=dimension of xccc3d1 ; 0 if no XC core correction is used
!!   otherwise, cplex*nfft
!!  occkq(mband*nkpt_rbz*nsppol)=occupation number for each band (often 2)
!!   at each k+q point of the reduced Brillouin zone.
!!  occ_rbz(mband*nkpt_rbz*nsppol)=occupation number for each band (often 2)
!!   at each k point of the reduced Brillouin zone.
!!  pertcase=fuill index of the perturbation
!!  phnons1(2,nfft**(1-1/nsym1),nspden/nsppol)=nonsymmorphic transl. phases,
!!   for RF symmetries
!!  ph1d(2,3*(2*mgfft+1)*natom)=one-dimensional structure factor information
!!  prtbbb=if 1, band-by-band decomposition (also dim of d2bbb)
!!  psps <type(pseudopotential_type)>=variables related to pseudopotentials
!!  qphon(3)=reduced coordinates for the phonon wavelength
!!  rhog(2,nfft)=array for Fourier transform of GS electron density
!!  rhor(nfft,nspden)=array for GS electron density in electrons/bohr**3.
!!  rprimd(3,3)=dimensional primitive translations in real space (bohr)
!!  symaf1(nsym1)=anti(ferromagnetic) part of symmetry operations
!!  symrc1(3,3,nsym1)=symmetry operations in reciprocal space
!!  symrl1(3,3,nsym1)=symmetry operations in real space in terms
!!   of primitive translations
!!  timrev=1 if time-reversal preserves the q wavevector; 0 otherwise.
!!  tnons1(3,nsym1)=nonsymmorphic translations for symmetry operations
!!  wffddk=struct info for ddk file
!!  wffnew=struct info for 1WF at exit, if mk1mem=0
!!  wffnow=struct info for 1WF at start, if mk1mem=0
!!  wfftgs=struct info for GS WF at start, if mkmem=0
!!  wfftkq=struct info for k+q GS WF at start, if mkqmem=0
!!  vpsp1(cplex*nfft)=first-order derivative of the ionic potential
!!  vtrial(nfft,nspden)=GS potential (Hartree).
!!  wtk_rbz(nkpt_rbz)=weight for each k point in the reduced Brillouin zone
!!  xccc3d1(cplex*n3xccc)=3D change in core charge density, see n3xccc
!!  xred(3,natom)=reduced dimensionless atomic coordinates
!!  ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k point
!!  ylm1(mpw1*mk1mem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k+q point
!!  ylmgr(mpw*mkmem,3,mpsang*mpsang*useylm)= gradients of real spherical harmonics at k
!!  ylmgr1(mpw1*mk1mem,3,mpsang*mpsang*useylm)= gradients of real spherical harmonics at k+q
!!
!! OUTPUT
!!  blkflg(3,mpert,3,mpert)=flags for each element of the 2DTE (=1 if computed)
!!  d2bbb(2,3,3,mpert,mband,mband*prtbbb)=band by band decomposition of some
!!       second order derivatives
!!  d2lo(2,mpert,3,mpert)=local contributions to the 2DTEs
!!  d2nl(2,mpert,3,mpert)=non-local contributions to the 2DTEs
!!  edocc=correction to 2nd-order total energy coming from changes of occupation
!!  eeig0=0th-order eigenenergies part of 2nd-order total energy
!!  ehart01=inhomogeneous 1st-order Hartree part of 2nd-order total energy
!!    for strain perturbation only (zero otherwise, and not used)
!!  ehart1=1st-order Hartree part of 2nd-order total energy
!!  eigen1(2*mband*mband*nkpt_rbz*nsppol)=array for holding eigenvalues (hartree)
!!  ek0=0th-order kinetic energy part of 2nd-order total energy.
!!  ek1=1st-order kinetic energy part of 2nd-order total energy.
!!  eloc0=0th-order local (psp+vxc+Hart) part of 2nd-order total energy
!!  elpsp1=1st-order local pseudopot. part of 2nd-order total energy.
!!  enl0=0th-order nonlocal pseudopot. part of 2nd-order total energy.
!!  enl1=1st-order nonlocal pseudopot. part of 2nd-order total energy.
!!  exc1=1st-order exchange-correlation part of 2nd-order total energy
!!  etotal=total energy (sum of 7 contributions) (hartree)
!!  resid(mband*nkpt_rbz*nsppol)=residuals for each band over all k points
!!   of the reduced Brillouin zone, and spins
!!  residm=maximum value from resid array (except for nbdbuf highest bands)
!!
!! SIDE EFFECTS
!!  cg1(2,mpw1*nspinor*mband*mk1mem*nsppol)=updated wavefunctions;  if mk1mem<=nkpt_rbz,
!!         these are kept in a disk file, see wffnow and wffnew.
!!  initialized= if 0 the initialization of the RF run is not yet finished
!!  mpi_enreg=informations about MPI parallelization
!!  rhog1(2,nfft)=array for Fourier transform of RF electron density
!!  rhor1(cplex*nfft,nspden)=array for RF electron density in electrons/bohr**3.
!!
!! TODO
!! Should be taken away of the list of arguments : ehart, ek, enl, enxc
!!
!! PARENTS
!!      loper3
!!
!! CHILDREN
!!      appdig,ebp3,eneres3,getcut,hartre,hartrestr,initberry3,int2char4,ioarr
!!      leave_new,leave_test,metric,mkvxc3,mkvxcstr3,newfermie1,newvtr3,nselt3
!!      nstdy3,qmatrix,rhofermi3,scprqt,status,timab,vtorho3,wffclose,wrtout
!!      xme_init
!!
!! SOURCE

#if defined HAVE_CONFIG_H
#include "config.h"
#endif

subroutine scfcv3(atindx,atindx1,blkflg,cg,cgq,cg1,cplex,cpus,&
&  densymop_rf,doccde_rbz,docckqde,dtfil,dtset,&
&  d2bbb,d2lo,d2nl,edocc,eeig0,eew,efrhar,efrkin,efrloc,efrnl,efrx1,efrx2,&
&  ehart,ehart01,ehart1,eigenq,eigen0,eigen1,eii,ek,ek0,ek1,eloc0,elpsp1,&
&  enl,enl0,enl1,enxc,etotal,exc1,fermie,fform,hdr,idir,indkpt1,&
&  indsy1,initialized,ipert,irrzon1,istwfk_rbz,&
&  kg,kg1,kpt_rbz,kxc,mgfft,mkmem,mkqmem,mk1mem,&
&  mpert,mpi_enreg,mpsang,mpw,mpw1,nattyp,nband_rbz,&
&  nfft,nkpt,nkpt_rbz,nkxc,npwarr,npwar1,nspden,nspinor,&
&  nsym1,n3xccc,occkq,occ_rbz,pertcase,phnons1,ph1d,&
&  prtbbb,psps,qphon,resid,residm,rhog,rhog1,&
&  rhor,rhor1,rprimd,symaf1,symrc1,symrl1,timrev,&
&  tnons1,wffddk,wffnew,wffnow,wfftgs,wfftkq,vpsp1,vtrial,&
&  wtk_rbz,xccc3d1,xred,ylm,ylm1,ylmgr,ylmgr1,ngfft)

 use defs_basis
 use defs_datatypes

!This section has been created automatically by the script Abilint (TD). Do not modify these by hand.
#ifdef HAVE_FORTRAN_INTERFACES
 use interfaces_01manage_mpi
 use interfaces_11util
 use interfaces_12geometry
 use interfaces_13io_mpi
 use interfaces_13recipspace
 use interfaces_13xc
 use interfaces_14iowfdenpot
 use interfaces_15common
 use interfaces_16response, except_this_one => scfcv3
 use interfaces_lib01hidempi
#endif
!End of the abilint section

 implicit none

!Arguments ------------------------------------
!no_abirules
!Needed for integer arrays
 type(dataset_type),intent(in) :: dtset
 type(pseudopotential_type),intent(in) :: psps
!---
 integer,intent(in) :: cplex,fform,idir,ipert,mgfft,mk1mem,mkmem,mkqmem
 integer,intent(in) :: mpert,mpsang,mpw,mpw1,n3xccc,nfft
 integer,intent(in) :: nkpt,nkpt_rbz,nkxc,nspden
 integer,intent(in) :: nsym1,pertcase,prtbbb,timrev
 integer,intent(inout) :: initialized,nspinor
! nfft**(1-1/nsym1) is 1 if nsym1==1, and nfft otherwise
 integer,intent(in) :: atindx(dtset%natom),atindx1(dtset%natom),ngfft(18)
 integer,intent(out) :: blkflg(3,mpert,3,mpert)
 integer,intent(in) :: indkpt1(nkpt_rbz),indsy1(4,nsym1,dtset%natom)
 integer,intent(in) :: irrzon1(nfft**(1-1/nsym1),2,nspden/dtset%nsppol)
 integer,intent(in) :: istwfk_rbz(nkpt_rbz)
 integer,intent(in) :: kg(3,mpw*mkmem),kg1(3,mpw1*mk1mem),nattyp(psps%ntypat)
 integer,intent(in) :: nband_rbz(nkpt_rbz*dtset%nsppol)
 integer,intent(in) :: npwar1(nkpt_rbz),npwarr(nkpt_rbz)
 integer,intent(in) :: symaf1(nsym1),symrc1(3,3,nsym1),symrl1(3,3,nsym1)
 real(dp),intent(in) :: cpus,eew,efrhar,efrkin,efrloc,efrnl,efrx1,efrx2,ehart
 real(dp),intent(out) :: edocc,eeig0,ehart01,ehart1,ek0,ek1,eloc0,elpsp1,enl0
 real(dp),intent(out) :: enl1,exc1,etotal,residm
 real(dp),intent(in) :: eii,ek,enl,enxc
 real(dp),intent(inout) :: fermie
 real(dp),intent(in) :: qphon(3)
! nfft**(1-1/nsym1) is 1 if nsym1==1, and nfft otherwise
 real(dp),intent(in) :: cg(2,mpw*nspinor*dtset%mband*mkmem*dtset%nsppol)
 real(dp),intent(inout) :: cg1(2,mpw1*nspinor*dtset%mband*mk1mem*dtset%nsppol)
 real(dp),intent(in) :: cgq(2,mpw1*nspinor*dtset%mband*mkqmem*dtset%nsppol)
 real(dp),intent(out) :: d2bbb(2,3,3,mpert,dtset%mband,dtset%mband*prtbbb)
 real(dp),intent(out) :: d2lo(2,3,mpert,3,mpert)
 real(dp),intent(out) :: d2nl(2,3,mpert,3,mpert)
 real(dp),intent(in) :: doccde_rbz(dtset%mband*nkpt_rbz*dtset%nsppol)
 real(dp),intent(in) :: docckqde(dtset%mband*nkpt_rbz*dtset%nsppol)
 real(dp),intent(in) :: eigen0(dtset%mband*nkpt_rbz*dtset%nsppol)
 real(dp),intent(out) :: eigen1(2*dtset%mband*dtset%mband*nkpt_rbz*dtset%nsppol)
 real(dp),intent(in) :: eigenq(dtset%mband*nkpt_rbz*dtset%nsppol)
 real(dp),intent(in) :: kpt_rbz(3,nkpt_rbz),kxc(nfft,nkxc)
 real(dp),intent(in) :: occ_rbz(dtset%mband*nkpt_rbz*dtset%nsppol)
 real(dp),intent(in) :: occkq(dtset%mband*nkpt_rbz*dtset%nsppol)
 real(dp),intent(in) :: ph1d(2,3*(2*mgfft+1)*dtset%natom)
 real(dp),intent(in) :: phnons1(2,nfft**(1-1/nsym1),nspden/dtset%nsppol)
 real(dp),intent(out) :: resid(dtset%mband*nkpt_rbz*nspden)
 real(dp),intent(in) :: rhog(2,nfft),rhor(nfft,nspden),rprimd(3,3)
 real(dp),intent(inout) :: rhog1(2,nfft),rhor1(cplex*nfft,nspden),vtrial(nfft,nspden),xred(3,dtset%natom)
 real(dp),intent(in) :: tnons1(3,nsym1),vpsp1(cplex*nfft)
 real(dp),intent(in) :: wtk_rbz(nkpt_rbz),xccc3d1(cplex*n3xccc)
 real(dp),intent(in) :: ylm(mpw*mkmem,mpsang*mpsang*psps%useylm)
 real(dp),intent(in) :: ylm1(mpw1*mk1mem,mpsang*mpsang*psps%useylm)
 real(dp),intent(in) :: ylmgr(mpw*mkmem,3,mpsang*mpsang*psps%useylm)
 real(dp),intent(in) :: ylmgr1(mpw1*mk1mem,3,mpsang*mpsang*psps%useylm)
 type(datafiles_type),intent(in) :: dtfil
 type(hdr_type),intent(inout) :: hdr
 type(dens_sym_operator_type),intent(in) :: densymop_rf
 type(MPI_type),intent(inout) :: mpi_enreg
 type(wffile_type),intent(inout) :: wffddk,wfftgs,wfftkq
 type(wffile_type),intent(inout) :: wffnow,wffnew

!Local variables-------------------------------
!scalars
 integer,parameter :: level=16,response=1
 integer :: accessfil,afford,choice,coordn,dbl_nnsclo,dielop,dielstrt,fformr
 integer :: fformv,fftalg,formeis,formeiw,i1,i2,i3,iatom,ierr,iexit,ifft,ii
 integer :: ikpt,index,iprcel,ipw,ipw1,ipw2,ir,iscf_mod,ispden,isppol,istep
 integer :: itypat,me,mgfftdiel,mvdum,n1,n2,n3,n_fftgr,n_index,nfftdiel,nfftot
 integer :: npwdiel,nstep,nsym,option,optxc,prtfor,quit,qzero,rdwr
 integer :: rdwrpaw
 real(dp) :: boxcut,deltae,diecut,diemix,diffor,eberry,ecut,ecutsus,elast
 real(dp) :: etotal_old,evar,fe1fixed,fermie1,gsqcut,maxfor,rho1_dn,rho1im_dn
 real(dp) :: rho1re_dn,ucvol,vres2,vxcavg
 logical :: ex,test_mpi
 character(len=4) :: tag
 character(len=500) :: message
 character(len=fnlen) :: fi1o,filapp,fildata,filfft,filkgs,kgnam
 type(efield_type) :: dtefield
 type(wffile_type) :: wfftmp
!arrays
 integer,allocatable :: i_rhor(:),i_vresid(:),i_vrespc(:),i_vtrial(:)
 real(dp) :: dielar(7),dummy2(6),favg(3),gmet(3,3),gprimd(3,3),k0(3)
 real(dp) :: kpt_diel(3),rmet(3,3),strsxc(6),tollist(12),tsec(2)
 real(dp),allocatable :: dielinv(:,:,:,:,:),dyfrx2(:,:,:),f_fftgr(:,:,:)
 real(dp),allocatable :: fcart(:,:),forold(:,:),grhf(:,:),grtn(:,:)
 real(dp),allocatable :: rhorfermi(:,:),susmat(:,:,:,:,:),vhartr01(:)
 real(dp),allocatable :: vhartr1(:),vresid1(:,:),vt1g(:,:),vtrial1(:,:)
 real(dp),allocatable :: vxc1(:,:),work(:)
 type(pawrhoij_type) :: pawrhoij(1)
!no_abirules
 integer,allocatable :: pwindall(:,:,:)
 real(dp),allocatable ::qmat(:,:,:,:,:,:)

! *********************************************************************

!DEBUG
!write(6,*)' scfcv3 : enter'
!ENDDEBUG

 call timab(120,1,tsec)
 call timab(154,1,tsec)

 call status(0,dtfil%filstat,iexit,level,'enter         ')

!Structured debugging if dtset%prtvol==-level
 if(dtset%prtvol==-level)then
  write(message,'(80a,a,a)') ('=',ii=1,80),ch10,' scfcv3: enter '
  call wrtout(06,message,'COLL')
 end if

!Init me
 call xme_init(mpi_enreg,me)

!
!  If dtset%accesswff == 2 set all array outputs to netcdf format
!
 accessfil = 0
 if (dtset%accesswff == 2) then
   accessfil = 1
 end if
 if (dtset%accesswff == 3) then
   accessfil = 3
 end if

 iprcel=dtset%iprcel
 iscf_mod=dtset%iscf
 formeiw=1
 formeis=0

!The value of iscf must be modified if ddk perturbation, see loper3.f
 if(ipert==dtset%natom+1)iscf_mod=-3

 ecut=dtset%ecut
 tollist(1)=dtset%tolmxf
 tollist(2)=dtset%tolwfr
 tollist(3)=dtset%toldff
 tollist(4)=dtset%toldfe
 tollist(6)=dtset%tolvrs
 tollist(7)=dtset%tolrff

 nstep=dtset%nstep

 quit=0 ; dbl_nnsclo=0 ; elast=zero
 n1=dtset%ngfft(1) ; n2=dtset%ngfft(2) ; n3=dtset%ngfft(3)

!This might be taken away later
 edocc=zero ; eeig0=zero ; ehart01=zero ; ehart1=zero ; ek0=zero ; ek1=zero
 eloc0=zero ; elpsp1=zero ; enl0=zero ; enl1=zero ; exc1=zero
 deltae=zero
 fermie1=zero

!Prepare the name of the _FFT file and the _KGS file
 filfft=trim(dtfil%filnam_ds(5))//'_FFT'
 filkgs=trim(dtfil%filnam_ds(5))//'_KGS'

 filapp=trim(dtfil%filnam_ds(5))

 if(mpi_enreg%paral_compil_kpt==1 .or. mpi_enreg%paral_compil_fft==1)then
!BEGIN TF_CHANGES
 call int2char4(me,tag)
!END TF_CHANGES
  filfft=trim(filfft)//'_P-'//tag
  filkgs=trim(filkgs)//'_P-'//tag
 end if

 call status(0,dtfil%filstat,iexit,level,'allocate      ')

 allocate(fcart(3,dtset%natom),forold(3,dtset%natom))
 allocate(vhartr1(cplex*nfft),vtrial1(cplex*nfft,nspden))

!Examine tolerance criteria, and eventually  print a line to the output
!file (with choice=1, the only non-dummy arguments of scprqt are
!nstep, tollist and iscf - still, diffor and vres2 are here initialized to 0)
 choice=1 ; prtfor=0
 diffor=0.0_dp ; vres2=0.0_dp
 call scprqt(choice,cpus,deltae,diffor,dtset,eigen0,&
&  etotal,favg,fcart,fermie,filapp,dtfil%filnam_ds(1),&
&  iscf_mod,istep,kpt_rbz,maxfor,&
&  mvdum,mpi_enreg,nband_rbz,nkpt_rbz,&
&  nstep,occ_rbz,0,prtfor,&
&  quit,vres2,resid,residm,response,&
&  tollist,psps%usepaw,vxcavg,wtk_rbz,xred)

 forold(:,:)=0.0_dp
 dielop=0
!Compute different geometric tensor, as well as ucvol, from rprimd
 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)

!Compute large sphere cut-off gsqcut
 k0(:)=0.0_dp
 call getcut(boxcut,ecut,gmet,gsqcut,dtset%iboxcut,6,k0,dtset%ngfft)

!These arrays are needed only in the self-consistent case
 if(nstep>0 .and. iscf_mod>0) then

  allocate(grhf(3,dtset%natom))

  if(iscf_mod==1) then

!  For iscf_mod==1, five additional vectors are needed
   n_fftgr=5 ; n_index=1
   allocate(i_rhor(n_index),i_vtrial(n_index),i_vresid(n_index),i_vrespc(n_index))
!  The index 1 is attributed to the old trial potential,
!  The new residual potential, and the new
!  preconditioned residual potential receive now a temporary index
   i_vtrial(1)=1 ; i_vresid(1)=2 ; i_vrespc(1)=3
!  The indices number 4 and 5 are attributed to work vectors.

  else if((iscf_mod==2).or.(iscf_mod==3)) then

!  For iscf_mod==2 or 3 , four additional vectors are needed.
   n_fftgr=4 ; n_index=2
   allocate(i_rhor(n_index),i_vtrial(n_index),i_vresid(n_index),i_vrespc(n_index))
!  The index number 1 is attributed to the old trial vector
!  The new residual potential, and the new and old preconditioned
!  residual potential, receive now a temporary index.
   i_vtrial(1)=1 ; i_vresid(1)=2 ; i_vrespc(1)=3 ; i_vrespc(2)=4

  else if((iscf_mod==5).or.(iscf_mod==6)) then

!  For iscf_mod==5 or 6, ten additional vectors are needed
   n_fftgr=10 ; n_index=3
   allocate(i_rhor(n_index),i_vtrial(n_index),i_vresid(n_index),i_vrespc(n_index))
!  The index number 6 is attributed to the search vector
!  The index number 1 is attributed to the old trial vector
!  Other indices are attributed now.
   i_vtrial(1)=1
   i_vresid(1)=2 ; i_vrespc(1)=3
   i_vresid(2)=4 ; i_vrespc(2)=5
   i_vresid(3)=7 ; i_vrespc(3)=8 ; i_rhor(2)=9 ; i_rhor(3)=10

  else if(iscf_mod==7) then

!  For iscf_mod==7, lot of additional vectors are needed
   n_fftgr=2+2*dtset%npulayit ; n_index=1+dtset%npulayit
   allocate(i_rhor(n_index),i_vtrial(n_index),i_vresid(n_index),i_vrespc(n_index))
!  The index number 1 is attributed to the old trial vector
!  The index number 2 is attributed to the old residual
!  The indices number 2 and 3 are attributed to two old precond. residuals
!  Other indices are attributed now.
   i_vrespc(dtset%npulayit+1)=2*dtset%npulayit+1; i_vresid(1)=2*dtset%npulayit+2
   do ii=1,dtset%npulayit
    i_vtrial(ii)=2*ii-1 ; i_vrespc(ii)=2*ii
   end do

  end if

! The big array containing functions defined on the fft grid is allocated now.
! Note however, that a zero value of mffmem will cause allocation with
! the third dimension set to zero ! In this case, another temporary
! will be used inside newvtr3.
!
  allocate(f_fftgr(cplex*nfft,nspden,n_fftgr*dtset%mffmem))

 end if

! Here, allocate arrays for computation of susceptibility and dielectric matrix
! or for TDDFT
 if( (nstep>0 .and. iscf_mod>0) .or. iscf_mod==-1 ) then

! Here, for TDDFT, artificially set iprcel . Also set a variable to reduce
! the memory needs.
  afford=1
  if(iscf_mod==-1) then
   iprcel=21
   afford=0
  end if

  npwdiel=1
  mgfftdiel=1
  nfftdiel=1

! Now, performs allocation
! CAUTION : the dimensions are still those of GS, except for phnonsdiel
  allocate(dielinv(2,npwdiel*afford,nspden,npwdiel,nspden))
  allocate(susmat(2,npwdiel*afford,nspden,npwdiel,nspden))

 end if

 call timab(154,2,tsec)


!Initialize Berry-phase related stuffs

if (dtset%berryopt == 4) then

 allocate(pwindall(max(mpw,mpw1)*mkmem,8,3))
 call initberry3(dtefield,dtfil,dtset,gmet,kg,kg1,dtset%mband,mkmem,mpi_enreg,&
&                mpw,mpw1,nkpt,npwarr,npwar1,dtset%nsppol,occ_rbz,pwindall,rprimd)


!calculate inverse of the overlap matrix
 allocate(qmat(2,dtefield%nband_occ,dtefield%nband_occ,nkpt,2,3))
 call qmatrix(cg,dtefield,qmat,mpw,mpw1,mkmem,dtset%mband,npwarr,nkpt,nspinor,dtset%nsppol,pwindall)

end if



!#####################################################################

!PERFORM ITERATIONS

!THIS TO BE UPDATED :
!Offer option of computing total energy with existing
!wavefunctions when nstep<=0, else do nstep iterations
!Note that for non-self-consistent calculations, this loop will be exited
!after the first call to vtorho

!Pass through the first routines even when nstep==0
 do istep=1,max(1,nstep)

  if (istep==1)then

   call status(istep,dtfil%filstat,iexit,level,'get vtrial1   ')

   call hartre(cplex,gmet,gsqcut,psps%usepaw,mpi_enreg,nfft,dtset%ngfft,dtset%qptn,rhog1,vhartr1)

!  section for hartree strain perturbation initialization
   if(ipert==dtset%natom+3 .or. ipert==dtset%natom+4)then
    allocate(vhartr01(cplex*nfft))
    call hartrestr(gmet,gprimd,gsqcut,idir,ipert,mpi_enreg,dtset%natom,nfft,dtset%ngfft,&
&    rhog,vhartr01)
    vhartr1(:)=vhartr1(:)+vhartr01(:)
    deallocate(vhartr01)
   end if

   allocate(vxc1(cplex*nfft,nspden))
   option=1

   if(ipert==dtset%natom+3 .or. ipert==dtset%natom+4) then
!   option for vxc strain perturbation
    call mkvxcstr3(cplex,gmet,gsqcut,idir,ipert,kxc,mpi_enreg,dtset%natom,nfft,dtset%ngfft,&
&    nkxc,nspden,n3xccc,option,dtset%qptn,rhor,rhor1,rprimd,vxc1,xccc3d1)
   else
    call mkvxc3(cplex,gmet,gsqcut,kxc,mpi_enreg,nfft,dtset%ngfft,nkxc,&
&    nspden,n3xccc,option,dtset%qptn,rhor1,rprimd,vxc1,xccc3d1)
   end if
   if(dtset%nsppol==1)then
    if(cplex==1)then
     do ir=1,nfft
      vtrial1(ir,1)=vpsp1(ir)+vhartr1(ir)+vxc1(ir,1)
     end do
    else
     do ir=1,nfft
      vtrial1(2*ir-1,1)=vpsp1(2*ir-1)+vhartr1(2*ir-1)+vxc1(2*ir-1,1)
      vtrial1(2*ir  ,1)=vpsp1(2*ir  )+vhartr1(2*ir  )+vxc1(2*ir  ,1)
     end do
    end if
   else
    if(cplex==1)then
     do ir=1,nfft
      vtrial1(ir,1)=vpsp1(ir)+vhartr1(ir)+vxc1(ir,1)
      vtrial1(ir,2)=vpsp1(ir)+vhartr1(ir)+vxc1(ir,2)
     end do
    else
     do ir=1,nfft
      vtrial1(2*ir-1,1)=vpsp1(2*ir-1)+vhartr1(2*ir-1)+vxc1(2*ir-1,1)
      vtrial1(2*ir  ,1)=vpsp1(2*ir  )+vhartr1(2*ir  )+vxc1(2*ir  ,1)
      vtrial1(2*ir-1,2)=vpsp1(2*ir-1)+vhartr1(2*ir-1)+vxc1(2*ir-1,2)
      vtrial1(2*ir  ,2)=vpsp1(2*ir  )+vhartr1(2*ir  )+vxc1(2*ir  ,2)
     end do
    end if
   end if
   deallocate(vxc1)

!  For Q=0 and metallic occupation, initialize quantities needed to
!   compute the first-order Fermi energy
   qzero=0
   if(qphon(1)**2+qphon(2)**2+qphon(3)**2 < 1.d-14) qzero=1

   if(qzero==1 .and. (dtset%occopt>=3 .and. dtset%occopt <=7) .and.&
&    (ipert<=dtset%natom .or. ipert==dtset%natom+3 .or. ipert==dtset%natom+4) .and. dtset%frzfermi==0) then
    allocate(rhorfermi(cplex*nfft,nspden))

    call rhofermi3(atindx,atindx1,cg,cgq,cplex,densymop_rf,&
&    doccde_rbz,docckqde,dtfil,dtset,&
&    edocc,eeig0,eigenq,eigen0,eigen1,&
&    fe1fixed,gmet,gprimd,gsqcut,hdr,idir,&
&    ipert,irrzon1,istwfk_rbz,kg,kg1,kpt_rbz,dtset%mband,mgfft,&
&    mkmem,mkqmem,mk1mem,mpi_enreg,mpsang,mpw,mpw1,&
&    dtset%natom,nattyp,nband_rbz,nfft,nkpt_rbz,npwarr,npwar1,nspden,nspinor,&
&    dtset%nsppol,nsym1,dtset%ntypat,occkq,occ_rbz,phnons1,&
&    ph1d,dtset%prtvol,psps,rhorfermi,rmet,symaf1,tnons1,ucvol,&
&    wfftgs,wfftkq,wtk_rbz,xred,ylm,ylm1,ylmgr1)

   end if !First-order fermi energy setup

!DEBUG
!  allocate(vt1g(2,nfft))
!   call symrhg(cplex,densymop_rf,irrzon1,nfft,dtset%ngfft,nspden,dtset%nsppol,&
!&   nsym1,phnons1,vt1g,vtrial1,symaf1)
!   call symrhg(cplex,densymop_rf,irrzon1,nfft,dtset%ngfft,nspden,dtset%nsppol,&
!&   nsym1,phnons1,vt1g,vtrial1,symaf1)
!  deallocate(vt1g)
!  write(6,*)' after symrhg(init)'
!  do ifft=1,nfft
!   i3=(ifft-1)/n1/n2
!   i2=(ifft-1-i3*n1*n2)/n1
!   i1=ifft-1-i3*n1*n2-i2*n1
!   if(i1==0 .and. i2==0)write(6,*)i1,i2,i3,vtrial1(ifft,1)
!  end do

!ENDDEBUG

! End the condition of istep==1
  end if

! No need to continue and call vtorho, when nstep==0
  if(nstep==0)exit

  call status(istep,dtfil%filstat,iexit,level,'loop istep    ')

  if(iscf_mod>0)then
   write(message, '(a,a,i4)' )ch10,' ITER STEP NUMBER  ',istep
   call wrtout(06,message,'COLL')
  end if

  call status(istep,dtfil%filstat,iexit,level,'call vtorho3  ')

!DEBUG
! write(6,*)' scfcv3 : before vtorho3, vtrial1(1,1)=',vtrial1(1,1)
!ENDDEBUG

  if(qzero==1 .and. (dtset%occopt>=3 .and. dtset%occopt <=7) .and.&
&   (ipert<=dtset%natom .or. ipert==dtset%natom+3 .or. ipert==dtset%natom+4) .and. dtset%frzfermi==0) then

!  For Q=0 and metallic occupation, calculate the first-order Fermi energy
   nfftot=n1*n2*n3
   call newfermie1(cplex,fermie1,fe1fixed,istep,&
&   mpi_enreg,nfft,nfftot,nspden,dtset%occopt,&
&    dtset%prtvol,rhorfermi,ucvol,vtrial1)
  end if

!#######################################################################
! Compute the density rhor1 (and rhog1) from the trial potential vtrial1

  call vtorho3(atindx,atindx1,cg,cgq,cg1,cplex,cpus,dbl_nnsclo,densymop_rf,&
&  doccde_rbz,docckqde,dtefield,&
&  dtfil,dtset,edocc,eeig0,eigenq,eigen0,eigen1,ek0,ek1,eloc0,&
&  enl0,enl1,fermie1,gmet,gprimd,gsqcut,hdr,idir,&
&  ipert,irrzon1,istep,istwfk_rbz,kg,kg1,kpt_rbz,dtset%mband,mgfft,&
&  mkmem,mkqmem,mk1mem,mpi_enreg,mpsang,mpw,mpw1,&
&  dtset%natom,nattyp,nband_rbz,nfft,nkpt_rbz,npwarr,npwar1,nspden,nspinor,&
&  dtset%nsppol,nsym1,dtset%ntypat,occkq,occ_rbz,phnons1,&
&  ph1d,dtset%prtvol,psps,pwindall,qmat,resid,residm,rhog1,rhor1,rmet,symaf1,tnons1,ucvol,&
&  wffddk,wffnew,wffnow,wfftgs,wfftkq,vtrial,vtrial1,wtk_rbz,xred,ylm,ylm1,ylmgr1)


!#######################################################################
  call status(istep,dtfil%filstat,iexit,level,'after vtorho3 ')

  if (dtset%berryopt == 4) then
!  calculate \Omega E \cdot P term 
   call  ebp3(cg,cg1,dtefield,eberry,dtset%mband,mkmem,&
&          mpw,mpw1,nkpt,npwarr,npwar1,dtset%nsppol,nspinor,pwindall,qmat)
  endif

! Rotate labels of disk files when wf i/o is used
  if (mk1mem==0) then
   wfftmp=wffnow ; wffnow=wffnew ; wffnew=wfftmp
  end if

! Skip out of step loop if non-SCF (completed)
! Indeed, nstep loops have been done inside vtorho
  if (iscf_mod<=0 .and. iscf_mod/=-3) exit

  allocate(vresid1(cplex*nfft,nspden))

!#######################################################################
! Compute the energy and the potential residual
! from the trial potential and the corresponding density.
  call eneres3(cplex,deltae,edocc,eeig0,ehart01,ehart1,ek0,ek1,&
&  elast,eloc0,elpsp1,&
&  enl0,enl1,exc1,evar,dtset%frzfermi,gsqcut,idir,ipert,kxc,&
&  mpi_enreg,dtset%natom,nfft,dtset%ngfft,&
&  nkxc,nspden,n3xccc,dtset%occopt,dtset%qptn,rhog,rhog1,rhor,rhor1,&
&  rprimd,dtset%tsmear,vhartr1,vpsp1,vresid1,vres2,vtrial1,xccc3d1)

if (dtset%berryopt == 4) then
  etotal=evar+eew+eii+efrhar+efrkin+efrloc+efrnl+efrx1+efrx2+2_dp*eberry
else
  etotal=evar+eew+eii+efrhar+efrkin+efrloc+efrnl+efrx1+efrx2
end if


!DEBUG
!  call prtene3(edocc,eeig0,eew,efrloc,efrnl,efrx1,efrx2,&
!&   ehart1,ek0,ek1,eloc0,elpsp1,enl0,enl1,exc1,6,ipert,dtset%natom)
!ENDDEBUG


!#######################################################################

!------Compute and print additional information, also check exit criteria.

  call timab(152,1,tsec)

! The next flag says whether the xred have to be changed in the current iteration

  call status(istep,dtfil%filstat,iexit,level,'print info    ')
  choice=2
  call scprqt(choice,cpus,deltae,diffor,dtset,eigen0,&
&  etotal,favg,fcart,fermie,filapp,dtfil%filnam_ds(1),&
&  iscf_mod,istep,kpt_rbz,maxfor,&
&  mvdum,mpi_enreg,nband_rbz,nkpt_rbz,&
&  nstep,occ_rbz,0,prtfor,&
&  quit,vres2,resid,residm,response,&
&  tollist,psps%usepaw,vxcavg,wtk_rbz,xred)

  call timab(152,2,tsec)

! If criteria in scprqt say to quit, then exit the loop over istep.
  if (quit==1) then
   deallocate(vresid1)
   exit
  end if

!#######################################################################
!------ Precondition the residual and forces, then determine the new vtrial
! (Warning : the (H)xc potential may have been subtracted from vtrial)


!DEBUG
!  write(6,*)' scfcv3 : before newvtr3, vtrial1'
!  do ifft=1,nfft
!   i3=(ifft-1)/n1/n2
!   i2=(ifft-1-i3*n1*n2)/n1
!   i1=ifft-1-i3*n1*n2-i2*n1
!   if(i1==0 .and. i2==0)write(6,*)i1,i2,i3,vtrial1(ifft,1)
!  end do
!  write(6,*)' scfcv3 : before newvtr3, vresid1'
!  do ifft=1,nfft
!   i3=(ifft-1)/n1/n2
!   i2=(ifft-1-i3*n1*n2)/n1
!   i1=ifft-1-i3*n1*n2-i2*n1
!   if(i1==0 .and. i2==0)write(6,*)i1,i2,i3,vresid1(ifft,1)
!  end do
!ENDDEBUG

  if(iscf_mod/=-3)then

!  Note that vresid1 and vtrial1 are called vresid and vtrial inside this routine
   dielar(1)=dtset%diecut
   dielar(2)=dtset%dielng
   dielar(3)=dtset%diemac
   dielar(4)=dtset%diemix
   dielar(5)=dtset%diegap
   dielar(6)=dtset%dielam
   call newvtr3(cplex,dbl_nnsclo,dielar,etotal,filfft,f_fftgr,&
&   initialized,iscf_mod,dtset%isecur,istep,i_rhor,&
&   i_vresid,i_vrespc,i_vtrial,dtset%mffmem,mpi_enreg,dtset%natom,nfft,dtset%ngfft,nspden,&
&   n_fftgr,n_index,pawrhoij,dtset%qptn,rhor1,rprimd,vresid1,vtrial1,xred)

!DEBUG
!  write(6,*)' scfcv3 : after newvtr3'
!  qphon(:)=zero
!  call filterpot(cplex,gmet,gsqcut,nfft,dtset%ngfft,nspden,qphon,vtrial1)
!  write(6,*)' scfcv3 : after newvtr3'

!  do ifft=1,nfft
!   i3=(ifft-1)/n1/n2
!   i2=(ifft-1-i3*n1*n2)/n1
!   i1=ifft-1-i3*n1*n2-i2*n1
!   if(i1==0 .and. i2==0)write(6,*)i1,i2,i3,vtrial1(ifft,1)
!  end do

!   allocate(vt1g(2,nfft))
!   call symrhg(cplex,densymop_rf,irrzon1,nfft,dtset%ngfft,nspden,dtset%nsppol,&
!&   nsym1,phnons1,rhog1,rhor1,symaf1)
!   deallocate(vt1g)

!   write(6,*)' scfcv3 : after symrhg3'
!   do ifft=1,nfft
!    i3=(ifft-1)/n1/n2
!    i2=(ifft-1-i3*n1*n2)/n1
!    i1=ifft-1-i3*n1*n2-i2*n1
!    if(i1==0 .and. i2==0)write(6,*)i1,i2,i3,vtrial1(ifft,1)
!   end do
!   stop

!ENDDEBUG

!  The initialisation of the RF run should be done when this point is reached
   initialized=1
  end if

!#######################################################################

  deallocate(vresid1)

! End minimization step loop here. Note that there are different "exit"
! instructions within the loop
 end do ! istep

!#####################################################################

! if(nstep==0) then (not implemented)
! end if

!#####################################################################

!DEBUG
!write(6,*)' scfcv3 : after istep loop, continue'
!stop
!ENDDEBUG

 call timab(160,1,tsec)
 call status(0,dtfil%filstat,iexit,level,'endloop istep ')

!Delete eventual _FFT file
 if(dtset%mffmem==0)then
  inquire (file=filfft,exist=ex)
  if(ex)then
   open(unit=tmp_unit,file=filfft,form='unformatted',status='old')
   close(unit=tmp_unit,status='DELETE')
  end if
 end if

!Eventually close the dot file, before calling nstdy3
 if(ipert==dtset%natom+2 .and. sum( (dtset%qptn(1:3)) **2 ) <= 1.0d-7 )then
  call WffClose(wffddk,ierr)
 end if

 call timab(160,2,tsec)
 call timab(147,1,tsec)

 if(ipert==dtset%natom+3 .or. ipert==dtset%natom+4) then
  call status(0,dtfil%filstat,iexit,level,'enter nselt3  ')
  call nselt3(atindx,atindx1,blkflg,cg,cg1,cplex,doccde_rbz,docckqde,&
&  d2bbb,d2lo,d2nl,ecut,dtset%ecutsm,dtset%effmass,eigen0,eigen1,fform,&
&  gmet,gprimd,gsqcut,idir,indkpt1,indsy1,&
&  ipert,iscf_mod,istep,istwfk_rbz,kg,kg1,kpt_rbz,kxc,dtset%mband,mgfft,&
&  mkmem,mk1mem,mpert,mpi_enreg,mpsang,mpw,mpw1,&
&  dtset%natom,nattyp,dtset%nband,nband_rbz,nfft,dtset%ngfft,&
&  nkpt,nkpt_rbz,nkxc,dtset%nline,dtset%nloalg,&
&  npwarr,npwar1,nspden,nspinor,dtset%nsppol,&
&  nsym1,dtset%ntypat,occkq,dtset%occopt,occ_rbz,dtset%ortalg,&
&  ph1d,dtset%prtbbb,dtset%prtvol,psps,dtset%qptn,rhog,rhog1,&
&  rhor,rhor1,rmet,rprimd,symrc1,dtset%tsmear,dtset%typat,ucvol,&
&  dtfil%unkg,dtfil%unkg1,&
&  wffnow,wfftgs,dtfil%unylm,dtfil%unylm1,&
&  dtfil%fnamewffddk,wtk_rbz,xred,ylm,ylm1,ylmgr,ylmgr1)
 end if
 if(ipert<=dtset%natom+4)then
  call status(0,dtfil%filstat,iexit,level,'enter nstdy3  ')
  call nstdy3(atindx,atindx1,blkflg,cg,cg1,cplex,doccde_rbz,docckqde,&
&  dtset,d2bbb,d2lo,d2nl,ecut,dtset%ecutsm,eigen0,eigen1,fform,&
&  gmet,gprimd,gsqcut,idir,indkpt1,indsy1,&
&  ipert,iscf_mod,istep,istwfk_rbz,kg,kg1,kpt_rbz,kxc,dtset%mband,mgfft,&
&  mkmem,mk1mem,mpert,mpi_enreg,mpsang,mpw,mpw1,&
&  dtset%natom,nattyp,dtset%nband,nband_rbz,nfft,dtset%ngfft,&
&  nkpt,nkpt_rbz,nkxc,dtset%nline,dtset%nloalg,&
&  npwarr,npwar1,nspden,nspinor,dtset%nsppol,&
&  nsym1,dtset%ntypat,occkq,dtset%occopt,occ_rbz,dtset%ortalg,&
&  ph1d,dtset%prtbbb,dtset%prtvol,psps,dtset%qptn,rhog1,&
&  rhor1,rmet,rprimd,symrc1,dtset%tsmear,dtset%typat,ucvol,&
&  dtfil%unkg,dtfil%unkg1,&
&  wffnow,wfftgs,dtfil%unylm,dtfil%unylm1,&
&  dtfil%fnamewffddk,wtk_rbz,xred,ylm,ylm1)
 end if
 call timab(147,2,tsec)
 call timab(160,1,tsec)

!If SCF convergence was not reached (for nstep>0),
!print a warning to the output file (non-dummy arguments: nstep,
!residm, diffor - infos from tollist have been saved inside )
 call status(0,dtfil%filstat,iexit,level,'call scprqt(en')
 choice=3
 call scprqt(choice,cpus,deltae,diffor,dtset,eigen0,&
& etotal,favg,fcart,fermie,filapp,dtfil%filnam_ds(1),&
& iscf_mod,istep,kpt_rbz,maxfor,&
& mvdum,mpi_enreg,nband_rbz,nkpt_rbz,&
& nstep,occ_rbz,0,prtfor,&
& quit,vres2,resid,residm,response,&
& tollist,psps%usepaw,vxcavg,wtk_rbz,xred)

!Optionally provide output of charge density and/or potential in real space,
!as well as analysis of geometrical factors (bond lengths and bond angles).
!Warnings :
! - core charge is excluded from the charge density;
! - the potential is the INPUT vtrial.
 if(me==0) then
  if (dtset%prtden>0) then
   call status(0,dtfil%filstat,iexit,level,'call ioarr-den')
   rdwr=2 ; fformr=52 ; rdwrpaw=0
   fildata=trim(dtfil%filnam_ds(4))//'_DEN'
   call appdig(pertcase,fildata,fi1o)
   call ioarr(accessfil,rhor1,etotal,fformr,fi1o,hdr,cplex*nfft,nspden,rdwr,rdwrpaw,ngfft)
  end if
  if (dtset%prtpot>0) then
   call status(0,dtfil%filstat,iexit,level,'call ioarr-pot')
   rdwr=2 ; fformv=102 ; rdwrpaw=0
   fildata=trim(dtfil%filnam_ds(4))//'_POT'
   call appdig(pertcase,fildata,fi1o)
   call ioarr(accessfil,vtrial1,fermie,fformv,fi1o,hdr,cplex*nfft,nspden,rdwr,rdwrpaw,ngfft)
  end if
 end if
 if(mpi_enreg%paral_compil_kpt==1 .or. mpi_enreg%paral_compil_fft==1)then
  call timab(61,1,tsec)
!BEGIN TF_CHANGES
  call leave_test(mpi_enreg)
!END TF_CHANGES
  call timab(61,2,tsec)
 end if

!DEBUG
! ii=0
! do ir=1,nfft
!  if(rhog1(1,ir)**2+rhog1(2,ir)**2>1.0D-15) then
!   ii=ii+1
!   if(ii>100) cycle
!   if(ipert==dtset%natom+4) then
!    write(6,'(a,i6,2d24.12)') 'scfcv3:rhog1',ir,rhog1(1,ir),rhog1(2,ir)
!   elseif(ipert==dtset%natom+3) then
!    write(6,'(a,i6,2d24.12)') 'scfcv3:rhog1',ir,rhog1(1,ir)-rhog(1,ir),&
!     rhog1(2,ir)-rhog(2,ir)
!   end if
!  end if
! end do
!ENDDEBUG


! Debugging : print the different parts of rhor1
! MPIWF Warning : this should not be parallelized over space, leave this debugging feature as such.
  if(dtset%prtvol==-level)then
   write(message,'(a)') '   ir       rhor(ir)     '
   call wrtout(06,message,'COLL')
   do ir=1,nfft
    if(ir<=11 .or. mod(ir,301)==0 )then
     write(message,'(i5,a,es13.6)')ir,' ',rhor(ir,1)
     call wrtout(06,message,'COLL')
     if(nspden==2)then
      write(message,'(a,es13.6)')'      ',rhor(ir,2)
      call wrtout(06,message,'COLL')
     end if
    end if
   end do
  end if

!Deallocate the arrays
 deallocate(fcart,forold,vtrial1,vhartr1)
 if(nstep>0 .and. iscf_mod>0) then
  deallocate(i_rhor,i_vtrial,i_vresid,i_vrespc)
  deallocate(f_fftgr,grhf)
 end if
 if( (nstep>0 .and. iscf_mod>0) .or. iscf_mod==-1 ) then
  deallocate(dielinv)
  deallocate(susmat) !! added by MM
 end if

if (dtset%berryopt == 4) then
 deallocate(pwindall,qmat)
 deallocate(dtefield%ikpt_dk,dtefield%cgindex,dtefield%idxkstr,dtefield%kgindex)
end if


 if(allocated(rhorfermi)) deallocate(rhorfermi)

!Structured debugging : if dtset%prtvol=-level, stop here.
 if(dtset%prtvol==-level)then
  write(message,'(a1,a,a1,a,i2,a)') ch10,' scfcv3: exit ',&
&   ch10,'  prtvol=-',level,', debugging mode => stop '
  call wrtout(06,message,'COLL')
  call leave_new('COLL')
 end if

 call status(0,dtfil%filstat,iexit,level,'exit          ')

 call timab(160,2,tsec)
 call timab(120,2,tsec)

!DEBUG
!write(6,*)' scfcv3: exit '
!ENDDEBUG

end subroutine scfcv3
!!***

Generated by  Doxygen 1.6.0   Back to index