TABLE OF CONTENTS


ABINIT/m_cgtools [ Modules ]

[ Top ] [ Modules ]

NAME

  m_cgtools

FUNCTION

 This module defines wrappers for BLAS routines. The arguments are stored
 using the "cg" convention, namely real array of shape cg(2,...)

COPYRIGHT

 Copyright (C) 1992-2018 ABINIT group (MG, MT, XG, DCA, GZ, FB, MVer, DCA, GMR, FF)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

NOTES

 1) The convention about names of interfaced routine is: cg_<name>,
    where <name> is equal to the name of the standard BLAS routine

 2) Blas routines are called without an explicit interface on purpose since

    a) The compiler should pass the base address of the array to the F77 BLAS

    b) Any compiler would complain about type mismatch (REAL,COMPLEX)
       if an explicit interface is given.

 3) The use of mpi_type is not allowed here. MPI parallelism should be handled in a generic
    way by passing the MPI communicator so that the caller can decide how to handle MPI.


ABINIT/subdiago [ Functions ]

[ Top ] [ Functions ]

NAME

 subdiago

FUNCTION

 This routine diagonalizes the Hamiltonian in the eigenfunction subspace

INPUTS

  icg=shift to be applied on the location of data in the array cg
  igsc=shift to be applied on the location of data in the array gsc
  istwf_k=input parameter that describes the storage of wfs
  mcg=second dimension of the cg array
  mgsc=second dimension of the gsc array
  nband_k=number of bands at this k point for that spin polarization
  npw_k=number of plane waves at this k point
  nspinor=number of spinorial components of the wavefunctions (on current proc)
  subham(nband_k*(nband_k+1))=Hamiltonian expressed in the WFs subspace
  subovl(nband_k*(nband_k+1)*use_subovl)=overlap matrix expressed in the WFs subspace
  use_subovl=1 if the overlap matrix is not identity in WFs subspace
  usepaw= 0 for non paw calculation; =1 for paw calculation
  me_g0=1 if this processors has G=0, 0 otherwise.

OUTPUT

  eig_k(nband_k)=array for holding eigenvalues (hartree)
  evec(2*nband_k,nband_k)=array for holding eigenvectors

SIDE EFFECTS

  cg(2,mcg)=wavefunctions
  gsc(2,mgsc)=<g|S|c> matrix elements (S=overlap)

PARENTS

      rayleigh_ritz,vtowfk

CHILDREN

      abi_xcopy,abi_xgemm,abi_xhpev,abi_xhpgv,cg_normev,hermit

SOURCE

4825 subroutine subdiago(cg,eig_k,evec,gsc,icg,igsc,istwf_k,&
4826 &                   mcg,mgsc,nband_k,npw_k,nspinor,paral_kgb,&
4827 &                   subham,subovl,use_subovl,usepaw,me_g0)
4828 
4829  use m_linalg_interfaces
4830  use m_abi_linalg
4831 
4832 !This section has been created automatically by the script Abilint (TD).
4833 !Do not modify the following lines by hand.
4834 #undef ABI_FUNC
4835 #define ABI_FUNC 'subdiago'
4836 !End of the abilint section
4837 
4838  implicit none
4839 
4840 !Arguments ------------------------------------
4841  integer,intent(in) :: icg,igsc,istwf_k,mcg,mgsc,nband_k,npw_k,me_g0
4842  integer,intent(in) :: nspinor,paral_kgb,use_subovl,usepaw
4843  real(dp),intent(inout) :: subham(nband_k*(nband_k+1)),subovl(nband_k*(nband_k+1)*use_subovl)
4844  real(dp),intent(out) :: eig_k(nband_k),evec(2*nband_k,nband_k)
4845  real(dp),intent(inout) :: cg(2,mcg),gsc(2,mgsc)
4846 
4847 !Local variables-------------------------------
4848  integer :: iband,ii,ierr,rvectsize,vectsize,use_slk
4849  character(len=500) :: message
4850  ! real(dp) :: tsec(2)
4851  real(dp),allocatable :: evec_tmp(:,:),subovl_tmp(:),subham_tmp(:)
4852  real(dp),allocatable :: work(:,:)
4853  real(dp),allocatable :: blockvectora(:,:),blockvectorb(:,:),blockvectorc(:,:)
4854 
4855 ! *********************************************************************
4856 
4857  if (paral_kgb<0) then
4858    MSG_BUG('paral_kgb should be positive ')
4859  end if
4860 
4861  ! 1 if Scalapack version is used.
4862  use_slk = paral_kgb
4863 
4864  rvectsize=npw_k*nspinor
4865  vectsize=2*rvectsize;if (me_g0==1) vectsize=vectsize-1
4866 
4867 !Impose Hermiticity on diagonal elements of subham (and subovl, if needed)
4868 ! MG FIXME: In these two calls we are aliasing the args
4869  call hermit(subham,subham,ierr,nband_k)
4870  if (use_subovl==1) then
4871    call hermit(subovl,subovl,ierr,nband_k)
4872  end if
4873 
4874 !Diagonalize the Hamitonian matrix
4875  if(istwf_k==2) then
4876    ABI_ALLOCATE(evec_tmp,(nband_k,nband_k))
4877    ABI_ALLOCATE(subham_tmp,(nband_k*(nband_k+1)/2))
4878    subham_tmp=subham(1:nband_k*(nband_k+1):2)
4879    evec_tmp=zero
4880    if (use_subovl==1) then
4881      ABI_ALLOCATE(subovl_tmp,(nband_k*(nband_k+1)/2))
4882      subovl_tmp=subovl(1:nband_k*(nband_k+1):2)
4883 !    TO DO: Not sure this one has been fully tested
4884      call abi_xhpgv(1,'V','U',nband_k, subham_tmp,subovl_tmp, eig_k,evec_tmp,istwf_k=istwf_k,use_slk=use_slk)
4885      ABI_DEALLOCATE(subovl_tmp)
4886    else
4887      call abi_xhpev('V','U',nband_k,subham_tmp,eig_k,evec_tmp,istwf_k=istwf_k,use_slk=use_slk)
4888    end if
4889    evec(:,:)=zero;evec(1:2*nband_k:2,:) =evec_tmp
4890    ABI_DEALLOCATE(evec_tmp)
4891    ABI_DEALLOCATE(subham_tmp)
4892  else
4893    if (use_subovl==1) then
4894      call abi_xhpgv(1,'V','U',nband_k,subham,subovl,eig_k,evec,istwf_k=istwf_k,use_slk=use_slk)
4895    else
4896      call abi_xhpev('V','U',nband_k,subham,eig_k,evec,istwf_k=istwf_k,use_slk=use_slk)
4897    end if
4898  end if
4899 
4900 !Normalize each eigenvector and set phase:
4901 !The problem with minus/plus signs might be present also if .not. use_subovl
4902 !if(use_subovl == 0) then
4903  call cg_normev(evec,nband_k,nband_k)
4904 !end if
4905 
4906  if(istwf_k==2)then
4907    do iband=1,nband_k
4908      do ii=1,nband_k
4909        if(abs(evec(2*ii,iband))>1.0d-10)then
4910          write(message,'(3a,2i0,2es16.6,a,a)')ch10,&
4911 &         ' subdiago: For istwf_k=2, observed the following element of evec :',ch10,&
4912 &         iband,ii,evec(2*ii-1,iband),evec(2*ii,iband),ch10,'  with a non-negligible imaginary part.'
4913          MSG_BUG(message)
4914        end if
4915      end do
4916    end do
4917  end if
4918 
4919 !=====================================================
4920 !Carry out rotation of bands C(G,n) according to evecs
4921 ! ZGEMM if istwfk==1, DGEMM if istwfk==2
4922 !=====================================================
4923  if (istwf_k==2) then
4924 
4925    ABI_STAT_ALLOCATE(blockvectora,(vectsize,nband_k), ierr)
4926    ABI_CHECK(ierr==0, "out-of-memory in blockvectora")
4927    ABI_STAT_ALLOCATE(blockvectorb,(nband_k,nband_k), ierr)
4928    ABI_CHECK(ierr==0, "out-of-memory in blockvectorb")
4929    ABI_STAT_ALLOCATE(blockvectorc,(vectsize,nband_k), ierr)
4930    ABI_CHECK(ierr==0, "out-of-memory in blockvectorc")
4931 
4932    do iband=1,nband_k
4933      if (me_g0 == 1) then
4934        call abi_xcopy(1,cg(1,cgindex_subd(iband)),1,blockvectora(1,iband),1)
4935        call abi_xcopy(rvectsize-1,cg(1,cgindex_subd(iband)+1),2,blockvectora(2,iband),1)
4936        call abi_xcopy(rvectsize-1,cg(2,cgindex_subd(iband)+1),2,blockvectora(rvectsize+1,iband),1)
4937      else
4938        call abi_xcopy(rvectsize,cg(1,cgindex_subd(iband)),2,blockvectora(1,iband),1)
4939        call abi_xcopy(rvectsize,cg(2,cgindex_subd(iband)),2,blockvectora(rvectsize+1,iband),1)
4940      end if
4941      call abi_xcopy(nband_k,evec(2*iband-1,1),2*nband_k,blockvectorb(iband,1),nband_k)
4942    end do
4943 
4944    call abi_xgemm('N','N',vectsize,nband_k,nband_k,&
4945 &   cone,blockvectora,vectsize,blockvectorb,nband_k,czero,blockvectorc,vectsize)
4946 
4947    do iband=1,nband_k
4948      if (me_g0 == 1) then
4949        call abi_xcopy(1,blockvectorc(1,iband),1,cg(1,cgindex_subd(iband)),1)
4950        call abi_xcopy(rvectsize-1,blockvectorc(2,iband),1,cg(1,cgindex_subd(iband)+1),2)
4951        call abi_xcopy(rvectsize-1,blockvectorc(rvectsize+1,iband),1,cg(2,cgindex_subd(iband)+1),2)
4952      else
4953        call abi_xcopy(rvectsize,blockvectorc(1,iband),1,cg(1,cgindex_subd(iband)),2)
4954        call abi_xcopy(rvectsize,blockvectorc(rvectsize+1,iband),1,cg(2,cgindex_subd(iband)),2)
4955      end if
4956    end do
4957 
4958 !  If paw, musb also rotate S.C(G,n):
4959    if (usepaw==1) then
4960 
4961      do iband=1,nband_k
4962        if (me_g0 == 1) then
4963          call abi_xcopy(1,gsc(1,gscindex_subd(iband)),1,blockvectora(1,iband),1)
4964          call abi_xcopy(rvectsize-1,gsc(1,gscindex_subd(iband)+1),2,blockvectora(2,iband),1)
4965          call abi_xcopy(rvectsize-1,gsc(2,gscindex_subd(iband)+1),2,blockvectora(rvectsize+1,iband),1)
4966        else
4967          call abi_xcopy(rvectsize  ,gsc(1,gscindex_subd(iband)),2,blockvectora(1,iband),1)
4968          call abi_xcopy(rvectsize  ,gsc(2,gscindex_subd(iband)),2,blockvectora(rvectsize+1,iband),1)
4969        end if
4970        call abi_xcopy(nband_k,evec(2*iband-1,1),2*nband_k,blockvectorb(iband,1),nband_k)
4971      end do
4972 
4973      call abi_xgemm('N','N',vectsize,nband_k,nband_k,&
4974 &     cone,blockvectora,vectsize,blockvectorb,nband_k,czero,blockvectorc,vectsize)
4975 
4976      do iband=1,nband_k
4977        if (me_g0 == 1) then
4978          call abi_xcopy(1,blockvectorc(1,iband),1,gsc(1,gscindex_subd(iband)),1)
4979          call abi_xcopy(rvectsize-1,blockvectorc(2,iband),1,gsc(1,gscindex_subd(iband)+1),2)
4980          call abi_xcopy(rvectsize-1,blockvectorc(rvectsize+1,iband),1,gsc(2,gscindex_subd(iband)+1),2)
4981        else
4982          call abi_xcopy(rvectsize,blockvectorc(1,iband),1,gsc(1,gscindex_subd(iband)),2)
4983          call abi_xcopy(rvectsize,blockvectorc(rvectsize+1,iband),1,gsc(2,gscindex_subd(iband)),2)
4984        end if
4985      end do
4986 
4987    end if
4988 
4989    ABI_DEALLOCATE(blockvectora)
4990    ABI_DEALLOCATE(blockvectorb)
4991    ABI_DEALLOCATE(blockvectorc)
4992 
4993  else
4994 
4995    ABI_STAT_ALLOCATE(work,(2,npw_k*nspinor*nband_k), ierr)
4996    ABI_CHECK(ierr==0, "out-of-memory in work")
4997 
4998 !  MG: Do not remove this initialization.
4999 !  telast_06 stops in fxphase on inca_debug and little_buda (very very strange, due to atlas?)
5000    work=zero
5001 
5002    call abi_xgemm('N','N',npw_k*nspinor,nband_k,nband_k,cone, &
5003 &   cg(:,icg+1:npw_k*nspinor*nband_k+icg),npw_k*nspinor, &
5004 &   evec,nband_k,czero,work,npw_k*nspinor,x_cplx=2)
5005 
5006    call abi_xcopy(npw_k*nspinor*nband_k,work(1,1),1,cg(1,1+icg),1,x_cplx=2)
5007 
5008 !  If paw, must also rotate S.C(G,n):
5009    if (usepaw==1) then
5010      call abi_xgemm('N','N',npw_k*nspinor,nband_k,nband_k,cone, &
5011 &     gsc(:,1+igsc:npw_k*nspinor*nband_k+igsc),npw_k*nspinor, &
5012 &     evec,nband_k,czero,work,npw_k*nspinor,x_cplx=2)
5013      call abi_xcopy(npw_k*nspinor*nband_k, work(1,1),1,gsc(1,1+igsc),1,x_cplx=2)
5014    end if
5015 
5016    ABI_DEALLOCATE(work)
5017  end if
5018 
5019  contains
5020 
5021    function cgindex_subd(iband)
5022 
5023 
5024 !This section has been created automatically by the script Abilint (TD).
5025 !Do not modify the following lines by hand.
5026 #undef ABI_FUNC
5027 #define ABI_FUNC 'cgindex_subd'
5028 !End of the abilint section
5029 
5030    integer :: iband,cgindex_subd
5031    cgindex_subd=npw_k*nspinor*(iband-1)+icg+1
5032  end function cgindex_subd
5033    function gscindex_subd(iband)
5034 
5035 
5036 !This section has been created automatically by the script Abilint (TD).
5037 !Do not modify the following lines by hand.
5038 #undef ABI_FUNC
5039 #define ABI_FUNC 'gscindex_subd'
5040 !End of the abilint section
5041 
5042    integer :: iband,gscindex_subd
5043    gscindex_subd=npw_k*nspinor*(iband-1)+igsc+1
5044  end function gscindex_subd
5045 
5046 end subroutine subdiago

m_cgtools/cg_addtorho [ Functions ]

[ Top ] [ m_cgtools ] [ Functions ]

NAME

  cg_addtorho

FUNCTION

  Add |ur|**2 to the ground-states density rho.
    rho = rho + weight_r * Re[ur]**2 + weight_i * Im[ur]**2

INPUTS

  nx,ny,nz=physical dimension of the FFT box.
  ldx,ldy,ldz=leading dimensions of the arrays.
  ndat=number of contributions to accumulate.
  weight_r=weight used for the accumulation of the density in real space
  weight_i=weight used for the accumulation of the density in real space
  ur(2,ldx,ldy,ldz*ndat)=wavefunctions in real space

SIDE EFFECTS

  rho(ldx,ldy,ldz) = contains the input density at input,
                  modified in input with the contribution gived by ur.

PARENTS

      fourwf,m_dfti,m_fftw3

CHILDREN

SOURCE

2466 subroutine cg_addtorho(nx,ny,nz,ldx,ldy,ldz,ndat,weight_r,weight_i,ur,rho)
2467 
2468 
2469 !This section has been created automatically by the script Abilint (TD).
2470 !Do not modify the following lines by hand.
2471 #undef ABI_FUNC
2472 #define ABI_FUNC 'cg_addtorho'
2473 !End of the abilint section
2474 
2475  implicit none
2476 
2477 !Arguments ------------------------------------
2478 !scalars
2479  integer,intent(in) :: nx,ny,nz,ldx,ldy,ldz,ndat
2480  real(dp),intent(in) :: weight_i,weight_r
2481 !arrays
2482  real(dp),intent(in) :: ur(2,ldx,ldy,ldz*ndat)
2483  real(dp),intent(inout) :: rho(ldx,ldy,ldz)
2484 
2485 !Local variables-------------------------------
2486 !scalars
2487  integer :: ix,iy,iz,idat,izdat
2488 
2489 ! *************************************************************************
2490 
2491  if (ndat==1) then
2492 !$OMP PARALLEL DO
2493    do iz=1,nz
2494      do iy=1,ny
2495        do ix=1,nx
2496          rho(ix,iy,iz) = rho(ix,iy,iz) + weight_r * ur(1,ix,iy,iz)**2 &
2497 &                                      + weight_i * ur(2,ix,iy,iz)**2
2498        end do
2499      end do
2500    end do
2501 
2502  else
2503 ! It would be nice to use $OMP PARALLEL DO PRIVATE(izdat) REDUCTION(+:rho)
2504 ! but it's risky as the private rho is allocated on the stack of the thread.
2505 !$OMP PARALLEL PRIVATE(izdat)
2506    do idat=1,ndat
2507 !$OMP DO
2508      do iz=1,nz
2509        izdat = iz + (idat-1)*ldz
2510        do iy=1,ny
2511          do ix=1,nx
2512            rho(ix,iy,iz) = rho(ix,iy,iz) + weight_r * ur(1,ix,iy,izdat)**2 &
2513 &                                        + weight_i * ur(2,ix,iy,izdat)**2
2514          end do
2515        end do
2516      end do
2517 !$OMP END DO NOWAIT
2518    end do
2519 !$OMP END PARALLEL
2520  end if
2521 
2522 end subroutine cg_addtorho

m_cgtools/cg_box2gsph [ Functions ]

[ Top ] [ m_cgtools ] [ Functions ]

NAME

  cg_box2gsph

FUNCTION

INPUTS

  nx,ny,nz=physical dimension of the FFT box.
  ldx,ldy,ldz=Logical dimensions of the arrays.
  ndat=number of data in iarrbox
  npw_k=Number of planewaves in the G-sphere.
  kg_k(3,npw_k)=Reduced coordinates of the G-vectoes.
  iarrbox(2,ldx,ldy,ldz*ndat)=Input arrays on the FFT box.
  [rscal] = Scaling factor

OUTPUT

  oarrsph(2,npw_k*ndat)=Data defined on the G-sphere.

PARENTS

      fourwf,m_dfti,m_fftw3

CHILDREN

SOURCE

2352 subroutine cg_box2gsph(nx,ny,nz,ldx,ldy,ldz,ndat,npw_k,kg_k,iarrbox,oarrsph,rscal)
2353 
2354 
2355 !This section has been created automatically by the script Abilint (TD).
2356 !Do not modify the following lines by hand.
2357 #undef ABI_FUNC
2358 #define ABI_FUNC 'cg_box2gsph'
2359 !End of the abilint section
2360 
2361  implicit none
2362 
2363 !Arguments ------------------------------------
2364 !scalars
2365  integer,intent(in) :: npw_k,nx,ny,nz,ldx,ldy,ldz,ndat
2366  real(dp),optional,intent(in) :: rscal
2367 !arrays
2368  integer,intent(in) :: kg_k(3,npw_k)
2369  real(dp),intent(in) :: iarrbox(2,ldx*ldy*ldz*ndat)
2370  real(dp),intent(out) :: oarrsph(2,npw_k*ndat)
2371 
2372 !Local variables-------------------------------
2373 !scalars
2374  integer :: ig,ix,iy,iz,idat,sph_pad,box_pad,ifft
2375 
2376 ! *************************************************************************
2377 
2378  if (.not. PRESENT(rscal)) then
2379    !
2380    if (ndat==1) then
2381 !$OMP PARALLEL DO PRIVATE(ix,iy,iz,ifft)
2382      do ig=1,npw_k
2383        ix=kg_k(1,ig); if (ix<0) ix=ix+nx; ix=ix+1
2384        iy=kg_k(2,ig); if (iy<0) iy=iy+ny; iy=iy+1
2385        iz=kg_k(3,ig); if (iz<0) iz=iz+nz; iz=iz+1
2386        ifft = ix + (iy-1)*ldx + (iz-1)*ldx*ldy
2387        oarrsph(1,ig) = iarrbox(1,ifft)
2388        oarrsph(2,ig) = iarrbox(2,ifft)
2389      end do
2390    else
2391 !$OMP PARALLEL DO PRIVATE(sph_pad,box_pad,ix,iy,iz,ifft)
2392      do idat=1,ndat
2393        sph_pad = (idat-1)*npw_k
2394        box_pad = (idat-1)*ldx*ldy*ldz
2395        do ig=1,npw_k
2396          ix=kg_k(1,ig); if (ix<0) ix=ix+nx; ix=ix+1
2397          iy=kg_k(2,ig); if (iy<0) iy=iy+ny; iy=iy+1
2398          iz=kg_k(3,ig); if (iz<0) iz=iz+nz; iz=iz+1
2399          ifft = ix + (iy-1)*ldx + (iz-1)*ldx*ldy
2400          oarrsph(1,ig+sph_pad) = iarrbox(1,ifft+box_pad)
2401          oarrsph(2,ig+sph_pad) = iarrbox(2,ifft+box_pad)
2402        end do
2403      end do
2404    end if
2405    !
2406  else
2407    if (ndat==1) then
2408 !$OMP PARALLEL DO PRIVATE(ix,iy,iz,ifft)
2409      do ig=1,npw_k
2410        ix=kg_k(1,ig); if (ix<0) ix=ix+nx; ix=ix+1
2411        iy=kg_k(2,ig); if (iy<0) iy=iy+ny; iy=iy+1
2412        iz=kg_k(3,ig); if (iz<0) iz=iz+nz; iz=iz+1
2413        ifft = ix + (iy-1)*ldx + (iz-1)*ldx*ldy
2414        oarrsph(1,ig) = iarrbox(1,ifft) * rscal
2415        oarrsph(2,ig) = iarrbox(2,ifft) * rscal
2416      end do
2417    else
2418 !$OMP PARALLEL DO PRIVATE(sph_pad,box_pad,ix,iy,iz,ifft)
2419      do idat=1,ndat
2420        sph_pad = (idat-1)*npw_k
2421        box_pad = (idat-1)*ldx*ldy*ldz
2422        do ig=1,npw_k
2423          ix=kg_k(1,ig); if (ix<0) ix=ix+nx; ix=ix+1
2424          iy=kg_k(2,ig); if (iy<0) iy=iy+ny; iy=iy+1
2425          iz=kg_k(3,ig); if (iz<0) iz=iz+nz; iz=iz+1
2426          ifft = ix + (iy-1)*ldx + (iz-1)*ldx*ldy
2427          oarrsph(1,ig+sph_pad) = iarrbox(1,ifft+box_pad) * rscal
2428          oarrsph(2,ig+sph_pad) = iarrbox(2,ifft+box_pad) * rscal
2429        end do
2430      end do
2431    end if
2432  end if
2433 
2434 end subroutine cg_box2gsph

m_cgtools/cg_dznrm2 [ Functions ]

[ Top ] [ m_cgtools ] [ Functions ]

NAME

  cg_dznrm2

FUNCTION

   returns the euclidean norm of a vector via the function name, so that
   DZNRM2 := sqrt( x**H*x )

INPUTS

  n = Specifies the number of elements in vector x.
  x(2*x) = Input array.

OUTPUT

PARENTS

SOURCE

631 function cg_dznrm2(n, x) result(res)
632 
633 
634 !This section has been created automatically by the script Abilint (TD).
635 !Do not modify the following lines by hand.
636 #undef ABI_FUNC
637 #define ABI_FUNC 'cg_dznrm2'
638 !End of the abilint section
639 
640  implicit none
641 
642 !Arguments ------------------------------------
643 !scalars
644  integer,intent(in) :: n
645  real(dp) :: res
646 !arrays
647  real(dp),intent(in) :: x(2*n)
648  real(dp),external :: dznrm2
649 
650 ! *************************************************************************
651 
652  res = dznrm2(n, x, 1)
653 
654 end function cg_dznrm2

m_cgtools/cg_envlop [ Functions ]

[ Top ] [ m_cgtools ] [ Functions ]

NAME

 cg_envlop

FUNCTION

 Multiply random number values in cg by envelope function to lower initial kinetic energy.
 Envelope  $\left( 1-\left( G/G_{\max }\right) ^2\right) ^{power}$ for |G|<= Gmax.
 Near G=0, little scaling, and goes to zero flatly near Gmax.Loop over perturbations

INPUTS

 cg(2,npw*nband)=initial random number wavefunctions
 ecut=kinetic energy cutoff in Ha
 gmet(3,3)=reciprocal space metric (bohr^-2)
 icgmod=shift to be given to the location of data in cg
 kg(3,npw)=reduced coordinates of G vectors in basis sphere
 kpoint(3)=reduced coordinates of k point
 mcg=maximum second dimension of cg (at least npw*nband*nspinor)
 nband=number of bands being considered
 npw=number of planewaves in basis sphere
 nspinor=number of spinorial components of the wavefunctions

OUTPUT

  cg(2,mcg)=revised values (not orthonormalized)

PARENTS

      wfconv

CHILDREN

SOURCE

3675 subroutine cg_envlop(cg,ecut,gmet,icgmod,kg,kpoint,mcg,nband,npw,nspinor)
3676 
3677 
3678 !This section has been created automatically by the script Abilint (TD).
3679 !Do not modify the following lines by hand.
3680 #undef ABI_FUNC
3681 #define ABI_FUNC 'cg_envlop'
3682 !End of the abilint section
3683 
3684  implicit none
3685 
3686 !Arguments ------------------------------------
3687 !scalars
3688  integer,intent(in) :: icgmod,mcg,nband,npw,nspinor
3689  real(dp),intent(in) :: ecut
3690 !arrays
3691  integer,intent(in) :: kg(3,npw)
3692  real(dp),intent(in) :: gmet(3,3),kpoint(3)
3693  real(dp),intent(inout) :: cg(2,mcg)
3694 
3695 !Local variables-------------------------------
3696 !scalars
3697  integer,parameter :: re=1,im=2,power=12
3698  integer :: i1,i2,i3,ig,igs,ispinor,nn,spad
3699  real(dp) :: cutoff,gs,kpgsqc
3700  !character(len=500) :: msg
3701 !arrays
3702  real(dp),allocatable :: cut_pws(:)
3703 
3704 ! *************************************************************************
3705 
3706 !$(k+G)^2$ cutoff from $(1/2)(2 Pi (k+G))^2 = ecut$
3707  kpgsqc=ecut/(2.0_dp*pi**2)
3708  cutoff = kpgsqc
3709 
3710  ABI_MALLOC(cut_pws,(npw))
3711 
3712 !Run through G vectors in basis
3713 !$OMP PARALLEL DO PRIVATE(gs)
3714  do ig=1,npw
3715    i1=kg(1,ig) ; i2=kg(2,ig) ; i3=kg(3,ig)
3716 !(k+G)^2 evaluated using metric and kpoint
3717    gs = gmet(1,1)*(kpoint(1)+dble(i1))**2+&
3718 &    gmet(2,2)*(kpoint(2)+dble(i2))**2+&
3719 &    gmet(3,3)*(kpoint(3)+dble(i3))**2+&
3720 &    2.0_dp*(gmet(2,1)*(kpoint(2)+dble(i2))*(kpoint(1)+dble(i1))+&
3721 &    gmet(3,2)*(kpoint(3)+dble(i3))*(kpoint(2)+dble(i2))+&
3722 &    gmet(1,3)*(kpoint(1)+dble(i1))*(kpoint(3)+dble(i3)))
3723    if (gs>cutoff) then
3724      cut_pws(ig) = zero
3725    else
3726      cut_pws(ig) = (1.0_dp-(gs/cutoff))**power
3727    end if
3728  end do
3729 
3730 !Run through bands (real and imaginary components)
3731 !$OMP PARALLEL DO PRIVATE(spad,igs)
3732  do nn=1,nband
3733    spad = (nn-1)*npw*nspinor+icgmod
3734    do ispinor=1,nspinor
3735      do ig=1,npw
3736        igs=ig+(ispinor-1)*npw
3737        cg(1,igs+spad) = cg(1,igs+spad) * cut_pws(ig)
3738        cg(2,igs+spad) = cg(2,igs+spad) * cut_pws(ig)
3739      end do
3740    end do
3741  end do
3742 
3743  ABI_FREE(cut_pws)
3744 
3745 end subroutine cg_envlop

