TABLE OF CONTENTS


ABINIT/m_outscfcv [ Modules ]

[ Top ] [ Modules ]

NAME

  m_outscfcv

FUNCTION

COPYRIGHT

  Copyright (C) 2005-2018 ABINIT group (XG)
  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_outscfcv
27 
28  implicit none
29 
30  private

ABINIT/outscfcv [ Functions ]

[ Top ] [ Functions ]

NAME

 outscfcv

FUNCTION

 Output routine for the scfcv.F90 routine

INPUTS

  atindx(natom)=index table for atoms (see gstate.f)
  atindx1(natom)=index table for atoms, inverse of atindx (see gstate.f)
  cg(2,mcg)=planewave coefficients of wavefunctions (see also side effects)
  compch_fft=compensation charge, from FFT grid
  compch_sph=compensation charge, from sphere
  cprj(natom,mcprj*usecprj)=<p_lmn|Cnk> coefficients for each WF |Cnk>
          and each |p_lmn> non-local projector. See also side effects
  dimcprj(natom*usecprj)=array of dimensions of array cprj (not ordered)
  dmatpawu= fixed occupation matrix of correlated orbitals (DFT+U or DMFT only)
  dtfil <type(datafiles_type)>=variables related to files
  dtset <type(dataset_type)>=all input variables for this dataset
  ecut=cut-off energy for plane wave basis sphere (Ha)
  eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree)
  elfr(nfft,nspden(+1))=electron localization function, real space.
   (+1) if spin-polarized in order to get total, spin up and spin down elf
  etotal=total energy
  gmet(3,3)=reciprocal space metric tensor in bohr**-2.
  gprimd(3,3)=dimensional reciprocal space primitive translations
  grhor(nfft,nspden,3)= gradient of electron density in electrons/bohr**4, real space
  hdr <type(hdr_type)>=the header of wf, den and pot files
  kg(3,mpw*mkmem)=reduced planewave coordinates.
  lrhor(nfft,nspden)= Laplacian of electron density in electrons/bohr**5, real space
  mband=maximum number of bands
  mcg=size of wave-functions array (cg) =mpw*my_nspinor*mband*mkmem*nsppol
  mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol
  mgfftc=maximum size of 1D FFTs for the PAW coarse grid
  mkmem =number of k points treated by this node.
  mpi_enreg=information about MPI parallelization
  mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
  mpw=maximum dimensioned size of npw.
  my_natom=number of atoms treated by current processor
  natom=number of atoms in cell.
  nattyp(ntypat)= # atoms of each type.
  nfft=(effective) number of FFT grid points (for this processor) (see NOTES at beginning of scfcv)
  ngfft(18)=contain all needed information about 3D FFT (see NOTES at beginning of scfcv)
  nhat(nfft,nspden*usepaw)= compensation charge density  (PAW)
  nkpt=number of k points.
  npwarr(nkpt)=number of planewaves in basis at this k point
  nspden=number of spin-density components
  nsppol=1 for unpolarized, 2 for spin-polarized
  nsym=number of symmetries in space group
  ntypat=number of types of atoms in unit cell.
  n3xccc=dimension of the xccc3d array (0 or nfft).
  occ(mband*nkpt*nsppol)=occupation number for each band (usually 2) for each k.
  paw_an(my_natom) <type(paw_an_type)>=paw arrays given on angular mesh
  pawang <type(pawang_type)>=paw angular mesh and related data
  pawfgr(natom) <type(pawfgr_type)>=fine grid parameters and related data
  pawfgrtab(my_natom) <type(pawfgrtab_type)> tables on PAW fine grid
  pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data
  pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data
  pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
  paw_ij(my_natom) <type(paw_ij_type)>=paw arrays given on (i,j) channels
            note:structure factors are given on the coarse grid for PAW
  prtvol=control print volume and debugging output
  psps <type(pseudopotential_type)>=variables related to pseudopotentials
  results_gs <type(results_gs_type)>=results (energy and its components,
     forces and its components, the stress tensor) of a ground-state computation
  rhor(nfft,nspden)=total electron density in electrons/bohr**3, real space.
  rprimd(3,3)=dimensional primitive translations for real space (bohr)
  taur(nfft,nspden)=total kinetic energy density in bohr**(-5), real space.
  ucvol=unit cell volume (bohr**3)
  usecprj=1 if cprj datastructure has been allocated
  vhartr(nfft)=Hartree potential
  vxc(nfft,nspden)=xc potential
  vtrial(nfft,nspden)=the trial potential = vxc + vpsp + vhartr, roughly speaking
  xccc3d(n3xccc)=3D core electron density for XC core correction (bohr^-3)
  xred(3,natom)=reduced dimensionless atomic coordinates

OUTPUT

  (only writing, printing)

SIDE EFFECTS

  electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation
  If prtwant==3 the following quantitities are updated using the unitary transformation
  defining the QP amplitudes in terms of the KS basis set:
   cg(2,mcg)=planewave coefficients of wavefunctions.
   cprj(natom,mcprj*usecpyj)=<p_lmn|Cnk> coefficients for each WF |Cnk> and each |p_lmn> non-local projector

NOTES

   The string passed to fftdatar_write (first argument) gives the name used to store the data in the netcdf file
   The function  varname_from_fname defined in the module m_hdr.F90 gives the mapping between the Abinit
   file extension and the netcdf name e.g. foo_VHXC.nc --> vxc
   This function is used in cut3d so that we can immediately select the data to analyze without having
   to prompt the user. Remember to update varname_from_fname if you add a new file or if you change the
   name of the variable.

PARENTS

      scfcv

CHILDREN

      bonds_lgth_angles,bound_deriv,calc_efg,calc_fc,calcdensph
      compute_coeff_plowannier,crystal_free,crystal_init,datafordmft,denfgr
      destroy_dmft,destroy_oper,destroy_plowannier,dos_calcnwrite,ebands_free
      ebands_init,ebands_interpolate_kpath,ebands_prtbltztrp,ebands_write
      epjdos_free,fatbands_ncwrite,fftdatar_write,free_my_atmtab
      get_my_atmtab,init_dmft,init_oper,init_plowannier,ioarr,mag_constr_e
      mlwfovlp,mlwfovlp_qp,multipoles_out,optics_paw,optics_paw_core
      optics_vloc,out1dm,outkss,outwant,partial_dos_fractions
      partial_dos_fractions_paw,pawmkaewf,pawprt,pawrhoij_copy
      pawrhoij_nullify,posdoppler,poslifetime,print_dmft,prt_cif,prtfatbands
      read_atomden,simpson_int,sort_dp,spline,splint,timab,wrtout,xmpi_sum
      xmpi_sum_master

