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

gstate.F90

!{\src2tex{textfont=tt}}
!!****f* ABINIT/gstate
!! NAME
!! gstate
!!
!! FUNCTION
!! Primary routine for conducting DFT calculations by CG minimization.
!!
!! COPYRIGHT
!! Copyright (C) 1998-2007 ABINIT group (DCA, XG, GMR, JYR, MKV, MT, FJ, MB)
!! 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
!!  codvsn=code version
!!  cpui=initial CPU time
!!  nspinor=number of spinorial components of the wavefunctions
!!  walli=initial wall clock time
!!
!! OUTPUT
!!  npwtot(nkpt) = total number of plane waves at each k point
!!  results_gs <type(results_gs_type)>=results (energy and its components,
!!   forces and its components, the stress tensor) of a ground-state computation
!!
!! SIDE EFFECTS
!!  acell(3)=unit cell length scales (bohr)
!!  dtfil <type(datafiles_type)>=variables related to files
!!  dtset <type(dataset_type)>=all input variables in this dataset
!!   | mband =maximum number of bands (IN)
!!   | mgfft =maximum single fft dimension (IN)
!!   | mkmem =maximum number of k points which can fit in core memory (IN)
!!   | mpw   =maximum number of planewaves in basis sphere (large number) (IN)
!!   | natom =number of atoms in unit cell (IN)
!!   | nfft  =(effective) number of FFT grid points (for this processor) (IN)
!!   | nkpt  =number of k points (IN)
!!   | nspden=number of spin-density components (IN)
!!   | nsppol=number of channels for spin-polarization (1 or 2) (IN)
!!   | nsym  =number of symmetry elements in space group
!!  iexit= exit flag
!!  mpi_enreg=MPI-parallelisation information (some already initialized,
!!   some others to be initialized here)
!!  occ(mband*nkpt*nsppol) = occupation number for each band and k
!!  pawang <type(pawang_type)>=paw angular mesh and related data
!!  pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data
!!  pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
!!  psps <type(pseudopotential_type)>=variables related to pseudopotentials
!!   Before entering the first time in gstate, a significant part of
!!   psps has been initialized :
!!   the integers dimekb,lmnmax,lnmax,mpssang,mpssoang,mpsso,mgrid,
!!     ntypat,n1xccc,usepaw,useylm, and the arrays dimensioned to npsp
!!   All the remaining components of psps are to be initialized in the call
!!   to pspini .
!!   The next time the code enters gstate, psps might be identical to the
!!   one of the previous dtset, in which case, no reinitialisation is scheduled
!!   in pspini.f .
!!  rprim(3,3)=dimensionless real space primitive translations
!!  vel(3,natom)=value of velocity
!!  xred(3,natom) = reduced atomic coordinates
!!
!! NOTES
!! USE OF FFT GRIDS:
!! =================
!! In case of PAW:
!! ---------------
!!    Two FFT grids are used:
!!    - A "coarse" FFT grid (defined by ecut)
!!      for the application of the Hamiltonian on the plane waves basis.
!!      It is defined by nfft, ngfft, mgfft, ...
!!      Hamiltonian, wave-functions, density related to WFs (rhor here), ...
!!      are expressed on this grid.
!!    - A "fine" FFT grid (defined) by ecutdg)
!!      for the computation of the density inside PAW spheres.
!!      It is defined by nfftf, ngfftf, mgfftf, ...
!!      Total density, potentials, ...
!!      are expressed on this grid.
!! In case of norm-conserving:
!! ---------------------------
!!    - Only the usual FFT grid (defined by ecut) is used.
!!      It is defined by nfft, ngfft, mgfft, ...
!!      For compatibility reasons, (nfftf,ngfftf,mgfftf)
!!      are set equal to (nfft,ngfft,mgfft) in that case.
!!
!! TODO
!! Not yet possible to use restartxf in parallel when localrdwf==0
!!
!! PARENTS
!!      driver,pstate
!!
!! CHILDREN
!!      blok8,brdmin,bstruct_clean,bstruct_init,chkexi,clnmpi_fft,clnmpi_gs
!!      clnup1,clnup2,delocint,diisrelax,energies_init,fconv,fixsym,fourdp
!!      getph,handle_ncerr,hdr_clean,hdr_init,hdr_update,indgrid,initberry
!!      initmpi_fft,initmpi_gs,initrhoij,initro,initylmg,int2char4,inwffil
!!      ioarr,ioddb8,kpgio,leave_new,mkrho,moldyn,move,mpi_comm_size,newocc
!!      outqmc,outwf,outxfhist,pawinit,pawpuinit,prtene,psddb8,pspini,scfcv
!!      setsym,setsymrhoij,setup1,setup2,status,timab,transgrid,wffclose
!!      wffdelete,wffopen,wffreadskiprec,wrtout,wvl_free_type_proj
!!      wvl_free_type_wfs,wvl_init_type_proj,wvl_init_type_wfs,wvl_mkrho
!!      wvl_setboxgeometry,xcomm_world,xme_init,xproc_max
!!
!! SOURCE

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

subroutine gstate(acell,codvsn,cpui,dtfil,dtset,iexit,&
& mpi_enreg,&
& npwtot,nspinor,&
& occ,pawang,pawrad,pawtab,psps,results_gs,rprim,vel,walli,xred)

 use defs_basis
 use defs_datatypes
#if defined HAVE_NETCDF
 use netcdf
#endif

!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_12ffts
 use interfaces_12geometry
 use interfaces_13io_mpi
 use interfaces_13iovars
 use interfaces_13paw
 use interfaces_13psp
 use interfaces_13recipspace
 use interfaces_14iowfdenpot
 use interfaces_14occeig
 use interfaces_14wvl_wfs
 use interfaces_15common
 use interfaces_16response
 use interfaces_18seqpar
 use interfaces_21drive, except_this_one => gstate
 use interfaces_lib01hidempi
#else
 use defs_berry
#endif
!End of the abilint section

 implicit none

!Arguments ------------------------------------
!scalars
 integer,intent(inout) :: iexit,nspinor
 real(dp),intent(in) :: cpui,walli
 character(len=6),intent(in) :: codvsn
 type(MPI_type),intent(inout) :: mpi_enreg
 type(datafiles_type),intent(inout) :: dtfil
 type(dataset_type),intent(inout) :: dtset
 type(pawang_type),intent(inout) :: pawang
 type(pseudopotential_type),intent(inout) :: psps
 type(results_gs_type),intent(out) :: results_gs
!arrays
 integer,intent(out) :: npwtot(dtset%nkpt)
 real(dp),intent(inout) :: acell(3),occ(dtset%mband*dtset%nkpt*dtset%nsppol)
 real(dp),intent(inout) :: rprim(3,3),vel(3,dtset%natom),xred(3,dtset%natom)
 type(pawrad_type),intent(inout) :: pawrad(psps%ntypat*psps%usepaw)
 type(pawtab_type),intent(inout) :: pawtab(psps%ntypat*psps%usepaw)

