TABLE OF CONTENTS
- ABINIT/m_vcoul
- m_vcoul/adapt_nmc
- m_vcoul/beigi_cylinder_limit
- m_vcoul/beigi_surface_limit
- m_vcoul/carrier_isz
- m_vcoul/cylinder_setup
- m_vcoul/gw_icutcoul_to_mode
- m_vcoul/gygi_baldereschi_isz
- m_vcoul/integratefaux
- m_vcoul/mc_free
- m_vcoul/mc_init
- m_vcoul/mc_integrate
- m_vcoul/mc_t
- m_vcoul/surface_setup
- m_vcoul/vcgen_free
- m_vcoul/vcgen_get_vc_sqrt
- m_vcoul/vcgen_init
- m_vcoul/vcgen_t
- m_vcoul/vcoul_free
- m_vcoul/vcoul_init
- m_vcoul/vcoul_plot
- m_vcoul/vcoul_print
- m_vcoul/vcoul_t
ABINIT/m_vcoul [ 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 ]
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 ]
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 ]
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