10
0
mirror of https://github.com/LCPQ/quantum_package synced 2024-06-25 22:52:15 +02:00

Corrected some bugs in MRCC

This commit is contained in:
Anthony Scemama 2016-11-14 23:57:23 +01:00
parent 9a9c5037bb
commit 80d0a9420e
4 changed files with 90 additions and 66 deletions

View File

@ -99,9 +99,9 @@ program fci_zmq
print *, 'N_states = ', N_states
do k=1,N_states
print *, 'State', k
print *, 'PT2 = ', pt2
print *, 'E = ', E_CI_before
print *, 'E+PT2 = ', E_CI_before+pt2
print *, 'PT2 = ', pt2(k)
print *, 'E = ', E_CI_before(k)
print *, 'E+PT2 = ', E_CI_before(k)+pt2(k)
print *, '-----'
enddo
call ezfio_set_cas_sd_zmq_energy_pt2(E_CI_before+pt2)

View File

@ -191,6 +191,15 @@ END_PROVIDER
end if
end do
if (a_col == at_row) then
t(:) = t(:) + 1.d0
endif
if (sum(dabs(t(:))) > 0.d0) then
wk = wk+1
A_ind_mwen(wk) = a_col
A_val_mwen(:,wk) = t(:)
endif
end do
if(wk /= 0) then

View File

