10
0
mirror of https://github.com/LCPQ/quantum_package synced 2024-12-25 05:43:47 +01:00

Forgot files

This commit is contained in:
Anthony Scemama 2018-02-09 18:23:05 +01:00
parent d8b94c0473
commit 0f75e345c0
8 changed files with 473 additions and 0 deletions

View File

@ -0,0 +1,32 @@
program davidson_slave
use f77_zmq
implicit none
integer(ZMQ_PTR), external :: new_zmq_to_qp_run_socket
integer(ZMQ_PTR) :: zmq_to_qp_run_socket
double precision :: energy(N_states_diag)
character*(64) :: state
call provide_everything
call switch_qp_run_to_master
call omp_set_nested(.True.)
zmq_context = f77_zmq_ctx_new ()
zmq_state = 'davidson'
state = 'Waiting'
zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
do
call wait_for_state(zmq_state,state)
if(trim(state) /= "davidson") exit
integer :: rc, i
print *, 'Davidson slave running'
call davidson_slave_tcp(i)
end do
end
subroutine provide_everything
PROVIDE mo_bielec_integrals_in_map psi_det_sorted_bit N_states_diag zmq_context ref_bitmask_energy
end subroutine

View File

@ -0,0 +1,16 @@
program diag_and_save
implicit none
read_wf = .True.
touch read_wf
call routine
end
subroutine routine
implicit none
call diagonalize_CI
print*,'N_det = ',N_det
call save_wavefunction_general(N_det,N_states_diag,psi_det_sorted,size(psi_coef_sorted,1),psi_coef_sorted)
end

View File

@ -0,0 +1,29 @@
program diag_and_save
implicit none
read_wf = .True.
touch read_wf
call routine
end
subroutine routine
implicit none
print*,'N_det = ',N_det
PROVIDE H_apply_buffer_allocated
if (s2_eig) then
call make_s2_eigenfunction
endif
call diagonalize_CI
integer :: igood_state
igood_state=1
double precision, allocatable :: psi_coef_tmp(:)
allocate(psi_coef_tmp(n_det))
integer :: i
do i = 1, N_det
psi_coef_tmp(i) = psi_coef(i,igood_state)
enddo
call save_wavefunction_general(N_det,1,psi_det,n_det,psi_coef_tmp)
deallocate(psi_coef_tmp)
end

View File

@ -0,0 +1,26 @@
program diag_and_save
implicit none
read_wf = .True.
touch read_wf
call routine
end
subroutine routine
implicit none
print*,'N_det = ',N_det
call diagonalize_CI
write(*,*)'Which state would you like to save ?'
integer :: igood_state
read(5,*)igood_state
double precision, allocatable :: psi_coef_tmp(:)
allocate(psi_coef_tmp(n_det))
integer :: i
do i = 1, N_det
psi_coef_tmp(i) = psi_coef(i,igood_state)
enddo
call save_wavefunction_general(N_det,1,psi_det,n_det,psi_coef_tmp)
deallocate(psi_coef_tmp)
end

View File

