TABLE OF CONTENTS


ABINIT/m_vcoul [ Modules ]

[ Top ] [ Modules ]

NAME

  m_vcoul

FUNCTION

  This module contains the definition of the vcoul_t as well
  as procedures to calculate the Coulomb interaction in reciprocal space
  taking into account a possible cutoff in real space.
  Procedures to deal with the singularity for q --> 0 are also provided.

COPYRIGHT

 Copyright (C) 1999-2022 ABINIT group (MG, FB)
 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 .

SOURCE

20 #if defined HAVE_CONFIG_H
21 #include "config.h"
22 #endif
23 
24 #include "abi_common.h"
25 
26 module m_vcoul
27 
28  use defs_basis
29  use m_abicore
30  use m_errors
31  use m_xmpi
32  use m_splines
33  use m_sort
34 
35  use m_fstrings,        only : sjoin, itoa
36  use m_special_funcs,   only : abi_derf
37  use m_bessel,          only : calck0
38  use m_io_tools,        only : open_file
39  use m_gwdefs,          only : GW_TOLQ0
40  use m_numeric_tools,   only : arth, geop, imin_loc, llsfit_svd, l2norm, OPERATOR(.x.), quadrature, isdiagmat
41  use m_hide_lapack,     only : matrginv
42  use m_geometry,        only : normv, metric
43  use m_qplusg,          only : cmod_qpg
44  use m_crystal,         only : crystal_t
45  use m_bz_mesh,         only : kmesh_t
46  use m_gsphere,         only : gsphere_t
47  use m_fftcore,         only : get_kg
48  use m_dtfil,           only : isfile
49 
50  ! Cut-off methods modules
51  use m_cutoff_sphere,   only : cutoff_sphere
52  use m_cutoff_slab,     only : cutoff_slab
53  use m_cutoff_cylinder, only : cutoff_cylinder
54 
55  implicit none
56 
57  public :: gw_icutcoul_to_mode
58  public :: carrier_isz
59 
60  private

m_vcoul/adapt_nmc [ Functions ]

[ Top ] [ m_vcoul ] [ Functions ]

NAME

 adapt_nmc

FUNCTION

 Empirical law to decrease the Monte Carlo sampling
 for large q+G, for which the accuracy is not an issue

m_vcoul/beigi_cylinder_limit [ Functions ]

[ Top ] [ m_vcoul ] [ Functions ]

NAME

 beigi_cylinder_limit

FUNCTION

SOURCE

1456 subroutine beigi_cylinder_limit(opt_cylinder, cryst, nqibz, nkbz, rcut, hcyl, boxcenter, pdir, i_sz)
1457 
1458 !Arguments ------------------------------------
1459  integer,intent(in) :: opt_cylinder, nqibz, nkbz, pdir(3)
1460  type(crystal_t),intent(in) :: cryst
1461  real(dp),intent(in) :: rcut, hcyl, boxcenter(3)
1462  real(dp),intent(out) :: i_sz
1463 
1464 !Local variables-------------------------------
1465  integer :: ii, iq, npar, npt, gamma_pt(3,1)
1466  real(dp) :: step, bz_plane, dx, integ, q0_vol, q0_volsph, b1(3),b2(3),b3(3)
1467  real(dp),allocatable :: cov(:,:),par(:),qfit(:,:),sigma(:),var(:), vcfit(:,:),xx(:),yy(:)
1468 
1469 ! *************************************************************************
1470 
1471  b1 = two_pi * cryst%gprimd(:,1); b2 = two_pi * cryst%gprimd(:,2); b3 = two_pi * cryst%gprimd(:,3)
1472 
1473  npt =100
1474  npar=8; gamma_pt = RESHAPE(([0, 0, 0]), [3, 1])
1475  ABI_MALLOC(qfit, (3, npt))
1476  ABI_MALLOC(vcfit, (1, npt))
1477  if (nqibz == 1) then
1478    ABI_ERROR("nqibz == 1 not supported when Beigi's method is used")
1479  endif
1480  qfit(:,:)=zero
1481  step=half/(npt * (nqibz-1))              ; qfit(3,:)=arth(tol6,step,npt)
1482  !step=(half/(nqibz-1)/tol6)**(one/npt) ; qfit(3,:)=geop(tol6,step,npt)
1483 
1484  do iq=1,npt
1485    call cutoff_cylinder(qfit(:,iq),1,gamma_pt,rcut,hcyl,pdir,boxcenter,&
1486                         Cryst%rprimd,vcfit(:,iq),opt_cylinder, xmpi_comm_self)
1487  end do
1488 
1489  ABI_MALLOC(xx, (npt))
1490  ABI_MALLOC(yy, (npt))
1491  ABI_MALLOC(sigma, (npt))
1492  ABI_MALLOC(par, (npar))
1493  ABI_MALLOC(var, (npar))
1494  ABI_MALLOC(cov, (npar, npar))
1495 
1496  do ii=1,npt
1497    xx(ii) = normv(qfit(:,ii), cryst%gmet, 'G')
1498  end do
1499  ABI_FREE(qfit)
1500  sigma=one ; yy(:)=vcfit(1,:)
1501  ABI_FREE(vcfit)
1502  !call llsfit_svd(xx,yy,sigma,npar,K0fit,chisq,par,var,cov,info)
1503  !do ii=1,npt
1504  ! write(99,*)xx(ii),yy(ii),DOT_PRODUCT(par,K0fit(xx(ii),npar))
1505  !end do
1506  bz_plane=l2norm(b1.x.b2)
1507  !integ=K0fit_int(xx(npt),par,npar)
1508  !write(std_out,*)' SVD fit : chi-square',chisq
1509  !write(std_out,*)' fit-parameters : ',par
1510  !write(std_out,*)' variance ',var
1511  !write(std_out,*)' bz_plane ',bz_plane
1512  !write(std_out,*)' SCD integ ',integ
1513  ! Here Im assuming homogeneous mesh
1514  dx=(xx(2)-xx(1))
1515  integ=yy(2)*dx*3.0/2.0
1516  do ii=3,npt-2
1517    integ=integ+yy(ii)*dx
1518  end do
1519  integ=integ+yy(npt-1)*dx*3.0/2.0
1520  !write(std_out,*)' simple integral',integ
1521  q0_volsph = (two_pi)**3 / (nkbz * cryst%ucvol)
1522  q0_vol=bz_plane*two*xx(npt)
1523  !write(std_out,*)' q0 sphere : ',q0_volsph,' q0_vol cyl ',q0_vol
1524  i_sz = bz_plane * two * integ / q0_vol
1525  !write(std_out,*)' spherical approximation ',four_pi*7.44*q0_volsph**(-two_thirds)
1526  !write(std_out,*)' Cylindrical cutoff value ',i_sz
1527  !i_sz=four_pi*7.44*q0_vol**(-two_thirds)
1528 
1529  ABI_FREE(xx)
1530  ABI_FREE(yy)
1531  ABI_FREE(sigma)
1532  ABI_FREE(par)
1533  ABI_FREE(var)
1534  ABI_FREE(cov)
1535 
1536 end subroutine beigi_cylinder_limit

m_vcoul/beigi_surface_limit [ Functions ]

[ Top ] [ m_vcoul ] [ Functions ]

NAME

 beigi_surface_limit

FUNCTION

SOURCE

1549 subroutine beigi_surface_limit(opt_slab, cryst, nqibz, nkbz, rcut, alpha, boxcenter, pdir, i_sz)
1550 
1551 !Arguments ------------------------------------
1552  integer,intent(in) :: opt_slab, nqibz, nkbz, pdir(3)
1553  type(crystal_t),intent(in) :: cryst
1554  real(dp),intent(in) :: rcut, alpha(3), boxcenter(3)
1555  real(dp),intent(out) :: i_sz
1556 
1557 !Local variables-------------------------------
1558  integer :: ii, npt, gamma_pt(3,1)
1559  real(dp) :: step, bz_plane, dx, integ, q0_vol, q0_volsph, b1(3),b2(3),b3(3)
1560  real(dp),allocatable :: qfit(:,:),sigma(:),vcfit(:,:),xx(:),yy(:), qcart(:,:)
1561 
1562 ! *************************************************************************
1563 
1564  b1 = two_pi * cryst%gprimd(:,1); b2 = two_pi * cryst%gprimd(:,2); b3 = two_pi * cryst%gprimd(:,3)
1565 
1566  gamma_pt=RESHAPE([0, 0, 0], [3, 1]) ! Gamma point
1567  npt=100 ! Number of points in 1D
1568  ABI_MALLOC(qfit, (3, npt))
1569  ABI_MALLOC(qcart, (3, npt))
1570  ABI_MALLOC(vcfit, (1, npt))
1571  if (nqibz == 1) then
1572    ABI_ERROR("nqibz == 1 not supported when Beigi's method is used")
1573  endif
1574  qfit(:,:)=zero
1575  qcart(:,:)=zero
1576  ! Size of the third vector
1577  bz_plane=l2norm(b3)
1578  q0_volsph=(two_pi)**3 / (nkbz * cryst%ucvol)
1579  ! radius that gives the same volume as q0_volsph
1580  ! Let's assume that c is perpendicular to the plane
1581  ! We also assume isotropic BZ around gamma
1582  step=sqrt((q0_volsph/bz_plane)/pi)/npt
1583 
1584  !step=half/(npt*(nqibz-1))
1585  ! Let's take qpoints along 1 line, the vcut does depend only on the norm
1586  qcart(1,:) = arth(tol6,step,npt)
1587 
1588  do ii=1,npt
1589    qfit(:,ii) = MATMUL(TRANSPOSE(Cryst%rprimd),qcart(:,ii)) / (2*pi)
1590    call cutoff_slab(qfit(:,ii), 1, gamma_pt, cryst%gprimd, rcut, &
1591                     boxcenter, pdir, alpha, vcfit(:,ii), opt_slab)
1592  end do
1593 
1594  ABI_MALLOC(xx, (npt))
1595  ABI_MALLOC(yy, (npt))
1596  ABI_MALLOC(sigma, (npt))
1597  do ii=1,npt
1598    !xx(ii)=qfit(1,:)
1599    xx(ii) = normv(qfit(:,ii), cryst%gmet, 'G')
1600  end do
1601  ABI_FREE(qfit)
1602  sigma=one
1603  yy(:)=vcfit(1,:)
1604  !yy(:)=one
1605  ABI_FREE(vcfit)
1606  dx=(xx(2)-xx(1))
1607  ! integ = \int dr r f(r)
1608  integ=xx(2)*yy(2)*dx*3.0/2.0
1609  do ii=3,npt-2
1610    integ=integ+xx(ii)*yy(ii)*dx
1611  end do
1612  integ=integ+xx(npt-1)*yy(npt-1)*dx*3.0/2.0
1613  !write(std_out,*)' simple integral',integ
1614  q0_vol=bz_plane*pi*xx(npt)**2
1615  !write(std_out,*)' q0 sphere : ',q0_volsph,' q0_vol cyl ',q0_vol
1616  i_sz=bz_plane*2*pi*integ/q0_vol
1617  !write(std_out,*)' spherical approximation ',four_pi*7.44*q0_volsph**(-two_thirds)
1618  !write(std_out,*)' Cylindrical cutoff value ',i_sz
1619  !i_sz=four_pi*7.44*q0_vol**(-two_thirds)
1620  ABI_FREE(xx)
1621  ABI_FREE(yy)
1622 
1623 end subroutine beigi_surface_limit

