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

PARENTS

      wf_mixing

CHILDREN

      dotprod_set_cgcprj,lincom_cgcprj,zpotrf,ztrsm

SOURCE

842  subroutine cgcprj_cholesky(atindx1,cg,cprj_k,dimcprj,icg,ikpt,isppol,istwf,mcg,mcprj,mkmem,&
843 &  mpi_enreg,natom,nattyp,nband,npw,nspinor,nsppol,ntypat,pawtab,usepaw)
844 
845 
846 !This section has been created automatically by the script Abilint (TD).
847 !Do not modify the following lines by hand.
848 #undef ABI_FUNC
849 #define ABI_FUNC 'cgcprj_cholesky'
850 !End of the abilint section
851 
852  implicit none
853 
854 !Arguments ------------------------------------
855 !scalars
856  integer,intent(in) :: icg,ikpt,isppol,istwf,mcg,mcprj,mkmem
857  integer,intent(in) :: natom,nband,npw,nspinor,nsppol,ntypat,usepaw
858 !arrays
859  integer, intent(in) :: atindx1(natom),dimcprj(natom),nattyp(ntypat)
860  real(dp), intent(inout) :: cg(2,mcg)
861  type(pawcprj_type),intent(inout) :: cprj_k(natom,mcprj)
862  type(MPI_type),intent(in) :: mpi_enreg
863  type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw)
864 
865 !Local variables ------------------------------
866 !scalars
867  integer :: hermitian,ierr,ii,inplace
868 !arrays
869  real(dp), allocatable :: dmn(:,:,:),smn(:,:,:)
870 
871 ! *************************************************************************
872 
873  ABI_ALLOCATE(smn,(2,nband,nband))
874  ABI_ALLOCATE(dmn,(2,nband,nband))
875 
876  hermitian=1
877  call dotprod_set_cgcprj(atindx1,cg,cg,cprj_k,cprj_k,dimcprj,hermitian,&
878 & 0,0,icg,icg,ikpt,isppol,istwf,nband,mcg,mcg,mcprj,mcprj,mkmem,&
879 & mpi_enreg,natom,nattyp,nband,nband,npw,nspinor,nsppol,ntypat,pawtab,smn,usepaw)
880 
881 !Cholesky factorization: O = U^H U with U upper triangle matrix.
882  call ZPOTRF('U',nband,smn,nband,ierr)
883 
884 !Solve X U = 1.
885  dmn=zero
886  do ii=1,nband
887    dmn(1,ii,ii)=one
888  end do
889  call ZTRSM('Right','Upper','Normal','Normal',nband,nband,cone,smn,nband,dmn,nband)
890 
891  inplace=1
892 !This call does not take into account the fact that X=dmn is an upper triangular matrix...
893 !The number of operations might be divided by two.
894  call lincom_cgcprj(dmn,cg,cprj_k,dimcprj,&
895 & icg,inplace,mcg,mcprj,natom,nband,nband,npw,nspinor,usepaw)
896 
897  ABI_DEALLOCATE(smn)
898  ABI_DEALLOCATE(dmn)
899 
900 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

PARENTS

      cgcprj_cholesky,wf_mixing

CHILDREN

      dotprod_g,pawcprj_alloc,pawcprj_free,pawcprj_get,zhpev

SOURCE

