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

Moves psi_cas in Dets

Browse Source
This commit is contained in:
Anthony Scemama 2015-04-09 22:32:25 +02:00
parent 99e63935d4
commit 0f2bdedb90
3 changed files with 70 additions and 191 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

27
src/MRCC/mrcc.irp.f
View File

@ -10,10 +10,6 @@ end
subroutine run
implicit none
print *, N_det
print *, N_det_cas
print *, N_det_sd
integer :: i,j
print *, 'CAS'
print *, '==='
@ -22,13 +18,13 @@ 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 *, 'xxx', 'Energy CAS+SD', ci_energy
! print *, 'SD'
! print *, '=='
! do i=1,N_det_non_cas
! print *, psi_non_cas_coefs(i,:)
! call debug_det(psi_non_cas(1,1,i),N_int)
! enddo
call write_double(6,ci_energy(1),"Initial CI energy")
end
subroutine run_mrcc_test
implicit none
@ -46,14 +42,7 @@ end
subroutine run_mrcc
implicit none
integer :: i,j
print *, 'MRCC'
print *, '===='
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
@ -71,5 +60,7 @@ subroutine run_mrcc
delta_E = dabs(E_new - E_old)
call write_double(6,ci_energy_dressed(1),"MRCC energy")
enddo
call ezfio_set_mrcc_energy(ci_energy_dressed(1))
! call save_wavefunction
end

89
src/MRCC/mrcc_dress.irp.f
View File

