TABLE OF CONTENTS
ABINIT/m_rttddft [ Modules ]
NAME
m_rttddft
FUNCTION
Contains various subroutines used in RT-TDDFT
COPYRIGHT
Copyright (C) 2021-2024 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
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_rttddft 23 24 use defs_basis 25 use defs_abitypes, only: MPI_type 26 use defs_datatypes, only: pseudopotential_type 27 28 use m_dtset, only: dataset_type 29 use m_energies, only: energies_type 30 use m_hamiltonian, only: init_hamiltonian, gs_hamiltonian_type 31 use m_kg, only: getph 32 use m_nonlop, only: nonlop 33 use m_paw_an, only: paw_an_reset_flags 34 use m_paw_correlations, only: setrhoijpbe0 35 use m_paw_denpot, only: pawdenpot 36 use m_pawdij, only: pawdij, symdij 37 use m_paw_ij, only: paw_ij_reset_flags 38 use m_paw_nhat, only: nhatgrid 39 use m_paw_tools, only: chkpawovlp 40 use m_profiling_abi, only: abimem_record 41 use m_rttddft_tdks, only: tdks_type 42 use m_specialmsg, only: wrtout 43 use m_setvtr, only: setvtr 44 use m_xmpi, only: xmpi_paral 45 46 implicit none 47 48 private
m_rttddft/rttddft_init_hamiltonian [ Functions ]
[ Top ] [ m_rttddft ] [ Functions ]
NAME
rttddft_init_hamiltonian
FUNCTION
Init/Update various quantities in order to set up the Hamiltonian
INPUTS
dtset <type(dataset_type)> = all input variables for this dataset energies <energies_type> = contains various contribution to the energy istep <integer> = step number mpi_enreg <MPI_type> = MPI-parallelisation information psps <type(pseudopotential_type)> = variables related to pseudopotentials tdks <type(tdks_type)> = Main RT-TDDFT object
OUTPUT
gs_hamk <type(gs_hamiltonian_type)> = Hamiltonian object
SOURCE
188 subroutine rttddft_init_hamiltonian(dtset, energies, gs_hamk, istep, mpi_enreg, psps, tdks) 189 190 implicit none 191 192 !Arguments ------------------------------------ 193 !scalars 194 integer, intent(in) :: istep 195 type(dataset_type), intent(inout) :: dtset 196 type(energies_type), intent(inout) :: energies 197 type(gs_hamiltonian_type), intent(out) :: gs_hamk 198 type(MPI_type), intent(inout) :: mpi_enreg 199 type(pseudopotential_type), intent(inout) :: psps 200 type(tdks_type), intent(inout) :: tdks 201 202 !Local variables------------------------------- 203 !scalars 204 character(len=500) :: msg 205 integer :: comm 206 integer,parameter :: cplex=1 207 integer,parameter :: ipert=0 208 integer :: initialized0 209 integer :: istep_mix 210 integer :: moved_atm_inside, moved_rhor 211 integer :: my_natom 212 integer :: nfftotf 213 integer :: nzlmopt 214 integer :: optene 215 integer :: option 216 integer :: nkxc, n1xccc, n3xccc 217 integer :: usecprj_local 218 logical :: calc_ewald 219 logical :: tfw_activated 220 real(dp) :: compch_sph 221 real(dp) :: vxcavg 222 !arrays 223 real(dp),allocatable :: grchempottn(:,:) 224 real(dp),allocatable :: grewtn(:,:) 225 real(dp),parameter :: k0(3)=(/zero,zero,zero/) 226 real(dp),allocatable :: kxc(:,:) 227 real(dp) :: strsxc(6) 228 real(dp) :: vpotzero(2) 229 230 ! *********************************************************************** 231 232 my_natom=mpi_enreg%my_natom 233 234 !** Set up the potential (calls setvtr) 235 !** The following steps have been gathered in the setvtr routine: 236 !** - get Ewald energy and Ewald forces 237 !** - compute local ionic pseudopotential vpsp 238 !** - possibly compute 3D core electron density xccc3d 239 !** - possibly compute 3D core kinetic energy density 240 !** - possibly compute vxc and vhartr 241 !** - set up vtrial 242 !** Only the local part of the potential is computed here 243 !FB: @MT Are the values of moved_atm_inside and moved_rhor correct? 244 optene = 4; nkxc=0; moved_atm_inside=0; moved_rhor=1 245 n1xccc=0;if (psps%n1xccc/=0) n1xccc=psps%n1xccc 246 n3xccc=0;if (psps%n1xccc/=0) n3xccc=tdks%pawfgr%nfft 247 strsxc(:)=zero 248 !FB: tfw_activated is a save variable in scfcv, should maybe check where it appears again 249 tfw_activated=.false. 250 if (dtset%tfkinfunc==12) tfw_activated=.true. 251 ABI_MALLOC(grchempottn,(3,dtset%natom)) 252 ABI_MALLOC(grewtn,(3,dtset%natom)) 253 ABI_MALLOC(kxc,(tdks%pawfgr%nfft,nkxc)) 254 calc_ewald = .false. 255 if (dtset%ionmov/=0 .or. istep == tdks%first_step) calc_ewald=.true. 256 !FB: Should we also add an option to avoid recomputing xccc3d as for Ewald? 257 call setvtr(tdks%atindx1,dtset,energies,tdks%gmet,tdks%gprimd,grchempottn, & 258 & grewtn,tdks%grvdw,tdks%gsqcut,istep,kxc,tdks%pawfgr%mgfft, & 259 & moved_atm_inside,moved_rhor,mpi_enreg,tdks%nattyp,tdks%pawfgr%nfft, & 260 & tdks%pawfgr%ngfft,tdks%ngrvdw,tdks%nhat,tdks%nhatgr,tdks%nhatgrdim, & 261 & nkxc,psps%ntypat,n1xccc,n3xccc,optene,tdks%pawang,tdks%pawrad, & 262 & tdks%pawrhoij,tdks%pawtab,tdks%ph1df,psps,tdks%rhog,tdks%rhor, & 263 & tdks%rmet,tdks%rprimd,strsxc,tdks%ucvol,tdks%usexcnhat,tdks%vhartr, & 264 & tdks%vpsp,tdks%vtrial,tdks%vxc,vxcavg,tdks%wvl,tdks%xccc3d,tdks%xred, & 265 & taur=tdks%taur,vxc_hybcomp=tdks%vxc_hybcomp,vxctau=tdks%vxctau, & 266 & add_tfw=tfw_activated,xcctau3d=tdks%xcctau3d,calc_ewald=calc_ewald) 267 ABI_FREE(grchempottn) 268 ABI_FREE(grewtn) 269 ABI_FREE(kxc) 270 271 ! set the zero of the potentials here 272 if(dtset%usepotzero==2) tdks%vpsp(:) = tdks%vpsp(:) + tdks%ecore / ( tdks%zion * tdks%ucvol ) 273 274 !** Update PAW quantities 275 !** Compute energies and potentials in the augmentation regions (spheres) 276 !** and pseudopotential strengths (Dij quantities) 277 if (psps%usepaw==1)then 278 !** Local exact exch.: impose occ. matrix if required 279 if (dtset%useexexch/=0) then 280 if (xmpi_paral==1.and.mpi_enreg%paral_hf==1) then 281 comm=mpi_enreg%comm_kpt 282 else 283 comm=mpi_enreg%comm_cell 284 end if 285 istep_mix=1; initialized0=0 286 call setrhoijpbe0(dtset,initialized0,istep,istep_mix, & 287 & comm,my_natom,dtset%natom,dtset%ntypat,tdks%pawrhoij,tdks%pawtab, & 288 & dtset%typat,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 289 end if 290 291 !** Computation of on-site densities/potentials/energies 292 !** Force the recomputation of on-site potentials and Dij 293 call paw_an_reset_flags(tdks%paw_an) 294 !FB: @MT Changed self_consistent to false here. Is this right? 295 call paw_ij_reset_flags(tdks%paw_ij,self_consistent=.false.) 296 option=0; compch_sph=-1.d5; nzlmopt=0 297 call pawdenpot(compch_sph,energies%e_paw,energies%e_pawdc,ipert, & 298 & dtset%ixc,my_natom,dtset%natom,dtset%nspden,psps%ntypat, & 299 & dtset%nucdipmom,nzlmopt,option,tdks%paw_an,tdks%paw_an, & 300 & tdks%paw_ij,tdks%pawang,dtset%pawprtvol,tdks%pawrad, & 301 & tdks%pawrhoij,dtset%pawspnorb,tdks%pawtab,dtset%pawxcdev, & 302 & dtset%spnorbscl,dtset%xclevel,dtset%xc_denpos, & 303 & dtset%xc_taupos,tdks%ucvol, & 304 & psps%znuclpsp,comm_atom=mpi_enreg%comm_atom, & 305 & mpi_atmtab=mpi_enreg%my_atmtab,vpotzero=vpotzero) 306 !Correct the average potential with the calculated constant vpotzero 307 !Correct the total energies accordingly 308 !vpotzero(1) = -beta/ucvol 309 !vpotzero(2) = -1/ucvol sum_ij rho_ij gamma_ij 310 write(msg,'(a,f14.6,2x,f14.6)') & 311 & ' average electrostatic smooth potential [Ha] , [eV]', & 312 & SUM(vpotzero(:)),SUM(vpotzero(:))*Ha_eV 313 call wrtout(std_out,msg,'COLL') 314 tdks%vtrial(:,:)=tdks%vtrial(:,:)+SUM(vpotzero(:)) 315 if(option/=1)then 316 !Fix the direct total energy (non-zero only for charged systems) 317 energies%e_paw=energies%e_paw-SUM(vpotzero(:))*dtset%cellcharge(1) 318 !Fix the double counting total energy accordingly (for both charged AND 319 !neutral systems) 320 energies%e_pawdc=energies%e_pawdc-SUM(vpotzero(:))*tdks%zion+ & 321 & vpotzero(2)*dtset%cellcharge(1) 322 end if 323 324 !** Dij computation 325 !FB: @MT fatvshift? 326 nfftotf=tdks%pawfgr%ngfft(1)*tdks%pawfgr%ngfft(2)*tdks%pawfgr%ngfft(3) 327 call pawdij(cplex,dtset%enunit,tdks%gprimd,ipert,my_natom,dtset%natom, & 328 & tdks%pawfgr%nfft,nfftotf,dtset%nspden,psps%ntypat,tdks%paw_an, & 329 & tdks%paw_ij,tdks%pawang,tdks%pawfgrtab,dtset%pawprtvol,tdks%pawrad, & 330 & tdks%pawrhoij,dtset%pawspnorb,tdks%pawtab,dtset%pawxcdev,k0, & 331 & dtset%spnorbscl,tdks%ucvol,dtset%cellcharge(1),tdks%vtrial, & 332 & tdks%vxc,tdks%xred,natvshift=dtset%natvshift,atvshift=dtset%atvshift, & 333 & fatvshift=one,comm_atom=mpi_enreg%comm_atom, & 334 & mpi_atmtab=mpi_enreg%my_atmtab,mpi_comm_grid=mpi_enreg%comm_fft, & 335 & nucdipmom=dtset%nucdipmom) 336 337 !Symetrize Dij 338 call symdij(tdks%gprimd,tdks%indsym,ipert,my_natom,dtset%natom,dtset%nsym, & 339 & psps%ntypat,0,tdks%paw_ij,tdks%pawang,dtset%pawprtvol, & 340 & tdks%pawtab,tdks%rprimd,dtset%symafm,tdks%symrec, & 341 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab) 342 end if 343 344 !** Initialize most of the Hamiltonian 345 !** Allocate all arrays and initialize quantities that do not depend on k and spin. 346 !FB: Should recompute cprj if ions have moved right? 347 usecprj_local=0; if (psps%usepaw==1 .and. dtset%ionmov==0) usecprj_local=1 348 call init_hamiltonian(gs_hamk,psps,tdks%pawtab,dtset%nspinor,dtset%nsppol,dtset%nspden,dtset%natom,dtset%typat, & 349 & tdks%xred,dtset%nfft,dtset%mgfft,dtset%ngfft,tdks%rprimd,dtset%nloalg,paw_ij=tdks%paw_ij, & 350 & ph1d=tdks%ph1d,usecprj=usecprj_local,comm_atom=mpi_enreg%comm_atom, & 351 & mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab,nucdipmom=dtset%nucdipmom, & 352 & gpu_option=dtset%gpu_option) 353 354 end subroutine rttddft_init_hamiltonian
m_rttddft/rttddft_setup_ele_step [ Functions ]
[ Top ] [ m_rttddft ] [ Functions ]
NAME
rttddft_setup_ele_step
FUNCTION
Init/Update various quantities needed before performing propagation of KS orbitals
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
76 subroutine rttddft_setup_ele_step(dtset, mpi_enreg, psps, tdks) 77 78 implicit none 79 80 !Arguments ------------------------------------ 81 !scalars 82 type(dataset_type), intent(inout) :: dtset 83 type(MPI_type), intent(inout) :: mpi_enreg 84 type(pseudopotential_type), intent(inout) :: psps 85 type(tdks_type), intent(inout) :: tdks 86 87 !Local variables------------------------------- 88 !scalars 89 integer :: forces_needed 90 integer :: my_natom 91 integer :: optcut, optgr0, optgr1, optgr2, optrad 92 integer :: stress_needed 93 !arrays 94 !real(dp),parameter :: k0(3)=(/zero,zero,zero/) 95 96 ! *********************************************************************** 97 98 my_natom=mpi_enreg%my_natom 99 100 !** Update various quantities that needs to be changed 101 !** after a change of xred during the nuclear step 102 103 !Compute large sphere G^2 cut-off (gsqcut) and box / sphere ratio 104 !FB: @MT Probably not needed? The box didn't change only nuclear pos.. 105 !if (psps%usepaw==1) then 106 ! call getcut(tdks%boxcut,dtset%pawecutdg,tdks%gmet,tdks%gsqcut,dtset%iboxcut, & 107 ! & std_out,k0,tdks%pawfgr%ngfft) 108 !else 109 ! call getcut(tdks%boxcut,dtset%ecut,tdks%gmet,tdks%gsqcut,dtset%iboxcut, & 110 ! & std_out,k0,tdks%pawfgr%ngfft) 111 !end if 112 113 !Compute structure factor phases (exp(2Pi i G.xred)) on coarse and fine grid 114 call getph(tdks%atindx,dtset%natom,tdks%pawfgr%ngfftc(1),tdks%pawfgr%ngfftc(2), & 115 & tdks%pawfgr%ngfftc(3),tdks%ph1d,tdks%xred) 116 if (psps%usepaw==1.and.tdks%pawfgr%usefinegrid==1) then 117 call getph(tdks%atindx,dtset%natom,tdks%pawfgr%ngfft(1),tdks%pawfgr%ngfft(2), & 118 & tdks%pawfgr%ngfft(3),tdks%ph1df,tdks%xred) 119 else 120 tdks%ph1df(:,:)=tdks%ph1d(:,:) 121 end if 122 123 !PAW specific 124 if (psps%usepaw==1) then 125 !Check for non-overlapping PAW spheres 126 call chkpawovlp(dtset%natom,psps%ntypat,dtset%pawovlp,tdks%pawtab,tdks%rmet, & 127 & dtset%typat,tdks%xred) 128 129 !Identify parts of the rectangular grid where the density has to be calculated 130 !FB: Needed? 131 optcut=0;optgr0=dtset%pawstgylm;optgr1=0;optgr2=0;optrad=1-dtset%pawstgylm 132 forces_needed=0 !FB TODO needs to be changed if Ehrenfest? 133 stress_needed=0 134 if ((forces_needed==1).or. & 135 & (dtset%xclevel==2.and.dtset%pawnhatxc>0.and.tdks%usexcnhat>0).or. & 136 & (dtset%positron/=0.and.forces_needed==2)) then 137 optgr1=dtset%pawstgylm; if (stress_needed==1) optrad=1; if (dtset%pawprtwf==1) optrad=1 138 end if 139 call nhatgrid(tdks%atindx1,tdks%gmet,my_natom,dtset%natom,tdks%nattyp, & 140 & tdks%pawfgr%ngfft,psps%ntypat,optcut,optgr0,optgr1, & 141 & optgr2,optrad,tdks%pawfgrtab,tdks%pawtab,tdks%rprimd, & 142 & dtset%typat, tdks%ucvol,tdks%xred,comm_atom=mpi_enreg%comm_atom, & 143 & mpi_atmtab=mpi_enreg%my_atmtab,comm_fft=mpi_enreg%comm_fft, & 144 & distribfft=mpi_enreg%distribfft) 145 endif 146 147 !!FB: @MT Needed? If yes, then don't forget to put it back in tdks_init/second_setup as well 148 !!if any nuclear dipoles are nonzero, compute the vector potential in real space (depends on 149 !!atomic position so should be done for nstep = 1 and for updated ion positions 150 !if ( any(abs(dtset%nucdipmom(:,:))>tol8) ) then 151 ! with_vectornd = 1 152 !else 153 ! with_vectornd = 0 154 !end if 155 !if(allocated(vectornd)) then 156 ! ABI_FREE(vectornd) 157 !end if 158 !ABI_MALLOC(vectornd,(with_vectornd*nfftf,3)) 159 !vectornd=zero 160 !if(with_vectornd .EQ. 1) then 161 ! call make_vectornd(1,gsqcut,psps%usepaw,mpi_enreg,dtset%natom,nfftf,ngfftf,dtset%nucdipmom,& 162 ! & rprimd,vectornd,xred) 163 164 end subroutine rttddft_setup_ele_step