10
0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-06-26 15:12:19 +02:00

added fci_tc_bi

This commit is contained in:
eginer 2022-10-05 17:59:23 +02:00
parent ad0203c959
commit 70516c8c05
17 changed files with 793 additions and 103 deletions

View File

@ -16,108 +16,6 @@ program bi_ort_ints
! call test_3_e_diag
! call test_3_e_diag_cycle1
! call test_3_e
call routine_test_one_int
! call routine_test_one_int
end
subroutine routine_test_one_int
implicit none
integer :: p,q,r,s,ii
integer :: i,j
i = 3
j = 5
double precision :: accu
double precision, allocatable :: vec(:)
integer, allocatable :: iorder(:)
allocate(vec(ao_num**4),iorder(ao_num**4))
accu = 0.d0
ii = 0
do p = 1, ao_num !
do q = 1, ao_num
do r = 1, ao_num
do s = 1, ao_num
!<ji | ji>
!
! j j i i
if(dabs(mo_l_coef(s,j) * mo_l_coef(q,i) * ao_two_e_tc_tot(s,r,q,p) * mo_r_coef(p,i) * mo_r_coef(r,j)).gt.10)then
write(33,'(3(F16.10,X),4(I3,X))')mo_l_coef(s,j) * mo_l_coef(q,i)* mo_r_coef(p,i) * mo_r_coef(r,j) , ao_two_e_tc_tot(s,r,q,p), mo_l_coef(s,j) * mo_l_coef(q,i) * ao_two_e_tc_tot(s,r,q,p) * mo_r_coef(p,i) * mo_r_coef(r,j) , s,q,p,r
endif
ii += 1
iorder(ii) = ii
vec(ii) = mo_l_coef(s,j) * mo_l_coef(q,i) * ao_two_e_tc_tot(s,r,q,p) * mo_r_coef(p,i) * mo_r_coef(r,j)
accu += mo_l_coef(s,j) * mo_l_coef(q,i) * ao_two_e_tc_tot(s,r,q,p) * mo_r_coef(p,i) * mo_r_coef(r,j)
enddo
enddo
enddo
enddo
call dsort(vec,iorder,ao_num**4)
accu = 0.d0
do i = 1, ao_num**4
accu += vec(i)
write(34,*)i,vec(i),accu
enddo
print*,'accu = ',accu
end
subroutine routine_twoe
implicit none
integer :: i,j,k,l
double precision :: old, get_mo_two_e_integral_tc_int
double precision :: ref,new, accu, contrib, bi_ortho_mo_ints
accu = 0.d0
print*,'Testing the bi ortho two e'
do j = 1, mo_num
do i = 1, mo_num
do l = 1, mo_num
do k = 1, mo_num
! mo_non_hermit_term(k,l,i,j) = <k l| V(r_12) |i j>
! ref = bi_ortho_mo_ints(k,l,i,j)
ref = bi_ortho_mo_ints(l,k,j,i)
new = mo_bi_ortho_tc_two_e(l,k,j,i)
! old = get_mo_two_e_integral_tc_int(k,l,i,j,mo_integrals_tc_int_map)
! old += mo_non_hermit_term(l,k,j,i)
contrib = dabs(ref - new)
if(dabs(ref).gt.1.d-10)then
if(contrib.gt.1.d-10)then
print*,k,l,i,j
print*,ref,new,contrib
endif
endif
accu += contrib
enddo
enddo
enddo
enddo
print*,'accu = ',accu/(dble(mo_num)**4)
end
subroutine routine_onee
implicit none
integer :: i,k
double precision :: ref,new,accu,contrib
print*,'Testing the bi ortho one e'
accu = 0.d0
do i = 1, mo_num
do k = 1, mo_num
ref = mo_bi_ortho_tc_one_e_slow(k,i)
new = mo_bi_ortho_tc_one_e(k,i)
contrib = dabs(ref - new)
if(dabs(ref).gt.1.d-10)then
if(contrib .gt. 1.d-10)then
print*,'i,k',i,k
print*,ref,new,contrib
endif
endif
accu += contrib
enddo
enddo
print*,'accu = ',accu/mo_num**2
end

17
src/fci_tc_bi/EZFIO.cfg Normal file
View File

