Line data Source code
1 : !!****m* ABINIT/m_rttddft_exponential
2 : !! NAME
3 : !! m_rttddft_exponential
4 : !!
5 : !! FUNCTION
6 : !! Contains subroutines to compute the exponential
7 : !! part of the propagator using various approximations
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_exponential
24 :
25 : use defs_basis
26 : use defs_abitypes, only: MPI_type
27 :
28 : use m_bandfft_kpt, only: bandfft_kpt, &
29 : & bandfft_kpt_get_ikpt
30 : use m_dtset, only: dataset_type
31 : use m_getghc, only: multithreaded_getghc
32 : use m_hamiltonian, only: gs_hamiltonian_type
33 : use m_invovl, only: apply_invovl
34 : use m_pawcprj, only: pawcprj_type, pawcprj_alloc, pawcprj_free
35 : use m_prep_kgb, only: prep_getghc, prep_index_wavef_bandpp
36 : use m_profiling_abi, only: abimem_record
37 : use m_rttddft_properties, only: rttddft_calc_eig, rttddft_calc_enl
38 : use m_xmpi, only: xmpi_alltoallv
39 :
40 : implicit none
41 :
42 : private
43 : !!***
44 :
45 : public :: rttddft_exp_taylor
46 : !!***
47 :
48 : contains
49 : !!***
50 :
51 : !!****f* m_rttddft_exponential/rttddft_exp_taylor
52 : !!
53 : !! NAME
54 : !! rttddft_propagator_exp
55 : !!
56 : !! FUNCTION
57 : !! Applies the propagator exp(-i*dt*H(*S^-1)) on the cg coeffs
58 : !! approximating the exponential using a Taylor exapansion
59 : !!
60 : !! INPUTS
61 : !! cg <real(npw*nspinor*nband)> = the wavefunction coefficients
62 : !! dtset <type(dataset_type)>=all input variables for this dataset
63 : !! ham_k <type(gs_hamiltonian_type)> = the Hamiltonian H
64 : !! mpi_enreg <MPI_type> = MPI-parallelisation information
65 : !! nband_k <integer> = number of bands
66 : !! npw_k <integer> = number of plane waves
67 : !! nspinor <integer> = dimension of spinors
68 : !!
69 : !! OUTPUT
70 : !! cg <real(npw*nspinor*nband)> = the new cg after application
71 : !! of the exponential propagator
72 : !! enl <real(bandpp)> = non local contribution to the energy in the
73 : !! NC case - optional
74 : !! eig <real(bandpp)> = eigenvalues - optional
75 : !!
76 : !! SOURCE
77 27380 : subroutine rttddft_exp_taylor(cg,dtset,ham_k,mpi_enreg,nband_k,npw_k,nspinor,enl,eig)
78 :
79 : !Arguments ------------------------------------
80 : !scalars
81 : integer, intent(in) :: nband_k
82 : integer, intent(in) :: npw_k
83 : integer, intent(in) :: nspinor
84 : type(dataset_type), intent(in) :: dtset
85 : type(gs_hamiltonian_type), intent(inout) :: ham_k
86 : type(MPI_type), intent(inout) :: mpi_enreg
87 : !arrays
88 : real(dp), target, intent(inout) :: cg(2,npw_k*nband_k*nspinor)
89 : real(dp), intent(out), optional :: enl(:)
90 : real(dp), pointer, intent(inout), optional :: eig(:)
91 :
92 : !Local variables-------------------------------
93 : !scalars
94 : integer :: cpopt
95 : integer :: iorder
96 : integer :: nband_t
97 : integer :: npw_t
98 : integer :: sij_opt
99 : integer :: tim_getghc = 5
100 : logical :: l_paw, l_eig, l_enl
101 : real(dp) :: dt
102 : !arrays
103 27380 : integer, allocatable :: index_wavef_band(:)
104 27380 : type(pawcprj_type), allocatable :: cwaveprj(:,:)
105 : real(dp), pointer :: cg_t(:,:) => null()
106 27380 : real(dp), target, allocatable :: cg_trans(:,:)
107 27380 : real(dp), allocatable :: ghc(:,:)
108 27380 : real(dp), allocatable :: gvnlxc_dummy(:,:)
109 27380 : real(dp), allocatable :: gsm1hc(:,:)
110 27380 : real(dp), allocatable :: gsc(:,:)
111 27380 : real(dp), allocatable :: tmp(:,:)
112 :
113 : ! ***********************************************************************
114 :
115 : !Check what properties we need to compute
116 27380 : l_eig = .false.; l_enl = .false.
117 27380 : if (present(eig)) then
118 27260 : if (associated(eig)) l_eig = .true.
119 : end if
120 27380 : if (present(enl)) l_enl = (size(enl)>0)
121 :
122 : !Transpose if paral_kgb
123 27380 : if (dtset%paral_kgb == 1 .and. mpi_enreg%nproc_band>1) then
124 6560 : call paral_kgb_transpose(cg,cg_trans,mpi_enreg,nband_t,npw_t,nspinor,1,index_wavef_band)
125 6560 : cg_t => cg_trans
126 : else
127 20820 : npw_t = npw_k
128 20820 : nband_t = nband_k
129 20820 : cg_t => cg
130 : end if
131 :
132 27380 : l_paw = (ham_k%usepaw == 1)
133 27380 : if(l_paw) then
134 724400 : ABI_MALLOC(cwaveprj, (ham_k%natom,nspinor*nband_t))
135 27020 : call pawcprj_alloc(cwaveprj,0,ham_k%dimcprj)
136 : else
137 360 : ABI_MALLOC(cwaveprj,(0,0))
138 : end if
139 27380 : cpopt = 0
140 :
141 27380 : dt = dtset%dtele !electronic timestep
142 :
143 82140 : ABI_MALLOC(ghc, (2, npw_t*nspinor*nband_t))
144 54760 : ABI_MALLOC(gsc, (2, npw_t*nspinor*nband_t))
145 54760 : ABI_MALLOC(gsm1hc,(2, npw_t*nspinor*nband_t))
146 27380 : ABI_MALLOC(gvnlxc_dummy, (0, 0))
147 :
148 : !*** Taylor expansion ***
149 54760 : ABI_MALLOC(tmp,(2, npw_t*nspinor*nband_t))
150 158948300 : tmp(:,:) = cg_t(:,:)
151 :
152 136900 : do iorder = 1, dtset%td_exp_order
153 :
154 : !** Apply Hamiltonian
155 109520 : sij_opt = 0
156 109520 : if (iorder == 1 .and. l_eig .and. l_paw) sij_opt = 1
157 109520 : if (dtset%paral_kgb /= 1) then
158 : call multithreaded_getghc(cpopt,tmp,cwaveprj,ghc,gsc,ham_k,gvnlxc_dummy,1.0_dp, &
159 83280 : & mpi_enreg,nband_t,dtset%prtvol,sij_opt,tim_getghc,0)
160 : else
161 : call prep_getghc(tmp,ham_k,gvnlxc_dummy,ghc,gsc,1.0_dp,nband_t,mpi_enreg, &
162 26240 : & dtset%prtvol,sij_opt,cpopt,cwaveprj,already_transposed=.true.)
163 : end if
164 :
165 : !** Also apply S^-1 in PAW case
166 109520 : if (l_paw) then
167 108080 : call apply_invovl(ham_k,ghc,gsm1hc,cwaveprj,npw_t,nband_t,mpi_enreg,nspinor,dtset%invovl_blksliced)
168 176904080 : tmp(1,:) = dt*gsm1hc(2,:)/real(iorder,dp)
169 176904080 : tmp(2,:) = -dt*gsm1hc(1,:)/real(iorder,dp)
170 : else
171 35100000 : tmp(1,:) = dt*ghc(2,:)/real(iorder,dp)
172 35100000 : tmp(2,:) = -dt*ghc(1,:)/real(iorder,dp)
173 : end if
174 :
175 : !** Compute properties if requested
176 109520 : if (iorder == 1) then
177 : !eigenvalues
178 27380 : if (l_eig) then
179 0 : if (l_paw) then
180 : call rttddft_calc_eig(cg_t,eig,ghc,ham_k%istwf_k,nband_t,npw_t,nspinor, &
181 0 : & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft,gsc=gsc)
182 : else
183 : call rttddft_calc_eig(cg_t,eig,ghc,ham_k%istwf_k,nband_t,npw_t,nspinor, &
184 0 : & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
185 : end if
186 : end if
187 :
188 : !non local PSP energy contribution in NC case
189 27380 : if (l_enl) then
190 320 : call rttddft_calc_enl(cg_t,enl,ham_k,nband_t,npw_t,nspinor,mpi_enreg)
191 : end if
192 : end if
193 :
194 : !** Update cg
195 635820580 : cg_t(:,:) = cg_t(:,:) + tmp(:,:)
196 :
197 : end do
198 :
199 27380 : ABI_FREE(tmp)
200 27380 : ABI_FREE(ghc)
201 27380 : ABI_FREE(gsc)
202 27380 : ABI_FREE(gsm1hc)
203 27380 : ABI_FREE(gvnlxc_dummy)
204 :
205 27380 : if(l_paw) call pawcprj_free(cwaveprj)
206 444980 : ABI_FREE(cwaveprj)
207 :
208 : !Transpose back if paral_kgb
209 27380 : if (dtset%paral_kgb == 1 .and. mpi_enreg%nproc_band > 1) then
210 6560 : call paral_kgb_transpose(cg,cg_trans,mpi_enreg,nband_t,npw_t,nspinor,-1,index_wavef_band)
211 : end if
212 :
213 27380 : nullify(cg_t)
214 :
215 27380 : end subroutine rttddft_exp_taylor
216 : !!***
217 :
218 : !!****f* m_rttddft_exponential/paral_kgb_transpose
219 : !!
220 : !! NAME
221 : !! paral_kgb_transpose
222 : !!
223 : !! FUNCTION
224 : !! if option = 1: Forward transpose
225 : !! Transpose cg_1 in linalg ((npw/npband),nband) distribution
226 : !! into cg_2 in fft (npw,bandpp) distribution
227 : !!
228 : !! if option = -1: Backward transpose
229 : !! Transpose back cg_2 in fft (npw,bandpp) distribution
230 : !! into cg_1 linakg (npw/npband),nband) distribution
231 : !!
232 : !! INPUTS
233 : !! if option = 1:
234 : !! cg_1 <real((npw/nband)*nspinor*nband)>
235 : !! if option = -1:
236 : !! cg_2 <real(npw*nspinor*bandpp)>
237 : !! nband_t <integer> = number of bands after forward transpose (bandpp)
238 : !! npw_t <integer> = number of pw after forward transpose (npw_k)
239 : !! dtset <type(dataset_type)> = all input variables for this dataset
240 : !! mpi_enreg <MPI_type> = MPI-parallelisation information
241 : !! nspinor <integer> = "dimension of spinors"
242 : !! option <integer> = option for forward or backward transpose
243 : !! index_wavef_band <integer> = order of the bands after transpose
244 : !!
245 : !! OUTPUT
246 : !! if option = 1 :
247 : !! cg_2 <real(npw*nspinor*bandpp)>
248 : !! nband_t <integer> = number of bands after forward transpose (bandpp)
249 : !! npw_t <integer> = number of pw after forward transpose (npw_k)
250 : !! if option = -1:
251 : !! cg_1 <real((npw/nband)*nspinor*nband)>
252 :
253 : !! SIDE EFFECTS
254 : !! if option = 1 :
255 : !! cg_2 has been allocated
256 : !! index_wavef_band has been allocated
257 : !! if option = -1:
258 : !! cg_2 has been deallocated
259 : !! index_wavef_band has been allocated
260 : !!
261 : !! SOURCE
262 13120 : subroutine paral_kgb_transpose(cg_1,cg_2,mpi_enreg,nband_t,npw_t,nspinor,option,index_wavef_band)
263 :
264 : !Arguments ------------------------------------
265 : !scalars
266 : integer, intent(inout) :: nband_t
267 : integer, intent(inout) :: npw_t
268 : integer, intent(in) :: nspinor
269 : integer, intent(in) :: option
270 : type(MPI_type), intent(inout) :: mpi_enreg
271 : !arrays
272 : integer, allocatable, intent(inout) :: index_wavef_band(:)
273 : real(dp), intent(inout) :: cg_1(:,:)
274 : real(dp), allocatable, intent(inout) :: cg_2(:,:)
275 :
276 : !Local variables-------------------------------
277 : !scalars
278 : integer :: bandpp
279 : integer :: ierr
280 : integer :: ikpt_this_proc
281 : !arrays
282 26240 : integer :: recvcountsloc(mpi_enreg%nproc_band)
283 26240 : integer :: rdisplsloc(mpi_enreg%nproc_band)
284 26240 : integer :: sendcountsloc(mpi_enreg%nproc_band)
285 26240 : integer :: sdisplsloc(mpi_enreg%nproc_band)
286 13120 : real(dp), allocatable :: cg_work(:,:)
287 :
288 : ! ***********************************************************************
289 :
290 : !Init useful MPI variables
291 13120 : bandpp = mpi_enreg%bandpp
292 13120 : ikpt_this_proc = bandfft_kpt_get_ikpt()
293 39360 : recvcountsloc = bandfft_kpt(ikpt_this_proc)%recvcounts*2*nspinor*bandpp
294 39360 : rdisplsloc = bandfft_kpt(ikpt_this_proc)%rdispls*2*nspinor*bandpp
295 39360 : sendcountsloc = bandfft_kpt(ikpt_this_proc)%sendcounts*2*nspinor
296 39360 : sdisplsloc = bandfft_kpt(ikpt_this_proc)%sdispls*2*nspinor
297 :
298 : !Forward transpose: cg_1 -> cg_2
299 13120 : if (option == 1) then
300 6560 : nband_t = bandpp
301 6560 : npw_t = bandfft_kpt(ikpt_this_proc)%ndatarecv
302 19680 : ABI_MALLOC(cg_2, (2, npw_t*nspinor*nband_t))
303 13120 : ABI_MALLOC(cg_work, (2, npw_t*nspinor*nband_t))
304 : !Transpose input cg_1 into cg_work
305 : call xmpi_alltoallv(cg_1,sendcountsloc,sdisplsloc,cg_work, &
306 6560 : & recvcountsloc,rdisplsloc,mpi_enreg%comm_band,ierr)
307 : !properly sort array according to bandd after alltoall
308 : call prep_index_wavef_bandpp(mpi_enreg%nproc_band,mpi_enreg%bandpp, &
309 : & nspinor,bandfft_kpt(ikpt_this_proc)%ndatarecv , &
310 : & bandfft_kpt(ikpt_this_proc)%recvcounts, &
311 6560 : & bandfft_kpt(ikpt_this_proc)%rdispls, index_wavef_band)
312 18759680 : cg_2(:,:) = cg_work(:,index_wavef_band)
313 6560 : ABI_FREE(cg_work)
314 : end if
315 :
316 : !Transpose back: cg_2 -> cg_1
317 13120 : if (option == -1) then
318 19680 : ABI_MALLOC(cg_work, (2, npw_t*nspinor*nband_t))
319 18759680 : cg_work(:,index_wavef_band) = cg_2(:,:)
320 : !Transpose cg_work to input cg_1
321 : call xmpi_alltoallv(cg_work,recvcountsloc,rdisplsloc,cg_1, &
322 6560 : & sendcountsloc,sdisplsloc,mpi_enreg%comm_band,ierr)
323 : !Free memory
324 6560 : ABI_FREE(cg_work)
325 6560 : ABI_FREE(cg_2)
326 6560 : ABI_FREE(index_wavef_band)
327 : end if
328 :
329 13120 : end subroutine paral_kgb_transpose
330 : !!***
331 :
332 : end module m_rttddft_exponential
333 : !!***
|
