From 6025113cbcc704f42c912506e27cdacda5f94b21 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Wed, 9 Dec 2015 10:18:12 +0100
Subject: [PATCH] Canonical orthonormalization

---
 src/MOGuess/mo_ortho_canonical.irp.f          | 42 +++++++++++++++++++
 src/MOGuess/pot_mo_ortho_canonical_ints.irp.f | 25 +++++++++++
 2 files changed, 67 insertions(+)
 create mode 100644 src/MOGuess/mo_ortho_canonical.irp.f
 create mode 100644 src/MOGuess/pot_mo_ortho_canonical_ints.irp.f

diff --git a/src/MOGuess/mo_ortho_canonical.irp.f b/src/MOGuess/mo_ortho_canonical.irp.f
new file mode 100644
index 00000000..a4027555
--- /dev/null
+++ b/src/MOGuess/mo_ortho_canonical.irp.f
@@ -0,0 +1,42 @@
+
+BEGIN_PROVIDER [double precision, ao_ortho_canonical_coef, (ao_num_align,ao_num)]
+  implicit none
+  BEGIN_DOC
+! matrix of the coefficients of the mos generated by the 
+! orthonormalization by the S^{-1/2} canonical transformation of the aos
+! ao_ortho_canonical_coef(i,j) = coefficient of the ith ao on the jth ao_ortho_canonical orbital
+  END_DOC
+  integer                        :: i,j,k,l
+  tmp_matrix(:,:) = 0.d0
+  do j=1, ao_num
+    tmp_matrix(j,j) = 1.d0
+  enddo
+  call ortho_canonical(ao_overlap,ao_num_align,ao_num,ao_ortho_canonical_coef,ao_num_align,mo_tot_num)
+  SOFT_TOUCH mo_tot_num
+END_PROVIDER
+
+BEGIN_PROVIDER [double precision, ao_ortho_canonical_overlap, (ao_num_align,ao_num)]
+  implicit none
+  BEGIN_DOC
+! overlap matrix of the ao_ortho_canonical
+! supposed to be the Identity
+  END_DOC
+  integer                        :: i,j,k,l
+  double precision               :: c
+  do j=1, ao_num
+    do i=1, ao_num
+      ao_ortho_canonical_overlap(i,j) = 0.d0
+    enddo
+  enddo
+  do k=1, ao_num
+    do j=1, ao_num
+      c = 0.d0
+      do l=1, ao_num
+        c +=  ao_ortho_canonical_coef(j,l) * ao_overlap(k,l)
+      enddo
+      do i=1, ao_num
+        ao_ortho_canonical_overlap(i,j) += ao_ortho_canonical_coef(i,k) * c
+      enddo
+    enddo
+  enddo
+END_PROVIDER
diff --git a/src/MOGuess/pot_mo_ortho_canonical_ints.irp.f b/src/MOGuess/pot_mo_ortho_canonical_ints.irp.f
new file mode 100644
index 00000000..9c61ebcd
--- /dev/null
+++ b/src/MOGuess/pot_mo_ortho_canonical_ints.irp.f
@@ -0,0 +1,25 @@
+BEGIN_PROVIDER [double precision, ao_ortho_canonical_nucl_elec_integral, (mo_tot_num_align,mo_tot_num)]
+ implicit none
+ integer :: i1,j1,i,j
+ double precision :: c_i1,c_j1
+
+ ao_ortho_canonical_nucl_elec_integral = 0.d0
+ !$OMP PARALLEL DO DEFAULT(none) &
+ !$OMP PRIVATE(i,j,i1,j1,c_j1,c_i1) &
+ !$OMP SHARED(mo_tot_num,ao_num,ao_ortho_canonical_coef, &
+ !$OMP   ao_ortho_canonical_nucl_elec_integral, ao_nucl_elec_integral)
+ do i = 1, mo_tot_num
+   do j = 1, mo_tot_num
+    do i1 = 1,ao_num
+     c_i1 = ao_ortho_canonical_coef(i1,i)
+     do j1 = 1,ao_num
+       c_j1 = c_i1*ao_ortho_canonical_coef(j1,j)
+       ao_ortho_canonical_nucl_elec_integral(j,i) = ao_ortho_canonical_nucl_elec_integral(j,i) + &
+                                           c_j1 * ao_nucl_elec_integral(j1,i1)
+     enddo
+    enddo
+   enddo
+ enddo
+ !$OMP END PARALLEL DO
+END_PROVIDER
+