m_vcoul/carrier_isz [ Functions ]

[ Top ] [ m_vcoul ] [ Functions ]

NAME

  carrier_isz

FUNCTION

SOURCE

1634 real(dp) function carrier_isz(cryst, nqbz, qbz, rcut, comm) result(i_sz)
1635 
1636 !Arguments ------------------------------------
1637  type(crystal_t),intent(in) :: cryst
1638  integer,intent(in) :: nqbz, comm
1639  real(dp), intent(in) :: qbz(3, nqbz), rcut
1640 
1641 !Local variables-------------------------------
1642  integer :: iq_bz
1643  real(dp) :: qbz_norm, bz_geometry_factor, qbz_cart(3), b1(3), b2(3), b3(3)
1644 
1645 !************************************************************************
1646 
1647  b1 = two_pi * cryst%gprimd(:,1); b2 = two_pi * cryst%gprimd(:,2); b3 = two_pi * cryst%gprimd(:,3)
1648 
1649  bz_geometry_factor = zero
1650  do iq_bz=1,nqbz
1651    qbz_cart(:) = qbz(1,iq_bz)*b1(:) + qbz(2,iq_bz)*b2(:) + qbz(3,iq_bz)*b3(:)
1652    qbz_norm = SQRT(SUM(qbz_cart(:)**2))
1653    if (qbz_norm > TOLQ0) bz_geometry_factor = bz_geometry_factor - faux(qbz(:,iq_bz), rcut, b1, b2, b3)
1654  end do
1655 
1656  bz_geometry_factor = bz_geometry_factor + integratefaux(rcut, cryst%gprimd, cryst%ucvol, comm) * nqbz
1657  i_sz = four_pi * bz_geometry_factor  ! Final result stored here
1658 
1659 end function carrier_isz

m_vcoul/cylinder_setup [ Functions ]

[ Top ] [ m_vcoul ] [ Functions ]

NAME

  cylinder_setup

FUNCTION

SOURCE

594 subroutine cylinder_setup(cryst, vcutgeo, hcyl, pdir, opt_cylinder)
595 
596  type(crystal_t),intent(in) :: cryst
597  real(dp),intent(in) :: vcutgeo(3)
598  real(dp),intent(out) :: hcyl
599  integer,intent(out) :: pdir(3), opt_cylinder
600 
601 !Local variables-------------------------------
602  integer :: ii
603  real(dp),parameter :: tol999 = 999.0
604  real(dp) :: check
605 
606 ! *************************************************************************
607 
608  ABI_CHECK(count(abs(vcutgeo) > tol6) == 1, 'Wrong cutgeo for cylinder')
609 
610  ! Beigi's method is the default one, i.e infinite cylinder of radius rcut.
611  ! Use negative values to use Rozzi's method with finite cylinder of extent hcyl.
612  opt_cylinder = 1; hcyl = zero; pdir(:) = 0
613  do ii=1,3
614    check = vcutgeo(ii)
615    if (abs(check) > tol6) then
616      pdir(ii) = 1
617      if (check < zero) then
618        ! use Rozzi's method.
619        hcyl = ABS(check) *SQRT(SUM(cryst%rprimd(:,ii)**2))
620        opt_cylinder = 2
621        ! Check to enter the infinite Rozzi treatment
622        if(vcutgeo(3) <= -tol999) hcyl = tol12
623      end if
624    end if
625  end do
626 
627  ABI_CHECK((count(pdir == 1) == 1), 'Wrong pdir for cylinder')
628  if (pdir(3) /= 1) then
629    ABI_ERROR("The cylinder must be along the z-axis")
630  end if
631 
632 end subroutine cylinder_setup

m_vcoul/gw_icutcoul_to_mode [ Functions ]

[ Top ] [ m_vcoul ] [ Functions ]

NAME

 gw_icutcoul_to_mode

FUNCTION

 Convert gw_icutcoul_to_mode to mode string.

SOURCE

259 subroutine gw_icutcoul_to_mode(gw_icutcoul, mode)
260 
261 !Arguments ------------------------------------
262  integer,intent(in) :: gw_icutcoul
263  character(len=*),intent(out) :: mode
264 
265 ! *************************************************************************
266 
267  mode = 'NONE'
268  if (gw_icutcoul == 0) mode = 'SPHERE'
269  if (gw_icutcoul == 1) mode = 'CYLINDER'
270  if (gw_icutcoul == 2) mode = 'SLAB'
271  if (gw_icutcoul == 3) mode = 'CRYSTAL'
272  if (gw_icutcoul == 4) mode = 'ERF'
273  if (gw_icutcoul == 5) mode = 'ERFC'
274  if (gw_icutcoul == 6) mode = 'AUXILIARY_FUNCTION'
275  if (gw_icutcoul == 7) mode = 'AUX_GB'
276  if (gw_icutcoul == 14) mode = 'MINIBZ-ERF'
277  if (gw_icutcoul == 15) mode = 'MINIBZ-ERFC'
278  if (gw_icutcoul == 16) mode = 'MINIBZ'
279 
280 end subroutine gw_icutcoul_to_mode

m_vcoul/gygi_baldereschi_isz [ Functions ]

[ Top ] [ m_vcoul ] [ Functions ]

NAME

  gygi_baldereschi_isz

FUNCTION

SOURCE

1670 real(dp) function gygi_baldereschi_isz(cryst, nqbz, qbz, ecut, ng, gvec) result(i_sz)
1671 
1672 !Arguments ------------------------------------
1673  type(crystal_t),intent(in) :: cryst
1674  integer,intent(in) :: nqbz, ng
1675  real(dp), intent(in) :: qbz(3, nqbz), ecut
1676  integer,intent(in) :: gvec(3,ng)
1677 
1678 !Local variables-------------------------------
1679  integer :: iq_bz, ig
1680  real(dp) :: bz_geometry_factor, intfauxgb, alfa, qpg2, qpg(3)
1681 
1682 !************************************************************************
1683 
1684  ! the choice of alfa (the width of the gaussian) is somehow empirical
1685  alfa = 150.0 / ecut
1686 
1687  bz_geometry_factor=zero
1688  do iq_bz=1,nqbz
1689    do ig = 1,ng
1690      qpg(:) = qbz(:,iq_bz) + gvec(:,ig)
1691      qpg2 = normv(qpg, cryst%gmet, 'G')**2
1692      if (qpg2 > TOLQ0) bz_geometry_factor = bz_geometry_factor - EXP(-alfa*qpg2)/qpg2
1693    end do
1694  end do
1695 
1696  intfauxgb = cryst%ucvol/four_pi/SQRT(0.5*two_pi*alfa)
1697  bz_geometry_factor = bz_geometry_factor + intfauxgb * nqbz
1698 
1699  i_sz = four_pi*bz_geometry_factor
1700 
1701 end function gygi_baldereschi_isz

m_vcoul/integratefaux [ Functions ]

[ Top ] [ m_vcoul ] [ Functions ]

NAME

  integratefaux

FUNCTION

SOURCE

690 real(dp) function integratefaux(rcut, gprimd, ucvol, comm)
691 
692  real(dp),intent(in) :: rcut, gprimd(3,3), ucvol
693  integer,intent(in) :: comm
694 
695 !Local variables-------------------------------
696  integer,parameter :: nref = 3, nq = 50
697  integer :: ierr,iq,iqx1,iqy1,iqz1,iqx2,iqy2,iqz2,miniqy1,maxiqy1,nqhalf, nprocs, my_rank
698  real(dp) :: invnq,invnq3,qq,weightq,weightxy,weightxyz
699  real(dp) :: qq1(3),qq2(3),bb4sinpiqq_2(3,nq),sin2piqq(nq),bb4sinpiqq2_2(3,0:nq),sin2piqq2(3,0:nq)
700  real(dp) :: b1(3), b2(3), b3(3), bb(3)
701  real(dp) :: b1b1,b2b2,b3b3,b1b2,b2b3,b3b1
702 
703 ! *************************************************************************
704 
705  ! nq is the number of sampling points along each axis for the numerical integration
706  ! nref is the area where the mesh is refined
707 
708  integratefaux = zero
709  invnq = one/DBLE(nq)
710  invnq3 = invnq**3
711  nqhalf = nq/2
712 
713  b1 = two_pi * gprimd(:,1); b2 = two_pi * gprimd(:,2); b3 = two_pi * gprimd(:,3)
714  b1b1 = dot_product(b1, b1); b2b2 = dot_product(b2, b2); b3b3 = dot_product(b3, b3)
715  bb(1) = b1b1; bb(2) = b2b2; bb(3) = b3b3
716  b1b2 = dot_product(b1, b2); b2b3 = dot_product(b2, b3); b3b1 = dot_product(b3, b1)
717 
718  nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
719 
720  ! In order to speed up the calculation, precompute the sines
721  do iq=1,nq
722    qq=DBLE(iq)*invnq-half
723    bb4sinpiqq_2(:,iq)=bb(:)*four*SIN(pi*qq)**2 ; sin2piqq(iq)=SIN(two_pi*qq)
724  end do
725 
726  do iqx1=1,nq
727    if (modulo(iqx1, nprocs) /= my_rank) cycle ! MPI parallelism
728    qq1(1)=DBLE(iqx1)*invnq-half
729    ! Here take advantage of the q <=> -q symmetry:
730    ! arrange the sampling of qx, qy space to avoid duplicating calculations. Need weights to do this ...
731    !do iqy1=1,nq
732    miniqy1 = nqhalf + 1; maxiqy1 = nq
733    if (iqx1 >= nqhalf) miniqy1 = nqhalf
734    if (iqx1 > nqhalf .and. iqx1 < nq) maxiqy1 = nq - 1
735 
736    do iqy1=miniqy1,maxiqy1
737      qq1(2) = DBLE(iqy1)*invnq - half
738      ! By default, the factor of two is for the q <=> -q symmetry
739      weightq = invnq3*two
740      ! But not all qx qy lines have a symmetric one ...
741      if( (iqx1 == nqhalf .or. iqx1 == nq) .and. (iqy1 == nqhalf .or. iqy1 == nq)) weightq = weightq*half
742 
743      do iqz1=1,nq
744        qq1(3) = DBLE(iqz1)*invnq - half
745 
746        ! Refine the mesh for the point close to the origin
747        if( abs(iqx1-nqhalf) <= nref .and. abs(iqy1-nqhalf) <= nref .and. abs(iqz1-nqhalf) <= nref ) then
748          ! Note that the set of point is symmetric around the central point, while weights are taken into account
749          do iq=0,nq
750            qq2(:) = qq1(:)+ (DBLE(iq)*invnq-half)*invnq
751            bb4sinpiqq2_2(:,iq) =bb(:)*four*SIN(pi*qq2(:))**2; sin2piqq2(:,iq)=SIN(two_pi*qq2(:))
752          end do
753          do iqx2=0,nq
754            qq2(1)=qq1(1) + (DBLE(iqx2)*invnq-half ) *invnq
755            do iqy2=0,nq
756              qq2(2)=qq1(2) + (DBLE(iqy2)*invnq-half ) *invnq
757              weightxy=invnq3*weightq
758              if (iqx2 == 0 .or. iqx2 == nq) weightxy = weightxy*half
759              if (iqy2 == 0 .or. iqy2 == nq) weightxy = weightxy*half
760              do iqz2=0,nq
761                qq2(3) = qq1(3) + (DBLE(iqz2)*invnq - half) * invnq
762                weightxyz = weightxy
763                if (iqz2 == 0 .or. iqz2 == nq) weightxyz = weightxy*half
764                !
765                ! Treat the remaining divergence in the origin as if it would be a spherical integration of 1/q^2
766                if (iqx1/=nqhalf .or. iqy1/=nqhalf .or. iqz1/=nqhalf .or. &
767                    iqx2/=nqhalf .or. iqy2/=nqhalf .or. iqz2/=nqhalf ) then
768                  !integratefaux=integratefaux+ faux(qq2, rcut, b1, b2, b3) *invnq**6
769                  integratefaux = integratefaux + &
770                    faux_fast(qq2, bb4sinpiqq2_2(1,iqx2), bb4sinpiqq2_2(2,iqy2), bb4sinpiqq2_2(3,iqz2), &
771                              sin2piqq2(1,iqx2), sin2piqq2(2,iqy2), sin2piqq2(3,iqz2), rcut, b1, b2, b3) * weightxyz
772                else
773                   integratefaux = integratefaux + 7.7955* ((two_pi)**3/ucvol*invnq3*invnq3 )**(-2./3.) *invnq3*invnq3
774                end if
775              end do
776            end do
777          end do
778        else
779         ! integratefaux=integratefaux+faux(qq1, rcut, b1, b2, b3)*invnq**3
780         integratefaux = integratefaux + &
781           faux_fast(qq1, bb4sinpiqq_2(1,iqx1), bb4sinpiqq_2(2,iqy1), bb4sinpiqq_2(3,iqz1), &
782                     sin2piqq(iqx1), sin2piqq(iqy1), sin2piqq(iqz1), rcut, b1, b2, b3) * weightq
783        end if
784      end do
785    end do
786  end do
787 
788  call xmpi_sum(integratefaux, comm, ierr)
789 
790 end function integratefaux

