TABLE OF CONTENTS


ABINIT/m_rttddft_properties [ Modules ]

[ Top ] [ Modules ]

NAME

  m_rttddft_properties

FUNCTION

  Contains most of the subroutines to compute
  properties (energy, occupations, eigenvalues..)

COPYRIGHT

  Copyright (C) 2021-2022 ABINIT group (FB)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

SOURCE

17 #if defined HAVE_CONFIG_H
18 #include "config.h"
19 #endif
20 
21 #include "abi_common.h"
22 
23 module m_rttddft_properties
24 
25  use defs_basis
26  use defs_abitypes,      only: MPI_type
27  use defs_datatypes,     only: pseudopotential_type
28 
29  use m_bandfft_kpt,      only: bandfft_kpt_type
30  use m_cgprj,            only: ctocprj
31  use m_cgtools,          only: dotprod_g
32  use m_dtset,            only: dataset_type
33  use m_energies,         only: energies_type
34  use m_fourier_interpol, only: transgrid
35  use m_hamiltonian,      only: gs_hamiltonian_type
36  use m_mkrho,            only: mkrho
37  use m_nonlop,           only: nonlop
38  use m_pawcprj,          only: pawcprj_type, pawcprj_alloc, pawcprj_get, &
39                              & pawcprj_free, pawcprj_mpi_allgather
40  use m_paw_mkrho,        only: pawmkrho
41  use m_paw_occupancies,  only: pawmkrhoij
42  use m_pawrhoij,         only: pawrhoij_type, pawrhoij_free, &
43                              & pawrhoij_alloc, pawrhoij_inquire_dim
44  use m_profiling_abi,    only: abimem_record
45  use m_rttddft_tdks,     only: tdks_type
46  use m_spacepar,         only: meanvalue_g
47  use m_xmpi,             only: xmpi_sum, xmpi_comm_rank
48 
49  implicit none
50 
51  private

m_rttddft_properties/rttddft_calc_density [ Functions ]

[ Top ] [ m_rttddft_properties ] [ Functions ]

NAME

  rttddft_calc_density

FUNCTION

  Compute electronic density (in 1/bohr^3) from the WF (cg coefficients)

INPUTS

  dtset <type(dataset_type)> = all input variables for this dataset
  mpi_enreg <MPI_type> = MPI-parallelisation information
  psps <type(pseudopotential_type)> = variables related to pseudopotentials
  tdks <type(tdks_type)> = Main RT-TDDFT object

OUTPUT

SOURCE

 82 subroutine rttddft_calc_density(dtset, mpi_enreg, psps, tdks)
 83 
 84  implicit none
 85 
 86  !Arguments ------------------------------------
 87  !scalars
 88  type(tdks_type),            intent(inout) :: tdks
 89  type(dataset_type),         intent(inout) :: dtset
 90  type(MPI_type),             intent(inout) :: mpi_enreg
 91  type(pseudopotential_type), intent(inout) :: psps
 92 
 93  !Local variables-------------------------------
 94  !scalars
 95  integer, parameter          :: cplex=1
 96  integer                     :: cplex_rhoij
 97  integer                     :: ipert, idir
 98  integer                     :: my_natom
 99  integer                     :: nspden_rhoij
