From 6985d4d5493a6204f95a8d1d1cccbccc80c12071 Mon Sep 17 00:00:00 2001 From: eginer Date: Fri, 12 Jul 2024 18:25:17 +0200 Subject: [PATCH] the casscf does not work with mo optimization ... --- src/casscf_cipsi/EZFIO.cfg | 6 + src/casscf_cipsi/NEED | 2 +- src/casscf_cipsi/casscf.irp.f | 166 +++++++++--------- .../class.irp.f | 0 4 files changed, 94 insertions(+), 80 deletions(-) rename src/{mo_optimization_utils => mo_optimization}/class.irp.f (100%) diff --git a/src/casscf_cipsi/EZFIO.cfg b/src/casscf_cipsi/EZFIO.cfg index 18e0b6b1..5b72d906 100644 --- a/src/casscf_cipsi/EZFIO.cfg +++ b/src/casscf_cipsi/EZFIO.cfg @@ -79,3 +79,9 @@ 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 + +[act_mos_opt] +type: logical +doc: If |true|, the active orbitals are also optimized variationally +interface: ezfio,provider,ocaml +default: False diff --git a/src/casscf_cipsi/NEED b/src/casscf_cipsi/NEED index 11d1a78c..32f5ae90 100644 --- a/src/casscf_cipsi/NEED +++ b/src/casscf_cipsi/NEED @@ -3,4 +3,4 @@ selectors_full generators_cas two_body_rdm dav_general_mat -mo_optimization +mo_optimization_utils diff --git a/src/casscf_cipsi/casscf.irp.f b/src/casscf_cipsi/casscf.irp.f index dc3e2245..b64a9d8f 100644 --- a/src/casscf_cipsi/casscf.irp.f +++ b/src/casscf_cipsi/casscf.irp.f @@ -46,94 +46,101 @@ 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 - pt2_max_before = pt2_max - - call write_time(6) - call write_int(6,iteration,'CAS-SCF iteration = ') - call write_double(6,energy,'State-average CAS-SCF energy = ') -! if(n_states == 1)then -! call ezfio_get_casscf_cipsi_energy_pt2(E_PT2) -! call ezfio_get_casscf_cipsi_energy(PT2) - double precision :: delta_E_istate, e_av - e_av = 0.d0 - do istate=1,N_states - e_av += state_average_weight(istate) * Ev(istate) - if(istate.gt.1)then - delta_E_istate = E_PT2(istate) - E_PT2(1) - write(*,'(A6,I2,A18,F16.10)')'state ',istate,' Delta E+PT2 = ',delta_E_istate - endif - write(*,'(A6,I2,A18,F16.10)')'state ',istate,' E + PT2 energy = ',E_PT2(istate) - write(*,'(A6,I2,A18,F16.10)')'state ',istate,' PT2 energy = ',PT2(istate) -! call write_double(6,E_PT2(istate),'E + PT2 energy = ') -! call write_double(6,PT2(istate),' PT2 = ') - enddo - call write_double(6,e_av,'State-average CAS-SCF energy bis = ') - call write_double(6,pt2_max,' PT2_MAX = ') +! if(act_mos_opt)then DOES NOT WORK +! call run_orb_opt_trust_v2 +! call run_stochastic_cipsi(Ev,PT2) ! 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 if (criterion_casscf == "e_pt2")then - 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(.not.small_active_space)then - if(adaptive_pt2_max)then - pt2_max = dabs(energy_improvement / (pt2_relative_error)) - pt2_max = min(pt2_max, pt2_max_before) - if(n_act_orb.ge.n_big_act_orb)then - pt2_max = max(pt2_max,pt2_min_casscf) - endif + if(.True.)then + 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 + pt2_max_before = pt2_max + + call write_time(6) + call write_int(6,iteration,'CAS-SCF iteration = ') + call write_double(6,energy,'State-average CAS-SCF energy = ') +!! if(n_states == 1)then +!! call ezfio_get_casscf_cipsi_energy_pt2(E_PT2) +!! call ezfio_get_casscf_cipsi_energy(PT2) + double precision :: delta_E_istate, e_av + e_av = 0.d0 + do istate=1,N_states + e_av += state_average_weight(istate) * Ev(istate) + if(istate.gt.1)then + delta_E_istate = E_PT2(istate) - E_PT2(1) + write(*,'(A6,I2,A18,F16.10)')'state ',istate,' Delta E+PT2 = ',delta_E_istate + endif + write(*,'(A6,I2,A18,F16.10)')'state ',istate,' E + PT2 energy = ',E_PT2(istate) + write(*,'(A6,I2,A18,F16.10)')'state ',istate,' PT2 energy = ',PT2(istate) +!! call write_double(6,E_PT2(istate),'E + PT2 energy = ') +!! call write_double(6,PT2(istate),' PT2 = ') + enddo + call write_double(6,e_av,'State-average CAS-SCF energy bis = ') + 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 if (criterion_casscf == "e_pt2")then + 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 - endif - print*,'' - call write_double(6,pt2_max, 'PT2_MAX for next iteration = ') - - mo_coef = NewOrbs - mo_occ = occnum - if(.not.converged)then - call save_mos - iteration += 1 - if(norm_grad_vec2.gt.0.01d0)then - N_det = N_states - else - N_det = max(N_det/8 ,N_states) - endif - psi_det = psi_det_sorted - psi_coef = psi_coef_sorted - read_wf = .True. - call clear_mo_map - SOFT_TOUCH mo_coef N_det psi_det psi_coef + ept2_before = E_PT2 if(.not.small_active_space)then if(adaptive_pt2_max)then - SOFT_TOUCH pt2_max + pt2_max = dabs(energy_improvement / (pt2_relative_error)) + pt2_max = min(pt2_max, pt2_max_before) + if(n_act_orb.ge.n_big_act_orb)then + pt2_max = max(pt2_max,pt2_min_casscf) + endif endif endif - if(iteration .gt. 3)then - state_following_casscf = state_following_casscf_cipsi_save - soft_touch state_following_casscf + print*,'' + call write_double(6,pt2_max, 'PT2_MAX for next iteration = ') + + mo_coef = NewOrbs + mo_occ = occnum + if(.not.converged)then + call save_mos + iteration += 1 + if(norm_grad_vec2.gt.0.01d0)then + N_det = N_states + else + N_det = max(N_det/8 ,N_states) + endif + psi_det = psi_det_sorted + psi_coef = psi_coef_sorted + read_wf = .True. + call clear_mo_map + SOFT_TOUCH mo_coef N_det psi_det psi_coef + if(.not.small_active_space)then + if(adaptive_pt2_max)then + SOFT_TOUCH pt2_max + endif + endif + if(iteration .gt. 3)then + state_following_casscf = state_following_casscf_cipsi_save + soft_touch state_following_casscf + endif endif endif - + enddo + if(.True.)then integer :: i print*,'Converged CASSCF ' print*,'--------------------------' @@ -153,6 +160,7 @@ subroutine run ! write(*,*)mcscf_fock_alpha_mo(i,i) enddo + endif end diff --git a/src/mo_optimization_utils/class.irp.f b/src/mo_optimization/class.irp.f similarity index 100% rename from src/mo_optimization_utils/class.irp.f rename to src/mo_optimization/class.irp.f