TABLE OF CONTENTS


ABINIT/cgcprj_cholesky [ Modules ]

[ Top ] [ Modules ]

NAME

  cgcprj_cholesky

FUNCTION

 Cholesky orthonormalization of the vectors stored in cg+cprj mode.

 This implementation is NOT band-parallelized
 Also, it is far of being optimal at the level of linear algebra

INPUTS

  atindx1(natom)=index table for atoms, inverse of atindx
  dimcprj(natom)=number of lmn components in the <p_{lmn}^i|\psi> for the i-th atom
  icg=shift in cg array to locate current k-point and spinpol
  ikpt=current k point index
  isppol=current spin polarization index
  istwf=input option parameter that describes the storage of wfs
  mcg=second dimension of cg array (mpw*nspinor*mband*mkmem*nsppol)
  mcprj=second dimension of cprj_k array
  mkmem=number of k points which can fit in memory
  mpi_enreg=information about MPI parallelization
  natom=number of atoms
  nattyp(ntypat)=number of atoms of each type in cell.
  nband=number of bands
  npw=number of planewaves in basis at this k point
  nspinor=number of spinor components
  nsppol=number of spin polarizations
  ntypat=number of types of atoms
  pawtab(dtset%ntypat*dtset%usepaw) <type(pawtab_type)>=paw tabulated starting data
  usepaw=1 if PAW is activated

SIDE EFFECTS

  cg(2,mcg)= plane wave wavefunction coefficients for the set of input wavefunctions (all k points and spinpol)
  cprj_k(natom,mcprj) <type(pawcprj_type)>= projected input wave functions <Proj_i|Cnk> with NL projectors for the specific k point and spinpol

SOURCE

786  subroutine cgcprj_cholesky(atindx1,cg,cprj_k,dimcprj,icg,ikpt,isppol,istwf,mcg,mcprj,mkmem,&
787 &  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  real(dp), allocatable :: dmn(:,:,:),smn(:,:,:)
805 
806 ! *************************************************************************
807 
808  ABI_MALLOC(smn,(2,nband,nband))
809  ABI_MALLOC(dmn,(2,nband,nband))
810 
811  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 & 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  call ZPOTRF('U',nband,smn,nband,ierr)
818 
819 !Solve X U = 1.
820  dmn=zero
821  do ii=1,nband
822    dmn(1,ii,ii)=one
823  end do
824  call ZTRSM('Right','Upper','Normal','Normal',nband,nband,cone,smn,nband,dmn,nband)
825 
826  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 & icg,inplace,mcg,mcprj,natom,nband,nband,npw,nspinor,usepaw)
831 
832  ABI_FREE(smn)
833  ABI_FREE(dmn)
834 
835 end subroutine cgcprj_cholesky

ABINIT/dotprod_set_cgcprj [ Functions ]

[ Top ] [ Functions ]

NAME

 dotprod_set_cgcprj

FUNCTION

 For one k point and spinpol, compute the matrix of scalar products between two sets of nband wavefunctions
 that are known in the cg+cprj representation
 Smn=<wf1m|wf2n>

 Can also treat the case of the computation of scalar products within one set of nband wavefunctions
 without recomputing already computed matrix elements (define hermitian=1)

 This implementation is NOT band-parallelized
 Also, it is far of being optimal at the level of linear algebra, and involves extra copying
 that are detrimental for performance...