@ -0,0 +1,17 @@
[energy]
type: double precision
doc: Calculated Selected |FCI| energy
interface: ezfio
size: (determinants.n_states)
[energy_pt2]
type: double precision
doc: Calculated |FCI| energy + |PT2|
interface: ezfio
size: (determinants.n_states)
[cipsi_tc]
type: character*(32)
doc: TODO
interface: ezfio,provider,ocaml
default: h_tc

3
src/fci_tc_bi/NEED Normal file
View File

@ -0,0 +1,3 @@
tc_bi_ortho
davidson_undressed
cipsi_tc_bi_ortho

12
src/fci_tc_bi/class.irp.f Normal file
View File

@ -0,0 +1,12 @@
BEGIN_PROVIDER [ logical, do_only_1h1p ]
&BEGIN_PROVIDER [ logical, do_only_cas ]
&BEGIN_PROVIDER [ logical, do_ddci ]
implicit none
BEGIN_DOC
! In the FCI case, all those are always false
END_DOC
do_only_1h1p = .False.
do_only_cas = .False.
do_ddci = .False.
END_PROVIDER

215
src/fci_tc_bi/copy_wf.irp.f Normal file
View File

@ -0,0 +1,215 @@
use bitmasks
subroutine copy_H_apply_buffer_to_wf_tc
use omp_lib
implicit none
BEGIN_DOC
! Copies the H_apply buffer to psi_coef.
! After calling this subroutine, N_det, psi_det and psi_coef need to be touched
END_DOC
integer(bit_kind), allocatable :: buffer_det(:,:,:)
double precision, allocatable :: buffer_r_coef(:,:), buffer_l_coef(:,:)
integer :: i,j,k
integer :: N_det_old
PROVIDE H_apply_buffer_allocated
ASSERT (N_int > 0)
ASSERT (N_det > 0)
allocate ( buffer_det(N_int,2,N_det), buffer_r_coef(N_det,N_states), buffer_l_coef(N_det,N_states) )
! Backup determinants
j=0
do i=1,N_det
! if (pruned(i)) cycle ! Pruned determinants
j+=1
ASSERT (sum(popcnt(psi_det(:,1,i))) == elec_alpha_num)
ASSERT (sum(popcnt(psi_det(:,2,i))) == elec_beta_num)
buffer_det(:,:,j) = psi_det(:,:,i)
enddo
N_det_old = j
! Backup coefficients
do k=1,N_states
j=0
do i=1,N_det
! if (pruned(i)) cycle ! Pruned determinants
j += 1
buffer_r_coef(j,k) = psi_r_coef_bi_ortho(i,k)
buffer_l_coef(j,k) = psi_l_coef_bi_ortho(i,k)
enddo
ASSERT ( j == N_det_old )
enddo
! Update N_det
N_det = N_det_old
do j=0,nproc-1
N_det = N_det + H_apply_buffer(j)%N_det
enddo
! Update array sizes
if (psi_det_size < N_det) then
psi_det_size = N_det
TOUCH psi_det_size
endif
! Restore backup in resized array
do i=1,N_det_old
psi_det(:,:,i) = buffer_det(:,:,i)
ASSERT (sum(popcnt(psi_det(:,1,i))) == elec_alpha_num)
ASSERT (sum(popcnt(psi_det(:,2,i))) == elec_beta_num )
enddo
do k=1,N_states
do i=1,N_det_old
psi_r_coef_bi_ortho(i,k) = buffer_r_coef(i,k)
psi_l_coef_bi_ortho(i,k) = buffer_l_coef(i,k)
enddo
enddo
! Copy new buffers
!$OMP PARALLEL DEFAULT(SHARED) &
!$OMP PRIVATE(j,k,i) FIRSTPRIVATE(N_det_old) &
!$OMP SHARED(N_int,H_apply_buffer,psi_det,psi_r_coef_bi_ortho,psi_l_coef_bi_ortho,N_states,psi_det_size)
j=0
!$ j=omp_get_thread_num()
do k=0,j-1
N_det_old += H_apply_buffer(k)%N_det
enddo
do i=1,H_apply_buffer(j)%N_det
do k=1,N_int
psi_det(k,1,i+N_det_old) = H_apply_buffer(j)%det(k,1,i)
psi_det(k,2,i+N_det_old) = H_apply_buffer(j)%det(k,2,i)
enddo
ASSERT (sum(popcnt(psi_det(:,1,i+N_det_old))) == elec_alpha_num)
ASSERT (sum(popcnt(psi_det(:,2,i+N_det_old))) == elec_beta_num )
enddo
do k=1,N_states
do i=1,H_apply_buffer(j)%N_det
psi_r_coef_bi_ortho(i+N_det_old,k) = H_apply_buffer(j)%coef(i,k)
psi_l_coef_bi_ortho(i+N_det_old,k) = 0.d0
enddo
enddo
!$OMP BARRIER
H_apply_buffer(j)%N_det = 0
!$OMP END PARALLEL
SOFT_TOUCH N_det psi_det psi_r_coef_bi_ortho psi_l_coef_bi_ortho
logical :: found_duplicates
call remove_duplicates_in_psi_det_tc(found_duplicates)
call bi_normalize(psi_l_coef_bi_ortho,psi_r_coef_bi_ortho,size(psi_l_coef_bi_ortho,1),N_det,N_states)
SOFT_TOUCH N_det psi_det psi_r_coef_bi_ortho psi_l_coef_bi_ortho
end
subroutine remove_duplicates_in_psi_det_tc(found_duplicates)
implicit none
logical, intent(out) :: found_duplicates
BEGIN_DOC
! Removes duplicate determinants in the wave function.
END_DOC
integer :: i,j,k
integer(bit_kind), allocatable :: bit_tmp(:)
logical,allocatable :: duplicate(:)
logical :: dup
allocate (duplicate(N_det), bit_tmp(N_det))
found_duplicates = .False.
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,k,dup)
!$OMP DO
do i=1,N_det
integer, external :: det_search_key
!$DIR FORCEINLINE
bit_tmp(i) = det_search_key(psi_det_sorted_bit_tc(1,1,i),N_int)
duplicate(i) = .False.
enddo
!$OMP END DO
!$OMP DO schedule(dynamic,1024)
do i=1,N_det-1
if (duplicate(i)) then
cycle
endif
j = i+1
do while (bit_tmp(j)==bit_tmp(i))
if (duplicate(j)) then
j = j+1
if (j > N_det) then
exit
else
cycle
endif
endif
dup = .True.
do k=1,N_int
if ( (psi_det_sorted_bit_tc(k,1,i) /= psi_det_sorted_bit_tc(k,1,j) ) &
.or. (psi_det_sorted_bit_tc(k,2,i) /= psi_det_sorted_bit_tc(k,2,j) ) ) then
dup = .False.
exit
endif
enddo
if (dup) then
duplicate(j) = .True.
found_duplicates = .True.
endif
j += 1
if (j > N_det) then
exit
endif
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
if (found_duplicates) then
k=0
do i=1,N_det
if (.not.duplicate(i)) then
k += 1
psi_det(:,:,k) = psi_det_sorted_bit_tc (:,:,i)
psi_r_coef_bi_ortho(k,:) = psi_r_coef_sorted_bit(i,:)
psi_l_coef_bi_ortho(k,:) = psi_l_coef_sorted_bit(i,:)
else
if (sum(abs(psi_r_coef_sorted_bit(i,:))) /= 0.d0 ) then
psi_r_coef_bi_ortho(k,:) = psi_r_coef_sorted_bit(i,:)
psi_l_coef_bi_ortho(k,:) = psi_l_coef_sorted_bit(i,:)
endif
endif
enddo
N_det = k
psi_det_sorted_bit_tc(:,:,1:N_det) = psi_det(:,:,1:N_det)
psi_r_coef_sorted_bit(1:N_det,:) = psi_r_coef_bi_ortho(1:N_det,:)
psi_l_coef_sorted_bit(1:N_det,:) = psi_l_coef_bi_ortho(1:N_det,:)
TOUCH N_det psi_det psi_det_sorted_bit_tc c0_weight psi_r_coef_sorted_bit psi_l_coef_sorted_bit
endif
psi_det = psi_det_sorted_tc
psi_r_coef_bi_ortho = psi_r_coef_sorted_bi_ortho
psi_l_coef_bi_ortho = psi_l_coef_sorted_bi_ortho
SOFT_TOUCH psi_det psi_r_coef_bi_ortho psi_l_coef_bi_ortho psi_det_sorted_bit_tc psi_r_coef_sorted_bit psi_l_coef_sorted_bit
deallocate (duplicate,bit_tmp)
end
BEGIN_PROVIDER [ integer(bit_kind), psi_det_sorted_bit_tc, (N_int,2,psi_det_size) ]
&BEGIN_PROVIDER [ double precision, psi_r_coef_sorted_bit, (N_det,N_states) ]
&BEGIN_PROVIDER [ double precision, psi_l_coef_sorted_bit, (N_det,N_states) ]
implicit none
BEGIN_DOC
! Determinants on which we apply $\langle i|H|psi \rangle$ for perturbation.
! They are sorted by determinants interpreted as integers. Useful
! to accelerate the search of a random determinant in the wave
! function.
END_DOC
call sort_dets_by_det_search_key(N_det, psi_det, psi_r_coef_bi_ortho, size(psi_r_coef_bi_ortho,1), &
psi_det_sorted_bit_tc, psi_r_coef_sorted_bit, N_states)
call sort_dets_by_det_search_key(N_det, psi_det, psi_l_coef_bi_ortho, size(psi_l_coef_bi_ortho,1), &
psi_det_sorted_bit_tc, psi_l_coef_sorted_bit, N_states)
END_PROVIDER

