Line data Source code
1 : !!****f* ABINIT/m_lobpcgwf_cprj
2 : !! NAME
3 : !! m_lobpcgwf_cprj
4 : !!
5 : !! FUNCTION
6 : !! This routine updates the whole wave functions at a given k-point,
7 : !! using the lobpcg method
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) 1998-2026 ABINIT group (LB)
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 : !! SOURCE
20 :
21 : #if defined HAVE_CONFIG_H
22 : #include "config.h"
23 : #endif
24 :
25 : #include "abi_common.h"
26 :
27 : module m_lobpcgwf_cprj
28 :
29 : use, intrinsic :: iso_c_binding
30 : use defs_basis
31 : use m_abicore
32 : use m_lobpcg
33 : use m_errors
34 : use m_time
35 : use m_xomp
36 : use m_fstrings
37 : use m_xg
38 : use m_xg_nonlop
39 : use m_xgTransposer
40 : use m_lobpcg2_cprj
41 : use m_dtset
42 :
43 : use defs_abitypes, only : mpi_type
44 : use m_hamiltonian, only : gs_hamiltonian_type
45 : use m_getghc, only : multithreaded_getghc
46 : use m_pawcprj, only : pawcprj_type
47 :
48 : implicit none
49 :
50 : private
51 :
52 : integer, parameter :: l_tim_getghc=5
53 :
54 : ! For use in getghc_gsc1
55 : integer, save :: l_prtvol
56 : type(mpi_type),pointer,save :: l_mpi_enreg
57 : type(gs_hamiltonian_type),pointer,save :: l_gs_hamk
58 :
59 : public :: lobpcgwf2_cprj
60 :
61 : contains
62 : !!***
63 :
64 4036 : subroutine lobpcgwf2_cprj(cg,dtset,eig,occ,enl_out,gs_hamk,isppol,ikpt,inonsc,istep,kinpw,mpi_enreg,&
65 4036 : & nband,npw,nspinor,prtvol,resid,nbdbuf,xg_nonlop)
66 :
67 :
68 : use m_cgtools, only : dotprod_g
69 :
70 : !Arguments ------------------------------------
71 : integer,intent(in) :: nband,npw,prtvol,nspinor
72 : integer,intent(in) :: isppol,ikpt,inonsc,istep,nbdbuf
73 : type(gs_hamiltonian_type),target,intent(inout) :: gs_hamk
74 : type(dataset_type) ,intent(in ) :: dtset
75 : type(mpi_type) ,target,intent(in) :: mpi_enreg
76 : real(dp) ,target,intent(inout) :: cg(2,nspinor*nband*npw)
77 : real(dp) ,intent(in ) :: kinpw(npw)
78 : real(dp) ,target,intent( out) :: resid(nband)
79 : real(dp) ,intent( out) :: enl_out(nband)
80 : real(dp) ,target,intent( out) :: eig(nband)
81 : real(dp) ,target,intent(in ) :: occ(nband)
82 : type(xg_nonlop_t) ,target,intent(in ) :: xg_nonlop
83 :
84 : !Local variables-------------------------------
85 :
86 : type(xgBlock_t) :: xgx0
87 : !type(xgBlock_t) :: cprj_xgx0
88 : type(xg_t) :: cprj_xgx0
89 : type(xgBlock_t) :: xgeigen
90 : type(xgBlock_t) :: xgresidu
91 : type(xgBlock_t) :: xgocc
92 : type(xgBlock_t) :: xgenl
93 : type(xgBlock_t) :: xg_precond,xg_kin
94 4036 : type(lobpcg_t) :: lobpcg
95 :
96 : integer :: space, space_cprj, blockdim, cprjdim, nband_cprj
97 : integer :: me_g0, me_g0_fft
98 :
99 : integer, parameter :: tim_lobpcgwf2 = 2030
100 : real(dp) :: tsec(2)
101 :
102 4036 : real(dp), allocatable :: occ_tmp(:)
103 : ! Important things for NC
104 4036 : real(dp), allocatable :: pcon(:),kin(:)
105 : ! real(dp), allocatable :: cprj_contiguous(:,:)
106 :
107 : ! *********************************************************************
108 :
109 4036 : call timab(tim_lobpcgwf2,1,tsec)
110 :
111 : ! Set module variables
112 4036 : l_prtvol = prtvol
113 4036 : l_mpi_enreg => mpi_enreg
114 4036 : l_gs_hamk => gs_hamk
115 :
116 4036 : cprjdim = xg_nonlop%cprjdim
117 :
118 : !Variables
119 4036 : blockdim=nband/dtset%nblock_lobpcg
120 4036 : if (blockdim/=mpi_enreg%nproc_band*mpi_enreg%bandpp) then
121 0 : ABI_ERROR('blockdim is not consistent with nproc_band and bandpp')
122 : end if
123 4036 : nband_cprj=nband/mpi_enreg%nproc_band
124 :
125 : !Depends on istwfk
126 4036 : if ( gs_hamk%istwf_k > 1 ) then ! Real only
127 : ! SPACE_CR mean that we have complex numbers but no re*im terms only re*re
128 : ! and im*im so that a vector of complex is consider as a long vector of real
129 : ! therefore the number of data is (2*npw*nspinor)*nband
130 : ! This space is completely equivalent to SPACE_R but will correctly set and
131 : ! get the array data into the xgBlock
132 2124 : space = SPACE_CR
133 : else ! complex
134 1912 : space = SPACE_C
135 : end if
136 4036 : space_cprj = xg_nonlop%space_cprj
137 :
138 : !For kinetic part of the Hamiltonian
139 12108 : ABI_MALLOC(kin,(npw))
140 4036 : call build_kin(kin,kinpw,npw)
141 4036 : call xgBlock_map_1d(xg_kin,kin,SPACE_R,npw)
142 :
143 : !For preconditionning
144 8072 : ABI_MALLOC(pcon,(npw))
145 4036 : call build_pcon(pcon,kinpw,npw)
146 :
147 : ! Local variables for lobpcg
148 4036 : me_g0 = -1
149 4036 : me_g0_fft = -1
150 4036 : if (space==SPACE_CR) then
151 2124 : me_g0 = 0
152 2124 : me_g0_fft = 0
153 2124 : if (gs_hamk%istwf_k == 2) then
154 339 : if (l_mpi_enreg%me_g0 == 1) me_g0 = 1
155 339 : if (l_mpi_enreg%me_g0_fft == 1) me_g0_fft = 1
156 : end if
157 : end if
158 4036 : call xgBlock_map(xgx0,cg,space,npw*nspinor,nband,l_mpi_enreg%comm_band,me_g0=me_g0)
159 :
160 4036 : call xgBlock_map_1d(xg_precond,pcon,SPACE_R,npw)
161 :
162 4036 : call xgBlock_map_1d(xgeigen,eig,SPACE_R,nband)
163 :
164 4036 : call xgBlock_map_1d(xgresidu,resid,SPACE_R,nband)
165 :
166 4036 : call xgBlock_map_1d(xgenl,enl_out,SPACE_R,nband)
167 :
168 4036 : call xg_init(cprj_xgx0,space_cprj,xg_nonlop%cprjdim,nband_cprj*nspinor,comm=l_mpi_enreg%comm_band)
169 :
170 : ! Occupancies in LOBPCG are used for convergence criteria only
171 4036 : if (dtset%nbdbuf==-101.and.nspinor==1.and.dtset%nsppol==1) then
172 0 : ABI_MALLOC(occ_tmp,(nband))
173 0 : occ_tmp(:) = half*occ(:)
174 0 : call xgBlock_map_1d(xgocc,occ_tmp,SPACE_R,nband,gpu_option=dtset%gpu_option)
175 : else
176 4036 : call xgBlock_map_1d(xgocc,occ,SPACE_R,nband,gpu_option=dtset%gpu_option)
177 : end if
178 :
179 : call lobpcg_init(lobpcg,mpi_enreg%bandpp,nband,npw*nspinor,cprjdim,blockdim,dtset%tolwfr_diago,dtset%nline,&
180 : space,space_cprj,l_mpi_enreg%comm_band,dtset%paral_kgb,xg_nonlop,l_mpi_enreg%comm_spinorfft,l_mpi_enreg%comm_band,&
181 4036 : me_g0,me_g0_fft)
182 :
183 : ! Run lobpcg
184 : call lobpcg_run_cprj(lobpcg,xgx0,cprj_xgx0%self,xg_getghc,xg_kin,xg_precond,xgeigen,xgocc,xgresidu,xgenl,&
185 4036 : prtvol,nspinor,isppol,ikpt,inonsc,istep,nbdbuf)
186 :
187 4036 : if (allocated(occ_tmp)) then
188 0 : ABI_FREE(occ_tmp)
189 : end if
190 : ! Free preconditionning since not needed anymore
191 4036 : ABI_FREE(pcon)
192 4036 : ABI_FREE(kin)
193 :
194 4036 : call xg_free(cprj_xgx0)
195 :
196 : ! Free lobpcg
197 4036 : call lobpcg_free(lobpcg)
198 :
199 4036 : call timab(tim_lobpcgwf2,2,tsec)
200 :
201 : DBG_EXIT("COLL")
202 :
203 8072 : end subroutine lobpcgwf2_cprj
204 : !!***
205 :
206 : !!****f* m_lobpcg/xg_getghc
207 : !! NAME
208 : !! xg_getghc
209 : !!
210 : !! FUNCTION
211 : !! This routine computes H|C> and possibly S|C> for a given wave function C.
212 : !! It acts as a driver for getghc, taken into account parallelism, multithreading, etc.
213 : !!
214 : !! SIDE EFFECTS
215 : !! X <type(xgBlock_t)>= memory block containing |C>
216 : !! AX <type(xgBlock_t)>= memory block containing H|C>
217 : !!
218 : !! SOURCE
219 : !
220 75096 : subroutine xg_getghc(X,AX)
221 :
222 : !Arguments ------------------------------------
223 : type(xgBlock_t), intent(inout) :: X
224 : type(xgBlock_t), intent(inout) :: AX
225 :
226 : !Local variables-------------------------------
227 : !scalars
228 : integer :: blockdim
229 : integer :: spacedim
230 : integer,parameter :: sij_opt=0,cpopt=-1,type_calc=1 ! Compute local part only
231 : ! integer :: iatom,iband,ispinor,cprj_index,cprj_rows,cprj_cols,ncpgr,nlmn
232 : real(dp) :: eval
233 1051344 : type(pawcprj_type) :: cprj_dum(l_gs_hamk%natom,1)
234 : !arrays
235 75096 : real(dp), pointer :: cg(:,:)
236 75096 : real(dp), pointer :: ghc(:,:)
237 : real(dp) :: gsc(1,1),gvnlxc(1,1)
238 :
239 : ! *********************************************************************
240 :
241 75096 : call xgBlock_getSize(X,spacedim,blockdim)
242 75096 : call xgBlock_check(X,AX)
243 :
244 75096 : call xgBlock_reverseMap(X,cg,rows=1,cols=spacedim*blockdim)
245 75096 : call xgBlock_reverseMap(AX,ghc,rows=1,cols=spacedim*blockdim)
246 :
247 : ! Apply only local part of the Hamiltonian
248 : call multithreaded_getghc(cpopt,cg,cprj_dum,ghc,gsc,&
249 75096 : l_gs_hamk,gvnlxc,eval,l_mpi_enreg,blockdim,l_prtvol,sij_opt,l_tim_getghc,type_calc)
250 :
251 150192 : end subroutine xg_getghc
252 : !!***
253 :
254 4036 : subroutine build_pcon(pcon,kinpw,npw)
255 :
256 : integer,intent(in) :: npw
257 : real(dp),intent(in) :: kinpw(:)
258 : real(dp),intent(out) :: pcon(:)
259 :
260 : integer :: ipw
261 :
262 : !$omp parallel do schedule(static), shared(pcon,kinpw)
263 269445 : do ipw=1,npw
264 269445 : if(kinpw(ipw)>huge(0.0_dp)*1.d-11) then
265 5848 : pcon(ipw)=0.d0
266 : else
267 : pcon(ipw) = (27+kinpw(ipw)*(18+kinpw(ipw)*(12+8*kinpw(ipw)))) &
268 259561 : & / (27+kinpw(ipw)*(18+kinpw(ipw)*(12+8*kinpw(ipw))) + 16*kinpw(ipw)**4)
269 : end if
270 : end do
271 :
272 4036 : end subroutine build_pcon
273 :
274 4036 : subroutine build_kin(kin,kinpw,npw)
275 :
276 : integer,intent(in) :: npw
277 : real(dp),intent(in) :: kinpw(:)
278 : real(dp),intent(out) :: kin(:)
279 :
280 : integer :: ipw
281 :
282 : !$omp parallel do schedule(static), shared(kin,kinpw)
283 269445 : do ipw=1,npw
284 269445 : if(kinpw(ipw)>huge(0.0_dp)*1.d-11) then
285 5848 : kin(ipw)=0.d0
286 : else
287 259561 : kin(ipw) = kinpw(ipw)
288 : end if
289 : end do
290 :
291 4036 : end subroutine build_kin
292 :
293 : end module m_lobpcgwf_cprj
294 : !!***
|
