TABLE OF CONTENTS


ABINIT/dsksta [ Functions ]

[ Top ] [ Functions ]

NAME

 dsksta

FUNCTION

 This routine evaluates the amount of disk space required by the _KSS file.

INPUTS

  dimlmn(natom*usepaw)=Number of nlm partial waves for each atom.
  ishm=Number of G-shells to be saved in _KSS file.
  mpsang=Max angular momentum +1 for pseudos.
  natom=Number of atoms in the unit cell.
  nbandkss=Number of desired bands to be saved in _KSS file
  nkpt=Number of k points.
  npwkss=Number of desired G-vectors to be saved in _KSS file.
  nspinor=Number of spinorial components.
  nsppol=Number of independent spin polarizations.
  ntypat=Number of type of atoms.
  nsym2=Number of symmetries in space group, without INV
  usepaw=1 if PAW.

OUTPUT

  Writes on standard output

PARENTS

      outkss

CHILDREN

      wrtout

SOURCE

2704 subroutine dsksta(ishm,usepaw,nbandkss,mpsang,natom,ntypat,npwkss,nkpt,nspinor,nsppol,nsym2,dimlmn)
2705 
2706 
2707 !This section has been created automatically by the script Abilint (TD).
2708 !Do not modify the following lines by hand.
2709 #undef ABI_FUNC
2710 #define ABI_FUNC 'dsksta'
2711 !End of the abilint section
2712 
2713  implicit none
2714 
2715 !Arguments ------------------------------------
2716 !scalars
2717  integer,intent(in) :: usepaw,ishm,nbandkss,mpsang,natom,ntypat,nkpt
2718  integer,intent(in) :: npwkss,nspinor,nsppol,nsym2
2719 !arrays
2720  integer,intent(in) :: dimlmn(natom*usepaw)
2721 
2722 !Local variables-------------------------------
2723 !scalars
2724  integer :: bsize_tot,bsize_hdr,bsize_kb,bsize_wf,bsize_cprj
2725  character(len=500) :: msg
2726 
2727 ! *********************************************************************
2728 
2729 !The Abinit header is not considered.
2730  bsize_hdr= 80*2 + & !title
2731 &5*4 + & !nsym2,nbandksseff,npwkss,ishm,mpsang
2732 &nsym2*9*4 + & !symrel2
2733 &nsym2*3*8 + & !tnons
2734 &npwkss*3*4 + & !gbig
2735 &ishm*4     !shlim
2736 
2737 !NOTE: vkb does not depend on nsppol, however the elements are written for each spin.
2738  bsize_kb=0
2739  if (usepaw==0) then
2740    bsize_kb= nsppol* &
2741 &   (         mpsang*ntypat       *8 + & !vkbsign
2742 &  2*(nkpt*mpsang*ntypat*npwkss*8)  & !vkbd,vkbd
2743 &  )
2744  end if
2745 
2746  bsize_wf= nsppol* &
2747 & ( nkpt*nbandkss                 *8 + & !energies
2748 &nkpt*nbandkss*nspinor*npwkss*2*8   & !wfg
2749 &)
2750 
2751 !For PAW add space required by projectors.
2752  bsize_cprj=0
2753  if (usepaw==1) then
2754    bsize_cprj=SUM(dimlmn(:))*(nsppol*nkpt*nspinor*nbandkss*2*8)
2755  end if
2756 
2757  bsize_tot = bsize_hdr + bsize_kb + bsize_wf + bsize_cprj
2758  write(msg,'(2a,f8.2,4a,4(a,f8.2,2a))')ch10,&
2759 & ' Total amount of disk space required by _KSS file = ',bsize_tot*b2Mb,' Mb.',ch10,&
2760 & '  Subdivided into : ',ch10,&
2761 & '  Header             = ',bsize_hdr *b2Mb,' Mb.',ch10,&
2762 & '  KB elements        = ',bsize_kb  *b2Mb,' Mb.',ch10,&
2763 & '  Wavefunctions (PW) = ',bsize_wf  *b2Mb,' Mb.',ch10,&
2764 & '  PAW projectors     = ',bsize_cprj*b2Mb,' Mb.',ch10
2765  call wrtout(std_out,msg,'COLL')
2766 
2767 end subroutine dsksta

ABINIT/m_io_kss [ Modules ]

[ Top ] [ Modules ]

NAME

  m_io_kss

FUNCTION

  This module contains procedured dealing with the IO of the KSS file.

COPYRIGHT

 Copyright (C) 1999-2018 ABINIT group (MG, MT, VO, AR, LR, RWG, MM, XG, RShaltaf)
 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 .

PARENTS

CHILDREN

SOURCE

22 #if defined HAVE_CONFIG_H
23 #include "config.h"
24 #endif
25 
26 #include "abi_common.h"
27 
28 MODULE m_io_kss
29 
30  use defs_basis
31  use defs_datatypes
32  use defs_abitypes
33  use m_abicore
34  use m_xmpi
35  use m_errors
36  use m_nctk
37 #ifdef HAVE_NETCDF
38  use netcdf
39 #endif
40  use m_hdr
41  use m_wfk
42  use m_cgtools
43  use m_hamiltonian
44  use m_electronpositron
45  use m_pawtab
46  use m_paw_ij
47  use m_pawcprj
48  use m_pawfgr
49 
50  use m_time,             only : timab
51  use m_io_tools,         only : open_file
52  use m_fstrings,         only : sjoin, itoa, strcat
53  use m_hide_lapack,      only : xheevx, xhegvx
54  use m_geometry,         only : metric, remove_inversion
55  use m_dtset,            only : dtset_copy, dtset_free
56  use m_mpinfo,           only : destroy_mpi_enreg, proc_distrb_cycle
57  use m_fftcore,          only : get_kg, sphere
58  use m_fft,              only : fftpac
59  use m_crystal ,         only : crystal_t
60  use m_crystal_io,       only : crystal_ncwrite
61  use m_gsphere,          only : table_gbig2kg, merge_and_sort_kg
62  use m_kg,               only : mkkin, mkkpg
63  use m_ksdiago,          only : ksdiago, init_ddiago_ctl, ddiago_ctl_type
64  use m_mkffnl,           only : mkffnl
65  use m_getghc,           only : getghc
66  use m_fourier_interpol, only : transgrid
67 
68  implicit none
69 
70  private
71 
72  public :: write_kss_header    ! Writes the header of the KSS file.
73  !private :: write_vkb         ! Writes the KB form factors and derivates on file for a single k-point.
74  public :: write_kss_wfgk      ! Write the Gamma-centered wavefunctions and energies on the KSS file for a single k-point.
75  public :: k2gamma_centered    ! Convert a set of wavefunctions from the k-centered to the gamma-centered basis set.
76  public :: make_gvec_kss       ! Build the list of G-vectors for the KSS file.
77  public :: gshgg_mkncwrite     ! Debugging tool used to build <G|H|G'> and dump data in netcdf format.
78  public :: outkss              ! Generate KSS file
79 
80 CONTAINS  !===========================================================

ABINIT/memkss [ Functions ]

[ Top ] [ Functions ]

NAME

 memkss

FUNCTION

 This routine evaluates the additional amount of memory required
 by routine 'outkss'.

INPUTS

  mband=maximum number of bands
  mgfft=maximum size of 1D FFTs
  mproj=maximum dimension for number of projection operators for each
   angular momentum for nonlocal pseudopotential
  mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
  mpw=maximum dimensioned size of npw.
  natom=number of atoms in cell.
  ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
  nkpt=number of k points.
  nspinor=number of spinorial components of the wavefunctions
  nsym=number of symmetries in space group
  ntypat=number of types of atoms in unit cell.

NOTES

 This routine is not available for paw calculations

PARENTS

      outkss

CHILDREN

      wrtout

SOURCE

2566 subroutine memkss(mband,mgfft,mproj,mpsang,mpw,natom,ngfft,nkpt,nspinor,nsym,ntypat)
2567 
2568 
2569 !This section has been created automatically by the script Abilint (TD).
2570 !Do not modify the following lines by hand.
2571 #undef ABI_FUNC
2572 #define ABI_FUNC 'memkss'
2573 !End of the abilint section
2574 
2575  implicit none
2576 
2577 !Arguments ------------------------------------
2578 !scalars
2579  integer,intent(in) :: mband,mgfft,mproj,mpsang,mpw,natom,nkpt,nspinor
2580  integer,intent(in) :: nsym,ntypat
2581 !arrays
2582  integer,intent(in) :: ngfft(18)
2583 
2584 !Local variables-------------------------------
2585 !scalars
2586  integer(i8b) :: isize,memsize
2587  character(len=500) :: msg
2588 
2589 ! *********************************************************************
2590 !
2591  isize=580+fnlen+4*(81+nkpt+9*nsym)+8*15    !non allocatable var.
2592  if(xmpi_paral==1)then
2593    isize=isize+4*4                           !kpt_distrb
2594  end if
2595  memsize=isize
2596  isize=isize+4*nkpt+12*mpw+20*nkpt*mpw      !nbasek,gbasek,cnormk,gcurr
2597  memsize=max(memsize,isize)
2598  if(xmpi_paral==1)then
2599    isize=isize+12*mpw*nkpt                   !ibuf1,ibuf2,rbuf1
2600    memsize=max(memsize,isize)
2601    isize=isize-12*mpw*nkpt                   !ibuf1,ibuf2,rbuf1
2602  end if
2603  isize=isize+40*mpw                         !gbase,cnorm
2604  memsize=max(memsize,isize)
2605  isize=isize-4*nkpt-20*mpw*nkpt             !nbasek,gbasek,cnormk
2606  isize=isize+4*mpw                          !insort
2607  memsize=max(memsize,isize)
2608  isize=isize-16*mpw                         !cnorm
2609  isize=isize+28*mpw+24*nsym                 !gbig,nshell,gshell
2610  memsize=max(memsize,isize)
2611  isize=isize+4*mpw                          !shlim
2612  memsize=max(memsize,isize)
2613  isize=isize-44*mpw-24*nsym                 !gcurr,gbase,gshell,insort,nshell
2614  isize=isize-4*mpw                          !shlim
2615  isize=isize+8*mpw*nspinor&
2616 & +16*mpw*nspinor*(mpw*nspinor+1)       !eigval,eigvec
2617  memsize=max(memsize,isize)
2618  isize=isize+8*mpw+8*ngfft(4)&
2619 & *ngfft(5)*ngfft(6)      !ts,vlocal
2620  memsize=max(memsize,isize)
2621  isize=isize+8*mgfft+4+28*mpw               !gbound,indpw_k,kg_k
2622  memsize=max(memsize,isize)
2623  isize=isize+8*natom&
2624 & +24*mpw*ntypat*mpsang*mproj      !phkxred,ffnl,kinpw
2625  memsize=max(memsize,isize)
2626  isize=isize+16*mpw*natom                   !ph3d
2627  memsize=max(memsize,isize)
2628  isize=isize+48*mpw*nspinor&
2629 & +8*mpw*nspinor*(mpw*nspinor+1)        !pwave,subghg,gvnlg
2630  if (nspinor==2)&
2631 & isize=isize+40*mpw*nspinor                !pwave_so,subghg_so
2632  memsize=max(memsize,isize)
2633  isize=isize+8*mpw*nspinor*(mpw*nspinor+1)  !ghg
2634  memsize=max(memsize,isize)
2635  isize=isize+8*ngfft(4)*ngfft(5)*ngfft(6)   !work
2636  memsize=max(memsize,isize)
2637  isize=isize-8*ngfft(4)*ngfft(5)*ngfft(6)   !work
2638  isize=isize-8*mgfft+4+28*mpw&              !gbound,indpw_k,kg_k
2639 &-8*natom-24*mpw*ntypat*mpsang*mproj&  !phkxred,ffnl,kinpw
2640 &-16*mpw*natom                        !ph3d
2641  isize=isize-48*mpw*nspinor&
2642 & -8*mpw*nspinor*(mpw*nspinor+1)        !pwave,subghg,gvnlg
2643  if (nspinor==2)&
2644 & isize=isize-40*mpw*nspinor                !pwave_so,subghg_so
2645 
2646  isize=isize+56*mpw*nspinor                 !cwork,rwork
2647  memsize=max(memsize,isize)
2648  isize=isize-56*mpw*nspinor                 !cwork,rwork
2649  isize=isize+112*mpw*nspinor                !cwork,rwork,iwork,ifail
2650  memsize=max(memsize,isize)
2651  isize=isize-112*mpw*nspinor                !cwork,rwork,iwork,ifail
2652  isize=isize-8*mpw*nspinor*(mpw*nspinor+1)  !ghg
2653  isize=isize+8*mband                        !occ_k
2654  memsize=max(memsize,isize)
2655  isize=isize-8*mband                        !occ_k
2656  isize=isize-8*mpw*nspinor&
2657 & -16*mpw*nspinor*(mpw*nspinor+1)       !eigval,eigvec
2658  isize=isize-32*mpw-8*ngfft(4)&
2659 & *ngfft(5)*ngfft(6)     !gbig,ts,vlocal
2660  if(xmpi_paral==1)then
2661    isize=isize-4*4                           !kpt_distrb
2662  end if
2663  isize=isize-580-fnlen-4*(81+nkpt+9*nsym)-8*15   !non allocatable var.
2664 !
2665  write(msg,'(2a,f8.2,a)')ch10,&
2666 & ' Additional amount of memory required by "outkss" routine=',memsize*b2Mb,' Mbytes.'
2667  call wrtout(std_out,msg,'COLL')
2668 !
2669 end subroutine memkss

m_io_kss/gshgg_mkncwrite [ Functions ]

[ Top ] [ m_io_kss ] [ Functions ]

NAME

 gshgg_mkncwrite

FUNCTION

  This routine builds <G|H|G'> matrix elements for all (G, G').
  starting from the knowledge of the local potential on the real-space FFT mesh.
  It can also perform the direct diagonalization of the Kohn-Sham Hamiltonian
  for a given k-point and spin. This a debugging tool

INPUTS

  kg(3,mpw*mkmem)=reduced planewave coordinates.
  ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k point
  kpoint(3)
  natom=number of atoms in cell.
  nfftf=(effective) number of FFT grid points in the dense FFT mesh (for this processor)
         (nfftf=nfft for norm-conserving potential runs)
  nspinor=number of spinorial components of the wavefunctions
  ntypat=number of types of atoms in unit cell.
  pawtab(psps%ntypat*psps%usepaw) <type(pawtab_type)>=paw tabulated starting data
  pawfgr<pawfgr_type>=fine grid parameters and related data
  Paw_ij(natom) <type(paw_ij_type)>=paw arrays given on (i,j) channels
  psps <type(pseudopotential_type)>=variables related to pseudopotentials
  rprimd(3,3)=dimensional primitive translations for real space (bohr)
  vtrial(nfftf,nspden)=the trial potential
  xred(3,natom)=reduced dimensionless atomic coordinates
  eigen(mband,nkpt,nsppol)=array for holding eigenvalues (hartree)
  ngfftc(18)=Info about 3D FFT for the coarse mesh, see ~abinit/doc/variables/vargs.htm#ngfft
  [Electronpositron] <electronpositron_type>=quantities for the electron-positron annihilation.

PARENTS

      scfcv

CHILDREN

      metric,mkffnl,mkkin

SOURCE

