LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-09-27 03:51:01 +02:00
Code Issues Releases Wiki Activity

Removed dead code

Browse Source
This commit is contained in:
Anthony Scemama 2018-12-18 15:45:40 +01:00
parent dea9dc9c36
commit 67e71e5931
21 changed files with 207 additions and 3320 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

9
TODO
View File

@ -1,7 +1,6 @@
* Virer tous les modules qui sont dans plugins
* Permettre aux utilisateurs de facilement deposer des plugins dans plugins via une commande
* Permettre de descendre plus bas dans l'arborescence de plugins pour permettre des `git clone` dans le
repertoire plugins
* Permettre de descendre plus bas dans l'arborescence de plugins pour permettre des `git clone` dans le repertoire plugins
* Mettre les fichiers de test dans le directory source
* Separer les integrales bielectroniques en AO et MO
* Creer une page web pas trop degueu et la mettre ici : http://lcpq.github.io/quantum_package
@ -39,5 +38,11 @@
* Programmers doc:
* Example : Simple Hartree-Fock program from scratch
* Examples : subroutine example_module
* Config file for Cray
* Noms des modules et fichiers en minuscule
* Implicit none obligatoire
* Doc des executables
* Il faut que le programme ait le meme nom que la subroutine

47
docs/source/users_guide/dft.rst Normal file
View File

