Line data Source code
1 : !!****f* ABINIT/m_slicewf_cprj
2 : !! NAME
3 : !! m_slicewf_cprj
4 : !!
5 : !! FUNCTION
6 : !! This module contains a routine updating the whole wave functions at a given k-point,
7 : !! using the Spectrum Slicing filtering method (2021 implementation using xG abstraction layer)
8 : !! for a given spin-polarization, from a fixed hamiltonian
9 : !! but might also simply compute eigenvectors and eigenvalues at this k point.
10 : !! it will also update the matrix elements of the hamiltonian.
11 : !!
12 : !! COPYRIGHT
13 : !! Copyright (C) 2023-2025 ABINIT group (LB,IML)
14 : !! This file is distributed under the terms of the
15 : !! gnu general public license, see ~abinit/COPYING
16 : !! or http://www.gnu.org/copyleft/gpl.txt .
17 : !! for the initials of contributors, see ~abinit/doc/developers/contributors.txt .
18 : !!
19 : !! PARENTS
20 : !! vtowfk
21 : !!
22 : !! CHILDREN
23 : !!
24 : !! SOURCE
25 :
26 : #if defined HAVE_CONFIG_H
27 : #include "config.h"
28 : #endif
29 :
30 : #include "abi_common.h"
31 :
32 : module m_slicewf_cprj
33 :
34 : use defs_abitypes
35 : use defs_basis
36 : use m_abicore
37 : use m_errors
38 : use m_fstrings
39 : use m_time
40 : use m_xg
41 : use m_xg_nonlop
42 : use m_slice_cprj
43 : use m_invovl
44 :
45 : use m_dtset, only : dataset_type
46 : use m_hamiltonian, only : gs_hamiltonian_type
47 : use m_pawcprj, only : pawcprj_type
48 : use m_getghc, only : multithreaded_getghc
49 :
50 : use m_xg
51 : use m_xgTransposer
52 :
53 : use iso_c_binding, only: c_associated,c_loc,c_ptr,c_f_pointer
54 :
55 : use m_xmpi
56 : use m_xomp
57 : #ifdef HAVE_OPENMP
58 : use omp_lib
59 : #endif
60 :
61 : implicit none
62 :
63 : private
64 :
65 : integer, parameter :: l_tim_getghc=5
66 :
67 : ! For use in getghc_gsc1
68 : integer, save :: l_prtvol
69 : type(mpi_type),pointer,save :: l_mpi_enreg
70 : type(gs_hamiltonian_type),pointer,save :: l_gs_hamk
71 :
72 : public :: slicewf_cprj
73 :
74 : contains
75 : !!***
76 :
77 : !!****f* m_slicewf/slicewf_cprj
78 : !! NAME
79 : !! slicewf_cprj
80 : !!
81 : !! FUNCTION
82 : !! This routine updates the whole wave functions set at a given k-point,
83 : !! using the Spectrum Slicing method (2021 version using xG abstraction layer)
84 : !!
85 : !! INPUTS
86 : !! dtset= input variables for this dataset
87 : !! mpi_enreg= MPI-parallelisation information
88 : !! nband= number of bands at this k point
89 : !! npw= number of plane waves at this k point
90 : !! nspinor= number of spinorial components of the wavefunctions
91 : !! prtvol= control print volume and debugging
92 : !!
93 : !! OUTPUT
94 : !! eig(nband)= eigenvalues (hartree) for all bands
95 : !! enl_out(nband)= contribution of each band to the nl part of energy
96 : !! resid(nband)= residuals for each band
97 : !!
98 : !! SIDE EFFECTS
99 : !! cg(2,npw*nspinor*nband)= planewave coefficients of wavefunctions
100 : !! gs_hamk <type(gs_hamiltonian_type)>=all data for the hamiltonian at k
101 : !!
102 : !! SOURCE
103 :
104 0 : subroutine slicewf_cprj(cg,dtset,eig,occ,enl_out,gs_hamk,mpi_enreg,&
105 0 : & nband,npw,nspinor,prtvol,resid,xg_nonlop)
106 :
107 : !Arguments ------------------------------------
108 : integer,intent(in) :: nband,npw,prtvol,nspinor
109 : type(gs_hamiltonian_type),target,intent(inout) :: gs_hamk
110 : type(dataset_type) ,intent(in ) :: dtset
111 : type(mpi_type) ,target,intent(in) :: mpi_enreg
112 : real(dp) ,target,intent(inout) :: cg(2,nspinor*nband*npw)
113 : real(dp) ,target,intent( out) :: resid(nband)
114 : real(dp) ,intent( out) :: enl_out(nband)
115 : real(dp) ,target,intent( out) :: eig(nband)
116 : real(dp) ,target,intent(in ) :: occ(nband)
117 : type(xg_nonlop_t) , intent(in ) :: xg_nonlop
118 :
119 : !Local variables-------------------------------
120 :
121 : type(xgBlock_t) :: xgx0
122 : type(xg_t) :: cprj_xgx0
123 : type(xgBlock_t) :: xgeigen
124 : type(xgBlock_t) :: xgocc
125 : type(xgBlock_t) :: xgresidu
126 : type(xgBlock_t) :: xgenl
127 : type(xgBlock_t) :: xg_kin
128 0 : type(slice_t) :: slice
129 :
130 : logical :: paw
131 :
132 : integer :: space, space_cprj, blockdim, cprjdim, nband_cprj
133 : integer :: me_g0,me_g0_fft
134 :
135 : integer, parameter :: tim_slicewf = 2170
136 : double precision :: tsec(2)
137 :
138 0 : real(dp), allocatable :: kin(:),occ_tmp(:)
139 :
140 : ! *********************************************************************
141 :
142 0 : call timab(tim_slicewf,1,tsec)
143 :
144 0 : paw = gs_hamk%usepaw==1
145 :
146 : ! Set module variables
147 0 : l_prtvol = prtvol
148 0 : l_mpi_enreg => mpi_enreg
149 0 : l_gs_hamk => gs_hamk
150 :
151 0 : cprjdim = xg_nonlop%cprjdim
152 :
153 : !Variables
154 0 : blockdim=mpi_enreg%nproc_band*mpi_enreg%bandpp
155 0 : nband_cprj=nband/mpi_enreg%nproc_band
156 :
157 : !Depends on istwfk
158 0 : if ( gs_hamk%istwf_k > 1 ) then ! Real only
159 : ! SPACE_CR mean that we have complex numbers but no re*im terms only re*re
160 : ! and im*im so that a vector of complex is consider as a long vector of real
161 : ! therefore the number of data is (2*npw*nspinor)*nband
162 : ! This space is completely equivalent to SPACE_R but will correctly set and
163 : ! get the array data into the xgBlock
164 0 : space = SPACE_CR
165 : else ! complex
166 0 : space = SPACE_C
167 : end if
168 0 : space_cprj = xg_nonlop%space_cprj
169 :
170 : !For kinetic part of the Hamiltonian
171 0 : ABI_MALLOC(kin,(gs_hamk%npw_fft_k))
172 0 : call build_kin(kin,gs_hamk%kinpw_k,gs_hamk%npw_fft_k)
173 0 : call xgBlock_map_1d(xg_kin,kin,SPACE_R,gs_hamk%npw_fft_k)
174 :
175 : ! Local variables for slice
176 0 : me_g0 = -1
177 0 : me_g0_fft = -1
178 0 : if (space==SPACE_CR) then
179 0 : me_g0 = 0
180 0 : me_g0_fft = 0
181 0 : if (gs_hamk%istwf_k == 2) then
182 0 : if (l_mpi_enreg%me_g0 == 1) me_g0 = 1
183 0 : if (l_mpi_enreg%me_g0_fft == 1) me_g0_fft = 1
184 : end if
185 : end if
186 0 : call xgBlock_map(xgx0,cg,space,npw*nspinor,nband,l_mpi_enreg%comm_band,me_g0=me_g0)
187 :
188 0 : call xgBlock_map_1d(xgeigen,eig,SPACE_R,nband)
189 :
190 0 : call xgBlock_map_1d(xgresidu,resid,SPACE_R,nband)
191 :
192 0 : call xgBlock_map_1d(xgenl,enl_out,SPACE_R,nband)
193 :
194 : ! Occupancies in chebyshev are used for convergence criteria only
195 0 : if (dtset%nbdbuf==-101.and.nspinor==1.and.dtset%nsppol==1) then
196 0 : ABI_MALLOC(occ_tmp,(nband))
197 0 : occ_tmp(:) = half*occ(:)
198 0 : call xgBlock_map_1d(xgocc,occ_tmp,SPACE_R,nband,gpu_option=dtset%gpu_option)
199 : else
200 0 : call xgBlock_map_1d(xgocc,occ,SPACE_R,nband,gpu_option=dtset%gpu_option)
201 : end if
202 :
203 : !call xg_cprj_copy(cprj_cwavef_bands,cprj_contiguous,space_cprj,nband_cprj,cprj_xgx0,&
204 : ! & xg_nonlop,l_mpi_enreg%comm_band,CPRJ_ALLOC)
205 0 : call xg_init(cprj_xgx0,space_cprj,xg_nonlop%cprjdim,nband_cprj*nspinor,comm=l_mpi_enreg%comm_band)
206 :
207 : call slice_init(slice,nband,npw*nspinor,cprjdim,dtset%tolwfr_diago,dtset%ecut, &
208 : & mpi_enreg%bandpp, dtset%nline, dtset%nbdbuf, space,space_cprj,1, &
209 : & l_mpi_enreg%comm_band,me_g0,paw,&
210 : & dtset%nslice,0.d0,dtset%paral_slice,0,&
211 0 : & xg_nonlop,me_g0_fft)
212 :
213 : ! Run slice
214 0 : call slice_run_cprj(slice,xgx0,cprj_xgx0%self,xg_getghc,xg_kin,xgeigen,xgocc,xgresidu,xgenl,nspinor)
215 :
216 0 : if (allocated(occ_tmp)) then
217 0 : ABI_FREE(occ_tmp)
218 : end if
219 0 : ABI_FREE(kin)
220 :
221 : ! call xg_cprj_copy(cprj_cwavef_bands,cprj_contiguous,space_cprj,nband_cprj,cprj_xgx0,&
222 : ! & xg_nonlop,l_mpi_enreg%comm_band,CPRJ_FREE)
223 0 : call xg_free(cprj_xgx0)
224 :
225 : ! Free slice
226 0 : call slice_free(slice)
227 :
228 0 : call timab(tim_slicewf,2,tsec)
229 :
230 : DBG_EXIT("COLL")
231 :
232 0 : end subroutine slicewf_cprj
233 : !!***
234 :
235 : !!****f* m_slice/xg_getghc
236 : !! NAME
237 : !! xg_getghc
238 : !!
239 : !! FUNCTION
240 : !! This routine computes H|C> and possibly S|C> for a given wave function C.
241 : !! It acts as a driver for getghc, taken into account parallelism, multithreading, etc.
242 : !!
243 : !! SIDE EFFECTS
244 : !! X <type(xgBlock_t)>= memory block containing |C>
245 : !! AX <type(xgBlock_t)>= memory block containing H|C>
246 : !!
247 : !! PARENTS
248 : !!
249 : !! CHILDREN
250 : !! xgBlock_getSize,xgBlock_reverseMap,xgBlock_scale,xgBlock_copy
251 : !! multithreaded_getghc
252 : !!
253 : !! SOURCE
254 : !
255 0 : subroutine xg_getghc(X,AX)
256 :
257 : use iso_c_binding
258 :
259 : !Arguments ------------------------------------
260 : type(xgBlock_t), intent(inout) :: X
261 : type(xgBlock_t), intent(inout) :: AX
262 :
263 : !Local variables-------------------------------
264 : !scalars
265 : integer :: blockdim
266 : integer :: spacedim
267 : integer,parameter :: sij_opt=0,cpopt=-1,type_calc=1 ! Compute local part only
268 : ! integer :: iatom,iband,ispinor,cprj_index,cprj_rows,cprj_cols,ncpgr,nlmn
269 : real(dp) :: eval
270 0 : type(pawcprj_type) :: cprj_dum(l_gs_hamk%natom,1)
271 : !arrays
272 0 : real(dp), pointer :: cg(:,:)
273 0 : real(dp), pointer :: ghc(:,:)
274 : real(dp) :: gsc(1,1),gvnlxc(1,1)
275 :
276 : ! *********************************************************************
277 :
278 0 : call xgBlock_getSize(X,spacedim,blockdim)
279 0 : call xgBlock_check(X,AX)
280 :
281 0 : call xgBlock_reverseMap(X,cg,rows=1,cols=spacedim*blockdim)
282 0 : call xgBlock_reverseMap(AX,ghc,rows=1,cols=spacedim*blockdim)
283 :
284 : ! Apply only local part of the Hamiltonian
285 : call multithreaded_getghc(cpopt,cg,cprj_dum,ghc,gsc,&
286 0 : l_gs_hamk,gvnlxc,eval,l_mpi_enreg,blockdim,l_prtvol,sij_opt,l_tim_getghc,type_calc)
287 :
288 0 : end subroutine xg_getghc
289 : !!***
290 :
291 0 : subroutine build_kin(kin,kinpw,npw)
292 :
293 : integer,intent(in) :: npw
294 : real(dp),intent(in) :: kinpw(:)
295 : real(dp),intent(out) :: kin(:)
296 :
297 : integer :: ipw
298 :
299 : !$omp parallel do schedule(static), shared(kin,kinpw)
300 0 : do ipw=1,npw
301 0 : if(kinpw(ipw)>huge(0.0_dp)*1.d-11) then
302 0 : kin(ipw)=0.d0
303 : else
304 0 : kin(ipw) = kinpw(ipw)
305 : end if
306 : end do
307 :
308 0 : end subroutine build_kin
309 :
310 : end module m_slicewf_cprj
311 : !!***
|
