LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-07-03 09:55:59 +02:00
Code Issues Releases Wiki Activity

Modifs of fobo-scf

Browse Source
This commit is contained in:
Emmanuel Giner 2016-07-16 16:09:50 +02:00
parent e50cfeee02
commit a0d5869054
39 changed files with 4424 additions and 369 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
config/ifort.cfg
View File

@ -18,7 +18,7 @@ IRPF90_FLAGS : --ninja --align=32
# 0 : Deactivate # 0 : Deactivate
# #
[OPTION] [OPTION]
MODE : DEBUG ; [ OPT | PROFILE | DEBUG ] : Chooses the section below MODE : OPT ; [ OPT | PROFILE | DEBUG ] : Chooses the section below
CACHE : 1 ; Enable cache_compile.py CACHE : 1 ; Enable cache_compile.py
OPENMP : 1 ; Append OpenMP flags OPENMP : 1 ; Append OpenMP flags

7
plugins/All_singles/H_apply.irp.f
View File

@ -8,6 +8,13 @@ s.unset_skip()
s.filter_only_1h1p() s.filter_only_1h1p()
print s print s
s = H_apply("just_1h_1p_singles",do_double_exc=False)
s.set_selection_pt2("epstein_nesbet_2x2")
s.unset_skip()
s.filter_only_1h1p()
print s
s = H_apply("just_mono",do_double_exc=False) s = H_apply("just_mono",do_double_exc=False)
s.set_selection_pt2("epstein_nesbet_2x2") s.set_selection_pt2("epstein_nesbet_2x2")
s.unset_skip() s.unset_skip()

2
plugins/All_singles/all_1h_1p.irp.f
View File

@ -49,7 +49,7 @@ subroutine routine
endif endif
call save_wavefunction call save_wavefunction
if(n_det_before == N_det)then if(n_det_before == N_det)then
selection_criterion = selection_criterion * 0.5d0 selection_criterion_factor = selection_criterion_factor * 0.5d0
endif endif
enddo enddo

76
plugins/All_singles/all_1h_1p_singles.irp.f Normal file
View File

@ -0,0 +1,76 @@
program restart_more_singles
BEGIN_DOC
! Generates and select single and double excitations of type 1h-1p
! on the top of a given restart wave function of type CAS
END_DOC
read_wf = .true.
touch read_wf
print*,'ref_bitmask_energy = ',ref_bitmask_energy
call routine
end
subroutine routine
implicit none
integer :: i,k
double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:),E_before(:)
integer :: N_st, degree
integer :: n_det_before
N_st = N_states
allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st),E_before(N_st))
i = 0
print*,'N_det = ',N_det
print*,'n_det_max = ',n_det_max
print*,'pt2_max = ',pt2_max
pt2=-1.d0
E_before = ref_bitmask_energy
do while (N_det < n_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max)
n_det_before = N_det
i += 1
print*,'-----------------------'
print*,'i = ',i
call H_apply_just_1h_1p_singles(pt2, norm_pert, H_pert_diag, N_st)
call diagonalize_CI
print*,'N_det = ',N_det
print*,'E = ',CI_energy(1)
print*,'pt2 = ',pt2(1)
print*,'E+PT2 = ',E_before + pt2(1)
E_before = CI_energy
if(N_states_diag.gt.1)then
print*,'Variational Energy difference'
do i = 2, N_st
print*,'Delta E = ',CI_energy(i) - CI_energy(1)
enddo
endif
if(N_states.gt.1)then
print*,'Variational + perturbative Energy difference'
do i = 2, N_st
print*,'Delta E = ',E_before(i)+ pt2(i) - (E_before(1) + pt2(1))
enddo
endif
call save_wavefunction
if(n_det_before == N_det)then
selection_criterion_factor = selection_criterion_factor * 0.5d0
endif
enddo
threshold_davidson = 1.d-10
soft_touch threshold_davidson davidson_criterion
call diagonalize_CI
if(N_states_diag.gt.1)then
print*,'Variational Energy difference'
do i = 2, N_st
print*,'Delta E = ',CI_energy(i) - CI_energy(1)
enddo
endif
if(N_states.gt.1)then
print*,'Variational + perturbative Energy difference'
do i = 2, N_st
print*,'Delta E = ',CI_energy(i)+ pt2(i) - (CI_energy(1) + pt2(1))
enddo
endif
call ezfio_set_all_singles_energy(CI_energy)
call save_wavefunction
deallocate(pt2,norm_pert)
end

31
plugins/DFT_Utils/grid_density.irp.f
View File

