TABLE OF CONTENTS


ABINIT/m_ppmodel [ Modules ]

[ Top ] [ Modules ]

NAME

 m_ppmodel

FUNCTION

  Module containing the definition of the ppmodel_t used to deal with
  the plasmonpole technique. Methods to operate on the object are also provided.

COPYRIGHT

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

PARENTS

CHILDREN

SOURCE

22 #if defined HAVE_CONFIG_H
23 #include "config.h"
24 #endif
25 
26 #include "abi_common.h"
27 
28 MODULE m_ppmodel
29 
30  use defs_basis
31  use defs_abitypes
32  use m_errors
33  use m_abicore
34  use m_array
35  use m_linalg_interfaces
36 
37  use m_fstrings,       only : sjoin, itoa
38  use m_hide_lapack,    only : xhegv
39  use m_gwdefs,         only : GW_Q0_DEFAULT, czero_gw
40  use m_crystal,        only : crystal_t
41  use m_bz_mesh,        only : kmesh_t, get_bz_item
42  use m_gsphere,        only : gsphere_t
43  use m_vcoul,          only : vcoul_t, cmod_qpg
44  use m_fft_mesh,       only : g2ifft
45  use m_fft,            only : fourdp
46  use m_mpinfo,         only : destroy_mpi_enreg, initmpi_seq
47 
48  implicit none
49 
50  private

m_ppmodel/calc_sig_ppm [ Functions ]

[ Top ] [ m_ppmodel ] [ Functions ]

NAME

 calc_sig_ppm

FUNCTION

  Calculate the contribution to self-energy operator using a plasmon-pole model.

INPUTS

  nomega=Number of frequencies.
  nspinor=Number of spinorial components.
  npwc=Number of G vectors in the plasmon pole.
  npwx=number of G vectors in rhotwgp, i.e. no. of G-vectors for the exchange part.
  theta_mu_minus_e0i= $\theta(\mu-\epsilon_{k-q,b1,s}), defines if the state is occupied or not
  zcut=Small imaginary part to avoid the divergence. (see related input variable)
  omegame0i(nomega)=Frequencies used to evaluate \Sigma_c ($\omega$ - $\epsilon_i)$
  otq(npwc,dm2_otq)=Plasmon pole parameters for this q-point.
  PPm<ppmodel_t>=structure gathering info on the Plasmon-pole technique.
     %model=type plasmon pole model
     %dm2_botsq= 1 if model==3, =npwc if model== 4, 1 for all the other cases
     %dm2_otq= 1 if model==3, =1    if model== 4, 1 for all the other cases
     %dm_eig=npwc if model=3, 0 otherwise
  botsq(npwc,dm2_botsq)=Plasmon pole parameters for this q-point.
  eig(dm_eig,dm_eig)=The eigvectors of the symmetrized inverse dielectric matrix for this q point
   (first index for G, second index for bands)
  rhotwgp(npwx)=oscillator matrix elements divided by |q+G| i.e.
    $\frac{\langle b1 k-q s | e^{-i(q+G)r | b2 k s \rangle}{|q+G|}$

OUTPUT

  sigcme(nomega) (to be described), only relevant if ppm3 or ppm4
  ket(npwc,nomega):
  === model==1,2 ====
    ket(G,omega) = Sum_G2       conjg(rhotw(G)) * Omega(G,G2) * rhotw(G2)
                            ---------------------------------------------------
                             2 omegatw(G,G2) (omega-E_i + omegatw(G,G2)(2f-1))

PARENTS

      calc_sigc_me

CHILDREN

SOURCE

2361 subroutine calc_sig_ppm(PPm,nspinor,npwc,nomega,rhotwgp,botsq,otq,&
2362 & omegame0i,zcut,theta_mu_minus_e0i,eig,npwx,ket,sigcme)
2363 
2364 
2365 !This section has been created automatically by the script Abilint (TD).
2366 !Do not modify the following lines by hand.
2367 #undef ABI_FUNC
2368 #define ABI_FUNC 'calc_sig_ppm'
2369 !End of the abilint section
2370 
2371  implicit none
2372 
2373 !Arguments ------------------------------------
2374 !scalars
2375  integer,intent(in) :: nomega,npwc,npwx,nspinor
2376  real(dp),intent(in) :: theta_mu_minus_e0i,zcut
2377  type(ppmodel_t),intent(in) :: PPm
2378 !arrays
2379  real(dp),intent(in) :: omegame0i(nomega)
2380  complex(gwpc),intent(in) :: botsq(npwc,PPm%dm2_botsq),eig(PPm%dm_eig,PPm%dm_eig),otq(npwc,PPm%dm2_otq)
2381  complex(gwpc),intent(in) :: rhotwgp(npwx*nspinor)
2382  complex(gwpc),intent(inout) :: ket(npwc*nspinor,nomega)
2383  complex(gwpc),intent(out) :: sigcme(nomega)
2384 
2385 !Local variables-------------------------------
2386 !scalars
2387  integer :: ig,igp,ii,ios,ispinor,spadc,spadx
2388  real(dp) :: den,ff,inv_den,omegame0i_io,otw,twofm1,twofm1_zcut
2389  complex(gwpc) :: ct,num,numf,rhotwgdp_igp
2390  logical :: fully_occupied,totally_empty
2391  !character(len=500) :: msg
2392 !arrays
2393  complex(gwpc),allocatable :: rhotwgdpcc(:)
2394 
2395 !*************************************************************************
2396 
2397  SELECT CASE (PPm%model)
2398  CASE (PPM_GODBY_NEEDS, PPM_HYBERTSEN_LOUIE)
2399    fully_occupied =(ABS(theta_mu_minus_e0i-one)<0.001)
2400    totally_empty  =(ABS(theta_mu_minus_e0i    )<0.001)
2401 
2402    do ispinor=1,nspinor
2403      spadx=(ispinor-1)*npwx; spadc=(ispinor-1)*npwc
2404 
2405      if (.not.totally_empty) then
2406        ! \Bomega^2_{G1G2}/\omegat_{G1G2} M_{G1,G2}. \theta(\mu-e_s) / (\omega+\omegat_{G1G2}-e_s-i\delta)
2407        twofm1_zcut=zcut
2408 !$omp parallel do private(omegame0i_io, rhotwgdp_igp, otw, num, den)
2409        do ios=1,nomega
2410          omegame0i_io=omegame0i(ios)
2411          do igp=1,npwc
2412            rhotwgdp_igp=rhotwgp(spadx+igp)
2413            do ig=1,npwc
2414              otw=DBLE(otq(ig,igp)) !in principle otw -> otw - ieta
2415              num = botsq(ig,igp)*rhotwgdp_igp
2416              den = omegame0i_io+otw
2417              if (den**2>zcut**2) then
2418                ket(spadc+ig,ios)=ket(spadc+ig,ios) + num/(den*otw)*theta_mu_minus_e0i
2419              else
2420                ket(spadc+ig,ios)=ket(spadc+ig,ios) + &
2421 &                num*CMPLX(den,twofm1_zcut)/((den**2+twofm1_zcut**2)*otw)*theta_mu_minus_e0i
2422              end if
2423            end do !ig
2424          end do !igp
2425        end do !ios
2426      end if !not totally empty
2427 
2428      if (.not.(fully_occupied)) then
2429        ! \Bomega^2_{G1G2}/\omegat_{G1G2} M_{G1,G2}. \theta(e_s-\mu) / (\omega-\omegat_{G1G2}-e_s+i\delta)
2430        twofm1_zcut=-zcut
2431 !$omp parallel do private(omegame0i_io, rhotwgdp_igp, otw, num, den)
2432        do ios=1,nomega
2433          omegame0i_io=omegame0i(ios)
2434          do igp=1,npwc
2435            rhotwgdp_igp=rhotwgp(spadx+igp)
2436            do ig=1,npwc
2437              otw=DBLE(otq(ig,igp)) !in principle otw -> otw - ieta
2438              num = botsq(ig,igp)*rhotwgdp_igp
2439              den=omegame0i_io-otw
2440              if (den**2>zcut**2) then
2441                ket(spadc+ig,ios)=ket(spadc+ig,ios) + num/(den*otw)*(one-theta_mu_minus_e0i)
2442              else
2443                ket(spadc+ig,ios)=ket(spadc+ig,ios) + &
2444 &                num*CMPLX(den,twofm1_zcut)/((den**2+twofm1_zcut**2)*otw)*(one-theta_mu_minus_e0i)
2445              end if
2446            end do !ig
2447          end do !igp
2448          !
2449        end do !ios
2450      end if !not fully occupied
2451 
2452    end do !ispinor
2453 
2454    ket=ket*half
2455 
2456  CASE (PPM_LINDEN_HORSH,PPM_ENGEL_FARID)
2457    ABI_CHECK(nspinor == 1, "nspinor/=1 not allowed")
2458 
2459    ! rho-twiddle(G) is formed, introduce rhotwgdpcc, for speed reason
2460    ABI_MALLOC(rhotwgdpcc,(npwx))
2461 
2462    ff=theta_mu_minus_e0i      ! occupation number f (include poles if ...)
2463    twofm1=two*ff-one          ! 2f-1
2464    twofm1_zcut=twofm1*zcut
2465    rhotwgdpcc(:)=CONJG(rhotwgp(:))
2466 
2467    do ios=1,nomega
2468      omegame0i_io=omegame0i(ios)
2469      ct=czero_gw
2470      do ii=1,npwc ! Loop over the DM bands
2471        num=czero_gw
2472 
2473        SELECT CASE (PPm%model)
2474        CASE (PPM_LINDEN_HORSH)
2475          ! Calculate \beta (eq. 106 pag 47)
2476          do ig=1,npwc
2477            num=num+rhotwgdpcc(ig)*eig(ig,ii)
2478          end do
2479          numf=num*CONJG(num) !MG this means that we cannot do SCGW
2480          numf=numf*botsq(ii,1)
2481 
2482        CASE (PPM_ENGEL_FARID)
2483          do ig=1,npwc
2484            num=num+rhotwgdpcc(ig)*botsq(ig,ii)
2485          end do
2486          numf=num*CONJG(num) !MG this means that we cannot do SCGW
2487 
2488        CASE DEFAULT
2489          MSG_ERROR("Wrong PPm%model")
2490        END SELECT
2491 
2492        !numf=num*CONJG(num) !MG this means that we cannot do SCGW
2493        !if (PPm%model==3) numf=numf*botsq(ii,1)
2494 
2495        otw=DBLE(otq(ii,1)) ! in principle otw -> otw - ieta
2496        den=omegame0i_io+otw*twofm1
2497 
2498        if (den**2>zcut**2) then
2499          inv_den=one/den
2500          ct=ct+numf*inv_den
2501        else
2502          inv_den=one/((den**2+twofm1_zcut**2))
2503          ct=ct+numf*CMPLX(den,twofm1_zcut)*inv_den
2504        end if
2505 
2506      end do !ii DM bands
2507      sigcme(ios)=ct*half
2508    end do !ios
2509    ABI_FREE(rhotwgdpcc)
2510 
2511  CASE DEFAULT
2512    MSG_BUG(sjoin('Wrong PPm%model:',itoa(PPm%model)))
2513  END SELECT
2514 
2515 end subroutine calc_sig_ppm

m_ppmodel/cppm1par [ Functions ]

[ Top ] [ m_ppmodel ] [ Functions ]

NAME

 cppm1par

FUNCTION

 Calculate the plasmon-pole parameters big-omega-twiddle-squared and omega-twiddle from
 epsilon-twiddle^-1 calculated for nomega (usually 2) frequencies omega=0 and omega=iE0.

INPUTS

  epsm1(npwc,npwc,nomega)=dielectric matrix at nomega frequencies.
  npwc=number of plane waves
  nomega=number of frequencies (usually 2)
  omega(nomega)=frequencies
  omegaplasma=input variable or Drude plasma frequency

OUTPUT

  bigomegatwsq(npwc,npwc)=parameter of the plasmon-pole model (see gwa.pdf file)
  omegatw(npwc,npwc)=parameter of the plasmon-pole model (see gwa.pdf file)

TODO

  Calculation can be done in place.

PARENTS

      m_ppmodel

CHILDREN

SOURCE

1339 subroutine cppm1par(npwc,nomega,omega,omegaplasma,epsm1,omegatw,bigomegatwsq)
1340 
1341 
1342 !This section has been created automatically by the script Abilint (TD).
1343 !Do not modify the following lines by hand.
1344 #undef ABI_FUNC
1345 #define ABI_FUNC 'cppm1par'
1346 !End of the abilint section
1347 
1348  implicit none
1349 
1350 !Arguments ------------------------------------
1351 !scalars
1352  integer,intent(in) :: nomega,npwc
1353  real(dp),intent(in) :: omegaplasma
1354 !arrays
1355  complex(gwpc),intent(in) :: epsm1(npwc,npwc,nomega)
1356  complex(dpc),intent(in) :: omega(nomega)
1357  complex(gwpc),intent(out) :: bigomegatwsq(npwc,npwc)
1358  complex(gwpc),intent(out) :: omegatw(npwc,npwc)
1359 
1360 !Local variables-------------------------------
1361 !scalars
1362  integer :: ig,igp,io,io0,ioe0
1363  real(dp) :: e0,minomega
1364  character(len=500) :: msg
1365  complex(gwpc) :: AA,omegatwsq,diff,ratio
1366  complex(gwpc) :: epsm1_io0,epsm1_ioe0
1367 
1368 ! *************************************************************************
1369 
1370  DBG_ENTER("COLL")
1371  !
1372  ! === Find omega=0 and omega=imag (closest to omegaplasma) to fit the ppm parameters ===
1373  minomega=1.0d-3; io0=0
1374  do io=1,nomega
1375    if (ABS(omega(io))<minomega) then
1376      io0=io; minomega=ABS(omega(io))
1377    end if
1378  end do
1379  ABI_CHECK(io0/=0,"omega=0 not found")
1380 
1381  minomega=1.0d-3; e0=200.0; ioe0=0
1382  do io=1,nomega
1383    if (REAL(omega(io))<minomega.and.AIMAG(omega(io))>minomega) then
1384      if (ABS(AIMAG(omega(io))-omegaplasma)<ABS(e0-omegaplasma)) then
1385        ioe0=io; e0=AIMAG(omega(io))
1386      end if
1387    end if
1388  end do
1389 
1390  write(msg,'(a,f9.4,a)')' Imaginary frequency for fit located at: ',e0*Ha_eV,' [eV] '
1391  call wrtout(std_out,msg,'COLL')
1392  ABI_CHECK(ioe0/=0,"Imaginary omega not found")
1393  !
1394  ! ================================================================
1395  ! === Calculate plasmon-pole A parameter A=epsilon^-1(0)-delta ===
1396  ! ================================================================
1397  do ig=1,npwc
1398    do igp=1,npwc
1399      epsm1_io0  = epsm1(ig,igp,io0)
1400      epsm1_ioe0 = epsm1(ig,igp,ioe0)
1401 
1402      AA=epsm1_io0
1403      if (ig==igp) AA=AA-one
1404 
1405      ! === Calculate plasmon-pole omega-twiddle-square parameter ===
1406      ! XG201009 Strangely, the next formula does not work with gcc43-debug
1407      ! omegatwsq=(AA/(epsm1_io0-epsm1_ioe0)-one)*e0**2
1408      ! This seems to be due to precision issue at the level of division by a complex whose norm squared
1409      ! is below the smallest representable number.
1410      ! After many trials, I have decided to shift the difference by a small number ... well, not so small ...
1411      ! for numerical issues
1412      diff=epsm1_io0-epsm1_ioe0
1413      diff=diff+cmplx(tol10,tol10)
1414      ratio=AA/diff
1415      omegatwsq=(ratio-cone)*e0**2
1416      !
1417      ! If omega-twiddle-squared is negative,set omega-twiddle-squared to 1.0 (a reasonable way of treating
1418      ! such terms, in which epsilon**-1 was originally increasing along this part of the imaginary axis)
1419      ! (note: originally these terms were ignored in Sigma; this was changed on 6 March 1990.)
1420 
1421      if (REAL(omegatwsq)<=zero) omegatwsq=one
1422      !
1423      ! Get omega-twiddle
1424      ! * Neglect the imag part (if any) in omega-twiddle-squared
1425      omegatw(ig,igp)=SQRT(REAL(omegatwsq))
1426      !
1427      ! Get big-omega-twiddle-squared=-omega-twiddle-squared AA
1428      bigomegatwsq(ig,igp)=-AA*omegatw(ig,igp)**2
1429 
1430    end do !igp
1431  end do !ig
1432 
1433  write(msg,'(2a,f15.12,2a,2i5,a)')ch10,&
1434 &  ' cppm1par : omega twiddle minval [eV]  = ',MINVAL(ABS(omegatw))*Ha_eV,ch10,&
1435 &  '            omega twiddle min location = ',MINLOC(ABS(omegatw)),ch10
1436  call wrtout(std_out,msg,'COLL')
1437 
1438  DBG_EXIT("COLL")
1439 
1440 end subroutine cppm1par

m_ppmodel/cppm2par [ Functions ]

[ Top ] [ m_ppmodel ] [ Functions ]

NAME

 cppm2par

FUNCTION

  Calculate plasmon-pole parameters of the Hybertsen and Louie model (PRB 34, 5390 (1986) [[cite:Hybertsen1986]])

INPUTS

  qpt(3)=The coordinates of the q-point in the IBZ.
  epsm1(npwc,npwc)=symmetrized inverse dielectric (static limit is used)
  gmet(3,3)=metric in reciprocal space
  ngfftf(18)=contain all needed information about the 3D fine FFT, see ~abinit/doc/variables/vargs.htm#ngfft
  npwc=number of plane waves in epsm1
  rhor(nfftf)=charge density on the real space FFT grid
  nfftf= total number of points in the fine FFT mesh  (for this processor)
  invalid_freq: what to do when PPM omega is found negative or imaginary
     0) drop it (default as specified in Hybersen-Louie original GW paper)
     1) set to 1 hartree
     2) set to infinity

