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