LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-11-12 17:13:54 +01:00
Code Issues Releases Wiki Activity

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

Browse Source
This commit is contained in:
Anthony Scemama 2014-10-08 12:00:45 +02:00
commit 45d9331cd7
16 changed files with 237 additions and 59 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
data/ezfio_defaults
View File

@ -17,8 +17,9 @@ cis_dressed
determinants
n_states 1
n_states_diag determinants_n_states
n_det_max_jacobi 5000
threshold_generators 0.999
threshold_generators 0.995
threshold_selectors 0.999
read_wf False

4
setup_environment.sh
View File

@ -45,8 +45,8 @@ make -C ocaml qp_create_ezfio_from_xyz.native
if [[ $? -ne 0 ]]
then
scripts/fetch_from_web.py "https://raw.github.com/hcarty/ocamlbrew/master/ocamlbrew-install" ocamlbrew-install.sh
cat < ocamlbrew-install.sh | env OCAMLBREW_FLAGS="-r" bash > ocaml_install.log
grep "source " install.log | grep "etc/ocamlbrew.bashrc" >> quantum_package.rc
cat < ocamlbrew-install.sh | env OCAMLBREW_FLAGS="-r" bash | tee ocaml_install.log
grep "source " ocaml_install.log | grep "etc/ocamlbrew.bashrc" >> quantum_package.rc
source quantum_package.rc
echo Y | opam install core async
fi

2
src/CISD_SC2_selected/H_apply.irp.f
View File

@ -4,7 +4,7 @@ from generate_h_apply import *
from perturbation import perturbations
s = H_apply("PT2",SingleRef=True)
s.set_perturbation("epstein_nesbet_sc2_projected")
s.set_perturbation("epstein_nesbet_sc2_no_projected")
print s
s = H_apply("PT2_en_sc2",SingleRef=True)

14
src/CISD_SC2_selected/cisd_sc2_selection.irp.f
View File

@ -12,8 +12,9 @@ program cisd_sc2_selected
pt2 = 1.d0
perturbation = "epstein_nesbet_sc2_projected"
E_old(1) = HF_energy
davidson_threshold = 1.d-8
davidson_threshold = 1.d-10
if (N_det > n_det_max_cisd_sc2) then
call diagonalize_CI_SC2
call save_wavefunction
@ -31,6 +32,8 @@ program cisd_sc2_selected
print *, '-----'
endif
integer :: i_count
i_count = 0
do while (N_det < n_det_max_cisd_sc2.and.maxval(abs(pt2(1:N_st))) > pt2_max)
print*,'----'
print*,''
@ -49,6 +52,13 @@ program cisd_sc2_selected
E_old(i) = CI_SC2_energy(i)
enddo
! print *, 'E corr = ', (E_old(1)) - HF_energy
if(dabs(E_old(i) - CI_SC2_energy(i) ).le.1.d-12)then
i_count += 1
selection_criterion_factor = selection_criterion_factor * 0.5d0
if(i_count > 5)then
exit
endif
endif
if (abort_all) then
exit
endif
@ -81,7 +91,7 @@ program cisd_sc2_selected
print *, 'PT2(SC2) = ', pt2(i)
print *, 'E(SC2) = ', CI_SC2_energy(i)
print *, 'E_before(SC2)+PT2(SC2) = ', CI_SC2_energy(i)+pt2(i)
print *, 'E_before(SC2)+PT2(SC2)_new = ', CI_SC2_energy(i)+pt2(i)*(1.d0+norm_pert)
print *, 'E_before(SC2)+PT2(SC2)_new = ', CI_SC2_energy(i)+pt2(i)* (1.d0 + norm_pert) - H_pert_diag(i)
print*,'greater coeficient of the state : ',dabs(psi_coef(imax,i))
call get_excitation_degree(ref_bitmask,psi_det(1,1,imax),degree,N_int)

2
src/Dets/SC2.irp.f
View File

@ -175,7 +175,7 @@ subroutine CISD_SC2(dets_in,u_in,energies,dim_in,sze,N_st,Nint,convergence)
enddo
call lapack_diag(eigenvalues,eigenvectors, &
H_matrix_tmp,size(H_matrix_all_dets,1),sze)
do j=1,min(N_states,sze)
do j=1,min(N_states_diag,sze)
do i=1,sze
u_in(i,j) = eigenvectors(i,j)
enddo

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

