LCPQ
/
qp2
mirror of https://github.com/QuantumPackage/qp2.git
9
1
Fork 0
Code Releases Activity

added cas_based_on_top

This commit is contained in:
Emmanuel Giner 2020-04-07 11:42:29 +02:00
parent 2fc294cd56
commit 87b4eab907
13 changed files with 779 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
src/cas_based_on_top/NEED Normal file
View File

@ -0,0 +1,2 @@
two_body_rdm
dft_utils_in_r

9
src/cas_based_on_top/README.rst Normal file
View File

@ -0,0 +1,9 @@
==============
on_top_density
==============
This plugin proposes different routines/providers to compute the on-top pair density of a CAS-based wave function.
This means that all determinants in psi_det must belong to an active-space.
As usual, see the file "example.irp.f" to get introduced to the main providers/routines.

125
src/cas_based_on_top/c_i_a_v_mos.irp.f Normal file
View File

@ -0,0 +1,125 @@
BEGIN_PROVIDER[double precision, core_mos_in_r_array, (n_core_orb,n_points_final_grid)]
&BEGIN_PROVIDER[double precision, core_mos_in_r_array_transp,(n_points_final_grid,n_core_orb)]
implicit none
BEGIN_DOC
! all COREE MOs on the grid points, arranged in two different ways
END_DOC
integer :: i,j,k
do i = 1, n_core_orb
j = list_core(i)
do k = 1, n_points_final_grid
core_mos_in_r_array_transp(k,i) = mos_in_r_array_transp(k,j)
enddo
enddo
do k = 1, n_points_final_grid
do i = 1, n_core_orb
core_mos_in_r_array(i,k) = core_mos_in_r_array_transp(k,i)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER[double precision, inact_mos_in_r_array, (n_inact_orb,n_points_final_grid)]
&BEGIN_PROVIDER[double precision, inact_mos_in_r_array_transp,(n_points_final_grid,n_inact_orb)]
implicit none
BEGIN_DOC
! all INACTIVE MOs on the grid points, arranged in two different ways
END_DOC
integer :: i,j,k
do i = 1, n_inact_orb
j = list_inact(i)
do k = 1, n_points_final_grid
inact_mos_in_r_array_transp(k,i) = mos_in_r_array_transp(k,j)
enddo
enddo
do k = 1, n_points_final_grid
do i = 1, n_inact_orb
inact_mos_in_r_array(i,k) = inact_mos_in_r_array_transp(k,i)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER[double precision, act_mos_in_r_array, (n_act_orb,n_points_final_grid)]
&BEGIN_PROVIDER[double precision, act_mos_in_r_array_transp,(n_points_final_grid,n_act_orb)]
implicit none
BEGIN_DOC
! all ACTIVE MOs on the grid points, arranged in two different ways
END_DOC
integer :: i,j,k
do i = 1, n_act_orb
j = list_act(i)
do k = 1, n_points_final_grid
act_mos_in_r_array_transp(k,i) = mos_in_r_array_transp(k,j)
enddo
enddo
do k = 1, n_points_final_grid
do i = 1, n_act_orb
act_mos_in_r_array(i,k) = act_mos_in_r_array_transp(k,i)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER[double precision, virt_mos_in_r_array, (n_virt_orb,n_points_final_grid)]
&BEGIN_PROVIDER[double precision, virt_mos_in_r_array_transp,(n_points_final_grid,n_virt_orb)]
implicit none
BEGIN_DOC
! all VIRTUAL MOs on the grid points, arranged in two different ways
END_DOC
integer :: i,j,k
do i = 1, n_virt_orb
j = list_virt(i)
do k = 1, n_points_final_grid
virt_mos_in_r_array_transp(k,i) = mos_in_r_array_transp(k,j)
enddo
enddo
do k = 1, n_points_final_grid
do i = 1, n_virt_orb
virt_mos_in_r_array(i,k) = virt_mos_in_r_array_transp(k,i)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER[double precision, core_inact_act_mos_in_r_array, (n_core_inact_act_orb,n_points_final_grid)]
&BEGIN_PROVIDER[double precision, core_inact_act_mos_in_r_array_transp,(n_points_final_grid,n_core_inact_act_orb)]
implicit none
integer :: i,j,k
do i = 1, n_core_inact_act_orb
j = list_core_inact_act(i)
do k = 1, n_points_final_grid
core_inact_act_mos_in_r_array_transp(k,i) = mos_in_r_array_transp(k,j)
enddo
enddo
do k = 1, n_points_final_grid
do i = 1, n_core_inact_act_orb
core_inact_act_mos_in_r_array(i,k) = core_inact_act_mos_in_r_array_transp(k,i)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER[double precision, core_inact_act_mos_grad_in_r_array, (3,n_core_inact_act_orb,n_points_final_grid)]
implicit none
integer :: i,j,k,l
do i = 1, n_core_inact_act_orb
j = list_core_inact_act(i)
do k = 1, n_points_final_grid
do l = 1, 3
core_inact_act_mos_grad_in_r_array(l,i,k) = mos_grad_in_r_array(j,k,l)
enddo
enddo
enddo
END_PROVIDER

