working on casscf_cipsi

src/casscf_cipsi/EZFIO.cfg

@@ -73,3 +73,9 @@ type: logical
 doc: If |true|, the pt2_max value in the CIPSI iterations will automatically adapt, otherwise it is fixed at the value given in the EZFIO folder
 interface: ezfio,provider,ocaml
 default: True
+
+[small_active_space]
+type: logical
+doc: If |true|, the pt2_max value in the CIPSI is set to 10-10 and will not change
+interface: ezfio,provider,ocaml
+default: False

src/casscf_cipsi/README.rst

@@ -3,3 +3,39 @@ casscf
 ======
 
 |CASSCF| program with the CIPSI algorithm.
+
+Example of inputs
+-----------------
+
+a) Small active space : standard CASSCF 
+---------------------------------------
+Let's do O2 (triplet) in aug-cc-pvdz with the following geometry (xyz format, Bohr units)
+3
+
+ O           0.0000000000        0.0000000000       -1.1408000000
+ O           0.0000000000        0.0000000000        1.1408000000
+
+# Create the ezfio folder 
+qp create_ezfio -b aug-cc-pvdz O2.xyz -m 3 -a -o O2_avdz
+
+# Start with an ROHF guess 
+qp run scf | tee ${EZFIO_FILE}.rohf.out
+
+# Get the ROHF energy for check 
+qp get hartree_fock energy # should be -149.4684509
+
+# Define the full valence active space: the two 1s are doubly occupied, the other 8 valence orbitals are active 
+# CASSCF(12e,10orb) 
+qp set_mo_class -c "[1-2]" -a "[3-10]" -v "[11-46]"
+
+# Specify that you want an near exact CASSCF, i.e. the CIPSI selection will stop at pt2_max = 10^-10
+qp set casscf_cipsi small_active_space True 
+# RUN THE CASSCF 
+qp run casscf | tee ${EZFIO_FILE}.casscf.out
+
+
+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 
+
+

src/casscf_cipsi/casscf.irp.f

@@ -8,17 +8,22 @@ program casscf
 !  touch no_vvvv_integrals
   n_det_max_full = 500
   touch n_det_max_full
-  pt2_relative_error = 0.04
+  if(small_active_space)then
+   pt2_relative_error = 0.00001
+  else
+   pt2_relative_error = 0.04
+  endif
   touch pt2_relative_error 
-!  call run_stochastic_cipsi
   call run
 end
 
 subroutine run
   implicit none
-  double precision               :: energy_old, energy, pt2_max_before, ept2_before,delta_E
+  double precision               :: energy_old, energy, pt2_max_before,delta_E
   logical                        :: converged,state_following_casscf_cipsi_save
-  integer                        :: iteration
+  integer                        :: iteration,istate
+  double precision, allocatable :: E_PT2(:), PT2(:), Ev(:), ept2_before(:)
+  allocate(E_PT2(N_states), PT2(N_states), Ev(N_states), ept2_before(N_states))
   converged = .False.
 
   energy = 0.d0
@@ -28,13 +33,19 @@ subroutine run
   state_following_casscf = .True.
   touch state_following_casscf
   ept2_before = 0.d0
-  if(adaptive_pt2_max)then
+  if(small_active_space)then
-   pt2_max = 0.005
+   pt2_max = 1.d-10
    SOFT_TOUCH pt2_max
+  else
+   if(adaptive_pt2_max)then
+    pt2_max = 0.005
+    SOFT_TOUCH pt2_max
+   endif
   endif
   do while (.not.converged)
     print*,'pt2_max = ',pt2_max
-    call run_stochastic_cipsi
+    call run_stochastic_cipsi(Ev,PT2)
+    E_PT2(1:N_states) = Ev(1:N_states) + PT2(1:N_states)
     energy_old = energy
     energy = eone+etwo+ecore
     pt2_max_before = pt2_max
@@ -42,15 +53,13 @@ subroutine run
     call write_time(6)
     call write_int(6,iteration,'CAS-SCF iteration = ')
     call write_double(6,energy,'CAS-SCF energy = ')