100  integer                     :: tim_mkrho
101  real(dp)                    :: compch_fft
102  !arrays
103  real(dp)                    :: qpt(3)
104  real(dp),allocatable        :: rhowfg(:,:), rhowfr(:,:)
105  type(pawrhoij_type),pointer :: pawrhoij_unsym(:)
106 
107 ! ***********************************************************************
108 
109  my_natom=mpi_enreg%my_natom
110 
111  tim_mkrho=1
112 
113  if (psps%usepaw==1) then
114 
115    ABI_MALLOC(rhowfg,(2,dtset%nfft))
116    ABI_MALLOC(rhowfr,(dtset%nfft,dtset%nspden))
117 
118    ! 1-Compute density from WFs (without compensation charge density nhat)
119    call mkrho(tdks%cg,dtset,tdks%gprimd,tdks%irrzon,tdks%kg,tdks%mcg,mpi_enreg, &
120             & tdks%npwarr,tdks%occ0,tdks%paw_dmft,tdks%phnons,rhowfg,rhowfr,    &
121             & tdks%rprimd,tim_mkrho,tdks%ucvol,tdks%wvl%den,tdks%wvl%wfs)
122 
123    ! 2-Compute cprj = <\psi_{n,k}|p_{i,j}>
124    call ctocprj(tdks%atindx,tdks%cg,1,tdks%cprj,tdks%gmet,tdks%gprimd,0,0,0,           &
125               & dtset%istwfk,tdks%kg,dtset%kptns,tdks%mcg,tdks%mcprj,dtset%mgfft,      &
126               & dtset%mkmem,mpi_enreg,psps%mpsang,dtset%mpw,dtset%natom,tdks%nattyp,   &
127               & dtset%nband,dtset%natom,dtset%ngfft,dtset%nkpt,dtset%nloalg,           &
128               & tdks%npwarr,dtset%nspinor,dtset%nsppol,dtset%nsppol,psps%ntypat,       &
129               & dtset%paral_kgb,tdks%ph1d,psps,tdks%rmet,dtset%typat,tdks%ucvol,       &
130               & tdks%unpaw,tdks%xred,tdks%ylm,tdks%ylmgr)
131 
132    !paral atom
133    if (my_natom/=dtset%natom) then
134      ABI_MALLOC(pawrhoij_unsym,(dtset%natom))
135      call pawrhoij_inquire_dim(cplex_rhoij=cplex_rhoij,nspden_rhoij=nspden_rhoij, &
136                              & nspden=dtset%nspden,spnorb=dtset%pawspnorb,        &
137                              & cpxocc=dtset%pawcpxocc)
138      call pawrhoij_alloc(pawrhoij_unsym,cplex_rhoij,nspden_rhoij,dtset%nspinor, &
139                        & dtset%nsppol,dtset%typat,pawtab=tdks%pawtab,use_rhoijp=0)
140    else
141        pawrhoij_unsym => tdks%pawrhoij
142    end if
143 
144    ! 3-Compute pawrhoij = \rho_{i,j} = \sum_{n,k}f_{n,k} \tilde{c}^{i,*}_{n,k} \tilde{c}^{j}_{n,k}
145    call pawmkrhoij(tdks%atindx,tdks%atindx1,tdks%cprj,tdks%dimcprj,dtset%istwfk,    &
146                  & dtset%kptopt,dtset%mband,tdks%mband_cprj,tdks%mcprj,dtset%mkmem, &
147                  & mpi_enreg,dtset%natom,dtset%nband,dtset%nkpt,dtset%nspinor,      &
148                  & dtset%nsppol,tdks%occ0,dtset%paral_kgb,tdks%paw_dmft,            &
149                  & pawrhoij_unsym,tdks%unpaw,dtset%usewvl,dtset%wtk)
150 
151    ! 4-Symetrize rhoij, compute nhat and add it to rhor
152    ! Note pawrhoij_unsym and pawrhoij are the same, which means that pawrhoij
153    ! cannot be distributed over different atomic sites.
154    ipert=0; idir=0; qpt(:)=zero; compch_fft=-1e-5_dp
155    tdks%nhat = zero
156    call pawmkrho(1,compch_fft,cplex,tdks%gprimd,idir,tdks%indsym,ipert,mpi_enreg, &
157                & my_natom,dtset%natom,dtset%nspden,dtset%nsym,dtset%ntypat,       &
158                & dtset%paral_kgb,tdks%pawang,tdks%pawfgr,tdks%pawfgrtab,          &
159                & dtset%pawprtvol,tdks%pawrhoij,pawrhoij_unsym,tdks%pawtab,qpt,    &
160                & rhowfg,rhowfr,tdks%rhor,tdks%rprimd,dtset%symafm,tdks%symrec,    &
161                & dtset%typat,tdks%ucvol,dtset%usewvl,tdks%xred,pawnhat=tdks%nhat, &
162                & rhog=tdks%rhog)
163 
164    ! 6-Take care of kinetic energy density
165    if(dtset%usekden==1)then
166      call mkrho(tdks%cg,dtset,tdks%gprimd,tdks%irrzon,tdks%kg,tdks%mcg,mpi_enreg, &
167               & tdks%npwarr,tdks%occ0,tdks%paw_dmft,tdks%phnons,rhowfg,rhowfr,     &
168               & tdks%rprimd,tim_mkrho,tdks%ucvol,tdks%wvl%den,tdks%wvl%wfs,option=1)
169      !FB: Useful?
170      call transgrid(1,mpi_enreg,dtset%nspden,+1,1,1,dtset%paral_kgb,tdks%pawfgr, &
171                   & rhowfg,tdks%taug,rhowfr,tdks%taur)
172    end if
173 
174    ABI_FREE(rhowfg)
175    ABI_FREE(rhowfr)
176 
177    if (my_natom/=dtset%natom) then
178      call pawrhoij_free(pawrhoij_unsym)
179      ABI_FREE(pawrhoij_unsym)
180    else
181       pawrhoij_unsym => NULL()
182    end if
183 
184  else
185 
186    ! 1-Compute density from WFs
187    call mkrho(tdks%cg,dtset,tdks%gprimd,tdks%irrzon,tdks%kg,tdks%mcg,mpi_enreg,   &
188             & tdks%npwarr,tdks%occ0,tdks%paw_dmft,tdks%phnons,tdks%rhog,tdks%rhor, &
189             & tdks%rprimd,tim_mkrho,tdks%ucvol,tdks%wvl%den,tdks%wvl%wfs)
190    ! 2-Take care of kinetic energy density
191    if(dtset%usekden==1)then
192      call mkrho(tdks%cg,dtset,tdks%gprimd,tdks%irrzon,tdks%kg,tdks%mcg,mpi_enreg,   &
193               & tdks%npwarr,tdks%occ0,tdks%paw_dmft,tdks%phnons,tdks%taug,tdks%taur, &
194               & tdks%rprimd,tim_mkrho,tdks%ucvol,tdks%wvl%den,tdks%wvl%wfs,option=1)
195    end if
196 
197  endif
198 
199 end subroutine rttddft_calc_density

