mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-12-22 12:23:43 +01:00
compiling after some cleaning
This commit is contained in:
eginer
parent
b749796d93
commit
366afb2933
1
plugins/local/old_delta_tc_qmc/NEED
Normal file
1
plugins/local/old_delta_tc_qmc/NEED
Normal file
@ -0,0 +1 @@
|
||||
|
||||
4
plugins/local/old_delta_tc_qmc/README.rst
Normal file
4
plugins/local/old_delta_tc_qmc/README.rst
Normal file
@ -0,0 +1,4 @@
|
||||
================
|
||||
old_delta_tc_qmc
|
||||
================
|
||||
|
||||
7
plugins/local/old_delta_tc_qmc/old_delta_tc_qmc.irp.f
Normal file
7
plugins/local/old_delta_tc_qmc/old_delta_tc_qmc.irp.f
Normal file
@ -0,0 +1,7 @@
|
||||
program old_delta_tc_qmc
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! TODO : Put the documentation of the program here
|
||||
END_DOC
|
||||
print *, 'Hello world'
|
||||
end
|
||||
198
plugins/local/slater_tc/h_mat_triple.irp.f
Normal file
198
plugins/local/slater_tc/h_mat_triple.irp.f
Normal file
@ -0,0 +1,198 @@
|
||||
subroutine H_tc_s2_u_0_with_pure_three(v_0, s_0, u_0, N_st, sze)
|
||||
BEGIN_DOC
|
||||
! Computes $v_0 = H^TC | u_0\rangle$ WITH PURE TRIPLE EXCITATION TERMS
|
||||
!
|
||||
! Assumes that the determinants are in psi_det
|
||||
!
|
||||
! istart, iend, ishift, istep are used in ZMQ parallelization.
|
||||
END_DOC
|
||||
|
||||
use bitmasks
|
||||
implicit none
|
||||
|
||||
integer, intent(in) :: N_st,sze
|
||||
double precision, intent(in) :: u_0(sze,N_st)
|
||||
double precision, intent(out) :: v_0(sze,N_st), s_0(sze,N_st)
|
||||
call H_tc_s2_u_0_opt(v_0, s_0, u_0, N_st, sze)
|
||||
integer :: i,j,degree,ist
|
||||
double precision :: hmono, htwoe, hthree, htot
|
||||
do i = 1, N_det
|
||||
do j = 1, N_det
|
||||
call get_excitation_degree(psi_det(1,1,i),psi_det(1,1,j),degree,N_int)
|
||||
if(degree .ne. 3)cycle
|
||||
call triple_htilde_mu_mat_fock_bi_ortho(N_int, psi_det(1,1,i), psi_det(1,1,j), hmono, htwoe, hthree, htot)
|
||||
do ist = 1, N_st
|
||||
v_0(i,ist) += htot * u_0(j,ist)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
end
|
||||
|
||||
subroutine H_tc_s2_u_0_with_pure_three_omp(v_0, s_0, u_0, N_st, sze)
|
||||
BEGIN_DOC
|
||||
! Computes $v_0 = H^TC | u_0\rangle$ WITH PURE TRIPLE EXCITATION TERMS
|
||||
!
|
||||
! Assumes that the determinants are in psi_det
|
||||
!
|
||||
! istart, iend, ishift, istep are used in ZMQ parallelization.
|
||||
END_DOC
|
||||
|
||||
use bitmasks
|
||||
implicit none
|
||||
|
||||
integer, intent(in) :: N_st,sze
|
||||
double precision, intent(in) :: u_0(sze,N_st)
|
||||
double precision, intent(out) :: v_0(sze,N_st), s_0(sze,N_st)
|
||||
call H_tc_s2_u_0_opt(v_0, s_0, u_0, N_st, sze)
|
||||
integer :: i,j,degree,ist
|
||||
double precision :: hmono, htwoe, hthree, htot
|
||||
!$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(dynamic,8) &
|
||||
!$OMP SHARED(N_st, N_det, N_int, psi_det, u_0, v_0) &
|
||||
!$OMP PRIVATE(ist, i, j, degree, hmono, htwoe, hthree,htot)
|
||||
do i = 1, N_det
|
||||
do j = 1, N_det
|
||||
call get_excitation_degree(psi_det(1,1,i),psi_det(1,1,j),degree,N_int)
|
||||
if(degree .ne. 3)cycle
|
||||
call triple_htilde_mu_mat_fock_bi_ortho(N_int, psi_det(1,1,i), psi_det(1,1,j), hmono, htwoe, hthree, htot)
|
||||
do ist = 1, N_st
|
||||
v_0(i,ist) += htot * u_0(j,ist)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END PARALLEL DO
|
||||
end
|
||||
|
||||
! ---
|
||||
|
||||
subroutine H_tc_s2_dagger_u_0_with_pure_three(v_0, s_0, u_0, N_st, sze)
|
||||
BEGIN_DOC
|
||||
! Computes $v_0 = (H^TC)^dagger | u_0\rangle$ WITH PURE TRIPLE EXCITATION TERMS
|
||||
!
|
||||
! Assumes that the determinants are in psi_det
|
||||
!
|
||||
! istart, iend, ishift, istep are used in ZMQ parallelization.
|
||||
END_DOC
|
||||
|
||||
use bitmasks
|
||||
implicit none
|
||||
|
||||
integer, intent(in) :: N_st,sze
|
||||
double precision, intent(in) :: u_0(sze,N_st)
|
||||
double precision, intent(out) :: v_0(sze,N_st), s_0(sze,N_st)
|
||||
call H_tc_s2_dagger_u_0_opt(v_0, s_0, u_0, N_st, sze)
|
||||
integer :: i,j,degree,ist
|
||||
double precision :: hmono, htwoe, hthree, htot
|
||||
do i = 1, N_det
|
||||
do j = 1, N_det
|
||||
call get_excitation_degree(psi_det(1,1,i),psi_det(1,1,j),degree,N_int)
|
||||
if(degree .ne. 3)cycle
|
||||
call triple_htilde_mu_mat_fock_bi_ortho(N_int, psi_det(1,1,j), psi_det(1,1,i), hmono, htwoe, hthree, htot)
|
||||
do ist = 1, N_st
|
||||
v_0(i,ist) += htot * u_0(j,ist)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
end
|
||||
|
||||
subroutine H_tc_s2_dagger_u_0_with_pure_three_omp(v_0, s_0, u_0, N_st, sze)
|
||||
BEGIN_DOC
|
||||
! Computes $v_0 = (H^TC)^dagger | u_0\rangle$ WITH PURE TRIPLE EXCITATION TERMS
|
||||
!
|
||||
! Assumes that the determinants are in psi_det
|
||||
!
|
||||
! istart, iend, ishift, istep are used in ZMQ parallelization.
|
||||
END_DOC
|
||||
|
||||
use bitmasks
|
||||
implicit none
|
||||
|
||||
integer, intent(in) :: N_st,sze
|
||||
double precision, intent(in) :: u_0(sze,N_st)
|
||||
double precision, intent(out) :: v_0(sze,N_st), s_0(sze,N_st)
|
||||
call H_tc_s2_dagger_u_0_opt(v_0, s_0, u_0, N_st, sze)
|
||||
integer :: i,j,degree,ist
|
||||
double precision :: hmono, htwoe, hthree, htot
|
||||
!$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(dynamic,8) &
|
||||
!$OMP SHARED(N_st, N_det, N_int, psi_det, u_0, v_0) &
|
||||
!$OMP PRIVATE(ist, i, j, degree, hmono, htwoe, hthree,htot)
|
||||
do i = 1, N_det
|
||||
do j = 1, N_det
|
||||
call get_excitation_degree(psi_det(1,1,i),psi_det(1,1,j),degree,N_int)
|
||||
if(degree .ne. 3)cycle
|
||||
call triple_htilde_mu_mat_fock_bi_ortho(N_int, psi_det(1,1,j), psi_det(1,1,i), hmono, htwoe, hthree, htot)
|
||||
do ist = 1, N_st
|
||||
v_0(i,ist) += htot * u_0(j,ist)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END PARALLEL DO
|
||||
end
|
||||
|
||||
! ---
|
||||
subroutine triple_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, hthree, htot)
|
||||
use bitmasks
|
||||
BEGIN_DOC
|
||||
! <key_j | H_tilde | key_i> for triple excitation
|
||||
!!
|
||||
!! WARNING !!
|
||||
!
|
||||
! Genuine triple excitations of the same spin are not yet implemented
|
||||
END_DOC
|
||||
implicit none
|
||||
integer(bit_kind), intent(in) :: key_j(N_int,2),key_i(N_int,2)
|
||||
integer, intent(in) :: Nint
|
||||
double precision, intent(out) :: hmono, htwoe, hthree, htot
|
||||
integer :: degree
|
||||
integer :: h1, p1, h2, p2, s1, s2, h3, p3, s3
|
||||
integer :: holes_array(100,2),particles_array(100,2),degree_array(2)
|
||||
double precision :: phase,sym_3_e_int_from_6_idx_tensor
|
||||
|
||||
hmono = 0.d0
|
||||
htwoe = 0.d0
|
||||
hthree = 0.d0
|
||||
htot = 0.d0
|
||||
call get_excitation_general(key_j, key_i, Nint,degree_array,holes_array, particles_array,phase)
|
||||
degree = degree_array(1) + degree_array(2)
|
||||
if(degree .ne. 3)return
|
||||
if(degree_array(1)==3.or.degree_array(2)==3)then
|
||||
if(degree_array(1) == 3)then
|
||||
h1 = holes_array(1,1)
|
||||
h2 = holes_array(2,1)
|
||||
h3 = holes_array(3,1)
|
||||
p1 = particles_array(1,1)
|
||||
p2 = particles_array(2,1)
|
||||
p3 = particles_array(3,1)
|
||||
else
|
||||
h1 = holes_array(1,2)
|
||||
h2 = holes_array(2,2)
|
||||
h3 = holes_array(3,2)
|
||||
p1 = particles_array(1,2)
|
||||
p2 = particles_array(2,2)
|
||||
p3 = particles_array(3,2)
|
||||
endif
|
||||
hthree = sym_3_e_int_from_6_idx_tensor(p3, p2, p1, h3, h2, h1)
|
||||
else
|
||||
if(degree_array(1) == 2.and.degree_array(2) == 1)then ! double alpha + single beta
|
||||
h1 = holes_array(1,1)
|
||||
h2 = holes_array(2,1)
|
||||
h3 = holes_array(1,2)
|
||||
p1 = particles_array(1,1)
|
||||
p2 = particles_array(2,1)
|
||||
p3 = particles_array(1,2)
|
||||
else if(degree_array(2) == 2 .and. degree_array(1) == 1)then ! double beta + single alpha
|
||||
h1 = holes_array(1,2)
|
||||
h2 = holes_array(2,2)
|
||||
h3 = holes_array(1,1)
|
||||
p1 = particles_array(1,2)
|
||||
p2 = particles_array(2,2)
|
||||
p3 = particles_array(1,1)
|
||||
else
|
||||
print*,'PB !!'
|
||||
stop
|
||||
endif
|
||||
hthree = three_body_ints_bi_ort(p3,p2,p1,h3,h2,h1) - three_body_ints_bi_ort(p3,p2,p1,h3,h1,h2)
|
||||
endif
|
||||
hthree *= phase
|
||||
htot = hthree
|
||||
end
|
||||
|
||||
@ -1,193 +0,0 @@
|
||||
subroutine get_excitation_general(key_i,key_j, Nint,degree_array,holes_array, particles_array,phase)
|
||||
use bitmasks
|
||||
BEGIN_DOC
|
||||
! returns the array, for each spin, of holes/particles between key_i and key_j
|
||||
!
|
||||
! with the following convention: a^+_{particle} a_{hole}|key_i> = |key_j>
|
||||
END_DOC
|
||||
include 'utils/constants.include.F'
|
||||
implicit none
|
||||
integer, intent(in) :: Nint
|
||||
integer(bit_kind), intent(in) :: key_j(Nint,2),key_i(Nint,2)
|
||||
integer, intent(out) :: holes_array(100,2),particles_array(100,2),degree_array(2)
|
||||
double precision, intent(out) :: phase
|
||||
integer :: ispin,k,i,pos
|
||||
integer(bit_kind) :: key_hole, key_particle
|
||||
integer(bit_kind) :: xorvec(N_int_max,2)
|
||||
holes_array = -1
|
||||
particles_array = -1
|
||||
degree_array = 0
|
||||
do i = 1, N_int
|
||||
xorvec(i,1) = xor( key_i(i,1), key_j(i,1))
|
||||
xorvec(i,2) = xor( key_i(i,2), key_j(i,2))
|
||||
degree_array(1) += popcnt(xorvec(i,1))
|
||||
degree_array(2) += popcnt(xorvec(i,2))
|
||||
enddo
|
||||
degree_array(1) = shiftr(degree_array(1),1)
|
||||
degree_array(2) = shiftr(degree_array(2),1)
|
||||
|
||||
do ispin = 1, 2
|
||||
k = 1
|
||||
!!! GETTING THE HOLES
|
||||
do i = 1, N_int
|
||||
key_hole = iand(xorvec(i,ispin),key_i(i,ispin))
|
||||
do while(key_hole .ne.0_bit_kind)
|
||||
pos = trailz(key_hole)
|
||||
holes_array(k,ispin) = 1+ bit_kind_size * (i-1) + pos
|
||||
key_hole = ibclr(key_hole,pos)
|
||||
k += 1
|
||||
if(k .gt.100)then
|
||||
print*,'WARNING in get_excitation_general'
|
||||
print*,'More than a 100-th excitation for spin ',ispin
|
||||
print*,'stoping ...'
|
||||
stop
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
do ispin = 1, 2
|
||||
k = 1
|
||||
!!! GETTING THE PARTICLES
|
||||
do i = 1, N_int
|
||||
key_particle = iand(xor(key_i(i,ispin),key_j(i,ispin)),key_j(i,ispin))
|
||||
do while(key_particle .ne.0_bit_kind)
|
||||
pos = trailz(key_particle)
|
||||
particles_array(k,ispin) = 1+ bit_kind_size * (i-1) + pos
|
||||
key_particle = ibclr(key_particle,pos)
|
||||
k += 1
|
||||
if(k .gt.100)then
|
||||
print*,'WARNING in get_excitation_general '
|
||||
print*,'More than a 100-th excitation for spin ',ispin
|
||||
print*,'stoping ...'
|
||||
stop
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
integer :: h,p, i_ok
|
||||
integer(bit_kind), allocatable :: det_i(:,:),det_ip(:,:)
|
||||
integer :: exc(0:2,2,2)
|
||||
double precision :: phase_tmp
|
||||
allocate(det_i(Nint,2),det_ip(N_int,2))
|
||||
det_i = key_i
|
||||
phase = 1.d0
|
||||
do ispin = 1, 2
|
||||
do i = 1, degree_array(ispin)
|
||||
h = holes_array(i,ispin)
|
||||
p = particles_array(i,ispin)
|
||||
det_ip = det_i
|
||||
call do_single_excitation(det_ip,h,p,ispin,i_ok)
|
||||
if(i_ok == -1)then
|
||||
print*,'excitation was not possible '
|
||||
stop
|
||||
endif
|
||||
call get_single_excitation(det_i,det_ip,exc,phase_tmp,Nint)
|
||||
phase *= phase_tmp
|
||||
det_i = det_ip
|
||||
enddo
|
||||
enddo
|
||||
|
||||
end
|
||||
|
||||
subroutine get_holes_general(key_i, key_j,Nint, holes_array)
|
||||
use bitmasks
|
||||
BEGIN_DOC
|
||||
! returns the array, per spin, of holes between key_i and key_j
|
||||
!
|
||||
! with the following convention: a_{hole}|key_i> --> |key_j>
|
||||
END_DOC
|
||||
implicit none
|
||||
integer, intent(in) :: Nint
|
||||
integer(bit_kind), intent(in) :: key_j(Nint,2),key_i(Nint,2)
|
||||
integer, intent(out) :: holes_array(100,2)
|
||||
integer(bit_kind) :: key_hole
|
||||
integer :: ispin,k,i,pos
|
||||
holes_array = -1
|
||||
do ispin = 1, 2
|
||||
k = 1
|
||||
do i = 1, N_int
|
||||
key_hole = iand(xor(key_i(i,ispin),key_j(i,ispin)),key_i(i,ispin))
|
||||
do while(key_hole .ne.0_bit_kind)
|
||||
pos = trailz(key_hole)
|
||||
holes_array(k,ispin) = 1+ bit_kind_size * (i-1) + pos
|
||||
key_hole = ibclr(key_hole,pos)
|
||||
k += 1
|
||||
if(k .gt.100)then
|
||||
print*,'WARNING in get_holes_general'
|
||||
print*,'More than a 100-th excitation for spin ',ispin
|
||||
print*,'stoping ...'
|
||||
stop
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
end
|
||||
|
||||
subroutine get_particles_general(key_i, key_j,Nint,particles_array)
|
||||
use bitmasks
|
||||
BEGIN_DOC
|
||||
! returns the array, per spin, of particles between key_i and key_j
|
||||
!
|
||||
! with the following convention: a^dagger_{particle}|key_i> --> |key_j>
|
||||
END_DOC
|
||||
implicit none
|
||||
integer, intent(in) :: Nint
|
||||
integer(bit_kind), intent(in) :: key_j(Nint,2),key_i(Nint,2)
|
||||
integer, intent(out) :: particles_array(100,2)
|
||||
integer(bit_kind) :: key_particle
|
||||
integer :: ispin,k,i,pos
|
||||
particles_array = -1
|
||||
do ispin = 1, 2
|
||||
k = 1
|
||||
do i = 1, N_int
|
||||
key_particle = iand(xor(key_i(i,ispin),key_j(i,ispin)),key_j(i,ispin))
|
||||
do while(key_particle .ne.0_bit_kind)
|
||||
pos = trailz(key_particle)
|
||||
particles_array(k,ispin) = 1+ bit_kind_size * (i-1) + pos
|
||||
key_particle = ibclr(key_particle,pos)
|
||||
k += 1
|
||||
if(k .gt.100)then
|
||||
print*,'WARNING in get_holes_general'
|
||||
print*,'More than a 100-th excitation for spin ',ispin
|
||||
print*,'Those are the two determinants'
|
||||
call debug_det(key_i, N_int)
|
||||
call debug_det(key_j, N_int)
|
||||
print*,'stoping ...'
|
||||
stop
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
end
|
||||
|
||||
subroutine get_phase_general(key_i,Nint,degree, holes_array, particles_array,phase)
|
||||
implicit none
|
||||
integer, intent(in) :: degree(2), Nint
|
||||
integer(bit_kind), intent(in) :: key_i(Nint,2)
|
||||
integer, intent(in) :: holes_array(100,2),particles_array(100,2)
|
||||
double precision, intent(out) :: phase
|
||||
integer :: i,ispin,h,p, i_ok
|
||||
integer(bit_kind), allocatable :: det_i(:,:),det_ip(:,:)
|
||||
integer :: exc(0:2,2,2)
|
||||
double precision :: phase_tmp
|
||||
allocate(det_i(Nint,2),det_ip(N_int,2))
|
||||
det_i = key_i
|
||||
phase = 1.d0
|
||||
do ispin = 1, 2
|
||||
do i = 1, degree(ispin)
|
||||
h = holes_array(i,ispin)
|
||||
p = particles_array(i,ispin)
|
||||
det_ip = det_i
|
||||
call do_single_excitation(det_ip,h,p,ispin,i_ok)
|
||||
if(i_ok == -1)then
|
||||
print*,'excitation was not possible '
|
||||
stop
|
||||
endif
|
||||
call get_single_excitation(det_i,det_ip,exc,phase_tmp,Nint)
|
||||
phase *= phase_tmp
|
||||
det_i = det_ip
|
||||
enddo
|
||||
enddo
|
||||
|
||||
end
|
||||
|
||||
@ -1,129 +0,0 @@
|
||||
program pt2_tc_cisd
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! TODO : Reads psi_det in the EZFIO folder and prints out the left- and right-eigenvectors together
|
||||
! with the energy. Saves the left-right wave functions at the end.
|
||||
!
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
|
||||
my_grid_becke = .True.
|
||||
PROVIDE tc_grid1_a tc_grid1_r
|
||||
my_n_pt_r_grid = tc_grid1_r
|
||||
my_n_pt_a_grid = tc_grid1_a
|
||||
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
|
||||
|
||||
read_wf = .True.
|
||||
touch read_wf
|
||||
|
||||
print*, ' nb of states = ', N_states
|
||||
print*, ' nb of det = ', N_det
|
||||
call routine_diag()
|
||||
|
||||
call routine
|
||||
end
|
||||
|
||||
subroutine routine
|
||||
implicit none
|
||||
integer :: i,h1,p1,h2,p2,s1,s2,degree
|
||||
double precision :: h0i,hi0,e00,ei,delta_e
|
||||
double precision :: norm,e_corr,coef,e_corr_pos,e_corr_neg,e_corr_abs
|
||||
|
||||
integer :: exc(0:2,2,2)
|
||||
double precision :: phase
|
||||
double precision :: eh1,ep1,eh2,ep2
|
||||
|
||||
norm = 0.d0
|
||||
e_corr = 0.d0
|
||||
e_corr_abs = 0.d0
|
||||
e_corr_pos = 0.d0
|
||||
e_corr_neg = 0.d0
|
||||
call htilde_mu_mat_bi_ortho_tot_slow(psi_det(1,1,1), psi_det(1,1,1), N_int, e00)
|
||||
do i = 2, N_det
|
||||
call htilde_mu_mat_bi_ortho_tot_slow(psi_det(1,1,i), psi_det(1,1,1), N_int, hi0)
|
||||
call htilde_mu_mat_bi_ortho_tot_slow(psi_det(1,1,1), psi_det(1,1,i), N_int, h0i)
|
||||
call htilde_mu_mat_bi_ortho_tot_slow(psi_det(1,1,i), psi_det(1,1,i), N_int, ei)
|
||||
call get_excitation_degree(psi_det(1,1,1), psi_det(1,1,i),degree,N_int)
|
||||
call get_excitation(psi_det(1,1,1), psi_det(1,1,i),exc,degree,phase,N_int)
|
||||
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
|
||||
eh1 = Fock_matrix_tc_diag_mo_tot(h1)
|
||||
ep1 = Fock_matrix_tc_diag_mo_tot(p1)
|
||||
delta_e = eh1 - ep1
|
||||
if (degree==2)then
|
||||
eh2 = Fock_matrix_tc_diag_mo_tot(h2)
|
||||
ep2 = Fock_matrix_tc_diag_mo_tot(p2)
|
||||
delta_e += eh2 - ep2
|
||||
endif
|
||||
! delta_e = e00 - ei
|
||||
coef = hi0/delta_e
|
||||
norm += coef*coef
|
||||
e_corr = coef* h0i
|
||||
if(e_corr.lt.0.d0)then
|
||||
e_corr_neg += e_corr
|
||||
elseif(e_corr.gt.0.d0)then
|
||||
e_corr_pos += e_corr
|
||||
endif
|
||||
e_corr_abs += dabs(e_corr)
|
||||
enddo
|
||||
print*,'e_corr_abs = ',e_corr_abs
|
||||
print*,'e_corr_pos = ',e_corr_pos
|
||||
print*,'e_corr_neg = ',e_corr_neg
|
||||
print*,'norm = ',dsqrt(norm)
|
||||
|
||||
end
|
||||
|
||||
subroutine routine_diag()
|
||||
|
||||
implicit none
|
||||
integer :: i, j, k
|
||||
double precision :: dE
|
||||
|
||||
! provide eigval_right_tc_bi_orth
|
||||
! provide overlap_bi_ortho
|
||||
! provide htilde_matrix_elmt_bi_ortho
|
||||
|
||||
if(N_states .eq. 1) then
|
||||
|
||||
print*,'eigval_right_tc_bi_orth = ',eigval_right_tc_bi_orth(1)
|
||||
print*,'e_tc_left_right = ',e_tc_left_right
|
||||
print*,'e_tilde_bi_orth_00 = ',e_tilde_bi_orth_00
|
||||
print*,'e_pt2_tc_bi_orth = ',e_pt2_tc_bi_orth
|
||||
print*,'e_pt2_tc_bi_orth_single = ',e_pt2_tc_bi_orth_single
|
||||
print*,'e_pt2_tc_bi_orth_double = ',e_pt2_tc_bi_orth_double
|
||||
print*,'***'
|
||||
print*,'e_corr_bi_orth = ',e_corr_bi_orth
|
||||
print*,'e_corr_bi_orth_proj = ',e_corr_bi_orth_proj
|
||||
print*,'e_corr_bi_orth_proj_abs = ',e_corr_bi_orth_proj_abs
|
||||
print*,'e_corr_single_bi_orth = ',e_corr_single_bi_orth
|
||||
print*,'e_corr_double_bi_orth = ',e_corr_double_bi_orth
|
||||
print*,'e_corr_single_bi_orth_abs = ',e_corr_single_bi_orth_abs
|
||||
print*,'e_corr_double_bi_orth_abs = ',e_corr_double_bi_orth_abs
|
||||
print*,'Left/right eigenvectors'
|
||||
do i = 1,N_det
|
||||
write(*,'(I5,X,(100(F12.7,X)))')i,leigvec_tc_bi_orth(i,1),reigvec_tc_bi_orth(i,1),leigvec_tc_bi_orth(i,1)*reigvec_tc_bi_orth(i,1)
|
||||
enddo
|
||||
|
||||
else
|
||||
|
||||
print*,'eigval_right_tc_bi_orth : '
|
||||
do i = 1, N_states
|
||||
print*, i, eigval_right_tc_bi_orth(i)
|
||||
enddo
|
||||
|
||||
print*,''
|
||||
print*,'******************************************************'
|
||||
print*,'TC Excitation energies (au) (eV)'
|
||||
do i = 2, N_states
|
||||
dE = eigval_right_tc_bi_orth(i) - eigval_right_tc_bi_orth(1)
|
||||
print*, i, dE, dE/0.0367502d0
|
||||
enddo
|
||||
print*,''
|
||||
|
||||
endif
|
||||
|
||||
end
|
||||
|
||||
|
||||
|
||||
@ -1,36 +0,0 @@
|
||||
|
||||
! ---
|
||||
|
||||
program tc_cisd_sc2
|
||||
|
||||
BEGIN_DOC
|
||||
! TODO : Put the documentation of the program here
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
|
||||
print *, 'Hello world'
|
||||
|
||||
my_grid_becke = .True.
|
||||
PROVIDE tc_grid1_a tc_grid1_r
|
||||
my_n_pt_r_grid = tc_grid1_r
|
||||
my_n_pt_a_grid = tc_grid1_a
|
||||
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
|
||||
|
||||
read_wf = .True.
|
||||
touch read_wf
|
||||
|
||||
call test
|
||||
|
||||
end
|
||||
|
||||
! ---
|
||||
|
||||
subroutine test()
|
||||
implicit none
|
||||
! double precision, allocatable :: dressing_dets(:),e_corr_dets(:)
|
||||
! allocate(dressing_dets(N_det),e_corr_dets(N_det))
|
||||
! e_corr_dets = 0.d0
|
||||
! call get_cisd_sc2_dressing(psi_det,e_corr_dets,N_det,dressing_dets)
|
||||
provide eigval_tc_cisd_sc2_bi_ortho
|
||||
end
|
||||
@ -1,145 +0,0 @@
|
||||
BEGIN_PROVIDER [ double precision, reigvec_tc_cisd_sc2_bi_ortho, (N_det,N_states)]
|
||||
&BEGIN_PROVIDER [ double precision, leigvec_tc_cisd_sc2_bi_ortho, (N_det,N_states)]
|
||||
&BEGIN_PROVIDER [ double precision, eigval_tc_cisd_sc2_bi_ortho, (N_states)]
|
||||
implicit none
|
||||
integer :: it,n_real,degree,i,istate
|
||||
double precision :: e_before, e_current,thr, hmono,htwoe,hthree,accu
|
||||
double precision, allocatable :: e_corr_dets(:),h0j(:), h_sc2(:,:), dressing_dets(:)
|
||||
double precision, allocatable :: leigvec_tc_bi_orth_tmp(:,:),reigvec_tc_bi_orth_tmp(:,:),eigval_right_tmp(:)
|
||||
allocate(leigvec_tc_bi_orth_tmp(N_det,N_det),reigvec_tc_bi_orth_tmp(N_det,N_det),eigval_right_tmp(N_det))
|
||||
allocate(e_corr_dets(N_det),h0j(N_det),h_sc2(N_det,N_det),dressing_dets(N_det))
|
||||
allocate(H_jj(N_det),vec_tmp(N_det,n_states_diag),eigval_tmp(N_states))
|
||||
dressing_dets = 0.d0
|
||||
do i = 1, N_det
|
||||
call htilde_mu_mat_bi_ortho_tot_slow(psi_det(1,1,i), psi_det(1,1,i), N_int, H_jj(i))
|
||||
call get_excitation_degree(HF_bitmask,psi_det(1,1,i),degree,N_int)
|
||||
if(degree == 1 .or. degree == 2)then
|
||||
call htilde_mu_mat_opt_bi_ortho(HF_bitmask,psi_det(1,1,i),N_int,hmono,htwoe,hthree,h0j(i))
|
||||
endif
|
||||
enddo
|
||||
reigvec_tc_bi_orth_tmp = 0.d0
|
||||
do i = 1, N_det
|
||||
reigvec_tc_bi_orth_tmp(i,1) = psi_r_coef_bi_ortho(i,1)
|
||||
enddo
|
||||
vec_tmp = 0.d0
|
||||
do istate = 1, N_states
|
||||
vec_tmp(:,istate) = reigvec_tc_bi_orth_tmp(:,istate)
|
||||
enddo
|
||||
do istate = N_states+1, n_states_diag
|
||||
vec_tmp(istate,istate) = 1.d0
|
||||
enddo
|
||||
print*,'Diagonalizing the TC CISD '
|
||||
call davidson_general_diag_dressed_ext_rout_nonsym_b1space(vec_tmp, H_jj, dressing_dets,eigval_tmp, N_det, n_states, n_states_diag, converged, htc_bi_ortho_calc_tdav_slow)
|
||||
do i = 1, N_det
|
||||
e_corr_dets(i) = reigvec_tc_bi_orth_tmp(i,1) * h0j(i)/reigvec_tc_bi_orth_tmp(1,1)
|
||||
enddo
|
||||
E_before = eigval_tmp(1)
|
||||
print*,'Starting from ',E_before
|
||||
|
||||
e_current = 10.d0
|
||||
thr = 1.d-5
|
||||
it = 0
|
||||
dressing_dets = 0.d0
|
||||
double precision, allocatable :: H_jj(:),vec_tmp(:,:),eigval_tmp(:)
|
||||
external htc_bi_ortho_calc_tdav_slow
|
||||
external htcdag_bi_ortho_calc_tdav_slow
|
||||
logical :: converged
|
||||
do while (dabs(E_before-E_current).gt.thr)
|
||||
it += 1
|
||||
E_before = E_current
|
||||
! h_sc2 = htilde_matrix_elmt_bi_ortho
|
||||
call get_cisd_sc2_dressing(psi_det,e_corr_dets,N_det,dressing_dets)
|
||||
do i = 1, N_det
|
||||
! print*,'dressing_dets(i) = ',dressing_dets(i)
|
||||
h_sc2(i,i) += dressing_dets(i)
|
||||
enddo
|
||||
print*,'********************'
|
||||
print*,'iteration ',it
|
||||
! call non_hrmt_real_diag(N_det,h_sc2,&
|
||||
! leigvec_tc_bi_orth_tmp,reigvec_tc_bi_orth_tmp,&
|
||||
! n_real,eigval_right_tmp)
|
||||
! print*,'eigval_right_tmp(1)',eigval_right_tmp(1)
|
||||
vec_tmp = 0.d0
|
||||
do istate = 1, N_states
|
||||
vec_tmp(:,istate) = reigvec_tc_bi_orth_tmp(:,istate)
|
||||
enddo
|
||||
do istate = N_states+1, n_states_diag
|
||||
vec_tmp(istate,istate) = 1.d0
|
||||
enddo
|
||||
call davidson_general_diag_dressed_ext_rout_nonsym_b1space(vec_tmp, H_jj, dressing_dets,eigval_tmp, N_det, n_states, n_states_diag, converged, htc_bi_ortho_calc_tdav_slow)
|
||||
print*,'outside Davidson'
|
||||
print*,'eigval_tmp(1) = ',eigval_tmp(1)
|
||||
do i = 1, N_det
|
||||
reigvec_tc_bi_orth_tmp(i,1) = vec_tmp(i,1)
|
||||
e_corr_dets(i) = reigvec_tc_bi_orth_tmp(i,1) * h0j(i)/reigvec_tc_bi_orth_tmp(1,1)
|
||||
enddo
|
||||
! E_current = eigval_right_tmp(1)
|
||||
E_current = eigval_tmp(1)
|
||||
print*,'it, E(SC)^2 = ',it,E_current
|
||||
enddo
|
||||
eigval_tc_cisd_sc2_bi_ortho(1:N_states) = eigval_right_tmp(1:N_states)
|
||||
reigvec_tc_cisd_sc2_bi_ortho(1:N_det,1:N_states) = reigvec_tc_bi_orth_tmp(1:N_det,1:N_states)
|
||||
leigvec_tc_cisd_sc2_bi_ortho(1:N_det,1:N_states) = leigvec_tc_bi_orth_tmp(1:N_det,1:N_states)
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
subroutine get_cisd_sc2_dressing(dets,e_corr_dets,ndet,dressing_dets)
|
||||
implicit none
|
||||
use bitmasks
|
||||
integer, intent(in) :: ndet
|
||||
integer(bit_kind), intent(in) :: dets(N_int,2,ndet)
|
||||
double precision, intent(in) :: e_corr_dets(ndet)
|
||||
double precision, intent(out) :: dressing_dets(ndet)
|
||||
integer, allocatable :: degree(:),hole(:,:),part(:,:),spin(:,:)
|
||||
integer(bit_kind), allocatable :: hole_part(:,:,:)
|
||||
integer :: i,j,k, exc(0:2,2,2),h1,p1,h2,p2,s1,s2
|
||||
integer(bit_kind) :: xorvec(2,N_int)
|
||||
|
||||
double precision :: phase
|
||||
dressing_dets = 0.d0
|
||||
allocate(degree(ndet),hole(2,ndet),part(2,ndet), spin(2,ndet),hole_part(N_int,2,ndet))
|
||||
do i = 2, ndet
|
||||
call get_excitation_degree(HF_bitmask,dets(1,1,i),degree(i),N_int)
|
||||
do j = 1, N_int
|
||||
hole_part(j,1,i) = xor( HF_bitmask(j,1), dets(j,1,i))
|
||||
hole_part(j,2,i) = xor( HF_bitmask(j,2), dets(j,2,i))
|
||||
enddo
|
||||
if(degree(i) == 1)then
|
||||
call get_single_excitation(HF_bitmask,psi_det(1,1,i),exc,phase,N_int)
|
||||
else if(degree(i) == 2)then
|
||||
call get_double_excitation(HF_bitmask,psi_det(1,1,i),exc,phase,N_int)
|
||||
endif
|
||||
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
|
||||
hole(1,i) = h1
|
||||
hole(2,i) = h2
|
||||
part(1,i) = p1
|
||||
part(2,i) = p2
|
||||
spin(1,i) = s1
|
||||
spin(2,i) = s2
|
||||
enddo
|
||||
|
||||
integer :: same
|
||||
if(elec_alpha_num+elec_beta_num<3)return
|
||||
do i = 2, ndet
|
||||
do j = i+1, ndet
|
||||
same = 0
|
||||
if(degree(i) == degree(j) .and. degree(i)==1)cycle
|
||||
do k = 1, N_int
|
||||
xorvec(k,1) = iand(hole_part(k,1,i),hole_part(k,1,j))
|
||||
xorvec(k,2) = iand(hole_part(k,2,i),hole_part(k,2,j))
|
||||
same += popcnt(xorvec(k,1)) + popcnt(xorvec(k,2))
|
||||
enddo
|
||||
! print*,'i,j',i,j
|
||||
! call debug_det(dets(1,1,i),N_int)
|
||||
! call debug_det(hole_part(1,1,i),N_int)
|
||||
! call debug_det(dets(1,1,j),N_int)
|
||||
! call debug_det(hole_part(1,1,j),N_int)
|
||||
! print*,'same = ',same
|
||||
if(same.eq.0)then
|
||||
dressing_dets(i) += e_corr_dets(j)
|
||||
dressing_dets(j) += e_corr_dets(i)
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
|
||||
end
|
||||
@ -1,170 +0,0 @@
|
||||
|
||||
! ---
|
||||
|
||||
program test_tc
|
||||
|
||||
implicit none
|
||||
|
||||
my_grid_becke = .True.
|
||||
PROVIDE tc_grid1_a tc_grid1_r
|
||||
my_n_pt_r_grid = tc_grid1_r
|
||||
my_n_pt_a_grid = tc_grid1_a
|
||||
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
|
||||
|
||||
read_wf = .True.
|
||||
touch read_wf
|
||||
|
||||
call provide_all_three_ints_bi_ortho()
|
||||
call routine_h_triple_left
|
||||
call routine_h_triple_right
|
||||
! call routine_test_s2_davidson
|
||||
|
||||
end
|
||||
|
||||
subroutine routine_h_triple_right
|
||||
implicit none
|
||||
logical :: do_right
|
||||
integer :: sze ,i, N_st, j
|
||||
double precision :: sij, accu_e, accu_s, accu_e_0, accu_s_0
|
||||
double precision, allocatable :: v_0_ref(:,:),u_0(:,:),s_0_ref(:,:)
|
||||
double precision, allocatable :: v_0_new(:,:),s_0_new(:,:)
|
||||
sze = N_det
|
||||
N_st = 1
|
||||
allocate(v_0_ref(N_det,1),u_0(N_det,1),s_0_ref(N_det,1),s_0_new(N_det,1),v_0_new(N_det,1))
|
||||
print*,'Checking first the Right '
|
||||
do i = 1, sze
|
||||
u_0(i,1) = psi_r_coef_bi_ortho(i,1)
|
||||
enddo
|
||||
double precision :: wall0,wall1
|
||||
call wall_time(wall0)
|
||||
call H_tc_s2_u_0_with_pure_three_omp(v_0_ref,s_0_ref, u_0,N_st,sze)
|
||||
call wall_time(wall1)
|
||||
print*,'time for omp',wall1 - wall0
|
||||
call wall_time(wall0)
|
||||
call H_tc_s2_u_0_with_pure_three(v_0_new, s_0_new, u_0, N_st, sze)
|
||||
call wall_time(wall1)
|
||||
print*,'time serial ',wall1 - wall0
|
||||
accu_e = 0.d0
|
||||
accu_s = 0.d0
|
||||
do i = 1, sze
|
||||
accu_e += dabs(v_0_ref(i,1) - v_0_new(i,1))
|
||||
accu_s += dabs(s_0_ref(i,1) - s_0_new(i,1))
|
||||
enddo
|
||||
print*,'accu_e = ',accu_e
|
||||
print*,'accu_s = ',accu_s
|
||||
|
||||
end
|
||||
|
||||
subroutine routine_h_triple_left
|
||||
implicit none
|
||||
logical :: do_right
|
||||
integer :: sze ,i, N_st, j
|
||||
double precision :: sij, accu_e, accu_s, accu_e_0, accu_s_0
|
||||
double precision, allocatable :: v_0_ref(:,:),u_0(:,:),s_0_ref(:,:)
|
||||
double precision, allocatable :: v_0_new(:,:),s_0_new(:,:)
|
||||
sze = N_det
|
||||
N_st = 1
|
||||
allocate(v_0_ref(N_det,1),u_0(N_det,1),s_0_ref(N_det,1),s_0_new(N_det,1),v_0_new(N_det,1))
|
||||
print*,'Checking the Left '
|
||||
do i = 1, sze
|
||||
u_0(i,1) = psi_l_coef_bi_ortho(i,1)
|
||||
enddo
|
||||
double precision :: wall0,wall1
|
||||
call wall_time(wall0)
|
||||
call H_tc_s2_dagger_u_0_with_pure_three_omp(v_0_ref,s_0_ref, u_0,N_st,sze)
|
||||
call wall_time(wall1)
|
||||
print*,'time for omp',wall1 - wall0
|
||||
call wall_time(wall0)
|
||||
call H_tc_s2_dagger_u_0_with_pure_three(v_0_new, s_0_new, u_0, N_st, sze)
|
||||
call wall_time(wall1)
|
||||
print*,'time serial ',wall1 - wall0
|
||||
accu_e = 0.d0
|
||||
accu_s = 0.d0
|
||||
do i = 1, sze
|
||||
accu_e += dabs(v_0_ref(i,1) - v_0_new(i,1))
|
||||
accu_s += dabs(s_0_ref(i,1) - s_0_new(i,1))
|
||||
enddo
|
||||
print*,'accu_e = ',accu_e
|
||||
print*,'accu_s = ',accu_s
|
||||
|
||||
end
|
||||
|
||||
|
||||
subroutine routine_test_s2_davidson
|
||||
implicit none
|
||||
double precision, allocatable :: H_jj(:),vec_tmp(:,:), energies(:) , s2(:)
|
||||
integer :: i,istate
|
||||
logical :: converged
|
||||
external H_tc_s2_dagger_u_0_opt
|
||||
external H_tc_s2_u_0_opt
|
||||
allocate(H_jj(N_det),vec_tmp(N_det,n_states_diag),energies(n_states_diag), s2(n_states_diag))
|
||||
do i = 1, N_det
|
||||
call htilde_mu_mat_bi_ortho_tot_slow(psi_det(1,1,i), psi_det(1,1,i), N_int, H_jj(i))
|
||||
enddo
|
||||
! Preparing the left-eigenvector
|
||||
print*,'Computing the left-eigenvector '
|
||||
vec_tmp = 0.d0
|
||||
do istate = 1, N_states
|
||||
vec_tmp(1:N_det,istate) = psi_l_coef_bi_ortho(1:N_det,istate)
|
||||
enddo
|
||||
do istate = N_states+1, n_states_diag
|
||||
vec_tmp(istate,istate) = 1.d0
|
||||
enddo
|
||||
do istate = 1, N_states
|
||||
leigvec_tc_bi_orth(1:N_det,istate) = vec_tmp(1:N_det,istate)
|
||||
enddo
|
||||
integer :: n_it_max
|
||||
n_it_max = 1
|
||||
call davidson_hs2_nonsym_b1space(vec_tmp, H_jj, s2, energies, N_det, n_states, n_states_diag, n_it_max, converged, H_tc_s2_dagger_u_0_opt)
|
||||
double precision, allocatable :: v_0_new(:,:),s_0_new(:,:)
|
||||
integer :: sze,N_st
|
||||
logical :: do_right
|
||||
sze = N_det
|
||||
N_st = 1
|
||||
do_right = .False.
|
||||
allocate(s_0_new(N_det,1),v_0_new(N_det,1))
|
||||
call H_tc_s2_u_0_nstates_openmp(v_0_new,s_0_new,vec_tmp,N_st,sze, do_right)
|
||||
double precision :: accu_e_0, accu_s_0
|
||||
accu_e_0 = 0.d0
|
||||
accu_s_0 = 0.d0
|
||||
do i = 1, sze
|
||||
accu_e_0 += v_0_new(i,1) * vec_tmp(i,1)
|
||||
accu_s_0 += s_0_new(i,1) * vec_tmp(i,1)
|
||||
enddo
|
||||
print*,'energies = ',energies
|
||||
print*,'s2 = ',s2
|
||||
print*,'accu_e_0',accu_e_0
|
||||
print*,'accu_s_0',accu_s_0
|
||||
|
||||
! Preparing the right-eigenvector
|
||||
print*,'Computing the right-eigenvector '
|
||||
vec_tmp = 0.d0
|
||||
do istate = 1, N_states
|
||||
vec_tmp(1:N_det,istate) = psi_r_coef_bi_ortho(1:N_det,istate)
|
||||
enddo
|
||||
do istate = N_states+1, n_states_diag
|
||||
vec_tmp(istate,istate) = 1.d0
|
||||
enddo
|
||||
do istate = 1, N_states
|
||||
leigvec_tc_bi_orth(1:N_det,istate) = vec_tmp(1:N_det,istate)
|
||||
enddo
|
||||
n_it_max = 1
|
||||
call davidson_hs2_nonsym_b1space(vec_tmp, H_jj, s2, energies, N_det, n_states, n_states_diag, n_it_max, converged, H_tc_s2_u_0_opt)
|
||||
sze = N_det
|
||||
N_st = 1
|
||||
do_right = .True.
|
||||
v_0_new = 0.d0
|
||||
s_0_new = 0.d0
|
||||
call H_tc_s2_u_0_nstates_openmp(v_0_new,s_0_new,vec_tmp,N_st,sze, do_right)
|
||||
accu_e_0 = 0.d0
|
||||
accu_s_0 = 0.d0
|
||||
do i = 1, sze
|
||||
accu_e_0 += v_0_new(i,1) * vec_tmp(i,1)
|
||||
accu_s_0 += s_0_new(i,1) * vec_tmp(i,1)
|
||||
enddo
|
||||
print*,'energies = ',energies
|
||||
print*,'s2 = ',s2
|
||||
print*,'accu_e_0',accu_e_0
|
||||
print*,'accu_s_0',accu_s_0
|
||||
|
||||
end
|
||||
Loading…
Reference in New Issue
Block a user