int1e.F90 Source File

module int1e_mod

  implicit none

  character(len=*), parameter :: module_name = "int1e_mod"

  private
  public int1e

contains

  subroutine int1e_C(c_handle) bind(C, name="int1e")
    use c_interop, only: oqp_handle_t, oqp_handle_get_info
    use types, only: information
    type(oqp_handle_t) :: c_handle
    type(information), pointer :: inf
    inf => oqp_handle_get_info(c_handle)
    call int1e(inf)
  end subroutine int1e_C

!> @brief Calculate the basic H, S, and T 1e-integrals
  subroutine int1e(infos)

    use types, only: information
    use oqp_tagarray_driver
    use precision, only: dp
    use io_constants, only: iw
    use constants, only: tol_int
    use int1, only: omp_hst
    use basis_tools, only: basis_set
    use printing, only: print_sym_labeled
    use messages, only: show_message, WITH_ABORT
    use strings, only: Cstring, fstring
    use physical_constants, only: BOHR_TO_ANGSTROM
    use printing, only: print_module_info
    use qmmm_mod, only: oqp_esp_qmmm
    use dk_scalar_mod, only: dk_scalar

    implicit none

    character(len=*), parameter :: subroutine_name = "int1e"

    type(information), target, intent(inout) :: infos
    type(basis_set), pointer :: basis

    real(kind=dp) :: tol
    integer :: i, nbf, nat, nbf2, dk

    ! tagarray
    real(kind=dp), contiguous, pointer :: &
      hcore(:), tmat(:), smat(:), Hqmmm(:), mm_potential(:)
    character(len=*), parameter :: tags_general(3) = (/ character(len=80) :: &
      OQP_SM, OQP_TM, OQP_Hcore /)
    character(len=*), parameter :: tags_stale(1) = (/ character(len=80) :: &
      OQP_QMAT /)
    character(len=*), parameter :: tags_qmmm(2) = (/ character(len=80) :: &
      OQP_mm_potential,  OQP_Hqmmm/)

    logical dbg
    dbg = .false.

    dk = infos%control%scal_rel
!   Files open:
!   LOG: Read and Write: Main output file
    open(unit=iw,  file=infos%log_filename, position="append")

!   Load basis set
    basis => infos%basis
    basis%atoms => infos%atoms
!
    call print_module_info('int1e','Computing H, S and T Matrices')

!   Print out the Cartesian Coordinates.
    write(iw, fmt="(&
              &/21X,32('=')&
              &/21X,a&
              &/21X,32('=')&
              &/8X,'ATOM     ZNUC',11X,'X',14X,'Y',14X,'Z'&
              &/6X,62('-'))") "Cartesian Coordinate in Angstrom"

    do i = 1, size(basis%atoms%zn(:))
       write(iw,'(7x,i4,5x,f4.1,3(x,f15.9))') &
               i, basis%atoms%zn(i), basis%atoms%xyz(1:3,i)*BOHR_TO_ANGSTROM
    end do

!   Allocate H, S and T matrices
    nbf2 = basis%nbf*(basis%nbf+1)/2
    nat = ubound(infos%atoms%zn,1)

!   The overlap matrix changes, so the cached Q = S^(-1/2) is stale
    call infos%dat%erase(tags_stale)

!   Allocate H, S and T and bind typed pointers (one call each). alloc_or_die
!   allocates or reuses in place, so the Python-set OQP_mm_potential tagarray is
!   preserved (only the stale SM/TM/Hcore/Q records are refreshed here).
    call infos%dat%alloc_or_die(OQP_SM,    (/ nbf2 /), smat,  description=OQP_SM_comment)
    call infos%dat%alloc_or_die(OQP_TM,    (/ nbf2 /), tmat,  description=OQP_TM_comment)
    call infos%dat%alloc_or_die(OQP_Hcore, (/ nbf2 /), hcore, description=OQP_Hcore_comment)

!   Create arrays of atomic coordinates and charges for one-electron code
    nbf = basis%nbf
    nat = ubound(infos%atoms%zn,1)

!   Compute conventional H, S, and T integrals
    tol = log(10.0d0)*tol_int
    call omp_hst(basis, infos%atoms%xyz, infos%atoms%zn - infos%basis%ecp_zn_num, hcore, smat, tmat,&
            logtol=tol, comm=infos%mpiinfo%comm, usempi=infos%mpiinfo%usempi)

!   Compute QM/MM interaction
!    if(infos%control%qmmm_flag) then
!       call infos%dat%reserve_data(OQP_Hqmmm, TA_TYPE_REAL64, nbf2, comment=OQP_Hqmmm_comment)
!       call data_has_tags(infos%dat, tags_qmmm, module_name, subroutine_name, WITH_ABORT)
!       call tagarray_get_data(infos%dat, OQP_Hqmmm, Hqmmm)
!       call tagarray_get_data(infos%dat, OQP_mm_potential, mm_potential)
!
!       write(iw,"(/1X,'  Computing ESP One Electron Integrals (QM/MM) '/)")
!       write(iw,"('External MM potential:'/)")
!       do i=1,nat
!          write(iw,"(i4,1X,f12.8)") i, mm_potential(i)
!       end do
!!   Compute QM/MM contribution to core Hamiltonian
!!       call oqp_esp_qmmm(infos, Hqmmm, mm_potential, smat, logtol=tol)
!!   Add QM/MM contribution to core Hamiltonian
!       hcore = hcore + hqmmm
!       write(iw,"(/1X,'  ... End of ESP One Electron Integrals ... '/)")
!    endif

!   Douglas-Kroll scalar-relativistic correction to the core Hamiltonian (SOC)
    if (dk.gt.0) call dk_scalar(infos)

    if (dbg) then
        if(infos%control%qmmm_flag) then
            write(iw,'(/"BARE NUCLEUS HAMILTONIAN INTEGRALS (H=T+V+QM/MM)")')
        else
            write(iw,'(/"BARE NUCLEUS HAMILTONIAN INTEGRALS (H=T+V)")')
        end if
        call print_sym_labeled(hcore,nbf, basis)

        write(iw,'(/"OVERLAP MATRIX")')
        call print_sym_labeled(Smat,nbf, basis)

        write(iw,'(/"KINETIC ENERGY INTEGRALS")')
        call print_sym_labeled(tmat,nbf, basis)

        if (infos%control%qmmm_flag) then
           write(iw,'(/"QM/MM HAMILTONIAN INTEGRALS")')
           call print_sym_labeled(Hqmmm,nbf, basis)
        end if
    end if

    write(iw,"(/1X,'...... End Of One Electron Integrals ......'/)")

    close(iw)

  end subroutine int1e

end module int1e_mod