mirror of
https://github.com/LCPQ/quantum_package
synced 2024-08-25 05:51:46 +02:00
Merge branch 'master' of github.com:LCPQ/quantum_package
This commit is contained in:
commit
d09b56c60e
@ -23,12 +23,15 @@ Documentation
|
||||
.. Do not edit this section. It was auto-generated from the
|
||||
.. NEEDED_MODULES file.
|
||||
|
||||
`fill_h_apply_buffer_cisd <http://github.com/LCPQ/quantum_package/tree/master/src/CISD/H_apply.irp.f#L6>`_
|
||||
`fill_h_apply_buffer_cisd <http://github.com/LCPQ/quantum_package/tree/master/src/CISD/H_apply.irp.f#L/subroutine fill_H_apply_buffer_cisd(n_selected,det_buffer,Nint)/;">`_
|
||||
Fill the H_apply buffer with determiants for CISD
|
||||
|
||||
`h_apply_cisd <http://github.com/LCPQ/quantum_package/tree/master/src/CISD/H_apply.irp.f#L43>`_
|
||||
`h_apply_cisd <http://github.com/LCPQ/quantum_package/tree/master/src/CISD/H_apply.irp.f#L/subroutine H_apply_cisd/;">`_
|
||||
Calls H_apply on the HF determinant and selects all connected single and double
|
||||
excitations (of the same symmetry).
|
||||
|
||||
`cisd <http://github.com/LCPQ/quantum_package/tree/master/src/CISD/cisd.irp.f#L/subroutine cisd/;">`_
|
||||
Undocumented
|
||||
|
||||
|
||||
|
||||
|
421
src/CISD/SC2.irp.f
Normal file
421
src/CISD/SC2.irp.f
Normal file
@ -0,0 +1,421 @@
|
||||
subroutine CISD_SC2(dets_in,u_in,energies,dim_in,sze,N_st,Nint)
|
||||
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)
|
||||
PROVIDE ref_bitmask_energy
|
||||
ASSERT (N_st > 0)
|
||||
ASSERT (sze > 0)
|
||||
ASSERT (Nint > 0)
|
||||
ASSERT (Nint == N_int)
|
||||
integer :: iter
|
||||
integer :: i,j,k,l,m
|
||||
logical :: converged
|
||||
double precision :: overlap(N_st,N_st)
|
||||
double precision :: u_dot_v, u_dot_u
|
||||
|
||||
integer :: degree,N_double,index_hf,index_double(sze)
|
||||
double precision :: hij_elec, e_corr_double,e_corr,diag_h_mat_elem,inv_c0
|
||||
double precision :: e_corr_array(sze),H_jj_ref(sze),H_jj_dressed(sze),hij_double(sze)
|
||||
double precision :: e_corr_double_before,accu,cpu_2,cpu_1
|
||||
integer :: i_ok
|
||||
|
||||
!$OMP PARALLEL DEFAULT(NONE) &
|
||||
!$OMP SHARED(sze,N_st, &
|
||||
!$OMP H_jj_ref,Nint,dets_in,u_in) &
|
||||
!$OMP PRIVATE(i)
|
||||
|
||||
!$OMP DO
|
||||
do i=1,sze
|
||||
H_jj_ref(i) = diag_h_mat_elem(dets_in(1,1,i),Nint)
|
||||
enddo
|
||||
!$OMP END DO NOWAIT
|
||||
!$OMP END PARALLEL
|
||||
|
||||
N_double = 0
|
||||
e_corr = 0.d0
|
||||
e_corr_double = 0.d0
|
||||
do i = 1, sze
|
||||
call get_excitation_degree(ref_bitmask,dets_in(1,1,i),degree,Nint)
|
||||
if(degree==0)then
|
||||
index_hf=i
|
||||
else if (degree == 2)then
|
||||
N_double += 1
|
||||
index_double(N_double) = i
|
||||
call i_H_j(ref_bitmask,dets_in(1,1,i),Nint,hij_elec)
|
||||
hij_double(N_double) = hij_elec
|
||||
e_corr_array(N_double) = u_in(i,1)* hij_elec
|
||||
e_corr_double += e_corr_array(N_double)
|
||||
e_corr += e_corr_array(N_double)
|
||||
index_double(N_double) = i
|
||||
else if (degree == 1)then
|
||||
call i_H_j(ref_bitmask,dets_in(1,1,i),Nint,hij_elec)
|
||||
print*,hij_elec
|
||||
e_corr += u_in(i,1)* hij_elec
|
||||
endif
|
||||
enddo
|
||||
inv_c0 = 1.d0/u_in(index_hf,1)
|
||||
do i = 1, N_double
|
||||
e_corr_array(i) = e_corr_array(i) * inv_c0
|
||||
enddo
|
||||
e_corr = e_corr * inv_c0
|
||||
e_corr_double = e_corr_double * inv_c0
|
||||
print*, 'E_corr = ',e_corr
|
||||
print*, 'E_corr_double = ', e_corr_double
|
||||
|
||||
converged = .False.
|
||||
e_corr_double_before = e_corr_double
|
||||
iter = 0
|
||||
do while (.not.converged)
|
||||
|
||||
iter +=1
|
||||
print*,'SC2 iteration : ',iter
|
||||
call cpu_time(cpu_1)
|
||||
do i=1,sze
|
||||
H_jj_dressed(i) = H_jj_ref(i)
|
||||
if (i==index_hf)cycle
|
||||
accu = 0.d0
|
||||
do j=1,N_double
|
||||
call repeat_excitation(dets_in(1,1,i),ref_bitmask,dets_in(1,1,index_double(j)),i_ok,Nint)
|
||||
if (i_ok==1)cycle! you check if the excitation is possible
|
||||
accu += e_corr_array(j)
|
||||
enddo
|
||||
H_jj_dressed(i) += accu
|
||||
enddo
|
||||
|
||||
call cpu_time(cpu_2)
|
||||
print*,'time for the excitations = ',cpu_2 - cpu_1
|
||||
print*,H_jj_ref(1),H_jj_ref(2)
|
||||
print*,H_jj_dressed(1),H_jj_dressed(2)
|
||||
print*,u_in(index_hf,1),u_in(index_double(1),1)
|
||||
call davidson_diag_hjj(dets_in,u_in,H_jj_dressed,energies,dim_in,sze,N_st,Nint)
|
||||
print*,u_in(index_hf,1),u_in(index_double(1),1)
|
||||
e_corr_double = 0.d0
|
||||
inv_c0 = 1.d0/u_in(index_hf,1)
|
||||
do i = 1, N_double
|
||||
e_corr_array(i) = u_in(index_double(i),1)*inv_c0 * hij_double(i)
|
||||
e_corr_double += e_corr_array(i)
|
||||
enddo
|
||||
print*,'E_corr = ',e_corr_double
|
||||
print*,'delta E_corr =',e_corr_double - e_corr_double_before
|
||||
converged = dabs(e_corr_double - e_corr_double_before) < 1.d-10
|
||||
if (converged) then
|
||||
exit
|
||||
endif
|
||||
e_corr_double_before = e_corr_double
|
||||
|
||||
enddo
|
||||
|
||||
|
||||
|
||||
end
|
||||
|
||||
subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint)
|
||||
use bitmasks
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Davidson diagonalization with specific diagonal elements of the H matrix
|
||||
!
|
||||
! H_jj : specific diagonal H matrix elements to diagonalize de Davidson
|
||||
!
|
||||
! 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) :: H_jj(dim_in)
|
||||
double precision, intent(inout) :: u_in(dim_in,N_st)
|
||||
double precision, intent(out) :: energies(N_st)
|
||||
|
||||
integer :: iter
|
||||
integer :: i,j,k,l,m
|
||||
logical :: converged
|
||||
|
||||
double precision :: overlap(N_st,N_st)
|
||||
double precision :: u_dot_v, u_dot_u
|
||||
|
||||
integer, allocatable :: kl_pairs(:,:)
|
||||
integer :: k_pairs, kl
|
||||
|
||||
integer :: iter2
|
||||
double precision, allocatable :: W(:,:,:), U(:,:,:), R(:,:)
|
||||
double precision, allocatable :: y(:,:,:,:), h(:,:,:,:), lambda(:)
|
||||
double precision :: diag_h_mat_elem
|
||||
double precision :: residual_norm(N_st)
|
||||
|
||||
PROVIDE ref_bitmask_energy
|
||||
|
||||
allocate( &
|
||||
kl_pairs(2,N_st*(N_st+1)/2), &
|
||||
W(sze,N_st,davidson_sze_max), &
|
||||
U(sze,N_st,davidson_sze_max), &
|
||||
R(sze,N_st), &
|
||||
h(N_st,davidson_sze_max,N_st,davidson_sze_max), &
|
||||
y(N_st,davidson_sze_max,N_st,davidson_sze_max), &
|
||||
lambda(N_st*davidson_sze_max))
|
||||
|
||||
ASSERT (N_st > 0)
|
||||
ASSERT (sze > 0)
|
||||
ASSERT (Nint > 0)
|
||||
ASSERT (Nint == N_int)
|
||||
|
||||
! Initialization
|
||||
! ==============
|
||||
|
||||
k_pairs=0
|
||||
do l=1,N_st
|
||||
do k=1,l
|
||||
k_pairs+=1
|
||||
kl_pairs(1,k_pairs) = k
|
||||
kl_pairs(2,k_pairs) = l
|
||||
enddo
|
||||
enddo
|
||||
|
||||
!$OMP PARALLEL DEFAULT(NONE) &
|
||||
!$OMP SHARED(U,sze,N_st,overlap,kl_pairs,k_pairs, &
|
||||
!$OMP Nint,dets_in,u_in) &
|
||||
!$OMP PRIVATE(k,l,kl,i)
|
||||
|
||||
|
||||
! Orthonormalize initial guess
|
||||
! ============================
|
||||
|
||||
!$OMP DO
|
||||
do kl=1,k_pairs
|
||||
k = kl_pairs(1,kl)
|
||||
l = kl_pairs(2,kl)
|
||||
if (k/=l) then
|
||||
overlap(k,l) = u_dot_v(U_in(1,k),U_in(1,l),sze)
|
||||
overlap(l,k) = overlap(k,l)
|
||||
else
|
||||
overlap(k,k) = u_dot_u(U_in(1,k),sze)
|
||||
endif
|
||||
enddo
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
call ortho_lowdin(overlap,size(overlap,1),N_st,U_in,size(U_in,1),sze)
|
||||
|
||||
! Davidson iterations
|
||||
! ===================
|
||||
|
||||
converged = .False.
|
||||
|
||||
do while (.not.converged)
|
||||
|
||||
!$OMP PARALLEL DEFAULT(NONE) &
|
||||
!$OMP PRIVATE(k,i) SHARED(U,u_in,sze,N_st)
|
||||
do k=1,N_st
|
||||
!$OMP DO
|
||||
do i=1,sze
|
||||
U(i,k,1) = u_in(i,k)
|
||||
enddo
|
||||
!$OMP END DO
|
||||
enddo
|
||||
!$OMP END PARALLEL
|
||||
|
||||
do iter=1,davidson_sze_max-1
|
||||
|
||||
! Compute W_k = H |u_k>
|
||||
! ----------------------
|
||||
|
||||
do k=1,N_st
|
||||
call H_u_0(W(1,k,iter),U(1,k,iter),H_jj,sze,dets_in,Nint)
|
||||
enddo
|
||||
|
||||
! Compute h_kl = <u_k | W_l> = <u_k| H |u_l>
|
||||
! -------------------------------------------
|
||||
|
||||
do l=1,N_st
|
||||
do k=1,N_st
|
||||
do iter2=1,iter-1
|
||||
h(k,iter2,l,iter) = u_dot_v(U(1,k,iter2),W(1,l,iter),sze)
|
||||
h(k,iter,l,iter2) = h(k,iter2,l,iter)
|
||||
enddo
|
||||
enddo
|
||||
do k=1,l
|
||||
h(k,iter,l,iter) = u_dot_v(U(1,k,iter),W(1,l,iter),sze)
|
||||
h(l,iter,k,iter) = h(k,iter,l,iter)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
! Diagonalize h
|
||||
! -------------
|
||||
call lapack_diag(lambda,y,h,N_st*davidson_sze_max,N_st*iter)
|
||||
|
||||
! Express eigenvectors of h in the determinant basis
|
||||
! --------------------------------------------------
|
||||
|
||||
! call dgemm ( 'N','N', sze, N_st*iter, N_st, &
|
||||
! 1.d0, U(1,1,1), size(U,1), y(1,1,1,1), size(y,1)*size(y,2), &
|
||||
! 0.d0, U(1,1,iter+1), size(U,1) )
|
||||
do k=1,N_st
|
||||
do i=1,sze
|
||||
U(i,k,iter+1) = 0.d0
|
||||
W(i,k,iter+1) = 0.d0
|
||||
do l=1,N_st
|
||||
do iter2=1,iter
|
||||
U(i,k,iter+1) = U(i,k,iter+1) + U(i,l,iter2)*y(l,iter2,k,1)
|
||||
W(i,k,iter+1) = W(i,k,iter+1) + W(i,l,iter2)*y(l,iter2,k,1)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
! Compute residual vector
|
||||
! -----------------------
|
||||
|
||||
do k=1,N_st
|
||||
do i=1,sze
|
||||
R(i,k) = lambda(k) * U(i,k,iter+1) - W(i,k,iter+1)
|
||||
enddo
|
||||
residual_norm(k) = u_dot_u(R(1,k),sze)
|
||||
enddo
|
||||
|
||||
print '(I3,15(F16.8,x))', iter, lambda(1:N_st) + nuclear_repulsion
|
||||
print '(3x,15(E16.5,x))', residual_norm(1:N_st)
|
||||
|
||||
converged = maxval(residual_norm) < 1.d-10
|
||||
if (converged) then
|
||||
exit
|
||||
endif
|
||||
|
||||
! Davidson step
|
||||
! -------------
|
||||
|
||||
do k=1,N_st
|
||||
do i=1,sze
|
||||
U(i,k,iter+1) = 1.d0/(lambda(k) - H_jj(i)) * R(i,k)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
! Gram-Schmidt
|
||||
! ------------
|
||||
|
||||
double precision :: c
|
||||
do k=1,N_st
|
||||
do iter2=1,iter
|
||||
do l=1,N_st
|
||||
c = u_dot_v(U(1,k,iter+1),U(1,l,iter2),sze)
|
||||
do i=1,sze
|
||||
U(i,k,iter+1) -= c * U(i,l,iter2)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
do l=1,k-1
|
||||
c = u_dot_v(U(1,k,iter+1),U(1,l,iter+1),sze)
|
||||
do i=1,sze
|
||||
U(i,k,iter+1) -= c * U(i,l,iter+1)
|
||||
enddo
|
||||
enddo
|
||||
call normalize( U(1,k,iter+1), sze )
|
||||
enddo
|
||||
enddo
|
||||
|
||||
if (.not.converged) then
|
||||
iter = davidson_sze_max-1
|
||||
endif
|
||||
|
||||
! Re-contract to u_in
|
||||
! -----------
|
||||
|
||||
do k=1,N_st
|
||||
energies(k) = lambda(k)
|
||||
do i=1,sze
|
||||
u_in(i,k) = 0.d0
|
||||
do iter2=1,iter
|
||||
do l=1,N_st
|
||||
u_in(i,k) += U(i,l,iter2)*y(l,iter2,k,1)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
enddo
|
||||
|
||||
deallocate ( &
|
||||
kl_pairs, &
|
||||
W, &
|
||||
U, &
|
||||
R, &
|
||||
h, &
|
||||
y, &
|
||||
lambda &
|
||||
)
|
||||
end
|
||||
|
||||
subroutine repeat_excitation(key_in,key_1,key_2,i_ok,Nint)
|
||||
use bitmasks
|
||||
implicit none
|
||||
integer(bit_kind), intent(in) :: key_in(Nint,2),key_1(Nint,2),key_2(Nint,2),Nint
|
||||
integer,intent(out):: i_ok
|
||||
integer :: ispin,i_hole,k_hole,j_hole,i_particl,k_particl,j_particl,i_trou,degree,exc(0:2,2,2)
|
||||
double precision :: phase
|
||||
i_ok = 1
|
||||
call get_excitation(key_1,key_2,exc,degree,phase,Nint)
|
||||
integer :: h1,p1,h2,p2,s1,s2
|
||||
if(degree==2)then
|
||||
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
|
||||
! first hole
|
||||
k_hole = ishft(h1-1,-5)+1
|
||||
j_hole = h1-ishft(k_hole-1,5)-1
|
||||
if(iand(key_in(k_hole,s1),ibset(0,j_hole)).eq.0)then
|
||||
i_ok = 0
|
||||
return
|
||||
endif
|
||||
|
||||
! second hole
|
||||
k_hole = ishft(h2-1,-5)+1
|
||||
j_hole = h2-ishft(k_hole-1,5)-1
|
||||
if(iand(key_in(k_hole,s2),ibset(0,j_hole)).eq.0)then
|
||||
i_ok = 0
|
||||
return
|
||||
endif
|
||||
|
||||
! first particle
|
||||
k_particl = ishft(p1-1,-5)+1
|
||||
j_particl = p1-ishft(k_particl-1,5)-1
|
||||
if(iand(key_in(k_particl,s1),ibset(0,j_particl)).ne.0)then
|
||||
i_ok = 0
|
||||
return
|
||||
endif
|
||||
|
||||
! second particle
|
||||
k_particl = ishft(p2-1,-5)+1
|
||||
j_particl = p2-ishft(k_particl-1,5)-1
|
||||
if(iand(key_in(k_particl,s2),ibset(0,j_particl)).ne.0)then
|
||||
i_ok = 0
|
||||
return
|
||||
endif
|
||||
return
|
||||
endif
|
||||
end
|
||||
|
@ -6,6 +6,8 @@ program cisd
|
||||
|
||||
double precision :: pt2(10), norm_pert(10), H_pert_diag
|
||||
|
||||
N_states = 3
|
||||
touch N_states
|
||||
call H_apply_cisd
|
||||
allocate(eigvalues(n_states),eigvectors(n_det,n_states))
|
||||
print *, 'N_det = ', N_det
|
||||
@ -20,7 +22,9 @@ program cisd
|
||||
print *, 'HF:', HF_energy
|
||||
print *, '---'
|
||||
do i = 1,1
|
||||
print *, 'energy(i) = ',eigvalues(i) + nuclear_repulsion
|
||||
print *, 'energy(i) = ',eigvalues(i) + nuclear_repulsion
|
||||
print *, 'E_corr = ',eigvalues(i) - ref_bitmask_energy
|
||||
enddo
|
||||
call CISD_SC2(psi_det,psi_coef,eigvalues,size(psi_coef,1),N_det,N_states,N_int)
|
||||
deallocate(eigvalues,eigvectors)
|
||||
end
|
||||
|
334
src/CIS_dressed/CIS_providers.irp.f
Normal file
334
src/CIS_dressed/CIS_providers.irp.f
Normal file
@ -0,0 +1,334 @@
|
||||
|
||||
use bitmasks
|
||||
BEGIN_PROVIDER [integer(bit_kind), psi_CIS,(N_int,2,size_psi_CIS)]
|
||||
&BEGIN_PROVIDER [integer, psi_CIS_holes,(size_psi_CIS)]
|
||||
&BEGIN_PROVIDER [integer, psi_CIS_particl,(size_psi_CIS)]
|
||||
&BEGIN_PROVIDER [integer, psi_CIS_spin,(size_psi_CIS)]
|
||||
&BEGIN_PROVIDER [integer, psi_CIS_adress,(n_core_cis+1:elec_alpha_num,elec_alpha_num+1:n_act_cis)]
|
||||
&BEGIN_PROVIDER [double precision, H_CIS,(size_psi_CIS,size_psi_CIS)]
|
||||
BEGIN_DOC
|
||||
!key of the CIS-matrix
|
||||
END_DOC
|
||||
|
||||
|
||||
implicit none
|
||||
integer :: a !control variable
|
||||
integer :: i,j,k,l !variables for going over the occupied (i,j) and virutal (k,l)
|
||||
integer :: key !key for CIS-matrix
|
||||
integer :: i_hole,j_hole,ispin,l_particle,k_particle
|
||||
double precision :: hij
|
||||
|
||||
|
||||
do a=1,N_int
|
||||
psi_CIS(a,1,1)=ref_bitmask(a,1)
|
||||
psi_CIS(a,2,1)=ref_bitmask(a,2)
|
||||
enddo
|
||||
|
||||
psi_CIS_holes(1) = 0
|
||||
psi_CIS_particl(1) = 0
|
||||
psi_CIS_spin(1) = 0
|
||||
|
||||
!loop on particles: create a particle in k
|
||||
do k=elec_alpha_num+1,n_act_cis
|
||||
|
||||
!loop on holes: destroy a particle in i
|
||||
do i=n_core_cis+1,elec_alpha_num
|
||||
|
||||
!alpha spin
|
||||
ispin=1
|
||||
|
||||
key=2*((k-elec_alpha_num-1)*(elec_alpha_num-n_core_cis) + i-n_core_cis) !index of such an excited determinant in the CIS WF
|
||||
psi_CIS_adress(i,k)=key
|
||||
|
||||
do a=1,N_int
|
||||
psi_CIS(a,1,key)=ref_bitmask(a,1)
|
||||
psi_CIS(a,2,key)=ref_bitmask(a,2)
|
||||
enddo
|
||||
|
||||
j_hole=ishft(i-1,-5)+1
|
||||
i_hole=i-ishft(j_hole-1,5)-1
|
||||
|
||||
psi_CIS(j_hole,ispin,key)=ibclr(psi_CIS(j_hole,ispin,key),i_hole)
|
||||
|
||||
l_particle=ishft(k-1,-5)+1
|
||||
k_particle=k-ishft(l_particle-1,5)-1
|
||||
|
||||
psi_CIS(l_particle,ispin,key)=ibset(psi_CIS(l_particle,ispin,key),k_particle)
|
||||
|
||||
psi_CIS_holes(key) = i
|
||||
psi_CIS_particl(key) = k
|
||||
psi_CIS_spin(key) = 1
|
||||
|
||||
!beta spin
|
||||
ispin=2
|
||||
|
||||
key=key+1
|
||||
|
||||
do a=1,N_int
|
||||
psi_CIS(a,1,key)=ref_bitmask(a,1)
|
||||
psi_CIS(a,2,key)=ref_bitmask(a,2)
|
||||
enddo
|
||||
|
||||
j_hole=ishft(i-1,-5)+1
|
||||
i_hole=i-ishft(j_hole-1,5)-1
|
||||
|
||||
psi_CIS(j_hole,ispin,key)=ibclr(psi_CIS(j_hole,ispin,key),i_hole)
|
||||
|
||||
l_particle=ishft(k-1,-5)+1
|
||||
k_particle=k-ishft(l_particle-1,5)-1
|
||||
psi_CIS_holes(key) = i
|
||||
psi_CIS_particl(key) = k
|
||||
psi_CIS_spin(key) = 2
|
||||
|
||||
psi_CIS(l_particle,ispin,key)=ibset(psi_CIS(l_particle,ispin,key),k_particle)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
!Building the CIS-matrix
|
||||
double precision :: diag_H_mat_elem
|
||||
|
||||
do key=1,size_psi_CIS
|
||||
H_CIS(key,key)=diag_H_mat_elem(psi_CIS(1,1,key),N_int)
|
||||
|
||||
do a=key+1,size_psi_CIS
|
||||
call i_H_j(psi_CIS(1,1,a),psi_CIS(1,1,key),N_int,hij)
|
||||
|
||||
H_CIS(key,a)=hij
|
||||
H_CIS(a,key)=hij
|
||||
enddo
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
BEGIN_PROVIDER[double precision, eigenvalues_CIS,(n_state_CIS)]
|
||||
&BEGIN_PROVIDER[double precision, coefs_CIS, (size_psi_CIS,n_state_CIS)]
|
||||
use bitmasks
|
||||
|
||||
BEGIN_DOC
|
||||
!the first states of the CIS matrix
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer :: i,j,k
|
||||
double precision :: eigvalues(size_psi_CIS),eigvectors(size_psi_CIS,size_psi_CIS)
|
||||
double precision :: coefs_tmp(size_psi_CIS)
|
||||
double precision :: s2
|
||||
|
||||
|
||||
!Diagonalisation of CIS-matrix
|
||||
call lapack_diag(eigvalues,eigvectors,H_CIS,size_psi_CIS,size_psi_CIS)
|
||||
|
||||
do i = 1,n_state_CIS
|
||||
eigenvalues_CIS(i) = eigvalues(i)
|
||||
|
||||
do k=1,size_psi_CIS
|
||||
|
||||
if (dabs(eigvectors(k,i)).ge.10.d-2) then
|
||||
write(11,*),'k,i,eigenvectors(k,i)=',k,i,eigvectors(k,i)
|
||||
write(11,*),'hole,particl,spin:',psi_CIS_holes(k),psi_CIS_particl(k),psi_CIS_spin(k)
|
||||
write(11,*),''
|
||||
endif
|
||||
coefs_tmp(k) = eigvectors(k,i)
|
||||
coefs_CIS(k,i)=eigvectors(k,i)
|
||||
enddo
|
||||
call get_s2_u0(psi_CIS,coefs_tmp,size_psi_CIS,size_psi_CIS,s2)
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
BEGIN_PROVIDER [double precision, eigenvalues_CIS_dress_D,(n_state_CIS)]
|
||||
&BEGIN_PROVIDER [double precision, s_2_CIS_dress_D,(n_state_CIS)]
|
||||
&BEGIN_PROVIDER [double precision, eigenvectors_CIS_dress_D,(size_psi_CIS,n_state_CIS)]
|
||||
&BEGIN_PROVIDER [double precision, overlap_D ]
|
||||
use bitmasks
|
||||
|
||||
BEGIN_DOC
|
||||
!The first states of the CIS matrix dressed by the doubles
|
||||
END_DOC
|
||||
implicit none
|
||||
double precision,allocatable :: delta_H_matrix_doub(:,:)
|
||||
double precision,allocatable :: eigvalues(:),eigvectors(:,:)
|
||||
double precision :: overlap,max_overlap,s2
|
||||
integer :: i_overlap,i,j,k
|
||||
allocate (delta_H_matrix_doub(size_psi_CIS,size_psi_CIS))
|
||||
allocate(eigvalues(size_psi_CIS),eigvectors(size_psi_CIS,size_psi_CIS))
|
||||
do i = 1,n_state_CIS
|
||||
call dress_by_doubles(eigenvalues_CIS(i),coefs_CIS(1,i),delta_H_matrix_doub,size_psi_CIS) !dressing of the Doubles
|
||||
do j = 1,size_psi_CIS
|
||||
do k = 1,size_psi_CIS
|
||||
delta_H_matrix_doub(j,k) += H_CIS(j,k)
|
||||
enddo
|
||||
enddo
|
||||
call lapack_diag(eigvalues,eigvectors,delta_H_matrix_doub,size_psi_CIS,size_psi_CIS)
|
||||
|
||||
! state following
|
||||
max_overlap = 0.d0
|
||||
do k = 1, size_psi_CIS
|
||||
overlap = 0.d0
|
||||
do j = 1,size_psi_CIS
|
||||
overlap += eigvectors(j,k)*coefs_CIS(j,i)
|
||||
enddo
|
||||
if(dabs(overlap).gt.max_overlap)then
|
||||
max_overlap = dabs(overlap)
|
||||
i_overlap = k
|
||||
endif
|
||||
! <CIS(i)|state(k)>
|
||||
enddo
|
||||
! print*,'overlap = ',max_overlap
|
||||
overlap_D=max_overlap
|
||||
do k = 1,size_psi_CIS
|
||||
eigenvectors_CIS_dress_D(k,i) = eigvectors(k,i_overlap)
|
||||
if (dabs(eigvectors(k,i_overlap)).ge.10.d-2) then
|
||||
write(12,*),'k,i,eigenvectors(k,i)=',k,i,eigvectors(k,i_overlap)
|
||||
write(12,*),'hole,particl,spin:',psi_CIS_holes(k),psi_CIS_particl(k),psi_CIS_spin(k)
|
||||
write(12,*),''
|
||||
endif
|
||||
enddo
|
||||
call get_s2_u0(psi_CIS,eigenvectors_CIS_dress_D(1,i),size_psi_CIS,size_psi_CIS,s2)
|
||||
s_2_CIS_dress_D(i) = s2
|
||||
eigenvalues_CIS_dress_D(i) = eigvalues(i_overlap)
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
|
||||
BEGIN_PROVIDER [double precision, eigenvalues_CIS_dress_D_dt,(n_state_CIS)]
|
||||
&BEGIN_PROVIDER [double precision, s_2_CIS_dress_D_dt,(n_state_CIS)]
|
||||
&BEGIN_PROVIDER [double precision, eigenvectors_CIS_dress_D_dt,(size_psi_CIS,n_state_CIS)
|
||||
&BEGIN_PROVIDER [double precision, overlap_Ddt]
|
||||
use bitmasks
|
||||
|
||||
BEGIN_DOC
|
||||
!The first states of the CIS matrix dressed by the doubles and the disconnected triples
|
||||
END_DOC
|
||||
implicit none
|
||||
double precision,allocatable :: delta_H_matrix_doub(:,:)
|
||||
double precision,allocatable :: eigvalues(:),eigvectors(:,:)
|
||||
double precision :: overlap,max_overlap,s2
|
||||
integer :: i_overlap,i,j,k
|
||||
allocate (delta_H_matrix_doub(size_psi_CIS,size_psi_CIS))
|
||||
allocate(eigvalues(size_psi_CIS),eigvectors(size_psi_CIS,size_psi_CIS))
|
||||
do i = 1,n_state_CIS
|
||||
call dress_by_doubles(eigenvalues_CIS(i),coefs_CIS(1,i),delta_H_matrix_doub,size_psi_CIS) !dressing of the Doubles
|
||||
|
||||
do j = 1,size_psi_CIS
|
||||
do k = 1,size_psi_CIS
|
||||
delta_H_matrix_doub(j,k) += H_CIS(j,k)
|
||||
enddo
|
||||
delta_H_matrix_doub(j,j) += dress_T_discon_array_CIS(j)
|
||||
enddo
|
||||
call lapack_diag(eigvalues,eigvectors,delta_H_matrix_doub,size_psi_CIS,size_psi_CIS)
|
||||
|
||||
! state following
|
||||
max_overlap = 0.d0
|
||||
do k = 1, size_psi_CIS
|
||||
overlap = 0.d0
|
||||
do j = 1,size_psi_CIS
|
||||
overlap += eigvectors(j,k)*coefs_CIS(j,i)
|
||||
enddo
|
||||
if(dabs(overlap).gt.max_overlap)then
|
||||
max_overlap = dabs(overlap)
|
||||
i_overlap = k
|
||||
endif
|
||||
! <CIS(i)|state(k)>
|
||||
enddo
|
||||
! print*,'overlap = ',max_overlap
|
||||
overlap_Ddt=max_overlap
|
||||
do k = 1,size_psi_CIS
|
||||
eigenvectors_CIS_dress_D_dt(k,i) = eigvectors(k,i_overlap)
|
||||
if (dabs(eigvectors(k,i_overlap)).ge.10.d-2) then
|
||||
write(13,*),'k,i,eigenvectors(k,i)=',k,i,eigvectors(k,i_overlap)
|
||||
write(13,*),'hole,particl,spin:',psi_CIS_holes(k),psi_CIS_particl(k),psi_CIS_spin(k)
|
||||
write(13,*),''
|
||||
endif
|
||||
enddo
|
||||
call get_s2_u0(psi_CIS,eigenvectors_CIS_dress_D_dt(1,i),size_psi_CIS,size_psi_CIS,s2)
|
||||
s_2_CIS_dress_D_dt(i) = s2
|
||||
eigenvalues_CIS_dress_D_dt(i) = eigvalues(i_overlap)
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
!BEGIN_PROVIDER [double precision, eigenvalues_CIS_dress_tot,(n_state_CIS)]
|
||||
!BEGIN_PROVIDER [double precision, s_2_CIS_dress_tot,(n_state_CIS)]
|
||||
!BEGIN_PROVIDER [double precision, eigenvectors_CIS_dress_tot,(size_psi_CIS,n_state_CIS)]
|
||||
!BEGIN_PROVIDER [double precision, overlap_tot]
|
||||
|
||||
!BEGIN_DOC
|
||||
!!The first states of the CIS matrix dressed by the doubles
|
||||
!END_DOC
|
||||
!implicit none
|
||||
!double precision,allocatable :: delta_H_matrix_doub(:,:)
|
||||
!double precision,allocatable :: eigvalues(:),eigvectors(:,:)
|
||||
!double precision,allocatable :: delta_H_trip(:,:)
|
||||
!double precision :: overlap,max_overlap,s2,average_eigvalue
|
||||
!integer :: i_overlap,i,j,k,m
|
||||
!allocate (delta_H_matrix_doub(size_psi_CIS,size_psi_CIS),delta_H_trip(size_psi_CIS,size_psi_CIS) )
|
||||
!allocate(eigvalues(size_psi_CIS),eigvectors(size_psi_CIS,size_psi_CIS))
|
||||
! do i = 1,n_state_CIS
|
||||
! call dress_by_doubles(eigenvalues_CIS(i),coefs_CIS(1,i),delta_H_matrix_doub,size_psi_CIS) !dressing of the Doubles
|
||||
! call dress_T_con(eigenvalues_CIS(i),delta_H_trip,size_psi_CIS)
|
||||
! do j = 1,size_psi_CIS
|
||||
! do k = 1,size_psi_CIS
|
||||
! delta_H_matrix_doub(j,k) += H_CIS(j,k)
|
||||
! delta_H_matrix_doub(j,k) += delta_H_trip(j,k)
|
||||
! enddo
|
||||
! delta_H_matrix_doub(j,j) += dress_T_discon_array_CIS(j)
|
||||
! enddo
|
||||
! call lapack_diag(eigvalues,eigvectors,delta_H_matrix_doub,size_psi_CIS,size_psi_CIS)
|
||||
! do m=1,n_state_CIS
|
||||
! write(12,*),'m,eigvalues(m)',m,eigvalues(m)
|
||||
! enddo
|
||||
! ! state following
|
||||
! max_overlap = 0.d0
|
||||
! do k = 1, size_psi_CIS
|
||||
! overlap = 0.d0
|
||||
! do j = 1,size_psi_CIS
|
||||
! overlap += eigvectors(j,k)*coefs_CIS(j,i)
|
||||
! enddo
|
||||
! if(dabs(overlap).gt.max_overlap)then
|
||||
! max_overlap = dabs(overlap)
|
||||
! i_overlap = k
|
||||
! endif
|
||||
! ! <CIS(i)|state(k)>
|
||||
! enddo
|
||||
!! print*,'overlap = ',max_overlap
|
||||
! overlap_tot=max_overlap
|
||||
! do k = 1,size_psi_CIS
|
||||
! eigenvectors_CIS_dress_tot(k,i) = eigvectors(k,i_overlap)
|
||||
! enddo
|
||||
! call get_s2_u0(psi_CIS,eigenvectors_CIS_dress_tot(1,i),size_psi_CIS,size_psi_CIS,s2)
|
||||
! s_2_CIS_dress_tot(i) = s2
|
||||
! eigenvalues_CIS_dress_tot(i) = eigvalues(i_overlap)
|
||||
! enddo
|
||||
|
||||
!END_PROVIDER
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
BEGIN_PROVIDER [double precision, diag_elements, (size_psi_CIS)]
|
||||
use bitmasks
|
||||
|
||||
|
||||
BEGIN_DOC
|
||||
!Array of the energy of the CIS determinants ordered in the CIS matrix
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
double precision :: hij
|
||||
integer :: i
|
||||
|
||||
do i = 1, size_psi_CIS
|
||||
call i_H_j(psi_CIS(1,1,i),psi_CIS(1,1,i),N_int,hij)
|
||||
diag_elements(i) = hij
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
1611
src/CIS_dressed/MP2.irp.f
Normal file
1611
src/CIS_dressed/MP2.irp.f
Normal file
File diff suppressed because it is too large
Load Diff
28
src/CIS_dressed/README.rst
Normal file
28
src/CIS_dressed/README.rst
Normal file
@ -0,0 +1,28 @@
|
||||
Needed Modules
|
||||
==============
|
||||
|
||||
.. Do not edit this section. It was auto-generated from the
|
||||
.. NEEDED_MODULES file.
|
||||
|
||||
* `AOs <http://github.com/LCPQ/quantum_package/tree/master/src/AOs>`_
|
||||
* `BiInts <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts>`_
|
||||
* `Bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask>`_
|
||||
* `Dets <http://github.com/LCPQ/quantum_package/tree/master/src/Dets>`_
|
||||
* `Electrons <http://github.com/LCPQ/quantum_package/tree/master/src/Electrons>`_
|
||||
* `Ezfio_files <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files>`_
|
||||
* `Hartree_Fock <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock>`_
|
||||
* `MonoInts <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts>`_
|
||||
* `MOs <http://github.com/LCPQ/quantum_package/tree/master/src/MOs>`_
|
||||
* `Nuclei <http://github.com/LCPQ/quantum_package/tree/master/src/Nuclei>`_
|
||||
* `Output <http://github.com/LCPQ/quantum_package/tree/master/src/Output>`_
|
||||
* `Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Utils>`_
|
||||
* `DensityMatrix <http://github.com/LCPQ/quantum_package/tree/master/src/DensityMatrix>`_
|
||||
|
||||
Documentation
|
||||
=============
|
||||
|
||||
.. Do not edit this section. It was auto-generated from the
|
||||
.. NEEDED_MODULES file.
|
||||
|
||||
|
||||
|
6
src/CIS_dressed/cis_dressed.ezfio_config
Normal file
6
src/CIS_dressed/cis_dressed.ezfio_config
Normal file
@ -0,0 +1,6 @@
|
||||
cis_dressed
|
||||
n_state_cis integer
|
||||
n_core_cis integer
|
||||
n_act_cis integer
|
||||
mp2_dressing logical
|
||||
standard_doubles logical
|
66
src/CIS_dressed/density_matrix_suroutine.irp.f
Normal file
66
src/CIS_dressed/density_matrix_suroutine.irp.f
Normal file
@ -0,0 +1,66 @@
|
||||
|
||||
subroutine get_dm_from_psi(dets_in,u_in,sze,dim_in,Nint,dm_alpha,dm_beta)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Alpha and beta one-body density matrix
|
||||
!
|
||||
! dets_in :: bitsrings corresponding to the determinants in the wave function
|
||||
!
|
||||
! u_in :: coefficients of the wave function
|
||||
!
|
||||
! sze :: number of determinants in the wave function
|
||||
!
|
||||
! dim_in :: physical dimension of the array u_in and dets_in
|
||||
!
|
||||
! Nint :: should be equal to N_int
|
||||
!
|
||||
! dm_alpha :: alpha one body density matrix
|
||||
!
|
||||
! dm_beta :: beta one body density matrix
|
||||
END_DOC
|
||||
use bitmasks
|
||||
integer, intent(in) :: sze,dim_in,Nint
|
||||
integer(bit_kind), intent(in) :: dets_in(Nint,2,dim_in)
|
||||
double precision, intent(in) :: u_in(dim_in)
|
||||
double precision, intent(out) :: dm_alpha(mo_tot_num,mo_tot_num)
|
||||
double precision, intent(out) :: dm_beta(mo_tot_num,mo_tot_num)
|
||||
|
||||
integer :: j,k,l
|
||||
integer :: occ(N_int*bit_kind_size,2)
|
||||
double precision :: ck, cl, ckl
|
||||
double precision :: phase
|
||||
integer :: h1,h2,p1,p2,s1,s2, degree
|
||||
integer :: exc(0:2,2,2),n_occ_alpha
|
||||
dm_alpha = 0.d0
|
||||
dm_beta = 0.d0
|
||||
|
||||
do k=1,sze
|
||||
call bitstring_to_list(dets_in(1,1,k), occ(1,1), n_occ_alpha, N_int)
|
||||
call bitstring_to_list(dets_in(1,2,k), occ(1,2), n_occ_alpha, N_int)
|
||||
ck = u_in(k)
|
||||
do l=1,elec_alpha_num
|
||||
j = occ(l,1)
|
||||
dm_alpha(j,j) += ck*ck
|
||||
enddo
|
||||
do l=1,elec_beta_num
|
||||
j = occ(l,2)
|
||||
dm_beta(j,j) += ck*ck
|
||||
enddo
|
||||
do l=1,k-1
|
||||
call get_excitation_degree(dets_in(1,1,k),dets_in(1,1,l),degree,N_int)
|
||||
if (degree /= 1) then
|
||||
cycle
|
||||
endif
|
||||
call get_mono_excitation(dets_in(1,1,k),dets_in(1,1,l),exc,phase,N_int)
|
||||
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
|
||||
ckl = ck * u_in(l) * phase
|
||||
if (s1==1) then
|
||||
dm_alpha(h1,p1) += ckl
|
||||
dm_alpha(p1,h1) += ckl
|
||||
else
|
||||
dm_beta(h1,p1) += ckl
|
||||
dm_beta(p1,h1) += ckl
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
end
|
92
src/CIS_dressed/natural_particl_orbitals.irp.f
Normal file
92
src/CIS_dressed/natural_particl_orbitals.irp.f
Normal file
@ -0,0 +1,92 @@
|
||||
BEGIN_PROVIDER [double precision, particle_natural_orb_CIS_properties, (6,n_state_cis)]
|
||||
&BEGIN_PROVIDER [double precision, CIS_states_properties, (6,n_state_cis)]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! properties of the natural orbital of the particle of the various n_state_cis eigenvectors of the CIS matrix
|
||||
!
|
||||
! You first build the density matrix of the one eigenvector and you take off the Hartree Fock density matrix
|
||||
!
|
||||
! particl(i,j)(state = k) == dm(i,j)(Hartree Fock) - dm(i,j)(state = k)
|
||||
!
|
||||
! you diagonalize particl(i,j) and the first eigenvector is the natural orbital corresponding to the particl
|
||||
!
|
||||
! that is specific to the excitation in the CIS state
|
||||
!
|
||||
! particle_natural_orb_CIS_properties(i,1) = <phi_i|x|phi_i>
|
||||
!
|
||||
! particle_natural_orb_CIS_properties(i,2) = <phi_i|y|phi_i>
|
||||
!
|
||||
! particle_natural_orb_CIS_properties(i,3) = <phi_i|z|phi_i>
|
||||
!
|
||||
! particle_natural_orb_CIS_properties(i,5) = <phi_i|x^2|phi_i>
|
||||
!
|
||||
! particle_natural_orb_CIS_properties(i,6) = <phi_i|y^2|phi_i>
|
||||
!
|
||||
! particle_natural_orb_CIS_properties(i,7) = <phi_i|z^2|phi_i>
|
||||
!
|
||||
! CIS_states_properties(i,1:6) = the same but for the hole state i
|
||||
END_DOC
|
||||
|
||||
integer :: i,j,k,l
|
||||
double precision :: dm_alpha(mo_tot_num,mo_tot_num)
|
||||
double precision :: dm_beta(mo_tot_num,mo_tot_num)
|
||||
double precision :: dm(mo_tot_num,mo_tot_num)
|
||||
double precision :: eigvalues(mo_tot_num)
|
||||
double precision :: eigvectors(mo_tot_num,mo_tot_num)
|
||||
double precision :: accu_n_elec,c_k
|
||||
|
||||
|
||||
do i = 1, n_state_cis
|
||||
print*,' state cis = ',i
|
||||
call get_dm_from_psi(psi_CIS,coefs_CIS(1,i),size_psi_CIS,size_psi_CIS,N_int,dm_alpha,dm_beta)
|
||||
|
||||
dm = dm_alpha + dm_beta
|
||||
call get_properties_from_density_matrix(dm,CIS_states_properties(1,i))
|
||||
dm = -dm
|
||||
do k = 1, elec_alpha_num
|
||||
dm(k,k) += 1.d0
|
||||
enddo
|
||||
do k = 1, elec_beta_num
|
||||
dm(k,k) += 1.d0
|
||||
enddo
|
||||
call lapack_diag(eigvalues,eigvectors,dm,mo_tot_num,mo_tot_num)
|
||||
accu_n_elec = 0.d0
|
||||
do k = 1, mo_tot_num
|
||||
accu_n_elec += eigvalues(k)
|
||||
enddo
|
||||
do k = 1, mo_tot_num
|
||||
do l = 1, mo_tot_num
|
||||
c_k = eigvectors(k,j) * eigvectors(l,j)
|
||||
particle_natural_orb_CIS_properties(1,i) += c_k * mo_dipole_x(k,l)
|
||||
particle_natural_orb_CIS_properties(2,i) += c_k * mo_dipole_y(k,l)
|
||||
particle_natural_orb_CIS_properties(3,i) += c_k * mo_dipole_z(k,l)
|
||||
particle_natural_orb_CIS_properties(4,i) += c_k * mo_spread_x(k,l)
|
||||
particle_natural_orb_CIS_properties(5,i) += c_k * mo_spread_y(k,l)
|
||||
particle_natural_orb_CIS_properties(6,i) += c_k * mo_spread_z(k,l)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
subroutine get_properties_from_density_matrix(dm,properties)
|
||||
implicit none
|
||||
double precision, intent(in) :: dm(mo_tot_num,mo_tot_num)
|
||||
double precision, intent(out) :: properties(6)
|
||||
integer :: k,l
|
||||
double precision :: c_k
|
||||
do k = 1, 6
|
||||
properties(k) = 0.d0
|
||||
enddo
|
||||
do k = 1, mo_tot_num
|
||||
do l = 1, mo_tot_num
|
||||
c_k = dm(k,l)
|
||||
properties(1) += c_k * mo_dipole_x(k,l)
|
||||
properties(2) += c_k * mo_dipole_y(k,l)
|
||||
properties(3) += c_k * mo_dipole_z(k,l)
|
||||
properties(4) += c_k * mo_spread_x(k,l)
|
||||
properties(5) += c_k * mo_spread_y(k,l)
|
||||
properties(6) += c_k * mo_spread_z(k,l)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
end
|
90
src/CIS_dressed/options.irp.f
Normal file
90
src/CIS_dressed/options.irp.f
Normal file
@ -0,0 +1,90 @@
|
||||
BEGIN_PROVIDER [ integer , n_state_cis ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Number of states asked for the CIS vector
|
||||
END_DOC
|
||||
|
||||
logical :: has
|
||||
PROVIDE ezfio_filename
|
||||
call ezfio_has_cis_dressed_n_state_cis(has)
|
||||
if (has) then
|
||||
call ezfio_get_cis_dressed_n_state_cis(n_state_cis)
|
||||
else
|
||||
n_state_cis = 5
|
||||
call ezfio_set_cis_dressed_n_state_cis(n_state_cis)
|
||||
endif
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ integer , n_core_cis]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Number of states asked for the CIS vector
|
||||
END_DOC
|
||||
|
||||
logical :: has
|
||||
PROVIDE ezfio_filename
|
||||
call ezfio_has_cis_dressed_n_core_cis(has)
|
||||
if (has) then
|
||||
call ezfio_get_cis_dressed_n_core_cis(n_core_cis)
|
||||
else
|
||||
n_core_cis = 0
|
||||
call ezfio_set_cis_dressed_n_core_cis(n_core_cis)
|
||||
endif
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ integer , n_act_cis]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Number of states asked for the CIS vector
|
||||
END_DOC
|
||||
|
||||
logical :: has
|
||||
PROVIDE ezfio_filename
|
||||
call ezfio_has_cis_dressed_n_act_cis(has)
|
||||
if (has) then
|
||||
call ezfio_get_cis_dressed_n_act_cis(n_act_cis)
|
||||
else
|
||||
n_act_cis = mo_tot_num
|
||||
call ezfio_set_cis_dressed_n_act_cis(n_act_cis)
|
||||
endif
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ logical , mp2_dressing]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Number of states asked for the CIS vector
|
||||
END_DOC
|
||||
|
||||
logical :: has
|
||||
PROVIDE ezfio_filename
|
||||
call ezfio_has_cis_dressed_mp2_dressing(has)
|
||||
if (has) then
|
||||
call ezfio_get_cis_dressed_mp2_dressing(mp2_dressing)
|
||||
else
|
||||
mp2_dressing = .False.
|
||||
call ezfio_set_cis_dressed_mp2_dressing(mp2_dressing)
|
||||
endif
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ logical , standard_doubles]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Number of states asked for the CIS vector
|
||||
END_DOC
|
||||
|
||||
logical :: has
|
||||
PROVIDE ezfio_filename
|
||||
call ezfio_has_cis_dressed_standard_doubles(has)
|
||||
if (has) then
|
||||
call ezfio_get_cis_dressed_standard_doubles(standard_doubles)
|
||||
else
|
||||
standard_doubles = .True.
|
||||
call ezfio_set_cis_dressed_standard_doubles(standard_doubles)
|
||||
endif
|
||||
|
||||
END_PROVIDER
|
||||
|
@ -8,6 +8,39 @@ Documentation
|
||||
.. Do not edit this section. It was auto-generated from the
|
||||
.. NEEDED_MODULES file.
|
||||
|
||||
`iunit_two_body_dm_aa <http://github.com/LCPQ/quantum_package/tree/master/src/DensityMatrix/density_matrix.irp.f#L/BEGIN_PROVIDER [ integer, iunit_two_body_dm_aa ]/;">`_
|
||||
Temporary files for 2-body dm calculation
|
||||
|
||||
`iunit_two_body_dm_ab <http://github.com/LCPQ/quantum_package/tree/master/src/DensityMatrix/density_matrix.irp.f#L/&BEGIN_PROVIDER [ integer, iunit_two_body_dm_ab ]/;">`_
|
||||
Temporary files for 2-body dm calculation
|
||||
|
||||
`iunit_two_body_dm_bb <http://github.com/LCPQ/quantum_package/tree/master/src/DensityMatrix/density_matrix.irp.f#L/&BEGIN_PROVIDER [ integer, iunit_two_body_dm_bb ]/;">`_
|
||||
Temporary files for 2-body dm calculation
|
||||
|
||||
`one_body_dm_a <http://github.com/LCPQ/quantum_package/tree/master/src/DensityMatrix/density_matrix.irp.f#L/BEGIN_PROVIDER [ double precision, one_body_dm_a, (mo_tot_num_align,mo_tot_num) ]/;">`_
|
||||
Alpha and beta one-body density matrix
|
||||
|
||||
`one_body_dm_b <http://github.com/LCPQ/quantum_package/tree/master/src/DensityMatrix/density_matrix.irp.f#L/&BEGIN_PROVIDER [ double precision, one_body_dm_b, (mo_tot_num_align,mo_tot_num) ]/;">`_
|
||||
Alpha and beta one-body density matrix
|
||||
|
||||
`two_body_dm_diag_aa <http://github.com/LCPQ/quantum_package/tree/master/src/DensityMatrix/density_matrix.irp.f#L/BEGIN_PROVIDER [ double precision, two_body_dm_diag_aa, (mo_tot_num_align,mo_tot_num)]/;">`_
|
||||
diagonal part of the two body density matrix
|
||||
|
||||
`two_body_dm_diag_ab <http://github.com/LCPQ/quantum_package/tree/master/src/DensityMatrix/density_matrix.irp.f#L/&BEGIN_PROVIDER [ double precision, two_body_dm_diag_ab, (mo_tot_num_align,mo_tot_num)]/;">`_
|
||||
diagonal part of the two body density matrix
|
||||
|
||||
`two_body_dm_diag_bb <http://github.com/LCPQ/quantum_package/tree/master/src/DensityMatrix/density_matrix.irp.f#L/&BEGIN_PROVIDER [ double precision, two_body_dm_diag_bb, (mo_tot_num_align,mo_tot_num)]/;">`_
|
||||
diagonal part of the two body density matrix
|
||||
|
||||
`det_coef_provider <http://github.com/LCPQ/quantum_package/tree/master/src/DensityMatrix/det_num.irp.f#L/&BEGIN_PROVIDER [ double precision , det_coef_provider, (det_num) ]/;">`_
|
||||
Undocumented
|
||||
|
||||
`det_num <http://github.com/LCPQ/quantum_package/tree/master/src/DensityMatrix/det_num.irp.f#L/BEGIN_PROVIDER [integer, det_num]/;">`_
|
||||
Undocumented
|
||||
|
||||
`det_provider <http://github.com/LCPQ/quantum_package/tree/master/src/DensityMatrix/det_num.irp.f#L/BEGIN_PROVIDER [ integer(bit_kind), det_provider, (N_int,2,det_num)]/;">`_
|
||||
Undocumented
|
||||
|
||||
|
||||
|
||||
Needed Modules
|
||||
@ -28,4 +61,5 @@ Needed Modules
|
||||
* `Nuclei <http://github.com/LCPQ/quantum_package/tree/master/src/Nuclei>`_
|
||||
* `Output <http://github.com/LCPQ/quantum_package/tree/master/src/Output>`_
|
||||
* `Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Utils>`_
|
||||
* `DensityMatrix <http://github.com/LCPQ/quantum_package/tree/master/src/DensityMatrix>`_
|
||||
|
||||
|
@ -43,14 +43,14 @@
|
||||
-0.316337207748718D-01 &
|
||||
/)
|
||||
|
||||
do i=1,10
|
||||
call write_bitstring( 6, det_provider(1,1,i), N_int )
|
||||
enddo
|
||||
print *, ''
|
||||
do i=1,10
|
||||
call write_bitstring( 6, det_provider(1,2,i), N_int )
|
||||
enddo
|
||||
print *, ''
|
||||
!do i=1,10
|
||||
! call write_bitstring( 6, det_provider(1,1,i), N_int )
|
||||
!enddo
|
||||
!print *, ''
|
||||
!do i=1,10
|
||||
! call write_bitstring( 6, det_provider(1,2,i), N_int )
|
||||
!enddo
|
||||
!print *, ''
|
||||
|
||||
|
||||
END_PROVIDER
|
||||
|
@ -275,3 +275,4 @@ subroutine davidson_diag(dets_in,u_in,energies,dim_in,sze,N_st,Nint)
|
||||
)
|
||||
end
|
||||
|
||||
|
||||
|
@ -1,11 +1,12 @@
|
||||
subroutine get_s2(key_i,key_j,phase,Nint)
|
||||
implicit none
|
||||
use bitmasks
|
||||
BEGIN_DOC
|
||||
! Returns <S^2>
|
||||
END_DOC
|
||||
integer, intent(in) :: Nint
|
||||
integer, intent(in) :: key_i(Nint,2)
|
||||
integer, intent(in) :: key_j(Nint,2)
|
||||
integer(bit_kind), intent(in) :: key_i(Nint,2)
|
||||
integer(bit_kind), intent(in) :: key_j(Nint,2)
|
||||
double precision, intent(out) :: phase
|
||||
integer :: exc(0:2,2,2)
|
||||
integer :: degree
|
||||
@ -32,3 +33,38 @@ subroutine get_s2(key_i,key_j,phase,Nint)
|
||||
end select
|
||||
end
|
||||
|
||||
BEGIN_PROVIDER [ double precision, S_z ]
|
||||
&BEGIN_PROVIDER [ double precision, S_z2_Sz ]
|
||||
implicit none
|
||||
|
||||
S_z = 0.5d0*dble(elec_alpha_num-elec_beta_num)
|
||||
S_z2_Sz = S_z*(S_z-1.d0)
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
subroutine get_s2_u0(psi_keys_tmp,psi_coefs_tmp,n,nmax,s2)
|
||||
implicit none
|
||||
use bitmasks
|
||||
integer(bit_kind), intent(in) :: psi_keys_tmp(N_int,2,nmax)
|
||||
integer, intent(in) :: n,nmax
|
||||
double precision, intent(in) :: psi_coefs_tmp(nmax)
|
||||
double precision, intent(out) :: s2
|
||||
integer :: i,j,l
|
||||
double precision :: s2_tmp
|
||||
s2 = S_z2_Sz
|
||||
!$OMP PARALLEL DO DEFAULT(NONE) &
|
||||
!$OMP PRIVATE(i,j,s2_tmp) SHARED(n,psi_coefs_tmp,psi_keys_tmp,N_int) &
|
||||
!$OMP REDUCTION(+:s2) SCHEDULE(dynamic)
|
||||
do i = 1, n
|
||||
call get_s2(psi_keys_tmp(1,1,i),psi_keys_tmp(1,1,i),s2_tmp,N_int)
|
||||
! print*,'s2_tmp = ',s2_tmp
|
||||
do j = 1, n
|
||||
call get_s2(psi_keys_tmp(1,1,i),psi_keys_tmp(1,1,j),s2_tmp,N_int)
|
||||
if (s2_tmp == 0.d0) cycle
|
||||
s2 += psi_coefs_tmp(i)*psi_coefs_tmp(j)*s2_tmp
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END PARALLEL DO
|
||||
end
|
||||
|
||||
|
@ -21,19 +21,19 @@ Documentation
|
||||
.. Do not edit this section. It was auto-generated from the
|
||||
.. NEEDED_MODULES file.
|
||||
|
||||
`fock_matrix_alpha_ao <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/Fock_matrix.irp.f#L83>`_
|
||||
`fock_matrix_alpha_ao <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/Fock_matrix.irp.f#L/BEGIN_PROVIDER [ double precision, Fock_matrix_alpha_ao, (ao_num_align, ao_num) ]/;">`_
|
||||
Alpha Fock matrix in AO basis set
|
||||
|
||||
`fock_matrix_alpha_mo <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/Fock_matrix.irp.f#L172>`_
|
||||
`fock_matrix_alpha_mo <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/Fock_matrix.irp.f#L/BEGIN_PROVIDER [ double precision, Fock_matrix_alpha_mo, (mo_tot_num_align,mo_tot_num) ]/;">`_
|
||||
Fock matrix on the MO basis
|
||||
|
||||
`fock_matrix_beta_ao <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/Fock_matrix.irp.f#L84>`_
|
||||
`fock_matrix_beta_ao <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/Fock_matrix.irp.f#L/&BEGIN_PROVIDER [ double precision, Fock_matrix_beta_ao, (ao_num_align, ao_num) ]/;">`_
|
||||
Alpha Fock matrix in AO basis set
|
||||
|
||||
`fock_matrix_beta_mo <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/Fock_matrix.irp.f#L192>`_
|
||||
`fock_matrix_beta_mo <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/Fock_matrix.irp.f#L/BEGIN_PROVIDER [ double precision, Fock_matrix_beta_mo, (mo_tot_num_align,mo_tot_num) ]/;">`_
|
||||
Fock matrix on the MO basis
|
||||
|
||||
`fock_matrix_diag_mo <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/Fock_matrix.irp.f#L2>`_
|
||||
`fock_matrix_diag_mo <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/Fock_matrix.irp.f#L/&BEGIN_PROVIDER [ double precision, Fock_matrix_diag_mo, (mo_tot_num)]/;">`_
|
||||
Fock matrix on the MO basis.
|
||||
For open shells, the ROHF Fock Matrix is
|
||||
.br
|
||||
@ -48,7 +48,7 @@ Documentation
|
||||
K = Fb - Fa
|
||||
.br
|
||||
|
||||
`fock_matrix_mo <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/Fock_matrix.irp.f#L1>`_
|
||||
`fock_matrix_mo <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/Fock_matrix.irp.f#L/BEGIN_PROVIDER [ double precision, Fock_matrix_mo, (mo_tot_num_align,mo_tot_num) ]/;">`_
|
||||
Fock matrix on the MO basis.
|
||||
For open shells, the ROHF Fock Matrix is
|
||||
.br
|
||||
@ -63,49 +63,49 @@ Documentation
|
||||
K = Fb - Fa
|
||||
.br
|
||||
|
||||
`hf_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/Fock_matrix.irp.f#L211>`_
|
||||
`hf_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/Fock_matrix.irp.f#L/BEGIN_PROVIDER [ double precision, HF_energy ]/;">`_
|
||||
Hartree-Fock energy
|
||||
|
||||
`hf_density_matrix_ao <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/HF_density_matrix_ao.irp.f#L46>`_
|
||||
`hf_density_matrix_ao <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/HF_density_matrix_ao.irp.f#L/BEGIN_PROVIDER [ double precision, HF_density_matrix_ao, (ao_num_align,ao_num) ]/;">`_
|
||||
Density matrix in the AO basis
|
||||
|
||||
`hf_density_matrix_ao_alpha <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/HF_density_matrix_ao.irp.f#L1>`_
|
||||
`hf_density_matrix_ao_alpha <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/HF_density_matrix_ao.irp.f#L/BEGIN_PROVIDER [ double precision, HF_density_matrix_ao_alpha, (ao_num_align,ao_num) ]/;">`_
|
||||
Alpha and Beta density matrix in the AO basis
|
||||
|
||||
`hf_density_matrix_ao_beta <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/HF_density_matrix_ao.irp.f#L2>`_
|
||||
`hf_density_matrix_ao_beta <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/HF_density_matrix_ao.irp.f#L/&BEGIN_PROVIDER [ double precision, HF_density_matrix_ao_beta, (ao_num_align,ao_num) ]/;">`_
|
||||
Alpha and Beta density matrix in the AO basis
|
||||
|
||||
`diagonal_fock_matrix_mo <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/diagonalize_fock.irp.f#L1>`_
|
||||
`diagonal_fock_matrix_mo <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/diagonalize_fock.irp.f#L/BEGIN_PROVIDER [ double precision, diagonal_Fock_matrix_mo, (mo_tot_num) ]/;">`_
|
||||
Diagonal Fock matrix in the MO basis
|
||||
|
||||
`eigenvectors_fock_matrix_mo <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/diagonalize_fock.irp.f#L2>`_
|
||||
`eigenvectors_fock_matrix_mo <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/diagonalize_fock.irp.f#L/&BEGIN_PROVIDER [ double precision, eigenvectors_Fock_matrix_mo, (ao_num_align,mo_tot_num) ]/;">`_
|
||||
Diagonal Fock matrix in the MO basis
|
||||
|
||||
`scf_iteration <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/mo_SCF_iterations.irp.f#L1>`_
|
||||
`scf_iteration <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/mo_SCF_iterations.irp.f#L/subroutine scf_iteration/;">`_
|
||||
Undocumented
|
||||
|
||||
`do_diis <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/options.irp.f#L41>`_
|
||||
`do_diis <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/options.irp.f#L/BEGIN_PROVIDER [ logical, do_DIIS ]/;">`_
|
||||
If True, compute integrals on the fly
|
||||
|
||||
`n_it_scf_max <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/options.irp.f#L22>`_
|
||||
`n_it_scf_max <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/options.irp.f#L/BEGIN_PROVIDER [ integer, n_it_scf_max]/;">`_
|
||||
Maximum number of SCF iterations
|
||||
|
||||
`thresh_scf <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/options.irp.f#L1>`_
|
||||
`thresh_scf <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/options.irp.f#L/BEGIN_PROVIDER [ double precision,thresh_SCF ]/;">`_
|
||||
Threshold on the convergence of the Hartree Fock energy
|
||||
|
||||
`bi_elec_ref_bitmask_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/ref_bitmask.irp.f#L5>`_
|
||||
`bi_elec_ref_bitmask_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/ref_bitmask.irp.f#L/&BEGIN_PROVIDER [ double precision, bi_elec_ref_bitmask_energy ]/;">`_
|
||||
Energy of the reference bitmask used in Slater rules
|
||||
|
||||
`kinetic_ref_bitmask_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/ref_bitmask.irp.f#L3>`_
|
||||
`kinetic_ref_bitmask_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/ref_bitmask.irp.f#L/&BEGIN_PROVIDER [ double precision, kinetic_ref_bitmask_energy ]/;">`_
|
||||
Energy of the reference bitmask used in Slater rules
|
||||
|
||||
`mono_elec_ref_bitmask_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/ref_bitmask.irp.f#L2>`_
|
||||
`mono_elec_ref_bitmask_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/ref_bitmask.irp.f#L/&BEGIN_PROVIDER [ double precision, mono_elec_ref_bitmask_energy ]/;">`_
|
||||
Energy of the reference bitmask used in Slater rules
|
||||
|
||||
`nucl_elec_ref_bitmask_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/ref_bitmask.irp.f#L4>`_
|
||||
`nucl_elec_ref_bitmask_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/ref_bitmask.irp.f#L/&BEGIN_PROVIDER [ double precision, nucl_elec_ref_bitmask_energy ]/;">`_
|
||||
Energy of the reference bitmask used in Slater rules
|
||||
|
||||
`ref_bitmask_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/ref_bitmask.irp.f#L1>`_
|
||||
`ref_bitmask_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/ref_bitmask.irp.f#L/BEGIN_PROVIDER [ double precision, ref_bitmask_energy ]/;">`_
|
||||
Energy of the reference bitmask used in Slater rules
|
||||
|
||||
|
||||
|
@ -1,29 +0,0 @@
|
||||
OPENMP =1
|
||||
PROFILE =0
|
||||
DEBUG = 0
|
||||
|
||||
IRPF90_FLAGS+= --align=32
|
||||
FC = ifort -g
|
||||
FCFLAGS=
|
||||
FCFLAGS+= -xHost
|
||||
#FCFLAGS+= -xAVX
|
||||
FCFLAGS+= -O2
|
||||
FCFLAGS+= -ip
|
||||
FCFLAGS+= -opt-prefetch
|
||||
FCFLAGS+= -ftz
|
||||
MKL=-mkl=parallel
|
||||
|
||||
ifeq ($(PROFILE),1)
|
||||
FC += -p -g
|
||||
CXX += -pg
|
||||
endif
|
||||
|
||||
ifeq ($(OPENMP),1)
|
||||
FC += -openmp
|
||||
CXX += -fopenmp
|
||||
endif
|
||||
|
||||
ifeq ($(DEBUG),1)
|
||||
FC += -C -traceback -fpe0
|
||||
#FCFLAGS =-O0
|
||||
endif
|
@ -17,96 +17,115 @@ Documentation
|
||||
.. Do not edit this section. It was auto-generated from the
|
||||
.. NEEDED_MODULES file.
|
||||
|
||||
`ao_mono_elec_integral <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/ao_mono_ints.irp.f#L122>`_
|
||||
`ao_mono_elec_integral <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/ao_mono_ints.irp.f#L/BEGIN_PROVIDER [ double precision, ao_mono_elec_integral,(ao_num_align,ao_num)]/;">`_
|
||||
array of the mono electronic hamiltonian on the AOs basis
|
||||
: sum of the kinetic and nuclear electronic potential
|
||||
|
||||
`ao_overlap <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/ao_mono_ints.irp.f#L1>`_
|
||||
`ao_overlap <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/ao_mono_ints.irp.f#L/BEGIN_PROVIDER [ double precision, ao_overlap,(ao_num_align,ao_num) ]/;">`_
|
||||
Overlap between atomic basis functions:
|
||||
:math:`\int \chi_i(r) \chi_j(r) dr)`
|
||||
|
||||
`ao_overlap_abs <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/ao_mono_ints.irp.f#L65>`_
|
||||
`ao_overlap_abs <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/ao_mono_ints.irp.f#L/BEGIN_PROVIDER [ double precision, ao_overlap_abs,(ao_num_align,ao_num) ]/;">`_
|
||||
Overlap between absolute value of atomic basis functions:
|
||||
:math:`\int |\chi_i(r)| |\chi_j(r)| dr)`
|
||||
|
||||
`ao_overlap_x <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/ao_mono_ints.irp.f#L2>`_
|
||||
`ao_overlap_x <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/ao_mono_ints.irp.f#L/&BEGIN_PROVIDER [ double precision, ao_overlap_x,(ao_num_align,ao_num) ]/;">`_
|
||||
Overlap between atomic basis functions:
|
||||
:math:`\int \chi_i(r) \chi_j(r) dr)`
|
||||
|
||||
`ao_overlap_y <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/ao_mono_ints.irp.f#L3>`_
|
||||
`ao_overlap_y <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/ao_mono_ints.irp.f#L/&BEGIN_PROVIDER [ double precision, ao_overlap_y,(ao_num_align,ao_num) ]/;">`_
|
||||
Overlap between atomic basis functions:
|
||||
:math:`\int \chi_i(r) \chi_j(r) dr)`
|
||||
|
||||
`ao_overlap_z <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/ao_mono_ints.irp.f#L4>`_
|
||||
`ao_overlap_z <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/ao_mono_ints.irp.f#L/&BEGIN_PROVIDER [ double precision, ao_overlap_z,(ao_num_align,ao_num) ]/;">`_
|
||||
Overlap between atomic basis functions:
|
||||
:math:`\int \chi_i(r) \chi_j(r) dr)`
|
||||
|
||||
`check_ortho <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/check_orthonormality.irp.f#L1>`_
|
||||
None
|
||||
`do_print <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/check_orthonormality.irp.f#L11>`_
|
||||
None
|
||||
`n_pt_max_i_x <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/dimensions.irp.f#L2>`_
|
||||
None
|
||||
`n_pt_max_integrals <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/dimensions.irp.f#L1>`_
|
||||
None
|
||||
`ao_deriv2_x <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/kin_ao_ints.irp.f#L1>`_
|
||||
`check_ortho <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/check_orthonormality.irp.f#L/subroutine check_ortho/;">`_
|
||||
Undocumented
|
||||
|
||||
`do_print <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/check_orthonormality.irp.f#L/subroutine do_print/;">`_
|
||||
Undocumented
|
||||
|
||||
`n_pt_max_i_x <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/dimensions.irp.f#L/&BEGIN_PROVIDER [ integer, n_pt_max_i_x]/;">`_
|
||||
Undocumented
|
||||
|
||||
`n_pt_max_integrals <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/dimensions.irp.f#L/BEGIN_PROVIDER [ integer, n_pt_max_integrals ]/;">`_
|
||||
Undocumented
|
||||
|
||||
`ao_deriv2_x <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/kin_ao_ints.irp.f#L/BEGIN_PROVIDER [ double precision, ao_deriv2_x,(ao_num_align,ao_num) ]/;">`_
|
||||
second derivatives matrix elements in the ao basis
|
||||
.. math::
|
||||
.br
|
||||
{\tt ao_deriv2_x} = \langle \chi_i(x,y,z) \frac{\partial^2}{\partial x^2} |\chi_j (x,y,z) \rangle
|
||||
|
||||
`ao_deriv2_y <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/kin_ao_ints.irp.f#L2>`_
|
||||
`ao_deriv2_y <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/kin_ao_ints.irp.f#L/&BEGIN_PROVIDER [ double precision, ao_deriv2_y,(ao_num_align,ao_num) ]/;">`_
|
||||
second derivatives matrix elements in the ao basis
|
||||
.. math::
|
||||
.br
|
||||
{\tt ao_deriv2_x} = \langle \chi_i(x,y,z) \frac{\partial^2}{\partial x^2} |\chi_j (x,y,z) \rangle
|
||||
|
||||
`ao_deriv2_z <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/kin_ao_ints.irp.f#L3>`_
|
||||
`ao_deriv2_z <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/kin_ao_ints.irp.f#L/&BEGIN_PROVIDER [ double precision, ao_deriv2_z,(ao_num_align,ao_num) ]/;">`_
|
||||
second derivatives matrix elements in the ao basis
|
||||
.. math::
|
||||
.br
|
||||
{\tt ao_deriv2_x} = \langle \chi_i(x,y,z) \frac{\partial^2}{\partial x^2} |\chi_j (x,y,z) \rangle
|
||||
|
||||
`ao_kinetic_integral <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/kin_ao_ints.irp.f#L125>`_
|
||||
`ao_kinetic_integral <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/kin_ao_ints.irp.f#L/BEGIN_PROVIDER [double precision, ao_kinetic_integral, (ao_num_align,ao_num)]/;">`_
|
||||
array of the priminitve basis kinetic integrals
|
||||
\langle \chi_i |\hat{T}| \chi_j \rangle
|
||||
|
||||
`mo_kinetic_integral <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/kin_mo_ints.irp.f#L1>`_
|
||||
None
|
||||
`mo_mono_elec_integral <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/mo_mono_ints.irp.f#L35>`_
|
||||
`mo_kinetic_integral <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/kin_mo_ints.irp.f#L/BEGIN_PROVIDER [double precision, mo_kinetic_integral, (mo_tot_num_align,mo_tot_num)]/;">`_
|
||||
Undocumented
|
||||
|
||||
`mo_mono_elec_integral <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/mo_mono_ints.irp.f#L/BEGIN_PROVIDER [ double precision, mo_mono_elec_integral,(mo_tot_num_align,mo_tot_num)]/;">`_
|
||||
array of the mono electronic hamiltonian on the MOs basis
|
||||
: sum of the kinetic and nuclear electronic potential
|
||||
|
||||
`mo_overlap <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/mo_mono_ints.irp.f#L1>`_
|
||||
None
|
||||
`orthonormalize_mos <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/orthonormalize.irp.f#L1>`_
|
||||
None
|
||||
`ao_nucl_elec_integral <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L1>`_
|
||||
`mo_overlap <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/mo_mono_ints.irp.f#L/BEGIN_PROVIDER [ double precision, mo_overlap,(mo_tot_num_align,mo_tot_num)]/;">`_
|
||||
Undocumented
|
||||
|
||||
`orthonormalize_mos <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/orthonormalize.irp.f#L/subroutine orthonormalize_mos/;">`_
|
||||
Undocumented
|
||||
|
||||
`ao_nucl_elec_integral <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L/BEGIN_PROVIDER [ double precision, ao_nucl_elec_integral, (ao_num_align,ao_num)]/;">`_
|
||||
interaction nuclear electron
|
||||
|
||||
`give_polynom_mult_center_mono_elec <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L157>`_
|
||||
None
|
||||
`i_x1_pol_mult_mono_elec <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L285>`_
|
||||
None
|
||||
`i_x2_pol_mult_mono_elec <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L357>`_
|
||||
None
|
||||
`int_gaus_pol <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L428>`_
|
||||
None
|
||||
`nai_pol_mult <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L82>`_
|
||||
None
|
||||
`v_e_n <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L409>`_
|
||||
None
|
||||
`v_phi <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L473>`_
|
||||
None
|
||||
`v_r <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L457>`_
|
||||
None
|
||||
`v_theta <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L486>`_
|
||||
None
|
||||
`wallis <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L502>`_
|
||||
None
|
||||
`mo_nucl_elec_integral <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_mo_ints.irp.f#L1>`_
|
||||
None
|
||||
`save_ortho_mos <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/save_ortho_mos.irp.f#L1>`_
|
||||
None
|
||||
`give_polynom_mult_center_mono_elec <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L/subroutine give_polynom_mult_center_mono_elec(A_center,B_center,alpha,beta,power_A,power_B,C_center,n_pt_in,d,n_pt_out)/;">`_
|
||||
Undocumented
|
||||
|
||||
`i_x1_pol_mult_mono_elec <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L/recursive subroutine I_x1_pol_mult_mono_elec(a,c,R1x,R1xp,R2x,d,nd,n_pt_in)/;">`_
|
||||
Undocumented
|
||||
|
||||
`i_x2_pol_mult_mono_elec <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L/recursive subroutine I_x2_pol_mult_mono_elec(c,R1x,R1xp,R2x,d,nd,dim)/;">`_
|
||||
Undocumented
|
||||
|
||||
`int_gaus_pol <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L/double precision function int_gaus_pol(alpha,n)/;">`_
|
||||
Undocumented
|
||||
|
||||
`nai_pol_mult <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L/double precision function NAI_pol_mult(A_center,B_center,power_A,power_B,alpha,beta,C_center,n_pt_in)/;">`_
|
||||
Undocumented
|
||||
|
||||
`v_e_n <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L/double precision function V_e_n(a_x,a_y,a_z,b_x,b_y,b_z,alpha,beta)/;">`_
|
||||
Undocumented
|
||||
|
||||
`v_phi <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L/double precision function V_phi(n,m)/;">`_
|
||||
Undocumented
|
||||
|
||||
`v_r <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L/double precision function V_r(n,alpha)/;">`_
|
||||
Undocumented
|
||||
|
||||
`v_theta <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L/double precision function V_theta(n,m)/;">`_
|
||||
Undocumented
|
||||
|
||||
`wallis <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L/double precision function Wallis(n)/;">`_
|
||||
Undocumented
|
||||
|
||||
`mo_nucl_elec_integral <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_mo_ints.irp.f#L/BEGIN_PROVIDER [double precision, mo_nucl_elec_integral, (mo_tot_num_align,mo_tot_num)]/;">`_
|
||||
Undocumented
|
||||
|
||||
`save_ortho_mos <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/save_ortho_mos.irp.f#L/subroutine save_ortho_mos/;">`_
|
||||
Undocumented
|
||||
|
||||
|
||||
|
||||
|
416
src/MonoInts/spread_dipole_ao.irp.f
Normal file
416
src/MonoInts/spread_dipole_ao.irp.f
Normal file
@ -0,0 +1,416 @@
|
||||
BEGIN_PROVIDER [ double precision, ao_spread_x, (ao_num_align,ao_num)]
|
||||
&BEGIN_PROVIDER [ double precision, ao_spread_y, (ao_num_align,ao_num)]
|
||||
&BEGIN_PROVIDER [ double precision, ao_spread_z, (ao_num_align,ao_num)]
|
||||
BEGIN_DOC
|
||||
! array of the integrals of AO_i * x^2 AO_j
|
||||
! array of the integrals of AO_i * y^2 AO_j
|
||||
! array of the integrals of AO_i * z^2 AO_j
|
||||
END_DOC
|
||||
implicit none
|
||||
integer :: i,j,n,l
|
||||
double precision :: f, tmp
|
||||
integer :: dim1
|
||||
double precision :: overlap, overlap_x, overlap_y, overlap_z
|
||||
double precision :: alpha, beta
|
||||
double precision :: A_center(3), B_center(3)
|
||||
integer :: power_A(3), power_B(3)
|
||||
double precision :: lower_exp_val, dx, c,accu_x,accu_y,accu_z
|
||||
dim1=500
|
||||
lower_exp_val = 40.d0
|
||||
ao_spread_x= 0.d0
|
||||
ao_spread_y= 0.d0
|
||||
ao_spread_z= 0.d0
|
||||
!$OMP PARALLEL DO SCHEDULE(GUIDED) &
|
||||
!$OMP DEFAULT(NONE) &
|
||||
!$OMP PRIVATE(A_center,B_center,power_A,power_B,&
|
||||
!$OMP overlap_x,overlap_y, overlap_z, overlap, &
|
||||
!$OMP alpha, beta,i,j,dx,tmp,c,accu_x,accu_y,accu_z) &
|
||||
!$OMP SHARED(nucl_coord,ao_power,ao_prim_num, &
|
||||
!$OMP ao_spread_x,ao_spread_y,ao_spread_z,ao_num,ao_coef_transp,ao_nucl, &
|
||||
!$OMP ao_expo_transp,dim1,lower_exp_val)
|
||||
do j=1,ao_num
|
||||
A_center(1) = nucl_coord( ao_nucl(j), 1 )
|
||||
A_center(2) = nucl_coord( ao_nucl(j), 2 )
|
||||
A_center(3) = nucl_coord( ao_nucl(j), 3 )
|
||||
power_A(1) = ao_power( j, 1 )
|
||||
power_A(2) = ao_power( j, 2 )
|
||||
power_A(3) = ao_power( j, 3 )
|
||||
!DEC$ VECTOR ALIGNED
|
||||
!DEC$ VECTOR ALWAYS
|
||||
do i= 1,ao_num
|
||||
B_center(1) = nucl_coord( ao_nucl(i), 1 )
|
||||
B_center(2) = nucl_coord( ao_nucl(i), 2 )
|
||||
B_center(3) = nucl_coord( ao_nucl(i), 3 )
|
||||
power_B(1) = ao_power( i, 1 )
|
||||
power_B(2) = ao_power( i, 2 )
|
||||
power_B(3) = ao_power( i, 3 )
|
||||
accu_x = 0.d0
|
||||
accu_y = 0.d0
|
||||
accu_z = 0.d0
|
||||
do n = 1,ao_prim_num(j)
|
||||
alpha = ao_expo_transp(n,j)
|
||||
!DEC$ VECTOR ALIGNED
|
||||
do l = 1, ao_prim_num(i)
|
||||
c = ao_coef_transp(n,j)*ao_coef_transp(l,i)
|
||||
beta = ao_expo_transp(l,i)
|
||||
call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim1)
|
||||
call overlap_bourrin_spread(A_center(1),B_center(1),alpha,beta,power_A(1),power_B(1),tmp,lower_exp_val,dx,dim1)
|
||||
accu_x += c*(tmp*overlap_y*overlap_z)
|
||||
call overlap_bourrin_spread(A_center(2),B_center(2),alpha,beta,power_A(2),power_B(2),tmp,lower_exp_val,dx,dim1)
|
||||
accu_y += c*(tmp*overlap_x*overlap_z)
|
||||
call overlap_bourrin_spread(A_center(3),B_center(3),alpha,beta,power_A(3),power_B(3),tmp,lower_exp_val,dx,dim1)
|
||||
accu_z += c*(tmp*overlap_y*overlap_x)
|
||||
enddo
|
||||
enddo
|
||||
ao_spread_x(i,j) = accu_x
|
||||
ao_spread_y(i,j) = accu_y
|
||||
ao_spread_z(i,j) = accu_z
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END PARALLEL DO
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ double precision, ao_dipole_x, (ao_num_align,ao_num)]
|
||||
&BEGIN_PROVIDER [ double precision, ao_dipole_y, (ao_num_align,ao_num)]
|
||||
&BEGIN_PROVIDER [ double precision, ao_dipole_z, (ao_num_align,ao_num)]
|
||||
BEGIN_DOC
|
||||
! array of the integrals of AO_i * x AO_j
|
||||
! array of the integrals of AO_i * y AO_j
|
||||
! array of the integrals of AO_i * z AO_j
|
||||
END_DOC
|
||||
implicit none
|
||||
integer :: i,j,n,l
|
||||
double precision :: f, tmp
|
||||
integer :: dim1
|
||||
double precision :: overlap, overlap_x, overlap_y, overlap_z,accu_x,accu_y,accu_z
|
||||
double precision :: alpha, beta
|
||||
double precision :: A_center(3), B_center(3)
|
||||
integer :: power_A(3), power_B(3)
|
||||
double precision :: lower_exp_val, dx, c
|
||||
dim1=500
|
||||
lower_exp_val = 40.d0
|
||||
ao_dipole_x= 0.d0
|
||||
ao_dipole_y= 0.d0
|
||||
ao_dipole_z= 0.d0
|
||||
!$OMP PARALLEL DO SCHEDULE(GUIDED) &
|
||||
!$OMP DEFAULT(NONE) &
|
||||
!$OMP PRIVATE(A_center,B_center,power_A,power_B,&
|
||||
!$OMP overlap_x,overlap_y, overlap_z, overlap, &
|
||||
!$OMP alpha, beta,i,j,dx,tmp,c,accu_x,accu_y,accu_z) &
|
||||
!$OMP SHARED(nucl_coord,ao_power,ao_prim_num, &
|
||||
!$OMP ao_dipole_x,ao_dipole_y,ao_dipole_z,ao_num,ao_coef_transp,ao_nucl, &
|
||||
!$OMP ao_expo_transp,dim1,lower_exp_val)
|
||||
do j=1,ao_num
|
||||
A_center(1) = nucl_coord( ao_nucl(j), 1 )
|
||||
A_center(2) = nucl_coord( ao_nucl(j), 2 )
|
||||
A_center(3) = nucl_coord( ao_nucl(j), 3 )
|
||||
power_A(1) = ao_power( j, 1 )
|
||||
power_A(2) = ao_power( j, 2 )
|
||||
power_A(3) = ao_power( j, 3 )
|
||||
!DEC$ VECTOR ALIGNED
|
||||
!DEC$ VECTOR ALWAYS
|
||||
do i= 1,ao_num
|
||||
B_center(1) = nucl_coord( ao_nucl(i), 1 )
|
||||
B_center(2) = nucl_coord( ao_nucl(i), 2 )
|
||||
B_center(3) = nucl_coord( ao_nucl(i), 3 )
|
||||
power_B(1) = ao_power( i, 1 )
|
||||
power_B(2) = ao_power( i, 2 )
|
||||
power_B(3) = ao_power( i, 3 )
|
||||
accu_x = 0.d0
|
||||
accu_y = 0.d0
|
||||
accu_z = 0.d0
|
||||
do n = 1,ao_prim_num(j)
|
||||
alpha = ao_expo_transp(n,j)
|
||||
!DEC$ VECTOR ALIGNED
|
||||
do l = 1, ao_prim_num(i)
|
||||
beta = ao_expo_transp(l,i)
|
||||
c = ao_coef_transp(l,i)*ao_coef_transp(n,j)
|
||||
call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim1)
|
||||
|
||||
call overlap_bourrin_dipole(A_center(1),B_center(1),alpha,beta,power_A(1),power_B(1),tmp,lower_exp_val,dx,dim1)
|
||||
accu_x = accu_x + c*(tmp*overlap_y*overlap_z)
|
||||
call overlap_bourrin_dipole(A_center(2),B_center(2),alpha,beta,power_A(2),power_B(2),tmp,lower_exp_val,dx,dim1)
|
||||
accu_y = accu_y + c*(tmp*overlap_x*overlap_z)
|
||||
call overlap_bourrin_dipole(A_center(3),B_center(3),alpha,beta,power_A(3),power_B(3),tmp,lower_exp_val,dx,dim1)
|
||||
accu_z = accu_z + c*(tmp*overlap_y*overlap_x)
|
||||
enddo
|
||||
enddo
|
||||
ao_dipole_x(i,j) = accu_x
|
||||
ao_dipole_y(i,j) = accu_y
|
||||
ao_dipole_z(i,j) = accu_z
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END PARALLEL DO
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ double precision, ao_deriv_1_x, (ao_num_align,ao_num)]
|
||||
&BEGIN_PROVIDER [ double precision, ao_deriv_1_y, (ao_num_align,ao_num)]
|
||||
&BEGIN_PROVIDER [ double precision, ao_deriv_1_z, (ao_num_align,ao_num)]
|
||||
BEGIN_DOC
|
||||
! array of the integrals of AO_i * d/dx AO_j
|
||||
! array of the integrals of AO_i * d/dy AO_j
|
||||
! array of the integrals of AO_i * d/dz AO_j
|
||||
END_DOC
|
||||
implicit none
|
||||
integer :: i,j,n,l
|
||||
double precision :: f, tmp
|
||||
integer :: dim1
|
||||
double precision :: overlap, overlap_x, overlap_y, overlap_z
|
||||
double precision :: alpha, beta
|
||||
double precision :: A_center(3), B_center(3)
|
||||
integer :: power_A(3), power_B(3)
|
||||
double precision :: lower_exp_val, dx, c,accu_x,accu_y,accu_z
|
||||
integer :: i_component
|
||||
dim1=500
|
||||
lower_exp_val = 40.d0
|
||||
ao_deriv_1_x= 0.d0
|
||||
ao_deriv_1_y= 0.d0
|
||||
ao_deriv_1_z= 0.d0
|
||||
!$OMP PARALLEL DO SCHEDULE(GUIDED) &
|
||||
!$OMP DEFAULT(NONE) &
|
||||
!$OMP PRIVATE(A_center,B_center,power_A,power_B,&
|
||||
!$OMP overlap_x,overlap_y, overlap_z, overlap, &
|
||||
!$OMP alpha, beta,i,j,dx,tmp,c,i_component,accu_x,accu_y,accu_z) &
|
||||
!$OMP SHARED(nucl_coord,ao_power,ao_prim_num, &
|
||||
!$OMP ao_deriv_1_x,ao_deriv_1_y,ao_deriv_1_z,ao_num,ao_coef_transp,ao_nucl, &
|
||||
!$OMP ao_expo_transp,dim1,lower_exp_val)
|
||||
do j=1,ao_num
|
||||
A_center(1) = nucl_coord( ao_nucl(j), 1 )
|
||||
A_center(2) = nucl_coord( ao_nucl(j), 2 )
|
||||
A_center(3) = nucl_coord( ao_nucl(j), 3 )
|
||||
power_A(1) = ao_power( j, 1 )
|
||||
power_A(2) = ao_power( j, 2 )
|
||||
power_A(3) = ao_power( j, 3 )
|
||||
!DEC$ VECTOR ALIGNED
|
||||
!DEC$ VECTOR ALWAYS
|
||||
do i= 1,ao_num
|
||||
B_center(1) = nucl_coord( ao_nucl(i), 1 )
|
||||
B_center(2) = nucl_coord( ao_nucl(i), 2 )
|
||||
B_center(3) = nucl_coord( ao_nucl(i), 3 )
|
||||
power_B(1) = ao_power( i, 1 )
|
||||
power_B(2) = ao_power( i, 2 )
|
||||
power_B(3) = ao_power( i, 3 )
|
||||
accu_x = 0.d0
|
||||
accu_y = 0.d0
|
||||
accu_z = 0.d0
|
||||
do n = 1,ao_prim_num(j)
|
||||
alpha = ao_expo_transp(n,j)
|
||||
!DEC$ VECTOR ALIGNED
|
||||
do l = 1, ao_prim_num(i)
|
||||
beta = ao_expo_transp(l,i)
|
||||
call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim1)
|
||||
c = ao_coef_transp(l,i) * ao_coef_transp(n,j)
|
||||
i_component = 1
|
||||
call overlap_bourrin_deriv_x(i_component,A_center,B_center,alpha,beta,power_A,power_B,dx,lower_exp_val,tmp,dim1)
|
||||
accu_x += c*(tmp*overlap_y*overlap_z)
|
||||
i_component = 2
|
||||
call overlap_bourrin_deriv_x(i_component,A_center,B_center,alpha,beta,power_A,power_B,dx,lower_exp_val,tmp,dim1)
|
||||
accu_y += c*(tmp*overlap_x*overlap_z)
|
||||
i_component = 3
|
||||
call overlap_bourrin_deriv_x(i_component,A_center,B_center,alpha,beta,power_A,power_B,dx,lower_exp_val,tmp,dim1)
|
||||
accu_z += c*(tmp*overlap_y*overlap_x)
|
||||
enddo
|
||||
enddo
|
||||
ao_deriv_1_x(i,j) = accu_x
|
||||
ao_deriv_1_y(i,j) = accu_y
|
||||
ao_deriv_1_z(i,j) = accu_z
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END PARALLEL DO
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
|
||||
subroutine overlap_bourrin_x_abs(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,lower_exp_val,dx,nx)
|
||||
implicit none
|
||||
! compute the following integral :
|
||||
! int [-infty ; +infty] of [(x-A_center)^(power_A) * (x-B_center)^power_B * exp(-alpha(x-A_center)^2) * exp(-beta(x-B_center)^2) ]
|
||||
integer :: i,j,k,l
|
||||
integer,intent(in) :: power_A,power_B
|
||||
double precision, intent(in) :: lower_exp_val
|
||||
double precision,intent(in) :: A_center, B_center,alpha,beta
|
||||
double precision, intent(out) :: overlap_x,dx
|
||||
integer, intent(in) :: nx
|
||||
double precision :: x_min,x_max,domain,x,power,factor,dist,p,p_inv,rho
|
||||
double precision :: P_center,pouet_timy
|
||||
if(power_A.lt.0.or.power_B.lt.0)then
|
||||
overlap_x = 0.d0
|
||||
dx = 0.d0
|
||||
return
|
||||
endif
|
||||
p = alpha + beta
|
||||
p_inv= 1.d0/p
|
||||
rho = alpha * beta * p_inv
|
||||
dist = (A_center - B_center)*(A_center - B_center)
|
||||
P_center = (alpha * A_center + beta * B_center) * p_inv
|
||||
factor = dexp(-rho * dist)
|
||||
pouet_timy = dsqrt(lower_exp_val/p)
|
||||
x_min = P_center - pouet_timy
|
||||
x_max = P_center + pouet_timy
|
||||
!print*,'xmin = ',x_min
|
||||
!print*,'xmax = ',x_max
|
||||
domain = x_max-x_min
|
||||
dx = domain/dble(nx)
|
||||
overlap_x = 0.d0
|
||||
x = x_min
|
||||
do i = 1, nx
|
||||
x += dx
|
||||
overlap_x += (power(power_A,x-A_center) * power(power_B,x-B_center)) * dexp(-p * (x-P_center)*(x-P_center))
|
||||
enddo
|
||||
overlap_x *= factor * dx
|
||||
end
|
||||
|
||||
subroutine overlap_bourrin_spread(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,lower_exp_val,dx,nx)
|
||||
! compute the following integral :
|
||||
! int [-infty ; +infty] of [(x-A_center)^(power_A) * (x-B_center)^power_B * exp(-alpha(x-A_center)^2) * exp(-beta(x-B_center)^2) * x ]
|
||||
! needed for the dipole and those things
|
||||
implicit none
|
||||
integer :: i,j,k,l
|
||||
integer,intent(in) :: power_A,power_B
|
||||
double precision, intent(in) :: lower_exp_val
|
||||
double precision,intent(in) :: A_center, B_center,alpha,beta
|
||||
double precision, intent(out) :: overlap_x,dx
|
||||
integer, intent(in) :: nx
|
||||
double precision :: x_min,x_max,domain,x,power,factor,dist,p,p_inv,rho
|
||||
double precision :: P_center,pouet_timy
|
||||
if(power_A.lt.0.or.power_B.lt.0)then
|
||||
overlap_x = 0.d0
|
||||
dx = 0.d0
|
||||
return
|
||||
endif
|
||||
p = alpha + beta
|
||||
p_inv= 1.d0/p
|
||||
rho = alpha * beta * p_inv
|
||||
dist = (A_center - B_center)*(A_center - B_center)
|
||||
P_center = (alpha * A_center + beta * B_center) * p_inv
|
||||
factor = dexp(-rho * dist)
|
||||
if(factor.lt.0.000001d0)then
|
||||
! print*,'factor = ',factor
|
||||
dx = 0.d0
|
||||
overlap_x = 0.d0
|
||||
return
|
||||
endif
|
||||
pouet_timy = dsqrt(lower_exp_val/p)
|
||||
x_min = P_center - pouet_timy
|
||||
x_max = P_center + pouet_timy
|
||||
domain = x_max-x_min
|
||||
dx = domain/dble(nx)
|
||||
overlap_x = 0.d0
|
||||
x = x_min
|
||||
do i = 1, nx
|
||||
x += dx
|
||||
overlap_x += power(power_A,x-A_center) * power(power_B,x-B_center) * dexp(-p * (x-P_center)*(x-P_center)) * x * x
|
||||
enddo
|
||||
overlap_x *= factor * dx
|
||||
|
||||
end
|
||||
|
||||
double precision function power(n,x)
|
||||
implicit none
|
||||
integer :: i,n
|
||||
double precision :: x,accu
|
||||
power = x**n
|
||||
return
|
||||
end
|
||||
|
||||
subroutine overlap_bourrin_dipole(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,lower_exp_val,dx,nx)
|
||||
! compute the following integral :
|
||||
! int [-infty ; +infty] of [(x-A_center)^(power_A) * (x-B_center)^power_B * exp(-alpha(x-A_center)^2) * exp(-beta(x-B_center)^2) * x ]
|
||||
! needed for the dipole and those things
|
||||
implicit none
|
||||
integer :: i,j,k,l
|
||||
integer,intent(in) :: power_A,power_B
|
||||
double precision, intent(in) :: lower_exp_val
|
||||
double precision,intent(in) :: A_center, B_center,alpha,beta
|
||||
double precision, intent(out) :: overlap_x,dx
|
||||
integer, intent(in) :: nx
|
||||
double precision :: x_min,x_max,domain,x,power,factor,dist,p,p_inv,rho
|
||||
double precision :: P_center
|
||||
if(power_A.lt.0.or.power_B.lt.0)then
|
||||
overlap_x = 0.d0
|
||||
dx = 0.d0
|
||||
return
|
||||
endif
|
||||
p = alpha + beta
|
||||
p_inv= 1.d0/p
|
||||
rho = alpha * beta * p_inv
|
||||
dist = (A_center - B_center)*(A_center - B_center)
|
||||
P_center = (alpha * A_center + beta * B_center) * p_inv
|
||||
factor = dexp(-rho * dist)
|
||||
double precision :: pouet_timy
|
||||
|
||||
pouet_timy = dsqrt(lower_exp_val/p)
|
||||
x_min = P_center - pouet_timy
|
||||
x_max = P_center + pouet_timy
|
||||
domain = x_max-x_min
|
||||
dx = domain/dble(nx)
|
||||
overlap_x = 0.d0
|
||||
x = x_min
|
||||
do i = 1, nx
|
||||
x += dx
|
||||
overlap_x += power(power_A,x-A_center) * power(power_B,x-B_center) * dexp(-p * (x-P_center)*(x-P_center)) * x
|
||||
enddo
|
||||
overlap_x *= factor * dx
|
||||
|
||||
end
|
||||
|
||||
subroutine overlap_bourrin_deriv_x(i_component,A_center,B_center,alpha,beta,power_A,power_B,dx,lower_exp_val,overlap_x,nx)
|
||||
implicit none
|
||||
integer :: i,j,k,l
|
||||
integer,intent(in) :: power_A(3),power_B(3),i_component
|
||||
double precision,intent(in) :: A_center(3), B_center(3),alpha,beta,lower_exp_val
|
||||
double precision, intent(out) :: overlap_x,dx
|
||||
integer, intent(in) :: nx
|
||||
double precision :: overlap_first, overlap_second
|
||||
! computes : <phi_i|d/dx|phi_j> = (a_x_i <phi_i_x|phi_j_x(a_x_j-1)> - 2 alpha <phi_i_x|phi_j_w(a_x_j+1)>)
|
||||
|
||||
call overlap_bourrin_x(A_center(i_component),B_center(i_component),alpha,beta,power_A(i_component)-1,power_B(i_component),overlap_first,lower_exp_val,dx,nx)
|
||||
call overlap_bourrin_x(A_center(i_component),B_center(i_component),alpha,beta,power_A(i_component)+1,power_B(i_component),overlap_second,lower_exp_val,dx,nx)
|
||||
overlap_x = (power_A(i_component) * overlap_first - 2.d0 * alpha * overlap_second)
|
||||
end
|
||||
|
||||
subroutine overlap_bourrin_x(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,lower_exp_val,dx,nx)
|
||||
implicit none
|
||||
! compute the following integral :
|
||||
! int [-infty ; +infty] of [(x-A_center)^(power_A) * (x-B_center)^power_B * exp(-alpha(x-A_center)^2) * exp(-beta(x-B_center)^2) ]
|
||||
integer :: i,j,k,l
|
||||
integer,intent(in) :: power_A,power_B
|
||||
double precision, intent(in) :: lower_exp_val
|
||||
double precision,intent(in) :: A_center, B_center,alpha,beta
|
||||
double precision, intent(out) :: overlap_x,dx
|
||||
integer, intent(in) :: nx
|
||||
double precision :: x_min,x_max,domain,x,power,factor,dist,p,p_inv,rho
|
||||
double precision :: P_center,pouet_timy
|
||||
if(power_A.lt.0.or.power_B.lt.0)then
|
||||
overlap_x = 0.d0
|
||||
dx = 0.d0
|
||||
return
|
||||
endif
|
||||
p = alpha + beta
|
||||
p_inv= 1.d0/p
|
||||
rho = alpha * beta * p_inv
|
||||
dist = (A_center - B_center)*(A_center - B_center)
|
||||
P_center = (alpha * A_center + beta * B_center) * p_inv
|
||||
factor = dexp(-rho * dist)
|
||||
if(factor.lt.0.000001d0)then
|
||||
dx = 0.d0
|
||||
overlap_x = 0.d0
|
||||
return
|
||||
endif
|
||||
|
||||
pouet_timy = dsqrt(lower_exp_val/p)
|
||||
x_min = P_center - pouet_timy
|
||||
x_max = P_center + pouet_timy
|
||||
domain = x_max-x_min
|
||||
dx = domain/dble(nx)
|
||||
overlap_x = 0.d0
|
||||
x = x_min
|
||||
do i = 1, nx
|
||||
x += dx
|
||||
overlap_x += power(power_A,x-A_center) * power(power_B,x-B_center) * dexp(-p * (x-P_center)*(x-P_center))
|
||||
enddo
|
||||
overlap_x *= factor * dx
|
||||
end
|
||||
|
101
src/MonoInts/spread_dipole_mo.irp.f
Normal file
101
src/MonoInts/spread_dipole_mo.irp.f
Normal file
@ -0,0 +1,101 @@
|
||||
BEGIN_PROVIDER [double precision, mo_dipole_x , (mo_tot_num_align,mo_tot_num)]
|
||||
&BEGIN_PROVIDER [double precision, mo_dipole_y , (mo_tot_num_align,mo_tot_num)]
|
||||
&BEGIN_PROVIDER [double precision, mo_dipole_z , (mo_tot_num_align,mo_tot_num)]
|
||||
BEGIN_DOC
|
||||
! array of the integrals of MO_i * x MO_j
|
||||
! array of the integrals of MO_i * y MO_j
|
||||
! array of the integrals of MO_i * z MO_j
|
||||
END_DOC
|
||||
implicit none
|
||||
integer :: i1,j1,i,j
|
||||
double precision :: c_i1,c_j1
|
||||
|
||||
mo_dipole_x = 0.d0
|
||||
mo_dipole_y = 0.d0
|
||||
mo_dipole_z = 0.d0
|
||||
!$OMP PARALLEL DO DEFAULT(none) &
|
||||
!$OMP PRIVATE(i,j,i1,j1,c_j1,c_i1) &
|
||||
!$OMP SHARED(mo_tot_num,ao_num,mo_coef, &
|
||||
!$OMP mo_dipole_x,mo_dipole_y,mo_dipole_z,ao_dipole_x,ao_dipole_y,ao_dipole_z)
|
||||
do i = 1, mo_tot_num
|
||||
do j = 1, mo_tot_num
|
||||
do i1 = 1,ao_num
|
||||
c_i1 = mo_coef(i1,i)
|
||||
do j1 = 1,ao_num
|
||||
c_j1 = c_i1*mo_coef(j1,j)
|
||||
mo_dipole_x(j,i) = mo_dipole_x(j,i) + c_j1 * ao_dipole_x(j1,i1)
|
||||
mo_dipole_y(j,i) = mo_dipole_y(j,i) + c_j1 * ao_dipole_y(j1,i1)
|
||||
mo_dipole_z(j,i) = mo_dipole_z(j,i) + c_j1 * ao_dipole_z(j1,i1)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END PARALLEL DO
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [double precision, mo_spread_x , (mo_tot_num_align,mo_tot_num)]
|
||||
&BEGIN_PROVIDER [double precision, mo_spread_y , (mo_tot_num_align,mo_tot_num)]
|
||||
&BEGIN_PROVIDER [double precision, mo_spread_z , (mo_tot_num_align,mo_tot_num)]
|
||||
BEGIN_DOC
|
||||
! array of the integrals of MO_i * x^2 MO_j
|
||||
! array of the integrals of MO_i * y^2 MO_j
|
||||
! array of the integrals of MO_i * z^2 MO_j
|
||||
END_DOC
|
||||
implicit none
|
||||
integer :: i1,j1,i,j
|
||||
double precision :: c_i1,c_j1
|
||||
|
||||
mo_nucl_elec_integral = 0.d0
|
||||
!$OMP PARALLEL DO DEFAULT(none) &
|
||||
!$OMP PRIVATE(i,j,i1,j1,c_j1,c_i1) &
|
||||
!$OMP SHARED(mo_tot_num,ao_num,mo_coef, &
|
||||
!$OMP mo_spread_x,mo_spread_y,mo_spread_z,ao_spread_x,ao_spread_y,ao_spread_z)
|
||||
do i = 1, mo_tot_num
|
||||
do j = 1, mo_tot_num
|
||||
do i1 = 1,ao_num
|
||||
c_i1 = mo_coef(i1,i)
|
||||
do j1 = 1,ao_num
|
||||
c_j1 = c_i1*mo_coef(j1,j)
|
||||
mo_spread_x(j,i) = mo_spread_x(j,i) + c_j1 * ao_spread_x(j1,i1)
|
||||
mo_spread_y(j,i) = mo_spread_y(j,i) + c_j1 * ao_spread_y(j1,i1)
|
||||
mo_spread_z(j,i) = mo_spread_z(j,i) + c_j1 * ao_spread_z(j1,i1)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END PARALLEL DO
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [double precision, mo_deriv_1_x , (mo_tot_num_align,mo_tot_num)]
|
||||
&BEGIN_PROVIDER [double precision, mo_deriv_1_y , (mo_tot_num_align,mo_tot_num)]
|
||||
&BEGIN_PROVIDER [double precision, mo_deriv_1_z , (mo_tot_num_align,mo_tot_num)]
|
||||
BEGIN_DOC
|
||||
! array of the integrals of MO_i * d/dx MO_j
|
||||
! array of the integrals of MO_i * d/dy MO_j
|
||||
! array of the integrals of MO_i * d/dz MO_j
|
||||
END_DOC
|
||||
implicit none
|
||||
integer :: i1,j1,i,j
|
||||
double precision :: c_i1,c_j1
|
||||
|
||||
mo_nucl_elec_integral = 0.d0
|
||||
!$OMP PARALLEL DO DEFAULT(none) &
|
||||
!$OMP PRIVATE(i,j,i1,j1,c_j1,c_i1) &
|
||||
!$OMP SHARED(mo_tot_num,ao_num,mo_coef, &
|
||||
!$OMP mo_deriv_1_x,mo_deriv_1_y,mo_deriv_1_z,ao_spread_x,ao_spread_y,ao_spread_z)
|
||||
do i = 1, mo_tot_num
|
||||
do j = 1, mo_tot_num
|
||||
do i1 = 1,ao_num
|
||||
c_i1 = mo_coef(i1,i)
|
||||
do j1 = 1,ao_num
|
||||
c_j1 = c_i1*mo_coef(j1,j)
|
||||
mo_deriv_1_x(j,i) = mo_deriv_1_x(j,i) + c_j1 * ao_spread_x(j1,i1)
|
||||
mo_deriv_1_y(j,i) = mo_deriv_1_y(j,i) + c_j1 * ao_spread_y(j1,i1)
|
||||
mo_deriv_1_z(j,i) = mo_deriv_1_z(j,i) + c_j1 * ao_spread_z(j1,i1)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END PARALLEL DO
|
||||
END_PROVIDER
|
||||
|
@ -129,4 +129,13 @@ Needed Modules
|
||||
* `Dets <http://github.com/LCPQ/quantum_package/tree/master/src/Dets>`_
|
||||
* `DensityMatrix <http://github.com/LCPQ/quantum_package/tree/master/src/DensityMatrix>`_
|
||||
* `CISD <http://github.com/LCPQ/quantum_package/tree/master/src/CISD>`_
|
||||
* `Perturbation <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation>`_
|
||||
|
||||
Documentation
|
||||
=============
|
||||
|
||||
.. Do not edit this section. It was auto-generated from the
|
||||
.. NEEDED_MODULES file.
|
||||
|
||||
|
||||
|
||||
|
@ -181,6 +181,13 @@ subroutine lapack_diag(eigvalues,eigvectors,H,nmax,n)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Diagonalize matrix H
|
||||
!
|
||||
! H is untouched between input and ouptut
|
||||
!
|
||||
! eigevalues(i) = ith lowest eigenvalue of the H matrix
|
||||
!
|
||||
! eigvectors(i,j) = <i|psi_j> where i is the basis function and psi_j is the j th eigenvector
|
||||
!
|
||||
END_DOC
|
||||
integer, intent(in) :: n,nmax
|
||||
double precision, intent(out) :: eigvectors(nmax,n)
|
||||
@ -189,7 +196,6 @@ subroutine lapack_diag(eigvalues,eigvectors,H,nmax,n)
|
||||
double precision,allocatable :: eigenvalues(:)
|
||||
double precision,allocatable :: work(:)
|
||||
double precision,allocatable :: A(:,:)
|
||||
!eigvectors(i,j) = <d_i|psi_j> where d_i is the basis function and psi_j is the j th eigenvector
|
||||
allocate(A(nmax,n),eigenvalues(nmax),work(4*nmax))
|
||||
integer :: LWORK, info, i,j,l,k
|
||||
A=H
|
||||
|
Loading…
Reference in New Issue
Block a user