LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2025-04-15 21:19:40 +02:00
Code Issues Releases Wiki Activity

Merge branch 'master' of github.com:LCPQ/quantum_package

Browse Source
This commit is contained in:
Anthony Scemama 2014-05-21 22:46:49 +02:00
commit 6b3c7c29b0
40 changed files with 4322 additions and 353 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
scripts/check_dependencies.sh
View File

@ -51,7 +51,7 @@ DEPS=$(unique_list $DEPS_LONG)
if [[ ! "$COMMAND_LINE" == "$DEPS" ]]
then
DEPS=$(check_dependencies.sh $DEPS)
DEPS=$(${QPACKAGE_ROOT}/scripts/check_dependencies.sh ${DEPS})
fi
echo $DEPS

95
scripts/generate_h_apply.py
View File

@ -10,14 +10,25 @@ subroutine
parameters
initialization
declarations
keys_work
keys_work_locked
keys_work_unlocked
finalization
""".split()
def new_dict(openmp=True):
s ={}
class H_apply(object):
def __init__(self,sub,openmp=True):
s = {}
for k in keywords:
s[k] = ""
s["subroutine"] = "H_apply_%s"%(sub)
self.openmp = openmp
if openmp:
s["subroutine"] += "_OpenMP"
self.selection_pt2 = None
self.perturbation = None
#s["omp_parallel"] = """!$OMP PARALLEL DEFAULT(NONE) &
s["omp_parallel"] = """!$OMP PARALLEL DEFAULT(SHARED) &
!$OMP PRIVATE(i,j,k,l,keys_out,hole,particle, &
@ -27,9 +38,9 @@ def new_dict(openmp=True):
!$OMP occ_hole_tmp,key_idx,i_b,j_b,key,N_elec_in_key_part_1,&
!$OMP N_elec_in_key_hole_1,N_elec_in_key_part_2, &
!$OMP N_elec_in_key_hole_2,ia_ja_pairs) &
!$OMP SHARED(key_in,N_int,elec_num_tab, &
!$OMP SHARED(key_in,N_int,elec_num_tab,mo_tot_num, &
!$OMP hole_1, particl_1, hole_2, particl_2, &
!$OMP lck,thresh,elec_alpha_num,E_ref)"""
!$OMP lck,thresh,elec_alpha_num)"""
s["omp_init_lock"] = "call omp_init_lock(lck)"
s["omp_set_lock"] = "call omp_set_lock(lck)"
s["omp_unset_lock"] = "call omp_unset_lock(lck)"
@ -50,16 +61,76 @@ def new_dict(openmp=True):
s["set_i_H_j_threshold"] = """
thresh = H_apply_threshold
"""
return s
self.data = s
def __setitem__(self,key,value):
self.data[key] = value
def __getitem__(self,key):
return self.data[key]
def create_h_apply(s):
def __repr__(self):
buffer = template
for key in s:
buffer = buffer.replace('$'+key, s[key])
print buffer
for key,value in self.data.items():
buffer = buffer.replace('$'+key, value)
return buffer
def set_perturbation(self,pert):
if self.perturbation is not None:
raise
self.perturbation = pert
if pert is not None:
self.data["parameters"] = ",sum_e_2_pert_in,sum_norm_pert_in,sum_H_pert_diag_in,N_st,Nint"
self.data["declarations"] = """
integer, intent(in) :: N_st,Nint
double precision, intent(inout) :: sum_e_2_pert_in(N_st)
double precision, intent(inout) :: sum_norm_pert_in(N_st)
double precision, intent(inout) :: sum_H_pert_diag_in
double precision :: sum_e_2_pert(N_st)
double precision :: sum_norm_pert(N_st)
double precision :: sum_H_pert_diag
double precision :: e_2_pert_buffer(N_st,size_max)
double precision :: coef_pert_buffer(N_st,size_max)
"""
self.data["size_max"] = "256"
self.data["initialization"] = """
sum_e_2_pert = sum_e_2_pert_in
sum_norm_pert = sum_norm_pert_in
sum_H_pert_diag = sum_H_pert_diag_in
"""
self.data["keys_work_unlocked"] += """
call perturb_buffer_%s(keys_out,key_idx,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert, &
sum_norm_pert,sum_H_pert_diag,N_st,Nint)
"""%(pert,)
self.data["finalization"] = """
sum_e_2_pert_in = sum_e_2_pert
sum_norm_pert_in = sum_norm_pert
sum_H_pert_diag_in = sum_H_pert_diag
"""
if self.openmp:
self.data["omp_set_lock"] = ""
self.data["omp_unset_lock"] = ""
self.data["omp_test_lock"] = ".False."
self.data["omp_parallel"] += """&
!$OMP SHARED(N_st,Nint) PRIVATE(e_2_pert_buffer,coef_pert_buffer) &
!$OMP REDUCTION(+:sum_e_2_pert, sum_norm_pert, sum_H_pert_diag)"""
def set_selection_pt2(self,pert):
if self.selection_pt2 is not None:
raise
self.set_perturbation(pert)
self.selection_pt2 = pert
if pert is not None:
self.data["size_max"] = str(1024*128)
self.data["keys_work_unlocked"] = """
e_2_pert_buffer = 0.d0
coef_pert_buffer = 0.d0
""" + self.data["keys_work_unlocked"]
self.data["keys_work_locked"] = """
call fill_H_apply_buffer_selection(key_idx,keys_out,e_2_pert_buffer,coef_pert_buffer,N_st,N_int)
"""
self.data["omp_test_lock"] = "omp_test_lock(lck)"
self.data["omp_set_lock"] = "call omp_set_lock(lck)"
self.data["omp_unset_lock"] = "call omp_unset_lock(lck)"

26
scripts/perturbation.py Executable file
View File

@ -0,0 +1,26 @@
#!/usr/bin/env python
import os
Pert_dir = os.environ["QPACKAGE_ROOT"]+"/src/Perturbation/"
perturbations = []
for filename in filter(lambda x: x.endswith(".irp.f"), os.listdir(Pert_dir)):
filename = Pert_dir+filename
file = open(filename,'r')
lines = file.readlines()
file.close()
for line in lines:
buffer = line.lower().lstrip().split()
if len(buffer) > 1:
if buffer[0] == "subroutine" and buffer[1].startswith("pt2_"):
p = (buffer[1].split('(')[0])[4:]
perturbations.append( p )
if __name__ == '__main__':
print 'Perturbations:'
for k in perturbations:
print '* ', k

10
src/CISD/H_apply.irp.f
View File

@ -54,14 +54,14 @@ subroutine H_apply_cisd
particle_mask(:,1) = iand(not(HF_bitmask(:,1)),full_ijkl_bitmask(:,1))
particle_mask(:,2) = iand(not(HF_bitmask(:,2)),full_ijkl_bitmask(:,2))
call H_apply_cisd_monoexc(HF_bitmask, &
call H_apply_cisd_OpenMP_monoexc(HF_bitmask, &
hole_mask, particle_mask)
call H_apply_cisd_diexc(HF_bitmask, &
call H_apply_cisd_OpenMP_diexc(HF_bitmask, &
hole_mask, particle_mask, &
hole_mask, particle_mask )
call copy_H_apply_buffer_to_wf
call copy_h_apply_buffer_to_wf
end

7
src/CISD/README.rst
View File

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

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

@ -3,6 +3,11 @@ program cisd
integer :: i,k
double precision, allocatable :: eigvalues(:),eigvectors(:,:)
PROVIDE ref_bitmask_energy
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
@ -17,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

12
src/CISD/h_apply.py
View File

@ -1,11 +1,11 @@
#!/usr/bin/env python
import generate_h_apply
from generate_h_apply import *
# H_apply
s = generate_h_apply.new_dict(openmp=True)
s["subroutine"] = "H_apply_cisd"
s["keys_work"] = "call fill_H_apply_buffer_cisd(key_idx,keys_out,N_int)"
generate_h_apply.create_h_apply(s)
s = H_apply("cisd",openmp=True)
s["keys_work"] += """
call fill_H_apply_buffer_cisd(key_idx,keys_out,N_int)
"""
print s

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

8
src/Hartree_Fock_AOs/README.rst → src/CIS_dressed/README.rst
View File

@ -1,8 +1,3 @@
===================
Hartree-Fock Module
===================
Needed Modules
==============
@ -12,13 +7,16 @@ Needed Modules
* `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
=============

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

22
src/Dets/H_apply.irp.f
View File

@ -33,11 +33,12 @@ subroutine resize_H_apply_buffer_det(new_size)
integer, intent(in) :: new_size
integer(bit_kind), allocatable :: buffer_det(:,:,:)
double precision, allocatable :: buffer_coef(:,:)
double precision, allocatable :: buffer_e2(:,:)
integer :: i,j,k
integer :: Ndet
ASSERT (new_size > 0)
allocate ( buffer_det(N_int,2,new_size), buffer_coef(new_size,N_states) )
allocate ( buffer_det(N_int,2,new_size), buffer_coef(new_size,N_states), buffer_e2(new_size,N_states) )
do i=1,min(new_size,H_apply_buffer_N_det)
do k=1,N_int
@ -48,9 +49,10 @@ subroutine resize_H_apply_buffer_det(new_size)
ASSERT (sum(popcnt(H_apply_buffer_det(:,2,i))) == elec_beta_num )
enddo
do k=1,N_states
do i=1,min(new_size,H_apply_buffer_N_det)
buffer_coef(i,k) = H_apply_buffer_coef(i,k)
enddo
do i=1,min(new_size,H_apply_buffer_N_det)
buffer_coef(i,k) = H_apply_buffer_coef(i,k)
buffer_e2(i,k) = H_apply_buffer_e2(i,k)
enddo
enddo
H_apply_buffer_size = new_size
@ -70,20 +72,23 @@ subroutine resize_H_apply_buffer_det(new_size)
do k=1,N_states
do i=1,H_apply_buffer_N_det
H_apply_buffer_coef(i,k) = buffer_coef(i,k)
H_apply_buffer_e2(i,k) = buffer_e2(i,k)
enddo
enddo
deallocate (buffer_det, buffer_coef)
SOFT_TOUCH H_apply_buffer_det H_apply_buffer_coef H_apply_buffer_N_det
deallocate (buffer_det, buffer_coef, buffer_e2)
SOFT_TOUCH H_apply_buffer_det H_apply_buffer_coef H_apply_buffer_N_det H_apply_buffer_e2
end
BEGIN_PROVIDER [ integer(bit_kind), H_apply_buffer_det,(N_int,2,H_apply_buffer_size) ]
&BEGIN_PROVIDER [ double precision, H_apply_buffer_coef,(H_apply_buffer_size,N_states) ]
&BEGIN_PROVIDER [ double precision, H_apply_buffer_e2,(H_apply_buffer_size,N_states) ]
&BEGIN_PROVIDER [ integer, H_apply_buffer_N_det ]
implicit none
BEGIN_DOC
! Buffer of determinants/coefficients for H_apply. Uninitialized. Filled by H_apply subroutines.
! Buffer of determinants/coefficients/perturbative energy for H_apply.
! Uninitialized. Filled by H_apply subroutines.
END_DOC
H_apply_buffer_N_det = 0
@ -148,8 +153,9 @@ subroutine copy_H_apply_buffer_to_wf
psi_coef(i+N_det_old,k) = H_apply_buffer_coef(i,k)
enddo
enddo
H_apply_buffer_N_det = 0
SOFT_TOUCH psi_det psi_coef
SOFT_TOUCH psi_det psi_coef H_apply_buffer_N_det H_apply_buffer_det H_apply_buffer_coef H_apply_buffer_e2
end

84
src/Dets/H_apply_template.f
View File

@ -2,6 +2,12 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2 $parame
use omp_lib
use bitmasks
implicit none
BEGIN_DOC
! Generate all double excitations of key_in using the bit masks of holes and
! particles.
! Assume N_int is already provided.
END_DOC
integer,parameter :: size_max = $size_max
$declarations
integer(omp_lock_kind) :: lck
integer(bit_kind),intent(in) :: key_in(N_int,2)
@ -21,21 +27,18 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2 $parame
integer :: N_elec_in_key_hole_1(2),N_elec_in_key_part_1(2)
integer :: N_elec_in_key_hole_2(2),N_elec_in_key_part_2(2)
integer,parameter :: size_max = $size_max
double precision :: hij_elec, mo_bielec_integral, thresh
integer, allocatable :: ia_ja_pairs(:,:,:)
double precision :: diag_H_mat_elem, E_ref
double precision :: diag_H_mat_elem
PROVIDE mo_integrals_map ref_bitmask_energy
PROVIDE mo_bielec_integrals_in_map
PROVIDE mo_integrals_map ref_bitmask_energy key_pattern_not_in_ref
PROVIDE mo_bielec_integrals_in_map reference_energy psi_ref_coef psi_ref
$set_i_H_j_threshold
$omp_init_lock
E_ref = diag_H_mat_elem(key_in,N_int)
$initialization
$omp_parallel
@ -153,8 +156,9 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2 $parame
hij_tab(key_idx) = hij_elec
ASSERT (key_idx <= size_max)
if (key_idx == size_max) then
$keys_work_unlocked
$omp_set_lock
$keys_work
$keys_work_locked
$omp_unset_lock
key_idx = 0
endif
@ -162,7 +166,8 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2 $parame
enddo
if (key_idx > ishft(size_max,-5)) then
if ($omp_test_lock) then
$keys_work
$keys_work_unlocked
$keys_work_locked
$omp_unset_lock
key_idx = 0
endif
@ -201,8 +206,9 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2 $parame
hij_tab(key_idx) = hij_elec
ASSERT (key_idx <= size_max)
if (key_idx == size_max) then
$keys_work_unlocked
$omp_set_lock
$keys_work
$keys_work_locked
$omp_unset_lock
key_idx = 0
endif
@ -210,7 +216,8 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2 $parame
enddo
if (key_idx > ishft(size_max,-5)) then
if ($omp_test_lock) then
$keys_work
$keys_work_locked
$keys_work_unlocked
$omp_unset_lock
key_idx = 0
endif
@ -219,8 +226,9 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2 $parame
enddo ! ii
$omp_enddo
enddo ! ispin
$keys_work_unlocked
$omp_set_lock
$keys_work
$keys_work_locked
$omp_unset_lock
deallocate (keys_out,hij_tab,ia_ja_pairs)
$omp_end_parallel
@ -233,6 +241,12 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1 $parameters )
use omp_lib
use bitmasks
implicit none
BEGIN_DOC
! Generate all single excitations of key_in using the bit masks of holes and
! particles.
! Assume N_int is already provided.
END_DOC
integer,parameter :: size_max = $size_max
$declarations
integer(omp_lock_kind) :: lck
integer(bit_kind),intent(in) :: key_in(N_int,2)
@ -252,21 +266,18 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1 $parameters )
integer :: N_elec_in_key_hole_1(2),N_elec_in_key_part_1(2)
integer :: N_elec_in_key_hole_2(2),N_elec_in_key_part_2(2)
integer,parameter :: size_max = $size_max
double precision :: hij_elec, thresh
integer, allocatable :: ia_ja_pairs(:,:,:)
double precision :: diag_H_mat_elem, E_ref
double precision :: diag_H_mat_elem
PROVIDE mo_integrals_map ref_bitmask_energy
PROVIDE mo_bielec_integrals_in_map
PROVIDE mo_integrals_map ref_bitmask_energy key_pattern_not_in_ref
PROVIDE mo_bielec_integrals_in_map reference_energy psi_ref_coef psi_ref
$set_i_H_j_threshold
$omp_init_lock
E_ref = diag_H_mat_elem(key_in,N_int)
$initialization
$omp_parallel
@ -305,7 +316,6 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1 $parameters )
integer(bit_kind) :: test(N_int,2)
double precision :: accu
accu = 0.d0
hij_elec = 0.d0
do ispin=1,2
other_spin = iand(ispin,1)+1
$omp_do
@ -319,33 +329,33 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1 $parameters )
k_a = ishft(j_a-1,-bit_kind_shift)+1
l_a = j_a-ishft(k_a-1,bit_kind_shift)-1
hole(k_a,ispin) = ibset(hole(k_a,ispin),l_a)
call i_H_j(hole,key_in,N_int,hij_elec)
if(dabs(hij_elec) .ge. thresh)then
key_idx += 1
do k=1,N_int
keys_out(k,1,key_idx) = hole(k,1)
keys_out(k,2,key_idx) = hole(k,2)
enddo
hij_tab(key_idx) = hij_elec
if (key_idx > ishft(size_max,-5)) then
if ($omp_test_lock) then
$keys_work
$omp_unset_lock
key_idx = 0
endif
endif
if (key_idx == size_max) then
$omp_set_lock
$keys_work
key_idx += 1
do k=1,N_int
keys_out(k,1,key_idx) = hole(k,1)
keys_out(k,2,key_idx) = hole(k,2)
enddo
hij_tab(key_idx) = hij_elec
if (key_idx > ishft(size_max,-5)) then
if ($omp_test_lock) then
$keys_work_unlocked
$keys_work_locked
$omp_unset_lock
key_idx = 0
endif
endif
if (key_idx == size_max) then
$keys_work_unlocked
$omp_set_lock
$keys_work_locked
$omp_unset_lock
key_idx = 0
endif
enddo ! ii
$omp_enddo
enddo ! ispin
$keys_work_unlocked
$omp_set_lock
$keys_work
$keys_work_locked
$omp_unset_lock
deallocate (keys_out,hij_tab,ia_ja_pairs)
$omp_end_parallel

18
src/Dets/README.rst
View File

@ -135,10 +135,10 @@ Documentation
`get_s2 <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/s2.irp.f#L1>`_
Returns <S^2>
`a_operator <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L840>`_
`a_operator <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L666>`_
Needed for diag_H_mat_elem
`ac_operator <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L885>`_
`ac_operator <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L711>`_
Needed for diag_H_mat_elem
`decode_exc <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L76>`_
@ -148,15 +148,9 @@ Documentation
s1,s2 : Spins (1:alpha, 2:beta)
degree : Degree of excitation
`diag_h_mat_elem <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L778>`_
`diag_h_mat_elem <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L604>`_
Computes <i|H|i>
`filter_connected <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L603>`_
Filters out the determinants that are not connected by H
`filter_connected_i_h_psi0 <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L687>`_
Undocumented
`get_double_excitation <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L141>`_
Returns the two excitation operators between two doubly excited determinants and the phase
@ -172,10 +166,10 @@ Documentation
`get_mono_excitation <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L274>`_
Returns the excitation operator between two singly excited determinants and the phase
`get_occ_from_key <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L933>`_
`get_occ_from_key <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L759>`_
Returns a list of occupation numbers from a bitstring
`h_u_0 <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L949>`_
`h_u_0 <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L775>`_
Computes v_0 = H|u_0>
.br
n : number of determinants
@ -185,7 +179,7 @@ Documentation
`i_h_j <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L355>`_
Returns <i|H|j> where i and j are determinants
`i_h_psim <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L491>`_
`i_h_psi <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L491>`_
Undocumented
`h_matrix_all_dets <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/utils.irp.f#L1>`_

256
src/Dets/connected_to_ref.irp.f Normal file
View File

@ -0,0 +1,256 @@
integer function connected_to_ref(key,keys,Nint,N_past_in,Ndet,thresh)
use bitmasks
implicit none
integer, intent(in) :: Nint, N_past_in, Ndet
integer(bit_kind), intent(in) :: keys(Nint,2,Ndet)
integer(bit_kind), intent(in) :: key(Nint,2)
double precision, intent(in) :: thresh
integer :: N_past
integer :: i, l
integer :: degree_x2
logical :: det_is_not_or_may_be_in_ref, t
double precision :: hij_elec
! output : 0 : not connected
! i : connected to determinant i of the past
! -i : is the ith determinant of the refernce wf keys
ASSERT (Nint > 0)
ASSERT (Nint == N_int)
connected_to_ref = 0
N_past = max(1,N_past_in)
if (Nint == 1) then
do i=N_past-1,1,-1
degree_x2 = popcnt(xor( key(1,1), keys(1,1,i))) + &
popcnt(xor( key(1,2), keys(1,2,i)))
if(degree_x2 == 0)then
connected_to_ref = -i
return
endif
if (degree_x2 > 5) then
cycle
else
call i_H_j(keys(1,1,i),key,N_int,hij_elec)
if(dabs(hij_elec).lt.thresh)cycle
connected_to_ref = i
return
endif
enddo
!DIR$ FORCEINLINE
t = det_is_not_or_may_be_in_ref(key,N_int)
if ( t ) then
return
endif
do i=N_past,Ndet
if ( (key(1,1) /= keys(1,1,i)).or. &
(key(1,2) /= keys(1,2,i)) ) then
cycle
endif
connected_to_ref = -i
return
enddo
return
else if (Nint==2) then
do i=N_past-1,1,-1
degree_x2 = popcnt(xor( key(1,1), keys(1,1,i))) + &
popcnt(xor( key(1,2), keys(1,2,i))) + &
popcnt(xor( key(2,1), keys(2,1,i))) + &
popcnt(xor( key(2,2), keys(2,2,i)))
if(degree_x2 == 0)then
connected_to_ref = -i
return
endif
if (degree_x2 > 5) then
cycle
else
call i_H_j(keys(1,1,i),key,N_int,hij_elec)
if(dabs(hij_elec).lt.thresh)cycle
connected_to_ref = i
return
endif
enddo
!DIR$ FORCEINLINE
t = det_is_not_or_may_be_in_ref(key,N_int)
if ( t ) then
return
endif
!DIR$ LOOP COUNT (1000)
do i=N_past+1,Ndet
if ( (key(1,1) /= keys(1,1,i)).or. &
(key(1,2) /= keys(1,2,i)).or. &
(key(2,1) /= keys(2,1,i)).or. &
(key(2,2) /= keys(2,2,i)) ) then
cycle
endif
connected_to_ref = -i
return
enddo
return
else if (Nint==3) then
do i=N_past-1,1,-1
degree_x2 = popcnt(xor( key(1,1), keys(1,1,i))) + &
popcnt(xor( key(1,2), keys(1,2,i))) + &
popcnt(xor( key(2,1), keys(2,1,i))) + &
popcnt(xor( key(2,2), keys(2,2,i))) + &
popcnt(xor( key(3,1), keys(3,1,i))) + &
popcnt(xor( key(3,2), keys(3,2,i)))
if(degree_x2 == 0)then
connected_to_ref = -i
return
endif
if (degree_x2 > 5) then
cycle
else
call i_H_j(keys(1,1,i),key,N_int,hij_elec)
if(dabs(hij_elec).lt.thresh)cycle
connected_to_ref = i
return
endif
enddo
!DIR$ FORCEINLINE
t = det_is_not_or_may_be_in_ref(key,N_int)
if ( t ) then
return
endif
do i=N_past+1,Ndet
if ( (key(1,1) /= keys(1,1,i)).or. &
(key(1,2) /= keys(1,2,i)).or. &
(key(2,1) /= keys(2,1,i)).or. &
(key(2,2) /= keys(2,2,i)).or. &
(key(3,1) /= keys(3,1,i)).or. &
(key(3,2) /= keys(3,2,i)) ) then
cycle
endif
connected_to_ref = -i
return
enddo
return
else
do i=N_past-1,1,-1
degree_x2 = popcnt(xor( key(1,1), keys(1,1,i))) + &
popcnt(xor( key(1,2), keys(1,2,i)))
!DEC$ LOOP COUNT MIN(3)
do l=2,Nint
degree_x2 = degree_x2 + popcnt(xor( key(l,1), keys(l,1,i))) +&
popcnt(xor( key(l,2), keys(l,2,i)))
enddo
if(degree_x2 == 0)then
connected_to_ref = -i
return
endif
if (degree_x2 > 5) then
cycle
else
call i_H_j(keys(1,1,i),key,N_int,hij_elec)
if(dabs(hij_elec).lt.thresh)cycle
connected_to_ref = i
return
endif
enddo
!DIR$ FORCEINLINE
t = det_is_not_or_may_be_in_ref(key,N_int)
if ( t ) then
return
endif
do i=N_past+1,Ndet
if ( (key(1,1) /= keys(1,1,i)).or. &
(key(1,2) /= keys(1,2,i)) ) then
cycle
else
connected_to_ref = -1
!DEC$ LOOP COUNT MIN(3)
do l=2,Nint
if ( (key(l,1) /= keys(l,1,i)).or. &
(key(l,2) /= keys(l,2,i)) ) then
connected_to_ref = 0
exit
endif
enddo
if (connected_to_ref /= 0) then
return
endif
endif
enddo
endif
end
logical function det_is_not_or_may_be_in_ref(key,Nint)
use bitmasks
implicit none
BEGIN_DOC
! If true, det is not in ref
! If false, det may be in ref
END_DOC
integer(bit_kind), intent(in) :: key(Nint,2), Nint
integer(bit_kind) :: key_int
integer*1 :: key_short(bit_kind)
!DIR$ ATTRIBUTES ALIGN : 32 :: key_short
equivalence (key_int,key_short)
integer :: i, ispin, k
det_is_not_or_may_be_in_ref = .False.
do ispin=1,2
do i=1,Nint
key_int = key(i,ispin)
do k=1,bit_kind
det_is_not_or_may_be_in_ref = &
det_is_not_or_may_be_in_ref .or. &
key_pattern_not_in_ref(key_short(k), i, ispin)
enddo
if(det_is_not_or_may_be_in_ref) then
return
endif
enddo
enddo
end
BEGIN_PROVIDER [ logical, key_pattern_not_in_ref, (-128:127,N_int,2) ]
use bitmasks
implicit none
BEGIN_DOC
! Min and max values of the integers of the keys of the reference
END_DOC
integer :: i, j, ispin
integer(bit_kind) :: key
integer*1 :: key_short(bit_kind)
equivalence (key,key_short)
integer :: idx, k
key_pattern_not_in_ref = .True.
do j=1,N_det
do ispin=1,2
do i=1,N_int
key = psi_det(i,ispin,j)
do k=1,bit_kind
key_pattern_not_in_ref( key_short(k), i, ispin ) = .False.
enddo
enddo
enddo
enddo
END_PROVIDER

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

@ -133,7 +133,6 @@ subroutine davidson_diag(dets_in,u_in,energies,dim_in,sze,N_st,Nint)
!$OMP END PARALLEL
do iter=1,davidson_sze_max-1
print *, 'iter = ',iter
! print *, '***************'
! do i=1,iter
@ -174,7 +173,6 @@ subroutine davidson_diag(dets_in,u_in,energies,dim_in,sze,N_st,Nint)
! -------------
call lapack_diag(lambda,y,h,N_st*davidson_sze_max,N_st*iter)
print *, lambda(1:4)
! Express eigenvectors of h in the determinant basis
! --------------------------------------------------
@ -203,13 +201,11 @@ subroutine davidson_diag(dets_in,u_in,energies,dim_in,sze,N_st,Nint)
enddo
residual_norm(k) = u_dot_u(R(1,k),sze)
enddo
print *, 'Lambda'
print *, lambda(1:N_st) + nuclear_repulsion
print *, 'Residual_norm'
print *, residual_norm(1:N_st)
print *, ''
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-5
converged = maxval(residual_norm) < 1.d-10
if (converged) then
exit
endif
@ -279,3 +275,4 @@ subroutine davidson_diag(dets_in,u_in,energies,dim_in,sze,N_st,Nint)
)
end

175
src/Dets/filter_connected.irp.f Normal file
View File

@ -0,0 +1,175 @@
subroutine filter_connected(key1,key2,Nint,sze,idx)
use bitmasks
implicit none
BEGIN_DOC
! Filters out the determinants that are not connected by H
END_DOC
integer, intent(in) :: Nint, sze
integer(bit_kind), intent(in) :: key1(Nint,2,sze)
integer(bit_kind), intent(in) :: key2(Nint,2)
integer, intent(out) :: idx(0:sze)
integer :: i,j,l
integer :: degree_x2
ASSERT (Nint > 0)
ASSERT (sze >= 0)
l=1
if (Nint==1) then
!DIR$ LOOP COUNT (1000)
do i=1,sze
degree_x2 = popcnt( xor( key1(1,1,i), key2(1,1))) &
+ popcnt( xor( key1(1,2,i), key2(1,2)))
if (degree_x2 < 5) then
idx(l) = i
l = l+1
endif
enddo
else if (Nint==2) then
!DIR$ LOOP COUNT (1000)
do i=1,sze
degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + &
popcnt(xor( key1(2,1,i), key2(2,1))) + &
popcnt(xor( key1(1,2,i), key2(1,2))) + &
popcnt(xor( key1(2,2,i), key2(2,2)))
if (degree_x2 < 5) then
idx(l) = i
l = l+1
endif
enddo
else if (Nint==3) then
!DIR$ LOOP COUNT (1000)
do i=1,sze
degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + &
popcnt(xor( key1(1,2,i), key2(1,2))) + &
popcnt(xor( key1(2,1,i), key2(2,1))) + &
popcnt(xor( key1(2,2,i), key2(2,2))) + &
popcnt(xor( key1(3,1,i), key2(3,1))) + &
popcnt(xor( key1(3,2,i), key2(3,2)))
if (degree_x2 < 5) then
idx(l) = i
l = l+1
endif
enddo
else
!DIR$ LOOP COUNT (1000)
do i=1,sze
degree_x2 = 0
!DEC$ LOOP COUNT MIN(4)
do j=1,Nint
degree_x2 = degree_x2+ popcnt(xor( key1(j,1,i), key2(j,1))) +&
popcnt(xor( key1(j,2,i), key2(j,2)))
if (degree_x2 > 4) then
exit
endif
enddo
if (degree_x2 <= 5) then
idx(l) = i
l = l+1
endif
enddo
endif
idx(0) = l-1
end
subroutine filter_connected_i_H_psi0(key1,key2,Nint,sze,idx)
use bitmasks
implicit none
integer, intent(in) :: Nint, sze
integer(bit_kind), intent(in) :: key1(Nint,2,sze)
integer(bit_kind), intent(in) :: key2(Nint,2)
integer, intent(out) :: idx(0:sze)
integer :: i,l
integer :: degree_x2
ASSERT (Nint > 0)
ASSERT (Nint == N_int)
ASSERT (sze > 0)
l=1
if (Nint==1) then
!DIR$ LOOP COUNT (1000)
do i=1,sze
degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + &
popcnt(xor( key1(1,2,i), key2(1,2)))
if (degree_x2 < 5) then
if(degree_x2 .ne. 0)then
idx(l) = i
l = l+1
endif
endif
enddo
else if (Nint==2) then
!DIR$ LOOP COUNT (1000)
do i=1,sze
degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + &
popcnt(xor( key1(2,1,i), key2(2,1))) + &
popcnt(xor( key1(1,2,i), key2(1,2))) + &
popcnt(xor( key1(2,2,i), key2(2,2)))
if (degree_x2 < 5) then
if(degree_x2 .ne. 0)then
idx(l) = i
l = l+1
endif
endif
enddo
else if (Nint==3) then
!DIR$ LOOP COUNT (1000)
do i=1,sze
degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + &
popcnt(xor( key1(1,2,i), key2(1,2))) + &
popcnt(xor( key1(2,1,i), key2(2,1))) + &
popcnt(xor( key1(2,2,i), key2(2,2))) + &
popcnt(xor( key1(3,1,i), key2(3,1))) + &
popcnt(xor( key1(3,2,i), key2(3,2)))
if (degree_x2 < 5) then
if(degree_x2 .ne. 0)then
idx(l) = i
l = l+1
endif
endif
enddo
else
!DIR$ LOOP COUNT (1000)
do i=1,sze
degree_x2 = 0
!DEC$ LOOP COUNT MIN(4)
do l=1,Nint
degree_x2 = degree_x2+ popcnt(xor( key1(l,1,i), key2(l,1))) +&
popcnt(xor( key1(l,2,i), key2(l,2)))
if (degree_x2 > 4) then
exit
endif
enddo
if (degree_x2 <= 5) then
if(degree_x2 .ne. 0)then
idx(l) = i
l = l+1
endif
endif
enddo
endif
idx(0) = l-1
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

221
src/Dets/slater_rules.irp.f
View File

@ -488,13 +488,13 @@ end
subroutine i_H_psim(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array)
subroutine i_H_psi(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array)
use bitmasks
implicit none
integer, intent(in) :: Nint, Ndet,Ndet_max,Nstate
integer, intent(in) :: keys(Nint,2,Ndet_max)
integer(bit_kind), intent(in) :: keys(Nint,2,Ndet)
integer(bit_kind), intent(in) :: key(Nint,2)
double precision, intent(in) :: coef(Ndet_max,Nstate)
integer, intent(in) :: key(Nint,2)
double precision, intent(out) :: i_H_psi_array(Nstate)
integer :: i, ii,j
@ -503,6 +503,11 @@ subroutine i_H_psim(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array)
double precision :: hij
integer :: idx(0:Ndet)
ASSERT (Nint > 0)
ASSERT (N_int == Nint)
ASSERT (Nstate > 0)
ASSERT (Ndet > 0)
ASSERT (Ndet_max >= Ndet)
i_H_psi_array = 0.d0
call filter_connected_i_H_psi0(keys,key,Nint,Ndet,idx)
do ii=1,idx(0)
@ -512,6 +517,7 @@ subroutine i_H_psim(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array)
do j = 1, Nstate
i_H_psi_array(j) = i_H_psi_array(j) + coef(i,j)*hij
enddo
print *, 'x', coef(i,1), hij, i_H_psi_array(1)
enddo
end
@ -600,180 +606,6 @@ end
subroutine filter_connected(key1,key2,Nint,sze,idx)
use bitmasks
implicit none
BEGIN_DOC
! Filters out the determinants that are not connected by H
END_DOC
integer, intent(in) :: Nint, sze
integer(bit_kind), intent(in) :: key1(Nint,2,sze)
integer(bit_kind), intent(in) :: key2(Nint,2)
integer, intent(out) :: idx(0:sze)
integer :: i,j,l
integer :: degree_x2
ASSERT (Nint > 0)
ASSERT (sze > 0)
l=1
if (Nint==1) then
!DIR$ LOOP COUNT (1000)
do i=1,sze
degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + &
popcnt(xor( key1(1,2,i), key2(1,2)))
if (degree_x2 < 5) then
idx(l) = i
l = l+1
endif
enddo
else if (Nint==2) then
!DIR$ LOOP COUNT (1000)
do i=1,sze
degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + &
popcnt(xor( key1(1,2,i), key2(1,2))) + &
popcnt(xor( key1(2,1,i), key2(2,1))) + &
popcnt(xor( key1(2,2,i), key2(2,2)))
if (degree_x2 < 5) then
idx(l) = i
l = l+1
endif
enddo
else if (Nint==3) then
!DIR$ LOOP COUNT (1000)
do i=1,sze
degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + &
popcnt(xor( key1(1,2,i), key2(1,2))) + &
popcnt(xor( key1(2,1,i), key2(2,1))) + &
popcnt(xor( key1(2,2,i), key2(2,2))) + &
popcnt(xor( key1(3,1,i), key2(3,1))) + &
popcnt(xor( key1(3,2,i), key2(3,2)))
if (degree_x2 < 5) then
idx(l) = i
l = l+1
endif
enddo
else
!DIR$ LOOP COUNT (1000)
do i=1,sze
degree_x2 = 0
!DEC$ LOOP COUNT MIN(4)
do j=1,Nint
degree_x2 = degree_x2+ popcnt(xor( key1(j,1,i), key2(j,1))) +&
popcnt(xor( key1(j,2,i), key2(j,2)))
if (degree_x2 > 4) then
exit
endif
enddo
if (degree_x2 <= 5) then
idx(l) = i
l = l+1
endif
enddo
endif
idx(0) = l-1
end
subroutine filter_connected_i_H_psi0(key1,key2,Nint,sze,idx)
use bitmasks
implicit none
integer, intent(in) :: Nint, sze
integer(bit_kind), intent(in) :: key1(Nint,2,sze)
integer(bit_kind), intent(in) :: key2(Nint,2)
integer, intent(out) :: idx(0:sze)
integer :: i,l
integer :: degree_x2
ASSERT (Nint > 0)
ASSERT (sze > 0)
l=1
if (Nint==1) then
!DIR$ LOOP COUNT (1000)
do i=1,sze
degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + &
popcnt(xor( key1(1,2,i), key2(1,2)))
if (degree_x2 < 5) then
if(degree_x2 .ne. 0)then
idx(l) = i
l = l+1
endif
endif
enddo
else if (Nint==2) then
!DIR$ LOOP COUNT (1000)
do i=1,sze
degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + &
popcnt(xor( key1(1,2,i), key2(1,2))) + &
popcnt(xor( key1(2,1,i), key2(2,1))) + &
popcnt(xor( key1(2,2,i), key2(2,2)))
if (degree_x2 < 5) then
if(degree_x2 .ne. 0)then
idx(l) = i
l = l+1
endif
endif
enddo
else if (Nint==3) then
!DIR$ LOOP COUNT (1000)
do i=1,sze
degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + &
popcnt(xor( key1(1,2,i), key2(1,2))) + &
popcnt(xor( key1(2,1,i), key2(2,1))) + &
popcnt(xor( key1(2,2,i), key2(2,2))) + &
popcnt(xor( key1(3,1,i), key2(3,1))) + &
popcnt(xor( key1(3,2,i), key2(3,2)))
if (degree_x2 < 5) then
if(degree_x2 .ne. 0)then
idx(l) = i
l = l+1
endif
endif
enddo
else
!DIR$ LOOP COUNT (1000)
do i=1,sze
degree_x2 = 0
!DEC$ LOOP COUNT MIN(4)
do l=1,Nint
degree_x2 = degree_x2+ popcnt(xor( key1(l,1,i), key2(l,1))) +&
popcnt(xor( key1(l,2,i), key2(l,2)))
if (degree_x2 > 4) then
exit
endif
enddo
if (degree_x2 <= 5) then
if(degree_x2 .ne. 0)then
idx(l) = i
l = l+1
endif
endif
enddo
endif
idx(0) = l-1
end
double precision function diag_H_mat_elem(det_in,Nint)
implicit none
@ -963,6 +795,7 @@ subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint)
integer(bit_kind),intent(in) :: keys_tmp(Nint,2,n)
integer, allocatable :: idx(:)
double precision :: hij
double precision, allocatable :: vt(:)
integer :: i,j,k,l, jj
integer :: i0, j0
ASSERT (Nint > 0)
@ -971,32 +804,34 @@ subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint)
PROVIDE ref_bitmask_energy
integer, parameter :: block_size = 157
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE(i,hij,j,k,idx,jj) SHARED(n,H_jj,u_0,keys_tmp,Nint)&
!$OMP PRIVATE(i,hij,j,k,idx,jj,vt) SHARED(n,H_jj,u_0,keys_tmp,Nint)&
!$OMP SHARED(v_0)
allocate(idx(0:block_size))
!$OMP DO SCHEDULE(static)
do i=1,n
v_0(i) = H_jj(i) * u_0(i)
enddo
!$OMP END DO
allocate(idx(0:n), vt(n))
Vt = 0.d0
!$OMP DO SCHEDULE(guided)
do i0=1,n,block_size
do j0=1,n,block_size
do i=i0,min(i0+block_size-1,n)
call filter_connected(keys_tmp(1,1,j0),keys_tmp(1,1,i),Nint,min(block_size,i-j0+1),idx)
do i=1,n
call filter_connected(keys_tmp(1,1,1),keys_tmp(1,1,i),Nint,i-1,idx)
do jj=1,idx(0)
j = idx(jj)+j0-1
if ( (j<i).and.(dabs(u_0(j)) > 1.d-8)) then
call i_H_j(keys_tmp(1,1,j),keys_tmp(1,1,i),Nint,hij)
v_0(i) = v_0(i) + hij*u_0(j)
v_0(j) = v_0(j) + hij*u_0(i)
endif
j = idx(jj)
call i_H_j(keys_tmp(1,1,j),keys_tmp(1,1,i),Nint,hij)
vt (i) = vt (i) + hij*u_0(j)
vt (j) = vt (j) + hij*u_0(i)
enddo
enddo
enddo
enddo
!$OMP END DO NOWAIT
deallocate(idx)
!$OMP END DO
!$OMP CRITICAL
do i=1,n
v_0(i) = v_0(i) + vt(i)
enddo
!$OMP END CRITICAL
deallocate(idx,vt)
!$OMP END PARALLEL
end

6
src/MOGuess/README.rst
View File

@ -23,7 +23,8 @@ Documentation
.. NEEDED_MODULES file.
`h_core_guess <http://github.com/LCPQ/quantum_package/tree/master/src/MOGuess/H_CORE_guess.irp.f#L1>`_
None
Undocumented
`ao_ortho_lowdin_coef <http://github.com/LCPQ/quantum_package/tree/master/src/MOGuess/mo_ortho_lowdin.irp.f#L2>`_
matrix of the coefficients of the mos generated by the
orthonormalization by the S^{-1/2} canonical transformation of the aos
@ -34,6 +35,7 @@ None
supposed to be the Identity
`ao_ortho_lowdin_nucl_elec_integral <http://github.com/LCPQ/quantum_package/tree/master/src/MOGuess/pot_mo_ortho_lowdin_ints.irp.f#L1>`_
None
Undocumented

29
src/Makefile.config.example
View File

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

167
src/MonoInts/README.rst
View File

@ -42,13 +42,17 @@ Documentation
: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
Undocumented
`do_print <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/check_orthonormality.irp.f#L11>`_
None
Undocumented
`n_pt_max_i_x <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/dimensions.irp.f#L2>`_
None
Undocumented
`n_pt_max_integrals <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/dimensions.irp.f#L1>`_
None
Undocumented
`ao_deriv2_x <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/kin_ao_ints.irp.f#L1>`_
second derivatives matrix elements in the ao basis
.. math::
@ -67,46 +71,169 @@ None
.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#L126>`_
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
Undocumented
`mo_mono_elec_integral <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/mo_mono_ints.irp.f#L35>`_
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
Undocumented
`orthonormalize_mos <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/orthonormalize.irp.f#L1>`_
None
Undocumented
`ao_nucl_elec_integral <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L1>`_
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
Undocumented
`i_x1_pol_mult_mono_elec <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L285>`_
None
Undocumented
`i_x2_pol_mult_mono_elec <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L357>`_
None
Undocumented
`int_gaus_pol <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L428>`_
None
Undocumented
`nai_pol_mult <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L82>`_
None
Undocumented
`v_e_n <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L409>`_
None
Undocumented
`v_phi <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L473>`_
None
Undocumented
`v_r <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L457>`_
None
Undocumented
`v_theta <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L486>`_
None
Undocumented
`wallis <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_ao_ints.irp.f#L502>`_
None
Undocumented
`mo_nucl_elec_integral <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/pot_mo_ints.irp.f#L1>`_
None
Undocumented
`save_ortho_mos <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/save_ortho_mos.irp.f#L1>`_
None
Undocumented
`ao_deriv_1_x <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/spread_dipole_ao.irp.f#L148>`_
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
`ao_deriv_1_y <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/spread_dipole_ao.irp.f#L149>`_
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
`ao_deriv_1_z <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/spread_dipole_ao.irp.f#L150>`_
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
`ao_dipole_x <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/spread_dipole_ao.irp.f#L75>`_
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
`ao_dipole_y <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/spread_dipole_ao.irp.f#L76>`_
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
`ao_dipole_z <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/spread_dipole_ao.irp.f#L77>`_
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
`ao_spread_x <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/spread_dipole_ao.irp.f#L1>`_
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
`ao_spread_y <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/spread_dipole_ao.irp.f#L2>`_
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
`ao_spread_z <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/spread_dipole_ao.irp.f#L3>`_
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
`overlap_bourrin_deriv_x <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/spread_dipole_ao.irp.f#L359>`_
Undocumented
`overlap_bourrin_dipole <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/spread_dipole_ao.irp.f#L318>`_
Undocumented
`overlap_bourrin_spread <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/spread_dipole_ao.irp.f#L265>`_
Undocumented
`overlap_bourrin_x <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/spread_dipole_ao.irp.f#L374>`_
Undocumented
`overlap_bourrin_x_abs <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/spread_dipole_ao.irp.f#L226>`_
Undocumented
`power <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/spread_dipole_ao.irp.f#L310>`_
Undocumented
`mo_deriv_1_x <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/spread_dipole_mo.irp.f#L69>`_
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
`mo_deriv_1_y <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/spread_dipole_mo.irp.f#L70>`_
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
`mo_deriv_1_z <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/spread_dipole_mo.irp.f#L71>`_
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
`mo_dipole_x <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/spread_dipole_mo.irp.f#L1>`_
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
`mo_dipole_y <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/spread_dipole_mo.irp.f#L2>`_
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
`mo_dipole_z <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/spread_dipole_mo.irp.f#L3>`_
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
`mo_spread_x <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/spread_dipole_mo.irp.f#L36>`_
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
`mo_spread_y <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/spread_dipole_mo.irp.f#L37>`_
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
`mo_spread_z <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts/spread_dipole_mo.irp.f#L38>`_
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

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

2
src/NEEDED_MODULES
View File

@ -1 +1 @@
AOs Bitmask Electrons Ezfio_files MOs Nuclei Output Utils Hartree_Fock BiInts MonoInts MOGuess Dets DensityMatrix CISD
AOs Bitmask Electrons Ezfio_files MOs Nuclei Output Utils Hartree_Fock BiInts MonoInts MOGuess Dets DensityMatrix CISD Perturbation Selection

6
src/Perturbation/ASSUMPTIONS.rst Normal file
View File

@ -0,0 +1,6 @@
* This is not allowed:
subroutine &
pt2_....

8
src/Perturbation/Makefile Normal file
View File

@ -0,0 +1,8 @@
default: all
# Define here all new external source files and objects.Don't forget to prefix the
# object files with IRPF90_temp/
SRC=perturbation_template.f
OBJ=
include $(QPACKAGE_ROOT)/src/Makefile.common

1
src/Perturbation/NEEDED_MODULES Normal file
View File

@ -0,0 +1 @@
AOs BiInts Bitmask Dets Electrons Ezfio_files Hartree_Fock MonoInts MOs Nuclei Output Utils

125
src/Perturbation/README.rst Normal file
View File

@ -0,0 +1,125 @@
===================
Perturbation Module
===================
All subroutines in `*.irp.f` starting with ``pt2_`` in the current directory are
perturbation computed using the routine ``i_H_psi``. Other cases are not allowed.
The arguments of the ``pt2_`` are always:
subroutine pt2_...( &
psi_ref, &
psi_ref_coefs, &
E_refs, &
det_pert, &
c_pert, &
e_2_pert, &
H_pert_diag, &
Nint, &
ndet, &
n_st )
integer, intent(in) :: Nint,ndet,n_st
integer(bit_kind), intent(in) :: psi_ref(Nint,2,ndet)
double precision , intent(in) :: psi_ref_coefs(ndet,n_st)
double precision , intent(in) :: E_refs(n_st)
integer(bit_kind), intent(in) :: det_pert(Nint,2)
double precision , intent(out) :: c_pert(n_st),e_2_pert(n_st),H_pert_diag
psi_ref
bitstring of the determinants present in the various n_st states
psi_ref_coefs
coefficients of the determinants on the various n_st states
E_refs
Energy of the various n_st states
det_pert
Perturber determinant
c_pert
Pertrubative coefficients for the various states
e_2_pert
Perturbative energetic contribution for the various states
H_pert_diag
Diagonal H matrix element of the perturber
Nint
Should be equal to N_int
Ndet
Number of determinants `i` in Psi on which we apply <det_pert|Hi>
N_st
Number of states
Assumptions
===========
.. Do not edit this section. It was auto-generated from the
.. NEEDED_MODULES file.
* This is not allowed:
subroutine &
pt2_....
Documentation
=============
.. Do not edit this section. It was auto-generated from the
.. NEEDED_MODULES file.
`pt2_epstein_nesbet <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/epstein_nesbet.irp.f#L1>`_
compute the standard Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
.br
for the various n_st states.
.br
c_pert(i) = <psi(i)|H|det_pert>/( E(i) - <det_pert|H|det_pert> )
.br
e_2_pert(i) = <psi(i)|H|det_pert>^2/( E(i) - <det_pert|H|det_pert> )
.br
`pt2_epstein_nesbet_2x2 <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/epstein_nesbet.irp.f#L33>`_
compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution
.br
for the various n_st states.
.br
e_2_pert(i) = 0.5 * (( <det_pert|H|det_pert> - E(i) ) - sqrt( ( <det_pert|H|det_pert> - E(i)) ^2 + 4 <psi(i)|H|det_pert>^2 )
.br
c_pert(i) = e_2_pert(i)/ <psi(i)|H|det_pert>
.br
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>`_

68
src/Perturbation/epstein_nesbet.irp.f Normal file
View File

@ -0,0 +1,68 @@
subroutine pt2_epstein_nesbet(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_st)
use bitmasks
implicit none
integer, intent(in) :: Nint,ndet,n_st
integer(bit_kind), intent(in) :: det_pert(Nint,2)
double precision , intent(out) :: c_pert(n_st),e_2_pert(n_st),H_pert_diag
double precision :: i_H_psi_array(N_st)
BEGIN_DOC
! compute the standard Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
!
! for the various n_st states.
!
! c_pert(i) = <psi(i)|H|det_pert>/( E(i) - <det_pert|H|det_pert> )
!
! e_2_pert(i) = <psi(i)|H|det_pert>^2/( E(i) - <det_pert|H|det_pert> )
!
END_DOC
integer :: i,j
double precision :: diag_H_mat_elem
ASSERT (Nint == N_int)
ASSERT (Nint > 0)
call i_H_psi(det_pert,psi_ref,psi_ref_coef,Nint,ndet,psi_ref_size,n_st,i_H_psi_array)
H_pert_diag = diag_H_mat_elem(det_pert,Nint)
do i =1,n_st
c_pert(i) = i_H_psi_array(i) / (reference_energy(i) - H_pert_diag)
e_2_pert(i) = c_pert(i) * i_H_psi_array(i)
enddo
end
subroutine pt2_epstein_nesbet_2x2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_st)
use bitmasks
implicit none
integer, intent(in) :: Nint,ndet,n_st
integer(bit_kind), intent(in) :: det_pert(Nint,2)
double precision , intent(out) :: c_pert(n_st),e_2_pert(n_st),H_pert_diag
double precision :: i_H_psi_array(N_st)
BEGIN_DOC
! compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution
!
! for the various n_st states.
!
! e_2_pert(i) = 0.5 * (( <det_pert|H|det_pert> - E(i) ) - sqrt( ( <det_pert|H|det_pert> - E(i)) ^2 + 4 <psi(i)|H|det_pert>^2 )
!
! c_pert(i) = e_2_pert(i)/ <psi(i)|H|det_pert>
!
END_DOC
integer :: i,j
double precision :: diag_H_mat_elem,delta_e
ASSERT (Nint == N_int)
ASSERT (Nint > 0)
print *, 'coefs'
print *, psi_ref_coef(1:N_det_ref,1)
print *, '-----'
call i_H_psi(det_pert,psi_ref,psi_ref_coef,Nint,N_det_ref,psi_ref_size,n_st,i_H_psi_array)
H_pert_diag = diag_H_mat_elem(det_pert,Nint)
do i =1,n_st
delta_e = H_pert_diag - reference_energy(i)
e_2_pert(i) = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * i_H_psi_array(i) * i_H_psi_array(i)))
c_pert(i) = e_2_pert(i)/i_H_psi_array(i)
enddo
print *, e_2_pert, delta_e , i_H_psi_array
end

13
src/Perturbation/perturbation.irp.f Normal file
View File

@ -0,0 +1,13 @@
BEGIN_SHELL [ /usr/bin/env python ]
from perturbation import perturbations
import os
filename = os.environ["QPACKAGE_ROOT"]+"/src/Perturbation/perturbation_template.f"
file = open(filename,'r')
template = file.read()
file.close()
for p in perturbations:
print template.replace("$PERT",p)
END_SHELL

47
src/Perturbation/perturbation_template.f Normal file
View File

@ -0,0 +1,47 @@
BEGIN_SHELL [ /usr/bin/env python ]
import perturbation
END_SHELL
subroutine perturb_buffer_$PERT(buffer,buffer_size,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert,sum_norm_pert,sum_H_pert_diag,N_st,Nint)
implicit none
BEGIN_DOC
! Applly pertubration ``$PERT`` to the buffer of determinants generated in the H_apply
! routine.
END_DOC
integer, intent(in) :: Nint, N_st, buffer_size
integer(bit_kind), intent(in) :: buffer(Nint,2,buffer_size)
double precision, intent(inout) :: sum_norm_pert(N_st),sum_e_2_pert(N_st)
double precision, intent(inout) :: coef_pert_buffer(N_st,buffer_size),e_2_pert_buffer(N_st,buffer_size),sum_H_pert_diag(N_st)
double precision :: c_pert(N_st), e_2_pert(N_st), H_pert_diag
integer :: i,k, c_ref
integer :: connected_to_ref
ASSERT (Nint > 0)
ASSERT (Nint == N_int)
ASSERT (buffer_size >= 0)
ASSERT (minval(sum_norm_pert) >= 0.d0)
ASSERT (N_st > 0)
do i = 1,buffer_size
c_ref = connected_to_ref(buffer(1,1,i),psi_det,Nint,N_det_ref,N_det,h_apply_threshold)
if (c_ref /= 0) then
cycle
endif
call pt2_$PERT(buffer(1,1,i), &
c_pert,e_2_pert,H_pert_diag,Nint,n_det_ref,n_st)
do k = 1,N_st
e_2_pert_buffer(k,i) = e_2_pert(k)
coef_pert_buffer(k,i) = c_pert(k)
sum_norm_pert(k) += c_pert(k) * c_pert(k)
sum_e_2_pert(k) += e_2_pert(k)
sum_H_pert_diag(k) += c_pert(k) * c_pert(k) * H_pert_diag
enddo
enddo
end

9
src/README.rst
View File

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

9
src/Utils/LinearAlgebra.irp.f
View File

@ -76,6 +76,7 @@ subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m)
!$OMP END DO NOWAIT
enddo
!$OMP BARRIER
!$OMP DO
do j=1,n
do i=1,m
@ -181,6 +182,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 +197,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