OUTPUT

  bigomegatwsq(npwc,npwc)= squared bare plasma frequencies
   \Omega^2_{G1 G2}(q) = 4\pi \frac {(q+G1).(q+G2)}/{|q+G1|^2} n(G1-G2)
  omegatw(npwc,npwc)= plasmon frequencies \tilde\omega_{G1 G2}(q) where:
  \tilde\omega^2_{G1 G2}(q) =
    \frac {\Omega^2_{G1 G2}(q)} {\delta_{G1 G2}-\tilde\epsilon^{-1}_{G1 G2} (q, \omega=0)}

PARENTS

      m_ppmodel

CHILDREN

SOURCE

1479 subroutine cppm2par(qpt,npwc,epsm1,ngfftf,gvec,gprimd,rhor,nfftf,gmet,bigomegatwsq,omegatw,invalid_freq)
1480 
1481 
1482 !This section has been created automatically by the script Abilint (TD).
1483 !Do not modify the following lines by hand.
1484 #undef ABI_FUNC
1485 #define ABI_FUNC 'cppm2par'
1486 !End of the abilint section
1487 
1488  implicit none
1489 
1490 !Arguments ------------------------------------
1491 !scalars
1492  integer,intent(in) :: npwc,nfftf,invalid_freq
1493 !arrays
1494  integer,intent(in) :: gvec(3,npwc)
1495  integer,intent(in) :: ngfftf(18)
1496  real(dp),intent(in) :: qpt(3),gmet(3,3),gprimd(3,3)
1497  real(dp),intent(in) :: rhor(nfftf)
1498  complex(gwpc),intent(in) :: epsm1(npwc,npwc)
1499  complex(gwpc),intent(out) :: bigomegatwsq(npwc,npwc)
1500  complex(gwpc),intent(out) :: omegatw(npwc,npwc)
1501 
1502 !Local variables-------------------------------
1503 !scalars
1504  integer,parameter :: paral_kgb0=0
1505  integer :: ig,igp,nimwp,ngfft1,ngfft2,ngfft3,gmgp_idx,ierr
1506  real(dp) :: lambda,phi,AA
1507  logical,parameter :: use_symmetrized=.TRUE.,check_imppf=.FALSE.
1508  character(len=500) :: msg
1509  type(MPI_type) :: MPI_enreg_seq
1510 !arrays
1511  real(dp) :: qlist(3,1)
1512  real(dp),allocatable :: tmp_rhor(:),qratio(:,:),qplusg(:),rhog_dp(:,:)
1513  complex(gwpc),allocatable :: omegatwsq(:,:)
1514  complex(gwpc),allocatable :: rhog(:),rhogg(:,:),temp(:,:)  !MG these should be double precision TODO
1515 !*************************************************************************
1516 
1517  call initmpi_seq(MPI_enreg_seq)
1518  call init_distribfft_seq(MPI_enreg_seq%distribfft,'c',ngfftf(2),ngfftf(3),'all')
1519  !
1520  ! === Calculate qratio(npwec,npvec) = (q+G).(q+Gp)/|q+G|^2 ===
1521  ABI_STAT_MALLOC(qratio,(npwc,npwc), ierr)
1522  ABI_CHECK(ierr==0, "out-of-memory in qratio")
1523 
1524  call cqratio(npwc,gvec,qpt,gmet,gprimd,qratio)
1525  !
1526  ! === Compute the density in G space rhor(R)--> rhog(G) ===
1527  ABI_MALLOC(rhog_dp,(2,nfftf))
1528  ABI_MALLOC(rhog,(nfftf))
1529  ngfft1=ngfftf(1)
1530  ngfft2=ngfftf(2)
1531  ngfft3=ngfftf(3)
1532 
1533  ABI_MALLOC(tmp_rhor,(nfftf))
1534  tmp_rhor=rhor ! To avoid having to use intent(inout).
1535  call fourdp(1,rhog_dp,tmp_rhor,-1,MPI_enreg_seq,nfftf,ngfftf,paral_kgb0,0)
1536  ABI_FREE(tmp_rhor)
1537 
1538  rhog(1:nfftf)=CMPLX(rhog_dp(1,1:nfftf),rhog_dp(2,1:nfftf))
1539  !
1540  ! Calculate the FFT index of each (G-Gp) vector and assign
1541  ! the value of the correspondent density simultaneously
1542  ABI_STAT_MALLOC(rhogg,(npwc,npwc), ierr)
1543  ABI_CHECK(ierr==0, "out of memory rhogg")
1544 
1545  ierr=0
1546  do ig=1,npwc
1547    do igp=1,npwc
1548      gmgp_idx = g2ifft(gvec(:,ig)-gvec(:,igp),ngfftf)
1549      if (gmgp_idx/=0) then
1550        rhogg(ig,igp)=rhog(gmgp_idx)
1551      else
1552        ierr=ierr+1
1553        rhogg(ig,igp)=czero
1554      end if
1555    end do
1556  end do
1557 
1558  if (ierr/=0) then
1559    write(msg,'(a,i4,3a)')&
1560 &   'Found ',ierr,' G1-G2 vectors falling outside the FFT box. ',ch10,&
1561 &   'Enlarge the FFT mesh to get rid of this problem. '
1562    MSG_WARNING(msg)
1563  end if
1564 
1565  rhogg=four_pi*rhogg
1566  ABI_FREE(rhog_dp)
1567  ABI_FREE(rhog)
1568  !
1569  ! Calculate GPP parameters
1570  ! unsymmetrized epsm1 -> epsm1=|q+Gp|/|q+G|*epsm1
1571  ABI_MALLOC(qplusg,(npwc))
1572  ABI_MALLOC(temp,(npwc,npwc))
1573  ABI_STAT_MALLOC(omegatwsq,(npwc,npwc), ierr)
1574  ABI_CHECK(ierr==0, 'out of memory in omegatwsq')
1575 
1576  temp = -epsm1(:,:)
1577  !
1578  ! RS still not obvious for me whether one shall use the symmetrized inverse DM or the unsymmetrized one
1579  ! the default here is to use the symmetrized one, I must discuss this with XG
1580  !
1581  ! MG it turns out that using the symmetrized inverse DM in the plasmon-pole
1582  ! equations give the same results for the squared plasmon frequencies omegatwsq while the
1583  ! squared bare plasma frequencies bigomegatwsq related to the symmetrized dielectric matrix
1584  ! are obtained multipling by |q+G1|/|q+G2|
1585  !
1586  if (.not.use_symmetrized) then
1587    qlist(:,1) = qpt
1588    call cmod_qpg(1,1,qlist,npwc,gvec,gprimd,qplusg) !MG TODO here take care of small q
1589    do ig=1,npwc
1590      do igp=1,npwc
1591        temp(ig,igp)=qplusg(igp)/qplusg(ig)*temp(ig,igp)
1592      end do
1593    end do
1594  end if
1595 
1596  nimwp=0
1597  do ig=1,npwc
1598    temp(ig,ig)=temp(ig,ig)+one
1599    do igp=1,npwc
1600      bigomegatwsq(ig,igp) = rhogg(ig,igp)*qratio(ig,igp)
1601      omegatwsq(ig,igp)=bigomegatwsq(ig,igp)/temp(ig,igp)
1602      !
1603      ! Set to an arbitrary value the omegawsq which become negative or imaginary
1604      ! in principle these correspond to cases where the imaginary part of epsm1 does not have
1605      ! a well defined peak. The imaginary part of epsm1 in these cases oscillates  with a small amplitude
1606      ! since the amplitude A_GGpr=-pi/2*bigomegatwsq/omegatw,
1607      ! it follows that bigomegatwsq shall be set to zero for these cases
1608      if ( REAL(omegatwsq(ig,igp))<= tol12 .or. AIMAG(omegatwsq(ig,igp))**2*tol12> REAL(omegatwsq(ig,igp))**2) then
1609        nimwp=nimwp+1
1610 
1611        if ( invalid_freq == 1 ) then
1612        ! set omegatwsq to 1 hartree
1613          omegatwsq(ig,igp)=cone
1614          AA = epsm1(ig,igp)
1615          if ( ig == igp ) AA = AA - one
1616          omegatw(ig,igp)=SQRT(REAL(omegatwsq(ig,igp)))
1617          bigomegatwsq(ig,igp)=-AA*omegatw(ig,igp)**2
1618        elseif ( invalid_freq == 2 ) then
1619        ! set omegatwsq to infinity
1620          omegatwsq(ig,igp)=cone/tol6
1621          AA = epsm1(ig,igp)
1622          if ( ig == igp ) AA = AA - one
1623          omegatw(ig,igp)=SQRT(REAL(omegatwsq(ig,igp)))
1624          bigomegatwsq(ig,igp)=-AA*omegatw(ig,igp)**2
1625        else
1626        ! simply ignore all cases of omegatw with imaginary values
1627          bigomegatwsq(ig,igp)=(0.,0.)
1628          omegatw(ig,igp)=(ten,0.)
1629        end if
1630        if (check_imppf) then
1631          write(msg,'(a,2i8)')' Imaginary plasmon frequency at : ',ig,igp
1632          call wrtout(std_out,msg,'COLL')
1633        end if
1634      else
1635        ! this part has been added to deal with systems without inversion symmetry
1636        ! this new implementation gives the same results as the previous one if
1637        ! omegatwsq is a pure real number and has the advantage of being an improved
1638        ! approach for systems without an inversion center.
1639        lambda=ABS(omegatwsq(ig,igp))
1640        phi=ATAN(AIMAG(omegatwsq(ig,igp))/REAL(omegatwsq(ig,igp)))
1641        omegatw(ig,igp)=SQRT(lambda/COS(phi))
1642        bigomegatwsq(ig,igp)=bigomegatwsq(ig,igp)*(1.-(0.,1.)*TAN(phi))
1643        ! Uncomment the following line and comment the previous to restore the old version.
1644        !omegatw(ig,igp)=sqrt(real(omegatwsq(ig,igp)))
1645      end if
1646    end do
1647  end do
1648 
1649  write(msg,'(a,3f12.6,a,i8,a,i8)')&
1650 &  ' at q-point : ',qpt,&
1651 &  ' number of imaginary plasmonpole frequencies = ',nimwp,' / ',npwc**2
1652  call wrtout(std_out,msg,'COLL')
1653 
1654  write(msg,'(2a,f12.8,2a,3i5)')ch10,&
1655 &  ' cppm2par : omega twiddle minval [eV]  = ',MINVAL(ABS(omegatw))*Ha_eV,ch10,&
1656 &  '            omega twiddle min location = ',MINLOC(ABS(omegatw))
1657  call wrtout(std_out,msg,'COLL')
1658 
1659  call destroy_mpi_enreg(MPI_enreg_seq)
1660 
1661  ABI_FREE(omegatwsq)
1662  ABI_FREE(rhogg)
1663  ABI_FREE(temp)
1664  ABI_FREE(qplusg)
1665  ABI_FREE(qratio)
1666 
1667 end subroutine cppm2par

m_ppmodel/cppm3par [ Functions ]

[ Top ] [ m_ppmodel ] [ Functions ]

NAME

 cppm3par

FUNCTION

 Calculate the plasmon-pole parameters using the von Linden-Horsh model (PRB 37, 8351, 1988) [[cite:vonderLinden1988]]
 (see also Pag 22 of Quasiparticle Calculations in Solids [[cite:Aulbur2001]].

INPUTS

 epsm1(npwc,npwc))= symmetrized inverse dielectric
 ngfftf(18)=contain all needed information about 3D fine FFT, see ~abinit/doc/variables/vargs.htm#ngfft
 npwc=number of plane waves in epsm1
 qratio=(q+G1).(q+G2)/(|q+G1|.|q+G2|)
 rhor(nfftf)=charge density on the real space FFT grid
 nfftf=number of points in the FFT grid (for this processor)
 gvec(3,npwc)= G vectors in reduced coordinates

OUTPUT

  omegatw(npwc,npwc)= plasmon pole positions
  bigomegatwsq(npwc,npwc)=(E_{q,ii}^{-1}-1)*omegatw
   where E^{-1} is the eigenvalue of the inverse dielectric matrix
  eigtot(npwc,npwc)=the eigvectors of the symmetrized inverse dielectric matrix
   (first index for G, second index for bands)

PARENTS

      m_ppmodel

CHILDREN

SOURCE

