TABLE OF CONTENTS


ABINIT/m_gsphere [ Modules ]

[ Top ] [ Modules ]

NAME

  m_gsphere

FUNCTION

   The Gsphere data type defines the set of G-vectors
   centered on Gamma used to describe (chi0|epsilon|W) in the GW code.
   Note that, unlike the kg_k arrays used for wavefunctions, here the
   G-vectors are ordered in shells (increasing length). Moreover
   the sphere can be enlarged to take into account umklapps for which
   one need the knowledge of several quantities at G-G0.

COPYRIGHT

 Copyright (C) 1999-2018 ABINIT group (MG, GMR, VO, LR, RWG, MT, XG)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .

PARENTS

SOURCE

24 #if defined HAVE_CONFIG_H
25 #include "config.h"
26 #endif
27 
28 #include "abi_common.h"
29 
30 MODULE m_gsphere
31 
32  use defs_basis
33  use m_abicore
34  use m_errors
35  use m_distribfft
36  use m_sort
37 
38  use defs_abitypes,   only : MPI_type
39  use m_fstrings,      only : sjoin, itoa
40  use m_numeric_tools, only : bisect
41  use m_geometry,      only : normv
42  use m_crystal,       only : crystal_t
43  use m_fftcore,       only : kpgsph, kgindex, sphereboundary
44  use m_mpinfo,        only : destroy_mpi_enreg, initmpi_seq
45 
46  implicit none
47 
48  private
49 
50 ! Low-level procedures.
51  public :: merge_and_sort_kg   ! Merges a set of k-centered G-spheres of cutoff ecut. Return a Gamma-centered G-spheres.
52  public :: table_gbig2kg       ! Associate the kg_k set of G-vectors with Gamma-centered G-sphere.
53  public :: get_irredg          ! Given a set of G vectors, find the set of G"s generating the others by symmetry.
54  public :: merge_kgirr         ! Merge a list of irreducible G vectors (see routine for more info)
55  public :: setshells           ! Set consistently the number of shells, the number of plane-waves,  and the energy cut-off
56  public :: kg_map              ! Compute the mapping between two lists of g-vectors.
57  public :: make_istwfk_table
58  public :: getkpgnorm          ! Compute the norms of the k+G vectors
59  public :: symg

m_gpshere/setshells [ Functions ]

[ Top ] [ Functions ]

NAME

 setshells

FUNCTION

 Set consistently the number of shells, the number of plane-waves, and the energy cut-off

INPUTS

  nsym=number of symmetry operations
  gmet(3,3)=metric tensor in reciprocal space
  gprimd(3,3)=dimensional primitive vectors in reciprocal space
  symrel(3,3,nsym)=symmetry operations in real space
  tag=suffix to account for the different possibilities for these variables (npw, ecut or nsh ..)
  ucvol=unit cell volume

OUTPUT

  (see side effects)

SIDE EFFECTS

  ecut,npw,nsh=one of them is an input, the two others are output
  ecut=cut-off energy for plane wave basis sphere (Ha)
  npw=number of plane waves
  nsh=number of shells

PARENTS

      invars2m,setup_screening,setup_sigma

CHILDREN

      gsph_free,gsph_init

SOURCE

1853 subroutine setshells(ecut,npw,nsh,nsym,gmet,gprimd,symrel,tag,ucvol)
1854 
1855 
1856 !This section has been created automatically by the script Abilint (TD).
1857 !Do not modify the following lines by hand.
1858 #undef ABI_FUNC
1859 #define ABI_FUNC 'setshells'
1860 !End of the abilint section
1861 
1862  implicit none
1863 
1864 !Arguments ------------------------------------
1865 !scalars
1866  integer,intent(in) :: nsym
1867  integer,intent(inout) :: npw,nsh
1868  real(dp),intent(in) :: ucvol
1869  real(dp),intent(inout) :: ecut
1870  character(len=*),intent(in) :: tag
1871 !arrays
1872  integer,intent(in) :: symrel(3,3,nsym)
1873  real(dp),intent(in) :: gmet(3,3),gprimd(3,3)
1874 
1875 !Local variables-------------------------------
1876 !scalars
1877  integer :: exchn2n3d,ifound,ig,ii,ish,isym,npw_found,npwave
1878  integer :: npwwrk,nsh_found,pad=50
1879  real(dp) :: ecut_found,ecut_trial,eps,scale=1.3_dp
1880  logical :: found
1881  character(len=500) :: msg
1882  type(MPI_type) :: MPI_enreg_seq
1883 !arrays
1884  integer :: geq(3)
1885  integer,allocatable :: gvec(:,:),gvec_sh(:,:),insort(:),npw_sh(:)
1886  real(dp) :: gctr(3)
1887  real(dp),allocatable :: gnorm(:),gnorm_sh(:)
1888 
1889 !******************************************************************
1890 
1891  DBG_ENTER("COLL")
1892 !
1893 !=== Check coherence of input variables ecut, npw, and nsh ===
1894 !1-> one at least should be non-null
1895  if (npw==0.and.nsh==0.and.ecut<=tol6) then
1896    write(msg,'(8a)')&
1897 &   'One of the three variables ecut',TRIM(tag),', npw',TRIM(tag),', or nsh',TRIM(tag),ch10,&
1898 &   'must be non-null. Returning.'
1899    MSG_COMMENT(msg)
1900    RETURN
1901  end if
1902 !2-> one and only one should be non-null
1903  if (npw/=0.and.nsh/=0) then
1904    write(msg,'(6a)')&
1905 &   'Only one of the two variables npw',TRIM(tag),' and nsh',TRIM(tag),ch10,&
1906 &   'can be non-null. Modify the value of one of these in input file.'
1907    MSG_ERROR(msg)
1908  end if
1909  if (ecut>tol6.and.npw/=0) then
1910    write(msg,'(6a)')&
1911 &   'Only one of the two variables ecut',TRIM(tag),' and npw',TRIM(tag),ch10,&
1912 &   'can be non-null. Modify the value of one of these in input file.'
1913    MSG_ERROR(msg)
1914  end if
1915  if (ecut>tol6.and.nsh/=0) then
1916    write(msg,'(6a)')&
1917 &   'Only one of the two variables ecut',TRIM(tag),' and nsh',TRIM(tag),ch10,&
1918 &   'can be non-null Action : modify the value of one of these in input file.'
1919    MSG_ERROR(msg)
1920  end if
1921 !
1922 !=== Calculates an upper bound for npw ===
1923 !* gctr is center of the g-vector sphere
1924  gctr(:)= [zero,zero,zero]
1925  if (ecut>tol6) then
1926 !  The average number of plane-waves in the cutoff sphere is given by:
1927 !  npwave = (2*ecut)**(3/2)*ucvol/(6*pi**2)
1928 !  The upper bound is calculated as npwwrk=int(scale * npwave) + pad
1929    npwave=NINT(ucvol*(two*ecut)**1.5_dp/(six*pi**2))
1930    npwwrk=NINT(DBLE(npwave)*scale)+pad
1931    ecut_trial=ecut
1932  else if (npw/=0) then
1933 !  npw is given in the input
1934    npwwrk=NINT(DBLE(npw)*scale)+pad
1935    ecut_trial=(six*pi**2*npw/ucvol)**two_thirds/two
1936  else
1937 !  If nsh is given in the input
1938    npwwrk=nsh*18+2*pad
1939    ecut_trial=(six*pi**2*nsh*18/ucvol)**two_thirds/two
1940  end if
1941 
1942  call initmpi_seq(MPI_enreg_seq)
1943 
1944  ABI_MALLOC(gvec,(3,npwwrk))
1945  ifound=0
1946  do while(ifound==0)
1947    !write(msg,'(a,f8.2)')' setshells : ecut_trial = ',ecut_trial
1948    !call wrtout(std_out,msg,'COLL')
1949    exchn2n3d=0 ! For the time being, no exchange of n2 and n3
1950 
1951    call kpgsph(ecut_trial,exchn2n3d,gmet,0,1,1,gvec,gctr,1,MPI_enreg_seq,npwwrk,npw_found)
1952 
1953    ABI_MALLOC(gnorm,(npw_found))
1954    ABI_MALLOC(insort,(npw_found))
1955 
1956    do ig=1,npw_found
1957      insort(ig)=ig
1958      gnorm(ig)=zero
1959      do ii=1,3
1960        gnorm(ig)=gnorm(ig)+(gvec(1,ig)*gprimd(ii,1)+&
1961 &       gvec(2,ig)*gprimd(ii,2)+&
1962 &       gvec(3,ig)*gprimd(ii,3))**2
1963      end do
1964    end do
1965    call sort_dp(npw_found,gnorm,insort,tol14)
1966    ABI_MALLOC(npw_sh,(npw_found))
1967    ABI_MALLOC(gnorm_sh,(npw_found))
1968    ABI_MALLOC(gvec_sh,(3,npw_found))
1969    npw_sh(:)=0
1970    gnorm_sh(:)=zero
1971    gvec_sh(:,:)=0
1972 !  Count the number of shells:
1973 !  (search for the G-vectors generating the others by symmetry)
1974    nsh_found=0
1975 
1976    do ig=1,npw_found
1977      eps=1.d-8*gnorm(ig)
1978      found=.FALSE.
1979      ish=1
1980      do while ((.not.found).and.(ish<=nsh_found))
1981        if (ABS(gnorm(ig)-gnorm_sh(ish))<=eps) then
1982          isym=1
1983          do while ((.not.found).and.(isym<=nsym))
1984            geq(:)=(symrel(1,:,isym)*gvec(1,insort(ig))+&
1985 &           symrel(2,:,isym)*gvec(2,insort(ig))+&
1986 &           symrel(3,:,isym)*gvec(3,insort(ig)))
1987 
1988            found=((geq(1)==gvec_sh(1,ish)).and.&
1989 &           (geq(2)==gvec_sh(2,ish)).and.&
1990 &           (geq(3)==gvec_sh(3,ish)))
1991            isym=isym+1
1992          end do
1993        end if
1994        ish=ish+1
1995      end do
1996      if (.not.found) then
1997        nsh_found=nsh_found+1
1998        gnorm_sh(nsh_found)=gnorm(ig)
1999        gvec_sh(:,nsh_found)=gvec(:,insort(ig))
2000        npw_sh(nsh_found)=1
2001      else
2002        ish=ish-1
2003        npw_sh(ish)=npw_sh(ish)+1
2004      end if
2005    end do
2006 
2007    ecut_found=two*pi**2*gnorm(npw_found)
2008 
2009    if(ecut>tol6) then
2010 !    ecut is given in the input
2011      if (ecut_found<ecut-0.1) then
2012        write(msg,'(3a,e14.6,9a,e14.6,3a)')&
2013 &       'The value ecut',TRIM(tag),'=',ecut,' given in the input file leads to',ch10,&
2014 &       'the same values for nsh',TRIM(tag),' and npw',TRIM(tag),' as ecut',TRIM(tag),'=',ecut_found,ch10,&
2015 &       'This value will be adopted for the calculation.',ch10
2016        MSG_WARNING(msg)
2017      end if
2018      ifound=1
2019    else if (npw/=0) then
2020 !    If npw is given in the input
2021      if (npw_found==npw) then
2022        ecut_found=two*pi**2*gnorm(npw_found)
2023        ifound=1
2024      else if (npw_found>npw) then
2025        npw_found=0
2026        nsh_found=0
2027        do while (npw_found<npw)
2028          nsh_found=nsh_found+1
2029          npw_found=npw_found+npw_sh(nsh_found)
2030        end do
2031 !      check that the shell is closed
2032        if(npw_found>npw) then
2033 !        shell not closed
2034          npw_found=npw_found-npw_sh(nsh_found)
2035          nsh_found=nsh_found-1
2036          do while (ABS(gnorm_sh(nsh_found)-gnorm_sh(nsh_found+1))<0.000001)
2037            npw_found=npw_found-npw_sh(nsh_found)
2038            nsh_found=nsh_found-1
2039          end do
2040          write(msg,'(3a,i6,5a,i6,3a)')&
2041 &         'The value npw',TRIM(tag),'=',npw,' given in the input file does not close the shell',ch10,&
2042 &         'The lower closed-shell is obtained for a value npw',TRIM(tag),'=',npw_found,ch10,&
2043 &         'This value will be adopted for the calculation.',ch10
2044          MSG_WARNING(msg)
2045        end if
2046        ecut_found=two*pi**2*gnorm(npw_found)
2047        ifound=1
2048      end if
2049    else if (nsh/=0) then
2050 !    If nsh is given in the input
2051      if (nsh_found==nsh) then
2052        ecut_found=two*pi**2*gnorm(npw_found)
2053        ifound=1
2054      else if (nsh_found>nsh) then
2055        npw_found=0
2056        nsh_found=0
2057        do ish=1,nsh
2058          npw_found=npw_found+npw_sh(ish)
2059          nsh_found=nsh_found+1
2060        end do
2061        if (ABS(gnorm_sh(nsh_found)-gnorm_sh(nsh_found+1))<0.000001) then
2062          do while (ABS(gnorm_sh(nsh_found)-gnorm_sh(nsh_found+1))<0.000001)
2063            nsh_found=nsh_found+1
2064            npw_found=npw_found+npw_sh(nsh_found)
2065          end do
2066          write(msg,'(3a,i6,5a,i6,3a)')&
2067 &         'The value nsh',TRIM(tag),'=',nsh,' given in the input file corresponds to the same',ch10,&
2068 &         'cut-off energy as for closed-shell upto nsh',TRIM(tag),'=',nsh_found,ch10,&
2069 &         'This value will be adopted for the calculation.',ch10
2070          MSG_WARNING(msg)
2071        end if
2072        ecut_found=two*pi**2*gnorm(npw_found)
2073        ifound=1
2074      end if
2075    end if
2076 
2077    if (ifound==0) then
2078      ecut_trial=1.1*ecut_trial
2079      ABI_FREE(gnorm)
2080      ABI_FREE(gnorm_sh)
2081      ABI_FREE(gvec_sh)
2082      ABI_FREE(insort)
2083      ABI_FREE(npw_sh)
2084    else
2085      ecut=ecut_found
2086      npw=npw_found
2087      nsh=nsh_found
2088    end if
2089 
2090  end do !while(ifound==0)
2091 
2092  call destroy_mpi_enreg(MPI_enreg_seq)
2093 
2094  ABI_FREE(gnorm)
2095  ABI_FREE(gnorm_sh)
2096  ABI_FREE(gvec)
2097  ABI_FREE(gvec_sh)
2098  ABI_FREE(insort)
2099  ABI_FREE(npw_sh)
2100 
2101  DBG_EXIT("COLL")
2102 
2103 end subroutine setshells