@ -12,7 +12,7 @@ BEGIN_PROVIDER [ integer, davidson_sze_max ]
! Max number of Davidson sizes
END_DOC
ASSERT (davidson_sze_max <= davidson_iter_max)
davidson_sze_max = 8*N_states
davidson_sze_max = 8*N_states_diag
END_PROVIDER
subroutine davidson_diag(dets_in,u_in,energies,dim_in,sze,N_st,Nint,iunit)

2
src/Dets/determinants.ezfio_config
View File

@ -4,6 +4,7 @@ determinants
mo_label character*(64)
n_det integer
n_states integer
n_states_diag integer
psi_coef double precision (determinants_n_det,determinants_n_states)
psi_det integer*8 (determinants_n_int*determinants_bit_kind/8,2,determinants_n_det)
n_det_max_jacobi integer
@ -13,4 +14,5 @@ determinants
det_occ integer (electrons_elec_alpha_num,determinants_det_num,2)
det_coef double precision (determinants_det_num)
read_wf logical
expected_s2 double precision

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

@ -150,7 +150,7 @@ subroutine read_dets(det,Nint,Ndet)
end
BEGIN_PROVIDER [ double precision, psi_coef, (psi_det_size,N_states) ]
BEGIN_PROVIDER [ double precision, psi_coef, (psi_det_size,N_states_diag) ]
implicit none
BEGIN_DOC
! The wave function coefficients. Initialized with Hartree-Fock if the EZFIO file
@ -161,6 +161,11 @@ BEGIN_PROVIDER [ double precision, psi_coef, (psi_det_size,N_states) ]
logical :: exists
double precision, allocatable :: psi_coef_read(:,:)
character*(64) :: label
psi_coef = 0.d0
do i=1,N_states_diag
psi_coef(i,i) = 1.d0
enddo
if (read_wf) then
call ezfio_has_determinants_psi_coef(exists)
@ -183,22 +188,8 @@ BEGIN_PROVIDER [ double precision, psi_coef, (psi_det_size,N_states) ]
enddo
deallocate(psi_coef_read)
else
psi_coef = 0.d0
do i=1,N_states
psi_coef(i,i) = 1.d0
enddo
endif
else
psi_coef = 0.d0
do i=1,N_states
psi_coef(i,i) = 1.d0
enddo
endif

40
src/Dets/diagonalize_CI.irp.f
View File

@ -9,13 +9,13 @@ BEGIN_PROVIDER [ character*(64), diag_algorithm ]
diag_algorithm = "Lapack"
endif
if (N_det < N_states) then
if (N_det < N_states_diag) then
diag_algorithm = "Lapack"
endif
END_PROVIDER
BEGIN_PROVIDER [ double precision, CI_energy, (N_states) ]
BEGIN_PROVIDER [ double precision, CI_energy, (N_states_diag) ]
implicit none
BEGIN_DOC
! N_states lowest eigenvalues of the CI matrix
@ -24,23 +24,25 @@ BEGIN_PROVIDER [ double precision, CI_energy, (N_states) ]
integer :: j
character*(8) :: st
call write_time(output_Dets)
do j=1,N_states
do j=1,N_states_diag
CI_energy(j) = CI_electronic_energy(j) + nuclear_repulsion
write(st,'(I4)') j
call write_double(output_Dets,CI_energy(j),'Energy of state '//trim(st))
call write_double(output_Dets,CI_eigenvectors_s2(j),'S^2 of state '//trim(st))
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, CI_electronic_energy, (N_states) ]
&BEGIN_PROVIDER [ double precision, CI_eigenvectors, (N_det,N_states) ]
BEGIN_PROVIDER [ double precision, CI_electronic_energy, (N_states_diag) ]
&BEGIN_PROVIDER [ double precision, CI_eigenvectors, (N_det,N_states_diag) ]
&BEGIN_PROVIDER [ double precision, CI_eigenvectors_s2, (N_states_diag) ]
implicit none
BEGIN_DOC
! Eigenvectors/values of the CI matrix
END_DOC
integer :: i,j
do j=1,N_states
do j=1,N_states_diag
do i=1,N_det
CI_eigenvectors(i,j) = psi_coef(i,j)
enddo
@ -49,7 +51,7 @@ END_PROVIDER
if (diag_algorithm == "Davidson") then
call davidson_diag(psi_det,CI_eigenvectors,CI_electronic_energy, &
size(CI_eigenvectors,1),N_det,N_states,N_int,output_Dets)
size(CI_eigenvectors,1),N_det,N_states_diag,N_int,output_Dets)
else if (diag_algorithm == "Lapack") then
@ -59,11 +61,21 @@ END_PROVIDER
call lapack_diag(eigenvalues,eigenvectors, &
H_matrix_all_dets,size(H_matrix_all_dets,1),N_det)
CI_electronic_energy(:) = 0.d0
do j=1,min(N_states,N_det)
do i=1,N_det
CI_eigenvectors(i,j) = eigenvectors(i,j)
enddo
CI_electronic_energy(j) = eigenvalues(j)
integer :: i_state
double precision :: s2
j=0
i_state = 0
do while(i_state.lt.min(N_states_diag,N_det))
j+=1
call get_s2_u0(psi_det,eigenvectors(1,j),N_det,N_det,s2)
if(dabs(s2-expected_s2).le.0.1d0)then
i_state += 1
do i=1,N_det
CI_eigenvectors(i,i_state) = eigenvectors(i,j)
enddo
CI_electronic_energy(i_state) = eigenvalues(j)
CI_eigenvectors_s2(i_state) = s2
endif
enddo
deallocate(eigenvectors,eigenvalues)
endif
@ -77,10 +89,10 @@ subroutine diagonalize_CI
! eigenstates of the CI matrix
END_DOC
integer :: i,j
do j=1,N_states
do j=1,N_states_diag
do i=1,N_det
psi_coef(i,j) = CI_eigenvectors(i,j)
enddo
enddo
SOFT_TOUCH psi_coef CI_electronic_energy CI_energy CI_eigenvectors
SOFT_TOUCH psi_coef CI_electronic_energy CI_energy CI_eigenvectors CI_eigenvectors_s2
end

