From 5c957cf1f25fedefc2afe054a824e0065ce97134 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 19 Mar 2015 21:14:52 +0100 Subject: [PATCH] Merged eginer-master --- data/ezfio_defaults | 13 + scripts/generate_h_apply.py | 39 +- src/BiInts/ao_bi_integrals.irp.f | 3 +- src/BiInts/map_integrals.irp.f | 3 - src/BiInts/mo_bi_integrals.irp.f | 1 + src/CAS_SD_selected/NEEDED_MODULES | 3 +- src/CID_SC2_selected/NEEDED_MODULES | 2 +- src/CID_SC2_selected/cid_sc2_selection.irp.f | 1 - src/CID_selected/NEEDED_MODULES | 2 +- src/CISD_SC2_selected/H_apply.irp.f | 4 + src/CISD_SC2_selected/NEEDED_MODULES | 2 +- .../cisd_sc2_selection.irp.f | 6 +- src/CISD_selected/NEEDED_MODULES | 2 +- src/CISD_selected/cisd.ezfio_config | 3 + src/CISD_selected/cisd_selection.irp.f | 8 +- src/CISD_selected/options.irp.f | 34 ++ src/DDCI_selected/NEEDED_MODULES | 3 +- src/Dets/H_apply_template.f | 18 +- src/Dets/SC2.irp.f | 4 +- src/Dets/connected_to_ref.irp.f | 108 ++++++ src/Dets/density_matrix.irp.f | 133 ++++++- src/Dets/determinants.ezfio_config | 1 + src/Dets/determinants.irp.f | 38 +- src/Dets/diagonalize_CI.irp.f | 2 - src/Dets/diagonalize_CI_SC2.irp.f | 3 +- src/Dets/diagonalize_CI_mono.irp.f | 72 ++++ src/Dets/excitations_utils.irp.f | 16 + src/Dets/guess_doublet.irp.f | 79 ++++ src/Dets/guess_singlet.irp.f | 44 +++ src/Dets/guess_triplet.irp.f | 48 +++ src/Dets/options.irp.f | 8 + src/Dets/slater_rules.irp.f | 186 ++++++++++ src/Full_CI/H_apply.irp.f | 27 ++ src/Full_CI/NEEDED_MODULES | 2 +- src/Full_CI/full_ci.ezfio_config | 13 +- src/Full_CI/full_ci.irp.f | 22 +- src/Full_CI/options.irp.f | 8 + src/Generators_full/generators.irp.f | 15 + .../ASSUMPTIONS.rst | 0 src/Generators_restart/Makefile | 8 + src/Generators_restart/NEEDED_MODULES | 1 + src/Generators_restart/README.rst | 4 + src/Generators_restart/generators.irp.f | 59 +++ src/MOs/utils.irp.f | 41 +++ src/MP2/NEEDED_MODULES | 2 +- src/Molden/print_mo.irp.f | 6 +- src/MonoInts/pot_ao_ints.irp.f | 68 +++- src/MonoInts/pot_mo_ints.irp.f | 39 +- src/MonoInts/spread_dipole_ao.irp.f | 6 + src/NEEDED_MODULES | 2 +- src/Nuclei/nuclei.irp.f | 18 +- src/Perturbation/NEEDED_MODULES | 2 +- src/Perturbation/delta_rho_perturbation.irp.f | 72 ++++ src/Perturbation/dipole_moment.irp.f | 69 ++++ src/Perturbation/exc_max.irp.f | 4 + src/Perturbation/pert_single.irp.f | 49 +++ src/Perturbation/perturbation_template.f | 54 ++- src/Perturbation/selection.irp.f | 4 +- src/Primitive_basis/NEEDED_MODULES | 1 - src/Primitive_basis/functions.irp.f | 18 - src/Primitive_basis/prim.irp.f | 175 --------- .../primitive_basis.ezfio_config | 5 - src/Properties/Makefile | 8 + src/Properties/NEEDED_MODULES | 1 + src/Properties/README.rst | 4 + src/Properties/average.irp.f | 25 ++ src/Properties/delta_rho.irp.f | 224 ++++++++++++ src/Properties/need.irp.f | 289 +++++++++++++++ src/Properties/options.irp.f | 13 + src/Properties/properties.ezfio_config | 2 + src/Properties/properties.irp.f | 51 +++ src/Properties/routines_test.irp.f | 89 +++++ .../slater_rules_mono_electronic.irp.f | 339 ++++++++++++++++++ src/Selectors_full/selectors.irp.f | 12 + src/Selectors_no_sorted/ASSUMPTIONS.rst | 0 src/Selectors_no_sorted/Makefile | 8 + src/Selectors_no_sorted/NEEDED_MODULES | 1 + src/Selectors_no_sorted/README.rst | 4 + .../e_corr_selectors.irp.f | 79 ++++ src/Selectors_no_sorted/selectors.irp.f | 156 ++++++++ 80 files changed, 2714 insertions(+), 274 deletions(-) create mode 100644 src/CISD_selected/cisd.ezfio_config create mode 100644 src/CISD_selected/options.irp.f create mode 100644 src/Dets/diagonalize_CI_mono.irp.f create mode 100644 src/Dets/excitations_utils.irp.f create mode 100644 src/Dets/guess_doublet.irp.f create mode 100644 src/Dets/guess_singlet.irp.f create mode 100644 src/Dets/guess_triplet.irp.f rename src/{Primitive_basis => Generators_restart}/ASSUMPTIONS.rst (100%) create mode 100644 src/Generators_restart/Makefile create mode 100644 src/Generators_restart/NEEDED_MODULES create mode 100644 src/Generators_restart/README.rst create mode 100644 src/Generators_restart/generators.irp.f create mode 100644 src/Perturbation/delta_rho_perturbation.irp.f create mode 100644 src/Perturbation/dipole_moment.irp.f create mode 100644 src/Perturbation/exc_max.irp.f create mode 100644 src/Perturbation/pert_single.irp.f delete mode 100644 src/Primitive_basis/NEEDED_MODULES delete mode 100644 src/Primitive_basis/functions.irp.f delete mode 100644 src/Primitive_basis/prim.irp.f delete mode 100644 src/Primitive_basis/primitive_basis.ezfio_config create mode 100644 src/Properties/Makefile create mode 100644 src/Properties/NEEDED_MODULES create mode 100644 src/Properties/README.rst create mode 100644 src/Properties/average.irp.f create mode 100644 src/Properties/delta_rho.irp.f create mode 100644 src/Properties/need.irp.f create mode 100644 src/Properties/options.irp.f create mode 100644 src/Properties/properties.ezfio_config create mode 100644 src/Properties/properties.irp.f create mode 100644 src/Properties/routines_test.irp.f create mode 100644 src/Properties/slater_rules_mono_electronic.irp.f create mode 100644 src/Selectors_no_sorted/ASSUMPTIONS.rst create mode 100644 src/Selectors_no_sorted/Makefile create mode 100644 src/Selectors_no_sorted/NEEDED_MODULES create mode 100644 src/Selectors_no_sorted/README.rst create mode 100644 src/Selectors_no_sorted/e_corr_selectors.irp.f create mode 100644 src/Selectors_no_sorted/selectors.irp.f diff --git a/data/ezfio_defaults b/data/ezfio_defaults index 21163b55..9d20a5b7 100644 --- a/data/ezfio_defaults +++ b/data/ezfio_defaults @@ -23,13 +23,20 @@ determinants threshold_selectors 0.999 read_wf False s2_eig False + only_single_double_dm False full_ci n_det_max_fci 10000 + n_det_max_fci_property 50000 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 + do_pt2_end False + cassd n_det_max_cassd 10000 pt2_max 1.e-4 @@ -39,8 +46,14 @@ hartree_fock n_it_scf_max 200 thresh_scf 1.e-10 +cisd_selected + n_det_max_cisd 10000 + pt2_max 1.e-4 + cisd_sc2_selected n_det_max_cisd_sc2 10000 pt2_max 1.e-4 do_pt2_end True +properties + z_one_point 3.9 diff --git a/scripts/generate_h_apply.py b/scripts/generate_h_apply.py index 8bde1d11..35f9223f 100755 --- a/scripts/generate_h_apply.py +++ b/scripts/generate_h_apply.py @@ -13,6 +13,8 @@ initialization declarations decls_main keys_work +do_double_excitations +check_double_excitation copy_buffer finalization generate_psi_guess @@ -23,7 +25,8 @@ deinit_thread skip init_main filter_integrals -filter2h2p +filterhole +filterparticle """.split() class H_apply(object): @@ -116,12 +119,23 @@ class H_apply(object): buffer = buffer.replace('$'+key, value) return buffer - def set_filter_2h_2p(self): - self["filter2h2p"] = """ -! ! DIR$ FORCEINLINE - if(is_a_two_holes_two_particles(key))cycle + def unset_double_excitations(self): + self["do_double_excitations"] = ".False." + self["check_double_excitation"] = """ + check_double_excitation = .False. + """ + def set_filter_holes(self): + self["filterhole"] = """ + if(iand(ibset(0_bit_kind,j),hole(k,other_spin)).eq.0_bit_kind ) then + cycle + endif + """ + def set_filter_particl(self): + self["filterparticle"] = """ + if(iand(ibset(0_bit_kind,j_a),hole(k_a,other_spin)).eq.0_bit_kind ) then + cycle + endif """ - def set_perturbation(self,pert): if self.perturbation is not None: @@ -165,9 +179,14 @@ class H_apply(object): PROVIDE CI_electronic_energy psi_selectors_coef psi_selectors E_corr_per_selectors psi_det_sorted_bit """ 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,N_int) - """%(pert,) + if(check_double_excitation)then + 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,N_int) + else + call perturb_buffer_by_mono_%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,N_int) + endif + """%(pert,pert) self.data["finalization"] = """ """ self.data["copy_buffer"] = "" @@ -233,7 +252,7 @@ class H_apply(object): if (s2_eig) then call make_s2_eigenfunction endif -! SOFT_TOUCH psi_det psi_coef N_det +! SOFT_TOUCH psi_det psi_coef N_det selection_criterion_min = min(selection_criterion_min, maxval(select_max))*0.1d0 selection_criterion = selection_criterion_min call write_double(output_Dets,selection_criterion,'Selection criterion') diff --git a/src/BiInts/ao_bi_integrals.irp.f b/src/BiInts/ao_bi_integrals.irp.f index 1cbd3189..0cde449b 100644 --- a/src/BiInts/ao_bi_integrals.irp.f +++ b/src/BiInts/ao_bi_integrals.irp.f @@ -354,6 +354,7 @@ BEGIN_PROVIDER [ logical, ao_bielec_integrals_in_map ] real :: map_mb if (read_ao_integrals) then integer :: load_ao_integrals + print*,'Reading the AO integrals' if (load_ao_integrals(trim(ezfio_filename)//'/work/ao_integrals.bin') == 0) then write(output_BiInts,*) 'AO integrals provided' ao_bielec_integrals_in_map = .True. @@ -373,7 +374,7 @@ BEGIN_PROVIDER [ logical, ao_bielec_integrals_in_map ] PROVIDE progress_bar call omp_init_lock(lock) lmax = ao_num*(ao_num+1)/2 - write(output_BiInts,*) 'providing the AO integrals' + write(output_BiInts,*) 'Providing the AO integrals' call wall_time(wall_0) call wall_time(wall_1) call cpu_time(cpu_1) diff --git a/src/BiInts/map_integrals.irp.f b/src/BiInts/map_integrals.irp.f index 82d3b49b..4f707e3a 100644 --- a/src/BiInts/map_integrals.irp.f +++ b/src/BiInts/map_integrals.irp.f @@ -276,7 +276,6 @@ double precision function get_mo_bielec_integral(i,j,k,l,map) implicit none BEGIN_DOC ! Returns one integral in the MO basis - ! i(1)j(1) 1/r12 k(2)l(2) END_DOC integer, intent(in) :: i,j,k,l integer*8 :: idx @@ -293,7 +292,6 @@ double precision function mo_bielec_integral(i,j,k,l) implicit none BEGIN_DOC ! Returns one integral in the MO basis - ! i(1)j(1) 1/r12 k(2)l(2) END_DOC integer, intent(in) :: i,j,k,l double precision :: get_mo_bielec_integral @@ -308,7 +306,6 @@ subroutine get_mo_bielec_integrals(j,k,l,sze,out_val,map) BEGIN_DOC ! Returns multiple integrals in the MO basis, all ! i for j,k,l fixed. - ! i(1)j(1) 1/r12 k(2)l(2) END_DOC integer, intent(in) :: j,k,l, sze real(integral_kind), intent(out) :: out_val(sze) diff --git a/src/BiInts/mo_bi_integrals.irp.f b/src/BiInts/mo_bi_integrals.irp.f index eb9e650d..a3bcc00c 100644 --- a/src/BiInts/mo_bi_integrals.irp.f +++ b/src/BiInts/mo_bi_integrals.irp.f @@ -28,6 +28,7 @@ BEGIN_PROVIDER [ logical, mo_bielec_integrals_in_map ] mo_bielec_integrals_in_map = .True. if (read_mo_integrals) then integer :: load_mo_integrals + print*,'Reading the MO integrals' if (load_mo_integrals(trim(ezfio_filename)//'/work/mo_integrals.bin') == 0) then write(output_BiInts,*) 'MO integrals provided' return diff --git a/src/CAS_SD_selected/NEEDED_MODULES b/src/CAS_SD_selected/NEEDED_MODULES index c24125e8..7d7c112f 100644 --- a/src/CAS_SD_selected/NEEDED_MODULES +++ b/src/CAS_SD_selected/NEEDED_MODULES @@ -1 +1,2 @@ -AOs BiInts Bitmask Dets Electrons Ezfio_files Generators_CAS Hartree_Fock MOGuess MonoInts MOs Nuclei Output Perturbation Selectors_full Utils +AOs BiInts Bitmask Dets Electrons Ezfio_files Generators_CAS Hartree_Fock MOGuess MonoInts MOs Nuclei Output Perturbation Properties Selectors_full Utils + diff --git a/src/CID_SC2_selected/NEEDED_MODULES b/src/CID_SC2_selected/NEEDED_MODULES index e372b444..d7bca283 100644 --- a/src/CID_SC2_selected/NEEDED_MODULES +++ b/src/CID_SC2_selected/NEEDED_MODULES @@ -1,2 +1,2 @@ -AOs BiInts Bitmask CISD CISD_selected Dets Electrons Ezfio_files Hartree_Fock MOGuess MonoInts MOs Nuclei Output Perturbation Selectors_full SingleRefMethod Utils +AOs BiInts Bitmask CISD CISD_selected Dets Electrons Ezfio_files Hartree_Fock MOGuess MonoInts MOs Nuclei Output Perturbation Properties Selectors_full SingleRefMethod Utils diff --git a/src/CID_SC2_selected/cid_sc2_selection.irp.f b/src/CID_SC2_selected/cid_sc2_selection.irp.f index 3610a0fe..d2390229 100644 --- a/src/CID_SC2_selected/cid_sc2_selection.irp.f +++ b/src/CID_SC2_selected/cid_sc2_selection.irp.f @@ -68,6 +68,5 @@ program cisd_sc2_selected print*,'degree of excitation of such determinant : ',degree enddo - print*,'coucou' deallocate(pt2,norm_pert,H_pert_diag) end diff --git a/src/CID_selected/NEEDED_MODULES b/src/CID_selected/NEEDED_MODULES index e0486dce..56946f0e 100644 --- a/src/CID_selected/NEEDED_MODULES +++ b/src/CID_selected/NEEDED_MODULES @@ -1,2 +1,2 @@ -AOs BiInts Bitmask CISD Dets Electrons Ezfio_files Hartree_Fock MOGuess MonoInts MOs Nuclei Output Perturbation Selectors_full SingleRefMethod Utils +AOs BiInts Bitmask CISD Dets Electrons Ezfio_files Hartree_Fock MOGuess MonoInts MOs Nuclei Output Perturbation Properties Selectors_full SingleRefMethod Utils diff --git a/src/CISD_SC2_selected/H_apply.irp.f b/src/CISD_SC2_selected/H_apply.irp.f index 85ecae95..76e3d95e 100644 --- a/src/CISD_SC2_selected/H_apply.irp.f +++ b/src/CISD_SC2_selected/H_apply.irp.f @@ -3,6 +3,10 @@ BEGIN_SHELL [ /usr/bin/env python ] from generate_h_apply import * from perturbation import perturbations +s = H_apply("SC2_selected",SingleRef=True) +s.set_selection_pt2("epstein_nesbet_sc2_no_projected") +print s + s = H_apply("PT2",SingleRef=True) s.set_perturbation("epstein_nesbet_sc2_no_projected") print s diff --git a/src/CISD_SC2_selected/NEEDED_MODULES b/src/CISD_SC2_selected/NEEDED_MODULES index e372b444..d7bca283 100644 --- a/src/CISD_SC2_selected/NEEDED_MODULES +++ b/src/CISD_SC2_selected/NEEDED_MODULES @@ -1,2 +1,2 @@ -AOs BiInts Bitmask CISD CISD_selected Dets Electrons Ezfio_files Hartree_Fock MOGuess MonoInts MOs Nuclei Output Perturbation Selectors_full SingleRefMethod Utils +AOs BiInts Bitmask CISD CISD_selected Dets Electrons Ezfio_files Hartree_Fock MOGuess MonoInts MOs Nuclei Output Perturbation Properties Selectors_full SingleRefMethod Utils diff --git a/src/CISD_SC2_selected/cisd_sc2_selection.irp.f b/src/CISD_SC2_selected/cisd_sc2_selection.irp.f index 290693b1..fb7de4cf 100644 --- a/src/CISD_SC2_selected/cisd_sc2_selection.irp.f +++ b/src/CISD_SC2_selected/cisd_sc2_selection.irp.f @@ -37,8 +37,7 @@ program cisd_sc2_selected do while (N_det < n_det_max_cisd_sc2.and.maxval(abs(pt2(1:N_st))) > pt2_max) print*,'----' print*,'' - call H_apply_cisd_selection(perturbation,pt2, norm_pert, H_pert_diag, N_st) -! soft_touch det_connections + call H_apply_SC2_selected(pt2, norm_pert, H_pert_diag, N_st) call diagonalize_CI_SC2 print *, 'N_det = ', N_det do i = 1, N_st @@ -51,13 +50,14 @@ program cisd_sc2_selected endif E_old(i) = CI_SC2_energy(i) enddo -! print *, 'E corr = ', (E_old(1)) - HF_energy if(dabs(E_old(i) - CI_SC2_energy(i) ).le.1.d-12)then i_count += 1 selection_criterion_factor = selection_criterion_factor * 0.5d0 if(i_count > 5)then exit endif + else + i_count = 0 endif if (abort_all) then exit diff --git a/src/CISD_selected/NEEDED_MODULES b/src/CISD_selected/NEEDED_MODULES index e0486dce..56946f0e 100644 --- a/src/CISD_selected/NEEDED_MODULES +++ b/src/CISD_selected/NEEDED_MODULES @@ -1,2 +1,2 @@ -AOs BiInts Bitmask CISD Dets Electrons Ezfio_files Hartree_Fock MOGuess MonoInts MOs Nuclei Output Perturbation Selectors_full SingleRefMethod Utils +AOs BiInts Bitmask CISD Dets Electrons Ezfio_files Hartree_Fock MOGuess MonoInts MOs Nuclei Output Perturbation Properties Selectors_full SingleRefMethod Utils diff --git a/src/CISD_selected/cisd.ezfio_config b/src/CISD_selected/cisd.ezfio_config new file mode 100644 index 00000000..be3854b1 --- /dev/null +++ b/src/CISD_selected/cisd.ezfio_config @@ -0,0 +1,3 @@ +cisd_selected + n_det_max_cisd integer + pt2_max double precision diff --git a/src/CISD_selected/cisd_selection.irp.f b/src/CISD_selected/cisd_selection.irp.f index 5d0b8c3a..c74d0f3a 100644 --- a/src/CISD_selected/cisd_selection.irp.f +++ b/src/CISD_selected/cisd_selection.irp.f @@ -18,11 +18,13 @@ program cisd print *, 'E = ', CI_energy(i) enddo E_old = CI_energy - do while (maxval(abs(pt2(1:N_st))) > 1.d-4) + do while (maxval(abs(pt2(1:N_st))) > pt2_max.and.n_det < n_det_max_cisd) print*,'----' print*,'' call H_apply_cisd_selection(perturbation,pt2, norm_pert, H_pert_diag, N_st) call diagonalize_CI + psi_det = psi_det_sorted + psi_coef = psi_coef_sorted print*,'N_det = ',N_det do i = 1, N_st print*,'state ',i @@ -36,5 +38,9 @@ program cisd exit endif enddo + N_det = min(N_det,n_det_max_cisd) + touch N_det psi_det psi_coef + call diagonalize_CI deallocate(pt2,norm_pert,H_pert_diag) + call save_wavefunction end diff --git a/src/CISD_selected/options.irp.f b/src/CISD_selected/options.irp.f new file mode 100644 index 00000000..53aa0ef6 --- /dev/null +++ b/src/CISD_selected/options.irp.f @@ -0,0 +1,34 @@ + BEGIN_PROVIDER [ integer, n_det_max_cisd ] + implicit none + BEGIN_DOC + ! Get n_det_max_cisd from EZFIO file + END_DOC + logical :: has_n_det_max_cisd + PROVIDE ezfio_filename + call ezfio_has_cisd_selected_n_det_max_cisd(has_n_det_max_cisd) + if (has_n_det_max_cisd) then + call ezfio_get_cisd_selected_n_det_max_cisd(n_det_max_cisd) + else + n_det_max_cisd = 30000 + call ezfio_set_cisd_selected_n_det_max_cisd(n_det_max_cisd) + endif + print*,'n_det_max_cisd = ',n_det_max_cisd + END_PROVIDER + + BEGIN_PROVIDER [ double precision , pt2_max ] + implicit none + BEGIN_DOC + ! Get pt2_max from EZFIO file + END_DOC + logical :: has_pt2_max + PROVIDE ezfio_filename + call ezfio_has_cisd_selected_pt2_max(has_pt2_max) + if (has_pt2_max) then + call ezfio_get_cisd_selected_pt2_max(pt2_max) + else + pt2_max = 1.d-9 + call ezfio_set_cisd_selected_pt2_max(pt2_max) + endif + print*,'pt2_max = ',pt2_max + END_PROVIDER + diff --git a/src/DDCI_selected/NEEDED_MODULES b/src/DDCI_selected/NEEDED_MODULES index c24125e8..7d7c112f 100644 --- a/src/DDCI_selected/NEEDED_MODULES +++ b/src/DDCI_selected/NEEDED_MODULES @@ -1 +1,2 @@ -AOs BiInts Bitmask Dets Electrons Ezfio_files Generators_CAS Hartree_Fock MOGuess MonoInts MOs Nuclei Output Perturbation Selectors_full Utils +AOs BiInts Bitmask Dets Electrons Ezfio_files Generators_CAS Hartree_Fock MOGuess MonoInts MOs Nuclei Output Perturbation Properties Selectors_full Utils + diff --git a/src/Dets/H_apply_template.f b/src/Dets/H_apply_template.f index 557a6325..6402751f 100644 --- a/src/Dets/H_apply_template.f +++ b/src/Dets/H_apply_template.f @@ -26,7 +26,6 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_gene integer :: N_elec_in_key_hole_2(2),N_elec_in_key_part_2(2) double precision :: mo_bielec_integral - logical :: is_a_two_holes_two_particles integer, allocatable :: ia_ja_pairs(:,:,:) integer, allocatable :: ib_jb_pairs(:,:) double precision :: diag_H_mat_elem @@ -36,6 +35,11 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_gene ifirst=1 endif + logical :: check_double_excitation + check_double_excitation = .True. + + + $initialization $omp_parallel @@ -163,7 +167,6 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_gene k = ishft(j_b-1,-bit_kind_shift)+1 l = j_b-ishft(k-1,bit_kind_shift)-1 key(k,other_spin) = ibset(key(k,other_spin),l) - $filter2h2p key_idx += 1 do k=1,N_int keys_out(k,1,key_idx) = key(k,1) @@ -212,7 +215,6 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_gene k = ishft(j_b-1,-bit_kind_shift)+1 l = j_b-ishft(k-1,bit_kind_shift)-1 key(k,ispin) = ibset(key(k,ispin),l) - $filter2h2p key_idx += 1 do k=1,N_int keys_out(k,1,key_idx) = key(k,1) @@ -270,12 +272,17 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,i_generator,iproc $param 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) - logical :: is_a_two_holes_two_particles integer, allocatable :: ia_ja_pairs(:,:,:) logical, allocatable :: array_pairs(:,:) double precision :: diag_H_mat_elem integer(omp_lock_kind), save :: lck, ifirst=0 + + logical :: check_double_excitation + check_double_excitation = .True. + $check_double_excitation + + if (ifirst == 0) then ifirst=1 !$ call omp_init_lock(lck) @@ -333,11 +340,12 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,i_generator,iproc $param hole = key_in k = ishft(i_a-1,-bit_kind_shift)+1 j = i_a-ishft(k-1,bit_kind_shift)-1 + $filterhole hole(k,ispin) = ibclr(hole(k,ispin),j) k_a = ishft(j_a-1,-bit_kind_shift)+1 l_a = j_a-ishft(k_a-1,bit_kind_shift)-1 + $filterparticle hole(k_a,ispin) = ibset(hole(k_a,ispin),l_a) - $filter2h2p key_idx += 1 do k=1,N_int keys_out(k,1,key_idx) = hole(k,1) diff --git a/src/Dets/SC2.irp.f b/src/Dets/SC2.irp.f index 7bd73dea..8a6c10d7 100644 --- a/src/Dets/SC2.irp.f +++ b/src/Dets/SC2.irp.f @@ -141,7 +141,7 @@ subroutine CISD_SC2(dets_in,u_in,energies,dim_in,sze,N_st,Nint,convergence) accu += e_corr_array(j) endif enddo - case default + case default do j=1,N_double call get_excitation_degree(dets_in(1,1,i),doubles(1,1,j),degree,N_int) if (degree<=degree_exc(i)) then @@ -177,7 +177,7 @@ subroutine CISD_SC2(dets_in,u_in,energies,dim_in,sze,N_st,Nint,convergence) else call davidson_diag_hjj(dets_in,u_in,H_jj_dressed,energies,dim_in,sze,N_st,Nint,output_Dets) endif - + e_corr_double = 0.d0 inv_c0 = 1.d0/u_in(index_hf,1) do i = 1, N_double diff --git a/src/Dets/connected_to_ref.irp.f b/src/Dets/connected_to_ref.irp.f index 83d7e297..2245b94c 100644 --- a/src/Dets/connected_to_ref.irp.f +++ b/src/Dets/connected_to_ref.irp.f @@ -225,6 +225,114 @@ integer function connected_to_ref(key,keys,Nint,N_past_in,Ndet) end + + +integer function connected_to_ref_by_mono(key,keys,Nint,N_past_in,Ndet) + use bitmasks + implicit none + integer, intent(in) :: Nint, N_past_in, Ndet + integer(bit_kind), intent(in) :: keys(Nint,2,Ndet) + integer(bit_kind), intent(in) :: key(Nint,2) + + integer :: N_past + integer :: i, l + integer :: degree_x2 + logical :: det_is_not_or_may_be_in_ref, t + double precision :: hij_elec + + ! output : 0 : not connected + ! i : connected to determinant i of the past + ! -i : is the ith determinant of the refernce wf keys + + ASSERT (Nint > 0) + ASSERT (Nint == N_int) + + connected_to_ref_by_mono = 0 + N_past = max(1,N_past_in) + if (Nint == 1) then + + do i=N_past-1,1,-1 + degree_x2 = popcnt(xor( key(1,1), keys(1,1,i))) + & + popcnt(xor( key(1,2), keys(1,2,i))) + if (degree_x2 > 3.and. degree_x2 <5) then + cycle + else if (degree_x2 == 4)then + cycle + else if(degree_x2 == 2)then + connected_to_ref_by_mono = i + return + endif + enddo + + return + + + else if (Nint==2) then + + do i=N_past-1,1,-1 + degree_x2 = popcnt(xor( key(1,1), keys(1,1,i))) + & + popcnt(xor( key(1,2), keys(1,2,i))) + & + popcnt(xor( key(2,1), keys(2,1,i))) + & + popcnt(xor( key(2,2), keys(2,2,i))) + if (degree_x2 > 3.and. degree_x2 <5) then + cycle + else if (degree_x2 == 4)then + cycle + else if(degree_x2 == 2)then + connected_to_ref_by_mono = i + return + endif + enddo + + return + + else if (Nint==3) then + + do i=N_past-1,1,-1 + degree_x2 = popcnt(xor( key(1,1), keys(1,1,i))) + & + popcnt(xor( key(1,2), keys(1,2,i))) + & + popcnt(xor( key(2,1), keys(2,1,i))) + & + popcnt(xor( key(2,2), keys(2,2,i))) + & + popcnt(xor( key(3,1), keys(3,1,i))) + & + popcnt(xor( key(3,2), keys(3,2,i))) + if (degree_x2 > 3.and. degree_x2 <5) then + cycle + else if (degree_x2 == 4)then + cycle + else if(degree_x2 == 2)then + connected_to_ref_by_mono = i + return + endif + enddo + + return + + else + + do i=N_past-1,1,-1 + degree_x2 = popcnt(xor( key(1,1), keys(1,1,i))) + & + popcnt(xor( key(1,2), keys(1,2,i))) + !DEC$ LOOP COUNT MIN(3) + do l=2,Nint + degree_x2 = degree_x2 + popcnt(xor( key(l,1), keys(l,1,i))) +& + popcnt(xor( key(l,2), keys(l,2,i))) + enddo + if (degree_x2 > 3.and. degree_x2 <5) then + cycle + else if (degree_x2 == 4)then + cycle + else if(degree_x2 == 2)then + connected_to_ref_by_mono = i + return + endif + enddo + + endif + +end + + + logical function det_is_not_or_may_be_in_ref(key,Nint) use bitmasks implicit none diff --git a/src/Dets/density_matrix.irp.f b/src/Dets/density_matrix.irp.f index 8fe608ab..f72b337c 100644 --- a/src/Dets/density_matrix.irp.f +++ b/src/Dets/density_matrix.irp.f @@ -13,13 +13,96 @@ integer :: exc(0:2,2,2),n_occ_alpha double precision, allocatable :: tmp_a(:,:), tmp_b(:,:) - one_body_dm_mo_alpha = 0.d0 - one_body_dm_mo_beta = 0.d0 + if(only_single_double_dm)then + print*,'ONLY DOUBLE DM' + one_body_dm_mo_alpha = one_body_single_double_dm_mo_alpha + one_body_dm_mo_beta = one_body_single_double_dm_mo_beta + else + one_body_dm_mo_alpha = 0.d0 + one_body_dm_mo_beta = 0.d0 + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(j,k,l,m,occ,ck, cl, ckl,phase,h1,h2,p1,p2,s1,s2, degree,exc, & + !$OMP tmp_a, tmp_b, n_occ_alpha)& + !$OMP SHARED(psi_det,psi_coef,N_int,N_states,state_average_weight,elec_alpha_num,& + !$OMP elec_beta_num,one_body_dm_mo_alpha,one_body_dm_mo_beta,N_det,mo_tot_num_align,& + !$OMP mo_tot_num) + allocate(tmp_a(mo_tot_num_align,mo_tot_num), tmp_b(mo_tot_num_align,mo_tot_num) ) + tmp_a = 0.d0 + tmp_b = 0.d0 + !$OMP DO SCHEDULE(dynamic) + do k=1,N_det + call bitstring_to_list(psi_det(1,1,k), occ(1,1), n_occ_alpha, N_int) + call bitstring_to_list(psi_det(1,2,k), occ(1,2), n_occ_alpha, N_int) + do m=1,N_states + ck = psi_coef(k,m)*psi_coef(k,m) * state_average_weight(m) + do l=1,elec_alpha_num + j = occ(l,1) + tmp_a(j,j) += ck + enddo + do l=1,elec_beta_num + j = occ(l,2) + tmp_b(j,j) += ck + enddo + enddo + do l=1,k-1 + call get_excitation_degree(psi_det(1,1,k),psi_det(1,1,l),degree,N_int) + if (degree /= 1) then + cycle + endif + call get_mono_excitation(psi_det(1,1,k),psi_det(1,1,l),exc,phase,N_int) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + do m=1,N_states + ckl = psi_coef(k,m) * psi_coef(l,m) * phase * state_average_weight(m) + if (s1==1) then + tmp_a(h1,p1) += ckl + tmp_a(p1,h1) += ckl + else + tmp_b(h1,p1) += ckl + tmp_b(p1,h1) += ckl + endif + enddo + enddo + enddo + !$OMP END DO NOWAIT + !$OMP CRITICAL + one_body_dm_mo_alpha = one_body_dm_mo_alpha + tmp_a + !$OMP END CRITICAL + !$OMP CRITICAL + one_body_dm_mo_beta = one_body_dm_mo_beta + tmp_b + !$OMP END CRITICAL + deallocate(tmp_a,tmp_b) + !$OMP BARRIER + !$OMP END PARALLEL + + endif +END_PROVIDER + + BEGIN_PROVIDER [ double precision, one_body_single_double_dm_mo_alpha, (mo_tot_num_align,mo_tot_num) ] +&BEGIN_PROVIDER [ double precision, one_body_single_double_dm_mo_beta, (mo_tot_num_align,mo_tot_num) ] + implicit none + BEGIN_DOC + ! Alpha and beta one-body density matrix for each state + END_DOC + + integer :: j,k,l,m + integer :: occ(N_int*bit_kind_size,2) + double precision :: ck, cl, ckl + double precision :: phase + integer :: h1,h2,p1,p2,s1,s2, degree + integer :: exc(0:2,2,2),n_occ_alpha + double precision, allocatable :: tmp_a(:,:), tmp_b(:,:) + integer :: degree_respect_to_HF_k + integer :: degree_respect_to_HF_l + + PROVIDE elec_alpha_num elec_beta_num + + one_body_single_double_dm_mo_alpha = 0.d0 + one_body_single_double_dm_mo_beta = 0.d0 !$OMP PARALLEL DEFAULT(NONE) & !$OMP PRIVATE(j,k,l,m,occ,ck, cl, ckl,phase,h1,h2,p1,p2,s1,s2, degree,exc, & - !$OMP tmp_a, tmp_b, n_occ_alpha)& - !$OMP SHARED(psi_det,psi_coef,N_int,N_states,state_average_weight,elec_alpha_num,& - !$OMP elec_beta_num,one_body_dm_mo_alpha,one_body_dm_mo_beta,N_det,mo_tot_num_align,& + !$OMP tmp_a, tmp_b, n_occ_alpha,degree_respect_to_HF_k,degree_respect_to_HF_l)& + !$OMP SHARED(ref_bitmask,psi_det,psi_coef,N_int,N_states,state_average_weight,elec_alpha_num,& + !$OMP elec_beta_num,one_body_single_double_dm_mo_alpha,one_body_single_double_dm_mo_beta,N_det,mo_tot_num_align,& !$OMP mo_tot_num) allocate(tmp_a(mo_tot_num_align,mo_tot_num), tmp_b(mo_tot_num_align,mo_tot_num) ) tmp_a = 0.d0 @@ -28,18 +111,26 @@ do k=1,N_det call bitstring_to_list(psi_det(1,1,k), occ(1,1), n_occ_alpha, N_int) call bitstring_to_list(psi_det(1,2,k), occ(1,2), n_occ_alpha, N_int) + call get_excitation_degree(ref_bitmask,psi_det(1,1,k),degree_respect_to_HF_k,N_int) + do m=1,N_states ck = psi_coef(k,m)*psi_coef(k,m) * state_average_weight(m) - do l=1,elec_alpha_num - j = occ(l,1) - tmp_a(j,j) += ck - enddo - do l=1,elec_beta_num - j = occ(l,2) - tmp_b(j,j) += ck - enddo + call get_excitation_degree(ref_bitmask,psi_det(1,1,k),degree_respect_to_HF_l,N_int) + if(degree_respect_to_HF_l.le.0)then + do l=1,elec_alpha_num + j = occ(l,1) + tmp_a(j,j) += ck + enddo + do l=1,elec_beta_num + j = occ(l,2) + tmp_b(j,j) += ck + enddo + endif enddo do l=1,k-1 + call get_excitation_degree(ref_bitmask,psi_det(1,1,l),degree_respect_to_HF_l,N_int) + if(degree_respect_to_HF_k.ne.0)cycle + if(degree_respect_to_HF_l.eq.2.and.degree_respect_to_HF_k.ne.2)cycle call get_excitation_degree(psi_det(1,1,k),psi_det(1,1,l),degree,N_int) if (degree /= 1) then cycle @@ -56,14 +147,14 @@ tmp_b(p1,h1) += ckl endif enddo + enddo enddo - enddo !$OMP END DO NOWAIT !$OMP CRITICAL - one_body_dm_mo_alpha = one_body_dm_mo_alpha + tmp_a + one_body_single_double_dm_mo_alpha = one_body_single_double_dm_mo_alpha + tmp_a !$OMP END CRITICAL !$OMP CRITICAL - one_body_dm_mo_beta = one_body_dm_mo_beta + tmp_b + one_body_single_double_dm_mo_beta = one_body_single_double_dm_mo_beta + tmp_b !$OMP END CRITICAL deallocate(tmp_a,tmp_b) !$OMP BARRIER @@ -78,6 +169,14 @@ BEGIN_PROVIDER [ double precision, one_body_dm_mo, (mo_tot_num_align,mo_tot_num) one_body_dm_mo = one_body_dm_mo_alpha + one_body_dm_mo_beta END_PROVIDER +BEGIN_PROVIDER [ double precision, one_body_spin_density_mo, (mo_tot_num_align,mo_tot_num) ] + implicit none + BEGIN_DOC + ! rho(alpha) - rho(beta) + END_DOC + one_body_spin_density_mo = one_body_dm_mo_alpha - one_body_dm_mo_beta +END_PROVIDER + subroutine set_natural_mos implicit none BEGIN_DOC @@ -100,7 +199,7 @@ subroutine save_natural_mos ! Save natural orbitals, obtained by diagonalization of the one-body density matrix in the MO basis END_DOC call set_natural_mos - call save_mos + call save_mos end diff --git a/src/Dets/determinants.ezfio_config b/src/Dets/determinants.ezfio_config index 38ff88c9..0937502a 100644 --- a/src/Dets/determinants.ezfio_config +++ b/src/Dets/determinants.ezfio_config @@ -16,4 +16,5 @@ determinants read_wf logical expected_s2 double precision s2_eig logical + only_single_double_dm logical diff --git a/src/Dets/determinants.irp.f b/src/Dets/determinants.irp.f index c75028b6..89335932 100644 --- a/src/Dets/determinants.irp.f +++ b/src/Dets/determinants.irp.f @@ -29,6 +29,20 @@ BEGIN_PROVIDER [ integer, N_det ] ASSERT (N_det > 0) END_PROVIDER +BEGIN_PROVIDER [integer, max_degree_exc] + implicit none + integer :: i,degree + max_degree_exc = 0 + BEGIN_DOC + ! Maximum degree of excitation in the wf + END_DOC + do i = 1, N_det + call get_excitation_degree(HF_bitmask,psi_det(1,1,i),degree,N_int) + if(degree.gt.max_degree_exc)then + max_degree_exc= degree + endif + enddo +END_PROVIDER BEGIN_PROVIDER [ integer, psi_det_size ] implicit none @@ -229,7 +243,7 @@ BEGIN_PROVIDER [ double precision, psi_coef, (psi_det_size,N_states_diag) ] do i=1,N_states_diag psi_coef(i,i) = 1.d0 enddo - + if (read_wf) then call ezfio_has_determinants_psi_coef(exists) if (exists) then @@ -729,6 +743,16 @@ subroutine save_wavefunction call save_wavefunction_general(N_det,N_states,psi_det_sorted,size(psi_coef_sorted,1),psi_coef_sorted) end + +subroutine save_wavefunction_unsorted + implicit none + use bitmasks + BEGIN_DOC +! Save the wave function into the EZFIO file + END_DOC + call save_wavefunction_general(N_det,N_states,psi_det,size(psi_coef,1),psi_coef) +end + subroutine save_wavefunction_general(ndet,nstates,psidet,dim_psicoef,psicoef) implicit none BEGIN_DOC @@ -792,11 +816,23 @@ subroutine save_wavefunction_general(ndet,nstates,psidet,dim_psicoef,psicoef) progress_bar(1) = 7 progress_value = dble(progress_bar(1)) allocate (psi_coef_save(ndet,nstates)) + double precision :: accu_norm(nstates) + accu_norm = 0.d0 do k=1,nstates do i=1,ndet + accu_norm(k) = accu_norm(k) + psicoef(i,k) * psicoef(i,k) psi_coef_save(i,k) = psicoef(i,k) enddo enddo + do k = 1, nstates + accu_norm(k) = 1.d0/dsqrt(accu_norm(k)) + enddo + do k=1,nstates + do i=1,ndet + psi_coef_save(i,k) = psi_coef_save(i,k) * accu_norm(k) + enddo + enddo + call ezfio_set_determinants_psi_coef(psi_coef_save) call write_int(output_dets,ndet,'Saved determinants') call stop_progress diff --git a/src/Dets/diagonalize_CI.irp.f b/src/Dets/diagonalize_CI.irp.f index 288f28c3..55612920 100644 --- a/src/Dets/diagonalize_CI.irp.f +++ b/src/Dets/diagonalize_CI.irp.f @@ -69,10 +69,8 @@ END_PROVIDER i_state = 0 do j=1,N_det call get_s2_u0(psi_det,eigenvectors(1,j),N_det,N_det,s2) -! print *, 'j = ',j,s2, expected_s2 if(dabs(s2-expected_s2).le.0.3d0)then i_state += 1 -! print *, 'i_state = ',i_state do i=1,N_det CI_eigenvectors(i,i_state) = eigenvectors(i,j) enddo diff --git a/src/Dets/diagonalize_CI_SC2.irp.f b/src/Dets/diagonalize_CI_SC2.irp.f index 81931a62..86ba72b9 100644 --- a/src/Dets/diagonalize_CI_SC2.irp.f +++ b/src/Dets/diagonalize_CI_SC2.irp.f @@ -35,7 +35,8 @@ END_PROVIDER do i=1,N_det CI_SC2_eigenvectors(i,j) = psi_coef(i,j) enddo - CI_SC2_electronic_energy(j) = CI_electronic_energy(j) +! TODO : check comment +! CI_SC2_electronic_energy(j) = CI_electronic_energy(j) enddo call CISD_SC2(psi_det,CI_SC2_eigenvectors,CI_SC2_electronic_energy, & diff --git a/src/Dets/diagonalize_CI_mono.irp.f b/src/Dets/diagonalize_CI_mono.irp.f new file mode 100644 index 00000000..a3c5b103 --- /dev/null +++ b/src/Dets/diagonalize_CI_mono.irp.f @@ -0,0 +1,72 @@ + BEGIN_PROVIDER [ double precision, CI_electronic_energy_mono, (N_states_diag) ] +&BEGIN_PROVIDER [ double precision, CI_eigenvectors_mono, (N_det,N_states_diag) ] +&BEGIN_PROVIDER [ double precision, CI_eigenvectors_s2_mono, (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_mono(i,j) = psi_coef(i,j) + enddo + enddo + + if (diag_algorithm == "Davidson") then + + call davidson_diag(psi_det,CI_eigenvectors_mono,CI_electronic_energy, & + size(CI_eigenvectors_mono,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_all_dets,1),N_det)) + allocate (eigenvalues(N_det)) + call lapack_diag(eigenvalues,eigenvectors, & + H_matrix_all_dets,size(H_matrix_all_dets,1),N_det) + CI_electronic_energy_mono(:) = 0.d0 + do i=1,N_det + CI_eigenvectors_mono(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 + print*,'j = ',j + print*,'e = ',eigenvalues(j) + print*,'c = ',dabs(eigenvectors(1,j)) + if(dabs(eigenvectors(1,j)).gt.0.9d0)then + i_state += 1 + do i=1,N_det + CI_eigenvectors_mono(i,i_state) = eigenvectors(i,j) + enddo + CI_electronic_energy_mono(i_state) = eigenvalues(j) + CI_eigenvectors_s2_mono(i_state) = s2 + endif + endif + if (i_state.ge.N_states_diag) then + exit + endif + enddo + deallocate(eigenvectors,eigenvalues) + endif + +END_PROVIDER + +subroutine diagonalize_CI_mono + 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_mono(i,j) + enddo + enddo + SOFT_TOUCH psi_coef CI_electronic_energy_mono CI_eigenvectors_mono CI_eigenvectors_s2_mono +end diff --git a/src/Dets/excitations_utils.irp.f b/src/Dets/excitations_utils.irp.f new file mode 100644 index 00000000..46e38b08 --- /dev/null +++ b/src/Dets/excitations_utils.irp.f @@ -0,0 +1,16 @@ +subroutine apply_mono(i_hole,i_particle,ispin_excit,key_in,Nint) + implicit none + integer, intent(in) :: i_hole,i_particle,ispin_excit,Nint + integer(bit_kind), intent(inout) :: key_in(Nint,2) + integer :: k,j + use bitmasks + ! hole + k = ishft(i_hole-1,-bit_kind_shift)+1 + j = i_hole-ishft(k-1,bit_kind_shift)-1 + key_in(k,ispin_excit) = ibclr(key_in(k,ispin_excit),j) + + k = ishft(i_particle-1,-bit_kind_shift)+1 + j = i_particle-ishft(k-1,bit_kind_shift)-1 + key_in(k,ispin_excit) = ibset(key_in(k,ispin_excit),j) + +end diff --git a/src/Dets/guess_doublet.irp.f b/src/Dets/guess_doublet.irp.f new file mode 100644 index 00000000..a44697c1 --- /dev/null +++ b/src/Dets/guess_doublet.irp.f @@ -0,0 +1,79 @@ +program put_gess + use bitmasks + implicit none + integer :: i,j,N_det_tmp,N_states_tmp + integer :: list(N_int*bit_kind_size,2) + integer(bit_kind) :: string(N_int,2) + integer(bit_kind) :: psi_det_tmp(N_int,2,3) + double precision :: psi_coef_tmp(3,1) + + integer :: iorb,jorb,korb + print*,'which open shells ?' + read(5,*)iorb,jorb,korb + print*,iorb,jorb,korb + N_states= 1 + N_det= 3 + + + list = 0 + list(1,1) = 1 + list(1,2) = 1 + list(2,1) = 2 + list(2,2) = 2 + list(3,1) = iorb + list(4,1) = jorb + list(3,2) = korb + print*,'passed' + call list_to_bitstring( string(1,1), list(1,1), elec_alpha_num, N_int) + print*,'passed' + call list_to_bitstring( string(1,2), list(1,2), elec_beta_num, N_int) + print*,'passed' + call print_det(string,N_int) + do j = 1,2 + do i = 1, N_int + psi_det(i,j,1) = string(i,j) + enddo + enddo + psi_coef(1,1) = 1.d0/dsqrt(3.d0) + + print*,'passed 1' + list = 0 + list(1,1) = 1 + list(1,2) = 1 + list(2,1) = 2 + list(2,2) = 2 + list(3,1) = iorb + list(4,1) = korb + list(3,2) = jorb + call list_to_bitstring( string(1,1), list(1,1), elec_alpha_num, N_int) + call list_to_bitstring( string(1,2), list(1,2), elec_beta_num, N_int) + call print_det(string,N_int) + do j = 1,2 + do i = 1, N_int + psi_det(i,j,2) = string(i,j) + enddo + enddo + psi_coef(2,1) = 1.d0/dsqrt(3.d0) + + print*,'passed 2' + list = 0 + list(1,1) = 1 + list(1,2) = 1 + list(2,1) = 2 + list(2,2) = 2 + list(3,1) = korb + list(4,1) = jorb + list(3,2) = iorb + call list_to_bitstring( string(1,1), list(1,1), elec_alpha_num, N_int) + call list_to_bitstring( string(1,2), list(1,2), elec_beta_num, N_int) + call print_det(string,N_int) + do j = 1,2 + do i = 1, N_int + psi_det(i,j,3) = string(i,j) + enddo + enddo + psi_coef(3,1) = 1.d0/dsqrt(3.d0) + print*,'passed 3' + + call save_wavefunction +end diff --git a/src/Dets/guess_singlet.irp.f b/src/Dets/guess_singlet.irp.f new file mode 100644 index 00000000..50f8dc4e --- /dev/null +++ b/src/Dets/guess_singlet.irp.f @@ -0,0 +1,44 @@ +program put_gess + use bitmasks + implicit none + integer :: i,j,N_det_tmp,N_states_tmp + integer :: list(N_int*bit_kind_size,2) + integer(bit_kind) :: string(N_int,2) + integer(bit_kind) :: psi_det_tmp(N_int,2,2) + double precision :: psi_coef_tmp(2,1) + + integer :: iorb,jorb + print*,'which open shells ?' + read(5,*)iorb,jorb + N_states= 1 + N_det= 2 + + + list = 0 + list(1,1) = iorb + list(1,2) = jorb + call list_to_bitstring( string(1,1), list(1,1), elec_alpha_num, N_int) + call list_to_bitstring( string(1,2), list(1,2), elec_beta_num, N_int) + call print_det(string,N_int) + do j = 1,2 + do i = 1, N_int + psi_det(i,j,1) = string(i,j) + enddo + enddo + psi_coef(1,1) = 1.d0/dsqrt(2.d0) + + list = 0 + list(1,1) = jorb + list(1,2) = iorb + call list_to_bitstring( string(1,1), list(1,1), elec_alpha_num, N_int) + call list_to_bitstring( string(1,2), list(1,2), elec_beta_num, N_int) + call print_det(string,N_int) + do j = 1,2 + do i = 1, N_int + psi_det(i,j,2) = string(i,j) + enddo + enddo + psi_coef(2,1) = 1.d0/dsqrt(2.d0) + + call save_wavefunction +end diff --git a/src/Dets/guess_triplet.irp.f b/src/Dets/guess_triplet.irp.f new file mode 100644 index 00000000..77f88c3e --- /dev/null +++ b/src/Dets/guess_triplet.irp.f @@ -0,0 +1,48 @@ +program put_gess + use bitmasks + implicit none + integer :: i,j,N_det_tmp,N_states_tmp + integer :: list(N_int*bit_kind_size,2) + integer(bit_kind) :: string(N_int,2) + integer(bit_kind) :: psi_det_tmp(N_int,2,2) + double precision :: psi_coef_tmp(2,1) + + integer :: iorb,jorb + print*,'which open shells ?' + read(5,*)iorb,jorb + N_states= 1 + N_det= 2 + print*,'iorb = ',iorb + print*,'jorb = ',jorb + + + list = 0 + list(1,1) = iorb + list(1,2) = jorb + string = 0 + call list_to_bitstring( string(1,1), list(1,1), elec_alpha_num, N_int) + call list_to_bitstring( string(1,2), list(1,2), elec_beta_num, N_int) + call print_det(string,N_int) + do j = 1,2 + do i = 1, N_int + psi_det(i,j,1) = string(i,j) + enddo + enddo + psi_coef(1,1) = 1.d0/dsqrt(2.d0) + + list = 0 + list(1,1) = jorb + list(1,2) = iorb + string = 0 + call list_to_bitstring( string(1,1), list(1,1), elec_alpha_num, N_int) + call list_to_bitstring( string(1,2), list(1,2), elec_beta_num, N_int) + call print_det(string,N_int) + do j = 1,2 + do i = 1, N_int + psi_det(i,j,2) = string(i,j) + enddo + enddo + psi_coef(2,1) = -1.d0/dsqrt(2.d0) + + call save_wavefunction +end diff --git a/src/Dets/options.irp.f b/src/Dets/options.irp.f index 0ff351ae..dda5c04a 100644 --- a/src/Dets/options.irp.f +++ b/src/Dets/options.irp.f @@ -22,6 +22,14 @@ T.set_ezfio_name( "read_wf" ) T.set_output ( "output_dets" ) print T +T.set_type ( "logical" ) +T.set_name ( "only_single_double_dm" ) +T.set_doc ( "If true, The One body DM is calculated with ignoring the Double<->Doubles extra diag elements" ) +T.set_ezfio_name( "only_single_double_dm" ) +T.set_output ( "output_dets" ) +print T + + T.set_name ( "s2_eig" ) T.set_doc ( "Force the wave function to be an eigenfunction of S^2" ) T.set_ezfio_name( "s2_eig" ) diff --git a/src/Dets/slater_rules.irp.f b/src/Dets/slater_rules.irp.f index f909bb24..7d431879 100644 --- a/src/Dets/slater_rules.irp.f +++ b/src/Dets/slater_rules.irp.f @@ -488,6 +488,146 @@ end + +subroutine i_H_j_verbose(key_i,key_j,Nint,hij,hmono,hdouble) + use bitmasks + implicit none + BEGIN_DOC + ! Returns where i and j are determinants + END_DOC + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2) + double precision, intent(out) :: hij,hmono,hdouble + + integer :: exc(0:2,2,2) + integer :: degree + double precision :: get_mo_bielec_integral + integer :: m,n,p,q + integer :: i,j,k + integer :: occ(Nint*bit_kind_size,2) + double precision :: diag_H_mat_elem, phase,phase_2 + integer :: n_occ_alpha, n_occ_beta + logical :: has_mipi(Nint*bit_kind_size) + double precision :: mipi(Nint*bit_kind_size), miip(Nint*bit_kind_size) + PROVIDE mo_bielec_integrals_in_map mo_integrals_map + + ASSERT (Nint > 0) + ASSERT (Nint == N_int) + ASSERT (sum(popcnt(key_i(:,1))) == elec_alpha_num) + ASSERT (sum(popcnt(key_i(:,2))) == elec_beta_num) + ASSERT (sum(popcnt(key_j(:,1))) == elec_alpha_num) + ASSERT (sum(popcnt(key_j(:,2))) == elec_beta_num) + + hij = 0.d0 + hmono = 0.d0 + hdouble = 0.d0 + !DEC$ FORCEINLINE + call get_excitation_degree(key_i,key_j,degree,Nint) + select case (degree) + case (2) + call get_double_excitation(key_i,key_j,exc,phase,Nint) + if (exc(0,1,1) == 1) then + ! Mono alpha, mono beta + hij = phase*get_mo_bielec_integral( & + exc(1,1,1), & + exc(1,1,2), & + exc(1,2,1), & + exc(1,2,2) ,mo_integrals_map) + else if (exc(0,1,1) == 2) then + ! Double alpha + hij = phase*(get_mo_bielec_integral( & + exc(1,1,1), & + exc(2,1,1), & + exc(1,2,1), & + exc(2,2,1) ,mo_integrals_map) - & + get_mo_bielec_integral( & + exc(1,1,1), & + exc(2,1,1), & + exc(2,2,1), & + exc(1,2,1) ,mo_integrals_map) ) + else if (exc(0,1,2) == 2) then + ! Double beta + hij = phase*(get_mo_bielec_integral( & + exc(1,1,2), & + exc(2,1,2), & + exc(1,2,2), & + exc(2,2,2) ,mo_integrals_map) - & + get_mo_bielec_integral( & + exc(1,1,2), & + exc(2,1,2), & + exc(2,2,2), & + exc(1,2,2) ,mo_integrals_map) ) + endif + case (1) + call get_mono_excitation(key_i,key_j,exc,phase,Nint) + call bitstring_to_list(key_i(1,1), occ(1,1), n_occ_alpha, Nint) + call bitstring_to_list(key_i(1,2), occ(1,2), n_occ_beta, Nint) + has_mipi = .False. + if (exc(0,1,1) == 1) then + ! Mono alpha + m = exc(1,1,1) + p = exc(1,2,1) + do k = 1, elec_alpha_num + i = occ(k,1) + if (.not.has_mipi(i)) then + mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) + miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map) + has_mipi(i) = .True. + endif + enddo + do k = 1, elec_beta_num + i = occ(k,2) + if (.not.has_mipi(i)) then + mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) + has_mipi(i) = .True. + endif + enddo + + do k = 1, elec_alpha_num + hdouble = hdouble + mipi(occ(k,1)) - miip(occ(k,1)) + enddo + do k = 1, elec_beta_num + hdouble = hdouble + mipi(occ(k,2)) + enddo + + else + ! Mono beta + m = exc(1,1,2) + p = exc(1,2,2) + do k = 1, elec_beta_num + i = occ(k,2) + if (.not.has_mipi(i)) then + mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) + miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map) + has_mipi(i) = .True. + endif + enddo + do k = 1, elec_alpha_num + i = occ(k,1) + if (.not.has_mipi(i)) then + mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) + has_mipi(i) = .True. + endif + enddo + + do k = 1, elec_alpha_num + hdouble = hdouble + mipi(occ(k,1)) + enddo + do k = 1, elec_beta_num + hdouble = hdouble + mipi(occ(k,2)) - miip(occ(k,2)) + enddo + + endif + hmono = mo_mono_elec_integral(m,p) + hij = phase*(hdouble + hmono) + + case (0) + hij = diag_H_mat_elem(key_i,Nint) + end select +end + + + subroutine i_H_psi(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array) use bitmasks implicit none @@ -523,6 +663,52 @@ subroutine i_H_psi(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array) enddo end +subroutine i_H_psi_sec_ord(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array,idx_interaction,interactions) + use bitmasks + implicit none + integer, intent(in) :: Nint, Ndet,Ndet_max,Nstate + integer(bit_kind), intent(in) :: keys(Nint,2,Ndet) + integer(bit_kind), intent(in) :: key(Nint,2) + double precision, intent(in) :: coef(Ndet_max,Nstate) + double precision, intent(out) :: i_H_psi_array(Nstate) + double precision, intent(out) :: interactions(Ndet) + integer,intent(out) :: idx_interaction(0:Ndet) + + integer :: i, ii,j + double precision :: phase + integer :: exc(0:2,2,2) + double precision :: hij + integer :: idx(0:Ndet),n_interact + BEGIN_DOC + ! for the various Nstates + END_DOC + + ASSERT (Nint > 0) + ASSERT (N_int == Nint) + ASSERT (Nstate > 0) + ASSERT (Ndet > 0) + ASSERT (Ndet_max >= Ndet) + i_H_psi_array = 0.d0 + call filter_connected_i_H_psi0(keys,key,Nint,Ndet,idx) + n_interact = 0 + do ii=1,idx(0) + i = idx(ii) + !DEC$ FORCEINLINE + call i_H_j(keys(1,1,i),key,Nint,hij) + if(dabs(hij).ge.1.d-8)then + if(i.ne.1)then + n_interact += 1 + interactions(n_interact) = hij + idx_interaction(n_interact) = i + endif + endif + do j = 1, Nstate + i_H_psi_array(j) = i_H_psi_array(j) + coef(i,j)*hij + enddo + enddo + idx_interaction(0) = n_interact +end + subroutine i_H_psi_SC2(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array,idx_repeat) use bitmasks diff --git a/src/Full_CI/H_apply.irp.f b/src/Full_CI/H_apply.irp.f index 76c046ae..deabdf01 100644 --- a/src/Full_CI/H_apply.irp.f +++ b/src/Full_CI/H_apply.irp.f @@ -10,5 +10,32 @@ s = H_apply("FCI_PT2") s.set_perturbation("epstein_nesbet_2x2") print s + + +s = H_apply("FCI_mono") +s.set_selection_pt2("epstein_nesbet_2x2") +s.unset_double_excitations() +print s + + +s = H_apply("select_mono_delta_rho") +s.unset_double_excitations() +s.set_selection_pt2("delta_rho_one_point") +print s + +s = H_apply("select_mono_di_delta_rho") +s.set_selection_pt2("delta_rho_one_point") +print s + +s = H_apply("pt2_mono_delta_rho") +s.unset_double_excitations() +s.set_perturbation("delta_rho_one_point") +print s + +s = H_apply("pt2_mono_di_delta_rho") +s.set_perturbation("delta_rho_one_point") +print s + + END_SHELL diff --git a/src/Full_CI/NEEDED_MODULES b/src/Full_CI/NEEDED_MODULES index f6002738..f2bb9ba7 100644 --- a/src/Full_CI/NEEDED_MODULES +++ b/src/Full_CI/NEEDED_MODULES @@ -1,2 +1,2 @@ -AOs BiInts Bitmask Dets Electrons Ezfio_files Generators_full Hartree_Fock MOGuess MonoInts MOs Nuclei Output Perturbation 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/Full_CI/full_ci.ezfio_config b/src/Full_CI/full_ci.ezfio_config index 65e9c393..a7e8a435 100644 --- a/src/Full_CI/full_ci.ezfio_config +++ b/src/Full_CI/full_ci.ezfio_config @@ -1,7 +1,8 @@ full_ci - n_det_max_fci integer - pt2_max double precision - do_pt2_end logical - energy double precision - energy_pt2 double precision - var_pt2_ratio double precision + n_det_max_fci integer + n_det_max_fci_property integer + pt2_max double precision + do_pt2_end logical + energy double precision + energy_pt2 double precision + var_pt2_ratio double precision diff --git a/src/Full_CI/full_ci.irp.f b/src/Full_CI/full_ci.irp.f index 6435ee3c..594cc474 100644 --- a/src/Full_CI/full_ci.irp.f +++ b/src/Full_CI/full_ci.irp.f @@ -27,8 +27,19 @@ program full_ci print *, 'E+PT2 = ', CI_energy+pt2 print *, '-----' endif + double precision :: i_H_psi_array(N_states),diag_H_mat_elem,h,i_O1_psi_array(N_states) + if(read_wf)then + call i_H_psi(psi_det(1,1,N_det),psi_det,psi_coef,N_int,N_det,psi_det_size,N_states,i_H_psi_array) + h = diag_H_mat_elem(psi_det(1,1,N_det),N_int) + selection_criterion = dabs(psi_coef(N_det,1) * (i_H_psi_array(1) - h * psi_coef(N_det,1))) * 0.1d0 + soft_touch selection_criterion + endif + + integer :: n_det_before + print*,'Beginning the selection ...' do while (N_det < n_det_max_fci.and.maxval(abs(pt2(1:N_st))) > pt2_max) + n_det_before = N_det call H_apply_FCI(pt2, norm_pert, H_pert_diag, N_st) PROVIDE psi_coef @@ -43,6 +54,9 @@ program full_ci endif call diagonalize_CI call save_wavefunction + if(n_det_before == N_det)then + selection_criterion = selection_criterion * 0.5d0 + endif print *, 'N_det = ', N_det print *, 'N_states = ', N_states print *, 'PT2 = ', pt2 @@ -55,16 +69,15 @@ program full_ci endif enddo N_det = min(n_det_max_fci,N_det) + touch N_det psi_det psi_coef + call diagonalize_CI if(do_pt2_end)then + print*,'Last iteration only to compute the PT2' threshold_selectors = 1.d0 threshold_generators = 0.999d0 - touch N_det psi_det psi_coef - call diagonalize_CI call H_apply_FCI_PT2(pt2, norm_pert, H_pert_diag, N_st) print *, 'Final step' -!! call remove_small_contributions -!! call diagonalize_CI print *, 'N_det = ', N_det print *, 'N_states = ', N_states print *, 'PT2 = ', pt2 @@ -73,5 +86,6 @@ program full_ci print *, '-----' call ezfio_set_full_ci_energy_pt2(CI_energy+pt2) endif + call save_wavefunction deallocate(pt2,norm_pert) end diff --git a/src/Full_CI/options.irp.f b/src/Full_CI/options.irp.f index 76ae6023..e276492b 100644 --- a/src/Full_CI/options.irp.f +++ b/src/Full_CI/options.irp.f @@ -9,6 +9,14 @@ T.set_ezfio_name( "N_det_max_fci" ) T.set_output ( "output_full_ci" ) print T +T.set_type ( "integer" ) +T.set_name ( "N_det_max_fci_property" ) +T.set_doc ( "Max number of determinants in the wave function when you select for a given property" ) +T.set_ezfio_dir ( "full_ci" ) +T.set_ezfio_name( "N_det_max_fci_property" ) +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" ) diff --git a/src/Generators_full/generators.irp.f b/src/Generators_full/generators.irp.f index 56e82935..58fefef7 100644 --- a/src/Generators_full/generators.irp.f +++ b/src/Generators_full/generators.irp.f @@ -50,6 +50,21 @@ BEGIN_PROVIDER [ integer(bit_kind), psi_generators, (N_int,2,psi_det_size) ] END_PROVIDER +BEGIN_PROVIDER [integer, degree_max_generators] + implicit none + BEGIN_DOC +! Max degree of excitation (respect to HF) of the generators + END_DOC + integer :: i,degree + degree_max_generators = 0 + do i = 1, N_det_generators + call get_excitation_degree(HF_bitmask,psi_generators(1,1,i),degree,N_int) + if(degree .gt. degree_max_generators)then + degree_max_generators = degree + endif + enddo +END_PROVIDER + BEGIN_PROVIDER [ integer, size_select_max] implicit none BEGIN_DOC diff --git a/src/Primitive_basis/ASSUMPTIONS.rst b/src/Generators_restart/ASSUMPTIONS.rst similarity index 100% rename from src/Primitive_basis/ASSUMPTIONS.rst rename to src/Generators_restart/ASSUMPTIONS.rst diff --git a/src/Generators_restart/Makefile b/src/Generators_restart/Makefile new file mode 100644 index 00000000..b2ea1de1 --- /dev/null +++ b/src/Generators_restart/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_restart/NEEDED_MODULES b/src/Generators_restart/NEEDED_MODULES new file mode 100644 index 00000000..26097b8b --- /dev/null +++ b/src/Generators_restart/NEEDED_MODULES @@ -0,0 +1 @@ +AOs BiInts Bitmask Dets Electrons Ezfio_files MonoInts MOs Nuclei Output Utils diff --git a/src/Generators_restart/README.rst b/src/Generators_restart/README.rst new file mode 100644 index 00000000..e7ab7045 --- /dev/null +++ b/src/Generators_restart/README.rst @@ -0,0 +1,4 @@ +========================= +Generators_restart Module +========================= + diff --git a/src/Generators_restart/generators.irp.f b/src/Generators_restart/generators.irp.f new file mode 100644 index 00000000..19a6a9fd --- /dev/null +++ b/src/Generators_restart/generators.irp.f @@ -0,0 +1,59 @@ +use bitmasks + +BEGIN_PROVIDER [ integer, N_det_generators ] + implicit none + BEGIN_DOC + ! Read the wave function + END_DOC + integer :: i + integer, save :: ifirst = 0 + double precision :: norm + read_wf = .True. + if(ifirst == 0)then + N_det_generators = N_det + ifirst = 1 + endif + call write_int(output_dets,N_det_generators,'Number of generators') +END_PROVIDER + + + BEGIN_PROVIDER [ integer(bit_kind), psi_generators, (N_int,2,N_det_generators) ] +&BEGIN_PROVIDER [ double precision, psi_generators_coef, (N_det_generators,N_states) + implicit none + BEGIN_DOC + ! read wf + ! + END_DOC + integer :: i, k + integer, save :: ifirst = 0 + if(ifirst == 0)then + do i=1,N_det_generators + do k=1,N_int + psi_generators(k,1,i) = psi_det(k,1,i) + psi_generators(k,2,i) = psi_det(k,2,i) + enddo + do k = 1, N_states + psi_generators_coef(i,k) = psi_coef(i,k) + enddo + enddo + ifirst = 1 + endif + +END_PROVIDER + +BEGIN_PROVIDER [ integer, size_select_max] + implicit none + BEGIN_DOC + ! Size of the select_max array + END_DOC + size_select_max = 10000 +END_PROVIDER + +BEGIN_PROVIDER [ double precision, select_max, (size_select_max) ] + implicit none + BEGIN_DOC + ! Memo to skip useless selectors + END_DOC + select_max = huge(1.d0) +END_PROVIDER + diff --git a/src/MOs/utils.irp.f b/src/MOs/utils.irp.f index 24765180..21d04c27 100644 --- a/src/MOs/utils.irp.f +++ b/src/MOs/utils.irp.f @@ -139,3 +139,44 @@ subroutine mo_as_eigvectors_of_mo_matrix_sort_by_observable(matrix,observable,n, mo_label = label SOFT_TOUCH mo_coef mo_label end + + +subroutine mo_sort_by_observable(observable,label) + implicit none + character*(64), intent(in) :: label + double precision, intent(in) :: observable(mo_tot_num) + + double precision, allocatable :: mo_coef_new(:,:),value(:) + integer,allocatable :: iorder(:) + + allocate(mo_coef_new(ao_num_align,mo_tot_num),value(mo_tot_num),iorder(mo_tot_num)) + print*,'allocate !' + mo_coef_new = mo_coef + + + do i = 1, mo_tot_num + iorder(i) = i + value(i) = observable(i) + enddo + integer :: i,j,k,index + print*,'sort ....' + call dsort(value,iorder,mo_tot_num) + do i = 1, mo_tot_num + index = iorder(i) + do j = 1, mo_tot_num + mo_coef(j,i) = mo_coef_new(j,index) + enddo + enddo + + write (output_mos,'(A)'), 'MOs are now **'//trim(label)//'**' + write (output_mos,'(A)'), '' + + + deallocate(mo_coef_new,value) +! call write_time(output_mos) + + mo_label = label + SOFT_TOUCH mo_coef mo_label +end + + diff --git a/src/MP2/NEEDED_MODULES b/src/MP2/NEEDED_MODULES index 424391da..98fbcc58 100644 --- a/src/MP2/NEEDED_MODULES +++ b/src/MP2/NEEDED_MODULES @@ -1,2 +1,2 @@ -AOs BiInts Bitmask Dets Electrons Ezfio_files Hartree_Fock MOGuess MonoInts MOs Nuclei Output Perturbation Selectors_full SingleRefMethod Utils +AOs BiInts Bitmask Dets Electrons Ezfio_files Hartree_Fock MOGuess MonoInts MOs Nuclei Output Perturbation Properties Selectors_full SingleRefMethod Utils diff --git a/src/Molden/print_mo.irp.f b/src/Molden/print_mo.irp.f index 30b14d92..9cec8fbd 100644 --- a/src/Molden/print_mo.irp.f +++ b/src/Molden/print_mo.irp.f @@ -91,7 +91,11 @@ subroutine write_Ao_basis(i_unit_output) ! write(i_unit_output,*),j,i_shell,i_ao!trim(character_shell) do k = 1, ao_prim_num(i_ao) i_prim +=1 - write(i_unit_output,'(4X,I3,3X,A1,6X,I2,6X,F16.7,2X,F16.12)')i_shell,character_shell,i_prim,ao_expo_unsorted(i_ao,k),ao_coef_unnormalized(i_ao,k) + if(i_prim.lt.100)then + write(i_unit_output,'(4X,I3,3X,A1,6X,I2,6X,F16.7,2X,F16.12)')i_shell,character_shell,i_prim,ao_expo_unsorted(i_ao,k),ao_coef_unnormalized(i_ao,k) + else + write(i_unit_output,'(4X,I3,3X,A1,5X,I3,6X,F16.7,2X,F16.12)')i_shell,character_shell,i_prim,ao_expo_unsorted(i_ao,k),ao_coef_unnormalized(i_ao,k) + endif enddo write(i_unit_output,*)'' enddo diff --git a/src/MonoInts/pot_ao_ints.irp.f b/src/MonoInts/pot_ao_ints.irp.f index d5d646ac..f430ace9 100644 --- a/src/MonoInts/pot_ao_ints.irp.f +++ b/src/MonoInts/pot_ao_ints.irp.f @@ -64,6 +64,73 @@ END_PROVIDER + BEGIN_PROVIDER [ double precision, ao_nucl_elec_integral_per_atom, (ao_num_align,ao_num,nucl_num)] + BEGIN_DOC +! ao_nucl_elec_integral_per_atom(i,j,k) = - +! where Rk is the geometry of the kth atom + END_DOC + implicit none + double precision :: alpha, beta, gama, delta + integer :: i_c,num_A,num_B + double precision :: A_center(3),B_center(3),C_center(3) + integer :: power_A(3),power_B(3) + integer :: i,j,k,l,n_pt_in,m + double precision ::overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult + integer :: nucl_numC + ! Important for OpenMP + + ao_nucl_elec_integral_per_atom = 0.d0 + + + do k = 1, nucl_num + C_center(1) = nucl_coord(k,1) + C_center(2) = nucl_coord(k,2) + C_center(3) = nucl_coord(k,3) + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i,j,l,m,alpha,beta,A_center,B_center,power_A,power_B, & + !$OMP num_A,num_B,c,n_pt_in) & + !$OMP SHARED (k,ao_num,ao_prim_num,ao_expo_transp,ao_power,ao_nucl,nucl_coord,ao_coef_transp, & + !$OMP n_pt_max_integrals,ao_nucl_elec_integral_per_atom,nucl_num,C_center) + n_pt_in = n_pt_max_integrals + !$OMP DO SCHEDULE (guided) + + double precision :: c + do j = 1, ao_num + power_A(1)= ao_power(j,1) + power_A(2)= ao_power(j,2) + power_A(3)= ao_power(j,3) + num_A = ao_nucl(j) + A_center(1) = nucl_coord(num_A,1) + A_center(2) = nucl_coord(num_A,2) + A_center(3) = nucl_coord(num_A,3) + do i = 1, ao_num + power_B(1)= ao_power(i,1) + power_B(2)= ao_power(i,2) + power_B(3)= ao_power(i,3) + num_B = ao_nucl(i) + B_center(1) = nucl_coord(num_B,1) + B_center(2) = nucl_coord(num_B,2) + B_center(3) = nucl_coord(num_B,3) + c = 0.d0 + do l=1,ao_prim_num(j) + alpha = ao_expo_transp(l,j) + do m=1,ao_prim_num(i) + beta = ao_expo_transp(m,i) + c = c + NAI_pol_mult(A_center,B_center,power_A,power_B,alpha,beta,C_center,n_pt_in) & + * ao_coef_transp(l,j)*ao_coef_transp(m,i) + enddo + enddo + ao_nucl_elec_integral_per_atom(i,j,k) = -c + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + enddo + +END_PROVIDER + + double precision function NAI_pol_mult(A_center,B_center,power_A,power_B,alpha,beta,C_center,n_pt_in) ! function that calculate the folowing integral : @@ -290,7 +357,6 @@ recursive subroutine I_x1_pol_mult_mono_elec(a,c,R1x,R1xp,R2x,d,nd,n_pt_in) ! print*,'nd_in = ',nd if( (a==0) .and. (c==0))then -! print*,'coucou !' nd = 0 d(0) = 1.d0 return diff --git a/src/MonoInts/pot_mo_ints.irp.f b/src/MonoInts/pot_mo_ints.irp.f index 6393e57e..50019abb 100644 --- a/src/MonoInts/pot_mo_ints.irp.f +++ b/src/MonoInts/pot_mo_ints.irp.f @@ -2,6 +2,9 @@ BEGIN_PROVIDER [double precision, mo_nucl_elec_integral, (mo_tot_num_align,mo_to implicit none integer :: i1,j1,i,j double precision :: c_i1,c_j1 + BEGIN_DOC +! interaction nuclear electron on the MO basis + END_DOC mo_nucl_elec_integral = 0.d0 !$OMP PARALLEL DO DEFAULT(none) & @@ -11,9 +14,9 @@ BEGIN_PROVIDER [double precision, mo_nucl_elec_integral, (mo_tot_num_align,mo_to do i = 1, mo_tot_num do j = 1, mo_tot_num do i1 = 1,ao_num - c_i1 = mo_coef(i1,i) ! + c_i1 = mo_coef(i1,i) do j1 = 1,ao_num - c_j1 = c_i1*mo_coef(j1,j) ! + c_j1 = c_i1*mo_coef(j1,j) mo_nucl_elec_integral(j,i) = mo_nucl_elec_integral(j,i) + & c_j1 * ao_nucl_elec_integral(j1,i1) enddo @@ -23,3 +26,35 @@ BEGIN_PROVIDER [double precision, mo_nucl_elec_integral, (mo_tot_num_align,mo_to !$OMP END PARALLEL DO END_PROVIDER + +BEGIN_PROVIDER [double precision, mo_nucl_elec_integral_per_atom, (mo_tot_num_align,mo_tot_num,nucl_num)] + implicit none + integer :: i1,j1,i,j,k + double precision :: c_i1,c_j1 + BEGIN_DOC +! mo_nucl_elec_integral_per_atom(i,j,k) = - +! where Rk is the geometry of the kth atom + END_DOC + + mo_nucl_elec_integral_per_atom = 0.d0 + do k = 1, nucl_num + !$OMP PARALLEL DO DEFAULT(none) & + !$OMP PRIVATE(i,j,i1,j1,c_j1,c_i1) & + !$OMP SHARED(mo_tot_num,ao_num,mo_coef, & + !$OMP mo_nucl_elec_integral_per_atom, ao_nucl_elec_integral_per_atom,k) + do i = 1, mo_tot_num + do j = 1, mo_tot_num + do i1 = 1,ao_num + c_i1 = mo_coef(i1,i) + do j1 = 1,ao_num + c_j1 = c_i1*mo_coef(j1,j) + mo_nucl_elec_integral_per_atom(j,i,k) = mo_nucl_elec_integral_per_atom(j,i,k) + & + c_j1 * ao_nucl_elec_integral_per_atom(j1,i1,k) + enddo + enddo + enddo + enddo + !$OMP END PARALLEL DO + enddo +END_PROVIDER + diff --git a/src/MonoInts/spread_dipole_ao.irp.f b/src/MonoInts/spread_dipole_ao.irp.f index 2628c6a0..c0d7c88e 100644 --- a/src/MonoInts/spread_dipole_ao.irp.f +++ b/src/MonoInts/spread_dipole_ao.irp.f @@ -339,6 +339,12 @@ end dist = (A_center - B_center)*(A_center - B_center) P_center = (alpha * A_center + beta * B_center) * p_inv factor = dexp(-rho * dist) + if(power_B == 0 .and. power_A ==0)then + double precision :: F_integral + overlap_x = P_center * F_integral(0,p) * factor + dx = 0.d0 + return + endif double precision :: pouet_timy pouet_timy = dsqrt(lower_exp_val/p) diff --git a/src/NEEDED_MODULES b/src/NEEDED_MODULES index b5489e54..830bb3d8 100644 --- a/src/NEEDED_MODULES +++ b/src/NEEDED_MODULES @@ -1 +1 @@ -AOs BiInts Bitmask Dets Electrons Ezfio_files Full_CI Generators_full Hartree_Fock MOGuess MonoInts MOs Nuclei Output Selectors_full Utils Molden FCIdump Generators_CAS CAS_SD_selected DDCI_selected +AOs BiInts Bitmask Dets Electrons Ezfio_files Full_CI Generators_full Hartree_Fock MOGuess MonoInts MOs Nuclei Output Selectors_full Utils Molden FCIdump Generators_CAS CAS_SD_selected diff --git a/src/Nuclei/nuclei.irp.f b/src/Nuclei/nuclei.irp.f index 9fdc9c4d..ec1fb7d4 100644 --- a/src/Nuclei/nuclei.irp.f +++ b/src/Nuclei/nuclei.irp.f @@ -168,7 +168,23 @@ END_PROVIDER END_PROVIDER -BEGIN_PROVIDER [ double precision, nuclear_repulsion ] +BEGIN_PROVIDER [ double precision, positive_charge_barycentre,(3)] + implicit none + BEGIN_DOC + ! Centroid of the positive charges + END_DOC + integer :: l + positive_charge_barycentre(1) = 0.d0 + positive_charge_barycentre(2) = 0.d0 + positive_charge_barycentre(3) = 0.d0 + do l = 1, nucl_num + positive_charge_barycentre(1) += nucl_charge(l) * nucl_coord(l,1) + positive_charge_barycentre(2) += nucl_charge(l) * nucl_coord(l,2) + positive_charge_barycentre(3) += nucl_charge(l) * nucl_coord(l,3) + enddo +END_PROVIDER + + BEGIN_PROVIDER [ double precision, nuclear_repulsion ] implicit none BEGIN_DOC ! Nuclear repulsion energy diff --git a/src/Perturbation/NEEDED_MODULES b/src/Perturbation/NEEDED_MODULES index 3f5f6001..4f7ab529 100644 --- a/src/Perturbation/NEEDED_MODULES +++ b/src/Perturbation/NEEDED_MODULES @@ -1,2 +1,2 @@ -AOs BiInts Bitmask Dets Electrons Ezfio_files Hartree_Fock MOGuess MonoInts MOs Nuclei Output Utils +AOs BiInts Bitmask Dets Electrons Ezfio_files Hartree_Fock MOGuess MonoInts MOs Nuclei Output Properties Utils diff --git a/src/Perturbation/delta_rho_perturbation.irp.f b/src/Perturbation/delta_rho_perturbation.irp.f new file mode 100644 index 00000000..d83eb9a8 --- /dev/null +++ b/src/Perturbation/delta_rho_perturbation.irp.f @@ -0,0 +1,72 @@ +subroutine pt2_delta_rho_one_point(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_st) + use bitmasks + implicit none + integer, intent(in) :: Nint,ndet,n_st + integer(bit_kind), intent(in) :: det_pert(Nint,2) + double precision , intent(out) :: c_pert(n_st),e_2_pert(n_st),H_pert_diag(N_st) + double precision :: i_O1_psi_array(N_st) + double precision :: i_H_psi_array(N_st) + + BEGIN_DOC + ! compute the perturbatibe contribution to the Integrated Spin density at z = z_one point of one determinant + ! + ! for the various n_st states, at various level of theory. + ! + ! c_pert(i) = /( - ) + ! + ! e_2_pert(i) = c_pert(i) * + ! + ! H_pert_diag(i) = c_pert(i)^2 * + ! + ! To get the contribution of the first order : + ! + ! = sum(over i) e_2_pert(i) + ! + ! To get the contribution of the diagonal elements of the second order : + ! + ! [ + + sum(over i) H_pert_diag(i) ] / [1. + sum(over i) c_pert(i) **2] + ! + END_DOC + + integer :: i,j + double precision :: diag_H_mat_elem,diag_o1_mat_elem_alpha_beta + integer :: exc(0:2,2,2) + integer :: degree + double precision :: phase,delta_e,h,oii,diag_o1_mat_elem + integer :: h1,h2,p1,p2,s1,s2 + ASSERT (Nint == N_int) + ASSERT (Nint > 0) + +! call get_excitation_degree(HF_bitmask,det_pert,degree,N_int) +! if(degree.gt.degree_max_generators+1)then +! H_pert_diag = 0.d0 +! e_2_pert = 0.d0 +! c_pert = 0.d0 +! return +! endif + call i_O1_psi_alpha_beta(mo_integrated_delta_rho_one_point,det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_O1_psi_array) + + 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) + oii = diag_O1_mat_elem_alpha_beta(mo_integrated_delta_rho_one_point,det_pert,N_int) + + + do i =1,N_st + if(CI_electronic_energy(i)>h.and.CI_electronic_energy(i).ne.0.d0)then + c_pert(i) = -1.d0 + e_2_pert(i) = selection_criterion*selection_criterion_factor*2.d0 + else 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_O1_psi_array(i)+i_O1_psi_array(i) ) + c_pert(i) * c_pert(i) * oii + H_pert_diag(i) = c_pert(i) * (i_O1_psi_array(i)+i_O1_psi_array(i) ) + else + c_pert(i) = -1.d0 + e_2_pert(i) = -dabs(i_H_psi_array(i)) + H_pert_diag(i) = c_pert(i) * i_O1_psi_array(i) + endif + enddo + + +end + diff --git a/src/Perturbation/dipole_moment.irp.f b/src/Perturbation/dipole_moment.irp.f new file mode 100644 index 00000000..ca09c31c --- /dev/null +++ b/src/Perturbation/dipole_moment.irp.f @@ -0,0 +1,69 @@ +subroutine pt2_dipole_moment_z(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_st) + use bitmasks + implicit none + integer, intent(in) :: Nint,ndet,n_st + integer(bit_kind), intent(in) :: det_pert(Nint,2) + double precision , intent(out) :: c_pert(n_st),e_2_pert(n_st),H_pert_diag(N_st) + double precision :: i_O1_psi_array(N_st) + double precision :: i_H_psi_array(N_st) + + BEGIN_DOC + ! compute the perturbatibe contribution to the dipole moment of one determinant + ! + ! for the various n_st states, at various level of theory. + ! + ! c_pert(i) = /( - ) + ! + ! e_2_pert(i) = c_pert(i) * + ! + ! H_pert_diag(i) = c_pert(i)^2 * + ! + ! To get the contribution of the first order : + ! + ! = sum(over i) e_2_pert(i) + ! + ! To get the contribution of the diagonal elements of the second order : + ! + ! [ + + sum(over i) H_pert_diag(i) ] / [1. + sum(over i) c_pert(i) **2] + ! + END_DOC + + integer :: i,j + double precision :: diag_H_mat_elem + integer :: exc(0:2,2,2) + integer :: degree + double precision :: phase,delta_e,h,oii,diag_o1_mat_elem + integer :: h1,h2,p1,p2,s1,s2 + ASSERT (Nint == N_int) + ASSERT (Nint > 0) + +! call get_excitation_degree(HF_bitmask,det_pert,degree,N_int) +! if(degree.gt.degree_max_generators+1)then +! H_pert_diag = 0.d0 +! e_2_pert = 0.d0 +! c_pert = 0.d0 +! return +! endif + + call i_O1_psi(mo_dipole_z,det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_O1_psi_array) + 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) + oii = diag_O1_mat_elem(mo_dipole_z,det_pert,N_int) + + + do i =1,N_st + if(CI_electronic_energy(i)>h.and.CI_electronic_energy(i).ne.0.d0)then + c_pert(i) = -1.d0 + e_2_pert(i) = selection_criterion*selection_criterion_factor*2.d0 + else 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_O1_psi_array(i)+i_O1_psi_array(i) ) + H_pert_diag(i) = e_2_pert(i) + c_pert(i) * c_pert(i) * oii + else + c_pert(i) = -1.d0 + e_2_pert(i) = -dabs(i_H_psi_array(i)) + H_pert_diag(i) = c_pert(i) * i_O1_psi_array(i) + endif + enddo +end + diff --git a/src/Perturbation/exc_max.irp.f b/src/Perturbation/exc_max.irp.f new file mode 100644 index 00000000..15e86728 --- /dev/null +++ b/src/Perturbation/exc_max.irp.f @@ -0,0 +1,4 @@ +BEGIN_PROVIDER [integer, max_exc_pert] + implicit none + max_exc_pert = 0 +END_PROVIDER diff --git a/src/Perturbation/pert_single.irp.f b/src/Perturbation/pert_single.irp.f new file mode 100644 index 00000000..d04ca7ca --- /dev/null +++ b/src/Perturbation/pert_single.irp.f @@ -0,0 +1,49 @@ +subroutine pt2_h_core(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st) + use bitmasks + implicit none + integer, intent(in) :: Nint,ndet,N_st + integer(bit_kind), intent(in) :: det_pert(Nint,2) + double precision , intent(out) :: c_pert(N_st),e_2_pert(N_st),H_pert_diag(N_st) + double precision :: i_H_psi_array(N_st) + + BEGIN_DOC + ! compute the standard Epstein-Nesbet perturbative first order coefficient and second order energetic contribution + ! + ! for the various N_st states. + ! + ! c_pert(i) = /( E(i) - ) + ! + ! e_2_pert(i) = ^2/( E(i) - ) + ! + END_DOC + + integer :: i,j + double precision :: diag_H_mat_elem, h + + ASSERT (Nint == N_int) + ASSERT (Nint > 0) + + integer :: exc(0:2,2,2) + integer :: degree + double precision :: phase + call get_excitation(ref_bitmask,det_pert,exc,degree,phase,N_int) + h = diag_H_mat_elem(det_pert,N_int) + print*,'delta E = ',h-ref_bitmask_energy + if(h1)then + c_pert = 0.d0 + e_2_pert = 0.d0 + H_pert_diag = 0.d0 + return + endif + integer :: h1,p1,h2,p2,s1,s2 + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + c_pert = phase * mo_mono_elec_integral(h1,p1) + e_2_pert = -dabs(mo_mono_elec_integral(h1,p1)+1.d0) + +end diff --git a/src/Perturbation/perturbation_template.f b/src/Perturbation/perturbation_template.f index a450edff..3627dd17 100644 --- a/src/Perturbation/perturbation_template.f +++ b/src/Perturbation/perturbation_template.f @@ -34,7 +34,59 @@ subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,c if (is_in_wavefunction(buffer(1,1,i),Nint,N_det)) then cycle endif - + + integer :: degree + call get_excitation_degree(HF_bitmask,buffer(1,1,i),degree,N_int) + call pt2_$PERT(buffer(1,1,i), & + c_pert,e_2_pert,H_pert_diag,Nint,N_det_selectors,n_st) + + do k = 1,N_st + e_2_pert_buffer(k,i) = e_2_pert(k) + 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) += H_pert_diag(k) + enddo + + enddo + +end + + +subroutine perturb_buffer_by_mono_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert,sum_norm_pert,sum_H_pert_diag,N_st,Nint) + implicit none + BEGIN_DOC + ! Applly pertubration ``$PERT`` to the buffer of determinants generated in the H_apply + ! routine. + END_DOC + + integer, intent(in) :: Nint, N_st, buffer_size, i_generator + 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(N_st) + integer :: i,k, c_ref + integer, external :: connected_to_ref_by_mono + logical, external :: is_in_wavefunction + + ASSERT (Nint > 0) + ASSERT (Nint == N_int) + ASSERT (buffer_size >= 0) + ASSERT (minval(sum_norm_pert) >= 0.d0) + ASSERT (N_st > 0) + do i = 1,buffer_size + + c_ref = connected_to_ref_by_mono(buffer(1,1,i),psi_generators,Nint,i_generator,N_det) + + if (c_ref /= 0) then + cycle + endif + + if (is_in_wavefunction(buffer(1,1,i),Nint,N_det)) then + cycle + endif + + integer :: degree call pt2_$PERT(buffer(1,1,i), & c_pert,e_2_pert,H_pert_diag,Nint,N_det_selectors,n_st) diff --git a/src/Perturbation/selection.irp.f b/src/Perturbation/selection.irp.f index 15d07ff0..4230293a 100644 --- a/src/Perturbation/selection.irp.f +++ b/src/Perturbation/selection.irp.f @@ -36,7 +36,7 @@ subroutine fill_H_apply_buffer_selection(n_selected,det_buffer,e_2_pert_buffer,c is_selected = .False. do j=1,N_st - s = -e_2_pert_buffer(j,i) + s = dabs(e_2_pert_buffer(j,i)) is_selected = s > selection_criterion*selection_criterion_factor .or. is_selected select_max_out = max(select_max_out,s) enddo @@ -78,7 +78,7 @@ end BEGIN_DOC ! Threshold to select determinants. Set by selection routines. END_DOC - selection_criterion = .1d0 + selection_criterion = 0.1d0 selection_criterion_factor = 0.01d0 selection_criterion_min = selection_criterion diff --git a/src/Primitive_basis/NEEDED_MODULES b/src/Primitive_basis/NEEDED_MODULES deleted file mode 100644 index 190f8c6e..00000000 --- a/src/Primitive_basis/NEEDED_MODULES +++ /dev/null @@ -1 +0,0 @@ -AOs Electrons Ezfio_files MOs Nuclei Output Utils diff --git a/src/Primitive_basis/functions.irp.f b/src/Primitive_basis/functions.irp.f deleted file mode 100644 index f58fe91e..00000000 --- a/src/Primitive_basis/functions.irp.f +++ /dev/null @@ -1,18 +0,0 @@ -double precision function mo_r(i,r) - implicit none - integer :: i - double precision :: r(3) - integer :: j - double precision :: x_center,y_center,z_center,r_center - mo_r = 0.d0 - do j = 1, prim_num - if(dabs(prim_mo_coef(j,i)).lt.1.d-10)cycle - x_center = r(1) - nucl_coord(prim_nucl(j),1) - y_center = r(2) - nucl_coord(prim_nucl(j),2) - z_center = r(3) - nucl_coord(prim_nucl(j),3) - r_center = x_center*x_center + y_center*y_center + z_center*z_center - mo_r = mo_r + prim_mo_coef(j,i) * (x_center ** (prim_power(j,1)) * y_center ** (prim_power(j,2)) * z_center ** (prim_power(j,3)) & - * dexp(-prim_expo(j) * r_center)) - enddo - -end diff --git a/src/Primitive_basis/prim.irp.f b/src/Primitive_basis/prim.irp.f deleted file mode 100644 index 28ae0dd6..00000000 --- a/src/Primitive_basis/prim.irp.f +++ /dev/null @@ -1,175 +0,0 @@ -BEGIN_PROVIDER [integer, prim_num] - implicit none - BEGIN_DOC -! Number of uniq primitive basis function - END_DOC - PROVIDE ezfio_filename - call ezfio_get_primitive_basis_prim_num(prim_num) - -END_PROVIDER - - BEGIN_PROVIDER [integer, prim_nucl, (prim_num)] - implicit none - BEGIN_DOC - ! Array of the nuclei on which are attached the primitives - END_DOC - PROVIDE ezfio_filename - integer, allocatable :: tmp(:) - allocate(tmp(prim_num)) - call ezfio_get_primitive_basis_prim_nucl(tmp) - prim_nucl = tmp - END_PROVIDER - - BEGIN_PROVIDER [integer, prim_power, (prim_num,3)] -&BEGIN_PROVIDER [integer, prim_l, (prim_num)] - implicit none - BEGIN_DOC - ! Array of the power of the primitives - ! Array of the L of the primitive - END_DOC - PROVIDE ezfio_filename - integer :: i,j - integer, allocatable :: tmp(:,:) - allocate(tmp(3,prim_num)) - call ezfio_get_primitive_basis_prim_power(tmp) - do i = 1, prim_num - do j = 1,3 - prim_power(i,j) = tmp(j,i) - enddo - enddo - do i = 1, prim_num - prim_l(i) = prim_power(i,1) + prim_power(i,2) + prim_power(i,3) - enddo - END_PROVIDER - - BEGIN_PROVIDER [double precision, prim_expo, (prim_num)] - implicit none - BEGIN_DOC - ! Array of the exponents of the primitives - END_DOC - PROVIDE ezfio_filename - double precision, allocatable :: tmp(:) - allocate(tmp(prim_num)) - call ezfio_get_primitive_basis_prim_expo(tmp) - prim_expo = tmp - END_PROVIDER - - BEGIN_PROVIDER [double precision, prim_overlap, (prim_num,prim_num)] - implicit none - BEGIN_DOC - ! Array of the overlap between the primitives - END_DOC - integer :: i,j - prim_overlap = 0.d0 - double precision :: overlap, overlap_x, overlap_y, overlap_z - double precision :: alpha, beta, c - double precision :: A_center(3), B_center(3) - integer :: power_A(3), power_B(3),dim1 - dim1=100 - do j = 1, prim_num - A_center(1) = nucl_coord( prim_nucl(j), 1 ) - A_center(2) = nucl_coord( prim_nucl(j), 2 ) - A_center(3) = nucl_coord( prim_nucl(j), 3 ) - power_A(1) = prim_power( j, 1 ) - power_A(2) = prim_power( j, 2 ) - power_A(3) = prim_power( j, 3 ) - alpha = prim_expo(j) - do i = 1, prim_num - B_center(1) = nucl_coord( prim_nucl(i), 1 ) - B_center(2) = nucl_coord( prim_nucl(i), 2 ) - B_center(3) = nucl_coord( prim_nucl(i), 3 ) - power_B(1) = prim_power( i, 1 ) - power_B(2) = prim_power( i, 2 ) - power_B(3) = prim_power( i, 3 ) - beta = prim_expo(i) - call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim1) - prim_overlap(i,j) = overlap - enddo - enddo - END_PROVIDER - - BEGIN_PROVIDER [ double precision, prim_ao_coef, (prim_num,ao_num)] - BEGIN_DOC - ! Developement of the primitive on the AO basis - ! prim_ao_coef(i,j) = - END_DOC - implicit none - integer :: i,j,k,l - - prim_ao_coef = 0.d0 - do i = 1, prim_num - do j = 1, ao_num - if(ao_nucl(j).ne.prim_nucl(i))cycle - if(ao_l(j) == prim_l(i))then - if(ao_power(j,1) == prim_power(i,1) .and. & - ao_power(j,2) == prim_power(i,2) .and. & - ao_power(j,3) == prim_power(i,3) )then - ! primitive of the same L than the AO - do k = 1, ao_prim_num(j) - if(ao_expo(j,k) == prim_expo(i))then - prim_ao_coef(i,j) += ao_coef(j,k) - endif - enddo - endif - endif - enddo - enddo - END_PROVIDER - - BEGIN_PROVIDER [ double precision, prim_ao_overlap, (ao_num,ao_num)] - implicit none - BEGIN_dOC - ! Array of the overlap of the AO basis, explicitely calculated on the primitive basis - END_DOC - integer :: i,j,k,l - do i= 1, ao_num - do j = 1,ao_num - prim_ao_overlap(i,j) = 0.d0 - do k = 1, prim_num - do l = 1, prim_num - prim_ao_overlap(i,j) += prim_ao_coef(k,i) * prim_overlap(k,l) * prim_ao_coef(l,j) - enddo - enddo - enddo - enddo - - END_PROVIDER - - BEGIN_PROVIDER [ double precision, prim_mo_coef, (prim_num,mo_tot_num)] - BEGIN_DOC - ! Developement of the primitive on the MO basis - ! prim_ao_coef(i,j) = - END_DOC - - integer :: i,j,k,l - double precision :: coef - prim_mo_coef = 0.d0 - do i = 1, mo_tot_num - do k = 1, ao_num - coef = mo_coef(k,i) - do l = 1, prim_num - prim_mo_coef(l,i) += coef * prim_ao_coef(l,k) - enddo - enddo - enddo - - END_PROVIDER - - BEGIN_PROVIDER [ double precision, prim_mo_overlap, (mo_tot_num,mo_tot_num)] - implicit none - BEGIN_dOC - ! Array of the overlap of the AO basis, explicitely calculated on the primitive basis - END_DOC - integer :: i,j,k,l - do i= 1, mo_tot_num - do j = 1,mo_tot_num - prim_mo_overlap(i,j) = 0.d0 - do k = 1, prim_num - do l = 1, prim_num - prim_mo_overlap(i,j) += prim_mo_coef(k,i) * prim_overlap(k,l) * prim_mo_coef(l,j) - enddo - enddo - enddo - enddo - - END_PROVIDER diff --git a/src/Primitive_basis/primitive_basis.ezfio_config b/src/Primitive_basis/primitive_basis.ezfio_config deleted file mode 100644 index 8988801c..00000000 --- a/src/Primitive_basis/primitive_basis.ezfio_config +++ /dev/null @@ -1,5 +0,0 @@ -primitive_basis - prim_num integer - prim_power integer (3,primitive_basis_prim_num) - prim_expo double precision (primitive_basis_prim_num) - prim_nucl integer (primitive_basis_prim_num) diff --git a/src/Properties/Makefile b/src/Properties/Makefile new file mode 100644 index 00000000..b2ea1de1 --- /dev/null +++ b/src/Properties/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/Properties/NEEDED_MODULES b/src/Properties/NEEDED_MODULES new file mode 100644 index 00000000..8cf33c90 --- /dev/null +++ b/src/Properties/NEEDED_MODULES @@ -0,0 +1 @@ +AOs BiInts Bitmask Dets Electrons Ezfio_files MonoInts MOs Nuclei Output Utils diff --git a/src/Properties/README.rst b/src/Properties/README.rst new file mode 100644 index 00000000..34565220 --- /dev/null +++ b/src/Properties/README.rst @@ -0,0 +1,4 @@ +================= +Properties Module +================= + diff --git a/src/Properties/average.irp.f b/src/Properties/average.irp.f new file mode 100644 index 00000000..b57c8ef2 --- /dev/null +++ b/src/Properties/average.irp.f @@ -0,0 +1,25 @@ +subroutine get_average(array,density,average) + implicit none + double precision, intent(in) :: array(mo_tot_num_align,mo_tot_num) + double precision, intent(in) :: density(mo_tot_num_align,mo_tot_num) + double precision, intent(out):: average + integer :: i,j + BEGIN_DOC +! computes the average value of a pure MONO ELECTRONIC OPERATOR +! whom integrals on the MO basis are stored in "array" +! and with the density is stored in "density" + END_DOC + average = 0.d0 + + !$OMP PARALLEL DO DEFAULT(none) & + !$OMP PRIVATE(i,j) & + !$OMP SHARED(mo_tot_num,array,density) & + !$OMP REDUCTION(+:average) + do i = 1, mo_tot_num + do j = 1, mo_tot_num + average += density(j,i) * array(j,i) + enddo + enddo + !$OMP END PARALLEL DO + +end diff --git a/src/Properties/delta_rho.irp.f b/src/Properties/delta_rho.irp.f new file mode 100644 index 00000000..3cfe136c --- /dev/null +++ b/src/Properties/delta_rho.irp.f @@ -0,0 +1,224 @@ + BEGIN_PROVIDER [integer, N_z_pts] +&BEGIN_PROVIDER [double precision, z_min] +&BEGIN_PROVIDER [double precision, z_max] +&BEGIN_PROVIDER [double precision, delta_z] + implicit none + z_min = -20.d0 + z_max = 20.d0 + delta_z = 0.1d0 + N_z_pts = (z_max - z_min)/delta_z + print*,'N_z_pts = ',N_z_pts + +END_PROVIDER + + +BEGIN_PROVIDER [double precision, integrated_delta_rho_all_points, (N_z_pts)] + BEGIN_DOC +! +! integrated_rho(alpha,z) - integrated_rho(beta,z) for all the z points +! chosen +! + END_DOC + implicit none + integer :: i,j,k,l,i_z,h + double precision :: z,function_integrated_delta_rho,c_k,c_j,n_i_h,accu + integrated_delta_rho_all_points = 0.d0 + !$OMP PARALLEL DO DEFAULT(none) & + !$OMP PRIVATE(i,h,j,k,c_j,c_k,n_i_h,i_z) & + !$OMP SHARED(mo_tot_num,ao_num,mo_coef, & + !$OMP ao_integrated_delta_rho_all_points,one_body_spin_density_mo,integrated_delta_rho_all_points,N_z_pts) + do i_z = 1, N_z_pts + do i = 1, mo_tot_num + do h = 1, mo_tot_num + n_i_h = one_body_spin_density_mo(i,h) + if(dabs(n_i_h).lt.1.d-10)cycle + do j = 1, ao_num + c_j = mo_coef(j,i) ! coefficient of the ith MO on the jth AO + do k = 1, ao_num + c_k = mo_coef(k,h) ! coefficient of the hth MO on the kth AO + integrated_delta_rho_all_points(i_z) += c_k * c_j * n_i_h * ao_integrated_delta_rho_all_points(j,k,i_z) + enddo + enddo + enddo + enddo + enddo + !$OMP END PARALLEL DO + + z = z_min + accu = 0.d0 + do i = 1, N_z_pts + accu += integrated_delta_rho_all_points(i) + write(i_unit_integrated_delta_rho,*)z,integrated_delta_rho_all_points(i),accu + z += delta_z + enddo + print*,'sum of integrated_delta_rho = ',accu + +END_PROVIDER + + + + +BEGIN_PROVIDER [ double precision, ao_integrated_delta_rho_all_points, (ao_num_align, ao_num, N_z_pts)] + BEGIN_DOC +! array of the overlap in x,y between the AO function and integrated between [z,z+dz] in the z axis +! for all the z points that are given (N_z_pts) + END_DOC + implicit none + integer :: i,j,n,l + double precision :: f,accu + integer :: dim1 + double precision :: overlap, overlap_x, overlap_y, overlap_z + double precision :: alpha, beta, c + double precision :: A_center(3), B_center(3) + integer :: power_A(3), power_B(3) + integer :: i_z + double precision :: z,SABpartial,accu_x,accu_y,accu_z + dim1=100 + z = z_min + do i_z = 1, N_z_pts + !$OMP PARALLEL DO DEFAULT(none) & + !$OMP PRIVATE(i,j,n,l,A_center,power_A,B_center,power_B,accu_z, & + !$OMP overlap_x,overlap_y,overlap_z,overlap,c,alpha,beta) & + !$OMP SHARED(ao_num,nucl_coord,ao_nucl,ao_power,ao_prim_num,ao_expo_transp,ao_coef_transp, & + !$OMP ao_integrated_delta_rho_all_points,N_z_pts,dim1,i_z,z,delta_z) + do j=1,ao_num + A_center(1) = nucl_coord( ao_nucl(j), 1 ) + A_center(2) = nucl_coord( ao_nucl(j), 2 ) + A_center(3) = nucl_coord( ao_nucl(j), 3 ) + power_A(1) = ao_power( j, 1 ) + power_A(2) = ao_power( j, 2 ) + power_A(3) = ao_power( j, 3 ) + do i= 1,ao_num + B_center(1) = nucl_coord( ao_nucl(i), 1 ) + B_center(2) = nucl_coord( ao_nucl(i), 2 ) + B_center(3) = nucl_coord( ao_nucl(i), 3 ) + power_B(1) = ao_power( i, 1 ) + power_B(2) = ao_power( i, 2 ) + power_B(3) = ao_power( i, 3 ) + + accu_z = 0.d0 + do n = 1,ao_prim_num(j) + alpha = ao_expo_transp(n,j) + do l = 1, ao_prim_num(i) + beta = ao_expo_transp(l,i) + call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim1) + + c = ao_coef_transp(n,j) * ao_coef_transp(l,i) + accu_z += c* overlap_x * overlap_y * SABpartial(z,z+delta_z,A_center,B_center,power_A,power_B,alpha,beta) + enddo + enddo + ao_integrated_delta_rho_all_points(i,j,i_z) = accu_z + enddo + enddo + !$OMP END PARALLEL DO + z += delta_z + enddo +END_PROVIDER + +BEGIN_PROVIDER [integer, i_unit_integrated_delta_rho] + implicit none + BEGIN_DOC +! fortran unit for the writing of the integrated delta_rho + END_DOC + integer :: getUnitAndOpen + character*(128) :: output_i_unit_integrated_delta_rho + output_i_unit_integrated_delta_rho=trim(ezfio_filename)//'/properties/delta_rho' + i_unit_integrated_delta_rho= getUnitAndOpen(output_i_unit_integrated_delta_rho,'w') + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, ao_integrated_delta_rho_one_point, (ao_num_align, ao_num )] + BEGIN_DOC +! array of the overlap in x,y between the AO function and integrated between [z,z+dz] in the z axis +! for one specific z point + END_DOC + implicit none + integer :: i,j,n,l + double precision :: f + integer :: dim1 + double precision :: overlap, overlap_x, overlap_y, overlap_z + double precision :: alpha, beta, c + double precision :: A_center(3), B_center(3) + integer :: power_A(3), power_B(3) + integer :: i_z + double precision :: z,SABpartial,accu_z + dim1=100 + z = z_one_point + !$OMP PARALLEL DO DEFAULT(none) & + !$OMP PRIVATE(i,j,n,l,A_center,power_A,B_center,power_B,accu_z, & + !$OMP overlap_x,overlap_y,overlap_z,overlap,c,alpha,beta) & + !$OMP SHARED(ao_num,nucl_coord,ao_nucl,ao_power,ao_prim_num,ao_expo_transp,ao_coef_transp, & + !$OMP ao_integrated_delta_rho_one_point,dim1,z,delta_z) + do j=1,ao_num + A_center(1) = nucl_coord( ao_nucl(j), 1 ) + A_center(2) = nucl_coord( ao_nucl(j), 2 ) + A_center(3) = nucl_coord( ao_nucl(j), 3 ) + power_A(1) = ao_power( j, 1 ) + power_A(2) = ao_power( j, 2 ) + power_A(3) = ao_power( j, 3 ) + do i= 1,ao_num + B_center(1) = nucl_coord( ao_nucl(i), 1 ) + B_center(2) = nucl_coord( ao_nucl(i), 2 ) + B_center(3) = nucl_coord( ao_nucl(i), 3 ) + power_B(1) = ao_power( i, 1 ) + power_B(2) = ao_power( i, 2 ) + power_B(3) = ao_power( i, 3 ) + + accu_z = 0.d0 + do n = 1,ao_prim_num(j) + alpha = ao_expo_transp(n,j) + do l = 1, ao_prim_num(i) + beta = ao_expo_transp(l,i) + call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim1) + + c = ao_coef_transp(n,j) * ao_coef_transp(l,i) + accu_z += c* overlap_x * overlap_y * SABpartial(z,z+delta_z,A_center,B_center,power_A,power_B,alpha,beta) + enddo + enddo + ao_integrated_delta_rho_one_point(i,j) = accu_z + enddo + enddo + !$OMP END PARALLEL DO +END_PROVIDER + +BEGIN_PROVIDER [double precision, mo_integrated_delta_rho_one_point, (mo_tot_num_align,mo_tot_num)] + BEGIN_DOC +! +! array of the integrals needed of integrated_rho(alpha,z) - integrated_rho(beta,z) for z = z_one_point +! on the MO basis +! + END_DOC + implicit none + integer :: i,j,k,l,i_z,h + double precision :: z,function_integrated_delta_rho,c_k,c_j + mo_integrated_delta_rho_one_point = 0.d0 + !$OMP PARALLEL DO DEFAULT(none) & + !$OMP PRIVATE(i,j,h,k,c_j,c_k) & + !$OMP SHARED(mo_tot_num,ao_num,mo_coef, & + !$OMP mo_integrated_delta_rho_one_point, ao_integrated_delta_rho_one_point) + do i = 1, mo_tot_num + do h = 1, mo_tot_num + do j = 1, ao_num + c_j = mo_coef(j,i) ! coefficient of the jth AO on the ith MO + do k = 1, ao_num + c_k = mo_coef(k,h) ! coefficient of the kth AO on the hth MO + mo_integrated_delta_rho_one_point(i,h) += c_k * c_j * ao_integrated_delta_rho_one_point(j,k) + enddo + enddo + enddo + enddo + !$OMP END PARALLEL DO +END_PROVIDER +BEGIN_PROVIDER [ double precision, integrated_delta_rho_one_point] + implicit none + BEGIN_DOC +! +! integral (x,y) and (z,z+delta_z) of rho(alpha) - rho(beta) +! on the MO basis +! + END_DOC + double precision :: average +call get_average(mo_integrated_delta_rho_one_point,one_body_spin_density_mo,average) + integrated_delta_rho_one_point = average +END_PROVIDER + diff --git a/src/Properties/need.irp.f b/src/Properties/need.irp.f new file mode 100644 index 00000000..eb4dfe34 --- /dev/null +++ b/src/Properties/need.irp.f @@ -0,0 +1,289 @@ + + double precision function SABpartial(zA,zB,A,B,nA,nB,gamA,gamB) + implicit double precision(a-h,o-z) + dimension nA(3),nB(3) + dimension A(3),B(3) + gamtot=gamA+gamB + SABpartial=1.d0 + + l=3 + u=gamA/gamtot*A(l)+gamB/gamtot*B(l) + arg=gamtot*u**2-gamA*A(l)**2-gamB*B(l)**2 + alpha=dexp(arg) + &/gamtot**((1.d0+dfloat(nA(l))+dfloat(nB(l)))/2.d0) + wA=dsqrt(gamtot)*(u-A(l)) + wB=dsqrt(gamtot)*(u-B(l)) + boundA=dsqrt(gamtot)*(zA-u) + boundB=dsqrt(gamtot)*(zB-u) + + accu=0.d0 + do n=0,nA(l) + do m=0,nB(l) + integ=nA(l)+nB(l)-n-m + accu=accu + & +wA**n*wB**m*binom(nA(l),n)*binom(nB(l),m) + & *(rinteg(integ,boundB)-rinteg(integ,boundA)) + enddo + enddo + SABpartial=SABpartial*accu*alpha + end + + double precision function rintgauss(n) + implicit double precision(a-h,o-z) + rintgauss=dsqrt(dacos(-1.d0)) + if(n.eq.0)return + if(n.eq.1)then + rintgauss=0.d0 + return + endif + if(iand(n,1).eq.1)then + rintgauss=0.d0 + return + endif + rintgauss=rintgauss/2.d0**(n/2) + rintgauss=rintgauss*ddfact2(n-1) + end + + double precision function rinteg(n,u) + implicit double precision(a-h,o-z) + include 'constants.F' +! pi=dacos(-1.d0) + ichange=1 + factor=1.d0 + if(u.lt.0.d0)then + u=-u + factor=(-1.d0)**(n+1) + ichange=-1 + endif + if(iand(n,1).eq.0)then + rinteg=0.d0 + do l=0,n-2,2 + prod=b_coef(l,u) + do k=l+2,n-2,2 + prod=prod*a_coef(k) + enddo + rinteg=rinteg+prod + enddo + prod=dsqrt(pi)/2.d0*erf0(u) + do k=0,n-2,2 + prod=prod*a_coef(k) + enddo + rinteg=rinteg+prod + endif + + if(iand(n,1).eq.1)then + rinteg=0.d0 + do l=1,n-2,2 + prod=b_coef(l,u) + do k=l+2,n-2,2 + prod=prod*a_coef(k) + enddo + rinteg=rinteg+prod + enddo + prod=0.5d0*(1.d0-dexp(-u**2)) + do k=1,n-2,2 + prod=prod*a_coef(k) + enddo + rinteg=rinteg+prod + endif + + rinteg=rinteg*factor + + if(ichange.eq.-1)u=-u + + end + +! +! +! +! +! +! +! +! + double precision function erf0(x) + implicit double precision (a-h,o-z) + if(x.lt.0.d0)then + erf0=-gammp(0.5d0,x**2) + else + erf0=gammp(0.5d0,x**2) + endif + end + + +! +! +! +! +! +! +! +! +! +! +! +! gcf +! gser +! +! +! + double precision function gammp(a,x) + implicit double precision (a-h,o-z) + if(x.lt.0..or.a.le.0.)pause + if(x.lt.a+1.)then + call gser(gammp,a,x,gln) + else + call gcf(gammcf,a,x,gln) + gammp=1.-gammcf + endif + return + end +! +! + + +! +! +! +! +! +! +! +! +! +! +! +! gammp +! +! +! + subroutine gser(gamser,a,x,gln) + implicit double precision (a-h,o-z) + parameter (itmax=100,eps=3.e-7) + gln=gammln(a) + if(x.le.0.)then + if(x.lt.0.)pause + gamser=0. + return + endif + ap=a + sum=1./a + del=sum + do 11 n=1,itmax + ap=ap+1. + del=del*x/ap + sum=sum+del + if(abs(del).lt.abs(sum)*eps)go to 1 +11 continue + pause 'a too large, itmax too small' +1 gamser=sum*exp(-x+a*log(x)-gln) + return + end +! + +! +! +! +! +! +! +! +! +! +! +! +! +! gammp +! +! +! + subroutine gcf(gammcf,a,x,gln) + implicit double precision (a-h,o-z) + parameter (itmax=100,eps=3.e-7) + gln=gammln(a) + gold=0. + a0=1. + a1=x + b0=0. + b1=1. + fac=1. + do 11 n=1,itmax + an=float(n) + ana=an-a + a0=(a1+a0*ana)*fac + b0=(b1+b0*ana)*fac + anf=an*fac + a1=x*a0+anf*a1 + b1=x*b0+anf*b1 + if(a1.ne.0.)then + fac=1./a1 + g=b1*fac + if(abs((g-gold)/g).lt.eps)go to 1 + gold=g + endif +11 continue + pause 'a too large, itmax too small' +1 gammcf=exp(-x+a*log(x)-gln)*g + return + end + +! +! + double precision function ddfact2(n) + implicit double precision(a-h,o-z) + if(iand(n,1).eq.0)stop 'error in ddfact2' + ddfact2=1.d0 + do i=1,n,2 + ddfact2=ddfact2*dfloat(i) + enddo + end + + double precision function a_coef(n) + implicit double precision(a-h,o-z) + a_coef=dfloat(n+1)/2.d0 + end + + double precision function b_coef(n,u) + implicit double precision(a-h,o-z) + b_coef=-0.5d0*u**(n+1)*dexp(-u**2) + end + +! +! +! +! +! +! +! +! + double precision function gammln(xx) + implicit double precision (a-h,o-z) + real*8 cof(6),stp,half,one,fpf,x,tmp,ser + data cof,stp/76.18009173d0,-86.50532033d0,24.01409822d0, + * -1.231739516d0,.120858003d-2,-.536382d-5,2.50662827465d0/ + data half,one,fpf/0.5d0,1.0d0,5.5d0/ + x=xx-one + tmp=x+fpf + tmp=(x+half)*log(tmp)-tmp + ser=one + do 11 j=1,6 + x=x+one + ser=ser+cof(j)/x +11 continue + gammln=tmp+log(stp*ser) + return + end +! +! diff --git a/src/Properties/options.irp.f b/src/Properties/options.irp.f new file mode 100644 index 00000000..0fd5a4c1 --- /dev/null +++ b/src/Properties/options.irp.f @@ -0,0 +1,13 @@ +BEGIN_SHELL [ /usr/bin/python ] +from ezfio_with_default import EZFIO_Provider +T = EZFIO_Provider() +T.set_type ( "double precision" ) +T.set_name ( "z_one_point" ) +T.set_doc ( "z point on which the integrated delta rho is calculated" ) +T.set_ezfio_dir ( "properties" ) +T.set_ezfio_name( "z_one_point" ) +T.set_output ( "output_full_ci" ) +print T + +END_SHELL + diff --git a/src/Properties/properties.ezfio_config b/src/Properties/properties.ezfio_config new file mode 100644 index 00000000..018b56d0 --- /dev/null +++ b/src/Properties/properties.ezfio_config @@ -0,0 +1,2 @@ +properties + z_one_point double precision diff --git a/src/Properties/properties.irp.f b/src/Properties/properties.irp.f new file mode 100644 index 00000000..aad310f9 --- /dev/null +++ b/src/Properties/properties.irp.f @@ -0,0 +1,51 @@ +BEGIN_PROVIDER [double precision, average_position,(3)] + implicit none + BEGIN_DOC + ! average_position(1) = + ! average_position(2) = + ! average_position(3) = + END_DOC + integer :: i,j + average_position = 0.d0 + !$OMP PARALLEL DO DEFAULT(none) & + !$OMP PRIVATE(i,j) & + !$OMP SHARED(mo_tot_num,mo_dipole_x,mo_dipole_y,mo_dipole_z,one_body_dm_mo) & + !$OMP REDUCTION(+:average_position) + do i = 1, mo_tot_num + do j = 1, mo_tot_num + average_position(1) += one_body_dm_mo(j,i) * mo_dipole_x(j,i) + average_position(2) += one_body_dm_mo(j,i) * mo_dipole_y(j,i) + average_position(3) += one_body_dm_mo(j,i) * mo_dipole_z(j,i) + enddo + enddo + !$OMP END PARALLEL DO +!call test_average_value(mo_dipole_z,average_position(3)) + +END_PROVIDER + + +BEGIN_PROVIDER [double precision, average_spread,(3)] + implicit none + BEGIN_DOC + ! average_spread(1) = + ! average_spread(2) = + ! average_spread(3) = + END_DOC + integer :: i,j + average_spread = 0.d0 + !$OMP PARALLEL DO DEFAULT(none) & + !$OMP PRIVATE(i,j) & + !$OMP SHARED(mo_tot_num,mo_spread_x,mo_spread_y,mo_spread_z,one_body_dm_mo) & + !$OMP REDUCTION(+:average_spread) + do i = 1, mo_tot_num + do j = 1, mo_tot_num + average_spread(1) += one_body_dm_mo(j,i) * mo_spread_x(j,i) + average_spread(2) += one_body_dm_mo(j,i) * mo_spread_y(j,i) + average_spread(3) += one_body_dm_mo(j,i) * mo_spread_z(j,i) + enddo + enddo + !$OMP END PARALLEL DO +!call test_average_value(mo_spread_z,average_spread(3)) + +END_PROVIDER + diff --git a/src/Properties/routines_test.irp.f b/src/Properties/routines_test.irp.f new file mode 100644 index 00000000..231c6f2d --- /dev/null +++ b/src/Properties/routines_test.irp.f @@ -0,0 +1,89 @@ + + +subroutine test_average_value(array,value) + implicit none + double precision, intent(in) :: array(mo_tot_num_align,mo_tot_num) + double precision, intent(in) :: value + double precision :: tmp,hij + integer :: i,j + tmp = 0.d0 + do i = 1, N_det + do j = 1, N_det + call i_O1_j(array,psi_det(1,1,i),psi_det(1,1,j),N_int,hij) + tmp+= hij * psi_coef(i,1) * psi_coef(j,1) + enddo + enddo + if(dabs(tmp-value).ge.1.d-8)then + print*,'subroutine test_average_value' + print*,'PB WITH AVERAGE VALUE !!!!!!!!!' + print*,' = ',tmp + print*,'value = ',value + stop + endif +end + +subroutine test_average_value_alpha_beta(array,value) + implicit none + double precision, intent(in) :: array(mo_tot_num_align,mo_tot_num) + double precision, intent(in) :: value + double precision :: tmp,hij + integer :: i,j + tmp = 0.d0 + do i = 1, N_det + do j = 1, N_det + call i_O1_j_alpha_beta(array,psi_det(1,1,i),psi_det(1,1,j),N_int,hij) + tmp+= hij * psi_coef(i,1) * psi_coef(j,1) + enddo + enddo + print*,'tmp = ',tmp + print*,'value = ',value + tmp = 0.d0 + do i = 1, N_det + call i_O1_psi_alpha_beta(array,psi_det(1,1,i),psi_det,psi_coef,N_int,N_det,psi_det_size,N_states,hij) + tmp+= hij * psi_coef(i,1) + enddo + print*,'tmp = ',tmp + print*,'value = ',value + if(dabs(tmp-value).ge.1.d-8)then + print*,'subroutine test_average_value' + print*,'PB WITH AVERAGE VALUE !!!!!!!!!' + print*,' = ',tmp + print*,'value = ',value + stop + endif +end + +subroutine test_dm(tmp_array) + use bitmasks + implicit none + double precision,intent(out) :: tmp_array(mo_tot_num,mo_tot_num) + double precision :: phase + integer :: i,j,ispin,k + integer :: degree + integer :: exc(0:2,2,2),h1,p1,h2,p2,s1,s2 + integer :: occ_det(N_int*bit_kind_size,2) + integer :: tmp + double precision :: accu + + tmp_array = 0.d0 + + do i = 1, N_det + call bitstring_to_list(psi_det(1,1,i), occ_det(1,1), tmp, N_int) + call bitstring_to_list(psi_det(1,2,i), occ_det(1,2), tmp, N_int) + do ispin = 1, 2 + do k = 1, elec_num_tab(ispin) + tmp_array(occ_det(k,ispin),occ_det(k,ispin)) += psi_coef(i,1) * psi_coef(i,1) + enddo + enddo + do j = i+1, N_det + call get_excitation_degree(psi_det(1,1,i),psi_det(1,1,j),degree,N_int) + if (degree.ne.1)cycle + call get_excitation(psi_det(1,1,i),psi_det(1,1,j),exc,degree,phase,N_int) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + tmp_array(h1,p1) += phase * psi_coef(i,1) * psi_coef(j,1) + tmp_array(p1,h1) = tmp_array(h1,p1) + enddo + enddo + + +end diff --git a/src/Properties/slater_rules_mono_electronic.irp.f b/src/Properties/slater_rules_mono_electronic.irp.f new file mode 100644 index 00000000..3e4d94bc --- /dev/null +++ b/src/Properties/slater_rules_mono_electronic.irp.f @@ -0,0 +1,339 @@ +subroutine i_O1_j(array,key_i,key_j,Nint,hij) + use bitmasks + implicit none + BEGIN_DOC + ! Returns where i and j are determinants + ! and O1 is a ONE BODY OPERATOR + ! array is the array of the mono electronic operator + ! on the MO basis + END_DOC + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2) + double precision, intent(out) :: hij + double precision, intent(in) :: array(mo_tot_num_align,mo_tot_num) + + integer :: exc(0:2,2,2) + integer :: degree + integer :: m,p + double precision :: diag_O1_mat_elem, phase + + ASSERT (Nint > 0) + ASSERT (Nint == N_int) + ASSERT (sum(popcnt(key_i(:,1))) == elec_alpha_num) + ASSERT (sum(popcnt(key_i(:,2))) == elec_beta_num) + ASSERT (sum(popcnt(key_j(:,1))) == elec_alpha_num) + ASSERT (sum(popcnt(key_j(:,2))) == elec_beta_num) + + hij = 0.d0 + !DEC$ FORCEINLINE + call get_excitation_degree(key_i,key_j,degree,Nint) + select case (degree) + case (2) + hij = 0.d0 + case (1) + call get_mono_excitation(key_i,key_j,exc,phase,Nint) + if (exc(0,1,1) == 1) then + ! Mono alpha + m = exc(1,1,1) + p = exc(1,2,1) + else + ! Mono beta + m = exc(1,1,2) + p = exc(1,2,2) + endif + hij = phase* array(m,p) + + case (0) + hij = diag_O1_mat_elem(array,key_i,Nint) + end select +end + + +subroutine i_O1_psi(array,key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array) + use bitmasks + implicit none + integer, intent(in) :: Nint, Ndet,Ndet_max,Nstate + double precision, intent(in) :: array(mo_tot_num_align,mo_tot_num) + integer(bit_kind), intent(in) :: keys(Nint,2,Ndet) + integer(bit_kind), intent(in) :: key(Nint,2) + double precision, intent(in) :: coef(Ndet_max,Nstate) + double precision, intent(out) :: i_H_psi_array(Nstate) + + integer :: i, ii,j + double precision :: phase + integer :: exc(0:2,2,2) + double precision :: hij + integer :: idx(0:Ndet) + BEGIN_DOC + ! for the various Nstates + ! and O1 is a ONE BODY OPERATOR + ! array is the array of the mono electronic operator + ! on the MO basis + END_DOC + + ASSERT (Nint > 0) + ASSERT (N_int == Nint) + ASSERT (Nstate > 0) + ASSERT (Ndet > 0) + ASSERT (Ndet_max >= Ndet) + i_H_psi_array = 0.d0 + call filter_connected_mono(keys,key,Nint,Ndet,idx) + do ii=1,idx(0) + i = idx(ii) + !DEC$ FORCEINLINE + call i_O1_j(array,keys(1,1,i),key,Nint,hij) + do j = 1, Nstate + i_H_psi_array(j) = i_H_psi_array(j) + coef(i,j)*hij + enddo + enddo +end + +double precision function diag_O1_mat_elem(array,det_in,Nint) + use bitmasks + implicit none + BEGIN_DOC + ! Computes + END_DOC + integer,intent(in) :: Nint + integer(bit_kind),intent(in) :: det_in(Nint,2) + double precision, intent(in) :: array(mo_tot_num_align,mo_tot_num) + + integer :: i, ispin,tmp + integer :: occ_det(Nint*bit_kind_size,2) + + ASSERT (Nint > 0) + ASSERT (sum(popcnt(det_in(:,1))) == elec_alpha_num) + ASSERT (sum(popcnt(det_in(:,2))) == elec_beta_num) + + call bitstring_to_list(det_in(1,1), occ_det(1,1), tmp, Nint) + call bitstring_to_list(det_in(1,2), occ_det(1,2), tmp, Nint) + diag_O1_mat_elem = 0.d0 + do ispin = 1, 2 + do i = 1, elec_num_tab(ispin) + diag_O1_mat_elem += array(occ_det(i,ispin),occ_det(i,ispin)) + enddo + enddo +end + + +subroutine i_O1_psi_alpha_beta(array,key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array) + use bitmasks + implicit none + integer, intent(in) :: Nint, Ndet,Ndet_max,Nstate + double precision, intent(in) :: array(mo_tot_num_align,mo_tot_num) + integer(bit_kind), intent(in) :: keys(Nint,2,Ndet) + integer(bit_kind), intent(in) :: key(Nint,2) + double precision, intent(in) :: coef(Ndet_max,Nstate) + double precision, intent(out) :: i_H_psi_array(Nstate) + + integer :: i, ii,j + double precision :: phase + integer :: exc(0:2,2,2) + double precision :: hij + integer :: idx(0:Ndet) + BEGIN_DOC + ! for the various Nstates + ! and O1 is a ONE BODY OPERATOR + ! array is the array of the mono electronic operator + ! on the MO basis + END_DOC + + ASSERT (Nint > 0) + ASSERT (N_int == Nint) + ASSERT (Nstate > 0) + ASSERT (Ndet > 0) + ASSERT (Ndet_max >= Ndet) + i_H_psi_array = 0.d0 + call filter_connected_mono(keys,key,Nint,Ndet,idx) + do ii=1,idx(0) + i = idx(ii) + !DEC$ FORCEINLINE + call i_O1_j_alpha_beta(array,keys(1,1,i),key,Nint,hij) + do j = 1, Nstate + i_H_psi_array(j) = i_H_psi_array(j) + coef(i,j)*hij + enddo + enddo +end + +subroutine i_O1_j_alpha_beta(array,key_i,key_j,Nint,hij) + use bitmasks + implicit none + BEGIN_DOC + ! Returns where i and j are determinants + ! and O1 is a ONE BODY OPERATOR + ! array is the array of the mono electronic operator + ! on the MO basis + END_DOC + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2) + double precision, intent(out) :: hij + double precision, intent(in) :: array(mo_tot_num_align,mo_tot_num) + + integer :: exc(0:2,2,2) + integer :: degree + integer :: m,p + double precision :: diag_O1_mat_elem_alpha_beta, phase + + ASSERT (Nint > 0) + ASSERT (Nint == N_int) + ASSERT (sum(popcnt(key_i(:,1))) == elec_alpha_num) + ASSERT (sum(popcnt(key_i(:,2))) == elec_beta_num) + ASSERT (sum(popcnt(key_j(:,1))) == elec_alpha_num) + ASSERT (sum(popcnt(key_j(:,2))) == elec_beta_num) + + hij = 0.d0 + !DEC$ FORCEINLINE + call get_excitation_degree(key_i,key_j,degree,Nint) + select case (degree) + case (2) + hij = 0.d0 + case (1) + call get_mono_excitation(key_i,key_j,exc,phase,Nint) + if (exc(0,1,1) == 1) then + ! Mono alpha + m = exc(1,1,1) + p = exc(1,2,1) + hij = phase* array(m,p) + else + ! Mono beta + m = exc(1,1,2) + p = exc(1,2,2) + hij = -phase* array(m,p) + endif + + case (0) + hij = diag_O1_mat_elem_alpha_beta(array,key_i,Nint) + end select +end + + +double precision function diag_O1_mat_elem_alpha_beta(array,det_in,Nint) + use bitmasks + implicit none + BEGIN_DOC + ! Computes + END_DOC + integer,intent(in) :: Nint + integer(bit_kind),intent(in) :: det_in(Nint,2) + double precision, intent(in) :: array(mo_tot_num_align,mo_tot_num) + + integer :: i, ispin,tmp + integer :: occ_det(Nint*bit_kind_size,2) + + ASSERT (Nint > 0) + ASSERT (sum(popcnt(det_in(:,1))) == elec_alpha_num) + ASSERT (sum(popcnt(det_in(:,2))) == elec_beta_num) + + call bitstring_to_list(det_in(1,1), occ_det(1,1), tmp, Nint) + call bitstring_to_list(det_in(1,2), occ_det(1,2), tmp, Nint) + diag_O1_mat_elem_alpha_beta = 0.d0 + ispin = 1 + do i = 1, elec_num_tab(ispin) + diag_O1_mat_elem_alpha_beta += array(occ_det(i,ispin),occ_det(i,ispin)) + enddo + ispin = 2 + do i = 1, elec_num_tab(ispin) + diag_O1_mat_elem_alpha_beta -= array(occ_det(i,ispin),occ_det(i,ispin)) + enddo +end + +subroutine filter_connected_mono(key1,key2,Nint,sze,idx) + use bitmasks + implicit none + BEGIN_DOC + ! Filters out the determinants that are not connected through PURE + ! + ! MONO EXCITATIONS OPERATORS (a^{\dagger}j a_i) + ! + ! returns the array idx which contains the index of the + ! + ! determinants in the array key1 that interact + ! + ! via some PURE MONO EXCITATIONS OPERATORS + ! + ! idx(0) is the number of determinants that interact with key1 + END_DOC + integer, intent(in) :: Nint, sze + integer(bit_kind), intent(in) :: key1(Nint,2,sze) + integer(bit_kind), intent(in) :: key2(Nint,2) + integer, intent(out) :: idx(0:sze) + + integer :: i,j,l + integer :: degree_x2 + + ASSERT (Nint > 0) + ASSERT (sze >= 0) + + l=1 + + if (Nint==1) then + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = popcnt( xor( key1(1,1,i), key2(1,1))) & + + popcnt( xor( key1(1,2,i), key2(1,2))) + if (degree_x2 > 3) then + cycle + else + idx(l) = i + l = l+1 + endif + enddo + + else if (Nint==2) then + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + & + popcnt(xor( key1(2,1,i), key2(2,1))) + & + popcnt(xor( key1(1,2,i), key2(1,2))) + & + popcnt(xor( key1(2,2,i), key2(2,2))) + if (degree_x2 > 3) then + cycle + else + idx(l) = i + l = l+1 + endif + enddo + + else if (Nint==3) then + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + & + popcnt(xor( key1(1,2,i), key2(1,2))) + & + popcnt(xor( key1(2,1,i), key2(2,1))) + & + popcnt(xor( key1(2,2,i), key2(2,2))) + & + popcnt(xor( key1(3,1,i), key2(3,1))) + & + popcnt(xor( key1(3,2,i), key2(3,2))) + if (degree_x2 > 3) then + cycle + else + idx(l) = i + l = l+1 + endif + enddo + + else + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = 0 + !DEC$ LOOP COUNT MIN(4) + do j=1,Nint + degree_x2 = degree_x2+ popcnt(xor( key1(j,1,i), key2(j,1))) +& + popcnt(xor( key1(j,2,i), key2(j,2))) + if (degree_x2 > 3) then + exit + endif + enddo + if (degree_x2 <= 3) then + idx(l) = i + l = l+1 + endif + enddo + + endif + idx(0) = l-1 +end + diff --git a/src/Selectors_full/selectors.irp.f b/src/Selectors_full/selectors.irp.f index 9de181eb..986241f5 100644 --- a/src/Selectors_full/selectors.irp.f +++ b/src/Selectors_full/selectors.irp.f @@ -60,6 +60,18 @@ END_PROVIDER enddo END_PROVIDER + BEGIN_PROVIDER [ double precision, psi_selectors_diag_h_mat, (psi_selectors_size) ] + implicit none + BEGIN_DOC + ! Diagonal elements of the H matrix for each selectors + END_DOC + integer :: i + double precision :: diag_H_mat_elem + do i = 1, N_det_selectors + psi_selectors_diag_h_mat(i) = diag_H_mat_elem(psi_selectors(1,1,i),N_int) + enddo + END_PROVIDER + BEGIN_PROVIDER [ integer(bit_kind), psi_selectors_ab, (N_int,2,psi_selectors_size) ] &BEGIN_PROVIDER [ double precision, psi_selectors_coef_ab, (psi_selectors_size,N_states) ] diff --git a/src/Selectors_no_sorted/ASSUMPTIONS.rst b/src/Selectors_no_sorted/ASSUMPTIONS.rst new file mode 100644 index 00000000..e69de29b diff --git a/src/Selectors_no_sorted/Makefile b/src/Selectors_no_sorted/Makefile new file mode 100644 index 00000000..b2ea1de1 --- /dev/null +++ b/src/Selectors_no_sorted/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/Selectors_no_sorted/NEEDED_MODULES b/src/Selectors_no_sorted/NEEDED_MODULES new file mode 100644 index 00000000..26097b8b --- /dev/null +++ b/src/Selectors_no_sorted/NEEDED_MODULES @@ -0,0 +1 @@ +AOs BiInts Bitmask Dets Electrons Ezfio_files MonoInts MOs Nuclei Output Utils diff --git a/src/Selectors_no_sorted/README.rst b/src/Selectors_no_sorted/README.rst new file mode 100644 index 00000000..79b78499 --- /dev/null +++ b/src/Selectors_no_sorted/README.rst @@ -0,0 +1,4 @@ +========================== +Selectors_no_sorted Module +========================== + diff --git a/src/Selectors_no_sorted/e_corr_selectors.irp.f b/src/Selectors_no_sorted/e_corr_selectors.irp.f new file mode 100644 index 00000000..952e1c23 --- /dev/null +++ b/src/Selectors_no_sorted/e_corr_selectors.irp.f @@ -0,0 +1,79 @@ + +use bitmasks + BEGIN_PROVIDER [integer, exc_degree_per_selectors, (N_det_selectors)] +&BEGIN_PROVIDER [integer, double_index_selectors, (N_det_selectors)] +&BEGIN_PROVIDER [integer, n_double_selectors] + implicit none + BEGIN_DOC + ! degree of excitation respect to Hartree Fock for the wave function + ! + ! for the all the selectors determinants + ! + ! double_index_selectors = list of the index of the double excitations + ! + ! n_double_selectors = number of double excitations in the selectors determinants + END_DOC + integer :: i,degree + n_double_selectors = 0 + do i = 1, N_det_selectors + call get_excitation_degree(psi_selectors(1,1,i),ref_bitmask,degree,N_int) + exc_degree_per_selectors(i) = degree + if(degree==2)then + n_double_selectors += 1 + double_index_selectors(n_double_selectors) =i + endif + enddo +END_PROVIDER + + BEGIN_PROVIDER[double precision, coef_hf_selector] + &BEGIN_PROVIDER[double precision, inv_selectors_coef_hf] + &BEGIN_PROVIDER[double precision, inv_selectors_coef_hf_squared] + &BEGIN_PROVIDER[double precision, E_corr_per_selectors, (N_det_selectors)] + &BEGIN_PROVIDER[double precision, i_H_HF_per_selectors, (N_det_selectors)] + &BEGIN_PROVIDER[double precision, Delta_E_per_selector, (N_det_selectors)] + &BEGIN_PROVIDER[double precision, E_corr_double_only ] + &BEGIN_PROVIDER[double precision, E_corr_second_order ] + implicit none + BEGIN_DOC + ! energy of correlation per determinant respect to the Hartree Fock determinant + ! + ! for the all the double excitations in the selectors determinants + ! + ! E_corr_per_selectors(i) = * c(D_i)/c(HF) if |D_i> is a double excitation + ! + ! E_corr_per_selectors(i) = -1000.d0 if it is not a double excitation + ! + ! coef_hf_selector = coefficient of the Hartree Fock determinant in the selectors determinants + END_DOC + PROVIDE ref_bitmask_energy psi_selectors ref_bitmask N_int psi_selectors + integer :: i,degree + double precision :: hij,diag_H_mat_elem + E_corr_double_only = 0.d0 + E_corr_second_order = 0.d0 + do i = 1, N_det_selectors + if(exc_degree_per_selectors(i)==2)then + call i_H_j(ref_bitmask,psi_selectors(1,1,i),N_int,hij) + i_H_HF_per_selectors(i) = hij + E_corr_per_selectors(i) = psi_selectors_coef(i,1) * hij + E_corr_double_only += E_corr_per_selectors(i) + E_corr_second_order += hij * hij /(ref_bitmask_energy - diag_H_mat_elem(psi_selectors(1,1,i),N_int)) + elseif(exc_degree_per_selectors(i) == 0)then + coef_hf_selector = psi_selectors_coef(i,1) + E_corr_per_selectors(i) = -1000.d0 + Delta_E_per_selector(i) = 0.d0 + else + E_corr_per_selectors(i) = -1000.d0 + endif + enddo + if (dabs(coef_hf_selector) > 1.d-8) then + inv_selectors_coef_hf = 1.d0/coef_hf_selector + inv_selectors_coef_hf_squared = inv_selectors_coef_hf * inv_selectors_coef_hf + else + inv_selectors_coef_hf = 0.d0 + inv_selectors_coef_hf_squared = 0.d0 + endif + do i = 1,n_double_selectors + E_corr_per_selectors(double_index_selectors(i)) *=inv_selectors_coef_hf + enddo + E_corr_double_only = E_corr_double_only * inv_selectors_coef_hf + END_PROVIDER diff --git a/src/Selectors_no_sorted/selectors.irp.f b/src/Selectors_no_sorted/selectors.irp.f new file mode 100644 index 00000000..d6c20804 --- /dev/null +++ b/src/Selectors_no_sorted/selectors.irp.f @@ -0,0 +1,156 @@ + +use bitmasks + + +BEGIN_PROVIDER [ integer, psi_selectors_size ] + implicit none + psi_selectors_size = psi_det_size +END_PROVIDER + +BEGIN_PROVIDER [ integer, N_det_selectors] + implicit none + BEGIN_DOC + ! For Single reference wave functions, the number of selectors is 1 : the + ! Hartree-Fock determinant + END_DOC + integer :: i + double precision :: norm + call write_time(output_dets) + norm = 0.d0 + N_det_selectors = N_det + N_det_selectors = max(N_det_selectors,1) + call write_int(output_dets,N_det_selectors,'Number of selectors') +END_PROVIDER + + + BEGIN_PROVIDER [ integer(bit_kind), psi_selectors, (N_int,2,psi_selectors_size) ] +&BEGIN_PROVIDER [ double precision, psi_selectors_coef, (psi_selectors_size,N_states) ] + implicit none + BEGIN_DOC + ! Determinants on which we apply for perturbation. + END_DOC + integer :: i,k + + do i=1,N_det_selectors + do k=1,N_int + psi_selectors(k,1,i) = psi_det(k,1,i) + psi_selectors(k,2,i) = psi_det(k,2,i) + enddo + enddo + do k=1,N_states + do i=1,N_det_selectors + psi_selectors_coef(i,k) = psi_coef(i,k) + enddo + enddo +END_PROVIDER + + BEGIN_PROVIDER [ double precision, psi_selectors_diag_h_mat, (psi_selectors_size) ] + implicit none + BEGIN_DOC + ! Diagonal elements of the H matrix for each selectors + END_DOC + integer :: i + double precision :: diag_H_mat_elem + do i = 1, N_det_selectors + psi_selectors_diag_h_mat(i) = diag_H_mat_elem(psi_selectors(1,1,i),N_int) + enddo + END_PROVIDER + + + BEGIN_PROVIDER [ integer(bit_kind), psi_selectors_ab, (N_int,2,psi_selectors_size) ] +&BEGIN_PROVIDER [ double precision, psi_selectors_coef_ab, (psi_selectors_size,N_states) ] +&BEGIN_PROVIDER [ integer, psi_selectors_next_ab, (2,psi_selectors_size) ] + implicit none + BEGIN_DOC + ! Determinants on which we apply . + ! They are sorted by the 3 highest electrons in the alpha part, + ! then by the 3 highest electrons in the beta part to accelerate + ! the research of connected determinants. + END_DOC + integer :: i,j,k + integer, allocatable :: iorder(:) + integer*8, allocatable :: bit_tmp(:) + integer*8, external :: det_search_key + + allocate ( iorder(N_det_selectors), bit_tmp(N_det_selectors) ) + + ! Sort alpha dets + ! --------------- + + integer(bit_kind) :: det_tmp(N_int) + + do i=1,N_det_selectors + iorder(i) = i + call int_of_3_highest_electrons(psi_selectors(1,1,i),bit_tmp(i),N_int) + enddo + call i8sort(bit_tmp,iorder,N_det_selectors) + !DIR$ IVDEP + do i=1,N_det_selectors + do j=1,N_int + psi_selectors_ab(j,1,i) = psi_selectors(j,1,iorder(i)) + psi_selectors_ab(j,2,i) = psi_selectors(j,2,iorder(i)) + enddo + do k=1,N_states + psi_coef_sorted_ab(i,k) = psi_selectors_coef(iorder(i),k) + enddo + enddo + + ! Find next alpha + ! --------------- + + integer :: next + + next = N_det_selectors+1 + psi_selectors_next_ab(1,N_det_selectors) = next + do i=N_det_selectors-1,1,-1 + if (bit_tmp(i) /= bit_tmp(i+1)) then + next = i+1 + endif + psi_selectors_next_ab(1,i) = next + enddo + + ! Sort beta dets + ! -------------- + + integer :: istart, iend + integer(bit_kind), allocatable :: psi_selectors_ab_temp (:,:) + + allocate ( psi_selectors_ab_temp (N_int,N_det_selectors) ) + do i=1,N_det_selectors + do j=1,N_int + psi_selectors_ab_temp(j,i) = psi_selectors_ab(j,2,i) + enddo + iorder(i) = i + call int_of_3_highest_electrons(psi_selectors_ab_temp(1,i),bit_tmp(i),N_int) + enddo + + istart=1 + do while ( istart