diff --git a/src/Dets/README.rst b/src/Dets/README.rst index 35c61523..09352c44 100644 --- a/src/Dets/README.rst +++ b/src/Dets/README.rst @@ -50,24 +50,24 @@ 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 `_ Copies the H_apply buffer to psi_coef. You need to touch psi_det, psi_coef and N_det after calling this function. -`fill_h_apply_buffer_no_selection `_ +`fill_h_apply_buffer_no_selection `_ Fill the H_apply buffer with determiants for CISD -`h_apply_buffer_allocated `_ +`h_apply_buffer_allocated `_ Buffer of determinants/coefficients/perturbative energy for H_apply. Uninitialized. Filled by H_apply subroutines. -`h_apply_threshold `_ +`h_apply_threshold `_ Theshold on | | -`resize_h_apply_buffer `_ +`resize_h_apply_buffer `_ Undocumented -`cisd_sc2 `_ +`cisd_sc2 `_ CISD+SC2 method :: take off all the disconnected terms of a CISD (selected or not) .br dets_in : bitmasks corresponding to determinants @@ -83,29 +83,29 @@ Documentation .br Initial guess vectors are not necessarily orthonormal -`repeat_excitation `_ +`repeat_excitation `_ Undocumented -`connected_to_ref `_ +`connected_to_ref `_ Undocumented -`det_is_not_or_may_be_in_ref `_ +`det_is_not_or_may_be_in_ref `_ If true, det is not in ref If false, det may be in ref -`is_in_wavefunction `_ +`is_in_wavefunction `_ Undocumented -`key_pattern_not_in_ref `_ +`key_pattern_not_in_ref `_ Min and max values of the integers of the keys of the reference -`davidson_converged `_ +`davidson_converged `_ True if the Davidson algorithm is converged -`davidson_criterion `_ +`davidson_criterion `_ Can be : [ energy | residual | both | wall_time | cpu_time | iterations ] -`davidson_diag `_ +`davidson_diag `_ Davidson diagonalization. .br dets_in : bitmasks corresponding to determinants @@ -123,7 +123,7 @@ Documentation .br Initial guess vectors are not necessarily orthonormal -`davidson_diag_hjj `_ +`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 @@ -143,114 +143,114 @@ 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 -`davidson_threshold `_ +`davidson_threshold `_ Can be : [ energy | residual | both | wall_time | cpu_time | iterations ] -`det_search_key `_ +`det_search_key `_ Return an integer*8 corresponding to a determinant index for searching -`n_det `_ +`n_det `_ Number of determinants in the wave function -`n_det_max_jacobi `_ +`n_det_max_jacobi `_ Maximum number of determinants diagonalized my jacobi -`n_states `_ +`n_states `_ Number of states to consider -`psi_average_norm_contrib `_ +`psi_average_norm_contrib `_ Contribution of determinants to the state-averaged density -`psi_average_norm_contrib_sorted `_ +`psi_average_norm_contrib_sorted `_ Wave function sorted by determinants contribution to the norm (state-averaged) -`psi_coef `_ +`psi_coef `_ The wave function coefficients. Initialized with Hartree-Fock if the EZFIO file is empty -`psi_coef_sorted `_ +`psi_coef_sorted `_ Wave function sorted by determinants contribution to the norm (state-averaged) -`psi_coef_sorted_bit `_ +`psi_coef_sorted_bit `_ Determinants on which we apply for perturbation. o They are sorted by determinants interpreted as integers. Useful to accelerate the search of a determinant -`psi_det `_ +`psi_det `_ The wave function determinants. Initialized with Hartree-Fock if the EZFIO file is empty -`psi_det_size `_ +`psi_det_size `_ Size of the psi_det/psi_coef arrays -`psi_det_sorted `_ +`psi_det_sorted `_ Wave function sorted by determinants contribution to the norm (state-averaged) -`psi_det_sorted_bit `_ +`psi_det_sorted_bit `_ Determinants on which we apply for perturbation. o They are sorted by determinants interpreted as integers. Useful to accelerate the search of a determinant -`read_dets `_ +`read_dets `_ Reads the determinants from the EZFIO file -`save_wavefunction `_ +`save_wavefunction `_ Save the wave function into the EZFIO file -`double_exc_bitmask `_ +`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. -`ci_eigenvectors `_ +`ci_eigenvectors `_ Eigenvectors/values of the CI matrix -`ci_electronic_energy `_ +`ci_electronic_energy `_ Eigenvectors/values of the CI matrix -`ci_energy `_ +`ci_energy `_ N_states lowest eigenvalues of the CI matrix -`diag_algorithm `_ +`diag_algorithm `_ Diagonalization algorithm (Davidson or Lapack) -`diagonalize_ci `_ +`diagonalize_ci `_ Replace the coefficients of the CI states by the coefficients of the eigenstates of the CI matrix -`ci_sc2_eigenvectors `_ +`ci_sc2_eigenvectors `_ Eigenvectors/values of the CI matrix -`ci_sc2_electronic_energy `_ +`ci_sc2_electronic_energy `_ Eigenvectors/values of the CI matrix -`ci_sc2_energy `_ +`ci_sc2_energy `_ N_states lowest eigenvalues of the CI matrix -`diagonalize_ci_sc2 `_ +`diagonalize_ci_sc2 `_ Replace the coefficients of the CI states by the coefficients of the eigenstates of the CI matrix -`filter_connected `_ +`filter_connected `_ Filters out the determinants that are not connected by H .br returns the array idx which contains the index of the @@ -261,7 +261,7 @@ Documentation .br idx(0) is the number of determinants that interact with key1 -`filter_connected_davidson `_ +`filter_connected_davidson `_ Filters out the determinants that are not connected by H .br returns the array idx which contains the index of the @@ -272,7 +272,7 @@ Documentation .br idx(0) is the number of determinants that interact with key1 -`filter_connected_i_h_psi0 `_ +`filter_connected_i_h_psi0 `_ returns the array idx which contains the index of the .br determinants in the array key1 that interact @@ -281,7 +281,7 @@ Documentation .br idx(0) is the number of determinants that interact with key1 -`filter_connected_i_h_psi0_sc2 `_ +`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 @@ -292,69 +292,69 @@ Documentation .br to repeat the excitations -`get_s2 `_ +`get_s2 `_ Returns -`get_s2_u0 `_ +`get_s2_u0 `_ Undocumented -`s_z `_ +`s_z `_ Undocumented -`s_z2_sz `_ +`s_z2_sz `_ Undocumented -`a_operator `_ +`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 -`det_connections `_ +`det_connections `_ .br -`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 `_ +`i_h_psi `_ for the various Nstates -`i_h_psi_sc2 `_ +`i_h_psi_sc2 `_ for the various Nstate .br returns in addition @@ -367,10 +367,10 @@ Documentation .br to repeat the excitations -`n_con_int `_ +`n_con_int `_ Number of integers to represent the connections between determinants -`h_matrix_all_dets `_ +`h_matrix_all_dets `_ H matrix on the basis of the slater deter;inants defined by psi_det diff --git a/src/Hartree_Fock/Fock_matrix.irp.f b/src/Hartree_Fock/Fock_matrix.irp.f index 42bed645..69aa7e26 100644 --- a/src/Hartree_Fock/Fock_matrix.irp.f +++ b/src/Hartree_Fock/Fock_matrix.irp.f @@ -217,3 +217,46 @@ BEGIN_PROVIDER [ double precision, HF_energy ] END_PROVIDER +BEGIN_PROVIDER [ double precision, Fock_matrix_ao, (ao_num_align, ao_num) ] + implicit none + BEGIN_DOC + ! Fock matrix in AO basis set + END_DOC + + if (elec_alpha_num == elec_beta_num) then + integer :: i,j + do j=1,ao_num + !DIR$ VECTOR ALIGNED + do i=1,ao_num_align + Fock_matrix_ao(i,j) = Fock_matrix_alpha_ao(i,j) + enddo + enddo + else + double precision, allocatable :: T(:,:), M(:,:) + ! F_ao = S C F_mo C^t S + allocate (T(mo_tot_num_align,mo_tot_num),M(ao_num_align,mo_tot_num)) + call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0, & + ao_overlap, size(ao_overlap,1), & + mo_coef, size(mo_coef,1), & + 0.d0, & + M, size(M,1)) + call dgemm('N','N', ao_num,mo_tot_num,mo_tot_num, 1.d0, & + M, size(M,1), & + Fock_matrix_mo, size(Fock_matrix_mo,1), & + 0.d0, & + T, size(T,1)) + call dgemm('N','T', mo_tot_num,ao_num,mo_tot_num, 1.d0, & + T, size(T,1), & + mo_coef, size(mo_coef,1), & + 0.d0, & + M, size(M,1)) + call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0, & + M, size(M,1), & + ao_overlap, size(ao_overlap,1), & + 0.d0, & + Fock_matrix_ao, size(Fock_matrix_ao,1)) + + deallocate(T) + endif +END_PROVIDER + diff --git a/src/Hartree_Fock/HF_density_matrix_ao.irp.f b/src/Hartree_Fock/HF_density_matrix_ao.irp.f index b0d58344..e85e4a6c 100644 --- a/src/Hartree_Fock/HF_density_matrix_ao.irp.f +++ b/src/Hartree_Fock/HF_density_matrix_ao.irp.f @@ -1,46 +1,27 @@ - BEGIN_PROVIDER [ double precision, HF_density_matrix_ao_alpha, (ao_num_align,ao_num) ] -&BEGIN_PROVIDER [ double precision, HF_density_matrix_ao_beta, (ao_num_align,ao_num) ] +BEGIN_PROVIDER [ double precision, HF_density_matrix_ao_alpha, (ao_num_align,ao_num) ] implicit none BEGIN_DOC - ! Alpha and Beta density matrix in the AO basis + ! Alpha density matrix in the AO basis END_DOC - integer :: i,j,k,l1,l2 - integer, allocatable :: mo_occ(:,:) - allocate ( mo_occ(elec_alpha_num,2) ) - call bitstring_to_list( HF_bitmask(1,1), mo_occ(1,1), j, N_int) - ASSERT ( j==elec_alpha_num ) + call dgemm('N','T',ao_num,ao_num,elec_alpha_num,1.d0, & + mo_coef, size(mo_coef,1), & + mo_coef, size(mo_coef,1), 0.d0, & + HF_density_matrix_ao_alpha, size(HF_density_matrix_ao_alpha,1)) + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, HF_density_matrix_ao_beta, (ao_num_align,ao_num) ] + implicit none + BEGIN_DOC + ! Beta density matrix in the AO basis + END_DOC - call bitstring_to_list( HF_bitmask(1,2), mo_occ(1,2), j, N_int) - ASSERT ( j==elec_beta_num ) - - do j=1,ao_num - !DIR$ VECTOR ALIGNED - do i=1,ao_num_align - HF_density_matrix_ao_alpha(i,j) = 0.d0 - HF_density_matrix_ao_beta (i,j) = 0.d0 - enddo - do k=1,elec_beta_num - l1 = mo_occ(k,1) - l2 = mo_occ(k,2) - !DIR$ VECTOR ALIGNED - do i=1,ao_num - HF_density_matrix_ao_alpha(i,j) = HF_density_matrix_ao_alpha(i,j) +& - mo_coef(i,l1) * mo_coef(j,l1) - HF_density_matrix_ao_beta (i,j) = HF_density_matrix_ao_beta (i,j) +& - mo_coef(i,l2) * mo_coef(j,l2) - enddo - enddo - do k=elec_beta_num+1,elec_alpha_num - l1 = mo_occ(k,1) - !DIR$ VECTOR ALIGNED - do i=1,ao_num - HF_density_matrix_ao_alpha(i,j) = HF_density_matrix_ao_alpha(i,j) +& - mo_coef(i,l1) * mo_coef(j,l1) - enddo - enddo - enddo - deallocate(mo_occ) + call dgemm('N','T',ao_num,ao_num,elec_beta_num,1.d0, & + mo_coef, size(mo_coef,1), & + mo_coef, size(mo_coef,1), 0.d0, & + HF_density_matrix_ao_beta, size(HF_density_matrix_ao_beta,1)) + END_PROVIDER BEGIN_PROVIDER [ double precision, HF_density_matrix_ao, (ao_num_align,ao_num) ] @@ -48,40 +29,13 @@ BEGIN_PROVIDER [ double precision, HF_density_matrix_ao, (ao_num_align,ao_num) ] BEGIN_DOC ! Density matrix in the AO basis END_DOC - integer :: i,j,k,l1,l2 - integer, allocatable :: mo_occ(:,:) + ASSERT (size(HF_density_matrix_ao,1) == size(HF_density_matrix_ao_alpha,1)) + if (elec_alpha_num== elec_beta_num) then + HF_density_matrix_ao = HF_density_matrix_ao_alpha + HF_density_matrix_ao_alpha + else + ASSERT (size(HF_density_matrix_ao,1) == size(HF_density_matrix_ao_beta ,1)) + HF_density_matrix_ao = HF_density_matrix_ao_alpha + HF_density_matrix_ao_beta + endif - allocate ( mo_occ(elec_alpha_num,2) ) - call bitstring_to_list( HF_bitmask(1,1), mo_occ(1,1), j, N_int) - ASSERT ( j==elec_alpha_num ) - - call bitstring_to_list( HF_bitmask(1,2), mo_occ(1,2), j, N_int) - ASSERT ( j==elec_beta_num ) - - do j=1,ao_num - !DIR$ VECTOR ALIGNED - do i=1,ao_num_align - HF_density_matrix_ao(i,j) = 0.d0 - enddo - do k=1,elec_beta_num - l1 = mo_occ(k,1) - l2 = mo_occ(k,2) - !DIR$ VECTOR ALIGNED - do i=1,ao_num - HF_density_matrix_ao(i,j) = HF_density_matrix_ao(i,j) + & - mo_coef(i,l1) * mo_coef(j,l1) + & - mo_coef(i,l2) * mo_coef(j,l2) - enddo - enddo - do k=elec_beta_num+1,elec_alpha_num - l1 = mo_occ(k,1) - !DIR$ VECTOR ALIGNED - do i=1,ao_num - HF_density_matrix_ao(i,j) = HF_density_matrix_ao(i,j) + & - mo_coef(i,l1) * mo_coef(j,l1) - enddo - enddo - enddo - deallocate(mo_occ) END_PROVIDER diff --git a/src/Hartree_Fock/README.rst b/src/Hartree_Fock/README.rst index 074928ad..7173d418 100644 --- a/src/Hartree_Fock/README.rst +++ b/src/Hartree_Fock/README.rst @@ -26,19 +26,22 @@ Documentation .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. -`fock_matrix_alpha_ao `_ +`fock_matrix_alpha_ao `_ Alpha Fock matrix in AO basis set -`fock_matrix_alpha_mo `_ +`fock_matrix_alpha_mo `_ Fock matrix on the MO basis -`fock_matrix_beta_ao `_ +`fock_matrix_ao `_ + Fock matrix in AO basis set + +`fock_matrix_beta_ao `_ Alpha Fock matrix in AO basis set -`fock_matrix_beta_mo `_ +`fock_matrix_beta_mo `_ Fock matrix on the MO basis -`fock_matrix_diag_mo `_ +`fock_matrix_diag_mo `_ Fock matrix on the MO basis. For open shells, the ROHF Fock Matrix is .br @@ -53,7 +56,7 @@ Documentation K = Fb - Fa .br -`fock_matrix_mo `_ +`fock_matrix_mo `_ Fock matrix on the MO basis. For open shells, the ROHF Fock Matrix is .br @@ -68,49 +71,49 @@ Documentation K = Fb - Fa .br -`hf_energy `_ +`hf_energy `_ Hartree-Fock energy -`hf_density_matrix_ao `_ +`hf_density_matrix_ao `_ Density matrix in the AO basis -`hf_density_matrix_ao_alpha `_ - Alpha and Beta density matrix in the AO basis +`hf_density_matrix_ao_alpha `_ + Alpha density matrix in the AO basis -`hf_density_matrix_ao_beta `_ - Alpha and Beta density matrix in the AO basis +`hf_density_matrix_ao_beta `_ + Beta density matrix in the AO basis -`diagonal_fock_matrix_mo `_ +`diagonal_fock_matrix_mo `_ Diagonal Fock matrix in the MO basis -`eigenvectors_fock_matrix_mo `_ +`eigenvectors_fock_matrix_mo `_ Diagonal Fock matrix in the MO basis -`scf_iteration `_ +`xcf_iteration `_ Undocumented -`do_diis `_ +`do_diis `_ If True, compute integrals on the fly -`n_it_scf_max `_ +`n_it_scf_max `_ Maximum number of SCF iterations -`thresh_scf `_ +`thresh_scf `_ Threshold on the convergence of the Hartree Fock energy -`bi_elec_ref_bitmask_energy `_ +`bi_elec_ref_bitmask_energy `_ Energy of the reference bitmask used in Slater rules -`kinetic_ref_bitmask_energy `_ +`kinetic_ref_bitmask_energy `_ Energy of the reference bitmask used in Slater rules -`mono_elec_ref_bitmask_energy `_ +`mono_elec_ref_bitmask_energy `_ Energy of the reference bitmask used in Slater rules -`nucl_elec_ref_bitmask_energy `_ +`nucl_elec_ref_bitmask_energy `_ Energy of the reference bitmask used in Slater rules -`ref_bitmask_energy `_ +`ref_bitmask_energy `_ Energy of the reference bitmask used in Slater rules diff --git a/src/Hartree_Fock/SCF.irp.f b/src/Hartree_Fock/SCF.irp.f new file mode 100644 index 00000000..cb6e8cbf --- /dev/null +++ b/src/Hartree_Fock/SCF.irp.f @@ -0,0 +1,95 @@ +BEGIN_PROVIDER [ integer, it_scf ] + implicit none + BEGIN_DOC + ! Number of the current SCF iteration + END_DOC + it_scf = 0 +END_PROVIDER + +BEGIN_PROVIDER [ double precision, SCF_density_matrices, (ao_num_align,ao_num,2,0:n_it_scf_max) ] + implicit none + BEGIN_DOC + ! Density matrices at every SCF iteration + END_DOC + SCF_density_matrices = 0.d0 +END_PROVIDER + +subroutine insert_new_SCF_density_matrix + implicit none + integer :: i,j + do j=1,ao_num + do i=1,ao_num + SCF_density_matrices(i,j,1,it_scf) = HF_density_matrix_ao_alpha(i,j) + SCF_density_matrices(i,j,1,0) += HF_density_matrix_ao_alpha(i,j) + SCF_density_matrices(i,j,2,it_scf) = HF_density_matrix_ao_beta(i,j) + SCF_density_matrices(i,j,2,0) += HF_density_matrix_ao_beta(i,j) + enddo + enddo +end + +subroutine Fock_mo_to_ao(FMO,LDFMO,FAO,LDFAO) + implicit none + integer, intent(in) :: LDFMO ! size(FMO,1) + integer, intent(in) :: LDFAO ! size(FAO,1) + double precision, intent(in) :: FMO(LDFMO,*) + double precision, intent(out) :: FAO(LDFAO,*) + + double precision, allocatable :: T(:,:), M(:,:) + ! F_ao = S C F_mo C^t S + allocate (T(mo_tot_num_align,mo_tot_num),M(ao_num_align,mo_tot_num)) + call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0, & + ao_overlap, size(ao_overlap,1), & + mo_coef, size(mo_coef,1), & + 0.d0, & + M, size(M,1)) + call dgemm('N','N', ao_num,mo_tot_num,mo_tot_num, 1.d0, & + M, size(M,1), & + FMO, size(FMO,1), & + 0.d0, & + T, size(T,1)) + call dgemm('N','T', mo_tot_num,ao_num,mo_tot_num, 1.d0, & + T, size(T,1), & + mo_coef, size(mo_coef,1), & + 0.d0, & + M, size(M,1)) + call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0, & + M, size(M,1), & + ao_overlap, size(ao_overlap,1), & + 0.d0, & + FAO, size(FAO,1)) + deallocate(T,M) +end + +subroutine DIIS_step + implicit none + integer :: i,j + double precision :: c + c = 1.d0/dble(it_scf) + do j=1,ao_num + do i=1,ao_num + HF_density_matrix_ao_alpha(i,j) = SCF_density_matrices(i,j,1,0) * c + HF_density_matrix_ao_beta (i,j) = SCF_density_matrices(i,j,2,0) * c + enddo + enddo + TOUCH HF_density_matrix_ao_alpha HF_density_matrix_ao_beta + +! call Fock_mo_to_ao(Fock_matrix_mo_alpha, size(Fock_matrix_mo_alpha,1), & +! Fock_matrix_alpha_ao, size(Fock_matrix_alpha_ao,1) ) +! call Fock_mo_to_ao(Fock_matrix_mo_beta, size(Fock_matrix_mo_beta,1), & +! Fock_matrix_beta_ao, size(Fock_matrix_beta_ao,1) ) +! SOFT_TOUCH Fock_matrix_alpha_ao Fock_matrix_beta_ao Fock_matrix_mo_alpha Fock_matrix_mo_beta +end + +subroutine scf_iteration + implicit none + integer :: i,j + do i=1,n_it_scf_max + it_scf += 1 + SOFT_TOUCH it_scf + mo_coef = eigenvectors_Fock_matrix_mo + TOUCH mo_coef + call insert_new_SCF_density_matrix + call DIIS_step + print *, HF_energy + enddo +end diff --git a/src/Hartree_Fock/diagonalize_fock.irp.f b/src/Hartree_Fock/diagonalize_fock.irp.f index c34a1415..c5faea9c 100644 --- a/src/Hartree_Fock/diagonalize_fock.irp.f +++ b/src/Hartree_Fock/diagonalize_fock.irp.f @@ -5,16 +5,52 @@ ! Diagonal Fock matrix in the MO basis END_DOC - double precision, allocatable :: R(:,:) - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: R + integer :: i,j + integer :: liwork, lwork, n, info + integer, allocatable :: iwork(:) + double precision, allocatable :: work(:), F(:,:), S(:,:) - allocate(R(mo_tot_num_align,mo_tot_num)) + allocate(F(ao_num_align,ao_num), S(ao_num_align,ao_num) ) + do j=1,ao_num + do i=1,ao_num + S(i,j) = ao_overlap(i,j) + F(i,j) = Fock_matrix_ao(i,j) + enddo + enddo + + n = ao_num + lwork = 1+6*n + 2*n*n + liwork = 3 + 5*n - call lapack_diag(diagonal_Fock_matrix_mo,R,Fock_matrix_mo,size(Fock_matrix_mo,1),& - mo_tot_num) - - call dgemm('N','N',ao_num,mo_tot_num,mo_tot_num,1.d0,mo_coef,size(mo_coef,1),& - R,size(R,1),0.d0,eigenvectors_Fock_matrix_mo,size(eigenvectors_Fock_matrix_mo,1)) - deallocate(R) + allocate(work(lwork), iwork(liwork) ) + + lwork = -1 + liwork = -1 + call dsygvd(1,'v','u',mo_tot_num,F,size(F,1),S,size(S,1),& + diagonal_Fock_matrix_mo, work, lwork, iwork, liwork, info) + + if (info /= 0) then + print *, irp_here//' failed' + stop 1 + endif + lwork = int(work(1)) + liwork = iwork(1) + deallocate(work,iwork) + allocate(work(lwork), iwork(liwork) ) + + call dsygvd(1,'v','u',mo_tot_num,F,size(F,1),S,size(S,1),& + diagonal_Fock_matrix_mo, work, lwork, iwork, liwork, info) + + if (info /= 0) then + print *, irp_here//' failed' + stop 1 + endif + do j=1,ao_num + do i=1,ao_num + eigenvectors_Fock_matrix_mo(i,j) = F(i,j) + enddo + enddo + + deallocate(work, iwork, F, S) END_PROVIDER diff --git a/src/Hartree_Fock/mo_SCF_iterations.irp.f b/src/Hartree_Fock/mo_SCF_iterations.irp.f index 01cab413..4854861f 100644 --- a/src/Hartree_Fock/mo_SCF_iterations.irp.f +++ b/src/Hartree_Fock/mo_SCF_iterations.irp.f @@ -1,4 +1,4 @@ -program scf_iteration +program xcf_iteration use bitmasks implicit none double precision :: SCF_energy_before,SCF_energy_after,diag_H_mat_elem,get_mo_bielec_integral @@ -32,6 +32,6 @@ program scf_iteration endif mo_label = "Canonical" TOUCH mo_label mo_coef - call save_mos +! call save_mos end diff --git a/src/MOGuess/README.rst b/src/MOGuess/README.rst index d8e72641..90785358 100644 --- a/src/MOGuess/README.rst +++ b/src/MOGuess/README.rst @@ -22,19 +22,19 @@ Documentation .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. -`h_core_guess `_ +`h_core_guess `_ Undocumented -`ao_ortho_lowdin_coef `_ +`ao_ortho_lowdin_coef `_ matrix of the coefficients of the mos generated by the orthonormalization by the S^{-1/2} canonical transformation of the aos ao_ortho_lowdin_coef(i,j) = coefficient of the ith ao on the jth ao_ortho_lowdin orbital -`ao_ortho_lowdin_overlap `_ +`ao_ortho_lowdin_overlap `_ overlap matrix of the ao_ortho_lowdin supposed to be the Identity -`ao_ortho_lowdin_nucl_elec_integral `_ +`ao_ortho_lowdin_nucl_elec_integral `_ Undocumented diff --git a/src/MOs/cholesky_mo.irp.f b/src/MOs/cholesky_mo.irp.f new file mode 100644 index 00000000..97b6abd2 --- /dev/null +++ b/src/MOs/cholesky_mo.irp.f @@ -0,0 +1,80 @@ +subroutine cholesky_mo(n,m,P,LDP,C,LDC,tol_in,rank) + implicit none + BEGIN_DOC +! Cholesky decomposition of AO Density matrix to +! generate MOs + END_DOC + integer, intent(in) :: n,m, LDC, LDP + double precision, intent(in) :: P(LDP,n) + double precision, intent(out) :: C(LDC,m) + double precision, intent(in) :: tol_in + integer, intent(out) :: rank + + integer :: info + integer :: i,k + integer :: ipiv(n) + double precision:: tol + double precision, allocatable :: W(:,:), work(:) + !DEC$ ATTRIBUTES ALIGN: 32 :: W + !DEC$ ATTRIBUTES ALIGN: 32 :: work + !DEC$ ATTRIBUTES ALIGN: 32 :: ipiv + + allocate(W(LDC,n),work(2*n)) + tol=tol_in + + info = 0 + do i=1,n + do k=1,i + W(i,k) = P(i,k) + enddo + do k=i+1,n + W(i,k) = 0. + enddo + enddo + call DPSTRF('L', n, W, LDC, ipiv, rank, tol, work, info ) + do i=1,n + do k=1,min(m,rank) + C(ipiv(i),k) = W(i,k) + enddo + enddo + + deallocate(W,work) +end + +BEGIN_PROVIDER [ double precision, mo_density_matrix, (mo_tot_num_align, mo_tot_num) ] + implicit none + BEGIN_DOC + ! Density matrix in MO basis + END_DOC + integer :: i,j,k + mo_density_matrix = 0.d0 + do k=1,mo_tot_num + if (mo_occ(k) == 0.d0) then + cycle + endif + do j=1,ao_num + do i=1,ao_num + mo_density_matrix(i,j) = mo_density_matrix(i,j) + & + mo_occ(k) * mo_coef(i,k) * mo_coef(j,k) + enddo + enddo + enddo +END_PROVIDER + +BEGIN_PROVIDER [ double precision, mo_density_matrix_virtual, (mo_tot_num_align, mo_tot_num) ] + implicit none + BEGIN_DOC + ! Density matrix in MO basis (virtual MOs) + END_DOC + integer :: i,j,k + mo_density_matrix_virtual = 0.d0 + do k=1,mo_tot_num + do j=1,ao_num + do i=1,ao_num + mo_density_matrix_virtual(i,j) = mo_density_matrix_virtual(i,j) + & + (2.d0-mo_occ(k)) * mo_coef(i,k) * mo_coef(j,k) + enddo + enddo + enddo +END_PROVIDER + diff --git a/src/MOs/mos.ezfio_config b/src/MOs/mos.ezfio_config index a0eda491..e34d1aab 100644 --- a/src/MOs/mos.ezfio_config +++ b/src/MOs/mos.ezfio_config @@ -1,5 +1,7 @@ mo_basis + ao_num integer mo_tot_num integer mo_coef double precision (ao_basis_ao_num,mo_basis_mo_tot_num) mo_label character*(64) + mo_occ double precision (mo_basis_mo_tot_num) diff --git a/src/MOs/mos.irp.f b/src/MOs/mos.irp.f index b081d8ce..7895ac55 100644 --- a/src/MOs/mos.irp.f +++ b/src/MOs/mos.irp.f @@ -75,3 +75,25 @@ BEGIN_PROVIDER [ double precision, mo_coef_transp, (mo_tot_num_align,ao_num) ] END_PROVIDER +BEGIN_PROVIDER [ double precision, mo_occ, (mo_tot_num) ] + implicit none + BEGIN_DOC + ! MO occupation numbers + END_DOC + PROVIDE ezfio_filename + logical :: exists + call ezfio_has_mo_basis_mo_occ(exists) + if (exists) then + call ezfio_get_mo_basis_mo_occ(mo_occ) + else + mo_occ = 0.d0 + integer :: i + do i=1,elec_beta_num + mo_occ(i) = 2.d0 + enddo + do i=elec_beta_num+1,elec_alpha_num + mo_occ(i) = 1.d0 + enddo + endif +END_PROVIDER + diff --git a/src/Perturbation/README.rst b/src/Perturbation/README.rst index 806ad292..aad1d0c0 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_moller_plesset `_ +`pt2_moller_plesset `_ compute the standard Moller-Plesset 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/(difference of orbital energies) .br -`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. @@ -102,7 +102,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. @@ -112,20 +112,46 @@ Documentation c_pert(i) = e_2_pert(i)/ .br -`fill_h_apply_buffer_selection `_ +`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 CI_electronic_energy(1) + .br + that could be repeated to this determinant. + .br + In addition, for the perturbative energetic contribution you have the standard second order + .br + e_2_pert = ^2/(Delta_E) + .br + and also the purely projected contribution + .br + H_pert_diag = c_pert + +`repeat_all_e_corr `_ + Undocumented + +`fill_h_apply_buffer_selection `_ Fill the H_apply buffer with determiants for the selection -`remove_small_contributions `_ +`remove_small_contributions `_ Remove determinants with small contributions. N_states is assumed to be provided. -`selection_criterion `_ +`selection_criterion `_ Threshold to select determinants. Set by selection routines. -`selection_criterion_factor `_ +`selection_criterion_factor `_ Threshold to select determinants. Set by selection routines. -`selection_criterion_min `_ +`selection_criterion_min `_ Threshold to select determinants. Set by selection routines.