From 2865ff8a701c005d3bdb3e0970f8db672b2489fd Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 19 Dec 2024 09:36:32 +0100 Subject: [PATCH 01/10] Fixed qp set_file autocompletion --- etc/qp.rc | 17 ++--------------- 1 file changed, 2 insertions(+), 15 deletions(-) diff --git a/etc/qp.rc b/etc/qp.rc index adabc404..4e206096 100644 --- a/etc/qp.rc +++ b/etc/qp.rc @@ -190,21 +190,8 @@ _qp_Complete() ;; esac;; set_file) - # Caching the search results to reduce repeated find calls - if [[ -z "$QP_FILE_CACHE" || "$CACHE_DIR" != "$PWD" ]]; then - CACHE_DIR="$PWD" - QP_FILE_CACHE=$(find . -type f -name .version -exec dirname {} \; | sed 's/\/\.version$//') - fi - - # Support for relative paths - prefix=$(dirname "${cur}") - if [[ "$prefix" != "." ]]; then - dirs=$(echo "$QP_FILE_CACHE" | grep "^$prefix") - else - dirs="$QP_FILE_CACHE" - fi - - COMPREPLY=( $(compgen -W "$dirs" -- "$cur") ) + compopt -o nospace + COMPREPLY=( $(compgen -d -- "$cur") ) return 0 ;; plugins) From b0191f2c72a3053abb6a9df644ea97274f3d11ae Mon Sep 17 00:00:00 2001 From: Yann Damour <77277447+Ydrnan@users.noreply.github.com> Date: Fri, 17 Jan 2025 04:24:27 +0100 Subject: [PATCH 02/10] Update qp_import_trexio.py If True the export under trexio format followed by the import under trexio format leads to a wrong energy --- scripts/qp_import_trexio.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/qp_import_trexio.py b/scripts/qp_import_trexio.py index 630ffb5e..fc76f8de 100755 --- a/scripts/qp_import_trexio.py +++ b/scripts/qp_import_trexio.py @@ -195,7 +195,7 @@ def write_ezfio(trexio_filename, filename): prim_factor = trexio.read_basis_prim_factor(trexio_file) for i,p in enumerate(prim_factor): coefficient[i] *= prim_factor[i] - ezfio.set_ao_basis_primitives_normalized(True) + ezfio.set_ao_basis_primitives_normalized(False) ezfio.set_basis_prim_coef(coefficient) elif basis_type.lower() == "numerical": From 09f6d338e88dfc66ecb598f1c719322c6ca3ca22 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 20 Jan 2025 17:56:13 +0100 Subject: [PATCH 03/10] provide all_mo_integrals --- src/davidson/diagonalization_h_dressed.irp.f | 2 +- .../diagonalization_hcsf_dressed.irp.f | 2 +- .../diagonalization_hs2_dressed.irp.f | 2 +- .../diagonalization_nonsym_h_dressed.irp.f | 2 +- src/determinants/h_apply_nozmq.template.f | 2 +- src/determinants/slater_rules.irp.f | 12 ++++----- src/determinants/slater_rules_wee_mono.irp.f | 4 +-- src/determinants/utils.irp.f | 9 +++---- src/mo_two_e_ints/map_integrals.irp.f | 25 ++++++++++++++----- 9 files changed, 36 insertions(+), 24 deletions(-) diff --git a/src/davidson/diagonalization_h_dressed.irp.f b/src/davidson/diagonalization_h_dressed.irp.f index 15bf256d..a7e501ea 100644 --- a/src/davidson/diagonalization_h_dressed.irp.f +++ b/src/davidson/diagonalization_h_dressed.irp.f @@ -31,7 +31,7 @@ subroutine davidson_diag_h(dets_in,u_in,dim_in,energies,sze,N_st,N_st_diag,Nint, ASSERT (sze > 0) ASSERT (Nint > 0) ASSERT (Nint == N_int) - PROVIDE mo_two_e_integrals_in_map + PROVIDE all_mo_integrals allocate(H_jj(sze)) H_jj(1) = diag_h_mat_elem(dets_in(1,1,1),Nint) diff --git a/src/davidson/diagonalization_hcsf_dressed.irp.f b/src/davidson/diagonalization_hcsf_dressed.irp.f index 656dd1d9..1fb24e08 100644 --- a/src/davidson/diagonalization_hcsf_dressed.irp.f +++ b/src/davidson/diagonalization_hcsf_dressed.irp.f @@ -30,7 +30,7 @@ subroutine davidson_diag_h_csf(dets_in,u_in,dim_in,energies,sze,sze_csf,N_st,N_s ASSERT (sze > 0) ASSERT (Nint > 0) ASSERT (Nint == N_int) - PROVIDE mo_two_e_integrals_in_map + PROVIDE all_mo_integrals allocate(H_jj(sze)) H_jj(1) = diag_h_mat_elem(dets_in(1,1,1),Nint) diff --git a/src/davidson/diagonalization_hs2_dressed.irp.f b/src/davidson/diagonalization_hs2_dressed.irp.f index d299f982..ce2cb63f 100644 --- a/src/davidson/diagonalization_hs2_dressed.irp.f +++ b/src/davidson/diagonalization_hs2_dressed.irp.f @@ -62,7 +62,7 @@ subroutine davidson_diag_hs2(dets_in,u_in,s2_out,dim_in,energies,sze,N_st,N_st_d ASSERT (sze > 0) ASSERT (Nint > 0) ASSERT (Nint == N_int) - PROVIDE mo_two_e_integrals_in_map + PROVIDE all_mo_integrals allocate(H_jj(sze)) H_jj(1) = diag_h_mat_elem(dets_in(1,1,1),Nint) diff --git a/src/davidson/diagonalization_nonsym_h_dressed.irp.f b/src/davidson/diagonalization_nonsym_h_dressed.irp.f index 86df3a19..3ac37f1d 100644 --- a/src/davidson/diagonalization_nonsym_h_dressed.irp.f +++ b/src/davidson/diagonalization_nonsym_h_dressed.irp.f @@ -42,7 +42,7 @@ subroutine davidson_diag_nonsym_h(dets_in, u_in, dim_in, energies, sze, N_st, N_ ASSERT (sze > 0) ASSERT (Nint > 0) ASSERT (Nint == N_int) - PROVIDE mo_two_e_integrals_in_map + PROVIDE all_mo_integrals allocate(H_jj(sze)) diff --git a/src/determinants/h_apply_nozmq.template.f b/src/determinants/h_apply_nozmq.template.f index bd261bbe..3463a818 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_two_e_integrals_in_map psi_det_generators psi_coef_generators + PROVIDE H_apply_buffer_allocated all_mo_integrals psi_det_generators psi_coef_generators call wall_time(wall_0) diff --git a/src/determinants/slater_rules.irp.f b/src/determinants/slater_rules.irp.f index 3a33a37d..5c170b9e 100644 --- a/src/determinants/slater_rules.irp.f +++ b/src/determinants/slater_rules.irp.f @@ -521,7 +521,7 @@ subroutine i_H_j_s2(key_i,key_j,Nint,hij,s2) integer :: occ(Nint*bit_kind_size,2) double precision :: diag_H_mat_elem, phase integer :: n_occ_ab(2) - PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals + PROVIDE all_mo_integrals ASSERT (Nint > 0) ASSERT (Nint == N_int) @@ -623,7 +623,7 @@ subroutine i_H_j(key_i,key_j,Nint,hij) integer :: occ(Nint*bit_kind_size,2) double precision :: diag_H_mat_elem, phase integer :: n_occ_ab(2) - PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals + PROVIDE all_mo_integrals ASSERT (Nint > 0) ASSERT (Nint == N_int) @@ -724,7 +724,7 @@ subroutine i_H_j_verbose(key_i,key_j,Nint,hij,hmono,hdouble,phase) integer :: n_occ_ab(2) logical :: has_mipi(Nint*bit_kind_size) double precision :: mipi(Nint*bit_kind_size), miip(Nint*bit_kind_size) - PROVIDE mo_two_e_integrals_in_map mo_integrals_map + PROVIDE all_mo_integrals ASSERT (Nint > 0) ASSERT (Nint == N_int) @@ -2227,7 +2227,7 @@ subroutine i_H_j_single_spin(key_i,key_j,Nint,spin,hij) integer :: exc(0:2,2) double precision :: phase - PROVIDE big_array_exchange_integrals mo_two_e_integrals_in_map + PROVIDE all_mo_integrals call get_single_excitation_spin(key_i(1,spin),key_j(1,spin),exc,phase,Nint) call get_single_excitation_from_fock(key_i,key_j,exc(1,1),exc(1,2),spin,phase,hij) @@ -2248,7 +2248,7 @@ subroutine i_H_j_double_spin(key_i,key_j,Nint,hij) double precision :: phase double precision, external :: get_two_e_integral - PROVIDE big_array_exchange_integrals mo_two_e_integrals_in_map + PROVIDE all_mo_integrals call get_double_excitation_spin(key_i,key_j,exc,phase,Nint) hij = phase*(get_two_e_integral( & exc(1,1), & @@ -2277,7 +2277,7 @@ subroutine i_H_j_double_alpha_beta(key_i,key_j,Nint,hij) double precision :: phase, phase2 double precision, external :: get_two_e_integral - PROVIDE big_array_exchange_integrals mo_two_e_integrals_in_map + PROVIDE all_mo_integrals call get_single_excitation_spin(key_i(1,1),key_j(1,1),exc(0,1,1),phase,Nint) call get_single_excitation_spin(key_i(1,2),key_j(1,2),exc(0,1,2),phase2,Nint) diff --git a/src/determinants/slater_rules_wee_mono.irp.f b/src/determinants/slater_rules_wee_mono.irp.f index 4c1c9330..b94927e3 100644 --- a/src/determinants/slater_rules_wee_mono.irp.f +++ b/src/determinants/slater_rules_wee_mono.irp.f @@ -13,7 +13,7 @@ subroutine i_Wee_j_single(key_i,key_j,Nint,spin,hij) integer :: exc(0:2,2) double precision :: phase - PROVIDE big_array_exchange_integrals mo_two_e_integrals_in_map + PROVIDE all_mo_integrals call get_single_excitation_spin(key_i(1,spin),key_j(1,spin),exc,phase,Nint) call single_excitation_wee(key_i,key_j,exc(1,1),exc(1,2),spin,phase,hij) @@ -285,7 +285,7 @@ subroutine i_H_j_two_e(key_i,key_j,Nint,hij) integer :: occ(Nint*bit_kind_size,2) double precision :: diag_H_mat_elem, phase,phase_2 integer :: n_occ_ab(2) - PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals ref_bitmask_two_e_energy + PROVIDE all_mo_integrals ref_bitmask_two_e_energy ASSERT (Nint > 0) ASSERT (Nint == N_int) diff --git a/src/determinants/utils.irp.f b/src/determinants/utils.irp.f index 7b75d985..07a6334b 100644 --- a/src/determinants/utils.irp.f +++ b/src/determinants/utils.irp.f @@ -4,12 +4,11 @@ BEGIN_PROVIDER [ double precision, H_matrix_all_dets,(N_det,N_det) ] BEGIN_DOC ! |H| matrix on the basis of the Slater determinants defined by psi_det END_DOC - integer :: i,j,k + integer :: i,j double precision :: hij - integer :: degree(N_det),idx(0:N_det) - call i_H_j(psi_det(1,1,1),psi_det(1,1,1),N_int,hij) + PROVIDE all_mo_integrals print*,'Providing the H_matrix_all_dets ...' - !$OMP PARALLEL DO SCHEDULE(GUIDED) DEFAULT(NONE) PRIVATE(i,j,hij,degree,idx,k) & + !$OMP PARALLEL DO SCHEDULE(GUIDED) DEFAULT(NONE) PRIVATE(i,j,hij) & !$OMP SHARED (N_det, psi_det, N_int,H_matrix_all_dets) do i =1,N_det do j = i, N_det @@ -31,7 +30,7 @@ BEGIN_PROVIDER [ double precision, H_matrix_diag_all_dets,(N_det) ] integer :: i double precision :: hij integer :: degree(N_det) - call i_H_j(psi_det(1,1,1),psi_det(1,1,1),N_int,hij) + PROVIDE all_mo_integrals !$OMP PARALLEL DO SCHEDULE(GUIDED) DEFAULT(NONE) PRIVATE(i,hij,degree) & !$OMP SHARED (N_det, psi_det, N_int,H_matrix_diag_all_dets) do i =1,N_det diff --git a/src/mo_two_e_ints/map_integrals.irp.f b/src/mo_two_e_ints/map_integrals.irp.f index eeb4279f..e6b2967a 100644 --- a/src/mo_two_e_ints/map_integrals.irp.f +++ b/src/mo_two_e_ints/map_integrals.irp.f @@ -1,5 +1,15 @@ use map_module +BEGIN_PROVIDER [ logical, all_mo_integrals ] + implicit none + BEGIN_DOC +! Used to provide everything needed before using MO integrals +! PROVIDE all_mo_integrals + END_DOC + PROVIDE mo_two_e_integrals_in_map mo_integrals_cache mo_two_e_integrals_jj_exchange mo_two_e_integrals_jj_anti mo_two_e_integrals_jj big_array_exchange_integrals big_array_coulomb_integrals +END_PROVIDER + + !! MO Map !! ====== @@ -35,20 +45,24 @@ end BEGIN_PROVIDER [ integer, mo_integrals_cache_min ] &BEGIN_PROVIDER [ integer, mo_integrals_cache_max ] &BEGIN_PROVIDER [ integer, mo_integrals_cache_size ] +&BEGIN_PROVIDER [ integer*8, mo_integrals_cache_size_8 ] implicit none BEGIN_DOC ! Min and max values of the MOs for which the integrals are in the cache END_DOC - mo_integrals_cache_size = 2**mo_integrals_cache_shift + mo_integrals_cache_size = shiftl(1,mo_integrals_cache_shift) + mo_integrals_cache_size_8 = shiftl(1_8, mo_integrals_cache_shift*4) + mo_integrals_cache_min = max(1,elec_alpha_num - (mo_integrals_cache_size/2 - 1) ) mo_integrals_cache_max = min(mo_num, mo_integrals_cache_min + mo_integrals_cache_size - 1) - print *, 'MO integrals cache: (', mo_integrals_cache_min, ', ', mo_integrals_cache_max, ')' + print *, 'MO integrals cache: (', mo_integrals_cache_min, ', ', mo_integrals_cache_max, '), ', & + shiftr(mo_integrals_cache_size_8, 17), 'MiB' END_PROVIDER -BEGIN_PROVIDER [ double precision, mo_integrals_cache, (0_8:(1_8*mo_integrals_cache_size)**4) ] +BEGIN_PROVIDER [ double precision, mo_integrals_cache, (0_8:mo_integrals_cache_size_8**4_8) ] implicit none BEGIN_DOC ! Cache of MO integrals for fast access @@ -67,8 +81,7 @@ BEGIN_PROVIDER [ double precision, mo_integrals_cache, (0_8:(1_8*mo_integrals_ca do k=mo_integrals_cache_min,mo_integrals_cache_max ii = int(l-mo_integrals_cache_min,8) ii = ior( shiftl(ii,mo_integrals_cache_shift), int(k-mo_integrals_cache_min,8)) - ii = shiftl(ii,mo_integrals_cache_shift) - ii = shiftl(ii,mo_integrals_cache_shift) + ii = shiftl(ii,2*mo_integrals_cache_shift) call dgemm('T','N', mo_integrals_cache_max-mo_integrals_cache_min+1, & mo_integrals_cache_max-mo_integrals_cache_min+1, & cholesky_mo_num, 1.d0, & @@ -328,7 +341,7 @@ double precision function mo_two_e_integral(i,j,k,l) END_DOC integer, intent(in) :: i,j,k,l double precision :: get_two_e_integral - PROVIDE mo_two_e_integrals_in_map mo_integrals_cache + PROVIDE all_mo_integrals !DIR$ FORCEINLINE mo_two_e_integral = get_two_e_integral(i,j,k,l,mo_integrals_map) return From 348f6067914983d8919644d464a0ad18e7af6dc4 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 20 Jan 2025 18:03:19 +0100 Subject: [PATCH 04/10] Provide mo...in_map -> all_mo_integrals --- src/casscf_cipsi/bielec.irp.f | 28 +-- src/casscf_cipsi/hessian.irp.f | 2 +- src/casscf_cipsi/hessian_old.irp.f | 52 ++-- src/ccsd/ccsd_space_orb_sub.irp.f | 2 +- src/cipsi/stochastic_cipsi.irp.f | 6 +- src/cipsi_utils/pt2_stoch_routines.irp.f | 2 +- src/cipsi_utils/slave_cipsi.irp.f | 2 +- src/dressing/dress_slave.irp.f | 2 +- src/dressing/dress_stoch_routines.irp.f | 2 +- src/fci/fci.irp.f | 5 +- src/fci/pt2.irp.f | 16 +- .../debug_gradient_list_opt.irp.f | 22 +- src/mo_optimization/debug_gradient_opt.irp.f | 22 +- .../debug_hessian_list_opt.irp.f | 32 +-- src/mo_optimization/debug_hessian_opt.irp.f | 34 +-- .../routine_opt_mos.irp.f | 6 +- .../run_orb_opt_trust_v2.irp.f | 38 +-- src/mo_two_e_ints/map_integrals.irp.f | 2 +- src/mu_of_r/f_hf_utils.irp.f | 80 +++--- src/mu_of_r/f_psi_i_a_v_utils.irp.f | 144 +++++------ src/mu_of_r/mu_of_r_conditions.irp.f | 82 +++---- src/tools/fcidump.irp.f | 2 +- src/tools/fcidump_pyscf.irp.f | 2 +- src/tools/four_idx_transform.irp.f | 2 +- src/trexio/export_trexio_routines.irp.f | 2 +- src/two_body_rdm/act_2_rdm.irp.f | 76 +++--- src/two_rdm_routines/davidson_like_2rdm.irp.f | 228 +++++++++--------- .../davidson_like_trans_2rdm.irp.f | 228 +++++++++--------- src/utils_cc/mo_integrals_cc.irp.f | 26 +- 29 files changed, 574 insertions(+), 573 deletions(-) diff --git a/src/casscf_cipsi/bielec.irp.f b/src/casscf_cipsi/bielec.irp.f index a4901985..cd96d25a 100644 --- a/src/casscf_cipsi/bielec.irp.f +++ b/src/casscf_cipsi/bielec.irp.f @@ -1,7 +1,7 @@ BEGIN_PROVIDER [real*8, bielec_PQxx_array, (mo_num, mo_num,n_core_inact_act_orb,n_core_inact_act_orb)] BEGIN_DOC - ! WARNING !!! Old version !!! NOT USED ANYMORE IN THE PROGRAM !!! TOO BIG TO BE STORED ON LARGE SYSTEMS !!! - ! + ! WARNING !!! Old version !!! NOT USED ANYMORE IN THE PROGRAM !!! TOO BIG TO BE STORED ON LARGE SYSTEMS !!! + ! ! Replaced by the Cholesky-based function bielec_PQxx ! ! bielec_PQxx_array : integral (pq|xx) with p,q arbitrary, x core or active @@ -13,10 +13,10 @@ BEGIN_PROVIDER [real*8, bielec_PQxx_array, (mo_num, mo_num,n_core_inact_act_orb, print*,'' print*,'Providing bielec_PQxx_array, WARNING IT CAN BE A VERY BIG ARRAY WHEN MO_NUM IS LARGE !!!' print*,'' - + bielec_PQxx_array(:,:,:,:) = 0.d0 - PROVIDE mo_two_e_integrals_in_map - + PROVIDE all_mo_integrals + !$OMP PARALLEL DEFAULT(NONE) & !$OMP PRIVATE(i,ii,j,jj,i3,j3) & !$OMP SHARED(n_core_inact_orb,list_core_inact,mo_num,bielec_PQxx_array, & @@ -61,8 +61,8 @@ END_PROVIDER BEGIN_PROVIDER [real*8, bielec_PxxQ_array, (mo_num,n_core_inact_act_orb,n_core_inact_act_orb, mo_num)] BEGIN_DOC - ! WARNING !!! Old version !!! NOT USED ANYMORE IN THE PROGRAM !!! TOO BIG TO BE STORED ON LARGE SYSTEMS !!! - ! + ! WARNING !!! Old version !!! NOT USED ANYMORE IN THE PROGRAM !!! TOO BIG TO BE STORED ON LARGE SYSTEMS !!! + ! ! Replaced by the Cholesky-based function bielec_PxxQ ! ! bielec_PxxQ_array : integral (px|xq) with p,q arbitrary, x core or active @@ -72,11 +72,11 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ_array, (mo_num,n_core_inact_act_orb,n_core_i integer :: i,j,ii,jj,p,q,i3,j3,t3,v3 double precision, allocatable :: integrals_array(:,:) real*8 :: mo_two_e_integral - + print*,'' print*,'Providing bielec_PxxQ_array, WARNING IT CAN BE A VERY BIG ARRAY WHEN MO_NUM IS LARGE !!!' print*,'' - PROVIDE mo_two_e_integrals_in_map + PROVIDE all_mo_integrals bielec_PxxQ_array = 0.d0 !$OMP PARALLEL DEFAULT(NONE) & @@ -85,7 +85,7 @@ BEGIN_PROVIDER [real*8, bielec_PxxQ_array, (mo_num,n_core_inact_act_orb,n_core_i !$OMP n_act_orb,mo_integrals_map,list_act) allocate(integrals_array(mo_num,mo_num)) - + !$OMP DO do i=1,n_core_inact_orb ii=list_core_inact(i) @@ -143,17 +143,17 @@ BEGIN_PROVIDER [real*8, bielecCI, (n_act_orb,n_act_orb,n_act_orb, mo_num)] BEGIN_DOC ! bielecCI : integrals (tu|vp) with p arbitrary, tuv active ! index p runs over the whole basis, t,u,v only over the active orbitals - ! + ! ! This array can be stored anyway. Ex: 50 active orbitals, 1500 MOs ==> 8x50^3x1500 = 1.5 Gb END_DOC implicit none integer :: i,j,k,p,t,u,v double precision, external :: mo_two_e_integral - double precision :: wall0, wall1 + double precision :: wall0, wall1 call wall_time(wall0) print*,'Providing bielecCI' - PROVIDE mo_two_e_integrals_in_map - + PROVIDE all_mo_integrals + !$OMP PARALLEL DO DEFAULT(NONE) & !$OMP PRIVATE(i,j,k,p,t,u,v) & !$OMP SHARED(mo_num,n_act_orb,list_act,bielecCI) diff --git a/src/casscf_cipsi/hessian.irp.f b/src/casscf_cipsi/hessian.irp.f index 9a7a9031..1ee073d2 100644 --- a/src/casscf_cipsi/hessian.irp.f +++ b/src/casscf_cipsi/hessian.irp.f @@ -380,7 +380,7 @@ BEGIN_PROVIDER [double precision, hessmat, (nMonoEx,nMonoEx)] ! c-v | X X ! a-v | X - provide mo_two_e_integrals_in_map + provide all_mo_integrals hessmat = 0.d0 diff --git a/src/casscf_cipsi/hessian_old.irp.f b/src/casscf_cipsi/hessian_old.irp.f index d17f1f0a..0f7c8682 100644 --- a/src/casscf_cipsi/hessian_old.irp.f +++ b/src/casscf_cipsi/hessian_old.irp.f @@ -13,18 +13,18 @@ BEGIN_PROVIDER [real*8, hessmat_old, (nMonoEx,nMonoEx)] integer :: jndx,jhole,jpart character*3 :: iexc,jexc real*8 :: res - + if (bavard) then write(6,*) ' providing Hessian matrix hessmat_old ' write(6,*) ' nMonoEx = ',nMonoEx endif - + do indx=1,nMonoEx do jndx=1,nMonoEx hessmat_old(indx,jndx)=0.D0 end do end do - + do indx=1,nMonoEx ihole=excit(1,indx) ipart=excit(2,indx) @@ -38,7 +38,7 @@ BEGIN_PROVIDER [real*8, hessmat_old, (nMonoEx,nMonoEx)] hessmat_old(jndx,indx)=res end do end do - + END_PROVIDER subroutine calc_hess_elem(ihole,ipart,jhole,jpart,res) @@ -75,9 +75,9 @@ subroutine calc_hess_elem(ihole,ipart,jhole,jpart,res) integer :: nu_rs_possible integer :: mu_pqrs_possible integer :: mu_rspq_possible - + res=0.D0 - + ! the terms <0|E E H |0> do mu=1,n_det ! get the string of the determinant @@ -174,10 +174,10 @@ subroutine calc_hess_elem(ihole,ipart,jhole,jpart,res) end if end do end do - + ! state-averaged Hessian res*=1.D0/dble(N_states) - + end subroutine calc_hess_elem BEGIN_PROVIDER [real*8, hessmat_peter, (nMonoEx,nMonoEx)] @@ -190,26 +190,26 @@ BEGIN_PROVIDER [real*8, hessmat_peter, (nMonoEx,nMonoEx)] END_DOC implicit none integer :: i,j,t,u,a,b,indx,jndx,bstart,ustart,indx_shift - + real*8 :: hessmat_itju real*8 :: hessmat_itja real*8 :: hessmat_itua real*8 :: hessmat_iajb real*8 :: hessmat_iatb real*8 :: hessmat_taub - + if (bavard) then write(6,*) ' providing Hessian matrix hessmat_peter ' write(6,*) ' nMonoEx = ',nMonoEx endif - provide mo_two_e_integrals_in_map - + provide all_mo_integrals + !$OMP PARALLEL DEFAULT(NONE) & !$OMP SHARED(hessmat_peter,n_core_inact_orb,n_act_orb,n_virt_orb,nMonoEx) & !$OMP PRIVATE(i,indx,jndx,j,ustart,t,u,a,bstart,indx_shift) !$OMP DO - ! (DOUBLY OCCUPIED ---> ACT ) + ! (DOUBLY OCCUPIED ---> ACT ) do i=1,n_core_inact_orb do t=1,n_act_orb indx = t + (i-1)*n_act_orb @@ -226,14 +226,14 @@ BEGIN_PROVIDER [real*8, hessmat_peter, (nMonoEx,nMonoEx)] jndx+=1 end do end do - ! (DOUBLY OCCUPIED ---> VIRTUAL) + ! (DOUBLY OCCUPIED ---> VIRTUAL) do j=1,n_core_inact_orb do a=1,n_virt_orb hessmat_peter(jndx,indx)=hessmat_itja(i,t,j,a) jndx+=1 end do end do - ! (ACTIVE ---> VIRTUAL) + ! (ACTIVE ---> VIRTUAL) do u=1,n_act_orb do a=1,n_virt_orb hessmat_peter(jndx,indx)=hessmat_itua(i,t,u,a) @@ -243,15 +243,15 @@ BEGIN_PROVIDER [real*8, hessmat_peter, (nMonoEx,nMonoEx)] end do end do !$OMP END DO NOWAIT - + indx_shift = n_core_inact_orb*n_act_orb !$OMP DO - ! (DOUBLY OCCUPIED ---> VIRTUAL) + ! (DOUBLY OCCUPIED ---> VIRTUAL) do a=1,n_virt_orb do i=1,n_core_inact_orb indx = a + (i-1)*n_virt_orb + indx_shift jndx=indx - ! (DOUBLY OCCUPIED ---> VIRTUAL) + ! (DOUBLY OCCUPIED ---> VIRTUAL) do j=i,n_core_inact_orb if (i.eq.j) then bstart=a @@ -263,7 +263,7 @@ BEGIN_PROVIDER [real*8, hessmat_peter, (nMonoEx,nMonoEx)] jndx+=1 end do end do - ! (ACT ---> VIRTUAL) + ! (ACT ---> VIRTUAL) do t=1,n_act_orb do b=1,n_virt_orb hessmat_peter(jndx,indx)=hessmat_iatb(i,a,t,b) @@ -273,15 +273,15 @@ BEGIN_PROVIDER [real*8, hessmat_peter, (nMonoEx,nMonoEx)] end do end do !$OMP END DO NOWAIT - + indx_shift += n_core_inact_orb*n_virt_orb - !$OMP DO - ! (ACT ---> VIRTUAL) + !$OMP DO + ! (ACT ---> VIRTUAL) do a=1,n_virt_orb do t=1,n_act_orb indx = a + (t-1)*n_virt_orb + indx_shift jndx=indx - ! (ACT ---> VIRTUAL) + ! (ACT ---> VIRTUAL) do u=t,n_act_orb if (t.eq.u) then bstart=a @@ -295,16 +295,16 @@ BEGIN_PROVIDER [real*8, hessmat_peter, (nMonoEx,nMonoEx)] end do end do end do - !$OMP END DO + !$OMP END DO !$OMP END PARALLEL - + do jndx=1,nMonoEx do indx=1,jndx-1 hessmat_peter(indx,jndx) = hessmat_peter(jndx,indx) enddo enddo - + END_PROVIDER diff --git a/src/ccsd/ccsd_space_orb_sub.irp.f b/src/ccsd/ccsd_space_orb_sub.irp.f index e4907f22..30f134fc 100644 --- a/src/ccsd/ccsd_space_orb_sub.irp.f +++ b/src/ccsd/ccsd_space_orb_sub.irp.f @@ -35,7 +35,7 @@ subroutine run_ccsd_space_orb PROVIDE cholesky_mo_transp FREE cholesky_ao else - PROVIDE mo_two_e_integrals_in_map + PROVIDE all_mo_integrals endif det = psi_det(:,:,cc_ref) diff --git a/src/cipsi/stochastic_cipsi.irp.f b/src/cipsi/stochastic_cipsi.irp.f index 289040f0..40d77ef5 100644 --- a/src/cipsi/stochastic_cipsi.irp.f +++ b/src/cipsi/stochastic_cipsi.irp.f @@ -1,11 +1,11 @@ -subroutine run_stochastic_cipsi(Ev,PT2) +subroutine run_stochastic_cipsi(Ev,PT2) use selection_types implicit none BEGIN_DOC ! Selected Full Configuration Interaction with Stochastic selection and PT2. END_DOC integer :: i,j,k - double precision, intent(out) :: Ev(N_states), PT2(N_states) + double precision, intent(out) :: Ev(N_states), PT2(N_states) double precision, allocatable :: zeros(:) integer :: to_select type(pt2_type) :: pt2_data, pt2_data_err @@ -14,7 +14,7 @@ subroutine run_stochastic_cipsi(Ev,PT2) double precision :: rss double precision, external :: memory_of_double - PROVIDE H_apply_buffer_allocated distributed_davidson mo_two_e_integrals_in_map + PROVIDE H_apply_buffer_allocated distributed_davidson all_mo_integrals threshold_generators = 1.d0 SOFT_TOUCH threshold_generators diff --git a/src/cipsi_utils/pt2_stoch_routines.irp.f b/src/cipsi_utils/pt2_stoch_routines.irp.f index 744c4006..144d052d 100644 --- a/src/cipsi_utils/pt2_stoch_routines.irp.f +++ b/src/cipsi_utils/pt2_stoch_routines.irp.f @@ -165,7 +165,7 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in) state_average_weight(pt2_stoch_istate) = 1.d0 TOUCH state_average_weight pt2_stoch_istate selection_weight - PROVIDE nproc pt2_F mo_two_e_integrals_in_map mo_one_e_integrals pt2_w + PROVIDE nproc pt2_F all_mo_integrals mo_one_e_integrals pt2_w PROVIDE psi_selectors pt2_u pt2_J pt2_R call new_parallel_job(zmq_to_qp_run_socket, zmq_socket_pull, 'pt2') diff --git a/src/cipsi_utils/slave_cipsi.irp.f b/src/cipsi_utils/slave_cipsi.irp.f index 3e778270..ae8d3472 100644 --- a/src/cipsi_utils/slave_cipsi.irp.f +++ b/src/cipsi_utils/slave_cipsi.irp.f @@ -14,7 +14,7 @@ subroutine run_slave_cipsi end subroutine provide_everything - PROVIDE H_apply_buffer_allocated mo_two_e_integrals_in_map psi_det_generators psi_coef_generators psi_det_sorted_bit psi_selectors n_det_generators n_states generators_bitmask zmq_context N_states_diag + PROVIDE H_apply_buffer_allocated all_mo_integrals psi_det_generators psi_coef_generators psi_det_sorted_bit psi_selectors n_det_generators n_states generators_bitmask zmq_context N_states_diag PROVIDE pt2_e0_denominator mo_num N_int ci_energy mpi_master zmq_state zmq_context PROVIDE psi_det psi_coef threshold_generators state_average_weight PROVIDE N_det_selectors pt2_stoch_istate N_det selection_weight pseudo_sym diff --git a/src/dressing/dress_slave.irp.f b/src/dressing/dress_slave.irp.f index 04e4f01b..e0274a8a 100644 --- a/src/dressing/dress_slave.irp.f +++ b/src/dressing/dress_slave.irp.f @@ -15,7 +15,7 @@ subroutine dress_slave end subroutine provide_everything - PROVIDE H_apply_buffer_allocated mo_two_e_integrals_in_map psi_det_generators psi_coef_generators psi_det_sorted_bit psi_selectors n_det_generators n_states generators_bitmask zmq_context + PROVIDE H_apply_buffer_allocated all_mo_integrals psi_det_generators psi_coef_generators psi_det_sorted_bit psi_selectors n_det_generators n_states generators_bitmask zmq_context end subroutine run_wf diff --git a/src/dressing/dress_stoch_routines.irp.f b/src/dressing/dress_stoch_routines.irp.f index 6b37fad1..cdbe2c4f 100644 --- a/src/dressing/dress_stoch_routines.irp.f +++ b/src/dressing/dress_stoch_routines.irp.f @@ -258,7 +258,7 @@ subroutine ZMQ_dress(E, dress, delta_out, delta_s2_out, relative_error) state_average_weight(dress_stoch_istate) = 1.d0 TOUCH state_average_weight dress_stoch_istate - provide nproc mo_two_e_integrals_in_map mo_one_e_integrals psi_selectors pt2_F pt2_N_teeth dress_M_m + provide nproc all_mo_integrals mo_one_e_integrals psi_selectors pt2_F pt2_N_teeth dress_M_m print *, '========== ================= ================= =================' print *, ' Samples Energy Stat. Error Seconds ' diff --git a/src/fci/fci.irp.f b/src/fci/fci.irp.f index 2059a53b..d3ecc45d 100644 --- a/src/fci/fci.irp.f +++ b/src/fci/fci.irp.f @@ -36,8 +36,9 @@ program fci ! END_DOC + PROVIDE all_mo_integrals if (.not.is_zmq_slave) then - PROVIDE psi_det psi_coef mo_two_e_integrals_in_map + PROVIDE psi_det psi_coef write(json_unit,json_array_open_fmt) 'fci' @@ -55,7 +56,7 @@ program fci call json_close else - PROVIDE mo_two_e_integrals_in_map pt2_min_parallel_tasks + PROVIDE pt2_min_parallel_tasks call run_slave_cipsi diff --git a/src/fci/pt2.irp.f b/src/fci/pt2.irp.f index 1c9f9dcd..53bf1699 100644 --- a/src/fci/pt2.irp.f +++ b/src/fci/pt2.irp.f @@ -11,7 +11,7 @@ program pt2 ! ! The main option for the |PT2| correction is the ! :option:`perturbation pt2_relative_error` which is the relative - ! stochastic error on the |PT2| to reach before stopping the + ! stochastic error on the |PT2| to reach before stopping the ! sampling. ! END_DOC @@ -19,7 +19,7 @@ program pt2 read_wf = .True. threshold_generators = 1.d0 SOFT_TOUCH read_wf threshold_generators - PROVIDE mo_two_e_integrals_in_map + PROVIDE all_mo_integrals PROVIDE psi_energy call run else @@ -32,22 +32,22 @@ subroutine run use selection_types integer :: i,j,k logical, external :: detEq - + type(pt2_type) :: pt2_data, pt2_data_err integer :: degree integer :: n_det_before, to_select double precision :: threshold_davidson_in - + double precision :: relative_error double precision, allocatable :: E_CI_before(:) - + allocate ( E_CI_before(N_states)) call pt2_alloc(pt2_data, N_states) call pt2_alloc(pt2_data_err, N_states) - + E_CI_before(:) = psi_energy(:) + nuclear_repulsion relative_error=PT2_relative_error - + if (do_pt2) then call ZMQ_pt2(psi_energy_with_nucl_rep, pt2_data, pt2_data_err, relative_error, 0) ! Stochastic PT2 else @@ -56,7 +56,7 @@ subroutine run call print_summary(psi_energy_with_nucl_rep(1:N_states), & pt2_data, pt2_data_err, N_det,N_configuration,N_states,psi_s2) - + call save_energy(E_CI_before, pt2_data % pt2) call pt2_dealloc(pt2_data) call pt2_dealloc(pt2_data_err) diff --git a/src/mo_optimization/debug_gradient_list_opt.irp.f b/src/mo_optimization/debug_gradient_list_opt.irp.f index 32cea90c..dcc731fc 100644 --- a/src/mo_optimization/debug_gradient_list_opt.irp.f +++ b/src/mo_optimization/debug_gradient_list_opt.irp.f @@ -19,7 +19,7 @@ program debug_gradient_list - + implicit none ! Variables @@ -30,12 +30,12 @@ program debug_gradient_list double precision :: threshold double precision :: max_error, max_elem, norm integer :: nb_error - + m = dim_list_act_orb - ! Definition of n + ! Definition of n n = m*(m-1)/2 - PROVIDE mo_two_e_integrals_in_map ! Verifier pour suppression + PROVIDE all_mo_integrals ! Verifier pour suppression ! Allocation allocate(v_grad(n), v_grad2(n)) @@ -44,15 +44,15 @@ program debug_gradient_list call diagonalize_ci ! Verifier pour suppression - ! Gradient + ! Gradient call gradient_list_opt(n,m,list_act,v_grad,max_elem,norm) call first_gradient_list_opt(n,m,list_act,v_grad2) - - + + v_grad = v_grad - v_grad2 nb_error = 0 - max_error = 0d0 - threshold = 1d-12 + max_error = 0d0 + threshold = 1d-12 do i = 1, n if (ABS(v_grad(i)) > threshold) then @@ -65,9 +65,9 @@ program debug_gradient_list endif enddo - + print*,'' - print*,'Check the gradient' + print*,'Check the gradient' print*,'Threshold:', threshold print*,'Nb error:', nb_error print*,'Max error:', max_error diff --git a/src/mo_optimization/debug_gradient_opt.irp.f b/src/mo_optimization/debug_gradient_opt.irp.f index 529a02b6..7aba985b 100644 --- a/src/mo_optimization/debug_gradient_opt.irp.f +++ b/src/mo_optimization/debug_gradient_opt.irp.f @@ -19,7 +19,7 @@ program debug_gradient - + implicit none ! Variables @@ -30,27 +30,27 @@ program debug_gradient double precision :: threshold double precision :: max_error, max_elem integer :: nb_error - - ! Definition of n + + ! Definition of n n = mo_num*(mo_num-1)/2 - PROVIDE mo_two_e_integrals_in_map ! Check for suppression + PROVIDE all_mo_integrals ! Check for suppression ! Allocation allocate(v_grad(n), v_grad2(n)) ! Calculation - call diagonalize_ci + call diagonalize_ci - ! Gradient + ! Gradient call first_gradient_opt(n,v_grad) call gradient_opt(n,v_grad2,max_elem) - + v_grad = v_grad - v_grad2 nb_error = 0 - max_error = 0d0 - threshold = 1d-12 + max_error = 0d0 + threshold = 1d-12 do i = 1, n if (ABS(v_grad(i)) > threshold) then @@ -63,9 +63,9 @@ program debug_gradient endif enddo - + print*,'' - print*,'Check the gradient' + print*,'Check the gradient' print*,'Threshold :', threshold print*,'Nb error :', nb_error print*,'Max error :', max_error diff --git a/src/mo_optimization/debug_hessian_list_opt.irp.f b/src/mo_optimization/debug_hessian_list_opt.irp.f index 65a7bcf3..417642b6 100644 --- a/src/mo_optimization/debug_hessian_list_opt.irp.f +++ b/src/mo_optimization/debug_hessian_list_opt.irp.f @@ -43,18 +43,18 @@ program debug_hessian_list_opt double precision :: max_error, max_error_H integer :: nb_error, nb_error_H double precision :: threshold - + m = dim_list_act_orb !mo_num - ! Definition of n + ! Definition of n n = m*(m-1)/2 - PROVIDE mo_two_e_integrals_in_map + PROVIDE all_mo_integrals - ! Hessian + ! Hessian if (optimization_method == 'full') then print*,'Use the full hessian matrix' - allocate(H(n,n),H2(n,n)) + allocate(H(n,n),H2(n,n)) allocate(h_f(m,m,m,m),h_f2(m,m,m,m)) call hessian_list_opt(n,m,list_act,H,h_f) @@ -65,7 +65,7 @@ program debug_hessian_list_opt h_f = h_f - h_f2 H = H - H2 max_error = 0d0 - nb_error = 0 + nb_error = 0 threshold = 1d-12 do l = 1, m @@ -99,7 +99,7 @@ program debug_hessian_list_opt endif enddo - enddo + enddo ! Deallocation deallocate(H, H2, h_f, h_f2) @@ -110,26 +110,26 @@ program debug_hessian_list_opt allocate(H(n,1),H2(n,1)) call diag_hessian_list_opt(n,m,list_act,H) call first_diag_hessian_list_opt(n,m,list_act,H2) - + H = H - H2 - + max_error_H = 0d0 nb_error_H = 0 - + do i = 1, n if (ABS(H(i,1)) > threshold) then print*, H(i,1) nb_error_H = nb_error_H + 1 - + if (ABS(H(i,1)) > ABS(max_error_H)) then max_error_H = H(i,1) endif - + endif enddo - + endif - + print*,'' if (optimization_method == 'full') then print*,'Check of the full hessian' @@ -140,8 +140,8 @@ program debug_hessian_list_opt else print*,'Check of the diagonal hessian' endif - + print*,'Nb error_H:', nb_error_H print*,'Max error_H:', max_error_H - + end program diff --git a/src/mo_optimization/debug_hessian_opt.irp.f b/src/mo_optimization/debug_hessian_opt.irp.f index 684a0da5..729fc2a3 100644 --- a/src/mo_optimization/debug_hessian_opt.irp.f +++ b/src/mo_optimization/debug_hessian_opt.irp.f @@ -36,20 +36,20 @@ program debug_hessian double precision :: max_error, max_error_H integer :: nb_error, nb_error_H double precision :: threshold - - ! Definition of n + + ! Definition of n n = mo_num*(mo_num-1)/2 - PROVIDE mo_two_e_integrals_in_map + PROVIDE all_mo_integrals ! Allocation - allocate(H(n,n),H2(n,n)) + allocate(H(n,n),H2(n,n)) allocate(h_f(mo_num,mo_num,mo_num,mo_num),h_f2(mo_num,mo_num,mo_num,mo_num)) ! Calculation - - ! Hessian - if (optimization_method == 'full') then + + ! Hessian + if (optimization_method == 'full') then print*,'Use the full hessian matrix' call hessian_opt(n,H,h_f) @@ -59,7 +59,7 @@ program debug_hessian h_f = h_f - h_f2 H = H - H2 max_error = 0d0 - nb_error = 0 + nb_error = 0 threshold = 1d-12 do l = 1, mo_num @@ -93,7 +93,7 @@ program debug_hessian endif enddo - enddo + enddo elseif (optimization_method == 'diag') then @@ -128,43 +128,43 @@ program debug_hessian enddo h=H-H2 - + max_error_H = 0d0 nb_error_H = 0 - + do j = 1, n do i = 1, n if (ABS(H(i,j)) > threshold) then print*, H(i,j) nb_error_H = nb_error_H + 1 - + if (ABS(H(i,j)) > ABS(max_error_H)) then max_error_H = H(i,j) endif - + endif enddo enddo - + else print*,'Unknown optimization_method, please select full, diag' call abort endif - + print*,'' if (optimization_method == 'full') then print*,'Check the full hessian' else print*,'Check the diagonal hessian' endif - + print*,'Threshold :', threshold print*,'Nb error :', nb_error print*,'Max error :', max_error print*,'' print*,'Nb error_H :', nb_error_H print*,'Max error_H :', max_error_H - + ! Deallocation deallocate(H,H2,h_f,h_f2) diff --git a/src/mo_optimization_utils/routine_opt_mos.irp.f b/src/mo_optimization_utils/routine_opt_mos.irp.f index fceba2c5..742e074b 100644 --- a/src/mo_optimization_utils/routine_opt_mos.irp.f +++ b/src/mo_optimization_utils/routine_opt_mos.irp.f @@ -9,7 +9,7 @@ subroutine run_optimization_mos_CIPSI logical :: not_converged character (len=100) :: filename - PROVIDE psi_det psi_coef mo_two_e_integrals_in_map ao_pseudo_integrals + PROVIDE psi_det psi_coef all_mo_integrals ao_pseudo_integrals allocate(Ev(N_states),PT2(N_states)) not_converged = .True. @@ -30,7 +30,7 @@ subroutine run_optimization_mos_CIPSI print*,'======================' print*,' Cipsi step:', nb_iter print*,'======================' - print*,'' + print*,'' print*,'********** cipsi step **********' ! cispi calculation call run_stochastic_cipsi(Ev,PT2) @@ -70,7 +70,7 @@ subroutine run_optimization_mos_CIPSI print*, 'The program will exit' exit endif - + ! To double the number of determinants in the wf N_det_max = int(dble(n_det * 2)*0.9) TOUCH N_det_max diff --git a/src/mo_optimization_utils/run_orb_opt_trust_v2.irp.f b/src/mo_optimization_utils/run_orb_opt_trust_v2.irp.f index e1431255..afa4da0a 100644 --- a/src/mo_optimization_utils/run_orb_opt_trust_v2.irp.f +++ b/src/mo_optimization_utils/run_orb_opt_trust_v2.irp.f @@ -93,7 +93,7 @@ subroutine run_orb_opt_trust_v2 integer,allocatable :: tmp_list(:), key(:) double precision, allocatable :: tmp_m_x(:,:),tmp_R(:,:), tmp_x(:), W(:,:), e_val(:) - PROVIDE mo_two_e_integrals_in_map ci_energy psi_det psi_coef + PROVIDE all_mo_integrals ci_energy psi_det psi_coef ! Allocation @@ -110,17 +110,17 @@ allocate(tmp_R(m,m), tmp_m_x(m,m), tmp_x(tmp_n)) allocate(e_val(tmp_n),key(tmp_n),v_grad(tmp_n)) ! Method -! There are three different methods : +! There are three different methods : ! - the "full" hessian, which uses all the elements of the hessian ! matrix" ! - the "diagonal" hessian, which uses only the diagonal elements of the ! hessian -! - without the hessian (hessian = identity matrix) +! - without the hessian (hessian = identity matrix) !Display the method print*, 'Method :', optimization_method -if (optimization_method == 'full') then +if (optimization_method == 'full') then print*, 'Full hessian' allocate(H(tmp_n,tmp_n), h_f(m,m,m,m),W(tmp_n,tmp_n)) tmp_n2 = tmp_n @@ -147,13 +147,13 @@ print*, 'Absolute value of the hessian:', absolute_eig ! - We diagonalize the hessian ! - We compute a step and loop to reduce the radius of the ! trust region (and the size of the step by the way) until the step is -! accepted -! - We repeat the process until the convergence +! accepted +! - We repeat the process until the convergence ! NB: the convergence criterion can be changed ! Loop until the convergence of the optimization -! call diagonalize_ci +! call diagonalize_ci !### Initialization ### nb_iter = 0 @@ -183,14 +183,14 @@ do while (not_converged) ! Full hessian call hessian_list_opt(tmp_n, m, tmp_list, H, h_f) - ! Diagonalization of the hessian + ! Diagonalization of the hessian call diagonalization_hessian(tmp_n, H, e_val, w) elseif (optimization_method == 'diag') then - ! Diagonal hessian + ! Diagonal hessian call diag_hessian_list_opt(tmp_n, m, tmp_list, H) else - ! Identity matrix + ! Identity matrix do tmp_i = 1, tmp_n H(tmp_i,1) = 1d0 enddo @@ -212,7 +212,7 @@ do while (not_converged) endif ! Init before the internal loop - cancel_step = .True. ! To enter in the loop just after + cancel_step = .True. ! To enter in the loop just after nb_cancel = 0 nb_sub_iter = 0 @@ -225,15 +225,15 @@ do while (not_converged) print*,'Max elem grad:', max_elem_grad print*,'-----------------------------' - ! Hessian,gradient,Criterion -> x - call trust_region_step_w_expected_e(tmp_n,tmp_n2,H,W,e_val,v_grad,prev_criterion,rho,nb_iter,delta,criterion_model,tmp_x,must_exit) + ! Hessian,gradient,Criterion -> x + call trust_region_step_w_expected_e(tmp_n,tmp_n2,H,W,e_val,v_grad,prev_criterion,rho,nb_iter,delta,criterion_model,tmp_x,must_exit) if (must_exit) then print*,'step_in_trust_region sends: Exit' exit endif - ! 1D tmp -> 2D tmp + ! 1D tmp -> 2D tmp call vec_to_mat_v2(tmp_n, m, tmp_x, tmp_m_x) ! Rotation matrix for the active MOs @@ -250,14 +250,14 @@ do while (not_converged) call sub_to_full_rotation_matrix(m, tmp_list, tmp_R, R) ! MO rotations - call apply_mo_rotation(R, prev_mos) + call apply_mo_rotation(R, prev_mos) ! Update of the energy before the diagonalization of the hamiltonian call clear_mo_map - TOUCH mo_coef psi_det psi_coef ci_energy two_e_dm_mo + TOUCH mo_coef psi_det psi_coef ci_energy two_e_dm_mo call state_average_energy(criterion) - ! Criterion -> step accepted or rejected + ! Criterion -> step accepted or rejected call trust_region_is_step_cancelled(nb_iter, prev_criterion, criterion, criterion_model, rho, cancel_step) ! Cancellation of the step if necessary @@ -283,7 +283,7 @@ do while (not_converged) ! To exit the external loop if must_exit = .True. if (must_exit) then exit - endif + endif ! Step accepted, nb iteration + 1 nb_iter = nb_iter + 1 @@ -295,7 +295,7 @@ do while (not_converged) endif if (nb_iter >= optimization_max_nb_iter) then print*,'Not converged: nb_iter >= optimization_max_nb_iter' - not_converged = .False. + not_converged = .False. endif if (.not. not_converged) then diff --git a/src/mo_two_e_ints/map_integrals.irp.f b/src/mo_two_e_ints/map_integrals.irp.f index e6b2967a..06300666 100644 --- a/src/mo_two_e_ints/map_integrals.irp.f +++ b/src/mo_two_e_ints/map_integrals.irp.f @@ -6,7 +6,7 @@ BEGIN_PROVIDER [ logical, all_mo_integrals ] ! Used to provide everything needed before using MO integrals ! PROVIDE all_mo_integrals END_DOC - PROVIDE mo_two_e_integrals_in_map mo_integrals_cache mo_two_e_integrals_jj_exchange mo_two_e_integrals_jj_anti mo_two_e_integrals_jj big_array_exchange_integrals big_array_coulomb_integrals + PROVIDE mo_two_e_integrals_in_map mo_integrals_cache mo_two_e_integrals_jj_exchange mo_two_e_integrals_jj_anti mo_two_e_integrals_jj big_array_exchange_integrals big_array_coulomb_integrals mo_one_e_integrals END_PROVIDER diff --git a/src/mu_of_r/f_hf_utils.irp.f b/src/mu_of_r/f_hf_utils.irp.f index 102b40a0..7c40fc80 100644 --- a/src/mu_of_r/f_hf_utils.irp.f +++ b/src/mu_of_r/f_hf_utils.irp.f @@ -1,41 +1,41 @@ BEGIN_PROVIDER [double precision, two_e_int_hf_f, (n_basis_orb,n_basis_orb,n_max_occ_val_orb_for_hf,n_max_occ_val_orb_for_hf)] implicit none BEGIN_DOC -! list of two-electron integrals (built with the MOs belonging to the \mathcal{B} space) -! -! needed to compute the function f_{HF}(r_1,r_2) +! list of two-electron integrals (built with the MOs belonging to the \mathcal{B} space) +! +! needed to compute the function f_{HF}(r_1,r_2) ! ! two_e_int_hf_f(j,i,n,m) = < j i | n m > where all orbitals belong to "list_basis" END_DOC integer :: orb_i,orb_j,i,j,orb_m,orb_n,m,n double precision :: get_two_e_integral - PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals - do orb_m = 1, n_max_occ_val_orb_for_hf! electron 1 + PROVIDE all_mo_integrals big_array_exchange_integrals + do orb_m = 1, n_max_occ_val_orb_for_hf! electron 1 m = list_valence_orb_for_hf(orb_m,1) - do orb_n = 1, n_max_occ_val_orb_for_hf! electron 2 + do orb_n = 1, n_max_occ_val_orb_for_hf! electron 2 n = list_valence_orb_for_hf(orb_n,1) - do orb_i = 1, n_basis_orb ! electron 1 + do orb_i = 1, n_basis_orb ! electron 1 i = list_basis(orb_i) - do orb_j = 1, n_basis_orb ! electron 2 + do orb_j = 1, n_basis_orb ! electron 2 j = list_basis(orb_j) ! 2 1 2 1 - two_e_int_hf_f(orb_j,orb_i,orb_n,orb_m) = get_two_e_integral(m,n,i,j,mo_integrals_map) + two_e_int_hf_f(orb_j,orb_i,orb_n,orb_m) = get_two_e_integral(m,n,i,j,mo_integrals_map) enddo enddo enddo enddo -END_PROVIDER +END_PROVIDER subroutine f_HF_valence_ab(r1,r2,f_HF_val_ab,two_bod_dens) implicit none BEGIN_DOC -! f_HF_val_ab(r1,r2) = function f_{\Psi^B}(X_1,X_2) of Eq. (22) of J. Chem. Phys. 149, 194301 (2018) -! -! for alpha beta spins and an HF wave function and excluding the "core" orbitals (see Eq. 16a of Phys.Chem.Lett.2019, 10, 2931 2937) -! -! two_bod_dens = on-top pair density of the HF wave function +! f_HF_val_ab(r1,r2) = function f_{\Psi^B}(X_1,X_2) of Eq. (22) of J. Chem. Phys. 149, 194301 (2018) ! -! < HF | wee_{\alpha\beta} | HF > = \int (r1,r2) f_HF_ab(r1,r2) excluding all contributions from "core" "electrons" +! for alpha beta spins and an HF wave function and excluding the "core" orbitals (see Eq. 16a of Phys.Chem.Lett.2019, 10, 2931 2937) +! +! two_bod_dens = on-top pair density of the HF wave function +! +! < HF | wee_{\alpha\beta} | HF > = \int (r1,r2) f_HF_ab(r1,r2) excluding all contributions from "core" "electrons" END_DOC double precision, intent(in) :: r1(3), r2(3) double precision, intent(out):: f_HF_val_ab,two_bod_dens @@ -50,8 +50,8 @@ subroutine f_HF_valence_ab(r1,r2,f_HF_val_ab,two_bod_dens) allocate(mos_array_valence_r1(n_basis_orb) , mos_array_valence_r2(n_basis_orb), mos_array_r1(mo_num), mos_array_r2(mo_num)) allocate(mos_array_valence_hf_r1(n_occ_val_orb_for_hf(1)) , mos_array_valence_hf_r2(n_occ_val_orb_for_hf(2)) ) ! You get all orbitals in r_1 and r_2 - call give_all_mos_at_r(r1,mos_array_r1) - call give_all_mos_at_r(r2,mos_array_r2) + call give_all_mos_at_r(r1,mos_array_r1) + call give_all_mos_at_r(r2,mos_array_r2) ! You extract the occupied ALPHA/BETA orbitals belonging to the space \mathcal{A} do i_m = 1, n_occ_val_orb_for_hf(1) mos_array_valence_hf_r1(i_m) = mos_array_r1(list_valence_orb_for_hf(i_m,1)) @@ -61,7 +61,7 @@ subroutine f_HF_valence_ab(r1,r2,f_HF_val_ab,two_bod_dens) enddo ! You extract the orbitals belonging to the space \mathcal{B} - do i_m = 1, n_basis_orb + do i_m = 1, n_basis_orb mos_array_valence_r1(i_m) = mos_array_r1(list_basis(i_m)) mos_array_valence_r2(i_m) = mos_array_r2(list_basis(i_m)) enddo @@ -70,24 +70,24 @@ subroutine f_HF_valence_ab(r1,r2,f_HF_val_ab,two_bod_dens) f_HF_val_ab = 0.d0 two_bod_dens = 0.d0 ! You browse all OCCUPIED ALPHA electrons in the \mathcal{A} space - do m = 1, n_occ_val_orb_for_hf(1)! electron 1 + do m = 1, n_occ_val_orb_for_hf(1)! electron 1 ! You browse all OCCUPIED BETA electrons in the \mathcal{A} space - do n = 1, n_occ_val_orb_for_hf(2)! electron 2 + do n = 1, n_occ_val_orb_for_hf(2)! electron 2 ! two_bod_dens(r_1,r_2) = n_alpha(r_1) * n_beta(r_2) - two_bod_dens += mos_array_valence_hf_r1(m) * mos_array_valence_hf_r1(m) * mos_array_valence_hf_r2(n) * mos_array_valence_hf_r2(n) - ! You browse all COUPLE OF ORBITALS in the \mathacal{B} space + two_bod_dens += mos_array_valence_hf_r1(m) * mos_array_valence_hf_r1(m) * mos_array_valence_hf_r2(n) * mos_array_valence_hf_r2(n) + ! You browse all COUPLE OF ORBITALS in the \mathacal{B} space do i = 1, n_basis_orb do j = 1, n_basis_orb ! 2 1 2 1 - f_HF_val_ab += two_e_int_hf_f(j,i,n,m) & - * mos_array_valence_r1(i) * mos_array_valence_hf_r1(m) & - * mos_array_valence_r2(j) * mos_array_valence_hf_r2(n) + f_HF_val_ab += two_e_int_hf_f(j,i,n,m) & + * mos_array_valence_r1(i) * mos_array_valence_hf_r1(m) & + * mos_array_valence_r2(j) * mos_array_valence_hf_r2(n) enddo enddo enddo enddo ! multiply by two to adapt to the N(N-1) normalization condition of the active two-rdm - f_HF_val_ab *= 2.d0 + f_HF_val_ab *= 2.d0 two_bod_dens *= 2.d0 end @@ -95,14 +95,14 @@ end subroutine integral_f_HF_valence_ab(r1,int_f_HF_val_ab) implicit none BEGIN_DOC -! in_f_HF_val_ab(r_1) = \int dr_2 f_{\Psi^B}(r_1,r_2) -! -! where f_{\Psi^B}(r_1,r_2) is defined by Eq. (22) of J. Chem. Phys. 149, 194301 (2018) +! in_f_HF_val_ab(r_1) = \int dr_2 f_{\Psi^B}(r_1,r_2) +! +! where f_{\Psi^B}(r_1,r_2) is defined by Eq. (22) of J. Chem. Phys. 149, 194301 (2018) ! ! for alpha beta spins and an HF wave function and excluding the "core" orbitals (see Eq. 16a of Phys.Chem.Lett.2019, 10, 2931 2937) ! -! Such function can be used to test if the f_HF_val_ab(r_1,r_2) is correctly built. -! +! Such function can be used to test if the f_HF_val_ab(r_1,r_2) is correctly built. +! ! < HF | wee_{\alpha\beta} | HF > = \int (r1) int_f_HF_val_ab(r_1) END_DOC double precision, intent(in) :: r1(3) @@ -114,28 +114,28 @@ subroutine integral_f_HF_valence_ab(r1,int_f_HF_val_ab) double precision, allocatable :: mos_array_valence_r1(:) double precision, allocatable :: mos_array_valence_hf_r1(:) double precision :: get_two_e_integral - call give_all_mos_at_r(r1,mos_array_r1) + call give_all_mos_at_r(r1,mos_array_r1) allocate(mos_array_valence_r1( n_basis_orb )) allocate(mos_array_valence_hf_r1( n_occ_val_orb_for_hf(1) ) ) do i_m = 1, n_occ_val_orb_for_hf(1) mos_array_valence_hf_r1(i_m) = mos_array_r1(list_valence_orb_for_hf(i_m,1)) enddo - do i_m = 1, n_basis_orb + do i_m = 1, n_basis_orb mos_array_valence_r1(i_m) = mos_array_r1(list_basis(i_m)) enddo int_f_HF_val_ab = 0.d0 ! You browse all OCCUPIED ALPHA electrons in the \mathcal{A} space - do m = 1, n_occ_val_orb_for_hf(1)! electron 1 + do m = 1, n_occ_val_orb_for_hf(1)! electron 1 ! You browse all OCCUPIED BETA electrons in the \mathcal{A} space - do n = 1, n_occ_val_orb_for_hf(2)! electron 2 - ! You browse all ORBITALS in the \mathacal{B} space + do n = 1, n_occ_val_orb_for_hf(2)! electron 2 + ! You browse all ORBITALS in the \mathacal{B} space do i = 1, n_basis_orb - ! due to integration in real-space and the use of orthonormal MOs, a Kronecker delta_jn shoes up + ! due to integration in real-space and the use of orthonormal MOs, a Kronecker delta_jn shoes up j = n ! 2 1 2 1 - int_f_HF_val_ab += two_e_int_hf_f(j,i,n,m) & - * mos_array_valence_r1(i) * mos_array_valence_hf_r1(m) + int_f_HF_val_ab += two_e_int_hf_f(j,i,n,m) & + * mos_array_valence_r1(i) * mos_array_valence_hf_r1(m) enddo enddo enddo diff --git a/src/mu_of_r/f_psi_i_a_v_utils.irp.f b/src/mu_of_r/f_psi_i_a_v_utils.irp.f index 69fa16ff..0d51155c 100644 --- a/src/mu_of_r/f_psi_i_a_v_utils.irp.f +++ b/src/mu_of_r/f_psi_i_a_v_utils.irp.f @@ -1,7 +1,7 @@ subroutine give_f_ii_val_ab(r1,r2,f_ii_val_ab,two_bod_dens) implicit none BEGIN_DOC -! contribution from purely inactive orbitals to f_{\Psi^B}(r_1,r_2) for a CAS wave function +! contribution from purely inactive orbitals to f_{\Psi^B}(r_1,r_2) for a CAS wave function END_DOC double precision, intent(in) :: r1(3),r2(3) double precision, intent(out):: f_ii_val_ab,two_bod_dens @@ -9,13 +9,13 @@ subroutine give_f_ii_val_ab(r1,r2,f_ii_val_ab,two_bod_dens) integer :: i_i,i_j double precision, allocatable :: mos_array_inact_r1(:),mos_array_inact_r2(:) double precision, allocatable :: mos_array_basis_r1(:),mos_array_basis_r2(:) - double precision, allocatable :: mos_array_r1(:) , mos_array_r2(:) + double precision, allocatable :: mos_array_r1(:) , mos_array_r2(:) double precision :: get_two_e_integral ! You get all orbitals in r_1 and r_2 allocate(mos_array_r1(mo_num) , mos_array_r2(mo_num) ) - call give_all_mos_at_r(r1,mos_array_r1) - call give_all_mos_at_r(r2,mos_array_r2) - ! You extract the inactive orbitals + call give_all_mos_at_r(r1,mos_array_r1) + call give_all_mos_at_r(r2,mos_array_r2) + ! You extract the inactive orbitals allocate(mos_array_inact_r1(n_inact_orb) , mos_array_inact_r2(n_inact_orb) ) do i_m = 1, n_inact_orb mos_array_inact_r1(i_m) = mos_array_r1(list_inact(i_m)) @@ -26,7 +26,7 @@ subroutine give_f_ii_val_ab(r1,r2,f_ii_val_ab,two_bod_dens) ! You extract the orbitals belonging to the space \mathcal{B} allocate(mos_array_basis_r1(n_basis_orb) , mos_array_basis_r2(n_basis_orb) ) - do i_m = 1, n_basis_orb + do i_m = 1, n_basis_orb mos_array_basis_r1(i_m) = mos_array_r1(list_basis(i_m)) mos_array_basis_r2(i_m) = mos_array_r2(list_basis(i_m)) enddo @@ -34,30 +34,30 @@ subroutine give_f_ii_val_ab(r1,r2,f_ii_val_ab,two_bod_dens) f_ii_val_ab = 0.d0 two_bod_dens = 0.d0 ! You browse all OCCUPIED ALPHA electrons in the \mathcal{A} space - do m = 1, n_inact_orb ! electron 1 + do m = 1, n_inact_orb ! electron 1 ! You browse all OCCUPIED BETA electrons in the \mathcal{A} space - do n = 1, n_inact_orb ! electron 2 + do n = 1, n_inact_orb ! electron 2 ! two_bod_dens(r_1,r_2) = n_alpha(r_1) * n_beta(r_2) - two_bod_dens += mos_array_inact_r1(m) * mos_array_inact_r1(m) * mos_array_inact_r2(n) * mos_array_inact_r2(n) - ! You browse all COUPLE OF ORBITALS in the \mathacal{B} space + two_bod_dens += mos_array_inact_r1(m) * mos_array_inact_r1(m) * mos_array_inact_r2(n) * mos_array_inact_r2(n) + ! You browse all COUPLE OF ORBITALS in the \mathacal{B} space do i = 1, n_basis_orb do j = 1, n_basis_orb ! 2 1 2 1 - f_ii_val_ab += two_e_int_ii_f(j,i,n,m) * mos_array_inact_r1(m) * mos_array_basis_r1(i) & - * mos_array_inact_r2(n) * mos_array_basis_r2(j) + f_ii_val_ab += two_e_int_ii_f(j,i,n,m) * mos_array_inact_r1(m) * mos_array_basis_r1(i) & + * mos_array_inact_r2(n) * mos_array_basis_r2(j) enddo enddo enddo enddo ! multiply by two to adapt to the N(N-1) normalization condition of the active two-rdm - f_ii_val_ab *= 2.d0 + f_ii_val_ab *= 2.d0 two_bod_dens *= 2.d0 end subroutine give_f_ia_val_ab(r1,r2,f_ia_val_ab,two_bod_dens,istate) BEGIN_DOC -! contribution from inactive and active orbitals to f_{\Psi^B}(r_1,r_2) for the "istate" state of a CAS wave function +! contribution from inactive and active orbitals to f_{\Psi^B}(r_1,r_2) for the "istate" state of a CAS wave function END_DOC implicit none integer, intent(in) :: istate @@ -65,7 +65,7 @@ subroutine give_f_ia_val_ab(r1,r2,f_ia_val_ab,two_bod_dens,istate) double precision, intent(out):: f_ia_val_ab,two_bod_dens integer :: i,orb_i,a,orb_a,n,m,b double precision :: rho - double precision, allocatable :: mos_array_r1(:) , mos_array_r2(:) + double precision, allocatable :: mos_array_r1(:) , mos_array_r2(:) double precision, allocatable :: mos_array_inact_r1(:),mos_array_inact_r2(:) double precision, allocatable :: mos_array_basis_r1(:),mos_array_basis_r2(:) double precision, allocatable :: mos_array_act_r1(:),mos_array_act_r2(:) @@ -75,10 +75,10 @@ subroutine give_f_ia_val_ab(r1,r2,f_ia_val_ab,two_bod_dens,istate) two_bod_dens = 0.d0 ! You get all orbitals in r_1 and r_2 allocate(mos_array_r1(mo_num) , mos_array_r2(mo_num) ) - call give_all_mos_at_r(r1,mos_array_r1) - call give_all_mos_at_r(r2,mos_array_r2) + call give_all_mos_at_r(r1,mos_array_r1) + call give_all_mos_at_r(r2,mos_array_r2) - ! You extract the inactive orbitals + ! You extract the inactive orbitals allocate( mos_array_inact_r1(n_inact_orb) , mos_array_inact_r2(n_inact_orb) ) do i = 1, n_inact_orb mos_array_inact_r1(i) = mos_array_r1(list_inact(i)) @@ -87,7 +87,7 @@ subroutine give_f_ia_val_ab(r1,r2,f_ia_val_ab,two_bod_dens,istate) mos_array_inact_r2(i) = mos_array_r2(list_inact(i)) enddo - ! You extract the active orbitals + ! You extract the active orbitals allocate( mos_array_act_r1(n_basis_orb) , mos_array_act_r2(n_basis_orb) ) do i= 1, n_act_orb mos_array_act_r1(i) = mos_array_r1(list_act(i)) @@ -109,11 +109,11 @@ subroutine give_f_ia_val_ab(r1,r2,f_ia_val_ab,two_bod_dens,istate) ! rho_tilde(i,a) = \sum_b rho(b,a) * phi_i(1) * phi_j(2) allocate(rho_tilde(n_inact_orb,n_act_orb)) two_bod_dens = 0.d0 - do a = 1, n_act_orb + do a = 1, n_act_orb do i = 1, n_inact_orb rho_tilde(i,a) = 0.d0 - do b = 1, n_act_orb - rho = one_e_act_dm_beta_mo_for_dft(b,a,istate) + one_e_act_dm_alpha_mo_for_dft(b,a,istate) + do b = 1, n_act_orb + rho = one_e_act_dm_beta_mo_for_dft(b,a,istate) + one_e_act_dm_alpha_mo_for_dft(b,a,istate) two_bod_dens += mos_array_inact_r1(i) * mos_array_inact_r1(i) * mos_array_act_r2(a) * mos_array_act_r2(b) * rho rho_tilde(i,a) += rho * mos_array_inact_r1(i) * mos_array_act_r2(b) enddo @@ -125,12 +125,12 @@ subroutine give_f_ia_val_ab(r1,r2,f_ia_val_ab,two_bod_dens,istate) allocate( v_tilde(n_act_orb,n_act_orb) ) allocate( integrals_array(mo_num,mo_num) ) v_tilde = 0.d0 - do a = 1, n_act_orb + do a = 1, n_act_orb orb_a = list_act(a) do i = 1, n_inact_orb v_tilde(i,a) = 0.d0 orb_i = list_inact(i) -! call get_mo_two_e_integrals_ij(orb_i,orb_a,mo_num,integrals_array,mo_integrals_map) +! call get_mo_two_e_integrals_ij(orb_i,orb_a,mo_num,integrals_array,mo_integrals_map) do m = 1, n_basis_orb do n = 1, n_basis_orb ! v_tilde(i,a) += integrals_array(n,m) * mos_array_basis_r2(n) * mos_array_basis_r1(m) @@ -146,14 +146,14 @@ subroutine give_f_ia_val_ab(r1,r2,f_ia_val_ab,two_bod_dens,istate) enddo enddo ! multiply by two to adapt to the N(N-1) normalization condition of the active two-rdm - f_ia_val_ab *= 2.d0 + f_ia_val_ab *= 2.d0 two_bod_dens *= 2.d0 end subroutine give_f_aa_val_ab(r1,r2,f_aa_val_ab,two_bod_dens,istate) BEGIN_DOC -! contribution from purely active orbitals to f_{\Psi^B}(r_1,r_2) for the "istate" state of a CAS wave function +! contribution from purely active orbitals to f_{\Psi^B}(r_1,r_2) for the "istate" state of a CAS wave function END_DOC implicit none integer, intent(in) :: istate @@ -161,7 +161,7 @@ subroutine give_f_aa_val_ab(r1,r2,f_aa_val_ab,two_bod_dens,istate) double precision, intent(out):: f_aa_val_ab,two_bod_dens integer :: i,orb_i,a,orb_a,n,m,b,c,d double precision :: rho - double precision, allocatable :: mos_array_r1(:) , mos_array_r2(:) + double precision, allocatable :: mos_array_r1(:) , mos_array_r2(:) double precision, allocatable :: mos_array_basis_r1(:),mos_array_basis_r2(:) double precision, allocatable :: mos_array_act_r1(:),mos_array_act_r2(:) double precision, allocatable :: integrals_array(:,:),rho_tilde(:,:),v_tilde(:,:) @@ -170,10 +170,10 @@ subroutine give_f_aa_val_ab(r1,r2,f_aa_val_ab,two_bod_dens,istate) two_bod_dens = 0.d0 ! You get all orbitals in r_1 and r_2 allocate(mos_array_r1(mo_num) , mos_array_r2(mo_num) ) - call give_all_mos_at_r(r1,mos_array_r1) - call give_all_mos_at_r(r2,mos_array_r2) + call give_all_mos_at_r(r1,mos_array_r1) + call give_all_mos_at_r(r2,mos_array_r2) - ! You extract the active orbitals + ! You extract the active orbitals allocate( mos_array_act_r1(n_basis_orb) , mos_array_act_r2(n_basis_orb) ) do i= 1, n_act_orb mos_array_act_r1(i) = mos_array_r1(list_act(i)) @@ -201,8 +201,8 @@ subroutine give_f_aa_val_ab(r1,r2,f_aa_val_ab,two_bod_dens,istate) do c = 1, n_act_orb ! 1 do d = 1, n_act_orb ! 2 rho = mos_array_act_r1(c) * mos_array_act_r2(d) * act_2_rdm_ab_mo(d,c,b,a,istate) - rho_tilde(b,a) += rho - two_bod_dens += rho * mos_array_act_r1(a) * mos_array_act_r2(b) + rho_tilde(b,a) += rho + two_bod_dens += rho * mos_array_act_r1(a) * mos_array_act_r2(b) enddo enddo enddo @@ -212,7 +212,7 @@ subroutine give_f_aa_val_ab(r1,r2,f_aa_val_ab,two_bod_dens,istate) ! v_tilde(i,a) = \sum_{m,n} phi_m(1) * phi_n(2) < i a | m n > allocate( v_tilde(n_act_orb,n_act_orb) ) v_tilde = 0.d0 - do a = 1, n_act_orb + do a = 1, n_act_orb do b = 1, n_act_orb v_tilde(b,a) = 0.d0 do m = 1, n_basis_orb @@ -235,92 +235,92 @@ end BEGIN_PROVIDER [double precision, two_e_int_aa_f, (n_basis_orb,n_basis_orb,n_act_orb,n_act_orb)] implicit none BEGIN_DOC -! list of two-electron integrals (built with the MOs belonging to the \mathcal{B} space) -! -! needed to compute the function f_{ii}(r_1,r_2) +! list of two-electron integrals (built with the MOs belonging to the \mathcal{B} space) +! +! needed to compute the function f_{ii}(r_1,r_2) ! ! two_e_int_aa_f(j,i,n,m) = < j i | n m > where all orbitals belong to "list_basis" END_DOC integer :: orb_i,orb_j,i,j,orb_m,orb_n,m,n double precision :: integrals_array(mo_num,mo_num),get_two_e_integral - PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals - do orb_m = 1, n_act_orb ! electron 1 + PROVIDE all_mo_integrals + do orb_m = 1, n_act_orb ! electron 1 m = list_act(orb_m) - do orb_n = 1, n_act_orb ! electron 2 + do orb_n = 1, n_act_orb ! electron 2 n = list_act(orb_n) - call get_mo_two_e_integrals_ij(m,n,mo_num,integrals_array,mo_integrals_map) - do orb_i = 1, n_basis_orb ! electron 1 + call get_mo_two_e_integrals_ij(m,n,mo_num,integrals_array,mo_integrals_map) + do orb_i = 1, n_basis_orb ! electron 1 i = list_basis(orb_i) - do orb_j = 1, n_basis_orb ! electron 2 + do orb_j = 1, n_basis_orb ! electron 2 j = list_basis(orb_j) ! 2 1 2 1 - two_e_int_aa_f(orb_j,orb_i,orb_n,orb_m) = get_two_e_integral(m,n,i,j,mo_integrals_map) -! two_e_int_aa_f(orb_j,orb_i,orb_n,orb_m) = integrals_array(j,i) + two_e_int_aa_f(orb_j,orb_i,orb_n,orb_m) = get_two_e_integral(m,n,i,j,mo_integrals_map) +! two_e_int_aa_f(orb_j,orb_i,orb_n,orb_m) = integrals_array(j,i) enddo enddo enddo enddo -END_PROVIDER +END_PROVIDER BEGIN_PROVIDER [double precision, two_e_int_ia_f, (n_basis_orb,n_basis_orb,n_inact_orb,n_act_orb)] implicit none BEGIN_DOC -! list of two-electron integrals (built with the MOs belonging to the \mathcal{B} space) -! -! needed to compute the function f_{ia}(r_1,r_2) +! list of two-electron integrals (built with the MOs belonging to the \mathcal{B} space) +! +! needed to compute the function f_{ia}(r_1,r_2) ! ! two_e_int_aa_f(j,i,n,m) = < j i | n m > where all orbitals belong to "list_basis" END_DOC integer :: orb_i,orb_j,i,j,orb_m,orb_n,m,n double precision :: integrals_array(mo_num,mo_num),get_two_e_integral - PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals - do orb_m = 1, n_act_orb ! electron 1 + PROVIDE all_mo_integrals + do orb_m = 1, n_act_orb ! electron 1 m = list_act(orb_m) - do orb_n = 1, n_inact_orb ! electron 2 + do orb_n = 1, n_inact_orb ! electron 2 n = list_inact(orb_n) - call get_mo_two_e_integrals_ij(m,n,mo_num,integrals_array,mo_integrals_map) - do orb_i = 1, n_basis_orb ! electron 1 + call get_mo_two_e_integrals_ij(m,n,mo_num,integrals_array,mo_integrals_map) + do orb_i = 1, n_basis_orb ! electron 1 i = list_basis(orb_i) - do orb_j = 1, n_basis_orb ! electron 2 + do orb_j = 1, n_basis_orb ! electron 2 j = list_basis(orb_j) ! 2 1 2 1 -! two_e_int_ia_f(orb_j,orb_i,orb_n,orb_m) = get_two_e_integral(m,n,i,j,mo_integrals_map) - two_e_int_ia_f(orb_j,orb_i,orb_n,orb_m) = integrals_array(j,i) +! two_e_int_ia_f(orb_j,orb_i,orb_n,orb_m) = get_two_e_integral(m,n,i,j,mo_integrals_map) + two_e_int_ia_f(orb_j,orb_i,orb_n,orb_m) = integrals_array(j,i) enddo enddo enddo enddo -END_PROVIDER +END_PROVIDER BEGIN_PROVIDER [double precision, two_e_int_ii_f, (n_basis_orb,n_basis_orb,n_inact_orb,n_inact_orb)] implicit none BEGIN_DOC -! list of two-electron integrals (built with the MOs belonging to the \mathcal{B} space) -! -! needed to compute the function f_{ii}(r_1,r_2) +! list of two-electron integrals (built with the MOs belonging to the \mathcal{B} space) +! +! needed to compute the function f_{ii}(r_1,r_2) ! ! two_e_int_ii_f(j,i,n,m) = < j i | n m > where all orbitals belong to "list_basis" END_DOC integer :: orb_i,orb_j,i,j,orb_m,orb_n,m,n double precision :: get_two_e_integral,integrals_array(mo_num,mo_num) - PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals - do orb_m = 1, n_inact_orb ! electron 1 + PROVIDE all_mo_integrals + do orb_m = 1, n_inact_orb ! electron 1 m = list_inact(orb_m) - do orb_n = 1, n_inact_orb ! electron 2 + do orb_n = 1, n_inact_orb ! electron 2 n = list_inact(orb_n) - call get_mo_two_e_integrals_ij(m,n,mo_num,integrals_array,mo_integrals_map) - do orb_i = 1, n_basis_orb ! electron 1 + call get_mo_two_e_integrals_ij(m,n,mo_num,integrals_array,mo_integrals_map) + do orb_i = 1, n_basis_orb ! electron 1 i = list_basis(orb_i) - do orb_j = 1, n_basis_orb ! electron 2 + do orb_j = 1, n_basis_orb ! electron 2 j = list_basis(orb_j) ! 2 1 2 1 -! two_e_int_ii_f(orb_j,orb_i,orb_n,orb_m) = get_two_e_integral(m,n,i,j,mo_integrals_map) - two_e_int_ii_f(orb_j,orb_i,orb_n,orb_m) = integrals_array(j,i) +! two_e_int_ii_f(orb_j,orb_i,orb_n,orb_m) = get_two_e_integral(m,n,i,j,mo_integrals_map) + two_e_int_ii_f(orb_j,orb_i,orb_n,orb_m) = integrals_array(j,i) enddo enddo enddo enddo -END_PROVIDER +END_PROVIDER subroutine give_mu_of_r_cas(r,istate,mu_of_r,f_psi,n2_psi) @@ -343,13 +343,13 @@ subroutine give_mu_of_r_cas(r,istate,mu_of_r,f_psi,n2_psi) ! active-active part of f_psi(r1,r2) call give_f_aa_val_ab(r,r,f_aa_val_ab,two_bod_dens_aa,istate) - f_psi = f_ii_val_ab + f_ia_val_ab + f_aa_val_ab + f_psi = f_ii_val_ab + f_ia_val_ab + f_aa_val_ab n2_psi = two_bod_dens_ii + two_bod_dens_ia + two_bod_dens_aa if(n2_psi.le.1.d-12.or.f_psi.le.0.d0.or.f_psi * n2_psi.lt.0.d0)then w_psi = 1.d+10 - else + else w_psi = f_psi / n2_psi endif mu_of_r = w_psi * sqpi * 0.5d0 - + end diff --git a/src/mu_of_r/mu_of_r_conditions.irp.f b/src/mu_of_r/mu_of_r_conditions.irp.f index f2bb7145..88dad8c3 100644 --- a/src/mu_of_r/mu_of_r_conditions.irp.f +++ b/src/mu_of_r/mu_of_r_conditions.irp.f @@ -1,12 +1,12 @@ BEGIN_PROVIDER [double precision, mu_of_r_prov, (n_points_final_grid,N_states) ] - implicit none + implicit none BEGIN_DOC - ! general variable for mu(r) + ! general variable for mu(r) ! - ! corresponds to Eq. (37) of J. Chem. Phys. 149, 194301 (2018) + ! corresponds to Eq. (37) of J. Chem. Phys. 149, 194301 (2018) ! - ! !!!!!! WARNING !!!!!! if no_core_density == .True. then all contributions from the core orbitals + ! !!!!!! WARNING !!!!!! if no_core_density == .True. then all contributions from the core orbitals ! ! in the two-body density matrix are excluded END_DOC @@ -31,7 +31,7 @@ mu_of_r_prov(ipoint,istate) = mu_of_r_hf_sparse(ipoint) else if(mu_of_r_potential.EQ."cas_full".or.mu_of_r_potential.EQ."cas_truncated".or.mu_of_r_potential.EQ."pure_act")then mu_of_r_prov(ipoint,istate) = mu_of_r_psi_cas(ipoint,istate) - else + else print*,'you requested the following mu_of_r_potential' print*,mu_of_r_potential print*,'which does not correspond to any of the options for such keyword' @@ -42,23 +42,23 @@ if (write_mu_of_r) then print*,'Writing mu(r) on disk ...' - call ezfio_set_mu_of_r_io_mu_of_r('Read') - call ezfio_set_mu_of_r_mu_of_r_disk(mu_of_r_prov) + call ezfio_set_mu_of_r_io_mu_of_r('Read') + call ezfio_set_mu_of_r_mu_of_r_disk(mu_of_r_prov) endif call wall_time(wall1) print*,'Time to provide mu_of_r = ',wall1-wall0 - END_PROVIDER + END_PROVIDER BEGIN_PROVIDER [double precision, mu_of_r_hf, (n_points_final_grid) ] - implicit none + implicit none BEGIN_DOC ! mu(r) computed with a HF wave function (assumes that HF MOs are stored in the EZFIO) ! ! corresponds to Eq. (37) of J. Chem. Phys. 149, 194301 (2018) but for \Psi^B = HF^B ! - ! !!!!!! WARNING !!!!!! if no_core_density == .True. then all contributions from the core orbitals + ! !!!!!! WARNING !!!!!! if no_core_density == .True. then all contributions from the core orbitals ! ! in the two-body density matrix are excluded END_DOC @@ -70,14 +70,14 @@ sqpi = dsqrt(dacos(-1.d0)) !$OMP PARALLEL DO & !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (ipoint,f_hf,on_top,w_hf) & - !$OMP ShARED (n_points_final_grid,mu_of_r_hf,f_hf_cholesky,on_top_hf_grid,sqpi) + !$OMP PRIVATE (ipoint,f_hf,on_top,w_hf) & + !$OMP ShARED (n_points_final_grid,mu_of_r_hf,f_hf_cholesky,on_top_hf_grid,sqpi) do ipoint = 1, n_points_final_grid f_hf = f_hf_cholesky(ipoint) on_top = on_top_hf_grid(ipoint) if(on_top.le.1.d-12.or.f_hf.le.0.d0.or.f_hf * on_top.lt.0.d0)then w_hf = 1.d+10 - else + else w_hf = f_hf / on_top endif mu_of_r_hf(ipoint) = w_hf * sqpi * 0.5d0 @@ -85,16 +85,16 @@ !$OMP END PARALLEL DO call wall_time(wall1) print*,'Time to provide mu_of_r_hf = ',wall1-wall0 - END_PROVIDER + END_PROVIDER BEGIN_PROVIDER [double precision, mu_of_r_hf_sparse, (n_points_final_grid) ] - implicit none + implicit none BEGIN_DOC ! mu(r) computed with a HF wave function (assumes that HF MOs are stored in the EZFIO) ! ! corresponds to Eq. (37) of J. Chem. Phys. 149, 194301 (2018) but for \Psi^B = HF^B ! - ! !!!!!! WARNING !!!!!! if no_core_density == .True. then all contributions from the core orbitals + ! !!!!!! WARNING !!!!!! if no_core_density == .True. then all contributions from the core orbitals ! ! in the two-body density matrix are excluded END_DOC @@ -106,14 +106,14 @@ PROVIDE f_hf_cholesky_sparse on_top_hf_grid !$OMP PARALLEL DO & !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (ipoint,f_hf,on_top,w_hf) & - !$OMP ShARED (n_points_final_grid,mu_of_r_hf_sparse,f_hf_cholesky_sparse,on_top_hf_grid,sqpi) + !$OMP PRIVATE (ipoint,f_hf,on_top,w_hf) & + !$OMP ShARED (n_points_final_grid,mu_of_r_hf_sparse,f_hf_cholesky_sparse,on_top_hf_grid,sqpi) do ipoint = 1, n_points_final_grid f_hf = f_hf_cholesky_sparse(ipoint) on_top = on_top_hf_grid(ipoint) if(on_top.le.1.d-12.or.f_hf.le.0.d0.or.f_hf * on_top.lt.0.d0)then w_hf = 1.d+10 - else + else w_hf = f_hf / on_top endif mu_of_r_hf_sparse(ipoint) = w_hf * sqpi * 0.5d0 @@ -121,36 +121,36 @@ !$OMP END PARALLEL DO call wall_time(wall1) print*,'Time to provide mu_of_r_hf_sparse = ',wall1-wall0 - END_PROVIDER + END_PROVIDER BEGIN_PROVIDER [double precision, mu_of_r_hf_old, (n_points_final_grid) ] - implicit none + implicit none BEGIN_DOC ! mu(r) computed with a HF wave function (assumes that HF MOs are stored in the EZFIO) ! ! corresponds to Eq. (37) of J. Chem. Phys. 149, 194301 (2018) but for \Psi^B = HF^B ! - ! !!!!!! WARNING !!!!!! if no_core_density == .True. then all contributions from the core orbitals + ! !!!!!! WARNING !!!!!! if no_core_density == .True. then all contributions from the core orbitals ! ! in the two-body density matrix are excluded END_DOC integer :: ipoint double precision :: wall0,wall1,f_hf,on_top,w_hf,sqpi - PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals + PROVIDE all_mo_integrals print*,'providing mu_of_r_hf_old ...' call wall_time(wall0) sqpi = dsqrt(dacos(-1.d0)) - provide f_psi_hf_ab + provide f_psi_hf_ab !$OMP PARALLEL DO & !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (ipoint,f_hf,on_top,w_hf) & - !$OMP ShARED (n_points_final_grid,mu_of_r_hf_old,f_psi_hf_ab,on_top_hf_mu_r,sqpi) + !$OMP PRIVATE (ipoint,f_hf,on_top,w_hf) & + !$OMP ShARED (n_points_final_grid,mu_of_r_hf_old,f_psi_hf_ab,on_top_hf_mu_r,sqpi) do ipoint = 1, n_points_final_grid f_hf = f_psi_hf_ab(ipoint) on_top = on_top_hf_mu_r(ipoint) if(on_top.le.1.d-12.or.f_hf.le.0.d0.or.f_hf * on_top.lt.0.d0)then w_hf = 1.d+10 - else + else w_hf = f_hf / on_top endif mu_of_r_hf_old(ipoint) = w_hf * sqpi * 0.5d0 @@ -158,17 +158,17 @@ !$OMP END PARALLEL DO call wall_time(wall1) print*,'Time to provide mu_of_r_hf_old = ',wall1-wall0 - END_PROVIDER + END_PROVIDER BEGIN_PROVIDER [double precision, mu_of_r_psi_cas, (n_points_final_grid,N_states) ] - implicit none + implicit none BEGIN_DOC ! mu(r) computed with a wave function developped in an active space ! - ! corresponds to Eq. (37) of J. Chem. Phys. 149, 194301 (2018) + ! corresponds to Eq. (37) of J. Chem. Phys. 149, 194301 (2018) ! - ! !!!!!! WARNING !!!!!! if no_core_density == .True. then all contributions from the core orbitals + ! !!!!!! WARNING !!!!!! if no_core_density == .True. then all contributions from the core orbitals ! ! in the one- and two-body density matrix are excluded END_DOC @@ -181,15 +181,15 @@ provide f_psi_cas_ab !$OMP PARALLEL DO & !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (ipoint,f_psi,on_top,w_psi,istate) & - !$OMP SHARED (n_points_final_grid,mu_of_r_psi_cas,f_psi_cas_ab,on_top_cas_mu_r,sqpi,N_states) + !$OMP PRIVATE (ipoint,f_psi,on_top,w_psi,istate) & + !$OMP SHARED (n_points_final_grid,mu_of_r_psi_cas,f_psi_cas_ab,on_top_cas_mu_r,sqpi,N_states) do istate = 1, N_states do ipoint = 1, n_points_final_grid - f_psi = f_psi_cas_ab(ipoint,istate) - on_top = on_top_cas_mu_r(ipoint,istate) + f_psi = f_psi_cas_ab(ipoint,istate) + on_top = on_top_cas_mu_r(ipoint,istate) if(on_top.le.1.d-12.or.f_psi.le.0.d0.or.f_psi * on_top.lt.0.d0)then w_psi = 1.d+10 - else + else w_psi = f_psi / on_top endif mu_of_r_psi_cas(ipoint,istate) = w_psi * sqpi * 0.5d0 @@ -198,15 +198,15 @@ !$OMP END PARALLEL DO call wall_time(wall1) print*,'Time to provide mu_of_r_psi_cas = ',wall1-wall0 - END_PROVIDER + END_PROVIDER BEGIN_PROVIDER [double precision, mu_average_prov, (N_states)] implicit none BEGIN_DOC - ! average value of mu(r) weighted with the total one-e density and divided by the number of electrons + ! average value of mu(r) weighted with the total one-e density and divided by the number of electrons ! - ! !!!!!! WARNING !!!!!! if no_core_density == .True. then all contributions from the core orbitals + ! !!!!!! WARNING !!!!!! if no_core_density == .True. then all contributions from the core orbitals ! ! in the one- and two-body density matrix are excluded END_DOC @@ -216,12 +216,12 @@ do istate = 1, N_states do ipoint = 1, n_points_final_grid weight =final_weight_at_r_vector(ipoint) - density = one_e_dm_and_grad_alpha_in_r(4,ipoint,istate) & + density = one_e_dm_and_grad_alpha_in_r(4,ipoint,istate) & + one_e_dm_and_grad_beta_in_r(4,ipoint,istate) if(mu_of_r_prov(ipoint,istate).gt.1.d+09)cycle mu_average_prov(istate) += mu_of_r_prov(ipoint,istate) * weight * density enddo mu_average_prov(istate) = mu_average_prov(istate) / elec_num_grid_becke(istate) enddo - END_PROVIDER + END_PROVIDER diff --git a/src/tools/fcidump.irp.f b/src/tools/fcidump.irp.f index bf4d07fb..df050218 100644 --- a/src/tools/fcidump.irp.f +++ b/src/tools/fcidump.irp.f @@ -41,7 +41,7 @@ program fcidump integer(key_kind), allocatable :: keys(:) double precision, allocatable :: values(:) integer(cache_map_size_kind) :: n_elements, n_elements_max - PROVIDE mo_two_e_integrals_in_map + PROVIDE all_mo_integrals double precision :: get_two_e_integral, integral diff --git a/src/tools/fcidump_pyscf.irp.f b/src/tools/fcidump_pyscf.irp.f index aaa552b4..e6962a60 100644 --- a/src/tools/fcidump_pyscf.irp.f +++ b/src/tools/fcidump_pyscf.irp.f @@ -41,7 +41,7 @@ program fcidump_pyscf integer(key_kind), allocatable :: keys(:) double precision, allocatable :: values(:) integer(cache_map_size_kind) :: n_elements, n_elements_max - PROVIDE mo_two_e_integrals_in_map + PROVIDE all_mo_integrals double precision :: get_two_e_integral, integral diff --git a/src/tools/four_idx_transform.irp.f b/src/tools/four_idx_transform.irp.f index fc6bface..cd1db0c0 100644 --- a/src/tools/four_idx_transform.irp.f +++ b/src/tools/four_idx_transform.irp.f @@ -18,6 +18,6 @@ program four_idx_transform io_mo_two_e_integrals = 'Write' SOFT_TOUCH io_mo_two_e_integrals if (.true.) then - PROVIDE mo_two_e_integrals_in_map + PROVIDE all_mo_integrals endif end diff --git a/src/trexio/export_trexio_routines.irp.f b/src/trexio/export_trexio_routines.irp.f index 5bc44880..0774bcd9 100644 --- a/src/trexio/export_trexio_routines.irp.f +++ b/src/trexio/export_trexio_routines.irp.f @@ -505,7 +505,7 @@ subroutine export_trexio(update,full_path) if (export_mo_two_e_ints) then print *, 'MO two-e integrals' - PROVIDE mo_two_e_integrals_in_map + PROVIDE all_mo_integrals double precision, external :: mo_two_e_integral diff --git a/src/two_body_rdm/act_2_rdm.irp.f b/src/two_body_rdm/act_2_rdm.irp.f index c550e991..9e2ea018 100644 --- a/src/two_body_rdm/act_2_rdm.irp.f +++ b/src/two_body_rdm/act_2_rdm.irp.f @@ -4,34 +4,34 @@ BEGIN_DOC ! 12 12 ! 1 2 1 2 == -! act_2_rdm_ab_mo(i,j,k,l,istate) = STATE SPECIFIC physicist notation for 2RDM of alpha/beta+beta/alpha electrons -! +! act_2_rdm_ab_mo(i,j,k,l,istate) = STATE SPECIFIC physicist notation for 2RDM of alpha/beta+beta/alpha electrons +! ! ! ! + ! -! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO AN ACTIVE SPACE DEFINED BY "list_act" +! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO AN ACTIVE SPACE DEFINED BY "list_act" ! ! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha}^{act} * N_{\beta}^{act} * 2 ! -! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" +! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" ! - END_DOC + END_DOC integer :: ispin double precision :: wall_1, wall_2 - character*(128) :: name_file + character*(128) :: name_file name_file = 'act_2_rdm_ab_mo' ! condition for alpha/beta spin print*,'' print*,'Providing act_2_rdm_ab_mo ' - ispin = 3 + ispin = 3 act_2_rdm_ab_mo = 0.d0 - provide mo_two_e_integrals_in_map + provide all_mo_integrals call wall_time(wall_1) if(read_two_body_rdm_ab)then print*,'Reading act_2_rdm_ab_mo from disk ...' call read_array_two_rdm(n_act_orb,N_states,act_2_rdm_ab_mo,name_file) - else + else call orb_range_2_rdm_openmp(act_2_rdm_ab_mo,n_act_orb,n_act_orb,list_act,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) endif if(write_two_body_rdm_ab)then @@ -42,36 +42,36 @@ call wall_time(wall_2) print*,'Wall time to provide act_2_rdm_ab_mo',wall_2 - wall_1 act_2_rdm_ab_mo *= 2.d0 - END_PROVIDER + END_PROVIDER BEGIN_PROVIDER [double precision, act_2_rdm_aa_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb,N_states)] implicit none BEGIN_DOC -! act_2_rdm_aa_mo(i,j,k,l,istate) = STATE SPECIFIC physicist notation for 2RDM of ALPHA/ALPHA electrons -! +! act_2_rdm_aa_mo(i,j,k,l,istate) = STATE SPECIFIC physicist notation for 2RDM of ALPHA/ALPHA electrons +! ! ! -! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO AN ACTIVE SPACE DEFINED BY "list_act" +! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO AN ACTIVE SPACE DEFINED BY "list_act" ! ! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha}^{act} * (N_{\alpha}^{act} - 1) ! -! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" - END_DOC +! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" + END_DOC integer :: ispin double precision :: wall_1, wall_2 ! condition for alpha/beta spin print*,'' print*,'Providing act_2_rdm_aa_mo ' - character*(128) :: name_file + character*(128) :: name_file name_file = 'act_2_rdm_aa_mo' - ispin = 1 + ispin = 1 act_2_rdm_aa_mo = 0.d0 call wall_time(wall_1) if(read_two_body_rdm_aa)then print*,'Reading act_2_rdm_aa_mo from disk ...' call read_array_two_rdm(n_act_orb,N_states,act_2_rdm_aa_mo,name_file) - else + else call orb_range_2_rdm_openmp(act_2_rdm_aa_mo,n_act_orb,n_act_orb,list_act,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) endif if(write_two_body_rdm_aa)then @@ -83,36 +83,36 @@ call wall_time(wall_2) print*,'Wall time to provide act_2_rdm_aa_mo',wall_2 - wall_1 act_2_rdm_aa_mo *= 2.d0 - END_PROVIDER - + END_PROVIDER + BEGIN_PROVIDER [double precision, act_2_rdm_bb_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb,N_states)] implicit none BEGIN_DOC -! act_2_rdm_bb_mo(i,j,k,l,istate) = STATE SPECIFIC physicist notation for 2RDM of BETA/BETA electrons -! +! act_2_rdm_bb_mo(i,j,k,l,istate) = STATE SPECIFIC physicist notation for 2RDM of BETA/BETA electrons +! ! ! -! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO AN ACTIVE SPACE DEFINED BY "list_act" +! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO AN ACTIVE SPACE DEFINED BY "list_act" ! ! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\beta}^{act} * (N_{\beta}^{act} - 1) ! -! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" - END_DOC +! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" + END_DOC integer :: ispin double precision :: wall_1, wall_2 ! condition for beta/beta spin print*,'' print*,'Providing act_2_rdm_bb_mo ' - character*(128) :: name_file + character*(128) :: name_file name_file = 'act_2_rdm_bb_mo' - ispin = 2 + ispin = 2 act_2_rdm_bb_mo = 0.d0 call wall_time(wall_1) if(read_two_body_rdm_bb)then print*,'Reading act_2_rdm_bb_mo from disk ...' call read_array_two_rdm(n_act_orb,N_states,act_2_rdm_bb_mo,name_file) - else + else call orb_range_2_rdm_openmp(act_2_rdm_bb_mo,n_act_orb,n_act_orb,list_act,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) endif if(write_two_body_rdm_bb)then @@ -124,35 +124,35 @@ call wall_time(wall_2) print*,'Wall time to provide act_2_rdm_bb_mo',wall_2 - wall_1 act_2_rdm_bb_mo *= 2.d0 - END_PROVIDER + END_PROVIDER BEGIN_PROVIDER [double precision, act_2_rdm_spin_trace_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb,N_states)] implicit none BEGIN_DOC -! act_2_rdm_spin_trace_mo(i,j,k,l,istate) = STATE SPECIFIC physicist notation for 2RDM -! +! act_2_rdm_spin_trace_mo(i,j,k,l,istate) = STATE SPECIFIC physicist notation for 2RDM +! ! \sum_{\sigma,\sigma'} ! -! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO AN ACTIVE SPACE DEFINED BY "list_act" +! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO AN ACTIVE SPACE DEFINED BY "list_act" ! ! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{elec}^{act} * (N_{elec}^{act} - 1) ! -! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" - END_DOC +! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" + END_DOC integer :: ispin double precision :: wall_1, wall_2 ! condition for beta/beta spin print*,'' print*,'Providing act_2_rdm_spin_trace_mo ' - character*(128) :: name_file + character*(128) :: name_file name_file = 'act_2_rdm_spin_trace_mo' - ispin = 4 + ispin = 4 act_2_rdm_spin_trace_mo = 0.d0 call wall_time(wall_1) if(read_two_body_rdm_spin_trace)then print*,'Reading act_2_rdm_spin_trace_mo from disk ...' call read_array_two_rdm(n_act_orb,N_states,act_2_rdm_spin_trace_mo,name_file) - else + else call orb_range_2_rdm_openmp(act_2_rdm_spin_trace_mo,n_act_orb,n_act_orb,list_act,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) endif if(write_two_body_rdm_spin_trace)then @@ -164,4 +164,4 @@ act_2_rdm_spin_trace_mo *= 2.d0 call wall_time(wall_2) print*,'Wall time to provide act_2_rdm_spin_trace_mo',wall_2 - wall_1 - END_PROVIDER + END_PROVIDER diff --git a/src/two_rdm_routines/davidson_like_2rdm.irp.f b/src/two_rdm_routines/davidson_like_2rdm.irp.f index ad7a3b21..09436663 100644 --- a/src/two_rdm_routines/davidson_like_2rdm.irp.f +++ b/src/two_rdm_routines/davidson_like_2rdm.irp.f @@ -2,9 +2,9 @@ subroutine orb_range_2_rdm_openmp(big_array,dim1,norb,list_orb,ispin,u_0,N_st,sz use bitmasks implicit none BEGIN_DOC - ! if ispin == 1 :: alpha/alpha 2rdm - ! == 2 :: beta /beta 2rdm - ! == 3 :: alpha/beta + beta/alpha 2rdm + ! if ispin == 1 :: alpha/alpha 2rdm + ! == 2 :: beta /beta 2rdm + ! == 3 :: alpha/beta + beta/alpha 2rdm ! == 4 :: spin traced 2rdm :: aa + bb + ab + ba ! ! Assumes that the determinants are in psi_det @@ -19,7 +19,7 @@ subroutine orb_range_2_rdm_openmp(big_array,dim1,norb,list_orb,ispin,u_0,N_st,sz integer :: k double precision, allocatable :: u_t(:,:) !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: u_t - PROVIDE mo_two_e_integrals_in_map + PROVIDE all_mo_integrals allocate(u_t(N_st,N_det)) do k=1,N_st call dset_order(u_0(1,k),psi_bilinear_matrix_order,N_det) @@ -30,14 +30,14 @@ subroutine orb_range_2_rdm_openmp(big_array,dim1,norb,list_orb,ispin,u_0,N_st,sz u_t, & size(u_t, 1), & N_det, N_st) - + call orb_range_2_rdm_openmp_work(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,1,N_det,0,1) deallocate(u_t) - + do k=1,N_st call dset_order(u_0(1,k),psi_bilinear_matrix_order_reverse,N_det) enddo - + end subroutine orb_range_2_rdm_openmp_work(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,istart,iend,ishift,istep) @@ -52,11 +52,11 @@ subroutine orb_range_2_rdm_openmp_work(big_array,dim1,norb,list_orb,ispin,u_t,N_ integer, intent(in) :: dim1,norb,list_orb(norb),ispin double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1) double precision, intent(in) :: u_t(N_st,N_det) - + integer :: k - + PROVIDE N_int - + select case (N_int) case (1) call orb_range_2_rdm_openmp_work_1(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,istart,iend,ishift,istep) @@ -70,9 +70,9 @@ subroutine orb_range_2_rdm_openmp_work(big_array,dim1,norb,list_orb,ispin,u_t,N_ call orb_range_2_rdm_openmp_work_N_int(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,istart,iend,ishift,istep) end select end - - - + + + BEGIN_TEMPLATE subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,istart,iend,ishift,istep) @@ -81,9 +81,9 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin implicit none BEGIN_DOC ! Computes the two rdm for the N_st vectors |u_t> - ! if ispin == 1 :: alpha/alpha 2rdm - ! == 2 :: beta /beta 2rdm - ! == 3 :: alpha/beta 2rdm + ! if ispin == 1 :: alpha/alpha 2rdm + ! == 2 :: beta /beta 2rdm + ! == 3 :: alpha/beta 2rdm ! == 4 :: spin traced 2rdm :: aa + bb + 0.5 (ab + ba)) ! The 2rdm will be computed only on the list of orbitals list_orb, which contains norb ! Default should be 1,N_det,0,1 for istart,iend,ishift,istep @@ -92,7 +92,7 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin double precision, intent(in) :: u_t(N_st,N_det) integer, intent(in) :: dim1,norb,list_orb(norb),ispin double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1) - + integer(omp_lock_kind) :: lock_2rdm integer :: i,j,k,l integer :: k_a, k_b, l_a, l_b @@ -116,7 +116,7 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin integer, allocatable :: keys(:,:) double precision, allocatable :: values(:,:) integer :: nkeys,sze_buff - alpha_alpha = .False. + alpha_alpha = .False. beta_beta = .False. alpha_beta = .False. spin_trace = .False. @@ -134,27 +134,27 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin stop endif - + PROVIDE N_int call list_to_bitstring( orb_bitmask, list_orb, norb, N_int) - sze_buff = 6 * norb + elec_alpha_num * elec_alpha_num * 60 - list_orb_reverse = -1000 + sze_buff = 6 * norb + elec_alpha_num * elec_alpha_num * 60 + list_orb_reverse = -1000 do i = 1, norb - list_orb_reverse(list_orb(i)) = i + list_orb_reverse(list_orb(i)) = i enddo maxab = max(N_det_alpha_unique, N_det_beta_unique)+1 allocate(idx0(maxab)) - + do i=1,maxab idx0(i) = i enddo call omp_init_lock(lock_2rdm) - + ! Prepare the array of all alpha single excitations ! ------------------------------------------------- - - PROVIDE N_int nthreads_davidson elec_alpha_num + + PROVIDE N_int nthreads_davidson elec_alpha_num !$OMP PARALLEL DEFAULT(NONE) NUM_THREADS(nthreads_davidson) & !$OMP SHARED(psi_bilinear_matrix_rows, N_det,lock_2rdm,& !$OMP psi_bilinear_matrix_columns, & @@ -166,7 +166,7 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin !$OMP psi_bilinear_matrix_order_transp_reverse, & !$OMP psi_bilinear_matrix_columns_loc, & !$OMP psi_bilinear_matrix_transp_rows_loc,elec_alpha_num, & - !$OMP istart, iend, istep, irp_here,list_orb_reverse, n_states, dim1, & + !$OMP istart, iend, istep, irp_here,list_orb_reverse, n_states, dim1, & !$OMP ishift, idx0, u_t, maxab, alpha_alpha,beta_beta,alpha_beta,spin_trace,ispin,big_array,sze_buff,orb_bitmask) & !$OMP PRIVATE(krow, kcol, tmp_det, spindet, k_a, k_b, i,c_1, & !$OMP lcol, lrow, l_a, l_b, & @@ -174,35 +174,35 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin !$OMP tmp_det2, idx, l, kcol_prev, & !$OMP singles_a, n_singles_a, singles_b, & !$OMP n_singles_b, nkeys, keys, values) - + ! Alpha/Beta double excitations ! ============================= nkeys = 0 - allocate( keys(4,sze_buff), values(n_st,sze_buff)) + allocate( keys(4,sze_buff), values(n_st,sze_buff)) allocate( buffer($N_int,maxab), & singles_a(maxab), & singles_b(maxab), & doubles(maxab), & idx(maxab)) - + kcol_prev=-1 - + ASSERT (iend <= N_det) ASSERT (istart > 0) ASSERT (istep > 0) - + !$OMP DO SCHEDULE(dynamic,64) do k_a=istart+ishift,iend,istep - + krow = psi_bilinear_matrix_rows(k_a) ASSERT (krow <= N_det_alpha_unique) - + kcol = psi_bilinear_matrix_columns(k_a) ASSERT (kcol <= N_det_beta_unique) - + tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow) tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol) - + if (kcol /= kcol_prev) then call get_all_spin_singles_$N_int( & psi_det_beta_unique, idx0, & @@ -210,58 +210,58 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin singles_b, n_singles_b) endif kcol_prev = kcol - + ! Loop over singly excited beta columns ! ------------------------------------- - + do i=1,n_singles_b lcol = singles_b(i) - + tmp_det2(1:$N_int,2) = psi_det_beta_unique(1:$N_int, lcol) - + l_a = psi_bilinear_matrix_columns_loc(lcol) ASSERT (l_a <= N_det) - + do j=1,psi_bilinear_matrix_columns_loc(lcol+1) - l_a lrow = psi_bilinear_matrix_rows(l_a) ASSERT (lrow <= N_det_alpha_unique) - + buffer(1:$N_int,j) = psi_det_alpha_unique(1:$N_int, lrow) - + ASSERT (l_a <= N_det) idx(j) = l_a l_a = l_a+1 enddo j = j-1 - + call get_all_spin_singles_$N_int( & buffer, idx, tmp_det(1,1), j, & singles_a, n_singles_a ) - + ! Loop over alpha singles ! ----------------------- - + if(alpha_beta.or.spin_trace)then do k = 1,n_singles_a l_a = singles_a(k) ASSERT (l_a <= N_det) - + lrow = psi_bilinear_matrix_rows(l_a) ASSERT (lrow <= N_det_alpha_unique) - + tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, lrow) ! print*,'nkeys before = ',nkeys do l= 1, N_states c_1(l) = u_t(l,l_a) * u_t(l,k_a) enddo if(alpha_beta)then - ! only ONE contribution + ! only ONE contribution if (nkeys+1 .ge. sze_buff) then call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) nkeys = 0 endif else if (spin_trace)then - ! TWO contributions + ! TWO contributions if (nkeys+2 .ge. sze_buff) then call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) nkeys = 0 @@ -271,42 +271,42 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin enddo endif - + call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) nkeys = 0 enddo - + enddo !$OMP END DO - + !$OMP DO SCHEDULE(dynamic,64) do k_a=istart+ishift,iend,istep - - + + ! Single and double alpha exitations ! =================================== - - + + ! Initial determinant is at k_a in alpha-major representation ! ----------------------------------------------------------------------- - + krow = psi_bilinear_matrix_rows(k_a) ASSERT (krow <= N_det_alpha_unique) - + kcol = psi_bilinear_matrix_columns(k_a) ASSERT (kcol <= N_det_beta_unique) - + tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow) tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol) - + ! Initial determinant is at k_b in beta-major representation ! ---------------------------------------------------------------------- - + k_b = psi_bilinear_matrix_order_transp_reverse(k_a) ASSERT (k_b <= N_det) - + spindet(1:$N_int) = tmp_det(1:$N_int,1) - + ! Loop inside the beta column to gather all the connected alphas lcol = psi_bilinear_matrix_columns(k_a) l_a = psi_bilinear_matrix_columns_loc(lcol) @@ -316,28 +316,28 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin if (lcol /= kcol) exit lrow = psi_bilinear_matrix_rows(l_a) ASSERT (lrow <= N_det_alpha_unique) - + buffer(1:$N_int,i) = psi_det_alpha_unique(1:$N_int, lrow) idx(i) = l_a l_a = l_a+1 enddo i = i-1 - + call get_all_spin_singles_and_doubles_$N_int( & buffer, idx, spindet, i, & singles_a, doubles, n_singles_a, n_doubles ) - + ! Compute Hij for all alpha singles ! ---------------------------------- - + tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol) do i=1,n_singles_a l_a = singles_a(i) ASSERT (l_a <= N_det) - + lrow = psi_bilinear_matrix_rows(l_a) ASSERT (lrow <= N_det_alpha_unique) - + tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, lrow) do l= 1, N_states c_1(l) = u_t(l,l_a) * u_t(l,k_a) @@ -356,23 +356,23 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin endif call orb_range_off_diag_single_to_all_states_aa_dm_buffer(tmp_det,tmp_det2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) endif - + enddo - + call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) nkeys = 0 - + ! Compute Hij for all alpha doubles ! ---------------------------------- - + if(alpha_alpha.or.spin_trace)then do i=1,n_doubles l_a = doubles(i) ASSERT (l_a <= N_det) - + lrow = psi_bilinear_matrix_rows(l_a) ASSERT (lrow <= N_det_alpha_unique) - + do l= 1, N_states c_1(l) = u_t(l,l_a) * u_t(l,k_a) enddo @@ -385,29 +385,29 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin endif call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) nkeys = 0 - - + + ! Single and double beta excitations ! ================================== - - + + ! Initial determinant is at k_a in alpha-major representation ! ----------------------------------------------------------------------- - + krow = psi_bilinear_matrix_rows(k_a) kcol = psi_bilinear_matrix_columns(k_a) - + tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow) tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol) - + spindet(1:$N_int) = tmp_det(1:$N_int,2) - + ! Initial determinant is at k_b in beta-major representation ! ----------------------------------------------------------------------- - + k_b = psi_bilinear_matrix_order_transp_reverse(k_a) ASSERT (k_b <= N_det) - + ! Loop inside the alpha row to gather all the connected betas lrow = psi_bilinear_matrix_transp_rows(k_b) l_b = psi_bilinear_matrix_transp_rows_loc(lrow) @@ -417,28 +417,28 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin if (lrow /= krow) exit lcol = psi_bilinear_matrix_transp_columns(l_b) ASSERT (lcol <= N_det_beta_unique) - + buffer(1:$N_int,i) = psi_det_beta_unique(1:$N_int, lcol) idx(i) = l_b l_b = l_b+1 enddo i = i-1 - + call get_all_spin_singles_and_doubles_$N_int( & buffer, idx, spindet, i, & singles_b, doubles, n_singles_b, n_doubles ) - + ! Compute Hij for all beta singles ! ---------------------------------- - + tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow) do i=1,n_singles_b l_b = singles_b(i) ASSERT (l_b <= N_det) - + lcol = psi_bilinear_matrix_transp_columns(l_b) ASSERT (lcol <= N_det_beta_unique) - + tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, lcol) l_a = psi_bilinear_matrix_transp_order(l_b) do l= 1, N_states @@ -461,18 +461,18 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin enddo call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) nkeys = 0 - + ! Compute Hij for all beta doubles ! ---------------------------------- - + if(beta_beta.or.spin_trace)then do i=1,n_doubles l_b = doubles(i) ASSERT (l_b <= N_det) - + lcol = psi_bilinear_matrix_transp_columns(l_b) ASSERT (lcol <= N_det_beta_unique) - + l_a = psi_bilinear_matrix_transp_order(l_b) do l= 1, N_states c_1(l) = u_t(l,l_a) * u_t(l,k_a) @@ -484,59 +484,59 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin call orb_range_off_diag_double_to_all_states_bb_dm_buffer(tmp_det(1,2),psi_det_beta_unique(1, lcol),c_1,N_st,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) ! print*,'to do orb_range_off_diag_double_to_2_rdm_bb_dm_buffer' ASSERT (l_a <= N_det) - + enddo endif call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) nkeys = 0 - - + + ! Diagonal contribution ! ===================== - - + + ! Initial determinant is at k_a in alpha-major representation ! ----------------------------------------------------------------------- - + krow = psi_bilinear_matrix_rows(k_a) ASSERT (krow <= N_det_alpha_unique) - + kcol = psi_bilinear_matrix_columns(k_a) ASSERT (kcol <= N_det_beta_unique) - + tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow) tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol) - + double precision, external :: diag_wee_mat_elem, diag_S_mat_elem - + double precision :: c_1(N_states) do l = 1, N_states c_1(l) = u_t(l,k_a) * u_t(l,k_a) enddo - + call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) nkeys = 0 call orb_range_diag_to_all_states_2_rdm_dm_buffer(tmp_det,c_1,N_states,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) nkeys = 0 - + end do !$OMP END DO deallocate(buffer, singles_a, singles_b, doubles, idx, keys, values) !$OMP END PARALLEL - + end - + SUBST [ N_int ] - + 1;; 2;; 3;; 4;; N_int;; - + END_TEMPLATE - + subroutine update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) use omp_lib @@ -545,7 +545,7 @@ subroutine update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,loc integer, intent(in) :: keys(4,nkeys) double precision, intent(in) :: values(n_st,nkeys) double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1,n_st) - + integer(omp_lock_kind),intent(inout):: lock_2rdm integer :: istate diff --git a/src/two_rdm_routines/davidson_like_trans_2rdm.irp.f b/src/two_rdm_routines/davidson_like_trans_2rdm.irp.f index 9e68a0e1..7f6b9459 100644 --- a/src/two_rdm_routines/davidson_like_trans_2rdm.irp.f +++ b/src/two_rdm_routines/davidson_like_trans_2rdm.irp.f @@ -2,12 +2,12 @@ subroutine orb_range_2_trans_rdm_openmp(big_array,dim1,norb,list_orb,ispin,u_0,N use bitmasks implicit none BEGIN_DOC - ! if ispin == 1 :: alpha/alpha 2_rdm - ! == 2 :: beta /beta 2_rdm - ! == 3 :: alpha/beta + beta/alpha 2trans_rdm + ! if ispin == 1 :: alpha/alpha 2_rdm + ! == 2 :: beta /beta 2_rdm + ! == 3 :: alpha/beta + beta/alpha 2trans_rdm ! == 4 :: spin traced 2_rdm :: aa + bb + ab + ba ! - ! notice that here it is the TRANSITION RDM THAT IS COMPUTED + ! notice that here it is the TRANSITION RDM THAT IS COMPUTED ! ! THE DIAGONAL PART IS THE USUAL ONE FOR A GIVEN STATE ! Assumes that the determinants are in psi_det @@ -22,7 +22,7 @@ subroutine orb_range_2_trans_rdm_openmp(big_array,dim1,norb,list_orb,ispin,u_0,N integer :: k double precision, allocatable :: u_t(:,:) !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: u_t - PROVIDE mo_two_e_integrals_in_map + PROVIDE all_mo_integrals allocate(u_t(N_st,N_det)) do k=1,N_st call dset_order(u_0(1,k),psi_bilinear_matrix_order,N_det) @@ -33,14 +33,14 @@ subroutine orb_range_2_trans_rdm_openmp(big_array,dim1,norb,list_orb,ispin,u_0,N u_t, & size(u_t, 1), & N_det, N_st) - + call orb_range_2_trans_rdm_openmp_work(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,1,N_det,0,1) deallocate(u_t) - + do k=1,N_st call dset_order(u_0(1,k),psi_bilinear_matrix_order_reverse,N_det) enddo - + end subroutine orb_range_2_trans_rdm_openmp_work(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,istart,iend,ishift,istep) @@ -55,11 +55,11 @@ subroutine orb_range_2_trans_rdm_openmp_work(big_array,dim1,norb,list_orb,ispin, integer, intent(in) :: dim1,norb,list_orb(norb),ispin double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1,N_st,N_st) double precision, intent(in) :: u_t(N_st,N_det) - + integer :: k - + PROVIDE N_int - + select case (N_int) case (1) call orb_range_2_trans_rdm_openmp_work_1(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,istart,iend,ishift,istep) @@ -73,8 +73,8 @@ subroutine orb_range_2_trans_rdm_openmp_work(big_array,dim1,norb,list_orb,ispin, call orb_range_2_trans_rdm_openmp_work_N_int(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,istart,iend,ishift,istep) end select end - - + + BEGIN_TEMPLATE subroutine orb_range_2_trans_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin,u_t,N_st,sze,istart,iend,ishift,istep) use bitmasks @@ -82,9 +82,9 @@ subroutine orb_range_2_trans_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb implicit none BEGIN_DOC ! Computes the two trans_rdm for the N_st vectors |u_t> - ! if ispin == 1 :: alpha/alpha 2trans_rdm - ! == 2 :: beta /beta 2trans_rdm - ! == 3 :: alpha/beta 2trans_rdm + ! if ispin == 1 :: alpha/alpha 2trans_rdm + ! == 2 :: beta /beta 2trans_rdm + ! == 3 :: alpha/beta 2trans_rdm ! == 4 :: spin traced 2trans_rdm :: aa + bb + 0.5 (ab + ba)) ! The 2trans_rdm will be computed only on the list of orbitals list_orb, which contains norb ! Default should be 1,N_det,0,1 for istart,iend,ishift,istep @@ -93,7 +93,7 @@ subroutine orb_range_2_trans_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb double precision, intent(in) :: u_t(N_st,N_det) integer, intent(in) :: dim1,norb,list_orb(norb),ispin double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1,N_st,N_st) - + integer(omp_lock_kind) :: lock_2trans_rdm integer :: i,j,k,l integer :: k_a, k_b, l_a, l_b @@ -118,7 +118,7 @@ subroutine orb_range_2_trans_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb double precision, allocatable :: values(:,:,:) integer :: nkeys,sze_buff integer :: ll - alpha_alpha = .False. + alpha_alpha = .False. beta_beta = .False. alpha_beta = .False. spin_trace = .False. @@ -136,27 +136,27 @@ subroutine orb_range_2_trans_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb stop endif - + PROVIDE N_int call list_to_bitstring( orb_bitmask, list_orb, norb, N_int) - sze_buff = 6 * norb + elec_alpha_num * elec_alpha_num * 60 - list_orb_reverse = -1000 + sze_buff = 6 * norb + elec_alpha_num * elec_alpha_num * 60 + list_orb_reverse = -1000 do i = 1, norb - list_orb_reverse(list_orb(i)) = i + list_orb_reverse(list_orb(i)) = i enddo maxab = max(N_det_alpha_unique, N_det_beta_unique)+1 allocate(idx0(maxab)) - + do i=1,maxab idx0(i) = i enddo call omp_init_lock(lock_2trans_rdm) - + ! Prepare the array of all alpha single excitations ! ------------------------------------------------- - - PROVIDE N_int nthreads_davidson elec_alpha_num + + PROVIDE N_int nthreads_davidson elec_alpha_num !$OMP PARALLEL DEFAULT(NONE) NUM_THREADS(nthreads_davidson) & !$OMP SHARED(psi_bilinear_matrix_rows, N_det,lock_2trans_rdm,& !$OMP psi_bilinear_matrix_columns, & @@ -168,7 +168,7 @@ subroutine orb_range_2_trans_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb !$OMP psi_bilinear_matrix_order_transp_reverse, & !$OMP psi_bilinear_matrix_columns_loc, & !$OMP psi_bilinear_matrix_transp_rows_loc,elec_alpha_num, & - !$OMP istart, iend, istep, irp_here,list_orb_reverse, n_states, dim1, & + !$OMP istart, iend, istep, irp_here,list_orb_reverse, n_states, dim1, & !$OMP ishift, idx0, u_t, maxab, alpha_alpha,beta_beta,alpha_beta,spin_trace,ispin,big_array,sze_buff,orb_bitmask) & !$OMP PRIVATE(krow, kcol, tmp_det, spindet, k_a, k_b, i,c_1, & !$OMP lcol, lrow, l_a, l_b, & @@ -176,35 +176,35 @@ subroutine orb_range_2_trans_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb !$OMP tmp_det2, idx, l, kcol_prev, & !$OMP singles_a, n_singles_a, singles_b, & !$OMP n_singles_b, nkeys, keys, values) - + ! Alpha/Beta double excitations ! ============================= nkeys = 0 - allocate( keys(4,sze_buff), values(n_st,n_st,sze_buff)) + allocate( keys(4,sze_buff), values(n_st,n_st,sze_buff)) allocate( buffer($N_int,maxab), & singles_a(maxab), & singles_b(maxab), & doubles(maxab), & idx(maxab)) - + kcol_prev=-1 - + ASSERT (iend <= N_det) ASSERT (istart > 0) ASSERT (istep > 0) - + !$OMP DO SCHEDULE(dynamic,64) do k_a=istart+ishift,iend,istep - + krow = psi_bilinear_matrix_rows(k_a) ASSERT (krow <= N_det_alpha_unique) - + kcol = psi_bilinear_matrix_columns(k_a) ASSERT (kcol <= N_det_beta_unique) - + tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow) tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol) - + if (kcol /= kcol_prev) then call get_all_spin_singles_$N_int( & psi_det_beta_unique, idx0, & @@ -212,45 +212,45 @@ subroutine orb_range_2_trans_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb singles_b, n_singles_b) endif kcol_prev = kcol - + ! Loop over singly excited beta columns ! ------------------------------------- - + do i=1,n_singles_b lcol = singles_b(i) - + tmp_det2(1:$N_int,2) = psi_det_beta_unique(1:$N_int, lcol) - + l_a = psi_bilinear_matrix_columns_loc(lcol) ASSERT (l_a <= N_det) - + do j=1,psi_bilinear_matrix_columns_loc(lcol+1) - l_a lrow = psi_bilinear_matrix_rows(l_a) ASSERT (lrow <= N_det_alpha_unique) - + buffer(1:$N_int,j) = psi_det_alpha_unique(1:$N_int, lrow) - + ASSERT (l_a <= N_det) idx(j) = l_a l_a = l_a+1 enddo j = j-1 - + call get_all_spin_singles_$N_int( & buffer, idx, tmp_det(1,1), j, & singles_a, n_singles_a ) - + ! Loop over alpha singles ! ----------------------- - + if(alpha_beta.or.spin_trace)then do k = 1,n_singles_a l_a = singles_a(k) ASSERT (l_a <= N_det) - + lrow = psi_bilinear_matrix_rows(l_a) ASSERT (lrow <= N_det_alpha_unique) - + tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, lrow) ! print*,'nkeys before = ',nkeys do ll = 1, N_states @@ -259,13 +259,13 @@ subroutine orb_range_2_trans_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb enddo enddo if(alpha_beta)then - ! only ONE contribution + ! only ONE contribution if (nkeys+1 .ge. sze_buff) then call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm) nkeys = 0 endif else if (spin_trace)then - ! TWO contributions + ! TWO contributions if (nkeys+2 .ge. sze_buff) then call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm) nkeys = 0 @@ -275,42 +275,42 @@ subroutine orb_range_2_trans_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb enddo endif - + call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm) nkeys = 0 enddo - + enddo !$OMP END DO - + !$OMP DO SCHEDULE(dynamic,64) do k_a=istart+ishift,iend,istep - - + + ! Single and double alpha exitations ! =================================== - - + + ! Initial determinant is at k_a in alpha-major representation ! ----------------------------------------------------------------------- - + krow = psi_bilinear_matrix_rows(k_a) ASSERT (krow <= N_det_alpha_unique) - + kcol = psi_bilinear_matrix_columns(k_a) ASSERT (kcol <= N_det_beta_unique) - + tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow) tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol) - + ! Initial determinant is at k_b in beta-major representation ! ---------------------------------------------------------------------- - + k_b = psi_bilinear_matrix_order_transp_reverse(k_a) ASSERT (k_b <= N_det) - + spindet(1:$N_int) = tmp_det(1:$N_int,1) - + ! Loop inside the beta column to gather all the connected alphas lcol = psi_bilinear_matrix_columns(k_a) l_a = psi_bilinear_matrix_columns_loc(lcol) @@ -320,28 +320,28 @@ subroutine orb_range_2_trans_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb if (lcol /= kcol) exit lrow = psi_bilinear_matrix_rows(l_a) ASSERT (lrow <= N_det_alpha_unique) - + buffer(1:$N_int,i) = psi_det_alpha_unique(1:$N_int, lrow) idx(i) = l_a l_a = l_a+1 enddo i = i-1 - + call get_all_spin_singles_and_doubles_$N_int( & buffer, idx, spindet, i, & singles_a, doubles, n_singles_a, n_doubles ) - + ! Compute Hij for all alpha singles ! ---------------------------------- - + tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol) do i=1,n_singles_a l_a = singles_a(i) ASSERT (l_a <= N_det) - + lrow = psi_bilinear_matrix_rows(l_a) ASSERT (lrow <= N_det_alpha_unique) - + tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, lrow) do ll= 1, N_states do l= 1, N_states @@ -362,23 +362,23 @@ subroutine orb_range_2_trans_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb endif call orb_range_off_diag_single_to_all_states_aa_trans_rdm_buffer(tmp_det,tmp_det2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) endif - + enddo - + call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm) nkeys = 0 - + ! Compute Hij for all alpha doubles ! ---------------------------------- - + if(alpha_alpha.or.spin_trace)then do i=1,n_doubles l_a = doubles(i) ASSERT (l_a <= N_det) - + lrow = psi_bilinear_matrix_rows(l_a) ASSERT (lrow <= N_det_alpha_unique) - + do ll= 1, N_states do l= 1, N_states c_1(l,ll) = u_t(ll,l_a) * u_t(l,k_a) @@ -393,29 +393,29 @@ subroutine orb_range_2_trans_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb endif call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm) nkeys = 0 - - + + ! Single and double beta excitations ! ================================== - - + + ! Initial determinant is at k_a in alpha-major representation ! ----------------------------------------------------------------------- - + krow = psi_bilinear_matrix_rows(k_a) kcol = psi_bilinear_matrix_columns(k_a) - + tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow) tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol) - + spindet(1:$N_int) = tmp_det(1:$N_int,2) - + ! Initial determinant is at k_b in beta-major representation ! ----------------------------------------------------------------------- - + k_b = psi_bilinear_matrix_order_transp_reverse(k_a) ASSERT (k_b <= N_det) - + ! Loop inside the alpha row to gather all the connected betas lrow = psi_bilinear_matrix_transp_rows(k_b) l_b = psi_bilinear_matrix_transp_rows_loc(lrow) @@ -425,28 +425,28 @@ subroutine orb_range_2_trans_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb if (lrow /= krow) exit lcol = psi_bilinear_matrix_transp_columns(l_b) ASSERT (lcol <= N_det_beta_unique) - + buffer(1:$N_int,i) = psi_det_beta_unique(1:$N_int, lcol) idx(i) = l_b l_b = l_b+1 enddo i = i-1 - + call get_all_spin_singles_and_doubles_$N_int( & buffer, idx, spindet, i, & singles_b, doubles, n_singles_b, n_doubles ) - + ! Compute Hij for all beta singles ! ---------------------------------- - + tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow) do i=1,n_singles_b l_b = singles_b(i) ASSERT (l_b <= N_det) - + lcol = psi_bilinear_matrix_transp_columns(l_b) ASSERT (lcol <= N_det_beta_unique) - + tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, lcol) l_a = psi_bilinear_matrix_transp_order(l_b) do ll= 1, N_states @@ -471,18 +471,18 @@ subroutine orb_range_2_trans_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb enddo call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm) nkeys = 0 - + ! Compute Hij for all beta doubles ! ---------------------------------- - + if(beta_beta.or.spin_trace)then do i=1,n_doubles l_b = doubles(i) ASSERT (l_b <= N_det) - + lcol = psi_bilinear_matrix_transp_columns(l_b) ASSERT (lcol <= N_det_beta_unique) - + l_a = psi_bilinear_matrix_transp_order(l_b) do ll= 1, N_states do l= 1, N_states @@ -496,61 +496,61 @@ subroutine orb_range_2_trans_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb call orb_range_off_diag_double_to_all_states_trans_rdm_bb_buffer(tmp_det(1,2),psi_det_beta_unique(1, lcol),c_1,N_st,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) ! print*,'to do orb_range_off_diag_double_to_2_trans_rdm_bb_dm_buffer' ASSERT (l_a <= N_det) - + enddo endif call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm) nkeys = 0 - - + + ! Diagonal contribution ! ===================== - - + + ! Initial determinant is at k_a in alpha-major representation ! ----------------------------------------------------------------------- - + krow = psi_bilinear_matrix_rows(k_a) ASSERT (krow <= N_det_alpha_unique) - + kcol = psi_bilinear_matrix_columns(k_a) ASSERT (kcol <= N_det_beta_unique) - + tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow) tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol) - + double precision, external :: diag_wee_mat_elem, diag_S_mat_elem - + double precision :: c_1(N_states,N_states) do ll = 1, N_states do l = 1, N_states c_1(l,ll) = u_t(ll,k_a) * u_t(l,k_a) enddo enddo - + call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm) nkeys = 0 call orb_range_diag_to_all_states_2_rdm_trans_buffer(tmp_det,c_1,N_states,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) call update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2trans_rdm) nkeys = 0 - + end do !$OMP END DO deallocate(buffer, singles_a, singles_b, doubles, idx, keys, values) !$OMP END PARALLEL - + end - + SUBST [ N_int ] - + 1;; 2;; 3;; 4;; N_int;; - + END_TEMPLATE - + subroutine update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) use omp_lib implicit none @@ -558,7 +558,7 @@ subroutine update_keys_values_n_states_trans(keys,values,nkeys,dim1,n_st,big_arr integer, intent(in) :: keys(4,nkeys) double precision, intent(in) :: values(n_st,n_st,nkeys) double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1,n_st,n_st) - + integer(omp_lock_kind),intent(inout):: lock_2rdm integer :: i,h1,h2,p1,p2,istate,jstate diff --git a/src/utils_cc/mo_integrals_cc.irp.f b/src/utils_cc/mo_integrals_cc.irp.f index b2b68d05..6f21c316 100644 --- a/src/utils_cc/mo_integrals_cc.irp.f +++ b/src/utils_cc/mo_integrals_cc.irp.f @@ -77,7 +77,7 @@ subroutine gen_v_space(n1,n2,n3,n4,list1,list2,list3,list4,v) else double precision :: get_two_e_integral - PROVIDE mo_two_e_integrals_in_map + PROVIDE all_mo_integrals !$OMP PARALLEL & !$OMP SHARED(n1,n2,n3,n4,list1,list2,list3,list4,v,mo_integrals_map) & @@ -161,7 +161,7 @@ BEGIN_PROVIDER [double precision, cc_space_v, (mo_num,mo_num,mo_num,mo_num)] integer :: i,j,k,l double precision :: get_two_e_integral - PROVIDE mo_two_e_integrals_in_map + PROVIDE all_mo_integrals !$OMP PARALLEL & !$OMP SHARED(cc_space_v,mo_num,mo_integrals_map) & @@ -194,7 +194,7 @@ BEGIN_PROVIDER [double precision, cc_space_v_oooo, (cc_nOa, cc_nOa, cc_nOa, cc_n integer :: i1, i2, i3, i4 integer :: n1, n2, n3, n4 - + n1 = size(cc_space_v_oooo,1) n2 = size(cc_space_v_oooo,2) n3 = size(cc_space_v_oooo,3) @@ -237,7 +237,7 @@ BEGIN_PROVIDER [double precision, cc_space_v_vooo, (cc_nVa, cc_nOa, cc_nOa, cc_n integer :: i1, i2, i3, i4 integer :: n1, n2, n3, n4 - + n1 = size(cc_space_v_vooo,1) n2 = size(cc_space_v_vooo,2) n3 = size(cc_space_v_vooo,3) @@ -281,7 +281,7 @@ BEGIN_PROVIDER [double precision, cc_space_v_ovoo, (cc_nOa, cc_nVa, cc_nOa, cc_n integer :: i1, i2, i3, i4 integer :: n1, n2, n3, n4 - + n1 = size(cc_space_v_ovoo,1) n2 = size(cc_space_v_ovoo,2) n3 = size(cc_space_v_ovoo,3) @@ -315,7 +315,7 @@ BEGIN_PROVIDER [double precision, cc_space_v_oovo, (cc_nOa, cc_nOa, cc_nVa, cc_n integer :: i1, i2, i3, i4 integer :: n1, n2, n3, n4 - + n1 = size(cc_space_v_oovo,1) n2 = size(cc_space_v_oovo,2) n3 = size(cc_space_v_oovo,3) @@ -349,7 +349,7 @@ BEGIN_PROVIDER [double precision, cc_space_v_ooov, (cc_nOa, cc_nOa, cc_nOa, cc_n integer :: i1, i2, i3, i4 integer :: n1, n2, n3, n4 - + n1 = size(cc_space_v_oovo,1) n2 = size(cc_space_v_oovo,2) n3 = size(cc_space_v_oovo,3) @@ -383,7 +383,7 @@ BEGIN_PROVIDER [double precision, cc_space_v_vvoo, (cc_nVa, cc_nVa, cc_nOa, cc_n integer :: i1, i2, i3, i4 integer :: n1, n2, n3, n4 - + n1 = size(cc_space_v_vvoo,1) n2 = size(cc_space_v_vvoo,2) n3 = size(cc_space_v_vvoo,3) @@ -426,7 +426,7 @@ BEGIN_PROVIDER [double precision, cc_space_v_vovo, (cc_nVa, cc_nOa, cc_nVa, cc_n integer :: i1, i2, i3, i4 integer :: n1, n2, n3, n4 - + n1 = size(cc_space_v_vovo,1) n2 = size(cc_space_v_vovo,2) n3 = size(cc_space_v_vovo,3) @@ -469,7 +469,7 @@ BEGIN_PROVIDER [double precision, cc_space_v_voov, (cc_nVa, cc_nOa, cc_nOa, cc_n integer :: i1, i2, i3, i4 integer :: n1, n2, n3, n4 - + n1 = size(cc_space_v_voov,1) n2 = size(cc_space_v_voov,2) n3 = size(cc_space_v_voov,3) @@ -503,7 +503,7 @@ BEGIN_PROVIDER [double precision, cc_space_v_ovvo, (cc_nOa, cc_nVa, cc_nVa, cc_n integer :: i1, i2, i3, i4 integer :: n1, n2, n3, n4 - + n1 = size(cc_space_v_ovvo,1) n2 = size(cc_space_v_ovvo,2) n3 = size(cc_space_v_ovvo,3) @@ -537,7 +537,7 @@ BEGIN_PROVIDER [double precision, cc_space_v_ovov, (cc_nOa, cc_nVa, cc_nOa, cc_n integer :: i1, i2, i3, i4 integer :: n1, n2, n3, n4 - + n1 = size(cc_space_v_ovov,1) n2 = size(cc_space_v_ovov,2) n3 = size(cc_space_v_ovov,3) @@ -571,7 +571,7 @@ BEGIN_PROVIDER [double precision, cc_space_v_oovv, (cc_nOa, cc_nOa, cc_nVa, cc_n integer :: i1, i2, i3, i4 integer :: n1, n2, n3, n4 - + n1 = size(cc_space_v_oovv,1) n2 = size(cc_space_v_oovv,2) n3 = size(cc_space_v_oovv,3) From a7bf04962be7f52fe292a07e239e6f1ac9742a23 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 20 Jan 2025 19:21:30 +0100 Subject: [PATCH 05/10] Fixed size of mo_integrals_cache --- src/mo_two_e_ints/map_integrals.irp.f | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mo_two_e_ints/map_integrals.irp.f b/src/mo_two_e_ints/map_integrals.irp.f index 06300666..9f485b79 100644 --- a/src/mo_two_e_ints/map_integrals.irp.f +++ b/src/mo_two_e_ints/map_integrals.irp.f @@ -62,7 +62,7 @@ end END_PROVIDER -BEGIN_PROVIDER [ double precision, mo_integrals_cache, (0_8:mo_integrals_cache_size_8**4_8) ] +BEGIN_PROVIDER [ double precision, mo_integrals_cache, (0_8:mo_integrals_cache_size_8) ] implicit none BEGIN_DOC ! Cache of MO integrals for fast access From 360aae7a104abfc086e269afd23c057702c06bfd Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 20 Jan 2025 19:35:24 +0100 Subject: [PATCH 06/10] Bug axes fixed --- src/scf_utils/roothaan_hall_scf.irp.f | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/scf_utils/roothaan_hall_scf.irp.f b/src/scf_utils/roothaan_hall_scf.irp.f index 947917af..274f08d6 100644 --- a/src/scf_utils/roothaan_hall_scf.irp.f +++ b/src/scf_utils/roothaan_hall_scf.irp.f @@ -222,13 +222,13 @@ END_DOC endif - ! Identify degenerate MOs and force them on the axes + ! Identify degenerate MOs and force them to be on the axes allocate(S(ao_num,ao_num)) i=1 do while (i Date: Tue, 21 Jan 2025 16:04:20 +0100 Subject: [PATCH 07/10] Fixed TREXIO norm problems --- ocaml/Input_ao_basis.ml | 44 ------------------------- scripts/qp_import_trexio.py | 10 +++--- src/ao_basis/EZFIO.cfg | 12 ------- src/ao_basis/aos.irp.f | 3 +- src/basis/EZFIO.cfg | 12 +++++++ src/basis/basis.irp.f | 17 ++++++++-- src/trexio/export_trexio_routines.irp.f | 27 ++++++--------- 7 files changed, 44 insertions(+), 81 deletions(-) diff --git a/ocaml/Input_ao_basis.ml b/ocaml/Input_ao_basis.ml index d9e28e04..343a4ae0 100644 --- a/ocaml/Input_ao_basis.ml +++ b/ocaml/Input_ao_basis.ml @@ -13,8 +13,6 @@ module Ao_basis : sig ao_coef : AO_coef.t array; ao_expo : AO_expo.t array; ao_cartesian : bool; - ao_normalized : bool; - primitives_normalized : bool; } [@@deriving sexp] ;; val read : unit -> t option @@ -36,8 +34,6 @@ end = struct ao_coef : AO_coef.t array; ao_expo : AO_expo.t array; ao_cartesian : bool; - ao_normalized : bool; - primitives_normalized : bool; } [@@deriving sexp] ;; @@ -133,24 +129,6 @@ end = struct Ezfio.get_ao_basis_ao_cartesian () ;; - let read_ao_normalized () = - if not (Ezfio.has_ao_basis_ao_normalized()) then - get_default "ao_normalized" - |> bool_of_string - |> Ezfio.set_ao_basis_ao_normalized - ; - Ezfio.get_ao_basis_ao_normalized () - ;; - - let read_primitives_normalized () = - if not (Ezfio.has_ao_basis_primitives_normalized()) then - get_default "primitives_normalized" - |> bool_of_string - |> Ezfio.set_ao_basis_primitives_normalized - ; - Ezfio.get_ao_basis_primitives_normalized () - ;; - let to_long_basis b = let ao_num = AO_number.to_int b.ao_num in let gto_array = Array.init (AO_number.to_int b.ao_num) @@ -213,8 +191,6 @@ end = struct ao_coef ; ao_expo ; ao_cartesian ; - ao_normalized ; - primitives_normalized ; } = b in write_md5 b ; @@ -247,8 +223,6 @@ end = struct ~rank:2 ~dim:[| ao_num ; 3 |] ~data:ao_power) ; Ezfio.set_ao_basis_ao_cartesian(ao_cartesian); - Ezfio.set_ao_basis_ao_normalized(ao_normalized); - Ezfio.set_ao_basis_primitives_normalized(primitives_normalized); let ao_coef = Array.to_list ao_coef @@ -280,8 +254,6 @@ end = struct ao_coef = read_ao_coef () ; ao_expo = read_ao_expo () ; ao_cartesian = read_ao_cartesian () ; - ao_normalized = read_ao_normalized () ; - primitives_normalized = read_primitives_normalized () ; } in to_md5 result @@ -392,8 +364,6 @@ end = struct { ao_basis = name ; ao_num ; ao_prim_num ; ao_prim_num_max ; ao_nucl ; ao_power ; ao_coef ; ao_expo ; ao_cartesian ; - ao_normalized = bool_of_string @@ get_default "ao_normalized"; - primitives_normalized = bool_of_string @@ get_default "primitives_normalized"; } ;; @@ -448,14 +418,6 @@ Cartesian coordinates (6d,10f,...) :: ao_cartesian = %s -Use normalized primitive functions :: - - primitives_normalized = %s - -Use normalized basis functions :: - - ao_normalized = %s - Basis set (read-only) :: %s @@ -469,8 +431,6 @@ Basis set (read-only) :: " (AO_basis_name.to_string b.ao_basis) (string_of_bool b.ao_cartesian) - (string_of_bool b.primitives_normalized) - (string_of_bool b.ao_normalized) (Basis.to_string short_basis |> String_ext.split ~on:'\n' |> list_map (fun x-> " "^x) @@ -507,8 +467,6 @@ ao_power = %s ao_coef = %s ao_expo = %s ao_cartesian = %s -ao_normalized = %s -primitives_normalized = %s md5 = %s " (AO_basis_name.to_string b.ao_basis) @@ -525,8 +483,6 @@ md5 = %s (b.ao_expo |> Array.to_list |> list_map AO_expo.to_string |> String.concat ", ") (b.ao_cartesian |> string_of_bool) - (b.ao_normalized |> string_of_bool) - (b.primitives_normalized |> string_of_bool) (to_md5 b |> MD5.to_string ) ;; diff --git a/scripts/qp_import_trexio.py b/scripts/qp_import_trexio.py index fc76f8de..ef4ac12a 100755 --- a/scripts/qp_import_trexio.py +++ b/scripts/qp_import_trexio.py @@ -193,9 +193,11 @@ def write_ezfio(trexio_filename, filename): shell_factor = trexio.read_basis_shell_factor(trexio_file) prim_factor = trexio.read_basis_prim_factor(trexio_file) - for i,p in enumerate(prim_factor): - coefficient[i] *= prim_factor[i] - ezfio.set_ao_basis_primitives_normalized(False) + ezfio.set_basis_prim_normalization_factor(prim_factor) + ezfio.set_basis_primitives_normalized(True) + ezfio.set_basis_ao_normalized(False) + for i, shell_idx in enumerate(shell_index): + coefficient[i] *= shell_factor[shell_idx] ezfio.set_basis_prim_coef(coefficient) elif basis_type.lower() == "numerical": @@ -391,7 +393,7 @@ def write_ezfio(trexio_filename, filename): # Renormalize MO coefs if needed if trexio.has_ao_normalization(trexio_file_cart): - ezfio.set_ao_basis_ao_normalized(False) + ezfio.set_basis_ao_normalized(False) norm = trexio.read_ao_normalization(trexio_file_cart) # for j in range(mo_num): # for i,f in enumerate(norm): diff --git a/src/ao_basis/EZFIO.cfg b/src/ao_basis/EZFIO.cfg index bd716383..c22f8029 100644 --- a/src/ao_basis/EZFIO.cfg +++ b/src/ao_basis/EZFIO.cfg @@ -55,18 +55,6 @@ doc: If |true|, use |AOs| in Cartesian coordinates (6d,10f,...) interface: ezfio, provider default: false -[ao_normalized] -type: logical -doc: Use normalized basis functions -interface: ezfio, provider -default: true - -[primitives_normalized] -type: logical -doc: Use normalized primitive functions -interface: ezfio, provider -default: true - [use_cgtos] type: logical doc: If true, use cgtos for AO integrals diff --git a/src/ao_basis/aos.irp.f b/src/ao_basis/aos.irp.f index 34853398..d718e935 100644 --- a/src/ao_basis/aos.irp.f +++ b/src/ao_basis/aos.irp.f @@ -82,10 +82,11 @@ END_PROVIDER enddo ao_coef_normalization_factor(i) = 1.d0/dsqrt(norm) - if (.not.ao_normalized) then + if (ao_normalized) then do j=1,ao_prim_num(i) ao_coef_normalized(i,j) = ao_coef_normalized(i,j) * ao_coef_normalization_factor(i) enddo + else ao_coef_normalization_factor(i) = 1.d0 endif enddo diff --git a/src/basis/EZFIO.cfg b/src/basis/EZFIO.cfg index a6864418..03e224e4 100644 --- a/src/basis/EZFIO.cfg +++ b/src/basis/EZFIO.cfg @@ -73,3 +73,15 @@ size: (basis.prim_num) interface: ezfio, provider +[primitives_normalized] +type: logical +doc: If true, assume primitive basis functions are normalized +interface: ezfio, provider, ocaml +default: true + +[ao_normalized] +type: logical +doc: If true, normalize the basis functions +interface: ezfio, provider, ocaml +default: false + diff --git a/src/basis/basis.irp.f b/src/basis/basis.irp.f index b750d75a..5374e5be 100644 --- a/src/basis/basis.irp.f +++ b/src/basis/basis.irp.f @@ -6,6 +6,11 @@ BEGIN_PROVIDER [ double precision, shell_normalization_factor , (shell_num) ] logical :: has PROVIDE ezfio_filename + if (.not.ao_normalized) then + shell_normalization_factor = 1.d0 + return + endif + if (mpi_master) then if (size(shell_normalization_factor) == 0) return @@ -70,6 +75,12 @@ BEGIN_PROVIDER [ double precision, prim_normalization_factor , (prim_num) ] logical :: has PROVIDE ezfio_filename + + if (.not.primitives_normalized) then + prim_normalization_factor(:) = 1.d0 + return + endif + if (mpi_master) then if (size(prim_normalization_factor) == 0) return @@ -95,9 +106,9 @@ BEGIN_PROVIDER [ double precision, prim_normalization_factor , (prim_num) ] do k=1, prim_num if (shell_index(k) /= i) cycle - call overlap_gaussian_xyz(C_A,C_A,prim_expo(k),prim_expo(k), & - powA,powA,overlap_x,overlap_y,overlap_z,norm,nz) - prim_normalization_factor(k) = 1.d0/dsqrt(norm) + call overlap_gaussian_xyz(C_A,C_A,prim_expo(k),prim_expo(k), & + powA,powA,overlap_x,overlap_y,overlap_z,norm,nz) + prim_normalization_factor(k) = 1.d0/dsqrt(norm) enddo enddo diff --git a/src/trexio/export_trexio_routines.irp.f b/src/trexio/export_trexio_routines.irp.f index 0774bcd9..c60b1aa0 100644 --- a/src/trexio/export_trexio_routines.irp.f +++ b/src/trexio/export_trexio_routines.irp.f @@ -271,11 +271,7 @@ subroutine export_trexio(update,full_path) call trexio_assert(rc, TREXIO_SUCCESS) allocate(factor(shell_num)) -! if (ao_normalized) then - factor(1:shell_num) = shell_normalization_factor(1:shell_num) -! else -! factor(1:shell_num) = 1.d0 -! endif + factor(1:shell_num) = shell_normalization_factor(1:shell_num) rc = trexio_write_basis_shell_factor(f(1), factor) call trexio_assert(rc, TREXIO_SUCCESS) @@ -291,11 +287,12 @@ subroutine export_trexio(update,full_path) call trexio_assert(rc, TREXIO_SUCCESS) allocate(factor(prim_num)) -! if (primitives_normalized) then + if (primitives_normalized) then factor(1:prim_num) = prim_normalization_factor(1:prim_num) -! else -! factor(1:prim_num) = 1.d0 -! endif + else + factor(1:prim_num) = 1.d0 + endif + rc = trexio_write_basis_prim_factor(f(1), factor) call trexio_assert(rc, TREXIO_SUCCESS) deallocate(factor) @@ -324,14 +321,10 @@ subroutine export_trexio(update,full_path) C_A(3) = 0.d0 allocate(factor(ao_num)) - if (ao_normalized) then - do i=1,ao_num - l = ao_first_of_shell(ao_shell(i)) - factor(i) = (ao_coef_normalized(i,1)+tiny(1.d0))/(ao_coef_normalized(l,1)+tiny(1.d0)) - enddo - else - factor(:) = 1.d0 - endif + do i=1,ao_num + l = ao_first_of_shell(ao_shell(i)) + factor(i) = (ao_coef_normalized(i,1)+tiny(1.d0))/(ao_coef_normalized(l,1)+tiny(1.d0)) + enddo rc = trexio_write_ao_normalization(f(1), factor) call trexio_assert(rc, TREXIO_SUCCESS) deallocate(factor) From 1a5f0b3e38c4f32f345c547f5bbd3a69af08e0f2 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 21 Jan 2025 17:48:37 +0100 Subject: [PATCH 08/10] Accelerated cholesky integrals in pt2 --- src/mo_two_e_ints/map_integrals.irp.f | 48 +++++++++++++++++++++------ 1 file changed, 38 insertions(+), 10 deletions(-) diff --git a/src/mo_two_e_ints/map_integrals.irp.f b/src/mo_two_e_ints/map_integrals.irp.f index 9f485b79..a1471a1c 100644 --- a/src/mo_two_e_ints/map_integrals.irp.f +++ b/src/mo_two_e_ints/map_integrals.irp.f @@ -143,7 +143,7 @@ double precision function get_two_e_integral(i,j,k,l,map) END_DOC integer, intent(in) :: i,j,k,l integer(key_kind) :: idx - integer :: ii + integer :: ii, kk type(map_type), intent(inout) :: map real(integral_kind) :: tmp @@ -177,10 +177,10 @@ double precision function get_two_e_integral(i,j,k,l,map) if (do_mo_cholesky) then - double precision, external :: ddot - get_two_e_integral = ddot(cholesky_mo_num, cholesky_mo_transp(1,i,k), 1, cholesky_mo_transp(1,j,l), 1) -! double precision, external :: get_from_mo_cholesky_cache -! get_two_e_integral = get_from_mo_cholesky_cache(i,j,k,l,.False.) + get_two_e_integral = 0.d0 + do kk=1,cholesky_mo_num + get_two_e_integral = get_two_e_integral + cholesky_mo_transp(kk,i,k)*cholesky_mo_transp(kk,i,l) + enddo else @@ -516,11 +516,39 @@ subroutine get_mo_two_e_integrals_exch_ii(k,l,sze,out_val,map) if (do_mo_cholesky) then - double precision, external :: ddot - do i=1,sze - out_val(i) = ddot(cholesky_mo_num, cholesky_mo_transp(1,i,k), 1, & - cholesky_mo_transp(1,i,l), 1) - enddo + if ( (k>=mo_integrals_cache_min).and.(k<=mo_integrals_cache_max).and. & + (l>=mo_integrals_cache_min).and.(l<=mo_integrals_cache_max) ) then + + integer :: kk + + do i=1,mo_integrals_cache_min-1 + out_val(i) = 0.d0 + do kk=1,cholesky_mo_num + out_val(i) = out_val(i) + cholesky_mo_transp(kk,i,k)*cholesky_mo_transp(kk,i,l) + enddo + enddo + + do i=mo_integrals_cache_min,mo_integrals_cache_max + out_val(i) = get_two_e_integral_cache(i,i,k,l) + enddo + + do i=mo_integrals_cache_max, sze + out_val(i) = 0.d0 + do kk=1,cholesky_mo_num + out_val(i) = out_val(i) + cholesky_mo_transp(kk,i,k)*cholesky_mo_transp(kk,i,l) + enddo + enddo + + else + + do i=1,sze + out_val(i) = 0.d0 + do kk=1,cholesky_mo_num + out_val(i) = out_val(i) + cholesky_mo_transp(kk,i,k)*cholesky_mo_transp(kk,i,l) + enddo + enddo + + endif else From db443e0a5af02d2fca237f58e3a78bf64560a3cd Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 21 Jan 2025 18:20:43 +0100 Subject: [PATCH 09/10] Put back ddot for AMD --- src/mo_two_e_ints/map_integrals.irp.f | 42 +++++++++++++++++---------- 1 file changed, 26 insertions(+), 16 deletions(-) diff --git a/src/mo_two_e_ints/map_integrals.irp.f b/src/mo_two_e_ints/map_integrals.irp.f index a1471a1c..76f169b4 100644 --- a/src/mo_two_e_ints/map_integrals.irp.f +++ b/src/mo_two_e_ints/map_integrals.irp.f @@ -177,10 +177,13 @@ double precision function get_two_e_integral(i,j,k,l,map) if (do_mo_cholesky) then - get_two_e_integral = 0.d0 - do kk=1,cholesky_mo_num - get_two_e_integral = get_two_e_integral + cholesky_mo_transp(kk,i,k)*cholesky_mo_transp(kk,i,l) - enddo + double precision, external :: ddot + get_two_e_integral = ddot(cholesky_mo_num, cholesky_mo_transp(1,i,k), 1, cholesky_mo_transp(1,j,l), 1) + +! get_two_e_integral = 0.d0 +! do kk=1,cholesky_mo_num +! get_two_e_integral = get_two_e_integral + cholesky_mo_transp(kk,i,k)*cholesky_mo_transp(kk,i,l) +! enddo else @@ -519,13 +522,16 @@ subroutine get_mo_two_e_integrals_exch_ii(k,l,sze,out_val,map) if ( (k>=mo_integrals_cache_min).and.(k<=mo_integrals_cache_max).and. & (l>=mo_integrals_cache_min).and.(l<=mo_integrals_cache_max) ) then + double precision, external :: ddot integer :: kk do i=1,mo_integrals_cache_min-1 - out_val(i) = 0.d0 - do kk=1,cholesky_mo_num - out_val(i) = out_val(i) + cholesky_mo_transp(kk,i,k)*cholesky_mo_transp(kk,i,l) - enddo + out_val(i) = ddot(cholesky_mo_num, cholesky_mo_transp(1,i,k), 1, & + cholesky_mo_transp(1,i,l), 1) +! out_val(i) = 0.d0 +! do kk=1,cholesky_mo_num +! out_val(i) = out_val(i) + cholesky_mo_transp(kk,i,k)*cholesky_mo_transp(kk,i,l) +! enddo enddo do i=mo_integrals_cache_min,mo_integrals_cache_max @@ -533,19 +539,23 @@ subroutine get_mo_two_e_integrals_exch_ii(k,l,sze,out_val,map) enddo do i=mo_integrals_cache_max, sze - out_val(i) = 0.d0 - do kk=1,cholesky_mo_num - out_val(i) = out_val(i) + cholesky_mo_transp(kk,i,k)*cholesky_mo_transp(kk,i,l) - enddo + out_val(i) = ddot(cholesky_mo_num, cholesky_mo_transp(1,i,k), 1, & + cholesky_mo_transp(1,i,l), 1) +! out_val(i) = 0.d0 +! do kk=1,cholesky_mo_num +! out_val(i) = out_val(i) + cholesky_mo_transp(kk,i,k)*cholesky_mo_transp(kk,i,l) +! enddo enddo else do i=1,sze - out_val(i) = 0.d0 - do kk=1,cholesky_mo_num - out_val(i) = out_val(i) + cholesky_mo_transp(kk,i,k)*cholesky_mo_transp(kk,i,l) - enddo + out_val(i) = ddot(cholesky_mo_num, cholesky_mo_transp(1,i,k), 1, & + cholesky_mo_transp(1,i,l), 1) +! out_val(i) = 0.d0 +! do kk=1,cholesky_mo_num +! out_val(i) = out_val(i) + cholesky_mo_transp(kk,i,k)*cholesky_mo_transp(kk,i,l) +! enddo enddo endif From ef234305e701ae50bf600127631fc9a8d84f4792 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 22 Jan 2025 11:06:53 +0100 Subject: [PATCH 10/10] Fixed ao_normalization. CP2K OK --- scripts/qp_import_trexio.py | 20 ++++++++------------ 1 file changed, 8 insertions(+), 12 deletions(-) diff --git a/scripts/qp_import_trexio.py b/scripts/qp_import_trexio.py index ef4ac12a..23f48eef 100755 --- a/scripts/qp_import_trexio.py +++ b/scripts/qp_import_trexio.py @@ -191,11 +191,12 @@ def write_ezfio(trexio_filename, filename): ezfio.set_basis_nucleus_shell_num(nucl_shell_num) - shell_factor = trexio.read_basis_shell_factor(trexio_file) prim_factor = trexio.read_basis_prim_factor(trexio_file) ezfio.set_basis_prim_normalization_factor(prim_factor) ezfio.set_basis_primitives_normalized(True) ezfio.set_basis_ao_normalized(False) + + shell_factor = trexio.read_basis_shell_factor(trexio_file) for i, shell_idx in enumerate(shell_index): coefficient[i] *= shell_factor[shell_idx] ezfio.set_basis_prim_coef(coefficient) @@ -249,7 +250,6 @@ def write_ezfio(trexio_filename, filename): ezfio.set_basis_shell_index([x+1 for x in shell_index]) ezfio.set_basis_nucleus_shell_num(nucl_shell_num) - shell_factor = trexio.read_basis_shell_factor(trexio_file) else: raise TypeError @@ -319,13 +319,18 @@ def write_ezfio(trexio_filename, filename): exponent.append(expo[i]) num_prim.append(num_prim0[i]) - print (len(coefficient), ao_num) assert (len(coefficient) == ao_num) + ezfio.set_ao_basis_ao_power(power_x + power_y + power_z) ezfio.set_ao_basis_ao_prim_num(num_prim) prim_num_max = max( [ len(x) for x in coefficient ] ) + ao_normalization = trexio.read_ao_normalization(trexio_file_cart) + for i, coef in enumerate(coefficient): + for j in range(len(coef)): + coef[j] *= ao_normalization[i] + for i in range(ao_num): coefficient[i] += [0. for j in range(len(coefficient[i]), prim_num_max)] exponent [i] += [0. for j in range(len(exponent[i]), prim_num_max)] @@ -340,7 +345,6 @@ def write_ezfio(trexio_filename, filename): coef.append(coefficient[j]) expo.append(exponent[j]) -# ezfio.set_ao_basis_ao_prim_num_max(prim_num_max) ezfio.set_ao_basis_ao_coef(coef) ezfio.set_ao_basis_ao_expo(expo) @@ -390,14 +394,6 @@ def write_ezfio(trexio_filename, filename): # Read coefs from temporary cartesian file created in the AO section MoMatrix = trexio.read_mo_coefficient(trexio_file_cart) - - # Renormalize MO coefs if needed - if trexio.has_ao_normalization(trexio_file_cart): - ezfio.set_basis_ao_normalized(False) - norm = trexio.read_ao_normalization(trexio_file_cart) -# for j in range(mo_num): -# for i,f in enumerate(norm): -# MoMatrix[i,j] *= f ezfio.set_mo_basis_mo_coef(MoMatrix) mo_occ = [ 0. for i in range(mo_num) ]