47
src/cas_based_on_top/cas_based_density.irp.f Normal file
View File

@ -0,0 +1,47 @@
program cas_based_density
implicit none
BEGIN_DOC
! TODO : Small example to use the different quantities in this plugin
END_DOC
!! You force QP2 to read the wave function in the EZFIO folder
!! It is assumed that all Slater determinants in the wave function
!! belongs to an active space defined by core, inactive and active list of orbitals
read_wf = .True.
touch read_wf
call routine_test_cas_based_density
end
subroutine routine_test_cas_based_density
implicit none
integer :: ipoint, istate
double precision :: accu_n_elec(N_states),accu_n_elec_2(N_states)
! PROVIDERS
accu_n_elec = 0.d0
do istate = 1, N_states
do ipoint = 1, n_points_final_grid
accu_n_elec(istate) += one_e_cas_total_density(ipoint,istate) * final_weight_at_r_vector(ipoint)
enddo
print*,'istate = ',istate
print*,'accu_n_elec = ',accu_n_elec(istate)
enddo
! ROUTINES
double precision :: r(3),core_dens,inact_dens,act_dens(2,N_states),total_cas_dens(N_states)
accu_n_elec_2 = 0.d0
do ipoint = 1, n_points_final_grid
r(:) = final_grid_points(:,ipoint)
call give_cas_density_in_r(core_dens,inact_dens,act_dens,total_cas_dens,r)
do istate = 1, N_states
accu_n_elec_2(istate) += total_cas_dens(istate) * final_weight_at_r_vector(ipoint)
enddo
enddo
do istate = 1, N_states
print*,'istate = ',istate
print*,'accu_n_elec = ',accu_n_elec_2(istate)
enddo
end

19
src/cas_based_on_top/cas_based_on_top.irp.f Normal file
View File

@ -0,0 +1,19 @@
program cas_based_on_top_density
implicit none
BEGIN_DOC
! TODO : Small example to use the different quantities in this plugin
END_DOC
!! You force QP2 to read the wave function in the EZFIO folder
!! It is assumed that all Slater determinants in the wave function
!! belongs to an active space defined by core, inactive and active list of orbitals
read_wf = .True.
touch read_wf
! call routine_test_cas_based_on_top_density
call routine
end
subroutine routine
implicit none
provide total_cas_on_top_density
end

99
src/cas_based_on_top/cas_dens_prov.irp.f Normal file
View File