m_rttddft_properties/rttddft_calc_eig [ Functions ]

[ Top ] [ m_rttddft_properties ] [ Functions ]

NAME

  rttddft_calc_eig

FUNCTION

  Computes eigenvalues from cg and ghc = <G|H|C> 
  and gsc = <G|S|C> if paw

INPUTS

  cg <real(2,npw*nspinor*nband)> = the wavefunction coefficients
  ghc <real(2,npw*nspinor*nband)> = <G|H|C>
  istwf_k <integer> = option describing the storage of wfs at k
  nband <integer> = number of bands
  npw <integer> = number of plane waves
  nspinor <integer> = dimension of spinors
  me_g0 <integer> = if set to 1 current proc contains G(0,0,0)
  comm <integer> = MPI communicator
  gsc <real(2,npw*nspinor*nband)> = <G|S|C> (optional - only in PAW)

OUTPUT

  eig <real(nband)> = the eigenvalues

SOURCE

307 subroutine rttddft_calc_eig(cg,eig,ghc,istwf_k,nband,npw,nspinor,me_g0,comm,gsc)
308 
309  implicit none
310 
311  !Arguments ------------------------------------
312  !scalars
313  integer,  intent(in)            :: istwf_k
314  integer,  intent(in)            :: nband
315  integer,  intent(in)            :: npw
316  integer,  intent(in)            :: nspinor
317  integer,  intent(in)            :: me_g0
318  integer,  intent(in)            :: comm
319  !arrays
320  real(dp), intent(in)            :: cg(2,npw*nspinor*nband)
321  real(dp), intent(out)           :: eig(nband)
322  real(dp), intent(in)            :: ghc(2,npw*nspinor*nband)
323  real(dp), intent(in), optional  :: gsc(2,npw*nspinor*nband)
324 
325  !Local variables-------------------------------
326  !scalars
327  integer   :: iband
328  integer   :: shift
329  real(dp)  :: dprod_r, dprod_i
330  !arrays
331 
332 ! ***********************************************************************
333 
334  do iband=1, nband
335     shift = npw*nspinor*(iband-1)
336     !Compute eigenvalues
337     call dotprod_g(eig(iband),dprod_i,istwf_k,npw*nspinor,1, &
338                  & ghc(:, shift+1:shift+npw*nspinor),        &
339                  & cg(:, shift+1:shift+npw*nspinor),         &
340                  & me_g0, comm)
341     if (present(gsc)) then
342        call dotprod_g(dprod_r,dprod_i,istwf_k,npw*nspinor,1, &
343                     & gsc(:, shift+1:shift+npw*nspinor),     &
344                     & cg(:, shift+1:shift+npw*nspinor),      &
345                     & me_g0, comm)
346        eig(iband) = eig(iband)/dprod_r
347     end if
348  end do
349 
350 end subroutine rttddft_calc_eig

