diff --git a/plugins/Full_CI/H_apply.irp.f b/plugins/Full_CI/H_apply.irp.f index 4f63e29e..a37e2165 100644 --- a/plugins/Full_CI/H_apply.irp.f +++ b/plugins/Full_CI/H_apply.irp.f @@ -5,24 +5,20 @@ from generate_h_apply import * s = H_apply("FCI") s.set_selection_pt2("epstein_nesbet_2x2") #s.set_selection_pt2("qdpt") -s.unset_skip() print s s = H_apply("FCI_PT2") s.set_perturbation("epstein_nesbet_2x2") #s.set_perturbation("qdpt") -s.unset_skip() s.unset_openmp() print s s = H_apply("FCI_no_selection") s.set_selection_pt2("dummy") -s.unset_skip() print s s = H_apply("FCI_mono") s.set_selection_pt2("epstein_nesbet_2x2") -s.unset_skip() s.unset_double_excitations() s.unset_openmp() print s diff --git a/plugins/Properties/print_hcc.irp.f b/plugins/Properties/print_hcc.irp.f index 45bca5e6..7f74f996 100644 --- a/plugins/Properties/print_hcc.irp.f +++ b/plugins/Properties/print_hcc.irp.f @@ -2,5 +2,16 @@ program print_hcc_main implicit none read_wf = .True. touch read_wf - call print_hcc +! call print_hcc + call routine end + + +subroutine routine + implicit none + integer :: i + do i = 1, mo_tot_num + write(*,'(1000(F16.10,X))')one_body_dm_mo_beta(i,:,1) + enddo +end + diff --git a/plugins/mrcepa0/energy.irp.f b/plugins/mrcepa0/energy.irp.f index 757a81e3..7eaf5562 100644 --- a/plugins/mrcepa0/energy.irp.f +++ b/plugins/mrcepa0/energy.irp.f @@ -13,6 +13,8 @@ BEGIN_PROVIDER [ double precision, mrcc_E0_denominator, (N_states) ] END_DOC if (initialize_mrcc_E0_denominator) then mrcc_E0_denominator(1:N_states) = psi_energy(1:N_states) +! mrcc_E0_denominator(1:N_states) = HF_energy - nuclear_repulsion +! mrcc_E0_denominator(1:N_states) = barycentric_electronic_energy(1:N_states) call write_double(6,mrcc_E0_denominator(1)+nuclear_repulsion, 'mrcc Energy denominator') else mrcc_E0_denominator = -huge(1.d0) diff --git a/plugins/mrcepa0/mrcc_stoch_routines.irp.f b/plugins/mrcepa0/mrcc_stoch_routines.irp.f index ac1573f3..1ad9b8da 100644 --- a/plugins/mrcepa0/mrcc_stoch_routines.irp.f +++ b/plugins/mrcepa0/mrcc_stoch_routines.irp.f @@ -11,6 +11,7 @@ subroutine ZMQ_mrcc(E, mrcc, delta, delta_s2, relative_error) implicit none character(len=64000) :: task + integer(ZMQ_PTR) :: zmq_to_qp_run_socket, zmq_socket_pull integer, external :: omp_get_thread_num double precision, intent(in) :: relative_error, E @@ -40,6 +41,7 @@ subroutine ZMQ_mrcc(E, mrcc, delta, delta_s2, relative_error) print *, ' Samples Energy Stat. Error Seconds ' print *, '========== ================= ================= =================' + call new_parallel_job(zmq_to_qp_run_socket,zmq_socket_pull, 'mrcc') integer, external :: zmq_put_psi @@ -63,7 +65,7 @@ subroutine ZMQ_mrcc(E, mrcc, delta, delta_s2, relative_error) ! do i=1,comb_teeth ! print *, "TOOTH", first_det_of_teeth(i+1) - first_det_of_teeth(i) ! end do - + integer(ZMQ_PTR), external :: new_zmq_to_qp_run_socket integer :: ipos ipos=1 @@ -75,6 +77,7 @@ subroutine ZMQ_mrcc(E, mrcc, delta, delta_s2, relative_error) if (add_task_to_taskserver(zmq_to_qp_run_socket,trim(task(1:ipos))) == -1) then stop 'Unable to add task to task server' endif + ipos=1 endif else @@ -105,6 +108,7 @@ subroutine ZMQ_mrcc(E, mrcc, delta, delta_s2, relative_error) i = omp_get_thread_num() if (i==0) then call mrcc_collector(zmq_socket_pull,E(mrcc_stoch_istate), relative_error, delta, delta_s2, mrcc) + else call mrcc_slave_inproc(i) endif @@ -123,14 +127,18 @@ subroutine mrcc_slave_inproc(i) end + subroutine mrcc_collector(zmq_socket_pull, E, relative_error, delta, delta_s2, mrcc) + use dress_types use f77_zmq use bitmasks implicit none + integer(ZMQ_PTR), intent(in) :: zmq_socket_pull + double precision, intent(in) :: relative_error, E double precision, intent(out) :: mrcc(N_states) double precision, allocatable :: cp(:,:,:,:) @@ -236,11 +244,13 @@ subroutine mrcc_collector(zmq_socket_pull, E, relative_error, delta, delta_s2, m end if end do + integer, external :: zmq_delete_tasks if (zmq_delete_tasks(zmq_to_qp_run_socket,zmq_socket_pull,task_id,ntask,more) == -1) then stop 'Unable to delete tasks' endif + time = omp_get_wtime() @@ -262,6 +272,7 @@ subroutine mrcc_collector(zmq_socket_pull, E, relative_error, delta, delta_s2, m !!!!!!!!!!!! double precision :: su, su2, eqt, avg, E0, val integer, external :: zmq_abort + su = 0d0 su2 = 0d0 @@ -284,6 +295,7 @@ subroutine mrcc_collector(zmq_socket_pull, E, relative_error, delta, delta_s2, m if ((dabs(eqt) < relative_error .and. cps_N(cur_cp) >= 30) .or. total_computed == N_det_generators) then ! Termination !print '(G10.3, 2X, F16.7, 2X, G16.3, 2X, F16.4, A20)', Nabove(tooth), avg+E, eqt, time-time0, '' + ! print *, "GREPME", cur_cp, E+E0+avg, eqt, time-time0, total_computed if (zmq_abort(zmq_to_qp_run_socket) == -1) then call sleep(1) @@ -296,6 +308,7 @@ subroutine mrcc_collector(zmq_socket_pull, E, relative_error, delta, delta_s2, m if (cur_cp > old_cur_cp) then old_cur_cp = cur_cp ! print *, "GREPME", cur_cp, E+E0+avg, eqt, time-time0, total_computed + !print '(G10.3, 2X, F16.7, 2X, G16.3, 2X, F16.4, A20)', Nabove(tooth), avg+E, eqt, time-time0, '' endif endif @@ -326,6 +339,8 @@ subroutine mrcc_collector(zmq_socket_pull, E, relative_error, delta, delta_s2, m mrcc(1) = E call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket) + + call end_zmq_pull_socket(zmq_socket_pull) end subroutine diff --git a/plugins/mrcepa0/mrcepa0_general.irp.f b/plugins/mrcepa0/mrcepa0_general.irp.f index 6a4a76cf..f58008a0 100644 --- a/plugins/mrcepa0/mrcepa0_general.irp.f +++ b/plugins/mrcepa0/mrcepa0_general.irp.f @@ -47,6 +47,10 @@ subroutine run(N_st,energy) enddo call diagonalize_ci_dressed(lambda) E_new = mrcc_e0_denominator(1) !sum(ci_energy_dressed(1:N_states)) + +! if (.true.) then +! provide delta_ij_mrcc_pouet +! endif delta_E = (E_new - E_old)/dble(N_states) print *, '' call write_double(6,thresh_mrcc,"thresh_mrcc") diff --git a/src/Determinants/density_matrix.irp.f b/src/Determinants/density_matrix.irp.f index b948182d..86c7fcd4 100644 --- a/src/Determinants/density_matrix.irp.f +++ b/src/Determinants/density_matrix.irp.f @@ -99,7 +99,7 @@ END_PROVIDER ! Alpha and beta one-body density matrix for each state END_DOC - integer :: j,k,l,m + integer :: j,k,l,m,k_a,k_b integer :: occ(N_int*bit_kind_size,2) double precision :: ck, cl, ckl double precision :: phase @@ -114,7 +114,7 @@ END_PROVIDER one_body_dm_mo_alpha = 0.d0 one_body_dm_mo_beta = 0.d0 !$OMP PARALLEL DEFAULT(NONE) & - !$OMP PRIVATE(j,k,l,m,occ,ck, cl, ckl,phase,h1,h2,p1,p2,s1,s2, degree,exc, & + !$OMP PRIVATE(j,k,k_a,k_b,l,m,occ,ck, cl, ckl,phase,h1,h2,p1,p2,s1,s2, degree,exc, & !$OMP tmp_a, tmp_b, n_occ, krow, kcol, lrow, lcol, tmp_det, tmp_det2)& !$OMP SHARED(psi_det,psi_coef,N_int,N_states,elec_alpha_num,& !$OMP elec_beta_num,one_body_dm_mo_alpha,one_body_dm_mo_beta,N_det,& @@ -124,28 +124,31 @@ END_PROVIDER !$OMP psi_bilinear_matrix_values, psi_bilinear_matrix_transp_values) allocate(tmp_a(mo_tot_num,mo_tot_num,N_states), tmp_b(mo_tot_num,mo_tot_num,N_states) ) tmp_a = 0.d0 - tmp_b = 0.d0 - !$OMP DO SCHEDULE(guided) - do k=1,N_det - krow = psi_bilinear_matrix_rows(k) - kcol = psi_bilinear_matrix_columns(k) - tmp_det(:,1) = psi_det_alpha_unique(:,krow) - tmp_det(:,2) = psi_det_beta_unique (:,kcol) + !$OMP DO SCHEDULE(dynamic,64) + do k_a=1,N_det + krow = psi_bilinear_matrix_rows(k_a) + ASSERT (krow <= N_det_alpha_unique) + + kcol = psi_bilinear_matrix_columns(k_a) + ASSERT (kcol <= N_det_beta_unique) + + tmp_det(1:N_int,1) = psi_det_alpha_unique(1:N_int,krow) + tmp_det(1:N_int,2) = psi_det_beta_unique (1:N_int,kcol) + + ! Diagonal part + ! ------------- + call bitstring_to_list_ab(tmp_det, occ, n_occ, N_int) do m=1,N_states - ck = psi_bilinear_matrix_values(k,m)*psi_bilinear_matrix_values(k,m) + ck = psi_bilinear_matrix_values(k_a,m)*psi_bilinear_matrix_values(k_a,m) do l=1,elec_alpha_num j = occ(l,1) tmp_a(j,j,m) += ck enddo - do l=1,elec_beta_num - j = occ(l,2) - tmp_b(j,j,m) += ck - enddo enddo - if (k == N_det) cycle - l = k+1 + if (k_a == N_det) cycle + l = k_a+1 lrow = psi_bilinear_matrix_rows(l) lcol = psi_bilinear_matrix_columns(l) ! Fix beta determinant, loop over alphas @@ -157,7 +160,7 @@ END_PROVIDER call get_mono_excitation_spin(tmp_det(1,1),tmp_det2,exc,phase,N_int) call decode_exc_spin(exc,h1,p1,h2,p2) do m=1,N_states - ckl = psi_bilinear_matrix_values(k,m)*psi_bilinear_matrix_values(l,m) * phase + ckl = psi_bilinear_matrix_values(k_a,m)*psi_bilinear_matrix_values(l,m) * phase tmp_a(h1,p1,m) += ckl tmp_a(p1,h1,m) += ckl enddo @@ -168,20 +171,52 @@ END_PROVIDER lcol = psi_bilinear_matrix_columns(l) enddo - l = psi_bilinear_matrix_order_reverse(k) - if (l == N_det) cycle - l = l+1 - ! Fix alpha determinant, loop over betas + enddo + !$OMP END DO NOWAIT + + !$OMP CRITICAL + one_body_dm_mo_alpha(:,:,:) = one_body_dm_mo_alpha(:,:,:) + tmp_a(:,:,:) + !$OMP END CRITICAL + deallocate(tmp_a) + + tmp_b = 0.d0 + !$OMP DO SCHEDULE(dynamic,64) + do k_b=1,N_det + krow = psi_bilinear_matrix_transp_rows(k_b) + ASSERT (krow <= N_det_alpha_unique) + + kcol = psi_bilinear_matrix_transp_columns(k_b) + ASSERT (kcol <= N_det_beta_unique) + + tmp_det(1:N_int,1) = psi_det_alpha_unique(1:N_int,krow) + tmp_det(1:N_int,2) = psi_det_beta_unique (1:N_int,kcol) + + ! Diagonal part + ! ------------- + + call bitstring_to_list_ab(tmp_det, occ, n_occ, N_int) + do m=1,N_states + ck = psi_bilinear_matrix_transp_values(k_b,m)*psi_bilinear_matrix_transp_values(k_b,m) + do l=1,elec_beta_num + j = occ(l,2) + tmp_b(j,j,m) += ck + enddo + enddo + + if (k_b == N_det) cycle + l = k_b+1 lrow = psi_bilinear_matrix_transp_rows(l) lcol = psi_bilinear_matrix_transp_columns(l) + ! Fix beta determinant, loop over alphas do while ( lrow == krow ) - tmp_det2(:) = psi_det_beta_unique (:, lcol) + tmp_det2(:) = psi_det_beta_unique(:, lcol) call get_excitation_degree_spin(tmp_det(1,2),tmp_det2,degree,N_int) if (degree == 1) then + exc = 0 call get_mono_excitation_spin(tmp_det(1,2),tmp_det2,exc,phase,N_int) call decode_exc_spin(exc,h1,p1,h2,p2) do m=1,N_states - ckl = psi_bilinear_matrix_values(k,m)*psi_bilinear_matrix_transp_values(l,m) * phase + ckl = psi_bilinear_matrix_transp_values(k_b,m)*psi_bilinear_matrix_transp_values(l,m) * phase tmp_b(h1,p1,m) += ckl tmp_b(p1,h1,m) += ckl enddo @@ -195,12 +230,10 @@ END_PROVIDER enddo !$OMP END DO NOWAIT !$OMP CRITICAL - one_body_dm_mo_alpha(:,:,:) = one_body_dm_mo_alpha(:,:,:) + tmp_a(:,:,:) - !$OMP END CRITICAL - !$OMP CRITICAL one_body_dm_mo_beta(:,:,:) = one_body_dm_mo_beta(:,:,:) + tmp_b(:,:,:) !$OMP END CRITICAL - deallocate(tmp_a,tmp_b) + + deallocate(tmp_b) !$OMP END PARALLEL END_PROVIDER diff --git a/src/Integrals_Bielec/mo_bi_integrals.irp.f b/src/Integrals_Bielec/mo_bi_integrals.irp.f index 479b373d..65594810 100644 --- a/src/Integrals_Bielec/mo_bi_integrals.irp.f +++ b/src/Integrals_Bielec/mo_bi_integrals.irp.f @@ -128,7 +128,7 @@ BEGIN_PROVIDER [ logical, mo_bielec_integrals_in_map ] mo_coef, size(mo_coef,1), & 1, 1, 1, 1, ao_num, ao_num, ao_num, ao_num, & 1, 1, 1, 1, mo_num, mo_num, mo_num, mo_num) -! + ! call four_index_transform_block(ao_integrals_map,mo_integrals_map, & ! mo_coef, size(mo_coef,1), & ! 1, 1, 1, 1, ao_num, ao_num, ao_num, ao_num, &