@ -0,0 +1,162 @@
program first_guess
use bitmasks
implicit none
BEGIN_DOC
! Select all the determinants with the lowest energy as a starting point.
END_DOC
integer :: i,j
double precision, allocatable :: orb_energy(:)
double precision :: E
integer, allocatable :: kept(:)
integer :: nelec_kept(2)
character :: occ_char, keep_char
PROVIDE H_apply_buffer_allocated psi_det
allocate (orb_energy(mo_tot_num), kept(0:mo_tot_num))
nelec_kept(1:2) = 0
kept(0) = 0
print *, 'Orbital energies'
print *, '================'
print *, ''
do i=1,mo_tot_num
keep_char = ' '
occ_char = '-'
orb_energy(i) = mo_mono_elec_integral(i,i)
do j=1,elec_beta_num
if (i==j) cycle
orb_energy(i) += mo_bielec_integral_jj_anti(i,j)
enddo
do j=1,elec_alpha_num
orb_energy(i) += mo_bielec_integral_jj(i,j)
enddo
if ( (orb_energy(i) > -.5d0).and.(orb_energy(i) < .1d0) ) then
kept(0) += 1
keep_char = 'X'
kept( kept(0) ) = i
if (i <= elec_beta_num) then
nelec_kept(2) += 1
endif
if (i <= elec_alpha_num) then
nelec_kept(1) += 1
endif
endif
if (i <= elec_alpha_num) then
if (i <= elec_beta_num) then
occ_char = '#'
else
occ_char = '+'
endif
endif
print '(I4, 3X, A, 3X, F10.6, 3X, A)', i, occ_char, orb_energy(i), keep_char
enddo
integer, allocatable :: list (:,:)
integer(bit_kind), allocatable :: string(:,:)
allocate ( list(N_int*bit_kind_size,2), string(N_int,2) )
string = ref_bitmask
call bitstring_to_list( string(1,1), list(1,1), elec_alpha_num, N_int)
call bitstring_to_list( string(1,2), list(1,2), elec_beta_num , N_int)
psi_det_alpha_unique(:,1) = string(:,1)
psi_det_beta_unique (:,1) = string(:,2)
N_det_alpha_unique = 1
N_det_beta_unique = 1
integer :: i1,i2,i3,i4,i5,i6,i7,i8,i9
psi_det_size = kept(0)**(nelec_kept(1)+nelec_kept(2))
print *, kept(0), nelec_kept(:)
call write_int(6,psi_det_size,'psi_det_size')
TOUCH psi_det_size
BEGIN_SHELL [ /usr/bin/python ]
template_alpha_ext = """
do %(i2)s = %(i1)s-1,1,-1
list(elec_alpha_num-%(i)d,1) = kept(%(i2)s)
call list_to_bitstring( string(1,1), list(1,1), elec_alpha_num, N_int)
"""
template_alpha = """
do %(i2)s = %(i1)s-1,1,-1
list(elec_alpha_num-%(i)d,1) = kept(%(i2)s)
call list_to_bitstring( string(1,1), list(1,1), elec_alpha_num, N_int)
N_det_alpha_unique += 1
psi_det_alpha_unique(:,N_det_alpha_unique) = string(:,1)
"""
template_beta_ext = """
do %(i2)s = %(i1)s-1,1,-1
list(elec_beta_num-%(i)d,2) = kept(%(i2)s)
call list_to_bitstring( string(1,2), list(1,2), elec_beta_num, N_int)
"""
template_beta = """
do %(i2)s = %(i1)s-1,1,-1
list(elec_beta_num-%(i)d,2) = kept(%(i2)s)
call list_to_bitstring( string(1,2), list(1,2), elec_beta_num, N_int)
N_det_beta_unique += 1
psi_det_beta_unique(:,N_det_beta_unique) = string(:,2)
"""
def write(template_ext,template,imax):
print "case(%d)"%(imax)
def aux(i2,i1,i,j):
if (i==imax-1):
print template%locals()
else:
print template_ext%locals()
i += 1
j -= 1
if (i != imax):
i1 = "i%d"%(i)
i2 = "i%d"%(i+1)
aux(i2,i1,i,j)
print "enddo"
i2 = "i1"
i1 = "kept(0)+1"
i = 0
aux(i2,i1,i,imax)
def main():
print """
select case (nelec_kept(1))
case(0)
continue
"""
for imax in range(1,10):
write(template_alpha_ext,template_alpha,imax)
print """
end select
select case (nelec_kept(2))
case(0)
continue
"""
for imax in range(1,10):
write(template_beta_ext,template_beta,imax)
print "end select"
main()
END_SHELL
TOUCH N_det_alpha_unique N_det_beta_unique psi_det_alpha_unique psi_det_beta_unique
call create_wf_of_psi_bilinear_matrix(.False.)
call diagonalize_ci
j= N_det
do i=1,N_det
if (psi_average_norm_contrib_sorted(i) < 1.d-6) then
j = i-1
exit
endif
! call debug_det(psi_det_sorted(1,1,i),N_int)
enddo
call save_wavefunction_general(j,N_states,psi_det_sorted,size(psi_coef_sorted,1),psi_coef_sorted)
deallocate(orb_energy, kept, list, string)
end

View File

@ -0,0 +1,10 @@
BEGIN_PROVIDER [ double precision, dressing_column_h, (N_det,N_states) ]
&BEGIN_PROVIDER [ double precision, dressing_column_s, (N_det,N_states) ]
implicit none
BEGIN_DOC
! Null dressing vectors
END_DOC
dressing_column_h(:,:) = 0.d0
dressing_column_s(:,:) = 0.d0
END_PROVIDER

View File