INPUTS

  atindx1(natom)=index table for atoms, inverse of atindx
  cg1(2,mcg1)= plane wave wavefunction coefficients for the first set of wavefunctions (all k points and spinpol)
  cg2(2,mcg2)= plane wave wavefunction coefficients for the second set of wavefunctions (all k points and spinpol)
  cprj1(natom,mcprj1) <type(pawcprj_type)>= projected input wave functions <Proj_i|Cnk> with NL projectors in the first set,
  cprj2(natom,mcprj2) <type(pawcprj_type)>= projected input wave functions <Proj_i|Cnk> with NL projectors in the second set,
  dimcprj(natom)=number of lmn components in the <p_{lmn}^i|\psi> for the i-th atom
  hermitian= if 1, consider that the Smn matrix is hermitian, and do not recompute already computed matrix elements.
  ibg1=shift in cprj1 array to locate current k-point.
    Might be 0, in which case cprj1 is not copied internally, which saves some time/space
  ibg2=shift in cprj2 array to locate current k-point.
    Might be 0, in which case cprj2 is not copied internally, which saves some time/space
  icg1=shift in cg1 array to locate current k-point
  icg2=shift in cg2 array to locate current k-point
  ikpt=current k point index
  isppol=current spin polarization index
  istwf=input option parameter that describes the storage of wfs
  mcg1=second dimension of cg1 array (mpw*nspinor*mband1*mkmem*nsppol)
  mcg2=second dimension of cg2 array (mpw*nspinor*mband2*mkmem*nsppol)
  mcprj1=second dimension of cprj1 array
  mcprj2=second dimension of cprj2 array
  mkmem=number of k points which can fit in memory
  mpi_enreg=information about MPI parallelization
  natom=number of atoms
  nattyp(ntypat)=number of atoms of each type in cell.
  nbd1=number of bands for the first set of wavefunctions
  nbd2=number of bands for the second set of wavefunctions
  npw=number of planewaves in basis at this k point
  nspinor=number of spinor components
  nsppol=number of spin polarizations
  ntypat=number of types of atoms
  pawtab(dtset%ntypat*dtset%usepaw) <type(pawtab_type)>=paw tabulated starting data
  usepaw=1 if PAW is activated

OUTPUT

  smn(2,nbd1,nbd2)=matrix of scalar products between the first set of wavefunctions and the second set of wavefunctions

SIDE EFFECTS

SOURCE