@ -0,0 +1,99 @@
BEGIN_PROVIDER [double precision, one_e_cas_total_density ,(n_points_final_grid,N_states) ]
implicit none
BEGIN_DOC
! one_e_cas_total_density = TOTAL DENSITY FOR a CAS wave function
!
! WARNING : if "no_core_density" == .True. then the core part of density is ignored
END_DOC
integer :: ipoint,i,j,istate
do istate = 1, N_states
do ipoint = 1, n_points_final_grid
one_e_cas_total_density(ipoint,istate) = one_e_act_density_alpha(ipoint,istate) + one_e_act_density_beta(ipoint,istate) &
+ 2.d0 * inact_density(ipoint)
if(.not.no_core_density)then !!! YOU ADD THE CORE DENSITY
one_e_cas_total_density(ipoint,istate) += 2.d0 * core_density(ipoint)
endif
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, one_e_act_density_alpha,(n_points_final_grid,N_states) ]
implicit none
BEGIN_DOC
! one_e_act_density_alpha = pure ACTIVE part of the DENSITY for ALPHA ELECTRONS
END_DOC
one_e_act_density_alpha = 0.d0
integer :: ipoint,i,j,istate
do istate = 1, N_states
do ipoint = 1, n_points_final_grid
do i = 1, n_act_orb
do j = 1, n_act_orb
one_e_act_density_alpha(ipoint,istate) += one_e_act_dm_alpha_mo_for_dft(j,i,istate) * act_mos_in_r_array(j,ipoint) * act_mos_in_r_array(i,ipoint)
enddo
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, one_e_act_density_beta,(n_points_final_grid,N_states) ]
implicit none
BEGIN_DOC
! one_e_act_density_beta = pure ACTIVE part of the DENSITY for BETA ELECTRONS
END_DOC
one_e_act_density_beta = 0.d0
integer :: ipoint,i,j,istate
do istate = 1, N_states
do ipoint = 1, n_points_final_grid
do i = 1, n_act_orb
do j = 1, n_act_orb
one_e_act_density_beta(ipoint,istate) += one_e_act_dm_beta_mo_for_dft(j,i,istate) * act_mos_in_r_array(j,ipoint) * act_mos_in_r_array(i,ipoint)
enddo
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, inact_density, (n_points_final_grid) ]
implicit none
BEGIN_DOC
! INACTIVE part of the density for alpha/beta.
!
! WARNING :: IF YOU NEED THE TOTAL DENSITY COMING FROM THE INACTIVE,
!
! YOU MUST MULTIPLY BY TWO
END_DOC
integer :: i,j
inact_density = 0.d0
do i = 1, n_points_final_grid
do j = 1, n_inact_orb
inact_density(i) += inact_mos_in_r_array(j,i) **2
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, core_density, (n_points_final_grid) ]
implicit none
BEGIN_DOC
! CORE part of the density for alpha/beta.
!
! WARNING :: IF YOU NEED THE TOTAL DENSITY COMING FROM THE CORE,
!
! YOU MUST MULTIPLY BY TWO
END_DOC
integer :: i,j
core_density = 0.d0
do i = 1, n_points_final_grid
do j = 1, n_core_orb
core_density(i) += core_mos_in_r_array(j,i) **2
enddo
enddo
END_PROVIDER

65
src/cas_based_on_top/cas_dens_rout.irp.f Normal file
View File

@ -0,0 +1,65 @@
subroutine give_cas_density_in_r(core_dens,inact_dens,act_dens,total_cas_dens,r)
implicit none
BEGIN_DOC
! returns the different component of the density at grid point r(3) for a CAS wave function
!
! core_dens : density coming from the CORE orbitals
!
! inact_dens : density coming from the INACT orbitals
!
! act_dens(1/2,1:N_states) : active part of the alpha/beta electrons for all states
!
! total_cas_dens : total density of the cas wave function
!
! WARNING : if "no_core_density" == .True. then the core part of density is ignored in total_cas_dens
END_DOC
double precision, intent(in) :: r(3)
double precision, intent(out) :: core_dens, inact_dens, act_dens(2,N_states), total_cas_dens(N_states)
double precision, allocatable :: mos_array(:),act_mos(:)
allocate(mos_array(mo_num))
call give_all_mos_at_r(r,mos_array)
integer :: i,iorb,j,jorb,istate
! core part of the density
core_dens = 0.d0
do i = 1, n_core_orb
iorb = list_core(i)
core_dens += mos_array(iorb)*mos_array(iorb)
enddo
core_dens = core_dens * 2.d0
! inactive part of the density
inact_dens = 0.d0
do i = 1, n_inact_orb
iorb = list_inact(i)
inact_dens += mos_array(iorb)*mos_array(iorb)
enddo
inact_dens = inact_dens * 2.d0
allocate(act_mos(n_act_orb))
do i = 1, n_act_orb
iorb = list_act(i)
act_mos(i) = mos_array(iorb)
enddo
! active part of the density for alpha/beta and all states
act_dens = 0.d0
do istate = 1, N_states
do i = 1, n_act_orb
do j = 1, n_act_orb
act_dens(1,istate) += one_e_act_dm_alpha_mo_for_dft(j,i,istate) * act_mos(j) * act_mos(i)
act_dens(2,istate) += one_e_act_dm_beta_mo_for_dft(j,i,istate) * act_mos(j) * act_mos(i)
enddo
enddo
enddo
! TOTAL density for all states
do istate = 1, N_states
total_cas_dens(istate) = inact_dens + act_dens(1,istate) + act_dens(2,istate)
if(.not.no_core_density)then !!! YOU ADD THE CORE DENSITY
total_cas_dens(istate) += core_dens
endif
enddo
end