1021 subroutine gshgg_mkncwrite(istep, dtset, dtfil, psps, hdr, pawtab, pawfgr, paw_ij, mpi_enreg, &
1022   rprimd, xred, eigen, npwarr, kg, ylm, ngfftc, nfftc, ngfftf, nfftf, vtrial,&
1023   electronpositron) ! Optional arguments
1024 
1025 
1026 !This section has been created automatically by the script Abilint (TD).
1027 !Do not modify the following lines by hand.
1028 #undef ABI_FUNC
1029 #define ABI_FUNC 'gshgg_mkncwrite'
1030 !End of the abilint section
1031 
1032  implicit none
1033 
1034 !Arguments ------------------------------------
1035 !scalars
1036  integer,intent(in) :: istep,nfftc,nfftf
1037  type(dataset_type),intent(in) :: dtset
1038  type(datafiles_type),intent(in) :: dtfil
1039  type(pseudopotential_type),intent(in) :: psps
1040  type(pawfgr_type),intent(in) :: pawfgr
1041  type(hdr_type),intent(in) :: hdr
1042  type(MPI_type),intent(inout) :: mpi_enreg
1043 !arrays
1044  integer,intent(in) :: ngfftc(18),ngfftf(18)
1045  integer, intent(in) :: kg(3,dtset%mpw*dtset%mkmem),npwarr(dtset%nkpt)
1046  real(dp),intent(in) :: rprimd(3,3)
1047  real(dp),intent(in) :: eigen(dtset%mband,dtset%nkpt,dtset%nsppol)
1048  real(dp),intent(inout) :: vtrial(nfftf,dtset%nspden)
1049  real(dp),intent(in) :: xred(3,dtset%natom)
1050  real(dp), intent(in) :: ylm(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm)
1051  type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw)
1052  type(paw_ij_type),intent(in) :: Paw_ij(dtset%natom*psps%usepaw)
1053  type(electronpositron_type),optional,pointer :: electronpositron
1054 
1055 !Local variables-------------------------------
1056 !scalars
1057  integer,parameter :: tim_getghc=4,master=0,fform=102
1058  integer :: comm,cpopt,dimffnl,ib,ider,idir,isppol,npw_k,npws,ierr
1059  integer :: ikg,istwf_k,ntypat,paral_kgb,ikpt,ilm,natom,nband_k,ncerr
1060  integer :: jj,n1,n2,n3,n4,n5,n6,negv,nkpg,nprocs,nspden,nspinor,mgfftc,ncid
1061  integer :: my_rank,ispden,ndat,type_calc,sij_opt,igsp2,cplex_ghg
1062  integer :: kg_varid,ghg_varid,pawgtg_varid,eighist_varid,eigdiago_varid
1063  real(dp) :: cfact,ucvol,lambda,size_mat
1064  character(len=50) :: jobz,range
1065  character(len=80) :: frmt1,frmt2
1066  character(len=10) :: stag(2)
1067  character(len=500) :: msg
1068  character(len=fnlen) :: path
1069  type(gs_hamiltonian_type) :: gs_hamk
1070 !arrays
1071  integer,allocatable :: kg_k(:,:)
1072  real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),rhodum(1),kpoint(3),ylmgr_dum(1,1,1)
1073  real(dp),allocatable :: eig_ene(:),eig_vec(:,:,:),ph3d(:,:,:),pwave(:,:)
1074  real(dp),allocatable :: ffnl(:,:,:,:),kinpw(:),kpg_k(:,:)
1075  real(dp),allocatable :: vlocal(:,:,:,:),ylm_k(:,:),vlocal_tmp(:,:,:)
1076  real(dp),allocatable :: ghc(:,:),gvnlc(:,:),gsc(:,:),ghg_mat(:,:,:),gtg_mat(:,:,:),cgrvtrial(:,:)
1077  type(pawcprj_type),allocatable :: cwaveprj(:,:)
1078 
1079 ! *********************************************************************
1080 
1081  DBG_ENTER("COLL")
1082  ABI_UNUSED(ngfftf(1))
1083 
1084  call wrtout(std_out, "Constructing H_{GG')(k,spin) and dumping matrices to file")
1085 
1086  comm = mpi_enreg%comm_cell; nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
1087  paral_kgb = mpi_enreg%paral_kgb
1088  ABI_CHECK(paral_kgb == 0 ,"paral_kgb == 1 not coded")
1089  ABI_CHECK(all(dtset%nband == dtset%mband) ,"nband must be constant")
1090 
1091  natom = dtset%natom; ntypat = dtset%ntypat
1092  nspden = dtset%nspden; nspinor = dtset%nspinor
1093  if (dtset%nsppol==1) stag=['          ','          ']
1094  if (dtset%nsppol==2) stag=['SPIN UP:  ','SPIN DOWN:']
1095 
1096  if (nprocs > 1 .and. .not. nctk_has_mpiio) then
1097    MSG_ERROR("Netcdf without MPI-IO support. Cannot produce HGG.nc file with nprocs > 1")
1098  end if
1099 
1100  path = strcat(dtfil%filnam_ds(4), "_HGG.nc")
1101 #ifdef HAVE_NETCDF
1102  if (istep == 1) then
1103    ! Open netcdf file, define dims and variables.
1104    NCF_CHECK(nctk_open_create(ncid, path, comm))
1105 
1106    ncerr = nctk_def_dims(ncid, [&
1107      nctkdim_t("complex", 2), nctkdim_t("max_number_of_coefficients", dtset%mpw),&
1108      nctkdim_t("max_number_of_states", dtset%mband),&
1109      nctkdim_t("number_of_kpoints", dtset%nkpt), nctkdim_t("number_of_spins", dtset%nsppol), &
1110      nctkdim_t("nstep", dtset%nstep) ])
1111    NCF_CHECK(ncerr)
1112 
1113    call wfk_ncdef_dims_vars(ncid, hdr, fform)
1114 
1115    NCF_CHECK(nctk_def_iscalars(ncid, ["last_step"]))
1116 
1117    ncerr = nctk_def_arrays(ncid, nctkarr_t('eigenvalues_history', "dp", &
1118      "max_number_of_states, number_of_kpoints, number_of_spins, nstep"))
1119    NCF_CHECK(ncerr)
1120 
1121    ncerr = nctk_def_arrays(ncid, nctkarr_t('eigenvalues_diago', "dp", &
1122      "max_number_of_states, number_of_kpoints, number_of_spins, nstep"))
1123    NCF_CHECK(ncerr)
1124 
1125    ncerr = nctk_def_arrays(ncid, nctkarr_t('ghg', "dp", &
1126      "complex, max_number_of_coefficients, max_number_of_coefficients, number_of_kpoints, number_of_spins, nstep"))
1127    NCF_CHECK(ncerr)
1128 
1129    if (psps%usepaw == 1) then
1130      ncerr = nctk_def_arrays(ncid, nctkarr_t('paw_gtg', "dp", &
1131        "complex, max_number_of_coefficients, max_number_of_coefficients, number_of_kpoints, number_of_spins, nstep"))
1132      NCF_CHECK(ncerr)
1133    end if
1134  else
1135    NCF_CHECK(nctk_open_modify(ncid, path, comm))
1136  end if
1137 
1138  ! Use individual IO (default) for the G-vectors [3, mpw, nkpt]
1139  NCF_CHECK(nctk_set_datamode(ncid))
1140  NCF_CHECK(nf90_inq_varid(ncid, "reduced_coordinates_of_plane_waves", kg_varid))
1141  NCF_CHECK(nf90_inq_varid(ncid, "eigenvalues_diago", eigdiago_varid))
1142  NCF_CHECK(nf90_inq_varid(ncid, "eigenvalues_history", eighist_varid))
1143  NCF_CHECK(nf90_inq_varid(ncid, "ghg", ghg_varid))
1144  if (psps%usepaw == 1) then
1145    NCF_CHECK(nf90_inq_varid(ncid, "paw_gtg", pawgtg_varid))
1146  end if
1147 
1148  NCF_CHECK(nctk_write_iscalars(ncid, ["last_step"], [istep]))
1149 #endif
1150 
1151  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
1152 
1153  ! The coarse FFT mesh.
1154  mgfftc = maxval(ngfftc(:3))
1155  n1=ngfftc(1); n2=ngfftc(2); n3=ngfftc(3)
1156  n4=ngfftc(4); n5=ngfftc(5); n6=ngfftc(6)
1157 
1158 !============================================
1159 !==== Initialize most of the Hamiltonian ====
1160 !============================================
1161 !1) Allocate all arrays and initialize quantities that do not depend on k and spin.
1162 !2) Perform the setup needed for the non-local factors:
1163 !* Norm-conserving: Constant kleimann-Bylander energies are copied from psps to gs_hamk.
1164 !* PAW: Initialize the overlap coefficients and allocate the Dij coefficients.
1165 
1166  if (PRESENT(Electronpositron)) then
1167    call init_hamiltonian(gs_hamk,psps,pawtab,nspinor,dtset%nsppol,nspden,natom,dtset%typat,xred,nfftc,&
1168 &   mgfftc,ngfftc,rprimd,dtset%nloalg,paw_ij=paw_ij,usecprj=0,electronpositron=electronpositron,&
1169 &   comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab)
1170  else
1171    call init_hamiltonian(gs_hamk,psps,pawtab,nspinor,dtset%nsppol,nspden,natom,dtset%typat,xred,nfftc,&
1172 &   mgfftc,ngfftc,rprimd,dtset%nloalg,paw_ij=paw_ij,usecprj=0,&
1173 &   comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab)
1174  end if
1175 
1176  ! Set up local potential vlocal with proper dimensioning, from vtrial.
1177  ! Select spin component of interest if nspden<=2 as nvloc==1, for nspden==4, nvloc==4
1178  ! option=2: vtrial(n1*n2*n3,ispden) --> vlocal(nd1,nd2,nd3) real case
1179  ABI_MALLOC(vlocal,(n4,n5,n6,gs_hamk%nvloc))
1180  ABI_MALLOC(kg_k,(3,dtset%mpw))
1181 
1182  do isppol=1,dtset%nsppol
1183    ikg = 0
1184 
1185    ! Set up local potential vlocal with proper dimensioning, from vtrial
1186    ! Also take into account the spin.
1187    if (nspden/=4)then
1188      if (psps%usepaw==0 .or. pawfgr%usefinegrid==0) then
1189        call fftpac(isppol,mpi_enreg,nspden,n1,n2,n3,n4,n5,n6,ngfftc,vtrial,vlocal,2)
1190      else
1191        ! Move from fine to coarse FFT mesh (PAW)
1192        ABI_MALLOC(cgrvtrial,(nfftc,nspden))
1193        call transgrid(1,mpi_enreg,nspden,-1,0,0,paral_kgb,pawfgr,rhodum,rhodum,cgrvtrial,vtrial)
1194        call fftpac(isppol,mpi_enreg,nspden,n1,n2,n3,n4,n5,n6,ngfftc,cgrvtrial,vlocal,2)
1195        ABI_FREE(cgrvtrial)
1196      end if
1197    else
1198      ABI_MALLOC(vlocal_tmp,(n4,n5,n6))
1199      if (psps%usepaw==0 .or. pawfgr%usefinegrid==0) then
1200        do ispden=1,nspden
1201          call fftpac(ispden,mpi_enreg,nspden,n1,n2,n3,n4,n5,n6,ngfftc,vtrial,vlocal_tmp,2)
1202          vlocal(:,:,:,ispden)=vlocal_tmp(:,:,:)
1203        end do
1204      else
1205        ! Move from fine to coarse FFT mesh (PAW)
1206        ABI_MALLOC(cgrvtrial,(nfftc,nspden))
1207        call transgrid(1,mpi_enreg,nspden,-1,0,0,paral_kgb,pawfgr,rhodum,rhodum,cgrvtrial,vtrial)
1208        do ispden=1,nspden
1209          call fftpac(ispden,mpi_enreg,nspden,n1,n2,n3,n4,n5,n6,ngfftc,cgrvtrial,vlocal_tmp,2)
1210          vlocal(:,:,:,ispden)=vlocal_tmp(:,:,:)
1211        end do
1212        ABI_FREE(cgrvtrial)
1213      end if
1214      ABI_FREE(vlocal_tmp)
1215    end if
1216 
1217    !Continue to initialize the Hamiltonian
1218    call load_spin_hamiltonian(gs_hamk,isppol,vlocal=vlocal,with_nonlocal=.true.)
1219 
1220    do ikpt=1,dtset%nkpt
1221      nband_k = dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
1222      istwf_k = dtset%istwfk(ikpt)
1223      npw_k = npwarr(ikpt)
1224      kpoint = dtset%kptns(:,ikpt)
1225 
1226      ! Skip the rest of the k-point loop
1227      if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,my_rank)) cycle
1228 
1229      ABI_MALLOC(ylm_k, (npw_k,psps%mpsang**2*psps%useylm))
1230 
1231      kg_k(:,1:npw_k) = kg(:,1+ikg:npw_k+ikg)
1232 
1233      if (psps%useylm==1) then
1234        do ilm=1,psps%mpsang**2
1235          ylm_k(1:npw_k,ilm) = ylm(1+ikg:npw_k+ikg,ilm)
1236        end do
1237      end if
1238 
1239      ! Set up remaining of the Hamiltonian
1240      ! Compute (1/2) (2 Pi)**2 (k+G)**2:
1241      ABI_ALLOCATE(kinpw,(npw_k))
1242      call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg_k,kinpw,kpoint,npw_k,0,0)
1243 
1244      ! Compute (k+G) vectors (only if useylm=1)
1245      nkpg=3*dtset%nloalg(3)
1246      ABI_ALLOCATE(kpg_k,(npw_k,nkpg))
1247      if (nkpg>0) then
1248        call mkkpg(kg_k,kpg_k,kpoint,nkpg,npw_k)
1249      end if
1250 
1251      ! Compute nonlocal form factors ffnl at all (k+G):
1252      ider=0; idir=0; dimffnl=1
1253      ABI_ALLOCATE(ffnl, (npw_k,dimffnl,psps%lmnmax,ntypat))
1254 
1255      call mkffnl(psps%dimekb,dimffnl,psps%ekb,ffnl,psps%ffspl,&
1256        gmet,gprimd,ider,idir,psps%indlmn,kg_k,kpg_k,kpoint,psps%lmnmax,&
1257        psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg,&
1258        npw_k,ntypat,psps%pspso,psps%qgrid_ff,rmet,&
1259        psps%usepaw,psps%useylm,ylm_k,ylmgr_dum)
1260 
1261      ABI_FREE(kpg_k)
1262      ABI_FREE(ylm_k)
1263 
1264 !    Load k-dependent part in the Hamiltonian datastructure
1265      ABI_ALLOCATE(ph3d,(2,npw_k,gs_hamk%matblk))
1266      call load_k_hamiltonian(gs_hamk,kpt_k=dtset%kptns(:,ikpt),npw_k=npw_k,istwf_k=istwf_k,&
1267 &                            kinpw_k=kinpw,kg_k=kg_k,ffnl_k=ffnl,ph3d_k=ph3d,&
1268 &                            compute_ph3d=.true.,compute_gbound=.true.)
1269 
1270      ! Prepare the call to getghc.
1271      ndat=1; lambda=zero; type_calc=0         ! For applying the whole Hamiltonian
1272      sij_opt=0; if (psps%usepaw==1) sij_opt=1 ! For PAW, <k+G|1+S|k+G"> is also needed.
1273      cpopt=-1                                 ! <p_lmn|in> (and derivatives) are computed here (and not saved)
1274 
1275      cplex_ghg=2; size_mat = cplex_ghg*(npw_k*nspinor)**2*dp*b2Mb
1276      ABI_STAT_MALLOC(ghg_mat,(cplex_ghg,npw_k*nspinor,npw_k*nspinor), ierr)
1277      write(msg,'(a,f0.3,a)')" Out-of-memory in ghg_mat. Memory required by the Hamiltonian matrix: ",size_mat," [Mb]."
1278      ABI_CHECK(ierr==0, msg)
1279      ghg_mat = huge(one)
1280 
1281      ABI_STAT_MALLOC(gtg_mat,(cplex_ghg,npw_k*nspinor,npw_k*nspinor*psps%usepaw), ierr)
1282      write(msg,'(a,f0.3,a)')" Out-of-memory in gtg_mat. Memory required by the PAW overlap operator: ",size_mat," [Mb]."
1283      ABI_CHECK(ierr==0, msg)
1284 
1285      ABI_DT_MALLOC(cwaveprj,(natom,nspinor*(1+cpopt)*gs_hamk%usepaw))
1286      if (cpopt==0) then
1287        call pawcprj_alloc(cwaveprj,0,gs_hamk%dimcprj)
1288      end if
1289 
1290      ABI_MALLOC(pwave,(2,npw_k*nspinor))
1291      pwave=zero ! Initialize plane-wave array:
1292 
1293      ABI_MALLOC(ghc  ,(2,npw_k*nspinor*ndat))
1294      ABI_MALLOC(gvnlc,(2,npw_k*nspinor*ndat))
1295      ABI_MALLOC(gsc  ,(2,npw_k*nspinor*ndat*(sij_opt+1)/2))
1296 
1297      if (dtset%prtvol > 0) call wrtout(std_out,' Calculating <G|H|G''> elements','PERS')
1298 
1299      ! Loop over the |beta,G''> component.
1300      ! Get <:|H|beta,G''> and <:|T_{PAW}|beta,G''>
1301      do igsp2=1,npw_k*nspinor
1302        ghc = zero
1303        pwave = zero
1304        pwave(1,igsp2)=one
1305 
1306        call getghc(cpopt,pwave,cwaveprj,ghc,gsc,gs_hamk,gvnlc,lambda,mpi_enreg,ndat,&
1307 &                  dtset%prtvol,sij_opt,tim_getghc,type_calc)
1308 
1309        ! Fill the upper triangle.
1310        ghg_mat(:,1:igsp2,igsp2) = ghc(:,1:igsp2)
1311        if (psps%usepaw==1) gtg_mat(:,1:igsp2,igsp2) = gsc(:,1:igsp2)
1312 
1313        pwave(1,igsp2)=zero ! Reset the |G,beta> component that has been treated.
1314      end do
1315 
1316      ABI_FREE(kinpw)
1317      ABI_FREE(ffnl)
1318      ABI_FREE(ph3d)
1319      ABI_FREE(pwave)
1320      ABI_FREE(ghc)
1321      ABI_FREE(gvnlc)
1322      ABI_FREE(gsc)
1323 
1324      if (psps%usepaw==1.and.cpopt==0) call pawcprj_free(cwaveprj)
1325      ABI_DT_FREE(cwaveprj)
1326 
1327      if (dtset%prtvol > 0) then
1328        !===========================================
1329        !=== Diagonalization of <G|H|G''> matrix ===
1330        !===========================================
1331        ABI_MALLOC(eig_ene,(nband_k))
1332        ABI_MALLOC(eig_vec,(cplex_ghg,npw_k*nspinor,nband_k))
1333 
1334        ! * Partial diagonalization
1335        write(msg,'(2a,3es16.8,3a,i5,a,i0)')ch10,&
1336        ' Begin partial diagonalization for kpt= ',kpoint,stag(isppol),ch10,&
1337        ' - Size of mat.=',npw_k*nspinor,' - # nband_k=',nband_k
1338        if (dtset%prtvol>0) call wrtout(std_out,msg,'PERS')
1339 
1340        jobz = "V"  ! Compute eigenvalues and eigenvectors.
1341        range = "I" ! the IL-th through IU-th eigenvalues will be found.
1342 
1343        if (psps%usepaw==0) then
1344          call xheevx(jobz,range,"Upper",cplex_ghg,npw_k*nspinor,ghg_mat,zero,zero,&
1345          1,nband_k,-tol8,negv,eig_ene,eig_vec,npw_k*nspinor)
1346        else
1347          call xhegvx(1,jobz,range,"Upper",cplex_ghg,npw_k*nspinor,ghg_mat,gtg_mat,zero,zero,&
1348          1,nband_k,-tol8,negv,eig_ene,eig_vec,npw_k*nspinor)
1349        end if
1350 
1351        ! Write eigenvalues.
1352        cfact=one !cfact=Ha_eV
1353        frmt1='(2x,8(1x,es9.2))' ; frmt2='(2x,8(1x,es9.2))'
1354        write(msg,'(a,3es16.8,3x,a)')' Eigenvalues in Ha for kpt= ',kpoint,stag(isppol)
1355        call wrtout(std_out,msg,'COLL')
1356 
1357        write(msg,frmt1)(eig_ene(ib)*cfact,ib=1,MIN(8,nband_k))
1358        call wrtout(std_out,msg,'COLL')
1359        if (nband_k>8) then
1360          do jj=9,nband_k,8
1361            write(msg,frmt2) (eig_ene(ib)*cfact,ib=jj,MIN(jj+7,nband_k))
1362            call wrtout(std_out,msg,'COLL')
1363          end do
1364        end if
1365 
1366 #ifdef HAVE_NETCDF
1367        ncerr = nf90_put_var(ncid, eigdiago_varid, eig_ene, start=[1,ikpt,isppol,istep], count=[nband_k,1,1,1])
1368        NCF_CHECK(ncerr)
1369 #endif
1370 
1371        ABI_FREE(eig_ene)
1372        ABI_FREE(eig_vec)
1373      end if
1374 
1375 #ifdef HAVE_NETCDF
1376      !write(std_out,*)"Writing H_GG slice for ikpt",ikpt
1377      if (isppol == 1) then
1378        ncerr = nf90_put_var(ncid, kg_varid, kg_k, start=[1,1,ikpt], count=[3,npw_k,1])
1379        NCF_CHECK(ncerr)
1380      end if
1381 
1382      ncerr = nf90_put_var(ncid, eighist_varid, eigen(:,ikpt,isppol), &
1383        start=[1,ikpt,isppol,istep], count=[nband_k,1,1,1])
1384      NCF_CHECK(ncerr)
1385 
1386      npws = nspinor * npw_k
1387      ncerr = nf90_put_var(ncid, ghg_varid, ghg_mat, start=[1,1,1,ikpt,isppol,istep], &
1388        count=[2,npws,npws,1,1,1])
1389      NCF_CHECK(ncerr)
1390 
1391      if (psps%usepaw == 1) then
1392        ncerr = nf90_put_var(ncid, pawgtg_varid, gtg_mat, start=[1,1,1,ikpt,isppol,istep], &
1393          count=[2,npws,npws,1,1,1])
1394        NCF_CHECK(ncerr)
1395      end if
1396 
1397      !write(666,*)"{istep: ",istep,", ikpt: ", ikpt, ", isppol: ", isppol, ", npws: ",npws,"}"
1398      !do jj=1,9
1399      !  write(666, "(18(es11.3))")ghg_mat(:,jj,:9)
1400      !  write(666, *)
1401      !end do
1402      !write(666, *) ghg_mat
1403 #endif
1404 
1405      ABI_FREE(ghg_mat)
1406      ABI_FREE(gtg_mat)
1407 
1408      ikg = ikg + npw_k
1409    end do ! ikpt
1410  end do ! isppol
1411 
1412  ABI_FREE(kg_k)
1413  ABI_FREE(vlocal)
1414 
1415  call destroy_hamiltonian(gs_hamk)
1416 
1417 #ifdef HAVE_NETCDF
1418  NCF_CHECK(nf90_close(ncid))
1419 #endif
1420 
1421  DBG_EXIT("COLL")
1422 
1423 end subroutine gshgg_mkncwrite

