prg_xlkernel_mod.F90

Go to the documentation of this file.

 
module prg_xlkernel_mod
 
  use bml
  use prg_kernelparser_mod
  use prg_extras_mod
 
  implicit none
 
  private
 
  integer, parameter :: dp = kind(1.0d0)
 
  type, public :: xlk_type
 
    character(20) :: kerneltype
 
    integer       :: verbose, nrank
 
    real(dp)      :: scalecoeff
 
  end type xlk_type
 
  public :: prg_parse_xlkernel, prg_rank1, prg_fermi
  public :: prg_v_kernel_fermi, prg_kernel_fermi_full
  private :: prg_mmult, prg_eig, prg_get_deriv_finite_temp, prg_inv
 
contains
 
  subroutine prg_parse_xlkernel(input,filename)
 
    type(xlk_type), intent(inout) :: input
    integer, parameter :: nkey_char = 1, nkey_int = 2, nkey_re = 1, nkey_log = 1
    character(len=*) :: filename
 
    !Library of keywords with the respective defaults.
    character(len=50), parameter :: keyvector_char(nkey_char) = [character(len=100) :: &
         'KernelType=']
    character(len=100) :: valvector_char(nkey_char) = [character(len=100) :: &
         'Full']
 
    character(len=50), parameter :: keyvector_int(nkey_int) = [character(len=50) :: &
         'Verbose=','Nrank=']
    integer :: valvector_int(nkey_int) = (/ &
         0, 1 /)
 
    character(len=50), parameter :: keyvector_re(nkey_re) = [character(len=50) :: &
         'ScaleCoeff=']
    real(dp) :: valvector_re(nkey_re) = (/&
         0.1 /)
 
    character(len=50), parameter :: keyvector_log(nkey_log) = [character(len=100) :: &
         'MixerON=']
    logical :: valvector_log(nkey_log) = (/&
         .true./)
 
    !Start and stop characters
    character(len=50), parameter :: startstop(2) = [character(len=50) :: &
         'XLKERNEL{', '}']
 
    call prg_parsing_kernel(keyvector_char,valvector_char&
         ,keyvector_int,valvector_int,keyvector_re,valvector_re,&
         keyvector_log,valvector_log,trim(filename),startstop)
 
    !Characters
    input%kerneltype = valvector_char(1)
 
    !Integers
    input%verbose = valvector_int(1)
    input%nrank = valvector_int(2)
 
    !Logicals
 
    !Reals
    input%scalecoeff = valvector_re(1)
 
  end subroutine prg_parse_xlkernel
 
  subroutine prg_fermi(D0,QQ,ee,gap,Fe_vec,mu0,H,Z,Nocc,T,OccErrLim,MaxIt,HDIM)
 
    integer, parameter             :: prec = 8
    integer(PREC), intent(in)      :: hdim, nocc
    real(prec), parameter          :: one = 1.d0, two = 2.d0, zero = 0.d0
    real(prec), parameter          :: kb = 8.61739d-5 ! eV/K, kB = 6.33366256e-6 Ry/K, kB = 3.166811429e-6 Ha/K
    real(prec), intent(in)         :: h(hdim,hdim), z(hdim,hdim)
    real(prec)                     :: h0(hdim,hdim), fe(hdim,hdim)
    real(prec), intent(out)        :: d0(hdim,hdim), qq(hdim,hdim), ee(hdim), fe_vec(hdim), gap
    real(prec)                     :: x(hdim,hdim)
    real(prec), intent(in)         :: t, occerrlim
    real(prec), intent(inout)      :: mu0
    integer(PREC), intent(in)      :: maxit
    real(prec)                     :: occerr, occ, docc, mu_0, beta, occ_i
    integer(PREC)                  :: i,iter
 
    !! Orthogonalize H with inverse overlap factor Z, where Z'SZ = I
    call prg_mmult(one,z,h,zero,x,'T','N',hdim)  !
    call prg_mmult(one,x,z,zero,h0,'N','N',hdim)   ! H0 = Z'*H*Z
    call eig(h0,qq,ee,'V',hdim)
    !! QQ eigenvectors of H0 and ee the eigenvalues
 
    occerr = one
    fe = zero
    beta = 1.d0/(kb*t)    ! Temp in Kelvin
    gap = ee(nocc+1)-ee(nocc)
    if (mu0 == 0.d0) then
      mu0 = 0.5d0*(ee(nocc)+ee(nocc+1)) ! Initial estimate of chemical potential
    endif
    iter = 0
    do while (occerr > 1e-9)
      iter = iter + 1
      occ = zero
      docc = zero
      do i = 1,hdim
        occ_i = one/(exp(beta*(ee(i)-mu0))+one)
        fe(i,i) = occ_i
        fe_vec(i) = occ_i
        occ = occ + occ_i
        docc = docc + beta*occ_i*(one-occ_i)
      enddo
      occerr = abs(nocc-occ)
      if (abs(occerr) > 1e-9) then
        mu0 = mu0 + (nocc-occ)/docc  ! Adjust occupation by shifting mu
      endif
      if(iter.gt.maxit) then
        occerr = zero
      endif
    enddo
    call prg_mmult(one,qq,fe,zero,x,'N','N',hdim)
    call prg_mmult(one,x,qq,zero,d0,'N','T',hdim)
 
  end subroutine prg_fermi
 
  subroutine prg_kernel_fermi_full(KK,JJ,D0,mu0,mu1,T,RX,RY,RZ,LBox,Hubbard_U,Element_Type,Nr_atoms, &
       MaxIt,eps, m,HDIM, Max_Nr_Neigh,Coulomb_acc,TIMERATIO,nnRx,nnRy,nnRz,nrnnlist,nnType,H_INDEX_START, &
       H_INDEX_END,H,S,Z,Nocc,Znuc,QQ,ee,Fe_vec)
 
    integer, parameter             :: prec = 8
    integer(PREC), intent(in)      :: nr_atoms, hdim, nocc,max_nr_neigh
    real(prec), parameter          :: one = 1.d0, two = 2.d0, zero = 0.d0
    real(prec), parameter          :: kb = 8.61739d-5 ! eV/K, kB = 6.33366256e-6 Ry/K, kB = 3.166811429e-6 Ha/K
    real(prec), intent(in)         :: coulomb_acc, timeratio
    real(prec)                     :: v(nr_atoms)
    real(prec), intent(in)         :: rx(nr_atoms), ry(nr_atoms), rz(nr_atoms), lbox(3)
    integer(PREC), intent(in)      :: h_index_start(nr_atoms), h_index_end(nr_atoms)
    real(prec), intent(in)         :: znuc(nr_atoms), hubbard_u(nr_atoms)
    real(prec), intent(in)         :: h(hdim,hdim), s(hdim,hdim), z(hdim,hdim)
    real(prec)                     :: h0(hdim,hdim), h1(hdim,hdim)
    real(prec), intent(inout)      :: d0(hdim,hdim)
    real(prec)                     :: x(hdim,hdim), xx(nr_atoms,nr_atoms)
    real(prec), intent(in)         :: t, eps
    real(prec), intent(inout)      :: mu0, mu1
    integer(PREC), intent(in)      :: m, maxit
    character(10), intent(in)      :: element_type(nr_atoms)
    integer(PREC), intent(in)      :: nrnnlist(nr_atoms), nntype(nr_atoms,max_nr_neigh)
    real(prec), intent(in)         :: nnrx(nr_atoms,max_nr_neigh)
    real(prec), intent(in)         :: nnry(nr_atoms,max_nr_neigh), nnrz(nr_atoms,max_nr_neigh)
    real(prec)                     :: coulomb_pot_real(nr_atoms), coulomb_pot(nr_atoms)
    real(prec)                     :: coulomb_pot_real_i, coulomb_pot_k(nr_atoms)
    real(prec)                     :: coulomb_pot_real_dq_v(nr_atoms), coulomb_pot_dq_v(nr_atoms)
    real(prec)                     :: coulomb_force_real_i(3),coulomb_force_k(3,nr_atoms)
    real(prec)                     :: d_dq_v(hdim,hdim), h_dq_v(hdim,hdim), dq_v(nr_atoms)
    real(prec)                     :: dq_dv(nr_atoms)
    real(prec), intent(out)        :: kk(nr_atoms,nr_atoms),jj(nr_atoms,nr_atoms)
    real(prec), intent(in)         :: qq(hdim,hdim), ee(hdim), fe_vec(hdim)
    integer(PREC)                  :: i,j,k, iter, mm
 
    xx = zero
    x = zero
    h1 = zero
    coulomb_pot_dq_v = zero
    coulomb_pot_k = zero
    dq_v = zero
    h_dq_v = zero
 
    do j = 1,nr_atoms  ! TRIVIAL OPEN_MP OR MPI PARALLELISM
      dq_v(j) = one
      do i = 1,nr_atoms
        call ewald_real_space(coulomb_pot_real_i,coulomb_force_real_i,i,rx,ry,rz,lbox, &
             dq_v,hubbard_u,element_type,nr_atoms,coulomb_acc,timeratio,nnrx,nnry,nnrz,nrnnlist,nntype,hdim,max_nr_neigh)
        coulomb_pot_real(i) = coulomb_pot_real_i
      enddo
      call ewald_k_space(coulomb_pot_k,coulomb_force_k,rx,ry,rz,lbox,dq_v,nr_atoms,coulomb_acc,timeratio,hdim,max_nr_neigh)
      coulomb_pot_dq_v = coulomb_pot_real+coulomb_pot_k
 
      h_dq_v = zero
      do i = 1,nr_atoms
        do k = h_index_start(i),h_index_end(i)
          h_dq_v(k,k) = hubbard_u(i)*dq_v(i) + coulomb_pot_dq_v(i)
        enddo
      enddo
      call prg_mmult(one,s,h_dq_v,zero,x,'N','N',hdim)
      call prg_mmult(one,h_dq_v,s,one,x,'N','T',hdim)
      h_dq_v =  (one/two)*x
 
      call prg_mmult(one,z,h,zero,x,'T','N',hdim)  !
      call prg_mmult(one,x,z,zero,h0,'N','N',hdim)   ! H0 = Z'*H*Z
 
      call prg_mmult(one,z,h_dq_v,zero,x,'T','N',hdim)  !
      call prg_mmult(one,x,z,zero,h1,'N','N',hdim)   ! H1 = Z'*H_dq_v*Z
 
      call get_deriv_finite_temp(d_dq_v,h0,h1,nocc,t,qq,ee,fe_vec,mu0,eps,hdim)
 
      call prg_mmult(two,z,d_dq_v,zero,x,'N','N',hdim)
      call prg_mmult(one,x,z,zero,d_dq_v,'N','T',hdim)
 
      call prg_mmult(one,d_dq_v,s,zero,x,'N','N',hdim)
      dq_dv = zero
      do i = 1, nr_atoms
        do k = h_index_start(i), h_index_end(i)
          dq_dv(i) = dq_dv(i) + x(k,k)
        enddo
        jj(i,j) = dq_dv(i)
      enddo
      dq_v = zero
    enddo
    xx = jj
    do i = 1,nr_atoms
      xx(i,i) = xx(i,i) - one
    enddo
    call inv(xx,kk,nr_atoms)
 
  end subroutine prg_kernel_fermi_full
 
  subroutine prg_v_kernel_fermi(D0,dq_dv,v,mu0,mu1,T,RX,RY,RZ,LBox,Hubbard_U,Element_Type,Nr_atoms, &
       MaxIt,eps, m,HDIM, Max_Nr_Neigh,Coulomb_acc,TIMERATIO,nnRx,nnRy,nnRz,nrnnlist,nnType,H_INDEX_START, &
       H_INDEX_END,H,S,Z,Nocc,Znuc,QQ,ee,Fe_vec)
 
    integer, parameter             :: prec = 8
    integer(PREC), intent(in)      :: nr_atoms, hdim, nocc,max_nr_neigh
    real(prec), parameter          :: one = 1.d0, two = 2.d0, zero = 0.d0
    real(prec), parameter          :: kb = 8.61739d-5 ! eV/K, kB = 6.33366256e-6 Ry/K, kB = 3.166811429e-6 Ha/K
    real(prec), intent(in)         :: coulomb_acc, timeratio, v(nr_atoms)
    real(prec), intent(in)         :: rx(nr_atoms), ry(nr_atoms), rz(nr_atoms), lbox(3)
    integer(PREC), intent(in)      :: h_index_start(nr_atoms), h_index_end(nr_atoms)
    real(prec), intent(in)         :: znuc(nr_atoms), hubbard_u(nr_atoms)
    real(prec), intent(in)         :: h(hdim,hdim), s(hdim,hdim), z(hdim,hdim)
    real(prec)                     :: h0(hdim,hdim), h1(hdim,hdim)
    real(prec), intent(inout)      :: d0(hdim,hdim)
    real(prec)                     :: x(hdim,hdim)
    real(prec), intent(in)         :: t, eps
    real(prec), intent(inout)      :: mu0, mu1
    integer(PREC), intent(in)      :: m, maxit
    character(10), intent(in)      :: element_type(nr_atoms)
    integer(PREC), intent(in)      :: nrnnlist(nr_atoms), nntype(nr_atoms,max_nr_neigh)
    real(prec), intent(in)         :: nnrx(nr_atoms,max_nr_neigh)
    real(prec), intent(in)         :: nnry(nr_atoms,max_nr_neigh), nnrz(nr_atoms,max_nr_neigh)
    real(prec)                     :: coulomb_pot_real(nr_atoms), coulomb_pot(nr_atoms)
    real(prec)                     :: coulomb_pot_real_i, coulomb_pot_k(nr_atoms)
    real(prec)                     :: coulomb_pot_real_dq_v(nr_atoms), coulomb_pot_dq_v(nr_atoms)
    real(prec)                     :: coulomb_force_real_i(3),coulomb_force_k(3,nr_atoms)
    real(prec)                     :: d_dq_v(hdim,hdim), h_dq_v(hdim,hdim), dq_v(nr_atoms)
    real(prec), intent(out)        :: dq_dv(nr_atoms)
    real(prec), intent(in)         :: qq(hdim,hdim), ee(hdim), fe_vec(hdim)
    integer(PREC)                  :: i,j,k, iter, mm
 
    dq_dv = zero
    dq_v = v/prg_norm2(v)
 
    do i = 1,nr_atoms
      call ewald_real_space(coulomb_pot_real_i,coulomb_force_real_i,i,rx,ry,rz,lbox, &
           dq_v,hubbard_u,element_type,nr_atoms,coulomb_acc,timeratio,nnrx,nnry,nnrz,nrnnlist,nntype,hdim,max_nr_neigh)
      coulomb_pot_real(i) = coulomb_pot_real_i
    enddo
    call ewald_k_space(coulomb_pot_k,coulomb_force_k,rx,ry,rz,lbox,dq_v,nr_atoms,coulomb_acc,timeratio,hdim,max_nr_neigh)
    coulomb_pot_dq_v = coulomb_pot_real+coulomb_pot_k
 
    h_dq_v = zero
    do i = 1,nr_atoms
      do k = h_index_start(i),h_index_end(i)
        h_dq_v(k,k) = hubbard_u(i)*dq_v(i) + coulomb_pot_dq_v(i)
      enddo
    enddo
    call prg_mmult(one,s,h_dq_v,zero,x,'N','N',hdim)
    call prg_mmult(one,h_dq_v,s,one,x,'N','T',hdim)
    h_dq_v =  (one/two)*x
 
    call prg_mmult(one,z,h,zero,x,'T','N',hdim)  !
    call prg_mmult(one,x,z,zero,h0,'N','N',hdim)   ! H0 = Z'*H*Z
 
    call prg_mmult(one,z,h_dq_v,zero,x,'T','N',hdim)  !
    call prg_mmult(one,x,z,zero,h1,'N','N',hdim)   ! H1 = Z'*H_dq_v*Z
 
    call get_deriv_finite_temp(d_dq_v,h0,h1,nocc,t,qq,ee,fe_vec,mu0,eps,hdim)
 
    call prg_mmult(two,z,d_dq_v,zero,x,'N','N',hdim)
    call prg_mmult(one,x,z,zero,d_dq_v,'N','T',hdim)
 
    call prg_mmult(one,d_dq_v,s,zero,x,'N','N',hdim)
    do i = 1, nr_atoms
      do k = h_index_start(i), h_index_end(i)
        dq_dv(i) = dq_dv(i) + x(k,k)
      enddo
    enddo
 
  end subroutine prg_v_kernel_fermi
 
  subroutine  prg_get_deriv_finite_temp(P1,H0,H1,Nocc,T,Q,ev,fe,mu0,eps,HDIM)
 
    integer, parameter      :: prec = 8
    integer(PREC), intent(in)       :: hdim, nocc
    real(prec), parameter   :: one = 1.d0, zero = 0.d0, two = 2.d0
    real(prec), intent(in)  :: h0(hdim,hdim), h1(hdim,hdim), q(hdim,hdim), ev(hdim), fe(hdim)
    real(prec)              :: x(hdim,hdim), yy(hdim,hdim), fermifun
    real(prec)              :: p02(hdim,hdim), ii(hdim,hdim), xi(hdim,hdim)
    real(prec)              :: dx1(hdim,hdim)
    real(prec)              :: id0(hdim,hdim), t12(hdim,hdim), ddt(hdim,hdim)
    real(prec), intent(out) :: p1(hdim,hdim) ! DM derivative
    real(prec), intent(in)  :: t, eps
    real(prec), intent(inout)  :: mu0
    real(prec)              ::  occerr, trdpdmu, trp0, trp1
    real(prec)              :: beta, cnst, kb, mu1, dh, dm, y, dy, trx, trddt
    integer(PREC)           :: n, cnt, i, j, k,l, ind, ind2, fact
 
    n = hdim
    kb = 8.61739e-5
    beta = 1.d0/(kb*t)     ! Temp in Kelvin
 
    ! T12 = F0_ort_Vecs'*F1_ort*F0_ort_Vecs;
    call prg_mmult(one,q,h1,zero,x,'T','N',hdim)
    call prg_mmult(one,x,q,zero,t12,'N','N',hdim)
 
    ! Divided differences matrix
    ddt = zero
    x = zero
    do i = 1,n
      do j = 1,n
        if (abs(ev(i)-ev(j)) < 1e-4) then
          ! divided difference from Taylor expansion
          y = (ev(i)+ev(j))/2.d0
          if (abs(beta*(y-mu0)) > 500.d0) then
            ddt(i,j) = 0.d0
          else
            ddt(i,j) = -beta*exp(beta*(y-mu0))/(exp(beta*(y-mu0))+1.d0)**2
          endif
        else
          ddt(i,j) = (fe(i)-fe(j))/(ev(i)-ev(j))
        endif
        x(i,j) = ddt(i,j)*t12(i,j)
      end do
    end do
    mu1 = zero
    trx = zero
    trddt = zero
    do i = 1,n
      trx = trx + x(i,i)
      trddt = trddt + ddt(i,i)
    end do
    mu1 = trx/trddt
    do i = 1,n
      x(i,i) = x(i,i)-ddt(i,i)*mu1
    end do
    call prg_mmult(one,q,x,zero,yy,'N','N',hdim)
    call prg_mmult(one,yy,q,zero,p1,'N','T',hdim)
 
  end subroutine prg_get_deriv_finite_temp
 
  subroutine prg_mmult(alpha,A,B,beta,C,TA,TB,HDIM)
 
    integer, parameter          :: prec = 8
    integer(PREC), intent(in)   :: hdim
    real(prec), intent(in)      :: a(hdim,hdim), b(hdim,hdim), alpha, beta
    real(prec), intent(inout)   :: c(hdim,hdim)
    character(1), intent(in)    :: ta, tb
 
    if (prec.eq.4) then
      call sgemm(ta, tb, hdim, hdim, hdim, alpha, &
           a, hdim, b, hdim, beta, c, hdim)
    else
      call dgemm(ta, tb, hdim, hdim, hdim, alpha, &
           a, hdim, b, hdim, beta, c, hdim)
    endif
 
  end subroutine prg_mmult
 
  subroutine prg_eig(A,Q,ee,type,HDIM)
 
    integer, parameter             :: prec = 8
    integer(PREC), intent(in)      :: hdim
    real(prec), intent(in)         :: a(hdim,hdim)
    real(prec), intent(out)        :: q(hdim,hdim), ee(hdim)
    character(1), intent(in)       :: type  ! 'N'(for eigenvaules only) or 'V' (otherwise)
 
    real(prec)    :: nono_evals(hdim), nono_work(1 + 6*hdim + 2*hdim*hdim)
    real(prec)    :: nonotmp(hdim,hdim), z(hdim,hdim)
    integer(PREC) :: nono_iwork(3+5*hdim)
    integer(PREC) :: i, j, info, nono_lwork
    real(prec)    :: invsqrt,err_check,numthresh
 
    nono_lwork = 1 + 6*hdim + 2*hdim*hdim
 
    q = a
    if (prec.eq.4) then
      call ssyev(type, 'U', hdim, q, hdim, ee, nono_work, &
           nono_lwork, info)
    else
      call dsyev('V', 'U', hdim, q, hdim, ee, nono_work, &
           nono_lwork, info)
    endif
 
  end subroutine prg_eig
 
  subroutine prg_inv(X,XI,HDIM)
 
    integer, parameter        :: prec = 8
    integer(PREC), intent(in) :: hdim
    integer(PREC)             :: lda, lwork, m, n, info, ipiv(hdim)
    real(prec), intent(in)    :: x(hdim,hdim)
    real(prec), intent(out)   :: xi(hdim,hdim)
    real(prec)                :: work(hdim*hdim)
 
    lda = hdim
    m = hdim
    n = hdim
    lwork = hdim*hdim
    xi = x
 
    if (prec.eq.4) then
      call sgetrf(m, n, xi, lda, ipiv, info)
      call sgetri(n, xi, n, ipiv, work, lwork, info)
    else
      call dgetrf(m, n, xi, lda, ipiv, info)
      call dgetri(n, xi, n, ipiv, work, lwork, info)
    endif
 
  end subroutine prg_inv
 
  subroutine prg_rank1(verbose)
 
    integer :: i,j,info,s,k,n
    integer :: nint, na, nb
    integer, intent(in) :: verbose
    real(dp), allocatable :: d(:),dl(:),dnewin(:)
    real(dp), allocatable :: dnewout(:)
    real(dp), allocatable :: coef(:,:),b(:),ipiv(:)
 
    write(*,*)"prt_rank1 Verbosity",verbose
 
  end subroutine prg_rank1
 
end module prg_xlkernel_mod