LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-12-23 04:43:50 +01:00
Code Issues Releases Wiki Activity

CIS(DT) works

Browse Source
This commit is contained in:
Manu 2014-05-19 18:35:56 +02:00
parent f55948890a
commit 3b18250977
20 changed files with 3351 additions and 85 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
src/CISD/README.rst
View File

@ -31,7 +31,7 @@ Documentation
excitations (of the same symmetry).
`cisd <http://github.com/LCPQ/quantum_package/tree/master/src/CISD/cisd.irp.f#L/subroutine cisd/;">`_
None
=======
Undocumented

421
src/CISD/SC2.irp.f Normal file
View 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

4
src/CISD/cisd.irp.f
View File

@ -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
@ -21,6 +23,8 @@ program cisd
print *, '---'
do i = 1,1
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
View 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
View File

File diff suppressed because it is too large Load Diff

28
src/CIS_dressed/README.rst Normal file
View 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
View 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
View 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
View 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
View 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

34
src/DensityMatrix/README.rst
View File

@ -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>`_

16
src/DensityMatrix/det_num.irp.f
View File

@ -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

1
src/Dets/davidson.irp.f
View File

@ -275,3 +275,4 @@ subroutine davidson_diag(dets_in,u_in,energies,dim_in,sze,N_st,Nint)
)
end

40
src/Dets/s2.irp.f
View File

@ -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

42
src/Hartree_Fock/README.rst
View File

@ -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

119
src/MonoInts/README.rst
View File

@ -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
View 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
View 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

1
src/README.rst
View File

@ -129,6 +129,7 @@ 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
=============

8
src/Utils/LinearAlgebra.irp.f
View 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