m_cgtools/cg_filter [ Functions ]

[ Top ] [ m_cgtools ] [ Functions ]

NAME

  cg_filter

FUNCTION

  Set all the elements of x to zero where mask is .TRUE.

INPUTS

  n=Specifies the number of elements in vectors x and y.
  mask(n)=Logical array.

SIDE EFFECTS

  x(n)=See description.

PARENTS

SOURCE

319 pure subroutine cg_filter(n, x, mask)
320 
321 
322 !This section has been created automatically by the script Abilint (TD).
323 !Do not modify the following lines by hand.
324 #undef ABI_FUNC
325 #define ABI_FUNC 'cg_filter'
326 !End of the abilint section
327 
328  implicit none
329 
330 !Arguments ------------------------------------
331 !scalars
332  integer,intent(in) :: n
333 !arrays
334  real(dp),intent(inout) :: x(2,n)
335  logical,intent(in) :: mask(n)
336 
337 ! *************************************************************************
338 
339  where (mask)
340    x(1,:) = zero
341    x(2,:) = zero
342  end where
343 
344 end subroutine cg_filter

m_cgtools/cg_from_reim [ Functions ]

[ Top ] [ m_cgtools ] [ Functions ]

NAME

  cg_from_reim

FUNCTION

INPUTS

PARENTS

CHILDREN

SOURCE

480 subroutine cg_from_reim(npw,ndat,reim,factor,cg)
481 
482 
483 !This section has been created automatically by the script Abilint (TD).
484 !Do not modify the following lines by hand.
485 #undef ABI_FUNC
486 #define ABI_FUNC 'cg_from_reim'
487 !End of the abilint section
488 
489  implicit none
490 
491 !Arguments ------------------------------------
492 !scalars
493  integer,intent(in) :: npw,ndat
494  real(dp),intent(in) :: factor
495 !arrays
496  real(dp),intent(in) :: reim(npw*ndat*2)
497  real(dp),intent(out) :: cg(2*npw*ndat)
498 
499 ! *************************************************************************
500 
501  ! UnPack real and imaginary part and multiply by scale factor if /= one.
502  ! Could use blocking but oh well
503  call dcopy(npw*ndat, reim(1), 1, cg(1), 2)
504  call dcopy(npw*ndat, reim(npw*ndat+1), 1, cg(2), 2)
505 
506  if (factor /= one) call dscal(2*npw*ndat,factor, cg(1), 1)
507 
508 end subroutine cg_from_reim

m_cgtools/cg_fromcplx [ Functions ]

[ Top ] [ m_cgtools ] [ Functions ]

NAME

  cg_fromcplx

FUNCTION

  Convert a complex array to a real array with (real,imag) part

INPUTS

  n = Specifies the number of elements in icplx and ocg.
  icplx(n)=Input complex array.

OUTPUT

  ocg(2*n)=Output array with real and imaginary part.

PARENTS

CHILDREN

SOURCE

265 subroutine cg_fromcplx(n,icplx,ocg)
266 
267 
268 !This section has been created automatically by the script Abilint (TD).
269 !Do not modify the following lines by hand.
270 #undef ABI_FUNC
271 #define ABI_FUNC 'cg_fromcplx'
272 !End of the abilint section
273 
274  implicit none
275 
276 !Arguments ------------------------------------
277 !scalars
278  integer,intent(in) :: n
279 !arrays
280  real(dp),intent(out) :: ocg(2*n)
281  complex(dpc),intent(in) :: icplx(n)
282 
283 !Local variables ------------------------------
284 !scalars
285  integer :: ii,idx
286 
287 ! *************************************************************************
288 
289 !$OMP PARALLEL DO PRIVATE(ii,idx)
290  do ii=1,n
291    idx = 2*ii-1
292    ocg(idx  ) = DBLE (icplx(ii))
293    ocg(idx+1) = AIMAG(icplx(ii))
294  end do
295 
296 end subroutine cg_fromcplx

m_cgtools/cg_getspin [ Functions ]

[ Top ] [ m_cgtools ] [ Functions ]

NAME

 cg_getspin

FUNCTION

  Sandwich a single wave function on the Pauli matrices

INPUTS

  npw_k = number of plane waves
  cgcband = coefficients of spinorial wave function

OUTPUT

  spin = 3-vector of spin components for this state
  cgcmat = outer spin product of spinorial wf with itself

PARENTS

      m_cut3d,partial_dos_fractions

CHILDREN

SOURCE

2103 subroutine  cg_getspin(cgcband, npw_k, spin, cgcmat)
2104 
2105 
2106 !This section has been created automatically by the script Abilint (TD).
2107 !Do not modify the following lines by hand.
2108 #undef ABI_FUNC
2109 #define ABI_FUNC 'cg_getspin'
2110 !End of the abilint section
2111 
2112  implicit none
2113 
2114 !Arguments ------------------------------------
2115 !scalars
2116  integer, intent(in) :: npw_k
2117  real(dp), intent(in) :: cgcband(2,2*npw_k)
2118  complex(dpc), intent(out),optional :: cgcmat(2,2)
2119  real(dp), intent(out) :: spin(3)
2120 
2121 !Local variables-------------------------------
2122 !scalars
2123  complex(dpc) :: cspin(0:3), cgcmat_(2,2)
2124 ! ***********************************************************************
2125 
2126 ! cgcmat_ = cgcband * cgcband^T*  i.e. 2x2 matrix of spin components (dpcomplex)
2127  cgcmat_ = czero
2128  call zgemm('n','c',2,2,npw_k,cone,cgcband,2,cgcband,2,czero,cgcmat_,2)
2129 
2130 ! spin(*)  = sum_{si sj pi} cgcband(si,pi)^* pauli_mat*(si,sj) cgcband(sj,pi)
2131  cspin(0) = cgcmat_(1,1)*pauli_mat(1,1,0) + cgcmat_(2,1)*pauli_mat(2,1,0) &
2132 &         + cgcmat_(1,2)*pauli_mat(1,2,0) + cgcmat_(2,2)*pauli_mat(2,2,0)
2133  cspin(1) = cgcmat_(1,1)*pauli_mat(1,1,1) + cgcmat_(2,1)*pauli_mat(2,1,1) &
2134 &         + cgcmat_(1,2)*pauli_mat(1,2,1) + cgcmat_(2,2)*pauli_mat(2,2,1)
2135  cspin(2) = cgcmat_(1,1)*pauli_mat(1,1,2) + cgcmat_(2,1)*pauli_mat(2,1,2) &
2136 &         + cgcmat_(1,2)*pauli_mat(1,2,2) + cgcmat_(2,2)*pauli_mat(2,2,2)
2137  cspin(3) = cgcmat_(1,1)*pauli_mat(1,1,3) + cgcmat_(2,1)*pauli_mat(2,1,3) &
2138 &         + cgcmat_(1,2)*pauli_mat(1,2,3) + cgcmat_(2,2)*pauli_mat(2,2,3)
2139 !write(std_out,*) 'cgmat: ', cgcmat_
2140 !write(std_out,*) 'real(spin): ', real(cspin)
2141 !write(std_out,*) 'aimag(spin): ', aimag(cspin)
2142 
2143  spin = real(cspin(1:3))
2144  if (present(cgcmat)) cgcmat = cgcmat_
2145 
2146 end subroutine cg_getspin

m_cgtools/cg_gsph2box [ Functions ]

[ Top ] [ m_cgtools ] [ Functions ]

NAME

 cg_gsph2box

FUNCTION

 Array iarrsph is defined in sphere with npw_k points. Insert iarrsph inside box
 of nx*ny*nz points to define array oarrbox for fft box. rest of oarrbox is filled with 0 s.

INPUTS

 iarrsph(2,npw_k*ndat)= contains values for npw_k G vectors in basis sphere
 ndat=number of FFT to perform.
 npw_k=number of G vectors in basis at this k point
 oarrbox(2,ldx*ldy*ldz*ndat) = fft box
 nx,ny,nz=physical dimension of the box (oarrbox)
 ldx,ldy,ldz=memory dimension of oarrbox
 kg_k(3,npw_k)=integer coordinates of G vectors in basis sphere
 istwf_k=option parameter that describes the storage of wfs

OUTPUT

   oarrbox(ldx*ldy*ldz*ndat)

NOTES

 If istwf_k differs from 1, then special storage modes must be taken
 into account, for symmetric wavefunctions coming from k=(0 0 0) or other
 special k points.

PARENTS

CHILDREN

SOURCE

2183 subroutine cg_gsph2box(nx,ny,nz,ldx,ldy,ldz,ndat,npw_k,istwf_k,kg_k,iarrsph,oarrbox)
2184 
2185 
2186 !This section has been created automatically by the script Abilint (TD).
2187 !Do not modify the following lines by hand.
2188 #undef ABI_FUNC
2189 #define ABI_FUNC 'cg_gsph2box'
2190 !End of the abilint section
2191 
2192  implicit none
2193 
2194 !Arguments ------------------------------------
2195 !scalars
2196  integer,intent(in) :: istwf_k,nx,ny,nz,ldx,ldy,ldz,ndat,npw_k
2197 !arrays
2198  integer,intent(in) :: kg_k(3,npw_k)
2199  real(dp),intent(in) :: iarrsph(2,npw_k*ndat)
2200  real(dp),intent(out) :: oarrbox(2,ldx*ldy*ldz*ndat)
2201 
2202 !Local variables-------------------------------
2203 !scalars
2204  integer,parameter :: me_g0=1
2205  integer :: ix,ixinv,iy,iyinv,iz,izinv,dat,ipw,npwmin,pad_box,pad_sph,ifft,ifft_inv,ldxyz
2206  character(len=500) :: msg
2207 !arrays
2208  integer,allocatable :: ixinver(:),iyinver(:),izinver(:)
2209 
2210 ! *************************************************************************
2211 
2212 !In the case of special k-points, invariant under time-reversal,
2213 !but not Gamma, initialize the inverse coordinates
2214 !Remember indeed that
2215 !u_k(G) = u_{k+G0}(G-G0); u_{-k}(G) = u_k(G)^*
2216 !and therefore:
2217 !u_{G0/2}(G) = u_{G0/2}(-G-G0)^*.
2218  if (istwf_k>=2) then
2219    ABI_MALLOC(ixinver,(nx))
2220    ABI_MALLOC(iyinver,(ny))
2221    ABI_MALLOC(izinver,(nz))
2222    if ( ANY(istwf_k==(/2,4,6,8/)) ) then
2223      ixinver(1)=1
2224      do ix=2,nx
2225        ixinver(ix)=nx+2-ix
2226      end do
2227    else
2228      do ix=1,nx
2229        ixinver(ix)=nx+1-ix
2230      end do
2231    end if
2232    if (istwf_k>=2 .and. istwf_k<=5) then
2233      iyinver(1)=1
2234      do iy=2,ny
2235        iyinver(iy)=ny+2-iy
2236      end do
2237    else
2238      do iy=1,ny
2239        iyinver(iy)=ny+1-iy
2240      end do
2241    end if
2242    if ( ANY(istwf_k==(/2,3,6,7/)) ) then
2243      izinver(1)=1
2244      do iz=2,nz
2245        izinver(iz)=nz+2-iz
2246      end do
2247    else
2248      do iz=1,nz
2249        izinver(iz)=nz+1-iz
2250      end do
2251    end if
2252  end if
2253 
2254  ldxyz = ldx*ldy*ldz
2255 
2256  if (istwf_k==1) then
2257 
2258 !$OMP PARALLEL DO PRIVATE(pad_sph,pad_box,ix,iy,iz,ifft)
2259    do dat=1,ndat
2260      pad_sph = (dat-1)*npw_k
2261      pad_box = (dat-1)*ldxyz
2262      oarrbox(:,1+pad_box:ldxyz+pad_box) = zero ! zero the sub-array
2263      do ipw=1,npw_k
2264        ix=kg_k(1,ipw); if (ix<0) ix=ix+nx; ix=ix+1
2265        iy=kg_k(2,ipw); if (iy<0) iy=iy+ny; iy=iy+1
2266        iz=kg_k(3,ipw); if (iz<0) iz=iz+nz; iz=iz+1
2267        ifft = ix + (iy-1)*ldx + (iz-1)*ldx*ldy
2268 #if defined __INTEL_COMPILER && defined HAVE_OPENMP
2269        if (ifft==0) then
2270          MSG_ERROR("prevent ifort+OMP from miscompiling this section on cronos")
2271        end if
2272 #endif
2273        oarrbox(1,ifft+pad_box) = iarrsph(1,ipw+pad_sph)
2274        oarrbox(2,ifft+pad_box) = iarrsph(2,ipw+pad_sph)
2275      end do
2276    end do
2277 
2278  else if (istwf_k>=2) then
2279    !
2280    npwmin=1
2281    if(istwf_k==2 .and. me_g0==1) then ! If gamma point, then oarrbox must be completed
2282      do dat=1,ndat
2283        pad_sph = (dat-1)*npw_k
2284        pad_box = (dat-1)*ldxyz
2285        oarrbox(1,1+pad_box) = iarrsph(1,1+pad_sph)
2286        oarrbox(2,1+pad_box) = zero
2287      end do
2288      npwmin=2
2289    end if
2290 
2291 !$OMP PARALLEL DO PRIVATE(pad_sph,pad_box,ix,iy,iz,ixinv,iyinv,izinv,ifft,ifft_inv)
2292    do dat=1,ndat
2293      pad_sph = (dat-1)*npw_k
2294      pad_box = (dat-1)*ldxyz
2295      oarrbox(:,npwmin+pad_box:ldxyz+pad_box) = zero
2296      do ipw=npwmin,npw_k
2297        ix=kg_k(1,ipw); if(ix<0)ix=ix+nx; ix=ix+1
2298        iy=kg_k(2,ipw); if(iy<0)iy=iy+ny; iy=iy+1
2299        iz=kg_k(3,ipw); if(iz<0)iz=iz+nz; iz=iz+1
2300        ifft = ix + (iy-1)*ldx + (iz-1)*ldx*ldy
2301        ! Construct the coordinates of -k-G
2302        ixinv=ixinver(ix); iyinv=iyinver(iy); izinv=izinver(iz)
2303        ifft_inv = ixinv + (iyinv-1)*ldx + (izinv-1)*ldx*ldy
2304 
2305        oarrbox(:,ifft    +pad_box) =  iarrsph(:,ipw+pad_sph)
2306        oarrbox(1,ifft_inv+pad_box) =  iarrsph(1,ipw+pad_sph)
2307        oarrbox(2,ifft_inv+pad_box) = -iarrsph(2,ipw+pad_sph)
2308      end do
2309    end do
2310    !
2311  else
2312    write(msg,'(a,i0)')"Wrong istwfk ",istwf_k
2313    MSG_ERROR(msg)
2314  end if
2315 
2316  if (istwf_k>=2) then
2317    ABI_FREE(ixinver)
2318    ABI_FREE(iyinver)
2319    ABI_FREE(izinver)
2320  end if
2321 
2322 end subroutine cg_gsph2box

m_cgtools/cg_normev [ Functions ]

[ Top ] [ m_cgtools ] [ Functions ]

NAME

 cg_normev

FUNCTION

 Normalize a set of nband eigenvectors of complex length npw
 (real length 2*npw) and set phases to make cg(i,i) real and positive.
 Near convergence, cg(i,j) approaches delta(i,j).

INPUTS

  cg(2*npw,nband)=unnormalized eigenvectors
  npw=dimension of cg as shown
  nband=number of eigenvectors and complex length thereof.

OUTPUT

  cg(2*npw,nband)=nband normalized eigenvectors

PARENTS

      subdiago

CHILDREN

SOURCE

3774 subroutine cg_normev(cg,npw,nband)
3775 
3776 
3777 !This section has been created automatically by the script Abilint (TD).
3778 !Do not modify the following lines by hand.
3779 #undef ABI_FUNC
3780 #define ABI_FUNC 'cg_normev'
3781 !End of the abilint section
3782 
3783  implicit none
3784 
3785 !Arguments ------------------------------------
3786 !scalars
3787  integer,intent(in) :: npw,nband
3788 !arrays
3789  real(dp),intent(inout) :: cg(2*npw,nband)
3790 
3791 !Local variables-------------------------------
3792 !scalars
3793  integer :: ii,jj
3794  real(dp) :: den,evim,evre,phim,phre,xnorm
3795  character(len=500) :: message
3796 
3797 ! *************************************************************************
3798 
3799 !Loop over vectors
3800  do ii=1,nband
3801    ! find norm
3802    xnorm=0.0d0
3803    do jj=1,2*npw
3804      xnorm=xnorm+cg(jj,ii)**2
3805    end do
3806 
3807    if((xnorm-one)**2>tol6)then
3808      write(message,'(6a,i6,a,es16.6,3a)' )ch10,&
3809 &     'normev: ',ch10,&
3810 &     'Starting xnorm should be close to one (tol is tol6).',ch10,&
3811 &     'However, for state number',ii,', xnorm=',xnorm,ch10,&
3812 &     'It might be that your LAPACK library has not been correctly installed.'
3813      MSG_BUG(message)
3814    end if
3815 
3816    xnorm=1.0d0/sqrt(xnorm)
3817 !  Set up phase
3818    phre=cg(2*ii-1,ii)
3819    phim=cg(2*ii,ii)
3820    if (phim/=0.0d0) then
3821      den=1.0d0/sqrt(phre**2+phim**2)
3822      phre=phre*xnorm*den
3823      phim=phim*xnorm*den
3824    else
3825 !    give xnorm the same sign as phre (negate if negative)
3826      phre=sign(xnorm,phre)
3827      phim=0.0d0
3828    end if
3829 !  normalize with phase change
3830    do jj=1,2*npw,2
3831      evre=cg(jj,ii)
3832      evim=cg(jj+1,ii)
3833      cg(jj,ii)=phre*evre+phim*evim
3834      cg(jj+1,ii)=phre*evim-phim*evre
3835    end do
3836  end do
3837 
3838 end subroutine cg_normev

m_cgtools/cg_precon [ Functions ]

[ Top ] [ m_cgtools ] [ Functions ]

NAME

 cg_precon

FUNCTION

 precondition $<G|(H-e_{n,k})|C_{n,k}>$

INPUTS

  $cg(2,npw)=<G|C_{n,k}>$.
  $eval=current band eigenvalue=<C_{n,k}|H|C_{n,k}>$.
  istwf_k=option parameter that describes the storage of wfs
  kinpw(npw)=(modified) kinetic energy for each plane wave (Hartree)
  nspinor=number of spinorial components of the wavefunctions
  $vect(2,npw)=<G|H|C_{n,k}>$.
  npw=number of planewaves at this k point.
  optekin= 1 if the kinetic energy used in preconditionning is modified
             according to Kresse, Furthmuller, PRB 54, 11169 (1996) [[cite:Kresse1996]]
           0 otherwise
  mg_g0=1 if the node treats G0.
  comm=MPI communicator

OUTPUT

  pcon(npw)=preconditionning matrix
  vect(2,npw*nspinor)=<G|(H-eval)|C_{n,k}>*(polynomial ratio)

PARENTS

      cgwf,dfpt_cgwf

CHILDREN

SOURCE

3876 subroutine cg_precon(cg,eval,istwf_k,kinpw,npw,nspinor,me_g0,optekin,pcon,vect,comm)
3877 
3878 
3879 !This section has been created automatically by the script Abilint (TD).
3880 !Do not modify the following lines by hand.
3881 #undef ABI_FUNC
3882 #define ABI_FUNC 'cg_precon'
3883 !End of the abilint section
3884 
3885  implicit none
3886 
3887 !Arguments ------------------------------------
3888 !scalars
3889  integer,intent(in) :: istwf_k,npw,nspinor,optekin,me_g0,comm
3890  real(dp),intent(in) :: eval
3891 !arrays
3892  real(dp),intent(in) :: cg(2,npw*nspinor),kinpw(npw)
3893  real(dp),intent(inout) :: pcon(npw),vect(2,npw*nspinor)
3894 
3895 !Local variables-------------------------------
3896 !scalars
3897  integer :: ierr,ig,igs,ipw1,ispinor
3898  real(dp) :: ek0,ek0_inv,fac,poly,xx
3899  character(len=500) :: message
3900 !arrays
3901  real(dp) :: tsec(2)
3902 
3903 ! *************************************************************************
3904 
3905 !Compute mean kinetic energy of band
3906  if(istwf_k==1)then
3907    ek0=zero
3908    do ispinor=1,nspinor
3909      igs=(ispinor-1)*npw
3910 !    !$OMP PARALLEL DO PRIVATE(ig) REDUCTION(+:ek0) SHARED(cg,igs,kinpw,npw)
3911      do ig=1+igs,npw+igs
3912        if(kinpw(ig-igs)<huge(0.0_dp)*1.d-11)then
3913          ek0=ek0+kinpw(ig-igs)*(cg(1,ig)**2+cg(2,ig)**2)
3914        end if
3915      end do
3916    end do
3917 
3918  else if (istwf_k>=2)then
3919    if (istwf_k==2 .and. me_g0 == 1)then
3920      ek0=zero ; ipw1=2
3921      if(kinpw(1)<huge(0.0_dp)*1.d-11)ek0=0.5_dp*kinpw(1)*cg(1,1)**2
3922    else
3923      ek0=zero ; ipw1=1
3924    end if
3925    do ispinor=1,nspinor
3926      igs=(ispinor-1)*npw
3927 !    !$OMP PARALLEL DO PRIVATE(ig) REDUCTION(+:ek0) SHARED(cg,ipw1,kinpw,npw)
3928      do ig=ipw1+igs,npw+igs
3929        if(kinpw(ig)<huge(0.0_dp)*1.d-11)then
3930          ek0=ek0+kinpw(ig)*(cg(1,ig)**2+cg(2,ig)**2)
3931        end if
3932      end do
3933    end do
3934    ek0=2.0_dp*ek0
3935  end if
3936 
3937  call timab(48,1,tsec)
3938  call xmpi_sum(ek0,comm,ierr)
3939  call timab(48,2,tsec)
3940 
3941  if(ek0<1.0d-10)then
3942    write(message,'(3a)')&
3943 &   'The mean kinetic energy of a wavefunction vanishes.',ch10,&
3944 &   'It is reset to 0.1Ha.'
3945    MSG_WARNING(message)
3946    ek0=0.1_dp
3947  end if
3948 
3949  if (optekin==1) then
3950    ek0_inv=2.0_dp/(3._dp*ek0)
3951  else
3952    ek0_inv=1.0_dp/ek0
3953  end if
3954 
3955 !Carry out preconditioning
3956  do ispinor=1,nspinor
3957    igs=(ispinor-1)*npw
3958 !$OMP PARALLEL DO PRIVATE(fac,ig,poly,xx) SHARED(cg,ek0_inv,eval,kinpw,igs,npw,vect,pcon)
3959    do ig=1+igs,npw+igs
3960      if(kinpw(ig-igs)<huge(0.0_dp)*1.d-11)then
3961        xx=kinpw(ig-igs)*ek0_inv
3962 !      Teter polynomial ratio
3963        poly=27._dp+xx*(18._dp+xx*(12._dp+xx*8._dp))
3964        fac=poly/(poly+16._dp*xx**4)
3965        if (optekin==1) fac=two*fac
3966        pcon(ig-igs)=fac
3967        vect(1,ig)=( vect(1,ig)-eval*cg(1,ig) )*fac
3968        vect(2,ig)=( vect(2,ig)-eval*cg(2,ig) )*fac
3969      else
3970        pcon(ig-igs)=zero
3971        vect(1,ig)=zero
3972        vect(2,ig)=zero
3973      end if
3974    end do
3975  end do
3976 
3977 end subroutine cg_precon

m_cgtools/cg_precon_block [ Functions ]

[ Top ] [ m_cgtools ] [ Functions ]

NAME

 cg_precon_block

FUNCTION

 precondition $<G|(H-e_{n,k})|C_{n,k}>$ for a block of band (band-FFT parallelisation)
 in the case of real WFs (istwfk/=1)

INPUTS

  blocksize= size of blocks of bands
  $cg(vectsize,blocksize)=<G|C_{n,k}> for a block of bands$.
  $eval(blocksize,blocksize)=current block of bands eigenvalues=<C_{n,k}|H|C_{n,k}>$.
  $ghc(vectsize,blocksize)=<G|H|C_{n,k}> for a block of bands$.
  iterationnumber=number of iterative minimizations in LOBPCG
  kinpw(npw)=(modified) kinetic energy for each plane wave (Hartree)
  nspinor=number of spinorial components of the wavefunctions (on current proc)
  $vect(vectsize,blocksize)=<G|H|C_{n,k}> for a block of bands$.
  npw=number of planewaves at this k point.
  optekin= 1 if the kinetic energy used in preconditionning is modified
             according to Kresse, Furthmuller, PRB 54, 11169 (1996) [[cite:Kresse1996]]
           0 otherwise
  optpcon= 0 the TPA preconditionning matrix does not depend on band
           1 the TPA preconditionning matrix (not modified)
           2 the TPA preconditionning matrix is independant of iterationnumber
  vectsize= size of vectors
  mg_g0=1 if this node has Gamma, 0 otherwise.

OUTPUT

  vect(2,npw)=<g|(h-eval)|c_{n,k}>*(polynomial ratio)

SIDE EFFECTS

  pcon(npw,blocksize)=preconditionning matrix
            input  if optpcon=0,2 and iterationnumber/=1
            output if optpcon=0,2 and iterationnumber==1

PARENTS

      m_lobpcg

CHILDREN

SOURCE

