From fe3c2bb875f822eda66af4c411b239c72e94338e Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 3 Apr 2015 14:26:14 +0200 Subject: [PATCH 01/12] Minor changes --- scripts/ezfio_with_default.py | 1 + src/CISD/cisd.irp.f | 1 + src/MRCC/H_apply.irp.f | 10 +- src/MRCC/mrcc_dress.irp.f | 197 ++++++++++++++-------------------- src/MRCC/mrcc_utils.irp.f | 132 +++++++++++++++++++++++ 5 files changed, 223 insertions(+), 118 deletions(-) diff --git a/scripts/ezfio_with_default.py b/scripts/ezfio_with_default.py index 1208f4b7..63666abd 100755 --- a/scripts/ezfio_with_default.py +++ b/scripts/ezfio_with_default.py @@ -94,6 +94,7 @@ END_PROVIDER file = open(filename,'r') lines = file.readlines() file.close() + k=-1 # Search directory for k,line in enumerate(lines): if line[0] != ' ': diff --git a/src/CISD/cisd.irp.f b/src/CISD/cisd.irp.f index 0b6bc9fd..6d310d95 100644 --- a/src/CISD/cisd.irp.f +++ b/src/CISD/cisd.irp.f @@ -13,6 +13,7 @@ program cisd print *, 'E_corr = ',CI_electronic_energy(i) - ref_bitmask_energy enddo + call save_wavefunction ! call CISD_SC2(psi_det,psi_coef,eigvalues,size(psi_coef,1),N_det,N_states,N_int) ! do i = 1, N_states ! print*,'eigvalues(i) = ',eigvalues(i) diff --git a/src/MRCC/H_apply.irp.f b/src/MRCC/H_apply.irp.f index 4f284112..a8519c25 100644 --- a/src/MRCC/H_apply.irp.f +++ b/src/MRCC/H_apply.irp.f @@ -2,13 +2,13 @@ use bitmasks BEGIN_SHELL [ /usr/bin/env python ] from generate_h_apply import * -s = H_apply("mrcc") +s = H_apply("mrcc_simple") 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["keys_work"] = "call mrcc_dress_simple(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"] += """ @@ -19,7 +19,13 @@ s.data["finalization"] = "" s.data["copy_buffer"] = "" s.data["generate_psi_guess"] = "" s.data["size_max"] = "256" +print s + + + +s.data["subroutine"] = "H_apply_mrcc" +s.data["keys_work"] = "call mrcc_dress(delta_ij_sd_,Ndet_sd,i_generator,key_idx,keys_out,N_int,iproc)" print s END_SHELL diff --git a/src/MRCC/mrcc_dress.irp.f b/src/MRCC/mrcc_dress.irp.f index e5c8d233..d9c8978c 100644 --- a/src/MRCC/mrcc_dress.irp.f +++ b/src/MRCC/mrcc_dress.irp.f @@ -1,3 +1,84 @@ +subroutine mrcc_dress_simple(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,m + integer :: new_size + logical :: is_in_wavefunction + integer :: degree(psi_det_size) + integer :: idx(0:psi_det_size) + 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_det_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 + + ! Compute / (E0 - Haa) + double precision :: hka, haa + double precision :: haj + double precision :: f(N_states) + + 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) + 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 + enddo + enddo + enddo +end + + + + + + + + subroutine mrcc_dress(delta_ij_sd_,Ndet_sd,i_generator,n_selected,det_buffer,Nint,iproc) use bitmasks implicit none @@ -72,119 +153,3 @@ subroutine mrcc_dress(delta_ij_sd_,Ndet_sd,i_generator,n_selected,det_buffer,Nin 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 - 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 2c5dc811..01e566ae 100644 --- a/src/MRCC/mrcc_utils.irp.f +++ b/src/MRCC/mrcc_utils.irp.f @@ -148,3 +148,135 @@ BEGIN_PROVIDER [ double precision, lambda_mrcc, (psi_det_size,n_states) ] enddo END_PROVIDER + + +BEGIN_PROVIDER [ character*(32), dressing_type ] + implicit none + BEGIN_DOC + ! [ Simple | MRCC ] + END_DOC + dressing_type = "MRCC" +END_PROVIDER + +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 + if (dressing_type == "MRCC") then + call H_apply_mrcc(delta_ij_sd,N_det_sd) + else if (dressing_type == "Simple") then + call H_apply_mrcc_simple(delta_ij_sd,N_det_sd) + else + print *, irp_here + stop 'dressing' + endif +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 + 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 + From dd1a51af635bf8da51e8a22405acb991f87297d0 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 9 Apr 2015 10:57:15 +0200 Subject: [PATCH 02/12] CAS-SD and CAS-SD selected --- .../ASSUMPTIONS.rst | 0 src/{CAS_SD_selected => CAS_SD}/H_apply.irp.f | 7 +- src/{CAS_SD_selected => CAS_SD}/Makefile | 0 .../NEEDED_MODULES | 0 src/{CAS_SD_selected => CAS_SD}/README.rst | 0 src/{CAS_SD_selected => CAS_SD}/cas_sd.irp.f | 4 +- src/CAS_SD/cas_sd_selected.irp.f | 86 +++++++++++++++++++ src/{CAS_SD_selected => CAS_SD}/options.irp.f | 0 8 files changed, 93 insertions(+), 4 deletions(-) rename src/{CAS_SD_selected => CAS_SD}/ASSUMPTIONS.rst (100%) rename src/{CAS_SD_selected => CAS_SD}/H_apply.irp.f (68%) rename src/{CAS_SD_selected => CAS_SD}/Makefile (100%) rename src/{CAS_SD_selected => CAS_SD}/NEEDED_MODULES (100%) rename src/{CAS_SD_selected => CAS_SD}/README.rst (100%) rename src/{CAS_SD_selected => CAS_SD}/cas_sd.irp.f (96%) create mode 100644 src/CAS_SD/cas_sd_selected.irp.f rename src/{CAS_SD_selected => CAS_SD}/options.irp.f (100%) diff --git a/src/CAS_SD_selected/ASSUMPTIONS.rst b/src/CAS_SD/ASSUMPTIONS.rst similarity index 100% rename from src/CAS_SD_selected/ASSUMPTIONS.rst rename to src/CAS_SD/ASSUMPTIONS.rst diff --git a/src/CAS_SD_selected/H_apply.irp.f b/src/CAS_SD/H_apply.irp.f similarity index 68% rename from src/CAS_SD_selected/H_apply.irp.f rename to src/CAS_SD/H_apply.irp.f index 76c046ae..85b2f1c6 100644 --- a/src/CAS_SD_selected/H_apply.irp.f +++ b/src/CAS_SD/H_apply.irp.f @@ -2,11 +2,14 @@ use bitmasks BEGIN_SHELL [ /usr/bin/env python ] from generate_h_apply import * -s = H_apply("FCI") +s = H_apply("CAS_SD") +print s + +s = H_apply("CAS_SD_selected") s.set_selection_pt2("epstein_nesbet_2x2") print s -s = H_apply("FCI_PT2") +s = H_apply("CAS_SD_PT2") s.set_perturbation("epstein_nesbet_2x2") print s diff --git a/src/CAS_SD_selected/Makefile b/src/CAS_SD/Makefile similarity index 100% rename from src/CAS_SD_selected/Makefile rename to src/CAS_SD/Makefile diff --git a/src/CAS_SD_selected/NEEDED_MODULES b/src/CAS_SD/NEEDED_MODULES similarity index 100% rename from src/CAS_SD_selected/NEEDED_MODULES rename to src/CAS_SD/NEEDED_MODULES diff --git a/src/CAS_SD_selected/README.rst b/src/CAS_SD/README.rst similarity index 100% rename from src/CAS_SD_selected/README.rst rename to src/CAS_SD/README.rst diff --git a/src/CAS_SD_selected/cas_sd.irp.f b/src/CAS_SD/cas_sd.irp.f similarity index 96% rename from src/CAS_SD_selected/cas_sd.irp.f rename to src/CAS_SD/cas_sd.irp.f index e5e0497a..bf772b1c 100644 --- a/src/CAS_SD_selected/cas_sd.irp.f +++ b/src/CAS_SD/cas_sd.irp.f @@ -29,7 +29,7 @@ program full_ci 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) + call H_apply_CAS_SD(pt2, norm_pert, H_pert_diag, N_st) PROVIDE psi_coef PROVIDE psi_det @@ -60,7 +60,7 @@ program full_ci integer :: exc_max, degree_min exc_max = 0 print *, 'CAS determinants : ', N_det_generators - do i=1,N_det_generators + do i=1,min(N_det_generators,10) 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) exc_max = max(exc_max,degree) diff --git a/src/CAS_SD/cas_sd_selected.irp.f b/src/CAS_SD/cas_sd_selected.irp.f new file mode 100644 index 00000000..1f1e0d24 --- /dev/null +++ b/src/CAS_SD/cas_sd_selected.irp.f @@ -0,0 +1,86 @@ +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_CAS_SD_selected(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 + + ! Check that it is a CAS-SD + logical :: in_cas + integer :: exc_max, degree_min + exc_max = 0 + print *, 'CAS determinants : ', N_det_generators + do i=1,min(N_det_generators,10) + 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) + exc_max = max(exc_max,degree) + enddo + 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. + 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) + degree_min = min(degree_min,degree) + 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/CAS_SD_selected/options.irp.f b/src/CAS_SD/options.irp.f similarity index 100% rename from src/CAS_SD_selected/options.irp.f rename to src/CAS_SD/options.irp.f From 4ea60695ea90aded67b015161828f7d7ec6cb3bb Mon Sep 17 00:00:00 2001 From: Thomas Applencourt Date: Thu, 9 Apr 2015 14:03:24 +0200 Subject: [PATCH 03/12] Add Patch requirement for OPEM --- INSTALL.md | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/INSTALL.md b/INSTALL.md index 712e3d8f..47319c9c 100644 --- a/INSTALL.md +++ b/INSTALL.md @@ -6,8 +6,9 @@ * m4 * GNU make * Fortran compiler (ifort or gfortran are tested) -* Python >= 2.7 +* Python >= 2.6 * Bash +* Patch (for opam) ## Standard installation From 102bbb0b4fe4c8d057103c28dd6145da79894595 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 9 Apr 2015 21:46:28 +0200 Subject: [PATCH 04/12] Lots of minor fixes --- data/ezfio_defaults | 6 + scripts/ezfio_interface.py | 6 +- scripts/ezfio_with_default.py | 7 +- scripts/generate_h_apply.py | 2 +- scripts/qp_convert.py | 203 ++++++++-------------------- src/BiInts/map_integrals.irp.f | 6 +- src/Bitmask/bitmasks.irp.f | 2 - src/CAS_SD/README.rst | 6 +- src/CAS_SD/cas_sd.irp.f | 10 +- src/CAS_SD/cas_sd_selected.irp.f | 10 +- src/Dets/H_apply_template.f | 25 ++-- src/Full_CI/EZFIO.cfg | 14 +- src/Generators_CAS/generators.irp.f | 1 - src/MRCC/NEEDED_MODULES | 2 +- src/NEEDED_MODULES | 2 +- src/Utils/LinearAlgebra.irp.f | 3 + 16 files changed, 121 insertions(+), 184 deletions(-) diff --git a/data/ezfio_defaults b/data/ezfio_defaults index 5da4d813..347e3751 100644 --- a/data/ezfio_defaults +++ b/data/ezfio_defaults @@ -32,6 +32,12 @@ full_ci do_pt2_end True var_pt2_ratio 0.75 +cas_sd + n_det_max_cas_sd 100000 + pt2_max 1.e-4 + do_pt2_end True + var_pt2_ratio 0.75 + all_singles n_det_max_fci 50000 pt2_max 1.e-8 diff --git a/scripts/ezfio_interface.py b/scripts/ezfio_interface.py index 28203b52..7c196cf8 100755 --- a/scripts/ezfio_interface.py +++ b/scripts/ezfio_interface.py @@ -176,7 +176,7 @@ def get_dict_config_file(config_file_path, module_lower): d[pvd][option] = d_default[option] # If interface is output we need a default value information - if d[pvd]["interface"] == "output": + if d[pvd]["interface"] == "input": try: d[pvd]["default"] = config_file.get(section, "default") except ConfigParser.NoOptionError: @@ -210,8 +210,8 @@ def create_ezfio_provider(dict_ezfio_cfg): ez_p.set_ezfio_dir(dict_info['ezfio_dir']) ez_p.set_ezfio_name(dict_info['ezfio_name']) ez_p.set_default(dict_info['default']) - ez_p.set_output("output_%s" % dict_info['ezfio_dir']) + dict_code_provider[provider_name] = str(ez_p) return dict_code_provider @@ -239,7 +239,7 @@ def save_ezfio_provider(path_head, dict_code_provider): "! from file {0}/EZFIO.cfg\n".format(path_head) + \ "\n" for provider_name, code in dict_code_provider.iteritems(): - output += code + "\n" + output += str(code) + "\n" if output != old_output: with open(path, "w") as f: diff --git a/scripts/ezfio_with_default.py b/scripts/ezfio_with_default.py index 63666abd..75fef0b6 100755 --- a/scripts/ezfio_with_default.py +++ b/scripts/ezfio_with_default.py @@ -112,16 +112,19 @@ END_PROVIDER break v = buffer[1] name = self.name + true = True + false= False try: v_eval = eval(v) + except: + v = "call ezfio_get_%(v)s(%(name)s)"%locals() + else: if type(v_eval) == bool: v = '.%s.'%(v) elif type(v_eval) == float: v = v.replace('e','d') v = v.replace('E','D') v = "%(name)s = %(v)s"%locals() - except: - v = "call ezfio_get_%(v)s(%(name)s)"%locals() self.default = v diff --git a/scripts/generate_h_apply.py b/scripts/generate_h_apply.py index ad6a57ee..280c9f72 100755 --- a/scripts/generate_h_apply.py +++ b/scripts/generate_h_apply.py @@ -53,7 +53,7 @@ class H_apply(object): !$OMP N_elec_in_key_hole_2,ia_ja_pairs) & !$OMP SHARED(key_in,N_int,elec_num_tab,mo_tot_num, & !$OMP hole_1, particl_1, hole_2, particl_2, & - !$OMP elec_alpha_num,i_generator)""" + !$OMP elec_alpha_num,i_generator) FIRSTPRIVATE(iproc)""" s["omp_end_parallel"] = "!$OMP END PARALLEL" s["omp_master"] = "!$OMP MASTER" s["omp_end_master"] = "!$OMP END MASTER" diff --git a/scripts/qp_convert.py b/scripts/qp_convert.py index 54a6d2da..ab008e9e 100755 --- a/scripts/qp_convert.py +++ b/scripts/qp_convert.py @@ -53,16 +53,68 @@ def write_ezfioFile(res,filename): basis = res.uncontracted_basis geom = res.geometry + res.clean_contractions() + # AO Basis + import string + at = [] + num_prim = [] + magnetic_number = [] + angular_number = [] + power_x = [] + power_y = [] + power_z = [] + coefficient = [] + exponent = [] + res.convert_to_cartesian() + for b in res.basis: + c = b.center + for i,atom in enumerate(res.geometry): + if atom.coord == c: + at.append(i+1) + num_prim.append(len(b.prim)) + s = b.sym + power_x.append( string.count(s,"x") ) + power_y.append( string.count(s,"y") ) + power_z.append( string.count(s,"z") ) + coefficient.append( b.coef ) + exponent.append( [ p.expo for p in b.prim ] ) + ezfio.set_ao_basis_ao_num(len(res.basis)) + ezfio.set_ao_basis_ao_nucl(at) + ezfio.set_ao_basis_ao_prim_num(num_prim) + ezfio.set_ao_basis_ao_power(power_x+power_y+power_z) + prim_num_max = ezfio.get_ao_basis_ao_prim_num_max() + len_res_basis = len(res.basis) + for i in range(len(res.basis)): + coefficient[i] += [ 0. for j in range(len(coefficient[i]),prim_num_max) ] + exponent[i] += [ 0. for j in range(len(exponent[i]),prim_num_max) ] + coefficient = reduce(lambda x, y: x+y, coefficient, []) + exponent = reduce(lambda x, y: x+y, exponent, []) + coef = [] + expo = [] + for i in range(prim_num_max): + for j in range(i,len(coefficient),prim_num_max): + coef.append ( coefficient[j] ) + expo.append ( exponent[j] ) + ezfio.set_ao_basis_ao_coef(coef) + ezfio.set_ao_basis_ao_expo(expo) + ezfio.set_ao_basis_ao_basis("Read by resultsFile") + + # MO MoTag = res.determinants_mo_type ezfio.set_mo_basis_mo_label('Orthonormalized') MO_type = MoTag - allMOs = res.uncontracted_mo_sets[MO_type] + allMOs = res.mo_sets[MO_type] - closed = [ (allMOs[i].eigenvalue,i) for i in res.closed_mos ] - active = [ (allMOs[i].eigenvalue,i) for i in res.active_mos ] - virtual =[ (allMOs[i].eigenvalue,i) for i in res.virtual_mos ] + try: + closed = [ (allMOs[i].eigenvalue,i) for i in res.closed_mos ] + active = [ (allMOs[i].eigenvalue,i) for i in res.active_mos ] + virtual =[ (allMOs[i].eigenvalue,i) for i in res.virtual_mos ] + except: + closed = [] + virtual = [] + active = [ (allMOs[i].eigenvalue,i) for i in range(len(allMOs)) ] # closed.sort() # active.sort() @@ -111,117 +163,6 @@ def write_ezfioFile(res,filename): while len(MoMatrix) < len(MOs[0].vector)**2: MoMatrix.append(0.) - ezfio.set_mo_basis_mo_tot_num(mo_tot_num) - ezfio.set_mo_basis_mo_occ(OccNum) - - - res.clean_contractions() - # AO Basis - import string - at = [] - num_prim = [] - magnetic_number = [] - angular_number = [] - power_x = [] - power_y = [] - power_z = [] - coefficient = [] - exponent = [] - res.convert_to_cartesian() - for b in res.basis: - c = b.center - for i,atom in enumerate(res.geometry): - if atom.coord == c: - at.append(i+1) - num_prim.append(len(b.prim)) - s = b.sym - power_x.append( string.count(s,"x") ) - power_y.append( string.count(s,"y") ) - power_z.append( string.count(s,"z") ) - coefficient.append( b.coef ) - exponent.append( [ p.expo for p in b.prim ] ) - ezfio.set_ao_basis_ao_num(len(res.basis)) - ezfio.set_ao_basis_ao_nucl(at) - ezfio.set_ao_basis_ao_prim_num(num_prim) - ezfio.set_ao_basis_ao_power(power_x+power_y+power_z) - prim_num_max = ezfio.get_ao_basis_ao_prim_num_max() - len_res_basis = len(res.basis) - for i in range(len(res.basis)): - coefficient[i] += [ 0. for j in range(len(coefficient[i]),prim_num_max) ] - exponent[i] += [ 0. for j in range(len(exponent[i]),prim_num_max) ] - coefficient = reduce(lambda x, y: x+y, coefficient, []) - exponent = reduce(lambda x, y: x+y, exponent, []) - coef = [] - expo = [] - for i in range(prim_num_max): - for j in range(i,len(coefficient),prim_num_max): - coef.append ( coefficient[j] ) - expo.append ( exponent[j] ) - ezfio.set_ao_basis_ao_coef(coef) - ezfio.set_ao_basis_ao_expo(expo) - ezfio.set_ao_basis_ao_basis("Read by resultsFile") - -# Apply threshold to determinants - if len(res.determinants) == 1: - sorted_determinants = [ (-1.,1.,res.determinants[0]) ] - else: - sorted_determinants = [] - for i,j in zip(res.det_coefficients[0],res.determinants): - sorted_determinants.append((-abs(i),i,j)) - sorted_determinants.sort() - norm = 0.0 - for length, (a,b,c) in enumerate(sorted_determinants): - if -a < det_threshold: - length -=1 - break - norm += a**2 - norm = sqrt(norm) - length += 1 - for i in xrange(length): - a = sorted_determinants[i] - sorted_determinants[i] = (a[0],a[1]/norm,a[2]) - sorted_determinants = sorted_determinants[:length] - -# MOs - mo_tot_num = len(res.mo_sets[MoTag]) - closed_mos = res.closed_mos - active_mos = res.active_mos - virtual_mos = res.virtual_mos - to_remove = [] - to_add = [] - for i in active_mos: - found = False - for (a,b,c) in sorted_determinants: - if i in c['alpha']+c['beta']: - found = True - break - if not found: - to_remove.append(i) - to_add.append(i) - virtual_mos = to_add + virtual_mos - for i in active_mos: - always = True - for (a,b,c) in sorted_determinants: - if not (i in c['alpha'] and i in c['beta']): - always = False - break - if always: - to_remove.append(i) - closed_mos.append(i) - for i in to_remove: - active_mos.remove(i) - - - MOindices = closed_mos + active_mos + virtual_mos - while len(MOindices) < mo_tot_num: - MOindices.append(len(MOindices)) - MOmap = list(MOindices) - for i in range(len(MOindices)): - MOmap[i] = MOindices.index(i) - - - ezfio.set_mo_basis_mo_tot_num(mo_tot_num) - mo = [] for i in MOindices: mo.append(res.mo_sets[MoTag][i]) @@ -234,35 +175,11 @@ def write_ezfioFile(res,filename): while len(mo) < mo_tot_num: mo.append(newmo) Energies = [ m.eigenvalue for m in mo ] - - if res.occ_num is not None: - OccNum = [] - for i in MOindices: - OccNum.append(res.occ_num[MoTag][i]) - - while len(OccNum) < mo_tot_num: - OccNum.append(0.) - ezfio.set_mo_basis_mo_occ(OccNum) - - cls = [ "v" for i in mo ] - for i in closed_mos: - cls[MOmap[i]] = 'c' - for i in active_mos: - cls[MOmap[i]] = 'a' - - sym0 = [ i.sym for i in res.mo_sets[MoTag] ] - sym = [ i.sym for i in res.mo_sets[MoTag] ] - for i in xrange(len(sym)): - sym[MOmap[i]] = sym0[i] - - MoMatrix = [] - for m in mo: - for coef in m.vector: - MoMatrix.append(coef) - while len(MoMatrix) < len(mo[0].vector)**2: - MoMatrix.append(0.) + + ezfio.set_mo_basis_mo_tot_num(mo_tot_num) + ezfio.set_mo_basis_mo_occ(OccNum) ezfio.set_mo_basis_mo_coef(MoMatrix) - del MoMatrix + diff --git a/src/BiInts/map_integrals.irp.f b/src/BiInts/map_integrals.irp.f index ae50207e..3ebbf6fe 100644 --- a/src/BiInts/map_integrals.irp.f +++ b/src/BiInts/map_integrals.irp.f @@ -368,11 +368,11 @@ subroutine get_mo_bielec_integrals_existing_ik(j,l,sze,out_array,map) enddo logical :: integral_is_in_map - if (cache_key_kind == 8) then + if (key_kind == 8) then call i8radix_sort(hash,iorder,kk,-1) - else if (cache_key_kind == 4) then + else if (key_kind == 4) then call iradix_sort(hash,iorder,kk,-1) - else if (cache_key_kind == 2) then + else if (key_kind == 2) then call i2radix_sort(hash,iorder,kk,-1) endif diff --git a/src/Bitmask/bitmasks.irp.f b/src/Bitmask/bitmasks.irp.f index e0e04614..2d044ca5 100644 --- a/src/Bitmask/bitmasks.irp.f +++ b/src/Bitmask/bitmasks.irp.f @@ -123,7 +123,6 @@ BEGIN_PROVIDER [ integer(bit_kind), generators_bitmask, (N_int,2,6,N_generators_ call ezfio_has_bitmasks_generators(exists) if (exists) then - print*,'EXIST !!' call ezfio_get_bitmasks_generators(generators_bitmask) else integer :: k, ispin @@ -181,7 +180,6 @@ BEGIN_PROVIDER [ integer(bit_kind), cas_bitmask, (N_int,2,N_cas_bitmask) ] PROVIDE ezfio_filename call ezfio_has_bitmasks_cas(exists) - print*,'exists = ',exists if (exists) then call ezfio_get_bitmasks_cas(cas_bitmask) else diff --git a/src/CAS_SD/README.rst b/src/CAS_SD/README.rst index f66fe229..71f828ab 100644 --- a/src/CAS_SD/README.rst +++ b/src/CAS_SD/README.rst @@ -2,4 +2,8 @@ CAS_SD_selected Module ====================== -Selected CAS + SD module +Selected CAS + SD module. + +1) Set the different MO classes using the ``qp_set_mo_class`` command +2) Run the selected CAS+SD program + diff --git a/src/CAS_SD/cas_sd.irp.f b/src/CAS_SD/cas_sd.irp.f index bf772b1c..144eec88 100644 --- a/src/CAS_SD/cas_sd.irp.f +++ b/src/CAS_SD/cas_sd.irp.f @@ -11,12 +11,12 @@ program full_ci pt2 = 1.d0 diag_algorithm = "Lapack" - if (N_det > n_det_max_fci) then + if (N_det > n_det_max_cas_sd) then call diagonalize_CI call save_wavefunction psi_det = psi_det_sorted psi_coef = psi_coef_sorted - N_det = n_det_max_fci + N_det = n_det_max_cas_sd soft_touch N_det psi_det psi_coef call diagonalize_CI call save_wavefunction @@ -28,17 +28,17 @@ program full_ci print *, '-----' endif - do while (N_det < n_det_max_fci.and.maxval(abs(pt2(1:N_st))) > pt2_max) + do while (N_det < n_det_max_cas_sd.and.maxval(abs(pt2(1:N_st))) > pt2_max) call H_apply_CAS_SD(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 + if (N_det > n_det_max_cas_sd) then psi_det = psi_det_sorted psi_coef = psi_coef_sorted - N_det = n_det_max_fci + N_det = n_det_max_cas_sd soft_touch N_det psi_det psi_coef endif call diagonalize_CI diff --git a/src/CAS_SD/cas_sd_selected.irp.f b/src/CAS_SD/cas_sd_selected.irp.f index 1f1e0d24..fc9d4dd2 100644 --- a/src/CAS_SD/cas_sd_selected.irp.f +++ b/src/CAS_SD/cas_sd_selected.irp.f @@ -11,12 +11,12 @@ program full_ci pt2 = 1.d0 diag_algorithm = "Lapack" - if (N_det > n_det_max_fci) then + if (N_det > n_det_max_cas_sd) then call diagonalize_CI call save_wavefunction psi_det = psi_det_sorted psi_coef = psi_coef_sorted - N_det = n_det_max_fci + N_det = n_det_max_cas_sd soft_touch N_det psi_det psi_coef call diagonalize_CI call save_wavefunction @@ -28,17 +28,17 @@ program full_ci print *, '-----' endif - do while (N_det < n_det_max_fci.and.maxval(abs(pt2(1:N_st))) > pt2_max) + do while (N_det < n_det_max_cas_sd.and.maxval(abs(pt2(1:N_st))) > pt2_max) call H_apply_CAS_SD_selected(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 + if (N_det > n_det_max_cas_sd) then psi_det = psi_det_sorted psi_coef = psi_coef_sorted - N_det = n_det_max_fci + N_det = n_det_max_cas_sd soft_touch N_det psi_det psi_coef endif call diagonalize_CI diff --git a/src/Dets/H_apply_template.f b/src/Dets/H_apply_template.f index e44562e7..a9a282ae 100644 --- a/src/Dets/H_apply_template.f +++ b/src/Dets/H_apply_template.f @@ -1,4 +1,4 @@ -subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_generator, iproc $parameters ) +subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_generator, iproc_in $parameters ) use omp_lib use bitmasks implicit none @@ -14,7 +14,7 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_gene integer(bit_kind),allocatable :: keys_out(:,:,:) integer(bit_kind), intent(in) :: hole_1(N_int,2), particl_1(N_int,2) integer(bit_kind), intent(in) :: hole_2(N_int,2), particl_2(N_int,2) - integer, intent(in) :: iproc + integer, intent(in) :: iproc_in integer(bit_kind), allocatable :: hole_save(:,:) integer(bit_kind), allocatable :: key(:,:),hole(:,:), particle(:,:) integer(bit_kind), allocatable :: hole_tmp(:,:), particle_tmp(:,:) @@ -30,6 +30,7 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_gene integer, allocatable :: ia_ja_pairs(:,:,:) integer, allocatable :: ib_jb_pairs(:,:) double precision :: diag_H_mat_elem + integer :: iproc integer(omp_lock_kind), save :: lck, ifirst=0 if (ifirst == 0) then !$ call omp_init_lock(lck) @@ -38,12 +39,13 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_gene logical :: check_double_excitation check_double_excitation = .True. - + iproc = iproc_in $initialization $omp_parallel +!$ iproc = omp_get_thread_num() allocate (keys_out(N_int,2,size_max), hole_save(N_int,2), & key(N_int,2),hole(N_int,2), particle(N_int,2), hole_tmp(N_int,2),& particle_tmp(N_int,2), occ_particle(N_int*bit_kind_size,2), & @@ -248,7 +250,7 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_gene $finalization end -subroutine $subroutine_monoexc(key_in, hole_1,particl_1,i_generator,iproc $parameters ) +subroutine $subroutine_monoexc(key_in, hole_1,particl_1,i_generator,iproc_in $parameters ) use omp_lib use bitmasks implicit none @@ -262,7 +264,7 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,i_generator,iproc $param integer ,intent(in) :: i_generator integer(bit_kind),intent(in) :: key_in(N_int,2) integer(bit_kind),intent(in) :: hole_1(N_int,2), particl_1(N_int,2) - integer, intent(in) :: iproc + integer, intent(in) :: iproc_in integer(bit_kind),allocatable :: keys_out(:,:,:) integer(bit_kind),allocatable :: hole_save(:,:) integer(bit_kind),allocatable :: key(:,:),hole(:,:), particle(:,:) @@ -281,8 +283,11 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,i_generator,iproc $param logical, allocatable :: array_pairs(:,:) double precision :: diag_H_mat_elem integer(omp_lock_kind), save :: lck, ifirst=0 + integer :: iproc logical :: check_double_excitation + iproc = iproc_in + check_double_excitation = .True. $check_double_excitation @@ -295,6 +300,7 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,i_generator,iproc $param $initialization $omp_parallel +!$ iproc = omp_get_thread_num() allocate (keys_out(N_int,2,size_max), hole_save(N_int,2), & key(N_int,2),hole(N_int,2), particle(N_int,2), hole_tmp(N_int,2),& particle_tmp(N_int,2), occ_particle(N_int*bit_kind_size,2), & @@ -396,7 +402,8 @@ subroutine $subroutine($params_main) integer :: iproc $initialization - PROVIDE H_apply_buffer_allocated mo_bielec_integrals_in_map + PROVIDE H_apply_buffer_allocated mo_bielec_integrals_in_map psi_det_generators psi_coef_generators + nmax = mod( N_det_generators,nproc ) @@ -406,6 +413,7 @@ subroutine $subroutine($params_main) call wall_time(wall_0) + iproc = 0 allocate( mask(N_int,2,6) ) do i_generator=1,nmax @@ -443,12 +451,12 @@ subroutine $subroutine($params_main) 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) + i_generator, iproc $params_post) endif if($do_mono_excitations)then 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) + i_generator, iproc $params_post) endif call wall_time(wall_1) $printout_always @@ -463,7 +471,6 @@ subroutine $subroutine($params_main) !$OMP PARALLEL DEFAULT(SHARED) & !$OMP PRIVATE(i_generator,wall_1,wall_0,ispin,k,mask,iproc) call wall_time(wall_0) - iproc = 0 !$ iproc = omp_get_thread_num() allocate( mask(N_int,2,6) ) !$OMP DO SCHEDULE(dynamic,1) diff --git a/src/Full_CI/EZFIO.cfg b/src/Full_CI/EZFIO.cfg index febd0530..f660a386 100644 --- a/src/Full_CI/EZFIO.cfg +++ b/src/Full_CI/EZFIO.cfg @@ -1,42 +1,42 @@ [N_det_max_fci] type: Det_number_max doc: Max number of determinants in the wave function -interface: output +interface: input default: 10000 [N_det_max_fci_property] type: Det_number_max doc: Max number of determinants in the wave function when you select for a given property -interface: output +interface: input default: 10000 [do_pt2_end] type: logical doc: If true, compute the PT2 at the end of the selection -interface: output +interface: input default: true [PT2_max] type: PT2_energy doc: The selection process stops when the largest PT2 (for all the state is lower than pt2_max in absolute value -interface: output +interface: input default: 0.0001 [var_pt2_ratio] type: Normalized_float doc: The selection process stops when the energy ratio variational/(variational+PT2) is equal to var_pt2_ratio -interface: output +interface: input default: 0.75 [energy] type: double precision doc: "Calculated Full CI energy" -interface: input +interface: output [energy_pt2] type: double precision doc: "Calculated Full CI energy" -interface: input +interface: output diff --git a/src/Generators_CAS/generators.irp.f b/src/Generators_CAS/generators.irp.f index 5dcfc1b0..511877a0 100644 --- a/src/Generators_CAS/generators.irp.f +++ b/src/Generators_CAS/generators.irp.f @@ -62,7 +62,6 @@ END_PROVIDER 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/MRCC/NEEDED_MODULES b/src/MRCC/NEEDED_MODULES index 3711786f..31ed6f4e 100644 --- a/src/MRCC/NEEDED_MODULES +++ b/src/MRCC/NEEDED_MODULES @@ -1,2 +1,2 @@ -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 +AOs BiInts Bitmask CAS_SD Dets Electrons Ezfio_files Generators_CAS Generators_full Hartree_Fock MOGuess MonoInts MOs Nuclei Output Perturbation Properties Selectors_full Utils diff --git a/src/NEEDED_MODULES b/src/NEEDED_MODULES index 7c84ef42..8f2f7a20 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 MRCC +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 DDCI_selected MRCC diff --git a/src/Utils/LinearAlgebra.irp.f b/src/Utils/LinearAlgebra.irp.f index 7d173246..a7462f94 100644 --- a/src/Utils/LinearAlgebra.irp.f +++ b/src/Utils/LinearAlgebra.irp.f @@ -451,3 +451,6 @@ subroutine set_zero_extra_diag(i1,i2,matrix,lda,m) end + + + From 99e63935d4ca6decd515c420a317fc0d2548d9a9 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 9 Apr 2015 21:46:37 +0200 Subject: [PATCH 05/12] MRCC works --- src/CAS_SD/options.irp.f | 32 ------- src/MRCC/README.rst | 88 +++++++++++++++-- src/MRCC/mrcc.irp.f | 70 ++++++++++---- src/MRCC/mrcc_dress.irp.f | 193 +++++++++++++++++++++++++------------- src/MRCC/mrcc_utils.irp.f | 46 +++++---- 5 files changed, 284 insertions(+), 145 deletions(-) delete mode 100644 src/CAS_SD/options.irp.f diff --git a/src/CAS_SD/options.irp.f b/src/CAS_SD/options.irp.f deleted file mode 100644 index 553fbe9f..00000000 --- a/src/CAS_SD/options.irp.f +++ /dev/null @@ -1,32 +0,0 @@ -BEGIN_SHELL [ /usr/bin/python ] -from ezfio_with_default import EZFIO_Provider -T = EZFIO_Provider() -T.set_type ( "integer" ) -T.set_name ( "N_det_max_fci" ) -T.set_doc ( "Max number of determinants in the wave function" ) -T.set_ezfio_dir ( "full_ci" ) -T.set_ezfio_name( "N_det_max_fci" ) -T.set_output ( "output_full_ci" ) -print T - -T.set_type ( "logical" ) -T.set_name ( "do_pt2_end" ) -T.set_doc ( "If true, compute the PT2 at the end of the selection" ) -T.set_ezfio_name( "do_pt2_end" ) -print T - -T.set_type ( "double precision" ) -T.set_name ( "pt2_max" ) -T.set_doc ( """The selection process stops when the largest PT2 (for all the states) -is lower than pt2_max in absolute value""" ) -T.set_ezfio_name( "pt2_max" ) -print T - -T.set_type ( "double precision" ) -T.set_name ( "var_pt2_ratio" ) -T.set_doc ( """The selection process stops when the energy ratio variational/(variational+PT2) -is equal to var_pt2_ratio""" ) -T.set_ezfio_name( "var_pt2_ratio" ) -print T -END_SHELL - diff --git a/src/MRCC/README.rst b/src/MRCC/README.rst index 1fbcccdb..3720884a 100644 --- a/src/MRCC/README.rst +++ b/src/MRCC/README.rst @@ -1,14 +1,82 @@ --4.142795384334731 +E(CAS+SD) = -3.93510180989509 + + -0.726935163332 + ++++---- + ++++---- + + 0.685482292841 + +++-+--- + +++-+--- + + 0.0409791961003 + ++-+-+-- + ++-+-+-- + +c2*c3/c1 = -0.038642391672 + + + -0.996413146877 + ++++---- + ++++---- + + -0.0598366139154 + ++-+-+-- + +++-+--- + + -0.0598366139154 + +++-+--- + ++-+-+-- + + +======= + + 0.706616635698 + +++-+--- + +++-+--- + + -0.692594178414 + ++++---- + ++++---- + + -0.0915716740265 + ++-+-+-- + +++-+--- + + -0.0915716740265 + +++-+--- + ++-+-+-- + + 0.0461418323107 + ++-+-+-- + ++-+-+-- + + -0.0458957789452 + ++--++-- + ++--++-- + + +---- + + 0.713102867186 + ++++---- + ++++---- + + -0.688067688244 + +++-+--- + +++-+--- + + 0.089262694488 + +++-+--- + ++-+-+-- + + 0.089262694488 + ++-+-+-- + +++-+--- + + -0.0459510603899 + ++-+-+-- + ++-+-+-- -4.695183071437694E-002 - Determinant 64 - --------------------------------------------- - 000000000000002E|000000000000002E - |-+++-+----------------------------------------------------------| - |-+++-+----------------------------------------------------------| -CAS-SD: -4.14214374069306 - : -4.14230904320551 -E0 = -11.5634986758976 diff --git a/src/MRCC/mrcc.irp.f b/src/MRCC/mrcc.irp.f index 663b9f5a..48587b6d 100644 --- a/src/MRCC/mrcc.irp.f +++ b/src/MRCC/mrcc.irp.f @@ -3,6 +3,8 @@ program mrcc read_wf = .True. TOUCH read_wf call run + call run_mrcc +! call run_mrcc_test end subroutine run @@ -12,8 +14,7 @@ subroutine run print *, N_det_cas print *, N_det_sd -! call update_generators - integer :: i + integer :: i,j print *, 'CAS' print *, '===' do i=1,N_det_cas @@ -21,23 +22,54 @@ subroutine run 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 -! + 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 *, 'xxx', 'Energy CAS+SD', ci_energy +end +subroutine run_mrcc_test + implicit none + integer :: i,j + double precision :: pt2 + pt2 = 0.d0 + do j=1,N_det + do i=1,N_det + pt2 += psi_coef(i,1)*psi_coef(j,1) * delta_ij(i,j,1) + enddo + enddo + print *, ci_energy(1) + print *, ci_energy(1)+pt2 +end +subroutine run_mrcc + implicit none + integer :: i,j print *, 'MRCC' print *, '====' - 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 + print *, '' + print *, 'CAS+SD energy : ', ci_energy_dressed(:) + print *, '' + +! call diagonalize_ci_dressed +! call save_wavefunction_unsorted + double precision :: E_new, E_old, delta_e + integer :: iteration + E_new = 0.d0 + delta_E = 1.d0 + iteration = 0 + do while (delta_E > 1.d-8) + iteration += 1 + print *, '===========================' + print *, 'MRCC Iteration', iteration + print *, '===========================' + print *, '' + E_old = sum(ci_energy_dressed) + call diagonalize_ci_dressed + E_new = sum(ci_energy_dressed) + delta_E = dabs(E_new - E_old) + call write_double(6,ci_energy_dressed(1),"MRCC energy") + enddo + end diff --git a/src/MRCC/mrcc_dress.irp.f b/src/MRCC/mrcc_dress.irp.f index d9c8978c..2e913b02 100644 --- a/src/MRCC/mrcc_dress.irp.f +++ b/src/MRCC/mrcc_dress.irp.f @@ -1,3 +1,120 @@ +subroutine mrcc_dress(delta_ij_,Ndet,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 + double precision, intent(inout) :: delta_ij_(Ndet,Ndet,*) + + integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected) + integer :: i,j,k,l,m + logical :: is_in_wavefunction + integer :: degree_alpha(psi_det_size) + integer :: degree_I(psi_det_size) + integer :: idx_I(0:psi_det_size) + integer :: idx_alpha(0:psi_det_size) + logical :: good + + integer(bit_kind) :: tq(Nint,2,n_selected) + integer :: N_tq, c_ref ,degree + integer :: connected_to_ref + + call find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq) + + double precision :: hIk, hIl, hla, dIk(N_states), dka(N_states), dIa(N_states) + double precision :: haj, phase, phase2 + double precision :: f(N_states), ci_inv(N_states) + integer :: exc(0:2,2,2) + integer :: h1,h2,p1,p2,s1,s2 + integer(bit_kind):: tmp_det(Nint,2) + integer :: iint, ipos +! integer :: istate, i_sd, i_cas + + + + ! |I> + + ! |alpha> + do i=1,N_tq + call get_excitation_degree_vector(psi_sd,tq(1,1,i),degree_alpha,Nint,N_det_sd,idx_alpha) + + ! |I> + do j=1,N_det_cas + ! Find triples and quadruple grand parents + call get_excitation_degree(tq(1,1,i),psi_cas(1,1,j),degree,Nint) + if (degree > 4) then + cycle + endif + dIa(:) = 0.d0 + ! |alpha> + do k=1,idx_alpha(0) + call get_excitation_degree(psi_cas(1,1,j),psi_sd(1,1,idx_alpha(k)),degree,Nint) + if (degree > 2) then + cycle + endif + ! + ! + call i_h_j(psi_cas(1,1,j),psi_sd(1,1,idx_alpha(k)),Nint,hIk) + dIk(:) = hIk * lambda_mrcc(idx_alpha(k),:) + ! Exc(k -> alpha) + call get_excitation(psi_sd(1,1,idx_alpha(k)),tq(1,1,i),exc,degree,phase,Nint) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + tmp_det(1:Nint,1:2) = psi_cas(1,1,j) + ! Hole (see list_to_bitstring) + iint = ishft(h1-1,-bit_kind_shift) + 1 + ipos = h1-ishft((iint-1),bit_kind_shift)-1 + tmp_det(iint,s1) = ibclr(tmp_det(iint,s1),ipos) + + ! Particle + iint = ishft(p1-1,-bit_kind_shift) + 1 + ipos = p1-ishft((iint-1),bit_kind_shift)-1 + tmp_det(iint,s1) = ibset(tmp_det(iint,s1),ipos) + if (degree == 2) then + ! Hole (see list_to_bitstring) + iint = ishft(h2-1,-bit_kind_shift) + 1 + ipos = h2-ishft((iint-1),bit_kind_shift)-1 + tmp_det(iint,s2) = ibclr(tmp_det(iint,s2),ipos) + + ! Particle + iint = ishft(p2-1,-bit_kind_shift) + 1 + ipos = p2-ishft((iint-1),bit_kind_shift)-1 + tmp_det(iint,s2) = ibset(tmp_det(iint,s2),ipos) + endif + + dka(:) = 0.d0 + do l=k+1,idx_alpha(0) + call get_excitation_degree(tmp_det,psi_sd(1,1,idx_alpha(l)),degree,Nint) + if (degree == 0) then + call get_excitation(psi_cas(1,1,j),psi_sd(1,1,idx_alpha(l)),exc,degree,phase2,Nint) + call i_h_j(psi_cas(1,1,j),psi_sd(1,1,idx_alpha(l)),Nint,hIl) + dka(:) = hIl * lambda_mrcc(idx_alpha(l),:) * phase * phase2 + exit + endif + enddo + do l=1,N_states + dIa(l) += dka(l)*dIk(l) + enddo + enddo + ci_inv(1:N_states) = 1.d0/psi_cas_coefs(j,1:N_states) + do l=1,idx_alpha(0) + k = idx_alpha(l) + call i_h_j(tq(1,1,i),psi_sd(1,1,idx_alpha(l)),Nint,hla) + do m=1,N_states + delta_ij_(idx_sd(k),idx_cas(j),m) += dIa(m) * hla + delta_ij_(idx_cas(j),idx_sd(k),m) += dIa(m) * hla + delta_ij_(idx_cas(j),idx_cas(j),m) -= dIa(m) * hla * ci_inv(m) * psi_sd_coefs(k,m) + enddo + enddo + enddo + enddo +end + + + + + + + subroutine mrcc_dress_simple(delta_ij_sd_,Ndet_sd,i_generator,n_selected,det_buffer,Nint,iproc) use bitmasks implicit none @@ -18,35 +135,7 @@ subroutine mrcc_dress_simple(delta_ij_sd_,Ndet_sd,i_generator,n_selected,det_buf 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_det_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 + call find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq) ! Compute / (E0 - Haa) double precision :: hka, haa @@ -73,30 +162,22 @@ subroutine mrcc_dress_simple(delta_ij_sd_,Ndet_sd,i_generator,n_selected,det_buf end - - - - - - -subroutine mrcc_dress(delta_ij_sd_,Ndet_sd,i_generator,n_selected,det_buffer,Nint,iproc) +subroutine find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq) 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, intent(in) :: i_generator,n_selected, Nint integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected) integer :: i,j,k,m - integer :: new_size logical :: is_in_wavefunction integer :: degree(psi_det_size) integer :: idx(0:psi_det_size) logical :: good - integer(bit_kind) :: tq(Nint,2,n_selected) - integer :: N_tq, c_ref + integer(bit_kind), intent(out) :: tq(Nint,2,n_selected) + integer, intent(out) :: N_tq + integer :: c_ref integer :: connected_to_ref N_tq = 0 @@ -129,27 +210,11 @@ subroutine mrcc_dress(delta_ij_sd_,Ndet_sd,i_generator,n_selected,det_buffer,Nin endif enddo - ! Compute / (E0 - Haa) - double precision :: hka, haa - double precision :: haj - double precision :: f(N_states) - - 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) - 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 - enddo - enddo - enddo end + + + + + + diff --git a/src/MRCC/mrcc_utils.irp.f b/src/MRCC/mrcc_utils.irp.f index 01e566ae..a2b91645 100644 --- a/src/MRCC/mrcc_utils.irp.f +++ b/src/MRCC/mrcc_utils.irp.f @@ -164,14 +164,7 @@ BEGIN_PROVIDER [ double precision, delta_ij_sd, (N_det_sd, N_det_sd,N_states) ] ! Dressing matrix in SD basis END_DOC delta_ij_sd = 0.d0 - if (dressing_type == "MRCC") then - call H_apply_mrcc(delta_ij_sd,N_det_sd) - else if (dressing_type == "Simple") then - call H_apply_mrcc_simple(delta_ij_sd,N_det_sd) - else - print *, irp_here - stop 'dressing' - endif + call H_apply_mrcc_simple(delta_ij_sd,N_det_sd) END_PROVIDER BEGIN_PROVIDER [ double precision, delta_ij, (N_det,N_det,N_states) ] @@ -181,13 +174,17 @@ BEGIN_PROVIDER [ double precision, delta_ij, (N_det,N_det,N_states) ] 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) + if (dressing_type == "MRCC") then + call H_apply_mrcc(delta_ij,N_det) + else if (dressing_type == "Simple") then + 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 - enddo - enddo + endif END_PROVIDER BEGIN_PROVIDER [ double precision, h_matrix_dressed, (N_det,N_det) ] @@ -199,9 +196,6 @@ BEGIN_PROVIDER [ double precision, h_matrix_dressed, (N_det,N_det) ] 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 @@ -273,10 +267,22 @@ BEGIN_PROVIDER [ double precision, CI_energy_dressed, (N_states_diag) ] 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 +subroutine diagonalize_CI_dressed + implicit none + BEGIN_DOC +! Replace the coefficients of the CI states by the coefficients of the +! eigenstates of the CI matrix + END_DOC + integer :: i,j + do j=1,N_states_diag + do i=1,N_det + psi_coef(i,j) = CI_eigenvectors_dressed(i,j) + enddo + enddo + SOFT_TOUCH psi_coef + +end From 0f2bdedb902c14c64defa6e5e50577332599b4f3 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 9 Apr 2015 22:32:25 +0200 Subject: [PATCH 06/12] Moves psi_cas in Dets --- src/MRCC/mrcc.irp.f | 27 +++---- src/MRCC/mrcc_dress.irp.f | 89 +++++++++++++---------- src/MRCC/mrcc_utils.irp.f | 145 +++----------------------------------- 3 files changed, 70 insertions(+), 191 deletions(-) diff --git a/src/MRCC/mrcc.irp.f b/src/MRCC/mrcc.irp.f index 48587b6d..71b87ddb 100644 --- a/src/MRCC/mrcc.irp.f +++ b/src/MRCC/mrcc.irp.f @@ -10,10 +10,6 @@ end subroutine run implicit none - print *, N_det - print *, N_det_cas - print *, N_det_sd - integer :: i,j print *, 'CAS' print *, '===' @@ -22,13 +18,13 @@ subroutine run 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 - print *, 'xxx', 'Energy CAS+SD', ci_energy +! print *, 'SD' +! print *, '==' +! do i=1,N_det_non_cas +! print *, psi_non_cas_coefs(i,:) +! call debug_det(psi_non_cas(1,1,i),N_int) +! enddo + call write_double(6,ci_energy(1),"Initial CI energy") end subroutine run_mrcc_test implicit none @@ -46,14 +42,7 @@ end subroutine run_mrcc implicit none integer :: i,j - print *, 'MRCC' - print *, '====' - print *, '' - print *, 'CAS+SD energy : ', ci_energy_dressed(:) - print *, '' -! call diagonalize_ci_dressed -! call save_wavefunction_unsorted double precision :: E_new, E_old, delta_e integer :: iteration E_new = 0.d0 @@ -71,5 +60,7 @@ subroutine run_mrcc delta_E = dabs(E_new - E_old) call write_double(6,ci_energy_dressed(1),"MRCC energy") enddo + call ezfio_set_mrcc_energy(ci_energy_dressed(1)) +! call save_wavefunction end diff --git a/src/MRCC/mrcc_dress.irp.f b/src/MRCC/mrcc_dress.irp.f index 2e913b02..96f278ca 100644 --- a/src/MRCC/mrcc_dress.irp.f +++ b/src/MRCC/mrcc_dress.irp.f @@ -28,38 +28,44 @@ subroutine mrcc_dress(delta_ij_,Ndet,i_generator,n_selected,det_buffer,Nint,ipro integer :: h1,h2,p1,p2,s1,s2 integer(bit_kind):: tmp_det(Nint,2) integer :: iint, ipos -! integer :: istate, i_sd, i_cas + integer :: i_state, k_sd, l_sd, i_I, i_alpha ! |I> ! |alpha> - do i=1,N_tq - call get_excitation_degree_vector(psi_sd,tq(1,1,i),degree_alpha,Nint,N_det_sd,idx_alpha) + do i_alpha=1,N_tq + call get_excitation_degree_vector(psi_non_cas,tq(1,1,i_alpha),degree_alpha,Nint,N_det_non_cas,idx_alpha) ! |I> - do j=1,N_det_cas + do i_I=1,N_det_cas ! Find triples and quadruple grand parents - call get_excitation_degree(tq(1,1,i),psi_cas(1,1,j),degree,Nint) + call get_excitation_degree(tq(1,1,i_alpha),psi_cas(1,1,i_I),degree,Nint) if (degree > 4) then cycle endif - dIa(:) = 0.d0 + + do i_state=1,N_states + dIa(i_state) = 0.d0 + enddo + ! |alpha> - do k=1,idx_alpha(0) - call get_excitation_degree(psi_cas(1,1,j),psi_sd(1,1,idx_alpha(k)),degree,Nint) + do k_sd=1,idx_alpha(0) + call get_excitation_degree(psi_cas(1,1,i_I),psi_non_cas(1,1,idx_alpha(k_sd)),degree,Nint) if (degree > 2) then cycle endif - ! + ! ! - call i_h_j(psi_cas(1,1,j),psi_sd(1,1,idx_alpha(k)),Nint,hIk) - dIk(:) = hIk * lambda_mrcc(idx_alpha(k),:) - ! Exc(k -> alpha) - call get_excitation(psi_sd(1,1,idx_alpha(k)),tq(1,1,i),exc,degree,phase,Nint) + call i_h_j(psi_cas(1,1,i_I),psi_non_cas(1,1,idx_alpha(k_sd)),Nint,hIk) + do i_state=1,N_states + dIk(i_state) = hIk * lambda_mrcc(idx_alpha(k_sd),i_state) + enddo + ! |l> = Exc(k -> alpha) |I> + call get_excitation(psi_non_cas(1,1,idx_alpha(k_sd)),tq(1,1,i_alpha),exc,degree,phase,Nint) call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) - tmp_det(1:Nint,1:2) = psi_cas(1,1,j) + tmp_det(1:Nint,1:2) = psi_cas(1,1,i_I) ! Hole (see list_to_bitstring) iint = ishft(h1-1,-bit_kind_shift) + 1 ipos = h1-ishft((iint-1),bit_kind_shift)-1 @@ -81,28 +87,37 @@ subroutine mrcc_dress(delta_ij_,Ndet,i_generator,n_selected,det_buffer,Nint,ipro tmp_det(iint,s2) = ibset(tmp_det(iint,s2),ipos) endif - dka(:) = 0.d0 - do l=k+1,idx_alpha(0) - call get_excitation_degree(tmp_det,psi_sd(1,1,idx_alpha(l)),degree,Nint) + ! + do i_state=1,N_states + dka(i_state) = 0.d0 + enddo + do l_sd=k_sd+1,idx_alpha(0) + call get_excitation_degree(tmp_det,psi_non_cas(1,1,idx_alpha(l_sd)),degree,Nint) if (degree == 0) then - call get_excitation(psi_cas(1,1,j),psi_sd(1,1,idx_alpha(l)),exc,degree,phase2,Nint) - call i_h_j(psi_cas(1,1,j),psi_sd(1,1,idx_alpha(l)),Nint,hIl) - dka(:) = hIl * lambda_mrcc(idx_alpha(l),:) * phase * phase2 + call get_excitation(psi_cas(1,1,i_I),psi_non_cas(1,1,idx_alpha(l_sd)),exc,degree,phase2,Nint) + call i_h_j(psi_cas(1,1,i_I),psi_non_cas(1,1,idx_alpha(l_sd)),Nint,hIl) + do i_state=1,N_states + dka(i_state) = hIl * lambda_mrcc(idx_alpha(l_sd),i_state) * phase * phase2 + enddo exit endif enddo - do l=1,N_states - dIa(l) += dka(l)*dIk(l) + do i_state=1,N_states + dIa(i_state) = dIa(i_state) + dIk(i_state) * dka(i_state) enddo enddo - ci_inv(1:N_states) = 1.d0/psi_cas_coefs(j,1:N_states) - do l=1,idx_alpha(0) - k = idx_alpha(l) - call i_h_j(tq(1,1,i),psi_sd(1,1,idx_alpha(l)),Nint,hla) + + do i_state=1,N_states + ci_inv(i_state) = 1.d0/psi_cas_coefs(i_I,i_state) + enddo + + do l_sd=1,idx_alpha(0) + k_sd = idx_alpha(l_sd) + call i_h_j(tq(1,1,i_alpha),psi_non_cas(1,1,idx_alpha(l_sd)),Nint,hla) do m=1,N_states - delta_ij_(idx_sd(k),idx_cas(j),m) += dIa(m) * hla - delta_ij_(idx_cas(j),idx_sd(k),m) += dIa(m) * hla - delta_ij_(idx_cas(j),idx_cas(j),m) -= dIa(m) * hla * ci_inv(m) * psi_sd_coefs(k,m) + delta_ij_(idx_non_cas(k_sd),idx_cas(i_I),m) += dIa(m) * hla + delta_ij_(idx_cas(i_I),idx_non_cas(k_sd),m) += dIa(m) * hla + delta_ij_(idx_cas(i_I),idx_cas(i_I),m) -= dIa(m) * hla * ci_inv(m) * psi_non_cas_coefs(k_sd,m) enddo enddo enddo @@ -115,13 +130,13 @@ end -subroutine mrcc_dress_simple(delta_ij_sd_,Ndet_sd,i_generator,n_selected,det_buffer,Nint,iproc) +subroutine mrcc_dress_simple(delta_ij_non_cas_,Ndet_non_cas,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, intent(in) :: Ndet_non_cas + double precision, intent(inout) :: delta_ij_non_cas_(Ndet_non_cas,Ndet_non_cas,*) integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected) integer :: i,j,k,m @@ -143,18 +158,18 @@ subroutine mrcc_dress_simple(delta_ij_sd_,Ndet_sd,i_generator,n_selected,det_buf double precision :: f(N_states) do i=1,N_tq - call get_excitation_degree_vector(psi_sd,tq(1,1,i),degree,Nint,Ndet_sd,idx) + call get_excitation_degree_vector(psi_non_cas,tq(1,1,i),degree,Nint,Ndet_non_cas,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) + call i_h_j(tq(1,1,i),psi_non_cas(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) + call i_h_j(tq(1,1,i),psi_non_cas(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) + delta_ij_non_cas_(idx(k), idx(j),m) += haj*hka* f(m) + delta_ij_non_cas_(idx(j), idx(k),m) += haj*hka* f(m) enddo enddo enddo diff --git a/src/MRCC/mrcc_utils.irp.f b/src/MRCC/mrcc_utils.irp.f index a2b91645..4a1fc1cc 100644 --- a/src/MRCC/mrcc_utils.irp.f +++ b/src/MRCC/mrcc_utils.irp.f @@ -1,130 +1,3 @@ - use bitmasks - -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_sd, (N_det) ] -&BEGIN_PROVIDER [ integer, N_det_sd] - implicit none - BEGIN_DOC - ! SD - END_DOC - integer :: i_sd,j,k - integer :: degree - logical :: in_cas - i_sd =0 - do k=1,N_det - in_cas = .False. - 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 - in_cas = .True. - 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 -END_PROVIDER BEGIN_PROVIDER [ double precision, lambda_mrcc, (psi_det_size,n_states) ] implicit none @@ -133,15 +6,15 @@ BEGIN_PROVIDER [ double precision, lambda_mrcc, (psi_det_size,n_states) ] 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, & + do i=1,N_det_non_cas + call i_h_psi(psi_non_cas(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) = psi_non_cas_coefs(i,k)/ihpsi(k) lambda_mrcc(i,k) = min( lambda_mrcc (i,k),0.d0 ) endif enddo @@ -158,13 +31,13 @@ BEGIN_PROVIDER [ character*(32), dressing_type ] dressing_type = "MRCC" END_PROVIDER -BEGIN_PROVIDER [ double precision, delta_ij_sd, (N_det_sd, N_det_sd,N_states) ] +BEGIN_PROVIDER [ double precision, delta_ij_non_cas, (N_det_non_cas, N_det_non_cas,N_states) ] implicit none BEGIN_DOC ! Dressing matrix in SD basis END_DOC - delta_ij_sd = 0.d0 - call H_apply_mrcc_simple(delta_ij_sd,N_det_sd) + delta_ij_non_cas = 0.d0 + call H_apply_mrcc_simple(delta_ij_non_cas,N_det_non_cas) END_PROVIDER BEGIN_PROVIDER [ double precision, delta_ij, (N_det,N_det,N_states) ] @@ -178,9 +51,9 @@ BEGIN_PROVIDER [ double precision, delta_ij, (N_det,N_det,N_states) ] call H_apply_mrcc(delta_ij,N_det) else if (dressing_type == "Simple") then 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) + do j=1,N_det_non_cas + do i=1,N_det_non_cas + delta_ij(idx_non_cas(i),idx_non_cas(j),m) = delta_ij_non_cas(i,j,m) enddo enddo enddo From b5b211aaea485c7688d48ae493c8a9e97025ee7a Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 9 Apr 2015 23:59:06 +0200 Subject: [PATCH 07/12] MRCC is parallelized --- src/Dets/determinants.irp.f | 37 ++++++++++++++----- src/MRCC/H_apply.irp.f | 2 +- src/MRCC/mrcc.irp.f | 6 +-- src/MRCC/mrcc_dress.irp.f | 74 ++++++++++++++++++++++++------------- src/MRCC/mrcc_utils.irp.f | 14 +++---- 5 files changed, 85 insertions(+), 48 deletions(-) diff --git a/src/Dets/determinants.irp.f b/src/Dets/determinants.irp.f index 0a0ec864..00e683fc 100644 --- a/src/Dets/determinants.irp.f +++ b/src/Dets/determinants.irp.f @@ -467,33 +467,52 @@ END_PROVIDER ! to accelerate the search of a random determinant in the wave ! function. END_DOC + + call sort_dets_by_det_search_key(N_det, psi_det, psi_coef, & + psi_det_sorted_bit, psi_coef_sorted_bit) +END_PROVIDER + +subroutine sort_dets_by_det_search_key(Ndet, det_in, coef_in, det_out, coef_out) + use bitmasks + implicit none + integer, intent(in) :: Ndet + integer(bit_kind), intent(in) :: det_in (N_int,2,psi_det_size) + double precision , intent(in) :: coef_in(psi_det_size,N_states) + integer(bit_kind), intent(out) :: det_out (N_int,2,psi_det_size) + double precision , intent(out) :: coef_out(psi_det_size,N_states) + BEGIN_DOC + ! Determinants are sorted are sorted according to their det_search_key. + ! Useful to accelerate the search of a random determinant in the wave + ! function. + END_DOC integer :: i,j,k integer, allocatable :: iorder(:) integer*8, allocatable :: bit_tmp(:) integer*8, external :: det_search_key - allocate ( iorder(N_det), bit_tmp(N_det) ) + allocate ( iorder(Ndet), bit_tmp(Ndet) ) - do i=1,N_det + do i=1,Ndet iorder(i) = i !$DIR FORCEINLINE - bit_tmp(i) = det_search_key(psi_det(1,1,i),N_int) + bit_tmp(i) = det_search_key(det_in(1,1,i),N_int) enddo - call i8sort(bit_tmp,iorder,N_det) + call i8sort(bit_tmp,iorder,Ndet) !DIR$ IVDEP - do i=1,N_det + do i=1,Ndet do j=1,N_int - psi_det_sorted_bit(j,1,i) = psi_det(j,1,iorder(i)) - psi_det_sorted_bit(j,2,i) = psi_det(j,2,iorder(i)) + det_out(j,1,i) = det_in(j,1,iorder(i)) + det_out(j,2,i) = det_in(j,2,iorder(i)) enddo do k=1,N_states - psi_coef_sorted_bit(i,k) = psi_coef(iorder(i),k) + coef_out(i,k) = coef_in(iorder(i),k) enddo enddo deallocate(iorder, bit_tmp) -END_PROVIDER +end + subroutine int_of_3_highest_electrons( det_in, res, Nint ) implicit none diff --git a/src/MRCC/H_apply.irp.f b/src/MRCC/H_apply.irp.f index a8519c25..cadaca32 100644 --- a/src/MRCC/H_apply.irp.f +++ b/src/MRCC/H_apply.irp.f @@ -18,7 +18,7 @@ s.data["decls_main"] += """ s.data["finalization"] = "" s.data["copy_buffer"] = "" s.data["generate_psi_guess"] = "" -s.data["size_max"] = "256" +s.data["size_max"] = "3072" print s diff --git a/src/MRCC/mrcc.irp.f b/src/MRCC/mrcc.irp.f index 71b87ddb..0ca9f169 100644 --- a/src/MRCC/mrcc.irp.f +++ b/src/MRCC/mrcc.irp.f @@ -14,14 +14,14 @@ subroutine run print *, 'CAS' print *, '===' do i=1,N_det_cas - print *, psi_cas_coefs(i,:) + print *, psi_cas_coef(i,:) call debug_det(psi_cas(1,1,i),N_int) enddo ! print *, 'SD' ! print *, '==' ! do i=1,N_det_non_cas -! print *, psi_non_cas_coefs(i,:) +! print *, psi_non_cas_coef(i,:) ! call debug_det(psi_non_cas(1,1,i),N_int) ! enddo call write_double(6,ci_energy(1),"Initial CI energy") @@ -56,9 +56,9 @@ subroutine run_mrcc print *, '' E_old = sum(ci_energy_dressed) call diagonalize_ci_dressed + call write_double(6,ci_energy_dressed(1),"MRCC energy") E_new = sum(ci_energy_dressed) delta_E = dabs(E_new - E_old) - call write_double(6,ci_energy_dressed(1),"MRCC energy") enddo call ezfio_set_mrcc_energy(ci_energy_dressed(1)) ! call save_wavefunction diff --git a/src/MRCC/mrcc_dress.irp.f b/src/MRCC/mrcc_dress.irp.f index 96f278ca..083cea5a 100644 --- a/src/MRCC/mrcc_dress.irp.f +++ b/src/MRCC/mrcc_dress.irp.f @@ -1,3 +1,17 @@ +use omp_lib + +BEGIN_PROVIDER [ integer(omp_lock_kind), psi_cas_lock, (psi_det_size) ] + implicit none + BEGIN_DOC + ! Locks on CAS determinants to fill delta_ij + END_DOC + integer :: i + do i=1,psi_det_size + call omp_init_lock( psi_cas_lock(i) ) + enddo + +END_PROVIDER + subroutine mrcc_dress(delta_ij_,Ndet,i_generator,n_selected,det_buffer,Nint,iproc) use bitmasks implicit none @@ -7,11 +21,8 @@ subroutine mrcc_dress(delta_ij_,Ndet,i_generator,n_selected,det_buffer,Nint,ipro double precision, intent(inout) :: delta_ij_(Ndet,Ndet,*) integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected) - integer :: i,j,k,l,m - logical :: is_in_wavefunction + integer :: i,j,k,l integer :: degree_alpha(psi_det_size) - integer :: degree_I(psi_det_size) - integer :: idx_I(0:psi_det_size) integer :: idx_alpha(0:psi_det_size) logical :: good @@ -19,18 +30,19 @@ subroutine mrcc_dress(delta_ij_,Ndet,i_generator,n_selected,det_buffer,Nint,ipro integer :: N_tq, c_ref ,degree integer :: connected_to_ref + double precision :: hIk, hla, hIl, dIk(N_states), dka(N_states), dIa(N_states) + double precision, allocatable :: dIa_hla(:,:) + double precision :: haj, phase, phase2 + double precision :: f(N_states), ci_inv(N_states) + integer :: exc(0:2,2,2) + integer :: h1,h2,p1,p2,s1,s2 + integer(bit_kind) :: tmp_det(Nint,2) + integer :: iint, ipos + integer :: i_state, k_sd, l_sd, i_I, i_alpha + call find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq) - double precision :: hIk, hIl, hla, dIk(N_states), dka(N_states), dIa(N_states) - double precision :: haj, phase, phase2 - double precision :: f(N_states), ci_inv(N_states) - integer :: exc(0:2,2,2) - integer :: h1,h2,p1,p2,s1,s2 - integer(bit_kind):: tmp_det(Nint,2) - integer :: iint, ipos - integer :: i_state, k_sd, l_sd, i_I, i_alpha - - + allocate (dIa_hla(N_states,Ndet)) ! |I> @@ -60,12 +72,15 @@ subroutine mrcc_dress(delta_ij_,Ndet,i_generator,n_selected,det_buffer,Nint,ipro ! call i_h_j(psi_cas(1,1,i_I),psi_non_cas(1,1,idx_alpha(k_sd)),Nint,hIk) do i_state=1,N_states - dIk(i_state) = hIk * lambda_mrcc(idx_alpha(k_sd),i_state) + dIk(i_state) = hIk * lambda_mrcc(i_state,idx_alpha(k_sd)) enddo ! |l> = Exc(k -> alpha) |I> call get_excitation(psi_non_cas(1,1,idx_alpha(k_sd)),tq(1,1,i_alpha),exc,degree,phase,Nint) call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) - tmp_det(1:Nint,1:2) = psi_cas(1,1,i_I) + do k=1,N_int + tmp_det(k,1) = psi_cas(k,1,i_I) + tmp_det(k,2) = psi_cas(k,2,i_I) + enddo ! Hole (see list_to_bitstring) iint = ishft(h1-1,-bit_kind_shift) + 1 ipos = h1-ishft((iint-1),bit_kind_shift)-1 @@ -75,7 +90,7 @@ subroutine mrcc_dress(delta_ij_,Ndet,i_generator,n_selected,det_buffer,Nint,ipro iint = ishft(p1-1,-bit_kind_shift) + 1 ipos = p1-ishft((iint-1),bit_kind_shift)-1 tmp_det(iint,s1) = ibset(tmp_det(iint,s1),ipos) - if (degree == 2) then + if (degree_alpha(k_sd) == 2) then ! Hole (see list_to_bitstring) iint = ishft(h2-1,-bit_kind_shift) + 1 ipos = h2-ishft((iint-1),bit_kind_shift)-1 @@ -97,31 +112,39 @@ subroutine mrcc_dress(delta_ij_,Ndet,i_generator,n_selected,det_buffer,Nint,ipro call get_excitation(psi_cas(1,1,i_I),psi_non_cas(1,1,idx_alpha(l_sd)),exc,degree,phase2,Nint) call i_h_j(psi_cas(1,1,i_I),psi_non_cas(1,1,idx_alpha(l_sd)),Nint,hIl) do i_state=1,N_states - dka(i_state) = hIl * lambda_mrcc(idx_alpha(l_sd),i_state) * phase * phase2 + dka(i_state) = hIl * lambda_mrcc(i_state,idx_alpha(l_sd)) * phase * phase2 enddo exit endif enddo do i_state=1,N_states - dIa(i_state) = dIa(i_state) + dIk(i_state) * dka(i_state) + dIa(i_state) = dIa(i_state) + dIk(i_state) * dka(i_state) enddo enddo do i_state=1,N_states - ci_inv(i_state) = 1.d0/psi_cas_coefs(i_I,i_state) + ci_inv(i_state) = 1.d0/psi_cas_coef(i_I,i_state) enddo - do l_sd=1,idx_alpha(0) k_sd = idx_alpha(l_sd) call i_h_j(tq(1,1,i_alpha),psi_non_cas(1,1,idx_alpha(l_sd)),Nint,hla) - do m=1,N_states - delta_ij_(idx_non_cas(k_sd),idx_cas(i_I),m) += dIa(m) * hla - delta_ij_(idx_cas(i_I),idx_non_cas(k_sd),m) += dIa(m) * hla - delta_ij_(idx_cas(i_I),idx_cas(i_I),m) -= dIa(m) * hla * ci_inv(m) * psi_non_cas_coefs(k_sd,m) + do i_state=1,N_states + dIa_hla(i_state,k_sd) = dIa(i_state) * hla enddo enddo + call omp_set_lock( psi_cas_lock(i_I) ) + do l_sd=1,idx_alpha(0) + k_sd = idx_alpha(l_sd) + do i_state=1,N_states + delta_ij_(idx_non_cas(k_sd),idx_cas(i_I),i_state) += dIa_hla(i_state,k_sd) + delta_ij_(idx_cas(i_I),idx_non_cas(k_sd),i_state) += dIa_hla(i_state,k_sd) + delta_ij_(idx_cas(i_I),idx_cas(i_I),i_state) -= dIa_hla(i_state,k_sd) * ci_inv(i_state) * psi_non_cas_coef(k_sd,i_state) + enddo + enddo + call omp_unset_lock( psi_cas_lock(i_I) ) enddo enddo + deallocate (dIa_hla) end @@ -141,7 +164,6 @@ subroutine mrcc_dress_simple(delta_ij_non_cas_,Ndet_non_cas,i_generator,n_select integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected) integer :: i,j,k,m integer :: new_size - logical :: is_in_wavefunction integer :: degree(psi_det_size) integer :: idx(0:psi_det_size) logical :: good diff --git a/src/MRCC/mrcc_utils.irp.f b/src/MRCC/mrcc_utils.irp.f index 4a1fc1cc..eca1a706 100644 --- a/src/MRCC/mrcc_utils.irp.f +++ b/src/MRCC/mrcc_utils.irp.f @@ -1,5 +1,5 @@ -BEGIN_PROVIDER [ double precision, lambda_mrcc, (psi_det_size,n_states) ] +BEGIN_PROVIDER [ double precision, lambda_mrcc, (N_states,psi_det_size) ] implicit none BEGIN_DOC ! cm/ @@ -7,16 +7,12 @@ BEGIN_PROVIDER [ double precision, lambda_mrcc, (psi_det_size,n_states) ] integer :: i,k double precision :: ihpsi(N_states) do i=1,N_det_non_cas - call i_h_psi(psi_non_cas(1,1,i), psi_cas, psi_cas_coefs, N_int, N_det_cas, & - size(psi_cas_coefs,1), n_states, ihpsi) + call i_h_psi(psi_non_cas(1,1,i), psi_cas, psi_cas_coef, N_int, N_det_cas, & + size(psi_cas_coef,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_non_cas_coefs(i,k)/ihpsi(k) - lambda_mrcc(i,k) = min( lambda_mrcc (i,k),0.d0 ) - endif + lambda_mrcc(k,i) = psi_non_cas_coef(i,k)/ihpsi(k) + lambda_mrcc(k,i) = min( lambda_mrcc (k,i),0.d0 ) enddo enddo END_PROVIDER From 8bb2c8c8981d452946deb91b83dd66d838b1e89f Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 10 Apr 2015 00:44:12 +0200 Subject: [PATCH 08/12] Removed debug in connected_to_ref.irp.f --- src/Dets/connected_to_ref.irp.f | 54 ++++++++++++++++----------------- src/MRCC/mrcc.irp.f | 3 +- 2 files changed, 29 insertions(+), 28 deletions(-) diff --git a/src/Dets/connected_to_ref.irp.f b/src/Dets/connected_to_ref.irp.f index 07331782..3c7eb581 100644 --- a/src/Dets/connected_to_ref.irp.f +++ b/src/Dets/connected_to_ref.irp.f @@ -95,9 +95,9 @@ integer function get_index_in_psi_det_sorted_bit(key,Nint) 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) @@ -116,39 +116,39 @@ integer function get_index_in_psi_det_sorted_bit(key,Nint) enddo if (is_in_wavefunction) then get_index_in_psi_det_sorted_bit = i - exit -! return +! exit + return endif endif i += 1 if (i > N_det) then - exit -! return +! exit + return endif enddo ! 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/mrcc.irp.f b/src/MRCC/mrcc.irp.f index 0ca9f169..2b879882 100644 --- a/src/MRCC/mrcc.irp.f +++ b/src/MRCC/mrcc.irp.f @@ -55,11 +55,12 @@ subroutine run_mrcc print *, '===========================' print *, '' E_old = sum(ci_energy_dressed) - call diagonalize_ci_dressed call write_double(6,ci_energy_dressed(1),"MRCC energy") + call diagonalize_ci_dressed E_new = sum(ci_energy_dressed) delta_E = dabs(E_new - E_old) enddo + call write_double(6,ci_energy_dressed(1),"Final MRCC energy") call ezfio_set_mrcc_energy(ci_energy_dressed(1)) ! call save_wavefunction From 17126c779e955244b319875eb742f126d6c99e94 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 10 Apr 2015 00:45:09 +0200 Subject: [PATCH 09/12] Forgt EZFIO.cfg file --- src/CAS_SD/EZFIO.cfg | 36 ++++++++++++++++++++++++++++++++++++ 1 file changed, 36 insertions(+) create mode 100644 src/CAS_SD/EZFIO.cfg diff --git a/src/CAS_SD/EZFIO.cfg b/src/CAS_SD/EZFIO.cfg new file mode 100644 index 00000000..5629e90b --- /dev/null +++ b/src/CAS_SD/EZFIO.cfg @@ -0,0 +1,36 @@ +[N_det_max_cas_sd] +type: Det_number_max +doc: Max number of determinants in the wave function +interface: input +default: 10000 + +[do_pt2_end] +type: logical +doc: If true, compute the PT2 at the end of the selection +interface: input +default: true + +[PT2_max] +type: PT2_energy +doc: The selection process stops when the largest PT2 (for all the state is lower + than pt2_max in absolute value +interface: input +default: 0.0001 + +[var_pt2_ratio] +type: Normalized_float +doc: The selection process stops when the energy ratio variational/(variational+PT2) + is equal to var_pt2_ratio +interface: input +default: 0.75 + +[energy] +type: double precision +doc: "Calculated CAS-SD energy" +interface: output + +[energy_pt2] +type: double precision +doc: "Calculated selected CAS-SD energy with PT2 correction" +interface: output + From 44e4cbe4b3a0cec8379909bacada9cfb966dab33 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 10 Apr 2015 10:37:07 +0200 Subject: [PATCH 10/12] Forgot file --- src/Dets/psi_cas.irp.f | 114 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 114 insertions(+) create mode 100644 src/Dets/psi_cas.irp.f diff --git a/src/Dets/psi_cas.irp.f b/src/Dets/psi_cas.irp.f new file mode 100644 index 00000000..299e5e8f --- /dev/null +++ b/src/Dets/psi_cas.irp.f @@ -0,0 +1,114 @@ +use bitmasks + + BEGIN_PROVIDER [ integer(bit_kind), psi_cas, (N_int,2,psi_det_size) ] +&BEGIN_PROVIDER [ double precision, psi_cas_coef, (psi_det_size,n_states) ] +&BEGIN_PROVIDER [ integer, idx_cas, (psi_det_size) ] +&BEGIN_PROVIDER [ integer, N_det_cas ] + implicit none + BEGIN_DOC + ! CAS wave function, defined from the application of the CAS bitmask on the + ! determinants. idx_cas gives the indice of the CAS determinant in psi_det. + END_DOC + integer :: i, k, l + logical :: good + 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 = N_det_cas+1 + do k=1,N_int + psi_cas(k,1,N_det_cas) = psi_det(k,1,i) + psi_cas(k,2,N_det_cas) = psi_det(k,2,i) + enddo + idx_cas(N_det_cas) = i + do k=1,N_states + psi_cas_coef(N_det_cas,k) = psi_coef(i,k) + enddo + endif + enddo + call write_int(output_dets,N_det_cas, 'Number of determinants in the CAS') + +END_PROVIDER + + + BEGIN_PROVIDER [ integer(bit_kind), psi_cas_sorted_bit, (N_int,2,psi_det_size) ] +&BEGIN_PROVIDER [ double precision, psi_cas_coef_sorted_bit, (psi_det_size,N_states) ] + implicit none + BEGIN_DOC + ! CAS determinants sorted to accelerate the search of a random determinant in the wave + ! function. + END_DOC + call sort_dets_by_det_search_key(N_det_cas, psi_cas, psi_cas_coef, & + psi_cas_sorted_bit, psi_cas_coef_sorted_bit) + +END_PROVIDER + + + + BEGIN_PROVIDER [ integer(bit_kind), psi_non_cas, (N_int,2,psi_det_size) ] +&BEGIN_PROVIDER [ double precision, psi_non_cas_coef, (psi_det_size,n_states) ] +&BEGIN_PROVIDER [ integer, idx_non_cas, (psi_det_size) ] +&BEGIN_PROVIDER [ integer, N_det_non_cas ] + implicit none + BEGIN_DOC + ! Set of determinants which are not part of the CAS, defined from the application + ! of the CAS bitmask on the determinants. + ! idx_non_cas gives the indice of the determinant in psi_det. + END_DOC + integer :: i_non_cas,j,k + integer :: degree + logical :: in_cas + i_non_cas =0 + do k=1,N_det + in_cas = .False. + 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 + in_cas = .True. + exit + endif + enddo + if (.not.in_cas) then + double precision :: hij + i_non_cas += 1 + do j=1,N_int + psi_non_cas(j,1,i_non_cas) = psi_det(j,1,k) + psi_non_cas(j,2,i_non_cas) = psi_det(j,2,k) + enddo + do j=1,N_states + psi_non_cas_coef(i_non_cas,j) = psi_coef(k,j) + enddo + idx_non_cas(i_non_cas) = k + endif + enddo + N_det_non_cas = i_non_cas +END_PROVIDER + + BEGIN_PROVIDER [ integer(bit_kind), psi_non_cas_sorted_bit, (N_int,2,psi_det_size) ] +&BEGIN_PROVIDER [ double precision, psi_non_cas_coef_sorted_bit, (psi_det_size,N_states) ] + implicit none + BEGIN_DOC + ! CAS determinants sorted to accelerate the search of a random determinant in the wave + ! function. + END_DOC + call sort_dets_by_det_search_key(N_det_cas, psi_non_cas, psi_non_cas_coef, & + psi_non_cas_sorted_bit, psi_non_cas_coef_sorted_bit) + +END_PROVIDER + + + + + From 6ac88cf7891a0f0a2f26de80d2cf359cb336f439 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 10 Apr 2015 11:06:33 +0200 Subject: [PATCH 11/12] Forgot file --- src/MRCC/EZFIO.cfg | 5 +++++ 1 file changed, 5 insertions(+) create mode 100644 src/MRCC/EZFIO.cfg diff --git a/src/MRCC/EZFIO.cfg b/src/MRCC/EZFIO.cfg new file mode 100644 index 00000000..d596b785 --- /dev/null +++ b/src/MRCC/EZFIO.cfg @@ -0,0 +1,5 @@ +[energy] +type: double precision +doc: "Calculated MRCC energy" +interface: output + From 47d913a9cee930c0bc125bdae3bb9f653b983690 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 10 Apr 2015 17:26:07 +0200 Subject: [PATCH 12/12] Fixed CAS-SD --- src/CAS_SD/cas_sd_selected.irp.f | 15 ++++++++------- src/MRCC/NEEDED_MODULES | 2 +- src/MRCC/mrcc_utils.irp.f | 8 ++++++-- 3 files changed, 15 insertions(+), 10 deletions(-) diff --git a/src/CAS_SD/cas_sd_selected.irp.f b/src/CAS_SD/cas_sd_selected.irp.f index fc9d4dd2..86cac969 100644 --- a/src/CAS_SD/cas_sd_selected.irp.f +++ b/src/CAS_SD/cas_sd_selected.irp.f @@ -8,6 +8,7 @@ program full_ci N_st = N_states allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st)) character*(64) :: perturbation + PROVIDE N_det_cas pt2 = 1.d0 diag_algorithm = "Lapack" @@ -59,21 +60,21 @@ program full_ci logical :: in_cas integer :: exc_max, degree_min exc_max = 0 - print *, 'CAS determinants : ', N_det_generators - do i=1,min(N_det_generators,10) - 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) + print *, 'CAS determinants : ', N_det_cas + do i=1,min(N_det_cas,10) + do k=i,N_det_cas + call get_excitation_degree(psi_cas(1,1,k),psi_cas(1,1,i),degree,N_int) exc_max = max(exc_max,degree) enddo - call debug_det(psi_det_generators(1,1,i),N_int) + call debug_det(psi_cas(1,1,i),N_int) print *, '' enddo 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) + do k=1,N_det_cas + call get_excitation_degree(psi_cas(1,1,k),psi_det(1,1,i),degree,N_int) degree_min = min(degree_min,degree) enddo if (degree_min > 2) then diff --git a/src/MRCC/NEEDED_MODULES b/src/MRCC/NEEDED_MODULES index 31ed6f4e..f2bb9ba7 100644 --- a/src/MRCC/NEEDED_MODULES +++ b/src/MRCC/NEEDED_MODULES @@ -1,2 +1,2 @@ -AOs BiInts Bitmask CAS_SD Dets Electrons Ezfio_files Generators_CAS Generators_full Hartree_Fock MOGuess MonoInts MOs Nuclei Output Perturbation Properties Selectors_full Utils +AOs BiInts Bitmask Dets Electrons Ezfio_files Generators_full Hartree_Fock MOGuess MonoInts MOs Nuclei Output Perturbation Properties Selectors_full Utils diff --git a/src/MRCC/mrcc_utils.irp.f b/src/MRCC/mrcc_utils.irp.f index eca1a706..d33b7902 100644 --- a/src/MRCC/mrcc_utils.irp.f +++ b/src/MRCC/mrcc_utils.irp.f @@ -11,8 +11,12 @@ BEGIN_PROVIDER [ double precision, lambda_mrcc, (N_states,psi_det_size) ] size(psi_cas_coef,1), n_states, ihpsi) double precision :: hij do k=1,N_states - lambda_mrcc(k,i) = psi_non_cas_coef(i,k)/ihpsi(k) - lambda_mrcc(k,i) = min( lambda_mrcc (k,i),0.d0 ) + if (dabs(ihpsi(k)) > 1.d-5) then + lambda_mrcc(k,i) = psi_non_cas_coef(i,k)/ihpsi(k) + lambda_mrcc(k,i) = min( lambda_mrcc (k,i),0.d0 ) + else + lambda_mrcc(k,i) = 0.d0 + endif enddo enddo END_PROVIDER