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