m_rttddft_properties/rttddft_calc_energy [ Functions ]

[ Top ] [ m_rttddft_properties ] [ Functions ]

NAME

  rttddft_calc_energy

FUNCTION

  Computes total energy

INPUTS

  dtset <type(dataset_type)> = all input variables for this dataset
  energies <energies_type> = contains the different contribution ot the total energy
  occ <real(nband*nkpt*nsspol)> = occupation numbers at time t

OUTPUT

  etotal <real(dp)> = the total energy

SOURCE

219 subroutine rttddft_calc_etot(dtset, energies, etotal, occ)
220 
221  implicit none
222 
223  !Arguments ------------------------------------
224  !scalars
225  type(dataset_type),  intent(inout) :: dtset
226  type(energies_type), intent(inout) :: energies
227  real(dp),            intent(out)   :: etotal
228  !arrays
229  real(dp),            intent(in)    :: occ(:)
230 
231  !Local variables-------------------------------
232  !scalars
233  real(dp) :: entropy
234  !arrays
235 
236 ! ***********************************************************************
237 
238  ! Compute electronic entropy
239  call rttddft_calc_ent(entropy, dtset, occ)
240 
241 !  When the finite-temperature VG broadening scheme is used,
242 !  the total entropy contribution "tsmear*entropy" has a meaning,
243 !  and gather the two last terms of Eq.8 of VG paper
244 !  Warning : might have to be changed for fixed moment calculations
245  if(dtset%occopt>=3 .and. dtset%occopt<=8) then
246    if (abs(dtset%tphysel) < tol10) then
247       energies%e_entropy = - dtset%tsmear * entropy
248    else
249       energies%e_entropy = - dtset%tphysel * entropy
250    end if
251  else
252     !FB - TODO: Might want to clean this ?!
253    energies%e_entropy = -entropy
254  end if
255 
256  etotal = energies%e_kinetic     &
257       & + energies%e_hartree     &  
258       & + energies%e_xc          &
259       & + energies%e_localpsp    &   
260       & + energies%e_corepsp     &
261       & + energies%e_entropy     &
262       & + energies%e_ewald       &
263       & + energies%e_vdw_dftd    &
264       & + energies%e_nlpsp_vfock &
265       & + energies%e_paw         
266 !FB: @MT Should one add the last e_paw contribution or not?
267 !FB: Seeems like all the other contributions are not relevant here @MT?
268 !     & + energies%e_chempot     &
269 !     & + energies%e_elecfield   &  
270 !     & + energies%e_magfield    &
271 !     & + energies%e_nucdip      &
272 !     & + energies%e_hybcomp_E0  &
273 !     & - energies%e_hybcomp_v0  &
274 !     & + energies%e_hybcomp_v   &
275 !     & + energies%e_constrained_dft
276 
277 !if (psps%usepaw==0) etotal = etotal + energies%e_nlpsp_vfock - energies%e_fock0
278 !if (psps%usepaw==1) etotal = etotal + energies%e_paw + energies%e_fock
279 
280 end subroutine rttddft_calc_etot

m_rttddft_properties/rttddft_calc_enl [ Functions ]

[ Top ] [ m_rttddft_properties ] [ Functions ]

NAME

  rttddft_calc_enl

FUNCTION

  Computes the NL part of energy in NC case

INPUTS

  cg <real(2,npw*nspinor*nband)> = the wavefunction coefficients
  ham_k <gs_hamiltonian_type> = hamiltonian at point k
  nband <integer> = number of bands
  npw <integer> = number of plane waves
  nspinor <integer> = dimension of spinors
  mpi_enreg <MPI_type> = MPI-parallelisation information

OUTPUT

  enl <real(nband)> = the non local part of the energy in NC case

SIDE EFFECTS

SOURCE