!Local variables-------------------------------
!Define file format for different type of files. Presently,
!only one file format is supported for each type of files, but this might
!change soon ...
!2   for wavefunction file, new format (version 2.0 and after)    (fform)
!52  for density rho(r)       (fformr)
!102 for potential V(r) file. (fformv)
!scalars
 integer,parameter :: fform=2,fformv=102,formeig=0,level=3,response=0
 integer,save :: nsym_old=-1
 integer :: accessfil,ask_accurate,bantot,blktyp,choice,fformr=52,fullinit
 integer :: gscase,iapp,iatom,icount,idir,ierr,ifft,ii,ilmn,index,initialized
 integer :: ionmov,ios,ir,iscf,ispden,isppol,itime,itimexit,itypat,ixfh,ixx
 integer :: lpawu,master,me,mgfftf,mpert,mpsang,msize,mu,mxfh,mygroup,nblok
 integer :: ncerr,ncid_hdr,nfftf,nfftftot,nfftot,normchoice,nproc,ntime,ntypat
 integer :: nxfh,nxfhr,openexit,option,optorth,prtvol,psp_gencond,pwind_alloc
 integer :: rdwr,rdwrpaw,restartxf,spaceworld,tim_mkrho,tmkmem,vrsddb
 real(dp) :: cpus,ecore,ecut_eff,ecutdg_eff,epulay,etot,fermie,gsqcut_eff
 real(dp) :: gsqcutc_eff,residm,rhosum,tolwfr,ucvol
 logical :: ex,od
 character(len=3) :: ipara
 character(len=4) :: tag
 character(len=500) :: message
 character(len=fnlen) :: ddbnm,dscrpt
 type(bandstructure_type) :: bstruct
 type(dens_sym_operator_type) :: densymop_gs
 type(efield_type) :: dtefield
 type(hdr_type) :: hdr
 type(pawfgr_type) :: pawfgr
 type(wffile_type) :: wff1,wffnew,wffnow
 type(wvl_data) :: wvl
!arrays
 integer :: ngfft(18),ngfftf(18)
 integer,allocatable :: atindx(:),atindx1(:),blkflg(:),indsym(:,:,:)
 integer,allocatable :: irrzon(:,:,:),kg(:,:),nattyp(:),npwarr(:),symrec(:,:,:)
 integer,pointer :: pwind(:,:,:)
 real(dp) :: blknrm(3),blkqpt(9),corstr(6),gmet(3,3),gprimd(3,3)
 real(dp) :: k0(3),rmet(3,3),rprimd(3,3),tsec(2)
 real(dp),allocatable :: amass(:),blkval(:,:),cg(:,:),doccde(:),dyfrx2(:,:,:)
 real(dp),allocatable :: eigen(:),ph1d(:,:),ph1df(:,:),phnons(:,:,:),resid(:)
 real(dp),allocatable :: rhog(:,:),rhor(:,:),rhowfg(:,:),rhowfr(:,:),start(:,:)
 real(dp),allocatable :: work(:),xfhist(:,:,:,:),xred_old(:,:)
 real(dp),allocatable :: ylm(:,:),ylmgr(:,:,:)
 real(dp),pointer :: pwnsfac(:,:)
 character(len=fnlen) :: tmpfil(7)
 type(pawrhoij_type),allocatable :: pawrhoij(:)
!no_abirules

! ***********************************************************************
!DEBUG
!write(6,*)' gstate : enter'
!write(6,*) 'nfft',dtset%nfft
!stop
!ENDDEBUG

 call timab(32,1,tsec)
 call timab(33,1,tsec)

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

  if (mpi_enreg%me == 0 .and. dtset%outputXML == 1) then
    ! gstate() will handle a dataset, so we output the dataSet markup.
    write(ab_xml_out, "(A)") '  <dataSet>'
    ! We output the variables of the dataset given in argument.
!    call outvarsXML()
  end if

!Set up mpi informations from the dataset
 if (mpi_enreg%parareel == 0) then
  mpi_enreg%paral_level=2
  call initmpi_gs(dtset,mpi_enreg)
 else
  mpi_enreg%paral_level=1
 end if

 call initmpi_fft(dtset,mpi_enreg)

!Define FFT grid(s) sizes (be careful !)
!See NOTES in the comments at the beginning of this file.
 ngfft(:)=dtset%ngfft(:)
 if (psps%usepaw==1) then
  if (dtset%pawecutdg >= 1.0000001_dp*dtset%ecut) then
   pawfgr%usefinegrid=1
   nfftf=dtset%nfftdg;mgfftf=dtset%mgfftdg;ngfftf(:)=dtset%ngfftdg(:)
   nfftot=ngfft(1)*ngfft(2)*ngfft(3)
   nfftftot=ngfftf(1)*ngfftf(2)*ngfftf(3)
   allocate(pawfgr%coatofin(nfftot),pawfgr%fintocoa(nfftftot))
   call indgrid(pawfgr%coatofin,pawfgr%fintocoa,nfftot,nfftftot,ngfft,ngfftf)
  else !this is a simple transfer, this can be done in parallel with only local info
   pawfgr%usefinegrid=0
   nfftf=dtset%nfft;mgfftf=dtset%mgfft;ngfftf(:)=dtset%ngfft(:)
   allocate(pawfgr%coatofin(dtset%nfft),pawfgr%fintocoa(dtset%nfft))
   do ii=1,dtset%nfft;pawfgr%coatofin(ii)=ii;pawfgr%fintocoa(ii)=ii;end do
  end if
  pawfgr%natom=dtset%natom
  pawfgr%nfftc=dtset%nfft;pawfgr%mgfftc=dtset%mgfft;pawfgr%ngfftc(:)=dtset%ngfft(:)
  pawfgr%nfft =nfftf     ;pawfgr%mgfft=mgfftf      ;pawfgr%ngfft(:)=ngfftf(:)
  ecutdg_eff = dtset%pawecutdg * (dtset%dilatmx)**2
  ecut_eff   = dtset%ecut      * (dtset%dilatmx)**2
 else
  pawfgr%usefinegrid=0
  nfftf=dtset%nfft;mgfftf=dtset%mgfft;ngfftf(:)=dtset%ngfft(:)
  allocate(pawfgr%coatofin(0),pawfgr%fintocoa(0))
  ecut_eff= dtset%ecut * (dtset%dilatmx)**2
  ecutdg_eff=ecut_eff
 end if

!
!  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

!Init spaceworld
 call xcomm_world(mpi_enreg,spaceworld)
 master =0
!Define me

!BEGIN TF_CHANGES

call xme_init(mpi_enreg,me)
!END TF_CHANGES

!Define nproc

 call xproc_max(nproc,ierr)

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

 ntime=dtset%ntime

!Option input variables
 ionmov   =dtset%ionmov
 iscf     =dtset%iscf
 restartxf=dtset%restartxf

!Create names for the temporary files based on dtfil%filnam_ds(5)
!by appending adequate string.
!'_WF1' -> dtfil%unwft1
!'_WF2' -> dtfil%unwft2
!'_KG' ->  dtfil%unkg
!'_DUM' -> tmp_unit (real dummy name)
!'_YLM' -> dtfil%unylm
!'_PAW' -> dtfil%unpaw
 tmpfil(1)=trim(dtfil%filnam_ds(5))//'_WF1'
 tmpfil(2)=trim(dtfil%filnam_ds(5))//'_WF2'
 tmpfil(3)=trim(dtfil%filnam_ds(5))//'_KG'
 tmpfil(4)=trim(dtfil%filnam_ds(5))//'_DUM'
 tmpfil(6)=trim(dtfil%filnam_ds(5))//'_YLM'
 tmpfil(7)=trim(dtfil%filnam_ds(5))//'_PAW'

 if(mpi_enreg%paral_compil_kpt==1)then