37
src/cas_based_on_top/cas_one_e_rdm.irp.f Normal file
View File

@ -0,0 +1,37 @@
BEGIN_PROVIDER [double precision, one_e_act_dm_beta_mo_for_dft, (n_act_orb,n_act_orb,N_states)]
implicit none
BEGIN_DOC
! one_e_act_dm_beta_mo_for_dft = pure ACTIVE part of the ONE ELECTRON REDUCED DENSITY MATRIX for the BETA ELECTRONS
END_DOC
integer :: i,j,ii,jj,istate
do istate = 1, N_states
do ii = 1, n_act_orb
i = list_act(ii)
do jj = 1, n_act_orb
j = list_act(jj)
one_e_act_dm_beta_mo_for_dft(jj,ii,istate) = one_e_dm_mo_beta_for_dft(j,i,istate)
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, one_e_act_dm_alpha_mo_for_dft, (n_act_orb,n_act_orb,N_states)]
implicit none
BEGIN_DOC
! one_e_act_dm_alpha_mo_for_dft = pure ACTIVE part of the ONE ELECTRON REDUCED DENSITY MATRIX for the ALPHA ELECTRONS
END_DOC
integer :: i,j,ii,jj,istate
do istate = 1, N_states
do ii = 1, n_act_orb
i = list_act(ii)
do jj = 1, n_act_orb
j = list_act(jj)
one_e_act_dm_alpha_mo_for_dft(jj,ii,istate) = one_e_dm_mo_alpha_for_dft(j,i,istate)
enddo
enddo
enddo
END_PROVIDER

68
src/cas_based_on_top/eff_spin_dens.irp.f Normal file
View File

