From a1415feab714faa2d40ef9bf667414f051e5268e Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 9 Jan 2018 17:34:15 +0100 Subject: [PATCH 01/14] Bug H_Core guess --- plugins/Hartree_Fock/SCF.irp.f | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/plugins/Hartree_Fock/SCF.irp.f b/plugins/Hartree_Fock/SCF.irp.f index 3d71d826..75a07145 100644 --- a/plugins/Hartree_Fock/SCF.irp.f +++ b/plugins/Hartree_Fock/SCF.irp.f @@ -23,7 +23,7 @@ subroutine create_guess mo_coef = ao_ortho_lowdin_coef TOUCH mo_coef mo_label = 'Guess' - call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral,size(mo_mono_elec_integral,1),size(mo_mono_elec_integral,2),mo_label,.false.) + call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral,size(mo_mono_elec_integral,1),size(mo_mono_elec_integral,2),mo_label,1,.false.) SOFT_TOUCH mo_coef mo_label else if (mo_guess_type == "Huckel") then call huckel_guess From 09ec85f5e9e96cbdee800099f77df577767aa468 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 10 Jan 2018 16:54:12 +0100 Subject: [PATCH 02/14] SCF --- plugins/Hartree_Fock/Roothaan_Hall_SCF.irp.f | 2 ++ 1 file changed, 2 insertions(+) diff --git a/plugins/Hartree_Fock/Roothaan_Hall_SCF.irp.f b/plugins/Hartree_Fock/Roothaan_Hall_SCF.irp.f index 860cc825..c66c8985 100644 --- a/plugins/Hartree_Fock/Roothaan_Hall_SCF.irp.f +++ b/plugins/Hartree_Fock/Roothaan_Hall_SCF.irp.f @@ -14,6 +14,8 @@ END_DOC integer :: i,j double precision, allocatable :: mo_coef_save(:,:) + + PROVIDE ao_md5 mo_occ level_shift allocate(mo_coef_save(ao_num,mo_tot_num), & Fock_matrix_DIIS (ao_num,ao_num,max_dim_DIIS), & From 9103c6cf52aed971d0386b701d8b9f32dd5385fd Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 10 Jan 2018 18:11:49 +0100 Subject: [PATCH 03/14] Fixed multi-state --- plugins/Full_CI_ZMQ/pt2_stoch_routines.irp.f | 15 ++++---- plugins/Generators_full/generators.irp.f | 11 ++---- src/Determinants/density_matrix.irp.f | 6 ++-- src/Determinants/determinants.irp.f | 38 +++----------------- 4 files changed, 18 insertions(+), 52 deletions(-) diff --git a/plugins/Full_CI_ZMQ/pt2_stoch_routines.irp.f b/plugins/Full_CI_ZMQ/pt2_stoch_routines.irp.f index 5848cec0..9d1c50d4 100644 --- a/plugins/Full_CI_ZMQ/pt2_stoch_routines.irp.f +++ b/plugins/Full_CI_ZMQ/pt2_stoch_routines.irp.f @@ -25,8 +25,8 @@ subroutine ZMQ_pt2(E, pt2,relative_error, absolute_error, error) double precision :: sumabove(comb_teeth), sum2above(comb_teeth), Nabove(comb_teeth) double precision, external :: omp_get_wtime + double precision :: state_average_weight_save(N_states), w(N_states) double precision :: time - double precision :: w(N_states) integer(ZMQ_PTR), external :: new_zmq_to_qp_run_socket if (N_det < max(10,N_states)) then @@ -35,18 +35,19 @@ subroutine ZMQ_pt2(E, pt2,relative_error, absolute_error, error) error(:) = 0.d0 else + state_average_weight_save(:) = state_average_weight(:) do pt2_stoch_istate=1,N_states SOFT_TOUCH pt2_stoch_istate - w(:) = 0.d0 - w(pt2_stoch_istate) = 1.d0 - call update_psi_average_norm_contrib(w) + state_average_weight(:) = 0.d0 + state_average_weight(pt2_stoch_istate) = 1.d0 + TOUCH state_average_weight allocate(pt2_detail(N_states,N_det_generators+1), comb(N_det_generators), computed(N_det_generators), tbc(0:size_tbc)) sumabove = 0d0 sum2above = 0d0 Nabove = 0d0 - provide nproc fragment_first fragment_count mo_bielec_integrals_in_map mo_mono_elec_integral pt2_weight psi_selectors + provide nproc fragment_first fragment_count mo_bielec_integrals_in_map mo_mono_elec_integral pt2_weight psi_selectors computed = .false. @@ -141,7 +142,9 @@ subroutine ZMQ_pt2(E, pt2,relative_error, absolute_error, error) deallocate(pt2_detail, comb, computed, tbc) enddo - FREE psi_average_norm_contrib pt2_stoch_istate + FREE pt2_stoch_istate + state_average_weight(:) = state_average_weight_save(:) + TOUCH state_average_weight endif do k=N_det+1,N_states pt2(k) = 0.d0 diff --git a/plugins/Generators_full/generators.irp.f b/plugins/Generators_full/generators.irp.f index a04065cf..4f2c715e 100644 --- a/plugins/Generators_full/generators.irp.f +++ b/plugins/Generators_full/generators.irp.f @@ -30,15 +30,8 @@ END_PROVIDER ! Hartree-Fock determinant END_DOC integer :: i, k - psi_coef_generators = 0.d0 - psi_det_generators = 0_bit_kind - do i=1,N_det_generators - do k=1,N_int - psi_det_generators(k,1,i) = psi_det_sorted(k,1,i) - psi_det_generators(k,2,i) = psi_det_sorted(k,2,i) - enddo - psi_coef_generators(i,:) = psi_coef_sorted(i,:) - enddo + psi_det_generators(1:N_int,1:2,1:N_det) = psi_det_sorted(1:N_int,1:2,1:N_det) + psi_coef_generators(1:N_det,1:N_states) = psi_coef_sorted(1:N_det,1:N_states) END_PROVIDER diff --git a/src/Determinants/density_matrix.irp.f b/src/Determinants/density_matrix.irp.f index be28183b..bd5f0741 100644 --- a/src/Determinants/density_matrix.irp.f +++ b/src/Determinants/density_matrix.irp.f @@ -368,13 +368,13 @@ BEGIN_PROVIDER [ double precision, state_average_weight, (N_states) ] END_DOC logical :: exists - state_average_weight = 1.d0 + state_average_weight(:) = 1.d0 call ezfio_has_determinants_state_average_weight(exists) if (exists) then call ezfio_get_determinants_state_average_weight(state_average_weight) endif - state_average_weight = state_average_weight+1.d-31 - state_average_weight = state_average_weight/(sum(state_average_weight)) + state_average_weight(:) = state_average_weight(:)+1.d-31 + state_average_weight(:) = state_average_weight(:)/(sum(state_average_weight(:))) END_PROVIDER diff --git a/src/Determinants/determinants.irp.f b/src/Determinants/determinants.irp.f index 2ef5dfac..8530fa64 100644 --- a/src/Determinants/determinants.irp.f +++ b/src/Determinants/determinants.irp.f @@ -225,34 +225,6 @@ BEGIN_PROVIDER [ double precision, psi_coef, (psi_det_size,N_states) ] END_PROVIDER -subroutine update_psi_average_norm_contrib(w) - implicit none - BEGIN_DOC - ! Compute psi_average_norm_contrib for different state average weights w(:) - END_DOC - double precision, intent(in) :: w(N_states) - double precision :: w0(N_states), f - w0(:) = w(:)/sum(w(:)) - - integer :: i,j,k - do i=1,N_det - psi_average_norm_contrib(i) = psi_coef(i,1)*psi_coef(i,1)*w(1) - enddo - do k=2,N_states - do i=1,N_det - psi_average_norm_contrib(i) = psi_average_norm_contrib(i) + & - psi_coef(i,k)*psi_coef(i,k)*w(k) - enddo - enddo - f = 1.d0/sum(psi_average_norm_contrib(1:N_det)) - do i=1,N_det - psi_average_norm_contrib(i) = psi_average_norm_contrib(i)*f - enddo - SOFT_TOUCH psi_average_norm_contrib - -end subroutine - - BEGIN_PROVIDER [ double precision, psi_average_norm_contrib, (psi_det_size) ] implicit none BEGIN_DOC @@ -260,14 +232,12 @@ BEGIN_PROVIDER [ double precision, psi_average_norm_contrib, (psi_det_size) ] END_DOC integer :: i,j,k double precision :: f - f = 1.d0/dble(N_states) - do i=1,N_det - psi_average_norm_contrib(i) = psi_coef(i,1)*psi_coef(i,1)*f - enddo - do k=2,N_states + + psi_average_norm_contrib(:) = 0.d0 + do k=1,N_states do i=1,N_det psi_average_norm_contrib(i) = psi_average_norm_contrib(i) + & - psi_coef(i,k)*psi_coef(i,k)*f + psi_coef(i,k)*psi_coef(i,k)*state_average_weight(k) enddo enddo f = 1.d0/sum(psi_average_norm_contrib(1:N_det)) From 2550047bb2378779294cd3e3aad9aa0e1cb40a72 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 10 Jan 2018 19:30:02 +0100 Subject: [PATCH 04/14] Fixed mrcc stoch --- plugins/mrcepa0/mrcc_stoch_routines.irp.f | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/plugins/mrcepa0/mrcc_stoch_routines.irp.f b/plugins/mrcepa0/mrcc_stoch_routines.irp.f index 78940d5e..d68b5137 100644 --- a/plugins/mrcepa0/mrcc_stoch_routines.irp.f +++ b/plugins/mrcepa0/mrcc_stoch_routines.irp.f @@ -32,17 +32,15 @@ subroutine ZMQ_mrcc(E, mrcc, delta, delta_s2, relative_error) double precision, external :: omp_get_wtime double precision :: time - double precision :: w(N_states) + state_average_weight(:) = 0.d0 + state_average_weight(mrcc_stoch_istate) = 1.d0 + TOUCH state_average_weight + provide nproc fragment_first fragment_count mo_bielec_integrals_in_map mo_mono_elec_integral mrcc_weight psi_selectors - w(:) = 0.d0 - w(mrcc_stoch_istate) = 1.d0 - call update_psi_average_norm_contrib(w) - - print *, '========== ================= ================= =================' From 5d0dbb304a4848041f2235ccc0654644e7beaf63 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 11 Jan 2018 10:08:11 +0100 Subject: [PATCH 05/14] Memo --- plugins/Full_CI_ZMQ/selection.irp.f | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/plugins/Full_CI_ZMQ/selection.irp.f b/plugins/Full_CI_ZMQ/selection.irp.f index 81ff5795..378e51c4 100644 --- a/plugins/Full_CI_ZMQ/selection.irp.f +++ b/plugins/Full_CI_ZMQ/selection.irp.f @@ -623,6 +623,9 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d Hii = diag_H_mat_elem_fock(psi_det_generators(1,1,i_generator),det,fock_diag_tmp,N_int) min_e_pert = 0d0 +! double precision :: hij +! call i_h_j(psi_det_generators(1,1,i_generator), det, N_int, hij) + do istate=1,N_states delta_E = E0(istate) - Hii val = mat(istate, p1, p2) + mat(istate, p1, p2) @@ -633,7 +636,7 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d e_pert = 0.5d0 * (tmp - delta_E) pt2(istate) = pt2(istate) + e_pert min_e_pert = min(e_pert,min_e_pert) -! ci(istate) = e_pert / mat(istate, p1, p2) +! ci(istate) = e_pert / hij end do if(min_e_pert <= buf%mini) then From d31cca38df6f9b78ec1b5e1444e2cb32f6fd45a6 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 12 Jan 2018 21:42:37 +0100 Subject: [PATCH 06/14] Fixes #223 --- plugins/GPI2/NEEDED_CHILDREN_MODULES | 1 - plugins/GPI2/README.rst | 14 ---- plugins/GPI2/gpi_test.irp.f | 13 ---- plugins/GPI2/utils.irp.f | 76 ------------------- .../read_integral/print_integrals_mo.irp.f | 6 +- plugins/read_integral/read_integrals_mo.irp.f | 3 +- 6 files changed, 5 insertions(+), 108 deletions(-) delete mode 100644 plugins/GPI2/NEEDED_CHILDREN_MODULES delete mode 100644 plugins/GPI2/README.rst delete mode 100644 plugins/GPI2/gpi_test.irp.f delete mode 100644 plugins/GPI2/utils.irp.f diff --git a/plugins/GPI2/NEEDED_CHILDREN_MODULES b/plugins/GPI2/NEEDED_CHILDREN_MODULES deleted file mode 100644 index aae89501..00000000 --- a/plugins/GPI2/NEEDED_CHILDREN_MODULES +++ /dev/null @@ -1 +0,0 @@ -Determinants diff --git a/plugins/GPI2/README.rst b/plugins/GPI2/README.rst deleted file mode 100644 index d6be4958..00000000 --- a/plugins/GPI2/README.rst +++ /dev/null @@ -1,14 +0,0 @@ -===== -GASPI -===== - -Providers for GASPI programs (with the GPI2 library). - -Needed Modules -============== -.. Do not edit this section It was auto-generated -.. by the `update_README.py` script. -Documentation -============= -.. Do not edit this section It was auto-generated -.. by the `update_README.py` script. diff --git a/plugins/GPI2/gpi_test.irp.f b/plugins/GPI2/gpi_test.irp.f deleted file mode 100644 index 1fc109e6..00000000 --- a/plugins/GPI2/gpi_test.irp.f +++ /dev/null @@ -1,13 +0,0 @@ -program test - double precision :: energy(N_states) - if (is_gaspi_master) then - energy = 1.d0 - else - energy = 0.d0 - endif - call broadcast_wf(energy) - print *, 'energy (1.d0) :', GASPI_rank, energy(1) - print *, 'coef :', GASPI_rank, psi_coef(1,1) - print *, 'det :', GASPI_rank, psi_det (1,1,1) - call gaspi_finalize -end diff --git a/plugins/GPI2/utils.irp.f b/plugins/GPI2/utils.irp.f deleted file mode 100644 index cfb17b75..00000000 --- a/plugins/GPI2/utils.irp.f +++ /dev/null @@ -1,76 +0,0 @@ - BEGIN_PROVIDER [ logical, GASPI_is_initialized ] -&BEGIN_PROVIDER [ logical, has_gaspi ] - implicit none - BEGIN_DOC -! This is true when GASPI_Init has been called - END_DOC - - has_gaspi = .False. - IRP_IF GASPI - use GASPI - integer(gaspi_return_t) :: res - res = gaspi_proc_init(GASPI_BLOCK) - if (res /= GASPI_SUCCESS) then - print *, res - print *, 'GASPI failed to initialize' - stop -1 - endif - has_gaspi = .True. - IRP_ENDIF - GASPI_is_initialized = .True. -END_PROVIDER - - - BEGIN_PROVIDER [ integer, GASPI_rank ] -&BEGIN_PROVIDER [ integer, GASPI_size ] -&BEGIN_PROVIDER [ logical, is_GASPI_master ] - implicit none - BEGIN_DOC -! Usual GASPI variables - END_DOC - - PROVIDE GASPI_is_initialized - - IRP_IF GASPI - use GASPI - integer(gaspi_return_t) :: res - integer(gaspi_rank_t) :: n - res = gaspi_proc_num(n) - GASPI_size = n - if (res /= GASPI_SUCCESS) then - print *, res - print *, 'Unable to get GASPI_size' - stop -1 - endif - res = gaspi_proc_rank(n) - GASPI_rank = n - if (res /= GASPI_SUCCESS) then - print *, res - print *, 'Unable to get GASPI_rank' - stop -1 - endif - is_GASPI_master = (GASPI_rank == 0) - IRP_ELSE - GASPI_rank = 0 - GASPI_size = 1 - is_GASPI_master = .True. - IRP_ENDIF - - -END_PROVIDER - -subroutine gaspi_finalize() - implicit none - PROVIDE GASPI_is_initialized - IRP_IF GASPI - use GASPI - integer(gaspi_return_t) :: res - res = gaspi_proc_term(GASPI_BLOCK) - if (res /= GASPI_SUCCESS) then - print *, res - print *, 'Unable to finalize GASPI' - stop -1 - endif - IRP_ENDIF -end subroutine - diff --git a/plugins/read_integral/print_integrals_mo.irp.f b/plugins/read_integral/print_integrals_mo.irp.f index 18795249..45745c13 100644 --- a/plugins/read_integral/print_integrals_mo.irp.f +++ b/plugins/read_integral/print_integrals_mo.irp.f @@ -45,13 +45,13 @@ program print_integrals do k=1,mo_tot_num do j=l,mo_tot_num do i=k,mo_tot_num - !if (i>=j) then + if (i>=j) then double precision :: get_mo_bielec_integral integral = get_mo_bielec_integral(i,j,k,l,mo_integrals_map) if (dabs(integral) > mo_integrals_threshold) then - write (iunit,'(4(I6,X),F20.15)') i,j,k,l, integral + write (iunit,'(4(I6,X),E25.15)') i,j,k,l, integral endif - !end if + end if enddo enddo enddo diff --git a/plugins/read_integral/read_integrals_mo.irp.f b/plugins/read_integral/read_integrals_mo.irp.f index 06db0ddf..a2d1cb6b 100644 --- a/plugins/read_integral/read_integrals_mo.irp.f +++ b/plugins/read_integral/read_integrals_mo.irp.f @@ -69,9 +69,10 @@ subroutine run 13 continue close(iunit) - call insert_into_mo_integrals_map(n_integrals,buffer_i,buffer_values,0.d0) + call map_append(mo_integrals_map, buffer_i, buffer_values, n_integrals) call map_sort(mo_integrals_map) + call map_unique(mo_integrals_map) call map_save_to_disk(trim(ezfio_filename)//'/work/mo_ints',mo_integrals_map) call ezfio_set_integrals_bielec_disk_access_mo_integrals("Read") From 02161df6d0078ef3c7a40d0433e5bb10c67b0a69 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sat, 13 Jan 2018 18:15:12 +0100 Subject: [PATCH 07/14] Fixes #223 --- plugins/read_integral/NEEDED_CHILDREN_MODULES | 2 +- .../read_integral/print_integrals_mo.irp.f | 14 ++++---- plugins/read_integral/read_integrals_mo.irp.f | 36 +++++++++++++++++++ 3 files changed, 43 insertions(+), 9 deletions(-) diff --git a/plugins/read_integral/NEEDED_CHILDREN_MODULES b/plugins/read_integral/NEEDED_CHILDREN_MODULES index e492a3ce..566762ba 100644 --- a/plugins/read_integral/NEEDED_CHILDREN_MODULES +++ b/plugins/read_integral/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -Integrals_Monoelec Integrals_Bielec +Integrals_Monoelec Integrals_Bielec Hartree_Fock diff --git a/plugins/read_integral/print_integrals_mo.irp.f b/plugins/read_integral/print_integrals_mo.irp.f index 45745c13..2381da52 100644 --- a/plugins/read_integral/print_integrals_mo.irp.f +++ b/plugins/read_integral/print_integrals_mo.irp.f @@ -44,14 +44,12 @@ program print_integrals do l=1,mo_tot_num do k=1,mo_tot_num do j=l,mo_tot_num - do i=k,mo_tot_num - if (i>=j) then - double precision :: get_mo_bielec_integral - integral = get_mo_bielec_integral(i,j,k,l,mo_integrals_map) - if (dabs(integral) > mo_integrals_threshold) then - write (iunit,'(4(I6,X),E25.15)') i,j,k,l, integral - endif - end if + do i=max(j,k),mo_tot_num + double precision :: get_mo_bielec_integral + integral = get_mo_bielec_integral(i,j,k,l,mo_integrals_map) + if (dabs(integral) > mo_integrals_threshold) then + write (iunit,'(4(I6,X),E25.15)') i,j,k,l, integral + endif enddo enddo enddo diff --git a/plugins/read_integral/read_integrals_mo.irp.f b/plugins/read_integral/read_integrals_mo.irp.f index a2d1cb6b..c021941c 100644 --- a/plugins/read_integral/read_integrals_mo.irp.f +++ b/plugins/read_integral/read_integrals_mo.irp.f @@ -5,8 +5,44 @@ program read_integrals ! - nuclear_mo ! - bielec_mo END_DOC + + integer :: iunit + integer :: getunitandopen + integer :: i,j,n + PROVIDE ezfio_filename call ezfio_set_integrals_monoelec_disk_access_mo_one_integrals("None") + + logical :: has + call ezfio_has_mo_basis_mo_tot_num(has) + if (.not.has) then + + iunit = getunitandopen('nuclear_mo','r') + n=0 + do + read (iunit,*,end=12) i + n = max(n,i) + enddo + 12 continue + close(iunit) + call ezfio_set_mo_basis_mo_tot_num(n) + + call ezfio_has_ao_basis_ao_num(has) + mo_label = "None" + if (has) then + call huckel_guess + else + call ezfio_set_ao_basis_ao_num(n) + double precision, allocatable :: X(:,:) + allocate (X(n,n)) + X = 0.d0 + do i=1,n + X(i,i) = 1.d0 + enddo + call ezfio_set_mo_basis_mo_coef(X) + call save_mos + endif + endif call run end From 35708de94457a19ec8111522ac91884e29c9379b Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 22 Jan 2018 01:19:03 +0100 Subject: [PATCH 08/14] Upgrade to OCaml 4.06 and Core 0.10 --- README.md | 6 +++-- configure | 28 +++++++++++++--------- install/scripts/install_ocaml.sh | 5 ++-- ocaml/TaskServer.ml | 11 +++++---- ocaml/qp_create_ezfio_from_xyz.ml | 2 +- ocaml/qp_create_guess.ml | 2 +- ocaml/qp_find_pi_space.ml | 2 +- ocaml/qp_print.ml | 2 +- ocaml/qp_run.ml | 2 +- ocaml/qp_set_mo_class.ml | 2 +- scripts/ezfio_interface/upgrade_1.0_2.0.sh | 27 --------------------- scripts/qp_upgrade_ocaml.sh | 19 +++++++++++++++ 12 files changed, 55 insertions(+), 53 deletions(-) delete mode 100755 scripts/ezfio_interface/upgrade_1.0_2.0.sh create mode 100755 scripts/qp_upgrade_ocaml.sh diff --git a/README.md b/README.md index a11c5713..52f949c3 100644 --- a/README.md +++ b/README.md @@ -1,7 +1,9 @@ ## IMPORTANT -If you have problems upgrading to the current version, consider re-installing everything from scratch including the OCaml compiler. -To do this, you will have to remove the `quantum_package` directory **and** the `$HOME/.opam` directory as well. +If you have problems upgrading to the current version, first try +`qp_upgrade_ocaml.sh`. If it fails, then consider re-installing everything from +scratch including the OCaml compiler. To do this, you will have to remove the +`quantum_package` directory **and** the `$HOME/.opam` directory as well. diff --git a/configure b/configure index 9b59b209..9f677e92 100755 --- a/configure +++ b/configure @@ -49,7 +49,7 @@ QP_ROOT_INSTALL = join(QP_ROOT, "install") os.environ["PATH"] = os.environ["PATH"] + ":" + QP_ROOT_BIN d_dependency = { - "ocaml": ["m4", "curl", "zlib", "patch", "gcc", "zeromq"], + "ocaml": ["m4", "curl", "zlib", "patch", "gcc", "zeromq", "gmp"], "m4": ["make"], "curl": ["make"], "zlib": ["gcc", "make"], @@ -67,7 +67,8 @@ d_dependency = { "ninja": ["g++", "python"], "make": [], "p_graphviz": ["python"], - "bats": [] + "bats": [], + "gmp" : ["make", "g++"] } from collections import namedtuple @@ -136,6 +137,11 @@ zeromq = Info( description=' ZeroMQ', default_path=join(QP_ROOT_LIB, "libzmq.a")) +gmp= Info( + url='https://gmplib.org/download/gmp/gmp-6.1.2.tar.bz2', + description=' The GNU Multiple Precision Arithmetic Library', + default_path=join(QP_ROOT_LIB, "libgmp.a")) + f77zmq = Info( url='{head}/zeromq/f77_zmq/{tail}'.format(**path_github), description=' F77-ZeroMQ', @@ -155,7 +161,7 @@ d_info = dict() for m in ["ocaml", "m4", "curl", "zlib", "patch", "irpf90", "docopt", "resultsFile", "ninja", "emsl", "ezfio", "p_graphviz", - "zeromq", "f77zmq", "bats"]: + "zeromq", "f77zmq", "bats", "gmp"]: exec ("d_info['{0}']={0}".format(m)) @@ -480,16 +486,16 @@ def create_ninja_and_rc(l_installed): 'export QP_PYTHON={0}'.format(":".join(l_python)), "", 'export IRPF90={0}'.format(path_irpf90.replace(QP_ROOT,"${QP_ROOT}")), 'export NINJA={0}'.format(path_ninja.replace(QP_ROOT,"${QP_ROOT}")), - 'function qp_append_export () {', - ' #Append path $2:${!1}. Add the semicolon only if ${!1} is defined', + 'function qp_prepend_export () {', + ' #Prepend path $2:${!1}. Add the semicolon only if ${!1} is defined', ' eval "value_1=\"\${$1}\""', - ' echo ${2}${value_1:+:${value_1}}', + ' echo ${value_1:+${value_1}:}${2}', '}', - 'export PYTHONPATH=$(qp_append_export "PYTHONPATH" "${QP_EZFIO}/Python":"${QP_PYTHON}")', - 'export PATH=$(qp_append_export "PATH" "${QP_PYTHON}":"${QP_ROOT}"/bin:"${QP_ROOT}"/ocaml)', - 'export LD_LIBRARY_PATH=$(qp_append_export "LD_LIBRARY_PATH" "${QP_ROOT}"/lib:"${QP_ROOT}"/lib64)', - 'export LIBRARY_PATH=$(qp_append_export "LIBRARY_PATH" "${QP_ROOT}"/lib:"${QP_ROOT}"/lib64)', - 'export C_INCLUDE_PATH=$(qp_append_export "C_INCLUDE_PATH" "${QP_ROOT}"/include)', + 'export PYTHONPATH=$(qp_prepend_export "PYTHONPATH" "${QP_EZFIO}/Python":"${QP_PYTHON}")', + 'export PATH=$(qp_prepend_export "PATH" "${QP_PYTHON}":"${QP_ROOT}"/bin:"${QP_ROOT}"/ocaml)', + 'export LD_LIBRARY_PATH=$(qp_prepend_export "LD_LIBRARY_PATH" "${QP_ROOT}"/lib:"${QP_ROOT}"/lib64)', + 'export LIBRARY_PATH=$(qp_prepend_export "LIBRARY_PATH" "${QP_ROOT}"/lib:"${QP_ROOT}"/lib64)', + 'export C_INCLUDE_PATH=$(qp_prepend_export "C_INCLUDE_PATH" "${QP_ROOT}"/include)', '', 'if [[ $SHELL == "bash" ]] ; then', ' source ${QP_ROOT}/install/EZFIO/Bash/ezfio.sh', diff --git a/install/scripts/install_ocaml.sh b/install/scripts/install_ocaml.sh index f322bd0b..9e8a2b25 100755 --- a/install/scripts/install_ocaml.sh +++ b/install/scripts/install_ocaml.sh @@ -5,11 +5,12 @@ QP_ROOT=$PWD cd - # Normal installation -PACKAGES="core.v0.9.1 cryptokit.1.10 ocamlfind sexplib.v0.9.1 ZMQ ppx_sexp_conv ppx_deriving" +PACKAGES="core.v0.10.0 cryptokit ocamlfind sexplib.v0.10.0 ZMQ ppx_sexp_conv ppx_deriving" # Needed for ZeroMQ export C_INCLUDE_PATH="${QP_ROOT}"/include:"${C_INCLUDE_PATH}" export LIBRARY_PATH="${QP_ROOT}"/lib:"${LIBRARY_PATH}" +export LDFLAGS="-L$QP_ROOT/lib" export LD_LIBRARY_PATH="${QP_ROOT}"/lib:"${LD_LIBRARY_PATH}" # return 0 if program version is equal or greater than check version @@ -64,7 +65,7 @@ fi cd Downloads || exit 1 chmod +x ocaml.sh || exit 1 -echo N | ./ocaml.sh ${QP_ROOT}/bin/ 4.04.2 || exit 1 +echo N | ./ocaml.sh ${QP_ROOT}/bin/ 4.06.0 || exit 1 ${QP_ROOT}/bin/opam config setup -a -q || exit 1 diff --git a/ocaml/TaskServer.ml b/ocaml/TaskServer.ml index 31d6ab3b..170e011a 100644 --- a/ocaml/TaskServer.ml +++ b/ocaml/TaskServer.ml @@ -1,6 +1,7 @@ open Core open Qptypes +module StringHashtbl = Hashtbl.Make(String) type pub_state = | Waiting @@ -28,7 +29,7 @@ type t = progress_bar : Progress_bar.t option ; running : bool; accepting_clients : bool; - data : (string, string) Hashtbl.t; + data : string StringHashtbl.t; } @@ -208,7 +209,7 @@ let end_job msg program_state rep_socket pair_socket = address_inproc = None; running = true; accepting_clients = false; - data = Hashtbl.create ~hashable:String.hashable (); + data = StringHashtbl.create (); } and wait n = @@ -592,7 +593,7 @@ let put_data msg rest_of_msg program_state rep_socket = in let success () = - Hashtbl.set program_state.data ~key ~data:value ; + StringHashtbl.set program_state.data ~key ~data:value ; Message.PutDataReply (Message.PutDataReply_msg.create ()) |> Message.to_string |> ZMQ.Socket.send rep_socket; @@ -622,7 +623,7 @@ let get_data msg program_state rep_socket = let success () = let value = - match Hashtbl.find program_state.data key with + match StringHashtbl.find program_state.data key with | Some value -> value | None -> "" in @@ -776,7 +777,7 @@ let run ~port = address_inproc = None; progress_bar = None ; accepting_clients = false; - data = Hashtbl.create ~hashable:String.hashable (); + data = StringHashtbl.create (); } in diff --git a/ocaml/qp_create_ezfio_from_xyz.ml b/ocaml/qp_create_ezfio_from_xyz.ml index 93c8c8ff..737052ee 100644 --- a/ocaml/qp_create_ezfio_from_xyz.ml +++ b/ocaml/qp_create_ezfio_from_xyz.ml @@ -665,7 +665,7 @@ let run ?o b au c d m p cart xyz_file = let command = - Command.basic + Command.basic_spec ~summary: "Quantum Package command" ~readme:(fun () -> " diff --git a/ocaml/qp_create_guess.ml b/ocaml/qp_create_guess.ml index b841c350..71a5b296 100644 --- a/ocaml/qp_create_guess.ml +++ b/ocaml/qp_create_guess.ml @@ -128,7 +128,7 @@ let spec = +> anon ("ezfio_file" %: string) let () = - Command.basic + Command.basic_spec ~summary: "Quantum Package command" ~readme:( fun () -> " Creates an open-shell multiplet initial guess\n\n" ) diff --git a/ocaml/qp_find_pi_space.ml b/ocaml/qp_find_pi_space.ml index 0f5f7365..dcd671ce 100644 --- a/ocaml/qp_find_pi_space.ml +++ b/ocaml/qp_find_pi_space.ml @@ -95,7 +95,7 @@ let spec = let command = - Command.basic + Command.basic_spec ~summary: "Quantum Package command" ~readme:(fun () -> "Find all the pi molecular orbitals to create a pi space. diff --git a/ocaml/qp_print.ml b/ocaml/qp_print.ml index ea52bd7f..efbdc01e 100644 --- a/ocaml/qp_print.ml +++ b/ocaml/qp_print.ml @@ -141,7 +141,7 @@ let run_o ~action ezfio_filename = ;; let command = - Command.basic + Command.basic_spec ~summary: "Quantum Package command" ~readme:(fun () -> " diff --git a/ocaml/qp_run.ml b/ocaml/qp_run.ml index f426932b..f3f0b14e 100644 --- a/ocaml/qp_run.ml +++ b/ocaml/qp_run.ml @@ -150,7 +150,7 @@ let spec = let () = - Command.basic + Command.basic_spec ~summary: "Quantum Package command" ~readme:( fun () -> " Executes a Quantum Package binary file among these:\n\n" diff --git a/ocaml/qp_set_mo_class.ml b/ocaml/qp_set_mo_class.ml index ef2cc977..e12a2d75 100644 --- a/ocaml/qp_set_mo_class.ml +++ b/ocaml/qp_set_mo_class.ml @@ -323,7 +323,7 @@ let spec = let command = - Command.basic + Command.basic_spec ~summary: "Quantum Package command" ~readme:(fun () -> "Set the orbital classes in an EZFIO directory diff --git a/scripts/ezfio_interface/upgrade_1.0_2.0.sh b/scripts/ezfio_interface/upgrade_1.0_2.0.sh deleted file mode 100755 index ec0ab770..00000000 --- a/scripts/ezfio_interface/upgrade_1.0_2.0.sh +++ /dev/null @@ -1,27 +0,0 @@ -#!/bin/bash -# Convert a old ezfio file (with option.irp.f ezfio_default) -# into a new EZFIO.cfg type - -# Hartree Fock -# Changin the case, don't know if is needed or not -mv $1/Hartree_Fock $1/hartree_fock 2> /dev/null - -mv $1/hartree_Fock/thresh_SCF $1/hartree_fock/thresh_scf 2> /dev/null - -# BiInts -mv $1/bi_integrals $1/bielect_integrals 2> /dev/null - -if [ -f $1/bielect_integrals/read_ao_integrals ]; then - if [ `cat $1/bielect_integrals/read_ao_integrals` -eq "True" ] - then - echo "Read" > $1/bielect_integrals/disk_access_ao_integrals - - elif [ `cat bielect_integrals/write_ao_integrals` -eq "True" ] - then - echo "Write" > $1/bielect_integrals/disk_access_ao_integrals - - else - echo "None" > $1/bielect_integrals/disk_access_ao_integrals - - fi -fi \ No newline at end of file diff --git a/scripts/qp_upgrade_ocaml.sh b/scripts/qp_upgrade_ocaml.sh new file mode 100755 index 00000000..ad2156e9 --- /dev/null +++ b/scripts/qp_upgrade_ocaml.sh @@ -0,0 +1,19 @@ +#!/bin/bash + +OCAML_VERSION="4.06.0" +PACKAGES="core.v0.10.0 cryptokit ocamlfind sexplib.v0.10.0 ZMQ ppx_sexp_conv ppx_deriving" + +if [[ -z ${QP_ROOT} ]] +then + print "The QP_ROOT environment variable is not set." + print "Please reload the quantum_package.rc file." + exit -1 +fi + +cd $QP_ROOT/ocaml +opam update +opam switch ${OCAML_VERSION} +eval `opam config env` +opam install -y ${PACKAGES} || echo "Upgrade failed. You can try running +configure ; $0" + From cee0b40463849c06aec515014407d3155ced6c9b Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 22 Jan 2018 01:46:46 +0100 Subject: [PATCH 10/14] Fixed qp_edit --- scripts/ezfio_interface/qp_edit_template | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/ezfio_interface/qp_edit_template b/scripts/ezfio_interface/qp_edit_template index 55a35f83..2e2c26c6 100644 --- a/scripts/ezfio_interface/qp_edit_template +++ b/scripts/ezfio_interface/qp_edit_template @@ -256,7 +256,7 @@ let spec = let command = - Command.basic + Command.basic_spec ~summary: "Quantum Package command" ~readme:(fun () -> " From 5f6349e7ac3c79d1d1a109f1cc8bfb41969119b0 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 8 Feb 2018 17:48:47 +0100 Subject: [PATCH 11/14] Generalized Davdison for dressed methods --- ocaml/Input_determinants_by_hand.ml | 12 +- plugins/CAS_SD_ZMQ/NEEDED_CHILDREN_MODULES | 2 +- plugins/CID/NEEDED_CHILDREN_MODULES | 2 +- plugins/CID_selected/NEEDED_CHILDREN_MODULES | 2 +- plugins/CIS/NEEDED_CHILDREN_MODULES | 2 +- plugins/CISD/NEEDED_CHILDREN_MODULES | 2 +- plugins/Casino/NEEDED_CHILDREN_MODULES | 2 +- plugins/DavidsonDressed/README.rst | 14 + .../diagonalization_hs2_dressed.irp.f | 103 +- .../DavidsonUndressed/NEEDED_CHILDREN_MODULES | 1 + plugins/DavidsonUndressed/README.rst | 14 + .../DavidsonUndressed}/davidson_slave.irp.f | 0 .../diag_restart_save_all_states.irp.f | 0 .../diag_restart_save_lowest_state.irp.f | 0 .../diag_restart_save_one_state.irp.f | 0 .../guess_lowest_state.irp.f | 0 .../print_H_matrix_restart.irp.f | 0 .../DavidsonUndressed}/print_energy.irp.f | 0 plugins/Full_CI_ZMQ/NEEDED_CHILDREN_MODULES | 2 +- plugins/MRCC_Utils/H_apply.irp.f | 10 +- plugins/MRCC_Utils/davidson.irp.f | 1128 ----------------- plugins/MRCC_Utils/mrcc_dress.irp.f | 5 +- plugins/MRCC_Utils/mrcc_utils.irp.f | 204 --- plugins/MRCC_Utils/multi_state.irp.f | 101 -- plugins/Properties/NEEDED_CHILDREN_MODULES | 2 +- plugins/Psiref_CAS/NEEDED_CHILDREN_MODULES | 2 +- .../UndressedMethod/NEEDED_CHILDREN_MODULES | 1 + plugins/UndressedMethod/README.rst | 14 + .../null_dressing_vector.irp.f | 10 + plugins/mrcepa0/NEEDED_CHILDREN_MODULES | 2 +- plugins/mrcepa0/dressing.irp.f | 424 ++----- plugins/mrcepa0/dressing_slave.irp.f | 52 +- plugins/mrcepa0/mrcepa0_general.irp.f | 8 +- plugins/mrcepa0/run_mrcc_slave.irp.f | 10 +- src/Davidson/NEEDED_CHILDREN_MODULES | 2 +- src/Davidson/diagonalize_CI.irp.f | 2 +- 36 files changed, 308 insertions(+), 1827 deletions(-) create mode 100644 plugins/DavidsonDressed/README.rst rename src/Davidson/diagonalization_hs2.irp.f => plugins/DavidsonDressed/diagonalization_hs2_dressed.irp.f (84%) create mode 100644 plugins/DavidsonUndressed/NEEDED_CHILDREN_MODULES create mode 100644 plugins/DavidsonUndressed/README.rst rename {src/Davidson => plugins/DavidsonUndressed}/davidson_slave.irp.f (100%) rename src/Davidson/diagonalize_restart_and_save_all_states.irp.f => plugins/DavidsonUndressed/diag_restart_save_all_states.irp.f (100%) rename src/Davidson/diagonalize_restart_and_save_lowest_state.irp.f => plugins/DavidsonUndressed/diag_restart_save_lowest_state.irp.f (100%) rename src/Davidson/diagonalize_restart_and_save_one_state.irp.f => plugins/DavidsonUndressed/diag_restart_save_one_state.irp.f (100%) rename {src/Davidson => plugins/DavidsonUndressed}/guess_lowest_state.irp.f (100%) rename {src/Davidson => plugins/DavidsonUndressed}/print_H_matrix_restart.irp.f (100%) rename {src/Davidson => plugins/DavidsonUndressed}/print_energy.irp.f (100%) delete mode 100644 plugins/MRCC_Utils/davidson.irp.f delete mode 100644 plugins/MRCC_Utils/multi_state.irp.f create mode 100644 plugins/UndressedMethod/NEEDED_CHILDREN_MODULES create mode 100644 plugins/UndressedMethod/README.rst create mode 100644 plugins/UndressedMethod/null_dressing_vector.irp.f diff --git a/ocaml/Input_determinants_by_hand.ml b/ocaml/Input_determinants_by_hand.ml index 90174e18..ecc68e70 100644 --- a/ocaml/Input_determinants_by_hand.ml +++ b/ocaml/Input_determinants_by_hand.ml @@ -93,8 +93,16 @@ end = struct ;; let write_n_states n = - States_number.to_int n - |> Ezfio.set_determinants_n_states + let n_states = + States_number.to_int n + in + Ezfio.set_determinants_n_states n_states; + let data = + Array.create n_states 1. + |> Array.to_list + in + Ezfio.ezfio_array_of_list ~rank:1 ~dim:[| n_states |] ~data + |> Ezfio.set_determinants_state_average_weight ;; let write_state_average_weight data = diff --git a/plugins/CAS_SD_ZMQ/NEEDED_CHILDREN_MODULES b/plugins/CAS_SD_ZMQ/NEEDED_CHILDREN_MODULES index ae599426..91dd3eff 100644 --- a/plugins/CAS_SD_ZMQ/NEEDED_CHILDREN_MODULES +++ b/plugins/CAS_SD_ZMQ/NEEDED_CHILDREN_MODULES @@ -1,2 +1,2 @@ -Generators_CAS Perturbation Selectors_CASSD ZMQ +Generators_CAS Perturbation Selectors_CASSD ZMQ DavidsonUndressed diff --git a/plugins/CID/NEEDED_CHILDREN_MODULES b/plugins/CID/NEEDED_CHILDREN_MODULES index 1632a44d..3272abe5 100644 --- a/plugins/CID/NEEDED_CHILDREN_MODULES +++ b/plugins/CID/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -Selectors_full SingleRefMethod Davidson +Selectors_full SingleRefMethod DavidsonUndressed diff --git a/plugins/CID_selected/NEEDED_CHILDREN_MODULES b/plugins/CID_selected/NEEDED_CHILDREN_MODULES index 1e0c52c2..ea9febd6 100644 --- a/plugins/CID_selected/NEEDED_CHILDREN_MODULES +++ b/plugins/CID_selected/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -Perturbation CID +Perturbation CID DavidsonUndressed diff --git a/plugins/CIS/NEEDED_CHILDREN_MODULES b/plugins/CIS/NEEDED_CHILDREN_MODULES index 1632a44d..3272abe5 100644 --- a/plugins/CIS/NEEDED_CHILDREN_MODULES +++ b/plugins/CIS/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -Selectors_full SingleRefMethod Davidson +Selectors_full SingleRefMethod DavidsonUndressed diff --git a/plugins/CISD/NEEDED_CHILDREN_MODULES b/plugins/CISD/NEEDED_CHILDREN_MODULES index 1632a44d..3272abe5 100644 --- a/plugins/CISD/NEEDED_CHILDREN_MODULES +++ b/plugins/CISD/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -Selectors_full SingleRefMethod Davidson +Selectors_full SingleRefMethod DavidsonUndressed diff --git a/plugins/Casino/NEEDED_CHILDREN_MODULES b/plugins/Casino/NEEDED_CHILDREN_MODULES index 34de8ddb..2a87d1c1 100644 --- a/plugins/Casino/NEEDED_CHILDREN_MODULES +++ b/plugins/Casino/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -Determinants Davidson +Determinants DavidsonUndressed diff --git a/plugins/DavidsonDressed/README.rst b/plugins/DavidsonDressed/README.rst new file mode 100644 index 00000000..ce9c78ba --- /dev/null +++ b/plugins/DavidsonDressed/README.rst @@ -0,0 +1,14 @@ +=============== +DavidsonDressed +=============== + +Davidson with single-column dressing + +Needed Modules +============== +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. +Documentation +============= +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. diff --git a/src/Davidson/diagonalization_hs2.irp.f b/plugins/DavidsonDressed/diagonalization_hs2_dressed.irp.f similarity index 84% rename from src/Davidson/diagonalization_hs2.irp.f rename to plugins/DavidsonDressed/diagonalization_hs2_dressed.irp.f index 2dfe468e..8a477b5a 100644 --- a/src/Davidson/diagonalization_hs2.irp.f +++ b/plugins/DavidsonDressed/diagonalization_hs2_dressed.irp.f @@ -1,4 +1,17 @@ -subroutine davidson_diag_hs2(dets_in,u_in,s2_out,dim_in,energies,sze,N_st,N_st_diag,Nint,iunit) +BEGIN_PROVIDER [ integer, dressed_column_idx, (N_states) ] + implicit none + BEGIN_DOC + ! Index of the dressed columns + END_DOC + integer :: i + double precision :: tmp + integer, external :: idamax + do i=1,N_states + dressed_column_idx(i) = idamax(size(psi_coef,1), psi_coef(1,i), 1) + enddo +END_PROVIDER + +subroutine davidson_diag_hs2(dets_in,u_in,s2_out,dim_in,energies,sze,N_st,N_st_diag,Nint,dressing_state) use bitmasks implicit none BEGIN_DOC @@ -15,41 +28,45 @@ subroutine davidson_diag_hs2(dets_in,u_in,s2_out,dim_in,energies,sze,N_st,N_st_d ! ! N_st : Number of eigenstates ! - ! iunit : Unit number for the I/O - ! ! Initial guess vectors are not necessarily orthonormal END_DOC - integer, intent(in) :: dim_in, sze, N_st, N_st_diag, Nint, iunit + integer, intent(in) :: dim_in, sze, N_st, N_st_diag, Nint integer(bit_kind), intent(in) :: dets_in(Nint,2,sze) double precision, intent(inout) :: u_in(dim_in,N_st_diag) double precision, intent(out) :: energies(N_st_diag), s2_out(N_st_diag) - double precision, allocatable :: H_jj(:) + integer, intent(in) :: dressing_state + double precision, allocatable :: H_jj(:), S2_jj(:) - double precision :: diag_H_mat_elem, diag_S_mat_elem + double precision, external :: diag_H_mat_elem, diag_S_mat_elem integer :: i ASSERT (N_st > 0) ASSERT (sze > 0) ASSERT (Nint > 0) ASSERT (Nint == N_int) PROVIDE mo_bielec_integrals_in_map - allocate(H_jj(sze) ) + allocate(H_jj(sze),S2_jj(sze)) + H_jj(1) = diag_h_mat_elem(dets_in(1,1,1),Nint) !$OMP PARALLEL DEFAULT(NONE) & !$OMP SHARED(sze,H_jj, dets_in,Nint) & !$OMP PRIVATE(i) !$OMP DO SCHEDULE(static) - do i=1,sze + do i=2,sze H_jj(i) = diag_H_mat_elem(dets_in(1,1,i),Nint) enddo !$OMP END DO !$OMP END PARALLEL - call davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_out,energies,dim_in,sze,N_st,N_st_diag,Nint,iunit) - deallocate (H_jj) + if (dressing_state > 0) then + H_jj(dressed_column_idx(dressing_state)) += dressing_column_h(dressed_column_idx(dressing_state),dressing_state) + endif + + call davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_out,energies,dim_in,sze,N_st,N_st_diag,Nint,dressing_state) + deallocate (H_jj,S2_jj) end -subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_st,N_st_diag,Nint,iunit) +subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_st,N_st_diag,Nint,dressing_state) use bitmasks implicit none BEGIN_DOC @@ -72,15 +89,13 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_ ! ! N_st_diag : Number of states in which H is diagonalized. Assumed > sze ! - ! iunit : Unit for the I/O - ! ! Initial guess vectors are not necessarily orthonormal END_DOC integer, intent(in) :: dim_in, sze, N_st, N_st_diag, Nint integer(bit_kind), intent(in) :: dets_in(Nint,2,sze) double precision, intent(in) :: H_jj(sze) + integer, intent(in) :: dressing_state double precision, intent(inout) :: s2_out(N_st_diag) - integer, intent(in) :: iunit double precision, intent(inout) :: u_in(dim_in,N_st_diag) double precision, intent(out) :: energies(N_st_diag) @@ -88,7 +103,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_ integer :: i,j,k,l,m logical :: converged - double precision :: u_dot_v, u_dot_u + double precision, external :: u_dot_v, u_dot_u integer :: k_pairs, kl @@ -101,7 +116,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_ character*(16384) :: write_buffer double precision :: to_print(3,N_st) double precision :: cpu, wall - integer :: shift, shift2, itermax + integer :: shift, shift2, itermax, istate double precision :: r1, r2 logical :: state_ok(N_st_diag*davidson_sze_max) include 'constants.include.F' @@ -117,35 +132,35 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_ PROVIDE nuclear_repulsion expected_s2 psi_bilinear_matrix_order psi_bilinear_matrix_order_reverse - call write_time(iunit) + call write_time(6) call wall_time(wall) call cpu_time(cpu) - write(iunit,'(A)') '' - write(iunit,'(A)') 'Davidson Diagonalization' - write(iunit,'(A)') '------------------------' - write(iunit,'(A)') '' - call write_int(iunit,N_st,'Number of states') - call write_int(iunit,N_st_diag,'Number of states in diagonalization') - call write_int(iunit,sze,'Number of determinants') + write(6,'(A)') '' + write(6,'(A)') 'Davidson Diagonalization' + write(6,'(A)') '------------------------' + write(6,'(A)') '' + call write_int(6,N_st,'Number of states') + call write_int(6,N_st_diag,'Number of states in diagonalization') + call write_int(6,sze,'Number of determinants') r1 = 8.d0*(3.d0*dble(sze*N_st_diag*itermax+5.d0*(N_st_diag*itermax)**2 & + 4.d0*(N_st_diag*itermax)+nproc*(4.d0*N_det_alpha_unique+2.d0*N_st_diag*sze)))/(1024.d0**3) - call write_double(iunit, r1, 'Memory(Gb)') - write(iunit,'(A)') '' + call write_double(6, r1, 'Memory(Gb)') + write(6,'(A)') '' write_buffer = '=====' do i=1,N_st write_buffer = trim(write_buffer)//' ================ =========== ===========' enddo - write(iunit,'(A)') write_buffer(1:6+41*N_states) + write(6,'(A)') write_buffer(1:6+41*N_states) write_buffer = 'Iter' do i=1,N_st write_buffer = trim(write_buffer)//' Energy S^2 Residual ' enddo - write(iunit,'(A)') write_buffer(1:6+41*N_states) + write(6,'(A)') write_buffer(1:6+41*N_states) write_buffer = '=====' do i=1,N_st write_buffer = trim(write_buffer)//' ================ =========== ===========' enddo - write(iunit,'(A)') write_buffer(1:6+41*N_states) + write(6,'(A)') write_buffer(1:6+41*N_states) allocate( & @@ -225,7 +240,21 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_ call H_S2_u_0_nstates_openmp(W(1,shift+1),S(1,shift+1),U(1,shift+1),N_st_diag,sze) endif - + if (dressing_state > 0) then + + do istate=1,N_st_diag + l = dressed_column_idx(dressing_state) + do i=1,sze + W(i,shift+istate) += dressing_column_h(i,dressing_state) * U(l,shift+istate) + S(i,shift+istate) += dressing_column_s(i,dressing_state) * U(l,shift+istate) + W(l,shift+istate) += dressing_column_h(i,dressing_state) * U(i,shift+istate) + S(l,shift+istate) += dressing_column_s(i,dressing_state) * U(i,shift+istate) + enddo + W(l,shift+istate) -= dressing_column_h(l,dressing_state) * U(l,shift+istate) + S(l,shift+istate) -= dressing_column_s(l,dressing_state) * U(l,shift+istate) + enddo + endif + ! Compute h_kl = = ! ------------------------------------------- @@ -399,7 +428,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_ endif enddo - write(iunit,'(1X,I3,1X,100(1X,F16.10,1X,F11.6,1X,E11.3))') iter, to_print(1:3,1:N_st) + write(6,'(1X,I3,1X,100(1X,F16.10,1X,F11.6,1X,E11.3))') iter, to_print(1:3,1:N_st) call davidson_converged(lambda,residual_norm,wall,iter,cpu,N_st,converged) do k=1,N_st if (residual_norm(k) > 1.e8) then @@ -429,9 +458,9 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_ do i=1,N_st write_buffer = trim(write_buffer)//' ================ =========== ===========' enddo - write(iunit,'(A)') trim(write_buffer) - write(iunit,'(A)') '' - call write_time(iunit) + write(6,'(A)') trim(write_buffer) + write(6,'(A)') '' + call write_time(6) deallocate ( & W, residual_norm, & @@ -443,6 +472,12 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_ ) end + + + + + + subroutine u_0_H_u_0(e_0,u_0,n,keys_tmp,Nint,N_st,sze) use bitmasks implicit none diff --git a/plugins/DavidsonUndressed/NEEDED_CHILDREN_MODULES b/plugins/DavidsonUndressed/NEEDED_CHILDREN_MODULES new file mode 100644 index 00000000..25180044 --- /dev/null +++ b/plugins/DavidsonUndressed/NEEDED_CHILDREN_MODULES @@ -0,0 +1 @@ +Davidson UndressedMethod diff --git a/plugins/DavidsonUndressed/README.rst b/plugins/DavidsonUndressed/README.rst new file mode 100644 index 00000000..e11d0703 --- /dev/null +++ b/plugins/DavidsonUndressed/README.rst @@ -0,0 +1,14 @@ +================= +DavidsonUndressed +================= + +Module for main files with undressed Davidson + +Needed Modules +============== +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. +Documentation +============= +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. diff --git a/src/Davidson/davidson_slave.irp.f b/plugins/DavidsonUndressed/davidson_slave.irp.f similarity index 100% rename from src/Davidson/davidson_slave.irp.f rename to plugins/DavidsonUndressed/davidson_slave.irp.f diff --git a/src/Davidson/diagonalize_restart_and_save_all_states.irp.f b/plugins/DavidsonUndressed/diag_restart_save_all_states.irp.f similarity index 100% rename from src/Davidson/diagonalize_restart_and_save_all_states.irp.f rename to plugins/DavidsonUndressed/diag_restart_save_all_states.irp.f diff --git a/src/Davidson/diagonalize_restart_and_save_lowest_state.irp.f b/plugins/DavidsonUndressed/diag_restart_save_lowest_state.irp.f similarity index 100% rename from src/Davidson/diagonalize_restart_and_save_lowest_state.irp.f rename to plugins/DavidsonUndressed/diag_restart_save_lowest_state.irp.f diff --git a/src/Davidson/diagonalize_restart_and_save_one_state.irp.f b/plugins/DavidsonUndressed/diag_restart_save_one_state.irp.f similarity index 100% rename from src/Davidson/diagonalize_restart_and_save_one_state.irp.f rename to plugins/DavidsonUndressed/diag_restart_save_one_state.irp.f diff --git a/src/Davidson/guess_lowest_state.irp.f b/plugins/DavidsonUndressed/guess_lowest_state.irp.f similarity index 100% rename from src/Davidson/guess_lowest_state.irp.f rename to plugins/DavidsonUndressed/guess_lowest_state.irp.f diff --git a/src/Davidson/print_H_matrix_restart.irp.f b/plugins/DavidsonUndressed/print_H_matrix_restart.irp.f similarity index 100% rename from src/Davidson/print_H_matrix_restart.irp.f rename to plugins/DavidsonUndressed/print_H_matrix_restart.irp.f diff --git a/src/Davidson/print_energy.irp.f b/plugins/DavidsonUndressed/print_energy.irp.f similarity index 100% rename from src/Davidson/print_energy.irp.f rename to plugins/DavidsonUndressed/print_energy.irp.f diff --git a/plugins/Full_CI_ZMQ/NEEDED_CHILDREN_MODULES b/plugins/Full_CI_ZMQ/NEEDED_CHILDREN_MODULES index 6736cc4e..cc81a88f 100644 --- a/plugins/Full_CI_ZMQ/NEEDED_CHILDREN_MODULES +++ b/plugins/Full_CI_ZMQ/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -Perturbation Selectors_full Generators_full ZMQ FourIdx MPI +Perturbation Selectors_full Generators_full ZMQ FourIdx MPI DavidsonUndressed diff --git a/plugins/MRCC_Utils/H_apply.irp.f b/plugins/MRCC_Utils/H_apply.irp.f index d8dfb62d..7fedd1a8 100644 --- a/plugins/MRCC_Utils/H_apply.irp.f +++ b/plugins/MRCC_Utils/H_apply.irp.f @@ -3,19 +3,17 @@ BEGIN_SHELL [ /usr/bin/env python ] from generate_h_apply import * s = H_apply("mrcc") -s.data["parameters"] = ", delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref" +s.data["parameters"] = ", delta_ij_, Nstates, Ndet_non_ref, Ndet_ref" s.data["declarations"] += """ integer, intent(in) :: Nstates, Ndet_ref, Ndet_non_ref double precision, intent(in) :: delta_ij_(Nstates, Ndet_non_ref, Ndet_ref) - double precision, intent(in) :: delta_ii_(Nstates, Ndet_ref) """ -s.data["keys_work"] = "call mrcc_dress(delta_ij_,delta_ii_,Nstates,Ndet_non_ref,Ndet_ref,i_generator,key_idx,keys_out,N_int,iproc,key_mask)" -s.data["params_post"] += ", delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref" -s.data["params_main"] += "delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref" +s.data["keys_work"] = "call mrcc_dress(delta_ij_,Nstates,Ndet_non_ref,Ndet_ref,i_generator,key_idx,keys_out,N_int,iproc,key_mask)" +s.data["params_post"] += ", delta_ij_, Nstates, Ndet_non_ref, Ndet_ref" +s.data["params_main"] += "delta_ij_, Nstates, Ndet_non_ref, Ndet_ref" s.data["decls_main"] += """ integer, intent(in) :: Ndet_ref, Ndet_non_ref, Nstates double precision, intent(in) :: delta_ij_(Nstates,Ndet_non_ref,Ndet_ref) - double precision, intent(in) :: delta_ii_(Nstates,Ndet_ref) """ s.data["finalization"] = "" s.data["copy_buffer"] = "" diff --git a/plugins/MRCC_Utils/davidson.irp.f b/plugins/MRCC_Utils/davidson.irp.f deleted file mode 100644 index 14b8d816..00000000 --- a/plugins/MRCC_Utils/davidson.irp.f +++ /dev/null @@ -1,1128 +0,0 @@ -subroutine davidson_diag_mrcc(dets_in,u_in,energies,dim_in,sze,N_st,N_st_diag,Nint,iunit,istate) - use bitmasks - implicit none - BEGIN_DOC - ! Davidson diagonalization. - ! - ! dets_in : bitmasks corresponding to determinants - ! - ! u_in : guess coefficients on the various states. Overwritten - ! on exit - ! - ! dim_in : leftmost dimension of u_in - ! - ! sze : Number of determinants - ! - ! N_st : Number of eigenstates - ! - ! iunit : Unit number for the I/O - ! - ! Initial guess vectors are not necessarily orthonormal - END_DOC - integer, intent(in) :: dim_in, sze, N_st, Nint, iunit, istate, N_st_diag - integer(bit_kind), intent(in) :: dets_in(Nint,2,sze) - double precision, intent(inout) :: u_in(dim_in,N_st_diag) - double precision, intent(out) :: energies(N_st_diag) - double precision, allocatable :: H_jj(:) - - double precision :: diag_h_mat_elem - integer :: i - ASSERT (N_st > 0) - ASSERT (N_st_diag >= N_st) - ASSERT (sze > 0) - ASSERT (Nint > 0) - ASSERT (Nint == N_int) - PROVIDE mo_bielec_integrals_in_map - allocate(H_jj(sze)) - - H_jj(1) = diag_h_mat_elem(dets_in(1,1,1),Nint) - !$OMP PARALLEL DEFAULT(NONE) & - !$OMP SHARED(sze,H_jj,N_det_ref,dets_in,Nint,istate,delta_ii,idx_ref) & - !$OMP PRIVATE(i) - !$OMP DO - do i=2,sze - H_jj(i) = diag_h_mat_elem(dets_in(1,1,i),Nint) - enddo - !$OMP END DO - !$OMP END PARALLEL - - do i=1,N_det_ref - H_jj(idx_ref(i)) += delta_ii(istate,i) - enddo - call davidson_diag_hjj_mrcc(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,N_st_diag,Nint,iunit,istate) - deallocate (H_jj) -end - -subroutine davidson_diag_hjj_mrcc(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,N_st_diag,Nint,iunit,istate) - use bitmasks - implicit none - BEGIN_DOC - ! Davidson diagonalization with specific diagonal elements of the H matrix - ! - ! H_jj : specific diagonal H matrix elements to diagonalize de Davidson - ! - ! dets_in : bitmasks corresponding to determinants - ! - ! u_in : guess coefficients on the various states. Overwritten - ! on exit - ! - ! dim_in : leftmost dimension of u_in - ! - ! sze : Number of determinants - ! - ! N_st : Number of eigenstates - ! - ! N_st_diag : Number of states in which H is diagonalized - ! - ! iunit : Unit for the I/O - ! - ! Initial guess vectors are not necessarily orthonormal - END_DOC - integer, intent(in) :: dim_in, sze, N_st, N_st_diag, Nint, istate - integer(bit_kind), intent(in) :: dets_in(Nint,2,sze) - double precision, intent(in) :: H_jj(sze) - integer, intent(in) :: iunit - double precision, intent(inout) :: u_in(dim_in,N_st_diag) - double precision, intent(out) :: energies(N_st_diag) - - integer :: iter - integer :: i,j,k,l,m - logical :: converged - - double precision, allocatable :: overlap(:,:) - double precision :: u_dot_v, u_dot_u - - integer :: k_pairs, kl - - integer :: iter2 - double precision, allocatable :: W(:,:,:), U(:,:,:), R(:,:) - double precision, allocatable :: y(:,:,:,:), h(:,:,:,:), lambda(:) - double precision, allocatable :: c(:), H_small(:,:) - double precision :: diag_h_mat_elem - double precision, allocatable :: residual_norm(:) - character*(16384) :: write_buffer - double precision :: to_print(2,N_st) - double precision :: cpu, wall - include 'constants.include.F' - - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: U, W, R, y, h, lambda - - PROVIDE nuclear_repulsion - - call write_time(iunit) - call wall_time(wall) - call cpu_time(cpu) - write(iunit,'(A)') '' - write(iunit,'(A)') 'Davidson Diagonalization' - write(iunit,'(A)') '------------------------' - write(iunit,'(A)') '' - call write_int(iunit,N_st,'Number of states') - call write_int(iunit,N_st_diag,'Number of states in diagonalization') - call write_int(iunit,sze,'Number of determinants') - call write_int(iunit,istate,'Using dressing for state ') - write(iunit,'(A)') '' - write_buffer = '===== ' - do i=1,N_st - write_buffer = trim(write_buffer)//' ================ ================' - enddo - write(iunit,'(A)') trim(write_buffer) - write_buffer = ' Iter' - do i=1,N_st - write_buffer = trim(write_buffer)//' Energy Residual' - enddo - write(iunit,'(A)') trim(write_buffer) - write_buffer = '===== ' - do i=1,N_st - write_buffer = trim(write_buffer)//' ================ ================' - enddo - write(iunit,'(A)') trim(write_buffer) - - allocate( & - W(sze,N_st_diag,davidson_sze_max), & - U(sze,N_st_diag,davidson_sze_max), & - R(sze,N_st_diag), & - h(N_st_diag,davidson_sze_max,N_st_diag,davidson_sze_max), & - y(N_st_diag,davidson_sze_max,N_st_diag,davidson_sze_max), & - residual_norm(N_st_diag), & - overlap(N_st_diag,N_st_diag), & - c(N_st_diag*davidson_sze_max), & - H_small(N_st_diag,N_st_diag), & - lambda(N_st_diag*davidson_sze_max)) - - ASSERT (N_st > 0) - ASSERT (N_st_diag >= N_st) - ASSERT (sze > 0) - ASSERT (Nint > 0) - ASSERT (Nint == N_int) - - ! Initialization - ! ============== - - - do k=1,N_st_diag - - if (k > N_st) then - do i=1,sze - double precision :: r1, r2 - call random_number(r1) - call random_number(r2) - u_in(i,k) = dsqrt(-2.d0*dlog(r1))*dcos(dtwo_pi*r2) - enddo - endif - - ! Gram-Schmidt - ! ------------ - call dgemv('T',sze,k-1,1.d0,u_in,size(u_in,1), & - u_in(1,k),1,0.d0,c,1) - call dgemv('N',sze,k-1,-1.d0,u_in,size(u_in,1), & - c,1,1.d0,u_in(1,k),1) - call normalize(u_in(1,k),sze) - enddo - - - - converged = .False. - do while (.not.converged) - - do k=1,N_st_diag - do i=1,sze - U(i,k,1) = u_in(i,k) - enddo - enddo - - do iter=1,davidson_sze_max-1 - - ! Compute |W_k> = \sum_i |i> - ! ----------------------------------------- - - call H_u_0_mrcc_nstates(W(1,1,iter),U(1,1,iter),H_jj,sze,dets_in,Nint,istate,N_st_diag,sze) - - - ! Compute h_kl = = - ! ------------------------------------------- - - - call dgemm('T','N', N_st_diag*iter, N_st_diag, sze, & - 1.d0, U, size(U,1), W(1,1,iter), size(W,1), & - 0.d0, h(1,1,1,iter), size(h,1)*size(h,2)) - - ! Diagonalize h - ! ------------- - call lapack_diag(lambda,y,h,N_st_diag*davidson_sze_max,N_st_diag*iter) - - ! Express eigenvectors of h in the determinant basis - ! -------------------------------------------------- - - do k=1,N_st_diag - do i=1,sze - U(i,k,iter+1) = 0.d0 - W(i,k,iter+1) = 0.d0 - enddo - enddo -! - call dgemm('N','N', sze, N_st_diag, N_st_diag*iter, & - 1.d0, U, size(U,1), y, size(y,1)*size(y,2), 0.d0, U(1,1,iter+1), size(U,1)) - call dgemm('N','N',sze,N_st_diag,N_st_diag*iter, & - 1.d0, W, size(W,1), y, size(y,1)*size(y,2), 0.d0, W(1,1,iter+1), size(W,1)) - - - ! Compute residual vector - ! ----------------------- - - do k=1,N_st_diag - do i=1,sze - R(i,k) = lambda(k) * U(i,k,iter+1) - W(i,k,iter+1) - enddo - if (k <= N_st) then - residual_norm(k) = u_dot_u(R(1,k),sze) - to_print(1,k) = lambda(k) + nuclear_repulsion - to_print(2,k) = residual_norm(k) - endif - enddo - - write(iunit,'(X,I3,X,100(X,F16.10,X,E16.6))') iter, to_print(:,1:N_st) - call davidson_converged(lambda,residual_norm,wall,iter,cpu,N_st,converged) - if (converged) then - exit - endif - - ! Davidson step - ! ------------- - - do k=1,N_st_diag - do i=1,sze - U(i,k,iter+1) = -1.d0/max(H_jj(i) - lambda(k),1.d-2) * R(i,k) - enddo - enddo - - ! Gram-Schmidt - ! ------------ - - do k=1,N_st_diag - - call dgemv('T',sze,N_st_diag*iter,1.d0,U,size(U,1), & - U(1,k,iter+1),1,0.d0,c,1) - call dgemv('N',sze,N_st_diag*iter,-1.d0,U,size(U,1), & - c,1,1.d0,U(1,k,iter+1),1) - - call dgemv('T',sze,k-1,1.d0,U(1,1,iter+1),size(U,1), & - U(1,k,iter+1),1,0.d0,c,1) - call dgemv('N',sze,k-1,-1.d0,U(1,1,iter+1),size(U,1), & - c,1,1.d0,U(1,k,iter+1),1) - - call normalize( U(1,k,iter+1), sze ) - enddo - - enddo - - if (.not.converged) then - iter = davidson_sze_max-1 - endif - - ! Re-contract to u_in - ! ----------- - - do k=1,N_st_diag - do i=1,sze - u_in(i,k) = 0.d0 - enddo - enddo - - call dgemm('N','N', sze, N_st_diag, N_st_diag*iter, 1.d0, & - U, size(U,1), y, N_st_diag*davidson_sze_max, & - 0.d0, u_in, size(u_in,1)) - - enddo - - do k=1,N_st_diag - energies(k) = lambda(k) - enddo - write_buffer = '===== ' - do i=1,N_st - write_buffer = trim(write_buffer)//' ================ ================' - enddo - write(iunit,'(A)') trim(write_buffer) - write(iunit,'(A)') '' - call write_time(iunit) - - deallocate ( & - W, residual_norm, & - U, overlap, & - R, c, & - h, & - y, & - lambda & - ) -end - - -subroutine u_0_H_u_0_mrcc_nstates(e_0,u_0,n,keys_tmp,Nint,istate,N_st,sze) - use bitmasks - implicit none - BEGIN_DOC - ! Computes e_0 = / - ! - ! n : number of determinants - ! - END_DOC - integer, intent(in) :: n,Nint,N_st,sze - double precision, intent(out) :: e_0(N_st) - double precision, intent(in) :: u_0(sze,N_st) - integer(bit_kind),intent(in) :: keys_tmp(Nint,2,n) - integer,intent(in) :: istate - - double precision, allocatable :: v_0(:,:), H_jj(:) - double precision :: u_dot_u,u_dot_v,diag_H_mat_elem - integer :: i,j - allocate(H_jj(n), v_0(sze,N_st)) - do i = 1, n - H_jj(i) = diag_H_mat_elem(keys_tmp(1,1,i),Nint) - enddo - - do i=1,N_det_ref - H_jj(idx_ref(i)) += delta_ii(istate,i) - enddo - - call H_u_0_mrcc_nstates(v_0,u_0,H_jj,n,keys_tmp,Nint,istate,N_st,sze) - do i=1,N_st - e_0(i) = u_dot_v(v_0(1,i),u_0(1,i),n)/u_dot_u(u_0(1,i),n) - enddo - deallocate(H_jj, v_0) -end - - -subroutine H_u_0_mrcc_nstates(v_0,u_0,H_jj,n,keys_tmp,Nint,istate_in,N_st,sze) - use bitmasks - implicit none - BEGIN_DOC - ! Computes v_0 = H|u_0> - ! - ! n : number of determinants - ! - ! H_jj : array of - END_DOC - integer, intent(in) :: n,Nint,istate_in,N_st,sze - double precision, intent(out) :: v_0(sze,N_st) - double precision, intent(in) :: u_0(sze,N_st) - double precision, intent(in) :: H_jj(n) - integer(bit_kind),intent(in) :: keys_tmp(Nint,2,n) - double precision :: hij - double precision, allocatable :: vt(:,:) - integer :: i,j,k,l, jj,ii - integer :: i0, j0 - integer(bit_kind) :: sorted_i(Nint) - - - integer,allocatable :: shortcut(:,:), sort_idx(:,:) - integer(bit_kind), allocatable :: sorted(:,:,:), version(:,:,:) - - - integer :: sh, sh2, ni, exa, ext, org_i, org_j, endi, pass, istate - - - ASSERT (Nint > 0) - ASSERT (Nint == N_int) - ASSERT (n>0) - PROVIDE ref_bitmask_energy - allocate (shortcut(0:n+1,2), sort_idx(n,2), sorted(Nint,n,2), version(Nint,n,2)) - v_0 = 0.d0 - - call sort_dets_ab_v(keys_tmp, sorted(1,1,1), sort_idx(1,1), shortcut(0,1), version(1,1,1), n, Nint) - call sort_dets_ba_v(keys_tmp, sorted(1,1,2), sort_idx(1,2), shortcut(0,2), version(1,1,2), n, Nint) - - !$OMP PARALLEL DEFAULT(NONE) & - !$OMP PRIVATE(i,hij,j,k,jj,vt,ii,sh,sh2,ni,exa,ext,org_i,org_j,endi,sorted_i,istate)& - !$OMP SHARED(n,H_jj,u_0,keys_tmp,Nint,v_0,sorted,shortcut,sort_idx,version,N_st,sze,& - !$OMP istate_in,delta_ij,N_det_ref,N_det_non_ref,idx_ref,idx_non_ref) - allocate(vt(sze,N_st)) - Vt = 0.d0 - - !$OMP DO SCHEDULE(static,1) - do sh=1,shortcut(0,1) - do sh2=sh,shortcut(0,1) - exa = 0 - do ni=1,Nint - exa = exa + popcnt(xor(version(ni,sh,1), version(ni,sh2,1))) - end do - if(exa > 2) then - cycle - end if - - do i=shortcut(sh,1),shortcut(sh+1,1)-1 - org_i = sort_idx(i,1) - if(sh==sh2) then - endi = i-1 - else - endi = shortcut(sh2+1,1)-1 - end if - do ni=1,Nint - sorted_i(ni) = sorted(ni,i,1) - enddo - - do j=shortcut(sh2,1),endi - org_j = sort_idx(j,1) - ext = exa - do ni=1,Nint - ext = ext + popcnt(xor(sorted_i(ni), sorted(ni,j,1))) - end do - if(ext <= 4) then - call i_H_j(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),Nint,hij) - do istate=1,N_st - vt (org_i,istate) = vt (org_i,istate) + hij*u_0(org_j,istate) - vt (org_j,istate) = vt (org_j,istate) + hij*u_0(org_i,istate) - enddo - endif - enddo - enddo - enddo - enddo - !$OMP END DO - - !$OMP DO SCHEDULE(static,1) - do sh=1,shortcut(0,2) - do i=shortcut(sh,2),shortcut(sh+1,2)-1 - org_i = sort_idx(i,2) - do j=shortcut(sh,2),i-1 - org_j = sort_idx(j,2) - ext = 0 - do ni=1,Nint - ext = ext + popcnt(xor(sorted(ni,i,2), sorted(ni,j,2))) - end do - if(ext == 4) then - call i_H_j(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),Nint,hij) - do istate=1,N_st - vt (org_i,istate) = vt (org_i,istate) + hij*u_0(org_j,istate) - vt (org_j,istate) = vt (org_j,istate) + hij*u_0(org_i,istate) - enddo - end if - end do - end do - enddo - !$OMP END DO - - !$OMP DO - do ii=1,n_det_ref - i = idx_ref(ii) - do jj = 1, n_det_non_ref - j = idx_non_ref(jj) - do istate=1,N_st - vt (i,istate) = vt (i,istate) + delta_ij(istate_in,jj,ii)*u_0(j,istate) - vt (j,istate) = vt (j,istate) + delta_ij(istate_in,jj,ii)*u_0(i,istate) - enddo - enddo - enddo - !$OMP END DO - - do istate=1,N_st - do i=n,1,-1 - !$OMP ATOMIC - v_0(i,istate) = v_0(i,istate) + vt(i,istate) - enddo - enddo - - deallocate(vt) - !$OMP END PARALLEL - - do istate=1,N_st - do i=1,n - v_0(i,istate) += H_jj(i) * u_0(i,istate) - enddo - enddo - deallocate (shortcut, sort_idx, sorted, version) - -end - - -subroutine davidson_diag_mrcc_hs2(dets_in,u_in,dim_in,energies,sze,N_st,N_st_diag,Nint,iunit,istate) - use bitmasks - implicit none - BEGIN_DOC - ! Davidson diagonalization. - ! - ! dets_in : bitmasks corresponding to determinants - ! - ! u_in : guess coefficients on the various states. Overwritten - ! on exit - ! - ! dim_in : leftmost dimension of u_in - ! - ! sze : Number of determinants - ! - ! N_st : Number of eigenstates - ! - ! iunit : Unit number for the I/O - ! - ! Initial guess vectors are not necessarily orthonormal - END_DOC - integer, intent(in) :: dim_in, sze, N_st, N_st_diag, Nint, iunit, istate - integer(bit_kind), intent(in) :: dets_in(Nint,2,sze) - double precision, intent(inout) :: u_in(dim_in,N_st_diag) - double precision, intent(out) :: energies(N_st_diag) - double precision, allocatable :: H_jj(:), S2_jj(:) - - double precision :: diag_h_mat_elem - integer :: i - ASSERT (N_st > 0) - ASSERT (sze > 0) - ASSERT (Nint > 0) - ASSERT (Nint == N_int) - PROVIDE mo_bielec_integrals_in_map - allocate(H_jj(sze), S2_jj(sze)) - - H_jj(1) = diag_h_mat_elem(dets_in(1,1,1),Nint) - call get_s2(dets_in(1,1,1),dets_in(1,1,1),Nint,S2_jj(1)) - !$OMP PARALLEL DEFAULT(NONE) & - !$OMP SHARED(sze,H_jj,S2_jj, dets_in,Nint,N_det_ref,delta_ii, & - !$OMP idx_ref, istate) & - !$OMP PRIVATE(i) - !$OMP DO - do i=2,sze - H_jj(i) = diag_h_mat_elem(dets_in(1,1,i),Nint) - call get_s2(dets_in(1,1,i),dets_in(1,1,i),Nint,S2_jj(i)) - enddo - !$OMP END DO - !$OMP END PARALLEL - - do i=1,N_det_ref - H_jj(idx_ref(i)) += delta_ii(istate,i) - enddo - - call davidson_diag_hjj_sjj_mrcc(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_st,N_st_diag,Nint,iunit,istate) - deallocate (H_jj,S2_jj) -end - - -subroutine davidson_diag_hjj_sjj_mrcc(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_st,N_st_diag,Nint,iunit,istate ) - use bitmasks - implicit none - BEGIN_DOC - ! Davidson diagonalization with specific diagonal elements of the H matrix - ! - ! H_jj : specific diagonal H matrix elements to diagonalize de Davidson - ! - ! S2_jj : specific diagonal S^2 matrix elements - ! - ! dets_in : bitmasks corresponding to determinants - ! - ! u_in : guess coefficients on the various states. Overwritten - ! on exit - ! - ! dim_in : leftmost dimension of u_in - ! - ! sze : Number of determinants - ! - ! N_st : Number of eigenstates - ! - ! N_st_diag : Number of states in which H is diagonalized. Assumed > sze - ! - ! iunit : Unit for the I/O - ! - ! Initial guess vectors are not necessarily orthonormal - END_DOC - integer, intent(in) :: dim_in, sze, N_st, N_st_diag, Nint, istate - integer(bit_kind), intent(in) :: dets_in(Nint,2,sze) - double precision, intent(in) :: H_jj(sze), S2_jj(sze) - integer, intent(in) :: iunit - double precision, intent(inout) :: u_in(dim_in,N_st_diag) - double precision, intent(out) :: energies(N_st_diag) - - integer :: iter - integer :: i,j,k,l,m - logical :: converged - - double precision :: u_dot_v, u_dot_u - - integer :: k_pairs, kl - - integer :: iter2 - double precision, allocatable :: W(:,:), U(:,:), S(:,:), overlap(:,:) - double precision, allocatable :: y(:,:), h(:,:), lambda(:), s2(:) - double precision, allocatable :: c(:), s_(:,:), s_tmp(:,:) - double precision :: diag_h_mat_elem - double precision, allocatable :: residual_norm(:) - character*(16384) :: write_buffer - double precision :: to_print(3,N_st) - double precision :: cpu, wall - integer :: shift, shift2, itermax - include 'constants.include.F' - - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: U, W, S, y, h, lambda - if (N_st_diag*3 > sze) then - print *, 'error in Davidson :' - print *, 'Increase n_det_max_jacobi to ', N_st_diag*3 - stop -1 - endif - - PROVIDE nuclear_repulsion - - call write_time(iunit) - call wall_time(wall) - call cpu_time(cpu) - write(iunit,'(A)') '' - write(iunit,'(A)') 'Davidson Diagonalization' - write(iunit,'(A)') '------------------------' - write(iunit,'(A)') '' - call write_int(iunit,N_st,'Number of states') - call write_int(iunit,N_st_diag,'Number of states in diagonalization') - call write_int(iunit,sze,'Number of determinants') - call write_int(iunit,istate,'Using dressing for state ') - - write(iunit,'(A)') '' - write_buffer = '===== ' - do i=1,N_st - write_buffer = trim(write_buffer)//' ================ =========== ===========' - enddo - write(iunit,'(A)') trim(write_buffer) - write_buffer = ' Iter' - do i=1,N_st - write_buffer = trim(write_buffer)//' Energy S^2 Residual ' - enddo - write(iunit,'(A)') trim(write_buffer) - write_buffer = '===== ' - do i=1,N_st - write_buffer = trim(write_buffer)//' ================ =========== ===========' - enddo - write(iunit,'(A)') trim(write_buffer) - - itermax = min(davidson_sze_max, sze/N_st_diag) - allocate( & - W(sze,N_st_diag*itermax), & - U(sze,N_st_diag*itermax), & - S(sze,N_st_diag*itermax), & - h(N_st_diag*itermax,N_st_diag*itermax), & - y(N_st_diag*itermax,N_st_diag*itermax), & - s_(N_st_diag*itermax,N_st_diag*itermax), & - s_tmp(N_st_diag*itermax,N_st_diag*itermax), & - residual_norm(N_st_diag), & - c(N_st_diag*itermax), & - s2(N_st_diag*itermax), & - overlap(N_st_diag*itermax,N_st_diag*itermax), & - lambda(N_st_diag*itermax)) - - h = 0.d0 - s_ = 0.d0 - s_tmp = 0.d0 - U = 0.d0 - W = 0.d0 - S = 0.d0 - y = 0.d0 - - - ASSERT (N_st > 0) - ASSERT (N_st_diag >= N_st) - ASSERT (sze > 0) - ASSERT (Nint > 0) - ASSERT (Nint == N_int) - - ! Davidson iterations - ! =================== - - converged = .False. - - double precision :: r1, r2 - do k=N_st+1,N_st_diag - u_in(k,k) = 10.d0 - do i=1,sze - call random_number(r1) - call random_number(r2) - r1 = dsqrt(-2.d0*dlog(r1)) - r2 = dtwo_pi*r2 - u_in(i,k) = r1*dcos(r2) - enddo - enddo - do k=1,N_st_diag - call normalize(u_in(1,k),sze) - enddo - - - do while (.not.converged) - - do k=1,N_st_diag - do i=1,sze - U(i,k) = u_in(i,k) - enddo - enddo - - do iter=1,davidson_sze_max-1 - - shift = N_st_diag*(iter-1) - shift2 = N_st_diag*iter - - call ortho_qr(U,size(U,1),sze,shift2) - - ! Compute |W_k> = \sum_i |i> - ! ----------------------------------------- - - call H_S2_u_0_mrcc_nstates(W(1,shift+1),S(1,shift+1),U(1,shift+1),H_jj,S2_jj,sze,dets_in,Nint,& - istate,N_st_diag,sze) - - - ! Compute h_kl = = - ! ------------------------------------------- - - - call dgemm('T','N', shift2, shift2, sze, & - 1.d0, U, size(U,1), W, size(W,1), & - 0.d0, h, size(h,1)) - - call dgemm('T','N', shift2, shift2, sze, & - 1.d0, U, size(U,1), S, size(S,1), & - 0.d0, s_, size(s_,1)) - -! ! Diagonalize S^2 -! ! --------------- -! -! call lapack_diag(s2,y,s_,size(s_,1),shift2) -! -! ! Rotate H in the basis of eigenfunctions of s2 -! ! --------------------------------------------- -! -! call dgemm('N','N',shift2,shift2,shift2, & -! 1.d0, h, size(h,1), y, size(y,1), & -! 0.d0, s_tmp, size(s_tmp,1)) -! -! call dgemm('T','N',shift2,shift2,shift2, & -! 1.d0, y, size(y,1), s_tmp, size(s_tmp,1), & -! 0.d0, h, size(h,1)) -! -! ! Damp interaction between different spin states -! ! ------------------------------------------------ -! -! do k=1,shift2 -! do l=1,shift2 -! if (dabs(s2(k) - s2(l)) > 1.d0) then -! h(k,l) = h(k,l)*(max(0.d0,1.d0 - dabs(s2(k) - s2(l)))) -! endif -! enddo -! enddo -! -! ! Rotate back H -! ! ------------- -! -! call dgemm('N','T',shift2,shift2,shift2, & -! 1.d0, h, size(h,1), y, size(y,1), & -! 0.d0, s_tmp, size(s_tmp,1)) -! -! call dgemm('N','N',shift2,shift2,shift2, & -! 1.d0, y, size(y,1), s_tmp, size(s_tmp,1), & -! 0.d0, h, size(h,1)) - - - ! Diagonalize h - ! ------------- - call lapack_diag(lambda,y,h,size(h,1),shift2) - - ! Compute S2 for each eigenvector - ! ------------------------------- - - call dgemm('N','N',shift2,shift2,shift2, & - 1.d0, s_, size(s_,1), y, size(y,1), & - 0.d0, s_tmp, size(s_tmp,1)) - - call dgemm('T','N',shift2,shift2,shift2, & - 1.d0, y, size(y,1), s_tmp, size(s_tmp,1), & - 0.d0, s_, size(s_,1)) - - do k=1,shift2 - s2(k) = s_(k,k) + S_z2_Sz - enddo - - if (s2_eig) then - logical :: state_ok(N_st_diag*davidson_sze_max) - do k=1,shift2 - state_ok(k) = (dabs(s2(k)-expected_s2) < 0.6d0) - enddo - else - do k=1,size(state_ok) - state_ok(k) = .True. - enddo - endif - - do k=1,shift2 - if (.not. state_ok(k)) then - do l=k+1,shift2 - if (state_ok(l)) then - call dswap(shift2, y(1,k), 1, y(1,l), 1) - call dswap(1, s2(k), 1, s2(l), 1) - call dswap(1, lambda(k), 1, lambda(l), 1) - state_ok(k) = .True. - state_ok(l) = .False. - exit - endif - enddo - endif - enddo - - if (state_following) then - - ! Compute overlap with U_in - ! ------------------------- - - integer :: order(N_st_diag) - double precision :: cmax - overlap = -1.d0 - do k=1,shift2 - do i=1,shift2 - overlap(k,i) = dabs(y(k,i)) - enddo - enddo - do k=1,N_st - cmax = -1.d0 - do i=1,N_st - if (overlap(i,k) > cmax) then - cmax = overlap(i,k) - order(k) = i - endif - enddo - do i=1,shift2 - overlap(order(k),i) = -1.d0 - enddo - enddo - overlap = y - do k=1,N_st - l = order(k) - if (k /= l) then - y(1:shift2,k) = overlap(1:shift2,l) - endif - enddo - do k=1,N_st - overlap(k,1) = lambda(k) - overlap(k,2) = s2(k) - enddo - do k=1,N_st - l = order(k) - if (k /= l) then - lambda(k) = overlap(l,1) - s2(k) = overlap(l,2) - endif - enddo - - endif - - - ! Express eigenvectors of h in the determinant basis - ! -------------------------------------------------- - - call dgemm('N','N', sze, N_st_diag, shift2, & - 1.d0, U, size(U,1), y, size(y,1), 0.d0, U(1,shift2+1), size(U,1)) - call dgemm('N','N', sze, N_st_diag, shift2, & - 1.d0, W, size(W,1), y, size(y,1), 0.d0, W(1,shift2+1), size(W,1)) - call dgemm('N','N', sze, N_st_diag, shift2, & - 1.d0, S, size(S,1), y, size(y,1), 0.d0, S(1,shift2+1), size(S,1)) - - ! Compute residual vector - ! ----------------------- - - do k=1,N_st_diag -! if (state_ok(k)) then - do i=1,sze - U(i,shift2+k) = (lambda(k) * U(i,shift2+k) - W(i,shift2+k) ) & - * (1.d0 + s2(k) * U(i,shift2+k) - S(i,shift2+k) - S_z2_Sz & - )/max(H_jj(i) - lambda (k),1.d-2) - enddo -! else -! ! Randomize components with bad -! do i=1,sze-2,2 -! call random_number(r1) -! call random_number(r2) -! r1 = dsqrt(-2.d0*dlog(r1)) -! r2 = dtwo_pi*r2 -! U(i,shift2+k) = r1*dcos(r2) -! U(i+1,shift2+k) = r1*dsin(r2) -! enddo -! do i=sze-2+1,sze -! call random_number(r1) -! call random_number(r2) -! r1 = dsqrt(-2.d0*dlog(r1)) -! r2 = dtwo_pi*r2 -! U(i,shift2+k) = r1*dcos(r2) -! enddo -! endif - - if (k <= N_st) then - residual_norm(k) = u_dot_u(U(1,shift2+k),sze) - to_print(1,k) = lambda(k) + nuclear_repulsion - to_print(2,k) = s2(k) - to_print(3,k) = residual_norm(k) - endif - enddo - - write(iunit,'(X,I3,X,100(X,F16.10,X,F11.6,X,E11.3))') iter, to_print(1:3,1:N_st) - call davidson_converged(lambda,residual_norm,wall,iter,cpu,N_st,converged) - do k=1,N_st - if (residual_norm(k) > 1.e8) then - print *, '' - stop 'Davidson failed' - endif - enddo - if (converged) then - exit - endif - - enddo - - ! Re-contract to u_in - ! ----------- - - call dgemm('N','N', sze, N_st_diag, shift2, & - 1.d0, U, size(U,1), y, size(y,1), 0.d0, u_in, size(u_in,1)) - - enddo - - do k=1,N_st_diag - energies(k) = lambda(k) - enddo - - write_buffer = '===== ' - do i=1,N_st - write_buffer = trim(write_buffer)//' ================ =========== ===========' - enddo - write(iunit,'(A)') trim(write_buffer) - write(iunit,'(A)') '' - call write_time(iunit) - - deallocate ( & - W, residual_norm, & - U, overlap, & - c, S, & - h, & - y, s_, s_tmp, & - lambda & - ) -end - - -subroutine H_S2_u_0_mrcc_nstates(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,istate_in,N_st,sze) - use bitmasks - implicit none - BEGIN_DOC - ! Computes v_0 = H|u_0> and s_0 = S^2 |u_0> - ! - ! n : number of determinants - ! - ! H_jj : array of - ! - ! S2_jj : array of - END_DOC - integer, intent(in) :: N_st,n,Nint, sze, istate_in - double precision, intent(out) :: v_0(sze,N_st), s_0(sze,N_st) - double precision, intent(in) :: u_0(sze,N_st) - double precision, intent(in) :: H_jj(n), S2_jj(n) - integer(bit_kind),intent(in) :: keys_tmp(Nint,2,n) - double precision :: hij,s2 - double precision, allocatable :: vt(:,:), ut(:,:), st(:,:) - integer :: i,j,k,l, jj,ii - integer :: i0, j0 - - integer, allocatable :: shortcut(:,:), sort_idx(:,:) - integer(bit_kind), allocatable :: sorted(:,:,:), version(:,:,:) - integer(bit_kind) :: sorted_i(Nint) - - integer :: sh, sh2, ni, exa, ext, org_i, org_j, endi, istate - - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: vt, ut - - ASSERT (Nint > 0) - ASSERT (Nint == N_int) - ASSERT (n>0) - PROVIDE ref_bitmask_energy - - allocate (shortcut(0:n+1,2), sort_idx(n,2), sorted(Nint,n,2), version(Nint,n,2)) - allocate(ut(N_st,n)) - - v_0 = 0.d0 - s_0 = 0.d0 - - do i=1,n - do istate=1,N_st - ut(istate,i) = u_0(i,istate) - enddo - enddo - - call sort_dets_ab_v(keys_tmp, sorted(1,1,1), sort_idx(1,1), shortcut(0,1), version(1,1,1), n, Nint) - call sort_dets_ba_v(keys_tmp, sorted(1,1,2), sort_idx(1,2), shortcut(0,2), version(1,1,2), n, Nint) - - PROVIDE delta_ij_s2 - !$OMP PARALLEL DEFAULT(NONE) & - !$OMP PRIVATE(i,hij,s2,j,k,jj,vt,st,ii,sh,sh2,ni,exa,ext,org_i,org_j,endi,sorted_i,istate)& - !$OMP SHARED(n,keys_tmp,ut,Nint,v_0,s_0,sorted,shortcut,sort_idx,version,N_st, & - !$OMP N_det_ref, idx_ref, N_det_non_ref, idx_non_ref, delta_ij, delta_ij_s2,istate_in) - allocate(vt(N_st,n),st(N_st,n)) - Vt = 0.d0 - St = 0.d0 - - !$OMP DO SCHEDULE(guided) - do sh=1,shortcut(0,1) - do sh2=sh,shortcut(0,1) - exa = 0 - do ni=1,Nint - exa = exa + popcnt(xor(version(ni,sh,1), version(ni,sh2,1))) - end do - if(exa > 2) then - cycle - end if - - do i=shortcut(sh,1),shortcut(sh+1,1)-1 - org_i = sort_idx(i,1) - if(sh==sh2) then - endi = i-1 - else - endi = shortcut(sh2+1,1)-1 - end if - do ni=1,Nint - sorted_i(ni) = sorted(ni,i,1) - enddo - - do j=shortcut(sh2,1),endi - org_j = sort_idx(j,1) - ext = exa - do ni=1,Nint - ext = ext + popcnt(xor(sorted_i(ni), sorted(ni,j,1))) - end do - if(ext <= 4) then - call i_h_j (keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),nint,hij) - call get_s2(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),nint,s2) - do istate=1,n_st - vt (istate,org_i) = vt (istate,org_i) + hij*ut(istate,org_j) - vt (istate,org_j) = vt (istate,org_j) + hij*ut(istate,org_i) - st (istate,org_i) = st (istate,org_i) + s2*ut(istate,org_j) - st (istate,org_j) = st (istate,org_j) + s2*ut(istate,org_i) - enddo - endif - enddo - enddo - enddo - enddo - !$OMP END DO - - !$OMP DO SCHEDULE(guided) - do sh=1,shortcut(0,2) - do i=shortcut(sh,2),shortcut(sh+1,2)-1 - org_i = sort_idx(i,2) - do j=shortcut(sh,2),i-1 - org_j = sort_idx(j,2) - ext = 0 - do ni=1,Nint - ext = ext + popcnt(xor(sorted(ni,i,2), sorted(ni,j,2))) - end do - if(ext == 4) then - call i_h_j (keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),nint,hij) - call get_s2(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),nint,s2) - do istate=1,n_st - vt (istate,org_i) = vt (istate,org_i) + hij*ut(istate,org_j) - vt (istate,org_j) = vt (istate,org_j) + hij*ut(istate,org_i) - st (istate,org_i) = st (istate,org_i) + s2*ut(istate,org_j) - st (istate,org_j) = st (istate,org_j) + s2*ut(istate,org_i) - enddo - end if - end do - end do - enddo - !$OMP END DO - -! -------------------------- -! Begin Specific to dressing -! -------------------------- - -!TODO : DRESSING 1 column - - !$OMP DO - do ii=1,n_det_ref - i = idx_ref(ii) - do jj = 1, n_det_non_ref - j = idx_non_ref(jj) - do istate=1,N_st - vt (istate,i) = vt (istate,i) + delta_ij(istate_in,jj,ii)*ut(istate,j) - vt (istate,j) = vt (istate,j) + delta_ij(istate_in,jj,ii)*ut(istate,i) - st (istate,i) = st (istate,i) + delta_ij_s2(istate_in,jj,ii)*ut(istate,j) - st (istate,j) = st (istate,j) + delta_ij_s2(istate_in,jj,ii)*ut(istate,i) - enddo - enddo - enddo - !$OMP END DO - -! ------------------------ -! End Specific to dressing -! ------------------------ - - do istate=1,N_st - do i=n,1,-1 - !$OMP ATOMIC - v_0(i,istate) = v_0(i,istate) + vt(istate,i) - !$OMP ATOMIC - s_0(i,istate) = s_0(i,istate) + st(istate,i) - enddo - enddo - - deallocate(vt,st) - !$OMP END PARALLEL - - do istate=1,N_st - do i=1,n - v_0(i,istate) = v_0(i,istate) + H_jj(i) * u_0(i,istate) - s_0(i,istate) = s_0(i,istate) + s2_jj(i)* u_0(i,istate) - enddo - enddo - deallocate (shortcut, sort_idx, sorted, version, ut) -end - diff --git a/plugins/MRCC_Utils/mrcc_dress.irp.f b/plugins/MRCC_Utils/mrcc_dress.irp.f index 5c2f5efc..41041288 100644 --- a/plugins/MRCC_Utils/mrcc_dress.irp.f +++ b/plugins/MRCC_Utils/mrcc_dress.irp.f @@ -14,14 +14,13 @@ BEGIN_PROVIDER [ integer(omp_lock_kind), psi_ref_lock, (psi_det_size) ] END_PROVIDER -subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_generator,n_selected,det_buffer,Nint,iproc,key_mask) +subroutine mrcc_dress(delta_ij_, Nstates, Ndet_non_ref, Ndet_ref,i_generator,n_selected,det_buffer,Nint,iproc,key_mask) use bitmasks implicit none integer, intent(in) :: i_generator,n_selected, Nint, iproc integer, intent(in) :: Nstates, Ndet_ref, Ndet_non_ref double precision, intent(inout) :: delta_ij_(Nstates,Ndet_non_ref,Ndet_ref) - double precision, intent(inout) :: delta_ii_(Nstates,Ndet_ref) integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected) integer :: i,j,k,l,m @@ -265,10 +264,8 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge do l_sd=1,idx_alpha(0) k_sd = idx_alpha(l_sd) delta_ij_(i_state,k_sd,i_I) = delta_ij_(i_state,k_sd,i_I) + dIa_hla(i_state,k_sd) - delta_ii_(i_state,i_I) = delta_ii_(i_state,i_I) - dIa_hla(i_state,k_sd) * ci_inv(i_state) * psi_non_ref_coef_transp(i_state,k_sd) enddo else - !delta_ii_(i_state,i_I) = 0.d0 do l_sd=1,idx_alpha(0) k_sd = idx_alpha(l_sd) delta_ij_(i_state,k_sd,i_I) = delta_ij_(i_state,k_sd,i_I) + 0.5d0 * dIa_hla(i_state,k_sd) diff --git a/plugins/MRCC_Utils/mrcc_utils.irp.f b/plugins/MRCC_Utils/mrcc_utils.irp.f index 6609790b..ad653c8c 100644 --- a/plugins/MRCC_Utils/mrcc_utils.irp.f +++ b/plugins/MRCC_Utils/mrcc_utils.irp.f @@ -139,210 +139,6 @@ BEGIN_PROVIDER [ double precision, hij_mrcc, (N_det_non_ref,N_det_ref) ] END_PROVIDER -BEGIN_PROVIDER [ double precision, h_matrix_dressed, (N_det,N_det,N_states) ] - implicit none - BEGIN_DOC - ! Dressed H with Delta_ij - END_DOC - integer :: i, j,istate,ii,jj - do istate = 1,N_states - do j=1,N_det - do i=1,N_det - h_matrix_dressed(i,j,istate) = h_matrix_all_dets(i,j) - enddo - enddo - do ii = 1, N_det_ref - i =idx_ref(ii) - h_matrix_dressed(i,i,istate) += delta_ii(istate,ii) - do jj = 1, N_det_non_ref - j =idx_non_ref(jj) - h_matrix_dressed(i,j,istate) += delta_ij(istate,jj,ii) - h_matrix_dressed(j,i,istate) += delta_ij(istate,jj,ii) - enddo - enddo - enddo -END_PROVIDER - - - BEGIN_PROVIDER [ double precision, CI_electronic_energy_dressed, (N_states_diag) ] -&BEGIN_PROVIDER [ double precision, CI_eigenvectors_dressed, (N_det,N_states_diag) ] -&BEGIN_PROVIDER [ double precision, CI_eigenvectors_s2_dressed, (N_states_diag) ] - implicit none - BEGIN_DOC - ! Eigenvectors/values of the dressed CI matrix - END_DOC - double precision :: ovrlp,u_dot_v - integer :: i_good_state - integer, allocatable :: index_good_state_array(:) - logical, allocatable :: good_state_array(:) - double precision, allocatable :: s2_values_tmp(:) - integer :: i_other_state - double precision, allocatable :: eigenvectors(:,:), eigenvalues(:) - integer :: i_state - double precision :: e_0 - integer :: i,j,k - double precision, allocatable :: s2_eigvalues(:) - double precision, allocatable :: e_array(:) - integer, allocatable :: iorder(:) - - integer :: mrcc_state - - do j=1,min(N_states,N_det) - do i=1,N_det - CI_eigenvectors_dressed(i,j) = psi_coef(i,j) - enddo - enddo - - if (diag_algorithm == "Davidson") then - - allocate (eigenvectors(size(CI_eigenvectors_dressed,1),size(CI_eigenvectors_dressed,2)),& - eigenvalues(size(CI_electronic_energy_dressed,1))) - do j=1,min(N_states,N_det) - do i=1,N_det - eigenvectors(i,j) = psi_coef(i,j) - enddo - enddo - do mrcc_state=1,N_states - do j=mrcc_state,min(N_states,N_det) - do i=1,N_det - eigenvectors(i,j) = psi_coef(i,j) - enddo - enddo - call davidson_diag_mrcc_HS2(psi_det,eigenvectors, & - size(eigenvectors,1), & - eigenvalues,N_det,N_states,N_states_diag,N_int, & - 6,mrcc_state) - CI_eigenvectors_dressed(1:N_det,mrcc_state) = eigenvectors(1:N_det,mrcc_state) - CI_electronic_energy_dressed(mrcc_state) = eigenvalues(mrcc_state) - enddo - do k=N_states+1,N_states_diag - CI_eigenvectors_dressed(1:N_det,k) = eigenvectors(1:N_det,k) - CI_electronic_energy_dressed(k) = eigenvalues(k) - enddo - call u_0_S2_u_0(CI_eigenvectors_s2_dressed,CI_eigenvectors_dressed,N_det,psi_det,N_int,& - N_states_diag,size(CI_eigenvectors_dressed,1)) - - deallocate (eigenvectors,eigenvalues) - - else if (diag_algorithm == "Lapack") then - - allocate (eigenvectors(size(H_matrix_dressed,1),N_det)) - allocate (eigenvalues(N_det)) - call lapack_diag(eigenvalues,eigenvectors, & - H_matrix_dressed,size(H_matrix_dressed,1),N_det) - CI_electronic_energy_dressed(:) = 0.d0 - if (s2_eig) then - i_state = 0 - allocate (s2_eigvalues(N_det)) - allocate(index_good_state_array(N_det),good_state_array(N_det)) - good_state_array = .False. - call u_0_S2_u_0(s2_eigvalues,eigenvectors,N_det,psi_det,N_int,& - N_det,size(eigenvectors,1)) - do j=1,N_det - ! Select at least n_states states with S^2 values closed to "expected_s2" - if(dabs(s2_eigvalues(j)-expected_s2).le.0.5d0)then - i_state += 1 - index_good_state_array(i_state) = j - good_state_array(j) = .True. - endif - if (i_state==N_states) then - exit - endif - enddo - if (i_state /= 0) then - ! Fill the first "i_state" states that have a correct S^2 value - do j = 1, i_state - do i=1,N_det - CI_eigenvectors_dressed(i,j) = eigenvectors(i,index_good_state_array(j)) - enddo - CI_electronic_energy_dressed(j) = eigenvalues(index_good_state_array(j)) - CI_eigenvectors_s2_dressed(j) = s2_eigvalues(index_good_state_array(j)) - enddo - i_other_state = 0 - do j = 1, N_det - if(good_state_array(j))cycle - i_other_state +=1 - if(i_state+i_other_state.gt.n_states_diag)then - exit - endif - do i=1,N_det - CI_eigenvectors_dressed(i,i_state+i_other_state) = eigenvectors(i,j) - enddo - CI_electronic_energy_dressed(i_state+i_other_state) = eigenvalues(j) - CI_eigenvectors_s2_dressed(i_state+i_other_state) = s2_eigvalues(i_state+i_other_state) - enddo - - else - print*,'' - print*,'!!!!!!!! WARNING !!!!!!!!!' - print*,' Within the ',N_det,'determinants selected' - print*,' and the ',N_states_diag,'states requested' - print*,' We did not find any state with S^2 values close to ',expected_s2 - print*,' We will then set the first N_states eigenvectors of the H matrix' - print*,' as the CI_eigenvectors_dressed' - print*,' You should consider more states and maybe ask for s2_eig to be .True. or just enlarge the CI space' - print*,'' - do j=1,min(N_states_diag,N_det) - do i=1,N_det - CI_eigenvectors_dressed(i,j) = eigenvectors(i,j) - enddo - CI_electronic_energy_dressed(j) = eigenvalues(j) - CI_eigenvectors_s2_dressed(j) = s2_eigvalues(j) - enddo - endif - deallocate(index_good_state_array,good_state_array) - deallocate(s2_eigvalues) - else - call u_0_S2_u_0(CI_eigenvectors_s2_dressed,eigenvectors,N_det,psi_det,N_int,& - min(N_det,N_states_diag),size(eigenvectors,1)) - ! Select the "N_states_diag" states of lowest energy - do j=1,min(N_det,N_states_diag) - do i=1,N_det - CI_eigenvectors_dressed(i,j) = eigenvectors(i,j) - enddo - CI_electronic_energy_dressed(j) = eigenvalues(j) - enddo - endif - deallocate(eigenvectors,eigenvalues) - endif - -END_PROVIDER - -BEGIN_PROVIDER [ double precision, CI_energy_dressed, (N_states_diag) ] - implicit none - BEGIN_DOC - ! N_states lowest eigenvalues of the dressed CI matrix - END_DOC - - integer :: j - character*(8) :: st - call write_time(6) - do j=1,min(N_det,N_states) - write(st,'(I4)') j - CI_energy_dressed(j) = CI_electronic_energy_dressed(j) + nuclear_repulsion - call write_double(6,CI_energy_dressed(j),'Energy of state '//trim(st)) - call write_double(6,CI_eigenvectors_s2_dressed(j),'S^2 of state '//trim(st)) - enddo - -END_PROVIDER - -subroutine diagonalize_CI_dressed(lambda) - implicit none - BEGIN_DOC -! Replace the coefficients of the CI states by the coefficients of the -! eigenstates of the CI matrix - END_DOC - double precision, intent(in) :: lambda - integer :: i,j - do j=1,N_states - do i=1,N_det - psi_coef(i,j) = lambda * CI_eigenvectors_dressed(i,j) + (1.d0 - lambda) * psi_coef(i,j) - enddo - call normalize(psi_coef(1,j), N_det) - enddo - SOFT_TOUCH psi_coef - -end logical function is_generable(det1, det2, Nint) diff --git a/plugins/MRCC_Utils/multi_state.irp.f b/plugins/MRCC_Utils/multi_state.irp.f deleted file mode 100644 index b4a2a3cb..00000000 --- a/plugins/MRCC_Utils/multi_state.irp.f +++ /dev/null @@ -1,101 +0,0 @@ -subroutine multi_state(CI_electronic_energy_dressed_,CI_eigenvectors_dressed_,LDA) - implicit none - BEGIN_DOC - ! Multi-state mixing - END_DOC - integer, intent(in) :: LDA - double precision, intent(inout) :: CI_electronic_energy_dressed_(N_states) - double precision, intent(inout) :: CI_eigenvectors_dressed_(LDA,N_states) - double precision, allocatable :: h(:,:,:), s(:,:), Psi(:,:), H_Psi(:,:,:), H_jj(:) - - allocate( h(N_states,N_states,0:N_states), s(N_states,N_states) ) - allocate( Psi(LDA,N_states), H_Psi(LDA,N_states,0:N_states) ) - allocate (H_jj(LDA) ) - -! e_0(i) = u_dot_v(v_0(1,i),u_0(1,i),n)/u_dot_u(u_0(1,i),n) - - integer :: i,j,k,istate - double precision :: U(N_states,N_states), Vt(N_states,N_states), D(N_states) - double precision, external :: diag_H_mat_elem - do istate=1,N_states - do i=1,N_det - H_jj(i) = diag_H_mat_elem(psi_det(1,1,i),N_int) - enddo - - do i=1,N_det_ref - H_jj(idx_ref(i)) += delta_ii(istate,i) - enddo - - do k=1,N_states - do i=1,N_det - Psi(i,k) = CI_eigenvectors_dressed_(i,k) - enddo - enddo - call H_u_0_mrcc_nstates(H_Psi(1,1,istate),Psi,H_jj,N_det,psi_det,N_int,istate,N_states,LDA) - - do k=1,N_states - do i=1,N_states - double precision, external :: u_dot_v - h(i,k,istate) = u_dot_v(Psi(1,i), H_Psi(1,k,istate), N_det) - enddo - enddo - enddo - - do k=1,N_states - do i=1,N_states - s(i,k) = u_dot_v(Psi(1,i), Psi(1,k), N_det) - enddo - enddo - - print *, s(:,:) - print *, '' - - h(:,:,0) = h(:,:,1) - do istate=2,N_states - U(:,:) = h(:,:,0) - call dgemm('N','N',N_states,N_states,N_states,1.d0,& - U, size(U,1), h(1,1,istate), size(h,1), 0.d0, & - h(1,1,0), size(Vt,1)) - enddo - - call svd(h(1,1,0), size(h,1), U, size(U,1), D, Vt, size(Vt,1), N_states, N_states) - do k=1,N_states - D(k) = D(k)**(1./dble(N_states)) - if (D(k) > 0.d0) then - D(k) = -D(k) - endif - enddo - - do j=1,N_states - do i=1,N_states - h(i,j,0) = 0.d0 - do k=1,N_states - h(i,j,0) += U(i,k) * D(k) * Vt(k,j) - enddo - enddo - enddo - - print *, h(:,:,0) - print *,'' - - integer :: LWORK, INFO - double precision, allocatable :: WORK(:) - LWORK=3*N_states - allocate (WORK(LWORK)) - call dsygv(1, 'V', 'U', N_states, h(1,1,0), size(h,1), s, size(s,1), D, WORK, LWORK, INFO) - deallocate(WORK) - - do j=1,N_states - do i=1,N_det - CI_eigenvectors_dressed_(i,j) = 0.d0 - do k=1,N_states - CI_eigenvectors_dressed_(i,j) += Psi(i,k) * h(k,j,0) - enddo - enddo - CI_electronic_energy_dressed_(j) = D(j) - enddo - - - deallocate (h,s, H_jj) - deallocate( Psi, H_Psi ) -end diff --git a/plugins/Properties/NEEDED_CHILDREN_MODULES b/plugins/Properties/NEEDED_CHILDREN_MODULES index 34de8ddb..2a87d1c1 100644 --- a/plugins/Properties/NEEDED_CHILDREN_MODULES +++ b/plugins/Properties/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -Determinants Davidson +Determinants DavidsonUndressed diff --git a/plugins/Psiref_CAS/NEEDED_CHILDREN_MODULES b/plugins/Psiref_CAS/NEEDED_CHILDREN_MODULES index 107c1643..22828878 100644 --- a/plugins/Psiref_CAS/NEEDED_CHILDREN_MODULES +++ b/plugins/Psiref_CAS/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -Psiref_Utils Davidson +Psiref_Utils DavidsonUndressed diff --git a/plugins/UndressedMethod/NEEDED_CHILDREN_MODULES b/plugins/UndressedMethod/NEEDED_CHILDREN_MODULES new file mode 100644 index 00000000..8b137891 --- /dev/null +++ b/plugins/UndressedMethod/NEEDED_CHILDREN_MODULES @@ -0,0 +1 @@ + diff --git a/plugins/UndressedMethod/README.rst b/plugins/UndressedMethod/README.rst new file mode 100644 index 00000000..1b739e5c --- /dev/null +++ b/plugins/UndressedMethod/README.rst @@ -0,0 +1,14 @@ +=============== +UndressedMethod +=============== + +Defines a null dressing vector + +Needed Modules +============== +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. +Documentation +============= +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. diff --git a/plugins/UndressedMethod/null_dressing_vector.irp.f b/plugins/UndressedMethod/null_dressing_vector.irp.f new file mode 100644 index 00000000..faffe964 --- /dev/null +++ b/plugins/UndressedMethod/null_dressing_vector.irp.f @@ -0,0 +1,10 @@ + BEGIN_PROVIDER [ double precision, dressing_column_h, (N_det,N_states) ] +&BEGIN_PROVIDER [ double precision, dressing_column_s, (N_det,N_states) ] + implicit none + BEGIN_DOC + ! Null dressing vectors + END_DOC + dressing_column_h(:,:) = 0.d0 + dressing_column_s(:,:) = 0.d0 +END_PROVIDER + diff --git a/plugins/mrcepa0/NEEDED_CHILDREN_MODULES b/plugins/mrcepa0/NEEDED_CHILDREN_MODULES index 8b6c5a18..fe8255d1 100644 --- a/plugins/mrcepa0/NEEDED_CHILDREN_MODULES +++ b/plugins/mrcepa0/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -Perturbation Selectors_full Generators_full Psiref_CAS MRCC_Utils ZMQ +Perturbation Selectors_full Generators_full Psiref_CAS MRCC_Utils ZMQ diff --git a/plugins/mrcepa0/dressing.irp.f b/plugins/mrcepa0/dressing.irp.f index 727bdba7..951c8c4c 100644 --- a/plugins/mrcepa0/dressing.irp.f +++ b/plugins/mrcepa0/dressing.irp.f @@ -74,10 +74,8 @@ BEGIN_PROVIDER [ double precision, mrcc_norm_acc, (0:N_det_non_ref, N_states) ] END_PROVIDER - BEGIN_PROVIDER [ double precision, delta_ij_mrcc_sto,(N_states,N_det_non_ref,N_det_ref) ] -&BEGIN_PROVIDER [ double precision, delta_ii_mrcc_sto, (N_states, N_det_ref) ] -&BEGIN_PROVIDER [ double precision, delta_ij_s2_mrcc_sto, (N_states,N_det_non_ref,N_det_ref) ] -&BEGIN_PROVIDER [ double precision, delta_ii_s2_mrcc_sto, (N_states, N_det_ref) ] + BEGIN_PROVIDER [ double precision, delta_ij_mrcc_sto,(N_states,N_det_non_ref) ] +&BEGIN_PROVIDER [ double precision, delta_ij_s2_mrcc_sto, (N_states,N_det_non_ref) ] use bitmasks implicit none integer :: gen, h, p, n, t, i, j, h1, h2, p1, p2, s1, s2, iproc @@ -94,10 +92,8 @@ END_PROVIDER read(*,*) n_in_teeth !n_in_teeth = 2 in_teeth_step = 1d0 / dfloat(n_in_teeth) -!double precision :: delta_ij_mrcc_tmp,(N_states,N_det_non_ref,N_det_ref) ] - !double precision :: delta_ii_mrcc_tmp, (N_states,N_det_ref) ] - !double precision :: delta_ij_s2_mrcc_tmp(N_states,N_det_non_ref,N_det_ref) - !double precision :: delta_ii_s2_mrcc_tmp(N_states, N_det_ref) + !double precision :: delta_ij_mrcc_tmp,(N_states,N_det_non_ref) + !double precision :: delta_ij_s2_mrcc_tmp(N_states,N_det_non_ref) coefs = 0d0 coefs(:mrcc_teeth(1,1)-1) = 1d0 @@ -144,15 +140,13 @@ END_PROVIDER delta_ij_mrcc_sto = 0d0 - delta_ii_mrcc_sto = 0d0 delta_ij_s2_mrcc_sto = 0d0 - delta_ii_s2_mrcc_sto = 0d0 PROVIDE dij provide hh_shortcut psi_det_size! lambda_mrcc !$OMP PARALLEL DO default(none) schedule(dynamic) & !$OMP shared(psi_ref, psi_non_ref, hh_exists, pp_exists, N_int, hh_shortcut) & - !$OMP shared(N_det_generators, coefs,N_det_non_ref, N_det_ref, delta_ii_mrcc_sto, delta_ij_mrcc_sto) & - !$OMP shared(contrib,psi_det_generators, delta_ii_s2_mrcc_sto, delta_ij_s2_mrcc_sto) & + !$OMP shared(N_det_generators, coefs,N_det_non_ref, delta_ij_mrcc_sto) & + !$OMP shared(contrib,psi_det_generators, delta_ij_s2_mrcc_sto) & !$OMP private(i,j,curnorm,myCoef, h, n, mask, omask, buf, ok, iproc) do gen= 1,N_det_generators if(coefs(gen) == 0d0) cycle @@ -174,8 +168,8 @@ END_PROVIDER end do n = n - 1 if(n /= 0) then - call mrcc_part_dress(delta_ij_mrcc_sto, delta_ii_mrcc_sto, delta_ij_s2_mrcc_sto, & - delta_ii_s2_mrcc_sto, gen,n,buf,N_int,omask,myCoef,contrib) + call mrcc_part_dress(delta_ij_mrcc_sto, delta_ij_s2_mrcc_sto, & + gen,n,buf,N_int,omask,myCoef,contrib) endif end do deallocate(buf) @@ -185,21 +179,17 @@ END_PROVIDER curnorm = 0d0 - do i=1,N_det_ref do j=1,N_det_non_ref - curnorm += delta_ij_mrcc_sto(1, j, i)**2 + curnorm += delta_ij_mrcc_sto(1,j)*delta_ij_mrcc_sto(1,j) end do - end do - print *, "NORM DELTA ", curnorm**0.5d0 + print *, "NORM DELTA ", dsqrt(curnorm) END_PROVIDER - BEGIN_PROVIDER [ double precision, delta_ij_cancel, (N_states,N_det_non_ref,N_det_ref) ] -&BEGIN_PROVIDER [ double precision, delta_ii_cancel, (N_states, N_det_ref) ] -&BEGIN_PROVIDER [ double precision, delta_ij_s2_cancel, (N_states,N_det_non_ref,N_det_ref) ] -&BEGIN_PROVIDER [ double precision, delta_ii_s2_cancel, (N_states, N_det_ref) ] + BEGIN_PROVIDER [ double precision, delta_ij_cancel, (N_states,N_det_non_ref) ] +&BEGIN_PROVIDER [ double precision, delta_ij_s2_cancel, (N_states,N_det_non_ref) ] use bitmasks implicit none @@ -216,15 +206,19 @@ END_PROVIDER integer(bit_kind) :: det_tmp(N_int, 2), det_tmp2(N_int, 2),inac, virt integer, external :: get_index_in_psi_det_sorted_bit, searchDet,detCmp logical, external :: is_in_wavefunction + double precision :: c0(N_states) provide dij delta_ij_cancel = 0d0 - delta_ii_cancel = 0d0 + + do i_state = 1, N_states + c0(i_state) = 1.d0/psi_coef(dressed_column_idx(i_state),i_state) + enddo do i=1,N_det_ref !$OMP PARALLEL DO default(shared) private(kk, k, blok, exc_Ik,det_tmp2,ok,deg,phase_Ik, l,ll) & - !$OMP private(contrib, contrib_s2, i_state) + !$OMP private(contrib, contrib_s2, i_state, c0) do kk = 1, nlink(i) k = det_cepa0_idx(linked(kk, i)) blok = blokMwen(kk, i) @@ -244,21 +238,10 @@ END_PROVIDER do i_state = 1, N_states contrib = (dij(j, l, i_state) - dij(i, k, i_state)) * delta_cas(i,j,i_state)! * Hla *phase_ia * phase_ik contrib_s2 = dij(j, l, i_state) - dij(i, k, i_state)! * Sla*phase_ia * phase_ik - if(dabs(psi_ref_coef(i,i_state)).ge.1.d-3) then - !$OMP ATOMIC - delta_ij_cancel(i_state,l,i) += contrib - !$OMP ATOMIC - delta_ij_s2_cancel(i_state,l,i) += contrib_s2 - !$OMP ATOMIC - delta_ii_cancel(i_state,i) -= contrib / psi_ref_coef(i, i_state) * psi_non_ref_coef(l,i_state) - !$OMP ATOMIC - delta_ii_s2_cancel(i_state,i) -= contrib_s2 / psi_ref_coef(i, i_state) * psi_non_ref_coef(l,i_state) - else - !$OMP ATOMIC - delta_ij_cancel(i_state,l,i) += contrib * 0.5d0 - !$OMP ATOMIC - delta_ij_s2_cancel(i_state,l,i) += contrib_s2 * 0.5d0 - endif + !$OMP ATOMIC + delta_ij_cancel(i_state,l) += contrib * psi_ref_coef(i,i_state) * c0(i_state) + !$OMP ATOMIC + delta_ij_s2_cancel(i_state,l) += contrib_s2* psi_ref_coef(i,i_state) * c0(i_state) end do end do end do @@ -268,10 +251,8 @@ END_PROVIDER END_PROVIDER - BEGIN_PROVIDER [ double precision, delta_ij_mrcc, (N_states,N_det_non_ref,N_det_ref) ] -&BEGIN_PROVIDER [ double precision, delta_ii_mrcc, (N_states, N_det_ref) ] -&BEGIN_PROVIDER [ double precision, delta_ij_s2_mrcc, (N_states,N_det_non_ref,N_det_ref) ] -&BEGIN_PROVIDER [ double precision, delta_ii_s2_mrcc, (N_states, N_det_ref) ] + BEGIN_PROVIDER [ double precision, delta_ij_mrcc, (N_states,N_det_non_ref) ] +&BEGIN_PROVIDER [ double precision, delta_ij_s2_mrcc, (N_states,N_det_non_ref) ] use bitmasks implicit none integer :: gen, h, p, n, t, i, h1, h2, p1, p2, s1, s2, iproc @@ -286,14 +267,12 @@ END_PROVIDER contrib = 0d0 delta_ij_mrcc = 0d0 - delta_ii_mrcc = 0d0 delta_ij_s2_mrcc = 0d0 - delta_ii_s2_mrcc = 0d0 !$OMP PARALLEL DO default(none) schedule(dynamic) & !$OMP shared(contrib,psi_det_generators, N_det_generators, hh_exists, pp_exists, N_int, hh_shortcut) & - !$OMP shared(N_det_non_ref, N_det_ref, delta_ii_mrcc, delta_ij_mrcc, delta_ii_s2_mrcc, delta_ij_s2_mrcc) & + !$OMP shared(N_det_non_ref, N_det_ref, delta_ij_mrcc, delta_ij_s2_mrcc) & !$OMP private(h, n, mask, omask, buf, ok, iproc) do gen= 1, N_det_generators allocate(buf(N_int, 2, N_det_non_ref)) @@ -313,7 +292,7 @@ END_PROVIDER n = n - 1 if(n /= 0) then - call mrcc_part_dress(delta_ij_mrcc, delta_ii_mrcc, delta_ij_s2_mrcc, delta_ii_s2_mrcc, gen,n,buf,N_int,omask,1d0,contrib) + call mrcc_part_dress(delta_ij_mrcc, delta_ij_s2_mrcc, gen,n,buf,N_int,omask,1d0,contrib) endif end do @@ -324,20 +303,18 @@ END_PROVIDER ! subroutine blit(b1, b2) -! double precision :: b1(N_states,N_det_non_ref,N_det_ref), b2(N_states,N_det_non_ref,N_det_ref) +! double precision :: b1(N_states,N_det_non_ref), b2(N_states,N_det_non_ref) ! b1 = b1 + b2 ! end subroutine -subroutine mrcc_part_dress(delta_ij_, delta_ii_,delta_ij_s2_, delta_ii_s2_,i_generator,n_selected,det_buffer,Nint,key_mask,coef,contrib) +subroutine mrcc_part_dress(delta_ij_, delta_ij_s2_, i_generator,n_selected,det_buffer,Nint,key_mask,coef,contrib) use bitmasks implicit none integer, intent(in) :: i_generator,n_selected, Nint - double precision, intent(inout) :: delta_ij_(N_states,N_det_non_ref,N_det_ref) - double precision, intent(inout) :: delta_ii_(N_states,N_det_ref) - double precision, intent(inout) :: delta_ij_s2_(N_states,N_det_non_ref,N_det_ref) - double precision, intent(inout) :: delta_ii_s2_(N_states,N_det_ref) + double precision, intent(inout) :: delta_ij_(N_states,N_det_non_ref) + double precision, intent(inout) :: delta_ij_s2_(N_states,N_det_non_ref) integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected) integer :: i,j,k,l,m @@ -399,6 +376,11 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,delta_ij_s2_, delta_ii_s2_,i_gen deallocate(microlist, idx_microlist) + double precision :: c0(N_states) + do i_state=1,N_states + c0(i_state) = 1.d0/psi_coef(dressed_column_idx(i_state),i_state) + enddo + allocate (dIa_hla(N_states,N_det_non_ref), dIa_sla(N_states,N_det_non_ref)) ! |I> @@ -436,8 +418,8 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,delta_ij_s2_, delta_ii_s2_,i_gen do i_alpha=1,N_tq - if(key_mask(1,1) /= 0) then - call getMobiles(tq(1,1,i_alpha), key_mask, mobiles, Nint) + if(key_mask(1,1) /= 0) then + call getMobiles(tq(1,1,i_alpha), key_mask, mobiles, Nint) if(N_microlist(mobiles(1)) < N_microlist(mobiles(2))) then smallerlist = mobiles(1) @@ -445,7 +427,7 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,delta_ij_s2_, delta_ii_s2_,i_gen smallerlist = mobiles(2) end if - + do l=0,N_microlist(smallerlist)-1 microlist_zero(:,:,ptr_microlist(1) + l) = microlist(:,:,ptr_microlist(smallerlist) + l) idx_microlist_zero(ptr_microlist(1) + l) = idx_microlist(ptr_microlist(smallerlist) + l) @@ -467,9 +449,9 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,delta_ij_s2_, delta_ii_s2_,i_gen k_sd = idx_alpha(l_sd) call i_h_j(tq(1,1,i_alpha),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,hij_cache(k_sd)) call get_s2(tq(1,1,i_alpha),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,sij_cache(k_sd)) - !if(sij_cache(k_sd) /= 0D0) PRINT *, "SIJ ", sij_cache(k_sd) + !if(sij_cache(k_sd) /= 0D0) PRINT *, "SIJ ", sij_cache(k_sd) enddo - + ! |I> do i_I=1,N_det_ref ! Find triples and quadruple grand parents @@ -484,12 +466,12 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,delta_ij_s2_, delta_ii_s2_,i_gen ! |alpha> do k_sd=1,idx_alpha(0) - + call get_excitation_degree(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(k_sd)),degree,Nint) if (degree > 2) then cycle endif - + ! ! |l> = Exc(k -> alpha) |I> @@ -499,7 +481,7 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,delta_ij_s2_, delta_ii_s2_,i_gen tmp_det(k,1) = psi_ref(k,1,i_I) tmp_det(k,2) = psi_ref(k,2,i_I) enddo - logical :: ok + logical :: ok call apply_excitation(psi_ref(1,1,i_I), exc, tmp_det, ok, Nint) do i_state=1,N_states @@ -510,7 +492,7 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,delta_ij_s2_, delta_ii_s2_,i_gen do i_state=1,N_states dka(i_state) = 0.d0 enddo - + if (ok) then do l_sd=k_sd+1,idx_alpha(0) call get_excitation_degree(tmp_det,psi_non_ref(1,1,idx_alpha(l_sd)),degree,Nint) @@ -522,40 +504,40 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,delta_ij_s2_, delta_ii_s2_,i_gen exit endif enddo - + else if (perturbative_triples) then - ! Linked - - hka = hij_cache(idx_alpha(k_sd)) - if (dabs(hka) > 1.d-12) then - call get_delta_e_dyall_general_mp(psi_ref(1,1,i_I),tq(1,1,i_alpha),Delta_E_inv) - - do i_state=1,N_states - ASSERT (Delta_E_inv(i_state) < 0.d0) - dka(i_state) = hka / Delta_E_inv(i_state) - enddo - endif - + ! Linked + + hka = hij_cache(idx_alpha(k_sd)) + if (dabs(hka) > 1.d-12) then + call get_delta_e_dyall_general_mp(psi_ref(1,1,i_I),tq(1,1,i_alpha),Delta_E_inv) + + do i_state=1,N_states + ASSERT (Delta_E_inv(i_state) < 0.d0) + dka(i_state) = hka / Delta_E_inv(i_state) + enddo + endif + endif - + if (perturbative_triples.and. (degree2 == 1) ) then - call i_h_j(psi_ref(1,1,i_I),tmp_det,Nint,hka) - hka = hij_cache(idx_alpha(k_sd)) - hka - if (dabs(hka) > 1.d-12) then - call get_delta_e_dyall_general_mp(psi_ref(1,1,i_I),tq(1,1,i_alpha),Delta_E_inv) - do i_state=1,N_states - ASSERT (Delta_E_inv(i_state) < 0.d0) - dka(i_state) = hka / Delta_E_inv(i_state) - enddo - endif - + call i_h_j(psi_ref(1,1,i_I),tmp_det,Nint,hka) + hka = hij_cache(idx_alpha(k_sd)) - hka + if (dabs(hka) > 1.d-12) then + call get_delta_e_dyall_general_mp(psi_ref(1,1,i_I),tq(1,1,i_alpha),Delta_E_inv) + do i_state=1,N_states + ASSERT (Delta_E_inv(i_state) < 0.d0) + dka(i_state) = hka / Delta_E_inv(i_state) + enddo + endif + endif - + do i_state=1,N_states dIa(i_state) = dIa(i_state) + dIk(i_state) * dka(i_state) enddo enddo - + do i_state=1,N_states ci_inv(i_state) = psi_ref_coef_inv(i_I,i_state) enddo @@ -569,39 +551,17 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,delta_ij_s2_, delta_ii_s2_,i_gen enddo enddo do i_state=1,N_states - if(dabs(psi_ref_coef(1,i_state)).ge.1.d-3)then - do l_sd=1,idx_alpha(0) - k_sd = idx_alpha(l_sd) - p1 = 1 - hdress = dIa_hla(i_state,k_sd) * psi_ref_coef(i_I,i_state) / psi_ref_coef(p1,i_state) - sdress = dIa_sla(i_state,k_sd) * psi_ref_coef(i_I,i_state) / psi_ref_coef(p1,i_state) - !$OMP ATOMIC - contrib(i_state) += hdress * psi_ref_coef(p1, i_state) * psi_non_ref_coef(k_sd, i_state) - !$OMP ATOMIC - delta_ij_(i_state,k_sd,p1) += hdress - !$OMP ATOMIC - !delta_ii_(i_state,i_I) = delta_ii_(i_state,i_I) - dIa_hla(i_state,k_sd) * ci_inv(i_state) * psi_non_ref_coef_transp(i_state,k_sd) - delta_ii_(i_state,p1) -= hdress / psi_ref_coef(p1,i_state) * psi_non_ref_coef_transp(i_state,k_sd) - !$OMP ATOMIC - delta_ij_s2_(i_state,k_sd,p1) += sdress - !$OMP ATOMIC - !delta_ii_s2_(i_state,i_I) = delta_ii_s2_(i_state,i_I) - dIa_sla(i_state,k_sd) * ci_inv(i_state) * psi_non_ref_coef_transp(i_state,k_sd) - delta_ii_s2_(i_state,p1) -= sdress / psi_ref_coef(p1,i_state) * psi_non_ref_coef_transp(i_state,k_sd) - enddo - else - !stop "dress with coef < 1d-3" - delta_ii_(i_state,1) = 0.d0 - do l_sd=1,idx_alpha(0) - k_sd = idx_alpha(l_sd) - p1 = 1 - hdress = dIa_hla(i_state,k_sd) * psi_ref_coef(i_I,i_state) / psi_ref_coef(p1,i_state) - sdress = dIa_sla(i_state,k_sd) * psi_ref_coef(i_I,i_state) / psi_ref_coef(p1,i_state) - !$OMP ATOMIC - delta_ij_(i_state,k_sd,p1) = delta_ij_(i_state,k_sd,p1) + 0.5d0*hdress - !$OMP ATOMIC - delta_ij_s2_(i_state,k_sd,p1) = delta_ij_s2_(i_state,k_sd,p1) + 0.5d0*sdress - enddo - endif + do l_sd=1,idx_alpha(0) + k_sd = idx_alpha(l_sd) + hdress = dIa_hla(i_state,k_sd) * psi_ref_coef(i_I,i_state) * c0(i_state) + sdress = dIa_sla(i_state,k_sd) * psi_ref_coef(i_I,i_state) * c0(i_state) + !$OMP ATOMIC + contrib(i_state) += hdress * psi_coef(dressed_column_idx(i_state), i_state) * psi_non_ref_coef(k_sd, i_state) + !$OMP ATOMIC + delta_ij_(i_state,k_sd) += hdress + !$OMP ATOMIC + delta_ij_s2_(i_state,k_sd) += sdress + enddo enddo enddo enddo @@ -611,15 +571,13 @@ end -subroutine mrcc_part_dress_1c(delta_ij_, delta_ii_,delta_ij_s2_, delta_ii_s2_,i_generator,n_selected,det_buffer,Nint,key_mask,contrib) +subroutine mrcc_part_dress_1c(delta_ij_, delta_ij_s2_, i_generator,n_selected,det_buffer,Nint,key_mask,contrib) use bitmasks implicit none integer, intent(in) :: i_generator,n_selected, Nint double precision, intent(inout) :: delta_ij_(N_states,N_det_non_ref) - double precision, intent(inout) :: delta_ii_(N_states) double precision, intent(inout) :: delta_ij_s2_(N_states,N_det_non_ref) - double precision, intent(inout) :: delta_ii_s2_(N_states) integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected) integer :: i,j,k,l,m @@ -715,6 +673,11 @@ subroutine mrcc_part_dress_1c(delta_ij_, delta_ii_,delta_ij_s2_, delta_ii_s2_,i_ end if end if + double precision :: c0(N_states) + do i_state=1,N_states + c0(i_state) = 1.d0/psi_coef(dressed_column_idx(i_state),i_state) + enddo + do i_alpha=1,N_tq if(key_mask(1,1) /= 0) then @@ -850,39 +813,17 @@ subroutine mrcc_part_dress_1c(delta_ij_, delta_ii_,delta_ij_s2_, delta_ii_s2_,i_ enddo enddo do i_state=1,N_states - if(dabs(psi_ref_coef(1,i_state)).ge.1.d-3)then do l_sd=1,idx_alpha(0) k_sd = idx_alpha(l_sd) - p1 = 1 - hdress = dIa_hla(i_state,k_sd) * psi_ref_coef(i_I,i_state) / psi_ref_coef(p1,i_state) - sdress = dIa_sla(i_state,k_sd) * psi_ref_coef(i_I,i_state) / psi_ref_coef(p1,i_state) + hdress = dIa_hla(i_state,k_sd) * psi_ref_coef(i_I,i_state) * c0(i_state) + sdress = dIa_sla(i_state,k_sd) * psi_ref_coef(i_I,i_state) * c0(i_state) !$OMP ATOMIC - contrib(i_state) += hdress * psi_ref_coef(p1, i_state) * psi_non_ref_coef(k_sd, i_state) + contrib(i_state) += hdress * psi_ref_coef(dressed_column_idx(i_state), i_state) * psi_non_ref_coef(k_sd, i_state) !$OMP ATOMIC delta_ij_(i_state,k_sd) += hdress !$OMP ATOMIC - !delta_ii_(i_state,i_I) = delta_ii_(i_state,i_I) - dIa_hla(i_state,k_sd) * ci_inv(i_state) * psi_non_ref_coef_transp(i_state,k_sd) - delta_ii_(i_state) -= hdress / psi_ref_coef(p1,i_state) * psi_non_ref_coef_transp(i_state,k_sd) - !$OMP ATOMIC delta_ij_s2_(i_state,k_sd) += sdress - !$OMP ATOMIC - !delta_ii_s2_(i_state,i_I) = delta_ii_s2_(i_state,i_I) - dIa_sla(i_state,k_sd) * ci_inv(i_state) * psi_non_ref_coef_transp(i_state,k_sd) - delta_ii_s2_(i_state) -= sdress / psi_ref_coef(p1,i_state) * psi_non_ref_coef_transp(i_state,k_sd) enddo - else - !stop "dress with coef < 1d-3" - delta_ii_(i_state) = 0.d0 - do l_sd=1,idx_alpha(0) - k_sd = idx_alpha(l_sd) - p1 = 1 - hdress = dIa_hla(i_state,k_sd) * psi_ref_coef(i_I,i_state) / psi_ref_coef(p1,i_state) - sdress = dIa_sla(i_state,k_sd) * psi_ref_coef(i_I,i_state) / psi_ref_coef(p1,i_state) - !$OMP ATOMIC - delta_ij_(i_state,k_sd) = delta_ij_(i_state,k_sd) + 0.5d0*hdress - !$OMP ATOMIC - delta_ij_s2_(i_state,k_sd) = delta_ij_s2_(i_state,k_sd) + 0.5d0*sdress - enddo - endif enddo enddo enddo @@ -900,10 +841,8 @@ end END_PROVIDER - BEGIN_PROVIDER [ double precision, delta_ij_mrcc_zmq, (N_states,N_det_non_ref,N_det_ref) ] -&BEGIN_PROVIDER [ double precision, delta_ii_mrcc_zmq, (N_states, N_det_ref) ] -&BEGIN_PROVIDER [ double precision, delta_ij_s2_mrcc_zmq, (N_states,N_det_non_ref,N_det_ref) ] -&BEGIN_PROVIDER [ double precision, delta_ii_s2_mrcc_zmq, (N_states, N_det_ref) ] + BEGIN_PROVIDER [ double precision, delta_ij_mrcc_zmq, (N_states,N_det_non_ref) ] +&BEGIN_PROVIDER [ double precision, delta_ij_s2_mrcc_zmq, (N_states,N_det_non_ref) ] use bitmasks implicit none @@ -917,9 +856,7 @@ end delta_ij_mrcc_zmq = 0d0 - delta_ii_mrcc_zmq = 0d0 delta_ij_s2_mrcc_zmq = 0d0 - delta_ii_s2_mrcc_zmq = 0d0 !call random_seed() E_CI_before = mrcc_E0_denominator(1) + nuclear_repulsion @@ -935,142 +872,67 @@ end call ZMQ_mrcc(E_CI_before, mrcc, delta_ij_mrcc_zmq, delta_ij_s2_mrcc_zmq, abs(relative_error)) mrcc_previous_E(:) = mrcc_E0_denominator(:) - do i=N_det_non_ref,1,-1 - delta_ii_mrcc_zmq(:,1) -= delta_ij_mrcc_zmq(:, i, 1) / psi_ref_coef(1,1) * psi_non_ref_coef(i, 1) - delta_ii_s2_mrcc_zmq(:,1) -= delta_ij_s2_mrcc_zmq(:, i, 1) / psi_ref_coef(1,1) * psi_non_ref_coef(i, 1) - end do END_PROVIDER - BEGIN_PROVIDER [ double precision, delta_ij, (N_states,N_det_non_ref,N_det_ref) ] -&BEGIN_PROVIDER [ double precision, delta_ii, (N_states, N_det_ref) ] -&BEGIN_PROVIDER [ double precision, delta_ij_s2, (N_states,N_det_non_ref,N_det_ref) ] -&BEGIN_PROVIDER [ double precision, delta_ii_s2, (N_states, N_det_ref) ] + BEGIN_PROVIDER [ double precision, delta_ij, (N_states,N_det_non_ref) ] +&BEGIN_PROVIDER [ double precision, delta_ij_s2, (N_states,N_det_non_ref) ] use bitmasks implicit none integer :: i, j, i_state - !mrmode : 1=mrcepa0, 2=mrsc2 add, 3=mrcc + !mrmode : 1=mrcepa0, 2=mrsc2 add, 3=mrcc, 4=stoch if(mrmode == 4) then - do i = 1, N_det_ref - do i_state = 1, N_states - delta_ii(i_state,i)= delta_ii_mrcc_sto(i_state,i) - delta_ii_s2(i_state,i)= delta_ii_s2_mrcc_sto(i_state,i) - enddo do j = 1, N_det_non_ref do i_state = 1, N_states - delta_ij(i_state,j,i) = delta_ij_mrcc_sto(i_state,j,i) - delta_ij_s2(i_state,j,i) = delta_ij_s2_mrcc_sto(i_state,j,i) + delta_ij(i_state,j) = delta_ij_mrcc_sto(i_state,j) + delta_ij_s2(i_state,j) = delta_ij_s2_mrcc_sto(i_state,j) enddo end do - end do ! else if(mrmode == 10) then -! do i = 1, N_det_ref -! do i_state = 1, N_states -! delta_ii(i_state,i)= delta_ii_mrsc2(i_state,i) -! delta_ii_s2(i_state,i)= delta_ii_s2_mrsc2(i_state,i) -! enddo ! do j = 1, N_det_non_ref ! do i_state = 1, N_states -! delta_ij(i_state,j,i) = delta_ij_mrsc2(i_state,j,i) -! delta_ij_s2(i_state,j,i) = delta_ij_s2_mrsc2(i_state,j,i) +! delta_ij(i_state,j) = delta_ij_mrsc2(i_state,j) +! delta_ij_s2(i_state,j) = delta_ij_s2_mrsc2(i_state,j) ! enddo ! end do -! end do else if(mrmode == 5) then - do i = 1, N_det_ref - do i_state = 1, N_states - delta_ii(i_state,i)= delta_ii_mrcc_zmq(i_state,i) - delta_ii_s2(i_state,i)= delta_ii_s2_mrcc_zmq(i_state,i) - enddo do j = 1, N_det_non_ref do i_state = 1, N_states - delta_ij(i_state,j,i) = delta_ij_mrcc_zmq(i_state,j,i) - delta_ij_s2(i_state,j,i) = delta_ij_s2_mrcc_zmq(i_state,j,i) + delta_ij(i_state,j) = delta_ij_mrcc_zmq(i_state,j) + delta_ij_s2(i_state,j) = delta_ij_s2_mrcc_zmq(i_state,j) enddo end do - end do else if(mrmode == 3) then - do i = 1, N_det_ref - do i_state = 1, N_states - delta_ii(i_state,i)= delta_ii_mrcc(i_state,i) - delta_ii_s2(i_state,i)= delta_ii_s2_mrcc(i_state,i) - enddo do j = 1, N_det_non_ref do i_state = 1, N_states - delta_ij(i_state,j,i) = delta_ij_mrcc(i_state,j,i) - delta_ij_s2(i_state,j,i) = delta_ij_s2_mrcc(i_state,j,i) + delta_ij(i_state,j) = delta_ij_mrcc(i_state,j) + delta_ij_s2(i_state,j) = delta_ij_s2_mrcc(i_state,j) enddo end do - end do - - ! =-=-= BEGIN STATE AVERAGE -! do i = 1, N_det_ref -! delta_ii(:,i)= delta_ii_mrcc(1,i) -! delta_ii_s2(:,i)= delta_ii_s2_mrcc(1,i) -! do i_state = 2, N_states -! delta_ii(:,i) += delta_ii_mrcc(i_state,i) -! delta_ii_s2(:,i) += delta_ii_s2_mrcc(i_state,i) -! enddo -! do j = 1, N_det_non_ref -! delta_ij(:,j,i) = delta_ij_mrcc(1,j,i) -! delta_ij_s2(:,j,i) = delta_ij_s2_mrcc(1,j,i) -! do i_state = 2, N_states -! delta_ij(:,j,i) += delta_ij_mrcc(i_state,j,i) -! delta_ij_s2(:,j,i) += delta_ij_s2_mrcc(i_state,j,i) -! enddo -! end do -! end do -! delta_ij = delta_ij * (1.d0/dble(N_states)) -! delta_ii = delta_ii * (1.d0/dble(N_states)) - ! =-=-= END STATE AVERAGE - ! - ! do i = 1, N_det_ref - ! delta_ii(i_state,i)= delta_mrcepa0_ii(i,i_state) - delta_sub_ii(i,i_state) - ! do j = 1, N_det_non_ref - ! delta_ij(i_state,j,i) = delta_mrcepa0_ij(i,j,i_state) - delta_sub_ij(i,j,i_state) - ! end do - ! end do else if(mrmode == 2) then - do i = 1, N_det_ref - do i_state = 1, N_states - delta_ii(i_state,i)= delta_ii_old(i_state,i) - delta_ii_s2(i_state,i)= delta_ii_s2_old(i_state,i) - enddo do j = 1, N_det_non_ref do i_state = 1, N_states - delta_ij(i_state,j,i) = delta_ij_old(i_state,j,i) - delta_ij_s2(i_state,j,i) = delta_ij_s2_old(i_state,j,i) + delta_ij(i_state,j) = delta_ij_old(i_state,j) + delta_ij_s2(i_state,j) = delta_ij_s2_old(i_state,j) enddo end do - end do else if(mrmode == 1) then - do i = 1, N_det_ref - do i_state = 1, N_states - delta_ii(i_state,i)= delta_mrcepa0_ii(i,i_state) - delta_ii_s2(i_state,i)= delta_mrcepa0_ii_s2(i,i_state) - enddo do j = 1, N_det_non_ref do i_state = 1, N_states - delta_ij(i_state,j,i) = delta_mrcepa0_ij(i,j,i_state) - delta_ij_s2(i_state,j,i) = delta_mrcepa0_ij_s2(i,j,i_state) + delta_ij(i_state,j) = delta_mrcepa0_ij(j,i_state) + delta_ij_s2(i_state,j) = delta_mrcepa0_ij_s2(j,i_state) enddo end do - end do else stop "invalid mrmode" end if !if(mrmode == 2 .or. mrmode == 3) then - ! do i = 1, N_det_ref - ! do i_state = 1, N_states - ! delta_ii(i_state,i) += delta_ii_cancel(i_state,i) - ! enddo ! do j = 1, N_det_non_ref ! do i_state = 1, N_states - ! delta_ij(i_state,j,i) += delta_ij_cancel(i_state,j,i) + ! delta_ij(i_state,j) += delta_ij_cancel(i_state,j) ! enddo ! end do - ! end do !end if END_PROVIDER @@ -1352,10 +1214,8 @@ subroutine getHP(a,h,p,Nint) end subroutine - BEGIN_PROVIDER [ double precision, delta_mrcepa0_ij, (N_det_ref,N_det_non_ref,N_states) ] -&BEGIN_PROVIDER [ double precision, delta_mrcepa0_ii, (N_det_ref,N_states) ] -&BEGIN_PROVIDER [ double precision, delta_mrcepa0_ij_s2, (N_det_ref,N_det_non_ref,N_states) ] -&BEGIN_PROVIDER [ double precision, delta_mrcepa0_ii_s2, (N_det_ref,N_states) ] + BEGIN_PROVIDER [ double precision, delta_mrcepa0_ij, (N_det_non_ref,N_states) ] +&BEGIN_PROVIDER [ double precision, delta_mrcepa0_ij_s2, (N_det_non_ref,N_states) ] use bitmasks implicit none @@ -1363,7 +1223,7 @@ end subroutine integer :: p1,p2,h1,h2,s1,s2, p1_,p2_,h1_,h2_,s1_,s2_, sortRefIdx(N_det_ref) logical :: ok double precision :: phase_iI, phase_Ik, phase_Jl, phase_IJ, phase_al, diI, hIi, hJi, delta_JI, dkI(1), HkI, ci_inv(1), dia_hla(1) - double precision :: contrib, contrib2, contrib_s2, contrib2_s2, HIIi, HJk, wall + double precision :: contrib, contrib_s2, HIIi, HJk, wall integer, dimension(0:2,2,2) :: exc_iI, exc_Ik, exc_IJ integer(bit_kind) :: det_tmp(N_int, 2), made_hole(N_int,2), made_particle(N_int,2), myActive(N_int,2) integer(bit_kind),allocatable :: sortRef(:,:,:) @@ -1385,20 +1245,23 @@ end subroutine idx_sorted_bit(get_index_in_psi_det_sorted_bit(psi_non_ref(1,1,i), N_int)) = i enddo + double precision :: c0(N_states) + do i_state=1,N_states + c0(i_state) = 1.d0/psi_coef(dressed_column_idx(i_state),i_state) + enddo + ! To provide everything contrib = dij(1, 1, 1) - delta_mrcepa0_ii(:,:) = 0d0 - delta_mrcepa0_ij(:,:,:) = 0d0 - delta_mrcepa0_ii_s2(:,:) = 0d0 - delta_mrcepa0_ij_s2(:,:,:) = 0d0 + delta_mrcepa0_ij(:,:) = 0d0 + delta_mrcepa0_ij_s2(:,:) = 0d0 do i_state = 1, N_states - !$OMP PARALLEL DO default(none) schedule(dynamic) shared(delta_mrcepa0_ij, delta_mrcepa0_ii, delta_mrcepa0_ij_s2, delta_mrcepa0_ii_s2) & - !$OMP private(m,i,II,J,k,degree,myActive,made_hole,made_particle,hjk,contrib,contrib2,contrib_s2,contrib2_s2) & + !$OMP PARALLEL DO default(none) schedule(dynamic) shared(delta_mrcepa0_ij, delta_mrcepa0_ij_s2) & + !$OMP private(m,i,II,J,k,degree,myActive,made_hole,made_particle,hjk,contrib,contrib_s2) & !$OMP shared(active_sorb, psi_non_ref, psi_non_ref_coef, psi_ref, psi_ref_coef, cepa0_shortcut, det_cepa0_active) & !$OMP shared(N_det_ref, N_det_non_ref,N_int,det_cepa0_idx,lambda_mrcc,det_ref_active, delta_cas, delta_cas_s2) & - !$OMP shared(notf,i_state, sortRef, sortRefIdx, dij) + !$OMP shared(notf,i_state, sortRef, sortRefIdx, dij,c0) do blok=1,cepa0_shortcut(0) do i=cepa0_shortcut(blok), cepa0_shortcut(blok+1)-1 do II=1,N_det_ref @@ -1438,23 +1301,12 @@ end subroutine !$OMP ATOMIC notf = notf+1 -! call i_h_j(psi_non_ref(1,1,det_cepa0_idx(k)),psi_ref(1,1,J),N_int,HJk) contrib = delta_cas(II, J, i_state)* dij(J, det_cepa0_idx(k), i_state) contrib_s2 = delta_cas_s2(II, J, i_state) * dij(J, det_cepa0_idx(k), i_state) - if(dabs(psi_ref_coef(J,i_state)).ge.1.d-3) then - contrib2 = contrib / psi_ref_coef(J, i_state) * psi_non_ref_coef(det_cepa0_idx(i),i_state) - contrib2_s2 = contrib_s2 / psi_ref_coef(J, i_state) * psi_non_ref_coef(det_cepa0_idx(i),i_state) - !$OMP ATOMIC - delta_mrcepa0_ii(J,i_state) -= contrib2 - delta_mrcepa0_ii_s2(J,i_state) -= contrib2_s2 - else - contrib = contrib * 0.5d0 - contrib_s2 = contrib_s2 * 0.5d0 - end if !$OMP ATOMIC - delta_mrcepa0_ij(J, det_cepa0_idx(i), i_state) += contrib - delta_mrcepa0_ij_s2(J, det_cepa0_idx(i), i_state) += contrib_s2 + delta_mrcepa0_ij(det_cepa0_idx(i), i_state) += contrib * c0(i_state) * psi_ref_coef(J,i_state) + delta_mrcepa0_ij_s2(det_cepa0_idx(i), i_state) += contrib_s2 * c0(i_state) * psi_ref_coef(J,i_state) end do kloop end do @@ -1469,8 +1321,7 @@ end subroutine END_PROVIDER - BEGIN_PROVIDER [ double precision, delta_sub_ij, (N_det_ref,N_det_non_ref,N_states) ] -&BEGIN_PROVIDER [ double precision, delta_sub_ii, (N_det_ref, N_states) ] +BEGIN_PROVIDER [ double precision, delta_sub_ij, (N_det_non_ref,N_states) ] use bitmasks implicit none @@ -1478,7 +1329,7 @@ END_PROVIDER integer :: p1,p2,h1,h2,s1,s2, p1_,p2_,h1_,h2_,s1_,s2_ logical :: ok double precision :: phase_Ji, phase_Ik, phase_Ii - double precision :: contrib, contrib2, delta_IJk, HJk, HIk, HIl + double precision :: contrib, delta_IJk, HJk, HIk, HIl integer, dimension(0:2,2,2) :: exc_Ik, exc_Ji, exc_Ii integer(bit_kind) :: det_tmp(N_int, 2), det_tmp2(N_int, 2) integer, allocatable :: idx_sorted_bit(:) @@ -1492,21 +1343,27 @@ END_PROVIDER do i=1,N_det_non_ref idx_sorted_bit(get_index_in_psi_det_sorted_bit(psi_non_ref(1,1,i), N_int)) = i enddo + + double precision :: c0(N_states) + do i_state=1,N_states + c0(i_state) = 1.d0/psi_coef(dressed_column_idx(i_state),i_state) + enddo + + do i_state = 1, N_states - delta_sub_ij(:,:,:) = 0d0 - delta_sub_ii(:,:) = 0d0 + delta_sub_ij(:,:) = 0d0 provide mo_bielec_integrals_in_map - !$OMP PARALLEL DO default(none) schedule(dynamic,10) shared(delta_sub_ij, delta_sub_ii) & + !$OMP PARALLEL DO default(none) schedule(dynamic,10) shared(delta_sub_ij) & !$OMP private(i, J, k, degree, degree2, l, deg, ni) & !$OMP private(p1,p2,h1,h2,s1,s2, p1_,p2_,h1_,h2_,s1_,s2_) & - !$OMP private(ok, phase_Ji, phase_Ik, phase_Ii, contrib2, contrib, delta_IJk, HJk, HIk, HIl, exc_Ik, exc_Ji, exc_Ii) & + !$OMP private(ok, phase_Ji, phase_Ik, phase_Ii, contrib, delta_IJk, HJk, HIk, HIl, exc_Ik, exc_Ji, exc_Ii) & !$OMP private(det_tmp, det_tmp2, II, blok) & !$OMP shared(idx_sorted_bit, N_det_non_ref, N_det_ref, N_int, psi_non_ref, psi_non_ref_coef, psi_ref, psi_ref_coef) & - !$OMP shared(i_state,lambda_mrcc, hf_bitmask, active_sorb) + !$OMP shared(i_state,lambda_mrcc, hf_bitmask, active_sorb,c0) do i=1,N_det_non_ref if(mod(i,1000) == 0) print *, i, "/", N_det_non_ref do J=1,N_det_ref @@ -1553,15 +1410,8 @@ END_PROVIDER call apply_excitation(psi_non_ref(1,1,i),exc_Ik,det_tmp,ok,N_int) if(ok) cycle contrib = delta_IJk * HIl * lambda_mrcc(i_state,l) - if(dabs(psi_ref_coef(II,i_state)).ge.1.d-3) then - contrib2 = contrib / psi_ref_coef(II, i_state) * psi_non_ref_coef(l,i_state) - !$OMP ATOMIC - delta_sub_ii(II,i_state) -= contrib2 - else - contrib = contrib * 0.5d0 - endif !$OMP ATOMIC - delta_sub_ij(II, i, i_state) += contrib + delta_sub_ij(i, i_state) += contrib* c0(i_state) * psi_ref_coef(II,i_state) end do end do end do diff --git a/plugins/mrcepa0/dressing_slave.irp.f b/plugins/mrcepa0/dressing_slave.irp.f index fa486101..b0c3a360 100644 --- a/plugins/mrcepa0/dressing_slave.irp.f +++ b/plugins/mrcepa0/dressing_slave.irp.f @@ -402,17 +402,15 @@ end -subroutine mrsc2_dressing_collector(zmq_socket_pull,delta_ii_,delta_ij_,delta_ii_s2_,delta_ij_s2_) +subroutine mrsc2_dressing_collector(zmq_socket_pull,delta_ij_,delta_ij_s2_) use f77_zmq implicit none BEGIN_DOC ! Collects results from the AO integral calculation END_DOC - double precision,intent(inout) :: delta_ij_(N_states,N_det_non_ref,N_det_ref) - double precision,intent(inout) :: delta_ii_(N_states,N_det_ref) - double precision,intent(inout) :: delta_ij_s2_(N_states,N_det_non_ref,N_det_ref) - double precision,intent(inout) :: delta_ii_s2_(N_states,N_det_ref) + double precision,intent(inout) :: delta_ij_(N_states,N_det_non_ref) + double precision,intent(inout) :: delta_ij_s2_(N_states,N_det_non_ref) integer(ZMQ_PTR), intent(in) :: zmq_socket_pull ! integer :: j,l @@ -431,15 +429,18 @@ subroutine mrsc2_dressing_collector(zmq_socket_pull,delta_ii_,delta_ij_,delta_ii integer :: I_i, J, l, i_state, n(2), kk integer,allocatable :: idx(:,:) - delta_ii_(:,:) = 0d0 - delta_ij_(:,:,:) = 0d0 - delta_ii_s2_(:,:) = 0d0 - delta_ij_s2_(:,:,:) = 0d0 + delta_ij_(:,:) = 0d0 + delta_ij_s2_(:,:) = 0d0 zmq_to_qp_run_socket = new_zmq_to_qp_run_socket() allocate ( delta(N_states,0:N_det_non_ref,2), delta_s2(N_states,0:N_det_non_ref,2) ) + double precision :: c0(N_states) + do i_state=1,N_states + c0(i_state) = 1.d0/psi_coef(dressed_column_idx(i_state),i_state) + enddo + allocate(idx(N_det_non_ref,2)) more = 1 do while (more == 1) @@ -449,34 +450,19 @@ subroutine mrsc2_dressing_collector(zmq_socket_pull,delta_ii_,delta_ij_,delta_ii do l=1, n(1) do i_state=1,N_states - delta_ij_(i_state,idx(l,1),i_I) += delta(i_state,l,1) - delta_ij_s2_(i_state,idx(l,1),i_I) += delta_s2(i_state,l,1) + delta_ij_(i_state,idx(l,1)) += delta(i_state,l,1) * psi_ref_coef(i_I,i_state) * c0(i_state) + delta_ij_s2_(i_state,idx(l,1)) += delta_s2(i_state,l,1) * psi_ref_coef(i_I,i_state) * c0(i_state) end do end do do l=1, n(2) do i_state=1,N_states - delta_ij_(i_state,idx(l,2),J) += delta(i_state,l,2) - delta_ij_s2_(i_state,idx(l,2),J) += delta_s2(i_state,l,2) + delta_ij_(i_state,idx(l,2)) += delta(i_state,l,2) * psi_ref_coef(J,i_state) * c0(i_state) + delta_ij_s2_(i_state,idx(l,2)) += delta_s2(i_state,l,2) * psi_ref_coef(J,i_state) * c0(i_state) end do end do - if(n(1) /= 0) then - do i_state=1,N_states - delta_ii_(i_state,i_I) += delta(i_state,0,1) - delta_ii_s2_(i_state,i_I) += delta_s2(i_state,0,1) - end do - end if - - if(n(2) /= 0) then - do i_state=1,N_states - delta_ii_(i_state,J) += delta(i_state,0,2) - delta_ii_s2_(i_state,J) += delta_s2(i_state,0,2) - end do - end if - - if (task_id /= 0) then integer, external :: zmq_delete_task if (zmq_delete_task(zmq_to_qp_run_socket,zmq_socket_pull,task_id,more) == -1) then @@ -495,10 +481,8 @@ end - BEGIN_PROVIDER [ double precision, delta_ij_old, (N_states,N_det_non_ref,N_det_ref) ] -&BEGIN_PROVIDER [ double precision, delta_ii_old, (N_states,N_det_ref) ] -&BEGIN_PROVIDER [ double precision, delta_ij_s2_old, (N_states,N_det_non_ref,N_det_ref) ] -&BEGIN_PROVIDER [ double precision, delta_ii_s2_old, (N_states,N_det_ref) ] + BEGIN_PROVIDER [ double precision, delta_ij_old, (N_states,N_det_non_ref) ] +&BEGIN_PROVIDER [ double precision, delta_ij_s2_old, (N_states,N_det_non_ref) ] implicit none integer :: i_state, i, i_I, J, k, kk, degree, degree2, m, l, deg, ni, m2 @@ -612,11 +596,11 @@ end print *, nzer, ntot, float(nzer) / float(ntot) provide nproc !$OMP PARALLEL DEFAULT(none) & - !$OMP SHARED(delta_ii_old,delta_ij_old,delta_ii_s2_old,delta_ij_s2_old,zmq_socket_pull)& + !$OMP SHARED(delta_ij_old,delta_ij_s2_old,zmq_socket_pull)& !$OMP PRIVATE(i) NUM_THREADS(nproc+1) i = omp_get_thread_num() if (i==0) then - call mrsc2_dressing_collector(zmq_socket_pull,delta_ii_old,delta_ij_old,delta_ii_s2_old,delta_ij_s2_old) + call mrsc2_dressing_collector(zmq_socket_pull,delta_ij_old,delta_ij_s2_old) else call mrsc2_dressing_slave_inproc(i) endif diff --git a/plugins/mrcepa0/mrcepa0_general.irp.f b/plugins/mrcepa0/mrcepa0_general.irp.f index f58008a0..3c390763 100644 --- a/plugins/mrcepa0/mrcepa0_general.irp.f +++ b/plugins/mrcepa0/mrcepa0_general.irp.f @@ -14,8 +14,6 @@ subroutine run(N_st,energy) integer :: n_it_mrcc_max double precision :: thresh_mrcc - double precision, allocatable :: lambda(:) - allocate (lambda(N_states)) thresh_mrcc = thresh_dressed_ci n_it_mrcc_max = n_it_max_dressed_ci @@ -34,7 +32,6 @@ subroutine run(N_st,energy) E_new = 0.d0 delta_E = 1.d0 iteration = 0 - lambda = 1.d0 do while (delta_E > thresh_mrcc) iteration += 1 print *, '===============================================' @@ -45,12 +42,9 @@ subroutine run(N_st,energy) do i=1,N_st call write_double(6,ci_energy_dressed(i),"Energy") enddo - call diagonalize_ci_dressed(lambda) + call diagonalize_ci_dressed E_new = mrcc_e0_denominator(1) !sum(ci_energy_dressed(1:N_states)) -! if (.true.) then -! provide delta_ij_mrcc_pouet -! endif delta_E = (E_new - E_old)/dble(N_states) print *, '' call write_double(6,thresh_mrcc,"thresh_mrcc") diff --git a/plugins/mrcepa0/run_mrcc_slave.irp.f b/plugins/mrcepa0/run_mrcc_slave.irp.f index 3b3cfe44..c2d871e0 100644 --- a/plugins/mrcepa0/run_mrcc_slave.irp.f +++ b/plugins/mrcepa0/run_mrcc_slave.irp.f @@ -35,17 +35,13 @@ subroutine run_mrcc_slave(thread,iproc,energy) integer(bit_kind) :: mask(N_int,2), omask(N_int,2) double precision,allocatable :: delta_ij_loc(:,:,:) - double precision,allocatable :: delta_ii_loc(:,:) !double precision,allocatable :: delta_ij_s2_loc(:,:,:) - !double precision,allocatable :: delta_ii_s2_loc(:,:) integer :: h,p,n logical :: ok double precision :: contrib(N_states) - allocate(delta_ij_loc(N_states,N_det_non_ref,2) & - ,delta_ii_loc(N_states,2))! & + allocate(delta_ij_loc(N_states,N_det_non_ref,2) ) !,delta_ij_s2_loc(N_states,N_det_non_ref,N_det_ref) & - !,delta_ii_s2_loc(N_states, N_det_ref)) allocate(abuf(N_int, 2, N_det_non_ref)) @@ -82,9 +78,7 @@ subroutine run_mrcc_slave(thread,iproc,energy) contrib = 0d0 i_generator = ind(i_i_generator) delta_ij_loc = 0d0 - delta_ii_loc = 0d0 !delta_ij_s2_loc = 0d0 - !delta_ii_s2_loc = 0d0 !call select_connected(i_generator,energy,mrcc_detail(1, i_i_generator),buf,subset) !!!!!!!!!!!!!!!!!!!!!! @@ -103,7 +97,7 @@ subroutine run_mrcc_slave(thread,iproc,energy) n = n - 1 if(n /= 0) then - call mrcc_part_dress_1c(delta_ij_loc(1,1,1), delta_ii_loc(1,1), delta_ij_loc(1,1,2), delta_ii_loc(1,2), & + call mrcc_part_dress_1c(delta_ij_loc(1,1,1), delta_ij_loc(1,1,2), & i_generator,n,abuf,N_int,omask,contrib) endif end do diff --git a/src/Davidson/NEEDED_CHILDREN_MODULES b/src/Davidson/NEEDED_CHILDREN_MODULES index aae89501..9361eccd 100644 --- a/src/Davidson/NEEDED_CHILDREN_MODULES +++ b/src/Davidson/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -Determinants +Determinants DavidsonDressed diff --git a/src/Davidson/diagonalize_CI.irp.f b/src/Davidson/diagonalize_CI.irp.f index c65d9763..39edad18 100644 --- a/src/Davidson/diagonalize_CI.irp.f +++ b/src/Davidson/diagonalize_CI.irp.f @@ -65,7 +65,7 @@ END_PROVIDER call davidson_diag_HS2(psi_det,CI_eigenvectors, CI_eigenvectors_s2, & size(CI_eigenvectors,1),CI_electronic_energy, & - N_det,min(N_det,N_states),min(N_det,N_states_diag),N_int,6) + N_det,min(N_det,N_states),min(N_det,N_states_diag),N_int,0) else if (diag_algorithm == "Lapack") then From e44c040434eeaf6d51d18074f64ea194bf0282bf Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 8 Feb 2018 22:16:58 +0100 Subject: [PATCH 12/14] Fixed travis --- .travis.yml | 1 + install/scripts/install_lapack.sh | 6 +++++- {plugins => src}/DavidsonDressed/README.rst | 0 .../DavidsonDressed/diagonalization_hs2_dressed.irp.f | 0 {plugins => src}/DavidsonUndressed/NEEDED_CHILDREN_MODULES | 0 {plugins => src}/DavidsonUndressed/README.rst | 0 {plugins => src}/DavidsonUndressed/davidson_slave.irp.f | 0 .../DavidsonUndressed/diag_restart_save_all_states.irp.f | 0 .../DavidsonUndressed/diag_restart_save_lowest_state.irp.f | 0 .../DavidsonUndressed/diag_restart_save_one_state.irp.f | 0 {plugins => src}/DavidsonUndressed/guess_lowest_state.irp.f | 0 .../DavidsonUndressed/print_H_matrix_restart.irp.f | 0 {plugins => src}/DavidsonUndressed/print_energy.irp.f | 0 13 files changed, 6 insertions(+), 1 deletion(-) rename {plugins => src}/DavidsonDressed/README.rst (100%) rename {plugins => src}/DavidsonDressed/diagonalization_hs2_dressed.irp.f (100%) rename {plugins => src}/DavidsonUndressed/NEEDED_CHILDREN_MODULES (100%) rename {plugins => src}/DavidsonUndressed/README.rst (100%) rename {plugins => src}/DavidsonUndressed/davidson_slave.irp.f (100%) rename {plugins => src}/DavidsonUndressed/diag_restart_save_all_states.irp.f (100%) rename {plugins => src}/DavidsonUndressed/diag_restart_save_lowest_state.irp.f (100%) rename {plugins => src}/DavidsonUndressed/diag_restart_save_one_state.irp.f (100%) rename {plugins => src}/DavidsonUndressed/guess_lowest_state.irp.f (100%) rename {plugins => src}/DavidsonUndressed/print_H_matrix_restart.irp.f (100%) rename {plugins => src}/DavidsonUndressed/print_energy.irp.f (100%) diff --git a/.travis.yml b/.travis.yml index dd28c132..10b19f7f 100644 --- a/.travis.yml +++ b/.travis.yml @@ -24,6 +24,7 @@ addons: cache: directories: - $HOME/.opam/ + - $HOME/lapack-release language: python python: diff --git a/install/scripts/install_lapack.sh b/install/scripts/install_lapack.sh index 25cbefc6..e23a14d8 100755 --- a/install/scripts/install_lapack.sh +++ b/install/scripts/install_lapack.sh @@ -1,6 +1,10 @@ #!/bin/bash -x -git clone https://github.com/Reference-LAPACK/lapack-release.git +if [[ ! -d lapack-release ]] +then + git clone https://github.com/Reference-LAPACK/lapack-release.git +fi + cd lapack-release cp make.inc.example make.inc make -j 8 diff --git a/plugins/DavidsonDressed/README.rst b/src/DavidsonDressed/README.rst similarity index 100% rename from plugins/DavidsonDressed/README.rst rename to src/DavidsonDressed/README.rst diff --git a/plugins/DavidsonDressed/diagonalization_hs2_dressed.irp.f b/src/DavidsonDressed/diagonalization_hs2_dressed.irp.f similarity index 100% rename from plugins/DavidsonDressed/diagonalization_hs2_dressed.irp.f rename to src/DavidsonDressed/diagonalization_hs2_dressed.irp.f diff --git a/plugins/DavidsonUndressed/NEEDED_CHILDREN_MODULES b/src/DavidsonUndressed/NEEDED_CHILDREN_MODULES similarity index 100% rename from plugins/DavidsonUndressed/NEEDED_CHILDREN_MODULES rename to src/DavidsonUndressed/NEEDED_CHILDREN_MODULES diff --git a/plugins/DavidsonUndressed/README.rst b/src/DavidsonUndressed/README.rst similarity index 100% rename from plugins/DavidsonUndressed/README.rst rename to src/DavidsonUndressed/README.rst diff --git a/plugins/DavidsonUndressed/davidson_slave.irp.f b/src/DavidsonUndressed/davidson_slave.irp.f similarity index 100% rename from plugins/DavidsonUndressed/davidson_slave.irp.f rename to src/DavidsonUndressed/davidson_slave.irp.f diff --git a/plugins/DavidsonUndressed/diag_restart_save_all_states.irp.f b/src/DavidsonUndressed/diag_restart_save_all_states.irp.f similarity index 100% rename from plugins/DavidsonUndressed/diag_restart_save_all_states.irp.f rename to src/DavidsonUndressed/diag_restart_save_all_states.irp.f diff --git a/plugins/DavidsonUndressed/diag_restart_save_lowest_state.irp.f b/src/DavidsonUndressed/diag_restart_save_lowest_state.irp.f similarity index 100% rename from plugins/DavidsonUndressed/diag_restart_save_lowest_state.irp.f rename to src/DavidsonUndressed/diag_restart_save_lowest_state.irp.f diff --git a/plugins/DavidsonUndressed/diag_restart_save_one_state.irp.f b/src/DavidsonUndressed/diag_restart_save_one_state.irp.f similarity index 100% rename from plugins/DavidsonUndressed/diag_restart_save_one_state.irp.f rename to src/DavidsonUndressed/diag_restart_save_one_state.irp.f diff --git a/plugins/DavidsonUndressed/guess_lowest_state.irp.f b/src/DavidsonUndressed/guess_lowest_state.irp.f similarity index 100% rename from plugins/DavidsonUndressed/guess_lowest_state.irp.f rename to src/DavidsonUndressed/guess_lowest_state.irp.f diff --git a/plugins/DavidsonUndressed/print_H_matrix_restart.irp.f b/src/DavidsonUndressed/print_H_matrix_restart.irp.f similarity index 100% rename from plugins/DavidsonUndressed/print_H_matrix_restart.irp.f rename to src/DavidsonUndressed/print_H_matrix_restart.irp.f diff --git a/plugins/DavidsonUndressed/print_energy.irp.f b/src/DavidsonUndressed/print_energy.irp.f similarity index 100% rename from plugins/DavidsonUndressed/print_energy.irp.f rename to src/DavidsonUndressed/print_energy.irp.f From 8e689236be5dada408876da89953dbede6a89894 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 8 Feb 2018 22:22:21 +0100 Subject: [PATCH 13/14] Forgot files --- src/DavidsonDressed/NEEDED_CHILDREN_MODULES | 1 + src/DavidsonDressed/diagonalize_CI.irp.f | 210 ++++++++++++++++++++ 2 files changed, 211 insertions(+) create mode 100644 src/DavidsonDressed/NEEDED_CHILDREN_MODULES create mode 100644 src/DavidsonDressed/diagonalize_CI.irp.f diff --git a/src/DavidsonDressed/NEEDED_CHILDREN_MODULES b/src/DavidsonDressed/NEEDED_CHILDREN_MODULES new file mode 100644 index 00000000..8b137891 --- /dev/null +++ b/src/DavidsonDressed/NEEDED_CHILDREN_MODULES @@ -0,0 +1 @@ + diff --git a/src/DavidsonDressed/diagonalize_CI.irp.f b/src/DavidsonDressed/diagonalize_CI.irp.f new file mode 100644 index 00000000..7b12bc1c --- /dev/null +++ b/src/DavidsonDressed/diagonalize_CI.irp.f @@ -0,0 +1,210 @@ +BEGIN_PROVIDER [ double precision, CI_energy_dressed, (N_states_diag) ] + implicit none + BEGIN_DOC + ! N_states lowest eigenvalues of the CI matrix + END_DOC + + integer :: j + character*(8) :: st + call write_time(6) + do j=1,min(N_det,N_states_diag) + CI_energy_dressed(j) = CI_electronic_energy_dressed(j) + nuclear_repulsion + enddo + do j=1,min(N_det,N_states) + write(st,'(I4)') j + call write_double(6,CI_energy_dressed(j),'Energy of state '//trim(st)) + call write_double(6,CI_eigenvectors_s2_dressed(j),'S^2 of state '//trim(st)) + enddo + +END_PROVIDER + + BEGIN_PROVIDER [ double precision, CI_electronic_energy_dressed, (N_states_diag) ] +&BEGIN_PROVIDER [ double precision, CI_eigenvectors_dressed, (N_det,N_states_diag) ] +&BEGIN_PROVIDER [ double precision, CI_eigenvectors_s2_dressed, (N_states_diag) ] + BEGIN_DOC + ! Eigenvectors/values of the CI matrix + END_DOC + implicit none + double precision :: ovrlp,u_dot_v + integer :: i_good_state + integer, allocatable :: index_good_state_array(:) + logical, allocatable :: good_state_array(:) + double precision, allocatable :: s2_values_tmp(:) + integer :: i_other_state + double precision, allocatable :: eigenvectors(:,:), eigenvectors_s2(:,:), eigenvalues(:) + integer :: i_state + double precision :: e_0 + integer :: i,j,k,mrcc_state + double precision, allocatable :: s2_eigvalues(:) + double precision, allocatable :: e_array(:) + integer, allocatable :: iorder(:) + + PROVIDE threshold_davidson nthreads_davidson + ! Guess values for the "N_states" states of the CI_eigenvectors_dressed + do j=1,min(N_states,N_det) + do i=1,N_det + CI_eigenvectors_dressed(i,j) = psi_coef(i,j) + enddo + enddo + + do j=min(N_states,N_det)+1,N_states_diag + do i=1,N_det + CI_eigenvectors_dressed(i,j) = 0.d0 + enddo + enddo + + if (diag_algorithm == "Davidson") then + + allocate (eigenvectors(size(CI_eigenvectors_dressed,1),size(CI_eigenvectors_dressed,2)),& + eigenvectors_s2(size(CI_eigenvectors_dressed,1),size(CI_eigenvectors_dressed,2)),& + eigenvalues(size(CI_electronic_energy_dressed,1))) + do j=1,min(N_states,N_det) + do i=1,N_det + eigenvectors(i,j) = psi_coef(i,j) + enddo + enddo + do mrcc_state=1,N_states + do j=mrcc_state,min(N_states,N_det) + do i=1,N_det + eigenvectors(i,j) = psi_coef(i,j) + enddo + enddo + call davidson_diag_HS2(psi_det,eigenvectors, eigenvectors_s2, & + size(eigenvectors,1), & + eigenvalues,N_det,min(N_det,N_states),min(N_det,N_states_diag),N_int,& + mrcc_state) + CI_eigenvectors_dressed(1:N_det,mrcc_state) = eigenvectors(1:N_det,mrcc_state) + CI_electronic_energy_dressed(mrcc_state) = eigenvalues(mrcc_state) + enddo + do k=N_states+1,N_states_diag + CI_eigenvectors_dressed(1:N_det,k) = eigenvectors(1:N_det,k) + CI_electronic_energy_dressed(k) = eigenvalues(k) + enddo + call u_0_S2_u_0(CI_eigenvectors_s2_dressed,CI_eigenvectors_dressed,N_det,psi_det,N_int,& + N_states_diag,size(CI_eigenvectors_dressed,1)) + + deallocate (eigenvectors,eigenvalues) + + + else if (diag_algorithm == "Lapack") then + + allocate (eigenvectors(size(H_matrix_dressed,1),N_det)) + allocate (eigenvalues(N_det)) + call lapack_diag(eigenvalues,eigenvectors, & + H_matrix_dressed,size(H_matrix_dressed,1),N_det) + CI_electronic_energy_dressed(:) = 0.d0 + if (s2_eig) then + i_state = 0 + allocate (s2_eigvalues(N_det)) + allocate(index_good_state_array(N_det),good_state_array(N_det)) + good_state_array = .False. + call u_0_S2_u_0(s2_eigvalues,eigenvectors,N_det,psi_det,N_int, & + N_det,size(eigenvectors,1)) + do j=1,N_det + ! Select at least n_states states with S^2 values closed to "expected_s2" + if(dabs(s2_eigvalues(j)-expected_s2).le.0.5d0)then + i_state +=1 + index_good_state_array(i_state) = j + good_state_array(j) = .True. + endif + if(i_state.eq.N_states) then + exit + endif + enddo + if(i_state .ne.0)then + ! Fill the first "i_state" states that have a correct S^2 value + do j = 1, i_state + do i=1,N_det + CI_eigenvectors_dressed(i,j) = eigenvectors(i,index_good_state_array(j)) + enddo + CI_electronic_energy_dressed(j) = eigenvalues(index_good_state_array(j)) + CI_eigenvectors_s2_dressed(j) = s2_eigvalues(index_good_state_array(j)) + enddo + i_other_state = 0 + do j = 1, N_det + if(good_state_array(j))cycle + i_other_state +=1 + if(i_state+i_other_state.gt.n_states_diag)then + exit + endif + do i=1,N_det + CI_eigenvectors_dressed(i,i_state+i_other_state) = eigenvectors(i,j) + enddo + CI_electronic_energy_dressed(i_state+i_other_state) = eigenvalues(j) + CI_eigenvectors_s2_dressed(i_state+i_other_state) = s2_eigvalues(i_state+i_other_state) + enddo + + else + print*,'' + print*,'!!!!!!!! WARNING !!!!!!!!!' + print*,' Within the ',N_det,'determinants selected' + print*,' and the ',N_states_diag,'states requested' + print*,' We did not find any state with S^2 values close to ',expected_s2 + print*,' We will then set the first N_states eigenvectors of the H matrix' + print*,' as the CI_eigenvectors_dressed' + print*,' You should consider more states and maybe ask for s2_eig to be .True. or just enlarge the CI space' + print*,'' + do j=1,min(N_states_diag,N_det) + do i=1,N_det + CI_eigenvectors_dressed(i,j) = eigenvectors(i,j) + enddo + CI_electronic_energy_dressed(j) = eigenvalues(j) + CI_eigenvectors_s2_dressed(j) = s2_eigvalues(j) + enddo + endif + deallocate(index_good_state_array,good_state_array) + deallocate(s2_eigvalues) + else + call u_0_S2_u_0(CI_eigenvectors_s2_dressed,eigenvectors,N_det,psi_det,N_int,& + min(N_det,N_states_diag),size(eigenvectors,1)) + ! Select the "N_states_diag" states of lowest energy + do j=1,min(N_det,N_states_diag) + do i=1,N_det + CI_eigenvectors_dressed(i,j) = eigenvectors(i,j) + enddo + CI_electronic_energy_dressed(j) = eigenvalues(j) + enddo + endif + deallocate(eigenvectors,eigenvalues) + endif + +END_PROVIDER + +subroutine diagonalize_CI_dressed + implicit none + BEGIN_DOC +! Replace the coefficients of the CI states by the coefficients of the +! eigenstates of the CI matrix + END_DOC + integer :: i,j + do j=1,N_states + do i=1,N_det + psi_coef(i,j) = CI_eigenvectors_dressed(i,j) + enddo + enddo + SOFT_TOUCH psi_coef +end + + + +BEGIN_PROVIDER [ double precision, h_matrix_dressed, (N_det,N_det,N_states) ] + implicit none + BEGIN_DOC + ! Dressed H with Delta_ij + END_DOC + integer :: i, j,istate,ii,jj + do istate = 1,N_states + do j=1,N_det + do i=1,N_det + h_matrix_dressed(i,j,istate) = h_matrix_all_dets(i,j) + enddo + enddo + i = dressed_column_idx(istate) + do j = 1, N_det + h_matrix_dressed(i,j,istate) += dressing_column_h(j,istate) + h_matrix_dressed(j,i,istate) += dressing_column_h(j,istate) + enddo + h_matrix_dressed(i,i,istate) -= dressing_column_h(i,istate) + enddo +END_PROVIDER + From 08c2a46f449e958dd0d2f19ac3900995edf7d2cf Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 8 Feb 2018 23:05:17 +0100 Subject: [PATCH 14/14] Fix travis --- install/scripts/install_lapack.sh | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/install/scripts/install_lapack.sh b/install/scripts/install_lapack.sh index e23a14d8..9ea47619 100755 --- a/install/scripts/install_lapack.sh +++ b/install/scripts/install_lapack.sh @@ -1,10 +1,6 @@ #!/bin/bash -x -if [[ ! -d lapack-release ]] -then - git clone https://github.com/Reference-LAPACK/lapack-release.git -fi - +git clone https://github.com/Reference-LAPACK/lapack-release.git || echo "Clone failed" cd lapack-release cp make.inc.example make.inc make -j 8