104 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 & 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  real(dp),allocatable :: cwavef1(:,:),cwavef2(:,:),proj(:,:,:)
128  real(dp),allocatable :: eigval(:),eigvec(:,:,:),matrx(:,:),zhpev1(:,:),zhpev2(:)
129  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  if(hermitian==1)then
141    if(nbd1/=nbd2)then
142      ABI_ERROR(' With hermitian==1, nb1 and nb2 must be equal ')
143    end if
144  end if
145 
146  ABI_MALLOC(cwavef1,(2,npw*nspinor))
147  ABI_MALLOC(cwavef2,(2,npw*nspinor))
148  if(usepaw==1) then
149    ABI_MALLOC(cprj1_k,(natom,nspinor*nbd1))
150    ABI_MALLOC(cprj2_k,(natom,nspinor*nbd2))
151    iorder=0 ! There is no change of ordering of cprj when copying wavefunctions
152  end if
153  if(usepaw==1 .and. ibg1/=0) then
154    call pawcprj_alloc(cprj1_k,cprj1(1,1)%ncpgr,dimcprj)
155  end if
156  if(usepaw==1 .and. ibg2/=0) then
157    call pawcprj_alloc(cprj2_k,cprj1(1,1)%ncpgr,dimcprj)
158  end if
159 
160  icgb1=icg1
161  do ibd1=1,nbd1
162 
163 !  Extract wavefunction information
164    do ig=1,npw*nspinor
165      cwavef1(1,ig)=cg1(1,ig+icgb1)
166      cwavef1(2,ig)=cg1(2,ig+icgb1)
167    end do
168    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 &     mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
172    end if
173 
174    icgb2=icg2
175    max_nbd2=nbd2
176    if(hermitian==1)max_nbd2=ibd1
177    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      do ig=1,npw*nspinor
184        cwavef2(1,ig)=cg2(1,ig+icgb2)
185        cwavef2(2,ig)=cg2(2,ig+icgb2)
186      end do
187 
188      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 &       mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
192      end if
193 
194 !    Calculate Smn=<cg1|cg2>
195      call dotprod_g(dotr,doti,istwf,npw*nspinor,2,cwavef1,cwavef2,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
196 
197      if(usepaw==1) then
198        ia =0
199        do itypat=1,ntypat
200          if(ibg1/=0 .and. ibg2/=0)then
201            do iat=1+ia,nattyp(itypat)+ia
202              do ilmn1=1,pawtab(itypat)%lmn_size
203                do ilmn2=1,ilmn1
204                  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 &                 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 &                 cprj1_k(iat,ibd1)%cp(2,ilmn1)*cprj2_k(iat,ibd2)%cp(1,ilmn2))
209                end do
210                do ilmn2=ilmn1+1,pawtab(itypat)%lmn_size
211                  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 &                 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 &                 cprj1_k(iat,ibd1)%cp(2,ilmn1)*cprj2_k(iat,ibd2)%cp(1,ilmn2))
216                end do
217              end do
218            end do
219          else if(ibg1/=0 .and. ibg2==0)then
220            do iat=1+ia,nattyp(itypat)+ia
221              do ilmn1=1,pawtab(itypat)%lmn_size
222                do ilmn2=1,ilmn1
223                  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 &                 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 &                 cprj1_k(iat,ibd1)%cp(2,ilmn1)*cprj2(iat,ibd2)%cp(1,ilmn2))
228                end do
229                do ilmn2=ilmn1+1,pawtab(itypat)%lmn_size
230                  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 &                 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 &                 cprj1_k(iat,ibd1)%cp(2,ilmn1)*cprj2(iat,ibd2)%cp(1,ilmn2))
235                end do
236              end do
237            end do
238          else if(ibg1==0 .and. ibg2/=0)then
239            do iat=1+ia,nattyp(itypat)+ia
240              do ilmn1=1,pawtab(itypat)%lmn_size
241                do ilmn2=1,ilmn1
242                  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 &                 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 &                 cprj1(iat,ibd1)%cp(2,ilmn1)*cprj2_k(iat,ibd2)%cp(1,ilmn2))
247                end do
248                do ilmn2=ilmn1+1,pawtab(itypat)%lmn_size
249                  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 &                 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 &                 cprj1(iat,ibd1)%cp(2,ilmn1)*cprj2_k(iat,ibd2)%cp(1,ilmn2))
254                end do
255              end do
256            end do
257          else if(ibg1==0 .and. ibg2==0)then
258            do iat=1+ia,nattyp(itypat)+ia
259              do ilmn1=1,pawtab(itypat)%lmn_size
260                do ilmn2=1,ilmn1
261                  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 &                 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 &                 cprj1(iat,ibd1)%cp(2,ilmn1)*cprj2(iat,ibd2)%cp(1,ilmn2))
266                end do
267                do ilmn2=ilmn1+1,pawtab(itypat)%lmn_size
268                  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 &                 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 &                 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          ia=ia+nattyp(itypat)
278        end do
279      end if
280      smn(1,ibd1,ibd2)=dotr
281      smn(2,ibd1,ibd2)=doti
282 !    End loop over bands ibd2
283      icgb2=icgb2+npw*nspinor
284 
285    end do
286 
287 !  End loop over bands ibd1
288    icgb1=icgb1+npw*nspinor
289  end do
290 
291 !Complete the matrix if hermitian
292  if(hermitian==1)then
293    do ibd1=1,nbd1-1
294      do ibd2=ibd1+1,nbd2
295        smn(1,ibd1,ibd2)= smn(1,ibd2,ibd1)
296        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  ABI_FREE(cwavef1)
364  ABI_FREE(cwavef2)
365  if(usepaw==1) then
366    if(ibg1/=0)then
367      call pawcprj_free(cprj1_k)
368    end if
369    if(ibg2/=0)then
370      call pawcprj_free(cprj2_k)
371    end if
372    ABI_FREE(cprj1_k)
373    ABI_FREE(cprj2_k)
374  end if
375 
376 end subroutine dotprod_set_cgcprj

ABINIT/dotprodm_sumdiag_cgcprj [ Functions ]

[ Top ] [ Functions ]

NAME

 dotprodm_sumdiag_cgcprj

FUNCTION

 For one k point and spinpol, compute the matrix of sum of diagonal scalar products
 between sets of nband wavefunctions that are known in the cg+cprj representation
 The sets of wavefunctions must be contained in a big array of sets of wavefunctions.
 Sij=Sum_m <wfm(set i)|wfm(set j)>

 Can also treat the case of the computation of scalar products within one set of wavefunctions

 This implementation is NOT band-parallelized
 Also, it is far of being optimal at the level of linear algebra, and involves extra copying
 that are detrimental for performance...

