mirror of
https://github.com/LCPQ/quantum_package
synced 2024-11-12 17:13:54 +01:00
Merge lpqlx139:/home/scemama/quantum_package
This commit is contained in:
commit
098a386411
@ -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
|
||||
|
||||
|
||||
@ -30,7 +30,5 @@ Documentation
|
||||
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#L1>`_
|
||||
None
|
||||
|
||||
|
||||
|
||||
@ -1,19 +1,23 @@
|
||||
program cisd
|
||||
implicit none
|
||||
integer :: i
|
||||
integer :: i,k
|
||||
double precision, allocatable :: eigvalues(:),eigvectors(:,:)
|
||||
PROVIDE ref_bitmask_energy
|
||||
call H_apply_cisd
|
||||
double precision, allocatable :: eigvalues(:),eigvectors(:,:)
|
||||
allocate(eigvalues(n_det),eigvectors(n_det,n_det))
|
||||
allocate(eigvalues(n_states),eigvectors(n_det,n_states))
|
||||
print *, 'N_det = ', N_det
|
||||
call lapack_diag(eigvalues,eigvectors,H_matrix_all_dets,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
|
||||
|
||||
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)
|
||||
|
||||
@ -26,7 +26,7 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2 $parame
|
||||
integer, allocatable :: ia_ja_pairs(:,:,:)
|
||||
double precision :: diag_H_mat_elem, E_ref
|
||||
|
||||
PROVIDE mo_integrals_map
|
||||
PROVIDE mo_integrals_map ref_bitmask_energy
|
||||
PROVIDE mo_bielec_integrals_in_map
|
||||
|
||||
$set_i_H_j_threshold
|
||||
@ -257,7 +257,7 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1 $parameters )
|
||||
integer, allocatable :: ia_ja_pairs(:,:,:)
|
||||
double precision :: diag_H_mat_elem, E_ref
|
||||
|
||||
PROVIDE mo_integrals_map
|
||||
PROVIDE mo_integrals_map ref_bitmask_energy
|
||||
PROVIDE mo_bielec_integrals_in_map
|
||||
|
||||
$set_i_H_j_threshold
|
||||
|
||||
@ -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#L840>`_
|
||||
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#L885>`_
|
||||
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,16 @@ 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#L778>`_
|
||||
Computes <i|H|i>
|
||||
|
||||
`filter_connected <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L602>`_
|
||||
`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>`_
|
||||
None
|
||||
`get_double_excitation <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L140>`_
|
||||
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
|
||||
|
||||
`get_excitation <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L30>`_
|
||||
@ -138,20 +166,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#L933>`_
|
||||
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#L949>`_
|
||||
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_psim <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
|
||||
|
||||
|
||||
281
src/Dets/davidson.irp.f
Normal file
281
src/Dets/davidson.irp.f
Normal file
@ -0,0 +1,281 @@
|
||||
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 *, 'iter = ',iter
|
||||
|
||||
! 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)
|
||||
|
||||
print *, lambda(1:4)
|
||||
! 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 *, 'Lambda'
|
||||
print *, lambda(1:N_st) + nuclear_repulsion
|
||||
print *, 'Residual_norm'
|
||||
print *, residual_norm(1:N_st)
|
||||
print *, ''
|
||||
|
||||
converged = maxval(residual_norm) < 1.d-5
|
||||
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
|
||||
|
||||
@ -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.
|
||||
@ -488,6 +489,7 @@ end
|
||||
|
||||
|
||||
subroutine i_H_psim(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,14 +600,14 @@ end
|
||||
|
||||
|
||||
|
||||
subroutine filter_connected(key1,key2,Nint,sze,sze_max,idx)
|
||||
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,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) :: idx(0:sze)
|
||||
|
||||
@ -615,7 +616,6 @@ subroutine filter_connected(key1,key2,Nint,sze,sze_max,idx)
|
||||
|
||||
ASSERT (Nint > 0)
|
||||
ASSERT (sze > 0)
|
||||
ASSERT (sze_max >= sze)
|
||||
|
||||
l=1
|
||||
|
||||
@ -684,11 +684,11 @@ subroutine filter_connected(key1,key2,Nint,sze,sze_max,idx)
|
||||
idx(0) = l-1
|
||||
end
|
||||
|
||||
subroutine filter_connected_i_H_psi0(key1,key2,Nint,sze,sze_max,idx)
|
||||
subroutine filter_connected_i_H_psi0(key1,key2,Nint,sze,idx)
|
||||
use bitmasks
|
||||
implicit none
|
||||
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) :: idx(0:sze)
|
||||
|
||||
@ -697,7 +697,6 @@ subroutine filter_connected_i_H_psi0(key1,key2,Nint,sze,sze_max,idx)
|
||||
|
||||
ASSERT (Nint > 0)
|
||||
ASSERT (sze > 0)
|
||||
ASSERT (sze_max >= sze)
|
||||
|
||||
l=1
|
||||
|
||||
@ -777,7 +776,6 @@ end
|
||||
|
||||
|
||||
double precision function diag_H_mat_elem(det_in,Nint)
|
||||
use bitmasks
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Computes <i|H|i>
|
||||
@ -947,3 +945,58 @@ 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
|
||||
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) 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
|
||||
!$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 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
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO NOWAIT
|
||||
deallocate(idx)
|
||||
!$OMP END PARALLEL
|
||||
end
|
||||
|
||||
|
||||
@ -87,7 +87,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,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