m_vcoul/mc_free [ Functions ]

[ Top ] [ m_vcoul ] [ Functions ]

NAME

 mc_free

FUNCTION

  Free dynamic memory

SOURCE

1434 subroutine mc_free(mc)
1435 
1436 !Arguments ------------------------------------
1437  class(mc_t),intent(inout) :: mc
1438 
1439 ! *************************************************************************
1440 
1441  ABI_SFREE(mc%qran)
1442 
1443 end subroutine mc_free

m_vcoul/mc_init [ Functions ]

[ Top ] [ m_vcoul ] [ Functions ]

NAME

 mc_init

FUNCTION

SOURCE

1203 subroutine mc_init(mc, rprimd, ucvol, gprimd, gmet, kptrlatt)
1204 
1205 !Arguments ------------------------------------
1206  class(mc_t),intent(out) :: mc
1207  real(dp),intent(in) :: rprimd(3,3), ucvol, gprimd(3,3), gmet(3,3)
1208  integer,intent(in) :: kptrlatt(3,3)
1209 
1210 !Local variables-------------------------------
1211  integer,parameter :: ncell=3
1212  integer :: nseed, i1,i2,i3,imc
1213  real(dp) :: lmin,vlength, ucvol_sc
1214  real(dp) :: rprimd_sc(3,3),gprimd_sc(3,3),gmet_sc(3,3),rmet_sc(3,3), qcart2red(3,3)
1215  real(dp) :: qtmp(3),qmin(3),qmin_cart(3)
1216  integer, allocatable :: seed(:)
1217 
1218 ! *************************************************************************
1219 
1220  mc%gmet = gmet
1221  mc%ucvol = ucvol
1222 
1223  ! Supercell defined by the k-mesh
1224  rprimd_sc(:,:) = MATMUL(rprimd, kptrlatt)
1225  call metric(gmet_sc, gprimd_sc, -1, rmet_sc, rprimd_sc, ucvol_sc)
1226 
1227  qcart2red(:,:) = two_pi * gprimd(:,:)
1228  call matrginv(qcart2red, 3, 3)
1229 
1230  ! Find the largest sphere inside the miniBZ
1231  ! in order to integrate the divergence analytically
1232  mc%q0sph = HUGE(one)
1233  do i1 = -ncell+1, ncell
1234    qtmp(1) = dble(i1) * 0.5_dp
1235    do i2 = -ncell+1, ncell
1236      qtmp(2) = dble(i2) * 0.5_dp
1237      do i3 = -ncell+1, ncell
1238        qtmp(3) = dble(i3) * 0.5_dp
1239        if (i1 == 0 .AND. i2 == 0 .AND. i3 == 0) cycle
1240        vlength = normv(qtmp, gmet_sc, 'G')
1241        if (vlength < mc%q0sph) mc%q0sph = vlength
1242      enddo
1243    enddo
1244  enddo
1245 
1246  ! Setup the random vectors for the Monte Carlo sampling of the miniBZ at q = 0
1247  mc%nmc_max = 2500000
1248  ABI_MALLOC(mc%qran,(3, mc%nmc_max))
1249  call random_seed(size=nseed)
1250  ABI_MALLOC(seed, (nseed))
1251  do i1=1,nseed
1252    seed(i1) = NINT(SQRT(DBLE(i1) * 103731))
1253  end do
1254  call random_seed(put=seed)
1255  call random_number(mc%qran)
1256  ABI_FREE(seed)
1257 
1258  ! Overide the first "random vector" with 0
1259  mc%qran(:,1) = zero
1260 
1261  ! Fold qran into the Wignez-Seitz cell around q = 0
1262  do imc=2,mc%nmc_max
1263    lmin = HUGE(one)
1264    do i1 = -ncell+1, ncell
1265      qtmp(1) = mc%qran(1,imc) + dble(i1)
1266      do i2 = -ncell+1, ncell
1267        qtmp(2) = mc%qran(2,imc) + dble(i2)
1268        do i3 = -ncell+1, ncell
1269          qtmp(3) = mc%qran(3,imc) + dble(i3)
1270          vlength = normv(qtmp, gmet_sc, 'G')
1271          if (vlength < lmin) then
1272            lmin = vlength
1273            ! Get the q-vector in cartesian coordinates
1274            qmin_cart(:) = two_pi * MATMUL( gprimd_sc(:,:) , qtmp )
1275            ! Transform it back to the reciprocal space
1276            qmin(:) = MATMUL( qcart2red , qmin_cart )
1277          end if
1278        enddo
1279      enddo
1280    enddo
1281 
1282    mc%qran(:,imc) = qmin(:)
1283  enddo
1284 
1285 end subroutine mc_init

m_vcoul/mc_integrate [ Functions ]

[ Top ] [ m_vcoul ] [ Functions ]

NAME

 mc_integrate

FUNCTION

SOURCE

1298 subroutine mc_integrate(mc, mode, qibz, ng, gvec, rcut2, nkbz, vcoul, comm)
1299 
1300 !Arguments ------------------------------------
1301  class(mc_t),intent(in) :: mc
1302  real(dp),intent(in) :: rcut2
1303  integer,intent(in) :: nkbz, ng, comm
1304  character(len=*),intent(in) :: mode
1305  real(dp),intent(in) :: qibz(3)
1306  integer,intent(in) :: gvec(3, ng)
1307  real(dp),intent(out) :: vcoul(ng)
1308 
1309 !Local variables-------------------------------
1310  integer,parameter :: master = 0
1311  integer :: ig, ig0, imc, nmc, my_rank, nprocs, ierr
1312  logical :: q_is_gamma
1313  real(dp)  :: qpg2, qpg(3)
1314 
1315 ! *************************************************************************
1316 
1317  my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
1318  q_is_gamma = all(abs(qibz) < tol16)
1319 
1320  ! Find index of G=0 in gvec.
1321  ig0 = -1
1322  do ig=1,ng
1323    if (all(gvec(:, ig) == 0)) then
1324       ig0 = ig; exit
1325    end if
1326  end do
1327  ABI_CHECK(ig0 /= -1, "Cannot find G=0 in gvec!")
1328 
1329  vcoul = zero
1330 
1331  select case (trim(mode))
1332  case('MINIBZ')
1333 
1334    do ig=1,ng
1335      if (mod(ig, nprocs) /= my_rank) cycle ! MPI parallelism.
1336      if (q_is_gamma .and. ig == ig0) cycle
1337      qpg(:) = qibz(:) + gvec(:,ig)
1338      qpg2 = normv(qpg, mc%gmet, 'G')**2
1339      nmc = adapt_nmc(mc%nmc_max, qpg2)
1340      do imc=1,nmc
1341        qpg(:) = qibz(:) + gvec(:,ig) + mc%qran(:,imc)
1342        qpg2 = normv(qpg, mc%gmet, 'G')**2
1343        vcoul(ig) = vcoul(ig) + four_pi / qpg2 / REAL(nmc, dp)
1344      end do
1345    end do ! ig
1346 
1347    if (q_is_gamma .and. my_rank == master) then
1348      ! Override ig0 component
1349      vcoul(ig0) = four_pi**2 * nkbz * mc%ucvol / ( 8.0_dp * pi**3 ) * mc%q0sph
1350      do imc=1,mc%nmc_max
1351        qpg(:) = qibz(:) + gvec(:,ig0) + mc%qran(:,imc)
1352        qpg2 = normv(qpg, mc%gmet, 'G')**2
1353        if (qpg2 > mc%q0sph ** 2) vcoul(ig0) = vcoul(ig0) + four_pi / qpg2 / REAL(mc%nmc_max, dp)
1354      end do
1355    end if
1356 
1357  case('MINIBZ-ERFC')
1358 
1359    do ig=1,ng
1360      if (mod(ig, nprocs) /= my_rank) cycle ! MPI parallelism.
1361      if (q_is_gamma .and. ig == ig0) cycle
1362      qpg(:) = qibz(:) + gvec(:,ig)
1363      qpg2 = normv(qpg, mc%gmet, 'G')**2
1364      nmc = adapt_nmc(mc%nmc_max, qpg2)
1365      do imc=1,nmc
1366        qpg(:) = qibz(:) + gvec(:,ig) + mc%qran(:,imc)
1367        qpg2 = normv(qpg, mc%gmet, 'G')**2
1368        vcoul(ig) = vcoul(ig) + four_pi / qpg2 / REAL(nmc,dp) * (  one - EXP( -0.25d0 * rcut2 * qpg2 ) )
1369      end do
1370    end do ! ig
1371 
1372    if (q_is_gamma .and. my_rank == master) then
1373      ! Override ig0 component
1374      vcoul(ig0) = four_pi**2 * nkbz * mc%ucvol / ( 8.0_dp * pi**3 ) &
1375         * ( mc%q0sph - SQRT(pi/rcut2) * abi_derf(0.5_dp*SQRT(rcut2)*mc%q0sph) )
1376      do imc=1,mc%nmc_max
1377        qpg(:) = qibz(:) + gvec(:,ig0) + mc%qran(:,imc)
1378        qpg2 = normv(qpg, mc%gmet, 'G')**2
1379        if (qpg2 > mc%q0sph**2) then
1380          vcoul(ig0) = vcoul(ig0) + four_pi / qpg2 / REAL(mc%nmc_max,dp) * (one - EXP( -0.25d0 * rcut2 * qpg2))
1381        end if
1382      end do
1383    end if
1384 
1385  case('MINIBZ-ERF')
1386 
1387    do ig=1,ng
1388      if (mod(ig, nprocs) /= my_rank) cycle ! MPI parallelism.
1389      if (q_is_gamma .and. ig == ig0) cycle
1390      qpg(:) = qibz(:) + gvec(:,ig)
1391      qpg2 = normv(qpg, mc%gmet, 'G')**2
1392      nmc = adapt_nmc(mc%nmc_max, qpg2)
1393      do imc=1,nmc
1394        qpg(:) = qibz(:) + gvec(:,ig) + mc%qran(:,imc)
1395        qpg2 = normv(qpg, mc%gmet, 'G')**2
1396        vcoul(ig) = vcoul(ig) + four_pi / qpg2 / REAL(nmc,dp) * EXP( -0.25d0 * rcut2 * qpg2 )
1397      end do
1398    end do ! ig
1399 
1400    if (q_is_gamma .and. my_rank == master) then
1401      ! Override ig=ig0 component
1402      vcoul(ig0) = four_pi**2 * nkbz * mc%ucvol / ( 8.0_dp * pi**3 ) * SQRT(pi/rcut2) * abi_derf(0.5_dp*SQRT(rcut2)*mc%q0sph)
1403 
1404      do imc=1,mc%nmc_max
1405        qpg(:) = qibz(:) + gvec(:,ig0) + mc%qran(:,imc)
1406        qpg2 = normv(qpg, mc%gmet, 'G')**2
1407        if (qpg2 > mc%q0sph**2) then
1408          vcoul(ig0) = vcoul(ig0) + four_pi / qpg2 / REAL(mc%nmc_max,dp) *  EXP( -0.25d0 * rcut2 * qpg2 )
1409        end if
1410      end do
1411    end if
1412 
1413  case default
1414    ABI_ERROR(sjoin("Invalid mode:", mode))
1415  end select
1416 
1417  ! Collect result on each MPI proc.
1418  call xmpi_sum(vcoul, comm, ierr)
1419 
1420 end subroutine mc_integrate

