prg_pulaycomponent_mod.F90

Go to the documentation of this file.

 
module prg_pulaycomponent_mod
 
  use bml
  use iso_c_binding
 
#ifdef USE_NVTX
  use prg_nvtx_mod
#endif
  
  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,&
       &verbose)
    !&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
    character(20) ::  bml_type
 
    if(verbose.eq.2) write(*,*)"In prg_PulayComponent0 ..."
 
    norb = bml_get_n(rho_bml)
    bml_type = bml_get_deep_type(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
      ! re-allocate will break the c pointer to pcm (for the c wrapper)
      !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)
      call bml_scale(0.0_dp, pcm_bml)
    endif
 
    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,partracex,partracey,partracez
#ifdef USE_OFFLOAD
    type(c_ptr)                        :: xtmp_bml_c_ptr, ytmp_bml_c_ptr, ztmp_bml_c_ptr
    integer :: ld
    real(c_double), pointer            :: xtmp_bml_ptr(:,:), ytmp_bml_ptr(:,:), ztmp_bml_ptr(:,:)
#endif
 
    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(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)
    call bml_multiply(dsx_bml,sihd_bml,xtmp_bml,1.0_dp,0.0_dp,threshold)
    call bml_copy_new(rho_bml,ytmp_bml)
    call bml_multiply(dsy_bml,sihd_bml,ytmp_bml,1.0_dp,0.0_dp,threshold)
    call bml_copy_new(rho_bml,ztmp_bml)
    call bml_multiply(dsz_bml,sihd_bml,ztmp_bml,1.0_dp,0.0_dp,threshold)
 
#ifdef USE_OFFLOAD
    xtmp_bml_c_ptr = bml_get_data_ptr_dense(xtmp_bml)
    ytmp_bml_c_ptr = bml_get_data_ptr_dense(ytmp_bml)
    ztmp_bml_c_ptr = bml_get_data_ptr_dense(ztmp_bml)
    ld = bml_get_ld_dense(xtmp_bml)
    call c_f_pointer(xtmp_bml_c_ptr,xtmp_bml_ptr,shape=[ld,norb])
    call c_f_pointer(ytmp_bml_c_ptr,ytmp_bml_ptr,shape=[ld,norb])
    call c_f_pointer(ztmp_bml_c_ptr,ztmp_bml_ptr,shape=[ld,norb])
    !$acc enter data copyin(FPUL(1:3,1:nats),hindex(1:2,1:nats))
    !$acc parallel loop deviceptr(Xtmp_bml_ptr,Ytmp_bml_ptr,Ztmp_bml_ptr) &
    !$acc private(i,j,partracex,partracey,partracez)
    do i = 1,nats
      partracex = 0.0_dp
      partracey = 0.0_dp
      partracez = 0.0_dp
      !$acc loop reduction(partracex,partracey,partracez)
      do j=hindex(1,i),hindex(2,i)
        partracex = partracex + xtmp_bml_ptr(j,j)
        partracey = partracey + ytmp_bml_ptr(j,j)
        partracez = partracez + ztmp_bml_ptr(j,j)
     enddo
     !$acc end loop
      fpul(1,i) = partracex;
      fpul(2,i) = partracey;
      fpul(3,i) = partracez;
 
   enddo
   !$acc end loop
   !$acc exit data copyout(FPUL(1:3,1:nats)) delete(hindex(1:2,1:nats))
#else
    allocate(diagxtmp(norb))
    call bml_get_diagonal(xtmp_bml,diagxtmp)
 
    allocate(diagytmp(norb))
    call bml_get_diagonal(ytmp_bml,diagytmp)
 
    allocate(diagztmp(norb))
    call bml_get_diagonal(ztmp_bml,diagztmp)
 
    !$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)
#endif
    call bml_deallocate(xtmp_bml)
    call bml_deallocate(ytmp_bml)
    call bml_deallocate(ztmp_bml)
    call bml_deallocate(sihd_bml)
 
  end subroutine prg_get_pulayforce
 
end module prg_pulaycomponent_mod