@ -77,18 +77,18 @@ BEGIN_PROVIDER [ double precision, hij_mrcc, (N_det_non_ref,N_det_ref) ]
END_PROVIDER
BEGIN_PROVIDER [ double precision, delta_ij, (N_states,N_det_non_ref,N_det_ref) ]
&BEGIN_PROVIDER [ double precision, delta_ii, (N_states,N_det_ref) ]
implicit none
BEGIN_DOC
! Dressing matrix in N_det basis
END_DOC
integer :: i,j,m
delta_ij = 0.d0
delta_ii = 0.d0
call H_apply_mrcc(delta_ij,delta_ii,N_states,N_det_non_ref,N_det_ref)
END_PROVIDER
! BEGIN_PROVIDER [ double precision, delta_ij, (N_states,N_det_non_ref,N_det_ref) ]
!&BEGIN_PROVIDER [ double precision, delta_ii, (N_states,N_det_ref) ]
! implicit none
! BEGIN_DOC
! ! Dressing matrix in N_det basis
! END_DOC
! integer :: i,j,m
! delta_ij = 0.d0
! delta_ii = 0.d0
! call H_apply_mrcc(delta_ij,delta_ii,N_states,N_det_non_ref,N_det_ref)
!
!END_PROVIDER
BEGIN_PROVIDER [ double precision, h_matrix_dressed, (N_det,N_det,N_states) ]
@ -139,7 +139,6 @@ END_PROVIDER
integer :: mrcc_state
mrcc_state = N_states
do j=1,min(N_states,N_det)
do i=1,N_det
CI_eigenvectors_dressed(i,j) = psi_coef(i,j)
@ -148,16 +147,28 @@ END_PROVIDER
if (diag_algorithm == "Davidson") then
! call davidson_diag_mrcc(psi_det,CI_eigenvectors_dressed,CI_electronic_energy_dressed,&
! size(CI_eigenvectors_dressed,1),N_det,N_states,N_states_diag,N_int,output_determinants,mrcc_state)
call davidson_diag_mrcc_HS2(psi_det,CI_eigenvectors_dressed,&
size(CI_eigenvectors_dressed,1), &
CI_electronic_energy_dressed,N_det,N_states,N_states_diag,N_int, &
output_determinants,mrcc_state)
allocate (eigenvectors(size(CI_eigenvectors_dressed,1),size(CI_eigenvectors_dressed,2)), &
eigenvalues(size(CI_electronic_energy_dressed,1)))
do mrcc_state=N_states,1,-1
do j=1,min(N_states,N_det)
do i=1,N_det
eigenvectors(i,j) = psi_coef(i,j)
enddo
enddo
call davidson_diag_mrcc_HS2(psi_det,eigenvectors,&
size(eigenvectors,1), &
eigenvalues,N_det,N_states,N_states_diag,N_int, &
output_determinants,mrcc_state)
CI_eigenvectors_dressed(1:N_det,mrcc_state) = eigenvectors(1:N_det,mrcc_state)
CI_electronic_energy_dressed(mrcc_state) = eigenvalues(mrcc_state)
enddo
do mrcc_state=N_states+1,N_states_diag
CI_eigenvectors_dressed(1:N_det,mrcc_state) = eigenvectors(1:N_det,mrcc_state)
CI_electronic_energy_dressed(mrcc_state) = eigenvalues(mrcc_state)
enddo
call u_0_S2_u_0(CI_eigenvectors_s2_dressed,CI_eigenvectors_dressed,N_det,psi_det,N_int,&
N_states_diag,size(CI_eigenvectors_dressed,1))
N_states_diag,size(CI_eigenvectors_dressed,1))
deallocate (eigenvectors,eigenvalues)
else if (diag_algorithm == "Lapack") then
@ -628,12 +639,12 @@ END_PROVIDER
double precision :: phase
double precision, allocatable :: rho_mrcc_init(:,:)
double precision, allocatable :: rho_mrcc_init(:)
integer :: a_coll, at_roww
print *, "TI", hh_nex, N_det_non_ref
allocate(rho_mrcc_init(N_det_non_ref, N_states))
allocate(rho_mrcc_init(N_det_non_ref))
allocate(x_new(hh_nex))
allocate(x(hh_nex), AtB(hh_nex))
x = 0d0
@ -644,9 +655,8 @@ END_PROVIDER
AtB(:) = 0.d0
!$OMP PARALLEL default(none) shared(k, psi_non_ref_coef, active_excitation_to_determinants_idx,&
!$OMP active_excitation_to_determinants_val, x, N_det_ref, hh_nex, N_det_non_ref) &
!$OMP private(at_row, a_col, t, i, j, r1, r2, wk, A_ind_mwen, A_val_mwen, a_coll, at_roww)&
!$OMP private(at_row, a_col, i, j, r1, r2, wk, A_ind_mwen, A_val_mwen, a_coll, at_roww)&
!$OMP shared(N_states,mrcc_col_shortcut, mrcc_N_col, AtB, mrcc_AtA_val, mrcc_AtA_ind, s, n_exc_active, active_pp_idx)
allocate(A_val_mwen(N_states,hh_nex), A_ind_mwen(hh_nex), t(N_states))
!$OMP DO schedule(dynamic, 100)
do at_roww = 1, n_exc_active ! hh_nex
@ -655,11 +665,11 @@ END_PROVIDER
AtB(at_row) = AtB(at_row) + psi_non_ref_coef(active_excitation_to_determinants_idx(i, at_roww), s) * active_excitation_to_determinants_val(s,i, at_roww)
end do
end do
!$OMP END DO NOWAIT
deallocate (A_ind_mwen, A_val_mwen)
!$OMP END DO
!$OMP END PARALLEL
x = 0d0
X(:) = 0d0
do a_coll = 1, n_exc_active
@ -669,10 +679,7 @@ END_PROVIDER
rho_mrcc_init = 0d0
!$OMP PARALLEL default(shared) &
!$OMP private(lref, hh, pp, II, myMask, myDet, ok, ind, phase)
allocate(lref(N_det_ref))
!$OMP DO schedule(static, 1)
do hh = 1, hh_shortcut(0)
do pp = hh_shortcut(hh), hh_shortcut(hh+1)-1
if(is_active_exc(pp)) cycle
@ -694,16 +701,14 @@ END_PROVIDER
X(pp) = AtB(pp) / X(pp)
do II=1,N_det_ref
if(lref(II) > 0) then
rho_mrcc_init(lref(II),s) = psi_ref_coef(II,s) * X(pp)
rho_mrcc_init(lref(II)) = psi_ref_coef(II,s) * X(pp)
else if(lref(II) < 0) then
rho_mrcc_init(-lref(II),s) = -psi_ref_coef(II,s) * X(pp)
rho_mrcc_init(-lref(II)) = -psi_ref_coef(II,s) * X(pp)
end if
end do
end do
end do
!$OMP END DO
deallocate(lref)
!$OMP END PARALLEL
x_new = x
@ -716,9 +721,9 @@ END_PROVIDER
!$OMP DO
do i=1,N_det_non_ref
rho_mrcc(i,s) = rho_mrcc_init(i,s)
rho_mrcc(i,s) = rho_mrcc_init(i)
enddo
!$OMP END DO
!$OMP END DO NOWAIT
!$OMP DO
do a_coll = 1, n_exc_active
@ -928,7 +933,7 @@ END_PROVIDER
! rho_mrcc now contains the product of the scaling factors and the
! normalization constant
dIj_unique(:size(X), s) = X(:)
dIj_unique(1:size(X), s) = X(1:size(X))
end do
END_PROVIDER

