Branch data 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 [ - + - + : 1755 : 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 : 1755 : real(dp),allocatable :: grchempottn(:,:)
202 : 1755 : real(dp),allocatable :: grewtn(:,:)
203 : : real(dp),parameter :: k0(3)=(/zero,zero,zero/)
204 [ - + ]: 1755 : real(dp),allocatable :: kxc(:,:)
205 : : real(dp) :: strsxc(6)
206 : : real(dp) :: vpotzero(2)
207 : :
208 : : ! ***********************************************************************
209 : :
210 : 1755 : my_natom=mpi_enreg%my_natom
211 : :
212 : 1755 : 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 : 1755 : optene = 4; nkxc=0; moved_atm_inside=0; moved_rhor=1
225 [ + - ]: 1755 : n1xccc=0;if (psps%n1xccc/=0) n1xccc=psps%n1xccc
226 [ + - ]: 1755 : n3xccc=0;if (psps%n1xccc/=0) n3xccc=tdks%pawfgr%nfft
227 : 1755 : strsxc(:)=zero
228 : : !FB: tfw_activated is a save variable in scfcv, should maybe check where it appears again
229 : 1755 : tfw_activated=.false.
230 [ - + ]: 1755 : if (dtset%tfkinfunc==12) tfw_activated=.true.
231 [ - + + - : 5265 : ABI_MALLOC(grchempottn,(3,dtset%natom))
+ - + - -
+ - + ]
232 [ - + + - : 3510 : ABI_MALLOC(grewtn,(3,dtset%natom))
- + - + ]
233 [ - + + - : 3510 : ABI_MALLOC(kxc,(tdks%pawfgr%nfft,nkxc))
- + - + ]
234 : 1755 : calc_ewald = .false.
235 [ + - + + ]: 1755 : 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 [ + - - - : 1755 : & add_tfw=tfw_activated,xcctau3d=tdks%xcctau3d,calc_ewald=calc_ewald)
- - + - -
- - - - +
- + ]
247 [ - + ]: 1755 : ABI_FREE(grchempottn)
248 [ - + ]: 1755 : ABI_FREE(grewtn)
249 [ - + ]: 1755 : ABI_FREE(kxc)
250 : :
251 : : ! Get electronic temperature from dtset
252 [ - + - - : 1755 : 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 [ - + - - ]: 1755 : 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 [ + + ]: 1755 : if (psps%usepaw==1)then
261 : : !** Local exact exch.: impose occ. matrix if required
262 [ - + ]: 1545 : 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 : 1545 : call paw_an_reset_flags(tdks%paw_an)
277 : : !FB: @MT Changed self_consistent to false here. Is this right?
278 : 1545 : call paw_ij_reset_flags(tdks%paw_ij,self_consistent=.false.)
279 : 1545 : 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 [ + - - - : 1545 : & 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 [ - + ]: 1545 : 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 [ + + + + : 24207420 : tdks%vtrial(:,:)=tdks%vtrial(:,:)+SUM(vpotzero(:))
+ + ]
301 : : if(option/=1)then
302 : : !Fix the direct total energy (non-zero only for charged systems)
303 [ + + ]: 4635 : 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 [ + + ]: 4635 : & vpotzero(2)*dtset%cellcharge(1)
308 : : end if
309 : :
310 : : !** Dij computation
311 : : !FB: @MT fatvshift?
312 : 1545 : 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 [ + - - - : 1545 : & 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 [ + - - - : 1545 : & 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 [ + + + - ]: 1755 : 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 [ + - - - : 1755 : & gpu_option=dtset%gpu_option)
- - - + ]
340 : :
341 : 1755 : end subroutine rttddft_init_hamiltonian
342 : : !!***
343 : :
344 : : end module m_rttddft
345 : : !!***
|