466 subroutine rttddft_calc_enl(cg,enl,ham_k,nband,npw,nspinor,mpi_enreg)
467 
468  implicit none
469 
470  !Arguments ------------------------------------
471  !scalars
472  integer,                   intent(in)    :: nband
473  integer,                   intent(in)    :: npw
474  integer,                   intent(in)    :: nspinor
475  type(gs_hamiltonian_type), intent(in)    :: ham_k
476  type(MPI_type),            intent(in)    :: mpi_enreg
477  !arrays
478  real(dp),                  intent(inout) :: cg(2,npw*nspinor*nband)
479  real(dp),                  intent(out)   :: enl(nband)
480 
481  !Local variables-------------------------------
482  !scalars
483  integer, parameter              :: choice=1
484  integer, parameter              :: cpopt=-1
485  integer, parameter              :: paw_opt=0
486  integer, parameter              :: signs=1
487  integer, parameter              :: tim_getghc = 5
488  !arrays
489  type(pawcprj_type), allocatable :: cprj_dummy(:,:)
490  real(dp),           allocatable :: eig_dummy(:)
491  real(dp),           allocatable :: gvnlxc_dummy(:,:)
492  real(dp),           allocatable :: gsc_dummy(:,:)
493 ! ***********************************************************************
494 
495  ABI_MALLOC(cprj_dummy,(ham_k%natom,0))
496  ABI_MALLOC(gsc_dummy,(0,0))
497  ABI_MALLOC(gvnlxc_dummy, (0, 0))
498  ABI_MALLOC(eig_dummy,(nband))
499  call nonlop(choice,cpopt,cprj_dummy,enl,ham_k,0,eig_dummy,mpi_enreg,nband, &
500             & 1,paw_opt,signs,gsc_dummy,tim_getghc,cg,gvnlxc_dummy)
501  ABI_FREE(cprj_dummy)
502  ABI_FREE(gsc_dummy)
503  ABI_FREE(eig_dummy)
504  ABI_FREE(gvnlxc_dummy)
505 
506 end subroutine rttddft_calc_enl

m_rttddft_properties/rttddft_calc_ent [ Functions ]

[ Top ] [ m_rttddft_properties ] [ Functions ]

NAME

  rttddft_calc_ent

FUNCTION

  Computes electronic entropy from occupation numbers f_{nk}
  S = -2 \sum_k \sum_n w(k) [ f_{nk}*ln(f_{nk}) + (1-f_{nk})*ln(1-f_{nk}) ]

INPUTS

  dtset <type(dataset_type)> = all input variables for this dataset
  occ <real(nband*nkpt*nsspol)> = occupation numbers at time t

OUTPUT

  entropy <real>

SOURCE

704 subroutine rttddft_calc_ent(entropy,dtset,occ)
705 
706  implicit none
707 
708  !Arguments ------------------------------------
709  !scalars
710  real(dp),           intent(out)   :: entropy
711  type(dataset_type), intent(inout) :: dtset
712  !arrays
713  real(dp),           intent(in)    :: occ(:)
714 
715  !Local variables-------------------------------
716  !scalars
717  integer  :: band_index
718  integer  :: iband, isppol, ikpt
719  integer  :: nband_k
720  real(dp) :: fnk
721 
722 ! ***********************************************************************
723 
724  entropy = zero
725  band_index=0
726 
727  do isppol = 1, dtset%nsppol
728    do ikpt = 1, dtset%nkpt
729       nband_k = dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
730       do iband = 1, nband_k
731          fnk = occ(iband+band_index)*0.5_dp
732          if( fnk>tol16 .and. (one-fnk)>tol16 ) then
733             entropy = entropy - two*dtset%wtk(ikpt)*(fnk*log(fnk)+(one-fnk)*log(one-fnk))
734          end if
735       end do
736       band_index = band_index + nband_k
737    end do
738  end do
739 
740 end subroutine rttddft_calc_ent

m_rttddft_properties/rttddft_calc_kin [ Functions ]

[ Top ] [ m_rttddft_properties ] [ Functions ]

NAME

  rttddft_calc_kin

FUNCTION

  Computes the NL part of energy in NC case

INPUTS

  cg <real(2,npw*nspinor*nband)> = the wavefunction coefficients
  dtset <type(dataset_type)> = all input variables for this dataset
  ham_k <gs_hamiltonian_type> = hamiltonian at point k
  nband <integer> = number of bands
  npw <integer> = number of plane waves
  nspinor <integer> = dimension of spinors
  occ0 <real(nband)> = initial occupations
  wk <real> = weight of associated kpt
  mpi_enreg <MPI_type> = MPI-parallelisation information
  bandfft_kpt <bandfft_kpt_type> = additional info on parallelisation
   for the associated kpt