View File

@ -0,0 +1,100 @@
subroutine diagonalize_CI_tc_bi_ortho(ndet, E_tc,norm,pt2_data,print_pt2)
use selection_types
implicit none
integer, intent(inout) :: ndet ! number of determinants from before
double precision, intent(inout) :: E_tc,norm ! E and norm from previous wave function
type(pt2_type) , intent(in) :: pt2_data ! PT2 from previous wave function
logical, intent(in) :: print_pt2
BEGIN_DOC
! Replace the coefficients of the CI states by the coefficients of the
! eigenstates of the CI matrix
END_DOC
integer :: i,j
double precision :: pt2_tmp,pt1_norm,rpt2_tmp,abs_pt2
pt2_tmp = pt2_data % pt2(1)
abs_pt2 = pt2_data % variance(1)
pt1_norm = pt2_data % overlap(1,1)
rpt2_tmp = pt2_tmp/(1.d0 + pt1_norm)
print*,'*****'
print*,'New wave function information'
print*,'N_det tc = ',N_det
print*,'norm_ground_left_right_bi_orth = ',norm_ground_left_right_bi_orth
print*,'eigval_right_tc = ',eigval_right_tc_bi_orth(1)
print*,'Ndet, E_tc = ',N_det,eigval_right_tc_bi_orth(1)
print*,'*****'
if(print_pt2)then
print*,'*****'
print*,'previous wave function info'
print*,'norm(before) = ',norm
print*,'E(before) = ',E_tc
print*,'PT1 norm = ',dsqrt(pt1_norm)
print*,'PT2 = ',pt2_tmp
print*,'rPT2 = ',rpt2_tmp
print*,'|PT2| = ',abs_pt2
print*,'Positive PT2 = ',(pt2_tmp + abs_pt2)*0.5d0
print*,'Negative PT2 = ',(pt2_tmp - abs_pt2)*0.5d0
print*,'E(before) + PT2 = ',E_tc + pt2_tmp/norm
print*,'E(before) +rPT2 = ',E_tc + rpt2_tmp/norm
write(*,'(A28,X,I10,X,100(F16.8,X))')'Ndet,E,E+PT2,E+RPT2,|PT2|=',ndet,E_tc ,E_tc + pt2_tmp/norm,E_tc + rpt2_tmp/norm,abs_pt2
print*,'*****'
endif
E_tc = eigval_right_tc_bi_orth(1)
norm = norm_ground_left_right_bi_orth
ndet = N_det
do j=1,N_states
do i=1,N_det
psi_l_coef_bi_ortho(i,j) = leigvec_tc_bi_orth(i,j)
psi_r_coef_bi_ortho(i,j) = reigvec_tc_bi_orth(i,j)
psi_coef(i,j) = dabs(psi_l_coef_bi_ortho(i,j) * psi_r_coef_bi_ortho(i,j))
enddo
enddo
SOFT_TOUCH eigval_left_tc_bi_orth eigval_right_tc_bi_orth leigvec_tc_bi_orth reigvec_tc_bi_orth norm_ground_left_right_bi_orth psi_coef psi_l_coef_bi_ortho psi_r_coef_bi_ortho
call save_tc_bi_ortho_wavefunction
end
subroutine print_CI_dressed(ndet, E_tc,norm,pt2_data,print_pt2)
use selection_types
implicit none
integer, intent(inout) :: ndet ! number of determinants from before
double precision, intent(inout) :: E_tc,norm ! E and norm from previous wave function
type(pt2_type) , intent(in) :: pt2_data ! PT2 from previous wave function
logical, intent(in) :: print_pt2
BEGIN_DOC
! Replace the coefficients of the CI states by the coefficients of the
! eigenstates of the CI matrix
END_DOC
integer :: i,j
print*,'*****'
print*,'New wave function information'
print*,'N_det tc = ',N_det
print*,'norm_ground_left_right_bi_orth = ',norm_ground_left_right_bi_orth
print*,'eigval_right_tc = ',eigval_right_tc_bi_orth(1)
print*,'Ndet, E_tc = ',N_det,eigval_right_tc_bi_orth(1)
print*,'*****'
if(print_pt2)then
print*,'*****'
print*,'previous wave function info'
print*,'norm(before) = ',norm
print*,'E(before) = ',E_tc
print*,'PT1 norm = ',dsqrt(pt2_data % overlap(1,1))
print*,'E(before) + PT2 = ',E_tc + (pt2_data % pt2(1))/norm
print*,'PT2 = ',pt2_data % pt2(1)
print*,'Ndet, E_tc, E+PT2 = ',ndet,E_tc,E_tc + (pt2_data % pt2(1))/norm,dsqrt(pt2_data % overlap(1,1))
print*,'*****'
endif
E_tc = eigval_right_tc_bi_orth(1)
norm = norm_ground_left_right_bi_orth
ndet = N_det
do j=1,N_states
do i=1,N_det
psi_coef(i,j) = reigvec_tc_bi_orth(i,j)
enddo
enddo
SOFT_TOUCH eigval_left_tc_bi_orth eigval_right_tc_bi_orth leigvec_tc_bi_orth norm_ground_left_right_bi_orth psi_coef reigvec_tc_bi_orth
end

