diff --git a/plugins/MRPT_Utils/H_apply.irp.f b/plugins/MRPT_Utils/H_apply.irp.f new file mode 100644 index 00000000..18bfc86f --- /dev/null +++ b/plugins/MRPT_Utils/H_apply.irp.f @@ -0,0 +1,26 @@ +use bitmasks +BEGIN_SHELL [ /usr/bin/env python ] +from generate_h_apply import * + +s = H_apply("mrpt") +s.data["parameters"] = ", delta_ij_, Ndet" +s.data["declarations"] += """ + integer, intent(in) :: Ndet + double precision, intent(in) :: delta_ij_(Ndet,Ndet,*) +""" +s.data["keys_work"] = "call mrpt_dress(delta_ij_,Ndet,i_generator,key_idx,keys_out,N_int,iproc,key_mask)" +s.data["params_post"] += ", delta_ij_, Ndet" +s.data["params_main"] += "delta_ij_, Ndet" +s.data["decls_main"] += """ + integer, intent(in) :: Ndet + double precision, intent(in) :: delta_ij_(Ndet,Ndet,*) +""" +s.data["finalization"] = "" +s.data["copy_buffer"] = "" +s.data["generate_psi_guess"] = "" +s.data["size_max"] = "3072" +print s + + +END_SHELL + diff --git a/plugins/MRPT_Utils/MRPT_Utils.main.irp.f b/plugins/MRPT_Utils/MRPT_Utils.main.irp.f index c8140f70..1e2fe8a8 100644 --- a/plugins/MRPT_Utils/MRPT_Utils.main.irp.f +++ b/plugins/MRPT_Utils/MRPT_Utils.main.irp.f @@ -11,8 +11,13 @@ end subroutine routine_3 implicit none !provide fock_virt_total_spin_trace - provide energy_cas_dyall - print*, 'nuclear_reuplsion = ',nuclear_repulsion + provide delta_ij + + print *, 'N_det = ', N_det + print *, 'N_states = ', N_states + print *, 'PT2 = ', second_order_pt_new(1) + print *, 'E = ', CI_energy + print *, 'E+PT2 = ', CI_energy+second_order_pt_new(1) end diff --git a/plugins/MRPT_Utils/excitations_cas.irp.f b/plugins/MRPT_Utils/excitations_cas.irp.f new file mode 100644 index 00000000..acd70e3e --- /dev/null +++ b/plugins/MRPT_Utils/excitations_cas.irp.f @@ -0,0 +1,390 @@ +subroutine apply_exc_to_psi(orb,hole_particle,spin_exc, & + norm_out,psi_in_out,psi_in_out_coef, ndet,dim_psi_in,dim_psi_coef,N_states_in) + use bitmasks + implicit none + integer, intent(in) :: orb, hole_particle,spin_exc,N_states_in,ndet,dim_psi_in,dim_psi_coef + double precision, intent(out) :: norm_out(N_states_in) + integer(bit_kind), intent(inout) :: psi_in_out(N_int,2,dim_psi_in) + double precision, intent(inout) :: psi_in_out_coef(dim_psi_coef,N_states_in) + BEGIN_DOC + ! apply a contracted excitation to psi_in_out whose coefficients + ! are psi_in_out_coef + ! hole_particle = 1 ===> creation of an electron in psi_in_out + ! = -1 ===> annhilation of an electron in psi_in_out + ! orb ===> is the index of orbital where you want wether to create or + ! annhilate an electron + ! spin_exc ===> is the spin of the electron (1 == alpha) (2 == beta) + ! the wave function gets out normalized to unity + ! + ! norm_out is the sum of the squared of the coefficients + ! on which the excitation has been possible + END_DOC + + integer :: elec_num_tab_local(2) + elec_num_tab_local = 0 + do i = 1, ndet + if( psi_in_out_coef (i,1) .ne. 0.d0)then + do j = 1, N_int + elec_num_tab_local(1) += popcnt(psi_in_out(j,1,i)) + elec_num_tab_local(2) += popcnt(psi_in_out(j,2,i)) + enddo + exit + endif + enddo + integer :: i,j,accu_elec + if(hole_particle == 1)then + do i = 1, ndet + call set_bit_to_integer(orb,psi_in_out(1,spin_exc,i),N_int) + accu_elec = 0 + do j = 1, N_int + accu_elec += popcnt(psi_in_out(j,spin_exc,i)) + enddo + if(accu_elec .ne. elec_num_tab_local(spin_exc)+1)then + do j = 1, N_int + psi_in_out(j,1,i) = 0_bit_kind + psi_in_out(j,2,i) = 0_bit_kind + enddo + do j = 1, N_states_in + psi_in_out_coef(i,j) = 0.d0 + enddo + endif + enddo + else if (hole_particle == -1)then + do i = 1, ndet + call clear_bit_to_integer(orb,psi_in_out(1,spin_exc,i),N_int) + accu_elec = 0 + do j = 1, N_int + accu_elec += popcnt(psi_in_out(j,spin_exc,i)) + enddo + if(accu_elec .ne. elec_num_tab_local(spin_exc)-1)then + do j = 1, N_int + psi_in_out(j,1,i) = 0_bit_kind + psi_in_out(j,2,i) = 0_bit_kind + enddo + do j = 1, N_states_in + psi_in_out_coef(i,j) = 0.d0 + enddo + endif + enddo + endif + norm_out = 0.d0 + double precision :: norm_factor + do j = 1, N_states_in + do i = 1, ndet + norm_out(j) += psi_in_out_coef(i,j) * psi_in_out_coef(i,j) + enddo + if(norm_out(j).le.1.d-10)then + norm_factor = 0.d0 + else + norm_factor = 1.d0/(dsqrt(norm_out(j))) + endif + do i = 1, ndet + psi_in_out_coef(i,j) = psi_in_out_coef(i,j) * norm_factor + enddo + enddo +end + + +double precision function diag_H_mat_elem_no_elec_check(det_in,Nint) + implicit none + BEGIN_DOC + ! Computes + END_DOC + integer,intent(in) :: Nint + integer(bit_kind),intent(in) :: det_in(Nint,2) + + integer :: i, j, iorb, jorb + integer :: occ(Nint*bit_kind_size,2) + integer :: elec_num_tab_local(2) + + diag_H_mat_elem_no_elec_check = 0.d0 + call bitstring_to_list(det_in(1,1), occ(1,1), elec_num_tab_local(1), N_int) + call bitstring_to_list(det_in(1,2), occ(1,2), elec_num_tab_local(2), N_int) + ! alpha - alpha + do i = 1, elec_num_tab_local(1) + iorb = occ(i,1) + diag_H_mat_elem_no_elec_check += mo_mono_elec_integral(iorb,iorb) + do j = i+1, elec_num_tab_local(1) + jorb = occ(j,1) + diag_H_mat_elem_no_elec_check += mo_bielec_integral_jj_anti(jorb,iorb) + enddo + enddo + + ! beta - beta + do i = 1, elec_num_tab_local(2) + iorb = occ(i,2) + diag_H_mat_elem_no_elec_check += mo_mono_elec_integral(iorb,iorb) + do j = i+1, elec_num_tab_local(2) + jorb = occ(j,2) + diag_H_mat_elem_no_elec_check += mo_bielec_integral_jj_anti(jorb,iorb) + enddo + enddo + + + ! alpha - beta + do i = 1, elec_num_tab_local(2) + iorb = occ(i,2) + do j = 1, elec_num_tab_local(1) + jorb = occ(j,1) + diag_H_mat_elem_no_elec_check += mo_bielec_integral_jj(jorb,iorb) + enddo + enddo + +end + +subroutine a_operator_no_check(iorb,ispin,key,hjj,Nint,na,nb) + use bitmasks + implicit none + BEGIN_DOC + ! Needed for diag_H_mat_elem + END_DOC + integer, intent(in) :: iorb, ispin, Nint + integer, intent(inout) :: na, nb + integer(bit_kind), intent(inout) :: key(Nint,2) + double precision, intent(inout) :: hjj + + integer :: occ(Nint*bit_kind_size,2) + integer :: other_spin + integer :: k,l,i + integer :: tmp(2) + + ASSERT (iorb > 0) + ASSERT (ispin > 0) + ASSERT (ispin < 3) + ASSERT (Nint > 0) + + k = ishft(iorb-1,-bit_kind_shift)+1 + ASSERT (k > 0) + l = iorb - ishft(k-1,bit_kind_shift)-1 + key(k,ispin) = ibclr(key(k,ispin),l) + other_spin = iand(ispin,1)+1 + + !DIR$ FORCEINLINE + call bitstring_to_list_ab(key, occ, tmp, Nint) + na = na-1 + + hjj = hjj - mo_mono_elec_integral(iorb,iorb) + + ! Same spin + do i=1,na + hjj = hjj - mo_bielec_integral_jj_anti(occ(i,ispin),iorb) + enddo + + ! Opposite spin + do i=1,nb + hjj = hjj - mo_bielec_integral_jj(occ(i,other_spin),iorb) + enddo + +end + + +subroutine ac_operator_no_check(iorb,ispin,key,hjj,Nint,na,nb) + use bitmasks + implicit none + BEGIN_DOC + ! Needed for diag_H_mat_elem + END_DOC + integer, intent(in) :: iorb, ispin, Nint + integer, intent(inout) :: na, nb + integer(bit_kind), intent(inout) :: key(Nint,2) + double precision, intent(inout) :: hjj + + integer :: occ(Nint*bit_kind_size,2) + integer :: other_spin + integer :: k,l,i + + ASSERT (iorb > 0) + ASSERT (ispin > 0) + ASSERT (ispin < 3) + ASSERT (Nint > 0) + + integer :: tmp(2) + !DIR$ FORCEINLINE + call bitstring_to_list_ab(key, occ, tmp, Nint) + + k = ishft(iorb-1,-bit_kind_shift)+1 + ASSERT (k > 0) + l = iorb - ishft(k-1,bit_kind_shift)-1 + key(k,ispin) = ibset(key(k,ispin),l) + other_spin = iand(ispin,1)+1 + + hjj = hjj + mo_mono_elec_integral(iorb,iorb) + + print*,'na.nb = ',na,nb + ! Same spin + do i=1,na + hjj = hjj + mo_bielec_integral_jj_anti(occ(i,ispin),iorb) + enddo + + ! Opposite spin + do i=1,nb + hjj = hjj + mo_bielec_integral_jj(occ(i,other_spin),iorb) + enddo + na = na+1 +end + + +subroutine i_H_j_dyall(key_i,key_j,Nint,hij) + use bitmasks + implicit none + BEGIN_DOC + ! Returns where i and j are determinants + END_DOC + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2) + double precision, intent(out) :: hij + + integer :: exc(0:2,2,2) + integer :: degree + double precision :: get_mo_bielec_integral_schwartz + integer :: m,n,p,q + integer :: i,j,k + integer :: occ(Nint*bit_kind_size,2) + double precision :: diag_H_mat_elem_no_elec_check, phase,phase_2 + integer :: n_occ_ab(2) + logical :: has_mipi(Nint*bit_kind_size) + double precision :: mipi(Nint*bit_kind_size), miip(Nint*bit_kind_size) + PROVIDE mo_bielec_integrals_in_map mo_integrals_map + + ASSERT (Nint > 0) + ASSERT (Nint == N_int) + + hij = 0.d0 + !DIR$ FORCEINLINE + call get_excitation_degree(key_i,key_j,degree,Nint) + select case (degree) + case (2) + call get_double_excitation(key_i,key_j,exc,phase,Nint) + if (exc(0,1,1) == 1) then + ! Mono alpha, mono beta + hij = phase*get_mo_bielec_integral_schwartz( & + exc(1,1,1), & + exc(1,1,2), & + exc(1,2,1), & + exc(1,2,2) ,mo_integrals_map) + else if (exc(0,1,1) == 2) then + ! Double alpha + hij = phase*(get_mo_bielec_integral_schwartz( & + exc(1,1,1), & + exc(2,1,1), & + exc(1,2,1), & + exc(2,2,1) ,mo_integrals_map) - & + get_mo_bielec_integral_schwartz( & + exc(1,1,1), & + exc(2,1,1), & + exc(2,2,1), & + exc(1,2,1) ,mo_integrals_map) ) + else if (exc(0,1,2) == 2) then + ! Double beta + hij = phase*(get_mo_bielec_integral_schwartz( & + exc(1,1,2), & + exc(2,1,2), & + exc(1,2,2), & + exc(2,2,2) ,mo_integrals_map) - & + get_mo_bielec_integral_schwartz( & + exc(1,1,2), & + exc(2,1,2), & + exc(2,2,2), & + exc(1,2,2) ,mo_integrals_map) ) + endif + case (1) + call get_mono_excitation(key_i,key_j,exc,phase,Nint) + !DIR$ FORCEINLINE + call bitstring_to_list_ab(key_i, occ, n_occ_ab, Nint) + has_mipi = .False. + if (exc(0,1,1) == 1) then + ! Mono alpha + m = exc(1,1,1) + p = exc(1,2,1) + do k = 1, n_occ_ab(1) + i = occ(k,1) + if (.not.has_mipi(i)) then + mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) + miip(i) = get_mo_bielec_integral_schwartz(m,i,i,p,mo_integrals_map) + has_mipi(i) = .True. + endif + enddo + do k = 1, n_occ_ab(2) + i = occ(k,2) + if (.not.has_mipi(i)) then + mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) + has_mipi(i) = .True. + endif + enddo + + do k = 1, n_occ_ab(1) + hij = hij + mipi(occ(k,1)) - miip(occ(k,1)) + enddo + do k = 1, n_occ_ab(2) + hij = hij + mipi(occ(k,2)) + enddo + + else + ! Mono beta + m = exc(1,1,2) + p = exc(1,2,2) + do k = 1, n_occ_ab(2) + i = occ(k,2) + if (.not.has_mipi(i)) then + mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) + miip(i) = get_mo_bielec_integral_schwartz(m,i,i,p,mo_integrals_map) + has_mipi(i) = .True. + endif + enddo + do k = 1, n_occ_ab(1) + i = occ(k,1) + if (.not.has_mipi(i)) then + mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) + has_mipi(i) = .True. + endif + enddo + + do k = 1, n_occ_ab(1) + hij = hij + mipi(occ(k,1)) + enddo + do k = 1, n_occ_ab(2) + hij = hij + mipi(occ(k,2)) - miip(occ(k,2)) + enddo + + endif + hij = phase*(hij + mo_mono_elec_integral(m,p) + fock_operator_active_from_core_inact(m,p) ) + + case (0) + hij = diag_H_mat_elem_no_elec_check(key_i,Nint) + end select +end + + +subroutine u0_H_dyall_u0(energies,psi_in,psi_in_coef,ndet,dim_psi_in,dim_psi_coef,N_states_in,state_target) + use bitmasks + implicit none + integer, intent(in) :: N_states_in,ndet,dim_psi_in,dim_psi_coef,state_target + integer(bit_kind), intent(in) :: psi_in(N_int,2,dim_psi_in) + double precision, intent(in) :: psi_in_coef(dim_psi_coef,N_states_in) + double precision, intent(out) :: energies(N_states_in) + + integer :: i,j + double precision :: hij,accu + energies = 0.d0 + accu = 0.d0 + double precision, allocatable :: psi_coef_tmp(:) + allocate(psi_coef_tmp(ndet)) + + do i = 1, ndet + psi_coef_tmp(i) = psi_in_coef(i,state_target) + enddo + + double precision :: hij_bis + do i = 1, ndet + if(psi_coef_tmp(i)==0.d0)cycle + do j = 1, ndet + if(psi_coef_tmp(j)==0.d0)cycle + call i_H_j_dyall(psi_in(1,1,i),psi_in(1,1,j),N_int,hij) +! call i_H_j(psi_in(1,1,i),psi_in(1,1,j),N_int,hij_bis) +! print*, hij_bis,hij + accu += psi_coef_tmp(i) * psi_coef_tmp(j) * hij + enddo + enddo + energies(state_target) = accu + deallocate(psi_coef_tmp) +end diff --git a/plugins/MRPT_Utils/mrpt_dress.irp.f b/plugins/MRPT_Utils/mrpt_dress.irp.f new file mode 100644 index 00000000..f0b58d24 --- /dev/null +++ b/plugins/MRPT_Utils/mrpt_dress.irp.f @@ -0,0 +1,161 @@ +use omp_lib +use bitmasks + +BEGIN_PROVIDER [ integer(omp_lock_kind), psi_ref_bis_lock, (psi_det_size) ] + implicit none + BEGIN_DOC + ! Locks on ref determinants to fill delta_ij + END_DOC + integer :: i + do i=1,psi_det_size + call omp_init_lock( psi_ref_bis_lock(i) ) + enddo + +END_PROVIDER + + +subroutine mrpt_dress(delta_ij_, Ndet,i_generator,n_selected,det_buffer,Nint,iproc,key_mask) + use bitmasks + implicit none + + integer, intent(in) :: i_generator,n_selected, Nint, iproc + integer, intent(in) :: Ndet + integer(bit_kind),intent(in) :: key_mask(Nint, 2) + integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected) + double precision, intent(inout) :: delta_ij_(Ndet,Ndet,*) + + + integer :: i,j,k,l + integer :: idx_alpha(0:psi_det_size) + integer :: degree_alpha(psi_det_size) + logical :: fullMatch + + double precision :: delta_e_array(psi_det_size) + double precision :: hij_array(psi_det_size) + + integer(bit_kind) :: tq(Nint,2,n_selected) + integer :: N_tq + + double precision :: hialpha + integer :: i_state, i_alpha + + integer(bit_kind),allocatable :: miniList(:,:,:) + integer,allocatable :: idx_miniList(:) + integer :: N_miniList, leng + double precision :: delta_e_final,hij_tmp + integer :: index_i,index_j + + + leng = max(N_det_generators, N_det) + allocate(miniList(Nint, 2, leng), idx_miniList(leng)) + + !create_minilist_find_previous(key_mask, fullList, miniList, N_fullList, N_miniList, fullMatch, Nint) + call create_minilist_find_previous(key_mask, psi_det_generators, miniList, i_generator-1, N_miniList, fullMatch, Nint) + + if(fullMatch) then + return + end if + + + call find_connections_previous(i_generator,n_selected,det_buffer,Nint,tq,N_tq,miniList,N_minilist) + + if(N_tq > 0) then + call create_minilist(key_mask, psi_det, miniList, idx_miniList, N_det, N_minilist, Nint) + end if + + + do i_alpha=1,N_tq + call get_excitation_degree_vector(miniList,tq(1,1,i_alpha),degree_alpha,Nint,N_minilist,idx_alpha) + + do j=1,idx_alpha(0) + idx_alpha(j) = idx_miniList(idx_alpha(j)) + enddo + + do i = 1,idx_alpha(0) + index_i = idx_alpha(i) + call get_delta_e_dyall(psi_det(1,1,index_i),tq(1,1,i_alpha),delta_e_final) + call i_h_j(tq(1,1,i_alpha),psi_det(1,1,index_i),Nint,hialpha) + delta_e_array(index_i) = 1.d0/delta_e_final + hij_array(index_i) = hialpha + enddo + + do i=1,idx_alpha(0) + index_i = idx_alpha(i) + hij_tmp = hij_array(index_i) + call omp_set_lock( psi_ref_bis_lock(index_i) ) + do j = 1, idx_alpha(0) + index_j = idx_alpha(j) + do i_state=1,N_states + delta_ij_(index_i,index_j,i_state) += hij_array(index_j) * hij_tmp * delta_e_array(index_j) + enddo + enddo + call omp_unset_lock( psi_ref_bis_lock(index_i)) + enddo + enddo + deallocate(miniList, idx_miniList) +end + + + + BEGIN_PROVIDER [ integer(bit_kind), gen_det_sorted, (N_int,2,N_det_generators,2) ] +&BEGIN_PROVIDER [ integer, gen_det_shortcut, (0:N_det_generators,2) ] +&BEGIN_PROVIDER [ integer, gen_det_version, (N_int, N_det_generators,2) ] +&BEGIN_PROVIDER [ integer, gen_det_idx, (N_det_generators,2) ] + gen_det_sorted(:,:,:,1) = psi_det_generators(:,:,:N_det_generators) + gen_det_sorted(:,:,:,2) = psi_det_generators(:,:,:N_det_generators) + call sort_dets_ab_v(gen_det_sorted(:,:,:,1), gen_det_idx(:,1), gen_det_shortcut(0:,1), gen_det_version(:,:,1), N_det_generators, N_int) + call sort_dets_ba_v(gen_det_sorted(:,:,:,2), gen_det_idx(:,2), gen_det_shortcut(0:,2), gen_det_version(:,:,2), N_det_generators, N_int) +END_PROVIDER + + +subroutine find_connections_previous(i_generator,n_selected,det_buffer,Nint,tq,N_tq,miniList,N_miniList) + + use bitmasks + implicit none + + integer, intent(in) :: i_generator,n_selected, Nint + + integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected) + integer :: i,j,k,m + logical :: is_in_wavefunction + integer :: degree(psi_det_size) + integer :: idx(0:psi_det_size) + logical :: good + + integer(bit_kind), intent(out) :: tq(Nint,2,n_selected) + integer, intent(out) :: N_tq + + + integer :: nt,ni + logical, external :: is_connected_to + + + integer(bit_kind),intent(in) :: miniList(Nint,2,N_det_generators) + integer,intent(in) :: N_miniList + + + + N_tq = 0 + + + i_loop : do i=1,N_selected + if(is_connected_to(det_buffer(1,1,i), miniList, Nint, N_miniList)) then + cycle + end if + + if (.not. is_in_wavefunction(det_buffer(1,1,i),Nint,N_det)) then + N_tq += 1 + do k=1,N_int + tq(k,1,N_tq) = det_buffer(k,1,i) + tq(k,2,N_tq) = det_buffer(k,2,i) + enddo + endif + enddo i_loop +end + + + + + + + diff --git a/plugins/MRPT_Utils/mrpt_utils.irp.f b/plugins/MRPT_Utils/mrpt_utils.irp.f new file mode 100644 index 00000000..38c33377 --- /dev/null +++ b/plugins/MRPT_Utils/mrpt_utils.irp.f @@ -0,0 +1,22 @@ + + BEGIN_PROVIDER [ double precision, delta_ij, (N_det,N_det,N_states) ] +&BEGIN_PROVIDER [ double precision, second_order_pt_new, (N_states) ] + implicit none + BEGIN_DOC + ! Dressing matrix in N_det basis + END_DOC + integer :: i,j,m + delta_ij = 0.d0 + call H_apply_mrpt(delta_ij,N_det) + double precision :: accu + accu = 0.d0 + do i = 1, N_det + do j = 1, N_det + accu += delta_ij(i,j,1) * psi_coef(i,1) * psi_coef(j,1) + enddo + write(*,'(1000(F16.10,x))')delta_ij(i,:,:) + enddo + print*, 'accu = ',accu + second_order_pt_new(1) = accu +END_PROVIDER + diff --git a/plugins/Properties/print_spin_density.irp.f b/plugins/Properties/print_spin_density.irp.f index 9daa6fb7..b9cbe4e8 100644 --- a/plugins/Properties/print_spin_density.irp.f +++ b/plugins/Properties/print_spin_density.irp.f @@ -14,7 +14,7 @@ subroutine routine double precision, allocatable :: aos_array(:) allocate(aos_array(ao_num)) r = 0.d0 - r(3) = z_min + r(1) = z_min do i = 1, N_z_pts call give_all_aos_at_r(r,aos_array) accu = 0.d0 @@ -28,8 +28,8 @@ subroutine routine 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 + r(1) += delta_z + write(33,'(100(f16.10,X))')r(1),accu,accu_alpha,accu_beta enddo diff --git a/src/Determinants/NEEDED_CHILDREN_MODULES b/src/Determinants/NEEDED_CHILDREN_MODULES index 5505ce78..566762ba 100644 --- a/src/Determinants/NEEDED_CHILDREN_MODULES +++ b/src/Determinants/NEEDED_CHILDREN_MODULES @@ -1 +1 @@ -Integrals_Monoelec Integrals_Bielec \ No newline at end of file +Integrals_Monoelec Integrals_Bielec Hartree_Fock diff --git a/src/MO_Basis/print_mo_in_space.irp.f b/src/MO_Basis/print_mo_in_space.irp.f index 5a2bc297..a5a324ed 100644 --- a/src/MO_Basis/print_mo_in_space.irp.f +++ b/src/MO_Basis/print_mo_in_space.irp.f @@ -35,7 +35,7 @@ program pouet 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(3), mos_array(17), mos_array(23) !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