energy Module Reference

Functions/Subroutines
subroutine	total_e ()
subroutine	kinetic (Etot)
subroutine	hartree (Etot)
subroutine	atomic (Etot)
subroutine	ewald (Etot)
subroutine	xc (Etot)

Function/Subroutine Documentation

total_e()

subroutine energy::total_e

Definition at line 30 of file energy.F90.

    !
    use constant, only : htr2ev
    !
    real(8) :: Etot
    !
    etot = 0.0d0
    !
    write(*,*) "  Energy per u.c. [eV]"
    !
    call kinetic(etot)
    call hartree(etot)
    call atomic(etot)
    call ewald(etot)
    call xc(etot)
    !
    write(*,*) "    Total energy : ", etot*htr2ev
    !

References atomic(), ewald(), hartree(), constant::htr2ev, kinetic(), and xc().

Referenced by pwdft().

Here is the call graph for this function:

Here is the caller graph for this function:

kinetic()

subroutine energy::kinetic

(

real(8), intent(inout)

Etot

)

Definition at line 50 of file energy.F90.

    !
    use constant, only : htr2ev
    use atm_spec, only : bvec
    use k_point, only : nk, kvec, kgrd
    use kohn_sham, only : nbnd, ef, eval, evec
    use gvec, only : g_wf
    use libtetrabz, only : libtetrabz_occ
    !
    real(8),intent(inout) :: Etot
    !
    integer :: ik, ibnd, ipw
    real(8) :: occ(nbnd,nk), kgv(3), gg05(g_wf%npw), Ekin
    !
    eval(1:nbnd,1:nk) = eval(1:nbnd,1:nk) - ef
    call libtetrabz_occ(2,bvec,nbnd,kgrd,eval,kgrd,occ)
    eval(1:nbnd,1:nk) = eval(1:nbnd,1:nk) + ef
    !
    ekin = 0.0d0
    do ik = 1, nk
       !
       do ipw = 1, g_wf%npw
          kgv(1:3) = kvec(1:3,ik) + matmul(bvec(1:3,1:3), dble(g_wf%mill(1:3,g_wf%map(ipw))))
          gg05(ipw) = 0.5d0 * dot_product(kgv,kgv)
       end do
       !
       do ibnd = 1, nbnd
          ekin = ekin + occ(ibnd,ik) &
          &           * dble(sum(conjg(evec(1:g_wf%npw,ibnd,ik)) &
          &                          * evec(1:g_wf%npw,ibnd,ik)  &
          &                          * gg05(1:g_wf%npw        )  ))
       end do
    end do
    !
    ! Spin
    !
    ekin = ekin * 2.0d0
    !
    write(*,*) "    Kinetic energy : ", ekin*htr2ev
    etot = etot + ekin
    !

References atm_spec::bvec, kohn_sham::ef, kohn_sham::eval, kohn_sham::evec, gvec::g_wf, constant::htr2ev, k_point::kgrd, k_point::kvec, kohn_sham::nbnd, and k_point::nk.

Referenced by total_e().

Here is the caller graph for this function:

hartree()

subroutine energy::hartree

(

real(8), intent(inout)

Etot

)

Definition at line 93 of file energy.F90.

    !
    use constant, only : pi, htr2ev
    use gvec, only : g_rh
    use atm_spec, only : bvec, vcell
    use fftw_wrapper, only : fft_r2g
    use rho_v, only : rho
    !
    real(8),intent(inout) :: Etot
    !
    integer :: ir
    real(8) :: g3(3), EH
    complex(8) :: rhog(g_rh%nr)
    !
    call fft_r2g(rho, rhog)
    !
    ! G = 0 : Compensation to the ionic potential 
    !
    eh = 0.0d0
    do ir = 2, g_rh%nr
       g3(1:3) = matmul(bvec(1:3,1:3), dble(g_rh%mill(1:3,ir)))
       eh = eh + 4.0d0 * pi * dble(conjg(rhog(ir))*rhog(ir)) &
       &       / dot_product(g3(1:3),g3(1:3))
    end do
    eh = eh * 0.5d0 * vcell
    !
    write(*,*) "    Hartree energy : ", eh*htr2ev
    etot = etot + eh
    !

References atm_spec::bvec, fftw_wrapper::fft_r2g(), gvec::g_rh, constant::htr2ev, constant::pi, rho_v::rho, and atm_spec::vcell.

Referenced by total_e().

Here is the call graph for this function:

Here is the caller graph for this function:

atomic()

subroutine energy::atomic

(

real(8), intent(inout)

Etot

)

Definition at line 124 of file energy.F90.

    !
    use constant, only : htr2ev
    use rho_v, only : rho, vps
    use gvec, only : g_rh
    use atm_spec, only : vcell
    !
    real(8),intent(inout) :: Etot
    !
    real(8) :: Eps
    !
    eps = dot_product(rho(1:g_rh%nr), vps(1:g_rh%nr)) * vcell / dble(g_rh%nr)
    !
    write(*,*) "    rho*V : ", eps*htr2ev
    etot = etot + eps
    !

References gvec::g_rh, constant::htr2ev, rho_v::rho, atm_spec::vcell, and rho_v::vps.

Referenced by total_e().

Here is the caller graph for this function:

ewald()

subroutine energy::ewald

(

real(8), intent(inout)

Etot

)