4023 subroutine cg_precon_block(cg,eval,blocksize,iterationnumber,kinpw,&
4024 & npw,nspinor,me_g0,optekin,optpcon,pcon,ghc,vect,vectsize,comm)
4025 
4026 
4027 !This section has been created automatically by the script Abilint (TD).
4028 !Do not modify the following lines by hand.
4029 #undef ABI_FUNC
4030 #define ABI_FUNC 'cg_precon_block'
4031 !End of the abilint section
4032 
4033  implicit none
4034 
4035 !Arguments ------------------------------------
4036 !scalars
4037  integer,intent(in) :: blocksize,iterationnumber,npw,nspinor,optekin
4038  integer,intent(in) :: optpcon,vectsize,me_g0,comm
4039 !arrays
4040  real(dp),intent(in) :: cg(vectsize,blocksize),eval(blocksize,blocksize)
4041  real(dp),intent(in) :: ghc(vectsize,blocksize),kinpw(npw)
4042  real(dp),intent(inout) :: pcon(npw,blocksize),vect(vectsize,blocksize)
4043 
4044 !Local variables-------------------------------
4045 !scalars
4046  integer :: iblocksize,ierr,ig,igs,ipw1,ispinor
4047  real(dp) :: fac,poly,xx
4048  character(len=500) :: message
4049 !arrays
4050  real(dp) :: tsec(2)
4051  real(dp),allocatable :: ek0(:),ek0_inv(:)
4052 
4053 ! *************************************************************************
4054 
4055  call timab(536,1,tsec)
4056 
4057 !In this case, the Teter, Allan and Payne preconditioner is approximated:
4058 !the factor xx=Ekin(G) and no more Ekin(G)/Ekin(iband)
4059  if (optpcon==0) then
4060    do ispinor=1,nspinor
4061      igs=(ispinor-1)*npw
4062      if (me_g0 == 1) then
4063        do ig=1+igs,1+igs !g=0
4064          if (iterationnumber==1) then
4065            if(kinpw(ig-igs)<huge(0.0_dp)*1.d-11)then
4066              xx=kinpw(ig-igs)
4067 !            teter polynomial ratio
4068              poly=27._dp+xx*(18._dp+xx*(12._dp+xx*8._dp))
4069              fac=poly/(poly+16._dp*xx**4)
4070              if (optekin==1) fac=two*fac
4071              pcon(ig-igs,1)=fac
4072              do iblocksize=1,blocksize
4073                vect(ig,iblocksize)=(ghc(ig,iblocksize)-&
4074 &               eval(iblocksize,iblocksize)*cg(ig,iblocksize))*pcon(ig-igs,1)
4075              end do
4076            else
4077              pcon(ig-igs,1)=zero
4078              vect(ig,:)=0.0_dp
4079            end if
4080          else
4081            do iblocksize=1,blocksize
4082              vect(ig,iblocksize)=(ghc(ig,iblocksize)-&
4083 &             eval(iblocksize,iblocksize)*cg(ig,iblocksize))*pcon(ig-igs,1)
4084            end do
4085          end if
4086        end do
4087        do ig=2+igs,npw+igs
4088          if (iterationnumber==1) then
4089            if(kinpw(ig-igs)<huge(0.0_dp)*1.d-11)then
4090              xx=kinpw(ig-igs)
4091 !            teter polynomial ratio
4092              poly=27._dp+xx*(18._dp+xx*(12._dp+xx*8._dp))
4093              fac=poly/(poly+16._dp*xx**4)
4094              if (optekin==1) fac=two*fac
4095              pcon(ig-igs,1)=fac
4096              do iblocksize=1,blocksize
4097                vect(ig,iblocksize)=(ghc(ig,iblocksize)-&
4098 &               eval(iblocksize,iblocksize)*cg(ig,iblocksize))*pcon(ig-igs,1)
4099                vect(ig+npw-1,iblocksize)=(ghc(ig+npw-1,iblocksize)-&
4100 &               eval(iblocksize,iblocksize)*cg(ig+npw-1,iblocksize))*pcon(ig-igs,1)
4101              end do
4102            else
4103              pcon(ig-igs,1)=zero
4104              vect(ig,:)=zero
4105              vect(ig+npw-1,:)=zero
4106            end if
4107          else
4108            do iblocksize=1,blocksize
4109              vect(ig,iblocksize)=(ghc(ig,iblocksize)-&
4110 &             eval(iblocksize,iblocksize)*cg(ig,iblocksize))*pcon(ig-igs,1)
4111              vect(ig+npw-1,iblocksize)=(ghc(ig+npw-1,iblocksize)-&
4112 &             eval(iblocksize,iblocksize)*cg(ig+npw-1,iblocksize))*pcon(ig-igs,1)
4113            end do
4114          end if
4115        end do
4116      else
4117        do ig=1+igs,npw+igs
4118          if (iterationnumber==1) then
4119            if(kinpw(ig-igs)<huge(0.0_dp)*1.d-11)then
4120              xx=kinpw(ig-igs)
4121 !            teter polynomial ratio
4122              poly=27._dp+xx*(18._dp+xx*(12._dp+xx*8._dp))
4123              fac=poly/(poly+16._dp*xx**4)
4124              if (optekin==1) fac=two*fac
4125              pcon(ig-igs,1)=fac
4126              do iblocksize=1,blocksize
4127                vect(ig,iblocksize)=(ghc(ig,iblocksize)-&
4128 &               eval(iblocksize,iblocksize)*cg(ig,iblocksize))*pcon(ig-igs,1)
4129                vect(ig+npw,iblocksize)=(ghc(ig+npw,iblocksize)-&
4130 &               eval(iblocksize,iblocksize)*cg(ig+npw,iblocksize))*pcon(ig-igs,1)
4131              end do
4132            else
4133              pcon(ig-igs,:)=zero
4134              vect(ig,:)=zero
4135              vect(ig+npw,:)=zero
4136            end if
4137          else
4138            do iblocksize=1,blocksize
4139              vect(ig,iblocksize)=(ghc(ig,iblocksize)-&
4140 &             eval(iblocksize,iblocksize)*cg(ig,iblocksize))*pcon(ig-igs,1)
4141              vect(ig+npw,iblocksize)=(ghc(ig+npw,iblocksize)-&
4142 &             eval(iblocksize,iblocksize)*cg(ig+npw,iblocksize))*pcon(ig-igs,1)
4143            end do
4144          end if
4145        end do
4146      end if
4147    end do
4148 
4149  else if (optpcon>0) then
4150 !  Compute mean kinetic energy of all bands
4151    ABI_ALLOCATE(ek0,(blocksize))
4152    ABI_ALLOCATE(ek0_inv,(blocksize))
4153    if (iterationnumber==1.or.optpcon==1) then
4154      do iblocksize=1,blocksize
4155        if (me_g0 == 1)then
4156          ek0(iblocksize)=0.0_dp ; ipw1=2
4157          if(kinpw(1)<huge(0.0_dp)*1.d-11)ek0(iblocksize)=0.5_dp*kinpw(1)*cg(1,iblocksize)**2
4158          do ig=ipw1,npw
4159            if(kinpw(ig)<huge(0.0_dp)*1.d-11)then
4160              ek0(iblocksize)=ek0(iblocksize)+&
4161 &             kinpw(ig)*(cg(ig,iblocksize)**2+cg(ig+npw-1,iblocksize)**2)
4162            end if
4163          end do
4164        else
4165          ek0(iblocksize)=0.0_dp ; ipw1=1
4166          do ig=ipw1,npw
4167            if(kinpw(ig)<huge(0.0_dp)*1.d-11)then
4168              ek0(iblocksize)=ek0(iblocksize)+&
4169 &             kinpw(ig)*(cg(ig,iblocksize)**2+cg(ig+npw,iblocksize)**2)
4170            end if
4171          end do
4172        end if
4173      end do
4174 
4175      call xmpi_sum(ek0,comm,ierr)
4176 
4177      do iblocksize=1,blocksize
4178        if(ek0(iblocksize)<1.0d-10)then
4179          write(message, '(4a)' )ch10,&
4180 &         'cg_precon_block: the mean kinetic energy of a wavefunction vanishes.',ch10,&
4181 &         'it is reset to 0.1ha.'
4182          MSG_WARNING(message)
4183          ek0(iblocksize)=0.1_dp
4184        end if
4185      end do
4186      if (optekin==1) then
4187        ek0_inv(:)=2.0_dp/(3._dp*ek0(:))
4188      else
4189        ek0_inv(:)=1.0_dp/ek0(:)
4190      end if
4191    end if !iterationnumber==1.or.optpcon==1
4192 
4193 !  Carry out preconditioning
4194    do iblocksize=1,blocksize
4195      do ispinor=1,nspinor
4196        igs=(ispinor-1)*npw
4197        if (me_g0 == 1) then
4198          do ig=1+igs,1+igs !g=0
4199            if (iterationnumber==1.or.optpcon==1) then
4200              if(kinpw(ig-igs)<huge(0.0_dp)*1.d-11)then
4201                xx=kinpw(ig-igs)*ek0_inv(iblocksize)
4202 !              teter polynomial ratio
4203                poly=27._dp+xx*(18._dp+xx*(12._dp+xx*8._dp))
4204                fac=poly/(poly+16._dp*xx**4)
4205                if (optekin==1) fac=two*fac
4206                pcon(ig-igs,iblocksize)=fac
4207                vect(ig,iblocksize)=(ghc(ig,iblocksize)-&
4208 &               eval(iblocksize,iblocksize)*cg(ig,iblocksize))*fac
4209              else
4210                pcon(ig-igs,iblocksize)=zero
4211                vect(ig,iblocksize)=0.0_dp
4212              end if
4213            else
4214              vect(ig,iblocksize)=(ghc(ig,iblocksize)-&
4215 &             eval(iblocksize,iblocksize)*cg(ig,iblocksize))*pcon(ig-igs,iblocksize)
4216            end if
4217          end do
4218          do ig=2+igs,npw+igs
4219            if (iterationnumber==1.or.optpcon==1) then
4220              if(kinpw(ig-igs)<huge(0.0_dp)*1.d-11)then
4221                xx=kinpw(ig-igs)*ek0_inv(iblocksize)
4222 !              teter polynomial ratio
4223                poly=27._dp+xx*(18._dp+xx*(12._dp+xx*8._dp))
4224                fac=poly/(poly+16._dp*xx**4)
4225                if (optekin==1) fac=two*fac
4226                pcon(ig-igs,iblocksize)=fac
4227                vect(ig,iblocksize)=(ghc(ig,iblocksize)-&
4228 &               eval(iblocksize,iblocksize)*cg(ig,iblocksize))*fac
4229                vect(ig+npw-1,iblocksize)=(ghc(ig+npw-1,iblocksize)-&
4230 &               eval(iblocksize,iblocksize)*cg(ig+npw-1,iblocksize))*fac
4231              else
4232                pcon(ig-igs,iblocksize)=zero
4233                vect(ig,iblocksize)=zero
4234                vect(ig+npw-1,iblocksize)=zero
4235              end if
4236            else
4237              vect(ig,iblocksize)=(ghc(ig,iblocksize)-&
4238 &             eval(iblocksize,iblocksize)*cg(ig,iblocksize))*pcon(ig-igs,iblocksize)
4239              vect(ig+npw-1,iblocksize)=(ghc(ig+npw-1,iblocksize)-&
4240 &             eval(iblocksize,iblocksize)*cg(ig+npw-1,iblocksize))*pcon(ig-igs,iblocksize)
4241            end if
4242          end do
4243        else
4244          do ig=1+igs,npw+igs
4245            if (iterationnumber==1.or.optpcon==1) then
4246              if(kinpw(ig-igs)<huge(0.0_dp)*1.d-11)then
4247                xx=kinpw(ig-igs)*ek0_inv(iblocksize)
4248 !              teter polynomial ratio
4249                poly=27._dp+xx*(18._dp+xx*(12._dp+xx*8._dp))
4250                fac=poly/(poly+16._dp*xx**4)
4251                if (optekin==1) fac=two*fac
4252                pcon(ig-igs,iblocksize)=fac
4253                vect(ig,iblocksize)=(ghc(ig,iblocksize)-&
4254 &               eval(iblocksize,iblocksize)*cg(ig,iblocksize))*fac
4255                vect(ig+npw,iblocksize)=(ghc(ig+npw,iblocksize)-&
4256 &               eval(iblocksize,iblocksize)*cg(ig+npw,iblocksize))*fac
4257              else
4258                pcon(ig-igs,iblocksize)=zero
4259                vect(ig,iblocksize)=zero
4260                vect(ig+npw,iblocksize)=zero
4261              end if
4262            else
4263              vect(ig,iblocksize)=(ghc(ig,iblocksize)-&
4264 &             eval(iblocksize,iblocksize)*cg(ig,iblocksize))*pcon(ig-igs,iblocksize)
4265              vect(ig+npw,iblocksize)=(ghc(ig+npw,iblocksize)-&
4266 &             eval(iblocksize,iblocksize)*cg(ig+npw,iblocksize))*pcon(ig-igs,iblocksize)
4267            end if
4268          end do
4269        end if
4270      end do
4271    end do
4272    ABI_DEALLOCATE(ek0)
4273    ABI_DEALLOCATE(ek0_inv)
4274  end if !optpcon
4275 
4276  call timab(536,2,tsec)
4277 
4278 end subroutine cg_precon_block

m_cgtools/cg_real_zdotc [ Functions ]

[ Top ] [ m_cgtools ] [ Functions ]

NAME

  cg_real_zdotc

FUNCTION

   Perform a vector-vector operation defined as res = REAL (\Sigma (conjg(x)*y)) where x and y are n-element vectors.

INPUTS

  n = Specifies the number of elements in vector x and y
  x,y = Input arrays.

OUTPUT

  res=Real part of the scalar product.

PARENTS

SOURCE

742 function cg_real_zdotc(n,x,y) result(res)
743 
744 
745 !This section has been created automatically by the script Abilint (TD).
746 !Do not modify the following lines by hand.
747 #undef ABI_FUNC
748 #define ABI_FUNC 'cg_real_zdotc'
749 !End of the abilint section
750 
751  implicit none
752 
753 !Arguments ------------------------------------
754 !scalars
755  integer,intent(in) :: n
756 !arrays
757  real(dp),intent(in) :: x(2,n)
758  real(dp),intent(in) :: y(2,n)
759  real(dp) :: res
760 
761 !Local variables-------------------------------
762  real(dp),external :: ddot
763 
764 ! *************************************************************************
765 
766  res = ddot(2*n,x,1,y,1)
767 
768 end function cg_real_zdotc

m_cgtools/cg_setaug_zero [ Functions ]

[ Top ] [ m_cgtools ] [ Functions ]

NAME

  cg_setaug_zero

FUNCTION

  Set to zero all elements of the array that are not in the FFT box.

INPUTS

 nx,ny,nz=physical dimensions of the FFT box
 ldx,ldy,ldx=memory dimension of arr
 ndat=number of FFTs

 SIDE EFFECT
  arr(2,ldx,ldy,ldz*ndat)= all entries in the augmented region are set to zero

PARENTS

SOURCE

368 pure subroutine cg_setaug_zero(cplex,nx,ny,nz,ldx,ldy,ldz,ndat,arr)
369 
370 
371 !This section has been created automatically by the script Abilint (TD).
372 !Do not modify the following lines by hand.
373 #undef ABI_FUNC
374 #define ABI_FUNC 'cg_setaug_zero'
375 !End of the abilint section
376 
377  implicit none
378 
379 !Arguments ------------------------------------
380 !scalars
381  integer,intent(in) :: cplex,nx,ny,nz,ldx,ldy,ldz,ndat
382 !arrays
383  real(dp),intent(inout) :: arr(cplex,ldx,ldy,ldz*ndat)
384 
385 !Local variables-------------------------------
386  integer :: iy,iz,dat,padat
387 
388 ! *************************************************************************
389 
390  if (nx /= ldx) then
391    do iz=1,ldz*ndat
392      do iy=1,ldy
393        arr(:,nx+1:ldx,iy,iz) = zero
394      end do
395    end do
396  end if
397 
398  if (ny /= ldy) then
399    do iz=1,ldz*ndat
400      arr(:,:,ny+1:ldy,iz) = zero
401    end do
402  end if
403 
404  if (nz /= ldz) then
405    do dat=1,ndat
406      padat = ldz*(dat-1)
407      do iz=nz+1,ldz
408        arr(:,:,:,iz+padat) = zero
409      end do
410    end do
411  end if
412 
413 end subroutine cg_setaug_zero

m_cgtools/cg_setval [ Functions ]

[ Top ] [ m_cgtools ] [ Functions ]

NAME

  cg_setval

FUNCTION

  Set cg=alpha.

INPUTS

OUTPUT

PARENTS

CHILDREN

SOURCE

150 subroutine cg_setval(n,cg,alpha)
151 
152 
153 !This section has been created automatically by the script Abilint (TD).
154 !Do not modify the following lines by hand.
155 #undef ABI_FUNC
156 #define ABI_FUNC 'cg_setval'
157 !End of the abilint section
158 
159  implicit none
160 
161 !Arguments ------------------------------------
162 !scalars
163  integer,intent(in) :: n
164 !arrays
165  real(dp),optional,intent(in) :: alpha(2)
166  real(dp),intent(inout) :: cg(2,n)
167 
168 ! *************************************************************************
169 
170  if (PRESENT(alpha)) then
171 !$OMP PARALLEL
172 !$OMP WORKSHARE
173    cg(1,:)=alpha(1)
174    cg(2,:)=alpha(2)
175 !$OMP END WORKSHARE
176 !$OMP END PARALLEL
177  else
178 !$OMP PARALLEL
179 !$OMP WORKSHARE
180    cg(:,:)=zero
181 !$OMP END WORKSHARE
182 !$OMP END PARALLEL
183  end if
184 
185 end subroutine cg_setval

m_cgtools/cg_to_reim [ Functions ]

[ Top ] [ m_cgtools ] [ Functions ]

NAME

  cg_to_reim

FUNCTION

INPUTS

PARENTS

CHILDREN

SOURCE

432 subroutine cg_to_reim(npw,ndat,cg,factor,reim)
433 
434 
435 !This section has been created automatically by the script Abilint (TD).
436 !Do not modify the following lines by hand.
437 #undef ABI_FUNC
438 #define ABI_FUNC 'cg_to_reim'
439 !End of the abilint section
440 
441  implicit none
442 
443 !Arguments ------------------------------------
444 !scalars
445  integer,intent(in) :: npw,ndat
446  real(dp),intent(in) :: factor
447 !arrays
448  real(dp),intent(in) :: cg(2*npw*ndat)
449  real(dp),intent(out) :: reim(npw*ndat*2)
450 
451 ! *************************************************************************
452 
453  ! Pack real and imaginary part of the wavefunctions.
454  ! and multiply by scale factor if factor /= one. Could block but oh well
455  call dcopy(npw*ndat, cg(1), 2, reim(1), 1)
456  if (factor /= one) call dscal(npw*ndat,factor,reim(1),1)
457 
458  call dcopy(npw*ndat, cg(2), 2, reim(npw*ndat+1), 1)
459  if (factor /= one) call dscal(npw*ndat,factor,reim(npw*ndat+1),1)
460 
461 end subroutine cg_to_reim

m_cgtools/cg_tocplx [ Functions ]

[ Top ] [ m_cgtools ] [ Functions ]

NAME

  cg_tocplx

FUNCTION

  Convert a real array with (real,imag) part to complex.

INPUTS

  n = Specifies the number of elements in cg and ocplx
  cg(2*n)=Input array with real and imaginary part.

OUTPUT

  ocplx(n)=Output complex array.

PARENTS

CHILDREN

SOURCE

210 subroutine cg_tocplx(n, cg, ocplx)
211 
212 
213 !This section has been created automatically by the script Abilint (TD).
214 !Do not modify the following lines by hand.
215 #undef ABI_FUNC
216 #define ABI_FUNC 'cg_tocplx'
217 !End of the abilint section
218 
219  implicit none
220 
221 !Arguments ------------------------------------
222 !scalars
223  integer,intent(in) :: n
224 !arrays
225  real(dp),intent(in) :: cg(2*n)
226  complex(dpc),intent(out) :: ocplx(n)
227 
228 !Local variables ------------------------------
229 !scalars
230  integer :: ii,idx
231 
232 ! *************************************************************************
233 
234 !$OMP PARALLEL DO PRIVATE(ii,idx)
235  do ii=1,n
236    idx = 2*ii-1
237    ocplx(ii) = DCMPLX(cg(idx),cg(idx+1))
238  end do
239 
240 end subroutine cg_tocplx

m_cgtools/cg_vlocpsi [ Functions ]

[ Top ] [ m_cgtools ] [ Functions ]

NAME

  cg_vlocpsi

FUNCTION

  Apply the local part of the potentatil to the wavefunction in real space.

INPUTS

  nx,ny,nz=physical dimension of the FFT box.
  ldx,ldy,ldz=leading dimensions of the arrays.
  ndat=number of wavefunctions.
  cplex=  1 if vloc is real, 2 for complex
  vloc(cplex*ldx,ldy,ldz)=Local potential on the FFT box.

SIDE EFFECTS

  ur(2,ldx,ldy,ldz*ndat)=
    Input = wavefunctions in real space.
    Output= vloc |ur>

PARENTS

      m_dfti,m_fftw3

CHILDREN

SOURCE

2553 subroutine cg_vlocpsi(nx,ny,nz,ldx,ldy,ldz,ndat,cplex,vloc,ur)
2554 
2555 
2556 !This section has been created automatically by the script Abilint (TD).
2557 !Do not modify the following lines by hand.
2558 #undef ABI_FUNC
2559 #define ABI_FUNC 'cg_vlocpsi'
2560 !End of the abilint section
2561 
2562  implicit none
2563 
2564 !Arguments ------------------------------------
2565 !scalars
2566  integer,intent(in) :: nx,ny,nz,ldx,ldy,ldz,ndat,cplex
2567 !arrays
2568  real(dp),intent(in) :: vloc(cplex*ldx,ldy,ldz)
2569  real(dp),intent(inout) :: ur(2,ldx,ldy,ldz*ndat)
2570 
2571 !Local variables-------------------------------
2572 !scalars
2573  integer :: idat,ix,iy,iz,padat
2574  real(dp) :: fim,fre
2575 
2576 ! *************************************************************************
2577 
2578  if (cplex==1) then
2579    !
2580    if (ndat==1) then
2581 !$OMP PARALLEL DO
2582      do iz=1,nz
2583        do iy=1,ny
2584          do ix=1,nx
2585            ur(1,ix,iy,iz) = vloc(ix,iy,iz) * ur(1,ix,iy,iz)
2586            ur(2,ix,iy,iz) = vloc(ix,iy,iz) * ur(2,ix,iy,iz)
2587          end do
2588        end do
2589      end do
2590      !
2591    else
2592      !
2593 !$OMP PARALLEL DO PRIVATE(padat)
2594      do idat=1,ndat
2595        padat = ldz*(idat-1)
2596        do iz=1,nz
2597          do iy=1,ny
2598            do ix=1,nx
2599              ur(1,ix,iy,iz+padat) = vloc(ix,iy,iz) * ur(1,ix,iy,iz+padat)
2600              ur(2,ix,iy,iz+padat) = vloc(ix,iy,iz) * ur(2,ix,iy,iz+padat)
2601            end do
2602          end do
2603        end do
2604      end do
2605      !
2606    end if
2607    !
2608  else if (cplex==2)then
2609    !
2610    if (ndat==1) then
2611 !$OMP PARALLEL DO PRIVATE(fre,fim)
2612      do iz=1,nz
2613        do iy=1,ny
2614          do ix=1,nx
2615            fre = ur(1,ix,iy,iz)
2616            fim = ur(2,ix,iy,iz)
2617            ur(1,ix,iy,iz) = vloc(2*ix-1,iy,iz)*fre - vloc(2*ix,iy,iz)*fim
2618            ur(2,ix,iy,iz) = vloc(2*ix-1,iy,iz)*fim + vloc(2*ix,iy,iz)*fre
2619          end do
2620        end do
2621      end do
2622    else
2623 !$OMP PARALLEL DO PRIVATE(padat,fre,fim)
2624      do idat=1,ndat
2625        padat = ldz*(idat-1)
2626        do iz=1,nz
2627          do iy=1,ny
2628            do ix=1,nx
2629              fre = ur(1,ix,iy,iz+padat)
2630              fim = ur(2,ix,iy,iz+padat)
2631              ur(1,ix,iy,iz+padat) = vloc(2*ix-1,iy,iz)*fre - vloc(2*ix,iy,iz)*fim
2632              ur(2,ix,iy,iz+padat) = vloc(2*ix-1,iy,iz)*fim + vloc(2*ix,iy,iz)*fre
2633            end do
2634          end do
2635        end do
2636      end do
2637    end if
2638    !
2639  else
2640    ur = huge(one)
2641    !MSG_BUG("Wrong cplex")
2642  end if
2643 
2644 end subroutine cg_vlocpsi

m_cgtools/cg_zaxpby [ Functions ]

[ Top ] [ m_cgtools ] [ Functions ]

NAME

  cg_zaxpby

FUNCTION

  Scales two vectors, adds them to one another and stores result in the vector.
  y := a*x + b*y

INPUTS

 n = the number of elements in vectors x and y.
 a = Specifies the scalar a.
 x = Array.
 b = Specifies the scalar b.
 y = Array

OUTPUT

 y Contains the updated vector y.

PARENTS

CHILDREN

SOURCE

920 subroutine cg_zaxpby(n,a,x,b,y)
921 
922 
923 !This section has been created automatically by the script Abilint (TD).
924 !Do not modify the following lines by hand.
925 #undef ABI_FUNC
926 #define ABI_FUNC 'cg_zaxpby'
927 !End of the abilint section
928 
929  implicit none
930 
931 !Arguments ------------------------------------
932 !scalars
933  integer,intent(in) :: n
934  real(dp),intent(in) :: a(2),b(2)
935 !arrays
936  real(dp),intent(in) :: x(2*n)
937  real(dp),intent(inout) :: y(2*n)
938 
939 ! *************************************************************************
940 
941 #ifdef HAVE_LINALG_AXPBY
942  call zaxpby(n, a, x, 1, b, y, 1)
943 #else
944  call zscal(n, b, y, 1)
945  call zaxpy(n, a, x, 1, y,1)
946 #endif
947 
948 end subroutine cg_zaxpby

m_cgtools/cg_zaxpy [ Functions ]

[ Top ] [ m_cgtools ] [ Functions ]

NAME

  cg_zaxpy

FUNCTION

  Computes y = alpha*x + y

INPUTS

  n = Specifies the number of elements in vectors x and y.
  alpha = Specifies the scalar alpha.
  x = Array

SIDE EFFECTS

  y = Array. In output, y contains the updated vector.

PARENTS

      cgwf,dfpt_cgwf,lapackprof,rf2_init

CHILDREN

SOURCE

861 subroutine cg_zaxpy(n,alpha,x,y)
862 
863 
864 !This section has been created automatically by the script Abilint (TD).
865 !Do not modify the following lines by hand.
866 #undef ABI_FUNC
867 #define ABI_FUNC 'cg_zaxpy'
868 !End of the abilint section
869 
870  implicit none
871 
872 !Arguments ------------------------------------
873 !scalars
874  integer,intent(in) :: n
875  real(dp),intent(in) :: alpha(2)
876 !arrays
877  real(dp),intent(in) :: x(2*n)
878  real(dp),intent(inout) :: y(2*n)
879 
880 !local variables
881 ! integer :: ii
882 
883 ! *************************************************************************
884 
885  if (alpha(2) == zero) then
886    call daxpy(2*n,alpha(1),x,1,y,1)
887  else
888    call zaxpy(n,alpha,x,1,y,1)
889  end if
890 
891 end subroutine cg_zaxpy

m_cgtools/cg_zcopy [ Functions ]

[ Top ] [ m_cgtools ] [ Functions ]

NAME

  cg_zcopy

FUNCTION

  Perform y = x, where x and y are vectors.

INPUTS

  n = Specifies the number of elements in vectors x and y.
  x = Input Array

OUTPUT

  y = In output, y contains a copy of the values of x.

PARENTS

      cgwf,corrmetalwf1,dfpt_cgwf,dfpt_mkrho,dfpt_vtowfk,lapackprof,m_iowf

CHILDREN

SOURCE

534 subroutine cg_zcopy(n, x, y)
535 
536 
537 !This section has been created automatically by the script Abilint (TD).
538 !Do not modify the following lines by hand.
539 #undef ABI_FUNC
540 #define ABI_FUNC 'cg_zcopy'
541 !End of the abilint section
542 
543  implicit none
544 
545 !Arguments ------------------------------------
546 !scalars
547  integer,intent(in) :: n
548 !arrays
549  real(dp),intent(in) :: x(2*n)
550  real(dp),intent(out) :: y(2*n)
551 
552 ! *************************************************************************
553 
554  call zcopy(n,x,1,y,1)
555 
556 end subroutine cg_zcopy

