LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2025-01-10 13:08:23 +01:00
Code Issues Releases Wiki Activity

Davidson multistates works in a very standard version

Browse Source
This commit is contained in:
Anthony Scemama 2014-05-14 15:40:40 +02:00
parent 59c37a79ff
commit e1d2a79b3b
11 changed files with 755 additions and 66 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
src/CISD/README.rst
View File

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

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

@ -1,19 +1,18 @@
program cisd
implicit none
integer :: i
double precision, allocatable :: eigvalues(:),eigvectors(:,:)
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)
psi_coef = psi_coef - 1.d-4
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
View 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
View 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

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

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

72
src/Dets/README.rst
View File

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

261
src/Dets/davidson.irp.f Normal file
View File

@ -0,0 +1,261 @@
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 = 10
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
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)
allocate( &
kl_pairs(2,N_st*(N_st+1)/2), &
H_jj(sze), &
W(sze,N_st), &
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)
! Conventions:
! i,j : 1,sze
! k,l : 1,N_st
!$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 NOWAIT
enddo
!$OMP END PARALLEL
! Orthonormalize initial guess
! ============================
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 DO DEFAULT(NONE) &
!$OMP SHARED(U,sze,N_st,overlap,kl_pairs,k_pairs) &
!$OMP PRIVATE(k,l,kl)
do kl=1,k_pairs
k = kl_pairs(1,kl)
l = kl_pairs(2,kl)
if (k==l) then
overlap(k,k) = u_dot_u(U(1,k,1),sze)
endif
overlap(k,l) = u_dot_v(U(1,k,1),U(1,l,1),sze)
overlap(l,k) = overlap(k,l)
enddo
!$OMP END PARALLEL DO
call ortho_lowdin(overlap,size(overlap,1),N_st,U,size(U,1),sze)
!$OMP PARALLEL DO DEFAULT(NONE) &
!$OMP PRIVATE(i) &
!$OMP SHARED(H_jj,Nint,dets_in,sze)
do i=1,sze
H_jj(i) = diag_h_mat_elem(dets_in(1,1,i),Nint)
enddo
!$OMP END PARALLEL DO
! Davidson iterations
! ===================
converged = .False.
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),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 iter2=1,iter-1
do k=1,N_st
h(k,iter2,l,iter) = u_dot_v(U(1,k,iter2),W(1,l),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),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
! --------------------------------------------------
!TODO dgemm
do k=1,N_st
do i=1,sze
U(i,k,iter+1) = 0.d0
do iter2=1,iter
do l=1,N_st
U(i,k,iter+1) += U(i,l,iter2)*y(l,iter2,k,1)
enddo
enddo
enddo
enddo
! Compute residual vector
! -----------------------
do k=1,N_st
call H_u_0(W(1,k),U(1,k,iter+1),H_jj,sze,dets_in,Nint)
enddo
do k=1,N_st
do i=1,sze
R(i,k) = lambda(k) * U(i,k,iter+1) - W(i,k)
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-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
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
deallocate ( &
kl_pairs, &
H_jj, &
W, &
U, &
R, &
h, &
y, &
lambda &
)
end

16
src/Dets/determinants.irp.f
View File

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

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

@ -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,46 @@ 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
ASSERT (Nint > 0)
ASSERT (Nint == N_int)
ASSERT (n>0)
!$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:n))
!$OMP DO SCHEDULE(dynamic)
do i=1,n
v_0(i) = H_jj(i) * u_0(i)
call filter_connected(keys_tmp,keys_tmp(1,1,i),Nint,n,idx)
do jj=1,idx(0)
j = idx(jj)
if (j/=i) 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)
endif
enddo
enddo
!$OMP END DO
deallocate(idx)
!$OMP END PARALLEL
end

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

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

6
src/Utils/util.irp.f
View File

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