24
src/Dets/diagonalize_CI_SC2.irp.f
View File

@ -1,13 +1,13 @@
BEGIN_PROVIDER [ double precision, CI_SC2_energy, (N_states) ]
BEGIN_PROVIDER [ double precision, CI_SC2_energy, (N_states_diag) ]
implicit none
BEGIN_DOC
! N_states lowest eigenvalues of the CI matrix
! N_states_diag lowest eigenvalues of the CI matrix
END_DOC
integer :: j
character*(8) :: st
call write_time(output_Dets)
do j=1,N_states
do j=1,N_states_diag
CI_SC2_energy(j) = CI_SC2_electronic_energy(j) + nuclear_repulsion
write(st,'(I4)') j
call write_double(output_Dets,CI_SC2_energy(j),'Energy of state '//trim(st))
@ -20,38 +20,38 @@ END_PROVIDER
BEGIN_DOC
! convergence of the correlation energy of SC2 iterations
END_DOC
threshold_convergence_SC2 = 1.d-8
threshold_convergence_SC2 = 1.d-10
END_PROVIDER
BEGIN_PROVIDER [ double precision, CI_SC2_electronic_energy, (N_states) ]
&BEGIN_PROVIDER [ double precision, CI_SC2_eigenvectors, (N_det,N_states) ]
BEGIN_PROVIDER [ double precision, CI_SC2_electronic_energy, (N_states_diag) ]
&BEGIN_PROVIDER [ double precision, CI_SC2_eigenvectors, (N_det,N_states_diag) ]
implicit none
BEGIN_DOC
! Eigenvectors/values of the CI matrix
END_DOC
integer :: i,j
do j=1,N_states
do j=1,N_states_diag
do i=1,N_det
! CI_SC2_eigenvectors(i,j) = psi_coef(i,j)
CI_SC2_eigenvectors(i,j) = CI_eigenvectors(i,j)
CI_SC2_eigenvectors(i,j) = psi_coef(i,j)
! CI_SC2_eigenvectors(i,j) = CI_eigenvectors(i,j)
enddo
CI_SC2_electronic_energy(j) = CI_electronic_energy(j)
enddo
double precision :: convergence
call CISD_SC2(psi_det,CI_SC2_eigenvectors,CI_SC2_electronic_energy, &
size(CI_SC2_eigenvectors,1),N_det,N_states,N_int,threshold_convergence_SC2)
size(CI_SC2_eigenvectors,1),N_det,N_states_diag,N_int,threshold_convergence_SC2)
END_PROVIDER
subroutine diagonalize_CI_SC2
implicit none
BEGIN_DOC
! Replace the coefficients of the CI states by the coefficients of the
! Replace the coefficients of the CI states_diag by the coefficients of the
! eigenstates of the CI matrix
END_DOC
integer :: i,j
do j=1,N_states
do j=1,N_states_diag
do i=1,N_det
psi_coef(i,j) = CI_SC2_eigenvectors(i,j)
enddo

11
src/Dets/options.irp.f
View File