m_io_kss/k2gamma_centered [ Functions ]

[ Top ] [ m_io_kss ] [ Functions ]

NAME

  k2gamma_centered

FUNCTION

  Helper function to translate a set of wavefunctions from the k-centered G-sphere
  to the Gamma-centered G-sphere used for GW calculations.

INPUTS

  npw_k=Number of planewaves in the k-centered basis set.
  kss_npw=Number of planewaves in the Gamma-centered G-sphere.
  nspinor=Number of spinorial component.
  nbandksseff=Number of bands in input-output arrays.
  [icg]=Shift to be used when accessing the cg array. 0 if not specified (usually k_index).
  [eig_vec(2,npw_k*nspinor,nbandksseff)]=wavefunctions defined on the k-centered G-sphere.
  [cg(2,ikg+1:ikg+npw_k*nspinor*nbandksseff)]=wavefunctions defined on the k-centered G-sphere.
  ngfft(18)=Info on the FFT.
  MPI_enreg<MPI_type>=Structure gathering info on the parallelization.
  istwf_k
  ecut
  gbig(3,kss_npw)
  kg_k(3,npw_k)
  gmet(3,3)
  kpoint(3)

OUTPUT

  wfg(2,kss_npw*nspinor,nbandksseff)=Wavefunctions in the Gamma-centered representation.

NOTES

  1) icg is used only if cg is present.
  2) cg and eig_vec are mutually exclusive. One and only one can be passed to the routine.

PARENTS

      outkss

CHILDREN

      metric,mkffnl,mkkin

SOURCE

654 subroutine k2gamma_centered(kpoint,npw_k,istwf_k,ecut,kg_k,kss_npw,nspinor,nbandksseff,ngfft,gmet,&
655 &  MPI_enreg,gbig,ug,icg,cg,eig_vec)
656 
657 
658 !This section has been created automatically by the script Abilint (TD).
659 !Do not modify the following lines by hand.
660 #undef ABI_FUNC
661 #define ABI_FUNC 'k2gamma_centered'
662 !End of the abilint section
663 
664  implicit none
665 
666 !Arguments ------------------------------------
667 !scalars
668  integer,intent(in) :: nbandksseff,nspinor,kss_npw,npw_k,istwf_k
669  integer,optional,intent(in) :: icg
670  real(dp),intent(in) :: ecut
671  type(MPI_type),intent(inout) :: MPI_enreg
672 !arrays
673  integer,intent(in) :: gbig(3,kss_npw)
674  integer,intent(in) :: kg_k(3,npw_k)
675  integer,intent(in) :: ngfft(18)
676  real(dp),intent(in) :: gmet(3,3),kpoint(3)
677  real(dp),intent(out) :: ug(2,kss_npw*nspinor,nbandksseff)
678  real(dp),optional,intent(in) :: eig_vec(2,npw_k*nspinor,nbandksseff)
679  real(dp),optional,intent(in) :: cg(:,:)
680 
681 !Local variables-------------------------------
682 !scalars
683  integer,parameter :: tobox=1,tosph=-1
684  integer :: band,ispinor,spinor_shift2,spinor_shift1,ig,my_icg,ierr
685  integer :: n1,n2,n3,n4,n5,n6,ndat,full_npw_k,ii
686  character(len=500) :: msg
687 !arrays
688  integer :: identity(3,3)=RESHAPE((/1,0,0,0,1,0,0,0,1/),(/3,3/))
689  integer :: no_shift(3)=(/0,0,0/)
690  integer,allocatable :: trsl(:),full_kg_k(:,:)
691  real(dp),allocatable :: cfft(:,:,:,:)
692  real(dp),allocatable :: full_cg(:,:),tmp_cg(:,:)
693 
694 ! *********************************************************************
695 
696  if (PRESENT(cg).and.PRESENT(eig_vec)) then
697    MSG_ERROR("Both cg and eig_vec are present!")
698  end if
699 
700 ! Mapping between the gamma-centered basis set and the k-centered one.
701 ! trsl(ig)=npw_k+1 if vector ig is not inside the k-centered G-sphere.
702  ABI_ALLOCATE(trsl,(kss_npw))
703 
704  n1=ngfft(1); n2=ngfft(2); n3=ngfft(3)
705  n4=ngfft(4); n5=ngfft(5); n6=ngfft(6)
706 
707  if (istwf_k==1) then ! Full k-centered G-sphere.
708    call table_gbig2kg(npw_k,kg_k,kss_npw,gbig,trsl,ierr)
709    if (ierr/=0.and.(kss_npw>=npw_k)) then
710      MSG_ERROR(' The set of G vectors is inconsistent')
711    end if
712 
713  else  ! Calculate full kg with istwf_k=1 then do the mapping.
714    call get_kg(kpoint,1,ecut,gmet,full_npw_k,full_kg_k)
715 
716    call table_gbig2kg(full_npw_k,full_kg_k,kss_npw,gbig,trsl,ierr)
717    if (ierr/=0.and.(kss_npw>=npw_k)) then
718      MSG_ERROR(' The set of G vectors is inconsistent')
719    end if
720  end if
721  !
722  ! Branching, depending on optional arguments.
723  if (PRESENT(cg)) then
724    my_icg=0; if (PRESENT(icg)) my_icg=icg
725 
726    SELECT CASE (istwf_k)
727 
728    CASE (1)
729      do band=1,nbandksseff
730        do ispinor=1,nspinor
731          spinor_shift1=(ispinor-1)*kss_npw
732          spinor_shift2=(ispinor-1)*npw_k
733          do ig=1,kss_npw ! Retrieve the correct components
734            if (trsl(ig)<=npw_k) then
735              ug(:,ig+spinor_shift1,band)=cg(:,trsl(ig)+spinor_shift2+(band-1)*npw_k*nspinor+my_icg)
736            else
737              ug(:,ig+spinor_shift1,band)=zero
738            end if
739          end do
740        end do
741      end do
742 
743    CASE (2:9)
744 
745      ABI_CHECK(nspinor==1,"nspinor/=1!")
746      !
747      ! Convert input wfs from reduced to full G-sphere.
748      ndat=1
749      ABI_ALLOCATE(cfft,(2,n4,n5,n6*ndat))
750      ABI_ALLOCATE(full_cg,(2,full_npw_k*ndat))
751      ABI_ALLOCATE(tmp_cg,(2,npw_k*ndat))
752 
753      !write(std_out,*)"npw_k, full_kg_k",npw_k,full_npw_k
754 
755      do band=1,nbandksseff
756        ii = (band-1)*npw_k
757        tmp_cg = cg(:,my_icg+ii+1:my_icg+ii+npw_k)
758        !write(776,*)"band= ",band,tmp_cg !cg(1:,my_icg+1+ii:my_icg+ii+npw_k)
759 
760        call sphere(tmp_cg,ndat,npw_k,cfft,n1,n2,n3,n4,n5,n6,kg_k,istwf_k,tobox,MPI_enreg%me_g0,no_shift,identity,one)
761 
762        call sphere(full_cg,ndat,full_npw_k,cfft,n1,n2,n3,n4,n5,n6,full_kg_k,1,tosph,MPI_enreg%me_g0,no_shift,identity,one)
763        !write(777,*)"band= ",band,full_cg(:,:)
764 
765        do ig=1,kss_npw ! Retrieve the correct components
766          if (trsl(ig)<=full_npw_k) then
767            ug(:,ig,band)=full_cg(:,trsl(ig))
768          else
769            ug(:,ig,band)=zero
770          end if
771        end do
772      end do !band
773 
774      ABI_DEALLOCATE(cfft)
775      ABI_DEALLOCATE(tmp_cg)
776      ABI_DEALLOCATE(full_cg)
777 
778    CASE DEFAULT
779      MSG_BUG("Wrong istwf_k")
780    END SELECT
781 
782  else if (PRESENT(eig_vec)) then
783 
784    SELECT CASE (istwf_k)
785 
786    CASE (1)
787      do band=1,nbandksseff
788        do ispinor=1,nspinor
789          spinor_shift1=(ispinor-1)*kss_npw
790          spinor_shift2=(ispinor-1)*npw_k
791          do ig=1,kss_npw ! Retrieve the correct components
792            if (trsl(ig)<=npw_k) then
793              ug(:,ig+spinor_shift1,band)=eig_vec(:,trsl(ig)+spinor_shift2,band)
794            else
795              ug(:,ig+spinor_shift1,band)=zero
796            end if
797          end do
798        end do
799      end do
800 
801    CASE DEFAULT
802      write(msg,'(a,i0)')" Unsupported value for istwf_k: ",istwf_k
803      MSG_ERROR(msg)
804    END SELECT
805 
806  else
807    MSG_ERROR("neither cg not eig_vec are in input")
808  end if
809 
810  ABI_DEALLOCATE(trsl)
811  if (allocated(full_kg_k))  then
812    ABI_DEALLOCATE(full_kg_k)
813  end if
814 
815 end subroutine k2gamma_centered

m_io_kss/kss_calc_vkb [ Functions ]

[ Top ] [ m_io_kss ] [ Functions ]

NAME

  kss_calc_vkb

FUNCTION

  This routine calculates the Kleynman-Bylander form factors and its derivatives
  needed for the evaluation of the matrix elements of the dipole operator <phi1|r|phi2>.

INPUTS

  npw_k=Number of plane waves for this k-point.
  Psps<pseudopotential_type>=Structured datatype gathering information on the pseudopotentials.
  kg_k(3,npw_k)=Reduced coordinates of the G-vectors.
  kpoint(3)=The k-point in reduced coordinates.
  rprimd(3,3)=dimensional primitive translations for real space (bohr)

OUTPUT

  vkb (npw_k,Psps%ntypat,Psps%mpsang)=KB form factors.
  vkbd(npw_k,Psps%ntypat,Psps%mpsang)=KB form factor derivatives.
  vkbsign(Psps%mpsang,Psps%ntypat)   =KS dyadic sign.

NOTES

  This piece of code has been extracted from outkss.F90. The implementation is consistent
  with the KSS file formata (Fortran version) but it presents two design flaws.

   1) Pseudo with more that one projector per l-channel are not supported.
   2) Ordering of dimensions in vkb and vkbd is not optimal. We are not programming C!!!

TODO

  *) Spinorial case is not implemented.

PARENTS

      m_io_kss

CHILDREN

      metric,mkffnl,mkkin

SOURCE

