mirror of
https://github.com/LCPQ/quantum_package
synced 2024-11-12 17:13:54 +01:00
merged
mend
This commit is contained in:
Manu
commit
f55948890a
@ -14,10 +14,17 @@ keys_work
|
||||
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.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, &
|
||||
@ -29,7 +36,7 @@ def new_dict(openmp=True):
|
||||
!$OMP N_elec_in_key_hole_2,ia_ja_pairs) &
|
||||
!$OMP SHARED(key_in,N_int,elec_num_tab, &
|
||||
!$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 +57,54 @@ 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):
|
||||
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
|
||||
"""
|
||||
self.data["size_max"] = "256"
|
||||
self.data["initialization"] = """
|
||||
E_ref = diag_H_mat_elem(key_in,N_int)
|
||||
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"] += """
|
||||
call perturb_buffer_%s(keys_out,key_idx,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_test_lock"] = ".False."
|
||||
self.data["omp_parallel"] += """&
|
||||
!$OMP SHARED(N_st) &
|
||||
!$OMP REDUCTION(+:sum_e_2_pert, sum_norm_pert, sum_H_pert_diag)"""
|
||||
|
||||
|
||||
|
||||
|
||||
26
scripts/perturbation.py
Executable file
26
scripts/perturbation.py
Executable 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
|
||||
@ -126,8 +126,12 @@ def update_documentation(data):
|
||||
inside = line.startswith(".SH Description")
|
||||
else:
|
||||
if line.startswith(".SH"):
|
||||
return "".join(result)
|
||||
result.append(" "+line.strip()+"\n")
|
||||
break
|
||||
result.append(" "+line.strip())
|
||||
|
||||
if result == []:
|
||||
result = [" Undocumented"]
|
||||
return "\n".join(result)+'\n'
|
||||
|
||||
|
||||
|
||||
|
||||
@ -85,8 +85,10 @@ Documentation
|
||||
Number of primitives per atomic orbital
|
||||
|
||||
`ao_prim_num_max <http://github.com/LCPQ/quantum_package/tree/master/src/AOs/aos.irp.f#L176>`_
|
||||
None
|
||||
Undocumented
|
||||
|
||||
`ao_prim_num_max_align <http://github.com/LCPQ/quantum_package/tree/master/src/AOs/aos.irp.f#L177>`_
|
||||
None
|
||||
Undocumented
|
||||
|
||||
|
||||
|
||||
|
||||
33
src/AOs/tests/Makefile
Normal file
33
src/AOs/tests/Makefile
Normal file
@ -0,0 +1,33 @@
|
||||
OPENMP =1
|
||||
PROFILE =0
|
||||
DEBUG = 0
|
||||
|
||||
IRPF90+= -I tests
|
||||
|
||||
REF_FILES=$(subst %.irp.f, %.ref, $(wildcard *.irp.f))
|
||||
|
||||
.PHONY: clean executables serial_tests parallel_tests
|
||||
|
||||
all: clean executables serial_tests parallel_tests
|
||||
|
||||
parallel_tests: $(REF_FILES)
|
||||
@echo ; echo " ---- Running parallel tests ----" ; echo
|
||||
@OMP_NUM_THREADS=10 ${QPACKAGE_ROOT}/scripts/run_tests.py
|
||||
|
||||
serial_tests: $(REF_FILES)
|
||||
@echo ; echo " ---- Running serial tests ----" ; echo
|
||||
@OMP_NUM_THREADS=1 ${QPACKAGE_ROOT}/scripts/run_tests.py
|
||||
|
||||
executables: $(wildcard *.irp.f) veryclean
|
||||
$(MAKE) -C ..
|
||||
|
||||
%.ref: $(wildcard $(QPACKAGE_ROOT)/data/inputs/*.md5) executables
|
||||
$(QPACKAGE_ROOT)/scripts/create_test_ref.sh $*
|
||||
|
||||
clean:
|
||||
$(MAKE) -C .. clean
|
||||
|
||||
veryclean:
|
||||
$(MAKE) -C .. veryclean
|
||||
|
||||
|
||||
@ -107,7 +107,8 @@ Documentation
|
||||
AO integrals
|
||||
|
||||
`bielec_integrals_index <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/map_integrals.irp.f#L17>`_
|
||||
None
|
||||
Undocumented
|
||||
|
||||
`clear_ao_map <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts/map_integrals.irp.f#L128>`_
|
||||
Frees the memory of the AO map
|
||||
|
||||
|
||||
@ -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
|
||||
|
||||
|
||||
|
||||
@ -32,5 +32,6 @@ Documentation
|
||||
|
||||
`cisd <http://github.com/LCPQ/quantum_package/tree/master/src/CISD/cisd.irp.f#L/subroutine cisd/;">`_
|
||||
None
|
||||
=======
|
||||
|
||||
|
||||
|
||||
@ -1,19 +1,26 @@
|
||||
program cisd
|
||||
implicit none
|
||||
integer :: i
|
||||
call H_apply_cisd
|
||||
double precision, allocatable :: eigvalues(:),eigvectors(:,:)
|
||||
allocate(eigvalues(n_det),eigvectors(n_det,n_det))
|
||||
print *, 'N_det = ', N_det
|
||||
call lapack_diag(eigvalues,eigvectors,H_matrix_all_dets,n_det,n_det)
|
||||
integer :: i,k
|
||||
double precision, allocatable :: eigvalues(:),eigvectors(:,:)
|
||||
PROVIDE ref_bitmask_energy
|
||||
|
||||
double precision :: pt2(10), norm_pert(10), H_pert_diag
|
||||
|
||||
call H_apply_cisd
|
||||
allocate(eigvalues(n_states),eigvectors(n_det,n_states))
|
||||
print *, 'N_det = ', N_det
|
||||
print *, 'N_states = ', N_states
|
||||
psi_coef = - 1.d-4
|
||||
do k=1,N_states
|
||||
psi_coef(k,k) = 1.d0
|
||||
enddo
|
||||
call davidson_diag(psi_det,psi_coef,eigvalues,size(psi_coef,1),N_det,N_states,N_int)
|
||||
|
||||
! print *, H_matrix_all_dets
|
||||
print *, '---'
|
||||
print *, 'HF:', HF_energy
|
||||
print *, '---'
|
||||
do i = 1,3
|
||||
do i = 1,1
|
||||
print *, 'energy(i) = ',eigvalues(i) + nuclear_repulsion
|
||||
enddo
|
||||
! print *, eigvectors(:,1)
|
||||
deallocate(eigvalues,eigvectors)
|
||||
end
|
||||
|
||||
@ -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
|
||||
|
||||
|
||||
|
||||
259
src/DensityMatrix/density_matrix.irp.f
Normal file
259
src/DensityMatrix/density_matrix.irp.f
Normal file
@ -0,0 +1,259 @@
|
||||
use bitmasks
|
||||
BEGIN_PROVIDER [ integer, iunit_two_body_dm_aa ]
|
||||
&BEGIN_PROVIDER [ integer, iunit_two_body_dm_ab ]
|
||||
&BEGIN_PROVIDER [ integer, iunit_two_body_dm_bb ]
|
||||
implicit none
|
||||
use bitmasks
|
||||
BEGIN_DOC
|
||||
! Temporary files for 2-body dm calculation
|
||||
END_DOC
|
||||
integer :: getUnitAndOpen
|
||||
|
||||
iunit_two_body_dm_aa = getUnitAndOpen(trim(ezfio_filename)//'/work/two_body_aa.tmp','w')
|
||||
iunit_two_body_dm_ab = getUnitAndOpen(trim(ezfio_filename)//'/work/two_body_ab.tmp','w')
|
||||
iunit_two_body_dm_bb = getUnitAndOpen(trim(ezfio_filename)//'/work/two_body_bb.tmp','w')
|
||||
! Compute two body DM in file
|
||||
integer :: k,l,degree, idx,i
|
||||
integer :: exc(0:2,2,2),n_occ_alpha
|
||||
double precision :: phase, coef
|
||||
integer :: h1,h2,p1,p2,s1,s2
|
||||
double precision :: ck, cl
|
||||
character*(128), parameter :: f = '(i8,4(x,i5),x,d16.8)'
|
||||
do k=1,det_num
|
||||
ck = (det_coef_provider(k)+det_coef_provider(k))
|
||||
do l=1,k-1
|
||||
cl = det_coef_provider(l)
|
||||
call get_excitation_degree(det_provider(1,1,k),det_provider(1,1,l),degree,N_int)
|
||||
if (degree == 2) then
|
||||
call get_double_excitation(det_provider(1,1,k),det_provider(1,1,l),exc,phase,N_int)
|
||||
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
|
||||
call bielec_integrals_index(h1,h2,p1,p2,idx)
|
||||
ckl = phase*ck*cl
|
||||
select case (s1+s2)
|
||||
case(2) ! alpha alpha
|
||||
write(iunit_two_body_dm_aa,f) idx, h1,h2,p1,p2, ckl
|
||||
call bielec_integrals_index(h1,h2,p2,p1,idx)
|
||||
write(iunit_two_body_dm_aa,f) idx, h1,h2,p2,p1, -ckl
|
||||
case(3) ! alpha beta
|
||||
write(iunit_two_body_dm_ab,f) idx, h1,h2,p1,p2, ckl
|
||||
case(4) ! beta beta
|
||||
write(iunit_two_body_dm_bb,f) idx, h1,h2,p1,p2, ckl
|
||||
call bielec_integrals_index(h1,h2,p2,p1,idx)
|
||||
write(iunit_two_body_dm_bb,f) idx, h1,h2,p2,p1, -ckl
|
||||
end select
|
||||
else if (degree == 1) then
|
||||
call get_mono_excitation(det_provider(1,1,k),det_provider(1,1,l),exc,phase,N_int)
|
||||
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
|
||||
double precision :: ckl
|
||||
ckl = phase*ck*cl
|
||||
call bitstring_to_list(det_provider(1,1,k), occ(1,1), n_occ_alpha, N_int)
|
||||
call bitstring_to_list(det_provider(1,2,k), occ(1,2), n_occ_alpha, N_int)
|
||||
select case (s1)
|
||||
case (1) ! Alpha single excitation
|
||||
integer :: occ(N_int*bit_kind_size,2)
|
||||
do i = 1, elec_alpha_num
|
||||
p2=occ(i,1)
|
||||
h2=p2
|
||||
call bielec_integrals_index(h1,h2,p1,p2,idx)
|
||||
write(iunit_two_body_dm_aa,f) idx, h1,h2,p1,p2, ckl
|
||||
call bielec_integrals_index(h1,h2,p2,p1,idx)
|
||||
write(iunit_two_body_dm_aa,f) idx, h1,h2,p2,p1, -ckl
|
||||
enddo
|
||||
do i = 1, elec_beta_num
|
||||
p2=occ(i,2)
|
||||
h2=p2
|
||||
call bielec_integrals_index(h1,h2,p1,p2,idx)
|
||||
write(iunit_two_body_dm_ab,f) idx, h1,h2,p1,p2, ckl
|
||||
enddo
|
||||
case (2) ! Beta single excitation
|
||||
do i = 1, elec_alpha_num
|
||||
p2=occ(i,1)
|
||||
h2=p2
|
||||
call bielec_integrals_index(h1,h2,p1,p2,idx)
|
||||
write(iunit_two_body_dm_ab,f) idx, h1,h2,p1,p2, ckl
|
||||
enddo
|
||||
do i = 1, elec_beta_num
|
||||
p2=occ(i,2)
|
||||
h2=p2
|
||||
call bielec_integrals_index(h1,h2,p1,p2,idx)
|
||||
write(iunit_two_body_dm_bb,f) idx, h1,h2,p1,p2, ckl
|
||||
call bielec_integrals_index(h1,h2,p2,p1,idx)
|
||||
write(iunit_two_body_dm_bb,f) idx, h1,h2,p2,p1, -ckl
|
||||
enddo
|
||||
end select
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
! Sort file
|
||||
! Merge coefs
|
||||
|
||||
close(iunit_two_body_dm_aa)
|
||||
close(iunit_two_body_dm_ab)
|
||||
close(iunit_two_body_dm_bb)
|
||||
character*(128) :: filename
|
||||
filename = trim(ezfio_filename)//'/work/two_body_aa.tmp'
|
||||
call system('sort -n '//trim(filename)//' > '//trim(filename)//'2 ; cp '//trim(filename)//'2 '//trim(filename))
|
||||
filename = trim(ezfio_filename)//'/work/two_body_ab.tmp'
|
||||
call system('sort -n '//trim(filename)//' > '//trim(filename)//'2 ; cp '//trim(filename)//'2 '//trim(filename))
|
||||
filename = trim(ezfio_filename)//'/work/two_body_bb.tmp'
|
||||
call system('sort -n '//trim(filename)//' > '//trim(filename)//'2 ; cp '//trim(filename)//'2 '//trim(filename))
|
||||
iunit_two_body_dm_aa = getUnitAndOpen(trim(ezfio_filename)//'/work/two_body_aa.tmp','r')
|
||||
iunit_two_body_dm_ab = getUnitAndOpen(trim(ezfio_filename)//'/work/two_body_ab.tmp','r')
|
||||
iunit_two_body_dm_bb = getUnitAndOpen(trim(ezfio_filename)//'/work/two_body_bb.tmp','r')
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
BEGIN_TEMPLATE
|
||||
|
||||
BEGIN_PROVIDER [ integer, size_two_body_dm_$AA ]
|
||||
implicit none
|
||||
use bitmasks
|
||||
BEGIN_DOC
|
||||
! Size of the two body $ALPHA density matrix
|
||||
END_DOC
|
||||
integer *8 :: key, key_old
|
||||
rewind(iunit_two_body_dm_$AA)
|
||||
size_two_body_dm_$AA = 0
|
||||
key = 0_8
|
||||
key_old = key
|
||||
do while (.True.)
|
||||
read(iunit_two_body_dm_$AA,*,END=99) key
|
||||
if (key /= key_old) then
|
||||
size_two_body_dm_$AA += 1
|
||||
key_old = key
|
||||
endif
|
||||
end do
|
||||
99 continue
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ integer, two_body_dm_index_$AA, (4,size_two_body_dm_$AA) ]
|
||||
&BEGIN_PROVIDER [ double precision, two_body_dm_value_$AA, (size_two_body_dm_$AA) ]
|
||||
implicit none
|
||||
use bitmasks
|
||||
BEGIN_DOC
|
||||
! Two body $ALPHA density matrix
|
||||
END_DOC
|
||||
rewind(iunit_two_body_dm_$AA)
|
||||
integer *8 :: key, key_old
|
||||
integer :: ii, i,j,k,l
|
||||
double precision :: c
|
||||
key = 0_8
|
||||
key_old = key
|
||||
ii = 0
|
||||
do while (.True.)
|
||||
read(iunit_two_body_dm_$AA,*,END=99) key, i,j,k,l, c
|
||||
if (key /= key_old) then
|
||||
ii += 1
|
||||
two_body_dm_index_$AA(1,ii) = i
|
||||
two_body_dm_index_$AA(2,ii) = j
|
||||
two_body_dm_index_$AA(3,ii) = k
|
||||
two_body_dm_index_$AA(4,ii) = l
|
||||
two_body_dm_value_$AA(ii) = 0.d0
|
||||
key_old = key
|
||||
endif
|
||||
two_body_dm_value_$AA(ii) += c
|
||||
enddo
|
||||
99 continue
|
||||
close(iunit_two_body_dm_$AA, status='DELETE')
|
||||
END_PROVIDER
|
||||
|
||||
SUBST [ AA, ALPHA ]
|
||||
|
||||
aa ; alpha-alpha ;;
|
||||
ab ; alpha-beta ;;
|
||||
bb ; beta-beta ;;
|
||||
|
||||
END_TEMPLATE
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ double precision, two_body_dm_diag_aa, (mo_tot_num_align,mo_tot_num)]
|
||||
&BEGIN_PROVIDER [ double precision, two_body_dm_diag_bb, (mo_tot_num_align,mo_tot_num)]
|
||||
&BEGIN_PROVIDER [ double precision, two_body_dm_diag_ab, (mo_tot_num_align,mo_tot_num)]
|
||||
implicit none
|
||||
use bitmasks
|
||||
BEGIN_DOC
|
||||
! diagonal part of the two body density matrix
|
||||
END_DOC
|
||||
integer :: i,j,k,e1,e2
|
||||
integer :: occ(N_int*bit_kind_size,2)
|
||||
double precision :: ck
|
||||
integer :: n_occ_alpha
|
||||
two_body_dm_diag_aa=0.d0
|
||||
two_body_dm_diag_ab=0.d0
|
||||
two_body_dm_diag_bb=0.d0
|
||||
do k = 1, det_num
|
||||
call bitstring_to_list(det_provider(1,1,k), occ(1,1), n_occ_alpha, N_int)
|
||||
call bitstring_to_list(det_provider(1,2,k), occ(1,2), n_occ_alpha, N_int)
|
||||
ck = det_coef_provider(k) * det_coef_provider(k)
|
||||
do i = 1,elec_alpha_num
|
||||
e1=occ(i,1)
|
||||
do j = 1,elec_alpha_num
|
||||
e2=occ(j,1)
|
||||
! alpha-alpha
|
||||
two_body_dm_diag_aa(e1,e2) = two_body_dm_diag_aa(e1,e2) + ck
|
||||
enddo
|
||||
do j = 1,elec_beta_num
|
||||
e2=occ(j,2)
|
||||
! alpha-beta
|
||||
two_body_dm_diag_ab(e1,e2) = two_body_dm_diag_ab(e1,e2) + ck
|
||||
enddo
|
||||
enddo
|
||||
do i = 1,elec_beta_num
|
||||
e1=occ(i,2)
|
||||
do j = 1,elec_beta_num
|
||||
e2=occ(j,2)
|
||||
! beta-beta
|
||||
two_body_dm_diag_bb(e1,e2) = two_body_dm_diag_bb(e1,e2) + ck
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ double precision, one_body_dm_a, (mo_tot_num_align,mo_tot_num) ]
|
||||
&BEGIN_PROVIDER [ double precision, one_body_dm_b, (mo_tot_num_align,mo_tot_num) ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Alpha and beta one-body density matrix
|
||||
END_DOC
|
||||
|
||||
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
|
||||
one_body_dm_a = 0.d0
|
||||
one_body_dm_b = 0.d0
|
||||
|
||||
do k=1,det_num
|
||||
call bitstring_to_list(det_provider(1,1,k), occ(1,1), n_occ_alpha, N_int)
|
||||
call bitstring_to_list(det_provider(1,2,k), occ(1,2), n_occ_alpha, N_int)
|
||||
ck = det_coef_provider(k)
|
||||
do l=1,elec_alpha_num
|
||||
j = occ(l,1)
|
||||
one_body_dm_a(j,j) += ck*ck
|
||||
enddo
|
||||
do l=1,elec_beta_num
|
||||
j = occ(l,2)
|
||||
one_body_dm_b(j,j) += ck*ck
|
||||
enddo
|
||||
do l=1,k-1
|
||||
call get_excitation_degree(det_provider(1,1,k),det_provider(1,1,l),degree,N_int)
|
||||
if (degree /= 1) then
|
||||
cycle
|
||||
endif
|
||||
call get_mono_excitation(det_provider(1,1,k),det_provider(1,1,l),exc,phase,N_int)
|
||||
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
|
||||
ckl = ck * det_coef_provider(l) * phase
|
||||
if (s1==1) then
|
||||
one_body_dm_a(h1,p1) += ckl
|
||||
one_body_dm_a(p1,h1) += ckl
|
||||
else
|
||||
one_body_dm_b(h1,p1) += ckl
|
||||
one_body_dm_b(p1,h1) += ckl
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
END_PROVIDER
|
||||
|
||||
56
src/DensityMatrix/det_num.irp.f
Normal file
56
src/DensityMatrix/det_num.irp.f
Normal file
@ -0,0 +1,56 @@
|
||||
use bitmasks
|
||||
|
||||
BEGIN_PROVIDER [integer, det_num]
|
||||
det_num = 10
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ integer(bit_kind), det_provider, (N_int,2,det_num)]
|
||||
&BEGIN_PROVIDER [ double precision , det_coef_provider, (det_num) ]
|
||||
use bitmasks
|
||||
implicit none
|
||||
integer :: i
|
||||
det_provider = 0
|
||||
det_provider(1,1,1 ) = #001f ! 0000 0000 0001 1111
|
||||
det_provider(1,1,2 ) = #003b ! 0000 0000 0011 1011
|
||||
det_provider(1,1,3 ) = #008f ! 0000 0000 1000 1111
|
||||
det_provider(1,1,4 ) = #0057 ! 0000 0000 0101 0111
|
||||
det_provider(1,1,5 ) = #100f ! 0001 0000 0000 1111
|
||||
det_provider(1,1,6 ) = #001f ! 0000 0000 0001 1111
|
||||
det_provider(1,1,7 ) = #003b ! 0000 0000 0011 1011
|
||||
det_provider(1,1,8 ) = #00c7 ! 0000 0000 1100 0111
|
||||
det_provider(1,1,9 ) = #00ab ! 0000 0000 1010 1011
|
||||
det_provider(1,1,10) = #0073 ! 0000 0000 0111 0011
|
||||
det_provider(1,2,1 ) = #0007 ! 0000 0000 0001 0111
|
||||
det_provider(1,2,2 ) = #0023 ! 0000 0000 0010 0011
|
||||
det_provider(1,2,3 ) = #0023 ! 0000 0000 0010 0011
|
||||
det_provider(1,2,4 ) = #0023 ! 0000 0000 0010 0011
|
||||
det_provider(1,2,5 ) = #0015 ! 0000 0000 0001 0101
|
||||
det_provider(1,2,6 ) = #000d ! 0000 0000 0000 1101
|
||||
det_provider(1,2,7 ) = #0007 ! 0000 0000 0000 0111
|
||||
det_provider(1,2,8 ) = #0007 ! 0000 0000 0000 0111
|
||||
det_provider(1,2,9 ) = #0007 ! 0000 0000 0000 0111
|
||||
det_provider(1,2,10) = #0007 ! 0000 0000 0000 0111
|
||||
det_coef_provider = (/ &
|
||||
0.993536117982429D+00, &
|
||||
-0.556089064313864D-01, &
|
||||
0.403074722590178D-01, &
|
||||
0.403074717461626D-01, &
|
||||
-0.340290975461932D-01, &
|
||||
-0.340290958781670D-01, &
|
||||
-0.333949939765448D-01, &
|
||||
0.333418373363987D-01, &
|
||||
-0.316337211787351D-01, &
|
||||
-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 *, ''
|
||||
|
||||
|
||||
END_PROVIDER
|
||||
@ -120,6 +120,10 @@ subroutine copy_H_apply_buffer_to_wf
|
||||
N_det = N_det + H_apply_buffer_N_det
|
||||
TOUCH N_det
|
||||
|
||||
if (psi_det_size < N_det) then
|
||||
psi_det_size = N_det
|
||||
TOUCH psi_det_size
|
||||
endif
|
||||
do i=1,N_det_old
|
||||
do k=1,N_int
|
||||
psi_det(k,1,i) = buffer_det(k,1,i)
|
||||
|
||||
@ -2,6 +2,11 @@ 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
|
||||
$declarations
|
||||
integer(omp_lock_kind) :: lck
|
||||
integer(bit_kind),intent(in) :: key_in(N_int,2)
|
||||
@ -24,9 +29,9 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2 $parame
|
||||
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
|
||||
PROVIDE mo_integrals_map ref_bitmask_energy
|
||||
PROVIDE mo_bielec_integrals_in_map
|
||||
|
||||
$set_i_H_j_threshold
|
||||
@ -34,8 +39,6 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2 $parame
|
||||
$omp_init_lock
|
||||
|
||||
|
||||
E_ref = diag_H_mat_elem(key_in,N_int)
|
||||
|
||||
$initialization
|
||||
|
||||
$omp_parallel
|
||||
@ -233,6 +236,11 @@ 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
|
||||
$declarations
|
||||
integer(omp_lock_kind) :: lck
|
||||
integer(bit_kind),intent(in) :: key_in(N_int,2)
|
||||
@ -255,9 +263,9 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1 $parameters )
|
||||
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
|
||||
PROVIDE mo_integrals_map ref_bitmask_energy
|
||||
PROVIDE mo_bielec_integrals_in_map
|
||||
|
||||
$set_i_H_j_threshold
|
||||
@ -265,8 +273,6 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1 $parameters )
|
||||
$omp_init_lock
|
||||
|
||||
|
||||
E_ref = diag_H_mat_elem(key_in,N_int)
|
||||
|
||||
$initialization
|
||||
|
||||
$omp_parallel
|
||||
|
||||
@ -51,7 +51,8 @@ Documentation
|
||||
.. NEEDED_MODULES file.
|
||||
|
||||
`copy_h_apply_buffer_to_wf <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/H_apply.irp.f#L93>`_
|
||||
None
|
||||
Undocumented
|
||||
|
||||
`h_apply_buffer_coef <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/H_apply.irp.f#L82>`_
|
||||
Buffer of determinants/coefficients for H_apply. Uninitialized. Filled by H_apply subroutines.
|
||||
|
||||
@ -68,23 +69,49 @@ None
|
||||
Theshold on | <Di|H|Dj> |
|
||||
|
||||
`resize_h_apply_buffer_det <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/H_apply.irp.f#L31>`_
|
||||
None
|
||||
Undocumented
|
||||
|
||||
`davidson_diag <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/davidson.irp.f#L18>`_
|
||||
Davidson diagonalization.
|
||||
.br
|
||||
dets_in : bitmasks corresponding to determinants
|
||||
.br
|
||||
u_in : guess coefficients on the various states. Overwritten
|
||||
on exit
|
||||
.br
|
||||
dim_in : leftmost dimension of u_in
|
||||
.br
|
||||
sze : Number of determinants
|
||||
.br
|
||||
N_st : Number of eigenstates
|
||||
.br
|
||||
Initial guess vectors are not necessarily orthonormal
|
||||
|
||||
`davidson_iter_max <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/davidson.irp.f#L1>`_
|
||||
Max number of Davidson iterations
|
||||
|
||||
`davidson_sze_max <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/davidson.irp.f#L9>`_
|
||||
Max number of Davidson sizes
|
||||
|
||||
`n_det <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L11>`_
|
||||
Number of determinants in the wave function
|
||||
|
||||
`n_det_generators <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L47>`_
|
||||
`n_det_generators <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L55>`_
|
||||
Number of generator determinants in the wave function
|
||||
|
||||
`n_states <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L3>`_
|
||||
Number of states to consider
|
||||
|
||||
`psi_coef <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L20>`_
|
||||
`psi_coef <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L28>`_
|
||||
The wave function. Initialized with Hartree-Fock
|
||||
|
||||
`psi_det <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L19>`_
|
||||
`psi_det <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L27>`_
|
||||
The wave function. Initialized with Hartree-Fock
|
||||
|
||||
`psi_generators <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L55>`_
|
||||
`psi_det_size <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L19>`_
|
||||
Size of the psi_det/psi_coef arrays
|
||||
|
||||
`psi_generators <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L63>`_
|
||||
Determinants on which H is applied
|
||||
|
||||
`double_exc_bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants_bitmasks.irp.f#L40>`_
|
||||
@ -108,10 +135,10 @@ None
|
||||
`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#L842>`_
|
||||
`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#L887>`_
|
||||
`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>`_
|
||||
@ -121,15 +148,10 @@ None
|
||||
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#L779>`_
|
||||
`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#L602>`_
|
||||
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>`_
|
||||
None
|
||||
`get_double_excitation <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L140>`_
|
||||
`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
|
||||
|
||||
`get_excitation <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L30>`_
|
||||
@ -138,20 +160,28 @@ None
|
||||
`get_excitation_degree <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L1>`_
|
||||
Returns the excitation degree between two determinants
|
||||
|
||||
`get_excitation_degree_vector <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L518>`_
|
||||
`get_excitation_degree_vector <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L520>`_
|
||||
Applies get_excitation_degree to an array of determinants
|
||||
|
||||
`get_mono_excitation <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L273>`_
|
||||
`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#L935>`_
|
||||
`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
|
||||
|
||||
`i_h_j <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L354>`_
|
||||
`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
|
||||
.br
|
||||
H_jj : array of <j|H|j>
|
||||
|
||||
`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#L490>`_
|
||||
None
|
||||
`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>`_
|
||||
H matrix on the basis of the slater deter;inants defined by psi_det
|
||||
|
||||
|
||||
254
src/Dets/connected_to_ref.irp.f
Normal file
254
src/Dets/connected_to_ref.irp.f
Normal file
@ -0,0 +1,254 @@
|
||||
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(ishft(Nint,-1),2,Ndet)
|
||||
integer(bit_kind), intent(in) :: key(ishft(Nint,-1),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 == 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)
|
||||
implicit none
|
||||
! If true, det is not in ref
|
||||
! If false, det may be in ref
|
||||
|
||||
integer, intent(in) :: key(Nint,2), Nint
|
||||
integer :: key_int
|
||||
integer*1 :: key_short(4)
|
||||
!DIR$ ATTRIBUTES ALIGN : 32 :: key_short
|
||||
equivalence (key_int,key_short)
|
||||
|
||||
integer :: i, ispin
|
||||
|
||||
det_is_not_or_may_be_in_ref = .False.
|
||||
do ispin=1,2
|
||||
do i=1,Nint
|
||||
key_int = key(i,ispin)
|
||||
if ( &
|
||||
key_pattern_not_in_ref(key_short(1), i, ispin) &
|
||||
.or.key_pattern_not_in_ref(key_short(2), i, ispin) &
|
||||
.or.key_pattern_not_in_ref(key_short(3), i, ispin) &
|
||||
.or.key_pattern_not_in_ref(key_short(4), i, ispin) &
|
||||
) then
|
||||
det_is_not_or_may_be_in_ref = .True.
|
||||
return
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
|
||||
end
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ logical, key_pattern_not_in_ref, (-128:127,N_int,2) ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Min and max values of the integers of the keys of the reference
|
||||
END_DOC
|
||||
|
||||
integer :: i, j, ispin
|
||||
integer :: key
|
||||
integer*1 :: key_short(4)
|
||||
equivalence (key,key_short)
|
||||
integer :: idx
|
||||
|
||||
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)
|
||||
key_pattern_not_in_ref( key_short(1), i, ispin ) = .False.
|
||||
key_pattern_not_in_ref( key_short(2), i, ispin ) = .False.
|
||||
key_pattern_not_in_ref( key_short(3), i, ispin ) = .False.
|
||||
key_pattern_not_in_ref( key_short(4), i, ispin ) = .False.
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
277
src/Dets/davidson.irp.f
Normal file
277
src/Dets/davidson.irp.f
Normal file
@ -0,0 +1,277 @@
|
||||
BEGIN_PROVIDER [ integer, davidson_iter_max]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Max number of Davidson iterations
|
||||
END_DOC
|
||||
davidson_iter_max = 100
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ integer, davidson_sze_max]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Max number of Davidson sizes
|
||||
END_DOC
|
||||
ASSERT (davidson_sze_max <= davidson_iter_max)
|
||||
davidson_sze_max = 8
|
||||
END_PROVIDER
|
||||
|
||||
subroutine davidson_diag(dets_in,u_in,energies,dim_in,sze,N_st,Nint)
|
||||
use bitmasks
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Davidson diagonalization.
|
||||
!
|
||||
! 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)
|
||||
|
||||
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(:,:,:), H_jj(:), 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), &
|
||||
H_jj(sze), &
|
||||
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 H_jj,Nint,dets_in,u_in) &
|
||||
!$OMP PRIVATE(k,l,kl,i)
|
||||
|
||||
!$OMP DO
|
||||
do i=1,sze
|
||||
H_jj(i) = diag_h_mat_elem(dets_in(1,1,i),Nint)
|
||||
enddo
|
||||
!$OMP END DO NOWAIT
|
||||
|
||||
! 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
|
||||
|
||||
! print *, '***************'
|
||||
! do i=1,iter
|
||||
! do k=1,N_st
|
||||
! do j=1,iter
|
||||
! do l=1,N_st
|
||||
! print '(4(I4,X),F16.8)', i,j,k,l, u_dot_v(U(1,k,i),U(1,l,j),sze)
|
||||
! enddo
|
||||
! enddo
|
||||
! enddo
|
||||
! enddo
|
||||
! print *, '***************'
|
||||
|
||||
! 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, &
|
||||
H_jj, &
|
||||
W, &
|
||||
U, &
|
||||
R, &
|
||||
h, &
|
||||
y, &
|
||||
lambda &
|
||||
)
|
||||
end
|
||||
|
||||
@ -13,11 +13,19 @@ BEGIN_PROVIDER [ integer, N_det ]
|
||||
BEGIN_DOC
|
||||
! Number of determinants in the wave function
|
||||
END_DOC
|
||||
N_det = max(1,N_states)
|
||||
N_det = 1
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ integer(bit_kind), psi_det, (N_int,2,N_det) ]
|
||||
&BEGIN_PROVIDER [ double precision, psi_coef, (N_det,N_states) ]
|
||||
BEGIN_PROVIDER [ integer, psi_det_size ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Size of the psi_det/psi_coef arrays
|
||||
END_DOC
|
||||
psi_det_size = 1000
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ integer(bit_kind), psi_det, (N_int,2,psi_det_size) ]
|
||||
&BEGIN_PROVIDER [ double precision, psi_coef, (psi_det_size,N_states) ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! The wave function. Initialized with Hartree-Fock
|
||||
@ -52,7 +60,7 @@ BEGIN_PROVIDER [ integer, N_det_generators ]
|
||||
N_det_generators = N_det
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ integer(bit_kind), psi_generators, (N_int,2,N_det) ]
|
||||
BEGIN_PROVIDER [ integer(bit_kind), psi_generators, (N_int,2,psi_det_size) ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Determinants on which H is applied
|
||||
|
||||
174
src/Dets/filter_connected.irp.f
Normal file
174
src/Dets/filter_connected.irp.f
Normal file
@ -0,0 +1,174 @@
|
||||
|
||||
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 (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
|
||||
|
||||
@ -74,6 +74,7 @@ subroutine get_excitation(det1,det2,exc,degree,phase,Nint)
|
||||
end
|
||||
|
||||
subroutine decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
|
||||
use bitmasks
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Decodes the exc arrays returned by get_excitation.
|
||||
@ -487,7 +488,8 @@ 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)
|
||||
@ -515,14 +517,14 @@ end
|
||||
|
||||
|
||||
|
||||
subroutine get_excitation_degree_vector(key1,key2,degree,Nint,sze,sze_max,idx)
|
||||
subroutine get_excitation_degree_vector(key1,key2,degree,Nint,sze,idx)
|
||||
use bitmasks
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Applies get_excitation_degree to an array of determinants
|
||||
END_DOC
|
||||
integer, intent(in) :: Nint, sze,sze_max
|
||||
integer(bit_kind), intent(in) :: key1(Nint,2,sze_max)
|
||||
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) :: degree(sze)
|
||||
integer, intent(out) :: idx(0:sze)
|
||||
@ -531,7 +533,6 @@ subroutine get_excitation_degree_vector(key1,key2,degree,Nint,sze,sze_max,idx)
|
||||
|
||||
ASSERT (Nint > 0)
|
||||
ASSERT (sze > 0)
|
||||
ASSERT (sze_max >= sze)
|
||||
|
||||
l=1
|
||||
if (Nint==1) then
|
||||
@ -599,185 +600,8 @@ end
|
||||
|
||||
|
||||
|
||||
subroutine filter_connected(key1,key2,Nint,sze,sze_max,idx)
|
||||
use bitmasks
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Filters out the determinants that are not connected by H
|
||||
END_DOC
|
||||
integer, intent(in) :: Nint, sze,sze_max
|
||||
integer(bit_kind), intent(in) :: key1(Nint,2,sze_max)
|
||||
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)
|
||||
ASSERT (sze_max >= sze)
|
||||
|
||||
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,sze_max,idx)
|
||||
use bitmasks
|
||||
implicit none
|
||||
integer, intent(in) :: Nint, sze,sze_max
|
||||
integer(bit_kind), intent(in) :: key1(Nint,2,sze_max)
|
||||
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)
|
||||
ASSERT (sze_max >= sze)
|
||||
|
||||
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)
|
||||
use bitmasks
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Computes <i|H|i>
|
||||
@ -947,3 +771,61 @@ subroutine get_occ_from_key(key,occ,Nint)
|
||||
call bitstring_to_list(key(1,2), occ(1,2), tmp, Nint)
|
||||
|
||||
end
|
||||
|
||||
subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint)
|
||||
use bitmasks
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Computes v_0 = H|u_0>
|
||||
!
|
||||
! n : number of determinants
|
||||
!
|
||||
! H_jj : array of <j|H|j>
|
||||
END_DOC
|
||||
integer, intent(in) :: n,Nint
|
||||
double precision, intent(out) :: v_0(n)
|
||||
double precision, intent(in) :: u_0(n)
|
||||
double precision, intent(in) :: H_jj(n)
|
||||
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)
|
||||
ASSERT (Nint == N_int)
|
||||
ASSERT (n>0)
|
||||
PROVIDE ref_bitmask_energy
|
||||
integer, parameter :: block_size = 157
|
||||
!$OMP PARALLEL DEFAULT(NONE) &
|
||||
!$OMP PRIVATE(i,hij,j,k,idx,jj,vt) SHARED(n,H_jj,u_0,keys_tmp,Nint)&
|
||||
!$OMP SHARED(v_0)
|
||||
!$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 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)
|
||||
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
|
||||
!$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
|
||||
|
||||
|
||||
|
||||
|
||||
@ -82,7 +82,8 @@ Documentation
|
||||
Diagonal Fock matrix in the MO basis
|
||||
|
||||
`scf_iteration <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/mo_SCF_iterations.irp.f#L1>`_
|
||||
None
|
||||
Undocumented
|
||||
|
||||
`do_diis <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/options.irp.f#L41>`_
|
||||
If True, compute integrals on the fly
|
||||
|
||||
|
||||
@ -54,8 +54,10 @@ Documentation
|
||||
Aligned variable for dimensioning of arrays
|
||||
|
||||
`mo_as_eigvectors_of_mo_matrix <http://github.com/LCPQ/quantum_package/tree/master/src/MOs/utils.irp.f#L21>`_
|
||||
None
|
||||
Undocumented
|
||||
|
||||
`save_mos <http://github.com/LCPQ/quantum_package/tree/master/src/MOs/utils.irp.f#L1>`_
|
||||
None
|
||||
Undocumented
|
||||
|
||||
|
||||
|
||||
|
||||
@ -85,11 +85,6 @@ $(info -----------------------------------------------)
|
||||
endif
|
||||
|
||||
|
||||
# Update the README.rst file for GitHub, and git add it
|
||||
|
||||
$(shell update_README.py)
|
||||
|
||||
|
||||
# Define the Makefile common variables and rules
|
||||
|
||||
EZFIO_DIR=$(QPACKAGE_ROOT)/EZFIO
|
||||
@ -116,7 +111,8 @@ LIB+=$(EZFIO) $(MKL)
|
||||
IRPF90+=$(patsubst %, -I %, $(INCLUDE_DIRS)) $(IRPF90_FLAGS)
|
||||
|
||||
irpf90.make: $(filter-out IRPF90_temp/%, $(wildcard */*.irp.f)) $(wildcard *.irp.f) $(wildcard *.inc.f) Makefile $(EZFIO) NEEDED_MODULES $(wildcard *.py)
|
||||
$(IRPF90)
|
||||
- $(IRPF90)
|
||||
- update_README.py
|
||||
|
||||
Makefile.depend: Makefile
|
||||
$(QPACKAGE_ROOT)/scripts/create_Makefile_depend.sh
|
||||
|
||||
@ -11,6 +11,7 @@
|
||||
double precision :: A_center(3), B_center(3)
|
||||
integer :: power_A(3), power_B(3)
|
||||
double precision :: d_a_2,d_2
|
||||
PROVIDE all_utils
|
||||
dim1=100
|
||||
BEGIN_DOC
|
||||
! second derivatives matrix elements in the ao basis
|
||||
|
||||
@ -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
|
||||
|
||||
105
src/Perturbation/README.rst
Normal file
105
src/Perturbation/README.rst
Normal file
@ -0,0 +1,105 @@
|
||||
===================
|
||||
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.
|
||||
|
||||
|
||||
|
||||
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>`_
|
||||
|
||||
13
src/Perturbation/perturbation.irp.f
Normal file
13
src/Perturbation/perturbation.irp.f
Normal 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
|
||||
@ -1,10 +1,24 @@
|
||||
subroutine ortho_lowdin(overlap,lda,n,C,ldc,m)
|
||||
subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Compute U.S^-1/2 canonical orthogonalization
|
||||
! Compute C_new=C_old.S^-1/2 canonical orthogonalization.
|
||||
!
|
||||
! overlap : overlap matrix
|
||||
!
|
||||
! LDA : leftmost dimension of overlap array
|
||||
!
|
||||
! N : Overlap matrix is NxN (array is (LDA,N) )
|
||||
!
|
||||
! C : Coefficients of the vectors to orthogonalize. On exit,
|
||||
! orthogonal vectors
|
||||
!
|
||||
! LDC : leftmost dimension of C
|
||||
!
|
||||
! m : Coefficients matrix is MxN, ( array is (LDC,N) )
|
||||
!
|
||||
END_DOC
|
||||
|
||||
integer, intent(in) :: lda, ldc, n, m
|
||||
integer, intent(in) :: LDA, ldc, n, m
|
||||
double precision, intent(in) :: overlap(lda,n)
|
||||
double precision, intent(inout) :: C(ldc,n)
|
||||
double precision :: U(ldc,n)
|
||||
@ -34,37 +48,45 @@ subroutine ortho_lowdin(overlap,lda,n,C,ldc,m)
|
||||
stop
|
||||
endif
|
||||
|
||||
|
||||
!$OMP PARALLEL DEFAULT(NONE) &
|
||||
!$OMP SHARED(S_half,U,D,Vt,n,C,m) &
|
||||
!$OMP PRIVATE(i,j,k)
|
||||
|
||||
!$OMP DO
|
||||
do i=1,n
|
||||
if ( D(i) < 1.d-6 ) then
|
||||
D(i) = 0.d0
|
||||
else
|
||||
D(i) = 1.d0/dsqrt(D(i))
|
||||
endif
|
||||
enddo
|
||||
|
||||
S_half = 0.d0
|
||||
do k=1,n
|
||||
do j=1,n
|
||||
do i=1,n
|
||||
S_half(i,j) += U(i,k)*D(k)*Vt(k,j)
|
||||
enddo
|
||||
S_half(j,i) = 0.d0
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
|
||||
do k=1,n
|
||||
!$OMP DO
|
||||
do j=1,n
|
||||
do i=1,n
|
||||
S_half(i,j) = S_half(i,j) + U(i,k)*D(k)*Vt(k,j)
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO NOWAIT
|
||||
enddo
|
||||
|
||||
!$OMP DO
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
U(i,j) = C(i,j)
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
|
||||
C = 0.d0
|
||||
do j=1,n
|
||||
do i=1,m
|
||||
do k=1,n
|
||||
C(i,j) += U(i,k)*S_half(k,j)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END PARALLEL
|
||||
|
||||
call dgemm('N','N',m,n,n,1.d0,U,size(U,1),S_half,size(S_half,1),0.d0,C,size(C,1))
|
||||
|
||||
end
|
||||
|
||||
@ -171,6 +193,17 @@ subroutine lapack_diag(eigvalues,eigvectors,H,nmax,n)
|
||||
allocate(A(nmax,n),eigenvalues(nmax),work(4*nmax))
|
||||
integer :: LWORK, info, i,j,l,k
|
||||
A=H
|
||||
|
||||
! if (n<30) then
|
||||
! do i=1,n
|
||||
! do j=1,n
|
||||
! print *, j,i, H(j,i)
|
||||
! enddo
|
||||
! print *, '---'
|
||||
! enddo
|
||||
! print *, '---'
|
||||
! endif
|
||||
|
||||
LWORK = 4*nmax
|
||||
call dsyev( 'V', 'U', n, A, nmax, eigenvalues, work, LWORK, info )
|
||||
if (info < 0) then
|
||||
|
||||
@ -482,7 +482,7 @@ subroutine cache_map_get_interval(map, key, value, ibegin, iend, idx)
|
||||
if (idx > 0) then
|
||||
value = map%value(idx)
|
||||
else
|
||||
value = 0.
|
||||
value = 0._integral_kind
|
||||
endif
|
||||
end
|
||||
|
||||
|
||||
@ -214,7 +214,6 @@ double precision function u_dot_v(u,v,sze)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Compute <u|v>
|
||||
! u and v are expected to be aligned in memory.
|
||||
END_DOC
|
||||
integer, intent(in) :: sze
|
||||
double precision, intent(in) :: u(sze),v(sze)
|
||||
@ -227,14 +226,10 @@ double precision function u_dot_v(u,v,sze)
|
||||
t3 = t2+t2
|
||||
t4 = t3+t2
|
||||
u_dot_v = 0.d0
|
||||
!DIR$ VECTOR ALWAYS
|
||||
!DIR$ VECTOR ALIGNED
|
||||
do i=1,t2
|
||||
u_dot_v = u_dot_v + u(t1+i)*v(t1+i) + u(t2+i)*v(t2+i) + &
|
||||
u(t3+i)*v(t3+i) + u(t4+i)*v(t4+i)
|
||||
enddo
|
||||
!DIR$ VECTOR ALWAYS
|
||||
!DIR$ VECTOR ALIGNED
|
||||
do i=t4+t2+1,sze
|
||||
u_dot_v = u_dot_v + u(i)*v(i)
|
||||
enddo
|
||||
@ -245,7 +240,6 @@ double precision function u_dot_u(u,sze)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Compute <u|u>
|
||||
! u is expected to be aligned in memory.
|
||||
END_DOC
|
||||
integer, intent(in) :: sze
|
||||
double precision, intent(in) :: u(sze)
|
||||
@ -259,12 +253,16 @@ double precision function u_dot_u(u,sze)
|
||||
t3 = t2+t2
|
||||
t4 = t3+t2
|
||||
u_dot_u = 0.d0
|
||||
do i=1,t2
|
||||
u_dot_u = u_dot_u + u(t1+i)*u(t1+i) + u(t2+i)*u(t2+i) + &
|
||||
u(t3+i)*u(t3+i) + u(t4+i)*u(t4+i)
|
||||
enddo
|
||||
do i=t4+t2+1,sze
|
||||
u_dot_u = u_dot_u+u(i)*u(i)
|
||||
! do i=1,t2
|
||||
! u_dot_u = u_dot_u + u(t1+i)*u(t1+i) + u(t2+i)*u(t2+i) + &
|
||||
! u(t3+i)*u(t3+i) + u(t4+i)*u(t4+i)
|
||||
! enddo
|
||||
! do i=t4+t2+1,sze
|
||||
! u_dot_u = u_dot_u+u(i)*u(i)
|
||||
! enddo
|
||||
|
||||
do i=1,sze
|
||||
u_dot_u = u_dot_u + u(i)*u(i)
|
||||
enddo
|
||||
|
||||
end
|
||||
|
||||
Loading…
Reference in New Issue
Block a user