INPUTS

  atindx1(natom)=index table for atoms, inverse of atindx
  cg_set(2,mcg,mset)= plane wave wavefunction coefficients for several sets of wavefunctions (all k points and spins)
  cprj_set(natom,mcprj,mset) <type(pawcprj_type)>= projected input wave functions <Proj_i|Cnk>
    with NL projectors in the different sets
  dimcprj(natom)=number of lmn components in the <p_{lmn}^i|\psi> for the i-th atom
  ibg=shift in cprj_set array to locate current k-point
  icg=shift in cg_set array to locate current k-point
  ikpt=current k point index
  isppol=current spin polarization index
  istwf=input option parameter that describes the storage of wfs
  mband=maximum number of bands (used in the dimensioning of cprj_set)
  mcg=second dimension of cg array (mpw*nspinor*mband*mkmem*nsppol)
  mcprj=second dimension of cprj array
  mkmem=number of k points which can fit in memory
  mpi_enreg=information about MPI parallelization
  mset=third dimension of cg_set and cprj_set, maximum number of sets
  natom=number of atoms
  nattyp(ntypat)=number of atoms of each type in cell.
  nbd=number of bands for each set of wavefunctions
  npw=number of planewaves in basis at this k point
  nset1=number of sets of wavefunctions to be considered in the left side of the scalar products
  nset2=number of sets of wavefunctions to be considered in the right side of the scalar products
  nspinor=number of spinor components
  nsppol=number of spin polarizations
  ntypat=number of types of atoms
  pawtab(dtset%ntypat*dtset%usepaw) <type(pawtab_type)>=paw tabulated starting data
  shift_set1=shift that defines the first set of wavefunctions to be considered in the left side of the scalar products
  shift_set2=shift that defines the first set of wavefunctions to be considered in the right side of the scalar products
  usepaw=1 if PAW is activated

OUTPUT

  smn(2,nset1,nset2)=matrix of sum of diagonal scalar products between the first set
    of wavefunctions and the second set of wavefunctions

SIDE EFFECTS

SOURCE