View File

@ -0,0 +1,85 @@
program fci
implicit none
BEGIN_DOC
! Selected Full Configuration Interaction with stochastic selection
! and PT2.
!
! This program performs a |CIPSI|-like selected |CI| using a
! stochastic scheme for both the selection of the important Slater
! determinants and the computation of the |PT2| correction. This
! |CIPSI|-like algorithm will be performed for the lowest states of
! the variational space (see :option:`determinants n_states`). The
! |FCI| program will stop when reaching at least one the two following
! conditions:
!
! * number of Slater determinants > :option:`determinants n_det_max`
! * abs(|PT2|) less than :option:`perturbation pt2_max`
!
! The following other options can be of interest:
!
! :option:`determinants read_wf`
! When set to |false|, the program starts with a ROHF-like Slater
! determinant as a guess wave function. When set to |true|, the
! program starts with the wave function(s) stored in the |EZFIO|
! directory as guess wave function(s).
!
! :option:`determinants s2_eig`
! When set to |true|, the selection will systematically add all the
! necessary Slater determinants in order to have a pure spin wave
! function with an |S^2| value corresponding to
! :option:`determinants expected_s2`.
!
! For excited states calculations, it is recommended to start with
! :ref:`cis` or :ref:`cisd` guess wave functions, eventually in
! a restricted set of |MOs|, and to set :option:`determinants s2_eig`
! to |true|.
!
END_DOC
my_grid_becke = .True.
my_n_pt_r_grid = 30
my_n_pt_a_grid = 50
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
pruning = -1.d0
touch pruning
! pt2_relative_error = 0.01d0
! touch pt2_relative_error
call run_cipsi_tc
end
subroutine run_cipsi_tc
implicit none
if (.not.is_zmq_slave) then
PROVIDE psi_det psi_coef mo_bi_ortho_tc_two_e mo_bi_ortho_tc_one_e
if(elec_alpha_num+elec_beta_num.ge.3)then
if(three_body_h_tc)then
call provide_all_three_ints_bi_ortho
endif
endif
! ---
if (do_pt2) then
call run_stochastic_cipsi
else
call run_cipsi
endif
else
PROVIDE mo_bi_ortho_tc_one_e mo_bi_ortho_tc_two_e pt2_min_parallel_tasks
if(elec_alpha_num+elec_beta_num.ge.3)then
if(three_body_h_tc)then
call provide_all_three_ints_bi_ortho
endif
endif
! ---
call run_slave_cipsi
endif
end

