From 8718f5fd35ded4a420a34fce9d18108d0ba962d3 Mon Sep 17 00:00:00 2001 From: Manu Date: Mon, 6 Oct 2014 15:49:16 +0200 Subject: [PATCH 1/3] minor changes --- src/CISD_SC2_selected/H_apply.irp.f | 2 +- .../cisd_sc2_selection.irp.f | 14 +++- src/Dets/determinants.ezfio_config | 1 + src/Dets/diagonalize_CI.irp.f | 21 ++++-- src/Dets/diagonalize_CI_SC2.irp.f | 6 +- src/Dets/options.irp.f | 1 + src/Dets/s2.irp.f | 15 ++++ src/Perturbation/pert_sc2.irp.f | 74 ++++++++++++++++++- src/Utils/LinearAlgebra.irp.f | 66 +++++++++++++++++ src/Utils/one_e_integration.irp.f | 8 +- src/Utils/util.irp.f | 8 +- 11 files changed, 194 insertions(+), 22 deletions(-) diff --git a/src/CISD_SC2_selected/H_apply.irp.f b/src/CISD_SC2_selected/H_apply.irp.f index c40d6d4a..85ecae95 100644 --- a/src/CISD_SC2_selected/H_apply.irp.f +++ b/src/CISD_SC2_selected/H_apply.irp.f @@ -4,7 +4,7 @@ from generate_h_apply import * from perturbation import perturbations s = H_apply("PT2",SingleRef=True) -s.set_perturbation("epstein_nesbet_sc2_projected") +s.set_perturbation("epstein_nesbet_sc2_no_projected") print s s = H_apply("PT2_en_sc2",SingleRef=True) diff --git a/src/CISD_SC2_selected/cisd_sc2_selection.irp.f b/src/CISD_SC2_selected/cisd_sc2_selection.irp.f index acce82f8..b8770155 100644 --- a/src/CISD_SC2_selected/cisd_sc2_selection.irp.f +++ b/src/CISD_SC2_selected/cisd_sc2_selection.irp.f @@ -12,8 +12,9 @@ program cisd_sc2_selected pt2 = 1.d0 perturbation = "epstein_nesbet_sc2_projected" + E_old(1) = HF_energy - davidson_threshold = 1.d-8 + davidson_threshold = 1.d-10 if (N_det > n_det_max_cisd_sc2) then call diagonalize_CI_SC2 call save_wavefunction @@ -31,6 +32,8 @@ program cisd_sc2_selected print *, '-----' endif + integer :: i_count + i_count = 0 do while (N_det < n_det_max_cisd_sc2.and.maxval(abs(pt2(1:N_st))) > pt2_max) print*,'----' print*,'' @@ -49,6 +52,13 @@ program cisd_sc2_selected E_old(i) = CI_SC2_energy(i) enddo ! print *, 'E corr = ', (E_old(1)) - HF_energy + if(dabs(E_old(i) - CI_SC2_energy(i) ).le.1.d-12)then + i_count += 1 + selection_criterion_factor = selection_criterion_factor * 0.5d0 + if(i_count > 5)then + exit + endif + endif if (abort_all) then exit endif @@ -81,7 +91,7 @@ program cisd_sc2_selected print *, 'PT2(SC2) = ', pt2(i) print *, 'E(SC2) = ', CI_SC2_energy(i) print *, 'E_before(SC2)+PT2(SC2) = ', CI_SC2_energy(i)+pt2(i) - print *, 'E_before(SC2)+PT2(SC2)_new = ', CI_SC2_energy(i)+pt2(i)*(1.d0+norm_pert) + print *, 'E_before(SC2)+PT2(SC2)_new = ', CI_SC2_energy(i)+pt2(i)* (1.d0 + norm_pert) - H_pert_diag(i) print*,'greater coeficient of the state : ',dabs(psi_coef(imax,i)) call get_excitation_degree(ref_bitmask,psi_det(1,1,imax),degree,N_int) diff --git a/src/Dets/determinants.ezfio_config b/src/Dets/determinants.ezfio_config index 47fe30e8..ff97ca36 100644 --- a/src/Dets/determinants.ezfio_config +++ b/src/Dets/determinants.ezfio_config @@ -13,4 +13,5 @@ determinants det_occ integer (electrons_elec_alpha_num,determinants_det_num,2) det_coef double precision (determinants_det_num) read_wf logical + expected_s2 double precision diff --git a/src/Dets/diagonalize_CI.irp.f b/src/Dets/diagonalize_CI.irp.f index f9cc3ef0..1329b581 100644 --- a/src/Dets/diagonalize_CI.irp.f +++ b/src/Dets/diagonalize_CI.irp.f @@ -59,12 +59,21 @@ END_PROVIDER call lapack_diag(eigenvalues,eigenvectors, & H_matrix_all_dets,size(H_matrix_all_dets,1),N_det) CI_electronic_energy(:) = 0.d0 - do j=1,min(N_states,N_det) - do i=1,N_det - CI_eigenvectors(i,j) = eigenvectors(i,j) - enddo - CI_electronic_energy(j) = eigenvalues(j) - enddo + integer :: i_state + double precision :: s2 + j=0 + i_state = 0 +! do while(i_state.lt.N_states) + j+=1 +! call get_s2_u0(psi_det,eigenvectors(1,j),N_det,N_det,s2) +! if(dabs(s2-expected_s2).le.0.1d0)then + i_state += 1 + do i=1,N_det + CI_eigenvectors(i,i_state) = eigenvectors(i,j) + enddo + CI_electronic_energy(i_state) = eigenvalues(j) +! endif +! enddo deallocate(eigenvectors,eigenvalues) endif diff --git a/src/Dets/diagonalize_CI_SC2.irp.f b/src/Dets/diagonalize_CI_SC2.irp.f index 987586b2..412dc2ae 100644 --- a/src/Dets/diagonalize_CI_SC2.irp.f +++ b/src/Dets/diagonalize_CI_SC2.irp.f @@ -20,7 +20,7 @@ END_PROVIDER BEGIN_DOC ! convergence of the correlation energy of SC2 iterations END_DOC - threshold_convergence_SC2 = 1.d-8 + threshold_convergence_SC2 = 1.d-10 END_PROVIDER BEGIN_PROVIDER [ double precision, CI_SC2_electronic_energy, (N_states) ] @@ -33,8 +33,8 @@ END_PROVIDER do j=1,N_states do i=1,N_det -! CI_SC2_eigenvectors(i,j) = psi_coef(i,j) - CI_SC2_eigenvectors(i,j) = CI_eigenvectors(i,j) + CI_SC2_eigenvectors(i,j) = psi_coef(i,j) +! CI_SC2_eigenvectors(i,j) = CI_eigenvectors(i,j) enddo CI_SC2_electronic_energy(j) = CI_electronic_energy(j) enddo diff --git a/src/Dets/options.irp.f b/src/Dets/options.irp.f index 2a8ffe47..df2ba1b5 100644 --- a/src/Dets/options.irp.f +++ b/src/Dets/options.irp.f @@ -20,6 +20,7 @@ T.set_doc ( "If true, read the wave function from the EZFIO file" ) T.set_ezfio_name( "read_wf" ) T.set_output ( "output_dets" ) print T + END_SHELL diff --git a/src/Dets/s2.irp.f b/src/Dets/s2.irp.f index b218e7cb..6993675c 100644 --- a/src/Dets/s2.irp.f +++ b/src/Dets/s2.irp.f @@ -42,6 +42,21 @@ BEGIN_PROVIDER [ double precision, S_z ] END_PROVIDER +BEGIN_PROVIDER [ double precision, expected_s2] + implicit none + PROVIDE ezfio_filename + logical :: has_expected_s2 + + call ezfio_has_determinants_expected_s2(has_expected_s2) + if (has_expected_s2) then + call ezfio_get_determinants_expected_s2(expected_s2) + else + expected_s2 = elec_alpha_num - elec_beta_num + 0.5d0 * ((elec_alpha_num - elec_beta_num)**2*0.5d0 - (elec_alpha_num-elec_beta_num)) + call ezfio_set_determinants_expected_s2(expected_s2) + endif + +END_PROVIDER + subroutine get_s2_u0(psi_keys_tmp,psi_coefs_tmp,n,nmax,s2) implicit none diff --git a/src/Perturbation/pert_sc2.irp.f b/src/Perturbation/pert_sc2.irp.f index e5094a2a..f0fc8a0f 100644 --- a/src/Perturbation/pert_sc2.irp.f +++ b/src/Perturbation/pert_sc2.irp.f @@ -72,8 +72,8 @@ subroutine pt2_epstein_nesbet_SC2_projected(det_pert,c_pert,e_2_pert,H_pert_diag enddo if(degree==4)then ! - H_pert_diag(1) = e_2_pert(1) - e_2_pert_fonda = H_pert_diag(1) + e_2_pert_fonda = e_2_pert(1) + H_pert_diag(1) = e_2_pert(1) * c_pert(1) * c_pert(1) do i = 1, N_st do j = 1, idx_repeat(0) e_2_pert(i) += e_2_pert_fonda * psi_selectors_coef(idx_repeat(j),i) * psi_selectors_coef(idx_repeat(j),i) @@ -83,6 +83,76 @@ subroutine pt2_epstein_nesbet_SC2_projected(det_pert,c_pert,e_2_pert,H_pert_diag end + +subroutine pt2_epstein_nesbet_SC2_no_projected(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st) + use bitmasks + implicit none + integer, intent(in) :: Nint,ndet,N_st + integer(bit_kind), intent(in) :: det_pert(Nint,2) + double precision , intent(out) :: c_pert(N_st),e_2_pert(N_st),H_pert_diag(N_st) + double precision :: i_H_psi_array(N_st) + integer :: idx_repeat(0:ndet) + + BEGIN_DOC + ! compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution + ! + ! for the various N_st states, + ! + ! but with the correction in the denominator + ! + ! comming from the interaction of that determinant with all the others determinants + ! + ! that can be repeated by repeating all the double excitations + ! + ! : you repeat all the correlation energy already taken into account in CI_electronic_energy(1) + ! + ! that could be repeated to this determinant. + ! + ! In addition, for the perturbative energetic contribution you have the standard second order + ! + ! e_2_pert = ^2/(Delta_E) + ! + ! and also the purely projected contribution + ! + ! H_pert_diag = c_pert + END_DOC + + integer :: i,j,degree,l + double precision :: diag_H_mat_elem,accu_e_corr,hij,h0j,h,delta_E + double precision :: repeat_all_e_corr,accu_e_corr_tmp,e_2_pert_fonda + + ASSERT (Nint == N_int) + ASSERT (Nint > 0) + + call i_H_psi_SC2(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array,idx_repeat) + accu_e_corr = 0.d0 + !$IVDEP + do i = 1, idx_repeat(0) + accu_e_corr = accu_e_corr + E_corr_per_selectors(idx_repeat(i)) + enddo + h = diag_H_mat_elem(det_pert,Nint) + accu_e_corr + delta_E = 1.d0/(CI_SC2_electronic_energy(1) - h) + + + c_pert(1) = i_H_psi_array(1) /(CI_SC2_electronic_energy(1) - h) + e_2_pert(1) = i_H_psi_array(1) * c_pert(1) + + do i =2,N_st + H_pert_diag(i) = h + if (dabs(CI_SC2_electronic_energy(i) - h) > 1.d-6) then + c_pert(i) = i_H_psi_array(i) / (-dabs(CI_SC2_electronic_energy(i) - h)) + e_2_pert(i) = (c_pert(i) * i_H_psi_array(i)) + else + c_pert(i) = i_H_psi_array(i) + e_2_pert(i) = -dabs(i_H_psi_array(i)) + endif + enddo +end + + + + + double precision function repeat_all_e_corr(key_in) implicit none integer(bit_kind), intent(in) :: key_in(N_int,2) diff --git a/src/Utils/LinearAlgebra.irp.f b/src/Utils/LinearAlgebra.irp.f index 9f2336db..b29654b1 100644 --- a/src/Utils/LinearAlgebra.irp.f +++ b/src/Utils/LinearAlgebra.irp.f @@ -307,6 +307,72 @@ subroutine lapack_diag(eigvalues,eigvectors,H,nmax,n) deallocate(A,eigenvalues) end +subroutine lapack_diag_s2(eigvalues,eigvectors,H,nmax,n) + implicit none + BEGIN_DOC + ! Diagonalize matrix H + ! + ! H is untouched between input and ouptut + ! + ! eigevalues(i) = ith lowest eigenvalue of the H matrix + ! + ! eigvectors(i,j) = where i is the basis function and psi_j is the j th eigenvector + ! + END_DOC + integer, intent(in) :: n,nmax + double precision, intent(out) :: eigvectors(nmax,n) + double precision, intent(out) :: eigvalues(n) + double precision, intent(in) :: H(nmax,n) + double precision,allocatable :: eigenvalues(:) + double precision,allocatable :: work(:) + double precision,allocatable :: A(:,:) + integer :: lwork, info, i,j,l,k, liwork + + allocate(A(nmax,n),eigenvalues(n)) +! print*,'Diagonalization by jacobi' +! print*,'n = ',n + + A=H + lwork = 2*n*n + 6*n+ 1 + allocate (work(lwork)) + + lwork = -1 + call DSYEV( 'V', 'U', n, A, nmax, eigenvalues, work, lwork, & + info ) + if (info < 0) then + print *, irp_here, ': DSYEV: the ',-info,'-th argument had an illegal value' + stop 2 + endif + lwork = int( work( 1 ) ) + deallocate (work) + + allocate (work(lwork)) + call DSYEV( 'V', 'U', n, A, nmax, eigenvalues, work, lwork, & + info ) + deallocate(work) + + if (info < 0) then + print *, irp_here, ': DSYEV: the ',-info,'-th argument had an illegal value' + stop 2 + else if( info > 0 ) then + write(*,*)'DSYEV Failed' + stop 1 + end if + + eigvectors = 0.d0 + eigvalues = 0.d0 + do j = 1, n + eigvalues(j) = eigenvalues(j) + do i = 1, n + eigvectors(i,j) = A(i,j) + enddo + enddo + deallocate(A,eigenvalues) +end + + + + subroutine lapack_partial_diag(eigvalues,eigvectors,H,nmax,n,n_st) implicit none BEGIN_DOC diff --git a/src/Utils/one_e_integration.irp.f b/src/Utils/one_e_integration.irp.f index d7e52d8e..273b4f48 100644 --- a/src/Utils/one_e_integration.irp.f +++ b/src/Utils/one_e_integration.irp.f @@ -15,10 +15,10 @@ double precision function overlap_gaussian_x(A_center,B_center,alpha,beta,power_ call give_explicit_poly_and_gaussian_x(P_new,P_center,p,fact_p,iorder_p,alpha,& beta,power_A,power_B,A_center,B_center,dim) - if(fact_p.lt.0.000001d0)then - overlap_gaussian_x = 0.d0 - return - endif +! if(fact_p.lt.0.000001d0)then +! overlap_gaussian_x = 0.d0 +! return +! endif overlap_gaussian_x = 0.d0 integer :: i diff --git a/src/Utils/util.irp.f b/src/Utils/util.irp.f index 7d81d49b..ab7ce9db 100644 --- a/src/Utils/util.irp.f +++ b/src/Utils/util.irp.f @@ -49,17 +49,17 @@ double precision function binom_func(i,j) end - BEGIN_PROVIDER [ double precision, binom, (0:20,0:20) ] -&BEGIN_PROVIDER [ double precision, binom_transp, (0:20,0:20) ] + BEGIN_PROVIDER [ double precision, binom, (0:40,0:40) ] +&BEGIN_PROVIDER [ double precision, binom_transp, (0:40,0:40) ] implicit none BEGIN_DOC ! Binomial coefficients END_DOC integer :: k,l double precision :: fact, f - do k=0,20 + do k=0,40 f = fact(k) - do l=0,20 + do l=0,40 binom(k,l) = f/(fact(l)*fact(k-l)) binom_transp(l,k) = binom(k,l) enddo From 5e7644e1664c3a2a7eb97d6bced0de7fb589babe Mon Sep 17 00:00:00 2001 From: Manu Date: Tue, 7 Oct 2014 16:02:42 +0200 Subject: [PATCH 2/3] Bug in setup_environment corrected --- setup_environment.sh | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/setup_environment.sh b/setup_environment.sh index 5b2d8342..ae0471c9 100755 --- a/setup_environment.sh +++ b/setup_environment.sh @@ -45,8 +45,8 @@ make -C ocaml qp_create_ezfio_from_xyz.native if [[ $? -ne 0 ]] then scripts/fetch_from_web.py "https://raw.github.com/hcarty/ocamlbrew/master/ocamlbrew-install" ocamlbrew-install.sh - cat < ocamlbrew-install.sh | env OCAMLBREW_FLAGS="-r" bash > ocaml_install.log - grep "source " install.log | grep "etc/ocamlbrew.bashrc" >> quantum_package.rc + cat < ocamlbrew-install.sh | env OCAMLBREW_FLAGS="-r" bash | tee ocaml_install.log + grep "source " ocaml_install.log | grep "etc/ocamlbrew.bashrc" >> quantum_package.rc source quantum_package.rc echo Y | opam install core async fi From 590cc9574b315c2d3598d6748d92bbf93a5f8a59 Mon Sep 17 00:00:00 2001 From: Manu Date: Tue, 7 Oct 2014 16:23:19 +0200 Subject: [PATCH 3/3] Working on S^2 --- data/ezfio_defaults | 3 ++- src/Dets/SC2.irp.f | 2 +- src/Dets/davidson.irp.f | 2 +- src/Dets/determinants.ezfio_config | 1 + src/Dets/determinants.irp.f | 21 ++++++-------------- src/Dets/diagonalize_CI.irp.f | 31 ++++++++++++++++-------------- src/Dets/diagonalize_CI_SC2.irp.f | 18 ++++++++--------- src/Dets/options.irp.f | 10 ++++++++++ src/Dets/s2.irp.f | 2 +- 9 files changed, 48 insertions(+), 42 deletions(-) diff --git a/data/ezfio_defaults b/data/ezfio_defaults index 1de27eff..be7e7605 100644 --- a/data/ezfio_defaults +++ b/data/ezfio_defaults @@ -17,8 +17,9 @@ cis_dressed determinants n_states 1 + n_states_diag determinants_n_states n_det_max_jacobi 5000 - threshold_generators 0.999 + threshold_generators 0.995 threshold_selectors 0.999 read_wf False diff --git a/src/Dets/SC2.irp.f b/src/Dets/SC2.irp.f index aa2f49d1..921cf896 100644 --- a/src/Dets/SC2.irp.f +++ b/src/Dets/SC2.irp.f @@ -175,7 +175,7 @@ subroutine CISD_SC2(dets_in,u_in,energies,dim_in,sze,N_st,Nint,convergence) enddo call lapack_diag(eigenvalues,eigenvectors, & H_matrix_tmp,size(H_matrix_all_dets,1),sze) - do j=1,min(N_states,sze) + do j=1,min(N_states_diag,sze) do i=1,sze u_in(i,j) = eigenvectors(i,j) enddo diff --git a/src/Dets/davidson.irp.f b/src/Dets/davidson.irp.f index a79e9123..1efb70ff 100644 --- a/src/Dets/davidson.irp.f +++ b/src/Dets/davidson.irp.f @@ -12,7 +12,7 @@ BEGIN_PROVIDER [ integer, davidson_sze_max ] ! Max number of Davidson sizes END_DOC ASSERT (davidson_sze_max <= davidson_iter_max) - davidson_sze_max = 8*N_states + davidson_sze_max = 8*N_states_diag END_PROVIDER subroutine davidson_diag(dets_in,u_in,energies,dim_in,sze,N_st,Nint,iunit) diff --git a/src/Dets/determinants.ezfio_config b/src/Dets/determinants.ezfio_config index ff97ca36..4641db0f 100644 --- a/src/Dets/determinants.ezfio_config +++ b/src/Dets/determinants.ezfio_config @@ -4,6 +4,7 @@ determinants mo_label character*(64) n_det integer n_states integer + n_states_diag integer psi_coef double precision (determinants_n_det,determinants_n_states) psi_det integer*8 (determinants_n_int*determinants_bit_kind/8,2,determinants_n_det) n_det_max_jacobi integer diff --git a/src/Dets/determinants.irp.f b/src/Dets/determinants.irp.f index f7ec2dfc..1b52fafd 100644 --- a/src/Dets/determinants.irp.f +++ b/src/Dets/determinants.irp.f @@ -150,7 +150,7 @@ subroutine read_dets(det,Nint,Ndet) end -BEGIN_PROVIDER [ double precision, psi_coef, (psi_det_size,N_states) ] +BEGIN_PROVIDER [ double precision, psi_coef, (psi_det_size,N_states_diag) ] implicit none BEGIN_DOC ! The wave function coefficients. Initialized with Hartree-Fock if the EZFIO file @@ -161,6 +161,11 @@ BEGIN_PROVIDER [ double precision, psi_coef, (psi_det_size,N_states) ] logical :: exists double precision, allocatable :: psi_coef_read(:,:) character*(64) :: label + + psi_coef = 0.d0 + do i=1,N_states_diag + psi_coef(i,i) = 1.d0 + enddo if (read_wf) then call ezfio_has_determinants_psi_coef(exists) @@ -183,22 +188,8 @@ BEGIN_PROVIDER [ double precision, psi_coef, (psi_det_size,N_states) ] enddo deallocate(psi_coef_read) - else - - psi_coef = 0.d0 - do i=1,N_states - psi_coef(i,i) = 1.d0 - enddo - endif - else - - psi_coef = 0.d0 - do i=1,N_states - psi_coef(i,i) = 1.d0 - enddo - endif diff --git a/src/Dets/diagonalize_CI.irp.f b/src/Dets/diagonalize_CI.irp.f index 1329b581..f3a352f3 100644 --- a/src/Dets/diagonalize_CI.irp.f +++ b/src/Dets/diagonalize_CI.irp.f @@ -9,13 +9,13 @@ BEGIN_PROVIDER [ character*(64), diag_algorithm ] diag_algorithm = "Lapack" endif - if (N_det < N_states) then + if (N_det < N_states_diag) then diag_algorithm = "Lapack" endif END_PROVIDER -BEGIN_PROVIDER [ double precision, CI_energy, (N_states) ] +BEGIN_PROVIDER [ double precision, CI_energy, (N_states_diag) ] implicit none BEGIN_DOC ! N_states lowest eigenvalues of the CI matrix @@ -24,23 +24,25 @@ BEGIN_PROVIDER [ double precision, CI_energy, (N_states) ] integer :: j character*(8) :: st call write_time(output_Dets) - do j=1,N_states + do j=1,N_states_diag CI_energy(j) = CI_electronic_energy(j) + nuclear_repulsion write(st,'(I4)') j call write_double(output_Dets,CI_energy(j),'Energy of state '//trim(st)) + call write_double(output_Dets,CI_eigenvectors_s2(j),'S^2 of state '//trim(st)) enddo END_PROVIDER - BEGIN_PROVIDER [ double precision, CI_electronic_energy, (N_states) ] -&BEGIN_PROVIDER [ double precision, CI_eigenvectors, (N_det,N_states) ] + BEGIN_PROVIDER [ double precision, CI_electronic_energy, (N_states_diag) ] +&BEGIN_PROVIDER [ double precision, CI_eigenvectors, (N_det,N_states_diag) ] +&BEGIN_PROVIDER [ double precision, CI_eigenvectors_s2, (N_states_diag) ] implicit none BEGIN_DOC ! Eigenvectors/values of the CI matrix END_DOC integer :: i,j - do j=1,N_states + do j=1,N_states_diag do i=1,N_det CI_eigenvectors(i,j) = psi_coef(i,j) enddo @@ -49,7 +51,7 @@ END_PROVIDER if (diag_algorithm == "Davidson") then call davidson_diag(psi_det,CI_eigenvectors,CI_electronic_energy, & - size(CI_eigenvectors,1),N_det,N_states,N_int,output_Dets) + size(CI_eigenvectors,1),N_det,N_states_diag,N_int,output_Dets) else if (diag_algorithm == "Lapack") then @@ -63,17 +65,18 @@ END_PROVIDER double precision :: s2 j=0 i_state = 0 -! do while(i_state.lt.N_states) + do while(i_state.lt.min(N_states_diag,N_det)) j+=1 -! call get_s2_u0(psi_det,eigenvectors(1,j),N_det,N_det,s2) -! if(dabs(s2-expected_s2).le.0.1d0)then + call get_s2_u0(psi_det,eigenvectors(1,j),N_det,N_det,s2) + if(dabs(s2-expected_s2).le.0.1d0)then i_state += 1 do i=1,N_det CI_eigenvectors(i,i_state) = eigenvectors(i,j) enddo CI_electronic_energy(i_state) = eigenvalues(j) -! endif -! enddo + CI_eigenvectors_s2(i_state) = s2 + endif + enddo deallocate(eigenvectors,eigenvalues) endif @@ -86,10 +89,10 @@ subroutine diagonalize_CI ! eigenstates of the CI matrix END_DOC integer :: i,j - do j=1,N_states + do j=1,N_states_diag do i=1,N_det psi_coef(i,j) = CI_eigenvectors(i,j) enddo enddo - SOFT_TOUCH psi_coef CI_electronic_energy CI_energy CI_eigenvectors + SOFT_TOUCH psi_coef CI_electronic_energy CI_energy CI_eigenvectors CI_eigenvectors_s2 end diff --git a/src/Dets/diagonalize_CI_SC2.irp.f b/src/Dets/diagonalize_CI_SC2.irp.f index 412dc2ae..02252826 100644 --- a/src/Dets/diagonalize_CI_SC2.irp.f +++ b/src/Dets/diagonalize_CI_SC2.irp.f @@ -1,13 +1,13 @@ -BEGIN_PROVIDER [ double precision, CI_SC2_energy, (N_states) ] +BEGIN_PROVIDER [ double precision, CI_SC2_energy, (N_states_diag) ] implicit none BEGIN_DOC - ! N_states lowest eigenvalues of the CI matrix + ! N_states_diag lowest eigenvalues of the CI matrix END_DOC integer :: j character*(8) :: st call write_time(output_Dets) - do j=1,N_states + do j=1,N_states_diag CI_SC2_energy(j) = CI_SC2_electronic_energy(j) + nuclear_repulsion write(st,'(I4)') j call write_double(output_Dets,CI_SC2_energy(j),'Energy of state '//trim(st)) @@ -23,15 +23,15 @@ END_PROVIDER threshold_convergence_SC2 = 1.d-10 END_PROVIDER - BEGIN_PROVIDER [ double precision, CI_SC2_electronic_energy, (N_states) ] -&BEGIN_PROVIDER [ double precision, CI_SC2_eigenvectors, (N_det,N_states) ] + BEGIN_PROVIDER [ double precision, CI_SC2_electronic_energy, (N_states_diag) ] +&BEGIN_PROVIDER [ double precision, CI_SC2_eigenvectors, (N_det,N_states_diag) ] implicit none BEGIN_DOC ! Eigenvectors/values of the CI matrix END_DOC integer :: i,j - do j=1,N_states + do j=1,N_states_diag do i=1,N_det CI_SC2_eigenvectors(i,j) = psi_coef(i,j) ! CI_SC2_eigenvectors(i,j) = CI_eigenvectors(i,j) @@ -41,17 +41,17 @@ END_PROVIDER double precision :: convergence call CISD_SC2(psi_det,CI_SC2_eigenvectors,CI_SC2_electronic_energy, & - size(CI_SC2_eigenvectors,1),N_det,N_states,N_int,threshold_convergence_SC2) + size(CI_SC2_eigenvectors,1),N_det,N_states_diag,N_int,threshold_convergence_SC2) END_PROVIDER subroutine diagonalize_CI_SC2 implicit none BEGIN_DOC -! Replace the coefficients of the CI states by the coefficients of the +! Replace the coefficients of the CI states_diag by the coefficients of the ! eigenstates of the CI matrix END_DOC integer :: i,j - do j=1,N_states + do j=1,N_states_diag do i=1,N_det psi_coef(i,j) = CI_SC2_eigenvectors(i,j) enddo diff --git a/src/Dets/options.irp.f b/src/Dets/options.irp.f index df2ba1b5..c918e2c2 100644 --- a/src/Dets/options.irp.f +++ b/src/Dets/options.irp.f @@ -9,6 +9,16 @@ T.set_ezfio_name( "N_states" ) T.set_output ( "output_dets" ) print T + +T.set_type ( "integer" ) +T.set_name ( "N_states_diag" ) +T.set_doc ( "Number of states to consider for the diagonalization " ) +T.set_ezfio_dir ( "determinants" ) +T.set_ezfio_name( "N_states_diag" ) +T.set_output ( "output_dets" ) +print T + + T.set_name ( "N_det_max_jacobi" ) T.set_doc ( "Maximum number of determinants diagonalized by Jacobi" ) T.set_ezfio_name( "N_det_max_jacobi" ) diff --git a/src/Dets/s2.irp.f b/src/Dets/s2.irp.f index 6993675c..18c3dfaf 100644 --- a/src/Dets/s2.irp.f +++ b/src/Dets/s2.irp.f @@ -52,7 +52,7 @@ BEGIN_PROVIDER [ double precision, expected_s2] call ezfio_get_determinants_expected_s2(expected_s2) else expected_s2 = elec_alpha_num - elec_beta_num + 0.5d0 * ((elec_alpha_num - elec_beta_num)**2*0.5d0 - (elec_alpha_num-elec_beta_num)) - call ezfio_set_determinants_expected_s2(expected_s2) +! call ezfio_set_determinants_expected_s2(expected_s2) endif END_PROVIDER