diff --git a/scripts/generate_h_apply.py b/scripts/generate_h_apply.py index 4a5ad695..eeac4669 100755 --- a/scripts/generate_h_apply.py +++ b/scripts/generate_h_apply.py @@ -16,6 +16,8 @@ keys_work copy_buffer finalization generate_psi_guess +init_thread +deinit_thread """.split() class H_apply(object): @@ -42,7 +44,7 @@ class H_apply(object): !$OMP N_elec_in_key_hole_2,ia_ja_pairs,iproc) & !$OMP SHARED(key_in,N_int,elec_num_tab,mo_tot_num, & !$OMP hole_1, particl_1, hole_2, particl_2, & - !$OMP thresh,elec_alpha_num,i_generator)""" + !$OMP elec_alpha_num,i_generator)""" s["omp_end_parallel"] = "!$OMP END PARALLEL" s["omp_master"] = "!$OMP MASTER" s["omp_end_master"] = "!$OMP END MASTER" @@ -54,7 +56,7 @@ class H_apply(object): s["generate_psi_guess"] = """ ! Sort H_jj to find the N_states lowest states - integer :: i,k + integer :: i integer, allocatable :: iorder(:) double precision, allocatable :: H_jj(:) double precision, external :: diag_h_mat_elem @@ -81,9 +83,6 @@ class H_apply(object): for k in s: s[k] = "" s["size_max"] = str(1024*128) - s["set_i_H_j_threshold"] = """ - thresh = H_apply_threshold - """ s["copy_buffer"] = "call copy_h_apply_buffer_to_wf" self.data = s @@ -109,28 +108,41 @@ class H_apply(object): integer, intent(in) :: N_st,Nint double precision, intent(inout) :: sum_e_2_pert_in(N_st) double precision, intent(inout) :: sum_norm_pert_in(N_st) - double precision, intent(inout) :: sum_H_pert_diag_in + double precision, intent(inout) :: sum_H_pert_diag_in(N_st) double precision :: sum_e_2_pert(N_st) double precision :: sum_norm_pert(N_st) - double precision :: sum_H_pert_diag - double precision :: e_2_pert_buffer(N_st,size_max) - double precision :: coef_pert_buffer(N_st,size_max) + double precision :: sum_H_pert_diag(N_st) + double precision, allocatable :: e_2_pert_buffer(:,:) + double precision, allocatable :: coef_pert_buffer(:,:) + ASSERT (Nint == N_int) + """ + self.data["init_thread"] = """ + allocate (e_2_pert_buffer(N_st,size_max), coef_pert_buffer(N_st,size_max)) + do k=1,N_st + sum_e_2_pert(k) = 0.d0 + sum_norm_pert(k) = 0.d0 + sum_H_pert_diag(k) = 0.d0 + enddo + """ + self.data["deinit_thread"] = """ + !$OMP CRITICAL + do k=1,N_st + sum_e_2_pert_in(k) = sum_e_2_pert_in(k) + sum_e_2_pert(k) + sum_norm_pert_in(k) = sum_norm_pert_in(k) + sum_norm_pert(k) + sum_H_pert_diag_in(k) = sum_H_pert_diag_in(k) + sum_H_pert_diag(k) + enddo + !$OMP END CRITICAL + deallocate (e_2_pert_buffer, coef_pert_buffer) """ self.data["size_max"] = "256" self.data["initialization"] = """ - sum_e_2_pert = sum_e_2_pert_in - sum_norm_pert = sum_norm_pert_in - sum_H_pert_diag = sum_H_pert_diag_in PROVIDE CI_electronic_energy psi_selectors_coef psi_selectors """ self.data["keys_work"] = """ call perturb_buffer_%s(i_generator,keys_out,key_idx,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert, & - sum_norm_pert,sum_H_pert_diag,N_st,Nint) + sum_norm_pert,sum_H_pert_diag,N_st,N_int) """%(pert,) self.data["finalization"] = """ - sum_e_2_pert_in = sum_e_2_pert - sum_norm_pert_in = sum_norm_pert - sum_H_pert_diag_in = sum_H_pert_diag """ self.data["copy_buffer"] = "" self.data["generate_psi_guess"] = "" @@ -140,17 +152,19 @@ class H_apply(object): self.data["decls_main"] = """ integer, intent(in) :: N_st double precision, intent(inout):: pt2(N_st) double precision, intent(inout):: norm_pert(N_st) - double precision, intent(inout):: H_pert_diag + double precision, intent(inout):: H_pert_diag(N_st) PROVIDE CI_electronic_energy N_det_generators key_pattern_not_in_ref - pt2 = 0.d0 - norm_pert = 0.d0 - H_pert_diag = 0.d0 + do k=1,N_st + pt2(k) = 0.d0 + norm_pert(k) = 0.d0 + H_pert_diag(k) = 0.d0 + enddo """ if self.openmp: self.data["omp_parallel"] += """& - !$OMP SHARED(N_st,Nint) PRIVATE(e_2_pert_buffer,coef_pert_buffer) & - !$OMP REDUCTION(+:sum_e_2_pert, sum_norm_pert, sum_H_pert_diag)""" + !$OMP SHARED(N_st) PRIVATE(e_2_pert_buffer,coef_pert_buffer) & + !$OMP PRIVATE(sum_e_2_pert, sum_norm_pert, sum_H_pert_diag)""" def set_selection_pt2(self,pert): if self.selection_pt2 is not None: @@ -163,7 +177,7 @@ class H_apply(object): call copy_h_apply_buffer_to_wf selection_criterion_min = selection_criterion_min*0.1d0 selection_criterion = selection_criterion_min - call remove_small_contributions + !call remove_small_contributions """ self.data["keys_work"] = """ e_2_pert_buffer = 0.d0 diff --git a/setup_environment.sh b/setup_environment.sh index 9147d7e6..62b579ed 100755 --- a/setup_environment.sh +++ b/setup_environment.sh @@ -6,8 +6,6 @@ QPACKAGE_ROOT=${PWD} -IRPF90=$(which irpf90) - if [[ -z ${IRPF90} ]] ; then make irpf90 diff --git a/src/CISD/SC2.irp.f b/src/CISD/SC2.irp.f index efd85202..35f19e2b 100644 --- a/src/CISD/SC2.irp.f +++ b/src/CISD/SC2.irp.f @@ -143,7 +143,8 @@ end subroutine repeat_excitation(key_in,key_1,key_2,i_ok,Nint) use bitmasks implicit none - integer(bit_kind), intent(in) :: key_in(Nint,2),key_1(Nint,2),key_2(Nint,2),Nint + integer(bit_kind), intent(in) :: key_in(Nint,2),key_1(Nint,2),key_2(Nint,2) + integer :: Nint integer,intent(out) :: i_ok integer :: ispin,i_hole,k_hole,j_hole,i_particl,k_particl,j_particl,i_trou,degree,exc(0:2,2,2) double precision :: phase diff --git a/src/CISD_selected/H_apply.irp.f b/src/CISD_selected/H_apply.irp.f index 9b761963..1a01b53e 100644 --- a/src/CISD_selected/H_apply.irp.f +++ b/src/CISD_selected/H_apply.irp.f @@ -14,8 +14,7 @@ subroutine H_apply_cisd_selection(perturbation,pt2, norm_pert, H_pert_diag, N_st implicit none character*(64), intent(in) :: perturbation integer, intent(in) :: N_st - double precision, intent(inout):: pt2(N_st), norm_pert(N_st) - double precision,intent(inout) :: H_pert_diag + double precision, intent(inout):: pt2(N_st), norm_pert(N_st), H_pert_diag(N_st) BEGIN_SHELL [ /usr/bin/env python ] from perturbation import perturbations diff --git a/src/CISD_selected/README.rst b/src/CISD_selected/README.rst index 95499a8f..4fed26e5 100644 --- a/src/CISD_selected/README.rst +++ b/src/CISD_selected/README.rst @@ -8,6 +8,12 @@ Documentation .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. +`h_apply_cisd_selection `_ + Undocumented + +`cisd `_ + Undocumented + Needed Modules diff --git a/src/CISD_selected/cisd_selection.irp.f b/src/CISD_selected/cisd_selection.irp.f index 40691065..6d88f61e 100644 --- a/src/CISD_selected/cisd_selection.irp.f +++ b/src/CISD_selected/cisd_selection.irp.f @@ -3,27 +3,25 @@ program cisd integer :: i,k - double precision, allocatable :: pt2(:), norm_pert(:) - double precision :: H_pert_diag, E_old + double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:) integer :: N_st, iter character*(64) :: perturbation N_st = N_states - allocate (pt2(N_st), norm_pert(N_st)) + allocate (pt2(N_st), norm_pert(N_st), H_pert_diag(N_st)) pt2 = 1.d0 perturbation = "epstein_nesbet" do while (maxval(abs(pt2(1:N_st))) > 1.d-6) - E_old = CI_energy(1) call H_apply_cisd_selection(perturbation,pt2, norm_pert, H_pert_diag, N_st) call diagonalize_CI print *, 'N_det = ', N_det print *, 'N_states = ', N_states print *, 'PT2 = ', pt2 - print *, 'E = ', E_old - print *, 'E+PT2 = ', E_old+pt2 + print *, 'E = ', CI_energy + print *, 'E+PT2 = ', CI_energy+pt2 if (abort_all) then exit endif enddo - deallocate(pt2,norm_pert) + deallocate(pt2,norm_pert,H_pert_diag) end diff --git a/src/Dets/H_apply.irp.f b/src/Dets/H_apply.irp.f index 3a4922fb..be94a7a0 100644 --- a/src/Dets/H_apply.irp.f +++ b/src/Dets/H_apply.irp.f @@ -19,7 +19,7 @@ BEGIN_PROVIDER [ logical, H_apply_buffer_allocated ] ! Uninitialized. Filled by H_apply subroutines. END_DOC integer :: iproc, sze - sze = 10 + sze = 100 if (.not.associated(H_apply_buffer)) then allocate(H_apply_buffer(0:nproc-1)) iproc = 0 diff --git a/src/Dets/H_apply_template.f b/src/Dets/H_apply_template.f index adeff13d..3b6f1e35 100644 --- a/src/Dets/H_apply_template.f +++ b/src/Dets/H_apply_template.f @@ -14,31 +14,33 @@ 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(bit_kind) :: hole_save(N_int,2) - integer(bit_kind) :: key(N_int,2),hole(N_int,2), particle(N_int,2) - integer(bit_kind) :: hole_tmp(N_int,2), particle_tmp(N_int,2) + integer(bit_kind), allocatable :: hole_save(:,:) + integer(bit_kind), allocatable :: key(:,:),hole(:,:), particle(:,:) + integer(bit_kind), allocatable :: hole_tmp(:,:), particle_tmp(:,:) integer :: ii,i,jj,j,k,ispin,l - integer :: occ_particle(N_int*bit_kind_size,2) - integer :: occ_hole(N_int*bit_kind_size,2) - integer :: occ_particle_tmp(N_int*bit_kind_size,2) - integer :: occ_hole_tmp(N_int*bit_kind_size,2) + integer, allocatable :: occ_particle(:,:), occ_hole(:,:) + integer, allocatable :: occ_particle_tmp(:,:), occ_hole_tmp(:,:) integer :: kk,pp,other_spin,key_idx integer :: N_elec_in_key_hole_1(2),N_elec_in_key_part_1(2) integer :: N_elec_in_key_hole_2(2),N_elec_in_key_part_2(2) - double precision :: mo_bielec_integral, thresh + double precision :: mo_bielec_integral integer, allocatable :: ia_ja_pairs(:,:,:) double precision :: diag_H_mat_elem integer :: iproc - - $set_i_H_j_threshold + PROVIDE H_apply_threshold $initialization iproc = 0 $omp_parallel !$ iproc = omp_get_thread_num() - allocate (keys_out(N_int,2,size_max)) + 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), & + occ_hole(N_int*bit_kind_size,2), occ_particle_tmp(N_int*bit_kind_size,2),& + occ_hole_tmp(N_int*bit_kind_size,2)) + $init_thread !print*,'key_in !!' !call print_key(key_in) @@ -142,7 +144,6 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_gene j_b = occ_particle_tmp(jjj,other_spin) ASSERT (j_b > 0) ASSERT (j_b <= mo_tot_num) - if(dabs( mo_bielec_integral(j_a,j_b,i_a,i_b)) 0) ASSERT (j_b <= mo_tot_num) if (j_b <= j_a) cycle - if(dabs( mo_bielec_integral(j_a,j_b,i_b,i_a)) 5) then cycle + else + connected_to_ref = i + return endif enddo @@ -64,6 +67,9 @@ integer function connected_to_ref(key,keys,Nint,N_past_in,Ndet,thresh) endif if (degree_x2 > 5) then cycle + else + connected_to_ref = i + return endif enddo @@ -101,6 +107,9 @@ integer function connected_to_ref(key,keys,Nint,N_past_in,Ndet,thresh) endif if (degree_x2 > 5) then cycle + else + connected_to_ref = i + return endif enddo @@ -140,6 +149,9 @@ integer function connected_to_ref(key,keys,Nint,N_past_in,Ndet,thresh) endif if (degree_x2 > 5) then cycle + else + connected_to_ref = i + return endif enddo diff --git a/src/Dets/davidson.irp.f b/src/Dets/davidson.irp.f index 26fa1195..2bc500d2 100644 --- a/src/Dets/davidson.irp.f +++ b/src/Dets/davidson.irp.f @@ -374,8 +374,8 @@ end BEGIN_DOC ! Can be : [ energy | residual | both | wall_time | cpu_time | iterations ] END_DOC - davidson_criterion = 'both' - davidson_threshold = 1.d-8 + davidson_criterion = 'residual' + davidson_threshold = 1.d-6 END_PROVIDER subroutine davidson_converged(energy,residual,wall,iterations,cpu,N_st,converged) diff --git a/src/Dets/determinants.irp.f b/src/Dets/determinants.irp.f index 266898da..b33bf664 100644 --- a/src/Dets/determinants.irp.f +++ b/src/Dets/determinants.irp.f @@ -57,18 +57,18 @@ END_PROVIDER ifirst = 1 psi_det = 0_bit_kind psi_coef = 0.d0 - - integer :: i - do i=1,N_int - psi_det(i,1,1) = HF_bitmask(i,1) - psi_det(i,2,1) = HF_bitmask(i,2) - enddo - - do i=1,N_states - psi_coef(i,i) = 1.d0 - enddo endif + integer :: i + do i=1,N_int + psi_det(i,1,1) = HF_bitmask(i,1) + psi_det(i,2,1) = HF_bitmask(i,2) + enddo + + do i=1,N_states + psi_coef(i,i) = 1.d0 + enddo + END_PROVIDER @@ -80,3 +80,54 @@ BEGIN_PROVIDER [ integer, N_det_reference ] N_det_reference = N_det ASSERT (N_det_reference > 0) END_PROVIDER + +BEGIN_PROVIDER [ double precision, psi_average_norm_contrib, (N_det) ] + implicit none + BEGIN_DOC + ! Contribution of determinants to the state-averaged density + END_DOC + integer :: i,j,k + double precision :: f + f = 1.d0/dble(N_states) + do i=1,N_det + psi_average_norm_contrib(i) = psi_coef(i,1)*psi_coef(i,1)*f + enddo + do k=2,N_states + do i=1,N_det + psi_average_norm_contrib(i) = psi_average_norm_contrib(i) + & + psi_coef(i,k)*psi_coef(i,k)*f + enddo + enddo +END_PROVIDER + + BEGIN_PROVIDER [ integer(bit_kind), psi_det_sorted, (N_int,2,N_det) ] +&BEGIN_PROVIDER [ double precision, psi_coef_sorted, (N_det,N_states) ] +&BEGIN_PROVIDER [ double precision, psi_average_norm_contrib_sorted, (N_det) ] + implicit none + BEGIN_DOC + ! Wave function sorted by determinants (state-averaged) + END_DOC + integer :: i,j,k + integer, allocatable :: iorder(:) + allocate ( iorder(N_det) ) + do i=1,N_det + psi_average_norm_contrib_sorted(i) = -psi_average_norm_contrib(i) + iorder(i) = i + enddo + call dsort(psi_average_norm_contrib_sorted,iorder,N_det) + !DIR$ IVDEP + do i=1,N_det + do j=1,N_int + psi_det_sorted(j,1,i) = psi_det(j,1,iorder(i)) + psi_det_sorted(j,2,i) = psi_det(j,2,iorder(i)) + enddo + do k=1,N_states + psi_coef_sorted(i,k) = psi_coef(iorder(i),k) + enddo + psi_average_norm_contrib_sorted(i) = -psi_average_norm_contrib_sorted(i) + enddo + + deallocate(iorder) + +END_PROVIDER + diff --git a/src/Dets/diagonalize_CI.irp.f b/src/Dets/diagonalize_CI.irp.f index 1f1ad1cf..0dc39448 100644 --- a/src/Dets/diagonalize_CI.irp.f +++ b/src/Dets/diagonalize_CI.irp.f @@ -22,8 +22,12 @@ BEGIN_PROVIDER [ double precision, CI_energy, (N_states) ] END_DOC integer :: j + character*(8) :: st + call write_time(output_Dets) do j=1,N_states CI_energy(j) = CI_electronic_energy(j) + nuclear_repulsion + write(st,'(I4)') j + call write_double(output_Dets,CI_energy(j),'Energy of state '//trim(st)) enddo END_PROVIDER diff --git a/src/Full_CI/ASSUMPTIONS.rst b/src/Full_CI/ASSUMPTIONS.rst new file mode 100644 index 00000000..e69de29b diff --git a/src/Full_CI/H_apply.irp.f b/src/Full_CI/H_apply.irp.f new file mode 100644 index 00000000..c9b44cbd --- /dev/null +++ b/src/Full_CI/H_apply.irp.f @@ -0,0 +1,10 @@ +use bitmasks +BEGIN_SHELL [ /usr/bin/env python ] +from generate_h_apply import * + +s = H_apply("FCI",openmp=True) +s.set_selection_pt2("epstein_nesbet_2x2") +print s + +END_SHELL + diff --git a/src/Full_CI/Makefile b/src/Full_CI/Makefile new file mode 100644 index 00000000..b2ea1de1 --- /dev/null +++ b/src/Full_CI/Makefile @@ -0,0 +1,8 @@ +default: all + +# 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/Full_CI/NEEDED_MODULES b/src/Full_CI/NEEDED_MODULES new file mode 100644 index 00000000..8a732901 --- /dev/null +++ b/src/Full_CI/NEEDED_MODULES @@ -0,0 +1 @@ +AOs BiInts Bitmask Dets Electrons Ezfio_files Generators_full Hartree_Fock MonoInts MOs Nuclei Output Selectors_full Utils Perturbation diff --git a/src/Full_CI/README.rst b/src/Full_CI/README.rst new file mode 100644 index 00000000..5f0ae69c --- /dev/null +++ b/src/Full_CI/README.rst @@ -0,0 +1,39 @@ +============== +Full_CI Module +============== + +Performs a perturbatively selected Full-CI. + +Documentation +============= + +.. Do not edit this section. It was auto-generated from the +.. NEEDED_MODULES file. + +`cisd `_ + Undocumented + + + +Needed Modules +============== + +.. Do not edit this section. It was auto-generated from the +.. NEEDED_MODULES file. + +* `AOs `_ +* `BiInts `_ +* `Bitmask `_ +* `Dets `_ +* `Electrons `_ +* `Ezfio_files `_ +* `Generators_full `_ +* `Hartree_Fock `_ +* `MonoInts `_ +* `MOs `_ +* `Nuclei `_ +* `Output `_ +* `Selectors_full `_ +* `Utils `_ +* `Perturbation `_ + diff --git a/src/Full_CI/full_ci.irp.f b/src/Full_CI/full_ci.irp.f new file mode 100644 index 00000000..72eec40e --- /dev/null +++ b/src/Full_CI/full_ci.irp.f @@ -0,0 +1,28 @@ +program cisd + implicit none + integer :: i,k + + + double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:) + integer :: N_st, degree + character*(64) :: perturbation + N_st = N_states + allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st)) + + pt2 = 1.d0 + diag_algorithm = "Lapack" + do while (maxval(abs(pt2(1:N_st))) > 1.d-4) + call H_apply_FCI(pt2, norm_pert, H_pert_diag, N_st) + call diagonalize_CI + print *, 'N_det = ', N_det + print *, 'N_states = ', N_states + print *, 'PT2 = ', pt2 + print *, 'E = ', CI_energy + print *, 'E+PT2 = ', CI_energy+pt2 + print *, '-----' + if (abort_all) then + exit + endif + enddo + deallocate(pt2,norm_pert) +end diff --git a/src/Generators_full/ASSUMPTIONS.rst b/src/Generators_full/ASSUMPTIONS.rst new file mode 100644 index 00000000..e69de29b diff --git a/src/Generators_full/Makefile b/src/Generators_full/Makefile new file mode 100644 index 00000000..b2ea1de1 --- /dev/null +++ b/src/Generators_full/Makefile @@ -0,0 +1,8 @@ +default: all + +# 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/Generators_full/NEEDED_MODULES b/src/Generators_full/NEEDED_MODULES new file mode 100644 index 00000000..f6551dd9 --- /dev/null +++ b/src/Generators_full/NEEDED_MODULES @@ -0,0 +1 @@ +AOs BiInts Bitmask Dets Electrons Ezfio_files Hartree_Fock MonoInts MOs Nuclei Output Utils diff --git a/src/Generators_full/README.rst b/src/Generators_full/README.rst new file mode 100644 index 00000000..e3ba2312 --- /dev/null +++ b/src/Generators_full/README.rst @@ -0,0 +1,7 @@ +====================== +Generators_full Module +====================== + +All the determinants of the wave function are generators. In this way, the Full CI +space is explored. + diff --git a/src/Generators_full/generators.irp.f b/src/Generators_full/generators.irp.f new file mode 100644 index 00000000..00930ae6 --- /dev/null +++ b/src/Generators_full/generators.irp.f @@ -0,0 +1,44 @@ +use bitmasks + +BEGIN_PROVIDER [ double precision, generators_threshold ] + implicit none + BEGIN_DOC + ! Percentage of the norm of the state-averaged wave function to + ! consider for the generators + END_DOC + generators_threshold = 0.97d0 +END_PROVIDER + +BEGIN_PROVIDER [ integer, N_det_generators ] + implicit none + BEGIN_DOC + ! For Single reference wave functions, the number of generators is 1 : the + ! Hartree-Fock determinant + END_DOC + integer :: i + double precision :: norm + call write_time(output_dets) + norm = 0.d0 + N_det_generators = N_det + do i=1,N_det + norm = norm + psi_average_norm_contrib_sorted(i) + if (norm > generators_threshold) then + N_det_generators = i-1 + exit + endif + enddo + N_det_generators = max(N_det_generators,1) + 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) ] + implicit none + BEGIN_DOC + ! For Single reference wave functions, the generator is the + ! Hartree-Fock determinant + END_DOC + psi_generators = psi_det_sorted + +END_PROVIDER + + diff --git a/src/Generators_full/tests/Makefile b/src/Generators_full/tests/Makefile new file mode 100644 index 00000000..77bd84ba --- /dev/null +++ b/src/Generators_full/tests/Makefile @@ -0,0 +1,33 @@ +OPENMP =1 +PROFILE =0 +DEBUG = 0 + +IRPF90+= -I tests + +REF_FILES=$(subst %.irp.f, %.ref, $(wildcard *.irp.f)) + +.PHONY: clean executables serial_tests parallel_tests + +all: clean executables serial_tests parallel_tests + +parallel_tests: $(REF_FILES) + @echo ; echo " ---- Running parallel tests ----" ; echo + @OMP_NUM_THREADS=10 ${QPACKAGE_ROOT}/scripts/run_tests.py + +serial_tests: $(REF_FILES) + @echo ; echo " ---- Running serial tests ----" ; echo + @OMP_NUM_THREADS=1 ${QPACKAGE_ROOT}/scripts/run_tests.py + +executables: $(wildcard *.irp.f) veryclean + $(MAKE) -C .. + +%.ref: $(wildcard $(QPACKAGE_ROOT)/data/inputs/*.md5) executables + $(QPACKAGE_ROOT)/scripts/create_test_ref.sh $* + +clean: + $(MAKE) -C .. clean + +veryclean: + $(MAKE) -C .. veryclean + + diff --git a/src/NEEDED_MODULES b/src/NEEDED_MODULES index 725707f3..6151539c 100644 --- a/src/NEEDED_MODULES +++ b/src/NEEDED_MODULES @@ -1 +1 @@ -AOs Bitmask Electrons Ezfio_files MOs Nuclei Output Utils Hartree_Fock BiInts MonoInts MOGuess Dets DensityMatrix CISD Perturbation SingleRefMethod CISD_selected Selectors_full MP2 +AOs Bitmask Electrons Ezfio_files MOs Nuclei Output Utils Hartree_Fock BiInts MonoInts MOGuess Dets DensityMatrix CISD Perturbation SingleRefMethod CISD_selected Selectors_full MP2 Generators_full Full_CI diff --git a/src/Perturbation/Moller_plesset.irp.f b/src/Perturbation/Moller_plesset.irp.f index bca48358..7435f70c 100644 --- a/src/Perturbation/Moller_plesset.irp.f +++ b/src/Perturbation/Moller_plesset.irp.f @@ -3,7 +3,7 @@ subroutine pt2_moller_plesset(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_s implicit none integer, intent(in) :: Nint,ndet,n_st integer(bit_kind), intent(in) :: det_pert(Nint,2) - double precision , intent(out) :: c_pert(n_st),e_2_pert(n_st),H_pert_diag + double precision , intent(out) :: c_pert(n_st),e_2_pert(n_st),H_pert_diag(N_st) double precision :: i_H_psi_array(N_st) BEGIN_DOC @@ -21,7 +21,7 @@ subroutine pt2_moller_plesset(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_s double precision :: diag_H_mat_elem integer :: exc(0:2,2,2) integer :: degree - double precision :: phase,delta_e + double precision :: phase,delta_e,h integer :: h1,h2,p1,p2,s1,s2 ASSERT (Nint == N_int) ASSERT (Nint > 0) @@ -32,8 +32,9 @@ subroutine pt2_moller_plesset(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_s delta_e = 1.d0/delta_e call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det,psi_selectors_size,n_st,i_H_psi_array) - H_pert_diag = diag_H_mat_elem(det_pert,Nint) + h = diag_H_mat_elem(det_pert,Nint) do i =1,n_st + H_pert_diag(i) = h c_pert(i) = i_H_psi_array(i) *delta_e e_2_pert(i) = c_pert(i) * i_H_psi_array(i) enddo diff --git a/src/Perturbation/epstein_nesbet.irp.f b/src/Perturbation/epstein_nesbet.irp.f index b609cd02..d6a903c7 100644 --- a/src/Perturbation/epstein_nesbet.irp.f +++ b/src/Perturbation/epstein_nesbet.irp.f @@ -1,15 +1,15 @@ -subroutine pt2_epstein_nesbet(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_st) +subroutine pt2_epstein_nesbet(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st) use bitmasks implicit none - integer, intent(in) :: Nint,ndet,n_st + integer, intent(in) :: Nint,ndet,N_st integer(bit_kind), intent(in) :: det_pert(Nint,2) - double precision , intent(out) :: c_pert(n_st),e_2_pert(n_st),H_pert_diag + double precision , intent(out) :: c_pert(N_st),e_2_pert(N_st),H_pert_diag(N_st) double precision :: i_H_psi_array(N_st) BEGIN_DOC ! compute the standard Epstein-Nesbet perturbative first order coefficient and second order energetic contribution ! - ! for the various n_st states. + ! for the various N_st states. ! ! c_pert(i) = /( E(i) - ) ! @@ -18,14 +18,15 @@ subroutine pt2_epstein_nesbet(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_s END_DOC integer :: i,j - double precision :: diag_H_mat_elem + double precision :: diag_H_mat_elem, h ASSERT (Nint == N_int) ASSERT (Nint > 0) call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array) - H_pert_diag = diag_H_mat_elem(det_pert,Nint) - do i =1,n_st - if (dabs(CI_electronic_energy(i) - H_pert_diag) > 1.d-6) then - c_pert(i) = i_H_psi_array(i) / (CI_electronic_energy(i) - H_pert_diag) + h = diag_H_mat_elem(det_pert,Nint) + do i =1,N_st + H_pert_diag(i) = h + if (dabs(CI_electronic_energy(i) - h) > 1.d-6) then + c_pert(i) = i_H_psi_array(i) / (CI_electronic_energy(i) - h) e_2_pert(i) = c_pert(i) * i_H_psi_array(i) else c_pert(i) = 0.d0 @@ -35,18 +36,18 @@ subroutine pt2_epstein_nesbet(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_s end -subroutine pt2_epstein_nesbet_2x2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_st) +subroutine pt2_epstein_nesbet_2x2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st) use bitmasks implicit none - integer, intent(in) :: Nint,ndet,n_st + integer, intent(in) :: Nint,ndet,N_st integer(bit_kind), intent(in) :: det_pert(Nint,2) - double precision , intent(out) :: c_pert(n_st),e_2_pert(n_st),H_pert_diag + double precision , intent(out) :: c_pert(N_st),e_2_pert(N_st),H_pert_diag(N_st) double precision :: i_H_psi_array(N_st) BEGIN_DOC ! compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution ! - ! for the various n_st states. + ! for the various N_st states. ! ! e_2_pert(i) = 0.5 * (( - E(i) ) - sqrt( ( - E(i)) ^2 + 4 ^2 ) ! @@ -55,13 +56,14 @@ subroutine pt2_epstein_nesbet_2x2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet END_DOC integer :: i,j - double precision :: diag_H_mat_elem,delta_e + double precision :: diag_H_mat_elem,delta_e, h ASSERT (Nint == N_int) ASSERT (Nint > 0) - call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,n_st,i_H_psi_array) - H_pert_diag = diag_H_mat_elem(det_pert,Nint) - do i =1,n_st - delta_e = H_pert_diag - CI_electronic_energy(i) + call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array) + h = diag_H_mat_elem(det_pert,Nint) + do i =1,N_st + H_pert_diag(i) = h + delta_e = h - CI_electronic_energy(i) e_2_pert(i) = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * i_H_psi_array(i) * i_H_psi_array(i))) if (dabs(i_H_psi_array(i)) > 1.d-6) then c_pert(i) = e_2_pert(i)/i_H_psi_array(i) @@ -72,19 +74,19 @@ subroutine pt2_epstein_nesbet_2x2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet end -subroutine pt2_epstein_nesbet_SC2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_st) +subroutine pt2_epstein_nesbet_SC2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st) use bitmasks implicit none - integer, intent(in) :: Nint,ndet,n_st + integer, intent(in) :: Nint,ndet,N_st integer(bit_kind), intent(in) :: det_pert(Nint,2) - double precision , intent(out) :: c_pert(n_st),e_2_pert(n_st),H_pert_diag + double precision , intent(out) :: c_pert(N_st),e_2_pert(N_st),H_pert_diag(N_st) double precision :: i_H_psi_array(N_st) integer :: idx_repeat(0:ndet) BEGIN_DOC ! compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution ! - ! for the various n_st states, + ! for the various N_st states, ! ! but with the correction in the denominator ! @@ -105,40 +107,41 @@ subroutine pt2_epstein_nesbet_SC2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet END_DOC integer :: i,j - double precision :: diag_H_mat_elem,accu_e_corr,hij + double precision :: diag_H_mat_elem,accu_e_corr,hij,h ASSERT (Nint == N_int) ASSERT (Nint > 0) - call i_H_psi_SC2(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,n_st,i_H_psi_array,idx_repeat) + call i_H_psi_SC2(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array,idx_repeat) accu_e_corr = 0.d0 do i = 1, idx_repeat(0) call i_H_j(psi_selectors(1,1,idx_repeat(i)),det_pert,Nint,hij) accu_e_corr = accu_e_corr + hij * psi_selectors_coef(idx_repeat(i),1) enddo accu_e_corr = accu_e_corr / psi_selectors_coef(1,1) - H_pert_diag = diag_H_mat_elem(det_pert,Nint) + accu_e_corr - do i =1,n_st + h = diag_H_mat_elem(det_pert,Nint) + accu_e_corr + do i =1,N_st + H_pert_diag(i) = h e_2_pert(i) = c_pert(i) * i_H_psi_array(i) - if (dabs(CI_electronic_energy(i) - H_pert_diag) > 1.d-6) then - c_pert(i) = i_H_psi_array(i) / (CI_electronic_energy(i) - H_pert_diag) + if (dabs(CI_electronic_energy(i) - h) > 1.d-6) then + c_pert(i) = i_H_psi_array(i) / (CI_electronic_energy(i) - h) else c_pert(i) = 0.d0 endif enddo end -subroutine pt2_epstein_nesbet_2x2_SC2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_st) +subroutine pt2_epstein_nesbet_2x2_SC2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st) use bitmasks implicit none - integer, intent(in) :: Nint,ndet,n_st + integer, intent(in) :: Nint,ndet,N_st integer(bit_kind), intent(in) :: det_pert(Nint,2) - double precision , intent(out) :: c_pert(n_st),e_2_pert(n_st),H_pert_diag + double precision , intent(out) :: c_pert(N_st),e_2_pert(N_st),H_pert_diag(N_st) double precision :: i_H_psi_array(N_st) integer :: idx_repeat(0:ndet) BEGIN_DOC ! compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution ! - ! for the various n_st states. + ! for the various N_st states. ! ! but with the correction in the denominator ! @@ -159,19 +162,20 @@ subroutine pt2_epstein_nesbet_2x2_SC2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint, END_DOC integer :: i,j - double precision :: diag_H_mat_elem,accu_e_corr,hij,delta_e + double precision :: diag_H_mat_elem,accu_e_corr,hij,delta_e,h ASSERT (Nint == N_int) ASSERT (Nint > 0) - call i_H_psi_SC2(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,n_st,i_H_psi_array,idx_repeat) + call i_H_psi_SC2(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array,idx_repeat) accu_e_corr = 0.d0 do i = 1, idx_repeat(0) call i_H_j(psi_selectors(1,1,idx_repeat(i)),det_pert,Nint,hij) accu_e_corr = accu_e_corr + hij * psi_selectors_coef(idx_repeat(i),1) enddo accu_e_corr = accu_e_corr / psi_selectors_coef(1,1) - H_pert_diag = diag_H_mat_elem(det_pert,Nint) + accu_e_corr - do i =1,n_st - delta_e = H_pert_diag - CI_electronic_energy(i) + h = diag_H_mat_elem(det_pert,Nint) + accu_e_corr + do i =1,N_st + H_pert_diag(i) = h + delta_e = h - CI_electronic_energy(i) e_2_pert(i) = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * i_H_psi_array(i) * i_H_psi_array(i))) if (dabs(i_H_psi_array(i)) > 1.d-6) then c_pert(i) = e_2_pert(i)/i_H_psi_array(i) @@ -182,19 +186,19 @@ subroutine pt2_epstein_nesbet_2x2_SC2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint, end -subroutine pt2_epstein_nesbet_SC2_projected(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_st) +subroutine pt2_epstein_nesbet_SC2_projected(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st) use bitmasks implicit none - integer, intent(in) :: Nint,ndet,n_st + integer, intent(in) :: Nint,ndet,N_st integer(bit_kind), intent(in) :: det_pert(Nint,2) - double precision , intent(out) :: c_pert(n_st),e_2_pert(n_st),H_pert_diag + double precision , intent(out) :: c_pert(N_st),e_2_pert(N_st),H_pert_diag(N_st) double precision :: i_H_psi_array(N_st) integer :: idx_repeat(0:ndet) BEGIN_DOC ! compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution ! - ! for the various n_st states, + ! for the various N_st states, ! ! but with the correction in the denominator ! @@ -219,10 +223,10 @@ subroutine pt2_epstein_nesbet_SC2_projected(det_pert,c_pert,e_2_pert,H_pert_diag END_DOC integer :: i,j - double precision :: diag_H_mat_elem,accu_e_corr,hij,h0j + double precision :: diag_H_mat_elem,accu_e_corr,hij,h0j,h ASSERT (Nint == N_int) ASSERT (Nint > 0) - call i_H_psi_SC2(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,n_st,i_H_psi_array,idx_repeat) + call i_H_psi_SC2(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array,idx_repeat) accu_e_corr = 0.d0 call i_H_j(ref_bitmask,det_pert,Nint,h0j) do i = 1, idx_repeat(0) @@ -230,13 +234,14 @@ subroutine pt2_epstein_nesbet_SC2_projected(det_pert,c_pert,e_2_pert,H_pert_diag accu_e_corr = accu_e_corr + hij * psi_selectors_coef(idx_repeat(i),1) enddo accu_e_corr = accu_e_corr / psi_selectors_coef(1,1) - H_pert_diag = diag_H_mat_elem(det_pert,Nint) + accu_e_corr + h = diag_H_mat_elem(det_pert,Nint) + accu_e_corr - c_pert(1) = 1.d0/psi_selectors_coef(1,1) * i_H_psi_array(1) / (CI_electronic_energy(i) - H_pert_diag) + c_pert(1) = 1.d0/psi_selectors_coef(1,1) * i_H_psi_array(1) / (CI_electronic_energy(i) - h) e_2_pert(1) = c_pert(i) * h0j - do i =2,n_st - if (dabs(CI_electronic_energy(i) - H_pert_diag) > 1.d-6) then - c_pert(i) = i_H_psi_array(i) / (CI_electronic_energy(i) - H_pert_diag) + do i =2,N_st + H_pert_diag(i) = h + if (dabs(CI_electronic_energy(i) - h) > 1.d-6) then + c_pert(i) = i_H_psi_array(i) / (CI_electronic_energy(i) - h) else c_pert(i) = 0.d0 endif diff --git a/src/Perturbation/perturbation_template.f b/src/Perturbation/perturbation_template.f index 55aa7f56..e3033820 100644 --- a/src/Perturbation/perturbation_template.f +++ b/src/Perturbation/perturbation_template.f @@ -13,7 +13,7 @@ subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,c integer(bit_kind), intent(in) :: buffer(Nint,2,buffer_size) double precision, intent(inout) :: sum_norm_pert(N_st),sum_e_2_pert(N_st) double precision, intent(inout) :: coef_pert_buffer(N_st,buffer_size),e_2_pert_buffer(N_st,buffer_size),sum_H_pert_diag(N_st) - double precision :: c_pert(N_st), e_2_pert(N_st), H_pert_diag + double precision :: c_pert(N_st), e_2_pert(N_st), H_pert_diag(N_st) integer :: i,k, c_ref integer :: connected_to_ref @@ -38,7 +38,7 @@ subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,c coef_pert_buffer(k,i) = c_pert(k) sum_norm_pert(k) += c_pert(k) * c_pert(k) sum_e_2_pert(k) += e_2_pert(k) - sum_H_pert_diag(k) += c_pert(k) * c_pert(k) * H_pert_diag + sum_H_pert_diag(k) += c_pert(k) * c_pert(k) * H_pert_diag(k) enddo enddo diff --git a/src/Perturbation/selection.irp.f b/src/Perturbation/selection.irp.f index 25f6cf37..203dccec 100644 --- a/src/Perturbation/selection.irp.f +++ b/src/Perturbation/selection.irp.f @@ -31,13 +31,12 @@ subroutine fill_H_apply_buffer_selection(n_selected,det_buffer,e_2_pert_buffer,c l=H_apply_buffer(iproc)%N_det do i=1,n_selected - s = 0.d0 + is_selected = .False. do j=1,N_st - s -= e_2_pert_buffer(j,i) + s = -e_2_pert_buffer(j,i) + is_selected = s > selection_criterion*selection_criterion_factor .or. is_selected enddo - ASSERT (s>=-1.d-8) - is_selected = s > selection_criterion * selection_criterion_factor if (is_selected) then l = l+1 @@ -60,8 +59,10 @@ subroutine fill_H_apply_buffer_selection(n_selected,det_buffer,e_2_pert_buffer,c ASSERT (sum(popcnt(h_apply_buffer(iproc)%det(:,1,i)) )== elec_alpha_num) ASSERT (sum(popcnt(h_apply_buffer(iproc)%det(:,2,i))) == elec_beta_num) enddo - selection_criterion = smax - selection_criterion_min = smin + !$OMP CRITICAL + selection_criterion = max(selection_criterion,smax) + selection_criterion_min = min(selection_criterion_min,smin) + !$OMP END CRITICAL end BEGIN_PROVIDER [ double precision, selection_criterion ] @@ -71,7 +72,7 @@ end BEGIN_DOC ! Threshold to select determinants. Set by selection routines. END_DOC - selection_criterion = .1d0 + selection_criterion = 1.d0 selection_criterion_factor = 0.01d0 selection_criterion_min = selection_criterion @@ -80,33 +81,52 @@ END_PROVIDER subroutine remove_small_contributions implicit none BEGIN_DOC -! Remove determinants with small contributions +! Remove determinants with small contributions. N_states is assumed to be +! provided. END_DOC integer :: i,j,k, N_removed - logical keep - N_removed = 0 - do i=N_det,1,-1 - keep = .False. + logical, allocatable :: keep(:) + double precision :: i_H_psi_array(N_states) + allocate (keep(N_det)) + call diagonalize_CI + do i=1,N_det + keep(i) = .True. + enddo + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(i,j,i_H_psi_array) & + !$OMP SHARED(k,psi_det_sorted,psi_coef_sorted,N_int,N_det,psi_det_size,N_states, & + !$OMP selection_criterion_min,keep,N_det_generators) & + !$OMP REDUCTION(+:N_removed) + !$OMP DO + do i=2*N_det_generators+1, N_det + call i_H_psi(psi_det_sorted(1,1,i),psi_det_sorted,psi_coef_sorted,N_int,min(N_det,2*N_det_generators),psi_det_size,N_states,i_H_psi_array) + keep(i) = .False. do j=1,N_states - keep = keep .or. (dabs(psi_coef(i,j)) > selection_criterion_min) + keep(i) = keep(i) .or. (-(psi_coef_sorted(i,j)*i_H_psi_array(j)) > selection_criterion_min) enddo - if (.not.keep) then - do k=i+1,N_det - do j=1,N_int - psi_det(j,1,k-1) = psi_det(j,1,k) - psi_det(j,2,k-1) = psi_det(j,2,k) - enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + N_removed = 0 + k = 0 + do i=1, N_det + if (keep(i)) then + k += 1 + do j=1,N_int + psi_det(j,1,k) = psi_det_sorted(j,1,i) + psi_det(j,2,k) = psi_det_sorted(j,2,i) enddo do j=1,N_states - do k=i+1,N_det - psi_coef(k-1,j) = psi_coef(k,j) - enddo + psi_coef(k,j) = psi_coef_sorted(i,j) enddo + else N_removed += 1 endif enddo + deallocate(keep) if (N_removed > 0) then - N_det -= N_removed + N_det = N_det - N_removed + SOFT_TOUCH N_det psi_det psi_coef call write_int(output_dets,N_removed, 'Removed determinants') endif end