1703 subroutine cppm3par(qpt,npwc,epsm1,ngfftf,gvec,gprimd,rhor,nfftf,bigomegatwsq,omegatw,eigtot)
1704 
1705 
1706 !This section has been created automatically by the script Abilint (TD).
1707 !Do not modify the following lines by hand.
1708 #undef ABI_FUNC
1709 #define ABI_FUNC 'cppm3par'
1710 !End of the abilint section
1711 
1712  implicit none
1713 
1714 !Arguments ------------------------------------
1715 !scalars
1716  integer,intent(in) :: nfftf,npwc
1717 !arrays
1718  integer,intent(in) :: gvec(3,npwc),ngfftf(18)
1719  real(dp),intent(in) :: qpt(3),gprimd(3,3),rhor(nfftf)
1720  complex(gwpc),intent(in) :: epsm1(npwc,npwc)
1721  complex(gwpc),intent(out) :: bigomegatwsq(npwc,1),eigtot(npwc,npwc)
1722  complex(gwpc),intent(out) :: omegatw(npwc)
1723 
1724 !Local variables-------------------------------
1725 !TODO these should be dp
1726 !scalars
1727  integer,parameter :: paral_kgb0=0
1728  integer :: idx,ierr,ig,igp,ii,jj,ngfft1,ngfft2,ngfft3,gmgp_idx
1729  real(dp) :: num,qpg_dot_qpgp
1730  complex(dpc) :: conjg_eig
1731  logical :: qiszero
1732  character(len=500) :: msg
1733  type(MPI_type) :: MPI_enreg_seq
1734 !arrays
1735  real(dp) :: b1(3),b2(3),b3(3),gppq(3),gpq(3),qlist(3,1)
1736  real(dp),allocatable :: eigval(:),qplusg(:),rhog_dp(:,:),zhpev2(:),tmp_rhor(:)
1737  complex(dpc),allocatable :: eigvec(:,:),matr(:),mm(:,:),rhog(:),rhogg(:,:)
1738  complex(dpc),allocatable :: zhpev1(:),zz(:)
1739 
1740 !*************************************************************************
1741 
1742  ! Fake MPI_type for the sequential part.
1743  call initmpi_seq(MPI_enreg_seq)
1744  call init_distribfft_seq(MPI_enreg_seq%distribfft,'c',ngfftf(2),ngfftf(3),'all')
1745 
1746  qiszero = (ALL(ABS(qpt)<1.0e-3))
1747 
1748  b1=two_pi*gprimd(:,1)
1749  b2=two_pi*gprimd(:,2)
1750  b3=two_pi*gprimd(:,3)
1751 
1752  ngfft1=ngfftf(1)
1753  ngfft2=ngfftf(2)
1754  ngfft3=ngfftf(3)
1755 
1756  ABI_MALLOC(rhog_dp,(2,nfftf))
1757  ABI_MALLOC(rhog,(nfftf))
1758 
1759  ABI_STAT_MALLOC(rhogg,(npwc,npwc), ierr)
1760  ABI_CHECK(ierr==0, "out of memory in rhogg")
1761  !
1762  ! === Compute the density in G space rhog(r)--> rho(G) ===
1763  ! FIXME this has to be fixed, rho(G) should be passed instead of doing FFT for each q
1764 
1765  ABI_MALLOC(tmp_rhor,(nfftf))
1766  tmp_rhor=rhor ! To avoid having to use intent(inout).
1767  call fourdp(1,rhog_dp,tmp_rhor,-1,MPI_enreg_seq,nfftf,ngfftf,paral_kgb0,0)
1768  ABI_FREE(tmp_rhor)
1769 
1770  rhog(1:nfftf)=CMPLX(rhog_dp(1,1:nfftf),rhog_dp(2,1:nfftf))
1771  !
1772  ! Calculate the FFT index of each (G-Gp) vector and assign the value
1773  ! of the correspondent density simultaneously
1774  ierr=0
1775  do ig=1,npwc
1776    do igp=1,npwc
1777      gmgp_idx = g2ifft(gvec(:,ig)-gvec(:,igp),ngfftf)
1778      if (gmgp_idx/=0) then
1779        rhogg(ig,igp)=rhog(gmgp_idx)
1780      else
1781        ierr=ierr+1
1782        rhogg(ig,igp)=czero
1783      end if
1784    end do
1785  end do
1786 
1787  if (ierr/=0) then
1788    write(msg,'(a,i4,3a)')&
1789 &   'Found ',ierr,' G1-G2 vectors falling outside the FFT box. ',ch10,&
1790 &   'Enlarge the FFT mesh to get rid of this problem. '
1791    MSG_WARNING(msg)
1792  end if
1793  !
1794  ! mm(G,Gp) = (q+G) \cdot (q+Gp) n(G-Gp)
1795  ABI_STAT_MALLOC(mm,(npwc,npwc), ierr)
1796  ABI_CHECK(ierr==0, 'out of memory in mm')
1797 
1798  do ig=1,npwc
1799    if (qiszero) then
1800      ! To be discussed with Riad, here we should use the small q
1801      ! to be consistent and consider the limit q-->0
1802      gpq(:)=gvec(:,ig)
1803    else
1804      gpq(:)=gvec(:,ig)+qpt
1805    end if
1806    do igp=1,npwc
1807      if (qiszero) then
1808        gppq(:)=gvec(:,igp)
1809      else
1810        gppq(:)=gvec(:,igp)+qpt
1811      end if
1812      qpg_dot_qpgp=zero
1813      do ii=1,3
1814        qpg_dot_qpgp=qpg_dot_qpgp+&
1815 &        ( gpq(1)*b1(ii) +gpq(2)*b2(ii) +gpq(3)*b3(ii))*&
1816 &        (gppq(1)*b1(ii)+gppq(2)*b2(ii)+gppq(3)*b3(ii))
1817      end do
1818      mm(ig,igp)=rhogg(ig,igp)*qpg_dot_qpgp
1819    end do !igp
1820  end do !ig
1821 
1822  ABI_FREE(rhog_dp)
1823  ABI_FREE(rhog)
1824  ! === Now we have rhogg,rho0 ===
1825  !
1826  ! Calculate the dielectric matrix eigenvalues and vectors
1827  ! Use only the static epsm1 i.e., only the w=0 part (eps(:,:,1,:))
1828  ABI_MALLOC(eigval,(npwc))
1829 
1830  ABI_STAT_MALLOC(eigvec,(npwc,npwc), ierr)
1831  ABI_CHECK(ierr==0, 'eigvec out of memory')
1832 
1833  ABI_MALLOC(zz,(npwc))
1834  zz=czero
1835 
1836  ABI_MALLOC(qplusg,(npwc))
1837 
1838  ! Store the susceptibility matrix in upper mode before calling zhpev.
1839  ABI_STAT_MALLOC(matr,(npwc*(npwc+1)/2), ierr)
1840  ABI_CHECK(ierr==0, 'matr of memory')
1841 
1842  idx=1
1843  do ii=1,npwc
1844    do jj=1,ii
1845      matr(idx)=epsm1(jj,ii); idx=idx+1
1846    end do
1847  end do
1848 
1849  ABI_MALLOC(zhpev2,(3*npwc-2))
1850  ABI_MALLOC(zhpev1,(2*npwc-1))
1851 
1852  call ZHPEV('V','U',npwc,matr,eigval,eigvec,npwc,zhpev1,zhpev2,ierr)
1853  ABI_FREE(matr)
1854  ABI_FREE(zhpev2)
1855  ABI_FREE(zhpev1)
1856 
1857  if (ierr<0) then
1858    write (msg,'(2a,i4,a)')&
1859 &    ' Failed to calculate the eigenvalues and eigenvectors of the dielectric matrix ',ch10,&
1860 &    ierr*(-1),'-th argument in the matrix has an illegal value. '
1861    MSG_ERROR(msg)
1862  end if
1863 
1864  if (ierr>0) then
1865    write(msg,'(3a,i4,2a)')&
1866 &    ' Failed to calculate the eigenvalues and eigenvectors of the dielectric matrix ',ch10,&
1867 &    ' the algorithm failed to converge; ierr = ', ierr,ch10,&
1868 &    ' off-diagonal elements of an intermediate tridiagonal form did not converge to zero. '
1869    MSG_ERROR(msg)
1870  end if
1871  !
1872  ! Calculate the PPM parameters and the eigenpotentials needed for
1873  ! the calculation of the generalized overlap matrix
1874  ! Note: the eigenpotentials has to be calculated on the FFT (G-Gp) index
1875  !
1876  ! Save eigenvectors of \tilde\epsilon^{-1}
1877  ! MG well it is better to save \Theta otherwise
1878  ! we have to calculare \Theta for each band, spin, k-point but oh well
1879  eigtot=eigvec
1880 
1881  qlist(:,1) = qpt
1882  call cmod_qpg(1,1,qlist,npwc,gvec,gprimd,qplusg) !MG TODO here take care of small q
1883  !
1884  ! Basic Equation:
1885  !
1886  ! \Theta_{q,ii}(G)=\Psi_{q,ii}(G)/|q+G|
1887  ! where \Psi_{q,ii}(G) is the eigenvector of \tilde\epsilon^{-1}
1888 
1889  ! \tilde\omega_{ii,q}^2= 4\pi (1-eigenval(ii,q)))
1890  ! \sum_{G,Gp} \Theta^*_{q,ii}(G) (q+G)\cdot(q+Gp) n(G-Gp) \Theta_{q,ii}(Gp)
1891 
1892  do ii=1,npwc !DM band
1893    ! Calculate \Theta_{q,ii}(G)
1894    ! why the first element is not modified? if the problem is the small value of qplusg(1)
1895    ! we could multiply by sqrt(mod((q+G)(q+G'))) and then add the sing at the end
1896    if (qiszero)then
1897      eigvec(2:,ii)=eigvec(2:,ii)/qplusg(2:)
1898    else
1899      eigvec(:,ii)=eigvec(:,ii)/qplusg(:)
1900    end if
1901    do ig=1,npwc
1902      conjg_eig=CONJG(eigvec(ig,ii))
1903      do igp=1,npwc
1904        if(qiszero .and. ig==1 .and. igp==1)then
1905          zz(ii)=zz(ii)+conjg_eig*rhogg(ig,igp)*eigvec(igp,ii)
1906        else
1907          zz(ii)=zz(ii)+conjg_eig*mm(ig,igp)*eigvec(igp,ii)
1908        end if
1909      end do
1910    end do
1911 
1912    num=one-eigval(ii)
1913    if (num<=zero) then
1914 !    here I think we should set bigomegatwsq=0 and omegatw to an arbitrary value
1915 !    maybe we can output a warning TO BE discussed with Riad
1916      if (ABS(num)<1.0d-4) then
1917        num=1.0d-5
1918      else
1919        MSG_ERROR("One or more imaginary plasmon pole energies")
1920      end if
1921    end if
1922 
1923    omegatw(ii)=SQRT(4*pi*REAL(zz(ii))/num)
1924    ! this should be \alpha = 2\pi omegatw * (1-eigenval)
1925    ! MG check this, in the review I found a factor 2\pi, maybe it is reintroduced later
1926    bigomegatwsq(ii,1)=num*omegatw(ii)
1927  end do
1928 
1929  ABI_FREE(rhogg)
1930  ABI_FREE(mm)
1931  ABI_FREE(eigval)
1932  ABI_FREE(zz)
1933  ABI_FREE(eigvec)
1934  ABI_FREE(qplusg)
1935 
1936  call destroy_mpi_enreg(MPI_enreg_seq)
1937 
1938  write(msg,'(2a,f12.8,2a,3i5)')ch10,&
1939 &  ' cppm3par : omega twiddle minval [eV]  = ',MINVAL(ABS(omegatw))*Ha_eV,ch10,&
1940 &  '            omega twiddle min location = ',MINLOC(ABS(omegatw))
1941  call wrtout(std_out,msg,'COLL')
1942 
1943 end subroutine cppm3par

m_ppmodel/cppm4par [ Functions ]

[ Top ] [ m_ppmodel ] [ Functions ]

NAME

 cppm4par

FUNCTION

 Calculate the plasmon-pole parameters using Engel and Farid model (PRB47,15931,1993) [[cite:Engel1993]].
 See also Quasiparticle Calculations in Solids [[cite:Aulbur2001]] p. 23

INPUTS

  qpt(3)=Reduced coordinates of the q-point.
  epsm1(npwc,npwc)=symmetrized inverse dielectric matrix.
  gprimd(3,3)=dimensional primitive translations for reciprocal space ($\textrm{bohr}^{-1}$)
  ngfftf(18)=contain all needed information about 3D fine FFT, see ~abinit/doc/variables/vargs.htm#ngfft
  npwc=number of plane waves in epsm1
  rhor(nfftf)=charge density on the real space FFT grid
  gvec(3,npwc)=G vectors in reduced coordinated

OUTPUT

  bigomegatwsq(npwc,npwc)=plasmon-pole strength
  omegatw(npwc)=plasmon-pole frequencies

PARENTS

      m_ppmodel

CHILDREN

SOURCE

1976 subroutine cppm4par(qpt,npwc,epsm1,ngfftf,gvec,gprimd,rhor,nfftf,bigomegatwsq,omegatw)
1977 
1978 
1979 !This section has been created automatically by the script Abilint (TD).
1980 !Do not modify the following lines by hand.
1981 #undef ABI_FUNC
1982 #define ABI_FUNC 'cppm4par'
1983 !End of the abilint section
1984 
1985  implicit none
1986 
1987 !Arguments ------------------------------------
1988 !scalars
1989  integer,intent(in) :: nfftf,npwc
1990 !arrays
1991  integer,intent(in) :: gvec(3,npwc),ngfftf(18)
1992  real(dp),intent(in) :: gprimd(3,3),qpt(3)
1993  real(dp),intent(in) :: rhor(nfftf)
1994  complex(gwpc),intent(in) :: epsm1(npwc,npwc)
1995  complex(gwpc),intent(out) :: bigomegatwsq(npwc,npwc),omegatw(npwc)
1996 
1997 !Local variables-------------------------------
1998 !scalars
1999  integer,parameter :: paral_kgb0=0
2000  integer :: ierr,ig,igp,ii,ngfft1,ngfft2,ngfft3,gmgp_idx
2001  real(dp) :: qpg_dot_qpgp
2002  character(len=500) :: msg
2003  character(len=80) :: bar
2004  type(MPI_type) :: MPI_enreg_seq
2005 !arrays
2006  real(dp) :: b1(3),b2(3),b3(3),gppq(3),gpq(3),qlist(3,1)
2007  real(dp),allocatable :: eigval(:),qplusg(:),rhog_dp(:,:)
2008  real(dp),allocatable :: tmp_rhor(:)
2009  complex(dpc),allocatable :: chi(:,:),chitmp(:,:),chitmps(:,:),eigvec(:,:)
2010  complex(dpc),allocatable :: mm(:,:),mtemp(:,:),rhog(:)
2011  complex(dpc),allocatable :: rhogg(:,:),tmp1(:),zz2(:,:)
2012 
2013 !*************************************************************************
2014 
2015  DBG_ENTER("COLL")
2016 
2017  call initmpi_seq(MPI_enreg_seq)
2018  call init_distribfft_seq(MPI_enreg_seq%distribfft,'c',ngfftf(2),ngfftf(3),'all')
2019 
2020  b1=two_pi*gprimd(:,1)
2021  b2=two_pi*gprimd(:,2)
2022  b3=two_pi*gprimd(:,3)
2023  !
2024  ! === Calculate density in G space rhog(G) ===
2025  ABI_MALLOC(rhog_dp,(2,nfftf))
2026  ABI_MALLOC(rhog,(nfftf))
2027 
2028  ABI_STAT_MALLOC(rhogg,(npwc,npwc), ierr)
2029  ABI_CHECK(ierr==0, 'rhogg out of memory')
2030  !
2031  ! Conduct FFT tho(r)-->rhog(G)
2032  ! FIXME this has to be fixed, rho(G) should be passed instead of doing FFT for each q
2033  ABI_MALLOC(tmp_rhor,(nfftf))
2034  tmp_rhor=rhor ! To avoid having to use intent(inout).
2035  call fourdp(1,rhog_dp,tmp_rhor,-1,MPI_enreg_seq,nfftf,ngfftf,paral_kgb0,0)
2036  ABI_FREE(tmp_rhor)
2037 
2038  rhog(1:nfftf)=CMPLX(rhog_dp(1,1:nfftf),rhog_dp(2,1:nfftf))
2039  ABI_FREE(rhog_dp)
2040  !
2041  ! Calculate the FFT index of each (G-Gp) vector and assign the value
2042  ! of the correspondent density simultaneously
2043  ngfft1=ngfftf(1)
2044  ngfft2=ngfftf(2)
2045  ngfft3=ngfftf(3)
2046 
2047  ierr=0
2048  do ig=1,npwc
2049    do igp=1,npwc
2050      gmgp_idx = g2ifft(gvec(:,ig)-gvec(:,igp),ngfftf)
2051      if (gmgp_idx/=0) then
2052        rhogg(ig,igp)=rhog(gmgp_idx)
2053      else
2054        ierr=ierr+1
2055        rhogg(ig,igp)=czero
2056      end if
2057    end do
2058  end do
2059 
2060  if (ierr/=0) then
2061    write(msg,'(a,i0,3a)')&
2062 &   'Found ',ierr,' G1-G2 vectors falling outside the FFT box. ',ch10,&
2063 &   'Enlarge the FFT mesh to get rid of this problem. '
2064    MSG_WARNING(msg)
2065  end if
2066 
2067  ABI_FREE(rhog)
2068  !
2069  ! Now we have rhogg, calculate the M matrix (q+G1).(q+G2) n(G1-G2)
2070  ABI_STAT_MALLOC(mm,(npwc,npwc), ierr)
2071  ABI_CHECK(ierr==0, 'mm out of memory')
2072 
2073  do ig=1,npwc
2074    gpq(:)=gvec(:,ig)+qpt
2075    do igp=1,npwc
2076      gppq(:)=gvec(:,igp)+qpt
2077      qpg_dot_qpgp=zero
2078      do ii=1,3
2079        qpg_dot_qpgp=qpg_dot_qpgp+&
2080 &        ( gpq(1)*b1(ii) +gpq(2)*b2(ii) +gpq(3)*b3(ii))*&
2081 &        (gppq(1)*b1(ii)+gppq(2)*b2(ii)+gppq(3)*b3(ii))
2082      end do
2083      mm(ig,igp)=rhogg(ig,igp)*qpg_dot_qpgp
2084    end do !igp
2085  end do !ig
2086 
2087  !MG TODO too much memory in chi, we can do all this stuff inside a loop
2088  ABI_STAT_MALLOC(chitmp,(npwc,npwc), ierr)
2089  ABI_CHECK(ierr==0, 'chitmp out of memory')
2090 
2091  ABI_STAT_MALLOC(chi,(npwc,npwc), ierr)
2092  ABI_CHECK(ierr==0, 'chi out of memory')
2093 
2094  ABI_MALLOC(qplusg,(npwc))
2095  !
2096  ! Extract the full polarizability from \tilde \epsilon^{-1}
2097  ! \tilde\epsilon^{-1}_{G1 G2} = \delta_{G1 G2} + 4\pi \frac{\chi_{G1 G2}}{|q+G1| |q+G2|}
2098  chitmp(:,:)=epsm1(:,:)
2099 
2100  qlist(:,1) = qpt
2101  call cmod_qpg(1,1,qlist,npwc,gvec,gprimd,qplusg) !MG TODO here take care of small q
2102 
2103  do ig=1,npwc
2104    chitmp(ig,ig)=chitmp(ig,ig)-one
2105  end do
2106  do ig=1,npwc
2107    do igp=1,npwc
2108      chi(ig,igp)=chitmp(ig,igp)*qplusg(ig)*qplusg(igp)/four_pi
2109    end do
2110  end do
2111  ABI_FREE(chitmp)
2112  !
2113  ! === Solve chi*X = Lambda M*X where Lambda=-1/em(q)**2 ===
2114  ABI_MALLOC(eigval,(npwc))
2115 
2116  ABI_STAT_MALLOC(eigvec,(npwc,npwc), ierr)
2117  ABI_CHECK(ierr==0, "eigvec out-of-memory")
2118 
2119  ABI_STAT_MALLOC(mtemp,(npwc,npwc), ierr)
2120  ABI_CHECK(ierr==0, "mtemp out-of-memory")
2121 
2122  ABI_STAT_MALLOC(chitmps,(npwc,npwc), ierr)
2123  ABI_CHECK(ierr==0, "chitmps out-of-memory")
2124  !
2125  ! Copy chi and mm into working arrays
2126  chitmps(:,:)=chi(:,:)
2127  mtemp(:,:)=mm(:,:)
2128 
2129  call xhegv(1,"Vectors","Upper",npwc,chitmps,mtemp,eigval)
2130  !
2131  ! Eigenvectors are normalized as : X_i^* M X_j = \delta_{ij}
2132  eigvec(:,:)=chitmps(:,:)
2133  ABI_FREE(mtemp)
2134  ABI_FREE(chitmps)
2135  !
2136  ! === Calculate the plasmon pole parameters ===
2137  ABI_MALLOC(tmp1,(npwc))
2138 
2139  ABI_STAT_MALLOC(zz2,(npwc,npwc), ierr)
2140  ABI_CHECK(ierr==0, "zz2 out-of-memory")
2141  !
2142  ! good check:
2143  ! the lowest plasmon energy on gamma should be
2144  ! close to experimental plasma energy within an error of 10 %
2145  ! this error can be reduced further if one includes the non local
2146  ! commutators in the calculation of the polarizability at q==0
2147  zz2(:,:)=(0.0,0.0)
2148 
2149  qlist(:,1) = qpt
2150  call cmod_qpg(1,1,qlist,npwc,gvec,gprimd,qplusg) !MG TODO here take care of small q
2151 
2152  do ii=1,npwc
2153    !
2154    ! keeping in mind that the above matrix is negative definite
2155    ! we might have a small problem with the eigval that correspond to large G vectors
2156    ! i.e. DM band index, where the eigevalues become very small with
2157    ! possibility of being small positive numbers (due to numerical problems)
2158    ! thus as a caution one can use the following condition
2159    ! this will not affect the result since such a huge plasmon energy give almost zero
2160    ! contribution to the self correlation energy
2161 
2162    if (eigval(ii)>=zero) then
2163      eigval(ii) = -1.0d-4
2164      if (eigval(ii)>1.0d-3) then
2165        eigval(ii) = -1.0d-22
2166        write(msg,'(a,i6,a,es16.6)')&
2167 &        ' Imaginary plasmon pole eigenenergy, eigenvector number ',ii,' with eigval',eigval(ii),ch10
2168        MSG_ERROR(msg)
2169      end if
2170    end if
2171    !
2172    ! === Save plasmon energies ===
2173    omegatw(ii)=SQRT(-1/eigval(ii))
2174    !
2175    ! Calculate and save scaled plasmon-pole eigenvectors
2176    ! defined as \sqrt{4\pi} \frac{Mx}{\sqrt{\tilde\omega} |q+G|}
2177    tmp1(:)=eigvec(:,ii)
2178 
2179    do ig=1,npwc
2180      do igp=1,npwc
2181        zz2(ig,ii)=zz2(ig,ii)+mm(ig,igp)*tmp1(igp) ! z--->y
2182      end do
2183      bigomegatwsq(ig,ii)=SQRT(four_pi)*zz2(ig,ii)/SQRT(omegatw(ii))
2184      bigomegatwsq(ig,ii)=bigomegatwsq(ig,ii)/qplusg(ig)
2185    end do
2186 
2187  end do
2188 
2189  ABI_FREE(tmp1)
2190  ABI_FREE(eigvec)
2191  ABI_FREE(eigval)
2192  ABI_FREE(zz2)
2193 
2194  call destroy_mpi_enreg(MPI_enreg_seq)
2195 
2196  ABI_FREE(qplusg)
2197  ABI_FREE(chi)
2198  ABI_FREE(rhogg)
2199  ABI_FREE(mm)
2200 
2201  bar=REPEAT('-',80)
2202  write(msg,'(3a)')bar,ch10,' plasmon energies vs q vector shown for lowest 10 bands                 '
2203  call wrtout(std_out,msg,'COLL')
2204  write(msg,'(2x,5x,10f7.3)')(REAL(omegatw(ig))*Ha_eV, ig=1,10)
2205  call wrtout(std_out,msg,'COLL')
2206  write(msg,'(a)')bar
2207  call wrtout(std_out,msg,'COLL')
2208 
2209  write(msg,'(2a,f12.8,2a,3i5)')ch10,&
2210 &  ' cppm4par : omega twiddle minval [eV]  = ',MINVAL(ABS(omegatw))*Ha_eV,ch10,&
2211 &  '            omega twiddle min location = ',MINLOC(ABS(omegatw))
2212  call wrtout(std_out,msg,'COLL')
2213 
2214  DBG_EXIT("COLL")
2215 
2216 end subroutine cppm4par

m_ppmodel/cqratio [ Functions ]

[ Top ] [ m_ppmodel ] [ Functions ]

NAME

 cqratio

FUNCTION

  Calculate qratio(G,Gp,q)= (q+G)\cdot(q+Gp) / |q+G|^2 needed for Hybertsen-Louie and Plasmonpole model

INPUTS

  npwc=number of planewaves considered (used for the correlation part)
  gvec(3,npwc)=reduced coordinates of the plane waves
  q(3)=coordinates of q points
  gmet(3,3)=metric in reciprocal space
  gprimd(3,3)=reciprocal lattice vectors

OUTPUT

  qratio(npwc,npwc)=(q+G).(q+Gp)

PARENTS

      m_ppmodel,mrgscr

CHILDREN

SOURCE

2245 subroutine cqratio(npwc,gvec,q,gmet,gprimd,qratio)
2246 
2247 
2248 !This section has been created automatically by the script Abilint (TD).
2249 !Do not modify the following lines by hand.
2250 #undef ABI_FUNC
2251 #define ABI_FUNC 'cqratio'
2252 !End of the abilint section
2253 
2254  implicit none
2255 
2256 !Arguments ------------------------------------
2257 !scalars
2258  integer,intent(in) :: npwc
2259 !arrays
2260  integer,intent(in) :: gvec(3,npwc)
2261  real(dp),intent(in) :: gmet(3,3),gprimd(3,3),q(3)
2262  real(dp),intent(out) :: qratio(npwc,npwc)
2263 
2264 !Local variables ------------------------------
2265 !scalars
2266  integer :: ig,igp,ii
2267  real(dp) :: qpg_dot_qpgp
2268 !arrays
2269  real(dp) :: b1(3),b2(3),b3(3),gppq(3),gpq(3),norm(npwc)
2270 
2271 !************************************************************************
2272 
2273  b1=two_pi*gprimd(:,1); b2=two_pi*gprimd(:,2); b3=two_pi*gprimd(:,3)
2274 
2275  norm(:)=zero; qratio=zero
2276 
2277  !FIXME this loops have to be rewritten!!!!
2278  do ig=1,npwc
2279    gpq(:)=gvec(:,ig)+q
2280    norm(ig)=two_pi*SQRT(DOT_PRODUCT(gpq,MATMUL(gmet,gpq)))
2281 !  norm(ig)=normv(gpq,gmet,'g')
2282  end do
2283  do ig=1,npwc
2284    gpq(:)=gvec(:,ig)+q
2285    do igp=1,npwc
2286      gppq(:)=gvec(:,igp)+q
2287      qpg_dot_qpgp=zero
2288 !    qpg_dot_qpgp=vdotw(gpq,gppq,gmet,'g')
2289      do ii=1,3
2290        qpg_dot_qpgp=qpg_dot_qpgp+&
2291 &        ( gpq(1)*b1(ii) +  gpq(2)*b2(ii) + gpq(3)*b3(ii))*&
2292 &        (gppq(1)*b1(ii) + gppq(2)*b2(ii) +gppq(3)*b3(ii))
2293      end do
2294 
2295 !    Now calculate qratio = (q+G).(q+Gp)/|q+G|^2
2296 !    when |q+G|^2 and (q+G).(q+Gp) are both zero set (q+G).(q+Gp)/|q+G|^2 = 1
2297 !    when |q+G|^2 is zero and |q+Gp| is not zero set (q+G).(q+Gp)/|q+G|^2 = 0
2298      if (norm(ig)<0.001) then
2299        if (norm(igp)<0.001) then     ! Case q=0 and G=Gp=0
2300          qratio(ig,igp)=one
2301        else                          ! Case q=0 and G=0 and Gp !=0
2302          qratio(ig,igp)=zero
2303        end if
2304      else if (norm(igp)<0.001) then  ! Case q=0 and G= !0 and Gp=0
2305        qratio(ig,igp)=zero
2306      else
2307        qratio(ig,igp)=qpg_dot_qpgp/norm(ig)**2
2308      end if
2309 
2310    end do
2311  end do
2312 
2313 end subroutine cqratio

m_ppmodel/get_ppm_eigenvalues [ Functions ]

[ Top ] [ m_ppmodel ] [ Functions ]

NAME

  get_ppm_eigenvalues

FUNCTION

  Constructs the inverse dielectri matrix starting from the plasmon-pole
  parameters and calculates the frequency-dependent eigenvalues for each
  of the nomega frequencies specifies in the array omega.

INPUTS

OUTPUT

SIDE EFFECTS

PARENTS

      mrgscr

CHILDREN

SOURCE

1186 subroutine get_ppm_eigenvalues(PPm,iqibz,zcut,nomega,omega,Vcp,eigenvalues)
1187 
1188 
1189 !This section has been created automatically by the script Abilint (TD).
1190 !Do not modify the following lines by hand.
1191 #undef ABI_FUNC
1192 #define ABI_FUNC 'get_ppm_eigenvalues'
1193 !End of the abilint section
1194 
1195  implicit none
1196 
1197 !Arguments ------------------------------------
1198 !scalars
1199  integer,intent(in) :: iqibz,nomega
1200  type(ppmodel_t),intent(in) :: PPm
1201  type(vcoul_t),intent(in) :: Vcp
1202  real(dp),intent(in) :: zcut
1203 !arrays
1204  complex(dpc),intent(in) :: omega(nomega)
1205  complex(dpc),intent(out) :: eigenvalues(PPm%npwc,nomega)
1206 
1207 !Local variables-------------------------------
1208 !scalars
1209  integer :: info,lwork,negw,ig1,ig2,idx,sdim,iomega,ierr
1210  character(len=500) :: msg
1211 !arrays
1212  real(dp),allocatable :: ww(:),rwork(:)
1213  complex(dpc),allocatable :: work(:),Adpp(:),eigvec(:,:),wwc(:),vs(:,:),Afull(:,:)
1214  complex(dpc),allocatable :: em1q(:,:,:)
1215  logical,allocatable :: bwork(:)
1216  logical :: sortcplx !BUG in abilint
1217 
1218 ! *************************************************************************
1219 
1220  ABI_CHECK(PPm%mqmem/=0,'mqmem==0 not implemented')
1221 
1222  ABI_MALLOC(em1q,(PPm%npwc,PPm%npwc,nomega))
1223 
1224  call getem1_from_PPm(PPm,PPm%npwc,iqibz,zcut,nomega,omega,Vcp,em1q)
1225 
1226  do iomega=1,nomega
1227    !
1228    if (ABS(REAL(omega(iomega)))>0.00001) then
1229      !if (.TRUE.) then
1230      ! === Eigenvalues for a generic complex matrix ===
1231 
1232      lwork=4*2*PPm%npwc
1233      ABI_MALLOC(wwc,(PPm%npwc))
1234      ABI_MALLOC(work,(lwork))
1235      ABI_MALLOC(rwork,(PPm%npwc))
1236      ABI_MALLOC(bwork,(PPm%npwc))
1237      ABI_MALLOC(vs,(PPm%npwc,PPm%npwc))
1238      ABI_MALLOC(Afull,(PPm%npwc,PPm%npwc))
1239      Afull=em1q(:,:,iomega)
1240 
1241      !for the moment no sort, maybe here I should sort using the real part?
1242      call ZGEES('V','N',sortcplx,PPm%npwc,Afull,PPm%npwc,sdim,wwc,vs,PPm%npwc,work,lwork,rwork,bwork,info)
1243      if (info/=0) then
1244       write(msg,'(2a,i10)')' get_ppm_eigenvalues : Error in ZGEES, diagonalizing complex matrix, info = ',info
1245       call wrtout(std_out,msg,'COLL')
1246      end if
1247 
1248      eigenvalues(:,iomega)=wwc(:)
1249 
1250      ABI_FREE(wwc)
1251      ABI_FREE(work)
1252      ABI_FREE(rwork)
1253      ABI_FREE(bwork)
1254      ABI_FREE(vs)
1255      ABI_FREE(Afull)
1256 
1257    else
1258      ! === Hermitian Case ===
1259      lwork=2*PPm%npwc-1
1260      ABI_MALLOC(ww,(PPm%npwc))
1261      ABI_MALLOC(work,(lwork))
1262      ABI_MALLOC(rwork,(3*PPm%npwc-2))
1263      ABI_MALLOC(eigvec,(PPm%npwc,PPm%npwc))
1264 
1265      ABI_STAT_MALLOC(Adpp,(PPm%npwc*(PPm%npwc+1)/2), ierr)
1266      ABI_CHECK(ierr==0, "out of memory in Adpp")
1267 
1268      write(std_out,*) 'in hermitian'
1269 
1270      idx=0
1271      do ig2=1,PPm%npwc
1272        do ig1=1,ig2
1273          idx=idx+1
1274          Adpp(idx)=em1q(ig1,ig2,iomega)
1275        end do
1276      end do
1277 
1278      ! For the moment we require also the eigenvectors.
1279      call ZHPEV('V','U',PPm%npwc,Adpp,ww,eigvec,PPm%npwc,work,rwork,info)
1280 
1281      if (info/=0) then
1282        write(msg,'(2a,i10)')' get_ppm_eigenvalues : Error diagonalizing matrix, info = ',info
1283        call wrtout(std_out,msg,'COLL')
1284      end if
1285      negw = (COUNT((REAL(ww)<tol6)))
1286      if (negw/=0) then
1287        write(msg,'(a,i0,a,i0,a,f8.4)')&
1288 &       'Found negative eigenvalues. No. ',negw,' at iqibz= ',iqibz,' minval= ',MINVAL(REAL(ww))
1289         MSG_WARNING(msg)
1290      end if
1291 
1292      eigenvalues(:,iomega)=ww(:)
1293 
1294      ABI_FREE(ww)
1295      ABI_FREE(work)
1296      ABI_FREE(rwork)
1297      ABI_FREE(eigvec)
1298      ABI_FREE(Adpp)
1299    end if
1300    !
1301  end do !iomega
1302 
1303  ABI_FREE(em1q)
1304 
1305 end subroutine get_ppm_eigenvalues

m_ppmodel/getem1_from_ppm [ Functions ]

[ Top ] [ m_ppmodel ] [ Functions ]

NAME

  getem1_from_ppm

FUNCTION

  Calculate the symmetrized inverse dielectric matrix from the parameters of the plasmon-pole model.

INPUTS

OUTPUT

SIDE EFFECTS

PARENTS

      m_ppmodel

CHILDREN

SOURCE

 910 subroutine getem1_from_ppm(PPm,mpwc,iqibz,zcut,nomega,omega,Vcp,em1q,only_ig1,only_ig2)
 911 
 912 
 913 !This section has been created automatically by the script Abilint (TD).
 914 !Do not modify the following lines by hand.
 915 #undef ABI_FUNC
 916 #define ABI_FUNC 'getem1_from_ppm'
 917 !End of the abilint section
 918 
 919  implicit none
 920 
 921 !Arguments ------------------------------------
 922 !scalars
 923  integer,intent(in) :: mpwc,iqibz,nomega
 924  type(ppmodel_t),intent(in) :: PPm
 925  type(vcoul_t),intent(in) :: Vcp
 926  real(dp),intent(in) :: zcut
 927  integer,optional,intent(in) :: only_ig1,only_ig2
 928 !arrays
 929  complex(dpc),intent(in) :: omega(nomega)
 930  complex(dpc),intent(out) :: em1q(mpwc,mpwc,nomega)
 931 
 932 !Local variables-------------------------------
 933 !scalars
 934  integer :: ig1,ig2,io,idm,ig1_min,ig2_min,ig1_max,ig2_max
 935  real(dp) :: den
 936  complex(dpc) :: qpg1,qpg2,ug1,ug2
 937  complex(dpc) :: delta,em1ggp,otw,zzpq,yg1,yg2,bot1,bot2,chig1g2
 938  !character(len=500) :: msg
 939 
 940 ! *************************************************************************
 941 
 942  ABI_CHECK(PPm%mqmem/=0,'mqmem==0 not implemented')
 943 
 944  !TODO zcut should be an entry in PPm
 945  delta=CMPLX(zero,zcut)
 946 
 947  ! To save memory, a particular combination of
 948  ! ig1 and ig2 can be selected
 949  ig1_min = 1
 950  ig2_min = 1
 951  ig1_max = PPm%npwc
 952  ig2_max = PPm%npwc
 953  if (present(only_ig1)) then
 954    ig1_min = only_ig1
 955    ig1_max = only_ig1
 956  end if
 957  if (present(only_ig2)) then
 958    ig2_min = only_ig2
 959    ig2_max = only_ig2
 960  end if
 961 
 962  select case (PPm%model)
 963  case (PPM_GODBY_NEEDS, PPM_HYBERTSEN_LOUIE)
 964    do io=1,nomega
 965      !
 966      do ig2=ig2_min,ig2_max
 967        do ig1=ig1_min,ig1_max
 968         !den = omega(io)**2-REAL(PPm%omegatw(iqibz)%vals(ig1,ig2)**2)
 969         !if (den**2<zcut**2) den = omega(io)**2-REAL( (PPm%omegatw(iqibz)%vals(ig1,ig2)-delta)**2 )
 970         den = omega(io)**2 - REAL( (PPm%omegatw(iqibz)%vals(ig1,ig2)-delta)**2 )
 971         em1ggp = PPm%bigomegatwsq(iqibz)%vals(ig1,ig2)/den
 972         if (ig1==ig2) em1ggp=em1ggp+one
 973         em1q(ig1,ig2,io)=em1ggp
 974         !em1q(ig1,ig2,io)=em1ggp*Vcp%vc_sqrt(ig1,iqibz)*Vcp%vc_sqrt(ig2,iqibz)
 975        end do
 976      end do
 977      !
 978    end do !io
 979 
 980  case (PPM_LINDEN_HORSH)
 981    !TODO Check coefficients
 982    do io=1,nomega
 983      !
 984      do ig2=ig2_min,ig2_max
 985        do ig1=ig1_min,ig1_max
 986          !
 987          em1ggp=czero
 988          do idm=1,PPm%npwc
 989            !den=omega(io)**2-(PPm%omegatw(iqibz)%vals(ig1,ig2)-delta)**2
 990            !em1w(io)=em1w(io)+eigvec(ig1,idm,iqibz)*conjg(eigvec(ig2,idm,iqibz))*bigomegatwsq(ig1,ig2,iqibz)/den
 991            ug1 = PPm%eigpot(iqibz)%vals(ig1,idm)
 992            ug2 = PPm%eigpot(iqibz)%vals(ig2,idm)
 993            otw = PPm%bigomegatwsq(iqibz)%vals(idm,1)*PPm%omegatw(iqibz)%vals(idm,1)
 994            zzpq=PPm%bigomegatwsq(iqibz)%vals(idm,1)
 995            den=half*REAL(zzpq*otw*( one/(omega(io)-otw+delta) - one/(omega(io)+otw-delta) ))
 996            em1ggp=em1ggp+ug1*CONJG(ug2)*den
 997            !eigenvalues(idm,io)=one + half*REAL(zzpq*otw*( one/(omega(io)-otw+delta) - one/(omega(io)+otw-delta) ))
 998          end do
 999          if (ig2==ig1) em1ggp=em1ggp+one
1000          em1q(ig1,ig2,io)=em1ggp
1001        end do !ig1
1002      end do !ig2
1003      !
1004    end do !iomega
1005 
1006  case (PPM_ENGEL_FARID)
1007    ! Make e^-1
1008    do io=1,nomega
1009      !
1010      do ig2=ig2_min,ig2_max
1011        qpg2=one/Vcp%vc_sqrt(ig2,iqibz)
1012        do ig1=ig1_min,ig1_max
1013          qpg1=one/Vcp%vc_sqrt(ig1,iqibz)
1014 
1015          chig1g2=czero
1016          do idm=1,PPm%npwc
1017            otw =PPm%omegatw(iqibz)%vals(idm,1)
1018            bot1=PPm%bigomegatwsq(iqibz)%vals(ig1,idm)
1019            bot2=PPm%bigomegatwsq(iqibz)%vals(ig2,idm)
1020            yg1=SQRT(otw/four_pi)*qpg1*bot1
1021            yg2=SQRT(otw/four_pi)*qpg2*bot2
1022            chig1g2=chig1g2 + yg1*CONJG(yg2)/(omega(io)**2-(otw-delta)**2)
1023          end do
1024 
1025          em1ggp=four_pi*chig1g2/(qpg1*qpg2)
1026          if (ig1==ig2) em1ggp=em1ggp+one
1027          em1q(ig1,ig2,io)=em1ggp !*Vcp%vc_sqrt(ig1,iqibz)*Vcp%vc_sqrt(ig2,iqibz)
1028        end do !ig1
1029      end do !ig2
1030      !
1031    end do !iomega
1032 
1033  case default
1034    MSG_BUG(sjoin('Wrong PPm%model:',itoa(PPm%model)))
1035  end select
1036 
1037 end subroutine getem1_from_ppm

m_ppmodel/getem1_from_ppm_one_ggp [ Functions ]

[ Top ] [ m_ppmodel ] [ Functions ]

NAME

  getem1_from_ppm

FUNCTION

  Same as getem1_from_ppm, but does it for a single set of G,G' vectors

INPUTS

OUTPUT

SIDE EFFECTS

PARENTS

      m_ppmodel,mrgscr

CHILDREN

SOURCE

1062 subroutine getem1_from_ppm_one_ggp(PPm,iqibz,zcut,nomega,omega,Vcp,em1q,ig1,ig2)
1063 
1064 
1065 !This section has been created automatically by the script Abilint (TD).
1066 !Do not modify the following lines by hand.
1067 #undef ABI_FUNC
1068 #define ABI_FUNC 'getem1_from_ppm_one_ggp'
1069 !End of the abilint section
1070 
1071  implicit none
1072 
1073 !Arguments ------------------------------------
1074 !scalars
1075  integer,intent(in) :: iqibz,nomega
1076  type(ppmodel_t),intent(in) :: PPm
1077  type(vcoul_t),intent(in) :: Vcp
1078  real(dp),intent(in) :: zcut
1079  integer, intent(in) :: ig1,ig2
1080 !arrays
1081  complex(dpc),intent(in) :: omega(nomega)
1082  complex(dpc),intent(out) :: em1q(nomega)
1083 
1084 !Local variables-------------------------------
1085 !scalars
1086  integer :: io,idm !,ig1_min,ig2_min,ig2_max
1087  real(dp) :: den
1088  complex(dpc) :: qpg1,qpg2,ug1,ug2
1089  complex(dpc) :: delta,em1ggp,otw,zzpq,yg1,yg2,bot1,bot2,chig1g2
1090  !character(len=500) :: msg
1091 
1092 ! *************************************************************************
1093 
1094  ABI_CHECK(PPm%mqmem/=0,'mqmem==0 not implemented')
1095 
1096  !TODO zcut should be an entry in PPm
1097  delta=CMPLX(zero,zcut)
1098 
1099  select case (PPm%model)
1100 
1101  case (PPM_GODBY_NEEDS, PPM_HYBERTSEN_LOUIE)
1102    do io=1,nomega
1103      !den = omega(io)**2-REAL(PPm%omegatw(iqibz)%vals(ig1,ig2)**2)
1104      !if (den**2<zcut**2) den = omega(io)**2-REAL( (PPm%omegatw(iqibz)%value(ig1,ig2)-delta)**2 )
1105      den = omega(io)**2-REAL( (PPm%omegatw(iqibz)%vals(ig1,ig2)-delta)**2 )
1106      em1ggp = PPm%bigomegatwsq(iqibz)%vals(ig1,ig2)/den
1107      if (ig1==ig2) em1ggp=em1ggp+one
1108      em1q(io)=em1ggp
1109      !em1q(io)=em1ggp*Vcp%vc_sqrt(ig1,iqibz)*Vcp%vc_sqrt(ig2,iqibz)
1110    end do !io
1111 
1112  case (PPM_LINDEN_HORSH)
1113    !TODO Check coefficients
1114    do io=1,nomega
1115      em1ggp=czero
1116      do idm=1,PPm%npwc
1117        !den=omega(io)**2-(PPm%omegatw(iqibz)%vals(ig1,ig2)-delta)**2
1118        !em1w(io)=em1w(io)+eigvec(ig1,idm,iqibz)*conjg(eigvec(ig2,idm,iqibz))*bigomegatwsq(ig1,ig2,iqibz)/den
1119        ug1 =PPm%eigpot(iqibz)%vals(ig1,idm)
1120        ug2 =PPm%eigpot(iqibz)%vals(ig2,idm)
1121        otw =PPm%bigomegatwsq(iqibz)%vals(idm,1)*PPm%omegatw(iqibz)%vals(idm,1)
1122        zzpq=PPm%bigomegatwsq(iqibz)%vals(idm,1)
1123        den=half*REAL(zzpq*otw*( one/(omega(io)-otw+delta) - one/(omega(io)+otw-delta) ))
1124        em1ggp=em1ggp+ug1*CONJG(ug2)*den
1125        !eigenvalues(idm,io)=one + half*REAL(zzpq*otw*( one/(omega(io)-otw+delta) - one/(omega(io)+otw-delta) ))
1126      end do
1127 
1128      if (ig2==ig1) em1ggp=em1ggp+one
1129      em1q(io)=em1ggp
1130    end do !iomega
1131 
1132  case (PPM_ENGEL_FARID)
1133    ! Make e^-1
1134    do io=1,nomega
1135      !
1136      qpg2=one/Vcp%vc_sqrt(ig2,iqibz)
1137      qpg1=one/Vcp%vc_sqrt(ig1,iqibz)
1138 
1139      chig1g2=czero
1140      do idm=1,PPm%npwc
1141        otw =PPm%omegatw(iqibz)%vals(idm,1)
1142        bot1=PPm%bigomegatwsq(iqibz)%vals(ig1,idm)
1143        bot2=PPm%bigomegatwsq(iqibz)%vals(ig2,idm)
1144        yg1=SQRT(otw/four_pi)*qpg1*bot1
1145        yg2=SQRT(otw/four_pi)*qpg2*bot2
1146        chig1g2=chig1g2 + yg1*CONJG(yg2)/(omega(io)**2-(otw-delta)**2)
1147      end do
1148 
1149      em1ggp=four_pi*chig1g2/(qpg1*qpg2)
1150      if (ig1==ig2) em1ggp=em1ggp+one
1151      em1q(io)=em1ggp !*Vcp%vc_sqrt(ig1,iqibz)*Vcp%vc_sqrt(ig2,iqibz)
1152 
1153    end do !iomega
1154 
1155  case default
1156    MSG_BUG(sjoin('Wrong PPm%model:',itoa(PPm%model)))
1157  end select
1158 
1159 end subroutine getem1_from_ppm_one_ggp

m_ppmodel/new_setup_ppmodel [ Functions ]

[ Top ] [ m_ppmodel ] [ Functions ]

NAME

 new_setup_ppmodel

FUNCTION

  Initialize some values of several arrays of the PPm datastructure
  that are used in case of plasmonpole calculations
  Just a wrapper around different plasmonpole routines.

INPUTS

  Cryst<crystal_t>=Info on the unit cell and crystal symmetries.
  Qmesh<kmesh_t>=the q-mesh used for the inverse dielectric matrix
  iq_ibz=Index of the q-point in the BZ.
  npwe=number of G vectors for the correlation part
  nomega=number of frequencies in $\epsilon^{-1}$
  omega=frequencies in epsm1_ggw
  epsm1_ggw(npwe,npwe,nomega)=the inverse dielctric matrix
  ngfftf(18)=contain all needed information about the 3D fine FFT mesh, see ~abinit/doc/variables/vargs.htm#ngfft
  gprimd(3,3)=dimensional primitive translations for reciprocal space ($\textrm{bohr}^{-1}$)
  nfftf=the number of points in the FFT mesh (for this processor)
  rhor_tot(nfftf)=the total charge in real space
  PPm<ppmodel_t>:

SIDE EFFECTS

  PPm<ppmodel_t>:
  == if ppmodel 1 or 2 ==
   %omegatw and %bigomegatwsq
  == if ppmodel 3 ==
   %omegatw, %bigomegatwsq and %eigpot
  == if ppmodel 4 ==
   %omegatw and %bigomegatwsq

NOTES

 * FFT parallelism not implemented.
 * TODO: rhor_tot should be replaced by rhog_tot

PARENTS

      m_ppmodel,m_screen

CHILDREN

SOURCE

2719 subroutine new_setup_ppmodel(PPm,iq_ibz,Cryst,Qmesh,npwe,nomega,omega,epsm1_ggw,nfftf,gvec,ngfftf,rhor_tot)
2720 
2721 
2722 !This section has been created automatically by the script Abilint (TD).
2723 !Do not modify the following lines by hand.
2724 #undef ABI_FUNC
2725 #define ABI_FUNC 'new_setup_ppmodel'
2726 !End of the abilint section
2727 
2728  implicit none
2729 
2730 !Arguments ------------------------------------
2731 !scalars
2732  integer,intent(in) :: nfftf,npwe,nomega,iq_ibz
2733  type(kmesh_t),intent(in) :: Qmesh
2734  type(crystal_t),intent(in) :: Cryst
2735  type(ppmodel_t),intent(inout) :: PPm
2736 !arrays
2737  integer,intent(in) :: gvec(3,npwe),ngfftf(18)
2738  real(dp),intent(in) :: rhor_tot(nfftf)
2739  complex(dpc),intent(in) :: omega(nomega)
2740  complex(gwpc),intent(in) :: epsm1_ggw(npwe,npwe,nomega)
2741 
2742 !Local variables-------------------------------
2743 !scalars
2744  real(dp) :: n_at_G_zero
2745  character(len=500) :: msg
2746 !scalars
2747  real(dp) :: qpt(3)
2748 ! *************************************************************************
2749 
2750  DBG_ENTER("COLL")
2751 
2752  !@ppmodel_t
2753  !
2754  if (PPm%has_q(iq_ibz) /= PPM_TAB_ALLOCATED) then
2755    MSG_ERROR("ppmodel tables are not allocated")
2756  end if
2757 
2758  qpt = Qmesh%ibz(:,iq_ibz)
2759  PPm%has_q(iq_ibz) = PPM_TAB_STORED
2760  !
2761  ! Calculate plasmonpole parameters
2762  ! TODO ppmodel==1 by default, should be set to 0 if AC and CD
2763  SELECT CASE (PPm%model)
2764 
2765  CASE (PPM_NONE)
2766    MSG_COMMENT('Skipping Plasmonpole model calculation')
2767 
2768  CASE (PPM_GODBY_NEEDS)
2769    ! Note: the q-dependency enters only through epsilon^-1.
2770    call cppm1par(npwe,nomega,omega,PPm%drude_plsmf,epsm1_ggw,&
2771 &     PPm%omegatw(iq_ibz)%vals,PPm%bigomegatwsq(iq_ibz)%vals)
2772 
2773  CASE (PPM_HYBERTSEN_LOUIE)
2774    call cppm2par(qpt,npwe,epsm1_ggw(:,:,1),ngfftf,gvec,Cryst%gprimd,rhor_tot,nfftf,Cryst%gmet,&
2775 &    PPm%bigomegatwsq(iq_ibz)%vals,PPm%omegatw(iq_ibz)%vals,PPm%invalid_freq)
2776    !
2777    ! Quick-and-dirty change of the plasma frequency. Never executed in standard runs.
2778    if (PPm%force_plsmf>tol6) then ! Integrate the real-space density
2779       n_at_G_zero = SUM(rhor_tot(:))/nfftf
2780       ! Change the prefactor
2781       write(msg,'(2(a,es16.8))') 'Forced ppmfreq:',PPm%force_plsmf*Ha_eV,' nelec/ucvol:',n_at_G_zero
2782       MSG_WARNING(msg)
2783       PPm%force_plsmf = (PPm%force_plsmf**2)/(four_pi*n_at_G_zero)
2784       PPm%bigomegatwsq(iq_ibz)%vals = PPm%force_plsmf * PPm%bigomegatwsq(iq_ibz)%vals
2785       PPm%omegatw(iq_ibz)%vals      = PPm%force_plsmf * PPm%omegatw(iq_ibz)%vals
2786       write(msg,'(a,es16.8)') 'Plasma frequency forced in HL ppmodel, new prefactor is:',PPm%force_plsmf
2787       MSG_WARNING(msg)
2788    end if
2789 
2790  CASE (PPM_LINDEN_HORSH)
2791    ! TODO Check better double precision, this routine is in a messy state
2792    call cppm3par(qpt,npwe,epsm1_ggw(:,:,1),ngfftf,gvec,Cryst%gprimd,rhor_tot,nfftf,&
2793 &                PPm%bigomegatwsq(iq_ibz)%vals,PPm%omegatw(iq_ibz)%vals(:,1),PPm%eigpot(iq_ibz)%vals)
2794 
2795  CASE (PPM_ENGEL_FARID)  ! TODO Check better double precision, this routine is in a messy state
2796    if ((ALL(ABS(qpt)<1.0e-3))) qpt = GW_Q0_DEFAULT ! FIXME
2797 
2798    call cppm4par(qpt,npwe,epsm1_ggw(:,:,1),ngfftf,gvec,Cryst%gprimd,rhor_tot,nfftf,&
2799 &    PPm%bigomegatwsq(iq_ibz)%vals,PPm%omegatw(iq_ibz)%vals(:,1))
2800 
2801  CASE DEFAULT
2802    MSG_BUG(sjoin('Wrong PPm%model:',itoa(PPm%model)))
2803  END SELECT
2804 
2805  DBG_EXIT("COLL")
2806 
2807 end subroutine new_setup_ppmodel

m_ppmodel/ppm_free [ Functions ]

[ Top ] [ m_ppmodel ] [ Functions ]

NAME

  ppm_free

FUNCTION

  Deallocate all associated pointers defined in a variable of type ppmodel_t.

SIDE EFFECTS

  PPm<ppmodel_t>=All dynamic memory is released.

PARENTS

      m_screen,mrgscr,sigma

CHILDREN

SOURCE

397 subroutine ppm_free(PPm)
398 
399 
400 !This section has been created automatically by the script Abilint (TD).
401 !Do not modify the following lines by hand.
402 #undef ABI_FUNC
403 #define ABI_FUNC 'ppm_free'
404 !End of the abilint section
405 
406  implicit none
407 
408 !Arguments ------------------------------------
409  type(ppmodel_t),intent(inout) :: PPm
410 
411 !Local variables-------------------------------
412 !scalars
413  integer :: dim_q,iq_ibz
414 
415 ! *********************************************************************
416 
417  !@ppmodel_t
418 
419  ! Be careful here to avoid dangling pointers.
420  if (associated(PPm%bigomegatwsq_qbz)) then
421    select case (PPm%bigomegatwsq_qbz_stat)
422    case (PPM_ISALLOCATED)
423      call array_free(PPm%bigomegatwsq_qbz)
424    case (PPM_ISPOINTER)
425      nullify(PPm%bigomegatwsq_qbz)
426    end select
427  end if
428 
429  if (associated(PPm%omegatw_qbz)) then
430    select case (PPm%omegatw_qbz_stat)
431    case (PPM_ISALLOCATED)
432      call array_free(PPm%omegatw_qbz)
433    case (PPM_ISPOINTER)
434      nullify(PPm%omegatw_qbz)
435    end select
436  end if
437 
438  if (associated(PPm%eigpot_qbz)) then
439    select case (PPm%eigpot_qbz_stat)
440    case (PPM_ISALLOCATED)
441      call array_free(PPm%eigpot_qbz)
442    case (PPM_ISPOINTER)
443      nullify(PPm%eigpot_qbz)
444    end select
445  end if
446 
447  dim_q=PPm%nqibz; if (PPm%mqmem==0) dim_q=1
448 
449 #if 0
450  do iq_ibz=1,dim_q
451    call ppm_table_free(PPm,iq_ibz)
452  end do
453  ABI_FREE(PPm%bigomegatwsq)
454  ABI_FREE(PPm%omegatw)
455  ABI_FREE(PPm%eigpot)
456 
457 #else
458  if (allocated(PPm%bigomegatwsq)) then
459    do iq_ibz=1,dim_q
460      call array_free(PPm%bigomegatwsq(iq_ibz))
461    end do
462    ABI_DT_FREE(PPm%bigomegatwsq)
463  end if
464  !
465  if (allocated(PPm%omegatw)) then
466    do iq_ibz=1,dim_q
467      call array_free(PPm%omegatw(iq_ibz))
468    end do
469    ABI_DT_FREE(PPm%omegatw)
470  end if
471  !
472  if (allocated(PPm%eigpot)) then
473    do iq_ibz=1,dim_q
474      call array_free(PPm%eigpot(iq_ibz))
475    end do
476    ABI_DT_FREE(PPm%eigpot)
477  end if
478 #endif
479 
480  ! logical flags must be deallocated here.
481  if (allocated(PPm%keep_q)) then
482    ABI_FREE(PPm%keep_q)
483  end if
484  if (allocated(PPm%has_q)) then
485    ABI_FREE(PPm%has_q)
486  end if
487 
488 end subroutine ppm_free

m_ppmodel/ppm_get_qbz [ Functions ]

[ Top ] [ m_ppmodel ] [ Functions ]

NAME

  ppm_get_qbz

FUNCTION

  Calculates the plasmonpole matrix elements in the full BZ zone.

INPUTS

  PPm<ppmodel_t>=data type containing information on the plasmonpole technique
  Gsph<gsphere_t>=data related to the G-sphere
    %grottb
    %phmSGt
  Qmesh<kmesh_t>=Info on the q-mesh
    %nbz=number if q-points in the BZ
    %tab(nbz)=index of the symmeric q-point in the IBZ, for each point in the BZ
    %tabo(nbz)=the operation that rotates q_ibz onto \pm q_bz (depending on tabi)
    %tabi(nbz)=-1 if time-reversal has to be considered, 1 otherwise
  iq_bz=Index of the q-point in the BZ where PPmodel parameters have to be symmetrized

OUTPUT

  botsq
  otq
  eig (only if PPm%ppmodel==3)

NOTES

  In the present implementation we are not considering a possible umklapp vector G0.
  In this case,indeed, the equation is different since we have to consider G-G0.
  There is however a check in sigma

  * Remember the symmetry properties of \tilde\espilon^{-1}
    If q_bz=Sq_ibz+G0:

    $\epsilon^{-1}_{SG1-G0,SG2-G0}(q_bz) = e^{+iS(G2-G1).\tau}\epsilon^{-1}_{G1,G2)}(q)

    If time-reversal symmetry can be used then :
    $\epsilon^{-1}_{G1,G2}(-q_bz) = e^{+i(G1-G2).\tau}\epsilon^{-1}_{-S^{-1}(G1+Go),-S^{-1}(G2+G0)}^*(q)

 * Notice that eig is used only if PPm%model==3

PARENTS

      calc_sigc_me,m_ppmodel

CHILDREN

SOURCE

228 subroutine ppm_get_qbz(PPm,Gsph,Qmesh,iq_bz,botsq,otq,eig)
229 
230 
231 !This section has been created automatically by the script Abilint (TD).
232 !Do not modify the following lines by hand.
233 #undef ABI_FUNC
234 #define ABI_FUNC 'ppm_get_qbz'
235 !End of the abilint section
236 
237  implicit none
238 
239 !Arguments ------------------------------------
240 !scalars
241  integer,intent(in) :: iq_bz
242  type(ppmodel_t),target,intent(inout) :: PPm
243  type(gsphere_t),target,intent(in) :: Gsph
244  type(kmesh_t),intent(in) :: Qmesh
245  complex(gwpc),intent(out) :: botsq(:,:),otq(:,:),eig(:,:)
246 
247 !Local variables-------------------------------
248 !scalars
249  integer :: ii,jj,iq_ibz,itim_q,isym_q,iq_curr,isg1,isg2
250  !character(len=500) :: msg
251 !arrays
252  integer, ABI_CONTIGUOUS pointer :: grottb(:)
253  complex(gwpc),ABI_CONTIGUOUS pointer :: phsgt(:),bigomegatwsq(:,:),omegatw(:,:)
254 ! *********************************************************************
255 
256  !@ppmodel_t
257  ! Save the index of the q-point for checking purpose.
258  PPm%iq_bz = iq_bz
259 
260  ! Here there is a problem with the small q, still cannot use BZ methods
261  iq_ibz=Qmesh%tab(iq_bz)
262  isym_q=Qmesh%tabo(iq_bz)
263  itim_q=(3-Qmesh%tabi(iq_bz))/2
264 
265  !call get_bz_item(Qmesh,iq_bz,qbz,iq_ibz,isym_q,itim_q,isirred=q_isirred)
266 
267  iq_curr=iq_ibz; if (PPm%mqmem==0) iq_curr=1
268 
269  grottb => Gsph%rottb (1:PPm%npwc,itim_q,isym_q)
270  phsgt  => Gsph%phmSGt(1:PPm%npwc,isym_q)
271 
272  bigomegatwsq => PPm%bigomegatwsq(iq_curr)%vals
273  omegatw      => PPm%omegatw(iq_curr)%vals
274 
275  ! Symmetrize the PPM parameters
276  SELECT CASE (PPm%model)
277  CASE (PPM_GODBY_NEEDS, PPM_HYBERTSEN_LOUIE)
278    ! Plasmon pole frequencies otq are invariant under symmetry
279 !$omp parallel do private(isg1, isg2)
280    do jj=1,PPm%npwc
281      isg2 = grottb(jj)
282      do ii=1,PPm%npwc
283        isg1 = grottb(ii)
284        botsq(isg1,isg2) = bigomegatwsq(ii,jj)*phsgt(ii)*CONJG(phsgt(jj))
285        otq  (isg1,isg2) = omegatw(ii,jj)
286      end do
287    end do
288 
289  CASE (PPM_LINDEN_HORSH)
290    ! * For notations see pag 22 of Quasiparticle Calculations in solid (Aulbur et al)
291    !  If q_bz=Sq_ibz+G0 then:
292    !
293    ! $\omega^2_{ii}(q_bz) = \omega^2_{ii}(q)$        (otq array)
294    ! $\alpha_{ii}(q_bz)   = \alpha_{ii}(q)$          (botq array
295    ! $\Phi_{SG-G0}(q_bz)  = \Phi_{G}(q) e^{-iSG.t}$  (eigenvectors of e^{-1}, eig array)
296    !
297    do ii=1,PPm%npwc ! DM bands index
298      botsq(ii,1) = bigomegatwsq(ii,1)
299      otq  (ii,1) = omegatw     (ii,1)
300      do jj=1,PPm%npwc
301        eig(grottb(jj),ii) = PPm%eigpot(iq_curr)%vals(jj,ii) * phsgt(jj)
302      end do
303    end do
304    if (itim_q==2) eig=CONJG(eig) ! Time-reversal
305 
306  CASE (PPM_ENGEL_FARID)
307    ! * For notations see pag 23 of Quasiparticle Calculations in solid (Aulbur et al)
308    ! If q_bz=Sq_ibz+G0 then:
309    !
310    ! $\omega^2_{ii}(q_bz) = \omega^2_{ii}(q)$        (otq array)
311    ! $y_{SG-G0}(q_bz)     = y_{G}(q) e^{-iSG.t}$     (y=Lx)
312    !
313    do ii=1,PPm%npwc ! DM bands index
314      otq(ii,1) = omegatw(ii,1)
315      do jj=1,PPm%npwc
316        botsq(grottb(jj),ii) = bigomegatwsq(jj,ii)*phsgt(jj)
317      end do
318    end do
319 
320  CASE DEFAULT
321    MSG_BUG(sjoin('Wrong PPm%model:',itoa(PPm%model)))
322  END SELECT
323  !
324  ! * Take into account time-reversal symmetry.
325  if (itim_q==2) then
326 !$omp parallel workshare
327    botsq=CONJG(botsq)
328 !$omp end parallel workshare
329  end if
330 
331 end subroutine ppm_get_qbz

m_ppmodel/ppm_init [ Functions ]

[ Top ] [ m_ppmodel ] [ Functions ]

NAME

  ppm_init

FUNCTION

  Initialize dimensions and other useful variables related to the PPmodel

INPUTS

 ppmodel=
 drude_plsmf=

SIDE EFFECTS

  PPm<ppmodel_t>=Arrays needed to store the ppmodel parameters are allocated with
   proper dimensions according to the plasmon-pole model.

PARENTS

      m_screen,mrgscr,sigma

CHILDREN

SOURCE

623 subroutine ppm_init(PPm,mqmem,nqibz,npwe,ppmodel,drude_plsmf,invalid_freq)
624 
625 
626 !This section has been created automatically by the script Abilint (TD).
627 !Do not modify the following lines by hand.
628 #undef ABI_FUNC
629 #define ABI_FUNC 'ppm_init'
630 !End of the abilint section
631 
632  implicit none
633 
634 !Arguments ------------------------------------
635  integer,intent(in) :: mqmem,nqibz,npwe,ppmodel,invalid_freq
636  real(dp),intent(in) :: drude_plsmf
637  type(ppmodel_t),intent(out) :: PPm
638 
639 !Local variables-------------------------------
640 !scalars
641  integer :: dim_q,iq_ibz,ierr
642  logical :: ltest
643  !character(len=500) :: msg
644 ! *********************************************************************
645 
646  DBG_ENTER("COLL")
647 
648  !@ppmodel_t
649  call ppm_nullify(PPm)
650 
651  PPm%invalid_freq=invalid_freq
652  PPm%nqibz = nqibz
653  PPm%mqmem = mqmem
654  ltest = (PPm%mqmem==0.or.PPm%mqmem==PPm%nqibz)
655 ! write(std_out,'(a,I0)') ' PPm%mqmem = ',PPm%mqmem
656 ! write(std_out,'(a,I0)') ' PPm%nqibz = ',PPm%nqibz
657 ! ABI_CHECK(ltest,'Wrong value for mqmem')
658 
659  PPm%npwc        = npwe
660  PPm%model       = ppmodel
661  PPm%drude_plsmf = drude_plsmf
662  PPm%userho      = 0
663  if (ANY(ppmodel==(/PPM_HYBERTSEN_LOUIE, PPM_LINDEN_HORSH, PPM_ENGEL_FARID/))) PPm%userho=1
664 
665  ABI_MALLOC(PPm%keep_q,(nqibz))
666  PPm%keep_q=.FALSE.; if (PPm%mqmem>0) PPm%keep_q=.TRUE.
667 
668  ABI_MALLOC(PPm%has_q,(nqibz))
669  PPm%has_q=PPM_NOTAB
670  !
671  ! Full q-mesh is stored or out-of-memory solution.
672  dim_q=PPm%nqibz; if (PPm%mqmem==0) dim_q=1
673 
674  ABI_DT_MALLOC(PPm%bigomegatwsq, (dim_q))
675  ABI_DT_MALLOC(PPm%omegatw,(dim_q))
676  ABI_DT_MALLOC(PPm%eigpot,(dim_q))
677 
678  SELECT CASE (PPm%model)
679 
680  CASE (PPM_NONE)
681    MSG_WARNING("Called with ppmodel==0")
682    PPm%dm2_botsq = 0
683    PPm%dm2_otq   = 0
684    PPm%dm_eig    = 0
685    RETURN
686 
687  CASE (PPM_GODBY_NEEDS, PPM_HYBERTSEN_LOUIE)
688    PPm%dm2_botsq = PPm%npwc
689    PPm%dm2_otq   = PPm%npwc
690    PPm%dm_eig    = 1 ! Should be set to 0, but g95 doesnt like zero-sized arrays
691 
692  CASE (PPM_LINDEN_HORSH)
693    PPm%dm2_botsq = 1
694    PPm%dm2_otq   = 1
695    PPm%dm_eig    = PPm%npwc
696 
697  CASE (PPM_ENGEL_FARID)
698    PPm%dm2_botsq = PPm%npwc
699    PPm%dm2_otq   = 1
700    PPm%dm_eig    = 1 ! Should be set to 0, but g95 doesnt like zero-sized arrays
701 
702  CASE DEFAULT
703    MSG_BUG(sjoin('Wrong PPm%model:',itoa(PPm%model)))
704  END SELECT
705  !
706  ! Allocate tables depending on the value of keep_q.
707  do iq_ibz=1,dim_q
708    !%if (keep_q(iq_ibz)) then
709 #if 1
710    ABI_STAT_MALLOC(PPm%bigomegatwsq(iq_ibz)%vals, (PPm%npwc,PPm%dm2_botsq), ierr)
711    ABI_CHECK(ierr==0, "out of memory bigomegatwsq")
712    !
713    ABI_STAT_MALLOC(PPm%omegatw(iq_ibz)%vals, (PPm%npwc,PPm%dm2_otq), ierr)
714    ABI_CHECK(ierr==0, "out of memory omegatw")
715    !
716    ABI_STAT_MALLOC(PPm%eigpot(iq_ibz)%vals, (PPm%dm_eig,PPm%dm_eig), ierr)
717    ABI_CHECK(ierr==0, "out of memory eigpot")
718    !%endif
719 #else
720    call ppm_mallocq(PPm,iq_ibz)
721 #endif
722    !%end if
723  end do
724 
725  DBG_EXIT("COLL")
726 
727 end subroutine ppm_init

m_ppmodel/ppm_mallocq [ Functions ]

[ Top ] [ m_ppmodel ] [ Functions ]

NAME

  ppm_mallocq

FUNCTION

  Allocate the ppmodel tables for the selected q-point in the IBZ.

 INPUT
  iq_ibz=Index of the q-point in the IBZ.

SIDE EFFECTS

  PPm<ppmodel_t>=PPm tables are allocated.

PARENTS

      m_ppmodel

CHILDREN

SOURCE

513 subroutine ppm_mallocq(PPm,iq_ibz)
514 
515 
516 !This section has been created automatically by the script Abilint (TD).
517 !Do not modify the following lines by hand.
518 #undef ABI_FUNC
519 #define ABI_FUNC 'ppm_mallocq'
520 !End of the abilint section
521 
522  implicit none
523 
524 !Arguments ------------------------------------
525  integer,intent(in) :: iq_ibz
526  type(ppmodel_t),intent(inout) :: PPm
527 
528 ! *********************************************************************
529 
530  !@ppmodel_t
531  ABI_MALLOC(PPm%bigomegatwsq(iq_ibz)%vals, (PPm%npwc,PPm%dm2_botsq))
532 
533  ABI_MALLOC(PPm%omegatw(iq_ibz)%vals, (PPm%npwc,PPm%dm2_otq))
534 
535  ABI_MALLOC(PPm%eigpot(iq_ibz)%vals, (PPm%dm_eig,PPm%dm_eig))
536 
537  PPm%has_q(iq_ibz) = PPM_TAB_ALLOCATED
538 
539 end subroutine ppm_mallocq

m_ppmodel/ppm_nullify [ Functions ]

[ Top ] [ m_ppmodel ] [ Functions ]

NAME

  ppm_nullify

FUNCTION

  Nullify dynamic entities in a ppmodel_t object.

INPUTS

OUTPUT

PARENTS

      m_ppmodel,m_screen

CHILDREN

SOURCE

354 subroutine ppm_nullify(PPm)
355 
356 
357 !This section has been created automatically by the script Abilint (TD).
358 !Do not modify the following lines by hand.
359 #undef ABI_FUNC
360 #define ABI_FUNC 'ppm_nullify'
361 !End of the abilint section
362 
363  implicit none
364 
365 !Arguments ------------------------------------
366  type(ppmodel_t),intent(inout) :: PPm
367 ! *********************************************************************
368 
369  !@ppmodel_t
370  ! types
371  nullify(Ppm%bigomegatwsq_qbz)
372  nullify(Ppm%omegatw_qbz)
373  nullify(Ppm%eigpot_qbz)
374 
375 end subroutine ppm_nullify

m_ppmodel/ppm_symmetrizer [ Functions ]

[ Top ] [ m_ppmodel ] [ Functions ]

NAME

  ppm_symmetrizer

FUNCTION

  Symmetrize the plasmonpole matrix elements in the full BZ zone.

INPUTS

  iq_bz=Index of the q-point in the BZ where the PPmodel parameters are wanted.
  Gsph<gsphere_t>=data related to the G-sphere.
  Cryst<crystal_t>=Info on the unit cell and crystal symmetries.
  Qmesh<kmesh_t>=the q-mesh used for the inverse dielectric matrix
  iq_ibz=Index of the q-point in the BZ.
  npwe=number of G vectors for the correlation part
  nomega=number of frequencies in $\epsilon^{-1}$
  omega=frequencies in epsm1_ggw
  epsm1_ggw(npwe,npwe,nomega)=the inverse dielctric matrix
  ngfftf(18)=contain all needed information about the 3D fine FFT mesh, see ~abinit/doc/variables/vargs.htm#ngfft
  gprimd(3,3)=dimensional primitive translations for reciprocal space ($\textrm{bohr}^{-1}$)
  nfftf=the number of points in the FFT mesh (for this processor)
  rhor_tot(nfftf)=the total charge in real space

SIDE EFFECTS

  PPm<ppmodel_t>=data type containing information on the plasmonpole technique.
  Internal tables are modified so that they (point|store) the plasmon-pole parameters
  for the specified q-point in the BZ.

PARENTS

      m_screen

CHILDREN

SOURCE

2554 subroutine ppm_symmetrizer(PPm,iq_bz,Cryst,Qmesh,Gsph,npwe,nomega,omega,epsm1_ggw,nfftf,ngfftf,rhor_tot)
2555 
2556 
2557 !This section has been created automatically by the script Abilint (TD).
2558 !Do not modify the following lines by hand.
2559 #undef ABI_FUNC
2560 #define ABI_FUNC 'ppm_symmetrizer'
2561 !End of the abilint section
2562 
2563  implicit none
2564 
2565 !Arguments ------------------------------------
2566 !scalars
2567  integer,intent(in) :: nfftf,npwe,nomega,iq_bz
2568  type(crystal_t),intent(in) :: Cryst
2569  type(gsphere_t),intent(in) :: Gsph
2570  type(kmesh_t),intent(in) :: Qmesh
2571  type(ppmodel_t),target,intent(inout) :: PPm
2572 !arrays
2573  integer,intent(in) :: ngfftf(18)
2574  real(dp),intent(in) :: rhor_tot(nfftf)
2575  complex(dpc),intent(in) :: omega(nomega)
2576  complex(gwpc),intent(in) :: epsm1_ggw(npwe,npwe,nomega)
2577 
2578 !Local variables-------------------------------
2579 !scalars
2580  integer :: iq_ibz,itim_q,isym_q,iq_curr,ierr
2581  logical :: q_isirred
2582  !character(len=500) :: msg
2583 !arrays
2584  real(dp) :: qbz(3)
2585 ! *********************************************************************
2586 
2587  !@ppmodel_t
2588  !
2589  ! Save the index of the q-point in the BZ for checking purpose.
2590  PPm%iq_bz = iq_bz
2591  !
2592  ! Here there is a problem with the small q, still cannot use BZ methods
2593  !iq_ibz=Qmesh%tab(iq_bz)
2594  !isym_q=Qmesh%tabo(iq_bz)
2595  !itim_q=(3-Qmesh%tabi(iq_bz))/2
2596 
2597  call get_bz_item(Qmesh,iq_bz,qbz,iq_ibz,isym_q,itim_q,isirred=q_isirred)
2598  iq_curr=iq_ibz; if (PPm%mqmem==0) iq_curr=1
2599  !
2600  ! =======================================================
2601  ! ==== Branching for in-core or out-of-core solution ====
2602  ! =======================================================
2603  !
2604  if (PPm%has_q(iq_ibz)==PPM_NOTAB) then
2605    ! Allocate the tables for this q_ibz
2606    call ppm_mallocq(PPm,iq_ibz)
2607  end if
2608 
2609  if (PPm%has_q(iq_ibz)==PPM_TAB_ALLOCATED) then
2610    ! Calculate the ppmodel tables for this q_ibz
2611    call new_setup_ppmodel(PPm,iq_ibz,Cryst,Qmesh,npwe,nomega,omega,epsm1_ggw,nfftf,Gsph%gvec,ngfftf,rhor_tot) !Optional
2612  end if
2613 
2614  if (q_isirred) then
2615    ! Symmetrization is not needed.
2616    if (PPm%bigomegatwsq_qbz_stat==PPM_ISALLOCATED)  then
2617      ABI_FREE(PPm%bigomegatwsq_qbz%vals)
2618    end if
2619    if (PPm%omegatw_qbz_stat==PPM_ISALLOCATED)  then
2620      ABI_FREE(PPm%omegatw_qbz%vals)
2621    end if
2622    if (PPm%eigpot_qbz_stat==PPM_ISALLOCATED)  then
2623      ABI_FREE(PPm%eigpot_qbz%vals)
2624    end if
2625 
2626    ! Point the data in memory and change the status.
2627    PPm%bigomegatwsq_qbz => PPm%bigomegatwsq(iq_ibz)
2628    PPm%bigomegatwsq_qbz_stat = PPM_ISPOINTER
2629 
2630    PPm%omegatw_qbz => PPm%omegatw(iq_ibz)
2631    PPm%omegatw_qbz_stat = PPM_ISPOINTER
2632 
2633    PPm%eigpot_qbz => PPm%eigpot(iq_ibz)
2634    PPm%eigpot_qbz_stat = PPM_ISPOINTER
2635  else
2636    ! Allocate memory if not done yet.
2637    if (PPm%bigomegatwsq_qbz_stat==PPM_ISPOINTER) then
2638       nullify(PPm%bigomegatwsq_qbz)
2639       ABI_STAT_MALLOC(PPm%bigomegatwsq_qbz%vals, (PPm%npwc,PPm%dm2_botsq), ierr)
2640       ABI_CHECK(ierr==0, "out of memory bigomegatwsq")
2641       PPm%bigomegatwsq_qbz_stat=PPM_ISALLOCATED
2642    end if
2643 
2644    if (PPm%omegatw_qbz_stat==PPM_ISPOINTER) then
2645       nullify(PPm%omegatw_qbz)
2646       ABI_STAT_MALLOC(PPm%omegatw_qbz%vals,(PPm%npwc,PPm%dm2_otq), ierr)
2647       ABI_CHECK(ierr==0, "out of memory omegatw")
2648       PPm%omegatw_qbz_stat=PPM_ISALLOCATED
2649    end if
2650 
2651    if (PPm%eigpot_qbz_stat==PPM_ISPOINTER) then
2652      nullify(PPm%eigpot_qbz)
2653      ABI_STAT_MALLOC(PPm%eigpot_qbz%vals,(PPm%dm_eig,PPm%dm_eig), ierr)
2654      ABI_CHECK(ierr==0, "out of memory eigpot")
2655      PPm%eigpot_qbz_stat=PPM_ISALLOCATED
2656    end if
2657    !
2658    ! Calculate new table for this q-point in the BZ.
2659    ! Beware: Dimensions should not change.
2660    !botsq => PPm%bigomegatwsq_qbz%vals
2661    !otq   => PPm%omegatw_qbz%vals
2662    !eig   => PPm%eigpot_qbz%vals
2663 
2664    call ppm_get_qbz(PPm,Gsph,Qmesh,iq_bz,PPm%bigomegatwsq_qbz%vals,&
2665 &    PPm%omegatw_qbz%vals,PPm%eigpot_qbz%vals)
2666 
2667    ! Release the table in the IBZ if required.
2668    if (.not.PPm%keep_q(iq_ibz)) call ppm_table_free(PPm,iq_ibz)
2669  end if
2670 
2671 end subroutine ppm_symmetrizer

m_ppmodel/ppm_table_free [ Functions ]

[ Top ] [ m_ppmodel ] [ Functions ]

NAME

  ppm_table_free

FUNCTION

  Free the ppmodel tables for the selected q-point in the IBZ.

 INPUT
  iq_ibz=Index of the q-point in the IBZ.

SIDE EFFECTS

  PPm<ppmodel_t>=PPm tables are deallocated.

PARENTS

      m_ppmodel

CHILDREN

SOURCE

564 subroutine ppm_table_free(PPm,iq_ibz)
565 
566 
567 !This section has been created automatically by the script Abilint (TD).
568 !Do not modify the following lines by hand.
569 #undef ABI_FUNC
570 #define ABI_FUNC 'ppm_table_free'
571 !End of the abilint section
572 
573  implicit none
574 
575 !Arguments ------------------------------------
576  integer,intent(in) :: iq_ibz
577  type(ppmodel_t),intent(inout) :: PPm
578 
579 ! *********************************************************************
580 
581  !@ppmodel_t
582  if (allocated(PPm%bigomegatwsq)) then
583    call array_free(PPm%bigomegatwsq(iq_ibz))
584  end if
585 
586  if (allocated(PPm%omegatw)) then
587    call array_free(PPm%omegatw(iq_ibz))
588  end if
589 
590  if (allocated(PPm%eigpot)) then
591    call array_free(PPm%eigpot(iq_ibz))
592  end if
593 
594  PPm%has_q(iq_ibz) = PPM_NOTAB
595 
596 end subroutine ppm_table_free

m_ppmodel/ppm_times_ket [ Functions ]

[ Top ] [ m_ppmodel ] [ Functions ]

NAME

 ppm_times_ket

FUNCTION

  Calculate the contribution to self-energy operator using a plasmon-pole model.

INPUTS

  nomega=Number of frequencies.
  nspinor=Number of spinorial components.
  npwc=Number of G vectors in the plasmon pole.
  npwx=number of G vectors in rhotwgp, i.e. no. of G-vectors for the exchange part.
  theta_mu_minus_e0i= $\theta(\mu-\epsilon_{k-q,b1,s}), defines if the state is occupied or not
  zcut=Small imaginary part to avoid the divergence. (see related input variable)
  omegame0i(nomega)=Frequencies used to evaluate \Sigma_c ($\omega$ - $\epsilon_i)$
  PPm<ppmodel_t>=structure gathering info on the Plasmon-pole technique.
     %model=type plasmon pole model
     %dm2_botsq= 1 if model==3, =npwc if model== 4, 1 for all the other cases
     %dm2_otq= 1 if model==3, =1    if model== 4, 1 for all the other cases
     %dm_eig=npwc if model=3, 0 otherwise
     %botsq(npwc,dm2_botsq)=Plasmon pole parameters for this q-point.
     %eig(dm_eig,dm_eig)=The eigvectors of the symmetrized inverse dielectric matrix for this q point
       (first index for G, second index for bands)
     %otq(npwc,dm2_otq)=Plasmon pole parameters for this q-point.
  rhotwgp(npwx)=oscillator matrix elements divided by |q+G| i.e.
    $\frac{\langle b1 k-q s | e^{-i(q+G)r | b2 k s \rangle}{|q+G|}$

OUTPUT

  sigcme(nomega) (to be described), only relevant if ppm3 or ppm4
  ket(npwc,nomega):
   === model==1,2 ====

   ket(G,omega) = Sum_G2       conjg(rhotw(G)) * Omega(G,G2) * rhotw(G2)
                          ---------------------------------------------------
                            2 omegatw(G,G2) (omega-E_i + omegatw(G,G2)(2f-1))

PARENTS

CHILDREN

SOURCE

2855 subroutine ppm_times_ket(PPm,nspinor,npwc,nomega,rhotwgp,omegame0i,zcut,theta_mu_minus_e0i,npwx,ket,sigcme)
2856 
2857 
2858 !This section has been created automatically by the script Abilint (TD).
2859 !Do not modify the following lines by hand.
2860 #undef ABI_FUNC
2861 #define ABI_FUNC 'ppm_times_ket'
2862 !End of the abilint section
2863 
2864  implicit none
2865 
2866 !Arguments ------------------------------------
2867 !scalars
2868  integer,intent(in) :: nomega,npwc,npwx,nspinor
2869  real(dp),intent(in) :: theta_mu_minus_e0i,zcut
2870  type(ppmodel_t),intent(in) :: PPm
2871 !arrays
2872  real(dp),intent(in) :: omegame0i(nomega)
2873  complex(gwpc),intent(in) :: rhotwgp(npwx*nspinor)
2874  complex(gwpc),intent(inout) :: ket(npwc*nspinor,nomega)
2875  complex(gwpc),intent(out) :: sigcme(nomega)
2876 
2877 !Local variables-------------------------------
2878 !scalars
2879  integer :: ig,igp,ii,ios,ispinor,spadc,spadx
2880  real(dp) :: den,ff,inv_den,omegame0i_io,otw,twofm1,twofm1_zcut
2881  complex(gwpc) :: ct,num,numf,rhotwgdp_igp
2882  logical :: fully_occupied,totally_empty
2883 !arrays
2884  complex(gwpc),allocatable :: rhotwgdpcc(:)
2885  complex(gwpc), ABI_CONTIGUOUS pointer :: botsq(:,:),eig(:,:),otq(:,:)
2886 
2887 !*************************************************************************
2888 
2889  SELECT CASE (PPm%model)
2890  CASE (PPM_GODBY_NEEDS, PPM_HYBERTSEN_LOUIE)
2891    botsq  => PPm%bigomegatwsq_qbz%vals !(1:npwc,1:PPm%dm2_botsq)
2892    otq    => PPm%omegatw_qbz%vals      !(1:npwc,1:PPm%dm2_otq)
2893    fully_occupied =(ABS(theta_mu_minus_e0i-one)<0.001)
2894    totally_empty  =(ABS(theta_mu_minus_e0i    )<0.001)
2895 
2896    do ispinor=1,nspinor
2897      spadx=(ispinor-1)*npwx; spadc=(ispinor-1)*npwc
2898 
2899      if (.not.totally_empty) then
2900        ! \Bomega^2_{G1G2}/\omegat_{G1G2} M_{G1,G2}. \theta(\mu-e_s) / (\omega+\omegat_{G1G2}-e_s-i\delta)
2901        twofm1_zcut=zcut
2902 !$omp parallel do private(omegame0i_io, rhotwgdp_igp, otw, num, den)
2903        do ios=1,nomega
2904          omegame0i_io=omegame0i(ios)
2905          do igp=1,npwc
2906            rhotwgdp_igp=rhotwgp(spadx+igp)
2907            do ig=1,npwc
2908              otw=DBLE(otq(ig,igp)) !in principle otw -> otw - ieta
2909              num = botsq(ig,igp)*rhotwgdp_igp
2910              den = omegame0i_io+otw
2911              if (den**2>zcut**2) then
2912                ket(spadc+ig,ios)=ket(spadc+ig,ios) + num/(den*otw)*theta_mu_minus_e0i
2913              else
2914                ket(spadc+ig,ios)=ket(spadc+ig,ios) + &
2915 &                num*CMPLX(den,twofm1_zcut)/((den**2+twofm1_zcut**2)*otw)*theta_mu_minus_e0i
2916              end if
2917            end do !ig
2918          end do !igp
2919        end do !ios
2920      end if !not totally empty
2921 
2922      if (.not.(fully_occupied)) then
2923        ! \Bomega^2_{G1G2}/\omegat_{G1G2} M_{G1,G2}. \theta(e_s-\mu) / (\omega-\omegat_{G1G2}-e_s+i\delta)
2924        twofm1_zcut=-zcut
2925 !$omp parallel do private(omegame0i_io, rhotwgdp_igp, otw, num, den)
2926        do ios=1,nomega
2927          omegame0i_io=omegame0i(ios)
2928          do igp=1,npwc
2929            rhotwgdp_igp=rhotwgp(spadx+igp)
2930            do ig=1,npwc
2931              otw=DBLE(otq(ig,igp)) !in principle otw -> otw - ieta
2932              num = botsq(ig,igp)*rhotwgdp_igp
2933              den=omegame0i_io-otw
2934              if (den**2>zcut**2) then
2935                ket(spadc+ig,ios)=ket(spadc+ig,ios) + num/(den*otw)*(one-theta_mu_minus_e0i)
2936              else
2937                ket(spadc+ig,ios)=ket(spadc+ig,ios) + &
2938 &               num*CMPLX(den,twofm1_zcut)/((den**2+twofm1_zcut**2)*otw)*(one-theta_mu_minus_e0i)
2939              end if
2940            end do !ig
2941          end do !igp
2942        end do !ios
2943      end if !not fully occupied
2944 
2945    end do !ispinor
2946 
2947 !$omp parallel workshare
2948    ket=ket*half
2949 !$omp end parallel workshare
2950 
2951  CASE (PPM_LINDEN_HORSH,PPM_ENGEL_FARID)
2952    ABI_CHECK(nspinor == 1, "nspinor/=1 not allowed")
2953 
2954    botsq => PPm%bigomegatwsq_qbz%vals !(1:npwc,1:PPm%dm2_botsq)
2955    otq   => PPm%omegatw_qbz%vals      !(1:npwc,1:PPm%dm2_otq)
2956    eig   => PPm%eigpot_qbz%vals       !(1:PPm%dm_eig,1:PPm%dm_eig)
2957 
2958    ! rho-twiddle(G) is formed, introduce rhotwgdpcc, for speed reason
2959    ABI_MALLOC(rhotwgdpcc,(npwx))
2960 
2961    ff=theta_mu_minus_e0i      ! occupation number f (include poles if ...)
2962    twofm1=two*ff-one          ! 2f-1
2963    twofm1_zcut=twofm1*zcut
2964    rhotwgdpcc(:)=CONJG(rhotwgp(:))
2965 
2966    do ios=1,nomega
2967      omegame0i_io=omegame0i(ios)
2968      ct=czero_gw
2969      do ii=1,npwc ! Loop over the DM bands
2970        num=czero_gw
2971 
2972        SELECT CASE (PPm%model)
2973        CASE (PPM_LINDEN_HORSH)
2974          ! Calculate \beta (eq. 106 pag 47)
2975          do ig=1,npwc
2976            num=num+rhotwgdpcc(ig)*eig(ig,ii)
2977          end do
2978          numf=num*CONJG(num) !MG this means that we cannot do SCGW
2979          numf=numf*botsq(ii,1)
2980 
2981        CASE (PPM_ENGEL_FARID)
2982          do ig=1,npwc
2983            num=num+rhotwgdpcc(ig)*botsq(ig,ii)
2984          end do
2985          numf=num*CONJG(num) !MG this means that we cannot do SCGW
2986 
2987        CASE DEFAULT
2988          MSG_ERROR("Wrong PPm%model")
2989        END SELECT
2990 
2991        !numf=num*CONJG(num) !MG this means that we cannot do SCGW
2992        !if (PPm%model==3) numf=numf*botsq(ii,1)
2993 
2994        otw=DBLE(otq(ii,1)) ! in principle otw -> otw - ieta
2995        den=omegame0i_io+otw*twofm1
2996 
2997        if (den**2>zcut**2) then
2998          inv_den=one/den
2999          ct=ct+numf*inv_den
3000        else
3001          inv_den=one/((den**2+twofm1_zcut**2))
3002          ct=ct+numf*CMPLX(den,twofm1_zcut)*inv_den
3003        end if
3004 
3005      end do !ii DM bands
3006      sigcme(ios)=ct*half
3007    end do !ios
3008    ABI_FREE(rhotwgdpcc)
3009 
3010  CASE DEFAULT
3011    MSG_BUG(sjoin('Wrong PPm%model:',itoa(PPm%model)))
3012  END SELECT
3013 
3014 end subroutine ppm_times_ket

m_ppmodel/ppmodel_t [ Types ]

[ Top ] [ m_ppmodel ] [ Types ]

NAME

 ppmodel_t

FUNCTION

  For the GW part of ABINIT, this structure datatype gathers all
  the information on the Plasmonpole technique used in the calculations

NOTES

  If you modify this datatype, please check there there is no creation/destruction/copy routine,
  declared in another part of ABINIT, that might need to take into account your modification.
  Procedures that should be modified to reflect the changes done in the datatype have the tag @ppmodel_t.

SOURCE

 85  type,public :: ppmodel_t
 86 
 87    !integers
 88    integer :: dm2_botsq
 89    ! =npwc if ppmodel=1,2; =1 if ppmodel=3,4
 90 
 91    integer :: dm_eig
 92    ! =npwc if ppmodel=3;   =0 if ppmodel=1,2,4
 93 
 94    integer :: dm2_otq
 95    ! =npwc if ppmodel=1,2; =1 if ppmodel=3,4
 96 
 97    integer :: invalid_freq
 98    ! what to do when PPM frequencies are invalid
 99 
100    integer :: model
101    ! The type of Plasmonpole model
102 
103    integer :: mqmem
104    ! =nqibz if in-core solutions, =0 for out-of-core for which the last dimension in the PPm arrays has size 1.
105 
106    integer :: nqibz
107    ! Number of q-points in the IBZ
108 
109    integer :: npwc
110    ! Number of G vectors in $\tilde \epsilon $
111 
112    integer :: userho
113    ! 1 if the ppmodel requires rho(G).
114 
115    integer :: iq_bz=0
116    ! The index of the q-point in the BZ that is referenced by the internal pointer.
117 
118    integer :: bigomegatwsq_qbz_stat = PPM_ISPOINTER
119    ! Status of bigomegatwsq.
120 
121    integer :: omegatw_qbz_stat = PPM_ISPOINTER
122    ! Status of omegatw_qbz.
123 
124    integer :: eigpot_qbz_stat = PPM_ISPOINTER
125    ! Status of new_eigpot_qbz_stat.
126 
127    real(dp) :: drude_plsmf
128    ! Drude plasma frequency
129 
130    real(dp) :: force_plsmf=zero
131    ! Force plasma frequency to that set by ppmfreq (not used if set zero).
132 
133    ! arrays
134    logical,allocatable :: keep_q(:)
135    ! .TRUE. if the ppmodel tables for this q in the IBZ are kept in memory.
136 
137    integer,allocatable :: has_q(:)
138    ! Flag defining the status of the tables for the different q. See the PPM_TAB flags.
139 
140    type(array2_gwpc_t),pointer :: bigomegatwsq_qbz   => null()
141    ! (Points|Stores) the symmetrized plasmon pole parameters $\tilde\Omega^2_{G Gp}(q_bz)$.
142 
143    type(array2_gwpc_t),pointer :: omegatw_qbz  => null()
144    ! (Points|Stores) the symmetrized plasmon pole parameters $\tilde\omega_{G Gp}(q_bz)$.
145 
146    type(array2_gwpc_t),pointer :: eigpot_qbz   => null()
147    ! (Points|Stores) the eigvectors of the symmetrized inverse dielectric matrix
148 
149    type(array2_gwpc_t),allocatable :: bigomegatwsq(:)
150    ! bigomegatwsq(nqibz)%value(npwc,dm2_botsq)
151    ! Plasmon pole parameters $\tilde\Omega^2_{G Gp}(q)$.
152 
153    type(array2_gwpc_t),allocatable :: omegatw(:)
154    ! omegatw(nqibz)%value(npwc,dm2_otq)
155    ! Plasmon pole parameters $\tilde\omega_{G Gp}(q)$.
156 
157    type(array2_gwpc_t),allocatable :: eigpot(:)
158    ! eigpot(nqibz)%value(dm_eig,dm_eig)
159    ! Eigvectors of the symmetrized inverse dielectric matrix
160 
161  end type ppmodel_t
162 
163 
164  public :: ppm_get_qbz              ! Symmetrize the PPm parameters in the BZ.
165  public :: ppm_nullify              ! Nullify all pointers
166  public :: ppm_init                 ! Initialize dimensions and pointers
167  public :: ppm_free                 ! Destruction method.
168  public :: setup_ppmodel            ! Main Driver
169  public :: new_setup_ppmodel        ! Main Driver
170  public :: getem1_from_PPm          ! Reconstruct e^{-1}(w) from PPm.
171  public :: getem1_from_PPm_one_ggp  ! Reconstruct e^{-1}(w) from PPm for one G,G' pair
172  public :: get_ppm_eigenvalues
173  public :: calc_sig_ppm             ! Matrix elements of the self-energy with ppmodel.
174  public :: ppm_times_ket            ! Matrix elements of the self-energy with ppmodel.
175  public :: ppm_symmetrizer
176  public :: cqratio

m_ppmodel/setup_ppmodel [ Functions ]

[ Top ] [ m_ppmodel ] [ Functions ]

NAME

 setup_ppmodel

FUNCTION

  Initialize some values of several arrays of the PPm datastructure
  that are used in case of plasmonpole calculations
  This is a wrapper around different plasmonpole routines.

INPUTS

  Cryst<crystal_t>=Info on the unit cell and crystal symmetries.
  Qmesh<kmesh_t>=the q-mesh used for the inverse dielectric matrix
    %nibz=number of irreducible q-points
    %ibz(3,%nibz)=the irred q-point
  npwe=number of G vectors for the correlation part
  nomega=number of frequencies in $\epsilon^{-1}$
  omega=frequencies in epsm1
  epsm1=the inverse dielctric matrix
  ngfftf(18)=contain all needed information about the 3D fine FFT mesh, see ~abinit/doc/variables/vargs.htm#ngfft
  gprimd(3,3)=dimensional primitive translations for reciprocal space ($\textrm{bohr}^{-1}$)
  nfftf=the number of points in the FFT mesh (for this processor)
  rhor_tot(nfftf)=the total charge in real space
  PPm<ppmodel_t>:
    %ppmodel=the type of  plasmonpole model

OUTPUT

SIDE EFFECTS

  PPm<ppmodel_t>:
  == if ppmodel 1 or 2 ==
   %omegatw and %bigomegatwsq
  == if ppmodel 3 ==
   %omegatw, %bigomegatwsq and %eigpot
  == if ppmodel 4 ==
   %omegatw and %bigomegatwsq

NOTES

 * FFT parallelism not implemented.
 * TODO: rhor_tot should be replaced by rhog_tot

PARENTS

      calc_sigc_me,mrgscr,sigma

CHILDREN

SOURCE

779 subroutine setup_ppmodel(PPm,Cryst,Qmesh,npwe,nomega,omega,epsm1,nfftf,gvec,ngfftf,rhor_tot,&
780 & iqiA) !Optional
781 
782 
783 !This section has been created automatically by the script Abilint (TD).
784 !Do not modify the following lines by hand.
785 #undef ABI_FUNC
786 #define ABI_FUNC 'setup_ppmodel'
787 !End of the abilint section
788 
789  implicit none
790 
791 !Arguments ------------------------------------
792 !scalars
793  integer,intent(in) :: nfftf,npwe,nomega
794  integer,intent(in),optional :: iqiA
795  type(kmesh_t),intent(in) :: Qmesh
796  type(crystal_t),intent(in) :: Cryst
797  type(ppmodel_t),intent(inout) :: PPm
798 !arrays
799  integer,intent(in) :: gvec(3,npwe),ngfftf(18)
800  real(dp),intent(in) :: rhor_tot(nfftf)
801  complex(dpc),intent(in) :: omega(nomega)
802  complex(gwpc),intent(in) :: epsm1(:,:,:,:)
803 
804 !Local variables-------------------------------
805 !scalars
806  integer :: nqiA,iq_ibz
807  real(dp) :: n_at_G_zero
808  logical :: single_q
809  character(len=500) :: msg
810 !scalars
811  real(dp) :: qpt(3)
812 ! *************************************************************************
813 
814  DBG_ENTER("COLL")
815 
816  !@ppmodel_t
817  !
818  ! === if iqiA is present, then consider only one qpoint to save memory ===
819  ! * This means the object has been already initialized
820  nqiA=Qmesh%nibz; single_q=.FALSE.
821  if (PRESENT(iqiA)) then
822    nqiA=1; single_q=.TRUE.
823  end if
824  !
825  ! Allocate plasmonpole parameters
826  ! TODO ppmodel==1 by default, should be set to 0 if AC and CD
827  SELECT CASE (PPm%model)
828 
829  CASE (PPM_NONE)
830    MSG_COMMENT(' Skipping Plasmompole model calculation')
831 
832  CASE (PPM_GODBY_NEEDS)
833    ! * Note: the q-dependency enters only through epsilon^-1.
834    do iq_ibz=1,nqiA
835      call cppm1par(npwe,nomega,omega,PPm%drude_plsmf,&
836 &       epsm1(:,:,:,iq_ibz),PPm%omegatw(iq_ibz)%vals,PPm%bigomegatwsq(iq_ibz)%vals)
837    end do
838 
839  CASE (PPM_HYBERTSEN_LOUIE)
840    do iq_ibz=1,nqiA
841      qpt = Qmesh%ibz(:,iq_ibz); if (single_q) qpt=Qmesh%ibz(:,iqiA)
842 
843      call cppm2par(qpt,npwe,epsm1(:,:,1,iq_ibz),ngfftf,gvec,Cryst%gprimd,rhor_tot,nfftf,Cryst%gmet,&
844 &      PPm%bigomegatwsq(iq_ibz)%vals,PPm%omegatw(iq_ibz)%vals,PPm%invalid_freq)
845    end do
846    !
847    ! Quick-and-dirty change of the plasma frequency. Never executed in standard runs.
848    if (PPm%force_plsmf>tol6) then ! Integrate the real-space density
849       n_at_G_zero = SUM(rhor_tot(:))/nfftf
850       ! Change the prefactor
851       write(msg,'(2(a,es16.8))') 'Forced ppmfreq:',PPm%force_plsmf*Ha_eV,' nelec/ucvol:',n_at_G_zero
852       MSG_WARNING(msg)
853       PPm%force_plsmf = (PPm%force_plsmf**2)/(four_pi*n_at_G_zero)
854       do iq_ibz=1,PPm%nqibz
855         PPm%bigomegatwsq(iq_ibz)%vals = PPm%force_plsmf * PPm%bigomegatwsq(iq_ibz)%vals
856         PPm%omegatw(iq_ibz)%vals      = PPm%force_plsmf * PPm%omegatw(iq_ibz)%vals
857       end do
858       write(msg,'(a,es16.8)') 'Plasma frequency forced in HL ppmodel, new prefactor is:',PPm%force_plsmf
859       MSG_WARNING(msg)
860    end if
861 
862  CASE (PPM_LINDEN_HORSH) ! TODO Check better double precision, this routine is in a messy state
863    do iq_ibz=1,nqiA
864      qpt = Qmesh%ibz(:,iq_ibz); if (single_q) qpt=Qmesh%ibz(:,iqiA)
865 
866      call cppm3par(qpt,npwe,epsm1(:,:,1,iq_ibz),ngfftf,gvec,Cryst%gprimd,rhor_tot,nfftf,&
867 &      PPm%bigomegatwsq(iq_ibz)%vals,PPm%omegatw(iq_ibz)%vals(:,1),PPm%eigpot(iq_ibz)%vals)
868    end do
869 
870  CASE (PPM_ENGEL_FARID)  ! TODO Check better double precision, this routine is in a messy state
871    do iq_ibz=1,nqiA
872      qpt = Qmesh%ibz(:,iq_ibz); if (single_q) qpt=Qmesh%ibz(:,iqiA)
873      if ((ALL(ABS(qpt)<1.0e-3))) qpt = GW_Q0_DEFAULT ! FIXME
874 
875      call cppm4par(qpt,npwe,epsm1(:,:,1,iq_ibz),ngfftf,gvec,Cryst%gprimd,rhor_tot,nfftf,&
876 &      PPm%bigomegatwsq(iq_ibz)%vals,PPm%omegatw(iq_ibz)%vals(:,1))
877    end do
878 
879  CASE DEFAULT
880    MSG_BUG(sjoin('Wrong PPm%model:',itoa(PPm%model)))
881  END SELECT
882 
883  DBG_EXIT("COLL")
884 
885 end subroutine setup_ppmodel