mirror of
https://github.com/LCPQ/quantum_package
synced 2024-10-02 06:20:59 +02:00
MRCC works
This commit is contained in:
parent
102bbb0b4f
commit
99e63935d4
@ -1,32 +0,0 @@
|
||||
BEGIN_SHELL [ /usr/bin/python ]
|
||||
from ezfio_with_default import EZFIO_Provider
|
||||
T = EZFIO_Provider()
|
||||
T.set_type ( "integer" )
|
||||
T.set_name ( "N_det_max_fci" )
|
||||
T.set_doc ( "Max number of determinants in the wave function" )
|
||||
T.set_ezfio_dir ( "full_ci" )
|
||||
T.set_ezfio_name( "N_det_max_fci" )
|
||||
T.set_output ( "output_full_ci" )
|
||||
print T
|
||||
|
||||
T.set_type ( "logical" )
|
||||
T.set_name ( "do_pt2_end" )
|
||||
T.set_doc ( "If true, compute the PT2 at the end of the selection" )
|
||||
T.set_ezfio_name( "do_pt2_end" )
|
||||
print T
|
||||
|
||||
T.set_type ( "double precision" )
|
||||
T.set_name ( "pt2_max" )
|
||||
T.set_doc ( """The selection process stops when the largest PT2 (for all the states)
|
||||
is lower than pt2_max in absolute value""" )
|
||||
T.set_ezfio_name( "pt2_max" )
|
||||
print T
|
||||
|
||||
T.set_type ( "double precision" )
|
||||
T.set_name ( "var_pt2_ratio" )
|
||||
T.set_doc ( """The selection process stops when the energy ratio variational/(variational+PT2)
|
||||
is equal to var_pt2_ratio""" )
|
||||
T.set_ezfio_name( "var_pt2_ratio" )
|
||||
print T
|
||||
END_SHELL
|
||||
|
||||
@ -1,14 +1,82 @@
|
||||
-4.142795384334731
|
||||
E(CAS+SD) = -3.93510180989509
|
||||
|
||||
-0.726935163332
|
||||
++++----
|
||||
++++----
|
||||
|
||||
0.685482292841
|
||||
+++-+---
|
||||
+++-+---
|
||||
|
||||
0.0409791961003
|
||||
++-+-+--
|
||||
++-+-+--
|
||||
|
||||
c2*c3/c1 = -0.038642391672
|
||||
|
||||
|
||||
-0.996413146877
|
||||
++++----
|
||||
++++----
|
||||
|
||||
-0.0598366139154
|
||||
++-+-+--
|
||||
+++-+---
|
||||
|
||||
-0.0598366139154
|
||||
+++-+---
|
||||
++-+-+--
|
||||
|
||||
|
||||
=======
|
||||
|
||||
0.706616635698
|
||||
+++-+---
|
||||
+++-+---
|
||||
|
||||
-0.692594178414
|
||||
++++----
|
||||
++++----
|
||||
|
||||
-0.0915716740265
|
||||
++-+-+--
|
||||
+++-+---
|
||||
|
||||
-0.0915716740265
|
||||
+++-+---
|
||||
++-+-+--
|
||||
|
||||
0.0461418323107
|
||||
++-+-+--
|
||||
++-+-+--
|
||||
|
||||
-0.0458957789452
|
||||
++--++--
|
||||
++--++--
|
||||
|
||||
|
||||
----
|
||||
|
||||
0.713102867186
|
||||
++++----
|
||||
++++----
|
||||
|
||||
-0.688067688244
|
||||
+++-+---
|
||||
+++-+---
|
||||
|
||||
0.089262694488
|
||||
+++-+---
|
||||
++-+-+--
|
||||
|
||||
0.089262694488
|
||||
++-+-+--
|
||||
+++-+---
|
||||
|
||||
-0.0459510603899
|
||||
++-+-+--
|
||||
++-+-+--
|
||||
|
||||
4.695183071437694E-002
|
||||
Determinant 64
|
||||
---------------------------------------------
|
||||
000000000000002E|000000000000002E
|
||||
|-+++-+----------------------------------------------------------|
|
||||
|-+++-+----------------------------------------------------------|
|
||||
|
||||
CAS-SD: -4.14214374069306
|
||||
: -4.14230904320551
|
||||
|
||||
E0 = -11.5634986758976
|
||||
|
||||
|
||||
@ -3,6 +3,8 @@ program mrcc
|
||||
read_wf = .True.
|
||||
TOUCH read_wf
|
||||
call run
|
||||
call run_mrcc
|
||||
! call run_mrcc_test
|
||||
end
|
||||
|
||||
subroutine run
|
||||
@ -12,8 +14,7 @@ subroutine run
|
||||
print *, N_det_cas
|
||||
print *, N_det_sd
|
||||
|
||||
! call update_generators
|
||||
integer :: i
|
||||
integer :: i,j
|
||||
print *, 'CAS'
|
||||
print *, '==='
|
||||
do i=1,N_det_cas
|
||||
@ -21,23 +22,54 @@ subroutine run
|
||||
call debug_det(psi_cas(1,1,i),N_int)
|
||||
enddo
|
||||
|
||||
! print *, 'SD'
|
||||
! print *, '=='
|
||||
! do i=1,N_det_sd
|
||||
! print *, psi_sd_coefs(i,:)
|
||||
! call debug_det(psi_sd(1,1,i),N_int)
|
||||
! enddo
|
||||
!
|
||||
print *, 'SD'
|
||||
print *, '=='
|
||||
do i=1,N_det_sd
|
||||
print *, psi_sd_coefs(i,:)
|
||||
call debug_det(psi_sd(1,1,i),N_int)
|
||||
enddo
|
||||
print *, 'xxx', 'Energy CAS+SD', ci_energy
|
||||
end
|
||||
subroutine run_mrcc_test
|
||||
implicit none
|
||||
integer :: i,j
|
||||
double precision :: pt2
|
||||
pt2 = 0.d0
|
||||
do j=1,N_det
|
||||
do i=1,N_det
|
||||
pt2 += psi_coef(i,1)*psi_coef(j,1) * delta_ij(i,j,1)
|
||||
enddo
|
||||
enddo
|
||||
print *, ci_energy(1)
|
||||
print *, ci_energy(1)+pt2
|
||||
end
|
||||
subroutine run_mrcc
|
||||
implicit none
|
||||
integer :: i,j
|
||||
print *, 'MRCC'
|
||||
print *, '===='
|
||||
print *, ci_energy(:)
|
||||
print *, h_matrix_all_dets(3,3), delta_ij(3,3,1)
|
||||
print *, h_matrix_all_dets(3,3), delta_ij(3,3,1)
|
||||
print *, ci_energy_dressed(:)
|
||||
! print *, 'max', maxval(delta_ij_sd)
|
||||
! print *, 'min', minval(delta_ij_sd)
|
||||
!
|
||||
! do i=1,N_det
|
||||
! print '(10(F10.6,X))', delta_ij(i,1:N_det,1)
|
||||
! enddo
|
||||
print *, ''
|
||||
print *, 'CAS+SD energy : ', ci_energy_dressed(:)
|
||||
print *, ''
|
||||
|
||||
! call diagonalize_ci_dressed
|
||||
! call save_wavefunction_unsorted
|
||||
double precision :: E_new, E_old, delta_e
|
||||
integer :: iteration
|
||||
E_new = 0.d0
|
||||
delta_E = 1.d0
|
||||
iteration = 0
|
||||
do while (delta_E > 1.d-8)
|
||||
iteration += 1
|
||||
print *, '==========================='
|
||||
print *, 'MRCC Iteration', iteration
|
||||
print *, '==========================='
|
||||
print *, ''
|
||||
E_old = sum(ci_energy_dressed)
|
||||
call diagonalize_ci_dressed
|
||||
E_new = sum(ci_energy_dressed)
|
||||
delta_E = dabs(E_new - E_old)
|
||||
call write_double(6,ci_energy_dressed(1),"MRCC energy")
|
||||
enddo
|
||||
|
||||
end
|
||||
|
||||
@ -1,3 +1,120 @@
|
||||
subroutine mrcc_dress(delta_ij_,Ndet,i_generator,n_selected,det_buffer,Nint,iproc)
|
||||
use bitmasks
|
||||
implicit none
|
||||
|
||||
integer, intent(in) :: i_generator,n_selected, Nint, iproc
|
||||
integer, intent(in) :: Ndet
|
||||
double precision, intent(inout) :: delta_ij_(Ndet,Ndet,*)
|
||||
|
||||
integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected)
|
||||
integer :: i,j,k,l,m
|
||||
logical :: is_in_wavefunction
|
||||
integer :: degree_alpha(psi_det_size)
|
||||
integer :: degree_I(psi_det_size)
|
||||
integer :: idx_I(0:psi_det_size)
|
||||
integer :: idx_alpha(0:psi_det_size)
|
||||
logical :: good
|
||||
|
||||
integer(bit_kind) :: tq(Nint,2,n_selected)
|
||||
integer :: N_tq, c_ref ,degree
|
||||
integer :: connected_to_ref
|
||||
|
||||
call find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq)
|
||||
|
||||
double precision :: hIk, hIl, hla, dIk(N_states), dka(N_states), dIa(N_states)
|
||||
double precision :: haj, phase, phase2
|
||||
double precision :: f(N_states), ci_inv(N_states)
|
||||
integer :: exc(0:2,2,2)
|
||||
integer :: h1,h2,p1,p2,s1,s2
|
||||
integer(bit_kind):: tmp_det(Nint,2)
|
||||
integer :: iint, ipos
|
||||
! integer :: istate, i_sd, i_cas
|
||||
|
||||
|
||||
|
||||
! |I>
|
||||
|
||||
! |alpha>
|
||||
do i=1,N_tq
|
||||
call get_excitation_degree_vector(psi_sd,tq(1,1,i),degree_alpha,Nint,N_det_sd,idx_alpha)
|
||||
|
||||
! |I>
|
||||
do j=1,N_det_cas
|
||||
! Find triples and quadruple grand parents
|
||||
call get_excitation_degree(tq(1,1,i),psi_cas(1,1,j),degree,Nint)
|
||||
if (degree > 4) then
|
||||
cycle
|
||||
endif
|
||||
dIa(:) = 0.d0
|
||||
! <I| <> |alpha>
|
||||
do k=1,idx_alpha(0)
|
||||
call get_excitation_degree(psi_cas(1,1,j),psi_sd(1,1,idx_alpha(k)),degree,Nint)
|
||||
if (degree > 2) then
|
||||
cycle
|
||||
endif
|
||||
! <I| k |alpha>
|
||||
! <I|H|k>
|
||||
call i_h_j(psi_cas(1,1,j),psi_sd(1,1,idx_alpha(k)),Nint,hIk)
|
||||
dIk(:) = hIk * lambda_mrcc(idx_alpha(k),:)
|
||||
! Exc(k -> alpha)
|
||||
call get_excitation(psi_sd(1,1,idx_alpha(k)),tq(1,1,i),exc,degree,phase,Nint)
|
||||
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
|
||||
tmp_det(1:Nint,1:2) = psi_cas(1,1,j)
|
||||
! Hole (see list_to_bitstring)
|
||||
iint = ishft(h1-1,-bit_kind_shift) + 1
|
||||
ipos = h1-ishft((iint-1),bit_kind_shift)-1
|
||||
tmp_det(iint,s1) = ibclr(tmp_det(iint,s1),ipos)
|
||||
|
||||
! Particle
|
||||
iint = ishft(p1-1,-bit_kind_shift) + 1
|
||||
ipos = p1-ishft((iint-1),bit_kind_shift)-1
|
||||
tmp_det(iint,s1) = ibset(tmp_det(iint,s1),ipos)
|
||||
if (degree == 2) then
|
||||
! Hole (see list_to_bitstring)
|
||||
iint = ishft(h2-1,-bit_kind_shift) + 1
|
||||
ipos = h2-ishft((iint-1),bit_kind_shift)-1
|
||||
tmp_det(iint,s2) = ibclr(tmp_det(iint,s2),ipos)
|
||||
|
||||
! Particle
|
||||
iint = ishft(p2-1,-bit_kind_shift) + 1
|
||||
ipos = p2-ishft((iint-1),bit_kind_shift)-1
|
||||
tmp_det(iint,s2) = ibset(tmp_det(iint,s2),ipos)
|
||||
endif
|
||||
|
||||
dka(:) = 0.d0
|
||||
do l=k+1,idx_alpha(0)
|
||||
call get_excitation_degree(tmp_det,psi_sd(1,1,idx_alpha(l)),degree,Nint)
|
||||
if (degree == 0) then
|
||||
call get_excitation(psi_cas(1,1,j),psi_sd(1,1,idx_alpha(l)),exc,degree,phase2,Nint)
|
||||
call i_h_j(psi_cas(1,1,j),psi_sd(1,1,idx_alpha(l)),Nint,hIl)
|
||||
dka(:) = hIl * lambda_mrcc(idx_alpha(l),:) * phase * phase2
|
||||
exit
|
||||
endif
|
||||
enddo
|
||||
do l=1,N_states
|
||||
dIa(l) += dka(l)*dIk(l)
|
||||
enddo
|
||||
enddo
|
||||
ci_inv(1:N_states) = 1.d0/psi_cas_coefs(j,1:N_states)
|
||||
do l=1,idx_alpha(0)
|
||||
k = idx_alpha(l)
|
||||
call i_h_j(tq(1,1,i),psi_sd(1,1,idx_alpha(l)),Nint,hla)
|
||||
do m=1,N_states
|
||||
delta_ij_(idx_sd(k),idx_cas(j),m) += dIa(m) * hla
|
||||
delta_ij_(idx_cas(j),idx_sd(k),m) += dIa(m) * hla
|
||||
delta_ij_(idx_cas(j),idx_cas(j),m) -= dIa(m) * hla * ci_inv(m) * psi_sd_coefs(k,m)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
end
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
subroutine mrcc_dress_simple(delta_ij_sd_,Ndet_sd,i_generator,n_selected,det_buffer,Nint,iproc)
|
||||
use bitmasks
|
||||
implicit none
|
||||
@ -18,35 +135,7 @@ subroutine mrcc_dress_simple(delta_ij_sd_,Ndet_sd,i_generator,n_selected,det_buf
|
||||
integer :: N_tq, c_ref
|
||||
integer :: connected_to_ref
|
||||
|
||||
N_tq = 0
|
||||
do i=1,N_selected
|
||||
c_ref = connected_to_ref(det_buffer(1,1,i),psi_det_generators,Nint, &
|
||||
i_generator,N_det_generators)
|
||||
|
||||
if (c_ref /= 0) then
|
||||
cycle
|
||||
endif
|
||||
|
||||
! Select determinants that are triple or quadruple excitations
|
||||
! from the CAS
|
||||
good = .True.
|
||||
call get_excitation_degree_vector(psi_cas,det_buffer(1,1,i),degree,Nint,N_det_cas,idx)
|
||||
do k=1,idx(0)
|
||||
if (degree(k) < 3) then
|
||||
good = .False.
|
||||
exit
|
||||
endif
|
||||
enddo
|
||||
if (good) then
|
||||
if (.not. is_in_wavefunction(det_buffer(1,1,i),Nint,N_det)) then
|
||||
N_tq += 1
|
||||
do k=1,N_int
|
||||
tq(k,1,N_tq) = det_buffer(k,1,i)
|
||||
tq(k,2,N_tq) = det_buffer(k,2,i)
|
||||
enddo
|
||||
endif
|
||||
endif
|
||||
enddo
|
||||
call find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq)
|
||||
|
||||
! Compute <k|H|a><a|H|j> / (E0 - Haa)
|
||||
double precision :: hka, haa
|
||||
@ -73,30 +162,22 @@ subroutine mrcc_dress_simple(delta_ij_sd_,Ndet_sd,i_generator,n_selected,det_buf
|
||||
end
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
subroutine mrcc_dress(delta_ij_sd_,Ndet_sd,i_generator,n_selected,det_buffer,Nint,iproc)
|
||||
subroutine find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq)
|
||||
use bitmasks
|
||||
implicit none
|
||||
|
||||
integer, intent(in) :: i_generator,n_selected, Nint, iproc
|
||||
integer, intent(in) :: Ndet_sd
|
||||
double precision, intent(inout) :: delta_ij_sd_(Ndet_sd,Ndet_sd,*)
|
||||
integer, intent(in) :: i_generator,n_selected, Nint
|
||||
|
||||
integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected)
|
||||
integer :: i,j,k,m
|
||||
integer :: new_size
|
||||
logical :: is_in_wavefunction
|
||||
integer :: degree(psi_det_size)
|
||||
integer :: idx(0:psi_det_size)
|
||||
logical :: good
|
||||
|
||||
integer(bit_kind) :: tq(Nint,2,n_selected)
|
||||
integer :: N_tq, c_ref
|
||||
integer(bit_kind), intent(out) :: tq(Nint,2,n_selected)
|
||||
integer, intent(out) :: N_tq
|
||||
integer :: c_ref
|
||||
integer :: connected_to_ref
|
||||
|
||||
N_tq = 0
|
||||
@ -129,27 +210,11 @@ subroutine mrcc_dress(delta_ij_sd_,Ndet_sd,i_generator,n_selected,det_buffer,Nin
|
||||
endif
|
||||
enddo
|
||||
|
||||
! Compute <k|H|a><a|H|j> / (E0 - Haa)
|
||||
double precision :: hka, haa
|
||||
double precision :: haj
|
||||
double precision :: f(N_states)
|
||||
|
||||
do i=1,N_tq
|
||||
call get_excitation_degree_vector(psi_sd,tq(1,1,i),degree,Nint,Ndet_sd,idx)
|
||||
call i_h_j(tq(1,1,i),tq(1,1,i),Nint,haa)
|
||||
do m=1,N_states
|
||||
f(m) = 1.d0/(ci_electronic_energy(m)-haa)
|
||||
enddo
|
||||
do k=1,idx(0)
|
||||
call i_h_j(tq(1,1,i),psi_sd(1,1,idx(k)),Nint,hka)
|
||||
do j=k,idx(0)
|
||||
call i_h_j(tq(1,1,i),psi_sd(1,1,idx(j)),Nint,haj)
|
||||
do m=1,N_states
|
||||
delta_ij_sd_(idx(k), idx(j),m) += haj*hka* f(m)
|
||||
delta_ij_sd_(idx(j), idx(k),m) += haj*hka* f(m)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
end
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
@ -164,14 +164,7 @@ BEGIN_PROVIDER [ double precision, delta_ij_sd, (N_det_sd, N_det_sd,N_states) ]
|
||||
! Dressing matrix in SD basis
|
||||
END_DOC
|
||||
delta_ij_sd = 0.d0
|
||||
if (dressing_type == "MRCC") then
|
||||
call H_apply_mrcc(delta_ij_sd,N_det_sd)
|
||||
else if (dressing_type == "Simple") then
|
||||
call H_apply_mrcc_simple(delta_ij_sd,N_det_sd)
|
||||
else
|
||||
print *, irp_here
|
||||
stop 'dressing'
|
||||
endif
|
||||
call H_apply_mrcc_simple(delta_ij_sd,N_det_sd)
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ double precision, delta_ij, (N_det,N_det,N_states) ]
|
||||
@ -181,13 +174,17 @@ BEGIN_PROVIDER [ double precision, delta_ij, (N_det,N_det,N_states) ]
|
||||
END_DOC
|
||||
integer :: i,j,m
|
||||
delta_ij = 0.d0
|
||||
do m=1,N_states
|
||||
do j=1,N_det_sd
|
||||
do i=1,N_det_sd
|
||||
delta_ij(idx_sd(i),idx_sd(j),m) = delta_ij_sd(i,j,m)
|
||||
if (dressing_type == "MRCC") then
|
||||
call H_apply_mrcc(delta_ij,N_det)
|
||||
else if (dressing_type == "Simple") then
|
||||
do m=1,N_states
|
||||
do j=1,N_det_sd
|
||||
do i=1,N_det_sd
|
||||
delta_ij(idx_sd(i),idx_sd(j),m) = delta_ij_sd(i,j,m)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
endif
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ double precision, h_matrix_dressed, (N_det,N_det) ]
|
||||
@ -199,9 +196,6 @@ BEGIN_PROVIDER [ double precision, h_matrix_dressed, (N_det,N_det) ]
|
||||
do j=1,N_det
|
||||
do i=1,N_det
|
||||
h_matrix_dressed(i,j) = h_matrix_all_dets(i,j) + delta_ij(i,j,1)
|
||||
if (i==j) then
|
||||
print *, i, delta_ij(i,j,1), h_matrix_all_dets(i,j)
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
|
||||
@ -273,10 +267,22 @@ BEGIN_PROVIDER [ double precision, CI_energy_dressed, (N_states_diag) ]
|
||||
call write_time(output_Dets)
|
||||
do j=1,N_states_diag
|
||||
CI_energy_dressed(j) = CI_electronic_energy_dressed(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
|
||||
|
||||
subroutine diagonalize_CI_dressed
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Replace the coefficients of the CI states by the coefficients of the
|
||||
! eigenstates of the CI matrix
|
||||
END_DOC
|
||||
integer :: i,j
|
||||
do j=1,N_states_diag
|
||||
do i=1,N_det
|
||||
psi_coef(i,j) = CI_eigenvectors_dressed(i,j)
|
||||
enddo
|
||||
enddo
|
||||
SOFT_TOUCH psi_coef
|
||||
|
||||
end
|
||||
|
||||
Loading…
Reference in New Issue
Block a user