View File

@ -0,0 +1,51 @@
use bitmasks
BEGIN_PROVIDER [ integer, N_det_generators ]
implicit none
BEGIN_DOC
! For Single reference wave functions, the number of generators is 1 : the
! Hartree-Fock determinant
END_DOC
integer :: i
double precision :: norm
call write_time(6)
norm = 1.d0
N_det_generators = N_det
do i=1,N_det
norm = norm - psi_average_norm_contrib_sorted_tc(i)
if (norm - 1.d-10 < 1.d0 - threshold_generators) then
N_det_generators = i
exit
endif
enddo
N_det_generators = max(N_det_generators,1)
call write_int(6,N_det_generators,'Number of generators')
END_PROVIDER
BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators, (N_int,2,psi_det_size) ]
&BEGIN_PROVIDER [ double precision, psi_coef_generators, (psi_det_size,N_states) ]
implicit none
BEGIN_DOC
! For Single reference wave functions, the generator is the
! Hartree-Fock determinant
END_DOC
psi_det_generators(1:N_int,1:2,1:N_det) = psi_det_sorted_tc(1:N_int,1:2,1:N_det)
psi_coef_generators(1:N_det,1:N_states) = psi_coef_sorted_tc(1:N_det,1:N_states)
END_PROVIDER
BEGIN_PROVIDER [ integer(bit_kind), psi_det_sorted_tc_gen, (N_int,2,psi_det_size) ]
&BEGIN_PROVIDER [ double precision, psi_coef_sorted_tc_gen, (psi_det_size,N_states) ]
&BEGIN_PROVIDER [ integer, psi_det_sorted_tc_gen_order, (psi_det_size) ]
implicit none
BEGIN_DOC
! For Single reference wave functions, the generator is the
! Hartree-Fock determinant
END_DOC
psi_det_sorted_tc_gen = psi_det_sorted_tc
psi_coef_sorted_tc_gen = psi_coef_sorted_tc
psi_det_sorted_tc_gen_order = psi_det_sorted_tc_order
END_PROVIDER

