From 70516c8c05f625085a72dc9a95031c8d06fd2aec Mon Sep 17 00:00:00 2001 From: eginer Date: Wed, 5 Oct 2022 17:59:23 +0200 Subject: [PATCH] added fci_tc_bi --- src/bi_ort_ints/bi_ort_ints.irp.f | 104 +-------- src/fci_tc_bi/EZFIO.cfg | 17 ++ src/fci_tc_bi/NEED | 3 + src/fci_tc_bi/class.irp.f | 12 + src/fci_tc_bi/copy_wf.irp.f | 215 ++++++++++++++++++ src/fci_tc_bi/diagonalize_ci.irp.f | 100 ++++++++ src/fci_tc_bi/fci_tc_bi_ortho.irp.f | 85 +++++++ src/fci_tc_bi/generators.irp.f | 51 +++++ src/fci_tc_bi/save_energy.irp.f | 9 + src/fci_tc_bi/scripts_fci_tc/CH2.xyz | 6 + src/fci_tc_bi/scripts_fci_tc/FH.xyz | 5 + .../scripts_fci_tc/extract_tables.sh | 16 ++ src/fci_tc_bi/scripts_fci_tc/h2o.sh | 41 ++++ src/fci_tc_bi/scripts_fci_tc/h2o.xyz | 6 + src/fci_tc_bi/scripts_fci_tc/script.sh | 31 +++ src/fci_tc_bi/selectors.irp.f | 92 ++++++++ src/fci_tc_bi/zmq.irp.f | 103 +++++++++ 17 files changed, 793 insertions(+), 103 deletions(-) create mode 100644 src/fci_tc_bi/EZFIO.cfg create mode 100644 src/fci_tc_bi/NEED create mode 100644 src/fci_tc_bi/class.irp.f create mode 100644 src/fci_tc_bi/copy_wf.irp.f create mode 100644 src/fci_tc_bi/diagonalize_ci.irp.f create mode 100644 src/fci_tc_bi/fci_tc_bi_ortho.irp.f create mode 100644 src/fci_tc_bi/generators.irp.f create mode 100644 src/fci_tc_bi/save_energy.irp.f create mode 100644 src/fci_tc_bi/scripts_fci_tc/CH2.xyz create mode 100644 src/fci_tc_bi/scripts_fci_tc/FH.xyz create mode 100755 src/fci_tc_bi/scripts_fci_tc/extract_tables.sh create mode 100644 src/fci_tc_bi/scripts_fci_tc/h2o.sh create mode 100644 src/fci_tc_bi/scripts_fci_tc/h2o.xyz create mode 100755 src/fci_tc_bi/scripts_fci_tc/script.sh create mode 100644 src/fci_tc_bi/selectors.irp.f create mode 100644 src/fci_tc_bi/zmq.irp.f diff --git a/src/bi_ort_ints/bi_ort_ints.irp.f b/src/bi_ort_ints/bi_ort_ints.irp.f index 6884ff38..bb894b44 100644 --- a/src/bi_ort_ints/bi_ort_ints.irp.f +++ b/src/bi_ort_ints/bi_ort_ints.irp.f @@ -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 - ! - ! - ! 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) = -! 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 - - - - diff --git a/src/fci_tc_bi/EZFIO.cfg b/src/fci_tc_bi/EZFIO.cfg new file mode 100644 index 00000000..a2552c74 --- /dev/null +++ b/src/fci_tc_bi/EZFIO.cfg @@ -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 diff --git a/src/fci_tc_bi/NEED b/src/fci_tc_bi/NEED new file mode 100644 index 00000000..000b0deb --- /dev/null +++ b/src/fci_tc_bi/NEED @@ -0,0 +1,3 @@ +tc_bi_ortho +davidson_undressed +cipsi_tc_bi_ortho diff --git a/src/fci_tc_bi/class.irp.f b/src/fci_tc_bi/class.irp.f new file mode 100644 index 00000000..b4a68ac2 --- /dev/null +++ b/src/fci_tc_bi/class.irp.f @@ -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 + diff --git a/src/fci_tc_bi/copy_wf.irp.f b/src/fci_tc_bi/copy_wf.irp.f new file mode 100644 index 00000000..cdb1ead8 --- /dev/null +++ b/src/fci_tc_bi/copy_wf.irp.f @@ -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 diff --git a/src/fci_tc_bi/diagonalize_ci.irp.f b/src/fci_tc_bi/diagonalize_ci.irp.f new file mode 100644 index 00000000..56c561ac --- /dev/null +++ b/src/fci_tc_bi/diagonalize_ci.irp.f @@ -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 + diff --git a/src/fci_tc_bi/fci_tc_bi_ortho.irp.f b/src/fci_tc_bi/fci_tc_bi_ortho.irp.f new file mode 100644 index 00000000..84ac8166 --- /dev/null +++ b/src/fci_tc_bi/fci_tc_bi_ortho.irp.f @@ -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 diff --git a/src/fci_tc_bi/generators.irp.f b/src/fci_tc_bi/generators.irp.f new file mode 100644 index 00000000..55c0cbb9 --- /dev/null +++ b/src/fci_tc_bi/generators.irp.f @@ -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 + + diff --git a/src/fci_tc_bi/save_energy.irp.f b/src/fci_tc_bi/save_energy.irp.f new file mode 100644 index 00000000..7c41d00f --- /dev/null +++ b/src/fci_tc_bi/save_energy.irp.f @@ -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 diff --git a/src/fci_tc_bi/scripts_fci_tc/CH2.xyz b/src/fci_tc_bi/scripts_fci_tc/CH2.xyz new file mode 100644 index 00000000..9fa57f4b --- /dev/null +++ b/src/fci_tc_bi/scripts_fci_tc/CH2.xyz @@ -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 + diff --git a/src/fci_tc_bi/scripts_fci_tc/FH.xyz b/src/fci_tc_bi/scripts_fci_tc/FH.xyz new file mode 100644 index 00000000..9a1563f4 --- /dev/null +++ b/src/fci_tc_bi/scripts_fci_tc/FH.xyz @@ -0,0 +1,5 @@ +2 + +H 0.000000 0.000000 -0.825120 +F 0.000000 0.000000 0.091680 + diff --git a/src/fci_tc_bi/scripts_fci_tc/extract_tables.sh b/src/fci_tc_bi/scripts_fci_tc/extract_tables.sh new file mode 100755 index 00000000..a585884e --- /dev/null +++ b/src/fci_tc_bi/scripts_fci_tc/extract_tables.sh @@ -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* diff --git a/src/fci_tc_bi/scripts_fci_tc/h2o.sh b/src/fci_tc_bi/scripts_fci_tc/h2o.sh new file mode 100644 index 00000000..d0afca30 --- /dev/null +++ b/src/fci_tc_bi/scripts_fci_tc/h2o.sh @@ -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 diff --git a/src/fci_tc_bi/scripts_fci_tc/h2o.xyz b/src/fci_tc_bi/scripts_fci_tc/h2o.xyz new file mode 100644 index 00000000..dee51ffc --- /dev/null +++ b/src/fci_tc_bi/scripts_fci_tc/h2o.xyz @@ -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 + diff --git a/src/fci_tc_bi/scripts_fci_tc/script.sh b/src/fci_tc_bi/scripts_fci_tc/script.sh new file mode 100755 index 00000000..58585658 --- /dev/null +++ b/src/fci_tc_bi/scripts_fci_tc/script.sh @@ -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 ../ + diff --git a/src/fci_tc_bi/selectors.irp.f b/src/fci_tc_bi/selectors.irp.f new file mode 100644 index 00000000..734c8ed0 --- /dev/null +++ b/src/fci_tc_bi/selectors.irp.f @@ -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 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 + diff --git a/src/fci_tc_bi/zmq.irp.f b/src/fci_tc_bi/zmq.irp.f new file mode 100644 index 00000000..cb2df483 --- /dev/null +++ b/src/fci_tc_bi/zmq.irp.f @@ -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 +