TABLE OF CONTENTS


ABINIT/m_gstate [ Modules ]

[ Top ] [ Modules ]

NAME

  m_gstate

FUNCTION

COPYRIGHT

  Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, JYR, MKV, MT, FJ, MB, DJA)
  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

 20 #if defined HAVE_CONFIG_H
 21 #include "config.h"
 22 #endif
 23 
 24 #include "abi_common.h"
 25 
 26 module m_gstate
 27 
 28  use defs_basis
 29  use defs_datatypes
 30  use defs_abitypes
 31  use defs_rectypes
 32  use m_errors
 33  use m_xmpi
 34  use m_abicore
 35  use libxc_functionals
 36  use m_exit
 37  use m_crystal
 38  use m_crystal_io
 39  use m_scf_history
 40  use m_abimover
 41  use m_wffile
 42  use m_rec
 43  use m_efield
 44  use m_ddb
 45  use m_bandfft_kpt
 46  use m_invovl
 47  use m_gemm_nonlop
 48  use m_tetrahedron
 49  use m_wfk
 50  use m_nctk
 51  use m_hdr
 52  use m_ebands
 53 
 54  use m_time,             only : timab
 55  use m_symtk,            only : matr3inv
 56  use m_io_tools,         only : open_file
 57  use m_occ,              only : newocc, getnel
 58  use m_ddb_hdr,          only : ddb_hdr_type, ddb_hdr_init, ddb_hdr_free, ddb_hdr_open_write
 59  use m_fstrings,         only : strcat, sjoin
 60  use m_geometry,         only : fixsym, mkradim, metric
 61  use m_kpts,             only : tetra_from_kptrlatt
 62  use m_kg,               only : kpgio, getph
 63  use m_fft,              only : fourdp
 64  use m_pawang,           only : pawang_type
 65  use m_pawrad,           only : pawrad_type
 66  use m_pawtab,           only : pawtab_type
 67  use m_pawcprj,          only : pawcprj_type,pawcprj_free,pawcprj_alloc, pawcprj_getdim
 68  use m_pawfgr,           only : pawfgr_type, pawfgr_init, pawfgr_destroy
 69  use m_abi2big,          only : wvl_occ_abi2big, wvl_setngfft, wvl_setBoxGeometry
 70  use m_energies,         only : energies_type, energies_init
 71  use m_args_gs,          only : args_gs_type
 72  use m_results_gs,       only : results_gs_type
 73  use m_pawrhoij,         only : pawrhoij_type, pawrhoij_copy, pawrhoij_free
 74  use m_paw_dmft,         only : init_sc_dmft,destroy_sc_dmft,print_sc_dmft,paw_dmft_type,readocc_dmft
 75  use m_paw_sphharm,      only : setsym_ylm
 76  use m_paw_init,         only : pawinit,paw_gencond
 77  use m_paw_occupancies,  only : initrhoij
 78  use m_paw_correlations, only : pawpuxinit
 79  use m_orbmag,           only : initorbmag,destroy_orbmag,orbmag_type
 80  use m_paw_uj,           only : pawuj_ini,pawuj_free,pawuj_det
 81  use m_data4entropyDMFT, only : data4entropyDMFT_t, data4entropyDMFT_init, data4entropyDMFT_destroy
 82  use m_electronpositron, only : electronpositron_type,init_electronpositron,destroy_electronpositron, &
 83                                 electronpositron_calctype
 84  use m_scfcv,            only : scfcv_t, scfcv_init, scfcv_destroy, scfcv_run
 85  use m_dtfil,            only : dtfil_init_time, status
 86  use m_jellium,          only : jellium
 87  use m_iowf,             only : outwf
 88  use m_outqmc,           only : outqmc
 89  use m_ioarr,            only : ioarr,read_rhor
 90  use m_inwffil,          only : inwffil
 91  use m_spacepar,         only : setsym
 92  use m_mkrho,            only : mkrho, initro, prtrhomxmn
 93  use m_initylmg,         only : initylmg
 94  use m_pspini,           only : pspini
 95  use m_mover,            only : mover
 96  use m_mpinfo,           only : proc_distrb_cycle
 97  use m_common,           only : setup1, prteigrs, prtene
 98  use m_fourier_interpol, only : transgrid
 99  use m_psolver,          only : psolver_kernel
100  use m_paw2wvl,          only : paw2wvl, wvl_paw_free
101  use m_berryphase_new,   only : init_e_field_vars,prtefield
102  use m_wvl_wfs,          only : wvl_wfs_set, wvl_wfs_free, wvl_wfs_lr_copy
103  use m_wvl_rho,          only : wvl_initro, wvl_mkrho
104  use m_wvl_descr_psp,    only : wvl_descr_psp_set, wvl_descr_free, wvl_descr_atoms_set, wvl_descr_atoms_set_sym
105  use m_wvl_denspot,      only : wvl_denspot_set, wvl_denspot_free
106  use m_wvl_projectors,   only : wvl_projectors_set, wvl_projectors_free
107 
108 #if defined HAVE_GPU_CUDA
109  use m_alloc_hamilt_gpu, only : alloc_hamilt_gpu, dealloc_hamilt_gpu
110 #endif
111 
112  use defs_wvltypes,      only : wvl_data,coulomb_operator,wvl_wf_type
113 #if defined HAVE_BIGDFT
114  use BigDFT_API,         only : wvl_timing => timing,xc_init,xc_end,XC_MIXED,XC_ABINIT,&
115 &                               local_potential_dimensions,nullify_gaussian_basis, &
116 &                               copy_coulomb_operator,deallocate_coulomb_operator
117 #else
118  use defs_wvltypes,      only : coulomb_operator
119 #endif
120 #if defined HAVE_LOTF
121  use defs_param_lotf,    only : lotfparam_init
122 #endif
123 
124  implicit none
125 
126  private

ABINIT/outxfhist [ Functions ]

[ Top ] [ Functions ]

NAME

 outxfhist

FUNCTION

  read/write xfhist