! This is the parallel case : the index of the processor must be appended
  call int2char4(mpi_enreg%me,tag)
  ixx=1
  if (mpi_enreg%paral_compil_mpio == 1 .and. dtset%accesswff == 1 ) ixx=3
  do ii=ixx,7
   tmpfil(ii)=trim(tmpfil(ii))//'_P-'//tag
  end do
 end if

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

 initialized=0
 ecore=zero

 results_gs%grewtn(:,:)=zero
 call energies_init(results_gs%energies)
 results_gs%pel(1:3)   =zero

!Set up for iterations
 allocate(amass(dtset%natom))
 call setup1(acell,amass,dtset%amu,bantot,&
& ecutdg_eff,ecut_eff,gmet,gprimd,gsqcut_eff,gsqcutc_eff,dtset%iboxcut,dtset%intxc,ionmov,&
& dtset%natom,dtset%nband,ngfftf,ngfft,dtset%nkpt,dtset%nqpt,dtset%nsppol,dtset%nsym,psps%ntypat,&
& dtset%qptn,response,rmet,rprim,rprimd,dtset%typat,ucvol,psps%usepaw)
 call status(0,dtfil%filstat,iexit,level,'call kpgio    ')
!Set up the basis sphere of planewaves
 allocate(kg(3,dtset%mpw*dtset%mkmem),npwarr(dtset%nkpt))
 call kpgio(ecut_eff,dtset%exchn2n3,gmet,dtset%istwfk,kg,tmpfil(3),dtset%kptns,dtset%mkmem,&
& dtset%nband,dtset%nkpt,'PERS',mpi_enreg,dtset%mpw,npwarr,npwtot,dtset%nsppol,dtfil%unkg)

!Set up the Ylm for each k point
 allocate(ylm(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm))
 allocate(ylmgr(dtset%mpw*dtset%mkmem,3,psps%mpsang*psps%mpsang*psps%useylm))
 if (psps%useylm==1) then
  if(dtset%mkmem==0) open(dtfil%unylm,file=tmpfil(6),form='unformatted',status='unknown')
  call status(0,dtfil%filstat,iexit,level,'call initylmg ')
  option=0;if (dtset%prtstm==0.and.iscf>0) option=1
  call initylmg(gprimd,kg,dtset%kptns,dtset%mkmem,mpi_enreg,psps%mpsang,dtset%mpw,dtset%nband,dtset%nkpt,&
&               npwarr,dtset%nsppol,option,rprimd,dtfil%unkg,dtfil%unylm,ylm,ylmgr)
 end if

 call timab(33,2,tsec)


!Open and read pseudopotential files
 call status(0,dtfil%filstat,iexit,level,'call pspini   ')
 call pspini(dtset,ecore,psp_gencond,gsqcutc_eff,gsqcut_eff,level,&
&            pawrad,pawtab,psps,rprimd)

 call timab(33,1,tsec)

 ! In case of isolated computations, ecore must set to zero
 ! because its contribution is counted in the ewald energy
 ! as the ion-ion interaction.
 if (dtset%icoultrtmt == 1) then
  ecore = 0._dp
 end if

  ! WVL - Now that psp data are available, we compute rprimd, acell...
  !       from the atomic positions.
  if (dtset%usewvl == 1) then
    call wvl_setBoxGeometry(acell, dtset, mpi_enreg, psps, rprimd, xred)
    rprim(:, :)    = reshape((/ &
         & real(1., dp), real(0., dp), real(0., dp), 
          real(0., dp), real(1., dp), real(0., dp), &
          real(0., dp), real(0., dp), real(1., dp) /), (/ 3, 3 /))
    nfftf          = dtset%nfft
    mgfftf         = dtset%mgfft
    ngfftf(:)      = dtset%ngfft(:)
  end if

!Initialize band structure datatype
 allocate(doccde(bantot),eigen(bantot))
 doccde(:)=zero ; eigen(:)=zero
 call bstruct_init(bantot,bstruct,doccde,eigen,dtset%istwfk,dtset%kptns,&
& dtset%nband,dtset%nkpt,npwarr,dtset%nsppol,occ,dtset%wtk)
 deallocate(doccde,eigen)

!DEBUG
!write(6,*)' gstate : return for test memory leak '
!return
!ENDDEBUG

!Initialize PAW atomic occupancies
 if (psps%usepaw==1) then
  allocate(pawrhoij(dtset%natom))
  call initrhoij(psps%indlmn,psps%lmnmax,dtset%lpawu,dtset%natom,dtset%nspden,&
&                dtset%nsppol,dtset%ntypat,pawrhoij,pawtab,dtset%spinat,dtset%typat)
 end if

!Initialize header
 gscase=0
 call hdr_init(bstruct,codvsn,dtset,hdr,pawtab,gscase,psps)

!Update header, with evolving variables, when available
!Here, rprimd, xred and occ are available
 etot=hdr%etot ; fermie=hdr%fermie ; residm=hdr%residm
 call hdr_update(bantot,etot,fermie,hdr,dtset%natom,&
&                residm,rprimd,occ,pawrhoij,psps%usepaw,xred)

!Clean band structure datatype (should use it more in the future !)
 call bstruct_clean(bstruct)

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

 allocate(cg(2,dtset%mpw*nspinor*dtset%mband*dtset%mkmem*dtset%nsppol))
 allocate(eigen(dtset%mband*dtset%nkpt*dtset%nsppol))
 allocate(resid(dtset%mband*dtset%nkpt*dtset%nsppol))
 eigen(:)=0.0_dp ; resid(:)=0.0_dp
! mpi_enreg%paralbd=0 ; ask_accurate=0
 ask_accurate=0
  ! WVL - Branching, allocating wavefunctions as wavelets.
  if (dtset%usewvl == 1) then
    ! Create access arrays for wavefunctions and allocate them.
    call wvl_init_type_wfs(dtset, mpi_enreg, psps, rprimd, wvl%wfs, xred)
    ! Create access arrays for projectors and allocate them.
    ! Compute projectors from each atom.
    call wvl_init_type_proj(dtset, mpi_enreg, wvl%projectors, psps, rprimd, xred)
  end if

!XG 020711 : dtfil should not be reinitialized here !!!
 if (mpi_enreg%parareel == 1) then
  if (mpi_enreg%ipara > 0 ) then
   if (mpi_enreg%jpara == 0) then
    dtfil%ireadwf = 0
   else
    dtfil%ireadwf = 0
     if (mpi_enreg%paral_compil_mpio == 1 .and. dtset%accesswff == 1 ) then
      dtfil%fnamewffk=trim(dtfil%filnam_ds(4))//'_WFK'
     else
      if (mpi_enreg%ipara < 11) write(ipara,'(i1)')mpi_enreg%ipara-1
      if (mpi_enreg%ipara >= 11) write(ipara,'(i2)')mpi_enreg%ipara-1
      if (mpi_enreg%ipara >= 101) write(ipara,'(i3)')mpi_enreg%ipara-1
      dtfil%fnamewffk=trim(dtfil%filnam_ds(4))//'_WFK_'//ipara
     end if
   end if
  else
   dtfil%ireadwf = 0
  end if
 end if