m_vcoul/mc_t [ Types ]

[ Top ] [ m_vcoul ] [ Types ]

NAME

  mc_t

FUNCTION

 Mimicking the BerkeleyGW technique
 A Monte-Carlo sampling of each miniBZ surrounding each (q+G) point
 However:
    - extended to multiple shifts
    - with an adaptative number of MonteCarlo sampling points

SOURCE

169 type, public :: mc_t
170 
171   integer :: nmc_max = -1
172 
173   real(dp) :: q0sph = -one
174 
175   real(dp) :: ucvol = -one
176 
177   real(dp) :: gmet(3,3) = -one
178 
179   real(dp),allocatable :: qran(:,:)
180   ! (3, nmc_max)
181 
182 contains
183   procedure :: init => mc_init
184   procedure :: integrate => mc_integrate
185   procedure :: free => mc_free
186 end type mc_t

m_vcoul/surface_setup [ Functions ]

[ Top ] [ m_vcoul ] [ Functions ]

NAME

  surface_setup

FUNCTION

SOURCE

643 subroutine surface_setup(cryst, vcutgeo, alpha, rcut, pdir, opt_slab)
644 
645  type(crystal_t),intent(in) :: cryst
646  real(dp),intent(in) :: vcutgeo(3)
647  real(dp),intent(out) :: alpha(3)
648  real(dp),intent(inout) :: rcut
649  integer,intent(out) :: pdir(3), opt_slab
650 
651 !Local variables-------------------------------
652  integer :: ii
653  real(dp) :: check
654  character(len=500) :: msg
655 
656 ! *************************************************************************
657 
658  ABI_CHECK(count(vcutgeo /= zero) == 2, "Wrong vcutgeo")
659 
660  ! Default is Beigi's method.
661  opt_slab = 1; if (any(vcutgeo < zero)) opt_slab = 2
662  pdir(:) = zero; alpha(:)=zero
663  do ii=1,3
664    check = vcutgeo(ii)
665    if (abs(check) > zero) then
666      ! Use Rozzi's method with a finite surface along x-y
667      pdir(ii) = 1
668      if (check < zero) alpha(ii) = normv(check * cryst%rprimd(:,ii), cryst%rmet, 'R')
669    end if
670  end do
671 
672  ! In Beigi's method, the surface must be along x-y and R must be L_Z/2.
673  if (opt_slab == 1) then
674    msg = "2D Beigi method, the periodicity must be in the x-y plane. Modify vcutgeo and/or your geometry."
675    ABI_CHECK(all(pdir == [1, 1, 0]), msg)
676    rcut = half*SQRT(DOT_PRODUCT(cryst%rprimd(:,3), cryst%rprimd(:,3)))
677  end if
678 
679 end subroutine surface_setup

m_vcoul/vcgen_free [ Functions ]

[ Top ] [ m_vcoul ] [ Functions ]

NAME

 vcgen_free

FUNCTION

  Free dynamic memory.

SOURCE

1938 subroutine vcgen_free(vcgen)
1939 
1940 !Arguments ------------------------------------
1941  class(vcgen_t),intent(inout) :: vcgen
1942 
1943 ! *************************************************************************
1944 
1945  call vcgen%mc%free()
1946 
1947 end subroutine vcgen_free

m_vcoul/vcgen_get_vc_sqrt [ Functions ]

[ Top ] [ m_vcoul ] [ Functions ]

NAME

 vcgen_get_vc_sqrt

FUNCTION

  Compute sqrt(vc(q, g))

SOURCE

1851 subroutine vcgen_get_vc_sqrt(vcgen, qpt, npw, gvec, q0, cryst, vc_sqrt, comm)
1852 
1853 !Arguments ------------------------------------
1854  class(vcgen_t),intent(in) :: vcgen
1855  real(dp),intent(in) :: qpt(3), q0(3)
1856  integer,intent(in) :: npw, gvec(3,npw), comm
1857  type(crystal_t),intent(in) :: cryst
1858  complex(gwpc),intent(out) :: vc_sqrt(npw)
1859 
1860 !Local variables-------------------------------
1861  integer :: ig, ig0
1862  real(dp) :: rcut2
1863  logical :: q_is_gamma
1864  real(dp),allocatable :: vcoul(:)
1865 
1866 ! *************************************************************************
1867 
1868  q_is_gamma = normv(qpt, cryst%gmet, "G") < GW_TOLQ0
1869 
1870  ! Find index of G=0 in gvec.
1871  ig0 = -1
1872  do ig=1,npw
1873    if (all(gvec(:,ig) == 0)) then
1874     ig0 = ig; exit
1875    end if
1876  end do
1877  ABI_CHECK(ig0 /= -1, "Cannot find G=0 in gvec!")
1878 
1879  ABI_MALLOC(vcoul, (npw))
1880 
1881  select case (trim(vcgen%mode))
1882  case ('MINIBZ', 'MINIBZ-ERFC', 'MINIBZ-ERF')
1883    rcut2 = vcgen%rcut**2
1884    call vcgen%mc%integrate(vcgen%mode, qpt, npw, gvec, rcut2, vcgen%nkbz, vcoul, comm)
1885 
1886  case ('SPHERE')
1887    call cutoff_sphere(qpt, npw, gvec, cryst%gmet, vcgen%rcut, vcoul)
1888 
1889  case ('CYLINDER')
1890    call cutoff_cylinder(qpt, npw, gvec, vcgen%rcut, vcgen%hcyl, vcgen%pdir, &
1891                         vcgen%boxcenter, cryst%rprimd, vcoul, vcgen%opt_cylinder, comm)
1892 
1893  case ('SLAB')
1894    call cutoff_slab(qpt, npw, gvec, cryst%gprimd, vcgen%rcut, &
1895                     vcgen%boxcenter, vcgen%pdir, vcgen%alpha, vcoul, vcgen%opt_slab)
1896 
1897  case ('CRYSTAL', 'AUXILIARY_FUNCTION', "AUX_GB", "ERF", "ERFC")
1898    ! Compute |q+G| with special treatment of (q=0, g=0).
1899    do ig=1,npw
1900      !if (q_is_gamma) then
1901      if (q_is_gamma .and. ig == ig0) then
1902        vcoul(ig) = normv(q0 + gvec(:,ig), cryst%gmet, "G")
1903      else
1904        vcoul(ig) = normv(qpt + gvec(:,ig), cryst%gmet, "G")
1905      end if
1906    end do
1907 
1908    if (vcgen%mode == "ERF") then
1909      vcoul(:)  = four_pi/(vcoul(:)**2) *  EXP( -0.25d0 * (vcgen%rcut*vcoul(:))**2 )
1910    else if (vcgen%mode == "ERFC") then
1911      vcoul(:) = four_pi/(vcoul(:)**2) * ( one - EXP( -0.25d0 * (vcgen%rcut*vcoul(:))**2 ) )
1912    else
1913      vcoul = four_pi/vcoul**2
1914    end if
1915 
1916  case default
1917    ABI_BUG(sjoin('Unsupported cutoff mode:', vcgen%mode))
1918  end select
1919 
1920  ! Store final results in complex array as Rozzi's cutoff can give real negative values
1921  vc_sqrt = SQRT(CMPLX(vcoul, zero))
1922  ABI_FREE(vcoul)
1923 
1924 end subroutine vcgen_get_vc_sqrt

m_vcoul/vcgen_init [ Functions ]

[ Top ] [ m_vcoul ] [ Functions ]

NAME

 vcgen_init

FUNCTION

SOURCE