m_gsphere/get_irredg [ Functions ]

[ Top ] [ m_gsphere ] [ Functions ]

NAME

 get_irredg

FUNCTION

  Given a set of reciprocal lattice vectors, find the set of G"s generating the others by symmetry.

INPUTS

  nsym=number of symmetry operations
  pinv=-1 if time-reversal can be used, 1 otherwise
  npw_k=number of G vectors (for this k-point, as the set of G is k-centered)
  gcurr(3,npw_k)=the list of G vectors
  gprimd(3,3)=dimensional primitive translations for reciprocal space ($\textrm{bohr}^{-1}$)
  symrec(3,3,nsym)=symmetry operations in terms of reciprocal space primitive translations.

OUTPUT

  nbasek=number of irreducible G vectors found
  cnormk(npw_k)=first nbasek elements are the norm of each irreducible G-vector
  gbasek(3,npw_k)=first nbasek elements are the irreducible G vectors

NOTES

  The search can be optimized by looping over shells. See m_skw for a faster algo

PARENTS

      m_gsphere,m_skw

CHILDREN

      gsph_free,gsph_init

SOURCE

1634 subroutine get_irredg(npw_k,nsym,pinv,gprimd,symrec,gcurr,nbasek,gbasek,cnormk)
1635 
1636 
1637 !This section has been created automatically by the script Abilint (TD).
1638 !Do not modify the following lines by hand.
1639 #undef ABI_FUNC
1640 #define ABI_FUNC 'get_irredg'
1641 !End of the abilint section
1642 
1643  implicit none
1644 
1645 !Arguments ------------------------------------
1646 !scalars
1647  integer,intent(in) :: npw_k,nsym,pinv
1648  integer,intent(out) :: nbasek
1649 !arrays
1650  integer,intent(in) :: gcurr(3,npw_k),symrec(3,3,nsym)
1651  integer,intent(out) :: gbasek(3,npw_k)
1652  real(dp),intent(in) :: gprimd(3,3)
1653  real(dp),intent(out) :: cnormk(npw_k)
1654 
1655 !Local variables-------------------------------
1656 !scalars
1657  integer :: ig,irr,isym,jj
1658  real(dp) :: eps,norm
1659  logical :: found
1660 !arrays
1661  integer :: gbas(3),gcur(3),geq(3)
1662  real(dp) :: gcar(3)
1663 
1664 ! *************************************************************************
1665 
1666  DBG_ENTER("COLL")
1667 
1668  if (pinv/=1.and.pinv/=-1) then
1669    MSG_BUG(sjoin('pinv should be -1 or 1, however, pinv =', itoa(pinv)))
1670  end if
1671 
1672  ! Zero irred G vectors found, zeroing output arrays.
1673  nbasek = 0; cnormk(:) = zero; gbasek(:,:) = 0
1674 
1675  do ig=1,npw_k
1676    gcur(:) = gcurr(:,ig); norm = zero
1677    do jj=1,3
1678      gcar(jj)=gcur(1)*gprimd(jj,1)+gcur(2)*gprimd(jj,2)+gcur(3)*gprimd(jj,3)
1679      norm=norm+gcar(jj)**2
1680    end do
1681    eps = tol8 * norm; found = .False.; irr = 1
1682    do while (.not.found .and. irr <= nbasek)  ! This loop can be optimized by looping inside the shell.
1683      if (abs(norm - cnormk(irr)) <= eps) then
1684        gbas(:) = gbasek(:,irr); isym = 1
1685        do while (.not.found .and. isym <= nsym)
1686          geq(:) = matmul(symrec(:,:,isym),gcur)
1687          found = all(geq(:) == gbas(:))
1688          if (pinv == -1) found = (found .or. all(geq == -gbas)) ! For time-reversal
1689          isym = isym + 1
1690        end do
1691      end if
1692      irr = irr + 1
1693    end do
1694    if (.not. found) then
1695      nbasek = nbasek + 1; cnormk(nbasek) = norm; gbasek(:,nbasek) = gcur(:)
1696    end if
1697  end do
1698 
1699  DBG_EXIT("COLL")
1700 
1701 end subroutine get_irredg

m_gsphere/getfullg [ Functions ]

[ Top ] [ m_gsphere ] [ Functions ]

NAME

 getfullg

FUNCTION

  Reconstruct a G-sphere starting from a set of irreducible lattice vectors

INPUTS

  pinv=-1 if time-reversal can be used, 1 otherwise
  nsym=number of symmetry operations
  sizepw=Max expected number of G vectors in the shere
  symrec(3,3,nsym)=symmetry operation in reciprocal space
  nbase=number of irreducible G vectors
  gbase(3,nbase)=irreducible G-vectors
  cnorm(nbase)=norm of the irreducible G vectors (supposed not yet sorted)

OUTPUT

  maxpw=Number of G vectors found
  gbig(3,sizepw)=G vectors in the sphere packed in the first maxpw columns
  shlim(nbase)=number of G vectors within each shell
  ierr= Exit status, if /=0 the number of G vectors found exceeds sizepw

SIDE EFFECTS

NOTES

  cnorm is a bit redundant since it can be calculated from gbase. However this procedure
  is called by outkss in which cnorm is already calculated and we dont want to do it twice

PARENTS

      m_gsphere

CHILDREN

      gsph_free,gsph_init

SOURCE

1481 subroutine getfullg(nbase,nsym,pinv,sizepw,gbase,symrec,cnorm,maxpw,gbig,shlim,ierr)
1482 
1483 
1484 !This section has been created automatically by the script Abilint (TD).
1485 !Do not modify the following lines by hand.
1486 #undef ABI_FUNC
1487 #define ABI_FUNC 'getfullg'
1488 !End of the abilint section
1489 
1490  implicit none
1491 
1492 !Arguments ------------------------------------
1493 !scalars
1494  integer,intent(in) :: nbase,nsym,pinv,sizepw
1495  integer,intent(out) :: ierr,maxpw
1496 !arrays
1497  integer,intent(in) :: gbase(3,nbase),symrec(3,3,nsym)
1498  integer,intent(out) :: gbig(3,sizepw),shlim(nbase)
1499  real(dp),intent(inout) :: cnorm(nbase) !sort_dp can change cnorm
1500 
1501 !Local variables-------------------------------
1502 !scalars
1503  integer :: ibase,ig,ilim,ish,isym,itim
1504  logical :: found
1505  character(len=500) :: msg
1506 !arrays
1507  integer :: gcur(3),geq(3)
1508  integer,allocatable :: gshell(:,:),insort(:),nshell(:)
1509 
1510 ! *************************************************************************
1511 
1512  if (pinv/=1.and.pinv/=-1) then
1513    write(msg,'(a,i6)')&
1514 &   ' The argument pinv should be -1 or 1, however, pinv =',pinv
1515    MSG_BUG(msg)
1516  end if
1517  !
1518  ! === Reorder base g-vectors in order of increasing module ===
1519  ABI_MALLOC(insort,(nbase))
1520  do ibase=1,nbase
1521    insort(ibase)=ibase
1522  end do
1523  call sort_dp(nbase,cnorm,insort,tol14)
1524  !
1525  ! === Generate all stars of G-vectors ===
1526  ! Star of G is the set of all symetrical images of the vector
1527  ! gshell contains the symmetrical G at fixed gbase. No need to add an additional dimension
1528  ! or initialize to zero the array inside the loop over nbase as we loop over (ish<=nshell(ibase))
1529  ABI_MALLOC(nshell,(nbase))
1530  ABI_MALLOC(gshell,(3,2*nsym))
1531  !
1532  ! === Start with zero number of G vectors found ===
1533  maxpw=0 ; ierr=0
1534  do ibase=1,nbase
1535    !
1536    ! === Loop over all different modules of G ===
1537    ! * Start with zero G vectors found in this star
1538    nshell(ibase)=0
1539    gcur(:)=gbase(:,insort(ibase))
1540    !
1541    !  === Loop over symmetries ===
1542    do isym=1,nsym
1543      do itim=pinv,1,2
1544        geq(:)=itim*MATMUL(symrec(:,:,isym),gcur)
1545        !
1546        ! * Search for symetric of g and eventually add it:
1547        found=.FALSE. ; ish=1
1548        do while ((.not.found).and. (ish<=nshell(ibase)))
1549          found=ALL(geq(:)==gshell(:,ish))
1550          ish=ish+1
1551        end do
1552        if (.not.found) then
1553          nshell(ibase)=nshell(ibase)+1
1554          gshell(:,nshell(ibase))=geq(:)
1555        end if
1556      end do
1557    end do
1558    !
1559    ! * Was sizepw large enough?
1560    if ((maxpw+nshell(ibase))>sizepw) then
1561      write(msg,'(a,i6,2a)')&
1562 &     ' Number of G in sphere exceeds maximum allowed value =',sizepw,ch10,&
1563 &     ' check the value of sizepw in calling routine '
1564      MSG_WARNING(msg)
1565      ierr=1; RETURN
1566    end if
1567    !
1568    ! === Store this shell of Gs in a big array (gbig) ===
1569    do ig=1,nshell(ibase)
1570      gbig(:,ig+maxpw)=gshell(:,ig)
1571    end do
1572    maxpw=maxpw+nshell(ibase)
1573  end do ! ibase
1574  !
1575  ! === Compute number of G"s within each shell ===
1576  ilim=0
1577  do ibase=1,nbase
1578    ilim=ilim+nshell(ibase)
1579    shlim(ibase)=ilim
1580  end do
1581  !
1582  ! === Print out shell limits ===
1583  write(msg,'(3a)')&
1584 & ' Shells found:',ch10,&
1585 & ' number of shell    number of G vectors      cut-off energy [Ha] '
1586  call wrtout(std_out,msg,'COLL')
1587 
1588  do ibase=1,nbase
1589    write(msg,'(12x,i4,17x,i6,12x,f8.3)')ibase,shlim(ibase),two*pi**2*cnorm(ibase)
1590    call wrtout(std_out,msg,'COLL')
1591  end do
1592  write(msg,'(a)')ch10
1593  call wrtout(std_out,msg,'COLL')
1594  ABI_FREE(gshell)
1595  ABI_FREE(insort)
1596  ABI_FREE(nshell)
1597 
1598 end subroutine getfullg

m_gsphere/getkpgnorm [ Functions ]

[ Top ] [ m_gsphere ] [ Functions ]

NAME

 getkpgnorm

FUNCTION

  compute the norms of the k+G vectors

INPUTS

  gprimd(3,3)=metric tensor
  kg_k(3,npw_k)= G vectors, in reduced coordinates
  kpt(3)=k vector, in reduced coordinates
  npw_k=size of the G-vector set

OUTPUT

  kpgnorm(npw_k)=norms of the k+G vectors

PARENTS

      m_cut3d,partial_dos_fractions

CHILDREN

SOURCE

