diff --git a/scripts/generate_h_apply.py b/scripts/generate_h_apply.py index a9d1c966..a57f869f 100755 --- a/scripts/generate_h_apply.py +++ b/scripts/generate_h_apply.py @@ -27,7 +27,7 @@ filter_integrals class H_apply(object): - def __init__(self,sub,SingleRef=False): + def __init__(self,sub,SingleRef=False,do_mono_exc=True, do_double_exc=True): s = {} for k in keywords: s[k] = "" @@ -56,6 +56,9 @@ class H_apply(object): s["omp_do"] = "!$OMP DO SCHEDULE (static)" s["omp_enddo"] = "!$OMP ENDDO NOWAIT" + d = { True : '.True.', False : '.False.'} + s["do_mono_excitations"] = d[do_mono_exc] + s["do_double_excitations"] = d[do_double_exc] s["keys_work"] += "call fill_H_apply_buffer_no_selection(key_idx,keys_out,N_int,iproc)" s["filter_integrals"] = "array_pairs = .True." diff --git a/src/AOs/aos.irp.f b/src/AOs/aos.irp.f index 4b3a74c0..396d9aab 100644 --- a/src/AOs/aos.irp.f +++ b/src/AOs/aos.irp.f @@ -93,46 +93,6 @@ END_PROVIDER enddo END_PROVIDER -BEGIN_PROVIDER [ double precision, ao_overlap, (ao_num_align,ao_num) ] - implicit none - BEGIN_DOC -! matrix of the overlap for tha AOs -! S(i,j) = overlap between the ith and the jth atomic basis function - END_DOC - integer :: i,j,k,l,nz,num_i,num_j,powA(3),powB(3) - double precision :: accu,overlap_x,overlap_y,overlap_z,A_center(3),B_center(3),norm - nz=100 - do i = 1, ao_num - num_i = ao_nucl(i) - powA(1) = ao_power(i,1) - powA(2) = ao_power(i,2) - powA(3) = ao_power(i,3) - A_center(1)=nucl_coord(num_i,1) - A_center(2)=nucl_coord(num_i,2) - A_center(3)=nucl_coord(num_i,3) - do j = i, ao_num - num_j = ao_nucl(j) - powB(1) = ao_power(j,1) - powB(2) = ao_power(j,2) - powB(3) = ao_power(j,3) - B_center(1)=nucl_coord(num_j,1) - B_center(2)=nucl_coord(num_j,2) - B_center(3)=nucl_coord(num_j,3) - accu = 0.d0 - do k = 1, ao_prim_num(i) - do l = 1, ao_prim_num(j) - call overlap_gaussian_xyz(A_center,B_center,ao_expo(i,k),ao_expo(j,l),powA,powB,overlap_x,overlap_y,overlap_z,norm,nz) - accu = accu + norm * ao_coef(i,k) * ao_coef(j,l) - enddo - enddo - ao_overlap(i,j) = accu - ao_overlap(j,i) = accu - enddo - enddo - -END_PROVIDER - - BEGIN_PROVIDER [ double precision, ao_coef_transp, (ao_prim_num_max_align,ao_num) ] &BEGIN_PROVIDER [ double precision, ao_expo_transp, (ao_prim_num_max_align,ao_num) ] implicit none diff --git a/src/CID/H_apply.irp.f b/src/CID/H_apply.irp.f index 0df1da38..6a1bb536 100644 --- a/src/CID/H_apply.irp.f +++ b/src/CID/H_apply.irp.f @@ -3,7 +3,7 @@ BEGIN_SHELL [ /usr/bin/env python ] from generate_h_apply import H_apply -H = H_apply("cisd") +H = H_apply("cisd",do_double_exc=True,do_mono_exc=False) print H END_SHELL diff --git a/src/CID/README.rst b/src/CID/README.rst index 003b51fa..5a0c5026 100644 --- a/src/CID/README.rst +++ b/src/CID/README.rst @@ -36,7 +36,7 @@ Documentation .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. -`cisd `_ +`cisd `_ Undocumented diff --git a/src/CID/do_mono_double.irp.f b/src/CID/do_mono_double.irp.f deleted file mode 100644 index f211879e..00000000 --- a/src/CID/do_mono_double.irp.f +++ /dev/null @@ -1,19 +0,0 @@ -BEGIN_PROVIDER [logical, do_double_excitations] - implicit none - BEGIN_DOC - ! if True then the double excitations are performed in the calculation - ! always true in the CISD - END_DOC - do_double_excitations = .True. - -END_PROVIDER - -BEGIN_PROVIDER [logical, do_mono_excitations] - implicit none - BEGIN_DOC - ! if True then the mono excitations are performed in the calculation - ! always true in the CISD - END_DOC - do_mono_excitations = .False. - -END_PROVIDER diff --git a/src/CID_SC2_selected/H_apply.irp.f b/src/CID_SC2_selected/H_apply.irp.f index 79668af7..e284ab43 100644 --- a/src/CID_SC2_selected/H_apply.irp.f +++ b/src/CID_SC2_selected/H_apply.irp.f @@ -3,7 +3,7 @@ BEGIN_SHELL [ /usr/bin/env python ] from generate_h_apply import * from perturbation import perturbations -s = H_apply("PT2",SingleRef=True) +s = H_apply("PT2",SingleRef=True,do_mono_exc=False,do_double_exc=True) s.set_perturbation("epstein_nesbet_sc2_projected") print s END_SHELL diff --git a/src/CID_SC2_selected/README.rst b/src/CID_SC2_selected/README.rst index b6206850..8688fd41 100644 --- a/src/CID_SC2_selected/README.rst +++ b/src/CID_SC2_selected/README.rst @@ -8,7 +8,7 @@ Documentation .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. -`cisd_sc2_selected `_ +`cisd_sc2_selected `_ Undocumented diff --git a/src/CID_SC2_selected/do_mono_double.irp.f b/src/CID_SC2_selected/do_mono_double.irp.f deleted file mode 100644 index f211879e..00000000 --- a/src/CID_SC2_selected/do_mono_double.irp.f +++ /dev/null @@ -1,19 +0,0 @@ -BEGIN_PROVIDER [logical, do_double_excitations] - implicit none - BEGIN_DOC - ! if True then the double excitations are performed in the calculation - ! always true in the CISD - END_DOC - do_double_excitations = .True. - -END_PROVIDER - -BEGIN_PROVIDER [logical, do_mono_excitations] - implicit none - BEGIN_DOC - ! if True then the mono excitations are performed in the calculation - ! always true in the CISD - END_DOC - do_mono_excitations = .False. - -END_PROVIDER diff --git a/src/CID_selected/H_apply.irp.f b/src/CID_selected/H_apply.irp.f index 91dfb9fc..e3afaa9d 100644 --- a/src/CID_selected/H_apply.irp.f +++ b/src/CID_selected/H_apply.irp.f @@ -4,7 +4,7 @@ from generate_h_apply import * from perturbation import perturbations for perturbation in perturbations: - s = H_apply("cisd_selection_"+perturbation) + s = H_apply("cisd_selection_"+perturbation,do_mono_exc=False) s.set_selection_pt2(perturbation) print s END_SHELL diff --git a/src/CID_selected/README.rst b/src/CID_selected/README.rst index 4fed26e5..8c64bc25 100644 --- a/src/CID_selected/README.rst +++ b/src/CID_selected/README.rst @@ -8,10 +8,10 @@ Documentation .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. -`h_apply_cisd_selection `_ +`h_apply_cisd_selection `_ Undocumented -`cisd `_ +`cisd `_ Undocumented diff --git a/src/CID_selected/do_mono_double.irp.f b/src/CID_selected/do_mono_double.irp.f deleted file mode 100644 index f211879e..00000000 --- a/src/CID_selected/do_mono_double.irp.f +++ /dev/null @@ -1,19 +0,0 @@ -BEGIN_PROVIDER [logical, do_double_excitations] - implicit none - BEGIN_DOC - ! if True then the double excitations are performed in the calculation - ! always true in the CISD - END_DOC - do_double_excitations = .True. - -END_PROVIDER - -BEGIN_PROVIDER [logical, do_mono_excitations] - implicit none - BEGIN_DOC - ! if True then the mono excitations are performed in the calculation - ! always true in the CISD - END_DOC - do_mono_excitations = .False. - -END_PROVIDER diff --git a/src/CISD/do_mono_double.irp.f b/src/CISD/do_mono_double.irp.f deleted file mode 100644 index 8a07c292..00000000 --- a/src/CISD/do_mono_double.irp.f +++ /dev/null @@ -1,19 +0,0 @@ -BEGIN_PROVIDER [logical, do_double_excitations] - implicit none - BEGIN_DOC - ! if True then the double excitations are performed in the calculation - ! always true in the CISD - END_DOC - do_double_excitations = .True. - -END_PROVIDER - -BEGIN_PROVIDER [logical, do_mono_excitations] - implicit none - BEGIN_DOC - ! if True then the mono excitations are performed in the calculation - ! always true in the CISD - END_DOC - do_mono_excitations = .True. - -END_PROVIDER diff --git a/src/CISD_SC2_selected/do_mono_double.irp.f b/src/CISD_SC2_selected/do_mono_double.irp.f deleted file mode 100644 index 8a07c292..00000000 --- a/src/CISD_SC2_selected/do_mono_double.irp.f +++ /dev/null @@ -1,19 +0,0 @@ -BEGIN_PROVIDER [logical, do_double_excitations] - implicit none - BEGIN_DOC - ! if True then the double excitations are performed in the calculation - ! always true in the CISD - END_DOC - do_double_excitations = .True. - -END_PROVIDER - -BEGIN_PROVIDER [logical, do_mono_excitations] - implicit none - BEGIN_DOC - ! if True then the mono excitations are performed in the calculation - ! always true in the CISD - END_DOC - do_mono_excitations = .True. - -END_PROVIDER diff --git a/src/CISD_selected/do_mono_double.irp.f b/src/CISD_selected/do_mono_double.irp.f deleted file mode 100644 index 8a07c292..00000000 --- a/src/CISD_selected/do_mono_double.irp.f +++ /dev/null @@ -1,19 +0,0 @@ -BEGIN_PROVIDER [logical, do_double_excitations] - implicit none - BEGIN_DOC - ! if True then the double excitations are performed in the calculation - ! always true in the CISD - END_DOC - do_double_excitations = .True. - -END_PROVIDER - -BEGIN_PROVIDER [logical, do_mono_excitations] - implicit none - BEGIN_DOC - ! if True then the mono excitations are performed in the calculation - ! always true in the CISD - END_DOC - do_mono_excitations = .True. - -END_PROVIDER diff --git a/src/DensityMatrix/density_matrix.irp.f b/src/DensityMatrix/density_matrix.irp.f index 5afcb523..3fc1d6f9 100644 --- a/src/DensityMatrix/density_matrix.irp.f +++ b/src/DensityMatrix/density_matrix.irp.f @@ -209,51 +209,3 @@ END_PROVIDER enddo enddo END_PROVIDER - - BEGIN_PROVIDER [ double precision, one_body_dm_a, (mo_tot_num_align,mo_tot_num) ] -&BEGIN_PROVIDER [ double precision, one_body_dm_b, (mo_tot_num_align,mo_tot_num) ] - implicit none - BEGIN_DOC - ! Alpha and beta one-body density matrix - END_DOC - - integer :: j,k,l - integer :: occ(N_int*bit_kind_size,2) - double precision :: ck, cl, ckl - double precision :: phase - integer :: h1,h2,p1,p2,s1,s2, degree - integer :: exc(0:2,2,2),n_occ_alpha - one_body_dm_a = 0.d0 - one_body_dm_b = 0.d0 - - do k=1,det_num - call bitstring_to_list(det_provider(1,1,k), occ(1,1), n_occ_alpha, N_int) - call bitstring_to_list(det_provider(1,2,k), occ(1,2), n_occ_alpha, N_int) - ck = det_coef_provider(k) - do l=1,elec_alpha_num - j = occ(l,1) - one_body_dm_a(j,j) += ck*ck - enddo - do l=1,elec_beta_num - j = occ(l,2) - one_body_dm_b(j,j) += ck*ck - enddo - do l=1,k-1 - call get_excitation_degree(det_provider(1,1,k),det_provider(1,1,l),degree,N_int) - if (degree /= 1) then - cycle - endif - call get_mono_excitation(det_provider(1,1,k),det_provider(1,1,l),exc,phase,N_int) - call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) - ckl = ck * det_coef_provider(l) * phase - if (s1==1) then - one_body_dm_a(h1,p1) += ckl - one_body_dm_a(p1,h1) += ckl - else - one_body_dm_b(h1,p1) += ckl - one_body_dm_b(p1,h1) += ckl - endif - enddo - enddo -END_PROVIDER - diff --git a/src/Dets/H_apply_template.f b/src/Dets/H_apply_template.f index 773fbd92..1cb5d4e7 100644 --- a/src/Dets/H_apply_template.f +++ b/src/Dets/H_apply_template.f @@ -389,7 +389,7 @@ subroutine $subroutine($params_main) !$OMP PARALLEL DEFAULT(SHARED) & !$OMP PRIVATE(i_generator,wall_2,ispin,k,mask) allocate( mask(N_int,2,6) ) - !$OMP DO SCHEDULE(guided) + !$OMP DO SCHEDULE(dynamic,4) do i_generator=1,nmax if (abort_here) then cycle @@ -420,13 +420,13 @@ subroutine $subroutine($params_main) enddo enddo - if(do_double_excitations)then + if($do_double_excitations)then call $subroutine_diexc(psi_generators(1,1,i_generator), & mask(1,1,d_hole1), mask(1,1,d_part1), & mask(1,1,d_hole2), mask(1,1,d_part2), & i_generator $params_post) endif - if(do_mono_excitations)then + if($do_mono_excitations)then call $subroutine_monoexc(psi_generators(1,1,i_generator), & mask(1,1,s_hole ), mask(1,1,s_part ), & i_generator $params_post) @@ -475,13 +475,13 @@ subroutine $subroutine($params_main) not(psi_generators(k,ispin,i_generator)) ) enddo enddo - if(do_double_excitations)then + if($do_double_excitations)then call $subroutine_diexc(psi_generators(1,1,i_generator), & mask(1,1,d_hole1), mask(1,1,d_part1), & mask(1,1,d_hole2), mask(1,1,d_part2), & i_generator $params_post) endif - if(do_mono_excitations)then + if($do_mono_excitations)then call $subroutine_monoexc(psi_generators(1,1,i_generator), & mask(1,1,s_hole ), mask(1,1,s_part ), & i_generator $params_post) diff --git a/src/Dets/README.rst b/src/Dets/README.rst index 35c61523..4a4e8ff2 100644 --- a/src/Dets/README.rst +++ b/src/Dets/README.rst @@ -50,11 +50,11 @@ 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 `_ @@ -64,7 +64,7 @@ Documentation `h_apply_threshold `_ Theshold on | | -`resize_h_apply_buffer `_ +`resize_h_apply_buffer `_ Undocumented `cisd_sc2 `_ @@ -152,55 +152,55 @@ Documentation `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 `_ Number of determinants in the wave function -`n_det_max_jacobi `_ +`n_det_max_jacobi `_ Maximum number of determinants diagonalized my jacobi `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 `_ @@ -237,19 +237,22 @@ Documentation 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 `_ 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 +`threshold_convergence_sc2 `_ + convergence of the correlation energy of SC2 iterations + `filter_connected `_ Filters out the determinants that are not connected by H .br diff --git a/src/Dets/density_matrix.irp.f b/src/Dets/density_matrix.irp.f new file mode 100644 index 00000000..e5a907fb --- /dev/null +++ b/src/Dets/density_matrix.irp.f @@ -0,0 +1,106 @@ + BEGIN_PROVIDER [ double precision, one_body_dm_mo_alpha, (mo_tot_num_align,mo_tot_num) ] +&BEGIN_PROVIDER [ double precision, one_body_dm_mo_beta, (mo_tot_num_align,mo_tot_num) ] + implicit none + BEGIN_DOC + ! Alpha and beta one-body density matrix for each state + END_DOC + + integer :: j,k,l,m + integer :: occ(N_int*bit_kind_size,2) + double precision :: ck, cl, ckl + double precision :: phase + integer :: h1,h2,p1,p2,s1,s2, degree + integer :: exc(0:2,2,2),n_occ_alpha + double precision, allocatable :: tmp_a(:,:), tmp_b(:,:) + + one_body_dm_mo_alpha = 0.d0 + one_body_dm_mo_beta = 0.d0 + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(j,k,l,m,occ,ck, cl, ckl,phase,h1,h2,p1,p2,s1,s2, degree,exc, & + !$OMP tmp_a, tmp_b, n_occ_alpha)& + !$OMP SHARED(psi_det,psi_coef,N_int,N_states,state_average_weight,elec_alpha_num,& + !$OMP elec_beta_num,one_body_dm_mo_alpha,one_body_dm_mo_beta,N_det,mo_tot_num_align,& + !$OMP mo_tot_num) + allocate(tmp_a(mo_tot_num_align,mo_tot_num), tmp_b(mo_tot_num_align,mo_tot_num) ) + tmp_a = 0.d0 + tmp_b = 0.d0 + !$OMP DO SCHEDULE(dynamic) + do k=1,N_det + call bitstring_to_list(psi_det(1,1,k), occ(1,1), n_occ_alpha, N_int) + call bitstring_to_list(psi_det(1,2,k), occ(1,2), n_occ_alpha, N_int) + do m=1,N_states + ck = psi_coef(k,m)*psi_coef(k,m) * state_average_weight(m) + do l=1,elec_alpha_num + j = occ(l,1) + tmp_a(j,j) += ck + enddo + do l=1,elec_beta_num + j = occ(l,2) + tmp_b(j,j) += ck + enddo + enddo + do l=1,k-1 + call get_excitation_degree(psi_det(1,1,k),psi_det(1,1,l),degree,N_int) + if (degree /= 1) then + cycle + endif + call get_mono_excitation(psi_det(1,1,k),psi_det(1,1,l),exc,phase,N_int) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + do m=1,N_states + ckl = psi_coef(k,m) * psi_coef(l,m) * phase * state_average_weight(m) + if (s1==1) then + tmp_a(h1,p1) += ckl + tmp_a(p1,h1) += ckl + else + tmp_b(h1,p1) += ckl + tmp_b(p1,h1) += ckl + endif + enddo + enddo + enddo + !$OMP END DO NOWAIT + !$OMP CRITICAL + one_body_dm_mo_alpha = one_body_dm_mo_alpha + tmp_a + !$OMP END CRITICAL + !$OMP CRITICAL + one_body_dm_mo_beta = one_body_dm_mo_beta + tmp_b + !$OMP END CRITICAL + deallocate(tmp_a,tmp_b) + !$OMP BARRIER + !$OMP END PARALLEL +END_PROVIDER + +BEGIN_PROVIDER [ double precision, one_body_dm_mo, (mo_tot_num_align,mo_tot_num) ] + implicit none + BEGIN_DOC + ! One-body density matrix + END_DOC + one_body_dm_mo = one_body_dm_mo_alpha + one_body_dm_mo_beta +END_PROVIDER + +subroutine save_natural_mos + implicit none + BEGIN_DOC + ! Save natural orbitals, obtained by diagonalization of the one-body density matrix in the MO basis + END_DOC + character*(64) :: label + double precision, allocatable :: tmp(:,:) + allocate(tmp(size(one_body_dm_mo,1),size(one_body_dm_mo,2))) + + ! Negation to have the occupied MOs first after the diagonalization + tmp = -one_body_dm_mo + label = "Natural" + call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label) + deallocate(tmp) + call save_mos + +end + +BEGIN_PROVIDER [ double precision, state_average_weight, (N_states) ] + implicit none + BEGIN_DOC + ! Weights in the state-average calculation of the density matrix + END_DOC + state_average_weight = 1.d0/dble(N_states) +END_PROVIDER + diff --git a/src/Dets/save_natorb.irp.f b/src/Dets/save_natorb.irp.f new file mode 100644 index 00000000..dd3d7684 --- /dev/null +++ b/src/Dets/save_natorb.irp.f @@ -0,0 +1,4 @@ +program save_natorb + call save_natural_mos +end + diff --git a/src/Full_CI/full_ci.irp.f b/src/Full_CI/full_ci.irp.f index 532d3788..49a25b76 100644 --- a/src/Full_CI/full_ci.irp.f +++ b/src/Full_CI/full_ci.irp.f @@ -25,14 +25,14 @@ program cisd exit endif enddo - print *, 'Final step' - call remove_small_contributions - call diagonalize_CI - print *, 'N_det = ', N_det - print *, 'N_states = ', N_states - print *, 'PT2 = ', pt2 - print *, 'E = ', CI_energy - print *, 'E+PT2 = ', CI_energy+pt2 - print *, '-----' +! print *, 'Final step' +! call remove_small_contributions +! call diagonalize_CI +! print *, 'N_det = ', N_det +! print *, 'N_states = ', N_states +! print *, 'PT2 = ', pt2 +! print *, 'E = ', CI_energy +! print *, 'E+PT2 = ', CI_energy+pt2 +! print *, '-----' deallocate(pt2,norm_pert) end diff --git a/src/Generators_full/generators.irp.f b/src/Generators_full/generators.irp.f index 3f3fba86..9d8761d2 100644 --- a/src/Generators_full/generators.irp.f +++ b/src/Generators_full/generators.irp.f @@ -47,7 +47,13 @@ BEGIN_PROVIDER [ integer(bit_kind), psi_generators, (N_int,2,psi_det_size) ] ! For Single reference wave functions, the generator is the ! Hartree-Fock determinant END_DOC - psi_generators = psi_det_sorted + integer :: i, k + do i=1,N_det + do k=1,N_int + psi_generators(k,1,i) = psi_det_sorted(k,1,i) + psi_generators(k,2,i) = psi_det_sorted(k,2,i) + enddo + enddo END_PROVIDER 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..341c6284 100644 --- a/src/Hartree_Fock/README.rst +++ b/src/Hartree_Fock/README.rst @@ -32,6 +32,9 @@ Documentation `fock_matrix_alpha_mo `_ Fock matrix on the MO basis +`fock_matrix_ao `_ + Fock matrix in AO basis set + `fock_matrix_beta_ao `_ Alpha Fock matrix in AO basis set @@ -71,14 +74,35 @@ Documentation `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 + 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 + +`fock_mo_to_ao `_ + Undocumented + +`insert_new_scf_density_matrix `_ + Undocumented + +`it_scf `_ + Number of the current SCF iteration + +`scf_density_matrices `_ + Density matrices at every SCF iteration + +`scf_energies `_ + Density matrices at every SCF iteration + +`scf_interpolation_step `_ + Undocumented + +`scf_iterations `_ + Undocumented `diagonal_fock_matrix_mo `_ Diagonal Fock matrix in the MO basis @@ -86,13 +110,16 @@ Documentation `eigenvectors_fock_matrix_mo `_ Diagonal Fock matrix in the MO basis -`scf_iteration `_ +`run `_ Undocumented -`do_diis `_ +`scf `_ + Undocumented + +`do_diis `_ If True, compute integrals on the fly -`n_it_scf_max `_ +`n_it_scf_max `_ Maximum number of SCF iterations `thresh_scf `_ diff --git a/src/Hartree_Fock/SCF.irp.f b/src/Hartree_Fock/SCF.irp.f new file mode 100644 index 00000000..650aa365 --- /dev/null +++ b/src/Hartree_Fock/SCF.irp.f @@ -0,0 +1,108 @@ +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,n_it_scf_max) ] +&BEGIN_PROVIDER [ double precision, SCF_energies, (n_it_scf_max) ] + implicit none + BEGIN_DOC + ! Density matrices at every SCF iteration + END_DOC + SCF_density_matrices = 0.d0 + SCF_energies = 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,2,it_scf) = HF_density_matrix_ao_beta(i,j) + enddo + enddo + SCF_energies(it_scf) = HF_energy +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 SCF_interpolation_step + implicit none + integer :: i,j + double precision :: c + + if (it_scf == 1) then + return + endif + call random_number(c) + c = c*0.5d0 + do j=1,ao_num + do i=1,ao_num + HF_density_matrix_ao_alpha(i,j) = c*SCF_density_matrices(i,j,1,it_scf-1)+SCF_density_matrices(i,j,1,it_scf) * (1.d0 - c) + HF_density_matrix_ao_beta (i,j) = c*SCF_density_matrices(i,j,2,it_scf-1)+SCF_density_matrices(i,j,2,it_scf) * (1.d0 - 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_iterations + 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 + print *, HF_energy + if (SCF_energies(it_scf)>SCF_energies(it_scf-1)) then + call SCF_interpolation_step + endif + if (it_scf>1 ) then + if (dabs(SCF_energies(it_scf)-SCF_energies(it_scf-1)) < thresh_SCF) then + exit + endif + endif + 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..55c42178 100644 --- a/src/Hartree_Fock/mo_SCF_iterations.irp.f +++ b/src/Hartree_Fock/mo_SCF_iterations.irp.f @@ -1,35 +1,21 @@ -program scf_iteration + +program scf + call orthonormalize_mos + call run +end + +subroutine run + use bitmasks implicit none double precision :: SCF_energy_before,SCF_energy_after,diag_H_mat_elem,get_mo_bielec_integral double precision :: E0 - integer :: i_it - + integer :: i_it, i, j, k + E0 = HF_energy - i_it = 0 - n_it_scf_max = 10 - SCF_energy_before = huge(1.d0) - SCF_energy_after = E0 - print *, E0 - mo_label = "Canonical" - thresh_SCF = 1.d-10 - do while (i_it < 40 .and. dabs(SCF_energy_before - SCF_energy_after) > thresh_SCF) - SCF_energy_before = SCF_energy_after - mo_coef = eigenvectors_Fock_matrix_mo - TOUCH mo_coef mo_label - call clear_mo_map - SCF_energy_after = HF_energy - print*,SCF_energy_after, dabs(SCF_energy_before - SCF_energy_after) - i_it +=1 - if(i_it > n_it_scf_max)exit - enddo - if (i_it >= n_it_scf_max) then - stop 'Failed' - endif - if (SCF_energy_after - E0 > thresh_SCF) then - stop 'Failed' - endif + thresh_SCF = 1.d-10 + call damping_SCF mo_label = "Canonical" TOUCH mo_label mo_coef call save_mos 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/NEEDED_MODULES b/src/MOs/NEEDED_MODULES index 6ec1892c..9bbf1b74 100644 --- a/src/MOs/NEEDED_MODULES +++ b/src/MOs/NEEDED_MODULES @@ -1 +1 @@ -AOs Ezfio_files Nuclei Output Utils +AOs Ezfio_files Nuclei Output Utils Electrons diff --git a/src/MOs/README.rst b/src/MOs/README.rst index 659db5d5..f27425bc 100644 --- a/src/MOs/README.rst +++ b/src/MOs/README.rst @@ -27,6 +27,7 @@ Needed Modules * `Nuclei `_ * `Output `_ * `Utils `_ +* `Electrons `_ Documentation ============= @@ -34,12 +35,19 @@ Documentation .. Do not edit this section. It was auto-generated from the .. NEEDED_MODULES file. +`cholesky_mo `_ + Cholesky decomposition of MO Density matrix to + generate MOs + +`mo_density_matrix `_ + Density matrix in MO basis + `mo_coef `_ Molecular orbital coefficients on AO basis set mo_coef(i,j) = coefficient of the ith ao on the jth mo mo_label : Label characterizing the MOS (local, canonical, natural, etc) -`mo_coef_transp `_ +`mo_coef_transp `_ Molecular orbital coefficients on AO basis set `mo_label `_ @@ -47,13 +55,16 @@ Documentation mo_coef(i,j) = coefficient of the ith ao on the jth mo mo_label : Label characterizing the MOS (local, canonical, natural, etc) +`mo_occ `_ + MO occupation numbers + `mo_tot_num `_ Total number of molecular orbitals and the size of the keys corresponding `mo_tot_num_align `_ Aligned variable for dimensioning of arrays -`mo_as_eigvectors_of_mo_matrix `_ +`mo_as_eigvectors_of_mo_matrix `_ Undocumented `save_mos `_ 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..738d2c2a 100644 --- a/src/MOs/mos.irp.f +++ b/src/MOs/mos.irp.f @@ -32,7 +32,7 @@ END_PROVIDER !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: buffer logical :: exists PROVIDE ezfio_filename - + !Label call ezfio_has_mo_basis_mo_label(exists) if (exists) then @@ -40,7 +40,7 @@ END_PROVIDER else mo_label = 'no_label' endif - + ! Coefs allocate(buffer(ao_num,mo_tot_num)) buffer = 0.d0 @@ -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/MOs/utils.irp.f b/src/MOs/utils.irp.f index 186f5efc..5e7984ff 100644 --- a/src/MOs/utils.irp.f +++ b/src/MOs/utils.irp.f @@ -14,6 +14,7 @@ subroutine save_mos enddo enddo call ezfio_set_mo_basis_mo_coef(buffer) + call ezfio_set_mo_basis_mo_occ(mo_occ) deallocate (buffer) end @@ -27,21 +28,31 @@ subroutine mo_as_eigvectors_of_mo_matrix(matrix,n,m,label) double precision, allocatable :: mo_coef_new(:,:), R(:,:),eigvalues(:) !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: mo_coef_new, R + call write_time(output_mos) if (m /= mo_tot_num) then print *, irp_here, ': Error : m/= mo_tot_num' + stop 1 endif - allocate(R(n,m)) - allocate(mo_coef_new(ao_num_align,m),eigvalues(m)) + allocate(R(n,m),mo_coef_new(ao_num_align,m),eigvalues(m)) mo_coef_new = mo_coef call lapack_diag(eigvalues,R,matrix,size(matrix,1),size(matrix,2)) integer :: i + write (output_mos,'(A)'), 'MOs are now **'//trim(label)//'**' + write (output_mos,'(A)'), '' + write (output_mos,'(A)'), 'Eigenvalues' + write (output_mos,'(A)'), '-----------' + write (output_mos,'(A)'), '' + write (output_mos,'(A)'), '======== ================' do i = 1, m - print*,'eigvalues(i) = ',eigvalues(i) + write (output_mos,'(I8,X,F16.10)'), i,eigvalues(i) enddo + write (output_mos,'(A)'), '======== ================' + write (output_mos,'(A)'), '' call dgemm('N','N',ao_num,m,m,1.d0,mo_coef_new,size(mo_coef_new,1),R,size(R,1),0.d0,mo_coef,size(mo_coef,1)) deallocate(mo_coef_new,R,eigvalues) + call write_time(output_mos) mo_label = label SOFT_TOUCH mo_coef mo_label diff --git a/src/MonoInts/NEEDED_MODULES b/src/MonoInts/NEEDED_MODULES index 56587081..190f8c6e 100644 --- a/src/MonoInts/NEEDED_MODULES +++ b/src/MonoInts/NEEDED_MODULES @@ -1 +1 @@ -AOs Ezfio_files MOs Nuclei Output Utils +AOs Electrons Ezfio_files MOs Nuclei Output Utils diff --git a/src/MonoInts/README.rst b/src/MonoInts/README.rst index a9bd07bc..d10cee82 100644 --- a/src/MonoInts/README.rst +++ b/src/MonoInts/README.rst @@ -5,6 +5,7 @@ Needed Modules .. NEEDED_MODULES file. * `AOs `_ +* `Electrons `_ * `Ezfio_files `_ * `MOs `_ * `Nuclei `_ 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.