OUTPUT

  kin <real(nband)> = the non local part of the energy in NC case

SIDE EFFECTS

SOURCE

380 subroutine rttddft_calc_kin(kin,cg,dtset,ham_k,nband,npw,nspinor,occ0,wk,mpi_enreg,bandfft)
381 
382  implicit none
383 
384  !Arguments ------------------------------------
385  !scalars
386  integer,                   intent(in)    :: nband
387  integer,                   intent(in)    :: npw
388  integer,                   intent(in)    :: nspinor
389  real(dp),                  intent(in)    :: wk
390  type(dataset_type),        intent(inout) :: dtset
391  type(gs_hamiltonian_type), intent(in)    :: ham_k
392  type(MPI_type),            intent(in)    :: mpi_enreg
393  type(bandfft_kpt_type),    intent(in)    :: bandfft
394  !arrays
395  real(dp),                  intent(in)    :: cg(2,npw*nspinor*nband)
396  real(dp),                  intent(inout) :: kin
397  real(dp),                  intent(in)    :: occ0(nband)
398 
399  !Local variables-------------------------------
400  !scalars
401  integer  :: displ
402  integer  :: iband, ipw
403  integer  :: jpw
404  integer  :: me_bandfft
405  integer  :: shift
406  real(dp) :: ar
407  !arrays
408 
409 ! ***********************************************************************
410 
411  if (dtset%paral_kgb /= 1) then 
412    displ = 0
413  else
414    me_bandfft = xmpi_comm_rank(mpi_enreg%comm_bandspinorfft) 
415    displ = bandfft%rdispls(me_bandfft+1)
416  end if
417 
418  do iband=1, nband
419    if (abs(occ0(iband))>tol8) then 
420       shift = npw*nspinor*(iband-1)
421       !FB: meanvalue_g does the mpi_sum over the bands inside, that's not very efficient since 
422       !FB: we could do it only once at the end
423       !FB: From Lucas: meanvalue_g seems slow
424       !FB: It maybe useful not to use meanvalue_g at all here
425       !call meanvalue_g(ar,ham_k%kinpw_k(1+displ:displ+npw*nspinor),0,1,mpi_enreg,npw,nspinor, &
426       !               & cg(:,1+shift:shift+npw*nspinor),cg(:,1+shift:shift+npw*nspinor),0)
427       ar = zero
428       do ipw = 1, npw
429          ar = ar + ham_k%kinpw_k(displ+ipw)*(cg(1,shift+ipw)*cg(1,shift+ipw)+cg(2,shift+ipw)*cg(2,shift+ipw))
430       end do
431       if(nspinor==2)then
432          do ipw = 1+npw, 2*npw
433             jpw = ipw - npw
434             ar = ar + ham_k%kinpw_k(displ+jpw)*(cg(1,shift+ipw)*cg(1,shift+ipw)+cg(2,shift+ipw)*cg(2,shift+ipw))
435          end do
436       end if
437       kin = kin + wk*occ0(iband)*ar
438    end if
439  end do
440 
441 end subroutine rttddft_calc_kin

m_rttddft_properties/rttddft_calc_occ [ Functions ]

[ Top ] [ m_rttddft_properties ] [ Functions ]

NAME

  rttddft_calc_occ

FUNCTION

  Computes occupations at time t from cg(t), cg0 and occ0
  In NC:
    f_{n,k}(t) = \sum_{m} f_{m,k}(0) <\phi_m(0)|\phi_n(t)>
  In PAW:
    f_{n,k}(t) = \sum_{m} f_{m,k}(0) <\phi_m(0)|S|\phi_n(t)>
               = \sum_{m} f_{m,k}(0) [ <\phi_m(0)|\phi_n(t)> +
                 \sum_{i} <\phi_m(0)|p_{i}>\sum_jS_{i,j}<p_{j}|\phi_n(t)> ]

