From 1b9a75f4886a4112920da3c4f611e19bf35cae12 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Mon, 12 Feb 2024 18:18:53 +0100
Subject: [PATCH 1/2] Fixed pseudo-inverse (extrapolations)

---
 src/mol_properties/EZFIO.cfg                  |  7 +++
 .../print_e_components.irp.f                  |  0
 src/mol_properties/print_mol_properties.irp.f |  7 ++-
 src/utils/linear_algebra.irp.f                | 44 +++++++++----------
 4 files changed, 34 insertions(+), 24 deletions(-)
 rename src/{two_body_rdm => mol_properties}/print_e_components.irp.f (100%)

diff --git a/src/mol_properties/EZFIO.cfg b/src/mol_properties/EZFIO.cfg
index 35a095fb..3ddba227 100644
--- a/src/mol_properties/EZFIO.cfg
+++ b/src/mol_properties/EZFIO.cfg
@@ -21,3 +21,10 @@ type: logical
 doc: If true and N_states > 1, the oscillator strength will be computed
 interface: ezfio,provider,ocaml
 default: false
+
+[calc_energy_components]
+type: logical
+doc: If true, the components of the energy (1e, 2e, kinetic) will be computed
+interface: ezfio,provider,ocaml
+default: false
+
diff --git a/src/two_body_rdm/print_e_components.irp.f b/src/mol_properties/print_e_components.irp.f
similarity index 100%
rename from src/two_body_rdm/print_e_components.irp.f
rename to src/mol_properties/print_e_components.irp.f
diff --git a/src/mol_properties/print_mol_properties.irp.f b/src/mol_properties/print_mol_properties.irp.f
index 3753a3dd..00ccb826 100644
--- a/src/mol_properties/print_mol_properties.irp.f
+++ b/src/mol_properties/print_mol_properties.irp.f
@@ -6,6 +6,11 @@ subroutine print_mol_properties()
   ! Run the propertie calculations
   END_DOC
 
+  ! Energy components
+  if (calc_energy_components) then
+    call print_energy_components
+  endif
+
   ! Electric dipole moment
   if (calc_dipole_moment) then
     call print_dipole_moment
@@ -18,7 +23,7 @@ subroutine print_mol_properties()
 
   ! Oscillator strength
   if (calc_osc_str .and. N_states > 1) then
-    call print_oscillator_strength 
+    call print_oscillator_strength
   endif
 
 end
diff --git a/src/utils/linear_algebra.irp.f b/src/utils/linear_algebra.irp.f
index c9d0be72..76b280b7 100644
--- a/src/utils/linear_algebra.irp.f
+++ b/src/utils/linear_algebra.irp.f
@@ -1377,31 +1377,29 @@ subroutine get_pseudo_inverse(A, LDA, m, n, C, LDC, cutoff)
     enddo
   endif
 
-  print*, ' n_svd = ', n_svd
-
-  !$OMP PARALLEL       &
-  !$OMP DEFAULT (NONE) &
-  !$OMP PRIVATE (i, j) &
-  !$OMP SHARED (n, n_svd, D, Vt)
-  !$OMP DO
-  do j = 1, n
-    do i = 1, n_svd
-      Vt(i,j) = D(i) * Vt(i,j)
-    enddo
-  enddo
-  !$OMP END DO
-  !$OMP END PARALLEL
-
-  call dgemm("N", "N", m, n, n_svd, 1.d0, U, m, Vt, n, 0.d0, C, LDC)
-
-!  C = 0.d0
-!  do i=1,m
-!    do j=1,n
-!      do k=1,n
-!        C(j,i) = C(j,i) + U(i,k) * D(k) * Vt(k,j)
-!      enddo
+!  !$OMP PARALLEL       &
+!  !$OMP DEFAULT (NONE) &
+!  !$OMP PRIVATE (i, j) &
+!  !$OMP SHARED (n, n_svd, D, Vt)
+!  !$OMP DO
+!  do j = 1, n
+!    do i = 1, n_svd
+!      Vt(i,j) = D(i) * Vt(i,j)
 !    enddo
 !  enddo
+!  !$OMP END DO
+!  !$OMP END PARALLEL
+
+!  call dgemm('N', 'N', n, m, n_svd, 1.d0, Vt, size(Vt,1), U, size(U,1), 0.d0, C, size(C,1))
+
+  C = 0.d0
+  do i=1,m
+    do j=1,n
+      do k=1,n_svd
+        C(j,i) = C(j,i) + U(i,k) * D(k) * Vt(k,j)
+      enddo
+    enddo
+  enddo
 
   deallocate(U,D,Vt,work,A_tmp)
 

From d619c621fcd00e35d8bb8c2f956486044366a6d0 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Mon, 12 Feb 2024 18:21:59 +0100
Subject: [PATCH 2/2] DGEMM in pseudo-inverse

---
 src/utils/linear_algebra.irp.f | 42 +++++++++++++++++-----------------
 1 file changed, 21 insertions(+), 21 deletions(-)

diff --git a/src/utils/linear_algebra.irp.f b/src/utils/linear_algebra.irp.f
index 76b280b7..2db47092 100644
--- a/src/utils/linear_algebra.irp.f
+++ b/src/utils/linear_algebra.irp.f
@@ -1377,29 +1377,29 @@ subroutine get_pseudo_inverse(A, LDA, m, n, C, LDC, cutoff)
     enddo
   endif
 
-!  !$OMP PARALLEL       &
-!  !$OMP DEFAULT (NONE) &
-!  !$OMP PRIVATE (i, j) &
-!  !$OMP SHARED (n, n_svd, D, Vt)
-!  !$OMP DO
-!  do j = 1, n
-!    do i = 1, n_svd
-!      Vt(i,j) = D(i) * Vt(i,j)
-!    enddo
-!  enddo
-!  !$OMP END DO
-!  !$OMP END PARALLEL
-
-!  call dgemm('N', 'N', n, m, n_svd, 1.d0, Vt, size(Vt,1), U, size(U,1), 0.d0, C, size(C,1))
-
-  C = 0.d0
-  do i=1,m
-    do j=1,n
-      do k=1,n_svd
-        C(j,i) = C(j,i) + U(i,k) * D(k) * Vt(k,j)
-      enddo
+  !$OMP PARALLEL       &
+  !$OMP DEFAULT (NONE) &
+  !$OMP PRIVATE (i, j) &
+  !$OMP SHARED (n, n_svd, D, Vt)
+  !$OMP DO
+  do j = 1, n
+    do i = 1, n_svd
+      Vt(i,j) = D(i) * Vt(i,j)
     enddo
   enddo
+  !$OMP END DO
+  !$OMP END PARALLEL
+
+  call dgemm('T', 'T', n, m, n_svd, 1.d0, Vt, size(Vt,1), U, size(U,1), 0.d0, C, size(C,1))
+
+!  C = 0.d0
+!  do i=1,m
+!    do j=1,n
+!      do k=1,n_svd
+!        C(j,i) = C(j,i) + U(i,k) * D(k) * Vt(k,j)
+!      enddo
+!    enddo
+!  enddo
 
   deallocate(U,D,Vt,work,A_tmp)