!Initialize wavefunctions.
!Warning : ideally, results_gs%fermie and results_gs%residm
!should not be initialized here. One might make them separate variables.

 wff1%unwff=dtfil%unwff1
 optorth=1   !if (psps%usepaw==1) optorth=0
 if(psps%usepaw==1 .and. dtfil%ireadwf==1)optorth=0
 call inwffil(ask_accurate,cg,dtset,dtset%ecut,ecut_eff,eigen,dtset%exchn2n3,&
& formeig,gmet,hdr,dtfil%ireadwf,dtset%istwfk,kg,dtset%kptns,&
& dtset%localrdwf,dtset%mband,dtset%mkmem,mpi_enreg,&
& dtset%mpw,dtset%nband,ngfft,dtset%nkpt,npwarr,dtset%nspden,nspinor,dtset%nsppol,dtset%nsym,occ,&
& optorth,psps,prtvol,rprimd,dtset%symafm,dtset%symrel,dtset%tnons,&
& dtfil%unkg,wff1,wffnow,dtfil%unwff1,dtfil%unwft1,&
& dtfil%fnamewffk,wvl%wfs,tmpfil(1))

 if (psps%usepaw==1.and.dtfil%ireadwf==1)then
  do iatom=1,dtset%natom
   pawrhoij(iatom)%nspden=hdr%pawrhoij(iatom)%nspden
   pawrhoij(iatom)%lmn2_size=hdr%pawrhoij(iatom)%lmn2_size
   do ilmn=1,pawrhoij(iatom)%lmn2_size
    do ispden=1,pawrhoij(iatom)%nspden
     pawrhoij(iatom)%rhoijp(ilmn,ispden)=hdr%pawrhoij(iatom)%rhoijp(ilmn,ispden)+zero
     pawrhoij(iatom)%rhoijselect(ilmn,ispden)=hdr%pawrhoij(iatom)%rhoijselect(ilmn,ispden)+0
    end do
   end do
   do ispden=1,pawrhoij(iatom)%nspden
    pawrhoij(iatom)%nrhoijsel(ispden)=hdr%pawrhoij(iatom)%nrhoijsel(ispden)+0
   end do
  end do
! Has to update header again (because pawrhoij has changed)
  call hdr_update(bantot,etot,fermie,hdr,dtset%natom,&
&                 residm,rprimd,occ,pawrhoij,psps%usepaw,xred)
 end if

!DEBUG
!write(6,*)' gstate : stop for test memory leak '
!call hdr_clean(hdr)
!return
!ENDDEBUG

!Initialize xf history (should be put in inwffil)
 nxfh=0
 if(restartxf>=1 .and. dtfil%ireadwf==1)then

! Should exchange the data about history in parallel localrdwf==0
  if(mpi_enreg%paral_compil_kpt==1 .and. dtset%localrdwf==0)then
   write(message, '(a,a,a,a,a,a)' )ch10,&
&   ' gstate : BUG -',ch10,&
&   '  It is not yet possible to use non-zero restartxf,',ch10,&
&   '  in parallel, when localrdwf=0. Sorry for this ...'
   call wrtout(6,message,'COLL')
   call leave_new('COLL')
  end if

  allocate(xfhist(3,dtset%natom+4,2,0))
  call outxfhist(nxfh,dtset%natom,mxfh,xfhist,2,wff1,ios)
  deallocate(xfhist)

  if(ios>0)then
   write(message, '(a,a,a,a,a,a)' )ch10,&
&   ' gstate : BUG -',ch10,&
&   '  An error occurred reading the input wavefunction file,',ch10,&
&   '  with restartxf=1.'
   call wrtout(6,message,'COLL')
   call leave_new('COLL')
  else if(ios==0)then
   write(message, '(a,a,i4,a)' )ch10,&
&   ' gstate : reading',nxfh,' (x,f) history pairs from input wf file.'
   call wrtout(6,message,'COLL')
   call wrtout(ab_out,message,'COLL')
  end if
!WARNING : should check that restartxf is not negative
!WARNING : should check that restartxf /= only when dtfil%ireadwf is activated
 end if

!Allocate the xf history array : takes into account the existing
!pairs, minus those that will be discarded, then those that will
!be computed, governed by ntime, and some additional pairs
!(needed when it will be possible to use xfhist for move.f)
 mxfh=(nxfh-restartxf+1)+ntime+5
 if(mpi_enreg%parareel==1)mxfh=mxfh+500  ! XG020711 : why this value ?
 allocate(xfhist(3,dtset%natom+4,2,mxfh))
!WARNING : should check that the number of atoms in the wf file and natom
!are the same

!Initialize the xf history array
 if(nxfh>=restartxf .and. nxfh>0)then
! Eventually skip some of the previous history
  if(restartxf>=2)then
   do ixfh=1,restartxf-1
    call WffReadSkipRec(ios,1,wff1)
   end do
  end if

! Read and store the relevant history
  nxfhr=nxfh-restartxf+1
  call outxfhist(nxfhr,dtset%natom,mxfh,xfhist,3,wff1,ios)
 end if

!Close wff1, if it was ever opened (in inwffil)
 if (dtfil%ireadwf==1) then
  call WffClose(wff1,ierr)
 end if

!Initialize second wavefunction file if needed
 if(dtset%mkmem==0 .and. dtset%nstep/=0) then
  write(message, '(a,i4,a,a)' )&
&  ' gstate about to open unit',dtfil%unwft2,' for file=',trim(tmpfil(2))
  call wrtout(06,message,'PERS')

#if defined HAVE_NETCDF
     if(dtset%accesswff==2) then
    !  Create empty netCDF file
        ncerr = nf90_create(path=trim(tmpfil(2)), cmode=NF90_CLOBBER, ncid=ncid_hdr)
        call handle_ncerr(ncerr," create netcdf wavefunction file")
        ncerr = nf90_close(ncid_hdr)
        call handle_ncerr(ncerr," close netcdf wavefunction file")
     else if(dtset%accesswff==3) then
        write (std_out,*) "FIXME: ETSF I/O support in gstate"
     end if
#endif


  call WffOpen(dtset%accesswff,spaceworld,tmpfil(2),ierr,wffnew,master,me,dtfil%unwft2)
 end if

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

!Further setup
 allocate(start(3,dtset%natom))
 call setup2(dtset%dedlnn,dtset%ecut,results_gs%energies%e_pulay,iscf, &
           & dtset%istwfk,dtset%natom,dtset%nkpt,npwtot,start,ucvol,dtset%wtk,xred)

!Allocation for forces and atomic positions
 allocate(xred_old(3,dtset%natom))

!Do symmetry stuff only for nsym>1
 nfftot=ngfft(1)*ngfft(2)*ngfft(3)
 allocate(irrzon(nfftot**(1-1/dtset%nsym),2,dtset%nspden/dtset%nsppol))
 allocate(phnons(2,nfftot**(1-1/dtset%nsym),dtset%nspden/dtset%nsppol))
 irrzon(:,:,:)=0
 allocate(indsym(4,dtset%nsym,dtset%natom),symrec(3,3,dtset%nsym))

 if (dtset%nsym>1) then

  call status(0,dtfil%filstat,iexit,level,'call setsym   ')
  call setsym(densymop_gs,indsym,irrzon,iscf,dtset%natom,&
&  nfftot,ngfft,dtset%nspden,dtset%nsppol,dtset%nsym,&
&  phnons,dtset%symafm,symrec,dtset%symrel,dtset%tnons,dtset%typat,xred)

! Make sure dtset%iatfix does not break symmetry
  call status(0,dtfil%filstat,iexit,level,'call fixsym   ')
  call fixsym(dtset%iatfix,indsym,dtset%natom,dtset%nsym)

 else

! The symrec array is used by initberry even in case nsym = 1
  symrec(:,:,1) = 0
  symrec(1,1,1) = 1 ; symrec(2,2,1) = 1 ; symrec(3,3,1) = 1

 end if