m_cgtools/cg_zdotc [ Functions ]

[ Top ] [ m_cgtools ] [ Functions ]

NAME

  cg_zdotc

FUNCTION

   Perform a vector-vector operation defined as res = \Sigma (conjg(x)*y) where x and y are n-element vectors.

INPUTS

  n = Specifies the number of elements in vector x and y
  x,y = Input arrays.

OUTPUT

  res(2)=Real and Imaginary part of the scalar product.

PARENTS

SOURCE

676 function cg_zdotc(n,x,y) result(res)
677 
678 
679 !This section has been created automatically by the script Abilint (TD).
680 !Do not modify the following lines by hand.
681 #undef ABI_FUNC
682 #define ABI_FUNC 'cg_zdotc'
683 !End of the abilint section
684 
685  implicit none
686 
687 !Arguments ------------------------------------
688 !scalars
689  integer,intent(in) :: n
690 !arrays
691  real(dp),intent(in) :: x(2,n)
692  real(dp),intent(in) :: y(2,n)
693  real(dp) :: res(2)
694 
695 !Local variables-------------------------------
696 #ifdef HAVE_LINALG_ZDOTC_BUG
697  integer :: ii
698 #endif
699  complex(dpc) :: cres
700  complex(dpc),external :: zdotc
701 
702 ! *************************************************************************
703 
704 #ifdef HAVE_LINALG_ZDOTC_BUG
705  ! Workaround for veclib on MacOSx
706  res = zero
707 !$OMP PARALLEL DO PRIVATE(ii) REDUCTION(+:res)
708  do ii=1,n
709    res(1) = res(1) + x(1,ii)*y(1,ii) + x(2,ii)*y(2,ii)
710    res(2) = res(2) + x(1,ii)*y(2,ii) - x(2,ii)*y(1,ii)
711  end do
712 
713 #else
714  cres = zdotc(n, x, 1, y, 1)
715  res(1) = REAL(cres)
716  res(2) = AIMAG(cres)
717 #endif
718 
719 end function cg_zdotc

m_cgtools/cg_zdotu [ Functions ]

[ Top ] [ m_cgtools ] [ Functions ]

NAME

  cg_zdotu

FUNCTION

   Perform a vector-vector operation defined as res = \Sigma (x*y) where x and y are n-element vectors.
   Note that x is unconjugated.

INPUTS

  n = Specifies the number of elements in vector x and y
  x,y = Input arrays.

OUTPUT

  res(2)=Real and Imaginary part of the scalar product.

PARENTS

SOURCE

792 function cg_zdotu(n, x, y) result(res)
793 
794 
795 !This section has been created automatically by the script Abilint (TD).
796 !Do not modify the following lines by hand.
797 #undef ABI_FUNC
798 #define ABI_FUNC 'cg_zdotu'
799 !End of the abilint section
800 
801  implicit none
802 
803 !Arguments ------------------------------------
804 !scalars
805  integer,intent(in) :: n
806 !arrays
807  real(dp),intent(in) :: x(2,n)
808  real(dp),intent(in) :: y(2,n)
809  real(dp) :: res(2)
810 
811 !Local variables-------------------------------
812 #ifdef HAVE_LINALG_ZDOTU_BUG
813  integer :: ii
814 #endif
815  complex(dpc) :: cres
816  complex(dpc),external :: zdotu
817 
818 ! *************************************************************************
819 
820 #ifdef HAVE_LINALG_ZDOTU_BUG
821  ! Workaround for veclib on MacOSx
822  res = zero
823 !$OMP PARALLEL DO PRIVATE(ii) REDUCTION(+:res)
824  do ii=1,n
825    res(1) = res(1) + x(1,ii)*y(1,ii) - x(2,ii)*y(2,ii)
826    res(2) = res(2) + x(1,ii)*y(2,ii) + x(2,ii)*y(1,ii)
827  end do
828 #else
829  cres = zdotu(n, x, 1, y, 1)
830  res(1) = REAL(cres)
831  res(2) = AIMAG(cres)
832 #endif
833 
834 end function cg_zdotu

m_cgtools/cg_zgemm [ Functions ]

[ Top ] [ m_cgtools ] [ Functions ]

NAME

  cg_zgemm

FUNCTION

  The ?gemm routines perform a matrix-matrix operation with general matrices.
  The operation is defined as C := alpha*op(A)*op(B) + beta*C,
  where:

  op(x) is one of op(x) = x, or op(x) = x', or op(x) = conjg(x'),

  alpha and beta are scalars,
  A, B and C are matrices:
  op(A) is an m-by-k matrix,
  op(B) is a k-by-n matrix,
  C is an m-by-n matrix.

INPUTS

OUTPUT

PARENTS

      lapackprof,m_cgtools

CHILDREN

SOURCE

1061 subroutine cg_zgemm(transa,transb,npws,ncola,ncolb,cg_a,cg_b,cg_c,alpha,beta)
1062 
1063 
1064 !This section has been created automatically by the script Abilint (TD).
1065 !Do not modify the following lines by hand.
1066 #undef ABI_FUNC
1067 #define ABI_FUNC 'cg_zgemm'
1068 !End of the abilint section
1069 
1070  implicit none
1071 
1072 !Arguments ------------------------------------
1073 !scalars
1074  integer,intent(in) :: npws,ncola,ncolb
1075  real(dp),optional,intent(in) :: alpha(2),beta(2)
1076  character(len=1),intent(in) :: transa,transb
1077 !arrays
1078  real(dp),intent(in) :: cg_a(2,npws*ncola)
1079  real(dp),intent(in) :: cg_b(2,npws*ncolb)
1080  real(dp),intent(inout) :: cg_c(2,*)
1081 
1082 !Local variables-------------------------------
1083 !scalars
1084  integer :: mm,nn,kk,lda,ldb,ldc
1085  real(dp) :: my_alpha(2),my_beta(2)
1086 
1087 ! *************************************************************************
1088 
1089  lda = npws
1090  ldb = npws
1091 
1092  mm  = npws
1093  nn  = ncolb
1094  kk  = ncola
1095 
1096  if (toupper(transa) /= 'N') then
1097    mm = ncola
1098    kk = npws
1099  end if
1100  if (toupper(transb) /= 'N') nn = npws
1101 
1102  ldc = mm
1103 
1104  my_alpha = cg_cone;  if (PRESENT(alpha)) my_alpha = alpha
1105  my_beta  = cg_czero; if (PRESENT(beta))  my_beta  = beta
1106 
1107  call ZGEMM(transa,transb,mm,nn,kk,my_alpha,cg_a,lda,cg_b,ldb,my_beta,cg_c,ldc)
1108 
1109 end subroutine cg_zgemm

m_cgtools/cg_zgemv [ Functions ]

[ Top ] [ m_cgtools ] [ Functions ]

NAME

  cg_zgemv

FUNCTION

 The ?gemv routines perform a matrix-vector operation defined as

 y := alpha*A*x + beta*y,
 or
 y := alpha*A'*x + beta*y,
 or
 y := alpha*conjg(A')*x + beta*y,

 where: alpha and beta are scalars, x and y are vectors, A is an m-by-n matrix.

INPUTS

OUTPUT

PARENTS

      lapackprof,m_cgtools

CHILDREN

SOURCE

 979 subroutine cg_zgemv(trans,nrows,ncols,cgmat,vec,matvec,alpha,beta)
 980 
 981 
 982 !This section has been created automatically by the script Abilint (TD).
 983 !Do not modify the following lines by hand.
 984 #undef ABI_FUNC
 985 #define ABI_FUNC 'cg_zgemv'
 986 !End of the abilint section
 987 
 988  implicit none
 989 
 990 !Arguments ------------------------------------
 991 !scalars
 992  integer,intent(in) :: nrows,ncols
 993  real(dp),optional,intent(in) :: alpha(2),beta(2)
 994  character(len=1),intent(in) :: trans
 995 !arrays
 996  real(dp),intent(in) :: cgmat(2,nrows*ncols)
 997  real(dp),intent(in) :: vec(2,*)
 998  real(dp),intent(inout) :: matvec(2,*)
 999 
1000 !Local variables-------------------------------
1001 !scalars
1002  integer :: mm,nn,kk,lda,ldb,ldc
1003  real(dp) :: my_alpha(2),my_beta(2)
1004 
1005 ! *************************************************************************
1006 
1007  lda = nrows
1008  mm  = nrows
1009  nn  = 1
1010  kk  = ncols
1011 
1012  if (toupper(trans) /= 'N') then
1013    mm = ncols
1014    kk = nrows
1015  end if
1016 
1017  ldb = kk
1018  ldc = mm
1019 
1020  my_alpha = cg_cone;  if (PRESENT(alpha)) my_alpha = alpha
1021  my_beta  = cg_czero; if (PRESENT(beta))  my_beta  = beta
1022 
1023  call ZGEMM(trans,"N",mm,nn,kk,my_alpha,cgmat,lda,vec,ldb,my_beta,matvec,ldc)
1024  ! ZGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
1025 
1026  !call ZGEMV(trans,mm,nn,my_alpha,cgmat,lda,vec,1,my_beta,matvec,1)
1027 
1028 end subroutine cg_zgemv

m_cgtools/cg_zscal [ Functions ]

[ Top ] [ m_cgtools ] [ Functions ]

NAME

  cg_zscal

FUNCTION

  Perform x = a*x

INPUTS

  n = Specifies the number of elements in vector x.
  a(2)= The scalar a. If a(2) is zero, x = a*x is computed via zdscal

OUTPUT

  x = Updated vector.

PARENTS

      cgwf,m_cgtools

CHILDREN

SOURCE

582 subroutine cg_zscal(n, a, x)
583 
584 
585 !This section has been created automatically by the script Abilint (TD).
586 !Do not modify the following lines by hand.
587 #undef ABI_FUNC
588 #define ABI_FUNC 'cg_zscal'
589 !End of the abilint section
590 
591  implicit none
592 
593 !Arguments ------------------------------------
594 !scalars
595  integer,intent(in) :: n
596  real(dp),intent(in) :: a(2)
597 !arrays
598  real(dp),intent(inout) :: x(2*n)
599 
600 ! *************************************************************************
601 
602  if (a(2) == zero) then
603    call zdscal(n, a, x, 1)
604  else
605    call zscal(n, a, x, 1)
606  end if
607 
608 end subroutine cg_zscal

m_cgtools/cgnc_cholesky [ Functions ]

[ Top ] [ m_cgtools ] [ Functions ]

NAME

  cgnc_cholesky

FUNCTION

  Cholesky orthonormalization of the vectors stored in cgblock (version optimized for NC wavefunctions).

INPUTS

  npws=Size of each vector (usually npw*nspinor)
  nband=Number of band in cgblock
  istwfk=Storage mode for the wavefunctions. 1 for standard full mode
  me_g0=1 if this node has G=0.
  comm_pw=MPI communicator for the planewave group. Set to xmpi_comm_self for sequential mode.

SIDE EFFECTS

  cgblock(2*npws*nband)
    input: Input set of vectors.
    output: Orthonormalized set.

PARENTS

      pw_orthon

CHILDREN

SOURCE

2675 subroutine cgnc_cholesky(npws,nband,cgblock,istwfk,me_g0,comm_pw,use_gemm)
2676 
2677 
2678 !This section has been created automatically by the script Abilint (TD).
2679 !Do not modify the following lines by hand.
2680 #undef ABI_FUNC
2681 #define ABI_FUNC 'cgnc_cholesky'
2682 !End of the abilint section
2683 
2684  implicit none
2685 
2686 !Arguments ------------------------------------
2687 !scalars
2688  integer,intent(in) :: npws,nband,istwfk
2689  integer,intent(in) :: comm_pw,me_g0
2690  logical,optional,intent(in) :: use_gemm
2691 !arrays
2692  real(dp),intent(inout) :: cgblock(2*npws*nband)
2693 
2694 !Local variables ------------------------------
2695 !scalars
2696  integer :: ierr,b1,b2
2697 #ifdef DEBUG_MODE
2698  integer :: ptr
2699 #endif
2700  logical :: my_usegemm
2701  character(len=500) :: msg
2702 !arrays
2703  real(dp) :: rcg0(nband)
2704  real(dp),allocatable :: rovlp(:,:)
2705  complex(dpc),allocatable :: cf_ovlp(:,:)
2706 
2707 ! *************************************************************************
2708 
2709 #ifdef DEBUG_MODE
2710  if (istwfk==2) then
2711    ierr = 0
2712    do b1=1,nband
2713      ptr = 2 + 2*(b1-1)*npws
2714      if (ABS(cgblock(ptr)) > zero) then
2715        ierr = ierr + 1
2716        write(msg,'(a,i0,es13.6)')" Input b1, Im u(g=0) should be zero ",b1,cgblock(ptr)
2717        call wrtout(std_out,msg,"COLL")
2718        !cgblock(ptr) = zero
2719      end if
2720    end do
2721    ABI_CHECK(ierr==0,"Non zero imag part")
2722  end if
2723 #endif
2724 
2725  my_usegemm=.FALSE.; if (PRESENT(use_gemm)) my_usegemm = use_gemm
2726 
2727  ABI_MALLOC(cf_ovlp,(nband,nband))
2728  !
2729  ! 1) Calculate O_ij = <phi_i|phi_j>
2730  if (my_usegemm) then
2731    call ABI_ZGEMM("Conjugate","Normal",nband,nband,npws,cone,cgblock,npws,cgblock,npws,czero,cf_ovlp,nband)
2732  else
2733    call ZHERK("U","C",nband,npws,one,cgblock,npws,zero,cf_ovlp,nband)
2734  end if
2735 
2736  if (istwfk==1) then
2737    !
2738    ! Sum the overlap if PW are distributed.
2739    if (comm_pw /= xmpi_comm_self) then
2740      call xmpi_sum(cf_ovlp,comm_pw,ierr)
2741    end if
2742    !
2743    ! 2) Cholesky factorization: O = U^H U with U upper triangle matrix.
2744    call ZPOTRF('U',nband,cf_ovlp,nband,ierr)
2745 
2746    if (ierr/=0)  then
2747      write(msg,'(a,i0)')' ZPOTRF returned info = ',ierr
2748      MSG_ERROR(msg)
2749    end if
2750 
2751  else
2752    ! overlap is real. Note that nspinor is always 1 in this case.
2753    ABI_MALLOC(rovlp,(nband,nband))
2754    rovlp = two * REAL(cf_ovlp)
2755 
2756    if (istwfk==2 .and. me_g0==1) then
2757      ! Extract the real part at G=0 and subtract its contribution to the overlap.
2758      call dcopy(nband,cgblock,2*npws,rcg0,1)
2759      do b2=1,nband
2760        do b1=1,b2
2761          rovlp(b1,b2) = rovlp(b1,b2) - rcg0(b1)*rcg0(b2)
2762        end do
2763      end do
2764    end if
2765 
2766    ! Sum the overlap if PW are distributed.
2767    if (comm_pw /= xmpi_comm_self) then
2768      call xmpi_sum(rovlp,comm_pw,ierr)
2769    end if
2770    !
2771    ! 2) Cholesky factorization: O = U^H U with U upper triangle matrix.
2772    call DPOTRF('U',nband,rovlp,nband,ierr)
2773 
2774    if (ierr/=0)  then
2775      write(msg,'(a,i0)')' DPOTRF returned info = ',ierr
2776      MSG_ERROR(msg)
2777    end if
2778 
2779    cf_ovlp = DCMPLX(rovlp)
2780    ABI_FREE(rovlp)
2781  end if
2782  !
2783  ! 3) Solve X U = cgblock. On exit cgblock is orthonormalized.
2784  call ZTRSM('Right','Upper','Normal','Normal',npws,nband,cone,cf_ovlp,nband,cgblock,npws)
2785 
2786 #ifdef DEBUG_MODE
2787  if (istwfk==2) then
2788    ierr = 0
2789    do b1=1,nband
2790      ptr = 2 + 2*(b1-1)*npws
2791      if (ABS(cgblock(ptr)) > zero) then
2792        ierr = ierr + 1
2793        write(msg,'(a,i0,es13.6)')" Output b1, Im u(g=0) should be zero ",b1,cgblock(ptr)
2794      end if
2795    end do
2796    ABI_CHECK(ierr==0,"Non zero imag part")
2797  end if
2798 #endif
2799 
2800  ABI_FREE(cf_ovlp)
2801 
2802 end subroutine cgnc_cholesky

m_cgtools/cgnc_gramschmidt [ Functions ]

[ Top ] [ m_cgtools ] [ Functions ]

NAME

  cgnc_grortho

FUNCTION

  Gram-Schmidt orthonormalization of the vectors stored in cgblock

INPUTS

  npws=Size of each vector (usually npw*nspinor)
  nband=Number of band in cgblock
  istwfk=Storage mode for the wavefunctions. 1 for standard full mode
  me_g0=1 if this node has G=0.
  comm_pw=MPI communicator for the planewave group. Set to xmpi_comm_self for sequential mode.

SIDE EFFECTS

  cgblock(2*npws*nband)
    input: Input set of vectors.
    output: Orthonormalized set.

PARENTS

SOURCE

3147 subroutine cgnc_gramschmidt(npws,nband,cgblock,istwfk,me_g0,comm_pw)
3148 
3149 
3150 !This section has been created automatically by the script Abilint (TD).
3151 !Do not modify the following lines by hand.
3152 #undef ABI_FUNC
3153 #define ABI_FUNC 'cgnc_gramschmidt'
3154 !End of the abilint section
3155 
3156  implicit none
3157 
3158 !Arguments ------------------------------------
3159 !scalars
3160  integer,intent(in) :: npws,nband,istwfk
3161  integer,intent(in) :: comm_pw,me_g0
3162 !arrays
3163  real(dp),intent(inout) :: cgblock(2*npws*nband)
3164 
3165 !Local variables ------------------------------
3166 !scalars
3167  integer :: b1,nb2,opt
3168  logical :: normalize
3169 
3170 ! *************************************************************************
3171 
3172  ! Normalize the first vector.
3173  call cgnc_normalize(npws,1,cgblock(1),istwfk,me_g0,comm_pw)
3174  if (nband == 1) RETURN
3175 
3176  ! Orthogonaluze b1 wrt to the bands in [1,b1-1].
3177  normalize = .TRUE.
3178  do b1=2,nband
3179     opt = 1 + 2*npws*(b1-1)
3180     nb2=b1-1
3181     call cgnc_gsortho(npws,nb2,cgblock(1),1,cgblock(opt),istwfk,normalize,me_g0,comm_pw)
3182  end do
3183 
3184 end subroutine cgnc_gramschmidt

m_cgtools/cgnc_gsortho [ Functions ]

[ Top ] [ m_cgtools ] [ Functions ]

NAME

  cgnc_gsortho

FUNCTION

INPUTS

  npws=Size of each vector (usually npw*nspinor)
  nband1=Number of vectors in icg1
  nband1=Number of vectors in cg2
  comm_pw=MPI communicator.

SIDE EFFECTS

  cg2(2*npws*nband2)
  icg1(2*npws*nband1)
    input: Input set of vectors.
    output: Orthonormalized set.

PARENTS

      m_cgtools

CHILDREN

SOURCE

3050 subroutine cgnc_gsortho(npws,nband1,icg1,nband2,iocg2,istwfk,normalize,me_g0,comm_pw)
3051 
3052 
3053 !This section has been created automatically by the script Abilint (TD).
3054 !Do not modify the following lines by hand.
3055 #undef ABI_FUNC
3056 #define ABI_FUNC 'cgnc_gsortho'
3057 !End of the abilint section
3058 
3059  implicit none
3060 
3061 !Arguments ------------------------------------
3062 !scalars
3063  integer,intent(in) :: npws,nband1,nband2,istwfk,me_g0
3064  integer,optional,intent(in) :: comm_pw
3065  logical,intent(in) :: normalize
3066 !arrays
3067  real(dp),intent(in) :: icg1(2*npws*nband1)
3068  real(dp),intent(inout) :: iocg2(2*npws*nband2)
3069 
3070 !Local variables ------------------------------
3071 !scalars
3072  integer :: ierr,b1,b2
3073 !arrays
3074  real(dp) :: r_icg1(nband1),r_iocg2(nband2)
3075  real(dp),allocatable :: proj(:,:,:)
3076 
3077 ! *************************************************************************
3078 
3079  ABI_MALLOC(proj,(2,nband1,nband2))
3080  !proj = zero
3081 
3082  ! 1) Calculate <cg1|cg2>
3083  call cg_zgemm("C","N",npws,nband1,nband2,icg1,iocg2,proj)
3084 
3085  if (istwfk>1) then
3086    ! nspinor is always 1 in this case.
3087    ! Account for the missing G and set the imaginary part to zero since wavefunctions are real.
3088    proj(1,:,:) = two * proj(1,:,:)
3089    proj(2,:,:) = zero
3090    !
3091    if (istwfk==2 .and. me_g0==1) then
3092      ! Extract the real part at G=0 and subtract its contribution.
3093      call dcopy(nband1,icg1, 2*npws,r_icg1, 1)
3094      call dcopy(nband2,iocg2,2*npws,r_iocg2,1)
3095      do b2=1,nband2
3096        do b1=1,nband1
3097          proj(1,b1,b2) = proj(1,b1,b2) - r_icg1(b1) * r_iocg2(b2)
3098        end do
3099      end do
3100    end if
3101    !
3102  end if
3103  !
3104  ! This is for the MPI version
3105  if (comm_pw /= xmpi_comm_self) then
3106    call xmpi_sum(proj,comm_pw,ierr)
3107  end if
3108 
3109  ! 2) cg2 = cg2 - <cg1|cg2> cg1
3110  call cg_zgemm("N","N",npws,nband1,nband2,icg1,proj,iocg2,alpha=-cg_cone,beta=cg_cone)
3111 
3112  ABI_FREE(proj)
3113 
3114  ! 3) Normalize iocg2 if required.
3115  if (normalize) then
3116    call cgnc_normalize(npws,nband2,iocg2,istwfk,me_g0,comm_pw)
3117  end if
3118 
3119 end subroutine cgnc_gsortho

m_cgtools/cgnc_normalize [ Functions ]

[ Top ] [ m_cgtools ] [ Functions ]

NAME

  cgnc_normalize

FUNCTION

INPUTS

  npws=Size of each vector (usually npw*nspinor)
  nband=Number of vectors in icg1

SIDE EFFECTS

PARENTS

      m_cgtools

CHILDREN

SOURCE

2950 subroutine cgnc_normalize(npws,nband,cg,istwfk,me_g0,comm_pw)
2951 
2952 
2953 !This section has been created automatically by the script Abilint (TD).
2954 !Do not modify the following lines by hand.
2955 #undef ABI_FUNC
2956 #define ABI_FUNC 'cgnc_normalize'
2957 !End of the abilint section
2958 
2959  implicit none
2960 
2961 !Arguments ------------------------------------
2962 !scalars
2963  integer,intent(in) :: npws,nband,istwfk,me_g0,comm_pw
2964 !arrays
2965  real(dp),intent(inout) :: cg(2*npws*nband)
2966 
2967 !Local variables ------------------------------
2968 !scalars
2969  integer :: ptr,ierr,band
2970  character(len=500) :: msg
2971 !arrays
2972  real(dp) :: norm(nband),alpha(2)
2973 
2974 ! *************************************************************************
2975 
2976 !$OMP PARALLEL DO PRIVATE(ptr) IF (nband > 1)
2977  do band=1,nband
2978    ptr = 1 + 2*npws*(band-1)
2979    norm(band) = cg_dznrm2(npws, cg(ptr))
2980    norm(band) = norm(band) ** 2
2981    !norm(band) = cg_real_zdotc(npws, cg(ptr), cg(ptr))
2982  end do
2983 
2984  if (istwfk>1) then
2985    norm = two * norm
2986    if (istwfk==2 .and. me_g0==1) then
2987 !$OMP PARALLEL DO PRIVATE(ptr) IF (nband >1)
2988      do band=1,nband
2989        ptr = 1 + 2*npws*(band-1)
2990        norm(band) = norm(band) - cg(ptr)**2
2991      end do
2992    end if
2993  end if
2994 
2995  if (comm_pw /= xmpi_comm_self) then
2996    call xmpi_sum(norm,comm_pw,ierr)
2997  end if
2998 
2999  ierr = 0
3000  do band=1,nband
3001    if (norm(band) > zero) then
3002      norm(band) = SQRT(norm(band))
3003    else
3004      ierr = ierr + 1
3005    end if
3006  end do
3007 
3008  if (ierr/=0) then
3009    write(msg,'(a,i0,a)')" Found ",ierr," vectors with norm <= zero!"
3010    MSG_ERROR(msg)
3011  end if
3012 
3013 !$OMP PARALLEL DO PRIVATE(ptr,alpha) IF (nband > 1)
3014  do band=1,nband
3015    ptr = 1 + 2*npws*(band-1)
3016    alpha = (/one/norm(band), zero/)
3017    call cg_zscal(npws,alpha,cg(ptr))
3018  end do
3019 
3020 end subroutine cgnc_normalize

m_cgtools/cgpaw_cholesky [ Functions ]

[ Top ] [ m_cgtools ] [ Functions ]

NAME

  cgpaw_cholesky

FUNCTION

  Cholesky orthonormalization of the vectors stored in cgblock. (version for PAW wavefunctions).

INPUTS

  npws=Size of each vector (usually npw*nspinor)
  nband=Number of band in cgblock and gsc
  istwfk=Storage mode for the wavefunctions. 1 for standard full mode
  me_g0=1 if this node has G=0.
  comm_pw=MPI communicator for the planewave group. Set to xmpi_comm_self for sequential mode.

SIDE EFFECTS

  cgblock(2*npws*nband)
    input: Input set of vectors |C>, S|C>
    output: Orthonormalized set such as  <C|S|C> = 1
  gsc(2*npws*nband)
   destroyed in output.

PARENTS

      pw_orthon

CHILDREN

SOURCE

