TABLE OF CONTENTS
- ABINIT/m_kxc
- m_kxc/kxc_ADA
- m_kxc/kxc_alda
- m_kxc/kxc_driver
- m_kxc/kxc_eok
- m_kxc/kxc_local
- m_kxc/kxc_pgg
- m_kxc/kxc_rpa
ABINIT/m_kxc [ 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-2024 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)
SOURCE
34 #if defined HAVE_CONFIG_H 35 #include "config.h" 36 #endif 37 38 #include "abi_common.h" 39 40 MODULE m_kxc 41 42 use defs_basis 43 use m_abicore 44 use m_errors 45 use m_xmpi 46 use m_crystal 47 use m_distribfft 48 use m_xcdata 49 use libxc_functionals 50 use m_dtset 51 52 use defs_abitypes, only : MPI_type 53 use m_io_tools, only : open_file 54 use m_pptools, only : printxsf 55 use m_numeric_tools, only : hermitianize 56 use m_fft_mesh, only : g2ifft 57 use m_fft, only : fourdp_6d, fourdp 58 use m_mpinfo, only : destroy_mpi_enreg, initmpi_seq 59 use m_spacepar, only : hartre 60 use m_rhotoxc, only : rhotoxc 61 use m_dfpt_mkvxc, only : dfpt_mkvxc 62 63 implicit none 64 65 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
SOURCE
1232 subroutine kxc_ADA(Dtset,Cryst,ixc,ngfft,nfft,nspden,rhor,& 1233 & npw,nqibz,qibz,fxc_ADA,gvec,comm,kappa_init,dbg_mode) 1234 1235 !Arguments ------------------------------------ 1236 !scalars 1237 integer,intent(in) :: ixc,nfft,nspden,npw,comm 1238 real(dp),intent(in),optional :: kappa_init 1239 logical,optional,intent(in) :: dbg_mode 1240 type(crystal_t),intent(in) :: Cryst 1241 type(dataset_type),intent(in) :: Dtset 1242 !arrays 1243 integer,intent(in) :: gvec(3,npw),ngfft(18) 1244 integer,intent(in) :: nqibz 1245 real(dp),intent(in) :: rhor(nfft,nspden) 1246 real(dp),intent(in) :: qibz(3,nqibz) 1247 complex(gwpc),intent(out) :: fxc_ADA(npw,npw,nqibz) 1248 1249 !Local variables ------------------------------ 1250 !scalars 1251 integer :: i1,i2,i3,ig,igp,ir,irp,n3xccc,ngfft1,ngfft2,izero !,isp 1252 integer :: ngfft3,nkxc,option,ikxc,ierr,nproc 1253 integer :: nk3xc,igrid,iqbz,my_rank,master,unt_dmp,gmgp_idx 1254 logical :: non_magnetic_xc 1255 real(dp) :: el_temp,enxc,gsqcut,ucvol !,rs,Kx,Kc 1256 real(dp) :: vxcavg,kappa,abs_qpg_sq,abs_qpgp_sq 1257 real(dp) :: difx,dify,difz,inv_kappa_sq 1258 character(len=500) :: msg,fname 1259 type(MPI_type) :: MPI_enreg_seq 1260 type(xcdata_type) :: xcdata 1261 !arrays 1262 real(dp) :: qpg(3),qpgp(3),qphon(3),strsxc(6),q_point(3),dum(0) 1263 real(dp),parameter :: dummyvgeo(3)=zero 1264 real(dp),allocatable :: kxcr(:,:) 1265 real(dp),allocatable :: rhog(:,:),vhartr(:),vxclda(:,:) 1266 real(dp),allocatable :: xccc3d(:),my_rhor(:,:) 1267 real(dp),allocatable :: my_kxcg(:,:) 1268 real(dp),allocatable :: rhotilder(:,:) 1269 complex(gwpc),allocatable :: my_fxc_ADA_ggpq(:,:,:) 1270 complex(gwpc),allocatable :: FT_fxc_ADA_ggpq(:,:,:),dummy(:,:) 1271 real(dp),allocatable :: rvec(:,:),my_fxc_ADA_rrp(:,:) 1272 real(dp) :: rmrp(3),abs_rmrp 1273 integer :: n1,n2,n3,ig_idx_fft(npw) 1274 1275 ! ************************************************************************ 1276 1277 ABI_CHECK(Dtset%nsppol==1,'nsppol/=1 not coded') 1278 ABI_CHECK(nspden==1,'nsppol/=1 not coded') 1279 ABI_CHECK(nfft==PRODUCT(ngfft(1:3)),"mismatch in nfftot") 1280 1281 !Fake MPI_type for the sequential part. 1282 call initmpi_seq(MPI_enreg_seq) 1283 1284 my_rank = xmpi_comm_rank(comm) 1285 nproc = xmpi_comm_size(comm) 1286 master=0 1287 1288 write(msg,'(a,i3)') ' kxc_ADA: calculating exchange-correlation kernel using ixc = ',ixc 1289 call wrtout(std_out,msg,'COLL') 1290 call wrtout(std_out,' kxc_ADA: using smeared density','COLL') 1291 1292 if (.not.present(kappa_init)) then 1293 kappa = 2.1_dp 1294 else 1295 kappa = kappa_init 1296 end if 1297 write(msg,'(a,F10.3)') ' kxc_ADA: inverse smearing length, kappa = ',kappa 1298 call wrtout(std_out,msg,'COLL') 1299 inv_kappa_sq = one/(kappa*kappa) 1300 1301 call xcdata_init(xcdata,dtset=dtset,intxc=0,ixc=ixc,nspden=nspden) 1302 1303 if (ALL(xcdata%xclevel/=(/1,2/))) then 1304 write(msg,'(a,i0)')"Unsupported xclevel = ",xcdata%xclevel 1305 ABI_ERROR(msg) 1306 end if 1307 1308 ngfft1=ngfft(1) 1309 ngfft2=ngfft(2) 1310 ngfft3=ngfft(3) 1311 1312 non_magnetic_xc=(abs(dtset%usepawu)==4.or.dtset%usepawu==14) 1313 1314 if (ixc>=1.and.ixc<11) then ! LDA case 1315 nkxc= 2*min(Dtset%nspden,2)-1 ! 1 or 3 1316 else ! GGA case 1317 nkxc=12*min(Dtset%nspden,2)-5 ! 7 or 19 1318 ABI_CHECK(dtset%xclevel==2,"Functional should be GGA") 1319 ABI_ERROR("GGA functional not implemented for ADA vertex") 1320 end if 1321 1322 ABI_MALLOC(kxcr,(nfft,nkxc)) 1323 1324 !gsqcut and rhog are zeroed because they are not used by rhotoxc if 1<=ixc<=16 and option=0 1325 gsqcut=zero 1326 1327 ABI_MALLOC(rhog,(2,nfft)) 1328 ABI_MALLOC(vhartr,(nfft)) 1329 rhog(:,:)=zero 1330 !MG FIXME this is the 3D core electron density for XC core correction (bohr^-3) 1331 !should implement the non linear core correction 1332 n3xccc=0 1333 ABI_MALLOC(xccc3d,(n3xccc)) 1334 ABI_MALLOC(vxclda,(nfft,nspden)) 1335 1336 option=2 ! 2 for Hxc and kxcr (no paramagnetic part if nspden=1) 1337 qphon(:)=zero 1338 1339 !to be adjusted for the call to rhotoxc 1340 nk3xc=1 1341 1342 !Compute the kernel. 1343 izero=0 1344 1345 !DEBUG print density 1346 if (present(dbg_mode)) then 1347 if (dbg_mode.and.my_rank==master) then 1348 fname = 'xc_ADA_den.xsf' 1349 if (open_file(fname,msg,newunit=unt_dmp,status='unknown',form='formatted') /= 0) then 1350 ABI_ERROR(msg) 1351 end if 1352 call printxsf(ngfft1,ngfft2,ngfft3,rhor(:,1),Cryst%rprimd,(/zero,zero,zero/),& 1353 & Cryst%natom,Cryst%ntypat,Cryst%typat,Cryst%xcart,Cryst%znucl,unt_dmp,0) 1354 close(unt_dmp) 1355 end if 1356 end if 1357 !DEBUG 1358 1359 !Calculate the smeared density 1360 ABI_MALLOC(my_rhor,(nfft,nspden)) 1361 ABI_MALLOC(rhotilder,(nfft,nspden)) 1362 ucvol = Cryst%ucvol 1363 my_rhor = rhor 1364 !do isp = 1,nsppol 1365 !call calc_smeared_density(my_rhor(:,isp),1,rhotilder(:,isp),nfft,ngfft,npw,& 1366 !& gvec,Cryst%gprimd,Cryst%ucvol,MPI_enreg_seq,paral_kgb0,kappa_in=kappa) 1367 !my_rhor(:,isp) = rhotilder(:,isp) 1368 !end do 1369 1370 !DEBUG print smeared density 1371 if (present(dbg_mode)) then 1372 if (dbg_mode.and.my_rank==master) then 1373 fname = 'xc_ADA_smeared_den.xsf' 1374 if (open_file(fname,msg,newunit=unt_dmp,status='unknown',form='formatted') /= 0) then 1375 ABI_ERROR(msg) 1376 end if 1377 call printxsf(ngfft1,ngfft2,ngfft3,my_rhor(:,1),Cryst%rprimd,(/zero,zero,zero/),& 1378 & Cryst%natom,Cryst%ntypat,Cryst%typat,Cryst%xcart,Cryst%znucl,unt_dmp,0) 1379 close(unt_dmp) 1380 end if 1381 end if 1382 !DEBUG 1383 1384 ! Reinitialize the libxc module with the overriden values 1385 if (dtset%ixc<0) then 1386 call libxc_functionals_end() 1387 end if 1388 if (ixc<0) then 1389 el_temp=merge(Dtset%tphysel,Dtset%tsmear,Dtset%tphysel>tol8.and.Dtset%occopt/=3.and.Dtset%occopt/=9) 1390 call libxc_functionals_init(ixc,Dtset%nspden,el_temp=el_temp,xc_tb09_c=Dtset%xc_tb09_c) 1391 end if 1392 1393 call hartre(1,gsqcut,3,izero,MPI_enreg_seq,nfft,ngfft,1,zero,rhog,Cryst%rprimd,dummyvgeo,vhartr) 1394 call rhotoxc(enxc,kxcr,MPI_enreg_seq,nfft,ngfft,& 1395 & dum,0,dum,0,nkxc,nk3xc,non_magnetic_xc,& 1396 & n3xccc,option,my_rhor,Cryst%rprimd,& 1397 & strsxc,1,vxclda,vxcavg,xccc3d,xcdata,vhartr=vhartr) 1398 1399 !Check for extreme (NaN) values 1400 !do ir=1,nfft 1401 !if (isnan(kxcr(ir,1))) kxcr(ir,1) = HUGE(kxcr(ir,1)) 1402 !end do 1403 1404 !DEBUG test with direct way of calculating Kxc 1405 !do i1=1,nfft 1406 !rs = (three/(four_pi*my_rhor(i1,1)))**third 1407 !Kx = 16._dp/27._dp*0.3141592653589793e1_dp*(rs**2)*(-0.4581652_dp) 1408 ! 1409 !Kc = -0.4e1_dp / 0.9e1_dp * 0.3141592654e1_dp * rs ** 4 & 1410 !* (0.207271333333333333333333333333e-1_dp * & 1411 !(-0.177442658629204480000000e3_dp * rs - 0.17565190511219200000000e2_dp & 1412 !* sqrt(rs) - 0.1332650665120000e2_dp * rs ** 2 & 1413 !- 0.51031691247948928000000e2_dp * rs ** (0.3e1_dp / 0.2e1_dp)) & 1414 !* rs ** (-0.3e1_dp / 0.2e1_dp) / (rs + 0.37274400e1_dp * sqrt(rs) & 1415 !+ 0.129352000e2_dp) ** 2 / (-sqrt(rs) - 0.1049800_dp) & 1416 !+ 0.518178333333333333333333333333e-2_dp * rs ** (-0.3e1_dp / 0.2e1_dp) & 1417 !* (0.617071835390850041282140897280e3_dp * sqrt(rs) & 1418 !+ 0.659369347307557491857191871552e5_dp * rs ** 2 + & 1419 !0.700403648491298930017835369562e5_dp * rs ** (0.3e1_dp / 0.2e1_dp) & 1420 !+ 0.398437532951539263722720308167e5_dp * rs ** (0.5e1_dp / 0.2e1_dp) & 1421 !+ 0.368852071032531998953472000000e4_dp * rs ** (0.7e1_dp / 0.2e1_dp) & 1422 !+ 0.5330602660480000e2_dp * rs ** (0.9e1_dp / 0.2e1_dp) & 1423 !+ 0.143783940386264738593799346176e5_dp * rs ** 3 & 1424 !+ 0.124672564145568409213848436081e5_dp * rs & 1425 !+ 0.557398029956167136000000e3_dp * rs ** 4) & 1426 !/ (rs + 0.37274400e1_dp * sqrt(rs) + 0.129352000e2_dp) ** 4 & 1427 !/ (-sqrt(rs) - 0.1049800_dp) ** 2) 1428 !kxcr(i1,1) = Kx + Kc 1429 !end do 1430 !END DEBUG 1431 1432 !DEBUG print Kxc 1433 if (present(dbg_mode)) then 1434 if (dbg_mode.and.my_rank==master) then 1435 fname = 'xc_ADA_Kxc.xsf' 1436 if (open_file(fname,msg,newunit=unt_dmp,status='unknown',form='formatted') /= 0) then 1437 ABI_ERROR(msg) 1438 end if 1439 call printxsf(ngfft1,ngfft2,ngfft3,kxcr(:,1),Cryst%rprimd,(/zero,zero,zero/),& 1440 & Cryst%natom,Cryst%ntypat,Cryst%typat,Cryst%xcart,Cryst%znucl,unt_dmp,0) 1441 close(unt_dmp) 1442 end if 1443 end if 1444 !DEBUG 1445 1446 ABI_FREE(xccc3d) 1447 ABI_FREE(vxclda) 1448 ABI_FREE(vhartr) 1449 1450 ABI_MALLOC(my_kxcg,(2,nfft)) 1451 1452 do ikxc=1,nkxc 1453 call fourdp(1,my_kxcg,kxcr(:,ikxc),-1,MPI_enreg_seq,nfft,1,ngfft,0) 1454 ! kxcg(:,ikxc)=CMPLX(my_kxcg(1,:),my_kxcg(2,:)) 1455 end do 1456 !TODO Check symmetry of kxcg 1457 1458 !set up ADA vertex 1459 ABI_MALLOC(my_fxc_ADA_ggpq,(npw,npw,nqibz)) 1460 my_fxc_ADA_ggpq = czero 1461 !Calculate f_xc(R,R')=(kappa^2/2)K_xc[\tilde{n}](G-G') 1462 !x(1/(kappa^2+|q+G|^2) + 1/(kappa^2+|q+G'|^2 1463 !First get G vectors and indices 1464 1465 ierr=0 1466 do iqbz=1,nqibz 1467 q_point(:) = qibz(:,iqbz) 1468 ! 1469 do ig=1,npw 1470 do igp=1,npw 1471 ! Calculate |q+G| and |q+G'| 1472 qpg(:) = two_pi*MATMUL(Cryst%gprimd,q_point(:)+gvec(:,ig)) 1473 qpgp(:) = two_pi*MATMUL(Cryst%gprimd,q_point(:)+gvec(:,igp)) 1474 abs_qpg_sq = 1.0_dp/(1.0_dp+dot_product(qpg,qpg)*inv_kappa_sq) 1475 abs_qpgp_sq = 1.0_dp/(1.0_dp+dot_product(qpgp,qpgp)*inv_kappa_sq) 1476 1477 gmgp_idx = g2ifft(gvec(:,ig)-gvec(:,igp),ngfft) 1478 if (gmgp_idx>0) then 1479 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) 1480 else 1481 ierr=ierr+1 1482 my_fxc_ADA_ggpq(ig,igp,iqbz) = czero 1483 end if 1484 end do 1485 end do 1486 if (ierr/=0) then 1487 write(msg,'(a,i4,3a)')& 1488 & ' Found ',ierr,' G1-G2 vectors falling outside the FFT box. ',ch10,& 1489 & ' Enlarge the FFT mesh to get rid of this problem. ' 1490 ABI_WARNING(msg) 1491 end if 1492 end do 1493 1494 fxc_ADA = my_fxc_ADA_ggpq 1495 1496 !do iqbz=1,nqibz 1497 !call hermitianize(my_fxc_ADA_ggpq(:,:,iqbz),"All") 1498 !end do 1499 1500 1501 !DEBUG check symmetry 1502 if (.FALSE.) then 1503 ! do iqbz=1,nkptgw 1504 ! do ig=1,npw 1505 ! do igp=ig,npw 1506 ! if (ABS(REAL(fxc_ADA(ig,igp,iqbz))-REAL(fxc_ADA(igp,ig,iqbz)))>tol15.OR.& 1507 ! ABS(AIMAG(fxc_ADA(ig,igp,iqbz))-AIMAG(-fxc_ADA(igp,ig,iqbz)))>tol15) then 1508 ! write(std_out,*) 'Elements:' 1509 ! write(std_out,*) 'fxc_ADA(ig,igp,iqbz):',ig,igp,iqbz,fxc_ADA(ig,igp,iqbz) 1510 ! write(std_out,*) 'fxc_ADA(igp,ig,iqbz):',igp,ig,iqbz,fxc_ADA(igp,ig,iqbz) 1511 ! ABI_ERROR('fxc_ADA not symmetric') 1512 ! end if 1513 ! end do 1514 ! end do 1515 ! end do 1516 1517 ! write(std_out,*)"kxcr(r=0)",kxcr(1,1) 1518 ! write(std_out,*)"my_kxg(G=0)",my_kxcg(:,1) 1519 ! write(std_out,*)"SUM kxcr/nfft ",SUM(kxcr(:,1))/nfft 1520 ! write(std_out,*)"SUM my_kxg ",SUM(kxcg(:,1)) 1521 1522 ! DEBUG Check FT to real space 1523 ! The real-space expression is: 1524 ! f_xc(R,R')=(1/2)(kappa^2/(4*Pi)) 1525 ! \{K_xc[\tilde{n(R)}]+K_xc[\tilde{n(R')}]\} 1526 ! x exp(-kappa|R-R'|)/|R-R'| 1527 ABI_MALLOC(my_fxc_ADA_rrp,(nfft,nfft)) 1528 ABI_MALLOC(FT_fxc_ADA_ggpq,(npw,npw,nqibz)) 1529 ABI_MALLOC(rvec,(3,nfft)) 1530 ABI_MALLOC(dummy,(nfft,nfft)) 1531 my_fxc_ADA_rrp=zero; FT_fxc_ADA_ggpq=czero; dummy=czero; rvec=zero 1532 1533 ! First find coordinates of real-space fft points 1534 igrid = 0 1535 ngfft1 = ngfft(1) 1536 ngfft2 = ngfft(2) 1537 ngfft3 = ngfft(3) 1538 do i3=0,ngfft3-1 1539 difz=dble(i3)/dble(ngfft3) 1540 do i2=0,ngfft2-1 1541 dify=dble(i2)/dble(ngfft2) 1542 do i1=0,ngfft1-1 1543 difx=dble(i1)/dble(ngfft1) 1544 igrid = igrid + 1 1545 rvec(1,igrid)=difx*Cryst%rprimd(1,1)+dify*Cryst%rprimd(1,2)+difz*Cryst%rprimd(1,3) 1546 rvec(2,igrid)=difx*Cryst%rprimd(2,1)+dify*Cryst%rprimd(2,2)+difz*Cryst%rprimd(2,3) 1547 rvec(3,igrid)=difx*Cryst%rprimd(3,1)+dify*Cryst%rprimd(3,2)+difz*Cryst%rprimd(3,3) 1548 end do 1549 end do 1550 end do 1551 if (igrid/=nfft) then 1552 ABI_ERROR('kxc_ADA: igrid not equal to nfft') 1553 end if 1554 1555 ! Construct kernel in real space 1556 do ir=1,nfft 1557 do irp=ir,nfft 1558 rmrp(:) = rvec(:,ir)-rvec(:,irp) 1559 abs_rmrp = sqrt(dot_product(rmrp,rmrp)) 1560 my_fxc_ADA_rrp(ir,irp) = eighth*kappa*kappa*piinv* & 1561 (kxcr(ir,1)+kxcr(irp,1))* & 1562 EXP(-kappa*abs_rmrp)/(abs_rmrp+1.e-3_dp) 1563 ! (a small convergence factor is introduced 1564 ! to avoid a singularity) 1565 my_fxc_ADA_rrp(irp,ir) = my_fxc_ADA_rrp(ir,irp) 1566 end do 1567 end do 1568 1569 ! Find FFT index for all G 1570 n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3) 1571 ! Use the following indexing (N means ngfft of the adequate direction) 1572 ! 0 1 2 3 ... N/2 -(N-1)/2 ... -1 <= kg 1573 ! 1 2 3 4 ....N/2+1 N/2+2 ... N <= index 1574 do ig=1,npw 1575 i1=modulo(gvec(1,ig),n1) 1576 i2=modulo(gvec(2,ig),n2) 1577 i3=modulo(gvec(3,ig),n3) 1578 ig_idx_fft(ig)=i1+1+n1*(i2+n2*i3) 1579 end do 1580 ! FT kernel to reciprocal space for each q 1581 do iqbz=1,nqibz 1582 dummy = CMPLX(my_fxc_ADA_rrp,0.0_dp) 1583 ! Multiply with q-point phase factors exp(-iq.r)*f_xc(r,r')*exp(iq.r') 1584 do ir=1,nfft 1585 do irp=1,nfft 1586 ! Calculate q (variables defined for other purposes 1587 ! are being re-used as dummy variables) 1588 q_point(:) = qibz(:,iqbz) 1589 qpg(:) = two_pi*MATMUL(Cryst%gprimd,q_point(:)) 1590 abs_qpg_sq = dot_product(qpg(:),rvec(:,ir)) 1591 abs_qpgp_sq = dot_product(qpg(:),rvec(:,irp)) 1592 dummy(ir,irp) = EXP(-j_dpc*abs_qpg_sq)* & 1593 & dummy(ir,irp)* & 1594 & EXP(j_dpc*abs_qpgp_sq) 1595 end do 1596 end do 1597 call fourdp_6d(2,dummy,-1,MPI_enreg_seq,nfft,ngfft, 0) 1598 do ig=1,npw 1599 do igp=1,npw 1600 FT_fxc_ADA_ggpq(ig,igp,iqbz) = dummy(ig_idx_fft(ig),ig_idx_fft(igp)) 1601 end do 1602 end do 1603 1604 ! Output 1605 msg='' 1606 if (iqbz<10) write(msg,'(a,i1,a)') './debug_fxc_ADA_q',iqbz,'.dat' 1607 if ((iqbz>9).and.(iqbz<100)) write(msg,'(a,i2,a)') './debug_fxc_ADA_q',iqbz,'.dat' 1608 if ((iqbz>99).and.(iqbz<1000)) write(msg,'(a,i3,a)') './debug_fxc_ADA_q',iqbz,'.dat' 1609 1610 !open(777,file=TRIM(msg),STATUS='REPLACE') 1611 !do igp=1,npw 1612 ! do ig=1,npw 1613 ! write(777,*) ig,igp,REAL(my_fxc_ADA_ggpq(ig,igp,iqbz)),AIMAG(my_fxc_ADA_ggpq(ig,igp,iqbz)), & 1614 ! REAL(FT_fxc_ADA_ggpq(ig,igp,iqbz)),AIMAG(FT_fxc_ADA_ggpq(ig,igp,iqbz)), & 1615 ! ABS(ABS(my_fxc_ADA_ggpq(ig,igp,iqbz))-ABS(FT_fxc_ADA_ggpq(ig,igp,iqbz))) 1616 ! end do 1617 ! write(777,*) '' 1618 !end do 1619 !close(777) 1620 1621 end do ! iqbz 1622 1623 ABI_ERROR('Stopping in kxc_ADA for debugging') 1624 1625 ABI_FREE(rvec) 1626 ABI_FREE(my_fxc_ADA_rrp) 1627 ABI_FREE(FT_fxc_ADA_ggpq) 1628 1629 if (xcdata%xclevel==2) then 1630 ABI_ERROR(" GGA not implemented for kxc_ADA") 1631 end if !xclevel==2 1632 1633 end if ! Debugging section 1634 1635 ! Revert libxc module to the original settings 1636 if (ixc<0) then 1637 call libxc_functionals_end() 1638 end if 1639 if (dtset%ixc<0) then 1640 el_temp=merge(dtset%tphysel,dtset%tsmear,dtset%tphysel>tol8.and.dtset%occopt/=3.and.dtset%occopt/=9) 1641 call libxc_functionals_init(dtset%ixc,dtset%nspden,el_temp=el_temp,xc_tb09_c=dtset%xc_tb09_c) 1642 end if 1643 1644 call destroy_mpi_enreg(MPI_enreg_seq) 1645 ABI_FREE(my_kxcg) 1646 ABI_FREE(my_rhor) 1647 ABI_FREE(rhotilder) 1648 ABI_FREE(rhog) 1649 ABI_FREE(kxcr) 1650 1651 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.
SOURCE
375 subroutine kxc_alda(dtset,ixc,kxcg,mpi_enreg,nfft,ngfft,nspden,option,rhor,rhocut,rprimd) 376 377 !Arguments ------------------------------------------------------------- 378 !scalars 379 integer,intent(in) :: ixc,nfft,nspden,option 380 real(dp),intent(in) :: rhocut 381 type(MPI_type),intent(in) :: mpi_enreg 382 type(dataset_type),intent(in) :: dtset 383 !arrays 384 integer,intent(in) :: ngfft(18) 385 real(dp),intent(in) :: rhor(nfft,2*nspden-1),rprimd(3,3) 386 real(dp),intent(out) :: kxcg(2,nfft,*) 387 388 !Local variables ------------------------------------------------------- 389 !No improved xc quadrature. 390 !No core correction. 391 !Dummy here. 392 !For debugging purposes (see tests below): 393 !integer :: i1,i2,i3,k1,n1,n2,n3 394 !real(dp) :: kx,rho,rhomax,ftest 395 !scalars 396 integer :: ifft,ikxc,isp,n3xccc,ncut,nk3xc,nkxc,optionrhoxc,tim_fourdp 397 logical :: non_magnetic_xc 398 real(dp),parameter :: gsqcut=1._dp 399 real(dp) :: el_temp,enxc,rhocuttot,rhomin,vxcavg 400 character(len=500) :: message 401 type(xcdata_type) :: xcdata 402 !arrays 403 real(dp) :: strsxc(6) 404 real(dp) :: dum(0) 405 real(dp),parameter :: dummyvgeo(3)=zero 406 real(dp),allocatable :: kxcr(:,:),rhog(:,:),rhorcut(:,:),vhartree(:) 407 real(dp),allocatable :: vxc(:,:),xccc3d(:) 408 409 !*********************************************************************** 410 !For debugging purposes (see tests below): 411 412 !ftest(i1,n1,k1) = 0._dp+1._dp*cos(k1*two_pi*float(i1)/float(n1)) 413 414 !*********************************************************************** 415 416 !Check input parameters. 417 418 if (nspden > 2) then 419 message = ' kxc_alda does not work yet for nspden > 2.' 420 ABI_ERROR(message) 421 end if 422 423 !Allocate memory. 424 425 ABI_MALLOC(rhorcut,(nfft,nspden)) 426 ABI_MALLOC(rhog,(2,nfft)) 427 ABI_MALLOC(vhartree,(nfft)) 428 ABI_MALLOC(vxc,(nfft,nspden)) 429 430 call xcdata_init(xcdata,dtset=dtset,intxc=0,ixc=ixc,nspden=nspden) 431 432 non_magnetic_xc=(dtset%usepaw==1.and.mod(abs(dtset%usepawu),10)==4) 433 434 ! Reinitialize the libxc module with the overriden values 435 if (dtset%ixc<0) then 436 call libxc_functionals_end() 437 end if 438 if (ixc<0) then 439 el_temp=merge(dtset%tphysel,dtset%tsmear,dtset%tphysel>tol8.and.dtset%occopt/=3.and.dtset%occopt/=9) 440 call libxc_functionals_init(ixc,dtset%nspden,el_temp=el_temp,xc_tb09_c=dtset%xc_tb09_c) 441 end if 442 443 !to be adjusted for the call to rhotoxc 444 nk3xc=1 445 446 !Cut-off the density. 447 448 rhorcut(:,:) = rhor(:,:) 449 450 do isp = 1,nspden 451 452 rhomin = maxval(rhorcut(:,isp))*rhocut 453 454 ncut = 0 455 rhocuttot = 0._dp 456 457 do ifft = 1,nfft 458 if (rhorcut(ifft,isp) < rhomin) then 459 ncut = ncut+1 460 rhocuttot = rhocuttot+rhorcut(ifft,isp) 461 rhorcut(ifft,isp) = rhomin 462 end if 463 end do 464 465 if (ncut > 0) then 466 write (message,'(a,es10.3,3a,i1,a,i6,a,f6.3,3a,f6.3,a)') & 467 & 'rhocut = ',rhocut,'.',ch10,& 468 & 'For isp = ',isp,' the density was cut-off at ',ncut,' (',100._dp*float(ncut)/float(ifft),'%) grid points.',ch10,& 469 & 'These points account for ',100._dp*rhocuttot/sum(rhor(:,isp)),'% of the total density.' 470 ABI_WARNING(message) 471 end if 472 473 end do 474 475 !Calculate the AL(S)DA kernel in real space. 476 477 rhog(:,:) = 0._dp !We do not need the Hartree potential. 478 tim_fourdp=0 479 480 if ((option == 1).or.((option == 2).and.(nspden == 2))) then 481 482 nkxc = 2*nspden-1 483 n3xccc=0 484 ABI_MALLOC(kxcr,(nfft,nkxc)) 485 ABI_MALLOC(xccc3d,(n3xccc)) 486 487 optionrhoxc = 2 !See rhotoxc.f 488 489 call hartre(1,gsqcut,3,0,mpi_enreg,nfft,ngfft,1,zero,rhog,rprimd,dummyvgeo,vhartree) 490 call rhotoxc(enxc,kxcr,mpi_enreg,nfft,ngfft,dum,0,dum,0,nkxc,nk3xc,non_magnetic_xc,n3xccc,& 491 & optionrhoxc,rhorcut,rprimd,strsxc,1,vxc,vxcavg,xccc3d,xcdata,vhartr=vhartree) 492 493 ! DEBUG 494 ! fx for tests. 495 ! write (std_out,'(a)') ' kxc_alda: Using exchange-only kernel for tests.' 496 ! rhomin = minval(rhor(:,1)) 497 ! rhomax = maxval(rhor(:,1)) 498 ! write (std_out,'(a,es12.5,a,es12.5)') ' kxc_alda: rhomin = ',rhomin,' rhomax = ',rhomax 499 ! write (std_out,'(a)') ' kxc_alda: loping below 0.2*rhomax.' 500 ! kx = (3._dp/4._dp)*((3._dp/pi)**(1._dp/3._dp)) 501 ! do ifft = 1,nfft 502 ! rho = rhor(ifft,1) 503 ! rho = max(rho,0.2_dp*rhomax) 504 ! kxcr(ifft,1) = -(4._dp/9._dp)*kx*(rho**(-2._dp/3._dp)) 505 ! write (std_out,'(i4,a,es12.5)') ifft,': ',kxcr(ifft,1) 506 ! end do 507 ! write (std_out,'(a,es12.5)') 'kxcrmin: ',minval(kxcr(:,1)) 508 ! write (std_out,'(a,es12.5)') 'kxcrmax: ',maxval(kxcr(:,1)) 509 ! ENDDEBUG 510 511 ! DEBUG 512 ! test kernel. 513 ! write(std_out,'(a)') ' kxc_alda: Using test kernel for tests.' 514 ! n1 = ngfft(1) ; n2 = ngfft(2) ; n3 = ngfft(3) 515 ! do i3 = 0,n3-1 516 ! do i2 = 0,n2-1 517 ! do i1 = 0,n1-1 518 ! ifft = i1+n1*(i2+n2*i3)+1 519 ! kxcr(ifft,1) = ftest(i1,n1,1)*ftest(i2,n2,2)*ftest(i3,n3,3) 520 ! end do 521 ! end do 522 ! end do 523 ! ENDDEBUG 524 525 ! Calculate the Fourier transform of the AL(S)DA kernel. 526 527 if (option == 1) then 528 do ikxc = 1,nkxc 529 call fourdp(1,kxcg(:,:,ikxc),kxcr(:,ikxc),-1,mpi_enreg,nfft,1,ngfft,tim_fourdp) 530 end do 531 else 532 call fourdp(1,kxcg(:,:,1),kxcr(:,2),-1,mpi_enreg,nfft,1,ngfft,tim_fourdp) 533 end if 534 535 else if ((option == 2).and.(nspden == 1)) then 536 537 nkxc = 2 538 n3xccc=0 539 ABI_MALLOC(kxcr,(nfft,nkxc)) 540 ABI_MALLOC(xccc3d,(n3xccc)) 541 542 optionrhoxc = -2 !See rhotoxc.f 543 544 call hartre(1,gsqcut,3,0,mpi_enreg,nfft,ngfft,1,zero,rhog,rprimd,dummyvgeo,vhartree) 545 call rhotoxc(enxc,kxcr,mpi_enreg,nfft,ngfft,dum,0,dum,0,nkxc,nk3xc,non_magnetic_xc,n3xccc,& 546 & optionrhoxc,rhorcut,rprimd,strsxc,1,vxc,vxcavg,xccc3d,xcdata,vhartr=vhartree) 547 548 kxcr(:,2) = 0.5_dp*kxcr(:,2) 549 550 ! Calculate the Fourier transform of the up-down channel of the AL(S)DA kernel. 551 552 call fourdp(1,kxcg(:,:,1),kxcr(:,2),-1,mpi_enreg,nfft,1,ngfft,tim_fourdp) 553 554 else 555 write (message,'(4a,i0)')' Invalid option = ',option 556 ABI_ERROR(message) 557 end if 558 559 !DEBUG 560 !write(std_out,*) ' kxc_alda: Exc = ',enxc 561 !write(std_out,*) ' kxc_alda: <Vxc> = ',vxcavg 562 !ENDDEBUG 563 564 ! Revert libxc module to the original settings 565 if (ixc<0) then 566 call libxc_functionals_end() 567 end if 568 if (dtset%ixc<0) then 569 el_temp=merge(dtset%tphysel,dtset%tsmear,dtset%tphysel>tol8.and.dtset%occopt/=3.and.dtset%occopt/=9) 570 call libxc_functionals_init(dtset%ixc,dtset%nspden,el_temp=el_temp,xc_tb09_c=dtset%xc_tb09_c) 571 end if 572 573 !Free memory. 574 ABI_FREE(rhorcut) 575 ABI_FREE(rhog) 576 ABI_FREE(vhartree) 577 ABI_FREE(vxc) 578 ABI_FREE(kxcr) 579 ABI_FREE(xccc3d) 580 581 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
SOURCE
985 subroutine kxc_driver(Dtset,Cryst,ixc,ngfft,nfft_tot,nspden,rhor,npw,dim_kxcg,kxcg,gvec,comm,dbg_mode) 986 987 !Arguments ------------------------------------ 988 !scalars 989 integer,intent(in) :: ixc,npw,nfft_tot,nspden,dim_kxcg,comm 990 logical,optional,intent(in) :: dbg_mode 991 type(crystal_t),intent(in) :: Cryst 992 type(dataset_type),intent(in) :: Dtset 993 !arrays 994 integer,intent(in) :: gvec(3,npw),ngfft(18) 995 real(dp),intent(in) :: rhor(nfft_tot,nspden) 996 complex(gwpc),intent(out) :: kxcg(nfft_tot,dim_kxcg) 997 998 !Local variables ------------------------------ 999 !scalars 1000 integer :: cplex,i1,i2,i3,ig,igp,iq,ir,n3xccc,ngfft1,ngfft2,izero 1001 integer :: ngfft3,nkxc,option,ikxc,nk3xc,my_rank,master,unt_dmp 1002 logical :: non_magnetic_xc 1003 real(dp) :: el_temp,enxc,expo,gpqx,gpqy,gpqz,gsqcut,vxcavg 1004 character(len=500) :: msg,fname 1005 type(xcdata_type) :: xcdata 1006 type(MPI_type) :: MPI_enreg_seq 1007 !arrays 1008 real(dp) :: qphon(3),strsxc(6),dum(0) 1009 real(dp),parameter :: dummyvgeo(3)=zero 1010 real(dp),allocatable :: kxcpw_g(:,:),kxcr(:,:),phas(:,:,:) 1011 real(dp),allocatable :: rhog(:,:),vhartr(:),kxcpw_r(:,:),vxclda(:,:) 1012 real(dp),allocatable :: xccc3d(:),xx(:,:) 1013 real(dp),allocatable :: my_kxcg(:,:) 1014 1015 !************************************************************************ 1016 1017 ABI_CHECK(Dtset%nsppol==1,'nsppol/=1 not coded') 1018 ABI_CHECK(Dtset%nspden==1,'nsppol/=1 not coded') 1019 ABI_CHECK(nfft_tot==PRODUCT(ngfft(1:3)),"mismatch in nfftot") 1020 1021 !Fake MPI_type for the sequential part. 1022 call initmpi_seq(MPI_enreg_seq) 1023 call init_distribfft_seq(MPI_enreg_seq%distribfft,'c',ngfft(2),ngfft(3),'all') 1024 my_rank = xmpi_comm_rank(comm) 1025 master=0 1026 1027 write(msg,'(a,i3)') ' kxc_driver: calculating exchange-correlation kernel using ixc = ',ixc 1028 call wrtout(std_out,msg,'COLL') 1029 1030 call xcdata_init(xcdata,dtset=Dtset,intxc=0,ixc=ixc,nspden=nspden) 1031 1032 if (ALL(xcdata%xclevel/=(/1,2/))) then 1033 write(msg,'(a,i0)')"Unsupported xclevel = ",xcdata%xclevel 1034 ABI_ERROR(msg) 1035 end if 1036 1037 ngfft1=ngfft(1) 1038 ngfft2=ngfft(2) 1039 ngfft3=ngfft(3) 1040 1041 non_magnetic_xc=(dtset%usepaw==1.and.mod(abs(dtset%usepawu),10)==4) 1042 1043 if (ixc>=1.and.ixc<11) then ! LDA case 1044 nkxc= 2*min(nspden,2)-1 ! 1 or 3 1045 else ! GGA case 1046 nkxc=12*min(nspden,2)-5 ! 7 or 19 1047 ABI_CHECK(dtset%xclevel==2,"Functional should be GGA") 1048 ABI_ERROR("GGA functional not tested") 1049 end if 1050 1051 ABI_MALLOC(kxcr,(nfft_tot,nkxc)) 1052 1053 !gsqcut and rhog are zeroed because they are not used by rhotoxc if 1<=ixc<=16 and option=0 1054 gsqcut=zero 1055 1056 ABI_MALLOC(rhog,(2,nfft_tot)) 1057 ABI_MALLOC(vhartr,(nfft_tot)) 1058 rhog(:,:)=zero 1059 !MG FIXME this is the 3D core electron density for XC core correction (bohr^-3) 1060 !should implement the non linear core correction 1061 n3xccc=0 1062 ABI_MALLOC(xccc3d,(n3xccc)) 1063 ABI_MALLOC(vxclda,(nfft_tot,nspden)) 1064 1065 option=2 ! 2 for Hxc and kxcr (no paramagnetic part if nspden=1) 1066 qphon =zero 1067 1068 !to be adjusted for the call to rhotoxc 1069 nk3xc=1 1070 izero=0 1071 1072 ! Reinitialize the libxc module with the overriden values 1073 if (dtset%ixc<0) then 1074 call libxc_functionals_end() 1075 end if 1076 if (ixc<0) then 1077 el_temp=merge(Dtset%tphysel,Dtset%tsmear,Dtset%tphysel>tol8.and.Dtset%occopt/=3.and.Dtset%occopt/=9) 1078 call libxc_functionals_init(ixc,Dtset%nspden,el_temp=el_temp,xc_tb09_c=Dtset%xc_tb09_c) 1079 end if 1080 1081 call hartre(1,gsqcut,3,izero,MPI_enreg_seq,nfft_tot,ngfft,1,zero,rhog,Cryst%rprimd,dummyvgeo,vhartr) 1082 1083 !Compute the kernel. 1084 call rhotoxc(enxc,kxcr,MPI_enreg_seq,nfft_tot,ngfft,& 1085 & dum,0,dum,0,nkxc,nk3xc,non_magnetic_xc,& 1086 & n3xccc,option,rhor,Cryst%rprimd,& 1087 & strsxc,1,vxclda,vxcavg,xccc3d,xcdata,vhartr=vhartr) 1088 1089 ABI_FREE(rhog) 1090 ABI_FREE(vhartr) 1091 !DEBUG print Kxc 1092 if (present(dbg_mode)) then 1093 if (dbg_mode .and. my_rank==master) then 1094 fname = 'xc_Kxc.xsf' 1095 if (open_file(fname,msg,newunit=unt_dmp,status='unknown',form='formatted') /= 0) then 1096 ABI_ERROR(msg) 1097 end if 1098 call printxsf(ngfft1,ngfft2,ngfft3,kxcr(:,1),Cryst%rprimd,(/zero,zero,zero/),& 1099 & Cryst%natom,Cryst%ntypat,Cryst%typat,Cryst%xcart,Cryst%znucl,unt_dmp,0) 1100 close(unt_dmp) 1101 end if 1102 end if 1103 !DEBUG 1104 1105 ABI_FREE(xccc3d) 1106 ABI_FREE(vxclda) 1107 1108 ABI_MALLOC(my_kxcg,(2,nfft_tot)) 1109 1110 do ikxc=1,nkxc 1111 call fourdp(1,my_kxcg,kxcr(:,ikxc),-1,MPI_enreg_seq,nfft_tot,1,ngfft,0) 1112 kxcg(:,ikxc)=CMPLX(my_kxcg(1,:),my_kxcg(2,:)) 1113 end do 1114 1115 !write(std_out,*)"kxcr(r=0)",kxcr(1,1) 1116 !write(std_out,*)"my_kxg(G=0)",my_kxcg(:,1) 1117 !write(std_out,*)"SUM kxcr/nfft_tot ",SUM(kxcr(:,1))/nfft_tot 1118 !write(std_out,*)"SUM my_kxg ",SUM(kxcg(:,1)) 1119 1120 ABI_FREE(my_kxcg) 1121 1122 !MG this part is never executed, but one should use dfpt_mkvxc for the GGA kernel. 1123 if (xcdata%xclevel==2) then 1124 ABI_ERROR("check GGA implementation") 1125 cplex=2 1126 ABI_MALLOC(phas,(cplex*nfft_tot,npw,nspden)) 1127 ABI_MALLOC(kxcpw_r,(cplex*nfft_tot,nspden)) 1128 ABI_MALLOC(xx,(3,nfft_tot)) 1129 ABI_MALLOC(kxcpw_g,(2,nfft_tot)) 1130 1131 kxcg = czero 1132 1133 ! find the coordinates for all r in the FFT grid 1134 ir=0 1135 do i3=1,ngfft3 1136 do i2=1,ngfft2 1137 do i1=1,ngfft1 1138 ir=ir+1 1139 xx(1,ir)=dble((i1-1))/ngfft1 1140 xx(2,ir)=dble((i2-1))/ngfft2 1141 xx(3,ir)=dble((i3-1))/ngfft3 1142 end do 1143 end do 1144 end do 1145 1146 do iq=1,1 1147 1148 ! Calculate at once exp(i(G+q).r), for all possible q,G,r 1149 do ig=1,npw 1150 gpqx=dble(gvec(1,ig)) 1151 gpqy=dble(gvec(2,ig)) 1152 gpqz=dble(gvec(3,ig)) 1153 do ir=1,nfft_tot 1154 expo=gpqx*xx(1,ir)+gpqy*xx(2,ir)+gpqz*xx(3,ir) 1155 phas(2*ir-1,ig,1)= cos(two_pi*expo) 1156 phas(2*ir,ig,1) = sin(two_pi*expo) 1157 end do 1158 end do 1159 1160 ! 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} $ 1161 1162 do igp=1,npw 1163 1164 kxcpw_r(:,:)=zero 1165 1166 call dfpt_mkvxc(cplex,ixc,kxcr,MPI_enreg_seq,nfft_tot,ngfft,dum,0,dum,0,nkxc,non_magnetic_xc,& 1167 & nspden,n3xccc,option,qphon(:),phas(:,igp,:),Cryst%rprimd,1,kxcpw_r,xccc3d) 1168 1169 ! FFT the first index to --> to G space 1170 call fourdp(cplex,kxcpw_g(:,:),kxcpw_r(:,1),-1,MPI_enreg_seq,nfft_tot,1,ngfft,0) 1171 1172 ! kxcg(:,igp,iq)=CMPLX(kxcpw_g(1,igfft(:)),kxcpw_g(2,igfft(:))) 1173 ! kxcg(:,igp)=CMPLX(kxcpw_g(1,igfft(:)),kxcpw_g(2,igfft(:))) 1174 1175 end do ! igp 1176 end do ! iq 1177 1178 ABI_FREE(phas) 1179 ABI_FREE(kxcpw_r) 1180 ABI_FREE(xx) 1181 ABI_FREE(kxcpw_g) 1182 end if !xclevel==2 1183 1184 ! Revert libxc module to the original settings 1185 if (ixc<0) then 1186 call libxc_functionals_end() 1187 end if 1188 if (dtset%ixc<0) then 1189 el_temp=merge(dtset%tphysel,dtset%tsmear,dtset%tphysel>tol8.and.dtset%occopt/=3.and.dtset%occopt/=9) 1190 call libxc_functionals_init(dtset%ixc,dtset%nspden,el_temp=el_temp,xc_tb09_c=dtset%xc_tb09_c) 1191 end if 1192 1193 call destroy_mpi_enreg(MPI_enreg_seq) 1194 ABI_FREE(kxcr) 1195 1196 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
SOURCE
836 subroutine kxc_eok(ixceok,kxcg,mpi_enreg,nfft,ngfft,nspden,rhor,rhocut) 837 838 !Arguments ------------------------------------------------------------- 839 !scalars 840 integer,intent(in) :: ixceok,nfft,nspden 841 real(dp),intent(in) :: rhocut 842 type(MPI_type),intent(in) :: mpi_enreg 843 !arrays 844 integer,intent(in) :: ngfft(18) 845 real(dp),intent(in) :: rhor(nfft,2*nspden-1) 846 real(dp),intent(out) :: kxcg(2,nfft,2*nspden-1) 847 848 !Local variables ------------------------------------------------------- 849 !Maximum value allowed for rs. 850 !scalars 851 integer :: ifft,ikxc,ncut,nkxc,nlop,tim_fourdp 852 real(dp),parameter :: rslim=50._dp,dummyvgeo(3)=zero 853 real(dp) :: a2,a3,a4,rho,rhocuttot,rhomin,rs 854 character(len=500) :: message 855 !arrays 856 real(dp),allocatable :: kxcr(:,:) 857 858 !*********************************************************************** 859 860 !Check input parameters. 861 862 if (nspden > 1) then 863 message = ' kxc_eok does not work yet for nspden > 1.' 864 ABI_ERROR(message) 865 end if 866 867 !Values of a2, a3 and a4 for case 1 868 a2 = 0.0_dp 869 a3 = 0.0_dp 870 a4 = 0.0_dp 871 872 select case (ixceok) 873 case (1) 874 a2 = -0.51887_dp 875 a3 = 4.9359d-03 876 a4 = -5.9603d-05 877 case (2) 878 a2 = -0.50044_dp 879 a3 = 4.9653d-03 880 a4 = -3.3660d-05 881 case default 882 message = ' kxc_eok: ixceok /= 1 (linear EOK) or 2 (non-linear EOK).' 883 ABI_ERROR(message) 884 end select 885 886 !Allocate memory. 887 888 nkxc = 2*nspden-1 889 890 ABI_MALLOC(kxcr,(nfft,nkxc)) 891 892 !Calculate the energy optimized kernel in real space. 893 894 nlop = 0 895 896 rhomin = rhocut*maxval(rhor(:,:)) 897 898 ncut = 0 899 rhocuttot = 0._dp 900 901 do ifft = 1,nfft 902 903 rho = rhor(ifft,1) 904 905 if (rho < rhomin) then 906 ncut = ncut+1 907 rhocuttot = rhocuttot+rho 908 rho = rhomin 909 end if 910 911 rs = (3._dp/(4._dp*pi*rho))**(1._dp/3._dp) 912 913 if (rs > rslim) then 914 rs = rslim 915 nlop = nlop+1 916 end if 917 918 kxcr(ifft,1) = a2*rs**2+a3*rs**3+a4*rs**4 919 920 end do 921 922 if (ncut > 0) then 923 write (message,'(a,es10.3,3a,i1,a,i6,a,f6.3,3a,f6.3,a)') & 924 & 'rhocut = ',rhocut,'.',ch10,& 925 & 'For isp = ',1,' the density was cut-off at ',ncut,' (',100._dp*float(ncut)/float(ifft),'%) grid points.',ch10,& 926 & 'These points account for ',100._dp*rhocuttot/sum(rhor(:,1)),'% of the total density.' 927 ABI_WARNING(message) 928 end if 929 930 if (nlop > 0) then 931 write (message,'(a,f6.2,a,i6,a,f6.3,a)') & 932 & 'rs still exceeds ',rslim,' Bohr at ',nlop,' (',100._dp*float(nlop)/float(ifft),'%) grid points (after cut-off).' 933 ABI_WARNING(message) 934 end if 935 936 !Calculate the Fourier transform of the energy optimized kernel. 937 tim_fourdp=0 938 do ikxc = 1,nkxc 939 call fourdp(1,kxcg(:,:,ikxc),kxcr(:,ikxc),-1,mpi_enreg,nfft,1,ngfft,tim_fourdp) 940 end do 941 942 !Free memory. 943 944 ABI_FREE(kxcr) 945 946 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.
SOURCE
187 subroutine kxc_local(ispxc,kg_diel,kxc,kxcg,nfft,ngfft,npwdiel,nspden,option) 188 189 !Arguments ------------------------------------------------------------- 190 !scalars 191 integer,intent(in) :: ispxc,nfft,npwdiel,nspden,option 192 !arrays 193 integer,intent(in) :: kg_diel(3,npwdiel),ngfft(18) 194 real(dp),intent(in) :: kxcg(2,nfft) 195 real(dp),intent(out) :: kxc(2,npwdiel,nspden,npwdiel,nspden) 196 197 !Local variables ------------------------------------------------------- 198 !For debbuging purposes: 199 !real(dp) :: c1,c2,c3 200 !scalars 201 integer :: i1,i2,i3,ifft,ipw1,ipw2,ipwstart,isp1,isp2,j1,j2,j3,k1,k2,k3,n1,n2 202 integer :: n3 203 logical :: ok 204 character(len=500) :: message 205 206 !*********************************************************************** 207 208 !Check input parameters. 209 210 if (nspden > 2) then 211 message =' kxc_local does not work yet for nspden > 2.' 212 ABI_ERROR(message) 213 end if 214 215 isp1 = 1 216 isp2 = 1 217 ok = .true. 218 if (nspden == 1) then 219 select case (ispxc) 220 case (1) 221 isp1 = 1 222 isp2 = 1 223 case default 224 ok = .false. 225 end select 226 else 227 select case (ispxc) 228 case (1) 229 isp1 = 1 230 isp2 = 1 231 case (2) 232 isp1 = 1 233 isp2 = 2 234 case (3) 235 isp1 = 2 236 isp2 = 2 237 case default 238 ok = .false. 239 end select 240 end if 241 242 if (.not.ok) then 243 write (message,'(2(a,i0))')' The input ispxc = ',ispxc,' is not compatible with nspden = ',nspden 244 ABI_BUG(message) 245 end if 246 247 if (option == 0) then 248 ipwstart = 2 249 kxc(:,1,isp1,:,isp2) = 0._dp 250 kxc(:,:,isp1,1,isp2) = 0._dp 251 else 252 ipwstart = 1 253 end if 254 255 !Calculate the xc matrix. 256 257 n1 = ngfft(1) ; n2 = ngfft(2) ; n3 = ngfft(3) 258 259 do ipw2 = ipwstart,npwdiel 260 261 j1 = kg_diel(1,ipw2) ; j2 = kg_diel(2,ipw2) ; j3 = kg_diel(3,ipw2) 262 263 ! Fill the diagonal. 264 265 kxc(:,ipw2,isp1,ipw2,isp2) = kxcg(:,1) 266 267 ! Fill the off-diagonal elements. 268 269 do ipw1 = ipw2+1,npwdiel 270 271 i1 = kg_diel(1,ipw1) ; i2 = kg_diel(2,ipw1) ; i3 = kg_diel(3,ipw1) 272 273 ! Compute the difference between G vectors. 274 ! The use of two mod calls handles both i1-j1 >= n1 AND i1-j1 < 0. 275 276 k1 = mod(n1+mod(i1-j1,n1),n1) 277 k2 = mod(n2+mod(i2-j2,n2),n2) 278 k3 = mod(n3+mod(i3-j3,n3),n3) 279 280 ifft = k1+n1*(k2+n2*k3)+1 281 282 kxc(1,ipw1,isp1,ipw2,isp2) = kxcg(1,ifft) 283 kxc(2,ipw1,isp1,ipw2,isp2) = kxcg(2,ifft) 284 285 kxc(1,ipw2,isp1,ipw1,isp2) = kxcg(1,ifft) 286 kxc(2,ipw2,isp1,ipw1,isp2) = -kxcg(2,ifft) 287 288 end do 289 290 end do 291 292 !If needed, copy the up-down to the down-up spin channel. 293 294 if (ispxc == 2) then 295 do ipw2 = 1,npwdiel 296 do ipw1 = 1,npwdiel 297 kxc(1,ipw2,isp2,ipw1,isp1) = kxc(1,ipw1,isp1,ipw2,isp2) 298 kxc(2,ipw2,isp2,ipw1,isp1) = -kxc(2,ipw1,isp1,ipw2,isp2) 299 end do 300 end do 301 end if 302 303 !DEBUG 304 !See kxc_alda.f, "test kernel" DEBUG section. 305 !do ipw2 = 1,npwdiel 306 !j1 = kg_diel(1,ipw2) ; j2 = kg_diel(2,ipw2) ; j3 = kg_diel(3,ipw2) 307 !do ipw1 = ipw2+1,npwdiel 308 !i1 = kg_diel(1,ipw1) ; i2 = kg_diel(2,ipw1) ; i3 = kg_diel(3,ipw1) 309 !k1 = mod(n1+mod(i1-j1,n1),n1) 310 !k2 = mod(n2+mod(i2-j2,n2),n2) 311 !k3 = mod(n3+mod(i3-j3,n3),n3) 312 !ifft = k1+n1*(k2+n2*k3)+1 313 !c1 = 0._dp ; c2 = 0._dp ; c3 = 0._dp 314 !if (i1-j1 == 0) c1 = c1+0.0_dp 315 !if (i2-j2 == 0) c2 = c2+0.0_dp 316 !if (i3-j3 == 0) c3 = c3+0.0_dp 317 !if (i1-j1 == 1) c1 = c1+0.5_dp 318 !if (i2-j2 == 2) c2 = c2+0.5_dp 319 !if (i3-j3 == 3) c3 = c3+0.5_dp 320 !if (i1-j1 == -1) c1 = c1+0.5_dp 321 !if (i2-j2 == -2) c2 = c2+0.5_dp 322 !if (i3-j3 == -3) c3 = c3+0.5_dp 323 !if ((abs(kxcg(1,ifft)-c1*c2*c3) > tol10).or.(abs(kxcg(2,ifft)) > tol10)) then 324 !write (std_out,*) ' i1 i2 i3 ifft: ',i1,i2,i3,ifft 325 !write (std_out,*) ' exp.: ',c1*c2*c3,' got: ',kxcg(:,ifft) 326 !end if 327 !end do 328 !end do 329 !ENDDEBUG 330 331 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.
SOURCE
623 subroutine kxc_pgg(gmet,kg,khxcg,npw,rcut_coulomb,susmat,ucvol) 624 625 !Arguments ------------------------------------ 626 !scalars 627 integer,intent(in) :: npw 628 real(dp),intent(in) :: rcut_coulomb,ucvol 629 !arrays 630 integer,intent(in) :: kg(3,npw) 631 real(dp),intent(in) :: gmet(3,3),susmat(2,npw,npw) 632 real(dp),intent(out) :: khxcg(2,npw,npw) 633 634 !Local variables------------------------------- 635 !scalars 636 integer :: i1,i2,i3,ig,ii,ikg11,ikg12,ikg13,ikg21,ikg22,ikg23,ipw1,ipw2 637 integer :: j1,j2,j3,jg,jj 638 real(dp),parameter :: diffgsq=1.d-2 639 real(dp) :: kg_red1,kg_red2,kg_red3,gsquar,tpisq 640 !arrays 641 integer :: kgmax(3) 642 integer,allocatable :: index_g_inv(:,:,:),jgarr(:) 643 real(dp),allocatable :: gsq(:),sumg(:),vcoul(:) 644 645 ! ************************************************************************* 646 647 !DEBUG 648 !write(std_out,*) '%kxc_pgg: enter' 649 !write(std_out,*) 'npw=',npw 650 !ENDDEBUG 651 652 !tpisq is (2 Pi) **2: 653 tpisq=(two_pi)**2 654 655 kgmax(:)=0 656 do ipw1=1,npw 657 do jj=1,3 658 kgmax(jj)=max( kg(jj,ipw1), kgmax(jj) ) 659 end do 660 end do 661 662 !DEBUG 663 !write(std_out,*) 'kgmax:',kgmax(1:3) 664 !ENDDEBUG 665 666 !Perform allocations 667 ABI_MALLOC(index_g_inv,(-2*kgmax(1):2*kgmax(1),-2*kgmax(2):2*kgmax(2),-2*kgmax(3):2*kgmax(3))) 668 ABI_MALLOC(jgarr,(npw)) 669 ABI_MALLOC(gsq,(npw)) 670 ABI_MALLOC(sumg,(2)) 671 ABI_MALLOC(vcoul,(npw)) 672 673 !DEBUG 674 !write(std_out,*) '%kxc_pg: creating plane wave index and coulomb potential' 675 !ENDDEBUG 676 index_g_inv(:,:,:)=0 677 do ipw1=1,npw 678 index_g_inv(kg(1,ipw1),kg(2,ipw1),kg(3,ipw1))=ipw1 679 680 ! DEBUG 681 ! write(std_out,'(i5,2x,3i3,2x,i4)') ipw1,kg(1,ipw1),kg(2,ipw1),kg(3,ipw1) 682 ! ENDDEBUG 683 684 kg_red1=dble(kg(1,ipw1)) 685 kg_red2=dble(kg(2,ipw1)) 686 kg_red3=dble(kg(3,ipw1)) 687 gsquar=tpisq*(gmet(1,1)*kg_red1**2+gmet(2,2)*kg_red2**2+gmet(3,3)*kg_red3**2 & 688 & +2.0_dp*( (gmet(1,2)*kg_red2+gmet(1,3)*kg_red3)* kg_red1 + & 689 & gmet(2,3)*kg_red2*kg_red3) ) 690 ! Distinguish G=0 from other elements 691 if(gsquar > 1.0d-12)then 692 vcoul(ipw1)=four_pi/gsquar*(1._dp-cos(sqrt(gsquar)*rcut_coulomb)) 693 else 694 vcoul(ipw1)=four_pi*0.5_dp*rcut_coulomb**2 695 end if 696 697 end do 698 699 !DEBUG 700 !write(std_out,*) '%kxc_pg: starting convolution integral' 701 !ENDDEBUG 702 703 !loop over G1,G2 components of the density matrix 704 do ipw2=1,npw 705 ikg21=kg(1,ipw2) 706 ikg22=kg(2,ipw2) 707 ikg23=kg(3,ipw2) 708 709 do ii=1,npw 710 j1=ikg21-kg(1,ii) 711 j2=ikg22-kg(2,ii) 712 j3=ikg23-kg(3,ii) 713 jgarr(ii)=index_g_inv(j1,j2,j3) 714 end do 715 716 do ipw1=1,ipw2 717 ikg11=kg(1,ipw1) 718 ikg12=kg(2,ipw1) 719 ikg13=kg(3,ipw1) 720 721 ! do the convolution integral 722 sumg(:)=0._dp 723 do ii=1,npw 724 725 if( jgarr(ii) /= 0 ) then 726 727 i1=ikg11-kg(1,ii) 728 i2=ikg12-kg(2,ii) 729 i3=ikg13-kg(3,ii) 730 731 ! j1=ikg21-kg(1,ii) 732 ! j2=ikg22-kg(2,ii) 733 ! j3=ikg23-kg(3,ii) 734 735 ig=index_g_inv(i1,i2,i3) 736 ! jg=index_g_inv(j1,j2,j3) 737 738 if( ig /= 0 ) then 739 740 jg=jgarr(ii) 741 742 ! DEBUG 743 ! write(std_out,'(i5,2x,3i3,1x,3i3,2x,2i4)') ii,i1,i2,i3,& 744 ! & kg(1,jg),kg(2,jg),kg(3,jg),& 745 ! & ig,jg 746 ! ENDDEBUG 747 748 sumg(1)=sumg(1)+susmat(1,ig,jg)*vcoul(ii) 749 sumg(2)=sumg(2)+susmat(2,ig,jg)*vcoul(ii) 750 751 end if 752 753 end if 754 755 end do 756 khxcg(:,ipw1,ipw2)=-sumg(:)*ucvol 757 758 ! DEBUG 759 ! if(ipw1==ipw2) write(std_out,'(2i4,2(1x,es14.6))') ipw1,ipw2,khxcg(1,ipw1,ipw1),vcoul(ipw1) 760 ! write(std_out,'(2i4,3(1x,es14.6))') ipw1,ipw2,khxcg(1:2,ipw1,ipw2),vcoul(ipw1) 761 ! ENDDEBUG 762 763 end do 764 end do 765 766 !DEBUG 767 !verify hermiticity, note: ipw1 loop above must end at npw 768 !write(std_out,*) '%kxc_pgg: check hermiticity of pgg kernel' 769 !do ipw2=1,npw,max(2,npw/10) 770 !do ipw1=ipw2,npw,max(2,npw/10) 771 !write(std_out,'(2i4,2(1x,es14.6))') ipw1,ipw2,& 772 !& khxcg(1,ipw1,ipw2)-khxcg(1,ipw2,ipw1),& 773 !& khxcg(2,ipw1,ipw2)+khxcg(2,ipw2,ipw1) 774 !end do 775 !end do 776 !ENDDEBUG 777 778 !Impose hermiticity 779 write(std_out,*) '%kxc_pg: imposing hermiticity' 780 do ipw2=1,npw 781 do ipw1=ipw2+1,npw 782 khxcg(1,ipw1,ipw2)= khxcg(1,ipw2,ipw1) 783 khxcg(2,ipw1,ipw2)=-khxcg(2,ipw2,ipw1) 784 end do 785 end do 786 787 !DEBUG 788 !write(std_out,'(a10,2(1x,es20.12))') '+ex_pgg? ', 0.5_dp*khxcg(1,1,1)/ucvol 789 !ENDDEBUG 790 791 ABI_FREE(index_g_inv) 792 ABI_FREE(jgarr) 793 ABI_FREE(gsq) 794 ABI_FREE(sumg) 795 ABI_FREE(vcoul) 796 797 !DEBUG 798 !write(std_out,*) '%kxc_pgg: done' 799 !ENDDEBUG 800 801 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.
SOURCE
104 subroutine kxc_rpa(gsq,krpa,npw,option,rcut_coulomb) 105 106 !Arguments ------------------------------------------------------------- 107 !scalars 108 integer,intent(in) :: npw,option 109 real(dp),intent(in) :: rcut_coulomb 110 !arrays 111 real(dp),intent(in) :: gsq(npw) 112 real(dp),intent(out) :: krpa(npw) 113 114 !Local variables ------------------------------------------------------- 115 !scalars 116 integer :: ipw 117 118 !*********************************************************************** 119 120 if (option == 0) then 121 122 ! Compute the bare Hartree kernel. 123 124 do ipw = 1,npw 125 126 if (gsq(ipw) > tol12) then 127 krpa(ipw) = four_pi/gsq(ipw) 128 else 129 krpa(ipw) = 0._dp 130 end if 131 132 end do 133 134 else 135 136 ! Compute the Hartree kernel with a cut-off in real space beyond rcut_coulomb: 137 138 do ipw = 1,npw 139 140 if (gsq(ipw) > tol12) then 141 krpa(ipw) = (four_pi/gsq(ipw))*(1._dp-cos(sqrt(gsq(ipw))*rcut_coulomb)) 142 else 143 krpa(ipw) = two_pi*rcut_coulomb**2 144 end if 145 146 end do 147 148 end if 149 150 end subroutine kxc_rpa