SOURCE

 152 subroutine outscfcv(atindx1,cg,compch_fft,compch_sph,cprj,dimcprj,dmatpawu,dtfil,dtset,&
 153 & ecut,eigen,electronpositron,elfr,etotal,gmet,gprimd,grhor,hdr,kg,&
 154 & lrhor,mband,mcg,mcprj,mgfftc,mkmem,mpi_enreg,mpsang,mpw,my_natom,natom,&
 155 & nattyp,nfft,ngfft,nhat,nkpt,npwarr,nspden,nsppol,nsym,ntypat,n3xccc,occ,&
 156 & paw_dmft,pawang,pawfgr,pawfgrtab,pawrad,pawrhoij,pawtab,paw_an,paw_ij,&
 157 & prtvol,psps,results_gs,rhor,rprimd,&
 158 & taur,ucvol,usecprj,vhartr,vpsp,vtrial,vxc,wvl_den,xccc3d,xred)
 159 
 160  use defs_basis
 161  use defs_datatypes
 162  use defs_wvltypes
 163  use defs_abitypes
 164  use m_abicore
 165  use m_sort
 166  use m_efield
 167  use m_errors
 168  use m_xmpi
 169  use m_mpinfo
 170 #ifdef HAVE_NETCDF
 171  use netcdf
 172 #endif
 173  use m_abi_etsf
 174  use m_nctk
 175  use m_hdr
 176  use m_plowannier
 177  use m_splines
 178  use m_ebands
 179 
 180  use m_time,             only : timab
 181  use m_io_tools,         only : open_file
 182  use m_fstrings,         only : strcat, endswith
 183  use m_geometry,         only : bonds_lgth_angles
 184  use m_electronpositron, only : electronpositron_type,electronpositron_calctype
 185  use m_oper,             only : oper_type,init_oper,destroy_oper
 186  use m_crystal,          only : crystal_init, crystal_free, crystal_t, prt_cif
 187  use m_crystal_io,       only : crystal_ncwrite
 188  use m_results_gs,       only : results_gs_type, results_gs_ncwrite
 189  use m_ioarr,            only : ioarr, fftdatar_write
 190  use m_nucprop,          only : calc_efg,calc_fc
 191  use m_outwant,          only : outwant
 192  use m_pawang,           only : pawang_type
 193  use m_pawrad,           only : pawrad_type, simp_gen, bound_deriv
 194  use m_pawtab,           only : pawtab_type
 195  use m_paw_an,           only : paw_an_type
 196  use m_paw_ij,           only : paw_ij_type
 197  use m_paw_mkrho,        only : denfgr
 198  use m_pawfgrtab,        only : pawfgrtab_type
 199  use m_pawrhoij,         only : pawrhoij_type, pawrhoij_nullify, pawrhoij_copy, pawrhoij_free
 200  use m_pawcprj,          only : pawcprj_type
 201  use m_pawfgr,           only : pawfgr_type
 202  use m_paw_dmft,         only : paw_dmft_type,init_dmft,destroy_dmft,print_dmft
 203  use m_paw_optics,       only : optics_paw,optics_paw_core
 204  use m_paw_tools,        only : pawprt
 205  use m_numeric_tools,    only : simpson_int
 206  use m_epjdos,           only : dos_calcnwrite, partial_dos_fractions, partial_dos_fractions_paw, &
 207                                 epjdos_t, epjdos_new, epjdos_free, prtfatbands, fatbands_ncwrite
 208  use m_paral_atom,       only : get_my_atmtab, free_my_atmtab
 209  use m_io_kss,           only : outkss
 210  use m_multipoles,       only : multipoles_out, out1dm
 211  use m_mlwfovlp_qp,      only : mlwfovlp_qp
 212  use m_paw_mkaewf,       only : pawmkaewf
 213  use m_dens,             only : mag_constr_e, calcdensph
 214  use m_mlwfovlp,         only : mlwfovlp
 215  use m_datafordmft,      only : datafordmft
 216  use m_mkrho,            only : read_atomden
 217  use m_positron,         only : poslifetime, posdoppler
 218  use m_optics_vloc,      only : optics_vloc
 219 
 220 !This section has been created automatically by the script Abilint (TD).
 221 !Do not modify the following lines by hand.
 222 #undef ABI_FUNC
 223 #define ABI_FUNC 'outscfcv'
 224 !End of the abilint section
 225 
 226  implicit none
 227 
 228 !Arguments ------------------------------------
 229 !scalars
 230  integer,intent(in) :: mband,mcg,mcprj,mgfftc,mkmem,mpsang,mpw,n3xccc,my_natom,natom,nfft
 231  integer,intent(in) :: nkpt,nspden,nsppol,nsym,ntypat,prtvol,usecprj
 232  real(dp),intent(in) :: compch_fft,compch_sph,ecut,ucvol
 233  real(dp),intent(inout) :: etotal
 234  type(electronpositron_type),pointer :: electronpositron
 235  type(MPI_type),intent(inout) :: mpi_enreg
 236  type(datafiles_type),intent(in) :: dtfil
 237  type(dataset_type),intent(in) :: dtset
 238  type(hdr_type),intent(inout) :: hdr
 239  type(paw_dmft_type), intent(inout)  :: paw_dmft
 240  type(pawang_type),intent(in) :: pawang
 241  type(pawfgr_type),intent(in) :: pawfgr
 242  type(pseudopotential_type),intent(in) :: psps
 243  type(results_gs_type),intent(in) :: results_gs
 244  type(wvl_denspot_type), intent(in) :: wvl_den
 245 !arrays
 246  integer,intent(in) :: atindx1(natom),dimcprj(natom*usecprj)
 247  integer,intent(in) :: kg(3,mpw*mkmem),nattyp(ntypat),ngfft(18),npwarr(nkpt)
 248  real(dp),intent(in) :: dmatpawu(:,:,:,:),eigen(mband*nkpt*nsppol)
 249  real(dp),intent(in) :: gmet(3,3),gprimd(3,3)
 250  real(dp),intent(in) :: occ(mband*nkpt*nsppol)
 251  real(dp),intent(in) :: rprimd(3,3),vhartr(nfft),xccc3d(n3xccc)
 252  real(dp),intent(in) :: vpsp(nfft)
 253  real(dp),intent(inout) :: cg(2,mcg)
 254  real(dp),intent(inout) :: nhat(nfft,nspden*psps%usepaw)
 255  real(dp),intent(inout),target :: rhor(nfft,nspden),vtrial(nfft,nspden)
 256  real(dp),intent(inout) :: vxc(nfft,nspden),xred(3,natom)
 257  real(dp),pointer :: elfr(:,:),grhor(:,:,:),lrhor(:,:),taur(:,:)
 258  type(pawcprj_type),intent(inout) :: cprj(natom,mcprj*usecprj)
 259  type(paw_an_type),intent(inout) :: paw_an(my_natom*psps%usepaw)
 260  type(pawfgrtab_type),intent(in) :: pawfgrtab(my_natom*psps%usepaw)
 261  type(paw_ij_type),intent(inout) :: paw_ij(my_natom*psps%usepaw)
 262  type(pawrad_type),intent(in) :: pawrad(psps%ntypat*psps%usepaw)
 263  type(pawrhoij_type),target,intent(inout) :: pawrhoij(my_natom*psps%usepaw)
 264  type(pawtab_type),intent(inout) :: pawtab(ntypat*psps%usepaw)
 265 
 266 !Local variables-------------------------------
 267 !scalars
 268  integer,parameter :: master=0,cplex1=1,fform_den=52,rdwr2=2,rdwrpaw0=0
 269  integer :: bantot,fform,collect,timrev
 270  integer :: accessfil,coordn
 271  integer :: ii,ierr,ifft,ikpt,ispden,isppol,itypat
 272  integer :: me_fft,n1,n2,n3
 273  integer :: ifgd, iatom, iatom_tot,nradint
 274  integer :: me,my_natom_tmp
 275  integer :: occopt
 276  integer :: prtnabla
 277  integer :: pawprtden
 278  integer :: iband,nocc,spacecomm,comm_fft,tmp_unt,nfft_tot
 279  integer :: my_comm_atom
 280 #ifdef HAVE_NETCDF
 281  integer :: ncid
 282 #endif
 283  real(dp) :: norm,occ_norm,unocc_norm
 284  real(dp) :: rate_dum,rate_dum2
 285  real(dp) :: yp1, ypn, dr
 286  character(len=500) :: message
 287  character(len=fnlen) :: fname
 288 !arrays
 289  integer, allocatable :: isort(:)
 290  integer, pointer :: my_atmtab(:)
 291  real(dp) :: tsec(2),nt_ntone_norm(nspden)
 292  real(dp),allocatable :: eigen2(:)
 293  real(dp),allocatable :: elfr_down(:,:),elfr_up(:,:)
 294  real(dp),allocatable :: rhor_paw(:,:),rhor_paw_core(:,:),rhor_paw_val(:,:),vpaw(:,:),vwork(:,:)
 295  real(dp),allocatable :: rhor_n_one(:,:),rhor_nt_one(:,:),ps_norms(:,:,:)
 296  real(dp), allocatable :: doccde(:)
 297  real(dp), allocatable :: vh1spl(:)
 298  real(dp), allocatable :: vh1_interp(:)
 299  real(dp), allocatable :: vh1_integ(:)
 300  real(dp), allocatable :: vh1_corrector(:)
 301  real(dp), allocatable :: radii(:)
 302  real(dp), ABI_CONTIGUOUS pointer :: rho_ptr(:,:)
 303  type(pawrhoij_type) :: pawrhoij_dum(0)
 304  type(pawrhoij_type),pointer :: pawrhoij_all(:)
 305  logical :: remove_inv
 306  logical :: paral_atom, paral_fft, my_atmtab_allocated
 307  real(dp) :: e_fermie
 308  type(oper_type) :: lda_occup
 309  type(crystal_t) :: crystal
 310  type(ebands_t) :: ebands
 311  type(epjdos_t) :: dos
 312  type(plowannier_type) :: wan
 313 
 314 ! *************************************************************************
 315 
 316  DBG_ENTER("COLL")
 317  call timab(950,1,tsec) ! outscfcv
 318 
 319  if ((usecprj==0.or.mcprj==0).and.psps%usepaw==1.and. &
 320 & (dtset%prtwant==2.or.dtset%prtwant==3.or.dtset%prtnabla>0.or.dtset%prtdos==3 &
 321 & .or.dtset%kssform==3.or.dtset%pawfatbnd>0.or.dtset%pawprtwf>0)) then
 322    write (message,'(5a)')&
 323 &   'cprj datastructure must be allocated',ch10,&
 324 &   'with options prtwant=2,3, prtnabla>0, prtdos>3, kssform==3, pawfatbnd>0, pawprtwf>0',ch10,&
 325 &   'Action: change pawusecp input keyword.'
 326    MSG_ERROR(message)
 327  end if
 328 
 329 !Initialize two objects to facilitate the propagation of info.
 330 !These objects should used more frequently, actually they should
 331 !become basic objects used in abinit.
 332 
 333 !Crystalline structure.
 334  remove_inv=.false.
 335  if(dtset%nspden==4 .and. dtset%usedmft==1) remove_inv=.true. ! MG: why this?
 336 
 337  timrev = 2; if (any(dtset%kptopt == [3, 4])) timrev= 1
 338  call crystal_init(dtset%amu_orig(:,1),crystal,dtset%spgroup,natom,dtset%npsp,ntypat, &
 339 & dtset%nsym,rprimd,dtset%typat,xred,dtset%ziontypat,dtset%znucl,timrev,&
 340 & dtset%nspden==2.and.dtset%nsppol==1,remove_inv,hdr%title,&
 341 & dtset%symrel,dtset%tnons,dtset%symafm)
 342 
 343 !Electron band energies.
 344  bantot= dtset%mband*dtset%nkpt*dtset%nsppol
 345  ABI_MALLOC(doccde,(bantot))
 346  doccde=zero
 347 
 348  call ebands_init(bantot,ebands,dtset%nelect,doccde,eigen,hdr%istwfk,hdr%kptns,hdr%nband,&
 349 & hdr%nkpt,hdr%npwarr,hdr%nsppol,hdr%nspinor,hdr%tphysel,hdr%tsmear,hdr%occopt,hdr%occ,hdr%wtk,&
 350 & hdr%charge, hdr%kptopt, hdr%kptrlatt_orig, hdr%nshiftk_orig, hdr%shiftk_orig, &
 351 & hdr%kptrlatt, hdr%nshiftk, hdr%shiftk)
 352 
 353  ABI_FREE(doccde)
 354 
 355  ebands%fermie  = results_gs%energies%e_fermie
 356  e_fermie = results_gs%energies%e_fermie
 357  ebands%entropy = results_gs%energies%entropy
 358  !write(std_out,*)"ebands%efermi in outscfcv",ebands%fermie
 359  !write(std_out,*)"results_gs%energies%e_fermie in outscfcv",e_fermie
 360  !write(std_out,*)"results_gs%fermie in outscfcv",results_gs%fermie
 361  !write(std_out,*)"hdr%efermi in outscfcv",hdr%fermie
 362 
 363  ! Parameters for MPI-FFT
 364  n1 = ngfft(1); n2 = ngfft(2); n3 = ngfft(3); nfft_tot = product(ngfft(1:3))
 365  comm_fft = mpi_enreg%comm_fft
 366  me_fft = xmpi_comm_rank(comm_fft)
 367  paral_fft = (mpi_enreg%paral_kgb==1)
 368 
 369  spacecomm = mpi_enreg%comm_cell
 370  me = xmpi_comm_rank(spacecomm)
 371 
 372  paral_atom=(my_natom/=natom)
 373  my_comm_atom = mpi_enreg%comm_atom
 374  nullify(my_atmtab)
 375  if (paral_atom) then
 376    call get_my_atmtab(mpi_enreg%comm_atom, my_atmtab, my_atmtab_allocated, paral_atom,natom,my_natom_ref=my_natom)
 377  else
 378    ABI_ALLOCATE(my_atmtab, (natom))
 379    my_atmtab = (/ (iatom, iatom=1, natom) /)
 380    my_atmtab_allocated = .true.
 381  end if
 382 
 383 !wannier interface
 384  call timab(951,1,tsec)
 385  if (dtset%prtwant==2) then
 386 
 387    call mlwfovlp(crystal, ebands, hdr, atindx1,cg,cprj,dtset,dtfil,eigen,gprimd,kg,&
 388 &   mband,mcg,mcprj,mgfftc,mkmem,mpi_enreg,mpw,natom,&
 389 &   nattyp,nfft,ngfft,nkpt,npwarr,nsppol,ntypat,occ,&
 390 &   pawang,pawrad,pawtab,prtvol,psps,rprimd,ucvol,xred)
 391 
 392  else if (dtset%prtwant==3) then
 393 
 394 !  Convert cg and eigen to GW quasiparticle wave functions and eigenvalues in mlwfovlp_qp
 395    ABI_ALLOCATE(eigen2,(mband*nkpt*nsppol))
 396    eigen2=eigen
 397 
 398    call mlwfovlp_qp(cg,cprj,dtset,dtfil,eigen2,mband,mcg,mcprj,mkmem,mpw,natom,&
 399 &   nkpt,npwarr,nspden,nsppol,ntypat,Hdr,pawtab,rprimd,MPI_enreg)
 400 
 401 !  Call Wannier90
 402    call mlwfovlp(crystal, ebands, hdr, atindx1,cg,cprj,dtset,dtfil,eigen2,gprimd,kg,&
 403 &   mband,mcg,mcprj,mgfftc,mkmem,mpi_enreg,mpw,natom,&
 404 &   nattyp,nfft,ngfft,nkpt,npwarr,nsppol,ntypat,occ,&
 405 &   pawang,pawrad,pawtab,prtvol,psps,rprimd,ucvol,xred)
 406 
 407 !  this is the old implementation, risky due to unpredictable size effects
 408 !  now eigen is not overwritten, one should use other ways to print the GW corrections
 409 !  eigen=eigen2
 410    ABI_DEALLOCATE(eigen2)
 411  end if !prtwant
 412  call timab(951,2,tsec)
 413 
 414  occopt=dtset%occopt
 415 
 416  prtnabla=dtset%prtnabla
 417  pawprtden=dtset%prtden-1
 418 
 419  call timab(952,1,tsec)
 420 
 421  spacecomm=mpi_enreg%comm_cell; me=xmpi_comm_rank(spacecomm)
 422  comm_fft=mpi_enreg%comm_fft
 423  paral_atom=(my_natom/=natom)
 424 
 425 !Warnings :
 426 !- core charge is excluded from the charge density;
 427 !- the potential is the INPUT vtrial.
 428 
 429  if (iwrite_fftdatar(mpi_enreg) .and. dtset%usewvl==0) then
 430 
 431    ! output the density.
 432    if (dtset%prtden/=0) then
 433      if (dtset%positron/=1) rho_ptr => rhor
 434      if (dtset%positron==1) rho_ptr => electronpositron%rhor_ep
 435      call fftdatar_write("density",dtfil%fnameabo_app_den,dtset%iomode,hdr,&
 436      crystal,ngfft,cplex1,nfft,nspden,rho_ptr,mpi_enreg,ebands=ebands)
 437 
 438      if (dtset%positron/=0) then
 439        if (dtset%positron/=1) rho_ptr => electronpositron%rhor_ep
 440        if (dtset%positron==1) rho_ptr => rhor
 441        fname = trim(dtfil%fnameabo_app_den)//'_POSITRON'
 442        if (dtset%iomode == IO_MODE_ETSF) fname = strcat(fname, ".nc")
 443        call fftdatar_write("positron_density",fname,dtset%iomode,hdr,&
 444        crystal,ngfft,cplex1,nfft,nspden,rho_ptr,mpi_enreg,ebands=ebands)
 445      end if
 446    end if
 447 
 448  else if (dtset%usewvl == 1 .and. dtset%prtden /= 0) then
 449    !if iomode == 2 then set all outputs to netcdf format
 450    !if iomode == 3 then set all outputs to ETSF format
 451    accessfil = 0
 452    if (dtset%iomode == IO_MODE_ETSF) accessfil = 3
 453    if (dtset%iomode == IO_MODE_MPI) accessfil = 4
 454    fform = fform_den
 455     ! Write wavelet DEN. Note however that this should be delegate to separated Bigdft routines.
 456     ! a lot of stuff written in outscf does not make sense if usewvl==0
 457    call ioarr(accessfil,rhor,dtset,etotal,fform,dtfil%fnameabo_app_den, &
 458    hdr,mpi_enreg,ngfft,cplex1,nfft,pawrhoij_dum,rdwr2,rdwrpaw0,wvl_den)
 459  end if ! if master
 460 
 461 !! MS - Printing of PAWDEN parallellised and several possible options included
 462 !We output the total electron density in the PAW case
 463 !this requires removing nhat from rhor and making PAW on-site corrections
 464  if (pawprtden>0 .and. psps%usepaw==1) then
 465 !  pawprtden 1 --> output PAW valence density
 466 !  "     2 --> output PAW valence+core density
 467 !  "     3 --> output core, valence and full atomic protodensity
 468 !  "     4 --> options 1+3
 469 !  "     5 --> options 2+3
 470 !  "     6 --> output all individual PAW density contributions
 471    if (pawprtden/=3) then ! calc PAW valence density
 472      ABI_ALLOCATE(rhor_paw,(pawfgr%nfft,nspden))
 473      ABI_ALLOCATE(rhor_n_one,(pawfgr%nfft,nspden))
 474      ABI_ALLOCATE(rhor_nt_one,(pawfgr%nfft,nspden))
 475 !    If the communicator used for denfgr is kpt_comm, it is not compatible with paral_atom
 476      if (mpi_enreg%paral_kgb==0.and.my_natom/=natom) then
 477        my_natom_tmp=natom
 478        ABI_DATATYPE_ALLOCATE(pawrhoij_all,(natom))
 479        call pawrhoij_nullify(pawrhoij_all)
 480        call pawrhoij_copy(pawrhoij,pawrhoij_all,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,&
 481 &       keep_cplex=.false.,keep_itypat=.false.,keep_nspden=.false.)
 482      else
 483        my_natom_tmp=my_natom
 484        pawrhoij_all => pawrhoij
 485      end if
 486      if (pawprtden/=6) then
 487        call denfgr(atindx1,gmet,comm_fft,my_natom_tmp,natom,nattyp,ngfft,nhat,dtset%nspinor,nsppol,nspden,&
 488 &       ntypat,pawfgr,pawrad,pawrhoij_all,pawtab,prtvol,rhor,rhor_paw,rhor_n_one,&
 489 &       rhor_nt_one,rprimd,dtset%typat,ucvol,xred,&
 490 &       comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
 491      else
 492        call denfgr(atindx1,gmet,comm_fft,my_natom_tmp,natom,nattyp,ngfft,nhat,dtset%nspinor,nsppol,nspden,&
 493 &       ntypat,pawfgr,pawrad,pawrhoij_all,pawtab,prtvol,rhor,rhor_paw,rhor_n_one,&
 494 &       rhor_nt_one,rprimd,dtset%typat,ucvol,xred,&
 495 &       abs_n_tilde_nt_diff=nt_ntone_norm,znucl=dtset%znucl,&
 496 &       comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
 497      end if
 498      if (mpi_enreg%paral_kgb==0.and.my_natom/=natom) then
 499        call pawrhoij_free(pawrhoij_all)
 500        ABI_DATATYPE_DEALLOCATE(pawrhoij_all)
 501      end if
 502 
 503      if (prtvol>9) then  ! Check normalisation
 504        norm = SUM(rhor_paw(:,1))*ucvol/PRODUCT(pawfgr%ngfft(1:3))
 505        call xmpi_sum(norm,comm_fft,ierr)
 506        write(message,'(a,F8.4)') '  PAWDEN - NORM OF DENSITY: ',norm
 507        call wrtout(std_out,message,'COLL')
 508      end if
 509    end if
 510 
 511    if (pawprtden>1.AND.pawprtden<6) then ! We will need the core density
 512      ABI_ALLOCATE(rhor_paw_core,(pawfgr%nfft,nspden))
 513      call read_atomden(mpi_enreg,natom,pawfgr%nfft,pawfgr%ngfft,nspden,ntypat,rhor_paw_core,&
 514 &     dtset%typat,rprimd,xred,prtvol,file_prefix='core   ')
 515 
 516      if (prtvol>9) then  ! Check normalisation
 517        norm = SUM(rhor_paw_core(:,1))*ucvol/PRODUCT(pawfgr%ngfft(1:3))
 518        call xmpi_sum(norm,comm_fft,ierr)
 519        write(message,'(a,F8.4)') '  ATMDEN - NORM OF CORE DENSITY: ', norm
 520        call wrtout(std_out,message,'COLL')
 521      end if
 522    end if
 523 
 524    if (pawprtden>2.AND.pawprtden<6) then ! We will need the valence protodensity
 525      ABI_ALLOCATE(rhor_paw_val,(pawfgr%nfft,nspden))
 526      call read_atomden(mpi_enreg,natom,pawfgr%nfft,pawfgr%ngfft,nspden,ntypat,rhor_paw_val,&
 527 &     dtset%typat,rprimd,xred,prtvol,file_prefix='valence')
 528 
 529      if (prtvol>9) then ! Check normalisation
 530        norm = SUM(rhor_paw_val(:,1))*ucvol/PRODUCT(pawfgr%ngfft(1:3))
 531        call xmpi_sum(norm,comm_fft,ierr)
 532        write(message,'(a,F8.4)') '  ATMDEN - NORM OF VALENCE PROTODENSITY: ', norm
 533        call wrtout(std_out,message,'COLL')
 534      end if
 535    end if
 536 
 537    if (iwrite_fftdatar(mpi_enreg)) then
 538      if (pawprtden/=3) then
 539        if (pawprtden==2.or.pawprtden==5) rhor_paw = rhor_paw + rhor_paw_core
 540 !      PAWDEN
 541        call fftdatar_write("pawrhor",dtfil%fnameabo_app_pawden,dtset%iomode,hdr,&
 542        crystal,ngfft,cplex1,nfft,nspden,rhor_paw,mpi_enreg,ebands=ebands)
 543      end if
 544 
 545      if (pawprtden>2.AND.pawprtden<6) then
 546        ! ATMDEN_CORE
 547        call fftdatar_write("pawrhor_core",dtfil%fnameabo_app_atmden_core,dtset%iomode,hdr,&
 548        crystal,ngfft,cplex1,nfft,nspden,rhor_paw_core,mpi_enreg,ebands=ebands)
 549 
 550        ! valence protodensity. ATMDEN_VAL
 551        call fftdatar_write("pawrhor_val",dtfil%fnameabo_app_atmden_val,dtset%iomode,hdr,&
 552        crystal,ngfft,cplex1,nfft,nspden,rhor_paw_val,mpi_enreg,ebands=ebands)
 553 
 554        ! full protodensity. ATMDEN_FULL
 555        rhor_paw_val = rhor_paw_val + rhor_paw_core
 556        call fftdatar_write("pawrhor_full",dtfil%fnameabo_app_atmden_full,dtset%iomode,hdr,&
 557        crystal,ngfft,cplex1,nfft,nspden,rhor_paw_val,mpi_enreg,ebands=ebands)
 558      end if
 559 
 560      if (pawprtden==6) then ! Print all individual contributions to the density
 561        ! N_TILDE - N_HAT
 562        ! Use rhor_paw_val as temporary array
 563        if (.not.allocated(rhor_paw_val))  then
 564          ABI_ALLOCATE(rhor_paw_val,(pawfgr%nfft,nspden))
 565        end if
 566        rhor_paw_val = rhor - nhat
 567 
 568        call fftdatar_write("pawrhor_ntilde_minus_nhat",dtfil%fnameabo_app_n_tilde,dtset%iomode,hdr,&
 569        crystal,ngfft,cplex1,nfft,nspden,rhor_paw_val,mpi_enreg,ebands=ebands)
 570 
 571 !      N_ONSITE
 572        call fftdatar_write("pawrhor_n_one",dtfil%fnameabo_app_n_one,dtset%iomode,hdr,&
 573        crystal,ngfft,cplex1,nfft,nspden,rhor_n_one,mpi_enreg,ebands=ebands)
 574 
 575 !      N_TILDE_ONSITE
 576        call fftdatar_write("pawrhor_nt_one",dtfil%fnameabo_app_nt_one,dtset%iomode,hdr,&
 577        crystal,ngfft,cplex1,nfft,nspden,rhor_nt_one,mpi_enreg,ebands=ebands)
 578 
 579      end if ! All indivdual density cont.
 580    end if ! if master
 581 
 582    if (allocated(rhor_paw))  then
 583      ABI_DEALLOCATE(rhor_paw)
 584    end if
 585    if (allocated(rhor_paw_core))  then
 586      ABI_DEALLOCATE(rhor_paw_core)
 587    end if
 588    if (allocated(rhor_paw_val))  then
 589      ABI_DEALLOCATE(rhor_paw_val)
 590    end if
 591    if (allocated(rhor_n_one))  then
 592      ABI_DEALLOCATE(rhor_n_one)
 593    end if
 594    if (allocated(rhor_nt_one))  then
 595      ABI_DEALLOCATE(rhor_nt_one)
 596    end if
 597 
 598  end if ! if paw+pawprtden
 599 
 600  call timab(952,2,tsec)
 601  call timab(953,1,tsec)
 602 
 603  ! Output of the GSR file (except when we are inside mover)
 604 #ifdef HAVE_NETCDF
 605  if (me == master .and. dtset%prtgsr==1 .and. dtset%usewvl == 0) then
 606    !.and. (dtset%ionmov /= 0 .or. dtset%optcell /= 0)) then
 607    fname = strcat(dtfil%filnam_ds(4), "_GSR.nc")
 608 
 609    ! Write density, crystal, band structure energies.
 610    NCF_CHECK(nctk_open_create(ncid, fname, xmpi_comm_self))
 611    NCF_CHECK(hdr_ncwrite(hdr, ncid, fform_den, nc_define=.True.))
 612    NCF_CHECK(crystal_ncwrite(crystal, ncid))
 613    NCF_CHECK(ebands_ncwrite(ebands, ncid))
 614    ! Add energy, forces, stresses
 615    NCF_CHECK(results_gs_ncwrite(results_gs, ncid, dtset%ecut, dtset%pawecutdg))
 616    NCF_CHECK(nf90_close(ncid))
 617  end if
 618 #endif
 619 
 620  ! Output of VCLMB file
 621  ! The PAW correction has to be computed here (all processors contribute)
 622  if (psps%usepaw > 0 .AND. dtset%prtvclmb>0) then
 623    nradint = 1000 ! radial integration grid density
 624    ABI_ALLOCATE(vpaw,(nfft,nspden))
 625    vpaw(:,:)=zero
 626    if (me == master .and. my_natom > 0) then
 627      if (paw_an(1)%cplex > 1) then
 628        MSG_WARNING('cplex = 2 : complex hartree potential in PAW spheres. This is not coded yet. Imag part ignored')
 629      end if
 630    end if
 631 
 632    do ispden=1,nspden
 633      ! for points inside spheres, replace with full AE hartree potential.
 634      ! In principle the correction could be more subtle (not spherical)
 635      do iatom=1,my_natom
 636        iatom_tot=iatom;if (paral_atom) iatom_tot=mpi_enreg%my_atmtab(iatom)
 637        itypat=dtset%typat(iatom_tot)
 638 
 639        ABI_ALLOCATE(vh1spl,(paw_an(iatom)%mesh_size))
 640        ABI_ALLOCATE(vh1_corrector,(paw_an(iatom)%mesh_size))
 641        ABI_ALLOCATE(vh1_interp,(pawfgrtab(iatom)%nfgd))
 642        ABI_ALLOCATE(radii,(pawfgrtab(iatom)%nfgd))
 643        ABI_ALLOCATE(isort,(pawfgrtab(iatom)%nfgd))
 644        ! vh1 vht1 contain the spherical first moments of the Hartree potentials, so re-divide by Y_00 = sqrt(four_pi)
 645        vh1_corrector(:) = (paw_an(iatom)%vh1(:,1,ispden)-paw_an(iatom)%vht1(:,1,ispden)) / sqrt(four_pi)
 646 
 647        ! get end point derivatives
 648        call bound_deriv(vh1_corrector, pawrad(itypat), pawrad(itypat)%mesh_size, yp1, ypn)
 649        ! spline the vh1 function
 650        ! NB for second argument of vh1: only first moment lm_size appears to be used
 651        ! NB2: vh1 can in principle be complex - not sure what to do with the imaginary part. Ignored for now.
 652        call spline(pawrad(itypat)%rad, vh1_corrector, paw_an(iatom)%mesh_size, yp1, ypn, vh1spl)
 653 
 654        do ifgd = 1, pawfgrtab(iatom)%nfgd
 655          ! get radii for this point
 656          isort(ifgd) = ifgd
 657          radii(ifgd) = sqrt(sum(pawfgrtab(iatom)%rfgd(:,ifgd)**2))
 658        end do
 659 
 660        if (pawfgrtab(iatom)%nfgd/=0) then
 661        ! spline interpolate the vh1 value for current radii
 662          call sort_dp(pawfgrtab(iatom)%nfgd, radii, isort, tol12)
 663          call splint(pawrad(itypat)%mesh_size, pawrad(itypat)%rad, &
 664 &         vh1_corrector, vh1spl, pawfgrtab(iatom)%nfgd, radii,  vh1_interp, ierr)
 665        end if
 666 
 667        norm=SUM(vh1_interp)*ucvol/PRODUCT(ngfft(1:3))
 668        call xmpi_sum(norm,comm_fft,ierr)
 669        write(message,'(a,i6,a,E20.10)') ' sum of Hartree correction term on fft grid of atom : ', iatom, &
 670 &       ' = ', norm
 671        call wrtout(std_out,message,'COLL')
 672 
 673        if (pawfgrtab(iatom)%nfgd/=0) then
 674          vpaw(pawfgrtab(iatom)%ifftsph(isort(1:pawfgrtab(iatom)%nfgd)),ispden) = &
 675 &         vpaw(pawfgrtab(iatom)%ifftsph(isort(1:pawfgrtab(iatom)%nfgd)),ispden) + &
 676 &         vh1_interp(1:pawfgrtab(iatom)%nfgd)
 677        end if
 678 
 679        ! get integral of correction term in whole sphere
 680        ABI_DEALLOCATE(radii)
 681        ABI_DEALLOCATE(vh1_interp)
 682 
 683        ABI_ALLOCATE(radii,(nradint))
 684        ABI_ALLOCATE(vh1_interp,(nradint))
 685 
 686        ABI_ALLOCATE(vh1_integ,(nradint))
 687        dr = pawrad(itypat)%rad(paw_an(iatom)%mesh_size) / dble(nradint)
 688        do ifgd = 1, nradint
 689          radii(ifgd) = dble(ifgd-1)*dr
 690        end do
 691 
 692        ! spline interpolate the vh1 value for current radii
 693        call splint(pawrad(itypat)%mesh_size, pawrad(itypat)%rad, &
 694 &       vh1_corrector, vh1spl, nradint, radii,  vh1_interp, ierr)
 695 
 696        do ifgd = 1, nradint
 697          vh1_interp(ifgd) = vh1_interp(ifgd)*radii(ifgd)**2
 698        end do
 699 
 700        call simpson_int(nradint, dr, vh1_interp, vh1_integ)
 701        write(message,'(a,i6,a,E20.10)') ' integral of Hartree correction term in sphere of atom: ', iatom, &
 702 &       ' = ', vh1_integ(nradint)*four*pi
 703        call wrtout(std_out,message,'COLL')
 704 
 705        ABI_DEALLOCATE(vh1spl)
 706        ABI_DEALLOCATE(vh1_corrector)
 707        ABI_DEALLOCATE(vh1_interp)
 708        ABI_DEALLOCATE(vh1_integ)
 709        ABI_DEALLOCATE(radii)
 710        ABI_DEALLOCATE(isort)
 711      end do ! iatom
 712    end do !ispden
 713    call xmpi_sum_master(vpaw,master,mpi_enreg%comm_atom,ierr)
 714    if (.not.iwrite_fftdatar(mpi_enreg)) then
 715      ABI_DEALLOCATE(vpaw)
 716    end if
 717  end if ! if paw - add all electron vhartree in spheres
 718 
 719  if (iwrite_fftdatar(mpi_enreg)) then
 720 
 721    ! output the electron localization function ELF
 722    if (dtset%prtelf/=0) then
 723      call fftdatar_write("elfr",dtfil%fnameabo_app_elf,dtset%iomode,hdr,&
 724      crystal,ngfft,cplex1,nfft,nspden,elfr,mpi_enreg,ebands=ebands)
 725 
 726      if (nspden==2)then
 727        ABI_ALLOCATE(elfr_up,(nfft,nspden))
 728        elfr_up(:,:) = zero
 729        do ifft=1,nfft
 730          elfr_up(ifft,1) = elfr(ifft,2)
 731        end do
 732 !      ELF_UP
 733        call fftdatar_write("elfr_up",dtfil%fnameabo_app_elf_up,dtset%iomode,hdr,&
 734        crystal,ngfft,cplex1,nfft,nspden,elfr_up,mpi_enreg,ebands=ebands)
 735 
 736        ABI_ALLOCATE(elfr_down,(nfft,nspden))
 737        elfr_down(:,:) = zero
 738        do ifft=1,nfft
 739          elfr_down(ifft,1) = elfr(ifft,3)
 740        end do
 741 !      ELF_DOWN'
 742        call fftdatar_write("elfr_down",dtfil%fnameabo_app_elf_down,dtset%iomode,hdr,&
 743        crystal,ngfft,cplex1,nfft,nspden,elfr_down,mpi_enreg,ebands=ebands)
 744 
 745        ABI_DEALLOCATE(elfr_up)
 746        ABI_DEALLOCATE(elfr_down)
 747      end if
 748    end if
 749 
 750    call timab(953,2,tsec)
 751    call timab(954,1,tsec)
 752 
 753 !  We output the gradient of density
 754    if (dtset%prtgden/=0) then
 755 
 756      call fftdatar_write("grhor_1",dtfil%fnameabo_app_gden1,dtset%iomode,hdr,&
 757      crystal,ngfft,cplex1,nfft,nspden,grhor(:,:,1),mpi_enreg,ebands=ebands)
 758 
 759      call fftdatar_write("grhor_2",dtfil%fnameabo_app_gden2,dtset%iomode,hdr,&
 760      crystal,ngfft,cplex1,nfft,nspden,grhor(:,:,2),mpi_enreg,ebands=ebands)
 761 
 762      call fftdatar_write("grhor_3",dtfil%fnameabo_app_gden3,dtset%iomode,hdr,&
 763      crystal,ngfft,cplex1,nfft,nspden,grhor(:,:,3),mpi_enreg,ebands=ebands)
 764    end if
 765 
 766 !  We output the total kinetic energy density KDEN
 767    if (dtset%prtkden/=0) then
 768      call fftdatar_write("kinedr",dtfil%fnameabo_app_kden,dtset%iomode,hdr,&
 769      crystal,ngfft,cplex1,nfft,nspden,taur,mpi_enreg,ebands=ebands)
 770    end if
 771 
 772 !  We output the Laplacian of density
 773    if (dtset%prtlden/=0) then
 774      call fftdatar_write("laprhor",dtfil%fnameabo_app_lden,dtset%iomode,hdr,&
 775      crystal,ngfft,cplex1,nfft,nspden,lrhor,mpi_enreg,ebands=ebands)
 776    end if
 777 
 778    call timab(954,2,tsec)
 779    call timab(955,1,tsec)
 780 
 781    call timab(955,2,tsec)
 782    call timab(956,1,tsec)
 783 
 784 !  POT
 785    if (dtset%prtpot>0) then
 786      call fftdatar_write("vtrial",dtfil%fnameabo_app_pot,dtset%iomode,hdr,&
 787      crystal,ngfft,cplex1,nfft,nspden,vtrial,mpi_enreg,ebands=ebands)
 788    end if
 789 
 790    call timab(956,2,tsec)
 791    call timab(957,1,tsec)
 792 
 793    if (dtset%prtgeo>0) then
 794      coordn=dtset%prtgeo
 795      call bonds_lgth_angles(coordn,dtfil%fnameabo_app_geo,natom,psps%ntypat,&
 796 &     rprimd,dtset%typat,xred,dtset%znucl)
 797    end if
 798 
 799    if (dtset%prtcif > 0) then
 800      call prt_cif(dtset%brvltt, dtfil%fnameabo_app_cif, natom, dtset%nsym, dtset%ntypat, rprimd, &
 801 &     dtset%spgaxor, dtset%spgroup, dtset%spgorig, dtset%symrel, dtset%tnons, dtset%typat, xred, dtset%znucl)
 802    end if
 803 
 804    call timab(957,2,tsec)
 805    call timab(958,1,tsec)
 806 
 807 !  STM
 808    if (dtset%prtstm>0) then
 809      call fftdatar_write("stm",dtfil%fnameabo_app_stm,dtset%iomode,hdr,&
 810      crystal,ngfft,cplex1,nfft,nspden,rhor,mpi_enreg,ebands=ebands)
 811    end if
 812 
 813    if (dtset%prt1dm>0) then
 814      call out1dm(dtfil%fnameabo_app_1dm,mpi_enreg,natom,nfft,ngfft,nspden,psps%ntypat,&
 815 &     rhor,rprimd,dtset%typat,ucvol,vtrial,xred,dtset%znucl)
 816    end if
 817 
 818 !  VHA
 819    if (dtset%prtvha>0) then
 820      ABI_ALLOCATE(vwork,(nfft,nspden))
 821      do ispden=1,nspden
 822        vwork(:,ispden)=vhartr(:)
 823      end do
 824 
 825      call fftdatar_write("vhartree",dtfil%fnameabo_app_vha,dtset%iomode,hdr,&
 826      crystal,ngfft,cplex1,nfft,nspden,vwork,mpi_enreg,ebands=ebands)
 827 
 828      ABI_DEALLOCATE(vwork)
 829    end if
 830 
 831 !  VPSP
 832    if (dtset%prtvpsp>0) then
 833      ABI_ALLOCATE(vwork,(nfft,nspden))
 834      do ispden=1,nspden
 835        vwork(:,ispden)=vpsp(:)
 836      end do
 837 
 838      call fftdatar_write("vpsp",dtfil%fnameabo_app_vpsp,dtset%iomode,hdr,&
 839      crystal,ngfft,cplex1,nfft,nspden,vwork,mpi_enreg,ebands=ebands)
 840 
 841      ABI_DEALLOCATE(vwork)
 842    end if
 843 
 844 ! VCouLoMB
 845    if (dtset%prtvclmb>0) then
 846 
 847      ABI_ALLOCATE(vwork,(nfft,nspden))
 848      do ispden=1,nspden
 849        vwork(:,ispden)=vpsp(:)+vhartr(:)
 850      end do
 851      if (psps%usepaw==1) then
 852        do ispden=1,nspden
 853          vwork(:,ispden)=vwork(:,ispden)+vpaw(:,ispden)
 854        end do
 855        ABI_DEALLOCATE(vpaw)
 856      end if
 857 
 858      call fftdatar_write("vhartree_vloc",dtfil%fnameabo_app_vclmb,dtset%iomode,hdr,&
 859 &     crystal,ngfft,cplex1,nfft,nspden,vwork,mpi_enreg,ebands=ebands)
 860 
 861 !TODO: find out why this combination of calls with fftdatar_write then out1dm fails on buda with 4 mpi-fft procs (npkpt 1).
 862 !      For the moment comment it out. Only DS2 of mpiio test 27 fails
 863 !     call out1dm(dtfil%fnameabo_app_vclmb_1dm,mpi_enreg,natom,nfft,ngfft,nspden,psps%ntypat,&
 864 !&         rhor,rprimd,dtset%typat,ucvol,vwork,xred,dtset%znucl)
 865 
 866 ! TODO: add TEM phase with CE = (2 pi / lambda) (E+E0)/(E(E+2E0)) from p.49 of RE Dunin Borkowski 2004 encyclopedia of nanoscience volume 3 pp 41-99
 867 !   where E is energy of electron, E0 rest mass, lambda the relativistic wavelength
 868 !   values of CE at 200 300 and 1000 kV:  7.29e6  6.53e6   5.39e6 rad / V / m
 869 !   vertical integral of vclmb * c / ngfft(3) / cross sectional area factor (= sin(gamma))
 870 !      * 0.5291772083e-10*27.2113834 to get to SI
 871 !      * CE factor above
 872 !   should be done for each plane perpendicular to the axes...
 873      ABI_DEALLOCATE(vwork)
 874    end if ! prtvclmb
 875 
 876 
 877 !  VHXC
 878    if (dtset%prtvhxc>0) then
 879      ABI_ALLOCATE(vwork,(nfft,nspden))
 880      do ispden=1,nspden
 881        vwork(:,ispden)=vhartr(:)+vxc(:,ispden)
 882      end do
 883 
 884      call fftdatar_write("vhxc",dtfil%fnameabo_app_vhxc,dtset%iomode,hdr,&
 885 &     crystal,ngfft,cplex1,nfft,nspden,vwork,mpi_enreg,ebands=ebands)
 886 
 887      ABI_DEALLOCATE(vwork)
 888    end if
 889 
 890 !  VXC
 891    if (dtset%prtvxc>0) then
 892      call fftdatar_write("exchange_correlation_potential",dtfil%fnameabo_app_vxc,dtset%iomode,hdr,&
 893      crystal,ngfft,cplex1,nfft,nspden,vxc,mpi_enreg,ebands=ebands)
 894    end if
 895 
 896    call timab(958,2,tsec)
 897  end if ! if iwrite_fftdatar
 898 
 899  call timab(959,1,tsec)
 900 
 901 !Generate DOS using the tetrahedron method or using Gaussians
 902 !FIXME: Should centralize all calculations of DOS here in outscfcv
 903  if (dtset%prtdos>=2.or.dtset%pawfatbnd>0) then
 904    dos = epjdos_new(dtset, psps, pawtab)
 905 
 906    if (dos%partial_dos_flag>=1 .or. dos%fatbands_flag==1)then
 907      ! Generate fractions for partial DOSs if needed partial_dos 1,2,3,4  give different decompositions
 908      collect = 1 !; if (psps%usepaw==1 .and. dos%partial_dos_flag /= 2) collect = 0
 909      if ((psps%usepaw==0.or.dtset%pawprtdos/=2) .and. dos%partial_dos_flag>=1) then
 910        call partial_dos_fractions(dos,crystal,dtset,eigen,occ,npwarr,kg,cg,mcg,collect,mpi_enreg)
 911      end if
 912 
 913      if (psps%usepaw==1 .and. dos%partial_dos_flag /= 2) then
 914 !      TODO: update partial_dos_fractions_paw for extra atoms - no PAW contribution normally, but check bounds and so on.
 915        call partial_dos_fractions_paw(dos,cprj,dimcprj,dtset,mcprj,mkmem,mpi_enreg,pawrad,pawtab)
 916      end if
 917 
 918    else
 919      dos%fractions(:,:,:,1)=one
 920    end if
 921 
 922 !  Here, print out fatbands for the k-points given in file appended _FATBANDS
 923    if (me == master .and. dtset%pawfatbnd>0 .and. dos%fatbands_flag==1) then
 924      call prtfatbands(dos,dtset,ebands,dtfil%fnameabo_app_fatbands,dtset%pawfatbnd,pawtab)
 925    end if
 926 
 927 !  Here, computation and output of DOS and partial DOS  _DOS
 928    if (dos%fatbands_flag == 0) then
 929      if (dos%prtdos /= 4) then
 930        call dos_calcnwrite(dos,dtset,crystal,ebands,dtfil%fnameabo_app_dos,spacecomm)
 931      end if
 932    end if
 933 
 934 #ifdef HAVE_NETCDF
 935    ! Write netcdf file with dos% results.
 936    if (me == master) then
 937      fname = trim(dtfil%filnam_ds(4))//'_FATBANDS.nc'
 938      NCF_CHECK(nctk_open_create(ncid, fname, xmpi_comm_self))
 939      call fatbands_ncwrite(dos, crystal, ebands, hdr, dtset, psps, pawtab, ncid)
 940      NCF_CHECK(nf90_close(ncid))
 941    end if
 942 #endif
 943 
 944    call epjdos_free(dos)
 945  end if ! prtdos > 1
 946 
 947  call timab(959,2,tsec)
 948  call timab(960,1,tsec)
 949 
 950 !Output of integrated density inside atomic spheres
 951  if (dtset%prtdensph==1.and.dtset%usewvl==0)then
 952    call calcdensph(gmet,mpi_enreg,natom,nfft,ngfft,nspden,&
 953 &   ntypat,ab_out,dtset%ratsph,rhor,rprimd,dtset%typat,ucvol,xred,1,cplex1)
 954  end if
 955 
 956  call timab(960,2,tsec)
 957 
 958  if (dtset%magconon /= 0) then
 959 !  calculate final value of terms for magnetic constraint: "energy" term, lagrange multiplier term, and atomic contributions
 960    call mag_constr_e(dtset%magconon,dtset%magcon_lambda,mpi_enreg,&
 961 &   natom,nfft,ngfft,nspden,ntypat,dtset%ratsph,rhor,rprimd,dtset%spinat,dtset%typat,xred)
 962  end if
 963 
 964  call timab(961,1,tsec)
 965 
 966 !If PAW, provide additional outputs
 967  if (psps%usepaw==1) then
 968 !  Output of compensation charge
 969    if (dtset%nstep>0.or.dtfil%ireadwf/=0) then
 970      write(message, '(4a)' )ch10,' PAW TEST:',ch10,&
 971 &     ' ==== Compensation charge inside spheres ============'
 972      if (compch_sph>-1.d4.and.compch_fft>-1.d4) &
 973 &     write(message, '(3a)' ) trim(message),ch10,' The following values must be close to each other ...'
 974      if (compch_sph>-1.d4) write(message, '(3a,f22.15)' ) trim(message),ch10,&
 975 &     ' Compensation charge over spherical meshes = ',compch_sph
 976      if (compch_fft>-1.d4) then
 977        if (pawfgr%usefinegrid==1) then
 978          write(message, '(3a,f22.15)' ) trim(message),ch10,&
 979 &         ' Compensation charge over fine fft grid    = ',compch_fft
 980        else
 981          write(message, '(3a,f22.15)' ) trim(message),ch10,&
 982 &         ' Compensation charge over fft grid         = ',compch_fft
 983        end if
 984      end if
 985      call wrtout(ab_out,message,'COLL')
 986      call wrtout(std_out,message,'COLL')
 987    end if
 988 !  Output of pseudopotential strength Dij and augmentation occupancies Rhoij
 989    call pawprt(dtset,my_natom,paw_ij,pawrhoij,pawtab,&
 990 &   comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,&
 991 &   electronpositron=electronpositron)
 992  end if
 993 
 994  call timab(961,2,tsec)
 995  call timab(962,1,tsec)
 996 
 997 
 998 !PAW + output for optical conductivity   _OPT and _OPT2
 999  if (psps%usepaw==1.and.prtnabla>0) then
1000    if (prtnabla==1.or.prtnabla==2) then
1001      call optics_paw(atindx1,cg,cprj,dimcprj,dtfil,dtset,eigen,gprimd,hdr,kg,&
1002 &     mband,mcg,mcprj,mkmem,mpi_enreg,mpsang,mpw,natom,nkpt,npwarr,nsppol,pawrad,pawtab)
1003    end if
1004    if (prtnabla==2.or.prtnabla==3) then
1005      call optics_paw_core(atindx1,cprj,dimcprj,dtfil,dtset,eigen,psps%filpsp,hdr,&
1006 &     mband,mcprj,mkmem,mpi_enreg,mpsang,natom,nkpt,nsppol,pawrad,pawtab)
1007    end if
1008  end if
1009  if (prtnabla<0) then
1010    ! TODO: This routine is not tested but it's used in production.
1011    call optics_vloc(cg,dtfil,dtset,eigen,gprimd,hdr,kg,&
1012 &   mband,mcg,mkmem,mpi_enreg,mpw,nkpt,npwarr,nsppol)
1013  end if
1014 
1015  call timab(962,2,tsec)
1016  call timab(963,1,tsec)
1017 
1018 !Optionally provide output for AE wavefunctions (only for PAW)
1019  if (psps%usepaw==1 .and. dtset%pawprtwf==1) then
1020    ABI_ALLOCATE(ps_norms,(nsppol,nkpt,mband))
1021 
1022    call pawmkaewf(dtset,crystal,ebands,my_natom,mpw,mband,mcg,mcprj,nkpt,mkmem,nsppol,Dtset%nband,&
1023 &   Dtset%istwfk,npwarr,Dtset%kptns,Dtset%ngfftdg,kg,dimcprj,pawfgrtab,&
1024 &   Pawrad,Pawtab,Hdr,Dtfil,cg,Cprj,&
1025 &   MPI_enreg,ierr,pseudo_norms=ps_norms,set_k=dtset%pawprt_k,set_band=dtset%pawprt_b,&
1026 &   comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
1027 
1028    if (dtset%pawprt_b==0) then
1029      fname = strcat(dtfil%filnam_ds(4), '_PAWSTAT')
1030      if (open_file(fname, message,newunit=tmp_unt,status='unknown',form='formatted') /= 0) then
1031        MSG_ERROR(message)
1032      end if
1033      write(tmp_unt,'(5a)') '# This file contains the statistics on the cancellation of',ch10,&
1034 &     '# the onsite pseudo component of the all-electron wavefunction',ch10,&
1035 &     '# with the plane wave part'
1036      ii = 0
1037      do isppol=1,nsppol
1038        write(tmp_unt,'(a,i0)') '# isppol = ',isppol
1039        do ikpt=1,nkpt
1040          write(tmp_unt,'(a,i0)') '# ikpt = ',ikpt
1041          write(tmp_unt,'(a)') '#    band      norm'
1042          occ_norm = zero; unocc_norm = zero; nocc = 0
1043          do iband=1,dtset%nband(ikpt + (isppol-1)*nkpt)
1044            ii = ii + 1
1045            write(tmp_unt,'(i8,ES16.6)') iband,ps_norms(isppol,ikpt,iband)
1046            if (abs(occ(ii)) <= tol16) then
1047              unocc_norm = unocc_norm + ps_norms(isppol,ikpt,iband)
1048            else
1049              occ_norm = occ_norm + ps_norms(isppol,ikpt,iband)
1050              nocc = nocc + 1
1051            end if
1052          end do
1053          if(mband/=nocc)then
1054            write(tmp_unt,'(2(a,ES16.6))') '# occ average: ',occ_norm/real(nocc),&
1055 &           ' unocc average: ',unocc_norm/real(mband-nocc)
1056          else
1057            write(tmp_unt,'(2(a,ES16.6))') '# occ average: ',occ_norm/real(nocc)
1058          end if
1059        end do
1060      end do
1061      close(tmp_unt)
1062    end if
1063    ABI_DEALLOCATE(ps_norms)
1064  end if
1065 
1066  call timab(963,2,tsec)
1067  if(dtset%plowan_compute>0) then
1068    write(message,'(2a,i3)') ch10,&
1069 &   ' ====================================================================================== '
1070    call wrtout(std_out,message,'COLL')
1071    call wrtout(ab_out,message,'COLL')
1072    write(message,'(2a,i3)') ch10,&
1073 &   ' == Start computation of Projected Local Orbitals Wannier functions == ',dtset%nbandkss
1074    call wrtout(std_out,message,'COLL')
1075    call wrtout(ab_out,message,'COLL')
1076 
1077 !  ==  compute psichi
1078 
1079    call init_plowannier(dtset,wan)
1080    call compute_coeff_plowannier(crystal,cprj,dimcprj,dtset,eigen,e_fermie,&
1081 &   mpi_enreg,occ,wan,pawtab,psps,usecprj,dtfil%unpaw,pawrad,dtfil)
1082    call destroy_plowannier(wan)
1083  end if
1084 
1085 !Optionally provide output for the GW part of ABINIT
1086  if (dtset%nbandkss/=0) then
1087    ! Use DMFT to compute wannier function for cRPA calculation.
1088    if(dtset%usedmft==1) then
1089      write(message,'(2a,i3)') ch10,&
1090 &     '  Warning: Psichi are renormalized in datafordmft because nbandkss is used',dtset%nbandkss
1091      call wrtout(std_out,message,'COLL')
1092      call init_dmft(dmatpawu,dtset,e_fermie,dtfil%fnameabo_app,dtset%nspinor,paw_dmft,pawtab,psps,dtset%typat)
1093      call print_dmft(paw_dmft,dtset%pawprtvol)
1094 
1095 !    ==  compute psichi
1096      call init_oper(paw_dmft,lda_occup)
1097 
1098      call datafordmft(crystal,cprj,dimcprj,dtset,eigen,e_fermie,&
1099 &     lda_occup,dtset%mband,dtset%mband,dtset%mkmem,mpi_enreg,&
1100 &     dtset%nkpt,dtset%nspinor,dtset%nsppol,occ,&
1101 &     paw_dmft,paw_ij,pawang,pawtab,psps,usecprj,dtfil%unpaw,dtset%nbandkss)
1102 
1103      call destroy_dmft(paw_dmft)
1104      call destroy_oper(lda_occup)
1105    end if
1106 
1107    call timab(964,1,tsec) ! outscfcv(outkss)
1108 
1109    call outkss(crystal,dtfil,dtset,ecut,gmet,gprimd,hdr,&
1110 &   dtset%kssform,mband,mcg,mcprj,mgfftc,mkmem,mpi_enreg,mpsang,mpw,natom,natom,&
1111 &   nfft,nkpt,npwarr,nspden,nsppol,nsym,psps%ntypat,occ,pawtab,pawfgr,paw_ij,&
1112 &   prtvol,psps,rprimd,vtrial,xred,cg,usecprj,cprj,eigen,ierr)
1113    call timab(964,2,tsec) ! outscfcv(outkss)
1114    if (ierr/=0) then
1115      MSG_WARNING("outkss returned a non zero status error, check log")
1116    end if
1117  end if
1118 
1119  if (electronpositron_calctype(electronpositron)/=0) then
1120 
1121 !  Optionally provide output for  positron life time calculation
1122    call timab(965,1,tsec)
1123    call poslifetime(dtset,electronpositron,gprimd,my_natom,&
1124 &   mpi_enreg,n3xccc,nfft,ngfft,nhat,1,pawang,&
1125 &   pawrad,pawrhoij,pawtab,rate_dum,rate_dum2,&
1126 &   rhor,ucvol,xccc3d)
1127    call timab(965,2,tsec)
1128 
1129 !  Optionally provide output for momentum distribution of annihilation radiation
1130    if (dtset%posdoppler>0) then
1131      call posdoppler(cg,cprj,crystal,dimcprj,dtfil,dtset,electronpositron,psps%filpsp,&
1132 &     kg,mcg,mcprj,mpi_enreg,my_natom,n3xccc,nfft,ngfft,nhat,npwarr,&
1133 &     occ,pawang,pawrad,pawrhoij,pawtab,rhor,xccc3d)
1134    end if
1135  end if
1136 
1137 !Optionally provide output for WanT
1138  if (dtset%prtwant==1) then
1139    call timab(966,1,tsec)
1140    ! WARNING: mpi_enreg not used --> MPI is not supported
1141    call outwant(dtset,eigen,cg,kg,npwarr,mband,mcg,nkpt,nsppol,mkmem,mpw,dtset%prtwant)
1142    call timab(966,2,tsec)
1143  end if
1144 
1145 !Optionally provide output for electric field gradient calculation
1146  if (dtset%prtefg > 0) then
1147    call timab(967,1,tsec)
1148    call calc_efg(mpi_enreg,my_natom,natom,nfft,ngfft,nspden,dtset%nsym,ntypat,&
1149 &   dtset%paral_kgb,paw_an,pawang,pawrad,pawrhoij,pawtab,&
1150 &   dtset%ptcharge,dtset%prtefg,dtset%quadmom,rhor,rprimd,dtset%symrel,&
1151 &   dtset%tnons,dtset%typat,ucvol,psps%usepaw,xred,psps%zionpsp,&
1152 &   comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
1153    call timab(967,2,tsec)
1154  end if
1155 
1156 !Optionally provide output for Fermi-contact term at nuclear positions
1157  if (dtset%prtfc > 0) then
1158    call timab(967,1,tsec)
1159    call calc_fc(my_natom,natom,nspden,ntypat,pawrad,pawrhoij,pawtab,dtset%typat,psps%usepaw,&
1160 &   comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
1161    call timab(967,2,tsec)
1162  end if
1163 
1164  ! Output electron bands.
1165  if (me == master .and. dtset%tfkinfunc==0) then
1166    if (size(dtset%kptbounds, dim=2) > 0) then
1167      call ebands_write(ebands, dtset%prtebands, dtfil%filnam_ds(4), kptbounds=dtset%kptbounds)
1168    else
1169      call ebands_write(ebands, dtset%prtebands, dtfil%filnam_ds(4))
1170    end if
1171  end if
1172 
1173 !Optionally provide Xcrysden output for the Fermi surface (Only master writes)
1174  if (me == master .and. dtset%prtfsurf == 1) then
1175    if (ebands_write_bxsf(ebands,crystal,dtfil%fnameabo_app_bxsf) /= 0) then
1176      message = "Cannot produce BXSF file with Fermi surface, see log file for more info"
1177      MSG_WARNING(message)
1178      call wrtout(ab_out, message)
1179    end if
1180  end if ! prtfsurf==1
1181 
1182 !output nesting factor for Fermi surface (requires ph_nqpath)
1183  if (me == master .and. dtset%prtnest>0 .and. dtset%ph_nqpath > 0) then
1184    call timab(968,1,tsec)
1185    ierr = ebands_write_nesting(ebands,crystal,dtfil%fnameabo_app_nesting,dtset%prtnest,&
1186 &   dtset%tsmear,dtset%fermie_nest,dtset%ph_qpath(:,1:dtset%ph_nqpath),message)
1187    if (ierr /=0) then
1188      MSG_WARNING(message)
1189      call wrtout(ab_out,message,'COLL')
1190    end if
1191    call timab(968,2,tsec)
1192  end if ! prtnest=1
1193 
1194  call timab(969,1,tsec)
1195 
1196  if (dtset%prtdipole == 1) then
1197    call multipoles_out(rhor,mpi_enreg,natom,nfft,ngfft,dtset%nspden,dtset%ntypat,rprimd,&
1198 &   dtset%typat,ucvol,ab_out,xred,dtset%ziontypat)
1199  end if
1200 
1201  ! BoltzTraP output files in GENEric format
1202  if (dtset%prtbltztrp == 1 .and. me==master) then
1203    call ebands_prtbltztrp(ebands, crystal, dtfil%filnam_ds(4))
1204  end if
1205 
1206  ! Band structure interpolation from eigenvalues computed on the k-mesh.
1207  if (nint(dtset%einterp(1)) /= 0) then
1208    call ebands_interpolate_kpath(ebands, dtset, crystal, [0,0], dtfil%filnam_ds(4), spacecomm)
1209  end if
1210 
1211  call crystal_free(crystal)
1212  call ebands_free(ebands)
1213 
1214 !Destroy atom table used for parallelism
1215  call free_my_atmtab(my_atmtab,my_atmtab_allocated)
1216 
1217  call timab(969,2,tsec)
1218  call timab(950,2,tsec) ! outscfcv
1219 
1220  DBG_EXIT("COLL")
1221 
1222 end subroutine outscfcv