INPUTS

  cg <real(2,npw*nspinor*nband)> = the wavefunction coefficients
  cg0 <real(2,npw*nspinor*nband)> = the initial wavefunction coefficients
  dtset <type(dataset_type)> = all input variables for this dataset
  ham_k <gs_hamiltonian_type> = hamiltonian at point k
  ikpt <integer> = indice of the considered k-point
  ibg <integer> = indice of the considered k-point for cprj
  isppol <integer> = indice of the considered spin-polarization
  mpi_enreg <MPI_type> = MPI-parallelisation information
  nband_k <integer> = number of bands
  npw_k <integer> = number of plane waves
  nspinor <integer> = dimension of spinors
  occ0 <real(nband)> = initial occupations
  tdks <type(tdks_type)> = Main RT-TDDFT object

OUTPUT

  occ <real(nband)> = the occupations at time t

SOURCE

542 subroutine rttddft_calc_occ(cg,cg0,dtset,ham_k,ikpt,ibg,isppol,mpi_enreg,nband_k,npw_k,nspinor,occ,occ0,tdks)
543 
544  implicit none
545 
546  !Arguments ------------------------------------
547  !scalars
548  integer,                   intent(in)    :: ikpt
549  integer,                   intent(in)    :: ibg
550  integer,                   intent(in)    :: isppol
551  integer,                   intent(in)    :: nband_k
552  integer,                   intent(in)    :: npw_k
553  integer,                   intent(in)    :: nspinor
554  type(dataset_type),        intent(inout) :: dtset
555  type(gs_hamiltonian_type), intent(inout) :: ham_k
556  type(MPI_type),            intent(inout) :: mpi_enreg
557  type(tdks_type), target,   intent(inout) :: tdks
558  !arrays
559  real(dp),                  intent(inout) :: cg(2,npw_k*nspinor*nband_k)
560  real(dp),                  intent(inout) :: cg0(2,npw_k*nspinor*nband_k)
561  real(dp),                  intent(out)   :: occ(nband_k)
562  real(dp),                  intent(in)    :: occ0(nband_k)
563 
564  !Local variables-------------------------------
565  !scalars
566  integer                     :: iband, jband, ierr, ind
567  integer                     :: nband_cprj_k
568  integer                     :: natom
569  integer                     :: shift
570  !parameters for nonlop
571  integer, parameter          :: cpopt = 2
572  integer, parameter          :: choice = 1
573  integer, parameter          :: idir = 1
574  integer, parameter          :: nnlout = 1
575  integer, parameter          :: paw_opt = 3
576  integer, parameter          :: signs = 1
577  integer, parameter          :: tim_nonlop = 16
578  logical                     :: cprj_paral_band
579  !arrays
580  real(dp)                    :: csc(2*nband_k)
581  real(dp),allocatable        :: gsc(:,:),gvnlxc(:,:)
582  real(dp),allocatable        :: enlout(:),enlout_im(:)
583  type(pawcprj_type), pointer :: cprj(:,:), cprj_k(:,:)
584  type(pawcprj_type), pointer :: cprj0(:,:), cprj0_k(:,:)
585 
586 ! ***********************************************************************
587 
588  !Prepare cprj in PAW case
589  cprj_paral_band=.false.
590  if (ham_k%usepaw == 1) then
591     ! Determine if cprj datastructure is distributed over bands
592     cprj_paral_band=(tdks%mband_cprj<dtset%mband)
593     nband_cprj_k=nband_k; if (cprj_paral_band) nband_cprj_k=mpi_enreg%bandpp
594     natom = dtset%natom
595     ! Extract the right cprj for this k-point
596     if (dtset%mkmem*dtset%nsppol/=1) then
597        ABI_MALLOC(cprj_k,(natom,nspinor*nband_cprj_k))
598        ABI_MALLOC(cprj0_k,(natom,nspinor*nband_cprj_k))
599        call pawcprj_alloc(cprj_k,0,tdks%dimcprj)
600        call pawcprj_alloc(cprj0_k,0,tdks%dimcprj)
601        call pawcprj_get(tdks%atindx1,cprj_k,tdks%cprj,natom,1,ibg,ikpt,0,isppol,     &
602                       & tdks%mband_cprj,dtset%mkmem,natom,nband_cprj_k,nband_cprj_k, &
603                       & nspinor,dtset%nsppol,tdks%unpaw,mpicomm=mpi_enreg%comm_kpt,  &
604                       & proc_distrb=mpi_enreg%proc_distrb)
605        call pawcprj_get(tdks%atindx1,cprj0_k,tdks%cprj0,natom,1,ibg,ikpt,0,isppol,   &
606                       & tdks%mband_cprj,dtset%mkmem,natom,nband_cprj_k,nband_cprj_k, &
607                       & nspinor,dtset%nsppol,tdks%unpaw,mpicomm=mpi_enreg%comm_kpt,  &
608                       & proc_distrb=mpi_enreg%proc_distrb)
609     else
610        cprj_k => tdks%cprj
611        cprj0_k => tdks%cprj0
612     end if
613     ! If cprj are distributed over bands, gather them (because we need to mix bands)
614     if (cprj_paral_band) then
615        ABI_MALLOC(cprj,(natom,nspinor*nband_k))
616        ABI_MALLOC(cprj0,(natom,nspinor*nband_k))
617        call pawcprj_alloc(cprj,0,tdks%dimcprj)
618        call pawcprj_alloc(cprj0,0,tdks%dimcprj)
619        call pawcprj_mpi_allgather(cprj_k,cprj,natom,nspinor*nband_cprj_k,mpi_enreg%bandpp, &
620                                 & tdks%dimcprj,0,mpi_enreg%nproc_band,mpi_enreg%comm_band, &
621                                 & ierr,rank_ordered=.false.)
622        call pawcprj_mpi_allgather(cprj0_k,cprj0,natom,nspinor*nband_cprj_k,mpi_enreg%bandpp, &
623                                 & tdks%dimcprj,0,mpi_enreg%nproc_band,mpi_enreg%comm_band,   &
624                                 & ierr,rank_ordered=.false.)
625     else
626        cprj => cprj_k
627        cprj0 => cprj0_k
628     end if
629 
630    !allocate necessary arrays for nonlop
631    ABI_MALLOC(gsc,(0,0))
632    ABI_MALLOC(gvnlxc,(0,0))
633    ABI_MALLOC(enlout,(nband_k))
634    ABI_MALLOC(enlout_im,(nband_k))
635  end if
636 
637  csc = zero
638  do iband = 1, nband_k
639    !* 1 - Compute csc = <cg0|cg>
640    shift = npw_k*nspinor*(iband-1)
641    call zgemv('C',npw_k*nspinor,nband_k,cone,cg0,npw_k*nspinor, &
642             & cg(:,shift+1:shift+npw_k*nspinor),1,czero,csc,1)
643    !If band parallel then reduce csc
644    if (mpi_enreg%nproc_band > 1) then
645       call xmpi_sum(csc,mpi_enreg%comm_bandfft,ierr)
646    end if
647    !* 2 - If PAW, add the additional term \sum_i cprj_i \sum_j S_{i,j} cprj_j
648    if (ham_k%usepaw == 1) then
649       call nonlop(choice,cpopt,cprj(:,iband:iband+(nspinor-1)),enlout,     &
650                 & ham_k,idir,(/zero/),mpi_enreg,1,nnlout,paw_opt,signs,    &
651                 & gsc,tim_nonlop,cg(:,shift+1:shift+npw_k*nspinor),gvnlxc, &
652                 & cprjin_left=cprj0,enlout_im=enlout_im,ndat_left=nband_k)
653       do jband = 1, nband_k
654          csc(2*jband-1) = csc(2*jband-1) + enlout(jband)
655          csc(2*jband)   = csc(2*jband)   + enlout_im(jband)
656       end do
657    end if
658    !* 3 - Calc occupations from csc and occ0
659    do jband = 1, nband_k
660       ind = 2*jband
661       occ(iband) = occ(iband) + occ0(jband)*(csc(ind-1)**2+csc(ind)**2)
662    end do
663  end do
664 
665  if (ham_k%usepaw == 1) then
666    ABI_FREE(gsc)
667    ABI_FREE(gvnlxc)
668    ABI_FREE(enlout)
669    ABI_FREE(enlout_im)
670    if (cprj_paral_band) then
671       call pawcprj_free(cprj)
672       ABI_FREE(cprj)
673       call pawcprj_free(cprj0)
674       ABI_FREE(cprj0)
675    end if
676    if (dtset%mkmem*dtset%nsppol/=1) then
677       call pawcprj_free(cprj_k)
678       ABI_FREE(cprj_k)
679       call pawcprj_free(cprj0_k)
680       ABI_FREE(cprj0_k)
681    end if
682  end if
683 
684 end subroutine rttddft_calc_occ