prg_pulaycomponent_mod.F90

Go to the documentation of this file.

 
module prg_pulaycomponent_mod
 
  use bml
 
  implicit none
 
  private
 
  integer,parameter :: dp = kind(1.0d0)
 
  public :: prg_pulaycomponent0, prg_pulaycomponentt, prg_get_pulayforce
 
contains
 
  subroutine prg_pulaycomponent0(rho_bml,ham_bml,pcm_bml,threshold,M,&
       &bml_type,verbose)
 
    implicit none
 
    type(bml_matrix_t), intent(in) :: rho_bml
    type(bml_matrix_t), intent(in) :: ham_bml
    type(bml_matrix_t) :: aux_bml
    type(bml_matrix_t), intent(inout) :: pcm_bml
    integer, intent(in) :: m
    integer :: norb, verbose
    real(dp), intent(in) :: threshold
    character(20), intent(in) :: bml_type
 
    if(verbose.eq.2) write(*,*)"In prg_PulayComponent0 ..."
 
    norb = bml_get_n(rho_bml)
 
    call bml_zero_matrix(bml_type,bml_element_real,dp,norb ,norb,aux_bml)
 
    if(bml_get_n(pcm_bml).le.0)then !If pcm is not allocated
      call bml_zero_matrix(bml_type,bml_element_real,dp,norb,norb,pcm_bml)
    else
      call bml_deallocate(pcm_bml) !If pcm is allocated then we set it to 0
      call bml_zero_matrix(bml_type,bml_element_real,dp,norb,norb,pcm_bml)
    endif
 
    call bml_zero_matrix(bml_type,bml_element_real,dp,norb ,norb,pcm_bml)
 
    call bml_multiply(rho_bml, ham_bml, aux_bml, 1.0d0, 1.0d0,threshold) !D*H
 
    call bml_multiply(aux_bml,rho_bml,pcm_bml, 1.0d0, 1.0d0,threshold)  !(D*H)*D
 
    call bml_scale(0.5_dp, pcm_bml)
 
    call bml_deallocate(aux_bml)
 
  end subroutine prg_pulaycomponent0
 
  subroutine prg_pulaycomponentt(rho_bml,ham_bml,zmat_bml,pcm_bml,threshold &
       &,M,bml_type,verbose)
 
    implicit none
 
    type(bml_matrix_t), intent(in) :: rho_bml
    type(bml_matrix_t), intent(in) :: ham_bml
    type(bml_matrix_t), intent(in) :: zmat_bml
    type(bml_matrix_t) :: aux_bml
    type(bml_matrix_t) :: aux1_bml
    type(bml_matrix_t), intent(inout) :: pcm_bml
    integer, intent(in) :: m
    integer ::  norb, verbose
    real(dp), intent(in) :: threshold
    character(20), intent(in) :: bml_type
 
    if(verbose.eq.2) write(*,*)"In prg_PulayComponentT ..."
 
    norb = bml_get_n(rho_bml)
 
    if(bml_get_n(pcm_bml).le.0)then !If pcm is not allocated.
      call bml_zero_matrix(bml_type,bml_element_real,dp,norb,norb,pcm_bml)
    else
      call bml_deallocate(pcm_bml) !If pcm is allocated then we set it to 0.
      call bml_zero_matrix(bml_type,bml_element_real,dp,norb,norb,pcm_bml)
    endif
 
    call bml_zero_matrix(bml_type,bml_element_real,dp,norb ,norb,aux_bml)
    call bml_zero_matrix(bml_type,bml_element_real,dp,norb ,norb,aux1_bml)
 
    if(bml_get_n(zmat_bml).lt.0)then
      stop 'ERROR: zmat_bml not allocated'
    endif
 
    call bml_multiply(zmat_bml,zmat_bml,aux_bml,1.0d0,1.0d0,threshold)  !aux=Z*Z
 
    call bml_multiply(aux_bml,ham_bml,pcm_bml,1.0d0,1.0d0,threshold)  !Z*Z*H
 
    call bml_multiply(pcm_bml,rho_bml,aux1_bml,1.0d0,1.0d0,threshold) !aux1=Z*Z*H*D
 
    call bml_scale(0.0_dp, pcm_bml)
 
    call bml_multiply(ham_bml,aux_bml,pcm_bml,1.0d0,1.0d0,threshold)  !H*Z*Z
 
    call bml_scale(0.0_dp, aux_bml)  !We set aux to 0 to recicle the variable
 
    call bml_multiply(rho_bml,pcm_bml,aux_bml,1.0d0,1.0d0,threshold) !D*H*Z*Z
 
    call bml_add_deprecated(0.5d0,aux_bml,0.5d0,aux1_bml)
 
    call bml_copy(aux_bml,pcm_bml)
 
    call bml_deallocate(aux_bml)
    call bml_deallocate(aux1_bml)
 
 
  end subroutine prg_pulaycomponentt
 
 
  subroutine prg_get_pulayforce(nats,zmat_bml,ham_bml,rho_bml,&
       dSx_bml,dSy_bml,dSz_bml,hindex,FPUL,threshold)
    implicit none
    real(dp), allocatable, intent(inout) :: fpul(:,:)
    integer, intent(in) :: nats
    type(bml_matrix_t), intent(in)  ::  dsx_bml, dsy_bml, dsz_bml
    type(bml_matrix_t), intent(in)  ::  rho_bml, ham_bml, zmat_bml
    integer, intent(in)                ::  hindex(:,:)
    integer :: i_a, i_b, i , j, norb
    real(dp), intent(in) :: threshold
    type(bml_matrix_t)  :: xtmp_bml, ytmp_bml, ztmp_bml, sihd_bml
    type(bml_matrix_t)  :: aux_bml
    real(dp), allocatable :: diagxtmp(:), diagytmp(:), diagztmp(:)
    real(dp) :: partrace
 
    write(*,*)"In prg_get_pulayforce ..."
 
    if(.not.allocated(fpul))then
      allocate(fpul(3,nats))
    endif
 
    fpul = 0.0_dp
 
    norb = bml_get_n(rho_bml)
 
    !SIHD = 2*Z*Z'*H*D;
    call bml_copy_new(zmat_bml,sihd_bml)
    call bml_copy_new(zmat_bml,aux_bml)
    call bml_transpose(zmat_bml,aux_bml) !Z'
    call bml_multiply(zmat_bml,aux_bml,sihd_bml,2.0_dp,0.0_dp,threshold) !2*Z*Z'
    call bml_multiply(sihd_bml,ham_bml,aux_bml,1.0_dp,0.0_dp,threshold) !2*Z*Z'*H
    call bml_multiply(aux_bml,rho_bml,sihd_bml,1.0_dp,0.0_dp,threshold) !2*Z*Z'*D
    call bml_deallocate(aux_bml)
 
    call bml_copy_new(rho_bml,xtmp_bml)
    allocate(diagxtmp(norb))
    call bml_multiply(dsx_bml,sihd_bml,xtmp_bml,1.0_dp,0.0_dp,threshold)
    call bml_get_diagonal(xtmp_bml,diagxtmp)
    call bml_deallocate(xtmp_bml)
 
    call bml_copy_new(rho_bml,ytmp_bml)
    allocate(diagytmp(norb))
    call bml_multiply(dsy_bml,sihd_bml,ytmp_bml,1.0_dp,0.0_dp,threshold)
    call bml_get_diagonal(ytmp_bml,diagytmp)
    call bml_deallocate(ytmp_bml)
 
    call bml_copy_new(rho_bml,ztmp_bml)
    allocate(diagztmp(norb))
    call bml_multiply(dsz_bml,sihd_bml,ztmp_bml,1.0_dp,0.0_dp,threshold)
    call bml_get_diagonal(ztmp_bml,diagztmp)
    call bml_deallocate(ztmp_bml)
 
    call bml_deallocate(sihd_bml)
 
    !$omp parallel do default(none) private(i) &
    !$omp private(I_A,I_B,j,partrace) &
    !$omp shared(hindex,diagxtmp,diagytmp,diagztmp,FPUL,nats)
    do i = 1,nats
      i_a = hindex(1,i);
      i_b = hindex(2,i);
 
      partrace = 0.0_dp
      do j=i_a,i_b
        partrace = partrace + diagxtmp(j)
      enddo
      fpul(1,i) = partrace;
 
      partrace = 0.0_dp
      do j=i_a,i_b
        partrace = partrace + diagytmp(j)
      enddo
      fpul(2,i) = partrace;
 
      partrace = 0.0_dp
      do j=i_a,i_b
        partrace = partrace + diagztmp(j)
      enddo
      fpul(3,i) = partrace;
 
    enddo
    !$omp end parallel do
 
    deallocate(diagxtmp)
    deallocate(diagytmp)
    deallocate(diagztmp)
 
  end subroutine prg_get_pulayforce
 
end module prg_pulaycomponent_mod