added new stopping criterion in CASSCF

devel/casscf/EZFIO.cfg

@@ -11,23 +11,23 @@ interface: ezfio
 size: (determinants.n_states)
 
 [state_following_casscf]
-type: logical                                                                                                                                
+type: logical
 doc: If |true|, the CASSCF will try to follow the guess CI vector and orbitals
 interface: ezfio,provider,ocaml
 default: False
 
 
 [diag_hess_cas]
-type: logical                                                                                                                                
+type: logical
 doc: If |true|, only the DIAGONAL part of the hessian is retained for the CASSCF
 interface: ezfio,provider,ocaml
 default: False
 
 [hess_cv_cv]
-type: logical                                                                                                                                
-doc: If |true|, the core-virtual - core-virtual part of the hessian is computed 
+type: logical
+doc: If |true|, the core-virtual - core-virtual part of the hessian is computed
 interface: ezfio,provider,ocaml
-default: True 
+default: True
 
 
 [level_shift_casscf]
@@ -41,11 +41,17 @@ default: 0.005
 type: logical
 doc: If true, the two-rdm are computed with a fast algo
 interface: ezfio,provider,ocaml
-default: True 
+default: True
 
 [criterion_casscf]
 type: character*(32)
-doc: choice of the criterion for the convergence of the casscf: can be  energy  or  gradients  
+doc: choice of the criterion for the convergence of the casscf: can be  energy  or  gradients  or e_pt2
 interface: ezfio, provider, ocaml
-default: energy
+default: e_pt2
+
+[thresh_casscf]
+type: Threshold
+doc: Threshold on the convergence of the CASCF energy.
+interface: ezfio,provider,ocaml
+default: 1.e-06
 

devel/casscf/casscf.irp.f

@@ -10,13 +10,15 @@ program casscf
   SOFT_TOUCH pt2_max
   n_det_max_full = 500
   touch n_det_max_full
+  pt2_relative_error = 0.02
+  touch pt2_relative_error 
   call run_stochastic_cipsi
   call run
 end
 
 subroutine run
   implicit none
-  double precision               :: energy_old, energy
+  double precision               :: energy_old, energy, pt2_max_before, ept2_before,delta_E
   logical                        :: converged,state_following_casscf_save
   integer                        :: iteration
   converged = .False.
@@ -27,31 +29,50 @@ subroutine run
   state_following_casscf_save = state_following_casscf
   state_following_casscf = .True.
   touch state_following_casscf
+  ept2_before = 0.d0
   do while (.not.converged)
     call run_stochastic_cipsi
     energy_old = energy
     energy = eone+etwo+ecore
+    pt2_max_before = pt2_max
 
     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
+     double precision :: E_PT2, PT2
+     call ezfio_get_casscf_energy_pt2(E_PT2)
+     call ezfio_get_casscf_energy(PT2)
+     PT2 -= E_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
+
+    print*,''
     call write_double(6,norm_grad_vec2,'Norm of gradients = ')
     call write_double(6,norm_grad_vec2_tab(1), '     Core-active  gradients = ')
     call write_double(6,norm_grad_vec2_tab(2), '     Core-virtual gradients = ')
     call write_double(6,norm_grad_vec2_tab(3), '   Active-virtual gradients = ')
+    print*,''
     call write_double(6,energy_improvement, 'Predicted energy improvement = ')
 
     if(criterion_casscf == "energy")then
      converged = dabs(energy_improvement) < thresh_scf
     else if (criterion_casscf == "gradients")then
      converged = norm_grad_vec2 < thresh_scf
-    else 
-     converged = dabs(energy_improvement) < thresh_scf
+    else if (criterion_casscf == "e_pt2")then
+     delta_E = dabs(E_PT2 - ept2_before)
+     converged = dabs(delta_E) < thresh_casscf
     endif
-    pt2_max = dabs(energy_improvement / pt2_relative_error)
+    ept2_before = E_PT2
+    pt2_max = dabs(energy_improvement / (pt2_relative_error))
+    pt2_max = min(pt2_max, pt2_max_before)
+    print*,''
+    call write_double(6,pt2_max, 'PT2_MAX for next iteration = ')
 
     mo_coef = NewOrbs
-    mo_occ  = occnum  
+    mo_occ  = occnum
     call save_mos
     if(.not.converged)then
      iteration += 1
@@ -60,9 +81,9 @@ subroutine run
      psi_coef = psi_coef_sorted
      read_wf = .True.
      call clear_mo_map
-     SOFT_TOUCH mo_coef N_det pt2_max  psi_det psi_coef 
+     SOFT_TOUCH mo_coef N_det pt2_max  psi_det psi_coef
      if(iteration .gt. 3)then
-      state_following_casscf = state_following_casscf_save 
+      state_following_casscf = state_following_casscf_save
       soft_touch state_following_casscf
      endif
     endif
@@ -70,3 +91,5 @@ subroutine run
   enddo
 
 end
+
+