View File

@ -0,0 +1,9 @@
subroutine save_energy(E,pt2)
implicit none
BEGIN_DOC
! Saves the energy in |EZFIO|.
END_DOC
double precision, intent(in) :: E(N_states), pt2(N_states)
call ezfio_set_fci_tc_energy(E(1:N_states))
call ezfio_set_fci_tc_energy_pt2(E(1:N_states)+pt2(1:N_states))
end

View File

@ -0,0 +1,6 @@
3
C 6.000000 0.000000 0.000000 0.173480
H 1.000000 0.000000 -0.861500 -0.520430
H 1.000000 0.000000 0.861500 -0.520430

View File

@ -0,0 +1,5 @@
2
H 0.000000 0.000000 -0.825120
F 0.000000 0.000000 0.091680

View File

@ -0,0 +1,16 @@
input=h2o
basis=dz
EZFIO=${input}_${basis}_bi_ortho
file=${EZFIO}.tc_fci.out
grep "Ndet,E,E+PT2,E+RPT2,|PT2|=" ${file} | cut -d "=" -f 2 > data_${EZFIO}
file=${EZFIO}.tc_fci_normal_order.out
grep "Ndet,E,E+PT2,E+RPT2=" ${file} | cut -d "=" -f 2 > data_${EZFIO}_normal
#EZFIO=${input}_${basis}_ortho
#file=${EZFIO}.tc_fci.out
#grep "Ndet, E_tc, E+PT2 =" ${file} | cut -d "=" -f 2 > data_${EZFIO}
#file=${EZFIO}.tc_fci_normal_order.out
#grep "Ndet, E_tc, E+PT2 =" ${file} | cut -d "=" -f 2 > data_${EZFIO}_normal
#zip data_${input}_${basis}.zip data*

View File