1466 subroutine kss_calc_vkb(Psps,kpoint,npw_k,kg_k,rprimd,vkbsign,vkb,vkbd)
1467 
1468 
1469 !This section has been created automatically by the script Abilint (TD).
1470 !Do not modify the following lines by hand.
1471 #undef ABI_FUNC
1472 #define ABI_FUNC 'kss_calc_vkb'
1473 !End of the abilint section
1474 
1475  implicit none
1476 
1477 !Arguments ------------------------------------
1478 !scalars
1479  integer,intent(in) :: npw_k
1480  type(Pseudopotential_type),intent(in) :: Psps
1481 !arrays
1482  integer,intent(in) :: kg_k(3,npw_k)
1483  real(dp),intent(in) :: kpoint(3),rprimd(3,3)
1484  real(dp),intent(out) :: vkb (npw_k,Psps%ntypat,Psps%mpsang)
1485  real(dp),intent(out) :: vkbd(npw_k,Psps%ntypat,Psps%mpsang)
1486  real(dp),intent(out) :: vkbsign(Psps%mpsang,Psps%ntypat)
1487 
1488 !Local variables ------------------------------
1489 !scalars
1490  integer :: dimffnl,ider,idir,itypat,nkpg,il0,in
1491  integer :: il,ilmn,ig,is
1492  real(dp) :: ucvol,effmass_free,ecutsm,ecut
1493 !arrays
1494  real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3)
1495  real(dp),allocatable :: ffnl(:,:,:,:),kpg_dum(:,:),modkplusg(:)
1496  real(dp),allocatable :: ylm(:,:),ylm_gr(:,:,:),ylm_k(:,:)
1497 
1498 ! *************************************************************************
1499 
1500  DBG_ENTER("COLL")
1501 
1502  ABI_CHECK(Psps%usepaw==0,"You should not be here!")
1503  ABI_CHECK(Psps%useylm==0,"useylm/=0 not considered!")
1504 
1505  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
1506  !
1507  ! === Save KB dyadic sign (integer-valued) ===
1508  vkbsign=zero
1509  do itypat=1,Psps%ntypat
1510    il0=0
1511    do ilmn=1,Psps%lmnmax
1512      il=1+Psps%indlmn(1,ilmn,itypat)
1513      in=Psps%indlmn(3,ilmn,itypat)
1514      if (il/=il0 .and. in==1) then
1515        il0=il
1516        vkbsign(il,itypat)=DSIGN(one,Psps%ekb(ilmn,itypat))
1517      end if
1518    end do
1519  end do
1520 
1521  ! === Allocate KB form factor and derivative wrt k+G ===
1522  ! * Here we do not use correct ordering for dimensions
1523 
1524  ider=1; dimffnl=2 ! To retrieve the first derivative.
1525  idir=0; nkpg=0
1526  !
1527  ! Quantities used only if useylm==1
1528  ABI_MALLOC(ylm,(npw_k,Psps%mpsang**2*Psps%useylm))
1529  ABI_MALLOC(ylm_gr,(npw_k,3+6*(ider/2),Psps%mpsang**2*Psps%useylm))
1530  ABI_MALLOC(ylm_k,(npw_k,Psps%mpsang**2*Psps%useylm))
1531  ABI_MALLOC(kpg_dum,(npw_k,nkpg))
1532 
1533  ABI_MALLOC(ffnl,(npw_k,dimffnl,Psps%lmnmax,Psps%ntypat))
1534 
1535  call mkffnl(Psps%dimekb,dimffnl,Psps%ekb,ffnl,Psps%ffspl,gmet,gprimd,ider,idir,Psps%indlmn,&
1536    kg_k,kpg_dum,kpoint,Psps%lmnmax,Psps%lnmax,Psps%mpsang,Psps%mqgrid_ff,nkpg,npw_k,&
1537    Psps%ntypat,Psps%pspso,Psps%qgrid_ff,rmet,Psps%usepaw,Psps%useylm,ylm_k,ylm_gr)
1538 
1539  ABI_FREE(kpg_dum)
1540  ABI_FREE(ylm)
1541  ABI_FREE(ylm_gr)
1542  ABI_FREE(ylm_k)
1543 
1544  ABI_MALLOC(modkplusg,(npw_k))
1545 
1546  effmass_free=one; ecutsm=zero; ecut=HUGE(one)
1547 ! call mkkin(ecut,ecutsm,effmass_free,gmet,kg_k,modkplusg,kpoint,npw_k)
1548  call mkkin(ecut,ecutsm,effmass_free,gmet,kg_k,modkplusg,kpoint,npw_k,0,0)
1549  modkplusg(:)=SQRT(half/pi**2*modkplusg(:))
1550  modkplusg(:)=MAX(modkplusg(:),tol10)
1551 
1552  !do ig=1,npw_k
1553  ! kpg(:)= kpoint(:)+kg_k(:,ig)
1554  ! modkplusg(ig) = normv(kpg,gmet,"G")
1555  !end do
1556 
1557  ! Calculate matrix elements.
1558  vkb=zero; vkbd=zero
1559 
1560  do is=1,Psps%ntypat
1561    il0=0
1562    do ilmn=1,Psps%lmnmax
1563      il=1+Psps%indlmn(1,ilmn,is)
1564      in=Psps%indlmn(3,ilmn,is)
1565      if ((il/=il0).and.(in==1)) then
1566        il0=il
1567        if (ABS(Psps%ekb(ilmn,is))>1.0d-10) then
1568          if (il==1) then
1569            vkb (1:npw_k,is,il) = ffnl(:,1,ilmn,is)
1570            vkbd(1:npw_k,is,il) = ffnl(:,2,ilmn,is)*modkplusg(:)/two_pi
1571          else if (il==2) then
1572            vkb(1:npw_k,is,il)  = ffnl(:,1,ilmn,is)*modkplusg(:)
1573            do ig=1,npw_k
1574              vkbd(ig,is,il) = ((ffnl(ig,2,ilmn,is)*modkplusg(ig)*modkplusg(ig))+&
1575               ffnl(ig,1,ilmn,is) )/two_pi
1576            end do
1577          else if (il==3) then
1578            vkb (1:npw_k,is,il) =  ffnl(:,1,ilmn,is)*modkplusg(:)**2
1579            vkbd(1:npw_k,is,il) = (ffnl(:,2,ilmn,is)*modkplusg(:)**3+&
1580             2*ffnl(:,1,ilmn,is)*modkplusg(:) )/two_pi
1581          else if (il==4) then
1582            vkb (1:npw_k,is,il) =  ffnl(:,1,ilmn,is)*modkplusg(:)**3
1583            vkbd(1:npw_k,is,il) = (ffnl(:,2,ilmn,is)*modkplusg(:)**4+&
1584             3*ffnl(:,1,ilmn,is)*modkplusg(:)**2 )/two_pi
1585          end if
1586          vkb (:,is,il) = SQRT(4*pi/ucvol*(2*il-1)*ABS(Psps%ekb(ilmn,is)))*vkb (:,is,il)
1587          vkbd(:,is,il) = SQRT(4*pi/ucvol*(2*il-1)*ABS(Psps%ekb(ilmn,is)))*vkbd(:,is,il)
1588        else
1589          vkb (:,is,il)=zero
1590          vkbd(:,is,il)=zero
1591        end if
1592      end if
1593    end do
1594  end do
1595 
1596  ABI_FREE(ffnl)
1597  ABI_FREE(modkplusg)
1598 
1599  DBG_EXIT("COLL")
1600 
1601 end subroutine kss_calc_vkb

m_io_kss/make_gvec_kss [ Functions ]

[ Top ] [ m_io_kss ] [ Functions ]

NAME

 make_gvec_kss

FUNCTION

   Build the list of G-vectors using the KSS convention.

INPUTS

  nkpt=Number of k-points.
  nsym=Number of symmetries.
  prtvol=Verbosity option.
  symmorphi=
    0 : Old (Obsolete) implementation => Suppress inversion from symmetries list
    1 : Use input symrel, tnons.
  ecut_eff=Effective cutoff
  symrel(3,3,nsym)= Symmetry operation in real space.
  tnons(3,nsym)=Fractional translations
  kptns(3,nkpt)=K-points in reduced coordinates.

OUTPUT

  npwkss = Input: Initial guess for the number of G-vectors required. Use 0 to have the
           full list of G-vectors that form a closed shell.
           Output: Actual number of G-vectors that form a set of closed shells
  gvec_kss(:,:) = Input: null pointer. Output: gvec_kss(3,npwkss), list of G-vectors (closed shells)
  ierr=Status error

PARENTS

      gwls_hamiltonian,setup_screening,setup_sigma

CHILDREN

      metric,mkffnl,mkkin

SOURCE

854 subroutine make_gvec_kss(nkpt,kptns,ecut_eff,symmorphi,nsym,symrel,tnons,gprimd,prtvol,npwkss,gvec_kss,ierr)
855 
856 
857 !This section has been created automatically by the script Abilint (TD).
858 !Do not modify the following lines by hand.
859 #undef ABI_FUNC
860 #define ABI_FUNC 'make_gvec_kss'
861 !End of the abilint section
862 
863  implicit none
864 
865 !Arguments ------------------------------------
866 !scalars
867  integer,intent(in) :: nkpt,nsym,prtvol,symmorphi
868  integer,intent(out) :: ierr
869  integer,intent(inout) :: npwkss
870  real(dp),intent(in) :: ecut_eff
871 !arrays
872  integer,intent(in) :: symrel(3,3,nsym)
873  integer,pointer :: gvec_kss(:,:)
874  real(dp),intent(in) :: tnons(3,nsym),kptns(3,nkpt)
875  real(dp),intent(in) :: gprimd(3,3)
876 
877 !Local variables-------------------------------
878 !scalars
879  integer :: ii,ishm,maxpw,nbase
880  integer :: nrst1,nrst2,nsym2,pinv
881  integer,pointer :: gbig(:,:)
882  character(len=500) :: msg
883 !arrays
884  integer,pointer :: symrel2(:,:,:),shlim(:)
885  real(dp),pointer :: tnons2(:,:)
886 ! *********************************************************************
887 
888  ierr = 0
889  write(msg,'(2a)')ch10,' Sorting g-vecs for an output of states on an unique "big" PW basis.'
890  call wrtout(std_out,msg,'COLL')
891 
892  !ecut_eff = ecut * Dtset%dilatmx**2  ! Use ecut_eff instead of ecut_eff since otherwise
893  !
894  !============================================================
895  !=== Prepare set containing all G-vectors sorted by stars ===
896  !============================================================
897  !
898  !=== Analyze symmetry operations ===
899  if (symmorphi==0) then  ! Old (Obsolete) implementation: Suppress inversion from symmetries list:
900    nullify(symrel2,tnons2)
901    call remove_inversion(nsym,symrel,tnons,nsym2,symrel2,tnons2,pinv)
902    if (ANY(ABS(tnons2(:,1:nsym2))>tol8)) then
903      write(msg,'(3a)')&
904 &     ' Non-symmorphic operations still remain in the symmetries list ',ch10,&
905 &     ' Program does not stop but _KSS file will not be created...'
906      MSG_WARNING(msg)
907      ierr=ierr+1 ; RETURN
908    end if
909  else if (symmorphi==1) then
910 !  If in the input file symmorphi==1 all the symmetry operations are retained:
911 !  both identity and inversion (if any) as well as non-symmorphic operations.
912    nsym2=nsym ; pinv=1
913    ABI_MALLOC(symrel2,(3,3,nsym))
914    ABI_MALLOC(tnons2,(3,nsym))
915    symrel2(:,:,:)=symrel(:,:,1:nsym)
916    tnons2(:,:)   =tnons(:,1:nsym)
917  else
918    write(msg,'(a,i4,3a)')&
919 &   ' symmorphi = ',symmorphi,' while it must be 0 or 1',ch10,&
920 &   ' Program does not stop but KSS file will not be created...'
921    MSG_WARNING(msg)
922    ierr=ierr+1 ; RETURN
923  end if
924  !
925  !===================================================================
926  !==== Merge the set of k-centered G-spheres into a big set gbig ====
927  !===================================================================
928  !* Vectors in gbig are ordered by shells
929  !
930  nullify(gbig,shlim)
931  call merge_and_sort_kg(nkpt,kptns,ecut_eff,nsym2,pinv,symrel2,gprimd,gbig,prtvol,shlim_p=shlim)
932 
933  nbase = SIZE(shlim)   ! Number of independent G in the big sphere.
934  maxpw = shlim(nbase)  ! Total number of G"s in the big sphere.
935  !
936  ! * Determine optimal number of bands and G"s to be written.
937  !npwkss=Dtset%npwkss
938  if ((npwkss==0).or.(npwkss>=maxpw)) then
939    npwkss=maxpw
940    write(msg,'(5a)')&
941 &   ' Since the number of g''s to be written on file',ch10,&
942 &   ' was 0 or too large, it has been set to the max. value.,',ch10,&
943 &   ' computed from the union of the sets of G vectors for the different k-points.'
944    call wrtout(std_out,msg,'COLL')
945  end if
946 
947  ishm=0
948  do ii=1,nbase
949    if (shlim(ii)<=npwkss) then
950      ishm=ii
951    else
952      EXIT
953    end if
954  end do
955  !ishm=bisect(shlim,npwkss)
956 
957  if (shlim(ishm)/=npwkss) then
958    nrst1=shlim(ishm)
959    nrst2=MIN0(shlim(MIN0(ishm+1,nbase)),maxpw)
960    if (IABS(npwkss-nrst2)<IABS(npwkss-nrst1)) nrst1=nrst2
961    npwkss=nrst1
962    if (shlim(ishm)<npwkss) ishm=ishm+1
963    write(msg,'(3a)')&
964 &   ' The number of G''s to be written on file is not a whole number of stars ',ch10,&
965 &   ' the program set it to the nearest star limit.'
966    call wrtout(std_out,msg,'COLL')
967  end if
968 
969  write(msg,'(a,i5)')' Number of G-vectors is: ',npwkss
970  call wrtout(std_out,msg,'COLL')
971 
972  ABI_MALLOC(gvec_kss,(3,npwkss))
973  gvec_kss = gbig(:,1:npwkss)
974 
975  ABI_FREE(gbig)
976  ABI_FREE(symrel2)
977  ABI_FREE(tnons2)
978  ABI_FREE(shlim)
979 
980 end subroutine make_gvec_kss

m_io_kss/outkss [ Functions ]

[ Top ] [ m_io_kss ] [ Functions ]

NAME

 outkss

FUNCTION

  This routine creates an output file containing the Kohn-Sham electronic Structure
  for a large number of eigenstates (energies and eigen-functions).
  The resulting file (_KSS) is needed for a GW post-treatment.

 The routine drives the following operations:
  - Re-ordering G-vectors according to stars (sets of Gs related by symmetry operations).
    A set of g for all k-points is created.
  - Creating and opening the output "_KSS'" file
  - Printing out output file header information...
 ... and, for each k-point:
    According to 'kssform', either
      - Re-computing <G|H|G_prim> matrix elements for all (G, G_prim).
        Diagonalizing H in the plane-wave basis.
   or - Taking eigenvalues/vectors from congugate-gradient ones.
  - Writing out eigenvalues and eigenvectors.

INPUTS

  cg(2,mcg)=planewave coefficients of wavefunctions.
  usecprj=1 if cprj datastructure has been allocated (ONLY PAW)
  Cprj(natom,mcprj*usecprj) <type(pawcprj_type)>=
    projected input wave functions <Proj_i|Cnk> with all NL projectors (only for PAW)
    NOTE that Cprj are unsorted, see ctoprj.F90
  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)
  gmet(3,3)=reciprocal space metric tensor in bohr**-2.
  gprimd(3,3)=dimensional reciprocal space primitive translations
  Hdr <type(hdr_type)>=the header of wf, den and pot files
  kssform=govern the Kohn-Sham Structure file format
  mband=maximum number of bands
  mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
  mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol
  mgfft=maximum size of 1D FFTs
  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.
  nfft=(effective) number of FFT grid points (for this processor)
  nkpt=number of k points.
  npwarr(nkpt)=number of planewaves in basis at this k point
  nsppol=1 for unpolarized, 2 for spin-polarized
  nspden=number of density components
  nsym=number of symmetries in space group
  ntypat=number of types of atoms in unit cell.
  occ(mband*nkpt*nsppol)=occupation number for each band (usually 2) for each k.
  Pawtab(Psps%ntypat*Psps%usepaw) <type(pawtab_type)>=paw tabulated starting data
  Pawfgr<pawfgr_type>=fine grid parameters and related data
  prtvol=control print volume and debugging output
  Psps <type(pseudopotential_type)>=variables related to pseudopotentials
  rprimd(3,3)=dimensional primitive translations for real space (bohr)
  vtrial(nfft,nspden)=the trial potential
  xred(3,natom)=reduced dimensionless atomic coordinates

OUTPUT

  Output is written on file.
  ierr=Status error.

NOTES

 * This routine is maintained for legacy reasons. Abinit8 is not able to read KSS files
   anymore hence KSS files should be used only to interface Abinit with external codes
   that are still using the old KSS format.

 * The routine can be time consuming (in particular when computing
   <G|H|G_prim> elements for all (G, G_prim)) (kssform=1).
   So, it is recommended to call it once per run...

 * The IO code is not parallelized and this represents a serious bottleneck when np is large.

 * when kssform==1, the routine RE-computes all Hamiltonian terms.
   So it is equivalent to an additional electronic SC cycle.
   (This has no effect is convergence was reach...
   If not, eigenvalues/vectors may differs from the congugaste gradient ones)

 *  The KB form factors and derivatives are not calculated correctly if there are
    pseudos with more than one projector in an angular momentum channel.

 * In the ETSF output format (Dtset%iomode == 3), the complete symmetry set
   is output. So, if reading programs need only the symmorphic symmetries, they
   will need to remove themselves the non-symmorphic ones.

 * There exists two file formats:
    kssform==1 diagonalized file _KSS in real(dp) is generated.
    kssform==3 same as kssform=1 but the wavefunctions are not diagonalized
               (they are taken from conjugate-gradient ones)
    Old kssform=0 and kssform=2 are obsolete and no longer available

 TESTS
 * ETSF_IO output is tested in tests/etsf_io/t02.

PARENTS

      outscfcv

CHILDREN

      wrtout

SOURCE