-    if(n_states == 1)then
+!    if(n_states == 1)then
-     double precision :: E_PT2, PT2
+!     call ezfio_get_casscf_cipsi_energy_pt2(E_PT2)
-     call ezfio_get_casscf_cipsi_energy_pt2(E_PT2)
+!     call ezfio_get_casscf_cipsi_energy(PT2)
-     call ezfio_get_casscf_cipsi_energy(PT2)
+     call write_double(6,E_PT2(1:N_states),'E + PT2 energy = ')
-     PT2 -= E_PT2
+     call write_double(6,PT2(1:N_states),'  PT2          = ')
-     call write_double(6,E_PT2,'E + PT2 energy = ')
-     call write_double(6,PT2,'  PT2          = ')
      call write_double(6,pt2_max,' PT2_MAX       = ')
-    endif
+!    endif
 
     print*,''
     call write_double(6,norm_grad_vec2,'Norm of gradients = ')
@@ -65,15 +74,20 @@ subroutine run
     else if (criterion_casscf == "gradients")then
      converged = norm_grad_vec2 < thresh_scf
     else if (criterion_casscf == "e_pt2")then
-     delta_E = dabs(E_PT2 - ept2_before)
+     delta_E = 0.d0
+     do istate = 1, N_states
+      delta_E += dabs(E_PT2(istate) - ept2_before(istate))
+     enddo
      converged = dabs(delta_E) < thresh_casscf
     endif
     ept2_before = E_PT2
-    if(adaptive_pt2_max)then
+    if(.not.small_active_space)then
-     pt2_max = dabs(energy_improvement / (pt2_relative_error))
+     if(adaptive_pt2_max)then
-     pt2_max = min(pt2_max, pt2_max_before)
+      pt2_max = dabs(energy_improvement / (pt2_relative_error))
-     if(n_act_orb.ge.n_big_act_orb)then
+      pt2_max = min(pt2_max, pt2_max_before)
-      pt2_max = max(pt2_max,pt2_min_casscf)
+      if(n_act_orb.ge.n_big_act_orb)then
+       pt2_max = max(pt2_max,pt2_min_casscf)
+      endif
      endif
     endif
     print*,''
@@ -94,8 +108,10 @@ subroutine run
      read_wf = .True.
      call clear_mo_map
      SOFT_TOUCH mo_coef N_det psi_det psi_coef
-     if(adaptive_pt2_max)then
+     if(.not.small_active_space)then
-       SOFT_TOUCH pt2_max  
+      if(adaptive_pt2_max)then
+        SOFT_TOUCH pt2_max  
+      endif
      endif
      if(iteration .gt. 3)then
       state_following_casscf = state_following_casscf_cipsi_save
@@ -104,6 +120,27 @@ subroutine run
     endif
 
   enddo
+     integer :: i
+!    print*,'Converged CASSCF '
+!    print*,'--------------------------'
+!    write(6,*) ' occupation numbers of orbitals '
+!    do i=1,mo_num
+!      write(6,*) i,occnum(i)
+!    end do
+!
+!     write(6,*)
+!     write(6,*) ' the diagonal of the inactive effective Fock matrix '
+!     write(6,'(5(i3,F12.5))') (i,Fipq(i,i),i=1,mo_num)
+!     write(6,*)
+  print*,'Fock ROHF '
+  do i = 1, ao_num
+   write(33,*)fock_matrix_ao_alpha(i,1:ao_num)
+  enddo
+  print*,'Fock MCSCF'
+  do i = 1, ao_num
+   write(34,*)mcscf_fock_alpha(i,1:ao_num)
+  enddo
+
 
 end
 

src/casscf_cipsi/densities.irp.f

@@ -17,6 +17,35 @@ BEGIN_PROVIDER [real*8, D0tu, (n_act_orb,n_act_orb) ]
   
 END_PROVIDER
 