View File

@ -317,43 +317,53 @@ end
&BEGIN_PROVIDER [ double precision, delta_ii, (N_states, N_det_ref) ]
use bitmasks
implicit none
integer :: i, j, i_state
integer :: i, j, i_state
!mrmode : 1=mrcepa0, 2=mrsc2 add, 3=mrcc
do i_state = 1, N_states
if(mrmode == 3) then
if(mrmode == 3) then
do i = 1, N_det_ref
delta_ii(i_state,i)= delta_ii_mrcc(i_state,i)
do i_state = 1, N_states
delta_ii(i_state,i)= delta_ii_mrcc(i_state,i)
enddo
do j = 1, N_det_non_ref
delta_ij(i_state,j,i) = delta_ij_mrcc(i_state,j,i)
do i_state = 1, N_states
delta_ij(i_state,j,i) = delta_ij_mrcc(i_state,j,i)
enddo
end do
end do
!
! do i = 1, N_det_ref
! delta_ii(i_state,i)= delta_mrcepa0_ii(i,i_state) - delta_sub_ii(i,i_state)
! do j = 1, N_det_non_ref
! delta_ij(i_state,j,i) = delta_mrcepa0_ij(i,j,i_state) - delta_sub_ij(i,j,i_state)
! end do
! end do
else if(mrmode == 2) then
do i = 1, N_det_ref
!
! do i = 1, N_det_ref
! delta_ii(i_state,i)= delta_mrcepa0_ii(i,i_state) - delta_sub_ii(i,i_state)
! do j = 1, N_det_non_ref
! delta_ij(i_state,j,i) = delta_mrcepa0_ij(i,j,i_state) - delta_sub_ij(i,j,i_state)
! end do
! end do
else if(mrmode == 2) then
do i = 1, N_det_ref
do i_state = 1, N_states
delta_ii(i_state,i)= delta_ii_old(i_state,i)
do j = 1, N_det_non_ref
enddo
do j = 1, N_det_non_ref
do i_state = 1, N_states
delta_ij(i_state,j,i) = delta_ij_old(i_state,j,i)
end do
enddo
end do
else if(mrmode == 1) then
do i = 1, N_det_ref
end do
else if(mrmode == 1) then
do i = 1, N_det_ref
do i_state = 1, N_states
delta_ii(i_state,i)= delta_mrcepa0_ii(i,i_state)
do j = 1, N_det_non_ref
enddo
do j = 1, N_det_non_ref
do i_state = 1, N_states
delta_ij(i_state,j,i) = delta_mrcepa0_ij(i,j,i_state)
end do
enddo
end do
else
stop "invalid mrmode"
end if
end do
end do
else
stop "invalid mrmode"
end if
END_PROVIDER