2835 subroutine cgpaw_cholesky(npws,nband,cgblock,gsc,istwfk,me_g0,comm_pw)
2836 
2837 
2838 !This section has been created automatically by the script Abilint (TD).
2839 !Do not modify the following lines by hand.
2840 #undef ABI_FUNC
2841 #define ABI_FUNC 'cgpaw_cholesky'
2842 !End of the abilint section
2843 
2844  implicit none
2845 
2846 !Arguments ------------------------------------
2847 !scalars
2848  integer,intent(in) :: npws,nband,istwfk
2849  integer,intent(in) :: me_g0,comm_pw
2850 !arrays
2851  real(dp),intent(inout) :: cgblock(2*npws*nband)
2852  real(dp),intent(inout) :: gsc(2*npws*nband)
2853 
2854 !Local variables ------------------------------
2855 !scalars
2856  integer :: ierr,b1,b2
2857  character(len=500) :: msg
2858 !arrays
2859  real(dp),allocatable :: rovlp(:,:)
2860  real(dp) :: rcg0(nband),rg0sc(nband)
2861  complex(dpc),allocatable :: cf_ovlp(:,:)
2862 
2863 ! *************************************************************************
2864 
2865  ! 1) Calculate O_ij =  <phi_i|S|phi_j>
2866  ABI_MALLOC(cf_ovlp,(nband,nband))
2867 
2868  call ABI_ZGEMM("C","N",nband,nband,npws,cone,cgblock,npws,gsc,npws,czero,cf_ovlp,nband)
2869 
2870  if (istwfk==1) then
2871    !
2872    ! Sum the overlap if PW are distributed.
2873    if (comm_pw /= xmpi_comm_self) then
2874      call xmpi_sum(cf_ovlp,comm_pw,ierr)
2875    end if
2876    !
2877    ! 2) Cholesky factorization: O = U^H U with U upper triangle matrix.
2878    call ZPOTRF('U',nband,cf_ovlp,nband,ierr)
2879 
2880    if (ierr/=0)  then
2881      write(msg,'(a,i0)')' ZPOTRF returned info= ',ierr
2882      MSG_ERROR(msg)
2883    end if
2884 
2885  else
2886    ! overlap is real. Note that nspinor is always 1 in this case.
2887    ABI_MALLOC(rovlp,(nband,nband))
2888    rovlp = two * REAL(cf_ovlp)
2889 
2890    if (istwfk==2 .and. me_g0==1) then
2891      ! Extract the real part at G=0 and subtract its contribution to the overlap.
2892      call dcopy(nband,cgblock,2*npws,rcg0,1)
2893      call dcopy(nband,gsc,2*npws,rg0sc,1)
2894      do b2=1,nband
2895        do b1=1,b2
2896         rovlp(b1,b2) = rovlp(b1,b2) - rcg0(b1)*rg0sc(b2)
2897        end do
2898      end do
2899    end if
2900    !
2901    ! Sum the overlap if PW are distributed.
2902    if (comm_pw /= xmpi_comm_self) then
2903      call xmpi_sum(rovlp,comm_pw,ierr)
2904    end if
2905   !
2906    ! 2) Cholesky factorization: O = U^H U with U upper triangle matrix.
2907    call DPOTRF('U',nband,rovlp,nband,ierr)
2908 
2909    if (ierr/=0)  then
2910      write(msg,'(a,i0)')' DPOTRF returned info= ',ierr
2911      MSG_ERROR(msg)
2912    end if
2913 
2914    cf_ovlp = DCMPLX(rovlp)
2915    ABI_FREE(rovlp)
2916  end if
2917  !
2918  ! 3) Solve X U = cgblock.
2919  call ZTRSM('Right','Upper','Normal','Normal',npws,nband,cone,cf_ovlp,nband,cgblock,npws)
2920 
2921  ! 4) Solve Y U = gsc. On exit <cgblock|gsc> = 1
2922  call ZTRSM('Right','Upper','Normal','Normal',npws,nband,cone,cf_ovlp,nband,gsc,npws)
2923 
2924  ABI_FREE(cf_ovlp)
2925 
2926 end subroutine cgpaw_cholesky

m_cgtools/cgpaw_gramschmidt [ Functions ]

[ Top ] [ m_cgtools ] [ Functions ]

NAME

  cgpaw_grortho

FUNCTION

  Gram-Schmidt orthonormalization of the vectors stored in cg

INPUTS

  npws=Size of each vector (usually npw*nspinor)
  nband=Number of bands in cg
  istwfk=Storage mode for the wavefunctions. 1 for standard full mode
  me_g0=1 if this node has G=0.
  comm_pw=MPI communicator for the planewave group. Set to xmpi_comm_self for sequential mode.

SIDE EFFECTS

  cg(2*npws*nband), gsc(2*npws*nband)
    input: Input set of vectors.
    output: Orthonormalized set.

PARENTS

SOURCE

3423 subroutine cgpaw_gramschmidt(npws,nband,cg,gsc,istwfk,me_g0,comm_pw)
3424 
3425 
3426 !This section has been created automatically by the script Abilint (TD).
3427 !Do not modify the following lines by hand.
3428 #undef ABI_FUNC
3429 #define ABI_FUNC 'cgpaw_gramschmidt'
3430 !End of the abilint section
3431 
3432  implicit none
3433 
3434 !Arguments ------------------------------------
3435 !scalars
3436  integer,intent(in) :: npws,nband,istwfk,comm_pw,me_g0
3437 !arrays
3438  real(dp),intent(inout) :: cg(2*npws*nband),gsc(2*npws*nband)
3439 
3440 !Local variables ------------------------------
3441 !scalars
3442  integer :: b1,nb2,opt
3443  logical :: normalize
3444 
3445 ! *************************************************************************
3446 
3447  ! Normalize the first vector.
3448  call cgpaw_normalize(npws,1,cg(1),gsc(1),istwfk,me_g0,comm_pw)
3449  if (nband == 1) RETURN
3450 
3451  ! Orthogonalize b1 wrt to the bands in [1,b1-1].
3452  normalize = .TRUE.
3453  do b1=2,nband
3454     opt = 1 + 2*npws*(b1-1)
3455     nb2=b1-1
3456     call cgpaw_gsortho(npws,nb2,cg(1),gsc(1),1,cg(opt),gsc(opt),istwfk,normalize,me_g0,comm_pw)
3457  end do
3458 
3459 end subroutine cgpaw_gramschmidt

m_cgtools/cgpaw_gsortho [ Functions ]

[ Top ] [ m_cgtools ] [ Functions ]

NAME

  cgpaw_gsortho

FUNCTION

  This routine uses the Gram-Schmidt method to orthogonalize a set of PAW wavefunctions.
  with respect to an input block of states.

INPUTS

  npws=Size of each vector (usually npw*nspinor)
  nband1=Number of vectors in the input block icg1
  icg1(2*npws*nband1)=Input block of vectors.
  igsc1(2*npws*nband1)= S|C> for C in icg1.
  nband2=Number of vectors to orthogonalize
  normalize=True if output wavefunction must be normalized.
  istwfk=Storage mode for the wavefunctions. 1 for standard full mode
  me_g0=1 if this node has G=0.
  comm_pw=MPI communicator for the planewave group. Set to xmpi_comm_self for sequential mode.

SIDE EFFECTS

  iocg2(2*npws*nband2), iogsc2(2*npws*nband1)
    input: set of |C> and S|C> wher |C> is the set of states to orthogonalize
    output: Orthonormalized set.

PARENTS

      m_cgtools

CHILDREN

SOURCE

3324 subroutine cgpaw_gsortho(npws,nband1,icg1,igsc1,nband2,iocg2,iogsc2,istwfk,normalize,me_g0,comm_pw)
3325 
3326 
3327 !This section has been created automatically by the script Abilint (TD).
3328 !Do not modify the following lines by hand.
3329 #undef ABI_FUNC
3330 #define ABI_FUNC 'cgpaw_gsortho'
3331 !End of the abilint section
3332 
3333  implicit none
3334 
3335 !Arguments ------------------------------------
3336 !scalars
3337  integer,intent(in) :: npws,nband1,nband2,istwfk,me_g0
3338  integer,optional,intent(in) :: comm_pw
3339  logical,intent(in) :: normalize
3340 !arrays
3341  real(dp),intent(in) :: icg1(2*npws*nband1),igsc1(2*npws*nband1)
3342  real(dp),intent(inout) :: iocg2(2*npws*nband2),iogsc2(2*npws*nband2)
3343 
3344 !Local variables ------------------------------
3345 !scalars
3346  integer :: ierr,b1,b2
3347 !arrays
3348  real(dp) :: r_icg1(nband1),r_iocg2(nband2)
3349  real(dp),allocatable :: proj(:,:,:)
3350 
3351 ! *************************************************************************
3352 
3353  ABI_MALLOC(proj,(2,nband1,nband2))
3354 
3355  ! 1) Calculate <cg1|cg2>
3356  call cg_zgemm("C","N",npws,nband1,nband2,igsc1,iocg2,proj)
3357 
3358  if (istwfk>1) then
3359    ! nspinor is always 1 in this case.
3360    ! Account for the missing G and set the imaginary part to zero since wavefunctions are real.
3361    proj(1,:,:) = two * proj(1,:,:)
3362    proj(2,:,:) = zero
3363    !
3364    if (istwfk==2 .and. me_g0==1) then
3365      ! Extract the real part at G=0 and subtract its contribution.
3366      call dcopy(nband1,igsc1,2*npws,r_icg1, 1)
3367      call dcopy(nband2,iocg2,2*npws,r_iocg2,1)
3368      do b2=1,nband2
3369        do b1=1,nband1
3370          proj(1,b1,b2) = proj(1,b1,b2) - r_icg1(b1) * r_iocg2(b2)
3371        end do
3372      end do
3373    end if
3374    !
3375  end if
3376  !
3377  ! This is for the MPI version
3378  if (comm_pw /= xmpi_comm_self) then
3379    call xmpi_sum(proj,comm_pw,ierr)
3380  end if
3381 
3382  ! 2)
3383  !   cg2 = cg2 - <Scg1|cg2> cg1
3384  ! S cg2 = S cg2 - <Scg1|cg2> S cg1
3385  call cg_zgemm("N","N",npws,nband1,nband2,icg1,proj,iocg2,alpha=-cg_cone,beta=cg_cone)
3386  call cg_zgemm("N","N",npws,nband1,nband2,igsc1,proj,iogsc2,alpha=-cg_cone,beta=cg_cone)
3387 
3388  ABI_FREE(proj)
3389 
3390  ! 3) Normalize iocg2 and iogsc2 if required.
3391  if (normalize) then
3392    call cgpaw_normalize(npws,nband2,iocg2,iogsc2,istwfk,me_g0,comm_pw)
3393  end if
3394 
3395 end subroutine cgpaw_gsortho

m_cgtools/cgpaw_normalize [ Functions ]

[ Top ] [ m_cgtools ] [ Functions ]

NAME

  cgpaw_normalize

FUNCTION

  Normalize a set of PAW pseudo wavefunctions.

INPUTS

  npws=Size of each vector (usually npw*nspinor)
  nband=Number of band in cgblock and gsc
  istwfk=Storage mode for the wavefunctions. 1 for standard full mode
  me_g0=1 if this node has G=0.
  comm_pw=MPI communicator for the planewave group. Set to xmpi_comm_self for sequential mode.

SIDE EFFECTS

  cg(2*npws*nband)
    input: Input set of vectors |C>
    output: Normalized set such as  <C|S|C> = 1
  gsc(2*npws*nband)
    input: Input set of vectors S|C>
    output: New S|C> compute with the new |C>

PARENTS

      m_cgtools

CHILDREN

SOURCE

3218 subroutine cgpaw_normalize(npws,nband,cg,gsc,istwfk,me_g0,comm_pw)
3219 
3220 
3221 !This section has been created automatically by the script Abilint (TD).
3222 !Do not modify the following lines by hand.
3223 #undef ABI_FUNC
3224 #define ABI_FUNC 'cgpaw_normalize'
3225 !End of the abilint section
3226 
3227  implicit none
3228 
3229 !Arguments ------------------------------------
3230 !scalars
3231  integer,intent(in) :: npws,nband,istwfk,me_g0,comm_pw
3232 !arrays
3233  real(dp),intent(inout) :: cg(2*npws*nband),gsc(2*npws*nband)
3234 
3235 !Local variables ------------------------------
3236 !scalars
3237  integer :: ptr,ierr,band
3238  character(len=500) :: msg
3239 !arrays
3240  real(dp) :: norm(nband),alpha(2)
3241 
3242 ! *************************************************************************
3243 
3244 !$OMP PARALLEL DO PRIVATE(ptr) IF (nband > 1)
3245  do band=1,nband
3246    ptr = 1 + 2*npws*(band-1)
3247    norm(band) = cg_real_zdotc(npws,gsc(ptr),cg(ptr))
3248  end do
3249 
3250  if (istwfk>1) then
3251    norm = two * norm
3252    if (istwfk==2 .and. me_g0==1) then
3253 !$OMP PARALLEL DO PRIVATE(ptr) IF (nband > 1)
3254      do band=1,nband
3255        ptr = 1 + 2*npws*(band-1)
3256        norm(band) = norm(band) - gsc(ptr) * cg(ptr)
3257      end do
3258    end if
3259  end if
3260 
3261  if (comm_pw /= xmpi_comm_self) then
3262    call xmpi_sum(norm,comm_pw,ierr)
3263  end if
3264 
3265  ierr = 0
3266  do band=1,nband
3267    if (norm(band) > zero) then
3268      norm(band) = SQRT(norm(band))
3269    else
3270      ierr = ierr + 1
3271    end if
3272  end do
3273 
3274  if (ierr/=0) then
3275    write(msg,'(a,i0,a)')" Found ",ierr," vectors with norm <= zero!"
3276    MSG_ERROR(msg)
3277  end if
3278 
3279  ! Scale |C> and S|C>.
3280 !$OMP PARALLEL DO PRIVATE(ptr,alpha) IF (nband > 1)
3281  do band=1,nband
3282    ptr = 1 + 2*npws*(band-1)
3283    alpha = (/one/norm(band), zero/)
3284    call cg_zscal(npws,alpha,cg(ptr))
3285    call cg_zscal(npws,alpha,gsc(ptr))
3286  end do
3287 
3288 end subroutine cgpaw_normalize

m_cgtools/cz_zprecon_block [ Functions ]

[ Top ] [ m_cgtools ] [ Functions ]

NAME

 cg_zprecon_block

FUNCTION

 precondition $<G|(H-e_{n,k})|C_{n,k}>$ for a block of band (band-FFT parallelisation)

INPUTS

  blocksize= size of blocks of bands
  $cg(vectsize,blocksize)=<G|C_{n,k}> for a block of bands$.
  $eval(blocksize,blocksize)=current block of bands eigenvalues=<C_{n,k}|H|C_{n,k}>$.
  $ghc(vectsize,blocksize)=<G|H|C_{n,k}> for a block of bands$.
  iterationnumber=number of iterative minimizations in LOBPCG
  kinpw(npw)=(modified) kinetic energy for each plane wave (Hartree)
  nspinor=number of spinorial components of the wavefunctions (on current proc)
  $vect(vectsize,blocksize)=<G|H|C_{n,k}> for a block of bands$.
  npw=number of planewaves at this k point.
  optekin= 1 if the kinetic energy used in preconditionning is modified
             according to Kresse, Furthmuller, PRB 54, 11169 (1996) [[cite:Kresse1996]]
           0 otherwise
  optpcon= 0 the TPA preconditionning matrix does not depend on band
           1 the TPA preconditionning matrix (not modified)
           2 the TPA preconditionning matrix is independant of iterationnumber
  vectsize= size of vectors
  comm=MPI communicator.

OUTPUT

  vect(2,npw)=<g|(h-eval)|c_{n,k}>*(polynomial ratio)

SIDE EFFECTS

  pcon(npw,blocksize)=preconditionning matrix
            input  if optpcon=0,2 and iterationnumber/=1
            output if optpcon=0,2 and iterationnumber==1

PARENTS

      m_lobpcg

CHILDREN

SOURCE

4323 subroutine cg_zprecon_block(cg,eval,blocksize,iterationnumber,kinpw,&
4324 &  npw,nspinor,optekin,optpcon,pcon,ghc,vect,vectsize,comm)
4325 
4326 
4327 !This section has been created automatically by the script Abilint (TD).
4328 !Do not modify the following lines by hand.
4329 #undef ABI_FUNC
4330 #define ABI_FUNC 'cg_zprecon_block'
4331 !End of the abilint section
4332 
4333  implicit none
4334 
4335 !Arguments ------------------------------------
4336 !scalars
4337  integer,intent(in) :: blocksize,iterationnumber,npw,nspinor,optekin
4338  integer,intent(in) :: optpcon,vectsize,comm
4339 !arrays
4340  real(dp),intent(in) :: kinpw(npw)
4341  real(dp),intent(inout) :: pcon(npw,blocksize)
4342  complex(dpc),intent(in) :: cg(vectsize,blocksize),eval(blocksize,blocksize)
4343  complex(dpc),intent(in) :: ghc(vectsize,blocksize)
4344  complex(dpc),intent(inout) :: vect(vectsize,blocksize)
4345 
4346 !Local variables-------------------------------
4347 !scalars
4348  integer :: iblocksize,ierr,ig,igs,ispinor
4349  real(dp) :: fac,poly,xx
4350  !character(len=500) :: message
4351 !arrays
4352  real(dp) :: tsec(2)
4353  real(dp),allocatable :: ek0(:),ek0_inv(:)
4354 
4355 ! *************************************************************************
4356 
4357  call timab(536,1,tsec)
4358 
4359 !In this case, the Teter, Allan and Payne preconditioner is approximated:
4360 !the factor xx=Ekin(G) and no more Ekin(G)/Ekin(iband)
4361  if (optpcon==0) then
4362    do ispinor=1,nspinor
4363      igs=(ispinor-1)*npw
4364      do ig=1+igs,npw+igs
4365        if (iterationnumber==1) then
4366          if(kinpw(ig-igs)<huge(0.0_dp)*1.d-11)then
4367            xx=kinpw(ig-igs)
4368 !          teter polynomial ratio
4369            poly=27._dp+xx*(18._dp+xx*(12._dp+xx*8._dp))
4370            fac=poly/(poly+16._dp*xx**4)
4371            if (optekin==1) fac=two*fac
4372            pcon(ig-igs,1)=fac
4373            do iblocksize=1,blocksize
4374              vect(ig,iblocksize)=(ghc(ig,iblocksize)-eval(iblocksize,iblocksize)*cg(ig,iblocksize))*pcon(ig-igs,1)
4375            end do
4376          else
4377            pcon(ig-igs,1)=zero
4378            vect(ig,:)=dcmplx(0.0_dp,0.0_dp)
4379          end if
4380        else
4381          do iblocksize=1,blocksize
4382            vect(ig,iblocksize)=(ghc(ig,iblocksize)-eval(iblocksize,iblocksize)*cg(ig,iblocksize))*pcon(ig-igs,1)
4383          end do
4384        end if
4385      end do
4386    end do
4387 
4388  else if (optpcon>0) then
4389 !  Compute mean kinetic energy of all bands
4390    ABI_ALLOCATE(ek0,(blocksize))
4391    ABI_ALLOCATE(ek0_inv,(blocksize))
4392    if (iterationnumber==1.or.optpcon==1) then
4393      do iblocksize=1,blocksize
4394        ek0(iblocksize)=0.0_dp
4395        do ispinor=1,nspinor
4396          igs=(ispinor-1)*npw
4397          do ig=1+igs,npw+igs
4398            if(kinpw(ig-igs)<huge(0.0_dp)*1.d-11)then
4399              ek0(iblocksize)=ek0(iblocksize)+kinpw(ig-igs)*&
4400 &             (real(cg(ig,iblocksize))**2+aimag(cg(ig,iblocksize))**2)
4401            end if
4402          end do
4403        end do
4404      end do
4405 
4406      call xmpi_sum(ek0,comm,ierr)
4407 
4408      do iblocksize=1,blocksize
4409        if(ek0(iblocksize)<1.0d-10)then
4410          MSG_WARNING('the mean kinetic energy of a wavefunction vanishes. it is reset to 0.1ha.')
4411          ek0(iblocksize)=0.1_dp
4412        end if
4413      end do
4414      if (optekin==1) then
4415        ek0_inv(:)=2.0_dp/(3._dp*ek0(:))
4416      else
4417        ek0_inv(:)=1.0_dp/ek0(:)
4418      end if
4419    end if !iterationnumber==1.or.optpcon==1
4420 
4421 !  Carry out preconditioning
4422    do iblocksize=1,blocksize
4423      do ispinor=1,nspinor
4424        igs=(ispinor-1)*npw
4425        do ig=1+igs,npw+igs
4426          if (iterationnumber==1.or.optpcon==1) then
4427            if(kinpw(ig-igs)<huge(0.0_dp)*1.d-11)then
4428              xx=kinpw(ig-igs)*ek0_inv(iblocksize)
4429 !            teter polynomial ratio
4430              poly=27._dp+xx*(18._dp+xx*(12._dp+xx*8._dp))
4431              fac=poly/(poly+16._dp*xx**4)
4432              if (optekin==1) fac=two*fac
4433              pcon(ig-igs,iblocksize)=fac
4434              vect(ig,iblocksize)=(ghc(ig,iblocksize)-&
4435 &             eval(iblocksize,iblocksize)*cg(ig,iblocksize))*pcon(ig-igs,iblocksize)
4436            else
4437              pcon(ig-igs,iblocksize)=zero
4438              vect(ig,iblocksize)=dcmplx(0.0_dp,0.0_dp)
4439            end if
4440          else
4441            vect(ig,iblocksize)=(ghc(ig,iblocksize)-&
4442 &           eval(iblocksize,iblocksize)*cg(ig,iblocksize))*pcon(ig-igs,iblocksize)
4443          end if
4444        end do
4445      end do
4446    end do
4447    ABI_DEALLOCATE(ek0)
4448    ABI_DEALLOCATE(ek0_inv)
4449  end if !optpcon
4450 
4451  call timab(536,2,tsec)
4452 
4453 end subroutine cg_zprecon_block

m_cgtools/dotprod_g [ Functions ]

[ Top ] [ m_cgtools ] [ Functions ]

NAME

 dotprod_g

FUNCTION

 Compute scalar product <vec1|vect2> of complex vectors vect1 and vect2 (can be the same)
 Take into account the storage mode of the vectors (istwf_k)
 If option=1, compute only real part, if option=2 compute also imaginary part.
 If the number of calls to the dot product scales quadratically
 with the volume of the system, it is preferable not to
 call the present routine, but but to write a specially
 optimized routine, that will avoid many branches related to
 the existence of 'option' and 'istwf_k'.

INPUTS

  istwf_k=option parameter that describes the storage of wfs
  vect1(2,npw)=first vector (one should take its complex conjugate)
  vect2(2,npw)=second vector
  npw= (effective) number of planewaves at this k point (including spinorial level)
  option= 1 if only real part to be computed,
          2 if both real and imaginary.
          3 if in case istwf_k==1 must compute real and imaginary parts,
               but if  istwf_k >1 must compute only real part
  me_g0=1 if this processor treats G=0, 0 otherwise
  comm=MPI communicator used to reduce the results.

OUTPUT

  $doti=\Im ( <vect1|vect2> )$ , output only if option=2 and eventually option=3.
  $dotr=\Re ( <vect1|vect2> )$

PARENTS

      cgwf,chebfi,corrmetalwf1,d2frnl,dfpt_cgwf,dfpt_nsteltwf,dfpt_nstpaw
      dfpt_nstwf,dfpt_vtowfk,dfpt_wfkfermi,dfptnl_resp,dotprod_set_cgcprj
      dotprodm_sumdiag_cgcprj,eig2stern,extrapwf,fock2ACE,fock_ACE_getghc
      fock_getghc,m_efmas,m_gkk,m_phgamma,m_phpi,m_rf2,m_sigmaph,mkresi
      nonlop_gpu,nonlop_test,rf2_init

CHILDREN

SOURCE

1308 subroutine dotprod_g(dotr,doti,istwf_k,npw,option,vect1,vect2,me_g0,comm)
1309 
1310 
1311 !This section has been created automatically by the script Abilint (TD).
1312 !Do not modify the following lines by hand.
1313 #undef ABI_FUNC
1314 #define ABI_FUNC 'dotprod_g'
1315 !End of the abilint section
1316 
1317  implicit none
1318 
1319 !Arguments ------------------------------------
1320 !scalars
1321  integer,intent(in) :: istwf_k,npw,option,me_g0,comm
1322  real(dp),intent(out) :: doti,dotr
1323 !arrays
1324  real(dp),intent(in) :: vect1(2,npw),vect2(2,npw)
1325 
1326 !Local variables-------------------------------
1327 !scalars
1328  integer :: ierr
1329 !arrays
1330  real(dp) :: dotarr(2)
1331 
1332 ! *************************************************************************
1333 
1334  if (istwf_k==1) then ! General k-point
1335 
1336    if(option==1)then
1337      dotr = cg_real_zdotc(npw,vect1,vect2)
1338    else
1339      dotarr = cg_zdotc(npw,vect1,vect2)
1340      dotr = dotarr(1)
1341      doti = dotarr(2)
1342    end if
1343 
1344  else if (istwf_k==2 .and. me_g0==1) then ! Gamma k-point and I have G=0
1345 
1346    dotr=half*vect1(1,1)*vect2(1,1)
1347    dotr = dotr + cg_real_zdotc(npw-1,vect1(1,2),vect2(1,2))
1348    dotr = two*dotr
1349    if (option==2) doti=zero
1350 
1351  else ! Other TR k-points
1352 
1353    dotr = cg_real_zdotc(npw,vect1,vect2)
1354    dotr=two*dotr
1355    if (option==2) doti=zero
1356  end if
1357 
1358 !Reduction in case of parallelism
1359  if (xmpi_comm_size(comm)>1) then
1360    if (option==1.or.istwf_k/=1) then
1361      call xmpi_sum(dotr,comm,ierr)
1362    else
1363      dotarr(1)=dotr ; dotarr(2)=doti
1364      call xmpi_sum(dotarr,comm,ierr)
1365      dotr=dotarr(1) ; doti=dotarr(2)
1366    end if
1367  end if
1368 
1369 end subroutine dotprod_g

m_cgtools/dotprod_v [ Functions ]

[ Top ] [ m_cgtools ] [ Functions ]

NAME

 dotprod_v

FUNCTION

 Compute dot product of two potentials (integral over FFT grid), to obtain
 a square residual-like quantity (so the sum of product of values
 is NOT divided by the number of FFT points, and NOT multiplied by the primitive cell volume).
 Take into account the spin components of the potentials (nspden), and sum over them.

INPUTS

  cplex=if 1, real space functions on FFT grid are REAL, if 2, COMPLEX
  nfft= (effective) number of FFT grid points (for this processor)
  nspden=number of spin-density components
  opt_storage: 0, if potentials are stored as V^up-up, V^dn-dn, Re[V^up-dn], Im[V^up-dn]
               1, if potentials are stored as V, B_x, B_y, Bz  (B=magn. field)
  pot1(cplex*nfft,nspden)=first real space potential on FFT grid
  pot2(cplex*nfft,nspden)=second real space potential on FFT grid
  comm=MPI communicator in which results will be reduced.

OUTPUT

  dotr= value of the dot product

PARENTS

      m_epjdos

CHILDREN

SOURCE

1571 subroutine dotprod_v(cplex,dotr,nfft,nspden,opt_storage,pot1,pot2,comm)
1572 
1573 
1574 !This section has been created automatically by the script Abilint (TD).
1575 !Do not modify the following lines by hand.
1576 #undef ABI_FUNC
1577 #define ABI_FUNC 'dotprod_v'
1578 !End of the abilint section
1579 
1580  implicit none
1581 
1582 !Arguments ------------------------------------
1583 !scalars
1584  integer,intent(in) :: cplex,nfft,nspden,opt_storage,comm
1585  real(dp),intent(out) :: dotr
1586 !arrays
1587  real(dp),intent(in) :: pot1(cplex*nfft,nspden),pot2(cplex*nfft,nspden)
1588 
1589 !Local variables-------------------------------
1590 !scalars
1591  integer :: ierr,ifft,ispden
1592  real(dp) :: ar
1593 !arrays
1594 
1595 ! *************************************************************************
1596 
1597 !Real or complex inputs are coded
1598 
1599  dotr=zero
1600 !$OMP PARALLEL DO COLLAPSE(2) REDUCTION(+:dotr)
1601  do ispden=1,min(nspden,2)
1602    do ifft=1,cplex*nfft
1603      dotr =dotr + pot1(ifft,ispden)*pot2(ifft,ispden)
1604    end do
1605  end do
1606 
1607  if (nspden==4) then
1608    ar=zero
1609 !$OMP PARALLEL DO COLLAPSE(2) REDUCTION(+:ar)
1610    do ispden=3,4
1611      do ifft=1,cplex*nfft
1612        ar = ar + pot1(ifft,ispden)*pot2(ifft,ispden)
1613      end do
1614    end do
1615 
1616    if (opt_storage==0) then
1617      if (cplex==1) then
1618        dotr = dotr+two*ar
1619      else
1620        dotr = dotr+ar
1621      end if
1622    else
1623      dotr = half*(dotr+ar)
1624    end if
1625  end if
1626 
1627 !MPIWF reduction (addition) on dotr is needed here
1628  if (xmpi_comm_size(comm)>1) then
1629    call xmpi_sum(dotr,comm,ierr)
1630  end if
1631 
1632 end subroutine dotprod_v

