Branch data Line data Source code
1 : : !!****m* ABINIT/m_rttddft_propagators
2 : : !! NAME
3 : : !! m_rttddft_propagators
4 : : !!
5 : : !! FUNCTION
6 : : !! Contains various propagators for the KS orbitals
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_propagators
23 : :
24 : : use defs_basis
25 : : use defs_abitypes, only: MPI_type
26 : : use defs_datatypes, only: pseudopotential_type
27 : :
28 : : use m_bandfft_kpt, only: bandfft_kpt, bandfft_kpt_type, &
29 : : & bandfft_kpt_set_ikpt, &
30 : : & prep_bandfft_tabs
31 : : use m_dtset, only: dataset_type
32 : : use m_energies, only: energies_type
33 : : use m_errors, only: msg_hndl
34 : : use m_gemm_nonlop_projectors, only: set_gemm_nonlop_ikpt
35 : : use m_hamiltonian, only: gs_hamiltonian_type, gspot_transgrid_and_pack
36 : : use m_invovl, only: make_invovl
37 : : use m_kg, only: mkkin, mkkpg
38 : : use m_mkffnl, only: mkffnl
39 : : use m_mpinfo, only: proc_distrb_cycle
40 : : use m_profiling_abi, only: abimem_record
41 : : use m_rttddft, only: rttddft_init_hamiltonian
42 : : use m_rttddft_exponential, only: rttddft_exp_taylor
43 : : use m_rttddft_properties, only: rttddft_calc_density, &
44 : : & rttddft_calc_occ, &
45 : : & rttddft_calc_kin
46 : : use m_rttddft_tdks, only: tdks_type
47 : : use m_specialmsg, only: wrtout
48 : : use m_xmpi, only: xmpi_comm_rank, xmpi_sum, xmpi_max
49 : :
50 : : implicit none
51 : :
52 : : private
53 : : !!***
54 : :
55 : : public :: rttddft_propagator_er
56 : : public :: rttddft_propagator_emr
57 : : !!***
58 : :
59 : : contains
60 : : !!***
61 : :
62 : : !!****f* m_rttddft/rttddft_propagator_er
63 : : !!
64 : : !! NAME
65 : : !! rttddft_propagator_er
66 : : !!
67 : : !! FUNCTION
68 : : !! Main subroutine to propagate the KS orbitals using
69 : : !! the Exponential Rule (ER) propagator.
70 : : !!
71 : : !! INPUTS
72 : : !! dtset <type(dataset_type)> = all input variables for this dataset
73 : : !! ham_k <type(gs_hamiltonian_type)> = Hamiltonian object
74 : : !! istep <integer> = step number
75 : : !! mpi_enreg <MPI_type> = MPI-parallelisation information
76 : : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
77 : : !! tdks <type(tdks_type)> = the tdks object to initialize
78 : : !! calc_properties <logical> = logical governing the computation of
79 : : !! some properties (energies, occupations, eigenvalues..)
80 : : !!
81 : : !! NOTES
82 : : !! Other propagators such as the Exponential Midpoint Rule (EMR)
83 : : !! should usually be prefered over this one since the ER propagator
84 : : !! alone violates time reversal symmetry. Using this propagator with
85 : : !! the exponential approximated by Taylor expansion of order 1 leads
86 : : !! to the famous Euler method which is fast and simple but unstable
87 : : !! and thus insufficient for RT-TDDFT.
88 : : !!
89 : : !! SOURCE
90 : 1755 : subroutine rttddft_propagator_er(dtset, istep, mpi_enreg, psps, tdks, calc_properties)
91 : :
92 : : !Arguments ------------------------------------
93 : : !scalars
94 : : integer, intent(in) :: istep
95 : : logical, optional, intent(in) :: calc_properties
96 : : type(dataset_type), intent(inout) :: dtset
97 : : type(MPI_type), intent(inout) :: mpi_enreg
98 : : type(pseudopotential_type), intent(inout) :: psps
99 : : type(tdks_type), target, intent(inout) :: tdks
100 : :
101 : : !Local variables-------------------------------
102 : : !scalars
103 : : integer :: bdtot_index
104 : : integer :: calc_forces
105 : : integer :: dimffnl
106 : : integer :: iband
107 : : integer :: ibg, icg
108 : : integer :: ider, idir
109 : : integer :: ierr, ilm
110 : : integer :: ikpt, ikpt_loc, ikg
111 : : integer :: isppol
112 : : integer :: istwf_k
113 : : integer :: me_distrb
114 : : integer :: me_bandfft
115 : : integer :: my_ikpt, my_nspinor
116 : : integer :: nband_k, nband_k_mem
117 : : integer :: npw_k, nkpg
118 : : integer :: shift
119 : : integer :: spaceComm_distrb
120 : : integer :: n4, n5, n6
121 : : logical :: with_vxctau
122 : : logical :: lcalc_properties
123 : : type(energies_type) :: energies
124 : : type(bandfft_kpt_type),pointer :: my_bandfft_kpt => null()
125 : 1755 : type(gs_hamiltonian_type) :: ham_k
126 : : !arrays
127 : 1755 : integer, allocatable :: kg_k(:,:)
128 : : real(dp), pointer :: cg(:,:) => null()
129 : : real(dp), pointer :: cg0(:,:) => null()
130 : 1755 : real(dp), allocatable :: enl(:)
131 : : real(dp), pointer :: eig(:) => null()
132 : 1755 : real(dp), allocatable :: ffnl(:,:,:,:)
133 : 1755 : real(dp), allocatable :: kpg_k(:,:)
134 : : real(dp) :: kpoint(3)
135 : : real(dp) :: kpa(3)
136 : 1755 : real(dp), allocatable :: kinpw(:)
137 : : real(dp), pointer :: occ(:) => null()
138 : : real(dp), pointer :: occ0(:) => null()
139 : 1755 : real(dp), allocatable :: ph3d(:,:,:)
140 : 1755 : real(dp), allocatable :: vlocal(:,:,:,:)
141 : 1755 : real(dp), allocatable :: vxctaulocal(:,:,:,:,:)
142 : 1755 : real(dp), allocatable :: ylm_k(:,:)
143 : : logical :: lproperties(4)
144 : :
145 : : ! ***********************************************************************
146 : :
147 : : !Init MPI
148 : 1755 : spaceComm_distrb=mpi_enreg%comm_cell
149 : 960 : if (mpi_enreg%paral_kgb==1) spaceComm_distrb=mpi_enreg%comm_kpt
150 : 1755 : me_distrb=xmpi_comm_rank(spaceComm_distrb)
151 : :
152 : : !Do we calculate properties?
153 : : !Governed by lproperties:
154 : : ! lproperties(1) = compute energy contributions (kinetic)
155 : : ! lproperties(2) = NL energy contribution in NC case
156 : : ! lproperties(3) = eigenvalues
157 : : ! lproperties(4) = occupations
158 : 1755 : lproperties(:) = .false.
159 : 1755 : lcalc_properties = .false.
160 [ + + ]: 1755 : if (present(calc_properties)) then
161 : 1725 : lcalc_properties = calc_properties
162 [ + - ]: 1725 : if (lcalc_properties) then
163 : : !compute energy contributions
164 : 1725 : lproperties(1) = .true.
165 : : !Init to zero different energies
166 : 1725 : call energies%init()
167 : 1725 : energies%entropy_ks=tdks%energies%entropy_ks
168 : 1725 : energies%e_corepsp=tdks%energies%e_corepsp
169 : 1725 : energies%e_ewald=tdks%energies%e_ewald
170 : : !including NL part in NC case?
171 [ + + ]: 1725 : if (dtset%usepaw == 0) then
172 : 200 : lproperties(2) = .true.
173 [ - + + - : 600 : ABI_MALLOC(enl,(mpi_enreg%bandpp))
+ - + - -
+ - + ]
174 : : else
175 [ - + ]: 1525 : ABI_MALLOC(enl,(0))
176 : : end if
177 : : !other properties
178 : : ! eigenvalues
179 [ + - - + ]: 1725 : if (dtset%prteig /= 0 .or. dtset%prtdos /= 0) then
180 [ # # ]: 0 : if (mod(istep-1,dtset%td_prtstr) == 0) then
181 : 0 : lproperties(3) = .true.
182 [ # # ]: 0 : tdks%eigen(:) = zero
183 : : end if
184 : : end if
185 : : ! occupations
186 : 1725 : lproperties(4) = .true.
187 [ + + ]: 395485 : tdks%occ(:) = zero
188 : : end if
189 : : end if
190 : :
191 : : !Set "vtrial" and initialize the Hamiltonian
192 : 1755 : call rttddft_init_hamiltonian(dtset,energies,ham_k,istep,mpi_enreg,psps,tdks)
193 : :
194 : 1755 : my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
195 : 1755 : n4=dtset%ngfft(4); n5=dtset%ngfft(5); n6=dtset%ngfft(6)
196 [ - + + - : 10530 : ABI_MALLOC(vlocal,(n4,n5,n6,ham_k%nvloc))
- + + - -
+ + - - +
+ - + - +
- - + - +
- + ]
197 [ + + ]: 7020 : with_vxctau=(size(tdks%vxctau)>0)
198 [ - + ]: 1755 : if(with_vxctau) then
199 [ # # # # : 0 : ABI_MALLOC(vxctaulocal,(n4,n5,n6,ham_k%nvloc,4))
# # # # #
# # # # #
# # # # #
# # # # #
# # # # ]
200 : : end if
201 [ + - ]: 1755 : if (dtset%ionmov/=0) then
202 : : calc_forces=1
203 : : else
204 : 1755 : calc_forces=0
205 : : end if
206 : :
207 : : !FB: @MT Needed?
208 : : !has_vectornd = (with_vectornd .EQ. 1)
209 : : !if(has_vectornd) then
210 : : ! ABI_MALLOC(vectornd_pac,(n4,n5,n6,ham_k%nvloc,3))
211 : : ! vectornd_pac=zero
212 : : !end if
213 : :
214 : 1755 : icg=0; ibg=0
215 : 1755 : bdtot_index=0
216 : :
217 : : !*** LOOP OVER SPINS
218 [ + + ]: 3510 : do isppol=1,dtset%nsppol
219 : :
220 : 1755 : ikpt_loc=0
221 : 1755 : ikg=0
222 : :
223 : : ! Set up local potential vlocal on the coarse FFT mesh from vtrial taking into account the spin.
224 : : ! Also, continue to initialize the Hamiltonian.
225 : : call gspot_transgrid_and_pack(isppol, psps%usepaw, dtset%paral_kgb, dtset%nfft, dtset%ngfft, tdks%nfftf, &
226 : 1755 : & dtset%nspden, ham_k%nvloc, 1, tdks%pawfgr, mpi_enreg, tdks%vtrial, vlocal)
227 : 1755 : call ham_k%load_spin(isppol, vlocal=vlocal, with_nonlocal=.true.)
228 : :
229 [ - + ]: 1755 : if (with_vxctau) then
230 : : call gspot_transgrid_and_pack(isppol, psps%usepaw, dtset%paral_kgb, dtset%nfft, dtset%ngfft, tdks%nfftf, &
231 : 0 : & dtset%nspden, ham_k%nvloc, 4, tdks%pawfgr, mpi_enreg, tdks%vxctau, vxctaulocal)
232 : 0 : call ham_k%load_spin(isppol, vxctaulocal=vxctaulocal)
233 : : end if
234 : :
235 : : !FB: @MT Needed?
236 : : ! ! if vectornd is present, set it up for addition to ham_k similarly to how it's done for
237 : : ! ! vtrial. Note that it must be done for the three Cartesian directions. Also, the following
238 : : ! ! code assumes explicitly and implicitly that nvloc = 1. This should eventually be generalized.
239 : : ! if(has_vectornd) then
240 : : ! do idir = 1, 3
241 : : ! ABI_MALLOC(cgrvtrial,(dtset%nfft,dtset%nspden))
242 : : ! call transgrid(1,mpi_enreg,dtset%nspden,-1,0,0,dtset%paral_kgb,pawfgr,rhodum,rhodum,cgrvtrial,vectornd(:,idir))
243 : : ! call fftpac(isppol,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,cgrvtrial,vectornd_pac(:,:,:,1,idir),2)
244 : : ! ABI_FREE(cgrvtrial)
245 : : ! end do
246 : : ! call ham_k%load_spin(isppol, vectornd=vectornd_pac)
247 : : ! end if
248 : :
249 : : !*** BIG FAT k POINT LOOP
250 : 1755 : ikpt = 0
251 [ + + ]: 50570 : do while (ikpt_loc < dtset%nkpt)
252 : :
253 : 47060 : ikpt_loc = ikpt_loc + 1
254 : 47060 : ikpt = ikpt_loc
255 : 47060 : my_ikpt = mpi_enreg%my_kpttab(ikpt)
256 : :
257 : 47060 : nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
258 [ + + ]: 47060 : if (mpi_enreg%paral_kgb==1) then
259 : 26240 : nband_k_mem=mpi_enreg%bandpp
260 : : else
261 : : nband_k_mem=nband_k
262 : : end if
263 : 47060 : istwf_k=dtset%istwfk(ikpt)
264 : 47060 : npw_k=tdks%npwarr(ikpt)
265 : :
266 [ + + ]: 47060 : if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) then
267 : 19680 : bdtot_index=bdtot_index+nband_k
268 : 19680 : cycle
269 : : end if
270 : :
271 [ + + ]: 27380 : if (mpi_enreg%paral_kgb==1) my_bandfft_kpt => bandfft_kpt(my_ikpt)
272 : 27380 : call bandfft_kpt_set_ikpt(ikpt,mpi_enreg)
273 : :
274 [ - + + - : 82140 : ABI_MALLOC(kg_k,(3,npw_k))
+ - + - -
+ - + -
+ ]
275 [ - + + - : 109520 : ABI_MALLOC(ylm_k,(npw_k,psps%mpsang*psps%mpsang*psps%useylm))
+ + + - +
- + + - +
- + - + ]
276 [ + + + + ]: 22474460 : kg_k(:,1:npw_k)=tdks%kg(:,1+ikg:npw_k+ikg)
277 [ + + ]: 27380 : if (psps%useylm==1) then
278 [ + + ]: 136300 : do ilm=1,psps%mpsang*psps%mpsang
279 [ + + ]: 21448840 : ylm_k(1:npw_k,ilm)=tdks%ylm(1+ikg:npw_k+ikg,ilm)
280 : : end do
281 : : end if
282 : :
283 : : !** Set up the remaining k-dependent part of the Hamiltonian
284 [ + + ]: 109520 : kpoint(:)=dtset%kptns(:,ikpt) !k-point
285 [ + + ]: 109520 : kpa(:)=tdks%tdef%kpa(:,ikpt) !k+A
286 : :
287 : : ! Kinetic energy
288 [ - + + - : 82140 : ABI_MALLOC(kinpw,(npw_k))
+ - + - -
+ - + -
+ ]
289 : 27380 : call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,tdks%gmet,kg_k,kinpw,kpoint,npw_k,0,0,tdks%tdef%vecpot_red)
290 : :
291 : : ! Compute (k+G) vectors (only if useylm=1)
292 : 27380 : nkpg=3*calc_forces*dtset%nloalg(3)
293 [ - + + - : 109520 : ABI_MALLOC(kpg_k,(npw_k,nkpg))
+ - - - +
- - + - +
- + - + ]
294 : : ! Compute nonlocal form factors ffnl at all (k+G):
295 : 27380 : ider=0;idir=0;dimffnl=1
296 [ - + + - : 136900 : ABI_MALLOC(ffnl,(npw_k,dimffnl,psps%lmnmax,psps%ntypat))
- + + - -
+ + - + -
+ - - + -
+ - + ]
297 [ + + + + : 27380 : if (mpi_enreg%paral_kgb/=1 .or. istep<=tdks%first_step .or. tdks%tdef%ef_type/=0) then
+ + ]
298 [ - + ]: 27236 : if (nkpg > 0) call mkkpg(kg_k,kpg_k,kpa,nkpg,npw_k)
299 : : call mkffnl(psps%dimekb,dimffnl,psps%ekb,ffnl,psps%ffspl,tdks%gmet,tdks%gprimd, &
300 : : & ider,idir,psps%indlmn,kg_k,kpg_k,kpa,psps%lmnmax,psps%lnmax, &
301 : : & psps%mpsang,psps%mqgrid_ff,nkpg,npw_k,psps%ntypat,psps%pspso, &
302 : 27236 : & psps%qgrid_ff,tdks%rmet,psps%usepaw,psps%useylm,ylm_k,tdks%ylmgr)
303 : : end if
304 : :
305 : : !** Load k-dependent part in the Hamiltonian datastructure
306 : : !** - Compute 3D phase factors
307 : : !** - Prepare various tabs in case of band-FFT parallelism
308 : : !** - Load k-dependent quantities in the Hamiltonian
309 [ - + + - : 109520 : ABI_MALLOC(ph3d,(2,npw_k,ham_k%matblk))
- + + - +
- + - - +
- + - + ]
310 : : call ham_k%load_k(kpt_k=kpa,istwf_k=istwf_k,npw_k=npw_k,kinpw_k=kinpw,kg_k=kg_k,kpg_k=kpg_k, &
311 : 27380 : & ffnl_k=ffnl,ph3d_k=ph3d,compute_ph3d=.true.,compute_gbound=.true.)
312 : :
313 : : !** Load band-FFT tabs (transposed k-dependent arrays)
314 [ + + ]: 27380 : if (mpi_enreg%paral_kgb==1) then
315 [ + + + + ]: 6560 : if (istep<=tdks%first_step .or. tdks%tdef%ef_type/=0) then
316 : 6416 : call prep_bandfft_tabs(ham_k,ikpt,dtset%mkmem,mpi_enreg)
317 : : end if
318 : : call ham_k%load_k(npw_fft_k=my_bandfft_kpt%ndatarecv, &
319 : : & gbound_k =my_bandfft_kpt%gbound, &
320 : : & kinpw_k =my_bandfft_kpt%kinpw_gather, &
321 : : & kg_k =my_bandfft_kpt%kg_k_gather, &
322 : : & kpg_k =my_bandfft_kpt%kpg_k_gather, &
323 : : & ffnl_k =my_bandfft_kpt%ffnl_gather, &
324 : 6560 : & ph3d_k =my_bandfft_kpt%ph3d_gather)
325 : : end if
326 : :
327 : : !** Build inverse of overlap matrix
328 [ + + + + : 27380 : if(psps%usepaw == 1 .and. (istep <= tdks%first_step .or. tdks%tdef%ef_type/=0)) then
+ + ]
329 : 16164 : call make_invovl(ham_k,dimffnl,ffnl,ph3d,mpi_enreg)
330 : : end if
331 : :
332 : : ! Setup gemm_nonlop
333 [ - + ]: 27380 : if (tdks%gemm_nonlop_use_gemm) then
334 : : call set_gemm_nonlop_ikpt(my_ikpt,ham_k%npw_fft_k,ham_k%istwf_k,ham_k%indlmn,&
335 : 0 : & ham_k%ntypat,ham_k%nattyp,ham_k%gpu_option)
336 : : end if
337 : :
338 : : !** Compute the exp[-i(S^{-1})H]*cg using Taylor expansion to approximate the exponential
339 : 27380 : cg => tdks%cg(:,icg+1:icg+nband_k*npw_k*my_nspinor)
340 : : ! Compute properties "on-the-fly" if required
341 [ + + ]: 27380 : if (lcalc_properties) then
342 : 27260 : cg0 => tdks%cg0(:,icg+1:icg+nband_k*npw_k*my_nspinor)
343 : 27260 : occ => tdks%occ(bdtot_index+1:bdtot_index+nband_k)
344 : 27260 : occ0 => tdks%occ0(bdtot_index+1:bdtot_index+nband_k)
345 [ + + ]: 27260 : if (dtset%paral_kgb /= 1) then
346 : : shift = bdtot_index
347 : : else
348 : 6560 : me_bandfft = xmpi_comm_rank(mpi_enreg%comm_band)
349 : 6560 : shift = bdtot_index+me_bandfft*mpi_enreg%bandpp
350 : : end if
351 : : ! kinetic energy
352 [ + - ]: 27260 : if (lproperties(1)) then
353 : : call rttddft_calc_kin(energies%e_kinetic,cg,dtset,ham_k,nband_k,npw_k,my_nspinor, &
354 [ - + - + ]: 27260 : & occ0,dtset%wtk(ikpt),mpi_enreg,my_bandfft_kpt)
355 : : end if
356 : : ! for eigenvalues
357 [ - + ]: 27260 : if (lproperties(3)) then
358 : 0 : eig => tdks%eigen(1+shift:nband_k_mem+shift)
359 : : end if
360 : : ! occupations
361 [ + - ]: 27260 : if (lproperties(4)) then
362 : : !computed at istep-1 like energies and eigenvalues
363 : : call rttddft_calc_occ(cg,cg0,dtset,ham_k,ikpt,ibg,isppol,mpi_enreg, &
364 [ - + - + : 27260 : & nband_k,npw_k,my_nspinor,occ,occ0,tdks)
- + - + ]
365 : : end if
366 : :
367 : : !Propagate cg and compute the requested properties
368 [ - + ]: 27260 : call rttddft_exp_taylor(cg,dtset,ham_k,mpi_enreg,nband_k,npw_k,my_nspinor,enl=enl,eig=eig)
369 : :
370 : : !Finish computing NL PSP part for NC PSP
371 [ + + ]: 27260 : if (lproperties(2)) then
372 [ + + ]: 6080 : do iband = 1, mpi_enreg%bandpp
373 : 6080 : energies%e_nlpsp_vfock=energies%e_nlpsp_vfock+dtset%wtk(ikpt)*tdks%occ0(shift+iband)*enl(iband)
374 : : end do
375 : : end if
376 : : else
377 : : !Propagate cg only
378 [ - + ]: 120 : call rttddft_exp_taylor(cg,dtset,ham_k,mpi_enreg,nband_k,npw_k,my_nspinor)
379 : : end if
380 : :
381 [ - + ]: 27380 : ABI_FREE(kg_k)
382 [ - + ]: 27380 : ABI_FREE(ylm_k)
383 [ - + ]: 27380 : ABI_FREE(kpg_k)
384 [ - + ]: 27380 : ABI_FREE(kinpw)
385 [ - + ]: 27380 : ABI_FREE(ffnl)
386 [ - + ]: 27380 : ABI_FREE(ph3d)
387 : :
388 : 27380 : bdtot_index = bdtot_index+nband_k
389 : :
390 : : !** Also shift array memory if dtset%mkmem/=0
391 [ - + ]: 29135 : if (dtset%mkmem/=0) then
392 : 27380 : ibg=ibg+nband_k_mem
393 : 27380 : icg=icg+npw_k*my_nspinor*nband_k
394 : 27380 : ikg=ikg+npw_k
395 : : end if
396 : :
397 : : end do !nkpt
398 : :
399 : : end do !nsppol
400 : :
401 : : ! Free memory
402 [ - + ]: 1755 : ABI_FREE(vlocal)
403 [ - + ]: 1755 : if(dtset%usekden/=0) then
404 [ # # ]: 0 : ABI_FREE(vxctaulocal)
405 : : end if
406 : 1755 : call ham_k%free()
407 [ + + ]: 1755 : if (allocated(enl)) then
408 : 1725 : ABI_FREE(enl)
409 : : end if
410 : :
411 : : !Keep the computed energies in memory
412 [ + + ]: 1755 : if (lcalc_properties) then
413 : 1725 : call xmpi_sum(energies%e_kinetic,mpi_enreg%comm_kptband,ierr)
414 : 1725 : call xmpi_sum(energies%e_nlpsp_vfock,mpi_enreg%comm_kptband,ierr)
415 : 1725 : call energies%copy(tdks%energies)
416 [ - + ]: 1725 : if (lproperties(3)) call xmpi_sum(tdks%eigen,mpi_enreg%comm_kptband,ierr)
417 [ + - ]: 1725 : if (lproperties(4)) call xmpi_sum(tdks%occ,mpi_enreg%comm_kpt,ierr)
418 : : end if
419 : :
420 [ + + - + : 3510 : end subroutine rttddft_propagator_er
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + ]
421 : : !!***
422 : :
423 : : !!****f* m_rttddft/rttddft_propagator_emr
424 : : !!
425 : : !! NAME
426 : : !! rttddft_propagator_emr
427 : : !!
428 : : !! FUNCTION
429 : : !! Main subroutine to propagate the KS orbitals using
430 : : !! the Exponential Midpoint Rule (EMR) propagator
431 : : !!
432 : : !! INPUTS
433 : : !! dtset <type(dataset_type)>=all input variables for this dataset
434 : : !! ham_k <type(gs_hamiltonian_type)> = Hamiltonian object
435 : : !! istep <integer> = step number
436 : : !! mpi_enreg <MPI_type> = MPI-parallelisation information
437 : : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
438 : : !! tdks <type(tdks_type)> = the tdks object to initialize
439 : : !!
440 : : !! OUTPUT
441 : : !!
442 : : !! NOTES
443 : : !! This propagator is time reversible
444 : : !! (if H(t+dt/2) and the exponential are computed exactly).
445 : : !!
446 : : !! SOURCE
447 : 20 : subroutine rttddft_propagator_emr(dtset, istep, mpi_enreg, psps, tdks)
448 : :
449 : : !Arguments ------------------------------------
450 : : !scalars
451 : : integer, intent(in) :: istep
452 : : type(dataset_type), intent(inout) :: dtset
453 : : type(MPI_type), intent(inout) :: mpi_enreg
454 : : type(pseudopotential_type), intent(inout) :: psps
455 : : type(tdks_type), intent(inout) :: tdks
456 : :
457 : : !Local variables-------------------------------
458 : : integer :: me
459 : : !scalars
460 : : character(len=500) :: msg
461 : : integer :: ics
462 : : integer :: ierr
463 : : logical :: lconv
464 : : !arrays
465 : 40 : real(dp) :: cg(SIZE(tdks%cg(:,1)),SIZE(tdks%cg(1,:)))
466 : 40 : real(dp) :: diff(SIZE(tdks%cg(:,1)),SIZE(tdks%cg(1,:)))
467 : : real(dp) :: max_diff(2)
468 : :
469 : : ! ***********************************************************************
470 : :
471 [ + + + + ]: 5276180 : cg(:,:) = tdks%cg(:,:) !Psi(t)
472 : :
473 : : !** Predictor step
474 : : ! predict psi(t+dt) using ER propagator
475 : 20 : call rttddft_propagator_er(dtset,istep,mpi_enreg,psps,tdks,calc_properties=.true.)
476 : : ! for convergence check
477 [ + + + + ]: 5276180 : diff = tdks%cg
478 : : ! estimate psi(t+dt/2) = (psi(t)+psi(t+dt))/2
479 [ + + + + ]: 5276180 : tdks%cg(:,:) = 0.5_dp*(tdks%cg(:,:)+cg(:,:))
480 : : ! update electric field and vector potential to t+dt/2
481 : : ! We decide here not to update the induced vector potential to t+dt/2 for simplicity
482 : : call tdks%tdef%update(dtset,mpi_enreg,(istep-0.5_dp)*tdks%dt,tdks%rprimd,tdks%gprimd,tdks%kg, &
483 : 20 : & psps%mpsang,tdks%npwarr,tdks%ylm,tdks%ylmgr,tdks%current,update_vecpot_ind=.false.)
484 : : ! calc associated density at t+dt/2
485 : 20 : call rttddft_calc_density(dtset,mpi_enreg,psps,tdks)
486 : : ! go back to time t ..
487 [ + + + + ]: 5276180 : tdks%cg(:,:) = cg(:,:)
488 : : ! .. and evolve psi(t) using the EMR propagator with the estimated density at t+dt/2
489 : 20 : call rttddft_propagator_er(dtset,istep,mpi_enreg,psps,tdks)
490 : :
491 : : ! check convergence
492 [ + + + + ]: 5276180 : diff = abs(diff-tdks%cg)
493 : 20 : me = xmpi_comm_rank(mpi_enreg%comm_world)
494 [ + - - + : 1758760 : call xmpi_max(maxval(diff(1,:)),max_diff(1),mpi_enreg%comm_world,ierr)
- - + + ]
495 [ + - - + : 1758760 : call xmpi_max(maxval(diff(2,:)),max_diff(2),mpi_enreg%comm_world,ierr)
- - + + ]
496 [ + + - + ]: 20 : lconv = (max_diff(1) < dtset%td_scthr .and. max_diff(2) < dtset%td_scthr)
497 : 20 : ics = 0
498 [ + - ]: 20 : if (mpi_enreg%me == 0) then
499 : 20 : write(msg,'(a,a,i3,a,3(es8.2,1x),l1,a)') ch10, 'SC Step', ics, ' - ', max_diff(1), max_diff(2), &
500 : 40 : & dtset%td_scthr, lconv, ch10
501 [ + - ]: 20 : if (do_write_log) call wrtout(std_out,msg)
502 : : end if
503 [ + + ]: 20 : if (.not. lconv) then
504 : : !** Corrector steps
505 [ + - ]: 10 : do ics = 1, dtset%td_scnmax
506 : : ! for convergence check
507 [ + + + + ]: 1451530 : diff = tdks%cg
508 : : ! estimate psi(t+dt/2) = (psi(t)+psi(t+dt))/2
509 [ + + + + ]: 1451530 : tdks%cg(:,:) = 0.5_dp*(tdks%cg(:,:)+cg(:,:))
510 : : ! calc associated density at t+dt/2
511 : 10 : call rttddft_calc_density(dtset,mpi_enreg,psps,tdks)
512 : : ! Go back to time t ..
513 [ + + + + ]: 1451530 : tdks%cg(:,:) = cg(:,:)
514 : : ! .. and evolve psi(t) using estimated density at t+dt/2
515 : 10 : call rttddft_propagator_er(dtset,istep,mpi_enreg,psps,tdks)
516 : : ! check convergence
517 [ + + + + ]: 1451530 : diff = abs(diff-tdks%cg)
518 : 10 : me = xmpi_comm_rank(mpi_enreg%comm_world)
519 [ + - - + : 483860 : call xmpi_max(maxval(diff(1,:)),max_diff(1),mpi_enreg%comm_world,ierr)
- - + + ]
520 [ + - - + : 483860 : call xmpi_max(maxval(diff(2,:)),max_diff(2),mpi_enreg%comm_world,ierr)
- - + + ]
521 [ + - - + ]: 10 : lconv = (max_diff(1) < dtset%td_scthr .and. max_diff(2) < dtset%td_scthr)
522 [ + - ]: 10 : if (mpi_enreg%me == 0) then
523 : 10 : write(msg,'(a,a,i3,a,3(es8.2,1x),l1,a)') ch10, 'SC Step', ics, ' - ', max_diff(1), max_diff(2), &
524 : 20 : & dtset%td_scthr, lconv, ch10
525 [ + - ]: 10 : if (do_write_log) call wrtout(std_out,msg)
526 : : end if
527 [ - + ]: 10 : if (lconv) exit
528 : : end do
529 : : end if
530 : :
531 [ + - ]: 20 : if (lconv) then
532 : 20 : write(msg,'(a,i4,a)') "Converged after ", ics, " self-consistent corrector steps"
533 : 20 : call wrtout(ab_out,msg)
534 [ + - ]: 20 : if (do_write_log) call wrtout(std_out,msg)
535 : : else
536 : 0 : write(msg,'(a)') "Reached maximum number of corrector steps before convergence!"
537 : 0 : ABI_WARNING(msg)
538 : 0 : call wrtout(ab_out,msg)
539 : : end if
540 : :
541 : 20 : end subroutine rttddft_propagator_emr
542 : : !!***
543 : :
544 : : end module m_rttddft_propagators
545 : : !!***
|
