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 b795223f..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 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/Hartree_Fock/README.rst b/src/Hartree_Fock/README.rst index a9d70bd4..341c6284 100644 --- a/src/Hartree_Fock/README.rst +++ b/src/Hartree_Fock/README.rst @@ -26,22 +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_ao `_ +`fock_matrix_ao `_ Fock matrix in AO basis set -`fock_matrix_beta_ao `_ +`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 @@ -56,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 @@ -71,70 +71,73 @@ 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 `_ +`hf_density_matrix_ao_alpha `_ Alpha density matrix in the AO basis -`hf_density_matrix_ao_beta `_ +`hf_density_matrix_ao_beta `_ Beta density matrix in the AO basis -`fock_mo_to_ao `_ +`fock_mo_to_ao `_ Undocumented -`insert_new_scf_density_matrix `_ +`insert_new_scf_density_matrix `_ Undocumented -`it_scf `_ +`it_scf `_ Number of the current SCF iteration -`scf_density_matrices `_ +`scf_density_matrices `_ Density matrices at every SCF iteration -`scf_energies `_ +`scf_energies `_ Density matrices at every SCF iteration -`scf_interpolation_step `_ +`scf_interpolation_step `_ Undocumented -`scf_iterations `_ +`scf_iterations `_ Undocumented -`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 -`xcf_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 `_ +`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 index 895ca546..650aa365 100644 --- a/src/Hartree_Fock/SCF.irp.f +++ b/src/Hartree_Fock/SCF.irp.f @@ -70,10 +70,11 @@ subroutine SCF_interpolation_step 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)+SCF_density_matrices(i,j,1,it_scf-1) * (1.d0 - c) - HF_density_matrix_ao_beta (i,j) = c*SCF_density_matrices(i,j,2,it_scf)+SCF_density_matrices(i,j,2,it_scf-1) * (1.d0 - c) + 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 diff --git a/src/Hartree_Fock/mo_SCF_iterations.irp.f b/src/Hartree_Fock/mo_SCF_iterations.irp.f index 2698e010..55c42178 100644 --- a/src/Hartree_Fock/mo_SCF_iterations.irp.f +++ b/src/Hartree_Fock/mo_SCF_iterations.irp.f @@ -1,4 +1,11 @@ -program xcf_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 @@ -8,7 +15,7 @@ program xcf_iteration E0 = HF_energy thresh_SCF = 1.d-10 - call scf_iterations + call damping_SCF mo_label = "Canonical" TOUCH mo_label mo_coef call save_mos diff --git a/src/MOs/utils.irp.f b/src/MOs/utils.irp.f index ee04ac4b..5e7984ff 100644 --- a/src/MOs/utils.irp.f +++ b/src/MOs/utils.irp.f @@ -28,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