m_cgtools/dotprod_vn [ Functions ]

[ Top ] [ m_cgtools ] [ Functions ]

NAME

 dotprod_vn

FUNCTION

 Compute dot product of potential and density (integral over FFT grid), to obtain
 an energy-like quantity (so the usual dotproduct is divided
 by the number of FFT points, and multiplied by the primitive cell volume).
 Take into account the spin components of the density and potentials (nspden),
 and sum correctly over them. Note that the storage of densities and
 potentials is different : for potential, one stores the matrix components,
 while for the density, one stores the trace, and then, either the
 up component (if nspden=2), or the magnetization vector (if nspden=4).

INPUTS

  cplex=if 1, real space functions on FFT grid are REAL, if 2, COMPLEX
  dens(cplex*nfft,nspden)=real space density on FFT grid
  mpi_enreg=information about MPI parallelization
  nfft= (effective) number of FFT grid points (for this processor)
  nfftot= total number of FFT grid points
  nspden=number of spin-density components
  option= if 1, only the real part is computed
          if 2, both real and imaginary parts are computed  (not yet coded)
  pot(cplex*nfft,nspden)=real space potential on FFT grid
                 (will be complex conjugated if cplex=2 and option=2)
  ucvol=unit cell volume (Bohr**3)

OUTPUT

  doti= imaginary part of the dot product, output only if option=2 (and cplex=2).
  dotr= real part

NOTES

  Concerning storage when nspden=4:
   cplex=1:
     V is stored as : V^11, V^22, Re[V^12], Im[V^12] (complex, hermitian)
     N is stored as : n, m_x, m_y, m_z               (real)
   cplex=2:
     V is stored as : V^11, V^22, V^12, i.V^21 (complex)
     N is stored as : n, m_x, m_y, mz          (complex)

PARENTS

      dfpt_dyxc1,dfpt_eltfrxc,dfpt_nselt,dfpt_nstdy,dfpt_nstpaw,dfpt_rhotov
      dfptnl_loop,energy,newfermie1,odamix,prcrskerker2,rhotov,rhotoxc,setvtr

CHILDREN

      xmpi_sum

SOURCE

1684 subroutine dotprod_vn(cplex,dens,dotr,doti,nfft,nfftot,nspden,option,pot,ucvol, &
1685     mpi_comm_sphgrid)  ! Optional
1686 
1687 
1688 !This section has been created automatically by the script Abilint (TD).
1689 !Do not modify the following lines by hand.
1690 #undef ABI_FUNC
1691 #define ABI_FUNC 'dotprod_vn'
1692 !End of the abilint section
1693 
1694  implicit none
1695 
1696 !Arguments ------------------------------------
1697 !scalars
1698  integer,intent(in) :: cplex,nfft,nfftot,nspden,option
1699  integer,intent(in),optional :: mpi_comm_sphgrid
1700  real(dp),intent(in) :: ucvol
1701  real(dp),intent(out) :: doti,dotr
1702 !arrays
1703  real(dp),intent(in) :: dens(cplex*nfft,nspden),pot(cplex*nfft,nspden)
1704 
1705 !Local variables-------------------------------
1706 !scalars
1707  integer  :: ierr,ifft,jfft
1708  real(dp) :: dim11,dim12,dim21,dim22,dim_dn,dim_up,dre11,dre12,dre21,dre22
1709  real(dp) :: dre_dn,dre_up,factor,nproc_sphgrid,pim11,pim12,pim21,pim22,pim_dn,pim_up,pre11
1710  real(dp) :: pre12,pre21,pre22,pre_dn,pre_up
1711  real(dp) :: bx_re,bx_im,by_re,by_im,bz_re,bz_im,v0_re,v0_im
1712 !arrays
1713  real(dp) :: buffer2(2)
1714 
1715 ! *************************************************************************
1716 
1717 !Real or complex inputs are coded
1718  DBG_CHECK(ANY(cplex==(/1,2/)),"Wrong cplex")
1719  DBG_CHECK(ANY(nspden==(/1,2,4/)),"Wrong nspden")
1720 
1721 !Real or complex output are coded
1722  DBG_CHECK(ANY(option==(/1,2/)),"Wrong option")
1723 
1724  dotr=zero; doti=zero
1725 
1726  if(nspden==1)then
1727 
1728    if(option==1 .or. cplex==1 )then
1729 !$OMP PARALLEL DO REDUCTION(+:dotr)
1730      do ifft=1,cplex*nfft
1731        dotr=dotr + pot(ifft,1)*dens(ifft,1)
1732      end do
1733 !    dotr = ddot(cplex*nfft,pot,1,dens,1)
1734 
1735    else  ! option==2 and cplex==2 : one builds the imaginary part, from complex den/pot
1736 
1737 !$OMP PARALLEL DO PRIVATE(jfft) REDUCTION(+:dotr,doti)
1738      do ifft=1,nfft
1739        jfft=2*ifft
1740        dotr=dotr + pot(jfft-1,1)*dens(jfft-1,1) + pot(jfft,1)*dens(jfft  ,1)
1741        doti=doti + pot(jfft-1,1)*dens(jfft  ,1) - pot(jfft,1)*dens(jfft-1,1)
1742      end do
1743 
1744    end if
1745 
1746  else if(nspden==2)then
1747 
1748    if(option==1 .or. cplex==1 )then
1749 !$OMP PARALLEL DO REDUCTION(+:dotr)
1750      do ifft=1,cplex*nfft
1751        dotr=dotr + pot(ifft,1)* dens(ifft,2)     &    ! This is the spin up contribution
1752 &      + pot(ifft,2)*(dens(ifft,1)-dens(ifft,2))      ! This is the spin down contribution
1753      end do
1754 
1755    else ! option==2 and cplex==2 : one builds the imaginary part, from complex den/pot
1756 
1757 !$OMP PARALLEL DO DEFAULT(PRIVATE) SHARED(nfft,dens,pot) REDUCTION(+:dotr,doti)
1758      do ifft=1,nfft
1759 
1760        jfft=2*ifft
1761        dre_up=dens(jfft-1,2)
1762        dim_up=dens(jfft  ,2)
1763        dre_dn=dens(jfft-1,1)-dre_up
1764        dim_dn=dens(jfft  ,1)-dim_up
1765        pre_up=pot(jfft-1,1)
1766        pim_up=pot(jfft  ,1)
1767        pre_dn=pot(jfft-1,2)
1768        pim_dn=pot(jfft  ,2)
1769 
1770        dotr=dotr + pre_up * dre_up &
1771 &       + pim_up * dim_up &
1772 &       + pre_dn * dre_dn &
1773 &       + pim_dn * dim_dn
1774        doti=doti + pre_up * dim_up &
1775 &       - pim_up * dre_up &
1776 &       + pre_dn * dim_dn &
1777 &       - pim_dn * dre_dn
1778 
1779      end do
1780    end if
1781 
1782  else if(nspden==4)then
1783 !  \rho{\alpha,\beta} V^{\alpha,\beta} =
1784 !  rho*(V^{11}+V^{22})/2$
1785 !  + m_x Re(V^{12})- m_y Im{V^{12}}+ m_z(V^{11}-V^{22})/2
1786    if (cplex==1) then
1787 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(nfft,dens,pot) REDUCTION(+:dotr)
1788      do ifft=1,nfft
1789        dotr=dotr + &
1790 &       (pot(ifft,1) + pot(ifft,2))*half*dens(ifft,1) &   ! This is the density contrib
1791 &      + pot(ifft,3)      *dens(ifft,2) &   ! This is the m_x contrib
1792 &      - pot(ifft,4)      *dens(ifft,3) &   ! This is the m_y contrib
1793 &      +(pot(ifft,1) - pot(ifft,2))*half*dens(ifft,4)     ! This is the m_z contrib
1794      end do
1795    else ! cplex=2
1796 !    Note concerning storage when cplex=2:
1797 !    V is stored as : v^11, v^22, V^12, i.V^21 (each are complex)
1798 !    N is stored as : n, m_x, m_y, mZ          (each are complex)
1799      if (option==1) then
1800 !$OMP PARALLEL DO DEFAULT(PRIVATE) SHARED(nfft,dens,pot) REDUCTION(+:dotr)
1801        do ifft=1,nfft
1802          jfft=2*ifft
1803          dre11=half*(dens(jfft-1,1)+dens(jfft-1,4))
1804          dim11=half*(dens(jfft  ,1)+dens(jfft-1,4))
1805          dre22=half*(dens(jfft-1,1)-dens(jfft-1,4))
1806          dim22=half*(dens(jfft  ,1)-dens(jfft-1,4))
1807          dre12=half*(dens(jfft-1,2)+dens(jfft  ,3))
1808          dim12=half*(dens(jfft  ,2)-dens(jfft-1,3))
1809          dre21=half*(dens(jfft-1,2)-dens(jfft  ,3))
1810          dim21=half*(dens(jfft  ,2)+dens(jfft-1,3))
1811          pre11= pot(jfft-1,1)
1812          pim11= pot(jfft  ,1)
1813          pre22= pot(jfft-1,2)
1814          pim22= pot(jfft  ,2)
1815          pre12= pot(jfft-1,3)
1816          pim12= pot(jfft  ,3)
1817          pre21= pot(jfft  ,4)
1818          pim21=-pot(jfft-1,4)
1819 
1820          v0_re=half*(pre11+pre22)
1821          v0_im=half*(pim11+pim22)
1822          bx_re=half*(pre12+pre21)
1823          bx_im=half*(pim12+pim21)
1824          by_re=half*(-pim12+pim21)
1825          by_im=half*(pre12-pre21)
1826          bz_re=half*(pre11-pre22)
1827          bz_im=half*(pim11-pim22)
1828          dotr=dotr+v0_re * dens(jfft-1,1)&
1829 &         + v0_im * dens(jfft  ,1) &
1830 &         + bx_re * dens(jfft-1,2) &
1831 &         + bx_im * dens(jfft  ,2) &
1832 &         + by_re * dens(jfft-1,3) &
1833 &         + by_im * dens(jfft  ,3) &
1834 &         + bz_re * dens(jfft-1,4) &
1835 &         + bz_im * dens(jfft  ,4)
1836 !         dotr=dotr + pre11 * dre11 &
1837 !&         + pim11 * dim11 &
1838 !&         + pre22 * dre22 &
1839 !&         + pim22 * dim22 &
1840 !&         + pre12 * dre12 &
1841 !&         + pim12 * dim12 &
1842 !&         + pre21 * dre21 &
1843 !&         + pim21 * dim21
1844        end do
1845      else ! option=2
1846 !$OMP PARALLEL DO DEFAULT(PRIVATE) SHARED(nfft,dens,pot) REDUCTION(+:dotr,doti)
1847        do ifft=1,nfft
1848          jfft=2*ifft
1849          dre11=half*(dens(jfft-1,1)+dens(jfft-1,4))
1850          dim11=half*(dens(jfft  ,1)+dens(jfft-1,4))
1851          dre22=half*(dens(jfft-1,1)-dens(jfft-1,4))
1852          dim22=half*(dens(jfft  ,1)-dens(jfft-1,4))
1853          dre12=half*(dens(jfft-1,2)+dens(jfft  ,3))
1854          dim12=half*(dens(jfft  ,2)-dens(jfft-1,3))
1855          dre21=half*(dens(jfft-1,2)-dens(jfft  ,3))
1856          dim21=half*(dens(jfft  ,2)+dens(jfft-1,3))
1857          pre11= pot(jfft-1,1)
1858          pim11= pot(jfft  ,1)
1859          pre22= pot(jfft-1,2)
1860          pim22= pot(jfft  ,2)
1861          pre12= pot(jfft-1,3)
1862          pim12= pot(jfft  ,3)
1863          pre21= pot(jfft  ,4)
1864          pim21=-pot(jfft-1,4)
1865          dotr=dotr + pre11 * dre11 &
1866 &         + pim11 * dim11 &
1867 &         + pre22 * dre22 &
1868 &         + pim22 * dim22 &
1869 &         + pre12 * dre12 &
1870 &         + pim12 * dim12 &
1871 &         + pre21 * dre21 &
1872 &         + pim21 * dim21
1873          doti=doti + pre11 * dim11 &
1874 &         - pim11 * dre11 &
1875 &         + pre22 * dim22 &
1876 &         - pim22 * dre22 &
1877 &         + pre12 * dim12 &
1878 &         - pim12 * dre12 &
1879 &         + pre21 * dim21 &
1880 &         - pim21 * dre21
1881        end do
1882      end if ! option
1883    end if ! cplex
1884  end if ! nspden
1885 
1886  factor=ucvol/dble(nfftot)
1887  dotr=factor*dotr
1888  doti=factor*doti
1889 
1890 !MPIWF reduction (addition) on dotr, doti is needed here
1891  if(present(mpi_comm_sphgrid)) then
1892    nproc_sphgrid=xmpi_comm_size(mpi_comm_sphgrid)
1893    if(nproc_sphgrid>1) then
1894      buffer2(1)=dotr
1895      buffer2(2)=doti
1896      call xmpi_sum(buffer2,mpi_comm_sphgrid,ierr)
1897      dotr=buffer2(1)
1898      doti=buffer2(2)
1899    end if
1900  end if
1901 
1902 end subroutine dotprod_vn

m_cgtools/fxphas_seq [ Functions ]

[ Top ] [ m_cgtools ] [ Functions ]

NAME

 fxphas_seq

FUNCTION

 Fix phase of all bands. Keep normalization but maximize real part
 (minimize imag part). Also fix the sign of real part
 by setting the first non-zero element to be positive.

 This version has been stripped of all the mpi_enreg junk by MJV

INPUTS

  cg(2,mcg)= contains the wavefunction |c> coefficients.
  gsc(2,mgsc)= if useoverlap==1, contains the S|c> coefficients,
               where S is an overlap matrix.
  icg=shift to be applied on the location of data in the array cg
  igsc=shift to be applied on the location of data in the array gsc
  istwfk=input option parameter that describes the storage of wfs
    (set to 1 if usual complex vectors)
  mcg=size of second dimension of cg
  mgsc=size of second dimension of gsc
  nband_k=number of bands
  npw_k=number of planewaves
  useoverlap=describe the overlap of wavefunctions:
               0: no overlap (S=Identi0,ty_matrix)
               1: wavefunctions are overlapping

OUTPUT

  cg(2,mcg)=same array with altered phase.
  gsc(2,mgsc)= same array with altered phase.

NOTES

 When the sign of the real part was fixed (modif v3.1.3g.6), the
 test Tv3#5 , dataset 5, behaved differently than previously.
 This should be cleared up.

PARENTS

      m_dynmat,rayleigh_ritz

CHILDREN

SOURCE

4500 subroutine fxphas_seq(cg,gsc,icg,igsc,istwfk,mcg,mgsc,nband_k,npw_k,useoverlap)
4501 
4502 
4503 !This section has been created automatically by the script Abilint (TD).
4504 !Do not modify the following lines by hand.
4505 #undef ABI_FUNC
4506 #define ABI_FUNC 'fxphas_seq'
4507 !End of the abilint section
4508 
4509  implicit none
4510 
4511 !Arguments ------------------------------------
4512 !scalars
4513  integer,intent(in) :: icg,igsc,istwfk,mcg,mgsc,nband_k,npw_k,useoverlap
4514 !arrays
4515  real(dp),intent(inout) :: cg(2,mcg),gsc(2,mgsc*useoverlap)
4516 
4517 !Local variables-------------------------------
4518 !scalars
4519  integer :: iband,ii,indx
4520  real(dp) :: cim,cre,gscim,gscre,quotient,root1,root2,saa,sab,sbb,theta
4521  real(dp) :: thppi,xx,yy
4522  character(len=500) :: message
4523 !arrays
4524  real(dp),allocatable :: cimb(:),creb(:),saab(:),sabb(:),sbbb(:) !,sarr(:,:)
4525 
4526 ! *************************************************************************
4527 
4528 !The general case, where a complex phase indeterminacy is present
4529  if(istwfk==1)then
4530 
4531    ABI_ALLOCATE(cimb,(nband_k))
4532    ABI_ALLOCATE(creb,(nband_k))
4533    ABI_ALLOCATE(saab,(nband_k))
4534    ABI_ALLOCATE(sabb,(nband_k))
4535    ABI_ALLOCATE(sbbb,(nband_k))
4536    cimb(:)=zero ; creb(:)=zero
4537 
4538 !  Loop over bands
4539 !  TODO: MG store saa arrays in sarr(3,nband_k) to reduce false sharing.
4540    do iband=1,nband_k
4541      indx=icg+(iband-1)*npw_k
4542 
4543 !    Compute several sums over Re, Im parts of c
4544      saa=0.0_dp ; sbb=0.0_dp ; sab=0.0_dp
4545      do ii=1+indx,npw_k+indx
4546        saa=saa+cg(1,ii)*cg(1,ii)
4547        sbb=sbb+cg(2,ii)*cg(2,ii)
4548        sab=sab+cg(1,ii)*cg(2,ii)
4549      end do
4550      saab(iband)=saa
4551      sbbb(iband)=sbb
4552      sabb(iband)=sab
4553    end do ! iband
4554 
4555 
4556    do iband=1,nband_k
4557 
4558      indx=icg+(iband-1)*npw_k
4559 
4560      saa=saab(iband)
4561      sbb=sbbb(iband)
4562      sab=sabb(iband)
4563 
4564 !    Get phase angle theta
4565      if (sbb+saa>tol8)then
4566        if(abs(sbb-saa)>tol8*(sbb+saa) .or. 2*abs(sab)>tol8*(sbb+saa))then
4567          if (abs(sbb-saa)>tol8*abs(sab)) then
4568            quotient=sab/(sbb-saa)
4569            theta=0.5_dp*atan(2.0_dp*quotient)
4570          else
4571 !          Taylor expansion of the atan in terms of inverse of its argument. Correct up to 1/x2, included.
4572            theta=0.25_dp*(pi-(sbb-saa)/sab)
4573          end if
4574 !        Check roots to get theta for max Re part
4575          root1=cos(theta)**2*saa+sin(theta)**2*sbb-2.0_dp*cos(theta)*sin(theta)*sab
4576          thppi=theta+0.5_dp*pi
4577          root2=cos(thppi)**2*saa+sin(thppi)**2*sbb-2.0_dp*cos(thppi)*sin(thppi)*sab
4578          if (root2>root1) theta=thppi
4579        else
4580 !        The real part vector and the imaginary part vector are orthogonal, and of same norm. Strong indeterminacy.
4581 !        Will determine the first non-zero coefficient, and fix its phase
4582          do ii=1+indx,npw_k+indx
4583            cre=cg(1,ii)
4584            cim=cg(2,ii)
4585            if(cre**2+cim**2>tol8**2*(saa+sbb))then
4586              if(cre**2>tol8**2**cim**2)then
4587                theta=atan(cim/cre)
4588              else
4589 !              Taylor expansion of the atan in terms of inverse of its argument. Correct up to 1/x2, included.
4590                theta=pi/2-cre/cim
4591              end if
4592              exit
4593            end if
4594          end do
4595        end if
4596      else
4597        write(message,'(a,i0,5a)')&
4598 &       'The eigenvector with band ',iband,' has zero norm.',ch10,&
4599 &       'This usually happens when the number of bands (nband) is comparable to the number of planewaves (mpw)',ch10,&
4600 &       'Action: Check the parameters of the calculation. If nband ~ mpw, then decrease nband or, alternatively, increase ecut'
4601        MSG_ERROR(message)
4602      end if
4603 
4604      xx=cos(theta)
4605      yy=sin(theta)
4606 
4607 !    Here, set the first non-zero element to be positive
4608      do ii=1+indx,npw_k+indx
4609        cre=cg(1,ii)
4610        cim=cg(2,ii)
4611        cre=xx*cre-yy*cim
4612        if(abs(cre)>tol8)exit
4613      end do
4614      if(cre<zero)then
4615        xx=-xx ; yy=-yy
4616      end if
4617 
4618      creb(iband)=xx
4619      cimb(iband)=yy
4620 
4621    end do
4622 
4623    do iband=1,nband_k
4624 
4625      indx=icg+(iband-1)*npw_k
4626 
4627      xx=creb(iband)
4628      yy=cimb(iband)
4629      do ii=1+indx,npw_k+indx
4630        cre=cg(1,ii)
4631        cim=cg(2,ii)
4632        cg(1,ii)=xx*cre-yy*cim
4633        cg(2,ii)=xx*cim+yy*cre
4634      end do
4635 
4636 !    Alter phase of array S|cg>
4637      if (useoverlap==1) then
4638        indx=igsc+(iband-1)*npw_k
4639        do ii=1+indx,npw_k+indx
4640          gscre=gsc(1,ii)
4641          gscim=gsc(2,ii)
4642          gsc(1,ii)=xx*gscre-yy*gscim
4643          gsc(2,ii)=xx*gscim+yy*gscre
4644        end do
4645      end if
4646 
4647    end do ! iband
4648 
4649    ABI_DEALLOCATE(cimb)
4650    ABI_DEALLOCATE(creb)
4651    ABI_DEALLOCATE(saab)
4652    ABI_DEALLOCATE(sabb)
4653    ABI_DEALLOCATE(sbbb)
4654 
4655 !  ====================================================================
4656 
4657 !  Storages that take into account the time-reversal symmetry : the freedom is only a sign freedom
4658  else  ! if istwfk/=1
4659 
4660    ABI_ALLOCATE(creb,(nband_k))
4661    creb(:)=zero
4662 !  Loop over bands
4663    do iband=1,nband_k
4664 
4665      indx=icg+(iband-1)*npw_k
4666 
4667 !    Here, set the first non-zero real element to be positive
4668      do ii=1+indx,npw_k+indx
4669        cre=cg(1,ii)
4670        if(abs(cre)>tol8)exit
4671      end do
4672      creb(iband)=cre
4673 
4674    end do ! iband
4675 
4676    do iband=1,nband_k
4677 
4678      cre=creb(iband)
4679      if(cre<zero)then
4680        indx=icg+(iband-1)*npw_k
4681        do ii=1+indx,npw_k+indx
4682          cg(1,ii)=-cg(1,ii)
4683          cg(2,ii)=-cg(2,ii)
4684        end do
4685        if(useoverlap==1)then
4686          indx=igsc+(iband-1)*npw_k
4687          do ii=1+indx,npw_k+indx
4688            gsc(1,ii)=-gsc(1,ii)
4689            gsc(2,ii)=-gsc(2,ii)
4690          end do
4691        end if
4692      end if
4693 
4694    end do ! iband
4695 
4696    ABI_DEALLOCATE(creb)
4697 
4698  end if ! istwfk
4699 
4700 end subroutine fxphas_seq

m_cgtools/matrixelmt_g [ Functions ]

[ Top ] [ m_cgtools ] [ Functions ]

NAME

 matrixelmt_g

FUNCTION

  Compute a matrix element of two wavefunctions, in reciprocal space,
  for an operator that is diagonal in reciprocal space: <wf1|op|wf2>
  For the time being, only spin-independent operators are treated.

INPUTS

  diag(npw)=diagonal operator (real, spin-independent !)
  istwf_k=storage mode of the vectors
  needimag=0 if the imaginary part is not needed ; 1 if the imaginary part is needed
  npw=number of planewaves of the first vector
  nspinor=number of spinor components
  vect1(2,npw*nspinor)=first vector
  vect2(2,npw*nspinor)=second vector
  comm_fft=MPI communicator for the FFT
  me_g0=1 if this processors treats the G=0 component.

OUTPUT

  ai=imaginary part of the matrix element
  ar=real part of the matrix element

PARENTS

      dfpt_vtowfk

CHILDREN

SOURCE

1405 subroutine matrixelmt_g(ai,ar,diag,istwf_k,needimag,npw,nspinor,vect1,vect2,me_g0,comm_fft)
1406 
1407 
1408 !This section has been created automatically by the script Abilint (TD).
1409 !Do not modify the following lines by hand.
1410 #undef ABI_FUNC
1411 #define ABI_FUNC 'matrixelmt_g'
1412 !End of the abilint section
1413 
1414  implicit none
1415 
1416 !Arguments ------------------------------------
1417 !scalars
1418  integer,intent(in) :: istwf_k,needimag,npw,nspinor,me_g0,comm_fft
1419  real(dp),intent(out) :: ai,ar
1420 !arrays
1421  real(dp),intent(in) :: diag(npw),vect1(2,npw*nspinor),vect2(2,npw*nspinor)
1422 
1423 !Local variables-------------------------------
1424 !scalars
1425  integer :: i1,ierr,ipw
1426  character(len=500) :: message
1427 !arrays
1428  real(dp) :: buffer2(2)
1429  !real(dp),allocatable :: re_prod(:), im_prod(:)
1430 
1431 ! *************************************************************************
1432 
1433  if (nspinor==2 .and. istwf_k/=1) then
1434    write(message,'(a,a,a,i6,a,i6)')&
1435 &   'When istwf_k/=1, nspinor must be 1,',ch10,&
1436 &   'however, nspinor=',nspinor,', and istwf_k=',istwf_k
1437    MSG_BUG(message)
1438  end if
1439 
1440 #if 0
1441  !TODO
1442  ABI_MALLOC(re_prod,(npw*nspinor))
1443  do ipw=1,npw*nspinor
1444   re_prod(ipw) = vect1(1,ipw)*vect2(1,ipw) + vect1(2,ipw)*vect2(2,ipw)
1445  end do
1446 
1447  if (needimag == 1) then
1448    ABI_MALLOC(im_prod,(npw*nspinor))
1449    do ipw=1,npw*nspinor
1450      im_prod(ipw) = vect1(1,ipw)*vect2(2,ipw) - vect1(2,ipw)*vect2(1,ipw)
1451    end do
1452  end if
1453 #endif
1454 
1455  ar=zero
1456  if(needimag==1)ai=zero
1457 
1458 !Normal storage mode
1459  if(istwf_k==1)then
1460 
1461 !  Need only real part
1462    if(needimag==0)then
1463 
1464      do ipw=1,npw
1465        ar=ar+diag(ipw)*(vect1(1,ipw)*vect2(1,ipw)+vect1(2,ipw)*vect2(2,ipw))
1466      end do
1467 
1468      if(nspinor==2)then
1469        do ipw=1+npw,2*npw
1470          ar=ar+diag(ipw-npw)*(vect1(1,ipw)*vect2(1,ipw)+vect1(2,ipw)*vect2(2,ipw))
1471        end do
1472      end if
1473 
1474    else ! Need also the imaginary part
1475 
1476      do ipw=1,npw
1477        ar=ar+diag(ipw)*(vect1(1,ipw)*vect2(1,ipw)+vect1(2,ipw)*vect2(2,ipw))
1478        ai=ai+diag(ipw)*(vect1(1,ipw)*vect2(2,ipw)-vect1(2,ipw)*vect2(1,ipw))
1479      end do
1480 
1481      if(nspinor==2)then
1482        do ipw=1+npw,2*npw
1483          ar=ar+diag(ipw-npw)*(vect1(1,ipw)*vect2(1,ipw)+vect1(2,ipw)*vect2(2,ipw))
1484          ai=ai+diag(ipw-npw)*(vect1(1,ipw)*vect2(2,ipw)-vect1(2,ipw)*vect2(1,ipw))
1485        end do
1486      end if
1487 
1488    end if ! needimag
1489 
1490  else if(istwf_k>=2)then
1491 
1492 !  XG030513 : MPIWF need to know which proc has G=0
1493 
1494    i1=1
1495    if(istwf_k==2 .and. me_g0==1)then
1496      ar=half*diag(1)*vect1(1,1)*vect2(1,1) ; i1=2
1497    end if
1498 
1499 !  Need only real part
1500    if(needimag==0)then
1501 
1502      do ipw=i1,npw
1503        ar=ar+diag(ipw)*(vect1(1,ipw)*vect2(1,ipw)+vect1(2,ipw)*vect2(2,ipw))
1504      end do
1505      ar=two*ar
1506 
1507    else ! Need also the imaginary part
1508 
1509      do ipw=i1,npw
1510        ar=ar+diag(ipw)*(vect1(1,ipw)*vect2(1,ipw)+vect1(2,ipw)*vect2(2,ipw))
1511        ai=ai+diag(ipw)*(vect1(1,ipw)*vect2(2,ipw)-vect1(2,ipw)*vect2(1,ipw))
1512      end do
1513      ar=two*ar ; ai=two*ai
1514 
1515    end if
1516 
1517  end if ! istwf_k
1518 
1519 #if 0
1520  ABI_FREE(re_prod)
1521  if (needimag == 1) then
1522    ABI_FREE(im_prod)
1523  end if
1524 #endif
1525 
1526 !MPIWF need to make reduction on ar and ai .
1527  if (xmpi_comm_size(comm_fft)>1) then
1528    buffer2(1)=ai
1529    buffer2(2)=ar
1530    call xmpi_sum(buffer2,comm_fft ,ierr)
1531    ai=buffer2(1)
1532    ar=buffer2(2)
1533  end if
1534 
1535 end subroutine matrixelmt_g

