diff --git a/README.md b/README.md index 3caef9a0..c9e1b12d 100644 --- a/README.md +++ b/README.md @@ -137,6 +137,10 @@ interface: ezfio #FAQ +### Opam error: cryptokit + +You need to install `gmp-dev`. + ### Error: ezfio_* is already defined. #### Why ? diff --git a/configure b/configure index 128d7126..cf8d1e03 100755 --- a/configure +++ b/configure @@ -102,7 +102,7 @@ curl = Info( default_path=join(QP_ROOT_BIN, "curl")) zlib = Info( - url='http://zlib.net/zlib-1.2.8.tar.gz', + url='http://www.zlib.net/zlib-1.2.10.tar.gz', description=' zlib', default_path=join(QP_ROOT_LIB, "libz.a")) @@ -496,15 +496,21 @@ def create_ninja_and_rc(l_installed): 'export LIBRARY_PATH="${QP_ROOT}"/lib:"${LIBRARY_PATH}"', 'export C_INCLUDE_PATH="${C_INCLUDE_PATH}":"${QP_ROOT}"/include', '', - 'source ${QP_ROOT}/install/EZFIO/Bash/ezfio.sh', "", - 'source ${HOME}/.opam/opam-init/init.sh > /dev/null 2> /dev/null || true', + 'source ${QP_ROOT}/install/EZFIO/Bash/ezfio.sh', '', '# Choose the correct network interface', '# export QP_NIC=ib0', '# export QP_NIC=eth0', - "" + '' ] + qp_opam_root = os.getenv('OPAMROOT') + if not qp_opam_root: + qp_opam_root = '${HOME}/.opam' + l_rc.append('export QP_OPAM={0}'.format(qp_opam_root)) + l_rc.append('source ${QP_OPAM}/opam-init/init.sh > /dev/null 2> /dev/null || true') + l_rc.append('') + path = join(QP_ROOT, "quantum_package.rc") with open(path, "w+") as f: f.write("\n".join(l_rc)) diff --git a/ocaml/Basis.ml b/ocaml/Basis.ml index 869fb132..797d53f2 100644 --- a/ocaml/Basis.ml +++ b/ocaml/Basis.ml @@ -36,9 +36,11 @@ let read_element in_channel at_number element = -let to_string_general ~fmt ~atom_sep b = +let to_string_general ~fmt ~atom_sep ?ele_array b = let new_nucleus n = - Printf.sprintf "Atom %d" n + match ele_array with + | None -> Printf.sprintf "Atom %d" n + | Some x -> Printf.sprintf "%s" (Element.to_string x.(n-1)) in let rec do_work accu current_nucleus = function | [] -> List.rev accu @@ -56,12 +58,12 @@ let to_string_general ~fmt ~atom_sep b = do_work [new_nucleus 1] 1 b |> String.concat ~sep:"\n" -let to_string_gamess = - to_string_general ~fmt:Gto.Gamess ~atom_sep:"" +let to_string_gamess ?ele_array = + to_string_general ?ele_array ~fmt:Gto.Gamess ~atom_sep:"" -let to_string_gaussian b = +let to_string_gaussian ?ele_array b = String.concat ~sep:"\n" - [ to_string_general ~fmt:Gto.Gaussian ~atom_sep:"****" b ; "****" ] + [ to_string_general ?ele_array ~fmt:Gto.Gaussian ~atom_sep:"****" b ; "****" ] let to_string ?(fmt=Gto.Gamess) = match fmt with diff --git a/ocaml/Basis.mli b/ocaml/Basis.mli index 249c14f9..41ddc184 100644 --- a/ocaml/Basis.mli +++ b/ocaml/Basis.mli @@ -14,7 +14,7 @@ val read_element : in_channel -> Nucl_number.t -> Element.t -> (Gto.t * Nucl_number.t) list (** Convert the basis to a string *) -val to_string : ?fmt:Gto.fmt -> (Gto.t * Nucl_number.t) list -> string +val to_string : ?fmt:Gto.fmt -> ?ele_array:Element.t array -> (Gto.t * Nucl_number.t) list -> string (** Convert the basis to an MD5 hash *) val to_md5 : (Gto.t * Nucl_number.t) list -> MD5.t diff --git a/plugins/All_singles/README.rst b/plugins/All_singles/README.rst index d3888edc..8836ddd6 100644 --- a/plugins/All_singles/README.rst +++ b/plugins/All_singles/README.rst @@ -15,6 +15,7 @@ Needed Modules * `Properties `_ * `Selectors_no_sorted `_ * `Utils `_ +* `Davidson `_ Documentation ============= diff --git a/plugins/All_singles/tree_dependency.png b/plugins/All_singles/tree_dependency.png new file mode 100644 index 00000000..e69de29b diff --git a/plugins/CAS_SD/.gitignore b/plugins/CAS_SD/.gitignore index 380d6cbf..57b1926f 100644 --- a/plugins/CAS_SD/.gitignore +++ b/plugins/CAS_SD/.gitignore @@ -3,6 +3,7 @@ .ninja_log AO_Basis Bitmask +Davidson Determinants Electrons Ezfio_files diff --git a/plugins/CAS_SD/README.rst b/plugins/CAS_SD/README.rst index 11f5d4cc..20ffa64f 100644 --- a/plugins/CAS_SD/README.rst +++ b/plugins/CAS_SD/README.rst @@ -107,6 +107,7 @@ Needed Modules * `Perturbation `_ * `Selectors_full `_ * `Generators_CAS `_ +* `Davidson `_ Documentation ============= @@ -193,31 +194,6 @@ h_apply_cas_s_selected_monoexc Assume N_int is already provided. -h_apply_cas_s_selected_no_skip - Calls H_apply on the HF determinant and selects all connected single and double - excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script. - - -h_apply_cas_s_selected_no_skip_diexc - Undocumented - - -h_apply_cas_s_selected_no_skip_diexcorg - Generate all double excitations of key_in using the bit masks of holes and - particles. - Assume N_int is already provided. - - -h_apply_cas_s_selected_no_skip_diexcp - Undocumented - - -h_apply_cas_s_selected_no_skip_monoexc - Generate all single excitations of key_in using the bit masks of holes and - particles. - Assume N_int is already provided. - - h_apply_cas_sd Calls H_apply on the HF determinant and selects all connected single and double excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script. diff --git a/plugins/Full_CI/.gitignore b/plugins/Full_CI/.gitignore index 674f56da..70d637ea 100644 --- a/plugins/Full_CI/.gitignore +++ b/plugins/Full_CI/.gitignore @@ -3,6 +3,7 @@ .ninja_log AO_Basis Bitmask +Davidson Determinants Electrons Ezfio_files @@ -28,7 +29,6 @@ full_ci full_ci_no_skip irpf90.make irpf90_entities -micro_pt2 tags target_pt2 var_pt2_ratio \ No newline at end of file diff --git a/plugins/Full_CI/README.rst b/plugins/Full_CI/README.rst index 750db44c..77a0bd64 100644 --- a/plugins/Full_CI/README.rst +++ b/plugins/Full_CI/README.rst @@ -16,6 +16,7 @@ Needed Modules * `Perturbation `_ * `Selectors_full `_ * `Generators_full `_ +* `Davidson `_ Documentation ============= @@ -77,6 +78,31 @@ h_apply_fci_monoexc Assume N_int is already provided. +h_apply_fci_no_selection + Calls H_apply on the HF determinant and selects all connected single and double + excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script. + + +h_apply_fci_no_selection_diexc + Undocumented + + +h_apply_fci_no_selection_diexcorg + Generate all double excitations of key_in using the bit masks of holes and + particles. + Assume N_int is already provided. + + +h_apply_fci_no_selection_diexcp + Undocumented + + +h_apply_fci_no_selection_monoexc + Generate all single excitations of key_in using the bit masks of holes and + particles. + Assume N_int is already provided. + + h_apply_fci_no_skip Calls H_apply on the HF determinant and selects all connected single and double excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script. @@ -144,118 +170,6 @@ h_apply_fci_pt2_slave_tcp Computes a buffer over the network -h_apply_pt2_mono_delta_rho - Calls H_apply on the HF determinant and selects all connected single and double - excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script. - - -h_apply_pt2_mono_delta_rho_diexc - Undocumented - - -h_apply_pt2_mono_delta_rho_diexcorg - Generate all double excitations of key_in using the bit masks of holes and - particles. - Assume N_int is already provided. - - -h_apply_pt2_mono_delta_rho_diexcp - Undocumented - - -h_apply_pt2_mono_delta_rho_monoexc - Generate all single excitations of key_in using the bit masks of holes and - particles. - Assume N_int is already provided. - - -h_apply_pt2_mono_di_delta_rho - Calls H_apply on the HF determinant and selects all connected single and double - excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script. - - -h_apply_pt2_mono_di_delta_rho_diexc - Undocumented - - -h_apply_pt2_mono_di_delta_rho_diexcorg - Generate all double excitations of key_in using the bit masks of holes and - particles. - Assume N_int is already provided. - - -h_apply_pt2_mono_di_delta_rho_diexcp - Undocumented - - -h_apply_pt2_mono_di_delta_rho_monoexc - Generate all single excitations of key_in using the bit masks of holes and - particles. - Assume N_int is already provided. - - -h_apply_select_mono_delta_rho - Calls H_apply on the HF determinant and selects all connected single and double - excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script. - - -h_apply_select_mono_delta_rho_diexc - Undocumented - - -h_apply_select_mono_delta_rho_diexcorg - Generate all double excitations of key_in using the bit masks of holes and - particles. - Assume N_int is already provided. - - -h_apply_select_mono_delta_rho_diexcp - Undocumented - - -h_apply_select_mono_delta_rho_monoexc - Generate all single excitations of key_in using the bit masks of holes and - particles. - Assume N_int is already provided. - - -h_apply_select_mono_di_delta_rho - Calls H_apply on the HF determinant and selects all connected single and double - excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script. - - -h_apply_select_mono_di_delta_rho_diexc - Undocumented - - -h_apply_select_mono_di_delta_rho_diexcorg - Generate all double excitations of key_in using the bit masks of holes and - particles. - Assume N_int is already provided. - - -h_apply_select_mono_di_delta_rho_diexcp - Undocumented - - -h_apply_select_mono_di_delta_rho_monoexc - Generate all single excitations of key_in using the bit masks of holes and - particles. - Assume N_int is already provided. - - -`micro_pt2 `_ - Helper program to compute the PT2 in distributed mode. - - -`provide_everything `_ - Undocumented - - -`run_wf `_ - Undocumented - - `var_pt2_ratio_run `_ Undocumented diff --git a/plugins/Full_CI_ZMQ/README.rst b/plugins/Full_CI_ZMQ/README.rst new file mode 100644 index 00000000..d1677a7d --- /dev/null +++ b/plugins/Full_CI_ZMQ/README.rst @@ -0,0 +1,461 @@ +Needed Modules +============== +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. + + +.. image:: tree_dependency.png + +* `Perturbation `_ +* `Selectors_full `_ +* `Generators_full `_ +* `ZMQ `_ +* `Full_CI `_ + +Documentation +============= +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. + + +`add_task_to_taskserver `_ + Get a task from the task server + + +`add_to_selection_buffer `_ + Undocumented + + +`assert `_ + Undocumented + + +`connect_to_taskserver `_ + Connect to the task server and obtain the worker ID + + +`create_selection_buffer `_ + Undocumented + + +`disconnect_from_taskserver `_ + Disconnect from the task server + + +`end_parallel_job `_ + End a new parallel job with name 'name'. The slave tasks execute subroutine 'slave' + + +`end_zmq_pair_socket `_ + Terminate socket on which the results are sent. + + +`end_zmq_pull_socket `_ + Terminate socket on which the results are sent. + + +`end_zmq_push_socket `_ + Terminate socket on which the results are sent. + + +`end_zmq_sub_socket `_ + Terminate socket on which the results are sent. + + +`end_zmq_to_qp_run_socket `_ + Terminate the socket from the application to qp_run + + +`fci_zmq `_ + Undocumented + + +`fill_buffer_double `_ + Undocumented + + +`fill_buffer_single `_ + Undocumented + + +`full_ci `_ + Undocumented + + +`get_d0 `_ + Undocumented + + +`get_d1 `_ + Undocumented + + +`get_d2 `_ + Undocumented + + +`get_m0 `_ + Undocumented + + +`get_m1 `_ + Undocumented + + +`get_m2 `_ + Undocumented + + +`get_mask_phase `_ + Undocumented + + +`get_phase_bi `_ + Undocumented + + +`get_task_from_taskserver `_ + Get a task from the task server + + +h_apply_fci + Calls H_apply on the HF determinant and selects all connected single and double + excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script. + + +h_apply_fci_diexc + Undocumented + + +h_apply_fci_diexcorg + Generate all double excitations of key_in using the bit masks of holes and + particles. + Assume N_int is already provided. + + +h_apply_fci_diexcp + Undocumented + + +h_apply_fci_mono + Calls H_apply on the HF determinant and selects all connected single and double + excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script. + + +h_apply_fci_mono_diexc + Undocumented + + +h_apply_fci_mono_diexcorg + Generate all double excitations of key_in using the bit masks of holes and + particles. + Assume N_int is already provided. + + +h_apply_fci_mono_diexcp + Undocumented + + +h_apply_fci_mono_monoexc + Generate all single excitations of key_in using the bit masks of holes and + particles. + Assume N_int is already provided. + + +h_apply_fci_monoexc + Generate all single excitations of key_in using the bit masks of holes and + particles. + Assume N_int is already provided. + + +h_apply_fci_no_selection + Calls H_apply on the HF determinant and selects all connected single and double + excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script. + + +h_apply_fci_no_selection_diexc + Undocumented + + +h_apply_fci_no_selection_diexcorg + Generate all double excitations of key_in using the bit masks of holes and + particles. + Assume N_int is already provided. + + +h_apply_fci_no_selection_diexcp + Undocumented + + +h_apply_fci_no_selection_monoexc + Generate all single excitations of key_in using the bit masks of holes and + particles. + Assume N_int is already provided. + + +h_apply_fci_no_skip + Calls H_apply on the HF determinant and selects all connected single and double + excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script. + + +h_apply_fci_no_skip_diexc + Undocumented + + +h_apply_fci_no_skip_diexcorg + Generate all double excitations of key_in using the bit masks of holes and + particles. + Assume N_int is already provided. + + +h_apply_fci_no_skip_diexcp + Undocumented + + +h_apply_fci_no_skip_monoexc + Generate all single excitations of key_in using the bit masks of holes and + particles. + Assume N_int is already provided. + + +h_apply_fci_pt2 + Calls H_apply on the HF determinant and selects all connected single and double + excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script. + + +h_apply_fci_pt2_collector + Collects results from the selection in an array of generators + + +h_apply_fci_pt2_diexc + Undocumented + + +h_apply_fci_pt2_diexcorg + Generate all double excitations of key_in using the bit masks of holes and + particles. + Assume N_int is already provided. + + +h_apply_fci_pt2_diexcp + Undocumented + + +h_apply_fci_pt2_monoexc + Generate all single excitations of key_in using the bit masks of holes and + particles. + Assume N_int is already provided. + + +h_apply_fci_pt2_slave + Calls H_apply on the HF determinant and selects all connected single and double + excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script. + + +h_apply_fci_pt2_slave_inproc + Computes a buffer using threads + + +h_apply_fci_pt2_slave_tcp + Computes a buffer over the network + + +`integral8 `_ + Undocumented + + +`new_parallel_job `_ + Start a new parallel job with name 'name'. The slave tasks execute subroutine 'slave' + + +`new_zmq_pair_socket `_ + Socket on which the collector and the main communicate + + +`new_zmq_pull_socket `_ + Socket on which the results are sent. If thread is 1, use inproc + + +`new_zmq_push_socket `_ + Socket on which the results are sent. If thread is 1, use inproc + + +`new_zmq_sub_socket `_ + Socket to read the state published by the Task server + + +`new_zmq_to_qp_run_socket `_ + Socket on which the qp_run process replies + + +`past_d1 `_ + Undocumented + + +`past_d2 `_ + Undocumented + + +`provide_everything `_ + Undocumented + + +`psi_phasemask `_ + Undocumented + + +`pull_selection_results `_ + Undocumented + + +`push_selection_results `_ + Undocumented + + +`qp_run_address `_ + Address of the qp_run socket + Example : tcp://130.120.229.139:12345 + + +`reset_zmq_addresses `_ + Socket which pulls the results (2) + + +`run_selection_slave `_ + Undocumented + + +`run_wf `_ + Undocumented + + +`select_connected `_ + Undocumented + + +`select_doubles `_ + Undocumented + + +`select_singles `_ + Select determinants connected to i_det by H + + +`selection_collector `_ + Undocumented + + +`selection_slave `_ + Helper program to compute the PT2 in distributed mode. + + +`selection_slave_inproc `_ + Undocumented + + +`selection_slave_tcp `_ + Undocumented + + +`sort_selection_buffer `_ + Undocumented + + +`splash_p `_ + Undocumented + + +`splash_pq `_ + Undocumented + + +`spot_hasbeen `_ + Undocumented + + +`spot_isinwf `_ + Undocumented + + +`switch_qp_run_to_master `_ + Address of the master qp_run socket + Example : tcp://130.120.229.139:12345 + + +`task_done_to_taskserver `_ + Get a task from the task server + + +`update_energy `_ + Update energy when it is received from ZMQ + + +`var_pt2_ratio_run `_ + Undocumented + + +`wait_for_next_state `_ + Undocumented + + +`wait_for_state `_ + Wait for the ZMQ state to be ready + + +`wait_for_states `_ + Wait for the ZMQ state to be ready + + +`zmq_context `_ + Context for the ZeroMQ library + + +`zmq_delete_task `_ + When a task is done, it has to be removed from the list of tasks on the qp_run + queue. This guarantees that the results have been received in the pull. + + +`zmq_port `_ + Return the value of the ZMQ port from the corresponding integer + + +`zmq_port_start `_ + Address of the qp_run socket + Example : tcp://130.120.229.139:12345 + + +`zmq_selection `_ + Undocumented + + +`zmq_set_running `_ + Set the job to Running in QP-run + + +`zmq_socket_pair_inproc_address `_ + Socket which pulls the results (2) + + +`zmq_socket_pull_inproc_address `_ + Socket which pulls the results (2) + + +`zmq_socket_pull_tcp_address `_ + Socket which pulls the results (2) + + +`zmq_socket_push_inproc_address `_ + Socket which pulls the results (2) + + +`zmq_socket_push_tcp_address `_ + Socket which pulls the results (2) + + +`zmq_socket_sub_tcp_address `_ + Socket which pulls the results (2) + + +`zmq_state `_ + Threads executing work through the ZeroMQ interface + diff --git a/plugins/Full_CI_ZMQ/fci_zmq.irp.f b/plugins/Full_CI_ZMQ/fci_zmq.irp.f index ae0d7989..087a5579 100644 --- a/plugins/Full_CI_ZMQ/fci_zmq.irp.f +++ b/plugins/Full_CI_ZMQ/fci_zmq.irp.f @@ -96,11 +96,14 @@ program fci_zmq if(do_pt2_end)then print*,'Last iteration only to compute the PT2' - threshold_selectors = max(threshold_selectors,threshold_selectors_pt2) - threshold_generators = max(threshold_generators,threshold_generators_pt2) - TOUCH threshold_selectors threshold_generators + !threshold_selectors = max(threshold_selectors,threshold_selectors_pt2) + !threshold_generators = max(threshold_generators,threshold_generators_pt2) + !TOUCH threshold_selectors threshold_generators + threshold_selectors = 1.d0 + threshold_generators = 1d0 E_CI_before(1:N_states) = CI_energy(1:N_states) - call ZMQ_selection(0, pt2) + !call ZMQ_selection(0, pt2)! pour non-stochastic + call ZMQ_pt2(pt2) print *, 'Final step' print *, 'N_det = ', N_det print *, 'N_states = ', N_states @@ -119,6 +122,145 @@ program fci_zmq end +subroutine ZMQ_pt2(pt2) + use f77_zmq + use selection_types + + implicit none + + character*(512) :: task + integer(ZMQ_PTR) :: zmq_to_qp_run_socket + type(selection_buffer) :: b + integer, external :: omp_get_thread_num + double precision, intent(out) :: pt2(N_states) + + + double precision, allocatable :: pt2_detail(:,:), comb(:) + logical, allocatable :: computed(:) + integer, allocatable :: tbc(:) + integer :: i, j, Ncomb, generator_per_task, i_generator_end + integer, external :: pt2_find + + double precision :: sumabove(comb_teeth), sum2above(comb_teeth), Nabove(comb_teeth) + double precision, external :: omp_get_wtime + double precision :: time0, time + + allocate(pt2_detail(N_states, N_det_generators), comb(100000), computed(N_det_generators), tbc(0:N_det_generators)) + provide nproc + + !call random_seed() + + computed = .false. + tbc(0) = first_det_of_comb - 1 + do i=1, tbc(0) + tbc(i) = i + computed(i) = .true. + end do + + pt2_detail = 0d0 + + time0 = omp_get_wtime() + print *, "grep - time - avg - err - n_combs" + do while(.true.) + + call new_parallel_job(zmq_to_qp_run_socket,"pt2") + call zmq_put_psi(zmq_to_qp_run_socket,1,ci_electronic_energy,size(ci_electronic_energy)) + call zmq_set_running(zmq_to_qp_run_socket) + call create_selection_buffer(1, 1*2, b) + + + + + + + call get_carlo_workbatch(1d-3, computed, comb, Ncomb, tbc) + generator_per_task = 1 ! tbc(0)/300 + 1 + print *, 'TASKS REVERSED' + !do i=1,tbc(0),generator_per_task + do i=1,tbc(0) ! generator_per_task + i_generator_end = min(i+generator_per_task-1, tbc(0)) + !print *, "TASK", (i_generator_end-i+1), tbc(i:i_generator_end) + if(i > 10) then + integer :: zero + zero = 0 + write(task,*) (i_generator_end-i+1), zero, tbc(i:i_generator_end) + call add_task_to_taskserver(zmq_to_qp_run_socket,task) + else + do j=1,8 + write(task,*) (i_generator_end-i+1), j, tbc(i:i_generator_end) + call add_task_to_taskserver(zmq_to_qp_run_socket,task) + end do + end if + end do + + print *, "tasked" + !$OMP PARALLEL DEFAULT(shared) SHARED(b, pt2) PRIVATE(i) NUM_THREADS(nproc+1) + i = omp_get_thread_num() + if (i==0) then + call pt2_collector(b, pt2_detail) + else + call pt2_slave_inproc(i) + endif + !$OMP END PARALLEL + call end_parallel_job(zmq_to_qp_run_socket, 'pt2') + double precision :: E0, avg, eqt + call do_carlo(tbc, Ncomb, comb, pt2_detail, sumabove, sum2above, Nabove) + !END LOOP? + integer :: tooth + call get_first_tooth(computed, tooth) + !print *, "TOOTH ", tooth + + !!! ASSERT + !do i=1,first_det_of_teeth(tooth) + ! if(not(computed(i))) stop "deter non calc" + !end do + !logical :: ok + !ok = .false. + !do i=first_det_of_teeth(tooth), first_det_of_teeth(tooth+1) + ! if(not(computed(i))) ok = .true. + !end do + !if(not(ok)) stop "not OK..." + !!!!! + double precision :: prop + if(Nabove(tooth) >= 30) then + E0 = sum(pt2_detail(1,:first_det_of_teeth(tooth)-1)) + prop = ((1d0 - dfloat(comb_teeth - tooth + 1) * comb_step) - cweight(first_det_of_teeth(tooth)-1)) + prop = prop / weight(first_det_of_teeth(tooth)) + E0 += pt2_detail(1,first_det_of_teeth(tooth)) * prop + avg = E0 + (sumabove(tooth) / Nabove(tooth)) + eqt = sqrt(1d0 / (Nabove(tooth)-1) * abs(sum2above(tooth) / Nabove(tooth) - (sumabove(tooth)/Nabove(tooth))**2)) + time = omp_get_wtime() + print "(A, 5(E15.7))", "PT2stoch ", time - time0, avg, eqt, Nabove(tooth) + else + print *, Nabove(tooth), "< 30 combs" + end if + tbc(0) = 0 + end do + + pt2 = 0d0 +end subroutine + + +subroutine do_carlo(tbc, Ncomb, comb, pt2_detail, sumabove, sum2above, Nabove) + integer, intent(in) :: tbc(0:N_det_generators), Ncomb + double precision, intent(in) :: comb(Ncomb), pt2_detail(N_states, N_det_generators) + double precision, intent(inout) :: sumabove(comb_teeth), sum2above(comb_teeth), Nabove(comb_teeth) + integer :: i, dets(comb_teeth) + double precision :: myVal, myVal2 + + + do i=1,Ncomb + call get_comb(comb(i), dets) + myVal = 0d0 + myVal2 = 0d0 + do j=comb_teeth,1,-1 + myVal += pt2_detail(1, dets(j)) / weight(dets(j)) * comb_step + sumabove(j) += myVal + sum2above(j) += myVal**2 + Nabove(j) += 1 + end do + end do +end subroutine subroutine ZMQ_selection(N_in, pt2) @@ -136,43 +278,38 @@ subroutine ZMQ_selection(N_in, pt2) double precision, intent(out) :: pt2(N_states) - if (.True.) then - PROVIDE pt2_e0_denominator - N = max(N_in,1) - provide nproc - call new_parallel_job(zmq_to_qp_run_socket,"selection") - call zmq_put_psi(zmq_to_qp_run_socket,1,pt2_e0_denominator,size(pt2_e0_denominator)) - call zmq_set_running(zmq_to_qp_run_socket) - call create_selection_buffer(N, N*2, b) - endif + N = max(N_in,1) + provide nproc + provide ci_electronic_energy + call new_parallel_job(zmq_to_qp_run_socket,"selection") + call zmq_put_psi(zmq_to_qp_run_socket,1,ci_electronic_energy,size(ci_electronic_energy)) + call zmq_set_running(zmq_to_qp_run_socket) + call create_selection_buffer(N, N*2, b) integer :: i_generator, i_generator_start, i_generator_max, step ! step = int(max(1.,10*elec_num/mo_tot_num) step = int(5000000.d0 / dble(N_int * N_states * elec_num * elec_num * mo_tot_num * mo_tot_num )) step = max(1,step) - do i= 1, N_det_generators,step - i_generator_start = i - i_generator_max = min(i+step-1,N_det_generators) + do i= N_det_generators, 1, -step + i_generator_start = max(i-step+1,1) + i_generator_max = i write(task,*) i_generator_start, i_generator_max, 1, N call add_task_to_taskserver(zmq_to_qp_run_socket,task) end do - !$OMP PARALLEL DEFAULT(shared) SHARED(b, pt2) PRIVATE(i) NUM_THREADS(nproc+1) - i = omp_get_thread_num() - if (i==0) then - call selection_collector(b, pt2) - else - call selection_slave_inproc(i) - endif + !$OMP PARALLEL DEFAULT(none) SHARED(b, pt2) PRIVATE(i) NUM_THREADS(nproc+1) shared(ci_electronic_energy_is_built, n_det_generators_is_built, n_states_is_built, n_int_is_built, nproc_is_built) + i = omp_get_thread_num() + if (i==0) then + call selection_collector(b, pt2) + else + call selection_slave_inproc(i) + endif !$OMP END PARALLEL - call end_parallel_job(zmq_to_qp_run_socket, 'selection') + call end_parallel_job(zmq_to_qp_run_socket, 'selection') if (N_in > 0) then call fill_H_apply_buffer_no_selection(b%cur,b%det,N_int,0) !!! PAS DE ROBIN call copy_H_apply_buffer_to_wf() - if (s2_eig) then - call make_s2_eigenfunction - endif endif end subroutine @@ -181,9 +318,75 @@ subroutine selection_slave_inproc(i) implicit none integer, intent(in) :: i - call run_selection_slave(1,i,pt2_e0_denominator) + call run_selection_slave(1,i,ci_electronic_energy) end +subroutine pt2_slave_inproc(i) + implicit none + integer, intent(in) :: i + + call run_pt2_slave(1,i,ci_electronic_energy) +end + +subroutine pt2_collector(b, pt2_detail) + use f77_zmq + use selection_types + use bitmasks + implicit none + + + type(selection_buffer), intent(inout) :: b + double precision, intent(inout) :: pt2_detail(N_states, N_det) + double precision :: pt2_mwen(N_states, N_det) + integer(ZMQ_PTR),external :: new_zmq_to_qp_run_socket + integer(ZMQ_PTR) :: zmq_to_qp_run_socket + + integer(ZMQ_PTR), external :: new_zmq_pull_socket + integer(ZMQ_PTR) :: zmq_socket_pull + + integer :: msg_size, rc, more + integer :: acc, i, j, robin, N, ntask + double precision, allocatable :: val(:) + integer(bit_kind), allocatable :: det(:,:,:) + integer, allocatable :: task_id(:) + integer :: done, Nindex + integer, allocatable :: index(:) + real :: time, time0 + zmq_to_qp_run_socket = new_zmq_to_qp_run_socket() + zmq_socket_pull = new_zmq_pull_socket() + allocate(val(b%N), det(N_int, 2, b%N), task_id(N_det), index(N_det)) + done = 0 + more = 1 + !pt2_detail = 0d0 + call CPU_TIME(time0) + do while (more == 1) + call pull_pt2_results(zmq_socket_pull, Nindex, index, pt2_mwen, task_id, ntask) + do i=1,Nindex + pt2_detail(:, index(i)) += pt2_mwen(:,i) + end do + + !do i=1, N + ! call add_to_selection_buffer(b, det(1,1,i), val(i)) + !end do + + do i=1, ntask + if(task_id(i) == 0) then + print *, "Error in collector" + endif + call zmq_delete_task(zmq_to_qp_run_socket,zmq_socket_pull,task_id(i),more) + end do + done += ntask + call CPU_TIME(time) +! print *, "DONE" , done, time - time0 + end do + + + call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket) + call end_zmq_pull_socket(zmq_socket_pull) + call sort_selection_buffer(b) +end subroutine + + subroutine selection_collector(b, pt2) use f77_zmq use selection_types @@ -238,3 +441,194 @@ subroutine selection_collector(b, pt2) call sort_selection_buffer(b) end subroutine + + +integer function pt2_find(v, w) + implicit none + double precision :: v, w(N_det) + integer :: i,l,h + + l = 0 + h = N_det-1 + + do while(h >= l) + i = (h+l)/2 + if(w(i+1) > v) then + h = i-1 + else + l = i+1 + end if + end do + pt2_find = l+1 +end function + + +BEGIN_PROVIDER [ integer, comb_teeth ] + implicit none + comb_teeth = 100 +END_PROVIDER + + + +subroutine get_first_tooth(computed, first_teeth) + implicit none + logical, intent(in) :: computed(N_det_generators) + integer, intent(out) :: first_teeth + integer :: i + + first_teeth = 1 + do i=first_det_of_comb, N_det_generators + if(not(computed(i))) then + first_teeth = i + exit + end if + end do + + do i=comb_teeth, 1, -1 + if(first_det_of_teeth(i) < first_teeth) then + first_teeth = i + exit + end if + end do +end subroutine + + +subroutine get_carlo_workbatch(maxWorkload, computed, comb, Ncomb, tbc) + implicit none + double precision, intent(in) :: maxWorkload + double precision, intent(out) :: comb(N_det_generators) + integer, intent(inout) :: tbc(0:N_det_generators) + integer, intent(out) :: Ncomb + logical, intent(inout) :: computed(N_det_generators) + integer :: i, dets(comb_teeth) + double precision :: myWorkload + + myWorkload = 0d0 + + do i=1,size(comb) + call RANDOM_NUMBER(comb(i)) + comb(i) = comb(i) * comb_step + call add_comb(comb(i), computed, tbc, myWorkload) + Ncomb = i + if(myWorkload > maxWorkload .and. i >= 50) exit + end do + call reorder_tbc(tbc) +end subroutine + + +subroutine reorder_tbc(tbc) + implicit none + integer, intent(inout) :: tbc(0:N_det_generators) + logical, allocatable :: ltbc(:) + integer :: i, ci + + allocate(ltbc(N_det_generators)) + ltbc = .false. + do i=1,tbc(0) + ltbc(tbc(i)) = .true. + end do + + ci = 0 + do i=1,N_det_generators + if(ltbc(i)) then + ci += 1 + tbc(ci) = i + end if + end do +end subroutine + + +subroutine get_comb(stato, dets) + implicit none + double precision, intent(in) :: stato + integer, intent(out) :: dets(comb_teeth) + double precision :: curs + integer :: j + integer, external :: pt2_find + + curs = 1d0 - stato + do j = comb_teeth, 1, -1 + dets(j) = pt2_find(curs, cweight) + curs -= comb_step + end do +end subroutine + + +subroutine add_comb(comb, computed, tbc, workload) + implicit none + double precision, intent(in) :: comb + logical, intent(inout) :: computed(N_det_generators) + double precision, intent(inout) :: workload + integer, intent(inout) :: tbc(0:N_det_generators) + integer :: i, dets(comb_teeth) + + call get_comb(comb, dets) + + do i = 1, comb_teeth + if(not(computed(dets(i)))) then + tbc(0) += 1 + tbc(tbc(0)) = dets(i) + workload += comb_workload(dets(i)) + computed(dets(i)) = .true. + end if + end do +end subroutine + + + + BEGIN_PROVIDER [ double precision, weight, (N_det_generators) ] +&BEGIN_PROVIDER [ double precision, cweight, (N_det_generators) ] +&BEGIN_PROVIDER [ double precision, comb_workload, (N_det_generators) ] +&BEGIN_PROVIDER [ double precision, comb_step ] +&BEGIN_PROVIDER [ integer, first_det_of_teeth, (comb_teeth+1) ] +&BEGIN_PROVIDER [ integer, first_det_of_comb ] + implicit none + integer :: i + double precision :: norm_left, stato + integer, external :: pt2_find + + weight(1) = psi_coef_generators(1,1)**2 + cweight(1) = psi_coef_generators(1,1)**2 + + do i=2,N_det_generators + weight(i) = psi_coef_generators(i,1)**2 + cweight(i) = cweight(i-1) + psi_coef_generators(i,1)**2 + end do + + weight = weight / cweight(N_det_generators) + cweight = cweight / cweight(N_det_generators) + comb_workload = 1d0 / dfloat(N_det_generators) + + norm_left = 1d0 + + comb_step = 1d0/dfloat(comb_teeth) + do i=1,N_det_generators + if(weight(i)/norm_left < comb_step/2d0) then + first_det_of_comb = i + exit + end if + norm_left -= weight(i) + end do + + comb_step = 1d0 / dfloat(comb_teeth) * (1d0 - cweight(first_det_of_comb-1)) + + stato = 1d0 - comb_step! + 1d-5 + do i=comb_teeth, 1, -1 + first_det_of_teeth(i) = pt2_find(stato, cweight) + stato -= comb_step + end do + first_det_of_teeth(comb_teeth+1) = N_det_generators + 1 + first_det_of_teeth(1) = first_det_of_comb + if(first_det_of_teeth(1) /= first_det_of_comb) stop "comb provider" + +END_PROVIDER + + + + + + + + + + diff --git a/plugins/Full_CI_ZMQ/pt2_slave.irp.f b/plugins/Full_CI_ZMQ/pt2_slave.irp.f new file mode 100644 index 00000000..91c3db63 --- /dev/null +++ b/plugins/Full_CI_ZMQ/pt2_slave.irp.f @@ -0,0 +1,93 @@ +program pt2_slave + implicit none + BEGIN_DOC +! Helper program to compute the PT2 in distributed mode. + END_DOC + + read_wf = .False. + SOFT_TOUCH read_wf + call provide_everything + call switch_qp_run_to_master + call run_wf +end + +subroutine provide_everything + PROVIDE H_apply_buffer_allocated mo_bielec_integrals_in_map psi_det_generators psi_coef_generators psi_det_sorted_bit psi_selectors n_det_generators n_states generators_bitmask zmq_context +! PROVIDE ci_electronic_energy mo_tot_num N_int +end + +subroutine run_wf + use f77_zmq + implicit none + + integer(ZMQ_PTR), external :: new_zmq_to_qp_run_socket + integer(ZMQ_PTR) :: zmq_to_qp_run_socket + double precision :: energy(N_states_diag) + character*(64) :: states(1) + integer :: rc, i + + call provide_everything + + zmq_context = f77_zmq_ctx_new () + states(1) = 'pt2' + + zmq_to_qp_run_socket = new_zmq_to_qp_run_socket() + + do + + call wait_for_states(states,zmq_state,1) + + if(trim(zmq_state) == 'Stopped') then + + exit + + else if (trim(zmq_state) == 'pt2') then + + ! Selection + ! --------- + + print *, 'PT2' + call zmq_get_psi(zmq_to_qp_run_socket,1,energy,N_states_diag) + + !$OMP PARALLEL PRIVATE(i) + i = omp_get_thread_num() + call pt2_slave_tcp(i, energy) + !$OMP END PARALLEL + print *, 'PT2 done' + + endif + + end do +end + +subroutine update_energy(energy) + implicit none + double precision, intent(in) :: energy(N_states_diag) + BEGIN_DOC +! Update energy when it is received from ZMQ + END_DOC + integer :: j,k + do j=1,N_states + do k=1,N_det + CI_eigenvectors(k,j) = psi_coef(k,j) + enddo + enddo + call u_0_S2_u_0(CI_eigenvectors_s2,CI_eigenvectors,N_det,psi_det,N_int) + if (.True.) then + do k=1,size(ci_electronic_energy) + ci_electronic_energy(k) = energy(k) + enddo + TOUCH ci_electronic_energy CI_eigenvectors_s2 CI_eigenvectors + endif + + call write_double(6,ci_energy,'Energy') +end + +subroutine pt2_slave_tcp(i,energy) + implicit none + double precision, intent(in) :: energy(N_states_diag) + integer, intent(in) :: i + + call run_pt2_slave(0,i,energy) +end + diff --git a/plugins/Full_CI_ZMQ/run_pt2_slave.irp.f b/plugins/Full_CI_ZMQ/run_pt2_slave.irp.f new file mode 100644 index 00000000..949a6d28 --- /dev/null +++ b/plugins/Full_CI_ZMQ/run_pt2_slave.irp.f @@ -0,0 +1,168 @@ + +subroutine run_pt2_slave(thread,iproc,energy) + use f77_zmq + use selection_types + implicit none + + double precision, intent(in) :: energy(N_states_diag) + integer, intent(in) :: thread, iproc + integer :: rc, i + + integer :: worker_id, task_id(1), ctask, ltask + character*(1000000) :: task + + integer(ZMQ_PTR),external :: new_zmq_to_qp_run_socket + integer(ZMQ_PTR) :: zmq_to_qp_run_socket + + integer(ZMQ_PTR), external :: new_zmq_push_socket + integer(ZMQ_PTR) :: zmq_socket_push + + type(selection_buffer) :: buf, buf2 + logical :: done + + double precision :: pt2(N_states) + double precision,allocatable :: pt2_detail(:,:) + integer,allocatable :: index(:) + integer :: Nindex + + allocate(pt2_detail(N_states, N_det), index(N_det)) + zmq_to_qp_run_socket = new_zmq_to_qp_run_socket() + zmq_socket_push = new_zmq_push_socket(thread) + call connect_to_taskserver(zmq_to_qp_run_socket,worker_id,thread) + if(worker_id == -1) then + print *, "WORKER -1" + call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket) + call end_zmq_push_socket(zmq_socket_push,thread) + return + end if + buf%N = 0 + ctask = 1 + pt2 = 0d0 + pt2_detail = 0d0 + do + call get_task_from_taskserver(zmq_to_qp_run_socket,worker_id, task_id(ctask), task) + + done = task_id(ctask) == 0 + if (done) then + ctask = ctask - 1 + else + integer :: i_generator, i_i_generator, N, subset + read (task,*) Nindex + read (task,*) Nindex, subset, index(:Nindex) + + !!!!! + N=1 + !!!!! + if(buf%N == 0) then + ! Only first time + call create_selection_buffer(N, N*2, buf) + call create_selection_buffer(N, N*3, buf2) + else + if(N /= buf%N) stop "N changed... wtf man??" + end if + do i_i_generator=1, Nindex + i_generator = index(i_i_generator) + call select_connected(i_generator,energy,pt2_detail(1, i_i_generator),buf,subset) + pt2(:) += pt2_detail(:, i_generator) + enddo + endif + + if(done .or. ctask == size(task_id)) then + if(buf%N == 0 .and. ctask > 0) stop "uninitialized selection_buffer" + do i=1, ctask + call task_done_to_taskserver(zmq_to_qp_run_socket,worker_id,task_id(i)) + end do + if(ctask > 0) then + call push_pt2_results(zmq_socket_push, Nindex, index, pt2_detail, task_id(1), ctask) + !print *, "pushed ", index(:Nindex) + do i=1,buf%cur + call add_to_selection_buffer(buf2, buf%det(1,1,i), buf%val(i)) + enddo + call sort_selection_buffer(buf2) + buf%mini = buf2%mini + pt2 = 0d0 + pt2_detail(:,:Nindex) = 0d0 + buf%cur = 0 + end if + ctask = 0 + end if + + if(done) exit + ctask = ctask + 1 + end do + call disconnect_from_taskserver(zmq_to_qp_run_socket,zmq_socket_push,worker_id) + call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket) + call end_zmq_push_socket(zmq_socket_push,thread) +end subroutine + + +subroutine push_pt2_results(zmq_socket_push, N, index, pt2_detail, task_id, ntask) + use f77_zmq + use selection_types + implicit none + + integer(ZMQ_PTR), intent(in) :: zmq_socket_push + double precision, intent(in) :: pt2_detail(N_states, N_det) + integer, intent(in) :: ntask, N, index(N), task_id(*) + integer :: rc + + + rc = f77_zmq_send( zmq_socket_push, N, 4, ZMQ_SNDMORE) + if(rc /= 4) stop "push" + + rc = f77_zmq_send( zmq_socket_push, index, 4*N, ZMQ_SNDMORE) + if(rc /= 4*N) stop "push" + + + rc = f77_zmq_send( zmq_socket_push, pt2_detail, 8*N_states*N, ZMQ_SNDMORE) + if(rc /= 8*N_states*N) stop "push" + + rc = f77_zmq_send( zmq_socket_push, ntask, 4, ZMQ_SNDMORE) + if(rc /= 4) stop "push" + + rc = f77_zmq_send( zmq_socket_push, task_id(1), ntask*4, 0) + if(rc /= 4*ntask) stop "push" + +! Activate is zmq_socket_push is a REQ +! rc = f77_zmq_recv( zmq_socket_push, task_id(1), ntask*4, 0) +end subroutine + + +subroutine pull_pt2_results(zmq_socket_pull, N, index, pt2_detail, task_id, ntask) + use f77_zmq + use selection_types + implicit none + integer(ZMQ_PTR), intent(in) :: zmq_socket_pull + double precision, intent(inout) :: pt2_detail(N_states, N_det) + integer, intent(out) :: index(N_det) + integer, intent(out) :: N, ntask, task_id(*) + integer :: rc, rn, i + + rc = f77_zmq_recv( zmq_socket_pull, N, 4, 0) + if(rc /= 4) stop "pull" + + rc = f77_zmq_recv( zmq_socket_pull, index, 4*N, 0) + if(rc /= 4*N) stop "pull" + + rc = f77_zmq_recv( zmq_socket_pull, pt2_detail, N_states*8*N, 0) + if(rc /= 8*N_states*N) stop "pull" + + rc = f77_zmq_recv( zmq_socket_pull, ntask, 4, 0) + if(rc /= 4) stop "pull" + + rc = f77_zmq_recv( zmq_socket_pull, task_id(1), ntask*4, 0) + if(rc /= 4*ntask) stop "pull" + +! Activate is zmq_socket_pull is a REP +! rc = f77_zmq_send( zmq_socket_pull, task_id(1), ntask*4, 0) +end subroutine + + +BEGIN_PROVIDER [ double precision, pt2_workload, (N_det) ] + integer :: i + do i=1,N_det + pt2_workload(:) = dfloat(N_det - i + 1)**2 + end do + pt2_workload = pt2_workload / sum(pt2_workload) +END_PROVIDER + diff --git a/plugins/Full_CI_ZMQ/run_selection_slave.irp.f b/plugins/Full_CI_ZMQ/run_selection_slave.irp.f index dfaee629..7d48e5c0 100644 --- a/plugins/Full_CI_ZMQ/run_selection_slave.irp.f +++ b/plugins/Full_CI_ZMQ/run_selection_slave.irp.f @@ -53,7 +53,7 @@ subroutine run_selection_slave(thread,iproc,energy) !print *, "psi_selectors_coef ", psi_selectors_coef(N_det_selectors-5:N_det_selectors, 1) !call debug_det(psi_selectors(1,1,N_det_selectors), N_int) do i_generator=i_generator_start,i_generator_max,step - call select_connected(i_generator,energy,pt2,buf) + call select_connected(i_generator,energy,pt2,buf,0) enddo endif diff --git a/plugins/Full_CI_ZMQ/selection.irp.f b/plugins/Full_CI_ZMQ/selection.irp.f index 3f351004..2f524022 100644 --- a/plugins/Full_CI_ZMQ/selection.irp.f +++ b/plugins/Full_CI_ZMQ/selection.irp.f @@ -67,11 +67,11 @@ subroutine get_mask_phase(det, phasemask) end subroutine -subroutine select_connected(i_generator,E0,pt2,b) +subroutine select_connected(i_generator,E0,pt2,b,subset) use bitmasks use selection_types implicit none - integer, intent(in) :: i_generator + integer, intent(in) :: i_generator, subset type(selection_buffer), intent(inout) :: b double precision, intent(inout) :: pt2(N_states) integer :: k,l @@ -90,8 +90,7 @@ subroutine select_connected(i_generator,E0,pt2,b) particle_mask(k,2) = iand(generators_bitmask(k,2,s_part,l), not(psi_det_generators(k,2,i_generator)) ) enddo - call select_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,pt2,b) - call select_singles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,pt2,b) + call select_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,pt2,b,subset) enddo end subroutine @@ -116,164 +115,6 @@ end subroutine -! Selection single -! ---------------- - -subroutine select_singles(i_gen,hole_mask,particle_mask,fock_diag_tmp,E0,pt2,buf) - use bitmasks - use selection_types - implicit none - BEGIN_DOC -! Select determinants connected to i_det by H - END_DOC - integer, intent(in) :: i_gen - integer(bit_kind), intent(in) :: hole_mask(N_int,2), particle_mask(N_int,2) - double precision, intent(in) :: fock_diag_tmp(mo_tot_num) - double precision, intent(in) :: E0(N_states) - double precision, intent(inout) :: pt2(N_states) - type(selection_buffer), intent(inout) :: buf - - double precision :: vect(N_states, mo_tot_num) - logical :: bannedOrb(mo_tot_num) - integer :: i, j, k - integer :: h1,h2,s1,s2,i1,i2,ib,sp - integer(bit_kind) :: hole(N_int,2), particle(N_int,2), mask(N_int, 2) - logical :: fullMatch, ok - - - do k=1,N_int - hole (k,1) = iand(psi_det_generators(k,1,i_gen), hole_mask(k,1)) - hole (k,2) = iand(psi_det_generators(k,2,i_gen), hole_mask(k,2)) - particle(k,1) = iand(not(psi_det_generators(k,1,i_gen)), particle_mask(k,1)) - particle(k,2) = iand(not(psi_det_generators(k,2,i_gen)), particle_mask(k,2)) - enddo - - ! Create lists of holes and particles - ! ----------------------------------- - - integer :: N_holes(2), N_particles(2) - integer :: hole_list(N_int*bit_kind_size,2) - integer :: particle_list(N_int*bit_kind_size,2) - - call bitstring_to_list_ab(hole , hole_list , N_holes , N_int) - call bitstring_to_list_ab(particle, particle_list, N_particles, N_int) - - do sp=1,2 - do i=1, N_holes(sp) - h1 = hole_list(i,sp) - call apply_hole(psi_det_generators(1,1,i_gen), sp, h1, mask, ok, N_int) - bannedOrb = .true. - do j=1,N_particles(sp) - bannedOrb(particle_list(j, sp)) = .false. - end do - call spot_hasBeen(mask, sp, psi_det_sorted, i_gen, N_det, bannedOrb, fullMatch) - if(fullMatch) cycle - vect = 0d0 - call splash_p(mask, sp, psi_selectors(1,1,i_gen), psi_phasemask(1,1,i_gen), psi_selectors_coef_transp(1,i_gen), N_det_selectors - i_gen + 1, bannedOrb, vect) - call fill_buffer_single(i_gen, sp, h1, bannedOrb, fock_diag_tmp, E0, pt2, vect, buf) - end do - enddo -end subroutine - - -subroutine fill_buffer_single(i_generator, sp, h1, bannedOrb, fock_diag_tmp, E0, pt2, vect, buf) - use bitmasks - use selection_types - implicit none - - integer, intent(in) :: i_generator, sp, h1 - double precision, intent(in) :: vect(N_states, mo_tot_num) - logical, intent(in) :: bannedOrb(mo_tot_num) - double precision, intent(in) :: fock_diag_tmp(mo_tot_num) - double precision, intent(in) :: E0(N_states) - double precision, intent(inout) :: pt2(N_states) - type(selection_buffer), intent(inout) :: buf - logical :: ok - integer :: s1, s2, p1, p2, ib, istate - integer(bit_kind) :: mask(N_int, 2), det(N_int, 2) - double precision :: e_pert, delta_E, val, Hii, max_e_pert, tmp - double precision, external :: diag_H_mat_elem_fock - - - call apply_hole(psi_det_generators(1,1,i_generator), sp, h1, mask, ok, N_int) - - do p1=1,mo_tot_num - if(bannedOrb(p1)) cycle - if(vect(1, p1) == 0d0) cycle - call apply_particle(mask, sp, p1, det, ok, N_int) - - - Hii = diag_H_mat_elem_fock(psi_det_generators(1,1,i_generator),det,fock_diag_tmp,N_int) - max_e_pert = 0d0 - - do istate=1,N_states - val = vect(istate, p1) + vect(istate, p1) - delta_E = E0(istate) - Hii - tmp = dsqrt(delta_E * delta_E + val * val) - if (delta_E < 0.d0) then - tmp = -tmp - endif - e_pert = 0.5d0 * ( tmp - delta_E) - pt2(istate) += e_pert - if(dabs(e_pert) > dabs(max_e_pert)) max_e_pert = e_pert - end do - - if(dabs(max_e_pert) > buf%mini) call add_to_selection_buffer(buf, det, max_e_pert) - end do -end subroutine - - -subroutine splash_p(mask, sp, det, phasemask, coefs, N_sel, bannedOrb, vect) - use bitmasks - implicit none - - integer(bit_kind),intent(in) :: mask(N_int, 2), det(N_int,2,N_sel) - integer(1), intent(in) :: phasemask(N_int*bit_kind_size, 2, N_sel) - double precision, intent(in) :: coefs(N_states, N_sel) - integer, intent(in) :: sp, N_sel - logical, intent(inout) :: bannedOrb(mo_tot_num) - double precision, intent(inout) :: vect(N_states, mo_tot_num) - - integer :: i, j, h(0:2,2), p(0:3,2), nt - integer(bit_kind) :: perMask(N_int, 2), mobMask(N_int, 2), negMask(N_int, 2) - - do i=1,N_int - negMask(i,1) = not(mask(i,1)) - negMask(i,2) = not(mask(i,2)) - end do - - do i=1, N_sel - nt = 0 - do j=1,N_int - mobMask(j,1) = iand(negMask(j,1), det(j,1,i)) - mobMask(j,2) = iand(negMask(j,2), det(j,2,i)) - nt += popcnt(mobMask(j, 1)) + popcnt(mobMask(j, 2)) - end do - - if(nt > 3) cycle - - do j=1,N_int - perMask(j,1) = iand(mask(j,1), not(det(j,1,i))) - perMask(j,2) = iand(mask(j,2), not(det(j,2,i))) - end do - - call bitstring_to_list(perMask(1,1), h(1,1), h(0,1), N_int) - call bitstring_to_list(perMask(1,2), h(1,2), h(0,2), N_int) - - call bitstring_to_list(mobMask(1,1), p(1,1), p(0,1), N_int) - call bitstring_to_list(mobMask(1,2), p(1,2), p(0,2), N_int) - - if(nt == 3) then - call get_m2(det(1,1,i), phasemask(1,1,i), bannedOrb, vect, mask, h, p, sp, coefs(1, i)) - else if(nt == 2) then - call get_m1(det(1,1,i), phasemask(1,1,i), bannedOrb, vect, mask, h, p, sp, coefs(1, i)) - else - call get_m0(det(1,1,i), phasemask(1,1,i), bannedOrb, vect, mask, h, p, sp, coefs(1, i)) - end if - end do -end subroutine - - subroutine get_m2(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs) use bitmasks implicit none @@ -420,67 +261,12 @@ subroutine get_m0(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs) end do end subroutine - -subroutine spot_hasBeen(mask, sp, det, i_gen, N, banned, fullMatch) - use bitmasks - implicit none - - integer(bit_kind),intent(in) :: mask(N_int, 2), det(N_int, 2, N) - integer, intent(in) :: i_gen, N, sp - logical, intent(inout) :: banned(mo_tot_num) - logical, intent(out) :: fullMatch - - - integer :: i, j, na, nb, list(3), nt - integer(bit_kind) :: myMask(N_int, 2), negMask(N_int, 2) - - fullMatch = .false. - - do i=1,N_int - negMask(i,1) = not(mask(i,1)) - negMask(i,2) = not(mask(i,2)) - end do - - genl : do i=1, N - nt = 0 - - do j=1, N_int - myMask(j, 1) = iand(det(j, 1, i), negMask(j, 1)) - myMask(j, 2) = iand(det(j, 2, i), negMask(j, 2)) - nt += popcnt(myMask(j, 1)) + popcnt(myMask(j, 2)) - end do - - if(nt > 3) cycle - - if(nt <= 2 .and. i < i_gen) then - fullMatch = .true. - return - end if - - call bitstring_to_list(myMask(1,sp), list(1), na, N_int) - - if(nt == 3 .and. i < i_gen) then - do j=1,na - banned(list(j)) = .true. - end do - else if(nt == 1 .and. na == 1) then - banned(list(1)) = .true. - end if - end do genl -end subroutine - - - - -! Selection double -! ---------------- - -subroutine select_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,pt2,buf) +subroutine select_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,pt2,buf,subset) use bitmasks use selection_types implicit none - integer, intent(in) :: i_generator + integer, intent(in) :: i_generator, subset integer(bit_kind), intent(in) :: hole_mask(N_int,2), particle_mask(N_int,2) double precision, intent(in) :: fock_diag_tmp(mo_tot_num) double precision, intent(in) :: E0(N_states) @@ -496,6 +282,13 @@ subroutine select_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,p integer,allocatable :: preinteresting(:), prefullinteresting(:), interesting(:), fullinteresting(:) integer(bit_kind), allocatable :: minilist(:, :, :), fullminilist(:, :, :) + logical :: monoAdo, monoBdo; + integer :: maskInd + maskInd = -1 + + monoAdo = .true. + monoBdo = .true. + allocate(minilist(N_int, 2, N_det_selectors), fullminilist(N_int, 2, N_det)) allocate(preinteresting(0:N_det_selectors), prefullinteresting(0:N_det), interesting(0:N_det_selectors), fullinteresting(0:N_det)) @@ -544,6 +337,19 @@ subroutine select_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,p do s1=1,2 do i1=N_holes(s1),1,-1 ! Generate low excitations first + !if(subset /= 0 .and. mod(maskInd, 10) /= (subset-1)) then + ! maskInd += 1 + ! cycle + !end if + maskInd += 1 + + + + + if(subset == 0 .or. mod(maskInd, 8) == (subset-1)) then + + + h1 = hole_list(i1,s1) call apply_hole(psi_det_generators(1,1,i_generator), s1,h1, pmask, ok, N_int) @@ -592,43 +398,65 @@ subroutine select_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,p end if end do + + + + + end if + + + do s2=s1,2 sp = s1 + if(s1 /= s2) sp = 3 ib = 1 if(s1 == s2) ib = i1+1 + monoAdo = .true. do i2=N_holes(s2),ib,-1 ! Generate low excitations first - - h2 = hole_list(i2,s2) - call apply_hole(pmask, s2,h2, mask, ok, N_int) - logical :: banned(mo_tot_num, mo_tot_num,2) logical :: bannedOrb(mo_tot_num, 2) - banned = .false. - - call spot_isinwf(mask, fullminilist, i_generator, fullinteresting(0), banned, fullMatch, fullinteresting) - - if(fullMatch) cycle - - bannedOrb(1:mo_tot_num, 1:2) = .true. - do s3=1,2 - do i=1,N_particles(s3) - bannedOrb(particle_list(i,s3), s3) = .false. + if(subset == 0 .or. mod(maskInd, 8) == (subset-1)) then + h2 = hole_list(i2,s2) + call apply_hole(pmask, s2,h2, mask, ok, N_int) + banned = .false. + bannedOrb(1:mo_tot_num, 1:2) = .true. + do s3=1,2 + do i=1,N_particles(s3) + bannedOrb(particle_list(i,s3), s3) = .false. + enddo enddo - enddo - - mat = 0d0 - call splash_pq(mask, sp, minilist, i_generator, interesting(0), bannedOrb, banned, mat, interesting) - call fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2, mat, buf) + + if(s1 /= s2) then + if(monoBdo) then + bannedOrb(h1,s1) = .false. + end if + if(monoAdo) then + bannedOrb(h2,s2) = .false. + monoAdo = .false. + end if + end if + + + call spot_isinwf(mask, fullminilist, i_generator, fullinteresting(0), banned, fullMatch, fullinteresting) + if(fullMatch) cycle + + mat = 0d0 + call splash_pq(mask, sp, minilist, i_generator, interesting(0), bannedOrb, banned, mat, interesting) + + call fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2, mat, buf) + end if enddo + if(s1 /= s2) monoBdo = .false. enddo enddo enddo end subroutine + subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2, mat, buf) use bitmasks use selection_types diff --git a/plugins/Full_CI_ZMQ/target_pt2_ratio_zmq.irp.f b/plugins/Full_CI_ZMQ/target_pt2_ratio_zmq.irp.f new file mode 100644 index 00000000..77bbab03 --- /dev/null +++ b/plugins/Full_CI_ZMQ/target_pt2_ratio_zmq.irp.f @@ -0,0 +1,105 @@ +program fci_zmq + implicit none + integer :: i,j,k + logical, external :: detEq + + double precision, allocatable :: pt2(:) + integer :: Nmin, Nmax + integer :: n_det_before, to_select + double precision :: threshold_davidson_in, ratio, E_ref + + double precision, allocatable :: psi_coef_ref(:,:) + integer(bit_kind), allocatable :: psi_det_ref(:,:,:) + + + allocate (pt2(N_states)) + + pt2 = 1.d0 + threshold_davidson_in = threshold_davidson + threshold_davidson = threshold_davidson_in * 100.d0 + SOFT_TOUCH threshold_davidson + + ! Stopping criterion is the PT2max + + double precision :: E_CI_before(N_states) + do while (dabs(pt2(1)) > pt2_max) + print *, 'N_det = ', N_det + print *, 'N_states = ', N_states + do k=1, N_states + print*,'State ',k + print *, 'PT2 = ', pt2(k) + print *, 'E = ', CI_energy(k) + print *, 'E(before)+PT2 = ', E_CI_before(k)+pt2(k) + enddo + print *, '-----' + E_CI_before(1:N_states) = CI_energy(1:N_states) + call ezfio_set_full_ci_zmq_energy(CI_energy(1)) + + n_det_before = N_det + to_select = N_det + to_select = max(64-to_select, to_select) + call ZMQ_selection(to_select, pt2) + + PROVIDE psi_coef + PROVIDE psi_det + PROVIDE psi_det_sorted + + call diagonalize_CI + call save_wavefunction + call ezfio_set_full_ci_zmq_energy(CI_energy(1)) + enddo + + threshold_selectors = max(threshold_selectors,threshold_selectors_pt2) + threshold_generators = max(threshold_generators,threshold_generators_pt2) + threshold_davidson = threshold_davidson_in + TOUCH threshold_selectors threshold_generators threshold_davidson + call diagonalize_CI + call ZMQ_selection(0, pt2) + + E_ref = CI_energy(1) + pt2(1) + print *, 'Est FCI = ', E_ref + + Nmax = N_det + Nmin = 2 + allocate (psi_coef_ref(size(psi_coef_sorted,1),size(psi_coef_sorted,2))) + allocate (psi_det_ref(N_int,2,size(psi_det_sorted,3))) + psi_coef_ref = psi_coef_sorted + psi_det_ref = psi_det_sorted + psi_det = psi_det_sorted + psi_coef = psi_coef_sorted + TOUCH psi_coef psi_det + do while (Nmax-Nmin > 1) + psi_coef = psi_coef_ref + psi_det = psi_det_ref + TOUCH psi_det psi_coef + call diagonalize_CI + ratio = (CI_energy(1) - HF_energy) / (E_ref - HF_energy) + if (ratio < var_pt2_ratio) then + Nmin = N_det + else + Nmax = N_det + endif + N_det = Nmin + (Nmax-Nmin)/2 + print *, '-----' + print *, 'Det min, Det max: ', Nmin, Nmax + print *, 'Ratio : ', ratio, ' ~ ', var_pt2_ratio + print *, 'N_det = ', N_det + print *, 'E = ', CI_energy(1) + enddo + call ZMQ_selection(0, pt2) + print *, '------' + print *, 'HF_energy = ', HF_energy + print *, 'Est FCI = ', E_ref + print *, 'E = ', CI_energy(1) + print *, 'PT2 = ', pt2(1) + print *, 'E+PT2 = ', CI_energy(1)+pt2(1) + + E_CI_before(1:N_states) = CI_energy(1:N_states) + call save_wavefunction + call ezfio_set_full_ci_zmq_energy(CI_energy(1)) + call ezfio_set_full_ci_zmq_energy_pt2(E_CI_before(1)+pt2(1)) +end + + + + diff --git a/plugins/Full_CI_ZMQ/target_pt2_zmq.irp.f b/plugins/Full_CI_ZMQ/target_pt2_zmq.irp.f new file mode 100644 index 00000000..52f825f1 --- /dev/null +++ b/plugins/Full_CI_ZMQ/target_pt2_zmq.irp.f @@ -0,0 +1,95 @@ +program fci_zmq + implicit none + integer :: i,j,k + logical, external :: detEq + + double precision, allocatable :: pt2(:) + integer :: Nmin, Nmax + integer :: n_det_before, to_select + double precision :: threshold_davidson_in, ratio, E_ref, pt2_ratio + + allocate (pt2(N_states)) + + pt2 = 1.d0 + threshold_davidson_in = threshold_davidson + threshold_davidson = threshold_davidson_in * 100.d0 + SOFT_TOUCH threshold_davidson + + double precision :: E_CI_before(N_states) + do while (dabs(pt2(1)) > pt2_max) + print *, 'N_det = ', N_det + print *, 'N_states = ', N_states + do k=1, N_states + print*,'State ',k + print *, 'PT2 = ', pt2(k) + print *, 'E = ', CI_energy(k) + print *, 'E(before)+PT2 = ', E_CI_before(k)+pt2(k) + enddo + print *, '-----' + E_CI_before(1:N_states) = CI_energy(1:N_states) + call ezfio_set_full_ci_zmq_energy(CI_energy(1)) + + n_det_before = N_det + to_select = N_det + to_select = max(64-to_select, to_select) + call ZMQ_selection(to_select, pt2) + + PROVIDE psi_coef + PROVIDE psi_det + PROVIDE psi_det_sorted + + call diagonalize_CI + call save_wavefunction + call ezfio_set_full_ci_zmq_energy(CI_energy(1)) + enddo + + threshold_selectors = max(threshold_selectors,threshold_selectors_pt2) + threshold_generators = max(threshold_generators,threshold_generators_pt2) + threshold_davidson = threshold_davidson_in + TOUCH threshold_selectors threshold_generators threshold_davidson + call diagonalize_CI + call ZMQ_selection(0, pt2) + + E_ref = CI_energy(1) + pt2(1) + pt2_ratio = (E_ref + pt2_max - HF_energy) / (E_ref - HF_energy) + print *, 'Est FCI = ', E_ref + + Nmax = N_det + Nmin = N_det/8 + do while (Nmax-Nmin > 1) + call diagonalize_CI + ratio = (CI_energy(1) - HF_energy) / (E_ref - HF_energy) + psi_det = psi_det_sorted + psi_coef = psi_coef_sorted + TOUCH psi_coef psi_det + if (ratio < pt2_ratio) then + Nmin = N_det + to_select = (Nmax-Nmin)/2 + call ZMQ_selection(to_select, pt2) + else + Nmax = N_det + N_det = Nmin + (Nmax-Nmin)/2 + endif + print *, '-----' + print *, 'Det min, Det max: ', Nmin, Nmax + print *, 'Ratio : ', ratio, ' ~ ', pt2_ratio + print *, 'HF_energy = ', HF_energy + print *, 'Est FCI = ', E_ref + print *, 'N_det = ', N_det + print *, 'E = ', CI_energy(1) + print *, 'PT2 = ', pt2(1) + enddo + call ZMQ_selection(0, pt2) + print *, '------' + print *, 'E = ', CI_energy(1) + print *, 'PT2 = ', pt2(1) + + E_CI_before(1:N_states) = CI_energy(1:N_states) + call save_wavefunction + call ezfio_set_full_ci_zmq_energy(CI_energy(1)) + call ezfio_set_full_ci_zmq_energy_pt2(E_CI_before(1)+pt2(1)) +end + + + + diff --git a/plugins/Full_CI_ZMQ/tree_dependency.png b/plugins/Full_CI_ZMQ/tree_dependency.png new file mode 100644 index 00000000..e69de29b diff --git a/plugins/Full_CI_ZMQ/zmq_selection.irp.f b/plugins/Full_CI_ZMQ/zmq_selection.irp.f new file mode 100644 index 00000000..75992273 --- /dev/null +++ b/plugins/Full_CI_ZMQ/zmq_selection.irp.f @@ -0,0 +1,117 @@ +subroutine ZMQ_selection(N_in, pt2) + use f77_zmq + use selection_types + + implicit none + + character*(512) :: task + integer(ZMQ_PTR) :: zmq_to_qp_run_socket + integer, intent(in) :: N_in + type(selection_buffer) :: b + integer :: i, N + integer, external :: omp_get_thread_num + double precision, intent(out) :: pt2(N_states) + + + if (.True.) then + PROVIDE pt2_e0_denominator + N = max(N_in,1) + provide nproc + call new_parallel_job(zmq_to_qp_run_socket,"selection") + call zmq_put_psi(zmq_to_qp_run_socket,1,pt2_e0_denominator,size(pt2_e0_denominator)) + call zmq_set_running(zmq_to_qp_run_socket) + call create_selection_buffer(N, N*2, b) + endif + + integer :: i_generator, i_generator_start, i_generator_max, step +! step = int(max(1.,10*elec_num/mo_tot_num) + + step = int(5000000.d0 / dble(N_int * N_states * elec_num * elec_num * mo_tot_num * mo_tot_num )) + step = max(1,step) + do i= 1, N_det_generators,step + i_generator_start = i + i_generator_max = min(i+step-1,N_det_generators) + write(task,*) i_generator_start, i_generator_max, 1, N + call add_task_to_taskserver(zmq_to_qp_run_socket,task) + end do + + !$OMP PARALLEL DEFAULT(shared) SHARED(b, pt2) PRIVATE(i) NUM_THREADS(nproc+1) + i = omp_get_thread_num() + if (i==0) then + call selection_collector(b, pt2) + else + call selection_slave_inproc(i) + endif + !$OMP END PARALLEL + call end_parallel_job(zmq_to_qp_run_socket, 'selection') + if (N_in > 0) then + call fill_H_apply_buffer_no_selection(b%cur,b%det,N_int,0) !!! PAS DE ROBIN + call copy_H_apply_buffer_to_wf() + if (s2_eig) then + call make_s2_eigenfunction + endif + endif +end subroutine + + +subroutine selection_slave_inproc(i) + implicit none + integer, intent(in) :: i + + call run_selection_slave(1,i,pt2_e0_denominator) +end + +subroutine selection_collector(b, pt2) + use f77_zmq + use selection_types + use bitmasks + implicit none + + + type(selection_buffer), intent(inout) :: b + double precision, intent(out) :: pt2(N_states) + double precision :: pt2_mwen(N_states) + integer(ZMQ_PTR),external :: new_zmq_to_qp_run_socket + integer(ZMQ_PTR) :: zmq_to_qp_run_socket + + integer(ZMQ_PTR), external :: new_zmq_pull_socket + integer(ZMQ_PTR) :: zmq_socket_pull + + integer :: msg_size, rc, more + integer :: acc, i, j, robin, N, ntask + double precision, allocatable :: val(:) + integer(bit_kind), allocatable :: det(:,:,:) + integer, allocatable :: task_id(:) + integer :: done + real :: time, time0 + zmq_to_qp_run_socket = new_zmq_to_qp_run_socket() + zmq_socket_pull = new_zmq_pull_socket() + allocate(val(b%N), det(N_int, 2, b%N), task_id(N_det)) + done = 0 + more = 1 + pt2(:) = 0d0 + call CPU_TIME(time0) + do while (more == 1) + call pull_selection_results(zmq_socket_pull, pt2_mwen, val(1), det(1,1,1), N, task_id, ntask) + pt2 += pt2_mwen + do i=1, N + call add_to_selection_buffer(b, det(1,1,i), val(i)) + end do + + do i=1, ntask + if(task_id(i) == 0) then + print *, "Error in collector" + endif + call zmq_delete_task(zmq_to_qp_run_socket,zmq_socket_pull,task_id(i),more) + end do + done += ntask + call CPU_TIME(time) +! print *, "DONE" , done, time - time0 + end do + + + call end_zmq_to_qp_run_socket(zmq_to_qp_run_socket) + call end_zmq_pull_socket(zmq_socket_pull) + call sort_selection_buffer(b) +end subroutine + diff --git a/plugins/Generators_full/README.rst b/plugins/Generators_full/README.rst index c30193a2..d1fc68ec 100644 --- a/plugins/Generators_full/README.rst +++ b/plugins/Generators_full/README.rst @@ -33,7 +33,7 @@ Documentation .. by the `update_README.py` script. -`degree_max_generators `_ +`degree_max_generators `_ Max degree of excitation (respect to HF) of the generators @@ -52,10 +52,10 @@ Documentation Hartree-Fock determinant -`select_max `_ +`select_max `_ Memo to skip useless selectors -`size_select_max `_ +`size_select_max `_ Size of the select_max array diff --git a/plugins/Generators_restart/tree_dependency.png b/plugins/Generators_restart/tree_dependency.png new file mode 100644 index 00000000..e69de29b diff --git a/plugins/Hartree_Fock/README.rst b/plugins/Hartree_Fock/README.rst index 77521b94..2e329163 100644 --- a/plugins/Hartree_Fock/README.rst +++ b/plugins/Hartree_Fock/README.rst @@ -67,11 +67,11 @@ Documentation Alpha Fock matrix in AO basis set -`fock_matrix_alpha_mo `_ +`fock_matrix_alpha_mo `_ Fock matrix on the MO basis -`fock_matrix_ao `_ +`fock_matrix_ao `_ Fock matrix in AO basis set @@ -79,7 +79,7 @@ Documentation Alpha Fock matrix in AO basis set -`fock_matrix_beta_mo `_ +`fock_matrix_beta_mo `_ Fock matrix on the MO basis @@ -115,7 +115,7 @@ Documentation .br -`fock_mo_to_ao `_ +`fock_mo_to_ao `_ Undocumented @@ -135,7 +135,7 @@ Documentation S^-1 Beta density matrix in the AO basis x S^-1 -`hf_energy `_ +`hf_energy `_ Hartree-Fock energy diff --git a/plugins/MRCC_Utils/.gitignore b/plugins/MRCC_Utils/.gitignore index 4c65ce66..7a0dd517 100644 --- a/plugins/MRCC_Utils/.gitignore +++ b/plugins/MRCC_Utils/.gitignore @@ -3,6 +3,7 @@ .ninja_log AO_Basis Bitmask +Davidson Determinants Electrons Ezfio_files diff --git a/plugins/MRCC_Utils/README.rst b/plugins/MRCC_Utils/README.rst index 39b5684c..ae041734 100644 --- a/plugins/MRCC_Utils/README.rst +++ b/plugins/MRCC_Utils/README.rst @@ -36,11 +36,19 @@ Documentation Compute 1st dimension such that it is aligned for vectorization. -`apply_rotation `_ +`apply_hole_local `_ + Undocumented + + +`apply_particle_local `_ + Undocumented + + +`apply_rotation `_ Apply the rotation found by find_rotation -`approx_dble `_ +`approx_dble `_ Undocumented @@ -63,23 +71,23 @@ Documentation Binomial coefficients -`ci_eigenvectors_dressed `_ - Eigenvectors/values of the CI matrix +`ci_eigenvectors_dressed `_ + Eigenvectors/values of the dressed CI matrix -`ci_eigenvectors_s2_dressed `_ - Eigenvectors/values of the CI matrix +`ci_eigenvectors_s2_dressed `_ + Eigenvectors/values of the dressed CI matrix -`ci_electronic_energy_dressed `_ - Eigenvectors/values of the CI matrix +`ci_electronic_energy_dressed `_ + Eigenvectors/values of the dressed CI matrix -`ci_energy_dressed `_ +`ci_energy_dressed `_ N_states lowest eigenvalues of the dressed CI matrix -`davidson_diag_hjj_mrcc `_ +`davidson_diag_hjj_mrcc `_ Davidson diagonalization with specific diagonal elements of the H matrix .br H_jj : specific diagonal H matrix elements to diagonalize de Davidson @@ -95,12 +103,39 @@ Documentation .br N_st : Number of eigenstates .br + N_st_diag : Number of states in which H is diagonalized + .br iunit : Unit for the I/O .br Initial guess vectors are not necessarily orthonormal -`davidson_diag_mrcc `_ +`davidson_diag_hjj_sjj_mrcc `_ + Davidson diagonalization with specific diagonal elements of the H matrix + .br + H_jj : specific diagonal H matrix elements to diagonalize de Davidson + .br + S2_jj : specific diagonal S^2 matrix elements + .br + dets_in : bitmasks corresponding to determinants + .br + u_in : guess coefficients on the various states. Overwritten + on exit + .br + dim_in : leftmost dimension of u_in + .br + sze : Number of determinants + .br + N_st : Number of eigenstates + .br + N_st_diag : Number of states in which H is diagonalized. Assumed > sze + .br + iunit : Unit for the I/O + .br + Initial guess vectors are not necessarily orthonormal + + +`davidson_diag_mrcc `_ Davidson diagonalization. .br dets_in : bitmasks corresponding to determinants @@ -119,19 +154,38 @@ Documentation Initial guess vectors are not necessarily orthonormal -`dble_fact `_ +`davidson_diag_mrcc_hs2 `_ + Davidson diagonalization. + .br + dets_in : bitmasks corresponding to determinants + .br + u_in : guess coefficients on the various states. Overwritten + on exit + .br + dim_in : leftmost dimension of u_in + .br + sze : Number of determinants + .br + N_st : Number of eigenstates + .br + iunit : Unit number for the I/O + .br + Initial guess vectors are not necessarily orthonormal + + +`dble_fact `_ Undocumented -`dble_fact_even `_ +`dble_fact_even `_ n!! -`dble_fact_odd `_ +`dble_fact_odd `_ n!! -`dble_logfact `_ +`dble_logfact `_ n!! @@ -139,19 +193,23 @@ Documentation Undocumented -`delta_ii `_ - Dressing matrix in N_det basis +`dec_exc `_ + Undocumented -`delta_ij `_ - Dressing matrix in N_det basis - - -`diagonalize_ci_dressed `_ +`diagonalize_ci_dressed `_ Replace the coefficients of the CI states by the coefficients of the eigenstates of the CI matrix +`dij `_ + Undocumented + + +`dij_unique `_ + Undocumented + + `dset_order `_ array A has already been sorted, and iorder has contains the new order of elements of A. This subroutine changes the order of x to match the new order of A. @@ -170,10 +228,26 @@ Documentation contains the new order of the elements. +`dtranspose `_ + Transpose input matrix A into output matrix B + + `erf0 `_ Undocumented +`exc_inf `_ + Undocumented + + +`exccmp `_ + Undocumented + + +`exceq `_ + Undocumented + + `f_integral `_ function that calculates the following integral \int_{\-infty}^{+\infty} x^n \exp(-p x^2) dx @@ -183,19 +257,19 @@ Documentation n! -`fact_inv `_ +`fact_inv `_ 1/n! -`find_rotation `_ +`find_rotation `_ Find A.C = B -`find_triples_and_quadruples `_ +`find_triples_and_quadruples `_ Undocumented -`find_triples_and_quadruples_micro `_ +`find_triples_and_quadruples_micro `_ Undocumented @@ -221,7 +295,15 @@ Documentation Undocumented -`get_pseudo_inverse `_ +`get_dij `_ + Undocumented + + +`get_dij_index `_ + Undocumented + + +`get_pseudo_inverse `_ Find C = A^-1 @@ -306,11 +388,63 @@ h_apply_mrcc_pt2_monoexc Assume N_int is already provided. -`h_matrix_dressed `_ +h_apply_mrcepa_pt2 + Calls H_apply on the HF determinant and selects all connected single and double + excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script. + + +h_apply_mrcepa_pt2_collector + Collects results from the selection in an array of generators + + +h_apply_mrcepa_pt2_diexc + Undocumented + + +h_apply_mrcepa_pt2_diexcorg + Generate all double excitations of key_in using the bit masks of holes and + particles. + Assume N_int is already provided. + + +h_apply_mrcepa_pt2_diexcp + Undocumented + + +h_apply_mrcepa_pt2_monoexc + Generate all single excitations of key_in using the bit masks of holes and + particles. + Assume N_int is already provided. + + +h_apply_mrcepa_pt2_slave + Calls H_apply on the HF determinant and selects all connected single and double + excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script. + + +h_apply_mrcepa_pt2_slave_inproc + Computes a buffer using threads + + +h_apply_mrcepa_pt2_slave_tcp + Computes a buffer over the network + + +`h_matrix_dressed `_ Dressed H with Delta_ij -`h_u_0_mrcc `_ +`h_s2_u_0_mrcc_nstates `_ + Computes v_0 = H|u_0> and s_0 = S^2 |u_0> + .br + n : number of determinants + .br + H_jj : array of + .br + S2_jj : array of + + +`h_u_0_mrcc_nstates `_ Computes v_0 = H|u_0> .br n : number of determinants @@ -392,7 +526,15 @@ h_apply_mrcc_pt2_monoexc Hermite polynomial -`hij_mrcc `_ +`hh_exists `_ + Undocumented + + +`hh_shortcut `_ + Undocumented + + +`hij_mrcc `_ < ref | H | Non-ref > matrix @@ -523,7 +665,7 @@ h_apply_mrcc_pt2_monoexc to be in integer*8 format -`inv_int `_ +`inv_int `_ 1/i @@ -541,6 +683,10 @@ h_apply_mrcc_pt2_monoexc iradix should be -1 in input. +`is_generable `_ + Undocumented + + `iset_order `_ array A has already been sorted, and iorder has contains the new order of elements of A. This subroutine changes the order of x to match the new order of A. @@ -559,15 +705,19 @@ h_apply_mrcc_pt2_monoexc contains the new order of the elements. -`lambda_mrcc `_ +`lambda_mrcc `_ cm/ or perturbative 1/Delta_E(m) -`lambda_mrcc_pt2 `_ +`lambda_mrcc_kept `_ cm/ or perturbative 1/Delta_E(m) -`lapack_diag `_ +`lambda_mrcc_pt2 `_ + cm/ or perturbative 1/Delta_E(m) + + +`lapack_diag `_ Diagonalize matrix H .br H is untouched between input and ouptut @@ -578,7 +728,7 @@ h_apply_mrcc_pt2_monoexc .br -`lapack_diag_s2 `_ +`lapack_diag_s2 `_ Diagonalize matrix H .br H is untouched between input and ouptut @@ -589,7 +739,7 @@ h_apply_mrcc_pt2_monoexc .br -`lapack_diagd `_ +`lapack_diagd `_ Diagonalize matrix H .br H is untouched between input and ouptut @@ -600,7 +750,7 @@ h_apply_mrcc_pt2_monoexc .br -`lapack_partial_diag `_ +`lapack_partial_diag `_ Diagonalize matrix H .br H is untouched between input and ouptut @@ -611,19 +761,27 @@ h_apply_mrcc_pt2_monoexc .br -`logfact `_ +`logfact `_ n! -`lowercase `_ +`lowercase `_ Transform to lower case +`map_load_from_disk `_ + Undocumented + + +`map_save_to_disk `_ + Undocumented + + `mrcc_dress `_ Undocumented -`mrcc_iterations `_ +`mrmode `_ Undocumented @@ -632,12 +790,24 @@ h_apply_mrcc_pt2_monoexc D(t) =! D(t) +( B(t)*C(t)) -`normalize `_ +`n_ex_exists `_ + Undocumented + + +`n_hh_exists `_ + Undocumented + + +`n_pp_exists `_ + Undocumented + + +`normalize `_ Normalizes vector u u is expected to be aligned in memory. -`nproc `_ +`nproc `_ Number of current OpenMP threads @@ -659,7 +829,7 @@ h_apply_mrcc_pt2_monoexc .br -`ortho_lowdin `_ +`ortho_lowdin `_ Compute C_new=C_old.S^-1/2 orthogonalization. .br overlap : overlap matrix @@ -677,6 +847,19 @@ h_apply_mrcc_pt2_monoexc .br +`ortho_qr `_ + Orthogonalization using Q.R factorization + .br + A : matrix to orthogonalize + .br + LDA : leftmost dimension of A + .br + n : Number of rows of A + .br + m : Number of columns of A + .br + + `overlap_a_b_c `_ Undocumented @@ -707,6 +890,10 @@ h_apply_mrcc_pt2_monoexc Undocumented +`pp_exists `_ + Undocumented + + `progress_active `_ Current status for displaying progress bars. Global variable. @@ -727,6 +914,14 @@ h_apply_mrcc_pt2_monoexc Current status for displaying progress bars. Global variable. +`psi_non_ref_sorted `_ + Undocumented + + +`psi_non_ref_sorted_idx `_ + Undocumented + + `psi_ref_lock `_ Locks on ref determinants to fill delta_ij @@ -735,6 +930,10 @@ h_apply_mrcc_pt2_monoexc Recenter two polynomials +`rho_mrcc `_ + Undocumented + + `rint `_ .. math:: .br @@ -762,10 +961,6 @@ h_apply_mrcc_pt2_monoexc Undocumented -`run_mrcc `_ - Undocumented - - `run_progress `_ Display a progress bar with documentation of what is happening @@ -774,7 +969,15 @@ h_apply_mrcc_pt2_monoexc Undocumented -`set_generators_bitmasks_as_holes_and_particles `_ +`searchdet `_ + Undocumented + + +`searchexc `_ + Undocumented + + +`set_generators_bitmasks_as_holes_and_particles `_ Undocumented @@ -790,7 +993,7 @@ h_apply_mrcc_pt2_monoexc to be in integer*8 format -`set_zero_extra_diag `_ +`set_zero_extra_diag `_ Undocumented @@ -800,6 +1003,14 @@ h_apply_mrcc_pt2_monoexc contains the new order of the elements. +`sort_det `_ + Undocumented + + +`sort_exc `_ + Undocumented + + `start_progress `_ Starts the progress bar @@ -817,18 +1028,37 @@ h_apply_mrcc_pt2_monoexc .br -`u_dot_u `_ +`tamise_exc `_ + Uncodumented : TODO + + +`transpose `_ + Transpose input matrix A into output matrix B + + +`u_0_h_u_0_mrcc_nstates `_ + Computes e_0 = / + .br + n : number of determinants + .br + + +`u_dot_u `_ Compute -`u_dot_v `_ +`u_dot_v `_ Compute -`wall_time `_ +`unsortedsearchdet `_ + Undocumented + + +`wall_time `_ The equivalent of cpu_time, but for the wall time. -`write_git_log `_ +`write_git_log `_ Write the last git commit in file iunit. diff --git a/plugins/MRCC_Utils/mrcc_utils.irp.f b/plugins/MRCC_Utils/mrcc_utils.irp.f index 5a5fb656..d6b9cc79 100644 --- a/plugins/MRCC_Utils/mrcc_utils.irp.f +++ b/plugins/MRCC_Utils/mrcc_utils.irp.f @@ -33,6 +33,7 @@ END_PROVIDER if (ihpsi_current(k) == 0.d0) then ihpsi_current(k) = 1.d-32 endif +! lambda_mrcc(k,i) = psi_non_ref_coef(i,k)/ihpsi_current(k) lambda_mrcc(k,i) = min(-1.d-32,psi_non_ref_coef(i,k)/ihpsi_current(k) ) lambda_pert = 1.d0 / (psi_ref_energy_diagonalized(k)-hii) if (lambda_pert / lambda_mrcc(k,i) < 0.5d0) then @@ -77,19 +78,6 @@ BEGIN_PROVIDER [ double precision, hij_mrcc, (N_det_non_ref,N_det_ref) ] END_PROVIDER -! BEGIN_PROVIDER [ double precision, delta_ij, (N_states,N_det_non_ref,N_det_ref) ] -!&BEGIN_PROVIDER [ double precision, delta_ii, (N_states,N_det_ref) ] -! implicit none -! BEGIN_DOC -! ! Dressing matrix in N_det basis -! END_DOC -! integer :: i,j,m -! delta_ij = 0.d0 -! delta_ii = 0.d0 -! call H_apply_mrcc(delta_ij,delta_ii,N_states,N_det_non_ref,N_det_ref) -! -!END_PROVIDER - BEGIN_PROVIDER [ double precision, h_matrix_dressed, (N_det,N_det,N_states) ] implicit none @@ -174,43 +162,6 @@ END_PROVIDER call u_0_S2_u_0(CI_eigenvectors_s2_dressed,CI_eigenvectors_dressed,N_det,psi_det,N_int,& N_states_diag,size(CI_eigenvectors_dressed,1)) -! double precision :: u_dot_u -! double precision, allocatable :: h(:,:,:), s(:,:) -! allocate (h(N_states,N_states,N_states), s(N_states,N_states)) -! do i=1,N_states -! do j=1,N_states -! s(i,j) = u_dot_v(CI_eigenvectors_dressed(1,i),CI_eigenvectors_dressed(1,j),N_det) -! print *, 'S(',i,',',j,')', s(i,j) -! enddo -! enddo -! -! do i=1,N_states -! h(i,i) = CI_electronic_energy_dressed(i) -! do j=i+1,N_states -! h(j,i) = (CI_electronic_energy_dressed(j)-CI_electronic_energy_dressed(i)) * s(i,j) -! h(i,j) = -h(j,i) -! print *, 'h(',i,',',i,')', h(i,j) -! enddo -! print *, 'h(',i,',',i,')', h(i,i) -! enddo -! call lapack_diag(eigenvalues,eigenvectors, h,size(h,1),N_states) -! do i=1,N_states -! CI_electronic_energy_dressed(i) = eigenvalues(i) -! do j=1,N_states -! h(i,j) = eigenvectors(i,j) -! enddo -! enddo -! do k=1,N_states -! eigenvectors(1:N_det,k) = 0.d0 -! do i=1,N_states -! eigenvectors(1:N_det,k) += CI_eigenvectors_dressed(1:N_det,k) * h(k,i) -! enddo -! enddo -! deallocate(h,s) -! - -! call multi_state(CI_electronic_energy_dressed,CI_eigenvectors_dressed,size(CI_eigenvectors_dressed,1)) - deallocate (eigenvectors,eigenvalues) else if (diag_algorithm == "Lapack") then @@ -790,29 +741,8 @@ END_PROVIDER end do dIj_unique(1:size(X), s) = X(1:size(X)) -! double precision, external :: ddot -! if (ddot (size(X), dIj_unique, 1, X, 1) < 0.d0) then -! dIj_unique(1:size(X),s) = -X(1:size(X)) -! endif - enddo - ! Adjust phase of dIj_unique - -! double precision :: snorm -! X = 0.d0 -! snorm = 0.d0 -! do s=1,N_states -! norm = 0.d0 -! do i=1,N_det_non_ref -! norm = norm + psi_non_ref_coef(i,s)*psi_non_ref_coef(i,s) -! enddo -! norm = dsqrt(norm) -! X(1:size(X)) = X(1:size(X)) + dIj_unique(1:size(X),s) * norm -! snorm += norm -! enddo -! X = X/snorm - do s=1,N_states do a_coll=1,n_exc_active @@ -821,7 +751,6 @@ END_PROVIDER i = active_excitation_to_determinants_idx(j,a_coll) if (i==0) exit rho_mrcc(i,s) = rho_mrcc(i,s) + active_excitation_to_determinants_val(s,j,a_coll) * dIj_unique(a_col,s) -! rho_mrcc(i,s) = rho_mrcc(i,s) + active_excitation_to_determinants_val(s,j,a_coll) * X(a_col) enddo end do diff --git a/plugins/MRPT_Utils/ezfio_interface.irp.f b/plugins/MRPT_Utils/ezfio_interface.irp.f index 6bd8931d..18fb453e 100644 --- a/plugins/MRPT_Utils/ezfio_interface.irp.f +++ b/plugins/MRPT_Utils/ezfio_interface.irp.f @@ -1,6 +1,6 @@ ! DO NOT MODIFY BY HAND ! Created by $QP_ROOT/scripts/ezfio_interface/ei_handler.py -! from file /home/scemama/quantum_package/src/MRPT_Utils/EZFIO.cfg +! from file /home/garniron/quantum_package/src/MRPT_Utils/EZFIO.cfg BEGIN_PROVIDER [ logical, do_third_order_1h1p ] diff --git a/plugins/Perturbation/README.rst b/plugins/Perturbation/README.rst index 810a58e1..1657e079 100644 --- a/plugins/Perturbation/README.rst +++ b/plugins/Perturbation/README.rst @@ -88,6 +88,7 @@ Needed Modules * `Properties `_ * `Hartree_Fock `_ +* `Davidson `_ Documentation ============= @@ -107,13 +108,13 @@ Documentation Undocumented -perturb_buffer_by_mono_delta_rho_one_point - Applly pertubration ``delta_rho_one_point`` to the buffer of determinants generated in the H_apply +perturb_buffer_by_mono_dipole_moment_z + Applly pertubration ``dipole_moment_z`` to the buffer of determinants generated in the H_apply routine. -perturb_buffer_by_mono_dipole_moment_z - Applly pertubration ``dipole_moment_z`` to the buffer of determinants generated in the H_apply +perturb_buffer_by_mono_dummy + Applly pertubration ``dummy`` to the buffer of determinants generated in the H_apply routine. @@ -152,13 +153,13 @@ perturb_buffer_by_mono_moller_plesset routine. -perturb_buffer_delta_rho_one_point - Applly pertubration ``delta_rho_one_point`` to the buffer of determinants generated in the H_apply +perturb_buffer_dipole_moment_z + Applly pertubration ``dipole_moment_z`` to the buffer of determinants generated in the H_apply routine. -perturb_buffer_dipole_moment_z - Applly pertubration ``dipole_moment_z`` to the buffer of determinants generated in the H_apply +perturb_buffer_dummy + Applly pertubration ``dummy`` to the buffer of determinants generated in the H_apply routine. @@ -197,27 +198,6 @@ perturb_buffer_moller_plesset routine. -`pt2_delta_rho_one_point `_ - compute the perturbatibe contribution to the Integrated Spin density at z = z_one point of one determinant - .br - for the various n_st states, at various level of theory. - .br - c_pert(i) = /( - ) - .br - e_2_pert(i) = c_pert(i) * - .br - H_pert_diag(i) = c_pert(i)^2 * - .br - To get the contribution of the first order : - .br - = sum(over i) e_2_pert(i) - .br - To get the contribution of the diagonal elements of the second order : - .br - [ + + sum(over i) H_pert_diag(i) ] / [1. + sum(over i) c_pert(i) **2] - .br - - `pt2_dipole_moment_z `_ compute the perturbatibe contribution to the dipole moment of one determinant .br @@ -239,7 +219,11 @@ perturb_buffer_moller_plesset .br -`pt2_epstein_nesbet `_ +`pt2_dummy `_ + Dummy perturbation to add all connected determinants. + + +`pt2_epstein_nesbet `_ compute the standard Epstein-Nesbet perturbative first order coefficient and second order energetic contribution .br for the various N_st states. @@ -250,7 +234,7 @@ perturb_buffer_moller_plesset .br -`pt2_epstein_nesbet_2x2 `_ +`pt2_epstein_nesbet_2x2 `_ compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution .br for the various N_st states. @@ -261,7 +245,7 @@ perturb_buffer_moller_plesset .br -`pt2_epstein_nesbet_sc2 `_ +`pt2_epstein_nesbet_sc2 `_ compute the standard Epstein-Nesbet perturbative first order coefficient and second order energetic contribution .br for the various N_st states, but with the CISD_SC2 energies and coefficients @@ -272,7 +256,7 @@ perturb_buffer_moller_plesset .br -`pt2_epstein_nesbet_sc2_no_projected `_ +`pt2_epstein_nesbet_sc2_no_projected `_ compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution .br for the various N_st states, @@ -296,7 +280,7 @@ perturb_buffer_moller_plesset H_pert_diag = c_pert -`pt2_epstein_nesbet_sc2_projected `_ +`pt2_epstein_nesbet_sc2_projected `_ compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution .br for the various N_st states, @@ -331,12 +315,12 @@ perturb_buffer_moller_plesset .br -`pt2_max `_ +`pt2_max `_ The selection process stops when the largest PT2 (for all the state) is lower than pt2_max in absolute value -`pt2_moller_plesset `_ +`pt2_moller_plesset `_ compute the standard Moller-Plesset perturbative first order coefficient and second order energetic contribution .br for the various n_st states. @@ -368,7 +352,7 @@ perturb_buffer_moller_plesset Threshold to select determinants. Set by selection routines. -`var_pt2_ratio `_ +`var_pt2_ratio `_ The selection process stops when the energy ratio variational/(variational+PT2) is equal to var_pt2_ratio diff --git a/plugins/Psiref_CAS/.gitignore b/plugins/Psiref_CAS/.gitignore index 69ebdc69..d79d94d9 100644 --- a/plugins/Psiref_CAS/.gitignore +++ b/plugins/Psiref_CAS/.gitignore @@ -3,6 +3,7 @@ .ninja_log AO_Basis Bitmask +Davidson Determinants Electrons Ezfio_files diff --git a/plugins/Psiref_CAS/README.rst b/plugins/Psiref_CAS/README.rst index 5d511317..a217e36c 100644 --- a/plugins/Psiref_CAS/README.rst +++ b/plugins/Psiref_CAS/README.rst @@ -58,6 +58,7 @@ Needed Modules .. image:: tree_dependency.png * `Psiref_Utils `_ +* `Davidson `_ Documentation ============= diff --git a/plugins/Psiref_CAS/psi_ref.irp.f b/plugins/Psiref_CAS/psi_ref.irp.f index d3b6c28f..ab9e6943 100644 --- a/plugins/Psiref_CAS/psi_ref.irp.f +++ b/plugins/Psiref_CAS/psi_ref.irp.f @@ -67,3 +67,36 @@ END_PROVIDER END_PROVIDER + BEGIN_PROVIDER [double precision, norm_psi_ref, (N_states)] +&BEGIN_PROVIDER [double precision, inv_norm_psi_ref, (N_states)] + implicit none + integer :: i,j + norm_psi_ref = 0.d0 + do j = 1, N_states + do i = 1, N_det_ref + norm_psi_ref(j) += psi_ref_coef(i,j) * psi_ref_coef(i,j) + enddo + inv_norm_psi_ref(j) = 1.d0/(dsqrt(norm_psi_Ref(j))) + enddo + + END_PROVIDER + + BEGIN_PROVIDER [double precision, psi_ref_coef_interm_norm, (N_det_ref,N_states)] + implicit none + integer :: i,j + do j = 1, N_states + do i = 1, N_det_ref + psi_ref_coef_interm_norm(i,j) = inv_norm_psi_ref(j) * psi_ref_coef(i,j) + enddo + enddo + END_PROVIDER + + BEGIN_PROVIDER [double precision, psi_non_ref_coef_interm_norm, (N_det_non_ref,N_states)] + implicit none + integer :: i,j + do j = 1, N_states + do i = 1, N_det_non_ref + psi_non_ref_coef_interm_norm(i,j) = psi_non_ref_coef(i,j) * inv_norm_psi_ref(j) + enddo + enddo + END_PROVIDER diff --git a/plugins/Psiref_Utils/README.rst b/plugins/Psiref_Utils/README.rst index 35232d23..2ceb6b98 100644 --- a/plugins/Psiref_Utils/README.rst +++ b/plugins/Psiref_Utils/README.rst @@ -154,11 +154,11 @@ Documentation Compute 1st dimension such that it is aligned for vectorization. -`apply_rotation `_ +`apply_rotation `_ Apply the rotation found by find_rotation -`approx_dble `_ +`approx_dble `_ Undocumented @@ -181,19 +181,19 @@ Documentation Binomial coefficients -`dble_fact `_ +`dble_fact `_ Undocumented -`dble_fact_even `_ +`dble_fact_even `_ n!! -`dble_fact_odd `_ +`dble_fact_odd `_ n!! -`dble_logfact `_ +`dble_logfact `_ n!! @@ -219,6 +219,10 @@ Documentation contains the new order of the elements. +`dtranspose `_ + Transpose input matrix A into output matrix B + + `erf0 `_ Undocumented @@ -236,11 +240,11 @@ Documentation n! -`fact_inv `_ +`fact_inv `_ 1/n! -`find_rotation `_ +`find_rotation `_ Find A.C = B @@ -270,7 +274,7 @@ Documentation Returns the index of the determinant in the ``psi_ref_sorted_bit`` array -`get_pseudo_inverse `_ +`get_pseudo_inverse `_ Find C = A^-1 @@ -531,7 +535,7 @@ Documentation to be in integer*8 format -`inv_int `_ +`inv_int `_ 1/i @@ -571,7 +575,7 @@ Documentation contains the new order of the elements. -`lapack_diag `_ +`lapack_diag `_ Diagonalize matrix H .br H is untouched between input and ouptut @@ -582,7 +586,7 @@ Documentation .br -`lapack_diag_s2 `_ +`lapack_diag_s2 `_ Diagonalize matrix H .br H is untouched between input and ouptut @@ -593,7 +597,7 @@ Documentation .br -`lapack_diagd `_ +`lapack_diagd `_ Diagonalize matrix H .br H is untouched between input and ouptut @@ -604,7 +608,7 @@ Documentation .br -`lapack_partial_diag `_ +`lapack_partial_diag `_ Diagonalize matrix H .br H is untouched between input and ouptut @@ -615,14 +619,22 @@ Documentation .br -`logfact `_ +`logfact `_ n! -`lowercase `_ +`lowercase `_ Transform to lower case +`map_load_from_disk `_ + Undocumented + + +`map_save_to_disk `_ + Undocumented + + `multiply_poly `_ Multiply two polynomials D(t) =! D(t) +( B(t)*C(t)) @@ -635,12 +647,12 @@ Documentation idx_non_ref_rev gives the reverse. -`normalize `_ +`normalize `_ Normalizes vector u u is expected to be aligned in memory. -`nproc `_ +`nproc `_ Number of current OpenMP threads @@ -662,7 +674,7 @@ Documentation .br -`ortho_lowdin `_ +`ortho_lowdin `_ Compute C_new=C_old.S^-1/2 orthogonalization. .br overlap : overlap matrix @@ -680,6 +692,19 @@ Documentation .br +`ortho_qr `_ + Orthogonalization using Q.R factorization + .br + A : matrix to orthogonalize + .br + LDA : leftmost dimension of A + .br + n : Number of rows of A + .br + m : Number of columns of A + .br + + `overlap_a_b_c `_ Undocumented @@ -860,7 +885,7 @@ Documentation to be in integer*8 format -`set_zero_extra_diag `_ +`set_zero_extra_diag `_ Undocumented @@ -887,18 +912,22 @@ Documentation .br -`u_dot_u `_ +`transpose `_ + Transpose input matrix A into output matrix B + + +`u_dot_u `_ Compute -`u_dot_v `_ +`u_dot_v `_ Compute -`wall_time `_ +`wall_time `_ The equivalent of cpu_time, but for the wall time. -`write_git_log `_ +`write_git_log `_ Write the last git commit in file iunit. diff --git a/plugins/Psiref_threshold/psi_ref.irp.f b/plugins/Psiref_threshold/psi_ref.irp.f index ee69ef5c..62321140 100644 --- a/plugins/Psiref_threshold/psi_ref.irp.f +++ b/plugins/Psiref_threshold/psi_ref.irp.f @@ -1,5 +1,44 @@ use bitmasks +! BEGIN_PROVIDER [ integer(bit_kind), psi_ref, (N_int,2,psi_det_size) ] +!&BEGIN_PROVIDER [ double precision, psi_ref_coef, (psi_det_size,n_states) ] +!&BEGIN_PROVIDER [ integer, idx_ref, (psi_det_size) ] +!&BEGIN_PROVIDER [ integer, N_det_ref ] +! implicit none +! BEGIN_DOC +! ! Reference wave function, defined as determinants with amplitudes > 0.05 +! ! idx_ref gives the indice of the ref determinant in psi_det. +! END_DOC +! integer :: i, k, l +! logical :: good +! double precision, parameter :: threshold=0.01d0 +! double precision :: t(N_states) +! N_det_ref = 0 +! do l = 1, N_states +! t(l) = threshold * abs_psi_coef_max(l) +! enddo +! do i=1,N_det +! good = .False. +! do l=1, N_states +! psi_ref_coef(i,l) = 0.d0 +! good = good.or.(dabs(psi_coef(i,l)) > t(l)) +! enddo +! if (good) then +! N_det_ref = N_det_ref+1 +! do k=1,N_int +! psi_ref(k,1,N_det_ref) = psi_det(k,1,i) +! psi_ref(k,2,N_det_ref) = psi_det(k,2,i) +! enddo +! idx_ref(N_det_ref) = i +! do k=1,N_states +! psi_ref_coef(N_det_ref,k) = psi_coef(i,k) +! enddo +! endif +! enddo +! call write_int(output_determinants,N_det_ref, 'Number of determinants in the reference') +! +!END_PROVIDER + BEGIN_PROVIDER [ integer(bit_kind), psi_ref, (N_int,2,psi_det_size) ] &BEGIN_PROVIDER [ double precision, psi_ref_coef, (psi_det_size,n_states) ] &BEGIN_PROVIDER [ integer, idx_ref, (psi_det_size) ] @@ -10,30 +49,16 @@ use bitmasks ! idx_ref gives the indice of the ref determinant in psi_det. END_DOC integer :: i, k, l - logical :: good - double precision, parameter :: threshold=0.05d0 - double precision :: t(N_states) - N_det_ref = 0 - do l = 1, N_states - t(l) = threshold * abs_psi_coef_max(l) - enddo - do i=1,N_det - good = .False. - do l=1, N_states - psi_ref_coef(i,l) = 0.d0 - good = good.or.(dabs(psi_coef(i,l)) > t(l)) + double precision, parameter :: threshold=0.01d0 + + call find_reference(threshold, N_det_ref, idx_ref) + do l=1,N_states + do i=1,N_det_ref + psi_ref_coef(i,l) = psi_coef(idx_ref(i), l) enddo - if (good) then - N_det_ref = N_det_ref+1 - do k=1,N_int - psi_ref(k,1,N_det_ref) = psi_det(k,1,i) - psi_ref(k,2,N_det_ref) = psi_det(k,2,i) - enddo - idx_ref(N_det_ref) = i - do k=1,N_states - psi_ref_coef(N_det_ref,k) = psi_coef(i,k) - enddo - endif + enddo + do i=1,N_det_ref + psi_ref(:,:,i) = psi_det(:,:,idx_ref(i)) enddo call write_int(output_determinants,N_det_ref, 'Number of determinants in the reference') diff --git a/plugins/QmcChem/e_curve_qmc.irp.f b/plugins/QmcChem/e_curve_qmc.irp.f index d45624a0..169db84e 100644 --- a/plugins/QmcChem/e_curve_qmc.irp.f +++ b/plugins/QmcChem/e_curve_qmc.irp.f @@ -5,6 +5,8 @@ program e_curve double precision :: norm, E, hij, num, ci, cj integer, allocatable :: iorder(:) double precision , allocatable :: norm_sort(:) + PROVIDE mo_bielec_integrals_in_map + nab = n_det_alpha_unique+n_det_beta_unique allocate ( norm_sort(0:nab), iorder(0:nab) ) diff --git a/plugins/Selectors_full/README.rst b/plugins/Selectors_full/README.rst index 393e9421..fc264fc1 100644 --- a/plugins/Selectors_full/README.rst +++ b/plugins/Selectors_full/README.rst @@ -161,15 +161,19 @@ Documentation n_double_selectors = number of double excitations in the selectors determinants -`psi_selectors `_ +`psi_selectors `_ Determinants on which we apply for perturbation. -`psi_selectors_coef `_ +`psi_selectors_coef `_ Determinants on which we apply for perturbation. -`psi_selectors_diag_h_mat `_ +`psi_selectors_coef_transp `_ + Transposed psi_selectors + + +`psi_selectors_diag_h_mat `_ Diagonal elements of the H matrix for each selectors @@ -177,7 +181,7 @@ Documentation Undocumented -`zmq_get_psi `_ +`zmq_get_psi `_ Get the wave function from the qp_run scheduler diff --git a/plugins/Selectors_no_sorted/tree_dependency.png b/plugins/Selectors_no_sorted/tree_dependency.png new file mode 100644 index 00000000..e69de29b diff --git a/plugins/analyze_wf/NEEDED_CHILDREN_MODULES b/plugins/analyze_wf/NEEDED_CHILDREN_MODULES new file mode 100644 index 00000000..aae89501 --- /dev/null +++ b/plugins/analyze_wf/NEEDED_CHILDREN_MODULES @@ -0,0 +1 @@ +Determinants diff --git a/plugins/analyze_wf/README.rst b/plugins/analyze_wf/README.rst new file mode 100644 index 00000000..179e407d --- /dev/null +++ b/plugins/analyze_wf/README.rst @@ -0,0 +1,12 @@ +========== +analyze_wf +========== + +Needed Modules +============== +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. +Documentation +============= +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. diff --git a/plugins/analyze_wf/analyze_wf.irp.f b/plugins/analyze_wf/analyze_wf.irp.f new file mode 100644 index 00000000..6d8bffcf --- /dev/null +++ b/plugins/analyze_wf/analyze_wf.irp.f @@ -0,0 +1,70 @@ +program analyze_wf + implicit none + BEGIN_DOC +! Wave function analyzis + END_DOC + read_wf = .True. + SOFT_TOUCH read_wf + call run() +end + +subroutine run + implicit none + integer :: istate, i + integer :: class(0:mo_tot_num,5) + double precision :: occupation(mo_tot_num) + + write(*,'(A)') 'MO Occupation' + write(*,'(A)') '=============' + write(*,'(A)') '' + do istate=1,N_states + call get_occupation_from_dets(occupation,1) + write(*,'(A)') '' + write(*,'(A,I3)'), 'State ', istate + write(*,'(A)') '---------------' + write(*,'(A)') '' + write (*,'(A)') '======== ================' + class = 0 + do i=1,mo_tot_num + write (*,'(I8,X,F16.10)') i, occupation(i) + if (occupation(i) > 1.999d0) then + class(0,1) += 1 + class( class(0,1), 1) = i + else if (occupation(i) > 1.95d0) then + class(0,2) += 1 + class( class(0,2), 2) = i + else if (occupation(i) < 0.001d0) then + class(0,5) += 1 + class( class(0,5), 5) = i + else if (occupation(i) < 0.01d0) then + class(0,4) += 1 + class( class(0,4), 4) = i + else + class(0,3) += 1 + class( class(0,3), 3) = i + endif + enddo + write (*,'(A)') '======== ================' + write (*,'(A)') '' + + write (*,'(A)') 'Suggested classes' + write (*,'(A)') '-----------------' + write (*,'(A)') '' + write (*,'(A)') 'Core :' + write (*,*) (class(i,1), ',', i=1,class(0,1)) + write (*,*) '' + write (*,'(A)') 'Inactive :' + write (*,*) (class(i,2), ',', i=1,class(0,2)) + write (*,'(A)') '' + write (*,'(A)') 'Active :' + write (*,*) (class(i,3), ',', i=1,class(0,3)) + write (*,'(A)') '' + write (*,'(A)') 'Virtual :' + write (*,*) (class(i,4), ',', i=1,class(0,4)) + write (*,'(A)') '' + write (*,'(A)') 'Deleted :' + write (*,*) (class(i,5), ',', i=1,class(0,5)) + write (*,'(A)') '' + enddo + +end diff --git a/plugins/analyze_wf/occupation.irp.f b/plugins/analyze_wf/occupation.irp.f new file mode 100644 index 00000000..d426dc14 --- /dev/null +++ b/plugins/analyze_wf/occupation.irp.f @@ -0,0 +1,23 @@ +subroutine get_occupation_from_dets(occupation, istate) + implicit none + double precision, intent(out) :: occupation(mo_tot_num) + integer, intent(in) :: istate + BEGIN_DOC + ! Returns the average occupation of the MOs + END_DOC + integer :: i,j, ispin + integer :: list(N_int*bit_kind_size,2) + integer :: n_elements(2) + double precision :: c + + occupation = 0.d0 + do i=1,N_det + c = psi_coef(i,istate)*psi_coef(i,istate) + call bitstring_to_list_ab(psi_det(1,1,i), list, n_elements, N_int) + do ispin=1,2 + do j=1,n_elements(ispin) + occupation( list(j,ispin) ) += c + enddo + enddo + enddo +end diff --git a/plugins/mrcc_selected/EZFIO.cfg b/plugins/mrcc_selected/EZFIO.cfg new file mode 100644 index 00000000..b64637e6 --- /dev/null +++ b/plugins/mrcc_selected/EZFIO.cfg @@ -0,0 +1,33 @@ +[lambda_type] +type: Positive_int +doc: lambda type +interface: ezfio,provider,ocaml +default: 0 + +[energy] +type: double precision +doc: Calculated energy +interface: ezfio + +[energy_pt2] +type: double precision +doc: Calculated energy with PT2 contribution +interface: ezfio + +[energy] +type: double precision +doc: Calculated energy +interface: ezfio + +[thresh_dressed_ci] +type: Threshold +doc: Threshold on the convergence of the dressed CI energy +interface: ezfio,provider,ocaml +default: 1.e-5 + +[n_it_max_dressed_ci] +type: Strictly_positive_int +doc: Maximum number of dressed CI iterations +interface: ezfio,provider,ocaml +default: 10 + diff --git a/plugins/mrcc_selected/NEEDED_CHILDREN_MODULES b/plugins/mrcc_selected/NEEDED_CHILDREN_MODULES new file mode 100644 index 00000000..ea28c761 --- /dev/null +++ b/plugins/mrcc_selected/NEEDED_CHILDREN_MODULES @@ -0,0 +1 @@ +Perturbation Selectors_full Generators_full Psiref_threshold MRCC_Utils ZMQ diff --git a/plugins/mrcc_selected/README.rst b/plugins/mrcc_selected/README.rst new file mode 100644 index 00000000..997d005e --- /dev/null +++ b/plugins/mrcc_selected/README.rst @@ -0,0 +1,12 @@ +======= +mrcepa0 +======= + +Needed Modules +============== +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. +Documentation +============= +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. diff --git a/plugins/mrcc_selected/dressing.irp.f b/plugins/mrcc_selected/dressing.irp.f index c772e2aa..23fedcee 100644 --- a/plugins/mrcc_selected/dressing.irp.f +++ b/plugins/mrcc_selected/dressing.irp.f @@ -534,63 +534,9 @@ END_PROVIDER END_PROVIDER -! BEGIN_PROVIDER [ double precision, delta_cas, (N_det_ref, N_det_ref, N_states) ] -! use bitmasks -! implicit none -! integer :: i,j,k -! double precision :: Hjk, Hki, Hij, pre(N_det_ref), wall -! integer :: i_state, degree, npre, ipre(N_det_ref), npres(N_det_ref) -! -! ! provide lambda_mrcc -! npres = 0 -! delta_cas = 0d0 -! call wall_time(wall) -! print *, "dcas ", wall -! do i_state = 1, N_states -! !!$OMP PARALLEL DO default(none) schedule(dynamic) private(pre,npre,ipre,j,k,Hjk,Hki,degree) shared(npres,lambda_mrcc,i_state, N_det_non_ref,psi_ref, psi_non_ref,N_int,delta_cas,N_det_ref) -! do k=1,N_det_non_ref -! if(lambda_mrcc(i_state, k) == 0d0) cycle -! npre = 0 -! do i=1,N_det_ref -! call i_h_j(psi_non_ref(1,1,k),psi_ref(1,1,i), N_int,Hki) -! if(Hki /= 0d0) then -! !!$OMP ATOMIC -! npres(i) += 1 -! npre += 1 -! ipre(npre) = i -! pre(npre) = Hki -! end if -! end do -! -! -! do i=1,npre -! do j=1,i -! !!$OMP ATOMIC -! delta_cas(ipre(i),ipre(j),i_state) += pre(i) * pre(j) * lambda_mrcc(i_state, k) -! end do -! end do -! end do -! !!$OMP END PARALLEL DO -! npre=0 -! do i=1,N_det_ref -! npre += npres(i) -! end do -! !stop -! do i=1,N_det_ref -! do j=1,i -! delta_cas(j,i,i_state) = delta_cas(i,j,i_state) -! end do -! end do -! end do -! -! call wall_time(wall) -! print *, "dcas", wall -! ! stop -! END_PROVIDER - - BEGIN_PROVIDER [ double precision, delta_cas, (N_det_ref, N_det_ref, N_states) ] -&BEGIN_PROVIDER [ double precision, delta_cas_s2, (N_det_ref, N_det_ref, N_states) ] + BEGIN_PROVIDER [ double precision, delta_ref, (N_det_ref, N_det_ref, N_states) ] +&BEGIN_PROVIDER [ double precision, delta_ref_s2, (N_det_ref, N_det_ref, N_states) ] use bitmasks implicit none integer :: i,j,k @@ -600,22 +546,22 @@ END_PROVIDER provide lambda_mrcc dIj do i_state = 1, N_states - !$OMP PARALLEL DO default(none) schedule(dynamic) private(j,k,Sjk,Hjk,Hki,degree) shared(lambda_mrcc,i_state, N_det_non_ref,psi_ref, psi_non_ref,N_int,delta_cas,delta_cas_s2,N_det_ref,dij) + !$OMP PARALLEL DO default(none) schedule(dynamic) private(j,k,Sjk,Hjk,Hki,degree) shared(lambda_mrcc,i_state, N_det_non_ref,psi_ref, psi_non_ref,N_int,delta_ref,delta_ref_s2,N_det_ref,dij) do i=1,N_det_ref do j=1,i call get_excitation_degree(psi_ref(1,1,i), psi_ref(1,1,j), degree, N_int) - delta_cas(i,j,i_state) = 0d0 - delta_cas_s2(i,j,i_state) = 0d0 + delta_ref(i,j,i_state) = 0d0 + delta_ref_s2(i,j,i_state) = 0d0 do k=1,N_det_non_ref call i_h_j(psi_ref(1,1,j), psi_non_ref(1,1,k),N_int,Hjk) call get_s2(psi_ref(1,1,j), psi_non_ref(1,1,k),N_int,Sjk) - delta_cas(i,j,i_state) += Hjk * dij(i, k, i_state) ! * Hki * lambda_mrcc(i_state, k) - delta_cas_s2(i,j,i_state) += Sjk * dij(i, k, i_state) ! * Ski * lambda_mrcc(i_state, k) + delta_ref(i,j,i_state) += Hjk * dij(i, k, i_state) ! * Hki * lambda_mrcc(i_state, k) + delta_ref_s2(i,j,i_state) += Sjk * dij(i, k, i_state) ! * Ski * lambda_mrcc(i_state, k) end do - delta_cas(j,i,i_state) = delta_cas(i,j,i_state) - delta_cas_s2(j,i,i_state) = delta_cas_s2(i,j,i_state) + delta_ref(j,i,i_state) = delta_ref(i,j,i_state) + delta_ref_s2(j,i,i_state) = delta_ref_s2(i,j,i_state) end do end do !$OMP END PARALLEL DO @@ -739,7 +685,7 @@ end function !$OMP PARALLEL DO default(none) schedule(dynamic) shared(delta_mrcepa0_ij, delta_mrcepa0_ii, delta_mrcepa0_ij_s2, delta_mrcepa0_ii_s2) & !$OMP private(m,i,II,J,k,degree,myActive,made_hole,made_particle,hjk,contrib,contrib2,contrib_s2,contrib2_s2) & !$OMP shared(active_sorb, psi_non_ref, psi_non_ref_coef, psi_ref, psi_ref_coef, cepa0_shortcut, det_cepa0_active) & - !$OMP shared(N_det_ref, N_det_non_ref,N_int,det_cepa0_idx,lambda_mrcc,det_ref_active, delta_cas, delta_cas_s2) & + !$OMP shared(N_det_ref, N_det_non_ref,N_int,det_cepa0_idx,lambda_mrcc,det_ref_active, delta_ref, delta_ref_s2) & !$OMP shared(notf,i_state, sortRef, sortRefIdx, dij) do blok=1,cepa0_shortcut(0) do i=cepa0_shortcut(blok), cepa0_shortcut(blok+1)-1 @@ -781,8 +727,8 @@ end function notf = notf+1 ! call i_h_j(psi_non_ref(1,1,det_cepa0_idx(k)),psi_ref(1,1,J),N_int,HJk) - contrib = delta_cas(II, J, i_state) * dij(J, det_cepa0_idx(k), i_state) - contrib_s2 = delta_cas_s2(II, J, i_state) * dij(J, det_cepa0_idx(k), i_state) + contrib = delta_ref(II, J, i_state) * dij(J, det_cepa0_idx(k), i_state) + contrib_s2 = delta_ref_s2(II, J, i_state) * dij(J, det_cepa0_idx(k), i_state) if(dabs(psi_ref_coef(J,i_state)).ge.1.d-3) then contrib2 = contrib / psi_ref_coef(J, i_state) * psi_non_ref_coef(det_cepa0_idx(i),i_state) @@ -828,7 +774,7 @@ END_PROVIDER integer :: II, blok - provide delta_cas lambda_mrcc + provide delta_ref lambda_mrcc allocate(idx_sorted_bit(N_det)) idx_sorted_bit(:) = -1 do i=1,N_det_non_ref diff --git a/plugins/mrcc_selected/mrcc_selected.irp.f b/plugins/mrcc_selected/mrcc_selected.irp.f index 91592e62..b64f968d 100644 --- a/plugins/mrcc_selected/mrcc_selected.irp.f +++ b/plugins/mrcc_selected/mrcc_selected.irp.f @@ -8,7 +8,6 @@ program mrsc2sub read_wf = .True. SOFT_TOUCH read_wf - call print_cas_coefs call set_generators_bitmasks_as_holes_and_particles call run(N_states,energy) if(do_pt2_end)then diff --git a/plugins/mrcc_selected/mrcepa0_general.irp.f b/plugins/mrcc_selected/mrcepa0_general.irp.f index e3a2d1f5..812aeef0 100644 --- a/plugins/mrcc_selected/mrcepa0_general.irp.f +++ b/plugins/mrcc_selected/mrcepa0_general.irp.f @@ -60,16 +60,17 @@ subroutine run(N_st,energy) end -subroutine print_cas_coefs +subroutine print_ref_coefs implicit none integer :: i,j - print *, 'CAS' - print *, '===' - do i=1,N_det_cas - print *, (psi_cas_coef(i,j), j=1,N_states) - call debug_det(psi_cas(1,1,i),N_int) + print *, 'Reference' + print *, '=========' + do i=1,N_det_ref + print *, (psi_ref_coef(i,j), j=1,N_states) + call debug_det(psi_ref(1,1,i),N_int) enddo + print *, '' call write_double(6,ci_energy(1),"Initial CI energy") end @@ -202,7 +203,7 @@ subroutine run_pt2(N_st,energy) print*,'Last iteration only to compute the PT2' - N_det_generators = N_det_cas + N_det_generators = N_det_ref N_det_selectors = N_det_non_ref do i=1,N_det_generators diff --git a/plugins/mrcepa0/README.rst b/plugins/mrcepa0/README.rst index 997d005e..9e66ca0d 100644 --- a/plugins/mrcepa0/README.rst +++ b/plugins/mrcepa0/README.rst @@ -6,7 +6,203 @@ Needed Modules ============== .. Do not edit this section It was auto-generated .. by the `update_README.py` script. + + +.. image:: tree_dependency.png + +* `Perturbation `_ +* `Selectors_full `_ +* `Generators_full `_ +* `Psiref_CAS `_ +* `MRCC_Utils `_ +* `ZMQ `_ + Documentation ============= .. Do not edit this section It was auto-generated .. by the `update_README.py` script. + + +`active_sorb `_ + Undocumented + + +`blokmwen `_ + Undocumented + + +`cepa0_shortcut `_ + Undocumented + + +`child_num `_ + Undocumented + + +`delta_cas `_ + Undocumented + + +`delta_ii `_ + Undocumented + + +`delta_ii_mrcc `_ + Undocumented + + +`delta_ii_old `_ + Undocumented + + +`delta_ij `_ + Undocumented + + +`delta_ij_mrcc `_ + Undocumented + + +`delta_ij_old `_ + Undocumented + + +`delta_mrcepa0_ii `_ + Undocumented + + +`delta_mrcepa0_ij `_ + Undocumented + + +`delta_sub_ii `_ + Undocumented + + +`delta_sub_ij `_ + Undocumented + + +`det_cepa0 `_ + Undocumented + + +`det_cepa0_active `_ + Undocumented + + +`det_cepa0_idx `_ + Undocumented + + +`det_ref_active `_ + Undocumented + + +`filter_tq `_ + Undocumented + + +`filter_tq_micro `_ + Undocumented + + +`gethp `_ + Undocumented + + +`h_ `_ + Undocumented + + +`hp `_ + Undocumented + + +`isincassd `_ + Undocumented + + +`lambda_type `_ + lambda type + + +`linked `_ + Undocumented + + +`mrcc_part_dress `_ + Undocumented + + +`mrcepa0 `_ + Undocumented + + +`mrsc2 `_ + Undocumented + + +`mrsc2_dressing_collector `_ + Collects results from the AO integral calculation + + +`mrsc2_dressing_slave `_ + Task for parallel MR-SC2 + + +`mrsc2_dressing_slave_inproc `_ + Task for parallel MR-SC2 + + +`mrsc2_dressing_slave_tcp `_ + Task for parallel MR-SC2 + + +`mrsc2sub `_ + Undocumented + + +`n_it_max_dressed_ci `_ + Maximum number of dressed CI iterations + + +`nlink `_ + Undocumented + + +`print_cas_coefs `_ + Undocumented + + +`pull_mrsc2_results `_ + Push integrals in the push socket + + +`push_mrsc2_results `_ + Push integrals in the push socket + + +`run `_ + Undocumented + + +`run_pt2 `_ + Undocumented + + +`run_pt2_old `_ + Undocumented + + +`searchance `_ + Undocumented + + +`set_det_bit `_ + Undocumented + + +`thresh_dressed_ci `_ + Threshold on the convergence of the dressed CI energy + diff --git a/plugins/mrcepa0/tree_dependency.png b/plugins/mrcepa0/tree_dependency.png new file mode 100644 index 00000000..e69de29b diff --git a/plugins/mrsc2_no_amp/NEEDED_CHILDREN_MODULES b/plugins/mrsc2_no_amp/NEEDED_CHILDREN_MODULES new file mode 100644 index 00000000..f04fe3b0 --- /dev/null +++ b/plugins/mrsc2_no_amp/NEEDED_CHILDREN_MODULES @@ -0,0 +1 @@ +Psiref_CAS Determinants Davidson diff --git a/plugins/mrsc2_no_amp/README.rst b/plugins/mrsc2_no_amp/README.rst new file mode 100644 index 00000000..b24848f7 --- /dev/null +++ b/plugins/mrsc2_no_amp/README.rst @@ -0,0 +1,12 @@ +============ +mrsc2_no_amp +============ + +Needed Modules +============== +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. +Documentation +============= +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. diff --git a/plugins/mrsc2_no_amp/mrsc2_no_amp.irp.f b/plugins/mrsc2_no_amp/mrsc2_no_amp.irp.f new file mode 100644 index 00000000..b8b021e8 --- /dev/null +++ b/plugins/mrsc2_no_amp/mrsc2_no_amp.irp.f @@ -0,0 +1,78 @@ + BEGIN_PROVIDER [double precision, CI_eigenvectors_sc2_no_amp, (N_det,N_states_diag)] +&BEGIN_PROVIDER [double precision, CI_eigenvectors_s2_sc2_no_amp, (N_states_diag)] +&BEGIN_PROVIDER [double precision, CI_electronic_energy_sc2_no_amp, (N_states_diag)] + implicit none + integer :: i,j,k,l + integer, allocatable :: idx(:) + double precision, allocatable :: e_corr(:,:) + double precision, allocatable :: accu(:) + double precision, allocatable :: ihpsi_current(:) + double precision, allocatable :: H_jj(:),H_jj_total(:),S2_jj(:) + allocate(e_corr(N_det_non_ref,N_states),ihpsi_current(N_states),accu(N_states),H_jj(N_det_non_ref),idx(0:N_det_non_ref)) + allocate(H_jj_total(N_det),S2_jj(N_det)) + accu = 0.d0 + do i = 1, N_det_non_ref + call i_h_psi(psi_non_ref(1,1,i), psi_ref, psi_ref_coef_interm_norm, N_int, N_det_ref,& + size(psi_ref_coef_interm_norm,1), N_states,ihpsi_current) + do j = 1, N_states + e_corr(i,j) = psi_non_ref_coef_interm_norm(i,j) * ihpsi_current(j) + accu(j) += e_corr(i,j) + enddo + enddo + double precision :: hjj,diag_h_mat_elem + do i = 1, N_det_non_ref + call filter_not_connected(psi_non_ref,psi_non_ref(1,1,i),N_int,N_det_non_ref,idx) + H_jj(i) = 0.d0 + do j = 1, idx(0) + H_jj(i) += e_corr(idx(j),1) + enddo + enddo + do i=1,N_Det + H_jj_total(i) = diag_h_mat_elem(psi_det(1,1,i),N_int) + call get_s2(psi_det(1,1,i),psi_det(1,1,i),N_int,S2_jj(i)) + enddo + do i=1, N_det_non_ref + H_jj_total(idx_non_ref(i)) += H_jj(i) + enddo + + + call davidson_diag_hjj_sjj(psi_det,CI_eigenvectors_sc2_no_amp,H_jj_total,S2_jj,CI_electronic_energy_sc2_no_amp,size(CI_eigenvectors_sc2_no_amp,1),N_Det,N_states,N_states_diag,N_int,6) + do i=1,N_states_diag + CI_eigenvectors_s2_sc2_no_amp(i) = S2_jj(i) + enddo + + deallocate(e_corr,ihpsi_current,accu,H_jj,idx,H_jj_total,s2_jj) +END_PROVIDER + +BEGIN_PROVIDER [ double precision, CI_energy_sc2_no_amp, (N_states_diag) ] + implicit none + BEGIN_DOC + ! N_states lowest eigenvalues of the CI matrix + END_DOC + + integer :: j + character*(8) :: st + call write_time(output_determinants) + do j=1,min(N_det,N_states_diag) + CI_energy_sc2_no_amp(j) = CI_electronic_energy_sc2_no_amp(j) + nuclear_repulsion + enddo + do j=1,min(N_det,N_states) + write(st,'(I4)') j + call write_double(output_determinants,CI_energy_sc2_no_amp(j),'Energy of state '//trim(st)) + call write_double(output_determinants,CI_eigenvectors_s2_sc2_no_amp(j),'S^2 of state '//trim(st)) + enddo + +END_PROVIDER + +subroutine diagonalize_CI_sc2_no_amp + implicit none + integer :: i,j + do j=1,N_states + do i=1,N_det + psi_coef(i,j) = CI_eigenvectors_sc2_no_amp(i,j) + enddo + enddo + SOFT_TOUCH ci_eigenvectors_s2_sc2_no_amp ci_eigenvectors_sc2_no_amp ci_electronic_energy_sc2_no_amp ci_energy_sc2_no_amp psi_coef + +end + diff --git a/plugins/mrsc2_no_amp/sc2_no_amp.irp.f b/plugins/mrsc2_no_amp/sc2_no_amp.irp.f new file mode 100644 index 00000000..622d7236 --- /dev/null +++ b/plugins/mrsc2_no_amp/sc2_no_amp.irp.f @@ -0,0 +1,9 @@ +program pouet + implicit none + integer :: i + do i = 1, 10 + call diagonalize_CI_sc2_no_amp + TOUCH psi_coef + enddo + +end diff --git a/plugins/qmcpack/qp_convert_qmcpack_to_ezfio.py b/plugins/qmcpack/qp_convert_qmcpack_to_ezfio.py index e911af28..a1f47ccd 100755 --- a/plugins/qmcpack/qp_convert_qmcpack_to_ezfio.py +++ b/plugins/qmcpack/qp_convert_qmcpack_to_ezfio.py @@ -364,10 +364,6 @@ for line_raw in det_without_header.split("\n"): try: float(line) except ValueError: - - print line_raw.strip(), len(line_raw.strip()) - print l_order_mo, len(l_order_mo) - line_order = [line_raw[i] for i in l_order_mo] line= "".join([d_rep[x] if x in d_rep else x for x in line_raw]) diff --git a/scripts/ezfio_interface/qp_convert_output_to_ezfio.py b/scripts/ezfio_interface/qp_convert_output_to_ezfio.py index ff7ad225..1c541a21 100755 --- a/scripts/ezfio_interface/qp_convert_output_to_ezfio.py +++ b/scripts/ezfio_interface/qp_convert_output_to_ezfio.py @@ -12,12 +12,10 @@ Option: """ - import sys import os from functools import reduce - # ~#~#~#~#~#~#~#~ # # Add to the path # # ~#~#~#~#~#~#~#~ # @@ -29,9 +27,9 @@ except: print "Error: QP_ROOT environment variable not found." sys.exit(1) else: - sys.path = [QP_ROOT + "/install/EZFIO/Python", - QP_ROOT + "/resultsFile", - QP_ROOT + "/scripts"] + sys.path + sys.path = [ QP_ROOT + "/install/EZFIO/Python", + QP_ROOT + "/resultsFile", + QP_ROOT + "/scripts"] + sys.path # ~#~#~#~#~#~ # # I m p o r t # @@ -39,7 +37,6 @@ else: from ezfio import ezfio - try: from resultsFile import * except: @@ -254,7 +251,7 @@ def write_ezfio(res, filename): for coef in m.vector: MoMatrix.append(coef) - while len(MoMatrix) < len(MOs[0].vector) ** 2: + while len(MoMatrix) < len(MOs[0].vector)**2: MoMatrix.append(0.) # ~#~#~#~#~ # @@ -273,7 +270,130 @@ def write_ezfio(res, filename): # \_| |___/\___|\__,_|\__,_|\___/ # - ezfio.set_pseudo_do_pseudo(False) + # INPUT + # {% for lanel,zcore, l_block in l_atom $} + # #local l_block l=0} + # {label} GEN {zcore} {len(l_block)-1 #lmax_block} + # {% for l_param in l_block%} + # {len(l_param) # list of parameter aka n_max_bock_max(n)} + # {% for coef,n,zeta for l_param} + # {coef,n, zeta} + + # OUTPUT + + # Local are 1 array padded by max(n_max_block) when l == 0 (output:k_loc_max) + # v_k[n-2][atom] = value + + #No Local are 2 array padded with max of lmax_block when l!=0 (output:lmax+1) and max(n_max_block)whem l !=0 (kmax) + # v_kl[l][n-2][atom] = value + + def pad(array, size, value=0): + new_array = array + for add in xrange(len(array), size): + new_array.append(value) + + return new_array + + def parse_str(pseudo_str): + '''Return 4d array atom,l,n, attribute (attribute is coef, n, zeta)''' + matrix = [] + array_l_max_block = [] + array_z_remove = [] + + for block in [b for b in pseudo_str.split('\n\n') if b]: + #First element is header, the rest are l_param + array_party = [i for i in re.split(r"\n\d+\n", block) if i] + + z_remove, l_max_block = map(int, array_party[0].split()[-2:]) + array_l_max_block.append(l_max_block) + array_z_remove.append(z_remove) + + matrix.append([[coef_n_zeta.split()[1:] for coef_n_zeta in l.split('\n')] for l in array_party[1:]]) + + return (matrix, array_l_max_block, array_z_remove) + + def get_local_stuff(matrix): + + matrix_local_unpad = [atom[0] for atom in matrix] + k_loc_max = max(len(i) for i in matrix_local_unpad) + + matrix_local = [ pad(ll, k_loc_max, [0., 2, 0.]) for ll in matrix_local_unpad] + + m_coef = [[float(i[0]) for i in atom] for atom in matrix_local] + m_n = [[int(i[1]) - 2 for i in atom] for atom in matrix_local] + m_zeta = [[float(i[2]) for i in atom] for atom in matrix_local] + return (k_loc_max, m_coef, m_n, m_zeta) + + def get_non_local_stuff(matrix): + + matrix_unlocal_unpad = [atom[1:] for atom in matrix] + l_max_block = max(len(i) for i in matrix_unlocal_unpad) + k_max = max([len(item) for row in matrix_unlocal_unpad for item in row]) + + matrix_unlocal_semipaded = [[pad(item, k_max, [0., 2, 0.]) for item in row] for row in matrix_unlocal_unpad] + + empty_row = [[0., 2, 0.] for k in range(l_max_block)] + matrix_unlocal = [ pad(ll, l_max_block, empty_row) for ll in matrix_unlocal_semipaded ] + + m_coef_noloc = [[[float(k[0]) for k in j] for j in i] for i in matrix_unlocal] + m_n_noloc = [[[int(k[1]) - 2 for k in j] for j in i] for i in matrix_unlocal] + m_zeta_noloc = [[[float(k[2]) for k in j] for j in i] for i in matrix_unlocal] + + return (l_max_block, k_max, m_coef_noloc, m_n_noloc, m_zeta_noloc) + + try: + pseudo_str = res_file.get_pseudo() + matrix, array_l_max_block, array_z_remove = parse_str(pseudo_str) + except: + ezfio.set_pseudo_do_pseudo(False) + else: + ezfio.set_pseudo_do_pseudo(True) + + # ~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~ # + # Z _ e f f , a l p h a / b e t a _ e l e c # + # ~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~ # + + ezfio.pseudo_charge_remove = array_z_remove + ezfio.nuclei_nucl_charge = [ + i - j for i, j in zip(ezfio.nuclei_nucl_charge, array_z_remove) + ] + + import math + num_elec = sum(ezfio.nuclei_nucl_charge) + + ezfio.electrons_elec_alpha_num = int(math.ceil(num_elec / 2.)) + ezfio.electrons_elec_beta_num = int(math.floor(num_elec / 2.)) + + # Change all the array 'cause EZFIO + # v_kl (v, l) => v_kl(l,v) + # v_kl => zip(*_v_kl) + # [[7.0, 79.74474797, -49.45159098], [1.0, 5.41040609, -4.60151975]] + # [(7.0, 1.0), (79.74474797, 5.41040609), (-49.45159098, -4.60151975)] + + # ~#~#~#~#~ # + # L o c a l # + # ~#~#~#~#~ # + + klocmax, m_coef, m_n, m_zeta = get_local_stuff(matrix) + ezfio.pseudo_pseudo_klocmax = klocmax + + ezfio.pseudo_pseudo_v_k = zip(*m_coef) + ezfio.pseudo_pseudo_n_k = zip(*m_n) + ezfio.pseudo_pseudo_dz_k = zip(*m_zeta) + + # ~#~#~#~#~#~#~#~#~ # + # N o n _ L o c a l # + # ~#~#~#~#~#~#~#~#~ # + + l_max_block, k_max, m_coef_noloc, m_n_noloc, m_zeta_noloc = get_non_local_stuff( + matrix) + + ezfio.pseudo_pseudo_lmax = l_max_block - 1 + ezfio.pseudo_pseudo_kmax = k_max + + ezfio.pseudo_pseudo_v_kl = zip(*m_coef_noloc) + ezfio.pseudo_pseudo_n_kl = zip(*m_n_noloc) + ezfio.pseudo_pseudo_dz_kl = zip(*m_zeta_noloc) def get_full_path(file_path): @@ -282,6 +402,7 @@ def get_full_path(file_path): file_path = os.path.abspath(file_path) return file_path + if __name__ == '__main__': arguments = docopt(__doc__) diff --git a/src/AO_Basis/README.rst b/src/AO_Basis/README.rst index ae9acdf0..d67a3a63 100644 --- a/src/AO_Basis/README.rst +++ b/src/AO_Basis/README.rst @@ -133,7 +133,7 @@ Documentation :math:`\int \chi_i(r) \chi_j(r) dr)` -`ao_overlap_abs `_ +`ao_overlap_abs `_ Overlap between absolute value of atomic basis functions: :math:`\int |\chi_i(r)| |\chi_j(r)| dr)` diff --git a/src/Davidson/README.rst b/src/Davidson/README.rst new file mode 100644 index 00000000..15e9b46a --- /dev/null +++ b/src/Davidson/README.rst @@ -0,0 +1,322 @@ +Needed Modules +============== +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. + + +.. image:: tree_dependency.png + +* `Determinants `_ + +Documentation +============= +.. Do not edit this section It was auto-generated +.. by the `update_README.py` script. + + +`ci_eigenvectors `_ + Eigenvectors/values of the CI matrix + + +`ci_eigenvectors_mono `_ + Eigenvectors/values of the CI matrix + + +`ci_eigenvectors_s2 `_ + Eigenvectors/values of the CI matrix + + +`ci_eigenvectors_s2_mono `_ + Eigenvectors/values of the CI matrix + + +`ci_electronic_energy `_ + Eigenvectors/values of the CI matrix + + +`ci_electronic_energy_mono `_ + Eigenvectors/values of the CI matrix + + +`ci_energy `_ + N_states lowest eigenvalues of the CI matrix + + +`dav_det `_ + Temporary arrays for parallel davidson + .br + Touched in davidson_miniserver_get + + +`dav_size `_ + Size of the arrays for Davidson + .br + Touched in davidson_miniserver_get + + +`dav_ut `_ + Temporary arrays for parallel davidson + .br + Touched in davidson_miniserver_get + + +`davidson_add_task `_ + Undocumented + + +`davidson_collect `_ + Undocumented + + +`davidson_collector `_ + Undocumented + + +`davidson_converged `_ + True if the Davidson algorithm is converged + + +`davidson_criterion `_ + Can be : [ energy | residual | both | wall_time | cpu_time | iterations ] + + +`davidson_diag `_ + Davidson diagonalization. + .br + dets_in : bitmasks corresponding to determinants + .br + u_in : guess coefficients on the various states. Overwritten + on exit + .br + dim_in : leftmost dimension of u_in + .br + sze : Number of determinants + .br + N_st : Number of eigenstates + .br + iunit : Unit number for the I/O + .br + Initial guess vectors are not necessarily orthonormal + + +`davidson_diag_hjj `_ + Davidson diagonalization with specific diagonal elements of the H matrix + .br + H_jj : specific diagonal H matrix elements to diagonalize de Davidson + .br + dets_in : bitmasks corresponding to determinants + .br + u_in : guess coefficients on the various states. Overwritten + on exit + .br + dim_in : leftmost dimension of u_in + .br + sze : Number of determinants + .br + N_st : Number of eigenstates + .br + N_st_diag : Number of states in which H is diagonalized + .br + iunit : Unit for the I/O + .br + Initial guess vectors are not necessarily orthonormal + + +`davidson_diag_hjj_sjj `_ + Davidson diagonalization with specific diagonal elements of the H matrix + .br + H_jj : specific diagonal H matrix elements to diagonalize de Davidson + .br + S2_jj : specific diagonal S^2 matrix elements + .br + dets_in : bitmasks corresponding to determinants + .br + u_in : guess coefficients on the various states. Overwritten + on exit + .br + dim_in : leftmost dimension of u_in + .br + sze : Number of determinants + .br + N_st : Number of eigenstates + .br + N_st_diag : Number of states in which H is diagonalized. Assumed > sze + .br + iunit : Unit for the I/O + .br + Initial guess vectors are not necessarily orthonormal + + +`davidson_diag_hs2 `_ + Davidson diagonalization. + .br + dets_in : bitmasks corresponding to determinants + .br + u_in : guess coefficients on the various states. Overwritten + on exit + .br + dim_in : leftmost dimension of u_in + .br + sze : Number of determinants + .br + N_st : Number of eigenstates + .br + iunit : Unit number for the I/O + .br + Initial guess vectors are not necessarily orthonormal + + +`davidson_init `_ + Undocumented + + +`davidson_iter_max `_ + Max number of Davidson iterations + + +`davidson_miniserver_end `_ + Undocumented + + +`davidson_miniserver_get `_ + Undocumented + + +`davidson_miniserver_run `_ + Undocumented + + +`davidson_process `_ + Undocumented + + +`davidson_pull_results `_ + Undocumented + + +`davidson_push_results `_ + Undocumented + + +`davidson_run `_ + Undocumented + + +`davidson_run_slave `_ + Undocumented + + +`davidson_slave `_ + Undocumented + + +`davidson_slave_inproc `_ + Undocumented + + +`davidson_slave_tcp `_ + Undocumented + + +`davidson_slave_work `_ + Undocumented + + +`davidson_sze_max `_ + Max number of Davidson sizes + + +`det_inf `_ + Ordering function for determinants + + +`diagonalize_ci `_ + Replace the coefficients of the CI states by the coefficients of the + eigenstates of the CI matrix + + +`diagonalize_ci_mono `_ + Replace the coefficients of the CI states by the coefficients of the + eigenstates of the CI matrix + + +`first_guess `_ + Select all the determinants with the lowest energy as a starting point. + + +`h_s2_u_0_nstates `_ + Computes v_0 = H|u_0> and s_0 = S^2 |u_0> + .br + n : number of determinants + .br + H_jj : array of + .br + S2_jj : array of + + +`h_u_0_nstates `_ + Computes v_0 = H|u_0> + .br + n : number of determinants + .br + H_jj : array of + + +`max_blocksize `_ + Undocumented + + +`n_states_diag `_ + n_states_diag + + +`provide_everything `_ + Undocumented + + +`psi_energy `_ + Energy of the current wave function + + +`shortcut_ `_ + Undocumented + + +`sort_dets_ab `_ + Uncodumented : TODO + + +`sort_dets_ab_v `_ + Uncodumented : TODO + + +`sort_dets_ba_v `_ + Uncodumented : TODO + + +`sort_idx_ `_ + Undocumented + + +`sorted_ `_ + Undocumented + + +`tamiser `_ + Uncodumented : TODO + + +`threshold_davidson `_ + Thresholds of Davidson's algorithm + + +`u_0_h_u_0 `_ + Computes e_0 = / + .br + n : number of determinants + .br + + +`version_ `_ + Undocumented + diff --git a/src/Davidson/find_reference.irp.f b/src/Davidson/find_reference.irp.f new file mode 100644 index 00000000..0cafd739 --- /dev/null +++ b/src/Davidson/find_reference.irp.f @@ -0,0 +1,41 @@ +subroutine find_reference(thresh,n_ref,result) + implicit none + double precision, intent(in) :: thresh + integer, intent(out) :: result(N_det),n_ref + integer :: i,j,istate + double precision :: i_H_psi_array(1), E0, hii, norm + double precision :: de + integer(bit_kind), allocatable :: psi_ref_(:,:,:) + double precision, allocatable :: psi_ref_coef_(:,:) + + allocate(psi_ref_coef_(N_det,1), psi_ref_(N_int,2,N_det)) + n_ref = 1 + result(1) = 1 + istate = 1 + psi_ref_coef_(1,1) = psi_coef(1,istate) + psi_ref_(:,:,1) = psi_det(:,:,1) + norm = psi_ref_coef_(1,1) * psi_ref_coef_(1,1) + call u_0_H_u_0(E0,psi_ref_coef_,n_ref,psi_ref_,N_int,1,size(psi_ref_coef_,1)) + print *, '' + print *, 'Reference determinants' + print *, '======================' + print *, '' + print *, n_ref, ': E0 = ', E0 + nuclear_repulsion + call debug_det(psi_ref_(1,1,n_ref),N_int) + do i=2,N_det + call i_h_psi(psi_det(1,1,i),psi_ref_(1,1,1),psi_ref_coef_(1,istate),N_int, & + n_ref,size(psi_ref_coef_,1),1,i_H_psi_array) + call i_H_j(psi_det(1,1,i),psi_det(1,1,i),N_int,hii) + de = i_H_psi_array(istate)**2 / (E0 - hii) + if (dabs(de) > thresh) then + n_ref += 1 + result(n_ref) = i + psi_ref_(:,:,n_ref) = psi_det(:,:,i) + psi_ref_coef_(n_ref,1) = psi_coef(i,istate) + call u_0_H_u_0(E0,psi_ref_coef_,n_ref,psi_ref_,N_int,1,size(psi_ref_coef_,1)) + print *, n_ref, ': E0 = ', E0 + nuclear_repulsion + call debug_det(psi_ref_(1,1,n_ref),N_int) + endif + enddo +end + diff --git a/src/Davidson/tree_dependency.png b/src/Davidson/tree_dependency.png new file mode 100644 index 00000000..e69de29b diff --git a/src/Davidson/u0Hu0.irp.f b/src/Davidson/u0Hu0.irp.f index 2589f0b3..9e76bc92 100644 --- a/src/Davidson/u0Hu0.irp.f +++ b/src/Davidson/u0Hu0.irp.f @@ -32,20 +32,20 @@ subroutine H_u_0_nstates(v_0,u_0,H_jj,n,keys_tmp,Nint,N_st,sze_8) use bitmasks implicit none BEGIN_DOC - ! Computes v_0 = H|u_0> + ! Computes v_0 = H|u_0> ! ! n : number of determinants ! ! H_jj : array of + ! END_DOC integer, intent(in) :: N_st,n,Nint, sze_8 double precision, intent(out) :: v_0(sze_8,N_st) double precision, intent(in) :: u_0(sze_8,N_st) double precision, intent(in) :: H_jj(n) integer(bit_kind),intent(in) :: keys_tmp(Nint,2,n) - double precision :: hij - double precision, allocatable :: vt(:,:) - double precision, allocatable :: ut(:,:) + double precision :: hij,s2 + double precision, allocatable :: vt(:,:), ut(:,:), st(:,:) integer :: i,j,k,l, jj,ii integer :: i0, j0 @@ -57,38 +57,79 @@ subroutine H_u_0_nstates(v_0,u_0,H_jj,n,keys_tmp,Nint,N_st,sze_8) integer :: N_st_8 integer, external :: align_double - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: vt, ut + integer :: blockb, blockb2, istep + double precision :: ave_workload, workload, target_workload_inv + + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: vt, ut, st N_st_8 = align_double(N_st) ASSERT (Nint > 0) ASSERT (Nint == N_int) ASSERT (n>0) - PROVIDE ref_bitmask_energy + PROVIDE ref_bitmask_energy allocate (shortcut(0:n+1,2), sort_idx(n,2), sorted(Nint,n,2), version(Nint,n,2)) - allocate(ut(N_st_8,n)) + allocate( ut(N_st_8,n)) v_0 = 0.d0 - do i=1,n - do istate=1,N_st - ut(istate,i) = u_0(i,istate) - enddo - enddo - call sort_dets_ab_v(keys_tmp, sorted(1,1,1), sort_idx(1,1), shortcut(0,1), version(1,1,1), n, Nint) call sort_dets_ba_v(keys_tmp, sorted(1,1,2), sort_idx(1,2), shortcut(0,2), version(1,1,2), n, Nint) - + !$OMP PARALLEL DEFAULT(NONE) & - !$OMP PRIVATE(i,hij,j,k,jj,vt,ii,sh,sh2,ni,exa,ext,org_i,org_j,endi,sorted_i,istate)& - !$OMP SHARED(n,H_jj,keys_tmp,ut,Nint,v_0,sorted,shortcut,sort_idx,version,N_st,N_st_8) - allocate(vt(N_st_8,n)) + !$OMP PRIVATE(i,hij,s2,j,k,jj,vt,st,ii,sh,sh2,ni,exa,ext,org_i,org_j,endi,sorted_i,istate)& + !$OMP SHARED(n,keys_tmp,ut,Nint,u_0,v_0,sorted,shortcut,sort_idx,version,N_st,N_st_8) + allocate(vt(N_st_8,n),st(N_st_8,n)) Vt = 0.d0 - + St = 0.d0 + + !$OMP DO + do i=1,n + do istate=1,N_st + ut(istate,i) = u_0(sort_idx(i,2),istate) + enddo + enddo + !$OMP END DO + + !$OMP DO SCHEDULE(dynamic) + do sh=1,shortcut(0,2) + do i=shortcut(sh,2),shortcut(sh+1,2)-1 + org_i = sort_idx(i,2) + do j=shortcut(sh,2),shortcut(sh+1,2)-1 + org_j = sort_idx(j,2) + ext = popcnt(xor(sorted(1,i,2), sorted(1,j,2))) + if (ext > 4) cycle + do ni=2,Nint + ext = ext + popcnt(xor(sorted(ni,i,2), sorted(ni,j,2))) + if (ext > 4) exit + end do + if(ext == 4) then + call i_h_j (keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),nint,hij) + call get_s2(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),nint,s2) + do istate=1,n_st + vt (istate,org_i) = vt (istate,org_i) + hij*ut(istate,j) + st (istate,org_i) = st (istate,org_i) + s2*ut(istate,j) + enddo + end if + end do + end do + enddo + !$OMP END DO + + !$OMP DO + do i=1,n + do istate=1,N_st + ut(istate,i) = u_0(sort_idx(i,1),istate) + enddo + enddo + !$OMP END DO + !$OMP DO SCHEDULE(dynamic) do sh=1,shortcut(0,1) - do sh2=sh,shortcut(0,1) + do sh2=1,shortcut(0,1) + if (sh==sh2) cycle + exa = 0 do ni=1,Nint exa = exa + popcnt(xor(version(ni,sh,1), version(ni,sh2,1))) @@ -96,78 +137,126 @@ subroutine H_u_0_nstates(v_0,u_0,H_jj,n,keys_tmp,Nint,N_st,sze_8) if(exa > 2) then cycle end if - + do i=shortcut(sh,1),shortcut(sh+1,1)-1 org_i = sort_idx(i,1) - if(sh==sh2) then - endi = i-1 - else - endi = shortcut(sh2+1,1)-1 - end if do ni=1,Nint sorted_i(ni) = sorted(ni,i,1) enddo - - do j=shortcut(sh2,1),endi - org_j = sort_idx(j,1) - ext = exa - do ni=1,Nint + + do j=shortcut(sh2,1),shortcut(sh2+1,1)-1 + ext = exa + popcnt(xor(sorted_i(1), sorted(1,j,1))) + if (ext > 4) cycle + do ni=2,Nint ext = ext + popcnt(xor(sorted_i(ni), sorted(ni,j,1))) + if (ext > 4) exit end do if(ext <= 4) then - call i_H_j(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),Nint,hij) - do istate=1,N_st - vt (istate,org_i) = vt (istate,org_i) + hij*ut(istate,org_j) - vt (istate,org_j) = vt (istate,org_j) + hij*ut(istate,org_i) - enddo + org_j = sort_idx(j,1) + call i_h_j (keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),nint,hij) + if (hij /= 0.d0) then + do istate=1,n_st + vt (istate,org_i) = vt (istate,org_i) + hij*ut(istate,j) + enddo + endif + if (ext /= 2) then + call get_s2(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),nint,s2) + if (s2 /= 0.d0) then + do istate=1,n_st + st (istate,org_i) = st (istate,org_i) + s2*ut(istate,j) + enddo + endif + endif endif enddo + + enddo + enddo + + exa = 0 + + do i=shortcut(sh,1),shortcut(sh+1,1)-1 + org_i = sort_idx(i,1) + do ni=1,Nint + sorted_i(ni) = sorted(ni,i,1) + enddo + + do j=shortcut(sh,1),i-1 + ext = exa + popcnt(xor(sorted_i(1), sorted(1,j,1))) + if (ext > 4) cycle + do ni=2,Nint + ext = ext + popcnt(xor(sorted_i(ni), sorted(ni,j,1))) + if (ext > 4) exit + end do + if(ext <= 4) then + org_j = sort_idx(j,1) + call i_h_j (keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),nint,hij) + if (hij /= 0.d0) then + do istate=1,n_st + vt (istate,org_i) = vt (istate,org_i) + hij*ut(istate,j) + enddo + endif + if (ext /= 2) then + call get_s2(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),nint,s2) + if (s2 /= 0.d0) then + do istate=1,n_st + st (istate,org_i) = st (istate,org_i) + s2*ut(istate,j) + enddo + endif + endif + endif + enddo + + do j=i+1,shortcut(sh+1,1)-1 + if (i==j) cycle + ext = exa + popcnt(xor(sorted_i(1), sorted(1,j,1))) + if (ext > 4) cycle + do ni=2,Nint + ext = ext + popcnt(xor(sorted_i(ni), sorted(ni,j,1))) + if (ext > 4) exit + end do + if(ext <= 4) then + org_j = sort_idx(j,1) + call i_h_j (keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),nint,hij) + if (hij /= 0.d0) then + do istate=1,n_st + vt (istate,org_i) = vt (istate,org_i) + hij*ut(istate,j) + enddo + endif + if (ext /= 2) then + call get_s2(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),nint,s2) + if (s2 /= 0.d0) then + do istate=1,n_st + st (istate,org_i) = st (istate,org_i) + s2*ut(istate,j) + enddo + endif + endif + endif enddo enddo enddo - !$OMP END DO NOWAIT - - !$OMP DO SCHEDULE(dynamic) - do sh=1,shortcut(0,2) - do i=shortcut(sh,2),shortcut(sh+1,2)-1 - org_i = sort_idx(i,2) - do j=shortcut(sh,2),i-1 - org_j = sort_idx(j,2) - ext = 0 - do ni=1,Nint - ext = ext + popcnt(xor(sorted(ni,i,2), sorted(ni,j,2))) - end do - if(ext == 4) then - call i_H_j(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),Nint,hij) - do istate=1,N_st - vt (istate,org_i) = vt (istate,org_i) + hij*ut(istate,org_j) - vt (istate,org_j) = vt (istate,org_j) + hij*ut(istate,org_i) - enddo - end if - end do - end do - enddo - !$OMP END DO NOWAIT - - !$OMP CRITICAL + !$OMP END DO + + !$OMP CRITICAL (u0Hu0) do istate=1,N_st - do i=n,1,-1 + do i=1,n v_0(i,istate) = v_0(i,istate) + vt(istate,i) enddo enddo - !$OMP END CRITICAL + !$OMP END CRITICAL (u0Hu0) - deallocate(vt) + deallocate(vt,st) !$OMP END PARALLEL - + do istate=1,N_st do i=1,n - v_0(i,istate) += H_jj(i) * u_0(i,istate) + v_0(i,istate) = v_0(i,istate) + H_jj(i) * u_0(i,istate) enddo enddo deallocate (shortcut, sort_idx, sorted, version, ut) end + BEGIN_PROVIDER [ double precision, psi_energy, (N_states) ] implicit none BEGIN_DOC @@ -323,35 +412,38 @@ subroutine H_S2_u_0_nstates(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,N_st,sze_8) PROVIDE ref_bitmask_energy allocate (shortcut(0:n+1,2), sort_idx(n,2), sorted(Nint,n,2), version(Nint,n,2)) - allocate(ut(N_st_8,n)) + allocate( ut(N_st_8,n)) v_0 = 0.d0 s_0 = 0.d0 - do i=1,n - do istate=1,N_st - ut(istate,i) = u_0(i,istate) - enddo - enddo - call sort_dets_ab_v(keys_tmp, sorted(1,1,1), sort_idx(1,1), shortcut(0,1), version(1,1,1), n, Nint) call sort_dets_ba_v(keys_tmp, sorted(1,1,2), sort_idx(1,2), shortcut(0,2), version(1,1,2), n, Nint) !$OMP PARALLEL DEFAULT(NONE) & !$OMP PRIVATE(i,hij,s2,j,k,jj,vt,st,ii,sh,sh2,ni,exa,ext,org_i,org_j,endi,sorted_i,istate)& - !$OMP SHARED(n,keys_tmp,ut,Nint,v_0,s_0,sorted,shortcut,sort_idx,version,N_st,N_st_8) + !$OMP SHARED(n,keys_tmp,ut,Nint,u_0,v_0,s_0,sorted,shortcut,sort_idx,version,N_st,N_st_8) allocate(vt(N_st_8,n),st(N_st_8,n)) Vt = 0.d0 St = 0.d0 + !$OMP DO + do i=1,n + do istate=1,N_st + ut(istate,i) = u_0(sort_idx(i,2),istate) + enddo + enddo + !$OMP END DO + !$OMP DO SCHEDULE(dynamic) do sh=1,shortcut(0,2) do i=shortcut(sh,2),shortcut(sh+1,2)-1 org_i = sort_idx(i,2) - do j=shortcut(sh,2),i-1 + do j=shortcut(sh,2),shortcut(sh+1,2)-1 org_j = sort_idx(j,2) - ext = 0 - do ni=1,Nint + ext = popcnt(xor(sorted(1,i,2), sorted(1,j,2))) + if (ext > 4) cycle + do ni=2,Nint ext = ext + popcnt(xor(sorted(ni,i,2), sorted(ni,j,2))) if (ext > 4) exit end do @@ -359,19 +451,28 @@ subroutine H_S2_u_0_nstates(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,N_st,sze_8) call i_h_j (keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),nint,hij) call get_s2(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),nint,s2) do istate=1,n_st - vt (istate,org_i) = vt (istate,org_i) + hij*ut(istate,org_j) - vt (istate,org_j) = vt (istate,org_j) + hij*ut(istate,org_i) - st (istate,org_i) = st (istate,org_i) + s2*ut(istate,org_j) - st (istate,org_j) = st (istate,org_j) + s2*ut(istate,org_i) + vt (istate,org_i) = vt (istate,org_i) + hij*ut(istate,j) + st (istate,org_i) = st (istate,org_i) + s2*ut(istate,j) enddo end if end do end do enddo - !$OMP END DO NOWAIT + !$OMP END DO + + !$OMP DO + do i=1,n + do istate=1,N_st + ut(istate,i) = u_0(sort_idx(i,1),istate) + enddo + enddo + !$OMP END DO + + !$OMP DO SCHEDULE(dynamic) do sh=1,shortcut(0,1) - !$OMP DO SCHEDULE(dynamic) - do sh2=sh,shortcut(0,1) + do sh2=1,shortcut(0,1) + if (sh==sh2) cycle + exa = 0 do ni=1,Nint exa = exa + popcnt(xor(version(ni,sh,1), version(ni,sh2,1))) @@ -382,18 +483,14 @@ subroutine H_S2_u_0_nstates(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,N_st,sze_8) do i=shortcut(sh,1),shortcut(sh+1,1)-1 org_i = sort_idx(i,1) - if(sh==sh2) then - endi = i-1 - else - endi = shortcut(sh2+1,1)-1 - end if do ni=1,Nint sorted_i(ni) = sorted(ni,i,1) enddo - do j=shortcut(sh2,1),endi - ext = exa - do ni=1,Nint + do j=shortcut(sh2,1),shortcut(sh2+1,1)-1 + ext = exa + popcnt(xor(sorted_i(1), sorted(1,j,1))) + if (ext > 4) cycle + do ni=2,Nint ext = ext + popcnt(xor(sorted_i(ni), sorted(ni,j,1))) if (ext > 4) exit end do @@ -402,25 +499,86 @@ subroutine H_S2_u_0_nstates(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,N_st,sze_8) call i_h_j (keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),nint,hij) if (hij /= 0.d0) then do istate=1,n_st - vt (istate,org_i) = vt (istate,org_i) + hij*ut(istate,org_j) - vt (istate,org_j) = vt (istate,org_j) + hij*ut(istate,org_i) + vt (istate,org_i) = vt (istate,org_i) + hij*ut(istate,j) enddo endif if (ext /= 2) then call get_s2(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),nint,s2) if (s2 /= 0.d0) then do istate=1,n_st - st (istate,org_i) = st (istate,org_i) + s2*ut(istate,org_j) - st (istate,org_j) = st (istate,org_j) + s2*ut(istate,org_i) + st (istate,org_i) = st (istate,org_i) + s2*ut(istate,j) enddo endif endif endif enddo + + enddo + enddo + + exa = 0 + + do i=shortcut(sh,1),shortcut(sh+1,1)-1 + org_i = sort_idx(i,1) + do ni=1,Nint + sorted_i(ni) = sorted(ni,i,1) + enddo + + do j=shortcut(sh,1),i-1 + ext = exa + popcnt(xor(sorted_i(1), sorted(1,j,1))) + if (ext > 4) cycle + do ni=2,Nint + ext = ext + popcnt(xor(sorted_i(ni), sorted(ni,j,1))) + if (ext > 4) exit + end do + if(ext <= 4) then + org_j = sort_idx(j,1) + call i_h_j (keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),nint,hij) + if (hij /= 0.d0) then + do istate=1,n_st + vt (istate,org_i) = vt (istate,org_i) + hij*ut(istate,j) + enddo + endif + if (ext /= 2) then + call get_s2(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),nint,s2) + if (s2 /= 0.d0) then + do istate=1,n_st + st (istate,org_i) = st (istate,org_i) + s2*ut(istate,j) + enddo + endif + endif + endif + enddo + + do j=i+1,shortcut(sh+1,1)-1 + if (i==j) cycle + ext = exa + popcnt(xor(sorted_i(1), sorted(1,j,1))) + if (ext > 4) cycle + do ni=2,Nint + ext = ext + popcnt(xor(sorted_i(ni), sorted(ni,j,1))) + if (ext > 4) exit + end do + if(ext <= 4) then + org_j = sort_idx(j,1) + call i_h_j (keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),nint,hij) + if (hij /= 0.d0) then + do istate=1,n_st + vt (istate,org_i) = vt (istate,org_i) + hij*ut(istate,j) + enddo + endif + if (ext /= 2) then + call get_s2(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),nint,s2) + if (s2 /= 0.d0) then + do istate=1,n_st + st (istate,org_i) = st (istate,org_i) + s2*ut(istate,j) + enddo + endif + endif + endif enddo enddo - !$OMP END DO NOWAIT enddo + !$OMP END DO !$OMP CRITICAL (u0Hu0) do istate=1,N_st diff --git a/src/Determinants/Fock_diag.irp.f b/src/Determinants/Fock_diag.irp.f index a99bbcad..01393fe1 100644 --- a/src/Determinants/Fock_diag.irp.f +++ b/src/Determinants/Fock_diag.irp.f @@ -19,6 +19,15 @@ subroutine build_fock_tmp(fock_diag_tmp,det_ref,Nint) fock_diag_tmp = 0.d0 E0 = 0.d0 + if (Ne(1) /= elec_alpha_num) then + print *, 'Error in build_fock_tmp (alpha)', Ne(1), Ne(2) + stop -1 + endif + if (Ne(2) /= elec_beta_num) then + print *, 'Error in build_fock_tmp (beta)', Ne(1), Ne(2) + stop -1 + endif + ! Occupied MOs do ii=1,elec_alpha_num i = occ(ii,1) diff --git a/src/Determinants/README.rst b/src/Determinants/README.rst index c6685945..9ad0f1a3 100644 --- a/src/Determinants/README.rst +++ b/src/Determinants/README.rst @@ -15,23 +15,31 @@ Documentation .. by the `update_README.py` script. -`a_operator `_ +`a_operator `_ Needed for diag_H_mat_elem -`abs_psi_coef_max `_ +`abs_psi_coef_max `_ Max and min values of the coefficients -`abs_psi_coef_min `_ +`abs_psi_coef_min `_ Max and min values of the coefficients -`ac_operator `_ +`ac_operator `_ Needed for diag_H_mat_elem -`apply_excitation `_ +`apply_excitation `_ + Undocumented + + +`apply_hole `_ + Undocumented + + +`apply_holes `_ Undocumented @@ -39,16 +47,24 @@ Documentation Undocumented +`apply_particle `_ + Undocumented + + +`apply_particles `_ + Undocumented + + `bi_elec_ref_bitmask_energy `_ Energy of the reference bitmask used in Slater rules -`bitstring_to_list_ab `_ +`bitstring_to_list_ab `_ Gives the inidices(+1) of the bits set to 1 in the bit string For alpha/beta determinants -`bitstring_to_list_ab_old `_ +`bitstring_to_list_ab_old `_ Gives the inidices(+1) of the bits set to 1 in the bit string For alpha/beta determinants @@ -58,72 +74,15 @@ Documentation determinant. F_00 is = E0. -`ci_eigenvectors `_ - Eigenvectors/values of the CI matrix - - -`ci_eigenvectors_mono `_ - Eigenvectors/values of the CI matrix - - -`ci_eigenvectors_s2 `_ - Eigenvectors/values of the CI matrix - - -`ci_eigenvectors_s2_mono `_ - Eigenvectors/values of the CI matrix - - -`ci_electronic_energy `_ - Eigenvectors/values of the CI matrix - - -`ci_electronic_energy_mono `_ - Eigenvectors/values of the CI matrix - - -`ci_energy `_ - N_states lowest eigenvalues of the CI matrix - - -`ci_sc2_eigenvectors `_ - Eigenvectors/values of the CI matrix - - -`ci_sc2_electronic_energy `_ - Eigenvectors/values of the CI matrix - - -`ci_sc2_energy `_ - N_states_diag lowest eigenvalues of the CI matrix - - `cisd `_ Undocumented -`cisd_sc2 `_ - CISD+SC2 method :: take off all the disconnected terms of a CISD (selected or not) - .br - dets_in : bitmasks corresponding to determinants - .br - u_in : guess coefficients on the various states. Overwritten - on exit - .br - dim_in : leftmost dimension of u_in - .br - sze : Number of determinants - .br - N_st : Number of eigenstates - .br - Initial guess vectors are not necessarily orthonormal - - -`connected_to_ref `_ +`connected_to_ref `_ Undocumented -`connected_to_ref_by_mono `_ +`connected_to_ref_by_mono `_ Undocumented @@ -136,11 +95,11 @@ Documentation Undocumented -`create_minilist `_ +`create_minilist `_ Undocumented -`create_minilist_find_previous `_ +`create_minilist_find_previous `_ Undocumented @@ -149,62 +108,6 @@ Documentation of alpha and beta determinants -`davidson_converged `_ - True if the Davidson algorithm is converged - - -`davidson_criterion `_ - Can be : [ energy | residual | both | wall_time | cpu_time | iterations ] - - -`davidson_diag `_ - Davidson diagonalization. - .br - dets_in : bitmasks corresponding to determinants - .br - u_in : guess coefficients on the various states. Overwritten - on exit - .br - dim_in : leftmost dimension of u_in - .br - sze : Number of determinants - .br - N_st : Number of eigenstates - .br - iunit : Unit number for the I/O - .br - Initial guess vectors are not necessarily orthonormal - - -`davidson_diag_hjj `_ - Davidson diagonalization with specific diagonal elements of the H matrix - .br - H_jj : specific diagonal H matrix elements to diagonalize de Davidson - .br - dets_in : bitmasks corresponding to determinants - .br - u_in : guess coefficients on the various states. Overwritten - on exit - .br - dim_in : leftmost dimension of u_in - .br - sze : Number of determinants - .br - N_st : Number of eigenstates - .br - iunit : Unit for the I/O - .br - Initial guess vectors are not necessarily orthonormal - - -`davidson_iter_max `_ - Max number of Davidson iterations - - -`davidson_sze_max `_ - Max number of Davidson sizes - - `decode_exc `_ Decodes the exc arrays returned by get_excitation. h1,h2 : Holes @@ -213,6 +116,14 @@ Documentation degree : Degree of excitation +`decode_exc_int2 `_ + Decodes the exc arrays returned by get_excitation. + h1,h2 : Holes + p1,p2 : Particles + s1,s2 : Spins (1:alpha, 2:beta) + degree : Degree of excitation + + `det_alpha_norm `_ Norm of the alpha and beta spin determinants in the wave function: .br @@ -225,15 +136,11 @@ Documentation ||Da||_i \sum_j C_{ij}**2 -`det_coef `_ +`det_coef `_ det_coef -`det_inf `_ - Undocumented - - -`det_occ `_ +`det_occ `_ det_occ @@ -245,44 +152,29 @@ Documentation Transform a determinant to an occupation pattern -`diag_algorithm `_ +`detcmp `_ + Undocumented + + +`deteq `_ + Undocumented + + +`diag_algorithm `_ Diagonalization algorithm (Davidson or Lapack) -`diag_h_elements_sc2 `_ - Eigenvectors/values of the CI matrix - - -`diag_h_mat_elem `_ +`diag_h_mat_elem `_ Computes -`diag_h_mat_elem_fock `_ +`diag_h_mat_elem_fock `_ Computes when i is at most a double excitation from a reference. -`diagonalize_ci `_ - Replace the coefficients of the CI states by the coefficients of the - eigenstates of the CI matrix - - -`diagonalize_ci_mono `_ - Replace the coefficients of the CI states by the coefficients of the - eigenstates of the CI matrix - - -`diagonalize_ci_sc2 `_ - Replace the coefficients of the CI states_diag by the coefficients of the - eigenstates of the CI matrix - - -`diagonalize_s2 `_ - Diagonalize the S^2 operator within the n_states_diag states required. Notice : the vectors are sorted by increasing S^2 values. - - -`diagonalize_s2_betweenstates `_ - You enter with nstates vectors in psi_coefs_inout that may be coupled by S^2 +`diagonalize_s2_betweenstates `_ + You enter with nstates vectors in u_0 that may be coupled by S^2 The subroutine diagonalize the S^2 operator in the basis of these states. The vectors that you obtain in output are no more coupled by S^2, which does not necessary mean that they are eigenfunction of S^2. @@ -349,7 +241,7 @@ Documentation idx(0) is the number of determinants that interact with key1 -`filter_connected_i_h_psi0 `_ +`filter_connected_i_h_psi0 `_ returns the array idx which contains the index of the .br determinants in the array key1 that interact @@ -359,7 +251,7 @@ Documentation idx(0) is the number of determinants that interact with key1 -`filter_connected_i_h_psi0_sc2 `_ +`filter_connected_i_h_psi0_sc2 `_ standard filter_connected_i_H_psi but returns in addition .br the array of the index of the non connected determinants to key1 @@ -371,18 +263,22 @@ Documentation to repeat the excitations -`first_guess `_ - Select all the determinants with the lowest energy as a starting point. +`flip_generators `_ + Undocumented `generate_all_alpha_beta_det_products `_ Create a wave function from all possible alpha x beta determinants -`get_double_excitation `_ +`get_double_excitation `_ Returns the two excitation operators between two doubly excited determinants and the phase +`get_double_excitation_phase `_ + Undocumented + + `get_excitation `_ Returns the excitation operators between two determinants and the phase @@ -391,7 +287,7 @@ Documentation Returns the excitation degree between two determinants -`get_excitation_degree_vector `_ +`get_excitation_degree_vector `_ Applies get_excitation_degree to an array of determinants @@ -407,27 +303,23 @@ Documentation Returns the index of the determinant in the ``psi_det_sorted_bit`` array -`get_mono_excitation `_ +`get_mono_excitation `_ Returns the excitation operator between two singly excited determinants and the phase -`get_occ_from_key `_ +`get_occ_from_key `_ Returns a list of occupation numbers from a bitstring +`get_phase `_ + Returns the phase between key1 and key2 + + `get_s2 `_ Returns -`get_s2_u0 `_ - Undocumented - - -`get_s2_u0_old `_ - Undocumented - - -`get_uj_s2_ui `_ +`get_uj_s2_ui `_ returns the matrix elements of S^2 "s2(i,j)" between the "nstates" states psi_coefs_tmp(:,i) and psi_coefs_tmp(:,j) @@ -458,27 +350,19 @@ Documentation Undocumented -`h_u_0 `_ - Computes v_0 = H|u_0> - .br - n : number of determinants - .br - H_jj : array of - - -`i_h_j `_ +`i_h_j `_ Returns where i and j are determinants -`i_h_j_phase_out `_ +`i_h_j_phase_out `_ Returns where i and j are determinants -`i_h_j_verbose `_ +`i_h_j_verbose `_ Returns where i and j are determinants -`i_h_psi `_ +`i_h_psi `_ Computes = \sum_J c_J . .br Uses filter_connected_i_H_psi0 to get all the |J> to which |i> @@ -487,14 +371,14 @@ Documentation minilists -`i_h_psi_minilist `_ +`i_h_psi_minilist `_ Computes = \sum_J c_J . .br Uses filter_connected_i_H_psi0 to get all the |J> to which |i> is connected. The |J> are searched in short pre-computed lists. -`i_h_psi_sc2 `_ +`i_h_psi_sc2 `_ for the various Nstate .br returns in addition @@ -508,7 +392,7 @@ Documentation to repeat the excitations -`i_h_psi_sc2_verbose `_ +`i_h_psi_sc2_verbose `_ for the various Nstate .br returns in addition @@ -522,10 +406,17 @@ Documentation to repeat the excitations -`i_h_psi_sec_ord `_ +`i_h_psi_sec_ord `_ for the various Nstates +`i_s2_psi_minilist `_ + Computes = \sum_J c_J . + .br + Uses filter_connected_i_H_psi0 to get all the |J> to which |i> + is connected. The |J> are searched in short pre-computed lists. + + `idx_cas `_ CAS wave function, defined from the application of the CAS bitmask on the determinants. idx_cas gives the indice of the CAS determinant in psi_det. @@ -537,11 +428,15 @@ Documentation idx_non_cas gives the indice of the determinant in psi_det. -`is_connected_to `_ +`is_connected_to `_ Undocumented -`is_connected_to_by_mono `_ +`is_connected_to_by_mono `_ + Undocumented + + +`is_generable_cassd `_ Undocumented @@ -557,7 +452,7 @@ Documentation Undocumented -`max_degree_exc `_ +`max_degree_exc `_ Maximum degree of excitation in the wf @@ -573,7 +468,7 @@ Documentation Undocumented -`n_det `_ +`n_det `_ Number of determinants in the wave function @@ -598,7 +493,7 @@ Documentation Maximum number of determinants diagonalized by Jacobi -`n_det_max_property `_ +`n_det_max_property `_ Max number of determinants in the wave function when you select for a given property @@ -630,10 +525,6 @@ Documentation Number of states to consider -`n_states_diag `_ - Number of states to consider for the diagonalization - - `neutral_no_hund_in_couple `_ n_couples is the number of couples of orbitals to be checked couples(i,1) = first orbital of the ith couple @@ -696,15 +587,15 @@ Documentation rho(alpha) - rho(beta) -`only_single_double_dm `_ +`only_single_double_dm `_ If true, The One body DM is calculated with ignoring the Double<->Doubles extra diag elements -`psi_average_norm_contrib `_ +`psi_average_norm_contrib `_ Contribution of determinants to the state-averaged density -`psi_average_norm_contrib_sorted `_ +`psi_average_norm_contrib_sorted `_ Wave function sorted by determinants contribution to the norm (state-averaged) @@ -756,7 +647,7 @@ Documentation function. -`psi_coef `_ +`psi_coef `_ The wave function coefficients. Initialized with Hartree-Fock if the EZFIO file is empty @@ -765,26 +656,26 @@ Documentation Undocumented -`psi_coef_max `_ +`psi_coef_max `_ Max and min values of the coefficients -`psi_coef_min `_ +`psi_coef_min `_ Max and min values of the coefficients -`psi_coef_sorted `_ +`psi_coef_sorted `_ Wave function sorted by determinants contribution to the norm (state-averaged) -`psi_coef_sorted_bit `_ +`psi_coef_sorted_bit `_ Determinants on which we apply for perturbation. They are sorted by determinants interpreted as integers. Useful to accelerate the search of a random determinant in the wave function. -`psi_det `_ +`psi_det `_ The wave function determinants. Initialized with Hartree-Fock if the EZFIO file is empty @@ -805,15 +696,15 @@ Documentation Unique beta determinants -`psi_det_size `_ +`psi_det_size `_ Size of the psi_det/psi_coef arrays -`psi_det_sorted `_ +`psi_det_sorted `_ Wave function sorted by determinants contribution to the norm (state-averaged) -`psi_det_sorted_bit `_ +`psi_det_sorted_bit `_ Determinants on which we apply for perturbation. They are sorted by determinants interpreted as integers. Useful to accelerate the search of a random determinant in the wave @@ -860,7 +751,7 @@ Documentation Undocumented -`read_dets `_ +`read_dets `_ Reads the determinants from the EZFIO file @@ -885,11 +776,25 @@ Documentation be set before calling this function. -`s2_eig `_ +`s2_eig `_ Force the wave function to be an eigenfunction of S^2 -`s2_values `_ +`s2_u_0 `_ + Computes v_0 = S^2|u_0> + .br + n : number of determinants + .br + + +`s2_u_0_nstates `_ + Computes v_0 = S^2|u_0> + .br + n : number of determinants + .br + + +`s2_values `_ array of the averaged values of the S^2 operator on the various states @@ -913,23 +818,23 @@ Documentation Save natural orbitals, obtained by diagonalization of the one-body density matrix in the MO basis -`save_ref_determinant `_ +`save_ref_determinant `_ Undocumented -`save_wavefunction `_ +`save_wavefunction `_ Save the wave function into the EZFIO file -`save_wavefunction_general `_ +`save_wavefunction_general `_ Save the wave function into the EZFIO file -`save_wavefunction_specified `_ +`save_wavefunction_specified `_ Save the wave function into the EZFIO file -`save_wavefunction_unsorted `_ +`save_wavefunction_unsorted `_ Save the wave function into the EZFIO file @@ -947,49 +852,25 @@ Documentation for a given couple of hole/particle excitations i. -`sort_dets_ab `_ - Uncodumented : TODO - - -`sort_dets_ab_v `_ - Uncodumented : TODO - - -`sort_dets_ba_v `_ - Uncodumented : TODO - - -`sort_dets_by_det_search_key `_ +`sort_dets_by_det_search_key `_ Determinants are sorted are sorted according to their det_search_key. Useful to accelerate the search of a random determinant in the wave function. `spin_det_search_key `_ - Return an integer*8 corresponding to a determinant index for searching + Return an integer(8) corresponding to a determinant index for searching `state_average_weight `_ Weights in the state-average calculation of the density matrix -`tamiser `_ - Uncodumented : TODO - - -`target_energy `_ +`target_energy `_ Energy that should be obtained when truncating the wave function (optional) -`threshold_convergence_sc2 `_ - convergence of the correlation energy of SC2 iterations - - -`threshold_davidson `_ - Thresholds of Davidson's algorithm - - -`threshold_generators `_ +`threshold_generators `_ Thresholds on generators (fraction of the norm) @@ -997,8 +878,8 @@ Documentation Thresholds on selectors (fraction of the norm) -`u0_h_u_0 `_ - Computes e_0 = / +`u_0_s2_u_0 `_ + Computes e_0 = / .br n : number of determinants .br diff --git a/src/Determinants/density_matrix.irp.f b/src/Determinants/density_matrix.irp.f index 118bbdf7..ed2f49bd 100644 --- a/src/Determinants/density_matrix.irp.f +++ b/src/Determinants/density_matrix.irp.f @@ -194,7 +194,6 @@ subroutine set_natural_mos double precision, allocatable :: tmp(:,:) label = "Natural" -! call mo_as_eigvectors_of_mo_matrix(one_body_dm_mo,size(one_body_dm_mo,1),mo_tot_num,label,-1) call mo_as_svd_vectors_of_mo_matrix(one_body_dm_mo,size(one_body_dm_mo,1),mo_tot_num,mo_tot_num,label) end diff --git a/src/Determinants/filter_connected.irp.f b/src/Determinants/filter_connected.irp.f index da333b1e..b76540f7 100644 --- a/src/Determinants/filter_connected.irp.f +++ b/src/Determinants/filter_connected.irp.f @@ -1,4 +1,102 @@ +subroutine filter_not_connected(key1,key2,Nint,sze,idx) + use bitmasks + implicit none + BEGIN_DOC + ! Returns the array idx which contains the index of the + ! + ! determinants in the array key1 that DO NOT interact + ! + ! via the H operator with key2. + ! + ! idx(0) is the number of determinants that DO NOT interact with key1 + END_DOC + integer, intent(in) :: Nint, sze + integer(bit_kind), intent(in) :: key1(Nint,2,sze) + integer(bit_kind), intent(in) :: key2(Nint,2) + integer, intent(out) :: idx(0:sze) + + integer :: i,j,l + integer :: degree_x2 + + + ASSERT (Nint > 0) + ASSERT (sze >= 0) + + l=1 + + if (Nint==1) then + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = popcnt( xor( key1(1,1,i), key2(1,1))) & + + popcnt( xor( key1(1,2,i), key2(1,2))) + if (degree_x2 > 4) then + idx(l) = i + l = l+1 + else + cycle + endif + enddo + + else if (Nint==2) then + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + & + popcnt(xor( key1(2,1,i), key2(2,1))) + & + popcnt(xor( key1(1,2,i), key2(1,2))) + & + popcnt(xor( key1(2,2,i), key2(2,2))) + if (degree_x2 > 4) then + idx(l) = i + l = l+1 + else + cycle + endif + enddo + + else if (Nint==3) then + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + & + popcnt(xor( key1(1,2,i), key2(1,2))) + & + popcnt(xor( key1(2,1,i), key2(2,1))) + & + popcnt(xor( key1(2,2,i), key2(2,2))) + & + popcnt(xor( key1(3,1,i), key2(3,1))) + & + popcnt(xor( key1(3,2,i), key2(3,2))) + if (degree_x2 > 4) then + idx(l) = i + l = l+1 + else + cycle + endif + enddo + + else + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = 0 + !DEC$ LOOP COUNT MIN(4) + do j=1,Nint + degree_x2 = degree_x2+ popcnt(xor( key1(j,1,i), key2(j,1))) +& + popcnt(xor( key1(j,2,i), key2(j,2))) + if (degree_x2 > 4) then + idx(l) = i + l = l+1 + endif + enddo + if (degree_x2 <= 5) then + exit + endif + enddo + + endif + idx(0) = l-1 +end + + subroutine filter_connected(key1,key2,Nint,sze,idx) use bitmasks implicit none diff --git a/src/Determinants/slater_rules.irp.f b/src/Determinants/slater_rules.irp.f index 789dc93c..7556e2b9 100644 --- a/src/Determinants/slater_rules.irp.f +++ b/src/Determinants/slater_rules.irp.f @@ -1041,13 +1041,15 @@ subroutine i_H_psi(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array) double precision :: phase integer :: exc(0:2,2,2) double precision :: hij - integer :: idx(0:Ndet) + integer, allocatable :: idx(:) ASSERT (Nint > 0) ASSERT (N_int == Nint) ASSERT (Nstate > 0) ASSERT (Ndet > 0) ASSERT (Ndet_max >= Ndet) + allocate(idx(0:Ndet)) + i_H_psi_array = 0.d0 call filter_connected_i_H_psi0(keys,key,Nint,Ndet,idx) @@ -1089,7 +1091,7 @@ subroutine i_H_psi_minilist(key,keys,idx_key,N_minilist,coef,Nint,Ndet,Ndet_max, double precision :: phase integer :: exc(0:2,2,2) double precision :: hij - integer :: idx(0:Ndet) + integer, allocatable :: idx(:) BEGIN_DOC ! Computes = \sum_J c_J . ! @@ -1102,6 +1104,7 @@ subroutine i_H_psi_minilist(key,keys,idx_key,N_minilist,coef,Nint,Ndet,Ndet_max, ASSERT (Nstate > 0) ASSERT (Ndet > 0) ASSERT (Ndet_max >= Ndet) + allocate(idx(0:Ndet)) i_H_psi_array = 0.d0 call filter_connected_i_H_psi0(keys,key,Nint,N_minilist,idx) @@ -1148,7 +1151,8 @@ subroutine i_H_psi_sec_ord(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array double precision :: phase integer :: exc(0:2,2,2) double precision :: hij - integer :: idx(0:Ndet),n_interact + integer,allocatable :: idx(:) + integer :: n_interact BEGIN_DOC ! for the various Nstates END_DOC @@ -1158,6 +1162,7 @@ subroutine i_H_psi_sec_ord(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array ASSERT (Nstate > 0) ASSERT (Ndet > 0) ASSERT (Ndet_max >= Ndet) + allocate(idx(0:Ndet)) i_H_psi_array = 0.d0 call filter_connected_i_H_psi0(keys,key,Nint,Ndet,idx) n_interact = 0 @@ -1207,7 +1212,7 @@ subroutine i_H_psi_SC2(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array,idx double precision :: phase integer :: exc(0:2,2,2) double precision :: hij - integer :: idx(0:Ndet) + integer,allocatable :: idx(:) ASSERT (Nint > 0) ASSERT (N_int == Nint) @@ -1215,6 +1220,7 @@ subroutine i_H_psi_SC2(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array,idx ASSERT (Ndet > 0) ASSERT (Ndet_max >= Ndet) i_H_psi_array = 0.d0 + allocate(idx(0:Ndet)) call filter_connected_i_H_psi0_SC2(keys,key,Nint,Ndet,idx,idx_repeat) do ii=1,idx(0) i = idx(ii) @@ -1254,7 +1260,7 @@ subroutine i_H_psi_SC2_verbose(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_a double precision :: phase integer :: exc(0:2,2,2) double precision :: hij - integer :: idx(0:Ndet) + integer,allocatable :: idx(:) ASSERT (Nint > 0) ASSERT (N_int == Nint) @@ -1262,6 +1268,7 @@ subroutine i_H_psi_SC2_verbose(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_a ASSERT (Ndet > 0) ASSERT (Ndet_max >= Ndet) i_H_psi_array = 0.d0 + allocate(idx(0:Ndet)) call filter_connected_i_H_psi0_SC2(keys,key,Nint,Ndet,idx,idx_repeat) print*,'--------' do ii=1,idx(0) diff --git a/src/Electrons/README.rst b/src/Electrons/README.rst index d1c342b5..484617bb 100644 --- a/src/Electrons/README.rst +++ b/src/Electrons/README.rst @@ -44,7 +44,7 @@ Documentation .. by the `update_README.py` script. -`elec_alpha_num `_ +`elec_alpha_num `_ Numbers of electrons alpha ("up") diff --git a/src/Ezfio_files/README.rst b/src/Ezfio_files/README.rst index 6b494339..ad87e7f5 100644 --- a/src/Ezfio_files/README.rst +++ b/src/Ezfio_files/README.rst @@ -219,6 +219,10 @@ output_cas_sd Initial CPU and wall times when printing in the output files +output_davidson + Output file for Davidson + + output_determinants Output file for Determinants @@ -235,6 +239,10 @@ output_full_ci Output file for Full_CI +output_full_ci_zmq + Output file for Full_CI_ZMQ + + output_generators_cas Output file for Generators_CAS @@ -267,14 +275,14 @@ output_moguess Output file for MOGuess -output_mrcc_cassd - Output file for MRCC_CASSD - - output_mrcc_utils Output file for MRCC_Utils +output_mrcepa0 + Output file for mrcepa0 + + output_nuclei Output file for Nuclei diff --git a/src/Integrals_Bielec/README.rst b/src/Integrals_Bielec/README.rst index 98fbbb92..f6644db4 100644 --- a/src/Integrals_Bielec/README.rst +++ b/src/Integrals_Bielec/README.rst @@ -45,7 +45,7 @@ Documentation .. by the `update_README.py` script. -`add_integrals_to_map `_ +`add_integrals_to_map `_ Adds integrals to tha MO map according to some bitmask @@ -54,7 +54,7 @@ Documentation i(r1) j(r1) 1/r12 k(r2) l(r2) -`ao_bielec_integral_schwartz `_ +`ao_bielec_integral_schwartz `_ Needed to compute Schwartz inequalities @@ -68,7 +68,7 @@ Documentation i(r1) j(r2) 1/r12 k(r1) l(r2) -`ao_bielec_integrals_in_map_collector `_ +`ao_bielec_integrals_in_map_collector `_ Collects results from the AO integral calculation @@ -84,11 +84,23 @@ Documentation Computes a buffer of integrals. i is the ID of the current thread. +`ao_integrals_cache `_ + Cache of AO integrals for fast access + + +`ao_integrals_cache_max `_ + Min and max values of the AOs for which the integrals are in the cache + + +`ao_integrals_cache_min `_ + Min and max values of the AOs for which the integrals are in the cache + + `ao_integrals_map `_ AO integrals -`ao_integrals_threshold `_ +`ao_integrals_threshold `_ If || < ao_integrals_threshold then is zero @@ -108,11 +120,11 @@ Documentation Undocumented -`clear_ao_map `_ +`clear_ao_map `_ Frees the memory of the AO map -`clear_mo_map `_ +`clear_mo_map `_ Frees the memory of the MO map @@ -120,15 +132,15 @@ Documentation Compute AO 1/r12 integrals for all i and fixed j,k,l -`compute_ao_integrals_jl `_ +`compute_ao_integrals_jl `_ Parallel client for AO integrals -`disk_access_ao_integrals `_ +`disk_access_ao_integrals `_ Read/Write AO integrals from/to disk [ Write | Read | None ] -`disk_access_mo_integrals `_ +`disk_access_mo_integrals `_ Read/Write MO integrals from/to disk [ Write | Read | None ] @@ -136,15 +148,15 @@ Documentation Compute integrals on the fly -`dump_ao_integrals `_ +`dump_ao_integrals `_ Save to disk the $ao integrals -`dump_mo_integrals `_ +`dump_mo_integrals `_ Save to disk the $ao integrals -`eri `_ +`eri `_ ATOMIC PRIMTIVE bielectronic integral between the 4 primitives :: primitive_1 = x1**(a_x) y1**(a_y) z1**(a_z) exp(-alpha * r1**2) primitive_2 = x1**(b_x) y1**(b_y) z1**(b_z) exp(- beta * r1**2) @@ -166,148 +178,156 @@ Documentation t_w(i,2,k) = t(i) -`general_primitive_integral `_ +`general_primitive_integral `_ Computes the integral where p,q,r,s are Gaussian primitives -`get_ao_bielec_integral `_ +`get_ao_bielec_integral `_ Gets one AO bi-electronic integral from the AO map -`get_ao_bielec_integrals `_ +`get_ao_bielec_integrals `_ Gets multiple AO bi-electronic integral from the AO map . All i are retrieved for j,k,l fixed. -`get_ao_bielec_integrals_non_zero `_ +`get_ao_bielec_integrals_non_zero `_ Gets multiple AO bi-electronic integral from the AO map . All non-zero i are retrieved for j,k,l fixed. -`get_ao_map_size `_ +`get_ao_map_size `_ Returns the number of elements in the AO map -`get_mo_bielec_integral `_ +`get_mo_bielec_integral `_ Returns one integral in the MO basis -`get_mo_bielec_integral_schwartz `_ +`get_mo_bielec_integral_schwartz `_ Returns one integral in the MO basis -`get_mo_bielec_integrals `_ +`get_mo_bielec_integrals `_ Returns multiple integrals in the MO basis, all i for j,k,l fixed. -`get_mo_bielec_integrals_ij `_ +`get_mo_bielec_integrals_ij `_ Returns multiple integrals in the MO basis, all i(1)j(2) 1/r12 k(1)l(2) i, j for k,l fixed. -`get_mo_map_size `_ +`get_mo_map_size `_ Return the number of elements in the MO map -`give_polynom_mult_center_x `_ +`give_polynom_mult_center_x `_ subroutine that returns the explicit polynom in term of the "t" variable of the following polynomw : I_x1(a_x, d_x,p,q) * I_x1(a_y, d_y,p,q) * I_x1(a_z, d_z,p,q) -`i_x1_new `_ +`i_x1_new `_ recursive function involved in the bielectronic integral -`i_x1_pol_mult `_ +`i_x1_pol_mult `_ recursive function involved in the bielectronic integral -`i_x1_pol_mult_a1 `_ +`i_x1_pol_mult_a1 `_ recursive function involved in the bielectronic integral -`i_x1_pol_mult_a2 `_ +`i_x1_pol_mult_a2 `_ recursive function involved in the bielectronic integral -`i_x1_pol_mult_recurs `_ +`i_x1_pol_mult_recurs `_ recursive function involved in the bielectronic integral -`i_x2_new `_ +`i_x2_new `_ recursive function involved in the bielectronic integral -`i_x2_pol_mult `_ +`i_x2_pol_mult `_ recursive function involved in the bielectronic integral -`insert_into_ao_integrals_map `_ +`insert_into_ao_integrals_map `_ Create new entry into AO map -`insert_into_mo_integrals_map `_ +`insert_into_mo_integrals_map `_ Create new entry into MO map, or accumulate in an existing entry -`integrale_new `_ +`integrale_new `_ calculate the integral of the polynom :: I_x1(a_x+b_x, c_x+d_x,p,q) * I_x1(a_y+b_y, c_y+d_y,p,q) * I_x1(a_z+b_z, c_z+d_z,p,q) between ( 0 ; 1) -`load_ao_integrals `_ +`load_ao_integrals `_ Read from disk the $ao integrals -`load_mo_integrals `_ +`load_mo_integrals `_ Read from disk the $ao integrals -`mo_bielec_integral `_ +`mo_bielec_integral `_ Returns one integral in the MO basis -`mo_bielec_integral_jj `_ +`mo_bielec_integral_jj `_ mo_bielec_integral_jj(i,j) = J_ij mo_bielec_integral_jj_exchange(i,j) = K_ij mo_bielec_integral_jj_anti(i,j) = J_ij - K_ij -`mo_bielec_integral_jj_anti `_ +`mo_bielec_integral_jj_anti `_ mo_bielec_integral_jj(i,j) = J_ij mo_bielec_integral_jj_exchange(i,j) = K_ij mo_bielec_integral_jj_anti(i,j) = J_ij - K_ij -`mo_bielec_integral_jj_anti_from_ao `_ +`mo_bielec_integral_jj_anti_from_ao `_ mo_bielec_integral_jj_from_ao(i,j) = J_ij mo_bielec_integral_jj_exchange_from_ao(i,j) = J_ij mo_bielec_integral_jj_anti_from_ao(i,j) = J_ij - K_ij -`mo_bielec_integral_jj_exchange `_ +`mo_bielec_integral_jj_exchange `_ mo_bielec_integral_jj(i,j) = J_ij mo_bielec_integral_jj_exchange(i,j) = K_ij mo_bielec_integral_jj_anti(i,j) = J_ij - K_ij -`mo_bielec_integral_jj_exchange_from_ao `_ +`mo_bielec_integral_jj_exchange_from_ao `_ mo_bielec_integral_jj_from_ao(i,j) = J_ij mo_bielec_integral_jj_exchange_from_ao(i,j) = J_ij mo_bielec_integral_jj_anti_from_ao(i,j) = J_ij - K_ij -`mo_bielec_integral_jj_from_ao `_ +`mo_bielec_integral_jj_from_ao `_ mo_bielec_integral_jj_from_ao(i,j) = J_ij mo_bielec_integral_jj_exchange_from_ao(i,j) = J_ij mo_bielec_integral_jj_anti_from_ao(i,j) = J_ij - K_ij -`mo_bielec_integral_schwartz `_ +`mo_bielec_integral_mipi `_ + and - . Indices are (i,m,p) + + +`mo_bielec_integral_mipi_anti `_ + and - . Indices are (i,m,p) + + +`mo_bielec_integral_schwartz `_ Needed to compute Schwartz inequalities @@ -319,11 +339,23 @@ Documentation Computes an unique index for i,j,k,l integrals -`mo_integrals_map `_ +`mo_integrals_cache `_ + Cache of MO integrals for fast access + + +`mo_integrals_cache_max `_ + Min and max values of the MOs for which the integrals are in the cache + + +`mo_integrals_cache_min `_ + Min and max values of the MOs for which the integrals are in the cache + + +`mo_integrals_map `_ MO integrals -`mo_integrals_threshold `_ +`mo_integrals_threshold `_ If || < ao_integrals_threshold then is zero @@ -331,20 +363,16 @@ Documentation Aligned n_pt_max_integrals -`n_pt_sup `_ +`n_pt_sup `_ Returns the upper boundary of the degree of the polynomial involved in the bielctronic integral : Ix(a_x,b_x,c_x,d_x) * Iy(a_y,b_y,c_y,d_y) * Iz(a_z,b_z,c_z,d_z) -`provide_all_mo_integrals `_ +`provide_all_mo_integrals `_ Undocumented -`pull_integrals `_ - How the collector pulls the computed integrals - - `push_integrals `_ Push integrals in the push socket diff --git a/src/Integrals_Monoelec/README.rst b/src/Integrals_Monoelec/README.rst index d92cea0a..7e926fd5 100644 --- a/src/Integrals_Monoelec/README.rst +++ b/src/Integrals_Monoelec/README.rst @@ -102,7 +102,7 @@ Documentation interaction nuclear electron -`ao_nucl_elec_integral_per_atom `_ +`ao_nucl_elec_integral_per_atom `_ ao_nucl_elec_integral_per_atom(i,j,k) = - where Rk is the geometry of the kth atom @@ -115,7 +115,7 @@ Documentation Local pseudo-potential -`ao_pseudo_integral_non_local `_ +`ao_pseudo_integral_non_local `_ Local pseudo-potential @@ -153,19 +153,19 @@ Documentation Undocumented -`give_polynom_mult_center_mono_elec `_ +`give_polynom_mult_center_mono_elec `_ Undocumented -`i_x1_pol_mult_mono_elec `_ +`i_x1_pol_mult_mono_elec `_ Undocumented -`i_x2_pol_mult_mono_elec `_ +`i_x2_pol_mult_mono_elec `_ Undocumented -`int_gaus_pol `_ +`int_gaus_pol `_ Undocumented @@ -200,7 +200,7 @@ Documentation interaction nuclear electron on the MO basis -`mo_nucl_elec_integral_per_atom `_ +`mo_nucl_elec_integral_per_atom `_ mo_nucl_elec_integral_per_atom(i,j,k) = - where Rk is the geometry of the kth atom @@ -227,7 +227,7 @@ Documentation array of the integrals of MO_i * z^2 MO_j -`nai_pol_mult `_ +`nai_pol_mult `_ Undocumented @@ -259,27 +259,27 @@ Documentation Undocumented -`pseudo_dz_k_transp `_ +`pseudo_dz_k_transp `_ Transposed arrays for pseudopotentials -`pseudo_dz_kl_transp `_ +`pseudo_dz_kl_transp `_ Transposed arrays for pseudopotentials -`pseudo_n_k_transp `_ +`pseudo_n_k_transp `_ Transposed arrays for pseudopotentials -`pseudo_n_kl_transp `_ +`pseudo_n_kl_transp `_ Transposed arrays for pseudopotentials -`pseudo_v_k_transp `_ +`pseudo_v_k_transp `_ Transposed arrays for pseudopotentials -`pseudo_v_kl_transp `_ +`pseudo_v_kl_transp `_ Transposed arrays for pseudopotentials @@ -299,23 +299,23 @@ Documentation Undocumented -`v_e_n `_ +`v_e_n `_ Undocumented -`v_phi `_ +`v_phi `_ Undocumented -`v_r `_ +`v_r `_ Undocumented -`v_theta `_ +`v_theta `_ Undocumented -`wallis `_ +`wallis `_ Undocumented diff --git a/src/Nuclei/README.rst b/src/Nuclei/README.rst index bf7e6f52..356b8e9e 100644 --- a/src/Nuclei/README.rst +++ b/src/Nuclei/README.rst @@ -38,7 +38,7 @@ Documentation Array of the name of element, sorted by nuclear charge (integer) -`nucl_charge `_ +`nucl_charge `_ Nuclear charges diff --git a/src/Utils/README.rst b/src/Utils/README.rst index 03ec80f5..902a5250 100644 --- a/src/Utils/README.rst +++ b/src/Utils/README.rst @@ -28,11 +28,11 @@ Documentation Compute 1st dimension such that it is aligned for vectorization. -`apply_rotation `_ +`apply_rotation `_ Apply the rotation found by find_rotation -`approx_dble `_ +`approx_dble `_ Undocumented @@ -55,19 +55,19 @@ Documentation Binomial coefficients -`dble_fact `_ +`dble_fact `_ Undocumented -`dble_fact_even `_ +`dble_fact_even `_ n!! -`dble_fact_odd `_ +`dble_fact_odd `_ n!! -`dble_logfact `_ +`dble_logfact `_ n!! @@ -93,6 +93,10 @@ Documentation contains the new order of the elements. +`dtranspose `_ + Transpose input matrix A into output matrix B + + `erf0 `_ Undocumented @@ -106,11 +110,11 @@ Documentation n! -`fact_inv `_ +`fact_inv `_ 1/n! -`find_rotation `_ +`find_rotation `_ Find A.C = B @@ -136,7 +140,7 @@ Documentation Undocumented -`get_pseudo_inverse `_ +`get_pseudo_inverse `_ Find C = A^-1 @@ -372,7 +376,7 @@ Documentation to be in integer*8 format -`inv_int `_ +`inv_int `_ 1/i @@ -408,7 +412,7 @@ Documentation contains the new order of the elements. -`lapack_diag `_ +`lapack_diag `_ Diagonalize matrix H .br H is untouched between input and ouptut @@ -419,7 +423,7 @@ Documentation .br -`lapack_diag_s2 `_ +`lapack_diag_s2 `_ Diagonalize matrix H .br H is untouched between input and ouptut @@ -430,7 +434,7 @@ Documentation .br -`lapack_diagd `_ +`lapack_diagd `_ Diagonalize matrix H .br H is untouched between input and ouptut @@ -441,7 +445,7 @@ Documentation .br -`lapack_partial_diag `_ +`lapack_partial_diag `_ Diagonalize matrix H .br H is untouched between input and ouptut @@ -452,25 +456,33 @@ Documentation .br -`logfact `_ +`logfact `_ n! -`lowercase `_ +`lowercase `_ Transform to lower case +`map_load_from_disk `_ + Undocumented + + +`map_save_to_disk `_ + Undocumented + + `multiply_poly `_ Multiply two polynomials D(t) =! D(t) +( B(t)*C(t)) -`normalize `_ +`normalize `_ Normalizes vector u u is expected to be aligned in memory. -`nproc `_ +`nproc `_ Number of current OpenMP threads @@ -492,7 +504,7 @@ Documentation .br -`ortho_lowdin `_ +`ortho_lowdin `_ Compute C_new=C_old.S^-1/2 orthogonalization. .br overlap : overlap matrix @@ -510,6 +522,19 @@ Documentation .br +`ortho_qr `_ + Orthogonalization using Q.R factorization + .br + A : matrix to orthogonalize + .br + LDA : leftmost dimension of A + .br + n : Number of rows of A + .br + m : Number of columns of A + .br + + `overlap_a_b_c `_ Undocumented @@ -607,7 +632,7 @@ Documentation to be in integer*8 format -`set_zero_extra_diag `_ +`set_zero_extra_diag `_ Undocumented @@ -634,18 +659,22 @@ Documentation .br -`u_dot_u `_ +`transpose `_ + Transpose input matrix A into output matrix B + + +`u_dot_u `_ Compute -`u_dot_v `_ +`u_dot_v `_ Compute -`wall_time `_ +`wall_time `_ The equivalent of cpu_time, but for the wall time. -`write_git_log `_ +`write_git_log `_ Write the last git commit in file iunit. diff --git a/src/ZMQ/README.rst b/src/ZMQ/README.rst index 187af23e..b73dc42d 100644 --- a/src/ZMQ/README.rst +++ b/src/ZMQ/README.rst @@ -21,59 +21,67 @@ Documentation .. by the `update_README.py` script. -`add_task_to_taskserver `_ +`add_task_to_taskserver `_ Get a task from the task server -`connect_to_taskserver `_ +`connect_to_taskserver `_ Connect to the task server and obtain the worker ID -`disconnect_from_taskserver `_ +`disconnect_from_taskserver `_ Disconnect from the task server -`end_parallel_job `_ +`end_parallel_job `_ End a new parallel job with name 'name'. The slave tasks execute subroutine 'slave' -`end_zmq_pair_socket `_ +`end_zmq_pair_socket `_ Terminate socket on which the results are sent. -`end_zmq_pull_socket `_ +`end_zmq_pull_socket `_ Terminate socket on which the results are sent. -`end_zmq_push_socket `_ +`end_zmq_push_socket `_ Terminate socket on which the results are sent. -`end_zmq_to_qp_run_socket `_ +`end_zmq_sub_socket `_ + Terminate socket on which the results are sent. + + +`end_zmq_to_qp_run_socket `_ Terminate the socket from the application to qp_run -`get_task_from_taskserver `_ +`get_task_from_taskserver `_ Get a task from the task server -`new_parallel_job `_ +`new_parallel_job `_ Start a new parallel job with name 'name'. The slave tasks execute subroutine 'slave' -`new_zmq_pair_socket `_ +`new_zmq_pair_socket `_ Socket on which the collector and the main communicate -`new_zmq_pull_socket `_ +`new_zmq_pull_socket `_ Socket on which the results are sent. If thread is 1, use inproc -`new_zmq_push_socket `_ +`new_zmq_push_socket `_ Socket on which the results are sent. If thread is 1, use inproc -`new_zmq_to_qp_run_socket `_ +`new_zmq_sub_socket `_ + Socket to read the state published by the Task server + + +`new_zmq_to_qp_run_socket `_ Socket on which the qp_run process replies @@ -82,29 +90,41 @@ Documentation Example : tcp://130.120.229.139:12345 -`reset_zmq_addresses `_ - Undocumented +`reset_zmq_addresses `_ + Socket which pulls the results (2) -`switch_qp_run_to_master `_ +`switch_qp_run_to_master `_ Address of the master qp_run socket Example : tcp://130.120.229.139:12345 -`task_done_to_taskserver `_ +`task_done_to_taskserver `_ Get a task from the task server +`wait_for_next_state `_ + Undocumented + + +`wait_for_state `_ + Wait for the ZMQ state to be ready + + +`wait_for_states `_ + Wait for the ZMQ state to be ready + + `zmq_context `_ Context for the ZeroMQ library -`zmq_delete_task `_ +`zmq_delete_task `_ When a task is done, it has to be removed from the list of tasks on the qp_run queue. This guarantees that the results have been received in the pull. -`zmq_port `_ +`zmq_port `_ Return the value of the ZMQ port from the corresponding integer @@ -113,6 +133,10 @@ Documentation Example : tcp://130.120.229.139:12345 +`zmq_set_running `_ + Set the job to Running in QP-run + + `zmq_socket_pair_inproc_address `_ Socket which pulls the results (2) @@ -133,6 +157,10 @@ Documentation Socket which pulls the results (2) -`zmq_state `_ +`zmq_socket_sub_tcp_address `_ + Socket which pulls the results (2) + + +`zmq_state `_ Threads executing work through the ZeroMQ interface diff --git a/src/ZMQ/tree_dependency.png b/src/ZMQ/tree_dependency.png new file mode 100644 index 00000000..e69de29b diff --git a/tests/bats/cassd.bats b/tests/bats/cassd.bats index 8cbc4899..2a8fabc2 100644 --- a/tests/bats/cassd.bats +++ b/tests/bats/cassd.bats @@ -13,7 +13,7 @@ source $QP_ROOT/tests/bats/common.bats.sh qp_set_mo_class $INPUT -core "[1]" -inact "[2,5]" -act "[3,4,6,7]" -virt "[8-24]" qp_run cassd_zmq $INPUT energy="$(ezfio get cas_sd_zmq energy_pt2)" - eq $energy -76.2311422169983 5.E-5 + eq $energy -76.231084536315 5.E-5 ezfio set determinants n_det_max 2048 ezfio set determinants read_wf True @@ -21,6 +21,6 @@ source $QP_ROOT/tests/bats/common.bats.sh qp_run cassd_zmq $INPUT ezfio set determinants read_wf False energy="$(ezfio get cas_sd_zmq energy)" - eq $energy -76.2300888408526 2.E-5 + eq $energy -76.2300887947446 2.E-5 } diff --git a/tests/bats/mrcepa0.bats b/tests/bats/mrcepa0.bats index e4b45f65..6bca8b7e 100644 --- a/tests/bats/mrcepa0.bats +++ b/tests/bats/mrcepa0.bats @@ -16,7 +16,7 @@ source $QP_ROOT/tests/bats/common.bats.sh ezfio set mrcepa0 n_it_max_dressed_ci 3 qp_run $EXE $INPUT energy="$(ezfio get mrcepa0 energy_pt2)" - eq $energy -76.23857580469 1.e-4 + eq $energy -76.2385617521816 1.e-4 } @test "MRCC H2O cc-pVDZ" { @@ -32,7 +32,7 @@ source $QP_ROOT/tests/bats/common.bats.sh ezfio set mrcepa0 n_it_max_dressed_ci 3 qp_run $EXE $INPUT energy="$(ezfio get mrcepa0 energy_pt2)" - eq $energy -76.238510023275 2.e-4 + eq $energy -76.2385052514433 2.e-4 } @test "MRSC2 H2O cc-pVDZ" { @@ -48,7 +48,7 @@ source $QP_ROOT/tests/bats/common.bats.sh ezfio set mrcepa0 n_it_max_dressed_ci 3 qp_run $EXE $INPUT energy="$(ezfio get mrcepa0 energy_pt2)" - eq $energy -76.2357889658142 2.e-4 + eq $energy -76.235786994991 2.e-4 } @test "MRCEPA0 H2O cc-pVDZ" { @@ -64,6 +64,6 @@ source $QP_ROOT/tests/bats/common.bats.sh ezfio set mrcepa0 n_it_max_dressed_ci 3 qp_run $EXE $INPUT energy="$(ezfio get mrcepa0 energy_pt2)" - eq $energy -76.2417748223423 2.e-4 + eq $energy -76.2417725924747 2.e-4 } diff --git a/tests/run_tests.sh b/tests/run_tests.sh index 367937f5..3ac452ad 100755 --- a/tests/run_tests.sh +++ b/tests/run_tests.sh @@ -1,16 +1,14 @@ -#!/bin/bash +#!/bin/bash -e LIST=" - convert.bats hf.bats -foboci.bats pseudo.bats fci.bats cassd.bats mrcepa0.bats - " +#foboci.bats export QP_PREFIX="timeout -s 9 600" @@ -30,10 +28,9 @@ do if [[ "$1" == "-v" ]] then echo "Verbose mode" - ./bats_to_sh.py $BATS_FILE | bash + ./bats_to_sh.py $BATS_FILE | bash else bats $BATS_FILE fi done -