@ -1,11 +1,11 @@
BEGIN_PROVIDER [integer, n_points_angular_grid] BEGIN_PROVIDER [integer, n_points_angular_grid]
implicit none implicit none
n_points_angular_grid = 51 n_points_angular_grid = 50
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [integer, n_points_radial_grid] BEGIN_PROVIDER [integer, n_points_radial_grid]
implicit none implicit none
n_points_radial_grid = 1000 n_points_radial_grid = 10000
END_PROVIDER END_PROVIDER
@ -15,8 +15,28 @@ END_PROVIDER
BEGIN_DOC BEGIN_DOC
! weights and grid points for the integration on the angular variables on ! weights and grid points for the integration on the angular variables on
! the unit sphere centered on (0,0,0) ! the unit sphere centered on (0,0,0)
! According to the LEBEDEV scheme
END_DOC END_DOC
call cal_quad(n_points_angular_grid, angular_quadrature_points,weights_angular_points) call cal_quad(n_points_angular_grid, angular_quadrature_points,weights_angular_points)
include 'constants.include.F'
integer :: i
double precision :: accu
double precision :: degre_rad
!degre_rad = 180.d0/pi
!accu = 0.d0
!do i = 1, n_points_integration_angular_lebedev
! accu += weights_angular_integration_lebedev(i)
! weights_angular_points(i) = weights_angular_integration_lebedev(i) * 2.d0 * pi
! angular_quadrature_points(i,1) = dcos ( degre_rad * theta_angular_integration_lebedev(i)) &
! * dsin ( degre_rad * phi_angular_integration_lebedev(i))
! angular_quadrature_points(i,2) = dsin ( degre_rad * theta_angular_integration_lebedev(i)) &
! * dsin ( degre_rad * phi_angular_integration_lebedev(i))
! angular_quadrature_points(i,3) = dcos ( degre_rad * phi_angular_integration_lebedev(i))
!enddo
!print*,'ANGULAR'
!print*,''
!print*,'accu = ',accu
!ASSERT( dabs(accu - 1.D0) < 1.d-10)
END_PROVIDER END_PROVIDER
@ -55,7 +75,7 @@ BEGIN_PROVIDER [double precision, grid_points_per_atom, (3,n_points_angular_grid
x_ref = nucl_coord(i,1) x_ref = nucl_coord(i,1)
y_ref = nucl_coord(i,2) y_ref = nucl_coord(i,2)
z_ref = nucl_coord(i,3) z_ref = nucl_coord(i,3)
do j = 1, n_points_radial_grid do j = 1, n_points_radial_grid-1
double precision :: x,r double precision :: x,r
x = grid_points_radial(j) ! x value for the mapping of the [0, +\infty] to [0,1] x = grid_points_radial(j) ! x value for the mapping of the [0, +\infty] to [0,1]
r = knowles_function(alpha_knowles(int(nucl_charge(i))),m_knowles,x) ! value of the radial coordinate for the integration r = knowles_function(alpha_knowles(int(nucl_charge(i))),m_knowles,x) ! value of the radial coordinate for the integration
@ -81,7 +101,7 @@ BEGIN_PROVIDER [double precision, weight_functions_at_grid_points, (n_points_ang
double precision :: tmp_array(nucl_num) double precision :: tmp_array(nucl_num)
! run over all points in space ! run over all points in space
do j = 1, nucl_num ! that are referred to each atom do j = 1, nucl_num ! that are referred to each atom
do k = 1, n_points_radial_grid !for each radial grid attached to the "jth" atom do k = 1, n_points_radial_grid -1 !for each radial grid attached to the "jth" atom
do l = 1, n_points_angular_grid ! for each angular point attached to the "jth" atom do l = 1, n_points_angular_grid ! for each angular point attached to the "jth" atom
r(1) = grid_points_per_atom(1,l,k,j) r(1) = grid_points_per_atom(1,l,k,j)
r(2) = grid_points_per_atom(2,l,k,j) r(2) = grid_points_per_atom(2,l,k,j)
@ -95,6 +115,7 @@ BEGIN_PROVIDER [double precision, weight_functions_at_grid_points, (n_points_ang
enddo enddo
accu = 1.d0/accu accu = 1.d0/accu
weight_functions_at_grid_points(l,k,j) = tmp_array(j) * accu weight_functions_at_grid_points(l,k,j) = tmp_array(j) * accu
! print*,weight_functions_at_grid_points(l,k,j)
enddo enddo
enddo enddo
enddo enddo
@ -110,7 +131,7 @@ END_PROVIDER
double precision :: r(3) double precision :: r(3)
double precision :: aos_array(ao_num),mos_array(mo_tot_num) double precision :: aos_array(ao_num),mos_array(mo_tot_num)
do j = 1, nucl_num do j = 1, nucl_num
do k = 1, n_points_radial_grid do k = 1, n_points_radial_grid -1
do l = 1, n_points_angular_grid do l = 1, n_points_angular_grid
one_body_dm_mo_alpha_at_grid_points(l,k,j) = 0.d0 one_body_dm_mo_alpha_at_grid_points(l,k,j) = 0.d0
one_body_dm_mo_beta_at_grid_points(l,k,j) = 0.d0 one_body_dm_mo_beta_at_grid_points(l,k,j) = 0.d0

26
plugins/DFT_Utils/integration_3d.irp.f
View File

@ -4,12 +4,18 @@ double precision function step_function_becke(x)
double precision :: f_function_becke double precision :: f_function_becke
integer :: i,n_max_becke integer :: i,n_max_becke
!if(x.lt.-1.d0)then
! step_function_becke = 0.d0
!else if (x .gt.1)then
! step_function_becke = 0.d0
!else
step_function_becke = f_function_becke(x) step_function_becke = f_function_becke(x)
n_max_becke = 3 !!n_max_becke = 1
do i = 1, n_max_becke do i = 1, 4
step_function_becke = f_function_becke(step_function_becke) step_function_becke = f_function_becke(step_function_becke)
enddo enddo
step_function_becke = 0.5d0*(1.d0 - step_function_becke) step_function_becke = 0.5d0*(1.d0 - step_function_becke)
!endif
end end
double precision function f_function_becke(x) double precision function f_function_becke(x)
@ -46,19 +52,3 @@ double precision function cell_function_becke(r,atom_number)
enddo enddo
end end
double precision function weight_function_becke(r,atom_number)
implicit none
double precision, intent(in) :: r(3)
integer, intent(in) :: atom_number
BEGIN_DOC
! atom_number :: atom on which the weight function of Becke (1988, JCP,88(4))
! r(1:3) :: x,y,z coordinantes of the current point
END_DOC
double precision :: cell_function_becke,accu
integer :: j
accu = 0.d0
do j = 1, nucl_num
accu += cell_function_becke(r,j)
enddo
weight_function_becke = cell_function_becke(r,atom_number)/accu
end

3
plugins/DFT_Utils/integration_radial.irp.f
View File

@ -16,7 +16,7 @@
do j = 1, nucl_num do j = 1, nucl_num
integral_density_alpha_knowles_becke_per_atom(j) = 0.d0 integral_density_alpha_knowles_becke_per_atom(j) = 0.d0
integral_density_beta_knowles_becke_per_atom(j) = 0.d0 integral_density_beta_knowles_becke_per_atom(j) = 0.d0
do i = 1, n_points_radial_grid do i = 1, n_points_radial_grid-1
! Angular integration over the solid angle Omega for a FIXED angular coordinate "r" ! Angular integration over the solid angle Omega for a FIXED angular coordinate "r"
f_average_angular_alpha = 0.d0 f_average_angular_alpha = 0.d0
f_average_angular_beta = 0.d0 f_average_angular_beta = 0.d0
@ -48,7 +48,6 @@ END_PROVIDER
double precision, intent(in) :: alpha,x double precision, intent(in) :: alpha,x
integer, intent(in) :: m integer, intent(in) :: m
knowles_function = -alpha * dlog(1.d0-x**m) knowles_function = -alpha * dlog(1.d0-x**m)
!knowles_function = 1.d0
end end
double precision function derivative_knowles_function(alpha,m,x) double precision function derivative_knowles_function(alpha,m,x)

699
plugins/FOBOCI/SC2_1h1p.irp.f Normal file
View File

@ -0,0 +1,699 @@
subroutine dressing_1h1p(dets_in,u_in,diag_H_elements,dim_in,sze,N_st,Nint,convergence)
use bitmasks
implicit none
BEGIN_DOC
! CISD+SC2 method :: take off all the disconnected terms of a ROHF+1h1p (selected or not)
!
! dets_in : bitmasks corresponding to determinants
!
! u_in : guess coefficients on the various states. Overwritten
! on exit
!
! dim_in : leftmost dimension of u_in
!
! sze : Number of determinants
!
! N_st : Number of eigenstates
!
! Initial guess vectors are not necessarily orthonormal
END_DOC
integer, intent(in) :: dim_in, sze, N_st, Nint
integer(bit_kind), intent(in) :: dets_in(Nint,2,sze)
double precision, intent(inout) :: u_in(dim_in,N_st)
double precision, intent(out) :: diag_H_elements(dim_in)
double precision, intent(in) :: convergence
integer :: i,j,k,l
integer :: n_singles
integer :: index_singles(sze),hole_particles_singles(sze,3)
integer :: n_doubles
integer :: index_doubles(sze),hole_particles_doubles(sze,2)
integer :: index_hf
double precision :: e_corr_singles(mo_tot_num,2)
double precision :: e_corr_doubles(mo_tot_num)
double precision :: e_corr_singles_total(2)
double precision :: e_corr_doubles_1h1p
integer :: exc(0:2,2,2),degree
integer :: h1,h2,p1,p2,s1,s2
integer :: other_spin(2)
double precision :: phase
integer(bit_kind) :: key_tmp(N_int,2)
integer :: i_ok
double precision :: phase_single_double,phase_double_hf,get_mo_bielec_integral_schwartz
double precision :: hij,c_ref,contrib
integer :: iorb
other_spin(1) = 2
other_spin(2) = 1
n_singles = 0
n_doubles = 0
do i = 1,sze
call get_excitation(ref_bitmask,dets_in(1,1,i),exc,degree,phase,N_int)
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
call i_H_j(dets_in(1,1,i),dets_in(1,1,i),N_int,hij)
diag_H_elements(i) = hij
if(degree == 0)then
index_hf = i
else if (degree == 1)then
n_singles +=1
index_singles(n_singles) = i
! h1 = inactive orbital of the hole
hole_particles_singles(n_singles,1) = h1
! p1 = virtual orbital of the particle
hole_particles_singles(n_singles,2) = p1
! s1 = spin of the electron excited
hole_particles_singles(n_singles,3) = s1
else if (degree == 2)then
n_doubles +=1
index_doubles(n_doubles) = i
! h1 = inactive orbital of the hole (beta of course)
hole_particles_doubles(n_doubles,1) = h1
! p1 = virtual orbital of the particle (alpha of course)
hole_particles_doubles(n_doubles,2) = p2
else
print*,'PB !! found out other thing than a single or double'
print*,'stopping ..'
stop
endif
enddo
e_corr_singles = 0.d0
e_corr_doubles = 0.d0
e_corr_singles_total = 0.d0
e_corr_doubles_1h1p = 0.d0
c_ref = 1.d0/u_in(index_hf,1)
print*,'c_ref = ',c_ref
do i = 1,sze
call get_excitation(ref_bitmask,dets_in(1,1,i),exc,degree,phase,N_int)
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
call i_H_j(ref_bitmask,dets_in(1,1,i),N_int,hij)
contrib = hij * u_in(i,1) * c_ref
if (degree == 1)then
e_corr_singles(h1,s1) += contrib
e_corr_singles(p1,s1) += contrib
e_corr_singles_total(s1)+= contrib
else if (degree == 2)then
e_corr_doubles_1h1p += contrib
e_corr_doubles(h1) += contrib
e_corr_doubles(p2) += contrib
endif
enddo
print*,'e_corr_singles alpha = ',e_corr_singles_total(1)
print*,'e_corr_singles beta = ',e_corr_singles_total(2)
print*,'e_corr_doubles_1h1p = ',e_corr_doubles_1h1p
! repeat all the correlation energy on the singles
do i = 1,n_singles
! you can repeat all the correlation energy of the single excitation of the other spin
diag_H_elements(index_singles(i)) += e_corr_singles_total(other_spin(hole_particles_singles(i,3)))
! you can repeat all the correlation energy of the single excitation of the same spin
do j = 1, n_inact_orb
iorb = list_inact(j)
! except the one of the hole
if(iorb == hole_particles_singles(i,1))cycle
! ispin = hole_particles_singles(i,3)
diag_H_elements(index_singles(i)) += e_corr_singles(iorb,hole_particles_singles(i,3))
enddo
! also exclude all the energy coming from the virtual orbital
diag_H_elements(index_singles(i)) -= e_corr_singles(hole_particles_singles(i,2),hole_particles_singles(i,3))
! If it is a single excitation alpha, you can repeat :
! +) all the double excitation 1h1p, appart the part involving the virtual orbital "r"
! If it is a single excitation alpha, you can repeat :
! +) all the double excitation 1h1p, appart the part involving the inactive orbital "i"
diag_H_elements(index_singles(i)) += e_corr_doubles_1h1p
if(hole_particles_singles(i,3) == 1)then ! alpha single excitation
diag_H_elements(index_singles(i)) -= e_corr_doubles(hole_particles_singles(i,2))
else ! beta single exctitation
diag_H_elements(index_singles(i)) -= e_corr_doubles(hole_particles_singles(i,1))
endif
enddo
! repeat all the correlation energy on the doubles
! as all the doubles involve the active space, you cannot repeat any of them one on another
do i = 1, n_doubles
! on a given double, you can repeat all the correlation energy of the singles alpha
do j = 1, n_inact_orb
iorb = list_inact(j)
! ispin = hole_particles_singles(i,3)
diag_H_elements(index_doubles(i)) += e_corr_singles(iorb,1)
enddo
! except the part involving the virtual orbital "hole_particles_doubles(i,2)"
diag_H_elements(index_doubles(i)) -= e_corr_singles(hole_particles_doubles(i,2),1)
! on a given double, you can repeat all the correlation energy of the singles beta
do j = 1, n_inact_orb
iorb = list_inact(j)
! except the one of the hole
if(iorb == hole_particles_doubles(i,1))cycle
! ispin = hole_particles_singles(i,3)
diag_H_elements(index_doubles(i)) += e_corr_singles(iorb,2)
enddo
enddo
! Taking into account the connected part of the 2h2p on the HF determinant
! 1/2 \sum_{ir,js} c_{ir}^{sigma} c_{js}^{sigma}
! diag_H_elements(index_hf) += total_corr_e_2h2p
c_ref = c_ref * c_ref
print*,'diag_H_elements(index_hf) = ',diag_H_elements(index_hf)
do i = 1, n_singles
! start on the single excitation "|i>"
h1 = hole_particles_singles(i,1)
p1 = hole_particles_singles(i,2)
do j = 1, n_singles
do k = 1, N_int
key_tmp(k,1) = dets_in(k,1,index_singles(i))
key_tmp(k,2) = dets_in(k,2,index_singles(i))
enddo
h2 = hole_particles_singles(j,1)
p2 = hole_particles_singles(j,2)
call do_mono_excitation(key_tmp,h2,p2,hole_particles_singles(j,3),i_ok)
! apply the excitation operator from the single excitation "|j>"
if(i_ok .ne. 1)cycle
double precision :: phase_ref_other_single,diag_H_mat_elem,hijj,contrib_e2,coef_1
call get_excitation(key_tmp,dets_in(1,1,index_singles(i)),exc,degree,phase_single_double,N_int)
call get_excitation(ref_bitmask,dets_in(1,1,index_singles(j)),exc,degree,phase_ref_other_single,N_int)
call i_H_j(ref_bitmask,key_tmp,N_int,hij)
diag_H_elements(index_hf) += u_in(index_singles(i),1) * u_in(index_singles(j),1) * c_ref * hij &
* phase_single_double * phase_ref_other_single
enddo
enddo
print*,'diag_H_elements(index_hf) = ',diag_H_elements(index_hf)
end
subroutine dressing_1h1p_full(dets_in,u_in,H_matrix,dim_in,sze,N_st,Nint,convergence)
use bitmasks
implicit none
BEGIN_DOC
! CISD+SC2 method :: take off all the disconnected terms of a ROHF+1h1p (selected or not)
!
! dets_in : bitmasks corresponding to determinants
!
! u_in : guess coefficients on the various states. Overwritten
! on exit
!
! dim_in : leftmost dimension of u_in
!
! sze : Number of determinants
!
! N_st : Number of eigenstates
!
! Initial guess vectors are not necessarily orthonormal
END_DOC
integer, intent(in) :: dim_in, sze, N_st, Nint
integer(bit_kind), intent(in) :: dets_in(Nint,2,sze)
double precision, intent(in) :: u_in(dim_in,N_st)
double precision, intent(inout) :: H_matrix(sze,sze)
double precision, intent(in) :: convergence
integer :: i,j,k,l
integer :: n_singles
integer :: index_singles(sze),hole_particles_singles(sze,3)
integer :: n_doubles
integer :: index_doubles(sze),hole_particles_doubles(sze,2)
integer :: index_hf
double precision :: e_corr_singles(mo_tot_num,2)
double precision :: e_corr_doubles(mo_tot_num)
double precision :: e_corr_singles_total(2)
double precision :: e_corr_doubles_1h1p
integer :: exc(0:2,2,2),degree
integer :: h1,h2,p1,p2,s1,s2
integer :: other_spin(2)
double precision :: phase
integer(bit_kind) :: key_tmp(N_int,2)
integer :: i_ok
double precision :: phase_single_double,phase_double_hf,get_mo_bielec_integral_schwartz
double precision :: hij,c_ref,contrib
integer :: iorb
other_spin(1) = 2
other_spin(2) = 1
n_singles = 0
n_doubles = 0
do i = 1,sze
call get_excitation(ref_bitmask,dets_in(1,1,i),exc,degree,phase,N_int)
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
if(degree == 0)then
index_hf = i
else if (degree == 1)then
n_singles +=1
index_singles(n_singles) = i
! h1 = inactive orbital of the hole
hole_particles_singles(n_singles,1) = h1
! p1 = virtual orbital of the particle
hole_particles_singles(n_singles,2) = p1
! s1 = spin of the electron excited
hole_particles_singles(n_singles,3) = s1
else if (degree == 2)then
n_doubles +=1
index_doubles(n_doubles) = i
! h1 = inactive orbital of the hole (beta of course)
hole_particles_doubles(n_doubles,1) = h1
! p1 = virtual orbital of the particle (alpha of course)
hole_particles_doubles(n_doubles,2) = p2
else
print*,'PB !! found out other thing than a single or double'
print*,'stopping ..'
stop
endif
enddo
double precision, allocatable :: dressing_H_mat_elem(:)
allocate(dressing_H_mat_elem(N_det))
logical :: lmct
dressing_H_mat_elem = 0.d0
call dress_diag_elem_2h2p(dressing_H_mat_elem,N_det)
lmct = .False.
call dress_diag_elem_2h1p(dressing_H_mat_elem,N_det,lmct,1000)
lmct = .true.
call dress_diag_elem_1h2p(dressing_H_mat_elem,N_det,lmct,1000)
do i = 1, N_det
H_matrix(i,i) += dressing_H_mat_elem(i)
enddo
e_corr_singles = 0.d0
e_corr_doubles = 0.d0
e_corr_singles_total = 0.d0
e_corr_doubles_1h1p = 0.d0
c_ref = 1.d0/u_in(index_hf,1)
print*,'c_ref = ',c_ref
do i = 1,sze
call get_excitation(ref_bitmask,dets_in(1,1,i),exc,degree,phase,N_int)
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
call i_H_j(ref_bitmask,dets_in(1,1,i),N_int,hij)
contrib = hij * u_in(i,1) * c_ref
if (degree == 1)then
e_corr_singles(h1,s1) += contrib
e_corr_singles(p1,s1) += contrib
e_corr_singles_total(s1)+= contrib
else if (degree == 2)then
e_corr_doubles_1h1p += contrib
e_corr_doubles(h1) += contrib
e_corr_doubles(p2) += contrib
endif
enddo
print*,'e_corr_singles alpha = ',e_corr_singles_total(1)
print*,'e_corr_singles beta = ',e_corr_singles_total(2)
print*,'e_corr_doubles_1h1p = ',e_corr_doubles_1h1p
! repeat all the correlation energy on the singles
! do i = 1,n_singles
! ! you can repeat all the correlation energy of the single excitation of the other spin
! H_matrix(index_singles(i),index_singles(i)) += e_corr_singles_total(other_spin(hole_particles_singles(i,3)))
! ! you can repeat all the correlation energy of the single excitation of the same spin
! do j = 1, n_inact_orb
! iorb = list_inact(j)
! ! except the one of the hole
! if(iorb == hole_particles_singles(i,1))cycle
! ! ispin = hole_particles_singles(i,3)
! H_matrix(index_singles(i),index_singles(i)) += e_corr_singles(iorb,hole_particles_singles(i,3))
! enddo
! ! also exclude all the energy coming from the virtual orbital
! H_matrix(index_singles(i),index_singles(i)) -= e_corr_singles(hole_particles_singles(i,2),hole_particles_singles(i,3))
!
! ! If it is a single excitation alpha, you can repeat :
! ! +) all the double excitation 1h1p, appart the part involving the virtual orbital "r"
! ! If it is a single excitation alpha, you can repeat :
! ! +) all the double excitation 1h1p, appart the part involving the inactive orbital "i"
! H_matrix(index_singles(i),index_singles(i)) += e_corr_doubles_1h1p
! if(hole_particles_singles(i,3) == 1)then ! alpha single excitation
! H_matrix(index_singles(i),index_singles(i)) -= e_corr_doubles(hole_particles_singles(i,2))
! else ! beta single exctitation
! H_matrix(index_singles(i),index_singles(i)) -= e_corr_doubles(hole_particles_singles(i,1))
! endif
! enddo
! ! repeat all the correlation energy on the doubles
! ! as all the doubles involve the active space, you cannot repeat any of them one on another
! do i = 1, n_doubles
! ! on a given double, you can repeat all the correlation energy of the singles alpha
! do j = 1, n_inact_orb
! iorb = list_inact(j)
! ! ispin = hole_particles_singles(i,3)
! H_matrix(index_doubles(i),index_doubles(i)) += e_corr_singles(iorb,1)
! enddo
! ! except the part involving the virtual orbital "hole_particles_doubles(i,2)"
! H_matrix(index_doubles(i),index_doubles(i)) -= e_corr_singles(hole_particles_doubles(i,2),1)
! ! on a given double, you can repeat all the correlation energy of the singles beta
! do j = 1, n_inact_orb
! iorb = list_inact(j)
! ! except the one of the hole
! if(iorb == hole_particles_doubles(i,1))cycle
! ! ispin = hole_particles_singles(i,3)
! H_matrix(index_doubles(i),index_doubles(i)) += e_corr_singles(iorb,2)
! enddo
! enddo
! Taking into account the connected part of the 2h2p on the HF determinant
! 1/2 \sum_{ir,js} c_{ir}^{sigma} c_{js}^{sigma}
! H_matrix(index_hf) += total_corr_e_2h2p
print*,'H_matrix(index_hf,index_hf) = ',H_matrix(index_hf,index_hf)
do i = 1, n_singles
! start on the single excitation "|i>"
h1 = hole_particles_singles(i,1)
p1 = hole_particles_singles(i,2)
print*,'i = ',i
do j = i+1, n_singles
do k = 1, N_int
key_tmp(k,1) = dets_in(k,1,index_singles(i))
key_tmp(k,2) = dets_in(k,2,index_singles(i))
enddo
h2 = hole_particles_singles(j,1)
p2 = hole_particles_singles(j,2)
call do_mono_excitation(key_tmp,h2,p2,hole_particles_singles(j,3),i_ok)
! apply the excitation operator from the single excitation "|j>"
if(i_ok .ne. 1)cycle
double precision :: H_array(sze),diag_H_mat_elem,hjj
do k = 1, sze
call get_excitation_degree(dets_in(1,1,k),key_tmp,degree,N_int)
H_array(k) = 0.d0
if(degree > 2)cycle
call i_H_j(dets_in(1,1,k),key_tmp,N_int,hij)
H_array(k) = hij
enddo
hjj = 1.d0/(ref_bitmask_energy - diag_H_mat_elem(key_tmp,N_int))
! contrib_e2 = 0.5d0 * (delta_e + dsqrt(delta_e * delta_e + 4.d0 * hij * hij))
do l = 2, sze
! pause
H_matrix(l,l) += H_array(l) * H_array(l) * hjj
! H_matrix(1,l) += H_array(1) * H_array(l) * hjj
! H_matrix(l,1) += H_array(1) * H_array(l) * hjj
enddo
enddo
enddo
print*,'H_matrix(index_hf,index_hf) = ',H_matrix(index_hf,index_hf)
end
subroutine SC2_1h1p_full(dets_in,u_in,energies,H_matrix,dim_in,sze,N_st,Nint,convergence)
use bitmasks
implicit none
BEGIN_DOC
! CISD+SC2 method :: take off all the disconnected terms of a CISD (selected or not)
!
! dets_in : bitmasks corresponding to determinants
!
! u_in : guess coefficients on the various states. Overwritten
! on exit
!
! dim_in : leftmost dimension of u_in
!
! sze : Number of determinants
!
! N_st : Number of eigenstates
!
! Initial guess vectors are not necessarily orthonormal
END_DOC
integer, intent(in) :: dim_in, sze, N_st, Nint
integer(bit_kind), intent(in) :: dets_in(Nint,2,sze)
double precision, intent(inout) :: u_in(dim_in,N_st)
double precision, intent(out) :: energies(N_st)
double precision, intent(out) :: H_matrix(sze,sze)
double precision, intent(in) :: convergence
integer :: i,j,iter
print*,'sze = ',sze
do iter = 1, 1
! if(sze<=N_det_max_jacobi)then
double precision, allocatable :: eigenvectors(:,:), eigenvalues(:),H_matrix_tmp(:,:)
allocate (H_matrix_tmp(size(H_matrix_all_dets,1),sze),eigenvalues(sze),eigenvectors(size(H_matrix_all_dets,1),sze))
H_matrix_tmp = 0.d0
call dressing_1h1p_full(dets_in,u_in,H_matrix_tmp,dim_in,sze,N_st,Nint,convergence)
do j=1,sze
do i=1,sze
H_matrix_tmp(i,j) += H_matrix_all_dets(i,j)
enddo
enddo
print*,'passed the dressing'
call lapack_diag(eigenvalues,eigenvectors, &
H_matrix_tmp,size(H_matrix_all_dets,1),sze)
do j=1,min(N_states_diag,sze)
do i=1,sze
u_in(i,j) = eigenvectors(i,j)
enddo
energies(j) = eigenvalues(j)
enddo
deallocate (H_matrix_tmp, eigenvalues, eigenvectors)
! else
! call davidson_diag_hjj(dets_in,u_in,diag_H_elements,energies,dim_in,sze,N_st,Nint,output_determinants)
! endif
print*,'E = ',energies(1) + nuclear_repulsion
enddo
end
subroutine SC2_1h1p(dets_in,u_in,energies,diag_H_elements,dim_in,sze,N_st,Nint,convergence)
use bitmasks
implicit none
BEGIN_DOC
! CISD+SC2 method :: take off all the disconnected terms of a CISD (selected or not)
!
! dets_in : bitmasks corresponding to determinants
!
! u_in : guess coefficients on the various states. Overwritten
! on exit
!
! dim_in : leftmost dimension of u_in
!
! sze : Number of determinants
!
! N_st : Number of eigenstates
!
! Initial guess vectors are not necessarily orthonormal
END_DOC
integer, intent(in) :: dim_in, sze, N_st, Nint
integer(bit_kind), intent(in) :: dets_in(Nint,2,sze)
double precision, intent(inout) :: u_in(dim_in,N_st)
double precision, intent(out) :: energies(N_st)
double precision, intent(out) :: diag_H_elements(dim_in)
double precision, intent(in) :: convergence
integer :: i,j,iter
do iter = 1, 1
call dressing_1h1p(dets_in,u_in,diag_H_elements,dim_in,sze,N_st,Nint,convergence)
if(sze<=N_det_max_jacobi)then
double precision, allocatable :: eigenvectors(:,:), eigenvalues(:),H_matrix_tmp(:,:)
allocate (H_matrix_tmp(size(H_matrix_all_dets,1),sze),eigenvalues(sze),eigenvectors(size(H_matrix_all_dets,1),sze))
do j=1,sze
do i=1,sze
H_matrix_tmp(i,j) = H_matrix_all_dets(i,j)
enddo
enddo
do i = 1,sze
H_matrix_tmp(i,i) = diag_H_elements(i)
enddo
call lapack_diag(eigenvalues,eigenvectors, &
H_matrix_tmp,size(H_matrix_all_dets,1),sze)
do j=1,min(N_states_diag,sze)
do i=1,sze
u_in(i,j) = eigenvectors(i,j)
enddo
energies(j) = eigenvalues(j)
enddo
deallocate (H_matrix_tmp, eigenvalues, eigenvectors)
else
call davidson_diag_hjj(dets_in,u_in,diag_H_elements,energies,dim_in,sze,N_st,Nint,output_determinants)
endif
print*,'E = ',energies(1) + nuclear_repulsion
enddo
end
subroutine density_matrix_1h1p(dets_in,u_in,density_matrix_alpha,density_matrix_beta,norm,dim_in,sze,N_st,Nint)
use bitmasks
implicit none
BEGIN_DOC
! CISD+SC2 method :: take off all the disconnected terms of a ROHF+1h1p (selected or not)
!
! dets_in : bitmasks corresponding to determinants
!
! u_in : guess coefficients on the various states. Overwritten
! on exit
!
! dim_in : leftmost dimension of u_in
!
! sze : Number of determinants
!
! N_st : Number of eigenstates
!
! Initial guess vectors are not necessarily orthonormal
END_DOC
integer, intent(in) :: dim_in, sze, N_st, Nint
integer(bit_kind), intent(in) :: dets_in(Nint,2,sze)
double precision, intent(inout) :: u_in(dim_in,N_st)
double precision, intent(inout) :: density_matrix_alpha(mo_tot_num_align,mo_tot_num)
double precision, intent(inout) :: density_matrix_beta(mo_tot_num_align,mo_tot_num)
double precision, intent(inout) :: norm
integer :: i,j,k,l
integer :: n_singles
integer :: index_singles(sze),hole_particles_singles(sze,3)
integer :: n_doubles
integer :: index_doubles(sze),hole_particles_doubles(sze,2)
integer :: index_hf
integer :: exc(0:2,2,2),degree
integer :: h1,h2,p1,p2,s1,s2
integer :: other_spin(2)
double precision :: phase
integer(bit_kind) :: key_tmp(N_int,2)
integer :: i_ok
double precision :: phase_single_double,phase_double_hf,get_mo_bielec_integral_schwartz
double precision :: hij,c_ref,contrib
integer :: iorb
other_spin(1) = 2
other_spin(2) = 1
n_singles = 0
n_doubles = 0
norm = 0.d0
do i = 1,sze
call get_excitation(ref_bitmask,dets_in(1,1,i),exc,degree,phase,N_int)
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
norm += u_in(i,1)* u_in(i,1)
if(degree == 0)then
index_hf = i
c_ref = 1.d0/psi_coef(i,1)
else if (degree == 1)then
n_singles +=1
index_singles(n_singles) = i
! h1 = inactive orbital of the hole
hole_particles_singles(n_singles,1) = h1
! p1 = virtual orbital of the particle
hole_particles_singles(n_singles,2) = p1
! s1 = spin of the electron excited
hole_particles_singles(n_singles,3) = s1
else if (degree == 2)then
n_doubles +=1
index_doubles(n_doubles) = i
! h1 = inactive orbital of the hole (beta of course)
hole_particles_doubles(n_doubles,1) = h1
! p1 = virtual orbital of the particle (alpha of course)
hole_particles_doubles(n_doubles,2) = p2
else
print*,'PB !! found out other thing than a single or double'
print*,'stopping ..'
stop
endif
enddo
print*,'norm = ',norm
! Taking into account the connected part of the 2h2p on the HF determinant
! 1/2 \sum_{ir,js} c_{ir}^{sigma} c_{js}^{sigma}
do i = 1, n_singles
! start on the single excitation "|i>"
h1 = hole_particles_singles(i,1)
p1 = hole_particles_singles(i,2)
do j = 1, n_singles
do k = 1, N_int
key_tmp(k,1) = dets_in(k,1,index_singles(i))
key_tmp(k,2) = dets_in(k,2,index_singles(i))
enddo
h2 = hole_particles_singles(j,1)
p2 = hole_particles_singles(j,2)
call do_mono_excitation(key_tmp,h2,p2,hole_particles_singles(j,3),i_ok)
! apply the excitation operator from the single excitation "|j>"
if(i_ok .ne. 1)cycle
double precision :: coef_ijrs,phase_other_single_ref
integer :: occ(N_int*bit_kind_size,2),n_occ(2)
call get_excitation(key_tmp,dets_in(1,1,index_singles(i)),exc,degree,phase_single_double,N_int)
call get_excitation(ref_bitmask,dets_in(1,1,index_singles(j)),exc,degree,phase_other_single_ref,N_int)
call get_excitation(key_tmp,dets_in(1,1,index_singles(j)),exc,degree,phase_other_single_ref,N_int)
coef_ijrs = u_in(index_singles(i),1) * u_in(index_singles(j),1) * c_ref * c_ref &
* phase_single_double * phase_other_single_ref
call bitstring_to_list_ab(key_tmp, occ, n_occ, N_int)
do k=1,elec_alpha_num
l = occ(k,1)
density_matrix_alpha(l,l) += coef_ijrs*coef_ijrs
enddo
do k=1,elec_beta_num
l = occ(k,1)
density_matrix_beta(l,l) += coef_ijrs*coef_ijrs
enddo
norm += coef_ijrs* coef_ijrs
if(hole_particles_singles(j,3) == 1)then ! single alpha
density_matrix_alpha(h2,p2) += coef_ijrs * phase_single_double * u_in(index_singles(i),1) * c_ref
density_matrix_alpha(p2,h2) += coef_ijrs * phase_single_double * u_in(index_singles(i),1) * c_ref
else
density_matrix_beta(h2,p2) += coef_ijrs * phase_single_double * u_in(index_singles(i),1) * c_ref
density_matrix_beta(p2,h2) += coef_ijrs * phase_single_double * u_in(index_singles(i),1) * c_ref
endif
enddo
enddo
do i = 1, n_doubles
! start on the double excitation "|i>"
h1 = hole_particles_doubles(i,1)
p1 = hole_particles_doubles(i,2)
do j = 1, n_singles
do k = 1, N_int
key_tmp(k,1) = dets_in(k,1,index_doubles(i))
key_tmp(k,2) = dets_in(k,2,index_doubles(i))
enddo
h2 = hole_particles_singles(j,1)
p2 = hole_particles_singles(j,2)
call do_mono_excitation(key_tmp,h2,p2,hole_particles_singles(j,3),i_ok)
! apply the excitation operator from the single excitation "|j>"
if(i_ok .ne. 1)cycle
double precision :: coef_ijrs_kv,phase_double_triple
call get_excitation(key_tmp,dets_in(1,1,index_singles(i)),exc,degree,phase_double_triple,N_int)
call get_excitation(ref_bitmask,dets_in(1,1,index_singles(j)),exc,degree,phase_other_single_ref,N_int)
call get_excitation(key_tmp,dets_in(1,1,index_singles(j)),exc,degree,phase_other_single_ref,N_int)
coef_ijrs_kv = u_in(index_doubles(i),1) * u_in(index_singles(j),1) * c_ref * c_ref &
* phase_double_triple * phase_other_single_ref
call bitstring_to_list_ab(key_tmp, occ, n_occ, N_int)
do k=1,elec_alpha_num
l = occ(k,1)
density_matrix_alpha(l,l) += coef_ijrs_kv*coef_ijrs_kv
enddo
do k=1,elec_beta_num
l = occ(k,1)
density_matrix_beta(l,l) += coef_ijrs_kv*coef_ijrs_kv
enddo
norm += coef_ijrs_kv* coef_ijrs_kv
if(hole_particles_singles(j,3) == 1)then ! single alpha
density_matrix_alpha(h2,p2) += coef_ijrs_kv * phase_double_triple * u_in(index_doubles(i),1) * c_ref
density_matrix_alpha(p2,h2) += coef_ijrs_kv * phase_double_triple * u_in(index_doubles(i),1) * c_ref
else
density_matrix_beta(h2,p2) += coef_ijrs_kv * phase_double_triple * u_in(index_doubles(i),1) * c_ref
density_matrix_beta(p2,h2) += coef_ijrs_kv * phase_double_triple * u_in(index_doubles(i),1) * c_ref
endif
enddo
enddo
print*,'norm = ',norm
norm = 1.d0/norm
do i = 1, mo_tot_num
do j = 1, mo_tot_num
density_matrix_alpha(i,j) *= norm
density_matrix_beta(i,j) *= norm
enddo
enddo
coef_ijrs = 0.d0
do i = 1, mo_tot_num
coef_ijrs += density_matrix_beta(i,i) + density_matrix_beta(i,i)
enddo
print*,'accu = ',coef_ijrs
end

2
plugins/FOBOCI/fobo_scf.irp.f
View File

@ -1,6 +1,6 @@
program foboscf program foboscf
implicit none implicit none
!call run_prepare call run_prepare
no_oa_or_av_opt = .True. no_oa_or_av_opt = .True.
touch no_oa_or_av_opt touch no_oa_or_av_opt
call routine_fobo_scf call routine_fobo_scf

4
plugins/FOBOCI/foboci_lmct_mlct_threshold_old.irp.f
View File

@ -81,6 +81,8 @@ subroutine FOBOCI_lmct_mlct_old_thr
call set_bitmask_particl_as_input(reunion_of_bitmask) call set_bitmask_particl_as_input(reunion_of_bitmask)
call set_bitmask_hole_as_input(reunion_of_bitmask) call set_bitmask_hole_as_input(reunion_of_bitmask)
call all_single call all_single
call make_s2_eigenfunction
call diagonalize_ci
! if(dressing_2h2p)then ! if(dressing_2h2p)then
! call diag_dressed_2h2p_hamiltonian_and_update_psi_det(i_hole_osoci,lmct) ! call diag_dressed_2h2p_hamiltonian_and_update_psi_det(i_hole_osoci,lmct)
! endif ! endif
@ -193,6 +195,8 @@ subroutine FOBOCI_lmct_mlct_old_thr
! call all_single_split(psi_det_generators,psi_coef_generators,N_det_generators,dressing_matrix) ! call all_single_split(psi_det_generators,psi_coef_generators,N_det_generators,dressing_matrix)
! call diag_dressed_matrix_and_set_to_psi_det(psi_det_generators,N_det_generators,dressing_matrix) ! call diag_dressed_matrix_and_set_to_psi_det(psi_det_generators,N_det_generators,dressing_matrix)
call all_single call all_single
call make_s2_eigenfunction
call diagonalize_ci
! if(dressing_2h2p)then ! if(dressing_2h2p)then
! call diag_dressed_2h2p_hamiltonian_and_update_psi_det(i_particl_osoci,lmct) ! call diag_dressed_2h2p_hamiltonian_and_update_psi_det(i_particl_osoci,lmct)
! endif ! endif

6
plugins/Properties/delta_rho.irp.f
View File

@ -3,9 +3,9 @@
&BEGIN_PROVIDER [double precision, z_max] &BEGIN_PROVIDER [double precision, z_max]
&BEGIN_PROVIDER [double precision, delta_z] &BEGIN_PROVIDER [double precision, delta_z]
implicit none implicit none
z_min = -20.d0 z_min = 0.d0
z_max = 20.d0 z_max = 10.d0
delta_z = 0.1d0 delta_z = 0.005d0
N_z_pts = (z_max - z_min)/delta_z N_z_pts = (z_max - z_min)/delta_z
print*,'N_z_pts = ',N_z_pts print*,'N_z_pts = ',N_z_pts

15
plugins/Properties/mulliken.irp.f
View File

@ -14,13 +14,16 @@ BEGIN_PROVIDER [double precision, spin_population, (ao_num_align,ao_num)]
enddo enddo
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [double precision, spin_population_angular_momentum, (0:ao_l_max)] BEGIN_PROVIDER [double precision, spin_population_angular_momentum, (0:ao_l_max)]
&BEGIN_PROVIDER [double precision, spin_population_angular_momentum_per_atom, (0:ao_l_max,nucl_num)]
implicit none implicit none
integer :: i integer :: i
double precision :: accu double precision :: accu
spin_population_angular_momentum = 0.d0 spin_population_angular_momentum = 0.d0
spin_population_angular_momentum_per_atom = 0.d0
do i = 1, ao_num do i = 1, ao_num
spin_population_angular_momentum(ao_l(i)) += spin_gross_orbital_product(i) spin_population_angular_momentum(ao_l(i)) += spin_gross_orbital_product(i)
spin_population_angular_momentum_per_atom(ao_l(i),ao_nucl(i)) += spin_gross_orbital_product(i)
enddo enddo
END_PROVIDER END_PROVIDER
@ -133,6 +136,16 @@ subroutine print_mulliken_sd
print*,' ',trim(l_to_charater(i)),spin_population_angular_momentum(i) print*,' ',trim(l_to_charater(i)),spin_population_angular_momentum(i)
print*,'sum = ',accu print*,'sum = ',accu
enddo enddo
print*,'Angular momentum analysis per atom'
print*,'Angular momentum analysis'
do j = 1,nucl_num
accu = 0.d0
do i = 0, ao_l_max
accu += spin_population_angular_momentum_per_atom(i,j)
write(*,'(XX,I3,XX,A4,X,A4,X,F10.7)')j,trim(element_name(int(nucl_charge(j)))),trim(l_to_charater(i)),spin_population_angular_momentum_per_atom(i,j)
print*,'sum = ',accu
enddo
enddo
end end

36
plugins/Properties/print_spin_density.irp.f Normal file
View File

@ -0,0 +1,36 @@
program print_sd
implicit none
read_wf = .True.
touch read_wf
call routine
end
subroutine routine
implicit none
integer :: i,j,k
double precision :: z
double precision :: r(3),accu,accu_alpha,accu_beta,tmp
double precision, allocatable :: aos_array(:)
allocate(aos_array(ao_num))
r = 0.d0
r(3) = z_min
do i = 1, N_z_pts
call give_all_aos_at_r(r,aos_array)
accu = 0.d0
accu_alpha = 0.d0
accu_beta = 0.d0
do j = 1, ao_num
do k = 1, ao_num
tmp = aos_array(k) * aos_array(j)
accu += one_body_spin_density_ao(k,j) * tmp
accu_alpha += one_body_dm_ao_alpha(k,j) * tmp
accu_beta += one_body_dm_ao_beta(k,j) * tmp
enddo
enddo
r(3) += delta_z
write(33,'(100(f16.10,X))')r(3),accu,accu_alpha,accu_beta
enddo
end

11
plugins/Properties/provide_deltarho.irp.f Normal file
View File

@ -0,0 +1,11 @@
program pouet
implicit none
read_wf = .True.
touch read_wf
call routine
end
subroutine routine
implicit none
provide integrated_delta_rho_all_points
end

46
plugins/Properties/test_two_body_dm.irp.f
View File

@ -7,28 +7,67 @@ end
subroutine routine subroutine routine
implicit none implicit none
integer :: i,j,k,l integer :: i,j,k,l
integer :: h1,p1,h2,p2,s1,s2
double precision :: accu,get_two_body_dm_ab_map_element,get_mo_bielec_integral_schwartz double precision :: accu,get_two_body_dm_ab_map_element,get_mo_bielec_integral_schwartz
accu = 0.d0 accu = 0.d0
! Diag part of the two body dm ! Diag part of the core two body dm
do i = 1, n_core_orb
h1 = list_core(i)
do j = 1, n_core_orb
h2 = list_core(j)
accu += two_body_dm_ab_diag_core(j,i) * mo_bielec_integral_jj(h1,h2)
enddo
enddo
! Diag part of the active two body dm
do i = 1, n_act_orb do i = 1, n_act_orb
h1 = list_act(i)
do j = 1, n_act_orb do j = 1, n_act_orb
accu += two_body_dm_ab_diag(i,j) * mo_bielec_integral_jj(i,j) h2 = list_act(j)
accu += two_body_dm_ab_diag_act(j,i) * mo_bielec_integral_jj(h1,h2)
enddo
enddo
! Diag part of the active <-> core two body dm
do i = 1, n_act_orb
h1 = list_act(i)
do j = 1, n_core_orb
h2 = list_core(j)
accu += two_body_dm_diag_core_act(j,i) * mo_bielec_integral_jj(h1,h2)
enddo enddo
enddo enddo
print*,'BI ELECTRONIC = ',accu print*,'BI ELECTRONIC = ',accu
double precision :: accu_extra_diag double precision :: accu_extra_diag
accu_extra_diag = 0.d0 accu_extra_diag = 0.d0
! purely active part of the two body dm
do l = 1, n_act_orb ! p2 do l = 1, n_act_orb ! p2
p2 = list_act(l)
do k = 1, n_act_orb ! h2 do k = 1, n_act_orb ! h2
h2 = list_act(k)
do j = 1, n_act_orb ! p1 do j = 1, n_act_orb ! p1
p1 = list_act(j)
do i = 1,n_act_orb ! h1 do i = 1,n_act_orb ! h1
accu_extra_diag += two_body_dm_ab_big_array(i,j,k,l) * get_mo_bielec_integral_schwartz(i,k,j,l,mo_integrals_map) h1 = list_act(i)
accu_extra_diag += two_body_dm_ab_big_array_act(i,j,k,l) * get_mo_bielec_integral_schwartz(h1,h2,p1,p2,mo_integrals_map)
enddo enddo
enddo enddo
enddo enddo
enddo enddo
! core <-> active part of the two body dm
do l = 1, n_act_orb ! p1
p1 = list_act(l)
do k = 1, n_act_orb ! h1
h1 = list_act(k)
do i = 1,n_core_orb ! h2
h2 = list_core(i)
accu_extra_diag += two_body_dm_ab_big_array_core_act(i,k,l) * get_mo_bielec_integral_schwartz(h1,h2,p1,h2,mo_integrals_map)
enddo
enddo
enddo
print*,'extra_diag = ',accu_extra_diag print*,'extra_diag = ',accu_extra_diag
double precision :: average_mono double precision :: average_mono
call get_average(mo_mono_elec_integral,one_body_dm_mo,average_mono) call get_average(mo_mono_elec_integral,one_body_dm_mo,average_mono)
@ -41,7 +80,6 @@ subroutine routine
print*,'<Psi| H |Psi> = ',e_0 + nuclear_repulsion print*,'<Psi| H |Psi> = ',e_0 + nuclear_repulsion
integer :: degree integer :: degree
integer :: exc(0:2,2,2) integer :: exc(0:2,2,2)
integer :: h1,h2,p1,p2,s1,s2
double precision :: phase double precision :: phase
integer :: n_elements integer :: n_elements
n_elements = 0 n_elements = 0

283
plugins/loc_cele/loc_cele.irp.f
View File

@ -92,7 +92,7 @@
nrot(1) = 64 ! number of orbitals to be localized nrot(1) = 6 ! number of orbitals to be localized
integer :: index_rot(1000,1) integer :: index_rot(1000,1)
@ -101,261 +101,72 @@
cmoref = 0.d0 cmoref = 0.d0
irot = 0 irot = 0
! H2 molecule for the mixed localization do i=1,nrot(1)
do i=1,64 irot(i,1) = 19+i
irot(i,1) = i+2
enddo enddo
do i=1,17 ! ESATRIENE with 3 bonding and anti bonding orbitals
cmoref(i+1,i,1)=1.d0 ! First bonding orbital for esa
enddo
cmoref(19,19-1,1)=1.d0
cmoref(20,19-1,1)=-1.d0
cmoref(19,20-1,1)=-1.d0
cmoref(20,20-1,1)=-1.d0
cmoref(21,20-1,1)=2.d0
cmoref(22,21-1,1)=1.d0
cmoref(23,22-1,1)=1.d0
cmoref(24,23-1,1)=1.d0
cmoref(25,24-1,1)=1.d0
cmoref(26,24-1,1)=-1.d0
cmoref(25,25-1,1)=-1.d0
cmoref(26,25-1,1)=-1.d0
cmoref(27,25-1,1)=2.d0
cmoref(28,26-1,1)=1.d0
cmoref(29,27-1,1)=1.d0
cmoref(30,28-1,1)=1.d0
cmoref(31,29-1,1)=1.d0
cmoref(32,29-1,1)=-1.d0
cmoref(31,30-1,1)=-1.d0
cmoref(32,30-1,1)=-1.d0
cmoref(33,30-1,1)=2.d0
cmoref(34,31-1,1)=1.d0
cmoref(35,32-1,1)=1.d0
cmoref(36,33-1,1)=1.d0
do i=33,49
cmoref(i+5,i,1)= 1.d0
enddo
cmoref(55,52-2,1)=1.d0
cmoref(56,52-2,1)=-1.d0
cmoref(55,53-2,1)=-1.d0
cmoref(56,53-2,1)=-1.d0
cmoref(57,53-2,1)=2.d0
cmoref(58,54-2,1)=1.d0
cmoref(59,55-2,1)=1.d0
cmoref(60,56-2,1)=1.d0
cmoref(61,57-2,1)=1.d0
cmoref(62,57-2,1)=-1.d0
cmoref(61,58-2,1)=-1.d0
cmoref(62,58-2,1)=-1.d0
cmoref(63,58-2,1)=2.d0
cmoref(64,59-2,1)=1.d0
cmoref(65,60-2,1)=1.d0
cmoref(66,61-2,1)=1.d0
cmoref(67,62-2,1)=1.d0
cmoref(68,62-2,1)=-1.d0
cmoref(67,63-2,1)=-1.d0
cmoref(68,63-2,1)=-1.d0
cmoref(69,63-2,1)=2.d0
cmoref(70,64-2,1)=1.d0
cmoref(71,65-2,1)=1.d0
cmoref(72,66-2,1)=1.d0
! H2 molecule
! do i=1,66
! irot(i,1) = i
! enddo
!
! do i=1,18
! cmoref(i,i,1)=1.d0
! enddo
! cmoref(19,19,1)=1.d0
! cmoref(20,19,1)=-1.d0
! cmoref(19,20,1)=-1.d0
! cmoref(20,20,1)=-1.d0
! cmoref(21,20,1)=2.d0
! cmoref(22,21,1)=1.d0
! cmoref(23,22,1)=1.d0
! cmoref(24,23,1)=1.d0
!
!
! cmoref(25,24,1)=1.d0
! cmoref(26,24,1)=-1.d0
! cmoref(25,25,1)=-1.d0
! cmoref(26,25,1)=-1.d0
! cmoref(27,25,1)=2.d0
! cmoref(28,26,1)=1.d0
! cmoref(29,27,1)=1.d0
! cmoref(30,28,1)=1.d0
!
! cmoref(31,29,1)=1.d0
! cmoref(32,29,1)=-1.d0
! cmoref(31,30,1)=-1.d0
! cmoref(32,30,1)=-1.d0
! cmoref(33,30,1)=2.d0
! cmoref(34,31,1)=1.d0
! cmoref(35,32,1)=1.d0
! cmoref(36,33,1)=1.d0
!
! do i=34,51
! cmoref(i+3,i,1)= 1.d0
! enddo
!
! cmoref(55,52,1)=1.d0
! cmoref(56,52,1)=-1.d0
! cmoref(55,53,1)=-1.d0
! cmoref(56,53,1)=-1.d0
! cmoref(57,53,1)=2.d0
! cmoref(58,54,1)=1.d0
! cmoref(59,55,1)=1.d0
! cmoref(60,56,1)=1.d0
!
! cmoref(61,57,1)=1.d0
! cmoref(62,57,1)=-1.d0
! cmoref(61,58,1)=-1.d0
! cmoref(62,58,1)=-1.d0
! cmoref(63,58,1)=2.d0
! cmoref(64,59,1)=1.d0
! cmoref(65,60,1)=1.d0
! cmoref(66,61,1)=1.d0
!
! cmoref(67,62,1)=1.d0
! cmoref(68,62,1)=-1.d0
! cmoref(67,63,1)=-1.d0
! cmoref(68,63,1)=-1.d0
! cmoref(69,63,1)=2.d0
! cmoref(70,64,1)=1.d0
! cmoref(71,65,1)=1.d0
! cmoref(72,66,1)=1.d0
! H atom
! do i=1,33
! irot(i,1) = i
! enddo
!
! do i=1,18
! cmoref(i,i,1)=1.d0
! enddo
! cmoref(19,19,1)=1.d0
! cmoref(20,19,1)=-1.d0
! cmoref(19,20,1)=-1.d0
! cmoref(20,20,1)=-1.d0
! cmoref(21,20,1)=2.d0
! cmoref(22,21,1)=1.d0
! cmoref(23,22,1)=1.d0
! cmoref(24,23,1)=1.d0
! cmoref(25,24,1)=1.d0
! cmoref(26,24,1)=-1.d0
! cmoref(25,25,1)=-1.d0
! cmoref(26,25,1)=-1.d0
! cmoref(27,25,1)=2.d0
! cmoref(28,26,1)=1.d0
! cmoref(29,27,1)=1.d0
! cmoref(30,28,1)=1.d0
!
! cmoref(31,29,1)=1.d0
! cmoref(32,29,1)=-1.d0
! cmoref(31,30,1)=-1.d0
! cmoref(32,30,1)=-1.d0
! cmoref(33,30,1)=2.d0
! cmoref(34,31,1)=1.d0
! cmoref(35,32,1)=1.d0
! cmoref(36,33,1)=1.d0
! Definition of the index of the MO to be rotated
! irot(2,1) = 21 ! the first mo to be rotated is the 21 th MO
! irot(3,1) = 22 ! etc....
! irot(4,1) = 23 !
! irot(5,1) = 24 !
! irot(6,1) = 25 !
!N2
! irot(1,1) = 5
! irot(2,1) = 6
! irot(3,1) = 7
! irot(4,1) = 8
! irot(5,1) = 9
! irot(6,1) = 10
!
! cmoref(5,1,1) = 1.d0 !
! cmoref(6,2,1) = 1.d0 !
! cmoref(7,3,1) = 1.d0 !
! cmoref(40,4,1) = 1.d0 !
! cmoref(41,5,1) = 1.d0 !
! cmoref(42,6,1) = 1.d0 !
!END N2
!HEXATRIENE
! irot(1,1) = 20
! irot(2,1) = 21
! irot(3,1) = 22
! irot(4,1) = 23
! irot(5,1) = 24
! irot(6,1) = 25
!
! cmoref(7,1,1) = 1.d0 ! ! cmoref(7,1,1) = 1.d0 !
! cmoref(26,1,1) = 1.d0 ! ! cmoref(26,1,1) = 1.d0 !
! Second bonding orbital for esa
! cmoref(45,2,1) = 1.d0 ! ! cmoref(45,2,1) = 1.d0 !
! cmoref(64,2,1) = 1.d0 ! ! cmoref(64,2,1) = 1.d0 !
! Third bonding orbital for esa
! cmoref(83,3,1) = 1.d0 ! ! cmoref(83,3,1) = 1.d0 !
! cmoref(102,3,1) = 1.d0 ! ! cmoref(102,3,1) = 1.d0 !
! First anti bonding orbital for esa
! cmoref(7,4,1) = 1.d0 ! ! cmoref(7,4,1) = 1.d0 !
! cmoref(26,4,1) = -1.d0 ! ! cmoref(26,4,1) = -1.d0 !
! Second anti bonding orbital for esa
! cmoref(45,5,1) = 1.d0 ! ! cmoref(45,5,1) = 1.d0 !
! cmoref(64,5,1) = -1.d0 ! ! cmoref(64,5,1) = -1.d0 !
! Third anti bonding orbital for esa
! cmoref(83,6,1) = 1.d0 ! ! cmoref(83,6,1) = 1.d0 !
! cmoref(102,6,1) = -1.d0 ! ! cmoref(102,6,1) = -1.d0 !
!END HEXATRIENE
!!!!H2 H2 CAS ! ESATRIENE with 2 bonding and anti bonding orbitals
! irot(1,1) = 1 ! AND 2 radical orbitals
! irot(2,1) = 2 ! First radical orbital
! ! cmoref(7,1,1) = 1.d0 !
! cmoref(1,1,1) = 1.d0 ! First bonding orbital
! cmoref(37,2,1) = 1.d0 ! cmoref(26,2,1) = 1.d0 !
!END H2 ! cmoref(45,2,1) = 1.d0 !
!!!! LOCALIZATION ON THE BASIS FUNCTIONS ! Second bonding orbital
! do i = 1, nrot(1) ! cmoref(64,3,1) = 1.d0 !
! irot(i,1) = i ! cmoref(83,3,1) = 1.d0 !
! cmoref(i,i,1) = 1.d0 ! Second radical orbital for esa
! enddo ! cmoref(102,4,1) = 1.d0 !
! First anti bonding orbital for esa
! cmoref(26,5,1) = 1.d0 !
! cmoref(45,5,1) =-1.d0 !
! Second anti bonding orbital for esa
! cmoref(64,6,1) = 1.d0 !
! cmoref(83,6,1) =-1.d0 !
! ESATRIENE with 1 central bonding and anti bonding orbitals
! AND 4 radical orbitals
! First radical orbital
cmoref(7,1,1) = 1.d0 !
! Second radical orbital
cmoref(26,2,1) = 1.d0 !
! First bonding orbital
cmoref(45,3,1) = 1.d0 !
cmoref(64,3,1) = 1.d0 !
! Third radical orbital for esa
cmoref(83,4,1) = 1.d0 !
! Fourth radical orbital for esa
cmoref(102,5,1) = 1.d0 !
! First anti bonding orbital
cmoref(45,6,1) = 1.d0 !
cmoref(64,6,1) =-1.d0 !
!END BASISLOC
! do i = 1, nrot(1)
! irot(i,1) = 4+i
! enddo
do i = 1, nrot(1) do i = 1, nrot(1)
print*,'irot(i,1) = ',irot(i,1) print*,'irot(i,1) = ',irot(i,1)
enddo enddo
! pause
! you define the guess vectors that you want
! the new MO to be close to
! cmore(i,j,1) = < AO_i | guess_vector_MO(j) >
! i goes from 1 to ao_num
! j goes from 1 to nrot(1)
! Here you must go to the GAMESS output file
! where the AOs are listed and explicited
! From the basis of this knowledge you can build your
! own guess vectors for the MOs
! The new MOs are provided in output
! in the same order than the guess MOs
! do i = 1, nrot(1)
! j = 5+(i-1)*15
! cmoref(j,i,1) = 0.2d0
! cmoref(j+3,i,1) = 0.12d0
! print*,'j = ',j
! enddo
! pause

84
src/Bitmask/bitmasks.irp.f
View File

@ -37,6 +37,30 @@ BEGIN_PROVIDER [ integer(bit_kind), full_ijkl_bitmask_4, (N_int,4) ]
enddo enddo
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ integer(bit_kind), core_inact_act_bitmask_4, (N_int,4) ]
implicit none
integer :: i
do i=1,N_int
core_inact_act_bitmask_4(i,1) = reunion_of_core_inact_act_bitmask(i,1)
core_inact_act_bitmask_4(i,2) = reunion_of_core_inact_act_bitmask(i,1)
core_inact_act_bitmask_4(i,3) = reunion_of_core_inact_act_bitmask(i,1)
core_inact_act_bitmask_4(i,4) = reunion_of_core_inact_act_bitmask(i,1)
enddo
END_PROVIDER
BEGIN_PROVIDER [ integer(bit_kind), virt_bitmask_4, (N_int,4) ]
implicit none
integer :: i
do i=1,N_int
virt_bitmask_4(i,1) = virt_bitmask(i,1)
virt_bitmask_4(i,2) = virt_bitmask(i,1)
virt_bitmask_4(i,3) = virt_bitmask(i,1)
virt_bitmask_4(i,4) = virt_bitmask(i,1)
enddo
END_PROVIDER
BEGIN_PROVIDER [ integer(bit_kind), HF_bitmask, (N_int,2)] BEGIN_PROVIDER [ integer(bit_kind), HF_bitmask, (N_int,2)]
implicit none implicit none
@ -369,11 +393,19 @@ END_PROVIDER
BEGIN_PROVIDER [ integer, list_inact, (n_inact_orb)] BEGIN_PROVIDER [ integer, list_inact, (n_inact_orb)]
&BEGIN_PROVIDER [ integer, list_virt, (n_virt_orb)] &BEGIN_PROVIDER [ integer, list_virt, (n_virt_orb)]
&BEGIN_PROVIDER [ integer, list_inact_reverse, (mo_tot_num)]
&BEGIN_PROVIDER [ integer, list_virt_reverse, (mo_tot_num)]
BEGIN_DOC BEGIN_DOC
! list_inact : List of the inactive orbitals which are supposed to be doubly excited ! list_inact : List of the inactive orbitals which are supposed to be doubly excited
! in post CAS methods ! in post CAS methods
! list_virt : List of vritual orbitals which are supposed to be recieve electrons ! list_virt : List of vritual orbitals which are supposed to be recieve electrons
! in post CAS methods ! in post CAS methods
! list_inact_reverse : reverse list of inactive orbitals
! list_inact_reverse(i) = 0 ::> not an inactive
! list_inact_reverse(i) = k ::> IS the kth inactive
! list_virt_reverse : reverse list of virtual orbitals
! list_virt_reverse(i) = 0 ::> not an virtual
! list_virt_reverse(i) = k ::> IS the kth virtual
END_DOC END_DOC
implicit none implicit none
integer :: occ_inact(N_int*bit_kind_size) integer :: occ_inact(N_int*bit_kind_size)
@ -381,15 +413,20 @@ END_PROVIDER
occ_inact = 0 occ_inact = 0
call bitstring_to_list(inact_bitmask(1,1), occ_inact(1), itest, N_int) call bitstring_to_list(inact_bitmask(1,1), occ_inact(1), itest, N_int)
ASSERT(itest==n_inact_orb) ASSERT(itest==n_inact_orb)
list_inact_reverse = 0
do i = 1, n_inact_orb do i = 1, n_inact_orb
list_inact(i) = occ_inact(i) list_inact(i) = occ_inact(i)
list_inact_reverse(occ_inact(i)) = i
enddo enddo
occ_inact = 0 occ_inact = 0
call bitstring_to_list(virt_bitmask(1,1), occ_inact(1), itest, N_int) call bitstring_to_list(virt_bitmask(1,1), occ_inact(1), itest, N_int)
ASSERT(itest==n_virt_orb) ASSERT(itest==n_virt_orb)
list_virt_reverse = 0
do i = 1, n_virt_orb do i = 1, n_virt_orb
list_virt(i) = occ_inact(i) list_virt(i) = occ_inact(i)
list_virt_reverse(occ_inact(i)) = i
enddo enddo
END_PROVIDER END_PROVIDER
@ -397,7 +434,7 @@ END_PROVIDER
BEGIN_PROVIDER [ integer(bit_kind), reunion_of_core_inact_bitmask, (N_int,2)] BEGIN_PROVIDER [ integer(bit_kind), reunion_of_core_inact_bitmask, (N_int,2)]
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! Reunion of the inactive, active and virtual bitmasks ! Reunion of the core and inactive and virtual bitmasks
END_DOC END_DOC
integer :: i,j integer :: i,j
do i = 1, N_int do i = 1, N_int
@ -407,6 +444,20 @@ END_PROVIDER
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [integer(bit_kind), reunion_of_core_inact_act_bitmask, (N_int,2)]
implicit none
BEGIN_DOC
! Reunion of the core, inactive and active bitmasks
END_DOC
integer :: i,j
do i = 1, N_int
reunion_of_core_inact_act_bitmask(i,1) = ior(reunion_of_core_inact_bitmask(i,1),cas_bitmask(i,1,1))
reunion_of_core_inact_act_bitmask(i,2) = ior(reunion_of_core_inact_bitmask(i,2),cas_bitmask(i,1,1))
enddo
END_PROVIDER
BEGIN_PROVIDER [ integer(bit_kind), reunion_of_bitmask, (N_int,2)] BEGIN_PROVIDER [ integer(bit_kind), reunion_of_bitmask, (N_int,2)]
@ -435,6 +486,7 @@ END_PROVIDER
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [integer, list_core, (n_core_orb)] BEGIN_PROVIDER [integer, list_core, (n_core_orb)]
&BEGIN_PROVIDER [integer, list_core_reverse, (mo_tot_num)]
BEGIN_DOC BEGIN_DOC
! List of the core orbitals that are never excited in post CAS method ! List of the core orbitals that are never excited in post CAS method
END_DOC END_DOC
@ -444,8 +496,10 @@ END_PROVIDER
occ_core = 0 occ_core = 0
call bitstring_to_list(core_bitmask(1,1), occ_core(1), itest, N_int) call bitstring_to_list(core_bitmask(1,1), occ_core(1), itest, N_int)
ASSERT(itest==n_core_orb) ASSERT(itest==n_core_orb)
list_core_reverse = 0
do i = 1, n_core_orb do i = 1, n_core_orb
list_core(i) = occ_core(i) list_core(i) = occ_core(i)
list_core_reverse(occ_core(i)) = i
enddo enddo
END_PROVIDER END_PROVIDER
@ -497,11 +551,17 @@ BEGIN_PROVIDER [ integer, n_act_orb]
do i = 1, N_int do i = 1, N_int
n_act_orb += popcnt(cas_bitmask(i,1,1)) n_act_orb += popcnt(cas_bitmask(i,1,1))
enddo enddo
print*,'n_act_orb = ',n_act_orb
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [integer, list_act, (n_act_orb)] BEGIN_PROVIDER [integer, list_act, (n_act_orb)]
&BEGIN_PROVIDER [integer, list_act_reverse, (mo_tot_num)]
BEGIN_DOC BEGIN_DOC
! list of active orbitals ! list_act(i) = index of the ith active orbital
!
! list_act_reverse : reverse list of active orbitals
! list_act_reverse(i) = 0 ::> not an active
! list_act_reverse(i) = k ::> IS the kth active orbital
END_DOC END_DOC
implicit none implicit none
integer :: occ_act(N_int*bit_kind_size) integer :: occ_act(N_int*bit_kind_size)
@ -509,10 +569,11 @@ BEGIN_PROVIDER [integer, list_act, (n_act_orb)]
occ_act = 0 occ_act = 0
call bitstring_to_list(cas_bitmask(1,1,1), occ_act(1), itest, N_int) call bitstring_to_list(cas_bitmask(1,1,1), occ_act(1), itest, N_int)
ASSERT(itest==n_act_orb) ASSERT(itest==n_act_orb)
list_act_reverse = 0
do i = 1, n_act_orb do i = 1, n_act_orb
list_act(i) = occ_act(i) list_act(i) = occ_act(i)
list_act_reverse(occ_act(i)) = i
enddo enddo
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [integer(bit_kind), closed_shell_ref_bitmask, (N_int,2)] BEGIN_PROVIDER [integer(bit_kind), closed_shell_ref_bitmask, (N_int,2)]
@ -538,3 +599,18 @@ END_PROVIDER
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [integer, n_core_orb_allocate]
implicit none
n_core_orb_allocate = max(n_core_orb,1)
END_PROVIDER
BEGIN_PROVIDER [integer, n_inact_orb_allocate]
implicit none
n_inact_orb_allocate = max(n_inact_orb,1)
END_PROVIDER
BEGIN_PROVIDER [integer, n_virt_orb_allocate]
implicit none
n_virt_orb_allocate = max(n_virt_orb,1)
END_PROVIDER

12
src/Determinants/density_matrix.irp.f
View File

@ -238,17 +238,19 @@ END_PROVIDER
END_DOC END_DOC
implicit none implicit none
integer :: i,j,k,l integer :: i,j,k,l
double precision :: dm_mo double precision :: mo_alpha,mo_beta
one_body_spin_density_ao = 0.d0 one_body_dm_ao_alpha = 0.d0
one_body_dm_ao_beta = 0.d0
do k = 1, ao_num do k = 1, ao_num
do l = 1, ao_num do l = 1, ao_num
do i = 1, mo_tot_num do i = 1, mo_tot_num
do j = 1, mo_tot_num do j = 1, mo_tot_num
dm_mo = one_body_dm_mo_alpha(j,i) mo_alpha = one_body_dm_mo_alpha(j,i)
mo_beta = one_body_dm_mo_beta(j,i)
! if(dabs(dm_mo).le.1.d-10)cycle ! if(dabs(dm_mo).le.1.d-10)cycle
one_body_dm_ao_alpha(l,k) += mo_coef(k,i) * mo_coef(l,j) * dm_mo one_body_dm_ao_alpha(l,k) += mo_coef(k,i) * mo_coef(l,j) * mo_alpha
one_body_dm_ao_beta(l,k) += mo_coef(k,i) * mo_coef(l,j) * dm_mo one_body_dm_ao_beta(l,k) += mo_coef(k,i) * mo_coef(l,j) * mo_beta
enddo enddo
enddo enddo

16
src/Determinants/diagonalize_restart_and_save_all_states.irp.f Normal file
View File

@ -0,0 +1,16 @@
program diag_and_save
implicit none
read_wf = .True.
touch read_wf
call routine
end
subroutine routine
implicit none
call diagonalize_CI
print*,'N_det = ',N_det
call save_wavefunction_general(N_det,N_states_diag,psi_det_sorted,size(psi_coef_sorted,1),psi_coef_sorted)
end

25
src/Determinants/diagonalize_restart_and_save_lowest_state.irp.f Normal file
View File

@ -0,0 +1,25 @@
program diag_and_save
implicit none
read_wf = .True.
touch read_wf
call routine
end
subroutine routine
implicit none
print*,'N_det = ',N_det
call diagonalize_CI
integer :: igood_state
igood_state=1
double precision, allocatable :: psi_coef_tmp(:)
allocate(psi_coef_tmp(n_det))
integer :: i
do i = 1, N_det
psi_coef_tmp(i) = psi_coef(i,igood_state)
enddo
call save_wavefunction_general(N_det,1,psi_det,n_det,psi_coef_tmp)
deallocate(psi_coef_tmp)
end

26
src/Determinants/diagonalize_restart_and_save_one_states.irp.f Normal file
View File

@ -0,0 +1,26 @@
program diag_and_save
implicit none
read_wf = .True.
touch read_wf
call routine
end
subroutine routine
implicit none
print*,'N_det = ',N_det
call diagonalize_CI
write(*,*)'Which state would you like to save ?'
integer :: igood_state
read(5,*)igood_state
double precision, allocatable :: psi_coef_tmp(:)
allocate(psi_coef_tmp(n_det))
integer :: i
do i = 1, N_det
psi_coef_tmp(i) = psi_coef(i,igood_state)
enddo
call save_wavefunction_general(N_det,1,psi_det,n_det,psi_coef_tmp)
deallocate(psi_coef_tmp)
end

27
src/Determinants/diagonalize_restart_and_save_two_states.irp.f Normal file
View File

@ -0,0 +1,27 @@
program diag_and_save
implicit none
read_wf = .True.
touch read_wf
call routine
end
subroutine routine
implicit none
integer :: igood_state_1,igood_state_2
double precision, allocatable :: psi_coef_tmp(:,:)
integer :: i
print*,'N_det = ',N_det
!call diagonalize_CI
write(*,*)'Which couple of states would you like to save ?'
read(5,*)igood_state_1,igood_state_2
allocate(psi_coef_tmp(n_det,2))
do i = 1, N_det
psi_coef_tmp(i,1) = psi_coef(i,igood_state_1)
psi_coef_tmp(i,2) = psi_coef(i,igood_state_2)
enddo
call save_wavefunction_general(N_det,2,psi_det,n_det,psi_coef_tmp)
deallocate(psi_coef_tmp)
end

2
src/Determinants/occ_pattern.irp.f
View File

@ -256,7 +256,7 @@ subroutine make_s2_eigenfunction
integer :: N_det_new integer :: N_det_new
integer, parameter :: bufsze = 1000 integer, parameter :: bufsze = 1000
logical, external :: is_in_wavefunction logical, external :: is_in_wavefunction
return ! return
! !TODO DEBUG ! !TODO DEBUG
! do i=1,N_det ! do i=1,N_det

179
src/Determinants/print_H_matrix_restart.irp.f Normal file
View File

@ -0,0 +1,179 @@
program print_H_matrix_restart
implicit none
read_wf = .True.
touch read_wf
call routine
end
subroutine routine
use bitmasks
implicit none
integer :: i,j
integer, allocatable :: H_matrix_degree(:,:)
double precision, allocatable :: H_matrix_phase(:,:)
integer :: degree
integer(bit_kind), allocatable :: keys_tmp(:,:,:)
allocate(keys_tmp(N_int,2,N_det))
do i = 1, N_det
print*,''
call debug_det(psi_det(1,1,i),N_int)
do j = 1, N_int
keys_tmp(j,1,i) = psi_det(j,1,i)
keys_tmp(j,2,i) = psi_det(j,2,i)
enddo
enddo
if(N_det.ge.10000)then
print*,'Warning !!!'
print*,'Number of determinants is ',N_det
print*,'It means that the H matrix will be enormous !'
print*,'stoppping ..'
stop
endif
print*,''
print*,'Determinants '
do i = 1, N_det
enddo
allocate(H_matrix_degree(N_det,N_det),H_matrix_phase(N_det,N_det))
integer :: exc(0:2,2,2)
double precision :: phase
do i = 1, N_det
do j = i, N_det
call get_excitation_degree(psi_det(1,1,i),psi_det(1,1,j),degree,N_int)
H_matrix_degree(i,j) = degree
H_matrix_degree(j,i) = degree
phase = 0.d0
if(degree==1.or.degree==2)then
call get_excitation(psi_det(1,1,i),psi_det(1,1,j),exc,degree,phase,N_int)
endif
H_matrix_phase(i,j) = phase
H_matrix_phase(j,i) = phase
enddo
enddo
print*,'H matrix '
double precision :: ref_h_matrix,s2
ref_h_matrix = H_matrix_all_dets(1,1)
print*,'HF like determinant energy = ',ref_bitmask_energy+nuclear_repulsion
print*,'Ref element of H_matrix = ',ref_h_matrix+nuclear_repulsion
print*,'Printing the H matrix ...'
print*,''
print*,''
!do i = 1, N_det
! H_matrix_all_dets(i,i) -= ref_h_matrix
!enddo
do i = 1, N_det
H_matrix_all_dets(i,i) += nuclear_repulsion
enddo
!do i = 5,N_det
! H_matrix_all_dets(i,3) = 0.d0
! H_matrix_all_dets(3,i) = 0.d0
! H_matrix_all_dets(i,4) = 0.d0
! H_matrix_all_dets(4,i) = 0.d0
!enddo
do i = 1, N_det
write(*,'(I3,X,A3,1000(F16.7))')i,' | ',H_matrix_all_dets(i,:)
enddo
print*,''
print*,''
print*,''
print*,'Printing the degree of excitations within the H matrix'
print*,''
print*,''
do i = 1, N_det
write(*,'(I3,X,A3,X,1000(I1,X))')i,' | ',H_matrix_degree(i,:)
enddo
print*,''
print*,''
print*,'Printing the phase of the Hamiltonian matrix elements '
print*,''
print*,''
do i = 1, N_det
write(*,'(I3,X,A3,X,1000(F3.0,X))')i,' | ',H_matrix_phase(i,:)
enddo
print*,''
double precision, allocatable :: eigenvectors(:,:), eigenvalues(:)
double precision, allocatable :: s2_eigvalues(:)
allocate (eigenvectors(size(H_matrix_all_dets,1),N_det))
allocate (eigenvalues(N_det))
call lapack_diag(eigenvalues,eigenvectors, &
H_matrix_all_dets,size(H_matrix_all_dets,1),N_det)
print*,'Two first eigenvectors '
do j = 1, n_states
!do j = 1, 1
print*,'State ',j
call get_s2_u0(keys_tmp,eigenvectors(1,j),N_det,size(eigenvectors,1),s2)
print*,'s2 = ',s2
print*,'e = ',eigenvalues(j)
print*,'coefs : '
do i = 1, N_det
print*,'i = ',i,eigenvectors(i,j)
enddo
if(j>1)then
print*,'Delta E(H) = ',eigenvalues(1) - eigenvalues(j)
print*,'Delta E(eV) = ',(eigenvalues(1) - eigenvalues(j))*27.2114d0
endif
enddo
double precision :: get_mo_bielec_integral_schwartz,k_a_iv,k_b_iv
integer :: h1,p1,h2,p2
h1 = 10
p1 = 16
h2 = 14
p2 = 14
!h1 = 1
!p1 = 4
!h2 = 2
!p2 = 2
k_a_iv = get_mo_bielec_integral_schwartz(h1,h2,p2,p1,mo_integrals_map)
h2 = 15
p2 = 15
k_b_iv = get_mo_bielec_integral_schwartz(h1,h2,p2,p1,mo_integrals_map)
print*,'k_a_iv = ',k_a_iv
print*,'k_b_iv = ',k_b_iv
double precision :: k_av,k_bv,k_ai,k_bi
h1 = 16
p1 = 14
h2 = 14
p2 = 16
k_av = get_mo_bielec_integral_schwartz(h1,h2,p1,p2,mo_integrals_map)
h1 = 16
p1 = 15
h2 = 15
p2 = 16
k_bv = get_mo_bielec_integral_schwartz(h1,h2,p1,p2,mo_integrals_map)
h1 = 10
p1 = 14
h2 = 14
p2 = 10
k_ai = get_mo_bielec_integral_schwartz(h1,h2,p1,p2,mo_integrals_map)
h1 = 10
p1 = 15
h2 = 15
p2 = 10
k_bi = get_mo_bielec_integral_schwartz(h1,h2,p1,p2,mo_integrals_map)
print*,'k_av, k_bv = ',k_av,k_bv
print*,'k_ai, k_bi = ',k_ai,k_bi
double precision :: k_iv
h1 = 10
p1 = 16
h2 = 16
p2 = 10
k_iv = get_mo_bielec_integral_schwartz(h1,h2,p1,p2,mo_integrals_map)
print*,'k_iv = ',k_iv
end

11
src/Determinants/print_bitmask.irp.f Normal file
View File

@ -0,0 +1,11 @@
program print_bitmask
implicit none
print*,'core'
call debug_det(core_bitmask,N_int)
print*,'inact'
call debug_det(inact_bitmask,N_int)
print*,'virt'
call debug_det(virt_bitmask,N_int)
end

36
src/Determinants/print_holes_particles.irp.f Normal file
View File

@ -0,0 +1,36 @@
program pouet
implicit none
read_wf = .True.
touch read_wf
call routine
end
subroutine routine
implicit none
integer :: i,j,number_of_holes,number_of_particles
integer :: n_h,n_p
do i = 1, N_det
n_h = number_of_holes(psi_det(1,1,i))
n_p = number_of_particles(psi_det(1,1,i))
if(n_h == 0 .and. n_p == 0)then
print*,'CAS'
else if(n_h == 1 .and. n_p ==0)then
print*,'1h'
else if(n_h == 0 .and. n_p ==1)then
print*,'1p'
else if(n_h == 1 .and. n_p ==1)then
print*,'1h1p'
else if(n_h == 2 .and. n_p ==1)then
print*,'2h1p'
else if(n_h == 1 .and. n_p ==2)then
print*,'1h2p'
else
print*,'PB !! '
call debug_det(psi_det(1,1,i), N_int)
stop
endif
enddo
end

71
src/Determinants/print_wf.irp.f Normal file
View File

@ -0,0 +1,71 @@
program printwf
implicit none
read_wf = .True.
touch read_wf
print*,'ref_bitmask_energy = ',ref_bitmask_energy
call routine
end
subroutine routine
implicit none
integer :: i
integer :: degree
double precision :: hij
integer :: exc(0:2,2,2)
double precision :: phase
integer :: h1,p1,h2,p2,s1,s2
double precision :: get_mo_bielec_integral_schwartz
double precision :: norm_mono_a,norm_mono_b
norm_mono_a = 0.d0
norm_mono_b = 0.d0
do i = 1, min(500,N_det)
print*,''
print*,'i = ',i
call debug_det(psi_det(1,1,i),N_int)
call get_excitation_degree(psi_det(1,1,i),psi_det(1,1,1),degree,N_int)
print*,'degree = ',degree
if(degree == 0)then
print*,'Reference determinant '
else
call i_H_j(psi_det(1,1,i),psi_det(1,1,1),N_int,hij)
call get_excitation(psi_det(1,1,1),psi_det(1,1,i),exc,degree,phase,N_int)
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
print*,'phase = ',phase
if(degree == 1)then
print*,'s1',s1
print*,'h1,p1 = ',h1,p1
if(s1 == 1)then
norm_mono_a += dabs(psi_coef(i,1)/psi_coef(1,1))
else
norm_mono_b += dabs(psi_coef(i,1)/psi_coef(1,1))
endif
print*,'< h | Ka| p > = ',get_mo_bielec_integral_schwartz(h1,list_act(1),list_act(1),p1,mo_integrals_map)
double precision :: hmono,hdouble
call i_H_j_verbose(psi_det(1,1,1),psi_det(1,1,i),N_int,hij,hmono,hdouble)
print*,'hmono = ',hmono
print*,'hdouble = ',hdouble
print*,'hmono+hdouble = ',hmono+hdouble
print*,'hij = ',hij
else
print*,'s1',s1
print*,'h1,p1 = ',h1,p1
print*,'s2',s2
print*,'h2,p2 = ',h2,p2
print*,'< h | Ka| p > = ',get_mo_bielec_integral_schwartz(h1,h2,p1,p2,mo_integrals_map)
endif
print*,'<Ref| H |D_I> = ',hij
endif
print*,'amplitude = ',psi_coef(i,1)/psi_coef(1,1)
enddo
print*,''
print*,''
print*,''
print*,'mono alpha = ',norm_mono_a
print*,'mono beta = ',norm_mono_b
end

50
src/Determinants/save_only_singles.irp.f Normal file
View File

@ -0,0 +1,50 @@
program save_only_singles
implicit none
read_wf = .True.
touch read_wf
call routine
end
subroutine routine
implicit none
integer :: i,j,k,l
use bitmasks
integer :: n_det_restart,degree
integer(bit_kind),allocatable :: psi_det_tmp(:,:,:)
double precision ,allocatable :: psi_coef_tmp(:,:),accu(:)
integer, allocatable :: index_restart(:)
allocate(index_restart(N_det))
N_det_restart = 0
do i = 1, N_det
call get_excitation_degree(psi_det(1,1,1),psi_det(1,1,i),degree,N_int)
if(degree == 0 .or. degree==1)then
N_det_restart +=1
index_restart(N_det_restart) = i
cycle
endif
enddo
allocate (psi_det_tmp(N_int,2,N_det_restart),psi_coef_tmp(N_det_restart,N_states),accu(N_states))
accu = 0.d0
do i = 1, N_det_restart
do j = 1, N_int
psi_det_tmp(j,1,i) = psi_det(j,1,index_restart(i))
psi_det_tmp(j,2,i) = psi_det(j,2,index_restart(i))
enddo
do j = 1,N_states
psi_coef_tmp(i,j) = psi_coef(index_restart(i),j)
accu(j) += psi_coef_tmp(i,j) * psi_coef_tmp(i,j)
enddo
enddo
do j = 1, N_states
accu(j) = 1.d0/dsqrt(accu(j))
enddo
do j = 1,N_states
do i = 1, N_det_restart
psi_coef_tmp(i,j) = psi_coef_tmp(i,j) * accu(j)
enddo
enddo
call save_wavefunction_general(N_det_restart,N_states,psi_det_tmp,N_det_restart,psi_coef_tmp)
deallocate (psi_det_tmp,psi_coef_tmp,accu,index_restart)
end

15
src/Determinants/test_3d.irp.f Normal file
View File

@ -0,0 +1,15 @@
program test_3d
implicit none
integer :: i,npt
double precision :: dx,domain,x_min,x,step_function_becke
domain = 5.d0
npt = 100
dx = domain/dble(npt)
x_min = -0.5d0 * domain
x = x_min
do i = 1, npt
write(33,*)x,step_function_becke(x)
x += dx
enddo
end

18
src/Determinants/test_two_body.irp.f Normal file
View File

@ -0,0 +1,18 @@
program test
implicit none
read_wf = .True.
touch read_wf
call routine
end
subroutine routine
implicit none
integer :: i,j,k,l
do i = 1, n_act_orb
do j = 1, n_act_orb
do k = 1, n_act_orb
enddo
enddo
enddo
end

21
src/Determinants/truncate_wf.irp.f
View File

@ -1,18 +1,11 @@
program cisd program s2_eig_restart
implicit none implicit none
integer :: i,k read_wf = .True.
call routine
end
double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:) subroutine routine
integer :: N_st, degree implicit none
N_det=10000 call make_s2_eigenfunction
do i=1,N_det
do k=1,N_int
psi_det(k,1,i) = psi_det_sorted(k,1,i)
psi_det(k,2,i) = psi_det_sorted(k,2,i)
enddo
psi_coef(i,:) = psi_coef_sorted(i,:)
enddo
TOUCH psi_det psi_coef psi_det_sorted psi_coef_sorted psi_average_norm_contrib_sorted N_det TOUCH psi_det psi_coef psi_det_sorted psi_coef_sorted psi_average_norm_contrib_sorted N_det
call save_wavefunction call save_wavefunction
end end

