From 348f6067914983d8919644d464a0ad18e7af6dc4 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 20 Jan 2025 18:03:19 +0100 Subject: [PATCH] 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)