@ -0,0 +1,68 @@
BEGIN_PROVIDER [double precision, effective_spin_dm, (n_points_final_grid,N_states) ]
&BEGIN_PROVIDER [double precision, grad_effective_spin_dm, (3,n_points_final_grid,N_states) ]
implicit none
BEGIN_DOC
! effective_spin_dm(r_i) = \sqrt( n(r)^2 - 4 * ontop(r) )
! effective spin density obtained from the total density and on-top pair density
! see equation (6) of Phys. Chem. Chem. Phys., 2015, 17, 22412--22422 | 22413
END_DOC
provide total_cas_on_top_density
integer :: i_point,i_state,i
double precision :: n2,m2,thr
thr = 1.d-14
effective_spin_dm = 0.d0
grad_effective_spin_dm = 0.d0
do i_state = 1, N_states
do i_point = 1, n_points_final_grid
n2 = (one_e_dm_and_grad_alpha_in_r(4,i_point,i_state) + one_e_dm_and_grad_beta_in_r(4,i_point,i_state))
! density squared
n2 = n2 * n2
if(n2 - 4.D0 * total_cas_on_top_density(i_point,i_state).gt.thr)then
effective_spin_dm(i_point,i_state) = dsqrt(n2 - 4.D0 * total_cas_on_top_density(i_point,i_state))
if(isnan(effective_spin_dm(i_point,i_state)))then
print*,'isnan(effective_spin_dm(i_point,i_state)'
stop
endif
m2 = effective_spin_dm(i_point,i_state)
m2 = 0.5d0 / m2 ! 1/(2 * sqrt(n(r)^2 - 4 * ontop(r)) )
do i = 1, 3
grad_effective_spin_dm(i,i_point,i_state) = m2 * ( one_e_stuff_for_pbe(i,i_point,i_state) - 4.d0 * grad_total_cas_on_top_density(i,i_point,i_state) )
enddo
else
effective_spin_dm(i_point,i_state) = 0.d0
grad_effective_spin_dm(:,i_point,i_state) = 0.d0
endif
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, effective_alpha_dm, (n_points_final_grid,N_states) ]
&BEGIN_PROVIDER [double precision, effective_beta_dm, (n_points_final_grid,N_states) ]
&BEGIN_PROVIDER [double precision, grad_effective_alpha_dm, (3,n_points_final_grid,N_states) ]
&BEGIN_PROVIDER [double precision, grad_effective_beta_dm, (3,n_points_final_grid,N_states) ]
implicit none
BEGIN_DOC
! effective_alpha_dm(r_i) = 1/2 * (effective_spin_dm(r_i) + n(r_i))
! effective_beta_dm(r_i) = 1/2 * (-effective_spin_dm(r_i) + n(r_i))
END_DOC
provide total_cas_on_top_density
integer :: i_point,i_state,i
double precision :: n,grad_n
do i_state = 1, N_states
do i_point = 1, n_points_final_grid
n = (one_e_dm_and_grad_alpha_in_r(4,i_point,i_state) + one_e_dm_and_grad_beta_in_r(4,i_point,i_state))
effective_alpha_dm(i_point,i_state) = 0.5d0 * (n + effective_spin_dm(i_point,i_state))
effective_beta_dm(i_point,i_state) = 0.5d0 * (n - effective_spin_dm(i_point,i_state))
do i = 1, 3
grad_n = (one_e_dm_and_grad_alpha_in_r(i,i_point,i_state) + one_e_dm_and_grad_beta_in_r(i,i_point,i_state))
grad_effective_alpha_dm(i,i_point,i_state) = 0.5d0 * (grad_n + grad_effective_spin_dm(i,i_point,i_state) )
grad_effective_beta_dm(i,i_point,i_state) = 0.5d0 * (grad_n - grad_effective_spin_dm(i,i_point,i_state) )
enddo
enddo
enddo
END_PROVIDER

41
src/cas_based_on_top/example.irp.f Normal file
View File

@ -0,0 +1,41 @@
subroutine write_on_top_in_real_space
implicit none
BEGIN_DOC
! This routines is a simple example of how to plot the on-top pair density on a simple 1D grid
END_DOC
double precision :: zmax,dz,r(3),on_top_in_r,total_density,zcenter,dist
integer :: nz,i,istate
character*(128) :: output
integer :: i_unit_output,getUnitAndOpen
PROVIDE ezfio_filename
output=trim(ezfio_filename)//'.on_top'
print*,'output = ',trim(output)
i_unit_output = getUnitAndOpen(output,'w')
zmax = 2.0d0
print*,'nucl_coord(1,3) = ',nucl_coord(1,3)
print*,'nucl_coord(2,3) = ',nucl_coord(2,3)
dist = dabs(nucl_coord(1,3) - nucl_coord(2,3))
zmax += dist
zcenter = (nucl_coord(1,3) + nucl_coord(2,3))*0.5d0
print*,'zcenter = ',zcenter
print*,'zmax = ',zmax
nz = 1000
dz = zmax / dble(nz)
r(:) = 0.d0
r(3) = zcenter -zmax * 0.5d0
print*,'r(3) = ',r(3)
istate = 1
do i = 1, nz
call give_on_top_in_r_one_state(r,istate,on_top_in_r)
call give_cas_density_in_r(r,total_density)
write(i_unit_output,*)r(3),on_top_in_r,total_density
r(3) += dz
enddo
end

76
src/cas_based_on_top/on_top_cas_prov.irp.f Normal file
View File

@ -0,0 +1,76 @@
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

118
src/cas_based_on_top/on_top_cas_rout.irp.f Normal file
View File

