electric_moments.F90 Source File

module electric_moments_mod

  implicit none

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

  private
  public electric_moments
  public electric_moments_excited
  public electric_dipole_au_C

contains

  subroutine electric_dipole_au_C(c_handle, dipole) bind(C, name="electric_dipole_au")
    use c_interop, only: oqp_handle_t, oqp_handle_get_info
    use iso_c_binding, only: c_double
    use types, only: information
    type(oqp_handle_t) :: c_handle
    real(c_double), intent(out) :: dipole(3)
    type(information), pointer :: inf
    inf => oqp_handle_get_info(c_handle)
    call electric_dipole_au(inf, dipole)
  end subroutine electric_dipole_au_C

  subroutine electric_dipole_au(infos, dipole)
    use oqp_tagarray_driver
    use basis_tools, only: basis_set
    use messages, only: show_message, with_abort
    use types, only: information
    use mathlib, only: traceprod_sym_packed
    use int1, only: multipole_integrals
    type(information), target, intent(inout) :: infos
    real(kind=8), intent(out) :: dipole(3)

    integer :: nbf, nbf2, ok, nat, i
    logical :: urohf
    type(basis_set), pointer :: basis
    real(kind=8), allocatable :: mints(:,:)
    real(kind=8) :: origin(3), dr(3), z
    real(kind=8), contiguous, pointer :: dmat_a(:), dmat_b(:)
    integer(4) :: status

    urohf = infos%control%scftype == 2 .or. infos%control%scftype == 3
    basis => infos%basis
    basis%atoms => infos%atoms
    nbf = basis%nbf
    nbf2 = nbf*(nbf+1)/2
    nat = ubound(basis%atoms%zn,1)
    allocate(mints(nbf2,19), source=0.0d0, stat=ok)
    if (ok /= 0) call show_message('Cannot allocate memory', WITH_ABORT)

    call tagarray_get_data(infos%dat, OQP_DM_A, dmat_a, status)
    call check_status(status, module_name, "electric_dipole_au", OQP_DM_A)
    if (urohf) then
      call tagarray_get_data(infos%dat, OQP_DM_B, dmat_b, status)
      call check_status(status, module_name, "electric_dipole_au", OQP_DM_B)
    endif

    origin = 0.0d0
    call multipole_integrals(basis, mints, origin, 3)

    dipole = 0.0d0
    do i = 1, nat
      dr = infos%atoms%xyz(:,i)
      z = infos%atoms%zn(i) - infos%basis%ecp_zn_num(i)
      dipole = dipole + z * dr
    end do
    do i = 1, 3
      dipole(i) = dipole(i) - traceprod_sym_packed(mints(:,i), dmat_a, nbf)
    end do
    if (urohf) then
      do i = 1, 3
        dipole(i) = dipole(i) - traceprod_sym_packed(mints(:,i), dmat_b, nbf)
      end do
    end if

    deallocate(mints)
  end subroutine electric_dipole_au

  subroutine electric_moments_C(c_handle) bind(C, name="electric_moments")
    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 electric_moments(inf)
  end subroutine electric_moments_c

  subroutine electric_moments_excited_C(c_handle) bind(C, name="electric_moments_excited")
    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 electric_moments_excited(inf)

  end subroutine electric_moments_excited_C


  subroutine electric_moments(infos)
    use io_constants, only: iw
    use oqp_tagarray_driver
    use basis_tools, only: basis_set
    use messages, only: show_message, with_abort
    use types, only: information
    use strings, only: Cstring, fstring
    use mathlib, only: traceprod_sym_packed, triangular_to_full
    use int1, only: multipole_integrals
    use physical_constants, only: AU_TO_DEBYE, AU_TO_BUCK, AU_TO_OCT
    use xyz_order

    implicit none

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

    type(information), target, intent(inout) :: infos

    integer :: nbf, nbf2, ok
    logical :: urohf
    type(basis_set), pointer :: basis
    real(kind=8), allocatable :: mints(:,:)
    real(kind=8) :: com(3), dip(3), qxyz(6), quad(3,3), dr(3), z
    real(kind=8) :: oxyz(10), oct(10)
    integer :: nat, i

    real(kind=8), contiguous, pointer :: dmat_a(:), dmat_b(:)
    integer(4) :: status

    urohf = infos%control%scftype == 2 .or. infos%control%scftype == 3