+ BEGIN_PROVIDER [double precision, D0tu_alpha_ao, (ao_num, ao_num)]
+&BEGIN_PROVIDER [double precision, D0tu_beta_ao, (ao_num, ao_num)]
+ implicit none
+ integer :: i,ii,j,u,t,uu,tt
+ double precision, allocatable :: D0_tmp_alpha(:,:),D0_tmp_beta(:,:)
+ allocate(D0_tmp_alpha(mo_num, mo_num),D0_tmp_beta(mo_num, mo_num))
+ D0_tmp_beta = 0.d0
+ D0_tmp_alpha = 0.d0
+ do i = 1, n_core_inact_orb
+  ii = list_core_inact(i)
+  D0_tmp_alpha(ii,ii) = 1.d0
+  D0_tmp_beta(ii,ii) = 1.d0
+ enddo
+ print*,'Diagonal elements of the 1RDM in the active space'
+ do u=1,n_act_orb
+   uu = list_act(u)
+   print*,uu,one_e_dm_mo_alpha_average(uu,uu),one_e_dm_mo_beta_average(uu,uu)
+   do t=1,n_act_orb
+    tt = list_act(t)
+    D0_tmp_alpha(tt,uu) = one_e_dm_mo_alpha_average(tt,uu)
+    D0_tmp_beta(tt,uu) = one_e_dm_mo_beta_average(tt,uu)
+   enddo
+ enddo
+
+ call mo_to_ao_no_overlap(D0_tmp_alpha,mo_num,D0tu_alpha_ao,ao_num)
+ call mo_to_ao_no_overlap(D0_tmp_beta,mo_num,D0tu_beta_ao,ao_num)
+
+END_PROVIDER 
+
 BEGIN_PROVIDER [real*8, P0tuvx, (n_act_orb,n_act_orb,n_act_orb,n_act_orb) ]
    BEGIN_DOC
    ! The second-order density matrix in the basis of the starting MOs ONLY IN THE RANGE OF ACTIVE MOS

src/casscf_cipsi/mcscf_fock.irp.f

@@ -77,4 +77,15 @@ BEGIN_PROVIDER [real*8, Fapq, (mo_num,mo_num) ]
    
 END_PROVIDER
  
- 
+ BEGIN_PROVIDER [ double precision, mcscf_fock_alpha, (ao_num, ao_num)] 
+&BEGIN_PROVIDER [ double precision, mcscf_fock_beta, (ao_num, ao_num)] 
+ implicit none
+ BEGIN_DOC
+  ! mcscf_fock_alpha are set to usual Fock like operator but computed with the MCSCF densities 
+ END_DOC
+ SCF_density_matrix_ao_alpha = D0tu_alpha_ao
+ SCF_density_matrix_ao_beta = D0tu_beta_ao
+ soft_touch SCF_density_matrix_ao_alpha SCF_density_matrix_ao_beta 
+ mcscf_fock_beta = fock_matrix_ao_beta
+ mcscf_fock_alpha = fock_matrix_ao_alpha
+END_PROVIDER 

src/cipsi/stochastic_cipsi.irp.f

@@ -1,10 +1,11 @@
-subroutine run_stochastic_cipsi
+subroutine run_stochastic_cipsi(Ev,PT2) 
   use selection_types
   implicit none
   BEGIN_DOC
 ! Selected Full Configuration Interaction with Stochastic selection and PT2.
   END_DOC
   integer                        :: i,j,k
+  double precision, intent(out)  :: Ev(N_states), PT2(N_states) 
   double precision, allocatable  :: zeros(:)
   integer                        :: to_select
   type(pt2_type)                 :: pt2_data, pt2_data_err
@@ -139,6 +140,8 @@ subroutine run_stochastic_cipsi
     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/fci/fci.irp.f

@@ -41,8 +41,10 @@ program fci
 
     write(json_unit,json_array_open_fmt) 'fci'
 
+    double precision, allocatable :: Ev(:),PT2(:)
+    allocate(Ev(N_states), PT2(N_state))
     if (do_pt2) then
-      call run_stochastic_cipsi
+      call run_stochastic_cipsi(Ev,PT2)
     else
       call run_cipsi
     endif

src/mo_optimization/cipsi_orb_opt.irp.f

@@ -11,11 +11,13 @@ subroutine run_optimization
   implicit none
 
   double precision :: e_cipsi, e_opt, delta_e
+  double precision, allocatable :: Ev(:),PT2(:)
   integer :: nb_iter,i
   logical :: not_converged
   character (len=100) :: filename
 
   PROVIDE psi_det psi_coef mo_two_e_integrals_in_map ao_pseudo_integrals
+  allocate(Ev(N_states),PT2(N_states))
 
   not_converged = .True.
   nb_iter = 0
@@ -38,7 +40,7 @@ subroutine run_optimization
       print*,'' 
       print*,'********** cipsi step **********'
       ! cispi calculation
-      call run_stochastic_cipsi
+      call run_stochastic_cipsi(Ev,PT2)
 
       ! State average energy after the cipsi step
       call state_average_energy(e_cipsi)