@ -0,0 +1,118 @@
subroutine give_core_inact_act_density_in_r(r,mos_array,core_density_in_r,inact_density_in_r,act_density_in_r, total_density)
implicit none
double precision, intent(in) :: r(3),mos_array(mo_num)
double precision, intent(out):: core_density_in_r,inact_density_in_r,act_density_in_r(2,N_states),total_density(N_states)
BEGIN_DOC
! core, inactive and active part of the density for alpha/beta electrons
!
! the density coming from the core and inactive are the same for alpha/beta electrons
!
! act_density(1/2, i) = alpha/beta density for the ith state
!
! total_density(i) = 2 * (core_density_in_r+inact_density_in_r) + act_density_in_r(1,i) + act_density_in_r(2,i)
END_DOC
integer :: i,j,istate
core_density_in_r = 0.d0
do i = 1, n_core_orb
j = list_core(i)
core_density_in_r += mos_array(j) * mos_array(j)
enddo
inact_density_in_r = 0.d0
do i = 1, n_inact_orb
j = list_inact(i)
inact_density_in_r += mos_array(j) * mos_array(j)
enddo
double precision, allocatable :: act_mos(:)
double precision :: tmp
allocate(act_mos(n_act_orb))
do i = 1, n_act_orb
j = list_act(i)
act_mos(i) = mos_array(j)
enddo
act_density_in_r = 0.d0
do istate = 1, N_states
do i = 1, n_act_orb
do j = 1, n_act_orb
tmp = act_mos(i) * act_mos(j)
act_density_in_r(1,istate) += tmp * one_e_act_dm_alpha_mo_for_dft(j,i,istate)
act_density_in_r(2,istate) += tmp * one_e_act_dm_beta_mo_for_dft(j,i,istate)
enddo
enddo
total_density(istate) = 2.d0 * (core_density_in_r + inact_density_in_r) + act_density_in_r(1,istate) + act_density_in_r(2,istate)
enddo
end
subroutine give_active_on_top_in_r_one_state(r,istate,mos_array,act_on_top)
implicit none
BEGIN_DOC
! gives the purely active on-top pair density for a given state
END_DOC
integer, intent(in) :: istate
double precision, intent(in) :: r(3),mos_array(mo_num)
double precision, intent(out) :: act_on_top
double precision :: phi_i,phi_j,phi_k,phi_l
integer :: i,j,k,l
double precision, allocatable :: act_mos(:)
double precision :: tmp
allocate(act_mos(n_act_orb))
do i = 1, n_act_orb
j = list_act(i)
act_mos(i) = mos_array(j)
enddo
act_on_top = 0.d0
do l = 1, n_act_orb
phi_l = act_mos(l)
do k = 1, n_act_orb
phi_k = act_mos(k)
do j = 1, n_act_orb
phi_j = act_mos(j)
tmp = phi_l * phi_k * phi_j
do i = 1, n_act_orb
phi_i = act_mos(i)
! 1 2 1 2
act_on_top += act_2_rdm_ab_mo(i,j,k,l,istate) * tmp * phi_i
enddo
enddo
enddo
enddo
end
subroutine give_on_top_in_r_one_state(r,istate,on_top_in_r)
implicit none
integer, intent(in) :: istate
double precision, intent(in) :: r(3)
double precision, intent(out) :: on_top_in_r
BEGIN_DOC
! on top pair density in r for the state istate a CAS-BASED wf
!
! note that if no_core_density .EQ. .True., all core contributions are excluded
END_DOC
double precision, allocatable :: mos_array(:)
provide act_2_rdm_ab_mo one_e_act_dm_alpha_mo_for_dft one_e_act_dm_beta_mo_for_dft
allocate(mos_array(mo_num))
call give_all_mos_at_r(r,mos_array)
double precision :: core_density_in_r, inact_density_in_r, act_density_in_r(2,N_states), total_density(N_states)
double precision :: act_on_top,core_inact_dm
! getting the different part of the density in r
call give_core_inact_act_density_in_r(r,mos_array,core_density_in_r,inact_density_in_r,act_density_in_r, total_density)
! getting the purely active part of the density in r
call give_active_on_top_in_r_one_state(r,istate,mos_array,act_on_top)
if(no_core_density) then
core_inact_dm = inact_density_in_r
else
core_inact_dm = core_density_in_r + inact_density_in_r
endif
on_top_in_r = act_on_top + core_inact_dm * (act_density_in_r(1,istate) + act_density_in_r(2,istate)) + core_inact_dm*core_inact_dm
end