367
src/Determinants/two_body_dm_map.irp.f
View File

@ -193,33 +193,141 @@ subroutine add_values_to_two_body_dm_map(mask_ijkl)
end end
BEGIN_PROVIDER [double precision, two_body_dm_ab_diag, (mo_tot_num, mo_tot_num)] BEGIN_PROVIDER [double precision, two_body_dm_ab_diag_act, (n_act_orb, n_act_orb)]
&BEGIN_PROVIDER [double precision, two_body_dm_ab_diag_inact, (n_inact_orb_allocate, n_inact_orb_allocate)]
&BEGIN_PROVIDER [double precision, two_body_dm_ab_diag_core, (n_core_orb_allocate, n_core_orb_allocate)]
&BEGIN_PROVIDER [double precision, two_body_dm_ab_diag_all, (mo_tot_num, mo_tot_num)]
&BEGIN_PROVIDER [double precision, two_body_dm_diag_core_a_act_b, (n_core_orb_allocate,n_act_orb)]
&BEGIN_PROVIDER [double precision, two_body_dm_diag_core_b_act_a, (n_core_orb_allocate,n_act_orb)]
&BEGIN_PROVIDER [double precision, two_body_dm_diag_core_act, (n_core_orb_allocate,n_act_orb)]
implicit none implicit none
use bitmasks
integer :: i,j,k,l,m integer :: i,j,k,l,m
integer :: occ(N_int*bit_kind_size,2) integer :: occ(N_int*bit_kind_size,2)
integer :: n_occ_ab(2) integer :: n_occ_ab(2)
integer :: occ_act(N_int*bit_kind_size,2)
integer :: n_occ_ab_act(2)
integer :: occ_core(N_int*bit_kind_size,2)
integer :: n_occ_ab_core(2)
double precision :: contrib double precision :: contrib
BEGIN_DOC BEGIN_DOC
! two_body_dm_ab_diag(k,m) = <\Psi | n_(k\alpha) n_(m\beta) | \Psi> ! two_body_dm_ab_diag_all(k,m) = <\Psi | n_(k\alpha) n_(m\beta) | \Psi>
! two_body_dm_ab_diag_act(k,m) is restricted to the active orbitals :
! orbital k = list_act(k)
! two_body_dm_ab_diag_inact(k,m) is restricted to the inactive orbitals :
! orbital k = list_inact(k)
! two_body_dm_ab_diag_core(k,m) is restricted to the core orbitals :
! orbital k = list_core(k)
! two_body_dm_ab_diag_core_b_act_a(k,m) represents the core beta <-> active alpha part of the two body dm
! orbital k = list_core(k)
! orbital m = list_act(m)
! two_body_dm_ab_diag_core_a_act_b(k,m) represents the core alpha <-> active beta part of the two body dm
! orbital k = list_core(k)
! orbital m = list_act(m)
! two_body_dm_ab_diag_core_act(k,m) represents the core<->active part of the diagonal two body dm
! when we traced on the spin
! orbital k = list_core(k)
! orbital m = list_act(m)
END_DOC END_DOC
integer(bit_kind) :: key_tmp_core(N_int,2)
integer(bit_kind) :: key_tmp_act(N_int,2)
two_body_dm_ab_diag = 0.d0 two_body_dm_ab_diag_all = 0.d0
two_body_dm_ab_diag_act = 0.d0
two_body_dm_ab_diag_core = 0.d0
two_body_dm_ab_diag_inact = 0.d0
two_body_dm_diag_core_a_act_b = 0.d0
two_body_dm_diag_core_b_act_a = 0.d0
two_body_dm_diag_core_act = 0.d0
do i = 1, N_det ! i == |I> do i = 1, N_det ! i == |I>
call bitstring_to_list_ab(psi_det(1,1,i), occ, n_occ_ab, N_int) ! Full diagonal part of the two body dm
contrib = psi_coef(i,1)**2 contrib = psi_coef(i,1)**2
call bitstring_to_list_ab(psi_det(1,1,i), occ, n_occ_ab, N_int)
do j = 1, elec_beta_num do j = 1, elec_beta_num
k = occ(j,2) k = occ(j,2)
do l = 1, elec_alpha_num do l = 1, elec_alpha_num
m = occ(l,1) m = occ(l,1)
two_body_dm_ab_diag(k,m) += 0.5d0 * contrib two_body_dm_ab_diag_all(k,m) += 0.5d0 * contrib
two_body_dm_ab_diag(m,k) += 0.5d0 * contrib two_body_dm_ab_diag_all(m,k) += 0.5d0 * contrib
enddo enddo
enddo enddo
! ACTIVE PART of the diagonal part of the two body dm
do j = 1, N_int
key_tmp_act(j,1) = psi_det(j,1,i)
key_tmp_act(j,2) = psi_det(j,2,i)
enddo
do j = 1, N_int
key_tmp_act(j,1) = iand(key_tmp_act(j,1),cas_bitmask(j,1,1))
key_tmp_act(j,2) = iand(key_tmp_act(j,2),cas_bitmask(j,1,1))
enddo
call bitstring_to_list_ab(key_tmp_act, occ_act, n_occ_ab_act, N_int)
do j = 1,n_occ_ab_act(2)
k = list_act_reverse(occ_act(j,2))
do l = 1, n_occ_ab_act(1)
m = list_act_reverse(occ_act(l,1))
two_body_dm_ab_diag_act(k,m) += 0.5d0 * contrib
two_body_dm_ab_diag_act(m,k) += 0.5d0 * contrib
enddo
enddo
! CORE PART of the diagonal part of the two body dm
do j = 1, N_int
key_tmp_core(j,1) = psi_det(j,1,i)
key_tmp_core(j,2) = psi_det(j,2,i)
enddo
do j = 1, N_int
key_tmp_core(j,1) = iand(key_tmp_core(j,1),core_bitmask(j,1))
key_tmp_core(j,2) = iand(key_tmp_core(j,2),core_bitmask(j,1))
enddo
call bitstring_to_list_ab(key_tmp_core, occ_core, n_occ_ab_core, N_int)
do j = 1,n_occ_ab_core(2)
k = list_core_reverse(occ_core(j,2))
do l = 1, n_occ_ab_core(1)
m = list_core_reverse(occ_core(l,1))
two_body_dm_ab_diag_core(k,m) += 0.5d0 * contrib
two_body_dm_ab_diag_core(m,k) += 0.5d0 * contrib
enddo
enddo
! ACT<->CORE PART
! alpha electron in active space
do j = 1,n_occ_ab_act(1)
k = list_act_reverse(occ_act(j,1))
! beta electron in core space
do l = 1, n_occ_ab_core(2)
m = list_core_reverse(occ_core(l,2))
! The fact that you have 1 * contrib and not 0.5 * contrib
! takes into account the following symmetry :
! 0.5 * <n_k n_m> + 0.5 * <n_m n_k>
two_body_dm_diag_core_b_act_a(m,k) += contrib
enddo
enddo
! beta electron in active space
do j = 1,n_occ_ab_act(2)
k = list_act_reverse(occ_act(j,2))
! alpha electron in core space
do l = 1, n_occ_ab_core(1)
m = list_core_reverse(occ_core(l,1))
! The fact that you have 1 * contrib and not 0.5 * contrib
! takes into account the following symmetry :
! 0.5 * <n_k n_m> + 0.5 * <n_m n_k>
two_body_dm_diag_core_a_act_b(m,k) += contrib
enddo
enddo
enddo
do j = 1, n_core_orb
do l = 1, n_act_orb
two_body_dm_diag_core_act(j,l) = two_body_dm_diag_core_b_act_a(j,l) + two_body_dm_diag_core_a_act_b(j,l)
enddo
enddo enddo
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [double precision, two_body_dm_ab_big_array, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)] BEGIN_PROVIDER [double precision, two_body_dm_ab_big_array_act, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)]
&BEGIN_PROVIDER [double precision, two_body_dm_ab_big_array_core_act, (n_core_orb_allocate,n_act_orb,n_act_orb)]
implicit none implicit none
use bitmasks
integer :: i,j,k,l,m integer :: i,j,k,l,m
integer :: degree integer :: degree
PROVIDE mo_coef psi_coef psi_det PROVIDE mo_coef psi_coef psi_det
@ -229,54 +337,108 @@ BEGIN_PROVIDER [double precision, two_body_dm_ab_big_array, (n_act_orb,n_act_orb
double precision :: contrib double precision :: contrib
integer :: occ(N_int*bit_kind_size,2) integer :: occ(N_int*bit_kind_size,2)
integer :: n_occ_ab(2) integer :: n_occ_ab(2)
two_body_dm_ab_big_array = 0.d0 integer :: occ_core(N_int*bit_kind_size,2)
integer :: n_occ_ab_core(2)
integer(bit_kind) :: key_tmp_i(N_int,2)
integer(bit_kind) :: key_tmp_i_core(N_int,2)
integer(bit_kind) :: key_tmp_j(N_int,2)
two_body_dm_ab_big_array_act = 0.d0
two_body_dm_ab_big_array_core_act = 0.d0
BEGIN_DOC BEGIN_DOC
! The alpha-beta energy can be computed thanks to ! two_body_dm_ab_big_array_act = Purely active part of the two body density matrix
! sum_{h1,p1,h2,p2} two_body_dm_ab_big_array(h1,p1,h2,p2) * (h1p1|h2p2) ! two_body_dm_ab_big_array_act_core takes only into account the single excitation
! within the active space that adds terms in the act <-> core two body dm
! two_body_dm_ab_big_array_act_core(i,j,k) = < a^\dagger_i n_k a_j >
! with i,j in the ACTIVE SPACE
! with k in the CORE SPACE
!
! The alpha-beta extra diagonal energy FOR WF DEFINED AS AN APPROXIMATION OF A CAS can be computed thanks to
! sum_{h1,p1,h2,p2} two_body_dm_ab_big_array_act(h1,p1,h2,p2) * (h1p1|h2p2)
! + sum_{h1,p1,h2,p2} two_body_dm_ab_big_array_core_act(h1,p1,h2,p2) * (h1p1|h2p2)
END_DOC END_DOC
do i = 1, N_det ! i == |I> do i = 1, N_det ! i == |I>
call bitstring_to_list_ab(psi_det(1,1,i), occ, n_occ_ab, N_int) ! active part of psi_det(i)
do j = 1, N_int
key_tmp_i(j,1) = psi_det(j,1,i)
key_tmp_i(j,2) = psi_det(j,2,i)
key_tmp_i_core(j,1) = psi_det(j,1,i)
key_tmp_i_core(j,2) = psi_det(j,2,i)
enddo
do j = 1, N_int
key_tmp_i(j,1) = iand(key_tmp_i(j,1),cas_bitmask(j,1,1))
key_tmp_i(j,2) = iand(key_tmp_i(j,2),cas_bitmask(j,1,1))
enddo
do j = 1, N_int
key_tmp_i_core(j,1) = iand(key_tmp_i_core(j,1),core_bitmask(j,1))
key_tmp_i_core(j,2) = iand(key_tmp_i_core(j,2),core_bitmask(j,1))
enddo
call bitstring_to_list_ab(key_tmp_i_core, occ_core, n_occ_ab_core, N_int)
call bitstring_to_list_ab(key_tmp_i, occ, n_occ_ab, N_int)
do j = i+1, N_det ! j == <J| do j = i+1, N_det ! j == <J|
call get_excitation_degree(psi_det(1,1,i),psi_det(1,1,j),degree,N_int) ! active part of psi_det(j)
do k = 1, N_int
key_tmp_j(k,1) = psi_det(k,1,j)
key_tmp_j(k,2) = psi_det(k,2,j)
enddo
do k = 1, N_int
key_tmp_j(k,1) = iand(key_tmp_j(k,1),cas_bitmask(k,1,1))
key_tmp_j(k,2) = iand(key_tmp_j(k,2),cas_bitmask(k,1,1))
enddo
! control if the two determinants are connected by
! at most a double excitation WITHIN THE ACTIVE SPACE
call get_excitation_degree(key_tmp_i,key_tmp_j,degree,N_int)
if(degree>2)cycle if(degree>2)cycle
! if it is the case, then compute the hamiltonian matrix element with the proper phase
call get_excitation(psi_det(1,1,i),psi_det(1,1,j),exc,degree,phase,N_int) call get_excitation(psi_det(1,1,i),psi_det(1,1,j),exc,degree,phase,N_int)
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
! if(i==3.or.j==3)then
! print*,'i,j = ',i,j
! call debug_det(psi_det(1,1,i),N_int)
! call debug_det(psi_det(1,1,j),N_int)
! print*,degree,s1,s2
! print*,h1,p1,h2,p2
! print*,phase
! pause
! endif
contrib = 0.5d0 * psi_coef(i,1) * psi_coef(j,1) * phase contrib = 0.5d0 * psi_coef(i,1) * psi_coef(j,1) * phase
! print*,'coucou'
! print*,'i,j = ',i,j
! print*,'contrib = ',contrib
! print*,h1,p1,h2,p2
! print*,'s1,s2',s1,s2
! call debug_det(psi_det(1,1,i),N_int)
! call debug_det(psi_det(1,1,j),N_int)
! pause
if(degree==2)then ! case of the DOUBLE EXCITATIONS ************************************ if(degree==2)then ! case of the DOUBLE EXCITATIONS ************************************
if(s1==s2)cycle ! Only the alpha/beta two body density matrix if(s1==s2)cycle ! Only the alpha/beta two body density matrix
! <J| a^{\dagger}_{p1 s1} a^{\dagger}_{p2 s2} a_{h2 s2} a_{h1 s1} |I> * c_I * c_J ! <J| a^{\dagger}_{p1 s1} a^{\dagger}_{p2 s2} a_{h2 s2} a_{h1 s1} |I> * c_I * c_J
call insert_into_two_body_dm_big_array( two_body_dm_ab_big_array,n_act_orb,n_act_orb,n_act_orb,n_act_orb,contrib,h1,p1,h2,p2) h1 = list_act_reverse(h1)
h2 = list_act_reverse(h2)
p1 = list_act_reverse(p1)
p2 = list_act_reverse(p2)
call insert_into_two_body_dm_big_array( two_body_dm_ab_big_array_act,n_act_orb,n_act_orb,n_act_orb,n_act_orb,contrib,h1,p1,h2,p2)
else if(degree==1)then! case of the SINGLE EXCITATIONS *************************************************** else if(degree==1)then! case of the SINGLE EXCITATIONS ***************************************************
print*,'h1 = ',h1
h1 = list_act_reverse(h1)
print*,'h1 = ',h1
print*,'p1 = ',p1
p1 = list_act_reverse(p1)
print*,'p1 = ',p1
if(s1==1)then ! Mono alpha : if(s1==1)then ! Mono alpha :
do k = 1, elec_beta_num ! purely active part of the extra diagonal two body dm
m = occ(k,2) do k = 1, n_occ_ab(2)
m = list_act_reverse(occ(k,2))
! <J| a^{\dagger}_{p1 \alpha} \hat{n}_{m \beta} a_{h1 \alpha} |I> * c_I * c_J ! <J| a^{\dagger}_{p1 \alpha} \hat{n}_{m \beta} a_{h1 \alpha} |I> * c_I * c_J
call insert_into_two_body_dm_big_array( two_body_dm_ab_big_array,n_act_orb,n_act_orb,n_act_orb,n_act_orb,contrib,h1,p1,m,m) call insert_into_two_body_dm_big_array( two_body_dm_ab_big_array_act,n_act_orb,n_act_orb,n_act_orb,n_act_orb,contrib,h1,p1,m,m)
enddo
! core <-> active part of the extra diagonal two body dm
do k = 1, n_occ_ab_core(2)
m = list_core_reverse(occ_core(k,2))
! <J| a^{\dagger}_{p1 \alpha} \hat{n}_{m \beta} a_{h1 \alpha} |I> * c_I * c_J
two_body_dm_ab_big_array_core_act(m,h1,p1) += 2.d0 * contrib
two_body_dm_ab_big_array_core_act(m,p1,h1) += 2.d0 * contrib
enddo enddo
else ! Mono Beta : else ! Mono Beta :
do k = 1, elec_alpha_num ! purely active part of the extra diagonal two body dm
m = occ(k,1) do k = 1, n_occ_ab(1)
m = list_act_reverse(occ(k,1))
! <J| a^{\dagger}_{p1 \beta} \hat{n}_{m \alpha} a_{h1 \beta} |I> * c_I * c_J ! <J| a^{\dagger}_{p1 \beta} \hat{n}_{m \alpha} a_{h1 \beta} |I> * c_I * c_J
call insert_into_two_body_dm_big_array(two_body_dm_ab_big_array,n_act_orb,n_act_orb,n_act_orb,n_act_orb,contrib,h1,p1,m,m) call insert_into_two_body_dm_big_array(two_body_dm_ab_big_array_act,n_act_orb,n_act_orb,n_act_orb,n_act_orb,contrib,h1,p1,m,m)
enddo
! core <-> active part of the extra diagonal two body dm
do k = 1, n_occ_ab_core(1)
m = list_core_reverse(occ_core(k,1))
! <J| a^{\dagger}_{p1 \alpha} \hat{n}_{m \beta} a_{h1 \alpha} |I> * c_I * c_J
two_body_dm_ab_big_array_core_act(m,h1,p1) += 2.d0 * contrib
two_body_dm_ab_big_array_core_act(m,p1,h1) += 2.d0 * contrib
enddo enddo
endif endif
@ -305,28 +467,37 @@ end
double precision function compute_extra_diag_two_body_dm_ab(r1,r2) double precision function compute_extra_diag_two_body_dm_ab(r1,r2)
implicit none implicit none
BEGIN_DOC
! compute the extra diagonal contribution to the alpha/bet two body density at r1, r2
END_DOC
double precision :: r1(3), r2(3)
double precision :: compute_extra_diag_two_body_dm_ab_act,compute_extra_diag_two_body_dm_ab_core_act
compute_extra_diag_two_body_dm_ab = compute_extra_diag_two_body_dm_ab_act(r1,r2)+compute_extra_diag_two_body_dm_ab_core_act(r1,r2)
end
double precision function compute_extra_diag_two_body_dm_ab_act(r1,r2)
implicit none
BEGIN_DOC
! compute the extra diagonal contribution to the two body density at r1, r2
! involving ONLY THE ACTIVE PART, which means that the four index of the excitations
! involved in the two body density matrix are ACTIVE
END_DOC
PROVIDE n_act_orb
double precision, intent(in) :: r1(3),r2(3) double precision, intent(in) :: r1(3),r2(3)
integer :: i,j,k,l integer :: i,j,k,l
double precision :: mos_array_r1(mo_tot_num),mos_array_r2(mo_tot_num) double precision :: mos_array_r1(n_act_orb),mos_array_r2(n_act_orb)
double precision :: contrib double precision :: contrib
compute_extra_diag_two_body_dm_ab = 0.d0 double precision :: contrib_tmp
!call give_all_mos_at_r(r1,mos_array_r1) !print*,'n_act_orb = ',n_act_orb
!call give_all_mos_at_r(r2,mos_array_r2) compute_extra_diag_two_body_dm_ab_act = 0.d0
call give_all_act_mos_at_r(r1,mos_array_r1) call give_all_act_mos_at_r(r1,mos_array_r1)
call give_all_act_mos_at_r(r2,mos_array_r2) call give_all_act_mos_at_r(r2,mos_array_r2)
do l = 1, n_act_orb ! p2 do l = 1, n_act_orb ! p2
do k = 1, n_act_orb ! h2 do k = 1, n_act_orb ! h2
do j = 1, n_act_orb ! p1 do j = 1, n_act_orb ! p1
do i = 1,n_act_orb ! h1 do i = 1,n_act_orb ! h1
double precision :: contrib_tmp
! print*,'i,j',i,j
! print*,mos_array_r1(i) , mos_array_r1(j)
! print*,'k,l',k,l
! print*,mos_array_r2(k) * mos_array_r2(l)
! print*,'gama = ',two_body_dm_ab_big_array(i,j,k,l)
! pause
contrib_tmp = mos_array_r1(i) * mos_array_r1(j) * mos_array_r2(k) * mos_array_r2(l) contrib_tmp = mos_array_r1(i) * mos_array_r1(j) * mos_array_r2(k) * mos_array_r2(l)
compute_extra_diag_two_body_dm_ab += two_body_dm_ab_big_array(i,j,k,l) * contrib_tmp compute_extra_diag_two_body_dm_ab_act += two_body_dm_ab_big_array_act(i,j,k,l) * contrib_tmp
enddo enddo
enddo enddo
enddo enddo
@ -334,13 +505,69 @@ double precision function compute_extra_diag_two_body_dm_ab(r1,r2)
end end
double precision function compute_diag_two_body_dm_ab(r1,r2) double precision function compute_extra_diag_two_body_dm_ab_core_act(r1,r2)
implicit none
BEGIN_DOC
! compute the extra diagonal contribution to the two body density at r1, r2
! involving ONLY THE ACTIVE PART, which means that the four index of the excitations
! involved in the two body density matrix are ACTIVE
END_DOC
double precision, intent(in) :: r1(3),r2(3)
integer :: i,j,k,l
double precision :: mos_array_act_r1(n_act_orb),mos_array_act_r2(n_act_orb)
double precision :: mos_array_core_r1(n_core_orb),mos_array_core_r2(n_core_orb)
double precision :: contrib_core_1,contrib_core_2
double precision :: contrib_act_1,contrib_act_2
double precision :: contrib_tmp
compute_extra_diag_two_body_dm_ab_core_act = 0.d0
call give_all_act_mos_at_r(r1,mos_array_act_r1)
call give_all_act_mos_at_r(r2,mos_array_act_r2)
call give_all_core_mos_at_r(r1,mos_array_core_r1)
call give_all_core_mos_at_r(r2,mos_array_core_r2)
do i = 1, n_act_orb ! h1
do j = 1, n_act_orb ! p1
contrib_act_1 = mos_array_act_r1(i) * mos_array_act_r1(j)
contrib_act_2 = mos_array_act_r2(i) * mos_array_act_r2(j)
do k = 1,n_core_orb ! h2
contrib_core_1 = mos_array_core_r1(k) * mos_array_core_r1(k)
contrib_core_2 = mos_array_core_r2(k) * mos_array_core_r2(k)
contrib_tmp = 0.5d0 * (contrib_act_1 * contrib_core_2 + contrib_act_2 * contrib_core_1)
compute_extra_diag_two_body_dm_ab_core_act += two_body_dm_ab_big_array_core_act(k,i,j) * contrib_tmp
enddo
enddo
enddo
end
double precision function compute_diag_two_body_dm_ab_core(r1,r2)
implicit none implicit none
double precision :: r1(3),r2(3) double precision :: r1(3),r2(3)
integer :: i,j,k,l integer :: i,j,k,l
double precision :: mos_array_r1(mo_tot_num),mos_array_r2(mo_tot_num) double precision :: mos_array_r1(n_core_orb_allocate),mos_array_r2(n_core_orb_allocate)
double precision :: contrib,contrib_tmp double precision :: contrib,contrib_tmp
compute_diag_two_body_dm_ab = 0.d0 compute_diag_two_body_dm_ab_core = 0.d0
call give_all_core_mos_at_r(r1,mos_array_r1)
call give_all_core_mos_at_r(r2,mos_array_r2)
do l = 1, n_core_orb !
contrib = mos_array_r2(l)*mos_array_r2(l)
! if(dabs(contrib).lt.threshld_two_bod_dm)cycle
do k = 1, n_core_orb !
contrib_tmp = contrib * mos_array_r1(k)*mos_array_r1(k)
! if(dabs(contrib).lt.threshld_two_bod_dm)cycle
compute_diag_two_body_dm_ab_core += two_body_dm_ab_diag_core(k,l) * contrib_tmp
enddo
enddo
end
double precision function compute_diag_two_body_dm_ab_act(r1,r2)
implicit none
double precision :: r1(3),r2(3)
integer :: i,j,k,l
double precision :: mos_array_r1(n_act_orb),mos_array_r2(n_act_orb)
double precision :: contrib,contrib_tmp
compute_diag_two_body_dm_ab_act = 0.d0
call give_all_act_mos_at_r(r1,mos_array_r1) call give_all_act_mos_at_r(r1,mos_array_r1)
call give_all_act_mos_at_r(r2,mos_array_r2) call give_all_act_mos_at_r(r2,mos_array_r2)
do l = 1, n_act_orb ! do l = 1, n_act_orb !
@ -349,10 +576,44 @@ double precision function compute_diag_two_body_dm_ab(r1,r2)
do k = 1, n_act_orb ! do k = 1, n_act_orb !
contrib_tmp = contrib * mos_array_r1(k)*mos_array_r1(k) contrib_tmp = contrib * mos_array_r1(k)*mos_array_r1(k)
! if(dabs(contrib).lt.threshld_two_bod_dm)cycle ! if(dabs(contrib).lt.threshld_two_bod_dm)cycle
compute_diag_two_body_dm_ab += two_body_dm_ab_diag(k,l) * contrib_tmp compute_diag_two_body_dm_ab_act += two_body_dm_ab_diag_act(k,l) * contrib_tmp
enddo enddo
enddo enddo
end end
double precision function compute_diag_two_body_dm_ab_core_act(r1,r2)
implicit none
double precision :: r1(3),r2(3)
integer :: i,j,k,l
double precision :: mos_array_core_r1(n_core_orb_allocate),mos_array_core_r2(n_core_orb_allocate)
double precision :: mos_array_act_r1(n_act_orb),mos_array_act_r2(n_act_orb)
double precision :: contrib_core_1,contrib_core_2
double precision :: contrib_act_1,contrib_act_2
double precision :: contrib_tmp
compute_diag_two_body_dm_ab_core_act = 0.d0
call give_all_act_mos_at_r(r1,mos_array_act_r1)
call give_all_act_mos_at_r(r2,mos_array_act_r2)
call give_all_core_mos_at_r(r1,mos_array_core_r1)
call give_all_core_mos_at_r(r2,mos_array_core_r2)
! if(dabs(contrib).lt.threshld_two_bod_dm)cycle
do k = 1, n_act_orb !
contrib_act_1 = mos_array_act_r1(k) * mos_array_act_r1(k)
contrib_act_2 = mos_array_act_r2(k) * mos_array_act_r2(k)
contrib_tmp = 0.5d0 * (contrib_act_1 * contrib_act_2 + contrib_act_2 * contrib_act_1)
! if(dabs(contrib).lt.threshld_two_bod_dm)cycle
do l = 1, n_core_orb !
contrib_core_1 = mos_array_core_r1(l) * mos_array_core_r1(l)
contrib_core_2 = mos_array_core_r2(l) * mos_array_core_r2(l)
compute_diag_two_body_dm_ab_core_act += two_body_dm_diag_core_act(l,k) * contrib_tmp
enddo
enddo
end
double precision function compute_diag_two_body_dm_ab(r1,r2)
implicit none
double precision,intent(in) :: r1(3),r2(3)
double precision :: compute_diag_two_body_dm_ab_act,compute_diag_two_body_dm_ab_core
double precision :: compute_diag_two_body_dm_ab_core_act
compute_diag_two_body_dm_ab = compute_diag_two_body_dm_ab_act(r1,r2)+compute_diag_two_body_dm_ab_core(r1,r2) &
+ compute_diag_two_body_dm_ab_core_act(r1,r2)
end

7
src/Integrals_Bielec/EZFIO.cfg
View File

@ -5,6 +5,13 @@ interface: ezfio,provider,ocaml
default: False default: False
ezfio_name: direct ezfio_name: direct
[no_vvvv_integrals]
type: logical
doc: If True, do not compute the bielectronic integrals when 4 indices are virtual
interface: ezfio,provider,ocaml
default: False
ezfio_name: None
[disk_access_mo_integrals] [disk_access_mo_integrals]
type: Disk_access type: Disk_access
doc: Read/Write MO integrals from/to disk [ Write | Read | None ] doc: Read/Write MO integrals from/to disk [ Write | Read | None ]

38
src/Integrals_Bielec/mo_bi_integrals.irp.f
View File

@ -21,6 +21,7 @@ end
BEGIN_PROVIDER [ logical, mo_bielec_integrals_in_map ] BEGIN_PROVIDER [ logical, mo_bielec_integrals_in_map ]
implicit none implicit none
integer(bit_kind) :: mask_ijkl(N_int,4)
BEGIN_DOC BEGIN_DOC
! If True, the map of MO bielectronic integrals is provided ! If True, the map of MO bielectronic integrals is provided
@ -36,7 +37,44 @@ BEGIN_PROVIDER [ logical, mo_bielec_integrals_in_map ]
endif endif
endif endif
if(no_vvvv_integrals)then
integer :: i,j,k,l
! (core+inact+act) ^ 4
do i = 1,N_int
mask_ijkl(i,1) = core_inact_act_bitmask_4(i,1)
mask_ijkl(i,2) = core_inact_act_bitmask_4(i,2)
mask_ijkl(i,3) = core_inact_act_bitmask_4(i,3)
mask_ijkl(i,4) = core_inact_act_bitmask_4(i,4)
enddo
call add_integrals_to_map(mask_ijkl)
! (core+inact+act) ^ 3 (virt) ^1
do i = 1,N_int
mask_ijkl(i,1) = core_inact_act_bitmask_4(i,1)
mask_ijkl(i,2) = core_inact_act_bitmask_4(i,2)
mask_ijkl(i,3) = core_inact_act_bitmask_4(i,3)
mask_ijkl(i,4) = virt_bitmask(i,1)
enddo
call add_integrals_to_map(mask_ijkl)
! (core+inact+act) ^ 2 (virt) ^2
do i = 1,N_int
mask_ijkl(i,1) = core_inact_act_bitmask_4(i,1)
mask_ijkl(i,2) = core_inact_act_bitmask_4(i,2)
mask_ijkl(i,3) = virt_bitmask(i,1)
mask_ijkl(i,4) = virt_bitmask(i,1)
enddo
call add_integrals_to_map(mask_ijkl)
! (core+inact+act) ^ 1 (virt) ^3
do i = 1,N_int
mask_ijkl(i,1) = core_inact_act_bitmask_4(i,1)
mask_ijkl(i,2) = virt_bitmask(i,1)
mask_ijkl(i,3) = virt_bitmask(i,1)
mask_ijkl(i,4) = virt_bitmask(i,1)
enddo
call add_integrals_to_map(mask_ijkl)
else
call add_integrals_to_map(full_ijkl_bitmask_4) call add_integrals_to_map(full_ijkl_bitmask_4)
endif
END_PROVIDER END_PROVIDER
subroutine add_integrals_to_map(mask_ijkl) subroutine add_integrals_to_map(mask_ijkl)

20
src/MO_Basis/mo_permutation.irp.f Normal file
View File

@ -0,0 +1,20 @@
program permut_mos
implicit none
integer :: mo1,mo2
integer :: i,j,k,l
double precision :: mo_coef_tmp(ao_num_align,2)
print*,'Which MOs would you like to change ?'
read(5,*)mo1,mo2
print*,''
do i= 1,ao_num
mo_coef_tmp(i,1) = mo_coef(i,mo1)
mo_coef_tmp(i,2) = mo_coef(i,mo2)
enddo
do i = 1,ao_num
mo_coef(i,mo1) = mo_coef_tmp(i,2)
mo_coef(i,mo2) = mo_coef_tmp(i,1)
enddo
touch mo_coef
call save_mos
end

53
src/MO_Basis/print_aos.irp.f Normal file
View File

@ -0,0 +1,53 @@
program pouet
implicit none
integer :: i,j,k
double precision :: r(3)
double precision, allocatable :: aos_array(:),mos_array(:),ao_ortho_array(:)
allocate(aos_array(ao_num),mos_array(mo_tot_num), ao_ortho_array(ao_num))
integer :: nx,ny
double precision :: interval_x
double precision :: xmin,xmax
double precision :: dx
double precision :: interval_y
double precision :: ymin,ymax
double precision :: dy
double precision :: val_max
!do i = 1, ao_num
! write(41,'(100(F16.10,X))'),ao_ortho_canonical_overlap(i,:)
!enddo
!stop
xmin = nucl_coord(1,1)-6.d0
xmax = nucl_coord(2,1)+6.d0
interval_x = xmax - xmin
!interval_x = nucl_dist(1,3)
nx = 500
dx = interval_x/dble(nx)
!dx = dabs(interval_x)/dble(nx) * 1.d0/sqrt(2.d0)
r = 0.d0
r(3) = xmin
!r(2) = nucl_coord(1,2)
!r(3) = nucl_coord(1,3)
!r(1) = nucl_coord(2,1)
!r(2) = 1.D0
!r(3) = nucl_coord(2,3)
double precision :: dr(3)
!dr = 0.d0
!dr(1) = -dx
!dr(3) = dx
do j = 1, nx+1
call give_all_mos_at_r(r,mos_array)
write(37,'(100(F16.10,X))') r(3),mos_array(1)*mos_array(1) , mos_array(2)*mos_array(2), mos_array(1)*mos_array(2)
write(38,'(100(F16.10,X))') r(3),mos_array(1), mos_array(2), mos_array(1)*mos_array(2)
! write(38,'(100(F16.10,X))') r(3),mos_array(10), mos_array(2) - 0.029916d0 * mos_array(10),mos_array(2) + 0.029916d0 * mos_array(10)
r(3) += dx
! r += dr
enddo
deallocate(aos_array,mos_array, ao_ortho_array)
end

50
src/MO_Basis/print_mo_in_space.irp.f Normal file
View File

@ -0,0 +1,50 @@
program pouet
implicit none
integer :: i,j,k
double precision :: r(3)
double precision, allocatable :: aos_array(:),mos_array(:),ao_ortho_array(:)
allocate(aos_array(ao_num),mos_array(mo_tot_num), ao_ortho_array(ao_num))
integer :: nx,ny
double precision :: interval_x
double precision :: xmin,xmax
double precision :: dx
double precision :: interval_y
double precision :: ymin,ymax
double precision :: dy
double precision :: val_max
!do i = 1, ao_num
! write(41,'(100(F16.10,X))'),ao_ortho_canonical_overlap(i,:)
!enddo
!stop
xmin = -4.d0
xmax = 4.d0
interval_x = xmax - xmin
nx = 100
dx = dabs(interval_x)/dble(nx)
r = 0.d0
!r(3) = xmin
r(1) = xmin
val_max = 0.d0
do j = 1, nx
! call give_all_aos_at_r(r,aos_array)
call give_all_mos_at_r(r,mos_array)
write(36,'(100(F16.10,X))') r(1), mos_array(1), mos_array(2), mos_array(1)* mos_array(2)
!write(36,'(100(F16.10,X))') r(1), mos_array(1), mos_array(2), mos_array(4)
!write(37,'(100(F16.10,X))') r(1),mos_array(1) * mos_array(2), mos_array(4)*mos_array(2)
! if(val_max.le.aos_array(1) * aos_array(2) )then
! val_max = aos_array(1) * aos_array(2)
! endif
r(1) += dx
! r(3) += dx
enddo
!write(40,'(100(F16.10,X))')nucl_coord(1,2),nucl_coord(1,3),val_max * 1.5d0
!write(41,'(100(F16.10,X))')nucl_coord(2,2),nucl_coord(2,3),val_max * 1.5d0
deallocate(aos_array,mos_array, ao_ortho_array)
end

112
src/Nuclei/atomic_radii.irp.f Normal file
View File

@ -0,0 +1,112 @@
BEGIN_PROVIDER [ double precision, slater_bragg_radii, (100)]
implicit none
BEGIN_DOC
! atomic radii in Angstrom defined in table I of JCP 41, 3199 (1964) Slater
! execpt for the Hydrogen atom where we took the value of Becke (1988, JCP)
END_DOC
slater_bragg_radii = 0.d0
slater_bragg_radii(1) = 0.35d0
slater_bragg_radii(2) = 0.35d0
slater_bragg_radii(3) = 1.45d0
slater_bragg_radii(4) = 1.05d0
slater_bragg_radii(5) = 0.85d0
slater_bragg_radii(6) = 0.70d0
slater_bragg_radii(7) = 0.65d0
slater_bragg_radii(8) = 0.60d0
slater_bragg_radii(9) = 0.50d0
slater_bragg_radii(10) = 0.45d0
slater_bragg_radii(11) = 1.80d0
slater_bragg_radii(12) = 1.70d0
slater_bragg_radii(13) = 1.50d0
slater_bragg_radii(14) = 1.25d0
slater_bragg_radii(15) = 1.10d0
slater_bragg_radii(16) = 1.00d0
slater_bragg_radii(17) = 1.00d0
slater_bragg_radii(18) = 1.00d0
slater_bragg_radii(19) = 2.20d0
slater_bragg_radii(20) = 1.80d0
slater_bragg_radii(21) = 1.60d0
slater_bragg_radii(22) = 1.40d0
slater_bragg_radii(23) = 1.34d0
slater_bragg_radii(24) = 1.40d0
slater_bragg_radii(25) = 1.40d0
slater_bragg_radii(26) = 1.40d0
slater_bragg_radii(27) = 1.35d0
slater_bragg_radii(28) = 1.35d0
slater_bragg_radii(29) = 1.35d0
slater_bragg_radii(30) = 1.35d0
slater_bragg_radii(31) = 1.30d0
slater_bragg_radii(32) = 1.25d0
slater_bragg_radii(33) = 1.15d0
slater_bragg_radii(34) = 1.15d0
slater_bragg_radii(35) = 1.15d0
slater_bragg_radii(36) = 1.15d0
END_PROVIDER
BEGIN_PROVIDER [double precision, slater_bragg_radii_ua, (100)]
implicit none
integer :: i
do i = 1, 100
slater_bragg_radii_ua(i) = slater_bragg_radii(i) * 1.889725989d0
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, slater_bragg_radii_per_atom, (nucl_num)]
implicit none
integer :: i
do i = 1, nucl_num
slater_bragg_radii_per_atom(i) = slater_bragg_radii(int(nucl_charge(i)))
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, slater_bragg_radii_per_atom_ua, (nucl_num)]
implicit none
integer :: i
do i = 1, nucl_num
slater_bragg_radii_per_atom_ua(i) = slater_bragg_radii_ua(int(nucl_charge(i)))
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, slater_bragg_type_inter_distance, (nucl_num, nucl_num)]
implicit none
integer :: i,j
double precision :: xhi_tmp,u_ij
slater_bragg_type_inter_distance = 0.d0
do i = 1, nucl_num
do j = i+1, nucl_num
xhi_tmp = slater_bragg_radii_per_atom(i) / slater_bragg_radii_per_atom(j)
u_ij = (xhi_tmp - 1.d0 ) / (xhi_tmp +1.d0)
slater_bragg_type_inter_distance(i,j) = u_ij / (u_ij * u_ij - 1.d0)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, slater_bragg_type_inter_distance_ua, (nucl_num, nucl_num)]
implicit none
integer :: i,j
double precision :: xhi_tmp,u_ij
slater_bragg_type_inter_distance_ua = 0.d0
do i = 1, nucl_num
do j = i+1, nucl_num
xhi_tmp = slater_bragg_radii_per_atom_ua(i) / slater_bragg_radii_per_atom_ua(j)
u_ij = (xhi_tmp - 1.d0 ) / (xhi_tmp +1.d0)
slater_bragg_type_inter_distance_ua(i,j) = u_ij / (u_ij * u_ij - 1.d0)
if(slater_bragg_type_inter_distance_ua(i,j).gt.0.5d0)then
slater_bragg_type_inter_distance_ua(i,j) = 0.5d0
else if( slater_bragg_type_inter_distance_ua(i,j) .le.-0.5d0)then
slater_bragg_type_inter_distance_ua(i,j) = -0.5d0
endif
enddo
enddo
END_PROVIDER

2264
src/Utils/angular_integration.irp.f Normal file
View File

File diff suppressed because it is too large Load Diff