From f2fef4a6dc37770d8c8842986c9b6f3d3488c159 Mon Sep 17 00:00:00 2001 From: Yann Garniron Date: Thu, 7 Apr 2016 17:54:13 +0200 Subject: [PATCH] added pt2 to mrsc2/mrcepa0 --- plugins/MRCC_Utils/mrcc_general.irp.f | 61 -------- plugins/MRCC_Utils/mrcc_utils.irp.f | 49 +++++++ plugins/mrcepa0/mrcepa0.irp.f | 31 ++--- plugins/mrcepa0/mrcepa0_general.irp.f | 192 +++++++++++++++++++------- plugins/mrcepa0/mrsc2.irp.f | 32 ++--- plugins/mrcepa0/mrsc2sub.irp.f | 33 ++--- 6 files changed, 228 insertions(+), 170 deletions(-) diff --git a/plugins/MRCC_Utils/mrcc_general.irp.f b/plugins/MRCC_Utils/mrcc_general.irp.f index a0a4c895..d356e4b9 100644 --- a/plugins/MRCC_Utils/mrcc_general.irp.f +++ b/plugins/MRCC_Utils/mrcc_general.irp.f @@ -1,61 +1,3 @@ -subroutine run_mrcc - implicit none - call set_generators_bitmasks_as_holes_and_particles - call mrcc_iterations -end - -subroutine mrcc_iterations - implicit none - - integer :: i,j - - double precision :: E_new, E_old, delta_e - integer :: iteration,i_oscillations - double precision :: E_past(4), lambda - E_new = 0.d0 - delta_E = 1.d0 - iteration = 0 - j = 1 - i_oscillations = 0 - lambda = 1.d0 - do while (delta_E > 1.d-7) - iteration += 1 - print *, '===========================' - print *, 'MRCC Iteration', iteration - print *, '===========================' - print *, '' - E_old = sum(ci_energy_dressed) - print *, iteration, ci_energy_dressed(1) - call write_double(6,ci_energy_dressed(1),"MRCC energy") - call diagonalize_ci_dressed(lambda) - E_new = sum(ci_energy_dressed) - delta_E = dabs(E_new - E_old) -! if (E_new > E_old) then -! lambda = lambda * 0.7d0 -! else -! lambda = min(1.d0, lambda * 1.1d0) -! endif -! print *, 'energy lambda ', lambda -! E_past(j) = E_new -! j +=1 - call save_wavefunction - if (iteration > 0) then - exit - endif - print*,'------------' - print*,'VECTOR' - do i = 1, N_det_ref - print*,'' - print*,'psi_ref_coef(i,1) = ',psi_ref_coef(i,1) - print*,'delta_ii(i,1) = ',delta_ii(i,1) - enddo - print*,'------------' - enddo - call write_double(6,ci_energy_dressed(1),"Final MRCC energy") - call ezfio_set_mrcc_cassd_energy(ci_energy_dressed(1)) - call save_wavefunction - -end subroutine set_generators_bitmasks_as_holes_and_particles implicit none @@ -82,7 +24,4 @@ subroutine set_generators_bitmasks_as_holes_and_particles enddo enddo touch generators_bitmask - - - end diff --git a/plugins/MRCC_Utils/mrcc_utils.irp.f b/plugins/MRCC_Utils/mrcc_utils.irp.f index c4f5c570..8445a9c6 100644 --- a/plugins/MRCC_Utils/mrcc_utils.irp.f +++ b/plugins/MRCC_Utils/mrcc_utils.irp.f @@ -1,3 +1,52 @@ +! +! BEGIN_PROVIDER [ double precision, lambda_mrcc, (N_states,psi_det_size) ] +! &BEGIN_PROVIDER [ integer, lambda_mrcc_pt2, (0:psi_det_size) ] +! implicit none +! BEGIN_DOC +! cm/ or perturbative 1/Delta_E(m) +! END_DOC +! integer :: i,k +! double precision :: ihpsi_current(N_states) +! integer :: i_pert_count +! double precision :: hii, lambda_pert +! lambda_mrcc_pt2(:) = 0d0 +! i_pert_count = 0 +! lambda_mrcc = 0.d0 +! +! do i=1,N_det_non_ref +! call i_h_psi(psi_non_ref(1,1,i), psi_ref, psi_ref_coef, N_int, N_det_ref, & +! size(psi_ref_coef,1), N_states,ihpsi_current) +! call i_H_j(psi_non_ref(1,1,i),psi_non_ref(1,1,i),N_int,hii) +! do k=1,N_states +! if (ihpsi_current(k) == 0.d0) then +! ihpsi_current(k) = 1.d-32 +! endif +! lambda_mrcc(k,i) = psi_non_ref_coef(i,k)/ihpsi_current(k) +! if ( dabs(psi_non_ref_coef(i,k)*ihpsi_current(k)) < 1.d-5 ) then +! i_pert_count += 1 +! lambda_mrcc(k,i) = 0.d0 +! lambda_pert = 1.d0 / (psi_ref_energy_diagonalized(k)-hii) +! if((ihpsi_current(k) * lambda_pert) < 0.5d0 * psi_non_ref_coef_restart(i,k) ) then +! lambda_mrcc(k,i) = 0.d0 +! endif +! endif +! double precision, parameter :: x = 2.d0 +! if (lambda_mrcc(k,i) > x) then +! lambda_mrcc(k,i) = x +! else if (lambda_mrcc(k,i) < -x) then +! lambda_mrcc(k,i) = -x +! endif +! enddo +! enddo +! +! print*,'N_det_non_ref = ',N_det_non_ref +! print*,'Number of ignored determinants = ',i_pert_count +! print*,'psi_coef_ref_ratio = ',psi_ref_coef(2,1)/psi_ref_coef(1,1) +! print*,'lambda min/max = ',maxval(dabs(lambda_mrcc)), minval(dabs(lambda_mrcc)) +! +! END_PROVIDER + + BEGIN_PROVIDER [ double precision, lambda_mrcc, (N_states,psi_det_size) ] &BEGIN_PROVIDER [ integer, lambda_mrcc_pt2, (0:psi_det_size) ] implicit none diff --git a/plugins/mrcepa0/mrcepa0.irp.f b/plugins/mrcepa0/mrcepa0.irp.f index 7877cdda..9473361b 100644 --- a/plugins/mrcepa0/mrcepa0.irp.f +++ b/plugins/mrcepa0/mrcepa0.irp.f @@ -1,26 +1,19 @@ program mrcepa0 implicit none + double precision, allocatable :: energy(:) + allocate (energy(N_states)) + !mrmode : 1=mrcepa0, 2=mrsc2 add, 3=mrsc2 sub mrmode = 1 - if (.not.read_wf) then - print *, 'read_wf has to be true.' - stop 1 - endif + + read_wf = .True. + SOFT_TOUCH read_wf call print_cas_coefs - call run_mrcepa0 -end - -subroutine print_cas_coefs - implicit none - - integer :: i,j - print *, 'CAS' - print *, '===' - do i=1,N_det_cas - print *, psi_cas_coef(i,:) - call debug_det(psi_cas(1,1,i),N_int) - enddo - - call write_double(6,ci_energy(1),"Initial CI energy") + call set_generators_bitmasks_as_holes_and_particles + call run(N_states,energy) + if(do_pt2_end)then + call run_pt2(N_states,energy) + endif + deallocate(energy) end diff --git a/plugins/mrcepa0/mrcepa0_general.irp.f b/plugins/mrcepa0/mrcepa0_general.irp.f index 525b70df..b3390577 100644 --- a/plugins/mrcepa0/mrcepa0_general.irp.f +++ b/plugins/mrcepa0/mrcepa0_general.irp.f @@ -1,28 +1,33 @@ -subroutine run_mrcepa0 - implicit none - call set_generators_bitmasks_as_holes_and_particles - call mrcepa0_iterations -end + BEGIN_PROVIDER [ integer, mrmode ] END_PROVIDER -subroutine mrcepa0_iterations + +subroutine run(N_st,energy) implicit none - integer :: i,j + integer, intent(in) :: N_st + double precision, intent(out) :: energy(N_st) + + integer :: i double precision :: E_new, E_old, delta_e - integer :: iteration,i_oscillations + integer :: iteration double precision :: E_past(4), lambda + + integer :: n_it_mrcc_max + double precision :: thresh_mrcc + + thresh_mrcc = 1d-7 + n_it_mrcc_max = 10 + E_new = 0.d0 delta_E = 1.d0 iteration = 0 - j = 1 - i_oscillations = 0 lambda = 1.d0 - do while (delta_E > 1.d-7) + do while (delta_E > thresh_mrcc) iteration += 1 print *, '===========================' print *, 'MRCEPA0 Iteration', iteration @@ -33,57 +38,142 @@ subroutine mrcepa0_iterations call diagonalize_ci_dressed(lambda) E_new = sum(ci_energy_dressed) delta_E = dabs(E_new - E_old) -! if (E_new > E_old) then -! lambda = lambda * 0.7d0 -! else -! lambda = min(1.d0, lambda * 1.1d0) -! endif -! print *, 'energy lambda ', lambda -! E_past(j) = E_new -! j +=1 call save_wavefunction - if (iteration > 10) then + call ezfio_set_mrcc_cassd_energy(ci_energy_dressed(1)) + if (iteration > n_it_mrcc_max) then exit endif - print*,'------------' - print*,'VECTOR' - do i = 1, N_det_ref - print*,'' - print*,'psi_ref_coef(i,1) = ',psi_ref_coef(i,1) - print*,'delta_ii(i,1) = ',delta_ii(i,1) - enddo - print*,'------------' enddo call write_double(6,ci_energy_dressed(1),"Final MRCEPA0 energy") - call ezfio_set_mrcc_cassd_energy(ci_energy_dressed(1)) - call save_wavefunction + energy(:) = ci_energy_dressed(:) + end -subroutine set_generators_bitmasks_as_holes_and_particles - implicit none - integer :: i,k - do k = 1, N_generators_bitmask - do i = 1, N_int - ! Pure single part - generators_bitmask(i,1,1,k) = holes_operators(i,1) ! holes for pure single exc alpha - generators_bitmask(i,1,2,k) = particles_operators(i,1) ! particles for pure single exc alpha - generators_bitmask(i,2,1,k) = holes_operators(i,2) ! holes for pure single exc beta - generators_bitmask(i,2,2,k) = particles_operators(i,2) ! particles for pure single exc beta - ! Double excitation - generators_bitmask(i,1,3,k) = holes_operators(i,1) ! holes for first single exc alpha - generators_bitmask(i,1,4,k) = particles_operators(i,1) ! particles for first single exc alpha - generators_bitmask(i,2,3,k) = holes_operators(i,2) ! holes for first single exc beta - generators_bitmask(i,2,4,k) = particles_operators(i,2) ! particles for first single exc beta +subroutine run_pt2(N_st,energy) + implicit none + integer :: i,j,k + double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:) + integer, intent(in) :: N_st + double precision, intent(in) :: energy(N_st) + allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st)) + pt2 = 0.d0 + !if(lambda_mrcc_pt2(0) == 0) return + + print*,'Last iteration only to compute the PT2' + threshold_selectors = 1.d0 + threshold_generators = 0.999d0 + + N_det_generators = lambda_mrcc_pt2(0) + do i=1,N_det_generators + j = lambda_mrcc_pt2(i) + do k=1,N_int + psi_det_generators(k,1,i) = psi_non_ref(k,1,j) + psi_det_generators(k,2,i) = psi_non_ref(k,2,j) + enddo + do k=1,N_st + psi_coef_generators(i,k) = psi_non_ref_coef(j,k) + enddo + enddo + SOFT_TOUCH N_det_generators psi_det_generators psi_coef_generators ci_eigenvectors_dressed ci_eigenvectors_s2_dressed ci_electronic_energy_dressed + + + call H_apply_mrcc_PT2(pt2, norm_pert, H_pert_diag, N_st) + print *, 'Final step' + print *, 'N_det = ', N_det + print *, 'N_states = ', N_states + print *, 'PT2 = ', pt2 + print *, 'E = ', energy + print *, 'E+PT2 = ', energy+pt2 + print *, '-----' + - generators_bitmask(i,1,5,k) = holes_operators(i,1) ! holes for second single exc alpha - generators_bitmask(i,1,6,k) = particles_operators(i,1) ! particles for second single exc alpha - generators_bitmask(i,2,5,k) = holes_operators(i,2) ! holes for second single exc beta - generators_bitmask(i,2,6,k) = particles_operators(i,2) ! particles for second single exc beta + call ezfio_set_full_ci_energy_pt2(energy+pt2) + deallocate(pt2,norm_pert) +end + + +subroutine print_cas_coefs + implicit none + + integer :: i,j + print *, 'CAS' + print *, '===' + do i=1,N_det_cas + print *, psi_cas_coef(i,:) + call debug_det(psi_cas(1,1,i),N_int) enddo - enddo - touch generators_bitmask + call write_double(6,ci_energy(1),"Initial CI energy") + end + + + + + + + + + +! subroutine run_mrcepa0 +! implicit none +! call set_generators_bitmasks_as_holes_and_particles +! call mrcepa0_iterations +! end +! + +! +! subroutine mrcepa0_iterations +! implicit none +! +! integer :: i,j +! +! double precision :: E_new, E_old, delta_e +! integer :: iteration,i_oscillations +! double precision :: E_past(4), lambda +! E_new = 0.d0 +! delta_E = 1.d0 +! iteration = 0 +! j = 1 +! i_oscillations = 0 +! lambda = 1.d0 +! do while (delta_E > 1.d-7) +! iteration += 1 +! print *, '===========================' +! print *, 'MRCEPA0 Iteration', iteration +! print *, '===========================' +! print *, '' +! E_old = sum(ci_energy_dressed) +! call write_double(6,ci_energy_dressed(1),"MRCEPA0 energy") +! call diagonalize_ci_dressed(lambda) +! E_new = sum(ci_energy_dressed) +! delta_E = dabs(E_new - E_old) +! ! if (E_new > E_old) then +! ! lambda = lambda * 0.7d0 +! ! else +! ! lambda = min(1.d0, lambda * 1.1d0) +! ! endif +! ! print *, 'energy lambda ', lambda +! ! E_past(j) = E_new +! ! j +=1 +! call save_wavefunction +! if (iteration > 10) then +! exit +! endif +! print*,'------------' +! print*,'VECTOR' +! do i = 1, N_det_ref +! print*,'' +! print*,'psi_ref_coef(i,1) = ',psi_ref_coef(i,1) +! print*,'delta_ii(i,1) = ',delta_ii(i,1) +! enddo +! print*,'------------' +! enddo +! call write_double(6,ci_energy_dressed(1),"Final MRCEPA0 energy") +! call ezfio_set_mrcc_cassd_energy(ci_energy_dressed(1)) +! call save_wavefunction +! end + diff --git a/plugins/mrcepa0/mrsc2.irp.f b/plugins/mrcepa0/mrsc2.irp.f index eb34adee..d4e1b1d4 100644 --- a/plugins/mrcepa0/mrsc2.irp.f +++ b/plugins/mrcepa0/mrsc2.irp.f @@ -1,26 +1,20 @@ -program mrcepa0 +program mrsc2 implicit none + double precision, allocatable :: energy(:) + allocate (energy(N_states)) + !mrmode : 1=mrcepa0, 2=mrsc2 add, 3=mrsc2 sub mrmode = 2 - if (.not.read_wf) then - print *, 'read_wf has to be true.' - stop 1 - endif + + read_wf = .True. + SOFT_TOUCH read_wf call print_cas_coefs - call run_mrcepa0 + call set_generators_bitmasks_as_holes_and_particles + call run(N_states,energy) + if(do_pt2_end)then + call run_pt2(N_states,energy) + endif + deallocate(energy) end -subroutine print_cas_coefs - implicit none - - integer :: i,j - print *, 'CAS' - print *, '===' - do i=1,N_det_cas - print *, psi_cas_coef(i,:) - call debug_det(psi_cas(1,1,i),N_int) - enddo - - call write_double(6,ci_energy(1),"Initial CI energy") -end diff --git a/plugins/mrcepa0/mrsc2sub.irp.f b/plugins/mrcepa0/mrsc2sub.irp.f index 524bbbd7..07a07c83 100644 --- a/plugins/mrcepa0/mrsc2sub.irp.f +++ b/plugins/mrcepa0/mrsc2sub.irp.f @@ -1,26 +1,19 @@ -program mrcepa0 +program mrsc2sub implicit none + double precision, allocatable :: energy(:) + allocate (energy(N_states)) + !mrmode : 1=mrcepa0, 2=mrsc2 add, 3=mrsc2 sub mrmode = 3 - if (.not.read_wf) then - print *, 'read_wf has to be true.' - stop 1 - endif + + read_wf = .True. + SOFT_TOUCH read_wf call print_cas_coefs - call run_mrcepa0 -end - -subroutine print_cas_coefs - implicit none - - integer :: i,j - print *, 'CAS' - print *, '===' - do i=1,N_det_cas - print *, psi_cas_coef(i,:) - call debug_det(psi_cas(1,1,i),N_int) - enddo - - call write_double(6,ci_energy(1),"Initial CI energy") + call set_generators_bitmasks_as_holes_and_particles + call run(N_states,energy) + if(do_pt2_end)then + call run_pt2(N_states,energy) + endif + deallocate(energy) end