@ -28,38 +28,44 @@ subroutine mrcc_dress(delta_ij_,Ndet,i_generator,n_selected,det_buffer,Nint,ipro
integer :: h1,h2,p1,p2,s1,s2
integer(bit_kind):: tmp_det(Nint,2)
integer :: iint, ipos
! integer :: istate, i_sd, i_cas
integer :: i_state, k_sd, l_sd, i_I, i_alpha
! |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)
do i_alpha=1,N_tq
call get_excitation_degree_vector(psi_non_cas,tq(1,1,i_alpha),degree_alpha,Nint,N_det_non_cas,idx_alpha)
! |I>
do j=1,N_det_cas
do i_I=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)
call get_excitation_degree(tq(1,1,i_alpha),psi_cas(1,1,i_I),degree,Nint)
if (degree > 4) then
cycle
endif
dIa(:) = 0.d0
do i_state=1,N_states
dIa(i_state) = 0.d0
enddo
! <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)
do k_sd=1,idx_alpha(0)
call get_excitation_degree(psi_cas(1,1,i_I),psi_non_cas(1,1,idx_alpha(k_sd)),degree,Nint)
if (degree > 2) then
cycle
endif
! <I| k |alpha>
! <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 i_h_j(psi_cas(1,1,i_I),psi_non_cas(1,1,idx_alpha(k_sd)),Nint,hIk)
do i_state=1,N_states
dIk(i_state) = hIk * lambda_mrcc(idx_alpha(k_sd),i_state)
enddo
! |l> = Exc(k -> alpha) |I>
call get_excitation(psi_non_cas(1,1,idx_alpha(k_sd)),tq(1,1,i_alpha),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)
tmp_det(1:Nint,1:2) = psi_cas(1,1,i_I)
! Hole (see list_to_bitstring)
iint = ishft(h1-1,-bit_kind_shift) + 1
ipos = h1-ishft((iint-1),bit_kind_shift)-1
@ -81,28 +87,37 @@ subroutine mrcc_dress(delta_ij_,Ndet,i_generator,n_selected,det_buffer,Nint,ipro
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)
! <I| \l/ |alpha>
do i_state=1,N_states
dka(i_state) = 0.d0
enddo
do l_sd=k_sd+1,idx_alpha(0)
call get_excitation_degree(tmp_det,psi_non_cas(1,1,idx_alpha(l_sd)),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
call get_excitation(psi_cas(1,1,i_I),psi_non_cas(1,1,idx_alpha(l_sd)),exc,degree,phase2,Nint)
call i_h_j(psi_cas(1,1,i_I),psi_non_cas(1,1,idx_alpha(l_sd)),Nint,hIl)
do i_state=1,N_states
dka(i_state) = hIl * lambda_mrcc(idx_alpha(l_sd),i_state) * phase * phase2
enddo
exit
endif
enddo
do l=1,N_states
dIa(l) += dka(l)*dIk(l)
do i_state=1,N_states
dIa(i_state) = dIa(i_state) + dIk(i_state) * dka(i_state)
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 i_state=1,N_states
ci_inv(i_state) = 1.d0/psi_cas_coefs(i_I,i_state)
enddo
do l_sd=1,idx_alpha(0)
k_sd = idx_alpha(l_sd)
call i_h_j(tq(1,1,i_alpha),psi_non_cas(1,1,idx_alpha(l_sd)),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)
delta_ij_(idx_non_cas(k_sd),idx_cas(i_I),m) += dIa(m) * hla
delta_ij_(idx_cas(i_I),idx_non_cas(k_sd),m) += dIa(m) * hla
delta_ij_(idx_cas(i_I),idx_cas(i_I),m) -= dIa(m) * hla * ci_inv(m) * psi_non_cas_coefs(k_sd,m)
enddo
enddo
enddo
@ -115,13 +130,13 @@ end
subroutine mrcc_dress_simple(delta_ij_sd_,Ndet_sd,i_generator,n_selected,det_buffer,Nint,iproc)
subroutine mrcc_dress_simple(delta_ij_non_cas_,Ndet_non_cas,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_sd
double precision, intent(inout) :: delta_ij_sd_(Ndet_sd,Ndet_sd,*)
integer, intent(in) :: Ndet_non_cas
double precision, intent(inout) :: delta_ij_non_cas_(Ndet_non_cas,Ndet_non_cas,*)
integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected)
integer :: i,j,k,m
@ -143,18 +158,18 @@ subroutine mrcc_dress_simple(delta_ij_sd_,Ndet_sd,i_generator,n_selected,det_buf
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 get_excitation_degree_vector(psi_non_cas,tq(1,1,i),degree,Nint,Ndet_non_cas,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)
call i_h_j(tq(1,1,i),psi_non_cas(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)
call i_h_j(tq(1,1,i),psi_non_cas(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)
delta_ij_non_cas_(idx(k), idx(j),m) += haj*hka* f(m)
delta_ij_non_cas_(idx(j), idx(k),m) += haj*hka* f(m)
enddo
enddo
enddo

145
src/MRCC/mrcc_utils.irp.f
View File

@ -1,130 +1,3 @@
use bitmasks
BEGIN_PROVIDER [ integer(bit_kind), cas_bitmask, (N_int,2,N_cas_bitmask) ]
implicit none
BEGIN_DOC
! Bitmasks for CAS reference determinants. (N_int, alpha/beta, CAS reference)
END_DOC
logical :: exists
integer :: i
PROVIDE ezfio_filename
call ezfio_has_bitmasks_cas(exists)
if (exists) then
call ezfio_get_bitmasks_cas(cas_bitmask)
else
do i=1,N_cas_bitmask
cas_bitmask(:,:,i) = iand(not(HF_bitmask(:,:)),full_ijkl_bitmask(:,:))
enddo
endif
END_PROVIDER
BEGIN_PROVIDER [ integer, N_det_cas ]
implicit none
BEGIN_DOC
! Number of generator detetrminants
END_DOC
integer :: i,k,l
logical :: good
call write_time(output_dets)
N_det_cas = 0
do i=1,N_det
do l=1,n_cas_bitmask
good = .True.
do k=1,N_int
good = good .and. ( &
iand(not(cas_bitmask(k,1,l)), psi_det(k,1,i)) == &
iand(not(cas_bitmask(k,1,l)), psi_det(k,1,1)) ) .and. ( &
iand(not(cas_bitmask(k,2,l)), psi_det(k,2,i)) == &
iand(not(cas_bitmask(k,2,l)), psi_det(k,2,1)) )
enddo
if (good) then
exit
endif
enddo
if (good) then
N_det_cas += 1
endif
enddo
N_det_cas = max(N_det_cas, 1)
call write_int(output_dets,N_det_cas, 'Number of determinants in the CAS')
END_PROVIDER
BEGIN_PROVIDER [ integer(bit_kind), psi_cas, (N_int,2,N_det_cas) ]
&BEGIN_PROVIDER [ double precision, psi_cas_coefs, (N_det_cas,n_states) ]
&BEGIN_PROVIDER [ integer, idx_cas, (N_det_cas) ]
implicit none
BEGIN_DOC
! For Single reference wave functions, the generator is the
! Hartree-Fock determinant
END_DOC
integer :: i, k, l, m
logical :: good
m=0
do i=1,N_det
do l=1,n_cas_bitmask
good = .True.
do k=1,N_int
good = good .and. ( &
iand(not(cas_bitmask(k,1,l)), psi_det(k,1,i)) == &
iand(not(cas_bitmask(k,1,l)), psi_det(k,1,1)) ) .and. ( &
iand(not(cas_bitmask(k,2,l)), psi_det(k,2,i)) == &
iand(not(cas_bitmask(k,2,l)), psi_det(k,2,1)) )
enddo
if (good) then
exit
endif
enddo
if (good) then
m = m+1
do k=1,N_int
psi_cas(k,1,m) = psi_det(k,1,i)
psi_cas(k,2,m) = psi_det(k,2,i)
enddo
idx_cas(m) = i
do k=1,N_states
psi_cas_coefs(m,k) = psi_coef(i,k)
enddo
endif
enddo
END_PROVIDER
BEGIN_PROVIDER [ integer(bit_kind), psi_sd, (N_int,2,N_det) ]
&BEGIN_PROVIDER [ double precision, psi_sd_coefs, (N_det,n_states) ]
&BEGIN_PROVIDER [ integer, idx_sd, (N_det) ]
&BEGIN_PROVIDER [ integer, N_det_sd]
implicit none
BEGIN_DOC
! SD
END_DOC
integer :: i_sd,j,k
integer :: degree
logical :: in_cas
i_sd =0
do k=1,N_det
in_cas = .False.
do j=1,N_det_cas
call get_excitation_degree(psi_cas(1,1,j), psi_det(1,1,k), degree, N_int)
if (degree == 0) then
in_cas = .True.
exit
endif
enddo
if (.not.in_cas) then
double precision :: hij
i_sd += 1
psi_sd(1:N_int,1:2,i_sd) = psi_det(1:N_int,1:2,k)
psi_sd_coefs(i_sd,1:N_states) = psi_coef(k,1:N_states)
idx_sd(i_sd) = k
endif
enddo
N_det_sd = i_sd
END_PROVIDER
BEGIN_PROVIDER [ double precision, lambda_mrcc, (psi_det_size,n_states) ]
implicit none
@ -133,15 +6,15 @@ BEGIN_PROVIDER [ double precision, lambda_mrcc, (psi_det_size,n_states) ]
END_DOC
integer :: i,k
double precision :: ihpsi(N_states)
do i=1,N_det_sd
call i_h_psi(psi_sd(1,1,i), psi_cas, psi_cas_coefs, N_int, N_det_cas, &
do i=1,N_det_non_cas
call i_h_psi(psi_non_cas(1,1,i), psi_cas, psi_cas_coefs, N_int, N_det_cas, &
size(psi_cas_coefs,1), n_states, ihpsi)
double precision :: hij
do k=1,N_states
if (dabs(ihpsi(k)) < 1.d-6) then
lambda_mrcc(i,k) = 0.d0
else
lambda_mrcc(i,k) = psi_sd_coefs(i,k)/ihpsi(k)
lambda_mrcc(i,k) = psi_non_cas_coefs(i,k)/ihpsi(k)
lambda_mrcc(i,k) = min( lambda_mrcc (i,k),0.d0 )
endif
enddo
@ -158,13 +31,13 @@ BEGIN_PROVIDER [ character*(32), dressing_type ]
dressing_type = "MRCC"
END_PROVIDER
BEGIN_PROVIDER [ double precision, delta_ij_sd, (N_det_sd, N_det_sd,N_states) ]
BEGIN_PROVIDER [ double precision, delta_ij_non_cas, (N_det_non_cas, N_det_non_cas,N_states) ]
implicit none
BEGIN_DOC
! Dressing matrix in SD basis
END_DOC
delta_ij_sd = 0.d0
call H_apply_mrcc_simple(delta_ij_sd,N_det_sd)
delta_ij_non_cas = 0.d0
call H_apply_mrcc_simple(delta_ij_non_cas,N_det_non_cas)
END_PROVIDER
BEGIN_PROVIDER [ double precision, delta_ij, (N_det,N_det,N_states) ]
@ -178,9 +51,9 @@ BEGIN_PROVIDER [ double precision, delta_ij, (N_det,N_det,N_states) ]
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)
do j=1,N_det_non_cas
do i=1,N_det_non_cas
delta_ij(idx_non_cas(i),idx_non_cas(j),m) = delta_ij_non_cas(i,j,m)
enddo
enddo
enddo