2505 subroutine getkpgnorm(gprimd,kpt,kg_k,kpgnorm,npw_k)
2506 
2507 
2508 !This section has been created automatically by the script Abilint (TD).
2509 !Do not modify the following lines by hand.
2510 #undef ABI_FUNC
2511 #define ABI_FUNC 'getkpgnorm'
2512 !End of the abilint section
2513 
2514  implicit none
2515 
2516 !Arguments ------------------------------------
2517 !scalars
2518  integer,intent(in) :: npw_k
2519 !arrays
2520  integer,intent(in) :: kg_k(3,npw_k)
2521  real(dp),intent(in) :: gprimd(3,3),kpt(3)
2522  real(dp),intent(out) :: kpgnorm(npw_k)
2523 
2524 !Local variables-------------------------------
2525 !scalars
2526  integer :: ipw
2527  real(dp) :: g11,g12,g13,g21,g22,g23,g31,g32,g33,k1,k2,k3,kpg1,kpg2,kpg3,rr,xx
2528  real(dp) :: yy,zz
2529 
2530 ! *************************************************************************
2531 
2532  k1=kpt(1) ; k2=kpt(2) ; k3=kpt(3)
2533  g11=gprimd(1,1)
2534  g12=gprimd(1,2)
2535  g13=gprimd(1,3)
2536  g21=gprimd(2,1)
2537  g22=gprimd(2,2)
2538  g23=gprimd(2,3)
2539  g31=gprimd(3,1)
2540  g32=gprimd(3,2)
2541  g33=gprimd(3,3)
2542 
2543 !Loop over all k+G
2544  do ipw=1,npw_k
2545 
2546 !  Load k+G
2547    kpg1=k1+dble(kg_k(1,ipw))
2548    kpg2=k2+dble(kg_k(2,ipw))
2549    kpg3=k3+dble(kg_k(3,ipw))
2550 
2551 !  Calculate module of k+G
2552    xx=g11*kpg1+g12*kpg2+g13*kpg3
2553    yy=g21*kpg1+g22*kpg2+g23*kpg3
2554    zz=g31*kpg1+g32*kpg2+g33*kpg3
2555    rr=sqrt(xx**2+yy**2+zz**2)
2556    kpgnorm(ipw) = rr
2557 
2558  end do ! ipw
2559 
2560 end subroutine getkpgnorm

m_gsphere/gsph_extend [ Functions ]

[ Top ] [ m_gsphere ] [ Functions ]

NAME

  gsph_extend

FUNCTION

  Construct a new gsphere_t with a larger cutoff energy
  while preserving the ordering of the first G-vectors stored in in_Gsph

INPUTS

OUTPUT

PARENTS

      setup_bse,setup_bse_interp

CHILDREN

      gsph_free,gsph_init

SOURCE

2403 subroutine gsph_extend(in_Gsph,Cryst,new_ecut,new_Gsph)
2404 
2405 
2406 !This section has been created automatically by the script Abilint (TD).
2407 !Do not modify the following lines by hand.
2408 #undef ABI_FUNC
2409 #define ABI_FUNC 'gsph_extend'
2410 !End of the abilint section
2411 
2412  implicit none
2413 
2414 !Arguments ------------------------------------
2415 !scalars
2416  real(dp),intent(in) :: new_ecut
2417  type(crystal_t),intent(in) :: Cryst
2418  type(gsphere_t),intent(in) :: in_Gsph
2419  type(gsphere_t),intent(out) :: new_Gsph
2420 !arrays
2421 
2422 !Local variables-------------------------------
2423 !scalars
2424  integer :: new_ng,in_ng,ig,ierr,sh
2425 !arrays
2426  integer,allocatable :: new_gvec(:,:)
2427 
2428 ! *********************************************************************
2429 
2430  call gsph_init(new_Gsph,Cryst,0,ecut=new_ecut)
2431 
2432  if (new_Gsph%ng > in_Gsph%ng) then
2433 
2434    new_ng = new_Gsph%ng
2435    in_ng  = in_Gsph%ng
2436 
2437    ierr = 0
2438    do ig=1,in_ng
2439      if (ANY(new_Gsph%gvec(:,ig) /= in_Gsph%gvec(:,ig)) ) then
2440        ierr = ierr + 1
2441        write(std_out,*)" new_gvec, in_gvec",ig,new_Gsph%gvec(:,ig),in_Gsph%gvec(:,ig)
2442      end if
2443    end do
2444 
2445    if (ierr==0) RETURN
2446 
2447    ierr = 0
2448    do sh=1,in_Gsph%nsh
2449      if ( new_Gsph%shlim(sh) /= in_Gsph%shlim(sh) .or. &
2450 &         ABS(new_Gsph%shlen(sh)-in_Gsph%shlen(sh)) > tol12 ) then
2451        ierr = ierr + 1
2452        write(std_out,*)"new_shlim, in_shlim",sh,new_Gsph%shlim(sh),in_Gsph%shlim(sh)
2453        write(std_out,*)"new_shlen, in_shlen",sh,new_Gsph%shlen(sh),in_Gsph%shlen(sh)
2454      end if
2455    end do
2456    ABI_CHECK(ierr==0,"Wrong shells")
2457 
2458    ABI_MALLOC(new_gvec,(3,new_ng))
2459    new_gvec = new_Gsph%gvec
2460    new_gvec(:,1:in_ng) = in_Gsph%gvec
2461 
2462    call gsph_free(new_Gsph)
2463    call gsph_init(new_Gsph,Cryst,new_ng,gvec=new_gvec)
2464    ABI_FREE(new_gvec)
2465 
2466  else
2467    ierr = 0
2468    do ig=1,MIN(new_Gsph%ng,in_Gsph%ng)
2469      if (ANY(new_Gsph%gvec(:,ig) /= in_Gsph%gvec(:,ig)) ) then
2470        ierr = ierr + 1
2471        write(std_out,*)" new_gvec, in_gvec",ig,new_Gsph%gvec(:,ig),in_Gsph%gvec(:,ig)
2472      end if
2473    end do
2474    ABI_CHECK(ierr==0,"Fatal error")
2475  end if
2476 
2477 end subroutine gsph_extend

m_gsphere/gsph_fft_tabs [ Functions ]

[ Top ] [ m_gsphere ] [ Functions ]

NAME

 gsph_fft_tabs

FUNCTION

INPUTS

  Gsph<gsphere_t>=Info on the G-sphere
  g0(3)
  mgfft=MAXVAL(ngfft(1:3))
  ngfftf(18)=Info on the FFT mesh.

OUTPUT

  use_padfft=1 if padded FFT can be used, 0 otherwise.
  gmg0_gbound(2*mgfft+8,2)=Tables for improved zero-padded FFTS. Calculated only if use_padfft==1
  gmg0_ifft(Gsph%ng)=Index of G-G0 in the FFT mesh defined by ngfft.

NOTES

  The routine will stop if any G-G0 happens to be outside the FFT box.

PARENTS

      calc_sigc_me,calc_sigx_me,cchi0,cchi0q0,cchi0q0_intraband,cohsex_me
      exc_build_block,exc_build_ham,prep_calc_ucrpa

CHILDREN

      gsph_free,gsph_init

SOURCE

515 subroutine gsph_fft_tabs(Gsph,g0,mgfft,ngfft,use_padfft,gmg0_gbound,gmg0_ifft)
516 
517 
518 !This section has been created automatically by the script Abilint (TD).
519 !Do not modify the following lines by hand.
520 #undef ABI_FUNC
521 #define ABI_FUNC 'gsph_fft_tabs'
522 !End of the abilint section
523 
524  implicit none
525 
526 !Arguments ------------------------------------
527 !scalars
528  integer,intent(in) :: mgfft
529  integer,intent(out) :: use_padfft
530  type(gsphere_t),intent(in) :: Gsph
531 !arrays
532  integer,intent(in) :: g0(3),ngfft(18)
533  integer,intent(out) :: gmg0_gbound(2*mgfft+8,2),gmg0_ifft(Gsph%ng)
534 
535 !Local variables-------------------------------
536 !scalars
537  integer :: ig,ng,ierr
538  character(len=500) :: msg
539  type(MPI_type) :: MPI_enreg_seq
540 !arrays
541  integer,allocatable :: gmg0(:,:)
542  logical,allocatable :: kg_mask(:)
543 
544 ! *************************************************************************
545 
546  if (mgfft/=MAXVAL(ngfft(1:3))) then
547    MSG_ERROR("mgfft/-MAXVAL(ngfft(1:3)")
548  end if
549 
550  ng = Gsph%ng
551 
552  ierr=0; use_padfft=0
553  ABI_MALLOC(gmg0,(3,ng))
554  do ig=1,ng
555    gmg0(:,ig) = Gsph%gvec(:,ig)-g0
556    ! Consider possible wrap around errors.
557    if ( ANY(gmg0(:,ig)>ngfft(1:3)/2) .or. ANY(gmg0(:,ig)<-(ngfft(1:3)-1)/2) ) then
558      !gmg0_ifft(ig,ig01+mg0(1)+1,ig02+mg0(2)+1,ig03+mg0(3)+1) = 0
559      write(std_out,*)" outside FFT box ",gmg0(:,ig)
560      ierr=ierr+1
561    end if
562    if (ALL(gmg0(:,ig) == 0)) use_padfft=1
563  end do
564 
565  if (ierr/=0) then
566    write(msg,'(a,i0,a)')'Found ',ierr,' G-G0 vectors falling outside the FFT box. This is not allowed '
567    MSG_ERROR(msg)
568  end if
569  !
570  ! Evaluate the tables needed for the padded FFT performed in rhotwg. Note that we have
571  ! to pass G-G0 to sphereboundary instead of G as we need FFT results on the shifted G-sphere,
572  ! If Gamma is not inside G-G0 one has to disable FFT padding as sphereboundary will give wrong tables.
573  if (use_padfft==1) then
574    call sphereboundary(gmg0_gbound,1,gmg0,mgfft,ng)
575  end if
576 
577  call initmpi_seq(MPI_enreg_seq) ! No FFT parallelism.
578  call init_distribfft_seq(MPI_enreg_seq%distribfft,'c',ngfft(2),ngfft(3),'all')
579 
580  ABI_MALLOC(kg_mask,(ng))
581  call kgindex(gmg0_ifft,gmg0,kg_mask,MPI_enreg_seq,ngfft,ng)
582 
583  ABI_CHECK(ALL(kg_mask),"FFT para not yet implemented")
584  ABI_FREE(kg_mask)
585 
586  ABI_DEALLOCATE(gmg0)
587  call destroy_mpi_enreg(MPI_enreg_seq)
588 
589 end subroutine gsph_fft_tabs

m_gsphere/gsph_free [ Functions ]

[ Top ] [ m_gsphere ] [ Functions ]

NAME

 gsph_free

FUNCTION

  Deallocate the memory in a gsphere_t data type.

INPUTS

   Gsph = datatype to be freed

PARENTS

      bethe_salpeter,cchi0,cchi0q0,gwls_hamiltonian,m_gsphere,mrgscr
      screening,sigma

CHILDREN

      gsph_free,gsph_init

SOURCE

809 subroutine gsph_free(Gsph)
810 
811 
812 !This section has been created automatically by the script Abilint (TD).
813 !Do not modify the following lines by hand.
814 #undef ABI_FUNC
815 #define ABI_FUNC 'gsph_free'
816 !End of the abilint section
817 
818  implicit none
819 
820 !Arguments ------------------------------------
821 !scalars
822  type(gsphere_t),intent(inout) :: Gsph
823 
824 ! *************************************************************************
825 
826  DBG_ENTER("COLL")
827 
828  !@gsphere_t
829 
830 ! integer arrays.
831  if (allocated(Gsph%g2sh)) then
832    ABI_FREE(Gsph%g2sh)
833  end if
834  if (allocated(Gsph%gvec)) then
835    ABI_FREE(Gsph%gvec)
836  end if
837  if (allocated(Gsph%g2mg)) then
838    ABI_FREE(Gsph%g2mg)
839  end if
840  if (allocated(Gsph%rottb)) then
841    ABI_FREE(Gsph%rottb)
842  end if
843  if (allocated(Gsph%rottbm1)) then
844    ABI_FREE(Gsph%rottbm1)
845  end if
846  if (allocated(Gsph%shlim)) then
847    ABI_FREE(Gsph%shlim)
848  end if
849 
850  if (allocated(Gsph%shlen)) then
851    ABI_FREE(Gsph%shlen)
852  end if
853 
854 ! complex arrays
855  if (allocated(Gsph%phmGt)) then
856    ABI_FREE(Gsph%phmGt)
857  end if
858  if (allocated(Gsph%phmSGt)) then
859    ABI_FREE(Gsph%phmSGt)
860  end if
861 
862  DBG_EXIT("COLL")
863 
864 end subroutine gsph_free

m_gsphere/gsph_g_idx [ Functions ]

[ Top ] [ m_gsphere ] [ Functions ]

NAME

 gsph_g_idx

FUNCTION

 Return the index of G in the sphere. zero if not in the sphere

INPUTS

  Gsph<gsphere_t>=Info on the G-sphere
  gg(3)=Reduced coordinates of the G-vector.

NOTES

  The function assumes that the G-vectors are ordered with increasing length.

PARENTS

SOURCE

887 pure function gsph_g_idx(Gsph,gg) result(g_idx)
888 
889 
890 !This section has been created automatically by the script Abilint (TD).
891 !Do not modify the following lines by hand.
892 #undef ABI_FUNC
893 #define ABI_FUNC 'gsph_g_idx'
894 !End of the abilint section
895 
896  implicit none
897 
898 !Arguments ------------------------------------
899 !scalars
900  type(gsphere_t),intent(in) :: Gsph
901  integer :: g_idx
902 !arrays
903  integer,intent(in) :: gg(3)
904 
905 !Local variables-------------------------------
906 !scalars
907  integer :: ishbsc,igs,ige
908  real(dp) :: glen
909  logical :: found
910 
911 ! *************************************************************************
912 
913  ! * Use shells and bisect to find the star and stop index thus avoiding the storage of a table (ig1,ig2)
914  glen = two_pi*SQRT(DOT_PRODUCT(gg,MATMUL(Gsph%gmet,gg)))
915 
916  ishbsc = bisect(Gsph%shlen,glen)
917  if ( ANY(ishbsc==(/0,Gsph%nsh/)) ) then ! glen out of range.
918    g_idx=0; RETURN
919  end if
920 
921  igs = Gsph%shlim(ishbsc)
922  ige = Gsph%shlim(MIN(ishbsc+2,Gsph%nsh+1))-1
923 
924  g_idx=igs-1; found=.FALSE.
925  do while (.not.found .and. g_idx<ige)
926    g_idx=g_idx+1
927    found=(ALL(Gsph%gvec(:,g_idx)==gg(:)))
928  end do
929  if (.not.found) g_idx=0
930 
931 end function gsph_g_idx

