TABLE OF CONTENTS
- ABINIT/m_rttddft_properties
- m_rttddft_properties/rttddft_calc_density
- m_rttddft_properties/rttddft_calc_eig
- m_rttddft_properties/rttddft_calc_energy
- m_rttddft_properties/rttddft_calc_enl
- m_rttddft_properties/rttddft_calc_ent
- m_rttddft_properties/rttddft_calc_kin
- m_rttddft_properties/rttddft_calc_occ
ABINIT/m_rttddft_properties [ 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