From 3beea8d2302ce4666d190116b71988976931e2af Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Thu, 12 Jun 2014 16:34:47 +0200
Subject: [PATCH] Added mo_occ and mo_density_matrix providers

---
 src/MOs/NEEDED_MODULES    |  2 +-
 src/MOs/README.rst        | 15 +++++++++++--
 src/MOs/cholesky_mo.irp.f | 44 +++++++++++++++++++++++++++++++++++++++
 src/MOs/mos.ezfio_config  |  1 +
 src/MOs/mos.irp.f         | 26 +++++++++++++++++++++--
 src/MOs/utils.irp.f       |  1 +
 6 files changed, 84 insertions(+), 5 deletions(-)

diff --git a/src/MOs/NEEDED_MODULES b/src/MOs/NEEDED_MODULES
index 6ec1892c..9bbf1b74 100644
--- a/src/MOs/NEEDED_MODULES
+++ b/src/MOs/NEEDED_MODULES
@@ -1 +1 @@
-AOs Ezfio_files Nuclei Output Utils
+AOs Ezfio_files Nuclei Output Utils Electrons
diff --git a/src/MOs/README.rst b/src/MOs/README.rst
index 659db5d5..f27425bc 100644
--- a/src/MOs/README.rst
+++ b/src/MOs/README.rst
@@ -27,6 +27,7 @@ Needed Modules
 * `Nuclei <http://github.com/LCPQ/quantum_package/tree/master/src/Nuclei>`_
 * `Output <http://github.com/LCPQ/quantum_package/tree/master/src/Output>`_
 * `Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Utils>`_
+* `Electrons <http://github.com/LCPQ/quantum_package/tree/master/src/Electrons>`_
 
 Documentation
 =============
@@ -34,12 +35,19 @@ Documentation
 .. Do not edit this section. It was auto-generated from the
 .. NEEDED_MODULES file.
 
+`cholesky_mo <http://github.com/LCPQ/quantum_package/tree/master/src/MOs/cholesky_mo.irp.f#L1>`_
+  Cholesky decomposition of MO Density matrix to
+  generate MOs
+
+`mo_density_matrix <http://github.com/LCPQ/quantum_package/tree/master/src/MOs/cholesky_mo.irp.f#L44>`_
+  Density matrix in MO basis
+
 `mo_coef <http://github.com/LCPQ/quantum_package/tree/master/src/MOs/mos.irp.f#L22>`_
   Molecular orbital coefficients on AO basis set
   mo_coef(i,j) = coefficient of the ith ao on the jth mo
   mo_label : Label characterizing the MOS (local, canonical, natural, etc)
 
-`mo_coef_transp <http://github.com/LCPQ/quantum_package/tree/master/src/MOs/mos.irp.f#L61>`_
+`mo_coef_transp <http://github.com/LCPQ/quantum_package/tree/master/src/MOs/mos.irp.f#L60>`_
   Molecular orbital coefficients on AO basis set
 
 `mo_label <http://github.com/LCPQ/quantum_package/tree/master/src/MOs/mos.irp.f#L23>`_
@@ -47,13 +55,16 @@ Documentation
   mo_coef(i,j) = coefficient of the ith ao on the jth mo
   mo_label : Label characterizing the MOS (local, canonical, natural, etc)
 
+`mo_occ <http://github.com/LCPQ/quantum_package/tree/master/src/MOs/mos.irp.f#L78>`_
+  MO occupation numbers
+
 `mo_tot_num <http://github.com/LCPQ/quantum_package/tree/master/src/MOs/mos.irp.f#L1>`_
   Total number of molecular orbitals and the size of the keys corresponding
 
 `mo_tot_num_align <http://github.com/LCPQ/quantum_package/tree/master/src/MOs/mos.irp.f#L12>`_
   Aligned variable for dimensioning of arrays
 
-`mo_as_eigvectors_of_mo_matrix <http://github.com/LCPQ/quantum_package/tree/master/src/MOs/utils.irp.f#L21>`_
+`mo_as_eigvectors_of_mo_matrix <http://github.com/LCPQ/quantum_package/tree/master/src/MOs/utils.irp.f#L22>`_
   Undocumented
 
 `save_mos <http://github.com/LCPQ/quantum_package/tree/master/src/MOs/utils.irp.f#L1>`_
diff --git a/src/MOs/cholesky_mo.irp.f b/src/MOs/cholesky_mo.irp.f
index d587d364..d3e04dfb 100644
--- a/src/MOs/cholesky_mo.irp.f
+++ b/src/MOs/cholesky_mo.irp.f
@@ -1,5 +1,9 @@
 subroutine cholesky_mo(n,m,P,C,LDC,tol_in,rank)
  implicit none
+ BEGIN_DOC
+! Cholesky decomposition of MO Density matrix to
+! generate MOs
+ END_DOC
  integer, intent(in) :: n,m, LDC
  double precision, intent(in) :: P(LDC,n)
  double precision, intent(out) :: C(LDC,m)
@@ -37,3 +41,43 @@ subroutine cholesky_mo(n,m,P,C,LDC,tol_in,rank)
  deallocate(W,work)
 end
 
+BEGIN_PROVIDER [ double precision, mo_density_matrix, (mo_tot_num_align, mo_tot_num) ]
+ implicit none
+ BEGIN_DOC
+ ! Density matrix in MO basis
+ END_DOC
+ integer :: i,j,k
+ mo_density_matrix = 0.d0
+ do k=1,mo_tot_num
+   if (mo_occ(k) == 0) then
+     cycle
+   endif
+   do j=1,ao_num
+     do i=1,ao_num
+       mo_density_matrix(i,j) = mo_density_matrix(i,j) + &
+         mo_occ(k) * mo_coef(i,k) * mo_coef(j,k)
+     enddo
+   enddo
+ enddo
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, mo_density_matrix_virtual, (mo_tot_num_align, mo_tot_num) ]
+ implicit none
+ BEGIN_DOC
+ ! Density matrix in MO basis (virtual MOs)
+ END_DOC
+ integer :: i,j,k
+ mo_density_matrix = 0.d0
+ do k=1,mo_tot_num
+   if (mo_occ(k) == 0) then
+     cycle
+   endif
+   do j=1,ao_num
+     do i=1,ao_num
+       mo_density_matrix(i,j) = mo_density_matrix(i,j) + &
+         (2.d0-mo_occ(k)) * mo_coef(i,k) * mo_coef(j,k)
+     enddo
+   enddo
+ enddo
+END_PROVIDER
+
diff --git a/src/MOs/mos.ezfio_config b/src/MOs/mos.ezfio_config
index a0eda491..b43f23b9 100644
--- a/src/MOs/mos.ezfio_config
+++ b/src/MOs/mos.ezfio_config
@@ -2,4 +2,5 @@ mo_basis
   mo_tot_num         integer
   mo_coef            double precision  (ao_basis_ao_num,mo_basis_mo_tot_num)
   mo_label           character*(64)
+  mo_occ             double precision  (mo_basis_mo_tot_num)
 
diff --git a/src/MOs/mos.irp.f b/src/MOs/mos.irp.f
index b081d8ce..738d2c2a 100644
--- a/src/MOs/mos.irp.f
+++ b/src/MOs/mos.irp.f
@@ -32,7 +32,7 @@ END_PROVIDER
   !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: buffer
   logical                        :: exists
   PROVIDE ezfio_filename
-
+  
   !Label
   call ezfio_has_mo_basis_mo_label(exists)
   if (exists) then
@@ -40,7 +40,7 @@ END_PROVIDER
   else
     mo_label = 'no_label'
   endif
-
+  
   ! Coefs
   allocate(buffer(ao_num,mo_tot_num))
   buffer = 0.d0
@@ -75,3 +75,25 @@ BEGIN_PROVIDER [ double precision, mo_coef_transp, (mo_tot_num_align,ao_num) ]
   
 END_PROVIDER
 
+BEGIN_PROVIDER [ double precision, mo_occ, (mo_tot_num) ]
+  implicit none
+  BEGIN_DOC
+  ! MO occupation numbers
+  END_DOC
+  PROVIDE ezfio_filename
+  logical :: exists
+  call ezfio_has_mo_basis_mo_occ(exists)
+  if (exists) then
+    call ezfio_get_mo_basis_mo_occ(mo_occ)
+  else
+    mo_occ = 0.d0
+    integer :: i
+    do i=1,elec_beta_num
+      mo_occ(i) = 2.d0
+    enddo
+    do i=elec_beta_num+1,elec_alpha_num
+      mo_occ(i) = 1.d0
+    enddo
+  endif
+END_PROVIDER
+
diff --git a/src/MOs/utils.irp.f b/src/MOs/utils.irp.f
index 186f5efc..ee04ac4b 100644
--- a/src/MOs/utils.irp.f
+++ b/src/MOs/utils.irp.f
@@ -14,6 +14,7 @@ subroutine save_mos
     enddo
   enddo
   call ezfio_set_mo_basis_mo_coef(buffer)
+  call ezfio_set_mo_basis_mo_occ(mo_occ)
   deallocate (buffer)
   
 end