1710 subroutine outkss(crystal,Dtfil,Dtset,ecut,gmet,gprimd,Hdr,&
1711 & kssform,mband,mcg,mcprj,mgfft,mkmem,MPI_enreg,mpsang,mpw,my_natom,natom,&
1712 & nfft,nkpt,npwarr,nspden,nsppol,nsym,ntypat,occ,Pawtab,Pawfgr,Paw_ij,&
1713 & prtvol,Psps,rprimd,vtrial,xred,cg,usecprj,Cprj,eigen,ierr)
1714 
1715  use m_linalg_interfaces
1716 
1717 !This section has been created automatically by the script Abilint (TD).
1718 !Do not modify the following lines by hand.
1719 #undef ABI_FUNC
1720 #define ABI_FUNC 'outkss'
1721 !End of the abilint section
1722 
1723  implicit none
1724 
1725 !Arguments ------------------------------------
1726 !scalars
1727  integer,intent(in) :: kssform,mband,mcg,mcprj,mgfft,mkmem,mpsang,mpw,my_natom,natom,usecprj
1728  integer,intent(in) :: nfft,nkpt,nsppol,nspden,nsym,ntypat,prtvol
1729  integer,intent(out) :: ierr
1730  real(dp),intent(in) :: ecut
1731  type(MPI_type),intent(inout) :: MPI_enreg
1732  type(Datafiles_type),intent(in) :: Dtfil
1733  type(Dataset_type),intent(in) :: Dtset
1734  type(Hdr_type),intent(inout) :: Hdr
1735  type(Pseudopotential_type),intent(in) :: Psps
1736  type(pawfgr_type), intent(in) :: Pawfgr
1737  type(crystal_t),intent(in) :: crystal
1738 !arrays
1739  integer,intent(in),target :: npwarr(nkpt)
1740  real(dp),intent(in) :: gmet(3,3),gprimd(3,3),occ(mband*nkpt*nsppol)
1741  real(dp),intent(in) :: rprimd(3,3)
1742  real(dp),intent(inout) :: vtrial(nfft,nspden)
1743  real(dp),intent(in) :: xred(3,natom)
1744  real(dp),intent(in) :: cg(2,mcg),eigen(mband*nkpt*nsppol)
1745  type(pawcprj_type),intent(in) :: Cprj(natom,mcprj*usecprj)
1746  type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat*Psps%usepaw)
1747  type(paw_ij_type),intent(inout),target :: Paw_ij(my_natom*Psps%usepaw)
1748 
1749 !Local variables-------------------------------
1750 !scalars
1751  integer,parameter :: tim_rwwf=0,bufnb=20
1752  integer :: untkss,onband_diago
1753  integer :: bdtot_index,i,ib,ibp,iomode
1754  integer :: ibsp,ibsp1,ibsp2,ibg,ig,ii,ikpt
1755  integer :: master,receiver,sender,spinor_shift1,shift
1756  integer :: ishm,ispinor,isppol,istwf_k,my_rank,j
1757  integer :: k_index,maxpw,mproj,n1,n2,n2dim,n3,n4,n5,n6,nband_k
1758  integer :: nbandkss_k,nbandksseff,nbase,nprocs,npw_k,onpw_k,npwkss
1759  integer :: nrst1,nrst2,nsym2,ntemp,pinv,sizepw,spaceComm,comm_self
1760  integer :: pad1,pad2
1761  integer :: bufrt,bufsz
1762  real(dp) :: cinf=1.0e24,csup=zero,einf=1.0e24,esup=zero
1763  real(dp) :: norm,cfact,ecut_eff
1764  logical :: do_diago,found,ltest,lhack
1765  logical,parameter :: skip_test_ortho=.FALSE.
1766  character(len=500) :: msg
1767  character(len=80) :: frmt1,frmt2
1768  character(len=10) :: stag(2)=(/'          ','          '/)
1769 !arrays
1770  integer :: nbandkssk(nkpt)
1771  integer,pointer :: symrel2(:,:,:)
1772  integer,pointer :: gbig(:,:)
1773  integer,pointer :: shlim(:)
1774  integer,allocatable :: kg_k(:,:)
1775  integer,allocatable :: dimlmn(:)
1776  integer :: nattyp_dum(0)
1777  real(dp) :: ovlp(2),kpoint(3),tsec(2)
1778  real(dp),pointer :: tnons2(:,:)
1779  real(dp),allocatable :: ene(:)
1780  real(dp),pointer :: eig_ene(:),eig_vec(:,:,:)
1781  real(dp),allocatable :: occ_k(:)
1782  real(dp),allocatable,target :: wfg(:,:,:)
1783  real(dp),ABI_CONTIGUOUS pointer :: ug1(:,:),ug2(:,:)
1784  type(pawcprj_type),allocatable :: Cprjnk_k(:,:)
1785  type(pawcprj_type),pointer :: Cprj_diago_k(:,:)
1786  type(ddiago_ctl_type) :: Diago_ctl
1787  type(paw_ij_type),pointer :: Paw_ij_all(:)
1788 ! *********************************************************************
1789 
1790  ABI_UNUSED(mkmem)
1791 
1792  DBG_ENTER("COLL")
1793 
1794  call timab(933,1,tsec) ! outkss
1795  call timab(934,1,tsec) ! outkss(Gsort+hd)
1796 
1797  spaceComm=MPI_enreg%comm_cell
1798  my_rank=xmpi_comm_rank(spaceComm)
1799  nprocs=xmpi_comm_size(spaceComm)
1800  master=0
1801 
1802  iomode = Dtset%iomode
1803  nullify(eig_ene)
1804  nullify(eig_vec)
1805  nullify(Cprj_diago_k)
1806 
1807  ! JB: Valgrind complains about non initialized value. Set to -1 so if an array
1808  ! should be allocated with this "unintialized value" it crashes
1809  onband_diago = -1
1810 
1811 !MG: since in seq case MPI_enreg%proc_distrb is not defined
1812 !we hack a bit the data type in order to get rid of MPI preprocessing options.
1813 !The previous status of %proc_distrb is restored before exiting.
1814 !Note that in case of seq run MPI_enreg%proc_distrb is nullified at the very beginning of abinit.F90
1815 !
1816 !FIXME this is a design flaw that should be solved: proc_distrb should always
1817 !be allocated and filled with my_rank in case of sequential run otherwise checks like
1818 !if (nprocs>1.and.MPI_enreg%proc_distrb(ii)==me) leads to SIGFAULT under gfortran.
1819 !as the second array is not allocated.
1820  lhack=.FALSE.
1821  if (nprocs==1) then
1822    ltest=allocated(MPI_enreg%proc_distrb)
1823    if (.not.ltest) then
1824      ABI_ALLOCATE(MPI_enreg%proc_distrb,(nkpt,mband,nsppol))
1825      MPI_enreg%proc_distrb=my_rank
1826      lhack=.TRUE.
1827    end if
1828    ltest=ALL(MPI_enreg%proc_distrb==my_rank)
1829    ABI_CHECK(ltest,'wrong values in %proc_distrb')
1830  end if
1831 !
1832 !============================
1833 !==== Perform some tests ====
1834 !============================
1835  ierr=0
1836 
1837  if (iomode==IO_MODE_ETSF) then
1838    write(msg,'(3a)')&
1839 &   'when iomode==3 in outkss, support for netcdf ',ch10,&
1840 &   'must be compiled. Use --enable-netcdf when configuring '
1841 #ifndef HAVE_NETCDF
1842    MSG_WARNING(msg)
1843    ierr = ierr + 1
1844 #endif
1845  end if
1846 
1847  if (kssform==3) then
1848    write(msg,'(a,70("="),4a)')ch10,ch10,&
1849 &   ' Calculating and writing out Kohn-Sham electronic Structure file',ch10, &
1850 &   ' Using conjugate gradient wavefunctions and energies (kssform=3)'
1851  else if (kssform==1) then
1852    write(msg,'(a,70("="),4a,i1,a)') ch10,ch10, &
1853 &   ' Calculating and writing out Kohn-Sham electronic Structure file',ch10, &
1854 &   ' Using diagonalized wavefunctions and energies (kssform=',kssform,')'
1855  else
1856    write(msg,'(a,i0,2a)')&
1857 &   " Unsupported value for kssform: ",kssform,ch10,&
1858 &   "  Program does not stop but _KSS file will not be created..."
1859    ierr=ierr+1
1860  end if
1861  call wrtout(std_out,msg,'COLL')
1862  call wrtout(ab_out,msg,'COLL')
1863 !
1864 !* Check whether nband is constant in metals
1865  if ( (Dtset%occopt>=2.and.Dtset%occopt<=8) .and. (ANY(Dtset%nband(1:nkpt*nsppol)/=Dtset%nband(1))) ) then
1866    write(msg,'(3a,i4,a,i3,a,i4,3a)')&
1867 &   ' The number of bands must be the same for all k-points ',ch10,&
1868 &   ' but nband(1)=',Dtset%nband(1),' is different of nband(',&
1869 &   ikpt+(isppol-1)*nkpt,')=',Dtset%nband(ikpt+(isppol-1)*nkpt),'.',ch10,&
1870 &   '  Program does not stop but _KSS file will not be created...'
1871    MSG_WARNING(msg)
1872    ierr=ierr+1
1873  end if
1874 !* istwfk must be 1 for each k-point
1875  if (ANY(Dtset%istwfk(1:nkpt)/=1).and.kssform/=3) then
1876    write(msg,'(7a)')&
1877 &   ' istwfk/=1 not allowed when kssform/=3 :',ch10,&
1878 &   ' States output not programmed for time-reversal symmetry.',ch10,&
1879 &   ' Action : change istwfk in input file (put it to 1 for all kpt).',ch10,&
1880 &   ' Program does not stop but _KSS file will not be created...'
1881    MSG_WARNING(msg)
1882    ierr=ierr+1
1883  end if
1884 !* Check spin-orbit
1885  if (Psps%mpssoang/=mpsang) then
1886    write(msg,'(3a)')&
1887 &   ' Variable mpspso should be 1 !',ch10,&
1888 &   ' Program does not stop but _KSS file will not be created...'
1889    MSG_WARNING(msg)
1890    ierr=ierr+1
1891  end if
1892 !* Check mproj
1893  mproj=MAXVAL(Psps%indlmn(3,:,:))
1894  if (mproj>1.and.Psps%usepaw==0) then ! TODO One has to derive the expression for [Vnl,r], in particular HGH and GTH psps
1895    write(msg,'(8a)')ch10,&
1896 &   ' outkss : COMMENT - ',ch10,&
1897 &   ' At least one NC pseudopotential has more that one projector per angular channel',ch10,&
1898 &   ' Note that inclvkb==0 should be used in screening, since the evaluation of the commutator',ch10,&
1899 &   ' for this particular case is not implemented yet'
1900    call wrtout(std_out,msg,'COLL')
1901    call wrtout(ab_out,msg,'COLL')
1902  end if
1903 !* Check max angular momentum
1904  if (MAXVAL(Psps%indlmn(1,:,:))+1 >= 5) then
1905    write(msg,'(3a)')&
1906 &   ' Pseudopotentials with f-projectors not implemented',ch10,&
1907 &   ' Program does not stop but _KSS file will not be created...'
1908    MSG_WARNING(msg)
1909    ierr=ierr+1
1910  end if
1911 !* Check useylm
1912  if (Psps%useylm/=0.and.Psps%usepaw==0) then
1913    write(msg,'(3a)')&
1914 &   ' The present version of outkss does not work with useylm/=0 !',ch10,&
1915 &   ' Program does not stop but _KSS file will not be created...'
1916    MSG_WARNING(msg)
1917    ierr=ierr+1
1918  end if
1919 !* Check PAW and kssform value
1920  if (Psps%usepaw/=0) then
1921    if (nprocs>1.and.kssform==1) then
1922      write(msg,'(3a)')&
1923 &     ' Parallel PAW with kssform=1, not yet allowed',ch10,&
1924 &     ' Program does not stop but _KSS file will not be created...'
1925      MSG_WARNING(msg)
1926      ierr=ierr+1
1927    end if
1928    if (kssform==3.and.usecprj/=1) then
1929      write(msg,'(3a)')&
1930 &     ' If PAW and kssform=3, usecprj must be 1',ch10,&
1931 &     ' Program does not stop but _KSS file will not be created...'
1932      MSG_WARNING(msg)
1933      ierr=ierr+1
1934    end if
1935  end if
1936 !* Check parallelization
1937  if (MPI_enreg%paralbd/=0) then
1938    write(msg,'(3a)')&
1939 &   ' outkss cannot be used with parallelization on bands (paralbd/=0) !',ch10,&
1940 &   ' Program does not stop but _KSS file will not be created...'
1941    MSG_WARNING(msg)
1942    ierr=ierr+1
1943  end if
1944  if (MPI_enreg%paral_spinor/=0) then
1945    write(msg,'(3a)')&
1946 &   ' outkss cannot be used yet with parallelization on nspinors !',ch10,&
1947 &   ' Program does not stop but _KSS file will not be created...'
1948    MSG_WARNING(msg)
1949    ierr=ierr+1
1950 
1951  endif
1952  if (ierr/=0) then
1953    write(msg,'(3a)')&
1954 &   ' outkss: Not allowed options found !',ch10,&
1955 &   ' Program does not stop but _KSS file will not be created...'
1956    call wrtout(std_out,msg,'COLL')
1957    call wrtout(ab_out,msg,'COLL')
1958    write(msg,'(a)')' outkss: see the log file for more information.'
1959    call wrtout(ab_out,msg,'COLL')
1960    RETURN ! Houston we have a problem!
1961  end if
1962 !
1963 !Estimate required memory in case of diagonalization.
1964 !TODO to be modified to take into account the case nsppol=2
1965  if (kssform/=3) then
1966    call memkss(mband,mgfft,mproj,Psps%mpssoang,mpw,natom,Dtset%ngfft,nkpt,dtset%nspinor,nsym,ntypat)
1967  end if
1968 !
1969 !=== Initialize some variables ===
1970  if (nsppol==2) stag(:)=(/'SPIN UP:  ','SPIN DOWN:'/)
1971  n1=Dtset%ngfft(1); n2=Dtset%ngfft(2); n3=Dtset%ngfft(3)
1972  n4=Dtset%ngfft(4); n5=Dtset%ngfft(5); n6=Dtset%ngfft(6)
1973  ecut_eff = ecut * Dtset%dilatmx**2  ! Use ecut_eff instead of ecut_eff since otherwise
1974 !one cannot restart from a previous density file
1975  sizepw=2*mpw ; do_diago=(kssform/=3)
1976  ABI_ALLOCATE(dimlmn,(natom*Psps%usepaw))
1977  if (Psps%usepaw==1) then
1978    call pawcprj_getdim(dimlmn,natom,nattyp_dum,ntypat,Dtset%typat,pawtab,'R')
1979  end if
1980 !
1981 !============================================================
1982 !=== Prepare set containing all G-vectors sorted by stars ===
1983 !============================================================
1984  write(msg,'(2a)')ch10,' Sorting g-vecs for an output of states on an unique "big" PW basis.'
1985  call wrtout(std_out,msg,'COLL')
1986 !
1987 !=== Analyze symmetry operations ===
1988  if (Dtset%symmorphi==0) then  ! Old (Obsolete) implementation: Suppress inversion from symmetries list:
1989    nullify(symrel2,tnons2)
1990    call remove_inversion(nsym,Dtset%symrel,Dtset%tnons,nsym2,symrel2,tnons2,pinv)
1991    if (ANY(ABS(tnons2(:,1:nsym2))>tol8)) then
1992      write(msg,'(3a)')&
1993 &     ' Non-symmorphic operations still remain in the symmetries list ',ch10,&
1994 &     ' Program does not stop but _KSS file will not be created...'
1995      MSG_WARNING(msg)
1996      ierr=ierr+1 ; RETURN
1997    end if
1998  else if (Dtset%symmorphi==1) then
1999 !  If in the input file symmorphi==1 all the symmetry operations are retained:
2000 !  both identity and inversion (if any) as well as non-symmorphic operations.
2001    nsym2=nsym ; pinv=1
2002    ABI_ALLOCATE(symrel2,(3,3,nsym))
2003    ABI_ALLOCATE(tnons2,(3,nsym))
2004    symrel2(:,:,:)=Dtset%symrel(:,:,1:nsym)
2005    tnons2(:,:)   =Dtset%tnons(:,1:nsym)
2006  else
2007    write(msg,'(a,i4,3a)')&
2008 &   ' symmorphi = ',Dtset%symmorphi,' while it must be 0 or 1',ch10,&
2009 &   ' Program does not stop but KSS file will not be created...'
2010    MSG_WARNING(msg)
2011    ierr=ierr+1 ; RETURN
2012  end if
2013 !
2014 !===================================================================
2015 !==== Merge the set of k-centered G-spheres into a big set gbig ====
2016 !===================================================================
2017 !* Vectors in gbig are ordered by shells
2018 !
2019  nullify(gbig,shlim)
2020  call merge_and_sort_kg(nkpt,Dtset%kptns,ecut_eff,nsym2,pinv,symrel2,gprimd,gbig,prtvol,shlim_p=shlim)
2021 
2022  nbase = SIZE(shlim)   ! Number of independent G in the big sphere.
2023  maxpw = shlim(nbase)  ! Total number of G"s in the big sphere.
2024 !
2025 !* Determine optimal number of bands and G"s to be written.
2026  npwkss=Dtset%npwkss
2027  if ((npwkss==0).or.(npwkss>=maxpw)) then
2028    npwkss=maxpw
2029    write(msg,'(5a)')&
2030 &   ' Since the number of g''s to be written on file',ch10,&
2031 &   ' was 0 or too large, it has been set to the max. value.,',ch10,&
2032 &   ' computed from the union of the sets of G vectors for the different k-points.'
2033    call wrtout(std_out,msg,'COLL')
2034  end if
2035 
2036  ishm=0
2037  do ii=1,nbase
2038    if (shlim(ii)<=npwkss) then
2039      ishm=ii
2040    else
2041      EXIT
2042    end if
2043  end do
2044 
2045  if (shlim(ishm)/=npwkss) then
2046    nrst1=shlim(ishm)
2047    nrst2=MIN0(shlim(MIN0(ishm+1,nbase)),maxpw)
2048    if (IABS(npwkss-nrst2)<IABS(npwkss-nrst1)) nrst1=nrst2
2049    npwkss=nrst1
2050    if (shlim(ishm)<npwkss) ishm=ishm+1
2051    write(msg,'(3a)')&
2052 &   ' The number of G''s to be written on file is not a whole number of stars ',ch10,&
2053 &   ' the program set it to the nearest star limit.'
2054    call wrtout(std_out,msg,'COLL')
2055  end if
2056 
2057  write(msg,'(a,i5)')' Number of g-vectors written on file is: ',npwkss
2058  call wrtout(std_out,msg,'COLL')
2059 !
2060 !=== Check on the number of stored bands ===
2061  if (do_diago) then
2062 
2063    if (Dtset%nbandkss==-1.or.Dtset%nbandkss>=maxpw) then
2064      nbandkssk(1:nkpt)=npwarr(1:nkpt)
2065      write(msg,'(6a)')ch10,&
2066 &     ' Since the number of bands to be computed was (-1) or',ch10,&
2067 &     ' too large, it has been set to the max. value. allowed for each k,',ch10,&
2068 &     ' thus, the minimum of the number of plane waves for each k point.'
2069      call wrtout(std_out,msg,'COLL')
2070    else
2071      nbandkssk(1:nkpt)=Dtset%nbandkss
2072      found=.FALSE.
2073      do ikpt=1,nkpt
2074        if (Dtset%nbandkss>npwarr(ikpt)) then
2075          nbandkssk(ikpt)=npwarr(ikpt)
2076          found=.TRUE.
2077        end if
2078      end do
2079      if (found) then
2080        write(msg,'(7a)')&
2081 &       ' The value choosen for the number of bands in file',ch10,&
2082 &       ' (nbandkss) was greater than at least one number of plane waves ',ch10,&
2083 &       ' for a given k-point (npw_k).',ch10,' It has been modified consequently.'
2084        MSG_WARNING(msg)
2085      end if
2086    end if
2087    found=.FALSE.
2088    do ikpt=1,nkpt
2089      if (nbandkssk(ikpt)>npwkss) then
2090        nbandkssk(ikpt)=npwkss
2091        found=.TRUE.
2092      end if
2093    end do
2094    if (found) then
2095      write(msg,'(5a)')&
2096 &     ' The number of bands to be computed (for one k) was',ch10,&
2097 &     ' greater than the number of g-vectors to be written.',ch10,&
2098 &     ' It has been modified consequently.'
2099      MSG_WARNING(msg)
2100    end if
2101    nbandksseff=MINVAL(nbandkssk)
2102 
2103  else ! .not. do_diago
2104    do ikpt=1,nkpt
2105      do isppol=1,nsppol
2106        nbandkssk(ikpt)=Dtset%nband(ikpt+(isppol-1)*nkpt)
2107      end do
2108    end do
2109    nbandksseff=MINVAL(nbandkssk)
2110    if (Dtset%nbandkss>0 .and. Dtset%nbandkss<nbandksseff) then
2111      write(msg,'(a,i5,a,i5,2a)')&
2112 &     ' Number of bands calculated=',nbandksseff,', greater than nbandkss=',Dtset%nbandkss,ch10,&
2113 &     ' will write nbandkss bands on the KSS file'
2114      MSG_COMMENT(msg)
2115      nbandksseff=Dtset%nbandkss
2116    end if
2117  end if
2118 
2119  write(msg,'(a,i5)')' Number of bands written on file is: ',nbandksseff
2120  call wrtout(std_out,msg,'COLL')
2121 
2122  found= ANY(nbandkssk(1:nkpt)<npwarr(1:nkpt))
2123 
2124  if (do_diago) then
2125    if (found) then
2126      write(msg,'(6a)')ch10,&
2127 &     ' Since the number of bands to be computed',ch10,&
2128 &     ' is less than the number of G-vectors found,',ch10,&
2129 &     ' the program will perform partial diagonalizations.'
2130    else
2131      write(msg,'(6a)')ch10,&
2132 &     ' Since the number of bands to be computed',ch10,&
2133 &     ' is equal to the nb of G-vectors found for each k-pt,',ch10,&
2134 &     ' the program will perform complete diagonalizations.'
2135    end if
2136    call wrtout(std_out,msg,'COLL')
2137  end if
2138 !
2139 !==========================================================================
2140 !=== Open KSS file for output, write header with dimensions and kb sign ===
2141 !==========================================================================
2142 !
2143 !* Output required disk space.
2144  call dsksta(ishm,Psps%usepaw,nbandksseff,mpsang,natom,ntypat,npwkss,nkpt,dtset%nspinor,nsppol,nsym2,dimlmn)
2145 
2146  if (my_rank==master) then
2147    call write_kss_header(dtfil%fnameabo_kss,npwkss,ishm,nbandksseff,mband,nsym2,symrel2,tnons2,occ,gbig,shlim,&
2148 &   crystal,Dtset,Hdr,Psps,iomode,untkss)
2149  end if
2150 
2151  ABI_DEALLOCATE(shlim)
2152 
2153  if (     do_diago) msg = ' Diagonalized eigenvalues'
2154  if (.not.do_diago) msg = ' Conjugate gradient eigenvalues'
2155  call wrtout(ab_out,msg,'COLL')
2156 
2157  if (Dtset%enunit==1) then
2158    msg='   k    eigenvalues [eV]'
2159  else
2160    msg='   k    eigenvalues [Hartree]'
2161  end if
2162  call wrtout(ab_out,msg,'COLL')
2163 !
2164 !=== In case of PAW distributed atomic sites, need to retrieve the full paw_ij%dij ===
2165  if (do_diago.and.Psps%usepaw==1.and.MPI_enreg%nproc_atom>1) then
2166    ABI_DATATYPE_ALLOCATE(Paw_ij_all,(dtset%natom))
2167    call paw_ij_gather(Paw_ij,Paw_ij_all,-1,MPI_enreg%comm_atom)
2168  else
2169    paw_ij_all => paw_ij
2170  end if
2171 
2172 
2173  call timab(934,2,tsec) ! outkss(Gsort+hd)
2174 !
2175 
2176  k_index=0; bdtot_index=0; ibg=0
2177 
2178  do isppol=1,nsppol ! Loop over spins
2179 !
2180    do ikpt=1,nkpt ! Loop over k-points.
2181      call timab(935,1,tsec) ! outkss(k-Loop)
2182 
2183      nband_k   =Dtset%nband(ikpt+(isppol-1)*nkpt)
2184      npw_k     =npwarr(ikpt)
2185      istwf_k   =Dtset%istwfk(ikpt)
2186      kpoint    =Dtset%kptns(:,ikpt)
2187      nbandkss_k=nbandkssk(ikpt)
2188 
2189      ! Get G-vectors, for this k-point.
2190      call get_kg(kpoint,istwf_k,ecut_eff,gmet,onpw_k,kg_k)
2191      ABI_CHECK(onpw_k==npw_k,"Mismatch in npw_k")
2192 !
2193 !    ============================================
2194 !    ==== Parallelism over k-points and spin ====
2195 !    ============================================
2196      if (MPI_enreg%proc_distrb(ikpt,1,isppol)==my_rank) then
2197 
2198        write(msg,'(2a,i3,3x,a)')ch10,' k-point ',ikpt,stag(isppol)
2199        call wrtout(std_out,msg,'PERS')
2200 
2201        if (do_diago) then ! Direct diagonalization of the KS Hamiltonian.
2202          if (associated(eig_ene))  then
2203            ABI_DEALLOCATE(eig_ene)
2204          end if
2205          if (associated(eig_vec))  then
2206            ABI_DEALLOCATE(eig_vec)
2207          end if
2208          comm_self = xmpi_comm_self
2209 
2210          call timab(936,1,tsec)
2211 
2212          call init_ddiago_ctl(Diago_ctl,"Vectors",isppol,dtset%nspinor,ecut_eff,Dtset%kptns(:,ikpt),Dtset%nloalg,gmet,&
2213 &         nband_k=nbandkssk(ikpt),effmass_free=Dtset%effmass_free,istwf_k=Dtset%istwfk(ikpt),prtvol=Dtset%prtvol)
2214 
2215          call ksdiago(Diago_ctl,nbandkssk(ikpt),Dtset%nfft,mgfft,Dtset%ngfft,natom,&
2216 &         Dtset%typat,nfft,dtset%nspinor,nspden,nsppol,Pawtab,Pawfgr,Paw_ij_all,&
2217 &         Psps,rprimd,vtrial,xred,onband_diago,eig_ene,eig_vec,Cprj_diago_k,comm_self,ierr)
2218 
2219          call timab(936,2,tsec)
2220        end if
2221 
2222      end if ! END of kpt+spin parallelism.
2223 !
2224 !    ===========================================================
2225 !    ==== Transfer data between master and the working proc ====
2226 !    ===========================================================
2227      call timab(937,1,tsec) !outkss(MPI_exch)
2228      ABI_DATATYPE_ALLOCATE(Cprjnk_k,(0,0))
2229      if (nprocs==1) then
2230 
2231        if (Psps%usepaw==1) then ! Copy projectors for this k-point
2232          n2dim=min(nbandksseff*dtset%nspinor,onband_diago)
2233          if (kssform==3) n2dim=nband_k*dtset%nspinor
2234          ABI_DATATYPE_DEALLOCATE(Cprjnk_k)
2235          ABI_DATATYPE_ALLOCATE(Cprjnk_k,(natom,n2dim))
2236          call pawcprj_alloc(Cprjnk_k,0,dimlmn)
2237          if (kssform==3) then
2238            call pawcprj_copy(Cprj(:,ibg+1:ibg+dtset%nspinor*nband_k),Cprjnk_k)
2239          else
2240           !MSG_WARNING("Here I have to use onband_diago") !FIXME
2241            call pawcprj_copy(Cprj_diago_k(:,1:n2dim),Cprjnk_k)
2242          end if
2243        end if
2244 
2245      else !parallel case
2246 
2247        receiver=master; sender=MPI_enreg%proc_distrb(ikpt,1,isppol)
2248 
2249        bufsz=nbandksseff/bufnb; bufrt=nbandksseff-bufnb*bufsz
2250 
2251        if (my_rank==receiver.or.my_rank==sender) then
2252 
2253          if (do_diago.and.(my_rank==receiver.and.my_rank/=sender)) then ! Alloc arrays if not done yet.
2254            ABI_ALLOCATE(eig_ene,(npw_k*dtset%nspinor))
2255            ABI_ALLOCATE(eig_vec,(2,npw_k*dtset%nspinor,nbandkssk(ikpt)))
2256          end if
2257 
2258          if (.not.do_diago) then
2259 
2260            ABI_ALLOCATE(eig_vec,(2,npw_k*dtset%nspinor,nbandkssk(ikpt)))
2261 
2262            if (my_rank==sender) then
2263              do ib=1,nbandksseff
2264                shift = k_index + (ib-1)*npw_k*dtset%nspinor
2265                do ig=1,npw_k*dtset%nspinor
2266                  eig_vec(:,ig,ib)=cg(:,ig+shift)
2267                end do
2268              end do
2269            end if
2270 !
2271 !          In case of PAW and kssform==3, retrieve matrix elements of the PAW projectors for this k-point
2272            if (Psps%usepaw==1) then
2273              n2dim=min(nbandksseff*dtset%nspinor,onband_diago)
2274              if (kssform==3) n2dim=nband_k*dtset%nspinor
2275              ABI_DATATYPE_DEALLOCATE(Cprjnk_k)
2276              ABI_DATATYPE_ALLOCATE(Cprjnk_k,(natom,n2dim))
2277              call pawcprj_alloc(Cprjnk_k,0,dimlmn)
2278              if (my_rank==sender) then
2279                if (kssform==3) then
2280                  call pawcprj_copy(Cprj(:,ibg+1:ibg+dtset%nspinor*nband_k),Cprjnk_k)
2281                else
2282                 !MSG_WARNING("Here I have to use onband_diago") !FIXME
2283                  call pawcprj_copy(Cprj_diago_k(:,1:n2dim),Cprjnk_k)
2284                end if
2285              end if
2286              if (sender/=receiver) then
2287                call pawcprj_mpi_exch(natom,n2dim,dimlmn,0,Cprjnk_k,Cprjnk_k,sender,receiver,spaceComm,ierr)
2288              end if
2289            end if ! usepaw
2290 
2291          else ! do_diago
2292            call xmpi_exch(eig_ene,nbandksseff,sender,eig_ene,receiver,spaceComm,ierr)
2293          end if
2294 
2295 !        Exchange eigenvectors.
2296          if (bufsz>0) then
2297            do i=0,bufnb-1
2298              call xmpi_exch(eig_vec(:,:,i*bufsz+1:(i+1)*bufsz),2*npw_k*dtset%nspinor*bufsz,&
2299 &             sender,eig_vec(:,:,i*bufsz+1:(i+1)*bufsz),receiver,spaceComm,ierr)
2300            end do
2301          end if
2302          if (bufrt>0) then
2303            call xmpi_exch(eig_vec(:,:,bufnb*bufsz+1:bufnb*bufsz+bufrt),2*npw_k*dtset%nspinor*bufrt,&
2304 &           sender,eig_vec(:,:,bufnb*bufsz+1:bufnb*bufsz+bufrt),receiver,spaceComm,ierr)
2305          end if
2306 
2307        end if
2308      end if !nprocs > 1
2309      call timab(937,2,tsec) !outkss(MPI_exch)
2310 
2311      call timab(938,1,tsec) !outkss(write)
2312 
2313      if (my_rank==master) then ! Prepare data for writing on disk.
2314        ABI_ALLOCATE(ene,(nbandksseff))
2315        ABI_ALLOCATE(wfg,(2,npwkss*dtset%nspinor,nbandksseff))
2316        ene=zero; wfg=zero
2317 
2318        if (.not.do_diago) then
2319          ene(1:nbandksseff)=eigen(1+bdtot_index:nbandksseff+bdtot_index)
2320 
2321          if (nprocs>1) then
2322            call k2gamma_centered(kpoint,npw_k,istwf_k,ecut_eff,kg_k,npwkss,dtset%nspinor,nbandksseff,Dtset%ngfft,gmet,&
2323 &           MPI_enreg,gbig,wfg,eig_vec=eig_vec)
2324          else
2325            call k2gamma_centered(kpoint,npw_k,istwf_k,ecut_eff,kg_k,npwkss,dtset%nspinor,nbandksseff,Dtset%ngfft,gmet,&
2326 &           MPI_enreg,gbig,wfg,icg=k_index,cg=cg)
2327          end if
2328 
2329        else ! Direct diagonalization.
2330          ene(1:nbandksseff)=eig_ene(1:nbandksseff)
2331 
2332 !        FIXME: receiver does not know Diago_ctl%npw_k
2333          call k2gamma_centered(kpoint,npw_k,istwf_k,ecut_eff,kg_k,npwkss,dtset%nspinor,nbandksseff,Dtset%ngfft,gmet,&
2334 &         MPI_enreg,gbig,wfg,eig_vec=eig_vec)
2335 
2336 !        * Check diagonalized eigenvalues with respect to conjugate gradient ones
2337          ntemp=MIN(nbandksseff,nband_k)
2338          if (ANY(ABS(ene(1:ntemp)-eigen(1+bdtot_index:ntemp+bdtot_index))>tol3)) then
2339            write(msg,'(3a)')&
2340 &           ' The diagonalized eigenvalues differ by more than 10^-3 Hartree',ch10,&
2341 &           ' with respect to the conjugated gradient values.'
2342            MSG_WARNING(msg)
2343          end if
2344        end if
2345 !
2346 !      * Write out energies
2347        if (Dtset%enunit==1) then
2348          cfact=Ha_eV ; frmt1='(i4,4x,9(1x,f7.2))' ; frmt2='(8x,9(1x,f7.2))'
2349          write(msg,'(a,i3,3x,a)')' Eigenvalues in eV for ikpt= ',ikpt,stag(isppol)
2350        else
2351          cfact=one   ; frmt1='(i4,4x,9(1x,f7.4))' ; frmt2='(8x,9(1x,f7.4))'
2352          write(msg,'(a,i3,3x,a)')' Eigenvalues in Hartree for ikpt= ',ikpt,stag(isppol)
2353        end if
2354        call wrtout(std_out,msg,'COLL')
2355 
2356        write(msg,frmt1)ikpt,(ene(ib)*cfact,ib=1,MIN(9,nbandksseff))
2357        call wrtout(std_out,msg,'COLL')
2358        call wrtout(ab_out,msg,'COLL')
2359 
2360        if (nbandksseff>9) then
2361          do j=10,nbandksseff,9
2362            write(msg,frmt2) (ene(ib)*cfact,ib=j,MIN(j+8,nbandksseff))
2363            call wrtout(std_out,msg,'COLL')
2364            call wrtout(ab_out,msg,'COLL')
2365          end do
2366        end if
2367 
2368        if (skip_test_ortho) then ! Set this if to FALSE to skip test below
2369          einf=one; esup=one; cinf=zero; csup=zero
2370        else
2371          !
2372          ! Test on the normalization of wavefunctions.
2373          ibsp=0
2374          do ib=1,nbandksseff
2375            norm=zero
2376            do ispinor=1,dtset%nspinor
2377              ibsp=ibsp+1
2378              spinor_shift1=(ispinor-1)*npwkss
2379              ug1 => wfg(:,1+spinor_shift1:npwkss+spinor_shift1,ib)
2380 
2381              !ovlp(1) =ddot(npwkss,ug1(1,:),1,ug1(1,:),1) + ddot(npwkss,ug1(2,:),1,ug1(2,:),1)
2382              ovlp(1) = cg_dznrm2(npwkss,ug1)
2383              ovlp(1) = ovlp(1)**2
2384              ovlp(2) = zero
2385              if (Psps%usepaw==1) ovlp = ovlp &
2386 &               + paw_overlap(Cprjnk_k(:,ibsp:ibsp),Cprjnk_k(:,ibsp:ibsp),Dtset%typat,Pawtab,&
2387 &                             spinor_comm=MPI_enreg%comm_spinor)
2388              norm = norm + DABS(ovlp(1))
2389            end do
2390            if (norm<einf) einf=norm
2391            if (norm>esup) esup=norm
2392          end do
2393 !
2394 !        Test on the orthogonalization of wavefunctions.
2395          do ib=1,nbandksseff
2396            pad1=(ib-1)*dtset%nspinor
2397            do ibp=ib+1,nbandksseff
2398              pad2=(ibp-1)*dtset%nspinor
2399              ovlp(:)=zero
2400              do ispinor=1,dtset%nspinor
2401                ibsp1=pad1+ispinor
2402                ibsp2=pad2+ispinor
2403                spinor_shift1=(ispinor-1)*npwkss
2404                ug1 => wfg(:,1+spinor_shift1:npwkss+spinor_shift1,ib )
2405                ug2 => wfg(:,1+spinor_shift1:npwkss+spinor_shift1,ibp)
2406 
2407                !ovlp(1)=ddot(npwkss,ug1(1,:),1,ug2(1,:),1) + ddot(npwkss,ug1(2,:),1,ug2(2,:),1)
2408                !ovlp(2)=ddot(npwkss,ug1(1,:),1,ug2(2,:),1) - ddot(npwkss,ug1(2,:),1,ug2(1,:),1)
2409                ovlp = cg_zdotc(npwkss,ug1,ug2)
2410 
2411                if (Psps%usepaw==1) ovlp= ovlp &
2412 &                 + paw_overlap(Cprjnk_k(:,ibsp1:ibsp1),Cprjnk_k(:,ibsp2:ibsp2),Dtset%typat,Pawtab,&
2413 &                               spinor_comm=MPI_enreg%comm_spinor)
2414              end do
2415              norm = DSQRT(ovlp(1)**2+ovlp(2)**2)
2416              if (norm<cinf) cinf=norm
2417              if (norm>csup) csup=norm
2418            end do
2419          end do
2420        end if
2421 
2422        write(msg,'(a,i3,3x,a)')' Writing out eigenvalues/vectors for ikpt=',ikpt,stag(isppol)
2423        call wrtout(std_out,msg,'COLL')
2424 !
2425 !      * Write occupation numbers on std_out.
2426        ABI_ALLOCATE(occ_k,(MAX(nband_k,nbandksseff)))
2427        occ_k(1:nband_k)=occ(1+bdtot_index:nband_k+bdtot_index)
2428        if (nband_k < nbandksseff) occ_k(nband_k+1:nbandksseff)=zero
2429 
2430        write(msg,'(a,i3,3x,a)')' Occupation numbers for ikpt=',ikpt,stag(isppol)
2431        call wrtout(std_out,msg,'COLL')
2432        write(msg,'(i4,4x,9(1x,f7.4))')ikpt,(occ_k(ib),ib=1,MIN(9,nbandksseff))
2433        call wrtout(std_out,msg,'COLL')
2434        if (nbandksseff>9) then
2435          do j=10,nbandksseff,9
2436            write(msg,'(8x,9(1x,f7.4))') (occ_k(ib),ib=j,MIN(j+8,nbandksseff))
2437            call wrtout(std_out,msg,'COLL')
2438          end do
2439        end if
2440 !
2441 !      =================================================================
2442 !      ==== Write wavefunctions, KB and PAW matrix elements on disk ====
2443 !      =================================================================
2444        call write_kss_wfgk(untkss,ikpt,isppol,kpoint,dtset%nspinor,npwkss,&
2445 &           nbandksseff,natom,Psps,ene,occ_k,rprimd,gbig,wfg,Cprjnk_k,iomode)
2446 
2447        ABI_DEALLOCATE(occ_k)
2448        ABI_DEALLOCATE(ene)
2449        ABI_DEALLOCATE(wfg)
2450 
2451      end if ! my_rank==master
2452      call timab(938,2,tsec) !outkss(write)
2453 
2454      if (my_rank==master.or.my_rank==MPI_enreg%proc_distrb(ikpt,1,isppol)) then
2455        if (associated(eig_ene))  then
2456          ABI_DEALLOCATE(eig_ene)
2457        end if
2458        if (associated(eig_vec))  then
2459          ABI_DEALLOCATE(eig_vec)
2460        end if
2461 
2462        if (Psps%usepaw==1) then
2463          call pawcprj_free(Cprjnk_k)
2464        end if
2465      end if
2466      ABI_DATATYPE_DEALLOCATE(Cprjnk_k)
2467 
2468      if (allocated(kg_k))  then
2469        ABI_DEALLOCATE(kg_k)
2470      end if
2471 
2472 !    if (MPI_enreg%paral_compil_kpt==1) then !cannot be used in seq run!
2473      if (.not.(proc_distrb_cycle(MPI_enreg%proc_distrb,ikpt,1,nband_k,isppol,my_rank))) then
2474        k_index=k_index+npw_k*nband_k*dtset%nspinor
2475        ibg=ibg+dtset%nspinor*nband_k
2476      end if
2477      bdtot_index=bdtot_index+nband_k
2478 
2479      call xmpi_barrier(spaceComm) ! FIXME this barrier is detrimental in the case of direct diago!
2480 
2481      call timab(935,2,tsec) !outkss(k-loop)
2482    end do ! ! End loop over k-points.
2483  end do ! spin
2484 
2485  write(msg,'(3a,f9.6,2a,f9.6,4a,f9.6,2a,f9.6,a)')&
2486 & ' Test on the normalization of the wavefunctions',ch10,&
2487 & '  min sum_G |a(n,k,G)| = ',einf,ch10,&
2488 & '  max sum_G |a(n,k,G)| = ',esup,ch10,&
2489 & ' Test on the orthogonalization of the wavefunctions',ch10,&
2490 & '  min sum_G a(n,k,G)a(n'',k,G) = ',cinf,ch10,&
2491 & '  max sum_G a(n,k,G)a(n'',k,G) = ',csup,ch10
2492  call wrtout(std_out,msg,'COLL')
2493  call wrtout(ab_out,msg,'COLL')
2494 
2495  ABI_DEALLOCATE(gbig)
2496  ABI_DEALLOCATE(symrel2)
2497  ABI_DEALLOCATE(tnons2)
2498  if (Psps%usepaw==1)  then
2499    ABI_DEALLOCATE(dimlmn)
2500    if (do_diago.and.MPI_enreg%nproc_atom>1) then
2501      ABI_DATATYPE_DEALLOCATE(Paw_ij_all)
2502    end if
2503  end if
2504 !
2505 !* Close file
2506  if (my_rank==master) then
2507    if (iomode==IO_MODE_FORTRAN) close(unit=untkss)
2508 #if defined HAVE_NETCDF
2509    if (iomode==IO_MODE_ETSF) then
2510      NCF_CHECK(nf90_close(untkss))
2511    end if
2512 #endif
2513  end if
2514 
2515  if (associated(Cprj_diago_k)) then
2516    call pawcprj_free(Cprj_diago_k)
2517    ABI_DATATYPE_DEALLOCATE(Cprj_diago_k)
2518  end if
2519 
2520  if (lhack)  then
2521    ABI_DEALLOCATE(MPI_enreg%proc_distrb)
2522  end if
2523 
2524  call wrtout(std_out, "outkss done", "COLL")
2525  call xmpi_barrier(spaceComm)
2526 
2527  DBG_EXIT("COLL")
2528  call timab(933,2,tsec) ! outkss
2529 
2530 contains

