From 14c43916ca549b3469b0424d86dbd37eddf83dcc Mon Sep 17 00:00:00 2001 From: Manu Date: Wed, 21 May 2014 22:41:26 +0200 Subject: [PATCH 1/4] Changes in /Hartree_Fock/README.srt --- src/AOs/README.rst | 24 +++---- src/BiInts/README.rst | 94 +++++++++++++-------------- src/Bitmask/bitmasks_routines.irp.f | 6 +- src/CISD/README.rst | 36 ++++++++++- src/CISD/SC2.irp.f | 27 ++++---- src/CISD/cisd.irp.f | 20 +++++- src/Hartree_Fock/README.rst | 3 - src/Nuclei/README.rst | 24 +++---- src/Utils/README.rst | 99 +++++++++++++++++------------ src/Utils/util.irp.f | 2 +- 10 files changed, 202 insertions(+), 133 deletions(-) diff --git a/src/AOs/README.rst b/src/AOs/README.rst index 378bf8ae..dfbef9ca 100644 --- a/src/AOs/README.rst +++ b/src/AOs/README.rst @@ -50,44 +50,44 @@ Documentation .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. -`ao_coef `_ +`ao_coef `_ Coefficients, exponents and powers of x,y and z ao_coef(i,j) = coefficient of the jth primitive on the ith ao -`ao_coef_transp `_ +`ao_coef_transp `_ Transposed ao_coef and ao_expo -`ao_expo `_ +`ao_expo `_ Coefficients, exponents and powers of x,y and z ao_coef(i,j) = coefficient of the jth primitive on the ith ao -`ao_expo_transp `_ +`ao_expo_transp `_ Transposed ao_coef and ao_expo -`ao_nucl `_ +`ao_nucl `_ Index of the nuclei on which the ao is centered -`ao_num `_ +`ao_num `_ Number of atomic orbitals -`ao_num_align `_ +`ao_num_align `_ Number of atomic orbitals -`ao_overlap `_ +`ao_overlap `_ matrix of the overlap for tha AOs S(i,j) = overlap between the ith and the jth atomic basis function -`ao_power `_ +`ao_power `_ Coefficients, exponents and powers of x,y and z ao_coef(i,j) = coefficient of the jth primitive on the ith ao -`ao_prim_num `_ +`ao_prim_num `_ Number of primitives per atomic orbital -`ao_prim_num_max `_ +`ao_prim_num_max `_ Undocumented -`ao_prim_num_max_align `_ +`ao_prim_num_max_align `_ Undocumented diff --git a/src/BiInts/README.rst b/src/BiInts/README.rst index 48fa022b..76d944b1 100644 --- a/src/BiInts/README.rst +++ b/src/BiInts/README.rst @@ -32,162 +32,162 @@ Documentation .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. -`ao_bielec_integral `_ +`ao_bielec_integral `_ integral of the AO basis or (ij|kl) i(r1) j(r1) 1/r12 k(r2) l(r2) -`ao_bielec_integral_schwartz `_ +`ao_bielec_integral_schwartz `_ Needed to compuet Schwartz inequalities -`ao_bielec_integrals_in_map `_ +`ao_bielec_integrals_in_map `_ Map of Atomic integrals i(r1) j(r2) 1/r12 k(r1) l(r2) -`compute_ao_bielec_integrals `_ +`compute_ao_bielec_integrals `_ Compute AO 1/r12 integrals for all i and fixed j,k,l -`eri `_ +`eri `_ ATOMIC PRIMTIVE bielectronic integral between the 4 primitives :: primitive_1 = x1**(a_x) y1**(a_y) z1**(a_z) exp(-alpha * r1**2) primitive_2 = x1**(b_x) y1**(b_y) z1**(b_z) exp(- beta * r1**2) primitive_3 = x2**(c_x) y2**(c_y) z2**(c_z) exp(-delta * r2**2) primitive_4 = x2**(d_x) y2**(d_y) z2**(d_z) exp(- gama * r2**2) -`general_primitive_integral `_ +`general_primitive_integral `_ Computes the integral where p,q,r,s are Gaussian primitives -`give_polynom_mult_center_x `_ +`give_polynom_mult_center_x `_ subroutine that returns the explicit polynom in term of the "t" variable of the following polynomw : I_x1(a_x, d_x,p,q) * I_x1(a_y, d_y,p,q) * I_x1(a_z, d_z,p,q) -`i_x1_new `_ +`i_x1_new `_ recursive function involved in the bielectronic integral -`i_x1_pol_mult `_ +`i_x1_pol_mult `_ recursive function involved in the bielectronic integral -`i_x1_pol_mult_a1 `_ +`i_x1_pol_mult_a1 `_ recursive function involved in the bielectronic integral -`i_x1_pol_mult_a2 `_ +`i_x1_pol_mult_a2 `_ recursive function involved in the bielectronic integral -`i_x1_pol_mult_recurs `_ +`i_x1_pol_mult_recurs `_ recursive function involved in the bielectronic integral -`i_x2_new `_ +`i_x2_new `_ recursive function involved in the bielectronic integral -`i_x2_pol_mult `_ +`i_x2_pol_mult `_ recursive function involved in the bielectronic integral -`integrale_new `_ +`integrale_new `_ calculate the integral of the polynom :: I_x1(a_x+b_x, c_x+d_x,p,q) * I_x1(a_y+b_y, c_y+d_y,p,q) * I_x1(a_z+b_z, c_z+d_z,p,q) between ( 0 ; 1) -`n_pt_sup `_ +`n_pt_sup `_ Returns the upper boundary of the degree of the polynom involved in the bielctronic integral : Ix(a_x,b_x,c_x,d_x) * Iy(a_y,b_y,c_y,d_y) * Iz(a_z,b_z,c_z,d_z) -`gauleg `_ +`gauleg `_ Gauss-Legendre -`gauleg_t2 `_ +`gauleg_t2 `_ t_w(i,1,k) = w(i) t_w(i,2,k) = t(i) -`gauleg_w `_ +`gauleg_w `_ t_w(i,1,k) = w(i) t_w(i,2,k) = t(i) -`ao_integrals_map `_ +`ao_integrals_map `_ AO integrals -`bielec_integrals_index `_ +`bielec_integrals_index `_ Undocumented -`clear_ao_map `_ +`clear_ao_map `_ Frees the memory of the AO map -`clear_mo_map `_ +`clear_mo_map `_ Frees the memory of the MO map -`get_ao_bielec_integral `_ +`get_ao_bielec_integral `_ Gets one AO bi-electronic integral from the AO map -`get_ao_bielec_integrals `_ +`get_ao_bielec_integrals `_ Gets multiple AO bi-electronic integral from the AO map . All i are retrieved for j,k,l fixed. -`get_ao_bielec_integrals_non_zero `_ +`get_ao_bielec_integrals_non_zero `_ Gets multiple AO bi-electronic integral from the AO map . All non-zero i are retrieved for j,k,l fixed. -`get_ao_map_size `_ +`get_ao_map_size `_ Returns the number of elements in the AO map -`get_mo_bielec_integral `_ +`get_mo_bielec_integral `_ Returns one integral in the MO basis -`get_mo_bielec_integrals `_ +`get_mo_bielec_integrals `_ Returns multiple integrals in the MO basis, all i for j,k,l fixed. -`get_mo_map_size `_ +`get_mo_map_size `_ Return the number of elements in the MO map -`insert_into_ao_integrals_map `_ +`insert_into_ao_integrals_map `_ Create new entry into AO map -`insert_into_mo_integrals_map `_ +`insert_into_mo_integrals_map `_ Create new entry into MO map, or accumulate in an existing entry -`mo_bielec_integral `_ +`mo_bielec_integral `_ Returns one integral in the MO basis -`mo_integrals_map `_ +`mo_integrals_map `_ MO integrals -`add_integrals_to_map `_ +`add_integrals_to_map `_ Adds integrals to tha MO map according to some bitmask -`mo_bielec_integral_jj `_ +`mo_bielec_integral_jj `_ Transform Bi-electronic integrals and -`mo_bielec_integral_jj_anti `_ +`mo_bielec_integral_jj_anti `_ Transform Bi-electronic integrals and -`mo_bielec_integral_jj_exchange `_ +`mo_bielec_integral_jj_exchange `_ Transform Bi-electronic integrals and -`mo_bielec_integrals_in_map `_ +`mo_bielec_integrals_in_map `_ If True, the map of MO bielectronic integrals is provided -`mo_bielec_integrals_index `_ +`mo_bielec_integrals_index `_ Computes an unique index for i,j,k,l integrals -`ao_integrals_threshold `_ +`ao_integrals_threshold `_ If < ao_integrals_threshold, = 0 -`do_direct_integrals `_ +`do_direct_integrals `_ If True, compute integrals on the fly -`mo_integrals_threshold `_ +`mo_integrals_threshold `_ If < mo_integrals_threshold, = 0 -`read_ao_integrals `_ +`read_ao_integrals `_ If true, read AO integrals in EZFIO -`read_mo_integrals `_ +`read_mo_integrals `_ If true, read MO integrals in EZFIO -`write_ao_integrals `_ +`write_ao_integrals `_ If true, write AO integrals in EZFIO -`write_mo_integrals `_ +`write_mo_integrals `_ If true, write MO integrals in EZFIO diff --git a/src/Bitmask/bitmasks_routines.irp.f b/src/Bitmask/bitmasks_routines.irp.f index cb5f98eb..6764c48b 100644 --- a/src/Bitmask/bitmasks_routines.irp.f +++ b/src/Bitmask/bitmasks_routines.irp.f @@ -119,9 +119,9 @@ subroutine debug_det(string,Nint) integer, intent(in) :: Nint integer(bit_kind), intent(in) :: string(Nint,2) character*(512) :: output(2) - call bitstring_to_hexa( output(1), string(1,1), Nint ) - call bitstring_to_hexa( output(2), string(1,2), Nint ) - print *, trim(output(1)) , '|', trim(output(2)) +! call bitstring_to_hexa( output(1), string(1,1), Nint ) +! call bitstring_to_hexa( output(2), string(1,2), Nint ) +! print *, trim(output(1)) , '|', trim(output(2)) call bitstring_to_str( output(1), string(1,1), N_int ) call bitstring_to_str( output(2), string(1,2), N_int ) diff --git a/src/CISD/README.rst b/src/CISD/README.rst index fbf1807f..59e077cd 100644 --- a/src/CISD/README.rst +++ b/src/CISD/README.rst @@ -30,7 +30,41 @@ Documentation Calls H_apply on the HF determinant and selects all connected single and double excitations (of the same symmetry). -`cisd `_ +`cisd_sc2 `_ + CISD+SC2 method :: take off all the disconnected terms of a CISD (selected or not) + .br + dets_in : bitmasks corresponding to determinants + .br + u_in : guess coefficients on the various states. Overwritten + on exit + .br + dim_in : leftmost dimension of u_in + .br + sze : Number of determinants + .br + N_st : Number of eigenstates + .br + Initial guess vectors are not necessarily orthonormal + +`davidson_diag_hjj `_ + Davidson diagonalization with specific diagonal elements of the H matrix + .br + H_jj : specific diagonal H matrix elements to diagonalize de Davidson + .br + dets_in : bitmasks corresponding to determinants + .br + u_in : guess coefficients on the various states. Overwritten + on exit + .br + dim_in : leftmost dimension of u_in + .br + sze : Number of determinants + .br + N_st : Number of eigenstates + .br + Initial guess vectors are not necessarily orthonormal + +`repeat_excitation `_ Undocumented diff --git a/src/CISD/SC2.irp.f b/src/CISD/SC2.irp.f index 330e3665..5a2c9dfd 100644 --- a/src/CISD/SC2.irp.f +++ b/src/CISD/SC2.irp.f @@ -36,6 +36,7 @@ subroutine CISD_SC2(dets_in,u_in,energies,dim_in,sze,N_st,Nint) double precision :: hij_elec, e_corr_double,e_corr,diag_h_mat_elem,inv_c0 double precision :: e_corr_array(sze),H_jj_ref(sze),H_jj_dressed(sze),hij_double(sze) double precision :: e_corr_double_before,accu,cpu_2,cpu_1 + integer :: degree_exc(sze) integer :: i_ok !$OMP PARALLEL DEFAULT(NONE) & @@ -55,6 +56,7 @@ subroutine CISD_SC2(dets_in,u_in,energies,dim_in,sze,N_st,Nint) e_corr_double = 0.d0 do i = 1, sze call get_excitation_degree(ref_bitmask,dets_in(1,1,i),degree,Nint) + degree_exc(i) = degree if(degree==0)then index_hf=i else if (degree == 2)then @@ -78,9 +80,6 @@ subroutine CISD_SC2(dets_in,u_in,energies,dim_in,sze,N_st,Nint) enddo e_corr = e_corr * inv_c0 e_corr_double = e_corr_double * inv_c0 - print*, 'E_corr = ',e_corr - print*, 'E_corr_double = ', e_corr_double - converged = .False. e_corr_double_before = e_corr_double iter = 0 @@ -93,19 +92,24 @@ subroutine CISD_SC2(dets_in,u_in,energies,dim_in,sze,N_st,Nint) H_jj_dressed(i) = H_jj_ref(i) if (i==index_hf)cycle accu = 0.d0 - do j=1,N_double - call repeat_excitation(dets_in(1,1,i),ref_bitmask,dets_in(1,1,index_double(j)),i_ok,Nint) - if (i_ok==1)cycle! you check if the excitation is possible - accu += e_corr_array(j) - enddo + if(degree_exc(i)==1)then + do j=1,N_double + call get_excitation_degree(dets_in(1,1,i),dets_in(1,1,index_double(j)),degree,N_int) + if (degree<=2)cycle + accu += e_corr_array(j) + enddo + else + do j=1,N_double + call get_excitation_degree(dets_in(1,1,i),dets_in(1,1,index_double(j)),degree,N_int) + if (degree<=3)cycle + accu += e_corr_array(j) + enddo + endif H_jj_dressed(i) += accu enddo call cpu_time(cpu_2) print*,'time for the excitations = ',cpu_2 - cpu_1 - print*,H_jj_ref(1),H_jj_ref(2) - print*,H_jj_dressed(1),H_jj_dressed(2) - print*,u_in(index_hf,1),u_in(index_double(1),1) call davidson_diag_hjj(dets_in,u_in,H_jj_dressed,energies,dim_in,sze,N_st,Nint) print*,u_in(index_hf,1),u_in(index_double(1),1) e_corr_double = 0.d0 @@ -418,4 +422,3 @@ subroutine repeat_excitation(key_in,key_1,key_2,i_ok,Nint) return endif end - diff --git a/src/CISD/cisd.irp.f b/src/CISD/cisd.irp.f index f3bb5488..29f5590a 100644 --- a/src/CISD/cisd.irp.f +++ b/src/CISD/cisd.irp.f @@ -5,18 +5,31 @@ program cisd PROVIDE ref_bitmask_energy double precision :: pt2(10), norm_pert(10), H_pert_diag + double precision,allocatable :: H_jj(:) + double precision :: diag_h_mat_elem + integer,allocatable :: iorder(:) - N_states = 3 + N_states = 6 touch N_states call H_apply_cisd allocate(eigvalues(n_states),eigvectors(n_det,n_states)) print *, 'N_det = ', N_det print *, 'N_states = ', N_states psi_coef = - 1.d-4 + allocate(H_jj(n_det),iorder(n_det)) + do i = 1, N_det + H_jj(i) = diag_h_mat_elem(psi_det(1,1,i),N_int) + iorder(i) = i + enddo + call dsort(H_jj,iorder,n_det) + do k=1,N_states - psi_coef(k,k) = 1.d0 + psi_coef(iorder(k),k) = 1.d0 enddo call davidson_diag(psi_det,psi_coef,eigvalues,size(psi_coef,1),N_det,N_states,N_int) + do i = 1, N_states + print*,'eigvalues(i) = ',eigvalues(i) + enddo print *, '---' print *, 'HF:', HF_energy @@ -26,5 +39,8 @@ program cisd print *, 'E_corr = ',eigvalues(i) - ref_bitmask_energy enddo call CISD_SC2(psi_det,psi_coef,eigvalues,size(psi_coef,1),N_det,N_states,N_int) + do i = 1, N_states + print*,'eigvalues(i) = ',eigvalues(i) + enddo deallocate(eigvalues,eigvectors) end diff --git a/src/Hartree_Fock/README.rst b/src/Hartree_Fock/README.rst index f6b194b8..c2c0840f 100644 --- a/src/Hartree_Fock/README.rst +++ b/src/Hartree_Fock/README.rst @@ -81,9 +81,6 @@ Documentation `eigenvectors_fock_matrix_mo `_ Diagonal Fock matrix in the MO basis -`scf_iteration `_ - Undocumented - `do_diis `_ If True, compute integrals on the fly diff --git a/src/Nuclei/README.rst b/src/Nuclei/README.rst index cf164529..048a7dd5 100644 --- a/src/Nuclei/README.rst +++ b/src/Nuclei/README.rst @@ -22,50 +22,50 @@ Documentation .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. -`nucl_charge `_ +`nucl_charge `_ Nuclear charges -`nucl_coord `_ +`nucl_coord `_ Nuclear coordinates in the format (:, {x,y,z}) -`nucl_coord_transp `_ +`nucl_coord_transp `_ Transposed array of nucl_coord -`nucl_dist `_ +`nucl_dist `_ nucl_dist : Nucleus-nucleus distances nucl_dist_2 : Nucleus-nucleus distances squared nucl_dist_vec : Nucleus-nucleus distances vectors -`nucl_dist_2 `_ +`nucl_dist_2 `_ nucl_dist : Nucleus-nucleus distances nucl_dist_2 : Nucleus-nucleus distances squared nucl_dist_vec : Nucleus-nucleus distances vectors -`nucl_dist_vec_x `_ +`nucl_dist_vec_x `_ nucl_dist : Nucleus-nucleus distances nucl_dist_2 : Nucleus-nucleus distances squared nucl_dist_vec : Nucleus-nucleus distances vectors -`nucl_dist_vec_y `_ +`nucl_dist_vec_y `_ nucl_dist : Nucleus-nucleus distances nucl_dist_2 : Nucleus-nucleus distances squared nucl_dist_vec : Nucleus-nucleus distances vectors -`nucl_dist_vec_z `_ +`nucl_dist_vec_z `_ nucl_dist : Nucleus-nucleus distances nucl_dist_2 : Nucleus-nucleus distances squared nucl_dist_vec : Nucleus-nucleus distances vectors -`nucl_label `_ +`nucl_label `_ Nuclear labels -`nucl_num `_ +`nucl_num `_ Number of nuclei -`nucl_num_aligned `_ +`nucl_num_aligned `_ Number of nuclei -`nuclear_repulsion `_ +`nuclear_repulsion `_ Nuclear repulsion energy diff --git a/src/Utils/README.rst b/src/Utils/README.rst index 7feee8ad..469a8ecc 100644 --- a/src/Utils/README.rst +++ b/src/Utils/README.rst @@ -10,42 +10,63 @@ Documentation .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. -`apply_rotation `_ +`apply_rotation `_ Apply the rotation found by find_rotation -`find_rotation `_ +`find_rotation `_ Find A.C = B -`get_pseudo_inverse `_ +`get_pseudo_inverse `_ Find C = A^-1 -`lapack_diag `_ +`lapack_diag `_ Diagonalize matrix H + .br + H is untouched between input and ouptut + .br + eigevalues(i) = ith lowest eigenvalue of the H matrix + .br + eigvectors(i,j) = where i is the basis function and psi_j is the j th eigenvector + .br -`ortho_lowdin `_ - Compute U.S^-1/2 canonical orthogonalization +`ortho_lowdin `_ + Compute C_new=C_old.S^-1/2 canonical orthogonalization. + .br + overlap : overlap matrix + .br + LDA : leftmost dimension of overlap array + .br + N : Overlap matrix is NxN (array is (LDA,N) ) + .br + C : Coefficients of the vectors to orthogonalize. On exit, + orthogonal vectors + .br + LDC : leftmost dimension of C + .br + m : Coefficients matrix is MxN, ( array is (LDC,N) ) + .br -`add_poly `_ +`add_poly `_ Add two polynomials D(t) =! D(t) +( B(t)+C(t)) -`add_poly_multiply `_ +`add_poly_multiply `_ Add a polynomial multiplied by a constant D(t) =! D(t) +( cst * B(t)) -`f_integral `_ +`f_integral `_ function that calculates the following integral \int_{\-infty}^{+\infty} x^n \exp(-p x^2) dx -`gaussian_product `_ +`gaussian_product `_ Gaussian product in 1D. e^{-a (x-x_A)^2} e^{-b (x-x_B)^2} = K_{ab}^x e^{-p (x-x_P)^2} -`gaussian_product_x `_ +`gaussian_product_x `_ Gaussian product in 1D. e^{-a (x-x_A)^2} e^{-b (x-x_B)^2} = K_{ab}^x e^{-p (x-x_P)^2} -`give_explicit_poly_and_gaussian `_ +`give_explicit_poly_and_gaussian `_ Transforms the product of (x-x_A)^a(1) (x-x_B)^b(1) (x-x_A)^a(2) (y-y_B)^b(2) (z-z_A)^a(3) (z-z_B)^b(3) exp(-(r-A)^2 alpha) exp(-(r-B)^2 beta) into @@ -53,105 +74,103 @@ Documentation * [ sum (l_y = 0,i_order(2)) P_new(l_y,2) * (y-P_center(2))^l_y ] exp (- p (y-P_center(2))^2 ) * [ sum (l_z = 0,i_order(3)) P_new(l_z,3) * (z-P_center(3))^l_z ] exp (- p (z-P_center(3))^2 ) -`give_explicit_poly_and_gaussian_x `_ +`give_explicit_poly_and_gaussian_x `_ Transform the product of (x-x_A)^a(1) (x-x_B)^b(1) (x-x_A)^a(2) (y-y_B)^b(2) (z-z_A)^a(3) (z-z_B)^b(3) exp(-(r-A)^2 alpha) exp(-(r-B)^2 beta) into fact_k (x-x_P)^iorder(1) (y-y_P)^iorder(2) (z-z_P)^iorder(3) exp(-p(r-P)^2) -`hermite `_ +`hermite `_ Hermite polynomial -`multiply_poly `_ +`multiply_poly `_ Multiply two polynomials D(t) =! D(t) +( B(t)*C(t)) -`recentered_poly2 `_ +`recentered_poly2 `_ Recenter two polynomials -`rint `_ +`rint `_ .. math:: .br \int_0^1 dx \exp(-p x^2) x^n .br -`rint1 `_ +`rint1 `_ Standard version of rint -`rint_large_n `_ +`rint_large_n `_ Version of rint for large values of n -`rint_sum `_ +`rint_sum `_ Needed for the calculation of two-electron integrals. -`overlap_gaussian_x `_ +`overlap_gaussian_x `_ .. math:: .br \sum_{-infty}^{+infty} (x-A_x)^ax (x-B_x)^bx exp(-alpha(x-A_x)^2) exp(-beta(x-B_X)^2) dx .br -`overlap_gaussian_xyz `_ +`overlap_gaussian_xyz `_ .. math:: .br S_x = \int (x-A_x)^{a_x} exp(-\alpha(x-A_x)^2) (x-B_x)^{b_x} exp(-beta(x-B_x)^2) dx \\ S = S_x S_y S_z .br -`overlap_x_abs `_ +`overlap_x_abs `_ .. math :: .br \int_{-infty}^{+infty} (x-A_center)^(power_A) * (x-B_center)^power_B * exp(-alpha(x-A_center)^2) * exp(-beta(x-B_center)^2) dx .br -`align_double `_ +`align_double `_ Compute 1st dimension such that it is aligned for vectorization. -`all_utils `_ +`all_utils `_ Dummy provider to provide all utils -`binom `_ +`binom `_ Binomial coefficients -`binom_func `_ +`binom_func `_ .. math :: .br \frac{i!}{j!(i-j)!} .br -`binom_transp `_ +`binom_transp `_ Binomial coefficients -`dble_fact `_ +`dble_fact `_ n!! -`fact `_ +`fact `_ n! -`fact_inv `_ +`fact_inv `_ 1/n! -`inv_int `_ +`inv_int `_ 1/i -`normalize `_ +`normalize `_ Normalizes vector u u is expected to be aligned in memory. -`nproc `_ +`nproc `_ Number of current OpenMP threads -`u_dot_u `_ +`u_dot_u `_ Compute - u is expected to be aligned in memory. -`u_dot_v `_ +`u_dot_v `_ Compute - u and v are expected to be aligned in memory. -`wall_time `_ +`wall_time `_ The equivalent of cpu_time, but for the wall time. -`write_git_log `_ +`write_git_log `_ Write the last git commit in file iunit. diff --git a/src/Utils/util.irp.f b/src/Utils/util.irp.f index e25148f1..47262616 100644 --- a/src/Utils/util.irp.f +++ b/src/Utils/util.irp.f @@ -174,7 +174,7 @@ BEGIN_PROVIDER [ double precision, inv_int, (128) ] ! 1/i END_DOC integer :: i - do i=1,size(inv_int) + do i=1,128 inv_int(i) = 1.d0/dble(i) enddo END_PROVIDER From 8860deb8b865163090c7b9f9650bb81455487e61 Mon Sep 17 00:00:00 2001 From: Manu Date: Thu, 22 May 2014 11:17:36 +0200 Subject: [PATCH 2/4] add Makefile.config.example --- src/Dets/README.rst | 159 +++++++++++++----- src/Dets/determinants.irp.f | 28 ++++ src/Dets/filter_connected.irp.f | 232 +++++++++++++++++++++++++- src/Dets/slater_rules.irp.f | 53 +++++- src/Makefile.config.example | 29 ++++ src/Perturbation/Moller_plesset.irp.f | 41 +++++ src/Perturbation/README.rst | 74 +++++++- src/Perturbation/epstein_nesbet.irp.f | 161 ++++++++++++++++++ 8 files changed, 733 insertions(+), 44 deletions(-) create mode 100644 src/Makefile.config.example create mode 100644 src/Perturbation/Moller_plesset.irp.f diff --git a/src/Dets/README.rst b/src/Dets/README.rst index ab562b17..bf16b740 100644 --- a/src/Dets/README.rst +++ b/src/Dets/README.rst @@ -50,28 +50,45 @@ Documentation .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. -`copy_h_apply_buffer_to_wf `_ +`copy_h_apply_buffer_to_wf `_ Undocumented -`h_apply_buffer_coef `_ - Buffer of determinants/coefficients for H_apply. Uninitialized. Filled by H_apply subroutines. +`h_apply_buffer_coef `_ + Buffer of determinants/coefficients/perturbative energy for H_apply. + Uninitialized. Filled by H_apply subroutines. -`h_apply_buffer_det `_ - Buffer of determinants/coefficients for H_apply. Uninitialized. Filled by H_apply subroutines. +`h_apply_buffer_det `_ + Buffer of determinants/coefficients/perturbative energy for H_apply. + Uninitialized. Filled by H_apply subroutines. -`h_apply_buffer_n_det `_ - Buffer of determinants/coefficients for H_apply. Uninitialized. Filled by H_apply subroutines. +`h_apply_buffer_e2 `_ + Buffer of determinants/coefficients/perturbative energy for H_apply. + Uninitialized. Filled by H_apply subroutines. -`h_apply_buffer_size `_ +`h_apply_buffer_n_det `_ + Buffer of determinants/coefficients/perturbative energy for H_apply. + Uninitialized. Filled by H_apply subroutines. + +`h_apply_buffer_size `_ Size of the H_apply buffer. -`h_apply_threshold `_ +`h_apply_threshold `_ Theshold on | | -`resize_h_apply_buffer_det `_ +`resize_h_apply_buffer_det `_ Undocumented -`davidson_diag `_ +`connected_to_ref `_ + Undocumented + +`det_is_not_or_may_be_in_ref `_ + If true, det is not in ref + If false, det may be in ref + +`key_pattern_not_in_ref `_ + Min and max values of the integers of the keys of the reference + +`davidson_diag `_ Davidson diagonalization. .br dets_in : bitmasks corresponding to determinants @@ -87,102 +104,164 @@ Documentation .br Initial guess vectors are not necessarily orthonormal -`davidson_iter_max `_ +`davidson_iter_max `_ Max number of Davidson iterations -`davidson_sze_max `_ +`davidson_sze_max `_ Max number of Davidson sizes -`n_det `_ +`n_det `_ Number of determinants in the wave function -`n_det_generators `_ +`n_det_generators `_ Number of generator determinants in the wave function -`n_states `_ +`n_states `_ Number of states to consider -`psi_coef `_ +`psi_coef `_ The wave function. Initialized with Hartree-Fock -`psi_det `_ +`psi_det `_ The wave function. Initialized with Hartree-Fock -`psi_det_size `_ +`psi_det_size `_ Size of the psi_det/psi_coef arrays -`psi_generators `_ +`psi_generators `_ Determinants on which H is applied -`double_exc_bitmask `_ +`psi_ref `_ + Determinants on which H is applied + +`psi_ref_coef `_ + Determinants on which H is applied + +`psi_ref_size `_ + Number of generator determinants in the wave function + +`double_exc_bitmask `_ double_exc_bitmask(:,1,i) is the bitmask for holes of excitation 1 double_exc_bitmask(:,2,i) is the bitmask for particles of excitation 1 double_exc_bitmask(:,3,i) is the bitmask for holes of excitation 2 double_exc_bitmask(:,4,i) is the bitmask for particles of excitation 2 for a given couple of hole/particle excitations i. -`n_double_exc_bitmasks `_ +`n_double_exc_bitmasks `_ Number of double excitation bitmasks -`n_single_exc_bitmasks `_ +`n_single_exc_bitmasks `_ Number of single excitation bitmasks -`single_exc_bitmask `_ +`single_exc_bitmask `_ single_exc_bitmask(:,1,i) is the bitmask for holes single_exc_bitmask(:,2,i) is the bitmask for particles for a given couple of hole/particle excitations i. -`get_s2 `_ +`filter_connected `_ + Filters out the determinants that are not connected by H + .br + returns the array idx which contains the index of the + .br + determinants in the array key1 that interact + .br + via the H operator with key2. + .br + idx(0) is the number of determinants that interact with key1 + +`filter_connected_i_h_psi0 `_ + returns the array idx which contains the index of the + .br + determinants in the array key1 that interact + .br + via the H operator with key2. + .br + idx(0) is the number of determinants that interact with key1 + +`filter_connected_i_h_psi0_sc2 `_ + standard filter_connected_i_H_psi but returns in addition + .br + the array of the index of the non connected determinants to key1 + .br + in order to know what double excitation can be repeated on key1 + .br + idx_repeat(0) is the number of determinants that can be used + .br + to repeat the excitations + +`get_s2 `_ Returns -`a_operator `_ +`get_s2_u0 `_ + Undocumented + +`s_z `_ + Undocumented + +`s_z2_sz `_ + Undocumented + +`a_operator `_ Needed for diag_H_mat_elem -`ac_operator `_ +`ac_operator `_ Needed for diag_H_mat_elem -`decode_exc `_ +`decode_exc `_ Decodes the exc arrays returned by get_excitation. h1,h2 : Holes p1,p2 : Particles s1,s2 : Spins (1:alpha, 2:beta) degree : Degree of excitation -`diag_h_mat_elem `_ +`diag_h_mat_elem `_ Computes -`get_double_excitation `_ +`get_double_excitation `_ Returns the two excitation operators between two doubly excited determinants and the phase -`get_excitation `_ +`get_excitation `_ Returns the excitation operators between two determinants and the phase -`get_excitation_degree `_ +`get_excitation_degree `_ Returns the excitation degree between two determinants -`get_excitation_degree_vector `_ +`get_excitation_degree_vector `_ Applies get_excitation_degree to an array of determinants -`get_mono_excitation `_ +`get_mono_excitation `_ Returns the excitation operator between two singly excited determinants and the phase -`get_occ_from_key `_ +`get_occ_from_key `_ Returns a list of occupation numbers from a bitstring -`h_u_0 `_ +`h_u_0 `_ Computes v_0 = H|u_0> .br n : number of determinants .br H_jj : array of -`i_h_j `_ +`i_h_j `_ Returns where i and j are determinants -`i_h_psi `_ - Undocumented +`i_h_psi `_ + for the various Nstate -`h_matrix_all_dets `_ +`i_h_psi_sc2 `_ + for the various Nstate + .br + returns in addition + .br + the array of the index of the non connected determinants to key1 + .br + in order to know what double excitation can be repeated on key1 + .br + idx_repeat(0) is the number of determinants that can be used + .br + to repeat the excitations + +`h_matrix_all_dets `_ H matrix on the basis of the slater deter;inants defined by psi_det diff --git a/src/Dets/determinants.irp.f b/src/Dets/determinants.irp.f index 1cb89644..b46fdaaf 100644 --- a/src/Dets/determinants.irp.f +++ b/src/Dets/determinants.irp.f @@ -75,3 +75,31 @@ BEGIN_PROVIDER [ integer(bit_kind), psi_generators, (N_int,2,psi_det_size) ] END_PROVIDER +BEGIN_PROVIDER [ integer, psi_ref_size] + implicit none + BEGIN_DOC + ! Number of generator determinants in the wave function + END_DOC + psi_det_size = N_det +END_PROVIDER + + BEGIN_PROVIDER [ integer(bit_kind), psi_ref, (N_int,2,psi_ref_size) ] +&BEGIN_PROVIDER [ double precision, psi_ref_coef, (psi_ref_size,N_states) ] + implicit none + BEGIN_DOC + ! Determinants on which H is applied + END_DOC + integer :: i,k + + do k = 1, psi_ref_size + do i=1,N_int + psi_ref(i,1,k) = psi_det(i,1,k) + psi_ref(i,2,k) = psi_det(i,1,k) + enddo + do i = 1, N_states + psi_ref_coef(k,i) = psi_coef(k,i) + enddo + enddo + +END_PROVIDER + diff --git a/src/Dets/filter_connected.irp.f b/src/Dets/filter_connected.irp.f index 0353c9c6..d6caa977 100644 --- a/src/Dets/filter_connected.irp.f +++ b/src/Dets/filter_connected.irp.f @@ -3,7 +3,16 @@ subroutine filter_connected(key1,key2,Nint,sze,idx) use bitmasks implicit none BEGIN_DOC -! Filters out the determinants that are not connected by H + ! Filters out the determinants that are not connected by H + ! + ! returns the array idx which contains the index of the + ! + ! determinants in the array key1 that interact + ! + ! via the H operator with key2. + ! + ! idx(0) is the number of determinants that interact with key1 + END_DOC END_DOC integer, intent(in) :: Nint, sze integer(bit_kind), intent(in) :: key1(Nint,2,sze) @@ -85,6 +94,15 @@ end subroutine filter_connected_i_H_psi0(key1,key2,Nint,sze,idx) use bitmasks + BEGIN_DOC + ! returns the array idx which contains the index of the + ! + ! determinants in the array key1 that interact + ! + ! via the H operator with key2. + ! + ! idx(0) is the number of determinants that interact with key1 + END_DOC implicit none integer, intent(in) :: Nint, sze integer(bit_kind), intent(in) :: key1(Nint,2,sze) @@ -173,3 +191,215 @@ subroutine filter_connected_i_H_psi0(key1,key2,Nint,sze,idx) idx(0) = l-1 end +subroutine filter_connected_i_H_psi0_SC2(key1,key2,Nint,sze,idx,idx_repeat) + use bitmasks + BEGIN_DOC + ! standard filter_connected_i_H_psi but returns in addition + ! + ! the array of the index of the non connected determinants to key1 + ! + ! in order to know what double excitation can be repeated on key1 + ! + ! idx_repeat(0) is the number of determinants that can be used + ! + ! to repeat the excitations + END_DOC + implicit none + integer, intent(in) :: Nint, sze + integer(bit_kind), intent(in) :: key1(Nint,2,sze) + integer(bit_kind), intent(in) :: key2(Nint,2) + integer, intent(out) :: idx(0:sze) + integer, intent(out) :: idx_repeat(0:sze) + + integer :: i,l,l_repeat + integer :: degree_x2 + + ASSERT (Nint > 0) + ASSERT (Nint == N_int) + ASSERT (sze > 0) + + l=1 + l_repeat=1 + call get_excitation_degree(ref_bitmask,key2,degree,Nint) + integer :: degree + + if(degree == 2)then + if (Nint==1) then + + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + & + popcnt(xor( key1(1,2,i), key2(1,2))) + if (degree_x2 < 5) then + if(degree_x2 .ne. 0)then + idx(l) = i + l = l+1 + endif + elseif(degree>6)then + idx_repeat(l_repeat) = i + l_repeat = l_repeat + 1 + endif + enddo + + else if (Nint==2) then + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + & + popcnt(xor( key1(2,1,i), key2(2,1))) + & + popcnt(xor( key1(1,2,i), key2(1,2))) + & + popcnt(xor( key1(2,2,i), key2(2,2))) + if (degree_x2 < 5) then + if(degree_x2 .ne. 0)then + idx(l) = i + l = l+1 + endif + elseif(degree>6)then + idx_repeat(l_repeat) = i + l_repeat = l_repeat + 1 + endif + enddo + + else if (Nint==3) then + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + & + popcnt(xor( key1(1,2,i), key2(1,2))) + & + popcnt(xor( key1(2,1,i), key2(2,1))) + & + popcnt(xor( key1(2,2,i), key2(2,2))) + & + popcnt(xor( key1(3,1,i), key2(3,1))) + & + popcnt(xor( key1(3,2,i), key2(3,2))) + if (degree_x2 < 5) then + if(degree_x2 .ne. 0)then + idx(l) = i + l = l+1 + endif + elseif(degree>6)then + idx_repeat(l_repeat) = i + l_repeat = l_repeat + 1 + endif + enddo + + else + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = 0 + !DEC$ LOOP COUNT MIN(4) + do l=1,Nint + degree_x2 = degree_x2+ popcnt(xor( key1(l,1,i), key2(l,1))) +& + popcnt(xor( key1(l,2,i), key2(l,2))) + if (degree_x2 > 4) then + exit + endif + enddo + if (degree_x2 <= 5) then + if(degree_x2 .ne. 0)then + idx(l) = i + l = l+1 + endif + elseif(degree>6)then + idx_repeat(l_repeat) = i + l_repeat = l_repeat + 1 + endif + enddo + + endif + elseif(degree==1)then + if (Nint==1) then + + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + & + popcnt(xor( key1(1,2,i), key2(1,2))) + if (degree_x2 < 5) then + if(degree_x2 .ne. 0)then + idx(l) = i + l = l+1 + endif + else + idx_repeat(l_repeat) = i + l_repeat = l_repeat + 1 + endif + enddo + + else if (Nint==2) then + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + & + popcnt(xor( key1(2,1,i), key2(2,1))) + & + popcnt(xor( key1(1,2,i), key2(1,2))) + & + popcnt(xor( key1(2,2,i), key2(2,2))) + if (degree_x2 < 5) then + if(degree_x2 .ne. 0)then + idx(l) = i + l = l+1 + endif + else + idx_repeat(l_repeat) = i + l_repeat = l_repeat + 1 + endif + enddo + + else if (Nint==3) then + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + & + popcnt(xor( key1(1,2,i), key2(1,2))) + & + popcnt(xor( key1(2,1,i), key2(2,1))) + & + popcnt(xor( key1(2,2,i), key2(2,2))) + & + popcnt(xor( key1(3,1,i), key2(3,1))) + & + popcnt(xor( key1(3,2,i), key2(3,2))) + if (degree_x2 < 5) then + if(degree_x2 .ne. 0)then + idx(l) = i + l = l+1 + endif + else + idx_repeat(l_repeat) = i + l_repeat = l_repeat + 1 + endif + enddo + + else + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = 0 + !DEC$ LOOP COUNT MIN(4) + do l=1,Nint + degree_x2 = degree_x2+ popcnt(xor( key1(l,1,i), key2(l,1))) +& + popcnt(xor( key1(l,2,i), key2(l,2))) + if (degree_x2 > 4) then + exit + endif + enddo + if (degree_x2 <= 5) then + if(degree_x2 .ne. 0)then + idx(l) = i + l = l+1 + endif + else + idx_repeat(l_repeat) = i + l_repeat = l_repeat + 1 + endif + enddo + + endif + + else + print*,'more than a double excitation, can not apply the ' + print*,'SC2 dressing of the diagonal element .....' + print*,'stop !!' + print*,'degree = ',degree + stop + endif + idx(0) = l-1 + idx_repeat(0) = l_repeat-1 +end + diff --git a/src/Dets/slater_rules.irp.f b/src/Dets/slater_rules.irp.f index 3f56eba3..3599b9a0 100644 --- a/src/Dets/slater_rules.irp.f +++ b/src/Dets/slater_rules.irp.f @@ -502,6 +502,9 @@ subroutine i_H_psi(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array) integer :: exc(0:2,2,2) double precision :: hij integer :: idx(0:Ndet) + BEGIN_DOC + ! for the various Nstate + END_DOC ASSERT (Nint > 0) ASSERT (N_int == Nint) @@ -517,7 +520,55 @@ subroutine i_H_psi(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array) do j = 1, Nstate i_H_psi_array(j) = i_H_psi_array(j) + coef(i,j)*hij enddo - print *, 'x', coef(i,1), hij, i_H_psi_array(1) +! print *, 'x', coef(i,1), hij, i_H_psi_array(1) + enddo +end + + +subroutine i_H_psi_SC2(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array,idx_repeat) + use bitmasks + BEGIN_DOC + ! for the various Nstate + ! + ! returns in addition + ! + ! the array of the index of the non connected determinants to key1 + ! + ! in order to know what double excitation can be repeated on key1 + ! + ! idx_repeat(0) is the number of determinants that can be used + ! + ! to repeat the excitations + END_DOC + implicit none + integer, intent(in) :: Nint, Ndet,Ndet_max,Nstate + integer(bit_kind), intent(in) :: keys(Nint,2,Ndet) + integer(bit_kind), intent(in) :: key(Nint,2) + double precision, intent(in) :: coef(Ndet_max,Nstate) + double precision, intent(out) :: i_H_psi_array(Nstate) + integer , intent(out) :: idx_repeat(0:Ndet) + + integer :: i, ii,j + double precision :: phase + integer :: exc(0:2,2,2) + double precision :: hij + integer :: idx(0:Ndet) + + ASSERT (Nint > 0) + ASSERT (N_int == Nint) + ASSERT (Nstate > 0) + ASSERT (Ndet > 0) + ASSERT (Ndet_max >= Ndet) + i_H_psi_array = 0.d0 + call filter_connected_i_H_psi0_SC2(keys,key,Nint,Ndet,idx,idx_repeat) + do ii=1,idx(0) + i = idx(ii) + !DEC$ FORCEINLINE + call i_H_j(keys(1,1,i),key,Nint,hij) + do j = 1, Nstate + i_H_psi_array(j) = i_H_psi_array(j) + coef(i,j)*hij + enddo +! print *, 'x', coef(i,1), hij, i_H_psi_array(1) enddo end diff --git a/src/Makefile.config.example b/src/Makefile.config.example new file mode 100644 index 00000000..f4cbb568 --- /dev/null +++ b/src/Makefile.config.example @@ -0,0 +1,29 @@ +OPENMP =1 +PROFILE =0 +DEBUG = 0 + +IRPF90_FLAGS+= --align=32 +FC = ifort -g +FCFLAGS= +FCFLAGS+= -xHost +#FCFLAGS+= -xAVX +FCFLAGS+= -O2 +FCFLAGS+= -ip +FCFLAGS+= -opt-prefetch +FCFLAGS+= -ftz +MKL=-mkl=parallel + +ifeq ($(PROFILE),1) +FC += -p -g +CXX += -pg +endif + +ifeq ($(OPENMP),1) +FC += -openmp +CXX += -fopenmp +endif + +ifeq ($(DEBUG),1) +FC += -C -traceback -fpe0 +#FCFLAGS =-O0 +endif diff --git a/src/Perturbation/Moller_plesset.irp.f b/src/Perturbation/Moller_plesset.irp.f new file mode 100644 index 00000000..4f99e9d2 --- /dev/null +++ b/src/Perturbation/Moller_plesset.irp.f @@ -0,0 +1,41 @@ +subroutine pt2_moller_plesset(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_st) + use bitmasks + implicit none + integer, intent(in) :: Nint,ndet,n_st + integer(bit_kind), intent(in) :: det_pert(Nint,2) + double precision , intent(out) :: c_pert(n_st),e_2_pert(n_st),H_pert_diag + double precision :: i_H_psi_array(N_st) + + BEGIN_DOC + ! compute the standard Moller-Plesset perturbative first order coefficient and second order energetic contribution + ! + ! for the various n_st states. + ! + ! c_pert(i) = /(difference of orbital energies) + ! + ! e_2_pert(i) = ^2/(difference of orbital energies) + ! + END_DOC + + integer :: i,j + double precision :: diag_H_mat_elem + integer :: exc(0:2,2,2) + integer :: degree + double precision :: phase,delta_e + integer :: h1,h2,p1,p2,s1,s2 + ASSERT (Nint == N_int) + ASSERT (Nint > 0) + call get_excitation(det_pert,ref_bitmask,exc,degree,phase,Nint) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + delta_e = Fock_matrix_diag_mo(h1) + Fock_matrix_diag_mo(h2) - & + Fock_matrix_diag_mo(p1) + Fock_matrix_diag_mo(p2) + delta_e = 1.d0/delta_e + + call i_H_psi(det_pert,psi_ref,psi_ref_coef,Nint,ndet,psi_ref_size,n_st,i_H_psi_array) + H_pert_diag = diag_H_mat_elem(det_pert,Nint) + do i =1,n_st + c_pert(i) = i_H_psi_array(i) *delta_e + e_2_pert(i) = c_pert(i) * i_H_psi_array(i) + enddo + +end diff --git a/src/Perturbation/README.rst b/src/Perturbation/README.rst index 83392b83..3acd8ef6 100644 --- a/src/Perturbation/README.rst +++ b/src/Perturbation/README.rst @@ -82,7 +82,7 @@ Documentation .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. -`pt2_epstein_nesbet `_ +`pt2_epstein_nesbet `_ compute the standard Epstein-Nesbet perturbative first order coefficient and second order energetic contribution .br for the various n_st states. @@ -92,7 +92,7 @@ Documentation e_2_pert(i) = ^2/( E(i) - ) .br -`pt2_epstein_nesbet_2x2 `_ +`pt2_epstein_nesbet_2x2 `_ compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution .br for the various n_st states. @@ -102,6 +102,76 @@ Documentation c_pert(i) = e_2_pert(i)/ .br +`pt2_epstein_nesbet_2x2_sc2 `_ + compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution + .br + for the various n_st states. + .br + but with the correction in the denominator + .br + comming from the interaction of that determinant with all the others determinants + .br + that can be repeated by repeating all the double excitations + .br + : you repeat all the correlation energy already taken into account in reference_energy(1) + .br + that could be repeated to this determinant. + .br + ---> + delta_e_corr + .br + e_2_pert(i) = 0.5 * (( - E(i) ) - sqrt( ( - E(i)) ^2 + 4 ^2 ) + .br + c_pert(i) = e_2_pert(i)/ + .br + +`pt2_epstein_nesbet_sc2 `_ + compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution + .br + for the various n_st states, + .br + but with the correction in the denominator + .br + comming from the interaction of that determinant with all the others determinants + .br + that can be repeated by repeating all the double excitations + .br + : you repeat all the correlation energy already taken into account in reference_energy(1) + .br + that could be repeated to this determinant. + .br + ---> + delta_e_corr + .br + c_pert(i) = /( E(i) - ( ) ) + .br + e_2_pert(i) = ^2/( E(i) - ( ) ) + .br + +`pt2_epstein_nesbet_sc2_projected `_ + compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution + .br + for the various n_st states, + .br + but with the correction in the denominator + .br + comming from the interaction of that determinant with all the others determinants + .br + that can be repeated by repeating all the double excitations + .br + : you repeat all the correlation energy already taken into account in reference_energy(1) + .br + that could be repeated to this determinant. + .br + BUT on the contrary with ""pt2_epstein_nesbet_SC2"", you compute the energy by projection + .br + ---> + delta_e_corr + .br + c_pert(1) = 1/c_HF /( E(i) - ( ) ) + .br + e_2_pert(1) = c_pert(1) + .br + NOTE :::: if you satisfy Brillouin Theorem, the singles don't contribute !! + .br + Needed Modules diff --git a/src/Perturbation/epstein_nesbet.irp.f b/src/Perturbation/epstein_nesbet.irp.f index 4c1e7c45..37572ae5 100644 --- a/src/Perturbation/epstein_nesbet.irp.f +++ b/src/Perturbation/epstein_nesbet.irp.f @@ -66,3 +66,164 @@ subroutine pt2_epstein_nesbet_2x2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet print *, e_2_pert, delta_e , i_H_psi_array end + +subroutine pt2_epstein_nesbet_SC2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_st) + use bitmasks + implicit none + integer, intent(in) :: Nint,ndet,n_st + integer(bit_kind), intent(in) :: det_pert(Nint,2) + double precision , intent(out) :: c_pert(n_st),e_2_pert(n_st),H_pert_diag + double precision :: i_H_psi_array(N_st) + integer :: idx_repeat(ndet) + + BEGIN_DOC + ! compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution + ! + ! for the various n_st states, + ! + ! but with the correction in the denominator + ! + ! comming from the interaction of that determinant with all the others determinants + ! + ! that can be repeated by repeating all the double excitations + ! + ! : you repeat all the correlation energy already taken into account in reference_energy(1) + ! + ! that could be repeated to this determinant. + ! + ! ---> + delta_e_corr + ! + ! c_pert(i) = /( E(i) - ( ) ) + ! + ! e_2_pert(i) = ^2/( E(i) - ( ) ) + ! + END_DOC + + integer :: i,j + double precision :: diag_H_mat_elem,accu_e_corr,hij + ASSERT (Nint == N_int) + ASSERT (Nint > 0) + call i_H_psi_SC2(det_pert,psi_ref,psi_ref_coef,Nint,ndet,psi_ref_size,n_st,i_H_psi_array,idx_repeat) + accu_e_corr = 0.d0 + do i = 1, idx_repeat(0) + call i_H_j(psi_ref(1,1,idx_repeat(i)),det_pert,Nint,hij) + accu_e_corr = accu_e_corr + hij * psi_ref_coef(idx_repeat(i)) + enddo + accu_e_corr = accu_e_corr / psi_ref_coef(1) + H_pert_diag = diag_H_mat_elem(det_pert,Nint) + accu_e_corr + do i =1,n_st + c_pert(i) = i_H_psi_array(i) / (reference_energy(i) - H_pert_diag) + e_2_pert(i) = c_pert(i) * i_H_psi_array(i) + enddo + +end +subroutine pt2_epstein_nesbet_2x2_SC2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_st) + use bitmasks + implicit none + integer, intent(in) :: Nint,ndet,n_st + integer(bit_kind), intent(in) :: det_pert(Nint,2) + double precision , intent(out) :: c_pert(n_st),e_2_pert(n_st),H_pert_diag + double precision :: i_H_psi_array(N_st) + integer :: idx_repeat(ndet) + + BEGIN_DOC + ! compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution + ! + ! for the various n_st states. + ! + ! but with the correction in the denominator + ! + ! comming from the interaction of that determinant with all the others determinants + ! + ! that can be repeated by repeating all the double excitations + ! + ! : you repeat all the correlation energy already taken into account in reference_energy(1) + ! + ! that could be repeated to this determinant. + ! + ! ---> + delta_e_corr + ! + ! e_2_pert(i) = 0.5 * (( - E(i) ) - sqrt( ( - E(i)) ^2 + 4 ^2 ) + ! + ! c_pert(i) = e_2_pert(i)/ + ! + END_DOC + + integer :: i,j + double precision :: diag_H_mat_elem,accu_e_corr,hij,delta_e + ASSERT (Nint == N_int) + ASSERT (Nint > 0) + call i_H_psi_SC2(det_pert,psi_ref,psi_ref_coef,Nint,ndet,psi_ref_size,n_st,i_H_psi_array,idx_repeat) + accu_e_corr = 0.d0 + do i = 1, idx_repeat(0) + call i_H_j(psi_ref(1,1,idx_repeat(i)),det_pert,Nint,hij) + accu_e_corr = accu_e_corr + hij * psi_ref_coef(idx_repeat(i)) + enddo + accu_e_corr = accu_e_corr / psi_ref_coef(1) + H_pert_diag = diag_H_mat_elem(det_pert,Nint) + accu_e_corr + do i =1,n_st + delta_e = H_pert_diag - reference_energy(i) + e_2_pert(i) = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * i_H_psi_array(i) * i_H_psi_array(i))) + c_pert(i) = e_2_pert(i)/i_H_psi_array(i) + enddo + +end + +subroutine pt2_epstein_nesbet_SC2_projected(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_st) + use bitmasks + implicit none + integer, intent(in) :: Nint,ndet,n_st + integer(bit_kind), intent(in) :: det_pert(Nint,2) + double precision , intent(out) :: c_pert(n_st),e_2_pert(n_st),H_pert_diag + double precision :: i_H_psi_array(N_st) + integer :: idx_repeat(ndet) + + BEGIN_DOC + ! compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution + ! + ! for the various n_st states, + ! + ! but with the correction in the denominator + ! + ! comming from the interaction of that determinant with all the others determinants + ! + ! that can be repeated by repeating all the double excitations + ! + ! : you repeat all the correlation energy already taken into account in reference_energy(1) + ! + ! that could be repeated to this determinant. + ! + ! BUT on the contrary with ""pt2_epstein_nesbet_SC2"", you compute the energy by projection + ! + ! ---> + delta_e_corr + ! + ! c_pert(1) = 1/c_HF /( E(i) - ( ) ) + ! + ! e_2_pert(1) = c_pert(1) + ! + ! NOTE :::: if you satisfy Brillouin Theorem, the singles don't contribute !! + ! + END_DOC + + integer :: i,j + double precision :: diag_H_mat_elem,accu_e_corr,hij,h0j + ASSERT (Nint == N_int) + ASSERT (Nint > 0) + call i_H_psi_SC2(det_pert,psi_ref,psi_ref_coef,Nint,ndet,psi_ref_size,n_st,i_H_psi_array,idx_repeat) + accu_e_corr = 0.d0 + call i_H_j(ref_bitmask,det_pert,Nint,h0j) + do i = 1, idx_repeat(0) + call i_H_j(psi_ref(1,1,idx_repeat(i)),det_pert,Nint,hij) + accu_e_corr = accu_e_corr + hij * psi_ref_coef(idx_repeat(i)) + enddo + accu_e_corr = accu_e_corr / psi_ref_coef(1) + H_pert_diag = diag_H_mat_elem(det_pert,Nint) + accu_e_corr + + c_pert(1) = 1.d0/psi_ref_coef(1) * i_H_psi_array(1) / (reference_energy(i) - H_pert_diag) + e_2_pert(1) = c_pert(i) * h0j + do i =2,n_st + c_pert(i) = i_H_psi_array(i) / (reference_energy(i) - H_pert_diag) + e_2_pert(i) = c_pert(i) * i_H_psi_array(i) + enddo + +end From 006e0affd08cd041001a3d581326df293c4825cb Mon Sep 17 00:00:00 2001 From: Manu Date: Thu, 22 May 2014 11:23:45 +0200 Subject: [PATCH 3/4] added NEEDED_MODULES Makefile for the CIS_dressed --- src/CIS_dressed/Makefile | 8 ++++++++ src/CIS_dressed/NEEDED_MODULES | 1 + 2 files changed, 9 insertions(+) create mode 100644 src/CIS_dressed/Makefile create mode 100644 src/CIS_dressed/NEEDED_MODULES diff --git a/src/CIS_dressed/Makefile b/src/CIS_dressed/Makefile new file mode 100644 index 00000000..b2ea1de1 --- /dev/null +++ b/src/CIS_dressed/Makefile @@ -0,0 +1,8 @@ +default: all + +# Define here all new external source files and objects.Don't forget to prefix the +# object files with IRPF90_temp/ +SRC= +OBJ= + +include $(QPACKAGE_ROOT)/src/Makefile.common diff --git a/src/CIS_dressed/NEEDED_MODULES b/src/CIS_dressed/NEEDED_MODULES new file mode 100644 index 00000000..e81cf251 --- /dev/null +++ b/src/CIS_dressed/NEEDED_MODULES @@ -0,0 +1 @@ +AOs BiInts Bitmask Dets Electrons Ezfio_files Hartree_Fock MonoInts MOs Nuclei Output Utils DensityMatrix From 1ffbeb616738a46ae18073832a935e9157334329 Mon Sep 17 00:00:00 2001 From: Manu Date: Thu, 22 May 2014 11:24:36 +0200 Subject: [PATCH 4/4] add ASSUMPTION.rst --- src/CIS_dressed/ASSUMPTIONS.rst | 0 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 src/CIS_dressed/ASSUMPTIONS.rst diff --git a/src/CIS_dressed/ASSUMPTIONS.rst b/src/CIS_dressed/ASSUMPTIONS.rst new file mode 100644 index 00000000..e69de29b