subroutine genr22(basis, ppairs, grotspd, shell_ids, flips, cutoffs, emu2)
implicit none
type(basis_set), intent(in) :: basis
type(int2_pair_storage), intent(in) :: ppairs
real(kind=dp), intent(inout) :: grotspd(*)
integer, intent(in) :: shell_ids(4)
integer, intent(out) :: flips(4)
type(int2_cutoffs_t), intent(in) :: cutoffs
real(kind=dp), optional :: emu2
real(kind=dp) :: a(3), b(3), c(3), d(3), p(3,3), t(3)
integer :: flips_inv(4)
integer :: angm(4)
integer :: shl_new(4)
integer :: itype, jtype
integer :: i
integer :: ji
real(kind=dp) :: acx, acz, cq, cqx, cqz, qx, qz
real(kind=dp) :: cosg, sing
real(kind=dp) :: rcd
real(kind=dp) :: tmp
real(kind=dp) :: x03, x04
integer :: nga, ngb, ngc, ngd
integer :: ppid_p, ppid_q, npp_p, npp_q
integer :: id1, id2
integer, parameter :: jtypes(*) = [ &
1, 2, 7, 2, 3, &
8, 7, 8, 10, 2, &
4, 9, 4, 5, 11, &
9, 11, 14, 7, 9, &
12, 9, 13, 15, 12, &
15, 17, 2, 4, 9, &
4, 5, 11, 9, 11, &
14, 3, 5, 13, 5, &
6, 16, 13, 16, 18, &
8, 11, 15, 11, 16, &
19, 15, 19, 20, 7, &
9, 12, 9, 13, 15, &
12, 15, 17, 8, 11, &
15, 11, 16, 19, 15, &
19, 20, 10, 14, 17, &
14, 18, 20, 17, 20, &
21]
type(rotaxis_data_t) :: rdat
angm = basis%am(shell_ids)
itype = 1 + angm(4) + 3*angm(3) + 9*angm(2) + 27*angm(1)
jtype = jtypes(itype)
flips = [1,2,3,4]
if (angm(1) > angm(2)) then
flips(1:2) = [2,1]
angm(1:2) = angm([2,1])
end if
if (angm(3) > angm(4)) then
flips(3:4) = [4,3]
angm(3:4) = angm([4,3])
end if
if (sum(angm(1:2)) > sum(angm(3:4)) .or. angm(2) > angm(4)) then
flips = flips([3,4,1,2])
angm = angm([3,4,1,2])
end if
flips_inv(flips) = [1, 2, 3, 4]
shl_new = shell_ids(flips_inv)
! If (la/=angm(1).or.lb/=angm(2)) print *, la, angm(1), lb, angm(2)
! Empty integral summation storage
call intclean(rdat,jtype)
if (present(emu2)) then
rdat%lrint = .true.
rdat%emu2 = emu2
end if
id1 = maxval(shl_new([1,2]))
id2 = minval(shl_new([1,2]))
npp_p = ppairs%ppid(1,id1*(id1-1)/2+id2)
ppid_p = ppairs%ppid(2,id1*(id1-1)/2+id2)
id1 = maxval(shl_new([3,4]))
id2 = minval(shl_new([3,4]))
npp_q = ppairs%ppid(1,id1*(id1-1)/2+id2)
ppid_q = ppairs%ppid(2,id1*(id1-1)/2+id2)
! If (npp_p == 0 .or. npp_q == 0) return
! Obtain information about shells: inew, knew, jnew, lnew
! Number of gaussians go into nga,... in common shllfo
! Shell angular quantum numbers la,... go into common shllfo
! Gaussian exponents go into arrays exa,exb,exc,exd in common shllfo
! Gaussian coefficients go into arrays csa,cpa,... in common shllfo
! Loop over gaussians in each shell
! First shell inew
! Coordinates of atoms associated with shells inew jnew knew and lnew
a = ppairs%p(:,ppid_p)-ppairs%pa(:, ppid_p)
b = ppairs%p(:,ppid_p)-ppairs%pb(:, ppid_p)
c = ppairs%p(:,ppid_q)-ppairs%pa(:, ppid_q)
d = ppairs%p(:,ppid_q)-ppairs%pb(:, ppid_q)
! Find direction cosines of penultimate axes from coordinates of ab
! P(1,1),p(1,2),... are direction cosines of axes at p. z-axis along ab
! T(1),t(2),t(3)... are direction cosines of axes at q. z-axis along cd
! Find direction cosines of ab and cd. these are local z-axes.
! If indeterminate take along space z-axis
p(:,3) = [0, 0, 1]
rdat%rab = ppairs%rab(ppid_p)
if (rdat%rab>0) then
p(:,3) = (b-a)*ppairs%uab(ppid_p)
end if
t = [0, 0, 1]
rcd = ppairs%rab(ppid_q)
if (rcd>0) then
t = (d-c)*ppairs%uab(ppid_q)
end if
! Find local y-axis as common perpendicular to ab and cd
! If indeterminate take perpendicular to ab and space z-axis
! If still indeterminate take perpendicular to ab and space x-axis
cosg = dot_product(t,p(:,3))
! Modified rotation testing.
! This fix cures the small angle problem.
p(1,2) = t(3)*p(2,3) - t(2)*p(3,3)
p(2,2) = t(1)*p(3,3) - t(3)*p(1,3)
p(3,2) = t(2)*p(1,3) - t(1)*p(2,3)
if (abs(cosg)>0.9) then
sing = norm2(p(:,2))
else
sing = sqrt(1-cosg*cosg)
end if
if (sing<1d-12) then
if (abs(p(1,3))<sqrt(0.5d0)) then
tmp = 1/sqrt(1-p(1,3)*p(1,3))
p(:,2) = [0d0, p(3,3), -p(2,3)] * tmp
else
tmp = 1/sqrt(1-p(3,3)*p(3,3))
p(:,2) = [p(2,3), -p(1,3), 0d0] * tmp
end if
else
p(:,2)= p(:,2)/sing
end if
p(1,1) = p(2,2)*p(3,3)-p(3,2)*p(2,3)
p(2,1) = p(3,2)*p(1,3)-p(1,2)*p(3,3)
p(3,1) = p(1,2)*p(2,3)-p(2,2)*p(1,3)
! Find coordinates of c relative to local axes at a
t = c - a
acx = dot_product(t,p(:,1))
rdat%acy = dot_product(t,p(:,2))
acz = dot_product(t,p(:,3))
! Set acy= 0 if close
if (abs(rdat%acy)<=acy_threshold) then
rdat%acy = 0.0_dp
rdat%acy2 = 0.0_dp
else
rdat%acy2 = rdat%acy*rdat%acy
end if
! Direction cosines of cd local axes with respect to ab local axes
! ( cosg, 0,-sing )
! ( 0, 1, 0 )
! ( sing, 0, cosg )
! Preliminary p loop
! Fill geompq with information about p in preliminary p-loop
rdat%cutoff = cutoffs%quartet_cutoff_squared
ji = 0
do i = 0, npp_p-1
if (abs(ppairs%k(ppid_p+i)*ppairs%ginv(ppid_p+i)) < cutoffs%quartet_cutoff) cycle
ji = ji+1
rdat%tx12(ji) = ppairs%g(ppid_p+i)
rdat%ty02(ji) = ppairs%alpha_b(ppid_p+i)*ppairs%ginv(ppid_p+i)*rdat%rab
rdat%sp(ji) = ppairs%k(ppid_p+i)*ppairs%ginv(ppid_p+i)
end do
rdat%ngangb = ji
! Begin q loop
do i = 0, npp_q-1
if (abs(ppairs%k(ppid_q+i)*ppairs%ginv(ppid_q+i)) < cutoffs%quartet_cutoff) cycle
x03 = ppairs%alpha_a(ppid_q+i)
x04 = ppairs%alpha_b(ppid_q+i)
rdat%x34 = ppairs%g(ppid_q+i)
rdat%x43 = ppairs%ginv(ppid_q+i)
rdat%y03 = x03*rdat%x43
rdat%y04 = x04*rdat%x43
rdat%sq = ppairs%k(ppid_q+i)*ppairs%ginv(ppid_q+i)
! Cqx = component of cq along penultimate x-axis
! Cqz = component of cq along penultimate z-axis
cq = rcd*rdat%y04
cqx = cq*sing
cqz = cq*cosg
! Find coordinates of q relative to axes at a
! Qpr is perpendicular from q to ab
rdat%aqx = acx+cqx
rdat%aqx2 = rdat%aqx*rdat%aqx
rdat%aqxy = rdat%aqx*rdat%acy
rdat%aqz = acz+cqz
rdat%qps = rdat%aqx2+rdat%acy2
! Use special fast routine for inner loops for 0000 ... 1111
call spdgen(jtype,rdat,1)
end do
qx = rcd*sing
qz = rcd*cosg
call mcdv_all(grotspd, rdat,qx, qz, jtype)
p = transpose(p)
call r30s1d(jtype, grotspd, p)
end subroutine genr22