corresponding_orbital_projection Subroutine

public subroutine corresponding_orbital_projection(vb, sba, va, ndoc, nact, nproj, nbf, l1co, l0)

Uses

  • proc~~corresponding_orbital_projection~~UsesGraph proc~corresponding_orbital_projection corresponding_orbital_projection module~eigen eigen proc~corresponding_orbital_projection->module~eigen module~mathlib_types mathlib_types module~eigen->module~mathlib_types module~messages messages module~eigen->module~messages module~oqp_linalg oqp_linalg module~eigen->module~oqp_linalg module~precision precision module~eigen->module~precision module~messages->module~precision comm_IOFILE comm_IOFILE module~messages->comm_IOFILE comm_PAR comm_PAR module~messages->comm_PAR module~io_constants io_constants module~messages->module~io_constants module~blas_wrap blas_wrap module~oqp_linalg->module~blas_wrap module~lapack_wrap lapack_wrap module~oqp_linalg->module~lapack_wrap iso_fortran_env iso_fortran_env module~precision->iso_fortran_env module~blas_wrap->module~mathlib_types module~blas_wrap->module~messages module~blas_wrap->module~precision module~lapack_wrap->module~mathlib_types module~lapack_wrap->module~messages module~lapack_wrap->module~precision

@brief Corresponding orbital projection @param[in] nproj number of orbitals from vb to projected onto @param[in] l0 number of the first orbitals from the va space to project vb to @param[in,out] va a orbitals on entry, replaced by corresponding orbitals on exit @param[in] vb b orbitals on entry, unchanged on exit @param[in] sba overlap integrals between the two bases,

Note

The corresponding orbitals in the b basis are not generated.

Note

H.F.King, R.E.Stanton, H.Kim, R.E.Wyatt, R.G.Parr J.Chem.Phys. 47, 1936-1941 (1967)

Arguments

Type IntentOptional Attributes Name
real(kind=dp) :: vb(l1co,nproj)
real(kind=dp) :: sba(l1co,nbf)
real(kind=dp) :: va(nbf,l0)
integer :: ndoc
integer :: nact
integer :: nproj
integer :: nbf
integer :: l1co
integer :: l0

Calls

proc~~corresponding_orbital_projection~~CallsGraph proc~corresponding_orbital_projection corresponding_orbital_projection interface~show_message show_message proc~corresponding_orbital_projection->interface~show_message proc~diag_symm_full diag_symm_full proc~corresponding_orbital_projection->proc~diag_symm_full proc~oqp_dgemm_i64 oqp_dgemm_i64 proc~corresponding_orbital_projection->proc~oqp_dgemm_i64 proc~oqp_dswap_i64 oqp_dswap_i64 proc~corresponding_orbital_projection->proc~oqp_dswap_i64 proc~oqp_dsyrk_i64 oqp_dsyrk_i64 proc~corresponding_orbital_projection->proc~oqp_dsyrk_i64 proc~schmd schmd proc~corresponding_orbital_projection->proc~schmd proc~diag_symm_full->interface~show_message dsyev dsyev proc~diag_symm_full->dsyev proc~oqp_dgemm_i64->interface~show_message dgemm dgemm proc~oqp_dgemm_i64->dgemm proc~oqp_dswap_i64->interface~show_message dswap dswap proc~oqp_dswap_i64->dswap proc~oqp_dsyrk_i64->interface~show_message dsyrk dsyrk proc~oqp_dsyrk_i64->dsyrk proc~schmd->interface~show_message proc~oqp_dgeqrf_i64 oqp_dgeqrf_i64 proc~schmd->proc~oqp_dgeqrf_i64 proc~oqp_dorgqr_i64 oqp_dorgqr_i64 proc~schmd->proc~oqp_dorgqr_i64 proc~oqp_dgeqrf_i64->interface~show_message dgeqrf dgeqrf proc~oqp_dgeqrf_i64->dgeqrf proc~oqp_dorgqr_i64->interface~show_message dorgqr dorgqr proc~oqp_dorgqr_i64->dorgqr

Called by