73
src/cas_based_on_top/on_top_grad.irp.f Normal file
View File

@ -0,0 +1,73 @@
subroutine give_on_top_gradient(ipoint,istate,ontop_grad)
implicit none
BEGIN_DOC
! on top pair density and its gradient evaluated at a given point of the grid
! ontop_grad(1:3) :: gradients of the on-top pair density
! ontop_grad(4) :: on-top pair density
END_DOC
double precision, intent(out) :: ontop_grad(4)
integer, intent(in) :: ipoint,istate
double precision :: phi_jkl,phi_ikl,phi_ijl,phi_ijk
integer :: i,j,k,l,m
ontop_grad = 0.d0
do l = 1, n_core_inact_act_orb
do k = 1, n_core_inact_act_orb
do j = 1, n_core_inact_act_orb
do i = 1, n_core_inact_act_orb
phi_jkl = core_inact_act_mos_in_r_array(j,ipoint) * core_inact_act_mos_in_r_array(k,ipoint) * core_inact_act_mos_in_r_array(l,ipoint)
phi_ikl = core_inact_act_mos_in_r_array(i,ipoint) * core_inact_act_mos_in_r_array(k,ipoint) * core_inact_act_mos_in_r_array(l,ipoint)
phi_ijl = core_inact_act_mos_in_r_array(i,ipoint) * core_inact_act_mos_in_r_array(j,ipoint) * core_inact_act_mos_in_r_array(l,ipoint)
phi_ijk = core_inact_act_mos_in_r_array(i,ipoint) * core_inact_act_mos_in_r_array(j,ipoint) * core_inact_act_mos_in_r_array(k,ipoint)
! 1 2 1 2
ontop_grad(4) += phi_ijk * core_inact_act_mos_in_r_array(l,ipoint) * full_occ_2_rdm_ab_mo(i,j,k,l,istate)
do m = 1,3
ontop_grad (m) += full_occ_2_rdm_ab_mo(i,j,k,l,istate) * &
( core_inact_act_mos_grad_in_r_array(m,i,ipoint) * phi_jkl + core_inact_act_mos_grad_in_r_array(m,j,ipoint) * phi_ikl + &
core_inact_act_mos_grad_in_r_array(m,k,ipoint) * phi_ijl + core_inact_act_mos_grad_in_r_array(m,l,ipoint) * phi_ijk )
enddo
enddo
enddo
enddo
enddo
end
BEGIN_PROVIDER [double precision, grad_total_cas_on_top_density,(4,n_points_final_grid,N_states) ]
&BEGIN_PROVIDER [double precision, wall_time_core_inact_act_on_top_of_r ]
implicit none
BEGIN_DOC
! grad_total_cas_on_top_density(1:3,ipoint,istate) : provider for the on top pair density gradient (x,y,z) for the point 'ipoint' and state 'istate'
!
! grad_total_cas_on_top_density(4,ipoint,istate) : on top pair density for the point 'ipoint' and state 'istate'
END_DOC
integer :: i_point,i_state,i
double precision :: wall_0,wall_1
double precision :: core_inact_act_on_top_of_r_from_provider,ontop_grad(4)
print*,'providing the core_inact_act_on_top_of_r'
i_point = 1
provide full_occ_2_rdm_ab_mo
i_state = 1
call give_on_top_gradient(i_point,i_state,ontop_grad)
call wall_time(wall_0)
!$OMP PARALLEL DO &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i_point,i_state,ontop_grad) &
!$OMP SHARED(grad_total_cas_on_top_density,n_points_final_grid,N_states)
do i_point = 1, n_points_final_grid
do i_state = 1, N_states
call give_on_top_gradient(i_point,i_state,ontop_grad)
do i = 1, 4
grad_total_cas_on_top_density(i,i_point,i_state) = ontop_grad(i)
enddo
enddo
enddo
!$OMP END PARALLEL DO
call wall_time(wall_1)
print*,'provided the core_inact_act_on_top_of_r'
print*,'Time to provide :',wall_1 - wall_0
wall_time_core_inact_act_on_top_of_r = wall_1 - wall_0
END_PROVIDER