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