m_gsphere/gsph_gmg_fftidx [ Functions ]

[ Top ] [ m_gsphere ] [ Functions ]

NAME

 gsph_gmg_fftidx

FUNCTION

 Return the index of G1-G2 in the FFT mesh defined by ngfft. zero if not found.

INPUTS

  Gsph<gsphere_t>=Info on the G-sphere
  ig1,ig2 index of g1 and g2 in the G-sphere.
  ngfft(18)=Info on the FFT mesh.

PARENTS

SOURCE

1023 pure function gsph_gmg_fftidx(Gsph,ig1,ig2,ngfft) result(fft_idx)
1024 
1025 
1026 !This section has been created automatically by the script Abilint (TD).
1027 !Do not modify the following lines by hand.
1028 #undef ABI_FUNC
1029 #define ABI_FUNC 'gsph_gmg_fftidx'
1030 !End of the abilint section
1031 
1032  implicit none
1033 
1034 !Arguments ------------------------------------
1035 !scalars
1036  type(gsphere_t),intent(in) :: Gsph
1037  integer,intent(in) :: ig1,ig2
1038  integer :: fft_idx
1039 !arrays
1040  integer,intent(in) :: ngfft(18)
1041 
1042 !Local variables-------------------------------
1043 !scalars
1044  integer :: id1,id2,id3
1045 !arrays
1046  integer :: g1mg2(3)
1047 
1048 ! *************************************************************************
1049 
1050  g1mg2(:)=Gsph%gvec(:,ig1)-Gsph%gvec(:,ig2)
1051 
1052  ! Make sure G1-G2 is still in the FFT mesh.
1053  ! MODULO wraps G1-G2 in the FFT box but the Fourier components are not periodic!
1054  if (ANY(g1mg2(:)>ngfft(1:3)/2) .or. ANY(g1mg2(:)<-(ngfft(1:3)-1)/2)) then
1055    fft_idx=0; RETURN
1056  end if
1057 
1058  id1=MODULO(g1mg2(1),ngfft(1))
1059  id2=MODULO(g1mg2(2),ngfft(2))
1060  id3=MODULO(g1mg2(3),ngfft(3))
1061  fft_idx= 1 + id1 + id2*ngfft(1) + id3*ngfft(1)*ngfft(2)
1062 
1063 end function gsph_gmg_fftidx

m_gsphere/gsph_gmg_idx [ Functions ]

[ Top ] [ m_gsphere ] [ Functions ]

NAME

 gsph_gmg_idx

FUNCTION

 Return the index of G1-G2 in the sphere. zero if not in the sphere

INPUTS

  Gsph<gsphere_t>=Info on the G-sphere
  ig1,ig2 index of g1 and g2 in the G-sphere.

NOTES

  The function assumes that the G-vectors are ordered with increasing length.

PARENTS

SOURCE

 954 pure function gsph_gmg_idx(Gsph,ig1,ig2) result(ig1mg2)
 955 
 956 
 957 !This section has been created automatically by the script Abilint (TD).
 958 !Do not modify the following lines by hand.
 959 #undef ABI_FUNC
 960 #define ABI_FUNC 'gsph_gmg_idx'
 961 !End of the abilint section
 962 
 963  implicit none
 964 
 965 !Arguments ------------------------------------
 966 !scalars
 967  type(gsphere_t),intent(in) :: Gsph
 968  integer,intent(in) :: ig1,ig2
 969  integer :: ig1mg2
 970 
 971 !Local variables-------------------------------
 972 !scalars
 973  integer :: ishbsc,igs,ige
 974  real(dp) :: difflen
 975  logical :: found
 976 !arrays
 977  integer :: g1mg2(3)
 978 
 979 ! *************************************************************************
 980 
 981  g1mg2 = Gsph%gvec(:,ig1)-Gsph%gvec(:,ig2)
 982 
 983  ! * Use shells and bisect to find the star and stop index thus avoiding the storage of a table (ig1,ig2)
 984  difflen = two_pi*SQRT(DOT_PRODUCT(g1mg2,MATMUL(Gsph%gmet,g1mg2)))
 985 
 986  ! FIXME It seems bisect is not portable, on my laptop test v5/t72 the number of skipped G-vectors is > 0
 987  ishbsc = bisect(Gsph%shlen,difflen)
 988  if ( ANY(ishbsc==(/0,Gsph%nsh/)) ) then ! difflen out of range.
 989    ig1mg2=0; RETURN
 990  end if
 991 
 992  igs = Gsph%shlim(ishbsc)
 993  ige = Gsph%shlim(MIN(ishbsc+2,Gsph%nsh+1))-1
 994 
 995  ig1mg2=igs-1; found=.FALSE.
 996  do while (.not.found .and. ig1mg2<ige)
 997    ig1mg2=ig1mg2+1
 998    found=(ALL(Gsph%gvec(:,ig1mg2)==g1mg2(:)))
 999  end do
1000  if (.not.found) ig1mg2=0
1001 
1002 end function gsph_gmg_idx

m_gsphere/gsph_in_fftbox [ Functions ]

[ Top ] [ m_gsphere ] [ Functions ]

NAME

 gsph_in_fftbox

FUNCTION

  Initialize the largest Gsphere contained in the FFT box.

INPUTS

  Cryst<crystal_t> = Info on unit cell and its symmetries.
  ngfft(18)=Info on the FFT box.

OUTPUT

  Gsph<gsphere_t>=Data type containing information related to the set of G vectors
   completetly initialized in output.

PARENTS

      cchi0,cchi0q0

CHILDREN

      gsph_free,gsph_init

SOURCE

617 subroutine gsph_in_fftbox(Gsph,Cryst,ngfft)
618 
619 
620 !This section has been created automatically by the script Abilint (TD).
621 !Do not modify the following lines by hand.
622 #undef ABI_FUNC
623 #define ABI_FUNC 'gsph_in_fftbox'
624 !End of the abilint section
625 
626  implicit none
627 
628 !Arguments ------------------------------------
629 !scalars
630  type(crystal_t),intent(in) :: Cryst
631  type(gsphere_t),intent(out) :: Gsph
632 !arrays
633  integer,intent(in) :: ngfft(18)
634 
635 !Local variables-------------------------------
636 !scalars
637  integer :: dir1,dir2,dir3,npw,ig,ist
638  real(dp) :: ecut,trial_ene
639 !arrays
640  integer :: n1_max(3),n2_max(3),n3_max(3),vec(3)
641  integer,allocatable :: gvec(:,:)
642 
643 !************************************************************************
644 
645  ! Find ecut for the largest G-sphere contained in the FFT box.
646  n1_max(1) = -(ngfft(1)-1)/2
647  n2_max(1) = -(ngfft(2)-1)/2
648  n3_max(1) = -(ngfft(3)-1)/2
649 
650  n1_max(2) = 0
651  n2_max(2) = 0
652  n3_max(2) = 0
653 
654  n1_max(3) = ngfft(1)/2
655  n2_max(3) = ngfft(2)/2
656  n3_max(3) = ngfft(3)/2
657 
658  ecut = HUGE(one)
659  do dir3=1,3
660    vec(3) = n1_max(dir3)
661    do dir2=1,3
662      vec(2) = n2_max(dir2)
663      do dir1=1,3
664        vec(1) = n1_max(dir1)
665        if (ANY(vec/=0)) then
666          trial_ene = half * normv(vec,Cryst%gmet,"G")**2
667          ecut = MIN(ecut,trial_ene)
668          !write(std_out,*)vec(:),trial_ene
669        end if
670      end do
671    end do
672  end do
673  !
674  ! Init sphere from ecut.
675  call gsph_init(Gsph,Cryst,0,ecut=ecut)
676  !
677  ! Make sure that Gsph does not contain G vectors outside the FFT box.
678  ! kpgsph might return G whose energy is larger than the input ecut.
679  npw = Gsph%ng
680  star_loop: do ist=1,Gsph%nsh-1
681    do ig=Gsph%shlim(ist),Gsph%shlim(ist+1)
682      if ( ANY(Gsph%gvec(:,ig)>ngfft(1:3)/2) .or. ANY(Gsph%gvec(:,ig)<-(ngfft(1:3)-1)/2) ) then
683        npw = Gsph%shlim(ist)-1  ! Gsph exceeds the FFT box. Only the shells up to npw will be used.
684        EXIT star_loop
685      end if
686    end do
687  end do star_loop
688 
689  if (npw<Gsph%ng) then
690    MSG_COMMENT("Have to reinit Gpshere")
691    ABI_MALLOC(gvec,(3,npw))
692    gvec =Gsph%gvec(:,1:npw)
693    call gsph_free(Gsph)
694    call gsph_init(Gsph,Cryst,npw,gvec=gvec)
695    ABI_FREE(gvec)
696  end if
697 
698 end subroutine gsph_in_fftbox

m_gsphere/gsph_init [ Functions ]

[ Top ] [ m_gsphere ] [ Functions ]

NAME

 gsph_init

FUNCTION

  Main creation method for the Gvectors data type

INPUTS

  Cryst<crystal_t> = Info on unit cell and its symmetries
     %nsym=number of symmetry operations
     %symrec(3,3,nsym)=symmetry operations in reciprocal space
     %tnons(3,nsym)=fractional translations
     %gmet(3,3)=reciprocal space metric (bohr**-2).
     %gprimd(3,3)=dimensional reciprocal space primitive translations
  ng=number of G vectors, needed only if gvec is passed.
  [gvec(3,ng)]=coordinates of G vectors
  [ecut]=Cutoff energy for G-sphere. gvec and ecut are mutually exclusive.

OUTPUT

  Gsph<gsphere_t>=Data type containing information related to the set of G vectors
   completetly initialized in output.

NOTES

  gvec are supposed to be ordered with increasing norm.

PARENTS

      gwls_hamiltonian,m_gsphere,mrgscr,setup_bse,setup_bse_interp
      setup_screening,setup_sigma

CHILDREN

      gsph_free,gsph_init

SOURCE