@ -0,0 +1,41 @@
#!/bin/bash
# This is a sample PBS script
# temps CPU a ajuster au calcul
#PBS -l cput=2000:00:00
#PBS -l nodes=1:ppn=16
# memoire a ajuster au calcul
#PBS -l vmem=100gb
# a changer
# Pour savoir sur quel noeud on est
#echo $HOSTNAME
# Startdir = ou sont les fichiers d'input, par defaut HOMEDIR
#
StartDir=$PBS_O_WORKDIR
echo $StartDir
#
# SCRATCHDIR = espace temporaire (local au noeud et a vider apres le calcul)
# NE PAS MODIFIER
ulimit -s unlimited
export SCRATCHDIR=/scratch/$USER/$PBS_JOBID
#
cd $StartDir
############################################################################
#### EXAMPLE OF SCRIPT TO RUN A CIPSI CALCULATION ON 5 STATES ON THE Ne^+ CATION
#### USING NATURAL ORBITALS OF A SMALL CIPSI AS MOS
#### ALL STATES WILL HAVE THE SAME SPIN SIMETRY : A DOUBLET
####### YOU PUT THE PATH TO YOUR
QP_ROOT=/home_lct/eginer/programs/qp2
source ${QP_ROOT}/quantum_package.rc
####### YOU LOAD SOME LIBRARIES
alias python3='/programmes/installation/Python/3.7.1/bin/python3'
type -a python3
export OMP_NUM_THREADS=16
module load intel2016_OMPI-V2
source ~/programs/qp2/quantum_package.rc
./script.sh h2o dz O 1

View File

@ -0,0 +1,6 @@
3
O 0.000000 0.000000 0.000000
H 0.000000 0.000000 0.957200
H -0.926627 0.000000 -0.239987

View File

@ -0,0 +1,31 @@
source /home_lct/eginer/qp2/quantum_package.rc
input=$1
basis=$2
atom=$3
mul=$4
EXPORT_OMP_NUM_THREADS=16
dir=${input}_${basis}
mkdir ${dir}
cp ${input}.xyz ${dir}/
cd $dir
EZFIO=${input}_${basis}_bi_ortho
qp create_ezfio -b "${atom}:cc-pcv${basis}|H:cc-pv${basis}" ${input}.xyz -m $mul -o $EZFIO
qp run scf
# Getting THE GOOD VALUE OF MU
qp run print_mu_av_tc | tee ${EZFIO_FILE}.mu_av.out
mu=`grep "average_mu_rs_c_lda =" ${EZFIO_FILE}.mu_av.out | cut -d "=" -f 2`
qp set ao_two_e_erf_ints mu_erf $mu
# Carrying the BI-ORTHO TC-SCF
qp run tc_scf | tee ${EZFIO_FILE}.tc_scf.out
# Three body terms without normal order
### THREE E TERMS FOR FCI
qp set tc_keywords three_body_h_tc True
qp set tc_keywords double_normal_ord False
qp set perturbation pt2_max 0.003
qp run fci_tc_bi_ortho | tee ${EZFIO_FILE}.tc_fci.out
# Three body terms with normal order
qp set tc_keywords double_normal_ord True
qp run fci_tc_bi_ortho | tee ${EZFIO_FILE}.tc_fci_normal_order.out
cd ../

View File

@ -0,0 +1,92 @@
use bitmasks
BEGIN_PROVIDER [ double precision, threshold_selectors ]
implicit none
BEGIN_DOC
! Thresholds on selectors (fraction of the square of the norm)
END_DOC
threshold_selectors = dsqrt(threshold_generators)
END_PROVIDER
BEGIN_PROVIDER [ integer, N_det_selectors]
implicit none
BEGIN_DOC
! For Single reference wave functions, the number of selectors is 1 : the
! Hartree-Fock determinant
END_DOC
integer :: i
double precision :: norm, norm_max
call write_time(6)
N_det_selectors = N_det
norm = 1.d0
do i=1,N_det
norm = norm - psi_average_norm_contrib_tc(i)
if (norm - 1.d-10 < 1.d0 - threshold_selectors) then
N_det_selectors = i
exit
endif
enddo
N_det_selectors = max(N_det_selectors,N_det_generators)
call write_int(6,N_det_selectors,'Number of selectors')
END_PROVIDER
BEGIN_PROVIDER [ integer(bit_kind), psi_selectors, (N_int,2,psi_selectors_size) ]
&BEGIN_PROVIDER [ double precision, psi_selectors_coef, (psi_selectors_size,N_states) ]
implicit none
BEGIN_DOC
! Determinants on which we apply <i|H|psi> for perturbation.
END_DOC
integer :: i,k
do i=1,N_det_selectors
do k=1,N_int
psi_selectors(k,1,i) = psi_det_sorted_tc(k,1,i)
psi_selectors(k,2,i) = psi_det_sorted_tc(k,2,i)
enddo
enddo
do k=1,N_states
do i=1,N_det_selectors
psi_selectors_coef(i,k) = psi_coef_sorted_tc_gen(i,k)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, psi_selectors_coef_transp, (N_states,psi_selectors_size) ]
implicit none
BEGIN_DOC
! Transposed psi_selectors
END_DOC
integer :: i,k
do i=1,N_det_selectors
do k=1,N_states
psi_selectors_coef_transp(k,i) = psi_selectors_coef(i,k)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, psi_selectors_rcoef_bi_orth_transp, (N_states, psi_det_size) ]
&BEGIN_PROVIDER [ double precision, psi_selectors_lcoef_bi_orth_transp, (N_states, psi_det_size) ]
implicit none
integer :: i, k
psi_selectors_rcoef_bi_orth_transp = 0.d0
psi_selectors_lcoef_bi_orth_transp = 0.d0
print*,'N_det,N_det_selectors',N_det,N_det_selectors
do i = 1, N_det_selectors
do k = 1, N_states
psi_selectors_rcoef_bi_orth_transp(k,i) = psi_r_coef_sorted_bi_ortho(i,k)
psi_selectors_lcoef_bi_orth_transp(k,i) = psi_l_coef_sorted_bi_ortho(i,k)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ integer, psi_selectors_size ]
implicit none
psi_selectors_size = psi_det_size
END_PROVIDER