435 subroutine dotprodm_sumdiag_cgcprj(atindx1,cg_set,cprj_set,dimcprj,&
436 & ibg,icg,ikpt,isppol,istwf,mband,mcg,mcprj,mkmem,&
437 & mpi_enreg,mset,natom,nattyp,nbd,npw,nset1,nset2,nspinor,nsppol,ntypat,&
438 & 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  real(dp),allocatable :: cwavef1(:,:),cwavef2(:,:)
461  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  ABI_MALLOC(cwavef1,(2,npw*nspinor))
471  ABI_MALLOC(cwavef2,(2,npw*nspinor))
472  if(usepaw==1) then
473    ABI_MALLOC(cprj1_k,(natom,nspinor*nbd))
474    ABI_MALLOC(cprj2_k,(natom,nspinor*nbd))
475    iorder=0 ! There is no change of ordering in the copy of wavefunctions
476    call pawcprj_alloc(cprj1_k,cprj_set(1,1,1)%ncpgr,dimcprj)
477  end if
478 
479  smn(:,:,:)=zero
480 
481  icgb=icg
482  do ibd=1,nbd
483 
484    do iset1=1,nset1
485 
486      ind_set1=iset1+shift_set1
487 
488 !    Extract wavefunction information
489      do ig=1,npw*nspinor
490        cwavef1(1,ig)=cg_set(1,ig+icgb,ind_set1)
491        cwavef1(2,ig)=cg_set(2,ig+icgb,ind_set1)
492      end do
493      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 &       mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
497      end if
498 
499      do iset2=1,nset2
500 
501        ind_set2=iset2+shift_set2
502        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        else if(ind_set1/=ind_set2)then
506 
507 !        Extract wavefunction information
508          do ig=1,npw*nspinor
509            cwavef2(1,ig)=cg_set(1,ig+icgb,ind_set2)
510            cwavef2(2,ig)=cg_set(2,ig+icgb,ind_set2)
511          end do
512 
513          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 &           mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
517          end if
518 
519 !        Calculate Smn=<cg1|cg2>
520          call dotprod_g(dotr,doti,istwf,npw*nspinor,2,cwavef1,cwavef2,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
521 
522          if(usepaw==1) then
523            ia =0
524            do itypat=1,ntypat
525              do iat=1+ia,nattyp(itypat)+ia
526                do ilmn1=1,pawtab(itypat)%lmn_size
527                  do ilmn2=1,ilmn1
528                    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 &                   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 &                   cprj1_k(iat,ibd)%cp(2,ilmn1)*cprj2_k(iat,ibd)%cp(1,ilmn2))
533                  end do
534                  do ilmn2=ilmn1+1,pawtab(itypat)%lmn_size
535                    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 &                   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 &                   cprj1_k(iat,ibd)%cp(2,ilmn1)*cprj2_k(iat,ibd)%cp(1,ilmn2))
540                  end do
541                end do
542              end do
543              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          call dotprod_g(dotr,doti,istwf,npw*nspinor,1,cwavef1,cwavef1,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
553 
554          if(usepaw==1) then
555            ia =0
556            do itypat=1,ntypat
557              do iat=1+ia,nattyp(itypat)+ia
558                do ilmn1=1,pawtab(itypat)%lmn_size
559                  do ilmn2=1,ilmn1
560                    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 &                   cprj1_k(iat,ibd)%cp(2,ilmn1)*cprj1_k(iat,ibd)%cp(2,ilmn2))
563                  end do
564                  do ilmn2=ilmn1+1,pawtab(itypat)%lmn_size
565                    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 &                   cprj1_k(iat,ibd)%cp(2,ilmn1)*cprj1_k(iat,ibd)%cp(2,ilmn2))
568                  end do
569                end do
570              end do
571              ia=ia+nattyp(itypat)
572            end do
573          end if ! usepaw
574          doti=zero
575 
576        end if ! Compare ind_set1 and ind_set2
577 
578        smn(1,iset1,iset2)=smn(1,iset1,iset2)+dotr
579        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    icgb=icgb+npw*nspinor
586  end do
587 
588 !Complete the matrix, using its hermitian property.
589  do iset1=1,nset1
590    ind_set1=iset1+shift_set1
591    do iset2=1,nset2
592      ind_set2=iset2+shift_set2
593      if(ind_set2<ind_set1 .and. ind_set2>shift_set1)then
594        smn(1,iset1,iset2)= smn(1,iset2,iset1)
595        smn(2,iset1,iset2)=-smn(2,iset2,iset1)
596      end if
597    end do
598  end do
599 
600  ABI_FREE(cwavef1)
601  ABI_FREE(cwavef2)
602  if(usepaw==1) then
603    call pawcprj_free(cprj1_k)
604    call pawcprj_free(cprj2_k)
605    ABI_FREE(cprj1_k)
606    ABI_FREE(cprj2_k)
607  end if
608 
609 end subroutine dotprodm_sumdiag_cgcprj

ABINIT/lincom_cgcprj [ Functions ]

[ Top ] [ Functions ]

NAME

 lincom_cgcprj

FUNCTION

 For one k point and spin, compute a set (size nband_out) of linear combinations of nband_in wavefunctions,
 that are known in the cg+cprj representation :
 cgout_n(:,:) <--- Sum_m [ cg_m(:,:) . alpha_mn ]
 cprjout_n(:,:) <--- Sum_m [ cprj_m(:,:) . alpha_mn ]
 If nband_out is smaller or equal to nband_in, the result might be in-place
 output in cg instead of cgout, and in cprj instead of cprjout).
 Otherwise, it is contained in the optional cgout+cprjout pair.

 In the present status, the cg and cgout relates to all the k points and spins, and rely on the icg index,
 while it is assumed that cprj and cprjout refer to the specific k point and spin.
 This is not coherent.
 THIS MIGHT BE CHANGED IN THE FUTURE !

 This implementation is NOT band-parallelized
 Also, it is far of being optimal at the level of linear algebra, and involves extra copying
 that are detrimental for performance...