308 subroutine gsph_init(Gsph,Cryst,ng,gvec,ecut)
309 
310 
311 !This section has been created automatically by the script Abilint (TD).
312 !Do not modify the following lines by hand.
313 #undef ABI_FUNC
314 #define ABI_FUNC 'gsph_init'
315 !End of the abilint section
316 
317  implicit none
318 
319 !Arguments ------------------------------------
320 !scalars
321  integer,intent(in) :: ng
322  real(dp),optional,intent(in) :: ecut
323  type(crystal_t),target,intent(in) :: Cryst
324  type(gsphere_t),intent(out) :: Gsph
325 !arrays
326  integer,optional,intent(in) :: gvec(3,ng)
327 !Local variables-------------------------------
328 !scalars
329  integer,parameter :: nkpt1=1
330  integer :: ig,isearch,img,ish,isym,nsh,nsym,timrev,pinv,g1,g2,g3,ss,ee
331  real(dp) :: eps,norm,norm_old,max_ecut,gsq
332 !arrays
333  real(dp),parameter :: k_gamma(3)=(/zero,zero,zero/)
334  integer :: sg(3),gsearch(3)
335  integer,allocatable :: shlim(:)
336  integer,pointer :: symrec(:,:,:),gvec_ptr(:,:)
337  real(dp) :: kptns1(3,nkpt1)
338  real(dp),allocatable :: shlen(:)
339  real(dp),pointer :: tnons(:,:)
340 
341 !************************************************************************
342 
343  DBG_ENTER("COLL")
344 
345  !@gsphere_t
346  ! === Copy info on symmetries ===
347  nsym   =  Cryst%nsym
348  timrev =  Cryst%timrev
349  symrec => Cryst%symrec
350  tnons  => Cryst%tnons
351  !
352  ! === Initialize the object ===
353  Gsph%istwfk = 1           ! Time reversal is not used here.
354  Gsph%nsym   = nsym
355  Gsph%timrev = timrev
356 
357  Gsph%gmet   = Cryst%gmet
358  Gsph%gprimd = Cryst%gprimd
359 
360  if (PRESENT(gvec)) then
361    if (PRESENT(ecut)) then
362      MSG_BUG("ecut cannot be present when gvec is used")
363    end if
364    Gsph%ng= ng
365    ABI_MALLOC(Gsph%gvec,(3,ng))
366    Gsph%gvec=gvec
367    !
368    ! Calculate cutoff energy.of the sphere.
369    max_ecut=-one
370    do ig=1,ng
371      g1=gvec(1,ig)
372      g2=gvec(2,ig)
373      g3=gvec(3,ig)
374      gsq=       Cryst%gmet(1,1)*g1**2+Cryst%gmet(2,2)*g2**2+Cryst%gmet(3,3)*g3**2+ &
375 &          two*(Cryst%gmet(1,2)*g1*g2+Cryst%gmet(1,3)*g1*g3+Cryst%gmet(2,3)*g2*g3)
376      max_ecut=MAX(max_ecut,gsq)
377    end do
378    max_ecut=two*max_ecut*pi**2
379    Gsph%ecut= max_ecut
380 
381  else
382    ! To be consistent with the previous implementation.
383    MSG_WARNING("Init from ecut has to be tested")
384    !call setshells(ecut,npw,nsh,nsym,Cryst%gmet,Cryst%gprimd,Cryst%symrel,tag,Cryst%ucvol)
385    Gsph%ecut = ecut
386    pinv=+1; kptns1(:,1)=k_gamma
387    call merge_and_sort_kg(nkpt1,kptns1,ecut,Cryst%nsym,pinv,Cryst%symrel,Cryst%gprimd,gvec_ptr,0)
388    Gsph%ng = SIZE(gvec_ptr,DIM=2)
389    ABI_MALLOC(Gsph%gvec, (3,Gsph%ng))
390    Gsph%gvec = gvec_ptr
391    ABI_FREE(gvec_ptr)
392  end if
393  !
394  ! === Calculate phase exp{-i2\pi G.\tau} ===
395  ABI_MALLOC(Gsph%phmGt,(Gsph%ng,nsym))
396  do ig=1,Gsph%ng
397    do isym=1,nsym
398     Gsph%phmGt(ig,isym)=EXP(-j_dpc*two_pi*DOT_PRODUCT(Gsph%gvec(:,ig),tnons(:,isym)))
399    end do
400  end do
401  !
402  ! === Calculate phase phsgt= exp{-i2\pi SG\cdot t} ===
403  ! TODO Here we can store only one of this arrays but I have to rewrite screeening!
404  ABI_MALLOC(Gsph%phmSGt,(Gsph%ng,nsym))
405  do ig=1,Gsph%ng
406    do isym=1,nsym
407      sg=MATMUL(symrec(:,:,isym),Gsph%gvec(:,ig))
408      Gsph%phmSGt(ig,isym)=EXP(-j_dpc*two_pi*DOT_PRODUCT(sg,tnons(:,isym)))
409    end do
410  end do
411  !
412  ! === Calculate number of shells and corresponding starting index ===
413  ! * Shells are useful to speed up search algorithms see e.g setup_G_rotation.
414  ! * The last shell ends at ng+1, thus gvec is supposed to be closed.
415 
416  ABI_CHECK(ALL(Gsph%gvec(1:3,1)==0),'First G must be 0')
417 
418  ABI_MALLOC(Gsph%g2sh,(Gsph%ng))
419  Gsph%g2sh(1)=1 ! This table is useful if we dont loop over shell
420 
421  ! For each shell, gives the index of the initial G-vector.
422  ABI_MALLOC(shlim,(Gsph%ng+1))
423  shlim(1)=1
424 
425  ! For each shell, gives the radius of the shell.
426  ABI_MALLOC(shlen,(Gsph%ng))
427  shlen(1)=zero
428 
429  nsh=1; norm_old=zero
430  do ig=2,Gsph%ng
431    norm=two_pi*SQRT(DOT_PRODUCT(Gsph%gvec(:,ig),MATMUL(Cryst%gmet,Gsph%gvec(:,ig))))
432    eps=norm*tol8
433    if (ABS(norm-norm_old)>eps) then
434      norm_old = norm; nsh = nsh + 1
435      shlim(nsh)=ig
436      shlen(nsh)=norm
437    end if
438    Gsph%g2sh(ig)=nsh
439  end do
440  shlim(nsh+1)=Gsph%ng+1
441 
442  ! === Save info on the shells ===
443  Gsph%nsh=nsh
444  ABI_MALLOC(Gsph%shlim,(nsh+1))
445  Gsph%shlim=shlim(1:nsh+1)
446  ABI_MALLOC(Gsph%shlen,(nsh  ))
447  Gsph%shlen=shlen(1:nsh)
448  ABI_FREE(shlim)
449  ABI_FREE(shlen)
450  !
451  ! === Calculate tables for rotated G"s ===
452  ABI_MALLOC(Gsph%rottb  ,(Gsph%ng,timrev,nsym))
453  ABI_MALLOC(Gsph%rottbm1,(Gsph%ng,timrev,nsym))
454 
455  call setup_G_rotation(nsym,symrec,timrev,Gsph%ng,Gsph%gvec,&
456 &  Gsph%g2sh,Gsph%nsh,Gsph%shlim,Gsph%rottb,Gsph%rottbm1)
457 
458  ! Store Mapping G --> -G
459  ! (we use a specialized table instead of rootb since rottb assumes time-reversal symmetry.
460  ABI_MALLOC(gsph%g2mg, (gsph%ng))
461 
462  do ig=1,gsph%ng
463    ish=gsph%g2sh(ig)
464    ss=gsph%shlim(ish); ee=gsph%shlim(ish+1)-1
465    gsearch = -gsph%gvec(:,ig)
466    img = 0
467    ! Loop on shells to speed up the search.
468    do isearch=ss,ee
469      if (all(gsph%gvec(:,isearch) == gsearch)) then
470        img = isearch; exit
471      end if
472    end do
473    if (img==0) MSG_ERROR("Cannot find -G in G-sphere!")
474    gsph%g2mg(ig) = img
475  end do
476 
477  !call print_gsphere(Gsph,unit=std_out,prtvol=1)
478 
479  DBG_EXIT("COLL")
480 
481 end subroutine gsph_init

m_gsphere/gsphere_t [ Types ]

[ Top ] [ m_gsphere ] [ Types ]

NAME

 gsphere_t

FUNCTION

 The gsphere_t data type contains information related to the set of G vectors
 used during a screening or a GW calculation, as well as symmetry tables relating
 these vectors. Presently the following quantities are stored

 1) The reduced coordinates of the G vectors (arrays gvec)
 2) Tables giving the correspondence between a G-vector and its rotated image
    through a symmetry operation in reciprocal space.
 3) List of the irreducible G pairs
 4) Tables giving, for each pair in the full reciprocal space, the corresponding
    irreducible pair as well as the symmetry operation in reciprocal space

 Note that, unlike the GS part, the basis set does not depend on the k-point.

NOTES

 To indicate the indices in the arrays grottb, grottbm1 we use the following notation :

  g defines the index of the reciprocal lattice vector in the array gvec
  s  indicates the index of the symmetry operation in reciprocal space
  i  can  be one or two. 1 is used to indicate the identity operator

SOURCE

 91  type,public :: gsphere_t
 92 
 93   integer :: ng
 94   ! Total number of G vectors in the sphere taking into account umklapps
 95   ! it defines the size of the array gvec and it accounts for possible umklapps for which
 96   ! one has to shift the sphere.
 97 
 98   ! TODO: The sphere should be enlarged in gpairs_init using mg0 and (ecut|ng) as input.
 99   ! For the time being we keep the old implementation.
100   !%integer :: ng_eff
101   ! Effective number of G vectors, i.e. the number of G in the smaller sphere without umklapps.
102   ! ng_eff<=ng and should be used to loop over the elements of (chi0|epsilon|W).
103   !
104   ! TODO: Add info on FFT including zero-padding algorithm.
105   ! table must be recalculated for each G0 and rho_tw_g should accept gpshere in input.
106 
107   integer :: nsh                   ! Number of shells
108   integer :: nsym                  ! The number of symmetry operations
109   integer :: timrev                ! 2 if time-reversal is used, 1 otherwise
110   integer :: istwfk=1              ! Storage mode. At present time-reversal is not used.
111 
112   !integer :: mg0(3)=0
113   ! For each reduced direction gives the max G0 component to account for umklapp processes.
114 
115   real(dp) :: ecut
116   ! Cutoff energy of the sphere.
117 
118   real(dp) :: gmet(3,3)
119   ! Reciprocal space metric ($\textrm{bohr}^{-2}$).
120 
121   real(dp) :: gprimd(3,3)
122   ! Dimensional reciprocal space primitive translations (bohr^{-1})
123 
124   integer,allocatable :: g2sh(:)
125   ! g2sh(ng)
126   ! For each G, it gives the index of the shell to which it belongs.
127 
128   integer,allocatable :: gvec(:,:)
129   ! gvec(3,ng)
130   ! Reduced coordinates of G vectors.
131 
132   integer,allocatable :: g2mg(:)
133   ! g2mg(ng)
134   ! Correspondence G --> -G
135 
136   integer,allocatable :: rottb(:,:,:)
137   ! rottb(ng,timrev,nsym)
138   ! rottb(G,I,S) is the index of (SI) G in the array gvec
139   ! where I is either the identity or the inversion.
140 
141   integer,allocatable :: rottbm1(:,:,:)
142   ! rottb(ng,timrev,nsym)
143   ! rottbm1(G,I,S) is the index of IS{^-1} G in the array gvec
144 
145   integer,allocatable :: shlim(:)
146   ! shlim(nsh+1)
147   ! Index of the first G vector in each shell, =ng+1 for nsh+1
148 
149   real(dp),allocatable :: shlen(:)
150   ! shlen(nsh)
151   ! Radius of each shell.
152 
153   !TODO switch to dpc
154   complex(gwpc),allocatable :: phmGt(:,:)
155   ! phmGt(ng,nsym)
156   ! Phase factor e^{-i2\pi(G.\tau)} where $\tau$ is the fractional translation associated to isym.
157 
158   complex(gwpc),allocatable :: phmSGt(:,:)
159   ! phmSGt(ng,nsym)
160   ! Phase factor e^{-i2\pi(SG.\tau)} where S is one of the symmetry properties in reciprocal space.
161 
162  end type gsphere_t
163 
164  public :: gsph_init          ! Initialize the G-sphere.
165  public :: gsph_fft_tabs      ! Returns useful tables for FFT (with or without padding).
166  public :: gsph_in_fftbox     ! Initialize the largest Gsphere contained in the FFT box.
167  public :: print_gsphere      ! Printout of basic dimensions.
168  public :: gsph_free          ! Free memory allocated in the object.
169  public :: gsph_g_idx         ! Returns the index of G from its reduced coordinates.
170  public :: gsph_gmg_idx       ! Returns the index of G1-G2 from their indeces
171  public :: gsph_gmg_fftidx    ! Returns the index of G1-G2 in the FFT mesh defined by ngfft.
172  public :: gsph_extend        ! Construct a new gsphere_t with a larger cutoff energy

m_gsphere/kg_map [ Functions ]

[ Top ] [ m_gsphere ] [ Functions ]

NAME

  kg_map

FUNCTION

  Compute the mapping between two lists of g-vectors.

INPUTS

   npw1, kg1(3,npw1)=First list of G-vectors
   npw2, kg2(3,npw2)=Second list of G-vectors

OUTPUT

   g2g1(npw2) = Mapping kg2 index --> kg1 index.
                Set to 0 if kg2(:,ig) not in kg1
   nmiss = Number of G-vectors in kg2 not found in kg1

PARENTS

      m_wfd

CHILDREN

      gsph_free,gsph_init

SOURCE

2132 subroutine kg_map(npw1,kg1,npw2,kg2,g2g1,nmiss)
2133 
2134 
2135 !This section has been created automatically by the script Abilint (TD).
2136 !Do not modify the following lines by hand.
2137 #undef ABI_FUNC
2138 #define ABI_FUNC 'kg_map'
2139 !End of the abilint section
2140 
2141  implicit none
2142 
2143 !Arguments ------------------------------------
2144 !scalars
2145  integer,intent(in) :: npw1,npw2
2146  integer,intent(out) :: nmiss
2147 !arrays
2148  integer,intent(in) :: kg1(3,npw1),kg2(3,npw2)
2149  integer,intent(out) :: g2g1(npw2)
2150 
2151 !Local variables ------------------------------
2152 !scalars
2153  integer :: ii,ipw,i1,i2,i3,n1,n2,n3
2154 !arrays
2155  integer :: gmax(3),g1_max(3),g2_max(3)
2156  integer,allocatable :: iwork(:,:,:)
2157 
2158 !************************************************************************
2159 
2160  g1_max = maxval(abs(kg1))
2161  g2_max = maxval(abs(kg2))
2162  do ii=1,3
2163    gmax(ii) = max(g1_max(ii), g2_max(ii))
2164  end do
2165  gmax = 2*gmax + 1
2166  n1 = gmax(1); n2 = gmax(2); n3 = gmax(3)
2167 
2168  ABI_MALLOC(iwork, (n1, n2, n3))
2169 
2170  ! Insert kg1 into work with extra 0 s around outside:
2171  iwork = 0
2172  do ipw=1,npw1
2173    i1 = kg1(1,ipw); if (i1<0) i1=i1+n1; i1=i1+1
2174    i2 = kg1(2,ipw); if (i2<0) i2=i2+n2; i2=i2+1
2175    i3 = kg1(3,ipw); if (i3<0) i3=i3+n3; i3=i3+1
2176    iwork(i1,i2,i3) = ipw
2177  end do
2178 
2179  g2g1 = 0; nmiss = 0
2180  do ipw=1,npw2
2181    i1 = kg2(1,ipw); if (i1<0) i1=i1+n1; i1=i1+1
2182    i2 = kg2(2,ipw); if (i2<0) i2=i2+n2; i2=i2+1
2183    i3 = kg2(3,ipw); if (i3<0) i3=i3+n3; i3=i3+1
2184    g2g1(ipw) = iwork(i1,i2,i3)
2185    if (g2g1(ipw) == 0) nmiss = nmiss + 1
2186  end do
2187 
2188  ABI_FREE(iwork)
2189 
2190 end subroutine kg_map

