diff --git a/plugins/Generators_CAS/Generators_full/.gitignore b/plugins/Generators_CAS/Generators_full/.gitignore new file mode 100644 index 00000000..8d85dede --- /dev/null +++ b/plugins/Generators_CAS/Generators_full/.gitignore @@ -0,0 +1,25 @@ +# Automatically created by /home/razoa/quantum_package/scripts/module/module_handler.py +IRPF90_temp +IRPF90_man +irpf90_entities +tags +irpf90.make +Makefile +Makefile.depend +build.ninja +.ninja_log +.ninja_deps +ezfio_interface.irp.f +Ezfio_files +Determinants +Integrals_Monoelec +MO_Basis +Utils +Pseudo +Bitmask +AO_Basis +Electrons +MOGuess +Nuclei +Hartree_Fock +Integrals_Bielec \ No newline at end of file diff --git a/plugins/Generators_CAS/Generators_full/NEEDED_CHILDREN_MODULES b/plugins/Generators_CAS/Generators_full/NEEDED_CHILDREN_MODULES new file mode 100644 index 00000000..54f54203 --- /dev/null +++ b/plugins/Generators_CAS/Generators_full/NEEDED_CHILDREN_MODULES @@ -0,0 +1 @@ +Determinants Hartree_Fock diff --git a/plugins/Generators_CAS/Generators_full/README.rst b/plugins/Generators_CAS/Generators_full/README.rst new file mode 100644 index 00000000..c30193a2 --- /dev/null +++ b/plugins/Generators_CAS/Generators_full/README.rst @@ -0,0 +1,61 @@ +====================== +Generators_full Module +====================== + +All the determinants of the wave function are generators. In this way, the Full CI +space is explored. + +Needed Modules +============== + +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. + +.. image:: tree_dependency.png + +* `Determinants `_ +* `Hartree_Fock `_ + +Needed Modules +============== +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. + + +.. image:: tree_dependency.png + +* `Determinants `_ +* `Hartree_Fock `_ + +Documentation +============= +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. + + +`degree_max_generators `_ + Max degree of excitation (respect to HF) of the generators + + +`n_det_generators `_ + For Single reference wave functions, the number of generators is 1 : the + Hartree-Fock determinant + + +`psi_coef_generators `_ + For Single reference wave functions, the generator is the + Hartree-Fock determinant + + +`psi_det_generators `_ + For Single reference wave functions, the generator is the + Hartree-Fock determinant + + +`select_max `_ + Memo to skip useless selectors + + +`size_select_max `_ + Size of the select_max array + diff --git a/plugins/Generators_CAS/Generators_full/generators.irp.f b/plugins/Generators_CAS/Generators_full/generators.irp.f new file mode 100644 index 00000000..eea5821b --- /dev/null +++ b/plugins/Generators_CAS/Generators_full/generators.irp.f @@ -0,0 +1,75 @@ +use bitmasks + +BEGIN_PROVIDER [ integer, N_det_generators ] + implicit none + BEGIN_DOC + ! For Single reference wave functions, the number of generators is 1 : the + ! Hartree-Fock determinant + END_DOC + integer :: i + double precision :: norm + call write_time(output_determinants) + norm = 0.d0 + N_det_generators = N_det + do i=1,N_det + norm = norm + psi_average_norm_contrib_sorted(i) + if (norm >= threshold_generators) then + N_det_generators = i + exit + endif + enddo + N_det_generators = max(N_det_generators,1) + call write_int(output_determinants,N_det_generators,'Number of generators') +END_PROVIDER + + BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators, (N_int,2,psi_det_size) ] +&BEGIN_PROVIDER [ double precision, psi_coef_generators, (psi_det_size,N_states) ] + implicit none + BEGIN_DOC + ! For Single reference wave functions, the generator is the + ! 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 + +END_PROVIDER + +BEGIN_PROVIDER [integer, degree_max_generators] + implicit none + BEGIN_DOC +! Max degree of excitation (respect to HF) of the generators + END_DOC + integer :: i,degree + degree_max_generators = 0 + do i = 1, N_det_generators + call get_excitation_degree(HF_bitmask,psi_det_generators(1,1,i),degree,N_int) + if(degree .gt. degree_max_generators)then + degree_max_generators = degree + endif + enddo +END_PROVIDER + +BEGIN_PROVIDER [ integer, size_select_max] + implicit none + BEGIN_DOC + ! Size of the select_max array + END_DOC + size_select_max = 10000 +END_PROVIDER + +BEGIN_PROVIDER [ double precision, select_max, (size_select_max) ] + implicit none + BEGIN_DOC + ! Memo to skip useless selectors + END_DOC + select_max = huge(1.d0) +END_PROVIDER + diff --git a/plugins/Generators_CAS/Generators_full/tree_dependency.png b/plugins/Generators_CAS/Generators_full/tree_dependency.png new file mode 100644 index 00000000..eed76866 Binary files /dev/null and b/plugins/Generators_CAS/Generators_full/tree_dependency.png differ diff --git a/plugins/Generators_CAS/generators.irp.f b/plugins/Generators_CAS/generators.irp.f index f47341de..10fbfaee 100644 --- a/plugins/Generators_CAS/generators.irp.f +++ b/plugins/Generators_CAS/generators.irp.f @@ -9,14 +9,14 @@ BEGIN_PROVIDER [ integer, N_det_generators ] logical :: good call write_time(output_determinants) N_det_generators = 0 - do i=1,N_det + do i=1,N_det_ref do l=1,n_cas_bitmask good = .True. do k=1,N_int good = good .and. ( & - iand(not(cas_bitmask(k,1,l)), psi_det(k,1,i)) == & + iand(not(cas_bitmask(k,1,l)), psi_ref(k,1,i)) == & iand(not(cas_bitmask(k,1,l)), HF_bitmask(k,1)) ) .and. ( & - iand(not(cas_bitmask(k,2,l)), psi_det(k,2,i)) == & + iand(not(cas_bitmask(k,2,l)), psi_ref(k,2,i)) == & iand(not(cas_bitmask(k,2,l)), HF_bitmask(k,2)) ) enddo if (good) then @@ -41,14 +41,14 @@ END_PROVIDER integer :: i, k, l, m logical :: good m=0 - do i=1,N_det + do i=1,N_det_ref do l=1,n_cas_bitmask good = .True. do k=1,N_int good = good .and. ( & - iand(not(cas_bitmask(k,1,l)), psi_det(k,1,i)) == & + iand(not(cas_bitmask(k,1,l)), psi_ref(k,1,i)) == & iand(not(cas_bitmask(k,1,l)), HF_bitmask(k,1)) .and. ( & - iand(not(cas_bitmask(k,2,l)), psi_det(k,2,i)) == & + iand(not(cas_bitmask(k,2,l)), psi_ref(k,2,i)) == & iand(not(cas_bitmask(k,2,l)), HF_bitmask(k,2) )) ) enddo if (good) then @@ -58,8 +58,8 @@ END_PROVIDER if (good) then m = m+1 do k=1,N_int - psi_det_generators(k,1,m) = psi_det(k,1,i) - psi_det_generators(k,2,m) = psi_det(k,2,i) + psi_det_generators(k,1,m) = psi_ref(k,1,i) + psi_det_generators(k,2,m) = psi_ref(k,2,i) enddo psi_coef_generators(m,:) = psi_coef(m,:) endif diff --git a/plugins/MRCC_Utils/NEEDED_CHILDREN_MODULES b/plugins/MRCC_Utils/NEEDED_CHILDREN_MODULES index 801d2f51..3dc21fd0 100644 --- a/plugins/MRCC_Utils/NEEDED_CHILDREN_MODULES +++ b/plugins/MRCC_Utils/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -Perturbation Selectors_full Generators_full Psiref_Utils Psiref_CAS +Perturbation Selectors_full Generators_full Psiref_Utils Psiref_CAS MRPT_Utils diff --git a/plugins/MRCC_Utils/mrcc_utils.irp.f b/plugins/MRCC_Utils/mrcc_utils.irp.f index 0540eed9..c1a277cf 100644 --- a/plugins/MRCC_Utils/mrcc_utils.irp.f +++ b/plugins/MRCC_Utils/mrcc_utils.irp.f @@ -629,6 +629,44 @@ END_PROVIDER call sort_det(psi_non_ref_sorted, psi_non_ref_sorted_idx, N_det_non_ref, N_int) END_PROVIDER + BEGIN_PROVIDER [ double precision, rho_mrpt, (N_det_non_ref, N_states) ] + implicit none + integer :: i, j, k + double precision :: coef_mrpt(N_States),coef_array(N_states),hij,delta_e(N_states) + double precision :: hij_array(N_det_Ref),delta_e_array(N_det_ref,N_states) + do i = 1, N_det_non_ref + do j = 1, N_det_ref + do k = 1, N_States + coef_array(j) = psi_ref_coef(j,k) + enddo + call i_h_j(psi_ref(1,1,j), psi_non_ref(1,1,i), N_int, Hij_array(j)) + call get_delta_e_dyall(psi_ref(1,1,j),psi_non_ref(1,1,i),coef_array,hij_array(j),delta_e) + print*,delta_e(:) + do k = 1, N_states + delta_e_Array(j,k) = delta_e(k) + enddo + enddo + do k = 1, N_states + do j = 1, N_det_Ref + coef_mrpt(k) += psi_ref_coef(j,k) * hij_array(j) / delta_e_array(j,k) + enddo + enddo + do k = 1, N_States + if(dabs(coef_mrpt(k)) .le.1.d-10)then + rho_mrpt(i,k) = 0.d0 + exit + endif + print*,k,psi_non_ref_coef(i,k) , coef_mrpt(k) + if(psi_non_ref_coef(i,k) / coef_mrpt(k) .lt.0d0)then + rho_mrpt(i,k) = 1.d0 + else + rho_mrpt(i,k) = psi_non_ref_coef(i,k) / coef_mrpt(k) + endif + enddo + enddo + + END_PROVIDER + BEGIN_PROVIDER [ double precision, dIj_unique, (hh_nex, N_states) ] &BEGIN_PROVIDER [ double precision, rho_mrcc, (N_det_non_ref, N_states) ] @@ -983,6 +1021,9 @@ double precision function get_dij_index(II, i, s, Nint) call get_phase(psi_ref(1,1,II), psi_non_ref(1,1,i), phase, N_int) get_dij_index = get_dij(psi_ref(1,1,II), psi_non_ref(1,1,i), s, Nint) * phase get_dij_index = get_dij_index + else if(lambda_type == 3) then + call i_h_j(psi_ref(1,1,II), psi_non_ref(1,1,i), Nint, HIi) + get_dij_index = HIi * rho_mrpt(i, s) end if end function diff --git a/plugins/MRPT/NEEDED_CHILDREN_MODULES b/plugins/MRPT/NEEDED_CHILDREN_MODULES index 7340c609..041b0136 100644 --- a/plugins/MRPT/NEEDED_CHILDREN_MODULES +++ b/plugins/MRPT/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -MRPT_Utils Selectors_full Generators_full +MRPT_Utils Selectors_full Psiref_CAS Generators_CAS diff --git a/src/MRPT_Utils/EZFIO.cfg b/plugins/MRPT_Utils/EZFIO.cfg similarity index 100% rename from src/MRPT_Utils/EZFIO.cfg rename to plugins/MRPT_Utils/EZFIO.cfg diff --git a/src/MRPT_Utils/H_apply.irp.f b/plugins/MRPT_Utils/H_apply.irp.f similarity index 100% rename from src/MRPT_Utils/H_apply.irp.f rename to plugins/MRPT_Utils/H_apply.irp.f diff --git a/src/MRPT_Utils/NEEDED_CHILDREN_MODULES b/plugins/MRPT_Utils/NEEDED_CHILDREN_MODULES similarity index 100% rename from src/MRPT_Utils/NEEDED_CHILDREN_MODULES rename to plugins/MRPT_Utils/NEEDED_CHILDREN_MODULES diff --git a/src/MRPT_Utils/README.rst b/plugins/MRPT_Utils/README.rst similarity index 100% rename from src/MRPT_Utils/README.rst rename to plugins/MRPT_Utils/README.rst diff --git a/src/MRPT_Utils/energies_cas.irp.f b/plugins/MRPT_Utils/energies_cas.irp.f similarity index 80% rename from src/MRPT_Utils/energies_cas.irp.f rename to plugins/MRPT_Utils/energies_cas.irp.f index 8f29717c..54e1a3f8 100644 --- a/src/MRPT_Utils/energies_cas.irp.f +++ b/plugins/MRPT_Utils/energies_cas.irp.f @@ -3,7 +3,7 @@ BEGIN_PROVIDER [ double precision, energy_cas_dyall, (N_states)] integer :: i double precision :: energies(N_states) do i = 1, N_states - call u0_H_dyall_u0(energies,psi_active,psi_coef,n_det,psi_det_size,psi_det_size,N_states,i) + call u0_H_dyall_u0(energies,psi_active,psi_ref_coef,n_det_ref,psi_det_size,psi_det_size,N_states,i) energy_cas_dyall(i) = energies(i) print*, 'energy_cas_dyall(i)', energy_cas_dyall(i) enddo @@ -15,7 +15,7 @@ BEGIN_PROVIDER [ double precision, energy_cas_dyall_no_exchange, (N_states)] integer :: i double precision :: energies(N_states) do i = 1, N_states - call u0_H_dyall_u0_no_exchange(energies,psi_active,psi_coef,n_det,psi_det_size,psi_det_size,N_states,i) + call u0_H_dyall_u0_no_exchange(energies,psi_active,psi_ref_coef,n_det_ref,psi_det_size,psi_det_size,N_states,i) energy_cas_dyall_no_exchange(i) = energies(i) print*, 'energy_cas_dyall(i)_no_exchange', energy_cas_dyall_no_exchange(i) enddo @@ -31,7 +31,7 @@ BEGIN_PROVIDER [ double precision, one_creat, (n_act_orb,2,N_states)] double precision :: norm_out(N_states) integer(bit_kind), allocatable :: psi_in_out(:,:,:) double precision, allocatable :: psi_in_out_coef(:,:) - allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states)) + allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det_ref,N_states)) use bitmasks integer :: iorb @@ -44,7 +44,7 @@ BEGIN_PROVIDER [ double precision, one_creat, (n_act_orb,2,N_states)] spin_exc = ispin do i = 1, n_det do j = 1, n_states - psi_in_out_coef(i,j) = psi_coef(i,j) + psi_in_out_coef(i,j) = psi_ref_coef(i,j) enddo do j = 1, N_int psi_in_out(j,1,i) = psi_active(j,1,i) @@ -53,8 +53,8 @@ BEGIN_PROVIDER [ double precision, one_creat, (n_act_orb,2,N_states)] enddo do state_target = 1,N_states call apply_exc_to_psi(orb,hole_particle,spin_exc, & - norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states) - call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target) + norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states) + call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target) one_creat(iorb,ispin,state_target) = energy_cas_dyall(state_target) - energies(state_target) enddo enddo @@ -72,7 +72,7 @@ BEGIN_PROVIDER [ double precision, one_anhil, (n_act_orb,2,N_states)] integer(bit_kind), allocatable :: psi_in_out(:,:,:) double precision, allocatable :: psi_in_out_coef(:,:) use bitmasks - allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states)) + allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det_ref,N_states)) integer :: iorb integer :: state_target @@ -84,7 +84,7 @@ BEGIN_PROVIDER [ double precision, one_anhil, (n_act_orb,2,N_states)] spin_exc = ispin do i = 1, n_det do j = 1, n_states - psi_in_out_coef(i,j) = psi_coef(i,j) + psi_in_out_coef(i,j) = psi_ref_coef(i,j) enddo do j = 1, N_int psi_in_out(j,1,i) = psi_active(j,1,i) @@ -93,8 +93,8 @@ BEGIN_PROVIDER [ double precision, one_anhil, (n_act_orb,2,N_states)] enddo do state_target = 1, N_states call apply_exc_to_psi(orb,hole_particle,spin_exc, & - norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states) - call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target) + norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states) + call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target) one_anhil(iorb,ispin,state_target) = energy_cas_dyall(state_target) - energies(state_target) enddo enddo @@ -113,7 +113,7 @@ BEGIN_PROVIDER [ double precision, two_creat, (n_act_orb,n_act_orb,2,2,N_states) integer(bit_kind), allocatable :: psi_in_out(:,:,:) double precision, allocatable :: psi_in_out_coef(:,:) use bitmasks - allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states)) + allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det_ref,N_states)) integer :: iorb,jorb integer :: state_target @@ -130,7 +130,7 @@ BEGIN_PROVIDER [ double precision, two_creat, (n_act_orb,n_act_orb,2,2,N_states) spin_exc_j = jspin do i = 1, n_det do j = 1, n_states - psi_in_out_coef(i,j) = psi_coef(i,j) + psi_in_out_coef(i,j) = psi_ref_coef(i,j) enddo do j = 1, N_int psi_in_out(j,1,i) = psi_active(j,1,i) @@ -139,10 +139,10 @@ BEGIN_PROVIDER [ double precision, two_creat, (n_act_orb,n_act_orb,2,2,N_states) enddo do state_target = 1 , N_states call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, & - norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states) + norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states) call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, & - norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states) - call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target) + norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states) + call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target) two_creat(iorb,jorb,ispin,jspin,state_target) = energy_cas_dyall(state_target) - energies(state_target) enddo enddo @@ -163,7 +163,7 @@ BEGIN_PROVIDER [ double precision, two_anhil, (n_act_orb,n_act_orb,2,2,N_states) integer(bit_kind), allocatable :: psi_in_out(:,:,:) double precision, allocatable :: psi_in_out_coef(:,:) use bitmasks - allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states)) + allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det_ref,N_states)) integer :: iorb,jorb integer :: state_target @@ -181,7 +181,7 @@ BEGIN_PROVIDER [ double precision, two_anhil, (n_act_orb,n_act_orb,2,2,N_states) spin_exc_j = jspin do i = 1, n_det do j = 1, n_states - psi_in_out_coef(i,j) = psi_coef(i,j) + psi_in_out_coef(i,j) = psi_ref_coef(i,j) enddo do j = 1, N_int psi_in_out(j,1,i) = psi_active(j,1,i) @@ -189,10 +189,10 @@ BEGIN_PROVIDER [ double precision, two_anhil, (n_act_orb,n_act_orb,2,2,N_states) enddo enddo call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, & - norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states) + norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states) call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, & - norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states) - call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target) + norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states) + call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target) two_anhil(iorb,jorb,ispin,jspin,state_target) = energy_cas_dyall(state_target) - energies(state_target) enddo enddo @@ -213,7 +213,7 @@ BEGIN_PROVIDER [ double precision, one_anhil_one_creat, (n_act_orb,n_act_orb,2,2 double precision, allocatable :: psi_in_out_coef(:,:) use bitmasks - allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states)) + allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det_ref,N_states)) integer :: iorb,jorb integer :: state_target double precision :: energies(n_states) @@ -229,7 +229,7 @@ BEGIN_PROVIDER [ double precision, one_anhil_one_creat, (n_act_orb,n_act_orb,2,2 spin_exc_j = jspin do i = 1, n_det do j = 1, n_states - psi_in_out_coef(i,j) = psi_coef(i,j) + psi_in_out_coef(i,j) = psi_ref_coef(i,j) enddo do j = 1, N_int psi_in_out(j,1,i) = psi_active(j,1,i) @@ -238,14 +238,14 @@ BEGIN_PROVIDER [ double precision, one_anhil_one_creat, (n_act_orb,n_act_orb,2,2 enddo do state_target = 1, N_states call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, & - norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states) + norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states) call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, & - norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states) + norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states) !if(orb_i == orb_j .and. ispin .ne. jspin)then - call u0_H_dyall_u0_no_exchange(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target) + call u0_H_dyall_u0_no_exchange(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target) one_anhil_one_creat(iorb,jorb,ispin,jspin,state_target) = energy_cas_dyall_no_exchange(state_target) - energies(state_target) !else - ! call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target) + ! call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target) ! one_anhil_one_creat(iorb,jorb,ispin,jspin,state_target) = energy_cas_dyall(state_target) - energies(state_target) !endif enddo @@ -268,7 +268,7 @@ BEGIN_PROVIDER [ double precision, two_anhil_one_creat, (n_act_orb,n_act_orb,n_a integer(bit_kind), allocatable :: psi_in_out(:,:,:) double precision, allocatable :: psi_in_out_coef(:,:) use bitmasks - allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states)) + allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det_ref,N_states)) integer :: iorb,jorb integer :: korb @@ -291,7 +291,7 @@ BEGIN_PROVIDER [ double precision, two_anhil_one_creat, (n_act_orb,n_act_orb,n_a spin_exc_k = kspin do i = 1, n_det do j = 1, n_states - psi_in_out_coef(i,j) = psi_coef(i,j) + psi_in_out_coef(i,j) = psi_ref_coef(i,j) enddo do j = 1, N_int psi_in_out(j,1,i) = psi_active(j,1,i) @@ -301,12 +301,12 @@ BEGIN_PROVIDER [ double precision, two_anhil_one_creat, (n_act_orb,n_act_orb,n_a do state_target = 1, N_states call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, & - norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states) + norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states) call apply_exc_to_psi(orb_k,hole_particle_k,spin_exc_k, & - norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states) + norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states) call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, & - norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states) - call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target) + norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states) + call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target) two_anhil_one_creat(iorb,jorb,korb,ispin,jspin,kspin,state_target) = energy_cas_dyall(state_target) - energies(state_target) enddo enddo @@ -330,7 +330,7 @@ BEGIN_PROVIDER [ double precision, two_creat_one_anhil, (n_act_orb,n_act_orb,n_a integer(bit_kind), allocatable :: psi_in_out(:,:,:) double precision, allocatable :: psi_in_out_coef(:,:) use bitmasks - allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states)) + allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det_ref,N_states)) integer :: iorb,jorb integer :: korb @@ -353,7 +353,7 @@ BEGIN_PROVIDER [ double precision, two_creat_one_anhil, (n_act_orb,n_act_orb,n_a spin_exc_k = kspin do i = 1, n_det do j = 1, n_states - psi_in_out_coef(i,j) = psi_coef(i,j) + psi_in_out_coef(i,j) = psi_ref_coef(i,j) enddo do j = 1, N_int psi_in_out(j,1,i) = psi_active(j,1,i) @@ -362,12 +362,12 @@ BEGIN_PROVIDER [ double precision, two_creat_one_anhil, (n_act_orb,n_act_orb,n_a enddo do state_target = 1, N_states call apply_exc_to_psi(orb_k,hole_particle_k,spin_exc_k, & - norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states) + norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states) call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, & - norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states) + norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states) call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, & - norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states) - call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target) + norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states) + call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target) two_creat_one_anhil(iorb,jorb,korb,ispin,jspin,kspin,state_target) = energy_cas_dyall(state_target) - energies(state_target) enddo enddo @@ -391,7 +391,7 @@ BEGIN_PROVIDER [ double precision, three_creat, (n_act_orb,n_act_orb,n_act_orb,2 integer(bit_kind), allocatable :: psi_in_out(:,:,:) double precision, allocatable :: psi_in_out_coef(:,:) use bitmasks - allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states)) + allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det_ref,N_states)) integer :: iorb,jorb integer :: korb @@ -414,7 +414,7 @@ BEGIN_PROVIDER [ double precision, three_creat, (n_act_orb,n_act_orb,n_act_orb,2 spin_exc_k = kspin do i = 1, n_det do j = 1, n_states - psi_in_out_coef(i,j) = psi_coef(i,j) + psi_in_out_coef(i,j) = psi_ref_coef(i,j) enddo do j = 1, N_int psi_in_out(j,1,i) = psi_active(j,1,i) @@ -423,12 +423,12 @@ BEGIN_PROVIDER [ double precision, three_creat, (n_act_orb,n_act_orb,n_act_orb,2 enddo do state_target = 1, N_states call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, & - norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states) + norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states) call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, & - norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states) + norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states) call apply_exc_to_psi(orb_k,hole_particle_k,spin_exc_k, & - norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states) - call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target) + norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states) + call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target) three_creat(iorb,jorb,korb,ispin,jspin,kspin,state_target) = energy_cas_dyall(state_target) - energies(state_target) enddo enddo @@ -452,7 +452,7 @@ BEGIN_PROVIDER [ double precision, three_anhil, (n_act_orb,n_act_orb,n_act_orb,2 integer(bit_kind), allocatable :: psi_in_out(:,:,:) double precision, allocatable :: psi_in_out_coef(:,:) use bitmasks - allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states)) + allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det_ref,N_states)) integer :: iorb,jorb integer :: korb @@ -475,7 +475,7 @@ BEGIN_PROVIDER [ double precision, three_anhil, (n_act_orb,n_act_orb,n_act_orb,2 spin_exc_k = kspin do i = 1, n_det do j = 1, n_states - psi_in_out_coef(i,j) = psi_coef(i,j) + psi_in_out_coef(i,j) = psi_ref_coef(i,j) enddo do j = 1, N_int psi_in_out(j,1,i) = psi_active(j,1,i) @@ -484,12 +484,12 @@ BEGIN_PROVIDER [ double precision, three_anhil, (n_act_orb,n_act_orb,n_act_orb,2 enddo do state_target = 1, N_states call apply_exc_to_psi(orb_i,hole_particle_i,spin_exc_i, & - norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states) + norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states) call apply_exc_to_psi(orb_j,hole_particle_j,spin_exc_j, & - norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states) + norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states) call apply_exc_to_psi(orb_k,hole_particle_k,spin_exc_k, & - norm_out,psi_in_out,psi_in_out_coef, n_det,n_det,n_det,N_states) - call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target) + norm_out,psi_in_out,psi_in_out_coef, n_det_ref,n_det_ref,n_det_ref,N_states) + call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target) three_anhil(iorb,jorb,korb,ispin,jspin,kspin,state_target) = energy_cas_dyall(state_target) - energies(state_target) enddo enddo @@ -515,7 +515,7 @@ END_PROVIDER integer(bit_kind), allocatable :: psi_in_out(:,:,:) double precision, allocatable :: psi_in_out_coef(:,:) use bitmasks - allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states)) + allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det_ref,N_states)) integer :: iorb,jorb,i_ok integer :: state_target @@ -543,8 +543,8 @@ END_PROVIDER enddo do i = 1, n_det do j = 1, N_int - psi_in_out(j,1,i) = psi_det(j,1,i) - psi_in_out(j,2,i) = psi_det(j,2,i) + psi_in_out(j,1,i) = psi_ref(j,1,i) + psi_in_out(j,2,i) = psi_ref(j,2,i) enddo call do_mono_excitation(psi_in_out(1,1,i),orb_i,orb_v,ispin,i_ok) if(i_ok.ne.1)then @@ -552,10 +552,10 @@ END_PROVIDER call debug_det(psi_in_out,N_int) print*, 'pb, i_ok ne 0 !!!' endif - call i_H_j(psi_in_out(1,1,i),psi_det(1,1,i),N_int,hij) + call i_H_j(psi_in_out(1,1,i),psi_ref(1,1,i),N_int,hij) do j = 1, n_states double precision :: coef,contrib - coef = psi_coef(i,j) !* psi_coef(i,j) + coef = psi_ref_coef(i,j) !* psi_ref_coef(i,j) psi_in_out_coef(i,j) = coef * hij norm(j,ispin) += psi_in_out_coef(i,j) * psi_in_out_coef(i,j) enddo @@ -571,7 +571,7 @@ END_PROVIDER norm(j,ispin) = 1.d0/dsqrt(norm(j,ispin)) endif enddo - do i = 1, N_det + do i = 1, N_det_ref do j = 1, N_states psi_in_out_coef(i,j) = psi_in_out_coef(i,j) * norm(j,ispin) norm_bis(j,ispin) += psi_in_out_coef(i,j) * psi_in_out_coef(i,j) @@ -584,8 +584,8 @@ END_PROVIDER do state_target = 1, N_states energies_alpha_beta(state_target, ispin) = 0.d0 if(norm(state_target,ispin) .ne. 0.d0 .and. dabs(norm_no_inv(state_target,ispin)) .gt. thresh_norm)then -! call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target) - call u0_H_dyall_u0_no_exchange(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target) +! call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target) + call u0_H_dyall_u0_no_exchange(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target) energies_alpha_beta(state_target, ispin) += energies(state_target) endif enddo @@ -616,7 +616,7 @@ BEGIN_PROVIDER [ double precision, one_anhil_inact, (n_inact_orb,n_act_orb,N_Sta integer(bit_kind), allocatable :: psi_in_out(:,:,:) double precision, allocatable :: psi_in_out_coef(:,:) use bitmasks - allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states)) + allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det_ref,N_states)) integer :: jorb,i_ok,aorb,orb_a integer :: state_target @@ -643,8 +643,8 @@ BEGIN_PROVIDER [ double precision, one_anhil_inact, (n_inact_orb,n_act_orb,N_Sta do ispin = 1,2 do i = 1, n_det do j = 1, N_int - psi_in_out(j,1,i) = psi_det(j,1,i) - psi_in_out(j,2,i) = psi_det(j,2,i) + psi_in_out(j,1,i) = psi_ref(j,1,i) + psi_in_out(j,2,i) = psi_ref(j,2,i) enddo call do_mono_excitation(psi_in_out(1,1,i),orb_i,orb_a,ispin,i_ok) if(i_ok.ne.1)then @@ -652,11 +652,11 @@ BEGIN_PROVIDER [ double precision, one_anhil_inact, (n_inact_orb,n_act_orb,N_Sta psi_in_out_coef(i,j) = 0.d0 enddo else - call i_H_j(psi_in_out(1,1,i),psi_det(1,1,i),N_int,hij) + call i_H_j(psi_in_out(1,1,i),psi_ref(1,1,i),N_int,hij) do j = 1, n_states double precision :: coef,contrib - coef = psi_coef(i,j) !* psi_coef(i,j) - psi_in_out_coef(i,j) = sign(coef,psi_coef(i,j)) * hij + coef = psi_ref_coef(i,j) !* psi_ref_coef(i,j) + psi_in_out_coef(i,j) = sign(coef,psi_ref_coef(i,j)) * hij norm(j,ispin) += psi_in_out_coef(i,j) * psi_in_out_coef(i,j) enddo endif @@ -671,7 +671,7 @@ BEGIN_PROVIDER [ double precision, one_anhil_inact, (n_inact_orb,n_act_orb,N_Sta endif enddo double precision :: norm_bis(N_states,2) - do i = 1, N_det + do i = 1, N_det_ref do j = 1, N_states psi_in_out_coef(i,j) = psi_in_out_coef(i,j) * norm(j,ispin) norm_bis(j,ispin) += psi_in_out_coef(i,j)* psi_in_out_coef(i,j) @@ -684,7 +684,7 @@ BEGIN_PROVIDER [ double precision, one_anhil_inact, (n_inact_orb,n_act_orb,N_Sta do state_target = 1, N_states energies_alpha_beta(state_target, ispin) = 0.d0 if(norm(state_target,ispin) .ne. 0.d0 .and. dabs(norm_no_inv(state_target,ispin)) .gt. thresh_norm)then - call u0_H_dyall_u0_no_exchange(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target) + call u0_H_dyall_u0_no_exchange(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target) energies_alpha_beta(state_target, ispin) += energies(state_target) endif enddo @@ -715,7 +715,7 @@ BEGIN_PROVIDER [ double precision, one_creat_virt, (n_act_orb,n_virt_orb,N_State integer(bit_kind), allocatable :: psi_in_out(:,:,:) double precision, allocatable :: psi_in_out_coef(:,:) use bitmasks - allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states)) + allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det_ref,N_states)) integer :: iorb,jorb,i_ok,aorb,orb_a integer :: state_target @@ -742,8 +742,8 @@ BEGIN_PROVIDER [ double precision, one_creat_virt, (n_act_orb,n_virt_orb,N_State do ispin = 1,2 do i = 1, n_det do j = 1, N_int - psi_in_out(j,1,i) = psi_det(j,1,i) - psi_in_out(j,2,i) = psi_det(j,2,i) + psi_in_out(j,1,i) = psi_ref(j,1,i) + psi_in_out(j,2,i) = psi_ref(j,2,i) enddo call do_mono_excitation(psi_in_out(1,1,i),orb_a,orb_v,ispin,i_ok) if(i_ok.ne.1)then @@ -751,11 +751,11 @@ BEGIN_PROVIDER [ double precision, one_creat_virt, (n_act_orb,n_virt_orb,N_State psi_in_out_coef(i,j) = 0.d0 enddo else - call i_H_j(psi_in_out(1,1,i),psi_det(1,1,i),N_int,hij) + call i_H_j(psi_in_out(1,1,i),psi_ref(1,1,i),N_int,hij) do j = 1, n_states double precision :: coef,contrib - coef = psi_coef(i,j) !* psi_coef(i,j) - psi_in_out_coef(i,j) = sign(coef,psi_coef(i,j)) * hij + coef = psi_ref_coef(i,j) !* psi_ref_coef(i,j) + psi_in_out_coef(i,j) = sign(coef,psi_ref_coef(i,j)) * hij norm(j,ispin) += psi_in_out_coef(i,j) * psi_in_out_coef(i,j) enddo endif @@ -770,7 +770,7 @@ BEGIN_PROVIDER [ double precision, one_creat_virt, (n_act_orb,n_virt_orb,N_State endif enddo double precision :: norm_bis(N_states,2) - do i = 1, N_det + do i = 1, N_det_ref do j = 1, N_states psi_in_out_coef(i,j) = psi_in_out_coef(i,j) * norm(j,ispin) norm_bis(j,ispin) += psi_in_out_coef(i,j)* psi_in_out_coef(i,j) @@ -783,8 +783,8 @@ BEGIN_PROVIDER [ double precision, one_creat_virt, (n_act_orb,n_virt_orb,N_State do state_target = 1, N_states energies_alpha_beta(state_target, ispin) = 0.d0 if(norm(state_target,ispin) .ne. 0.d0 .and. dabs(norm_no_inv(state_target,ispin)) .gt. thresh_norm)then -! call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target) - call u0_H_dyall_u0_no_exchange(energies,psi_in_out,psi_in_out_coef,n_det,n_det,n_det,N_states,state_target) +! call u0_H_dyall_u0(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target) + call u0_H_dyall_u0_no_exchange(energies,psi_in_out,psi_in_out_coef,n_det_ref,n_det_ref,n_det_ref,N_states,state_target) energies_alpha_beta(state_target, ispin) += energies(state_target) endif enddo @@ -812,38 +812,38 @@ END_PROVIDER subroutine give_singles_and_partial_doubles_1h1p_contrib(matrix_1h1p,e_corr_from_1h1p_singles) implicit none - double precision , intent(inout) :: matrix_1h1p(N_det,N_det,N_states) + double precision , intent(inout) :: matrix_1h1p(N_det_ref,N_det_ref,N_states) double precision , intent(out) :: e_corr_from_1h1p_singles(N_states) integer :: i,vorb,j integer :: ispin,jspin integer :: orb_i, hole_particle_i integer :: orb_v - double precision :: norm_out(N_states),diag_elem(N_det),interact_psi0(N_det) + double precision :: norm_out(N_states),diag_elem(N_det_ref),interact_psi0(N_det_ref) double precision :: delta_e_inact_virt(N_states) integer(bit_kind), allocatable :: psi_in_out(:,:,:) double precision, allocatable :: psi_in_out_coef(:,:) double precision, allocatable :: H_matrix(:,:),eigenvectors(:,:),eigenvalues(:),interact_cas(:,:) double precision, allocatable :: delta_e_det(:,:) use bitmasks - allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det,N_states),H_matrix(N_det+1,N_det+1)) - allocate (eigenvectors(size(H_matrix,1),N_det+1)) - allocate (eigenvalues(N_det+1),interact_cas(N_det,N_det)) - allocate (delta_e_det(N_det,N_det)) + allocate (psi_in_out(N_int,2,n_det),psi_in_out_coef(n_det_ref,N_states),H_matrix(N_det_ref+1,N_det_ref+1)) + allocate (eigenvectors(size(H_matrix,1),N_det_ref+1)) + allocate (eigenvalues(N_det_ref+1),interact_cas(N_det_ref,N_det_ref)) + allocate (delta_e_det(N_det_ref,N_det_ref)) integer :: iorb,jorb,i_ok integer :: state_target double precision :: energies(n_states) double precision :: hij double precision :: energies_alpha_beta(N_states,2) - double precision :: lamda_pt2(N_det) + double precision :: lamda_pt2(N_det_ref) double precision :: accu(N_states),norm - double precision :: amplitudes_alpha_beta(N_det,2) - double precision :: delta_e_alpha_beta(N_det,2) + double precision :: amplitudes_alpha_beta(N_det_ref,2) + double precision :: delta_e_alpha_beta(N_det_ref,2) double precision :: coef_array(N_states) - double precision :: coef_perturb(N_det) - double precision :: coef_perturb_bis(N_det) + double precision :: coef_perturb(N_det_ref) + double precision :: coef_perturb_bis(N_det_ref) do vorb = 1,n_virt_orb orb_v = list_virt(vorb) @@ -856,8 +856,8 @@ subroutine give_singles_and_partial_doubles_1h1p_contrib(matrix_1h1p,e_corr_from do ispin = 1,2 do i = 1, n_det do j = 1, N_int - psi_in_out(j,1,i) = psi_det(j,1,i) - psi_in_out(j,2,i) = psi_det(j,2,i) + psi_in_out(j,1,i) = psi_ref(j,1,i) + psi_in_out(j,2,i) = psi_ref(j,2,i) enddo call do_mono_excitation(psi_in_out(1,1,i),orb_i,orb_v,ispin,i_ok) if(i_ok.ne.1)then @@ -866,11 +866,11 @@ subroutine give_singles_and_partial_doubles_1h1p_contrib(matrix_1h1p,e_corr_from print*, 'pb, i_ok ne 0 !!!' endif interact_psi0(i) = 0.d0 - do j = 1 , N_det - call i_H_j(psi_in_out(1,1,i),psi_det(1,1,j),N_int,hij) - call get_delta_e_dyall(psi_det(1,1,j),psi_in_out(1,1,i),coef_array,hij,delta_e_det(i,j)) + do j = 1 , N_det_ref + call i_H_j(psi_in_out(1,1,i),psi_ref(1,1,j),N_int,hij) + call get_delta_e_dyall(psi_ref(1,1,j),psi_in_out(1,1,i),coef_array,hij,delta_e_det(i,j)) interact_cas(i,j) = hij - interact_psi0(i) += hij * psi_coef(j,1) + interact_psi0(i) += hij * psi_ref_coef(j,1) enddo do j = 1, N_int psi_in_out(j,1,i) = psi_active(j,1,i) @@ -882,27 +882,27 @@ subroutine give_singles_and_partial_doubles_1h1p_contrib(matrix_1h1p,e_corr_from do state_target = 1, N_states ! Building the Hamiltonian matrix H_matrix(1,1) = energy_cas_dyall(state_target) - do i = 1, N_det + do i = 1, N_det_ref ! interaction with psi0 - H_matrix(1,i+1) = interact_psi0(i)!* psi_coef(i,state_target) - H_matrix(i+1,1) = interact_psi0(i)!* psi_coef(i,state_target) + H_matrix(1,i+1) = interact_psi0(i)!* psi_ref_coef(i,state_target) + H_matrix(i+1,1) = interact_psi0(i)!* psi_ref_coef(i,state_target) ! diagonal elements H_matrix(i+1,i+1) = diag_elem(i) - delta_e_inact_virt(state_target) ! print*, 'H_matrix(i+1,i+1)',H_matrix(i+1,i+1) - do j = i+1, N_det + do j = i+1, N_det_ref call i_H_j_dyall(psi_in_out(1,1,i),psi_in_out(1,1,j),N_int,hij) H_matrix(i+1,j+1) = hij !0.d0 ! H_matrix(j+1,i+1) = hij !0.d0 ! enddo enddo - call lapack_diag(eigenvalues,eigenvectors,H_matrix,size(H_matrix,1),N_det+1) + call lapack_diag(eigenvalues,eigenvectors,H_matrix,size(H_matrix,1),N_det_ref+1) e_corr_from_1h1p_singles(state_target) += eigenvalues(1) - energy_cas_dyall(state_target) - do i = 1, N_det + do i = 1, N_det_ref psi_in_out_coef(i,state_target) = eigenvectors(i+1,1)/eigenvectors(1,1) coef_perturb(i) = 0.d0 - do j = 1, N_det - coef_perturb(i) += psi_coef(j,state_target) * interact_cas(i,j) *1.d0/delta_e_det(i,j) + do j = 1, N_det_ref + coef_perturb(i) += psi_ref_coef(j,state_target) * interact_cas(i,j) *1.d0/delta_e_det(i,j) enddo coef_perturb_bis(i) = interact_psi0(i) / (eigenvalues(1) - H_matrix(i+1,i+1)) if(dabs(interact_psi0(i)) .gt. 1.d-12)then @@ -913,22 +913,22 @@ subroutine give_singles_and_partial_doubles_1h1p_contrib(matrix_1h1p,e_corr_from enddo if(dabs(eigenvalues(1) - energy_cas_dyall(state_target)).gt.1.d-10)then print*, '' - do i = 1, N_det+1 + do i = 1, N_det_ref+1 write(*,'(100(F16.10))') H_matrix(i,:) enddo accu = 0.d0 - do i = 1, N_det + do i = 1, N_det_ref accu(state_target) += psi_in_out_coef(i,state_target) * interact_psi0(i) enddo print*, '' print*, 'e corr diagonal ',accu(state_target) accu = 0.d0 - do i = 1, N_det + do i = 1, N_det_ref accu(state_target) += coef_perturb(i) * interact_psi0(i) enddo print*, 'e corr perturb ',accu(state_target) accu = 0.d0 - do i = 1, N_det + do i = 1, N_det_ref accu(state_target) += coef_perturb_bis(i) * interact_psi0(i) enddo print*, 'e corr perturb EN',accu(state_target) @@ -941,10 +941,10 @@ subroutine give_singles_and_partial_doubles_1h1p_contrib(matrix_1h1p,e_corr_from write(*,'(100(F16.10,X))')coef_perturb_bis(:) endif integer :: k - do k = 1, N_det - do i = 1, N_det + do k = 1, N_det_ref + do i = 1, N_det_ref matrix_1h1p(i,i,state_target) += interact_cas(k,i) * interact_cas(k,i) * lamda_pt2(k) - do j = i+1, N_det + do j = i+1, N_det_ref matrix_1h1p(i,j,state_target) += interact_cas(k,i) * interact_cas(k,j) * lamda_pt2(k) matrix_1h1p(j,i,state_target) += interact_cas(k,i) * interact_cas(k,j) * lamda_pt2(k) enddo diff --git a/src/MRPT_Utils/excitations_cas.irp.f b/plugins/MRPT_Utils/excitations_cas.irp.f similarity index 100% rename from src/MRPT_Utils/excitations_cas.irp.f rename to plugins/MRPT_Utils/excitations_cas.irp.f diff --git a/plugins/MRPT_Utils/ezfio_interface.irp.f b/plugins/MRPT_Utils/ezfio_interface.irp.f new file mode 100644 index 00000000..ebe0bf52 --- /dev/null +++ b/plugins/MRPT_Utils/ezfio_interface.irp.f @@ -0,0 +1,42 @@ +! DO NOT MODIFY BY HAND +! Created by $QP_ROOT/scripts/ezfio_interface/ei_handler.py +! from file /home/giner/qp_bis/quantum_package/src/MRPT_Utils/EZFIO.cfg + + +BEGIN_PROVIDER [ logical, do_third_order_1h1p ] + implicit none + BEGIN_DOC +! If true, compute the third order contribution for the 1h1p + END_DOC + + logical :: has + PROVIDE ezfio_filename + + call ezfio_has_mrpt_utils_do_third_order_1h1p(has) + if (has) then + call ezfio_get_mrpt_utils_do_third_order_1h1p(do_third_order_1h1p) + else + print *, 'mrpt_utils/do_third_order_1h1p not found in EZFIO file' + stop 1 + endif + +END_PROVIDER + +BEGIN_PROVIDER [ logical, pure_state_specific_mrpt2 ] + implicit none + BEGIN_DOC +! If true, diagonalize the dressed matrix for each state and do a state following of the initial states + END_DOC + + logical :: has + PROVIDE ezfio_filename + + call ezfio_has_mrpt_utils_pure_state_specific_mrpt2(has) + if (has) then + call ezfio_get_mrpt_utils_pure_state_specific_mrpt2(pure_state_specific_mrpt2) + else + print *, 'mrpt_utils/pure_state_specific_mrpt2 not found in EZFIO file' + stop 1 + endif + +END_PROVIDER diff --git a/src/MRPT_Utils/fock_like_operators.irp.f b/plugins/MRPT_Utils/fock_like_operators.irp.f similarity index 100% rename from src/MRPT_Utils/fock_like_operators.irp.f rename to plugins/MRPT_Utils/fock_like_operators.irp.f diff --git a/src/MRPT_Utils/give_2h2p.irp.f b/plugins/MRPT_Utils/give_2h2p.irp.f similarity index 100% rename from src/MRPT_Utils/give_2h2p.irp.f rename to plugins/MRPT_Utils/give_2h2p.irp.f diff --git a/src/MRPT_Utils/mrpt_dress.irp.f b/plugins/MRPT_Utils/mrpt_dress.irp.f similarity index 78% rename from src/MRPT_Utils/mrpt_dress.irp.f rename to plugins/MRPT_Utils/mrpt_dress.irp.f index 60bb2b69..f5e7bd40 100644 --- a/src/MRPT_Utils/mrpt_dress.irp.f +++ b/plugins/MRPT_Utils/mrpt_dress.irp.f @@ -44,15 +44,15 @@ subroutine mrpt_dress(delta_ij_, Ndet,i_generator,n_selected,det_buffer,Nint,ip integer :: N_miniList, leng double precision :: delta_e(N_states),hij_tmp integer :: index_i,index_j - double precision :: phase_array(N_det),phase + double precision :: phase_array(N_det_ref),phase integer :: exc(0:2,2,2),degree - leng = max(N_det_generators, N_det) + leng = max(N_det_ref, N_det_ref) allocate(miniList(Nint, 2, leng), idx_miniList(leng)) !create_minilist_find_previous(key_mask, fullList, miniList, N_fullList, N_miniList, fullMatch, Nint) - call create_minilist_find_previous(key_mask, psi_det_generators, miniList, i_generator-1, N_miniList, fullMatch, Nint) + call create_minilist_find_previous(key_mask, psi_ref, miniList, i_generator-1, N_miniList, fullMatch, Nint) if(fullMatch) then return @@ -62,7 +62,7 @@ subroutine mrpt_dress(delta_ij_, Ndet,i_generator,n_selected,det_buffer,Nint,ip call find_connections_previous(i_generator,n_selected,det_buffer,Nint,tq,N_tq,miniList,N_minilist) if(N_tq > 0) then - call create_minilist(key_mask, psi_det, miniList, idx_miniList, N_det, N_minilist, Nint) + call create_minilist(key_mask, psi_ref, miniList, idx_miniList, N_det_ref, N_minilist, Nint) end if @@ -79,18 +79,18 @@ subroutine mrpt_dress(delta_ij_, Ndet,i_generator,n_selected,det_buffer,Nint,ip phase_array =0.d0 do i = 1,idx_alpha(0) index_i = idx_alpha(i) - call i_h_j(tq(1,1,i_alpha),psi_det(1,1,index_i),Nint,hialpha) + call i_h_j(tq(1,1,i_alpha),psi_ref(1,1,index_i),Nint,hialpha) double precision :: coef_array(N_states) do i_state = 1, N_states - coef_array(i_state) = psi_coef(index_i,i_state) + coef_array(i_state) = psi_ref_coef(index_i,i_state) enddo if(dabs(hialpha).le.1.d-10)then delta_e = 1.d+20 else - call get_delta_e_dyall(psi_det(1,1,index_i),tq(1,1,i_alpha),coef_array,hialpha,delta_e) + call get_delta_e_dyall(psi_ref(1,1,index_i),tq(1,1,i_alpha),coef_array,hialpha,delta_e) endif hij_array(index_i) = hialpha - call get_excitation(psi_det(1,1,index_i),tq(1,1,i_alpha),exc,degree,phase,N_int) + call get_excitation(psi_ref(1,1,index_i),tq(1,1,i_alpha),exc,degree,phase,N_int) ! phase_array(index_i) = phase do i_state = 1,N_states delta_e_inv_array(index_i,i_state) = 1.d0/delta_e(i_state) @@ -103,12 +103,12 @@ subroutine mrpt_dress(delta_ij_, Ndet,i_generator,n_selected,det_buffer,Nint,ip call omp_set_lock( psi_ref_bis_lock(index_i) ) do j = 1, idx_alpha(0) index_j = idx_alpha(j) -! call get_excitation(psi_det(1,1,index_i),psi_det(1,1,index_i),exc,degree,phase,N_int) +! call get_excitation(psi_ref(1,1,index_i),psi_ref(1,1,index_i),exc,degree,phase,N_int) ! if(index_j.ne.index_i)then ! if(phase_array(index_j) * phase_array(index_i) .ne. phase)then ! print*, phase_array(index_j) , phase_array(index_i) ,phase -! call debug_det(psi_det(1,1,index_i),N_int) -! call debug_det(psi_det(1,1,index_j),N_int) +! call debug_det(psi_ref(1,1,index_i),N_int) +! call debug_det(psi_ref(1,1,index_j),N_int) ! call debug_det(tq(1,1,i_alpha),N_int) ! stop ! endif @@ -126,14 +126,14 @@ end - BEGIN_PROVIDER [ integer(bit_kind), gen_det_sorted, (N_int,2,N_det_generators,2) ] -&BEGIN_PROVIDER [ integer, gen_det_shortcut, (0:N_det_generators,2) ] -&BEGIN_PROVIDER [ integer, gen_det_version, (N_int, N_det_generators,2) ] -&BEGIN_PROVIDER [ integer, gen_det_idx, (N_det_generators,2) ] - gen_det_sorted(:,:,:,1) = psi_det_generators(:,:,:N_det_generators) - gen_det_sorted(:,:,:,2) = psi_det_generators(:,:,:N_det_generators) - call sort_dets_ab_v(gen_det_sorted(:,:,:,1), gen_det_idx(:,1), gen_det_shortcut(0:,1), gen_det_version(:,:,1), N_det_generators, N_int) - call sort_dets_ba_v(gen_det_sorted(:,:,:,2), gen_det_idx(:,2), gen_det_shortcut(0:,2), gen_det_version(:,:,2), N_det_generators, N_int) + BEGIN_PROVIDER [ integer(bit_kind), gen_det_sorted, (N_int,2,N_det_ref,2) ] +&BEGIN_PROVIDER [ integer, gen_det_shortcut, (0:N_det_ref,2) ] +&BEGIN_PROVIDER [ integer, gen_det_version, (N_int, N_det_ref,2) ] +&BEGIN_PROVIDER [ integer, gen_det_idx, (N_det_ref,2) ] + gen_det_sorted(:,:,:,1) = psi_ref(:,:,:N_det_ref) + gen_det_sorted(:,:,:,2) = psi_ref(:,:,:N_det_ref) + call sort_dets_ab_v(gen_det_sorted(:,:,:,1), gen_det_idx(:,1), gen_det_shortcut(0:,1), gen_det_version(:,:,1), N_det_ref, N_int) + call sort_dets_ba_v(gen_det_sorted(:,:,:,2), gen_det_idx(:,2), gen_det_shortcut(0:,2), gen_det_version(:,:,2), N_det_ref, N_int) END_PROVIDER @@ -159,7 +159,7 @@ subroutine find_connections_previous(i_generator,n_selected,det_buffer,Nint,tq,N logical, external :: is_connected_to - integer(bit_kind),intent(in) :: miniList(Nint,2,N_det_generators) + integer(bit_kind),intent(in) :: miniList(Nint,2,N_det_ref) integer,intent(in) :: N_miniList @@ -172,7 +172,7 @@ subroutine find_connections_previous(i_generator,n_selected,det_buffer,Nint,tq,N cycle end if - if (.not. is_in_wavefunction(det_buffer(1,1,i),Nint,N_det)) then + if (.not. is_in_wavefunction(det_buffer(1,1,i),Nint,N_det_ref)) then N_tq += 1 do k=1,N_int tq(k,1,N_tq) = det_buffer(k,1,i) diff --git a/src/MRPT_Utils/mrpt_utils.irp.f b/plugins/MRPT_Utils/mrpt_utils.irp.f similarity index 72% rename from src/MRPT_Utils/mrpt_utils.irp.f rename to plugins/MRPT_Utils/mrpt_utils.irp.f index 8ac8e3e0..31013bb7 100644 --- a/src/MRPT_Utils/mrpt_utils.irp.f +++ b/plugins/MRPT_Utils/mrpt_utils.irp.f @@ -1,5 +1,5 @@ - BEGIN_PROVIDER [ double precision, delta_ij, (N_det,N_det,N_states) ] + BEGIN_PROVIDER [ double precision, delta_ij, (N_det_ref,N_det_ref,N_states) ] &BEGIN_PROVIDER [ double precision, second_order_pt_new, (N_states) ] &BEGIN_PROVIDER [ double precision, second_order_pt_new_1h, (N_states) ] &BEGIN_PROVIDER [ double precision, second_order_pt_new_1p, (N_states) ] @@ -11,7 +11,7 @@ &BEGIN_PROVIDER [ double precision, second_order_pt_new_2h2p, (N_states) ] implicit none BEGIN_DOC - ! Dressing matrix in N_det basis + ! Dressing matrix in N_det_ref basis END_DOC integer :: i,j,m integer :: i_state @@ -21,17 +21,17 @@ delta_ij = 0.d0 - allocate (delta_ij_tmp(N_det,N_det,N_states)) + allocate (delta_ij_tmp(N_det_ref,N_det_ref,N_states)) ! 1h delta_ij_tmp = 0.d0 - call H_apply_mrpt_1h(delta_ij_tmp,N_det) + call H_apply_mrpt_1h(delta_ij_tmp,N_det_ref) accu = 0.d0 do i_state = 1, N_states - do i = 1, N_det - do j = 1, N_det - accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) + do i = 1, N_det_ref + do j = 1, N_det_ref + accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_coef(j,i_state) delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) enddo enddo @@ -41,12 +41,12 @@ ! 1p delta_ij_tmp = 0.d0 - call H_apply_mrpt_1p(delta_ij_tmp,N_det) + call H_apply_mrpt_1p(delta_ij_tmp,N_det_ref) accu = 0.d0 do i_state = 1, N_states - do i = 1, N_det - do j = 1, N_det - accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) + do i = 1, N_det_ref + do j = 1, N_det_ref + accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_coef(j,i_state) delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) enddo enddo @@ -56,15 +56,15 @@ ! 1h1p delta_ij_tmp = 0.d0 - call H_apply_mrpt_1h1p(delta_ij_tmp,N_det) + call H_apply_mrpt_1h1p(delta_ij_tmp,N_det_ref) double precision :: e_corr_from_1h1p_singles(N_states) !call give_singles_and_partial_doubles_1h1p_contrib(delta_ij_tmp,e_corr_from_1h1p_singles) !call give_1h1p_only_doubles_spin_cross(delta_ij_tmp) accu = 0.d0 do i_state = 1, N_states - do i = 1, N_det - do j = 1, N_det - accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) + do i = 1, N_det_ref + do j = 1, N_det_ref + accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_coef(j,i_state) delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) enddo enddo @@ -78,9 +78,9 @@ call give_1h1p_sec_order_singles_contrib(delta_ij_tmp) accu = 0.d0 do i_state = 1, N_states - do i = 1, N_det - do j = 1, N_det - accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) + do i = 1, N_det_ref + do j = 1, N_det_ref + accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_coef(j,i_state) delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) enddo enddo @@ -91,12 +91,12 @@ ! 2h delta_ij_tmp = 0.d0 - call H_apply_mrpt_2h(delta_ij_tmp,N_det) + call H_apply_mrpt_2h(delta_ij_tmp,N_det_ref) accu = 0.d0 do i_state = 1, N_states - do i = 1, N_det - do j = 1, N_det - accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) + do i = 1, N_det_ref + do j = 1, N_det_ref + accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_coef(j,i_state) delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) enddo enddo @@ -106,12 +106,12 @@ ! 2p delta_ij_tmp = 0.d0 - call H_apply_mrpt_2p(delta_ij_tmp,N_det) + call H_apply_mrpt_2p(delta_ij_tmp,N_det_ref) accu = 0.d0 do i_state = 1, N_states - do i = 1, N_det - do j = 1, N_det - accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) + do i = 1, N_det_ref + do j = 1, N_det_ref + accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_coef(j,i_state) delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) enddo enddo @@ -122,12 +122,12 @@ ! 1h2p delta_ij_tmp = 0.d0 call give_1h2p_contrib(delta_ij_tmp) -!!!!call H_apply_mrpt_1h2p(delta_ij_tmp,N_det) +!!!!call H_apply_mrpt_1h2p(delta_ij_tmp,N_det_ref) accu = 0.d0 do i_state = 1, N_states - do i = 1, N_det - do j = 1, N_det - accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) + do i = 1, N_det_ref + do j = 1, N_det_ref + accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_coef(j,i_state) delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) enddo enddo @@ -138,12 +138,12 @@ ! 2h1p delta_ij_tmp = 0.d0 call give_2h1p_contrib(delta_ij_tmp) -!!!! call H_apply_mrpt_2h1p(delta_ij_tmp,N_det) + !!!!call H_apply_mrpt_2h1p(delta_ij_tmp,N_det_ref) accu = 0.d0 do i_state = 1, N_states - do i = 1, N_det - do j = 1, N_det - accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) + do i = 1, N_det_ref + do j = 1, N_det_ref + accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_coef(j,i_state) delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) enddo enddo @@ -153,12 +153,12 @@ ! 2h2p delta_ij_tmp = 0.d0 -!!!!!call H_apply_mrpt_2h2p(delta_ij_tmp,N_det) +!!!!!call H_apply_mrpt_2h2p(delta_ij_tmp,N_det_ref) accu = 0.d0 do i_state = 1, N_states - do i = 1, N_det - do j = 1, N_det - accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) + do i = 1, N_det_ref + do j = 1, N_det_ref + accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_coef(j,i_state) delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) enddo enddo @@ -169,7 +169,7 @@ double precision :: contrib_2h2p(N_states) call give_2h2p(contrib_2h2p) do i_state = 1, N_states - do i = 1, N_det + do i = 1, N_det_ref delta_ij(i,i,i_state) += contrib_2h2p(i_state) enddo second_order_pt_new_2h2p(i_state) = contrib_2h2p(i_state) @@ -180,10 +180,10 @@ ! total accu = 0.d0 do i_state = 1, N_states - do i = 1, N_det + do i = 1, N_det_ref ! write(*,'(1000(F16.10,x))')delta_ij(i,:,:) - do j = i_state, N_det - accu(i_state) += delta_ij(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) + do j = i_state, N_det_ref + accu(i_state) += delta_ij(j,i,i_state) * psi_ref_coef(i,i_state) * psi_ref_coef(j,i_state) enddo enddo second_order_pt_new(i_state) = accu(i_state) @@ -195,13 +195,15 @@ END_PROVIDER - BEGIN_PROVIDER [double precision, Hmatrix_dressed_pt2_new, (N_det,N_det,N_states)] + BEGIN_PROVIDER [double precision, Hmatrix_dressed_pt2_new, (N_det_ref,N_det_ref,N_states)] implicit none integer :: i,j,i_state + double precision :: hij do i_state = 1, N_states - do i = 1,N_det - do j = 1,N_det - Hmatrix_dressed_pt2_new(j,i,i_state) = H_matrix_all_dets(j,i) + delta_ij(j,i,i_state) + do i = 1,N_det_ref + do j = 1,N_det_ref + call i_h_j(psi_ref(1,1,j),psi_ref(1,1,i),N_int,hij) + Hmatrix_dressed_pt2_new(j,i,i_state) = hij + delta_ij(j,i,i_state) enddo enddo enddo @@ -209,13 +211,15 @@ END_PROVIDER - BEGIN_PROVIDER [double precision, Hmatrix_dressed_pt2_new_symmetrized, (N_det,N_det,N_states)] + BEGIN_PROVIDER [double precision, Hmatrix_dressed_pt2_new_symmetrized, (N_det_ref,N_det_ref,N_states)] implicit none integer :: i,j,i_state + double precision :: hij do i_state = 1, N_states - do i = 1,N_det - do j = i,N_det - Hmatrix_dressed_pt2_new_symmetrized(j,i,i_state) = H_matrix_all_dets(j,i) & + do i = 1,N_det_ref + do j = i,N_det_ref + call i_h_j(psi_ref(1,1,j),psi_ref(1,1,i),N_int,hij) + Hmatrix_dressed_pt2_new_symmetrized(j,i,i_state) = hij & + 0.5d0 * ( delta_ij(j,i,i_state) + delta_ij(i,j,i_state) ) Hmatrix_dressed_pt2_new_symmetrized(i,j,i_state) = Hmatrix_dressed_pt2_new_symmetrized(j,i,i_state) enddo @@ -224,7 +228,7 @@ END_PROVIDER END_PROVIDER BEGIN_PROVIDER [ double precision, CI_electronic_dressed_pt2_new_energy, (N_states_diag) ] - &BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_eigenvectors, (N_det,N_states_diag) ] + &BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_eigenvectors, (N_det_ref,N_states_diag) ] &BEGIN_PROVIDER [ double precision, CI_dressed_pt2_new_eigenvectors_s2, (N_states_diag) ] BEGIN_DOC ! Eigenvectors/values of the CI matrix @@ -243,18 +247,18 @@ END_PROVIDER double precision, allocatable :: s2_eigvalues(:) double precision, allocatable :: e_array(:) integer, allocatable :: iorder(:) - double precision :: overlap(N_det) + double precision :: overlap(N_det_ref) double precision, allocatable :: psi_tmp(:) ! Guess values for the "N_states_diag" states of the CI_dressed_pt2_new_eigenvectors - do j=1,min(N_states,N_det) - do i=1,N_det - CI_dressed_pt2_new_eigenvectors(i,j) = psi_coef(i,j) + do j=1,min(N_states,N_det_ref) + do i=1,N_det_ref + CI_dressed_pt2_new_eigenvectors(i,j) = psi_ref_coef(i,j) enddo enddo - do j=min(N_states,N_det)+1,N_states_diag - do i=1,N_det + do j=min(N_states,N_det_ref)+1,N_states_diag + do i=1,N_det_ref CI_dressed_pt2_new_eigenvectors(i,j) = 0.d0 enddo enddo @@ -265,33 +269,33 @@ END_PROVIDER stop else if (diag_algorithm == "Lapack") then - allocate (eigenvectors(N_det,N_det)) - allocate (eigenvalues(N_det)) + allocate (eigenvectors(N_det_ref,N_det_ref)) + allocate (eigenvalues(N_det_ref)) if(pure_state_specific_mrpt2)then - allocate (hmatrix_tmp(N_det,N_det)) - allocate (iorder(N_det)) - allocate (psi_tmp(N_det)) + allocate (hmatrix_tmp(N_det_ref,N_det_ref)) + allocate (iorder(N_det_ref)) + allocate (psi_tmp(N_det_ref)) print*,'' print*,'***************************' do i_state = 1, N_states !! Big loop over states print*,'' print*,'Diagonalizing with the dressing for state',i_state - do i = 1, N_det - do j = 1, N_det + do i = 1, N_det_ref + do j = 1, N_det_ref hmatrix_tmp(j,i) = Hmatrix_dressed_pt2_new_symmetrized(j,i,i_state) enddo enddo call lapack_diag(eigenvalues,eigenvectors, & - Hmatrix_dressed_pt2_new_symmetrized(1,1,1),N_det,N_det) + Hmatrix_dressed_pt2_new_symmetrized(1,1,1),N_det_ref,N_det_ref) write(*,'(A86)')'Looking for the most overlapping state within all eigenvectors of the dressed matrix' print*,'' print*,'Calculating the overlap for ...' - do i = 1, N_det + do i = 1, N_det_ref overlap(i) = 0.d0 iorder(i) = i print*,'eigenvector',i - do j = 1, N_det - overlap(i)+= psi_coef(j,i_state) * eigenvectors(j,i) + do j = 1, N_det_ref + overlap(i)+= psi_ref_coef(j,i_state) * eigenvectors(j,i) enddo overlap(i) = -dabs(overlap(i)) print*,'energy = ',eigenvalues(i) + nuclear_repulsion @@ -305,26 +309,26 @@ END_PROVIDER print*,'with the overlap of ',dabs(overlap(1)) print*,'and an energy of ',eigenvalues(iorder(1)) + nuclear_repulsion print*,'Calculating the S^2 value ...' - do i=1,N_det + do i=1,N_det_ref CI_dressed_pt2_new_eigenvectors(i,i_state) = eigenvectors(i,iorder(1)) psi_tmp(i) = eigenvectors(i,iorder(1)) enddo CI_electronic_dressed_pt2_new_energy(i_state) = eigenvalues(iorder(1)) - call u_0_S2_u_0(CI_dressed_pt2_new_eigenvectors_s2(i_state),psi_tmp,N_det,psi_det,N_int,1,N_det) + call u_0_S2_u_0(CI_dressed_pt2_new_eigenvectors_s2(i_state),psi_tmp,N_det_ref,psi_det,N_int,1,N_det_ref) print*,'S^2 = ', CI_dressed_pt2_new_eigenvectors_s2(i_state) enddo else call lapack_diag(eigenvalues,eigenvectors, & - Hmatrix_dressed_pt2_new_symmetrized(1,1,1),N_det,N_det) + Hmatrix_dressed_pt2_new_symmetrized(1,1,1),N_det_ref,N_det_ref) CI_electronic_dressed_pt2_new_energy(:) = 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)) + allocate (s2_eigvalues(N_det_ref)) + allocate(index_good_state_array(N_det_ref),good_state_array(N_det_ref)) 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 + call u_0_S2_u_0(s2_eigvalues,eigenvectors,N_det_ref,psi_det,N_int,& + N_det_ref,size(eigenvectors,1)) + do j=1,N_det_ref ! Select at least n_states states with S^2 values closed to "expected_s2" print*, eigenvalues(j)+nuclear_repulsion, s2_eigvalues(j) if(dabs(s2_eigvalues(j)-expected_s2).le.0.5d0)then @@ -339,20 +343,20 @@ END_PROVIDER 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 + do i=1,N_det_ref CI_dressed_pt2_new_eigenvectors(i,j) = eigenvectors(i,index_good_state_array(j)) enddo CI_electronic_dressed_pt2_new_energy(j) = eigenvalues(index_good_state_array(j)) CI_dressed_pt2_new_eigenvectors_s2(j) = s2_eigvalues(index_good_state_array(j)) enddo i_other_state = 0 - do j = 1, N_det + do j = 1, N_det_ref if(good_state_array(j))cycle i_other_state +=1 if(i_state+i_other_state.gt.n_states)then exit endif - do i=1,N_det + do i=1,N_det_ref CI_dressed_pt2_new_eigenvectors(i,i_state+i_other_state) = eigenvectors(i,j) enddo CI_electronic_dressed_pt2_new_energy(i_state+i_other_state) = eigenvalues(j) @@ -362,15 +366,15 @@ END_PROVIDER else print*,'' print*,'!!!!!!!! WARNING !!!!!!!!!' - print*,' Within the ',N_det,'determinants selected' + print*,' Within the ',N_det_ref,'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_dressed_pt2_new_eigenvectors' 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 + do j=1,min(N_states_diag,N_det_ref) + do i=1,N_det_ref CI_dressed_pt2_new_eigenvectors(i,j) = eigenvectors(i,j) enddo CI_electronic_dressed_pt2_new_energy(j) = eigenvalues(j) @@ -380,11 +384,11 @@ END_PROVIDER deallocate(index_good_state_array,good_state_array) deallocate(s2_eigvalues) else - call u_0_S2_u_0(CI_dressed_pt2_new_eigenvectors_s2,eigenvectors,N_det,psi_det,N_int,& - min(N_det,N_states_diag),size(eigenvectors,1)) + call u_0_S2_u_0(CI_dressed_pt2_new_eigenvectors_s2,eigenvectors,N_det_ref,psi_det,N_int,& + min(N_det_ref,N_states_diag),size(eigenvectors,1)) ! Select the "N_states_diag" states of lowest energy - do j=1,min(N_det,N_states) - do i=1,N_det + do j=1,min(N_det_ref,N_states) + do i=1,N_det_ref CI_dressed_pt2_new_eigenvectors(i,j) = eigenvectors(i,j) enddo CI_electronic_dressed_pt2_new_energy(j) = eigenvalues(j) diff --git a/src/MRPT_Utils/new_way.irp.f b/plugins/MRPT_Utils/new_way.irp.f similarity index 98% rename from src/MRPT_Utils/new_way.irp.f rename to plugins/MRPT_Utils/new_way.irp.f index 3624b7d3..a4bbe93a 100644 --- a/src/MRPT_Utils/new_way.irp.f +++ b/plugins/MRPT_Utils/new_way.irp.f @@ -38,8 +38,8 @@ subroutine give_2h1p_contrib(matrix_2h1p) active_int(a,2) = get_mo_bielec_integral(iorb,jorb,aorb,rorb,mo_integrals_map) ! exchange enddo - integer :: degree(N_det) - integer :: idx(0:N_det) + integer :: degree(N_det_Ref) + integer :: idx(0:N_det_Ref) double precision :: delta_e(n_act_orb,2,N_states) integer :: istate integer :: index_orb_act_mono(N_det,3) @@ -232,8 +232,8 @@ subroutine give_1h2p_contrib(matrix_1h2p) active_int(a,2) = get_mo_bielec_integral(iorb,aorb,vorb,rorb,mo_integrals_map) ! exchange enddo - integer :: degree(N_det) - integer :: idx(0:N_det) + integer :: degree(N_det_Ref) + integer :: idx(0:N_det_Ref) double precision :: delta_e(n_act_orb,2,N_states) integer :: istate integer :: index_orb_act_mono(N_det,3) @@ -413,8 +413,8 @@ subroutine give_1h1p_sec_order_singles_contrib(matrix_1h1p) double precision :: get_mo_bielec_integral double precision :: active_int(n_act_orb,2) double precision :: hij,phase - integer :: degree(N_det) - integer :: idx(0:N_det) + integer :: degree(N_det_Ref) + integer :: idx(0:N_det_Ref) integer :: istate double precision :: hja,delta_e_inact_virt(N_states) integer :: kspin,degree_scalar @@ -572,8 +572,8 @@ subroutine give_1p_sec_order_singles_contrib(matrix_1p) integer :: accu_elec double precision :: get_mo_bielec_integral double precision :: hij,phase - integer :: degree(N_det) - integer :: idx(0:N_det) + integer :: degree(N_det_Ref) + integer :: idx(0:N_det_Ref) integer :: istate double precision :: hja,delta_e_act_virt(N_states) integer :: kspin,degree_scalar @@ -715,8 +715,8 @@ subroutine give_1h1p_only_doubles_spin_cross(matrix_1h1p) double precision :: get_mo_bielec_integral double precision :: active_int(n_act_orb,2) double precision :: hij,phase - integer :: degree(N_det) - integer :: idx(0:N_det) + integer :: degree(N_det_Ref) + integer :: idx(0:N_det_Ref) integer :: istate double precision :: hja,delta_e_inact_virt(N_states) integer(bit_kind) :: pert_det(N_int,2,n_act_orb,n_act_orb,2) diff --git a/src/MRPT_Utils/new_way_second_order_coef.irp.f b/plugins/MRPT_Utils/new_way_second_order_coef.irp.f similarity index 99% rename from src/MRPT_Utils/new_way_second_order_coef.irp.f rename to plugins/MRPT_Utils/new_way_second_order_coef.irp.f index 4c12dbe1..676e14e9 100644 --- a/src/MRPT_Utils/new_way_second_order_coef.irp.f +++ b/plugins/MRPT_Utils/new_way_second_order_coef.irp.f @@ -44,8 +44,8 @@ subroutine give_2h1p_contrib_sec_order(matrix_2h1p) perturb_dets_phase(a,2,1) = -1000.d0 enddo - integer :: degree(N_det) - integer :: idx(0:N_det) + integer :: degree(N_det_Ref) + integer :: idx(0:N_det_Ref) double precision :: delta_e(n_act_orb,2,N_states) integer :: istate @@ -379,8 +379,8 @@ subroutine give_1h2p_contrib_sec_order(matrix_1h2p) double precision :: active_int(n_act_orb,2) double precision :: hij,phase double precision :: accu_contrib - integer :: degree(N_det) - integer :: idx(0:N_det) + integer :: degree(N_det_Ref) + integer :: idx(0:N_det_Ref) double precision :: delta_e(n_act_orb,2,N_states) integer :: istate integer :: index_orb_act_mono(N_det,6) diff --git a/src/MRPT_Utils/psi_active_prov.irp.f b/plugins/MRPT_Utils/psi_active_prov.irp.f similarity index 100% rename from src/MRPT_Utils/psi_active_prov.irp.f rename to plugins/MRPT_Utils/psi_active_prov.irp.f diff --git a/src/MRPT_Utils/second_order_new.irp.f b/plugins/MRPT_Utils/second_order_new.irp.f similarity index 99% rename from src/MRPT_Utils/second_order_new.irp.f rename to plugins/MRPT_Utils/second_order_new.irp.f index ba3b421b..2a61eece 100644 --- a/src/MRPT_Utils/second_order_new.irp.f +++ b/plugins/MRPT_Utils/second_order_new.irp.f @@ -22,8 +22,8 @@ subroutine give_1h2p_new(matrix_1h2p) double precision :: active_int(n_act_orb,2) double precision :: hij,phase double precision :: accu_contrib(N_states) - integer :: degree(N_det) - integer :: idx(0:N_det) + integer :: degree(N_det_Ref) + integer :: idx(0:N_det_Ref) double precision :: delta_e(n_act_orb,2,N_states) double precision :: delta_e_inv(n_act_orb,2,N_states) double precision :: delta_e_inactive_virt(N_states) @@ -502,8 +502,8 @@ subroutine give_2h1p_new(matrix_2h1p) double precision :: delta_e_inv(n_act_orb,2,N_states) double precision :: fock_operator_local(n_act_orb,n_act_orb,2) double precision :: delta_e_inactive_virt(N_states) - integer :: degree(N_det) - integer :: idx(0:N_det) + integer :: degree(N_det_Ref) + integer :: idx(0:N_det_Ref) double precision :: delta_e(n_act_orb,2,N_states) integer :: istate integer :: index_orb_act_mono(N_det,3) diff --git a/src/MRPT_Utils/second_order_new_2p.irp.f b/plugins/MRPT_Utils/second_order_new_2p.irp.f similarity index 99% rename from src/MRPT_Utils/second_order_new_2p.irp.f rename to plugins/MRPT_Utils/second_order_new_2p.irp.f index 11ae18da..d086b6c5 100644 --- a/src/MRPT_Utils/second_order_new_2p.irp.f +++ b/plugins/MRPT_Utils/second_order_new_2p.irp.f @@ -21,8 +21,8 @@ subroutine give_2p_new(matrix_2p) double precision :: active_int(n_act_orb,n_act_orb,2) double precision :: hij,phase double precision :: accu_contrib(N_states) - integer :: degree(N_det) - integer :: idx(0:N_det) + integer :: degree(N_det_Ref) + integer :: idx(0:N_det_Ref) double precision :: delta_e(n_act_orb,n_act_orb,2,2,N_states) double precision :: delta_e_inv(n_act_orb,n_act_orb,2,2,N_states) double precision :: delta_e_inactive_virt(N_states) diff --git a/src/MRPT_Utils/utils_bitmask.irp.f b/plugins/MRPT_Utils/utils_bitmask.irp.f similarity index 100% rename from src/MRPT_Utils/utils_bitmask.irp.f rename to plugins/MRPT_Utils/utils_bitmask.irp.f 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/scripts/generate_h_apply.py b/scripts/generate_h_apply.py index c7714e8a..5dd1e4f3 100755 --- a/scripts/generate_h_apply.py +++ b/scripts/generate_h_apply.py @@ -343,7 +343,7 @@ class H_apply(object): """ self.data["size_max"] = "8192" self.data["initialization"] = """ - PROVIDE psi_selectors_coef psi_selectors E_corr_per_selectors psi_det_sorted_bit +! PROVIDE psi_selectors_coef psi_selectors E_corr_per_selectors psi_det_sorted_bit """ if self.do_double_exc == True: self.data["keys_work"] = """ @@ -370,7 +370,7 @@ class H_apply(object): double precision, intent(inout):: norm_pert(N_st) double precision, intent(inout):: H_pert_diag(N_st) double precision :: delta_pt2(N_st), norm_psi(N_st), pt2_old(N_st) - PROVIDE N_det_generators +! PROVIDE N_det_generators do k=1,N_st pt2(k) = 0.d0 norm_pert(k) = 0.d0 @@ -478,7 +478,7 @@ class H_apply_zmq(H_apply): double precision, intent(inout):: norm_pert(N_st) double precision, intent(inout):: H_pert_diag(N_st) double precision :: delta_pt2(N_st), norm_psi(N_st), pt2_old(N_st) - PROVIDE N_det_generators +! PROVIDE N_det_generators do k=1,N_st pt2(k) = 0.d0 norm_pert(k) = 0.d0 diff --git a/src/Determinants/H_apply_nozmq.template.f b/src/Determinants/H_apply_nozmq.template.f index 0c319fe3..5550d9d1 100644 --- a/src/Determinants/H_apply_nozmq.template.f +++ b/src/Determinants/H_apply_nozmq.template.f @@ -17,7 +17,7 @@ subroutine $subroutine($params_main) double precision, allocatable :: fock_diag_tmp(:,:) $initialization - PROVIDE H_apply_buffer_allocated mo_bielec_integrals_in_map psi_det_generators psi_coef_generators + PROVIDE H_apply_buffer_allocated mo_bielec_integrals_in_map !psi_det_generators psi_coef_generators nmax = mod( N_det_generators,nproc ) diff --git a/src/Determinants/H_apply_zmq.template.f b/src/Determinants/H_apply_zmq.template.f index 59544b79..3e4c1867 100644 --- a/src/Determinants/H_apply_zmq.template.f +++ b/src/Determinants/H_apply_zmq.template.f @@ -20,7 +20,7 @@ subroutine $subroutine($params_main) double precision, allocatable :: fock_diag_tmp(:,:) $initialization - PROVIDE H_apply_buffer_allocated mo_bielec_integrals_in_map psi_det_generators psi_coef_generators +! PROVIDE H_apply_buffer_allocated mo_bielec_integrals_in_map psi_det_generators psi_coef_generators integer(ZMQ_PTR), external :: new_zmq_pair_socket integer(ZMQ_PTR) :: zmq_socket_pair