!Electric field: initialization
 if ((dtset%berryopt < 0).or.(dtset%berryopt == 4)) then
   nullify(pwind,pwnsfac)
   call initberry(dtefield,dtfil,dtset,gmet,kg,dtset%mband,dtset%mkmem,mpi_enreg,&
&              dtset%mpw,dtset%nkpt,npwarr,dtset%nsppol,dtset%nsym,occ,pwind,pwind_alloc,pwnsfac,rprimd,symrec)
 else
   pwind_alloc = 1
   allocate(pwind(pwind_alloc,2,3),pwnsfac(2,pwind_alloc))
 end if

!Timing for initialisation period
 call timab(33,2,tsec)
 call timab(34,1,tsec)

!Compute new occupation numbers, in case wavefunctions and eigenenergies
!were read from disk, occupation scheme is metallic (this excludes iscf=-1),
!and occupation numbers are required by iscf
 if( dtfil%ireadwf==1 .and. &
&    (dtset%occopt>=3.and.dtset%occopt<=7) .and. &
&    (iscf>0 .or. iscf==-3) ) then

  call status(0,dtfil%filstat,iexit,level,'call newocc   ')
  allocate(doccde(dtset%mband*dtset%nkpt*dtset%nsppol))
! Warning : ideally, results_gs%entropy should not be set up here XG 20011007
! Warning : ideally, results_gs%fermie should not be set up here XG 20011007
! Do not take into account the possible STM bias
  call newocc(doccde,eigen,results_gs%energies%entropy,&
&  results_gs%energies%e_fermie,&
&  dtset%fixmom,dtset%mband,dtset%nband,&
&  dtset%nelect,dtset%nkpt,nspinor,dtset%nsppol,occ,&
&  dtset%occopt,prtvol,zero,dtset%tphysel,dtset%tsmear,dtset%wtk)
  deallocate(doccde)

 else
! Warning : ideally, results_gs%entropy should not be set up here XG 20011007
  results_gs%energies%entropy=zero
 end if

!Generate an index table of atoms, in order for them to be used
!type after type.
 ntypat=psps%ntypat
 allocate(atindx(dtset%natom),atindx1(dtset%natom),nattyp(ntypat))
 index=1
 do itypat=1,ntypat
  nattyp(itypat)=0
  do iatom=1,dtset%natom
   if(dtset%typat(iatom)==itypat)then
    atindx(iatom)=index
    atindx1(index)=iatom
    index=index+1
    nattyp(itypat)=nattyp(itypat)+1
   end if
  end do
 end do

!Compute structure factor phases for current atomic pos:
 if (dtfil%ireadwf==0.or.psps%usepaw==1) then
  call status(0,dtfil%filstat,iexit,level,'call getph    ')
  allocate(ph1d(2,3*(2*dtset%mgfft+1)*dtset%natom))
  call getph(atindx,dtset%natom,dtset%ngfft(1),dtset%ngfft(2),dtset%ngfft(3),&
&            ph1d,xred)
 end if

!PAW: 1- Initialize values for several arrays unchanged during iterations
!     2- Initialize data for LDA+U
!     3- Eventually open temporary storage file
 if(psps%usepaw==1) then
! 1-
  if (psp_gencond==1) then
   call timab(553,1,tsec)
   call pawinit(psps%indlmn,dtset%pawlcutd,dtset%pawlmix,psps%lmnmax,psps%mpsang,psps%n1xccc,&
&       dtset%pawnphi,dtset%nsym,dtset%pawntheta,psps%ntypat,pawang,pawrad,pawtab,dtset%pawxcdev)
   call timab(553,2,tsec)
  end if
  if (psp_gencond==1.or.nsym_old/=dtset%nsym) then
   call setsymrhoij(gprimd,pawang%l_max-1,dtset%nsym,dtset%pawprtvol,&
&                   rprimd,dtset%symafm,symrec,pawang%zarot)
   nsym_old=dtset%nsym
  end if
  do itypat=1,ntypat
   pawtab(itypat)%usepawu=0
  end do
! 2-
!      Initialize and compute data for LDA+U
  if (dtset%usepawu>0) then
   if(dtset%nspden/=2) then
    write(message, '(a,a,i3,a,a,a,a)' )  ch10,&
&    ' gstate : nspden =',dtset%nspden,ch10,&
&    ' nspden.ne.2 not permitted in LDA+U',ch10,&
&    ' action: change your input file'
    call wrtout(06,message,'COLL')
    call leave_new('COLL')
   endif
   call pawpuinit(dtset%jpawu,dtset%lpawu,psps%indlmn,psps%lmnmax,ntypat,pawang,&
&    pawrad,pawtab,dtset%upawu,dtset%usepawu)
  end if
! 3-
  if(dtset%mkmem==0) then
   open(dtfil%unpaw,file=tmpfil(7),form='unformatted',status='unknown')
   rewind(unit=dtfil%unpaw)
  end if
 end if

!Get starting charge density : rhor as well as rhog
 allocate(rhog(2,nfftf),rhor(nfftf,dtset%nspden))
 if (iscf>0) then
  if(dtfil%ireadden/=0)then

   rdwr=1;rdwrpaw=psps%usepaw;if(dtfil%ireadwf/=0) rdwrpaw=0
   call ioarr(accessfil,rhor,results_gs%etotal,fformr,dtfil%fildensin,hdr,&
&              nfftf,dtset%nspden,rdwr,rdwrpaw,ngfft)
   if (rdwrpaw/=0) pawrhoij(1:dtset%natom)=hdr%pawrhoij(1:dtset%natom)
!  Compute up+down rho(G) by fft
   allocate(work(nfftf));work(:)=rhor(:,1)
   call fourdp(1,rhog,work,-1,mpi_enreg,nfftf,ngfftf,0)
   deallocate(work)

  else if(dtfil%ireadwf/=0)then

!  Obtain the charge density from wfs that were read previously
!  Be careful: in PAW, rho does not include the compensation
!              density (to be added in scfcv.F90) !
   call status(0,dtfil%filstat,iexit,level,'call mkrho    ')
!   tim_mkrho=1 ; mpi_enreg%paralbd=0
   tim_mkrho=1
   if (psps%usepaw==1) then
    allocate(rhowfg(2,dtset%nfft),rhowfr(dtset%nfft,dtset%nspden))
    call mkrho(cg,densymop_gs,dtset,irrzon,kg,&
&    mpi_enreg,npwarr,nspinor,occ,phnons,rhowfg,rhowfr,tim_mkrho,ucvol,&
&    dtfil%unkg,wffnow,wvl%wfs)
    call transgrid(mpi_enreg,dtset%nspden,+1,1,1,pawfgr,rhowfg,rhog,rhowfr,rhor)
    deallocate(rhowfg,rhowfr)
   else
    call mkrho(cg,densymop_gs,dtset,irrzon,kg,&
&    mpi_enreg,npwarr,nspinor,occ,phnons,rhog,rhor,tim_mkrho,ucvol,&
&    dtfil%unkg,wffnow,wvl%wfs)
   end if

  else if(dtfil%ireadwf==0)then

!  Crude, but realistic initialisation of the density
!  There is not point to compute it from random wavefunctions
!  except with wavelets.
   allocate(ph1df(2,3*(2*mgfftf+1)*dtset%natom))
   if (psps%usepaw==1) then
    call status(0,dtfil%filstat,iexit,level,'call getph    ')
    call getph(atindx,dtset%natom,ngfftf(1),ngfftf(2),ngfftf(3),&
&              ph1df,xred)
   else
    ph1df(:,:)=ph1d(:,:)
   end if
   call status(0,dtfil%filstat,iexit,level,'call initro   ')