m_gsphere/make_istwk_table [ Functions ]

[ Top ] [ m_gsphere ] [ Functions ]

NAME

 make_istwfk_table

FUNCTION

INPUTS

  ng1,ng2,ng3

OUTPUT

NOTES

   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

      gsph_free,gsph_init

SOURCE

2219 subroutine make_istwfk_table(istwf_k,ng1,ng2,ng3,ig1_inver,ig2_inver,ig3_inver)
2220 
2221 
2222 !This section has been created automatically by the script Abilint (TD).
2223 !Do not modify the following lines by hand.
2224 #undef ABI_FUNC
2225 #define ABI_FUNC 'make_istwfk_table'
2226 !End of the abilint section
2227 
2228  implicit none
2229 
2230 !Arguments ------------------------------------
2231 !scalars
2232  integer,intent(in) :: ng1,ng2,ng3,istwf_k
2233 !arrays
2234  integer,intent(out) :: ig1_inver(ng1),ig2_inver(ng2),ig3_inver(ng3)
2235 
2236 !Local variables ------------------------------
2237 !scalars
2238  integer :: i1,i2,i3
2239  character(len=500) :: msg
2240 
2241 !************************************************************************
2242 
2243 ! Initialize the inverse coordinates
2244  select case (istwf_k)
2245 
2246  case (1)
2247    ig1_inver(1)=1
2248    do i1=2,ng1
2249      ig1_inver(i1)=ng1+2-i1
2250    end do
2251    ig2_inver(1)=1
2252    do i2=2,ng2
2253      ig2_inver(i2)=ng2+2-i2
2254    end do
2255    ig3_inver(1)=1
2256    do i3=2,ng3
2257      ig3_inver(i3)=ng3+2-i3
2258    end do
2259 
2260  case (2:8)
2261    if (istwf_k==2 .or. istwf_k==4 .or. istwf_k==6 .or. istwf_k==8) then
2262      ig1_inver(1)=1
2263      do i1=2,ng1
2264        ig1_inver(i1)=ng1+2-i1
2265      end do
2266    else
2267      do i1=1,ng1
2268        ig1_inver(i1)=ng1+1-i1
2269      end do
2270    end if
2271    if (istwf_k>=2 .and. istwf_k<=5) then
2272      ig2_inver(1)=1
2273      do i2=2,ng2
2274        ig2_inver(i2)=ng2+2-i2
2275      end do
2276    else
2277      do i2=1,ng2
2278        ig2_inver(i2)=ng2+1-i2
2279      end do
2280    end if
2281    if (istwf_k==2 .or. istwf_k==3 .or. istwf_k==6 .or. istwf_k==7) then
2282      ig3_inver(1)=1
2283      do i3=2,ng3
2284        ig3_inver(i3)=ng3+2-i3
2285      end do
2286    else
2287      do i3=1,ng3
2288        ig3_inver(i3)=ng3+1-i3
2289      end do
2290    end if
2291 
2292  case default
2293    write(msg,'(a,i0)')" Wrong value for istwf_k: ",istwf_k
2294    MSG_ERROR(msg)
2295  end select
2296 
2297 end subroutine make_istwfk_table

m_gsphere/merge_and_sort_kg [ Functions ]

[ Top ] [ m_gsphere ] [ Functions ]

NAME

  merge_and_sort_kg

FUNCTION

  This routine merges a set of k-centered G-spheres of cutoff energy ecut and
  returns a Gamma-centered G-spheres. The elements in the final G-spheres are packed with increasing module.

INPUTS

  nkpt=Number of k-points
  kptns(3,nkpt)=The k-points in reduced coordinates defining the k-centered G-spheres.
  ecut=Cutoff energy for the k-centered G-spheres.
  nsym2=Number of symmetry operations.
  pinv=-1 if time-reversal can be used, 1 otherwise
  symrel2(3,3,nsym2)=symmetry operations in real space.
  gprimd(3,3)=dimensional primitive translations for reciprocal space ($\textrm{bohr}^{-1}$)
  prtvol=Flag defining the verbosity level.

SIDE EFFECTS

  gbig(:,:)
    in input : pointer to NULL
    in output: gbig(3,1:npw) contains the set of G-vectors ordered by shell obtained by
               merging the k-centered sphere.
  shlim_p(:)
    in input : pointer to NULL
    in output: shlim_p(nbase)=Cumulative number of G-vectors for each shell.
               where nbase is the number of irreducible G"s found.

PARENTS

      m_gsphere,m_io_kss,outkss,setup_sigma

CHILDREN

      gsph_free,gsph_init

SOURCE

1198 subroutine merge_and_sort_kg(nkpt,kptns,ecut,nsym2,pinv,symrel2,gprimd,gbig,prtvol,shlim_p)
1199 
1200 
1201 !This section has been created automatically by the script Abilint (TD).
1202 !Do not modify the following lines by hand.
1203 #undef ABI_FUNC
1204 #define ABI_FUNC 'merge_and_sort_kg'
1205 !End of the abilint section
1206 
1207  implicit none
1208 
1209 !Arguments ------------------------------------
1210 !scalars
1211  integer,intent(in) :: nkpt,nsym2,pinv,prtvol
1212  real(dp),intent(in) :: ecut
1213 !arrays
1214  integer,intent(in) :: symrel2(3,3,nsym2)
1215  real(dp),intent(in) :: kptns(3,nkpt),gprimd(3,3)
1216  integer,pointer :: gbig(:,:)
1217  integer,optional,pointer :: shlim_p(:)
1218 
1219 !Local variables-------------------------------
1220 !scalars
1221  integer,parameter :: mkmem_=1
1222  integer :: ikg,ig,ikpt,nbase,sizepw,in,maxpw,is,iinv,ish,ilim,mpw
1223  integer :: exchn2n3d,istwf_k,onpw_k,ierr,npw_k,ii,isym,sizeold
1224  logical :: found
1225  character(len=500) :: msg
1226  type(MPI_type) :: MPI_enreg_seq
1227 !arrays
1228  integer :: gcur(3),geq(3)
1229  integer :: symrec2t(3,3,nsym2)
1230  integer :: dum_kg(3,0)
1231  integer,allocatable :: gbase(:,:),gbasek(:,:,:)
1232  integer,allocatable :: gcurr(:,:),gshell(:,:),insort(:),gtmp(:,:)
1233  integer,allocatable :: nbasek(:),nshell(:),shlim(:)
1234  integer,allocatable :: npwarr(:)
1235  real(dp) :: kpoint(3),gmet(3,3)
1236  real(dp),allocatable :: cnorm(:),cnormk(:,:),ctmp(:)
1237 
1238 ! *********************************************************************
1239 
1240  ! * Fake MPI_type for the sequential part.
1241  ! This routine should not be parallelized as communicating gbig and other
1242  ! tables takes more time than recalculating them in sequential.
1243  call initmpi_seq(MPI_enreg_seq)
1244 
1245 !Compute reciprocal space metrics
1246  do ii=1,3
1247    gmet(ii,:)=gprimd(1,ii)*gprimd(1,:)+&
1248 &   gprimd(2,ii)*gprimd(2,:)+&
1249 &   gprimd(3,ii)*gprimd(3,:)
1250  end do
1251 
1252  ! * Here we use TRANSPOSE(symrel2) instead of the more intuitive symrel2^{-1t} for historical reasons
1253  ! It does not affect the results since in the code below we only check the module of G
1254  do isym=1,nsym2
1255    symrec2t(:,:,isym)=TRANSPOSE(symrel2(:,:,isym))
1256  end do
1257  !
1258  ! ==============================================
1259  ! ==== Find irreducible G-vectors at each k ====
1260  ! ==============================================
1261 
1262  ABI_MALLOC(npwarr,(nkpt))
1263  exchn2n3d=0; ikg=0
1264  do ikpt=1,nkpt
1265    kpoint=kptns(:,ikpt); istwf_k=1
1266    call kpgsph(ecut,exchn2n3d,gmet,ikg,0,istwf_k,dum_kg,kpoint,0,MPI_enreg_seq,0,npwarr(ikpt))
1267  end do
1268  mpw = MAXVAL(npwarr)
1269 
1270  ABI_MALLOC(nbasek,(nkpt))
1271  ABI_MALLOC(gbasek,(3,mpw,nkpt))
1272  ABI_MALLOC(cnormk,(mpw,nkpt))
1273  nbasek=0     ! # of irreducible G at each k.
1274  cnormk=zero  ! Norm of each irreducible G.
1275  gbasek=0     ! The set of irreducible G"s at each k.
1276 
1277  do ikpt=1,nkpt
1278 
1279    kpoint = kptns(:,ikpt)
1280    npw_k  = npwarr(ikpt)
1281 
1282    exchn2n3d=0; ikg=0; istwf_k=1
1283    ABI_MALLOC(gcurr,(3,npw_k))
1284    call kpgsph(ecut,exchn2n3d,gmet,ikg,0,istwf_k,gcurr,kpoint,mkmem_,MPI_enreg_seq,npw_k,onpw_k)
1285 
1286    if (ANY(gcurr(:,1)/=0)) then
1287      MSG_BUG("gcurr(:,1)/=0")
1288    end if
1289    !
1290    ! * Search for the G"s generating the others by symmetry.
1291    !  NB: Here we use symrec2t=TRANSPOSE(symrel2) for historical reasons, see note above
1292    call get_irredg(npw_k,nsym2,pinv,gprimd,symrec2t,gcurr,nbasek(ikpt),gbasek(:,:,ikpt),cnormk(:,ikpt))
1293 
1294    ABI_FREE(gcurr)
1295  end do
1296  !
1297  ! === Reduce info over k-points ===
1298  ! * Here symrec2t=TRANSPOSE(symrel2) for historical reasons, see note above
1299  sizepw=2*mpw
1300  ABI_MALLOC(gbase,(3,sizepw))
1301  ABI_MALLOC(cnorm,(sizepw))
1302  nbase=0                    ! # of irred G found.
1303 
1304  call merge_kgirr(nsym2,pinv,nkpt,mpw,sizepw,symrec2t,nbasek,cnormk,gbasek,nbase,gbase,cnorm,ierr)
1305  if (ierr/=0) then
1306    MSG_ERROR('merge_kgirr returned a non-zero status error')
1307  end if
1308 
1309  ABI_FREE(nbasek)
1310  ABI_FREE(cnormk)
1311  ABI_FREE(gbasek)
1312  !
1313  !=== Reorder base G-vectors in order of increasing module ===
1314  !
1315  !Generate all shells of G-vectors: star of a g==set of all symetrics of this g
1316  ABI_MALLOC(gshell,(3,2*nsym2))
1317  ABI_MALLOC(shlim,(nbase))
1318 
1319  ABI_MALLOC(gbig,(3,sizepw))
1320 !
1321 !TODO
1322 #if 0
1323 !* Here symrec2t=TRANSPOSE(symrel2) for historical reasons, see note above
1324  call getfullg(nbase,nsym2,pinv,sizepw,gbase,symrec2t,cnorm,maxpw,gbig,shlim,ierr)
1325  if (ierr/0) RETURN
1326 
1327 #else
1328  ABI_MALLOC(insort,(nbase))
1329  ABI_MALLOC(nshell,(nbase))
1330  do in=1,nbase
1331    insort(in)=in
1332  end do
1333  call sort_dp(nbase,cnorm,insort,tol14)
1334 !
1335 !Loop over all different modules of g''s (=shells):
1336  maxpw=0
1337  do in=1,nbase
1338    nshell(in)=0
1339    gcur(:)=gbase(:,insort(in))
1340 
1341    do is=1,nsym2 !  Loop over all symetries:
1342      do iinv=pinv,1,2
1343        geq(:)=iinv*(symrel2(1,:,is)*gcur(1)+symrel2(2,:,is)*gcur(2)+symrel2(3,:,is)*gcur(3))
1344 
1345        found=.FALSE.; ish=1
1346        do while ((.not.found) .and. (ish<=nshell(in))) ! Search for symetric of g and eventually add it:
1347          found=ALL(geq(:)==gshell(:,ish))
1348          ish=ish+1
1349        end do
1350        if (.not.found) then
1351          nshell(in)=nshell(in)+1
1352          gshell(:,nshell(in))=geq(:)
1353        end if
1354      end do
1355    end do
1356 
1357    if ((maxpw+nshell(in)) > sizepw) then
1358      ! We need to increase the size of the gbase, gbig and cnorm arrays while still keeping their content.
1359      ! This is done using two temporary arrays gtmp and ctmp
1360      MSG_WARNING("Had to reallocate gbase, gbig, cnorm")
1361      ABI_MALLOC(ctmp,(sizepw))
1362      ABI_MALLOC(gtmp,(3,sizepw))
1363      sizeold=sizepw
1364      sizepw=maxpw+nshell(in)
1365 
1366      ctmp(:)=cnorm(:)
1367      gtmp(:,:)=gbase(:,:)
1368 
1369      ABI_FREE(cnorm)
1370      ABI_MALLOC(cnorm,(sizepw))
1371      cnorm(1:sizeold)=ctmp(1:sizeold)
1372      cnorm(sizeold+1:sizepw)=zero
1373      ABI_FREE(ctmp)
1374 
1375 !    MG why this? gbase should not be changed!
1376      ABI_FREE(gbase)
1377      ABI_MALLOC(gbase,(3,sizepw))
1378      gbase(:,:sizeold)=gtmp(:,:sizeold)
1379      gbase(:,sizeold+1:sizepw)=0
1380      gtmp(:,:)=gbig(:,:)
1381 
1382      ABI_FREE(gbig)
1383      ABI_MALLOC(gbig,(3,sizepw))
1384      gbig(:,:sizeold)=gtmp(:,:sizeold)
1385      gbig(:,sizeold+1:sizepw)=0
1386      ABI_FREE(gtmp)
1387    end if
1388    !
1389    ! Store this shell of g''s in a big array of g (gbig):
1390    do ig=1,nshell(in)
1391      gbig(:,ig+maxpw)=gshell(:,ig)
1392    end do
1393    maxpw=maxpw+nshell(in)
1394 
1395  end do ! End loop over shells
1396  !
1397  ! * Compute shell limits
1398  ilim=0
1399  do in=1,nbase
1400    ilim=ilim+nshell(in)
1401    shlim(in)=ilim
1402  end do
1403 
1404  if (PRESENT(shlim_p)) then ! Return shlim_p
1405   ABI_MALLOC(shlim_p,(nbase))
1406   shlim_p = shlim
1407  end if
1408 
1409  ! Re-allocate gbig with correct sizes so that caller can inquire the size
1410  ABI_MALLOC(gtmp,(3,ilim))
1411  gtmp = gbig(:,1:ilim)
1412  ABI_FREE(gbig)
1413  ABI_MALLOC(gbig,(3,ilim))
1414  gbig=gtmp
1415  ABI_FREE(gtmp)
1416 
1417  if (prtvol>10) then ! Print out shell limits
1418    write(msg,'(3a)')&
1419 &    ' Shells found:',ch10,&
1420 &    ' number of shell    number of G vectors      cut-off energy [Ha} '
1421    call wrtout(std_out,msg,'COLL')
1422    do in=1,nbase
1423      write(msg,'(12x,i4,17x,i6,12x,f8.3)')in,shlim(in),2*pi**2*cnorm(in)
1424      call wrtout(std_out,msg,'COLL')
1425    end do
1426    call wrtout(std_out,ch10,'COLL')
1427  end if
1428 
1429  ABI_FREE(gshell)
1430  ABI_FREE(insort)
1431  ABI_FREE(nshell)
1432 #endif
1433 
1434  call destroy_mpi_enreg(MPI_enreg_seq)
1435  ABI_FREE(gbase)
1436  ABI_FREE(shlim)
1437  ABI_FREE(cnorm)
1438  ABI_FREE(npwarr)
1439 
1440 end subroutine merge_and_sort_kg