@ -0,0 +1,176 @@
program print_H_matrix_restart
implicit none
read_wf = .True.
touch read_wf
call routine
end
subroutine routine
use bitmasks
implicit none
integer :: i,j
integer, allocatable :: H_matrix_degree(:,:)
double precision, allocatable :: H_matrix_phase(:,:)
integer :: degree
integer(bit_kind), allocatable :: keys_tmp(:,:,:)
allocate(keys_tmp(N_int,2,N_det))
do i = 1, N_det
print*,''
call debug_det(psi_det(1,1,i),N_int)
do j = 1, N_int
keys_tmp(j,1,i) = psi_det(j,1,i)
keys_tmp(j,2,i) = psi_det(j,2,i)
enddo
enddo
if(N_det.ge.10000)then
print*,'Warning !!!'
print*,'Number of determinants is ',N_det
print*,'It means that the H matrix will be enormous !'
print*,'stoppping ..'
stop
endif
print*,''
print*,'Determinants '
do i = 1, N_det
enddo
allocate(H_matrix_degree(N_det,N_det),H_matrix_phase(N_det,N_det))
integer :: exc(0:2,2,2)
double precision :: phase
do i = 1, N_det
do j = i, N_det
call get_excitation_degree(psi_det(1,1,i),psi_det(1,1,j),degree,N_int)
H_matrix_degree(i,j) = degree
H_matrix_degree(j,i) = degree
phase = 0.d0
if(degree==1.or.degree==2)then
call get_excitation(psi_det(1,1,i),psi_det(1,1,j),exc,degree,phase,N_int)
endif
H_matrix_phase(i,j) = phase
H_matrix_phase(j,i) = phase
enddo
enddo
print*,'H matrix '
double precision :: ref_h_matrix,s2
ref_h_matrix = H_matrix_all_dets(1,1)
print*,'HF like determinant energy = ',ref_bitmask_energy+nuclear_repulsion
print*,'Ref element of H_matrix = ',ref_h_matrix+nuclear_repulsion
print*,'Printing the H matrix ...'
print*,''
print*,''
!do i = 1, N_det
! H_matrix_all_dets(i,i) -= ref_h_matrix
!enddo
do i = 1, N_det
H_matrix_all_dets(i,i) += nuclear_repulsion
enddo
!do i = 5,N_det
! H_matrix_all_dets(i,3) = 0.d0
! H_matrix_all_dets(3,i) = 0.d0
! H_matrix_all_dets(i,4) = 0.d0
! H_matrix_all_dets(4,i) = 0.d0
!enddo
do i = 1, N_det
write(*,'(I3,X,A3,1000(F16.7))')i,' | ',H_matrix_all_dets(i,:)
enddo
print*,''
print*,''
print*,''
print*,'Printing the degree of excitations within the H matrix'
print*,''
print*,''
do i = 1, N_det
write(*,'(I3,X,A3,X,1000(I1,X))')i,' | ',H_matrix_degree(i,:)
enddo
print*,''
print*,''
print*,'Printing the phase of the Hamiltonian matrix elements '
print*,''
print*,''
do i = 1, N_det
write(*,'(I3,X,A3,X,1000(F3.0,X))')i,' | ',H_matrix_phase(i,:)
enddo
print*,''
double precision, allocatable :: eigenvectors(:,:), eigenvalues(:)
double precision, allocatable :: s2_eigvalues(:)
allocate (eigenvectors(size(H_matrix_all_dets,1),N_det))
allocate (eigenvalues(N_det),s2_eigvalues(N_det))
call lapack_diag(eigenvalues,eigenvectors, &
H_matrix_all_dets,size(H_matrix_all_dets,1),N_det)
print*,'Two first eigenvectors '
call u_0_S2_u_0(s2_eigvalues,eigenvectors,n_det,keys_tmp,N_int,N_det,size(eigenvectors,1))
do j =1, N_states
print*,'s2 = ',s2_eigvalues(j)
print*,'e = ',eigenvalues(j)
print*,'coefs : '
do i = 1, N_det
print*,'i = ',i,eigenvectors(i,j)
enddo
if(j>1)then
print*,'Delta E(H) = ',eigenvalues(1) - eigenvalues(j)
print*,'Delta E(eV) = ',(eigenvalues(1) - eigenvalues(j))*27.2114d0
endif
enddo
double precision :: get_mo_bielec_integral,k_a_iv,k_b_iv
integer :: h1,p1,h2,p2
h1 = 10
p1 = 16
h2 = 14
p2 = 14
!h1 = 1
!p1 = 4
!h2 = 2
!p2 = 2
k_a_iv = get_mo_bielec_integral(h1,h2,p2,p1,mo_integrals_map)
h2 = 15
p2 = 15
k_b_iv = get_mo_bielec_integral(h1,h2,p2,p1,mo_integrals_map)
print*,'k_a_iv = ',k_a_iv
print*,'k_b_iv = ',k_b_iv
double precision :: k_av,k_bv,k_ai,k_bi
h1 = 16
p1 = 14
h2 = 14
p2 = 16
k_av = get_mo_bielec_integral(h1,h2,p1,p2,mo_integrals_map)
h1 = 16
p1 = 15
h2 = 15
p2 = 16
k_bv = get_mo_bielec_integral(h1,h2,p1,p2,mo_integrals_map)
h1 = 10
p1 = 14
h2 = 14
p2 = 10
k_ai = get_mo_bielec_integral(h1,h2,p1,p2,mo_integrals_map)
h1 = 10
p1 = 15
h2 = 15
p2 = 10
k_bi = get_mo_bielec_integral(h1,h2,p1,p2,mo_integrals_map)
print*,'k_av, k_bv = ',k_av,k_bv
print*,'k_ai, k_bi = ',k_ai,k_bi
double precision :: k_iv
h1 = 10
p1 = 16
h2 = 16
p2 = 10
k_iv = get_mo_bielec_integral(h1,h2,p1,p2,mo_integrals_map)
print*,'k_iv = ',k_iv
end

View File

@ -0,0 +1,22 @@
program print_energy
implicit none
read_wf = .true.
touch read_wf
call routine
end
subroutine routine
implicit none
integer :: i,j
double precision :: accu,hij
print*, 'psi_energy = ',psi_energy + nuclear_repulsion
accu = 0.d0
! do i = 1,N_det
! do j = 1,N_det
! call i_H_j(psi_det(1,1,j),psi_det(1,1,i),N_int,hij)
! accu += psi_coef(i,1) * psi_coef(j,1) * hij
! enddo
! enddo
! print*, 'accu = ',accu + nuclear_repulsion
end