proc~~corresponding_orbital_projection~~CalledByGraph proc~corresponding_orbital_projection corresponding_orbital_projection proc~huckel_guess huckel_guess proc~huckel_guess->proc~corresponding_orbital_projection proc~guess_huckel guess_huckel proc~guess_huckel->proc~huckel_guess proc~guess_huckel_c guess_huckel_C proc~guess_huckel_c->proc~guess_huckel

Source Code

 subroutine corresponding_orbital_projection(vb, sba, va, &
                   ndoc, nact, nproj, nbf, l1co, l0)

  use eigen, only: schmd, diag_symm_full

  implicit none

  real(kind=dp) :: vb(l1co,nproj), sba(l1co,nbf), va(nbf,l0)
  integer :: ndoc, nact, nproj, nbf, l1co, l0

  real(kind=dp), allocatable :: &
        vaco(:,:), rhs(:,:), d(:,:), tmp(:,:), scr(:)
  integer, allocatable :: iwrk(:)
  logical, allocatable :: unused(:)
  integer :: ok

  allocate(vaco(nbf,nbf), &
           rhs(nbf,nbf), &
           d(nbf,l0), &
           tmp(nbf,nbf), &
           scr(nbf), &
           source=0.0_dp, &
           stat=ok)
  allocate(iwrk(nbf), &
           source=0, &
           stat=ok)
  allocate(unused(nbf), &
           source=.false., &
           stat=ok)

! Vb^T * Sba
  call dgemm('t', 'n', nproj, nbf, l1co, &
             1.0_dp, vb, l1co, &
                     sba, l1co, &
             0.0_dp, rhs, nbf)

! 1. Project core space
  call project(va, rhs, 1, ndoc, l0, nbf)
! 2. Project active space
  call project(va, rhs, ndoc+1, ndoc+nact, l0, nbf)
! 3. Project virtual space
  call project(va, rhs, ndoc+nact+1, nproj, l0, nbf)

 contains

  subroutine project(va, rhs, minmo, maxmo, l0, nbf)

    implicit none

    real(kind=dp) :: va(nbf,*), rhs(nbf,*)
    integer :: minmo, maxmo, l0, nbf

    integer :: i, j, ok
    integer :: nrest, nmo

    nrest = l0-minmo+1
    nmo  = maxmo-minmo+1
    if (nrest <= 0) return
    if (nmo == 0) return

!   D = Vb^T * Sab * Va
    call dgemm('n','n', nmo, nrest, nbf, &
               1.0_dp, rhs(minmo,1), nbf, &
                       va(1,minmo),nbf, &
               0.0_dp, d, nbf)

!   Diagonalize D^T * D
!   Eigenvalues run 0 to 1, we want the highest ones first,
!   that is why we use factor -1
!   We don't need eigenvalues
    call dsyrk('u', 't', nrest, nmo, &
               -1.0_dp, d, nbf, &
                0.0_dp, vaco, nbf)
    call diag_symm_full(1, nrest, vaco, nbf, scr, ok)

!   The near zero overlap part of the space gets scrambled
!   in the above diagonalization, we can get a symmetry
!   adapted space by Gram-Schmidt instead.
    call schmd(vaco,nmo,nrest,nbf,scr)

!   Rotate va to corresponding orbital set, a'=a*v,
    call dgemm('n', 'n', nbf, nrest, nrest, &
               1.0_dp, va(1,minmo), nbf, &
                       vaco, nbf, &
               0.0_dp, tmp, nbf)

!   Copy projected orbitals back to end of va
    va(:,minmo:l0) = tmp(:,1:nrest)

!   Get overlap between `b` and the `a'` corresponding set,
!   D = Vb^T * Sab * VaCO
    call dgemm('t', 't', nrest, nmo, nbf, &
               1.0_dp, va(1,minmo), nbf, &
                       rhs(minmo,1), nbf, &
               0.0_dp, tmp, nbf)

!   Permute `a'` set to maximum overlap with `b` set
    unused(:nmo) = .true.
    do i = 1, nmo
      j = maxloc(abs(tmp(:nmo,i)), dim=1, mask=unused)
      iwrk(i) = j
      unused(j) = .false.
    end do

    call reorder_columns(va(1,minmo),iwrk,nmo,nbf)

  end subroutine

 end subroutine corresponding_orbital_projection