!   3. LOG: Write: Main output file
    open (unit=IW, file=infos%log_filename, position="append")

!   Load basis set
    basis => infos%basis
    basis%atoms => infos%atoms

!   Allocate memory
    nbf = basis%nbf
    nbf2 = nbf*(nbf+1)/2

    write(iw,'(2/)')
    write(iw,'(4x,a)') '========================'
    write(iw,'(4x,a)') 'Electric moment analysis'
    write(iw,'(4x,a)') '========================'
    call flush(iw)

    nat = ubound(basis%atoms%zn,1)
    allocate(mints(nbf2,19), source=0.0d0, stat=ok)
    if (ok /= 0) call show_message('Cannot allocate memory', WITH_ABORT)

    call tagarray_get_data(infos%dat, OQP_DM_A, dmat_a, status)
    call check_status(status, module_name, subroutine_name, OQP_DM_A)

    if (urohf) then
      call tagarray_get_data(infos%dat, OQP_DM_B, dmat_b, status)
      call check_status(status, module_name, subroutine_name, OQP_DM_B)
    endif

    com = 0
    do i = 1, nat
        com = com + basis%atoms%xyz(:,i)*basis%atoms%mass(i)
    end do
    com = com / sum(basis%atoms%mass)

    call multipole_integrals(basis, mints, com, 3)

    dip = 0
    qxyz = 0
    oxyz = 0

    do i = 1, nat
      dr = infos%atoms%xyz(:,i)-com
      z = infos%atoms%zn(i) - infos%basis%ecp_zn_num(i)

      dip = dip + z * dr

      qxyz(XX_) = qxyz(XX_) + z * dr(X__)*dr(X__)
      qxyz(YY_) = qxyz(YY_) + z * dr(Y__)*dr(Y__)
      qxyz(ZZ_) = qxyz(ZZ_) + z * dr(Z__)*dr(Z__)
      qxyz(XY_) = qxyz(XY_) + z * dr(X__)*dr(Y__)
      qxyz(XZ_) = qxyz(XZ_) + z * dr(X__)*dr(Z__)
      qxyz(YZ_) = qxyz(YZ_) + z * dr(Y__)*dr(Z__)

      oxyz(XXX) = oxyz(XXX) + z * dr(X__)*dr(X__)*dr(X__)
      oxyz(YYY) = oxyz(YYY) + z * dr(Y__)*dr(Y__)*dr(Y__)
      oxyz(ZZZ) = oxyz(ZZZ) + z * dr(Z__)*dr(Z__)*dr(Z__)
      oxyz(XXY) = oxyz(XXY) + z * dr(X__)*dr(X__)*dr(Y__)
      oxyz(XXZ) = oxyz(XXZ) + z * dr(X__)*dr(X__)*dr(Z__)
      oxyz(YYX) = oxyz(YYX) + z * dr(Y__)*dr(Y__)*dr(X__)
      oxyz(YYZ) = oxyz(YYZ) + z * dr(Y__)*dr(Y__)*dr(Z__)
      oxyz(ZZX) = oxyz(ZZX) + z * dr(Z__)*dr(Z__)*dr(X__)
      oxyz(ZZY) = oxyz(ZZY) + z * dr(Z__)*dr(Z__)*dr(Y__)
      oxyz(XYZ) = oxyz(XYZ) + z * dr(X__)*dr(Y__)*dr(Z__)
    end do

    do i = 1, 3
        dip(i) = dip(i) - traceprod_sym_packed(mints(:,i), dmat_a, nbf)
    end do

    do i = 1, 6
        qxyz(i) = qxyz(i) - traceprod_sym_packed(mints(:,3+i), dmat_a, nbf)
    end do

    do i = 1, 10
        oxyz(i) = oxyz(i) - traceprod_sym_packed(mints(:,9+i), dmat_a, nbf)
    end do

    if (urohf) then
      do i = 1, 3
        dip(i) = dip(i) - traceprod_sym_packed(mints(:,i), dmat_b, nbf)
      end do
      do i = 1, 6
        qxyz(i) = qxyz(i) - traceprod_sym_packed(mints(:,3+i), dmat_b, nbf)
      end do
    do i = 1, 10
        oxyz(i) = oxyz(i) - traceprod_sym_packed(mints(:,9+i), dmat_b, nbf)
    end do
    end if

    dip = dip*AU_TO_DEBYE

    ! Assemble quadrupole tensor:
    ! Q_ab = 1/2 * ( 3 * r_a * r_b - r^2 * \delta(a,b) )
    quad(X__,X__) = 0.5*(2*qxyz(XX_) - (qxyz(YY_)+qxyz(ZZ_)))
    quad(X__,Y__) = 0.5*(3*qxyz(XY_)                        )
    quad(X__,Z__) = 0.5*(3*qxyz(XZ_)                        )
    quad(Y__,Y__) = 0.5*(2*qxyz(YY_) - (qxyz(XX_)+qxyz(ZZ_)))
    quad(Y__,Z__) = 0.5*(3*qxyz(YZ_)                        )
    quad(Z__,Z__) = 0.5*(2*qxyz(ZZ_) - (qxyz(XX_)+qxyz(YY_)))
    call triangular_to_full(quad,3,'u')
    quad = quad*AU_TO_BUCK

    ! Assemble octopole tensor:
    ! O_abc = 1/2*(5*r_a*r_b*r_c - r^2*(r_a*\delta(b,c)+r_b*\delta(a,c)+r_c*\delta(a,b))
    oct(XXX) = 0.5*(2*oxyz(XXX) - 3*(oxyz(YYX)+oxyz(ZZX)))
    oct(YYY) = 0.5*(2*oxyz(YYY) - 3*(oxyz(XXY)+oxyz(ZZY)))
    oct(ZZZ) = 0.5*(2*oxyz(ZZZ) - 3*(oxyz(YYZ)+oxyz(XXZ)))
    oct(XXY) = 0.5*(4*oxyz(XXY) -   (oxyz(YYY)+oxyz(ZZY)))
    oct(XXZ) = 0.5*(4*oxyz(XXZ) -   (oxyz(YYZ)+oxyz(ZZZ)))
    oct(YYX) = 0.5*(4*oxyz(YYX) -   (oxyz(XXX)+oxyz(ZZX)))
    oct(YYZ) = 0.5*(4*oxyz(YYZ) -   (oxyz(XXZ)+oxyz(ZZZ)))
    oct(ZZX) = 0.5*(4*oxyz(ZZX) -   (oxyz(XXX)+oxyz(YYX)))
    oct(ZZY) = 0.5*(4*oxyz(ZZY) -   (oxyz(XXY)+oxyz(YYY)))
    oct(XYZ) = 0.5*(5*oxyz(XYZ)                          )

    oct = oct*AU_TO_OCT

    write(iw,'(/1x,a)') 'At point (Bohr):'
    write(iw,'(4x,40x,3(a8,7x))') 'X', 'Y', 'Z'
    write(iw,'(4x,40x,3f15.8,sp,f12.5)') com

    write(iw,'(/4x,a)') 'electric charge (a.u.):'
    write(iw,'(4x,34x,a6,sp,f12.5)') 'C_o', real(infos%mol_prop%charge)

    write(iw,'(/4x,a)') 'electric dipole (Debye):'
    write(iw,'(4x,40x,3(a8,7x),a10)') 'X', 'Y', 'Z', 'Norm'
    write(iw,'(4x,34x,a6,4f15.8)') 'D_o', dip, norm2(dip)

    write(iw,'(/4x,a)') 'electric quadrupole (Buckingham):'
    write(iw,'(4x,40x,3(a8,7x))') 'X', 'Y', 'Z'
    write(iw,'(4x,34x,a6,3F15.8)') 'Q_X', quad(:,1)
    write(iw,'(4x,34x,a6,3f15.8)') 'Q_Y', quad(:,2)
    write(iw,'(4x,34x,a6,3f15.8)') 'Q_Z', quad(:,3)

    write(iw,'(/4x,a)') 'electric octupole (Buckingham*Angstrom):'
    write(iw,'(4x,34x,a6,*(F15.8))') 'O_xxx', oct(XXX)
    write(iw,'(4x,34x,a6,*(F15.8))') 'O_xxy', oct(XXY)
    write(iw,'(4x,34x,a6,*(F15.8))') 'O_xxz', oct(XXZ)
    write(iw,'(4x,34x,a6,*(F15.8))') 'O_xyy', oct(YYX)
    write(iw,'(4x,34x,a6,*(F15.8))') 'O_xyz', oct(XYZ)
    write(iw,'(4x,34x,a6,*(F15.8))') 'O_xzz', oct(ZZX)
    write(iw,'(4x,34x,a6,*(F15.8))') 'O_yyy', oct(YYY)
    write(iw,'(4x,34x,a6,*(F15.8))') 'O_yyz', oct(YYZ)
    write(iw,'(4x,34x,a6,*(F15.8))') 'O_zzy', oct(ZZY)
    write(iw,'(4x,34x,a6,*(F15.8))') 'O_zzz', oct(ZZZ)

    close(iw)

  end subroutine electric_moments

  subroutine electric_moments_excited(infos)
    use io_constants, only: iw
    use oqp_tagarray_driver
    use basis_tools, only: basis_set
    use messages, only: show_message, with_abort
    use types, only: information
    use strings, only: Cstring, fstring
    use mathlib, only: traceprod_sym_packed, triangular_to_full
    use int1, only: multipole_integrals
    use physical_constants, only: AU_TO_DEBYE, AU_TO_BUCK, AU_TO_OCT
    use xyz_order

    implicit none

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

    type(information), target, intent(inout) :: infos

    integer :: nbf, nbf2, ok
    logical :: urohf
    type(basis_set), pointer :: basis
    real(kind=8), allocatable :: mints(:,:)
    real(kind=8), allocatable :: dens_ex(:)
    real(kind=8) :: com(3), dip(3), qxyz(6), quad(3,3), dr(3), z
    real(kind=8) :: oxyz(10), oct(10)
    integer :: nat, i, istate

    real(kind=8), contiguous, pointer :: dmat_a(:), dmat_b(:)
    real(kind=8), contiguous, pointer :: td_p(:,:)
    integer(4) :: status

    urohf = infos%control%scftype == 2 .or. infos%control%scftype == 3

    open (unit=IW, file=infos%log_filename, position="append")

    basis => infos%basis
    basis%atoms => infos%atoms

    nbf  = basis%nbf
    nbf2 = nbf*(nbf+1)/2

    istate = infos%tddft%target_state

    write(iw,'(2/)')
    write(iw,'(4x,a)')      '================================================='
    write(iw,'(4x,a,i4)')   'Electric moment analysis for excited state ', istate
    write(iw,'(4x,a)')      '     (using relaxed one-particle density)       '
    write(iw,'(4x,a)')      '================================================='
    call flush(iw)

    nat = ubound(basis%atoms%zn,1)

    allocate(mints(nbf2,19), source=0.0d0, stat=ok)
    if (ok /= 0) call show_message('Cannot allocate mints', WITH_ABORT)

    allocate(dens_ex(nbf2), source=0.0d0, stat=ok)
    if (ok /= 0) call show_message('Cannot allocate dens_ex', WITH_ABORT)


!   Ground-state density
    call tagarray_get_data(infos%dat, OQP_DM_A, dmat_a, status)
    call check_status(status, module_name, subroutine_name, OQP_DM_A)
    dens_ex = dmat_a

    if (urohf) then
      call tagarray_get_data(infos%dat, OQP_DM_B, dmat_b, status)
      call check_status(status, module_name, subroutine_name, OQP_DM_B)
      dens_ex = dens_ex + dmat_b
    end if

    call tagarray_get_data(infos%dat, OQP_TD_P, td_p, status)
    call check_status(status, module_name, subroutine_name, OQP_TD_P)
    dens_ex = dens_ex + td_p(:,1) + td_p(:,2)

    com = 0
    do i = 1, nat
      com = com + basis%atoms%xyz(:,i)*basis%atoms%mass(i)
    end do
    com = com / sum(basis%atoms%mass)

    call multipole_integrals(basis, mints, com, 3)

    dip  = 0
    qxyz = 0
    oxyz = 0

    do i = 1, nat
      dr = infos%atoms%xyz(:,i) - com
      z  = infos%atoms%zn(i) - infos%basis%ecp_zn_num(i)

      dip = dip + z * dr

      qxyz(XX_) = qxyz(XX_) + z * dr(X__)*dr(X__)
      qxyz(YY_) = qxyz(YY_) + z * dr(Y__)*dr(Y__)
      qxyz(ZZ_) = qxyz(ZZ_) + z * dr(Z__)*dr(Z__)
      qxyz(XY_) = qxyz(XY_) + z * dr(X__)*dr(Y__)
      qxyz(XZ_) = qxyz(XZ_) + z * dr(X__)*dr(Z__)
      qxyz(YZ_) = qxyz(YZ_) + z * dr(Y__)*dr(Z__)

      oxyz(XXX) = oxyz(XXX) + z * dr(X__)*dr(X__)*dr(X__)
      oxyz(YYY) = oxyz(YYY) + z * dr(Y__)*dr(Y__)*dr(Y__)
      oxyz(ZZZ) = oxyz(ZZZ) + z * dr(Z__)*dr(Z__)*dr(Z__)
      oxyz(XXY) = oxyz(XXY) + z * dr(X__)*dr(X__)*dr(Y__)
      oxyz(XXZ) = oxyz(XXZ) + z * dr(X__)*dr(X__)*dr(Z__)
      oxyz(YYX) = oxyz(YYX) + z * dr(Y__)*dr(Y__)*dr(X__)
      oxyz(YYZ) = oxyz(YYZ) + z * dr(Y__)*dr(Y__)*dr(Z__)
      oxyz(ZZX) = oxyz(ZZX) + z * dr(Z__)*dr(Z__)*dr(X__)
      oxyz(ZZY) = oxyz(ZZY) + z * dr(Z__)*dr(Z__)*dr(Y__)
      oxyz(XYZ) = oxyz(XYZ) + z * dr(X__)*dr(Y__)*dr(Z__)
    end do

    do i = 1, 3
      dip(i) = dip(i) - traceprod_sym_packed(mints(:,i), dens_ex, nbf)
    end do

    do i = 1, 6
      qxyz(i) = qxyz(i) - traceprod_sym_packed(mints(:,3+i), dens_ex, nbf)
    end do

    do i = 1, 10
      oxyz(i) = oxyz(i) - traceprod_sym_packed(mints(:,9+i), dens_ex, nbf)
    end do


    dip = dip*AU_TO_DEBYE

    quad(X__,X__) = 0.5*(2*qxyz(XX_) - (qxyz(YY_)+qxyz(ZZ_)))
    quad(X__,Y__) = 0.5*(3*qxyz(XY_)                        )
    quad(X__,Z__) = 0.5*(3*qxyz(XZ_)                        )
    quad(Y__,Y__) = 0.5*(2*qxyz(YY_) - (qxyz(XX_)+qxyz(ZZ_)))
    quad(Y__,Z__) = 0.5*(3*qxyz(YZ_)                        )
    quad(Z__,Z__) = 0.5*(2*qxyz(ZZ_) - (qxyz(XX_)+qxyz(YY_)))
    call triangular_to_full(quad,3,'u')
    quad = quad*AU_TO_BUCK

    oct(XXX) = 0.5*(2*oxyz(XXX) - 3*(oxyz(YYX)+oxyz(ZZX)))
    oct(YYY) = 0.5*(2*oxyz(YYY) - 3*(oxyz(XXY)+oxyz(ZZY)))
    oct(ZZZ) = 0.5*(2*oxyz(ZZZ) - 3*(oxyz(YYZ)+oxyz(XXZ)))
    oct(XXY) = 0.5*(4*oxyz(XXY) -   (oxyz(YYY)+oxyz(ZZY)))
    oct(XXZ) = 0.5*(4*oxyz(XXZ) -   (oxyz(YYZ)+oxyz(ZZZ)))
    oct(YYX) = 0.5*(4*oxyz(YYX) -   (oxyz(XXX)+oxyz(ZZX)))
    oct(YYZ) = 0.5*(4*oxyz(YYZ) -   (oxyz(XXZ)+oxyz(ZZZ)))
    oct(ZZX) = 0.5*(4*oxyz(ZZX) -   (oxyz(XXX)+oxyz(YYX)))
    oct(ZZY) = 0.5*(4*oxyz(ZZY) -   (oxyz(XXY)+oxyz(YYY)))
    oct(XYZ) = 0.5*(5*oxyz(XYZ)                          )

    oct = oct*AU_TO_OCT


    write(iw,'(/1x,a,i4)') 'Excited state index: ', istate
    write(iw,'( 1x,a)')    'Density used: relaxed (P_ground + dP_relaxed)'

    write(iw,'(/1x,a)') 'At point (Bohr):'
    write(iw,'(4x,40x,3(a8,7x))') 'X', 'Y', 'Z'
    write(iw,'(4x,40x,3f15.8,sp,f12.5)') com

    write(iw,'(/4x,a)') 'electric charge (a.u.):'
    write(iw,'(4x,34x,a6,sp,f12.5)') 'C_o', real(infos%mol_prop%charge)

    write(iw,'(/4x,a)') 'electric dipole (Debye) - excited state:'
    write(iw,'(4x,40x,3(a8,7x),a10)') 'X', 'Y', 'Z', 'Norm'
    write(iw,'(4x,34x,a6,4f15.8)') 'D_ex', dip, norm2(dip)

    write(iw,'(/4x,a)') 'electric quadrupole (Buckingham) - excited state:'
    write(iw,'(4x,40x,3(a8,7x))') 'X', 'Y', 'Z'
    write(iw,'(4x,34x,a6,3F15.8)') 'Q_X', quad(:,1)
    write(iw,'(4x,34x,a6,3f15.8)') 'Q_Y', quad(:,2)
    write(iw,'(4x,34x,a6,3f15.8)') 'Q_Z', quad(:,3)

    write(iw,'(/4x,a)') 'electric octupole (Buckingham*Angstrom) - excited state:'
    write(iw,'(4x,34x,a6,*(F15.8))') 'O_xxx', oct(XXX)
    write(iw,'(4x,34x,a6,*(F15.8))') 'O_xxy', oct(XXY)
    write(iw,'(4x,34x,a6,*(F15.8))') 'O_xxz', oct(XXZ)
    write(iw,'(4x,34x,a6,*(F15.8))') 'O_xyy', oct(YYX)
    write(iw,'(4x,34x,a6,*(F15.8))') 'O_xyz', oct(XYZ)
    write(iw,'(4x,34x,a6,*(F15.8))') 'O_xzz', oct(ZZX)
    write(iw,'(4x,34x,a6,*(F15.8))') 'O_yyy', oct(YYY)
    write(iw,'(4x,34x,a6,*(F15.8))') 'O_yyz', oct(YYZ)
    write(iw,'(4x,34x,a6,*(F15.8))') 'O_zzy', oct(ZZY)
    write(iw,'(4x,34x,a6,*(F15.8))') 'O_zzz', oct(ZZZ)

    close(iw)

    deallocate(mints, dens_ex)

  end subroutine electric_moments_excited

end module electric_moments_mod