m_cgtools/mean_fftr [ Functions ]

[ Top ] [ m_cgtools ] [ Functions ]

NAME

 mean_fftr

FUNCTION

  Compute the mean of an arraysp(nfft,nspden), over the FFT grid, for each component nspden,
  and return it in meansp(nspden).
  Take into account the spread of the array due to parallelism: the actual number of fft
  points is nfftot, but the number of points on this proc is nfft only.
  So : for ispden from 1 to nspden
       meansp(ispden) = sum(ifft=1,nfftot) arraysp(ifft,ispden) / nfftot

INPUTS

  arraysp(nfft,nspden)=the array whose average has to be computed
  nfft=number of FFT points stored by one proc
  nfftot=total number of FFT points
  nspden=number of spin-density components

OUTPUT

  meansp(nspden)=mean value for each nspden component

PARENTS

      fresid,newvtr,pawmknhat,prcref,prcref_PMA,psolver_rhohxc,rhohxcpositron
      rhotov,rhotoxc

CHILDREN

SOURCE

2032 subroutine mean_fftr(arraysp,meansp,nfft,nfftot,nspden,mpi_comm_sphgrid)
2033 
2034 
2035 !This section has been created automatically by the script Abilint (TD).
2036 !Do not modify the following lines by hand.
2037 #undef ABI_FUNC
2038 #define ABI_FUNC 'mean_fftr'
2039 !End of the abilint section
2040 
2041  implicit none
2042 
2043 !Arguments ------------------------------------
2044 !scalars
2045  integer,intent(in) :: nfft,nfftot,nspden
2046  integer,intent(in),optional:: mpi_comm_sphgrid
2047 !arrays
2048  real(dp),intent(in) :: arraysp(nfft,nspden)
2049  real(dp),intent(out) :: meansp(nspden)
2050 
2051 !Local variables-------------------------------
2052 !scalars
2053  integer :: ierr,ifft,ispden,nproc_sphgrid
2054  real(dp) :: invnfftot,tmean
2055 
2056 ! *************************************************************************
2057 
2058  invnfftot=one/(dble(nfftot))
2059 
2060  do ispden=1,nspden
2061    tmean=zero
2062 !$OMP PARALLEL DO REDUCTION(+:tmean)
2063    do ifft=1,nfft
2064      tmean=tmean+arraysp(ifft,ispden)
2065    end do
2066    meansp(ispden)=tmean*invnfftot
2067  end do
2068 
2069 !XG030514 : MPIWF The values of meansp(ispden) should
2070 !now be summed accross processors in the same WF group, and spread on all procs.
2071  if(present(mpi_comm_sphgrid)) then
2072    nproc_sphgrid=xmpi_comm_size(mpi_comm_sphgrid)
2073    if(nproc_sphgrid>1) then
2074      call xmpi_sum(meansp,nspden,mpi_comm_sphgrid,ierr)
2075    end if
2076  end if
2077 
2078 end subroutine mean_fftr

m_cgtools/overlap_g [ Functions ]

[ Top ] [ m_cgtools ] [ Functions ]

NAME

 overlap_g

FUNCTION

 Compute the scalar product between WF at two different k-points
 < u_{n,k1} | u_{n,k2}>

INPUTS

 mpw = maximum dimensioned size of npw
 npw_k1 = number of plane waves at k1
 npw_k2 = number of plane waves at k2
 nspinor = 1 for scalar, 2 for spinor wavefunctions
 pwind_k = array required to compute the scalar product (see initberry.f)
 vect1 = wavefunction at k1: | u_{n,k1} >
 vect2 = wavefunction at k1: | u_{n,k2} >

OUTPUT

 doti = imaginary part of the scalarproduct
 dotr = real part of the scalarproduct

NOTES

 In case a G-vector of the basis sphere of plane waves at k1
 does not belong to the basis sphere of plane waves at k2,
 pwind = 0. Therefore, the dimensions of vect1 &
 vect2 are (1:2,0:mpw) and the element (1:2,0) MUST be set to zero.

 The current implementation if not compatible with TR-symmetry (i.e. istwfk/=1) !

PARENTS

      dfptff_bec,dfptff_die,dfptff_ebp,dfptff_edie,dfptff_gbefd
      dfptff_gradberry,qmatrix,smatrix

CHILDREN

SOURCE

4740 subroutine overlap_g(doti,dotr,mpw,npw_k1,npw_k2,nspinor,pwind_k,vect1,vect2)
4741 
4742 
4743 !This section has been created automatically by the script Abilint (TD).
4744 !Do not modify the following lines by hand.
4745 #undef ABI_FUNC
4746 #define ABI_FUNC 'overlap_g'
4747 !End of the abilint section
4748 
4749  implicit none
4750 
4751 !Arguments ------------------------------------
4752 !scalars
4753  integer,intent(in) :: mpw,npw_k1,npw_k2,nspinor
4754  real(dp),intent(out) :: doti,dotr
4755 !arrays
4756  integer,intent(in) :: pwind_k(mpw)
4757  real(dp),intent(in) :: vect1(1:2,0:mpw*nspinor),vect2(1:2,0:mpw*nspinor)
4758 
4759 !Local variables-------------------------------
4760 !scalars
4761  integer :: ipw,ispinor,jpw,spnshft1,spnshft2
4762 
4763 ! *************************************************************************
4764 
4765 !Check if vect1(:,0) = 0 and vect2(:,0) = 0
4766  if ((abs(vect1(1,0)) > tol12).or.(abs(vect1(2,0)) > tol12).or. &
4767 & (abs(vect2(1,0)) > tol12).or.(abs(vect2(2,0)) > tol12)) then
4768    MSG_BUG('vect1(:,0) and/or vect2(:,0) are not equal to zero')
4769  end if
4770 
4771 !Compute the scalar product
4772  dotr = zero; doti = zero
4773  do ispinor = 1, nspinor
4774    spnshft1 = (ispinor-1)*npw_k1
4775    spnshft2 = (ispinor-1)*npw_k2
4776 !$OMP PARALLEL DO PRIVATE(jpw) REDUCTION(+:doti,dotr)
4777    do ipw = 1, npw_k1
4778      jpw = pwind_k(ipw)
4779      dotr = dotr + vect1(1,spnshft1+ipw)*vect2(1,spnshft2+jpw) + vect1(2,spnshft1+ipw)*vect2(2,spnshft2+jpw)
4780      doti = doti + vect1(1,spnshft1+ipw)*vect2(2,spnshft2+jpw) - vect1(2,spnshft1+ipw)*vect2(1,spnshft2+jpw)
4781    end do
4782  end do
4783 
4784 end subroutine overlap_g

m_cgtools/projbd [ Functions ]

[ Top ] [ m_cgtools ] [ Functions ]

NAME

 projbd

FUNCTION

 Project out vector "direc" onto the bands contained in "cg".
 if useoverlap==0
  New direc=direc-$sum_{j/=i} { <cg_{j}|direc>.|cg_{j}> }$
 if useoverlap==1 (use of overlap matrix S)
  New direc=direc-$sum_{j/=i} { <cg_{j}|S|direc>.|cg_{j}> }$
 (index i can be set to -1 to sum over all bands)

INPUTS

  cg(2,mcg)=wavefunction coefficients for ALL bands
  iband0=which particular band we are interested in ("i" in the above formula)
         Can be set to -1 to sum over all bands...
  icg=shift to be given to the location of the data in cg
  iscg=shift to be given to the location of the data in cg
  istwf_k=option parameter that describes the storage of wfs
  mcg=maximum size of second dimension of cg
  mscg=maximum size of second dimension of scg
  nband=number of bands
  npw=number of planewaves
  nspinor=number of spinorial components (on current proc)
  scg(2,mscg*useoverlap)=<G|S|band> for ALL bands,
                        where S is an overlap matrix
  scprod_io=0 if scprod array has to be computed; 1 if it is input (already in memory)
  tim_projbd=timing code of the calling subroutine(can be set to 0 if not attributed)
  useoverlap=describe the overlap of wavefunctions:
               0: no overlap (S=Identity_matrix)
               1: wavefunctions are overlapping
 me_g0=1 if this processors treats G=0, 0 otherwise.
 comm=MPI communicator (used if G vectors are distributed.

SIDE EFFECTS

  direc(2,npw)= input: vector to be orthogonalised with respect to cg (and S)
                output: vector that has been orthogonalized wrt cg (and S)

  scprod(2,nband)=scalar_product
        if useoverlap==0: scalar_product_i=$<cg_{j}|direc_{i}>$
        if useoverlap==1: scalar_product_i=$<cg_{j}|S|direc_{i}>$
    if scprod_io=0, scprod is output
    if scprod_io=1, scprod is input

NOTES

  1) MPIWF Might have to be recoded for efficient paralellism

  2) The new version employs BLAS2 routine so that the OMP parallelism is delegated to BLAS library.

  3) Note for PAW: ref.= PRB 73, 235101 (2006) [[cite:Audouze2006]], equations (71) and (72):
     in normal use, projbd applies P_c projector
     if cg and scg are inverted, projbd applies P_c+ projector

  4) cg_zgemv wraps ZGEMM whose implementation is more efficient, especially in the threaded case.

PARENTS

      cgwf,dfpt_cgwf,dfpt_nstpaw,getdc1,lapackprof

CHILDREN

SOURCE

3527 subroutine projbd(cg,direc,iband0,icg,iscg,istwf_k,mcg,mscg,nband,&
3528 &                 npw,nspinor,scg,scprod,scprod_io,tim_projbd,useoverlap,me_g0,comm)
3529 
3530 
3531 !This section has been created automatically by the script Abilint (TD).
3532 !Do not modify the following lines by hand.
3533 #undef ABI_FUNC
3534 #define ABI_FUNC 'projbd'
3535 !End of the abilint section
3536 
3537  implicit none
3538 
3539 !Arguments ------------------------------------
3540 !scalars
3541  integer,intent(in) :: iband0,icg,iscg,istwf_k,mcg,mscg,nband,npw,nspinor
3542  integer,intent(in) :: scprod_io,tim_projbd,useoverlap,me_g0,comm
3543 !arrays
3544  real(dp),intent(in) :: cg(2,mcg),scg(2,mscg*useoverlap)
3545  real(dp),intent(inout) :: direc(2,npw*nspinor)
3546  real(dp),intent(inout) :: scprod(2,nband)
3547 
3548 !Local variables-------------------------------
3549 !scalars
3550  integer :: nbandm,npw_sp,ierr
3551 !arrays
3552  real(dp) :: tsec(2),bkp_scprod(2),bkp_dirg0(2)
3553 
3554 ! *************************************************************************
3555 
3556  call timab(210+tim_projbd,1,tsec)
3557 
3558  npw_sp=npw*nspinor
3559 
3560  nbandm=nband
3561 
3562  if (istwf_k==1) then
3563 
3564    if (scprod_io==0) then
3565      if (useoverlap==1) then
3566        call cg_zgemv("C",npw_sp,nbandm,scg(1,iscg+1),direc,scprod)
3567      else
3568        call cg_zgemv("C",npw_sp,nbandm,cg(1,icg+1),  direc,scprod)
3569      end if
3570      call xmpi_sum(scprod,comm,ierr)
3571    end if
3572 
3573    if (iband0>0.and.iband0<=nbandm) then
3574      bkp_scprod = scprod(:,iband0)
3575      scprod(:,iband0) = zero
3576    end if
3577 
3578    call cg_zgemv("N",npw_sp,nbandm,cg(1,icg+1),scprod,direc,alpha=-cg_cone,beta=cg_cone)
3579 
3580    if (iband0>0.and.iband0<=nbandm) scprod(:,iband0) = bkp_scprod ! Restore previous value as scprod is output.
3581 
3582  else if (istwf_k>=2) then
3583    !
3584    !  u_{G0/2}(G) = u_{G0/2}(-G-G0)^*; k = G0/2
3585    !  hence:
3586    !  sum_G f*(G) g(G) = 2 REAL sum_G^{IZONE} w(G) f*(G)g(G)
3587    !  where
3588    !  w(G) = 1 except for k=0 and G=0 where w(G=0) = 1/2.
3589    !
3590    if (scprod_io==0) then
3591 
3592      if (useoverlap==1) then
3593 
3594        if (istwf_k==2 .and. me_g0==1) then
3595          bkp_dirg0 = direc(:,1)
3596          direc(1,1) = half * direc(1,1)
3597          direc(2,1) = zero
3598        end if
3599 
3600        call cg_zgemv("C",npw_sp,nbandm,scg(1,iscg+1),direc,scprod)
3601        scprod = two * scprod
3602        scprod(2,:) = zero
3603 
3604        if(istwf_k==2 .and. me_g0==1) direc(:,1) = bkp_dirg0
3605 
3606      else
3607 
3608        if (istwf_k==2 .and. me_g0==1) then
3609          bkp_dirg0 = direc(:,1)
3610          direc(1,1) = half * direc(1,1)
3611          direc(2,1) = zero
3612        end if
3613 
3614        call cg_zgemv("C",npw_sp,nbandm,cg(1,icg+1),direc,scprod)
3615        scprod = two * scprod
3616        scprod(2,:) = zero
3617 
3618        if (istwf_k==2 .and. me_g0==1) direc(:,1) = bkp_dirg0
3619      end if ! useoverlap
3620 
3621      call xmpi_sum(scprod,comm,ierr)
3622    end if
3623 
3624    if (iband0>0.and.iband0<=nbandm) then
3625      bkp_scprod = scprod(:,iband0)
3626      scprod(:,iband0) = zero
3627    end if
3628 
3629    call cg_zgemv("N",npw_sp,nbandm,cg(1,icg+1),scprod,direc,alpha=-cg_cone,beta=cg_cone)
3630 
3631    if (iband0>0.and.iband0<=nbandm) scprod(:,iband0) = bkp_scprod ! Restore previous value as scprod is output.
3632 
3633  end if ! Test on istwf_k
3634 
3635  call timab(210+tim_projbd,2,tsec)
3636 
3637 end subroutine projbd

m_cgtools/pw_orthon [ Functions ]

[ Top ] [ m_cgtools ] [ Functions ]

NAME

 pw_orthon

FUNCTION

 Normalize nvec complex vectors each of length nelem and then orthogonalize by modified Gram-Schmidt.
 Two orthogonality conditions are available:
  Simple orthogonality: ${<Vec_{i}|Vec_{j}>=Delta_ij}$
  Orthogonality with overlap S: ${<Vec_{i}|S|Vec_{j}>=Delta_ij}$

INPUTS

  icg=shift to be given to the location of the data in cg(=vecnm)
  igsc=shift to be given to the location of the data in gsc(=ovl_vecnm)
  istwf_k=option parameter that describes the storage of wfs
  mcg=maximum size of second dimension of cg(=vecnm)
  mgsc=maximum size of second dimension of gsc(=ovl_vecnm)
  nelem=number of complex elements in each vector
  nvec=number of vectors to be orthonormalized
  ortalgo= option for the choice of the algorithm
         -1: no orthogonalization (direct return)
          0 or 2: old algorithm (use of buffers)
          1: new algorithm (use of blas)
          3: new new algorithm (use of lapack without copy)
  useoverlap=select the orthogonality condition
               0: no overlap between vectors
               1: vectors are overlapping
  me_g0=1 if this processor has G=0, 0 otherwise
  comm=MPI communicator

SIDE EFFECTS

  vecnm= input: vectors to be orthonormalized; array of nvec column
                vectors,each of length nelem,shifted by icg
                This array is complex or else real(dp) of twice length
         output: orthonormalized set of vectors
  if (useoverlap==1) only:
    ovl_vecnm= input: product of overlap and input vectors:
                      S|vecnm>,where S is the overlap operator
               output: updated S|vecnm> according to vecnm

NOTES

 Note that each vector has an arbitrary phase which is not fixed in this routine.

 WARNING: not yet suited for nspinor=2 with istwfk/=1

PARENTS

      lapackprof,vtowfk,wfconv

CHILDREN

      abi_xcopy,abi_xorthonormalize,abi_xtrsm,cgnc_cholesky,cgnc_gramschmidt
      cgpaw_cholesky,cgpaw_gramschmidt,ortho_reim,timab,xmpi_sum

SOURCE

