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