115 subroutine dotprod_set_cgcprj(atindx1,cg1,cg2,cprj1,cprj2,dimcprj,hermitian,&
116 & ibg1,ibg2,icg1,icg2,ikpt,isppol,istwf,mband,mcg1,mcg2,mcprj1,mcprj2,mkmem,&
117 & mpi_enreg,natom,nattyp,nbd1,nbd2,npw,nspinor,nsppol,ntypat,pawtab,smn,usepaw)
118 
119 
120 !This section has been created automatically by the script Abilint (TD).
121 !Do not modify the following lines by hand.
122 #undef ABI_FUNC
123 #define ABI_FUNC 'dotprod_set_cgcprj'
124 !End of the abilint section
125 
126  implicit none
127 
128 !Arguments ------------------------------------
129 !scalars
130  integer, intent(in) :: hermitian,ibg1,ibg2,icg1,icg2,ikpt,isppol,istwf
131  integer, intent(in) :: mkmem,mband,mcg1,mcg2,mcprj1,mcprj2
132  integer, intent(in) :: natom,nbd1,nbd2,npw,nspinor,nsppol,ntypat,usepaw
133  type(MPI_type),intent(in) :: mpi_enreg
134 !arrays
135  integer, intent(in) :: atindx1(natom),dimcprj(natom),nattyp(ntypat)
136  real(dp), intent(in) :: cg1(2,mcg1),cg2(2,mcg2)
137  real(dp), intent(out) :: smn(2,nbd1,nbd2)
138  type(pawcprj_type),intent(in) :: cprj1(natom,mcprj1),cprj2(natom,mcprj2)
139  type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw)
140 
141 !Local variables-------------------------------
142 !scalars
143  integer :: ia,iat,itypat,ibd1,ibd2,icgb1,icgb2,ier,ig,ii,i1,i2,iorder
144  integer :: ilmn1,ilmn2,klmn,max_nbd2,nbd
145  real(dp) :: dotr,doti
146 !arrays
147  real(dp),allocatable :: cwavef1(:,:),cwavef2(:,:),proj(:,:,:)
148  real(dp),allocatable :: eigval(:),eigvec(:,:,:),matrx(:,:),zhpev1(:,:),zhpev2(:)
149  type(pawcprj_type),allocatable :: cprj1_k(:,:),cprj2_k(:,:)
150 
151 ! *************************************************************************
152 
153 !DEBUG
154 !write(std_out,*)' dotprod_set_cgcprj : enter '
155 !write(std_out,*)' dotprod_set_cgcprj : npw, nspinor=',npw,nspinor
156 !write(std_out,*)' dotprod_set_cgcprj : usepaw,nbd1,nbd2=',usepaw,nbd1,nbd2
157 !call flush(std_out)
158 !ENDDEBUG
159 
160  if(hermitian==1)then
161    if(nbd1/=nbd2)then
162      MSG_ERROR(' With hermitian==1, nb1 and nb2 must be equal ')
163    end if
164  end if
165 
166  ABI_ALLOCATE(cwavef1,(2,npw*nspinor))
167  ABI_ALLOCATE(cwavef2,(2,npw*nspinor))
168  if(usepaw==1) then
169    ABI_DATATYPE_ALLOCATE(cprj1_k,(natom,nspinor*nbd1))
170    ABI_DATATYPE_ALLOCATE(cprj2_k,(natom,nspinor*nbd2))
171    iorder=0 ! There is no change of ordering of cprj when copying wavefunctions
172  end if
173  if(usepaw==1 .and. ibg1/=0) then
174    call pawcprj_alloc(cprj1_k,cprj1(1,1)%ncpgr,dimcprj)
175  end if
176  if(usepaw==1 .and. ibg2/=0) then
177    call pawcprj_alloc(cprj2_k,cprj1(1,1)%ncpgr,dimcprj)
178  end if
179 
180  icgb1=icg1
181  do ibd1=1,nbd1
182 
183 !  Extract wavefunction information
184    do ig=1,npw*nspinor
185      cwavef1(1,ig)=cg1(1,ig+icgb1)
186      cwavef1(2,ig)=cg1(2,ig+icgb1)
187    end do
188    if(usepaw==1 .and. ibg1/=0) then
189      call pawcprj_get(atindx1,cprj1_k,cprj1,natom,1,ibg1,ikpt,iorder,isppol,mband,&
190 &     mkmem,natom,nbd1,nbd1,nspinor,nsppol,0,&
191 &     mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
192    end if
193 
194    icgb2=icg2
195    max_nbd2=nbd2
196    if(hermitian==1)max_nbd2=ibd1
197    do ibd2=1,max_nbd2
198 
199 !    XG171222 Note that this copy step, being inside the ibd1 loop, is quite detrimental.
200 !    It might be reduced by copying several cwavef2, and use a ZGEMM type of approach.
201 
202 !    Extract wavefunction information
203      do ig=1,npw*nspinor
204        cwavef2(1,ig)=cg2(1,ig+icgb2)
205        cwavef2(2,ig)=cg2(2,ig+icgb2)
206      end do
207 
208      if(usepaw==1 .and. ibg2/=0) then
209        call pawcprj_get(atindx1,cprj2_k,cprj2,natom,1,ibg2,ikpt,iorder,isppol,mband,&
210 &       mkmem,natom,nbd2,nbd2,nspinor,nsppol,0,&
211 &       mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
212      end if
213 
214 !    Calculate Smn=<cg1|cg2>
215      call dotprod_g(dotr,doti,istwf,npw*nspinor,2,cwavef1,cwavef2,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
216 
217      if(usepaw==1) then
218        ia =0
219        do itypat=1,ntypat
220          if(ibg1/=0 .and. ibg2/=0)then
221            do iat=1+ia,nattyp(itypat)+ia
222              do ilmn1=1,pawtab(itypat)%lmn_size
223                do ilmn2=1,ilmn1
224                  klmn=((ilmn1-1)*ilmn1)/2+ilmn2
225                  dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd1)%cp(1,ilmn1)*cprj2_k(iat,ibd2)%cp(1,ilmn2)+&
226 &                 cprj1_k(iat,ibd1)%cp(2,ilmn1)*cprj2_k(iat,ibd2)%cp(2,ilmn2))
227                  doti=doti+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd1)%cp(1,ilmn1)*cprj2_k(iat,ibd2)%cp(2,ilmn2)-&
228 &                 cprj1_k(iat,ibd1)%cp(2,ilmn1)*cprj2_k(iat,ibd2)%cp(1,ilmn2))
229                end do
230                do ilmn2=ilmn1+1,pawtab(itypat)%lmn_size
231                  klmn=((ilmn2-1)*ilmn2)/2+ilmn1
232                  dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd1)%cp(1,ilmn1)*cprj2_k(iat,ibd2)%cp(1,ilmn2)+&
233 &                 cprj1_k(iat,ibd1)%cp(2,ilmn1)*cprj2_k(iat,ibd2)%cp(2,ilmn2))
234                  doti=doti+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd1)%cp(1,ilmn1)*cprj2_k(iat,ibd2)%cp(2,ilmn2)-&
235 &                 cprj1_k(iat,ibd1)%cp(2,ilmn1)*cprj2_k(iat,ibd2)%cp(1,ilmn2))
236                end do
237              end do
238            end do
239          else if(ibg1/=0 .and. ibg2==0)then
240            do iat=1+ia,nattyp(itypat)+ia
241              do ilmn1=1,pawtab(itypat)%lmn_size
242                do ilmn2=1,ilmn1
243                  klmn=((ilmn1-1)*ilmn1)/2+ilmn2
244                  dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd1)%cp(1,ilmn1)*cprj2(iat,ibd2)%cp(1,ilmn2)+&
245 &                 cprj1_k(iat,ibd1)%cp(2,ilmn1)*cprj2(iat,ibd2)%cp(2,ilmn2))
246                  doti=doti+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd1)%cp(1,ilmn1)*cprj2(iat,ibd2)%cp(2,ilmn2)-&
247 &                 cprj1_k(iat,ibd1)%cp(2,ilmn1)*cprj2(iat,ibd2)%cp(1,ilmn2))
248                end do
249                do ilmn2=ilmn1+1,pawtab(itypat)%lmn_size
250                  klmn=((ilmn2-1)*ilmn2)/2+ilmn1
251                  dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd1)%cp(1,ilmn1)*cprj2(iat,ibd2)%cp(1,ilmn2)+&
252 &                 cprj1_k(iat,ibd1)%cp(2,ilmn1)*cprj2(iat,ibd2)%cp(2,ilmn2))
253                  doti=doti+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd1)%cp(1,ilmn1)*cprj2(iat,ibd2)%cp(2,ilmn2)-&
254 &                 cprj1_k(iat,ibd1)%cp(2,ilmn1)*cprj2(iat,ibd2)%cp(1,ilmn2))
255                end do
256              end do
257            end do
258          else if(ibg1==0 .and. ibg2/=0)then
259            do iat=1+ia,nattyp(itypat)+ia
260              do ilmn1=1,pawtab(itypat)%lmn_size
261                do ilmn2=1,ilmn1
262                  klmn=((ilmn1-1)*ilmn1)/2+ilmn2
263                  dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj1(iat,ibd1)%cp(1,ilmn1)*cprj2_k(iat,ibd2)%cp(1,ilmn2)+&
264 &                 cprj1(iat,ibd1)%cp(2,ilmn1)*cprj2_k(iat,ibd2)%cp(2,ilmn2))
265                  doti=doti+pawtab(itypat)%sij(klmn)*(cprj1(iat,ibd1)%cp(1,ilmn1)*cprj2_k(iat,ibd2)%cp(2,ilmn2)-&
266 &                 cprj1(iat,ibd1)%cp(2,ilmn1)*cprj2_k(iat,ibd2)%cp(1,ilmn2))
267                end do
268                do ilmn2=ilmn1+1,pawtab(itypat)%lmn_size
269                  klmn=((ilmn2-1)*ilmn2)/2+ilmn1
270                  dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj1(iat,ibd1)%cp(1,ilmn1)*cprj2_k(iat,ibd2)%cp(1,ilmn2)+&
271 &                 cprj1(iat,ibd1)%cp(2,ilmn1)*cprj2_k(iat,ibd2)%cp(2,ilmn2))
272                  doti=doti+pawtab(itypat)%sij(klmn)*(cprj1(iat,ibd1)%cp(1,ilmn1)*cprj2_k(iat,ibd2)%cp(2,ilmn2)-&
273 &                 cprj1(iat,ibd1)%cp(2,ilmn1)*cprj2_k(iat,ibd2)%cp(1,ilmn2))
274                end do
275              end do
276            end do
277          else if(ibg1==0 .and. ibg2==0)then
278            do iat=1+ia,nattyp(itypat)+ia
279              do ilmn1=1,pawtab(itypat)%lmn_size
280                do ilmn2=1,ilmn1
281                  klmn=((ilmn1-1)*ilmn1)/2+ilmn2
282                  dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj1(iat,ibd1)%cp(1,ilmn1)*cprj2(iat,ibd2)%cp(1,ilmn2)+&
283 &                 cprj1(iat,ibd1)%cp(2,ilmn1)*cprj2(iat,ibd2)%cp(2,ilmn2))
284                  doti=doti+pawtab(itypat)%sij(klmn)*(cprj1(iat,ibd1)%cp(1,ilmn1)*cprj2(iat,ibd2)%cp(2,ilmn2)-&
285 &                 cprj1(iat,ibd1)%cp(2,ilmn1)*cprj2(iat,ibd2)%cp(1,ilmn2))
286                end do
287                do ilmn2=ilmn1+1,pawtab(itypat)%lmn_size
288                  klmn=((ilmn2-1)*ilmn2)/2+ilmn1
289                  dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj1(iat,ibd1)%cp(1,ilmn1)*cprj2(iat,ibd2)%cp(1,ilmn2)+&
290 &                 cprj1(iat,ibd1)%cp(2,ilmn1)*cprj2(iat,ibd2)%cp(2,ilmn2))
291                  doti=doti+pawtab(itypat)%sij(klmn)*(cprj1(iat,ibd1)%cp(1,ilmn1)*cprj2(iat,ibd2)%cp(2,ilmn2)-&
292 &                 cprj1(iat,ibd1)%cp(2,ilmn1)*cprj2(iat,ibd2)%cp(1,ilmn2))
293                end do
294              end do
295            end do
296          end if
297          ia=ia+nattyp(itypat)
298        end do
299      end if
300      smn(1,ibd1,ibd2)=dotr
301      smn(2,ibd1,ibd2)=doti
302 !    End loop over bands ibd2
303      icgb2=icgb2+npw*nspinor
304 
305    end do
306 
307 !  End loop over bands ibd1
308    icgb1=icgb1+npw*nspinor
309  end do
310 
311 !Complete the matrix if hermitian
312  if(hermitian==1)then
313    do ibd1=1,nbd1-1
314      do ibd2=ibd1+1,nbd2
315        smn(1,ibd1,ibd2)= smn(1,ibd2,ibd1)
316        smn(2,ibd1,ibd2)=-smn(2,ibd2,ibd1)
317      end do
318    end do
319  end if
320 
321 !DEBUG
322 !write(std_out,*)' smn=',smn
323 !ENDDEBUG
324 
325 !====== Debugging section ==========
326  if(.false.)then
327 !DEBUG
328 !Compute the eigenvalues of the projector S herm(S) or herm(S) S, depending on which has lowest dimension.
329 !write(std_out,*)' dotprod_set_cgcprj : compute the projector matrix '
330    nbd=min(nbd1,nbd2)
331    ABI_ALLOCATE(proj,(2,nbd,nbd))
332    proj(:,:,:)=zero
333    if(nbd1<=nbd2)then
334      do ibd1=1,nbd1
335        do ibd2=1,nbd2
336          proj(1,:,ibd1)=proj(1,:,ibd1)+smn(1,:,ibd2)*smn(1,ibd1,ibd2)+smn(2,:,ibd2)*smn(2,ibd1,ibd2)
337          proj(2,:,ibd1)=proj(2,:,ibd1)-smn(1,:,ibd2)*smn(2,ibd1,ibd2)+smn(2,:,ibd2)*smn(1,ibd1,ibd2)
338        end do
339      end do
340    else
341      do ibd2=1,nbd2
342        do ibd1=1,nbd1
343          proj(1,:,ibd2)=proj(1,:,ibd2)+smn(1,ibd1,:)*smn(1,ibd1,ibd2)+smn(2,ibd1,:)*smn(2,ibd1,ibd2)
344          proj(2,:,ibd2)=proj(2,:,ibd2)+smn(1,ibd1,:)*smn(2,ibd1,ibd2)-smn(2,ibd1,:)*smn(1,ibd1,ibd2)
345        end do
346      end do
347    end if
348 
349 !write(std_out,*)' proj=',proj
350 
351 !write(std_out,*)' dotprod_set_cgcprj : compute the eigenvalues of the projector '
352    ABI_ALLOCATE(matrx,(2,(nbd*(nbd+1))/2))
353    ii=1
354    do i2=1,nbd
355      do i1=1,i2
356        matrx(1,ii)=proj(1,i1,i2)
357        matrx(2,ii)=proj(2,i1,i2)
358        ii=ii+1
359      end do
360    end do
361 
362    ABI_ALLOCATE(zhpev1,(2,2*nbd-1))
363    ABI_ALLOCATE(zhpev2,(3*nbd-2))
364    ABI_ALLOCATE(eigval,(nbd))
365    ABI_ALLOCATE(eigvec,(2,nbd,nbd))
366 
367    call ZHPEV ('V','U',nbd,matrx,eigval,eigvec,nbd,zhpev1,zhpev2,ier)
368 
369 !write(std_out,*)' eigval=',eigval
370 
371    ABI_DEALLOCATE(matrx)
372    ABI_DEALLOCATE(zhpev1)
373    ABI_DEALLOCATE(zhpev2)
374    ABI_DEALLOCATE(eigval)
375    ABI_DEALLOCATE(eigvec)
376 
377    ABI_DEALLOCATE(proj)
378 !stop
379 !ENDDEBUG
380  end if
381 !====== End of debugging section ==========
382 
383  ABI_DEALLOCATE(cwavef1)
384  ABI_DEALLOCATE(cwavef2)
385  if(usepaw==1) then
386    if(ibg1/=0)then
387      call pawcprj_free(cprj1_k)
388    end if
389    if(ibg2/=0)then
390      call pawcprj_free(cprj2_k)
391    end if
392    ABI_DATATYPE_DEALLOCATE(cprj1_k)
393    ABI_DATATYPE_DEALLOCATE(cprj2_k)
394  end if
395 
396 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

