prg_charges_mod.F90

Go to the documentation of this file.

 
module prg_charges_mod
 
  use prg_openfiles_mod
  use bml
  use prg_parallel_mod
  use prg_graph_mod
  use prg_system_mod
 
  implicit none
 
  private
 
  integer, parameter :: dp = kind(1.0d0)
 
  public :: prg_get_charges, prg_get_hscf, prg_get_hscf_v2
 
contains
 
  subroutine prg_get_charges(rho_bml,over_bml,hindex,charges,numel,spindex,mdimin,threshold)
    character(20)                        ::  bml_type, bml_mode
    integer                              ::  i, j, nats, norb
    integer                              ::  mdim
    integer, intent(in)                  ::  hindex(:,:), mdimin, spindex(:)
    real(dp)                             ::  znuc
    real(dp), allocatable                ::  rho_diag(:)
    real(dp), allocatable, intent(inout) ::  charges(:)
    real(dp), intent(in)                 ::  numel(:), threshold
    type(bml_matrix_t)                   ::  aux_bml
    type(bml_matrix_t), intent(inout)    ::  over_bml, rho_bml
 
    norb = bml_get_n(rho_bml)
    bml_type = bml_get_type(rho_bml)
    bml_mode = bml_get_distribution_mode(rho_bml)
    nats = size(hindex,dim=2)
 
    if(mdimin.lt.0)then
      mdim = norb
    else
      mdim = mdimin
    endif
 
    if(.not.allocated(charges)) allocate(charges(nats))
    allocate(rho_diag(norb))
 
    call bml_zero_matrix(bml_type,bml_element_real,dp,norb,mdim,aux_bml, &
         bml_mode)
 
    call bml_multiply(rho_bml,over_bml,aux_bml,1.0_dp,0.0_dp,threshold)
 
#ifdef DO_MPI
    if (getnranks() .gt. 1 .and. &
         bml_get_distribution_mode(aux_bml) == bml_dmode_distributed) then
      call prg_allgatherparallel(aux_bml)
    endif
#endif
 
    call bml_get_diagonal(aux_bml,rho_diag)
 
    do i = 1,nats
      znuc = numel(spindex(i))
      charges(i)=0.0_dp
      do j = hindex(1,i),hindex(2,i)
        charges(i) = charges(i) + rho_diag(j)
      enddo
      charges(i) = charges(i) - znuc
    enddo
 
    deallocate(rho_diag)
    call bml_deallocate(aux_bml)
 
  end subroutine prg_get_charges
 
  subroutine prg_get_hscf(ham0_bml,over_bml,ham_bml,spindex,hindex,hubbardu,charges,&
       coulomb_pot_r,coulomb_pot_k,mdimin,threshold)
    character(20)                      ::  bml_type
    integer                            ::  i, j, nats, norb, mdim
    integer, intent(in)                ::  hindex(:,:), mdimin, spindex(:)
    real(dp), allocatable              ::  coulomb_pot(:), diagonal(:)
    real(dp), intent(in)               ::  charges(:), coulomb_pot_k(:), coulomb_pot_r(:), hubbardu(:)
    real(dp), intent(in)               ::  threshold
    type(bml_matrix_t)                 ::  aux_bml
    type(bml_matrix_t), intent(in)     ::  ham0_bml, over_bml
    type(bml_matrix_t), intent(inout)  ::  ham_bml
 
    nats = size(coulomb_pot_r)
    allocate(coulomb_pot(nats))
 
    norb = bml_get_n(ham0_bml)
    if(hindex(2,nats) .ne. norb)then
      write(*,*)"ERROR in prg_get_hscf. hindex incompatible with ham0"
      stop
    endif
 
    if(mdimin.lt.0)then
      mdim = norb
    else
      mdim = mdimin
    endif
 
    allocate(diagonal(norb))
 
    call bml_copy_new(ham0_bml,ham_bml)
 
    bml_type = bml_get_type(ham_bml)
 
    coulomb_pot = coulomb_pot_k + coulomb_pot_r
 
    call bml_zero_matrix(bml_type,bml_element_real,dp,mdim,norb,aux_bml)
 
    do i = 1,nats
      do j = hindex(1,i),hindex(2,i)
        diagonal(j) = hubbardu(spindex(i))*charges(i) + coulomb_pot(i)
      enddo
    enddo
 
    call bml_set_diagonal(aux_bml,diagonal)
 
    call bml_multiply(over_bml,aux_bml,ham_bml,0.5_dp,1.0_dp,threshold) !  h = h + 0.5*s*h1
 
    call bml_multiply(aux_bml,over_bml,ham_bml,0.5_dp,1.0_dp,threshold) !  h = h + 0.5*h1*s
 
    deallocate(diagonal)
    deallocate(coulomb_pot)
    call bml_deallocate(aux_bml)
 
  end subroutine prg_get_hscf
 
 
  subroutine prg_get_hscf_v2(ham0_bml,over_bml,ham_bml,spindex,hindex,hubbardu,charges,&
       coulomb_pot_r,coulomb_pot_k,mdimin,threshold)
    character(20)                      ::  bml_type
    integer                            ::  i, j, nats, norb, mdim
    integer, intent(in)                ::  hindex(:,:), mdimin, spindex(:)
    real(dp), allocatable              ::  coulomb_pot(:), diagonal(:)
    real(dp), intent(in)               ::  charges(:), coulomb_pot_k(:),coulomb_pot_r(:), hubbardu(:)
    real(dp), intent(in)               ::  threshold
    type(bml_matrix_t)                 ::  aux_bml
    type(bml_matrix_t), intent(in)     ::  ham0_bml, over_bml
    type(bml_matrix_t), intent(inout)  ::  ham_bml
 
    nats = size(coulomb_pot_r)
    allocate(coulomb_pot(nats))
 
    norb = bml_get_n(ham0_bml)
    if(hindex(2,nats) .ne. norb)then
      write(*,*)"ERROR in prg_get_hscf. hindex incompatible with ham0"
      stop
    endif
 
    if(mdimin.lt.0)then
      mdim = norb
    else
      mdim = mdimin
    endif
 
    allocate(diagonal(norb))
 
    call bml_copy_new(ham0_bml,ham_bml)
 
    bml_type = bml_get_type(ham_bml)
 
    coulomb_pot = coulomb_pot_k + coulomb_pot_r
 
    call bml_zero_matrix(bml_type,bml_element_real,dp,mdim,norb,aux_bml)
 
    do i = 1,nats
      do j = hindex(1,i),hindex(2,i)
        diagonal(j) = hubbardu(spindex(i))*charges(i) + coulomb_pot(i)
      enddo
    enddo
 
    call bml_set_diagonal(aux_bml,diagonal)
    call bml_multiply(over_bml,aux_bml,ham_bml,0.5_dp,1.0_dp,threshold) !  h = h+ 0.5*s*h1
 
    call bml_multiply(aux_bml,over_bml,ham_bml,0.5_dp,1.0_dp,threshold) !  h = h+ 0.5*h1*s
 
    deallocate(diagonal)
    deallocate(coulomb_pot)
    call bml_deallocate(aux_bml)
 
  end subroutine prg_get_hscf_v2
 
 
end module prg_charges_mod