1714 subroutine vcgen_init(vcgen, cryst, kptrlatt, nkbz, nqibz, nqbz, qbz, rcut, gw_icutcoul, vcutgeo, ecut, comm)
1715 
1716 !Arguments ------------------------------------
1717  class(vcgen_t),intent(out) :: vcgen
1718  type(crystal_t),intent(in) :: cryst
1719  integer,intent(in) :: kptrlatt(3,3), nkbz, nqibz, nqbz, gw_icutcoul
1720  real(dp),intent(in) :: qbz(3,nqbz), rcut, ecut, vcutgeo(3)
1721  integer,intent(in) :: comm
1722 
1723 !Local variables-------------------------------
1724  integer,parameter :: istwfk1 = 1
1725  integer :: gvec0(3) ! npw_,
1726  real(dp) :: q0_vol, bz_geometry_factor, rcut2
1727  character(len=500) :: msg
1728  real(dp) :: vcoul0(1), q_gamma(3)
1729  !integer,allocatable :: gvec_(:,:)
1730 
1731 ! *************************************************************************
1732 
1733  ABI_UNUSED([ecut])
1734 
1735  ! Save dimension and other useful quantities in Vcp
1736  vcgen%rcut      = rcut                 ! Cutoff radius for cylinder.
1737  vcgen%hcyl      = zero                 ! Length of finite cylinder (Rozzi's method, default is Beigi).
1738  vcgen%boxcenter = zero                 ! Boxcenter at the moment is supposed to be at the origin.
1739  vcgen%vcutgeo   = vcutgeo(:)           ! Info on the orientation and extension of the cutoff region.
1740  vcgen%nkbz      = nkbz
1741 
1742  call gw_icutcoul_to_mode(gw_icutcoul, vcgen%mode)
1743  q_gamma = zero
1744  gvec0 = 0
1745 
1746  select case (trim(vcgen%mode))
1747  case ('MINIBZ', 'MINIBZ-ERFC', 'MINIBZ-ERF')
1748    call vcgen%mc%init(cryst%rprimd, cryst%ucvol, cryst%gprimd, cryst%gmet, kptrlatt)
1749    rcut2 = vcgen%rcut**2
1750    call vcgen%mc%integrate(vcgen%mode, q_gamma, 1, gvec0, rcut2, nkbz, vcoul0, xmpi_comm_self)
1751    ! Treat the limit q --> 0.
1752    vcgen%i_sz = vcoul0(1)
1753 
1754  case ('SPHERE')
1755    ! A non-positive value of rcut activates the recipe of Spencer & Alavi, PRB 77, 193110 (2008) [[cite:Spencer2008]].
1756    if (vcgen%rcut < tol12) then
1757      vcgen%rcut = (cryst%ucvol * nkbz * 3.d0 / four_pi) ** third
1758      write(msg,'(2a,2x,f8.4,a)')ch10,' Using calculated rcut: ',vcgen%rcut,' to have same volume as the BvK crystal'
1759      call wrtout(std_out, msg)
1760    end if
1761    vcgen%vcutgeo = zero
1762 
1763    ! Treat the limit q --> 0
1764    ! The small cube is approximated by a sphere, while vc(q=0) = 2piR**2.
1765    ! if a single q-point is used, the expression for the volume is exact.
1766    vcgen%i_sz = two_pi * vcgen%rcut**2
1767 
1768  case ('CYLINDER')
1769    call cylinder_setup(cryst, vcgen%vcutgeo, vcgen%hcyl, vcgen%pdir, vcgen%opt_cylinder)
1770 
1771    ! If Beigi, treat the limit q --> 0.
1772    if (vcgen%opt_cylinder == 1) then
1773      call beigi_cylinder_limit(vcgen%opt_cylinder, cryst, nqibz, nkbz, &
1774                                vcgen%rcut, vcgen%hcyl, vcgen%boxcenter, vcgen%pdir, vcgen%i_sz)
1775    else
1776      ! In Rozzi's method the lim q+G --> 0 is finite.
1777      call cutoff_cylinder(q_gamma, 1, gvec0, vcgen%rcut, vcgen%hcyl, vcgen%pdir,&
1778                           vcgen%boxcenter, cryst%rprimd, vcoul0, vcgen%opt_cylinder, xmpi_comm_self)
1779      vcgen%i_sz = vcoul0(1)
1780    end if
1781 
1782  case ('SLAB')
1783    call surface_setup(cryst, vcgen%vcutgeo, vcgen%alpha, vcgen%rcut, vcgen%pdir, vcgen%opt_slab)
1784 
1785    ! If Beigi, treat the limit q --> 0.
1786    if (vcgen%opt_slab == 1) then
1787      ! Integrate numerically in the plane close to 0
1788      call beigi_surface_limit(vcgen%opt_slab, cryst, nqibz, nkbz, vcgen%rcut, vcgen%alpha, &
1789                               vcgen%boxcenter, vcgen%pdir, vcgen%i_sz)
1790    else
1791      ! In Rozzi's method the lim q+G --> 0 is finite.
1792      call cutoff_slab(q_gamma, 1, gvec0, cryst%gprimd, vcgen%rcut, &
1793                       vcgen%boxcenter, vcgen%pdir, vcgen%alpha, vcoul0, vcgen%opt_slab)
1794      vcgen%i_sz = vcoul0(1)
1795    end if
1796 
1797  case ('CRYSTAL', 'AUXILIARY_FUNCTION', "AUX_GB")
1798 
1799    if (vcgen%mode == "CRYSTAL") then
1800      ! Analytic integration of 4pi/q^2 over the volume element:
1801      ! $4pi/V \int_V d^3q 1/q^2 =4pi bz_geometric_factor V^(-2/3)$
1802      ! i_sz=4*pi*bz_geometry_factor*q0_vol**(-two_thirds) where q0_vol= V_BZ/N_k
1803      ! bz_geometry_factor: sphere=7.79, fcc=7.44, sc=6.188, bcc=6.946, wz=5.255 (see gwa.pdf, appendix A.4)
1804      q0_vol = (two_pi) **3 / (nkbz*cryst%ucvol); bz_geometry_factor=zero
1805      vcgen%i_sz = four_pi*7.44*q0_vol**(-two_thirds)
1806 
1807    else if (vcgen%mode == "AUXILIARY_FUNCTION") then
1808      ! Numerical integration of the exact-exchange divergence through the
1809      ! auxiliary function of Carrier et al. PRB 75, 205126 (2007) [[cite:Carrier2007]].
1810      vcgen%i_sz = carrier_isz(cryst, nqbz, qbz, rcut, comm)
1811 
1812    else if (vcgen%mode == "AUX_GB") then
1813      ! We use the auxiliary function of a Gygi-Baldereschi variant [[cite:Gigy1986]]
1814      ABI_ERROR("AUX_GB not implemented in vcgen_init")
1815      !call get_kg(kk_bz, istwfk1, ecut, cryst%gmet, npw_, gvec_)
1816      !vcgen%i_sz = gygi_baldereschi_isz(cryst, nqbz, qbz, ecut, ng, gvec_)
1817      !ABI_FREE(gvec_)
1818 
1819    else
1820      ABI_ERROR(sjoin("Need treatment of 1/q^2 singularity! for mode", vcgen%mode))
1821    end if
1822 
1823  case ('ERF')
1824    vcgen%i_sz = carrier_isz(cryst, nqbz, qbz, rcut, xmpi_comm_self)
1825 
1826  case ('ERFC')
1827    ! === Treat 1/q^2 singularity ===
1828    ! * There is NO singularity in this case.
1829    vcgen%i_sz = pi * vcgen%rcut**2 ! Final result stored here
1830 
1831  case default
1832    ABI_BUG(sjoin('Unsupported cutoff mode:', vcgen%mode))
1833  end select
1834 
1835  !call wrtout(std_out, sjoin(" vcgen%i_sz", ftoa(vcgen%i_sz)))
1836 
1837 end subroutine vcgen_init

m_vcoul/vcgen_t [ Types ]

[ Top ] [ m_vcoul ] [ Types ]

NAME

  vcgen_t

FUNCTION

SOURCE

197 type, public :: vcgen_t
198 
199   integer :: nkbz = -1
200    ! Number of k-points in full BZ.
201 
202   integer :: opt_cylinder
203 
204   integer :: opt_slab
205 
206   real(dp) :: alpha(3) = -one
207    ! Lenght of the finite surface.
208 
209   real(dp) :: rcut = -one
210    ! Cutoff radius.
211 
212   real(dp) :: hcyl = -one
213    ! Length of the finite cylinder along the periodic dimension
214 
215   character(len=50) :: mode
216    ! String defining the cutoff mode
217 
218   integer :: pdir(3)
219    ! 1 if the system is periodic along this direction
220 
221   real(dp) :: boxcenter(3) = -1
222    ! 1 if the point in inside the cutoff region 0 otherwise
223    ! Reduced coordinates of the center of the box (input variable)
224 
225   real(dp) :: vcutgeo(3) = huge(one)
226    ! For each reduced direction gives the length of the finite system
227    ! 0 if the system is infinite along this direction.
228    ! negative values indicate that a finite size has to be used.
229 
230   real(dp) :: i_sz = huge(one)
231    ! Value of the integration of the Coulomb singularity 4\pi/V_BZ \int_BZ d^3q 1/q^2
232 
233   type(mc_t) :: mc
234    ! Monte carlo integrator.
235 
236 contains
237   procedure :: init => vcgen_init                  ! Initialize the object
238   procedure :: get_vc_sqrt => vcgen_get_vc_sqrt    ! Compute sqrt(vc(q,g))
239   procedure :: free => vcgen_free                  ! Free
240   !procedure :: print => vcgen_print
241 end type vcgen_t

m_vcoul/vcoul_free [ Functions ]

[ Top ] [ m_vcoul ] [ Functions ]

NAME

 vcoul_free

FUNCTION

  Free memory

SOURCE

1177 subroutine vcoul_free(Vcp)
1178 
1179 !Arguments ------------------------------------
1180  class(vcoul_t),intent(inout) :: Vcp
1181 
1182 ! *************************************************************************
1183 
1184  ABI_SFREE(Vcp%qibz)
1185  ABI_SFREE(Vcp%qlwl)
1186  ABI_SFREE(Vcp%vc_sqrt)
1187  ABI_SFREE(Vcp%vc_sqrt_resid)
1188  ABI_SFREE(Vcp%vcqlwl_sqrt)
1189 
1190 end subroutine vcoul_free

m_vcoul/vcoul_init [ Functions ]

[ Top ] [ m_vcoul ] [ Functions ]

NAME

 vcoul_init

FUNCTION

 Perform general check and initialize the data type containing information on the cutoff technique
 Note %vc_sqrt_resid and %i_sz_resid are simply initialized at the same value as %vc_sqrt and %i_sz

INPUTS

  Gsph=Info of the G sphere.
  Qmesh=Info on the q-point sampling.
  Kmesh=Info on the k-point sampling.
  rcut=Cutoff radius for the cylinder.
  gw_icutcoul=Option of the cutoff technique.
  vcutgeo(3)= Info on the orientation and extension of the cutoff region.
  ng=Number of G-vectors to be used to describe the Coulomb interaction
  nqlwl=Number of point around Gamma for treatment of long-wavelength limit
  qlwl(3,nqlwl)= The nqlwl "small" q-points
  comm=MPI communicator.

OUTPUT

  vcp=Datatype gathering information on the Coulomb interaction.

SOURCE

308 subroutine vcoul_init(vcp, Gsph, Cryst, Qmesh, Kmesh, rcut, gw_icutcoul, vcutgeo, ecut, ng, nqlwl, qlwl, comm)
309 
310 !Arguments ------------------------------------
311 !scalars
312  class(vcoul_t),intent(out) :: vcp
313  integer,intent(in) :: ng,nqlwl, gw_icutcoul, comm
314  real(dp),intent(in) :: rcut, ecut
315  type(kmesh_t),target,intent(in) :: Kmesh, Qmesh
316  type(gsphere_t),target,intent(in) :: Gsph
317  type(crystal_t),intent(in) :: Cryst
318 !arrays
319  real(dp),intent(in) :: qlwl(3,nqlwl),vcutgeo(3)
320 
321 !Local variables-------------------------------
322 !scalars
323  integer,parameter :: master=0
324  integer :: nqibz, nqbz, nkbz, iqlwl, iq_ibz
325  integer :: opt_cylinder,my_rank,nprocs, opt_slab
326  real(dp) :: bz_geometry_factor,q0_vol, rcut2
327  character(len=500) :: msg
328  type(mc_t) :: mc
329 !arrays
330  integer, contiguous, pointer :: gvec(:,:)
331  real(dp) :: a1(3),a2(3),a3(3),b1(3),b2(3),b3(3)
332  real(dp),allocatable :: vcoul(:,:),vcoul_lwl(:,:)
333  real(dp),contiguous, pointer :: qibz(:,:), qbz(:,:)
334 
335 ! *************************************************************************
336 
337  ! === Test if the first q-point is zero ===
338  ! FIXME this wont work if nqptdm/=0
339  !if (normv(Qmesh%ibz(:,1),gmet,'G') < GW_TOLQ0)) STOP 'vcoul_init, non zero first point '
340 
341  my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
342 
343  nqibz = qmesh%nibz; nqbz = qmesh%nbz
344  qibz => qmesh%ibz; qbz => qmesh%bz
345  nkbz = kmesh%nbz
346 
347  ! Save dimension and other useful quantities in vcp
348  vcp%ng        = ng                   ! Number of G-vectors in the Coulomb matrix elements.
349  vcp%nqibz     = nqibz                ! Number of irred q-point.
350  vcp%nqlwl     = nqlwl                ! Number of small q-directions to deal with singularity and non Analytic behavior.
351  vcp%rcut      = rcut                 ! Cutoff radius for cylinder.
352  vcp%hcyl      = zero                 ! Length of finite cylinder (Rozzi"s method, default is Beigi).
353  vcp%ucvol     = cryst%ucvol          ! Unit cell volume.
354  vcp%rprimd    = Cryst%rprimd(:,:)    ! Dimensional direct lattice.
355  vcp%boxcenter = zero                 ! Boxcenter at the moment is supposed to be at the origin.
356  vcp%vcutgeo   = vcutgeo(:)           ! Info on the orientation and extension of the cutoff region.
357 
358  gvec => Gsph%gvec
359 
360  call gw_icutcoul_to_mode(gw_icutcoul, vcp%mode)
361 
362  ABI_MALLOC(vcp%qibz, (3, nqibz))
363  vcp%qibz = Qmesh%ibz(:,:)
364  ABI_MALLOC(vcp%qlwl, (3, nqlwl))
365  vcp%qlwl = qlwl(:,:)
366 
367  ! ===============================================
368  ! == Calculation of the FT of the Coulomb term ==
369  ! ===============================================
370  a1 = cryst%rprimd(:,1); a2 = cryst%rprimd(:,2); a3 = cryst%rprimd(:,3)
371  b1 = two_pi * cryst%gprimd(:,1); b2 = two_pi * cryst%gprimd(:,2); b3 = two_pi * cryst%gprimd(:,3)
372 
373  ABI_MALLOC(vcoul    , (ng, nqibz))
374  ABI_MALLOC(vcoul_lwl, (ng, nqlwl))
375 
376  select case (trim(vcp%mode))
377  case ('MINIBZ', 'MINIBZ-ERFC', 'MINIBZ-ERF')
378    call mc%init(cryst%rprimd, cryst%ucvol, cryst%gprimd, cryst%gmet, kmesh%kptrlatt)
379 
380    rcut2 = vcp%rcut**2
381    do iq_ibz=1,nqibz
382      call mc%integrate(vcp%mode, qibz(:, iq_ibz), ng, gvec, rcut2, nkbz, vcoul(:, iq_ibz), comm)
383    end do
384 
385    ! Treat the limit q --> 0
386    vcp%i_sz = vcoul(1, 1)
387 
388    do iqlwl=1,nqlwl
389      call mc%integrate(vcp%mode, qlwl(:, iqlwl), ng, gvec, rcut2, nkbz, vcoul_lwl(:, iqlwl), comm)
390    end do
391 
392    call mc%free()
393 
394  case ('SPHERE')
395    ! A non-positive value of rcut activates the recipe of Spencer & Alavi, PRB 77, 193110 (2008) [[cite:Spencer2008]].
396    if (vcp%rcut < tol12) then
397      vcp%rcut = (cryst%ucvol * nkbz * 3.d0 / four_pi) ** third
398      write(msg,'(2a,2x,f8.4,a)')ch10,' Using calculated rcut: ',vcp%rcut,' to have same volume as the BvK crystal'
399      call wrtout(std_out, msg)
400    end if
401    vcp%vcutgeo = zero
402 
403    do iq_ibz=1,nqibz
404      call cutoff_sphere(qibz(:,iq_ibz), ng, gvec, cryst%gmet, vcp%rcut, vcoul(:,iq_ibz))
405    end do
406 
407    ! q-points for optical limit.
408    do iqlwl=1,nqlwl
409      call cutoff_sphere(qlwl(:,iqlwl), ng, gvec, cryst%gmet, vcp%rcut, vcoul_lwl(:,iqlwl))
410    end do
411 
412    ! Treat the limit q --> 0
413    ! The small cube is approximated by a sphere, while vc(q=0) = 2piR**2.
414    ! if a single q-point is used, the expression for the volume is exact.
415    vcp%i_sz = two_pi * vcp%rcut**2
416    call vcp%print(unit=ab_out)
417 
418  case ('CYLINDER')
419    call cylinder_setup(cryst, vcp%vcutgeo, vcp%hcyl, vcp%pdir, opt_cylinder)
420 
421    do iq_ibz=1,nqibz
422      call cutoff_cylinder(qibz(:,iq_ibz), ng, gvec, vcp%rcut, vcp%hcyl, vcp%pdir,&
423                           vcp%boxcenter, Cryst%rprimd, vcoul(:,iq_ibz), opt_cylinder, comm)
424    end do
425 
426    ! q-points for optical limit.
427    do iqlwl=1,nqlwl
428      call cutoff_cylinder(qlwl(:,iqlwl), ng, gvec, vcp%rcut, vcp%hcyl, vcp%pdir,&
429                           vcp%boxcenter, Cryst%rprimd, vcoul_lwl(:,iqlwl), opt_cylinder, comm)
430    end do
431 
432    ! If Beigi, treat the limit q --> 0.
433    if (opt_cylinder == 1) then
434      call beigi_cylinder_limit(opt_cylinder, cryst, nqibz, nkbz, vcp%rcut, vcp%hcyl, vcp%boxcenter, vcp%pdir, vcp%i_sz)
435    else
436      ! In Rozzi's method the lim q+G --> 0 is finite.
437      vcp%i_sz = vcoul(1,1)
438    end if
439 
440    call vcp%print(unit=ab_out)
441 
442  case ('SLAB')
443    call surface_setup(cryst, vcp%vcutgeo, vcp%alpha, vcp%rcut, vcp%pdir, opt_slab)
444 
445    do iq_ibz=1,nqibz
446      call cutoff_slab(qibz(:,iq_ibz), ng, gvec, cryst%gprimd, vcp%rcut, &
447                       vcp%boxcenter, vcp%pdir, vcp%alpha, vcoul(:,iq_ibz), opt_slab)
448    end do
449 
450    ! q-points for optical limit.
451    do iqlwl=1,nqlwl
452      call cutoff_slab(qlwl(:,iq_ibz), ng, gvec, cryst%gprimd, vcp%rcut, &
453                       vcp%boxcenter, vcp%pdir, vcp%alpha, vcoul_lwl(:,iqlwl), opt_slab)
454    end do
455 
456    ! If Beigi, treat the limit q --> 0.
457    if (opt_slab == 1) then
458      ! Integrate numerically in the plane close to 0
459      call beigi_surface_limit(opt_slab, cryst, nqibz, nkbz, vcp%rcut, vcp%alpha, &
460                               vcp%boxcenter, vcp%pdir, vcp%i_sz)
461    else
462      ! In Rozzi's method the lim q+G --> 0 is finite.
463      vcp%i_sz=vcoul(1,1)
464    end if
465 
466  case ('CRYSTAL', 'AUXILIARY_FUNCTION', "AUX_GB")
467    do iq_ibz=1,nqibz
468      call cmod_qpg(nqibz, iq_ibz, qibz, ng, gvec, cryst%gprimd, vcoul(:,iq_ibz))
469 
470      if (iq_ibz == 1) then
471        ! The singularity is treated using vcoul_lwl.
472        vcoul(1, iq_ibz) = zero
473        vcoul(2:,iq_ibz) = four_pi / vcoul(2:,iq_ibz)**2
474      else
475        vcoul(:,iq_ibz) = four_pi / vcoul(:,iq_ibz)**2
476      end if
477    end do ! iq_ibz
478 
479    ! q-points for optical limit.
480    do iqlwl=1,nqlwl
481      call cmod_qpg(nqlwl, iqlwl, qlwl, ng, gvec, cryst%gprimd, vcoul_lwl(:,iqlwl))
482    end do
483    vcoul_lwl = four_pi/vcoul_lwl**2
484 
485    ! Treatment of 1/q^2 singularity
486 
487    if (vcp%mode == "CRYSTAL") then
488      ! Analytic integration of 4pi/q^2 over the volume element:
489      ! $4pi/V \int_V d^3q 1/q^2 =4pi bz_geometric_factor V^(-2/3)$
490      ! i_sz=4*pi*bz_geometry_factor*q0_vol**(-two_thirds) where q0_vol= V_BZ/N_k
491      ! bz_geometry_factor: sphere=7.79, fcc=7.44, sc=6.188, bcc=6.946, wz=5.255 (see gwa.pdf, appendix A.4)
492      q0_vol = (two_pi) **3 / (nkbz*cryst%ucvol); bz_geometry_factor=zero
493      vcp%i_sz = four_pi*7.44*q0_vol**(-two_thirds)
494 
495    else if (vcp%mode == "AUXILIARY_FUNCTION") then
496      ! Numerical integration of the exact-exchange divergence through the
497      ! auxiliary function of Carrier et al. PRB 75, 205126 (2007) [[cite:Carrier2007]].
498      vcp%i_sz = carrier_isz(cryst, nqbz, qbz, rcut, comm)
499 
500    else if (vcp%mode == "AUX_GB") then
501      ! We use the auxiliary function of a Gygi-Baldereschi variant [[cite:Gigy1986]]
502      vcp%i_sz = gygi_baldereschi_isz(cryst, nqbz, qbz, ecut, ng, gvec)
503 
504    else
505      ABI_ERROR(sjoin("Need treatment of 1/q^2 singularity! for mode", vcp%mode))
506    end if
507 
508  case ('ERF')
509    ! Modified long-range only Coulomb interaction thanks to the error function:
510    ! * Vc = erf(r/rcut)/r
511    ! * The singularity is treated using vcoul_lwl.
512    do iq_ibz=1,nqibz
513      call cmod_qpg(nqibz, iq_ibz, qibz, ng, gvec, cryst%gprimd, vcoul(:,iq_ibz))
514 
515      ! The Fourier transform of the error function reads
516      if (iq_ibz == 1) then
517        vcoul(1, iq_ibz) = zero
518        vcoul(2:,iq_ibz) = four_pi/(vcoul(2:,iq_ibz)**2) *  EXP( -0.25d0 * (vcp%rcut*vcoul(2:,iq_ibz))**2 )
519      else
520        vcoul(:,iq_ibz)  = four_pi/(vcoul(:, iq_ibz)**2) *  EXP( -0.25d0 * (vcp%rcut*vcoul(: ,iq_ibz))**2 )
521      end if
522    end do
523 
524    ! q-points for optical limit.
525    do iqlwl=1,nqlwl
526      call cmod_qpg(nqlwl, iqlwl, qlwl, ng, gvec, cryst%gprimd, vcoul_lwl(:,iqlwl))
527    end do
528    vcoul_lwl = four_pi/(vcoul_lwl**2) *  EXP( -0.25d0 * (vcp%rcut*vcoul_lwl)**2 )
529 
530    ! === Treat 1/q^2 singularity ===
531    ! * We use the auxiliary function from PRB 75, 205126 (2007) [[cite:Carrier2007]]
532    vcp%i_sz = carrier_isz(cryst, nqbz, qbz, rcut, comm)
533 
534  case ('ERFC')
535    ! * Use a modified short-range only Coulomb interaction thanks to the complementary error function:
536    !   $ V_c = [1-erf(r/r_{cut})]/r $
537    ! * The Fourier transform of the error function reads
538    !   vcoul=four_pi/(vcoul**2) * ( 1.d0 - exp( -0.25d0 * (vcp%rcut*vcoul)**2 ) )
539    do iq_ibz=1,nqibz
540      call cmod_qpg(nqibz, iq_ibz, qibz,ng, gvec, cryst%gprimd, vcoul(:,iq_ibz))
541 
542      if (iq_ibz == 1) then
543        vcoul(1 ,iq_ibz) = zero
544        vcoul(2:,iq_ibz) = four_pi/(vcoul(2:,iq_ibz)**2) * ( one - EXP( -0.25d0 * (vcp%rcut*vcoul(2:,iq_ibz))**2 ) )
545      else
546        vcoul(:, iq_ibz) = four_pi/(vcoul(:, iq_ibz)**2) * ( one - EXP( -0.25d0 * (vcp%rcut*vcoul(:, iq_ibz))**2 ) )
547      end if
548    end do ! iq_ibz
549 
550    ! q-points for optical limit.
551    do iqlwl=1,nqlwl
552      call cmod_qpg(nqlwl, iqlwl, qlwl, ng, gvec, cryst%gprimd, vcoul_lwl(:,iqlwl))
553    end do
554    vcoul_lwl = four_pi/(vcoul_lwl**2) * ( one - EXP( -0.25d0 * (vcp%rcut*vcoul_lwl)**2 ) )
555 
556    ! === Treat 1/q^2 singularity ===
557    ! * There is NO singularity in this case.
558    vcp%i_sz = pi * vcp%rcut**2 ! Final result stored here
559 
560  case default
561    ABI_BUG(sjoin('Unsupported cutoff mode:', vcp%mode))
562  end select
563 
564  call wrtout(std_out, sjoin("vcp%i_sz", ftoa(vcp%i_sz)))
565  vcp%i_sz_resid = vcp%i_sz
566 
567  ! Store final results in complex array as Rozzi's cutoff can give real negative values
568  ABI_MALLOC(vcp%vc_sqrt, (ng, nqibz))
569  ABI_MALLOC(vcp%vc_sqrt_resid, (ng, nqibz))
570  vcp%vc_sqrt = CMPLX(vcoul, zero)
571  vcp%vc_sqrt = SQRT(vcp%vc_sqrt)
572  vcp%vc_sqrt_resid = vcp%vc_sqrt
573  ABI_FREE(vcoul)
574 
575  ABI_MALLOC(vcp%vcqlwl_sqrt, (ng, nqlwl))
576  vcp%vcqlwl_sqrt = CMPLX(vcoul_lwl, zero)
577  vcp%vcqlwl_sqrt = SQRT(vcp%vcqlwl_sqrt)
578  ABI_FREE(vcoul_lwl)
579 
580  call vcp%print(unit=std_out)
581 
582 end subroutine vcoul_init

m_vcoul/vcoul_plot [ Functions ]

[ Top ] [ m_vcoul ] [ Functions ]

NAME

 vcoul_plot

FUNCTION

 Plot vccut(q,G) as a function of |q+G|. Calculate also vc in real space.

INPUTS

OUTPUT

SOURCE

 889 subroutine vcoul_plot(Vcp, Qmesh, Gsph, ng, vc, comm)
 890 
 891 !Arguments ------------------------------------
 892 !scalars
 893  class(vcoul_t),intent(in) :: Vcp
 894  integer,intent(in) :: ng,comm
 895  type(kmesh_t),intent(in) :: Qmesh
 896  type(gsphere_t),intent(in) :: Gsph
 897 !arrays
 898  real(dp),intent(in) :: vc(ng,Qmesh%nibz)
 899 
 900 !Local variables-------------------------------
 901 !scalars
 902  integer,parameter :: master = 0
 903  integer :: icount,idx_Sm1G,ierr,ig,igs,ii,iq_bz,iq_ibz,iqg,ir,isym,itim
 904  integer :: my_start,my_stop,nqbz,nqibz,nr,ntasks,my_rank,unt
 905  real(dp) :: arg,fact,l1,l2,l3,lmax,step,tmp,vcft,vc_bare
 906  character(len=500) :: msg
 907  character(len=fnlen) :: filnam
 908 !arrays
 909  integer,allocatable :: insort(:)
 910  real(dp) :: b1(3),b2(3),b3(3),gmet(3,3),gprimd(3,3),qbz(3),qpgc(3)
 911  real(dp),allocatable :: qpg_mod(:),rr(:,:,:),vcr(:,:),vcr_cut(:,:)
 912 
 913 !************************************************************************
 914 
 915  if (trim(Vcp%mode) /= 'CYLINDER') RETURN
 916 
 917  my_rank = xmpi_comm_rank(comm)
 918 
 919  nqibz=Vcp%nqibz; nqbz=Qmesh%nbz
 920  gmet=Gsph%gmet; gprimd=Gsph%gprimd
 921 
 922  b1(:)=two_pi*gprimd(:,1)
 923  b2(:)=two_pi*gprimd(:,2)
 924  b3(:)=two_pi*gprimd(:,3)
 925 
 926  ! Compare in Fourier space the true Coulomb with the cutted one.
 927  if (my_rank == master) then
 928    ABI_MALLOC(insort, (nqibz * ng))
 929    ABI_MALLOC(qpg_mod, (nqibz * ng))
 930    iqg = 1
 931    do iq_ibz=1,nqibz
 932      do ig=1,ng
 933        qpg_mod(iqg) = normv(Qmesh%ibz(:,iq_ibz) + Gsph%gvec(:,ig), gmet,'g')
 934        insort(iqg) = iqg; iqg = iqg + 1
 935      end do
 936    end do
 937    call sort_dp(nqibz * ng, qpg_mod, insort, tol14)
 938 
 939    filnam='_VCoulFT_'
 940    call isfile(filnam, 'new')
 941    if (open_file(filnam, msg, newunit=unt, status='new', form='formatted') /= 0) then
 942      ABI_ERROR(msg)
 943    end if
 944    write(unt,'(a,i3,a,i6,a)')&
 945     '#   |q+G|       q-point (Tot no.',nqibz,')        Gvec (',ng,')     vc_bare(q,G)    vc_cutoff(q,G) '
 946 
 947    do iqg=1,nqibz*ng
 948      iq_ibz = (insort(iqg) - 1) / ng + 1
 949      ig = (insort(iqg)) - (iq_ibz-1) * ng
 950      vc_bare = zero
 951      if (qpg_mod(iqg) > tol16) vc_bare = four_pi / qpg_mod(iqg) ** 2
 952      write(unt,'(f12.6,2x,3f8.4,2x,3i6,2x,2es14.6)')&
 953        qpg_mod(iqg), Qmesh%ibz(:,iq_ibz), Gsph%gvec(:,ig), vc_bare, vc(ig, iq_ibz)
 954    end do
 955 
 956    close(unt)
 957    ABI_FREE(insort)
 958    ABI_FREE(qpg_mod)
 959  end if ! my_rank==master
 960 
 961  ! Fourier transform back to real space just to check cutoff implementation.
 962  ntasks= nqbz * ng
 963  call xmpi_split_work(ntasks, comm, my_start, my_stop)
 964 
 965  l1 = SQRT(SUM(Vcp%rprimd(:,1)**2))
 966  l2 = SQRT(SUM(Vcp%rprimd(:,2)**2))
 967  l3 = SQRT(SUM(Vcp%rprimd(:,3)**2))
 968 
 969  nr = 50
 970  lmax=MAX(l1,l2,l3) ; step=lmax/(nr-1)
 971  fact = one / (Vcp%ucvol * nqbz)
 972 
 973  ! numb coding
 974  ABI_CALLOC(rr, (3, nr, 3))
 975  do ii=1,3
 976    do ir=1,nr
 977      rr(ii,ir,ii)=(ir-1)*step
 978    end do
 979  end do
 980 
 981  ABI_CALLOC(vcr, (nr, 3))
 982  ABI_CALLOC(vcr_cut, (nr, 3))
 983 
 984  do iq_bz=1,nqbz
 985    call Qmesh%get_BZ_item(iq_bz, qbz, iq_ibz, isym, itim)
 986    if (ABS(qbz(1))<0.01) qbz(1)=zero
 987    if (ABS(qbz(2))<0.01) qbz(2)=zero
 988    if (ABS(qbz(3))<0.01) qbz(3)=zero
 989    igs=1; if (ALL(qbz(:)==zero)) igs=2
 990    do ig=igs,ng
 991      icount=ig+(iq_bz-1)*ng
 992      if (icount < my_start .or. icount > my_stop) CYCLE
 993      idx_Sm1G = Gsph%rottbm1(ig,itim,isym) ! IS{^-1}G
 994      vcft=vc(idx_Sm1G,iq_ibz)
 995      qpgc(:)=qbz(:)+Gsph%gvec(:,ig) ; qpgc(:)=b1(:)*qpgc(1)+b2(:)*qpgc(2)+b3(:)*qpgc(3)
 996      tmp=SQRT(DOT_PRODUCT(qpgc,qpgc)) ; tmp=tmp**2
 997      do ii=1,3
 998        do ir=1,nr
 999          arg=DOT_PRODUCT(rr(:,ir,ii),qpgc)
1000          vcr_cut(ir,ii)=vcr_cut(ir,ii) + vcft*COS(arg)
1001          vcr    (ir,ii)=vcr    (ir,ii) + four_pi/tmp*COS(arg)
1002        end do
1003      end do
1004    end do !ig
1005  end do !iq_ibz
1006 
1007  call xmpi_sum_master(vcr_cut,master,comm,ierr)
1008  call xmpi_sum_master(vcr    ,master,comm,ierr)
1009 
1010  if (my_rank == master) then
1011    filnam='_VCoulR_'
1012    call isfile(filnam, 'new')
1013    if (open_file(filnam,msg,newunit=unt,status='new',form='formatted') /= 0) then
1014      ABI_ERROR(msg)
1015    end if
1016    write(unt,'(a)')'# length  vc_bare(r)   vc_cut(r) '
1017    do ir=1,nr
1018      write(unt,'(7es18.6)')(ir-1)*step,(fact*vcr(ir,ii),fact*vcr_cut(ir,ii),ii=1,3)
1019    end do
1020    close(unt)
1021  end if
1022 
1023  ABI_FREE(rr)
1024  ABI_FREE(vcr)
1025  ABI_FREE(vcr_cut)
1026 
1027 end subroutine vcoul_plot

m_vcoul/vcoul_print [ Functions ]

[ Top ] [ m_vcoul ] [ Functions ]

NAME

 vcoul_print

FUNCTION

  Print the content of a Coulomb datatype.

INPUTS

  Vcp<vcoul_t>=The datatype whose content has to be printed.
  [unit]=The unit number for output
  [prtvol]=Verbosity level
  [mode_paral]=Either "COLL" or "PERS".

SOURCE

1047 subroutine vcoul_print(Vcp, unit, prtvol, mode_paral)
1048 
1049 !Arguments ------------------------------------
1050 !scalars
1051  class(vcoul_t),intent(in) :: Vcp
1052  integer,intent(in),optional :: prtvol,unit
1053  character(len=4),intent(in),optional :: mode_paral
1054 
1055 !Local variables-------------------------------
1056 !scalars
1057  integer :: ii,my_unt,my_prtvol,iqlwl
1058  character(len=4) :: my_mode
1059  character(len=500) :: msg
1060 
1061 ! *************************************************************************
1062 
1063  my_unt    =std_out; if (PRESENT(unit      )) my_unt   =unit
1064  my_mode   ='COLL' ; if (PRESENT(mode_paral)) my_mode  =mode_paral
1065  my_prtvol=0       ; if (PRESENT(prtvol    )) my_prtvol=prtvol
1066 
1067  select case (Vcp%mode)
1068 
1069  case ('MINIBZ')
1070    write(msg,'(3a)')ch10,' vcoul_init : cutoff-mode = ',trim(Vcp%mode)
1071    call wrtout(my_unt,msg,my_mode)
1072 
1073  case ('MINIBZ-ERF')
1074    write(msg,'(3a)')ch10,' vcoul_init : cutoff-mode = ',trim(Vcp%mode)
1075    call wrtout(my_unt,msg,my_mode)
1076    write(msg,'(5a,f10.4,3a,f10.2,3a,3f10.5,2a)')ch10,&
1077      ' === Error function cutoff === ',ch10,ch10,&
1078      '  Cutoff radius ......... ',Vcp%rcut,' [Bohr] ',ch10
1079    call wrtout(my_unt,msg,my_mode)
1080 
1081  case ('MINIBZ-ERFC')
1082    write(msg,'(3a)')ch10,' vcoul_init : cutoff-mode = ',trim(Vcp%mode)
1083    call wrtout(my_unt,msg,my_mode)
1084    write(msg,'(5a,f10.4,3a,f10.2,3a,3f10.5,2a)')ch10,&
1085      ' === Complement Error function cutoff === ',ch10,ch10,&
1086      '  Cutoff radius ......... ',Vcp%rcut,' [Bohr] ',ch10
1087    call wrtout(my_unt,msg,my_mode)
1088 
1089  case ('SPHERE')
1090    write(msg,'(5a,f10.4,3a,f10.2,3a,3f10.5,2a)')ch10,&
1091     ' === Spherical cutoff === ',ch10,ch10,&
1092     '  Cutoff radius ......... ',Vcp%rcut,' [Bohr] ',ch10,&
1093     '  Volume of the sphere .. ',four_pi/three*Vcp%rcut**3,' [Bohr^3] '
1094     !FB: This has no meaning here! &  '  Sphere centered at .... ',Vcp%boxcenter,' (r.l.u) ',ch10
1095     !MG It might be useful if the system is not centered on the origin because in this case the
1096     !   matrix elements of the Coulomb have to be multiplied by a phase depending on boxcenter.
1097     !   I still have to decide if it is useful to code this possibility and which variable use to
1098     !   define the center (boxcenter is used in the tddft part).
1099    call wrtout(my_unt,msg,my_mode)
1100 
1101  case ('CYLINDER')
1102    ii=imin_loc(ABS(Vcp%pdir-1))
1103    write(msg,'(5a,f10.4,3a,i2,2a,3f10.2,a)')ch10,&
1104      ' === Cylindrical cutoff === ',ch10,ch10,&
1105      '  Cutoff radius ............... ',Vcp%rcut,' [Bohr] ',ch10,&
1106      '  Axis parallel to direction... ',ii,ch10,&
1107      '  Passing through point ....... ',Vcp%boxcenter,' (r.l.u) '
1108    call wrtout(my_unt,msg,my_mode)
1109 
1110    write(msg,'(2a)')'  Infinite length  ....... ',ch10
1111    if (Vcp%hcyl/=zero) write(msg,'(a,f8.5,2a)')'  Finite length of ....... ',Vcp%hcyl,' [Bohr] ',ch10
1112    call wrtout(my_unt,msg,my_mode)
1113 
1114  CASE ('SLAB')
1115    write(msg,'(5a,f10.4,3a,3f10.2,2a)')ch10,&
1116      ' === Surface cutoff === ',ch10,ch10,&
1117      '  Cutoff radius .................... ',Vcp%rcut,' [Bohr] ',ch10,&
1118      '  Central plane passing through .... ',Vcp%boxcenter,' (r.l.u) ',ch10
1119    call wrtout(my_unt,msg,my_mode)
1120    !write(msg,'(a)')'  Infinite length  .......'
1121    !if (Vcp%hcyl/=zero) write(msg,'(a,f8.5,a)')'  Finite length of .......',Vcp%hcyl,' [Bohr] '
1122    !call wrtout(my_unt,msg,my_mode)
1123 
1124  case ('AUXILIARY_FUNCTION')
1125    write(msg,'(3a)')ch10,' vcoul_init : cutoff-mode = ',trim(Vcp%mode)
1126    call wrtout(my_unt,msg,my_mode)
1127 
1128  case ('AUX_GB')
1129    write(msg,'(3a)')ch10,' vcoul_init : cutoff-mode = ',trim(Vcp%mode)
1130    call wrtout(my_unt,msg,my_mode)
1131 
1132  case ('CRYSTAL')
1133    write(msg,'(3a)')ch10,' vcoul_init : cutoff-mode = ',trim(Vcp%mode)
1134    call wrtout(my_unt,msg,my_mode)
1135 
1136  case ('ERF')
1137    write(msg,'(5a,f10.4,3a,f10.2,3a,3f10.5,2a)')ch10,&
1138      ' === Error function cutoff === ',ch10,ch10,&
1139      '  Cutoff radius ......... ',Vcp%rcut,' [Bohr] ',ch10
1140    call wrtout(my_unt,msg,my_mode)
1141 
1142  case ('ERFC')
1143    write(msg,'(5a,f10.4,3a,f10.2,3a,3f10.5,2a)')ch10,&
1144      ' === Complement Error function cutoff === ',ch10,ch10,&
1145      '  Cutoff radius ......... ',Vcp%rcut,' [Bohr] ',ch10
1146    call wrtout(my_unt,msg,my_mode)
1147 
1148  case default
1149    ABI_BUG(sjoin('Unknown cutoff mode: ', Vcp%mode))
1150  end select
1151 
1152  if (Vcp%nqlwl > 0) then
1153    write(msg,'(a,i3)')" q-points for optical limit: ",Vcp%nqlwl
1154    call wrtout(my_unt,msg,my_mode)
1155    do iqlwl=1,Vcp%nqlwl
1156      write(msg,'(1x,i5,a,2x,3f12.6)') iqlwl,')',Vcp%qlwl(:,iqlwl)
1157      call wrtout(my_unt,msg,my_mode)
1158    end do
1159  end if
1160 
1161  !TODO add additional information
1162 
1163 end subroutine vcoul_print

m_vcoul/vcoul_t [ Types ]

[ Top ] [ m_vcoul ] [ Types ]

NAME

  vcoul_t

FUNCTION

  This data type contains the square root of the Fourier components of the Coulomb interaction
  calculated taking into account a possible cutoff. It also stores info on the particular geometry
  used for the cutoff as well as quantities required to deal with the Coulomb divergence for q --> 0.

SOURCE

 74 type,public :: vcoul_t
 75 
 76   integer :: ng = -1
 77    ! Number of G-vectors
 78 
 79   integer :: nqibz = -1
 80    ! Number of irreducible q-points
 81 
 82   integer :: nqlwl = -1
 83    ! Number of small q-points around Gamma
 84 
 85   real(dp) :: alpha(3) = -one
 86    ! Length of the finite slab
 87 
 88   real(dp) :: rcut = -one
 89    ! Cutoff radius.
 90 
 91   real(dp) :: i_sz = huge(one)
 92    ! Value of the integration of the Coulomb singularity 4\pi/V_BZ \int_BZ d^3q 1/q^2
 93 
 94   real(dp) :: i_sz_resid = huge(one)
 95    ! Residual difference between the i_sz in the sigma self-energy for exchange,
 96    ! and the i_sz already present in the generalized Kohn-Sham eigenenergies
 97    ! Initialized to the same value as i_sz
 98 
 99   real(dp) :: hcyl = -one
100    ! Length of the finite cylinder along the periodic dimension
101 
102   real(dp) :: ucvol = -one
103     ! Volume of the unit cell
104 
105   character(len=50) :: mode
106    ! String defining the cutoff mode, possible values are: sphere,cylinder,slab,crystal
107 
108   integer :: pdir(3)
109    ! 1 if the system is periodic along this direction
110 
111   real(dp) :: boxcenter(3) = -1
112    ! 1 if the point in inside the cutoff region 0 otherwise
113    ! Reduced coordinates of the center of the box (input variable)
114 
115   real(dp) :: vcutgeo(3) = huge(one)
116    ! For each reduced direction gives the length of the finite system
117    ! 0 if the system is infinite along this direction.
118    ! negative values indicate that a finite size has to be used.
119 
120   real(dp) :: rprimd(3,3) = zero
121    ! Lattice vectors in real space.
122 
123   real(dp),allocatable :: qibz(:,:)
124    ! (3, nqibz)
125    ! q-points in the IBZ.
126 
127   real(dp),allocatable :: qlwl(:,:)
128    ! (3, nqlwl)
129    ! q-points for the treatment of the Coulomb singularity.
130 
131   complex(gwpc),allocatable :: vc_sqrt(:,:)
132    ! (ng, nqibz)
133    ! Square root of the Coulomb interaction in reciprocal space.
134    ! complex-valued to allow for a possible cutoff (Rozzi's method)
135 
136   complex(gwpc),allocatable :: vcqlwl_sqrt(:,:)
137    ! (ng, nqlwl)
138    ! Square root of the Coulomb term calculated for small q-points
139 
140   complex(gwpc),allocatable :: vc_sqrt_resid(:,:)
141    ! (ng, nqibz)
142    ! Square root of the residual difference between the Coulomb interaction in the sigma self-energy for exchange,
143    ! and the Coulomb interaction already present in the generalized Kohn-Sham eigenenergies (when they come from an hybrid)
144    ! Given in reciprocal space. At the call to vcoul_init, it is simply initialized at the value of vc_sqrt(:,:),
145    ! and only later modified. A cutoff might be applied.
146 
147 contains
148    procedure :: init => vcoul_init    ! Main creation method.
149    procedure :: plot => vcoul_plot    ! Plot vc in real and reciprocal space.
150    procedure :: print => vcoul_print  ! Report info on the object.
151    procedure :: free => vcoul_free    ! Free memory
152 
153 end type vcoul_t