5102 subroutine pw_orthon(icg,igsc,istwf_k,mcg,mgsc,nelem,nvec,ortalgo,ovl_vecnm,useoverlap,vecnm,me_g0,comm)
5103 
5104  use m_abi_linalg
5105 
5106 !This section has been created automatically by the script Abilint (TD).
5107 !Do not modify the following lines by hand.
5108 #undef ABI_FUNC
5109 #define ABI_FUNC 'pw_orthon'
5110 !End of the abilint section
5111 
5112  implicit none
5113 
5114 !Arguments ------------------------------------
5115 !scalars
5116  integer,intent(in) :: icg,igsc,istwf_k,mcg,mgsc,nelem,nvec,ortalgo,useoverlap,me_g0,comm
5117 !arrays
5118  real(dp),intent(inout) :: ovl_vecnm(2,mgsc*useoverlap),vecnm(2,mcg)
5119 
5120 !Local variables-------------------------------
5121 !scalars
5122  integer :: ierr,ii,ii0,ii1,ii2,ivec,ivec2
5123  integer :: rvectsiz,vectsize,cg_idx,gsc_idx
5124  real(dp) :: doti,dotr,sum,xnorm
5125  character(len=500) :: msg
5126 !arrays
5127  integer :: cgindex(nvec),gscindex(nvec)
5128  real(dp) :: buffer2(2),tsec(2)
5129  real(dp),allocatable :: rblockvectorbx(:,:),rblockvectorx(:,:),rgramxbx(:,:)
5130  complex(dpc),allocatable :: cblockvectorbx(:,:),cblockvectorx(:,:)
5131  complex(dpc),allocatable :: cgramxbx(:,:)
5132 
5133 ! *************************************************************************
5134 
5135 #ifdef DEBUG_MODE
5136 !Make sure imaginary part at G=0 vanishes
5137  if (istwf_k==2) then
5138    do ivec=1,nvec
5139      if(abs(vecnm(2,1+nelem*(ivec-1)+icg))>zero)then
5140 !      if(abs(vecnm(2,1+nelem*(ivec-1)+icg))>tol16)then
5141        write(msg,'(2a,3i0,2es16.6,a,a)')&
5142 &       ' For istwf_k=2,observed the following element of vecnm :',ch10,&
5143 &       nelem,ivec,icg,vecnm(1:2,1+nelem*(ivec-1)+icg),ch10,&
5144 &       '  with a non-negligible imaginary part.'
5145        MSG_BUG(msg)
5146      end if
5147    end do
5148  end if
5149 #endif
5150 
5151 !Nothing to do if ortalgo=-1
5152  if(ortalgo==-1) return
5153 
5154  do ivec=1,nvec
5155    cgindex(ivec)=nelem*(ivec-1)+icg+1
5156    gscindex(ivec)=nelem*(ivec-1)+igsc+1
5157  end do
5158 
5159  if (ortalgo==3) then
5160 !  =========================
5161 !  First (new new) algorithm
5162 !  =========================
5163 !  NEW VERSION: avoid copies, use ZHERK for NC
5164    cg_idx = cgindex(1)
5165    if (useoverlap==1) then
5166      gsc_idx = gscindex(1)
5167      call cgpaw_cholesky(nelem,nvec,vecnm(1,cg_idx),ovl_vecnm(1,gsc_idx),istwf_k,me_g0,comm)
5168    else
5169      call cgnc_cholesky(nelem,nvec,vecnm(1,cg_idx),istwf_k,me_g0,comm,use_gemm=.FALSE.)
5170    end if
5171 
5172  else if(ortalgo==1) then
5173 !  =======================
5174 !  Second (new) algorithm
5175 !  =======================
5176 !  This first algorithm seems to be more efficient especially in the parallel band-FFT mode.
5177 
5178    if(istwf_k==1) then
5179      vectsize=nelem
5180      ABI_ALLOCATE(cgramxbx,(nvec,nvec))
5181      ABI_ALLOCATE(cblockvectorx,(vectsize,nvec))
5182      ABI_ALLOCATE(cblockvectorbx,(vectsize,nvec))
5183      call abi_xcopy(nvec*vectsize,vecnm(:,cgindex(1):cgindex(nvec)-1),1,cblockvectorx,1,x_cplx=2)
5184      if (useoverlap == 1) then
5185        call abi_xcopy(nvec*vectsize,ovl_vecnm(:,gscindex(1):gscindex(nvec)-1),1,cblockvectorbx,1,x_cplx=2)
5186      else
5187        call abi_xcopy(nvec*vectsize,vecnm(:,cgindex(1):cgindex(nvec)-1),1,cblockvectorbx,1,x_cplx=2)
5188      end if
5189      call abi_xorthonormalize(cblockvectorx,cblockvectorbx,nvec,comm,cgramxbx,vectsize)
5190      call abi_xcopy(nvec*vectsize,cblockvectorx,1,vecnm(:,cgindex(1):cgindex(nvec)-1),1,x_cplx=2)
5191      if (useoverlap == 1) then
5192        call abi_xtrsm('r','u','n','n',vectsize,nvec,cone,cgramxbx,nvec,cblockvectorbx,vectsize)
5193        call abi_xcopy(nvec*vectsize,cblockvectorbx,1,ovl_vecnm(:,gscindex(1):gscindex(nvec)-1),1,x_cplx=2)
5194      end if
5195      ABI_DEALLOCATE(cgramxbx)
5196      ABI_DEALLOCATE(cblockvectorx)
5197      ABI_DEALLOCATE(cblockvectorbx)
5198 
5199    else if(istwf_k==2) then
5200      ! Pack real and imaginary part of the wavefunctions.
5201      rvectsiz=nelem
5202      vectsize=2*nelem; if(me_g0==1) vectsize=vectsize-1
5203      ABI_ALLOCATE(rgramxbx,(nvec,nvec))
5204      ABI_ALLOCATE(rblockvectorx,(vectsize,nvec))
5205      ABI_ALLOCATE(rblockvectorbx,(vectsize,nvec))
5206      do ivec=1,nvec
5207        if (me_g0 == 1) then
5208          call abi_xcopy(1,vecnm(1,cgindex(ivec)),1,rblockvectorx (1,ivec),1)
5209          call abi_xcopy(rvectsiz-1,vecnm(1,cgindex(ivec)+1),2,rblockvectorx(2,ivec),1)
5210          call abi_xcopy(rvectsiz-1,vecnm(2,cgindex(ivec)+1),2,rblockvectorx(rvectsiz+1,ivec),1)
5211          if (useoverlap == 1) then
5212            call abi_xcopy(1,ovl_vecnm(1,gscindex(ivec)),1,rblockvectorbx(1,ivec),1)
5213            call abi_xcopy(rvectsiz-1,ovl_vecnm(1,gscindex(ivec)+1),2,rblockvectorbx(2,ivec),1)
5214            call abi_xcopy(rvectsiz-1,ovl_vecnm(2,gscindex(ivec)+1),2,rblockvectorbx(rvectsiz+1,ivec),1)
5215          else
5216            call abi_xcopy(1,vecnm(1,cgindex(ivec)),1,rblockvectorbx(1,ivec),1)
5217            call abi_xcopy(rvectsiz-1,vecnm(1,cgindex(ivec)+1),2,rblockvectorbx(2,ivec),1)
5218            call abi_xcopy(rvectsiz-1,vecnm(2,cgindex(ivec)+1),2,rblockvectorbx(rvectsiz+1,ivec),1)
5219          end if
5220          rblockvectorx (2:vectsize,ivec)=rblockvectorx (2:vectsize,ivec)*sqrt2
5221          rblockvectorbx(2:vectsize,ivec)=rblockvectorbx(2:vectsize,ivec)*sqrt2
5222        else
5223          call abi_xcopy(rvectsiz,vecnm(1,cgindex(ivec)),2,rblockvectorx(1,ivec),1)
5224          call abi_xcopy(rvectsiz,vecnm(2,cgindex(ivec)),2,rblockvectorx(rvectsiz+1,ivec),1)
5225          if (useoverlap == 1) then
5226            call abi_xcopy(rvectsiz,ovl_vecnm(1,gscindex(ivec)),2,rblockvectorbx(1,ivec),1)
5227            call abi_xcopy(rvectsiz,ovl_vecnm(2,gscindex(ivec)),2,rblockvectorbx(rvectsiz+1,ivec),1)
5228          else
5229            call abi_xcopy(rvectsiz,vecnm(1,cgindex(ivec)),2,rblockvectorbx(1,ivec),1)
5230            call abi_xcopy(rvectsiz,vecnm(2,cgindex(ivec)),2,rblockvectorbx(rvectsiz+1,ivec),1)
5231          end if
5232          rblockvectorx (1:vectsize,ivec)=rblockvectorx (1:vectsize,ivec)*sqrt2
5233          rblockvectorbx(1:vectsize,ivec)=rblockvectorbx(1:vectsize,ivec)*sqrt2
5234        end if
5235      end do
5236 
5237      call ortho_reim(rblockvectorx,rblockvectorbx,nvec,comm,rgramxbx,vectsize)
5238 
5239      do ivec=1,nvec
5240        ! Unpack results
5241        if (me_g0 == 1) then
5242          call abi_xcopy(1,rblockvectorx(1,ivec),1,vecnm(1,cgindex(ivec)),1)
5243          vecnm(2,cgindex(ivec))=zero
5244          rblockvectorx(2:vectsize,ivec)=rblockvectorx(2:vectsize,ivec)/sqrt2
5245          call abi_xcopy(rvectsiz-1,rblockvectorx(2,ivec),1,vecnm(1,cgindex(ivec)+1),2)
5246          call abi_xcopy(rvectsiz-1,rblockvectorx(rvectsiz+1,ivec),1,vecnm(2,cgindex(ivec)+1),2)
5247        else
5248          rblockvectorx(1:vectsize,ivec)=rblockvectorx(1:vectsize,ivec)/sqrt2
5249          call abi_xcopy(rvectsiz,rblockvectorx(1,ivec),1,vecnm(1,cgindex(ivec)),2)
5250          call abi_xcopy(rvectsiz,rblockvectorx(rvectsiz+1,ivec),1,vecnm(2,cgindex(ivec)),2)
5251        end if
5252 
5253        if(useoverlap == 1) then
5254          call abi_xtrsm('r','u','n','n',vectsize,nvec,one,rgramxbx,nvec,rblockvectorbx,vectsize)
5255          if (me_g0 == 1) then
5256            call abi_xcopy(1,rblockvectorbx(1,ivec),1,ovl_vecnm(1,gscindex(ivec)),1)
5257            ovl_vecnm(2,gscindex(ivec))=zero
5258            rblockvectorbx(2:vectsize,ivec)=rblockvectorbx(2:vectsize,ivec)/sqrt2
5259            call abi_xcopy(rvectsiz-1,rblockvectorbx(2,ivec),1,ovl_vecnm(1,gscindex(ivec)+1),2)
5260            call abi_xcopy(rvectsiz-1,rblockvectorbx(rvectsiz+1,ivec),1,ovl_vecnm(2,gscindex(ivec)+1),2)
5261          else
5262            rblockvectorbx(1:vectsize,ivec)=rblockvectorbx(1:vectsize,ivec)/sqrt2
5263            call abi_xcopy(rvectsiz,rblockvectorbx(1,ivec),1,ovl_vecnm(1,gscindex(ivec)),2)
5264            call abi_xcopy(rvectsiz,rblockvectorbx(rvectsiz+1,ivec),1,ovl_vecnm(2,gscindex(ivec)),2)
5265          end if
5266        end if
5267      end do
5268      ABI_DEALLOCATE(rgramxbx)
5269      ABI_DEALLOCATE(rblockvectorx)
5270      ABI_DEALLOCATE(rblockvectorbx)
5271    end if
5272 
5273  else if (ortalgo==4) then
5274 !  else if (ANY(ortalgo==(/0,2/))) then
5275 
5276    cg_idx = cgindex(1)
5277    if (useoverlap==0) then
5278      call cgnc_gramschmidt(nelem,nvec,vecnm(1,cg_idx),istwf_k,me_g0,comm)
5279    else
5280      gsc_idx = gscindex(1)
5281      call cgpaw_gramschmidt(nelem,nvec,vecnm(1,cg_idx),ovl_vecnm(1,gsc_idx),istwf_k,me_g0,comm)
5282    end if
5283 
5284  else if (ANY(ortalgo==(/0,2/))) then
5285 !  =======================
5286 !  Third (old) algorithm
5287 !  =======================
5288 
5289    do ivec=1,nvec
5290 !    Normalize each vecnm(n,m) in turn:
5291 
5292      if (useoverlap==1) then ! Using overlap S...
5293        if(istwf_k/=2)then
5294          sum=zero;ii0=1
5295        else
5296          if (me_g0 ==1) then
5297            sum=half*ovl_vecnm(1,1+nelem*(ivec-1)+igsc)*vecnm(1,1+nelem*(ivec-1)+icg)
5298            ii0=2
5299          else
5300            sum=zero;ii0=1
5301          end if
5302        end if
5303 !$OMP PARALLEL DO PRIVATE(ii) REDUCTION(+:sum) SHARED(icg,ivec,nelem,vecnm)
5304        do ii=ii0+nelem*(ivec-1),nelem*ivec
5305          sum=sum+vecnm(1,ii+icg)*ovl_vecnm(1,ii+igsc)+vecnm(2,ii+icg)*ovl_vecnm(2,ii+igsc)
5306        end do
5307 
5308      else ! Without overlap...
5309        if(istwf_k/=2)then
5310          sum=zero;ii0=1
5311        else
5312          if (me_g0 ==1) then
5313            sum=half*vecnm(1,1+nelem*(ivec-1)+icg)**2
5314            ii0=2
5315          else
5316            sum=zero;ii0=1
5317          end if
5318        end if
5319 !$OMP PARALLEL DO PRIVATE(ii) REDUCTION(+:sum) SHARED(icg,ivec,nelem,vecnm)
5320        do ii=ii0+nelem*(ivec-1)+icg,nelem*ivec+icg
5321          sum=sum+vecnm(1,ii)**2+vecnm(2,ii)**2
5322        end do
5323      end if
5324 
5325      call timab(48,1,tsec)
5326      call xmpi_sum(sum,comm,ierr)
5327      call timab(48,2,tsec)
5328 
5329      if(istwf_k>=2)sum=two*sum
5330      xnorm = sqrt(abs(sum)) ;  sum=1.0_dp/xnorm
5331 !$OMP PARALLEL DO PRIVATE(ii) SHARED(icg,ivec,nelem,sum,vecnm)
5332      do ii=1+nelem*(ivec-1)+icg,nelem*ivec+icg
5333        vecnm(1,ii)=vecnm(1,ii)*sum
5334        vecnm(2,ii)=vecnm(2,ii)*sum
5335      end do
5336      if (useoverlap==1) then
5337 !$OMP PARALLEL DO PRIVATE(ii) SHARED(icg,ivec,nelem,sum,ovl_vecnm)
5338        do ii=1+nelem*(ivec-1)+igsc,nelem*ivec+igsc
5339          ovl_vecnm(1,ii)=ovl_vecnm(1,ii)*sum
5340          ovl_vecnm(2,ii)=ovl_vecnm(2,ii)*sum
5341        end do
5342      end if
5343 
5344 !    Remove projection in all higher states.
5345      if (ivec<nvec) then
5346 
5347        if(istwf_k==1)then
5348 !        Cannot use time-reversal symmetry
5349 
5350          if (useoverlap==1) then ! Using overlap.
5351            do ivec2=ivec+1,nvec
5352 !            First compute scalar product
5353              dotr=zero ; doti=zero
5354              ii1=nelem*(ivec-1)+icg;ii2=nelem*(ivec2-1)+igsc
5355 !$OMP PARALLEL DO PRIVATE(ii) REDUCTION(+:doti,dotr) SHARED(ii1,ii2,nelem,vecnm)
5356              do ii=1,nelem
5357                dotr=dotr+vecnm(1,ii1+ii)*ovl_vecnm(1,ii2+ii)+vecnm(2,ii1+ii)*ovl_vecnm(2,ii2+ii)
5358                doti=doti+vecnm(1,ii1+ii)*ovl_vecnm(2,ii2+ii)-vecnm(2,ii1+ii)*ovl_vecnm(1,ii2+ii)
5359              end do
5360 
5361              call timab(48,1,tsec)
5362              buffer2(1)=doti;buffer2(2)=dotr
5363              call xmpi_sum(buffer2,comm,ierr)
5364              call timab(48,2,tsec)
5365              doti=buffer2(1)
5366              dotr=buffer2(2)
5367 
5368 !            Then subtract the appropriate amount of the lower state
5369              ii1=nelem*(ivec-1)+icg;ii2=nelem*(ivec2-1)+icg
5370 #ifdef FC_INTEL
5371 !            DIR$ ivdep
5372 #endif
5373 !$OMP PARALLEL DO PRIVATE(ii) SHARED(doti,dotr,ii1,ii2,nelem,vecnm)
5374              do ii=1,nelem
5375                vecnm(1,ii2+ii)=vecnm(1,ii2+ii)-dotr*vecnm(1,ii1+ii)+doti*vecnm(2,ii1+ii)
5376                vecnm(2,ii2+ii)=vecnm(2,ii2+ii)-doti*vecnm(1,ii1+ii)-dotr*vecnm(2,ii1+ii)
5377              end do
5378 
5379              ii1=nelem*(ivec-1)+igsc;ii2=nelem*(ivec2-1)+igsc
5380              do ii=1,nelem
5381                ovl_vecnm(1,ii2+ii)=ovl_vecnm(1,ii2+ii)&
5382 &               -dotr*ovl_vecnm(1,ii1+ii)&
5383 &               +doti*ovl_vecnm(2,ii1+ii)
5384                ovl_vecnm(2,ii2+ii)=ovl_vecnm(2,ii2+ii)&
5385                -doti*ovl_vecnm(1,ii1+ii)&
5386 &               -dotr*ovl_vecnm(2,ii1+ii)
5387              end do
5388            end do
5389          else
5390 !          ----- No overlap -----
5391            do ivec2=ivec+1,nvec
5392 !            First compute scalar product
5393              dotr=zero ; doti=zero
5394              ii1=nelem*(ivec-1)+icg;ii2=nelem*(ivec2-1)+icg
5395 !$OMP PARALLEL DO PRIVATE(ii) REDUCTION(+:doti,dotr) SHARED(ii1,ii2,nelem,vecnm)
5396              do ii=1,nelem
5397                dotr=dotr+vecnm(1,ii1+ii)*vecnm(1,ii2+ii)+&
5398 &               vecnm(2,ii1+ii)*vecnm(2,ii2+ii)
5399                doti=doti+vecnm(1,ii1+ii)*vecnm(2,ii2+ii)-&
5400 &               vecnm(2,ii1+ii)*vecnm(1,ii2+ii)
5401              end do
5402 !            Init mpi_comm
5403              buffer2(1)=doti
5404              buffer2(2)=dotr
5405              call timab(48,1,tsec)
5406              call xmpi_sum(buffer2,comm,ierr)
5407 !            call xmpi_sum(doti,spaceComm,ierr)
5408 !            call xmpi_sum(dotr,spaceComm,ierr)
5409              call timab(48,2,tsec)
5410              doti=buffer2(1)
5411              dotr=buffer2(2)
5412 
5413 !            Then subtract the appropriate amount of the lower state
5414 #ifdef FC_INTEL
5415 !            DIR$ ivdep
5416 #endif
5417 !$OMP PARALLEL DO PRIVATE(ii) SHARED(doti,dotr,ii1,ii2,nelem,vecnm)
5418              do ii=1,nelem
5419                vecnm(1,ii2+ii)=vecnm(1,ii2+ii)-dotr*vecnm(1,ii1+ii)+&
5420 &               doti*vecnm(2,ii1+ii)
5421                vecnm(2,ii2+ii)=vecnm(2,ii2+ii)-doti*vecnm(1,ii1+ii)-&
5422 &               dotr*vecnm(2,ii1+ii)
5423              end do
5424            end do
5425 
5426          end if  ! Test on useoverlap
5427 
5428        else if(istwf_k==2)then
5429 !        At gamma point use of time-reversal symmetry saves cpu time.
5430 
5431          if (useoverlap==1) then
5432 !          ----- Using overlap -----
5433            do ivec2=ivec+1,nvec
5434 !            First compute scalar product
5435              ii1=nelem*(ivec-1)+icg;ii2=nelem*(ivec2-1)+igsc
5436              if (me_g0 ==1) then
5437                dotr=half*vecnm(1,ii1+1)*ovl_vecnm(1,ii2+1)
5438 !              Avoid double counting G=0 contribution
5439 !              Imaginary part of vecnm at G=0 should be zero,so only take real part
5440 !$OMP PARALLEL DO PRIVATE(ii) REDUCTION(+:dotr) SHARED(ii1,ii2,nelem,vecnm)
5441                do ii=2,nelem
5442                  dotr=dotr+vecnm(1,ii1+ii)*ovl_vecnm(1,ii2+ii)+&
5443 &                 vecnm(2,ii1+ii)*ovl_vecnm(2,ii2+ii)
5444                end do
5445              else
5446                dotr=0._dp
5447 !$OMP PARALLEL DO PRIVATE(ii) REDUCTION(+:dotr) SHARED(ii1,ii2,nelem,vecnm)
5448                do ii=1,nelem
5449                  dotr=dotr+vecnm(1,ii1+ii)*ovl_vecnm(1,ii2+ii)+&
5450 &                 vecnm(2,ii1+ii)*ovl_vecnm(2,ii2+ii)
5451                end do
5452              end if
5453 
5454              dotr=two*dotr
5455 
5456              call timab(48,1,tsec)
5457              call xmpi_sum(dotr,comm,ierr)
5458              call timab(48,2,tsec)
5459 
5460 !            Then subtract the appropriate amount of the lower state
5461              ii1=nelem*(ivec-1)+icg;ii2=nelem*(ivec2-1)+icg
5462 #ifdef FC_INTEL
5463 !            DIR$ ivdep
5464 #endif
5465 !$OMP PARALLEL DO PRIVATE(ii) SHARED(dotr,ii1,ii2,nelem,vecnm)
5466              do ii=1,nelem
5467                vecnm(1,ii2+ii)=vecnm(1,ii2+ii)-dotr*vecnm(1,ii1+ii)
5468                vecnm(2,ii2+ii)=vecnm(2,ii2+ii)-dotr*vecnm(2,ii1+ii)
5469              end do
5470              ii1=nelem*(ivec-1)+igsc;ii2=nelem*(ivec2-1)+igsc
5471              do ii=1,nelem
5472                ovl_vecnm(1,ii2+ii)=ovl_vecnm(1,ii2+ii)-dotr*ovl_vecnm(1,ii1+ii)
5473                ovl_vecnm(2,ii2+ii)=ovl_vecnm(2,ii2+ii)-dotr*ovl_vecnm(2,ii1+ii)
5474              end do
5475            end do
5476          else
5477 !          ----- No overlap -----
5478            do ivec2=ivec+1,nvec
5479 !            First compute scalar product
5480              ii1=nelem*(ivec-1)+icg;ii2=nelem*(ivec2-1)+icg
5481              if (me_g0 ==1) then
5482 !              Avoid double counting G=0 contribution
5483 !              Imaginary part of vecnm at G=0 should be zero,so only take real part
5484                dotr=half*vecnm(1,ii1+1)*vecnm(1,ii2+1)
5485 !$OMP PARALLEL DO PRIVATE(ii) REDUCTION(+:dotr) SHARED(ii1,ii2,nelem,vecnm)
5486                do ii=2,nelem
5487                  dotr=dotr+vecnm(1,ii1+ii)*vecnm(1,ii2+ii)+vecnm(2,ii1+ii)*vecnm(2,ii2+ii)
5488                end do
5489              else
5490                dotr=0._dp
5491 !$OMP PARALLEL DO PRIVATE(ii) REDUCTION(+:dotr) SHARED(ii1,ii2,nelem,vecnm)
5492                do ii=1,nelem
5493                  dotr=dotr+vecnm(1,ii1+ii)*vecnm(1,ii2+ii)+vecnm(2,ii1+ii)*vecnm(2,ii2+ii)
5494                end do
5495              end if
5496              dotr=two*dotr
5497 
5498              call timab(48,1,tsec)
5499              call xmpi_sum(dotr,comm,ierr)
5500              call timab(48,2,tsec)
5501 
5502 !            Then subtract the appropriate amount of the lower state
5503 #ifdef FC_INTEL
5504 !            DIR$ ivdep
5505 #endif
5506 !$OMP PARALLEL DO PRIVATE(ii) SHARED(dotr,ii1,ii2,nelem,vecnm)
5507              do ii=1,nelem
5508                vecnm(1,ii2+ii)=vecnm(1,ii2+ii)-dotr*vecnm(1,ii1+ii)
5509                vecnm(2,ii2+ii)=vecnm(2,ii2+ii)-dotr*vecnm(2,ii1+ii)
5510              end do
5511            end do
5512          end if  ! Test on useoverlap
5513 
5514        else
5515 !        At other special points,use of time-reversal symmetry saves cpu time.
5516 
5517          if (useoverlap==1) then
5518 !          ----- Using overlap -----
5519            do ivec2=ivec+1,nvec
5520 !            First compute scalar product
5521              ii1=nelem*(ivec-1)+icg;ii2=nelem*(ivec2-1)+igsc
5522 !            Avoid double counting G=0 contribution
5523 !            Imaginary part of vecnm at G=0 should be zero,so only take real part
5524              dotr=zero
5525 !$OMP PARALLEL DO PRIVATE(ii) REDUCTION(+:dotr) SHARED(ii1,ii2,nelem,vecnm)
5526              do ii=1,nelem
5527                dotr=dotr+vecnm(1,ii1+ii)*ovl_vecnm(1,ii2+ii)+vecnm(2,ii1+ii)*ovl_vecnm(2,ii2+ii)
5528              end do
5529              dotr=two*dotr
5530 
5531              call timab(48,1,tsec)
5532              call xmpi_sum(dotr,comm,ierr)
5533              call timab(48,2,tsec)
5534 
5535 !            Then subtract the appropriate amount of the lower state
5536              ii1=nelem*(ivec-1)+icg;ii2=nelem*(ivec2-1)+icg
5537 #ifdef FC_INTEL
5538 !            DIR$ ivdep
5539 #endif
5540 !$OMP PARALLEL DO PRIVATE(ii) SHARED(dotr,ii1,ii2,nelem,vecnm)
5541              do ii=1,nelem
5542                vecnm(1,ii2+ii)=vecnm(1,ii2+ii)-dotr*vecnm(1,ii1+ii)
5543                vecnm(2,ii2+ii)=vecnm(2,ii2+ii)-dotr*vecnm(2,ii1+ii)
5544              end do
5545              ii1=nelem*(ivec-1)+igsc;ii2=nelem*(ivec2-1)+igsc
5546              do ii=1,nelem
5547                ovl_vecnm(1,ii2+ii)=ovl_vecnm(1,ii2+ii)-dotr*ovl_vecnm(1,ii1+ii)
5548                ovl_vecnm(2,ii2+ii)=ovl_vecnm(2,ii2+ii)-dotr*ovl_vecnm(2,ii1+ii)
5549              end do
5550            end do
5551          else
5552 !          ----- No overlap -----
5553            do ivec2=ivec+1,nvec
5554 !            First compute scalar product
5555              ii1=nelem*(ivec-1)+icg;ii2=nelem*(ivec2-1)+icg
5556 !            Avoid double counting G=0 contribution
5557 !            Imaginary part of vecnm at G=0 should be zero,so only take real part
5558              dotr=zero
5559 !$OMP PARALLEL DO PRIVATE(ii) REDUCTION(+:dotr) SHARED(ii1,ii2,nelem,vecnm)
5560              do ii=1,nelem
5561                dotr=dotr+vecnm(1,ii1+ii)*vecnm(1,ii2+ii)+vecnm(2,ii1+ii)*vecnm(2,ii2+ii)
5562              end do
5563              dotr=two*dotr
5564 
5565              call timab(48,1,tsec)
5566              call xmpi_sum(dotr,comm,ierr)
5567              call timab(48,2,tsec)
5568 
5569 !            Then subtract the appropriate amount of the lower state
5570 !$OMP PARALLEL DO PRIVATE(ii) SHARED(dotr,ii1,ii2,nelem,vecnm)
5571              do ii=1,nelem
5572                vecnm(1,ii2+ii)=vecnm(1,ii2+ii)-dotr*vecnm(1,ii1+ii)
5573                vecnm(2,ii2+ii)=vecnm(2,ii2+ii)-dotr*vecnm(2,ii1+ii)
5574              end do
5575            end do
5576          end if
5577 
5578 !        End use of time-reversal symmetry
5579        end if
5580 
5581      end if  ! Test on "ivec"
5582 
5583 !    end loop over vectors (or bands) with index ivec :
5584    end do
5585 
5586  else
5587    write(msg,'(a,i0)')"Wrong value for ortalgo: ",ortalgo
5588    MSG_ERROR(msg)
5589  end if ! End of the second algorithm
5590 
5591 end subroutine pw_orthon

m_cgtools/set_istwfk [ Functions ]

[ Top ] [ m_cgtools ] [ Functions ]

NAME

  set_istwfk

FUNCTION

  Returns the value of istwfk associated to the input k-point.

INPUTS

  kpoint(3)=The k-point in reduced coordinates.

OUTPUT

  istwfk= Integer flag internally used in the code to define the storage mode of the wavefunctions.
  It also define the algorithm used to apply an operator in reciprocal space as well as the FFT
  algorithm used to go from G- to r-space and vice versa.

   1 => time-reversal cannot be used
   2 => use time-reversal at the Gamma point.
   3 => use time-reversal symmetry for k=(1/2, 0 , 0 )
   4 => use time-reversal symmetry for k=( 0 , 0 ,1/2)
   5 => use time-reversal symmetry for k=(1/2, 0 ,1/2)
   6 => use time-reversal symmetry for k=( 0 ,1/2, 0 )
   7 => use time-reversal symmetry for k=(1/2,1/2, 0 )
   8 => use time-reversal symmetry for k=( 0 ,1/2,1/2)
   9 => use time-reversal symmetry for k=(1/2,1/2,1/2)

  Useful relations:
   u_k(G) = u_{k+G0}(G-G0); u_{-k}(G) = u_k(G)^*
  and therefore:
   u_{G0/2}(G) = u_{G0/2}(-G-G0)^*.

PARENTS

CHILDREN

SOURCE

1149 integer pure function set_istwfk(kpoint) result(istwfk)
1150 
1151 
1152 !This section has been created automatically by the script Abilint (TD).
1153 !Do not modify the following lines by hand.
1154 #undef ABI_FUNC
1155 #define ABI_FUNC 'set_istwfk'
1156 !End of the abilint section
1157 
1158  implicit none
1159 
1160 !Arguments ------------------------------------
1161  real(dp),intent(in) :: kpoint(3)
1162 
1163 !Local variables-------------------------------
1164 !scalars
1165  integer :: bit0,ii
1166 !arrays
1167  integer :: bit(3)
1168 
1169 ! *************************************************************************
1170 
1171  bit0=1
1172 
1173  do ii=1,3
1174    if (DABS(kpoint(ii))<tol10) then
1175      bit(ii)=0
1176    else if (DABS(kpoint(ii)-half)<tol10 ) then
1177      bit(ii)=1
1178    else
1179      bit0=0
1180    end if
1181  end do
1182 
1183  if (bit0==0) then
1184    istwfk=1
1185  else
1186    istwfk=2+bit(1)+4*bit(2)+2*bit(3) ! Note the inversion between bit(2) and bit(3)
1187  end if
1188 
1189 end function set_istwfk

m_cgtools/sqnorm_g [ Functions ]

[ Top ] [ m_cgtools ] [ Functions ]

NAME

 sqnorm_g

FUNCTION

 Compute the square of the norm of one complex vector vecti, in reciprocal space
 Take into account the storage mode of the vector (istwf_k)

INPUTS

  istwf_k=option parameter that describes the storage of wfs
  npwsp= (effective) number of planewaves at this k point.
  vect(2,npwsp)=the vector in reciprocal space (npw*nspinor, usually)
  me_g0=1 if this processors treats G=0, 0 otherwise.
  comm=MPI communicator for MPI sum.

OUTPUT

  dotr= <vect|vect>

PARENTS

      cgwf,dfpt_cgwf,dfpt_vtowfk,m_epjdos,mkresi,rf2_init

CHILDREN

SOURCE

1217 subroutine sqnorm_g(dotr,istwf_k,npwsp,vect,me_g0,comm)
1218 
1219 
1220 !This section has been created automatically by the script Abilint (TD).
1221 !Do not modify the following lines by hand.
1222 #undef ABI_FUNC
1223 #define ABI_FUNC 'sqnorm_g'
1224 !End of the abilint section
1225 
1226  implicit none
1227 
1228 !Arguments ------------------------------------
1229 !scalars
1230  integer,intent(in) :: istwf_k,npwsp,me_g0,comm
1231  real(dp),intent(out) :: dotr
1232 !arrays
1233  real(dp),intent(in) :: vect(2,npwsp)
1234 
1235 !Local variables-------------------------------
1236 !scalars
1237  integer :: ierr
1238 
1239 ! *************************************************************************
1240 
1241  if (istwf_k==1) then ! General k-point
1242    !dotr = cg_real_zdotc(npwsp,vect,vect)
1243    dotr = cg_dznrm2(npwsp, vect)
1244    dotr = dotr * dotr
1245 
1246  else
1247    if (istwf_k==2 .and. me_g0==1) then
1248      ! Gamma k-point and I have G=0
1249      dotr=half*vect(1,1)**2
1250      dotr = dotr + cg_real_zdotc(npwsp-1,vect(1,2),vect(1,2))
1251    else
1252      ! Other TR k-points
1253      dotr = cg_real_zdotc(npwsp,vect,vect)
1254    end if
1255    dotr=two*dotr
1256  end if
1257 
1258  if (xmpi_comm_size(comm)>1) then
1259    call xmpi_sum(dotr,comm,ierr)
1260  end if
1261 
1262 end subroutine sqnorm_g

m_cgtools/sqnorm_v [ Functions ]

[ Top ] [ m_cgtools ] [ Functions ]

NAME

 sqnorm_v

FUNCTION

 Compute square of the norm of a potential (integral over FFT grid), to obtain
 a square residual-like quantity (so the sum of product of values
 is NOT divided by the number of FFT points, and NOT multiplied by the primitive cell volume).
 Take into account the spin components of the potentials (nspden),
 and sum over them.

INPUTS

  cplex=if 1, real space function on FFT grid is REAL, if 2, COMPLEX
  nfft= (effective) number of FFT grid points (for this processor)
  nspden=number of spin-density components
  opt_storage: 0, if potential is stored as V^up-up, V^dn-dn, Re[V^up-dn], Im[V^up-dn]
               1, if potential is stored as V, B_x, B_y, Bz  (B=magn. field)
  pot(cplex*nfft,nspden)=real space potential on FFT grid

OUTPUT

  norm2= value of the square of the norm

PARENTS

      dfpt_rhotov,dfpt_vtorho,rhotov,vtorho

CHILDREN

SOURCE

1936 subroutine sqnorm_v(cplex,nfft,norm2,nspden,opt_storage,pot,mpi_comm_sphgrid)
1937 
1938 
1939 !This section has been created automatically by the script Abilint (TD).
1940 !Do not modify the following lines by hand.
1941 #undef ABI_FUNC
1942 #define ABI_FUNC 'sqnorm_v'
1943 !End of the abilint section
1944 
1945  implicit none
1946 
1947 !Arguments ------------------------------------
1948 !scalars
1949  integer,intent(in) :: cplex,nfft,nspden,opt_storage
1950  integer,intent(in),optional :: mpi_comm_sphgrid
1951  real(dp),intent(out) :: norm2
1952 !arrays
1953  real(dp),intent(in) :: pot(cplex*nfft,nspden)
1954 
1955 !Local variables-------------------------------
1956 !scalars
1957  integer :: ierr,ifft,ispden,nproc_sphgrid
1958  real(dp) :: ar
1959 
1960 ! *************************************************************************
1961 
1962 !Real or complex inputs are coded
1963 
1964  norm2=zero
1965  do ispden=1,min(nspden,2)
1966 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,ispden,nfft,pot) REDUCTION(+:norm2)
1967    do ifft=1,cplex*nfft
1968      norm2=norm2 + pot(ifft,ispden)**2
1969    end do
1970  end do
1971  if (nspden==4) then
1972    ar=zero
1973    do ispden=3,4
1974 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(cplex,ispden,nfft,pot) REDUCTION(+:ar)
1975      do ifft=1,cplex*nfft
1976        ar=ar + pot(ifft,ispden)**2
1977      end do
1978    end do
1979    if (opt_storage==0) then
1980      if (cplex==1) then
1981        norm2=norm2+two*ar
1982      else
1983        norm2=norm2+ar
1984      end if
1985    else
1986      norm2=half*(norm2+ar)
1987    end if
1988  end if
1989 
1990 !MPIWF reduction (addition) on norm2 is needed here
1991  if(present(mpi_comm_sphgrid)) then
1992    nproc_sphgrid=xmpi_comm_size(mpi_comm_sphgrid)
1993    if(nproc_sphgrid>1)then
1994      call xmpi_sum(norm2,mpi_comm_sphgrid,ierr)
1995    end if
1996  end if
1997 
1998 end subroutine sqnorm_v