@ -0,0 +1,47 @@
The DFT module uses Lebedev-Laikov grids, using the code distributed through CCL (http://www.ccl.net/).
.. code-block:: text
This subroutine is part of a set of subroutines that generate
Lebedev grids [1-6] for integration on a sphere. The original
C-code [1] was kindly provided by Dr. Dmitri N. Laikov and
translated into fortran by Dr. Christoph van Wuellen.
This subroutine was translated using a C to fortran77 conversion
tool written by Dr. Christoph van Wuellen.
Users of this code are asked to include reference [1] in their
publications, and in the user- and programmers-manuals
describing their codes.
This code was distributed through CCL (http://www.ccl.net/).
[1] V.I. Lebedev, and D.N. Laikov
"A quadrature formula for the sphere of the 131st
algebraic order of accuracy"
Doklady Mathematics, Vol. 59, No. 3, 1999, pp. 477-481.
[2] V.I. Lebedev
"A quadrature formula for the sphere of 59th algebraic
order of accuracy"
Russian Acad. Sci. Dokl. Math., Vol. 50, 1995, pp. 283-286.
[3] V.I. Lebedev, and A.L. Skorokhodov
"Quadrature formulas of orders 41, 47, and 53 for the sphere"
Russian Acad. Sci. Dokl. Math., Vol. 45, 1992, pp. 587-592.
[4] V.I. Lebedev
"Spherical quadrature formulas exact to orders 25-29"
Siberian Mathematical Journal, Vol. 18, 1977, pp. 99-107.
[5] V.I. Lebedev
"Quadratures on a sphere"
Computational Mathematics and Mathematical Physics, Vol. 16,
1976, pp. 10-24.
[6] V.I. Lebedev
"Values of the nodes and weights of ninth to seventeenth
order Gauss-Markov quadrature formulae invariant under the
octahedron group with inversion"
Computational Mathematics and Mathematical Physics, Vol. 15,
1975, pp. 44-51.

41
src/Davidson/find_reference.irp.f
View File

@ -1,41 +0,0 @@
subroutine find_reference(thresh,n_ref,result)
implicit none
double precision, intent(in) :: thresh
integer, intent(out) :: result(N_det),n_ref
integer :: i,j,istate
double precision :: i_H_psi_array(1), E0, hii, norm
double precision :: de
integer(bit_kind), allocatable :: psi_ref_(:,:,:)
double precision, allocatable :: psi_ref_coef_(:,:)
allocate(psi_ref_coef_(N_det,1), psi_ref_(N_int,2,N_det))
n_ref = 1
result(1) = 1
istate = 1
psi_ref_coef_(1,1) = psi_coef(1,istate)
psi_ref_(:,:,1) = psi_det(:,:,1)
norm = psi_ref_coef_(1,1) * psi_ref_coef_(1,1)
call u_0_H_u_0(E0,psi_ref_coef_,n_ref,psi_ref_,N_int,1,size(psi_ref_coef_,1))
print *, ''
print *, 'Reference determinants'
print *, '======================'
print *, ''
print *, n_ref, ': E0 = ', E0 + nuclear_repulsion
call debug_det(psi_ref_(1,1,n_ref),N_int)
do i=2,N_det
call i_h_psi(psi_det(1,1,i),psi_ref_(1,1,1),psi_ref_coef_(1,istate),N_int, &
n_ref,size(psi_ref_coef_,1),1,i_H_psi_array)
call i_H_j(psi_det(1,1,i),psi_det(1,1,i),N_int,hii)
de = i_H_psi_array(istate)**2 / (E0 - hii)
if (dabs(de) > thresh) then
n_ref += 1
result(n_ref) = i
psi_ref_(:,:,n_ref) = psi_det(:,:,i)
psi_ref_coef_(n_ref,1) = psi_coef(i,istate)
call u_0_H_u_0(E0,psi_ref_coef_,n_ref,psi_ref_,N_int,1,size(psi_ref_coef_,1))
print *, n_ref, ': E0 = ', E0 + nuclear_repulsion
call debug_det(psi_ref_(1,1,n_ref),N_int)
endif
enddo
end

12
src/Determinants/EZFIO.cfg
View File

@ -89,18 +89,6 @@ doc: Energy that should be obtained when truncating the wave function (optional)
type: Energy
default: 0.
[store_full_H_mat]
type: logical
doc: If |true|, the Davidson diagonalization is performed by storing the full |H| matrix up to n_det_max_stored. Be careful, it can cost a lot of memory but can also save a lot of CPU time
interface: ezfio,provider,ocaml
default: False
[n_det_max_stored]
type: Det_number_max
doc: Maximum number of determinants for which the full |H| matrix is stored. be careful, the memory requested scales as 10*n_det_max_stored**2. for instance, 90000 determinants represents a matrix of size 60 Gb.
interface: ezfio,provider,ocaml
default: 90000
[state_average_weight]
type: double precision
doc: Weight of the states in state-average calculations.

190
src/Determinants/create_excitations.irp.f
View File

@ -1,128 +1,90 @@
subroutine do_mono_excitation(key_in,i_hole,i_particle,ispin,i_ok)
implicit none
BEGIN_DOC
! Apply the mono excitation operator : a^{dager}_(i_particle) a_(i_hole) of spin = ispin
! on key_in
! ispin = 1 == alpha
! ispin = 2 == beta
! i_ok = 1 == the excitation is possible
! i_ok = -1 == the excitation is not possible
END_DOC
integer, intent(in) :: i_hole,i_particle,ispin
integer(bit_kind), intent(inout) :: key_in(N_int,2)
integer, intent(out) :: i_ok
integer :: k,j,i
use bitmasks
ASSERT (i_hole > 0 )
ASSERT (i_particle <= mo_tot_num)
i_ok = 1
! hole
k = shiftr(i_hole-1,bit_kind_shift)+1
j = i_hole-shiftl(k-1,bit_kind_shift)-1
key_in(k,ispin) = ibclr(key_in(k,ispin),j)
! particle
k = shiftr(i_particle-1,bit_kind_shift)+1
j = i_particle-shiftl(k-1,bit_kind_shift)-1
key_in(k,ispin) = ibset(key_in(k,ispin),j)
integer :: n_elec_tmp
n_elec_tmp = 0
do i = 1, N_int
n_elec_tmp += popcnt(key_in(i,1)) + popcnt(key_in(i,2))
enddo
if(n_elec_tmp .ne. elec_num)then
!print*, n_elec_tmp,elec_num
!call debug_det(key_in,N_int)
i_ok = -1
endif
end
subroutine do_spin_flip(key_in,i_flip,ispin,i_ok)
implicit none
BEGIN_DOC
! flip the spin ispin in the orbital i_flip
! on key_in
! ispin = 1 == alpha
! ispin = 2 == beta
! i_ok = 1 == the flip is possible
! i_ok = -1 == the flip is not possible
END_DOC
integer, intent(in) :: i_flip,ispin
integer(bit_kind), intent(inout) :: key_in(N_int,2)
integer, intent(out) :: i_ok
integer :: k,j,i
integer(bit_kind) :: key_tmp(N_int,2)
i_ok = -1
key_tmp = 0_bit_kind
k = shiftr(i_flip-1,bit_kind_shift)+1
j = i_flip-shiftl(k-1,bit_kind_shift)-1
key_tmp(k,1) = ibset(key_tmp(k,1),j)
integer :: other_spin(2)
other_spin(1) = 2
other_spin(2) = 1
if(popcnt(iand(key_tmp(k,1),key_in(k,ispin))) == 1 .and. popcnt(iand(key_tmp(k,1),key_in(k,other_spin(ispin)))) == 0 )then
! There is a spin "ispin" in the orbital i_flip AND There is no electron of opposit spin in the same orbital "i_flip"
key_in(k,ispin) = ibclr(key_in(k,ispin),j) ! destroy the electron ispin in the orbital i_flip
key_in(k,other_spin(ispin)) = ibset(key_in(k,other_spin(ispin)),j) ! create an electron of spin other_spin in the same orbital
implicit none
BEGIN_DOC
! Apply the mono excitation operator : a^{dager}_(i_particle) a_(i_hole) of spin = ispin
! on key_in
! ispin = 1 == alpha
! ispin = 2 == beta
! i_ok = 1 == the excitation is possible
! i_ok = -1 == the excitation is not possible
END_DOC
integer, intent(in) :: i_hole,i_particle,ispin
integer(bit_kind), intent(inout) :: key_in(N_int,2)
integer, intent(out) :: i_ok
integer :: k,j,i
use bitmasks
ASSERT (i_hole > 0 )
ASSERT (i_particle <= mo_tot_num)
i_ok = 1
else
return
endif
! hole
k = shiftr(i_hole-1,bit_kind_shift)+1
j = i_hole-shiftl(k-1,bit_kind_shift)-1
key_in(k,ispin) = ibclr(key_in(k,ispin),j)
! particle
k = shiftr(i_particle-1,bit_kind_shift)+1
j = i_particle-shiftl(k-1,bit_kind_shift)-1
key_in(k,ispin) = ibset(key_in(k,ispin),j)
integer :: n_elec_tmp
n_elec_tmp = 0
do i = 1, N_int
n_elec_tmp += popcnt(key_in(i,1)) + popcnt(key_in(i,2))
enddo
if(n_elec_tmp .ne. elec_num)then
!print*, n_elec_tmp,elec_num
!call debug_det(key_in,N_int)
i_ok = -1
endif
end
logical function is_spin_flip_possible(key_in,i_flip,ispin)
implicit none
BEGIN_DOC
! returns .True. if the spin-flip of spin ispin in the orbital i_flip is possible
! on key_in
END_DOC
integer, intent(in) :: i_flip,ispin
integer(bit_kind), intent(in) :: key_in(N_int,2)
integer :: k,j,i
integer(bit_kind) :: key_tmp(N_int,2)
is_spin_flip_possible = .False.
key_tmp = 0_bit_kind
k = shiftr(i_flip-1,bit_kind_shift)+1
j = i_flip-shiftl(k-1,bit_kind_shift)-1
key_tmp(k,1) = ibset(key_tmp(k,1),j)
integer :: other_spin(2)
other_spin(1) = 2
other_spin(2) = 1
if(popcnt(iand(key_tmp(k,1),key_in(k,ispin))) == 1 .and. popcnt(iand(key_tmp(k,1),key_in(k,other_spin(ispin)))) == 0 )then
! There is a spin "ispin" in the orbital i_flip AND There is no electron of opposit spin in the same orbital "i_flip"
is_spin_flip_possible = .True.
return
else
return
endif
implicit none
BEGIN_DOC
! returns .True. if the spin-flip of spin ispin in the orbital i_flip is possible
! on key_in
END_DOC
integer, intent(in) :: i_flip,ispin
integer(bit_kind), intent(in) :: key_in(N_int,2)
integer :: k,j,i
integer(bit_kind) :: key_tmp(N_int,2)
is_spin_flip_possible = .False.
key_tmp = 0_bit_kind
k = shiftr(i_flip-1,bit_kind_shift)+1
j = i_flip-shiftl(k-1,bit_kind_shift)-1
key_tmp(k,1) = ibset(key_tmp(k,1),j)
integer :: other_spin(2)
other_spin(1) = 2
other_spin(2) = 1
if(popcnt(iand(key_tmp(k,1),key_in(k,ispin))) == 1 .and. popcnt(iand(key_tmp(k,1),key_in(k,other_spin(ispin)))) == 0 )then
! There is a spin "ispin" in the orbital i_flip AND There is no electron of opposit spin in the same orbital "i_flip"
is_spin_flip_possible = .True.
return
else
return
endif
end
subroutine set_bit_to_integer(i_physical,key,Nint)
use bitmasks
implicit none
integer, intent(in) :: i_physical,Nint
integer(bit_kind), intent(inout) :: key(Nint)
integer :: k,j,i
k = shiftr(i_physical-1,bit_kind_shift)+1
j = i_physical-shiftl(k-1,bit_kind_shift)-1
key(k) = ibset(key(k),j)
use bitmasks
implicit none
integer, intent(in) :: i_physical,Nint
integer(bit_kind), intent(inout) :: key(Nint)
integer :: k,j,i
k = shiftr(i_physical-1,bit_kind_shift)+1
j = i_physical-shiftl(k-1,bit_kind_shift)-1
key(k) = ibset(key(k),j)
end
subroutine clear_bit_to_integer(i_physical,key,Nint)
use bitmasks
implicit none
integer, intent(in) :: i_physical,Nint
integer(bit_kind), intent(inout) :: key(Nint)
integer :: k,j,i
k = shiftr(i_physical-1,bit_kind_shift)+1
j = i_physical-shiftl(k-1,bit_kind_shift)-1
key(k) = ibclr(key(k),j)
use bitmasks
implicit none
integer, intent(in) :: i_physical,Nint
integer(bit_kind), intent(inout) :: key(Nint)
integer :: k,j,i
k = shiftr(i_physical-1,bit_kind_shift)+1
j = i_physical-shiftl(k-1,bit_kind_shift)-1
key(k) = ibclr(key(k),j)
end

18
src/Determinants/determinants.irp.f
View File

@ -308,24 +308,6 @@ END_PROVIDER
END_PROVIDER
subroutine flip_generators()
integer :: i,j,k
integer(bit_kind) :: detmp(N_int,2)
double precision :: tmp(N_states)
do i=1,N_det_generators/2
detmp(:,:) = psi_det_sorted(:,:,i)
tmp = psi_coef_sorted(i, :)
psi_det_sorted(:,:,i) = psi_det_sorted(:,:,N_det_generators+1-i)
psi_coef_sorted(i, :) = psi_coef_sorted(N_det_generators+1-i, :)
psi_det_sorted(:,:,N_det_generators+1-i) = detmp(:,:)
psi_coef_sorted(N_det_generators+1-i, :) = tmp
end do
TOUCH psi_det_sorted psi_coef_sorted psi_average_norm_contrib_sorted
end subroutine
BEGIN_PROVIDER [ integer(bit_kind), psi_det_sorted_bit, (N_int,2,psi_det_size) ]
&BEGIN_PROVIDER [ double precision, psi_coef_sorted_bit, (psi_det_size,N_states) ]
implicit none

223
src/Determinants/filter_connected.irp.f
View File

@ -198,12 +198,13 @@ end
subroutine getMobiles(key,key_mask, mobiles,Nint)
use bitmasks
implicit none
integer(bit_kind),intent(in) :: key(Nint,2), key_mask(Nint,2)
integer,intent(out) :: mobiles(2)
integer,intent(in) :: Nint
integer(bit_kind) :: mobileMask(Nint,2)
integer :: list(Nint*bit_kind_size), nel
integer :: list(Nint*bit_kind_size), nel,j
do j=1,Nint
mobileMask(j,1) = xor(key(j,1), key_mask(j,1))
@ -228,6 +229,7 @@ end subroutine
subroutine create_microlist(minilist, N_minilist, key_mask, microlist, idx_microlist, N_microlist, ptr_microlist, Nint)
use bitmasks
implicit none
integer, intent(in) :: Nint, N_minilist
integer(bit_kind), intent(in) :: minilist(Nint,2,N_minilist), key_mask(Nint,2)
@ -424,223 +426,4 @@ subroutine filter_connected_i_H_psi0(key1,key2,Nint,sze,idx)
idx(0) = l-1
end
subroutine filter_connected_i_H_psi0_SC2(key1,key2,Nint,sze,idx,idx_repeat)
use bitmasks
BEGIN_DOC
! standard filter_connected_i_H_psi but returns in addition
!
! the array of the index of the non connected determinants to key1
!
! in order to know what double excitation can be repeated on key1
!
! idx_repeat(0) is the number of determinants that can be used
!
! to repeat the excitations
END_DOC
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, intent(out) :: idx_repeat(0:sze)
integer :: i,l,l_repeat,m
integer :: degree_x2
ASSERT (Nint > 0)
ASSERT (Nint == N_int)
ASSERT (sze > 0)
integer :: degree
degree = popcnt(xor( ref_bitmask(1,1), key2(1,1))) + &
popcnt(xor( ref_bitmask(1,2), key2(1,2)))
!DIR$ NOUNROLL
do m=2,Nint
degree = degree+ popcnt(xor( ref_bitmask(m,1), key2(m,1))) + &
popcnt(xor( ref_bitmask(m,2), key2(m,2)))
enddo
degree = shiftr(degree,1)
l_repeat=1
l=1
if(degree == 2)then
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
elseif(degree_x2>6)then
idx_repeat(l_repeat) = i
l_repeat = l_repeat + 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
if(degree_x2 .ne. 0)then
idx(l) = i
l = l+1
endif
elseif(degree_x2>6)then
idx_repeat(l_repeat) = i
l_repeat = l_repeat + 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>6)then
idx_repeat(l_repeat) = i
l_repeat = l_repeat + 1
else 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
!DIR$ LOOP COUNT MIN(4)
do m=1,Nint
degree_x2 = degree_x2+ popcnt(xor( key1(m,1,i), key2(m,1))) +&
popcnt(xor( key1(m,2,i), key2(m,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
elseif(degree_x2>6)then
idx_repeat(l_repeat) = i
l_repeat = l_repeat + 1
endif
enddo
endif
elseif(degree==1)then
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
else
idx_repeat(l_repeat) = i
l_repeat = l_repeat + 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
if(degree_x2 .ne. 0)then
idx(l) = i
l = l+1
endif
else
idx_repeat(l_repeat) = i
l_repeat = l_repeat + 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
if(degree_x2 .ne. 0)then
idx(l) = i
l = l+1
endif
else
idx_repeat(l_repeat) = i
l_repeat = l_repeat + 1
endif
enddo
else
!DIR$ LOOP COUNT (1000)
do i=1,sze
degree_x2 = 0
!DIR$ LOOP COUNT MIN(4)
do m=1,Nint
degree_x2 = degree_x2+ popcnt(xor( key1(m,1,i), key2(m,1))) +&
popcnt(xor( key1(m,2,i), key2(m,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
else
idx_repeat(l_repeat) = i
l_repeat = l_repeat + 1
endif
enddo
endif
else
! print*,'more than a double excitation, can not apply the '
! print*,'SC2 dressing of the diagonal element .....'
! print*,'stop !!'
! print*,'degree = ',degree
! stop
idx(0) = 0
idx_repeat(0) = 0
endif
idx(0) = l-1
idx_repeat(0) = l_repeat-1
end

118
src/Determinants/s2.irp.f
View File

@ -328,124 +328,6 @@ subroutine get_uJ_s2_uI(psi_keys_tmp,psi_coefs_tmp,n,nmax_coefs,nmax_keys,s2,nst
enddo
end
subroutine diagonalize_s2_betweenstates(keys_tmp,u_0,n,nmax_keys,nmax_coefs,nstates,s2_eigvalues)
BEGIN_DOC
! You enter with nstates vectors in u_0 that may be coupled by S^2
! The subroutine diagonalize the S^2 operator in the basis of these states.
! The vectors that you obtain in output are no more coupled by S^2,
! which does not necessary mean that they are eigenfunction of S^2.
! n,nmax,nstates = number of determinants, physical dimension of the arrays and number of states
! keys_tmp = array of integer(bit_kind) that represents the determinants
! psi_coefs(i,j) = coeff of the ith determinant in the jth state
! VECTORS ARE SUPPOSED TO BE ORTHONORMAL IN INPUT
END_DOC
implicit none
use bitmasks
integer, intent(in) :: n,nmax_keys,nmax_coefs,nstates
integer(bit_kind), intent(in) :: keys_tmp(N_int,2,nmax_keys)
double precision, intent(inout) :: u_0(nmax_coefs,nstates)
double precision, intent(out) :: s2_eigvalues(nstates)
double precision,allocatable :: s2(:,:),overlap(:,:)
double precision, allocatable :: eigvectors(:,:,:)
integer :: i,j,k
double precision, allocatable :: psi_coefs_tmp(:,:)
double precision :: accu,coef_contract
double precision :: u_dot_u,u_dot_v
print*,''
print*,'*********************************************************************'
print*,'Cleaning the various vectors by diagonalization of the S^2 matrix ...'
print*,''
print*,'nstates = ',nstates
allocate(s2(nstates,nstates),overlap(nstates,nstates))
overlap = 0.d0
call dgemm('T','N',nstates,nstates,n, 1.d0, u_0, size(u_0,1), &
u_0, size(u_0,1), 0.d0, overlap, size(overlap,1))
call ortho_lowdin(overlap,size(overlap,1),nstates,u_0,size(u_0,1),n)
double precision, allocatable :: v_0(:,:)
allocate ( v_0(size(u_0,1),nstates) )
call S2_u_0_nstates(v_0,u_0,n,keys_tmp,N_int,nstates,size(u_0,1))
call dgemm('T','N',nstates,nstates,n, 1.d0, u_0, size(u_0,1), &
v_0, size(v_0,1), 0.d0, s2, size(s2,1))
print*,'S^2 matrix in the basis of the states considered'
do i = 1, nstates
write(*,'(100(F5.2,1X))')s2(i,:)
enddo
double precision :: accu_precision_diag,accu_precision_of_diag
accu_precision_diag = 0.d0
accu_precision_of_diag = 0.d0
do i = 1, nstates
! Do not combine states of the same spin symmetry
do j = i+1, nstates
if( dabs(s2(i,i) - s2(j,j)) .le.0.5d0) then
s2(i,j) = 0.d0
s2(j,i) = 0.d0
endif
enddo
! Do not rotate if the diagonal is correct
if( dabs(s2(i,i) - expected_s2).le.5.d-3) then
do j = 1, nstates
if (j==i) cycle
s2(i,j) = 0.d0
s2(j,i) = 0.d0
enddo
endif
enddo
print*,'Modified S^2 matrix that will be diagonalized'
do i = 1, nstates
write(*,'(10(F5.2,1X))')s2(i,:)
s2(i,i) = s2(i,i)
enddo
allocate(eigvectors(nstates,nstates,2))
! call svd(s2, size(s2,1), eigvectors, size(eigvectors,1), s2_eigvalues,&
! eigvectors(1,1,2), size(eigvectors,1), nstates, nstates)
call lapack_diagd(s2_eigvalues,eigvectors,s2,nstates,nstates)
print*,'Eigenvalues'
double precision :: t(nstates)
integer :: iorder(nstates)
do i = 1, nstates
t(i) = dabs(s2_eigvalues(i))
iorder(i) = i
enddo
call dsort(t,iorder,nstates)
do i = 1, nstates
s2_eigvalues(i) = t(i)
do j=1,nstates
eigvectors(j,i,2) = eigvectors(j,iorder(i),1)
enddo
print*,'s2 = ',s2_eigvalues(i)
enddo
allocate(psi_coefs_tmp(nmax_coefs,nstates))
psi_coefs_tmp = 0.d0
do j = 1, nstates
do k = 1, nstates
coef_contract = eigvectors(k,j,2) ! <phi_k|Psi_j>
do i = 1, n_det
psi_coefs_tmp(i,j) += u_0(i,k) * coef_contract
enddo
enddo
enddo
do j = 1, nstates
accu = 1.d0/u_dot_u(psi_coefs_tmp(1,j),n_det)
do i = 1, n_det
u_0(i,j) = psi_coefs_tmp(i,j) * accu
enddo
enddo
deallocate(s2,v_0,eigvectors,psi_coefs_tmp,overlap )
end
subroutine i_S2_psi_minilist(key,keys,idx_key,N_minilist,coef,Nint,Ndet,Ndet_max,Nstate,i_S2_psi_array)
use bitmasks

977
src/Determinants/slater_rules.irp.f
View File

File diff suppressed because it is too large Load Diff

618
src/Determinants/two_body_dm_map.irp.f
View File

@ -1,618 +0,0 @@
use map_module
BEGIN_PROVIDER [ type(map_type), two_body_dm_ab_map ]
implicit none
BEGIN_DOC
! Map of the two body density matrix elements for the alpha/beta elements
END_DOC
integer(key_kind) :: key_max
integer(map_size_kind) :: sze
use map_module
call bielec_integrals_index(mo_tot_num,mo_tot_num,mo_tot_num,mo_tot_num,key_max)
sze = key_max
call map_init(two_body_dm_ab_map,sze)
print*, 'two_body_dm_ab_map initialized'
END_PROVIDER
subroutine insert_into_two_body_dm_ab_map(n_product,buffer_i, buffer_values, thr)
use map_module
implicit none
BEGIN_DOC
! Create new entry into two_body_dm_ab_map, or accumulate in an existing entry
END_DOC
integer, intent(in) :: n_product
integer(key_kind), intent(inout) :: buffer_i(n_product)
real(integral_kind), intent(inout) :: buffer_values(n_product)
real(integral_kind), intent(in) :: thr
call map_update(two_body_dm_ab_map, buffer_i, buffer_values, n_product, thr)
end
double precision function get_two_body_dm_ab_map_element(i,j,k,l,map)
use map_module
implicit none
BEGIN_DOC
! Returns one value of the wo body density matrix \rho_{ijkl}^{\alpha \beta} defined as :
! \rho_{ijkl}^{\alpha \beta } = <\Psi|a^{\dagger}_{i\alpha} a^{\dagger}_{j\beta} a_{k\beta} a_{l\alpha}|\Psi>
END_DOC
PROVIDE two_body_dm_ab_map
integer, intent(in) :: i,j,k,l
integer(key_kind) :: idx
type(map_type), intent(inout) :: map
real(integral_kind) :: tmp
PROVIDE two_body_dm_in_map
!DIR$ FORCEINLINE
call bielec_integrals_index(i,j,k,l,idx)
!DIR$ FORCEINLINE
call map_get(two_body_dm_ab_map,idx,tmp)
get_two_body_dm_ab_map_element = dble(tmp)
end
subroutine get_get_two_body_dm_ab_map_elements(j,k,l,sze,out_val,map)
use map_module
implicit none
BEGIN_DOC
! Returns multiple elements of the \rho_{ijkl}^{\alpha \beta }, all
! i for j,k,l fixed.
END_DOC
integer, intent(in) :: j,k,l, sze
double precision, intent(out) :: out_val(sze)
type(map_type), intent(inout) :: map
integer :: i
integer(key_kind) :: hash(sze)
real(integral_kind) :: tmp_val(sze)
PROVIDE two_body_dm_in_map
do i=1,sze
!DIR$ FORCEINLINE
call bielec_integrals_index(i,j,k,l,hash(i))
enddo
if (key_kind == 8) then
call map_get_many(two_body_dm_ab_map, hash, out_val, sze)
else
call map_get_many(two_body_dm_ab_map, hash, tmp_val, sze)
! Conversion to double precision
do i=1,sze
out_val(i) = dble(tmp_val(i))
enddo
endif
end
BEGIN_PROVIDER [ logical, two_body_dm_in_map ]
implicit none
BEGIN_DOC
! If True, the map of the two body density matrix alpha/beta is provided
END_DOC
two_body_dm_in_map = .True.
call add_values_to_two_body_dm_map(full_ijkl_bitmask_4)
END_PROVIDER
subroutine add_values_to_two_body_dm_map(mask_ijkl)
use bitmasks
use map_module
implicit none
BEGIN_DOC
! Adds values to the map of two_body_dm according to some bitmask
END_DOC
integer(bit_kind), intent(in) :: mask_ijkl(N_int,4)
integer :: degree
PROVIDE mo_coef psi_coef psi_det
integer :: exc(0:2,2,2)
integer :: h1,h2,p1,p2,s1,s2
double precision :: phase
double precision :: contrib
integer(key_kind),allocatable :: buffer_i(:)
double precision ,allocatable :: buffer_value(:)
integer :: size_buffer
integer :: n_elements
integer :: occ(N_int*bit_kind_size,2)
integer :: n_occ_ab(2)
integer :: i,j,k,l,m
size_buffer = min(mo_tot_num*mo_tot_num*mo_tot_num,16000000)
allocate(buffer_i(size_buffer),buffer_value(size_buffer))
n_elements = 0
do i = 1, N_det ! i == |I>
call bitstring_to_list_ab(psi_det(1,1,i), occ, n_occ_ab, N_int)
do j = i+1, N_det ! j == <J|
call get_excitation_degree(psi_det(1,1,i),psi_det(1,1,j),degree,N_int)
if(degree>2)cycle
call get_excitation(psi_det(1,1,i),psi_det(1,1,j),exc,degree,phase,N_int)
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
if(degree==2)then ! case of the DOUBLE EXCITATIONS ************************************
if(s1==s2)cycle ! Only the alpha/beta two body density matrix
! <J| a^{\dagger}_{p1 s1} a^{\dagger}_{p2 s2} a_{h2 s2} a_{h1 s1} |I> * c_I * c_J
if(h1>p1)cycle
if(h2>p2)cycle
! if(s1.ne.1)cycle
n_elements += 1
contrib = psi_coef(i,1) * psi_coef(j,1) * phase
buffer_value(n_elements) = contrib
!DIR$ FORCEINLINE
call mo_bielec_integrals_index(h1,h2,p1,p2,buffer_i(n_elements))
! if (n_elements == size_buffer) then
! call insert_into_two_body_dm_ab_map(n_elements,buffer_i,buffer_value,&
! real(mo_integrals_threshold,integral_kind))
! n_elements = 0
! endif
else ! case of the SINGLE EXCITATIONS ***************************************************
cycle
! if(s1==1)then ! Mono alpha :
! do k = 1, elec_beta_num
! m = occ(k,2)
! n_elements += 1
! buffer_value(n_elements) = contrib
! ! <J| a^{\dagger}_{p1 \alpha} \hat{n}_{m \beta} a_{h1 \alpha} |I> * c_I * c_J
! call mo_bielec_integrals_index(h1,m,p1,m,buffer_i(n_elements))
! if (n_elements == size_buffer) then
! call insert_into_two_body_dm_ab_map(n_elements,buffer_i,buffer_value,&
! real(mo_integrals_threshold,integral_kind))
! n_elements = 0
! endif
! enddo
! else ! Mono Beta :
! do k = 1, elec_alpha_num
! m = occ(k,1)
! n_elements += 1
! buffer_value(n_elements) = contrib
! ! <J| a^{\dagger}_{p1 \beta} \hat{n}_{m \alpha} a_{h1 \beta} |I> * c_I * c_J
! call mo_bielec_integrals_index(h1,m,p1,m,buffer_i(n_elements))
! if (n_elements == size_buffer) then
! call insert_into_two_body_dm_ab_map(n_elements,buffer_i,buffer_value,&
! real(mo_integrals_threshold,integral_kind))
! n_elements = 0
! endif
! enddo
! endif
endif
enddo
enddo
print*,'n_elements = ',n_elements
call insert_into_two_body_dm_ab_map(n_elements,buffer_i,buffer_value,&
real(mo_integrals_threshold,integral_kind))
call map_merge(two_body_dm_ab_map)
deallocate(buffer_i,buffer_value)
end
BEGIN_PROVIDER [double precision, two_body_dm_ab_diag_act, (n_act_orb, n_act_orb)]
&BEGIN_PROVIDER [double precision, two_body_dm_ab_diag_inact, (n_inact_orb_allocate, n_inact_orb_allocate)]
&BEGIN_PROVIDER [double precision, two_body_dm_ab_diag_core, (n_core_orb_allocate, n_core_orb_allocate)]
&BEGIN_PROVIDER [double precision, two_body_dm_ab_diag_all, (mo_tot_num, mo_tot_num)]
&BEGIN_PROVIDER [double precision, two_body_dm_diag_core_a_act_b, (n_core_orb_allocate,n_act_orb)]
&BEGIN_PROVIDER [double precision, two_body_dm_diag_core_b_act_a, (n_core_orb_allocate,n_act_orb)]
&BEGIN_PROVIDER [double precision, two_body_dm_diag_core_act, (n_core_orb_allocate,n_act_orb)]
implicit none
use bitmasks
integer :: i,j,k,l,m
integer :: occ(N_int*bit_kind_size,2)
integer :: n_occ_ab(2)
integer :: occ_act(N_int*bit_kind_size,2)
integer :: n_occ_ab_act(2)
integer :: occ_core(N_int*bit_kind_size,2)
integer :: n_occ_ab_core(2)
double precision :: contrib
BEGIN_DOC
! two_body_dm_ab_diag_all(k,m) = <\Psi | n_(k\alpha) n_(m\beta) | \Psi>
! two_body_dm_ab_diag_act(k,m) is restricted to the active orbitals :
! orbital k = list_act(k)
! two_body_dm_ab_diag_inact(k,m) is restricted to the inactive orbitals :
! orbital k = list_inact(k)
! two_body_dm_ab_diag_core(k,m) is restricted to the core orbitals :
! orbital k = list_core(k)
! two_body_dm_ab_diag_core_b_act_a(k,m) represents the core beta <-> active alpha part of the two body dm
! orbital k = list_core(k)
! orbital m = list_act(m)
! two_body_dm_ab_diag_core_a_act_b(k,m) represents the core alpha <-> active beta part of the two body dm
! orbital k = list_core(k)
! orbital m = list_act(m)
! two_body_dm_ab_diag_core_act(k,m) represents the core<->active part of the diagonal two body dm
! when we traced on the spin
! orbital k = list_core(k)
! orbital m = list_act(m)
END_DOC
integer(bit_kind) :: key_tmp_core(N_int,2)
integer(bit_kind) :: key_tmp_act(N_int,2)
two_body_dm_ab_diag_all = 0.d0
two_body_dm_ab_diag_act = 0.d0
two_body_dm_ab_diag_core = 0.d0
two_body_dm_ab_diag_inact = 0.d0
two_body_dm_diag_core_a_act_b = 0.d0
two_body_dm_diag_core_b_act_a = 0.d0
two_body_dm_diag_core_act = 0.d0
do i = 1, N_det ! i == |I>
! Full diagonal part of the two body dm
contrib = psi_coef(i,1)**2
call bitstring_to_list_ab(psi_det(1,1,i), occ, n_occ_ab, N_int)
do j = 1, elec_beta_num
k = occ(j,2)
do l = 1, elec_alpha_num
m = occ(l,1)
two_body_dm_ab_diag_all(k,m) += 0.5d0 * contrib
two_body_dm_ab_diag_all(m,k) += 0.5d0 * contrib
enddo
enddo
! ACTIVE PART of the diagonal part of the two body dm
do j = 1, N_int
key_tmp_act(j,1) = psi_det(j,1,i)
key_tmp_act(j,2) = psi_det(j,2,i)
enddo
do j = 1, N_int
key_tmp_act(j,1) = iand(key_tmp_act(j,1),cas_bitmask(j,1,1))
key_tmp_act(j,2) = iand(key_tmp_act(j,2),cas_bitmask(j,1,1))
enddo
call bitstring_to_list_ab(key_tmp_act, occ_act, n_occ_ab_act, N_int)
do j = 1,n_occ_ab_act(2)
k = list_act_reverse(occ_act(j,2))
do l = 1, n_occ_ab_act(1)
m = list_act_reverse(occ_act(l,1))
two_body_dm_ab_diag_act(k,m) += 0.5d0 * contrib
two_body_dm_ab_diag_act(m,k) += 0.5d0 * contrib
enddo
enddo
! CORE PART of the diagonal part of the two body dm
do j = 1, N_int
key_tmp_core(j,1) = psi_det(j,1,i)
key_tmp_core(j,2) = psi_det(j,2,i)
enddo
do j = 1, N_int
key_tmp_core(j,1) = iand(key_tmp_core(j,1),core_bitmask(j,1))
key_tmp_core(j,2) = iand(key_tmp_core(j,2),core_bitmask(j,1))
enddo
call bitstring_to_list_ab(key_tmp_core, occ_core, n_occ_ab_core, N_int)
do j = 1,n_occ_ab_core(2)
k = list_core_reverse(occ_core(j,2))
do l = 1, n_occ_ab_core(1)
m = list_core_reverse(occ_core(l,1))
two_body_dm_ab_diag_core(k,m) += 0.5d0 * contrib
two_body_dm_ab_diag_core(m,k) += 0.5d0 * contrib
enddo
enddo
! ACT<->CORE PART
! alpha electron in active space
do j = 1,n_occ_ab_act(1)
k = list_act_reverse(occ_act(j,1))
! beta electron in core space
do l = 1, n_occ_ab_core(2)
m = list_core_reverse(occ_core(l,2))
! The fact that you have 1 * contrib and not 0.5 * contrib
! takes into account the following symmetry :
! 0.5 * <n_k n_m> + 0.5 * <n_m n_k>
two_body_dm_diag_core_b_act_a(m,k) += contrib
enddo
enddo
! beta electron in active space
do j = 1,n_occ_ab_act(2)
k = list_act_reverse(occ_act(j,2))
! alpha electron in core space
do l = 1, n_occ_ab_core(1)
m = list_core_reverse(occ_core(l,1))
! The fact that you have 1 * contrib and not 0.5 * contrib
! takes into account the following symmetry :
! 0.5 * <n_k n_m> + 0.5 * <n_m n_k>
two_body_dm_diag_core_a_act_b(m,k) += contrib
enddo
enddo
enddo
do j = 1, n_core_orb
do l = 1, n_act_orb
two_body_dm_diag_core_act(j,l) = two_body_dm_diag_core_b_act_a(j,l) + two_body_dm_diag_core_a_act_b(j,l)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, two_body_dm_ab_big_array_act, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)]
&BEGIN_PROVIDER [double precision, two_body_dm_ab_big_array_core_act, (n_core_orb_allocate,n_act_orb,n_act_orb)]
implicit none
use bitmasks
integer :: i,j,k,l,m
integer :: degree
PROVIDE mo_coef psi_coef psi_det
integer :: exc(0:2,2,2)
integer :: h1,h2,p1,p2,s1,s2
double precision :: phase
double precision :: contrib
integer :: occ(N_int*bit_kind_size,2)
integer :: n_occ_ab(2)
integer :: occ_core(N_int*bit_kind_size,2)
integer :: n_occ_ab_core(2)
integer(bit_kind) :: key_tmp_i(N_int,2)
integer(bit_kind) :: key_tmp_i_core(N_int,2)
integer(bit_kind) :: key_tmp_j(N_int,2)
two_body_dm_ab_big_array_act = 0.d0
two_body_dm_ab_big_array_core_act = 0.d0
BEGIN_DOC
! two_body_dm_ab_big_array_act = Purely active part of the two body density matrix
! two_body_dm_ab_big_array_act_core takes only into account the single excitation
! within the active space that adds terms in the act <-> core two body dm
! two_body_dm_ab_big_array_act_core(i,j,k) = < a^\dagger_i n_k a_j >
! with i,j in the ACTIVE SPACE
! with k in the CORE SPACE
!
! The alpha-beta extra diagonal energy FOR WF DEFINED AS AN APPROXIMATION OF A CAS can be computed thanks to
! sum_{h1,p1,h2,p2} two_body_dm_ab_big_array_act(h1,p1,h2,p2) * (h1p1|h2p2)
! + sum_{h1,p1,h2,p2} two_body_dm_ab_big_array_core_act(h1,p1,h2,p2) * (h1p1|h2p2)
END_DOC
do i = 1, N_det ! i == |I>
! active part of psi_det(i)
do j = 1, N_int
key_tmp_i(j,1) = psi_det(j,1,i)
key_tmp_i(j,2) = psi_det(j,2,i)
key_tmp_i_core(j,1) = psi_det(j,1,i)
key_tmp_i_core(j,2) = psi_det(j,2,i)
enddo
do j = 1, N_int
key_tmp_i(j,1) = iand(key_tmp_i(j,1),cas_bitmask(j,1,1))
key_tmp_i(j,2) = iand(key_tmp_i(j,2),cas_bitmask(j,1,1))
enddo
do j = 1, N_int
key_tmp_i_core(j,1) = iand(key_tmp_i_core(j,1),core_bitmask(j,1))
key_tmp_i_core(j,2) = iand(key_tmp_i_core(j,2),core_bitmask(j,1))
enddo
call bitstring_to_list_ab(key_tmp_i_core, occ_core, n_occ_ab_core, N_int)
call bitstring_to_list_ab(key_tmp_i, occ, n_occ_ab, N_int)
do j = i+1, N_det ! j == <J|
! active part of psi_det(j)
do k = 1, N_int
key_tmp_j(k,1) = psi_det(k,1,j)
key_tmp_j(k,2) = psi_det(k,2,j)
enddo
do k = 1, N_int
key_tmp_j(k,1) = iand(key_tmp_j(k,1),cas_bitmask(k,1,1))
key_tmp_j(k,2) = iand(key_tmp_j(k,2),cas_bitmask(k,1,1))
enddo
! control if the two determinants are connected by
! at most a double excitation WITHIN THE ACTIVE SPACE
call get_excitation_degree(key_tmp_i,key_tmp_j,degree,N_int)
if(degree>2)cycle
! if it is the case, then compute the hamiltonian matrix element with the proper phase
call get_excitation(psi_det(1,1,i),psi_det(1,1,j),exc,degree,phase,N_int)
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
contrib = 0.5d0 * psi_coef(i,1) * psi_coef(j,1) * phase
if(degree==2)then ! case of the DOUBLE EXCITATIONS ************************************
if(s1==s2)cycle ! Only the alpha/beta two body density matrix
! <J| a^{\dagger}_{p1 s1} a^{\dagger}_{p2 s2} a_{h2 s2} a_{h1 s1} |I> * c_I * c_J
h1 = list_act_reverse(h1)
h2 = list_act_reverse(h2)
p1 = list_act_reverse(p1)
p2 = list_act_reverse(p2)
call insert_into_two_body_dm_big_array( two_body_dm_ab_big_array_act,n_act_orb,n_act_orb,n_act_orb,n_act_orb,contrib,h1,p1,h2,p2)
else if(degree==1)then! case of the SINGLE EXCITATIONS ***************************************************
print*,'h1 = ',h1
h1 = list_act_reverse(h1)
print*,'h1 = ',h1
print*,'p1 = ',p1
p1 = list_act_reverse(p1)
print*,'p1 = ',p1
if(s1==1)then ! Mono alpha :
! purely active part of the extra diagonal two body dm
do k = 1, n_occ_ab(2)
m = list_act_reverse(occ(k,2))
! <J| a^{\dagger}_{p1 \alpha} \hat{n}_{m \beta} a_{h1 \alpha} |I> * c_I * c_J
call insert_into_two_body_dm_big_array( two_body_dm_ab_big_array_act,n_act_orb,n_act_orb,n_act_orb,n_act_orb,contrib,h1,p1,m,m)
enddo
! core <-> active part of the extra diagonal two body dm
do k = 1, n_occ_ab_core(2)
m = list_core_reverse(occ_core(k,2))
! <J| a^{\dagger}_{p1 \alpha} \hat{n}_{m \beta} a_{h1 \alpha} |I> * c_I * c_J
two_body_dm_ab_big_array_core_act(m,h1,p1) += 2.d0 * contrib
two_body_dm_ab_big_array_core_act(m,p1,h1) += 2.d0 * contrib
enddo
else ! Mono Beta :
! purely active part of the extra diagonal two body dm
do k = 1, n_occ_ab(1)
m = list_act_reverse(occ(k,1))
! <J| a^{\dagger}_{p1 \beta} \hat{n}_{m \alpha} a_{h1 \beta} |I> * c_I * c_J
call insert_into_two_body_dm_big_array(two_body_dm_ab_big_array_act,n_act_orb,n_act_orb,n_act_orb,n_act_orb,contrib,h1,p1,m,m)
enddo
! core <-> active part of the extra diagonal two body dm
do k = 1, n_occ_ab_core(1)
m = list_core_reverse(occ_core(k,1))
! <J| a^{\dagger}_{p1 \alpha} \hat{n}_{m \beta} a_{h1 \alpha} |I> * c_I * c_J
two_body_dm_ab_big_array_core_act(m,h1,p1) += 2.d0 * contrib
two_body_dm_ab_big_array_core_act(m,p1,h1) += 2.d0 * contrib
enddo
endif
endif
enddo
enddo
print*,'Big array for density matrix provided !'
END_PROVIDER
subroutine insert_into_two_body_dm_big_array(big_array,dim1,dim2,dim3,dim4,contrib,h1,p1,h2,p2)
implicit none
integer, intent(in) :: h1,p1,h2,p2
integer, intent(in) :: dim1,dim2,dim3,dim4
double precision, intent(inout) :: big_array(dim1,dim2,dim3,dim4)
double precision :: contrib
! Two spin symmetry
big_array(h1,p1,h2,p2) += contrib
big_array(h2,p2,h1,p1) += contrib
! Hermicity : hole-particle symmetry
big_array(p1,h1,p2,h2) += contrib
big_array(p2,h2,p1,h1) += contrib
end
double precision function compute_extra_diag_two_body_dm_ab(r1,r2)
implicit none
BEGIN_DOC
! compute the extra diagonal contribution to the alpha/bet two body density at r1, r2
END_DOC
double precision :: r1(3), r2(3)
double precision :: compute_extra_diag_two_body_dm_ab_act,compute_extra_diag_two_body_dm_ab_core_act
compute_extra_diag_two_body_dm_ab = compute_extra_diag_two_body_dm_ab_act(r1,r2)+compute_extra_diag_two_body_dm_ab_core_act(r1,r2)
end
double precision function compute_extra_diag_two_body_dm_ab_act(r1,r2)
implicit none
BEGIN_DOC
! compute the extra diagonal contribution to the two body density at r1, r2
! involving ONLY THE ACTIVE PART, which means that the four index of the excitations
! involved in the two body density matrix are ACTIVE
END_DOC
PROVIDE n_act_orb
double precision, intent(in) :: r1(3),r2(3)
integer :: i,j,k,l
double precision :: mos_array_r1(n_act_orb),mos_array_r2(n_act_orb)
double precision :: contrib
double precision :: contrib_tmp
!print*,'n_act_orb = ',n_act_orb
compute_extra_diag_two_body_dm_ab_act = 0.d0
call give_all_act_mos_at_r(r1,mos_array_r1)
call give_all_act_mos_at_r(r2,mos_array_r2)
do l = 1, n_act_orb ! p2
do k = 1, n_act_orb ! h2
do j = 1, n_act_orb ! p1
do i = 1,n_act_orb ! h1
contrib_tmp = mos_array_r1(i) * mos_array_r1(j) * mos_array_r2(k) * mos_array_r2(l)
compute_extra_diag_two_body_dm_ab_act += two_body_dm_ab_big_array_act(i,j,k,l) * contrib_tmp
enddo
enddo
enddo
enddo
end
double precision function compute_extra_diag_two_body_dm_ab_core_act(r1,r2)
implicit none
BEGIN_DOC
! compute the extra diagonal contribution to the two body density at r1, r2
! involving ONLY THE ACTIVE PART, which means that the four index of the excitations
! involved in the two body density matrix are ACTIVE
END_DOC
double precision, intent(in) :: r1(3),r2(3)
integer :: i,j,k,l
double precision :: mos_array_act_r1(n_act_orb),mos_array_act_r2(n_act_orb)
double precision :: mos_array_core_r1(n_core_orb),mos_array_core_r2(n_core_orb)
double precision :: contrib_core_1,contrib_core_2
double precision :: contrib_act_1,contrib_act_2
double precision :: contrib_tmp
compute_extra_diag_two_body_dm_ab_core_act = 0.d0
call give_all_act_mos_at_r(r1,mos_array_act_r1)
call give_all_act_mos_at_r(r2,mos_array_act_r2)
call give_all_core_mos_at_r(r1,mos_array_core_r1)
call give_all_core_mos_at_r(r2,mos_array_core_r2)
do i = 1, n_act_orb ! h1
do j = 1, n_act_orb ! p1
contrib_act_1 = mos_array_act_r1(i) * mos_array_act_r1(j)
contrib_act_2 = mos_array_act_r2(i) * mos_array_act_r2(j)
do k = 1,n_core_orb ! h2
contrib_core_1 = mos_array_core_r1(k) * mos_array_core_r1(k)
contrib_core_2 = mos_array_core_r2(k) * mos_array_core_r2(k)
contrib_tmp = 0.5d0 * (contrib_act_1 * contrib_core_2 + contrib_act_2 * contrib_core_1)
compute_extra_diag_two_body_dm_ab_core_act += two_body_dm_ab_big_array_core_act(k,i,j) * contrib_tmp
enddo
enddo
enddo
end
double precision function compute_diag_two_body_dm_ab_core(r1,r2)
implicit none
double precision :: r1(3),r2(3)
integer :: i,j,k,l
double precision :: mos_array_r1(n_core_orb_allocate),mos_array_r2(n_core_orb_allocate)
double precision :: contrib,contrib_tmp
compute_diag_two_body_dm_ab_core = 0.d0
call give_all_core_mos_at_r(r1,mos_array_r1)
call give_all_core_mos_at_r(r2,mos_array_r2)
do l = 1, n_core_orb !
contrib = mos_array_r2(l)*mos_array_r2(l)
! if(dabs(contrib).lt.threshld_two_bod_dm)cycle
do k = 1, n_core_orb !
contrib_tmp = contrib * mos_array_r1(k)*mos_array_r1(k)
! if(dabs(contrib).lt.threshld_two_bod_dm)cycle
compute_diag_two_body_dm_ab_core += two_body_dm_ab_diag_core(k,l) * contrib_tmp
enddo
enddo
end
double precision function compute_diag_two_body_dm_ab_act(r1,r2)
implicit none
double precision :: r1(3),r2(3)
integer :: i,j,k,l
double precision :: mos_array_r1(n_act_orb),mos_array_r2(n_act_orb)
double precision :: contrib,contrib_tmp
compute_diag_two_body_dm_ab_act = 0.d0
call give_all_act_mos_at_r(r1,mos_array_r1)
call give_all_act_mos_at_r(r2,mos_array_r2)
do l = 1, n_act_orb !
contrib = mos_array_r2(l)*mos_array_r2(l)
! if(dabs(contrib).lt.threshld_two_bod_dm)cycle
do k = 1, n_act_orb !
contrib_tmp = contrib * mos_array_r1(k)*mos_array_r1(k)
! if(dabs(contrib).lt.threshld_two_bod_dm)cycle
compute_diag_two_body_dm_ab_act += two_body_dm_ab_diag_act(k,l) * contrib_tmp
enddo
enddo
end
double precision function compute_diag_two_body_dm_ab_core_act(r1,r2)
implicit none
double precision :: r1(3),r2(3)
integer :: i,j,k,l
double precision :: mos_array_core_r1(n_core_orb_allocate),mos_array_core_r2(n_core_orb_allocate)
double precision :: mos_array_act_r1(n_act_orb),mos_array_act_r2(n_act_orb)
double precision :: contrib_core_1,contrib_core_2
double precision :: contrib_act_1,contrib_act_2
double precision :: contrib_tmp
compute_diag_two_body_dm_ab_core_act = 0.d0
call give_all_act_mos_at_r(r1,mos_array_act_r1)
call give_all_act_mos_at_r(r2,mos_array_act_r2)
call give_all_core_mos_at_r(r1,mos_array_core_r1)
call give_all_core_mos_at_r(r2,mos_array_core_r2)
! if(dabs(contrib).lt.threshld_two_bod_dm)cycle
do k = 1, n_act_orb !
contrib_act_1 = mos_array_act_r1(k) * mos_array_act_r1(k)
contrib_act_2 = mos_array_act_r2(k) * mos_array_act_r2(k)
contrib_tmp = 0.5d0 * (contrib_act_1 * contrib_act_2 + contrib_act_2 * contrib_act_1)
! if(dabs(contrib).lt.threshld_two_bod_dm)cycle
do l = 1, n_core_orb !
contrib_core_1 = mos_array_core_r1(l) * mos_array_core_r1(l)
contrib_core_2 = mos_array_core_r2(l) * mos_array_core_r2(l)
compute_diag_two_body_dm_ab_core_act += two_body_dm_diag_core_act(l,k) * contrib_tmp
enddo
enddo
end
double precision function compute_diag_two_body_dm_ab(r1,r2)
implicit none
double precision,intent(in) :: r1(3),r2(3)
double precision :: compute_diag_two_body_dm_ab_act,compute_diag_two_body_dm_ab_core
double precision :: compute_diag_two_body_dm_ab_core_act
compute_diag_two_body_dm_ab = compute_diag_two_body_dm_ab_act(r1,r2)+compute_diag_two_body_dm_ab_core(r1,r2) &
+ compute_diag_two_body_dm_ab_core_act(r1,r2)
end

286
src/Determinants/useful_for_ovb.irp.f
View File

@ -1,286 +0,0 @@
integer function n_open_shell(det_in,nint)
implicit none
use bitmasks
integer, intent(in) :: nint
integer(bit_kind), intent(in) :: det_in(nint,2)
integer :: i
n_open_shell = 0
do i=1,Nint
n_open_shell += popcnt(iand(xor(det_in(i,1),det_in(i,2)),det_in(i,1)))
enddo
end
integer function n_closed_shell(det_in,nint)
implicit none
use bitmasks
integer, intent(in) :: nint
integer(bit_kind), intent(in) :: det_in(nint,2)
integer :: i
n_closed_shell = 0
do i=1,Nint
n_closed_shell += popcnt(iand(det_in(i,1),det_in(i,2)))
enddo
end
integer function n_closed_shell_cas(det_in,nint)
implicit none
use bitmasks
integer, intent(in) :: nint
integer(bit_kind), intent(in) :: det_in(nint,2)
integer(bit_kind) :: det_tmp(nint,2)
integer :: i
n_closed_shell_cas = 0
do i=1,Nint
det_tmp(i,1) = xor(det_in(i,1),reunion_of_core_inact_bitmask(i,1))
det_tmp(i,2) = xor(det_in(i,2),reunion_of_core_inact_bitmask(i,2))
enddo
!call debug_det(det_tmp,nint)
do i=1,Nint
n_closed_shell_cas += popcnt(iand(det_tmp(i,1),det_tmp(i,2)))
enddo
end
subroutine doubly_occ_empty_in_couple(det_in,n_couples,couples,couples_out)
implicit none
use bitmasks
integer, intent(in) :: n_couples,couples(n_couples,2)
integer(bit_kind),intent(in) :: det_in(N_int,2)
logical, intent(out) :: couples_out(0:n_couples)
integer(bit_kind) :: det_tmp(N_int)
integer(bit_kind) :: det_tmp_bis(N_int)
BEGIN_DOC
! n_couples is the number of couples of orbitals to be checked
! couples(i,1) = first orbital of the ith couple
! couples(i,2) = second orbital of the ith couple
! returns the array couples_out
! couples_out(i) = .True. if det_in contains
! an orbital empty in the ith couple AND
! an orbital doubly occupied in the ith couple
END_DOC
integer :: i,j,k,l
! det_tmp tells you if the orbitals are occupied or not
do j = 1, N_int
det_tmp(j) = ior(det_in(j,1),det_in(j,2))
enddo
couples_out(0) = .False.
do i = 1, n_couples
do j = 1, N_int
det_tmp_bis(j) = 0_bit_kind
enddo
call set_bit_to_integer(couples(i,1),det_tmp_bis,N_int) ! first orb
call set_bit_to_integer(couples(i,2),det_tmp_bis,N_int) ! second orb
! det_tmp is zero except for the two orbitals of the couple
integer :: i_count
i_count = 0
do j = 1, N_int
i_count += popcnt(iand(det_tmp(j),det_tmp_bis(j))) ! check if the two orbitals are both occupied
enddo
if(i_count .ne. 1)then
couples_out(i) = .False.
cycle
endif
! test if orbital there are two electrons or not
i_count = 0
do j = 1, N_int
i_count += popcnt(iand(iand(det_in(j,1),det_in(j,2)),det_tmp_bis(j)))
enddo
if(i_count.ne.1)then
couples_out(i) = .False.
else
couples_out(i) = .True.
couples_out(0) = .True.
endif
enddo
end
subroutine give_index_of_doubly_occ_in_active_space(det_in,doubly_occupied_array)
implicit none
use bitmasks
integer(bit_kind), intent(in) :: det_in(N_int,2)
logical, intent(out) :: doubly_occupied_array(n_act_orb)
integer(bit_kind) :: det_tmp(N_int)
integer(bit_kind) :: det_tmp_bis(N_int)
BEGIN_DOC
END_DOC
integer :: i,j,k,l
! det_tmp tells you if the orbitals are occupied or not
do j = 1, N_int
det_tmp(j) = ior(det_in(j,1),det_in(j,2))
enddo
do i = 1, n_act_orb
do j = 1, N_int
det_tmp_bis(j) = 0_bit_kind
enddo
i_bite = list_act(i)
call set_bit_to_integer(i_bite,det_tmp_bis,N_int) ! act orb
! det_tmp is zero except for the orbital "ith" active orbital
integer :: i_count,i_bite
! test if orbital there are two electrons or not
i_count = 0
do j = 1, N_int
i_count += popcnt(iand(iand(det_in(j,1),det_in(j,2)),det_tmp_bis(j)))
enddo
if(i_count.ne.1)then
doubly_occupied_array(i) = .False.
else
doubly_occupied_array(i) = .True.
endif
enddo
end
subroutine doubly_occ_empty_in_couple_and_no_hund_elsewhere(det_in,n_couple_no_hund,couple_ion,couple_no_hund,is_ok)
implicit none
use bitmasks
integer, intent(in) :: n_couple_no_hund,couple_ion(2),couple_no_hund(n_couple_no_hund,2)
integer(bit_kind),intent(in) :: det_in(N_int,2)
logical, intent(out) :: is_ok
integer(bit_kind) :: det_tmp(N_int)
integer(bit_kind) :: det_tmp_bis(N_int)
BEGIN_DOC
! n_couples is the number of couples of orbitals to be checked
! couples(i,1) = first orbital of the ith couple
! couples(i,2) = second orbital of the ith couple
! returns the array couples_out
! couples_out(i) = .True. if det_in contains
! an orbital empty in the ith couple AND
! an orbital doubly occupied in the ith couple
END_DOC
integer :: i,j,k,l
! det_tmp tells you if the orbitals are occupied or not
do j = 1, N_int
det_tmp(j) = ior(det_in(j,1),det_in(j,2))
enddo
is_ok = .False.
do j = 1, N_int
det_tmp_bis(j) = 0_bit_kind
enddo
call set_bit_to_integer(couple_ion(1),det_tmp_bis,N_int) ! first orb
call set_bit_to_integer(couple_ion(2),det_tmp_bis,N_int) ! second orb
! det_tmp is zero except for the two orbitals of the couple
integer :: i_count
i_count = 0
do j = 1, N_int
i_count += popcnt(iand(det_tmp(j),det_tmp_bis(j))) ! check if the two orbitals are both occupied
enddo
if(i_count .ne. 1)then
is_ok = .False.
return
endif
! test if orbital there are two electrons or not
i_count = 0
do j = 1, N_int
i_count += popcnt(iand(iand(det_in(j,1),det_in(j,2)),det_tmp_bis(j)))
enddo
if(i_count.ne.1)then
is_ok = .False.
return
else
do i = 1, n_couple_no_hund
do j = 1, N_int
det_tmp_bis(j) = 0_bit_kind
enddo
call set_bit_to_integer(couple_no_hund (i,1),det_tmp_bis,N_int) ! first orb
call set_bit_to_integer(couple_no_hund (i,2),det_tmp_bis,N_int) ! second orb
! det_tmp_bis is zero except for the two orbitals of the couple
i_count = 0
do j = 1, N_int
i_count += popcnt(iand(det_tmp(j),det_tmp_bis(j))) ! check if the two orbitals are both occupied
enddo
if(i_count .ne. 2)then
is_ok = .False.
return
endif
! test if orbital there are one alpha and one beta
integer :: i_count_alpha,i_count_beta
i_count_alpha = 0
i_count_beta = 0
do j = 1, N_int
i_count_alpha += popcnt(iand(det_in(j,1),det_tmp_bis(j)))
i_count_beta += popcnt(iand(det_in(j,2),det_tmp_bis(j)))
enddo
if(i_count_alpha==1.and.i_count_beta==1)then
is_ok = .True.
else
is_ok = .False.
return
endif
enddo
is_ok = .True.
endif
end
subroutine neutral_no_hund_in_couple(det_in,n_couples,couples,couples_out)
implicit none
use bitmasks
integer, intent(in) :: n_couples,couples(n_couples,2)
integer(bit_kind),intent(in) :: det_in(N_int,2)
logical, intent(out) :: couples_out(0:n_couples)
integer(bit_kind) :: det_tmp(N_int)
integer(bit_kind) :: det_tmp_bis(N_int)
BEGIN_DOC
! n_couples is the number of couples of orbitals to be checked
! couples(i,1) = first orbital of the ith couple
! couples(i,2) = second orbital of the ith couple
! returns the array couples_out
! couples_out(i) = .True. if det_in contains
! an orbital empty in the ith couple AND
! an orbital doubly occupied in the ith couple
END_DOC
integer :: i,j,k,l
! det_tmp tells you if the orbitals are occupied or not
do j = 1, N_int
det_tmp(j) = ior(det_in(j,1),det_in(j,2))
enddo
couples_out(0) = .True.
do i = 1, n_couples
do j = 1, N_int
det_tmp_bis(j) = 0_bit_kind
enddo
call set_bit_to_integer(couples(i,1),det_tmp_bis,N_int) ! first orb
call set_bit_to_integer(couples(i,2),det_tmp_bis,N_int) ! second orb
! det_tmp_bis is zero except for the two orbitals of the couple
integer :: i_count
i_count = 0
do j = 1, N_int
i_count += popcnt(iand(det_tmp(j),det_tmp_bis(j))) ! check if the two orbitals are both occupied
enddo
if(i_count .ne. 2)then
couples_out(i) = .False.
cycle
endif
! test if orbital there are one alpha and one beta
integer :: i_count_alpha,i_count_beta
i_count_alpha = 0
i_count_beta = 0
do j = 1, N_int
i_count_alpha += popcnt(iand(det_in(j,1),det_tmp_bis(j)))
i_count_beta += popcnt(iand(det_in(j,2),det_tmp_bis(j)))
enddo
if(i_count_alpha==1.and.i_count_beta==1)then
couples_out(i) = .True.
else
couples_out(i) = .False.
endif
enddo
do i = 1, n_couples
if(.not.couples_out(i))then
couples_out(0) = .False.
endif
enddo
end

1
src/Hartree_Fock/SCF.irp.f
View File

@ -53,7 +53,6 @@ subroutine run
! Choose SCF algorithm
! call damping_SCF ! Deprecated routine
call Roothaan_Hall_SCF
end

131
src/Hartree_Fock/damping_SCF.irp.f
View File

@ -1,131 +0,0 @@
subroutine damping_SCF
implicit none
double precision :: E
double precision, allocatable :: D_alpha(:,:), D_beta(:,:)
double precision :: E_new
double precision, allocatable :: D_new_alpha(:,:), D_new_beta(:,:), F_new(:,:)
double precision, allocatable :: delta_alpha(:,:), delta_beta(:,:)
double precision :: lambda, E_half, a, b, delta_D, delta_E, E_min
integer :: i,j,k
logical :: saving
character :: save_char
allocate( &
D_alpha( ao_num, ao_num ), &
D_beta( ao_num, ao_num ), &
F_new( ao_num, ao_num ), &
D_new_alpha( ao_num, ao_num ), &
D_new_beta( ao_num, ao_num ), &
delta_alpha( ao_num, ao_num ), &
delta_beta( ao_num, ao_num ))
do j=1,ao_num
do i=1,ao_num
D_alpha(i,j) = HF_density_matrix_ao_alpha(i,j)
D_beta (i,j) = HF_density_matrix_ao_beta (i,j)
enddo
enddo
call write_time(6)
write(6,'(A4,1X,A16, 1X, A16, 1X, A16, 1X, A4 )') &
'====','================','================','================', '===='
write(6,'(A4,1X,A16, 1X, A16, 1X, A16, 1X, A4 )') &
' N ', 'Energy ', 'Energy diff ', 'Density diff ', 'Save'
write(6,'(A4,1X,A16, 1X, A16, 1X, A16, 1X, A4 )') &
'====','================','================','================', '===='
E = HF_energy + 1.d0
E_min = HF_energy
delta_D = 0.d0
do k=1,n_it_scf_max
delta_E = HF_energy - E
E = HF_energy
if ( (delta_E < 0.d0).and.(dabs(delta_E) < thresh_scf) ) then
exit
endif
saving = E < E_min
if (saving) then
call save_mos
save_char = 'X'
E_min = E
else
save_char = ' '
endif
write(6,'(I4,1X,F16.10, 1X, F16.10, 1X, F16.10, 3X, A )') &
k, E, delta_E, delta_D, save_char
D_alpha = HF_density_matrix_ao_alpha
D_beta = HF_density_matrix_ao_beta
mo_coef = eigenvectors_fock_matrix_mo
TOUCH mo_coef
D_new_alpha = HF_density_matrix_ao_alpha
D_new_beta = HF_density_matrix_ao_beta
F_new = Fock_matrix_ao
E_new = HF_energy
delta_alpha = D_new_alpha - D_alpha
delta_beta = D_new_beta - D_beta
lambda = .5d0
E_half = 0.d0
do while (E_half > E)
HF_density_matrix_ao_alpha = D_alpha + lambda * delta_alpha
HF_density_matrix_ao_beta = D_beta + lambda * delta_beta
TOUCH HF_density_matrix_ao_alpha HF_density_matrix_ao_beta
mo_coef = eigenvectors_fock_matrix_mo
TOUCH mo_coef
E_half = HF_energy
if ((E_half > E).and.(E_new < E)) then
lambda = 1.d0
exit
else if ((E_half > E).and.(lambda > 5.d-4)) then
lambda = 0.5d0 * lambda
E_new = E_half
else
exit
endif
enddo
a = (E_new + E - 2.d0*E_half)*2.d0
b = -E_new - 3.d0*E + 4.d0*E_half
lambda = -lambda*b/(a+1.d-16)
D_alpha = (1.d0-lambda) * D_alpha + lambda * D_new_alpha
D_beta = (1.d0-lambda) * D_beta + lambda * D_new_beta
delta_E = HF_energy - E
do j=1,ao_num
do i=1,ao_num
delta_D = delta_D + &
(D_alpha(i,j) - HF_density_matrix_ao_alpha(i,j))*(D_alpha(i,j) - HF_density_matrix_ao_alpha(i,j)) + &
(D_beta (i,j) - HF_density_matrix_ao_beta (i,j))*(D_beta (i,j) - HF_density_matrix_ao_beta (i,j))
enddo
enddo
delta_D = dsqrt(delta_D/dble(ao_num)**2)
HF_density_matrix_ao_alpha = D_alpha
HF_density_matrix_ao_beta = D_beta
TOUCH HF_density_matrix_ao_alpha HF_density_matrix_ao_beta
mo_coef = eigenvectors_fock_matrix_mo
TOUCH mo_coef
enddo
write(6,'(A4,1X,A16, 1X, A16, 1X, A16, 1X, A4 )') '====','================','================','================', '===='
write(6,*)
if(.not.no_oa_or_av_opt)then
call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo,size(Fock_matrix_mo,1),size(Fock_matrix_mo,2),mo_label,1,.true.)
endif
call write_double(6, E_min, 'Hartree-Fock energy')
call ezfio_set_hartree_fock_energy(E_min)
call write_time(6)
deallocate(D_alpha,D_beta,F_new,D_new_alpha,D_new_beta,delta_alpha,delta_beta)
end

46
src/Hartree_Fock/diagonalize_fock.irp.f
View File

@ -1,17 +1,19 @@
BEGIN_PROVIDER [ double precision, diagonal_Fock_matrix_mo, (mo_tot_num) ]
&BEGIN_PROVIDER [ double precision, eigenvectors_Fock_matrix_mo, (ao_num,mo_tot_num) ]
BEGIN_PROVIDER [ double precision, eigenvectors_Fock_matrix_mo, (ao_num,mo_tot_num) ]
implicit none
BEGIN_DOC
! Diagonal Fock matrix in the MO basis
! Eigenvector of the Fock matrix in the MO basis obtained with level shift.
END_DOC
integer :: i,j
integer :: liwork, lwork, n, info
integer, allocatable :: iwork(:)
double precision, allocatable :: work(:), F(:,:), S(:,:)
double precision, allocatable :: diag(:)
allocate( F(mo_tot_num,mo_tot_num) )
allocate (diag(mo_tot_num) )
do j=1,mo_tot_num
do i=1,mo_tot_num
F(i,j) = Fock_matrix_mo(i,j)
@ -40,8 +42,6 @@
endif
! Insert level shift here
do i = elec_beta_num+1, elec_alpha_num
F(i,i) += 0.5d0*level_shift
@ -62,8 +62,7 @@
liwork = -1
call dsyevd( 'V', 'U', mo_tot_num, F, &
size(F,1), diagonal_Fock_matrix_mo, &
work, lwork, iwork, liwork, info)
size(F,1), diag, work, lwork, iwork, liwork, info)
if (info /= 0) then
print *, irp_here//' DSYEVD failed : ', info
@ -77,15 +76,13 @@
allocate(work(lwork))
allocate(iwork(liwork) )
call dsyevd( 'V', 'U', mo_tot_num, F, &
size(F,1), diagonal_Fock_matrix_mo, &
work, lwork, iwork, liwork, info)
size(F,1), diag, work, lwork, iwork, liwork, info)
deallocate(iwork)
if (info /= 0) then
call dsyev( 'V', 'L', mo_tot_num, F, &
size(F,1), diagonal_Fock_matrix_mo, &
work, lwork, info)
size(F,1), diag, work, lwork, info)
if (info /= 0) then
print *, irp_here//' DSYEV failed : ', info
@ -96,33 +93,8 @@
call dgemm('N','N',ao_num,mo_tot_num,mo_tot_num, 1.d0, &
mo_coef, size(mo_coef,1), F, size(F,1), &
0.d0, eigenvectors_Fock_matrix_mo, size(eigenvectors_Fock_matrix_mo,1))
deallocate(work, F)
deallocate(work, F, diag)
END_PROVIDER
BEGIN_PROVIDER [double precision, diagonal_Fock_matrix_mo_sum, (mo_tot_num)]
implicit none
BEGIN_DOC
! diagonal element of the fock matrix calculated as the sum over all the interactions
! with all the electrons in the RHF determinant
! diagonal_Fock_matrix_mo_sum(i) = sum_{j=1, N_elec} 2 J_ij -K_ij
END_DOC
integer :: i,j
double precision :: accu
do j = 1,elec_alpha_num
accu = 0.d0
do i = 1, elec_alpha_num
accu += 2.d0 * mo_bielec_integral_jj_from_ao(i,j) - mo_bielec_integral_jj_exchange_from_ao(i,j)
enddo
diagonal_Fock_matrix_mo_sum(j) = accu + mo_mono_elec_integral(j,j)
enddo
do j = elec_alpha_num+1,mo_tot_num
accu = 0.d0
do i = 1, elec_alpha_num
accu += 2.d0 * mo_bielec_integral_jj_from_ao(i,j) - mo_bielec_integral_jj_exchange_from_ao(i,j)
enddo
diagonal_Fock_matrix_mo_sum(j) = accu + mo_mono_elec_integral(j,j)
enddo
END_PROVIDER

56
src/Integrals_Bielec/map_integrals.irp.f
View File

@ -642,62 +642,6 @@ subroutine dump_$ao_integrals(filename)
end
IRP_IF COARRAY
subroutine communicate_$ao_integrals()
use map_module
implicit none
BEGIN_DOC
! Communicate the $ao integrals with co-array
END_DOC
integer(cache_key_kind), pointer :: key(:)
real(integral_kind), pointer :: val(:)
integer*8 :: i,j, k, nmax
integer*8, save :: n[*]
integer :: copy_n
real(integral_kind), allocatable :: buffer_val(:)[:]
integer(cache_key_kind), allocatable :: buffer_key(:)[:]
real(integral_kind), allocatable :: copy_val(:)
integer(key_kind), allocatable :: copy_key(:)
n = 0_8
do i=0_8,$ao_integrals_map%map_size
n = max(n,$ao_integrals_map%map(i)%n_elements)
enddo
sync all
nmax = 0_8
do j=1,num_images()
nmax = max(nmax,n[j])
enddo
allocate( buffer_key(nmax)[*], buffer_val(nmax)[*])
allocate( copy_key(nmax), copy_val(nmax))
do i=0_8,$ao_integrals_map%map_size
key => $ao_integrals_map%map(i)%key
val => $ao_integrals_map%map(i)%value
n = $ao_integrals_map%map(i)%n_elements
do j=1,n
buffer_key(j) = key(j)
buffer_val(j) = val(j)
enddo
sync all
do j=1,num_images()
if (j /= this_image()) then
copy_n = n[j]
do k=1,copy_n
copy_val(k) = buffer_val(k)[j]
copy_key(k) = buffer_key(k)[j]
copy_key(k) = copy_key(k)+shiftl(i,map_shift)
enddo
call map_append($ao_integrals_map, copy_key, copy_val, copy_n )
endif
enddo
sync all
enddo
deallocate( buffer_key, buffer_val, copy_val, copy_key)
end
IRP_ENDIF
integer function load_$ao_integrals(filename)
implicit none

49
src/Integrals_Bielec/mo_bi_integrals.irp.f
View File

@ -157,49 +157,6 @@ BEGIN_PROVIDER [ logical, mo_bielec_integrals_in_map ]
END_PROVIDER
subroutine set_integrals_jj_into_map
use bitmasks
implicit none
integer :: i,j,n_integrals,i0,j0
double precision :: buffer_value(mo_tot_num * mo_tot_num)
integer(key_kind) :: buffer_i(mo_tot_num*mo_tot_num)
n_integrals = 0
do j0 = 1, n_virt_orb
j = list_virt(j0)
do i0 = j0, n_virt_orb
i = list_virt(i0)
n_integrals += 1
! mo_bielec_integral_jj_exchange(i,j) = mo_bielec_integral_vv_exchange_from_ao(i,j)
call mo_bielec_integrals_index(i,j,i,j,buffer_i(n_integrals))
buffer_value(n_integrals) = mo_bielec_integral_vv_from_ao(i,j)
enddo
enddo
call insert_into_mo_integrals_map(n_integrals,buffer_i,buffer_value,&
real(mo_integrals_threshold,integral_kind))
call map_merge(mo_integrals_map)
end
subroutine set_integrals_exchange_jj_into_map
use bitmasks
implicit none
integer :: i,j,n_integrals,i0,j0
double precision :: buffer_value(mo_tot_num * mo_tot_num)
integer(key_kind) :: buffer_i(mo_tot_num*mo_tot_num)
n_integrals = 0
do j0 = 1, n_virt_orb
j = list_virt(j0)
do i0 = j0+1, n_virt_orb
i = list_virt(i0)
n_integrals += 1
call mo_bielec_integrals_index(i,j,j,i,buffer_i(n_integrals))
buffer_value(n_integrals) = mo_bielec_integral_vv_exchange_from_ao(i,j)
enddo
enddo
call insert_into_mo_integrals_map(n_integrals,buffer_i,buffer_value,&
real(mo_integrals_threshold,integral_kind))
call map_merge(mo_integrals_map)
end
subroutine add_integrals_to_map(mask_ijkl)
use bitmasks
@ -1386,9 +1343,3 @@ subroutine clear_mo_map
end
subroutine provide_all_mo_integrals
implicit none
provide mo_integrals_map mo_bielec_integral_jj mo_bielec_integral_jj_anti
provide mo_bielec_integral_jj_exchange mo_bielec_integrals_in_map
end

91
src/MO_Basis/mos.irp.f
View File

@ -242,65 +242,42 @@ subroutine ao_to_mo(A_ao,LDA_ao,A_mo,LDA_mo)
deallocate(T)
end
subroutine mo_to_ao_no_overlap(A_mo,LDA_mo,A_ao,LDA_ao)
implicit none
BEGIN_DOC
! Transform A from the MO basis to the S^-1 AO basis
! Useful for density matrix
END_DOC
integer, intent(in) :: LDA_ao,LDA_mo
double precision, intent(in) :: A_mo(LDA_mo,mo_tot_num)
double precision, intent(out) :: A_ao(LDA_ao,ao_num)
double precision, allocatable :: T(:,:)
allocate ( T(mo_tot_num,ao_num) )
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T
call dgemm('N','T', mo_tot_num, ao_num, mo_tot_num, &
1.d0, A_mo,size(A_mo,1), &
mo_coef, size(mo_coef,1), &
0.d0, T, size(T,1))
call dgemm('N','N', ao_num, ao_num, mo_tot_num, &
1.d0, mo_coef,size(mo_coef,1), &
T, size(T,1), &
0.d0, A_ao, size(A_ao,1))
deallocate(T)
end
subroutine mix_mo_jk(j,k)
implicit none
integer, intent(in) :: j,k
integer :: i,i_plus,i_minus
BEGIN_DOC
! subroutine that rotates the jth MO with the kth MO
! to give two new MO's that are
! '+' = 1/sqrt(2) (|j> + |k>)
! '-' = 1/sqrt(2) (|j> - |k>)
! by convention, the '+' MO is in the lower index (min(j,k))
! by convention, the '-' MO is in the greater index (max(j,k))
END_DOC
double precision :: array_tmp(ao_num,2),dsqrt_2
if(j==k)then
print*,'You want to mix two orbitals that are the same !'
print*,'It does not make sense ... '
print*,'Stopping ...'
stop
endif
array_tmp = 0.d0
dsqrt_2 = 1.d0/dsqrt(2.d0)
do i = 1, ao_num
array_tmp(i,1) = dsqrt_2 * (mo_coef(i,j) + mo_coef(i,k))
array_tmp(i,2) = dsqrt_2 * (mo_coef(i,j) - mo_coef(i,k))
enddo
i_plus = min(j,k)
i_minus = max(j,k)
do i = 1, ao_num
mo_coef(i,i_plus) = array_tmp(i,1)
mo_coef(i,i_minus) = array_tmp(i,2)
enddo
implicit none
integer, intent(in) :: j,k
integer :: i,i_plus,i_minus
BEGIN_DOC
! Rotates the jth MO with the kth MO
! to give two new MO's that are
!
! '+' = 1/sqrt(2) (|j> + |k>)
!
! '-' = 1/sqrt(2) (|j> - |k>)
!
! by convention, the '+' MO is in the lower index (min(j,k))
! by convention, the '-' MO is in the larger index (max(j,k))
END_DOC
double precision :: array_tmp(ao_num,2),dsqrt_2
if(j==k)then
print*,'You want to mix two orbitals that are the same !'
print*,'It does not make sense ... '
print*,'Stopping ...'
stop
endif
array_tmp = 0.d0
dsqrt_2 = 1.d0/dsqrt(2.d0)
do i = 1, ao_num
array_tmp(i,1) = dsqrt_2 * (mo_coef(i,j) + mo_coef(i,k))
array_tmp(i,2) = dsqrt_2 * (mo_coef(i,j) - mo_coef(i,k))
enddo
i_plus = min(j,k)
i_minus = max(j,k)
do i = 1, ao_num
mo_coef(i,i_plus) = array_tmp(i,1)
mo_coef(i,i_minus) = array_tmp(i,2)
enddo
end
subroutine ao_ortho_cano_to_ao(A_ao,LDA_ao,A,LDA)

143
src/MO_Basis/utils.irp.f
View File

@ -216,147 +216,4 @@ subroutine mo_as_svd_vectors_of_mo_matrix_eig(matrix,lda,m,n,eig,label)
end
subroutine mo_as_eigvectors_of_mo_matrix_sort_by_observable(matrix,observable,n,m,label)
implicit none
integer,intent(in) :: n,m
character*(64), intent(in) :: label
double precision, intent(in) :: matrix(n,m),observable(n,n)
double precision, allocatable :: mo_coef_new(:,:), R(:,:),eigvalues(:),value(:)
integer,allocatable :: iorder(:)
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: mo_coef_new, R
call write_time(6)
if (m /= mo_tot_num) then
print *, irp_here, ': Error : m/= mo_tot_num'
stop 1
endif
allocate(R(n,m),mo_coef_new(ao_num,m),eigvalues(m),value(m),iorder(m))
mo_coef_new = mo_coef
call lapack_diag(eigvalues,R,matrix,size(matrix,1),size(matrix,2))
do i = 1, mo_tot_num
value(i) = 0.d0
iorder(i) = i
do j = 1, mo_tot_num
do k = 1, mo_tot_num
value(i) += R(k,i) * R(j,i) * observable(k,j)
enddo
enddo
! print*,'value(i) = ',i,value(i)
enddo
integer :: i,j,k,index
double precision :: R_tmp(m,m),obs(m)
print*,'sort ....'
call dsort(value,iorder,n)
do i = 1, mo_tot_num
index = iorder(i)
! print*,'iorder(i) = ',iorder(i)
obs(i) = eigvalues(iorder(i))
do j = 1, mo_tot_num
R_tmp(j,i) = R(j,index)
enddo
enddo
do i = 1, mo_tot_num
do j = 1, mo_tot_num
R(j,i) = R_tmp(j,i)
enddo
enddo
do i = 1, mo_tot_num
value(i) = 0.d0
do j = 1, mo_tot_num
do k = 1, mo_tot_num
value(i) += R(k,i) * R(j,i) * observable(k,j)
enddo
enddo
print*,'i = ',i
print*,'value(i) = ',value(i)
print*,'obs(i) = ',obs(i)
print*,''
enddo
write (6,'(A)') 'MOs are now **'//trim(label)//'**'
write (6,'(A)') ''
write (6,'(A)') 'Eigenvalues'
write (6,'(A)') '-----------'
write (6,'(A)') ''
write (6,'(A)') '======== ================'
do i = 1, m
write (6,'(I8,1X,F16.10)') i,eigvalues(i)
enddo
write (6,'(A)') '======== ================'
write (6,'(A)') ''
call dgemm('N','N',ao_num,m,m,1.d0,mo_coef_new,size(mo_coef_new,1),R,size(R,1),0.d0,mo_coef,size(mo_coef,1))
deallocate(mo_coef_new,R,eigvalues)
call write_time(6)
mo_label = label
SOFT_TOUCH mo_coef mo_label
end
subroutine mo_sort_by_observable(observable,label)
implicit none
character*(64), intent(in) :: label
double precision, intent(in) :: observable(mo_tot_num)
double precision, allocatable :: mo_coef_new(:,:),value(:)
integer,allocatable :: iorder(:)
allocate(mo_coef_new(ao_num,mo_tot_num),value(mo_tot_num),iorder(mo_tot_num))
print*,'allocate !'
mo_coef_new = mo_coef
do i = 1, mo_tot_num
iorder(i) = i
value(i) = observable(i)
enddo
integer :: i,j,k,index
print*,'sort ....'
call dsort(value,iorder,mo_tot_num)
do i = 1, mo_tot_num
index = iorder(i)
do j = 1, mo_tot_num
mo_coef(j,i) = mo_coef_new(j,index)
enddo
enddo
write (6,'(A)') 'MOs are now **'//trim(label)//'**'
write (6,'(A)') ''
deallocate(mo_coef_new,value)
! call write_time(6)
mo_label = label
SOFT_TOUCH mo_coef mo_label
end
subroutine give_all_mos_at_r(r,mos_array)
implicit none
double precision, intent(in) :: r(3)
double precision, intent(out) :: mos_array(mo_tot_num)
call give_specific_mos_at_r(r,mos_array, mo_coef)
end
subroutine give_specific_mos_at_r(r,mos_array, mo_coef_specific)
implicit none
double precision, intent(in) :: r(3)
double precision, intent(in) :: mo_coef_specific(ao_num, mo_tot_num)
double precision, intent(out) :: mos_array(mo_tot_num)
double precision :: aos_array(ao_num),accu
integer :: i,j
call give_all_aos_at_r(r,aos_array)
do i = 1, mo_tot_num
accu = 0.d0
do j = 1, ao_num
accu += mo_coef_specific(j,i) * aos_array(j)
enddo
mos_array(i) = accu
enddo
end

193
src/Perturbation/pt2_equations.irp.f
View File

@ -356,199 +356,6 @@ end
subroutine pt2_epstein_nesbet_SC2_projected ($arguments)
use bitmasks
implicit none
$declarations
BEGIN_DOC
! compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
!
! for the various N_st states,
!
! but with the correction in the denominator
!
! comming from the interaction of that determinant with all the others determinants
!
! that can be repeated by repeating all the double excitations
!
! : you repeat all the correlation energy already taken into account in electronic_energy(1)
!
! that could be repeated to this determinant.
!
! In addition, for the perturbative energetic contribution you have the standard second order
!
! e_2_pert = <psi_i|H|det_pert>^2/(Delta_E)
!
! and also the purely projected contribution
!
! H_pert_diag = <HF|H|det_pert> c_pert
END_DOC
double precision :: i_H_psi_array(N_st)
integer :: idx_repeat(0:ndet)
integer :: i,j,degree,l
double precision :: diag_H_mat_elem_fock,accu_e_corr,hij,h0j,h,delta_E
double precision :: repeat_all_e_corr,accu_e_corr_tmp,e_2_pert_fonda
ASSERT (Nint == N_int)
ASSERT (Nint > 0)
call i_H_psi_SC2(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array,idx_repeat)
accu_e_corr = 0.d0
!$IVDEP
do i = 1, idx_repeat(0)
accu_e_corr = accu_e_corr + E_corr_per_selectors(idx_repeat(i))
enddo
h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint)
h = h + accu_e_corr
delta_E = 1.d0/(electronic_energy(1) - h)
c_pert(1) = i_H_psi_array(1) /(electronic_energy(1) - h)
e_2_pert(1) = i_H_psi_array(1) * c_pert(1)
do i =2,N_st
H_pert_diag(i) = h
if (dabs(electronic_energy(i) - h) > 1.d-6) then
c_pert(i) = i_H_psi_array(i) / (-dabs(electronic_energy(i) - h))
e_2_pert(i) = (c_pert(i) * i_H_psi_array(i))
else
c_pert(i) = i_H_psi_array(i)
e_2_pert(i) = -dabs(i_H_psi_array(i))
endif
enddo
degree = popcnt(xor( ref_bitmask(1,1), det_pert(1,1))) + &
popcnt(xor( ref_bitmask(1,2), det_pert(1,2)))
!DIR$ NOUNROLL
do l=2,Nint
degree = degree+ popcnt(xor( ref_bitmask(l,1), det_pert(l,1))) + &
popcnt(xor( ref_bitmask(l,2), det_pert(l,2)))
enddo
if(degree==4)then
! <psi|delta_H|psi>
e_2_pert_fonda = e_2_pert(1)
H_pert_diag(1) = e_2_pert(1) * c_pert(1) * c_pert(1)
do i = 1, N_st
do j = 1, idx_repeat(0)
e_2_pert(i) += e_2_pert_fonda * psi_selectors_coef(idx_repeat(j),i) * psi_selectors_coef(idx_repeat(j),i)
enddo
enddo
endif
end
subroutine pt2_epstein_nesbet_SC2_no_projected ($arguments)
use bitmasks
implicit none
$declarations
BEGIN_DOC
! compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
!
! for the various N_st states,
!
! but with the correction in the denominator
!
! comming from the interaction of that determinant with all the others determinants
!
! that can be repeated by repeating all the double excitations
!
! : you repeat all the correlation energy already taken into account in electronic_energy(1)
!
! that could be repeated to this determinant.
!
! In addition, for the perturbative energetic contribution you have the standard second order
!
! e_2_pert = <psi_i|H|det_pert>^2/(Delta_E)
!
! and also the purely projected contribution
!
! H_pert_diag = <HF|H|det_pert> c_pert
END_DOC
double precision :: i_H_psi_array(N_st)
integer :: idx_repeat(0:ndet)
integer :: i,j,degree,l
double precision :: diag_H_mat_elem_fock,accu_e_corr,hij,h0j,h,delta_E
double precision :: repeat_all_e_corr,accu_e_corr_tmp,e_2_pert_fonda
ASSERT (Nint == N_int)
ASSERT (Nint > 0)
call i_H_psi_SC2(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array,idx_repeat)
accu_e_corr = 0.d0
!$IVDEP
do i = 1, idx_repeat(0)
accu_e_corr = accu_e_corr + E_corr_per_selectors(idx_repeat(i))
enddo
h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint)
h = h + accu_e_corr
delta_E = 1.d0/(electronic_energy(1) - h)
c_pert(1) = i_H_psi_array(1) /(electronic_energy(1) - h)
e_2_pert(1) = i_H_psi_array(1) * c_pert(1)
do i =2,N_st
H_pert_diag(i) = h
if (dabs(electronic_energy(i) - h) > 1.d-6) then
c_pert(i) = i_H_psi_array(i) / (-dabs(electronic_energy(i) - h))
e_2_pert(i) = (c_pert(i) * i_H_psi_array(i))
else
c_pert(i) = i_H_psi_array(i)
e_2_pert(i) = -dabs(i_H_psi_array(i))
endif
enddo
end
subroutine pt2_epstein_nesbet_sc2 ($arguments)
use bitmasks
implicit none
$declarations
BEGIN_DOC
! compute the standard Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
!
! for the various N_st states, but with the CISD_SC2 energies and coefficients
!
! c_pert(i) = <psi(i)|H|det_pert>/( E(i) - <det_pert|H|det_pert> )
!
! e_2_pert(i) = <psi(i)|H|det_pert>^2/( E(i) - <det_pert|H|det_pert> )
!
END_DOC
integer :: i,j
double precision :: i_H_psi_array(N_st)
double precision :: diag_H_mat_elem_fock, h
PROVIDE selection_criterion
ASSERT (Nint == N_int)
ASSERT (Nint > 0)
!call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array)
call i_H_psi_minilist(det_pert,minilist,idx_minilist,N_minilist,psi_selectors_coef,Nint,N_minilist,psi_selectors_size,N_st,i_H_psi_array)
h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint)
do i =1,N_st
if(electronic_energy(i)>h.and.electronic_energy(i).ne.0.d0)then
c_pert(i) = -1.d0
e_2_pert(i) = selection_criterion*selection_criterion_factor*2.d0
else if (dabs(electronic_energy(i) - h) > 1.d-6) then
c_pert(i) = i_H_psi_array(i) / (electronic_energy(i) - h)
H_pert_diag(i) = h*c_pert(i)*c_pert(i)
e_2_pert(i) = c_pert(i) * i_H_psi_array(i)
else
c_pert(i) = -1.d0
e_2_pert(i) = -dabs(i_H_psi_array(i))
H_pert_diag(i) = h
endif
enddo
end
subroutine pt2_dummy ($arguments)
use bitmasks

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

