From a0d58690546d906304978f83e8e360d66690e7d6 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Sat, 16 Jul 2016 16:09:50 +0200 Subject: [PATCH] Modifs of fobo-scf --- config/ifort.cfg | 2 +- plugins/All_singles/H_apply.irp.f | 7 + plugins/All_singles/all_1h_1p.irp.f | 2 +- plugins/All_singles/all_1h_1p_singles.irp.f | 76 + plugins/DFT_Utils/grid_density.irp.f | 31 +- plugins/DFT_Utils/integration_3d.irp.f | 34 +- plugins/DFT_Utils/integration_radial.irp.f | 3 +- plugins/FOBOCI/SC2_1h1p.irp.f | 699 +++++ plugins/FOBOCI/fobo_scf.irp.f | 2 +- .../foboci_lmct_mlct_threshold_old.irp.f | 4 + plugins/Properties/delta_rho.irp.f | 6 +- plugins/Properties/mulliken.irp.f | 15 +- plugins/Properties/print_spin_density.irp.f | 36 + plugins/Properties/provide_deltarho.irp.f | 11 + plugins/Properties/test_two_body_dm.irp.f | 46 +- plugins/loc_cele/loc_cele.irp.f | 309 +-- src/Bitmask/bitmasks.irp.f | 84 +- src/Determinants/density_matrix.irp.f | 12 +- ...gonalize_restart_and_save_all_states.irp.f | 16 + ...nalize_restart_and_save_lowest_state.irp.f | 25 + ...gonalize_restart_and_save_one_states.irp.f | 26 + ...gonalize_restart_and_save_two_states.irp.f | 27 + src/Determinants/occ_pattern.irp.f | 2 +- src/Determinants/print_H_matrix_restart.irp.f | 179 ++ src/Determinants/print_bitmask.irp.f | 11 + src/Determinants/print_holes_particles.irp.f | 36 + src/Determinants/print_wf.irp.f | 71 + src/Determinants/save_only_singles.irp.f | 50 + src/Determinants/test_3d.irp.f | 15 + src/Determinants/test_two_body.irp.f | 18 + src/Determinants/truncate_wf.irp.f | 23 +- src/Determinants/two_body_dm_map.irp.f | 369 ++- src/Integrals_Bielec/EZFIO.cfg | 7 + src/Integrals_Bielec/mo_bi_integrals.irp.f | 40 +- src/MO_Basis/mo_permutation.irp.f | 20 + src/MO_Basis/print_aos.irp.f | 53 + src/MO_Basis/print_mo_in_space.irp.f | 50 + src/Nuclei/atomic_radii.irp.f | 112 + src/Utils/angular_integration.irp.f | 2264 +++++++++++++++++ 39 files changed, 4424 insertions(+), 369 deletions(-) create mode 100644 plugins/All_singles/all_1h_1p_singles.irp.f create mode 100644 plugins/FOBOCI/SC2_1h1p.irp.f create mode 100644 plugins/Properties/print_spin_density.irp.f create mode 100644 plugins/Properties/provide_deltarho.irp.f create mode 100644 src/Determinants/diagonalize_restart_and_save_all_states.irp.f create mode 100644 src/Determinants/diagonalize_restart_and_save_lowest_state.irp.f create mode 100644 src/Determinants/diagonalize_restart_and_save_one_states.irp.f create mode 100644 src/Determinants/diagonalize_restart_and_save_two_states.irp.f create mode 100644 src/Determinants/print_H_matrix_restart.irp.f create mode 100644 src/Determinants/print_bitmask.irp.f create mode 100644 src/Determinants/print_holes_particles.irp.f create mode 100644 src/Determinants/print_wf.irp.f create mode 100644 src/Determinants/save_only_singles.irp.f create mode 100644 src/Determinants/test_3d.irp.f create mode 100644 src/Determinants/test_two_body.irp.f create mode 100644 src/MO_Basis/mo_permutation.irp.f create mode 100644 src/MO_Basis/print_aos.irp.f create mode 100644 src/MO_Basis/print_mo_in_space.irp.f create mode 100644 src/Nuclei/atomic_radii.irp.f create mode 100644 src/Utils/angular_integration.irp.f diff --git a/config/ifort.cfg b/config/ifort.cfg index a738a83c..da414912 100644 --- a/config/ifort.cfg +++ b/config/ifort.cfg @@ -18,7 +18,7 @@ IRPF90_FLAGS : --ninja --align=32 # 0 : Deactivate # [OPTION] -MODE : DEBUG ; [ OPT | PROFILE | DEBUG ] : Chooses the section below +MODE : OPT ; [ OPT | PROFILE | DEBUG ] : Chooses the section below CACHE : 1 ; Enable cache_compile.py OPENMP : 1 ; Append OpenMP flags diff --git a/plugins/All_singles/H_apply.irp.f b/plugins/All_singles/H_apply.irp.f index f34f003c..cb0976af 100644 --- a/plugins/All_singles/H_apply.irp.f +++ b/plugins/All_singles/H_apply.irp.f @@ -8,6 +8,13 @@ s.unset_skip() s.filter_only_1h1p() print s +s = H_apply("just_1h_1p_singles",do_double_exc=False) +s.set_selection_pt2("epstein_nesbet_2x2") +s.unset_skip() +s.filter_only_1h1p() +print s + + s = H_apply("just_mono",do_double_exc=False) s.set_selection_pt2("epstein_nesbet_2x2") s.unset_skip() diff --git a/plugins/All_singles/all_1h_1p.irp.f b/plugins/All_singles/all_1h_1p.irp.f index a2786248..7a3700b1 100644 --- a/plugins/All_singles/all_1h_1p.irp.f +++ b/plugins/All_singles/all_1h_1p.irp.f @@ -49,7 +49,7 @@ subroutine routine endif call save_wavefunction if(n_det_before == N_det)then - selection_criterion = selection_criterion * 0.5d0 + selection_criterion_factor = selection_criterion_factor * 0.5d0 endif enddo diff --git a/plugins/All_singles/all_1h_1p_singles.irp.f b/plugins/All_singles/all_1h_1p_singles.irp.f new file mode 100644 index 00000000..b76a14b3 --- /dev/null +++ b/plugins/All_singles/all_1h_1p_singles.irp.f @@ -0,0 +1,76 @@ +program restart_more_singles + BEGIN_DOC + ! Generates and select single and double excitations of type 1h-1p + ! on the top of a given restart wave function of type CAS + END_DOC + read_wf = .true. + touch read_wf + print*,'ref_bitmask_energy = ',ref_bitmask_energy + call routine + +end +subroutine routine + implicit none + integer :: i,k + double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:),E_before(:) + integer :: N_st, degree + integer :: n_det_before + N_st = N_states + allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st),E_before(N_st)) + i = 0 + print*,'N_det = ',N_det + print*,'n_det_max = ',n_det_max + print*,'pt2_max = ',pt2_max + pt2=-1.d0 + E_before = ref_bitmask_energy + do while (N_det < n_det_max.and.maxval(abs(pt2(1:N_st))) > pt2_max) + n_det_before = N_det + i += 1 + print*,'-----------------------' + print*,'i = ',i + call H_apply_just_1h_1p_singles(pt2, norm_pert, H_pert_diag, N_st) + call diagonalize_CI + print*,'N_det = ',N_det + print*,'E = ',CI_energy(1) + print*,'pt2 = ',pt2(1) + print*,'E+PT2 = ',E_before + pt2(1) + E_before = CI_energy + if(N_states_diag.gt.1)then + print*,'Variational Energy difference' + do i = 2, N_st + print*,'Delta E = ',CI_energy(i) - CI_energy(1) + enddo + endif + if(N_states.gt.1)then + print*,'Variational + perturbative Energy difference' + do i = 2, N_st + print*,'Delta E = ',E_before(i)+ pt2(i) - (E_before(1) + pt2(1)) + enddo + endif + call save_wavefunction + if(n_det_before == N_det)then + selection_criterion_factor = selection_criterion_factor * 0.5d0 + endif + + enddo + + threshold_davidson = 1.d-10 + soft_touch threshold_davidson davidson_criterion + call diagonalize_CI + if(N_states_diag.gt.1)then + print*,'Variational Energy difference' + do i = 2, N_st + print*,'Delta E = ',CI_energy(i) - CI_energy(1) + enddo + endif + if(N_states.gt.1)then + print*,'Variational + perturbative Energy difference' + do i = 2, N_st + print*,'Delta E = ',CI_energy(i)+ pt2(i) - (CI_energy(1) + pt2(1)) + enddo + endif + call ezfio_set_all_singles_energy(CI_energy) + + call save_wavefunction + deallocate(pt2,norm_pert) +end diff --git a/plugins/DFT_Utils/grid_density.irp.f b/plugins/DFT_Utils/grid_density.irp.f index ac3b702f..6071a18b 100644 --- a/plugins/DFT_Utils/grid_density.irp.f +++ b/plugins/DFT_Utils/grid_density.irp.f @@ -1,11 +1,11 @@ BEGIN_PROVIDER [integer, n_points_angular_grid] implicit none - n_points_angular_grid = 51 + n_points_angular_grid = 50 END_PROVIDER BEGIN_PROVIDER [integer, n_points_radial_grid] implicit none - n_points_radial_grid = 1000 + n_points_radial_grid = 10000 END_PROVIDER @@ -15,8 +15,28 @@ END_PROVIDER BEGIN_DOC ! weights and grid points for the integration on the angular variables on ! the unit sphere centered on (0,0,0) +! According to the LEBEDEV scheme END_DOC call cal_quad(n_points_angular_grid, angular_quadrature_points,weights_angular_points) + include 'constants.include.F' + integer :: i + double precision :: accu + double precision :: degre_rad +!degre_rad = 180.d0/pi +!accu = 0.d0 +!do i = 1, n_points_integration_angular_lebedev +! accu += weights_angular_integration_lebedev(i) +! weights_angular_points(i) = weights_angular_integration_lebedev(i) * 2.d0 * pi +! angular_quadrature_points(i,1) = dcos ( degre_rad * theta_angular_integration_lebedev(i)) & +! * dsin ( degre_rad * phi_angular_integration_lebedev(i)) +! angular_quadrature_points(i,2) = dsin ( degre_rad * theta_angular_integration_lebedev(i)) & +! * dsin ( degre_rad * phi_angular_integration_lebedev(i)) +! angular_quadrature_points(i,3) = dcos ( degre_rad * phi_angular_integration_lebedev(i)) +!enddo +!print*,'ANGULAR' +!print*,'' +!print*,'accu = ',accu +!ASSERT( dabs(accu - 1.D0) < 1.d-10) END_PROVIDER @@ -55,7 +75,7 @@ BEGIN_PROVIDER [double precision, grid_points_per_atom, (3,n_points_angular_grid x_ref = nucl_coord(i,1) y_ref = nucl_coord(i,2) z_ref = nucl_coord(i,3) - do j = 1, n_points_radial_grid + do j = 1, n_points_radial_grid-1 double precision :: x,r x = grid_points_radial(j) ! x value for the mapping of the [0, +\infty] to [0,1] r = knowles_function(alpha_knowles(int(nucl_charge(i))),m_knowles,x) ! value of the radial coordinate for the integration @@ -81,7 +101,7 @@ BEGIN_PROVIDER [double precision, weight_functions_at_grid_points, (n_points_ang double precision :: tmp_array(nucl_num) ! run over all points in space do j = 1, nucl_num ! that are referred to each atom - do k = 1, n_points_radial_grid !for each radial grid attached to the "jth" atom + do k = 1, n_points_radial_grid -1 !for each radial grid attached to the "jth" atom do l = 1, n_points_angular_grid ! for each angular point attached to the "jth" atom r(1) = grid_points_per_atom(1,l,k,j) r(2) = grid_points_per_atom(2,l,k,j) @@ -95,6 +115,7 @@ BEGIN_PROVIDER [double precision, weight_functions_at_grid_points, (n_points_ang enddo accu = 1.d0/accu weight_functions_at_grid_points(l,k,j) = tmp_array(j) * accu +! print*,weight_functions_at_grid_points(l,k,j) enddo enddo enddo @@ -110,7 +131,7 @@ END_PROVIDER double precision :: r(3) double precision :: aos_array(ao_num),mos_array(mo_tot_num) do j = 1, nucl_num - do k = 1, n_points_radial_grid + do k = 1, n_points_radial_grid -1 do l = 1, n_points_angular_grid one_body_dm_mo_alpha_at_grid_points(l,k,j) = 0.d0 one_body_dm_mo_beta_at_grid_points(l,k,j) = 0.d0 diff --git a/plugins/DFT_Utils/integration_3d.irp.f b/plugins/DFT_Utils/integration_3d.irp.f index e2c6d8d0..43eb1ab8 100644 --- a/plugins/DFT_Utils/integration_3d.irp.f +++ b/plugins/DFT_Utils/integration_3d.irp.f @@ -4,12 +4,18 @@ double precision function step_function_becke(x) double precision :: f_function_becke integer :: i,n_max_becke - step_function_becke = f_function_becke(x) - n_max_becke = 3 - do i = 1, n_max_becke - step_function_becke = f_function_becke(step_function_becke) - enddo - step_function_becke = 0.5d0*(1.d0 - step_function_becke) +!if(x.lt.-1.d0)then +! step_function_becke = 0.d0 +!else if (x .gt.1)then +! step_function_becke = 0.d0 +!else + step_function_becke = f_function_becke(x) +!!n_max_becke = 1 + do i = 1, 4 + step_function_becke = f_function_becke(step_function_becke) + enddo + step_function_becke = 0.5d0*(1.d0 - step_function_becke) +!endif end double precision function f_function_becke(x) @@ -46,19 +52,3 @@ double precision function cell_function_becke(r,atom_number) enddo end -double precision function weight_function_becke(r,atom_number) - implicit none - double precision, intent(in) :: r(3) - integer, intent(in) :: atom_number - BEGIN_DOC - ! atom_number :: atom on which the weight function of Becke (1988, JCP,88(4)) - ! r(1:3) :: x,y,z coordinantes of the current point - END_DOC - double precision :: cell_function_becke,accu - integer :: j - accu = 0.d0 - do j = 1, nucl_num - accu += cell_function_becke(r,j) - enddo - weight_function_becke = cell_function_becke(r,atom_number)/accu -end diff --git a/plugins/DFT_Utils/integration_radial.irp.f b/plugins/DFT_Utils/integration_radial.irp.f index 3670db14..4943783b 100644 --- a/plugins/DFT_Utils/integration_radial.irp.f +++ b/plugins/DFT_Utils/integration_radial.irp.f @@ -16,7 +16,7 @@ do j = 1, nucl_num integral_density_alpha_knowles_becke_per_atom(j) = 0.d0 integral_density_beta_knowles_becke_per_atom(j) = 0.d0 - do i = 1, n_points_radial_grid + do i = 1, n_points_radial_grid-1 ! Angular integration over the solid angle Omega for a FIXED angular coordinate "r" f_average_angular_alpha = 0.d0 f_average_angular_beta = 0.d0 @@ -48,7 +48,6 @@ END_PROVIDER double precision, intent(in) :: alpha,x integer, intent(in) :: m knowles_function = -alpha * dlog(1.d0-x**m) -!knowles_function = 1.d0 end double precision function derivative_knowles_function(alpha,m,x) diff --git a/plugins/FOBOCI/SC2_1h1p.irp.f b/plugins/FOBOCI/SC2_1h1p.irp.f new file mode 100644 index 00000000..d347c6e5 --- /dev/null +++ b/plugins/FOBOCI/SC2_1h1p.irp.f @@ -0,0 +1,699 @@ +subroutine dressing_1h1p(dets_in,u_in,diag_H_elements,dim_in,sze,N_st,Nint,convergence) + use bitmasks + implicit none + BEGIN_DOC + ! CISD+SC2 method :: take off all the disconnected terms of a ROHF+1h1p (selected or not) + ! + ! dets_in : bitmasks corresponding to determinants + ! + ! u_in : guess coefficients on the various states. Overwritten + ! on exit + ! + ! dim_in : leftmost dimension of u_in + ! + ! sze : Number of determinants + ! + ! N_st : Number of eigenstates + ! + ! Initial guess vectors are not necessarily orthonormal + END_DOC + integer, intent(in) :: dim_in, sze, N_st, Nint + integer(bit_kind), intent(in) :: dets_in(Nint,2,sze) + double precision, intent(inout) :: u_in(dim_in,N_st) + double precision, intent(out) :: diag_H_elements(dim_in) + double precision, intent(in) :: convergence + + integer :: i,j,k,l + integer :: n_singles + integer :: index_singles(sze),hole_particles_singles(sze,3) + integer :: n_doubles + integer :: index_doubles(sze),hole_particles_doubles(sze,2) + integer :: index_hf + double precision :: e_corr_singles(mo_tot_num,2) + double precision :: e_corr_doubles(mo_tot_num) + double precision :: e_corr_singles_total(2) + double precision :: e_corr_doubles_1h1p + + integer :: exc(0:2,2,2),degree + integer :: h1,h2,p1,p2,s1,s2 + integer :: other_spin(2) + double precision :: phase + integer(bit_kind) :: key_tmp(N_int,2) + integer :: i_ok + double precision :: phase_single_double,phase_double_hf,get_mo_bielec_integral_schwartz + double precision :: hij,c_ref,contrib + integer :: iorb + + other_spin(1) = 2 + other_spin(2) = 1 + + n_singles = 0 + n_doubles = 0 + do i = 1,sze + call get_excitation(ref_bitmask,dets_in(1,1,i),exc,degree,phase,N_int) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + call i_H_j(dets_in(1,1,i),dets_in(1,1,i),N_int,hij) + diag_H_elements(i) = hij + if(degree == 0)then + index_hf = i + else if (degree == 1)then + n_singles +=1 + index_singles(n_singles) = i + ! h1 = inactive orbital of the hole + hole_particles_singles(n_singles,1) = h1 + ! p1 = virtual orbital of the particle + hole_particles_singles(n_singles,2) = p1 + ! s1 = spin of the electron excited + hole_particles_singles(n_singles,3) = s1 + else if (degree == 2)then + n_doubles +=1 + index_doubles(n_doubles) = i + ! h1 = inactive orbital of the hole (beta of course) + hole_particles_doubles(n_doubles,1) = h1 + ! p1 = virtual orbital of the particle (alpha of course) + hole_particles_doubles(n_doubles,2) = p2 + else + print*,'PB !! found out other thing than a single or double' + print*,'stopping ..' + stop + endif + enddo + + e_corr_singles = 0.d0 + e_corr_doubles = 0.d0 + e_corr_singles_total = 0.d0 + e_corr_doubles_1h1p = 0.d0 + c_ref = 1.d0/u_in(index_hf,1) + print*,'c_ref = ',c_ref + do i = 1,sze + call get_excitation(ref_bitmask,dets_in(1,1,i),exc,degree,phase,N_int) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + call i_H_j(ref_bitmask,dets_in(1,1,i),N_int,hij) + contrib = hij * u_in(i,1) * c_ref + if (degree == 1)then + e_corr_singles(h1,s1) += contrib + e_corr_singles(p1,s1) += contrib + e_corr_singles_total(s1)+= contrib + else if (degree == 2)then + e_corr_doubles_1h1p += contrib + e_corr_doubles(h1) += contrib + e_corr_doubles(p2) += contrib + endif + enddo + print*,'e_corr_singles alpha = ',e_corr_singles_total(1) + print*,'e_corr_singles beta = ',e_corr_singles_total(2) + print*,'e_corr_doubles_1h1p = ',e_corr_doubles_1h1p + + ! repeat all the correlation energy on the singles + do i = 1,n_singles + ! you can repeat all the correlation energy of the single excitation of the other spin + diag_H_elements(index_singles(i)) += e_corr_singles_total(other_spin(hole_particles_singles(i,3))) + + ! you can repeat all the correlation energy of the single excitation of the same spin + do j = 1, n_inact_orb + iorb = list_inact(j) + ! except the one of the hole + if(iorb == hole_particles_singles(i,1))cycle + ! ispin = hole_particles_singles(i,3) + diag_H_elements(index_singles(i)) += e_corr_singles(iorb,hole_particles_singles(i,3)) + enddo + ! also exclude all the energy coming from the virtual orbital + diag_H_elements(index_singles(i)) -= e_corr_singles(hole_particles_singles(i,2),hole_particles_singles(i,3)) + + ! If it is a single excitation alpha, you can repeat : + ! +) all the double excitation 1h1p, appart the part involving the virtual orbital "r" + ! If it is a single excitation alpha, you can repeat : + ! +) all the double excitation 1h1p, appart the part involving the inactive orbital "i" + diag_H_elements(index_singles(i)) += e_corr_doubles_1h1p + if(hole_particles_singles(i,3) == 1)then ! alpha single excitation + diag_H_elements(index_singles(i)) -= e_corr_doubles(hole_particles_singles(i,2)) + else ! beta single exctitation + diag_H_elements(index_singles(i)) -= e_corr_doubles(hole_particles_singles(i,1)) + endif + enddo + + ! repeat all the correlation energy on the doubles + ! as all the doubles involve the active space, you cannot repeat any of them one on another + do i = 1, n_doubles + ! on a given double, you can repeat all the correlation energy of the singles alpha + do j = 1, n_inact_orb + iorb = list_inact(j) + ! ispin = hole_particles_singles(i,3) + diag_H_elements(index_doubles(i)) += e_corr_singles(iorb,1) + enddo + ! except the part involving the virtual orbital "hole_particles_doubles(i,2)" + diag_H_elements(index_doubles(i)) -= e_corr_singles(hole_particles_doubles(i,2),1) + ! on a given double, you can repeat all the correlation energy of the singles beta + do j = 1, n_inact_orb + iorb = list_inact(j) + ! except the one of the hole + if(iorb == hole_particles_doubles(i,1))cycle + ! ispin = hole_particles_singles(i,3) + diag_H_elements(index_doubles(i)) += e_corr_singles(iorb,2) + enddo + enddo + + + ! Taking into account the connected part of the 2h2p on the HF determinant + ! 1/2 \sum_{ir,js} c_{ir}^{sigma} c_{js}^{sigma} + +! diag_H_elements(index_hf) += total_corr_e_2h2p + c_ref = c_ref * c_ref + print*,'diag_H_elements(index_hf) = ',diag_H_elements(index_hf) + do i = 1, n_singles + ! start on the single excitation "|i>" + h1 = hole_particles_singles(i,1) + p1 = hole_particles_singles(i,2) + do j = 1, n_singles + do k = 1, N_int + key_tmp(k,1) = dets_in(k,1,index_singles(i)) + key_tmp(k,2) = dets_in(k,2,index_singles(i)) + enddo + h2 = hole_particles_singles(j,1) + p2 = hole_particles_singles(j,2) + call do_mono_excitation(key_tmp,h2,p2,hole_particles_singles(j,3),i_ok) + ! apply the excitation operator from the single excitation "|j>" + if(i_ok .ne. 1)cycle + double precision :: phase_ref_other_single,diag_H_mat_elem,hijj,contrib_e2,coef_1 + call get_excitation(key_tmp,dets_in(1,1,index_singles(i)),exc,degree,phase_single_double,N_int) + call get_excitation(ref_bitmask,dets_in(1,1,index_singles(j)),exc,degree,phase_ref_other_single,N_int) + call i_H_j(ref_bitmask,key_tmp,N_int,hij) + diag_H_elements(index_hf) += u_in(index_singles(i),1) * u_in(index_singles(j),1) * c_ref * hij & + * phase_single_double * phase_ref_other_single + enddo + enddo + print*,'diag_H_elements(index_hf) = ',diag_H_elements(index_hf) + +end + +subroutine dressing_1h1p_full(dets_in,u_in,H_matrix,dim_in,sze,N_st,Nint,convergence) + use bitmasks + implicit none + BEGIN_DOC + ! CISD+SC2 method :: take off all the disconnected terms of a ROHF+1h1p (selected or not) + ! + ! dets_in : bitmasks corresponding to determinants + ! + ! u_in : guess coefficients on the various states. Overwritten + ! on exit + ! + ! dim_in : leftmost dimension of u_in + ! + ! sze : Number of determinants + ! + ! N_st : Number of eigenstates + ! + ! Initial guess vectors are not necessarily orthonormal + END_DOC + integer, intent(in) :: dim_in, sze, N_st, Nint + integer(bit_kind), intent(in) :: dets_in(Nint,2,sze) + double precision, intent(in) :: u_in(dim_in,N_st) + double precision, intent(inout) :: H_matrix(sze,sze) + double precision, intent(in) :: convergence + + integer :: i,j,k,l + integer :: n_singles + integer :: index_singles(sze),hole_particles_singles(sze,3) + integer :: n_doubles + integer :: index_doubles(sze),hole_particles_doubles(sze,2) + integer :: index_hf + double precision :: e_corr_singles(mo_tot_num,2) + double precision :: e_corr_doubles(mo_tot_num) + double precision :: e_corr_singles_total(2) + double precision :: e_corr_doubles_1h1p + + integer :: exc(0:2,2,2),degree + integer :: h1,h2,p1,p2,s1,s2 + integer :: other_spin(2) + double precision :: phase + integer(bit_kind) :: key_tmp(N_int,2) + integer :: i_ok + double precision :: phase_single_double,phase_double_hf,get_mo_bielec_integral_schwartz + double precision :: hij,c_ref,contrib + integer :: iorb + + other_spin(1) = 2 + other_spin(2) = 1 + + n_singles = 0 + n_doubles = 0 + do i = 1,sze + call get_excitation(ref_bitmask,dets_in(1,1,i),exc,degree,phase,N_int) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + if(degree == 0)then + index_hf = i + else if (degree == 1)then + n_singles +=1 + index_singles(n_singles) = i + ! h1 = inactive orbital of the hole + hole_particles_singles(n_singles,1) = h1 + ! p1 = virtual orbital of the particle + hole_particles_singles(n_singles,2) = p1 + ! s1 = spin of the electron excited + hole_particles_singles(n_singles,3) = s1 + else if (degree == 2)then + n_doubles +=1 + index_doubles(n_doubles) = i + ! h1 = inactive orbital of the hole (beta of course) + hole_particles_doubles(n_doubles,1) = h1 + ! p1 = virtual orbital of the particle (alpha of course) + hole_particles_doubles(n_doubles,2) = p2 + else + print*,'PB !! found out other thing than a single or double' + print*,'stopping ..' + stop + endif + enddo + double precision, allocatable :: dressing_H_mat_elem(:) + allocate(dressing_H_mat_elem(N_det)) + logical :: lmct + dressing_H_mat_elem = 0.d0 + call dress_diag_elem_2h2p(dressing_H_mat_elem,N_det) + lmct = .False. + call dress_diag_elem_2h1p(dressing_H_mat_elem,N_det,lmct,1000) + lmct = .true. + call dress_diag_elem_1h2p(dressing_H_mat_elem,N_det,lmct,1000) + do i = 1, N_det + H_matrix(i,i) += dressing_H_mat_elem(i) + enddo + + e_corr_singles = 0.d0 + e_corr_doubles = 0.d0 + e_corr_singles_total = 0.d0 + e_corr_doubles_1h1p = 0.d0 + c_ref = 1.d0/u_in(index_hf,1) + print*,'c_ref = ',c_ref + do i = 1,sze + call get_excitation(ref_bitmask,dets_in(1,1,i),exc,degree,phase,N_int) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + call i_H_j(ref_bitmask,dets_in(1,1,i),N_int,hij) + contrib = hij * u_in(i,1) * c_ref + if (degree == 1)then + e_corr_singles(h1,s1) += contrib + e_corr_singles(p1,s1) += contrib + e_corr_singles_total(s1)+= contrib + else if (degree == 2)then + e_corr_doubles_1h1p += contrib + e_corr_doubles(h1) += contrib + e_corr_doubles(p2) += contrib + endif + enddo + print*,'e_corr_singles alpha = ',e_corr_singles_total(1) + print*,'e_corr_singles beta = ',e_corr_singles_total(2) + print*,'e_corr_doubles_1h1p = ',e_corr_doubles_1h1p + + + ! repeat all the correlation energy on the singles +! do i = 1,n_singles +! ! you can repeat all the correlation energy of the single excitation of the other spin +! H_matrix(index_singles(i),index_singles(i)) += e_corr_singles_total(other_spin(hole_particles_singles(i,3))) + +! ! you can repeat all the correlation energy of the single excitation of the same spin +! do j = 1, n_inact_orb +! iorb = list_inact(j) +! ! except the one of the hole +! if(iorb == hole_particles_singles(i,1))cycle +! ! ispin = hole_particles_singles(i,3) +! H_matrix(index_singles(i),index_singles(i)) += e_corr_singles(iorb,hole_particles_singles(i,3)) +! enddo +! ! also exclude all the energy coming from the virtual orbital +! H_matrix(index_singles(i),index_singles(i)) -= e_corr_singles(hole_particles_singles(i,2),hole_particles_singles(i,3)) +! +! ! If it is a single excitation alpha, you can repeat : +! ! +) all the double excitation 1h1p, appart the part involving the virtual orbital "r" +! ! If it is a single excitation alpha, you can repeat : +! ! +) all the double excitation 1h1p, appart the part involving the inactive orbital "i" +! H_matrix(index_singles(i),index_singles(i)) += e_corr_doubles_1h1p +! if(hole_particles_singles(i,3) == 1)then ! alpha single excitation +! H_matrix(index_singles(i),index_singles(i)) -= e_corr_doubles(hole_particles_singles(i,2)) +! else ! beta single exctitation +! H_matrix(index_singles(i),index_singles(i)) -= e_corr_doubles(hole_particles_singles(i,1)) +! endif +! enddo + +! ! repeat all the correlation energy on the doubles +! ! as all the doubles involve the active space, you cannot repeat any of them one on another +! do i = 1, n_doubles +! ! on a given double, you can repeat all the correlation energy of the singles alpha +! do j = 1, n_inact_orb +! iorb = list_inact(j) +! ! ispin = hole_particles_singles(i,3) +! H_matrix(index_doubles(i),index_doubles(i)) += e_corr_singles(iorb,1) +! enddo +! ! except the part involving the virtual orbital "hole_particles_doubles(i,2)" +! H_matrix(index_doubles(i),index_doubles(i)) -= e_corr_singles(hole_particles_doubles(i,2),1) +! ! on a given double, you can repeat all the correlation energy of the singles beta +! do j = 1, n_inact_orb +! iorb = list_inact(j) +! ! except the one of the hole +! if(iorb == hole_particles_doubles(i,1))cycle +! ! ispin = hole_particles_singles(i,3) +! H_matrix(index_doubles(i),index_doubles(i)) += e_corr_singles(iorb,2) +! enddo +! enddo + + + ! Taking into account the connected part of the 2h2p on the HF determinant + ! 1/2 \sum_{ir,js} c_{ir}^{sigma} c_{js}^{sigma} + +! H_matrix(index_hf) += total_corr_e_2h2p + print*,'H_matrix(index_hf,index_hf) = ',H_matrix(index_hf,index_hf) + do i = 1, n_singles + ! start on the single excitation "|i>" + h1 = hole_particles_singles(i,1) + p1 = hole_particles_singles(i,2) + print*,'i = ',i + do j = i+1, n_singles + do k = 1, N_int + key_tmp(k,1) = dets_in(k,1,index_singles(i)) + key_tmp(k,2) = dets_in(k,2,index_singles(i)) + enddo + h2 = hole_particles_singles(j,1) + p2 = hole_particles_singles(j,2) + call do_mono_excitation(key_tmp,h2,p2,hole_particles_singles(j,3),i_ok) + ! apply the excitation operator from the single excitation "|j>" + if(i_ok .ne. 1)cycle + double precision :: H_array(sze),diag_H_mat_elem,hjj + do k = 1, sze + call get_excitation_degree(dets_in(1,1,k),key_tmp,degree,N_int) + H_array(k) = 0.d0 + if(degree > 2)cycle + call i_H_j(dets_in(1,1,k),key_tmp,N_int,hij) + H_array(k) = hij + enddo + hjj = 1.d0/(ref_bitmask_energy - diag_H_mat_elem(key_tmp,N_int)) +! contrib_e2 = 0.5d0 * (delta_e + dsqrt(delta_e * delta_e + 4.d0 * hij * hij)) + do l = 2, sze +! pause + H_matrix(l,l) += H_array(l) * H_array(l) * hjj +! H_matrix(1,l) += H_array(1) * H_array(l) * hjj +! H_matrix(l,1) += H_array(1) * H_array(l) * hjj + enddo + enddo + enddo + print*,'H_matrix(index_hf,index_hf) = ',H_matrix(index_hf,index_hf) + +end + +subroutine SC2_1h1p_full(dets_in,u_in,energies,H_matrix,dim_in,sze,N_st,Nint,convergence) + use bitmasks + implicit none + BEGIN_DOC + ! CISD+SC2 method :: take off all the disconnected terms of a CISD (selected or not) + ! + ! dets_in : bitmasks corresponding to determinants + ! + ! u_in : guess coefficients on the various states. Overwritten + ! on exit + ! + ! dim_in : leftmost dimension of u_in + ! + ! sze : Number of determinants + ! + ! N_st : Number of eigenstates + ! + ! Initial guess vectors are not necessarily orthonormal + END_DOC + integer, intent(in) :: dim_in, sze, N_st, Nint + integer(bit_kind), intent(in) :: dets_in(Nint,2,sze) + double precision, intent(inout) :: u_in(dim_in,N_st) + double precision, intent(out) :: energies(N_st) + double precision, intent(out) :: H_matrix(sze,sze) + double precision, intent(in) :: convergence + integer :: i,j,iter + print*,'sze = ',sze + do iter = 1, 1 +! if(sze<=N_det_max_jacobi)then + double precision, allocatable :: eigenvectors(:,:), eigenvalues(:),H_matrix_tmp(:,:) + allocate (H_matrix_tmp(size(H_matrix_all_dets,1),sze),eigenvalues(sze),eigenvectors(size(H_matrix_all_dets,1),sze)) + H_matrix_tmp = 0.d0 + call dressing_1h1p_full(dets_in,u_in,H_matrix_tmp,dim_in,sze,N_st,Nint,convergence) + do j=1,sze + do i=1,sze + H_matrix_tmp(i,j) += H_matrix_all_dets(i,j) + enddo + enddo + print*,'passed the dressing' + call lapack_diag(eigenvalues,eigenvectors, & + H_matrix_tmp,size(H_matrix_all_dets,1),sze) + do j=1,min(N_states_diag,sze) + do i=1,sze + u_in(i,j) = eigenvectors(i,j) + enddo + energies(j) = eigenvalues(j) + enddo + deallocate (H_matrix_tmp, eigenvalues, eigenvectors) +! else +! call davidson_diag_hjj(dets_in,u_in,diag_H_elements,energies,dim_in,sze,N_st,Nint,output_determinants) +! endif + print*,'E = ',energies(1) + nuclear_repulsion + + enddo + + +end + + +subroutine SC2_1h1p(dets_in,u_in,energies,diag_H_elements,dim_in,sze,N_st,Nint,convergence) + use bitmasks + implicit none + BEGIN_DOC + ! CISD+SC2 method :: take off all the disconnected terms of a CISD (selected or not) + ! + ! dets_in : bitmasks corresponding to determinants + ! + ! u_in : guess coefficients on the various states. Overwritten + ! on exit + ! + ! dim_in : leftmost dimension of u_in + ! + ! sze : Number of determinants + ! + ! N_st : Number of eigenstates + ! + ! Initial guess vectors are not necessarily orthonormal + END_DOC + integer, intent(in) :: dim_in, sze, N_st, Nint + integer(bit_kind), intent(in) :: dets_in(Nint,2,sze) + double precision, intent(inout) :: u_in(dim_in,N_st) + double precision, intent(out) :: energies(N_st) + double precision, intent(out) :: diag_H_elements(dim_in) + double precision, intent(in) :: convergence + integer :: i,j,iter + do iter = 1, 1 + call dressing_1h1p(dets_in,u_in,diag_H_elements,dim_in,sze,N_st,Nint,convergence) + if(sze<=N_det_max_jacobi)then + double precision, allocatable :: eigenvectors(:,:), eigenvalues(:),H_matrix_tmp(:,:) + allocate (H_matrix_tmp(size(H_matrix_all_dets,1),sze),eigenvalues(sze),eigenvectors(size(H_matrix_all_dets,1),sze)) + do j=1,sze + do i=1,sze + H_matrix_tmp(i,j) = H_matrix_all_dets(i,j) + enddo + enddo + do i = 1,sze + H_matrix_tmp(i,i) = diag_H_elements(i) + enddo + call lapack_diag(eigenvalues,eigenvectors, & + H_matrix_tmp,size(H_matrix_all_dets,1),sze) + do j=1,min(N_states_diag,sze) + do i=1,sze + u_in(i,j) = eigenvectors(i,j) + enddo + energies(j) = eigenvalues(j) + enddo + deallocate (H_matrix_tmp, eigenvalues, eigenvectors) + else + call davidson_diag_hjj(dets_in,u_in,diag_H_elements,energies,dim_in,sze,N_st,Nint,output_determinants) + endif + print*,'E = ',energies(1) + nuclear_repulsion + + enddo + + +end + + +subroutine density_matrix_1h1p(dets_in,u_in,density_matrix_alpha,density_matrix_beta,norm,dim_in,sze,N_st,Nint) + use bitmasks + implicit none + BEGIN_DOC + ! CISD+SC2 method :: take off all the disconnected terms of a ROHF+1h1p (selected or not) + ! + ! dets_in : bitmasks corresponding to determinants + ! + ! u_in : guess coefficients on the various states. Overwritten + ! on exit + ! + ! dim_in : leftmost dimension of u_in + ! + ! sze : Number of determinants + ! + ! N_st : Number of eigenstates + ! + ! Initial guess vectors are not necessarily orthonormal + END_DOC + integer, intent(in) :: dim_in, sze, N_st, Nint + integer(bit_kind), intent(in) :: dets_in(Nint,2,sze) + double precision, intent(inout) :: u_in(dim_in,N_st) + double precision, intent(inout) :: density_matrix_alpha(mo_tot_num_align,mo_tot_num) + double precision, intent(inout) :: density_matrix_beta(mo_tot_num_align,mo_tot_num) + double precision, intent(inout) :: norm + + integer :: i,j,k,l + integer :: n_singles + integer :: index_singles(sze),hole_particles_singles(sze,3) + integer :: n_doubles + integer :: index_doubles(sze),hole_particles_doubles(sze,2) + integer :: index_hf + + integer :: exc(0:2,2,2),degree + integer :: h1,h2,p1,p2,s1,s2 + integer :: other_spin(2) + double precision :: phase + integer(bit_kind) :: key_tmp(N_int,2) + integer :: i_ok + double precision :: phase_single_double,phase_double_hf,get_mo_bielec_integral_schwartz + double precision :: hij,c_ref,contrib + integer :: iorb + + other_spin(1) = 2 + other_spin(2) = 1 + + n_singles = 0 + n_doubles = 0 + norm = 0.d0 + do i = 1,sze + call get_excitation(ref_bitmask,dets_in(1,1,i),exc,degree,phase,N_int) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + norm += u_in(i,1)* u_in(i,1) + if(degree == 0)then + index_hf = i + c_ref = 1.d0/psi_coef(i,1) + else if (degree == 1)then + n_singles +=1 + index_singles(n_singles) = i + ! h1 = inactive orbital of the hole + hole_particles_singles(n_singles,1) = h1 + ! p1 = virtual orbital of the particle + hole_particles_singles(n_singles,2) = p1 + ! s1 = spin of the electron excited + hole_particles_singles(n_singles,3) = s1 + else if (degree == 2)then + n_doubles +=1 + index_doubles(n_doubles) = i + ! h1 = inactive orbital of the hole (beta of course) + hole_particles_doubles(n_doubles,1) = h1 + ! p1 = virtual orbital of the particle (alpha of course) + hole_particles_doubles(n_doubles,2) = p2 + else + print*,'PB !! found out other thing than a single or double' + print*,'stopping ..' + stop + endif + enddo + print*,'norm = ',norm + + ! Taking into account the connected part of the 2h2p on the HF determinant + ! 1/2 \sum_{ir,js} c_{ir}^{sigma} c_{js}^{sigma} + + do i = 1, n_singles + ! start on the single excitation "|i>" + h1 = hole_particles_singles(i,1) + p1 = hole_particles_singles(i,2) + do j = 1, n_singles + do k = 1, N_int + key_tmp(k,1) = dets_in(k,1,index_singles(i)) + key_tmp(k,2) = dets_in(k,2,index_singles(i)) + enddo + h2 = hole_particles_singles(j,1) + p2 = hole_particles_singles(j,2) + call do_mono_excitation(key_tmp,h2,p2,hole_particles_singles(j,3),i_ok) + ! apply the excitation operator from the single excitation "|j>" + if(i_ok .ne. 1)cycle + double precision :: coef_ijrs,phase_other_single_ref + integer :: occ(N_int*bit_kind_size,2),n_occ(2) + call get_excitation(key_tmp,dets_in(1,1,index_singles(i)),exc,degree,phase_single_double,N_int) + call get_excitation(ref_bitmask,dets_in(1,1,index_singles(j)),exc,degree,phase_other_single_ref,N_int) + call get_excitation(key_tmp,dets_in(1,1,index_singles(j)),exc,degree,phase_other_single_ref,N_int) + coef_ijrs = u_in(index_singles(i),1) * u_in(index_singles(j),1) * c_ref * c_ref & + * phase_single_double * phase_other_single_ref + call bitstring_to_list_ab(key_tmp, occ, n_occ, N_int) + do k=1,elec_alpha_num + l = occ(k,1) + density_matrix_alpha(l,l) += coef_ijrs*coef_ijrs + enddo + do k=1,elec_beta_num + l = occ(k,1) + density_matrix_beta(l,l) += coef_ijrs*coef_ijrs + enddo + norm += coef_ijrs* coef_ijrs + if(hole_particles_singles(j,3) == 1)then ! single alpha + density_matrix_alpha(h2,p2) += coef_ijrs * phase_single_double * u_in(index_singles(i),1) * c_ref + density_matrix_alpha(p2,h2) += coef_ijrs * phase_single_double * u_in(index_singles(i),1) * c_ref + else + density_matrix_beta(h2,p2) += coef_ijrs * phase_single_double * u_in(index_singles(i),1) * c_ref + density_matrix_beta(p2,h2) += coef_ijrs * phase_single_double * u_in(index_singles(i),1) * c_ref + endif + enddo + enddo + + + do i = 1, n_doubles + ! start on the double excitation "|i>" + h1 = hole_particles_doubles(i,1) + p1 = hole_particles_doubles(i,2) + do j = 1, n_singles + do k = 1, N_int + key_tmp(k,1) = dets_in(k,1,index_doubles(i)) + key_tmp(k,2) = dets_in(k,2,index_doubles(i)) + enddo + h2 = hole_particles_singles(j,1) + p2 = hole_particles_singles(j,2) + call do_mono_excitation(key_tmp,h2,p2,hole_particles_singles(j,3),i_ok) + ! apply the excitation operator from the single excitation "|j>" + if(i_ok .ne. 1)cycle + double precision :: coef_ijrs_kv,phase_double_triple + call get_excitation(key_tmp,dets_in(1,1,index_singles(i)),exc,degree,phase_double_triple,N_int) + call get_excitation(ref_bitmask,dets_in(1,1,index_singles(j)),exc,degree,phase_other_single_ref,N_int) + call get_excitation(key_tmp,dets_in(1,1,index_singles(j)),exc,degree,phase_other_single_ref,N_int) + coef_ijrs_kv = u_in(index_doubles(i),1) * u_in(index_singles(j),1) * c_ref * c_ref & + * phase_double_triple * phase_other_single_ref + call bitstring_to_list_ab(key_tmp, occ, n_occ, N_int) + do k=1,elec_alpha_num + l = occ(k,1) + density_matrix_alpha(l,l) += coef_ijrs_kv*coef_ijrs_kv + enddo + do k=1,elec_beta_num + l = occ(k,1) + density_matrix_beta(l,l) += coef_ijrs_kv*coef_ijrs_kv + enddo + norm += coef_ijrs_kv* coef_ijrs_kv + if(hole_particles_singles(j,3) == 1)then ! single alpha + density_matrix_alpha(h2,p2) += coef_ijrs_kv * phase_double_triple * u_in(index_doubles(i),1) * c_ref + density_matrix_alpha(p2,h2) += coef_ijrs_kv * phase_double_triple * u_in(index_doubles(i),1) * c_ref + else + density_matrix_beta(h2,p2) += coef_ijrs_kv * phase_double_triple * u_in(index_doubles(i),1) * c_ref + density_matrix_beta(p2,h2) += coef_ijrs_kv * phase_double_triple * u_in(index_doubles(i),1) * c_ref + endif + enddo + enddo + + + + + print*,'norm = ',norm + norm = 1.d0/norm + do i = 1, mo_tot_num + do j = 1, mo_tot_num + density_matrix_alpha(i,j) *= norm + density_matrix_beta(i,j) *= norm + enddo + enddo + coef_ijrs = 0.d0 + do i = 1, mo_tot_num + coef_ijrs += density_matrix_beta(i,i) + density_matrix_beta(i,i) + enddo + print*,'accu = ',coef_ijrs + +end + diff --git a/plugins/FOBOCI/fobo_scf.irp.f b/plugins/FOBOCI/fobo_scf.irp.f index 0b0902b0..8656b633 100644 --- a/plugins/FOBOCI/fobo_scf.irp.f +++ b/plugins/FOBOCI/fobo_scf.irp.f @@ -1,6 +1,6 @@ program foboscf implicit none -!call run_prepare + call run_prepare no_oa_or_av_opt = .True. touch no_oa_or_av_opt call routine_fobo_scf diff --git a/plugins/FOBOCI/foboci_lmct_mlct_threshold_old.irp.f b/plugins/FOBOCI/foboci_lmct_mlct_threshold_old.irp.f index dc6519b8..e81b3fc1 100644 --- a/plugins/FOBOCI/foboci_lmct_mlct_threshold_old.irp.f +++ b/plugins/FOBOCI/foboci_lmct_mlct_threshold_old.irp.f @@ -81,6 +81,8 @@ subroutine FOBOCI_lmct_mlct_old_thr call set_bitmask_particl_as_input(reunion_of_bitmask) call set_bitmask_hole_as_input(reunion_of_bitmask) call all_single + call make_s2_eigenfunction + call diagonalize_ci ! if(dressing_2h2p)then ! call diag_dressed_2h2p_hamiltonian_and_update_psi_det(i_hole_osoci,lmct) ! endif @@ -193,6 +195,8 @@ subroutine FOBOCI_lmct_mlct_old_thr ! call all_single_split(psi_det_generators,psi_coef_generators,N_det_generators,dressing_matrix) ! call diag_dressed_matrix_and_set_to_psi_det(psi_det_generators,N_det_generators,dressing_matrix) call all_single + call make_s2_eigenfunction + call diagonalize_ci ! if(dressing_2h2p)then ! call diag_dressed_2h2p_hamiltonian_and_update_psi_det(i_particl_osoci,lmct) ! endif diff --git a/plugins/Properties/delta_rho.irp.f b/plugins/Properties/delta_rho.irp.f index 69894c38..7803ba3d 100644 --- a/plugins/Properties/delta_rho.irp.f +++ b/plugins/Properties/delta_rho.irp.f @@ -3,9 +3,9 @@ &BEGIN_PROVIDER [double precision, z_max] &BEGIN_PROVIDER [double precision, delta_z] implicit none - z_min = -20.d0 - z_max = 20.d0 - delta_z = 0.1d0 + z_min = 0.d0 + z_max = 10.d0 + delta_z = 0.005d0 N_z_pts = (z_max - z_min)/delta_z print*,'N_z_pts = ',N_z_pts diff --git a/plugins/Properties/mulliken.irp.f b/plugins/Properties/mulliken.irp.f index cc0a2f8e..deeb90bf 100644 --- a/plugins/Properties/mulliken.irp.f +++ b/plugins/Properties/mulliken.irp.f @@ -14,13 +14,16 @@ BEGIN_PROVIDER [double precision, spin_population, (ao_num_align,ao_num)] enddo END_PROVIDER -BEGIN_PROVIDER [double precision, spin_population_angular_momentum, (0:ao_l_max)] + BEGIN_PROVIDER [double precision, spin_population_angular_momentum, (0:ao_l_max)] +&BEGIN_PROVIDER [double precision, spin_population_angular_momentum_per_atom, (0:ao_l_max,nucl_num)] implicit none integer :: i double precision :: accu spin_population_angular_momentum = 0.d0 + spin_population_angular_momentum_per_atom = 0.d0 do i = 1, ao_num spin_population_angular_momentum(ao_l(i)) += spin_gross_orbital_product(i) + spin_population_angular_momentum_per_atom(ao_l(i),ao_nucl(i)) += spin_gross_orbital_product(i) enddo END_PROVIDER @@ -133,6 +136,16 @@ subroutine print_mulliken_sd print*,' ',trim(l_to_charater(i)),spin_population_angular_momentum(i) print*,'sum = ',accu enddo + print*,'Angular momentum analysis per atom' + print*,'Angular momentum analysis' + do j = 1,nucl_num + accu = 0.d0 + do i = 0, ao_l_max + accu += spin_population_angular_momentum_per_atom(i,j) + write(*,'(XX,I3,XX,A4,X,A4,X,F10.7)')j,trim(element_name(int(nucl_charge(j)))),trim(l_to_charater(i)),spin_population_angular_momentum_per_atom(i,j) + print*,'sum = ',accu + enddo + enddo end diff --git a/plugins/Properties/print_spin_density.irp.f b/plugins/Properties/print_spin_density.irp.f new file mode 100644 index 00000000..9daa6fb7 --- /dev/null +++ b/plugins/Properties/print_spin_density.irp.f @@ -0,0 +1,36 @@ +program print_sd + implicit none + read_wf = .True. + touch read_wf + call routine + +end + +subroutine routine + implicit none + integer :: i,j,k + double precision :: z + double precision :: r(3),accu,accu_alpha,accu_beta,tmp + double precision, allocatable :: aos_array(:) + allocate(aos_array(ao_num)) + r = 0.d0 + r(3) = z_min + do i = 1, N_z_pts + call give_all_aos_at_r(r,aos_array) + accu = 0.d0 + accu_alpha = 0.d0 + accu_beta = 0.d0 + do j = 1, ao_num + do k = 1, ao_num + tmp = aos_array(k) * aos_array(j) + accu += one_body_spin_density_ao(k,j) * tmp + accu_alpha += one_body_dm_ao_alpha(k,j) * tmp + accu_beta += one_body_dm_ao_beta(k,j) * tmp + enddo + enddo + r(3) += delta_z + write(33,'(100(f16.10,X))')r(3),accu,accu_alpha,accu_beta + enddo + + +end diff --git a/plugins/Properties/provide_deltarho.irp.f b/plugins/Properties/provide_deltarho.irp.f new file mode 100644 index 00000000..d576d622 --- /dev/null +++ b/plugins/Properties/provide_deltarho.irp.f @@ -0,0 +1,11 @@ +program pouet + implicit none + read_wf = .True. + touch read_wf + call routine +end + +subroutine routine + implicit none + provide integrated_delta_rho_all_points +end diff --git a/plugins/Properties/test_two_body_dm.irp.f b/plugins/Properties/test_two_body_dm.irp.f index 6fc02abf..f79aed86 100644 --- a/plugins/Properties/test_two_body_dm.irp.f +++ b/plugins/Properties/test_two_body_dm.irp.f @@ -7,28 +7,67 @@ end subroutine routine implicit none integer :: i,j,k,l + integer :: h1,p1,h2,p2,s1,s2 double precision :: accu,get_two_body_dm_ab_map_element,get_mo_bielec_integral_schwartz accu = 0.d0 - ! Diag part of the two body dm + ! Diag part of the core two body dm + do i = 1, n_core_orb + h1 = list_core(i) + do j = 1, n_core_orb + h2 = list_core(j) + accu += two_body_dm_ab_diag_core(j,i) * mo_bielec_integral_jj(h1,h2) + enddo + enddo + + ! Diag part of the active two body dm do i = 1, n_act_orb + h1 = list_act(i) do j = 1, n_act_orb - accu += two_body_dm_ab_diag(i,j) * mo_bielec_integral_jj(i,j) + h2 = list_act(j) + accu += two_body_dm_ab_diag_act(j,i) * mo_bielec_integral_jj(h1,h2) + enddo + enddo + + ! Diag part of the active <-> core two body dm + do i = 1, n_act_orb + h1 = list_act(i) + do j = 1, n_core_orb + h2 = list_core(j) + accu += two_body_dm_diag_core_act(j,i) * mo_bielec_integral_jj(h1,h2) enddo enddo print*,'BI ELECTRONIC = ',accu double precision :: accu_extra_diag accu_extra_diag = 0.d0 + ! purely active part of the two body dm do l = 1, n_act_orb ! p2 + p2 = list_act(l) do k = 1, n_act_orb ! h2 + h2 = list_act(k) do j = 1, n_act_orb ! p1 + p1 = list_act(j) do i = 1,n_act_orb ! h1 - accu_extra_diag += two_body_dm_ab_big_array(i,j,k,l) * get_mo_bielec_integral_schwartz(i,k,j,l,mo_integrals_map) + h1 = list_act(i) + accu_extra_diag += two_body_dm_ab_big_array_act(i,j,k,l) * get_mo_bielec_integral_schwartz(h1,h2,p1,p2,mo_integrals_map) enddo enddo enddo enddo + + ! core <-> active part of the two body dm + do l = 1, n_act_orb ! p1 + p1 = list_act(l) + do k = 1, n_act_orb ! h1 + h1 = list_act(k) + do i = 1,n_core_orb ! h2 + h2 = list_core(i) + accu_extra_diag += two_body_dm_ab_big_array_core_act(i,k,l) * get_mo_bielec_integral_schwartz(h1,h2,p1,h2,mo_integrals_map) + enddo + enddo + enddo + print*,'extra_diag = ',accu_extra_diag double precision :: average_mono call get_average(mo_mono_elec_integral,one_body_dm_mo,average_mono) @@ -41,7 +80,6 @@ subroutine routine print*,' = ',e_0 + nuclear_repulsion integer :: degree integer :: exc(0:2,2,2) - integer :: h1,h2,p1,p2,s1,s2 double precision :: phase integer :: n_elements n_elements = 0 diff --git a/plugins/loc_cele/loc_cele.irp.f b/plugins/loc_cele/loc_cele.irp.f index c9036aa1..8a110c05 100644 --- a/plugins/loc_cele/loc_cele.irp.f +++ b/plugins/loc_cele/loc_cele.irp.f @@ -92,7 +92,7 @@ - nrot(1) = 64 ! number of orbitals to be localized + nrot(1) = 6 ! number of orbitals to be localized integer :: index_rot(1000,1) @@ -101,261 +101,72 @@ cmoref = 0.d0 irot = 0 -! H2 molecule for the mixed localization - do i=1,64 - irot(i,1) = i+2 + do i=1,nrot(1) + irot(i,1) = 19+i enddo - do i=1,17 - cmoref(i+1,i,1)=1.d0 - enddo - cmoref(19,19-1,1)=1.d0 - cmoref(20,19-1,1)=-1.d0 - cmoref(19,20-1,1)=-1.d0 - cmoref(20,20-1,1)=-1.d0 - cmoref(21,20-1,1)=2.d0 - cmoref(22,21-1,1)=1.d0 - cmoref(23,22-1,1)=1.d0 - cmoref(24,23-1,1)=1.d0 +! ESATRIENE with 3 bonding and anti bonding orbitals +! First bonding orbital for esa +! cmoref(7,1,1) = 1.d0 ! +! cmoref(26,1,1) = 1.d0 ! +! Second bonding orbital for esa +! cmoref(45,2,1) = 1.d0 ! +! cmoref(64,2,1) = 1.d0 ! +! Third bonding orbital for esa +! cmoref(83,3,1) = 1.d0 ! +! cmoref(102,3,1) = 1.d0 ! + +! First anti bonding orbital for esa +! cmoref(7,4,1) = 1.d0 ! +! cmoref(26,4,1) = -1.d0 ! +! Second anti bonding orbital for esa +! cmoref(45,5,1) = 1.d0 ! +! cmoref(64,5,1) = -1.d0 ! +! Third anti bonding orbital for esa +! cmoref(83,6,1) = 1.d0 ! +! cmoref(102,6,1) = -1.d0 ! + +! ESATRIENE with 2 bonding and anti bonding orbitals +! AND 2 radical orbitals +! First radical orbital +! cmoref(7,1,1) = 1.d0 ! +! First bonding orbital +! cmoref(26,2,1) = 1.d0 ! +! cmoref(45,2,1) = 1.d0 ! +! Second bonding orbital +! cmoref(64,3,1) = 1.d0 ! +! cmoref(83,3,1) = 1.d0 ! +! Second radical orbital for esa +! cmoref(102,4,1) = 1.d0 ! + +! First anti bonding orbital for esa +! cmoref(26,5,1) = 1.d0 ! +! cmoref(45,5,1) =-1.d0 ! +! Second anti bonding orbital for esa +! cmoref(64,6,1) = 1.d0 ! +! cmoref(83,6,1) =-1.d0 ! + +! ESATRIENE with 1 central bonding and anti bonding orbitals +! AND 4 radical orbitals +! First radical orbital + cmoref(7,1,1) = 1.d0 ! +! Second radical orbital + cmoref(26,2,1) = 1.d0 ! +! First bonding orbital + cmoref(45,3,1) = 1.d0 ! + cmoref(64,3,1) = 1.d0 ! +! Third radical orbital for esa + cmoref(83,4,1) = 1.d0 ! +! Fourth radical orbital for esa + cmoref(102,5,1) = 1.d0 ! +! First anti bonding orbital + cmoref(45,6,1) = 1.d0 ! + cmoref(64,6,1) =-1.d0 ! - cmoref(25,24-1,1)=1.d0 - cmoref(26,24-1,1)=-1.d0 - cmoref(25,25-1,1)=-1.d0 - cmoref(26,25-1,1)=-1.d0 - cmoref(27,25-1,1)=2.d0 - cmoref(28,26-1,1)=1.d0 - cmoref(29,27-1,1)=1.d0 - cmoref(30,28-1,1)=1.d0 - - cmoref(31,29-1,1)=1.d0 - cmoref(32,29-1,1)=-1.d0 - cmoref(31,30-1,1)=-1.d0 - cmoref(32,30-1,1)=-1.d0 - cmoref(33,30-1,1)=2.d0 - cmoref(34,31-1,1)=1.d0 - cmoref(35,32-1,1)=1.d0 - cmoref(36,33-1,1)=1.d0 - - do i=33,49 - cmoref(i+5,i,1)= 1.d0 - enddo - - cmoref(55,52-2,1)=1.d0 - cmoref(56,52-2,1)=-1.d0 - cmoref(55,53-2,1)=-1.d0 - cmoref(56,53-2,1)=-1.d0 - cmoref(57,53-2,1)=2.d0 - cmoref(58,54-2,1)=1.d0 - cmoref(59,55-2,1)=1.d0 - cmoref(60,56-2,1)=1.d0 - - cmoref(61,57-2,1)=1.d0 - cmoref(62,57-2,1)=-1.d0 - cmoref(61,58-2,1)=-1.d0 - cmoref(62,58-2,1)=-1.d0 - cmoref(63,58-2,1)=2.d0 - cmoref(64,59-2,1)=1.d0 - cmoref(65,60-2,1)=1.d0 - cmoref(66,61-2,1)=1.d0 - - cmoref(67,62-2,1)=1.d0 - cmoref(68,62-2,1)=-1.d0 - cmoref(67,63-2,1)=-1.d0 - cmoref(68,63-2,1)=-1.d0 - cmoref(69,63-2,1)=2.d0 - cmoref(70,64-2,1)=1.d0 - cmoref(71,65-2,1)=1.d0 - cmoref(72,66-2,1)=1.d0 -! H2 molecule -! do i=1,66 -! irot(i,1) = i -! enddo -! -! do i=1,18 -! cmoref(i,i,1)=1.d0 -! enddo -! cmoref(19,19,1)=1.d0 -! cmoref(20,19,1)=-1.d0 -! cmoref(19,20,1)=-1.d0 -! cmoref(20,20,1)=-1.d0 -! cmoref(21,20,1)=2.d0 -! cmoref(22,21,1)=1.d0 -! cmoref(23,22,1)=1.d0 -! cmoref(24,23,1)=1.d0 -! -! -! cmoref(25,24,1)=1.d0 -! cmoref(26,24,1)=-1.d0 -! cmoref(25,25,1)=-1.d0 -! cmoref(26,25,1)=-1.d0 -! cmoref(27,25,1)=2.d0 -! cmoref(28,26,1)=1.d0 -! cmoref(29,27,1)=1.d0 -! cmoref(30,28,1)=1.d0 -! -! cmoref(31,29,1)=1.d0 -! cmoref(32,29,1)=-1.d0 -! cmoref(31,30,1)=-1.d0 -! cmoref(32,30,1)=-1.d0 -! cmoref(33,30,1)=2.d0 -! cmoref(34,31,1)=1.d0 -! cmoref(35,32,1)=1.d0 -! cmoref(36,33,1)=1.d0 -! -! do i=34,51 -! cmoref(i+3,i,1)= 1.d0 -! enddo -! -! cmoref(55,52,1)=1.d0 -! cmoref(56,52,1)=-1.d0 -! cmoref(55,53,1)=-1.d0 -! cmoref(56,53,1)=-1.d0 -! cmoref(57,53,1)=2.d0 -! cmoref(58,54,1)=1.d0 -! cmoref(59,55,1)=1.d0 -! cmoref(60,56,1)=1.d0 -! -! cmoref(61,57,1)=1.d0 -! cmoref(62,57,1)=-1.d0 -! cmoref(61,58,1)=-1.d0 -! cmoref(62,58,1)=-1.d0 -! cmoref(63,58,1)=2.d0 -! cmoref(64,59,1)=1.d0 -! cmoref(65,60,1)=1.d0 -! cmoref(66,61,1)=1.d0 -! -! cmoref(67,62,1)=1.d0 -! cmoref(68,62,1)=-1.d0 -! cmoref(67,63,1)=-1.d0 -! cmoref(68,63,1)=-1.d0 -! cmoref(69,63,1)=2.d0 -! cmoref(70,64,1)=1.d0 -! cmoref(71,65,1)=1.d0 -! cmoref(72,66,1)=1.d0 -! H atom -! do i=1,33 -! irot(i,1) = i -! enddo -! -! do i=1,18 -! cmoref(i,i,1)=1.d0 -! enddo -! cmoref(19,19,1)=1.d0 -! cmoref(20,19,1)=-1.d0 -! cmoref(19,20,1)=-1.d0 -! cmoref(20,20,1)=-1.d0 -! cmoref(21,20,1)=2.d0 -! cmoref(22,21,1)=1.d0 -! cmoref(23,22,1)=1.d0 -! cmoref(24,23,1)=1.d0 - - -! cmoref(25,24,1)=1.d0 -! cmoref(26,24,1)=-1.d0 -! cmoref(25,25,1)=-1.d0 -! cmoref(26,25,1)=-1.d0 -! cmoref(27,25,1)=2.d0 -! cmoref(28,26,1)=1.d0 -! cmoref(29,27,1)=1.d0 -! cmoref(30,28,1)=1.d0 -! -! cmoref(31,29,1)=1.d0 -! cmoref(32,29,1)=-1.d0 -! cmoref(31,30,1)=-1.d0 -! cmoref(32,30,1)=-1.d0 -! cmoref(33,30,1)=2.d0 -! cmoref(34,31,1)=1.d0 -! cmoref(35,32,1)=1.d0 -! cmoref(36,33,1)=1.d0 - - ! Definition of the index of the MO to be rotated -! irot(2,1) = 21 ! the first mo to be rotated is the 21 th MO -! irot(3,1) = 22 ! etc.... -! irot(4,1) = 23 ! -! irot(5,1) = 24 ! -! irot(6,1) = 25 ! - -!N2 -! irot(1,1) = 5 -! irot(2,1) = 6 -! irot(3,1) = 7 -! irot(4,1) = 8 -! irot(5,1) = 9 -! irot(6,1) = 10 -! -! cmoref(5,1,1) = 1.d0 ! -! cmoref(6,2,1) = 1.d0 ! -! cmoref(7,3,1) = 1.d0 ! -! cmoref(40,4,1) = 1.d0 ! -! cmoref(41,5,1) = 1.d0 ! -! cmoref(42,6,1) = 1.d0 ! -!END N2 - -!HEXATRIENE -! irot(1,1) = 20 -! irot(2,1) = 21 -! irot(3,1) = 22 -! irot(4,1) = 23 -! irot(5,1) = 24 -! irot(6,1) = 25 -! -! cmoref(7,1,1) = 1.d0 ! -! cmoref(26,1,1) = 1.d0 ! -! cmoref(45,2,1) = 1.d0 ! -! cmoref(64,2,1) = 1.d0 ! -! cmoref(83,3,1) = 1.d0 ! -! cmoref(102,3,1) = 1.d0 ! -! cmoref(7,4,1) = 1.d0 ! -! cmoref(26,4,1) = -1.d0 ! -! cmoref(45,5,1) = 1.d0 ! -! cmoref(64,5,1) = -1.d0 ! -! cmoref(83,6,1) = 1.d0 ! -! cmoref(102,6,1) = -1.d0 ! -!END HEXATRIENE - -!!!!H2 H2 CAS -! irot(1,1) = 1 -! irot(2,1) = 2 -! -! cmoref(1,1,1) = 1.d0 -! cmoref(37,2,1) = 1.d0 -!END H2 -!!!! LOCALIZATION ON THE BASIS FUNCTIONS -! do i = 1, nrot(1) -! irot(i,1) = i -! cmoref(i,i,1) = 1.d0 -! enddo - -!END BASISLOC - -! do i = 1, nrot(1) -! irot(i,1) = 4+i -! enddo do i = 1, nrot(1) print*,'irot(i,1) = ',irot(i,1) enddo -! pause - - ! you define the guess vectors that you want - ! the new MO to be close to - ! cmore(i,j,1) = < AO_i | guess_vector_MO(j) > - ! i goes from 1 to ao_num - ! j goes from 1 to nrot(1) - - ! Here you must go to the GAMESS output file - ! where the AOs are listed and explicited - ! From the basis of this knowledge you can build your - ! own guess vectors for the MOs - ! The new MOs are provided in output - ! in the same order than the guess MOs -! do i = 1, nrot(1) -! j = 5+(i-1)*15 -! cmoref(j,i,1) = 0.2d0 -! cmoref(j+3,i,1) = 0.12d0 -! print*,'j = ',j -! enddo -! pause diff --git a/src/Bitmask/bitmasks.irp.f b/src/Bitmask/bitmasks.irp.f index 7bb6e16e..f23d51e9 100644 --- a/src/Bitmask/bitmasks.irp.f +++ b/src/Bitmask/bitmasks.irp.f @@ -37,6 +37,30 @@ BEGIN_PROVIDER [ integer(bit_kind), full_ijkl_bitmask_4, (N_int,4) ] enddo END_PROVIDER +BEGIN_PROVIDER [ integer(bit_kind), core_inact_act_bitmask_4, (N_int,4) ] + implicit none + integer :: i + do i=1,N_int + core_inact_act_bitmask_4(i,1) = reunion_of_core_inact_act_bitmask(i,1) + core_inact_act_bitmask_4(i,2) = reunion_of_core_inact_act_bitmask(i,1) + core_inact_act_bitmask_4(i,3) = reunion_of_core_inact_act_bitmask(i,1) + core_inact_act_bitmask_4(i,4) = reunion_of_core_inact_act_bitmask(i,1) + enddo +END_PROVIDER + +BEGIN_PROVIDER [ integer(bit_kind), virt_bitmask_4, (N_int,4) ] + implicit none + integer :: i + do i=1,N_int + virt_bitmask_4(i,1) = virt_bitmask(i,1) + virt_bitmask_4(i,2) = virt_bitmask(i,1) + virt_bitmask_4(i,3) = virt_bitmask(i,1) + virt_bitmask_4(i,4) = virt_bitmask(i,1) + enddo +END_PROVIDER + + + BEGIN_PROVIDER [ integer(bit_kind), HF_bitmask, (N_int,2)] implicit none @@ -369,11 +393,19 @@ END_PROVIDER BEGIN_PROVIDER [ integer, list_inact, (n_inact_orb)] &BEGIN_PROVIDER [ integer, list_virt, (n_virt_orb)] + &BEGIN_PROVIDER [ integer, list_inact_reverse, (mo_tot_num)] + &BEGIN_PROVIDER [ integer, list_virt_reverse, (mo_tot_num)] BEGIN_DOC ! list_inact : List of the inactive orbitals which are supposed to be doubly excited ! in post CAS methods ! list_virt : List of vritual orbitals which are supposed to be recieve electrons ! in post CAS methods + ! list_inact_reverse : reverse list of inactive orbitals + ! list_inact_reverse(i) = 0 ::> not an inactive + ! list_inact_reverse(i) = k ::> IS the kth inactive + ! list_virt_reverse : reverse list of virtual orbitals + ! list_virt_reverse(i) = 0 ::> not an virtual + ! list_virt_reverse(i) = k ::> IS the kth virtual END_DOC implicit none integer :: occ_inact(N_int*bit_kind_size) @@ -381,15 +413,20 @@ END_PROVIDER occ_inact = 0 call bitstring_to_list(inact_bitmask(1,1), occ_inact(1), itest, N_int) ASSERT(itest==n_inact_orb) + list_inact_reverse = 0 do i = 1, n_inact_orb list_inact(i) = occ_inact(i) + list_inact_reverse(occ_inact(i)) = i enddo + occ_inact = 0 call bitstring_to_list(virt_bitmask(1,1), occ_inact(1), itest, N_int) ASSERT(itest==n_virt_orb) + list_virt_reverse = 0 do i = 1, n_virt_orb list_virt(i) = occ_inact(i) + list_virt_reverse(occ_inact(i)) = i enddo END_PROVIDER @@ -397,7 +434,7 @@ END_PROVIDER BEGIN_PROVIDER [ integer(bit_kind), reunion_of_core_inact_bitmask, (N_int,2)] implicit none BEGIN_DOC - ! Reunion of the inactive, active and virtual bitmasks + ! Reunion of the core and inactive and virtual bitmasks END_DOC integer :: i,j do i = 1, N_int @@ -407,6 +444,20 @@ END_PROVIDER END_PROVIDER + BEGIN_PROVIDER [integer(bit_kind), reunion_of_core_inact_act_bitmask, (N_int,2)] + implicit none + BEGIN_DOC + ! Reunion of the core, inactive and active bitmasks + END_DOC + integer :: i,j + + do i = 1, N_int + reunion_of_core_inact_act_bitmask(i,1) = ior(reunion_of_core_inact_bitmask(i,1),cas_bitmask(i,1,1)) + reunion_of_core_inact_act_bitmask(i,2) = ior(reunion_of_core_inact_bitmask(i,2),cas_bitmask(i,1,1)) + enddo + END_PROVIDER + + BEGIN_PROVIDER [ integer(bit_kind), reunion_of_bitmask, (N_int,2)] @@ -435,6 +486,7 @@ END_PROVIDER END_PROVIDER BEGIN_PROVIDER [integer, list_core, (n_core_orb)] +&BEGIN_PROVIDER [integer, list_core_reverse, (mo_tot_num)] BEGIN_DOC ! List of the core orbitals that are never excited in post CAS method END_DOC @@ -444,8 +496,10 @@ END_PROVIDER occ_core = 0 call bitstring_to_list(core_bitmask(1,1), occ_core(1), itest, N_int) ASSERT(itest==n_core_orb) + list_core_reverse = 0 do i = 1, n_core_orb list_core(i) = occ_core(i) + list_core_reverse(occ_core(i)) = i enddo END_PROVIDER @@ -497,11 +551,17 @@ BEGIN_PROVIDER [ integer, n_act_orb] do i = 1, N_int n_act_orb += popcnt(cas_bitmask(i,1,1)) enddo + print*,'n_act_orb = ',n_act_orb END_PROVIDER -BEGIN_PROVIDER [integer, list_act, (n_act_orb)] + BEGIN_PROVIDER [integer, list_act, (n_act_orb)] +&BEGIN_PROVIDER [integer, list_act_reverse, (mo_tot_num)] BEGIN_DOC - ! list of active orbitals + ! list_act(i) = index of the ith active orbital + ! + ! list_act_reverse : reverse list of active orbitals + ! list_act_reverse(i) = 0 ::> not an active + ! list_act_reverse(i) = k ::> IS the kth active orbital END_DOC implicit none integer :: occ_act(N_int*bit_kind_size) @@ -509,10 +569,11 @@ BEGIN_PROVIDER [integer, list_act, (n_act_orb)] occ_act = 0 call bitstring_to_list(cas_bitmask(1,1,1), occ_act(1), itest, N_int) ASSERT(itest==n_act_orb) + list_act_reverse = 0 do i = 1, n_act_orb list_act(i) = occ_act(i) + list_act_reverse(occ_act(i)) = i enddo - END_PROVIDER BEGIN_PROVIDER [integer(bit_kind), closed_shell_ref_bitmask, (N_int,2)] @@ -537,4 +598,19 @@ END_PROVIDER enddo END_PROVIDER + + BEGIN_PROVIDER [integer, n_core_orb_allocate] + implicit none + n_core_orb_allocate = max(n_core_orb,1) + END_PROVIDER + + BEGIN_PROVIDER [integer, n_inact_orb_allocate] + implicit none + n_inact_orb_allocate = max(n_inact_orb,1) + END_PROVIDER + + BEGIN_PROVIDER [integer, n_virt_orb_allocate] + implicit none + n_virt_orb_allocate = max(n_virt_orb,1) + END_PROVIDER diff --git a/src/Determinants/density_matrix.irp.f b/src/Determinants/density_matrix.irp.f index 62d09381..2253c33c 100644 --- a/src/Determinants/density_matrix.irp.f +++ b/src/Determinants/density_matrix.irp.f @@ -238,17 +238,19 @@ END_PROVIDER END_DOC implicit none integer :: i,j,k,l - double precision :: dm_mo + double precision :: mo_alpha,mo_beta - one_body_spin_density_ao = 0.d0 + one_body_dm_ao_alpha = 0.d0 + one_body_dm_ao_beta = 0.d0 do k = 1, ao_num do l = 1, ao_num do i = 1, mo_tot_num do j = 1, mo_tot_num - dm_mo = one_body_dm_mo_alpha(j,i) + mo_alpha = one_body_dm_mo_alpha(j,i) + mo_beta = one_body_dm_mo_beta(j,i) ! if(dabs(dm_mo).le.1.d-10)cycle - one_body_dm_ao_alpha(l,k) += mo_coef(k,i) * mo_coef(l,j) * dm_mo - one_body_dm_ao_beta(l,k) += mo_coef(k,i) * mo_coef(l,j) * dm_mo + one_body_dm_ao_alpha(l,k) += mo_coef(k,i) * mo_coef(l,j) * mo_alpha + one_body_dm_ao_beta(l,k) += mo_coef(k,i) * mo_coef(l,j) * mo_beta enddo enddo diff --git a/src/Determinants/diagonalize_restart_and_save_all_states.irp.f b/src/Determinants/diagonalize_restart_and_save_all_states.irp.f new file mode 100644 index 00000000..3bdc37c5 --- /dev/null +++ b/src/Determinants/diagonalize_restart_and_save_all_states.irp.f @@ -0,0 +1,16 @@ +program diag_and_save + implicit none + read_wf = .True. + touch read_wf + call routine +end + +subroutine routine + implicit none + call diagonalize_CI + print*,'N_det = ',N_det + call save_wavefunction_general(N_det,N_states_diag,psi_det_sorted,size(psi_coef_sorted,1),psi_coef_sorted) + + + +end diff --git a/src/Determinants/diagonalize_restart_and_save_lowest_state.irp.f b/src/Determinants/diagonalize_restart_and_save_lowest_state.irp.f new file mode 100644 index 00000000..11c98034 --- /dev/null +++ b/src/Determinants/diagonalize_restart_and_save_lowest_state.irp.f @@ -0,0 +1,25 @@ +program diag_and_save + implicit none + read_wf = .True. + touch read_wf + call routine +end + +subroutine routine + implicit none + print*,'N_det = ',N_det + call diagonalize_CI + integer :: igood_state + igood_state=1 + double precision, allocatable :: psi_coef_tmp(:) + allocate(psi_coef_tmp(n_det)) + integer :: i + do i = 1, N_det + psi_coef_tmp(i) = psi_coef(i,igood_state) + enddo + call save_wavefunction_general(N_det,1,psi_det,n_det,psi_coef_tmp) + deallocate(psi_coef_tmp) + + + +end diff --git a/src/Determinants/diagonalize_restart_and_save_one_states.irp.f b/src/Determinants/diagonalize_restart_and_save_one_states.irp.f new file mode 100644 index 00000000..c5f4e59d --- /dev/null +++ b/src/Determinants/diagonalize_restart_and_save_one_states.irp.f @@ -0,0 +1,26 @@ +program diag_and_save + implicit none + read_wf = .True. + touch read_wf + call routine +end + +subroutine routine + implicit none + print*,'N_det = ',N_det + call diagonalize_CI + write(*,*)'Which state would you like to save ?' + integer :: igood_state + read(5,*)igood_state + double precision, allocatable :: psi_coef_tmp(:) + allocate(psi_coef_tmp(n_det)) + integer :: i + do i = 1, N_det + psi_coef_tmp(i) = psi_coef(i,igood_state) + enddo + call save_wavefunction_general(N_det,1,psi_det,n_det,psi_coef_tmp) + deallocate(psi_coef_tmp) + + + +end diff --git a/src/Determinants/diagonalize_restart_and_save_two_states.irp.f b/src/Determinants/diagonalize_restart_and_save_two_states.irp.f new file mode 100644 index 00000000..97fed531 --- /dev/null +++ b/src/Determinants/diagonalize_restart_and_save_two_states.irp.f @@ -0,0 +1,27 @@ +program diag_and_save + implicit none + read_wf = .True. + touch read_wf + call routine +end + +subroutine routine + implicit none + integer :: igood_state_1,igood_state_2 + double precision, allocatable :: psi_coef_tmp(:,:) + integer :: i + print*,'N_det = ',N_det +!call diagonalize_CI + write(*,*)'Which couple of states would you like to save ?' + read(5,*)igood_state_1,igood_state_2 + allocate(psi_coef_tmp(n_det,2)) + do i = 1, N_det + psi_coef_tmp(i,1) = psi_coef(i,igood_state_1) + psi_coef_tmp(i,2) = psi_coef(i,igood_state_2) + enddo + call save_wavefunction_general(N_det,2,psi_det,n_det,psi_coef_tmp) + deallocate(psi_coef_tmp) + + + +end diff --git a/src/Determinants/occ_pattern.irp.f b/src/Determinants/occ_pattern.irp.f index aa059870..e2e12974 100644 --- a/src/Determinants/occ_pattern.irp.f +++ b/src/Determinants/occ_pattern.irp.f @@ -256,7 +256,7 @@ subroutine make_s2_eigenfunction integer :: N_det_new integer, parameter :: bufsze = 1000 logical, external :: is_in_wavefunction - return +! return ! !TODO DEBUG ! do i=1,N_det diff --git a/src/Determinants/print_H_matrix_restart.irp.f b/src/Determinants/print_H_matrix_restart.irp.f new file mode 100644 index 00000000..813f14d0 --- /dev/null +++ b/src/Determinants/print_H_matrix_restart.irp.f @@ -0,0 +1,179 @@ +program print_H_matrix_restart + implicit none + read_wf = .True. + touch read_wf + call routine + +end + +subroutine routine + use bitmasks + implicit none + integer :: i,j + integer, allocatable :: H_matrix_degree(:,:) + double precision, allocatable :: H_matrix_phase(:,:) + integer :: degree + integer(bit_kind), allocatable :: keys_tmp(:,:,:) + allocate(keys_tmp(N_int,2,N_det)) + do i = 1, N_det + print*,'' + call debug_det(psi_det(1,1,i),N_int) + do j = 1, N_int + keys_tmp(j,1,i) = psi_det(j,1,i) + keys_tmp(j,2,i) = psi_det(j,2,i) + enddo + enddo + if(N_det.ge.10000)then + print*,'Warning !!!' + print*,'Number of determinants is ',N_det + print*,'It means that the H matrix will be enormous !' + print*,'stoppping ..' + stop + endif + print*,'' + print*,'Determinants ' + do i = 1, N_det + enddo + allocate(H_matrix_degree(N_det,N_det),H_matrix_phase(N_det,N_det)) + integer :: exc(0:2,2,2) + double precision :: phase + do i = 1, N_det + do j = i, N_det + call get_excitation_degree(psi_det(1,1,i),psi_det(1,1,j),degree,N_int) + H_matrix_degree(i,j) = degree + H_matrix_degree(j,i) = degree + phase = 0.d0 + if(degree==1.or.degree==2)then + call get_excitation(psi_det(1,1,i),psi_det(1,1,j),exc,degree,phase,N_int) + endif + H_matrix_phase(i,j) = phase + H_matrix_phase(j,i) = phase + enddo + enddo + print*,'H matrix ' + double precision :: ref_h_matrix,s2 + ref_h_matrix = H_matrix_all_dets(1,1) + print*,'HF like determinant energy = ',ref_bitmask_energy+nuclear_repulsion + print*,'Ref element of H_matrix = ',ref_h_matrix+nuclear_repulsion + print*,'Printing the H matrix ...' + print*,'' + print*,'' +!do i = 1, N_det +! H_matrix_all_dets(i,i) -= ref_h_matrix +!enddo + + do i = 1, N_det + H_matrix_all_dets(i,i) += nuclear_repulsion + enddo + +!do i = 5,N_det +! H_matrix_all_dets(i,3) = 0.d0 +! H_matrix_all_dets(3,i) = 0.d0 +! H_matrix_all_dets(i,4) = 0.d0 +! H_matrix_all_dets(4,i) = 0.d0 +!enddo + + + + + + do i = 1, N_det + write(*,'(I3,X,A3,1000(F16.7))')i,' | ',H_matrix_all_dets(i,:) + enddo + + print*,'' + print*,'' + print*,'' + print*,'Printing the degree of excitations within the H matrix' + print*,'' + print*,'' + do i = 1, N_det + write(*,'(I3,X,A3,X,1000(I1,X))')i,' | ',H_matrix_degree(i,:) + enddo + + + print*,'' + print*,'' + print*,'Printing the phase of the Hamiltonian matrix elements ' + print*,'' + print*,'' + do i = 1, N_det + write(*,'(I3,X,A3,X,1000(F3.0,X))')i,' | ',H_matrix_phase(i,:) + enddo + print*,'' + + + double precision, allocatable :: eigenvectors(:,:), eigenvalues(:) + double precision, allocatable :: s2_eigvalues(:) + allocate (eigenvectors(size(H_matrix_all_dets,1),N_det)) + allocate (eigenvalues(N_det)) + call lapack_diag(eigenvalues,eigenvectors, & + H_matrix_all_dets,size(H_matrix_all_dets,1),N_det) + print*,'Two first eigenvectors ' + do j = 1, n_states +!do j = 1, 1 + print*,'State ',j + call get_s2_u0(keys_tmp,eigenvectors(1,j),N_det,size(eigenvectors,1),s2) + print*,'s2 = ',s2 + print*,'e = ',eigenvalues(j) + print*,'coefs : ' + do i = 1, N_det + print*,'i = ',i,eigenvectors(i,j) + enddo + + if(j>1)then + print*,'Delta E(H) = ',eigenvalues(1) - eigenvalues(j) + print*,'Delta E(eV) = ',(eigenvalues(1) - eigenvalues(j))*27.2114d0 + endif + enddo + double precision :: get_mo_bielec_integral_schwartz,k_a_iv,k_b_iv + integer :: h1,p1,h2,p2 + h1 = 10 + p1 = 16 + h2 = 14 + p2 = 14 +!h1 = 1 +!p1 = 4 +!h2 = 2 +!p2 = 2 + k_a_iv = get_mo_bielec_integral_schwartz(h1,h2,p2,p1,mo_integrals_map) + h2 = 15 + p2 = 15 + k_b_iv = get_mo_bielec_integral_schwartz(h1,h2,p2,p1,mo_integrals_map) + print*,'k_a_iv = ',k_a_iv + print*,'k_b_iv = ',k_b_iv + double precision :: k_av,k_bv,k_ai,k_bi + h1 = 16 + p1 = 14 + h2 = 14 + p2 = 16 + k_av = get_mo_bielec_integral_schwartz(h1,h2,p1,p2,mo_integrals_map) + h1 = 16 + p1 = 15 + h2 = 15 + p2 = 16 + k_bv = get_mo_bielec_integral_schwartz(h1,h2,p1,p2,mo_integrals_map) + + h1 = 10 + p1 = 14 + h2 = 14 + p2 = 10 + k_ai = get_mo_bielec_integral_schwartz(h1,h2,p1,p2,mo_integrals_map) + + h1 = 10 + p1 = 15 + h2 = 15 + p2 = 10 + k_bi = get_mo_bielec_integral_schwartz(h1,h2,p1,p2,mo_integrals_map) + + print*,'k_av, k_bv = ',k_av,k_bv + print*,'k_ai, k_bi = ',k_ai,k_bi + double precision :: k_iv + + h1 = 10 + p1 = 16 + h2 = 16 + p2 = 10 + k_iv = get_mo_bielec_integral_schwartz(h1,h2,p1,p2,mo_integrals_map) + print*,'k_iv = ',k_iv +end diff --git a/src/Determinants/print_bitmask.irp.f b/src/Determinants/print_bitmask.irp.f new file mode 100644 index 00000000..2f1c8f73 --- /dev/null +++ b/src/Determinants/print_bitmask.irp.f @@ -0,0 +1,11 @@ +program print_bitmask + implicit none + print*,'core' + call debug_det(core_bitmask,N_int) + print*,'inact' + call debug_det(inact_bitmask,N_int) + print*,'virt' + call debug_det(virt_bitmask,N_int) + + +end diff --git a/src/Determinants/print_holes_particles.irp.f b/src/Determinants/print_holes_particles.irp.f new file mode 100644 index 00000000..601015f7 --- /dev/null +++ b/src/Determinants/print_holes_particles.irp.f @@ -0,0 +1,36 @@ +program pouet + implicit none + read_wf = .True. + touch read_wf + call routine +end + +subroutine routine + implicit none + integer :: i,j,number_of_holes,number_of_particles + integer :: n_h,n_p + do i = 1, N_det + n_h = number_of_holes(psi_det(1,1,i)) + n_p = number_of_particles(psi_det(1,1,i)) + if(n_h == 0 .and. n_p == 0)then + print*,'CAS' + else if(n_h == 1 .and. n_p ==0)then + print*,'1h' + else if(n_h == 0 .and. n_p ==1)then + print*,'1p' + else if(n_h == 1 .and. n_p ==1)then + print*,'1h1p' + else if(n_h == 2 .and. n_p ==1)then + print*,'2h1p' + else if(n_h == 1 .and. n_p ==2)then + print*,'1h2p' + else + print*,'PB !! ' + call debug_det(psi_det(1,1,i), N_int) + stop + endif + enddo + + + +end diff --git a/src/Determinants/print_wf.irp.f b/src/Determinants/print_wf.irp.f new file mode 100644 index 00000000..e7ca5dc6 --- /dev/null +++ b/src/Determinants/print_wf.irp.f @@ -0,0 +1,71 @@ +program printwf + implicit none + read_wf = .True. + touch read_wf + print*,'ref_bitmask_energy = ',ref_bitmask_energy + call routine + +end + +subroutine routine + implicit none + integer :: i + integer :: degree + double precision :: hij + integer :: exc(0:2,2,2) + double precision :: phase + integer :: h1,p1,h2,p2,s1,s2 + double precision :: get_mo_bielec_integral_schwartz + double precision :: norm_mono_a,norm_mono_b + norm_mono_a = 0.d0 + norm_mono_b = 0.d0 + do i = 1, min(500,N_det) + print*,'' + print*,'i = ',i + call debug_det(psi_det(1,1,i),N_int) + call get_excitation_degree(psi_det(1,1,i),psi_det(1,1,1),degree,N_int) + print*,'degree = ',degree + if(degree == 0)then + print*,'Reference determinant ' + else + call i_H_j(psi_det(1,1,i),psi_det(1,1,1),N_int,hij) + call get_excitation(psi_det(1,1,1),psi_det(1,1,i),exc,degree,phase,N_int) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + print*,'phase = ',phase + if(degree == 1)then + print*,'s1',s1 + print*,'h1,p1 = ',h1,p1 + if(s1 == 1)then + norm_mono_a += dabs(psi_coef(i,1)/psi_coef(1,1)) + else + norm_mono_b += dabs(psi_coef(i,1)/psi_coef(1,1)) + endif + print*,'< h | Ka| p > = ',get_mo_bielec_integral_schwartz(h1,list_act(1),list_act(1),p1,mo_integrals_map) + double precision :: hmono,hdouble + call i_H_j_verbose(psi_det(1,1,1),psi_det(1,1,i),N_int,hij,hmono,hdouble) + print*,'hmono = ',hmono + print*,'hdouble = ',hdouble + print*,'hmono+hdouble = ',hmono+hdouble + print*,'hij = ',hij + else + print*,'s1',s1 + print*,'h1,p1 = ',h1,p1 + print*,'s2',s2 + print*,'h2,p2 = ',h2,p2 + print*,'< h | Ka| p > = ',get_mo_bielec_integral_schwartz(h1,h2,p1,p2,mo_integrals_map) + endif + + print*,' = ',hij + endif + print*,'amplitude = ',psi_coef(i,1)/psi_coef(1,1) + + enddo + + + print*,'' + print*,'' + print*,'' + print*,'mono alpha = ',norm_mono_a + print*,'mono beta = ',norm_mono_b + +end diff --git a/src/Determinants/save_only_singles.irp.f b/src/Determinants/save_only_singles.irp.f new file mode 100644 index 00000000..ae68a52c --- /dev/null +++ b/src/Determinants/save_only_singles.irp.f @@ -0,0 +1,50 @@ +program save_only_singles + implicit none + read_wf = .True. + touch read_wf + call routine +end + +subroutine routine + implicit none + integer :: i,j,k,l + use bitmasks + integer :: n_det_restart,degree + integer(bit_kind),allocatable :: psi_det_tmp(:,:,:) + double precision ,allocatable :: psi_coef_tmp(:,:),accu(:) + integer, allocatable :: index_restart(:) + allocate(index_restart(N_det)) + N_det_restart = 0 + do i = 1, N_det + call get_excitation_degree(psi_det(1,1,1),psi_det(1,1,i),degree,N_int) + if(degree == 0 .or. degree==1)then + N_det_restart +=1 + index_restart(N_det_restart) = i + cycle + endif + enddo + allocate (psi_det_tmp(N_int,2,N_det_restart),psi_coef_tmp(N_det_restart,N_states),accu(N_states)) + accu = 0.d0 + do i = 1, N_det_restart + do j = 1, N_int + psi_det_tmp(j,1,i) = psi_det(j,1,index_restart(i)) + psi_det_tmp(j,2,i) = psi_det(j,2,index_restart(i)) + enddo + do j = 1,N_states + psi_coef_tmp(i,j) = psi_coef(index_restart(i),j) + accu(j) += psi_coef_tmp(i,j) * psi_coef_tmp(i,j) + enddo + enddo + do j = 1, N_states + accu(j) = 1.d0/dsqrt(accu(j)) + enddo + do j = 1,N_states + do i = 1, N_det_restart + psi_coef_tmp(i,j) = psi_coef_tmp(i,j) * accu(j) + enddo + enddo + call save_wavefunction_general(N_det_restart,N_states,psi_det_tmp,N_det_restart,psi_coef_tmp) + + deallocate (psi_det_tmp,psi_coef_tmp,accu,index_restart) + +end diff --git a/src/Determinants/test_3d.irp.f b/src/Determinants/test_3d.irp.f new file mode 100644 index 00000000..748890da --- /dev/null +++ b/src/Determinants/test_3d.irp.f @@ -0,0 +1,15 @@ +program test_3d + implicit none + integer :: i,npt + double precision :: dx,domain,x_min,x,step_function_becke + domain = 5.d0 + npt = 100 + dx = domain/dble(npt) + x_min = -0.5d0 * domain + x = x_min + do i = 1, npt + write(33,*)x,step_function_becke(x) + x += dx + enddo + +end diff --git a/src/Determinants/test_two_body.irp.f b/src/Determinants/test_two_body.irp.f new file mode 100644 index 00000000..54c43c09 --- /dev/null +++ b/src/Determinants/test_two_body.irp.f @@ -0,0 +1,18 @@ +program test + implicit none + read_wf = .True. + touch read_wf + call routine +end + +subroutine routine + implicit none + integer :: i,j,k,l + do i = 1, n_act_orb + do j = 1, n_act_orb + do k = 1, n_act_orb + + enddo + enddo + enddo +end diff --git a/src/Determinants/truncate_wf.irp.f b/src/Determinants/truncate_wf.irp.f index 42340c71..aba16fa7 100644 --- a/src/Determinants/truncate_wf.irp.f +++ b/src/Determinants/truncate_wf.irp.f @@ -1,18 +1,11 @@ -program cisd - implicit none - integer :: i,k - - - double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:) - integer :: N_st, degree - N_det=10000 - do i=1,N_det - do k=1,N_int - psi_det(k,1,i) = psi_det_sorted(k,1,i) - psi_det(k,2,i) = psi_det_sorted(k,2,i) - enddo - psi_coef(i,:) = psi_coef_sorted(i,:) - enddo +program s2_eig_restart + implicit none + read_wf = .True. + call routine +end +subroutine routine + implicit none + call make_s2_eigenfunction TOUCH psi_det psi_coef psi_det_sorted psi_coef_sorted psi_average_norm_contrib_sorted N_det call save_wavefunction end diff --git a/src/Determinants/two_body_dm_map.irp.f b/src/Determinants/two_body_dm_map.irp.f index f88d6ea3..aa8f630b 100644 --- a/src/Determinants/two_body_dm_map.irp.f +++ b/src/Determinants/two_body_dm_map.irp.f @@ -193,33 +193,141 @@ subroutine add_values_to_two_body_dm_map(mask_ijkl) end -BEGIN_PROVIDER [double precision, two_body_dm_ab_diag, (mo_tot_num, mo_tot_num)] + BEGIN_PROVIDER [double precision, two_body_dm_ab_diag_act, (n_act_orb, n_act_orb)] +&BEGIN_PROVIDER [double precision, two_body_dm_ab_diag_inact, (n_inact_orb_allocate, n_inact_orb_allocate)] +&BEGIN_PROVIDER [double precision, two_body_dm_ab_diag_core, (n_core_orb_allocate, n_core_orb_allocate)] +&BEGIN_PROVIDER [double precision, two_body_dm_ab_diag_all, (mo_tot_num, mo_tot_num)] +&BEGIN_PROVIDER [double precision, two_body_dm_diag_core_a_act_b, (n_core_orb_allocate,n_act_orb)] +&BEGIN_PROVIDER [double precision, two_body_dm_diag_core_b_act_a, (n_core_orb_allocate,n_act_orb)] +&BEGIN_PROVIDER [double precision, two_body_dm_diag_core_act, (n_core_orb_allocate,n_act_orb)] implicit none + use bitmasks integer :: i,j,k,l,m integer :: occ(N_int*bit_kind_size,2) integer :: n_occ_ab(2) + integer :: occ_act(N_int*bit_kind_size,2) + integer :: n_occ_ab_act(2) + integer :: occ_core(N_int*bit_kind_size,2) + integer :: n_occ_ab_core(2) double precision :: contrib BEGIN_DOC - ! two_body_dm_ab_diag(k,m) = <\Psi | n_(k\alpha) n_(m\beta) | \Psi> + ! two_body_dm_ab_diag_all(k,m) = <\Psi | n_(k\alpha) n_(m\beta) | \Psi> + ! two_body_dm_ab_diag_act(k,m) is restricted to the active orbitals : + ! orbital k = list_act(k) + ! two_body_dm_ab_diag_inact(k,m) is restricted to the inactive orbitals : + ! orbital k = list_inact(k) + ! two_body_dm_ab_diag_core(k,m) is restricted to the core orbitals : + ! orbital k = list_core(k) + ! two_body_dm_ab_diag_core_b_act_a(k,m) represents the core beta <-> active alpha part of the two body dm + ! orbital k = list_core(k) + ! orbital m = list_act(m) + ! two_body_dm_ab_diag_core_a_act_b(k,m) represents the core alpha <-> active beta part of the two body dm + ! orbital k = list_core(k) + ! orbital m = list_act(m) + ! two_body_dm_ab_diag_core_act(k,m) represents the core<->active part of the diagonal two body dm + ! when we traced on the spin + ! orbital k = list_core(k) + ! orbital m = list_act(m) END_DOC + integer(bit_kind) :: key_tmp_core(N_int,2) + integer(bit_kind) :: key_tmp_act(N_int,2) - two_body_dm_ab_diag = 0.d0 + two_body_dm_ab_diag_all = 0.d0 + two_body_dm_ab_diag_act = 0.d0 + two_body_dm_ab_diag_core = 0.d0 + two_body_dm_ab_diag_inact = 0.d0 + two_body_dm_diag_core_a_act_b = 0.d0 + two_body_dm_diag_core_b_act_a = 0.d0 + two_body_dm_diag_core_act = 0.d0 do i = 1, N_det ! i == |I> - call bitstring_to_list_ab(psi_det(1,1,i), occ, n_occ_ab, N_int) + ! Full diagonal part of the two body dm contrib = psi_coef(i,1)**2 + call bitstring_to_list_ab(psi_det(1,1,i), occ, n_occ_ab, N_int) do j = 1, elec_beta_num k = occ(j,2) do l = 1, elec_alpha_num m = occ(l,1) - two_body_dm_ab_diag(k,m) += 0.5d0 * contrib - two_body_dm_ab_diag(m,k) += 0.5d0 * contrib + two_body_dm_ab_diag_all(k,m) += 0.5d0 * contrib + two_body_dm_ab_diag_all(m,k) += 0.5d0 * contrib enddo enddo + + ! ACTIVE PART of the diagonal part of the two body dm + do j = 1, N_int + key_tmp_act(j,1) = psi_det(j,1,i) + key_tmp_act(j,2) = psi_det(j,2,i) + enddo + do j = 1, N_int + key_tmp_act(j,1) = iand(key_tmp_act(j,1),cas_bitmask(j,1,1)) + key_tmp_act(j,2) = iand(key_tmp_act(j,2),cas_bitmask(j,1,1)) + enddo + call bitstring_to_list_ab(key_tmp_act, occ_act, n_occ_ab_act, N_int) + do j = 1,n_occ_ab_act(2) + k = list_act_reverse(occ_act(j,2)) + do l = 1, n_occ_ab_act(1) + m = list_act_reverse(occ_act(l,1)) + two_body_dm_ab_diag_act(k,m) += 0.5d0 * contrib + two_body_dm_ab_diag_act(m,k) += 0.5d0 * contrib + enddo + enddo + + ! CORE PART of the diagonal part of the two body dm + do j = 1, N_int + key_tmp_core(j,1) = psi_det(j,1,i) + key_tmp_core(j,2) = psi_det(j,2,i) + enddo + do j = 1, N_int + key_tmp_core(j,1) = iand(key_tmp_core(j,1),core_bitmask(j,1)) + key_tmp_core(j,2) = iand(key_tmp_core(j,2),core_bitmask(j,1)) + enddo + call bitstring_to_list_ab(key_tmp_core, occ_core, n_occ_ab_core, N_int) + do j = 1,n_occ_ab_core(2) + k = list_core_reverse(occ_core(j,2)) + do l = 1, n_occ_ab_core(1) + m = list_core_reverse(occ_core(l,1)) + two_body_dm_ab_diag_core(k,m) += 0.5d0 * contrib + two_body_dm_ab_diag_core(m,k) += 0.5d0 * contrib + enddo + enddo + + ! ACT<->CORE PART + ! alpha electron in active space + do j = 1,n_occ_ab_act(1) + k = list_act_reverse(occ_act(j,1)) + ! beta electron in core space + do l = 1, n_occ_ab_core(2) + m = list_core_reverse(occ_core(l,2)) + ! The fact that you have 1 * contrib and not 0.5 * contrib + ! takes into account the following symmetry : + ! 0.5 * + 0.5 * + two_body_dm_diag_core_b_act_a(m,k) += contrib + enddo + enddo + ! beta electron in active space + do j = 1,n_occ_ab_act(2) + k = list_act_reverse(occ_act(j,2)) + ! alpha electron in core space + do l = 1, n_occ_ab_core(1) + m = list_core_reverse(occ_core(l,1)) + ! The fact that you have 1 * contrib and not 0.5 * contrib + ! takes into account the following symmetry : + ! 0.5 * + 0.5 * + two_body_dm_diag_core_a_act_b(m,k) += contrib + enddo + enddo + enddo + + do j = 1, n_core_orb + do l = 1, n_act_orb + two_body_dm_diag_core_act(j,l) = two_body_dm_diag_core_b_act_a(j,l) + two_body_dm_diag_core_a_act_b(j,l) + enddo enddo END_PROVIDER -BEGIN_PROVIDER [double precision, two_body_dm_ab_big_array, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)] + BEGIN_PROVIDER [double precision, two_body_dm_ab_big_array_act, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)] +&BEGIN_PROVIDER [double precision, two_body_dm_ab_big_array_core_act, (n_core_orb_allocate,n_act_orb,n_act_orb)] implicit none + use bitmasks integer :: i,j,k,l,m integer :: degree PROVIDE mo_coef psi_coef psi_det @@ -229,54 +337,108 @@ BEGIN_PROVIDER [double precision, two_body_dm_ab_big_array, (n_act_orb,n_act_orb double precision :: contrib integer :: occ(N_int*bit_kind_size,2) integer :: n_occ_ab(2) - two_body_dm_ab_big_array = 0.d0 + integer :: occ_core(N_int*bit_kind_size,2) + integer :: n_occ_ab_core(2) + integer(bit_kind) :: key_tmp_i(N_int,2) + integer(bit_kind) :: key_tmp_i_core(N_int,2) + integer(bit_kind) :: key_tmp_j(N_int,2) + two_body_dm_ab_big_array_act = 0.d0 + two_body_dm_ab_big_array_core_act = 0.d0 BEGIN_DOC -! The alpha-beta energy can be computed thanks to -! sum_{h1,p1,h2,p2} two_body_dm_ab_big_array(h1,p1,h2,p2) * (h1p1|h2p2) +! two_body_dm_ab_big_array_act = Purely active part of the two body density matrix +! two_body_dm_ab_big_array_act_core takes only into account the single excitation +! within the active space that adds terms in the act <-> core two body dm +! two_body_dm_ab_big_array_act_core(i,j,k) = < a^\dagger_i n_k a_j > +! with i,j in the ACTIVE SPACE +! with k in the CORE SPACE +! +! The alpha-beta extra diagonal energy FOR WF DEFINED AS AN APPROXIMATION OF A CAS can be computed thanks to +! sum_{h1,p1,h2,p2} two_body_dm_ab_big_array_act(h1,p1,h2,p2) * (h1p1|h2p2) +! + sum_{h1,p1,h2,p2} two_body_dm_ab_big_array_core_act(h1,p1,h2,p2) * (h1p1|h2p2) END_DOC do i = 1, N_det ! i == |I> - call bitstring_to_list_ab(psi_det(1,1,i), occ, n_occ_ab, N_int) + ! active part of psi_det(i) + do j = 1, N_int + key_tmp_i(j,1) = psi_det(j,1,i) + key_tmp_i(j,2) = psi_det(j,2,i) + key_tmp_i_core(j,1) = psi_det(j,1,i) + key_tmp_i_core(j,2) = psi_det(j,2,i) + enddo + do j = 1, N_int + key_tmp_i(j,1) = iand(key_tmp_i(j,1),cas_bitmask(j,1,1)) + key_tmp_i(j,2) = iand(key_tmp_i(j,2),cas_bitmask(j,1,1)) + enddo + do j = 1, N_int + key_tmp_i_core(j,1) = iand(key_tmp_i_core(j,1),core_bitmask(j,1)) + key_tmp_i_core(j,2) = iand(key_tmp_i_core(j,2),core_bitmask(j,1)) + enddo + call bitstring_to_list_ab(key_tmp_i_core, occ_core, n_occ_ab_core, N_int) + call bitstring_to_list_ab(key_tmp_i, occ, n_occ_ab, N_int) do j = i+1, N_det ! j == 2)cycle + ! if it is the case, then compute the hamiltonian matrix element with the proper phase call get_excitation(psi_det(1,1,i),psi_det(1,1,j),exc,degree,phase,N_int) call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) -! if(i==3.or.j==3)then -! print*,'i,j = ',i,j -! call debug_det(psi_det(1,1,i),N_int) -! call debug_det(psi_det(1,1,j),N_int) -! print*,degree,s1,s2 -! print*,h1,p1,h2,p2 -! print*,phase -! pause -! endif contrib = 0.5d0 * psi_coef(i,1) * psi_coef(j,1) * phase -! print*,'coucou' -! print*,'i,j = ',i,j -! print*,'contrib = ',contrib -! print*,h1,p1,h2,p2 -! print*,'s1,s2',s1,s2 -! call debug_det(psi_det(1,1,i),N_int) -! call debug_det(psi_det(1,1,j),N_int) -! pause if(degree==2)then ! case of the DOUBLE EXCITATIONS ************************************ if(s1==s2)cycle ! Only the alpha/beta two body density matrix ! * c_I * c_J - call insert_into_two_body_dm_big_array( two_body_dm_ab_big_array,n_act_orb,n_act_orb,n_act_orb,n_act_orb,contrib,h1,p1,h2,p2) + h1 = list_act_reverse(h1) + h2 = list_act_reverse(h2) + p1 = list_act_reverse(p1) + p2 = list_act_reverse(p2) + call insert_into_two_body_dm_big_array( two_body_dm_ab_big_array_act,n_act_orb,n_act_orb,n_act_orb,n_act_orb,contrib,h1,p1,h2,p2) else if(degree==1)then! case of the SINGLE EXCITATIONS *************************************************** + print*,'h1 = ',h1 + h1 = list_act_reverse(h1) + print*,'h1 = ',h1 + print*,'p1 = ',p1 + p1 = list_act_reverse(p1) + print*,'p1 = ',p1 + if(s1==1)then ! Mono alpha : - do k = 1, elec_beta_num - m = occ(k,2) + ! purely active part of the extra diagonal two body dm + do k = 1, n_occ_ab(2) + m = list_act_reverse(occ(k,2)) ! * c_I * c_J - call insert_into_two_body_dm_big_array( two_body_dm_ab_big_array,n_act_orb,n_act_orb,n_act_orb,n_act_orb,contrib,h1,p1,m,m) + call insert_into_two_body_dm_big_array( two_body_dm_ab_big_array_act,n_act_orb,n_act_orb,n_act_orb,n_act_orb,contrib,h1,p1,m,m) + enddo + + ! core <-> active part of the extra diagonal two body dm + do k = 1, n_occ_ab_core(2) + m = list_core_reverse(occ_core(k,2)) + ! * c_I * c_J + two_body_dm_ab_big_array_core_act(m,h1,p1) += 2.d0 * contrib + two_body_dm_ab_big_array_core_act(m,p1,h1) += 2.d0 * contrib enddo else ! Mono Beta : - do k = 1, elec_alpha_num - m = occ(k,1) + ! purely active part of the extra diagonal two body dm + do k = 1, n_occ_ab(1) + m = list_act_reverse(occ(k,1)) ! * c_I * c_J - call insert_into_two_body_dm_big_array(two_body_dm_ab_big_array,n_act_orb,n_act_orb,n_act_orb,n_act_orb,contrib,h1,p1,m,m) + call insert_into_two_body_dm_big_array(two_body_dm_ab_big_array_act,n_act_orb,n_act_orb,n_act_orb,n_act_orb,contrib,h1,p1,m,m) + enddo + + ! core <-> active part of the extra diagonal two body dm + do k = 1, n_occ_ab_core(1) + m = list_core_reverse(occ_core(k,1)) + ! * c_I * c_J + two_body_dm_ab_big_array_core_act(m,h1,p1) += 2.d0 * contrib + two_body_dm_ab_big_array_core_act(m,p1,h1) += 2.d0 * contrib enddo endif @@ -303,30 +465,39 @@ subroutine insert_into_two_body_dm_big_array(big_array,dim1,dim2,dim3,dim4,contr end -double precision function compute_extra_diag_two_body_dm_ab(r1,r2) +double precision function compute_extra_diag_two_body_dm_ab(r1,r2) implicit none + BEGIN_DOC +! compute the extra diagonal contribution to the alpha/bet two body density at r1, r2 + END_DOC + double precision :: r1(3), r2(3) + double precision :: compute_extra_diag_two_body_dm_ab_act,compute_extra_diag_two_body_dm_ab_core_act + compute_extra_diag_two_body_dm_ab = compute_extra_diag_two_body_dm_ab_act(r1,r2)+compute_extra_diag_two_body_dm_ab_core_act(r1,r2) +end + +double precision function compute_extra_diag_two_body_dm_ab_act(r1,r2) + implicit none + BEGIN_DOC +! compute the extra diagonal contribution to the two body density at r1, r2 +! involving ONLY THE ACTIVE PART, which means that the four index of the excitations +! involved in the two body density matrix are ACTIVE + END_DOC + PROVIDE n_act_orb double precision, intent(in) :: r1(3),r2(3) integer :: i,j,k,l - double precision :: mos_array_r1(mo_tot_num),mos_array_r2(mo_tot_num) + double precision :: mos_array_r1(n_act_orb),mos_array_r2(n_act_orb) double precision :: contrib - compute_extra_diag_two_body_dm_ab = 0.d0 -!call give_all_mos_at_r(r1,mos_array_r1) -!call give_all_mos_at_r(r2,mos_array_r2) + double precision :: contrib_tmp +!print*,'n_act_orb = ',n_act_orb + compute_extra_diag_two_body_dm_ab_act = 0.d0 call give_all_act_mos_at_r(r1,mos_array_r1) call give_all_act_mos_at_r(r2,mos_array_r2) do l = 1, n_act_orb ! p2 do k = 1, n_act_orb ! h2 do j = 1, n_act_orb ! p1 do i = 1,n_act_orb ! h1 - double precision :: contrib_tmp -! print*,'i,j',i,j -! print*,mos_array_r1(i) , mos_array_r1(j) -! print*,'k,l',k,l -! print*,mos_array_r2(k) * mos_array_r2(l) -! print*,'gama = ',two_body_dm_ab_big_array(i,j,k,l) -! pause contrib_tmp = mos_array_r1(i) * mos_array_r1(j) * mos_array_r2(k) * mos_array_r2(l) - compute_extra_diag_two_body_dm_ab += two_body_dm_ab_big_array(i,j,k,l) * contrib_tmp + compute_extra_diag_two_body_dm_ab_act += two_body_dm_ab_big_array_act(i,j,k,l) * contrib_tmp enddo enddo enddo @@ -334,13 +505,69 @@ double precision function compute_extra_diag_two_body_dm_ab(r1,r2) end -double precision function compute_diag_two_body_dm_ab(r1,r2) +double precision function compute_extra_diag_two_body_dm_ab_core_act(r1,r2) + implicit none + BEGIN_DOC +! compute the extra diagonal contribution to the two body density at r1, r2 +! involving ONLY THE ACTIVE PART, which means that the four index of the excitations +! involved in the two body density matrix are ACTIVE + END_DOC + double precision, intent(in) :: r1(3),r2(3) + integer :: i,j,k,l + double precision :: mos_array_act_r1(n_act_orb),mos_array_act_r2(n_act_orb) + double precision :: mos_array_core_r1(n_core_orb),mos_array_core_r2(n_core_orb) + double precision :: contrib_core_1,contrib_core_2 + double precision :: contrib_act_1,contrib_act_2 + double precision :: contrib_tmp + compute_extra_diag_two_body_dm_ab_core_act = 0.d0 + call give_all_act_mos_at_r(r1,mos_array_act_r1) + call give_all_act_mos_at_r(r2,mos_array_act_r2) + call give_all_core_mos_at_r(r1,mos_array_core_r1) + call give_all_core_mos_at_r(r2,mos_array_core_r2) + do i = 1, n_act_orb ! h1 + do j = 1, n_act_orb ! p1 + contrib_act_1 = mos_array_act_r1(i) * mos_array_act_r1(j) + contrib_act_2 = mos_array_act_r2(i) * mos_array_act_r2(j) + do k = 1,n_core_orb ! h2 + contrib_core_1 = mos_array_core_r1(k) * mos_array_core_r1(k) + contrib_core_2 = mos_array_core_r2(k) * mos_array_core_r2(k) + contrib_tmp = 0.5d0 * (contrib_act_1 * contrib_core_2 + contrib_act_2 * contrib_core_1) + compute_extra_diag_two_body_dm_ab_core_act += two_body_dm_ab_big_array_core_act(k,i,j) * contrib_tmp + enddo + enddo + enddo + +end + +double precision function compute_diag_two_body_dm_ab_core(r1,r2) implicit none double precision :: r1(3),r2(3) integer :: i,j,k,l - double precision :: mos_array_r1(mo_tot_num),mos_array_r2(mo_tot_num) + double precision :: mos_array_r1(n_core_orb_allocate),mos_array_r2(n_core_orb_allocate) double precision :: contrib,contrib_tmp - compute_diag_two_body_dm_ab = 0.d0 + compute_diag_two_body_dm_ab_core = 0.d0 + call give_all_core_mos_at_r(r1,mos_array_r1) + call give_all_core_mos_at_r(r2,mos_array_r2) + do l = 1, n_core_orb ! + contrib = mos_array_r2(l)*mos_array_r2(l) +! if(dabs(contrib).lt.threshld_two_bod_dm)cycle + do k = 1, n_core_orb ! + contrib_tmp = contrib * mos_array_r1(k)*mos_array_r1(k) +! if(dabs(contrib).lt.threshld_two_bod_dm)cycle + compute_diag_two_body_dm_ab_core += two_body_dm_ab_diag_core(k,l) * contrib_tmp + enddo + enddo + +end + + +double precision function compute_diag_two_body_dm_ab_act(r1,r2) + implicit none + double precision :: r1(3),r2(3) + integer :: i,j,k,l + double precision :: mos_array_r1(n_act_orb),mos_array_r2(n_act_orb) + double precision :: contrib,contrib_tmp + compute_diag_two_body_dm_ab_act = 0.d0 call give_all_act_mos_at_r(r1,mos_array_r1) call give_all_act_mos_at_r(r2,mos_array_r2) do l = 1, n_act_orb ! @@ -349,10 +576,44 @@ double precision function compute_diag_two_body_dm_ab(r1,r2) do k = 1, n_act_orb ! contrib_tmp = contrib * mos_array_r1(k)*mos_array_r1(k) ! if(dabs(contrib).lt.threshld_two_bod_dm)cycle - compute_diag_two_body_dm_ab += two_body_dm_ab_diag(k,l) * contrib_tmp + compute_diag_two_body_dm_ab_act += two_body_dm_ab_diag_act(k,l) * contrib_tmp enddo enddo - end +double precision function compute_diag_two_body_dm_ab_core_act(r1,r2) + implicit none + double precision :: r1(3),r2(3) + integer :: i,j,k,l + double precision :: mos_array_core_r1(n_core_orb_allocate),mos_array_core_r2(n_core_orb_allocate) + double precision :: mos_array_act_r1(n_act_orb),mos_array_act_r2(n_act_orb) + double precision :: contrib_core_1,contrib_core_2 + double precision :: contrib_act_1,contrib_act_2 + double precision :: contrib_tmp + compute_diag_two_body_dm_ab_core_act = 0.d0 + call give_all_act_mos_at_r(r1,mos_array_act_r1) + call give_all_act_mos_at_r(r2,mos_array_act_r2) + call give_all_core_mos_at_r(r1,mos_array_core_r1) + call give_all_core_mos_at_r(r2,mos_array_core_r2) +! if(dabs(contrib).lt.threshld_two_bod_dm)cycle + do k = 1, n_act_orb ! + contrib_act_1 = mos_array_act_r1(k) * mos_array_act_r1(k) + contrib_act_2 = mos_array_act_r2(k) * mos_array_act_r2(k) + contrib_tmp = 0.5d0 * (contrib_act_1 * contrib_act_2 + contrib_act_2 * contrib_act_1) +! if(dabs(contrib).lt.threshld_two_bod_dm)cycle + do l = 1, n_core_orb ! + contrib_core_1 = mos_array_core_r1(l) * mos_array_core_r1(l) + contrib_core_2 = mos_array_core_r2(l) * mos_array_core_r2(l) + compute_diag_two_body_dm_ab_core_act += two_body_dm_diag_core_act(l,k) * contrib_tmp + enddo + enddo +end +double precision function compute_diag_two_body_dm_ab(r1,r2) + implicit none + double precision,intent(in) :: r1(3),r2(3) + double precision :: compute_diag_two_body_dm_ab_act,compute_diag_two_body_dm_ab_core + double precision :: compute_diag_two_body_dm_ab_core_act + compute_diag_two_body_dm_ab = compute_diag_two_body_dm_ab_act(r1,r2)+compute_diag_two_body_dm_ab_core(r1,r2) & + + compute_diag_two_body_dm_ab_core_act(r1,r2) +end diff --git a/src/Integrals_Bielec/EZFIO.cfg b/src/Integrals_Bielec/EZFIO.cfg index 3834b121..feed02c1 100644 --- a/src/Integrals_Bielec/EZFIO.cfg +++ b/src/Integrals_Bielec/EZFIO.cfg @@ -5,6 +5,13 @@ interface: ezfio,provider,ocaml default: False ezfio_name: direct +[no_vvvv_integrals] +type: logical +doc: If True, do not compute the bielectronic integrals when 4 indices are virtual +interface: ezfio,provider,ocaml +default: False +ezfio_name: None + [disk_access_mo_integrals] type: Disk_access doc: Read/Write MO integrals from/to disk [ Write | Read | None ] diff --git a/src/Integrals_Bielec/mo_bi_integrals.irp.f b/src/Integrals_Bielec/mo_bi_integrals.irp.f index 69ca0733..3557772d 100644 --- a/src/Integrals_Bielec/mo_bi_integrals.irp.f +++ b/src/Integrals_Bielec/mo_bi_integrals.irp.f @@ -21,6 +21,7 @@ end BEGIN_PROVIDER [ logical, mo_bielec_integrals_in_map ] implicit none + integer(bit_kind) :: mask_ijkl(N_int,4) BEGIN_DOC ! If True, the map of MO bielectronic integrals is provided @@ -36,7 +37,44 @@ BEGIN_PROVIDER [ logical, mo_bielec_integrals_in_map ] endif endif - call add_integrals_to_map(full_ijkl_bitmask_4) + if(no_vvvv_integrals)then + integer :: i,j,k,l + ! (core+inact+act) ^ 4 + do i = 1,N_int + mask_ijkl(i,1) = core_inact_act_bitmask_4(i,1) + mask_ijkl(i,2) = core_inact_act_bitmask_4(i,2) + mask_ijkl(i,3) = core_inact_act_bitmask_4(i,3) + mask_ijkl(i,4) = core_inact_act_bitmask_4(i,4) + enddo + call add_integrals_to_map(mask_ijkl) + ! (core+inact+act) ^ 3 (virt) ^1 + do i = 1,N_int + mask_ijkl(i,1) = core_inact_act_bitmask_4(i,1) + mask_ijkl(i,2) = core_inact_act_bitmask_4(i,2) + mask_ijkl(i,3) = core_inact_act_bitmask_4(i,3) + mask_ijkl(i,4) = virt_bitmask(i,1) + enddo + call add_integrals_to_map(mask_ijkl) + ! (core+inact+act) ^ 2 (virt) ^2 + do i = 1,N_int + mask_ijkl(i,1) = core_inact_act_bitmask_4(i,1) + mask_ijkl(i,2) = core_inact_act_bitmask_4(i,2) + mask_ijkl(i,3) = virt_bitmask(i,1) + mask_ijkl(i,4) = virt_bitmask(i,1) + enddo + call add_integrals_to_map(mask_ijkl) + ! (core+inact+act) ^ 1 (virt) ^3 + do i = 1,N_int + mask_ijkl(i,1) = core_inact_act_bitmask_4(i,1) + mask_ijkl(i,2) = virt_bitmask(i,1) + mask_ijkl(i,3) = virt_bitmask(i,1) + mask_ijkl(i,4) = virt_bitmask(i,1) + enddo + call add_integrals_to_map(mask_ijkl) + + else + call add_integrals_to_map(full_ijkl_bitmask_4) + endif END_PROVIDER subroutine add_integrals_to_map(mask_ijkl) diff --git a/src/MO_Basis/mo_permutation.irp.f b/src/MO_Basis/mo_permutation.irp.f new file mode 100644 index 00000000..72f132d7 --- /dev/null +++ b/src/MO_Basis/mo_permutation.irp.f @@ -0,0 +1,20 @@ +program permut_mos + implicit none + integer :: mo1,mo2 + integer :: i,j,k,l + double precision :: mo_coef_tmp(ao_num_align,2) + print*,'Which MOs would you like to change ?' + read(5,*)mo1,mo2 + print*,'' + do i= 1,ao_num + mo_coef_tmp(i,1) = mo_coef(i,mo1) + mo_coef_tmp(i,2) = mo_coef(i,mo2) + enddo + do i = 1,ao_num + mo_coef(i,mo1) = mo_coef_tmp(i,2) + mo_coef(i,mo2) = mo_coef_tmp(i,1) + enddo + touch mo_coef + call save_mos + +end diff --git a/src/MO_Basis/print_aos.irp.f b/src/MO_Basis/print_aos.irp.f new file mode 100644 index 00000000..f6b3bedf --- /dev/null +++ b/src/MO_Basis/print_aos.irp.f @@ -0,0 +1,53 @@ +program pouet + implicit none + integer :: i,j,k + double precision :: r(3) + double precision, allocatable :: aos_array(:),mos_array(:),ao_ortho_array(:) + allocate(aos_array(ao_num),mos_array(mo_tot_num), ao_ortho_array(ao_num)) + integer :: nx,ny + double precision :: interval_x + double precision :: xmin,xmax + double precision :: dx + + double precision :: interval_y + double precision :: ymin,ymax + double precision :: dy + + double precision :: val_max + +!do i = 1, ao_num +! write(41,'(100(F16.10,X))'),ao_ortho_canonical_overlap(i,:) +!enddo + +!stop + + + xmin = nucl_coord(1,1)-6.d0 + xmax = nucl_coord(2,1)+6.d0 + interval_x = xmax - xmin +!interval_x = nucl_dist(1,3) + nx = 500 + dx = interval_x/dble(nx) +!dx = dabs(interval_x)/dble(nx) * 1.d0/sqrt(2.d0) + + r = 0.d0 + r(3) = xmin +!r(2) = nucl_coord(1,2) +!r(3) = nucl_coord(1,3) +!r(1) = nucl_coord(2,1) +!r(2) = 1.D0 +!r(3) = nucl_coord(2,3) + double precision :: dr(3) +!dr = 0.d0 +!dr(1) = -dx +!dr(3) = dx + do j = 1, nx+1 + call give_all_mos_at_r(r,mos_array) + write(37,'(100(F16.10,X))') r(3),mos_array(1)*mos_array(1) , mos_array(2)*mos_array(2), mos_array(1)*mos_array(2) + write(38,'(100(F16.10,X))') r(3),mos_array(1), mos_array(2), mos_array(1)*mos_array(2) +! write(38,'(100(F16.10,X))') r(3),mos_array(10), mos_array(2) - 0.029916d0 * mos_array(10),mos_array(2) + 0.029916d0 * mos_array(10) + r(3) += dx +! r += dr + enddo + deallocate(aos_array,mos_array, ao_ortho_array) +end diff --git a/src/MO_Basis/print_mo_in_space.irp.f b/src/MO_Basis/print_mo_in_space.irp.f new file mode 100644 index 00000000..5a2bc297 --- /dev/null +++ b/src/MO_Basis/print_mo_in_space.irp.f @@ -0,0 +1,50 @@ +program pouet + implicit none + integer :: i,j,k + double precision :: r(3) + double precision, allocatable :: aos_array(:),mos_array(:),ao_ortho_array(:) + allocate(aos_array(ao_num),mos_array(mo_tot_num), ao_ortho_array(ao_num)) + integer :: nx,ny + double precision :: interval_x + double precision :: xmin,xmax + double precision :: dx + + double precision :: interval_y + double precision :: ymin,ymax + double precision :: dy + + double precision :: val_max + +!do i = 1, ao_num +! write(41,'(100(F16.10,X))'),ao_ortho_canonical_overlap(i,:) +!enddo + +!stop + + + xmin = -4.d0 + xmax = 4.d0 + interval_x = xmax - xmin + nx = 100 + dx = dabs(interval_x)/dble(nx) + + r = 0.d0 +!r(3) = xmin + r(1) = xmin + val_max = 0.d0 + do j = 1, nx +! call give_all_aos_at_r(r,aos_array) + call give_all_mos_at_r(r,mos_array) + write(36,'(100(F16.10,X))') r(1), mos_array(1), mos_array(2), mos_array(1)* mos_array(2) + !write(36,'(100(F16.10,X))') r(1), mos_array(1), mos_array(2), mos_array(4) + !write(37,'(100(F16.10,X))') r(1),mos_array(1) * mos_array(2), mos_array(4)*mos_array(2) +! if(val_max.le.aos_array(1) * aos_array(2) )then +! val_max = aos_array(1) * aos_array(2) +! endif + r(1) += dx +! r(3) += dx + enddo +!write(40,'(100(F16.10,X))')nucl_coord(1,2),nucl_coord(1,3),val_max * 1.5d0 +!write(41,'(100(F16.10,X))')nucl_coord(2,2),nucl_coord(2,3),val_max * 1.5d0 + deallocate(aos_array,mos_array, ao_ortho_array) +end diff --git a/src/Nuclei/atomic_radii.irp.f b/src/Nuclei/atomic_radii.irp.f new file mode 100644 index 00000000..7b04a97b --- /dev/null +++ b/src/Nuclei/atomic_radii.irp.f @@ -0,0 +1,112 @@ +BEGIN_PROVIDER [ double precision, slater_bragg_radii, (100)] + implicit none + BEGIN_DOC + ! atomic radii in Angstrom defined in table I of JCP 41, 3199 (1964) Slater + ! execpt for the Hydrogen atom where we took the value of Becke (1988, JCP) + END_DOC + + slater_bragg_radii = 0.d0 + + slater_bragg_radii(1) = 0.35d0 + slater_bragg_radii(2) = 0.35d0 + + slater_bragg_radii(3) = 1.45d0 + slater_bragg_radii(4) = 1.05d0 + + slater_bragg_radii(5) = 0.85d0 + slater_bragg_radii(6) = 0.70d0 + slater_bragg_radii(7) = 0.65d0 + slater_bragg_radii(8) = 0.60d0 + slater_bragg_radii(9) = 0.50d0 + slater_bragg_radii(10) = 0.45d0 + + slater_bragg_radii(11) = 1.80d0 + slater_bragg_radii(12) = 1.70d0 + + slater_bragg_radii(13) = 1.50d0 + slater_bragg_radii(14) = 1.25d0 + slater_bragg_radii(15) = 1.10d0 + slater_bragg_radii(16) = 1.00d0 + slater_bragg_radii(17) = 1.00d0 + slater_bragg_radii(18) = 1.00d0 + + slater_bragg_radii(19) = 2.20d0 + slater_bragg_radii(20) = 1.80d0 + + + slater_bragg_radii(21) = 1.60d0 + slater_bragg_radii(22) = 1.40d0 + slater_bragg_radii(23) = 1.34d0 + slater_bragg_radii(24) = 1.40d0 + slater_bragg_radii(25) = 1.40d0 + slater_bragg_radii(26) = 1.40d0 + slater_bragg_radii(27) = 1.35d0 + slater_bragg_radii(28) = 1.35d0 + slater_bragg_radii(29) = 1.35d0 + slater_bragg_radii(30) = 1.35d0 + + slater_bragg_radii(31) = 1.30d0 + slater_bragg_radii(32) = 1.25d0 + slater_bragg_radii(33) = 1.15d0 + slater_bragg_radii(34) = 1.15d0 + slater_bragg_radii(35) = 1.15d0 + slater_bragg_radii(36) = 1.15d0 + +END_PROVIDER + +BEGIN_PROVIDER [double precision, slater_bragg_radii_ua, (100)] + implicit none + integer :: i + do i = 1, 100 + slater_bragg_radii_ua(i) = slater_bragg_radii(i) * 1.889725989d0 + enddo +END_PROVIDER + +BEGIN_PROVIDER [double precision, slater_bragg_radii_per_atom, (nucl_num)] + implicit none + integer :: i + do i = 1, nucl_num + slater_bragg_radii_per_atom(i) = slater_bragg_radii(int(nucl_charge(i))) + enddo +END_PROVIDER + +BEGIN_PROVIDER [double precision, slater_bragg_radii_per_atom_ua, (nucl_num)] + implicit none + integer :: i + do i = 1, nucl_num + slater_bragg_radii_per_atom_ua(i) = slater_bragg_radii_ua(int(nucl_charge(i))) + enddo +END_PROVIDER + +BEGIN_PROVIDER [double precision, slater_bragg_type_inter_distance, (nucl_num, nucl_num)] + implicit none + integer :: i,j + double precision :: xhi_tmp,u_ij + slater_bragg_type_inter_distance = 0.d0 + do i = 1, nucl_num + do j = i+1, nucl_num + xhi_tmp = slater_bragg_radii_per_atom(i) / slater_bragg_radii_per_atom(j) + u_ij = (xhi_tmp - 1.d0 ) / (xhi_tmp +1.d0) + slater_bragg_type_inter_distance(i,j) = u_ij / (u_ij * u_ij - 1.d0) + enddo + enddo +END_PROVIDER + +BEGIN_PROVIDER [double precision, slater_bragg_type_inter_distance_ua, (nucl_num, nucl_num)] + implicit none + integer :: i,j + double precision :: xhi_tmp,u_ij + slater_bragg_type_inter_distance_ua = 0.d0 + do i = 1, nucl_num + do j = i+1, nucl_num + xhi_tmp = slater_bragg_radii_per_atom_ua(i) / slater_bragg_radii_per_atom_ua(j) + u_ij = (xhi_tmp - 1.d0 ) / (xhi_tmp +1.d0) + slater_bragg_type_inter_distance_ua(i,j) = u_ij / (u_ij * u_ij - 1.d0) + if(slater_bragg_type_inter_distance_ua(i,j).gt.0.5d0)then + slater_bragg_type_inter_distance_ua(i,j) = 0.5d0 + else if( slater_bragg_type_inter_distance_ua(i,j) .le.-0.5d0)then + slater_bragg_type_inter_distance_ua(i,j) = -0.5d0 + endif + enddo + enddo +END_PROVIDER diff --git a/src/Utils/angular_integration.irp.f b/src/Utils/angular_integration.irp.f new file mode 100644 index 00000000..1efd4abc --- /dev/null +++ b/src/Utils/angular_integration.irp.f @@ -0,0 +1,2264 @@ +BEGIN_PROVIDER [integer, degree_max_integration_lebedev] + BEGIN_DOC +! integrate correctly a polynom of order "degree_max_integration_lebedev" + ! needed for the angular integration according to LEBEDEV formulae + END_DOC + implicit none + degree_max_integration_lebedev= 15 + +END_PROVIDER + +BEGIN_PROVIDER [integer, n_points_integration_angular_lebedev] + BEGIN_DOC +! Number of points needed for the angular integral + END_DOC + implicit none + if (degree_max_integration_lebedev == 3)then + n_points_integration_angular_lebedev = 6 + else if (degree_max_integration_lebedev == 5)then + n_points_integration_angular_lebedev = 14 + else if (degree_max_integration_lebedev == 7)then + n_points_integration_angular_lebedev = 26 + else if (degree_max_integration_lebedev == 9)then + n_points_integration_angular_lebedev = 38 + else if (degree_max_integration_lebedev == 11)then + n_points_integration_angular_lebedev = 50 + else if (degree_max_integration_lebedev == 13)then + n_points_integration_angular_lebedev = 74 + else if (degree_max_integration_lebedev == 15)then + n_points_integration_angular_lebedev = 86 + else if (degree_max_integration_lebedev == 17)then + n_points_integration_angular_lebedev = 110 + else if (degree_max_integration_lebedev == 19)then + n_points_integration_angular_lebedev = 146 + else if (degree_max_integration_lebedev == 21)then + n_points_integration_angular_lebedev = 170 + endif + +END_PROVIDER + + BEGIN_PROVIDER [double precision, theta_angular_integration_lebedev, (n_points_integration_angular_lebedev)] +&BEGIN_PROVIDER [double precision, phi_angular_integration_lebedev, (n_points_integration_angular_lebedev)] +&BEGIN_PROVIDER [double precision, weights_angular_integration_lebedev, (n_points_integration_angular_lebedev)] + implicit none + BEGIN_DOC +! Theta phi values together with the weights values for the angular integration : +! integral [dphi,dtheta] f(x,y,z) = 4 * pi * sum (1