From 18b038e06d41aebc78563313fbbf0b4dc95d9cac Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 1 Apr 2015 13:22:00 +0200 Subject: [PATCH 1/5] Corrected little bugs --- src/Dets/connected_to_ref.irp.f | 1 + src/Perturbation/perturbation_template.f | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/src/Dets/connected_to_ref.irp.f b/src/Dets/connected_to_ref.irp.f index fd87cc3c..1b9e0c35 100644 --- a/src/Dets/connected_to_ref.irp.f +++ b/src/Dets/connected_to_ref.irp.f @@ -60,6 +60,7 @@ integer function get_index_in_psi_det_sorted_bit(key,Nint) integer*8, external :: det_search_key logical :: is_in_wavefunction + is_in_wavefunction = .False. get_index_in_psi_det_sorted_bit = 0 ibegin = 1 iend = N_det+1 diff --git a/src/Perturbation/perturbation_template.f b/src/Perturbation/perturbation_template.f index 3627dd17..e289d409 100644 --- a/src/Perturbation/perturbation_template.f +++ b/src/Perturbation/perturbation_template.f @@ -25,7 +25,7 @@ subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,c ASSERT (N_st > 0) do i = 1,buffer_size - c_ref = connected_to_ref(buffer(1,1,i),psi_generators,Nint,i_generator,N_det) + c_ref = connected_to_ref(buffer(1,1,i),psi_generators,Nint,i_generator,N_det_generators) if (c_ref /= 0) then cycle From ad75f3fd682f8d192366a4246ee664f671c91330 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 1 Apr 2015 13:23:02 +0200 Subject: [PATCH 2/5] Working on MRCC --- src/MRCC/ASSUMPTIONS.rst | 0 src/MRCC/H_apply.irp.f | 10 +++++ src/MRCC/Makefile | 6 +++ src/MRCC/NEEDED_MODULES | 1 + src/MRCC/README.rst | 4 ++ src/MRCC/cas_sd.irp.f | 57 +++++++++++++++++++++++++++++ src/MRCC/mrcc.irp.f | 39 ++++++++++++++++++++ src/MRCC/mrcc_dress.irp.f | 54 +++++++++++++++++++++++++++ src/MRCC/mrcc_utils.irp.f | 77 +++++++++++++++++++++++++++++++++++++++ 9 files changed, 248 insertions(+) create mode 100644 src/MRCC/ASSUMPTIONS.rst create mode 100644 src/MRCC/H_apply.irp.f create mode 100644 src/MRCC/Makefile create mode 100644 src/MRCC/NEEDED_MODULES create mode 100644 src/MRCC/README.rst create mode 100644 src/MRCC/cas_sd.irp.f create mode 100644 src/MRCC/mrcc.irp.f create mode 100644 src/MRCC/mrcc_dress.irp.f create mode 100644 src/MRCC/mrcc_utils.irp.f diff --git a/src/MRCC/ASSUMPTIONS.rst b/src/MRCC/ASSUMPTIONS.rst new file mode 100644 index 00000000..e69de29b diff --git a/src/MRCC/H_apply.irp.f b/src/MRCC/H_apply.irp.f new file mode 100644 index 00000000..757627e1 --- /dev/null +++ b/src/MRCC/H_apply.irp.f @@ -0,0 +1,10 @@ +use bitmasks +BEGIN_SHELL [ /usr/bin/env python ] +from generate_h_apply import * + +s = H_apply("mrcc") +s.data["keys_work"] = "call mrcc_dress(i_generator,key_idx,keys_out,N_int,iproc)" +print s + +END_SHELL + diff --git a/src/MRCC/Makefile b/src/MRCC/Makefile new file mode 100644 index 00000000..06dc50ff --- /dev/null +++ b/src/MRCC/Makefile @@ -0,0 +1,6 @@ +# Define here all new external source files and objects.Don't forget to prefix the +# object files with IRPF90_temp/ +SRC= +OBJ= + +include $(QPACKAGE_ROOT)/src/Makefile.common diff --git a/src/MRCC/NEEDED_MODULES b/src/MRCC/NEEDED_MODULES new file mode 100644 index 00000000..31a27888 --- /dev/null +++ b/src/MRCC/NEEDED_MODULES @@ -0,0 +1 @@ +AOs BiInts Bitmask CAS_SD_selected Dets Electrons Ezfio_files Generators_CAS Hartree_Fock MOGuess MonoInts MOs Nuclei Output Perturbation Properties Selectors_full Utils diff --git a/src/MRCC/README.rst b/src/MRCC/README.rst new file mode 100644 index 00000000..ebb762c1 --- /dev/null +++ b/src/MRCC/README.rst @@ -0,0 +1,4 @@ +======= + Module +======= + diff --git a/src/MRCC/cas_sd.irp.f b/src/MRCC/cas_sd.irp.f new file mode 100644 index 00000000..e3a8f6ec --- /dev/null +++ b/src/MRCC/cas_sd.irp.f @@ -0,0 +1,57 @@ +program full_ci + implicit none + integer :: i,k + + + double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:) + integer :: N_st, degree + N_st = N_states + allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st)) + character*(64) :: perturbation + + pt2 = 1.d0 + diag_algorithm = "Lapack" + if (N_det > n_det_max_fci) then + call diagonalize_CI + call save_wavefunction + psi_det = psi_det_sorted + psi_coef = psi_coef_sorted + N_det = n_det_max_fci + soft_touch N_det psi_det psi_coef + call diagonalize_CI + call save_wavefunction + print *, 'N_det = ', N_det + print *, 'N_states = ', N_states + print *, 'PT2 = ', pt2 + print *, 'E = ', CI_energy + print *, 'E+PT2 = ', CI_energy+pt2 + print *, '-----' + endif + + do while (N_det < n_det_max_fci.and.maxval(abs(pt2(1:N_st))) > pt2_max) + call H_apply_FCI(pt2, norm_pert, H_pert_diag, N_st) + + PROVIDE psi_coef + PROVIDE psi_det + PROVIDE psi_det_sorted + + if (N_det > n_det_max_fci) then + psi_det = psi_det_sorted + psi_coef = psi_coef_sorted + N_det = n_det_max_fci + soft_touch N_det psi_det psi_coef + endif + call diagonalize_CI + call save_wavefunction + print *, 'N_det = ', N_det + print *, 'N_states = ', N_states + print *, 'PT2 = ', pt2 + print *, 'E = ', CI_energy + print *, 'E+PT2 = ', CI_energy+pt2 + print *, '-----' + call ezfio_set_full_ci_energy(CI_energy) + if (abort_all) then + exit + endif + enddo +end diff --git a/src/MRCC/mrcc.irp.f b/src/MRCC/mrcc.irp.f new file mode 100644 index 00000000..4f2a6858 --- /dev/null +++ b/src/MRCC/mrcc.irp.f @@ -0,0 +1,39 @@ +program mrcc + implicit none + read_wf = .True. + TOUCH read_wf + + + print *, N_det + print *, N_det_cas + print *, N_det_sd +! psi_cas, (N_int,2,N_det_generators) ] +!psi_cas_coefs, (N_det_generators,n_states) ] +!psi_sd, (N_int,2,psi_det_size) ] +!psi_sd_coefs, (psi_det_size,n_states) ] + + call update_generators + integer :: i + print *, 'CAS' + print *, '===' + do i=1,N_det_cas + print *, psi_cas_coefs(i,:) + call debug_det(psi_cas(1,1,i),N_int) + enddo + + print *, 'SD' + print *, '==' + do i=1,N_det_sd + print *, psi_sd_coefs(i,:) + call debug_det(psi_sd(1,1,i),N_int) + enddo + + double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:) + integer :: N_st, degree + N_st = N_states + allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st)) + + print *, 'MRCC' + print *, '====' + call H_apply_mrcc(pt2, norm_pert, H_pert_diag, N_st) +end diff --git a/src/MRCC/mrcc_dress.irp.f b/src/MRCC/mrcc_dress.irp.f new file mode 100644 index 00000000..86ffc9f5 --- /dev/null +++ b/src/MRCC/mrcc_dress.irp.f @@ -0,0 +1,54 @@ +subroutine mrcc_dress(i_generator,n_selected,det_buffer,Nint,iproc) + use bitmasks + implicit none + + integer, intent(in) :: i_generator,n_selected, Nint, iproc + integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected) + integer :: i,j,k + integer :: new_size + logical :: is_in_wavefunction + double precision :: degree(N_det_cas) + integer :: idx(0:N_det_cas) + logical :: good + + integer(bit_kind) :: tq(Nint,2,n_selected) + integer :: N_tq, c_ref + integer :: connected_to_ref + + N_tq = 0 + do i=1,N_selected + + c_ref = connected_to_ref(det_buffer(1,1,i),psi_generators,Nint, & + i_generator,N_det_generators) + + if (c_ref /= 0) then + cycle + endif + + ! Select determinants that are triple or quadruple excitations + ! from the CAS + good = .True. + call get_excitation_degree_vector(psi_cas,det_buffer(1,1,i),degree,Nint,N_det_cas,idx) + do k=1,idx(0) + if (degree(k) < 3) then + good = .False. + exit + endif + enddo + if (good) then + if (.not. is_in_wavefunction(det_buffer(1,1,i),Nint,N_det)) then + N_tq += 1 + do k=1,N_int + tq(k,1,N_tq) = det_buffer(k,1,i) + tq(k,2,N_tq) = det_buffer(k,2,i) + enddo + endif + endif + enddo + + print *, N_tq + do i=1,N_tq + call debug_det(det_buffer(1,1,i),Nint) + enddo +end + diff --git a/src/MRCC/mrcc_utils.irp.f b/src/MRCC/mrcc_utils.irp.f new file mode 100644 index 00000000..327c220e --- /dev/null +++ b/src/MRCC/mrcc_utils.irp.f @@ -0,0 +1,77 @@ + use bitmasks + BEGIN_PROVIDER [ integer(bit_kind), psi_cas, (N_int,2,N_det_generators) ] +&BEGIN_PROVIDER [ double precision, psi_cas_coefs, (N_det_generators,n_states) ] +&BEGIN_PROVIDER [ integer(bit_kind), psi_sd, (N_int,2,N_det) ] +&BEGIN_PROVIDER [ double precision, psi_sd_coefs, (N_det,n_states) ] +&BEGIN_PROVIDER [ integer, idx_cas, (N_det_generators) ] +&BEGIN_PROVIDER [ integer, idx_sd, (N_det) ] +&BEGIN_PROVIDER [ integer, N_det_sd] +&BEGIN_PROVIDER [ integer, N_det_cas] + implicit none + BEGIN_DOC + ! SD + END_DOC + integer :: i_cas,i_sd,j,k + integer :: degree + logical :: in_cas + i_cas=0 + i_sd =0 + do k=1,N_det + in_cas = .False. + do j=1,n_det_generators + call get_excitation_degree(psi_generators(1,1,j), psi_det(1,1,k), degree, N_int) + if (degree == 0) then + i_cas += 1 + psi_cas(1:N_int,1:2,i_cas) = psi_det(1:N_int,1:2,k) + psi_cas_coefs(i_cas,1:N_states) = psi_coef(k,1:N_states) + in_cas = .True. + idx_cas(i_cas) = k + exit + endif + enddo + if (.not.in_cas) then + double precision :: hij + i_sd += 1 + psi_sd(1:N_int,1:2,i_sd) = psi_det(1:N_int,1:2,k) + psi_sd_coefs(i_sd,1:N_states) = psi_coef(k,1:N_states) + idx_sd(i_sd) = k + endif + enddo + N_det_sd = i_sd + N_det_cas = i_cas +END_PROVIDER + +BEGIN_PROVIDER [ double precision, lambda_mrcc, (psi_det_size,n_states) ] + implicit none + BEGIN_DOC + ! cm/ + END_DOC + integer :: i,k + double precision :: ihpsi(N_states) + do i=1,N_det_sd + call i_h_psi(psi_sd(1,1,i), psi_cas, psi_cas_coefs, N_int, N_det_cas, & + size(psi_cas_coefs,1), n_states, ihpsi) + double precision :: hij + do k=1,N_states + if (dabs(ihpsi(k)) < 1.d-6) then + lambda_mrcc(i,k) = 0.d0 + else + lambda_mrcc(i,k) = psi_sd_coefs(i,k)/ihpsi(k) + lambda_mrcc(i,k) = min( lambda_mrcc (i,k),0.d0 ) + endif + enddo + enddo +END_PROVIDER + +subroutine update_generators + implicit none + integer :: i,j,k + n_det_generators = N_det_sd + do k=1,N_det_sd + do j=1,2 + do i=1,N_int + psi_generators(i,j,k) = psi_sd(i,j,k) + enddo + enddo + enddo +end From a95adca2c702168f79873eeed3bebcb72a752e93 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 2 Apr 2015 10:13:33 +0200 Subject: [PATCH 3/5] Working on MRCC --- src/CAS_SD_selected/cas_sd.irp.f | 38 +++++++ src/Dets/connected_to_ref.irp.f | 60 +++++----- src/MRCC/H_apply.irp.f | 18 ++- src/MRCC/NEEDED_MODULES | 3 +- src/MRCC/README.rst | 16 ++- src/MRCC/cas_sd.irp.f | 57 ---------- src/MRCC/mrcc.irp.f | 40 ++++--- src/MRCC/mrcc_dress.irp.f | 188 +++++++++++++++++++++++++++++-- src/MRCC/mrcc_utils.irp.f | 125 +++++++++++++++----- src/NEEDED_MODULES | 2 +- 10 files changed, 403 insertions(+), 144 deletions(-) delete mode 100644 src/MRCC/cas_sd.irp.f diff --git a/src/CAS_SD_selected/cas_sd.irp.f b/src/CAS_SD_selected/cas_sd.irp.f index e3a8f6ec..2d115117 100644 --- a/src/CAS_SD_selected/cas_sd.irp.f +++ b/src/CAS_SD_selected/cas_sd.irp.f @@ -54,4 +54,42 @@ program full_ci exit endif enddo + + ! Check that it is a CAS-SD + logical :: in_cas + integer :: exc_max + exc_max = 0 + print *, 'CAS determinants' + do i=1,N_det_generators + do k=i,N_det_generators + call get_excitation_degree(psi_generators(1,1,k),psi_det(1,1,i),degree,N_int) + exc_max = max(exc_max,degree) + enddo + call debug_det(psi_generators(1,1,i),N_int) + print *, '' + enddo + print *, 'Max excitation degree in the CAS :', exc_max + do i=1,N_det + in_cas = .False. + do k=1,N_det_generators + call get_excitation_degree(psi_generators(1,1,k),psi_det(1,1,i),degree,N_int) + if (degree == 0) then + in_cas = .True. + exit + endif + enddo + if (.not.in_cas) then + do k=i,N_det + call get_excitation_degree(psi_det(1,1,k),psi_det(1,1,i),degree,N_int) + if (degree > exc_max+2) then + print *, 'Error : This is not a CAS-SD : ' + print *, 'CAS determinant:' + call debug_det(psi_det(1,1,i),N_int) + print *, 'Excited determinant:', degree + call debug_det(psi_det(1,1,k),N_int) + stop + endif + enddo + endif + enddo end diff --git a/src/Dets/connected_to_ref.irp.f b/src/Dets/connected_to_ref.irp.f index 1b9e0c35..07331782 100644 --- a/src/Dets/connected_to_ref.irp.f +++ b/src/Dets/connected_to_ref.irp.f @@ -94,9 +94,10 @@ integer function get_index_in_psi_det_sorted_bit(key,Nint) endif enddo i += 1 - if (i > N_det) then - return - endif + +! if (i > N_det) then +! return +! endif !DIR$ FORCEINLINE do while (det_search_key(psi_det_sorted_bit(1,1,i),Nint) == det_ref) @@ -105,50 +106,49 @@ integer function get_index_in_psi_det_sorted_bit(key,Nint) continue else is_in_wavefunction = .True. - !DEC$ LOOP COUNT MIN(3) + !DIR$ IVDEP + !DIR$ LOOP COUNT MIN(3) do l=2,Nint if ( (key(l,1) /= psi_det_sorted_bit(l,1,i)).or. & (key(l,2) /= psi_det_sorted_bit(l,2,i)) ) then is_in_wavefunction = .False. - exit endif enddo if (is_in_wavefunction) then + get_index_in_psi_det_sorted_bit = i exit +! return endif endif i += 1 if (i > N_det) then - return - ! exit + exit +! return endif enddo - if (is_in_wavefunction) then - get_index_in_psi_det_sorted_bit = i - endif ! DEBUG is_in_wf -! if (is_in_wavefunction) then -! degree = 1 -! do i=1,N_det -! integer :: degree -! call get_excitation_degree(key,psi_det(1,1,i),degree,N_int) -! if (degree == 0) then -! exit -! endif -! enddo -! if (degree /=0) then -! stop 'pouet 1' -! endif -! else -! do i=1,N_det -! call get_excitation_degree(key,psi_det(1,1,i),degree,N_int) -! if (degree == 0) then -! stop 'pouet 2' -! endif -! enddo -! endif + if (is_in_wavefunction) then + degree = 1 + do i=1,N_det + integer :: degree + call get_excitation_degree(key,psi_det(1,1,i),degree,N_int) + if (degree == 0) then + exit + endif + enddo + if (degree /=0) then + stop 'pouet 1' + endif + else + do i=1,N_det + call get_excitation_degree(key,psi_det(1,1,i),degree,N_int) + if (degree == 0) then + stop 'pouet 2' + endif + enddo + endif ! END DEBUG is_in_wf end diff --git a/src/MRCC/H_apply.irp.f b/src/MRCC/H_apply.irp.f index 757627e1..4f284112 100644 --- a/src/MRCC/H_apply.irp.f +++ b/src/MRCC/H_apply.irp.f @@ -3,7 +3,23 @@ BEGIN_SHELL [ /usr/bin/env python ] from generate_h_apply import * s = H_apply("mrcc") -s.data["keys_work"] = "call mrcc_dress(i_generator,key_idx,keys_out,N_int,iproc)" +s.data["parameters"] = ", delta_ij_sd_, Ndet_sd" +s.data["declarations"] += """ + integer, intent(in) :: Ndet_sd + double precision, intent(in) :: delta_ij_sd_(Ndet_sd,Ndet_sd,*) +""" +s.data["keys_work"] = "call mrcc_dress(delta_ij_sd_,Ndet_sd,i_generator,key_idx,keys_out,N_int,iproc)" +s.data["params_post"] += ", delta_ij_sd_, Ndet_sd" +s.data["params_main"] += "delta_ij_sd_, Ndet_sd" +s.data["decls_main"] += """ + integer, intent(in) :: Ndet_sd + double precision, intent(in) :: delta_ij_sd_(Ndet_sd,Ndet_sd,*) +""" +s.data["finalization"] = "" +s.data["copy_buffer"] = "" +s.data["generate_psi_guess"] = "" +s.data["size_max"] = "256" + print s END_SHELL diff --git a/src/MRCC/NEEDED_MODULES b/src/MRCC/NEEDED_MODULES index 31a27888..3711786f 100644 --- a/src/MRCC/NEEDED_MODULES +++ b/src/MRCC/NEEDED_MODULES @@ -1 +1,2 @@ -AOs BiInts Bitmask CAS_SD_selected Dets Electrons Ezfio_files Generators_CAS Hartree_Fock MOGuess MonoInts MOs Nuclei Output Perturbation Properties Selectors_full Utils +AOs BiInts Bitmask CAS_SD_selected Dets Electrons Ezfio_files Generators_CAS Generators_full Hartree_Fock MOGuess MonoInts MOs Nuclei Output Perturbation Properties Selectors_full Utils + diff --git a/src/MRCC/README.rst b/src/MRCC/README.rst index ebb762c1..1fbcccdb 100644 --- a/src/MRCC/README.rst +++ b/src/MRCC/README.rst @@ -1,4 +1,14 @@ -======= - Module -======= +-4.142795384334731 + +4.695183071437694E-002 + Determinant 64 + --------------------------------------------- + 000000000000002E|000000000000002E + |-+++-+----------------------------------------------------------| + |-+++-+----------------------------------------------------------| + +CAS-SD: -4.14214374069306 + : -4.14230904320551 + +E0 = -11.5634986758976 diff --git a/src/MRCC/cas_sd.irp.f b/src/MRCC/cas_sd.irp.f deleted file mode 100644 index e3a8f6ec..00000000 --- a/src/MRCC/cas_sd.irp.f +++ /dev/null @@ -1,57 +0,0 @@ -program full_ci - implicit none - integer :: i,k - - - double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:) - integer :: N_st, degree - N_st = N_states - allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st)) - character*(64) :: perturbation - - pt2 = 1.d0 - diag_algorithm = "Lapack" - if (N_det > n_det_max_fci) then - call diagonalize_CI - call save_wavefunction - psi_det = psi_det_sorted - psi_coef = psi_coef_sorted - N_det = n_det_max_fci - soft_touch N_det psi_det psi_coef - call diagonalize_CI - call save_wavefunction - print *, 'N_det = ', N_det - print *, 'N_states = ', N_states - print *, 'PT2 = ', pt2 - print *, 'E = ', CI_energy - print *, 'E+PT2 = ', CI_energy+pt2 - print *, '-----' - endif - - do while (N_det < n_det_max_fci.and.maxval(abs(pt2(1:N_st))) > pt2_max) - call H_apply_FCI(pt2, norm_pert, H_pert_diag, N_st) - - PROVIDE psi_coef - PROVIDE psi_det - PROVIDE psi_det_sorted - - if (N_det > n_det_max_fci) then - psi_det = psi_det_sorted - psi_coef = psi_coef_sorted - N_det = n_det_max_fci - soft_touch N_det psi_det psi_coef - endif - call diagonalize_CI - call save_wavefunction - print *, 'N_det = ', N_det - print *, 'N_states = ', N_states - print *, 'PT2 = ', pt2 - print *, 'E = ', CI_energy - print *, 'E+PT2 = ', CI_energy+pt2 - print *, '-----' - call ezfio_set_full_ci_energy(CI_energy) - if (abort_all) then - exit - endif - enddo -end diff --git a/src/MRCC/mrcc.irp.f b/src/MRCC/mrcc.irp.f index 4f2a6858..663b9f5a 100644 --- a/src/MRCC/mrcc.irp.f +++ b/src/MRCC/mrcc.irp.f @@ -2,17 +2,17 @@ program mrcc implicit none read_wf = .True. TOUCH read_wf + call run +end +subroutine run + implicit none print *, N_det print *, N_det_cas print *, N_det_sd -! psi_cas, (N_int,2,N_det_generators) ] -!psi_cas_coefs, (N_det_generators,n_states) ] -!psi_sd, (N_int,2,psi_det_size) ] -!psi_sd_coefs, (psi_det_size,n_states) ] - call update_generators +! call update_generators integer :: i print *, 'CAS' print *, '===' @@ -21,19 +21,23 @@ program mrcc call debug_det(psi_cas(1,1,i),N_int) enddo - print *, 'SD' - print *, '==' - do i=1,N_det_sd - print *, psi_sd_coefs(i,:) - call debug_det(psi_sd(1,1,i),N_int) - enddo - - double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:) - integer :: N_st, degree - N_st = N_states - allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st)) - +! print *, 'SD' +! print *, '==' +! do i=1,N_det_sd +! print *, psi_sd_coefs(i,:) +! call debug_det(psi_sd(1,1,i),N_int) +! enddo +! print *, 'MRCC' print *, '====' - call H_apply_mrcc(pt2, norm_pert, H_pert_diag, N_st) + print *, ci_energy(:) + print *, h_matrix_all_dets(3,3), delta_ij(3,3,1) + print *, h_matrix_all_dets(3,3), delta_ij(3,3,1) + print *, ci_energy_dressed(:) +! print *, 'max', maxval(delta_ij_sd) +! print *, 'min', minval(delta_ij_sd) +! +! do i=1,N_det +! print '(10(F10.6,X))', delta_ij(i,1:N_det,1) +! enddo end diff --git a/src/MRCC/mrcc_dress.irp.f b/src/MRCC/mrcc_dress.irp.f index 86ffc9f5..56257556 100644 --- a/src/MRCC/mrcc_dress.irp.f +++ b/src/MRCC/mrcc_dress.irp.f @@ -1,14 +1,17 @@ -subroutine mrcc_dress(i_generator,n_selected,det_buffer,Nint,iproc) +subroutine mrcc_dress(delta_ij_sd_,Ndet_sd,i_generator,n_selected,det_buffer,Nint,iproc) use bitmasks implicit none integer, intent(in) :: i_generator,n_selected, Nint, iproc + integer, intent(in) :: Ndet_sd + double precision, intent(inout) :: delta_ij_sd_(Ndet_sd,Ndet_sd,*) + integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected) - integer :: i,j,k + integer :: i,j,k,m integer :: new_size logical :: is_in_wavefunction - double precision :: degree(N_det_cas) - integer :: idx(0:N_det_cas) + integer :: degree(psi_det_size) + integer :: idx(0:psi_det_size) logical :: good integer(bit_kind) :: tq(Nint,2,n_selected) @@ -17,7 +20,6 @@ subroutine mrcc_dress(i_generator,n_selected,det_buffer,Nint,iproc) N_tq = 0 do i=1,N_selected - c_ref = connected_to_ref(det_buffer(1,1,i),psi_generators,Nint, & i_generator,N_det_generators) @@ -46,9 +48,181 @@ subroutine mrcc_dress(i_generator,n_selected,det_buffer,Nint,iproc) endif enddo - print *, N_tq + ! Compute / (E0 - Haa) + double precision :: hka, haa + double precision :: haj + double precision :: f(N_states) + + +! call i_h_j(psi_det(1,1,1), psi_det(1,1,64),Nint,hka) +! call debug_det(psi_det(1,1,1), N_int) +! call debug_det(psi_det(1,1,64), N_int) + double precision :: phase + integer :: exc(0:2,2,2) +! call get_excitation(psi_det(1,1,1),psi_det(1,1,64),exc,degree(1),phase,Nint) + integer :: h1, p1, h2, p2, s1, s2 +! call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) +! print *, hka +! print *, h1, p1, h2, p2 +! print *, s1, s2 +! pause + + + do i=1,N_tq - call debug_det(det_buffer(1,1,i),Nint) + call get_excitation_degree_vector(psi_sd,tq(1,1,i),degree,Nint,Ndet_sd,idx) + call i_h_j(tq(1,1,i),tq(1,1,i),Nint,haa) + do m=1,N_states + f(m) = 1.d0/(ci_electronic_energy(m)-haa) + enddo + do k=1,idx(0) + call i_h_j(tq(1,1,i),psi_sd(1,1,idx(k)),Nint,hka) + do j=k,idx(0) + call i_h_j(tq(1,1,i),psi_sd(1,1,idx(j)),Nint,haj) + do m=1,N_states + delta_ij_sd_(idx(k), idx(j),m) += haj*hka* f(m) + delta_ij_sd_(idx(j), idx(k),m) += haj*hka* f(m) + enddo + call get_excitation(tq(1,1,i),psi_sd(1,1,idx(j)),exc,degree(1),phase,Nint) + call decode_exc(exc,degree(1),h1,p1,h2,p2,s1,s2) + if ( (h1 == 1).and. & + (p1 == 6).and. & + (h2 == 1).and. & + (p2 == 6).and. & + (s1 == 1).and. & + (s2 == 2) ) then + call debug_det(tq(1,1,i), N_int) + call debug_det(psi_sd(1,1,idx(j)), N_int) + print *, haj + pause + endif + + + enddo + enddo enddo end +BEGIN_PROVIDER [ double precision, delta_ij_sd, (N_det_sd, N_det_sd,N_states) ] + implicit none + BEGIN_DOC + ! Dressing matrix in SD basis + END_DOC + delta_ij_sd = 0.d0 + call H_apply_mrcc(delta_ij_sd,N_det_sd) + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, delta_ij, (N_det,N_det,N_states) ] + implicit none + BEGIN_DOC + ! Dressing matrix in N_det basis + END_DOC + integer :: i,j,m + delta_ij = 0.d0 + do m=1,N_states + do j=1,N_det_sd + do i=1,N_det_sd + delta_ij(idx_sd(i),idx_sd(j),m) = delta_ij_sd(i,j,m) + enddo + enddo + enddo +END_PROVIDER + +BEGIN_PROVIDER [ double precision, h_matrix_dressed, (N_det,N_det) ] + implicit none + BEGIN_DOC + ! Dressed H with Delta_ij + END_DOC + integer :: i, j + do j=1,N_det + do i=1,N_det + h_matrix_dressed(i,j) = h_matrix_all_dets(i,j) + delta_ij(i,j,1) + if (i==j) then + print *, i, delta_ij(i,j,1), h_matrix_all_dets(i,j) + endif + enddo + enddo + +END_PROVIDER + + + BEGIN_PROVIDER [ double precision, CI_electronic_energy_dressed, (N_states_diag) ] +&BEGIN_PROVIDER [ double precision, CI_eigenvectors_dressed, (N_det,N_states_diag) ] +&BEGIN_PROVIDER [ double precision, CI_eigenvectors_s2_dressed, (N_states_diag) ] + implicit none + BEGIN_DOC + ! Eigenvectors/values of the CI matrix + END_DOC + integer :: i,j + + do j=1,N_states_diag + do i=1,N_det + CI_eigenvectors_dressed(i,j) = psi_coef(i,j) + enddo + enddo + + if (diag_algorithm == "Davidson") then + + stop 'use Lapack' +! call davidson_diag(psi_det,CI_eigenvectors_dressed,CI_electronic_energy_dressed, & +! size(CI_eigenvectors_dressed,1),N_det,N_states_diag,N_int,output_Dets) + + else if (diag_algorithm == "Lapack") then + + double precision, allocatable :: eigenvectors(:,:), eigenvalues(:) + allocate (eigenvectors(size(H_matrix_dressed,1),N_det)) + allocate (eigenvalues(N_det)) + call lapack_diag(eigenvalues,eigenvectors, & + H_matrix_dressed,size(H_matrix_dressed,1),N_det) + CI_electronic_energy_dressed(:) = 0.d0 + do i=1,N_det + CI_eigenvectors_dressed(i,1) = eigenvectors(i,1) + enddo + integer :: i_state + double precision :: s2 + i_state = 0 + do j=1,N_det + call get_s2_u0(psi_det,eigenvectors(1,j),N_det,N_det,s2) + if(dabs(s2-expected_s2).le.0.3d0)then + i_state += 1 + do i=1,N_det + CI_eigenvectors_dressed(i,i_state) = eigenvectors(i,j) + enddo + CI_electronic_energy_dressed(i_state) = eigenvalues(j) + CI_eigenvectors_s2_dressed(i_state) = s2 + endif + if (i_state.ge.N_states_diag) then + exit + endif + enddo +! if(i_state < min(N_states_diag,N_det))then +! print *, 'pb with the number of states' +! print *, 'i_state = ',i_state +! print *, 'N_states_diag ',N_states_diag +! print *,'stopping ...' +! stop +! endif + deallocate(eigenvectors,eigenvalues) + endif + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, CI_energy_dressed, (N_states_diag) ] + implicit none + BEGIN_DOC + ! N_states lowest eigenvalues of the dressed CI matrix + END_DOC + + integer :: j + character*(8) :: st + call write_time(output_Dets) + do j=1,N_states_diag + CI_energy_dressed(j) = CI_electronic_energy_dressed(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 + diff --git a/src/MRCC/mrcc_utils.irp.f b/src/MRCC/mrcc_utils.irp.f index 327c220e..2c5dc811 100644 --- a/src/MRCC/mrcc_utils.irp.f +++ b/src/MRCC/mrcc_utils.irp.f @@ -1,31 +1,117 @@ use bitmasks - BEGIN_PROVIDER [ integer(bit_kind), psi_cas, (N_int,2,N_det_generators) ] -&BEGIN_PROVIDER [ double precision, psi_cas_coefs, (N_det_generators,n_states) ] -&BEGIN_PROVIDER [ integer(bit_kind), psi_sd, (N_int,2,N_det) ] + +BEGIN_PROVIDER [ integer(bit_kind), cas_bitmask, (N_int,2,N_cas_bitmask) ] + implicit none + BEGIN_DOC + ! Bitmasks for CAS reference determinants. (N_int, alpha/beta, CAS reference) + END_DOC + logical :: exists + integer :: i + PROVIDE ezfio_filename + + call ezfio_has_bitmasks_cas(exists) + if (exists) then + call ezfio_get_bitmasks_cas(cas_bitmask) + else + do i=1,N_cas_bitmask + cas_bitmask(:,:,i) = iand(not(HF_bitmask(:,:)),full_ijkl_bitmask(:,:)) + enddo + endif + +END_PROVIDER + +BEGIN_PROVIDER [ integer, N_det_cas ] + implicit none + BEGIN_DOC + ! Number of generator detetrminants + END_DOC + integer :: i,k,l + logical :: good + call write_time(output_dets) + N_det_cas = 0 + do i=1,N_det + do l=1,n_cas_bitmask + good = .True. + do k=1,N_int + good = good .and. ( & + iand(not(cas_bitmask(k,1,l)), psi_det(k,1,i)) == & + iand(not(cas_bitmask(k,1,l)), psi_det(k,1,1)) ) .and. ( & + iand(not(cas_bitmask(k,2,l)), psi_det(k,2,i)) == & + iand(not(cas_bitmask(k,2,l)), psi_det(k,2,1)) ) + enddo + if (good) then + exit + endif + enddo + if (good) then + N_det_cas += 1 + endif + enddo + N_det_cas = max(N_det_cas, 1) + call write_int(output_dets,N_det_cas, 'Number of determinants in the CAS') +END_PROVIDER + + BEGIN_PROVIDER [ integer(bit_kind), psi_cas, (N_int,2,N_det_cas) ] +&BEGIN_PROVIDER [ double precision, psi_cas_coefs, (N_det_cas,n_states) ] +&BEGIN_PROVIDER [ integer, idx_cas, (N_det_cas) ] + implicit none + BEGIN_DOC + ! For Single reference wave functions, the generator is the + ! Hartree-Fock determinant + END_DOC + integer :: i, k, l, m + logical :: good + m=0 + do i=1,N_det + do l=1,n_cas_bitmask + good = .True. + do k=1,N_int + good = good .and. ( & + iand(not(cas_bitmask(k,1,l)), psi_det(k,1,i)) == & + iand(not(cas_bitmask(k,1,l)), psi_det(k,1,1)) ) .and. ( & + iand(not(cas_bitmask(k,2,l)), psi_det(k,2,i)) == & + iand(not(cas_bitmask(k,2,l)), psi_det(k,2,1)) ) + enddo + if (good) then + exit + endif + enddo + if (good) then + m = m+1 + do k=1,N_int + psi_cas(k,1,m) = psi_det(k,1,i) + psi_cas(k,2,m) = psi_det(k,2,i) + enddo + idx_cas(m) = i + do k=1,N_states + psi_cas_coefs(m,k) = psi_coef(i,k) + enddo + endif + enddo + +END_PROVIDER + + + + + BEGIN_PROVIDER [ integer(bit_kind), psi_sd, (N_int,2,N_det) ] &BEGIN_PROVIDER [ double precision, psi_sd_coefs, (N_det,n_states) ] -&BEGIN_PROVIDER [ integer, idx_cas, (N_det_generators) ] &BEGIN_PROVIDER [ integer, idx_sd, (N_det) ] &BEGIN_PROVIDER [ integer, N_det_sd] -&BEGIN_PROVIDER [ integer, N_det_cas] implicit none BEGIN_DOC ! SD END_DOC - integer :: i_cas,i_sd,j,k + integer :: i_sd,j,k integer :: degree logical :: in_cas - i_cas=0 i_sd =0 do k=1,N_det in_cas = .False. - do j=1,n_det_generators - call get_excitation_degree(psi_generators(1,1,j), psi_det(1,1,k), degree, N_int) + do j=1,N_det_cas + call get_excitation_degree(psi_cas(1,1,j), psi_det(1,1,k), degree, N_int) if (degree == 0) then - i_cas += 1 - psi_cas(1:N_int,1:2,i_cas) = psi_det(1:N_int,1:2,k) - psi_cas_coefs(i_cas,1:N_states) = psi_coef(k,1:N_states) in_cas = .True. - idx_cas(i_cas) = k exit endif enddo @@ -38,7 +124,6 @@ endif enddo N_det_sd = i_sd - N_det_cas = i_cas END_PROVIDER BEGIN_PROVIDER [ double precision, lambda_mrcc, (psi_det_size,n_states) ] @@ -63,15 +148,3 @@ BEGIN_PROVIDER [ double precision, lambda_mrcc, (psi_det_size,n_states) ] enddo END_PROVIDER -subroutine update_generators - implicit none - integer :: i,j,k - n_det_generators = N_det_sd - do k=1,N_det_sd - do j=1,2 - do i=1,N_int - psi_generators(i,j,k) = psi_sd(i,j,k) - enddo - enddo - enddo -end diff --git a/src/NEEDED_MODULES b/src/NEEDED_MODULES index 3df2ab88..7c84ef42 100644 --- a/src/NEEDED_MODULES +++ b/src/NEEDED_MODULES @@ -1 +1 @@ -AOs BiInts Bitmask CID CID_SC2_selected CID_selected CIS CISD CISD_selected CISD_SC2_selected Dets Electrons Ezfio_files Full_CI Generators_full Hartree_Fock MOGuess MonoInts MOs MP2 Nuclei Output Selectors_full Utils Molden FCIdump Generators_CAS CAS_SD_selected DDCI_selected +AOs BiInts Bitmask CID CID_SC2_selected CID_selected CIS CISD CISD_selected CISD_SC2_selected Dets Electrons Ezfio_files Full_CI Generators_full Hartree_Fock MOGuess MonoInts MOs MP2 Nuclei Output Selectors_full Utils Molden FCIdump Generators_CAS CAS_SD_selected DDCI_selected MRCC From c8c64a76ffbf411ac115b50b664b1bfcc26c0fc1 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 2 Apr 2015 10:26:45 +0200 Subject: [PATCH 4/5] Renamed psi_generators in psi_det_generators --- src/CAS_SD_selected/cas_sd.irp.f | 6 ++--- src/Dets/H_apply_template.f | 32 ++++++++++++------------ src/Generators_CAS/generators.irp.f | 8 +++--- src/Generators_full/generators.irp.f | 8 +++--- src/Generators_restart/generators.irp.f | 10 ++++---- src/MRCC/mrcc_dress.irp.f | 2 +- src/Perturbation/perturbation_template.f | 4 +-- src/SingleRefMethod/generators.irp.f | 17 ++++++------- 8 files changed, 42 insertions(+), 45 deletions(-) diff --git a/src/CAS_SD_selected/cas_sd.irp.f b/src/CAS_SD_selected/cas_sd.irp.f index 2d115117..e3721ca7 100644 --- a/src/CAS_SD_selected/cas_sd.irp.f +++ b/src/CAS_SD_selected/cas_sd.irp.f @@ -62,17 +62,17 @@ program full_ci print *, 'CAS determinants' do i=1,N_det_generators do k=i,N_det_generators - call get_excitation_degree(psi_generators(1,1,k),psi_det(1,1,i),degree,N_int) + call get_excitation_degree(psi_det_generators(1,1,k),psi_det_generators(1,1,i),degree,N_int) exc_max = max(exc_max,degree) enddo - call debug_det(psi_generators(1,1,i),N_int) + call debug_det(psi_det_generators(1,1,i),N_int) print *, '' enddo print *, 'Max excitation degree in the CAS :', exc_max do i=1,N_det in_cas = .False. do k=1,N_det_generators - call get_excitation_degree(psi_generators(1,1,k),psi_det(1,1,i),degree,N_int) + call get_excitation_degree(psi_det_generators(1,1,k),psi_det(1,1,i),degree,N_int) if (degree == 0) then in_cas = .True. exit diff --git a/src/Dets/H_apply_template.f b/src/Dets/H_apply_template.f index f4486c36..e44562e7 100644 --- a/src/Dets/H_apply_template.f +++ b/src/Dets/H_apply_template.f @@ -421,32 +421,32 @@ subroutine $subroutine($params_main) do k=1,N_int mask(k,ispin,s_hole) = & iand(generators_bitmask(k,ispin,s_hole,i_bitmask_gen), & - psi_generators(k,ispin,i_generator) ) + psi_det_generators(k,ispin,i_generator) ) mask(k,ispin,s_part) = & iand(generators_bitmask(k,ispin,s_part,i_bitmask_gen), & - not(psi_generators(k,ispin,i_generator)) ) + not(psi_det_generators(k,ispin,i_generator)) ) mask(k,ispin,d_hole1) = & iand(generators_bitmask(k,ispin,d_hole1,i_bitmask_gen), & - psi_generators(k,ispin,i_generator) ) + psi_det_generators(k,ispin,i_generator) ) mask(k,ispin,d_part1) = & iand(generators_bitmask(k,ispin,d_part1,i_bitmask_gen), & - not(psi_generators(k,ispin,i_generator)) ) + not(psi_det_generators(k,ispin,i_generator)) ) mask(k,ispin,d_hole2) = & iand(generators_bitmask(k,ispin,d_hole2,i_bitmask_gen), & - psi_generators(k,ispin,i_generator) ) + psi_det_generators(k,ispin,i_generator) ) mask(k,ispin,d_part2) = & iand(generators_bitmask(k,ispin,d_part2,i_bitmask_gen), & - not(psi_generators(k,ispin,i_generator)) ) + not(psi_det_generators(k,ispin,i_generator)) ) enddo enddo if($do_double_excitations)then - call $subroutine_diexc(psi_generators(1,1,i_generator), & + call $subroutine_diexc(psi_det_generators(1,1,i_generator), & mask(1,1,d_hole1), mask(1,1,d_part1), & mask(1,1,d_hole2), mask(1,1,d_part2), & i_generator, 0 $params_post) endif if($do_mono_excitations)then - call $subroutine_monoexc(psi_generators(1,1,i_generator), & + call $subroutine_monoexc(psi_det_generators(1,1,i_generator), & mask(1,1,s_hole ), mask(1,1,s_part ), & i_generator, 0 $params_post) endif @@ -481,33 +481,33 @@ subroutine $subroutine($params_main) do k=1,N_int mask(k,ispin,s_hole) = & iand(generators_bitmask(k,ispin,s_hole,i_bitmask_gen), & - psi_generators(k,ispin,i_generator) ) + psi_det_generators(k,ispin,i_generator) ) mask(k,ispin,s_part) = & iand(generators_bitmask(k,ispin,s_part,i_bitmask_gen), & - not(psi_generators(k,ispin,i_generator)) ) + not(psi_det_generators(k,ispin,i_generator)) ) mask(k,ispin,d_hole1) = & iand(generators_bitmask(k,ispin,d_hole1,i_bitmask_gen), & - psi_generators(k,ispin,i_generator) ) + psi_det_generators(k,ispin,i_generator) ) mask(k,ispin,d_part1) = & iand(generators_bitmask(k,ispin,d_part1,i_bitmask_gen), & - not(psi_generators(k,ispin,i_generator)) ) + not(psi_det_generators(k,ispin,i_generator)) ) mask(k,ispin,d_hole2) = & iand(generators_bitmask(k,ispin,d_hole2,i_bitmask_gen), & - psi_generators(k,ispin,i_generator) ) + psi_det_generators(k,ispin,i_generator) ) mask(k,ispin,d_part2) = & iand(generators_bitmask(k,ispin,d_part2,i_bitmask_gen), & - not (psi_generators(k,ispin,i_generator)) ) + not (psi_det_generators(k,ispin,i_generator)) ) enddo enddo if($do_double_excitations)then - call $subroutine_diexc(psi_generators(1,1,i_generator), & + call $subroutine_diexc(psi_det_generators(1,1,i_generator), & mask(1,1,d_hole1), mask(1,1,d_part1), & mask(1,1,d_hole2), mask(1,1,d_part2), & i_generator, iproc $params_post) endif if($do_mono_excitations)then - call $subroutine_monoexc(psi_generators(1,1,i_generator), & + call $subroutine_monoexc(psi_det_generators(1,1,i_generator), & mask(1,1,s_hole ), mask(1,1,s_part ), & i_generator, iproc $params_post) endif diff --git a/src/Generators_CAS/generators.irp.f b/src/Generators_CAS/generators.irp.f index 1ea98eee..55562248 100644 --- a/src/Generators_CAS/generators.irp.f +++ b/src/Generators_CAS/generators.irp.f @@ -31,7 +31,7 @@ BEGIN_PROVIDER [ integer, N_det_generators ] call write_int(output_dets,N_det_generators,'Number of generators') END_PROVIDER -BEGIN_PROVIDER [ integer(bit_kind), psi_generators, (N_int,2,psi_det_size) ] +BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators, (N_int,2,psi_det_size) ] implicit none BEGIN_DOC ! For Single reference wave functions, the generator is the @@ -57,10 +57,10 @@ BEGIN_PROVIDER [ integer(bit_kind), psi_generators, (N_int,2,psi_det_size) ] if (good) then m = m+1 do k=1,N_int - psi_generators(k,1,m) = psi_det(k,1,i) - psi_generators(k,2,m) = psi_det(k,2,i) + psi_det_generators(k,1,m) = psi_det(k,1,i) + psi_det_generators(k,2,m) = psi_det(k,2,i) enddo -! call debug_det(psi_generators(1,1,m),N_int) +! call debug_det(psi_det_generators(1,1,m),N_int) endif enddo diff --git a/src/Generators_full/generators.irp.f b/src/Generators_full/generators.irp.f index 58fefef7..ff2650b2 100644 --- a/src/Generators_full/generators.irp.f +++ b/src/Generators_full/generators.irp.f @@ -34,7 +34,7 @@ BEGIN_PROVIDER [ integer, N_det_generators ] call write_int(output_dets,N_det_generators,'Number of generators') END_PROVIDER -BEGIN_PROVIDER [ integer(bit_kind), psi_generators, (N_int,2,psi_det_size) ] +BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators, (N_int,2,psi_det_size) ] implicit none BEGIN_DOC ! For Single reference wave functions, the generator is the @@ -43,8 +43,8 @@ BEGIN_PROVIDER [ integer(bit_kind), psi_generators, (N_int,2,psi_det_size) ] integer :: i, k do i=1,N_det do k=1,N_int - psi_generators(k,1,i) = psi_det_sorted(k,1,i) - psi_generators(k,2,i) = psi_det_sorted(k,2,i) + psi_det_generators(k,1,i) = psi_det_sorted(k,1,i) + psi_det_generators(k,2,i) = psi_det_sorted(k,2,i) enddo enddo @@ -58,7 +58,7 @@ BEGIN_PROVIDER [integer, degree_max_generators] integer :: i,degree degree_max_generators = 0 do i = 1, N_det_generators - call get_excitation_degree(HF_bitmask,psi_generators(1,1,i),degree,N_int) + call get_excitation_degree(HF_bitmask,psi_det_generators(1,1,i),degree,N_int) if(degree .gt. degree_max_generators)then degree_max_generators = degree endif diff --git a/src/Generators_restart/generators.irp.f b/src/Generators_restart/generators.irp.f index 19a6a9fd..f3db53c3 100644 --- a/src/Generators_restart/generators.irp.f +++ b/src/Generators_restart/generators.irp.f @@ -17,8 +17,8 @@ BEGIN_PROVIDER [ integer, N_det_generators ] END_PROVIDER - BEGIN_PROVIDER [ integer(bit_kind), psi_generators, (N_int,2,N_det_generators) ] -&BEGIN_PROVIDER [ double precision, psi_generators_coef, (N_det_generators,N_states) + BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators, (N_int,2,N_det_generators) ] +&BEGIN_PROVIDER [ double precision, psi_coef_generators, (N_det_generators,N_states) ] implicit none BEGIN_DOC ! read wf @@ -29,11 +29,11 @@ END_PROVIDER if(ifirst == 0)then do i=1,N_det_generators do k=1,N_int - psi_generators(k,1,i) = psi_det(k,1,i) - psi_generators(k,2,i) = psi_det(k,2,i) + psi_det_generators(k,1,i) = psi_det(k,1,i) + psi_det_generators(k,2,i) = psi_det(k,2,i) enddo do k = 1, N_states - psi_generators_coef(i,k) = psi_coef(i,k) + psi_coef_generators(i,k) = psi_coef(i,k) enddo enddo ifirst = 1 diff --git a/src/MRCC/mrcc_dress.irp.f b/src/MRCC/mrcc_dress.irp.f index 56257556..460365a4 100644 --- a/src/MRCC/mrcc_dress.irp.f +++ b/src/MRCC/mrcc_dress.irp.f @@ -20,7 +20,7 @@ subroutine mrcc_dress(delta_ij_sd_,Ndet_sd,i_generator,n_selected,det_buffer,Nin N_tq = 0 do i=1,N_selected - c_ref = connected_to_ref(det_buffer(1,1,i),psi_generators,Nint, & + c_ref = connected_to_ref(det_buffer(1,1,i),psi_det_generators,Nint, & i_generator,N_det_generators) if (c_ref /= 0) then diff --git a/src/Perturbation/perturbation_template.f b/src/Perturbation/perturbation_template.f index e289d409..9f50b636 100644 --- a/src/Perturbation/perturbation_template.f +++ b/src/Perturbation/perturbation_template.f @@ -25,7 +25,7 @@ subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,c ASSERT (N_st > 0) do i = 1,buffer_size - c_ref = connected_to_ref(buffer(1,1,i),psi_generators,Nint,i_generator,N_det_generators) + c_ref = connected_to_ref(buffer(1,1,i),psi_det_generators,Nint,i_generator,N_det_generators) if (c_ref /= 0) then cycle @@ -76,7 +76,7 @@ subroutine perturb_buffer_by_mono_$PERT(i_generator,buffer,buffer_size,e_2_pert_ ASSERT (N_st > 0) do i = 1,buffer_size - c_ref = connected_to_ref_by_mono(buffer(1,1,i),psi_generators,Nint,i_generator,N_det) + c_ref = connected_to_ref_by_mono(buffer(1,1,i),psi_det_generators,Nint,i_generator,N_det) if (c_ref /= 0) then cycle diff --git a/src/SingleRefMethod/generators.irp.f b/src/SingleRefMethod/generators.irp.f index 6ef30fa1..eeb1da70 100644 --- a/src/SingleRefMethod/generators.irp.f +++ b/src/SingleRefMethod/generators.irp.f @@ -9,19 +9,20 @@ BEGIN_PROVIDER [ integer, N_det_generators ] N_det_generators = 1 END_PROVIDER -BEGIN_PROVIDER [ integer(bit_kind), psi_generators, (N_int,2,psi_det_size) ] + BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators, (N_int,2,psi_det_size) ] +&BEGIN_PROVIDER [ integer(bit_kind), psi_coef_generators, (1,N_states) ] implicit none BEGIN_DOC ! For Single reference wave functions, the generator is the ! Hartree-Fock determinant END_DOC - psi_generators = 0_bit_kind + psi_det_generators = 0_bit_kind integer :: i,j,k integer :: degree do i=1,N_int - psi_generators(i,1,1) = HF_bitmask(i,1) - psi_generators(i,2,1) = HF_bitmask(i,2) + psi_det_generators(i,1,1) = HF_bitmask(i,1) + psi_det_generators(i,2,1) = HF_bitmask(i,2) enddo do j=1,N_det @@ -32,12 +33,8 @@ BEGIN_PROVIDER [ integer(bit_kind), psi_generators, (N_int,2,psi_det_size) ] endif end do - do j=2,k - psi_generators(:,:,j) = psi_det(:,:,j-1) - enddo - do j=k+1,N_det - psi_generators(:,:,j) = psi_det(:,:,j) - enddo + psi_det_generators(:,:,1) = psi_det(:,:,j) + psi_coef_generators(1,:) = psi_coef_generators(j,:) END_PROVIDER From 2204bc4aba27dd157929714bbb7b0f78c61288eb Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 2 Apr 2015 11:40:16 +0200 Subject: [PATCH 5/5] Corrected CAS-SD bug --- scripts/generate_h_apply.py | 4 +-- src/CAS_SD_selected/cas_sd.irp.f | 27 ++++++------------ src/Generators_CAS/generators.irp.f | 14 +++++---- src/Generators_full/generators.irp.f | 4 ++- src/Generators_restart/generators.irp.f | 4 +-- src/MRCC/mrcc_dress.irp.f | 38 ------------------------- src/Makefile.common | 3 ++ src/SingleRefMethod/generators.irp.f | 2 +- 8 files changed, 28 insertions(+), 68 deletions(-) diff --git a/scripts/generate_h_apply.py b/scripts/generate_h_apply.py index 36bd48aa..ad6a57ee 100755 --- a/scripts/generate_h_apply.py +++ b/scripts/generate_h_apply.py @@ -220,7 +220,7 @@ class H_apply(object): self.data["printout_always"] = """ do k=1,N_st - norm_psi(k) = norm_psi(k) + psi_coef(i_generator,k)*psi_coef(i_generator,k) + norm_psi(k) = norm_psi(k) + psi_coef_generators(i_generator,k)*psi_coef_generators(i_generator,k) delta_pt2(k) = pt2(k) - pt2_old(k) enddo """ @@ -275,7 +275,7 @@ class H_apply(object): if (select_max(i_generator) < selection_criterion_min*selection_criterion_factor) then !$ call omp_set_lock(lck) do k=1,N_st - norm_psi(k) = norm_psi(k) + psi_coef(i_generator,k)*psi_coef(i_generator,k) + norm_psi(k) = norm_psi(k) + psi_coef_generators(i_generator,k)*psi_coef_generators(i_generator,k) ! delta_pt2(k) = 0.d0 ! pt2_old(k) = 0.d0 ! pt2(k) = select_max(i_generator) diff --git a/src/CAS_SD_selected/cas_sd.irp.f b/src/CAS_SD_selected/cas_sd.irp.f index e3721ca7..e5e0497a 100644 --- a/src/CAS_SD_selected/cas_sd.irp.f +++ b/src/CAS_SD_selected/cas_sd.irp.f @@ -57,9 +57,9 @@ program full_ci ! Check that it is a CAS-SD logical :: in_cas - integer :: exc_max + integer :: exc_max, degree_min exc_max = 0 - print *, 'CAS determinants' + print *, 'CAS determinants : ', N_det_generators do i=1,N_det_generators do k=i,N_det_generators call get_excitation_degree(psi_det_generators(1,1,k),psi_det_generators(1,1,i),degree,N_int) @@ -71,25 +71,16 @@ program full_ci print *, 'Max excitation degree in the CAS :', exc_max do i=1,N_det in_cas = .False. + degree_min = 1000 do k=1,N_det_generators call get_excitation_degree(psi_det_generators(1,1,k),psi_det(1,1,i),degree,N_int) - if (degree == 0) then - in_cas = .True. - exit - endif + degree_min = min(degree_min,degree) enddo - if (.not.in_cas) then - do k=i,N_det - call get_excitation_degree(psi_det(1,1,k),psi_det(1,1,i),degree,N_int) - if (degree > exc_max+2) then - print *, 'Error : This is not a CAS-SD : ' - print *, 'CAS determinant:' - call debug_det(psi_det(1,1,i),N_int) - print *, 'Excited determinant:', degree - call debug_det(psi_det(1,1,k),N_int) - stop - endif - enddo + if (degree_min > 2) then + print *, 'Error : This is not a CAS-SD : ' + print *, 'Excited determinant:', degree_min + call debug_det(psi_det(1,1,k),N_int) + stop endif enddo end diff --git a/src/Generators_CAS/generators.irp.f b/src/Generators_CAS/generators.irp.f index 55562248..5dcfc1b0 100644 --- a/src/Generators_CAS/generators.irp.f +++ b/src/Generators_CAS/generators.irp.f @@ -15,9 +15,9 @@ BEGIN_PROVIDER [ integer, N_det_generators ] do k=1,N_int good = good .and. ( & iand(not(cas_bitmask(k,1,l)), psi_det(k,1,i)) == & - iand(not(cas_bitmask(k,1,l)), psi_det(k,1,1)) ) .and. ( & + iand(not(cas_bitmask(k,1,l)), HF_bitmask(k,1)) ) .and. ( & iand(not(cas_bitmask(k,2,l)), psi_det(k,2,i)) == & - iand(not(cas_bitmask(k,2,l)), psi_det(k,2,1)) ) + iand(not(cas_bitmask(k,2,l)), HF_bitmask(k,2)) ) enddo if (good) then exit @@ -31,7 +31,8 @@ BEGIN_PROVIDER [ integer, N_det_generators ] call write_int(output_dets,N_det_generators,'Number of generators') END_PROVIDER -BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators, (N_int,2,psi_det_size) ] + BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators, (N_int,2,psi_det_size) ] +&BEGIN_PROVIDER [ double precision, psi_coef_generators, (psi_det_size,N_states) ] implicit none BEGIN_DOC ! For Single reference wave functions, the generator is the @@ -44,11 +45,11 @@ BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators, (N_int,2,psi_det_size) ] do l=1,n_cas_bitmask good = .True. do k=1,N_int - good = good .and. ( & + good = good .and. ( & iand(not(cas_bitmask(k,1,l)), psi_det(k,1,i)) == & - iand(not(cas_bitmask(k,1,l)), psi_det(k,1,1)) ) .and. ( & + iand(not(cas_bitmask(k,1,l)), HF_bitmask(k,1)) .and. ( & iand(not(cas_bitmask(k,2,l)), psi_det(k,2,i)) == & - iand(not(cas_bitmask(k,2,l)), psi_det(k,2,1)) ) + iand(not(cas_bitmask(k,2,l)), HF_bitmask(k,2) )) ) enddo if (good) then exit @@ -60,6 +61,7 @@ BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators, (N_int,2,psi_det_size) ] psi_det_generators(k,1,m) = psi_det(k,1,i) psi_det_generators(k,2,m) = psi_det(k,2,i) enddo + psi_coef_generators(m,:) = psi_coef(m,:) ! call debug_det(psi_det_generators(1,1,m),N_int) endif enddo diff --git a/src/Generators_full/generators.irp.f b/src/Generators_full/generators.irp.f index ff2650b2..4d261acd 100644 --- a/src/Generators_full/generators.irp.f +++ b/src/Generators_full/generators.irp.f @@ -34,7 +34,8 @@ BEGIN_PROVIDER [ integer, N_det_generators ] call write_int(output_dets,N_det_generators,'Number of generators') END_PROVIDER -BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators, (N_int,2,psi_det_size) ] + BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators, (N_int,2,psi_det_size) ] +&BEGIN_PROVIDER [ double precision, psi_coef_generators, (psi_det_size,N_states) ] implicit none BEGIN_DOC ! For Single reference wave functions, the generator is the @@ -46,6 +47,7 @@ BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators, (N_int,2,psi_det_size) ] psi_det_generators(k,1,i) = psi_det_sorted(k,1,i) psi_det_generators(k,2,i) = psi_det_sorted(k,2,i) enddo + psi_coef_generators(i,:) = psi_coef_sorted(i,:) enddo END_PROVIDER diff --git a/src/Generators_restart/generators.irp.f b/src/Generators_restart/generators.irp.f index f3db53c3..2e0bc375 100644 --- a/src/Generators_restart/generators.irp.f +++ b/src/Generators_restart/generators.irp.f @@ -17,8 +17,8 @@ BEGIN_PROVIDER [ integer, N_det_generators ] END_PROVIDER - BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators, (N_int,2,N_det_generators) ] -&BEGIN_PROVIDER [ double precision, psi_coef_generators, (N_det_generators,N_states) ] + BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators, (N_int,2,psi_det_size) ] +&BEGIN_PROVIDER [ double precision, psi_coef_generators, (psi_det_size,N_states) ] implicit none BEGIN_DOC ! read wf diff --git a/src/MRCC/mrcc_dress.irp.f b/src/MRCC/mrcc_dress.irp.f index 460365a4..e5c8d233 100644 --- a/src/MRCC/mrcc_dress.irp.f +++ b/src/MRCC/mrcc_dress.irp.f @@ -53,22 +53,6 @@ subroutine mrcc_dress(delta_ij_sd_,Ndet_sd,i_generator,n_selected,det_buffer,Nin double precision :: haj double precision :: f(N_states) - -! call i_h_j(psi_det(1,1,1), psi_det(1,1,64),Nint,hka) -! call debug_det(psi_det(1,1,1), N_int) -! call debug_det(psi_det(1,1,64), N_int) - double precision :: phase - integer :: exc(0:2,2,2) -! call get_excitation(psi_det(1,1,1),psi_det(1,1,64),exc,degree(1),phase,Nint) - integer :: h1, p1, h2, p2, s1, s2 -! call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) -! print *, hka -! print *, h1, p1, h2, p2 -! print *, s1, s2 -! pause - - - do i=1,N_tq call get_excitation_degree_vector(psi_sd,tq(1,1,i),degree,Nint,Ndet_sd,idx) call i_h_j(tq(1,1,i),tq(1,1,i),Nint,haa) @@ -83,21 +67,6 @@ subroutine mrcc_dress(delta_ij_sd_,Ndet_sd,i_generator,n_selected,det_buffer,Nin delta_ij_sd_(idx(k), idx(j),m) += haj*hka* f(m) delta_ij_sd_(idx(j), idx(k),m) += haj*hka* f(m) enddo - call get_excitation(tq(1,1,i),psi_sd(1,1,idx(j)),exc,degree(1),phase,Nint) - call decode_exc(exc,degree(1),h1,p1,h2,p2,s1,s2) - if ( (h1 == 1).and. & - (p1 == 6).and. & - (h2 == 1).and. & - (p2 == 6).and. & - (s1 == 1).and. & - (s2 == 2) ) then - call debug_det(tq(1,1,i), N_int) - call debug_det(psi_sd(1,1,idx(j)), N_int) - print *, haj - pause - endif - - enddo enddo enddo @@ -196,13 +165,6 @@ END_PROVIDER exit endif enddo -! if(i_state < min(N_states_diag,N_det))then -! print *, 'pb with the number of states' -! print *, 'i_state = ',i_state -! print *, 'N_states_diag ',N_states_diag -! print *,'stopping ...' -! stop -! endif deallocate(eigenvectors,eigenvalues) endif diff --git a/src/Makefile.common b/src/Makefile.common index 9ecc93c0..eddaa481 100644 --- a/src/Makefile.common +++ b/src/Makefile.common @@ -8,6 +8,9 @@ default: veryclean: clean_modules.sh +clean: + IN_MAKE=1 make clean + else # Called by scripts/build_module.sh default: all .gitignore diff --git a/src/SingleRefMethod/generators.irp.f b/src/SingleRefMethod/generators.irp.f index eeb1da70..ce71f996 100644 --- a/src/SingleRefMethod/generators.irp.f +++ b/src/SingleRefMethod/generators.irp.f @@ -10,7 +10,7 @@ BEGIN_PROVIDER [ integer, N_det_generators ] END_PROVIDER BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators, (N_int,2,psi_det_size) ] -&BEGIN_PROVIDER [ integer(bit_kind), psi_coef_generators, (1,N_states) ] +&BEGIN_PROVIDER [ double precision, psi_coef_generators, (psi_det_size,N_states) ] implicit none BEGIN_DOC ! For Single reference wave functions, the generator is the