PARENTS

      wf_mixing

CHILDREN

      dotprod_g,pawcprj_alloc,pawcprj_free,pawcprj_get

SOURCE

461 subroutine dotprodm_sumdiag_cgcprj(atindx1,cg_set,cprj_set,dimcprj,&
462 & ibg,icg,ikpt,isppol,istwf,mband,mcg,mcprj,mkmem,&
463 & mpi_enreg,mset,natom,nattyp,nbd,npw,nset1,nset2,nspinor,nsppol,ntypat,&
464 & shift_set1,shift_set2,pawtab,smn,usepaw)
465 
466 
467 !This section has been created automatically by the script Abilint (TD).
468 !Do not modify the following lines by hand.
469 #undef ABI_FUNC
470 #define ABI_FUNC 'dotprodm_sumdiag_cgcprj'
471 !End of the abilint section
472 
473  implicit none
474 
475 !Arguments ------------------------------------
476 !scalars
477  integer, intent(in) :: ibg,icg,ikpt,isppol,istwf
478  integer, intent(in) :: mband,mcg,mcprj,mkmem,mset
479  integer, intent(in) :: natom,nbd,npw,nset1,nset2,nspinor,nsppol,ntypat
480  integer, intent(in) :: shift_set1,shift_set2,usepaw
481  type(MPI_type),intent(in) :: mpi_enreg
482 !arrays
483  integer, intent(in) :: atindx1(natom),dimcprj(natom),nattyp(ntypat)
484  real(dp), intent(in) :: cg_set(2,mcg,mset)
485  real(dp), intent(out) :: smn(2,nset1,nset2)
486  type(pawcprj_type),intent(in) :: cprj_set(natom,mcprj,mset)
487  type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw)
488 
489 !Local variables-------------------------------
490 !scalars
491  integer :: ia,iat,itypat,ibd,icgb,ig,iorder
492  integer :: ilmn1,ilmn2,ind_set1,ind_set2,iset1,iset2,klmn
493  real(dp) :: dotr,doti
494 !arrays
495  real(dp),allocatable :: cwavef1(:,:),cwavef2(:,:)
496  type(pawcprj_type),allocatable :: cprj1_k(:,:),cprj2_k(:,:)
497 
498 ! *************************************************************************
499 
500 !DEBUG
501 !write(std_out,*)' dotprodm_sumdiag_cgcprj : enter '
502 !call flush(std_out)
503 !ENDDEBUG
504 
505  ABI_ALLOCATE(cwavef1,(2,npw*nspinor))
506  ABI_ALLOCATE(cwavef2,(2,npw*nspinor))
507  if(usepaw==1) then
508    ABI_DATATYPE_ALLOCATE(cprj1_k,(natom,nspinor*nbd))
509    ABI_DATATYPE_ALLOCATE(cprj2_k,(natom,nspinor*nbd))
510    iorder=0 ! There is no change of ordering in the copy of wavefunctions
511    call pawcprj_alloc(cprj1_k,cprj_set(1,1,1)%ncpgr,dimcprj)
512  end if
513 
514  smn(:,:,:)=zero
515 
516  icgb=icg
517  do ibd=1,nbd
518 
519    do iset1=1,nset1
520 
521      ind_set1=iset1+shift_set1
522 
523 !    Extract wavefunction information
524      do ig=1,npw*nspinor
525        cwavef1(1,ig)=cg_set(1,ig+icgb,ind_set1)
526        cwavef1(2,ig)=cg_set(2,ig+icgb,ind_set1)
527      end do
528      if(usepaw==1) then
529        call pawcprj_get(atindx1,cprj1_k,cprj_set(:,:,ind_set1),natom,1,ibg,ikpt,iorder,isppol,mband,&
530 &       mkmem,natom,nbd,nbd,nspinor,nsppol,0,&
531 &       mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
532      end if
533 
534      do iset2=1,nset2
535 
536        ind_set2=iset2+shift_set2
537        if(ind_set2<ind_set1 .and. ind_set2>shift_set1)then
538          continue ! These matrix elements have already been computed, the smn matrix will be completed later.
539 
540        else if(ind_set1/=ind_set2)then
541 
542 !        Extract wavefunction information
543          do ig=1,npw*nspinor
544            cwavef2(1,ig)=cg_set(1,ig+icgb,ind_set2)
545            cwavef2(2,ig)=cg_set(2,ig+icgb,ind_set2)
546          end do
547 
548          if(usepaw==1) then
549            call pawcprj_get(atindx1,cprj2_k,cprj_set(:,:,ind_set2),natom,1,ibg,ikpt,iorder,isppol,mband,&
550 &           mkmem,natom,nbd,nbd,nspinor,nsppol,0,&
551 &           mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
552          end if
553 
554 !        Calculate Smn=<cg1|cg2>
555          call dotprod_g(dotr,doti,istwf,npw*nspinor,2,cwavef1,cwavef2,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
556 
557          if(usepaw==1) then
558            ia =0
559            do itypat=1,ntypat
560              do iat=1+ia,nattyp(itypat)+ia
561                do ilmn1=1,pawtab(itypat)%lmn_size
562                  do ilmn2=1,ilmn1
563                    klmn=((ilmn1-1)*ilmn1)/2+ilmn2
564                    dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd)%cp(1,ilmn1)*cprj2_k(iat,ibd)%cp(1,ilmn2)+&
565 &                   cprj1_k(iat,ibd)%cp(2,ilmn1)*cprj2_k(iat,ibd)%cp(2,ilmn2))
566                    doti=doti+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd)%cp(1,ilmn1)*cprj2_k(iat,ibd)%cp(2,ilmn2)-&
567 &                   cprj1_k(iat,ibd)%cp(2,ilmn1)*cprj2_k(iat,ibd)%cp(1,ilmn2))
568                  end do
569                  do ilmn2=ilmn1+1,pawtab(itypat)%lmn_size
570                    klmn=((ilmn2-1)*ilmn2)/2+ilmn1
571                    dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd)%cp(1,ilmn1)*cprj2_k(iat,ibd)%cp(1,ilmn2)+&
572 &                   cprj1_k(iat,ibd)%cp(2,ilmn1)*cprj2_k(iat,ibd)%cp(2,ilmn2))
573                    doti=doti+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd)%cp(1,ilmn1)*cprj2_k(iat,ibd)%cp(2,ilmn2)-&
574 &                   cprj1_k(iat,ibd)%cp(2,ilmn1)*cprj2_k(iat,ibd)%cp(1,ilmn2))
575                  end do
576                end do
577              end do
578              ia=ia+nattyp(itypat)
579            end do
580          end if ! usepaw
581 
582 !      if(.false.)then
583        else
584 !        Diagonal part : no need to extract another wavefunction, and the scalar product must be real
585 
586 !        Calculate Smn=<cg1|cg1>
587          call dotprod_g(dotr,doti,istwf,npw*nspinor,1,cwavef1,cwavef1,mpi_enreg%me_g0,mpi_enreg%comm_spinorfft)
588 
589          if(usepaw==1) then
590            ia =0
591            do itypat=1,ntypat
592              do iat=1+ia,nattyp(itypat)+ia
593                do ilmn1=1,pawtab(itypat)%lmn_size
594                  do ilmn2=1,ilmn1
595                    klmn=((ilmn1-1)*ilmn1)/2+ilmn2
596                    dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd)%cp(1,ilmn1)*cprj1_k(iat,ibd)%cp(1,ilmn2)+&
597 &                   cprj1_k(iat,ibd)%cp(2,ilmn1)*cprj1_k(iat,ibd)%cp(2,ilmn2))
598                  end do
599                  do ilmn2=ilmn1+1,pawtab(itypat)%lmn_size
600                    klmn=((ilmn2-1)*ilmn2)/2+ilmn1
601                    dotr=dotr+pawtab(itypat)%sij(klmn)*(cprj1_k(iat,ibd)%cp(1,ilmn1)*cprj1_k(iat,ibd)%cp(1,ilmn2)+&
602 &                   cprj1_k(iat,ibd)%cp(2,ilmn1)*cprj1_k(iat,ibd)%cp(2,ilmn2))
603                  end do
604                end do
605              end do
606              ia=ia+nattyp(itypat)
607            end do
608          end if ! usepaw
609          doti=zero
610 
611        end if ! Compare ind_set1 and ind_set2
612 
613        smn(1,iset1,iset2)=smn(1,iset1,iset2)+dotr
614        smn(2,iset1,iset2)=smn(2,iset1,iset2)+doti
615 
616      end do ! iset2
617    end do ! iset1
618 
619 !  End loop over bands ibd
620    icgb=icgb+npw*nspinor
621  end do
622 
623 !Complete the matrix, using its hermitian property.
624  do iset1=1,nset1
625    ind_set1=iset1+shift_set1
626    do iset2=1,nset2
627      ind_set2=iset2+shift_set2
628      if(ind_set2<ind_set1 .and. ind_set2>shift_set1)then
629        smn(1,iset1,iset2)= smn(1,iset2,iset1)
630        smn(2,iset1,iset2)=-smn(2,iset2,iset1)
631      end if
632    end do
633  end do
634 
635  ABI_DEALLOCATE(cwavef1)
636  ABI_DEALLOCATE(cwavef2)
637  if(usepaw==1) then
638    call pawcprj_free(cprj1_k)
639    call pawcprj_free(cprj2_k)
640    ABI_DATATYPE_DEALLOCATE(cprj1_k)
641    ABI_DATATYPE_DEALLOCATE(cprj2_k)
642  end if
643 
644 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