INPUTS

  alpha_mn(2,nband_in,nband_out)=complex matrix of coefficients of the linear combinations to be computed
  dimcprj(natom)=number of lmn components in the <p_{lmn}^i|\psi> for the i-th atom
  icg=shift in cg array to locate current k-point and spinpol (for input, and possibly for in-place output)
  inplace= if 0, output in cgout and cprjout ; if 1, output in cg and cprj
  mcg=second dimension of cg array (mpw*nspinor*mband*mkmem*nsppol)
  mcprj=second dimension of cprj array
  natom=number of atoms
  nband_in=number of bands, size of the input set of wavefunctions
  nband_out=number of bands, size of the output set of wavefunctions (should be equal to nband_in if inplace==1)
  npw=number of planewaves in basis at this k point
  nspinor=number of spinor components
  usepaw=1 if PAW is activated
  [icgout= shift in cgout array to locate current k-point and spinpol (for output)]
  [mcgout=second dimension of cgout array (mpw*nspinor*mband*mkmem*nsppol)]
  [mcprjout=second dimension of cprjout array]

OUTPUT

  [cgout(2,mcgout)= plane wave wavefunction coefficients for the set of output wavefunctions]
  [cprjout(natom,mcprjout) <type(pawcprj_type)>= projected output wave functions <Proj_i|Cnk> with NL projectors]

SIDE EFFECTS

  (this quantities are input, and possibly updated output when inplace==1)
  cg(2,mcg)= plane wave wavefunction coefficients for the set of input wavefunctions (all k points and spinpol)
  cprj(natom,mcprj) <type(pawcprj_type)>= projected input wave functions <Proj_i|Cnk> with NL projectors

SOURCE

663  subroutine lincom_cgcprj(alpha_mn,cg,cprj,dimcprj,&
664 & icg,inplace,mcg,mcprj,natom,nband_in,nband_out,npw,nspinor,usepaw, &
665 & 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  real(dp),allocatable :: al(:,:),cgout_(:,:)
685  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  if(inplace==0)then
696    if(.not.present(cgout))then
697      ABI_ERROR(' inplace==0 while .not.present(cgout) is not permitted ')
698    end if
699    if(usepaw==1) then
700      if(.not.present(cprjout))then
701        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  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 & alpha_mn,nband_in,dcmplx(0._dp),cgout_,npw*nspinor)
712 
713  if(inplace==1)then
714    cg(:,icg+1:icg+npw*nspinor*nband_out)=cgout_
715  else
716    cgout(:,icgout+1:icgout+npw*nspinor*nband_out)=cgout_
717  end if
718  ABI_FREE(cgout_)
719 
720 !Take care of the cprj part
721  if(usepaw==1) then
722 
723    ABI_MALLOC(cprjout_,(natom,nspinor*nband_out))
724    call pawcprj_alloc(cprjout_,cprj(1,1)%ncpgr,dimcprj)
725    ABI_MALLOC(al,(2,nband_in))
726    do iband_out=1,nband_out
727      ii=(iband_out-1)*nspinor
728      do iband_in=1,nband_in
729        al(1,iband_in)=alpha_mn(1,iband_in,iband_out)
730        al(2,iband_in)=alpha_mn(2,iband_in,iband_out)
731      end do
732      call pawcprj_lincom(al,cprj,cprjout_(:,ii+1:ii+nspinor),nband_in)
733    end do
734    ABI_FREE(al)
735 
736    if(inplace==1)then
737      cprj=cprjout_
738    else
739      cprjout=cprjout_
740    end if
741    call pawcprj_free(cprjout_)
742    ABI_FREE(cprjout_)
743 
744  end if
745 
746 end subroutine lincom_cgcprj

ABINIT/m_cgcprj [ Modules ]

[ Top ] [ Modules ]

NAME

  m_cgcprj

FUNCTION

  Functions operating on wavefunctions in the cg+cprj representation.

COPYRIGHT

  Copyright (C) 2008-2022 ABINIT group (XG)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

SOURCE

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