m_io_kss/write_kss_header [ Functions ]

[ Top ] [ m_io_kss ] [ Functions ]

NAME

  write_kss_header

FUNCTION

  Write the header of the KSS file either using plain Fortran-IO or netcdf with ETSF-IO format.
  Returns the unit number to be used for further writing.
  It should be executed by master node only.

INPUTS

  filekss(len=fnlen)=The name of the KSS file.
  kss_npw=Number of planewaves used for the wavefunctions in the KSS files.
  ishm=Max number of shells written on file
  shlim(ishm)=The cumulative number of G"s in each shell.
  nbandksseff=Number of bands to be written.
  mband=The maximum number of bands treated by abinit.
  nsym2=Number of symmetry operations to be written on the header.
  symrel2(3,3,nsym2)=The symmetry operations in real space to be written.
  tnons2(3,nsym2)=The fractional translations associateed to symrel2.
  gbig(3,kss_npw)=The set of G-vectors for the KSS wavefunctions (Gamma-centered)
  Hdr<hdr_type>=The abinit header.
  Dtset <dataset_type>=all input variables for this dataset
  Psps<pseudopotential_type>=Structure gathering info on the pseudopotentials.
  iomode=Input variables specifying the fileformat. (0-->Fortran,3-->netcdf with ETSF-IO format).
  occ(mband*nkpt*nsppol)=The occupation factors.