Definition at line 142 of file energy.F90.

    !
    use constant, only : htr2ev, pi
    use atm_spec, only : avec, bvec, vcell, atm, spec, nat, nelec
    !
    real(8),intent(inout) :: Etot
    !
    integer :: nmax3(3), iat, jat, i1, i2, i3, ii
    real(8) :: Eew, eta, cutoff, dtau(3), ZZ, gv(3), g2, phase, norm, rv(3)
    !
    eew = 0.0d0
    !
    eta = 0.0d0
    do ii = 1, 3
       norm = sqrt(dot_product(bvec(1:3,ii), bvec(1:3,ii)))
       eta = eta + norm
    end do
    eta = eta / 3.0d0 / sqrt(2.0d0) / sqrt(2.0d0*pi)
    write(*,*) "      Ewald eta [Bohr^-1] : ", eta
    !
    ! Reciprocal space
    !
    cutoff = 10.0d0 * eta * 2.0d0
    do ii = 1, 3
       norm = sqrt(dot_product(avec(1:3,ii), avec(1:3,ii)))
       nmax3(ii) = ceiling(cutoff*norm/(2.0d0*pi))
    end do
    !
    write(*,*) "      Ewald grid (G) : ", nmax3(1:3)
    !
    do iat = 1, nat
       do jat = 1, nat
          !
          dtau(1:3) = atm(iat)%pos(1:3) - atm(jat)%pos(1:3)
          zz = spec(atm(iat)%ityp)%Zion * spec(atm(jat)%ityp)%Zion
          !
          do i3 = -nmax3(3), nmax3(3)
             do i2 = -nmax3(2), nmax3(2)
                do i1 = -nmax3(1), nmax3(1)
                   !
                   if(all((/i1,i2,i3/)==0)) cycle
                   !
                   gv(1:3) = matmul(bvec(1:3,1:3), dble((/i1,i2,i3/)))
                   g2 = dot_product(gv,gv)
                   phase = 2.0d0*pi*dot_product(dtau(1:3), dble((/i1,i2,i3/)))
                   !
                   eew = eew &
                   &   + 0.5d0 * zz * cos(phase) * 4.0d0 * pi * exp(-g2/(4.0d0*eta**2)) / (vcell*g2)
                   !
                end do
             end do
          end do
       end do
    end do
    !
    eew = eew - 0.5d0 * pi * nelec**2 / (vcell*eta**2)
    !
    ! Real space
    !
    cutoff = 10.0d0 / eta
    do ii = 1, 3
       norm = sqrt(dot_product(bvec(1:3,ii), bvec(1:3,ii)))
       nmax3(ii) = ceiling(cutoff*norm/(2.0d0*pi))
    end do
    !
    write(*,*) "      Ewald grid (R) : ", nmax3(1:3)
    !
    do iat = 1, nat
       do jat = 1, nat
          !
          dtau(1:3) = atm(iat)%pos(1:3) - atm(jat)%pos(1:3)
          zz = spec(atm(iat)%ityp)%Zion * spec(atm(jat)%ityp)%Zion
          !
          do i3 = -nmax3(3), nmax3(3)
             do i2 = -nmax3(2), nmax3(2)
                do i1 = -nmax3(1), nmax3(1)
                   !
                   if(all((/i1,i2,i3/)==0) .and. iat == jat) cycle
                   !
                   rv(1:3) = matmul(avec(1:3,1:3), dtau(1:3) + dble((/i1,i2,i3/)))
                   norm = sqrt(dot_product(rv(1:3), rv(1:3)))
                   !
                   eew = eew &
                   &   + 0.5d0 * zz * erfc(norm*eta) / norm
                   !
                end do
             end do
          end do
       end do
       !
       eew = eew - eta * spec(atm(iat)%ityp)%Zion**2 / sqrt(pi)
       !
    end do
    !
    write(*,*) "    Ewald energy : ", eew*htr2ev
    etot = etot + eew
    !

References atm_spec::atm, atm_spec::avec, atm_spec::bvec, constant::htr2ev, atm_spec::nat, atm_spec::nelec, constant::pi, atm_spec::spec, and atm_spec::vcell.

Referenced by total_e().

Here is the caller graph for this function:

xc()

subroutine energy::xc

(

real(8), intent(inout)

Etot

)

Definition at line 241 of file energy.F90.

    !
    use constant, only : htr2ev, pi
    use gvec, only : g_rh
    use rho_v, only : rho
    use atm_spec, only : vcell
    !
    real(8),intent(inout) :: Etot
    !
    integer :: ir
    real(8) :: rs, exc0, Exc
    !
    exc = 0.0d0
    do ir = 1, g_rh%nr
       !
       if(rho(ir) > 1.0d-10) then
          !
          rs = (3.0d0 / (4.0d0 * pi * rho(ir)))**(1.0d0/3.0d0)
          !
          if(rs < 1.0d0) then
             exc0 = -3.0d0 / (4.0d0*pi) * (9.0d0*pi/4)**(1.0d0/3.0d0) / rs &
             &    - 0.0480d0 + 0.031d0*log(rs) - 0.0116d0*rs + 0.0020d0*rs*log(rs)
          else
             exc0 = -3.0d0 / (4.0d0*pi) * (9.0d0*pi/4)**(1.0d0/3.0d0) / rs &
             &   -  0.1423d0 / (1.0d0 + 1.0529d0*sqrt(rs) + 0.3334d0*rs)
          end if
          !
          exc = exc + exc0 * rho(ir)
          !
       end if
       !
    end do
    exc = exc * vcell / dble(g_rh%nr)
    !
    write(*,*) "    XC energy : ", exc*htr2ev
    etot = etot + exc
    !

References gvec::g_rh, constant::htr2ev, constant::pi, rho_v::rho, and atm_spec::vcell.

Referenced by total_e().

Here is the caller graph for this function: