Added mo_occ and mo_density_matrix providers

src/MOs/NEEDED_MODULES

@@ -1 +1 @@
-AOs Ezfio_files Nuclei Output Utils
+AOs Ezfio_files Nuclei Output Utils Electrons

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>`_

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
+

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)
 

src/MOs/mos.irp.f

@@ -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
+

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