OUTPUT

  kss_unt=The unit number of the opened file.

SIDE EFFECTS

  The KSS Header is written on file.

PARENTS

      outkss

CHILDREN

      metric,mkffnl,mkkin

SOURCE

123 subroutine write_kss_header(filekss,kss_npw,ishm,nbandksseff,mband,nsym2,symrel2,tnons2,occ,gbig,shlim,&
124 &  crystal,Dtset,Hdr,Psps,iomode,kss_unt)
125 
126 
127 !This section has been created automatically by the script Abilint (TD).
128 !Do not modify the following lines by hand.
129 #undef ABI_FUNC
130 #define ABI_FUNC 'write_kss_header'
131 !End of the abilint section
132 
133  implicit none
134 
135 !Arguments ------------------------------------
136 !scalars
137  integer,intent(in) :: iomode,kss_npw,nbandksseff,ishm,nsym2,mband
138  integer,intent(out) :: kss_unt
139  character(len=fnlen),intent(in) :: filekss
140  type(crystal_t),intent(in) :: crystal
141  type(pseudopotential_type),intent(in) :: Psps
142  type(Hdr_type),intent(in) :: Hdr
143  type(Dataset_type),intent(in) :: Dtset
144 !arrays
145  integer,intent(in) :: symrel2(3,3,nsym2)
146  integer,intent(in) :: gbig(3,kss_npw),shlim(ishm)
147  real(dp),intent(in) :: tnons2(3,nsym2)
148  real(dp),intent(in) :: occ(mband*Dtset%nkpt*Dtset%nsppol)
149 
150 !Local variables-------------------------------
151 !scalars
152  integer :: nspinor,nsppol,nkpt,itypat,ierr
153  integer :: nb,isppol,ik,fform,ii,jj,kk,ig
154  integer :: il,il0,ilmn,in,ind1,ind2
155  character(len=80) :: title
156  character(len=500) :: msg
157  type(hdr_type) :: my_Hdr
158  type(dataset_type) :: Dtset_cpy
159 #ifdef HAVE_NETCDF
160  integer :: ncerr
161 #endif
162 !arrays
163  integer,allocatable :: vkbsign_int(:,:,:)
164  real(dp),allocatable :: vkbsign(:,:)
165 
166 ! *********************************************************************
167 
168  DBG_ENTER("COLL")
169 
170  nsppol = Dtset%nsppol
171  nkpt   = Dtset%nkpt
172  nspinor= Dtset%nspinor
173 
174  write(msg,'(3a)')ch10,' Opening file for KS structure output: ',TRIM(filekss)
175  call wrtout(std_out,msg,'COLL')
176 
177  write(msg,'(a,i6)') ' number of Gamma centered plane waves ',kss_npw
178  call wrtout(std_out,msg,'COLL')
179  call wrtout(ab_out,msg,'COLL')
180  write(msg,'(a,i6)') ' number of Gamma centered shells ',ishm
181  call wrtout(std_out,msg,'COLL')
182  call wrtout(ab_out,msg,'COLL')
183  write(msg,'(a,i6)') ' number of bands ',nbandksseff
184  call wrtout(std_out,msg,'COLL')
185  call wrtout(ab_out,msg,'COLL')
186  write(msg,'(a,i6)') ' maximum angular momentum components ',Psps%mpsang
187  call wrtout(std_out,msg,'COLL')
188  call wrtout(ab_out,msg,'COLL')
189  write(msg,'(a,i2,a)')' number of symmetry operations ',nsym2,' (without inversion)'
190  call wrtout(std_out,msg,'COLL')
191 
192 !Copy the header so that we can change some basic dimensions using the KSS values:
193 !(bantot, npwarr, nband) and the occupation factors
194 
195 !Note that nsym and symrel might have been changed this has to be fixed
196 !carefully in the next patch since in the new implementation symmorphy=0 should be dafault
197  call hdr_copy(Hdr,my_Hdr)
198 
199  my_Hdr%npwarr =kss_npw
200  my_Hdr%nband  =nbandksseff
201  my_hdr%mband = maxval(my_hdr%nband)
202  my_Hdr%bantot =nbandksseff*nkpt*nsppol
203 
204  my_Hdr%istwfk = 1  ! KSS file does not support istwfk/=1 even though the GS run
205                     ! can take advantage of time-reversal symmetry.
206 
207 !Copy the occ number in the new header with correct dimensions
208 !fill with zero the rest since mband can be < nbandksseff
209  !write(std_out,*)associated(my_Hdr%occ)
210  ABI_DEALLOCATE(my_Hdr%occ)
211  ABI_ALLOCATE(my_Hdr%occ,(my_Hdr%bantot))
212  !mband = MAXVAL(Hdr%nband)
213 
214  my_Hdr%occ=zero; nb=MIN(mband,nbandksseff)
215  do isppol=1,nsppol
216    do ik=1,nkpt
217      ind1=1+(ik-1)*nbandksseff+(isppol-1)*nkpt*nbandksseff
218      ind2=1+(ik-1)*mband      +(isppol-1)*nkpt*mband
219      my_Hdr%occ(ind1:ind1+nb-1) = occ(ind2:ind2+nb-1)
220    end do
221  end do
222 
223 !Change dimension in the local Dtset_cpy as well.
224  call dtset_copy(Dtset_cpy, Dtset)
225  Dtset_cpy%mpw   = kss_npw
226  Dtset_cpy%mband = nbandksseff
227 
228  fform=502
229 
230  SELECT CASE (iomode)
231 
232  CASE (IO_MODE_FORTRAN)
233 
234    if (open_file(filekss, msg, newunit=kss_unt, form="unformatted") /= 0) then
235      MSG_ERROR(msg)
236    end if
237 
238    call hdr_fort_write(my_hdr, kss_unt, fform, ierr)
239    ABI_CHECK(ierr == 0, "hdr_Fort_write returned ierr != 0")
240 
241    title='Results from ABINIT code';          write(kss_unt) title(1:80)
242    title='Ab-initio plane waves calculation'; write(kss_unt) title(1:80)
243 
244    write(kss_unt) nsym2,nbandksseff,kss_npw,ishm,Psps%mpsang ! To be modified to deal with more than one projector
245    write(kss_unt) (((symrel2(ii,jj,kk),ii=1,3),jj=1,3),kk=1,nsym2)
246    write(kss_unt) ((tnons2(ii,kk),ii=1,3),kk=1,nsym2)
247    write(kss_unt) ((gbig(ii,ig),ii=1,3),ig=1,kss_npw)
248    write(kss_unt) (shlim(in),in=1,ishm)
249 
250    ! Write vkbsign for NC pseudos with Fortran IO
251    ! MG FIXME: only one projector in each angular channel is treated.
252    ! Moreover the allocation is done in the wrong order for dimensions...
253    ! but if I change this code, compatibility with external codes is broken.
254    if (Psps%usepaw==0) then
255      ABI_ALLOCATE(vkbsign,(Psps%ntypat,Psps%mpsang))
256      vkbsign(:,:)=zero
257      do itypat=1,Psps%ntypat
258        il0=0
259        do ilmn=1,Psps%lmnmax
260          il=1+Psps%indlmn(1,ilmn,itypat)
261          in=Psps%indlmn(3,ilmn,itypat)
262          if (il/=il0 .and. in==1) then
263            il0=il
264            vkbsign(itypat,il)=DSIGN(one,Psps%ekb(ilmn,itypat))
265          end if
266        end do
267      end do
268      write(kss_unt) ((vkbsign(itypat,il),il=1,Psps%mpsang),itypat=1,Psps%ntypat)
269      ABI_DEALLOCATE(vkbsign)
270    end if
271 
272 #ifdef HAVE_NETCDF
273  CASE (IO_MODE_ETSF)
274 
275    ! Create file.
276    NCF_CHECK(nctk_open_create(kss_unt, nctk_ncify(filekss), xmpi_comm_self))
277 
278    ! Add additional info from abinit header.
279    NCF_CHECK(hdr_ncwrite(my_hdr, kss_unt, fform, nc_define=.True.))
280 
281    ! Add info on crystalline structure
282    ! FIXME: Check symmorphi trick and crystal%symrel!
283    ! We currently use the dataset symmetries, as defined in the Hdr structure
284    ! instead of the symmetries recomputed in outkss.
285    NCF_CHECK(crystal_ncwrite(crystal, kss_unt))
286 
287    ! Defined G-vectors and wavefunctions.
288    call wfk_ncdef_dims_vars(kss_unt, my_hdr, fform, iskss=.True.)
289    !call abi_etsf_init(Dtset_cpy, filekss, 4, .false., my_Hdr%lmn_size, Psps, Dummy_wfs)
290 
291    ! If NC pseudos, write vkbsign.
292    ! Here multi-projectors are supported, array is dimensioned according to etsf-io standard.
293    if (psps%usepaw == 0) then
294 
295      ! Define dims and variables needed for KB matrix elements.
296      ncerr = nctk_def_dims(kss_unt, [ &
297        nctkdim_t("max_number_of_angular_momenta", psps%mpsang), &
298        nctkdim_t("max_number_of_projectors", psps%mproj) &
299      ])
300      NCF_CHECK(ncerr)
301 
302      ncerr = nctk_def_arrays(kss_unt, [ &
303        nctkarr_t("kb_formfactor_sign", "int", &
304 &"max_number_of_projectors, max_number_of_angular_momenta, number_of_atom_species"), &
305        nctkarr_t("kb_formfactors", "dp", &
306 &"max_number_of_coefficients, number_of_kpoints, max_number_of_projectors,&
307 &max_number_of_angular_momenta, number_of_atom_species"), &
308        nctkarr_t("kb_formfactor_derivative", "dp", &
309 &"max_number_of_coefficients, number_of_kpoints, max_number_of_projectors,&
310 &max_number_of_angular_momenta, number_of_atom_species") &
311      ])
312      NCF_CHECK(ncerr)
313 
314      ABI_ALLOCATE(vkbsign_int, (psps%mproj, Psps%mpsang, Psps%ntypat))
315      vkbsign_int=0
316      do itypat=1,Psps%ntypat
317        do ilmn=1,Psps%lmnmax
318          il=1+Psps%indlmn(1,ilmn,itypat)
319          in=Psps%indlmn(3,ilmn,itypat)
320          vkbsign_int(in,il,itypat)=NINT(DSIGN(one,Psps%ekb(ilmn,itypat)))
321        end do
322      end do
323 
324      NCF_CHECK(nctk_set_datamode(kss_unt))
325 
326      ! Write KB sign here
327      NCF_CHECK(nf90_put_var(kss_unt, nctk_idname(kss_unt, "kb_formfactor_sign"), vkbsign_int))
328      ABI_DEALLOCATE(vkbsign_int)
329    end if
330 
331    NCF_CHECK(nctk_set_datamode(kss_unt))
332 #endif
333 
334  CASE DEFAULT
335    MSG_ERROR(sjoin("Unsupported value for iomode:", itoa(iomode)))
336  END SELECT
337 
338  call dtset_free(Dtset_cpy)
339  call hdr_free(my_Hdr)
340 
341  DBG_EXIT("COLL")
342 
343 end subroutine write_kss_header

m_io_kss/write_kss_wfgk [ Functions ]

[ Top ] [ m_io_kss ] [ Functions ]

NAME

  write_kss_wfgk

FUNCTION

  Write the Gamma-centered wavefunctions and energies on the KSS file for a single k-point.
  (Only the master node should call this routine).

INPUTS

  kss_unt=The unit number of the file
  ikpt=The index of the k-point
  isppol=The spin index.
  nspinor=number of spinorial components (on current proc)
  kss_npw=Number of planewaves used for the wavefunctions in the KSS files.
  npw_k=Number of plane-waves in the k-centered basis set.
  nbandksseff=Number of bands to be written.
  natom=Number of atoms.
  Psps<Pseudopotential_type>=Structure gathering pseudopotential data.
  kpoint(3)=The k-points in reduced coordinates.
  ene_k(nbandksseff)=Energies at this k-point
  occ_k(nbandksseff)=Occupation factors at this k-point.
  rprimd(3,3)=dimensional primitive translations for real space (bohr).
  gbig(3,kss_npw)=The set of G-vectors for the KSS wavefunctions (Gamma-centered)
  wfg(2,kss_npw*nspinor,nbandksseff)=The wavefunction Fourier coefficients.
  iomode=Input variables specifying the fileformat. (0-->Fortran,3--> netcdf with ETSF-IO format).

OUTPUT

  Only writing.

PARENTS

      outkss

CHILDREN

      metric,mkffnl,mkkin

SOURCE

516 subroutine write_kss_wfgk(kss_unt,ikpt,isppol,kpoint,nspinor,kss_npw,&
517 &          nbandksseff,natom,Psps,ene_k,occ_k,rprimd,gbig,wfg,Cprjnk_k,iomode)
518 
519 
520 !This section has been created automatically by the script Abilint (TD).
521 !Do not modify the following lines by hand.
522 #undef ABI_FUNC
523 #define ABI_FUNC 'write_kss_wfgk'
524 !End of the abilint section
525 
526  implicit none
527 
528 !Arguments ------------------------------------
529 !scalars
530  integer,intent(in) :: ikpt,isppol,iomode,kss_npw,nspinor,kss_unt,nbandksseff
531  integer,intent(in) :: natom
532  type(pseudopotential_type),intent(in) :: Psps
533 !arrays
534  integer,intent(in) :: gbig(3,kss_npw)
535  real(dp),intent(in) :: kpoint(3),rprimd(3,3)
536  real(dp),intent(in) :: ene_k(nbandksseff),occ_k(nbandksseff)
537  real(dp),intent(in) ::  wfg(2,kss_npw*nspinor,nbandksseff)
538  type(pawcprj_type),intent(in) :: Cprjnk_k(natom,nspinor*nbandksseff*Psps%usepaw)
539 
540 !Local variables-------------------------------
541 !scalars
542  integer :: ib,ibsp,ig,ispinor,iatom,ii !,ierr
543 #ifdef HAVE_NETCDF
544  integer :: kg_varid,cg_varid,ncerr
545  character(len=nctk_slen) :: kdep
546 #endif
547 
548 ! *********************************************************************
549 
550  ! Calculate and write KB form factors and derivative at this k-point.
551  if (Psps%usepaw==0) then
552    call write_vkb(kss_unt,ikpt,kpoint,kss_npw,gbig,rprimd,Psps,iomode)
553  end if
554 
555  ! ============================================================
556  ! ==== Write wavefunctions and PAW matrix elements on disk ====
557  ! ============================================================
558  SELECT CASE (iomode)
559 
560  CASE (IO_MODE_FORTRAN)
561    write(kss_unt) (ene_k(ib),ib=1,nbandksseff)
562 
563    ibsp=0
564    do ib=1,nbandksseff
565      write(kss_unt) (wfg(:,ig,ib),ig=1,kss_npw*nspinor)
566      if (Psps%usepaw==1) then ! Remember that cprj are unsorted.
567        do ispinor=1,nspinor
568          ibsp=ibsp+1
569          do iatom=1,natom
570            ii=Cprjnk_k(iatom,ibsp)%nlmn
571            write(kss_unt) (Cprjnk_k(iatom,ibsp)%cp(:,1:ii))
572          end do
573        end do
574      end if
575    end do
576 
577 #ifdef HAVE_NETCDF
578  CASE (IO_MODE_ETSF)
579    if (Psps%usepaw==1) then
580      MSG_WARNING("PAW output with ETSF-IO netcdf: cprj won't be written")
581    end if
582 
583    ! Write G-vectors (gbig because it's not k-dependent)
584    NCF_CHECK(nf90_inq_varid(kss_unt, "reduced_coordinates_of_plane_waves", kg_varid))
585    NCF_CHECK(nf90_get_att(kss_unt, kg_varid, "k_dependent", kdep))
586    if (kdep == "no") then
587      ncerr = nf90_put_var(kss_unt, kg_varid, gbig, start=[1,1], count=[3,kss_npw])
588    else
589      ncerr = nf90_put_var(kss_unt, kg_varid, gbig, start=[1,1,ikpt], count=[3,kss_npw,1])
590    end if
591    NCF_CHECK_MSG(ncerr, "putting gibg")
592 
593    ! Write wavefunctions
594    ! The coefficients_of_wavefunctions on file have shape [cplex, mpw, nspinor, mband, nkpt, nsppol]
595    NCF_CHECK(nf90_inq_varid(kss_unt, "coefficients_of_wavefunctions", cg_varid))
596    ncerr = nf90_put_var(kss_unt, cg_varid, wfg, start=[1,1,1,1,ikpt,isppol], &
597      count=[2,kss_npw,nspinor,nbandksseff,1,1])
598    NCF_CHECK_MSG(ncerr, "putting cg_k")
599 
600    ! Write eigenvalues and occupations
601    NCF_CHECK(nf90_put_var(kss_unt, nctk_idname(kss_unt, "eigenvalues"), ene_k, start=[1,ikpt,isppol]))
602    NCF_CHECK(nf90_put_var(kss_unt, nctk_idname(kss_unt, "occupations"), occ_k, start=[1,ikpt,isppol]))
603 #endif
604 
605  CASE DEFAULT
606    MSG_ERROR(sjoin("Unsupported iomode:", itoa(iomode)))
607  END SELECT
608 
609 end subroutine write_kss_wfgk

m_io_kss/write_vkb [ Functions ]

[ Top ] [ m_io_kss ] [ Functions ]

NAME

  write_vkb

FUNCTION

  Writes the KB form factors and derivates on file for a single k-point.
  Supports plain Fortran IO and netcdf with ETSF-IO format

INPUTS

  kss_unt=The unit number of the file
  ikpt=The index of the k-point
  kpoint(3)=The k-point in reduced coordinates.
  kss_npw=Number of planewaves used for the wavefunctions in the KSS files.
  npw_k=Number of planewaves at this k-point in the k-centered basis set used in abinit (ecut).
  ecut=cutoff energy used in abinit.
  rprimd(3,3)=dimensional primitive translations for real space (bohr).
  Psps<Pseudopotential_type>=Datatype gathering data on the Pseudopotentials.
  iomode=Input variables specifying the fileformat. (0-->Fortran,3-->netcdf with ETSF-IO).
  gbig(3,kss_npw)=Set of G-vectors used in the KSS file.

OUTPUT

  Only writing.

PARENTS

      m_io_kss

CHILDREN

      metric,mkffnl,mkkin

SOURCE

379 subroutine write_vkb(kss_unt,ikpt,kpoint,kss_npw,gbig,rprimd,Psps,iomode)
380 
381 
382 !This section has been created automatically by the script Abilint (TD).
383 !Do not modify the following lines by hand.
384 #undef ABI_FUNC
385 #define ABI_FUNC 'write_vkb'
386 !End of the abilint section
387 
388  implicit none
389 
390 !Arguments ------------------------------------
391 !scalars
392  integer,intent(in) :: ikpt,iomode,kss_npw,kss_unt
393  type(Pseudopotential_type),intent(in) :: Psps
394 !arrays
395  integer,intent(in) :: gbig(3,kss_npw)
396  real(dp),intent(in) :: kpoint(3),rprimd(3,3)
397 
398 !Local variables-------------------------------
399 !scalars
400  integer :: itypat,il,ig,mpsang,ntypat
401 !array
402  real(dp),allocatable :: vkb(:,:,:),vkbd(:,:,:)
403  real(dp),allocatable :: dum_vkbsign(:,:)
404 #ifdef HAVE_NETCDF
405  integer :: ncerr,varid
406  real(dp),allocatable,target :: vkb_tgt(:,:,:,:), vkbd_tgt(:,:,:,:)
407 #endif
408 
409 ! *********************************************************************
410 
411  mpsang = Psps%mpsang; ntypat = Psps%ntypat
412 
413  ABI_ALLOCATE(vkb ,(kss_npw,ntypat,mpsang))
414  ABI_ALLOCATE(vkbd,(kss_npw,ntypat,mpsang))
415  ABI_ALLOCATE(dum_vkbsign,(ntypat,mpsang))
416 
417  call kss_calc_vkb(Psps,kpoint,kss_npw,gbig,rprimd,dum_vkbsign,vkb,vkbd)
418  ABI_DEALLOCATE(dum_vkbsign)
419 
420  SELECT CASE (iomode)
421 
422  CASE (IO_MODE_FORTRAN)
423   do itypat=1,ntypat
424     do il=1,mpsang
425       write(kss_unt) (vkb (ig,itypat,il),ig=1,kss_npw)
426       write(kss_unt) (vkbd(ig,itypat,il),ig=1,kss_npw)
427     end do
428   end do
429 
430 #ifdef HAVE_NETCDF
431  CASE (IO_MODE_ETSF)
432    ABI_ALLOCATE(vkb_tgt ,(kss_npw,1,mpsang,ntypat))
433    ABI_ALLOCATE(vkbd_tgt,(kss_npw,1,mpsang,ntypat))
434    do itypat=1,ntypat
435      do il=1,mpsang
436        do ig=1,kss_npw
437          vkb_tgt (ig,1,il,itypat)=vkb (ig,itypat,il)
438          vkbd_tgt(ig,1,il,itypat)=vkbd(ig,itypat,il)
439        end do
440      end do
441    end do
442 
443    ! FIXME: Multiple projectors
444    !ABI_MALLOC(vkbsign, (psps%lnmax, cryst%ntypat))
445    !ABI_MALLOC(vkb, (npw, psps%lnmax, cryst%ntypat))
446    !ABI_MALLOC(vkbd, (npw, psps%lnmax, cryst%ntypat))
447    !call calc_vkb(cryst,psps,kpoint,npw,npw,gvec,vkbsign,vkb,vkbd)
448    !ABI_FREE(vkbsign)
449    !ABI_FREE(vkb)
450    !ABI_FREE(vkbd)
451 
452    ! The shape of the variable on disk is (Fortran API):
453    !  (max_number_of_coefficients, number_of_kpoints, max_number_of_projectors,
454    !   max_number_of_angular_momenta, number_of_atom_species)
455    varid = nctk_idname(kss_unt, "kb_formfactors")
456    ncerr = nf90_put_var(kss_unt, varid, vkb_tgt, start=[1,ikpt,1,1,1], count=[kss_npw,1,1,mpsang,ntypat])
457    NCF_CHECK(ncerr)
458 
459    varid = nctk_idname(kss_unt, "kb_formfactor_derivative")
460    ncerr = nf90_put_var(kss_unt, varid, vkbd_tgt, start=[1,ikpt,1,1,1], count=[kss_npw,1,1,mpsang,ntypat])
461    NCF_CHECK(ncerr)
462 
463    ABI_DEALLOCATE(vkb_tgt)
464    ABI_DEALLOCATE(vkbd_tgt)
465 #endif
466 
467  CASE DEFAULT
468    MSG_ERROR(sjoin("Unsupported value for iomode:", itoa(iomode)))
469  END SELECT
470 
471  ABI_DEALLOCATE(vkb)
472  ABI_DEALLOCATE(vkbd)
473 
474 end subroutine write_vkb