From 14f4a79a4bd663003ef277f1ec061c5fabf45c3b Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Wed, 21 Jul 2021 11:30:05 +0200 Subject: [PATCH] added new stopping criterion in CASSCF --- devel/casscf/EZFIO.cfg | 22 ++++++++++++++-------- devel/casscf/casscf.irp.f | 37 ++++++++++++++++++++++++++++++------- 2 files changed, 44 insertions(+), 15 deletions(-) diff --git a/devel/casscf/EZFIO.cfg b/devel/casscf/EZFIO.cfg index 95b032d..c8ce539 100644 --- a/devel/casscf/EZFIO.cfg +++ b/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 diff --git a/devel/casscf/casscf.irp.f b/devel/casscf/casscf.irp.f index 24f2f2f..03d2638 100644 --- a/devel/casscf/casscf.irp.f +++ b/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 + +