improved casscf and added README.rst

src/casscf_cipsi/README.rst

@@ -32,12 +32,16 @@ qp set_mo_class -c "[1-2]" -a "[3-10]" -v "[11-46]"
 qp set casscf_cipsi small_active_space True 
 # RUN THE CASSCF 
 qp run casscf | tee ${EZFIO_FILE}.casscf.out
+# you should find around -149.7243542
 
 
 b) Large active space : Exploit the selected CI in the active space 
 -------------------------------------------------------------------
-Let us start from the small active space calculation orbitals and add another shell of 
+#Let us start from the small active space calculation orbitals and add another 10 virtuals: CASSCF(12e,20orb)
+qp set_mo_class -c "[1-2]" -a "[3-20]" -v "[21-46]"
+# As this active space is larger, you unset the small_active_space feature 
+qp set casscf_cipsi small_active_space False
+# As it is a large active space, the energy convergence thereshold is set to be 0.0001
+qp run casscf | tee ${EZFIO_FILE}.casscf_large.out
+# you should find around -149.9046
 
-
-
-TODO : print FOCK MCSCF NEW in the MO BASIS AT THE END OF THE CASSCF 

src/casscf_cipsi/casscf.irp.f

@@ -11,6 +11,7 @@ program casscf
   if(small_active_space)then
    pt2_relative_error = 0.00001
   else
+   thresh_scf = 1.d-4
    pt2_relative_error = 0.04
   endif
   touch pt2_relative_error 
@@ -45,6 +46,7 @@ subroutine run
   do while (.not.converged)
     print*,'pt2_max = ',pt2_max
     call run_stochastic_cipsi(Ev,PT2)
+    print*,'Ev,PT2',Ev(1),PT2(1)
     E_PT2(1:N_states) = Ev(1:N_states) + PT2(1:N_states)
     energy_old = energy
     energy = eone+etwo+ecore

src/cipsi/stochastic_cipsi.irp.f

@@ -80,12 +80,14 @@ subroutine run_stochastic_cipsi(Ev,PT2)
     to_select = max(N_states_diag, to_select)
 
 
+    Ev(1:N_states) = psi_energy_with_nucl_rep(1:N_states)
     call pt2_dealloc(pt2_data)
     call pt2_dealloc(pt2_data_err)
     call pt2_alloc(pt2_data, N_states)
     call pt2_alloc(pt2_data_err, N_states)
     call ZMQ_pt2(psi_energy_with_nucl_rep,pt2_data,pt2_data_err,relative_error,to_select) ! Stochastic PT2 and selection
 
+    PT2(1:N_states) = pt2_data % pt2(1:N_states)
     correlation_energy_ratio = (psi_energy_with_nucl_rep(1) - hf_energy_ref)  /     &
                     (psi_energy_with_nucl_rep(1) + pt2_data % rpt2(1) - hf_energy_ref)
     correlation_energy_ratio = min(1.d0,correlation_energy_ratio)
@@ -140,8 +142,6 @@ subroutine run_stochastic_cipsi(Ev,PT2)
     call print_mol_properties()
     call write_cipsi_json(pt2_data,pt2_data_err)
   endif
-  Ev(1:N_states) = psi_energy_with_nucl_rep(1:N_states)
-  PT2(1:N_states) = pt2_data % pt2(1:N_states)
   call pt2_dealloc(pt2_data)
   call pt2_dealloc(pt2_data_err)
 

src/two_body_rdm/state_av_act_2rdm.irp.f

@@ -123,7 +123,7 @@
  state_av_act_2_rdm_spin_trace_mo = state_av_act_2_rdm_ab_mo & 
                                   + state_av_act_2_rdm_aa_mo & 
                                   + state_av_act_2_rdm_bb_mo   
-
+!
 ! call orb_range_2_rdm_state_av_openmp(state_av_act_2_rdm_spin_trace_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1))
 
  call wall_time(wall_2)

src/two_body_rdm/test_2_rdm.irp.f

@@ -2,7 +2,7 @@ program test_2_rdm
  implicit none
  read_wf = .True.
  touch read_wf 
-! call routine_active_only
+ call routine_active_only
  call routine_full_mos
 end