diff --git a/ocaml/Message.ml b/ocaml/Message.ml index f47d6cec..ee5ff80c 100644 --- a/ocaml/Message.ml +++ b/ocaml/Message.ml @@ -280,9 +280,8 @@ type t = let of_string s = let l = String.split ~on:' ' s - |> List.map ~f:String.strip - |> List.map ~f:String.lowercase |> List.filter ~f:(fun x -> (String.strip x) <> "") + |> List.map ~f:String.lowercase in match l with | "add_task" :: state :: task -> diff --git a/ocaml/Qpackage.ml b/ocaml/Qpackage.ml index 3656327f..bd0d34fc 100644 --- a/ocaml/Qpackage.ml +++ b/ocaml/Qpackage.ml @@ -55,7 +55,6 @@ let executables = lazy ( In_channel.input_lines in_channel |> List.map ~f:(fun x -> let e = String.split ~on:' ' x - |> List.map ~f:String.strip |> List.filter ~f:(fun x -> x <> "") in match e with diff --git a/ocaml/Queuing_system.ml b/ocaml/Queuing_system.ml index 7a927a60..7407d78d 100644 --- a/ocaml/Queuing_system.ml +++ b/ocaml/Queuing_system.ml @@ -29,12 +29,14 @@ let add_task ~task q = q.next_task_id in { q with - queued = q.queued @ [ task_id ] ; + queued = task_id :: q.queued ; tasks = Map.add q.tasks ~key:task_id ~data:task ; next_task_id = Id.Task.increment task_id ; }, task_id + + let add_client q = let client_id = q.next_client_id diff --git a/ocaml/TaskServer.ml b/ocaml/TaskServer.ml index 1f882540..61eec19f 100644 --- a/ocaml/TaskServer.ml +++ b/ocaml/TaskServer.ml @@ -80,7 +80,7 @@ let stop ~port = Message.Terminate (Message.Terminate_msg.create ()) |> Message.to_string - |> ZMQ.Socket.send req_socket ; + |> ZMQ.Socket.send ~block:false req_socket ; let msg = ZMQ.Socket.recv req_socket @@ -158,13 +158,13 @@ let run ~port = let terminate () = running := false; Message.to_string ok - |> ZMQ.Socket.send rep_socket + |> ZMQ.Socket.send ~block:false rep_socket and newjob x = q := Queuing_system.create (); job := Some x; Message.to_string ok - |> ZMQ.Socket.send rep_socket + |> ZMQ.Socket.send ~block:false rep_socket and connect state msg = let push_address = @@ -180,7 +180,7 @@ let run ~port = Message.ConnectReply (Message.ConnectReply_msg.create ~state ~client_id ~push_address) |> Message.to_string - |> ZMQ.Socket.send rep_socket + |> ZMQ.Socket.send ~block:false rep_socket and disconnect state msg = let s, c = @@ -199,7 +199,7 @@ let run ~port = Message.DisconnectReply (Message.DisconnectReply_msg.create ~state ~finished) |> Message.to_string - |> ZMQ.Socket.send rep_socket + |> ZMQ.Socket.send ~block:false rep_socket and add_task state msg = let s, task = @@ -207,12 +207,53 @@ let run ~port = msg.Message.AddTask_msg.task in assert (s = state); - let new_q, task_id = - Queuing_system.add_task ~task !q - in - q := new_q; Message.to_string ok - |> ZMQ.Socket.send rep_socket + |> ZMQ.Socket.send ~block:false rep_socket + ; + begin + match + String.split ~on:' ' msg.Message.AddTask_msg.task + |> List.filter ~f:(fun x -> x <> "") + with + | "triangle" :: str_l :: [] -> + begin + let l = + Int.of_string str_l + in + for j=1 to l + do + let task = + Printf.sprintf "%d %s" j str_l + in + let new_q, _ = + Queuing_system.add_task ~task !q + in + q := new_q + done + end + | "range" :: str_i :: str_j :: [] -> + begin + let i, j = + Int.of_string str_i, + Int.of_string str_j + in + for k=i to (j+1) + do + let task = + Int.to_string k + in + let new_q, task_id = + Queuing_system.add_task ~task !q + in + q := new_q + done + end + | _ -> + let new_q, task_id = + Queuing_system.add_task ~task !q + in + q := new_q + end and get_task state msg = let s, client_id = @@ -231,7 +272,7 @@ let run ~port = | _ -> Message.Terminate (Message.Terminate_msg.create ()) in Message.to_string reply - |> ZMQ.Socket.send rep_socket + |> ZMQ.Socket.send ~block:false rep_socket and task_done state msg = let s, client_id, task_id = @@ -245,12 +286,12 @@ let run ~port = in q := new_q; Message.to_string ok - |> ZMQ.Socket.send rep_socket + |> ZMQ.Socket.send ~block:false rep_socket and error msg = Message.Error (Message.Error_msg.create msg) |> Message.to_string - |> ZMQ.Socket.send rep_socket + |> ZMQ.Socket.send ~block:false rep_socket in if (polling.(0) = Some ZMQ.Poll.In) then diff --git a/plugins/Hartree_Fock/huckel.irp.f b/plugins/Hartree_Fock/huckel.irp.f index 8f61f0cf..103de83a 100644 --- a/plugins/Hartree_Fock/huckel.irp.f +++ b/plugins/Hartree_Fock/huckel.irp.f @@ -8,8 +8,6 @@ subroutine huckel_guess double precision :: c character*(64) :: label - mo_coef = ao_ortho_lowdin_coef - TOUCH mo_coef label = "Guess" call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral, & size(mo_mono_elec_integral,1), & diff --git a/plugins/MP2/H_apply.irp.f b/plugins/MP2/H_apply.irp.f index a79e3879..a5489149 100644 --- a/plugins/MP2/H_apply.irp.f +++ b/plugins/MP2/H_apply.irp.f @@ -8,7 +8,7 @@ s.set_perturbation("Moller_plesset") print s s = H_apply("mp2_selection") -s.set_selection_pt2("Moller_plesset") +s.set_selection_pt2("Moller_Plesset") print s END_SHELL diff --git a/plugins/Perturbation/pt2_equations.irp.f b/plugins/Perturbation/pt2_equations.irp.f index f49ee2ff..8e40d0fd 100644 --- a/plugins/Perturbation/pt2_equations.irp.f +++ b/plugins/Perturbation/pt2_equations.irp.f @@ -123,8 +123,8 @@ subroutine pt2_moller_plesset ($arguments) call get_excitation(ref_bitmask,det_pert,exc,degree,phase,Nint) if (degree == 2) then call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) - delta_e = Fock_matrix_diag_mo(h1) + Fock_matrix_diag_mo(h2) - & - (Fock_matrix_diag_mo(p1) + Fock_matrix_diag_mo(p2)) + delta_e = (Fock_matrix_diag_mo(h1) + Fock_matrix_diag_mo(h2)) - & + (Fock_matrix_diag_mo(p1) + Fock_matrix_diag_mo(p2)) delta_e = 1.d0/delta_e else if (degree == 1) then call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) @@ -134,8 +134,14 @@ subroutine pt2_moller_plesset ($arguments) delta_e = 0.d0 endif - call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det,psi_selectors_size,n_st,i_H_psi_array) - h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint) + if (delta_e /= 0.d0) then + ! call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det,psi_selectors_size,n_st,i_H_psi_array) + call i_H_psi_minilist(det_pert,minilist,idx_minilist,N_minilist,psi_selectors_coef,Nint,N_minilist,psi_selectors_size,N_st,i_H_psi_array) + h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint) + else + i_H_psi_array(:) = 0.d0 + h = 0.d0 + endif do i =1,n_st H_pert_diag(i) = h c_pert(i) = i_H_psi_array(i) *delta_e diff --git a/src/Determinants/H_apply.template.f b/src/Determinants/H_apply.template.f index 58ae8b08..e949c0d2 100644 --- a/src/Determinants/H_apply.template.f +++ b/src/Determinants/H_apply.template.f @@ -261,6 +261,7 @@ subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl ! Build array of the non-zero integrals of second excitation $filter_integrals + if (ispin == 1) then integer :: jjj @@ -269,7 +270,7 @@ subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl i_b = occ_hole_tmp(kk,other_spin) ASSERT (i_b > 0) ASSERT (i_b <= mo_tot_num) - do jjj=1,N_elec_in_key_part_2(other_spin) ! particule + do jjj=1,N_elec_in_key_part_2(other_spin) ! particle j_b = occ_particle_tmp(jjj,other_spin) ASSERT (j_b > 0) ASSERT (j_b <= mo_tot_num) diff --git a/src/Integrals_Bielec/ao_bi_integrals.irp.f b/src/Integrals_Bielec/ao_bi_integrals.irp.f index 6987d06b..8bae1a1e 100644 --- a/src/Integrals_Bielec/ao_bi_integrals.irp.f +++ b/src/Integrals_Bielec/ao_bi_integrals.irp.f @@ -371,19 +371,16 @@ BEGIN_PROVIDER [ logical, ao_bielec_integrals_in_map ] integer(ZMQ_PTR) :: zmq_to_qp_run_socket call new_parallel_job(zmq_to_qp_run_socket,'ao_integrals') + character*(32) :: task do l=1,ao_num - do j = 1, l - if (ao_overlap_abs(j,l) < ao_integrals_threshold) then - cycle - endif - write(task,*) j, l - call add_task_to_taskserver(zmq_to_qp_run_socket,task) - enddo + write(task,*) 'triangle', l + call add_task_to_taskserver(zmq_to_qp_run_socket,task) enddo + external :: ao_bielec_integrals_in_map_slave_inproc, ao_bielec_integrals_in_map_collector call new_parallel_threads(ao_bielec_integrals_in_map_slave_inproc, ao_bielec_integrals_in_map_collector) - + call end_parallel_job(zmq_to_qp_run_socket,'ao_integrals') print*, 'Sorting the map' diff --git a/src/MOGuess/H_CORE_guess.irp.f b/src/MOGuess/H_CORE_guess.irp.f index b65fe07d..d3e2eef9 100644 --- a/src/MOGuess/H_CORE_guess.irp.f +++ b/src/MOGuess/H_CORE_guess.irp.f @@ -5,8 +5,6 @@ program H_CORE_guess END_DOC implicit none character*(64) :: label - mo_coef = ao_ortho_lowdin_coef - TOUCH mo_coef label = "Guess" call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral, & size(mo_mono_elec_integral,1), & diff --git a/src/MOGuess/guess_overlap.irp.f b/src/MOGuess/guess_overlap.irp.f index c2f090e5..7d75e118 100644 --- a/src/MOGuess/guess_overlap.irp.f +++ b/src/MOGuess/guess_overlap.irp.f @@ -5,8 +5,6 @@ program guess_mimi implicit none character*(64) :: label - mo_coef = ao_ortho_lowdin_coef - TOUCH mo_coef label = "Guess" call mo_as_eigvectors_of_mo_matrix(ao_overlap, & size(ao_overlap,1), & diff --git a/src/MOGuess/h_core_guess_routine.irp.f b/src/MOGuess/h_core_guess_routine.irp.f index 605c7c8a..23899160 100644 --- a/src/MOGuess/h_core_guess_routine.irp.f +++ b/src/MOGuess/h_core_guess_routine.irp.f @@ -4,8 +4,6 @@ subroutine hcore_guess END_DOC implicit none character*(64) :: label - mo_coef = ao_ortho_lowdin_coef - TOUCH mo_coef label = "Guess" call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral, & size(mo_mono_elec_integral,1), & diff --git a/src/MOGuess/mo_ortho_lowdin.irp.f b/src/MOGuess/mo_ortho_lowdin.irp.f index c73bb553..90672b5e 100644 --- a/src/MOGuess/mo_ortho_lowdin.irp.f +++ b/src/MOGuess/mo_ortho_lowdin.irp.f @@ -1,4 +1,3 @@ - BEGIN_PROVIDER [double precision, ao_ortho_lowdin_coef, (ao_num_align,ao_num)] implicit none BEGIN_DOC @@ -8,13 +7,9 @@ BEGIN_PROVIDER [double precision, ao_ortho_lowdin_coef, (ao_num_align,ao_num)] END_DOC integer :: i,j,k,l double precision :: tmp_matrix(ao_num_align,ao_num),accu + tmp_matrix(:,:) = 0.d0 do j=1, ao_num - do i=1, ao_num - tmp_matrix(i,j) = 0.d0 - enddo - enddo - do i=1, ao_num - tmp_matrix(i,i) = 1.d0 + tmp_matrix(j,j) = 1.d0 enddo call ortho_lowdin(ao_overlap,ao_num_align,ao_num,tmp_matrix,ao_num_align,ao_num) do i=1, ao_num @@ -23,6 +18,7 @@ BEGIN_PROVIDER [double precision, ao_ortho_lowdin_coef, (ao_num_align,ao_num)] enddo enddo END_PROVIDER + BEGIN_PROVIDER [double precision, ao_ortho_lowdin_overlap, (ao_num_align,ao_num)] implicit none BEGIN_DOC @@ -40,7 +36,7 @@ BEGIN_PROVIDER [double precision, ao_ortho_lowdin_overlap, (ao_num_align,ao_num) do j=1, ao_num c = 0.d0 do l=1, ao_num - c = ao_ortho_lowdin_coef(j,l) * ao_overlap(k,l) + c += ao_ortho_lowdin_coef(j,l) * ao_overlap(k,l) enddo do i=1, ao_num ao_ortho_lowdin_overlap(i,j) += ao_ortho_lowdin_coef(i,k) * c diff --git a/src/MOGuess/mo_ortho_canonical.irp.f b/src/MO_Basis/ao_ortho_canonical.irp.f similarity index 59% rename from src/MOGuess/mo_ortho_canonical.irp.f rename to src/MO_Basis/ao_ortho_canonical.irp.f index a4027555..f9862abb 100644 --- a/src/MOGuess/mo_ortho_canonical.irp.f +++ b/src/MO_Basis/ao_ortho_canonical.irp.f @@ -1,21 +1,20 @@ - -BEGIN_PROVIDER [double precision, ao_ortho_canonical_coef, (ao_num_align,ao_num)] + BEGIN_PROVIDER [ double precision, ao_ortho_canonical_coef, (ao_num_align,ao_num)] +&BEGIN_PROVIDER [ integer, ao_ortho_canonical_num ] implicit none BEGIN_DOC ! matrix of the coefficients of the mos generated by the ! orthonormalization by the S^{-1/2} canonical transformation of the aos ! ao_ortho_canonical_coef(i,j) = coefficient of the ith ao on the jth ao_ortho_canonical orbital END_DOC - integer :: i,j,k,l - tmp_matrix(:,:) = 0.d0 - do j=1, ao_num - tmp_matrix(j,j) = 1.d0 + integer :: i + ao_ortho_canonical_coef(:,:) = 0.d0 + do i=1,ao_num + ao_ortho_canonical_coef(i,i) = 1.d0 enddo - call ortho_canonical(ao_overlap,ao_num_align,ao_num,ao_ortho_canonical_coef,ao_num_align,mo_tot_num) - SOFT_TOUCH mo_tot_num + call ortho_canonical(ao_overlap,ao_num_align,ao_num,ao_ortho_canonical_coef,ao_num_align,ao_ortho_canonical_num) END_PROVIDER -BEGIN_PROVIDER [double precision, ao_ortho_canonical_overlap, (ao_num_align,ao_num)] +BEGIN_PROVIDER [double precision, ao_ortho_canonical_overlap, (ao_ortho_canonical_num,ao_ortho_canonical_num)] implicit none BEGIN_DOC ! overlap matrix of the ao_ortho_canonical @@ -23,19 +22,19 @@ BEGIN_PROVIDER [double precision, ao_ortho_canonical_overlap, (ao_num_align,ao_n END_DOC integer :: i,j,k,l double precision :: c - do j=1, ao_num - do i=1, ao_num + do j=1, ao_ortho_canonical_num + do i=1, ao_ortho_canonical_num ao_ortho_canonical_overlap(i,j) = 0.d0 enddo enddo - do k=1, ao_num - do j=1, ao_num + do j=1, ao_ortho_canonical_num + do k=1, ao_num c = 0.d0 do l=1, ao_num - c += ao_ortho_canonical_coef(j,l) * ao_overlap(k,l) + c += ao_ortho_canonical_coef(l,j) * ao_overlap(l,k) enddo - do i=1, ao_num - ao_ortho_canonical_overlap(i,j) += ao_ortho_canonical_coef(i,k) * c + do i=1, ao_ortho_canonical_num + ao_ortho_canonical_overlap(i,j) += ao_ortho_canonical_coef(k,i) * c enddo enddo enddo diff --git a/src/MO_Basis/mos.irp.f b/src/MO_Basis/mos.irp.f index 016b48ad..5756b7de 100644 --- a/src/MO_Basis/mos.irp.f +++ b/src/MO_Basis/mos.irp.f @@ -9,8 +9,9 @@ BEGIN_PROVIDER [ integer, mo_tot_num ] if (exists) then call ezfio_get_mo_basis_mo_tot_num(mo_tot_num) else - mo_tot_num = ao_num + mo_tot_num = ao_ortho_canonical_num endif + call write_int(6,mo_tot_num,'mo_tot_num') ASSERT (mo_tot_num > 0) END_PROVIDER @@ -56,7 +57,14 @@ END_PROVIDER deallocate(buffer) else ! Orthonormalized AO basis - mo_coef = 0. + do i=1,mo_tot_num + do j=1,ao_num + mo_coef(j,i) = ao_ortho_canonical_coef(j,i) + enddo + do j=ao_num+1,ao_num_align + mo_coef(j,i) = 0.d0 + enddo + enddo endif END_PROVIDER diff --git a/src/Utils/LinearAlgebra.irp.f b/src/Utils/LinearAlgebra.irp.f index 14c1d8dd..a859a913 100644 --- a/src/Utils/LinearAlgebra.irp.f +++ b/src/Utils/LinearAlgebra.irp.f @@ -42,11 +42,81 @@ subroutine svd(A,LDA,U,LDU,D,Vt,LDVt,m,n) end +subroutine ortho_canonical(overlap,LDA,N,C,LDC,m) + implicit none + BEGIN_DOC + ! Compute C_new=C_old.U.s^-1/2 canonical orthogonalization. + ! + ! overlap : overlap matrix + ! + ! LDA : leftmost dimension of overlap array + ! + ! N : Overlap matrix is NxN (array is (LDA,N) ) + ! + ! C : Coefficients of the vectors to orthogonalize. On exit, + ! orthogonal vectors + ! + ! LDC : leftmost dimension of C + ! + ! m : Coefficients matrix is MxN, ( array is (LDC,N) ) + ! + END_DOC + + integer, intent(in) :: lda, ldc, n + integer, intent(out) :: m + double precision, intent(in) :: overlap(lda,n) + double precision, intent(inout) :: C(ldc,n) + double precision, allocatable :: U(:,:) + double precision, allocatable :: Vt(:,:) + double precision, allocatable :: D(:) + double precision, allocatable :: S_half(:,:) + !DEC$ ATTRIBUTES ALIGN : 64 :: U, Vt, D + integer :: info, i, j + + allocate (U(ldc,n), Vt(lda,n), D(n), S_half(lda,n)) + + call svd(overlap,lda,U,ldc,D,Vt,lda,n,n) + + m=n + do i=1,n + if ( D(i) >= 1.d-4 ) then + D(i) = 1.d0/dsqrt(D(i)) + else + m = i-1 + exit + endif + enddo + do i=m+1,n + D(i) = 0.d0 + enddo + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP SHARED(S_half,U,D,Vt,n,C,m) & + !$OMP PRIVATE(i,j) + + !$OMP DO + do j=1,n + do i=1,n + S_half(i,j) = U(i,j)*D(j) + enddo + do i=1,n + U(i,j) = C(i,j) + enddo + enddo + !$OMP END DO + + !$OMP END PARALLEL + + call dgemm('N','N',n,m,n,1.d0,U,size(U,1),S_half,size(S_half,1),0.d0,C,size(C,1)) + deallocate (U, Vt, D, S_half) + +end + subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m) implicit none BEGIN_DOC - ! Compute C_new=C_old.S^-1/2 canonical orthogonalization. + ! Compute C_new=C_old.S^-1/2 orthogonalization. ! ! overlap : overlap matrix ! @@ -81,7 +151,7 @@ subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m) !$OMP DO do i=1,n - if ( D(i) < 1.d-4 ) then + if ( D(i) < 1.d-5 ) then D(i) = 0.d0 else D(i) = 1.d0/dsqrt(D(i))