Branch data 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 : : !!***
|