#if defined MPI_FFT
   call MPI_COMM_SIZE(mpi_enreg%comm_fft,ir, itypat); write(6,*) ir,itypat
#endif
   if (dtset%usewvl == 0) then
     call initro(atindx,dtset%densty,gmet,gsqcut_eff,psps%usepaw,mgfftf,mpi_enreg,dtset%natom,nattyp,&
  &   nfftf,ngfftf,dtset%nspden,ntypat,ph1df,rhog,rhor,&
  &   dtset%spinat,ucvol,dtset%ziontypat,dtset%znucl)
   else
     call wvl_mkrho(dtset, mpi_enreg, occ, rhor, wvl%wfs)
   end if

  end if

 else if (iscf==-1.or.iscf==-2.or.iscf==-3) then

  call status(0,dtfil%filstat,iexit,level,'call ioarr    ')
! Read rho(r) from a disk file
  rdwr=1;rdwrpaw=psps%usepaw
! Note : results_gs%etotal is read here,
! and might serve in the tddft routine, but it is contrary to the
! intended use of results_gs ...
! Warning : should check the use of results_gs%fermie
! Warning : should check the use of results_gs%residm
! One might make them separate variables.

!DEBUG
!  write(6,*)' gstate : before ioarr, reading the density '
!ENDDEBUG

  call ioarr(accessfil,rhor,results_gs%etotal,fformr,dtfil%fildensin,hdr,&
&  nfftf,dtset%nspden,rdwr,rdwrpaw,ngfft)
! Compute up+down rho(G) by fft
  call status(0,dtfil%filstat,iexit,level,'call fourdp   ')
  allocate(work(nfftf))
  work(:)=rhor(:,1)
  call fourdp(1,rhog,work,-1,mpi_enreg,nfftf,ngfftf,0)
  deallocate(work)

 else

! Disallowed value for iscf
  write(message, '(a,a,a,a,i12,a)' )  ch10,&
&   ' gstate : BUG -',ch10,&
&   '  iscf has disallowed value=',iscf,'.'
  call wrtout(06,message,'COLL')
  call leave_new('COLL')

 end if

 if (dtfil%ireadwf==0.or.psps%usepaw==1) deallocate(ph1d)
 if (dtfil%ireadwf==0.and.dtfil%ireadden==0.and.iscf>0) deallocate(ph1df)

!Debugging : print the different parts of rhor
! MPIWF Warning : this should not be parallelized over space, leave this debugging feature as such.
 if(prtvol==-level)then
  write(message,'(a)') '   ir     rhor(ir)     '
  call wrtout(06,message,'COLL')
  do ir=1,nfftf
   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(dtset%nsppol==2)then
     write(message,'(a,es13.6)')'      ',rhor(ir,2)
     call wrtout(06,message,'COLL')
    end if
   end if
  end do
 end if

 call status(0,dtfil%filstat,iexit,level,'end gstate(1) ')

 if(prtvol==-level)then
  write(message,'(a1,a,a1,a,i1,a)') ch10,&
&   ' gstate : before scfcv, move or brdmin ',&
&   ch10,'  prtvol=-',level,', debugging mode => stop '
  call wrtout(06,message,'COLL')
  call leave_new('COLL')
 end if

 call timab(34,2,tsec)
!Check whether exiting was required by the user.
!If found then do not start minimization steps
!At this first call to chkexi, initialize cpus, if it
!is non-zero (which would mean that no action has to be taken)
!Should do this in driver ...
 cpus=dtset%cpus
 if(abs(cpus)>1.0d-5)cpus=cpus+cpui
 openexit=1 ; if(dtset%chkexit==0) openexit=0
 call chkexi(cpus,dtfil%filnam_ds(1),iexit,ab_out,mpi_enreg,openexit)
!If immediate exit, and wavefunctions were not read, must zero eigenvalues
 if (iexit/=0) then
  eigen(:)=zero
 end if
 if (iexit==0) then

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

! If atoms are not being moved, use scfcv directly; else
! call move or brdmin which in turn calls scfcv.

  call timab(35,1,tsec)

  write(message,'(a,80a)')ch10,('=',mu=1,80)
  call wrtout(ab_out,message,'COLL')
  call wrtout(06,message,'COLL')
  if (ionmov==0) then

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

!Should merge this call with the call for ionmov==4 and 5
   iapp=0
!   mpi_enreg%paralbd=0
   call scfcv(atindx,atindx1,cg,cpus,densymop_gs,dtefield,dtfil,dtset,&
&   ecore,eigen,hdr,iapp,indsym,initialized,&
&   irrzon,kg,mpi_enreg,nattyp,nfftf,npwarr,nspinor,occ,&
&   pawang,pawfgr,pawrad,pawrhoij,pawtab,&
&   phnons,psps,pwind,pwind_alloc,pwnsfac,resid,results_gs,rhog,rhor,rprimd,&
&   symrec,wffnew,wffnow,wvl,xred,xred_old,ylm,ylmgr)

  else if (ionmov==1) then
!  Conduct molecular dynamics, with or without viscous damping

   call status(0,dtfil%filstat,iexit,level,'call move     ')
!   mpi_enreg%paralbd=0
    call move(amass,atindx,atindx1,cg,cpus,densymop_gs,dtefield,dtfil,dtset,&
&   ecore,eigen,hdr,indsym,initialized,irrzon,&
&   kg,mpi_enreg,&
&   nattyp,nfftf,npwarr,nspinor,occ,&
&   pawang,pawfgr,pawrad,pawrhoij,pawtab,&
&   phnons,psps,pwind,pwind_alloc,pwnsfac,resid,results_gs,rhog,rhor,rprimd,&
&   symrec,wffnew,wffnow,vel,wvl,xred,xred_old,ylm,ylmgr)

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

!  Apply Broyden method for structural optimization, as
!  implemented by Jean-Christophe Charlier (May 1992)

   call status(0,dtfil%filstat,iexit,level,'call brdmin   ')
!   mpi_enreg%paralbd=0

   call brdmin(acell,atindx,atindx1,cg,cpus,densymop_gs,dtefield,dtfil,dtset,&
&   ecore,eigen,hdr,indsym,initialized,irrzon,&
&   kg,mpi_enreg,mxfh,&
&   nattyp,nfftf,npwarr,nspinor,nxfh,occ,&
&   pawang,pawfgr,pawrad,pawrhoij,pawtab,&
&   phnons,psps,pwind,pwind_alloc,pwnsfac,resid,results_gs,rhog,rhor,rprim,&
&   symrec,wffnew,wffnow,vel,wvl,xfhist,xred,xred_old,ylm,ylmgr)
!  call mkrdim(acell,rprim,rprimd)

  else if (ionmov==4 .or. ionmov==5) then

   do itime=1,ntime

    call status(itime,dtfil%filstat,iexit,level,'call scfcv(mv)')

    if(ionmov==4)then
     if(mod(itime,2)==1)then
      write(message, '(a,a,i3,a)' ) ch10,' STEP NUMBER ',itime,&
&      ' : OPTIMIZE ELECTRONS ------------------------------------'
     else
      write(message, '(a,a,i3,a)' ) ch10,' STEP NUMBER ',itime,&
&      ' : OPTIMIZE ELECTRONS AND IONS ---------------------------'
     end if
    else
     write(message, '(a,a,i3,a)' ) ch10,' STEP NUMBER ',itime,&
