Line data Source code
1 : !!****m* ABINIT/m_rttddft_properties
2 : !! NAME
3 : !! m_rttddft_properties
4 : !!
5 : !! FUNCTION
6 : !! Contains most of the subroutines to compute
7 : !! properties (energy, occupations, eigenvalues..)
8 : !!
9 : !! COPYRIGHT
10 : !! Copyright (C) 2021-2026 ABINIT group (FB)
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_properties
24 :
25 : use defs_basis
26 : use defs_abitypes, only: MPI_type
27 : use defs_datatypes, only: pseudopotential_type
28 :
29 : use m_bandfft_kpt, only: bandfft_kpt_type
30 : use m_cgprj, only: ctocprj
31 : use m_cgtools, only: dotprod_g
32 : use m_dtfil, only: datafiles_type
33 : use m_dtset, only: dataset_type
34 : use m_energies, only: energies_type
35 : use m_fourier_interpol, only: transgrid
36 : use m_hamiltonian, only: gs_hamiltonian_type
37 : use m_mkrho, only: mkrho
38 : use m_nonlop, only: nonlop
39 : use m_pawcprj, only: pawcprj_type, pawcprj_alloc, pawcprj_get, &
40 : & pawcprj_free, pawcprj_mpi_allgather
41 : use m_paw_optics, only: optics_paw
42 : use m_paw_mkrho, only: pawmkrho
43 : use m_paw_occupancies, only: pawmkrhoij
44 : use m_pawrhoij, only: pawrhoij_type, pawrhoij_free, &
45 : & pawrhoij_alloc, pawrhoij_inquire_dim
46 : use m_profiling_abi, only: abimem_record
47 : use m_rttddft_tdks, only: tdks_type
48 : use m_spacepar, only: meanvalue_g
49 : use m_xmpi, only: xmpi_sum, xmpi_comm_rank
50 : use m_dft_energy, only: entropy
51 :
52 : implicit none
53 :
54 : private
55 : !!***
56 :
57 : public :: rttddft_calc_density
58 : public :: rttddft_calc_etot
59 : public :: rttddft_calc_eig
60 : public :: rttddft_calc_enl
61 : public :: rttddft_calc_kin
62 : public :: rttddft_calc_occ
63 : public :: rttddft_calc_current
64 : !!***
65 :
66 : contains
67 : !!***
68 :
69 : !!****f* m_rttddft_properties/rttddft_calc_density
70 : !!
71 : !! NAME
72 : !! rttddft_calc_density
73 : !!
74 : !! FUNCTION
75 : !! Compute electronic density (in 1/bohr^3) from the WF (cg coefficients)
76 : !!
77 : !! INPUTS
78 : !! dtset <type(dataset_type)> = all input variables for this dataset
79 : !! mpi_enreg <MPI_type> = MPI-parallelisation information
80 : !! psps <type(pseudopotential_type)> = variables related to pseudopotentials
81 : !! tdks <type(tdks_type)> = Main RT-TDDFT object
82 : !!
83 : !! OUTPUT
84 : !!
85 : !! SOURCE
86 1805 : subroutine rttddft_calc_density(dtset, mpi_enreg, psps, tdks)
87 :
88 : !Arguments ------------------------------------
89 : !scalars
90 : type(tdks_type), intent(inout) :: tdks
91 : type(dataset_type), intent(in) :: dtset
92 : type(MPI_type), intent(inout) :: mpi_enreg
93 : type(pseudopotential_type), intent(in) :: psps
94 :
95 : !Local variables-------------------------------
96 : !scalars
97 : integer, parameter :: cplex=1
98 : integer :: cplex_rhoij
99 : integer :: ipert, idir
100 : integer :: my_natom
101 : integer :: nspden_rhoij
102 : integer :: tim_mkrho
103 : real(dp) :: compch_fft
104 : !arrays
105 : real(dp) :: qpt(3)
106 1805 : real(dp),allocatable :: rhowfg(:,:), rhowfr(:,:)
107 1805 : type(pawrhoij_type),pointer :: pawrhoij_unsym(:)
108 : logical :: printout
109 :
110 : ! ***********************************************************************
111 :
112 1805 : my_natom=mpi_enreg%my_natom
113 :
114 1805 : tim_mkrho=1
115 :
116 1805 : printout = (dtset%prtvol /= 0)
117 :
118 1805 : if (psps%usepaw==1) then
119 :
120 4725 : ABI_MALLOC(rhowfg,(2,dtset%nfft))
121 6300 : ABI_MALLOC(rhowfr,(dtset%nfft,dtset%nspden))
122 :
123 : ! 1-Compute density from WFs (without compensation charge density nhat)
124 : call mkrho(tdks%cg,dtset,tdks%gprimd,tdks%irrzon,tdks%kg,tdks%mcg,mpi_enreg, &
125 : & tdks%npwarr,tdks%occ0,tdks%paw_dmft,tdks%phnons,rhowfg,rhowfr, &
126 : & tdks%rprimd,tim_mkrho,tdks%ucvol,tdks%wvl%den,tdks%wvl%wfs, &
127 1575 : & printout=printout)
128 :
129 : ! 2-Compute cprj = <\psi_{n,k}|p_{i,j}>
130 : call ctocprj(tdks%atindx,tdks%cg,1,tdks%cprj,tdks%gmet,tdks%gprimd,0,0,0, &
131 : & dtset%istwfk,tdks%kg,tdks%tdef%kpa,tdks%mcg,tdks%mcprj,dtset%mgfft, &
132 : & dtset%mkmem,mpi_enreg,psps%mpsang,dtset%mpw,dtset%natom,tdks%nattyp, &
133 : & dtset%nband,dtset%natom,dtset%ngfft,dtset%nkpt,dtset%nloalg, &
134 : & tdks%npwarr,dtset%nspinor,dtset%nsppol,dtset%nsppol,psps%ntypat, &
135 : & dtset%paral_kgb,tdks%ph1d,psps,tdks%rmet,dtset%typat,tdks%ucvol, &
136 1575 : & tdks%unpaw,tdks%xred,tdks%ylm,tdks%ylmgr)
137 :
138 : !paral atom
139 1575 : if (my_natom/=dtset%natom) then
140 0 : ABI_MALLOC(pawrhoij_unsym,(dtset%natom))
141 : call pawrhoij_inquire_dim(cplex_rhoij=cplex_rhoij,nspden_rhoij=nspden_rhoij, &
142 : & nspden=dtset%nspden,spnorb=dtset%pawspnorb, &
143 0 : & cpxocc=dtset%pawcpxocc)
144 : call pawrhoij_alloc(pawrhoij_unsym,cplex_rhoij,nspden_rhoij,dtset%nspinor, &
145 0 : & dtset%nsppol,dtset%typat,pawtab=tdks%pawtab,use_rhoijp=0)
146 : else
147 1575 : pawrhoij_unsym => tdks%pawrhoij
148 : end if
149 :
150 : ! 3-Compute pawrhoij = \rho_{i,j} = \sum_{n,k}f_{n,k} \tilde{c}^{i,*}_{n,k} \tilde{c}^{j}_{n,k}
151 : call pawmkrhoij(tdks%atindx,tdks%atindx1,tdks%cprj,tdks%dimcprj_srt,dtset%istwfk, &
152 : & dtset%kptopt,dtset%mband,tdks%mband_cprj,tdks%mcprj,dtset%mkmem, &
153 : & mpi_enreg,dtset%natom,dtset%nband,dtset%nkpt,dtset%nspden, &
154 : & dtset%nspinor,dtset%nsppol,tdks%occ0,dtset%paral_kgb,tdks%paw_dmft, &
155 1575 : & pawrhoij_unsym,tdks%unpaw,dtset%usewvl,dtset%wtk)
156 :
157 : ! 4-Symetrize rhoij, compute nhat and add it to rhor
158 : ! Note pawrhoij_unsym and pawrhoij are the same, which means that pawrhoij
159 : ! cannot be distributed over different atomic sites.
160 1575 : ipert=0; idir=0; qpt(:)=zero; compch_fft=-1e-5_dp
161 24652719 : tdks%nhat = zero
162 : call pawmkrho(1,compch_fft,cplex,tdks%gprimd,idir,tdks%indsym,ipert,mpi_enreg, &
163 : & my_natom,dtset%natom,dtset%nspden,dtset%nsym,dtset%ntypat, &
164 : & dtset%paral_kgb,tdks%pawang,tdks%pawfgr,tdks%pawfgrtab, &
165 : & dtset%pawprtvol,tdks%pawrhoij,pawrhoij_unsym,tdks%pawtab,qpt, &
166 : & rhowfg,rhowfr,tdks%rhor,tdks%rprimd,dtset%symafm,tdks%symrec, &
167 : & dtset%typat,tdks%ucvol,dtset%usewvl,tdks%xred,pawnhat=tdks%nhat, &
168 1575 : & rhog=tdks%rhog)
169 :
170 : ! 5-Take care of kinetic energy density
171 1575 : if(dtset%usekden==1)then
172 : call mkrho(tdks%cg,dtset,tdks%gprimd,tdks%irrzon,tdks%kg,tdks%mcg,mpi_enreg, &
173 : & tdks%npwarr,tdks%occ0,tdks%paw_dmft,tdks%phnons,rhowfg,rhowfr, &
174 0 : & tdks%rprimd,tim_mkrho,tdks%ucvol,tdks%wvl%den,tdks%wvl%wfs,option=1)
175 :
176 : !FB: Useful?
177 : call transgrid(1,mpi_enreg,dtset%nspden,+1,1,1,dtset%paral_kgb,tdks%pawfgr, &
178 0 : & rhowfg,tdks%taug,rhowfr,tdks%taur)
179 : end if
180 :
181 1575 : ABI_FREE(rhowfg)
182 1575 : ABI_FREE(rhowfr)
183 :
184 1575 : if (my_natom/=dtset%natom) then
185 0 : call pawrhoij_free(pawrhoij_unsym)
186 0 : ABI_FREE(pawrhoij_unsym)
187 : else
188 : pawrhoij_unsym => NULL()
189 : end if
190 :
191 : else
192 :
193 : ! 1-Compute density from WFs
194 : call mkrho(tdks%cg,dtset,tdks%gprimd,tdks%irrzon,tdks%kg,tdks%mcg,mpi_enreg, &
195 : & tdks%npwarr,tdks%occ0,tdks%paw_dmft,tdks%phnons,tdks%rhog,tdks%rhor, &
196 230 : & tdks%rprimd,tim_mkrho,tdks%ucvol,tdks%wvl%den,tdks%wvl%wfs)
197 : ! 2-Take care of kinetic energy density
198 230 : if(dtset%usekden==1)then
199 : call mkrho(tdks%cg,dtset,tdks%gprimd,tdks%irrzon,tdks%kg,tdks%mcg,mpi_enreg, &
200 : & tdks%npwarr,tdks%occ0,tdks%paw_dmft,tdks%phnons,tdks%taug,tdks%taur, &
201 0 : & tdks%rprimd,tim_mkrho,tdks%ucvol,tdks%wvl%den,tdks%wvl%wfs,option=1)
202 : end if
203 :
204 : endif
205 :
206 1805 : end subroutine rttddft_calc_density
207 : !!***
208 :
209 : !!****f* m_rttddft_properties/rttddft_calc_etot
210 : !!
211 : !! NAME
212 : !! rttddft_calc_energy
213 : !!
214 : !! FUNCTION
215 : !! Computes total energy
216 : !!
217 : !! INPUTS
218 : !! dtset <type(dataset_type)> = all input variables for this dataset
219 : !! energies <energies_type> = contains the different contribution ot the total energy
220 : !! occ <real(nband*nkpt*nsspol)> = occupation numbers at time t
221 : !!
222 : !! OUTPUT
223 : !! etotal <real(dp)> = the total energy
224 : !!
225 : !! SOURCE
226 1725 : subroutine rttddft_calc_etot(dtset, energies, etotal, occ)
227 :
228 : !Arguments ------------------------------------
229 : !scalars
230 : type(dataset_type), intent(in) :: dtset
231 : type(energies_type), intent(inout) :: energies
232 : real(dp), intent(out) :: etotal
233 : !arrays
234 : real(dp), intent(in) :: occ(:)
235 :
236 : ! ***********************************************************************
237 :
238 : ! Compute electronic entropy
239 1725 : call rttddft_calc_ent(energies%entropy_ks, dtset, occ)
240 :
241 1725 : call entropy(dtset,energies)
242 :
243 : etotal = energies%e_kinetic &
244 : & + energies%e_hartree &
245 : & + energies%e_xc &
246 : & + energies%e_localpsp &
247 : & + energies%e_corepsp &
248 : & + energies%e_ewald &
249 : & + energies%e_vdw_dftd &
250 : & + energies%e_nlpsp_vfock &
251 1725 : & + energies%paw%epaw
252 : ! & + energies%e_entropy
253 : !FB: @MT Should one add the last e_paw contribution or not?
254 : !FB: Seeems like all the other contributions are not relevant here @MT?
255 : ! & + energies%e_chempot &
256 : ! & + energies%e_elecfield &
257 : ! & + energies%e_magfield &
258 : ! & + energies%e_nucdip &
259 : ! & + energies%e_hybcomp_E0 &
260 : ! & - energies%e_hybcomp_v0 &
261 : ! & + energies%e_hybcomp_v &
262 : ! & + energies%e_constrained_dft
263 :
264 : !if (psps%usepaw==0) etotal = etotal + energies%e_nlpsp_vfock - energies%e_fock0
265 : !if (psps%usepaw==1) etotal = etotal + energies%e_paw + energies%e_fock
266 :
267 1725 : end subroutine rttddft_calc_etot
268 : !!***
269 :
270 : !!****f* m_rttddft_properties/rttddft_calc_eig
271 : !!
272 : !! NAME
273 : !! rttddft_calc_eig
274 : !!
275 : !! FUNCTION
276 : !! Computes eigenvalues from cg and ghc = <G|H|C>
277 : !! and gsc = <G|S|C> if paw
278 : !!
279 : !! INPUTS
280 : !! cg <real(2,npw*nspinor*nband)> = the wavefunction coefficients
281 : !! ghc <real(2,npw*nspinor*nband)> = <G|H|C>
282 : !! istwf_k <integer> = option describing the storage of wfs at k
283 : !! nband <integer> = number of bands
284 : !! npw <integer> = number of plane waves
285 : !! nspinor <integer> = dimension of spinors
286 : !! me_g0 <integer> = if set to 1 current proc contains G(0,0,0)
287 : !! comm <integer> = MPI communicator
288 : !! gsc <real(2,npw*nspinor*nband)> = <G|S|C> (optional - only in PAW)
289 : !!
290 : !! OUTPUT
291 : !! eig <real(nband)> = the eigenvalues
292 : !!
293 : !! SOURCE
294 0 : subroutine rttddft_calc_eig(cg,eig,ghc,istwf_k,nband,npw,nspinor,me_g0,comm,gsc)
295 :
296 : !Arguments ------------------------------------
297 : !scalars
298 : integer, intent(in) :: istwf_k
299 : integer, intent(in) :: nband
300 : integer, intent(in) :: npw
301 : integer, intent(in) :: nspinor
302 : integer, intent(in) :: me_g0
303 : integer, intent(in) :: comm
304 : !arrays
305 : real(dp), intent(in) :: cg(2,npw*nspinor*nband)
306 : real(dp), intent(out) :: eig(nband)
307 : real(dp), intent(in) :: ghc(2,npw*nspinor*nband)
308 : real(dp), intent(in), optional :: gsc(2,npw*nspinor*nband)
309 :
310 : !Local variables-------------------------------
311 : !scalars
312 : integer :: iband
313 : integer :: shift
314 : real(dp) :: dprod_r, dprod_i
315 : !arrays
316 :
317 : ! ***********************************************************************
318 :
319 0 : do iband=1, nband
320 0 : shift = npw*nspinor*(iband-1)
321 : !Compute eigenvalues
322 : call dotprod_g(eig(iband),dprod_i,istwf_k,npw*nspinor,1, &
323 : & ghc(:, shift+1:shift+npw*nspinor), &
324 : & cg(:, shift+1:shift+npw*nspinor), &
325 0 : & me_g0, comm)
326 0 : if (present(gsc)) then
327 : call dotprod_g(dprod_r,dprod_i,istwf_k,npw*nspinor,1, &
328 : & gsc(:, shift+1:shift+npw*nspinor), &
329 : & cg(:, shift+1:shift+npw*nspinor), &
330 0 : & me_g0, comm)
331 0 : eig(iband) = eig(iband)/dprod_r
332 : end if
333 : end do
334 :
335 0 : end subroutine rttddft_calc_eig
336 : !!***
337 :
338 : !!****f* m_rttddft_properties/rttddft_calc_kin
339 : !!
340 : !! NAME
341 : !! rttddft_calc_kin
342 : !!
343 : !! FUNCTION
344 : !! Computes the NL part of energy in NC case
345 : !!
346 : !! INPUTS
347 : !! cg <real(2,npw*nspinor*nband)> = the wavefunction coefficients
348 : !! dtset <type(dataset_type)> = all input variables for this dataset
349 : !! ham_k <gs_hamiltonian_type> = hamiltonian at point k
350 : !! nband <integer> = number of bands
351 : !! npw <integer> = number of plane waves
352 : !! nspinor <integer> = dimension of spinors
353 : !! occ0 <real(nband)> = initial occupations
354 : !! wk <real> = weight of associated kpt
355 : !! mpi_enreg <MPI_type> = MPI-parallelisation information
356 : !! bandfft_kpt <bandfft_kpt_type> = additional info on parallelisation
357 : !! for the associated kpt
358 : !!
359 : !! OUTPUT
360 : !! kin <real(nband)> = the non local part of the energy in NC case
361 : !!
362 : !! SIDE EFFECTS
363 : !!
364 : !! SOURCE
365 27260 : subroutine rttddft_calc_kin(kin,cg,dtset,ham_k,nband,npw,nspinor,occ0,wk,mpi_enreg,bandfft)
366 :
367 : !Arguments ------------------------------------
368 : !scalars
369 : integer, intent(in) :: nband
370 : integer, intent(in) :: npw
371 : integer, intent(in) :: nspinor
372 : real(dp), intent(in) :: wk
373 : type(dataset_type), intent(in) :: dtset
374 : type(gs_hamiltonian_type), intent(in) :: ham_k
375 : type(MPI_type), intent(in) :: mpi_enreg
376 : type(bandfft_kpt_type), pointer, intent(in) :: bandfft
377 : !arrays
378 : real(dp), intent(in) :: cg(2,npw*nspinor*nband)
379 : real(dp), intent(inout) :: kin
380 : real(dp), intent(in) :: occ0(nband)
381 :
382 : !Local variables-------------------------------
383 : !scalars
384 : integer :: displ
385 : integer :: iband, ipw
386 : integer :: jpw
387 : integer :: me_bandfft
388 : integer :: shift
389 : real(dp) :: ar
390 : !arrays
391 :
392 : ! ***********************************************************************
393 :
394 27260 : if (dtset%paral_kgb /= 1) then
395 : displ = 0
396 : else
397 6560 : me_bandfft = xmpi_comm_rank(mpi_enreg%comm_bandspinorfft)
398 6560 : displ = bandfft%rdispls(me_bandfft+1)
399 : end if
400 :
401 255900 : do iband=1, nband
402 255900 : if (abs(occ0(iband))>tol8) then
403 115400 : shift = npw*nspinor*(iband-1)
404 : !FB: meanvalue_g does the mpi_sum over the bands inside, that's not very efficient since
405 : !FB: we could do it only once at the end
406 : !FB: From Lucas: meanvalue_g seems slow
407 : !FB: It maybe useful not to use meanvalue_g at all here
408 : !call meanvalue_g(ar,ham_k%kinpw_k(1+displ:displ+npw*nspinor),0,1,mpi_enreg,npw,nspinor, &
409 : ! & cg(:,1+shift:shift+npw*nspinor),cg(:,1+shift:shift+npw*nspinor),0)
410 115400 : ar = zero
411 24625780 : do ipw = 1, npw
412 24625780 : ar = ar + ham_k%kinpw_k(displ+ipw)*(cg(1,shift+ipw)*cg(1,shift+ipw)+cg(2,shift+ipw)*cg(2,shift+ipw))
413 : end do
414 115400 : if(nspinor==2)then
415 0 : do ipw = 1+npw, 2*npw
416 0 : jpw = ipw - npw
417 0 : ar = ar + ham_k%kinpw_k(displ+jpw)*(cg(1,shift+ipw)*cg(1,shift+ipw)+cg(2,shift+ipw)*cg(2,shift+ipw))
418 : end do
419 : end if
420 115400 : kin = kin + wk*occ0(iband)*ar
421 : end if
422 : end do
423 :
424 27260 : end subroutine rttddft_calc_kin
425 : !!***
426 :
427 : !!****f* m_rttddft_properties/rttddft_calc_enl
428 : !!
429 : !! NAME
430 : !! rttddft_calc_enl
431 : !!
432 : !! FUNCTION
433 : !! Computes the NL part of energy in NC case
434 : !!
435 : !! INPUTS
436 : !! cg <real(2,npw*nspinor*nband)> = the wavefunction coefficients
437 : !! ham_k <gs_hamiltonian_type> = hamiltonian at point k
438 : !! nband <integer> = number of bands
439 : !! npw <integer> = number of plane waves
440 : !! nspinor <integer> = dimension of spinors
441 : !! mpi_enreg <MPI_type> = MPI-parallelisation information
442 : !!
443 : !! OUTPUT
444 : !! enl <real(nband)> = the non local part of the energy in NC case
445 : !!
446 : !! SIDE EFFECTS
447 : !!
448 : !! SOURCE
449 320 : subroutine rttddft_calc_enl(cg,enl,ham_k,nband,npw,nspinor,mpi_enreg)
450 :
451 : !Arguments ------------------------------------
452 : !scalars
453 : integer, intent(in) :: nband
454 : integer, intent(in) :: npw
455 : integer, intent(in) :: nspinor
456 : type(gs_hamiltonian_type), intent(in) :: ham_k
457 : type(MPI_type), intent(in) :: mpi_enreg
458 : !arrays
459 : real(dp), intent(inout) :: cg(2,npw*nspinor*nband)
460 : real(dp), intent(out) :: enl(nband)
461 :
462 : !Local variables-------------------------------
463 : !scalars
464 : integer, parameter :: choice=1
465 : integer, parameter :: cpopt=-1
466 : integer, parameter :: paw_opt=0
467 : integer, parameter :: signs=1
468 : integer, parameter :: tim_getghc = 5
469 : !arrays
470 320 : type(pawcprj_type), allocatable :: cprj_dummy(:,:)
471 320 : real(dp), allocatable :: eig_dummy(:)
472 320 : real(dp), allocatable :: gvnlxc_dummy(:,:)
473 320 : real(dp), allocatable :: gsc_dummy(:,:)
474 : ! ***********************************************************************
475 :
476 640 : ABI_MALLOC(cprj_dummy,(ham_k%natom,0))
477 320 : ABI_MALLOC(gsc_dummy,(0,0))
478 320 : ABI_MALLOC(gvnlxc_dummy, (0, 0))
479 960 : ABI_MALLOC(eig_dummy,(nband))
480 : call nonlop(choice,cpopt,cprj_dummy,enl,ham_k,0,eig_dummy,mpi_enreg,nband, &
481 320 : & 1,paw_opt,signs,gsc_dummy,tim_getghc,cg,gvnlxc_dummy)
482 320 : ABI_FREE(cprj_dummy)
483 320 : ABI_FREE(gsc_dummy)
484 320 : ABI_FREE(eig_dummy)
485 320 : ABI_FREE(gvnlxc_dummy)
486 :
487 320 : end subroutine rttddft_calc_enl
488 : !!***
489 :
490 : !!****f* m_rttddft_properties/rttddft_calc_occ
491 : !!
492 : !! NAME
493 : !! rttddft_calc_occ
494 : !!
495 : !! FUNCTION
496 : !! Computes occupations at time t from cg(t), cg0 and occ0
497 : !! In NC:
498 : !! f_{n,k}(t) = \sum_{m} f_{m,k}(0) <\psi_m(0)|\psi_n(t)>
499 : !! In PAW:
500 : !! f_{n,k}(t) = \sum_{m} f_{m,k}(0) <\psi_m(0)|S|\psi_n(t)>
501 : !! = \sum_{m} f_{m,k}(0) [ <\psi_m(0)|\psi_n(t)> +
502 : !! \sum_{i} <\psi_m(0)|p_{i}>\sum_jS_{i,j}<p_{j}|\psi_n(t)> ]
503 : !!
504 : !! INPUTS
505 : !! cg <real(2,npw*nspinor*nband)> = the wavefunction coefficients
506 : !! cg0 <real(2,npw*nspinor*nband)> = the initial wavefunction coefficients
507 : !! dtset <type(dataset_type)> = all input variables for this dataset
508 : !! ham_k <gs_hamiltonian_type> = hamiltonian at point k
509 : !! ikpt <integer> = indice of the considered k-point
510 : !! ibg <integer> = indice of the considered k-point for cprj
511 : !! isppol <integer> = indice of the considered spin-polarization
512 : !! mpi_enreg <MPI_type> = MPI-parallelisation information
513 : !! nband_k <integer> = number of bands
514 : !! npw_k <integer> = number of plane waves
515 : !! nspinor <integer> = dimension of spinors
516 : !! occ0 <real(nband)> = initial occupations
517 : !! tdks <type(tdks_type)> = Main RT-TDDFT object
518 : !!
519 : !! OUTPUT
520 : !! occ <real(nband)> = the occupations at time t
521 : !!
522 : !! SOURCE
523 27260 : subroutine rttddft_calc_occ(cg,cg0,dtset,ham_k,ikpt,ibg,isppol,mpi_enreg,nband_k,npw_k,nspinor,occ,occ0,tdks)
524 :
525 : !Arguments ------------------------------------
526 : !scalars
527 : integer, intent(in) :: ikpt
528 : integer, intent(in) :: ibg
529 : integer, intent(in) :: isppol
530 : integer, intent(in) :: nband_k
531 : integer, intent(in) :: npw_k
532 : integer, intent(in) :: nspinor
533 : type(dataset_type), intent(in) :: dtset
534 : type(gs_hamiltonian_type), intent(in) :: ham_k
535 : type(MPI_type), intent(in) :: mpi_enreg
536 : type(tdks_type), target, intent(in) :: tdks
537 : !arrays
538 : real(dp), intent(inout) :: cg(2,npw_k*nspinor*nband_k)
539 : real(dp), intent(in) :: cg0(2,npw_k*nspinor*nband_k)
540 : real(dp), intent(out) :: occ(nband_k)
541 : real(dp), intent(in) :: occ0(nband_k)
542 :
543 : !Local variables-------------------------------
544 : !scalars
545 : integer :: iband, jband, ierr, ind
546 : integer :: nband_cprj_k
547 : integer :: natom
548 : integer :: shift
549 : !parameters for nonlop
550 : integer, parameter :: cpopt = 2
551 : integer, parameter :: choice = 1
552 : integer, parameter :: idir = 1
553 : integer, parameter :: nnlout = 1
554 : integer, parameter :: paw_opt = 3
555 : integer, parameter :: signs = 1
556 : integer, parameter :: tim_nonlop = 16
557 : logical :: cprj_paral_band
558 : !arrays
559 54520 : real(dp) :: csc(2*nband_k)
560 27260 : real(dp),allocatable :: gsc(:,:),gvnlxc(:,:)
561 27260 : real(dp),allocatable :: enlout(:),enlout_im(:)
562 27260 : type(pawcprj_type), pointer :: cprj(:,:), cprj_k(:,:)
563 27260 : type(pawcprj_type), pointer :: cprj0(:,:), cprj0_k(:,:)
564 :
565 : ! ***********************************************************************
566 :
567 : !Prepare cprj in PAW case
568 27260 : cprj_paral_band=.false.
569 27260 : if (ham_k%usepaw == 1) then
570 : ! Determine if cprj datastructure is distributed over bands
571 26940 : cprj_paral_band=(tdks%mband_cprj<dtset%mband)
572 26940 : nband_cprj_k=nband_k; if (cprj_paral_band) nband_cprj_k=mpi_enreg%bandpp
573 26940 : natom = dtset%natom
574 : ! Extract the right cprj for this k-point
575 26940 : if (dtset%mkmem*dtset%nsppol/=1) then
576 707280 : ABI_MALLOC(cprj_k,(natom,nspinor*nband_cprj_k))
577 680340 : ABI_MALLOC(cprj0_k,(natom,nspinor*nband_cprj_k))
578 26940 : call pawcprj_alloc(cprj_k,0,tdks%dimcprj)
579 26940 : call pawcprj_alloc(cprj0_k,0,tdks%dimcprj)
580 : call pawcprj_get(tdks%atindx1,cprj_k,tdks%cprj,natom,1,ibg,ikpt,0,isppol, &
581 : & tdks%mband_cprj,dtset%mkmem,natom,nband_cprj_k,nband_cprj_k, &
582 : & nspinor,dtset%nsppol,tdks%unpaw,mpicomm=mpi_enreg%comm_kpt, &
583 26940 : & proc_distrb=mpi_enreg%proc_distrb)
584 : call pawcprj_get(tdks%atindx1,cprj0_k,tdks%cprj0,natom,1,ibg,ikpt,0,isppol, &
585 : & tdks%mband_cprj,dtset%mkmem,natom,nband_cprj_k,nband_cprj_k, &
586 : & nspinor,dtset%nsppol,tdks%unpaw,mpicomm=mpi_enreg%comm_kpt, &
587 26940 : & proc_distrb=mpi_enreg%proc_distrb)
588 : else
589 0 : cprj_k => tdks%cprj
590 0 : cprj0_k => tdks%cprj0
591 : end if
592 : ! If cprj are distributed over bands, gather them (because we need to mix bands)
593 26940 : if (cprj_paral_band) then
594 179200 : ABI_MALLOC(cprj,(natom,nspinor*nband_k))
595 172800 : ABI_MALLOC(cprj0,(natom,nspinor*nband_k))
596 6400 : call pawcprj_alloc(cprj,0,tdks%dimcprj)
597 6400 : call pawcprj_alloc(cprj0,0,tdks%dimcprj)
598 : call pawcprj_mpi_allgather(cprj_k,cprj,natom,nspinor*nband_cprj_k,mpi_enreg%bandpp, &
599 : & tdks%dimcprj,0,mpi_enreg%nproc_band,mpi_enreg%comm_band, &
600 6400 : & ierr,rank_ordered=.false.)
601 : call pawcprj_mpi_allgather(cprj0_k,cprj0,natom,nspinor*nband_cprj_k,mpi_enreg%bandpp, &
602 : & tdks%dimcprj,0,mpi_enreg%nproc_band,mpi_enreg%comm_band, &
603 6400 : & ierr,rank_ordered=.false.)
604 : else
605 20540 : cprj => cprj_k
606 20540 : cprj0 => cprj0_k
607 : end if
608 :
609 : !allocate necessary arrays for nonlop
610 26940 : ABI_MALLOC(gsc,(0,0))
611 26940 : ABI_MALLOC(gvnlxc,(0,0))
612 80820 : ABI_MALLOC(enlout,(nband_k))
613 53880 : ABI_MALLOC(enlout_im,(nband_k))
614 : end if
615 :
616 484540 : csc = zero
617 255900 : do iband = 1, nband_k
618 : !* 1 - Compute csc = <cg0|cg>
619 228640 : shift = npw_k*nspinor*(iband-1)
620 : call zgemv('C',npw_k*nspinor,nband_k,cone,cg0,npw_k*nspinor, &
621 228640 : & cg(:,shift+1:shift+npw_k*nspinor),1,czero,csc,1)
622 : !If band parallel then reduce csc
623 228640 : if (mpi_enreg%nproc_band > 1) then
624 55040 : call xmpi_sum(csc,mpi_enreg%comm_bandfft,ierr)
625 : end if
626 : !* 2 - If PAW, add the additional term \sum_i cprj_i \sum_j S_{i,j} cprj_j
627 228640 : if (ham_k%usepaw == 1) then
628 : call nonlop(choice,cpopt,cprj(:,iband:iband+(nspinor-1)),enlout, &
629 : & ham_k,idir,(/zero/),mpi_enreg,1,nnlout,paw_opt,signs, &
630 : & gsc,tim_nonlop,cg(:,shift+1:shift+npw_k*nspinor),gvnlxc, &
631 220960 : & cprjin_left=cprj0,enlout_im=enlout_im,ndat_left=nband_k)
632 2217120 : do jband = 1, nband_k
633 1996160 : csc(2*jband-1) = csc(2*jband-1) + enlout(jband)
634 2217120 : csc(2*jband) = csc(2*jband) + enlout_im(jband)
635 : end do
636 : end if
637 : !* 3 - Calc occupations from csc and occ0
638 2436380 : do jband = 1, nband_k
639 2180480 : ind = 2*jband
640 2409120 : occ(iband) = occ(iband) + occ0(jband)*(csc(ind-1)**2+csc(ind)**2)
641 : end do
642 : end do
643 :
644 27260 : if (ham_k%usepaw == 1) then
645 26940 : ABI_FREE(gsc)
646 26940 : ABI_FREE(gvnlxc)
647 26940 : ABI_FREE(enlout)
648 26940 : ABI_FREE(enlout_im)
649 26940 : if (cprj_paral_band) then
650 6400 : call pawcprj_free(cprj)
651 6400 : ABI_FREE(cprj)
652 6400 : call pawcprj_free(cprj0)
653 6400 : ABI_FREE(cprj0)
654 : end if
655 26940 : if (dtset%mkmem*dtset%nsppol/=1) then
656 26940 : call pawcprj_free(cprj_k)
657 26940 : ABI_FREE(cprj_k)
658 26940 : call pawcprj_free(cprj0_k)
659 26940 : ABI_FREE(cprj0_k)
660 : end if
661 : end if
662 :
663 54520 : end subroutine rttddft_calc_occ
664 : !!***
665 :
666 : !!****f* m_rttddft_properties/rttddft_calc_ent
667 : !!
668 : !! NAME
669 : !! rttddft_calc_ent
670 : !!
671 : !! FUNCTION
672 : !! Computes electronic entropy from occupation numbers f_{nk}
673 : !! S = -2 \sum_k \sum_n w(k) [ f_{nk}*ln(f_{nk}) + (1-f_{nk})*ln(1-f_{nk}) ]
674 : !!
675 : !! INPUTS
676 : !! dtset <type(dataset_type)> = all input variables for this dataset
677 : !! occ <real(nband*nkpt*nsspol)> = occupation numbers at time t
678 : !!
679 : !! OUTPUT
680 : !! entropy <real>
681 : !!
682 : !! SOURCE
683 1725 : subroutine rttddft_calc_ent(entropy,dtset,occ)
684 :
685 : !Arguments ------------------------------------
686 : !scalars
687 : real(dp), intent(out) :: entropy
688 : type(dataset_type), intent(in) :: dtset
689 : !arrays
690 : real(dp), intent(in) :: occ(:)
691 :
692 : !Local variables-------------------------------
693 : !scalars
694 : integer :: band_index
695 : integer :: iband, isppol, ikpt
696 : integer :: nband_k
697 : real(dp) :: fnk
698 :
699 : ! ***********************************************************************
700 :
701 1725 : entropy = zero
702 1725 : band_index=0
703 :
704 3450 : do isppol = 1, dtset%nsppol
705 50390 : do ikpt = 1, dtset%nkpt
706 46940 : nband_k = dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
707 440700 : do iband = 1, nband_k
708 393760 : fnk = occ(iband+band_index)*0.5_dp
709 440700 : if( fnk>tol16 .and. (one-fnk)>tol16 ) then
710 352957 : entropy = entropy - two*dtset%wtk(ikpt)*(fnk*log(fnk)+(one-fnk)*log(one-fnk))
711 : end if
712 : end do
713 48665 : band_index = band_index + nband_k
714 : end do
715 : end do
716 :
717 1725 : end subroutine rttddft_calc_ent
718 : !!***
719 :
720 : !!****f* m_rttddft_properties/rttddft_calc_current
721 : !!
722 : !! NAME
723 : !! rttddft_calc_current
724 : !!
725 : !! FUNCTION
726 : !! Computes macroscopic current density
727 : !! In NC:
728 : !! J(t) = -1/Omega Im[\sum_{n,k} f_{nk}(0) <\tilde{psi}_{nk}|\nabla|\tilde{psi}_{nk}>] - A N_v/Omega
729 : !! + Gauge-dependent terms?
730 : !! In PAW:
731 : !! J(t) = -1/Omega Im[ \sum_{n,k} f_{nk}(0) <\tilde{psi}_{nk}|\nabla|\tilde{psi}_{nk}> +
732 : !! \sum_{aij} \rho_{aij} <\phi_{aij}|\nabla|\phi_{aij}> -
733 : !! \sum_{aij} \rho_{aij} <\tilde{\phi}_{aij}|\nabla|\tilde{\phi}_{aij}> ]
734 : !! - A/Omega \int \tilde{n}(r,t)dr
735 : !! INPUTS
736 : !! tdks <type(tdks_type)> = Main RT-TDDFT object
737 : !! dtset <type(dataset_type)> = all input variables for this dataset
738 : !! mpi_enreg <MPI_type> = MPI-parallelisation information
739 : !!
740 : !! OUTPUT
741 : !! current = the macroscopic current
742 : !!
743 : !! NOTES
744 : !! Im[<\tilde{psi}_{nk}|\nabla|\tilde{psi}_{nk}>] = <\tilde{psi}_{nk}|-i\nabla|\tilde{psi}_{nk}>
745 : !!
746 : !! SOURCE
747 1121 : subroutine rttddft_calc_current(tdks, dtset, dtfil, psps, mpi_enreg)
748 :
749 : !Arguments ------------------------------------
750 : !scalars
751 : type(tdks_type), target, intent(inout) :: tdks
752 : type(dataset_type), intent(in) :: dtset
753 : type(datafiles_type), intent(in) :: dtfil
754 : type(pseudopotential_type), intent(in) :: psps
755 : type(MPI_type), intent(in) :: mpi_enreg
756 :
757 : !Local variables-------------------------------
758 : !scalars
759 : integer :: iband, ikpt, isppol, bdtot_index, nband_k
760 : real(dp) :: current_k(3)
761 : !arrays
762 2242 : real(dp), target :: psinablapsi(2,3,dtset%mband,dtset%nkpt)
763 :
764 : ! ***********************************************************************
765 :
766 1121 : if (psps%usepaw==1) then
767 : ! 1 - Computes <\psi_{nk}|v|\psi_{nk}> = <\psi_{nk}|-i\nabla|\psi_{nk}> = Im[<\tilde{psi}_{nk}|\nabla|\tilde{psi}_{nk}>]
768 : call optics_paw(tdks%atindx1,tdks%cg,tdks%cprj,tdks%dimcprj,dtfil,dtset,tdks%eigen,tdks%gprimd,tdks%hdr, &
769 : & tdks%kg,dtset%mband,tdks%mcg,tdks%mcprj,dtset%mkmem,mpi_enreg,psps%mpsang,dtset%mpw, &
770 : & dtset%natom,dtset%nkpt,tdks%npwarr,dtset%nsppol,tdks%pawang,tdks%pawrad,tdks%pawrhoij, &
771 1121 : & tdks%pawtab,dtset%znucl,psinablapsi)
772 :
773 : ! 2 - Sum over bands and k-points
774 5605 : tdks%current = 0.0_dp
775 1121 : bdtot_index=0
776 : !Loop over spins
777 2242 : do isppol=1, dtset%nsppol
778 : !Loop over kpoints
779 36993 : do ikpt = 1, dtset%nkpt
780 35872 : current_k = 0.0_dp
781 35872 : nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
782 322848 : do iband = 1, nband_k
783 286976 : current_k(1) = current_k(1) + tdks%occ0(bdtot_index+iband)*psinablapsi(1,1,iband,ikpt)
784 286976 : current_k(2) = current_k(2) + tdks%occ0(bdtot_index+iband)*psinablapsi(1,2,iband,ikpt)
785 322848 : current_k(3) = current_k(3) + tdks%occ0(bdtot_index+iband)*psinablapsi(1,3,iband,ikpt)
786 : end do
787 143488 : tdks%current(:,isppol) = tdks%current(:,isppol) + dtset%wtk(ikpt)*current_k(:)
788 36993 : bdtot_index = bdtot_index + nband_k
789 : end do !nkpt
790 :
791 : ! 3 - Add last contribution from vector potential times integral of the density
792 : tdks%current(:,isppol) = tdks%current(:,isppol) + &
793 16092947 : & tdks%tdef%vecpot(:)*sum(tdks%rhor(:,isppol)-tdks%nhat(:,isppol))*tdks%ucvol/tdks%nfftf
794 :
795 5605 : tdks%current(:,isppol) = -tdks%current(:,isppol)/tdks%ucvol
796 :
797 : end do !nsspol
798 :
799 : end if
800 :
801 1121 : end subroutine rttddft_calc_current
802 : !!***
803 :
804 : end module m_rttddft_properties
805 : !!***
|
