From 274e903d3c5543a6d5446fb8b52271c7ff896246 Mon Sep 17 00:00:00 2001
From: eginer <giner.emmanuel@gmail.com>
Date: Fri, 24 Feb 2023 21:38:09 +0100
Subject: [PATCH] added spin density

---
 src/tc_bi_ortho/31.tc_bi_ortho.bats      | 49 ++++++++++++++++++++++++
 src/tc_bi_ortho/print_tc_spin_dens.irp.f | 16 ++++++++
 src/tc_bi_ortho/spin_mulliken.irp.f      | 14 +++++--
 src/tc_bi_ortho/test_spin_dens.irp.f     | 30 +++++++++++++++
 src/tc_keywords/EZFIO.cfg                |  6 +++
 5 files changed, 112 insertions(+), 3 deletions(-)
 create mode 100644 src/tc_bi_ortho/31.tc_bi_ortho.bats
 create mode 100644 src/tc_bi_ortho/print_tc_spin_dens.irp.f
 create mode 100644 src/tc_bi_ortho/test_spin_dens.irp.f

diff --git a/src/tc_bi_ortho/31.tc_bi_ortho.bats b/src/tc_bi_ortho/31.tc_bi_ortho.bats
new file mode 100644
index 00000000..f5b9d8c0
--- /dev/null
+++ b/src/tc_bi_ortho/31.tc_bi_ortho.bats
@@ -0,0 +1,49 @@
+#!/usr/bin/env bats
+
+source $QP_ROOT/tests/bats/common.bats.sh
+source $QP_ROOT/quantum_package.rc
+
+
+function run_Ne() {
+  qp set_file Ne_tc_scf 
+  qp run cisd 
+  qp run tc_bi_ortho | tee Ne_tc_scf.cisd_tc_bi_ortho.out  
+  eref=-128.77020441279302
+  energy="$(grep "eigval_right_tc_bi_orth =" Ne_tc_scf.cisd_tc_bi_ortho.out)"
+  eq $energy $eref 1e-6
+}
+
+
+@test "Ne" {
+ run_Ne 
+}
+
+
+function run_C() {
+  qp set_file C_tc_scf 
+  qp run cisd 
+  qp run tc_bi_ortho | tee C_tc_scf.cisd_tc_bi_ortho.out  
+  eref=-37.757536149952514
+  energy="$(grep "eigval_right_tc_bi_orth =" C_tc_scf.cisd_tc_bi_ortho.out)"
+  eq $energy $eref 1e-6
+}
+
+
+@test "C" {
+ run_C 
+}
+
+function run_O() {
+  qp set_file C_tc_scf 
+  qp run cisd 
+  qp run tc_bi_ortho | tee O_tc_scf.cisd_tc_bi_ortho.out  
+  eref=-74.908518517716161
+  energy="$(grep "eigval_right_tc_bi_orth =" O_tc_scf.cisd_tc_bi_ortho.out)"
+  eq $energy $eref 1e-6
+}
+
+
+@test "O" {
+ run_O 
+}
+
diff --git a/src/tc_bi_ortho/print_tc_spin_dens.irp.f b/src/tc_bi_ortho/print_tc_spin_dens.irp.f
new file mode 100644
index 00000000..8308140d
--- /dev/null
+++ b/src/tc_bi_ortho/print_tc_spin_dens.irp.f
@@ -0,0 +1,16 @@
+program test_spin_dens
+  implicit none
+  BEGIN_DOC
+! TODO : Reads psi_det in the EZFIO folder and prints out the left- and right-eigenvectors together with the energy. Saves the left-right wave functions at the end. 
+  END_DOC
+  print *, 'Hello world'
+  my_grid_becke = .True.
+  my_n_pt_r_grid = 30
+  my_n_pt_a_grid = 50
+  read_wf = .True.
+  touch read_wf
+  touch  my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
+ call tc_print_mulliken_sd
+! call test
+
+end
diff --git a/src/tc_bi_ortho/spin_mulliken.irp.f b/src/tc_bi_ortho/spin_mulliken.irp.f
index 922cc1b9..e225d299 100644
--- a/src/tc_bi_ortho/spin_mulliken.irp.f
+++ b/src/tc_bi_ortho/spin_mulliken.irp.f
@@ -7,13 +7,21 @@ BEGIN_PROVIDER [double precision, tc_spin_population, (ao_num,ao_num,N_states)]
 ! tc_spin_population(i,j) = rho_AO(alpha)(i,j) - rho_AO(beta)(i,j) * <AO_i|AO_j>
  END_DOC
  tc_spin_population = 0.d0
- do istate = 1, N_states
+ if(only_spin_tc_right)then
   do i = 1, ao_num
    do j = 1, ao_num
-    tc_spin_population(j,i,istate) = tc_spin_transition_matrix_ao(j,i,istate,istate) * ao_overlap(j,i)
+    tc_spin_population(j,i,1) = tc_spin_dens_right_only(j,i) * ao_overlap(j,i)
    enddo
   enddo
- enddo
+ else
+  do istate = 1, N_states
+   do i = 1, ao_num
+    do j = 1, ao_num
+     tc_spin_population(j,i,istate) = tc_spin_transition_matrix_ao(j,i,istate,istate) * ao_overlap(j,i)
+    enddo
+   enddo
+  enddo
+ endif
 END_PROVIDER
 
  BEGIN_PROVIDER [double precision, tc_spin_population_angular_momentum, (0:ao_l_max,N_states)]
diff --git a/src/tc_bi_ortho/test_spin_dens.irp.f b/src/tc_bi_ortho/test_spin_dens.irp.f
new file mode 100644
index 00000000..2c2f6e7c
--- /dev/null
+++ b/src/tc_bi_ortho/test_spin_dens.irp.f
@@ -0,0 +1,30 @@
+BEGIN_PROVIDER [ double precision, mo_r_coef_normalized, (ao_num,mo_num) ]
+ implicit none
+ integer :: i,j,k
+ double precision :: norm
+ do i = 1, mo_num
+  norm = 0.d0
+  do j = 1, ao_num
+   do k = 1, ao_num
+    norm += mo_r_coef(k,i) * mo_r_coef(j,i) * ao_overlap(k,j)
+   enddo
+  enddo
+  norm = 1.d0/dsqrt(norm)
+  do j = 1, ao_num
+   mo_r_coef_normalized(j,i) = mo_r_coef(j,i) * norm
+  enddo
+ enddo
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, tc_spin_dens_right_only, (ao_num, ao_num)]
+ implicit none
+ integer :: i,j,k
+ tc_spin_dens_right_only = 0.d0
+ do i = elec_beta_num+1, elec_alpha_num
+  do j = 1, ao_num
+   do k = 1, ao_num
+    tc_spin_dens_right_only(k,j) += mo_r_coef_normalized(k,i) * mo_r_coef_normalized(j,i)
+   enddo
+  enddo
+ enddo
+END_PROVIDER
diff --git a/src/tc_keywords/EZFIO.cfg b/src/tc_keywords/EZFIO.cfg
index 5d5477bc..4a1fcb9f 100644
--- a/src/tc_keywords/EZFIO.cfg
+++ b/src/tc_keywords/EZFIO.cfg
@@ -183,3 +183,9 @@ type: integer
 doc: If :: 1 then you compute the TC-PT2 the old way, :: 2 then you check with the new version but without three-body
 interface: ezfio,provider,ocaml
 default: -1
+
+[only_spin_tc_right]
+type: logical
+doc: If |true|, only the right part of WF is used to compute spin dens
+interface: ezfio,provider,ocaml
+default: False