&     ' : SIMPLE RELAXATION -------------------------------------'
    end if
    call wrtout(ab_out,message,'COLL')
    call wrtout(06,  message,'COLL')

!   In this case, iapp is simply itime
    iapp=itime
!    mpi_enreg%paralbd=0
    call scfcv(atindx,atindx1,cg,cpus,densymop_gs,dtefield,dtfil,dtset,ecore,&
&    eigen,hdr,iapp,indsym,initialized,irrzon,kg,mpi_enreg,&
&    nattyp,nfftf,npwarr,nspinor,occ,pawang,pawfgr,pawrad,pawrhoij,pawtab,&
&    phnons,psps,pwind,pwind_alloc,pwnsfac,resid,results_gs,rhog,rhor,rprimd,&
&    symrec,wffnew,wffnow,wvl,xred,xred_old,ylm,ylmgr)

    if(mod(itime,2)==1)then
!    When the SCF cycle dealt with electrons only,
!    check whether forces are below tolerance; if so, exit
!    from the itime loop
     itimexit=0 ; if(itime==ntime)itimexit=1
     call fconv(results_gs%fcart,dtset%iatfix,itimexit,itime,dtset%natom,&
&     ntime,0,1.0_dp,dtset%strtarget,results_gs%strten,dtset%tolmxf)
    end if
    if (itimexit/=0) exit

!   Check whether exiting was required by the user.
!   If found then beat a hasty exit from time steps
    if(dtset%chkexit==0) then
     openexit=0
    else
     openexit=1
    end if
    call chkexi(cpus,dtfil%filnam_ds(1),iexit,ab_out,mpi_enreg,openexit)
    if (iexit/=0) then
     iexit=0   ! In order not to exit of dataset loop automatically
     exit
    end if

   end do

  else if ( (ionmov>=6 .and. ionmov<=9) .or. ionmov==12) then

!  Molecular dynamics, using Verlet algorithm (ionmov=6)
!  or fake molecular dynamics for minimisation (ionmov=7)
!  or true molecular dynamics with Nose thermostat (ionmov=8)
!  or Langevin dynamics (ionmov=9) or Fei Zhang algorithm (ionmov=12)

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

   call moldyn(acell,amass,atindx,atindx1,cg,cpus,densymop_gs,dtefield,dtfil,&
&   dtset,ecore,eigen,hdr,indsym,initialized,&
&   irrzon,kg,mpi_enreg,mxfh,&
&   nattyp,nfftf,npwarr,nspinor,nxfh,occ,&
&   pawang,pawfgr,pawrad,pawrhoij,pawtab,&
&   phnons,psps,pwind,pwind_alloc,pwnsfac,resid,results_gs,rhog,rhor,rprim,&
&   symrec,wffnew,wffnow,vel,wvl,xfhist,xred,xred_old,ylm,ylmgr)
!  call mkrdim(acell,rprim,rprimd)

  else if (ionmov == 10) then

   call delocint(acell,atindx,atindx1,cg,cpus,densymop_gs,dtefield,dtfil,&
&   dtset,ecore,eigen,hdr,indsym,initialized,irrzon,&
&   kg,mpi_enreg,mxfh,&
&   nattyp,nfftf,npwarr,nspinor,nxfh,occ,&
&   pawang,pawfgr,pawrad,pawrhoij,pawtab,&
&   phnons,psps,pwind,pwind_alloc,pwnsfac,resid,results_gs,rhog,rhor,rprim,&
&   symrec,wffnew,wffnow,vel,wvl,xfhist,xred,xred_old,ylm,ylmgr)
!  call mkrdim(acell,rprim,rprimd)

  else if (ionmov == 20) then
    ! Ground state call.
    iapp = 0
    ! Ionic positions relaxation using DIIS. This algorithm is fast
    ! and converge to the nearest singular point (where gradient vanishes).
    ! This is a good algorithm to precisely tune saddle-points.
    call diisRelax(atindx, atindx1, cg, cpus, densymop_gs, dtefield, &
                 & dtfil, dtset, ecore, eigen, hdr, iapp, indsym, initialized, &
                 & irrzon, kg, mpi_enreg, nattyp, nfftf, npwarr, nspinor, occ, pawang, &
                 & pawfgr, pawrad, pawrhoij, pawtab, phnons, psps, pwind, pwind_alloc, pwnsfac, &
                 & resid, results_gs, rhog, rhor, rprimd, symrec, &
                 & wffnew, wffnow, wvl, xred, xred_old, ylm, ylmgr)

  else
!  Not an allowed option
   write(message, '(a,a,a,a,i12,a,a)' ) ch10,&
&   ' gstate : BUG -',ch10,&
&   '  Disallowed value for ionmov=',ionmov,ch10,&
&   '  Allowed values are 0 to 5.'
   call wrtout(06,message,'COLL')
   call leave_new('COLL')
  end if

  call timab(35,2,tsec)

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

!End of the check of hasty exit
 end if

 call timab(36,1,tsec)

 write(message, '(80a,a,a,a,a)' ) ('=',mu=1,80),ch10,ch10,&
&  ' ----iterations are completed or convergence reached----',&
&  ch10
 call wrtout(ab_out,message,'COLL')
 call wrtout(06,  message,'COLL')

!Close the unneeded temporary data files, if any.
!Other files are closed in clnup1.
 if (dtset%mkmem==0) then
  close (unit=dtfil%unkg,status='delete')
  if (psps%useylm==1) close (unit=dtfil%unylm,status='delete')
  if (psps%usepaw==1) close (unit=dtfil%unpaw,status='delete')
  call WffDelete(wffnew,ierr)
 end if

!Update the header, before using it
 call hdr_update(bantot,results_gs%etotal,results_gs%energies%e_fermie,hdr,dtset%natom,&
&                results_gs%residm,rprimd,occ,pawrhoij,psps%usepaw,xred)

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

  call outwf(cg,dtfil,dtset,eigen,dtfil%filnam_ds(4),hdr,kg,dtset%kptns,&
  & dtset%mband,dtset%mkmem,mpi_enreg,dtset%mpw,mxfh,dtset%natom,dtset%nband,dtset%nfft,ngfft,&
  & dtset%nkpt,npwarr,dtset%nqpt,nspinor,dtset%nsppol,dtset%nstep,nxfh,&
  & occ,resid,response,wffnow,wvl%wfs,xfhist)

 if(dtset%prtwf==2)then
  call outqmc(cg,dtset,eigen,gprimd,hdr,kg,dtset%mband,dtset%mkmem,mpi_enreg,dtset%mpw,dtset%nkpt,npwarr,&
&  nspinor,dtset%nsppol,occ,psps,results_gs)
 end if

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

 call clnup1(acell,dtset%dosdeltae,eigen,dtset%enunit,&
& results_gs%energies%e_fermie,dtfil%filnam_ds(4),&
& results_gs%fred,dtset%iatfix,iscf,dtset%kptns,dtset%kptopt,dtset%mband,dtset%mkmem,mpi_enreg,dtset%mpw,&
& dtset%natom,dtset%nband,nfftf,ngfftf,dtset%nkpt,dtset%nspden,nspinor,dtset%nsppol,dtset%nstep,&
& occ,dtset%occopt,dtset%prtdos,dtset%prteig,dtset%optforces,dtset%prtstm,prtvol,&
& resid,rhor,rprimd,dtset%tphysel,dtset%tsmear,results_gs%vxcavg,dtset%wtk,xred)

 if (iscf>0 .and. dtset%prtstm==0) then
   call status(0,dtfil%filstat,iexit,level,'call prtene   ')
   call prtene(dtset,results_gs%energies,ab_out,psps%usepaw)
 end if