PARENTS

      cgcprj_cholesky,wf_mixing

CHILDREN

      pawcprj_alloc,pawcprj_free,pawcprj_lincom,zgemm

SOURCE

704  subroutine lincom_cgcprj(alpha_mn,cg,cprj,dimcprj,&
705 & icg,inplace,mcg,mcprj,natom,nband_in,nband_out,npw,nspinor,usepaw, &
706 & cgout,cprjout,icgout) ! optional args
707 
708 
709 !This section has been created automatically by the script Abilint (TD).
710 !Do not modify the following lines by hand.
711 #undef ABI_FUNC
712 #define ABI_FUNC 'lincom_cgcprj'
713 !End of the abilint section
714 
715  implicit none
716 
717 !Arguments ------------------------------------
718 !scalars
719  integer, intent(in) :: icg,inplace,mcg,mcprj
720  integer, intent(in) :: natom,nband_in,nband_out,npw,nspinor,usepaw
721  integer, intent(in),optional :: icgout
722 !arrays
723  integer, intent(in) :: dimcprj(natom)
724  real(dp), intent(inout) :: cg(2,mcg)
725  real(dp), intent(in) :: alpha_mn(2,nband_in,nband_out)
726  real(dp), intent(out),optional :: cgout(:,:)
727  type(pawcprj_type),intent(inout) :: cprj(natom,mcprj)
728  type(pawcprj_type),intent(out),optional :: cprjout(:,:)
729 
730 !Local variables-------------------------------
731 !scalars
732  integer :: iband_in,iband_out,ii
733 !arrays
734  real(dp),allocatable :: al(:,:),cgout_(:,:)
735  type(pawcprj_type),allocatable :: cprjout_(:,:)
736 
737 ! *************************************************************************
738 
739 !DEBUG
740 !write(std_out,*)' lincom_cgcprj : enter '
741 !write(std_out,*)' lincom_cgcprj : npw, nspinor=',npw,nspinor
742 !write(std_out,*)' lincom_cgcprj : icgout=',icgout
743 !ENDDEBUG
744 
745  if(inplace==0)then
746    if(.not.present(cgout))then
747      MSG_ERROR(' inplace==0 while .not.present(cgout) is not permitted ')
748    end if
749    if(usepaw==1) then
750      if(.not.present(cprjout))then
751        MSG_ERROR(' inplace==0 and usepaw==1 while .not.present(cprjout) is not permitted ')
752      end if
753    end if
754  end if
755 
756 !Take care of the plane wave part
757  ABI_ALLOCATE(cgout_,(2,npw*nspinor*nband_out))
758 
759  call zgemm('N','N',npw*nspinor,nband_out,nband_in,dcmplx(1._dp), &
760 & cg(:,icg+1:icg+npw*nspinor*nband_in),npw*nspinor, &
761 & alpha_mn,nband_in,dcmplx(0._dp),cgout_,npw*nspinor)
762 
763  if(inplace==1)then
764    cg(:,icg+1:icg+npw*nspinor*nband_out)=cgout_
765  else
766    cgout(:,icgout+1:icgout+npw*nspinor*nband_out)=cgout_
767  end if
768  ABI_DEALLOCATE(cgout_)
769 
770 !Take care of the cprj part
771  if(usepaw==1) then
772 
773    ABI_DATATYPE_ALLOCATE(cprjout_,(natom,nspinor*nband_out))
774    call pawcprj_alloc(cprjout_,cprj(1,1)%ncpgr,dimcprj)
775    ABI_ALLOCATE(al,(2,nband_in))
776    do iband_out=1,nband_out
777      ii=(iband_out-1)*nspinor
778      do iband_in=1,nband_in
779        al(1,iband_in)=alpha_mn(1,iband_in,iband_out)
780        al(2,iband_in)=alpha_mn(2,iband_in,iband_out)
781      end do
782      call pawcprj_lincom(al,cprj,cprjout_(:,ii+1:ii+nspinor),nband_in)
783    end do
784    ABI_DEALLOCATE(al)
785 
786    if(inplace==1)then
787      cprj=cprjout_
788    else
789      cprjout=cprjout_
790    end if
791    call pawcprj_free(cprjout_)
792    ABI_DATATYPE_DEALLOCATE(cprjout_)
793 
794  end if
795 
796 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-2018 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 .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_cgcprj
28 
29  use defs_basis
30  use defs_abitypes
31  use m_cgtools
32  use m_errors
33  use m_xmpi
34 
35  use m_pawtab, only : pawtab_type
36  use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_get, pawcprj_free, pawcprj_lincom
37 
38  implicit none
39 
40  private