From 92294cf97372c8d1929654ec50349f2d96c50697 Mon Sep 17 00:00:00 2001
From: Kevin Gasperich <kgasperich@gmail.com>
Date: Tue, 17 Mar 2020 10:30:34 -0500
Subject: [PATCH] cleaner ao ortho canonical for kpts

---
 src/ao_one_e_ints/ao_ortho_cano_kpts.irp.f | 72 ++++++++++++++++++----
 1 file changed, 61 insertions(+), 11 deletions(-)

diff --git a/src/ao_one_e_ints/ao_ortho_cano_kpts.irp.f b/src/ao_one_e_ints/ao_ortho_cano_kpts.irp.f
index 3e3a371a..1cc4ba9d 100644
--- a/src/ao_one_e_ints/ao_ortho_cano_kpts.irp.f
+++ b/src/ao_one_e_ints/ao_ortho_cano_kpts.irp.f
@@ -1,19 +1,69 @@
 !todo: add kpts
 
-BEGIN_PROVIDER [ integer, ao_cart_to_sphe_num_per_kpt ]
+ BEGIN_PROVIDER [ complex*16, ao_cart_to_sphe_coef_kpts, (ao_num_per_kpt,ao_num_per_kpt)]
+&BEGIN_PROVIDER [ integer, ao_cart_to_sphe_num_per_kpt ]
   implicit none
-  ao_cart_to_sphe_num_per_kpt = ao_cart_to_sphe_num / kpt_num
-END_PROVIDER
+  BEGIN_DOC
+! Coefficients to go from cartesian to spherical coordinates in the current
+! basis set
+  END_DOC
+  integer :: i
+  integer, external              :: ao_power_index
+  integer                        :: ibegin,j,k
+  integer                        :: prev
+  prev = 0
+  ao_cart_to_sphe_coefi_kpts(:,:) = (0.d0,0.d0)
+  ! Assume order provided by ao_power_index
+  i = 1
+  ao_cart_to_sphe_num_per_kpt = 0
+  do while (i <= ao_num_per_kpt)
+    select case ( ao_l(i) )
+      case (0)
+        ao_cart_to_sphe_num_per_kpt += 1
+        ao_cart_to_sphe_coef_kpts(i,ao_cart_to_sphe_num_per_kpt) = (1.d0,0.d0)
+        i += 1
+      BEGIN_TEMPLATE
+      case ($SHELL)
+        if (ao_power(i,1) == $SHELL) then
+          do k=1,size(cart_to_sphe_$SHELL,2)
+            do j=1,size(cart_to_sphe_$SHELL,1)
+              ao_cart_to_sphe_coef_kpts(i+j-1,ao_cart_to_sphe_num_per_kpt+k) = dcmplx(cart_to_sphe_$SHELL(j,k),0.d0)
+            enddo
+          enddo
+          i += size(cart_to_sphe_$SHELL,1)
+          ao_cart_to_sphe_num_per_kpt += size(cart_to_sphe_$SHELL,2)
+        endif
+      SUBST [ SHELL ]
+        1;;
+        2;;
+        3;;
+        4;;
+        5;;
+        6;;
+        7;;
+        8;;
+        9;;
+      END_TEMPLATE
+      case default
+        stop 'Error in ao_cart_to_sphe_kpts : angular momentum too high'
+    end select
+  enddo
 
-BEGIN_PROVIDER [ complex*16, ao_cart_to_sphe_coef_kpts, (ao_num_per_kpt,ao_cart_to_sphe_num_per_kpt) ]
- implicit none
- BEGIN_DOC
- ! complex version of ao_cart_to_sphe_coef for one k-point
- END_DOC
- call zlacp2('A',ao_num_per_kpt,ao_cart_to_sphe_num_per_kpt, &
-        ao_cart_to_sphe_coef,size(ao_cart_to_sphe_coef,1), &
-        ao_cart_to_sphe_coef_kpts,size(ao_cart_to_sphe_coef_kpts,1))
 END_PROVIDER
+!BEGIN_PROVIDER [ integer, ao_cart_to_sphe_num_per_kpt ]
+!  implicit none
+!  ao_cart_to_sphe_num_per_kpt = ao_cart_to_sphe_num / kpt_num
+!END_PROVIDER
+!
+!BEGIN_PROVIDER [ complex*16, ao_cart_to_sphe_coef_kpts, (ao_num_per_kpt,ao_cart_to_sphe_num_per_kpt) ]
+! implicit none
+! BEGIN_DOC
+! ! complex version of ao_cart_to_sphe_coef for one k-point
+! END_DOC
+! call zlacp2('A',ao_num_per_kpt,ao_cart_to_sphe_num_per_kpt, &
+!        ao_cart_to_sphe_coef,size(ao_cart_to_sphe_coef,1), &
+!        ao_cart_to_sphe_coef_kpts,size(ao_cart_to_sphe_coef_kpts,1))
+!END_PROVIDER
 
 BEGIN_PROVIDER [ complex*16, ao_cart_to_sphe_overlap_kpts, (ao_cart_to_sphe_num_per_kpt,ao_cart_to_sphe_num_per_kpt,kpt_num) ]
   implicit none