TABLE OF CONTENTS


ABINIT/m_kxc [ Modules ]

[ Top ] [ Modules ]

NAME

 m_kxc

FUNCTION

 Helper functions to compute the XC kernel in reciprocal space.
 WARNING: At present (10/01/14) these routines are not tested
 since the ACFD code has been disabled.

COPYRIGHT

  Copyright (C) 1998-2018 ABINIT group (DCA, MF, XG, GMR, LSI, YMN, Rhaltaf, MS)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

 NOTES:
  libxc_functionals.F90 uses a global structure (funcs) to store the XC parameters.
  This structure is initialized in driver with the value of ixc specified by the user in the input file.
  In order to change the value of ixc at run-time, we have to reinitialize the global structure
  with the new value of ixc before computing XC quantities.
  Moreover one has to reinstate the old functional before returning so that the other routines
  will continue to used the previous ixc. This task can be accomplished with the following pseudocode

   ! Reinitialize the libxc module with the overriden values
   if (old_ixc<0)  call libxc_functionals_end()
   if (new_ixc<0) call libxc_functionals_init(new_ixc,nspden)
   ! Compute XC stuff here.
   ! Revert libxc module to the original settings
   if (new_ixc<0) call libxc_functionals_end()
   if (old_ixc<0) call libxc_functionals_init(old_ixc,nspden)

PARENTS

SOURCE

37 #if defined HAVE_CONFIG_H
38 #include "config.h"
39 #endif
40 
41 #include "abi_common.h"
42 
43 MODULE m_kxc
44 
45  use defs_basis
46  use m_abicore
47  use m_errors
48  use m_xmpi
49  use m_crystal
50  use m_distribfft
51  use m_xcdata
52  use libxc_functionals
53 
54  use defs_abitypes,   only : MPI_type, dataset_type
55  use m_dtset,         only : dtset_copy, dtset_free
56  use m_io_tools,      only : open_file
57  use m_pptools,       only : printxsf
58  use m_numeric_tools, only : hermitianize
59  use m_fft_mesh,      only : g2ifft
60  use m_fft,           only : fourdp_6d, fourdp
61  use m_mpinfo,        only : destroy_mpi_enreg, initmpi_seq
62  use m_spacepar,      only : hartre
63  use m_rhotoxc,       only : rhotoxc
64  use m_dfpt_mkvxc,    only : dfpt_mkvxc
65 
66  implicit none
67 
68  private

m_kxc/kxc_ADA [ Functions ]

[ Top ] [ m_kxc ] [ Functions ]

NAME

 kxc_ADA

FUNCTION

 Calculate exchange-correlation kernel in reciprocal space

INPUTS

 Dtset <type(dataset_type)>=all input variables in this dataset
 Cryst<crystal_t>=Info on the unit cell.
 ixc = choice for the exchange-correlation potential.
 ngfft(18)=contain all needed information about 3D FFT,
  see ~abinit/doc/variables/vargs.htm#ngfft
 nfft = total number of points on the FFT grid.
 rhor(nfft,nspden) = the charge density on the FFT grid.
  (total in first half and spin-up in second half if nsppol=2)
 npw: the size of kernel matrix
 dim_kxcg=dimension of the kernel.
 comm=MPI communicator.
 [dbg_mode]=Set it to .TRUE. to switch on the debug mode.

OUTPUT

  kxcg(nfft,dim_kxcg) = the exchange-correlation potential on the FFT grid.
  warning: the kernel is not divided by unit cell volume

NOTES

  No xc quadrature
  No nl core correction

PARENTS

      screening,sigma

CHILDREN

      destroy_mpi_enreg,fourdp,fourdp_6d,hartre,initmpi_seq
      libxc_functionals_end,libxc_functionals_init,printxsf,rhotoxc,wrtout
      xcdata_init

SOURCE

