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