103
src/fci_tc_bi/zmq.irp.f Normal file
View File

@ -0,0 +1,103 @@
BEGIN_TEMPLATE
integer function zmq_put_$X(zmq_to_qp_run_socket,worker_id)
use f77_zmq
implicit none
BEGIN_DOC
! Put $X on the qp_run scheduler
END_DOC
integer(ZMQ_PTR), intent(in) :: zmq_to_qp_run_socket
integer, intent(in) :: worker_id
integer :: rc
character*(256) :: msg
zmq_put_$X = 0
write(msg,'(A,1X,I8,1X,A200)') 'put_data '//trim(zmq_state), worker_id, '$X'
rc = f77_zmq_send(zmq_to_qp_run_socket,trim(msg),len(trim(msg)),ZMQ_SNDMORE)
if (rc /= len(trim(msg))) then
zmq_put_$X = -1
return
endif
rc = f77_zmq_send(zmq_to_qp_run_socket,$X,4,0)
if (rc /= 4) then
zmq_put_$X = -1
return
endif
rc = f77_zmq_recv(zmq_to_qp_run_socket,msg,len(msg),0)
if (msg(1:rc) /= 'put_data_reply ok') then
zmq_put_$X = -1
return
endif
end
integer function zmq_get_$X(zmq_to_qp_run_socket, worker_id)
use f77_zmq
implicit none
BEGIN_DOC
! Get $X from the qp_run scheduler
END_DOC
integer(ZMQ_PTR), intent(in) :: zmq_to_qp_run_socket
integer, intent(in) :: worker_id
integer :: rc
character*(256) :: msg
PROVIDE zmq_state
zmq_get_$X = 0
if (mpi_master) then
write(msg,'(A,1X,I8,1X,A200)') 'get_data '//trim(zmq_state), worker_id, '$X'
rc = f77_zmq_send(zmq_to_qp_run_socket,trim(msg),len(trim(msg)),0)
if (rc /= len(trim(msg))) then
zmq_get_$X = -1
go to 10
endif
rc = f77_zmq_recv(zmq_to_qp_run_socket,msg,len(msg),0)
if (msg(1:14) /= 'get_data_reply') then
zmq_get_$X = -1
go to 10
endif
rc = f77_zmq_recv(zmq_to_qp_run_socket,$X,4,0)
if (rc /= 4) then
zmq_get_$X = -1
go to 10
endif
endif
10 continue
IRP_IF MPI_DEBUG
print *, irp_here, mpi_rank
call MPI_BARRIER(MPI_COMM_WORLD, ierr)
IRP_ENDIF
IRP_IF MPI
include 'mpif.h'
integer :: ierr
call MPI_BCAST (zmq_get_$X, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
if (ierr /= MPI_SUCCESS) then
print *, irp_here//': Unable to broadcast N_det_generators'
stop -1
endif
call MPI_BCAST ($X, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
if (ierr /= MPI_SUCCESS) then
print *, irp_here//': Unable to broadcast N_det_generators'
stop -1
endif
IRP_ENDIF
end
SUBST [ X ]
N_det_generators ;;
N_det_selectors ;;
END_TEMPLATE