Line data Source code
1 : !!****m* ABINIT/m_rttddft_tdef
2 : !! NAME
3 : !! m_rttddft_tdef
4 : !!
5 : !! FUNCTION
6 : !! Contains definition of the tdef type
7 : !! related to time-dependent electric field
8 : !!
9 : !! COPYRIGHT
10 : !! Copyright (C) 2021-2026 ABINIT group (FBrieuc)
11 : !! This file is distributed under the terms of the
12 : !! GNU General Public License, see ~abinit/COPYING
13 : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : !!
15 : !! SOURCE
16 :
17 : #if defined HAVE_CONFIG_H
18 : #include "config.h"
19 : #endif
20 :
21 : #include "abi_common.h"
22 :
23 : module m_rttddft_tdef
24 :
25 : use defs_basis
26 : use defs_abitypes, only: MPI_type
27 :
28 : use m_dtset, only: dataset_type
29 : use m_errors, only: msg_hndl
30 : use m_initylmg, only: initylmg
31 : use m_profiling_abi, only: abimem_record
32 : use m_xmpi, only: xmpi_bcast
33 :
34 : implicit none
35 :
36 : private
37 : !!***
38 :
39 : !! NAME
40 : !! tdef_type: Time Dependent Electric Field type
41 : !! Object containing the TD Electric field
42 : !! for RT-TDDFT runs
43 : !!
44 : !! SOURCE
45 : type,public :: tdef_type
46 :
47 : integer :: ef_type !type of TD electric feld (Dirac or sin^2 pulse)
48 : real(dp) :: efield(3) !TD external elec. field perturbation
49 : real(dp) :: ef_ezero(3) !E_0 = |E_0|*polarization
50 : real(dp) :: vecpot(3) !total vector potential
51 : real(dp) :: vecpot_ext(3) !external vector potential
52 : real(dp) :: vecpot_ind(3,2) !induced vector potential at t and t-dt
53 : real(dp) :: vecpot_red(3) !total vector potential in reduced coord. (in reciprocal space)
54 : real(dp) :: ef_tzero !time at which elec field is switched on
55 : real(dp) :: ef_omega !angular freq. of TD elec field
56 : real(dp) :: ef_tau !time width of the pulse
57 : real(dp) :: ef_sin_a !useful constant for sin^2 pulse
58 : real(dp) :: ef_sin_b !useful constant for sin^2 pulse
59 : real(dp), allocatable :: kpa(:,:) ! contains kpts + A
60 : ! (in reduced coordinates in reciprocal space)
61 : logical :: induced_vecpot ! Add the vector potential induced by the current density
62 : ! in the Hamiltonian
63 :
64 : contains
65 :
66 : procedure :: init => tdef_init
67 : procedure :: update => tdef_update
68 : procedure :: restart => tdef_restart
69 :
70 : end type tdef_type
71 : !!***
72 :
73 : contains
74 : !!***
75 :
76 : !!****f* m_rttddft/tdef_init
77 : !!
78 : !! NAME
79 : !! tdef_init
80 : !!
81 : !! FUNCTION
82 : !! Update value of electric field and vector potential at time t
83 : !!
84 : !! INPUTS
85 : !! [tdef = tdef structure to update]
86 : !! td_ef_type = type of electric field (Dirac or sin^2 pulse)
87 : !! td_ef_pol = polarization
88 : !! td_ef_ezero = Amplitude (E_0 = |E_0|*polarization)
89 : !! td_ef_tzero = time at which the pulse is switched on
90 : !! td_ef_lambda = wavelength (for sin^2 pulse)
91 : !! td_ef_tau = time width of the pulse (for finite-width delta kick or sin^2 pulse)
92 : !! time = propagation time
93 : !! nkpt = number of kpoints
94 : !! kpts = kpoints array
95 : !!
96 : !! OUTPUT
97 : !! [tdef = updated tdef structure]
98 : !!
99 : !! SOURCE
100 50 : subroutine tdef_init(tdef,td_ef_type,td_ef_pol,td_ef_ezero,td_ef_tzero,td_ef_lambda,td_ef_tau,td_ef_induced_vecpot,nkpt,kpts)
101 :
102 : !Arguments ------------------------------------
103 : !scalars
104 : class(tdef_type), intent(inout) :: tdef
105 : integer, intent(in) :: td_ef_type, td_ef_induced_vecpot, nkpt
106 : real(dp), intent(in) :: td_ef_ezero, td_ef_tzero, td_ef_lambda, td_ef_tau
107 : !arrays
108 : real(dp), intent(in) :: kpts(:,:)
109 :
110 : !Local variables-------------------------------
111 : real(dp), intent(in) :: td_ef_pol(3)
112 :
113 : ! ***********************************************************************
114 :
115 50 : if (td_ef_type > 1 .or. td_ef_type < 0) then
116 0 : ABI_ERROR("Wrong value of td_ef_type")
117 : end if
118 :
119 50 : tdef%ef_type = td_ef_type
120 200 : tdef%ef_ezero = td_ef_pol*td_ef_ezero
121 50 : tdef%ef_tau = td_ef_tau
122 50 : tdef%ef_omega = 2.0_dp*pi*Speed_Light/td_ef_lambda !2*pi*f=2*pi*c/lambda
123 50 : tdef%ef_tzero = td_ef_tzero
124 50 : tdef%ef_sin_a = 2.0_dp*pi/td_ef_tau + tdef%ef_omega
125 50 : tdef%ef_sin_b = 2.0_dp*pi/td_ef_tau - tdef%ef_omega
126 50 : if (td_ef_induced_vecpot == 0) then
127 48 : tdef%induced_vecpot = .false.
128 2 : else if (td_ef_induced_vecpot == 1) then
129 2 : tdef%induced_vecpot = .true.
130 : else
131 0 : ABI_ERROR("Wrong value of td_ef_induced_vecpot")
132 : end if
133 :
134 200 : tdef%efield = zero
135 200 : tdef%vecpot = zero
136 200 : tdef%vecpot_ext = zero
137 450 : tdef%vecpot_ind = zero
138 200 : tdef%vecpot_red = zero
139 :
140 150 : ABI_MALLOC(tdef%kpa,(3,nkpt))
141 3732 : tdef%kpa = kpts
142 :
143 50 : end subroutine tdef_init
144 : !!***
145 :
146 : !!****f* m_rttddft/tdef_update
147 : !!
148 : !! NAME
149 : !! tdef_update
150 : !!
151 : !! FUNCTION
152 : !! Update value of electric field and vector potential at time t
153 : !!
154 : !! INPUTS
155 : !! [tdef = tdef structure to update]
156 : !! dtset = dataset structure
157 : !! mpi_enreg = MPI communicators structure
158 : !! time = propagation time
159 : !! rprimd = cell vectors (direct space)
160 : !! gprimd = cell vectors (reciprocal space)
161 : !! kg = kpoints in reciprocal space
162 : !! mpsang = 1+maximum angular momentum for nonlocal pseudopotential (required by initylmg)
163 : !! npwarr = array holding npw for each k point
164 : !! ylm = real spherical harmonics for each G and k point
165 : !! ylmgr = gradient of real spherical harmonics for each G and k point
166 : !! current = total current density
167 : !!
168 : !! OUTPUT
169 : !! [tdef = updated tdef structure]
170 : !!
171 : !! SOURCE
172 1795 : subroutine tdef_update(tdef,dtset,mpi_enreg,time,rprimd,gprimd,kg,mpsang,npwarr,ylm,ylmgr,current,update_vecpot_ind)
173 :
174 : !Arguments ------------------------------------
175 : !scalars
176 : class(tdef_type), intent(inout) :: tdef
177 : type(dataset_type), intent(inout) :: dtset
178 : type(MPI_type), intent(inout) :: mpi_enreg
179 : real(dp), intent(in) :: time
180 : real(dp), intent(in) :: rprimd(3,3)
181 : real(dp), intent(in) :: gprimd(3,3)
182 : integer, intent(in) :: kg(:,:)
183 : integer, intent(in) :: mpsang
184 : integer, intent(in) :: npwarr(:)
185 : real(dp), intent(in) :: current(:,:)
186 : real(dp), intent(out) :: ylm(:,:)
187 : real(dp), intent(out) :: ylmgr(:,:,:)
188 : logical, optional, intent(in) :: update_vecpot_ind
189 :
190 : !Local variables-------------------------------
191 : character(len=500) :: msg
192 : integer :: i
193 : logical :: lvecpot_ind
194 : real(dp) :: tmp(3), expt
195 :
196 : ! ***********************************************************************
197 :
198 : ! Update potential vector if time-dependent electric field perturbation is present
199 2916 : select case (tdef%ef_type)
200 : !No external field
201 : case (0)
202 : !Dirac pulse: vector potential is just an Heaviside function
203 : case (1)
204 1121 : if (abs(time-tdef%ef_tzero)<tol16) then
205 44 : tdef%efield(:) = tdef%ef_ezero(:)
206 : else
207 4440 : tdef%efield(:) = zero
208 : end if
209 1121 : if (time >= tdef%ef_tzero) then
210 4324 : tdef%vecpot_ext(:) = -tdef%ef_ezero(:)
211 : end if
212 : !"Finite" delta-kick pulse: Vector potential is a finite width sigmoid function
213 : case (2)
214 0 : expt = exp(-(time-tdef%ef_tzero)/tdef%ef_tau)
215 0 : tdef%efield(:) = tdef%ef_ezero(:)/tdef%ef_tau * expt/(1+expt)**2
216 0 : tdef%vecpot_ext(:) = -tdef%ef_ezero(:)/(1+expt)**2
217 : !Pulse with sin^2 shape:
218 : !E(t) = E0*cos(w*(t-t0))*sin^2(pi*(t-t0)/tau)
219 : !A(t) = -(E0/2w)*sin(w*(t-t0))+E0/(4*(2pi/tau+w))*sin((2pi/tau+w)*(t-t0))+E0/(4(2pi/tau-w))*sin((2pi/tau-w)*(t-t0))
220 : ! case(2)
221 : ! if (time >= tdef%ef_tzero+tdef%ef_tau) then
222 : ! tdef%efield(:) = zero
223 : ! else if (time >= tdef%ef_tzero) then
224 : ! t = time-tdef%ef_tzero
225 : ! tdef%efield(:) = tdef%ef_ezero*cos(tdef%ef_omega*t)*sin(pi*t/tdef%ef_tau)**2
226 : ! tdef%vecpot_ext(:) = tdef%ef_ezero*(-sin(tdef%ef_omega*t)/(two*tdef%ef_omega) &
227 : ! & +sin(tdef%ef_sin_a*t)/(four*tdef%ef_sin_a) &
228 : ! & +sin(tdef%ef_sin_b*t)/(four*tdef%ef_sin_b))
229 : ! end if
230 : case default
231 0 : write(msg,"(a)") "Unknown electric field type - check the value of td_ef_type"
232 1795 : ABI_ERROR(msg)
233 : end select
234 :
235 1795 : lvecpot_ind = .true.
236 1795 : if (present(update_vecpot_ind)) lvecpot_ind = update_vecpot_ind
237 :
238 : !Induced vector potential
239 : !Should deal with sppol?! How?
240 : !d^2A_ind/dt^2 = 4piJ(t)
241 : !A_ind(t+dt) = 2*A_ind(t) - A_ind(t-dt) + 4*pi*dt**2*J(t)
242 1795 : if (tdef%induced_vecpot) then
243 102 : if (lvecpot_ind) then
244 404 : tmp = tdef%vecpot_ind(:,2) ! t - 2dt
245 404 : tdef%vecpot_ind(:,2) = tdef%vecpot_ind(:,1) ! t - dt
246 404 : tdef%vecpot_ind(:,1) = 2*tdef%vecpot_ind(:,2) - tmp(:) + four*pi*(dtset%dtele**2)*current(:,1) ! t
247 : end if
248 408 : tdef%vecpot = tdef%vecpot_ext + tdef%vecpot_ind(:,1)
249 : else
250 6772 : tdef%vecpot = tdef%vecpot_ext
251 : end if
252 :
253 7180 : tdef%vecpot_red = matmul(transpose(rprimd),tdef%vecpot)
254 :
255 1795 : if (tdef%ef_type /= 0) then
256 : !Update the k+A grid used in cprojs
257 : !Divide by 2pi here seems necessary
258 4484 : do i = 1,3
259 112100 : tdef%kpa(i,:) = dtset%kptns(i,:) + tdef%vecpot_red(i)/(two*pi)
260 : end do
261 : ! update the spherical harmonics (computed at k+G+A)
262 : call initylmg(gprimd,kg,tdef%kpa,dtset%mkmem,mpi_enreg,mpsang,dtset%mpw,dtset%nband, &
263 1121 : & dtset%nkpt,npwarr,dtset%nsppol,0,rprimd,ylm,ylmgr)
264 : end if
265 :
266 1795 : end subroutine tdef_update
267 : !!***
268 :
269 : !!****f* m_rttddft/tdef_restart
270 : !!
271 : !! NAME
272 : !! tdef_restart
273 : !!
274 : !! FUNCTION
275 : !! Update some values for restart of calculation with TD electric field
276 : !! Essentially needed if we account for induced vector potential
277 : !!
278 : !! INPUTS
279 : !! [tdef = tdef structure to update]
280 : !! mpi_enreg = MPI communicators structure
281 : !! restart_unit = unit of restart file to read
282 : !!
283 : !! OUTPUT
284 : !! [tdef = updated tdef structure]
285 : !!
286 : !! SOURCE
287 20 : subroutine tdef_restart(tdef,mpi_enreg,restart_unit)
288 :
289 : !Arguments ------------------------------------
290 : !scalars
291 : class(tdef_type), intent(inout) :: tdef
292 : type(MPI_type), intent(inout) :: mpi_enreg
293 : integer, intent(in) :: restart_unit
294 :
295 : !Local variables-------------------------------
296 : integer :: ierr
297 :
298 : ! ***********************************************************************
299 :
300 20 : if (tdef%induced_vecpot) then
301 1 : if (mpi_enreg%me == 0) then
302 1 : read(restart_unit,*) tdef%vecpot_ind(:,2)
303 1 : read(restart_unit,*) tdef%vecpot_ind(:,1)
304 : end if
305 : end if
306 : !Send to all procs
307 20 : call xmpi_bcast(tdef%vecpot_ind,0,mpi_enreg%comm_world,ierr)
308 :
309 20 : end subroutine tdef_restart
310 : !!***
311 :
312 0 : end module m_rttddft_tdef
313 : !!***
|
