From ca4f8ebdca35acdb0fae181878b2eff5cc222144 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 22 Feb 2019 17:59:19 +0100 Subject: [PATCH] Fixed energies of non-expected s2 (#9) * Moved diag_algorithm in Davdison --- configure | 2 - ocaml/.gitignore | 2 + src/davidson/EZFIO.cfg | 6 ++ .../diagonalization_hs2_dressed.irp.f | 55 +++++++++++++++---- src/davidson/diagonalize_ci.irp.f | 39 +++++++++---- src/determinants/EZFIO.cfg | 6 -- src/determinants/determinants.irp.f | 17 ------ src/fci/40.fci.bats | 2 +- 8 files changed, 81 insertions(+), 48 deletions(-) diff --git a/configure b/configure index 767900cd..f8c59957 100755 --- a/configure +++ b/configure @@ -11,8 +11,6 @@ echo "QP_ROOT="$QP_ROOT # Check if the module to create new DFT functionals exists or not if [[ ! -d ${QP_ROOT}/src/new_functionals ]] ; then - echo "${QP_ROOT}/src/new_functionals does NOT exists !!" - echo "Creating a new one ..." cp -r ${QP_ROOT}/scripts/functionals/do_not_touch_func ${QP_ROOT}/src/new_functionals fi diff --git a/ocaml/.gitignore b/ocaml/.gitignore index 50a9344d..74e83c02 100644 --- a/ocaml/.gitignore +++ b/ocaml/.gitignore @@ -9,6 +9,7 @@ Input_ao_two_e_erf_ints.ml Input_ao_two_e_ints.ml Input_auto_generated.ml Input_becke_numerical_grid.ml +Input_champ.ml Input_davidson.ml Input_density_for_dft.ml Input_determinants.ml @@ -21,6 +22,7 @@ Input_nuclei.ml Input_perturbation.ml Input_pseudo.ml Input_scf_utils.ml +Input_variance.ml qp_create_ezfio qp_create_ezfio.native qp_edit diff --git a/src/davidson/EZFIO.cfg b/src/davidson/EZFIO.cfg index e1ad58ae..67d7786a 100644 --- a/src/davidson/EZFIO.cfg +++ b/src/davidson/EZFIO.cfg @@ -40,3 +40,9 @@ doc: If |true|, use filter out all vectors with bad |S^2| values default: True interface: ezfio,provider,ocaml +[n_det_max_full] +type: Det_number_max +doc: Maximum number of determinants where |H| is fully diagonalized +interface: ezfio,provider,ocaml +default: 1000 + diff --git a/src/davidson/diagonalization_hs2_dressed.irp.f b/src/davidson/diagonalization_hs2_dressed.irp.f index a9f07962..f70fe78b 100644 --- a/src/davidson/diagonalization_hs2_dressed.irp.f +++ b/src/davidson/diagonalization_hs2_dressed.irp.f @@ -1,3 +1,20 @@ +BEGIN_PROVIDER [ character*(64), diag_algorithm ] + implicit none + BEGIN_DOC + ! Diagonalization algorithm (Davidson or Lapack) + END_DOC + if (N_det > N_det_max_full) then + diag_algorithm = "Davidson" + else + diag_algorithm = "Lapack" + endif + + if (N_det < N_states) then + diag_algorithm = "Lapack" + endif +END_PROVIDER + + BEGIN_PROVIDER [ integer, dressed_column_idx, (N_states) ] implicit none BEGIN_DOC @@ -114,7 +131,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_ integer :: k_pairs, kl integer :: iter2, itertot - double precision, allocatable :: y(:,:), h(:,:), lambda(:), s2(:) + double precision, allocatable :: y(:,:), h(:,:), h_p(:,:), lambda(:), s2(:) real, allocatable :: y_s(:,:) double precision, allocatable :: s_(:,:), s_tmp(:,:) double precision :: diag_h_mat_elem @@ -264,6 +281,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_ ! Small h(N_st_diag*itermax,N_st_diag*itermax), & + h_p(N_st_diag*itermax,N_st_diag*itermax), & y(N_st_diag*itermax,N_st_diag*itermax), & s_(N_st_diag*itermax,N_st_diag*itermax), & s_tmp(N_st_diag*itermax,N_st_diag*itermax), & @@ -408,27 +426,44 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_ ! Compute h_kl = = ! ------------------------------------------- + call dgemm('T','N', shift2, shift2, sze, & + 1.d0, U, size(U,1), W, size(W,1), & + 0.d0, h, size(h_p,1)) + ! Penalty method ! -------------- if (s2_eig) then - h = s_ + h_p = s_ do k=1,shift2 - h(k,k) = h(k,k) + S_z2_Sz - expected_s2 + h_p(k,k) = h_p(k,k) + S_z2_Sz - expected_s2 enddo alpha = 0.1d0 + h_p = h + alpha*h_p else + h_p = h alpha = 0.d0 endif - call dgemm('T','N', shift2, shift2, sze, & - 1.d0, U, size(U,1), W, size(W,1), & - alpha , h, size(h,1)) + ! Diagonalize h_p + ! --------------- - ! Diagonalize h - ! ------------- + call lapack_diag(lambda,y,h_p,size(h_p,1),shift2) - call lapack_diag(lambda,y,h,size(h,1),shift2) + ! Compute Energy for each eigenvector + ! ----------------------------------- + + call dgemm('N','N',shift2,shift2,shift2, & + 1.d0, h, size(h,1), y, size(y,1), & + 0.d0, s_tmp, size(s_tmp,1)) + + call dgemm('T','N',shift2,shift2,shift2, & + 1.d0, y, size(y,1), s_tmp, size(s_tmp,1), & + 0.d0, h, size(h,1)) + + do k=1,shift2 + lambda(k) = h(k,k) + enddo ! Compute S2 for each eigenvector ! ------------------------------- @@ -441,8 +476,6 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_ 1.d0, y, size(y,1), s_tmp, size(s_tmp,1), & 0.d0, s_, size(s_,1)) - - do k=1,shift2 s2(k) = s_(k,k) + S_z2_Sz enddo diff --git a/src/davidson/diagonalize_ci.irp.f b/src/davidson/diagonalize_ci.irp.f index d4b3b48c..8339406f 100644 --- a/src/davidson/diagonalize_ci.irp.f +++ b/src/davidson/diagonalize_ci.irp.f @@ -117,17 +117,24 @@ END_PROVIDER good_state_array = .False. call u_0_S2_u_0(s2_eigvalues,eigenvectors,N_det,psi_det,N_int,& N_det,size(eigenvectors,1)) - do j=1,N_det - ! Select at least n_states states with S^2 values closed to "expected_s2" - if(dabs(s2_eigvalues(j)-expected_s2).le.0.5d0)then - i_state +=1 - index_good_state_array(i_state) = j - good_state_array(j) = .True. - endif - if(i_state.eq.N_states) then - exit - endif - enddo + if (only_expected_s2) then + do j=1,N_det + ! Select at least n_states states with S^2 values closed to "expected_s2" + if(dabs(s2_eigvalues(j)-expected_s2).le.0.5d0)then + i_state +=1 + index_good_state_array(i_state) = j + good_state_array(j) = .True. + endif + if(i_state.eq.N_states) then + exit + endif + enddo + else + do j=1,N_det + index_good_state_array(j) = j + good_state_array(j) = .True. + enddo + endif if(i_state .ne.0)then ! Fill the first "i_state" states that have a correct S^2 value do j = 1, i_state @@ -185,6 +192,16 @@ END_PROVIDER CI_electronic_energy(j) = eigenvalues(j) enddo endif + do k=1,N_states_diag + CI_electronic_energy(k) = 0.d0 + do j=1,N_det + do i=1,N_det + CI_electronic_energy(k) += & + CI_eigenvectors(i,k) * CI_eigenvectors(j,k) * & + H_matrix_all_dets(i,j) + enddo + enddo + enddo deallocate(eigenvectors,eigenvalues) endif diff --git a/src/determinants/EZFIO.cfg b/src/determinants/EZFIO.cfg index 79cf29eb..12a2c2d1 100644 --- a/src/determinants/EZFIO.cfg +++ b/src/determinants/EZFIO.cfg @@ -10,12 +10,6 @@ doc: Maximum number of determinants to be printed with the program print_wf interface: ezfio,provider,ocaml default: 10000 -[n_det_max_full] -type: Det_number_max -doc: Maximum number of determinants where |H| is fully diagonalized -interface: ezfio,provider,ocaml -default: 1000 - [n_states] type: States_number doc: Number of states to consider diff --git a/src/determinants/determinants.irp.f b/src/determinants/determinants.irp.f index f29939ad..fbb3a813 100644 --- a/src/determinants/determinants.irp.f +++ b/src/determinants/determinants.irp.f @@ -1,22 +1,5 @@ use bitmasks -BEGIN_PROVIDER [ character*(64), diag_algorithm ] - implicit none - BEGIN_DOC - ! Diagonalization algorithm (Davidson or Lapack) - END_DOC - if (N_det > N_det_max_full) then - diag_algorithm = "Davidson" - else - diag_algorithm = "Lapack" - endif - - if (N_det < N_states) then - diag_algorithm = "Lapack" - endif -END_PROVIDER - - BEGIN_PROVIDER [ integer, N_det ] implicit none BEGIN_DOC diff --git a/src/fci/40.fci.bats b/src/fci/40.fci.bats index 3d6a4858..7083f9fe 100644 --- a/src/fci/40.fci.bats +++ b/src/fci/40.fci.bats @@ -59,7 +59,7 @@ function run_stoch() { @test "H2O2" { # 12.9214s qp set_file h2o2.ezfio qp set_mo_class --core="[1-2]" --act="[3-24]" --del="[25-38]" - run -151.004888189874 2.e-5 + run -151.004888189874 4.e-5 } @test "HBO" { # 13.3144s