! Open the formatted derivative database file, and write the
! preliminary information
! In the // case, only one processor writes the energy and
! the gradients to the DDB

 if ((psps%usepaw == 0).and.(mpi_enreg%me==0).and.((iscf > 0).or.&
&    (dtset%berryopt == -1).or.(dtset%berryopt) == -3)) then

  call status(0,dtfil%filstat,iexit,level,'call ioddb8   ')
  vrsddb=010929
  dscrpt=' Note : temporary (transfer) database '
  choice=2
  ddbnm=trim(dtfil%filnam_ds(4))//'_DDB'
! tolwfr must be initialized here, but it is a dummy value
  tolwfr=1.0_dp
  call ioddb8 (choice,dscrpt,ddbnm,dtset%natom,dtset%mband,&
&   dtset%nkpt,dtset%nsym,psps%ntypat,dtfil%unddb,vrsddb,&
&   acell,dtset%amu,dtset%dilatmx,dtset%ecut,dtset%ecutsm,&
&   dtset%intxc,dtset%iscf,dtset%ixc,dtset%kpt,dtset%kptnrm,&
&   dtset%natom,dtset%nband,ngfft,dtset%nkpt,dtset%nspden,nspinor,&
&   dtset%nsppol,dtset%nsym,psps%ntypat,occ,dtset%occopt,&
&   rprim,dtset%sciss,dtset%spinat,dtset%symafm,dtset%symrel,&
&   dtset%tnons,tolwfr,dtset%tphysel,dtset%tsmear,&
&   dtset%typat,dtset%wtk,xred,psps%ziontypat,dtset%znucl)

  if (iscf > 0) then
    nblok = 2          ! 1st blok = energy, 2nd blok = gradients
  else
    nblok = 1
  end if
  fullinit = 0
  call psddb8 (choice,psps%dimekb,psps%ekb,fullinit,psps%indlmn,&
&   psps%lmnmax,psps%lnmax,nblok,&
&   psps%ntypat,dtfil%unddb,psps%pspso,psps%usepaw,psps%useylm,vrsddb)

  mpert = dtset%natom + 6   ; msize = 3*mpert
  allocate(blkflg(msize),blkval(2,msize))

  blkflg(:) = 0       ; blkval(:,:) = zero
  blkqpt(:) = zero    ; blknrm(:) = one

! Write total energy to the DDB
  if (iscf > 0) then
    blktyp = 0
    blkval(1,1) = results_gs%etotal
    blkflg(1) = 1
    call blok8(blkflg,blknrm,blkqpt,blktyp,blkval,choice,mpert,&
&      msize,dtset%natom,dtfil%unddb)
  end if

! Write gradients to the DDB
  blktyp = 4
  blkflg(:) = 0       ; blkval(:,:) = zero
  index = 0
  if (iscf > 0) then
    do iatom = 1, dtset%natom
      do idir = 1, 3
        index = index + 1
        blkflg(index) = 1
        blkval(1,index) = results_gs%fred(idir,iatom)
      end do
    end do
  end if

  index = 3*dtset%natom + 3
  if ((abs(dtset%berryopt) == 1).or.(abs(dtset%berryopt) == 3)) then
    do idir = 1, 3
      index = index + 1
      if (dtset%rfdir(idir) == 1) then
        blkflg(index) = 1
        blkval(1,index) = results_gs%pel(idir)
      end if
    end do
  end if

  index = 3*dtset%natom + 6
  if (iscf > 0) then
    blkflg(index+1:index+6) = 1
    blkval(1,index+1:index+6) = results_gs%strten(1:6)
  end if

  call blok8(blkflg,blknrm,blkqpt,blktyp,blkval,choice,mpert,&
&    msize,dtset%natom,dtfil%unddb)

  deallocate(blkflg,blkval)

! Close DDB
  close(dtfil%unddb)

 end if

 if (dtset%nstep>0 .and. dtset%prtstm==0 .and. dtset%positron==0) then
  call status(0,dtfil%filstat,iexit,level,'call clnup2   ')
  call clnup2(psps%n1xccc,results_gs%fred,results_gs%gresid,&
&  results_gs%grewtn,&
&  results_gs%grxc,iscf,dtset%natom,dtset%optforces,dtset%optstress,prtvol,start,&
&  results_gs%strten,results_gs%synlgr,psps%usepaw,xred)
 end if

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

!Deallocate arrays
 deallocate(amass,atindx,atindx1,cg,eigen,indsym)
 deallocate(irrzon,kg,npwarr,nattyp,phnons,resid)
 deallocate(rhog,rhor,start,symrec,xfhist,xred_old,ylm,ylmgr)
 deallocate(pawfgr%fintocoa,pawfgr%coatofin)
 if (psps%usepaw==1) then
  do iatom=1,dtset%natom
   deallocate(pawrhoij(iatom)%nrhoijsel,pawrhoij(iatom)%rhoijp,&
&             pawrhoij(iatom)%rhoijselect)
  end do
  deallocate(pawrhoij)
 end if

 if (dtset%usepawu>0) then
  do itypat=1,ntypat
   if(dtset%lpawu(itypat)/=-1) then
!    deallocate(pawtab(itypat)%vee)
    deallocate(pawtab(itypat)%lnproju)
    deallocate(pawtab(itypat)%phiphjint)
   end if
  end do
 end if
  ! WVL - Branching, deallocating wavefunctions as wavelets.
  if (dtset%usewvl == 1) then
    call wvl_free_type_wfs(wvl%wfs)
    call wvl_free_type_proj(wvl%projectors)
  end if

 if ((dtset%berryopt < 0).or.(dtset%berryopt == 4)) then
   deallocate(pwind,pwnsfac)
   deallocate(dtefield%ikpt_dk,dtefield%idxkstr)
   deallocate(dtefield%sflag,dtefield%cgindex,dtefield%kgindex)
   deallocate(dtefield%fkptns,dtefield%indkk_f2ibz,dtefield%i2fbz)
   if (mpi_enreg%paral_compil_kpt == 1) then
    deallocate(mpi_enreg%kptdstrb)
    if (dtset%berryopt == 4) then
     deallocate(mpi_enreg%kptdstrbi,dtefield%cgqindex,dtefield%nneigh)
    end if
   end if
 else
  deallocate(pwind,pwnsfac)
 end if
 if (dtset%berryopt == 4) deallocate(dtefield%smat)

!Clean the header
 call hdr_clean(hdr)

!DEBUG
!write(6,*)' gstate : return for test memory leak '
!return
!ENDDEBUG

  if (mpi_enreg%me == 0 .and. dtset%outputXML == 1) then
    ! The dataset given in argument has been treated, then
    ! we output its variables.
!    call outvarsXML()
    ! gstate() will handle a dataset, so we output the dataSet markup.
    write(ab_xml_out, "(A)") '  </dataSet>'
  end if

!Clean the MPI informations
 if (mpi_enreg%parareel == 0) then
   call clnmpi_gs(dtset,mpi_enreg)
 end if

 call clnmpi_fft(dtset,mpi_enreg)

 write(message, '(a,a)' ) ch10,' gstate : exiting '
 call wrtout(06,message,'COLL')

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

 call timab(36,2,tsec)
 call timab(32,2,tsec)

end subroutine gstate
!!***

Generated by  Doxygen 1.6.0   Back to index