@ -9,6 +9,16 @@ T.set_ezfio_name( "N_states" )
T.set_output ( "output_dets" )
print T
T.set_type ( "integer" )
T.set_name ( "N_states_diag" )
T.set_doc ( "Number of states to consider for the diagonalization " )
T.set_ezfio_dir ( "determinants" )
T.set_ezfio_name( "N_states_diag" )
T.set_output ( "output_dets" )
print T
T.set_name ( "N_det_max_jacobi" )
T.set_doc ( "Maximum number of determinants diagonalized by Jacobi" )
T.set_ezfio_name( "N_det_max_jacobi" )
@ -20,6 +30,7 @@ T.set_doc ( "If true, read the wave function from the EZFIO file" )
T.set_ezfio_name( "read_wf" )
T.set_output ( "output_dets" )
print T
END_SHELL

15
src/Dets/s2.irp.f
View File

@ -42,6 +42,21 @@ BEGIN_PROVIDER [ double precision, S_z ]
END_PROVIDER
BEGIN_PROVIDER [ double precision, expected_s2]
implicit none
PROVIDE ezfio_filename
logical :: has_expected_s2
call ezfio_has_determinants_expected_s2(has_expected_s2)
if (has_expected_s2) then
call ezfio_get_determinants_expected_s2(expected_s2)
else
expected_s2 = elec_alpha_num - elec_beta_num + 0.5d0 * ((elec_alpha_num - elec_beta_num)**2*0.5d0 - (elec_alpha_num-elec_beta_num))
! call ezfio_set_determinants_expected_s2(expected_s2)
endif
END_PROVIDER
subroutine get_s2_u0(psi_keys_tmp,psi_coefs_tmp,n,nmax,s2)
implicit none

74
src/Perturbation/pert_sc2.irp.f
View File

@ -72,8 +72,8 @@ subroutine pt2_epstein_nesbet_SC2_projected(det_pert,c_pert,e_2_pert,H_pert_diag
enddo
if(degree==4)then
! <psi|delta_H|psi>
H_pert_diag(1) = e_2_pert(1)
e_2_pert_fonda = H_pert_diag(1)
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)
@ -83,6 +83,76 @@ subroutine pt2_epstein_nesbet_SC2_projected(det_pert,c_pert,e_2_pert,H_pert_diag
end
subroutine pt2_epstein_nesbet_SC2_no_projected(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st)
use bitmasks
implicit none
integer, intent(in) :: Nint,ndet,N_st
integer(bit_kind), intent(in) :: det_pert(Nint,2)
double precision , intent(out) :: c_pert(N_st),e_2_pert(N_st),H_pert_diag(N_st)
double precision :: i_H_psi_array(N_st)
integer :: idx_repeat(0:ndet)
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 CI_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
integer :: i,j,degree,l
double precision :: diag_H_mat_elem,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(det_pert,Nint) + accu_e_corr
delta_E = 1.d0/(CI_SC2_electronic_energy(1) - h)
c_pert(1) = i_H_psi_array(1) /(CI_SC2_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(CI_SC2_electronic_energy(i) - h) > 1.d-6) then
c_pert(i) = i_H_psi_array(i) / (-dabs(CI_SC2_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
double precision function repeat_all_e_corr(key_in)
implicit none
integer(bit_kind), intent(in) :: key_in(N_int,2)

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

@ -307,6 +307,72 @@ 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

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

@ -15,10 +15,10 @@ double precision function overlap_gaussian_x(A_center,B_center,alpha,beta,power_
call give_explicit_poly_and_gaussian_x(P_new,P_center,p,fact_p,iorder_p,alpha,&
beta,power_A,power_B,A_center,B_center,dim)
if(fact_p.lt.0.000001d0)then
overlap_gaussian_x = 0.d0
return
endif
! if(fact_p.lt.0.000001d0)then
! overlap_gaussian_x = 0.d0
! return
! endif
overlap_gaussian_x = 0.d0
integer :: i

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

@ -49,17 +49,17 @@ double precision function binom_func(i,j)
end
BEGIN_PROVIDER [ double precision, binom, (0:20,0:20) ]
&BEGIN_PROVIDER [ double precision, binom_transp, (0:20,0:20) ]
BEGIN_PROVIDER [ double precision, binom, (0:40,0:40) ]
&BEGIN_PROVIDER [ double precision, binom_transp, (0:40,0:40) ]
implicit none
BEGIN_DOC
! Binomial coefficients
END_DOC
integer :: k,l
double precision :: fact, f
do k=0,20
do k=0,40
f = fact(k)
do l=0,20
do l=0,40
binom(k,l) = f/(fact(l)*fact(k-l))
binom_transp(l,k) = binom(k,l)
enddo