From b7474b4a608792568b6138e714ed2de74745e8ea Mon Sep 17 00:00:00 2001
From: Pierre-Francois Loos <pierrefrancois.loos@gmail.com>
Date: Sun, 15 Mar 2020 14:34:18 +0100
Subject: [PATCH] R and U density matrices

---
 src/eDFT/restricted_density_matrix.f90   | 35 ++++++++++++++++++
 src/eDFT/unrestricted_density_matrix.f90 | 47 ++++++++++++++++++++++++
 2 files changed, 82 insertions(+)
 create mode 100644 src/eDFT/restricted_density_matrix.f90
 create mode 100644 src/eDFT/unrestricted_density_matrix.f90

diff --git a/src/eDFT/restricted_density_matrix.f90 b/src/eDFT/restricted_density_matrix.f90
new file mode 100644
index 0000000..0f1ea16
--- /dev/null
+++ b/src/eDFT/restricted_density_matrix.f90
@@ -0,0 +1,35 @@
+subroutine restricted_density_matrix(nBas,nEns,nO,c,P)
+
+! Calculate density matrices
+
+  implicit none
+
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nEns
+  integer,intent(in)            :: nO
+  double precision,intent(in)   :: c(nBas,nBas)
+
+! Local variables
+
+  integer                       :: iEns
+
+! Output variables
+
+  double precision,intent(out)  :: P(nBas,nBas,nEns)
+
+! Ground state density matrix
+
+  iEns = 1
+  P(:,:,iEns) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO)))
+
+! Doubly-excited state density matrix
+
+  iEns = 2 
+  P(:,:,iEns) = 2d0*matmul(c(:,1:nO-1),transpose(c(:,1:nO-1))) &
+              + 2d0*matmul(c(:,nO+1:nO+1),transpose(c(:,nO+1:nO+1)))
+
+end subroutine restricted_density_matrix
diff --git a/src/eDFT/unrestricted_density_matrix.f90 b/src/eDFT/unrestricted_density_matrix.f90
new file mode 100644
index 0000000..0978142
--- /dev/null
+++ b/src/eDFT/unrestricted_density_matrix.f90
@@ -0,0 +1,47 @@
+subroutine unrestricted_density_matrix(nBas,nEns,nO,c,P)
+
+! Calculate density matrices
+
+  implicit none
+
+  include 'parameters.h'
+
+! Input variables
+
+  integer,intent(in)            :: nBas
+  integer,intent(in)            :: nEns
+  integer,intent(in)            :: nO(nspin)
+  double precision,intent(in)   :: c(nBas,nBas,nspin)
+
+! Local variables
+
+  integer                       :: ispin
+  integer                       :: iEns
+
+! Output variables
+
+  double precision,intent(out)  :: P(nBas,nBas,nspin,nEns)
+
+! Ground state density matrix
+
+  iEns = 1
+  do ispin=1,nspin
+    P(:,:,ispin,iEns) = matmul(c(:,1:nO(ispin),ispin),transpose(c(:,1:nO(ispin),ispin)))
+  end do
+
+! Singly-excited state density matrix
+
+  iEns = 2
+  P(:,:,1,iEns) = matmul(c(:,1:nO(1)-1,1),transpose(c(:,1:nO(1)-1,1))) &
+                + matmul(c(:,nO(1)+1:nO(1)+1,1),transpose(c(:,nO(1)+1:nO(1)+1,1)))
+  P(:,:,2,iEns) = matmul(c(:,1:nO(2),2),transpose(c(:,1:nO(2),2)))
+
+! Doubly-excited state density matrix
+
+  iEns = 3 
+  do ispin=1,nspin
+    P(:,:,ispin,iEns) = matmul(c(:,1:nO(ispin)-1,ispin),transpose(c(:,1:nO(ispin)-1,ispin))) &
+                      + matmul(c(:,nO(ispin)+1:nO(ispin)+1,ispin),transpose(c(:,nO(ispin)+1:nO(ispin)+1,ispin)))
+  end do
+
+end subroutine unrestricted_density_matrix