Line data Source code
1 : !!****m* ABINIT/m_cgcprj
2 : !! NAME
3 : !! m_cgcprj
4 : !!
5 : !! FUNCTION
6 : !! Functions operating on wavefunctions in the cg+cprj representation.
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2008-2026 ABINIT group (XG)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : module m_cgcprj
23 :
24 : use defs_basis
25 : use defs_abitypes
26 : use m_cgtools
27 : use m_errors
28 : use m_xmpi
29 :
30 : use m_pawtab, only : pawtab_type
31 : use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_get, pawcprj_free, pawcprj_lincom
32 :
33 : implicit none
34 :
35 : private
36 : !!***
37 :
38 : public :: dotprod_set_cgcprj
39 : public :: dotprodm_sumdiag_cgcprj
40 : public :: lincom_cgcprj
41 : public :: cgcprj_cholesky
42 : !!***
43 :
44 : contains
45 : !!***
46 :
47 : !!****f* ABINIT/dotprod_set_cgcprj
48 : !! NAME
49 : !! dotprod_set_cgcprj
50 : !!
51 : !! FUNCTION
52 : !! For one k point and spinpol, compute the matrix of scalar products between two sets of nband wavefunctions
53 : !! that are known in the cg+cprj representation
54 : !! Smn=<wf1m|wf2n>
55 : !!
56 : !! Can also treat the case of the computation of scalar products within one set of nband wavefunctions
57 : !! without recomputing already computed matrix elements (define hermitian=1)
58 : !!
59 : !! This implementation is NOT band-parallelized
60 : !! Also, it is far of being optimal at the level of linear algebra, and involves extra copying
61 : !! that are detrimental for performance...
62 : !!
63 : !! INPUTS
64 : !! atindx1(natom)=index table for atoms, inverse of atindx
65 : !! cg1(2,mcg1)= plane wave wavefunction coefficients for the first set of wavefunctions (all k points and spinpol)
66 : !! cg2(2,mcg2)= plane wave wavefunction coefficients for the second set of wavefunctions (all k points and spinpol)
67 : !! cprj1(natom,mcprj1) <type(pawcprj_type)>= projected input wave functions <Proj_i|Cnk> with NL projectors in the first set,
68 : !! cprj2(natom,mcprj2) <type(pawcprj_type)>= projected input wave functions <Proj_i|Cnk> with NL projectors in the second set,
69 : !! dimcprj(natom)=number of lmn components in the <p_{lmn}^i|\psi> for the i-th atom
70 : !! hermitian= if 1, consider that the Smn matrix is hermitian, and do not recompute already computed matrix elements.
71 : !! ibg1=shift in cprj1 array to locate current k-point.
72 : !! Might be 0, in which case cprj1 is not copied internally, which saves some time/space
73 : !! ibg2=shift in cprj2 array to locate current k-point.
74 : !! Might be 0, in which case cprj2 is not copied internally, which saves some time/space
75 : !! icg1=shift in cg1 array to locate current k-point
76 : !! icg2=shift in cg2 array to locate current k-point
77 : !! ikpt=current k point index
78 : !! isppol=current spin polarization index
79 : !! istwf=input option parameter that describes the storage of wfs
80 : !! mcg1=second dimension of cg1 array (mpw*nspinor*mband1*mkmem*nsppol)
81 : !! mcg2=second dimension of cg2 array (mpw*nspinor*mband2*mkmem*nsppol)
82 : !! mcprj1=second dimension of cprj1 array
83 : !! mcprj2=second dimension of cprj2 array
84 : !! mkmem=number of k points which can fit in memory
85 : !! mpi_enreg=information about MPI parallelization
86 : !! natom=number of atoms
87 : !! nattyp(ntypat)=number of atoms of each type in cell.
88 : !! nbd1=number of bands for the first set of wavefunctions
89 : !! nbd2=number of bands for the second set of wavefunctions
90 : !! npw=number of planewaves in basis at this k point
91 : !! nspinor=number of spinor components
92 : !! nsppol=number of spin polarizations
93 : !! ntypat=number of types of atoms
94 : !! pawtab(dtset%ntypat*dtset%usepaw) <type(pawtab_type)>=paw tabulated starting data
95 : !! usepaw=1 if PAW is activated
96 : !!
97 : !! OUTPUT
98 : !! smn(2,nbd1,nbd2)=matrix of scalar products between the first set of wavefunctions and the second set of wavefunctions
99 : !!
100 : !! SIDE EFFECTS
101 : !!
102 : !! SOURCE
103 :
104 282 : subroutine dotprod_set_cgcprj(atindx1,cg1,cg2,cprj1,cprj2,dimcprj,hermitian,&
105 : & ibg1,ibg2,icg1,icg2,ikpt,isppol,istwf,mband,mcg1,mcg2,mcprj1,mcprj2,mkmem,&
106 282 : & mpi_enreg,natom,nattyp,nbd1,nbd2,npw,nspinor,nsppol,ntypat,pawtab,smn,usepaw)
107 :
108 : !Arguments ------------------------------------
109 : !scalars
110 : integer, intent(in) :: hermitian,ibg1,ibg2,icg1,icg2,ikpt,isppol,istwf
111 : integer, intent(in) :: mkmem,mband,mcg1,mcg2,mcprj1,mcprj2
112 : integer, intent(in) :: natom,nbd1,nbd2,npw,nspinor,nsppol,ntypat,usepaw
113 : type(MPI_type),intent(in) :: mpi_enreg
114 : !arrays
115 : integer, intent(in) :: atindx1(natom),dimcprj(natom),nattyp(ntypat)
116 : real(dp), intent(in) :: cg1(2,mcg1),cg2(2,mcg2)
117 : real(dp), intent(out) :: smn(2,nbd1,nbd2)
118 : type(pawcprj_type),intent(in) :: cprj1(natom,mcprj1),cprj2(natom,mcprj2)
119 : type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw)
120 :
121 : !Local variables-------------------------------
122 : !scalars
123 : integer :: ia,iat,itypat,ibd1,ibd2,icgb1,icgb2,ier,ig,ii,i1,i2,iorder
124 : integer :: ilmn1,ilmn2,klmn,max_nbd2,nbd
125 : real(dp) :: dotr,doti
126 : !arrays
127 282 : real(dp),allocatable :: cwavef1(:,:),cwavef2(:,:),proj(:,:,:)
128 282 : real(dp),allocatable :: eigval(:),eigvec(:,:,:),matrx(:,:),zhpev1(:,:),zhpev2(:)
129 282 : type(pawcprj_type),allocatable :: cprj1_k(:,:),cprj2_k(:,:)
130 :
131 : ! *************************************************************************
132 :
133 : !DEBUG
134 : !write(std_out,*)' dotprod_set_cgcprj : enter '
135 : !write(std_out,*)' dotprod_set_cgcprj : npw, nspinor=',npw,nspinor
136 : !write(std_out,*)' dotprod_set_cgcprj : usepaw,nbd1,nbd2=',usepaw,nbd1,nbd2
137 : !call flush(std_out)
138 : !ENDDEBUG
139 :
140 282 : if(hermitian==1)then
141 141 : if(nbd1/=nbd2)then
142 0 : ABI_ERROR(' With hermitian==1, nb1 and nb2 must be equal ')
143 : end if
144 : end if
145 :
146 846 : ABI_MALLOC(cwavef1,(2,npw*nspinor))
147 564 : ABI_MALLOC(cwavef2,(2,npw*nspinor))
148 282 : if(usepaw==1) then
149 1920 : ABI_MALLOC(cprj1_k,(natom,nspinor*nbd1))
150 1920 : ABI_MALLOC(cprj2_k,(natom,nspinor*nbd2))
151 120 : iorder=0 ! There is no change of ordering of cprj when copying wavefunctions
152 : end if
153 282 : if(usepaw==1 .and. ibg1/=0) then
154 0 : call pawcprj_alloc(cprj1_k,cprj1(1,1)%ncpgr,dimcprj)
155 : end if
156 282 : if(usepaw==1 .and. ibg2/=0) then
157 0 : call pawcprj_alloc(cprj2_k,cprj1(1,1)%ncpgr,dimcprj)
158 : end if
159 :
160 282 : icgb1=icg1
161 1380 : do ibd1=1,nbd1
162 :
163 : ! Extract wavefunction information
164 197602 : do ig=1,npw*nspinor
165 196504 : cwavef1(1,ig)=cg1(1,ig+icgb1)
166 197602 : cwavef1(2,ig)=cg1(2,ig+icgb1)
167 : end do
168 1098 : if(usepaw==1 .and. ibg1/=0) then
169 : call pawcprj_get(atindx1,cprj1_k,cprj1,natom,1,ibg1,ikpt,iorder,isppol,mband,&
170 : & mkmem,natom,nbd1,nbd1,nspinor,nsppol,0,&
171 0 : & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
172 : end if
173 :
174 1098 : icgb2=icg2
175 1098 : max_nbd2=nbd2
176 1098 : if(hermitian==1)max_nbd2=ibd1
177 4644 : do ibd2=1,max_nbd2
178 :
179 : ! XG171222 Note that this copy step, being inside the ibd1 loop, is quite detrimental.
180 : ! It might be reduced by copying several cwavef2, and use a ZGEMM type of approach.
181 :
182 : ! Extract wavefunction information
183 574324 : do ig=1,npw*nspinor
184 570778 : cwavef2(1,ig)=cg2(1,ig+icgb2)
185 574324 : cwavef2(2,ig)=cg2(2,ig+icgb2)
186 : end do
187 :
188 3546 : if(usepaw==1 .and. ibg2/=0) then
189 : call pawcprj_get(atindx1,cprj2_k,cprj2,natom,1,ibg2,ikpt,iorder,isppol,mband,&
190 : & mkmem,natom,nbd2,nbd2,nspinor,nsppol,0,&
191 0 : & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
192 : end if
193 :
194 : ! Calculate Smn=<cg1|cg2>
195 3546 : call dotprod_g(dotr,doti,istwf,npw*nspinor,2,cwavef1,cwavef2,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
196 :
197 3546 : if(usepaw==1) then
198 : ia =0
199 3120 : do itypat=1,ntypat
200 1560 : if(ibg1/=0 .and. ibg2/=0)then
201 0 : do iat=1+ia,nattyp(itypat)+ia
202 0 : do ilmn1=1,pawtab(itypat)%lmn_size
203 0 : do ilmn2=1,ilmn1
204 0 : klmn=((ilmn1-1)*ilmn1)/2+ilmn2
205 : dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd1)%cp(1,ilmn1)*cprj2_k(iat,ibd2)%cp(1,ilmn2)+&
206 0 : & cprj1_k(iat,ibd1)%cp(2,ilmn1)*cprj2_k(iat,ibd2)%cp(2,ilmn2))
207 : doti=doti+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd1)%cp(1,ilmn1)*cprj2_k(iat,ibd2)%cp(2,ilmn2)-&
208 0 : & cprj1_k(iat,ibd1)%cp(2,ilmn1)*cprj2_k(iat,ibd2)%cp(1,ilmn2))
209 : end do
210 0 : do ilmn2=ilmn1+1,pawtab(itypat)%lmn_size
211 0 : klmn=((ilmn2-1)*ilmn2)/2+ilmn1
212 : dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd1)%cp(1,ilmn1)*cprj2_k(iat,ibd2)%cp(1,ilmn2)+&
213 0 : & cprj1_k(iat,ibd1)%cp(2,ilmn1)*cprj2_k(iat,ibd2)%cp(2,ilmn2))
214 : doti=doti+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd1)%cp(1,ilmn1)*cprj2_k(iat,ibd2)%cp(2,ilmn2)-&
215 0 : & cprj1_k(iat,ibd1)%cp(2,ilmn1)*cprj2_k(iat,ibd2)%cp(1,ilmn2))
216 : end do
217 : end do
218 : end do
219 1560 : else if(ibg1/=0 .and. ibg2==0)then
220 0 : do iat=1+ia,nattyp(itypat)+ia
221 0 : do ilmn1=1,pawtab(itypat)%lmn_size
222 0 : do ilmn2=1,ilmn1
223 0 : klmn=((ilmn1-1)*ilmn1)/2+ilmn2
224 : dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd1)%cp(1,ilmn1)*cprj2(iat,ibd2)%cp(1,ilmn2)+&
225 0 : & cprj1_k(iat,ibd1)%cp(2,ilmn1)*cprj2(iat,ibd2)%cp(2,ilmn2))
226 : doti=doti+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd1)%cp(1,ilmn1)*cprj2(iat,ibd2)%cp(2,ilmn2)-&
227 0 : & cprj1_k(iat,ibd1)%cp(2,ilmn1)*cprj2(iat,ibd2)%cp(1,ilmn2))
228 : end do
229 0 : do ilmn2=ilmn1+1,pawtab(itypat)%lmn_size
230 0 : klmn=((ilmn2-1)*ilmn2)/2+ilmn1
231 : dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd1)%cp(1,ilmn1)*cprj2(iat,ibd2)%cp(1,ilmn2)+&
232 0 : & cprj1_k(iat,ibd1)%cp(2,ilmn1)*cprj2(iat,ibd2)%cp(2,ilmn2))
233 : doti=doti+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd1)%cp(1,ilmn1)*cprj2(iat,ibd2)%cp(2,ilmn2)-&
234 0 : & cprj1_k(iat,ibd1)%cp(2,ilmn1)*cprj2(iat,ibd2)%cp(1,ilmn2))
235 : end do
236 : end do
237 : end do
238 1560 : else if(ibg1==0 .and. ibg2/=0)then
239 0 : do iat=1+ia,nattyp(itypat)+ia
240 0 : do ilmn1=1,pawtab(itypat)%lmn_size
241 0 : do ilmn2=1,ilmn1
242 0 : klmn=((ilmn1-1)*ilmn1)/2+ilmn2
243 : dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj1(iat,ibd1)%cp(1,ilmn1)*cprj2_k(iat,ibd2)%cp(1,ilmn2)+&
244 0 : & cprj1(iat,ibd1)%cp(2,ilmn1)*cprj2_k(iat,ibd2)%cp(2,ilmn2))
245 : doti=doti+pawtab(itypat)%sij(klmn)*(cprj1(iat,ibd1)%cp(1,ilmn1)*cprj2_k(iat,ibd2)%cp(2,ilmn2)-&
246 0 : & cprj1(iat,ibd1)%cp(2,ilmn1)*cprj2_k(iat,ibd2)%cp(1,ilmn2))
247 : end do
248 0 : do ilmn2=ilmn1+1,pawtab(itypat)%lmn_size
249 0 : klmn=((ilmn2-1)*ilmn2)/2+ilmn1
250 : dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj1(iat,ibd1)%cp(1,ilmn1)*cprj2_k(iat,ibd2)%cp(1,ilmn2)+&
251 0 : & cprj1(iat,ibd1)%cp(2,ilmn1)*cprj2_k(iat,ibd2)%cp(2,ilmn2))
252 : doti=doti+pawtab(itypat)%sij(klmn)*(cprj1(iat,ibd1)%cp(1,ilmn1)*cprj2_k(iat,ibd2)%cp(2,ilmn2)-&
253 0 : & cprj1(iat,ibd1)%cp(2,ilmn1)*cprj2_k(iat,ibd2)%cp(1,ilmn2))
254 : end do
255 : end do
256 : end do
257 1560 : else if(ibg1==0 .and. ibg2==0)then
258 4680 : do iat=1+ia,nattyp(itypat)+ia
259 29640 : do ilmn1=1,pawtab(itypat)%lmn_size
260 137280 : do ilmn2=1,ilmn1
261 112320 : klmn=((ilmn1-1)*ilmn1)/2+ilmn2
262 : dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj1(iat,ibd1)%cp(1,ilmn1)*cprj2(iat,ibd2)%cp(1,ilmn2)+&
263 112320 : & cprj1(iat,ibd1)%cp(2,ilmn1)*cprj2(iat,ibd2)%cp(2,ilmn2))
264 : doti=doti+pawtab(itypat)%sij(klmn)*(cprj1(iat,ibd1)%cp(1,ilmn1)*cprj2(iat,ibd2)%cp(2,ilmn2)-&
265 137280 : & cprj1(iat,ibd1)%cp(2,ilmn1)*cprj2(iat,ibd2)%cp(1,ilmn2))
266 : end do
267 115440 : do ilmn2=ilmn1+1,pawtab(itypat)%lmn_size
268 87360 : klmn=((ilmn2-1)*ilmn2)/2+ilmn1
269 : dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj1(iat,ibd1)%cp(1,ilmn1)*cprj2(iat,ibd2)%cp(1,ilmn2)+&
270 87360 : & cprj1(iat,ibd1)%cp(2,ilmn1)*cprj2(iat,ibd2)%cp(2,ilmn2))
271 : doti=doti+pawtab(itypat)%sij(klmn)*(cprj1(iat,ibd1)%cp(1,ilmn1)*cprj2(iat,ibd2)%cp(2,ilmn2)-&
272 112320 : & cprj1(iat,ibd1)%cp(2,ilmn1)*cprj2(iat,ibd2)%cp(1,ilmn2))
273 : end do
274 : end do
275 : end do
276 : end if
277 3120 : ia=ia+nattyp(itypat)
278 : end do
279 : end if
280 3546 : smn(1,ibd1,ibd2)=dotr
281 3546 : smn(2,ibd1,ibd2)=doti
282 : ! End loop over bands ibd2
283 4644 : icgb2=icgb2+npw*nspinor
284 :
285 : end do
286 :
287 : ! End loop over bands ibd1
288 1380 : icgb1=icgb1+npw*nspinor
289 : end do
290 :
291 : !Complete the matrix if hermitian
292 282 : if(hermitian==1)then
293 549 : do ibd1=1,nbd1-1
294 1365 : do ibd2=ibd1+1,nbd2
295 816 : smn(1,ibd1,ibd2)= smn(1,ibd2,ibd1)
296 1224 : smn(2,ibd1,ibd2)=-smn(2,ibd2,ibd1)
297 : end do
298 : end do
299 : end if
300 :
301 : !DEBUG
302 : !write(std_out,*)' smn=',smn
303 : !ENDDEBUG
304 :
305 : !====== Debugging section ==========
306 : if(.false.)then
307 : !DEBUG
308 : !Compute the eigenvalues of the projector S herm(S) or herm(S) S, depending on which has lowest dimension.
309 : !write(std_out,*)' dotprod_set_cgcprj : compute the projector matrix '
310 : nbd=min(nbd1,nbd2)
311 : ABI_MALLOC(proj,(2,nbd,nbd))
312 : proj(:,:,:)=zero
313 : if(nbd1<=nbd2)then
314 : do ibd1=1,nbd1
315 : do ibd2=1,nbd2
316 : proj(1,:,ibd1)=proj(1,:,ibd1)+smn(1,:,ibd2)*smn(1,ibd1,ibd2)+smn(2,:,ibd2)*smn(2,ibd1,ibd2)
317 : proj(2,:,ibd1)=proj(2,:,ibd1)-smn(1,:,ibd2)*smn(2,ibd1,ibd2)+smn(2,:,ibd2)*smn(1,ibd1,ibd2)
318 : end do
319 : end do
320 : else
321 : do ibd2=1,nbd2
322 : do ibd1=1,nbd1
323 : proj(1,:,ibd2)=proj(1,:,ibd2)+smn(1,ibd1,:)*smn(1,ibd1,ibd2)+smn(2,ibd1,:)*smn(2,ibd1,ibd2)
324 : proj(2,:,ibd2)=proj(2,:,ibd2)+smn(1,ibd1,:)*smn(2,ibd1,ibd2)-smn(2,ibd1,:)*smn(1,ibd1,ibd2)
325 : end do
326 : end do
327 : end if
328 :
329 : !write(std_out,*)' proj=',proj
330 :
331 : !write(std_out,*)' dotprod_set_cgcprj : compute the eigenvalues of the projector '
332 : ABI_MALLOC(matrx,(2,(nbd*(nbd+1))/2))
333 : ii=1
334 : do i2=1,nbd
335 : do i1=1,i2
336 : matrx(1,ii)=proj(1,i1,i2)
337 : matrx(2,ii)=proj(2,i1,i2)
338 : ii=ii+1
339 : end do
340 : end do
341 :
342 : ABI_MALLOC(zhpev1,(2,2*nbd-1))
343 : ABI_MALLOC(zhpev2,(3*nbd-2))
344 : ABI_MALLOC(eigval,(nbd))
345 : ABI_MALLOC(eigvec,(2,nbd,nbd))
346 :
347 : call ZHPEV ('V','U',nbd,matrx,eigval,eigvec,nbd,zhpev1,zhpev2,ier)
348 :
349 : !write(std_out,*)' eigval=',eigval
350 :
351 : ABI_FREE(matrx)
352 : ABI_FREE(zhpev1)
353 : ABI_FREE(zhpev2)
354 : ABI_FREE(eigval)
355 : ABI_FREE(eigvec)
356 :
357 : ABI_FREE(proj)
358 : !stop
359 : !ENDDEBUG
360 : end if
361 : !====== End of debugging section ==========
362 :
363 282 : ABI_FREE(cwavef1)
364 282 : ABI_FREE(cwavef2)
365 282 : if(usepaw==1) then
366 120 : if(ibg1/=0)then
367 0 : call pawcprj_free(cprj1_k)
368 : end if
369 120 : if(ibg2/=0)then
370 0 : call pawcprj_free(cprj2_k)
371 : end if
372 1080 : ABI_FREE(cprj1_k)
373 1080 : ABI_FREE(cprj2_k)
374 : end if
375 :
376 282 : end subroutine dotprod_set_cgcprj
377 : !!***
378 :
379 : !!****f* ABINIT/dotprodm_sumdiag_cgcprj
380 : !!
381 : !! NAME
382 : !! dotprodm_sumdiag_cgcprj
383 : !!
384 : !! FUNCTION
385 : !! For one k point and spinpol, compute the matrix of sum of diagonal scalar products
386 : !! between sets of nband wavefunctions that are known in the cg+cprj representation
387 : !! The sets of wavefunctions must be contained in a big array of sets of wavefunctions.
388 : !! Sij=Sum_m <wfm(set i)|wfm(set j)>
389 : !!
390 : !! Can also treat the case of the computation of scalar products within one set of wavefunctions
391 : !!
392 : !! This implementation is NOT band-parallelized
393 : !! Also, it is far of being optimal at the level of linear algebra, and involves extra copying
394 : !! that are detrimental for performance...
395 : !!
396 : !! INPUTS
397 : !! atindx1(natom)=index table for atoms, inverse of atindx
398 : !! cg_set(2,mcg,mset)= plane wave wavefunction coefficients for several sets of wavefunctions (all k points and spins)
399 : !! cprj_set(natom,mcprj,mset) <type(pawcprj_type)>= projected input wave functions <Proj_i|Cnk>
400 : !! with NL projectors in the different sets
401 : !! dimcprj(natom)=number of lmn components in the <p_{lmn}^i|\psi> for the i-th atom
402 : !! ibg=shift in cprj_set array to locate current k-point
403 : !! icg=shift in cg_set array to locate current k-point
404 : !! ikpt=current k point index
405 : !! isppol=current spin polarization index
406 : !! istwf=input option parameter that describes the storage of wfs
407 : !! mband=maximum number of bands (used in the dimensioning of cprj_set)
408 : !! mcg=second dimension of cg array (mpw*nspinor*mband*mkmem*nsppol)
409 : !! mcprj=second dimension of cprj array
410 : !! mkmem=number of k points which can fit in memory
411 : !! mpi_enreg=information about MPI parallelization
412 : !! mset=third dimension of cg_set and cprj_set, maximum number of sets
413 : !! natom=number of atoms
414 : !! nattyp(ntypat)=number of atoms of each type in cell.
415 : !! nbd=number of bands for each set of wavefunctions
416 : !! npw=number of planewaves in basis at this k point
417 : !! nset1=number of sets of wavefunctions to be considered in the left side of the scalar products
418 : !! nset2=number of sets of wavefunctions to be considered in the right side of the scalar products
419 : !! nspinor=number of spinor components
420 : !! nsppol=number of spin polarizations
421 : !! ntypat=number of types of atoms
422 : !! pawtab(dtset%ntypat*dtset%usepaw) <type(pawtab_type)>=paw tabulated starting data
423 : !! shift_set1=shift that defines the first set of wavefunctions to be considered in the left side of the scalar products
424 : !! shift_set2=shift that defines the first set of wavefunctions to be considered in the right side of the scalar products
425 : !! usepaw=1 if PAW is activated
426 : !!
427 : !! OUTPUT
428 : !! smn(2,nset1,nset2)=matrix of sum of diagonal scalar products between the first set
429 : !! of wavefunctions and the second set of wavefunctions
430 : !!
431 : !! SIDE EFFECTS
432 : !!
433 : !! SOURCE
434 :
435 41 : subroutine dotprodm_sumdiag_cgcprj(atindx1,cg_set,cprj_set,dimcprj,&
436 : & ibg,icg,ikpt,isppol,istwf,mband,mcg,mcprj,mkmem,&
437 41 : & mpi_enreg,mset,natom,nattyp,nbd,npw,nset1,nset2,nspinor,nsppol,ntypat,&
438 41 : & shift_set1,shift_set2,pawtab,smn,usepaw)
439 :
440 : !Arguments ------------------------------------
441 : !scalars
442 : integer, intent(in) :: ibg,icg,ikpt,isppol,istwf
443 : integer, intent(in) :: mband,mcg,mcprj,mkmem,mset
444 : integer, intent(in) :: natom,nbd,npw,nset1,nset2,nspinor,nsppol,ntypat
445 : integer, intent(in) :: shift_set1,shift_set2,usepaw
446 : type(MPI_type),intent(in) :: mpi_enreg
447 : !arrays
448 : integer, intent(in) :: atindx1(natom),dimcprj(natom),nattyp(ntypat)
449 : real(dp), intent(in) :: cg_set(2,mcg,mset)
450 : real(dp), intent(out) :: smn(2,nset1,nset2)
451 : type(pawcprj_type),intent(in) :: cprj_set(natom,mcprj,mset)
452 : type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw)
453 :
454 : !Local variables-------------------------------
455 : !scalars
456 : integer :: ia,iat,itypat,ibd,icgb,ig,iorder
457 : integer :: ilmn1,ilmn2,ind_set1,ind_set2,iset1,iset2,klmn
458 : real(dp) :: dotr,doti
459 : !arrays
460 41 : real(dp),allocatable :: cwavef1(:,:),cwavef2(:,:)
461 41 : type(pawcprj_type),allocatable :: cprj1_k(:,:),cprj2_k(:,:)
462 :
463 : ! *************************************************************************
464 :
465 : !DEBUG
466 : !write(std_out,*)' dotprodm_sumdiag_cgcprj : enter '
467 : !call flush(std_out)
468 : !ENDDEBUG
469 :
470 123 : ABI_MALLOC(cwavef1,(2,npw*nspinor))
471 82 : ABI_MALLOC(cwavef2,(2,npw*nspinor))
472 41 : if(usepaw==1) then
473 0 : ABI_MALLOC(cprj1_k,(natom,nspinor*nbd))
474 0 : ABI_MALLOC(cprj2_k,(natom,nspinor*nbd))
475 0 : iorder=0 ! There is no change of ordering in the copy of wavefunctions
476 0 : call pawcprj_alloc(cprj1_k,cprj_set(1,1,1)%ncpgr,dimcprj)
477 : end if
478 :
479 661 : smn(:,:,:)=zero
480 :
481 41 : icgb=icg
482 190 : do ibd=1,nbd
483 :
484 298 : do iset1=1,nset1
485 :
486 149 : ind_set1=iset1+shift_set1
487 :
488 : ! Extract wavefunction information
489 41581 : do ig=1,npw*nspinor
490 41432 : cwavef1(1,ig)=cg_set(1,ig+icgb,ind_set1)
491 41581 : cwavef1(2,ig)=cg_set(2,ig+icgb,ind_set1)
492 : end do
493 149 : if(usepaw==1) then
494 : call pawcprj_get(atindx1,cprj1_k,cprj_set(:,:,ind_set1),natom,1,ibg,ikpt,iorder,isppol,mband,&
495 : & mkmem,natom,nbd,nbd,nspinor,nsppol,0,&
496 0 : & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
497 : end if
498 :
499 873 : do iset2=1,nset2
500 :
501 575 : ind_set2=iset2+shift_set2
502 575 : if(ind_set2<ind_set1 .and. ind_set2>shift_set1)then
503 : continue ! These matrix elements have already been computed, the smn matrix will be completed later.
504 :
505 575 : else if(ind_set1/=ind_set2)then
506 :
507 : ! Extract wavefunction information
508 106714 : do ig=1,npw*nspinor
509 106288 : cwavef2(1,ig)=cg_set(1,ig+icgb,ind_set2)
510 106714 : cwavef2(2,ig)=cg_set(2,ig+icgb,ind_set2)
511 : end do
512 :
513 426 : if(usepaw==1) then
514 : call pawcprj_get(atindx1,cprj2_k,cprj_set(:,:,ind_set2),natom,1,ibg,ikpt,iorder,isppol,mband,&
515 : & mkmem,natom,nbd,nbd,nspinor,nsppol,0,&
516 0 : & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
517 : end if
518 :
519 : ! Calculate Smn=<cg1|cg2>
520 426 : call dotprod_g(dotr,doti,istwf,npw*nspinor,2,cwavef1,cwavef2,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
521 :
522 426 : if(usepaw==1) then
523 : ia =0
524 0 : do itypat=1,ntypat
525 0 : do iat=1+ia,nattyp(itypat)+ia
526 0 : do ilmn1=1,pawtab(itypat)%lmn_size
527 0 : do ilmn2=1,ilmn1
528 0 : klmn=((ilmn1-1)*ilmn1)/2+ilmn2
529 : dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd)%cp(1,ilmn1)*cprj2_k(iat,ibd)%cp(1,ilmn2)+&
530 0 : & cprj1_k(iat,ibd)%cp(2,ilmn1)*cprj2_k(iat,ibd)%cp(2,ilmn2))
531 : doti=doti+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd)%cp(1,ilmn1)*cprj2_k(iat,ibd)%cp(2,ilmn2)-&
532 0 : & cprj1_k(iat,ibd)%cp(2,ilmn1)*cprj2_k(iat,ibd)%cp(1,ilmn2))
533 : end do
534 0 : do ilmn2=ilmn1+1,pawtab(itypat)%lmn_size
535 0 : klmn=((ilmn2-1)*ilmn2)/2+ilmn1
536 : dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd)%cp(1,ilmn1)*cprj2_k(iat,ibd)%cp(1,ilmn2)+&
537 0 : & cprj1_k(iat,ibd)%cp(2,ilmn1)*cprj2_k(iat,ibd)%cp(2,ilmn2))
538 : doti=doti+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd)%cp(1,ilmn1)*cprj2_k(iat,ibd)%cp(2,ilmn2)-&
539 0 : & cprj1_k(iat,ibd)%cp(2,ilmn1)*cprj2_k(iat,ibd)%cp(1,ilmn2))
540 : end do
541 : end do
542 : end do
543 0 : ia=ia+nattyp(itypat)
544 : end do
545 : end if ! usepaw
546 :
547 : ! if(.false.)then
548 : else
549 : ! Diagonal part : no need to extract another wavefunction, and the scalar product must be real
550 :
551 : ! Calculate Smn=<cg1|cg1>
552 149 : call dotprod_g(dotr,doti,istwf,npw*nspinor,1,cwavef1,cwavef1,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
553 :
554 149 : if(usepaw==1) then
555 : ia =0
556 0 : do itypat=1,ntypat
557 0 : do iat=1+ia,nattyp(itypat)+ia
558 0 : do ilmn1=1,pawtab(itypat)%lmn_size
559 0 : do ilmn2=1,ilmn1
560 0 : klmn=((ilmn1-1)*ilmn1)/2+ilmn2
561 : dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd)%cp(1,ilmn1)*cprj1_k(iat,ibd)%cp(1,ilmn2)+&
562 0 : & cprj1_k(iat,ibd)%cp(2,ilmn1)*cprj1_k(iat,ibd)%cp(2,ilmn2))
563 : end do
564 0 : do ilmn2=ilmn1+1,pawtab(itypat)%lmn_size
565 0 : klmn=((ilmn2-1)*ilmn2)/2+ilmn1
566 : dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd)%cp(1,ilmn1)*cprj1_k(iat,ibd)%cp(1,ilmn2)+&
567 0 : & cprj1_k(iat,ibd)%cp(2,ilmn1)*cprj1_k(iat,ibd)%cp(2,ilmn2))
568 : end do
569 : end do
570 : end do
571 0 : ia=ia+nattyp(itypat)
572 : end do
573 : end if ! usepaw
574 149 : doti=zero
575 :
576 : end if ! Compare ind_set1 and ind_set2
577 :
578 575 : smn(1,iset1,iset2)=smn(1,iset1,iset2)+dotr
579 724 : smn(2,iset1,iset2)=smn(2,iset1,iset2)+doti
580 :
581 : end do ! iset2
582 : end do ! iset1
583 :
584 : ! End loop over bands ibd
585 190 : icgb=icgb+npw*nspinor
586 : end do
587 :
588 : !Complete the matrix, using its hermitian property.
589 82 : do iset1=1,nset1
590 41 : ind_set1=iset1+shift_set1
591 237 : do iset2=1,nset2
592 155 : ind_set2=iset2+shift_set2
593 196 : if(ind_set2<ind_set1 .and. ind_set2>shift_set1)then
594 0 : smn(1,iset1,iset2)= smn(1,iset2,iset1)
595 0 : smn(2,iset1,iset2)=-smn(2,iset2,iset1)
596 : end if
597 : end do
598 : end do
599 :
600 41 : ABI_FREE(cwavef1)
601 41 : ABI_FREE(cwavef2)
602 41 : if(usepaw==1) then
603 0 : call pawcprj_free(cprj1_k)
604 0 : call pawcprj_free(cprj2_k)
605 0 : ABI_FREE(cprj1_k)
606 0 : ABI_FREE(cprj2_k)
607 : end if
608 :
609 41 : end subroutine dotprodm_sumdiag_cgcprj
610 : !!***
611 :
612 : !!****f* ABINIT/lincom_cgcprj
613 : !!
614 : !! NAME
615 : !! lincom_cgcprj
616 : !!
617 : !! FUNCTION
618 : !! For one k point and spin, compute a set (size nband_out) of linear combinations of nband_in wavefunctions,
619 : !! that are known in the cg+cprj representation :
620 : !! cgout_n(:,:) <--- Sum_m [ cg_m(:,:) . alpha_mn ]
621 : !! cprjout_n(:,:) <--- Sum_m [ cprj_m(:,:) . alpha_mn ]
622 : !! If nband_out is smaller or equal to nband_in, the result might be in-place
623 : !! output in cg instead of cgout, and in cprj instead of cprjout).
624 : !! Otherwise, it is contained in the optional cgout+cprjout pair.
625 : !!
626 : !! In the present status, the cg and cgout relates to all the k points and spins, and rely on the icg index,
627 : !! while it is assumed that cprj and cprjout refer to the specific k point and spin.
628 : !! This is not coherent.
629 : !! THIS MIGHT BE CHANGED IN THE FUTURE !
630 : !!
631 : !! This implementation is NOT band-parallelized
632 : !! Also, it is far of being optimal at the level of linear algebra, and involves extra copying
633 : !! that are detrimental for performance...
634 : !!
635 : !! INPUTS
636 : !! alpha_mn(2,nband_in,nband_out)=complex matrix of coefficients of the linear combinations to be computed
637 : !! dimcprj(natom)=number of lmn components in the <p_{lmn}^i|\psi> for the i-th atom
638 : !! icg=shift in cg array to locate current k-point and spinpol (for input, and possibly for in-place output)
639 : !! inplace= if 0, output in cgout and cprjout ; if 1, output in cg and cprj
640 : !! mcg=second dimension of cg array (mpw*nspinor*mband*mkmem*nsppol)
641 : !! mcprj=second dimension of cprj array
642 : !! natom=number of atoms
643 : !! nband_in=number of bands, size of the input set of wavefunctions
644 : !! nband_out=number of bands, size of the output set of wavefunctions (should be equal to nband_in if inplace==1)
645 : !! npw=number of planewaves in basis at this k point
646 : !! nspinor=number of spinor components
647 : !! usepaw=1 if PAW is activated
648 : !! [icgout= shift in cgout array to locate current k-point and spinpol (for output)]
649 : !! [mcgout=second dimension of cgout array (mpw*nspinor*mband*mkmem*nsppol)]
650 : !! [mcprjout=second dimension of cprjout array]
651 : !!
652 : !! OUTPUT
653 : !! [cgout(2,mcgout)= plane wave wavefunction coefficients for the set of output wavefunctions]
654 : !! [cprjout(natom,mcprjout) <type(pawcprj_type)>= projected output wave functions <Proj_i|Cnk> with NL projectors]
655 : !!
656 : !! SIDE EFFECTS
657 : !! (this quantities are input, and possibly updated output when inplace==1)
658 : !! cg(2,mcg)= plane wave wavefunction coefficients for the set of input wavefunctions (all k points and spinpol)
659 : !! cprj(natom,mcprj) <type(pawcprj_type)>= projected input wave functions <Proj_i|Cnk> with NL projectors
660 : !!
661 : !! SOURCE
662 :
663 564 : subroutine lincom_cgcprj(alpha_mn,cg,cprj,dimcprj,&
664 : & icg,inplace,mcg,mcprj,natom,nband_in,nband_out,npw,nspinor,usepaw, &
665 282 : & cgout,cprjout,icgout) ! optional args
666 :
667 : !Arguments ------------------------------------
668 : !scalars
669 : integer, intent(in) :: icg,inplace,mcg,mcprj
670 : integer, intent(in) :: natom,nband_in,nband_out,npw,nspinor,usepaw
671 : integer, intent(in),optional :: icgout
672 : !arrays
673 : integer, intent(in) :: dimcprj(natom)
674 : real(dp), intent(inout) :: cg(2,mcg)
675 : real(dp), intent(in) :: alpha_mn(2,nband_in,nband_out)
676 : real(dp), intent(out),optional :: cgout(:,:)
677 : type(pawcprj_type),intent(inout) :: cprj(natom,mcprj)
678 : type(pawcprj_type),intent(inout),optional :: cprjout(:,:) ! ifort and others are buggy for optional intent(out) structured types
679 :
680 : !Local variables-------------------------------
681 : !scalars
682 : integer :: iband_in,iband_out,ii
683 : !arrays
684 282 : real(dp),allocatable :: al(:,:),cgout_(:,:)
685 282 : type(pawcprj_type),allocatable :: cprjout_(:,:)
686 :
687 : ! *************************************************************************
688 :
689 : !DEBUG
690 : !write(std_out,*)' lincom_cgcprj : enter '
691 : !write(std_out,*)' lincom_cgcprj : npw, nspinor=',npw,nspinor
692 : !write(std_out,*)' lincom_cgcprj : icgout=',icgout
693 : !ENDDEBUG
694 :
695 282 : if(inplace==0)then
696 36 : if(.not.present(cgout))then
697 0 : ABI_ERROR(' inplace==0 while .not.present(cgout) is not permitted ')
698 : end if
699 36 : if(usepaw==1) then
700 0 : if(.not.present(cprjout))then
701 0 : ABI_ERROR(' inplace==0 and usepaw==1 while .not.present(cprjout) is not permitted ')
702 : end if
703 : end if
704 : end if
705 :
706 : !Take care of the plane wave part
707 846 : ABI_MALLOC(cgout_,(2,npw*nspinor*nband_out))
708 :
709 : call zgemm('N','N',npw*nspinor,nband_out,nband_in,dcmplx(1._dp), &
710 : & cg(:,icg+1:icg+npw*nspinor*nband_in),npw*nspinor, &
711 282 : & alpha_mn,nband_in,dcmplx(0._dp),cgout_,npw*nspinor)
712 :
713 282 : if(inplace==1)then
714 483294 : cg(:,icg+1:icg+npw*nspinor*nband_out)=cgout_
715 : else
716 106500 : cgout(:,icgout+1:icgout+npw*nspinor*nband_out)=cgout_
717 : end if
718 282 : ABI_FREE(cgout_)
719 :
720 : !Take care of the cprj part
721 282 : if(usepaw==1) then
722 :
723 1920 : ABI_MALLOC(cprjout_,(natom,nspinor*nband_out))
724 120 : call pawcprj_alloc(cprjout_,cprj(1,1)%ncpgr,dimcprj)
725 360 : ABI_MALLOC(al,(2,nband_in))
726 600 : do iband_out=1,nband_out
727 480 : ii=(iband_out-1)*nspinor
728 2400 : do iband_in=1,nband_in
729 1920 : al(1,iband_in)=alpha_mn(1,iband_in,iband_out)
730 2400 : al(2,iband_in)=alpha_mn(2,iband_in,iband_out)
731 : end do
732 600 : call pawcprj_lincom(al,cprj,cprjout_(:,ii+1:ii+nspinor),nband_in)
733 : end do
734 120 : ABI_FREE(al)
735 :
736 120 : if(inplace==1)then
737 1560 : cprj=cprjout_
738 : else
739 0 : cprjout=cprjout_
740 : end if
741 120 : call pawcprj_free(cprjout_)
742 1080 : ABI_FREE(cprjout_)
743 :
744 : end if
745 :
746 282 : end subroutine lincom_cgcprj
747 : !!***
748 :
749 : !!****m* ABINIT/cgcprj_cholesky
750 : !! NAME
751 : !! cgcprj_cholesky
752 : !!
753 : !! FUNCTION
754 : !! Cholesky orthonormalization of the vectors stored in cg+cprj mode.
755 : !!
756 : !! This implementation is NOT band-parallelized
757 : !! Also, it is far of being optimal at the level of linear algebra
758 : !!
759 : !! INPUTS
760 : !! atindx1(natom)=index table for atoms, inverse of atindx
761 : !! dimcprj(natom)=number of lmn components in the <p_{lmn}^i|\psi> for the i-th atom
762 : !! icg=shift in cg array to locate current k-point and spinpol
763 : !! ikpt=current k point index
764 : !! isppol=current spin polarization index
765 : !! istwf=input option parameter that describes the storage of wfs
766 : !! mcg=second dimension of cg array (mpw*nspinor*mband*mkmem*nsppol)
767 : !! mcprj=second dimension of cprj_k array
768 : !! mkmem=number of k points which can fit in memory
769 : !! mpi_enreg=information about MPI parallelization
770 : !! natom=number of atoms
771 : !! nattyp(ntypat)=number of atoms of each type in cell.
772 : !! nband=number of bands
773 : !! npw=number of planewaves in basis at this k point
774 : !! nspinor=number of spinor components
775 : !! nsppol=number of spin polarizations
776 : !! ntypat=number of types of atoms
777 : !! pawtab(dtset%ntypat*dtset%usepaw) <type(pawtab_type)>=paw tabulated starting data
778 : !! usepaw=1 if PAW is activated
779 : !!
780 : !! SIDE EFFECTS
781 : !! cg(2,mcg)= plane wave wavefunction coefficients for the set of input wavefunctions (all k points and spinpol)
782 : !! cprj_k(natom,mcprj) <type(pawcprj_type)>= projected input wave functions <Proj_i|Cnk> with NL projectors for the specific k point and spinpol
783 : !!
784 : !! SOURCE
785 :
786 141 : subroutine cgcprj_cholesky(atindx1,cg,cprj_k,dimcprj,icg,ikpt,isppol,istwf,mcg,mcprj,mkmem,&
787 141 : & mpi_enreg,natom,nattyp,nband,npw,nspinor,nsppol,ntypat,pawtab,usepaw)
788 :
789 : !Arguments ------------------------------------
790 : !scalars
791 : integer,intent(in) :: icg,ikpt,isppol,istwf,mcg,mcprj,mkmem
792 : integer,intent(in) :: natom,nband,npw,nspinor,nsppol,ntypat,usepaw
793 : !arrays
794 : integer, intent(in) :: atindx1(natom),dimcprj(natom),nattyp(ntypat)
795 : real(dp), intent(inout) :: cg(2,mcg)
796 : type(pawcprj_type),intent(inout) :: cprj_k(natom,mcprj)
797 : type(MPI_type),intent(in) :: mpi_enreg
798 : type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw)
799 :
800 : !Local variables ------------------------------
801 : !scalars
802 : integer :: hermitian,ierr,ii,inplace
803 : !arrays
804 141 : real(dp), allocatable :: dmn(:,:,:),smn(:,:,:)
805 :
806 : ! *************************************************************************
807 :
808 564 : ABI_MALLOC(smn,(2,nband,nband))
809 423 : ABI_MALLOC(dmn,(2,nband,nband))
810 :
811 141 : hermitian=1
812 : call dotprod_set_cgcprj(atindx1,cg,cg,cprj_k,cprj_k,dimcprj,hermitian,&
813 : & 0,0,icg,icg,ikpt,isppol,istwf,nband,mcg,mcg,mcprj,mcprj,mkmem,&
814 141 : & mpi_enreg,natom,nattyp,nband,nband,npw,nspinor,nsppol,ntypat,pawtab,smn,usepaw)
815 :
816 : !Cholesky factorization: O = U^H U with U upper triangle matrix.
817 141 : call ZPOTRF('U',nband,smn,nband,ierr)
818 :
819 : !Solve X U = 1.
820 7233 : dmn=zero
821 690 : do ii=1,nband
822 690 : dmn(1,ii,ii)=one
823 : end do
824 141 : call ZTRSM('Right','Upper','Normal','Normal',nband,nband,cone,smn,nband,dmn,nband)
825 :
826 141 : inplace=1
827 : !This call does not take into account the fact that X=dmn is an upper triangular matrix...
828 : !The number of operations might be divided by two.
829 : call lincom_cgcprj(dmn,cg,cprj_k,dimcprj,&
830 141 : & icg,inplace,mcg,mcprj,natom,nband,nband,npw,nspinor,usepaw)
831 :
832 141 : ABI_FREE(smn)
833 141 : ABI_FREE(dmn)
834 :
835 141 : end subroutine cgcprj_cholesky
836 : !!***
837 :
838 : end module m_cgcprj
839 : !!***
|
