TABLE OF CONTENTS
ABINIT/m_rttddft_propagators [ Modules ]
NAME
m_rttddft_propagators
FUNCTION
Contains various propagators for the KS orbitals
COPYRIGHT
Copyright (C) 2021-2024 ABINIT group (FB) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
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, energies_init, energies_copy 33 use m_gemm_nonlop, only: make_gemm_nonlop 34 use m_hamiltonian, only: gs_hamiltonian_type, gspot_transgrid_and_pack 35 use m_invovl, only: make_invovl 36 use m_kg, only: mkkin, mkkpg 37 use m_mkffnl, only: mkffnl 38 use m_mpinfo, only: proc_distrb_cycle 39 use m_profiling_abi, only: abimem_record 40 use m_rttddft, only: rttddft_init_hamiltonian 41 use m_rttddft_exponential, only: rttddft_exp_taylor 42 use m_rttddft_properties, only: rttddft_calc_density, & 43 & rttddft_calc_occ, & 44 & rttddft_calc_kin 45 use m_rttddft_tdks, only: tdks_type 46 use m_specialmsg, only: wrtout 47 use m_xmpi, only: xmpi_comm_rank, xmpi_sum, xmpi_max 48 49 implicit none 50 51 private
m_rttddft/rttddft_propagator_emr [ Functions ]
[ Top ] [ m_rttddft ] [ Functions ]
NAME
rttddft_propagator_emr
FUNCTION
Main subroutine to propagate the KS orbitals using the Exponential Midpoint Rule (EMR) propagator
INPUTS
dtset <type(dataset_type)>=all input variables for this dataset ham_k <type(gs_hamiltonian_type)> = Hamiltonian object istep <integer> = step number mpi_enreg <MPI_type> = MPI-parallelisation information psps <type(pseudopotential_type)>=variables related to pseudopotentials tdks <type(tdks_type)> = the tdks object to initialize
OUTPUT
NOTES
This propagator is time reversible (if H(t+dt/2) and the exponential are computed exactly).
SOURCE
451 subroutine rttddft_propagator_emr(dtset, ham_k, istep, mpi_enreg, psps, tdks) 452 453 implicit none 454 455 !Arguments ------------------------------------ 456 !scalars 457 integer, intent(in) :: istep 458 type(dataset_type), intent(inout) :: dtset 459 type(gs_hamiltonian_type), intent(inout) :: ham_k 460 type(MPI_type), intent(inout) :: mpi_enreg 461 type(pseudopotential_type), intent(inout) :: psps 462 type(tdks_type), intent(inout) :: tdks 463 464 !Local variables------------------------------- 465 integer :: me 466 !scalars 467 character(len=500) :: msg 468 integer :: ics 469 integer :: ierr 470 logical :: lconv 471 !arrays 472 real(dp) :: cg(SIZE(tdks%cg(:,1)),SIZE(tdks%cg(1,:))) 473 real(dp) :: diff(SIZE(tdks%cg(:,1)),SIZE(tdks%cg(1,:))) 474 real(dp) :: max_diff(2) 475 476 ! *********************************************************************** 477 478 cg(:,:) = tdks%cg(:,:) !Psi(t) 479 480 !** Predictor step 481 ! predict psi(t+dt) using ER propagator 482 call rttddft_propagator_er(dtset,ham_k,istep,mpi_enreg,psps,tdks,calc_properties=.true.) 483 ! for convergence check 484 diff = tdks%cg 485 ! estimate psi(t+dt/2) = (psi(t)+psi(t+dt))/2 486 tdks%cg(:,:) = 0.5_dp*(tdks%cg(:,:)+cg(:,:)) 487 ! calc associated density at t+dt/2 488 call rttddft_calc_density(dtset,mpi_enreg,psps,tdks) 489 ! go back to time t .. 490 tdks%cg(:,:) = cg(:,:) 491 ! .. and evolve psi(t) using the EMR propagator with the estimated density at t+dt/2 492 call rttddft_propagator_er(dtset,ham_k,istep,mpi_enreg,psps,tdks) 493 494 ! check convergence 495 diff = abs(diff-tdks%cg) 496 me = xmpi_comm_rank(mpi_enreg%comm_world) 497 call xmpi_max(maxval(diff(1,:)),max_diff(1),mpi_enreg%comm_world,ierr) 498 call xmpi_max(maxval(diff(2,:)),max_diff(2),mpi_enreg%comm_world,ierr) 499 lconv = (max_diff(1) < dtset%td_scthr .and. max_diff(2) < dtset%td_scthr) 500 ics = 0 501 if (mpi_enreg%me == 0) then 502 write(msg,'(a,a,i3,a,3(es8.2,1x),l1,a)') ch10, 'SC Step', ics, ' - ', max_diff(1), max_diff(2), & 503 & dtset%td_scthr, lconv, ch10 504 if (do_write_log) call wrtout(std_out,msg) 505 end if 506 if (.not. lconv) then 507 !** Corrector steps 508 do ics = 1, dtset%td_scnmax 509 ! for convergence check 510 diff = tdks%cg 511 ! estimate psi(t+dt/2) = (psi(t)+psi(t+dt))/2 512 tdks%cg(:,:) = 0.5_dp*(tdks%cg(:,:)+cg(:,:)) 513 ! calc associated density at t+dt/2 514 call rttddft_calc_density(dtset,mpi_enreg,psps,tdks) 515 ! Go back to time t .. 516 tdks%cg(:,:) = cg(:,:) 517 ! .. and evolve psi(t) using estimated density at t+dt/2 518 call rttddft_propagator_er(dtset,ham_k,istep,mpi_enreg,psps,tdks) 519 ! check convergence 520 diff = abs(diff-tdks%cg) 521 me = xmpi_comm_rank(mpi_enreg%comm_world) 522 call xmpi_max(maxval(diff(1,:)),max_diff(1),mpi_enreg%comm_world,ierr) 523 call xmpi_max(maxval(diff(2,:)),max_diff(2),mpi_enreg%comm_world,ierr) 524 lconv = (max_diff(1) < dtset%td_scthr .and. max_diff(2) < dtset%td_scthr) 525 if (mpi_enreg%me == 0) then 526 write(msg,'(a,a,i3,a,3(es8.2,1x),l1,a)') ch10, 'SC Step', ics, ' - ', max_diff(1), max_diff(2), & 527 & dtset%td_scthr, lconv, ch10 528 if (do_write_log) call wrtout(std_out,msg) 529 end if 530 if (lconv) exit 531 end do 532 end if 533 534 if (lconv) then 535 write(msg,'(a,i4,a)') "Converged after ", ics, " self-consistent corrector steps" 536 call wrtout(ab_out,msg) 537 if (do_write_log) call wrtout(std_out,msg) 538 else 539 write(msg,'(a)') "Reached maximum number of corrector steps before convergence" 540 call wrtout(ab_out,msg) 541 if (do_write_log) call wrtout(std_out,msg) 542 end if 543 544 end subroutine rttddft_propagator_emr
m_rttddft/rttddft_propagator_er [ Functions ]
[ Top ] [ m_rttddft ] [ Functions ]
NAME
rttddft_propagator_er
FUNCTION
Main subroutine to propagate the KS orbitals using the Exponential Rule (ER) propagator.
INPUTS
dtset <type(dataset_type)> = all input variables for this dataset ham_k <type(gs_hamiltonian_type)> = Hamiltonian object istep <integer> = step number mpi_enreg <MPI_type> = MPI-parallelisation information psps <type(pseudopotential_type)>=variables related to pseudopotentials tdks <type(tdks_type)> = the tdks object to initialize calc_properties <logical> = logical governing the computation of some properties (energies, occupations, eigenvalues..)
NOTES
Other propagators such as the Exponential Midpoint Rule (EMR) should usually be prefered over this one since the ER propagator alone violates time reversal symmetry. Using this propagator with the exponential approximated by Taylor expansion of order 1 leads to the famous Euler method which is fast and simple but unstable and thus insufficient for RT-TDDFT.
SOURCE
89 subroutine rttddft_propagator_er(dtset, ham_k, istep, mpi_enreg, psps, tdks, calc_properties) 90 91 implicit none 92 93 !Arguments ------------------------------------ 94 !scalars 95 integer, intent(in) :: istep 96 logical, optional, intent(in) :: calc_properties 97 type(dataset_type), intent(inout) :: dtset 98 type(gs_hamiltonian_type), intent(inout) :: ham_k 99 type(MPI_type), intent(inout) :: mpi_enreg 100 type(pseudopotential_type), intent(inout) :: psps 101 type(tdks_type), target, intent(inout) :: tdks 102 103 !Local variables------------------------------- 104 !scalars 105 integer :: bdtot_index 106 integer :: calc_forces 107 integer :: dimffnl 108 integer :: gemm_nonlop_ikpt_this_proc_being_treated 109 integer :: iband 110 integer :: ibg, icg 111 integer :: ider, idir 112 integer :: ierr, ilm 113 integer :: ikpt, ikpt_loc, ikg 114 integer :: isppol 115 integer :: istwf_k 116 integer :: me_distrb 117 integer :: me_bandfft 118 integer :: my_ikpt, my_nspinor 119 integer :: nband_k, nband_k_mem 120 integer :: npw_k, nkpg 121 integer :: shift 122 integer :: spaceComm_distrb 123 integer :: n4, n5, n6 124 logical :: with_vxctau 125 logical :: lcalc_properties 126 type(energies_type) :: energies 127 type(bandfft_kpt_type),pointer :: my_bandfft_kpt => null() 128 !arrays 129 integer, allocatable :: kg_k(:,:) 130 real(dp), pointer :: cg(:,:) => null() 131 real(dp), pointer :: cg0(:,:) => null() 132 real(dp), allocatable :: enl(:) 133 real(dp), pointer :: eig(:) => null() 134 real(dp), allocatable :: ffnl(:,:,:,:) 135 real(dp), allocatable :: kpg_k(:,:) 136 real(dp) :: kpoint(3) 137 real(dp), allocatable :: kinpw(:) 138 real(dp), pointer :: occ(:) => null() 139 real(dp), pointer :: occ0(:) => null() 140 real(dp), allocatable :: ph3d(:,:,:) 141 real(dp), allocatable :: vlocal(:,:,:,:) 142 real(dp), allocatable :: vxctaulocal(:,:,:,:,:) 143 real(dp), allocatable :: ylm_k(:,:) 144 logical :: lproperties(4) 145 146 ! *********************************************************************** 147 148 !Init MPI 149 spaceComm_distrb=mpi_enreg%comm_cell 150 if (mpi_enreg%paral_kgb==1) spaceComm_distrb=mpi_enreg%comm_kpt 151 me_distrb=xmpi_comm_rank(spaceComm_distrb) 152 153 !Do we calculate properties? 154 !Governed by lproperties: 155 ! lproperties(1) = compute energy contributions (kinetic) 156 ! lproperties(2) = NL energy contribution in NC case 157 ! lproperties(3) = eigenvalues 158 ! lproperties(4) = occupations 159 lproperties(:) = .false. 160 lcalc_properties = .false. 161 if (present(calc_properties)) then 162 lcalc_properties = calc_properties 163 if (lcalc_properties) then 164 !compute energy contributions 165 lproperties(1) = .true. 166 !Init to zero different energies 167 call energies_init(energies) 168 energies%entropy=tdks%energies%entropy 169 energies%e_corepsp=tdks%energies%e_corepsp 170 energies%e_ewald=tdks%energies%e_ewald 171 !including NL part in NC case? 172 if (dtset%usepaw == 0) then 173 lproperties(2) = .true. 174 else 175 ABI_MALLOC(enl,(0)) 176 end if 177 !other properties 178 ! eigenvalues 179 if (dtset%prteig /= 0 .or. dtset%prtdos /= 0) then 180 if (mod(istep-1,dtset%td_prtstr) == 0) then 181 lproperties(3) = .true. 182 tdks%eigen(:) = zero 183 end if 184 end if 185 ! occupations 186 lproperties(4) = .true. 187 tdks%occ(:) = zero 188 end if 189 end if 190 191 !Set "vtrial" and initialize the Hamiltonian 192 call rttddft_init_hamiltonian(dtset,energies,ham_k,istep,mpi_enreg,psps,tdks) 193 194 my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor) 195 n4=dtset%ngfft(4); n5=dtset%ngfft(5); n6=dtset%ngfft(6) 196 ABI_MALLOC(vlocal,(n4,n5,n6,ham_k%nvloc)) 197 with_vxctau=(dtset%usekden/=0) 198 if(with_vxctau) then 199 ABI_MALLOC(vxctaulocal,(n4,n5,n6,ham_k%nvloc,4)) 200 end if 201 if (dtset%ionmov/=0) then 202 calc_forces=1 203 else 204 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 icg=0; ibg=0 215 bdtot_index=0 216 217 !*** LOOP OVER SPINS 218 do isppol=1,dtset%nsppol 219 220 ikpt_loc=0 221 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 & dtset%nspden, ham_k%nvloc, 1, tdks%pawfgr, mpi_enreg, tdks%vtrial, vlocal) 227 call ham_k%load_spin(isppol, vlocal=vlocal, with_nonlocal=.true.) 228 229 if (with_vxctau) then 230 call gspot_transgrid_and_pack(isppol, psps%usepaw, dtset%paral_kgb, dtset%nfft, dtset%ngfft, tdks%nfftf, & 231 & dtset%nspden, ham_k%nvloc, 4, tdks%pawfgr, mpi_enreg, tdks%vxctau, vxctaulocal) 232 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 ikpt = 0 251 do while (ikpt_loc < dtset%nkpt) 252 253 ikpt_loc = ikpt_loc + 1 254 ikpt = ikpt_loc 255 my_ikpt = mpi_enreg%my_kpttab(ikpt) 256 257 nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt) 258 if (mpi_enreg%paral_kgb==1) then 259 nband_k_mem=mpi_enreg%bandpp 260 else 261 nband_k_mem=nband_k 262 end if 263 istwf_k=dtset%istwfk(ikpt) 264 npw_k=tdks%npwarr(ikpt) 265 266 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) then 267 bdtot_index=bdtot_index+nband_k 268 cycle 269 end if 270 271 if (mpi_enreg%paral_kgb==1) my_bandfft_kpt => bandfft_kpt(my_ikpt) 272 call bandfft_kpt_set_ikpt(ikpt,mpi_enreg) 273 274 ABI_MALLOC(kg_k,(3,npw_k)) 275 ABI_MALLOC(ylm_k,(npw_k,psps%mpsang*psps%mpsang*psps%useylm)) 276 kg_k(:,1:npw_k)=tdks%kg(:,1+ikg:npw_k+ikg) 277 if (psps%useylm==1) then 278 do ilm=1,psps%mpsang*psps%mpsang 279 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 ! Kinetic energy - Compute (1/2) (2 Pi)**2 (k+G)**2: 285 kpoint(:)=dtset%kptns(:,ikpt) 286 ABI_MALLOC(kinpw,(npw_k)) 287 call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,tdks%gmet,kg_k,kinpw,kpoint,npw_k,0,0) 288 ! Compute (k+G) vectors (only if useylm=1) 289 nkpg=3*calc_forces*dtset%nloalg(3) 290 ABI_MALLOC(kpg_k,(npw_k,nkpg)) 291 if ((mpi_enreg%paral_kgb/=1.or.istep<=tdks%first_step).and.nkpg>0) call mkkpg(kg_k,kpg_k,kpoint,nkpg,npw_k) 292 ! Compute nonlocal form factors ffnl at all (k+G): 293 ider=0;idir=0;dimffnl=1 294 ABI_MALLOC(ffnl,(npw_k,dimffnl,psps%lmnmax,psps%ntypat)) 295 if (mpi_enreg%paral_kgb/=1.or.istep<=tdks%first_step) then 296 call mkffnl(psps%dimekb,dimffnl,psps%ekb,ffnl,psps%ffspl,tdks%gmet,tdks%gprimd, & 297 & ider,idir,psps%indlmn,kg_k,kpg_k,kpoint,psps%lmnmax,psps%lnmax, & 298 & psps%mpsang,psps%mqgrid_ff,nkpg,npw_k,psps%ntypat,psps%pspso, & 299 & psps%qgrid_ff,tdks%rmet,psps%usepaw,psps%useylm,ylm_k,tdks%ylmgr) 300 end if 301 302 !** Load k-dependent part in the Hamiltonian datastructure 303 !** - Compute 3D phase factors 304 !** - Prepare various tabs in case of band-FFT parallelism 305 !** - Load k-dependent quantities in the Hamiltonian 306 ABI_MALLOC(ph3d,(2,npw_k,ham_k%matblk)) 307 call ham_k%load_k(kpt_k=dtset%kptns(:,ikpt),istwf_k=istwf_k,npw_k=npw_k,kinpw_k=kinpw,kg_k=kg_k,kpg_k=kpg_k, & 308 & ffnl_k=ffnl,ph3d_k=ph3d,compute_ph3d=(mpi_enreg%paral_kgb/=1.or.istep<=tdks%first_step), & 309 & compute_gbound=(mpi_enreg%paral_kgb/=1)) 310 311 !** Load band-FFT tabs (transposed k-dependent arrays) 312 if (mpi_enreg%paral_kgb==1) then 313 if (istep<=tdks%first_step) call prep_bandfft_tabs(ham_k,ikpt,dtset%mkmem,mpi_enreg) 314 call ham_k%load_k(npw_fft_k =my_bandfft_kpt%ndatarecv, & 315 & gbound_k =my_bandfft_kpt%gbound, & 316 & kinpw_k =my_bandfft_kpt%kinpw_gather, & 317 & kg_k =my_bandfft_kpt%kg_k_gather, & 318 & kpg_k =my_bandfft_kpt%kpg_k_gather, & 319 & ffnl_k =my_bandfft_kpt%ffnl_gather, & 320 & ph3d_k =my_bandfft_kpt%ph3d_gather) 321 end if 322 323 !** Build inverse of overlap matrix 324 if(psps%usepaw == 1 .and. istep <= tdks%first_step) then 325 call make_invovl(ham_k, dimffnl, ffnl, ph3d, mpi_enreg) 326 end if 327 328 ! Setup gemm_nonlop 329 if (tdks%gemm_nonlop_use_gemm) then 330 !set the global variable indicating to gemm_nonlop where to get its data from 331 gemm_nonlop_ikpt_this_proc_being_treated = my_ikpt 332 if (istep <= tdks%first_step) then 333 !Init the arrays 334 call make_gemm_nonlop(my_ikpt,ham_k%npw_fft_k,ham_k%lmnmax, & 335 ham_k%ntypat, ham_k%indlmn, ham_k%nattyp, ham_k%istwf_k, & 336 ham_k%ucvol, ham_k%ffnl_k, & 337 ham_k%ph3d_k, ham_k%kpt_k, ham_k%kg_k, ham_k%kpg_k) 338 end if 339 end if 340 341 !** Compute the exp[(S^{-1})H]*cg using Taylor expansion to approximate the exponential 342 cg => tdks%cg(:,icg+1:icg+nband_k*npw_k*my_nspinor) 343 ! Compute properties "on-the-fly" if required 344 if (lcalc_properties) then 345 cg0 => tdks%cg0(:,icg+1:icg+nband_k*npw_k*my_nspinor) 346 occ => tdks%occ(bdtot_index+1:bdtot_index+nband_k) 347 occ0 => tdks%occ0(bdtot_index+1:bdtot_index+nband_k) 348 if (dtset%paral_kgb /= 1) then 349 shift = bdtot_index 350 else 351 me_bandfft = xmpi_comm_rank(mpi_enreg%comm_band) 352 shift = bdtot_index+me_bandfft*mpi_enreg%bandpp 353 end if 354 ! kinetic energy 355 if (lproperties(1)) then 356 call rttddft_calc_kin(energies%e_kinetic,cg,dtset,ham_k,nband_k,npw_k,my_nspinor, & 357 & occ0,dtset%wtk(ikpt),mpi_enreg,my_bandfft_kpt) 358 end if 359 ! for NL PSP part 360 if (lproperties(2)) then 361 ABI_MALLOC(enl,(mpi_enreg%bandpp)) 362 enl = zero 363 end if 364 ! for eigenvalues 365 if (lproperties(3)) then 366 eig => tdks%eigen(1+shift:nband_k_mem+shift) 367 end if 368 ! occupations 369 if (lproperties(4)) then 370 !note that occupations are computed at istep-1 like energies and eigenvalues 371 call rttddft_calc_occ(cg,cg0,dtset,ham_k,ikpt,ibg,isppol,mpi_enreg, & 372 & nband_k,npw_k,my_nspinor,occ,occ0,tdks) 373 end if 374 375 !Propagate cg and compute the requested properties 376 call rttddft_exp_taylor(cg,dtset,ham_k,mpi_enreg,nband_k,npw_k,my_nspinor,enl=enl,eig=eig) 377 378 !Finish computing NL PSP part for NC PSP 379 if (lproperties(2)) then 380 do iband = 1, mpi_enreg%bandpp 381 energies%e_nlpsp_vfock=energies%e_nlpsp_vfock+dtset%wtk(ikpt)*tdks%occ0(shift+iband)*enl(iband) 382 end do 383 ABI_FREE(enl) 384 end if 385 else 386 !Propagate cg only 387 call rttddft_exp_taylor(cg,dtset,ham_k,mpi_enreg,nband_k,npw_k,my_nspinor) 388 end if 389 390 ABI_FREE(kg_k) 391 ABI_FREE(ylm_k) 392 ABI_FREE(kpg_k) 393 ABI_FREE(kinpw) 394 ABI_FREE(ffnl) 395 ABI_FREE(ph3d) 396 397 bdtot_index = bdtot_index+nband_k 398 399 !** Also shift array memory if dtset%mkmem/=0 400 if (dtset%mkmem/=0) then 401 ibg=ibg+nband_k_mem 402 icg=icg+npw_k*my_nspinor*nband_k 403 ikg=ikg+npw_k 404 end if 405 406 end do !nkpt 407 408 end do !nsppol 409 410 ABI_FREE(vlocal) 411 if(dtset%usekden/=0) then 412 ABI_FREE(vxctaulocal) 413 end if 414 415 !Keep the computed energies in memory 416 if (lcalc_properties) then 417 call xmpi_sum(energies%e_kinetic,mpi_enreg%comm_kptband,ierr) 418 call xmpi_sum(energies%e_nlpsp_vfock,mpi_enreg%comm_kptband,ierr) 419 call energies_copy(energies,tdks%energies) 420 if (lproperties(3)) call xmpi_sum(tdks%eigen,mpi_enreg%comm_kptband,ierr) 421 if (lproperties(4)) call xmpi_sum(tdks%occ,mpi_enreg%comm_kpt,ierr) 422 end if 423 424 end subroutine rttddft_propagator_er