COPYRIGHT

 Copyright (C) 2003-2018 ABINIT group (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

  option =
   1: write
   2: read only nxfh
   3: read xfhist
  response =
   0: GS wavefunctions
   1: RF wavefunctions
  natom = number of atoms in unit cell
  mxfh = last dimension of the xfhist array

OUTPUT

  ios = error code returned by read operations

SIDE EFFECTS

  nxfh = actual number of (x,f) history pairs, see xfhist array
  wff2 = structured info for wavefunctions
  xfhist(3,natom+4,2,ab_xfh%mxfh) = (x,f) history array, also including
   rprim and stress

PARENTS

      gstate

CHILDREN

      xderiveread,xderiverrecend,xderiverrecinit,xderivewrecend
      xderivewrecinit,xderivewrite

SOURCE

2729 subroutine outxfhist(ab_xfh,natom,option,wff2,ios)
2730 
2731  use defs_basis
2732  use m_abicore
2733  use m_abimover
2734  use m_xmpi
2735  use m_wffile
2736  use m_errors
2737 #if defined HAVE_NETCDF
2738  use netcdf
2739 #endif
2740 
2741 !This section has been created automatically by the script Abilint (TD).
2742 !Do not modify the following lines by hand.
2743 #undef ABI_FUNC
2744 #define ABI_FUNC 'outxfhist'
2745 !End of the abilint section
2746 
2747  implicit none
2748 
2749 !Arguments ------------------------------------
2750  integer          ,intent(in)    :: natom,option
2751  integer          ,intent(out)   :: ios
2752  type(wffile_type),intent(inout)    :: wff2
2753  type(ab_xfh_type),intent(inout) :: ab_xfh
2754 
2755 !Local variables-------------------------------
2756  integer :: ierr,ixfh,ncid_hdr,spaceComm,xfdim2
2757  real(dp),allocatable :: xfhist_tmp(:)
2758  character(len=500) :: message
2759 !no_abirules
2760 #if defined HAVE_NETCDF
2761  integer :: ncerr
2762  integer :: nxfh_id, mxfh_id, xfdim2_id, dim2inout_id, dimr3_id,xfhist_id
2763  integer :: nxfh_tmp,mxfh_tmp,xfdim2_tmp,dim2inout_tmp
2764 #endif
2765 
2766 
2767 ! *************************************************************************
2768 
2769 !DEBUG
2770 !write(std_out,*)'outxfhist  : enter, option = ', option
2771 !ENDDEBUG
2772  ncid_hdr = wff2%unwff
2773  xfdim2 = natom+4
2774 
2775  ios = 0
2776 
2777 !### (Option=1) Write out content of all iterations
2778 !#####################################################################
2779  if ( option == 1 ) then
2780 
2781 !  Write the (x,f) history
2782    if (wff2%iomode == IO_MODE_FORTRAN) then
2783      write(unit=wff2%unwff)ab_xfh%nxfh
2784      do ixfh=1,ab_xfh%nxfh
2785        write(unit=wff2%unwff)ab_xfh%xfhist(:,:,:,ixfh)
2786      end do
2787 
2788    else if (wff2%iomode == IO_MODE_FORTRAN_MASTER) then
2789 !    FIXME: should copy the xfhist to other processors, and check that we are on the master to read in this case
2790 !    if node is master
2791      write(message, "(A,A,A,A)") ch10, " outxfhist: ERROR -", ch10, &
2792 &     'iomode == -1 (localrdwf ) has not been coded yet for xfhist rereading.'
2793      MSG_ERROR(message)
2794 
2795      write(unit=wff2%unwff)ab_xfh%nxfh
2796      do ixfh=1,ab_xfh%nxfh
2797        write(unit=wff2%unwff)ab_xfh%xfhist(:,:,:,ixfh)
2798      end do
2799 
2800 !    insert mpi broadcast here
2801 
2802    else if(wff2%iomode==IO_MODE_MPI)then
2803      ABI_ALLOCATE(xfhist_tmp,(3*(natom+4)*2))
2804      spaceComm=xmpi_comm_self
2805      call xderiveWRecInit(wff2,ierr)
2806      call xderiveWrite(wff2,ab_xfh%nxfh,ierr)
2807      call xderiveWRecEnd(wff2,ierr)
2808      do ixfh=1,ab_xfh%nxfh
2809        xfhist_tmp(:)=reshape(ab_xfh%xfhist(:,:,:,ixfh),(/3*(natom+4)*2/))
2810        call xderiveWRecInit(wff2,ierr)
2811        call xderiveWrite(wff2,xfhist_tmp,3*(natom+4)*2,spaceComm,ierr)
2812        call xderiveWRecEnd(wff2,ierr)
2813      end do
2814      ABI_DEALLOCATE(xfhist_tmp)
2815 
2816 #if defined HAVE_NETCDF
2817    else if (wff2%iomode == IO_MODE_NETCDF) then
2818 !    check if nxfh and xfhist are defined
2819      ncerr = nf90_inq_dimid(ncid=ncid_hdr,name="nxfh",dimid=nxfh_id)
2820 
2821      if (ncerr /= NF90_NOERR) then
2822 !      need to define everything
2823        ncerr = nf90_redef (ncid=ncid_hdr)
2824        NCF_CHECK_MSG(ncerr," outxfhist : going to define mode ")
2825 
2826        ncerr = nf90_def_dim(ncid=ncid_hdr,name="dim2inout",len=2,dimid=dim2inout_id)
2827        NCF_CHECK_MSG(ncerr," outxfhist : define dim2inout")
2828        ncerr = nf90_def_dim(ncid=ncid_hdr,name="mxfh",len=ab_xfh%mxfh,dimid=mxfh_id)
2829        NCF_CHECK_MSG(ncerr," outxfhist : define mxfh")
2830        ncerr = nf90_def_dim(ncid=ncid_hdr,name="nxfh",len=ab_xfh%nxfh,dimid=nxfh_id)
2831        NCF_CHECK_MSG(ncerr," outxfhist : define nxfh")
2832        ncerr = nf90_def_dim(ncid=ncid_hdr,name="xfdim2",len=xfdim2,dimid=xfdim2_id)
2833        NCF_CHECK_MSG(ncerr," outxfhist : define xfdim2")
2834 
2835        ncerr = nf90_inq_dimid(ncid=ncid_hdr,name="dimr3",dimid=dimr3_id)
2836        NCF_CHECK_MSG(ncerr," outxfhist : inquire dimr3")
2837 
2838 !      ab_xfh%xfhist(3,natom+4,2,ab_xfh%mxfh)
2839        ncerr = nf90_def_var(ncid=ncid_hdr,name="xfhist",xtype=NF90_DOUBLE,&
2840 &       dimids=(/dimr3_id,xfdim2_id,dim2inout_id,mxfh_id/),varid=xfhist_id)
2841        NCF_CHECK_MSG(ncerr," outxfhist : define xfhist")
2842 
2843 !      End define mode and go to data mode
2844        ncerr = nf90_enddef(ncid=ncid_hdr)
2845        NCF_CHECK_MSG(ncerr," outxfhist : enddef call ")
2846      else
2847 !      check that the dimensions are correct
2848        ncerr = nf90_inq_dimid(ncid=ncid_hdr,name="nxfh",dimid=nxfh_id)
2849        NCF_CHECK_MSG(ncerr," outxfhist : inquire nxfh")
2850        ncerr = nf90_Inquire_Dimension(ncid=ncid_hdr,dimid=nxfh_id,&
2851 &       len=nxfh_tmp)
2852        NCF_CHECK_MSG(ncerr,"  outxfhist : get nxfh")
2853        ncerr = nf90_inq_dimid(ncid=ncid_hdr,name="xfdim2",dimid=xfdim2_id)
2854        NCF_CHECK_MSG(ncerr," outxfhist : inquire xfdim2")
2855        ncerr = nf90_Inquire_Dimension(ncid=ncid_hdr,dimid=xfdim2_id,&
2856 &       len=xfdim2_tmp)
2857        NCF_CHECK_MSG(ncerr,"  outxfhist : get xfdim2")
2858        ncerr = nf90_inq_dimid(ncid=ncid_hdr,name="mxfh",dimid=mxfh_id)
2859        NCF_CHECK_MSG(ncerr," outxfhist : inquire mxfh")
2860        ncerr = nf90_Inquire_Dimension(ncid=ncid_hdr,dimid=mxfh_id,&
2861 &       len=mxfh_tmp)
2862        NCF_CHECK_MSG(ncerr,"  outxfhist : get mxfh")
2863        ncerr = nf90_inq_dimid(ncid=ncid_hdr,name="dim2inout",dimid=dim2inout_id)
2864        NCF_CHECK_MSG(ncerr," outxfhist : inquire dim2inout")
2865        ncerr = nf90_Inquire_Dimension(ncid=ncid_hdr,dimid=dim2inout_id,&
2866 &       len=dim2inout_tmp)
2867        NCF_CHECK_MSG(ncerr,"  outxfhist : get dim2inout")
2868 
2869        ncerr = nf90_inq_varid(ncid=ncid_hdr,name="xfhist",varid=xfhist_id)
2870        NCF_CHECK_MSG(ncerr," outxfhist : inquire xfhist")
2871 
2872        if (mxfh_tmp /= ab_xfh%mxfh .or. dim2inout_tmp /= 2 .or. xfdim2_tmp /= xfdim2) then
2873          write (message,"(A)") 'outxfhist : ERROR xfhist has bad dimensions in NetCDF file. Can not re-write it.'
2874          MSG_ERROR(message)
2875        end if
2876 
2877      end if
2878 
2879 !    Now fill the data
2880      ncerr = nf90_put_var(ncid=ncid_hdr,varid=xfhist_id,values=ab_xfh%xfhist)
2881      NCF_CHECK_MSG(ncerr," outxfhist : fill xfhist")
2882 
2883 !    end NETCDF definition ifdef
2884 #endif
2885    end if  ! end iomode if
2886 
2887 !  ### (Option=2) Read in number of iterations
2888 !  #####################################################################
2889  else if ( option == 2 ) then
2890 
2891    if (wff2%iomode == IO_MODE_FORTRAN) then
2892      read(unit=wff2%unwff,iostat=ios)ab_xfh%nxfh
2893 
2894    else if (wff2%iomode == IO_MODE_FORTRAN_MASTER) then
2895 !    FIXME: should copy the xfhist to other processors, and check that we are on the master to read in this case
2896 !    if node is master
2897      write(message, "(A,A,A,A)") ch10, " outxfhist: ERROR -", ch10, &
2898 &     'iomode == -1 (localrdwf ) has not been coded yet for xfhist rereading.'
2899      MSG_ERROR(message)
2900 
2901      read(unit=wff2%unwff,iostat=ios)ab_xfh%nxfh
2902 
2903    else if (wff2%iomode == IO_MODE_MPI) then
2904      call xderiveRRecInit(wff2,ierr)
2905      call xderiveRead(wff2,ab_xfh%nxfh,ierr)
2906      call xderiveRRecEnd(wff2,ierr)
2907 
2908 #if defined HAVE_NETCDF
2909    else if (wff2%iomode == IO_MODE_NETCDF) then
2910      ncerr = nf90_inq_dimid(ncid=ncid_hdr,name="nxfh",dimid=nxfh_id)
2911      NCF_CHECK_MSG(ncerr," outxfhist : inquire nxfh")
2912      ncerr = nf90_Inquire_Dimension(ncid=ncid_hdr,dimid=nxfh_id,&
2913 &     len=ab_xfh%nxfh)
2914      NCF_CHECK_MSG(ncerr,"  outxfhist : get nxfh")
2915 #endif
2916    end if
2917 
2918 !  ### (Option=3) Read in iteration content
2919 !  #####################################################################
2920  else if ( option == 3 ) then
2921    if (wff2%iomode == IO_MODE_FORTRAN) then
2922      do ixfh=1,ab_xfh%nxfhr
2923        read(unit=wff2%unwff,iostat=ios)ab_xfh%xfhist(:,:,:,ixfh)
2924      end do
2925    else if (wff2%iomode == IO_MODE_FORTRAN_MASTER) then
2926 !    FIXME: should copy the xfhist to other processors, and check that we are on the master to read in this case
2927 !    if node is master
2928      write(message, "(A,A,A,A)") ch10, " outxfhist: ERROR -", ch10, &
2929 &     'iomode == -1 (localrdwf ) has not been coded yet for xfhist rereading.'
2930      MSG_ERROR(message)
2931 
2932      do ixfh=1,ab_xfh%nxfhr
2933        read(unit=wff2%unwff,iostat=ios)ab_xfh%xfhist(:,:,:,ixfh)
2934      end do
2935 
2936    else if (wff2%iomode == IO_MODE_MPI) then
2937      ABI_ALLOCATE(xfhist_tmp,(3*(natom+4)*2))
2938      spaceComm=xmpi_comm_self
2939      do ixfh=1,ab_xfh%nxfhr
2940        call xderiveRRecInit(wff2,ierr)
2941        call xderiveRead(wff2,xfhist_tmp,3*(natom+4)*2,spaceComm,ierr)
2942        call xderiveRRecEnd(wff2,ierr)
2943        xfhist_tmp(:)=xfhist_tmp(:)
2944      end do
2945      ABI_DEALLOCATE(xfhist_tmp)
2946    end if
2947 
2948 !  FIXME: should this be inside the if not mpi as above for options 1 and 2?
2949 !  it is placed here because the netcdf read is a single operation
2950 #if defined HAVE_NETCDF
2951    if (wff2%iomode == IO_MODE_NETCDF) then
2952      ncerr = nf90_inq_dimid(ncid=ncid_hdr,name="nxfh",dimid=nxfh_id)
2953      NCF_CHECK_MSG(ncerr," outxfhist : inquire nxfh")
2954      ncerr = nf90_Inquire_Dimension(ncid=ncid_hdr,dimid=nxfh_id,&
2955 &     len=ab_xfh%nxfhr)
2956      NCF_CHECK_MSG(ncerr,"  outxfhist : get nxfh")
2957 
2958      ncerr = nf90_inq_varid(ncid=ncid_hdr,varid=xfhist_id,name="xfhist")
2959      NCF_CHECK_MSG(ncerr," outxfhist : inquire xfhist")
2960      ncerr = nf90_get_var(ncid=ncid_hdr,varid=xfhist_id,values=ab_xfh%xfhist,&
2961 &     start=(/1,1,1,1/),count=(/3,natom+4,2,ab_xfh%nxfhr/))
2962      NCF_CHECK_MSG(ncerr," outxfhist : read xfhist")
2963    end if
2964 #endif
2965 
2966  else
2967 !  write(std_out,*)' outxfhist : option ', option , ' not available '
2968    write(message, "(A,A,A,A,I3,A)") ch10, "outxfhist: ERROR -", ch10, &
2969 &   "option ", option, " not available."
2970    MSG_ERROR(message)
2971  end if
2972 
2973 end subroutine outxfhist

m_gstate/clnup1 [ Functions ]

[ Top ] [ m_gstate ] [ Functions ]

NAME

 clnup1

FUNCTION

 Perform "cleanup" at end of execution of gstate routine.

INPUTS

  acell(3)=length scales of primitive translations (bohr)
  dosdeltae=DOS delta of Energy
  dtset <type(dataset_type)>=all input variables in this dataset
  eigen(mband*nkpt*nsppol)=eigenvalues (hartree) for all bands
                           at each k point
  enunit=choice for units of output eigenvalues: 0=>hartree,
   1=> eV, 2=> hartree and eV
  fermie=fermi energy (Hartree)
  fnameabo_dos=filename of output DOS file
  fnameabo_eig=filename of output EIG file
  fred(3,natom)=d(E)/d(xred) (hartree)
  iatfix(3,natom)=0 if not fixed along specified direction,
                  1 if fixed
  iscf=parameter controlling scf or non-scf choice
  kptopt=option for the generation of k points
  kptns(3,nkpt)=k points in terms of recip primitive translations
  mband=maximum number of bands
  mpi_enreg=information about MPI parallelization
  natom=number of atoms in unit cell
  nband(nkpt*nsppol)=number of bands
  nfft=(effective) number of FFT grid points (for this processor)
  ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
  nkpt=number of k points
  nspden=number of spin-density components
  nsppol=1 for unpolarized, 2 for spin-polarized
  nstep=desired number of electron iteration steps
  occ(maxval(nband(:))*nkpt*nsppol)=occupancies for each band and k point
  occopt=option for occupancies
  prtdos= if == 1, will print the density of states
  prtfor= if >0, will print the forces
  prtstm= input variable prtstm
  prtvol=control print volume and debugging
  resid(mband*nkpt*nsppol)=squared residuals for each band and k point where
                     resid(n,k)=|<C(n,k)|(H-e(n,k))|C(n,k)>|^2
  rhor(nfft,nspden)=electron density (electrons/bohr^3)
  rprimd(3,3)=dimensional real space primitive translations (bohr)
  tphysel="physical" electronic temperature with FD occupations
  tsmear=smearing energy or temperature (if metal)
  vxcavg=average of vxc potential
  wtk(nkpt)=real(dp) array of k-point weights
  xred(3,natom)=reduced atomic coordinates

OUTPUT

  (only print and write to disk)

PARENTS

      gstate

CHILDREN

      getnel,metric,prteigrs,prtrhomxmn,prtxf,write_eig,wrtout

SOURCE

1946 subroutine clnup1(acell,dtset,eigen,fermie,&
1947   & fnameabo_dos,fnameabo_eig,fred,&
1948   & mpi_enreg,nfft,ngfft,occ,prtfor,&
1949   & resid,rhor,rprimd,vxcavg,xred)
1950 
1951 
1952 !This section has been created automatically by the script Abilint (TD).
1953 !Do not modify the following lines by hand.
1954 #undef ABI_FUNC
1955 #define ABI_FUNC 'clnup1'
1956 !End of the abilint section
1957 
1958  implicit none
1959 
1960 !Arguments ------------------------------------
1961 !scalars
1962  integer,intent(in) :: nfft
1963  integer,intent(in) :: prtfor
1964  real(dp),intent(in) :: fermie
1965  real(dp),intent(in) :: vxcavg
1966  character(len=*),intent(in) :: fnameabo_dos,fnameabo_eig
1967  type(dataset_type),intent(in) :: dtset
1968  type(MPI_type),intent(in) :: mpi_enreg
1969 !arrays
1970  integer,intent(in)  :: ngfft(18)
1971  real(dp),intent(in) :: acell(3)
1972  real(dp),intent(in) :: eigen(dtset%mband*dtset%nkpt*dtset%nsppol)
1973  real(dp),intent(in) :: fred(3,dtset%natom)
1974  real(dp),intent(in) :: resid(dtset%mband*dtset%nkpt*dtset%nsppol)
1975  real(dp),intent(in) :: rhor(nfft,dtset%nspden)
1976  real(dp),intent(in) :: rprimd(3,3)
1977  real(dp),intent(in) :: xred(3,dtset%natom)
1978  real(dp),intent(inout) :: occ(dtset%mband*dtset%nkpt*dtset%nsppol)
1979 
1980 !Local variables-------------------------------
1981 !scalars
1982  integer,parameter :: master=0
1983  integer :: comm,iatom,ii,iscf_dum,iwfrc,me,nnonsc,option,unitdos
1984  real(dp) :: entropy,grmax,grsum,maxocc,nelect,tolwf,ucvol
1985  real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3)
1986  character(len=500) :: message
1987  character(len=fnlen) filename
1988 !arrays
1989  real(dp),allocatable :: doccde(:)
1990 
1991 ! ****************************************************************
1992 
1993  comm=mpi_enreg%comm_cell; me=xmpi_comm_rank(comm)
1994 
1995  if(dtset%prtstm==0)then ! Write reduced coordinates xred
1996    write(message, '(a,i5,a)' )' reduced coordinates (array xred) for',dtset%natom,' atoms'
1997    call wrtout(ab_out,message,'COLL')
1998    do iatom=1,dtset%natom
1999      write(message, '(1x,3f20.12)' ) xred(:,iatom)
2000      call wrtout(ab_out,message,'COLL')
2001    end do
2002  end if
2003 
2004 !Write reduced gradients if iscf > 0 and dtset%nstep>0 and
2005  if (dtset%iscf>=0.and.dtset%nstep>0.and.dtset%prtstm==0) then
2006 
2007 !  Compute absolute maximum and root mean square value of gradients
2008    grmax=0.0_dp
2009    grsum=0.0_dp
2010    do iatom=1,dtset%natom
2011      do ii=1,3
2012 !      To be activated in v5.5
2013 !      grmax=max(grmax,abs(fred(ii,iatom)))
2014        grmax=max(grmax,fred(ii,iatom))
2015        grsum=grsum+fred(ii,iatom)**2
2016      end do
2017    end do
2018    grsum=sqrt(grsum/dble(3*dtset%natom))
2019 
2020    write(message, '(1x,a,1p,e12.4,a,e12.4,a)' )'rms dE/dt=',grsum,'; max dE/dt=',grmax,'; dE/dt below (all hartree)'
2021    call wrtout(ab_out,message,'COLL')
2022    do iatom=1,dtset%natom
2023      write(message, '(i5,1x,3f20.12)' ) iatom,fred(1:3,iatom)
2024      call wrtout(ab_out,message,'COLL')
2025    end do
2026 
2027  end if
2028 
2029  if(dtset%prtstm==0)then
2030 
2031 !  Compute and write out dimensional cartesian coords and forces:
2032    call wrtout(ab_out,' ','COLL')
2033 
2034 !  (only write forces if iscf > 0 and dtset%nstep>0)
2035    if (dtset%iscf<0.or.dtset%nstep<=0.or.prtfor==0) then
2036      iwfrc=0
2037    else
2038      iwfrc=1
2039    end if
2040 
2041    call prtxf(fred,dtset%iatfix,ab_out,iwfrc,dtset%natom,rprimd,xred)
2042 
2043 !  Write length scales
2044    write(message, '(1x,a,3f16.12,a)' )'length scales=',acell,' bohr'
2045    call wrtout(ab_out,message,'COLL')
2046    write(message, '(14x,a,3f16.12,a)' )'=',Bohr_Ang*acell(1:3),' angstroms'
2047    call wrtout(ab_out,message,'COLL')
2048 
2049  end if
2050 
2051  option=1; nnonsc=0; tolwf=0.0_dp
2052 
2053  if(dtset%iscf<0 .and. dtset%iscf/=-3)option=3
2054  iscf_dum=dtset%iscf
2055  if(dtset%nstep==0)iscf_dum=-1
2056 
2057  if(dtset%tfkinfunc==0)then
2058    if (me == master) then
2059      call prteigrs(eigen,dtset%enunit,fermie,fnameabo_eig,ab_out,&
2060 &     iscf_dum,dtset%kptns,dtset%kptopt,dtset%mband,&
2061 &     dtset%nband,dtset%nkpt,nnonsc,dtset%nsppol,occ,&
2062 &     dtset%occopt,option,dtset%prteig,dtset%prtvol,resid,tolwf,&
2063 &     vxcavg,dtset%wtk)
2064      call prteigrs(eigen,dtset%enunit,fermie,fnameabo_eig,std_out,&
2065 &     iscf_dum,dtset%kptns,dtset%kptopt,dtset%mband,&
2066 &     dtset%nband,dtset%nkpt,nnonsc,dtset%nsppol,occ,&
2067 &     dtset%occopt,option,dtset%prteig,dtset%prtvol,resid,tolwf,&
2068 &     vxcavg,dtset%wtk)
2069    end if
2070 
2071 #if defined HAVE_NETCDF
2072    if (dtset%prteig==1 .and. me == master) then
2073      filename=trim(fnameabo_eig)//'.nc'
2074      call write_eig(eigen,filename,dtset%kptns,dtset%mband,dtset%nband,dtset%nkpt,dtset%nsppol)
2075    end if
2076 #endif
2077  end if
2078 
2079 !Compute and print location of maximal and minimal density
2080  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
2081  call prtrhomxmn(std_out,mpi_enreg,nfft,ngfft,dtset%nspden,2,rhor,ucvol=ucvol)
2082  if( dtset%prtvol>1)then
2083    call prtrhomxmn(ab_out,mpi_enreg,nfft,ngfft,dtset%nspden,2,rhor,ucvol=ucvol)
2084  end if
2085 
2086 !If needed, print DOS (unitdos is closed in getnel, occ is not changed if option == 2
2087  if (dtset%prtdos==1 .and. me == master) then
2088    if (open_file(fnameabo_dos,message, newunit=unitdos, status='unknown', action="write", form='formatted') /= 0) then
2089      MSG_ERROR(message)
2090    end if
2091    rewind(unitdos)
2092    maxocc=two/(dtset%nspinor*dtset%nsppol)  ! Will not work in the fixed moment case
2093    option=2
2094    ABI_ALLOCATE(doccde,(dtset%mband*dtset%nkpt*dtset%nsppol))
2095    call getnel(doccde,dtset%dosdeltae,eigen,entropy,fermie,&
2096 &   maxocc,dtset%mband,dtset%nband,nelect,dtset%nkpt,&
2097 &   dtset%nsppol,occ,dtset%occopt,option,dtset%tphysel,&
2098 &   dtset%tsmear,unitdos,dtset%wtk)
2099    ABI_DEALLOCATE(doccde)
2100  end if
2101 
2102 end subroutine clnup1

m_gstate/clnup2 [ Functions ]

[ Top ] [ m_gstate ] [ Functions ]

NAME

 clnup2

FUNCTION

 Perform more "cleanup" after completion of iterations.
 This subroutine prints out more breakdown of force
 information, shifts of atomic positions, and stresses.

INPUTS

  fred(3,natom)=d(E_total)/d(xred) derivatives (hartree)
  grchempottn(3,natom)=d(E_chempot)/d(xred) derivatives (hartree)
  grewtn(3,natom)=d(E_Ewald)/d(xred) derivatives (hartree)
  grvdw(3,ngrvdw)=gradients of energy due to Van der Waals DFT-D2 dispersion (hartree)
  grxc(3,natom)=d(Exc)/d(xred) derivatives (0 without core charges)
  iscf=parameter controlling scf or non-scf iterations
  natom=number of atoms in unit cell
  ngrvdw=size of grvdw(:,:); can be 0 or natom according to dtset%vdw_xc
  n1xccc=dimension of xccc1d ; 0 if no XC core correction is used
  prtfor= >0 if forces have to be printed (0 otherwise)
  prtstr= >0 if stresses have to be printed (0 otherwise)
  prtvol=control print volume and debugging output
  start(3,natom)=starting coordinates in terms of real space
   primitive translations
  strten(6)=components of the stress tensor (hartree/bohr^3)
  synlgr(3,natom)=d(E_nlpsp)/d(xred) derivatives (hartree)
  xred(3,natom)=final coordinates in terms of primitive translations

OUTPUT

  (only print)

PARENTS

      gstate

CHILDREN

      wrtout

SOURCE

2322 subroutine clnup2(n1xccc,fred,grchempottn,gresid,grewtn,grvdw,grxc,iscf,natom,ngrvdw,&
2323 &                 prtfor,prtstr,prtvol,start,strten,synlgr,xred)
2324 
2325 
2326 !This section has been created automatically by the script Abilint (TD).
2327 !Do not modify the following lines by hand.
2328 #undef ABI_FUNC
2329 #define ABI_FUNC 'clnup2'
2330 !End of the abilint section
2331 
2332  implicit none
2333 
2334 !Arguments ------------------------------------
2335 !scalars
2336  integer,intent(in) :: iscf,n1xccc,natom,ngrvdw,prtfor,prtstr,prtvol
2337 !arrays
2338  real(dp),intent(in) :: fred(3,natom),grchempottn(3,natom),gresid(3,natom)
2339  real(dp),intent(in) :: grewtn(3,natom),grvdw(3,ngrvdw)
2340  real(dp),intent(in) :: grxc(3,natom),start(3,natom),strten(6),synlgr(3,natom)
2341  real(dp),intent(in) :: xred(3,natom)
2342 
2343 !Local variables-------------------------------
2344  character(len=*), parameter :: format01020 ="(i5,1x,3f20.12)"
2345 !scalars
2346  integer :: iatom,mu
2347  real(dp) :: devsqr,grchempot2
2348  character(len=500) :: message
2349 
2350 ! *************************************************************************
2351 !
2352 !DEBUG
2353 !write(std_out,*)' clnup2 : enter '
2354 !ENDDEBUG
2355 
2356 !Only print additional info for scf calculations
2357  if (iscf>=0) then
2358 
2359    if((prtvol>=10).and.(prtfor>0))then
2360 
2361      write(message, '(a,10x,a)' ) ch10,&
2362 &     '===> extra information on forces <==='
2363      call wrtout(ab_out,message,'COLL')
2364 
2365      write(message, '(a)' ) ' ewald contribution to reduced grads'
2366      call wrtout(ab_out,message,'COLL')
2367      do iatom=1,natom
2368        write(message,format01020) iatom,(grewtn(mu,iatom),mu=1,3)
2369        call wrtout(ab_out,message,'COLL')
2370      end do
2371 
2372      grchempot2=sum(grchempottn(:,:)**2)
2373      if(grchempot2>tol16)then
2374        write(message, '(a)' ) ' chemical potential contribution to reduced grads'
2375        call wrtout(ab_out,message,'COLL')
2376        do iatom=1,natom
2377          write(message,format01020) iatom,(grchempottn(mu,iatom),mu=1,3)
2378          call wrtout(ab_out,message,'COLL')
2379        end do
2380      end if
2381 
2382      write(message, '(a)' ) ' nonlocal contribution to red. grads'
2383      call wrtout(ab_out,message,'COLL')
2384      do iatom=1,natom
2385        write(message,format01020) iatom,(synlgr(mu,iatom),mu=1,3)
2386        call wrtout(ab_out,message,'COLL')
2387      end do
2388 
2389      write(message, '(a)' ) ' local psp contribution to red. grads'
2390      call wrtout(ab_out,message,'COLL')
2391      if (n1xccc/=0) then
2392        do iatom=1,natom
2393          write(message,format01020) iatom,fred(:,iatom)-&
2394 &         (grewtn(:,iatom)+grchempottn(:,iatom)+synlgr(:,iatom)+grxc(:,iatom)+gresid(:,iatom))
2395          call wrtout(ab_out,message,'COLL')
2396        end do
2397      else
2398        do iatom=1,natom
2399          write(message,format01020) iatom,fred(:,iatom)-&
2400 &         (grewtn(:,iatom)+grchempottn(:,iatom)+synlgr(:,iatom)+gresid(:,iatom))
2401          call wrtout(ab_out,message,'COLL')
2402        end do
2403      end if
2404 
2405      if (n1xccc/=0) then
2406        write(message, '(a)' ) ' core charge xc contribution to reduced grads'
2407        call wrtout(ab_out,message,'COLL')
2408        do iatom=1,natom
2409          write(message,format01020) iatom,(grxc(mu,iatom),mu=1,3)
2410          call wrtout(ab_out,message,'COLL')
2411        end do
2412      end if
2413 
2414      if (ngrvdw==natom) then
2415        write(message, '(a)' ) ' Van der Waals DFT-D contribution to reduced grads'
2416        call wrtout(ab_out,message,'COLL')
2417        do iatom=1,natom
2418          write(message,format01020) iatom,(grvdw(mu,iatom),mu=1,3)
2419          call wrtout(ab_out,message,'COLL')
2420        end do
2421      end if
2422 
2423      write(message, '(a)' ) ' residual contribution to red. grads'
2424      call wrtout(ab_out,message,'COLL')
2425      do iatom=1,natom
2426        write(message,format01020) iatom,(gresid(mu,iatom),mu=1,3)
2427        call wrtout(ab_out,message,'COLL')
2428      end do
2429 
2430    end if
2431 
2432 !  Compute mean squared deviation from starting coords
2433    devsqr=0.0_dp
2434    do iatom=1,natom
2435      do mu=1,3
2436        devsqr=devsqr+(xred(mu,iatom)-start(mu,iatom))**2
2437      end do
2438    end do
2439 
2440 !  When shift is nonnegligible then print values
2441    if (devsqr>1.d-14) then
2442      write(message, '(a,1p,e12.4,3x,a)' ) &
2443 &     ' rms coord change=',sqrt(devsqr/dble(3*natom)),&
2444 &     'atom, delta coord (reduced):'
2445      call wrtout(ab_out,message,'COLL')
2446      do iatom=1,natom
2447        write(message, '(1x,i5,2x,3f20.12)' ) iatom,&
2448 &       (xred(mu,iatom)-start(mu,iatom),mu=1,3)
2449        call wrtout(ab_out,message,'COLL')
2450      end do
2451    end if
2452 
2453 !  Write out stress results
2454    if (prtstr>0) then
2455      write(message, '(a,a)' ) ch10,&
2456 &     ' Cartesian components of stress tensor (hartree/bohr^3)'
2457      call wrtout(ab_out,message,'COLL')
2458      call wrtout(std_out,  message,'COLL')
2459 
2460      write(message, '(a,1p,e16.8,a,1p,e16.8)' ) &
2461 &     '  sigma(1 1)=',strten(1),'  sigma(3 2)=',strten(4)
2462      call wrtout(ab_out,message,'COLL')
2463      call wrtout(std_out,  message,'COLL')
2464      write(message, '(a,1p,e16.8,a,1p,e16.8)' ) &
2465 &     '  sigma(2 2)=',strten(2),'  sigma(3 1)=',strten(5)
2466      call wrtout(ab_out,message,'COLL')
2467      call wrtout(std_out,  message,'COLL')
2468      write(message, '(a,1p,e16.8,a,1p,e16.8)' ) &
2469 &     '  sigma(3 3)=',strten(3),'  sigma(2 1)=',strten(6)
2470      call wrtout(ab_out,message,'COLL')
2471      call wrtout(std_out,  message,'COLL')
2472 
2473 !    Also output the pressure (minus one third the trace of the stress
2474 !    tensor.
2475      write(message, '(a,a,es12.4,a)' ) ch10,&
2476 &     '-Cartesian components of stress tensor (GPa)         [Pressure=',&
2477 &     -(strten(1)+strten(2)+strten(3))*HaBohr3_GPa/3.0_dp,' GPa]'
2478 
2479      call wrtout(ab_out,message,'COLL')
2480      call wrtout(std_out,  message,'COLL')
2481 
2482      write(message, '(a,1p,e16.8,a,1p,e16.8)' ) &
2483 &     '- sigma(1 1)=',strten(1)*HaBohr3_GPa,&
2484 &     '  sigma(3 2)=',strten(4)*HaBohr3_GPa
2485      call wrtout(ab_out,message,'COLL')
2486      call wrtout(std_out,  message,'COLL')
2487      write(message, '(a,1p,e16.8,a,1p,e16.8)' ) &
2488 &     '- sigma(2 2)=',strten(2)*HaBohr3_GPa,&
2489 &     '  sigma(3 1)=',strten(5)*HaBohr3_GPa
2490      call wrtout(ab_out,message,'COLL')
2491      call wrtout(std_out,  message,'COLL')
2492      write(message, '(a,1p,e16.8,a,1p,e16.8)' ) &
2493 &     '- sigma(3 3)=',strten(3)*HaBohr3_GPa,&
2494 &     '  sigma(2 1)=',strten(6)*HaBohr3_GPa
2495      call wrtout(ab_out,message,'COLL')
2496      call wrtout(std_out,  message,'COLL')
2497    end if
2498 
2499 !  Last end if above refers to iscf > 0
2500  end if
2501 
2502 !DEBUG
2503 !write(std_out,*)' clnup2 : exit '
2504 !ENDDEBUG
2505 
2506 end subroutine clnup2

m_gstate/gstate [ Functions ]

[ Top ] [ m_gstate ] [ Functions ]

NAME

 gstate

FUNCTION

 Primary routine for conducting DFT calculations by CG minimization.

INPUTS

  args_gs<type(args_gs_type)>=various input arguments for the GS calculation
                              Possibly different from dtset
  codvsn=code version
  cpui=initial CPU 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 =number of k points treated by this processor (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
  initialized= 0 for the first GS calculation (not initialized), else 1
  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
  scf_history <type(scf_history_type)>=arrays obtained from previous SCF cycles
  vel(3,natom)=value of velocity
  vel_cell(3,3)=value of cell parameters velocity
  wvl <type(wvl_data)>=all wavelets data
  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.
 In case of wavelets:
 --------------------
    - Only the usual FFT grid (defined by wvl_crmult) is used.
      It is defined by nfft, ngfft, mgfft, ... This is strictly not
      an FFT grid since its dimensions are not suited for FFTs. They are
      defined by wvl_setngfft().
      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

      gstateimg

CHILDREN

      wrtout

SOURCE

 234 subroutine gstate(args_gs,acell,codvsn,cpui,dtfil,dtset,iexit,initialized,&
 235 &                 mpi_enreg,npwtot,occ,pawang,pawrad,pawtab,&
 236 &                 psps,results_gs,rprim,scf_history,vel,vel_cell,wvl,xred)
 237 
 238 
 239 !This section has been created automatically by the script Abilint (TD).
 240 !Do not modify the following lines by hand.
 241 #undef ABI_FUNC
 242 #define ABI_FUNC 'gstate'
 243 !End of the abilint section
 244 
 245  implicit none
 246 
 247 !Arguments ------------------------------------
 248 !scalars
 249  integer,intent(inout) :: iexit,initialized
 250  real(dp),intent(in) :: cpui
 251  character(len=6),intent(in) :: codvsn
 252  type(MPI_type),intent(inout) :: mpi_enreg
 253  type(args_gs_type),intent(in) :: args_gs
 254  type(datafiles_type),intent(inout) :: dtfil
 255  type(dataset_type),intent(inout) :: dtset
 256  type(pawang_type),intent(inout) :: pawang
 257  type(pseudopotential_type),intent(inout) :: psps
 258  type(results_gs_type),intent(inout) :: results_gs
 259  type(scf_history_type),target,intent(inout) :: scf_history
 260  type(wvl_data),intent(inout) :: wvl
 261 !arrays
 262  integer,intent(out) :: npwtot(dtset%nkpt)
 263  real(dp),intent(inout) :: acell(3),occ(dtset%mband*dtset%nkpt*dtset%nsppol)
 264  real(dp),intent(inout) :: rprim(3,3),vel(3,dtset%natom),vel_cell(3,3),xred(3,dtset%natom)
 265  type(pawrad_type),intent(inout) :: pawrad(psps%ntypat*psps%usepaw)
 266  type(pawtab_type),intent(inout) :: pawtab(psps%ntypat*psps%usepaw)
 267 
 268 !Local variables-------------------------------
 269 !Define file format for different type of files. Presently,
 270 !only one file format is supported for each type of files, but this might change soon ...
 271 !2   for wavefunction file, new format (version 2.0 and after)    (fform)   NOT USED
 272 !52  for density rho(r)       (fformr)
 273 !102 for potential V(r) file. (fformv)  NOT USED
 274 !scalars
 275  integer,parameter :: formeig=0,level=101,response=0 ,cplex1=1, master=0
 276  integer :: ndtpawuj=0  ! Cannot use parameter because scfargs points to this! Have to get rid of pointers to scalars!
 277 #if defined HAVE_BIGDFT
 278  integer :: icoulomb
 279 #endif
 280  integer :: accessfil,ask_accurate,bantot,choice,comm_psp,fform
 281  integer :: gnt_option,gscase,iatom,idir,ierr,ii,indx,jj,kk,ios,itypat
 282  integer :: ixfh,izero,mband_cprj,mcg,mcprj,me,mgfftf,mpert,msize,mu,my_natom,my_nspinor
 283  integer :: nblok,ncpgr,nfftf,nfftot,npwmin
 284  integer :: openexit,option,optorth,psp_gencond,conv_retcode
 285  integer :: pwind_alloc,rdwrpaw,comm,tim_mkrho,use_sc_dmft
 286  integer :: cnt,spin,band,ikpt,usecg,usecprj
 287  real(dp) :: cpus,ecore,ecut_eff,ecutdg_eff,etot,fermie
 288  real(dp) :: gsqcut_eff,gsqcut_shp,gsqcutc_eff,hyb_range_fock,residm,ucvol
 289  logical :: read_wf_or_den,has_to_init,call_pawinit,write_wfk
 290  logical :: wvlbigdft=.false.,wvl_debug=.false.
 291  character(len=500) :: message
 292  character(len=fnlen) :: ddbnm,dscrpt,filnam,wfkfull_path
 293  real(dp) :: fatvshift
 294  type(crystal_t) :: cryst
 295  type(ebands_t) :: bstruct,ebands
 296  type(t_tetrahedron) :: tetra
 297  type(efield_type) :: dtefield
 298  type(electronpositron_type),pointer :: electronpositron
 299  type(hdr_type) :: hdr,hdr_den
 300  type(macro_uj_type) :: dtpawuj(0)
 301  type(orbmag_type) :: dtorbmag
 302  type(paw_dmft_type) :: paw_dmft
 303  type(pawfgr_type) :: pawfgr
 304  type(recursion_type) ::rec_set
 305  type(wffile_type) :: wff1,wffnew,wffnow
 306  type(ab_xfh_type) :: ab_xfh
 307  type(ddb_type) :: ddb
 308  type(ddb_hdr_type) :: ddb_hdr
 309  type(scfcv_t) :: scfcv_args
 310 !arrays
 311  integer :: ngfft(18),ngfftf(18)
 312  integer,allocatable :: atindx(:),atindx1(:),indsym(:,:,:),dimcprj_srt(:)
 313  integer,allocatable :: irrzon(:,:,:),kg(:,:),nattyp(:),symrec(:,:,:)
 314  integer,allocatable,target :: npwarr(:)
 315  integer,pointer :: npwarr_(:),pwind(:,:,:)
 316  real(dp) :: efield_band(3),gmet(3,3),gmet_for_kg(3,3),gprimd(3,3),gprimd_for_kg(3,3)
 317  real(dp) :: rmet(3,3),rprimd(3,3),rprimd_for_kg(3,3),tsec(2)
 318  real(dp),allocatable :: doccde(:)
 319  real(dp),allocatable :: ph1df(:,:),phnons(:,:,:),resid(:),rhowfg(:,:)
 320  real(dp),allocatable :: rhowfr(:,:),spinat_dum(:,:),start(:,:),work(:)
 321  real(dp),allocatable :: ylm(:,:),ylmgr(:,:,:)
 322  real(dp),pointer :: cg(:,:),eigen(:),pwnsfac(:,:),rhog(:,:),rhor(:,:)
 323  real(dp),pointer :: taug(:,:),taur(:,:),xred_old(:,:)
 324  type(pawrhoij_type),pointer :: pawrhoij(:)
 325  type(coulomb_operator) :: kernel_dummy
 326  type(pawcprj_type),allocatable :: cprj(:,:)
 327 
 328 ! ***********************************************************************
 329 
 330  DBG_ENTER("COLL")
 331 
 332  call timab(32,1,tsec)
 333  call timab(33,3,tsec)
 334 
 335  call status(0,dtfil%filstat,iexit,level,'enter')
 336 
 337 !###########################################################
 338 !### 01. Initializations XML, MPI, WVL, etc
 339 
 340 !Init MPI data
 341  comm=mpi_enreg%comm_cell; me=xmpi_comm_rank(comm)
 342 
 343 !Set up MPI information from the dataset
 344  my_natom=mpi_enreg%my_natom
 345 
 346 !Set up information when wavelets are in use
 347  if (dtset%usewvl == 1) then
 348 
 349 !  If usewvl: wvlbigdft indicates that the BigDFT workflow will be followed
 350    wvlbigdft=(dtset%wvl_bigdft_comp==1)
 351 
 352 !  Default value, to be set-up elsewhere.
 353    wvl%descr%h(:) = dtset%wvl_hgrid
 354 
 355 #if defined HAVE_BIGDFT
 356    wvl%descr%paw%usepaw=psps%usepaw
 357    wvl%descr%paw%natom=dtset%natom
 358 #endif
 359 
 360 !  We set the atom-related internal wvl variables.
 361    call wvl_descr_atoms_set(acell, dtset%icoulomb, dtset%natom, &
 362 &   dtset%ntypat, dtset%typat, wvl%descr)
 363    if(dtset%usepaw==0) then
 364 !    nullify PAW proj_G in NC case:
 365 #if defined HAVE_BIGDFT
 366      ABI_DATATYPE_ALLOCATE(wvl%projectors%G,(dtset%ntypat))
 367      do itypat=1,dtset%ntypat
 368        call nullify_gaussian_basis(wvl%projectors%G(itypat))
 369      end do
 370 #endif
 371    end if
 372 
 373    wvl%descr%exctxpar = "OP2P"
 374  end if
 375 
 376  if (me == master .and. dtset%prtxml == 1) then
 377 !  gstate() will handle a dataset, so we output the dataSet markup.
 378    write(ab_xml_out, "(A)") '  <dataSet>'
 379  end if
 380 
 381 !Define FFT grid(s) sizes (be careful !)
 382 !See NOTES in the comments at the beginning of this file.
 383  call pawfgr_init(pawfgr,dtset,mgfftf,nfftf,ecut_eff,ecutdg_eff,ngfft,ngfftf)
 384 
 385 !Structured debugging if prtvol==-level
 386  if(dtset%prtvol==-level)then
 387    write(message,'(80a,a,a)')  ('=',ii=1,80),ch10,' gstate : enter , debug mode '
 388    call wrtout(std_out,message,'COLL')
 389  end if
 390 
 391 !###########################################################
 392 !### 02. Calls setup1, kpgio, initylmg
 393 
 394  ecore=zero
 395  results_gs%pel(1:3)   =zero
 396  results_gs%grchempottn(:,:)=zero
 397  results_gs%grewtn(:,:)=zero
 398 !MT Feb 2012: I dont know why but grvdw has to be allocated
 399 !when using BigDFT to ensure success on inca_gcc44_sdebug
 400  if ((dtset%vdw_xc>=5.and.dtset%vdw_xc<=7).or.dtset%usewvl==1) then
 401    results_gs%ngrvdw=dtset%natom
 402    if (allocated(results_gs%grvdw)) then
 403      ABI_DEALLOCATE(results_gs%grvdw)
 404    end if
 405    ABI_ALLOCATE(results_gs%grvdw,(3,dtset%natom))
 406    results_gs%grvdw(:,:)=zero
 407  end if
 408  call energies_init(results_gs%energies)
 409 
 410 !Set up for iterations
 411  call setup1(acell,bantot,dtset,&
 412 & ecutdg_eff,ecut_eff,gmet,gprimd,gsqcut_eff,gsqcutc_eff,&
 413 & dtset%natom,ngfftf,ngfft,dtset%nkpt,dtset%nsppol,&
 414 & response,rmet,rprim,rprimd,ucvol,psps%usepaw)
 415 
 416 !In some cases (e.g. getcell/=0), the plane wave vectors have
 417 ! to be generated from the original simulation cell
 418  rprimd_for_kg=rprimd
 419  if (dtset%getcell/=0.and.dtset%usewvl==0) rprimd_for_kg=args_gs%rprimd_orig
 420  call matr3inv(rprimd_for_kg,gprimd_for_kg)
 421  gmet_for_kg=matmul(transpose(gprimd_for_kg),gprimd_for_kg)
 422 
 423 !Set up the basis sphere of planewaves
 424  ABI_ALLOCATE(npwarr,(dtset%nkpt))
 425  if (dtset%usewvl == 0 .and. dtset%tfkinfunc /= 2) then
 426    ABI_ALLOCATE(kg,(3,dtset%mpw*dtset%mkmem))
 427    call kpgio(ecut_eff,dtset%exchn2n3d,gmet_for_kg,dtset%istwfk,kg, &
 428 &   dtset%kptns,dtset%mkmem,dtset%nband,dtset%nkpt,'PERS',mpi_enreg,&
 429 &   dtset%mpw,npwarr,npwtot,dtset%nsppol)
 430    call bandfft_kpt_init1(bandfft_kpt,dtset%istwfk,kg,dtset%mgfft,dtset%mkmem,mpi_enreg,&
 431 &   dtset%mpw,dtset%nband,dtset%nkpt,npwarr,dtset%nsppol)
 432  else
 433    ABI_ALLOCATE(kg,(0,0))
 434    npwarr(:) = 0
 435    npwtot(:) = 0
 436  end if
 437 
 438  if(dtset%wfoptalg == 1 .and. psps%usepaw == 1) then
 439    call init_invovl(dtset%nkpt)
 440  end if
 441 
 442  if(dtset%use_gemm_nonlop == 1 .and. dtset%use_gpu_cuda/=1) then
 443    ! set global variable
 444    gemm_nonlop_use_gemm = .true.
 445    call init_gemm_nonlop(dtset%nkpt)
 446  else
 447    gemm_nonlop_use_gemm = .false.
 448  end if
 449 
 450 !Set up the Ylm for each k point
 451  if ( dtset%tfkinfunc /= 2) then
 452    ABI_ALLOCATE(ylm,(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm))
 453    ABI_ALLOCATE(ylmgr,(dtset%mpw*dtset%mkmem,3,psps%mpsang*psps%mpsang*psps%useylm))
 454    if (psps%useylm==1) then
 455      option=0
 456      if (dtset%prtstm==0.and.dtset%iscf>0.and.dtset%positron/=1) option=1 ! compute gradients of YLM
 457      if (dtset%berryopt==4 .and. dtset%optstress /= 0 .and. psps%usepaw==1) option = 1 ! compute gradients of YLM
 458      call initylmg(gprimd,kg,dtset%kptns,dtset%mkmem,mpi_enreg,&
 459 &     psps%mpsang,dtset%mpw,dtset%nband,dtset%nkpt,&
 460 &     npwarr,dtset%nsppol,option,rprimd,ylm,ylmgr)
 461    end if
 462  else
 463    ABI_ALLOCATE(ylm,(0,0))
 464    ABI_ALLOCATE(ylmgr,(0,0,0))
 465  end if
 466 
 467 !SCF history management (allocate it at first call)
 468  if (initialized==0) then
 469 !  This call has to be done before any use of SCF history
 470    usecg=0
 471    if(dtset%extrapwf>0 .or. dtset%imgwfstor==1)usecg=1
 472    call scf_history_init(dtset,mpi_enreg,usecg,scf_history)
 473  end if
 474  has_to_init=(initialized==0.or.scf_history%history_size<0)
 475 
 476  call timab(33,2,tsec)
 477  call timab(701,3,tsec)
 478 
 479 !###########################################################
 480 !### 03. Calls pspini
 481 
 482 !Open and read pseudopotential files
 483  comm_psp=mpi_enreg%comm_cell;if (dtset%usewvl==1) comm_psp=mpi_enreg%comm_wvl
 484  if (dtset%nimage>1) psps%mixalch(:,:)=args_gs%mixalch(:,:) ! mixalch can evolve for some image algos
 485  call pspini(dtset,dtfil,ecore,psp_gencond,gsqcutc_eff,gsqcut_eff,&
 486 & pawrad,pawtab,psps,rprimd,comm_mpi=comm_psp)
 487 
 488  call timab(701,2,tsec)
 489  call timab(33,3,tsec)
 490  if (psps%usepaw==1.and.dtset%usefock==1)then
 491    do itypat=1,dtset%ntypat
 492      if(pawtab(itypat)%has_fock==0) then
 493        write(message, '(a)' )&
 494 &       'The PAW data file does not contain Fock information. Change the PAW data file.'
 495        MSG_BUG(message)
 496      end if
 497    end do
 498  end if
 499 !In case of isolated computations, ecore must set to zero
 500 !because its contribution is counted in the ewald energy as the ion-ion interaction.
 501  if (dtset%icoulomb == 1) ecore = zero
 502 
 503 !WVL - Now that psp data are available, we compute rprimd, acell... from the atomic positions.
 504  if (dtset%usewvl == 1) then
 505    call wvl_descr_psp_set(trim(dtfil%filnam_ds(3))//"_OCCUP",dtset%nsppol,psps,dtset%spinat,wvl%descr)
 506    call wvl_setBoxGeometry(dtset%prtvol, psps%gth_params%radii_cf, rprimd, xred, &
 507 &   wvl%descr, dtset%wvl_crmult, dtset%wvl_frmult)
 508    call mkradim(acell,rprim,rprimd)
 509    rprimd_for_kg=rprimd
 510    call wvl_denspot_set(wvl%den, psps%gth_params, dtset%ixc, dtset%natom, dtset%nsppol, rprimd, &
 511 &   wvl%descr, dtset%wvl_crmult, dtset%wvl_frmult, mpi_enreg%comm_wvl, xred)
 512 !  TODO: to be moved in a routine.
 513 #if defined HAVE_BIGDFT
 514    if (wvl%descr%atoms%astruct%geocode == "F") then
 515      icoulomb = 1
 516    else if (wvl%descr%atoms%astruct%geocode == "S") then
 517      icoulomb = 2
 518    else
 519      icoulomb = 0
 520    end if
 521 !  calculation of the Poisson kernel anticipated to reduce memory peak for small systems
 522    call psolver_kernel( wvl%den%denspot%dpbox%hgrids, 1, icoulomb, mpi_enreg%me_wvl, wvl%den%denspot%pkernel , &
 523 &   mpi_enreg%comm_wvl, wvl%den%denspot%dpbox%ndims, mpi_enreg%nproc_wvl, dtset%nscforder)
 524    nullify(wvl%den%denspot%pkernelseq%kernel)
 525    !call copy_coulomb_operator(wvl%den%denspot%pkernel,wvl%den%denspot%pkernelseq,ABI_FUNC)
 526 !  Associate the denspot distribution into mpi_enreg.
 527    mpi_enreg%nscatterarr  => wvl%den%denspot%dpbox%nscatterarr
 528    mpi_enreg%ngatherarr   => wvl%den%denspot%dpbox%ngatherarr
 529    mpi_enreg%ngfft3_ionic =  wvl%den%denspot%dpbox%n3pi
 530    call wvl_setngfft(mpi_enreg%me_wvl, dtset%mgfft, dtset%nfft, &
 531 &   dtset%ngfft, mpi_enreg%nproc_wvl, wvl%den%denspot%dpbox%ndims(1), &
 532 &   wvl%den%denspot%dpbox%ndims(2),wvl%den%denspot%dpbox%ndims(3),&
 533 &   wvl%den%denspot%dpbox%nscatterarr(mpi_enreg%me_wvl, 1))
 534 #endif
 535    nfftf     = dtset%nfft
 536    mgfftf    = dtset%mgfft
 537    ngfftf(:) = dtset%ngfft(:)
 538 !  Recalculate gprimd
 539    call matr3inv(rprimd,gprimd)
 540 !  PAW section
 541    if(psps%usepaw==1) then
 542 !    Reinitialize Pawfgr with new values of ngfft
 543      call pawfgr_destroy(pawfgr)
 544      call pawfgr_init(pawfgr,dtset,mgfftf,nfftf,ecut_eff,ecutdg_eff,ngfft,ngfftf)
 545 !    fill wvl objects from paw objects
 546 !    wvl%descr%paw%usepaw=dtset%usepaw
 547      call paw2wvl(pawtab,wvl%projectors,wvl%descr)
 548    end if
 549  else if (dtset%icoulomb /= 0) then
 550 #if defined HAVE_BIGDFT
 551    if (dtset%ixc < 0) then
 552      call xc_init(wvl%den%denspot%xc, dtset%ixc, XC_MIXED, dtset%nsppol)
 553    else
 554      call xc_init(wvl%den%denspot%xc, dtset%ixc, XC_ABINIT, dtset%nsppol)
 555    end if
 556 #endif
 557  end if
 558 
 559 !Initialize band structure datatype
 560  if (dtset%paral_kgb/=0) then     !  We decide to store total npw in bstruct,
 561    ABI_ALLOCATE(npwarr_,(dtset%nkpt))
 562    npwarr_(:)=npwarr(:)
 563    call xmpi_sum(npwarr_,mpi_enreg%comm_bandfft,ierr)
 564  else
 565    npwarr_ => npwarr
 566  end if
 567 
 568  bstruct = ebands_from_dtset(dtset, npwarr_)
 569 
 570  if (dtset%paral_kgb/=0)  then
 571    ABI_DEALLOCATE(npwarr_)
 572  end if
 573  nullify(npwarr_)
 574 
 575 !Initialize PAW atomic occupancies
 576  if (scf_history%history_size>=0) then
 577    pawrhoij => scf_history%pawrhoij_last
 578  else
 579    ABI_DATATYPE_ALLOCATE(pawrhoij,(my_natom*psps%usepaw))
 580  end if
 581  if (psps%usepaw==1.and.has_to_init) then
 582    call initrhoij(dtset%pawcpxocc,dtset%lexexch,&
 583 &   dtset%lpawu,my_natom,dtset%natom,dtset%nspden,dtset%nspinor,dtset%nsppol,&
 584 &   dtset%ntypat,pawrhoij,dtset%pawspnorb,pawtab,dtset%spinat,dtset%typat,&
 585 &   comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
 586  end if
 587 
 588 !Initialize header
 589  gscase=0
 590  call hdr_init(bstruct,codvsn,dtset,hdr,pawtab,gscase,psps,wvl%descr,&
 591 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
 592 
 593 !Clean band structure datatype (should use it more in the future !)
 594  call ebands_free(bstruct)
 595 
 596 !Update header, with evolving variables, when available
 597 !Here, rprimd, xred and occ are available
 598  etot=hdr%etot ; fermie=hdr%fermie ; residm=hdr%residm
 599  call hdr_update(hdr,bantot,etot,fermie,&
 600 & residm,rprimd,occ,pawrhoij,xred,args_gs%amu,&
 601 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
 602 
 603  ! PW basis set: test if the problem is ill-defined.
 604  if (dtset%usewvl == 0 .and. dtset%tfkinfunc /= 2) then
 605    npwmin=minval(hdr%npwarr(:))
 606    if (dtset%mband > npwmin) then
 607      ! No way we can solve the problem. Abort now!
 608      write(message,"(2(a,i0),4a)")&
 609 &     "Number of bands nband= ",dtset%mband," > number of planewaves npw= ",npwmin,ch10,&
 610 &     "The number of eigenvectors cannot be greater that the size of the Hamiltonian!",ch10,&
 611 &     "Action: decrease nband or, alternatively, increase ecut"
 612      if (dtset%ionmov/=23) then
 613        MSG_ERROR(message)
 614      else
 615        MSG_WARNING(message)
 616      end if
 617 
 618    else if (dtset%mband >= 0.9 * npwmin) then
 619      ! Warn the user
 620      write(message,"(a,i0,a,f6.1,4a)")&
 621 &     "Number of bands nband= ",dtset%mband," >= 0.9 * maximum number of planewaves= ",0.9*npwmin,ch10,&
 622 &     "The problem is ill-defined and the GS algorithm will show numerical instabilities!",ch10,&
 623 &     "Assume experienced user. Execution will continue."
 624      MSG_WARNING(message)
 625    end if
 626  end if
 627 
 628 !###########################################################
 629 !### 04. Symmetry operations when nsym>1
 630 
 631 !Do symmetry stuff only for nsym>1
 632  if (dtset%usewvl == 0) then
 633    nfftot=ngfft(1)*ngfft(2)*ngfft(3)
 634  else
 635 #if defined HAVE_BIGDFT
 636    nfftot=product(wvl%den%denspot%dpbox%ndims)
 637 #else
 638    BIGDFT_NOTENABLED_ERROR()
 639 #endif
 640  end if
 641  ABI_ALLOCATE(irrzon,(nfftot**(1-1/dtset%nsym),2,(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4)))
 642  ABI_ALLOCATE(phnons,(2,nfftot**(1-1/dtset%nsym),(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4)))
 643  ABI_ALLOCATE(indsym,(4,dtset%nsym,dtset%natom))
 644  ABI_ALLOCATE(symrec,(3,3,dtset%nsym))
 645  irrzon(:,:,:)=0
 646  phnons(:,:,:)=zero
 647  indsym(:,:,:)=0
 648  symrec(:,:,:)=0
 649 
 650  if (dtset%nsym>1) then
 651    call setsym(indsym,irrzon,dtset%iscf,dtset%natom,&
 652 &   nfftot,ngfft,dtset%nspden,dtset%nsppol,dtset%nsym,&
 653 &   phnons,dtset%symafm,symrec,dtset%symrel,dtset%tnons,dtset%typat,xred)
 654 
 655 !  Make sure dtset%iatfix does not break symmetry
 656    call fixsym(dtset%iatfix,indsym,dtset%natom,dtset%nsym)
 657  else
 658 !  The symrec array is used by initberry even in case nsym = 1
 659    symrec(:,:,1) = 0
 660    symrec(1,1,1) = 1 ; symrec(2,2,1) = 1 ; symrec(3,3,1) = 1
 661  end if
 662  if (dtset%usewvl == 1) then
 663    call wvl_descr_atoms_set_sym(wvl%descr, dtset%efield, irrzon, dtset%nsppol, &
 664 &   dtset%nsym, phnons, dtset%symafm, dtset%symrel, dtset%tnons, dtset%tolsym)
 665 #if defined HAVE_BIGDFT
 666    wvl%den%symObj = wvl%descr%atoms%astruct%sym%symObj
 667 #endif
 668  end if
 669 
 670 !###########################################################
 671 !### 05. Calls inwffil
 672 
 673  ! if paral_kgb == 0, it may happen that some processors are idle (no entry in proc_distrb)
 674  ! but mkmem == nkpt and this can cause integer overflow in mcg or allocation error.
 675  ! Here we count the number of states treated by the proc. if cnt == 0, mcg is then set to 0.
 676  cnt = 0
 677  do spin=1,dtset%nsppol
 678    do ikpt=1,dtset%nkpt
 679      do band=1,dtset%nband(ikpt + (spin-1) * dtset%nkpt)
 680        if (.not. proc_distrb_cycle(mpi_enreg%proc_distrb, ikpt, band, band, spin, mpi_enreg%me_kpt)) cnt = cnt + 1
 681      end do
 682    end do
 683  end do
 684 
 685  my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
 686  mcg=dtset%mpw*my_nspinor*dtset%mband*dtset%mkmem*dtset%nsppol
 687  if (cnt == 0) then
 688    mcg = 0
 689    write(message,"(2(a,i0))")"rank: ",mpi_enreg%me, "does not have wavefunctions to treat. Setting mcg to: ",mcg
 690    MSG_WARNING(message)
 691  end if
 692 
 693  if (dtset%usewvl == 0 .and. dtset%mpw > 0 .and. cnt /= 0)then
 694    if (my_nspinor*dtset%mband*dtset%mkmem*dtset%nsppol > floor(real(HUGE(0))/real(dtset%mpw) )) then
 695      write (message,'(9a)')&
 696 &     "Default integer is not wide enough to store the size of the wavefunction array (mcg).",ch10,&
 697 &     "This usually happens when paral_kgb == 0 and there are not enough procs to distribute kpts and spins",ch10,&
 698 &     "Action: if paral_kgb == 0, use nprocs = nkpt * nsppol to reduce the memory per node.",ch10,&
 699 &     "If this does not solve the problem, use paral_kgb 1 with nprocs > nkpt * nsppol and use npfft/npband/npspinor",ch10,&
 700 &     "to decrease the memory requirements. Consider also OpenMP threads."
 701      MSG_ERROR_NOSTOP(message,ii)
 702      write (message,'(5(a,i0), 2a)')&
 703 &     "my_nspinor: ",my_nspinor, ", mpw: ",dtset%mpw, ", mband: ",dtset%mband,&
 704 &     ", mkmem: ",dtset%mkmem, ", nsppol: ",dtset%nsppol,ch10,&
 705 &     'Note: Compiling with large int (int64) requires a full software stack (MPI/FFTW/BLAS...) compiled in int64 mode'
 706      MSG_ERROR(message)
 707    end if
 708  end if
 709 
 710  if (dtset%imgwfstor==1) then
 711    cg => scf_history%cg(:,:,1)
 712    eigen => scf_history%eigen(:,1)
 713  else
 714    ABI_STAT_ALLOCATE(cg,(2,mcg), ierr)
 715    ABI_CHECK(ierr==0, "out of memory in cg")
 716    ABI_ALLOCATE(eigen,(dtset%mband*dtset%nkpt*dtset%nsppol))
 717  end if
 718 
 719  ABI_ALLOCATE(resid,(dtset%mband*dtset%nkpt*dtset%nsppol))
 720  eigen(:)=zero ; resid(:)=zero
 721 !mpi_enreg%paralbd=0 ; ask_accurate=0
 722  ask_accurate=0
 723 
 724 !WVL - Branching, allocating wavefunctions as wavelets.
 725  if (dtset%usewvl == 1) then
 726    call wvl_wfs_lr_copy(wvl%wfs, wvl%descr)
 727 !  Create access arrays for wavefunctions and allocate wvl%wfs%psi (other arrays are left unallocated).
 728    call wvl_wfs_set(dtset%strprecon,dtset%spinmagntarget, dtset%kpt, mpi_enreg%me_wvl,&
 729 &   dtset%natom, sum(dtset%nband), &
 730 &   dtset%nkpt, mpi_enreg%nproc_wvl, dtset%nspinor, dtset%nsppol, dtset%nwfshist, occ, &
 731 &   psps, rprimd, wvl%wfs, dtset%wtk, wvl%descr, dtset%wvl_crmult, dtset%wvl_frmult, &
 732 &   xred)
 733 !  We transfer wavelets information to the hdr structure.
 734 #if defined HAVE_BIGDFT
 735    call local_potential_dimensions(mpi_enreg%me_wvl,wvl%wfs%ks%lzd,wvl%wfs%ks%orbs,wvl%den%denspot%xc,&
 736 &   wvl%den%denspot%dpbox%ngatherarr(0,1))
 737    hdr%nwvlarr(1) = wvl%wfs%ks%lzd%Glr%wfd%nvctr_c
 738    hdr%nwvlarr(2) = 7 * wvl%wfs%ks%lzd%Glr%wfd%nvctr_f
 739 #endif
 740 !  Create access arrays for projectors and allocate them.
 741 !  Compute projectors from each atom.
 742    call wvl_projectors_set(mpi_enreg%me_wvl, dtset%natom, wvl%projectors, psps, rprimd, &
 743 &   wvl%wfs, wvl%descr, dtset%wvl_frmult, xred)
 744  end if
 745 
 746  read_wf_or_den=(dtset%iscf<=0.or.dtfil%ireadwf/=0.or.(dtfil%ireadden/=0.and.dtset%positron<=0))
 747  read_wf_or_den=(read_wf_or_den.and.has_to_init)
 748 
 749 !RECURSION -  initialization
 750  if(has_to_init .and. dtset%userec==1) then
 751    call InitRec(dtset,mpi_enreg,rec_set,rmet,maxval(psps%indlmn(3,:,:)))
 752  end if
 753 
 754 !LOTF - initialization
 755 #if defined HAVE_LOTF
 756  if(has_to_init .and. dtset%ionmov==23) then
 757    call lotfparam_init(dtset%natom,dtset%lotf_version,1,&
 758 &   dtset%lotf_nitex,dtset%lotf_nneigx,&
 759 &   dtset%lotf_classic,1,1)
 760  end if
 761 #endif
 762 
 763 !Initialize wavefunctions.
 764  if(dtset%imgwfstor==1 .and. initialized==1)then
 765    cg(:,:)=scf_history%cg(:,:,1)
 766    eigen(:)=scf_history%eigen(:,1)
 767  else if(dtset%tfkinfunc /=2) then
 768 !if(dtset%tfkinfunc /=2) then
 769    wff1%unwff=dtfil%unwff1
 770    optorth=1   !if (psps%usepaw==1) optorth=0
 771    if(psps%usepaw==1 .and. dtfil%ireadwf==1)optorth=0
 772    hdr%rprimd=rprimd_for_kg ! We need the rprimd that was used to generate de G vectors
 773    call inwffil(ask_accurate,cg,dtset,dtset%ecut,ecut_eff,eigen,&
 774 &   dtset%exchn2n3d,formeig,hdr,dtfil%ireadwf,dtset%istwfk,kg,&
 775 &   dtset%kptns,dtset%localrdwf,dtset%mband,mcg,dtset%mkmem,mpi_enreg,&
 776 &   dtset%mpw,dtset%nband,ngfft,dtset%nkpt,npwarr,&
 777 &   dtset%nsppol,dtset%nsym,occ,optorth,dtset%symafm,&
 778 &   dtset%symrel,dtset%tnons,dtfil%unkg,wff1,wffnow,dtfil%unwff1,&
 779 &   dtfil%fnamewffk,wvl)
 780    hdr%rprimd=rprimd
 781  end if
 782 
 783  if (psps%usepaw==1.and.dtfil%ireadwf==1)then
 784    call pawrhoij_copy(hdr%pawrhoij,pawrhoij,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
 785 !  Has to update header again (because pawrhoij has changed)  -  MT 2007-10-22: Why ?
 786 !  call hdr_update(hdr,bantot,etot,fermie,residm,rprimd,occ,pawrhoij,xred,args_gs%amu, &
 787 !  &                  comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
 788  end if
 789 
 790 !###########################################################
 791 !### 06. Operations related to restartxf (Old version)
 792 
 793 !TODO: Remove ab_xfh
 794 !Initialize xf history (should be put in inwffil)
 795  ab_xfh%nxfh=0
 796  if(dtset%restartxf>=1 .and. dtfil%ireadwf==1)then
 797 
 798 !  Should exchange the data about history in parallel localrdwf==0
 799    if(xmpi_paral==1 .and. dtset%localrdwf==0)then
 800      write(message, '(a,a,a)' )&
 801 &     'It is not yet possible to use non-zero restartxf,',ch10,&
 802 &     'in parallel, when localrdwf=0. Sorry for this ...'
 803      MSG_BUG(message)
 804    end if
 805 
 806    ABI_ALLOCATE(ab_xfh%xfhist,(3,dtset%natom+4,2,0))
 807    call outxfhist(ab_xfh,dtset%natom,2,wff1,ios)
 808    ABI_DEALLOCATE(ab_xfh%xfhist)
 809 
 810    if(ios>0)then
 811      write(message,'(a,a,a)')&
 812 &     'An error occurred reading the input wavefunction file,',ch10,&
 813 &     'with restartxf=1.'
 814      MSG_ERROR(message)
 815    else if(ios==0)then
 816      write(message, '(a,a,i4,a)' )ch10,&
 817 &     ' gstate : reading',ab_xfh%nxfh,' (x,f) history pairs from input wf file.'
 818      call wrtout(std_out,message,'COLL')
 819      call wrtout(ab_out,message,'COLL')
 820    end if
 821 !  WARNING : should check that restartxf is not negative
 822 !  WARNING : should check that restartxf /= only when dtfil%ireadwf is activated
 823  end if
 824 
 825 !Allocate the xf history array : takes into account the existing
 826 !pairs, minus those that will be discarded, then those that will
 827 !be computed, governed by dtset%ntime, and some additional pairs
 828 !(needed when it will be possible to use xfhist for move.f)
 829  ab_xfh%mxfh=(ab_xfh%nxfh-dtset%restartxf+1)+dtset%ntime+5
 830  ABI_ALLOCATE(ab_xfh%xfhist,(3,dtset%natom+4,2,ab_xfh%mxfh))
 831  ab_xfh%xfhist(:,:,:,:) = zero
 832 !WARNING : should check that the number of atoms in the wf file and natom are the same
 833 
 834 !Initialize the xf history array
 835  if(ab_xfh%nxfh>=dtset%restartxf .and. ab_xfh%nxfh>0)then
 836 !  Eventually skip some of the previous history
 837    if(dtset%restartxf>=2)then
 838      do ixfh=1,dtset%restartxf-1
 839        call WffReadSkipRec(ios,1,wff1)
 840      end do
 841    end if
 842 
 843 !  Read and store the relevant history
 844    ab_xfh%nxfhr=ab_xfh%nxfh-dtset%restartxf+1
 845    call outxfhist(ab_xfh,dtset%natom,3,wff1,ios)
 846  end if
 847 
 848 !Close wff1, if it was ever opened (in inwffil)
 849  if (dtfil%ireadwf==1) then
 850    call WffClose(wff1,ierr)
 851  end if
 852 
 853 !###########################################################
 854 !### 07. Calls setup2
 855 
 856 !Further setup
 857  ABI_ALLOCATE(start,(3,dtset%natom))
 858  call setup2(dtset,npwtot,start,wvl%wfs,xred)
 859 
 860 !Allocation of previous atomic positions
 861  if (scf_history%history_size>=0) then
 862    xred_old => scf_history%xred_last
 863  else
 864    ABI_ALLOCATE(xred_old,(3,dtset%natom))
 865  end if
 866  if (has_to_init) xred_old=xred
 867 
 868 !Timing for initialisation period
 869  call timab(33,2,tsec)
 870  call timab(34,3,tsec)
 871 
 872 !###########################################################
 873 !### 08. Compute new occupation numbers
 874 
 875 !Compute new occupation numbers, in case wavefunctions and eigenenergies
 876 !were read from disk, occupation scheme is metallic (this excludes iscf=-1),
 877 !and occupation numbers are required by iscf
 878  if( dtfil%ireadwf==1 .and. &
 879 & (dtset%occopt>=3.and.dtset%occopt<=8) .and. &
 880 & (dtset%iscf>0 .or. dtset%iscf==-3) .and. dtset%positron/=1 ) then
 881 
 882    ABI_ALLOCATE(doccde,(dtset%mband*dtset%nkpt*dtset%nsppol))
 883 !  Warning : ideally, results_gs%entropy should not be set up here XG 20011007
 884 !  Do not take into account the possible STM bias
 885    call newocc(doccde,eigen,results_gs%energies%entropy,&
 886 &   results_gs%energies%e_fermie,&
 887 &   dtset%spinmagntarget,dtset%mband,dtset%nband,&
 888 &   dtset%nelect,dtset%nkpt,dtset%nspinor,dtset%nsppol,occ,&
 889 &   dtset%occopt,dtset%prtvol,zero,dtset%tphysel,dtset%tsmear,dtset%wtk)
 890    if (dtset%dmftcheck>=0.and.dtset%usedmft>=1.and.(sum(args_gs%upawu(:))>=tol8.or.  &
 891 &   sum(args_gs%jpawu(:))>tol8).and.dtset%dmft_entropy==0) results_gs%energies%entropy=zero
 892    ABI_DEALLOCATE(doccde)
 893 
 894 !  Transfer occupations to bigdft object:
 895    if(dtset%usewvl==1 .and. .not. wvlbigdft) then
 896      call wvl_occ_abi2big(dtset%mband,dtset%nkpt,dtset%nsppol,occ,1,wvl%wfs)
 897 !    call wvl_energies_abi2big(results_gs%energies,wvl%wfs,2)
 898    end if
 899 
 900  else
 901 !  Warning : ideally, results_gs%entropy should not be set up here XG 20011007
 902    results_gs%energies%entropy=zero
 903  end if
 904 
 905 !###########################################################
 906 !### 09. Generate an index table of atoms
 907 
 908 !Definition of atindx array
 909 !Generate an index table of atoms, in order for them to be used type after type.
 910  ABI_ALLOCATE(atindx,(dtset%natom))
 911  ABI_ALLOCATE(atindx1,(dtset%natom))
 912  ABI_ALLOCATE(nattyp,(psps%ntypat))
 913  indx=1
 914  do itypat=1,psps%ntypat
 915    nattyp(itypat)=0
 916    do iatom=1,dtset%natom
 917      if(dtset%typat(iatom)==itypat)then
 918        atindx(iatom)=indx
 919        atindx1(indx)=iatom
 920        indx=indx+1
 921        nattyp(itypat)=nattyp(itypat)+1
 922      end if
 923    end do
 924  end do
 925 
 926 !Compute structure factor phases for current atomic pos:
 927  if ((.not.read_wf_or_den).or.(scf_history%history_size>0.and.has_to_init)) then
 928    ABI_ALLOCATE(ph1df,(2,3*(2*mgfftf+1)*dtset%natom))
 929    call getph(atindx,dtset%natom,ngfftf(1),ngfftf(2),ngfftf(3),ph1df,xred)
 930  end if
 931 
 932 !Here allocation of GPU for vtorho calculations
 933 #if defined HAVE_GPU_CUDA
 934  if (dtset%use_gpu_cuda==1) then
 935    call alloc_hamilt_gpu(atindx1,dtset,gprimd,mpi_enreg,nattyp,npwarr,2,psps,dtset%use_gpu_cuda)
 936  end if
 937 #endif
 938 
 939 !###########################################################
 940 !### 10. PAW related operations
 941 
 942 !Initialize paw_dmft, even if neither dmft not paw are used
 943 !write(std_out,*) "dtset%usedmft",dtset%usedmft
 944  use_sc_dmft=dtset%usedmft
 945  if(dtset%paral_kgb>0) use_sc_dmft=0
 946  call init_sc_dmft(dtset%nbandkss,dtset%dmftbandi,dtset%dmftbandf,dtset%dmft_read_occnd,dtset%mband,&
 947 & dtset%nband,dtset%nkpt,dtset%nspden, &
 948 & dtset%nspinor,dtset%nsppol,occ,dtset%usedmft,paw_dmft,use_sc_dmft,dtset%dmft_solv,mpi_enreg)
 949  if (paw_dmft%use_dmft==1.and.me==0) then
 950    call readocc_dmft(paw_dmft,dtfil%filnam_ds(3),dtfil%filnam_ds(4))
 951  end if
 952  !Should be done inside init_sc_dmft
 953  if ( dtset%usedmft /= 0 ) then
 954    call data4entropyDMFT_init(paw_dmft%forentropyDMFT,&
 955    dtset%natom,&
 956    dtset%typat,&
 957    dtset%lpawu,&
 958    dtset%dmft_t2g==1, &
 959    args_gs%upawu,&
 960    args_gs%jpawu)
 961  end if
 962 !write(std_out,*) "paw_dmft%use_dmft",paw_dmft%use_dmft
 963 
 964 !PAW: 1- Initialize values for several arrays unchanged during iterations
 965 !2- Initialize data for LDA+U
 966 !3- Eventually open temporary storage file
 967  if(psps%usepaw==1) then
 968 !  1-
 969    gnt_option=1;if (dtset%pawxcdev==2.or.(dtset%pawxcdev==1.and.dtset%positron/=0)) gnt_option=2
 970 
 971 !  Test if we have to call pawinit
 972 !  Some gen-cond have to be added...
 973    call paw_gencond(Dtset,gnt_option,"test",call_pawinit)
 974 
 975    if (psp_gencond==1.or.call_pawinit) then
 976      call timab(553,1,tsec)
 977      gsqcut_shp=two*abs(dtset%diecut)*dtset%dilatmx**2/pi**2
 978      hyb_range_fock=zero;if (dtset%ixc<0) call libxc_functionals_get_hybridparams(hyb_range=hyb_range_fock)
 979      call pawinit(gnt_option,gsqcut_shp,hyb_range_fock,dtset%pawlcutd,dtset%pawlmix,&
 980 &     psps%mpsang,dtset%pawnphi,dtset%nsym,dtset%pawntheta,&
 981 &     pawang,pawrad,dtset%pawspnorb,pawtab,dtset%pawxcdev,dtset%xclevel,dtset%usepotzero)
 982 
 983      ! Update internal values
 984      call paw_gencond(Dtset,gnt_option,"save",call_pawinit)
 985      call timab(553,2,tsec)
 986 #if defined HAVE_BIGDFT
 987 !    In the PAW+WVL case, copy sij:
 988      if(dtset%usewvl==1) then
 989        do itypat=1,dtset%ntypat
 990          wvl%descr%paw%sij(1:pawtab(itypat)%lmn2_size,itypat)=pawtab(itypat)%sij(:)
 991        end do
 992      end if
 993 #endif
 994    end if
 995    psps%n1xccc=maxval(pawtab(1:psps%ntypat)%usetcore)
 996    call setsym_ylm(gprimd,pawang%l_max-1,dtset%nsym,dtset%pawprtvol,rprimd,symrec,pawang%zarot)
 997 
 998 !  2-Initialize and compute data for LDA+U, EXX, or LDA+DMFT
 999    pawtab(:)%usepawu=0
1000    pawtab(:)%useexexch=0
1001    pawtab(:)%exchmix=zero
1002    if(paw_dmft%use_dmft==1) then
1003      call print_sc_dmft(paw_dmft,dtset%pawprtvol)
1004    end if
1005    if (dtset%usepawu>0.or.dtset%useexexch>0.or.paw_dmft%use_dmft>0) then
1006      call pawpuxinit(dtset%dmatpuopt,dtset%exchmix,dtset%f4of2_sla,dtset%f6of2_sla,&
1007 &     args_gs%jpawu,dtset%lexexch,dtset%lpawu,dtset%ntypat,pawang,dtset%pawprtvol,&
1008 &     pawrad,pawtab,args_gs%upawu,dtset%usedmft,dtset%useexexch,dtset%usepawu,ucrpa=dtset%ucrpa)
1009    end if
1010  end if
1011 
1012 !###########################################################
1013 !### 11. Initialize (eventually) electron-positron data and
1014 !###     electric and magnetic field data
1015 
1016 !Initialize (eventually) electron-positron data
1017  nullify (electronpositron)
1018  if (dtset%positron/=0) then
1019    call init_electronpositron(dtfil%ireadwf,dtset,electronpositron,mpi_enreg,nfftf,pawrhoij,pawtab)
1020  end if
1021 
1022 !###########################################################
1023 ! Initialisation of cprj
1024  usecprj=0; mcprj=0;mband_cprj=0
1025  if (dtset%usepaw==1) then
1026    if (associated(electronpositron)) then
1027      if (dtset%positron/=0.and.electronpositron%dimcprj>0) usecprj=1
1028    end if
1029    if (dtset%prtnabla>0) usecprj=1
1030    if (dtset%extrapwf>0) usecprj=1
1031    if (dtset%pawfatbnd>0)usecprj=1
1032    if (dtset%prtdos==3)  usecprj=1
1033    if (dtset%usewvl==1)  usecprj=1
1034    if (dtset%nstep==0) usecprj=0
1035    if (dtset%usefock==1)  usecprj=1
1036  end if
1037  if (usecprj==0) then
1038    ABI_DATATYPE_ALLOCATE(cprj,(0,0))
1039  end if
1040  if (usecprj==1) then
1041    mband_cprj=dtset%mband;if (dtset%paral_kgb/=0) mband_cprj=mband_cprj/mpi_enreg%nproc_band
1042    mcprj=my_nspinor*mband_cprj*dtset%mkmem*dtset%nsppol
1043 !Was allocated above for valgrind sake so should always be true (safety)
1044    if (allocated(cprj)) then
1045      call pawcprj_free(cprj)
1046      ABI_DATATYPE_DEALLOCATE(cprj)
1047    end if
1048    ABI_DATATYPE_ALLOCATE(cprj,(dtset%natom,mcprj))
1049    ncpgr=0
1050    if (dtset%usefock==1) then
1051      if (dtset%optforces == 1) then
1052        ncpgr = 3
1053      end if
1054 !       if (dtset%optstress /= 0) then
1055 !         ncpgr = 6 ; ctocprj_choice = 3
1056 !       end if
1057    end if
1058    ABI_ALLOCATE(dimcprj_srt,(dtset%natom))
1059    call pawcprj_getdim(dimcprj_srt,dtset%natom,nattyp,dtset%ntypat,dtset%typat,pawtab,'O')
1060    call pawcprj_alloc(cprj,ncpgr,dimcprj_srt)
1061  end if
1062 
1063 
1064 !###########################################################
1065 !### 12. Operations dependent of iscf value
1066 
1067 !Get starting charge density : rhor as well as rhog
1068 !Also initialize the kinetic energy density
1069  if (scf_history%history_size>=0) then
1070    rhor => scf_history%rhor_last
1071    taur => scf_history%taur_last
1072  else
1073    ABI_ALLOCATE(rhor,(nfftf,dtset%nspden))
1074    ABI_ALLOCATE(taur,(nfftf,dtset%nspden*dtset%usekden))
1075  end if
1076  ABI_ALLOCATE(rhog,(2,nfftf))
1077  ABI_ALLOCATE(taug,(2,nfftf*dtset%usekden))
1078 
1079  if (has_to_init) then
1080    if (dtset%iscf>0 .or. (dtset%iscf==0 .and. dtset%usewvl==1 )) then ! .and. dtset%usepaw==1)) then
1081 
1082      if(dtfil%ireadden/=0.and.dtset%positron<=0)then
1083        ! Read density
1084        rdwrpaw=psps%usepaw; if(dtfil%ireadwf/=0) rdwrpaw=0
1085        if (dtset%usewvl==0) then
1086          call read_rhor(dtfil%fildensin, cplex1, dtset%nspden, nfftf, ngfftf, rdwrpaw, &
1087          mpi_enreg, rhor, hdr_den, pawrhoij, comm, check_hdr=hdr, allow_interp=.True.)
1088          results_gs%etotal = hdr_den%etot; call hdr_free(hdr_den)
1089        else
1090          fform=52 ; accessfil=0
1091          if (dtset%iomode == IO_MODE_MPI ) accessfil=4
1092          if (dtset%iomode == IO_MODE_ETSF) accessfil=3
1093          call ioarr(accessfil,rhor,dtset,results_gs%etotal,fform,dtfil%fildensin,hdr,&
1094 &         mpi_enreg,ngfftf,cplex1,nfftf,pawrhoij,1,rdwrpaw,wvl%den)
1095        end if
1096 
1097        if (rdwrpaw/=0) then
1098          call hdr_update(hdr,bantot,etot,fermie,residm,&
1099 &         rprimd,occ,pawrhoij,xred,args_gs%amu,&
1100 &         comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
1101        end if
1102 
1103        ! Read kinetic energy density
1104        if(dtfil%ireadkden/=0 .and. dtset%usekden==1 )then
1105          call read_rhor(dtfil%filkdensin, cplex1, dtset%nspden, nfftf, ngfftf, rdwrpaw, &
1106          mpi_enreg, taur, hdr_den, pawrhoij, comm, check_hdr=hdr, allow_interp=.True.)
1107          call hdr_free(hdr_den)
1108        end if
1109 
1110 !      Compute up+down rho(G) by fft
1111        ABI_ALLOCATE(work,(nfftf))
1112        work(:)=rhor(:,1)
1113        call fourdp(1,rhog,work,-1,mpi_enreg,nfftf,ngfftf,dtset%paral_kgb,0)
1114        if(dtset%usekden==1)then
1115          work(:)=taur(:,1)
1116          call fourdp(1,taug,work,-1,mpi_enreg,nfftf,ngfftf,dtset%paral_kgb,0)
1117        end if
1118        ABI_DEALLOCATE(work)
1119 
1120      else if(dtfil%ireadwf/=0)then
1121        izero=0
1122 !      Obtain the charge density from wfs that were read previously
1123 !      Be careful: in PAW, rho does not include the compensation
1124 !      density (to be added in scfcv.F90) !
1125 !      tim_mkrho=1 ; mpi_enreg%paralbd=0
1126        tim_mkrho=1
1127        if (psps%usepaw==1) then
1128          ABI_ALLOCATE(rhowfg,(2,dtset%nfft))
1129          ABI_ALLOCATE(rhowfr,(dtset%nfft,dtset%nspden))
1130 !        write(std_out,*) "mkrhogstate"
1131          call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,&
1132 &         mpi_enreg,npwarr,occ,paw_dmft,phnons,rhowfg,rhowfr,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs)
1133          call transgrid(1,mpi_enreg,dtset%nspden,+1,1,1,dtset%paral_kgb,pawfgr,rhowfg,rhog,rhowfr,rhor)
1134          ABI_DEALLOCATE(rhowfg)
1135          ABI_DEALLOCATE(rhowfr)
1136        else
1137          call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,&
1138 &         mpi_enreg,npwarr,occ,paw_dmft,phnons,rhog,rhor,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs)
1139          if(dtset%usekden==1)then
1140            call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,&
1141 &           mpi_enreg,npwarr,occ,paw_dmft,phnons,taug,taur,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs,option=1)
1142          end if
1143 
1144        end if
1145 
1146      else if(dtfil%ireadwf==0.and.dtset%positron/=1)then
1147 
1148 !      Crude, but realistic initialisation of the density
1149 !      There is not point to compute it from random wavefunctions.
1150        if (dtset%usewvl == 0) then
1151          call initro(atindx,dtset%densty,gmet,gsqcut_eff,psps%usepaw,&
1152 &         mgfftf,mpi_enreg,psps%mqgrid_vl,dtset%natom,nattyp,nfftf,&
1153 &         ngfftf,dtset%nspden,psps%ntypat,dtset%paral_kgb,psps,pawtab,ph1df,&
1154 &         psps%qgrid_vl,rhog,rhor,dtset%spinat,ucvol,psps%usepaw,&
1155 &         dtset%ziontypat,dtset%znucl)
1156 !        Update initialized density taking into account jellium slab
1157          if(dtset%jellslab/=0) then
1158            option=2
1159            ABI_ALLOCATE(work,(nfftf))
1160            call jellium(gmet,gsqcut_eff,mpi_enreg,nfftf,ngfftf,dtset%nspden,&
1161 &           option,dtset%paral_kgb,dtset%slabwsrad,rhog,rhor,rprimd,work,dtset%slabzbeg,dtset%slabzend)
1162            ABI_DEALLOCATE(work)
1163          end if ! of usejell
1164 !        Kinetic energy density initialized to zero (used only in metaGGAs ... )
1165          if(dtset%usekden==1)then
1166            taur=zero ; taug=zero
1167          end if
1168        else if (dtset%usewvl/=0) then
1169 !        skip for the moment for wvl+paw, not yet coded
1170          if(dtset%usepaw==0) then
1171            !Wavelet density corresponds exactly to the wavefunctions,
1172            !since wavefunctions are taken from diagonalisation of LCAO.
1173            call wvl_mkrho(dtset, irrzon, mpi_enreg, phnons, rhor, wvl%wfs, wvl%den)
1174          else !usepaw
1175 #if defined HAVE_BIGDFT
1176            call wvl_initro(atindx1,wvl%descr%atoms%astruct%geocode,wvl%descr%h,mpi_enreg%me_wvl,&
1177 &           dtset%natom,nattyp,nfftf,dtset%nspden,psps%ntypat,&
1178 &           wvl%descr%Glr%d%n1,wvl%descr%Glr%d%n1i,&
1179 &           wvl%descr%Glr%d%n2,wvl%descr%Glr%d%n2i,&
1180 &           wvl%descr%Glr%d%n3,&
1181 &           pawrad,pawtab,psps%gth_params%psppar,rhor,rprimd,&
1182 &           dtset%spinat,wvl%den,dtset%xc_denpos,xred,dtset%ziontypat)
1183            call wvl_mkrho(dtset, irrzon, mpi_enreg, phnons, rhor, wvl%wfs, wvl%den)
1184 #endif
1185          end if
1186        end if
1187 
1188      end if
1189 
1190    else if ((dtset%iscf==-1.or.dtset%iscf==-2.or.dtset%iscf==-3).and.dtset%positron<=0) then
1191 
1192      call status(0,dtfil%filstat,iexit,level,'call ioarr    ')
1193 !    Read rho(r) from a disk file
1194      rdwrpaw=psps%usepaw
1195 !    Note : results_gs%etotal is read here,
1196 !    and might serve in the tddft routine, but it is contrary to the intended use of results_gs ...
1197 !    Warning : should check the use of results_gs%e_fermie
1198 !    Warning : should check the use of results_gs%residm
1199 !    One might make them separate variables.
1200 
1201     ! Read density and get Fermi level from hdr_den
1202      call read_rhor(dtfil%fildensin, cplex1, dtset%nspden, nfftf, ngfftf, rdwrpaw, &
1203      mpi_enreg, rhor, hdr_den, pawrhoij, comm, check_hdr=hdr)
1204      results_gs%etotal = hdr_den%etot; results_gs%energies%e_fermie = hdr_den%fermie
1205      call hdr_free(hdr_den)
1206 
1207      if(dtfil%ireadkden/=0 .and. dtset%usekden==1)then
1208        call read_rhor(dtfil%filkdensin, cplex1, dtset%nspden, nfftf, ngfftf, rdwrpaw, &
1209        mpi_enreg, taur, hdr_den, pawrhoij, comm, check_hdr=hdr)
1210        call hdr_free(hdr_den)
1211      end if
1212 
1213 !    Compute up+down rho(G) by fft
1214      ABI_ALLOCATE(work,(nfftf))
1215      work(:)=rhor(:,1)
1216      call fourdp(1,rhog,work,-1,mpi_enreg,nfftf,ngfftf,dtset%paral_kgb,0)
1217      if(dtset%usekden==1)then
1218        work(:)=taur(:,1)
1219        call fourdp(1,taug,work,-1,mpi_enreg,nfftf,ngfftf,dtset%paral_kgb,0)
1220      end if
1221      ABI_DEALLOCATE(work)
1222 
1223    end if
1224  end if ! has_to_init
1225 
1226 !###########################################################
1227 !### 13. If needed, initialize SCF history variables
1228 
1229 !If needed, initialize atomic density in SCF history
1230  if (scf_history%history_size>0.and.has_to_init) then
1231 !  If rhor is an atomic density, just store it in history
1232    if (.not.read_wf_or_den) then
1233      scf_history%atmrho_last(:)=rhor(:,1)
1234    else
1235 !    If rhor is not an atomic density, has to compute rho_at(r)
1236      ABI_ALLOCATE(rhowfg,(2,nfftf))
1237      ABI_ALLOCATE(rhowfr,(nfftf,1))
1238      ABI_ALLOCATE(spinat_dum,(3,dtset%natom))
1239      spinat_dum=zero
1240      call initro(atindx,dtset%densty,gmet,gsqcut_eff,psps%usepaw,mgfftf,mpi_enreg,&
1241 &     psps%mqgrid_vl,dtset%natom,nattyp,nfftf,ngfftf,1,psps%ntypat,dtset%paral_kgb,psps,pawtab,&
1242 &     ph1df,psps%qgrid_vl,rhowfg,rhowfr,spinat_dum,ucvol,&
1243 &     psps%usepaw,dtset%ziontypat,dtset%znucl)
1244      scf_history%atmrho_last(:)=rhowfr(:,1)
1245      ABI_DEALLOCATE(rhowfg)
1246      ABI_DEALLOCATE(rhowfr)
1247      ABI_DEALLOCATE(spinat_dum)
1248    end if
1249  end if
1250 
1251  if ((.not.read_wf_or_den).or.(scf_history%history_size>0.and.has_to_init))  then
1252    ABI_DEALLOCATE(ph1df)
1253  end if
1254 
1255 !!Electric field: initialization stage
1256 !!further initialization and updates happen in scfcv.F90
1257  call init_e_field_vars(dtefield,dtset,gmet,gprimd,kg,&
1258 & mpi_enreg,npwarr,occ,pawang,pawrad,pawtab,psps,&
1259 & pwind,pwind_alloc,pwnsfac,rprimd,symrec,xred)
1260 
1261  !! orbital magnetization initialization
1262  dtorbmag%orbmag = dtset%orbmag
1263  if (dtorbmag%orbmag > 0) then
1264     call initorbmag(dtorbmag,dtset,gmet,gprimd,kg,mpi_enreg,npwarr,occ,&
1265 &                   pawtab,psps,pwind,pwind_alloc,pwnsfac,&
1266 &                   rprimd,symrec,xred)
1267  end if
1268 
1269 
1270  fatvshift=one
1271 
1272  if (dtset%usewvl == 1 .and. wvl_debug) then
1273 #if defined HAVE_BIGDFT
1274    call wvl_timing(me,'INIT','PR')
1275 #endif
1276  end if
1277 
1278 !Check whether exiting was required by the user. If found then do not start minimization steps
1279 !At this first call to chkexi, initialize cpus, if it
1280 !is non-zero (which would mean that no action has to be taken)
1281 !Should do this in driver ...
1282  cpus=dtset%cpus
1283  if(abs(cpus)>1.0d-5)cpus=cpus+cpui
1284  openexit=1 ; if(dtset%chkexit==0) openexit=0
1285  call exit_check(cpus,dtfil%filnam_ds(1),iexit,ab_out,mpi_enreg%comm_cell,openexit)
1286 
1287 !If immediate exit, and wavefunctions were not read, must zero eigenvalues
1288  if (iexit/=0) eigen(:)=zero
1289 
1290  call timab(34,2,tsec)
1291 
1292  conv_retcode = 0
1293 
1294  if (iexit==0) then
1295 
1296 !  ###########################################################
1297 !  ### 14. Move atoms and acell according to ionmov value
1298 
1299 
1300 !  Eventually symmetrize atomic coordinates over space group elements:
1301 !  call symmetrize_xred(indsym,dtset%natom,dtset%nsym,dtset%symrel,dtset%tnons,xred)
1302 
1303 !  If atoms are not being moved and U should not be determined,
1304 !  use scfcv directly; else
1305 !  call move, pawuj_drive or brdmin which in turn calls scfcv.
1306 
1307    call timab(35,3,tsec)
1308 
1309    call scfcv_init(scfcv_args,atindx,atindx1,cg,cprj,cpus,&
1310 &   args_gs%dmatpawu,dtefield,dtfil,dtorbmag,dtpawuj,dtset,ecore,eigen,hdr,&
1311 &   indsym,initialized,irrzon,kg,mcg,mcprj,mpi_enreg,my_natom,nattyp,ndtpawuj,&
1312 &   nfftf,npwarr,occ,pawang,pawfgr,pawrad,pawrhoij,&
1313 &   pawtab,phnons,psps,pwind,pwind_alloc,pwnsfac,rec_set,&
1314 &   resid,results_gs,scf_history,fatvshift,&
1315 &   symrec,taug,taur,wvl,ylm,ylmgr,paw_dmft,wffnew,wffnow)
1316 
1317    call dtfil_init_time(dtfil,0)
1318 
1319    write(message,'(a,80a)')ch10,('=',mu=1,80)
1320    call wrtout(ab_out,message,'COLL')
1321    call wrtout(std_out,message,'COLL')
1322 
1323    if (dtset%ionmov==0 .or. dtset%imgmov==6) then
1324 
1325 !    Should merge this call with the call for dtset%ionmov==4 and 5
1326 
1327      if (dtset%macro_uj==0) then
1328 
1329        call status(0,dtfil%filstat,iexit,level,'call scfcv_run')
1330        call scfcv_run(scfcv_args,electronpositron,rhog,rhor,rprimd,xred,xred_old,conv_retcode)
1331 
1332      else
1333 !      Conduct determination of U
1334 
1335        call pawuj_drive(scfcv_args,dtset,electronpositron,rhog,rhor,rprimd,xred,xred_old)
1336      end if
1337 
1338 !    ========================================
1339 !    New structure for geometry optimization
1340 !    ========================================
1341    else if (dtset%ionmov>50.or.dtset%ionmov<=27) then
1342 
1343      ! TODO: return conv_retcode
1344      call mover(scfcv_args,ab_xfh,acell,args_gs%amu,dtfil,&
1345 &     electronpositron,rhog,rhor,rprimd,vel,vel_cell,xred,xred_old)
1346 
1347 !    Compute rprim from rprimd and acell
1348      do kk=1,3
1349        do jj=1,3
1350          rprim(jj,kk)=rprimd(jj,kk)/acell(kk)
1351        end do
1352      end do
1353 
1354 !    =========================================
1355 !    New structure for geometry optimization
1356 !    =========================================
1357 
1358    else ! Not an allowed option
1359      write(message, '(a,i0,2a)' )&
1360 &     'Disallowed value for ionmov=',dtset%ionmov,ch10,&
1361 &     'Allowed values are: 1,2,3,4,5,6,7,8,9,10,11,12,13,14,20,21,22,23,24 and 30'
1362      MSG_BUG(message)
1363    end if
1364 
1365    call scfcv_destroy(scfcv_args)
1366 
1367    call timab(35,2,tsec)
1368 
1369 !  ###########################################################
1370 !  ### 15. Final operations and output for gstate
1371 
1372  end if !  End of the check of hasty exit
1373 
1374  call timab(36,3,tsec)
1375 
1376  write(message, '(80a,a,a,a,a)' ) ('=',mu=1,80),ch10,ch10,&
1377 & ' ----iterations are completed or convergence reached----',ch10
1378  call wrtout(ab_out,message,'COLL')
1379  call wrtout(std_out,  message,'COLL')
1380 
1381 !Mark this GS computation as done
1382  initialized=1
1383  if (dtset%usewvl == 1 .and. wvl_debug) then
1384 #if defined HAVE_BIGDFT
1385    call wvl_timing(me,'WFN_OPT','PR')
1386 #endif
1387  end if
1388 
1389 !Will be put here later.
1390 !! ! WVL - maybe compute the tail corrections to energy
1391 !! if (dtset%tl_radius > zero) then
1392 !!    ! Store xcart for each atom
1393 !!    allocate(xcart(3, dtset%natom))
1394 !!    call xred2xcart(dtset%natom, rprimd, xcart, xred)
1395 !!    ! Use the tails to improve energy precision.
1396 !!    call wvl_tail_corrections(dtset, results_gs%energies, results_gs%etotal, &
1397 !!         & mpi_enreg, occ, psps, vtrial, wvl, xcart)
1398 !!    deallocate(xcart)
1399 !! end if
1400 
1401 !Update the header, before using it
1402  call hdr_update(hdr,bantot,results_gs%etotal,results_gs%energies%e_fermie,&
1403 & results_gs%residm,rprimd,occ,pawrhoij,xred,args_gs%amu,&
1404 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
1405 
1406  ABI_ALLOCATE(doccde,(dtset%mband*dtset%nkpt*dtset%nsppol))
1407  doccde=zero
1408 
1409  call ebands_init(bantot,ebands,dtset%nelect,doccde,eigen,hdr%istwfk,hdr%kptns,hdr%nband,&
1410 & hdr%nkpt,hdr%npwarr,hdr%nsppol,hdr%nspinor,hdr%tphysel,hdr%tsmear,hdr%occopt,hdr%occ,hdr%wtk,&
1411 & hdr%charge, hdr%kptopt, hdr%kptrlatt_orig, hdr%nshiftk_orig, hdr%shiftk_orig, &
1412 & hdr%kptrlatt, hdr%nshiftk, hdr%shiftk)
1413 
1414  ebands%fermie = results_gs%energies%e_fermie
1415  ABI_DEALLOCATE(doccde)
1416  !write(std_out,*)"efermi after ebands_init",ebands%fermie
1417 
1418  ! Compute and print the gaps.
1419  call ebands_report_gap(ebands,header="Gap info",unit=std_out,mode_paral="COLL",gaps=results_gs%gaps)
1420 
1421  if(dtset%nqpt==0)filnam=dtfil%fnameabo_wfk
1422  if(dtset%nqpt==1)filnam=dtfil%fnameabo_wfq
1423 
1424  ! Write wavefunctions file only if convergence was not achieved.
1425  !write(std_out,*)"conv_retcode", conv_retcode
1426  write_wfk = .True.
1427  if (dtset%prtwf==-1 .and. conv_retcode == 0) then
1428    write_wfk = .False.
1429    message = "GS calculation converged with prtwf=-1 --> Skipping WFK file output"
1430    call wrtout(ab_out, message, "COLL")
1431    MSG_COMMENT(message)
1432  end if
1433 
1434 !To print out the WFs, need the rprimd that was used to generate the G vectors
1435  hdr%rprimd=rprimd_for_kg
1436 
1437  if (write_wfk) then
1438    call outwf(cg,dtset,psps,eigen,filnam,hdr,kg,dtset%kptns,&
1439 &   dtset%mband,mcg,dtset%mkmem,mpi_enreg,dtset%mpw,dtset%natom,&
1440 &   dtset%nband,dtset%nkpt,npwarr,dtset%nsppol,&
1441 &   occ,resid,response,dtfil%unwff2,wvl%wfs,wvl%descr)
1442    !SPr: add input variable managing the .vtk file OUTPUT (Please don't remove the next commented line)
1443    !call printmagvtk(mpi_enreg,cplex1,dtset%nspden,nfftf,ngfftf,rhor,rprimd,'DEN')
1444  end if
1445 
1446  if (dtset%prtwf==2) then
1447    call outqmc(cg,dtset,eigen,gprimd,hdr,kg,mcg,mpi_enreg,npwarr,occ,psps,results_gs)
1448  end if
1449 
1450 !Restore the original rprimd in hdr
1451  hdr%rprimd=rprimd
1452 
1453  ! Generate WFK in full BZ (needed by LOBSTER)
1454  if (me == master .and. dtset%prtwf == 1 .and. dtset%prtwf_full == 1 .and. dtset%nqpt == 0) then
1455    wfkfull_path = strcat(dtfil%filnam_ds(4), "_FULL_WFK")
1456    if (dtset%iomode == IO_MODE_ETSF) wfkfull_path = nctk_ncify(wfkfull_path)
1457    call wfk_tofullbz(filnam, dtset, psps, pawtab, wfkfull_path)
1458 
1459    ! Write tetrahedron tables.
1460    call crystal_from_hdr(cryst, hdr, 2)
1461    tetra = tetra_from_kptrlatt(cryst, dtset%kptopt, dtset%kptrlatt, dtset%nshiftk, &
1462    dtset%shiftk, dtset%nkpt, dtset%kptns, message, ierr)
1463    if (ierr == 0) then
1464      call tetra_write(tetra, dtset%nkpt, dtset%kptns, strcat(dtfil%filnam_ds(4), "_TETRA"))
1465    else
1466      MSG_WARNING(sjoin("Cannot produce TETRA file", ch10, message))
1467    end if
1468 
1469    call destroy_tetra(tetra)
1470    call crystal_free(cryst)
1471  end if
1472 
1473  call clnup1(acell,dtset,eigen,results_gs%energies%e_fermie,&
1474 & dtfil%fnameabo_dos,dtfil%fnameabo_eig,results_gs%fred,&
1475 & mpi_enreg,nfftf,ngfftf,occ,dtset%optforces,&
1476 & resid,rhor,rprimd,results_gs%vxcavg,xred)
1477 
1478  if ( (dtset%iscf>=0 .or. dtset%iscf==-3) .and. dtset%prtstm==0) then
1479    call prtene(dtset,results_gs%energies,ab_out,psps%usepaw)
1480  end if
1481 
1482 !write final electric field components HONG
1483 
1484  if (dtset%berryopt == 4 .or. dtset%berryopt == 6 .or. dtset%berryopt ==7 .or.  &
1485 & dtset%berryopt == 14 .or. dtset%berryopt == 16 .or. dtset%berryopt ==17 ) then ! output final elctric field data    !!HONG
1486    if (dtset%berryopt == 4) then
1487      write(message,'(a,a)')   ch10, 'Constant unreduced E calculation  - final values:'
1488    else if (dtset%berryopt == 6 ) then
1489      write(message,'(a,a)')   ch10, 'Constant unreduced D calculation  - final values:'
1490    else if (dtset%berryopt == 14) then
1491      write(message,'(a,a)')   ch10, 'Constant reduced ebar calculation  - final values:'
1492    else if (dtset%berryopt == 16 ) then
1493      write(message,'(a,a)')   ch10, 'Constant reduced d calculation  - final values:'
1494    else if (dtset%berryopt == 17) then
1495      write(message,'(a,a)')   ch10, 'Constant reduced ebar and d calculation  - final values:'
1496    end if
1497 
1498    call wrtout(ab_out,message,'COLL')
1499    call prtefield(dtset,dtefield,ab_out,rprimd)
1500 
1501    call wrtout(std_out,message,'COLL')
1502    call prtefield(dtset,dtefield,std_out,rprimd)
1503 
1504 !  To check if the final electric field is below the critical field
1505    do kk = 1, 3
1506      efield_band(kk) = abs(dtset%red_efieldbar(kk))*dtefield%nkstr(kk)
1507    end do
1508 !  eg = maxval(eg_dir)
1509 !  eg_ev = eg*Ha_eV
1510    write(message,'(a,a,a,a,a,a,a,a,f7.2,a,a)')ch10,&
1511 &   ' Please check: COMMENT - ',ch10,&
1512 &   '  As a rough estimate,',ch10,&
1513 &   '  to be below the critical field, the bandgap of your system',ch10,&
1514 &   '  should be larger than ',maxval(efield_band)*Ha_eV,' eV.',ch10
1515    call wrtout(ab_out,message,'COLL')
1516    call wrtout(std_out,message,'COLL')
1517 
1518    write(message,'(a)')  '--------------------------------------------------------------------------------'
1519    call wrtout(ab_out,message,'COLL')
1520    call wrtout(std_out,message,'COLL')
1521 
1522  end if
1523 
1524 !Open the formatted derivative database file, and write the preliminary information
1525 !In the // case, only one processor writes the energy and the gradients to the DDB
1526 
1527  if (me==0.and.dtset%nimage==1.and.((dtset%iscf > 0).or.&
1528 & (dtset%berryopt == -1).or.(dtset%berryopt) == -3)) then
1529 
1530    if (dtset%iscf > 0) then
1531      nblok = 2          ! 1st blok = energy, 2nd blok = gradients
1532    else
1533      nblok = 1
1534    end if
1535 
1536    dscrpt=' Note : temporary (transfer) database '
1537    ddbnm=trim(dtfil%filnam_ds(4))//'_DDB'
1538 
1539    call ddb_hdr_init(ddb_hdr,dtset,psps,pawtab,DDB_VERSION,dscrpt,&
1540 &   nblok,xred=xred,occ=occ,ngfft=ngfft)
1541 
1542    call ddb_hdr_open_write(ddb_hdr,ddbnm,dtfil%unddb,fullinit=0)
1543 
1544    call ddb_hdr_free(ddb_hdr)
1545 
1546    choice=2
1547    mpert = dtset%natom + 6 ; msize = 3*mpert
1548 
1549 !  create a ddb structure with just one blok
1550    call ddb_malloc(ddb,msize,1,dtset%natom,dtset%ntypat)
1551 
1552    ddb%flg = 0
1553    ddb%qpt = zero
1554    ddb%nrm = one
1555    ddb%val = zero
1556 
1557 !  Write total energy to the DDB
1558    if (dtset%iscf > 0) then
1559      ddb%typ(1) = 0
1560      ddb%val(1,1,1) = results_gs%etotal
1561      ddb%flg(1,1) = 1
1562      call ddb_write_blok(ddb,1,choice,dtset%mband,mpert,msize,dtset%nkpt,dtfil%unddb)
1563    end if
1564 
1565 !  Write gradients to the DDB
1566    ddb%typ = 4
1567    ddb%flg = 0
1568    ddb%val = zero
1569    indx = 0
1570    if (dtset%iscf > 0) then
1571      do iatom = 1, dtset%natom
1572        do idir = 1, 3
1573          indx = indx + 1
1574          ddb%flg(indx,1) = 1
1575          ddb%val(1,indx,1) = results_gs%fred(idir,iatom)
1576        end do
1577      end do
1578    end if
1579 
1580    indx = 3*dtset%natom + 3
1581    if ((abs(dtset%berryopt) == 1).or.(abs(dtset%berryopt) == 3)) then
1582      do idir = 1, 3
1583        indx = indx + 1
1584        if (dtset%rfdir(idir) == 1) then
1585          ddb%flg(indx,1) = 1
1586          ddb%val(1,indx,1) = results_gs%pel(idir)
1587        end if
1588      end do
1589    end if
1590 
1591    indx = 3*dtset%natom + 6
1592    if (dtset%iscf > 0) then
1593      ddb%flg(indx+1:indx+6,1) = 1
1594      ddb%val(1,indx+1:indx+6,1) = results_gs%strten(1:6)
1595    end if
1596 
1597    call ddb_write_blok(ddb,1,choice,dtset%mband,mpert,msize,dtset%nkpt,dtfil%unddb)
1598 
1599    call ddb_free(ddb)
1600 
1601 !  Close DDB
1602    close(dtfil%unddb)
1603  end if
1604 
1605  if (dtset%nstep>0 .and. dtset%prtstm==0 .and. dtset%positron/=1) then
1606    call clnup2(psps%n1xccc,results_gs%fred,results_gs%grchempottn,results_gs%gresid,&
1607 &   results_gs%grewtn,results_gs%grvdw,results_gs%grxc,dtset%iscf,dtset%natom,&
1608 &   results_gs%ngrvdw,dtset%optforces,dtset%optstress,dtset%prtvol,start,&
1609 &   results_gs%strten,results_gs%synlgr,xred)
1610  end if
1611 
1612  if(dtset%imgwfstor==1)then
1613    scf_history%cg(:,:,1)=cg(:,:)
1614    scf_history%eigen(:,1)=eigen(:)
1615  endif
1616 
1617 !Deallocate arrays
1618  ABI_DEALLOCATE(atindx)
1619  ABI_DEALLOCATE(atindx1)
1620  ABI_DEALLOCATE(indsym)
1621  ABI_DEALLOCATE(npwarr)
1622  ABI_DEALLOCATE(nattyp)
1623  ABI_DEALLOCATE(resid)
1624  ABI_DEALLOCATE(rhog)
1625  ABI_DEALLOCATE(start)
1626  ABI_DEALLOCATE(symrec)
1627  ABI_DEALLOCATE(taug)
1628  ABI_DEALLOCATE(ab_xfh%xfhist)
1629  call pawfgr_destroy(pawfgr)
1630 
1631  if(dtset%imgwfstor==0)then
1632    ABI_DEALLOCATE(cg)
1633    ABI_DEALLOCATE(eigen)
1634  else
1635    nullify(cg,eigen)
1636  endif
1637 
1638  if (dtset%usewvl == 0 .or. dtset%nsym <= 1) then
1639 !  In wavelet case, irrzon and phnons are deallocated by wavelet object.
1640    ABI_DEALLOCATE(irrzon)
1641    ABI_DEALLOCATE(phnons)
1642  end if
1643 
1644  ABI_DEALLOCATE(ylm)
1645  ABI_DEALLOCATE(ylmgr)
1646 
1647  if (scf_history%history_size<0) then
1648    if (psps%usepaw==1) then
1649      call pawrhoij_free(pawrhoij)
1650    end if
1651    ABI_DEALLOCATE(rhor)
1652    ABI_DEALLOCATE(taur)
1653    ABI_DATATYPE_DEALLOCATE(pawrhoij)
1654    ABI_DEALLOCATE(xred_old)
1655  else
1656    nullify(rhor,taur,pawrhoij,xred_old)
1657  end if
1658 
1659 !PAW+DMFT
1660  call destroy_sc_dmft(paw_dmft)
1661  ! This call should be done inside destroy_sc_dmft
1662  if ( dtset%usedmft /= 0 ) then
1663    call data4entropyDMFT_destroy(paw_dmft%forentropyDMFT)
1664  end if
1665 
1666 !Destroy electronpositron datastructure
1667  if (dtset%positron/=0) then
1668    call destroy_electronpositron(electronpositron)
1669  end if
1670 
1671 !Deallocating the basis set.
1672  if (dtset%usewvl == 1) then
1673    call wvl_projectors_free(wvl%projectors)
1674    call wvl_wfs_free(wvl%wfs)
1675    call wvl_descr_free(wvl%descr)
1676    call wvl_denspot_free(wvl%den)
1677    if(dtset%usepaw == 1) then
1678      call wvl_paw_free(wvl%descr)
1679    end if
1680  end if
1681 
1682  ABI_DEALLOCATE(kg)
1683 
1684  if (dtset%icoulomb /= 0) then
1685    call psolver_kernel((/ 0._dp, 0._dp, 0._dp /), 0, dtset%icoulomb, 0, kernel_dummy, &
1686 &   0, dtset%ngfft, 1, dtset%nscforder)
1687  end if
1688 
1689  if (associated(pwind)) then
1690    ABI_DEALLOCATE(pwind)
1691  end if
1692  if (associated(pwnsfac)) then
1693    ABI_DEALLOCATE(pwnsfac)
1694  end if
1695  if ((dtset%berryopt<0).or.&
1696 & (dtset%berryopt== 4.or.dtset%berryopt== 6.or.dtset%berryopt== 7.or.&
1697 & dtset%berryopt==14.or.dtset%berryopt==16.or.dtset%berryopt==17)) then
1698    if (xmpi_paral == 1) then
1699      ABI_DEALLOCATE(mpi_enreg%kptdstrb)
1700      if (dtset%berryopt== 4.or.dtset%berryopt== 6.or.dtset%berryopt== 7.or.&
1701 &     dtset%berryopt==14.or.dtset%berryopt==16.or.dtset%berryopt==17) then
1702        ABI_DEALLOCATE(mpi_enreg%kpt_loc2ibz_sp)
1703      end if
1704    end if
1705    if (allocated(mpi_enreg%kpt_loc2ibz_sp))  then
1706      ABI_DEALLOCATE(mpi_enreg%kpt_loc2ibz_sp)
1707    end if
1708    if (allocated(mpi_enreg%kpt_loc2fbz_sp)) then
1709      ABI_DEALLOCATE(mpi_enreg%kpt_loc2fbz_sp)
1710    end if
1711    if (allocated(mpi_enreg%mkmem)) then
1712      ABI_DEALLOCATE(mpi_enreg%mkmem)
1713    end if
1714  end if
1715  ! deallocate cprj
1716  if(usecprj==1) then
1717    ABI_DEALLOCATE(dimcprj_srt)
1718    call pawcprj_free(cprj)
1719  end if
1720  ABI_DATATYPE_DEALLOCATE(cprj)
1721 
1722  ! deallocate efield
1723  call destroy_efield(dtefield)
1724 
1725  ! deallocate orbmag
1726  call destroy_orbmag(dtorbmag)
1727 
1728 !deallocate Recursion
1729  if (dtset%userec == 1) then
1730    call CleanRec(rec_set)
1731  end if
1732 
1733 !Clean the header
1734  call hdr_free(hdr)
1735  call ebands_free(ebands)
1736 
1737  if (me == master .and. dtset%prtxml == 1) then
1738 !  The dataset given in argument has been treated, then we output its variables.
1739 !  call outvarsXML()
1740 !  gstate() will handle a dataset, so we output the dataSet markup.
1741    write(ab_xml_out, "(A)") '  </dataSet>'
1742  end if
1743 
1744  if (dtset%usewvl == 0 .and. dtset%tfkinfunc /= 2 .and. dtset%optdriver /= RUNL_GWLS) then
1745 !  Plane-wave case
1746    call bandfft_kpt_destroy_array(bandfft_kpt,mpi_enreg)
1747  end if
1748 
1749  if(dtset%wfoptalg == 1 .and. psps%usepaw == 1) then
1750    call destroy_invovl(dtset%nkpt)
1751  end if
1752 
1753  if(gemm_nonlop_use_gemm) then
1754    call destroy_gemm_nonlop(dtset%nkpt)
1755    gemm_nonlop_use_gemm = .false.
1756  end if
1757 
1758 !Eventually clean cuda runtime
1759 #if defined HAVE_GPU_CUDA
1760  if (dtset%use_gpu_cuda==1) then
1761    call dealloc_hamilt_gpu(2,dtset%use_gpu_cuda)
1762  end if
1763 #endif
1764 
1765  call status(0,dtfil%filstat,iexit,level,'exit')
1766 
1767  call timab(36,2,tsec)
1768  call timab(32,2,tsec)
1769 
1770  DBG_EXIT("COLL")
1771 end subroutine gstate

m_gstate/pawuj_drive [ Functions ]

[ Top ] [ m_gstate ] [ Functions ]

NAME

 pawuj_drive

FUNCTION

  Drive for automatic determination of U
  Relevant only in PAW+U context

INPUTS

  atindx(natom)=index table for atoms (see gstate.f)
  atindx1(natom)=index table for atoms, inverse of atindx (see gstate.f)
  cpus= cpu time limit in seconds
  dtfil <type(datafiles_type)>=variables related to files
  dtset <type(dataset_type)>=all input variables for this dataset
   | mband=maximum number of bands
   | mgfft=maximum size of 1D FFTs for the "coarse" grid (see NOTES below)
   | mkmem =number of k points treated by this node.
   | mpw=maximum dimensioned size of npw.
   | natom=number of atoms in cell.
   | nfft=(effective) number of FFT grid points (for this processor)
   |      for the "coarse" grid (see NOTES below)
   | nkpt=number of k points
   | nspden=number of spin-density components
   | nsppol=1 for unpolarized, 2 for spin-polarized
   | nsym=number of symmetry elements in space group
  ecore=core psp energy (part of total energy) (hartree)
  kg(3,mpw*mkmem)=reduced planewave coordinates.
  mpi_enreg=informations about MPI parallelization
  nattyp(ntypat)= # atoms of each type.
  npwarr(nkpt)=number of planewaves in basis at this k point
  nspinor=number of spinorial components of the wavefunctions
  psps <type(pseudopotential_type)>=variables related to pseudopotentials
   | mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
  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)
  ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k point
  ylmgr(mpw*mkmem,3,mpsang*mpsang*useylm)= gradients of real spherical harmonics

OUTPUT

  resid(mband*nkpt*nsppol)=residuals for each band over all k points and spins

SIDE EFFECTS

  cg(2,mpw*nspinor*mband*mkmem*nsppol)=updated wavefunctions.
  dtefield <type(efield_type)> = variables related to Berry phase calculations (see initberry.f)
  eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree)
  electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation
  indsym(4,nsym,natom)=indirect indexing array for atom labels
  initialized= if 0 the initialization of the gstate run is not yet finished
  irrzon(nfft**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data
  nfftf=(effective) number of FFT grid points (for this processor)
       for the "fine" grid (see NOTES below)
  occ(mband*nkpt*nsppol)=occupation number for each band (often 2) at each k point
  pawrhoij(natom*usepaw) <type(pawrhoij_type)>= -PAW only- atomic occupancies
  phnons(2,nfft**(1-1/nsym),(nspden/nsppol)-3*(nspden/4))=nonsymmorphic translation phases
  results_gs <type(results_gs_type)>=results (energy and its components,
   forces and its components, the stress tensor) of a ground-state computation
   (should be made a pure output quantity)
  rhog(2,nfftf)=array for Fourier transform of electron density
  rhor(nfftf,nspden)=array for electron density in el./bohr**3
  rprimd(3,3)=dimensional primitive translations in real space (bohr)
  scf_history <type(scf_history_type)>=arrays obtained from previous SCF cycles
  symrec(3,3,nsym)=symmetry operations in reciprocal space
  taug(2,nfftf*dtset%usekden)=array for Fourier transform of kinetic energy density
  taur(nfftf,nspden*dtset%usekden)=array for kinetic energy density
  wffnew,wffnow=struct info for wf disk files.
  wvl <type(wvl_data)>=all wavelets data.
  xred(3,natom)=reduced dimensionless atomic coordinates
  xred_old(3,natom)= at input, previous reduced dimensionless atomic coordinates
                     at output, current xred is transferred to xred_old

PARENTS

      gstate

CHILDREN

      pawuj_det,pawuj_free,pawuj_ini,scfcv_run

SOURCE

2590 subroutine pawuj_drive(scfcv_args, dtset,electronpositron,rhog,rhor,rprimd, xred,xred_old)
2591 
2592 
2593 !This section has been created automatically by the script Abilint (TD).
2594 !Do not modify the following lines by hand.
2595 #undef ABI_FUNC
2596 #define ABI_FUNC 'pawuj_drive'
2597 !End of the abilint section
2598 
2599  implicit none
2600 
2601 !Arguments ------------------------------------
2602 !scalars
2603  type(scfcv_t), intent(inout) :: scfcv_args
2604  type(dataset_type),intent(inout) :: dtset
2605  type(electronpositron_type),pointer :: electronpositron
2606  !type(wffile_type),intent(inout) :: wffnew,wffnow
2607 !arrays
2608  real(dp), intent(inout) :: rprimd(3,3)
2609  real(dp), pointer :: rhog(:,:),rhor(:,:)
2610  real(dp), intent(inout) :: xred(3,dtset%natom),xred_old(3,dtset%natom)
2611 
2612 !Local variables -------------------------
2613 !scalars
2614  integer,target :: ndtpawuj=4
2615  integer :: iuj,conv_retcode
2616  real(dp) :: ures
2617  !character(len=500) :: message
2618 !arrays
2619  real(dp),allocatable :: cgstart(:,:)
2620  type(macro_uj_type),allocatable,target :: dtpawuj(:)
2621 ! *********************************************************************
2622 
2623  DBG_ENTER("COLL")
2624 
2625  if (dtset%macro_uj==0) then
2626    MSG_BUG('Macro_uj must be set !')
2627  end if
2628 
2629  ABI_DATATYPE_ALLOCATE(dtpawuj,(0:ndtpawuj))
2630  ABI_ALLOCATE(cgstart,(2,scfcv_args%mcg))
2631 
2632 !DEBUG
2633 !write(std_out,*)'pawuj_drive: before ini dtpawuj(:)%iuj ', dtpawuj(:)%iuj
2634 !END DEBUG
2635  call pawuj_ini(dtpawuj,ndtpawuj)
2636 
2637  cgstart=scfcv_args%cg
2638  do iuj=1,ndtpawuj
2639 !  allocate(dtpawuj(iuj)%rprimd(3,3)) ! this has already been done in pawuj_ini
2640    dtpawuj(iuj)%macro_uj=dtset%macro_uj
2641    dtpawuj(iuj)%pawprtvol=dtset%pawprtvol
2642    dtpawuj(iuj)%diemix=dtset%diemix
2643    dtpawuj(iuj)%pawujat=dtset%pawujat
2644    dtpawuj(iuj)%nspden=dtset%nspden
2645    dtpawuj(iuj)%rprimd=dtset%rprimd_orig(1:3,1:3,1)
2646  end do
2647 
2648 !allocate(dtpawuj(0)%vsh(0,0),dtpawuj(0)%occ(0,0))
2649 
2650  do iuj=1,2
2651    if (iuj>1) scfcv_args%cg(:,:)=cgstart(:,:)
2652 
2653 !  DEBUG
2654 !  write(std_out,*)'drive_pawuj before count dtpawuj(:)%iuj ', dtpawuj(:)%iuj
2655 !  END DEBUG
2656 
2657    dtpawuj(iuj*2-1)%iuj=iuj*2-1
2658 
2659    scfcv_args%ndtpawuj=>ndtpawuj
2660    scfcv_args%dtpawuj=>dtpawuj
2661 
2662    !call scfcv_new(ab_scfcv_in,ab_scfcv_inout,dtset,electronpositron,&
2663 !&   paw_dmft,rhog,rhor,rprimd,wffnew,wffnow,xred,xred_old,conv_retcode)
2664    call scfcv_run(scfcv_args,electronpositron,rhog,rhor,rprimd,xred,xred_old,conv_retcode)
2665 
2666    scfcv_args%fatvshift=scfcv_args%fatvshift*(-one)
2667  end do
2668 
2669 !Calculate Hubbard U (or J)
2670  call pawuj_det(dtpawuj,ndtpawuj,trim(scfcv_args%dtfil%filnam_ds(4))//"_UJDET.nc",ures)
2671  dtset%upawu(dtset%typat(dtset%pawujat),1)=ures/Ha_eV
2672 
2673 !Deallocations
2674  do iuj=0,ndtpawuj
2675    call pawuj_free(dtpawuj(iuj))
2676  end do
2677 
2678  ABI_DATATYPE_DEALLOCATE(dtpawuj)
2679  ABI_DEALLOCATE(cgstart)
2680 
2681  DBG_EXIT("COLL")
2682 
2683 end subroutine pawuj_drive

m_gstate/prtxf [ Functions ]

[ Top ] [ m_gstate ] [ Functions ]

NAME

 prtxf

FUNCTION

 Compute and print out dimensional cartesian coordinates and forces.
 Note: for x=cartesian coordinates, t=reduced coordinates (xred),
 =>
  $ x= R t $
 =>
  $ x(1)=rprimd(1,1) t(1)+rprimd(2,1) t(2)+rprimd(3,1) t(3)$
  etc. Also $ t = (R^{-1}) x$ .
  To convert gradients, $d(E)/dx(n) = [d(E)/dt(m)] [dt(m)/dx(n)]$
  and $ dt(m)/dx(n) = (R^{-1})_{mn} = G_{nm}$ because G is the
  inverse transpose of R.  Finally then
  $d(E)/dx(n) = G_{nm} [d(E)/dt(m)]$.
  The vector $d(E)/dt(m)$ for each atom is input in fred
  (grad. wrt xred).

INPUTS

  fred(3,natom)=gradients of Etot (hartree) wrt xred(3,natom)
  iatfix(3,natom)=1 for each fixed atom along specified
  direction, else 0
  iout=unit number for output file
  iwfrc=controls force output: 0=> no forces output,
                               1=>forces out in eV/A and Ha/bohr,
                               2=>forces out in Ha/bohr
  natom=number of atoms in unit cell
  rprimd(3,3)=dimensional real space primitive translations (bohr)
  xred(3,natom)=relative coordinates of atoms (in terms of prim. transl.)

OUTPUT

  (data written to unit iout)

PARENTS

      clnup1

CHILDREN

      matr3inv,wrtout

SOURCE

2147 subroutine prtxf(fred,iatfix,iout,iwfrc,natom,rprimd,xred)
2148 
2149 
2150 !This section has been created automatically by the script Abilint (TD).
2151 !Do not modify the following lines by hand.
2152 #undef ABI_FUNC
2153 #define ABI_FUNC 'prtxf'
2154 !End of the abilint section
2155 
2156  implicit none
2157 
2158 !Arguments ------------------------------------
2159 !scalars
2160  integer,intent(in) :: iout,iwfrc,natom
2161 !arrays
2162  integer,intent(in) :: iatfix(3,natom)
2163  real(dp),intent(in) :: fred(3,natom),rprimd(3,3),xred(3,natom)
2164 
2165 !Local variables-------------------------------
2166 !scalars
2167  integer :: iatom,mu,unfixd
2168  real(dp) :: convt,fmax,frms
2169  character(len=15) :: format_line21
2170  character(len=15) :: format_line25
2171  character(len=15) :: format_line
2172  character(len=500) :: message
2173 !arrays
2174  real(dp) :: favg(3),favg_out(3),ff(3),gprimd(3,3),xx(3)
2175 
2176 ! ****************************************************************
2177 
2178  format_line21='(i5,1x,3f21.14)'
2179  format_line25='(i5,1x,3f25.14)'
2180 
2181 !Write cartesian coordinates in angstroms
2182  call wrtout(iout,' cartesian coordinates (angstrom) at end:','COLL')
2183  do iatom=1,natom
2184    format_line=format_line21
2185    do mu=1,3
2186      xx(mu)=(rprimd(mu,1)*xred(1,iatom)+&
2187 &     rprimd(mu,2)*xred(2,iatom)+&
2188 &     rprimd(mu,3)*xred(3,iatom))*Bohr_Ang
2189      if(xx(mu)>99999 .or. xx(mu)<-9999)format_line=format_line25
2190    end do
2191    write(message,format_line) iatom,xx
2192    call wrtout(iout,message,'COLL')
2193  end do
2194 
2195 !Optionally write cartesian forces in eV/Angstrom (also provide same in hartree/bohr)
2196  if (iwfrc/=0) then
2197 !  First, provide results in hartree/bohr
2198    write(message, '(a,a)' ) ch10,' cartesian forces (hartree/bohr) at end:'
2199    call wrtout(iout,message,'COLL')
2200    frms=zero
2201    fmax=zero
2202    favg(1)=zero
2203    favg(2)=zero
2204    favg(3)=zero
2205 !  To get cartesian forces from input gradients with respect to
2206 !  dimensionless coordinates xred, multiply by G and negate
2207 !  (see notes at top of this subroutine)
2208    call matr3inv(rprimd,gprimd)
2209 !  First compute (spurious) average force favg
2210    do iatom=1,natom
2211      do mu=1,3
2212        ff(mu)=-(gprimd(mu,1)*fred(1,iatom)+&
2213 &       gprimd(mu,2)*fred(2,iatom)+&
2214 &       gprimd(mu,3)*fred(3,iatom))
2215        favg(mu)=favg(mu)+ff(mu)
2216      end do
2217    end do
2218    favg(1) = favg(1)/dble(natom)
2219    favg(2) = favg(2)/dble(natom)
2220    favg(3) = favg(3)/dble(natom)
2221 
2222 !  Subtract off average force in what follows
2223 !  (avg is also subtracted off in carfor, called by loopcv,
2224 !  called by grad)
2225    unfixd=0
2226    do iatom=1,natom
2227      format_line=format_line21
2228      do mu=1,3
2229        ff(mu)=-(gprimd(mu,1)*fred(1,iatom)+&
2230 &       gprimd(mu,2)*fred(2,iatom)+&
2231 &       gprimd(mu,3)*fred(3,iatom))-favg(mu)
2232        if(ff(mu)>99999 .or. ff(mu)<-9999)format_line=format_line25
2233 !      For rms and max force, include only unfixed components
2234        if (iatfix(mu,iatom) /= 1) then
2235          unfixd=unfixd+1
2236          frms=frms+ff(mu)**2
2237          fmax=max(fmax,abs(ff(mu)))
2238        end if
2239      end do
2240      write(message, format_line) iatom,ff
2241      call wrtout(iout,message,'COLL')
2242    end do
2243    if ( unfixd /= 0 ) frms = sqrt(frms/dble(unfixd))
2244 
2245 !  The average force is obtained from the cancellation of numbers
2246 !  of typical size unity, so an absolute value lower
2247 !  than tol14 is meaningless for the output file.
2248    favg_out(:)=favg(:)
2249    if(abs(favg_out(1))<tol14)favg_out(1)=zero
2250    if(abs(favg_out(2))<tol14)favg_out(2)=zero
2251    if(abs(favg_out(3))<tol14)favg_out(3)=zero
2252 
2253    write(message, '(a,1p,2e14.7,1x,3e11.3,a)' )' frms,max,avg=',frms,fmax,favg_out(1:3),' h/b'
2254    call wrtout(iout,message,'COLL')
2255 
2256    if (iwfrc==1) then
2257 
2258      write(message, '(a,a)' )ch10,' cartesian forces (eV/Angstrom) at end:'
2259      call wrtout(iout,message,'COLL')
2260      convt=Ha_eV/Bohr_Ang
2261 
2262 !    Note: subtract off average force
2263      do iatom=1,natom
2264        format_line=format_line21
2265        do mu=1,3
2266          ff(mu)=(-(gprimd(mu,1)*fred(1,iatom)+&
2267 &         gprimd(mu,2)*fred(2,iatom)+&
2268 &         gprimd(mu,3)*fred(3,iatom))-favg(mu))*convt
2269          if(ff(mu)>99999 .or. ff(mu)<-9999)format_line=format_line25
2270        end do
2271        write(message, format_line) iatom,ff
2272        call wrtout(iout,message,'COLL')
2273      end do
2274      write(message, '(a,1p,2e14.7,1x,3e11.3,a)' )' frms,max,avg=',convt*frms,convt*fmax,convt*favg_out(1:3),' e/A'
2275      call wrtout(iout,message,'COLL')
2276 
2277    end if
2278  end if
2279 
2280 end subroutine prtxf

m_gstate/setup2 [ Functions ]

[ Top ] [ m_gstate ] [ Functions ]

NAME

 setup2

FUNCTION

 Call within main routine for setup of various arrays.

INPUTS

  dtset <type(dataset_type)>=all input variables in this dataset
   | ecut=kinetic energy cutoff for planewave basis (hartree)
   | natom=number of atoms in unit cell
   | nkpt=number of k points
   | wtk(nkpt)=integration weight associated with each k point
   | iscf=parameter controlling scf or non-scf choice
  npwtot(nkpt)=number of planewaves in basis and boundary at each k point
  xred(3,natom)=starting reduced atomic coordinates

OUTPUT

  start(3,natom)=copy of starting xred

PARENTS

      gstate

CHILDREN

      wrtout

SOURCE

1803 subroutine setup2(dtset,npwtot,start,wfs,xred)
1804 
1805 
1806 !This section has been created automatically by the script Abilint (TD).
1807 !Do not modify the following lines by hand.
1808 #undef ABI_FUNC
1809 #define ABI_FUNC 'setup2'
1810 !End of the abilint section
1811 
1812  implicit none
1813 
1814 !Arguments ------------------------------------
1815 !scalars
1816  type(dataset_type),intent(in) :: dtset
1817  type(wvl_wf_type),intent(in) :: wfs
1818 !arrays
1819  integer,intent(in) :: npwtot(dtset%nkpt)
1820  real(dp),intent(in) :: xred(3,dtset%natom)
1821  real(dp),intent(out) :: start(3,dtset%natom)
1822 
1823 !Local variables-------------------------------
1824 !scalars
1825  integer :: ikpt,npw
1826  real(dp) :: arith,geom,wtknrm
1827  character(len=500) :: message
1828 
1829 ! *************************************************************************
1830 
1831 !DEBUG
1832 !write(std_out,*)' setup2 : enter '
1833 !ENDDEBUG
1834 
1835    if (dtset%iscf>=0) then
1836 
1837 !  Copy coordinates into array start
1838      start(:,:)=xred(:,:)
1839 
1840      if (dtset%usewvl == 0) then
1841 !    Get average number of planewaves per k point:
1842 !    both arithmetic and GEOMETRIC averages are desired--
1843 !    need geometric average to use method of Francis and Payne,
1844 !    J. Phys.: Condens. Matter 2, 4395-4404 (1990) [[cite:Francis1990]].
1845 !    Also note: force k point wts to sum to 1 for this averaging.
1846 !    (wtk is not forced to add to 1 in a case with occopt=2)
1847        arith=zero
1848        geom=one
1849        wtknrm=zero
1850        do ikpt=1,dtset%nkpt
1851          npw=npwtot(ikpt)
1852          wtknrm=wtknrm+dtset%wtk(ikpt)
1853          arith=arith+npw*dtset%wtk(ikpt)
1854          geom=geom*npw**dtset%wtk(ikpt)
1855        end do
1856 
1857 !    Enforce normalization of weights to 1
1858        arith=arith/wtknrm
1859        geom=geom**(1.0_dp/wtknrm)
1860 
1861      end if
1862 
1863 !  Ensure portability of output thanks to tol8
1864      if (dtset%usewvl == 0) then
1865        write(message, '(a,2f12.3)' ) &
1866 &       '_setup2: Arith. and geom. avg. npw (full set) are',arith+tol8,geom
1867      else
1868 #if defined HAVE_BIGDFT
1869        write(message, '(a,2I8)' ) ' setup2: nwvl coarse and fine are', &
1870 &       wfs%ks%lzd%Glr%wfd%nvctr_c, wfs%ks%lzd%Glr%wfd%nvctr_f
1871 #endif
1872      end if
1873      call wrtout(ab_out,  message,'COLL')
1874      call wrtout(std_out, message,'COLL')
1875 
1876    end if
1877 
1878 #if !defined HAVE_BIGDFT
1879    if (.false.) write(std_out,*) wfs%ks
1880 #endif
1881 
1882  end subroutine setup2