From c413f076911568fbb7ef8f423431ab0f1907ed71 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 15 Dec 2020 23:46:19 +0100 Subject: [PATCH 1/9] Minor changes --- TODO.org | 36 ------------ ocaml/qptypes_generator.ml | 6 +- src/ao_two_e_ints/map_integrals.irp.f | 78 ------------------------- src/ao_two_e_ints/two_e_integrals.irp.f | 18 +++--- src/cipsi/selection.irp.f | 12 ++-- src/determinants/h_apply.template.f | 2 +- 6 files changed, 19 insertions(+), 133 deletions(-) delete mode 100644 TODO.org diff --git a/TODO.org b/TODO.org deleted file mode 100644 index 870855d2..00000000 --- a/TODO.org +++ /dev/null @@ -1,36 +0,0 @@ -det : dans la fonction d'onde vec -csf : dans la fonctions d'onde vec -conf_det : Nconf x Ndet'(conf) -conf_csf : Nconf x Ncsf'(conf) - -csf'(conf) : dans une conf (N_somo(conf) + N_domo(conf)) -det'(conf) : dans une conf (N_somo(conf) + N_domo(conf)) - -det''(conf) : probleme modele avec N_somo(conf) -csf''(conf) : probleme modele avec N_somo(conf) - -* [0/2] - - [ ] Precalculer dans la base des configurations - W_pq = \sum_{K} - K> configuration \in SOMO-SOMO SOMO-DOMO DOMO-VMO SOMO-VMO - W : matrice mono x mono - - [ ] Fonction conf -> - - [ ] Compute in conf basis - - [ ] Function hij_conf(i,j) -> matrix Ndet_ixM - - - - (1,Ndet) (Ndet,Ndet) (Ndet,1) -> E - - for i=1,Nconf - for j=1,Nconf - (i,Ndet') [(Ndet',Ndet'') (Ndet'',Ndet'') (Ndet'',Ndet')] (Ndet',j) -> E - (i,Ncsf') (Ncsf',Ncsf'') [(Ncsf'',Ndet'') (Ndet'',Ndet'') (Ndet'',Ncsf'')] (Ncsf'',Ncsf') (Ncsf',j) -> E - - < K | E_rs] | Psi> - - C.D^t - - K : toutes les mono - - = \sum_m diff --git a/ocaml/qptypes_generator.ml b/ocaml/qptypes_generator.ml index ce99fc78..a5ac22f2 100644 --- a/ocaml/qptypes_generator.ml +++ b/ocaml/qptypes_generator.ml @@ -220,16 +220,16 @@ end = struct type t = | EN | HF - | SOP + | CFG [@@deriving sexp] let to_string = function | EN -> \"EN\" | HF -> \"HF\" - | SOP -> \"SOP\" + | CFG -> \"CFG\" let of_string s = match (String.lowercase_ascii s) with - | \"sop\" -> SOP + | \"cfg\" -> CFG | \"en\" -> EN | \"hf\" -> HF | _ -> raise (Invalid_argument (\"Wrong Perturbation type : \"^s)) diff --git a/src/ao_two_e_ints/map_integrals.irp.f b/src/ao_two_e_ints/map_integrals.irp.f index c0ec9695..55b2d5e2 100644 --- a/src/ao_two_e_ints/map_integrals.irp.f +++ b/src/ao_two_e_ints/map_integrals.irp.f @@ -698,81 +698,3 @@ subroutine insert_into_ao_integrals_map(n_integrals,buffer_i, buffer_values) end -!subroutine dump_ao_integrals(filename) -! use map_module -! implicit none -! BEGIN_DOC -! ! Save to disk the |AO| integrals -! END_DOC -! character*(*), intent(in) :: filename -! integer(cache_key_kind), pointer :: key(:) -! real(integral_kind), pointer :: val(:) -! integer*8 :: i,j, n -! if (.not.mpi_master) then -! return -! endif -! call ezfio_set_work_empty(.False.) -! open(unit=66,file=filename,FORM='unformatted') -! write(66) integral_kind, key_kind -! write(66) ao_integrals_map%sorted, ao_integrals_map%map_size, & -! ao_integrals_map%n_elements -! do i=0_8,ao_integrals_map%map_size -! write(66) ao_integrals_map%map(i)%sorted, ao_integrals_map%map(i)%map_size,& -! ao_integrals_map%map(i)%n_elements -! enddo -! do i=0_8,ao_integrals_map%map_size -! key => ao_integrals_map%map(i)%key -! val => ao_integrals_map%map(i)%value -! n = ao_integrals_map%map(i)%n_elements -! write(66) (key(j), j=1,n), (val(j), j=1,n) -! enddo -! close(66) -! -!end - - -!integer function load_ao_integrals(filename) -! implicit none -! BEGIN_DOC -! ! Read from disk the |AO| integrals -! END_DOC -! character*(*), intent(in) :: filename -! integer*8 :: i -! integer(cache_key_kind), pointer :: key(:) -! real(integral_kind), pointer :: val(:) -! integer :: iknd, kknd -! integer*8 :: n, j -! load_ao_integrals = 1 -! open(unit=66,file=filename,FORM='unformatted',STATUS='UNKNOWN') -! read(66,err=98,end=98) iknd, kknd -! if (iknd /= integral_kind) then -! print *, 'Wrong integrals kind in file :', iknd -! stop 1 -! endif -! if (kknd /= key_kind) then -! print *, 'Wrong key kind in file :', kknd -! stop 1 -! endif -! read(66,err=98,end=98) ao_integrals_map%sorted, ao_integrals_map%map_size,& -! ao_integrals_map%n_elements -! do i=0_8, ao_integrals_map%map_size -! read(66,err=99,end=99) ao_integrals_map%map(i)%sorted, & -! ao_integrals_map%map(i)%map_size, ao_integrals_map%map(i)%n_elements -! call cache_map_reallocate(ao_integrals_map%map(i),ao_integrals_map%map(i)%map_size) -! enddo -! do i=0_8, ao_integrals_map%map_size -! key => ao_integrals_map%map(i)%key -! val => ao_integrals_map%map(i)%value -! n = ao_integrals_map%map(i)%n_elements -! read(66,err=99,end=99) (key(j), j=1,n), (val(j), j=1,n) -! enddo -! call map_sort(ao_integrals_map) -! load_ao_integrals = 0 -! return -! 99 continue -! call map_deinit(ao_integrals_map) -! 98 continue -! stop 'Problem reading ao_integrals_map file in work/' -! -!end -! diff --git a/src/ao_two_e_ints/two_e_integrals.irp.f b/src/ao_two_e_ints/two_e_integrals.irp.f index 8c6b3875..8032bd92 100644 --- a/src/ao_two_e_ints/two_e_integrals.irp.f +++ b/src/ao_two_e_ints/two_e_integrals.irp.f @@ -1162,17 +1162,17 @@ subroutine compute_ao_integrals_jl(j,l,n_integrals,buffer_i,buffer_value) ! Parallel client for AO integrals END_DOC - integer, intent(in) :: j,l - integer,intent(out) :: n_integrals - integer(key_kind),intent(out) :: buffer_i(ao_num*ao_num) + integer, intent(in) :: j,l + integer,intent(out) :: n_integrals + integer(key_kind),intent(out) :: buffer_i(ao_num*ao_num) real(integral_kind),intent(out) :: buffer_value(ao_num*ao_num) + logical, external :: ao_two_e_integral_zero - integer :: i,k - double precision :: ao_two_e_integral,cpu_1,cpu_2, wall_1, wall_2 - double precision :: integral, wall_0 - double precision :: thr - integer :: kk, m, j1, i1 - logical, external :: ao_two_e_integral_zero + integer :: i,k + double precision :: ao_two_e_integral,cpu_1,cpu_2, wall_1, wall_2 + double precision :: integral, wall_0 + double precision :: thr + integer :: kk, m, j1, i1 thr = ao_integrals_threshold diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f index b8cba9a5..03ed91cd 100644 --- a/src/cipsi/selection.irp.f +++ b/src/cipsi/selection.irp.f @@ -730,6 +730,12 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d if(bannedOrb(p2, s2)) cycle if(banned(p1,p2)) cycle + if(pseudo_sym)then + if(dabs(mat(1, p1, p2)).lt.thresh_sym)then + w = 0.d0 + endif + endif + val = maxval(abs(mat(1:N_states, p1, p2))) if( val == 0d0) cycle call apply_particles(mask, s1, p1, s2, p2, det, ok, N_int) @@ -918,12 +924,6 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d end do - if(pseudo_sym)then - if(dabs(mat(1, p1, p2)).lt.thresh_sym)then - w = 0.d0 - endif - endif - ! w = dble(n) * w diff --git a/src/determinants/h_apply.template.f b/src/determinants/h_apply.template.f index abdd6862..a6f9dadd 100644 --- a/src/determinants/h_apply.template.f +++ b/src/determinants/h_apply.template.f @@ -243,7 +243,7 @@ subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl ! Build array of the non-zero integrals of second excitation $filter_integrals - if (ispin == 1) then + if (ispin == 1) then integer :: jjj i=0 From 46767a43a38b75808b1651234f223bde04ded1b9 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 15 Dec 2020 23:51:55 +0100 Subject: [PATCH 2/9] Added AUTHORS file --- AUTHORS | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) create mode 100644 AUTHORS diff --git a/AUTHORS b/AUTHORS new file mode 100644 index 00000000..4e7cce4b --- /dev/null +++ b/AUTHORS @@ -0,0 +1,20 @@ +# If you contributed to this software, please make a pull request to add your +# name to this list (alphabetical order of the last name) + +- Thomas Applencourt +- Anouar Benali +- Michel Caffarel +- Grégoire David +- Anthony Ferté +- Yann Garniron +- Kevin Gasperich +- Vijay Gopal Chilkuri +- Emmanuel Giner +- Pierre-François Loos +- Jean-Paul Malrieu +- Julien Paquier +- Barthélémy Pradines +- Peter Reinhardt +- Anthony Scemama +- Julien Toulouse +- Mickaël Véril From b3e9c49514a01f45a48719aaf474cf6ff22c67ee Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 15 Dec 2020 23:52:17 +0100 Subject: [PATCH 3/9] Added RELEASE_NOTES.org --- RELEASE_NOTES.org | 81 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 81 insertions(+) create mode 100644 RELEASE_NOTES.org diff --git a/RELEASE_NOTES.org b/RELEASE_NOTES.org new file mode 100644 index 00000000..6e46d868 --- /dev/null +++ b/RELEASE_NOTES.org @@ -0,0 +1,81 @@ +#+TITLE: Quantum Package Release notes + +* Version 2.2 + +** New features + +** Changes + + - Python3 replaces Python2 + - Travis CI uses 3 jobs + - Moved Travis scripts into ~travis~ directory + - IRPF90 and EZFIO are now git submodules + - Now basis sets should be downloaded from basis-set-exchange website + - Added ~bse~ in the installable tools + - Documentation in ~src/README.rst~ + - Added two-body reduced density matrix + - Added basis set correction + - Added CAS-based on-top density functional + - Improve PT2 computation for excited-states: Mostly 2x2 + diagonalization, and some (n+1)x(n+1) diagonalizations + - Error bars for stochastic variance and norm of the perturbed wave function + - Improve PT2-matching for excited-states + - Compute the overlap of PT2 excited states + - Renamed SOP into CFG + - Improved parallelism in PT2 by splitting tasks + - Use max in multi-state PT2 instead of sum for the selection weight + - Added seniority + - Added excitation_max + - More tasks for distribueted Davidson + - Random guess vectors in Davidson have zeros to preserve symmetry + - Disk-based Davidson when too much memory is required + - Fixed bug in DIIS + - Fixed bug in molden (Au -> Angs) + +*** User interface + + - Added ~qp_basis~ script to install a basis set from the ~bse~ + command-line tool + - Introduced ~n_det_qp_edit~, ~psi_det_qp_edit~, and + ~psi_coef_qp_edit~ to accelerate the opening of qp_edit with + large wave functions + - Removed ~etc/ninja.rc~ + - Added flag to specify if the AOs are normalized + - Added flag to specify if the primitive Gaussians are normalized + - Added ~lin_dep_cutoff~, the cutoff for linear dependencies + - Davidson convergence threshold can be adapted from PT2 + - In ~density_for_dft~, ~no_core_density~ is now a logical + - Default for ~weight_selection~ has changed from 2 to 1 + - Nullify_small_elements in matrices to keep symmetry + - Default of density functional changed from LDA to PBE + - Added ~no_vvvv_integrals~ flag + - Added ~pt2_min_parallel_tasks~ to control parallelism in PT2 + - Added ~print_energy~ + - Added ~print_hamiltonian~ + - Added input for two body RDM + +*** Code + + - Many bug fixes + - Changed electron-nucleus from ~e_n~ to ~n_e~ in names of variables + - Changed ~occ_pattern~ to ~configuration~ + - Replaced ~List.map~ by a tail-recursive version ~Qputils.list_map~ + - Added possible imaginary part in OCaml MO coefficients + - Added ~qp_clean_source_files.sh~ to remove non-ascii characters + - Added flag ~is_periodic~ for periodic systems + - Possibilities to handle complex integrals and complex MOs + - Moved pseuodpotential integrals out of ~ao_one_e_integrals~ + - Removed Schwarz test and added logical functions + ~ao_two_e_integral_zero~ and ~ao_one_e_integral_zero~ + - Introduced type for ~pt2_data~ + - Banned excitations are used with far apart localized MOs + - S_z2_Sz is now included in S2 + - S^2 in single precision + - Added Shank function + - Added utilities for periodic calculations + + ao_one_e_integral_zero + banned_excitations + + + From 23a96f54aceb1397d76b18be279571417b6343de Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 22 Dec 2020 00:27:09 +0100 Subject: [PATCH 4/9] Multi-state Exc_max --- src/cipsi/selection.irp.f | 27 ++++++++++++++++++++------- src/determinants/determinants.irp.f | 20 ++++++++++++++++++++ 2 files changed, 40 insertions(+), 7 deletions(-) diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f index 03ed91cd..45b84a6a 100644 --- a/src/cipsi/selection.irp.f +++ b/src/cipsi/selection.irp.f @@ -773,22 +773,35 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d endif + integer :: degree + logical :: do_cycle if (excitation_max >= 0) then - integer :: degree - call get_excitation_degree(ref_bitmask(1,1),det(1,1),degree,N_int) - if (degree > excitation_max) cycle + do_cycle = .True. + do k=1,N_states + call get_excitation_degree(psi_det(1,1,dominant_det(k)),det(1,1),degree,N_int) + do_cycle = do_cycle .and. (degree > excitation_max) + enddo + if (do_cycle) cycle endif if (excitation_alpha_max >= 0) then - call get_excitation_degree_spin(ref_bitmask(1,1),det(1,1),degree,N_int) - if (degree > excitation_alpha_max) cycle + do_cycle = .True. + do k=1,N_states + call get_excitation_degree_spin(psi_det(1,1,dominant_det(k)),det(1,1),degree,N_int) + do_cycle = do_cycle .and. (degree > excitation_alpha_max) + enddo + if (do_cycle) cycle endif if (excitation_beta_max >= 0) then - call get_excitation_degree_spin(ref_bitmask(1,2),det(1,2),degree,N_int) - if (degree > excitation_beta_max) cycle + do_cycle = .True. + do k=1,N_states + call get_excitation_degree_spin(psi_det(1,2,dominant_det(k)),det(1,2),degree,N_int) + do_cycle = do_cycle .and. (degree > excitation_beta_max) + enddo + if (do_cycle) cycle endif diff --git a/src/determinants/determinants.irp.f b/src/determinants/determinants.irp.f index 5806634f..2a6057de 100644 --- a/src/determinants/determinants.irp.f +++ b/src/determinants/determinants.irp.f @@ -256,6 +256,26 @@ BEGIN_PROVIDER [ double precision, psi_average_norm_contrib, (psi_det_size) ] enddo END_PROVIDER +BEGIN_PROVIDER [ integer, dominant_det, (N_states) ] + implicit none + BEGIN_DOC + ! Determinant with the largest weight, for each state + END_DOC + integer :: i, k + double precision :: wmax, c + do k=1,N_states + wmax = 0.d0 + do i=1,N_det + c = psi_coef(i,k)*psi_coef(i,k) + if (c > wmax) then + dominant_det(k) = i + wmax = c + endif + enddo + enddo + +END_PROVIDER + !==============================================================================! From a982b0d196a3526429b14c80cb16e6daa1ed1ad1 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 8 Dec 2020 18:44:53 +0100 Subject: [PATCH 5/9] Global Replacement of 'occupation pattern' with 'configuration' --- src/cipsi/cipsi.irp.f | 6 +- src/cipsi/energy.irp.f | 2 +- src/cipsi/pert_rdm_providers.irp.f | 6 +- src/cipsi/pt2_stoch_routines.irp.f | 4 +- src/cipsi/selection.irp.f | 10 +- src/cipsi/selection_buffer.irp.f | 10 +- src/cipsi/slave_cipsi.irp.f | 4 +- src/cipsi/stochastic_cipsi.irp.f | 4 +- ...occ_pattern.irp.f => configurations.irp.f} | 156 +++++++++--------- src/determinants/connected_to_ref.irp.f | 4 +- src/determinants/prune_wf.irp.f | 8 +- src/fci/pt2.irp.f | 2 +- src/iterations/print_summary.irp.f | 6 +- src/perturbation/EZFIO.cfg | 2 +- 14 files changed, 112 insertions(+), 112 deletions(-) rename src/determinants/{occ_pattern.irp.f => configurations.irp.f} (66%) diff --git a/src/cipsi/cipsi.irp.f b/src/cipsi/cipsi.irp.f index 34b16ff3..112ebb56 100644 --- a/src/cipsi/cipsi.irp.f +++ b/src/cipsi/cipsi.irp.f @@ -102,7 +102,7 @@ subroutine run_cipsi call write_double(6,correlation_energy_ratio, 'Correlation ratio') call print_summary(psi_energy_with_nucl_rep, & - pt2_data, pt2_data_err, N_det,N_occ_pattern,N_states,psi_s2) + pt2_data, pt2_data_err, N_det,N_configuration,N_states,psi_s2) call save_energy(psi_energy_with_nucl_rep, pt2_data % pt2) @@ -144,13 +144,13 @@ subroutine run_cipsi SOFT_TOUCH threshold_generators endif print *, 'N_det = ', N_det - print *, 'N_sop = ', N_occ_pattern + print *, 'N_cfg = ', N_configuration print *, 'N_states = ', N_states print*, 'correlation_ratio = ', correlation_energy_ratio call save_energy(psi_energy_with_nucl_rep, pt2_data % pt2) call print_summary(psi_energy_with_nucl_rep(1:N_states), & - pt2_data, pt2_data_err, N_det,N_occ_pattern,N_states,psi_s2) + pt2_data, pt2_data_err, N_det,N_configuration,N_states,psi_s2) call save_iterations(psi_energy_with_nucl_rep(1:N_states),pt2_data % rpt2,N_det) call print_extrapolated_energy() endif diff --git a/src/cipsi/energy.irp.f b/src/cipsi/energy.irp.f index 1d8c6bf5..1f7cf122 100644 --- a/src/cipsi/energy.irp.f +++ b/src/cipsi/energy.irp.f @@ -22,7 +22,7 @@ BEGIN_PROVIDER [ double precision, pt2_E0_denominator, (N_states) ] enddo else if (h0_type == "Barycentric") then pt2_E0_denominator(1:N_states) = barycentric_electronic_energy(1:N_states) - else if (h0_type == "SOP") then + else if (h0_type == "CFG") then pt2_E0_denominator(1:N_states) = psi_energy(1:N_states) else print *, h0_type, ' not implemented' diff --git a/src/cipsi/pert_rdm_providers.irp.f b/src/cipsi/pert_rdm_providers.irp.f index caea57b2..eca8decc 100644 --- a/src/cipsi/pert_rdm_providers.irp.f +++ b/src/cipsi/pert_rdm_providers.irp.f @@ -71,9 +71,9 @@ subroutine fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fo call apply_holes(psi_det_generators(1,1,i_generator), s1, h1, s2, h2, mask, ok, N_int) E_shift = 0.d0 - if (h0_type == 'SOP') then - j = det_to_occ_pattern(i_generator) - E_shift = psi_det_Hii(i_generator) - psi_occ_pattern_Hii(j) + if (h0_type == 'CFG') then + j = det_to_configuration(i_generator) + E_shift = psi_det_Hii(i_generator) - psi_configuration_Hii(j) endif do p1=1,mo_num diff --git a/src/cipsi/pt2_stoch_routines.irp.f b/src/cipsi/pt2_stoch_routines.irp.f index 0966039b..428285c2 100644 --- a/src/cipsi/pt2_stoch_routines.irp.f +++ b/src/cipsi/pt2_stoch_routines.irp.f @@ -134,8 +134,8 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in) PROVIDE psi_bilinear_matrix_transp_order psi_selectors_coef_transp psi_det_sorted PROVIDE psi_det_hii selection_weight pseudo_sym - if (h0_type == 'SOP') then - PROVIDE psi_occ_pattern_hii det_to_occ_pattern + if (h0_type == 'CFG') then + PROVIDE psi_configuration_hii det_to_configuration endif if (N_det <= max(4,N_states) .or. pt2_N_teeth < 2) then diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f index 45b84a6a..bddf7095 100644 --- a/src/cipsi/selection.irp.f +++ b/src/cipsi/selection.irp.f @@ -695,9 +695,9 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d call apply_holes(psi_det_generators(1,1,i_generator), s1, h1, s2, h2, mask, ok, N_int) E_shift = 0.d0 - if (h0_type == 'SOP') then - j = det_to_occ_pattern(i_generator) - E_shift = psi_det_Hii(i_generator) - psi_occ_pattern_Hii(j) + if (h0_type == 'CFG') then + j = det_to_configuration(i_generator) + E_shift = psi_det_Hii(i_generator) - psi_configuration_Hii(j) endif do p1=1,mo_num @@ -810,8 +810,8 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d w = 0d0 ! integer(bit_kind) :: occ(N_int,2), n -! call occ_pattern_of_det(det,occ,N_int) -! call occ_pattern_to_dets_size(occ,n,elec_alpha_num,N_int) +! call configuration_of_det(det,occ,N_int) +! call configuration_to_dets_size(occ,n,elec_alpha_num,N_int) e_pert = 0.d0 coef = 0.d0 diff --git a/src/cipsi/selection_buffer.irp.f b/src/cipsi/selection_buffer.irp.f index 17b6e9a9..10132086 100644 --- a/src/cipsi/selection_buffer.irp.f +++ b/src/cipsi/selection_buffer.irp.f @@ -175,7 +175,7 @@ subroutine make_selection_buffer_s2(b) ! Sort integer, allocatable :: iorder(:) integer*8, allocatable :: bit_tmp(:) - integer*8, external :: occ_pattern_search_key + integer*8, external :: configuration_search_key integer(bit_kind), allocatable :: tmp_array(:,:,:) logical, allocatable :: duplicate(:) @@ -193,7 +193,7 @@ subroutine make_selection_buffer_s2(b) o(k,2,i) = iand(b%det(k,1,i), b%det(k,2,i)) enddo iorder(i) = i - bit_tmp(i) = occ_pattern_search_key(o(1,1,i),N_int) + bit_tmp(i) = configuration_search_key(o(1,1,i),N_int) enddo deallocate(b%det) @@ -279,7 +279,7 @@ subroutine make_selection_buffer_s2(b) ! Create determinants n_d = 0 do i=1,n_p - call occ_pattern_to_dets_size(o(1,1,i),sze,elec_alpha_num,N_int) + call configuration_to_dets_size(o(1,1,i),sze,elec_alpha_num,N_int) n_d = n_d + sze if (n_d > b%cur) then ! if (n_d - b%cur > b%cur - n_d + sze) then @@ -295,8 +295,8 @@ subroutine make_selection_buffer_s2(b) k=1 do i=1,n_p n=n_d - call occ_pattern_to_dets_size(o(1,1,i),n,elec_alpha_num,N_int) - call occ_pattern_to_dets(o(1,1,i),b%det(1,1,k),n,elec_alpha_num,N_int) + call configuration_to_dets_size(o(1,1,i),n,elec_alpha_num,N_int) + call configuration_to_dets(o(1,1,i),b%det(1,1,k),n,elec_alpha_num,N_int) do j=k,k+n-1 b%val(j) = val(i) enddo diff --git a/src/cipsi/slave_cipsi.irp.f b/src/cipsi/slave_cipsi.irp.f index 1dc3e784..4c5d8fca 100644 --- a/src/cipsi/slave_cipsi.irp.f +++ b/src/cipsi/slave_cipsi.irp.f @@ -296,8 +296,8 @@ subroutine run_slave_main print *, 'Number of threads', nproc_target endif - if (h0_type == 'SOP') then - PROVIDE det_to_occ_pattern + if (h0_type == 'CFG') then + PROVIDE det_to_configuration endif PROVIDE global_selection_buffer diff --git a/src/cipsi/stochastic_cipsi.irp.f b/src/cipsi/stochastic_cipsi.irp.f index c529795e..bd0c6a80 100644 --- a/src/cipsi/stochastic_cipsi.irp.f +++ b/src/cipsi/stochastic_cipsi.irp.f @@ -92,7 +92,7 @@ subroutine run_stochastic_cipsi call write_double(6,correlation_energy_ratio, 'Correlation ratio') call print_summary(psi_energy_with_nucl_rep, & - pt2_data, pt2_data_err, N_det,N_occ_pattern,N_states,psi_s2) + pt2_data, pt2_data_err, N_det,N_configuration,N_states,psi_s2) call save_energy(psi_energy_with_nucl_rep, pt2_data % pt2) @@ -131,7 +131,7 @@ subroutine run_stochastic_cipsi call save_energy(psi_energy_with_nucl_rep, pt2_data % pt2) call print_summary(psi_energy_with_nucl_rep, & - pt2_data , pt2_data_err, N_det, N_occ_pattern, N_states, psi_s2) + pt2_data , pt2_data_err, N_det, N_configuration, N_states, psi_s2) call save_iterations(psi_energy_with_nucl_rep(1:N_states),pt2_data % rpt2,N_det) call print_extrapolated_energy() endif diff --git a/src/determinants/occ_pattern.irp.f b/src/determinants/configurations.irp.f similarity index 66% rename from src/determinants/occ_pattern.irp.f rename to src/determinants/configurations.irp.f index d4d8c42b..c703a866 100644 --- a/src/determinants/occ_pattern.irp.f +++ b/src/determinants/configurations.irp.f @@ -1,9 +1,9 @@ use bitmasks -subroutine occ_pattern_of_det(d,o,Nint) +subroutine configuration_of_det(d,o,Nint) use bitmasks implicit none BEGIN_DOC - ! Transforms a determinant to an occupation pattern + ! Transforms a determinant to a configuration ! ! occ(:,1) : Single occupations ! @@ -23,11 +23,11 @@ subroutine occ_pattern_of_det(d,o,Nint) end -subroutine occ_pattern_to_dets_size(o,sze,n_alpha,Nint) +subroutine configuration_to_dets_size(o,sze,n_alpha,Nint) use bitmasks implicit none BEGIN_DOC -! Number of possible determinants for a given occ_pattern +! Number of possible determinants for a given configuration END_DOC integer ,intent(in) :: Nint, n_alpha integer(bit_kind),intent(in) :: o(Nint,2) @@ -49,15 +49,15 @@ subroutine occ_pattern_to_dets_size(o,sze,n_alpha,Nint) end -subroutine occ_pattern_to_dets(o,d,sze,n_alpha,Nint) +subroutine configuration_to_dets(o,d,sze,n_alpha,Nint) use bitmasks implicit none BEGIN_DOC - ! Generate all possible determinants for a given occ_pattern + ! Generate all possible determinants for a given configuration ! ! Input : - ! o : occupation pattern : (doubly occupied, singly occupied) - ! sze : Number of produced determinants, computed by `occ_pattern_to_dets_size` + ! o : configuration : (doubly occupied, singly occupied) + ! sze : Number of produced determinants, computed by `configuration_to_dets_size` ! n_alpha : Number of $\alpha$ electrons ! Nint : N_int ! @@ -184,32 +184,32 @@ subroutine occ_pattern_to_dets(o,d,sze,n_alpha,Nint) end - BEGIN_PROVIDER [ integer(bit_kind), psi_occ_pattern, (N_int,2,psi_det_size) ] -&BEGIN_PROVIDER [ integer, N_occ_pattern ] + BEGIN_PROVIDER [ integer(bit_kind), psi_configuration, (N_int,2,psi_det_size) ] +&BEGIN_PROVIDER [ integer, N_configuration ] implicit none BEGIN_DOC - ! Array of the occ_patterns present in the wave function. + ! Array of the configurations present in the wave function. ! - ! psi_occ_pattern(:,1,j) = j-th occ_pattern of the wave function : represents all the single occupations + ! psi_configuration(:,1,j) = j-th configuration of the wave function : represents all the single occupations ! - ! psi_occ_pattern(:,2,j) = j-th occ_pattern of the wave function : represents all the double occupations + ! psi_configuration(:,2,j) = j-th configuration of the wave function : represents all the double occupations ! - ! The occ patterns are sorted by :c:func:`occ_pattern_search_key` + ! The occ patterns are sorted by :c:func:`configuration_search_key` END_DOC integer :: i,j,k ! create do i = 1, N_det do k = 1, N_int - psi_occ_pattern(k,1,i) = ieor(psi_det(k,1,i),psi_det(k,2,i)) - psi_occ_pattern(k,2,i) = iand(psi_det(k,1,i),psi_det(k,2,i)) + psi_configuration(k,1,i) = ieor(psi_det(k,1,i),psi_det(k,2,i)) + psi_configuration(k,2,i) = iand(psi_det(k,1,i),psi_det(k,2,i)) enddo enddo ! Sort integer, allocatable :: iorder(:) integer*8, allocatable :: bit_tmp(:) - integer*8, external :: occ_pattern_search_key + integer*8, external :: configuration_search_key integer(bit_kind), allocatable :: tmp_array(:,:,:) logical,allocatable :: duplicate(:) logical :: dup @@ -219,7 +219,7 @@ end do i=1,N_det iorder(i) = i - bit_tmp(i) = occ_pattern_search_key(psi_occ_pattern(1,1,i),N_int) + bit_tmp(i) = configuration_search_key(psi_configuration(1,1,i),N_int) enddo call i8sort(bit_tmp,iorder,N_det) @@ -230,8 +230,8 @@ end !$OMP DO do i=1,N_det do k=1,N_int - tmp_array(k,1,i) = psi_occ_pattern(k,1,iorder(i)) - tmp_array(k,2,i) = psi_occ_pattern(k,2,iorder(i)) + tmp_array(k,1,i) = psi_configuration(k,1,iorder(i)) + tmp_array(k,2,i) = psi_configuration(k,2,iorder(i)) enddo duplicate(i) = .False. enddo @@ -273,36 +273,36 @@ end !$OMP END PARALLEL ! Copy filtered result - N_occ_pattern=0 + N_configuration=0 do i=1,N_det if (duplicate(i)) then cycle endif - N_occ_pattern += 1 + N_configuration += 1 do k=1,N_int - psi_occ_pattern(k,1,N_occ_pattern) = tmp_array(k,1,i) - psi_occ_pattern(k,2,N_occ_pattern) = tmp_array(k,2,i) + psi_configuration(k,1,N_configuration) = tmp_array(k,1,i) + psi_configuration(k,2,N_configuration) = tmp_array(k,2,i) enddo enddo !- Check ! print *, 'Checking for duplicates in occ pattern' -! do i=1,N_occ_pattern -! do j=i+1,N_occ_pattern +! do i=1,N_configuration +! do j=i+1,N_configuration ! duplicate(1) = .True. ! do k=1,N_int -! if (psi_occ_pattern(k,1,i) /= psi_occ_pattern(k,1,j)) then +! if (psi_configuration(k,1,i) /= psi_configuration(k,1,j)) then ! duplicate(1) = .False. ! exit ! endif -! if (psi_occ_pattern(k,2,i) /= psi_occ_pattern(k,2,j)) then +! if (psi_configuration(k,2,i) /= psi_configuration(k,2,j)) then ! duplicate(1) = .False. ! exit ! endif ! enddo ! if (duplicate(1)) then -! call debug_det(psi_occ_pattern(1,1,i),N_int) -! call debug_det(psi_occ_pattern(1,1,j),N_int) +! call debug_det(psi_configuration(1,1,i),N_int) +! call debug_det(psi_configuration(1,1,j),N_int) ! stop 'DUPLICATE' ! endif ! enddo @@ -313,21 +313,21 @@ end END_PROVIDER -BEGIN_PROVIDER [ integer, det_to_occ_pattern, (N_det) ] +BEGIN_PROVIDER [ integer, det_to_configuration, (N_det) ] implicit none BEGIN_DOC - ! Returns the index of the occupation pattern for each determinant + ! Returns the index of the configuration for each determinant END_DOC integer :: i,j,k,r,l integer*8 :: key integer(bit_kind) :: occ(N_int,2) logical :: found integer*8, allocatable :: bit_tmp(:) - integer*8, external :: occ_pattern_search_key + integer*8, external :: configuration_search_key - allocate(bit_tmp(N_occ_pattern)) - do i=1,N_occ_pattern - bit_tmp(i) = occ_pattern_search_key(psi_occ_pattern(1,1,i),N_int) + allocate(bit_tmp(N_configuration)) + do i=1,N_configuration + bit_tmp(i) = configuration_search_key(psi_configuration(1,1,i),N_int) enddo !$OMP PARALLEL DO DEFAULT(SHARED) & @@ -338,11 +338,11 @@ BEGIN_PROVIDER [ integer, det_to_occ_pattern, (N_det) ] occ(k,2) = iand(psi_det(k,1,i),psi_det(k,2,i)) enddo - key = occ_pattern_search_key(occ,N_int) + key = configuration_search_key(occ,N_int) ! TODO: Binary search l = 1 - r = N_occ_pattern + r = N_configuration ! do while(r-l > 32) ! j = shiftr(r+l,1) ! if (bit_tmp(j) < key) then @@ -354,14 +354,14 @@ BEGIN_PROVIDER [ integer, det_to_occ_pattern, (N_det) ] do j=l,r found = .True. do k=1,N_int - if ( (occ(k,1) /= psi_occ_pattern(k,1,j)) & - .or. (occ(k,2) /= psi_occ_pattern(k,2,j)) ) then + if ( (occ(k,1) /= psi_configuration(k,1,j)) & + .or. (occ(k,2) /= psi_configuration(k,2,j)) ) then found = .False. exit endif enddo if (found) then - det_to_occ_pattern(i) = j + det_to_configuration(i) = j exit endif enddo @@ -376,78 +376,78 @@ BEGIN_PROVIDER [ integer, det_to_occ_pattern, (N_det) ] END_PROVIDER -BEGIN_PROVIDER [ double precision, psi_occ_pattern_Hii, (N_occ_pattern) ] +BEGIN_PROVIDER [ double precision, psi_configuration_Hii, (N_configuration) ] implicit none BEGIN_DOC - ! $\langle I|H|I \rangle$ where $|I\rangle$ is an occupation pattern. + ! $\langle I|H|I \rangle$ where $|I\rangle$ is a configuration. ! This is the minimum $H_{ii}$, where the $|i\rangle$ are the ! determinants of $|I\rangle$. END_DOC integer :: j, i - psi_occ_pattern_Hii(:) = huge(1.d0) + psi_configuration_Hii(:) = huge(1.d0) do i=1,N_det - j = det_to_occ_pattern(i) - psi_occ_pattern_Hii(j) = min(psi_occ_pattern_Hii(j), psi_det_Hii(i)) + j = det_to_configuration(i) + psi_configuration_Hii(j) = min(psi_configuration_Hii(j), psi_det_Hii(i)) enddo END_PROVIDER -BEGIN_PROVIDER [ double precision, weight_occ_pattern, (N_occ_pattern,N_states) ] +BEGIN_PROVIDER [ double precision, weight_configuration, (N_configuration,N_states) ] implicit none BEGIN_DOC - ! Weight of the occupation patterns in the wave function + ! Weight of the configurations in the wave function END_DOC integer :: i,j,k - weight_occ_pattern = 0.d0 + weight_configuration = 0.d0 do i=1,N_det - j = det_to_occ_pattern(i) + j = det_to_configuration(i) do k=1,N_states - weight_occ_pattern(j,k) += psi_coef(i,k) * psi_coef(i,k) + weight_configuration(j,k) += psi_coef(i,k) * psi_coef(i,k) enddo enddo END_PROVIDER -BEGIN_PROVIDER [ double precision, weight_occ_pattern_average, (N_occ_pattern) ] +BEGIN_PROVIDER [ double precision, weight_configuration_average, (N_configuration) ] implicit none BEGIN_DOC - ! State-average weight of the occupation patterns in the wave function + ! State-average weight of the configurations in the wave function END_DOC integer :: i,j,k - weight_occ_pattern_average(:) = 0.d0 + weight_configuration_average(:) = 0.d0 do i=1,N_det - j = det_to_occ_pattern(i) + j = det_to_configuration(i) do k=1,N_states - weight_occ_pattern_average(j) += psi_coef(i,k) * psi_coef(i,k) * state_average_weight(k) + weight_configuration_average(j) += psi_coef(i,k) * psi_coef(i,k) * state_average_weight(k) enddo enddo END_PROVIDER - BEGIN_PROVIDER [ integer(bit_kind), psi_occ_pattern_sorted, (N_int,2,N_occ_pattern) ] -&BEGIN_PROVIDER [ double precision, weight_occ_pattern_average_sorted, (N_occ_pattern) ] -&BEGIN_PROVIDER [ integer, psi_occ_pattern_sorted_order, (N_occ_pattern) ] -&BEGIN_PROVIDER [ integer, psi_occ_pattern_sorted_order_reverse, (N_occ_pattern) ] + BEGIN_PROVIDER [ integer(bit_kind), psi_configuration_sorted, (N_int,2,N_configuration) ] +&BEGIN_PROVIDER [ double precision, weight_configuration_average_sorted, (N_configuration) ] +&BEGIN_PROVIDER [ integer, psi_configuration_sorted_order, (N_configuration) ] +&BEGIN_PROVIDER [ integer, psi_configuration_sorted_order_reverse, (N_configuration) ] implicit none BEGIN_DOC - ! Occupation patterns sorted by weight + ! Configurations sorted by weight END_DOC integer :: i,j,k integer, allocatable :: iorder(:) - allocate ( iorder(N_occ_pattern) ) - do i=1,N_occ_pattern - weight_occ_pattern_average_sorted(i) = -weight_occ_pattern_average(i) + allocate ( iorder(N_configuration) ) + do i=1,N_configuration + weight_configuration_average_sorted(i) = -weight_configuration_average(i) iorder(i) = i enddo - call dsort(weight_occ_pattern_average_sorted,iorder,N_occ_pattern) - do i=1,N_occ_pattern + call dsort(weight_configuration_average_sorted,iorder,N_configuration) + do i=1,N_configuration do j=1,N_int - psi_occ_pattern_sorted(j,1,i) = psi_occ_pattern(j,1,iorder(i)) - psi_occ_pattern_sorted(j,2,i) = psi_occ_pattern(j,2,iorder(i)) + psi_configuration_sorted(j,1,i) = psi_configuration(j,1,iorder(i)) + psi_configuration_sorted(j,2,i) = psi_configuration(j,2,iorder(i)) enddo - psi_occ_pattern_sorted_order(iorder(i)) = i - psi_occ_pattern_sorted_order_reverse(i) = iorder(i) - weight_occ_pattern_average_sorted(i) = -weight_occ_pattern_average_sorted(i) + psi_configuration_sorted_order(iorder(i)) = i + psi_configuration_sorted_order_reverse(i) = iorder(i) + weight_configuration_average_sorted(i) = -weight_configuration_average_sorted(i) enddo deallocate(iorder) @@ -466,27 +466,27 @@ subroutine make_s2_eigenfunction logical :: update update=.False. - call write_int(6,N_occ_pattern,'Number of occupation patterns') + call write_int(6,N_configuration,'Number of configurations') !$OMP PARALLEL DEFAULT(NONE) & - !$OMP SHARED(N_occ_pattern, psi_occ_pattern, elec_alpha_num,N_int,update) & + !$OMP SHARED(N_configuration, psi_configuration, elec_alpha_num,N_int,update) & !$OMP PRIVATE(s,ithread, d, det_buffer, smax, N_det_new,i,j,k) N_det_new = 0 - call occ_pattern_to_dets_size(psi_occ_pattern(1,1,1),s,elec_alpha_num,N_int) + call configuration_to_dets_size(psi_configuration(1,1,1),s,elec_alpha_num,N_int) allocate (d(N_int,2,s+64), det_buffer(N_int,2,bufsze) ) smax = s ithread=0 !$ ithread = omp_get_thread_num() !$OMP DO SCHEDULE (dynamic,1000) - do i=1,N_occ_pattern - call occ_pattern_to_dets_size(psi_occ_pattern(1,1,i),s,elec_alpha_num,N_int) + do i=1,N_configuration + call configuration_to_dets_size(psi_configuration(1,1,i),s,elec_alpha_num,N_int) s += 1 if (s > smax) then deallocate(d) allocate ( d(N_int,2,s+64) ) smax = s endif - call occ_pattern_to_dets(psi_occ_pattern(1,1,i),d,s,elec_alpha_num,N_int) + call configuration_to_dets(psi_configuration(1,1,i),d,s,elec_alpha_num,N_int) do j=1,s if ( is_in_wavefunction(d(1,1,j), N_int) ) then cycle @@ -511,7 +511,7 @@ subroutine make_s2_eigenfunction if (update) then call copy_H_apply_buffer_to_wf - TOUCH N_det psi_coef psi_det psi_occ_pattern N_occ_pattern + TOUCH N_det psi_coef psi_det psi_configuration N_configuration endif call write_time(6) diff --git a/src/determinants/connected_to_ref.irp.f b/src/determinants/connected_to_ref.irp.f index 086f19d1..a174659c 100644 --- a/src/determinants/connected_to_ref.irp.f +++ b/src/determinants/connected_to_ref.irp.f @@ -12,7 +12,7 @@ integer*8 function det_search_key(det,Nint) end -integer*8 function occ_pattern_search_key(det,Nint) +integer*8 function configuration_search_key(det,Nint) use bitmasks implicit none BEGIN_DOC @@ -22,7 +22,7 @@ integer*8 function occ_pattern_search_key(det,Nint) integer(bit_kind), intent(in) :: det(Nint,2) integer :: i i = shiftr(elec_alpha_num, bit_kind_shift)+1 - occ_pattern_search_key = int(shiftr(ior(det(i,1),det(i,2)),1)+sum(det),8) + configuration_search_key = int(shiftr(ior(det(i,1),det(i,2)),1)+sum(det),8) end diff --git a/src/determinants/prune_wf.irp.f b/src/determinants/prune_wf.irp.f index 7399945f..0bf4b246 100644 --- a/src/determinants/prune_wf.irp.f +++ b/src/determinants/prune_wf.irp.f @@ -10,16 +10,16 @@ BEGIN_PROVIDER [ logical, pruned, (N_det) ] return endif - integer :: i,j,k,ndet_new,nsop_max + integer :: i,j,k,ndet_new,ncfg_max double precision :: thr if (s2_eig) then - nsop_max = max(1,int ( dble(N_occ_pattern) * (1.d0 - pruning) + 0.5d0 )) + ncfg_max = max(1,int ( dble(N_configuration) * (1.d0 - pruning) + 0.5d0 )) do i=1,N_det - k = det_to_occ_pattern(i) - pruned(i) = psi_occ_pattern_sorted_order_reverse(k) > nsop_max + k = det_to_configuration(i) + pruned(i) = psi_configuration_sorted_order_reverse(k) > ncfg_max enddo else diff --git a/src/fci/pt2.irp.f b/src/fci/pt2.irp.f index eb11e28c..49112ec6 100644 --- a/src/fci/pt2.irp.f +++ b/src/fci/pt2.irp.f @@ -54,7 +54,7 @@ subroutine run endif call print_summary(psi_energy_with_nucl_rep(1:N_states), & - pt2_data, pt2_data_err, N_det,N_occ_pattern,N_states,psi_s2) + pt2_data, pt2_data_err, N_det,N_configuration,N_states,psi_s2) call save_energy(E_CI_before, pt2_data % pt2) call pt2_dealloc(pt2_data) diff --git a/src/iterations/print_summary.irp.f b/src/iterations/print_summary.irp.f index d04d8a93..4c7f5d48 100644 --- a/src/iterations/print_summary.irp.f +++ b/src/iterations/print_summary.irp.f @@ -1,11 +1,11 @@ -subroutine print_summary(e_,pt2_data,pt2_data_err,n_det_,n_occ_pattern_,n_st,s2_) +subroutine print_summary(e_,pt2_data,pt2_data_err,n_det_,n_configuration_,n_st,s2_) use selection_types implicit none BEGIN_DOC ! Print the extrapolated energy in the output END_DOC - integer, intent(in) :: n_det_, n_occ_pattern_, n_st + integer, intent(in) :: n_det_, n_configuration_, n_st double precision, intent(in) :: e_(n_st), s2_(n_st) type(pt2_type) , intent(in) :: pt2_data, pt2_data_err integer :: i, k @@ -57,7 +57,7 @@ subroutine print_summary(e_,pt2_data,pt2_data_err,n_det_,n_occ_pattern_,n_st,s2_ print *, 'N_det = ', N_det_ print *, 'N_states = ', n_st if (s2_eig) then - print *, 'N_sop = ', N_occ_pattern_ + print *, 'N_cfg = ', N_configuration_ endif print *, '' diff --git a/src/perturbation/EZFIO.cfg b/src/perturbation/EZFIO.cfg index f4bd8e65..5571f3bc 100644 --- a/src/perturbation/EZFIO.cfg +++ b/src/perturbation/EZFIO.cfg @@ -36,6 +36,6 @@ default: 1.00 [h0_type] type: character*(32) -doc: Type of denominator in PT2. [EN | SOP | HF] +doc: Type of denominator in PT2. [EN | CFG | HF] interface: ezfio,provider,ocaml default: EN From 239ba03231bf2680ea177f0417b21eeeec19f82e Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 22 Dec 2020 01:36:04 +0100 Subject: [PATCH 6/9] Fixed multi-state exc_max --- src/cipsi/selection.irp.f | 15 ++++----- src/determinants/configurations.irp.f | 44 +++++++++++++++++++++++++++ 2 files changed, 52 insertions(+), 7 deletions(-) diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f index bddf7095..12cd623a 100644 --- a/src/cipsi/selection.irp.f +++ b/src/cipsi/selection.irp.f @@ -679,6 +679,7 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d double precision, external :: diag_H_mat_elem_fock double precision :: E_shift double precision :: s_weight(N_states,N_states) + PROVIDE dominant_dets_of_cfgs N_dominant_dets_of_cfgs do jstate=1,N_states do istate=1,N_states s_weight(istate,jstate) = dsqrt(selection_weight(istate)*selection_weight(jstate)) @@ -777,9 +778,9 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d logical :: do_cycle if (excitation_max >= 0) then do_cycle = .True. - do k=1,N_states - call get_excitation_degree(psi_det(1,1,dominant_det(k)),det(1,1),degree,N_int) - do_cycle = do_cycle .and. (degree > excitation_max) + do k=1,N_dominant_dets_of_cfgs + call get_excitation_degree(dominant_dets_of_cfgs(1,1,k),det(1,1),degree,N_int) + do_cycle = do_cycle .and. (degree > excitation_max) enddo if (do_cycle) cycle endif @@ -787,8 +788,8 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d if (excitation_alpha_max >= 0) then do_cycle = .True. - do k=1,N_states - call get_excitation_degree_spin(psi_det(1,1,dominant_det(k)),det(1,1),degree,N_int) + do k=1,N_dominant_dets_of_cfgs + call get_excitation_degree(dominant_dets_of_cfgs(1,1,k),det(1,1),degree,N_int) do_cycle = do_cycle .and. (degree > excitation_alpha_max) enddo if (do_cycle) cycle @@ -797,8 +798,8 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d if (excitation_beta_max >= 0) then do_cycle = .True. - do k=1,N_states - call get_excitation_degree_spin(psi_det(1,2,dominant_det(k)),det(1,2),degree,N_int) + do k=1,N_dominant_dets_of_cfgs + call get_excitation_degree(dominant_dets_of_cfgs(1,1,k),det(1,1),degree,N_int) do_cycle = do_cycle .and. (degree > excitation_beta_max) enddo if (do_cycle) cycle diff --git a/src/determinants/configurations.irp.f b/src/determinants/configurations.irp.f index c703a866..5fb187e1 100644 --- a/src/determinants/configurations.irp.f +++ b/src/determinants/configurations.irp.f @@ -519,3 +519,47 @@ end +BEGIN_PROVIDER [ integer, dominant_cfg, (N_states) ] + implicit none + BEGIN_DOC + ! Configuration of the determinants with the largest weight, for each state + END_DOC + integer :: k + do k=1,N_states + dominant_cfg(k) = det_to_configuration(dominant_det(k)) + enddo +END_PROVIDER + + +BEGIN_PROVIDER [ integer, N_dominant_dets_of_cfgs ] + implicit none + BEGIN_DOC + ! Number of determinants in all dominant determinants + END_DOC + integer :: k, sze + + N_dominant_dets_of_cfgs = 0 + do k=1,N_states + call configuration_to_dets_size( & + psi_configuration(1,1,dominant_cfg(k)), & + sze, elec_alpha_num, N_int) + N_dominant_dets_of_cfgs += sze + enddo +END_PROVIDER + +BEGIN_PROVIDER [ integer(bit_kind), dominant_dets_of_cfgs, (N_int,2,N_dominant_dets_of_cfgs) ] + implicit none + BEGIN_DOC + ! Configuration of the determinants with the largest weight, for each state + END_DOC + integer :: i,k,sze + i=1 + do k=1,N_states + sze = N_dominant_dets_of_cfgs + call configuration_to_dets( & + psi_configuration(1,1,dominant_cfg(k)), & + dominant_dets_of_cfgs(1,1,i), & + sze,elec_alpha_num,N_int) + i += sze + enddo +END_PROVIDER From aeed0f7f3b9740d5b42595f3b6f54edd467119d9 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 23 Dec 2020 02:45:20 +0100 Subject: [PATCH 7/9] Fixed svd --- src/utils/linear_algebra.irp.f | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/src/utils/linear_algebra.irp.f b/src/utils/linear_algebra.irp.f index d858341d..c3b7be79 100644 --- a/src/utils/linear_algebra.irp.f +++ b/src/utils/linear_algebra.irp.f @@ -40,7 +40,7 @@ subroutine svd(A,LDA,U,LDU,D,Vt,LDVt,m,n) stop endif - do j=1,m + do j=1,min(m,n) do i=1,m if (dabs(U(i,j)) < 1.d-14) U(i,j) = 0.d0 enddo @@ -1119,7 +1119,12 @@ subroutine ortho_svd(A,LDA,m,n) double precision, allocatable :: U(:,:), D(:), Vt(:,:) allocate(U(m,n), D(n), Vt(n,n)) call SVD(A,LDA,U,size(U,1),D,Vt,size(Vt,1),m,n) - A(1:m,1:n) = U(1:m,1:n) + integer :: i,j + do j=1,n + do i=1,m + A(i,j) = U(i,j) + enddo + enddo deallocate(U,D, Vt) end From ebafb1b96880b4ffc0f674b78c02e516a7b7a3fe Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 23 Dec 2020 02:45:08 +0100 Subject: [PATCH 8/9] Fixing shiftedbk --- src/dressing/alpha_factory.irp.f | 16 +++++++--------- src/dressing/dress_slave.irp.f | 7 ++++++- 2 files changed, 13 insertions(+), 10 deletions(-) diff --git a/src/dressing/alpha_factory.irp.f b/src/dressing/alpha_factory.irp.f index 92af94d6..5eeeb1a6 100644 --- a/src/dressing/alpha_factory.irp.f +++ b/src/dressing/alpha_factory.irp.f @@ -14,9 +14,7 @@ subroutine alpha_callback(delta_ij_loc, i_generator, subset, csubset, iproc) integer(bit_kind) :: hole_mask(N_int,2), particle_mask(N_int,2) - do l=1,N_generators_bitmask - call generate_singles_and_doubles(delta_ij_loc,i_generator,l,subset,csubset,iproc) - enddo + call generate_singles_and_doubles(delta_ij_loc,i_generator,subset,csubset,iproc) end subroutine @@ -34,7 +32,7 @@ BEGIN_PROVIDER [ integer, psi_from_sorted_gen, (N_det) ] END_PROVIDER -subroutine generate_singles_and_doubles(delta_ij_loc, i_generator, bitmask_index, subset, csubset, iproc) +subroutine generate_singles_and_doubles(delta_ij_loc, i_generator, subset, csubset, iproc) use bitmasks implicit none BEGIN_DOC @@ -42,7 +40,7 @@ subroutine generate_singles_and_doubles(delta_ij_loc, i_generator, bitmask_index END_DOC double precision,intent(inout) :: delta_ij_loc(N_states,N_det,2) - integer, intent(in) :: i_generator, subset, csubset, bitmask_index + integer, intent(in) :: i_generator, subset, csubset integer, intent(in) :: iproc @@ -78,10 +76,10 @@ subroutine generate_singles_and_doubles(delta_ij_loc, i_generator, bitmask_index ! Masks adapted for MRCC do k=1,N_int - hole (k,1) = iand(psi_det_generators(k,1,i_generator), ior(generators_bitmask(k,1,s_hole,bitmask_index),generators_bitmask(k,1,s_part,bitmask_index) ) ) - hole (k,2) = iand(psi_det_generators(k,2,i_generator), ior(generators_bitmask(k,2,s_hole,bitmask_index),generators_bitmask(k,2,s_part,bitmask_index) ) ) - particle(k,1) = iand(not(psi_det_generators(k,1,i_generator)), ior(generators_bitmask(k,1,s_part,bitmask_index),generators_bitmask(k,1,s_hole,bitmask_index)) ) - particle(k,2) = iand(not(psi_det_generators(k,2,i_generator)), ior(generators_bitmask(k,2,s_part,bitmask_index),generators_bitmask(k,2,s_hole,bitmask_index)) ) + hole (k,1) = iand(psi_det_generators(k,1,i_generator), ior(generators_bitmask(k,1,s_hole),generators_bitmask(k,1,s_part) ) ) + hole (k,2) = iand(psi_det_generators(k,2,i_generator), ior(generators_bitmask(k,2,s_hole),generators_bitmask(k,2,s_part) ) ) + particle(k,1) = iand(not(psi_det_generators(k,1,i_generator)), ior(generators_bitmask(k,1,s_part),generators_bitmask(k,1,s_hole)) ) + particle(k,2) = iand(not(psi_det_generators(k,2,i_generator)), ior(generators_bitmask(k,2,s_part),generators_bitmask(k,2,s_hole)) ) enddo integer :: N_holes(2), N_particles(2) diff --git a/src/dressing/dress_slave.irp.f b/src/dressing/dress_slave.irp.f index 7401d0ba..04e4f01b 100644 --- a/src/dressing/dress_slave.irp.f +++ b/src/dressing/dress_slave.irp.f @@ -50,14 +50,19 @@ subroutine run_wf ! Dress ! --------- if (zmq_get_psi(zmq_to_qp_run_socket,1) == -1) cycle + if (zmq_get_dvector(zmq_to_qp_run_socket,1,'threshold_generators',threshold_generators,1) == -1) cycle + if (zmq_get_dvector(zmq_to_qp_run_socket,1,'energy',energy,N_states) == -1) cycle if (zmq_get_N_det_generators (zmq_to_qp_run_socket, 1) == -1) cycle if (zmq_get_N_det_selectors(zmq_to_qp_run_socket, 1) == -1) cycle if (zmq_get_dvector(zmq_to_qp_run_socket,1,'state_average_weight',state_average_weight,N_states) == -1) cycle + if (zmq_get_dvector(zmq_to_qp_run_socket,1,'selection_weight',selection_weight,N_states) == -1) cycle if (zmq_get_dvector(zmq_to_qp_run_socket,1,'energy',energy,N_states) == -1) cycle if (zmq_get_dvector(zmq_to_qp_run_socket,1,'dress_stoch_istate',tmp,1) == -1) cycle dress_stoch_istate = int(tmp) + pt2_e0_denominator(1:N_states) = energy(1:N_states) psi_energy(1:N_states) = energy(1:N_states) - TOUCH psi_energy dress_stoch_istate state_average_weight + TOUCH pt2_e0_denominator dress_stoch_istate state_average_weight threshold_generators selection_weight + PROVIDE psi_bilinear_matrix_columns_loc psi_det_alpha_unique psi_det_beta_unique PROVIDE psi_bilinear_matrix_rows psi_det_sorted_gen_order psi_bilinear_matrix_order From 6d33e6ce81eb658a389fa4b504c049c4e59f9069 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 23 Dec 2020 02:42:38 +0100 Subject: [PATCH 9/9] Moved selection_weight --- src/cipsi/selection.irp.f | 143 ------------------------ src/cipsi/selection_weight.irp.f | 150 ++++++++++++++++++++++++++ src/davidson/davidson_parallel.irp.f | 2 +- src/determinants/configurations.irp.f | 2 +- src/dressing/selection_weight.irp.f | 150 ++++++++++++++++++++++++++ src/utils/linear_algebra.irp.f | 6 +- 6 files changed, 307 insertions(+), 146 deletions(-) create mode 100644 src/cipsi/selection_weight.irp.f create mode 100644 src/dressing/selection_weight.irp.f diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f index 12cd623a..0edd3f63 100644 --- a/src/cipsi/selection.irp.f +++ b/src/cipsi/selection.irp.f @@ -1,148 +1,5 @@ use bitmasks -BEGIN_PROVIDER [ double precision, pt2_match_weight, (N_states) ] - implicit none - BEGIN_DOC - ! Weights adjusted along the selection to make the PT2 contributions - ! of each state coincide. - END_DOC - pt2_match_weight(:) = 1.d0 -END_PROVIDER - -BEGIN_PROVIDER [ double precision, variance_match_weight, (N_states) ] - implicit none - BEGIN_DOC - ! Weights adjusted along the selection to make the variances - ! of each state coincide. - END_DOC - variance_match_weight(:) = 1.d0 -END_PROVIDER - -subroutine update_pt2_and_variance_weights(pt2_data, N_st) - implicit none - use selection_types - BEGIN_DOC -! Updates the PT2- and Variance- matching weights. - END_DOC - integer, intent(in) :: N_st - type(pt2_type), intent(in) :: pt2_data - double precision :: pt2(N_st) - double precision :: variance(N_st) - - double precision :: avg, element, dt, x - integer :: k - integer, save :: i_iter=0 - integer, parameter :: i_itermax = 1 - double precision, allocatable, save :: memo_variance(:,:), memo_pt2(:,:) - - pt2(:) = pt2_data % pt2(:) - variance(:) = pt2_data % variance(:) - - if (i_iter == 0) then - allocate(memo_variance(N_st,i_itermax), memo_pt2(N_st,i_itermax)) - memo_pt2(:,:) = 1.d0 - memo_variance(:,:) = 1.d0 - endif - - i_iter = i_iter+1 - if (i_iter > i_itermax) then - i_iter = 1 - endif - - dt = 2.0d0 - - avg = sum(pt2(1:N_st)) / dble(N_st) - 1.d-32 ! Avoid future division by zero - do k=1,N_st - element = exp(dt*(pt2(k)/avg -1.d0)) - element = min(2.0d0 , element) - element = max(0.5d0 , element) - memo_pt2(k,i_iter) = element - pt2_match_weight(k) *= product(memo_pt2(k,:)) - enddo - - avg = sum(variance(1:N_st)) / dble(N_st) + 1.d-32 ! Avoid future division by zero - do k=1,N_st - element = exp(dt*(variance(k)/avg -1.d0)) - element = min(2.0d0 , element) - element = max(0.5d0 , element) - memo_variance(k,i_iter) = element - variance_match_weight(k) *= product(memo_variance(k,:)) - enddo - - if (N_det < 100) then - ! For tiny wave functions, weights are 1.d0 - pt2_match_weight(:) = 1.d0 - variance_match_weight(:) = 1.d0 - endif - - threshold_davidson_pt2 = min(1.d-6, & - max(threshold_davidson, 1.e-1 * PT2_relative_error * minval(abs(pt2(1:N_states)))) ) - - SOFT_TOUCH pt2_match_weight variance_match_weight threshold_davidson_pt2 -end - - -BEGIN_PROVIDER [ double precision, selection_weight, (N_states) ] - implicit none - BEGIN_DOC - ! Weights used in the selection criterion - END_DOC - select case (weight_selection) - - case (0) - print *, 'Using input weights in selection' - selection_weight(1:N_states) = c0_weight(1:N_states) * state_average_weight(1:N_states) - - case (1) - print *, 'Using 1/c_max^2 weight in selection' - selection_weight(1:N_states) = c0_weight(1:N_states) - - case (2) - print *, 'Using pt2-matching weight in selection' - selection_weight(1:N_states) = c0_weight(1:N_states) * pt2_match_weight(1:N_states) - print *, '# PT2 weight ', real(pt2_match_weight(:),4) - - case (3) - print *, 'Using variance-matching weight in selection' - selection_weight(1:N_states) = c0_weight(1:N_states) * variance_match_weight(1:N_states) - print *, '# var weight ', real(variance_match_weight(:),4) - - case (4) - print *, 'Using variance- and pt2-matching weights in selection' - selection_weight(1:N_states) = c0_weight(1:N_states) * sqrt(variance_match_weight(1:N_states) * pt2_match_weight(1:N_states)) - print *, '# PT2 weight ', real(pt2_match_weight(:),4) - print *, '# var weight ', real(variance_match_weight(:),4) - - case (5) - print *, 'Using variance-matching weight in selection' - selection_weight(1:N_states) = c0_weight(1:N_states) * variance_match_weight(1:N_states) - print *, '# var weight ', real(variance_match_weight(:),4) - - case (6) - print *, 'Using CI coefficient-based selection' - selection_weight(1:N_states) = c0_weight(1:N_states) - - case (7) - print *, 'Input weights multiplied by variance- and pt2-matching' - selection_weight(1:N_states) = c0_weight(1:N_states) * sqrt(variance_match_weight(1:N_states) * pt2_match_weight(1:N_states)) * state_average_weight(1:N_states) - print *, '# PT2 weight ', real(pt2_match_weight(:),4) - print *, '# var weight ', real(variance_match_weight(:),4) - - case (8) - print *, 'Input weights multiplied by pt2-matching' - selection_weight(1:N_states) = c0_weight(1:N_states) * pt2_match_weight(1:N_states) * state_average_weight(1:N_states) - print *, '# PT2 weight ', real(pt2_match_weight(:),4) - - case (9) - print *, 'Input weights multiplied by variance-matching' - selection_weight(1:N_states) = c0_weight(1:N_states) * variance_match_weight(1:N_states) * state_average_weight(1:N_states) - print *, '# var weight ', real(variance_match_weight(:),4) - - end select - print *, '# Total weight ', real(selection_weight(:),4) - -END_PROVIDER - subroutine get_mask_phase(det1, pm, Nint) use bitmasks implicit none diff --git a/src/cipsi/selection_weight.irp.f b/src/cipsi/selection_weight.irp.f new file mode 100644 index 00000000..24e2ace0 --- /dev/null +++ b/src/cipsi/selection_weight.irp.f @@ -0,0 +1,150 @@ +BEGIN_PROVIDER [ double precision, pt2_match_weight, (N_states) ] + implicit none + BEGIN_DOC + ! Weights adjusted along the selection to make the PT2 contributions + ! of each state coincide. + END_DOC + pt2_match_weight(:) = 1.d0 +END_PROVIDER + + + +BEGIN_PROVIDER [ double precision, variance_match_weight, (N_states) ] + implicit none + BEGIN_DOC + ! Weights adjusted along the selection to make the variances + ! of each state coincide. + END_DOC + variance_match_weight(:) = 1.d0 +END_PROVIDER + + + +subroutine update_pt2_and_variance_weights(pt2_data, N_st) + implicit none + use selection_types + BEGIN_DOC +! Updates the PT2- and Variance- matching weights. + END_DOC + integer, intent(in) :: N_st + type(pt2_type), intent(in) :: pt2_data + double precision :: pt2(N_st) + double precision :: variance(N_st) + + double precision :: avg, element, dt, x + integer :: k + integer, save :: i_iter=0 + integer, parameter :: i_itermax = 1 + double precision, allocatable, save :: memo_variance(:,:), memo_pt2(:,:) + + pt2(:) = pt2_data % pt2(:) + variance(:) = pt2_data % variance(:) + + if (i_iter == 0) then + allocate(memo_variance(N_st,i_itermax), memo_pt2(N_st,i_itermax)) + memo_pt2(:,:) = 1.d0 + memo_variance(:,:) = 1.d0 + endif + + i_iter = i_iter+1 + if (i_iter > i_itermax) then + i_iter = 1 + endif + + dt = 2.0d0 + + avg = sum(pt2(1:N_st)) / dble(N_st) - 1.d-32 ! Avoid future division by zero + do k=1,N_st + element = exp(dt*(pt2(k)/avg -1.d0)) + element = min(2.0d0 , element) + element = max(0.5d0 , element) + memo_pt2(k,i_iter) = element + pt2_match_weight(k) *= product(memo_pt2(k,:)) + enddo + + + avg = sum(variance(1:N_st)) / dble(N_st) + 1.d-32 ! Avoid future division by zero + do k=1,N_st + element = exp(dt*(variance(k)/avg -1.d0)) + element = min(2.0d0 , element) + element = max(0.5d0 , element) + memo_variance(k,i_iter) = element + variance_match_weight(k) *= product(memo_variance(k,:)) + enddo + + if (N_det < 100) then + ! For tiny wave functions, weights are 1.d0 + pt2_match_weight(:) = 1.d0 + variance_match_weight(:) = 1.d0 + endif + + threshold_davidson_pt2 = min(1.d-6, & + max(threshold_davidson, 1.e-1 * PT2_relative_error * minval(abs(pt2(1:N_states)))) ) + + SOFT_TOUCH pt2_match_weight variance_match_weight threshold_davidson_pt2 +end + + + + +BEGIN_PROVIDER [ double precision, selection_weight, (N_states) ] + implicit none + BEGIN_DOC + ! Weights used in the selection criterion + END_DOC + select case (weight_selection) + + case (0) + print *, 'Using input weights in selection' + selection_weight(1:N_states) = c0_weight(1:N_states) * state_average_weight(1:N_states) + + case (1) + print *, 'Using 1/c_max^2 weight in selection' + selection_weight(1:N_states) = c0_weight(1:N_states) + + case (2) + print *, 'Using pt2-matching weight in selection' + selection_weight(1:N_states) = c0_weight(1:N_states) * pt2_match_weight(1:N_states) + print *, '# PT2 weight ', real(pt2_match_weight(:),4) + + case (3) + print *, 'Using variance-matching weight in selection' + selection_weight(1:N_states) = c0_weight(1:N_states) * variance_match_weight(1:N_states) + print *, '# var weight ', real(variance_match_weight(:),4) + + case (4) + print *, 'Using variance- and pt2-matching weights in selection' + selection_weight(1:N_states) = c0_weight(1:N_states) * sqrt(variance_match_weight(1:N_states) * pt2_match_weight(1:N_states)) + print *, '# PT2 weight ', real(pt2_match_weight(:),4) + print *, '# var weight ', real(variance_match_weight(:),4) + + case (5) + print *, 'Using variance-matching weight in selection' + selection_weight(1:N_states) = c0_weight(1:N_states) * variance_match_weight(1:N_states) + print *, '# var weight ', real(variance_match_weight(:),4) + + case (6) + print *, 'Using CI coefficient-based selection' + selection_weight(1:N_states) = c0_weight(1:N_states) + + case (7) + print *, 'Input weights multiplied by variance- and pt2-matching' + selection_weight(1:N_states) = c0_weight(1:N_states) * sqrt(variance_match_weight(1:N_states) * pt2_match_weight(1:N_states)) * state_average_weight(1:N_states) + print *, '# PT2 weight ', real(pt2_match_weight(:),4) + print *, '# var weight ', real(variance_match_weight(:),4) + + case (8) + print *, 'Input weights multiplied by pt2-matching' + selection_weight(1:N_states) = c0_weight(1:N_states) * pt2_match_weight(1:N_states) * state_average_weight(1:N_states) + print *, '# PT2 weight ', real(pt2_match_weight(:),4) + + case (9) + print *, 'Input weights multiplied by variance-matching' + selection_weight(1:N_states) = c0_weight(1:N_states) * variance_match_weight(1:N_states) * state_average_weight(1:N_states) + print *, '# var weight ', real(variance_match_weight(:),4) + + end select + print *, '# Total weight ', real(selection_weight(:),4) + +END_PROVIDER + diff --git a/src/davidson/davidson_parallel.irp.f b/src/davidson/davidson_parallel.irp.f index aed81063..59318cd6 100644 --- a/src/davidson/davidson_parallel.irp.f +++ b/src/davidson/davidson_parallel.irp.f @@ -505,7 +505,7 @@ subroutine H_S2_u_0_nstates_zmq(v_0,s_0,u_0,N_st,sze) print *, irp_here, ': Failed in zmq_set_running' endif - call omp_set_nested(.True.) + call omp_set_max_active_levels(4) !$OMP PARALLEL DEFAULT(shared) NUM_THREADS(2) PRIVATE(ithread) ithread = omp_get_thread_num() if (ithread == 0 ) then diff --git a/src/determinants/configurations.irp.f b/src/determinants/configurations.irp.f index 5fb187e1..7e4e6de8 100644 --- a/src/determinants/configurations.irp.f +++ b/src/determinants/configurations.irp.f @@ -68,7 +68,7 @@ subroutine configuration_to_dets(o,d,sze,n_alpha,Nint) integer ,intent(in) :: Nint integer ,intent(in) :: n_alpha ! Number of alpha electrons integer ,intent(inout) :: sze ! Dimension of the output dets - integer(bit_kind),intent(in) :: o(Nint,2) ! Occ patters + integer(bit_kind),intent(in) :: o(Nint,2) ! Configurations integer(bit_kind),intent(out) :: d(Nint,2,sze) ! Output determinants integer :: i, k, n, ispin, ispin2 diff --git a/src/dressing/selection_weight.irp.f b/src/dressing/selection_weight.irp.f new file mode 100644 index 00000000..24e2ace0 --- /dev/null +++ b/src/dressing/selection_weight.irp.f @@ -0,0 +1,150 @@ +BEGIN_PROVIDER [ double precision, pt2_match_weight, (N_states) ] + implicit none + BEGIN_DOC + ! Weights adjusted along the selection to make the PT2 contributions + ! of each state coincide. + END_DOC + pt2_match_weight(:) = 1.d0 +END_PROVIDER + + + +BEGIN_PROVIDER [ double precision, variance_match_weight, (N_states) ] + implicit none + BEGIN_DOC + ! Weights adjusted along the selection to make the variances + ! of each state coincide. + END_DOC + variance_match_weight(:) = 1.d0 +END_PROVIDER + + + +subroutine update_pt2_and_variance_weights(pt2_data, N_st) + implicit none + use selection_types + BEGIN_DOC +! Updates the PT2- and Variance- matching weights. + END_DOC + integer, intent(in) :: N_st + type(pt2_type), intent(in) :: pt2_data + double precision :: pt2(N_st) + double precision :: variance(N_st) + + double precision :: avg, element, dt, x + integer :: k + integer, save :: i_iter=0 + integer, parameter :: i_itermax = 1 + double precision, allocatable, save :: memo_variance(:,:), memo_pt2(:,:) + + pt2(:) = pt2_data % pt2(:) + variance(:) = pt2_data % variance(:) + + if (i_iter == 0) then + allocate(memo_variance(N_st,i_itermax), memo_pt2(N_st,i_itermax)) + memo_pt2(:,:) = 1.d0 + memo_variance(:,:) = 1.d0 + endif + + i_iter = i_iter+1 + if (i_iter > i_itermax) then + i_iter = 1 + endif + + dt = 2.0d0 + + avg = sum(pt2(1:N_st)) / dble(N_st) - 1.d-32 ! Avoid future division by zero + do k=1,N_st + element = exp(dt*(pt2(k)/avg -1.d0)) + element = min(2.0d0 , element) + element = max(0.5d0 , element) + memo_pt2(k,i_iter) = element + pt2_match_weight(k) *= product(memo_pt2(k,:)) + enddo + + + avg = sum(variance(1:N_st)) / dble(N_st) + 1.d-32 ! Avoid future division by zero + do k=1,N_st + element = exp(dt*(variance(k)/avg -1.d0)) + element = min(2.0d0 , element) + element = max(0.5d0 , element) + memo_variance(k,i_iter) = element + variance_match_weight(k) *= product(memo_variance(k,:)) + enddo + + if (N_det < 100) then + ! For tiny wave functions, weights are 1.d0 + pt2_match_weight(:) = 1.d0 + variance_match_weight(:) = 1.d0 + endif + + threshold_davidson_pt2 = min(1.d-6, & + max(threshold_davidson, 1.e-1 * PT2_relative_error * minval(abs(pt2(1:N_states)))) ) + + SOFT_TOUCH pt2_match_weight variance_match_weight threshold_davidson_pt2 +end + + + + +BEGIN_PROVIDER [ double precision, selection_weight, (N_states) ] + implicit none + BEGIN_DOC + ! Weights used in the selection criterion + END_DOC + select case (weight_selection) + + case (0) + print *, 'Using input weights in selection' + selection_weight(1:N_states) = c0_weight(1:N_states) * state_average_weight(1:N_states) + + case (1) + print *, 'Using 1/c_max^2 weight in selection' + selection_weight(1:N_states) = c0_weight(1:N_states) + + case (2) + print *, 'Using pt2-matching weight in selection' + selection_weight(1:N_states) = c0_weight(1:N_states) * pt2_match_weight(1:N_states) + print *, '# PT2 weight ', real(pt2_match_weight(:),4) + + case (3) + print *, 'Using variance-matching weight in selection' + selection_weight(1:N_states) = c0_weight(1:N_states) * variance_match_weight(1:N_states) + print *, '# var weight ', real(variance_match_weight(:),4) + + case (4) + print *, 'Using variance- and pt2-matching weights in selection' + selection_weight(1:N_states) = c0_weight(1:N_states) * sqrt(variance_match_weight(1:N_states) * pt2_match_weight(1:N_states)) + print *, '# PT2 weight ', real(pt2_match_weight(:),4) + print *, '# var weight ', real(variance_match_weight(:),4) + + case (5) + print *, 'Using variance-matching weight in selection' + selection_weight(1:N_states) = c0_weight(1:N_states) * variance_match_weight(1:N_states) + print *, '# var weight ', real(variance_match_weight(:),4) + + case (6) + print *, 'Using CI coefficient-based selection' + selection_weight(1:N_states) = c0_weight(1:N_states) + + case (7) + print *, 'Input weights multiplied by variance- and pt2-matching' + selection_weight(1:N_states) = c0_weight(1:N_states) * sqrt(variance_match_weight(1:N_states) * pt2_match_weight(1:N_states)) * state_average_weight(1:N_states) + print *, '# PT2 weight ', real(pt2_match_weight(:),4) + print *, '# var weight ', real(variance_match_weight(:),4) + + case (8) + print *, 'Input weights multiplied by pt2-matching' + selection_weight(1:N_states) = c0_weight(1:N_states) * pt2_match_weight(1:N_states) * state_average_weight(1:N_states) + print *, '# PT2 weight ', real(pt2_match_weight(:),4) + + case (9) + print *, 'Input weights multiplied by variance-matching' + selection_weight(1:N_states) = c0_weight(1:N_states) * variance_match_weight(1:N_states) * state_average_weight(1:N_states) + print *, '# var weight ', real(variance_match_weight(:),4) + + end select + print *, '# Total weight ', real(selection_weight(:),4) + +END_PROVIDER + diff --git a/src/utils/linear_algebra.irp.f b/src/utils/linear_algebra.irp.f index c3b7be79..ab804983 100644 --- a/src/utils/linear_algebra.irp.f +++ b/src/utils/linear_algebra.irp.f @@ -19,7 +19,11 @@ subroutine svd(A,LDA,U,LDU,D,Vt,LDVt,m,n) double precision,allocatable :: A_tmp(:,:) allocate (A_tmp(LDA,n)) - A_tmp(:,:) = A(:,:) + do k=1,n + do i=1,m + A_tmp(i,k) = A(i,k) + enddo + enddo ! Find optimal size for temp arrays allocate(work(1))