@ -143,12 +143,12 @@ subroutine ortho_qr(A,LDA,m,n)
allocate (jpvt(n), tau(n), work(1))
LWORK=-1
call dgeqrf( m, n, A, LDA, TAU, WORK, LWORK, INFO )
call dgeqrf( m, n, A, LDA, TAU, WORK, LWORK, INFO )
LWORK=2*int(WORK(1))
deallocate(WORK)
allocate(WORK(LWORK))
call dgeqrf(m, n, A, LDA, TAU, WORK, LWORK, INFO )
call dorgqr(m, n, n, A, LDA, tau, WORK, LWORK, INFO)
call dgeqrf(m, n, A, LDA, TAU, WORK, LWORK, INFO )
call dorgqr(m, n, n, A, LDA, tau, WORK, LWORK, INFO)
deallocate(WORK,jpvt,tau)
end
@ -174,7 +174,7 @@ subroutine ortho_qr_unblocked(A,LDA,m,n)
double precision, allocatable :: tau(:), work(:)
allocate (jpvt(n), tau(n), work(n))
call dgeqr2( m, n, A, LDA, TAU, WORK, INFO )
call dgeqr2( m, n, A, LDA, TAU, WORK, INFO )
call dorg2r(m, n, n, A, LDA, tau, WORK, INFO)
deallocate(WORK,jpvt,tau)
end
@ -397,20 +397,20 @@ subroutine lapack_diagd(eigvalues,eigvectors,H,nmax,n)
integer ,allocatable :: iwork(:)
double precision,allocatable :: A(:,:)
integer :: lwork, info, i,j,l,k, liwork
allocate(A(nmax,n),eigenvalues(n))
! print*,'Diagonalization by jacobi'
! print*,'n = ',n
! print*,'Diagonalization by jacobi'
! print*,'n = ',n
A=H
lwork = max(1000,2*n*n + 6*n+ 1)
liwork = max(5*n + 3,1000)
allocate (work(lwork),iwork(liwork))
lwork = -1
liwork = -1
call DSYEVD( 'V', 'U', n, A, nmax, eigenvalues, work, lwork, &
iwork, liwork, info )
call DSYEVD( 'V', 'U', n, A, nmax, eigenvalues, work, lwork, &
iwork, liwork, info )
if (info < 0) then
print *, irp_here, ': DSYEVD: the ',-info,'-th argument had an illegal value'
stop 2
@ -418,20 +418,20 @@ subroutine lapack_diagd(eigvalues,eigvectors,H,nmax,n)
lwork = int( work( 1 ) )
liwork = iwork(1)
deallocate (work,iwork)
allocate (work(lwork),iwork(liwork))
call DSYEVD( 'V', 'U', n, A, nmax, eigenvalues, work, lwork, &
iwork, liwork, info )
call DSYEVD( 'V', 'U', n, A, nmax, eigenvalues, work, lwork, &
iwork, liwork, info )
deallocate(work,iwork)
if (info < 0) then
print *, irp_here, ': DSYEVD: the ',-info,'-th argument had an illegal value'
stop 2
else if( info > 0 ) then
write(*,*)'DSYEVD Failed'
stop 1
write(*,*)'DSYEVD Failed'
stop 1
end if
eigvectors = 0.d0
eigvalues = 0.d0
do j = 1, n
@ -465,8 +465,6 @@ subroutine lapack_diag(eigvalues,eigvectors,H,nmax,n)
integer :: lwork, info, i,j,l,k, liwork
allocate(A(nmax,n),eigenvalues(n))
! print*,'Diagonalization by jacobi'
! print*,'n = ',n
A=H
lwork = 2*n*n + 6*n+ 1
@ -511,165 +509,3 @@ subroutine lapack_diag(eigvalues,eigvectors,H,nmax,n)
deallocate(A,eigenvalues)
end
subroutine lapack_diag_s2(eigvalues,eigvectors,H,nmax,n)
implicit none
BEGIN_DOC
! Diagonalize matrix H
!
! H is untouched between input and ouptut
!
! eigevalues(i) = ith lowest eigenvalue of the H matrix
!
! eigvectors(i,j) = <i|psi_j> where i is the basis function and psi_j is the j th eigenvector
!
END_DOC
integer, intent(in) :: n,nmax
double precision, intent(out) :: eigvectors(nmax,n)
double precision, intent(out) :: eigvalues(n)
double precision, intent(in) :: H(nmax,n)
double precision,allocatable :: eigenvalues(:)
double precision,allocatable :: work(:)
double precision,allocatable :: A(:,:)
integer :: lwork, info, i,j,l,k, liwork
allocate(A(nmax,n),eigenvalues(n))
! print*,'Diagonalization by jacobi'
! print*,'n = ',n
A=H
lwork = 2*n*n + 6*n+ 1
allocate (work(lwork))
lwork = -1
call DSYEV( 'V', 'U', n, A, nmax, eigenvalues, work, lwork, &
info )
if (info < 0) then
print *, irp_here, ': DSYEV: the ',-info,'-th argument had an illegal value'
stop 2
endif
lwork = int( work( 1 ) )
deallocate (work)
allocate (work(lwork))
call DSYEV( 'V', 'U', n, A, nmax, eigenvalues, work, lwork, &
info )
deallocate(work)
if (info < 0) then
print *, irp_here, ': DSYEV: the ',-info,'-th argument had an illegal value'
stop 2
else if( info > 0 ) then
write(*,*)'DSYEV Failed'
stop 1
end if
eigvectors = 0.d0
eigvalues = 0.d0
do j = 1, n
eigvalues(j) = eigenvalues(j)
do i = 1, n
eigvectors(i,j) = A(i,j)
enddo
enddo
deallocate(A,eigenvalues)
end
subroutine lapack_partial_diag(eigvalues,eigvectors,H,nmax,n,n_st)
implicit none
BEGIN_DOC
! Diagonalize matrix H
!
! H is untouched between input and ouptut
!
! eigevalues(i) = ith lowest eigenvalue of the H matrix
!
! eigvectors(i,j) = <i|psi_j> where i is the basis function and psi_j is the j th eigenvector
!
END_DOC
integer, intent(in) :: n,nmax,n_st
double precision, intent(out) :: eigvectors(nmax,n)
double precision, intent(out) :: eigvalues(n)
double precision, intent(in) :: H(nmax,n)
double precision,allocatable :: work(:)
integer ,allocatable :: iwork(:), isuppz(:)
double precision,allocatable :: A(:,:)
integer :: lwork, info, i,j,l,k,m, liwork
allocate(A(nmax,n))
A=H
lwork = 2*n*n + 6*n+ 1
liwork = 5*n + 3
allocate (work(lwork),iwork(liwork),isuppz(2*N_st))
lwork = -1
liwork = -1
call DSYEVR( 'V', 'I', 'U', n, A, nmax, 0.d0, 0.d0, 1, n_st, 1.d-10, m, eigvalues, eigvectors, nmax, isuppz, work, lwork, &
iwork, liwork, info )
if (info < 0) then
print *, irp_here, ': DSYEVR: the ',-info,'-th argument had an illegal value'
stop 2
endif
lwork = int( work( 1 ) )
liwork = iwork(1)
deallocate (work,iwork)
allocate (work(lwork),iwork(liwork))
call DSYEVR( 'V', 'I', 'U', n, A, nmax, 0.d0, 0.d0, 1, n_st, 1.d-10, m, eigvalues, eigvectors, nmax, isuppz, work, lwork, &
iwork, liwork, info )
deallocate(work,iwork)
if (info < 0) then
print *, irp_here, ': DSYEVR: the ',-info,'-th argument had an illegal value'
stop 2
else if( info > 0 ) then
write(*,*)'DSYEVR Failed'
stop 1
end if
deallocate(A)
end
subroutine set_zero_extra_diag(i1,i2,matrix,lda,m)
implicit none
integer, intent(in) :: i1,i2,lda,m
double precision, intent(inout) :: matrix(lda,m)
integer :: i,j
do j=i1,i2
do i = 1,i1-1
matrix(i,j) = 0.d0
matrix(j,i) = 0.d0
enddo
enddo
do i = i2,i1
do j=i2+1,m
matrix(i,j) = 0.d0
matrix(j,i) = 0.d0
enddo
enddo
end
subroutine matrix_vector_product(u0,u1,matrix,sze,lda)
implicit none
BEGIN_DOC
! performs u1 += u0 * matrix
END_DOC
integer, intent(in) :: sze,lda
double precision, intent(in) :: u0(sze)
double precision, intent(inout) :: u1(sze)
double precision, intent(in) :: matrix(lda,sze)
integer :: i,j
integer :: incx,incy
incx = 1
incy = 1
call dsymv('U', sze, 1.d0, matrix, lda, u0, incx, 1.d0, u1, incy)
end