m_gsphere/merge_kgirr [ Functions ]

[ Top ] [ m_gsphere ] [ Functions ]

NAME

 merge_kgirr

FUNCTION

  Given a list of irreducible reciprocal vectors associated to different k-centered spheres,
  this subroutine finds the minimal set of G vectors needed to reconstruct the union of the spheres
  through symmetry operations.

INPUTS

 nsym=number of symmetry operations
 pinv=-1 if time-reversal can be used, 0 otherwise
 nkpt=number of k-points for k-centered spheres
 mpw=Max number of G vectors for each k-point
 sizepw=Max expected number of G vectors found
 symrec(3,3,nsym)=symmetries in reciprocal space given in reduced coordinates
 nbasek(nkpt)=number of irred G for each k-point
 cnormk(mpw,nkpt)=the norm of each k-centered G (only 1:nbase(ik)) is used
 gbasek(3,mpw,nkpt)

OUTPUT

 nbase=number of irreducible G needed to reconstruct the initial set of spheres
 gbase(3,sizepw)=irreducible G found in reciprocal coordinates
 cnorm(sizepw)=Norm of each irred G vector
 ierr= Exit status, if /=0 the number of G vectors found exceeds sizepw

PARENTS

      m_gsphere

CHILDREN

      gsph_free,gsph_init

SOURCE

1740 subroutine merge_kgirr(nsym,pinv,nkpt,mpw,sizepw,symrec,nbasek,cnormk,gbasek,nbase,gbase,cnorm,ierr)
1741 
1742 
1743 !This section has been created automatically by the script Abilint (TD).
1744 !Do not modify the following lines by hand.
1745 #undef ABI_FUNC
1746 #define ABI_FUNC 'merge_kgirr'
1747 !End of the abilint section
1748 
1749  implicit none
1750 
1751 !Arguments ------------------------------------
1752 !scalars
1753  integer,intent(in) :: mpw,nkpt,nsym,pinv,sizepw
1754  integer,intent(out) :: ierr,nbase
1755 !arrays
1756  integer,intent(in) :: gbasek(3,mpw,nkpt),nbasek(nkpt),symrec(3,3,nsym)
1757  integer,intent(inout) :: gbase(3,sizepw) !vz_i
1758  real(dp),intent(in) :: cnormk(mpw,nkpt)
1759  real(dp),intent(inout) :: cnorm(sizepw) !vz_i
1760 
1761 !Local variables-------------------------------
1762 !scalars
1763  integer :: ikpt,inb,irgk,isym
1764  real(dp) :: eps,norm
1765  logical :: found
1766  character(len=500) :: msg
1767 !arrays
1768  integer :: gbas(3),gcur(3),geq(3)
1769 
1770 ! *************************************************************************
1771 
1772  DBG_ENTER("COLL")
1773 
1774  if (pinv/=1.and.pinv/=-1) then
1775    write(msg,'(a,i6)')' The argument pinv should be -1 or 1, however, pinv =',pinv
1776    MSG_BUG(msg)
1777  end if
1778  !
1779  ! === Start with zero number of G found ===
1780  nbase=0 ; ierr=0
1781  do ikpt=1,nkpt
1782    do irgk=1,nbasek(ikpt)
1783      gcur(:)=gbasek(:,irgk,ikpt)
1784      norm=cnormk(irgk,ikpt) ; eps=tol8*norm
1785      found=.FALSE. ; inb=1
1786      do while ((.not.found).and.(inb<=nbase))
1787        if (ABS(norm-cnorm(inb))<=eps) then
1788          gbas(:)=gbase(:,inb)
1789          isym=1
1790          do while ((.not.found).and.(isym<=nsym))
1791            geq(:)=MATMUL(symrec(:,:,isym),gcur)
1792            found=ALL(geq(:)==gbas(:))
1793            if (pinv==-1) found= (found.or.ALL(geq(:)==-gbas(:)) ) ! For time-reversal
1794            isym=isym+1
1795          end do
1796        end if
1797        inb=inb+1
1798      end do
1799      if (.not.found) then
1800        ! === Add to the list ===
1801        nbase=nbase+1
1802        if (nbase>sizepw) then
1803          write(msg,'(2(a,i5),a)')&
1804 &         ' nbase (',nbase,') became greater than sizepw = ',sizepw,' returning ierr=1 '
1805          MSG_WARNING(msg)
1806          ierr=1; RETURN
1807        end if
1808        cnorm(nbase)=cnormk(irgk,ikpt)
1809        gbase(:,nbase)=gcur(:)
1810      end if
1811    end do
1812  end do
1813 
1814  DBG_EXIT("COLL")
1815 
1816 end subroutine merge_kgirr

m_gsphere/print_gsphere [ Functions ]

[ Top ] [ m_gsphere ] [ Functions ]

NAME

 print_gsphere

FUNCTION

  Print the content of a gvectors data type

INPUTS

  Gsph<gsphere_t>=Info on the G-sphere
  unit=the unit number for output
  prtvol = verbosity level
  mode_paral =either "COLL" or "PERS"

OUTPUT

  Only writing.

PARENTS

      cchi0q0,gwls_hamiltonian,setup_bse,setup_bse_interp

CHILDREN

      gsph_free,gsph_init

SOURCE

727 subroutine print_gsphere(Gsph,unit,prtvol,mode_paral)
728 
729 
730 !This section has been created automatically by the script Abilint (TD).
731 !Do not modify the following lines by hand.
732 #undef ABI_FUNC
733 #define ABI_FUNC 'print_gsphere'
734 !End of the abilint section
735 
736  implicit none
737 
738 !Arguments ------------------------------------
739 !scalars
740  integer,intent(in),optional :: prtvol,unit
741  character(len=4),intent(in),optional :: mode_paral
742  type(gsphere_t),intent(in) :: Gsph
743 
744 !Local variables-------------------------------
745 !scalars
746  integer :: ish,nsc,my_unt,my_prtvol
747  real(dp) :: fact,kin
748  character(len=4) :: my_mode
749  character(len=500) :: msg
750 
751 ! *************************************************************************
752 
753  my_unt    =std_out; if (PRESENT(unit      )) my_unt    =unit
754  my_prtvol=0       ; if (PRESENT(prtvol    )) my_prtvol=prtvol
755  my_mode   ='COLL' ; if (PRESENT(mode_paral)) my_mode   =mode_paral
756 
757  write(msg,'(3a,2(a,i8,a))')ch10,&
758 & ' ==== Info on the G-sphere ==== ',ch10,&
759 & '  Number of G vectors ... ',Gsph%ng,ch10,&
760 & '  Number of shells ...... ',Gsph%nsh,ch10
761  call wrtout(my_unt,msg,my_mode)
762 
763  SELECT CASE (Gsph%timrev)
764  CASE (1)
765    call wrtout(my_unt,' Time reversal symmetry cannot be used',my_mode)
766  CASE (2)
767    call wrtout(my_unt,' Time reversal symmetry is used',my_mode)
768  CASE DEFAULT
769    MSG_BUG("Wrong timrev")
770  END SELECT
771 
772  if (my_prtvol/=0) then
773    fact=half*two_pi**2
774    write(msg,'(a)')
775    call wrtout(my_unt,' Shell   Tot no. of Gs   Cutoff [Ha]',my_mode)
776    do ish=1,Gsph%nsh
777      nsc=Gsph%shlim(ish+1)-1
778      kin=half*Gsph%shlen(ish)**2
779      write(msg,'(2x,i4,10x,i6,5x,f8.3)')ish,nsc,kin
780      call wrtout(my_unt,msg,'COLL')
781    end do
782    call wrtout(my_unt,ch10,my_mode)
783  end if
784 
785 end subroutine print_gsphere

m_gsphere/prune_g1mg2 [ Functions ]

[ Top ] [ m_gsphere ] [ Functions ]

NAME

 prune_g1mg2

FUNCTION

 Given a list of G-vectors, evalute any possible difference G1-G2
 remove duplicated differences and report the list of inequivalent G-vectors.

INPUTS

  npw=Number of plane waves
  gvec(3,npw)= the reciprocal lattice vectors of the PW

OUTPUT

  ngdiff=Number of inequivalent differences G1-G2
  g1mg2(3,ngdiff)=The set of inequivalent G1-G2 vectors.

TODO

  Loop by shells to have better scaling.

PARENTS

CHILDREN

      gsph_free,gsph_init

SOURCE

1094 subroutine prune_g1mg2(npw,gvec,ngdiff,g1mg2)
1095 
1096 
1097 !This section has been created automatically by the script Abilint (TD).
1098 !Do not modify the following lines by hand.
1099 #undef ABI_FUNC
1100 #define ABI_FUNC 'prune_g1mg2'
1101 !End of the abilint section
1102 
1103  implicit none
1104 
1105 !Arguments ------------------------------------
1106 !scalars
1107  integer,intent(in) :: npw
1108  integer,intent(out) :: ngdiff
1109 !arrays
1110  integer,intent(in) :: gvec(3,npw)
1111  integer,allocatable,intent(out) :: g1mg2(:,:)
1112 
1113 !Local variables ------------------------------
1114 !scalars
1115  integer :: ii,ig1,ig2
1116  logical :: found
1117  character(len=500) :: msg
1118 !arrays
1119  integer :: gdiff(3)
1120  integer,allocatable :: g1mg2_tmp(:,:)
1121 
1122 !************************************************************************
1123 
1124  ABI_MALLOC(g1mg2_tmp,(3,9*npw))
1125 
1126  ngdiff=0
1127  do ig2=1,npw
1128    do ig1=1,npw
1129     gdiff = gvec(:,ig1) - gvec(:,ig2)
1130 
1131     found=.FALSE. ; ii=0
1132     do while (.not.found .and. ii<ngdiff)
1133       ii = ii+1
1134       found = ALL(gdiff == g1mg2_tmp(:,ii) )
1135     end do
1136 
1137     if (.not.found) then
1138       ngdiff = ngdiff + 1
1139       if (ngdiff > 9*npw) GOTO 100
1140       g1mg2_tmp(:,ngdiff) = gdiff
1141     end if
1142 
1143    end do
1144  end do
1145 
1146  ! * Save results
1147  ABI_MALLOC(g1mg2,(3,ngdiff))
1148  g1mg2 = g1mg2_tmp(:,1:ngdiff)
1149  ABI_FREE(g1mg2_tmp)
1150 
1151  RETURN
1152 
1153 100 continue
1154  write(msg,'(2(a,i6))')' ngdiff = ',ngdiff,' > 9*npw = ',9*npw
1155  MSG_BUG(msg)
1156 
1157 end subroutine prune_g1mg2