1331 subroutine kxc_ADA(Dtset,Cryst,ixc,ngfft,nfft,nspden,rhor,&
1332 &                  npw,nqibz,qibz,fxc_ADA,gvec,comm,kappa_init,dbg_mode)
1333 
1334 
1335 !This section has been created automatically by the script Abilint (TD).
1336 !Do not modify the following lines by hand.
1337 #undef ABI_FUNC
1338 #define ABI_FUNC 'kxc_ADA'
1339 !End of the abilint section
1340 
1341  implicit none
1342 
1343 !Arguments ------------------------------------
1344 !scalars
1345  integer,intent(in) :: ixc,nfft,nspden,npw,comm
1346  real(dp),intent(in),optional :: kappa_init
1347  logical,optional,intent(in) :: dbg_mode
1348  type(crystal_t),intent(in) :: Cryst
1349  type(dataset_type),intent(in) :: Dtset
1350 !arrays
1351  integer,intent(in) :: gvec(3,npw),ngfft(18)
1352  integer,intent(in) :: nqibz
1353  real(dp),intent(in) :: rhor(nfft,nspden)
1354  real(dp),intent(in) :: qibz(3,nqibz)
1355  complex(gwpc),intent(out) :: fxc_ADA(npw,npw,nqibz)
1356 
1357 !Local variables ------------------------------
1358 !scalars
1359  integer,parameter :: paral_kgb0=0
1360  integer :: i1,i2,i3,ig,igp,ir,irp,n3xccc,ngfft1,ngfft2,izero !,isp
1361  integer :: ngfft3,nkxc,option,ikxc,ierr,nproc
1362  integer :: nk3xc,igrid,iqbz,my_rank,master,unt_dmp,gmgp_idx
1363  real(dp) :: enxc,gsqcut,ucvol !,rs,Kx,Kc
1364  real(dp) :: vxcavg,kappa,abs_qpg_sq,abs_qpgp_sq
1365  real(dp) :: difx,dify,difz,inv_kappa_sq
1366  character(len=500) :: msg,fname
1367  type(MPI_type) :: MPI_enreg_seq
1368  type(xcdata_type) :: xcdata
1369 !arrays
1370  real(dp) :: qpg(3),qpgp(3),qphon(3),strsxc(6),q_point(3),dum(0)
1371  real(dp),allocatable :: kxcr(:,:)
1372  real(dp),allocatable :: rhog(:,:),vhartr(:),vxclda(:,:)
1373  real(dp),allocatable :: xccc3d(:),my_rhor(:,:)
1374  real(dp),allocatable :: my_kxcg(:,:)
1375  real(dp),allocatable :: rhotilder(:,:)
1376  complex(gwpc),allocatable :: my_fxc_ADA_ggpq(:,:,:)
1377  complex(gwpc),allocatable :: FT_fxc_ADA_ggpq(:,:,:),dummy(:,:)
1378  real(dp),allocatable :: rvec(:,:),my_fxc_ADA_rrp(:,:)
1379  real(dp) :: rmrp(3),abs_rmrp
1380  integer :: n1,n2,n3,ig_idx_fft(npw)
1381 
1382 ! ************************************************************************
1383 
1384  ABI_CHECK(Dtset%nsppol==1,'nsppol/=1 not coded')
1385  ABI_CHECK(nspden==1,'nsppol/=1 not coded')
1386  ABI_CHECK(nfft==PRODUCT(ngfft(1:3)),"mismatch in nfftot")
1387 
1388 !Fake MPI_type for the sequential part.
1389  call initmpi_seq(MPI_enreg_seq)
1390 
1391  my_rank = xmpi_comm_rank(comm)
1392  nproc   = xmpi_comm_size(comm)
1393  master=0
1394 
1395  write(msg,'(a,i3)') ' kxc_ADA: calculating exchange-correlation kernel using ixc = ',ixc
1396  call wrtout(std_out,msg,'COLL')
1397  call wrtout(std_out,' kxc_ADA: using smeared density','COLL')
1398 
1399  if (.not.present(kappa_init)) then
1400    kappa = 2.1_dp
1401  else
1402    kappa = kappa_init
1403  end if
1404  write(msg,'(a,F10.3)') ' kxc_ADA: inverse smearing length, kappa = ',kappa
1405  call wrtout(std_out,msg,'COLL')
1406  inv_kappa_sq = one/(kappa*kappa)
1407 
1408  call xcdata_init(xcdata,dtset=dtset,intxc=0,ixc=ixc,nspden=nspden)
1409 
1410  if (ALL(xcdata%xclevel/=(/1,2/))) then
1411    write(msg,'(a,i0)')"Unsupported xclevel = ",xcdata%xclevel
1412    MSG_ERROR(msg)
1413  end if
1414 
1415  ngfft1=ngfft(1)
1416  ngfft2=ngfft(2)
1417  ngfft3=ngfft(3)
1418 
1419  if (ixc>=1.and.ixc<11) then      ! LDA case
1420    nkxc= 2*min(Dtset%nspden,2)-1  ! 1 or 3
1421  else                             ! GGA case
1422    nkxc=12*min(Dtset%nspden,2)-5  ! 7 or 19
1423    ABI_CHECK(dtset%xclevel==2,"Functional should be GGA")
1424    MSG_ERROR("GGA functional not implemented for ADA vertex")
1425  end if
1426 
1427  ABI_MALLOC(kxcr,(nfft,nkxc))
1428 
1429 !gsqcut and rhog are zeroed because they are not used by rhotoxc if 1<=ixc<=16 and option=0
1430  gsqcut=zero
1431 
1432  ABI_MALLOC(rhog,(2,nfft))
1433  ABI_MALLOC(vhartr,(nfft))
1434  rhog(:,:)=zero
1435 !MG FIXME this is the 3D core electron density for XC core correction (bohr^-3)
1436 !should implement the non linear core correction
1437  n3xccc=0
1438  ABI_MALLOC(xccc3d,(n3xccc))
1439  ABI_MALLOC(vxclda,(nfft,nspden))
1440 
1441  option=2 ! 2 for Hxc and kxcr (no paramagnetic part if nspden=1)
1442  qphon(:)=zero
1443 
1444 !to be adjusted for the call to rhotoxc
1445  nk3xc=1
1446 
1447 !Compute the kernel.
1448  izero=0
1449 
1450 !DEBUG print density
1451  if (present(dbg_mode)) then
1452    if (dbg_mode.and.my_rank==master) then
1453      fname = 'xc_ADA_den.xsf'
1454      if (open_file(fname,msg,newunit=unt_dmp,status='unknown',form='formatted') /= 0) then
1455        MSG_ERROR(msg)
1456      end if
1457      call printxsf(ngfft1,ngfft2,ngfft3,rhor(:,1),Cryst%rprimd,(/zero,zero,zero/),&
1458 &     Cryst%natom,Cryst%ntypat,Cryst%typat,Cryst%xcart,Cryst%znucl,unt_dmp,0)
1459      close(unt_dmp)
1460    end if
1461  end if
1462 !DEBUG
1463 
1464 !Calculate the smeared density
1465  ABI_MALLOC(my_rhor,(nfft,nspden))
1466  ABI_MALLOC(rhotilder,(nfft,nspden))
1467  ucvol = Cryst%ucvol
1468  my_rhor = rhor
1469 !do isp = 1,nsppol
1470 !call calc_smeared_density(my_rhor(:,isp),1,rhotilder(:,isp),nfft,ngfft,npw,&
1471 !&   gvec,Cryst%gprimd,Cryst%ucvol,MPI_enreg_seq,paral_kgb0,kappa_in=kappa)
1472 !my_rhor(:,isp) = rhotilder(:,isp)
1473 !end do
1474 
1475 !DEBUG print smeared density
1476  if (present(dbg_mode)) then
1477    if (dbg_mode.and.my_rank==master) then
1478      fname = 'xc_ADA_smeared_den.xsf'
1479      if (open_file(fname,msg,newunit=unt_dmp,status='unknown',form='formatted') /= 0) then
1480        MSG_ERROR(msg)
1481      end if
1482      call printxsf(ngfft1,ngfft2,ngfft3,my_rhor(:,1),Cryst%rprimd,(/zero,zero,zero/),&
1483 &     Cryst%natom,Cryst%ntypat,Cryst%typat,Cryst%xcart,Cryst%znucl,unt_dmp,0)
1484      close(unt_dmp)
1485    end if
1486  end if
1487 !DEBUG
1488 
1489  ! Reinitialize the libxc module with the overriden values
1490  if (dtset%ixc<0) then
1491    call libxc_functionals_end()
1492  end if
1493  if (ixc<0) then
1494    call libxc_functionals_init(ixc,Dtset%nspden)
1495  end if
1496 
1497  call hartre(1,gsqcut,izero,MPI_enreg_seq,nfft,ngfft,dtset%paral_kgb,rhog,Cryst%rprimd,vhartr)
1498  call rhotoxc(enxc,kxcr,MPI_enreg_seq,nfft,ngfft,&
1499 & dum,0,dum,0,nkxc,nk3xc,dtset%usepawu==4.or.dtset%usepawu==14,&
1500 & n3xccc,option,dtset%paral_kgb,my_rhor,Cryst%rprimd,&
1501 & strsxc,1,vxclda,vxcavg,xccc3d,xcdata,vhartr=vhartr)
1502 
1503 !Check for extreme (NaN) values
1504 !do ir=1,nfft
1505 !if (isnan(kxcr(ir,1))) kxcr(ir,1) = HUGE(kxcr(ir,1))
1506 !end do
1507 
1508 !DEBUG test with direct way of calculating Kxc
1509 !do i1=1,nfft
1510 !rs = (three/(four_pi*my_rhor(i1,1)))**third
1511 !Kx = 16._dp/27._dp*0.3141592653589793e1_dp*(rs**2)*(-0.4581652_dp)
1512 !
1513 !Kc =  -0.4e1_dp / 0.9e1_dp * 0.3141592654e1_dp * rs ** 4 &
1514 !* (0.207271333333333333333333333333e-1_dp * &
1515 !(-0.177442658629204480000000e3_dp * rs - 0.17565190511219200000000e2_dp &
1516 !* sqrt(rs) - 0.1332650665120000e2_dp * rs ** 2 &
1517 !- 0.51031691247948928000000e2_dp * rs ** (0.3e1_dp / 0.2e1_dp)) &
1518 !* rs ** (-0.3e1_dp / 0.2e1_dp) / (rs + 0.37274400e1_dp * sqrt(rs) &
1519 !+ 0.129352000e2_dp) ** 2 / (-sqrt(rs) - 0.1049800_dp) &
1520 !+ 0.518178333333333333333333333333e-2_dp * rs ** (-0.3e1_dp / 0.2e1_dp) &
1521 !* (0.617071835390850041282140897280e3_dp * sqrt(rs) &
1522 !+ 0.659369347307557491857191871552e5_dp * rs ** 2 + &
1523 !0.700403648491298930017835369562e5_dp * rs ** (0.3e1_dp / 0.2e1_dp) &
1524 !+ 0.398437532951539263722720308167e5_dp * rs ** (0.5e1_dp / 0.2e1_dp) &
1525 !+ 0.368852071032531998953472000000e4_dp * rs ** (0.7e1_dp / 0.2e1_dp) &
1526 !+ 0.5330602660480000e2_dp * rs ** (0.9e1_dp / 0.2e1_dp) &
1527 !+ 0.143783940386264738593799346176e5_dp * rs ** 3 &
1528 !+ 0.124672564145568409213848436081e5_dp * rs &
1529 !+ 0.557398029956167136000000e3_dp * rs ** 4) &
1530 !/ (rs + 0.37274400e1_dp * sqrt(rs) + 0.129352000e2_dp) ** 4 &
1531 !/ (-sqrt(rs) - 0.1049800_dp) ** 2)
1532 !kxcr(i1,1) = Kx + Kc
1533 !end do
1534 !END DEBUG
1535 
1536 !DEBUG print Kxc
1537  if (present(dbg_mode)) then
1538    if (dbg_mode.and.my_rank==master) then
1539      fname = 'xc_ADA_Kxc.xsf'
1540      if (open_file(fname,msg,newunit=unt_dmp,status='unknown',form='formatted') /= 0) then
1541        MSG_ERROR(msg)
1542      end if
1543      call printxsf(ngfft1,ngfft2,ngfft3,kxcr(:,1),Cryst%rprimd,(/zero,zero,zero/),&
1544 &     Cryst%natom,Cryst%ntypat,Cryst%typat,Cryst%xcart,Cryst%znucl,unt_dmp,0)
1545      close(unt_dmp)
1546    end if
1547  end if
1548 !DEBUG
1549 
1550  ABI_FREE(xccc3d)
1551  ABI_FREE(vxclda)
1552  ABI_FREE(vhartr)
1553 
1554  ABI_MALLOC(my_kxcg,(2,nfft))
1555 
1556  do ikxc=1,nkxc
1557    call fourdp(1,my_kxcg,kxcr(:,ikxc),-1,MPI_enreg_seq,nfft,ngfft,paral_kgb0,0)
1558 !  kxcg(:,ikxc)=CMPLX(my_kxcg(1,:),my_kxcg(2,:))
1559  end do
1560 !TODO Check symmetry of kxcg
1561 
1562 !set up ADA vertex
1563  ABI_MALLOC(my_fxc_ADA_ggpq,(npw,npw,nqibz))
1564  my_fxc_ADA_ggpq = czero
1565 !Calculate f_xc(R,R')=(kappa^2/2)K_xc[\tilde{n}](G-G')
1566 !x(1/(kappa^2+|q+G|^2) + 1/(kappa^2+|q+G'|^2
1567 !First get G vectors and indices
1568 
1569  ierr=0
1570  do iqbz=1,nqibz
1571    q_point(:) = qibz(:,iqbz)
1572 !
1573    do ig=1,npw
1574      do igp=1,npw
1575 !      Calculate |q+G| and |q+G'|
1576        qpg(:) = two_pi*MATMUL(Cryst%gprimd,q_point(:)+gvec(:,ig))
1577        qpgp(:) = two_pi*MATMUL(Cryst%gprimd,q_point(:)+gvec(:,igp))
1578        abs_qpg_sq = 1.0_dp/(1.0_dp+dot_product(qpg,qpg)*inv_kappa_sq)
1579        abs_qpgp_sq = 1.0_dp/(1.0_dp+dot_product(qpgp,qpgp)*inv_kappa_sq)
1580 
1581        gmgp_idx = g2ifft(gvec(:,ig)-gvec(:,igp),ngfft)
1582        if (gmgp_idx>0) then
1583          my_fxc_ADA_ggpq(ig,igp,iqbz) = half*CMPLX(my_kxcg(1,gmgp_idx), my_kxcg(2,gmgp_idx))*(abs_qpg_sq+abs_qpgp_sq)
1584        else
1585          ierr=ierr+1
1586          my_fxc_ADA_ggpq(ig,igp,iqbz) = czero
1587        end if
1588      end do
1589    end do
1590    if (ierr/=0) then
1591      write(msg,'(a,i4,3a)')&
1592 &     ' Found ',ierr,' G1-G2 vectors falling outside the FFT box. ',ch10,&
1593 &     ' Enlarge the FFT mesh to get rid of this problem. '
1594      MSG_WARNING(msg)
1595    end if
1596  end do
1597 
1598  fxc_ADA = my_fxc_ADA_ggpq
1599 
1600 !do iqbz=1,nqibz
1601 !call hermitianize(my_fxc_ADA_ggpq(:,:,iqbz),"All")
1602 !end do
1603 
1604 
1605 !DEBUG check symmetry
1606  if (.FALSE.) then
1607 !  do iqbz=1,nkptgw
1608 !  do ig=1,npw
1609 !  do igp=ig,npw
1610 !  if (ABS(REAL(fxc_ADA(ig,igp,iqbz))-REAL(fxc_ADA(igp,ig,iqbz)))>tol15.OR.&
1611 !  ABS(AIMAG(fxc_ADA(ig,igp,iqbz))-AIMAG(-fxc_ADA(igp,ig,iqbz)))>tol15) then
1612 !  write(std_out,*) 'Elements:'
1613 !  write(std_out,*) 'fxc_ADA(ig,igp,iqbz):',ig,igp,iqbz,fxc_ADA(ig,igp,iqbz)
1614 !  write(std_out,*) 'fxc_ADA(igp,ig,iqbz):',igp,ig,iqbz,fxc_ADA(igp,ig,iqbz)
1615 !  MSG_ERROR('fxc_ADA not symmetric')
1616 !  end if
1617 !  end do
1618 !  end do
1619 !  end do
1620 
1621 !  write(std_out,*)"kxcr(r=0)",kxcr(1,1)
1622 !  write(std_out,*)"my_kxg(G=0)",my_kxcg(:,1)
1623 !  write(std_out,*)"SUM kxcr/nfft ",SUM(kxcr(:,1))/nfft
1624 !  write(std_out,*)"SUM my_kxg ",SUM(kxcg(:,1))
1625 
1626 !  DEBUG Check FT to real space
1627 !  The real-space expression is:
1628 !  f_xc(R,R')=(1/2)(kappa^2/(4*Pi))
1629 !  \{K_xc[\tilde{n(R)}]+K_xc[\tilde{n(R')}]\}
1630 !  x exp(-kappa|R-R'|)/|R-R'|
1631    ABI_MALLOC(my_fxc_ADA_rrp,(nfft,nfft))
1632    ABI_MALLOC(FT_fxc_ADA_ggpq,(npw,npw,nqibz))
1633    ABI_MALLOC(rvec,(3,nfft))
1634    ABI_MALLOC(dummy,(nfft,nfft))
1635    my_fxc_ADA_rrp=zero; FT_fxc_ADA_ggpq=czero; dummy=czero; rvec=zero
1636 
1637 !  First find coordinates of real-space fft points
1638    igrid = 0
1639    ngfft1 = ngfft(1)
1640    ngfft2 = ngfft(2)
1641    ngfft3 = ngfft(3)
1642    do i3=0,ngfft3-1
1643      difz=dble(i3)/dble(ngfft3)
1644      do i2=0,ngfft2-1
1645        dify=dble(i2)/dble(ngfft2)
1646        do i1=0,ngfft1-1
1647          difx=dble(i1)/dble(ngfft1)
1648          igrid = igrid + 1
1649          rvec(1,igrid)=difx*Cryst%rprimd(1,1)+dify*Cryst%rprimd(1,2)+difz*Cryst%rprimd(1,3)
1650          rvec(2,igrid)=difx*Cryst%rprimd(2,1)+dify*Cryst%rprimd(2,2)+difz*Cryst%rprimd(2,3)
1651          rvec(3,igrid)=difx*Cryst%rprimd(3,1)+dify*Cryst%rprimd(3,2)+difz*Cryst%rprimd(3,3)
1652        end do
1653      end do
1654    end do
1655    if (igrid/=nfft) then
1656      MSG_ERROR('kxc_ADA: igrid not equal to nfft')
1657    end if
1658 
1659 !  Construct kernel in real space
1660    do ir=1,nfft
1661      do irp=ir,nfft
1662        rmrp(:) = rvec(:,ir)-rvec(:,irp)
1663        abs_rmrp = sqrt(dot_product(rmrp,rmrp))
1664        my_fxc_ADA_rrp(ir,irp) = eighth*kappa*kappa*piinv* &
1665        (kxcr(ir,1)+kxcr(irp,1))* &
1666        EXP(-kappa*abs_rmrp)/(abs_rmrp+1.e-3_dp)
1667 !      (a small convergence factor is introduced
1668 !      to avoid a singularity)
1669        my_fxc_ADA_rrp(irp,ir) = my_fxc_ADA_rrp(ir,irp)
1670      end do
1671    end do
1672 
1673 !  Find FFT index for all G
1674    n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
1675 !  Use the following indexing (N means ngfft of the adequate direction)
1676 !  0 1 2 3 ... N/2    -(N-1)/2 ... -1    <= kg
1677 !  1 2 3 4 ....N/2+1  N/2+2    ...  N    <= index
1678    do ig=1,npw
1679      i1=modulo(gvec(1,ig),n1)
1680      i2=modulo(gvec(2,ig),n2)
1681      i3=modulo(gvec(3,ig),n3)
1682      ig_idx_fft(ig)=i1+1+n1*(i2+n2*i3)
1683    end do
1684 !  FT kernel to reciprocal space for each q
1685    do iqbz=1,nqibz
1686      dummy = CMPLX(my_fxc_ADA_rrp,0.0_dp)
1687 !    Multiply with q-point phase factors exp(-iq.r)*f_xc(r,r')*exp(iq.r')
1688      do ir=1,nfft
1689        do irp=1,nfft
1690 !        Calculate q (variables defined for other purposes
1691 !        are being re-used as dummy variables)
1692          q_point(:) = qibz(:,iqbz)
1693          qpg(:) = two_pi*MATMUL(Cryst%gprimd,q_point(:))
1694          abs_qpg_sq = dot_product(qpg(:),rvec(:,ir))
1695          abs_qpgp_sq = dot_product(qpg(:),rvec(:,irp))
1696          dummy(ir,irp) = EXP(-j_dpc*abs_qpg_sq)* &
1697 &         dummy(ir,irp)* &
1698 &         EXP(j_dpc*abs_qpgp_sq)
1699        end do
1700      end do
1701      call fourdp_6d(2,dummy,-1,MPI_enreg_seq,nfft,ngfft,paral_kgb0,0)
1702      do ig=1,npw
1703        do igp=1,npw
1704          FT_fxc_ADA_ggpq(ig,igp,iqbz) = dummy(ig_idx_fft(ig),ig_idx_fft(igp))
1705        end do
1706      end do
1707 
1708 !    Output
1709      msg=''
1710      if (iqbz<10) write(msg,'(a,i1,a)') './debug_fxc_ADA_q',iqbz,'.dat'
1711      if ((iqbz>9).and.(iqbz<100)) write(msg,'(a,i2,a)') './debug_fxc_ADA_q',iqbz,'.dat'
1712      if ((iqbz>99).and.(iqbz<1000)) write(msg,'(a,i3,a)') './debug_fxc_ADA_q',iqbz,'.dat'
1713 
1714      !open(777,file=TRIM(msg),STATUS='REPLACE')
1715      !do igp=1,npw
1716      !  do ig=1,npw
1717      !    write(777,*) ig,igp,REAL(my_fxc_ADA_ggpq(ig,igp,iqbz)),AIMAG(my_fxc_ADA_ggpq(ig,igp,iqbz)), &
1718      !    REAL(FT_fxc_ADA_ggpq(ig,igp,iqbz)),AIMAG(FT_fxc_ADA_ggpq(ig,igp,iqbz)), &
1719      !    ABS(ABS(my_fxc_ADA_ggpq(ig,igp,iqbz))-ABS(FT_fxc_ADA_ggpq(ig,igp,iqbz)))
1720      !  end do
1721      !  write(777,*) ''
1722      !end do
1723      !close(777)
1724 
1725    end do ! iqbz
1726 
1727    MSG_ERROR('Stopping in kxc_ADA for debugging')
1728 
1729    ABI_FREE(rvec)
1730    ABI_FREE(my_fxc_ADA_rrp)
1731    ABI_FREE(FT_fxc_ADA_ggpq)
1732 
1733    if (xcdata%xclevel==2) then
1734      MSG_ERROR(" GGA not implemented for kxc_ADA")
1735    end if !xclevel==2
1736 
1737  end if ! Debugging section
1738 
1739 ! Revert libxc module to the original settings
1740  if (ixc<0) then
1741    call libxc_functionals_end()
1742  end if
1743  if (dtset%ixc<0) then
1744    call libxc_functionals_init(dtset%ixc,dtset%nspden)
1745  end if
1746 
1747  call destroy_mpi_enreg(MPI_enreg_seq)
1748  ABI_FREE(my_kxcg)
1749  ABI_FREE(my_rhor)
1750  ABI_FREE(rhotilder)
1751  ABI_FREE(rhog)
1752  ABI_FREE(kxcr)
1753 
1754 end subroutine kxc_ADA

m_kxc/kxc_alda [ Functions ]

[ Top ] [ m_kxc ] [ Functions ]

NAME

 kxc_alda

FUNCTION

 If option = 1:
  Compute the AL(S)DA kernel in reciprocal space, on the FFT grid.
 If option = 2:
  Only computes the up-down channel of the AL(S)DA kernel, on the
  FFT grid, for use in the BPG kernel.

INPUTS

  dtset <type(dataset_type)>=all input variables in this dataset
  ixc = choice of exchange-correlation functional.
  mpi_enreg=information about MPI parallelization
  nfft = number of fft grid points.
  ngfft(1:3) = integer fft box dimensions, see getng for ngfft(4:8).
  nspden = number of spin-density components.
  option = 1 compute the AL(S)DA kernel in reciprocal space.
         = 2 only computes the up-down channel of the AL(S)DA kernel,
             for use in the BPG kernel.
  rhor(nfft,nspden) = electron density in real space in electrons/bohr**3
   (total in first half and spin-up in second half if nspden = 2).
  rhocut = cut-off density for the local kernels (ALDA, EOK),
           relative to max(rhor(:,:)).
  rprimd(3,3) = dimensional primitive translations for real space in Bohr.

OUTPUT

  kxcg(2,nfft,*) = the AL(S)DA kernel in reciprocal space, on the FFT grid
   (the third dimension is 2*nspden-1 if option = 1, and 1 if option = 2).

SIDE EFFECTS

WARNINGS

 Current restrictions are:
  a - Spin-polarized case not tested.

PARENTS

CHILDREN

      destroy_mpi_enreg,fourdp,fourdp_6d,hartre,initmpi_seq
      libxc_functionals_end,libxc_functionals_init,printxsf,rhotoxc,wrtout
      xcdata_init

SOURCE

417 subroutine kxc_alda(dtset,ixc,kxcg,mpi_enreg,nfft,ngfft,nspden,option,rhor,rhocut,rprimd)
418 
419 
420 !This section has been created automatically by the script Abilint (TD).
421 !Do not modify the following lines by hand.
422 #undef ABI_FUNC
423 #define ABI_FUNC 'kxc_alda'
424 !End of the abilint section
425 
426  implicit none
427 
428 !Arguments -------------------------------------------------------------
429 !scalars
430  integer,intent(in) :: ixc,nfft,nspden,option
431  real(dp),intent(in) :: rhocut
432  type(MPI_type),intent(in) :: mpi_enreg
433  type(dataset_type),intent(in) :: dtset
434 !arrays
435  integer,intent(in) :: ngfft(18)
436  real(dp),intent(in) :: rhor(nfft,2*nspden-1),rprimd(3,3)
437  real(dp),intent(out) :: kxcg(2,nfft,*)
438 
439 !Local variables -------------------------------------------------------
440 !No improved xc quadrature.
441 !No core correction.
442 !Dummy here.
443 !For debugging purposes (see tests below):
444 !integer :: i1,i2,i3,k1,n1,n2,n3
445 !real(dp) :: kx,rho,rhomax,ftest
446 !scalars
447  integer :: ifft,ikxc,isp,n3xccc,ncut,nk3xc,nkxc,optionrhoxc,tim_fourdp
448  real(dp),parameter :: gsqcut=1._dp
449  real(dp) :: enxc,rhocuttot,rhomin,vxcavg
450  character(len=500) :: message
451  type(xcdata_type) :: xcdata
452 !arrays
453  real(dp) :: strsxc(6)
454  real(dp) :: dum(0)
455  real(dp),allocatable :: kxcr(:,:),rhog(:,:),rhorcut(:,:),vhartree(:)
456  real(dp),allocatable :: vxc(:,:),xccc3d(:)
457 
458 !***********************************************************************
459 !For debugging purposes (see tests below):
460 
461 !ftest(i1,n1,k1) = 0._dp+1._dp*cos(k1*two_pi*float(i1)/float(n1))
462 
463 !***********************************************************************
464 
465 !Check input parameters.
466 
467  if (nspden > 2) then
468    message = ' kxc_alda does not work yet for nspden > 2.'
469    MSG_ERROR(message)
470  end if
471 
472 !Allocate memory.
473 
474  ABI_MALLOC(rhorcut,(nfft,nspden))
475  ABI_MALLOC(rhog,(2,nfft))
476  ABI_MALLOC(vhartree,(nfft))
477  ABI_MALLOC(vxc,(nfft,nspden))
478 
479  call xcdata_init(xcdata,dtset=dtset,intxc=0,ixc=ixc,nspden=nspden)
480 
481  ! Reinitialize the libxc module with the overriden values
482  if (dtset%ixc<0) then
483    call libxc_functionals_end()
484  end if
485  if (ixc<0) then
486    call libxc_functionals_init(ixc,dtset%nspden)
487  end if
488 
489 !to be adjusted for the call to rhotoxc
490  nk3xc=1
491 
492 !Cut-off the density.
493 
494  rhorcut(:,:) = rhor(:,:)
495 
496  do isp = 1,nspden
497 
498    rhomin = maxval(rhorcut(:,isp))*rhocut
499 
500    ncut = 0
501    rhocuttot = 0._dp
502 
503    do ifft = 1,nfft
504      if (rhorcut(ifft,isp) < rhomin) then
505        ncut = ncut+1
506        rhocuttot = rhocuttot+rhorcut(ifft,isp)
507        rhorcut(ifft,isp) = rhomin
508      end if
509    end do
510 
511    if (ncut > 0) then
512      write (message,'(a,es10.3,3a,i1,a,i6,a,f6.3,3a,f6.3,a)') &
513 &     'rhocut = ',rhocut,'.',ch10,&
514 &     'For isp = ',isp,' the density was cut-off at ',ncut,' (',100._dp*float(ncut)/float(ifft),'%) grid points.',ch10,&
515 &     'These points account for ',100._dp*rhocuttot/sum(rhor(:,isp)),'% of the total density.'
516      MSG_WARNING(message)
517    end if
518 
519  end do
520 
521 !Calculate the AL(S)DA kernel in real space.
522 
523  rhog(:,:) = 0._dp !We do not need the Hartree potential.
524  tim_fourdp=0
525 
526  if ((option == 1).or.((option == 2).and.(nspden == 2))) then
527 
528    nkxc = 2*nspden-1
529    n3xccc=0
530    ABI_MALLOC(kxcr,(nfft,nkxc))
531    ABI_MALLOC(xccc3d,(n3xccc))
532 
533    optionrhoxc = 2 !See rhotoxc.f
534 
535    call hartre(1,gsqcut,0,mpi_enreg,nfft,ngfft,dtset%paral_kgb,rhog,rprimd,vhartree)
536    call rhotoxc(enxc,kxcr,mpi_enreg,nfft,ngfft,dum,0,dum,0,nkxc,nk3xc,&
537 &   dtset%usepawu==4.or.dtset%usepawu==14,n3xccc,&
538 &   optionrhoxc,dtset%paral_kgb,rhorcut,rprimd,strsxc,1,vxc,vxcavg,xccc3d,xcdata,vhartr=vhartree)
539 
540 !  DEBUG
541 !  fx for tests.
542 !  write (std_out,'(a)') ' kxc_alda: Using exchange-only kernel for tests.'
543 !  rhomin = minval(rhor(:,1))
544 !  rhomax = maxval(rhor(:,1))
545 !  write (std_out,'(a,es12.5,a,es12.5)') ' kxc_alda: rhomin = ',rhomin,' rhomax = ',rhomax
546 !  write (std_out,'(a)') ' kxc_alda: loping below 0.2*rhomax.'
547 !  kx = (3._dp/4._dp)*((3._dp/pi)**(1._dp/3._dp))
548 !  do ifft = 1,nfft
549 !  rho = rhor(ifft,1)
550 !  rho = max(rho,0.2_dp*rhomax)
551 !  kxcr(ifft,1) = -(4._dp/9._dp)*kx*(rho**(-2._dp/3._dp))
552 !  write (std_out,'(i4,a,es12.5)') ifft,': ',kxcr(ifft,1)
553 !  end do
554 !  write (std_out,'(a,es12.5)') 'kxcrmin: ',minval(kxcr(:,1))
555 !  write (std_out,'(a,es12.5)') 'kxcrmax: ',maxval(kxcr(:,1))
556 !  ENDDEBUG
557 
558 !  DEBUG
559 !  test kernel.
560 !  write(std_out,'(a)') ' kxc_alda: Using test kernel for tests.'
561 !  n1 = ngfft(1) ; n2 = ngfft(2) ; n3 = ngfft(3)
562 !  do i3 = 0,n3-1
563 !  do i2 = 0,n2-1
564 !  do i1 = 0,n1-1
565 !  ifft = i1+n1*(i2+n2*i3)+1
566 !  kxcr(ifft,1) = ftest(i1,n1,1)*ftest(i2,n2,2)*ftest(i3,n3,3)
567 !  end do
568 !  end do
569 !  end do
570 !  ENDDEBUG
571 
572 !  Calculate the Fourier transform of the AL(S)DA kernel.
573 
574    if (option == 1) then
575      do ikxc = 1,nkxc
576        call fourdp(1,kxcg(:,:,ikxc),kxcr(:,ikxc),-1,mpi_enreg,nfft,ngfft,dtset%paral_kgb,tim_fourdp)
577      end do
578    else
579      call fourdp(1,kxcg(:,:,1),kxcr(:,2),-1,mpi_enreg,nfft,ngfft,dtset%paral_kgb,tim_fourdp)
580    end if
581 
582  else if ((option == 2).and.(nspden == 1)) then
583 
584    nkxc = 2
585    n3xccc=0
586    ABI_MALLOC(kxcr,(nfft,nkxc))
587    ABI_MALLOC(xccc3d,(n3xccc))
588 
589    optionrhoxc = -2 !See rhotoxc.f
590 
591    call hartre(1,gsqcut,0,mpi_enreg,nfft,ngfft,dtset%paral_kgb,rhog,rprimd,vhartree)
592    call rhotoxc(enxc,kxcr,mpi_enreg,nfft,ngfft,dum,0,dum,0,nkxc,nk3xc,&
593 &   dtset%usepawu==4.or.dtset%usepawu==14,n3xccc,&
594 &   optionrhoxc,dtset%paral_kgb,rhorcut,rprimd,strsxc,1,vxc,vxcavg,xccc3d,xcdata,vhartr=vhartree)
595 
596    kxcr(:,2) = 0.5_dp*kxcr(:,2)
597 
598 !  Calculate the Fourier transform of the up-down channel of the AL(S)DA kernel.
599 
600    call fourdp(1,kxcg(:,:,1),kxcr(:,2),-1,mpi_enreg,nfft,ngfft,dtset%paral_kgb,tim_fourdp)
601 
602  else
603    write (message,'(4a,i0)')'  Invalid option = ',option
604    MSG_ERROR(message)
605  end if
606 
607 !DEBUG
608 !write(std_out,*) ' kxc_alda:  Exc  = ',enxc
609 !write(std_out,*) ' kxc_alda: <Vxc> = ',vxcavg
610 !ENDDEBUG
611 
612 ! Revert libxc module to the original settings
613  if (ixc<0) then
614    call libxc_functionals_end()
615  end if
616  if (dtset%ixc<0) then
617    call libxc_functionals_init(dtset%ixc,dtset%nspden)
618  end if
619 
620 !Free memory.
621  ABI_FREE(rhorcut)
622  ABI_FREE(rhog)
623  ABI_FREE(vhartree)
624  ABI_FREE(vxc)
625  ABI_FREE(kxcr)
626  ABI_FREE(xccc3d)
627 
628 end subroutine kxc_alda

m_kxc/kxc_driver [ Functions ]

[ Top ] [ m_kxc ] [ Functions ]

NAME

 kxc_driver

FUNCTION

 Calculate the exchange-correlation kernel in reciprocal space.
 Require density in real space on the *full* FFT mesh

INPUTS

 Dtset<dataset_type>=all input variables in this dataset
 Cryst<crystal_t>=Info on the crystal structure.
 ixc = choice for the exchange-correlation potential.
 ngfft(18)=contain all needed information about 3D FFT,
  see ~abinit/doc/variables/vargs.htm#ngfft
 nfft_tot = Total number of points on the FFT grid.
 nspden=Number of independent spin densities.
 rhor(nfft_tot,nspden) = the charge density on the full FFT grid.
  (total in first half and spin-up in second half if nspden=2)
 npw: the size of kernel matrix
 dim_kxcg=dimension of the kernel.
 comm=MPI communicator.
 [dbg_mode]=Optional flag used to switch on the debug mode.

OUTPUT

  FIXME: Why are we using nfft_tot instead of the G-sphere
  kxcg(nfft_tot,dim_kxcg) = the exchange-correlation potential on the FFT grid.
  warning: the kernel is not divided by the unit cell volume

NOTES

  No xc quadrature
  No nl core correction

PARENTS

      screening,sigma

CHILDREN

      destroy_mpi_enreg,fourdp,fourdp_6d,hartre,initmpi_seq
      libxc_functionals_end,libxc_functionals_init,printxsf,rhotoxc,wrtout
      xcdata_init

SOURCE

1072 subroutine kxc_driver(Dtset,Cryst,ixc,ngfft,nfft_tot,nspden,rhor,npw,dim_kxcg,kxcg,gvec,comm,dbg_mode)
1073 
1074 
1075 !This section has been created automatically by the script Abilint (TD).
1076 !Do not modify the following lines by hand.
1077 #undef ABI_FUNC
1078 #define ABI_FUNC 'kxc_driver'
1079 !End of the abilint section
1080 
1081  implicit none
1082 
1083 !Arguments ------------------------------------
1084 !scalars
1085  integer,intent(in) :: ixc,npw,nfft_tot,nspden,dim_kxcg,comm
1086  logical,optional,intent(in) :: dbg_mode
1087  type(crystal_t),intent(in) :: Cryst
1088  type(dataset_type),intent(in) :: Dtset
1089 !arrays
1090  integer,intent(in) :: gvec(3,npw),ngfft(18)
1091  real(dp),intent(in) :: rhor(nfft_tot,nspden)
1092  complex(gwpc),intent(out) :: kxcg(nfft_tot,dim_kxcg)
1093 
1094 !Local variables ------------------------------
1095 !scalars
1096  integer,parameter :: paral_kgb0=0
1097  integer :: cplex,i1,i2,i3,ig,igp,iq,ir,n3xccc,ngfft1,ngfft2,izero
1098  integer :: ngfft3,nkxc,option,ikxc,nk3xc,my_rank,master,unt_dmp
1099  real(dp) :: enxc,expo,gpqx,gpqy,gpqz,gsqcut,vxcavg
1100  character(len=500) :: msg,fname
1101  type(xcdata_type) :: xcdata
1102  type(MPI_type) :: MPI_enreg_seq
1103 !arrays
1104  real(dp) :: qphon(3),strsxc(6),dum(0)
1105  real(dp),allocatable :: kxcpw_g(:,:),kxcr(:,:),phas(:,:,:)
1106  real(dp),allocatable :: rhog(:,:),vhartr(:),kxcpw_r(:,:),vxclda(:,:)
1107  real(dp),allocatable :: xccc3d(:),xx(:,:)
1108  real(dp),allocatable :: my_kxcg(:,:)
1109 
1110 !************************************************************************
1111 
1112  ABI_CHECK(Dtset%nsppol==1,'nsppol/=1 not coded')
1113  ABI_CHECK(Dtset%nspden==1,'nsppol/=1 not coded')
1114  ABI_CHECK(nfft_tot==PRODUCT(ngfft(1:3)),"mismatch in nfftot")
1115 
1116 !Fake MPI_type for the sequential part.
1117  call initmpi_seq(MPI_enreg_seq)
1118  call init_distribfft_seq(MPI_enreg_seq%distribfft,'c',ngfft(2),ngfft(3),'all')
1119  my_rank = xmpi_comm_rank(comm)
1120  master=0
1121 
1122  write(msg,'(a,i3)') ' kxc_driver: calculating exchange-correlation kernel using ixc = ',ixc
1123  call wrtout(std_out,msg,'COLL')
1124 
1125  call xcdata_init(xcdata,dtset=Dtset,intxc=0,ixc=ixc,nspden=nspden)
1126 
1127  if (ALL(xcdata%xclevel/=(/1,2/))) then
1128    write(msg,'(a,i0)')"Unsupported xclevel = ",xcdata%xclevel
1129    MSG_ERROR(msg)
1130  end if
1131 
1132  ngfft1=ngfft(1)
1133  ngfft2=ngfft(2)
1134  ngfft3=ngfft(3)
1135 
1136  if (ixc>=1.and.ixc<11) then ! LDA case
1137    nkxc= 2*min(nspden,2)-1   ! 1 or 3
1138  else                        ! GGA case
1139    nkxc=12*min(nspden,2)-5   ! 7 or 19
1140    ABI_CHECK(dtset%xclevel==2,"Functional should be GGA")
1141    MSG_ERROR("GGA functional not tested")
1142  end if
1143 
1144  ABI_MALLOC(kxcr,(nfft_tot,nkxc))
1145 
1146 !gsqcut and rhog are zeroed because they are not used by rhotoxc if 1<=ixc<=16 and option=0
1147  gsqcut=zero
1148 
1149  ABI_MALLOC(rhog,(2,nfft_tot))
1150  ABI_MALLOC(vhartr,(nfft_tot))
1151  rhog(:,:)=zero
1152 !MG FIXME this is the 3D core electron density for XC core correction (bohr^-3)
1153 !should implement the non linear core correction
1154  n3xccc=0
1155  ABI_MALLOC(xccc3d,(n3xccc))
1156  ABI_MALLOC(vxclda,(nfft_tot,nspden))
1157 
1158  option=2 ! 2 for Hxc and kxcr (no paramagnetic part if nspden=1)
1159  qphon =zero
1160 
1161 !to be adjusted for the call to rhotoxc
1162  nk3xc=1
1163  izero=0
1164 
1165  ! Reinitialize the libxc module with the overriden values
1166  if (dtset%ixc<0) then
1167    call libxc_functionals_end()
1168  end if
1169  if (ixc<0) then
1170    call libxc_functionals_init(ixc,Dtset%nspden)
1171  end if
1172 
1173  call hartre(1,gsqcut,izero,MPI_enreg_seq,nfft_tot,ngfft,dtset%paral_kgb,rhog,Cryst%rprimd,vhartr)
1174 
1175 !Compute the kernel.
1176  call rhotoxc(enxc,kxcr,MPI_enreg_seq,nfft_tot,ngfft,&
1177 & dum,0,dum,0,nkxc,nk3xc,dtset%usepawu==4.or.dtset%usepawu==14,&
1178 & n3xccc,option,Dtset%paral_kgb,rhor,Cryst%rprimd,&
1179 & strsxc,1,vxclda,vxcavg,xccc3d,xcdata,vhartr=vhartr)
1180 
1181  ABI_FREE(rhog)
1182  ABI_FREE(vhartr)
1183 !DEBUG print Kxc
1184  if (present(dbg_mode)) then
1185    if (dbg_mode .and. my_rank==master) then
1186      fname = 'xc_Kxc.xsf'
1187      if (open_file(fname,msg,newunit=unt_dmp,status='unknown',form='formatted') /= 0) then
1188        MSG_ERROR(msg)
1189      end if
1190      call printxsf(ngfft1,ngfft2,ngfft3,kxcr(:,1),Cryst%rprimd,(/zero,zero,zero/),&
1191 &     Cryst%natom,Cryst%ntypat,Cryst%typat,Cryst%xcart,Cryst%znucl,unt_dmp,0)
1192      close(unt_dmp)
1193    end if
1194  end if
1195 !DEBUG
1196 
1197  ABI_FREE(xccc3d)
1198  ABI_FREE(vxclda)
1199 
1200  ABI_MALLOC(my_kxcg,(2,nfft_tot))
1201 
1202  do ikxc=1,nkxc
1203    call fourdp(1,my_kxcg,kxcr(:,ikxc),-1,MPI_enreg_seq,nfft_tot,ngfft,paral_kgb0,0)
1204    kxcg(:,ikxc)=CMPLX(my_kxcg(1,:),my_kxcg(2,:))
1205  end do
1206 
1207 !write(std_out,*)"kxcr(r=0)",kxcr(1,1)
1208 !write(std_out,*)"my_kxg(G=0)",my_kxcg(:,1)
1209 !write(std_out,*)"SUM kxcr/nfft_tot ",SUM(kxcr(:,1))/nfft_tot
1210 !write(std_out,*)"SUM my_kxg ",SUM(kxcg(:,1))
1211 
1212  ABI_FREE(my_kxcg)
1213 
1214 !MG this part is never executed, but one should use dfpt_mkvxc for the GGA kernel.
1215  if (xcdata%xclevel==2) then
1216    MSG_ERROR("check GGA implementation")
1217    cplex=2
1218    ABI_MALLOC(phas,(cplex*nfft_tot,npw,nspden))
1219    ABI_MALLOC(kxcpw_r,(cplex*nfft_tot,nspden))
1220    ABI_MALLOC(xx,(3,nfft_tot))
1221    ABI_MALLOC(kxcpw_g,(2,nfft_tot))
1222 
1223    kxcg = czero
1224 
1225 !  find the coordinates for all r in the FFT grid
1226    ir=0
1227    do i3=1,ngfft3
1228      do i2=1,ngfft2
1229        do i1=1,ngfft1
1230          ir=ir+1
1231          xx(1,ir)=dble((i1-1))/ngfft1
1232          xx(2,ir)=dble((i2-1))/ngfft2
1233          xx(3,ir)=dble((i3-1))/ngfft3
1234        end do
1235      end do
1236    end do
1237 
1238    do iq=1,1
1239 
1240 !    Calculate at once exp(i(G+q).r), for all possible q,G,r
1241      do ig=1,npw
1242        gpqx=dble(gvec(1,ig))
1243        gpqy=dble(gvec(2,ig))
1244        gpqz=dble(gvec(3,ig))
1245        do ir=1,nfft_tot
1246          expo=gpqx*xx(1,ir)+gpqy*xx(2,ir)+gpqz*xx(3,ir)
1247          phas(2*ir-1,ig,1)= cos(two_pi*expo)
1248          phas(2*ir,ig,1) =  sin(two_pi*expo)
1249        end do
1250      end do
1251 
1252 !    Calculate $K(G,G'',q)=\frac{1}{nfft_tot}\sum_{r} exp(-i(q+G_{2}).r_{2} kxcr(r_{1}r_{2}) exp(i(q+G_{1}).r_{1} $
1253 
1254      do igp=1,npw
1255 
1256        kxcpw_r(:,:)=zero
1257 
1258        call dfpt_mkvxc(cplex,ixc,kxcr,MPI_enreg_seq,nfft_tot,ngfft,dum,0,dum,0,nkxc,nspden,&
1259 &       n3xccc,option,paral_kgb0,qphon(:),phas(:,igp,:),Cryst%rprimd,1,kxcpw_r,xccc3d)
1260 
1261 !      FFT the first index to --> to G space
1262        call fourdp(cplex,kxcpw_g(:,:),kxcpw_r(:,1),-1,MPI_enreg_seq,nfft_tot,ngfft,paral_kgb0,0)
1263 
1264 !      kxcg(:,igp,iq)=CMPLX(kxcpw_g(1,igfft(:)),kxcpw_g(2,igfft(:)))
1265 !      kxcg(:,igp)=CMPLX(kxcpw_g(1,igfft(:)),kxcpw_g(2,igfft(:)))
1266 
1267      end do ! igp
1268    end do ! iq
1269 
1270    ABI_FREE(phas)
1271    ABI_FREE(kxcpw_r)
1272    ABI_FREE(xx)
1273    ABI_FREE(kxcpw_g)
1274  end if !xclevel==2
1275 
1276 ! Revert libxc module to the original settings
1277  if (ixc<0) then
1278    call libxc_functionals_end()
1279  end if
1280  if (dtset%ixc<0) then
1281    call libxc_functionals_init(dtset%ixc,dtset%nspden)
1282  end if
1283 
1284  call destroy_mpi_enreg(MPI_enreg_seq)
1285  ABI_FREE(kxcr)
1286 
1287 end subroutine kxc_driver

m_kxc/kxc_eok [ Functions ]

[ Top ] [ m_kxc ] [ Functions ]

NAME

 kxc_eok

FUNCTION

  Compute the linear (ixceok = 1) or non-linear (ixceok = 2)
  energy optimized kernel of Dobson and Wang, in reciprocal
  space, on the FFT grid.
  See J. Dobson and J. Wang, Phys. Rev. B 62, 10038 (2000) [[cite:Dobson2000]].

INPUTS

  ixceok = 1 linear energy optimized kernel.
         = 2 non-linear energy optimized kernel.
  mpi_enreg=information about MPI parallelization
  nfft = number of fft grid points.
  ngfft(1:3) = integer fft box dimensions, see getng for ngfft(4:8).
  nspden = number of spin-density components.
  rhor(nfft,nspden) = electron density in real space in electrons/bohr**3
   (total in first half and spin-up in second half if nspden = 2).
  rhocut = cut-off density for the local kernels (ALDA, EOK),
           relative to max(rhor(:,:)).

OUTPUT

  kxcg(2,nfft,2*nspden-1) = the EOK kernel in reciprocal space, on the FFT grid.

SIDE EFFECTS

WARNINGS

PARENTS

CHILDREN

      destroy_mpi_enreg,fourdp,fourdp_6d,hartre,initmpi_seq
      libxc_functionals_end,libxc_functionals_init,printxsf,rhotoxc,wrtout
      xcdata_init

SOURCE

 906 subroutine kxc_eok(ixceok,kxcg,mpi_enreg,nfft,ngfft,nspden,paral_kgb,rhor,rhocut)
 907 
 908 
 909 !This section has been created automatically by the script Abilint (TD).
 910 !Do not modify the following lines by hand.
 911 #undef ABI_FUNC
 912 #define ABI_FUNC 'kxc_eok'
 913 !End of the abilint section
 914 
 915  implicit none
 916 
 917 !Arguments -------------------------------------------------------------
 918 !scalars
 919  integer,intent(in) :: ixceok,nfft,nspden,paral_kgb
 920  real(dp),intent(in) :: rhocut
 921  type(MPI_type),intent(in) :: mpi_enreg
 922 !arrays
 923  integer,intent(in) :: ngfft(18)
 924  real(dp),intent(in) :: rhor(nfft,2*nspden-1)
 925  real(dp),intent(out) :: kxcg(2,nfft,2*nspden-1)
 926 
 927 !Local variables -------------------------------------------------------
 928 !Maximum value allowed for rs.
 929 !scalars
 930  integer :: ifft,ikxc,ncut,nkxc,nlop,tim_fourdp
 931  real(dp),parameter :: rslim=50._dp
 932  real(dp) :: a2,a3,a4,rho,rhocuttot,rhomin,rs
 933  character(len=500) :: message
 934 !arrays
 935  real(dp),allocatable :: kxcr(:,:)
 936 
 937 !***********************************************************************
 938 
 939 !Check input parameters.
 940 
 941  if (nspden > 1) then
 942    message = ' kxc_eok does not work yet for nspden > 1.'
 943    MSG_ERROR(message)
 944  end if
 945 
 946 !Values of a2, a3 and a4 for case 1
 947  a2 = 0.0_dp
 948  a3 = 0.0_dp
 949  a4 = 0.0_dp
 950 
 951  select case (ixceok)
 952    case (1)
 953      a2 = -0.51887_dp
 954      a3 =  4.9359d-03
 955      a4 = -5.9603d-05
 956    case (2)
 957      a2 = -0.50044_dp
 958      a3 =  4.9653d-03
 959      a4 = -3.3660d-05
 960    case default
 961      message =  ' kxc_eok: ixceok /= 1 (linear EOK) or 2 (non-linear EOK).'
 962      MSG_ERROR(message)
 963  end select
 964 
 965 !Allocate memory.
 966 
 967  nkxc = 2*nspden-1
 968 
 969  ABI_MALLOC(kxcr,(nfft,nkxc))
 970 
 971 !Calculate the energy optimized kernel in real space.
 972 
 973  nlop = 0
 974 
 975  rhomin = rhocut*maxval(rhor(:,:))
 976 
 977  ncut = 0
 978  rhocuttot = 0._dp
 979 
 980  do ifft = 1,nfft
 981 
 982    rho = rhor(ifft,1)
 983 
 984    if (rho < rhomin) then
 985      ncut = ncut+1
 986      rhocuttot = rhocuttot+rho
 987      rho = rhomin
 988    end if
 989 
 990    rs = (3._dp/(4._dp*pi*rho))**(1._dp/3._dp)
 991 
 992    if (rs > rslim) then
 993      rs = rslim
 994      nlop = nlop+1
 995    end if
 996 
 997    kxcr(ifft,1) = a2*rs**2+a3*rs**3+a4*rs**4
 998 
 999  end do
1000 
1001  if (ncut > 0) then
1002    write (message,'(a,es10.3,3a,i1,a,i6,a,f6.3,3a,f6.3,a)') &
1003 &   'rhocut = ',rhocut,'.',ch10,&
1004 &   'For isp = ',1,' the density was cut-off at ',ncut,' (',100._dp*float(ncut)/float(ifft),'%) grid points.',ch10,&
1005 &   'These points account for ',100._dp*rhocuttot/sum(rhor(:,1)),'% of the total density.'
1006    MSG_WARNING(message)
1007  end if
1008 
1009  if (nlop > 0) then
1010    write (message,'(a,f6.2,a,i6,a,f6.3,a)') &
1011 &   'rs still exceeds ',rslim,' Bohr at ',nlop,' (',100._dp*float(nlop)/float(ifft),'%) grid points (after cut-off).'
1012    MSG_WARNING(message)
1013  end if
1014 
1015 !Calculate the Fourier transform of the energy optimized kernel.
1016  tim_fourdp=0
1017  do ikxc = 1,nkxc
1018    call fourdp(1,kxcg(:,:,ikxc),kxcr(:,ikxc),-1,mpi_enreg,nfft,ngfft,paral_kgb,tim_fourdp)
1019  end do
1020 
1021 !Free memory.
1022 
1023  ABI_FREE(kxcr)
1024 
1025 end subroutine kxc_eok

m_kxc/kxc_local [ Functions ]

[ Top ] [ m_kxc ] [ Functions ]

NAME

 kxc_local

FUNCTION

 In a planewave basis set, the matrix of a local xc kernel:
  $f_{\rm xc}(\vec{r},\vec{r}') = f(\vec{r})\delta(\vec{r}-\vec{r}')$
 is just:
  $f_{\rm xc}(\vec{G},\vec{G}') = f(\vec{G}-\vec{G}')$.
 This subroutine calculates the matrix of such a local xc kernel
 given $f(\vec{G})$ on the FFT grid.

INPUTS

  ispxc = 1 for the up-up spin channel.
        = 2 for the up-down (and down-up) spin channels.
        = 3 for the down-down spin channel.
  ispxc must be 1 if nspden = 1.
  kg_diel(3,npwdiel) = reduced planewave coordinates for the kxc matrix.
  kxcg(2,nfft) = $f(\vec{G})$ on the FFT grid.
  nfft = number of fft grid points.
  ngfft(1:3) = integer fft box dimensions, see getng for ngfft(4:8).
  npwdiel = number of planewaves for the susceptibility matrix.
  nspden = number of spin-density components.
  option = 0 do not compute the first row and column of the matrix of the
             xc kernel (which we assume to the G = 0 row and column).
        /= 0 compute the full matrix of the xc kernel.

OUTPUT

  kxc(2,npwdiel,nspden,npwdiel,nspden) = the matrix of the xc kernel.

PARENTS

CHILDREN

      destroy_mpi_enreg,fourdp,fourdp_6d,hartre,initmpi_seq
      libxc_functionals_end,libxc_functionals_init,printxsf,rhotoxc,wrtout
      xcdata_init

SOURCE

213 subroutine kxc_local(ispxc,kg_diel,kxc,kxcg,nfft,ngfft,npwdiel,nspden,option)
214 
215 
216 !This section has been created automatically by the script Abilint (TD).
217 !Do not modify the following lines by hand.
218 #undef ABI_FUNC
219 #define ABI_FUNC 'kxc_local'
220 !End of the abilint section
221 
222  implicit none
223 
224 !Arguments -------------------------------------------------------------
225 !scalars
226  integer,intent(in) :: ispxc,nfft,npwdiel,nspden,option
227 !arrays
228  integer,intent(in) :: kg_diel(3,npwdiel),ngfft(18)
229  real(dp),intent(in) :: kxcg(2,nfft)
230  real(dp),intent(out) :: kxc(2,npwdiel,nspden,npwdiel,nspden)
231 
232 !Local variables -------------------------------------------------------
233 !For debbuging purposes:
234 !real(dp) :: c1,c2,c3
235 !scalars
236  integer :: i1,i2,i3,ifft,ipw1,ipw2,ipwstart,isp1,isp2,j1,j2,j3,k1,k2,k3,n1,n2
237  integer :: n3
238  logical :: ok
239  character(len=500) :: message
240 
241 !***********************************************************************
242 
243 !Check input parameters.
244 
245  if (nspden > 2) then
246    message =' kxc_local does not work yet for nspden > 2.'
247    MSG_ERROR(message)
248  end if
249 
250  isp1 = 1
251  isp2 = 1
252  ok = .true.
253  if (nspden == 1) then
254    select case (ispxc)
255      case (1)
256        isp1 = 1
257        isp2 = 1
258      case default
259        ok = .false.
260    end select
261  else
262    select case (ispxc)
263      case (1)
264        isp1 = 1
265        isp2 = 1
266      case (2)
267        isp1 = 1
268        isp2 = 2
269      case (3)
270        isp1 = 2
271        isp2 = 2
272      case default
273        ok = .false.
274    end select
275  end if
276 
277  if (.not.ok) then
278    write (message,'(2(a,i0))')'  The input ispxc = ',ispxc,' is not compatible with nspden = ',nspden
279    MSG_BUG(message)
280  end if
281 
282  if (option == 0) then
283    ipwstart = 2
284    kxc(:,1,isp1,:,isp2) = 0._dp
285    kxc(:,:,isp1,1,isp2) = 0._dp
286  else
287    ipwstart = 1
288  end if
289 
290 !Calculate the xc matrix.
291 
292  n1 = ngfft(1) ; n2 = ngfft(2) ; n3 = ngfft(3)
293 
294  do ipw2 = ipwstart,npwdiel
295 
296    j1 = kg_diel(1,ipw2) ; j2 = kg_diel(2,ipw2) ; j3 = kg_diel(3,ipw2)
297 
298 !  Fill the diagonal.
299 
300    kxc(:,ipw2,isp1,ipw2,isp2) = kxcg(:,1)
301 
302 !  Fill the off-diagonal elements.
303 
304    do ipw1 = ipw2+1,npwdiel
305 
306      i1 = kg_diel(1,ipw1) ; i2 = kg_diel(2,ipw1) ; i3 = kg_diel(3,ipw1)
307 
308 !    Compute the difference between G vectors.
309 !    The use of two mod calls handles both i1-j1 >= n1 AND i1-j1 < 0.
310 
311      k1 = mod(n1+mod(i1-j1,n1),n1)
312      k2 = mod(n2+mod(i2-j2,n2),n2)
313      k3 = mod(n3+mod(i3-j3,n3),n3)
314 
315      ifft = k1+n1*(k2+n2*k3)+1
316 
317      kxc(1,ipw1,isp1,ipw2,isp2) =  kxcg(1,ifft)
318      kxc(2,ipw1,isp1,ipw2,isp2) =  kxcg(2,ifft)
319 
320      kxc(1,ipw2,isp1,ipw1,isp2) =  kxcg(1,ifft)
321      kxc(2,ipw2,isp1,ipw1,isp2) = -kxcg(2,ifft)
322 
323    end do
324 
325  end do
326 
327 !If needed, copy the up-down to the down-up spin channel.
328 
329  if (ispxc == 2) then
330    do ipw2 = 1,npwdiel
331      do ipw1 = 1,npwdiel
332        kxc(1,ipw2,isp2,ipw1,isp1) =  kxc(1,ipw1,isp1,ipw2,isp2)
333        kxc(2,ipw2,isp2,ipw1,isp1) = -kxc(2,ipw1,isp1,ipw2,isp2)
334      end do
335    end do
336  end if
337 
338 !DEBUG
339 !See kxc_alda.f, "test kernel" DEBUG section.
340 !do ipw2 = 1,npwdiel
341 !j1 = kg_diel(1,ipw2) ; j2 = kg_diel(2,ipw2) ; j3 = kg_diel(3,ipw2)
342 !do ipw1 = ipw2+1,npwdiel
343 !i1 = kg_diel(1,ipw1) ; i2 = kg_diel(2,ipw1) ; i3 = kg_diel(3,ipw1)
344 !k1 = mod(n1+mod(i1-j1,n1),n1)
345 !k2 = mod(n2+mod(i2-j2,n2),n2)
346 !k3 = mod(n3+mod(i3-j3,n3),n3)
347 !ifft = k1+n1*(k2+n2*k3)+1
348 !c1 = 0._dp ; c2 = 0._dp ; c3 = 0._dp
349 !if (i1-j1 ==  0) c1 = c1+0.0_dp
350 !if (i2-j2 ==  0) c2 = c2+0.0_dp
351 !if (i3-j3 ==  0) c3 = c3+0.0_dp
352 !if (i1-j1 ==  1) c1 = c1+0.5_dp
353 !if (i2-j2 ==  2) c2 = c2+0.5_dp
354 !if (i3-j3 ==  3) c3 = c3+0.5_dp
355 !if (i1-j1 == -1) c1 = c1+0.5_dp
356 !if (i2-j2 == -2) c2 = c2+0.5_dp
357 !if (i3-j3 == -3) c3 = c3+0.5_dp
358 !if ((abs(kxcg(1,ifft)-c1*c2*c3) > tol10).or.(abs(kxcg(2,ifft)) > tol10)) then
359 !write (std_out,*) ' i1 i2 i3 ifft: ',i1,i2,i3,ifft
360 !write (std_out,*) ' exp.: ',c1*c2*c3,' got: ',kxcg(:,ifft)
361 !end if
362 !end do
363 !end do
364 !ENDDEBUG
365 
366 end subroutine kxc_local

m_kxc/kxc_pgg [ Functions ]

[ Top ] [ m_kxc ] [ Functions ]

NAME

 kxc_pgg

FUNCTION

 Compute the PGG-exchange kernel in reciprocal space
 (Phys. Rev. Lett. 76, 1212 (1996) [[cite:Petersilka1996]]).

INPUTS

  gmet=reciprocal space metrix (bohr**-2)
  npw=number of plane waves
  rcut_coulomb=real space cutoff radius for Coulomb interaction (bohr)
  susmat(2,npw,npw)=density weighted squared density matrix in reciprocal space
  ucvol=unit cell volume (bohr**3)

OUTPUT

  khxcg(2,npwdiel,nspden,npwdiel,nspden)=PGG-exhange kernel in G space, at
       full interaction strength

NOTES

 The density weighted squared density matrix (actually the reduced=summed-over-spin
 density matrix) is convolved with the spherically cutoff Coulomb interaction.

WARNINGS

 a - 'rcut_coulomb' should be chosen consistently with cutoffs elsewhere,
     for instance dieltcel8.f
 b - applicable for spin-polarized case as well, through input 'susmat',
     but this has not been checked

TODO

 If simply the squared density matrix is input through 'susmat' the
 exchange energy is recovered as the zero-G component of the resulting 'khxcg'
 (then not the kernel of course). This could help to check convergence
 with respect to 'npw'. See +ex_pgg comment.

PARENTS

CHILDREN

      destroy_mpi_enreg,fourdp,fourdp_6d,hartre,initmpi_seq
      libxc_functionals_end,libxc_functionals_init,printxsf,rhotoxc,wrtout
      xcdata_init

SOURCE

677 subroutine kxc_pgg(gmet,kg,khxcg,npw,rcut_coulomb,susmat,ucvol)
678 
679 
680 !This section has been created automatically by the script Abilint (TD).
681 !Do not modify the following lines by hand.
682 #undef ABI_FUNC
683 #define ABI_FUNC 'kxc_pgg'
684 !End of the abilint section
685 
686  implicit none
687 
688 !Arguments ------------------------------------
689 !scalars
690  integer,intent(in) :: npw
691  real(dp),intent(in) :: rcut_coulomb,ucvol
692 !arrays
693  integer,intent(in) :: kg(3,npw)
694  real(dp),intent(in) :: gmet(3,3),susmat(2,npw,npw)
695  real(dp),intent(out) :: khxcg(2,npw,npw)
696 
697 !Local variables-------------------------------
698 !scalars
699  integer :: i1,i2,i3,ig,ii,ikg11,ikg12,ikg13,ikg21,ikg22,ikg23,ipw1,ipw2
700  integer :: j1,j2,j3,jg,jj
701  real(dp),parameter :: diffgsq=1.d-2
702  real(dp) :: gred1,gred2,gred3,gsquar,tpisq
703 !arrays
704  integer :: kgmax(3)
705  integer,allocatable :: index_g_inv(:,:,:),jgarr(:)
706  real(dp),allocatable :: gsq(:),sumg(:),vcoul(:)
707 
708 ! *************************************************************************
709 
710 !DEBUG
711 !write(std_out,*) '%kxc_pgg: enter'
712 !write(std_out,*) 'npw=',npw
713 !ENDDEBUG
714 
715 !tpisq is (2 Pi) **2:
716  tpisq=(two_pi)**2
717 
718  kgmax(:)=0
719  do ipw1=1,npw
720    do jj=1,3
721      kgmax(jj)=max( kg(jj,ipw1), kgmax(jj) )
722    end do
723  end do
724 
725 !DEBUG
726 !write(std_out,*) 'kgmax:',kgmax(1:3)
727 !ENDDEBUG
728 
729 !Perform allocations
730  ABI_MALLOC(index_g_inv,(-2*kgmax(1):2*kgmax(1),-2*kgmax(2):2*kgmax(2),-2*kgmax(3):2*kgmax(3)))
731  ABI_MALLOC(jgarr,(npw))
732  ABI_MALLOC(gsq,(npw))
733  ABI_MALLOC(sumg,(2))
734  ABI_MALLOC(vcoul,(npw))
735 
736 !DEBUG
737 !write(std_out,*) '%kxc_pg: creating plane wave index and coulomb potential'
738 !ENDDEBUG
739  index_g_inv(:,:,:)=0
740  do ipw1=1,npw
741    index_g_inv(kg(1,ipw1),kg(2,ipw1),kg(3,ipw1))=ipw1
742 
743 !  DEBUG
744 !  write(std_out,'(i5,2x,3i3,2x,i4)') ipw1,kg(1,ipw1),kg(2,ipw1),kg(3,ipw1)
745 !  ENDDEBUG
746 
747    gred1=dble(kg(1,ipw1))
748    gred2=dble(kg(2,ipw1))
749    gred3=dble(kg(3,ipw1))
750    gsquar=tpisq*(gmet(1,1)*gred1**2+gmet(2,2)*gred2**2+gmet(3,3)*gred3**2 &
751 &   +2.0_dp*( (gmet(1,2)*gred2+gmet(1,3)*gred3)* gred1 +      &
752 &   gmet(2,3)*gred2*gred3) )
753 !  Distinguish G=0 from other elements
754    if(gsquar > 1.0d-12)then
755      vcoul(ipw1)=four_pi/gsquar*(1._dp-cos(sqrt(gsquar)*rcut_coulomb))
756    else
757      vcoul(ipw1)=four_pi*0.5_dp*rcut_coulomb**2
758    end if
759 
760  end do
761 
762 !DEBUG
763 !write(std_out,*) '%kxc_pg: starting convolution integral'
764 !ENDDEBUG
765 
766 !loop over G1,G2 components of the density matrix
767  do ipw2=1,npw
768    ikg21=kg(1,ipw2)
769    ikg22=kg(2,ipw2)
770    ikg23=kg(3,ipw2)
771 
772    do ii=1,npw
773      j1=ikg21-kg(1,ii)
774      j2=ikg22-kg(2,ii)
775      j3=ikg23-kg(3,ii)
776      jgarr(ii)=index_g_inv(j1,j2,j3)
777    end do
778 
779    do ipw1=1,ipw2
780      ikg11=kg(1,ipw1)
781      ikg12=kg(2,ipw1)
782      ikg13=kg(3,ipw1)
783 
784 !    do the convolution integral
785      sumg(:)=0._dp
786      do ii=1,npw
787 
788        if( jgarr(ii) /= 0 ) then
789 
790          i1=ikg11-kg(1,ii)
791          i2=ikg12-kg(2,ii)
792          i3=ikg13-kg(3,ii)
793 
794 !        j1=ikg21-kg(1,ii)
795 !        j2=ikg22-kg(2,ii)
796 !        j3=ikg23-kg(3,ii)
797 
798          ig=index_g_inv(i1,i2,i3)
799 !        jg=index_g_inv(j1,j2,j3)
800 
801          if( ig /= 0 ) then
802 
803            jg=jgarr(ii)
804 
805 !          DEBUG
806 !          write(std_out,'(i5,2x,3i3,1x,3i3,2x,2i4)') ii,i1,i2,i3,&
807 !          &                                             kg(1,jg),kg(2,jg),kg(3,jg),&
808 !          &                                             ig,jg
809 !          ENDDEBUG
810 
811            sumg(1)=sumg(1)+susmat(1,ig,jg)*vcoul(ii)
812            sumg(2)=sumg(2)+susmat(2,ig,jg)*vcoul(ii)
813 
814          end if
815 
816        end if
817 
818      end do
819      khxcg(:,ipw1,ipw2)=-sumg(:)*ucvol
820 
821 !    DEBUG
822 !    if(ipw1==ipw2) write(std_out,'(2i4,2(1x,es14.6))') ipw1,ipw2,khxcg(1,ipw1,ipw1),vcoul(ipw1)
823 !    write(std_out,'(2i4,3(1x,es14.6))') ipw1,ipw2,khxcg(1:2,ipw1,ipw2),vcoul(ipw1)
824 !    ENDDEBUG
825 
826    end do
827  end do
828 
829 !DEBUG
830 !verify hermiticity, note: ipw1 loop above must end at npw
831 !write(std_out,*) '%kxc_pgg: check hermiticity of pgg kernel'
832 !do ipw2=1,npw,max(2,npw/10)
833 !do ipw1=ipw2,npw,max(2,npw/10)
834 !write(std_out,'(2i4,2(1x,es14.6))') ipw1,ipw2,&
835 !&   khxcg(1,ipw1,ipw2)-khxcg(1,ipw2,ipw1),&
836 !&   khxcg(2,ipw1,ipw2)+khxcg(2,ipw2,ipw1)
837 !end do
838 !end do
839 !ENDDEBUG
840 
841 !Impose hermiticity
842  write(std_out,*) '%kxc_pg: imposing hermiticity'
843  do ipw2=1,npw
844    do ipw1=ipw2+1,npw
845      khxcg(1,ipw1,ipw2)= khxcg(1,ipw2,ipw1)
846      khxcg(2,ipw1,ipw2)=-khxcg(2,ipw2,ipw1)
847    end do
848  end do
849 
850 !DEBUG
851 !write(std_out,'(a10,2(1x,es20.12))') '+ex_pgg? ', 0.5_dp*khxcg(1,1,1)/ucvol
852 !ENDDEBUG
853 
854  ABI_FREE(index_g_inv)
855  ABI_FREE(jgarr)
856  ABI_FREE(gsq)
857  ABI_FREE(sumg)
858  ABI_FREE(vcoul)
859 
860 !DEBUG
861 !write(std_out,*) '%kxc_pgg: done'
862 !ENDDEBUG
863 
864 end subroutine kxc_pgg

m_kxc/kxc_rpa [ Functions ]

[ Top ] [ m_kxc ] [ Functions ]

NAME

 kxc_rpa

FUNCTION

 Return the Hartree kernel:
  If option = 0, the bare Hartree kernel:
   krpa(ipw) = 4.0*pi/gsq(ipw) if gsq(ipw) /= 0.,
   krpa(ipw) = 0.0             if gsq(ipw) == 0. (1 <= ipw <= npw).
  If option /= 0, the Hartree kernel with a cut-off in real space beyond rcut_coulomb:
   krpa(ipw) = (4.0*pi/gsq(ipw))*(1.0-cos(sqrt(gsq(ipw))*rcut_coulomb)) if gsq(ipw) /= 0.,
   krpa(ipw) =  2.0*pi*rcut_coulomb**2                                  if gsq(ipw) == 0.

INPUTS

  gsq(npw) = the squared norm of the planewaves.
  npw = number of planewaves in the gsq array.
  option = 0 for the bare Hartree kernel, /=0 for the cut-off Hartree kernel.
  rcut_coulomb = real space cut-off radius for the Coulomb interaction in Bohr.

OUTPUT

  krpa(npw) = the Hartree kernel.

PARENTS

CHILDREN

      destroy_mpi_enreg,fourdp,fourdp_6d,hartre,initmpi_seq
      libxc_functionals_end,libxc_functionals_init,printxsf,rhotoxc,wrtout
      xcdata_init

SOURCE

114 subroutine kxc_rpa(gsq,krpa,npw,option,rcut_coulomb)
115 
116 
117 !This section has been created automatically by the script Abilint (TD).
118 !Do not modify the following lines by hand.
119 #undef ABI_FUNC
120 #define ABI_FUNC 'kxc_rpa'
121 !End of the abilint section
122 
123  implicit none
124 
125 !Arguments -------------------------------------------------------------
126 !scalars
127  integer,intent(in) :: npw,option
128  real(dp),intent(in) :: rcut_coulomb
129 !arrays
130  real(dp),intent(in) :: gsq(npw)
131  real(dp),intent(out) :: krpa(npw)
132 
133 !Local variables -------------------------------------------------------
134 !scalars
135  integer :: ipw
136 
137 !***********************************************************************
138 
139  if (option == 0) then
140 
141 !  Compute the bare Hartree kernel.
142 
143    do ipw = 1,npw
144 
145      if (gsq(ipw) > tol12) then
146        krpa(ipw) = four_pi/gsq(ipw)
147      else
148        krpa(ipw) = 0._dp
149      end if
150 
151    end do
152 
153  else
154 
155 !  Compute the Hartree kernel with a cut-off in real space beyond rcut_coulomb:
156 
157    do ipw = 1,npw
158 
159      if (gsq(ipw) > tol12) then
160        krpa(ipw) = (four_pi/gsq(ipw))*(1._dp-cos(sqrt(gsq(ipw))*rcut_coulomb))
161      else
162        krpa(ipw) = two_pi*rcut_coulomb**2
163      end if
164 
165    end do
166 
167  end if
168 
169 end subroutine kxc_rpa