78
src/Utils/one_e_integration.irp.f
View File

@ -32,84 +32,6 @@ double precision function overlap_gaussian_x(A_center,B_center,alpha,beta,power_
end
subroutine overlap_A_B_C(dim,alpha,beta,gama,a,b,A_center,B_center,Nucl_center,overlap)
implicit none
include 'constants.include.F'
integer, intent(in) :: dim
integer, intent(in) :: a(3),b(3) ! powers : (x-xa)**a_x = (x-A(1))**a(1)
double precision, intent(in) :: alpha, beta, gama ! exponents
double precision, intent(in) :: A_center(3) ! A center
double precision, intent(in) :: B_center (3) ! B center
double precision, intent(in) :: Nucl_center(3) ! B center
double precision, intent(out) :: overlap
double precision :: P_new(0:max_dim,3),P_center(3),fact_p,p
double precision :: F_integral_tab(0:max_dim)
integer :: iorder_p(3)
double precision :: overlap_x,overlap_z,overlap_y
call give_explicit_poly_and_gaussian_double(P_new,P_center,p,fact_p,iorder_p,alpha,beta,gama,a,b,A_center,B_center,Nucl_center,dim)
if(fact_p.lt.1d-10)then
! overlap_x = 0.d0
! overlap_y = 0.d0
! overlap_z = 0.d0
overlap = 0.d0
return
endif
integer :: nmax
double precision :: F_integral
nmax = maxval(iorder_p)
do i = 0,nmax
F_integral_tab(i) = F_integral(i,p)
enddo
overlap_x = P_new(0,1) * F_integral_tab(0)
overlap_y = P_new(0,2) * F_integral_tab(0)
overlap_z = P_new(0,3) * F_integral_tab(0)
integer :: i
do i = 1,iorder_p(1)
overlap_x += P_new(i,1) * F_integral_tab(i)
enddo
do i = 1,iorder_p(2)
overlap_y += P_new(i,2) * F_integral_tab(i)
enddo
do i = 1,iorder_p(3)
overlap_z += P_new(i,3) * F_integral_tab(i)
enddo
overlap = overlap_x * overlap_y * overlap_z * fact_p
!double precision :: overlap_x_1,overlap_y_1,overlap_z_1,overlap_1
!call test(alpha,beta,gama,a,b,A_center,B_center,Nucl_center,overlap_x_1,overlap_y_1,overlap_z_1,overlap_1)
!print*,'overlap_1 = ',overlap_1
!print*,'overlap = ',overlap
!if(dabs(overlap - overlap_1).ge.1.d-3)then
! print*,'power_A(1) = ',a(1)
! print*,'power_A(2) = ',a(2)
! print*,'power_A(3) = ',a(3)
! print*,'power_B(1) = ',b(1)
! print*,'power_B(2) = ',b(2)
! print*,'power_B(3) = ',b(3)
! print*,'alpha = ',alpha
! print*,'beta = ',beta
! print*,'gama = ',gama
! print*,'A_center(1) = ',A_center(1)
! print*,'A_center(2) = ',A_center(2)
! print*,'A_center(3) = ',A_center(3)
! print*,'B_center(1) = ',B_center(1)
! print*,'B_center(2) = ',B_center(2)
! print*,'B_center(3) = ',B_center(3)
! print*,'Nucl_center(1) = ',Nucl_center(1)
! print*,'Nucl_center(2) = ',Nucl_center(2)
! print*,'Nucl_center(3) = ',Nucl_center(3)
! print*,'overlap = ',overlap
! print*,'overlap_1=',overlap_1
! stop
!endif
end
subroutine overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,&
power_B,overlap_x,overlap_y,overlap_z,overlap,dim)
implicit none