m_gsphere/setup_G_rotation [ Functions ]

[ Top ] [ m_gsphere ] [ Functions ]

NAME

 setup_G_rotation

FUNCTION

 Set up tables indicating rotation of G-vectors.

INPUTS

 nsym=Number of symmetry operations.
 symrec(3,3,nsym)=Symmetry operations in reciprocal space.
 timrev=2 if time reversal can be used, 1 otherwise.
 npw=Number of planewaves in the sphere.
 gvec(3,npw)=Coordinates of plane waves, supposed to be ordered in increasing modulus
 g2sh(npw)=For each G, it gives the index of the shell to which it belongs.
 nsh=Number of shells
 shlim(nsh+1)=Index of the first G vector in each shell, =npw+1 for nsh+1

OUTPUT

  grottb  (npw,2,nsym)= grottb(G,I,S) is the index of (SI) G in the array gvec.
  grottbm1(npw,2,nsym)= index of IS^{-1} G.

 NOTES:
  I is either the identity or the inversion (time reversal in reciprocal space).
  S is one of the symmetry operation in reciprocal space belonging to the Space group.

PARENTS

      m_gsphere

CHILDREN

      gsph_free,gsph_init

SOURCE

211 subroutine setup_G_rotation(nsym,symrec,timrev,npw,gvec,g2sh,nsh,shlim,grottb,grottbm1)
212 
213 
214 !This section has been created automatically by the script Abilint (TD).
215 !Do not modify the following lines by hand.
216 #undef ABI_FUNC
217 #define ABI_FUNC 'setup_G_rotation'
218 !End of the abilint section
219 
220  implicit none
221 
222 !Arguments ------------------------------------
223 !scalars
224  integer,intent(in) :: npw,nsh,nsym,timrev
225 !arrays
226  integer,intent(in) :: g2sh(npw),gvec(3,npw),shlim(nsh+1),symrec(3,3,nsym)
227  integer,intent(inout) :: grottb  (npw,timrev,nsym)
228  integer,intent(inout) :: grottbm1(npw,timrev,nsym)
229 
230 !Local variables ------------------------------
231 !scalars
232  integer :: ee,ig1,ig2,ish1,isym,itim,ss
233  logical :: found
234  character(len=500) :: msg
235 !arrays
236  integer :: gbase(3),grot(3)
237 
238 !************************************************************************
239  !
240  ! === Set up G-rotation table ===
241  do ig1=1,npw
242    ish1=g2sh(ig1) ; ss=shlim(ish1) ; ee=shlim(ish1+1)-1
243    gbase(:)=gvec(:,ig1)
244 
245    do itim=1,timrev
246      do isym=1,nsym
247        grot=(3-2*itim)*MATMUL(symrec(:,:,isym),gbase)
248        found=.FALSE.
249        ! * Loop on the shell of ig1 to speed up the search.
250        do ig2=ss,ee
251          if (ALL(ABS(grot(:)-gvec(:,ig2))==0)) then
252            found=.TRUE.
253            grottb  (ig1,itim,isym)=ig2
254            grottbm1(ig2,itim,isym)=ig1
255          end if
256        end do
257        if (.not.found) then
258          write(msg,'(3a,i5,a,i5,1x,2(3i5,a),a,i3,a,i3)')&
259 &         'G-shell not closed',ch10,&
260 &         '  Initial G vector ',ig1,'/',npw,gbase(:),' Rotated G vector ',grot(:),ch10,&
261 &         '  Through sym ',isym,' and itim ',itim
262          MSG_ERROR(msg)
263        end if
264      end do
265    end do
266 
267  end do !ig1
268 
269 end subroutine setup_G_rotation

m_gsphere/symg [ Functions ]

[ Top ] [ m_gsphere ] [ Functions ]

NAME

 symg

FUNCTION

 Treat symmetries applied to the G vectors, in view of the application
 to symmetrization of the dielectric matrix.
 Generate a list of time-reversed G vectors, as well as a list
 of spatially-symmetric G vectors.

INPUTS

 kg_diel(3,npwdiel)=reduced planewave coordinates for the dielectric matrix.
 npwdiel=number of planewaves for the dielectric matrix
 nsym=number of symmetry
 symrel(3,3,nsym)=symmetry matrices in real space (integers)
 tnons(3,nsym)=reduced nonsymmorphic translations
 (symrel and tnons are in terms of real space primitive translations)

OUTPUT

 phdiel(2,npwdiel,nsym)=phase associated with point symmetries applied to G
 sym_g(npwdiel,nsym)=index list of symmetric G vectors
 (could save a bit of space by suppressing isym=1, since the
 corresponding symmetry is the identity)
 tmrev_g(npwdiel)=index list of inverted G vectors (time-reversed)

PARENTS

      suscep_stat

CHILDREN

SOURCE

2595 subroutine symg(kg_diel,npwdiel,nsym,phdiel,sym_g,symrel,tmrev_g,tnons)
2596 
2597 
2598 !This section has been created automatically by the script Abilint (TD).
2599 !Do not modify the following lines by hand.
2600 #undef ABI_FUNC
2601 #define ABI_FUNC 'symg'
2602 !End of the abilint section
2603 
2604  implicit none
2605 
2606 !Arguments ------------------------------------
2607 !scalars
2608  integer,intent(in) :: npwdiel,nsym
2609 !arrays
2610  integer,intent(in) :: kg_diel(3,npwdiel),symrel(3,3,nsym)
2611  integer,intent(out) :: sym_g(npwdiel,nsym),tmrev_g(npwdiel)
2612  real(dp),intent(in) :: tnons(3,nsym)
2613  real(dp),intent(out) :: phdiel(2,npwdiel,nsym)
2614 
2615 !Local variables-------------------------------
2616 !scalars
2617  integer :: g1,g2,g3,ipw,isym,j1,j2,j3,m1m,m1p,m2m,m2p,m3m,m3p,symmg,trevg
2618  real(dp) :: arg,tau1,tau2,tau3
2619  !character(len=500) :: message
2620 !arrays
2621  integer,allocatable :: grid(:,:,:)
2622 
2623 ! *************************************************************************
2624 
2625 !Determines maximal bounds of the zone spanned by the planewaves
2626  m1m=0 ; m2m=0 ; m3m=0 ; m1p=0 ; m2p=0 ; m3p=0
2627  do ipw=1,npwdiel
2628    g1=kg_diel(1,ipw)
2629    g2=kg_diel(2,ipw)
2630    g3=kg_diel(3,ipw)
2631    if(g1<m1m)m1m=g1 ; if(g1>m1p)m1p=g1
2632    if(g2<m2m)m2m=g2 ; if(g2>m2p)m2p=g2
2633    if(g3<m3m)m3m=g3 ; if(g3>m3p)m3p=g3
2634  end do
2635 
2636 !Set up grid, that associate to each point the index of the
2637 !corresponding planewave, if there is one
2638  ABI_ALLOCATE(grid,(m1m:m1p,m2m:m2p,m3m:m3p))
2639  grid(:,:,:)=0
2640  do ipw=1,npwdiel
2641    g1=kg_diel(1,ipw)
2642    g2=kg_diel(2,ipw)
2643    g3=kg_diel(3,ipw)
2644    grid(g1,g2,g3)=ipw
2645  end do
2646 
2647 !Set up tmrev_g and sym_g arrays
2648  do ipw=1,npwdiel
2649    g1=kg_diel(1,ipw)
2650    g2=kg_diel(2,ipw)
2651    g3=kg_diel(3,ipw)
2652 
2653 !  Treat first time-reversal symmetry
2654    trevg=grid(-g1,-g2,-g3)
2655    if(trevg==0)then
2656      MSG_BUG('Do not find the time-reversed symmetric of a G-vector.')
2657    end if
2658    tmrev_g(ipw)=trevg
2659 
2660 !  Treat now spatial symmetries
2661    do isym=1,nsym
2662 
2663 !    Get rotated G vector Gj for each symmetry element
2664 !    -- here we use the TRANSPOSE of symrel; assuming symrel expresses
2665 !    the rotation in real space, the transpose is then appropriate
2666 !    for G space symmetrization (according to Doug : see routine irrzg.f)
2667      j1=symrel(1,1,isym)*g1+&
2668 &     symrel(2,1,isym)*g2+symrel(3,1,isym)*g3
2669      j2=symrel(1,2,isym)*g1+&
2670 &     symrel(2,2,isym)*g2+symrel(3,2,isym)*g3
2671      j3=symrel(1,3,isym)*g1+&
2672 &     symrel(2,3,isym)*g2+symrel(3,3,isym)*g3
2673      symmg=grid(j1,j2,j3)
2674      if(symmg==0)then
2675        MSG_BUG('Do not find the spatially symmetric of a G-vector.')
2676      end if
2677      sym_g(ipw,isym)=symmg
2678 
2679 !    Get associated phase
2680      tau1=tnons(1,isym)
2681      tau2=tnons(2,isym)
2682      tau3=tnons(3,isym)
2683      if (abs(tau1)>tol12.or.abs(tau2)>tol12.or.abs(tau3)>tol12) then
2684 !      compute exp(-2*Pi*I*G dot tau) using original G
2685        arg=two_pi*(dble(g1)*tau1+dble(g2)*tau2+dble(g3)*tau3)
2686        phdiel(1,ipw,isym)=cos(arg)
2687        phdiel(2,ipw,isym)=-sin(arg)
2688      else
2689        phdiel(1,ipw,isym)=1._dp
2690        phdiel(2,ipw,isym)=0._dp
2691      end if
2692 
2693    end do
2694  end do
2695 
2696  ABI_DEALLOCATE(grid)
2697 
2698 end subroutine symg

m_gsphere/table_gbig2kg [ Functions ]

[ Top ] [ m_gsphere ] [ Functions ]

NAME

  table_gbig2kg

FUNCTION

  Associate the kg_k set of g-vectors with the big array of gbig
  The array gbig(3,maxpw) contains all g-vectors used for all k-points, in order of
  increasing shells. For a each k-point, the wave-functions are defined only on a particular set
  of g-vectors kg_k (included in gbig). This set is defined by array gamma2k:
  The array gamma2k(ig=1,maxpw) translates the index of the gbig (from 1 to maxpw) into the corresponding
  index in array kg_k. If gbig(ig) does not exist in kg_k, gamma2k(ig) contains npw_k+1.

INPUTS

  npw_k=Number of planewaves in the k-centered basis set
  kg_k(3,npw_k)=The k-centered basis set
  maxpw=Number of G in gbig
  gbig(3,maxpw)=The union of the G-spheres at different k-points.

OUTPUT

  ierr=Status error. It gives the number of G of kg_k not contained in gbig.
  gamma2k(maxpw)=Mapping gbig -> kg_k

PARENTS

      m_io_kss

CHILDREN

      gsph_free,gsph_init

SOURCE

2332 pure subroutine table_gbig2kg(npw_k,kg_k,maxpw,gbig,gamma2k,ierr)
2333 
2334 
2335 !This section has been created automatically by the script Abilint (TD).
2336 !Do not modify the following lines by hand.
2337 #undef ABI_FUNC
2338 #define ABI_FUNC 'table_gbig2kg'
2339 !End of the abilint section
2340 
2341  implicit none
2342 
2343 !Arguments ------------------------------------
2344 !scalars
2345  integer,intent(in) :: npw_k,maxpw
2346  integer,intent(out) :: ierr
2347 !arrays
2348  integer,intent(in) :: kg_k(3,npw_k)
2349  integer,intent(in) :: gbig(3,maxpw)
2350  integer,intent(out) :: gamma2k(maxpw)
2351 
2352 !Local variables-------------------------------
2353 !scalars
2354  integer :: ig,igp
2355  logical :: found
2356 !arrays
2357  integer :: gcur(3)
2358 
2359 ! *********************************************************************
2360 
2361  ierr=0
2362  gamma2k(:)=npw_k+1  ! Initialize array gamma2k
2363 
2364  do ig=1,npw_k       ! Loop over g-vectors, for this k point.
2365    gcur(:)=kg_k(:,ig)
2366    igp=0; found=.FALSE.
2367    do while ((.not.found) .and. igp<maxpw) ! Search selected vector in array gbig: TODO this part can be optimized
2368      igp=igp+1
2369      found=ALL(gcur(:)==gbig(:,igp))
2370    end do
2371    if (found) then ! Store it if found:
2372      gamma2k(igp)=ig
2373    else
2374      ierr=ierr+1
2375    end if
2376  end do
2377 
2378 end subroutine table_gbig2kg