qp2/src/cas_based_on_top/on_top_cas_prov.irp.f

77 lines
2.7 KiB
Fortran

subroutine act_on_top_on_grid_pt(ipoint,istate,pure_act_on_top_of_r)
implicit none
BEGIN_DOC
! act_on_top_on_grid_pt returns the purely ACTIVE part of the on top pair density
!
! at the grid point ipoint, for the state istate
END_DOC
integer, intent(in) :: ipoint,istate
double precision, intent(out) :: pure_act_on_top_of_r
double precision :: phi_i,phi_j,phi_k,phi_l
integer :: i,j,k,l
ASSERT (istate <= N_states)
pure_act_on_top_of_r = 0.d0
do l = 1, n_act_orb
phi_l = act_mos_in_r_array(l,ipoint)
do k = 1, n_act_orb
phi_k = act_mos_in_r_array(k,ipoint)
do j = 1, n_act_orb
phi_j = act_mos_in_r_array(j,ipoint)
do i = 1, n_act_orb
phi_i = act_mos_in_r_array(i,ipoint)
! 1 2 1 2
pure_act_on_top_of_r += act_2_rdm_ab_mo(i,j,k,l,istate) * phi_i * phi_j * phi_k * phi_l
enddo
enddo
enddo
enddo
end
BEGIN_PROVIDER [double precision, total_cas_on_top_density,(n_points_final_grid,N_states) ]
implicit none
BEGIN_DOC
! on top pair density :: n2(r,r) at each of the Becke's grid point of a CAS-BASED wf
!
! Contains all core/inact/act contribution.
!
! !!!!! WARNING !!!!! If no_core_density then you REMOVE ALL CONTRIBUTIONS COMING FROM THE CORE ORBITALS
END_DOC
integer :: i_point,istate
double precision :: wall_0,wall_1,core_inact_dm,pure_act_on_top_of_r
logical :: no_core
print*,'providing the total_cas_on_top_density'
! for parallelization
provide inact_density core_density one_e_act_density_beta one_e_act_density_alpha act_mos_in_r_array
i_point = 1
istate = 1
call act_on_top_on_grid_pt(i_point,istate,pure_act_on_top_of_r)
call wall_time(wall_0)
if(no_core_density)then
print*,'USING THE VALENCE ONLY TWO BODY DENSITY'
endif
do istate = 1, N_states
!$OMP PARALLEL DO &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i_point,core_inact_dm,pure_act_on_top_of_r) &
!$OMP SHARED(total_cas_on_top_density,n_points_final_grid,inact_density,core_density,one_e_act_density_beta,one_e_act_density_alpha,no_core_density,istate)
do i_point = 1, n_points_final_grid
call act_on_top_on_grid_pt(i_point,istate,pure_act_on_top_of_r)
if(no_core_density) then
core_inact_dm = inact_density(i_point)
else
core_inact_dm = (inact_density(i_point) + core_density(i_point))
endif
total_cas_on_top_density(i_point,istate) = pure_act_on_top_of_r + core_inact_dm * (one_e_act_density_beta(i_point,istate) + one_e_act_density_alpha(i_point,istate)) + core_inact_dm*core_inact_dm
enddo
!$OMP END PARALLEL DO
enddo
call wall_time(wall_1)
print*,'provided the total_cas_on_top_density'
print*,'Time to provide :',wall_1 - wall_0
END_PROVIDER