diff --git a/.gitignore b/.gitignore index 096e385b..085e4f74 100644 --- a/.gitignore +++ b/.gitignore @@ -5,7 +5,10 @@ build.ninja .ninja_deps bin/ lib/ +lib64/ +libexec/ config/qp_create_ninja.pickle src/*/.gitignore ezfio_interface.irp.f share +*.swp diff --git a/RELEASE_NOTES.org b/RELEASE_NOTES.org index 7724d1d1..d9135cdf 100644 --- a/RELEASE_NOTES.org +++ b/RELEASE_NOTES.org @@ -4,90 +4,100 @@ ** Changes - - Python3 replaces Python2 - - Travis CI uses 3 jobs - - Moved Travis scripts into ~travis~ directory - - IRPF90 and EZFIO are now git submodules - - Now basis sets should be downloaded from basis-set-exchange website - - Added ~bse~ in the installable tools - - Documentation in ~src/README.rst~ - - Added two-body reduced density matrix - - Added basis set correction - - Added CAS-based on-top density functional - - Improve PT2 computation for excited-states: Mostly 2x2 - diagonalization, and some (n+1)x(n+1) diagonalizations - - Error bars for stochastic variance and norm of the perturbed wave function - - Improve PT2-matching for excited-states - - Compute the overlap of PT2 excited states - - Renamed SOP into CFG - - Improved parallelism in PT2 by splitting tasks - - Use max in multi-state PT2 instead of sum for the selection weight - - Added seniority - - Added excitation_max - - More tasks for distribueted Davidson - - Random guess vectors in Davidson have zeros to preserve symmetry - - Disk-based Davidson when too much memory is required - - Fixed bug in DIIS - - Fixed bug in molden (Au -> Angs) - - Fixed bug with non-contiguous MOs in active space and deleter MOs - - Complete network-free installation - - Fixed bug in selection when computing full PT2 - - Updated version of f77-zmq + - Python3 replaces Python2 + - Travis CI uses 3 jobs + - Moved Travis scripts into ~travis~ directory + - IRPF90 and EZFIO are now git submodules + - Now basis sets should be downloaded from basis-set-exchange website + - Added ~bse~ in the installable tools + - Documentation in ~src/README.rst~ + - Added two-body reduced density matrix + - Added basis set correction + - Added GTOs with complex exponent + - Added many types of integrals + - Added CAS-based on-top density functional + - Improve PT2 computation for excited-states: Mostly 2x2 + diagonalization, and some (n+1)x(n+1) diagonalizations + - Error bars for stochastic variance and norm of the perturbed wave function + - Improve PT2-matching for excited-states + - Compute the overlap of PT2 excited states + - Renamed SOP into CFG + - Improved parallelism in PT2 by splitting tasks + - Use max in multi-state PT2 instead of sum for the selection weight + - Added seniority + - Added excitation_max + - More tasks for distribueted Davidson + - Random guess vectors in Davidson have zeros to preserve symmetry + - Disk-based Davidson when too much memory is required + - Fixed bug in DIIS + - Fixed bug in molden (Au -> Angs) + - Fixed bug with non-contiguous MOs in active space and deleter MOs + - Complete network-free installation + - Fixed bug in selection when computing full PT2 + - Updated version of f77-zmq + - Added transcorrelated SCF + - Added transcorrelated CIPSI + - Started to introduce shells in AOs + - Added ECMD UEG functional + - Introduced DFT-based basis set correction + - General davidson algorithm -*** User interface +** User interface - - Added ~qp_basis~ script to install a basis set from the ~bse~ - command-line tool - - Introduced ~n_det_qp_edit~, ~psi_det_qp_edit~, and - ~psi_coef_qp_edit~ to accelerate the opening of qp_edit with - large wave functions - - Removed ~etc/ninja.rc~ - - Added flag to specify if the AOs are normalized - - Added flag to specify if the primitive Gaussians are normalized - - Added ~lin_dep_cutoff~, the cutoff for linear dependencies - - Davidson convergence threshold can be adapted from PT2 - - In ~density_for_dft~, ~no_core_density~ is now a logical - - Default for ~weight_selection~ has changed from 2 to 1 - - Nullify_small_elements in matrices to keep symmetry - - Default of density functional changed from LDA to PBE - - Added ~no_vvvv_integrals~ flag - - Added ~pt2_min_parallel_tasks~ to control parallelism in PT2 - - Added ~print_energy~ - - Added ~print_hamiltonian~ - - Added input for two body RDM - - Added keyword ~save_wf_after_selection~ - - Added a ~restore_symm~ flag to enforce the restoration of - symmetry in matrices - - qp_export_as_tgz exports also plugin codes - - Added a basis module containing basis set information - - Added qp_run truncate_wf + - Added ~qp_basis~ script to install a basis set from the ~bse~ + command-line tool + - Introduced ~n_det_qp_edit~, ~psi_det_qp_edit~, and + ~psi_coef_qp_edit~ to accelerate the opening of qp_edit with + large wave functions + - Removed ~etc/ninja.rc~ + - Added flag to specify if the AOs are normalized + - Added flag to specify if the primitive Gaussians are normalized + - Added ~lin_dep_cutoff~, the cutoff for linear dependencies + - Davidson convergence threshold can be adapted from PT2 + - In ~density_for_dft~, ~no_core_density~ is now a logical + - Default for ~weight_selection~ has changed from 2 to 1 + - Nullify_small_elements in matrices to keep symmetry + - Default of density functional changed from LDA to PBE + - Added ~no_vvvv_integrals~ flag + - Added ~pt2_min_parallel_tasks~ to control parallelism in PT2 + - Added ~print_energy~ + - Added ~print_hamiltonian~ + - Added input for two body RDM + - Added keyword ~save_wf_after_selection~ + - Added a ~restore_symm~ flag to enforce the restoration of + symmetry in matrices + - qp_export_as_tgz exports also plugin codes + - Added a basis module containing basis set information + - Added qp_run truncate_wf -*** Code +** Code - - Many bug fixes - - Changed electron-nucleus from ~e_n~ to ~n_e~ in names of variables - - Changed ~occ_pattern~ to ~configuration~ - - Replaced ~List.map~ by a tail-recursive version ~Qputils.list_map~ - - Added possible imaginary part in OCaml MO coefficients - - Added ~qp_clean_source_files.sh~ to remove non-ascii characters - - Added flag ~is_periodic~ for periodic systems - - Possibilities to handle complex integrals and complex MOs - - Moved pseuodpotential integrals out of ~ao_one_e_integrals~ - - Removed Schwarz test and added logical functions - ~ao_two_e_integral_zero~ and ~ao_one_e_integral_zero~ - - Introduced type for ~pt2_data~ - - Banned excitations are used with far apart localized MOs - - S_z2_Sz is now included in S2 - - S^2 in single precision - - Added Shank function - - Added utilities for periodic calculations - - Added ~V_ne_psi_energy~ - - Added ~h_core_guess~ routine - - Fixed Laplacians in real space (indices) - - Added LIB file to add extra libs in plugin - - Using Intel IPP for sorting when using Intel compiler - - Removed parallelism in sorting - - Compute banned_excitations from exchange integrals to accelerate with local MOs + - Many bug fixes + - Changed electron-nucleus from ~e_n~ to ~n_e~ in names of variables + - Changed ~occ_pattern~ to ~configuration~ + - Replaced ~List.map~ by a tail-recursive version ~Qputils.list_map~ + - Added possible imaginary part in OCaml MO coefficients + - Added ~qp_clean_source_files.sh~ to remove non-ascii characters + - Added flag ~is_periodic~ for periodic systems + - Possibilities to handle complex integrals and complex MOs + - Moved pseuodpotential integrals out of ~ao_one_e_integrals~ + - Removed Schwarz test and added logical functions + ~ao_two_e_integral_zero~ and ~ao_one_e_integral_zero~ + - Introduced type for ~pt2_data~ + - Banned excitations are used with far apart localized MOs + - S_z2_Sz is now included in S2 + - S^2 in single precision + - Added Shank function + - Added utilities for periodic calculations + - Added ~V_ne_psi_energy~ + - Added ~h_core_guess~ routine + - Fixed Laplacians in real space (indices) + - Added LIB file to add extra libs in plugin + - Using Intel IPP for sorting when using Intel compiler + - Removed parallelism in sorting + - Compute banned_excitations from exchange integrals to accelerate with local MOs + - Updated OCaml for 4.13 + diff --git a/configure b/configure index 7d063b1d..79cd7119 100755 --- a/configure +++ b/configure @@ -246,7 +246,7 @@ EOF execute << EOF cd "\${QP_ROOT}"/external tar --gunzip --extract --file qp2-dependencies/bse-v0.8.11.tar.gz - pip install -e basis_set_exchange-* + python3 -m pip install -e basis_set_exchange-* EOF elif [[ ${PACKAGE} = zlib ]] ; then @@ -303,6 +303,25 @@ fi ZEROMQ=$(find_lib -lzmq) if [[ ${ZEROMQ} = $(not_found) ]] ; then + + MAKE=$(find_exe make) + if [[ ${MAKE} = $(not_found) ]] ; then + error "make is not installed." + fail + fi + + M4=$(find_exe autoreconf) + if [[ ${M4} = $(not_found) ]] ; then + error "autoreconf is not installed." + fail + fi + + M4=$(find_exe m4) + if [[ ${M4} = $(not_found) ]] ; then + error "m4 preprocesssor is not installed." + fail + fi + error "ZeroMQ (zeromq) is not installed." fail fi diff --git a/etc/qp.rc b/etc/qp.rc index 064ca3f7..c56661c7 100644 --- a/etc/qp.rc +++ b/etc/qp.rc @@ -80,6 +80,8 @@ function qp() if [[ -d $NAME ]] ; then [[ -d $EZFIO_FILE ]] && ezfio unset_file ezfio set_file $NAME + else + qp_create_ezfio -h | more fi unset _ARGS ;; diff --git a/external/.gitignore b/external/.gitignore index 676c79b7..241e560d 100644 --- a/external/.gitignore +++ b/external/.gitignore @@ -1,2 +1,2 @@ -#* +* diff --git a/external/qp2-dependencies b/external/qp2-dependencies index 242151e0..90ee61f5 160000 --- a/external/qp2-dependencies +++ b/external/qp2-dependencies @@ -1 +1 @@ -Subproject commit 242151e03d1d6bf042387226431d82d35845686a +Subproject commit 90ee61f5041c7c94a0c605625a264860292813a0 diff --git a/include/.gitignore b/include/.gitignore index afcb2de6..72e8ffc0 100644 --- a/include/.gitignore +++ b/include/.gitignore @@ -1,8 +1 @@ -zmq.h -gmp.h -zconf.h -zconf.h -zlib.h -zmq_utils.h -f77_zmq_free.h -f77_zmq.h +* diff --git a/ocaml/Command_line.ml b/ocaml/Command_line.ml index 1dd57892..602315c6 100644 --- a/ocaml/Command_line.ml +++ b/ocaml/Command_line.ml @@ -1,3 +1,5 @@ +exception Error of string + type short_opt = char type long_opt = string type optional = Mandatory | Optional @@ -181,15 +183,16 @@ let set_specs specs_in = Getopt.parse_cmdline cmd_specs (fun x -> anon_args := !anon_args @ [x]); if show_help () then - (help () ; exit 0); + help () + else + (* Check that all mandatory arguments are set *) + List.filter (fun x -> x.short <> ' ' && x.opt = Mandatory) !specs + |> List.iter (fun x -> + match get x.long with + | Some _ -> () + | None -> raise (Error ("--"^x.long^" option is missing.")) + ) - (* Check that all mandatory arguments are set *) - List.filter (fun x -> x.short <> ' ' && x.opt = Mandatory) !specs - |> List.iter (fun x -> - match get x.long with - | Some _ -> () - | None -> failwith ("Error: --"^x.long^" option is missing.") - ) ;; diff --git a/ocaml/Command_line.mli b/ocaml/Command_line.mli index 9f6e7022..5ad4ee08 100644 --- a/ocaml/Command_line.mli +++ b/ocaml/Command_line.mli @@ -59,6 +59,8 @@ let () = *) +exception Error of string + type short_opt = char type long_opt = string diff --git a/ocaml/Input_ao_two_e_eff_pot.ml b/ocaml/Input_ao_two_e_eff_pot.ml deleted file mode 100644 index e4e2c059..00000000 --- a/ocaml/Input_ao_two_e_eff_pot.ml +++ /dev/null @@ -1,113 +0,0 @@ -(* =~=~ *) -(* Init *) -(* =~=~ *) - -open Qptypes;; -open Qputils;; -open Sexplib.Std;; - -module Ao_two_e_eff_pot : sig -(* Generate type *) - type t = - { - adjoint_tc_h : bool; - grad_squared : bool; - } [@@deriving sexp] - ;; - val read : unit -> t option - val write : t-> unit - val to_string : t -> string - val to_rst : t -> Rst_string.t - val of_rst : Rst_string.t -> t option -end = struct -(* Generate type *) - type t = - { - adjoint_tc_h : bool; - grad_squared : bool; - } [@@deriving sexp] - ;; - - let get_default = Qpackage.get_ezfio_default "ao_two_e_eff_pot";; - -(* =~=~=~=~=~=~==~=~=~=~=~=~ *) -(* Generate Special Function *) -(* =~=~=~==~=~~=~=~=~=~=~=~=~ *) - -(* Read snippet for adjoint_tc_h *) - let read_adjoint_tc_h () = - if not (Ezfio.has_ao_two_e_eff_pot_adjoint_tc_h ()) then - get_default "adjoint_tc_h" - |> bool_of_string - |> Ezfio.set_ao_two_e_eff_pot_adjoint_tc_h - ; - Ezfio.get_ao_two_e_eff_pot_adjoint_tc_h () - ;; -(* Write snippet for adjoint_tc_h *) - let write_adjoint_tc_h = - Ezfio.set_ao_two_e_eff_pot_adjoint_tc_h - ;; - -(* Read snippet for grad_squared *) - let read_grad_squared () = - if not (Ezfio.has_ao_two_e_eff_pot_grad_squared ()) then - get_default "grad_squared" - |> bool_of_string - |> Ezfio.set_ao_two_e_eff_pot_grad_squared - ; - Ezfio.get_ao_two_e_eff_pot_grad_squared () - ;; -(* Write snippet for grad_squared *) - let write_grad_squared = - Ezfio.set_ao_two_e_eff_pot_grad_squared - ;; - -(* =~=~=~=~=~=~=~=~=~=~=~=~ *) -(* Generate Global Function *) -(* =~=~=~=~=~=~=~=~=~=~=~=~ *) - -(* Read all *) - let read() = - Some - { - adjoint_tc_h = read_adjoint_tc_h (); - grad_squared = read_grad_squared (); - } - ;; -(* Write all *) - let write{ - adjoint_tc_h; - grad_squared; - } = - write_adjoint_tc_h adjoint_tc_h; - write_grad_squared grad_squared; - ;; -(* to_string*) - let to_string b = - Printf.sprintf " - adjoint_tc_h = %s - grad_squared = %s - " - (string_of_bool b.adjoint_tc_h) - (string_of_bool b.grad_squared) - ;; -(* to_rst*) - let to_rst b = - Printf.sprintf " - If |true|, you compute the adjoint of the transcorrelated Hamiltonian :: - - adjoint_tc_h = %s - - If |true|, you compute also the square of the gradient of the correlation factor :: - - grad_squared = %s - - " - (string_of_bool b.adjoint_tc_h) - (string_of_bool b.grad_squared) - |> Rst_string.of_string - ;; - include Generic_input_of_rst;; - let of_rst = of_rst t_of_sexp;; - -end \ No newline at end of file diff --git a/ocaml/Input_bi_ortho_mos.ml b/ocaml/Input_bi_ortho_mos.ml deleted file mode 100644 index 5523a589..00000000 --- a/ocaml/Input_bi_ortho_mos.ml +++ /dev/null @@ -1,87 +0,0 @@ -(* =~=~ *) -(* Init *) -(* =~=~ *) - -open Qptypes;; -open Qputils;; -open Sexplib.Std;; - -module Bi_ortho_mos : sig -(* Generate type *) - type t = - { - bi_ortho : bool; - } [@@deriving sexp] - ;; - val read : unit -> t option - val write : t-> unit - val to_string : t -> string - val to_rst : t -> Rst_string.t - val of_rst : Rst_string.t -> t option -end = struct -(* Generate type *) - type t = - { - bi_ortho : bool; - } [@@deriving sexp] - ;; - - let get_default = Qpackage.get_ezfio_default "bi_ortho_mos";; - -(* =~=~=~=~=~=~==~=~=~=~=~=~ *) -(* Generate Special Function *) -(* =~=~=~==~=~~=~=~=~=~=~=~=~ *) - -(* Read snippet for bi_ortho *) - let read_bi_ortho () = - if not (Ezfio.has_bi_ortho_mos_bi_ortho ()) then - get_default "bi_ortho" - |> bool_of_string - |> Ezfio.set_bi_ortho_mos_bi_ortho - ; - Ezfio.get_bi_ortho_mos_bi_ortho () - ;; -(* Write snippet for bi_ortho *) - let write_bi_ortho = - Ezfio.set_bi_ortho_mos_bi_ortho - ;; - -(* =~=~=~=~=~=~=~=~=~=~=~=~ *) -(* Generate Global Function *) -(* =~=~=~=~=~=~=~=~=~=~=~=~ *) - -(* Read all *) - let read() = - Some - { - bi_ortho = read_bi_ortho (); - } - ;; -(* Write all *) - let write{ - bi_ortho; - } = - write_bi_ortho bi_ortho; - ;; -(* to_string*) - let to_string b = - Printf.sprintf " - bi_ortho = %s - " - (string_of_bool b.bi_ortho) - ;; -(* to_rst*) - let to_rst b = - Printf.sprintf " - If |true|, the MO basis is assumed to be bi-orthonormal :: - - bi_ortho = %s - - " - (string_of_bool b.bi_ortho) - |> Rst_string.of_string - ;; - include Generic_input_of_rst;; - let of_rst = of_rst t_of_sexp;; - -end \ No newline at end of file diff --git a/ocaml/Input_cassd.ml b/ocaml/Input_cassd.ml deleted file mode 100644 index 03416f42..00000000 --- a/ocaml/Input_cassd.ml +++ /dev/null @@ -1,113 +0,0 @@ -(* =~=~ *) -(* Init *) -(* =~=~ *) - -open Qptypes;; -open Qputils;; -open Sexplib.Std;; - -module Cassd : sig -(* Generate type *) - type t = - { - do_ddci : bool; - do_only_1h1p : bool; - } [@@deriving sexp] - ;; - val read : unit -> t option - val write : t-> unit - val to_string : t -> string - val to_rst : t -> Rst_string.t - val of_rst : Rst_string.t -> t option -end = struct -(* Generate type *) - type t = - { - do_ddci : bool; - do_only_1h1p : bool; - } [@@deriving sexp] - ;; - - let get_default = Qpackage.get_ezfio_default "cassd";; - -(* =~=~=~=~=~=~==~=~=~=~=~=~ *) -(* Generate Special Function *) -(* =~=~=~==~=~~=~=~=~=~=~=~=~ *) - -(* Read snippet for do_ddci *) - let read_do_ddci () = - if not (Ezfio.has_cassd_do_ddci ()) then - get_default "do_ddci" - |> bool_of_string - |> Ezfio.set_cassd_do_ddci - ; - Ezfio.get_cassd_do_ddci () - ;; -(* Write snippet for do_ddci *) - let write_do_ddci = - Ezfio.set_cassd_do_ddci - ;; - -(* Read snippet for do_only_1h1p *) - let read_do_only_1h1p () = - if not (Ezfio.has_cassd_do_only_1h1p ()) then - get_default "do_only_1h1p" - |> bool_of_string - |> Ezfio.set_cassd_do_only_1h1p - ; - Ezfio.get_cassd_do_only_1h1p () - ;; -(* Write snippet for do_only_1h1p *) - let write_do_only_1h1p = - Ezfio.set_cassd_do_only_1h1p - ;; - -(* =~=~=~=~=~=~=~=~=~=~=~=~ *) -(* Generate Global Function *) -(* =~=~=~=~=~=~=~=~=~=~=~=~ *) - -(* Read all *) - let read() = - Some - { - do_ddci = read_do_ddci (); - do_only_1h1p = read_do_only_1h1p (); - } - ;; -(* Write all *) - let write{ - do_ddci; - do_only_1h1p; - } = - write_do_ddci do_ddci; - write_do_only_1h1p do_only_1h1p; - ;; -(* to_string*) - let to_string b = - Printf.sprintf " - do_ddci = %s - do_only_1h1p = %s - " - (string_of_bool b.do_ddci) - (string_of_bool b.do_only_1h1p) - ;; -(* to_rst*) - let to_rst b = - Printf.sprintf " - If true, remove purely inactive double excitations :: - - do_ddci = %s - - If true, do only one hole/one particle excitations :: - - do_only_1h1p = %s - - " - (string_of_bool b.do_ddci) - (string_of_bool b.do_only_1h1p) - |> Rst_string.of_string - ;; - include Generic_input_of_rst;; - let of_rst = of_rst t_of_sexp;; - -end \ No newline at end of file diff --git a/ocaml/Input_cipsi_deb.ml b/ocaml/Input_cipsi_deb.ml deleted file mode 100644 index 9849b0e2..00000000 --- a/ocaml/Input_cipsi_deb.ml +++ /dev/null @@ -1,243 +0,0 @@ -(* =~=~ *) -(* Init *) -(* =~=~ *) - -open Qptypes;; -open Qputils;; -open Sexplib.Std;; - -module Cipsi_deb : sig -(* Generate type *) - type t = - { - pert_2rdm : bool; - save_wf_after_selection : bool; - seniority_max : int; - excitation_ref : int; - excitation_max : int; - excitation_alpha_max : int; - excitation_beta_max : int; - } [@@deriving sexp] - ;; - val read : unit -> t option - val write : t-> unit - val to_string : t -> string - val to_rst : t -> Rst_string.t - val of_rst : Rst_string.t -> t option -end = struct -(* Generate type *) - type t = - { - pert_2rdm : bool; - save_wf_after_selection : bool; - seniority_max : int; - excitation_ref : int; - excitation_max : int; - excitation_alpha_max : int; - excitation_beta_max : int; - } [@@deriving sexp] - ;; - - let get_default = Qpackage.get_ezfio_default "cipsi_deb";; - -(* =~=~=~=~=~=~==~=~=~=~=~=~ *) -(* Generate Special Function *) -(* =~=~=~==~=~~=~=~=~=~=~=~=~ *) - -(* Read snippet for excitation_alpha_max *) - let read_excitation_alpha_max () = - if not (Ezfio.has_cipsi_deb_excitation_alpha_max ()) then - get_default "excitation_alpha_max" - |> int_of_string - |> Ezfio.set_cipsi_deb_excitation_alpha_max - ; - Ezfio.get_cipsi_deb_excitation_alpha_max () - ;; -(* Write snippet for excitation_alpha_max *) - let write_excitation_alpha_max = - Ezfio.set_cipsi_deb_excitation_alpha_max - ;; - -(* Read snippet for excitation_beta_max *) - let read_excitation_beta_max () = - if not (Ezfio.has_cipsi_deb_excitation_beta_max ()) then - get_default "excitation_beta_max" - |> int_of_string - |> Ezfio.set_cipsi_deb_excitation_beta_max - ; - Ezfio.get_cipsi_deb_excitation_beta_max () - ;; -(* Write snippet for excitation_beta_max *) - let write_excitation_beta_max = - Ezfio.set_cipsi_deb_excitation_beta_max - ;; - -(* Read snippet for excitation_max *) - let read_excitation_max () = - if not (Ezfio.has_cipsi_deb_excitation_max ()) then - get_default "excitation_max" - |> int_of_string - |> Ezfio.set_cipsi_deb_excitation_max - ; - Ezfio.get_cipsi_deb_excitation_max () - ;; -(* Write snippet for excitation_max *) - let write_excitation_max = - Ezfio.set_cipsi_deb_excitation_max - ;; - -(* Read snippet for excitation_ref *) - let read_excitation_ref () = - if not (Ezfio.has_cipsi_deb_excitation_ref ()) then - get_default "excitation_ref" - |> int_of_string - |> Ezfio.set_cipsi_deb_excitation_ref - ; - Ezfio.get_cipsi_deb_excitation_ref () - ;; -(* Write snippet for excitation_ref *) - let write_excitation_ref = - Ezfio.set_cipsi_deb_excitation_ref - ;; - -(* Read snippet for pert_2rdm *) - let read_pert_2rdm () = - if not (Ezfio.has_cipsi_deb_pert_2rdm ()) then - get_default "pert_2rdm" - |> bool_of_string - |> Ezfio.set_cipsi_deb_pert_2rdm - ; - Ezfio.get_cipsi_deb_pert_2rdm () - ;; -(* Write snippet for pert_2rdm *) - let write_pert_2rdm = - Ezfio.set_cipsi_deb_pert_2rdm - ;; - -(* Read snippet for save_wf_after_selection *) - let read_save_wf_after_selection () = - if not (Ezfio.has_cipsi_deb_save_wf_after_selection ()) then - get_default "save_wf_after_selection" - |> bool_of_string - |> Ezfio.set_cipsi_deb_save_wf_after_selection - ; - Ezfio.get_cipsi_deb_save_wf_after_selection () - ;; -(* Write snippet for save_wf_after_selection *) - let write_save_wf_after_selection = - Ezfio.set_cipsi_deb_save_wf_after_selection - ;; - -(* Read snippet for seniority_max *) - let read_seniority_max () = - if not (Ezfio.has_cipsi_deb_seniority_max ()) then - get_default "seniority_max" - |> int_of_string - |> Ezfio.set_cipsi_deb_seniority_max - ; - Ezfio.get_cipsi_deb_seniority_max () - ;; -(* Write snippet for seniority_max *) - let write_seniority_max = - Ezfio.set_cipsi_deb_seniority_max - ;; - -(* =~=~=~=~=~=~=~=~=~=~=~=~ *) -(* Generate Global Function *) -(* =~=~=~=~=~=~=~=~=~=~=~=~ *) - -(* Read all *) - let read() = - Some - { - pert_2rdm = read_pert_2rdm (); - save_wf_after_selection = read_save_wf_after_selection (); - seniority_max = read_seniority_max (); - excitation_ref = read_excitation_ref (); - excitation_max = read_excitation_max (); - excitation_alpha_max = read_excitation_alpha_max (); - excitation_beta_max = read_excitation_beta_max (); - } - ;; -(* Write all *) - let write{ - pert_2rdm; - save_wf_after_selection; - seniority_max; - excitation_ref; - excitation_max; - excitation_alpha_max; - excitation_beta_max; - } = - write_pert_2rdm pert_2rdm; - write_save_wf_after_selection save_wf_after_selection; - write_seniority_max seniority_max; - write_excitation_ref excitation_ref; - write_excitation_max excitation_max; - write_excitation_alpha_max excitation_alpha_max; - write_excitation_beta_max excitation_beta_max; - ;; -(* to_string*) - let to_string b = - Printf.sprintf " - pert_2rdm = %s - save_wf_after_selection = %s - seniority_max = %s - excitation_ref = %s - excitation_max = %s - excitation_alpha_max = %s - excitation_beta_max = %s - " - (string_of_bool b.pert_2rdm) - (string_of_bool b.save_wf_after_selection) - (string_of_int b.seniority_max) - (string_of_int b.excitation_ref) - (string_of_int b.excitation_max) - (string_of_int b.excitation_alpha_max) - (string_of_int b.excitation_beta_max) - ;; -(* to_rst*) - let to_rst b = - Printf.sprintf " - If true, computes the one- and two-body rdms with perturbation theory :: - - pert_2rdm = %s - - If true, saves the wave function after the selection, before the diagonalization :: - - save_wf_after_selection = %s - - Maximum number of allowed open shells. Using -1 selects all determinants :: - - seniority_max = %s - - 1: Hartree-Fock determinant, 2:All determinants of the dominant configuration :: - - excitation_ref = %s - - Maximum number of excitation with respect to the Hartree-Fock determinant. Using -1 selects all determinants :: - - excitation_max = %s - - Maximum number of excitation for alpha determinants with respect to the Hartree-Fock determinant. Using -1 selects all determinants :: - - excitation_alpha_max = %s - - Maximum number of excitation for beta determinants with respect to the Hartree-Fock determinant. Using -1 selects all determinants :: - - excitation_beta_max = %s - - " - (string_of_bool b.pert_2rdm) - (string_of_bool b.save_wf_after_selection) - (string_of_int b.seniority_max) - (string_of_int b.excitation_ref) - (string_of_int b.excitation_max) - (string_of_int b.excitation_alpha_max) - (string_of_int b.excitation_beta_max) - |> Rst_string.of_string - ;; - include Generic_input_of_rst;; - let of_rst = of_rst t_of_sexp;; - -end \ No newline at end of file diff --git a/ocaml/Input_tc_h_clean.ml b/ocaml/Input_tc_h_clean.ml deleted file mode 100644 index 2fd145fa..00000000 --- a/ocaml/Input_tc_h_clean.ml +++ /dev/null @@ -1,351 +0,0 @@ -(* =~=~ *) -(* Init *) -(* =~=~ *) - -open Qptypes;; -open Qputils;; -open Sexplib.Std;; - -module Tc_h_clean : sig -(* Generate type *) - type t = - { - read_rl_eigv : bool; - comp_left_eigv : bool; - three_body_h_tc : bool; - pure_three_body_h_tc : bool; - double_normal_ord : bool; - core_tc_op : bool; - full_tc_h_solver : bool; - thresh_it_dav : Threshold.t; - max_it_dav : int; - thresh_psi_r : Threshold.t; - thresh_psi_r_norm : bool; - } [@@deriving sexp] - ;; - val read : unit -> t option - val write : t-> unit - val to_string : t -> string - val to_rst : t -> Rst_string.t - val of_rst : Rst_string.t -> t option -end = struct -(* Generate type *) - type t = - { - read_rl_eigv : bool; - comp_left_eigv : bool; - three_body_h_tc : bool; - pure_three_body_h_tc : bool; - double_normal_ord : bool; - core_tc_op : bool; - full_tc_h_solver : bool; - thresh_it_dav : Threshold.t; - max_it_dav : int; - thresh_psi_r : Threshold.t; - thresh_psi_r_norm : bool; - } [@@deriving sexp] - ;; - - let get_default = Qpackage.get_ezfio_default "tc_h_clean";; - -(* =~=~=~=~=~=~==~=~=~=~=~=~ *) -(* Generate Special Function *) -(* =~=~=~==~=~~=~=~=~=~=~=~=~ *) - -(* Read snippet for comp_left_eigv *) - let read_comp_left_eigv () = - if not (Ezfio.has_tc_h_clean_comp_left_eigv ()) then - get_default "comp_left_eigv" - |> bool_of_string - |> Ezfio.set_tc_h_clean_comp_left_eigv - ; - Ezfio.get_tc_h_clean_comp_left_eigv () - ;; -(* Write snippet for comp_left_eigv *) - let write_comp_left_eigv = - Ezfio.set_tc_h_clean_comp_left_eigv - ;; - -(* Read snippet for core_tc_op *) - let read_core_tc_op () = - if not (Ezfio.has_tc_h_clean_core_tc_op ()) then - get_default "core_tc_op" - |> bool_of_string - |> Ezfio.set_tc_h_clean_core_tc_op - ; - Ezfio.get_tc_h_clean_core_tc_op () - ;; -(* Write snippet for core_tc_op *) - let write_core_tc_op = - Ezfio.set_tc_h_clean_core_tc_op - ;; - -(* Read snippet for double_normal_ord *) - let read_double_normal_ord () = - if not (Ezfio.has_tc_h_clean_double_normal_ord ()) then - get_default "double_normal_ord" - |> bool_of_string - |> Ezfio.set_tc_h_clean_double_normal_ord - ; - Ezfio.get_tc_h_clean_double_normal_ord () - ;; -(* Write snippet for double_normal_ord *) - let write_double_normal_ord = - Ezfio.set_tc_h_clean_double_normal_ord - ;; - -(* Read snippet for full_tc_h_solver *) - let read_full_tc_h_solver () = - if not (Ezfio.has_tc_h_clean_full_tc_h_solver ()) then - get_default "full_tc_h_solver" - |> bool_of_string - |> Ezfio.set_tc_h_clean_full_tc_h_solver - ; - Ezfio.get_tc_h_clean_full_tc_h_solver () - ;; -(* Write snippet for full_tc_h_solver *) - let write_full_tc_h_solver = - Ezfio.set_tc_h_clean_full_tc_h_solver - ;; - -(* Read snippet for max_it_dav *) - let read_max_it_dav () = - if not (Ezfio.has_tc_h_clean_max_it_dav ()) then - get_default "max_it_dav" - |> int_of_string - |> Ezfio.set_tc_h_clean_max_it_dav - ; - Ezfio.get_tc_h_clean_max_it_dav () - ;; -(* Write snippet for max_it_dav *) - let write_max_it_dav = - Ezfio.set_tc_h_clean_max_it_dav - ;; - -(* Read snippet for pure_three_body_h_tc *) - let read_pure_three_body_h_tc () = - if not (Ezfio.has_tc_h_clean_pure_three_body_h_tc ()) then - get_default "pure_three_body_h_tc" - |> bool_of_string - |> Ezfio.set_tc_h_clean_pure_three_body_h_tc - ; - Ezfio.get_tc_h_clean_pure_three_body_h_tc () - ;; -(* Write snippet for pure_three_body_h_tc *) - let write_pure_three_body_h_tc = - Ezfio.set_tc_h_clean_pure_three_body_h_tc - ;; - -(* Read snippet for read_rl_eigv *) - let read_read_rl_eigv () = - if not (Ezfio.has_tc_h_clean_read_rl_eigv ()) then - get_default "read_rl_eigv" - |> bool_of_string - |> Ezfio.set_tc_h_clean_read_rl_eigv - ; - Ezfio.get_tc_h_clean_read_rl_eigv () - ;; -(* Write snippet for read_rl_eigv *) - let write_read_rl_eigv = - Ezfio.set_tc_h_clean_read_rl_eigv - ;; - -(* Read snippet for three_body_h_tc *) - let read_three_body_h_tc () = - if not (Ezfio.has_tc_h_clean_three_body_h_tc ()) then - get_default "three_body_h_tc" - |> bool_of_string - |> Ezfio.set_tc_h_clean_three_body_h_tc - ; - Ezfio.get_tc_h_clean_three_body_h_tc () - ;; -(* Write snippet for three_body_h_tc *) - let write_three_body_h_tc = - Ezfio.set_tc_h_clean_three_body_h_tc - ;; - -(* Read snippet for thresh_it_dav *) - let read_thresh_it_dav () = - if not (Ezfio.has_tc_h_clean_thresh_it_dav ()) then - get_default "thresh_it_dav" - |> float_of_string - |> Ezfio.set_tc_h_clean_thresh_it_dav - ; - Ezfio.get_tc_h_clean_thresh_it_dav () - |> Threshold.of_float - ;; -(* Write snippet for thresh_it_dav *) - let write_thresh_it_dav var = - Threshold.to_float var - |> Ezfio.set_tc_h_clean_thresh_it_dav - ;; - -(* Read snippet for thresh_psi_r *) - let read_thresh_psi_r () = - if not (Ezfio.has_tc_h_clean_thresh_psi_r ()) then - get_default "thresh_psi_r" - |> float_of_string - |> Ezfio.set_tc_h_clean_thresh_psi_r - ; - Ezfio.get_tc_h_clean_thresh_psi_r () - |> Threshold.of_float - ;; -(* Write snippet for thresh_psi_r *) - let write_thresh_psi_r var = - Threshold.to_float var - |> Ezfio.set_tc_h_clean_thresh_psi_r - ;; - -(* Read snippet for thresh_psi_r_norm *) - let read_thresh_psi_r_norm () = - if not (Ezfio.has_tc_h_clean_thresh_psi_r_norm ()) then - get_default "thresh_psi_r_norm" - |> bool_of_string - |> Ezfio.set_tc_h_clean_thresh_psi_r_norm - ; - Ezfio.get_tc_h_clean_thresh_psi_r_norm () - ;; -(* Write snippet for thresh_psi_r_norm *) - let write_thresh_psi_r_norm = - Ezfio.set_tc_h_clean_thresh_psi_r_norm - ;; - -(* =~=~=~=~=~=~=~=~=~=~=~=~ *) -(* Generate Global Function *) -(* =~=~=~=~=~=~=~=~=~=~=~=~ *) - -(* Read all *) - let read() = - Some - { - read_rl_eigv = read_read_rl_eigv (); - comp_left_eigv = read_comp_left_eigv (); - three_body_h_tc = read_three_body_h_tc (); - pure_three_body_h_tc = read_pure_three_body_h_tc (); - double_normal_ord = read_double_normal_ord (); - core_tc_op = read_core_tc_op (); - full_tc_h_solver = read_full_tc_h_solver (); - thresh_it_dav = read_thresh_it_dav (); - max_it_dav = read_max_it_dav (); - thresh_psi_r = read_thresh_psi_r (); - thresh_psi_r_norm = read_thresh_psi_r_norm (); - } - ;; -(* Write all *) - let write{ - read_rl_eigv; - comp_left_eigv; - three_body_h_tc; - pure_three_body_h_tc; - double_normal_ord; - core_tc_op; - full_tc_h_solver; - thresh_it_dav; - max_it_dav; - thresh_psi_r; - thresh_psi_r_norm; - } = - write_read_rl_eigv read_rl_eigv; - write_comp_left_eigv comp_left_eigv; - write_three_body_h_tc three_body_h_tc; - write_pure_three_body_h_tc pure_three_body_h_tc; - write_double_normal_ord double_normal_ord; - write_core_tc_op core_tc_op; - write_full_tc_h_solver full_tc_h_solver; - write_thresh_it_dav thresh_it_dav; - write_max_it_dav max_it_dav; - write_thresh_psi_r thresh_psi_r; - write_thresh_psi_r_norm thresh_psi_r_norm; - ;; -(* to_string*) - let to_string b = - Printf.sprintf " - read_rl_eigv = %s - comp_left_eigv = %s - three_body_h_tc = %s - pure_three_body_h_tc = %s - double_normal_ord = %s - core_tc_op = %s - full_tc_h_solver = %s - thresh_it_dav = %s - max_it_dav = %s - thresh_psi_r = %s - thresh_psi_r_norm = %s - " - (string_of_bool b.read_rl_eigv) - (string_of_bool b.comp_left_eigv) - (string_of_bool b.three_body_h_tc) - (string_of_bool b.pure_three_body_h_tc) - (string_of_bool b.double_normal_ord) - (string_of_bool b.core_tc_op) - (string_of_bool b.full_tc_h_solver) - (Threshold.to_string b.thresh_it_dav) - (string_of_int b.max_it_dav) - (Threshold.to_string b.thresh_psi_r) - (string_of_bool b.thresh_psi_r_norm) - ;; -(* to_rst*) - let to_rst b = - Printf.sprintf " - If |true|, read the right/left eigenvectors from ezfio :: - - read_rl_eigv = %s - - If |true|, computes also the left-eigenvector :: - - comp_left_eigv = %s - - If |true|, three-body terms are included :: - - three_body_h_tc = %s - - If |true|, pure triple excitation three-body terms are included :: - - pure_three_body_h_tc = %s - - If |true|, contracted double excitation three-body terms are included :: - - double_normal_ord = %s - - If |true|, takes the usual Hamiltonian for core orbitals (assumed to be doubly occupied) :: - - core_tc_op = %s - - If |true|, you diagonalize the full TC H matrix :: - - full_tc_h_solver = %s - - Thresholds on the energy for iterative Davidson used in TC :: - - thresh_it_dav = %s - - nb max of iteration in Davidson used in TC :: - - max_it_dav = %s - - Thresholds on the coefficients of the right-eigenvector. Used for PT2 computation. :: - - thresh_psi_r = %s - - If |true|, you prune the WF to compute the PT1 coef based on the norm. If False, the pruning is done through the amplitude on the right-coefficient. :: - - thresh_psi_r_norm = %s - - " - (string_of_bool b.read_rl_eigv) - (string_of_bool b.comp_left_eigv) - (string_of_bool b.three_body_h_tc) - (string_of_bool b.pure_three_body_h_tc) - (string_of_bool b.double_normal_ord) - (string_of_bool b.core_tc_op) - (string_of_bool b.full_tc_h_solver) - (Threshold.to_string b.thresh_it_dav) - (string_of_int b.max_it_dav) - (Threshold.to_string b.thresh_psi_r) - (string_of_bool b.thresh_psi_r_norm) - |> Rst_string.of_string - ;; - include Generic_input_of_rst;; - let of_rst = of_rst t_of_sexp;; - -end \ No newline at end of file diff --git a/ocaml/Input_tc_scf.ml b/ocaml/Input_tc_scf.ml deleted file mode 100644 index 2a709716..00000000 --- a/ocaml/Input_tc_scf.ml +++ /dev/null @@ -1,143 +0,0 @@ -(* =~=~ *) -(* Init *) -(* =~=~ *) - -open Qptypes;; -open Qputils;; -open Sexplib.Std;; - -module Tc_scf : sig -(* Generate type *) - type t = - { - bi_ortho : bool; - thresh_tcscf : Threshold.t; - n_it_tcscf_max : Strictly_positive_int.t; - } [@@deriving sexp] - ;; - val read : unit -> t option - val write : t-> unit - val to_string : t -> string - val to_rst : t -> Rst_string.t - val of_rst : Rst_string.t -> t option -end = struct -(* Generate type *) - type t = - { - bi_ortho : bool; - thresh_tcscf : Threshold.t; - n_it_tcscf_max : Strictly_positive_int.t; - } [@@deriving sexp] - ;; - - let get_default = Qpackage.get_ezfio_default "tc_scf";; - -(* =~=~=~=~=~=~==~=~=~=~=~=~ *) -(* Generate Special Function *) -(* =~=~=~==~=~~=~=~=~=~=~=~=~ *) - -(* Read snippet for bi_ortho *) - let read_bi_ortho () = - if not (Ezfio.has_tc_scf_bi_ortho ()) then - get_default "bi_ortho" - |> bool_of_string - |> Ezfio.set_tc_scf_bi_ortho - ; - Ezfio.get_tc_scf_bi_ortho () - ;; -(* Write snippet for bi_ortho *) - let write_bi_ortho = - Ezfio.set_tc_scf_bi_ortho - ;; - -(* Read snippet for n_it_tcscf_max *) - let read_n_it_tcscf_max () = - if not (Ezfio.has_tc_scf_n_it_tcscf_max ()) then - get_default "n_it_tcscf_max" - |> int_of_string - |> Ezfio.set_tc_scf_n_it_tcscf_max - ; - Ezfio.get_tc_scf_n_it_tcscf_max () - |> Strictly_positive_int.of_int - ;; -(* Write snippet for n_it_tcscf_max *) - let write_n_it_tcscf_max var = - Strictly_positive_int.to_int var - |> Ezfio.set_tc_scf_n_it_tcscf_max - ;; - -(* Read snippet for thresh_tcscf *) - let read_thresh_tcscf () = - if not (Ezfio.has_tc_scf_thresh_tcscf ()) then - get_default "thresh_tcscf" - |> float_of_string - |> Ezfio.set_tc_scf_thresh_tcscf - ; - Ezfio.get_tc_scf_thresh_tcscf () - |> Threshold.of_float - ;; -(* Write snippet for thresh_tcscf *) - let write_thresh_tcscf var = - Threshold.to_float var - |> Ezfio.set_tc_scf_thresh_tcscf - ;; - -(* =~=~=~=~=~=~=~=~=~=~=~=~ *) -(* Generate Global Function *) -(* =~=~=~=~=~=~=~=~=~=~=~=~ *) - -(* Read all *) - let read() = - Some - { - bi_ortho = read_bi_ortho (); - thresh_tcscf = read_thresh_tcscf (); - n_it_tcscf_max = read_n_it_tcscf_max (); - } - ;; -(* Write all *) - let write{ - bi_ortho; - thresh_tcscf; - n_it_tcscf_max; - } = - write_bi_ortho bi_ortho; - write_thresh_tcscf thresh_tcscf; - write_n_it_tcscf_max n_it_tcscf_max; - ;; -(* to_string*) - let to_string b = - Printf.sprintf " - bi_ortho = %s - thresh_tcscf = %s - n_it_tcscf_max = %s - " - (string_of_bool b.bi_ortho) - (Threshold.to_string b.thresh_tcscf) - (Strictly_positive_int.to_string b.n_it_tcscf_max) - ;; -(* to_rst*) - let to_rst b = - Printf.sprintf " - If |true|, the MO basis is assumed to be bi-orthonormal :: - - bi_ortho = %s - - Threshold on the convergence of the Hartree Fock energy. :: - - thresh_tcscf = %s - - Maximum number of SCF iterations :: - - n_it_tcscf_max = %s - - " - (string_of_bool b.bi_ortho) - (Threshold.to_string b.thresh_tcscf) - (Strictly_positive_int.to_string b.n_it_tcscf_max) - |> Rst_string.of_string - ;; - include Generic_input_of_rst;; - let of_rst = of_rst t_of_sexp;; - -end \ No newline at end of file diff --git a/ocaml/Makefile b/ocaml/Makefile index 40d292fe..8853a991 100644 --- a/ocaml/Makefile +++ b/ocaml/Makefile @@ -29,7 +29,7 @@ tests: $(ALL_TESTS) .gitignore: $(MLFILES) $(MLIFILES) @for i in .gitignore ezfio.ml element_create_db Qptypes.ml Git.ml *.byte *.native _build $(ALL_EXE) $(ALL_TESTS) \ $(patsubst %.ml,%,$(wildcard test_*.ml)) $(patsubst %.ml,%,$(wildcard qp_*.ml)) \ - $(shell grep Input Input_auto_generated.ml | awk '{print $$2 ".ml"}') \ + Input_*.ml \ qp_edit.ml qp_edit qp_edit.native Input_auto_generated.ml;\ do \ echo $$i ; \ diff --git a/ocaml/Molecule.ml b/ocaml/Molecule.ml index 9b01ac3a..603244c8 100644 --- a/ocaml/Molecule.ml +++ b/ocaml/Molecule.ml @@ -101,7 +101,7 @@ let to_string_general ~f m = |> String.concat "\n" let to_string = - to_string_general ~f:(fun x -> Atom.to_string Units.Angstrom x) + to_string_general ~f:(fun x -> Atom.to_string ~units:Units.Angstrom x) let to_xyz = to_string_general ~f:Atom.to_xyz @@ -113,7 +113,7 @@ let of_xyz_string s = let l = String_ext.split s ~on:'\n' |> List.filter (fun x -> x <> "") - |> list_map (fun x -> Atom.of_string units x) + |> list_map (fun x -> Atom.of_string ~units x) in let ne = ( get_charge { nuclei=l ; diff --git a/ocaml/Progress_bar.ml b/ocaml/Progress_bar.ml index bc720b95..2cd3d19e 100644 --- a/ocaml/Progress_bar.ml +++ b/ocaml/Progress_bar.ml @@ -10,7 +10,7 @@ type t = next : float; } -let init ?(bar_length=20) ?(start_value=0.) ?(end_value=1.) ~title = +let init ?(bar_length=20) ?(start_value=0.) ?(end_value=1.) title = { title ; start_value ; end_value ; bar_length ; cur_value=start_value ; init_time= Unix.time () ; dirty = false ; next = Unix.time () } diff --git a/ocaml/Qputils.ml b/ocaml/Qputils.ml index 270e069f..752a65a0 100644 --- a/ocaml/Qputils.ml +++ b/ocaml/Qputils.ml @@ -56,3 +56,7 @@ let string_of_string s = s let list_map f l = List.rev_map f l |> List.rev + +let socket_convert socket = + ((Obj.magic (Obj.repr socket)) : [ `Xsub ] Zmq.Socket.t ) + diff --git a/ocaml/TaskServer.ml b/ocaml/TaskServer.ml index ad827316..5d9d2416 100644 --- a/ocaml/TaskServer.ml +++ b/ocaml/TaskServer.ml @@ -155,7 +155,7 @@ let new_job msg program_state rep_socket pair_socket = ~start_value:0. ~end_value:1. ~bar_length:20 - ~title:(Message.State.to_string state) + (Message.State.to_string state) in let result = diff --git a/ocaml/qp_create_ezfio.ml b/ocaml/qp_create_ezfio.ml index d6c8d66c..4583b118 100644 --- a/ocaml/qp_create_ezfio.ml +++ b/ocaml/qp_create_ezfio.ml @@ -677,12 +677,15 @@ let run ?o b au c d m p cart xyz_file = let () = + try ( let open Command_line in begin "Creates an EZFIO directory from a standard xyz file or from a z-matrix file in Gaussian format. The basis set is defined as a single string if all the atoms are taken from the same basis set, otherwise specific elements can be defined as follows: + -b \"cc-pcvdz | H:cc-pvdz | C:6-31g\" -b \"cc-pvtz | 1,H:sto-3g | 3,H:6-31g\" + If a file with the same name as the basis set exists, this file will be read. Otherwise, the basis set is obtained from the database. " |> set_description_doc ; set_header_doc (Sys.argv.(0) ^ " - Quantum Package command"); @@ -732,7 +735,7 @@ If a file with the same name as the basis set exists, this file will be read. O let basis = match Command_line.get "basis" with - | None -> assert false + | None -> "" | Some x -> x in @@ -771,10 +774,14 @@ If a file with the same name as the basis set exists, this file will be read. O let xyz_filename = match Command_line.anon_args () with - | [x] -> x - | _ -> (Command_line.help () ; failwith "input file is missing") + | [] -> failwith "input file is missing" + | x::_ -> x in run ?o:output basis au charge dummy multiplicity pseudo cart xyz_filename + ) + with + | Failure txt -> Printf.eprintf "Fatal error: %s\n%!" txt + | Command_line.Error txt -> Printf.eprintf "Command line error: %s\n%!" txt diff --git a/ocaml/qp_tunnel.ml b/ocaml/qp_tunnel.ml index 84e50eb5..6885db73 100644 --- a/ocaml/qp_tunnel.ml +++ b/ocaml/qp_tunnel.ml @@ -2,7 +2,7 @@ open Qputils open Qptypes type ezfio_or_address = EZFIO of string | ADDRESS of string -type req_or_sub = REQ | SUB +type req_or_sub = REQ | SUB let localport = 42379 @@ -29,7 +29,7 @@ let () = end; let arg = - let x = + let x = match Command_line.anon_args () with | [x] -> x | _ -> begin @@ -44,7 +44,7 @@ let () = in - let localhost = + let localhost = Lazy.force TaskServer.ip_address in @@ -52,28 +52,28 @@ let () = let long_address = match arg with | ADDRESS x -> x - | EZFIO x -> - let ic = + | EZFIO x -> + let ic = Filename.concat (Qpackage.ezfio_work x) "qp_run_address" |> open_in in - let result = + let result = input_line ic |> String.trim in close_in ic; result in - + let protocol, address, port = match String.split_on_char ':' long_address with | t :: a :: p :: [] -> t, a, int_of_string p - | _ -> failwith @@ + | _ -> failwith @@ Printf.sprintf "%s : Malformed address" long_address in - let zmq_context = + let zmq_context = Zmq.Context.create () in @@ -105,10 +105,10 @@ let () = let create_socket sock_type bind_or_connect addr = - let socket = + let socket = Zmq.Socket.create zmq_context sock_type in - let () = + let () = try bind_or_connect socket addr with @@ -131,37 +131,64 @@ let () = Sys.set_signal Sys.sigint handler; - let new_thread req_or_sub addr_in addr_out = + let new_thread_req addr_in addr_out = let socket_in, socket_out = - match req_or_sub with - | REQ -> create_socket Zmq.Socket.router Zmq.Socket.bind addr_in, create_socket Zmq.Socket.dealer Zmq.Socket.connect addr_out - | SUB -> - create_socket Zmq.Socket.sub Zmq.Socket.connect addr_in, - create_socket Zmq.Socket.pub Zmq.Socket.bind addr_out in - if req_or_sub = SUB then - Zmq.Socket.subscribe socket_in ""; - - - let action_in = - match req_or_sub with - | REQ -> (fun () -> Zmq.Socket.recv_all socket_in |> Zmq.Socket.send_all socket_out) - | SUB -> (fun () -> Zmq.Socket.recv_all socket_in |> Zmq.Socket.send_all socket_out) + let action_in = + fun () -> Zmq.Socket.recv_all socket_in |> Zmq.Socket.send_all socket_out in - let action_out = - match req_or_sub with - | REQ -> (fun () -> Zmq.Socket.recv_all socket_out |> Zmq.Socket.send_all socket_in ) - | SUB -> (fun () -> () ) + let action_out = + fun () -> Zmq.Socket.recv_all socket_out |> Zmq.Socket.send_all socket_in in let pollitem = Zmq.Poll.mask_of - [| (socket_in, Zmq.Poll.In) ; (socket_out, Zmq.Poll.In) |] + [| (socket_convert socket_in, Zmq.Poll.In) ; (socket_convert socket_out, Zmq.Poll.In) |] + in + + while !run_status do + + let polling = + Zmq.Poll.poll ~timeout:1000 pollitem + in + + match polling with + | [| Some Zmq.Poll.In ; Some Zmq.Poll.In |] -> ( action_out () ; action_in () ) + | [| _ ; Some Zmq.Poll.In |] -> action_out () + | [| Some Zmq.Poll.In ; _ |] -> action_in () + | _ -> () + done; + + Zmq.Socket.close socket_in; + Zmq.Socket.close socket_out; + in + + let new_thread_sub addr_in addr_out = + let socket_in, socket_out = + create_socket Zmq.Socket.sub Zmq.Socket.connect addr_in, + create_socket Zmq.Socket.pub Zmq.Socket.bind addr_out + in + + Zmq.Socket.subscribe socket_in ""; + + + + let action_in = + fun () -> Zmq.Socket.recv_all socket_in |> Zmq.Socket.send_all socket_out + in + + let action_out = + fun () -> () + in + + let pollitem = + Zmq.Poll.mask_of + [| (socket_convert socket_in, Zmq.Poll.In) ; (socket_convert socket_out, Zmq.Poll.In) |] in @@ -173,8 +200,8 @@ let () = match polling with | [| Some Zmq.Poll.In ; Some Zmq.Poll.In |] -> ( action_out () ; action_in () ) - | [| _ ; Some Zmq.Poll.In |] -> action_out () - | [| Some Zmq.Poll.In ; _ |] -> action_in () + | [| _ ; Some Zmq.Poll.In |] -> action_out () + | [| Some Zmq.Poll.In ; _ |] -> action_in () | _ -> () done; @@ -193,8 +220,8 @@ let () = Printf.sprintf "tcp://*:%d" localport in - let f () = - new_thread REQ addr_in addr_out + let f () = + new_thread_req addr_in addr_out in (Thread.create f) () @@ -211,8 +238,8 @@ let () = Printf.sprintf "tcp://*:%d" (localport+2) in - let f () = - new_thread REQ addr_in addr_out + let f () = + new_thread_req addr_in addr_out in (Thread.create f) () in @@ -227,8 +254,8 @@ let () = Printf.sprintf "tcp://*:%d" (localport+1) in - let f () = - new_thread SUB addr_in addr_out + let f () = + new_thread_sub addr_in addr_out in (Thread.create f) () in @@ -236,7 +263,7 @@ let () = let input_thread = - let f () = + let f () = let addr_out = match arg with | EZFIO _ -> None @@ -248,22 +275,22 @@ let () = Printf.sprintf "tcp://*:%d" (localport+9) in - let socket_in = + let socket_in = create_socket Zmq.Socket.rep Zmq.Socket.bind addr_in in let socket_out = - match addr_out with + match addr_out with | Some addr_out -> Some ( create_socket Zmq.Socket.req Zmq.Socket.connect addr_out) | None -> None in - let temp_file = + let temp_file = Filename.temp_file "qp_tunnel" ".tar.gz" in - let get_ezfio_filename () = + let get_ezfio_filename () = match arg with | EZFIO x -> x | ADDRESS _ -> @@ -277,9 +304,9 @@ let () = end in - let get_input () = + let get_input () = match arg with - | EZFIO x -> + | EZFIO x -> begin Printf.sprintf "tar --exclude=\"*.gz.*\" -zcf %s %s" temp_file x |> Sys.command |> ignore; @@ -291,11 +318,11 @@ let () = in ignore @@ Unix.lseek fd 0 Unix.SEEK_SET ; let bstr = - Unix.map_file fd Bigarray.char + Unix.map_file fd Bigarray.char Bigarray.c_layout false [| len |] |> Bigarray.array1_of_genarray in - let result = + let result = String.init len (fun i -> bstr.{i}) ; in Unix.close fd; @@ -313,7 +340,7 @@ let () = end in - let () = + let () = match socket_out with | None -> () | Some socket_out -> @@ -329,7 +356,7 @@ let () = | ADDRESS _ -> begin Printf.printf "Getting input... %!"; - let ezfio_filename = + let ezfio_filename = get_ezfio_filename () in Printf.printf "%s%!" ezfio_filename; @@ -343,7 +370,7 @@ let () = |> Sys.command |> ignore ; let oc = Filename.concat (Qpackage.ezfio_work ezfio_filename) "qp_run_address" - |> open_out + |> open_out in Printf.fprintf oc "tcp://%s:%d\n" localhost localport; close_out oc; @@ -359,9 +386,9 @@ let () = let action () = match Zmq.Socket.recv socket_in with | "get_input" -> get_input () - |> Zmq.Socket.send socket_in + |> Zmq.Socket.send socket_in | "get_ezfio_filename" -> get_ezfio_filename () - |> Zmq.Socket.send socket_in + |> Zmq.Socket.send socket_in | "test" -> Zmq.Socket.send socket_in "OK" | x -> Printf.sprintf "Message '%s' not understood" x |> Zmq.Socket.send socket_in @@ -372,7 +399,7 @@ On remote hosts, create ssh tunnel using: ssh -L %d:%s:%d -L %d:%s:%d -L %d:%s:%d -L %d:%s:%d %s & Or from this host connect to clients using: ssh -R %d:localhost:%d -R %d:localhost:%d -R %d:localhost:%d -R %d:localhost:%d & -%!" +%!" (port ) localhost (localport ) (port+1) localhost (localport+1) (port+2) localhost (localport+2) @@ -392,12 +419,12 @@ Or from this host connect to clients using: match polling.(0) with | Some Zmq.Poll.In -> action () | None -> () - | Some Zmq.Poll.In_out + | Some Zmq.Poll.In_out | Some Zmq.Poll.Out -> () done; - let () = + let () = match socket_out with | Some socket_out -> Zmq.Socket.close socket_out | None -> () @@ -415,7 +442,7 @@ Or from this host connect to clients using: Thread.join ocaml_thread; Zmq.Context.terminate zmq_context; Printf.printf "qp_tunnel exited properly.\n" - + diff --git a/scripts/.gitignore b/scripts/.gitignore index b44ac5a2..103b3ae9 100644 --- a/scripts/.gitignore +++ b/scripts/.gitignore @@ -2,3 +2,6 @@ *.pyo docopt.py resultsFile/ +verif_omp/a.out +src/*/Makefile +src/*/*/ diff --git a/scripts/compilation/qp_create_ninja b/scripts/compilation/qp_create_ninja index 0f70f4c4..7df3c62d 100755 --- a/scripts/compilation/qp_create_ninja +++ b/scripts/compilation/qp_create_ninja @@ -99,9 +99,20 @@ def ninja_create_env_variable(pwd_config_file): l_string = ["builddir = {0}".format(os.path.dirname(ROOT_BUILD_NINJA)), ""] + for flag in ["FC", "FCFLAGS", "IRPF90", "IRPF90_FLAGS"]: str_ = "{0} = {1}".format(flag, get_compilation_option(pwd_config_file, flag)) + for directory in [real_join(QP_SRC, m) for m in sorted(os.listdir(QP_SRC))]: + includefile = real_join(directory, flag) + try: + content = "" + with open(includefile,'r') as f: + content = f.read() + str_ += " "+content + except IOError: + pass + l_string.append(str_) lib_lapack = get_compilation_option(pwd_config_file, "LAPACK_LIB") @@ -110,17 +121,19 @@ def ninja_create_env_variable(pwd_config_file): str_lib = " ".join([lib_lapack, EZFIO_LIB, ZMQ_LIB, LIB, lib_usr]) # Read all LIB files in modules - libfile = "LIB" - try: - content = "" - with open(libfile,'r') as f: - content = f.read() - str_lib += " "+content - except IOError: - pass + for directory in [real_join(QP_SRC, m) for m in sorted(os.listdir(QP_SRC))]: + libfile = real_join(directory, "LIB") + try: + content = "" + with open(libfile,'r') as f: + content = f.read() + str_lib += " "+content + except IOError: + pass l_string.append("LIB = {0} ".format(str_lib)) + l_string.append("CONFIG_FILE = {0}".format(pwd_config_file)) l_string.append("") diff --git a/src/.gitignore b/src/.gitignore new file mode 100644 index 00000000..6353c21a --- /dev/null +++ b/src/.gitignore @@ -0,0 +1,11 @@ +* +!README.rst +!*/ +*/* +!*/*.* +*/*.o +*/build.ninja +*/ezfio_interface.irp.f +*/.gitignore +*/*.swp + diff --git a/src/ao_basis/EZFIO.cfg b/src/ao_basis/EZFIO.cfg index e1a7d268..dd61b1be 100644 --- a/src/ao_basis/EZFIO.cfg +++ b/src/ao_basis/EZFIO.cfg @@ -67,4 +67,3 @@ doc: Use normalized primitive functions interface: ezfio, provider default: true - diff --git a/src/ao_basis/spherical_to_cartesian.irp.f b/src/ao_basis/spherical_to_cartesian.irp.f index 33a3bc89..336161f8 100644 --- a/src/ao_basis/spherical_to_cartesian.irp.f +++ b/src/ao_basis/spherical_to_cartesian.irp.f @@ -1,7 +1,7 @@ ! Spherical to cartesian transformation matrix obtained with ! Horton (http://theochem.github.com/horton/, 2015) -! First index is the index of the carteisan AO, obtained by ao_power_index +! First index is the index of the cartesian AO, obtained by ao_power_index ! Second index is the index of the spherical AO BEGIN_PROVIDER [ double precision, cart_to_sphe_0, (1,1) ] diff --git a/src/ao_many_one_e_ints/grad2_jmu_manu.irp.f b/src/ao_many_one_e_ints/grad2_jmu_manu.irp.f new file mode 100644 index 00000000..c25d8055 --- /dev/null +++ b/src/ao_many_one_e_ints/grad2_jmu_manu.irp.f @@ -0,0 +1,284 @@ + +BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test, (ao_num, ao_num, n_points_final_grid)] + + BEGIN_DOC + ! + ! int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2)^2 u_12^mu [\grad_1 u_12^mu] + ! + END_DOC + + implicit none + integer :: i, j, ipoint, i_1s, i_fit + double precision :: r(3), int_fit, expo_fit, coef_fit, coef_tmp + double precision :: coef, beta, B_center(3), dist + double precision :: alpha_1s, alpha_1s_inv, centr_1s(3), expo_coef_1s, tmp + double precision :: wall0, wall1 + double precision, external :: NAI_pol_mult_erf_ao_with1s + double precision :: j12_mu_r12,int_j1b + double precision :: sigma_ij,dist_ij_ipoint,dsqpi_3_2 + double precision :: beta_ij,center_ij_1s(3),factor_ij_1s + dsqpi_3_2 = (dacos(-1.d0))**(3/2) + + provide mu_erf final_grid_points j1b_pen ao_overlap_abs List_comb_thr_b3_cent + call wall_time(wall0) + + + int2_u_grad1u_j1b2_test = 0.d0 + + !$OMP PARALLEL DEFAULT (NONE) & + !$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, & + !$OMP coef_fit, expo_fit, int_fit, tmp, alpha_1s, dist, & + !$OMP beta_ij,center_ij_1s,factor_ij_1s, & + !$OMP int_j1b,alpha_1s_inv, centr_1s, expo_coef_1s, coef_tmp) & + !$OMP SHARED (n_points_final_grid, ao_num, List_comb_b3_size_thr, & + !$OMP final_grid_points, ng_fit_jast, & + !$OMP expo_gauss_j_mu_1_erf, coef_gauss_j_mu_1_erf, & + !$OMP ao_prod_dist_grid, ao_prod_sigma, ao_overlap_abs_grid,ao_prod_center,dsqpi_3_2, & + !$OMP List_comb_thr_b3_coef, List_comb_thr_b3_expo, ao_abs_comb_b3_j1b, & + !$OMP List_comb_thr_b3_cent, int2_u_grad1u_j1b2_test) + !$OMP DO + do ipoint = 1, n_points_final_grid + do i = 1, ao_num + do j = i, ao_num + if(dabs(ao_overlap_abs_grid(j,i)).lt.1.d-10)cycle + r(1) = final_grid_points(1,ipoint) + r(2) = final_grid_points(2,ipoint) + r(3) = final_grid_points(3,ipoint) + + tmp = 0.d0 + do i_1s = 1, List_comb_b3_size_thr(j,i) + + coef = List_comb_thr_b3_coef (i_1s,j,i) + beta = List_comb_thr_b3_expo (i_1s,j,i) + int_j1b = ao_abs_comb_b3_j1b(i_1s,j,i) + if(dabs(coef)*dabs(int_j1b).lt.1.d-10)cycle + B_center(1) = List_comb_thr_b3_cent(1,i_1s,j,i) + B_center(2) = List_comb_thr_b3_cent(2,i_1s,j,i) + B_center(3) = List_comb_thr_b3_cent(3,i_1s,j,i) + dist = (B_center(1) - r(1)) * (B_center(1) - r(1)) & + + (B_center(2) - r(2)) * (B_center(2) - r(2)) & + + (B_center(3) - r(3)) * (B_center(3) - r(3)) + + do i_fit = 1, ng_fit_jast + + expo_fit = expo_gauss_j_mu_1_erf(i_fit) + call gaussian_product(expo_fit,r,beta,B_center,factor_ij_1s,beta_ij,center_ij_1s) +! if(factor_ij_1s*dabs(coef*int_j1b)*dsqpi_3_2*beta_ij**(-3/2).lt.1.d-15)cycle + coef_fit = coef_gauss_j_mu_1_erf(i_fit) + + alpha_1s = beta + expo_fit + alpha_1s_inv = 1.d0 / alpha_1s + centr_1s(1) = alpha_1s_inv * (beta * B_center(1) + expo_fit * r(1)) + centr_1s(2) = alpha_1s_inv * (beta * B_center(2) + expo_fit * r(2)) + centr_1s(3) = alpha_1s_inv * (beta * B_center(3) + expo_fit * r(3)) + + expo_coef_1s = beta * expo_fit * alpha_1s_inv * dist + if(expo_coef_1s .gt. 20.d0) cycle + coef_tmp = coef * coef_fit * dexp(-expo_coef_1s) + if(dabs(coef_tmp) .lt. 1d-08) cycle + + int_fit = NAI_pol_mult_erf_ao_with1s(i, j, alpha_1s, centr_1s, 1.d+9, r) + + tmp += coef_tmp * int_fit + enddo + enddo + + int2_u_grad1u_j1b2_test(j,i,ipoint) = tmp + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + do ipoint = 1, n_points_final_grid + do i = 2, ao_num + do j = 1, i-1 + int2_u_grad1u_j1b2_test(j,i,ipoint) = int2_u_grad1u_j1b2_test(i,j,ipoint) + enddo + enddo + enddo + + call wall_time(wall1) + print*, ' wall time for int2_u_grad1u_j1b2_test', wall1 - wall0 + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test_no_v, (ao_num, ao_num, n_points_final_grid)] + + BEGIN_DOC + ! + ! -\frac{1}{4} x int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2)^2 [1 - erf(mu r12)]^2 + ! + END_DOC + + implicit none + integer :: i, j, ipoint, i_1s, i_fit + double precision :: r(3), expo_fit, coef_fit + double precision :: coef, beta, B_center(3) + double precision :: tmp + double precision :: wall0, wall1 + + double precision, allocatable :: int_fit_v(:) + double precision, external :: overlap_gauss_r12_ao_with1s + double precision :: factor_ij_1s,beta_ij,center_ij_1s(3),int_j1b,int_gauss,dsqpi_3_2 + dsqpi_3_2 = (dacos(-1.d0))**(3/2) + + provide mu_erf final_grid_points_transp j1b_pen List_comb_thr_b3_coef + call wall_time(wall0) + + int2_grad1u2_grad2u2_j1b2_test_no_v(:,:,:) = 0.d0 + + !$OMP PARALLEL DEFAULT (NONE) & + !$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center,& + !$OMP coef_fit, expo_fit, int_fit_v, tmp,int_gauss,int_j1b,factor_ij_1s,beta_ij,center_ij_1s) & + !$OMP SHARED (n_points_final_grid, ao_num, final_grid_points,List_comb_b3_size_thr,& + !$OMP final_grid_points_transp, ng_fit_jast, & + !$OMP expo_gauss_1_erf_x_2, coef_gauss_1_erf_x_2, & + !$OMP List_comb_thr_b3_coef, List_comb_thr_b3_expo, & + !$OMP List_comb_thr_b3_cent, int2_grad1u2_grad2u2_j1b2_test_no_v, ao_abs_comb_b3_j1b,& + !$OMP ao_overlap_abs,dsqpi_3_2) + !$OMP DO SCHEDULE(dynamic) + do ipoint = 1, n_points_final_grid + r(1) = final_grid_points(1,ipoint) + r(2) = final_grid_points(2,ipoint) + r(3) = final_grid_points(3,ipoint) + do i = 1, ao_num + do j = i, ao_num + if(ao_overlap_abs(j,i) .lt. 1.d-12) then + cycle + endif + + do i_1s = 1, List_comb_b3_size_thr(j,i) + + coef = List_comb_thr_b3_coef (i_1s,j,i) + beta = List_comb_thr_b3_expo (i_1s,j,i) + int_j1b = ao_abs_comb_b3_j1b(i_1s,j,i) +! if(dabs(coef)*dabs(int_j1b).lt.1.d-15)cycle + B_center(1) = List_comb_thr_b3_cent(1,i_1s,j,i) + B_center(2) = List_comb_thr_b3_cent(2,i_1s,j,i) + B_center(3) = List_comb_thr_b3_cent(3,i_1s,j,i) + + do i_fit = 1, ng_fit_jast + + expo_fit = expo_gauss_1_erf_x_2(i_fit) +! call gaussian_product(expo_fit,r,beta,B_center,factor_ij_1s,beta_ij,center_ij_1s) + coef_fit = -0.25d0 * coef_gauss_1_erf_x_2(i_fit) * coef +! if(dabs(coef_fit)*factor_ij_1s*dabs(int_j1b).lt.1.d-15)cycle + +! call overlap_gauss_r12_ao_with1s_v(B_center, beta, final_grid_points_transp, & +! expo_fit, i, j, int_fit_v, n_points_final_grid) + int_gauss = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j) + + int2_grad1u2_grad2u2_j1b2_test_no_v(j,i,ipoint) += coef_fit * int_gauss + + enddo + enddo + enddo + enddo + enddo + + !$OMP END DO + !$OMP END PARALLEL + + do ipoint = 1, n_points_final_grid + do i = 1, ao_num + do j = 1, i-1 + int2_grad1u2_grad2u2_j1b2_test_no_v(j,i,ipoint) = int2_grad1u2_grad2u2_j1b2_test_no_v(i,j,ipoint) + enddo + enddo + enddo + + call wall_time(wall1) + print*, ' wall time for int2_grad1u2_grad2u2_j1b2_test_no_v', wall1 - wall0 + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test, (ao_num, ao_num, n_points_final_grid)] +! +! BEGIN_DOC +! ! +! ! -\frac{1}{4} x int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2)^2 [1 - erf(mu r12)]^2 +! ! +! END_DOC +! + implicit none + integer :: i, j, ipoint, i_1s, i_fit + double precision :: r(3), expo_fit, coef_fit + double precision :: coef, beta, B_center(3) + double precision :: tmp + double precision :: wall0, wall1 + + double precision, allocatable :: int_fit_v(:) + double precision, external :: overlap_gauss_r12_ao_with1s + + provide mu_erf final_grid_points_transp j1b_pen + call wall_time(wall0) + + double precision :: int_j1b + int2_grad1u2_grad2u2_j1b2_test(:,:,:) = 0.d0 +! + !$OMP PARALLEL DEFAULT (NONE) & + !$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center,& + !$OMP coef_fit, expo_fit, int_fit_v, tmp,int_j1b) & + !$OMP SHARED (n_points_final_grid, ao_num, List_comb_b3_size_thr,& + !$OMP final_grid_points_transp, ng_fit_jast, & + !$OMP expo_gauss_1_erf_x_2, coef_gauss_1_erf_x_2, & + !$OMP List_comb_thr_b3_coef, List_comb_thr_b3_expo, & + !$OMP List_comb_thr_b3_cent, int2_grad1u2_grad2u2_j1b2_test,& + !$OMP ao_abs_comb_b3_j1b,ao_overlap_abs) +! + allocate(int_fit_v(n_points_final_grid)) + !$OMP DO SCHEDULE(dynamic) + do i = 1, ao_num + do j = i, ao_num + + if(ao_overlap_abs(j,i) .lt. 1.d-12) then + cycle + endif + + do i_1s = 1, List_comb_b3_size_thr(j,i) + + coef = List_comb_thr_b3_coef (i_1s,j,i) + beta = List_comb_thr_b3_expo (i_1s,j,i) + int_j1b = ao_abs_comb_b3_j1b(i_1s,j,i) +! if(dabs(coef)*dabs(int_j1b).lt.1.d-15)cycle + B_center(1) = List_comb_thr_b3_cent(1,i_1s,j,i) + B_center(2) = List_comb_thr_b3_cent(2,i_1s,j,i) + B_center(3) = List_comb_thr_b3_cent(3,i_1s,j,i) + + do i_fit = 1, ng_fit_jast + + expo_fit = expo_gauss_1_erf_x_2(i_fit) + coef_fit = -0.25d0 * coef_gauss_1_erf_x_2(i_fit) * coef + + call overlap_gauss_r12_ao_with1s_v(B_center, beta, final_grid_points_transp, size(final_grid_points_transp,1),& + expo_fit, i, j, int_fit_v, size(int_fit_v,1),n_points_final_grid) + + do ipoint = 1, n_points_final_grid + int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint) += coef_fit * int_fit_v(ipoint) + enddo + + enddo + + enddo + enddo + enddo + !$OMP END DO + deallocate(int_fit_v) + !$OMP END PARALLEL + + do ipoint = 1, n_points_final_grid + do i = 2, ao_num + do j = 1, i-1 + int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint) = int2_grad1u2_grad2u2_j1b2_test(i,j,ipoint) + enddo + enddo + enddo + + call wall_time(wall1) + print*, ' wall time for int2_grad1u2_grad2u2_j1b2_test', wall1 - wall0 + +END_PROVIDER + diff --git a/src/ao_many_one_e_ints/grad2_jmu_modif.irp.f b/src/ao_many_one_e_ints/grad2_jmu_modif.irp.f index b7fe234f..872bfaef 100644 --- a/src/ao_many_one_e_ints/grad2_jmu_modif.irp.f +++ b/src/ao_many_one_e_ints/grad2_jmu_modif.irp.f @@ -351,7 +351,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2, (ao_num, ao_num, n_points ! --- int_fit = NAI_pol_mult_erf_ao_with1s(i, j, expo_fit, r, 1.d+9, r) - if(dabs(int_fit) .lt. 1d-10) cycle +! if(dabs(int_fit) .lt. 1d-10) cycle tmp += coef_fit * int_fit @@ -375,9 +375,10 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2, (ao_num, ao_num, n_points centr_1s(3) = alpha_1s_inv * (beta * B_center(3) + expo_fit * r(3)) expo_coef_1s = beta * expo_fit * alpha_1s_inv * dist +! if(expo_coef_1s .gt. 80.d0) cycle coef_tmp = coef * coef_fit * dexp(-expo_coef_1s) - if(dabs(coef_tmp) .lt. 1d-10) cycle - +! if(dabs(coef_tmp) .lt. 1d-10) cycle + int_fit = NAI_pol_mult_erf_ao_with1s(i, j, alpha_1s, centr_1s, 1.d+9, r) tmp += coef_tmp * int_fit diff --git a/src/ao_many_one_e_ints/grad_lapl_jmu_manu.irp.f b/src/ao_many_one_e_ints/grad_lapl_jmu_manu.irp.f new file mode 100644 index 00000000..1b457d68 --- /dev/null +++ b/src/ao_many_one_e_ints/grad_lapl_jmu_manu.irp.f @@ -0,0 +1,287 @@ + +! --- + +BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_num, n_points_final_grid)] + + BEGIN_DOC + ! + ! int dr phi_i(r) phi_j(r) 1s_j1b(r) (erf(mu(R) |r - R| - 1) / |r - R| + ! + END_DOC + + implicit none + integer :: i, j, ipoint, i_1s + double precision :: r(3), int_mu, int_coulomb + double precision :: coef, beta, B_center(3) + double precision :: tmp,int_j1b + double precision :: wall0, wall1 + double precision, external :: NAI_pol_mult_erf_ao_with1s + double precision :: sigma_ij,dist_ij_ipoint,dsqpi_3_2 + dsqpi_3_2 = (dacos(-1.d0))**(3/2) + provide mu_erf final_grid_points j1b_pen + call wall_time(wall0) + + v_ij_erf_rk_cst_mu_j1b_test = 0.d0 + + !$OMP PARALLEL DEFAULT (NONE) & + !$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, int_mu, int_coulomb, tmp, int_j1b)& + !$OMP SHARED (n_points_final_grid, ao_num, List_comb_b2_size_thr, final_grid_points, & + !$OMP List_comb_thr_b2_coef, List_comb_thr_b2_expo, List_comb_thr_b2_cent,ao_abs_comb_b2_j1b, & + !$OMP v_ij_erf_rk_cst_mu_j1b_test, mu_erf, & + !$OMP ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2) + !$OMP DO + !do ipoint = 1, 10 + do ipoint = 1, n_points_final_grid + r(1) = final_grid_points(1,ipoint) + r(2) = final_grid_points(2,ipoint) + r(3) = final_grid_points(3,ipoint) + + do i = 1, ao_num + do j = i, ao_num + if(dabs(ao_overlap_abs_grid(j,i)).lt.1.d-20)cycle + + tmp = 0.d0 + do i_1s = 1, List_comb_b2_size_thr(j,i) + + coef = List_comb_thr_b2_coef (i_1s,j,i) + beta = List_comb_thr_b2_expo (i_1s,j,i) + int_j1b = ao_abs_comb_b2_j1b(i_1s,j,i) + if(dabs(coef)*dabs(int_j1b).lt.1.d-10)cycle + B_center(1) = List_comb_thr_b2_cent(1,i_1s,j,i) + B_center(2) = List_comb_thr_b2_cent(2,i_1s,j,i) + B_center(3) = List_comb_thr_b2_cent(3,i_1s,j,i) + + int_mu = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r) + int_coulomb = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r) + + tmp += coef * (int_mu - int_coulomb) + enddo + + v_ij_erf_rk_cst_mu_j1b_test(j,i,ipoint) = tmp + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + do ipoint = 1, n_points_final_grid + do i = 2, ao_num + do j = 1, i-1 + v_ij_erf_rk_cst_mu_j1b_test(j,i,ipoint) = v_ij_erf_rk_cst_mu_j1b_test(i,j,ipoint) + enddo + enddo + enddo + + call wall_time(wall1) + print*, ' wall time for v_ij_erf_rk_cst_mu_j1b_test', wall1 - wall0 + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_num, n_points_final_grid, 3)] + + BEGIN_DOC + ! int dr x phi_i(r) phi_j(r) 1s_j1b(r) (erf(mu(R) |r - R|) - 1)/|r - R| + END_DOC + + implicit none + integer :: i, j, ipoint + double precision :: wall0, wall1 + + call wall_time(wall0) + + do ipoint = 1, n_points_final_grid + do i = 1, ao_num + do j = 1, ao_num + x_v_ij_erf_rk_cst_mu_j1b_test(j,i,ipoint,1) = x_v_ij_erf_rk_cst_mu_tmp_j1b(1,j,i,ipoint) + x_v_ij_erf_rk_cst_mu_j1b_test(j,i,ipoint,2) = x_v_ij_erf_rk_cst_mu_tmp_j1b(2,j,i,ipoint) + x_v_ij_erf_rk_cst_mu_j1b_test(j,i,ipoint,3) = x_v_ij_erf_rk_cst_mu_tmp_j1b(3,j,i,ipoint) + enddo + enddo + enddo + + call wall_time(wall1) + print*, ' wall time for x_v_ij_erf_rk_cst_mu_j1b_test', wall1 - wall0 + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_tmp_j1b_test, (3, ao_num, ao_num, n_points_final_grid)] + + BEGIN_DOC + ! int dr x phi_i(r) phi_j(r) 1s_j1b(r) (erf(mu(R) |r - R|) - 1)/|r - R| + END_DOC + + implicit none + integer :: i, j, ipoint, i_1s + double precision :: coef, beta, B_center(3), r(3), ints(3), ints_coulomb(3) + double precision :: tmp_x, tmp_y, tmp_z + double precision :: wall0, wall1 + double precision :: sigma_ij,dist_ij_ipoint,dsqpi_3_2,int_j1b + dsqpi_3_2 = (dacos(-1.d0))**(3/2) + + call wall_time(wall0) + + x_v_ij_erf_rk_cst_mu_tmp_j1b_test = 0.d0 + + !$OMP PARALLEL DEFAULT (NONE) & + !$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, ints, ints_coulomb, & + !$OMP int_j1b, tmp_x, tmp_y, tmp_z) & + !$OMP SHARED (n_points_final_grid, ao_num, List_comb_b2_size_thr, final_grid_points,& + !$OMP List_comb_thr_b2_coef, List_comb_thr_b2_expo, List_comb_thr_b2_cent, & + !$OMP x_v_ij_erf_rk_cst_mu_tmp_j1b_test, mu_erf,ao_abs_comb_b2_j1b, & + !$OMP ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2) + !$OMP DO + !do ipoint = 1, 10 + do ipoint = 1, n_points_final_grid + r(1) = final_grid_points(1,ipoint) + r(2) = final_grid_points(2,ipoint) + r(3) = final_grid_points(3,ipoint) + + do i = 1, ao_num + do j = i, ao_num + if(dabs(ao_overlap_abs_grid(j,i)).lt.1.d-20)cycle + + tmp_x = 0.d0 + tmp_y = 0.d0 + tmp_z = 0.d0 + do i_1s = 1, List_comb_b2_size_thr(j,i) + + coef = List_comb_thr_b2_coef (i_1s,j,i) + beta = List_comb_thr_b2_expo (i_1s,j,i) + int_j1b = ao_abs_comb_b2_j1b(i_1s,j,i) + if(dabs(coef)*dabs(int_j1b).lt.1.d-10)cycle + B_center(1) = List_comb_thr_b2_cent(1,i_1s,j,i) + B_center(2) = List_comb_thr_b2_cent(2,i_1s,j,i) + B_center(3) = List_comb_thr_b2_cent(3,i_1s,j,i) + + call NAI_pol_x_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r, ints ) + call NAI_pol_x_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r, ints_coulomb) + + tmp_x += coef * (ints(1) - ints_coulomb(1)) + tmp_y += coef * (ints(2) - ints_coulomb(2)) + tmp_z += coef * (ints(3) - ints_coulomb(3)) + enddo + + x_v_ij_erf_rk_cst_mu_tmp_j1b_test(1,j,i,ipoint) = tmp_x + x_v_ij_erf_rk_cst_mu_tmp_j1b_test(2,j,i,ipoint) = tmp_y + x_v_ij_erf_rk_cst_mu_tmp_j1b_test(3,j,i,ipoint) = tmp_z + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + do ipoint = 1, n_points_final_grid + do i = 2, ao_num + do j = 1, i-1 + x_v_ij_erf_rk_cst_mu_tmp_j1b_test(1,j,i,ipoint) = x_v_ij_erf_rk_cst_mu_tmp_j1b_test(1,i,j,ipoint) + x_v_ij_erf_rk_cst_mu_tmp_j1b_test(2,j,i,ipoint) = x_v_ij_erf_rk_cst_mu_tmp_j1b_test(2,i,j,ipoint) + x_v_ij_erf_rk_cst_mu_tmp_j1b_test(3,j,i,ipoint) = x_v_ij_erf_rk_cst_mu_tmp_j1b_test(3,i,j,ipoint) + enddo + enddo + enddo + + call wall_time(wall1) + print*, ' wall time for x_v_ij_erf_rk_cst_mu_tmp_j1b_test', wall1 - wall0 + +END_PROVIDER + +! --- + +! TODO analytically +BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_test, (ao_num, ao_num, n_points_final_grid)] + + BEGIN_DOC + ! + ! int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2) u(mu, r12) + ! + END_DOC + + implicit none + integer :: i, j, ipoint, i_1s, i_fit + double precision :: r(3), int_fit, expo_fit, coef_fit + double precision :: coef, beta, B_center(3) + double precision :: tmp + double precision :: wall0, wall1 + + double precision, external :: overlap_gauss_r12_ao_with1s + double precision :: sigma_ij,dist_ij_ipoint,dsqpi_3_2,int_j1b + dsqpi_3_2 = (dacos(-1.d0))**(3/2) + + provide mu_erf final_grid_points j1b_pen + call wall_time(wall0) + + v_ij_u_cst_mu_j1b_test = 0.d0 + + !$OMP PARALLEL DEFAULT (NONE) & + !$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, & + !$OMP beta_ij_u, factor_ij_1s_u, center_ij_1s_u, & + !$OMP coef_fit, expo_fit, int_fit, tmp,coeftot,int_j1b) & + !$OMP SHARED (n_points_final_grid, ao_num, & + !$OMP final_grid_points, ng_fit_jast, & + !$OMP expo_gauss_j_mu_x, coef_gauss_j_mu_x, & + !$OMP List_comb_thr_b2_coef, List_comb_thr_b2_expo,List_comb_b2_size_thr, & + !$OMP List_comb_thr_b2_cent, v_ij_u_cst_mu_j1b_test,ao_abs_comb_b2_j1b, & + !$OMP ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2) + !$OMP DO + !do ipoint = 1, 10 + do ipoint = 1, n_points_final_grid + r(1) = final_grid_points(1,ipoint) + r(2) = final_grid_points(2,ipoint) + r(3) = final_grid_points(3,ipoint) + + do i = 1, ao_num + do j = i, ao_num + if(dabs(ao_overlap_abs_grid(j,i)).lt.1.d-20)cycle + + tmp = 0.d0 + do i_1s = 1, List_comb_b2_size_thr(j,i) + + coef = List_comb_thr_b2_coef (i_1s,j,i) + beta = List_comb_thr_b2_expo (i_1s,j,i) + int_j1b = ao_abs_comb_b2_j1b(i_1s,j,i) + if(dabs(coef)*dabs(int_j1b).lt.1.d-10)cycle + B_center(1) = List_comb_thr_b2_cent(1,i_1s,j,i) + B_center(2) = List_comb_thr_b2_cent(2,i_1s,j,i) + B_center(3) = List_comb_thr_b2_cent(3,i_1s,j,i) + + do i_fit = 1, ng_fit_jast + + expo_fit = expo_gauss_j_mu_x(i_fit) + coef_fit = coef_gauss_j_mu_x(i_fit) + coeftot = coef * coef_fit + if(dabs(coeftot).lt.1.d-15)cycle + double precision :: beta_ij_u, factor_ij_1s_u, center_ij_1s_u(3),coeftot + call gaussian_product(beta,B_center,expo_fit,r,factor_ij_1s_u,beta_ij_u,center_ij_1s_u) + if(factor_ij_1s_u*ao_overlap_abs_grid(j,i).lt.1.d-15)cycle + int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j) + + tmp += coef * coef_fit * int_fit + enddo + enddo + + v_ij_u_cst_mu_j1b_test(j,i,ipoint) = tmp + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + do ipoint = 1, n_points_final_grid + do i = 2, ao_num + do j = 1, i-1 + v_ij_u_cst_mu_j1b_test(j,i,ipoint) = v_ij_u_cst_mu_j1b_test(i,j,ipoint) + enddo + enddo + enddo + + call wall_time(wall1) + print*, ' wall time for v_ij_u_cst_mu_j1b_test', wall1 - wall0 + +END_PROVIDER + +! --- + diff --git a/src/ao_many_one_e_ints/listj1b.irp.f b/src/ao_many_one_e_ints/listj1b.irp.f index 0b40170c..c41b312d 100644 --- a/src/ao_many_one_e_ints/listj1b.irp.f +++ b/src/ao_many_one_e_ints/listj1b.irp.f @@ -168,7 +168,6 @@ END_PROVIDER do j = 1, nucl_num tmp_alphaj = dble(List_all_comb_b3(j,i)) * j1b_pen(j) - !print*, List_all_comb_b3(j,i), j1b_pen(j) List_all_comb_b3_expo(i) += tmp_alphaj List_all_comb_b3_cent(1,i) += tmp_alphaj * nucl_coord(j,1) List_all_comb_b3_cent(2,i) += tmp_alphaj * nucl_coord(j,2) diff --git a/src/ao_many_one_e_ints/listj1b_sorted.irp.f b/src/ao_many_one_e_ints/listj1b_sorted.irp.f new file mode 100644 index 00000000..606664f8 --- /dev/null +++ b/src/ao_many_one_e_ints/listj1b_sorted.irp.f @@ -0,0 +1,191 @@ + + BEGIN_PROVIDER [ integer, List_comb_b2_size_thr, (ao_num, ao_num)] +&BEGIN_PROVIDER [ integer, max_List_comb_b2_size_thr] + implicit none + integer :: i_1s,i,j,ipoint + double precision :: coef,beta,center(3),int_j1b,thr + double precision :: r(3),weight,dist + thr = 1.d-10 + List_comb_b2_size_thr = 0 + do i = 1, ao_num + do j = i, ao_num + do i_1s = 1, List_all_comb_b2_size + coef = List_all_comb_b2_coef (i_1s) + if(dabs(coef).lt.1.d-10)cycle + beta = List_all_comb_b2_expo (i_1s) + beta = max(beta,1.d-10) + center(1:3) = List_all_comb_b2_cent(1:3,i_1s) + int_j1b = 0.d0 + do ipoint = 1, n_points_final_grid + r(1:3) = final_grid_points(1:3,ipoint) + weight = final_weight_at_r_vector(ipoint) + dist = ( center(1) - r(1) )*( center(1) - r(1) ) + dist += ( center(2) - r(2) )*( center(2) - r(2) ) + dist += ( center(3) - r(3) )*( center(3) - r(3) ) + int_j1b += dabs(aos_in_r_array_transp(ipoint,i) * aos_in_r_array_transp(ipoint,j))*dexp(-beta*dist) * weight + enddo + if(dabs(coef)*dabs(int_j1b).gt.thr)then + List_comb_b2_size_thr(j,i) += 1 + endif + enddo + enddo + enddo + do i = 1, ao_num + do j = 1, i-1 + List_comb_b2_size_thr(j,i) = List_comb_b2_size_thr(i,j) + enddo + enddo + integer :: list(ao_num) + do i = 1, ao_num + list(i) = maxval(List_comb_b2_size_thr(:,i)) + enddo + max_List_comb_b2_size_thr = maxval(list) + +END_PROVIDER + + BEGIN_PROVIDER [ double precision, List_comb_thr_b2_coef, ( max_List_comb_b2_size_thr,ao_num, ao_num )] +&BEGIN_PROVIDER [ double precision, List_comb_thr_b2_expo, ( max_List_comb_b2_size_thr,ao_num, ao_num )] +&BEGIN_PROVIDER [ double precision, List_comb_thr_b2_cent, (3, max_List_comb_b2_size_thr,ao_num, ao_num )] +&BEGIN_PROVIDER [ double precision, ao_abs_comb_b2_j1b, ( max_List_comb_b2_size_thr ,ao_num, ao_num)] + implicit none + integer :: i_1s,i,j,ipoint,icount + double precision :: coef,beta,center(3),int_j1b,thr + double precision :: r(3),weight,dist + thr = 1.d-10 + ao_abs_comb_b2_j1b = 10000000.d0 + do i = 1, ao_num + do j = i, ao_num + icount = 0 + do i_1s = 1, List_all_comb_b2_size + coef = List_all_comb_b2_coef (i_1s) + if(dabs(coef).lt.1.d-10)cycle + beta = List_all_comb_b2_expo (i_1s) + center(1:3) = List_all_comb_b2_cent(1:3,i_1s) + int_j1b = 0.d0 + do ipoint = 1, n_points_final_grid + r(1:3) = final_grid_points(1:3,ipoint) + weight = final_weight_at_r_vector(ipoint) + dist = ( center(1) - r(1) )*( center(1) - r(1) ) + dist += ( center(2) - r(2) )*( center(2) - r(2) ) + dist += ( center(3) - r(3) )*( center(3) - r(3) ) + int_j1b += dabs(aos_in_r_array_transp(ipoint,i) * aos_in_r_array_transp(ipoint,j))*dexp(-beta*dist) * weight + enddo + if(dabs(coef)*dabs(int_j1b).gt.thr)then + icount += 1 + List_comb_thr_b2_coef(icount,j,i) = coef + List_comb_thr_b2_expo(icount,j,i) = beta + List_comb_thr_b2_cent(1:3,icount,j,i) = center(1:3) + ao_abs_comb_b2_j1b(icount,j,i) = int_j1b + endif + enddo + enddo + enddo + + do i = 1, ao_num + do j = 1, i-1 + do icount = 1, List_comb_b2_size_thr(j,i) + List_comb_thr_b2_coef(icount,j,i) = List_comb_thr_b2_coef(icount,i,j) + List_comb_thr_b2_expo(icount,j,i) = List_comb_thr_b2_expo(icount,i,j) + List_comb_thr_b2_cent(1:3,icount,j,i) = List_comb_thr_b2_cent(1:3,icount,i,j) + enddo + enddo + enddo + +END_PROVIDER + + + BEGIN_PROVIDER [ integer, List_comb_b3_size_thr, (ao_num, ao_num)] +&BEGIN_PROVIDER [ integer, max_List_comb_b3_size_thr] + implicit none + integer :: i_1s,i,j,ipoint + double precision :: coef,beta,center(3),int_j1b,thr + double precision :: r(3),weight,dist + thr = 1.d-10 + List_comb_b3_size_thr = 0 + do i = 1, ao_num + do j = 1, ao_num + do i_1s = 1, List_all_comb_b3_size + coef = List_all_comb_b3_coef (i_1s) + beta = List_all_comb_b3_expo (i_1s) + center(1:3) = List_all_comb_b3_cent(1:3,i_1s) + if(dabs(coef).lt.thr)cycle + int_j1b = 0.d0 + do ipoint = 1, n_points_final_grid + r(1:3) = final_grid_points(1:3,ipoint) + weight = final_weight_at_r_vector(ipoint) + dist = ( center(1) - r(1) )*( center(1) - r(1) ) + dist += ( center(2) - r(2) )*( center(2) - r(2) ) + dist += ( center(3) - r(3) )*( center(3) - r(3) ) + int_j1b += dabs(aos_in_r_array_transp(ipoint,i) * aos_in_r_array_transp(ipoint,j))*dexp(-beta*dist) * weight + enddo + if(dabs(coef)*dabs(int_j1b).gt.thr)then + List_comb_b3_size_thr(j,i) += 1 + endif + enddo + enddo + enddo +! do i = 1, ao_num +! do j = 1, i-1 +! List_comb_b3_size_thr(j,i) = List_comb_b3_size_thr(i,j) +! enddo +! enddo + integer :: list(ao_num) + do i = 1, ao_num + list(i) = maxval(List_comb_b3_size_thr(:,i)) + enddo + max_List_comb_b3_size_thr = maxval(list) + print*,'max_List_comb_b3_size_thr = ',max_List_comb_b3_size_thr + +END_PROVIDER + + BEGIN_PROVIDER [ double precision, List_comb_thr_b3_coef, ( max_List_comb_b3_size_thr,ao_num, ao_num )] +&BEGIN_PROVIDER [ double precision, List_comb_thr_b3_expo, ( max_List_comb_b3_size_thr,ao_num, ao_num )] +&BEGIN_PROVIDER [ double precision, List_comb_thr_b3_cent, (3, max_List_comb_b3_size_thr,ao_num, ao_num )] +&BEGIN_PROVIDER [ double precision, ao_abs_comb_b3_j1b, ( max_List_comb_b3_size_thr ,ao_num, ao_num)] + implicit none + integer :: i_1s,i,j,ipoint,icount + double precision :: coef,beta,center(3),int_j1b,thr + double precision :: r(3),weight,dist + thr = 1.d-10 + ao_abs_comb_b3_j1b = 10000000.d0 + do i = 1, ao_num + do j = 1, ao_num + icount = 0 + do i_1s = 1, List_all_comb_b3_size + coef = List_all_comb_b3_coef (i_1s) + beta = List_all_comb_b3_expo (i_1s) + beta = max(beta,1.d-10) + center(1:3) = List_all_comb_b3_cent(1:3,i_1s) + if(dabs(coef).lt.thr)cycle + int_j1b = 0.d0 + do ipoint = 1, n_points_final_grid + r(1:3) = final_grid_points(1:3,ipoint) + weight = final_weight_at_r_vector(ipoint) + dist = ( center(1) - r(1) )*( center(1) - r(1) ) + dist += ( center(2) - r(2) )*( center(2) - r(2) ) + dist += ( center(3) - r(3) )*( center(3) - r(3) ) + int_j1b += dabs(aos_in_r_array_transp(ipoint,i) * aos_in_r_array_transp(ipoint,j))*dexp(-beta*dist) * weight + enddo + if(dabs(coef)*dabs(int_j1b).gt.thr)then + icount += 1 + List_comb_thr_b3_coef(icount,j,i) = coef + List_comb_thr_b3_expo(icount,j,i) = beta + List_comb_thr_b3_cent(1:3,icount,j,i) = center(1:3) + ao_abs_comb_b3_j1b(icount,j,i) = int_j1b + endif + enddo + enddo + enddo + +! do i = 1, ao_num +! do j = 1, i-1 +! do icount = 1, List_comb_b3_size_thr(j,i) +! List_comb_thr_b3_coef(icount,j,i) = List_comb_thr_b3_coef(icount,i,j) +! List_comb_thr_b3_expo(icount,j,i) = List_comb_thr_b3_expo(icount,i,j) +! List_comb_thr_b3_cent(1:3,icount,j,i) = List_comb_thr_b3_cent(1:3,icount,i,j) +! enddo +! enddo +! enddo + +END_PROVIDER + diff --git a/src/ao_one_e_ints/pot_ao_ints.irp.f b/src/ao_one_e_ints/pot_ao_ints.irp.f index 2b6a4d05..dc19f6c7 100644 --- a/src/ao_one_e_ints/pot_ao_ints.irp.f +++ b/src/ao_one_e_ints/pot_ao_ints.irp.f @@ -38,11 +38,6 @@ BEGIN_PROVIDER [ double precision, ao_integrals_n_e, (ao_num,ao_num)] ao_integrals_n_e = 0.d0 - ! _ - ! /| / |_) - ! | / | \ - ! - !$OMP PARALLEL & !$OMP DEFAULT (NONE) & !$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,C_center,power_A,power_B,& @@ -106,7 +101,7 @@ BEGIN_PROVIDER [ double precision, ao_integrals_n_e, (ao_num,ao_num)] endif - IF(DO_PSEUDO) THEN + IF(do_pseudo) THEN ao_integrals_n_e += ao_pseudo_integrals ENDIF @@ -288,8 +283,6 @@ double precision function NAI_pol_mult(A_center,B_center,power_A,power_B,alpha,b ! sum of integrals of type : int {t,[0,1]} exp-(rho.(P-Q)^2 * t^2) * t^i do i =0 ,n_pt_out,2 accu += d(i) * rint(i/2,const) - -! print *, i/2, const, d(i), rint(shiftr(i, 1), const) enddo NAI_pol_mult = accu * coeff diff --git a/src/ao_one_e_ints/pot_ao_pseudo_ints.irp.f b/src/ao_one_e_ints/pot_ao_pseudo_ints.irp.f index 24f43311..1fca6acc 100644 --- a/src/ao_one_e_ints/pot_ao_pseudo_ints.irp.f +++ b/src/ao_one_e_ints/pot_ao_pseudo_ints.irp.f @@ -42,7 +42,7 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integrals_local, (ao_num,ao_num)] double precision :: wall_1, wall_2, wall_0 integer :: thread_num - integer :: omp_get_thread_num + integer, external :: omp_get_thread_num double precision :: c double precision :: Z @@ -169,7 +169,7 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integrals_local, (ao_num,ao_num)] integer :: power_A(3),power_B(3) integer :: i,j,k,l,m double precision :: Vloc, Vpseudo - integer :: omp_get_thread_num + integer, external :: omp_get_thread_num double precision :: wall_1, wall_2, wall_0 integer :: thread_num diff --git a/src/ao_one_e_ints/pseudopot.f90 b/src/ao_one_e_ints/pseudopot.f90 index bad641ab..7321dff7 100644 --- a/src/ao_one_e_ints/pseudopot.f90 +++ b/src/ao_one_e_ints/pseudopot.f90 @@ -1237,7 +1237,7 @@ end integer nptsgridmax,nptsgrid,ik double precision p,q,r,s parameter(nptsgridmax=50) - double precision :: coefs_pseudo(nptsgridmax),ptsgrid(nptsgridmax,3) + double precision coefs_pseudo(nptsgridmax),ptsgrid(nptsgridmax,3) common/pseudos/coefs_pseudo,ptsgrid p=1.d0/dsqrt(2.d0) @@ -1869,7 +1869,7 @@ double precision function int_prod_bessel(l,gam,n,m,a,b,arg) qk = dble(q) two_qkmp1 = 2.d0*(qk+mk)+1.d0 do k=0,q-1 - if (s_q_k < 1.d-32) then + if (s_q_k < 1.d-20) then s_q_k = 0.d0 exit endif diff --git a/src/ao_two_e_ints/map_integrals.irp.f b/src/ao_two_e_ints/map_integrals.irp.f index 0d34d95e..0ce6b87e 100644 --- a/src/ao_two_e_ints/map_integrals.irp.f +++ b/src/ao_two_e_ints/map_integrals.irp.f @@ -142,7 +142,7 @@ subroutine ao_idx2_sq(i,j,ij) ij=i*i endif end - + subroutine idx2_tri_int(i,j,ij) implicit none integer, intent(in) :: i,j @@ -152,7 +152,7 @@ subroutine idx2_tri_int(i,j,ij) q = min(i,j) ij = q+ishft(p*p-p,-1) end - + subroutine ao_idx2_tri_key(i,j,ij) use map_module implicit none @@ -163,8 +163,8 @@ subroutine ao_idx2_tri_key(i,j,ij) q = min(i,j) ij = q+ishft(p*p-p,-1) end - -subroutine two_e_integrals_index_2fold(i,j,k,l,i1) + +subroutine two_e_integrals_index_2fold(i,j,k,l,i1) use map_module implicit none integer, intent(in) :: i,j,k,l @@ -176,7 +176,7 @@ subroutine two_e_integrals_index_2fold(i,j,k,l,i1) call ao_idx2_tri_key(ik,jl,i1) end -subroutine ao_idx2_sq_rev(i,k,ik) +subroutine ao_idx2_sq_rev(i,k,ik) BEGIN_DOC ! reverse square compound index END_DOC @@ -399,7 +399,7 @@ BEGIN_PROVIDER [ complex*16, ao_integrals_cache_periodic, (0:64*64*64*64) ] tmp_im = 0.d0 integral = dcmplx(tmp_re,tmp_im) endif - + ii = l-ao_integrals_cache_min ii = ior( shiftl(ii,6), k-ao_integrals_cache_min) ii = ior( shiftl(ii,6), j-ao_integrals_cache_min) @@ -474,7 +474,7 @@ subroutine get_ao_two_e_integrals(j,k,l,sze,out_val) BEGIN_DOC ! Gets multiple AO bi-electronic integral from the AO map . ! All i are retrieved for j,k,l fixed. - ! physicist convention : + ! physicist convention : END_DOC implicit none integer, intent(in) :: j,k,l, sze @@ -483,7 +483,7 @@ subroutine get_ao_two_e_integrals(j,k,l,sze,out_val) integer :: i integer(key_kind) :: hash logical, external :: ao_one_e_integral_zero - PROVIDE ao_two_e_integrals_in_map ao_integrals_map + PROVIDE ao_two_e_integrals_in_map ao_integrals_map if (ao_one_e_integral_zero(j,l)) then out_val = 0.d0 @@ -503,7 +503,7 @@ subroutine get_ao_two_e_integrals_periodic(j,k,l,sze,out_val) BEGIN_DOC ! Gets multiple AO bi-electronic integral from the AO map . ! All i are retrieved for j,k,l fixed. - ! physicist convention : + ! physicist convention : END_DOC implicit none integer, intent(in) :: j,k,l, sze diff --git a/src/ao_two_e_ints/test_cosgtos_1e.irp.f b/src/ao_two_e_ints/test_cosgtos_1e.irp.f deleted file mode 100644 index 9c1a7215..00000000 --- a/src/ao_two_e_ints/test_cosgtos_1e.irp.f +++ /dev/null @@ -1,191 +0,0 @@ - -! --- - -program test_cosgtos - - implicit none - integer :: i, j - - call init_expo() - -! call test_coef() - call test_1e_kin() - call test_1e_coul() - - i = 1 - j = 1 -! call test_1e_coul_real(i, j) -! call test_1e_coul_cpx (i, j) - -end - -! --- - -subroutine init_expo() - - implicit none - - integer :: i, j - double precision, allocatable :: expo_im(:,:) - - allocate(expo_im(ao_num, ao_prim_num_max)) - - do j = 1, ao_prim_num_max - do i = 1, ao_num - ao_expoim_cosgtos(i,j) = 0.d0 - enddo - enddo - - call ezfio_set_cosgtos_ao_int_ao_expoim_cosgtos(expo_im) - - deallocate(expo_im) - -end subroutine init_expo - -! --- - -subroutine test_coef() - - implicit none - - integer :: i, j - double precision :: coef, coef_gtos, coef_cosgtos - double precision :: delta, accu_abs - - print*, ' check coefs' - - accu_abs = 0.d0 - accu_abs = 0.d0 - do i = 1, ao_num - do j = 1, ao_prim_num(i) - - coef = ao_coef(i,j) - coef_gtos = 1.d0 * ao_coef_normalized_ordered_transp(j,i) - coef_cosgtos = 2.d0 * ao_coef_norm_ord_transp_cosgtos (j,i) - - delta = dabs(coef_gtos - coef_cosgtos) - accu_abs += delta - - if(delta .gt. 1.d-10) then - print*, ' problem on: ' - print*, i, j - print*, coef_gtos, coef_cosgtos, delta - print*, coef - stop - endif - - enddo - enddo - - print*, 'accu_abs = ', accu_abs - -end subroutine test_coef - -! --- - -subroutine test_1e_kin() - - implicit none - - integer :: i, j - double precision :: integral_gtos, integral_cosgtos - double precision :: delta, accu_abs - - print*, ' check kin 1e integrals' - - accu_abs = 0.d0 - accu_abs = 0.d0 - - do j = 1, ao_num - do i = 1, ao_num - - integral_gtos = ao_kinetic_integrals (i,j) - integral_cosgtos = ao_kinetic_integrals_cosgtos(i,j) - - - delta = dabs(integral_gtos - integral_cosgtos) - accu_abs += delta - - if(delta .gt. 1.d-7) then - print*, ' problem on: ' - print*, i, j - print*, integral_gtos, integral_cosgtos, delta - !stop - endif - - enddo - enddo - - print*,'accu_abs = ', accu_abs - -end subroutine test_1e_kin - -! --- - -subroutine test_1e_coul() - - implicit none - - integer :: i, j - double precision :: integral_gtos, integral_cosgtos - double precision :: delta, accu_abs - - print*, ' check Coulomb 1e integrals' - - accu_abs = 0.d0 - accu_abs = 0.d0 - - do j = 1, ao_num - do i = 1, ao_num - - integral_gtos = ao_integrals_n_e (i,j) - integral_cosgtos = ao_integrals_n_e_cosgtos(i,j) - - delta = dabs(integral_gtos - integral_cosgtos) - accu_abs += delta - - if(delta .gt. 1.d-7) then - print*, ' problem on: ' - print*, i, j - print*, integral_gtos, integral_cosgtos, delta - !stop - endif - - enddo - enddo - - print*,'accu_abs = ', accu_abs - -end subroutine test_1e_coul - -! --- - -subroutine test_1e_coul_cpx(i, j) - - implicit none - - integer, intent(in) :: i, j - double precision :: integral - - integral = ao_integrals_n_e_cosgtos(i,j) - - print*, ' cpx Coulomb 1e integrals', integral - -end subroutine test_1e_coul_cpx - -! --- - -subroutine test_1e_coul_real(i, j) - - implicit none - - integer, intent(in) :: i, j - double precision :: integral - - integral = ao_integrals_n_e(i,j) - - print*, ' real Coulomb 1e integrals', integral - -end subroutine test_1e_coul_real - -! --- diff --git a/src/ao_two_e_ints/test_cosgtos_2e.irp.f b/src/ao_two_e_ints/test_cosgtos_2e.irp.f deleted file mode 100644 index de991dd1..00000000 --- a/src/ao_two_e_ints/test_cosgtos_2e.irp.f +++ /dev/null @@ -1,165 +0,0 @@ - -! --- - -program test_cosgtos - - implicit none - integer :: iao, jao, kao, lao - - call init_expo() - -! call test_coef() - call test_2e() - - iao = 1 - jao = 1 - kao = 1 - lao = 21 -! call test_2e_cpx (iao, jao, kao, lao) -! call test_2e_real(iao, jao, kao, lao) - -end - -! --- - -subroutine init_expo() - - implicit none - - integer :: i, j - double precision, allocatable :: expo_im(:,:) - - allocate(expo_im(ao_num, ao_prim_num_max)) - - do j = 1, ao_prim_num_max - do i = 1, ao_num - ao_expoim_cosgtos(i,j) = 0.d0 - enddo - enddo - - call ezfio_set_cosgtos_ao_int_ao_expoim_cosgtos(expo_im) - - deallocate(expo_im) - -end subroutine init_expo - -! --- - -subroutine test_coef() - - implicit none - - integer :: i, j - double precision :: coef, coef_gtos, coef_cosgtos - double precision :: delta, accu_abs - - print*, ' check coefs' - - accu_abs = 0.d0 - accu_abs = 0.d0 - do i = 1, ao_num - do j = 1, ao_prim_num(i) - - coef = ao_coef(i,j) - coef_gtos = 1.d0 * ao_coef_normalized_ordered_transp(j,i) - coef_cosgtos = 2.d0 * ao_coef_norm_ord_transp_cosgtos (j,i) - - delta = dabs(coef_gtos - coef_cosgtos) - accu_abs += delta - - if(delta .gt. 1.d-10) then - print*, ' problem on: ' - print*, i, j - print*, coef_gtos, coef_cosgtos, delta - print*, coef - stop - endif - - enddo - enddo - - print*, 'accu_abs = ', accu_abs - -end subroutine test_coef - - -! --- - -subroutine test_2e() - - implicit none - - integer :: iao, jao, kao, lao - double precision :: integral_gtos, integral_cosgtos - double precision :: delta, accu_abs - - double precision :: ao_two_e_integral, ao_two_e_integral_cosgtos - - print*, ' check integrals' - - accu_abs = 0.d0 - accu_abs = 0.d0 - - ! iao = 1 - ! jao = 1 - ! kao = 1 - ! lao = 24 - - do iao = 1, ao_num ! r1 - do jao = 1, ao_num ! r2 - do kao = 1, ao_num ! r1 - do lao = 1, ao_num ! r2 - - integral_gtos = ao_two_e_integral (iao, kao, jao, lao) - integral_cosgtos = ao_two_e_integral_cosgtos(iao, kao, jao, lao) - - delta = dabs(integral_gtos - integral_cosgtos) - accu_abs += delta - - if(delta .gt. 1.d-7) then - print*, ' problem on: ' - print*, iao, jao, kao, lao - print*, integral_gtos, integral_cosgtos, delta - !stop - endif - - enddo - enddo - enddo - enddo - - print*,'accu_abs = ', accu_abs - -end subroutine test_2e - -! --- - -subroutine test_2e_cpx(iao, jao, kao, lao) - - implicit none - integer, intent(in) :: iao, jao, kao, lao - double precision :: integral - double precision :: ao_two_e_integral_cosgtos - - integral = ao_two_e_integral_cosgtos(iao, kao, jao, lao) - print *, ' cosgtos: ', integral - -end subroutine test_2e_cpx - -! --- - -subroutine test_2e_real(iao, jao, kao, lao) - - implicit none - integer, intent(in) :: iao, jao, kao, lao - double precision :: integral - double precision :: ao_two_e_integral - - integral = ao_two_e_integral(iao, kao, jao, lao) - print *, ' gtos: ', integral - -end subroutine test_2e_real - -! --- - - diff --git a/src/ao_two_e_ints/two_e_integrals.irp.f b/src/ao_two_e_ints/two_e_integrals.irp.f index 80b4af2e..82ffbc90 100644 --- a/src/ao_two_e_ints/two_e_integrals.irp.f +++ b/src/ao_two_e_ints/two_e_integrals.irp.f @@ -603,10 +603,7 @@ double precision function general_primitive_integral(dim, & !DIR$ FORCEINLINE call multiply_poly(d_poly ,n_pt_tmp ,Iz_pol,n_Iz,d1,n_pt_out) double precision :: rint_sum - accu = accu + rint_sum(n_pt_out,const,d1) -! print *, n_pt_out, d1(0:n_pt_out) -! print *, accu general_primitive_integral = fact_p * fact_q * accu *pi_5_2*p_inv*q_inv/dsqrt(p+q) end @@ -871,15 +868,6 @@ subroutine give_polynom_mult_center_x(P_center,Q_center,a_x,d_x,p,q,n_pt_in,pq_i !DIR$ FORCEINLINE call I_x1_pol_mult(a_x,d_x,B10,B01,B00,C00,D00,d,n_pt1,n_pt_in) n_pt_out = n_pt1 - -! print *, ' ' -! print *, a_x, d_x -! print *, B10, B01, B00, C00, D00 -! print *, n_pt1, d(0:n_pt1) -! print *, ' ' - - - if(n_pt1<0)then n_pt_out = -1 do i = 0,n_pt_in diff --git a/src/bitmask/bitmasks_routines.irp.f b/src/bitmask/bitmasks_routines.irp.f index c34d54dc..9c6f4f0c 100644 --- a/src/bitmask/bitmasks_routines.irp.f +++ b/src/bitmask/bitmasks_routines.irp.f @@ -268,6 +268,21 @@ subroutine print_spindet(string,Nint) end +subroutine print_det_one_dimension(string,Nint) + use bitmasks + implicit none + BEGIN_DOC + ! Subroutine to print the content of a determinant using the '+-' notation + END_DOC + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: string(Nint) + character*(2048) :: output(1) + + call bitstring_to_str( output(1), string, Nint ) + print *, trim(output(1)) + +end + logical function is_integer_in_string(bite,string,Nint) use bitmasks implicit none diff --git a/src/cipsi/EZFIO.cfg b/src/cipsi/EZFIO.cfg index 19b45ac1..e01359c5 100644 --- a/src/cipsi/EZFIO.cfg +++ b/src/cipsi/EZFIO.cfg @@ -1,9 +1,3 @@ -[pert_2rdm] -type: logical -doc: If true, computes the one- and two-body rdms with perturbation theory -interface: ezfio,provider,ocaml -default: False - [save_wf_after_selection] type: logical doc: If true, saves the wave function after the selection, before the diagonalization @@ -40,3 +34,9 @@ doc: Maximum number of excitation for beta determinants with respect to the Hart interface: ezfio,ocaml,provider default: -1 +[twice_hierarchy_max] +type: integer +doc: Twice the maximum hierarchy parameter (excitation degree plus half the seniority number). Using -1 selects all determinants +interface: ezfio,ocaml,provider +default: -1 + diff --git a/src/cipsi/NEED b/src/cipsi/NEED index bfbc559a..85d01f79 100644 --- a/src/cipsi/NEED +++ b/src/cipsi/NEED @@ -2,5 +2,4 @@ perturbation zmq mpi iterations -two_body_rdm csf diff --git a/src/cipsi/pt2_stoch_routines.irp.f b/src/cipsi/pt2_stoch_routines.irp.f index b366a268..8d436c19 100644 --- a/src/cipsi/pt2_stoch_routines.irp.f +++ b/src/cipsi/pt2_stoch_routines.irp.f @@ -133,7 +133,7 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in) PROVIDE psi_bilinear_matrix_transp_order psi_selectors_coef_transp psi_det_sorted PROVIDE psi_det_hii selection_weight pseudo_sym PROVIDE n_act_orb n_inact_orb n_core_orb n_virt_orb n_del_orb seniority_max - PROVIDE pert_2rdm excitation_beta_max excitation_alpha_max excitation_max + PROVIDE excitation_beta_max excitation_alpha_max excitation_max if (h0_type == 'CFG') then PROVIDE psi_configuration_hii det_to_configuration @@ -288,12 +288,13 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in) call write_int(6,nproc_target,'Number of threads for PT2') call write_double(6,mem,'Memory (Gb)') - call omp_set_max_active_levels(1) + call set_multiple_levels_omp(.False.) +! call omp_set_max_active_levels(1) - print '(A)', '========== ======================= ===================== ===================== ===========' - print '(A)', ' Samples Energy Variance Norm^2 Seconds' - print '(A)', '========== ======================= ===================== ===================== ===========' + print '(A)', '========== ==================== ================ ================ ================ ============= ===========' + print '(A)', ' Samples Energy PT2 Variance Norm^2 Convergence Seconds' + print '(A)', '========== ==================== ================ ================ ================ ============= ===========' PROVIDE global_selection_buffer @@ -315,14 +316,16 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in) endif !$OMP END PARALLEL call end_parallel_job(zmq_to_qp_run_socket, zmq_socket_pull, 'pt2') - call omp_set_max_active_levels(8) + call set_multiple_levels_omp(.True.) +! call omp_set_max_active_levels(8) - print '(A)', '========== ======================= ===================== ===================== ===========' + print '(A)', '========== ==================== ================ ================ ================ ============= ===========' - do k=1,N_states - pt2_overlap(pt2_stoch_istate,k) = pt2_data % overlap(k,pt2_stoch_istate) - enddo - SOFT_TOUCH pt2_overlap + + do k=1,N_states + pt2_overlap(pt2_stoch_istate,k) = pt2_data % overlap(k,pt2_stoch_istate) + enddo + SOFT_TOUCH pt2_overlap enddo FREE pt2_stoch_istate @@ -415,6 +418,17 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, pt2_data_ double precision :: rss double precision, external :: memory_of_double, memory_of_int + character(len=20) :: format_str1, str_error1, format_str2, str_error2 + character(len=20) :: format_str3, str_error3, format_str4, str_error4 + character(len=20) :: format_value1, format_value2, format_value3, format_value4 + character(len=20) :: str_value1, str_value2, str_value3, str_value4 + character(len=20) :: str_conv + double precision :: value1, value2, value3, value4 + double precision :: error1, error2, error3, error4 + integer :: size1,size2,size3,size4 + + double precision :: conv_crit + sending =.False. rss = memory_of_int(pt2_n_tasks_max*2+N_det_generators*2) @@ -524,28 +538,74 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, pt2_data_ ! 1/(N-1.5) : see Brugger, The American Statistician (23) 4 p. 32 (1969) if(c > 2) then eqt = dabs((pt2_data_S2(t) % pt2(pt2_stoch_istate) / c) - (pt2_data_S(t) % pt2(pt2_stoch_istate)/c)**2) ! dabs for numerical stability - eqt = sqrt(eqt / (dble(c) - 1.5d0)) + eqt = dsqrt(eqt / (dble(c) - 1.5d0)) pt2_data_err % pt2(pt2_stoch_istate) = eqt eqt = dabs((pt2_data_S2(t) % variance(pt2_stoch_istate) / c) - (pt2_data_S(t) % variance(pt2_stoch_istate)/c)**2) ! dabs for numerical stability - eqt = sqrt(eqt / (dble(c) - 1.5d0)) + eqt = dsqrt(eqt / (dble(c) - 1.5d0)) pt2_data_err % variance(pt2_stoch_istate) = eqt eqta(:) = dabs((pt2_data_S2(t) % overlap(:,pt2_stoch_istate) / c) - (pt2_data_S(t) % overlap(:,pt2_stoch_istate)/c)**2) ! dabs for numerical stability - eqta(:) = sqrt(eqta(:) / (dble(c) - 1.5d0)) + eqta(:) = dsqrt(eqta(:) / (dble(c) - 1.5d0)) pt2_data_err % overlap(:,pt2_stoch_istate) = eqta(:) if ((time - time1 > 1.d0) .or. (n==N_det_generators)) then time1 = time - print '(I10, X, F12.6, X, G10.3, X, F10.6, X, G10.3, X, F10.6, X, G10.3, X, F10.4)', c, & - pt2_data % pt2(pt2_stoch_istate) +E, & - pt2_data_err % pt2(pt2_stoch_istate), & - pt2_data % variance(pt2_stoch_istate), & - pt2_data_err % variance(pt2_stoch_istate), & - pt2_data % overlap(pt2_stoch_istate,pt2_stoch_istate), & - pt2_data_err % overlap(pt2_stoch_istate,pt2_stoch_istate), & - time-time0 + + value1 = pt2_data % pt2(pt2_stoch_istate) + E + error1 = pt2_data_err % pt2(pt2_stoch_istate) + value2 = pt2_data % pt2(pt2_stoch_istate) + error2 = pt2_data_err % pt2(pt2_stoch_istate) + value3 = pt2_data % variance(pt2_stoch_istate) + error3 = pt2_data_err % variance(pt2_stoch_istate) + value4 = pt2_data % overlap(pt2_stoch_istate,pt2_stoch_istate) + error4 = pt2_data_err % overlap(pt2_stoch_istate,pt2_stoch_istate) + + ! Max size of the values (FX.Y) with X=size + size1 = 15 + size2 = 12 + size3 = 12 + size4 = 12 + + ! To generate the format: number(error) + call format_w_error(value1,error1,size1,8,format_value1,str_error1) + call format_w_error(value2,error2,size2,8,format_value2,str_error2) + call format_w_error(value3,error3,size3,8,format_value3,str_error3) + call format_w_error(value4,error4,size4,8,format_value4,str_error4) + + ! value > string with the right format + write(str_value1,'('//format_value1//')') value1 + write(str_value2,'('//format_value2//')') value2 + write(str_value3,'('//format_value3//')') value3 + write(str_value4,'('//format_value4//')') value4 + + ! Convergence criterion + conv_crit = dabs(pt2_data_err % pt2(pt2_stoch_istate)) / & + (1.d-20 + dabs(pt2_data % pt2(pt2_stoch_istate)) ) + write(str_conv,'(G10.3)') conv_crit + + write(*,'(I10,X,X,A20,X,A16,X,A16,X,A16,X,A12,X,F10.1)') c,& + adjustl(adjustr(str_value1)//'('//str_error1(1:1)//')'),& + adjustl(adjustr(str_value2)//'('//str_error2(1:1)//')'),& + adjustl(adjustr(str_value3)//'('//str_error3(1:1)//')'),& + adjustl(adjustr(str_value4)//'('//str_error4(1:1)//')'),& + adjustl(str_conv),& + time-time0 + + ! Old print + !print '(I10, X, F12.6, X, G10.3, X, F10.6, X, G10.3, X, F10.6, X, G10.3, X, F10.1,1pE16.6,1pE16.6)', c, & + ! pt2_data % pt2(pt2_stoch_istate) +E, & + ! pt2_data_err % pt2(pt2_stoch_istate), & + ! pt2_data % variance(pt2_stoch_istate), & + ! pt2_data_err % variance(pt2_stoch_istate), & + ! pt2_data % overlap(pt2_stoch_istate,pt2_stoch_istate), & + ! pt2_data_err % overlap(pt2_stoch_istate,pt2_stoch_istate), & + ! time-time0, & + ! pt2_data % pt2(pt2_stoch_istate), & + ! dabs(pt2_data_err % pt2(pt2_stoch_istate)) / & + ! (1.d-20 + dabs(pt2_data % pt2(pt2_stoch_istate)) ) + if (stop_now .or. ( & (do_exit .and. (dabs(pt2_data_err % pt2(pt2_stoch_istate)) / & (1.d-20 + dabs(pt2_data % pt2(pt2_stoch_istate)) ) <= relative_error))) ) then @@ -576,11 +636,11 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, pt2_data_ endif do i=1,n_tasks if(index(i).gt.size(pt2_data_I,1).or.index(i).lt.1)then - print*,'PB !!!' - print*,'If you see this, send a bug report with the following content' - print*,irp_here - print*,'i,index(i),size(pt2_data_I,1) = ',i,index(i),size(pt2_data_I,1) - stop -1 + print*,'PB !!!' + print*,'If you see this, send a bug report with the following content' + print*,irp_here + print*,'i,index(i),size(pt2_data_I,1) = ',i,index(i),size(pt2_data_I,1) + stop -1 endif call pt2_add(pt2_data_I(index(i)),1.d0,pt2_data_task(i)) f(index(i)) -= 1 @@ -843,7 +903,7 @@ END_PROVIDER do t=1, pt2_N_teeth tooth_width = tilde_cW(pt2_n_0(t+1)) - tilde_cW(pt2_n_0(t)) if (tooth_width == 0.d0) then - tooth_width = sum(tilde_w(pt2_n_0(t):pt2_n_0(t+1))) + tooth_width = max(1.d-15,sum(tilde_w(pt2_n_0(t):pt2_n_0(t+1)))) endif ASSERT(tooth_width > 0.d0) do i=pt2_n_0(t)+1, pt2_n_0(t+1) diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f index 00c9dc4d..62d7c52c 100644 --- a/src/cipsi/selection.irp.f +++ b/src/cipsi/selection.irp.f @@ -195,7 +195,10 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d integer :: l_a, nmax, idx integer, allocatable :: indices(:), exc_degree(:), iorder(:) - double precision, parameter :: norm_thr = 1.d-16 + + ! Removed to avoid introducing determinants already presents in the wf + !double precision, parameter :: norm_thr = 1.d-16 + allocate (indices(N_det), & exc_degree(max(N_det_alpha_unique,N_det_beta_unique))) @@ -215,10 +218,11 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d i = psi_bilinear_matrix_rows(l_a) if (nt + exc_degree(i) <= 4) then idx = psi_det_sorted_order(psi_bilinear_matrix_order(l_a)) - if (psi_average_norm_contrib_sorted(idx) > norm_thr) then + ! Removed to avoid introducing determinants already presents in the wf + !if (psi_average_norm_contrib_sorted(idx) > norm_thr) then indices(k) = idx k=k+1 - endif + !endif endif enddo enddo @@ -242,10 +246,11 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d idx = psi_det_sorted_order( & psi_bilinear_matrix_order( & psi_bilinear_matrix_transp_order(l_a))) - if (psi_average_norm_contrib_sorted(idx) > norm_thr) then + ! Removed to avoid introducing determinants already presents in the wf + !if (psi_average_norm_contrib_sorted(idx) > norm_thr) then indices(k) = idx k=k+1 - endif + !endif endif enddo enddo @@ -464,27 +469,21 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d allocate (fullminilist (N_int, 2, fullinteresting(0)), & minilist (N_int, 2, interesting(0)) ) - if(pert_2rdm)then - allocate(coef_fullminilist_rev(N_states,fullinteresting(0))) - do i=1,fullinteresting(0) - do j = 1, N_states - coef_fullminilist_rev(j,i) = psi_coef_sorted(fullinteresting(i),j) - enddo - enddo - endif +! if(pert_2rdm)then +! allocate(coef_fullminilist_rev(N_states,fullinteresting(0))) +! do i=1,fullinteresting(0) +! do j = 1, N_states +! coef_fullminilist_rev(j,i) = psi_coef_sorted(fullinteresting(i),j) +! enddo +! enddo +! endif do i=1,fullinteresting(0) - do k=1,N_int - fullminilist(k,1,i) = psi_det_sorted(k,1,fullinteresting(i)) - fullminilist(k,2,i) = psi_det_sorted(k,2,fullinteresting(i)) - enddo + fullminilist(:,:,i) = psi_det_sorted(:,:,fullinteresting(i)) enddo do i=1,interesting(0) - do k=1,N_int - minilist(k,1,i) = psi_det_sorted(k,1,interesting(i)) - minilist(k,2,i) = psi_det_sorted(k,2,interesting(i)) - enddo + minilist(:,:,i) = psi_det_sorted(:,:,interesting(i)) enddo do s2=s1,2 @@ -531,19 +530,19 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d call splash_pq(mask, sp, minilist, i_generator, interesting(0), bannedOrb, banned, mat, interesting) - if(.not.pert_2rdm)then +! if(.not.pert_2rdm)then call fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf) - else - call fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf,fullminilist, coef_fullminilist_rev, fullinteresting(0)) - endif +! else +! call fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf,fullminilist, coef_fullminilist_rev, fullinteresting(0)) +! endif end if enddo if(s1 /= s2) monoBdo = .false. enddo deallocate(fullminilist,minilist) - if(pert_2rdm)then - deallocate(coef_fullminilist_rev) - endif +! if(pert_2rdm)then +! deallocate(coef_fullminilist_rev) +! endif enddo enddo deallocate(preinteresting, prefullinteresting, interesting, fullinteresting) @@ -713,6 +712,25 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d if (do_cycle) cycle endif + if (twice_hierarchy_max >= 0) then + s = 0 + do k=1,N_int + s = s + popcnt(ieor(det(k,1),det(k,2))) + enddo + if ( mod(s,2)>0 ) stop 'For now, hierarchy CI is defined only for an even number of electrons' + if (excitation_ref == 1) then + call get_excitation_degree(HF_bitmask,det(1,1),degree,N_int) + else if (excitation_ref == 2) then + stop 'For now, hierarchy CI is defined only for a single reference determinant' +! do k=1,N_dominant_dets_of_cfgs +! call get_excitation_degree(dominant_dets_of_cfgs(1,1,k),det(1,1),degree,N_int) +! enddo + endif + integer :: twice_hierarchy + twice_hierarchy = degree + s/2 + if (twice_hierarchy > twice_hierarchy_max) cycle + endif + Hii = diag_H_mat_elem_fock(psi_det_generators(1,1,i_generator),det,fock_diag_tmp,N_int) w = 0d0 @@ -837,8 +855,28 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d endif end select + + ! To force the inclusion of determinants with a positive pt2 contribution + if (e_pert(istate) > 1d-8) then + w = -huge(1.0) + endif + end do +!!!BEGIN_DEBUG +! ! To check if the pt2 is taking determinants already in the wf +! if (is_in_wavefunction(det(N_int,1),N_int)) then +! logical, external :: is_in_wavefunction +! print*, 'A determinant contributing to the pt2 is already in' +! print*, 'the wave function:' +! call print_det(det(N_int,1),N_int) +! print*,'contribution to the pt2 for the states:', e_pert(:) +! print*,'error in the filtering in' +! print*, 'cipsi/selection.irp.f sub: selecte_singles_and_doubles' +! print*, 'abort' +! call abort +! endif +!!!END_DEBUG integer(bit_kind) :: occ(N_int,2), n if (h0_type == 'CFG') then @@ -1559,7 +1597,7 @@ subroutine bitstring_to_list_in_selection( string, list, n_elements, Nint) use bitmasks implicit none BEGIN_DOC - ! Gives the inidices(+1) of the bits set to 1 in the bit string + ! Gives the indices(+1) of the bits set to 1 in the bit string END_DOC integer, intent(in) :: Nint integer(bit_kind), intent(in) :: string(Nint) diff --git a/src/cipsi/selection_buffer.irp.f b/src/cipsi/selection_buffer.irp.f index 10132086..1f743e0e 100644 --- a/src/cipsi/selection_buffer.irp.f +++ b/src/cipsi/selection_buffer.irp.f @@ -60,6 +60,7 @@ subroutine add_to_selection_buffer(b, det, val) b%val(b%cur) = val if(b%cur == size(b%val)) then call sort_selection_buffer(b) + b%cur = b%cur-1 end if end if end subroutine @@ -86,43 +87,56 @@ subroutine merge_selection_buffers(b1, b2) double precision :: rss double precision, external :: memory_of_double sze = max(size(b1%val), size(b2%val)) - rss = memory_of_double(sze) + 2*N_int*memory_of_double(sze) - call check_mem(rss,irp_here) +! rss = memory_of_double(sze) + 2*N_int*memory_of_double(sze) +! call check_mem(rss,irp_here) allocate(val(sze), detmp(N_int, 2, sze)) i1=1 i2=1 - do i=1,nmwen - if ( (i1 > b1%cur).and.(i2 > b2%cur) ) then - exit - else if (i1 > b1%cur) then - val(i) = b2%val(i2) - detmp(1:N_int,1,i) = b2%det(1:N_int,1,i2) - detmp(1:N_int,2,i) = b2%det(1:N_int,2,i2) - i2=i2+1 - else if (i2 > b2%cur) then - val(i) = b1%val(i1) - detmp(1:N_int,1,i) = b1%det(1:N_int,1,i1) - detmp(1:N_int,2,i) = b1%det(1:N_int,2,i1) - i1=i1+1 - else - if (b1%val(i1) <= b2%val(i2)) then - val(i) = b1%val(i1) - detmp(1:N_int,1,i) = b1%det(1:N_int,1,i1) - detmp(1:N_int,2,i) = b1%det(1:N_int,2,i1) - i1=i1+1 + + select case (N_int) +BEGIN_TEMPLATE + case $case + do i=1,nmwen + if ( (i1 > b1%cur).and.(i2 > b2%cur) ) then + exit + else if (i1 > b1%cur) then + val(i) = b2%val(i2) + detmp(1:$N_int,1,i) = b2%det(1:$N_int,1,i2) + detmp(1:$N_int,2,i) = b2%det(1:$N_int,2,i2) + i2=i2+1 + else if (i2 > b2%cur) then + val(i) = b1%val(i1) + detmp(1:$N_int,1,i) = b1%det(1:$N_int,1,i1) + detmp(1:$N_int,2,i) = b1%det(1:$N_int,2,i1) + i1=i1+1 else - val(i) = b2%val(i2) - detmp(1:N_int,1,i) = b2%det(1:N_int,1,i2) - detmp(1:N_int,2,i) = b2%det(1:N_int,2,i2) - i2=i2+1 + if (b1%val(i1) <= b2%val(i2)) then + val(i) = b1%val(i1) + detmp(1:$N_int,1,i) = b1%det(1:$N_int,1,i1) + detmp(1:$N_int,2,i) = b1%det(1:$N_int,2,i1) + i1=i1+1 + else + val(i) = b2%val(i2) + detmp(1:$N_int,1,i) = b2%det(1:$N_int,1,i2) + detmp(1:$N_int,2,i) = b2%det(1:$N_int,2,i2) + i2=i2+1 + endif endif - endif - enddo + enddo + do i=nmwen+1,b2%N + val(i) = 0.d0 +! detmp(1:$N_int,1,i) = 0_bit_kind +! detmp(1:$N_int,2,i) = 0_bit_kind + enddo +SUBST [ case, N_int ] +(1); 1;; +(2); 2;; +(3); 3;; +(4); 4;; +default; N_int;; +END_TEMPLATE + end select deallocate(b2%det, b2%val) - do i=nmwen+1,b2%N - val(i) = 0.d0 - detmp(1:N_int,1:2,i) = 0_bit_kind - enddo b2%det => detmp b2%val => val b2%mini = min(b2%mini,b2%val(b2%N)) @@ -144,8 +158,8 @@ subroutine sort_selection_buffer(b) double precision :: rss double precision, external :: memory_of_double, memory_of_int - rss = memory_of_int(b%cur) + 2*N_int*memory_of_double(size(b%det,3)) - call check_mem(rss,irp_here) +! rss = memory_of_int(b%cur) + 2*N_int*memory_of_double(size(b%det,3)) +! call check_mem(rss,irp_here) allocate(iorder(b%cur), detmp(N_int, 2, size(b%det,3))) do i=1,b%cur iorder(i) = i @@ -225,14 +239,14 @@ subroutine make_selection_buffer_s2(b) endif dup = .True. do k=1,N_int - if ( (tmp_array(k,1,i) /= tmp_array(k,1,j)) & - .or. (tmp_array(k,2,i) /= tmp_array(k,2,j)) ) then + if ( (tmp_array(k,1,i) /= tmp_array(k,1,j)) .or. & + (tmp_array(k,2,i) /= tmp_array(k,2,j)) ) then dup = .False. exit endif enddo if (dup) then - val(i) = max(val(i), val(j)) + val(i) = min(val(i), val(j)) duplicate(j) = .True. endif j+=1 @@ -282,9 +296,6 @@ subroutine make_selection_buffer_s2(b) call configuration_to_dets_size(o(1,1,i),sze,elec_alpha_num,N_int) n_d = n_d + sze if (n_d > b%cur) then -! if (n_d - b%cur > b%cur - n_d + sze) then -! n_d = n_d - sze -! endif exit endif enddo @@ -329,10 +340,11 @@ subroutine remove_duplicates_in_selection_buffer(b) integer(bit_kind), allocatable :: tmp_array(:,:,:) logical, allocatable :: duplicate(:) - n_d = b%cur logical :: found_duplicates double precision :: rss double precision, external :: memory_of_double + + n_d = b%cur rss = (4*N_int+4)*memory_of_double(n_d) call check_mem(rss,irp_here) diff --git a/src/cipsi/selection_weight.irp.f b/src/cipsi/selection_weight.irp.f index 3c09e59a..756c65a1 100644 --- a/src/cipsi/selection_weight.irp.f +++ b/src/cipsi/selection_weight.irp.f @@ -38,11 +38,11 @@ subroutine update_pt2_and_variance_weights(pt2_data, N_st) avg = sum(pt2(1:N_st)) / dble(N_st) + 1.d-32 ! Avoid future division by zero - dt = 8.d0 !* selection_factor + dt = 4.d0 !* selection_factor do k=1,N_st - element = exp(dt*(pt2(k)/avg - 1.d0)) - element = min(2.0d0 , element) - element = max(0.5d0 , element) + element = pt2(k) !exp(dt*(pt2(k)/avg - 1.d0)) +! element = min(2.0d0 , element) +! element = max(0.5d0 , element) pt2_match_weight(k) *= element enddo @@ -50,9 +50,9 @@ subroutine update_pt2_and_variance_weights(pt2_data, N_st) avg = sum(variance(1:N_st)) / dble(N_st) + 1.d-32 ! Avoid future division by zero do k=1,N_st - element = exp(dt*(variance(k)/avg -1.d0)) - element = min(2.0d0 , element) - element = max(0.5d0 , element) + element = variance(k) ! exp(dt*(variance(k)/avg -1.d0)) +! element = min(2.0d0 , element) +! element = max(0.5d0 , element) variance_match_weight(k) *= element enddo @@ -62,6 +62,9 @@ subroutine update_pt2_and_variance_weights(pt2_data, N_st) variance_match_weight(:) = 1.d0 endif + pt2_match_weight(:) = pt2_match_weight(:)/sum(pt2_match_weight(:)) + variance_match_weight(:) = variance_match_weight(:)/sum(variance_match_weight(:)) + threshold_davidson_pt2 = min(1.d-6, & max(threshold_davidson, 1.e-1 * PT2_relative_error * minval(abs(pt2(1:N_states)))) ) @@ -87,7 +90,7 @@ BEGIN_PROVIDER [ double precision, selection_weight, (N_states) ] selection_weight(1:N_states) = c0_weight(1:N_states) case (2) - print *, 'Using pt2-matching weight in selection' + print *, 'Using PT2-matching weight in selection' selection_weight(1:N_states) = c0_weight(1:N_states) * pt2_match_weight(1:N_states) print *, '# PT2 weight ', real(pt2_match_weight(:),4) @@ -97,7 +100,7 @@ BEGIN_PROVIDER [ double precision, selection_weight, (N_states) ] print *, '# var weight ', real(variance_match_weight(:),4) case (4) - print *, 'Using variance- and pt2-matching weights in selection' + print *, 'Using variance- and PT2-matching weights in selection' selection_weight(1:N_states) = c0_weight(1:N_states) * sqrt(variance_match_weight(1:N_states) * pt2_match_weight(1:N_states)) print *, '# PT2 weight ', real(pt2_match_weight(:),4) print *, '# var weight ', real(variance_match_weight(:),4) @@ -112,7 +115,7 @@ BEGIN_PROVIDER [ double precision, selection_weight, (N_states) ] selection_weight(1:N_states) = c0_weight(1:N_states) case (7) - print *, 'Input weights multiplied by variance- and pt2-matching' + print *, 'Input weights multiplied by variance- and PT2-matching' selection_weight(1:N_states) = c0_weight(1:N_states) * sqrt(variance_match_weight(1:N_states) * pt2_match_weight(1:N_states)) * state_average_weight(1:N_states) print *, '# PT2 weight ', real(pt2_match_weight(:),4) print *, '# var weight ', real(variance_match_weight(:),4) @@ -128,6 +131,7 @@ BEGIN_PROVIDER [ double precision, selection_weight, (N_states) ] print *, '# var weight ', real(variance_match_weight(:),4) end select + selection_weight(:) = selection_weight(:)/sum(selection_weight(:)) print *, '# Total weight ', real(selection_weight(:),4) END_PROVIDER diff --git a/src/cipsi/slave_cipsi.irp.f b/src/cipsi/slave_cipsi.irp.f index 510c667b..ddfc050e 100644 --- a/src/cipsi/slave_cipsi.irp.f +++ b/src/cipsi/slave_cipsi.irp.f @@ -4,7 +4,7 @@ subroutine run_slave_cipsi ! Helper program for distributed parallelism END_DOC - call omp_set_max_active_levels(1) + call set_multiple_levels_omp(.False.) distributed_davidson = .False. read_wf = .False. SOFT_TOUCH read_wf distributed_davidson @@ -171,9 +171,9 @@ subroutine run_slave_main call write_double(6,(t1-t0),'Broadcast time') !--- - call omp_set_max_active_levels(8) + call set_multiple_levels_omp(.True.) call davidson_slave_tcp(0) - call omp_set_max_active_levels(1) + call set_multiple_levels_omp(.False.) print *, mpi_rank, ': Davidson done' !--- diff --git a/src/cipsi/zmq_selection.irp.f b/src/cipsi/zmq_selection.irp.f index 58630709..1bfe87c0 100644 --- a/src/cipsi/zmq_selection.irp.f +++ b/src/cipsi/zmq_selection.irp.f @@ -22,7 +22,7 @@ subroutine ZMQ_selection(N_in, pt2_data) PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns PROVIDE psi_bilinear_matrix_transp_order selection_weight pseudo_sym PROVIDE n_act_orb n_inact_orb n_core_orb n_virt_orb n_del_orb seniority_max - PROVIDE pert_2rdm excitation_beta_max excitation_alpha_max excitation_max + PROVIDE excitation_beta_max excitation_alpha_max excitation_max call new_parallel_job(zmq_to_qp_run_socket,zmq_socket_pull,'selection') diff --git a/src/cis/cis.irp.f b/src/cis/cis.irp.f index ab2294ad..2b16a5f7 100644 --- a/src/cis/cis.irp.f +++ b/src/cis/cis.irp.f @@ -62,6 +62,7 @@ subroutine run else call H_apply_cis endif + print*,'' print *, 'N_det = ', N_det print*,'******************************' print *, 'Energies of the states:' @@ -69,11 +70,13 @@ subroutine run print *, i, CI_energy(i) enddo if (N_states > 1) then - print*,'******************************' - print*,'Excitation energies ' + print*,'' + print*,'******************************************************' + print*,'Excitation energies (au) (eV)' do i = 2, N_states - print*, i ,CI_energy(i) - CI_energy(1) + print*, i ,CI_energy(i) - CI_energy(1), (CI_energy(i) - CI_energy(1)) * ha_to_ev enddo + print*,'' endif call ezfio_set_cis_energy(CI_energy) diff --git a/src/cisd/30.cisd.bats b/src/cisd/30.cisd.bats index f4e8018b..69b862b0 100644 --- a/src/cisd/30.cisd.bats +++ b/src/cisd/30.cisd.bats @@ -77,7 +77,7 @@ function run() { [[ -n $TRAVIS ]] && skip qp set_file ch4.ezfio qp set_mo_class --core="[1]" --act="[2-30]" --del="[31-59]" - run -40.2403962667047 -39.843315 + run -40.2403962667047 -39.8433221754964 } @test "SiH3" { # 20.2202s 1.38648m diff --git a/src/cisd/cisd.irp.f b/src/cisd/cisd.irp.f index 6c55e2ff..fca3b10e 100644 --- a/src/cisd/cisd.irp.f +++ b/src/cisd/cisd.irp.f @@ -63,7 +63,7 @@ subroutine run endif psi_coef = ci_eigenvectors SOFT_TOUCH psi_coef - call save_wavefunction + call save_wavefunction_truncated(save_threshold) call ezfio_set_cisd_energy(CI_energy) do i = 1,N_states diff --git a/src/csf/configurations.irp.f b/src/csf/configurations.irp.f index c11a49a4..aebf53d9 100644 --- a/src/csf/configurations.irp.f +++ b/src/csf/configurations.irp.f @@ -856,6 +856,7 @@ end subroutine !end subroutine ! BEGIN_PROVIDER [ integer, psi_configuration_to_psi_det, (2,N_configuration) ] +&BEGIN_PROVIDER [ integer, psi_configuration_n_det, (N_configuration) ] &BEGIN_PROVIDER [ integer, psi_configuration_to_psi_det_data, (N_det) ] implicit none @@ -944,6 +945,29 @@ end subroutine enddo deallocate(dets, old_order) + integer :: ndet_conf + do i = 1, N_configuration + ndet_conf = psi_configuration_to_psi_det(2,i) - psi_configuration_to_psi_det(1,i) + 1 + psi_configuration_n_det(i) = ndet_conf + enddo END_PROVIDER + +BEGIN_PROVIDER [ integer, n_elec_alpha_for_psi_configuration, (N_configuration)] + implicit none + integer :: i,j,k,l + integer(bit_kind) :: det_tmp(N_int,2),det_alpha(N_int) + n_elec_alpha_for_psi_configuration = 0 + do i = 1, N_configuration + j = psi_configuration_to_psi_det(2,i) + det_tmp(:,:) = psi_det(:,:,j) + k = 0 + do l = 1, N_int + det_alpha(N_int) = iand(det_tmp(l,1),psi_configuration(l,1,i)) + k += popcnt(det_alpha(l)) + enddo + n_elec_alpha_for_psi_configuration(i) = k + enddo + +END_PROVIDER diff --git a/src/csf/conversion.irp.f b/src/csf/conversion.irp.f index bdd8b327..494c3bfa 100644 --- a/src/csf/conversion.irp.f +++ b/src/csf/conversion.irp.f @@ -10,7 +10,6 @@ BEGIN_PROVIDER [ double precision, psi_csf_coef, (N_csf, N_states) ] call convertWFfromDETtoCSF(N_states, buffer, psi_csf_coef) END_PROVIDER - subroutine convertWFfromDETtoCSF(N_st,psi_coef_det_in, psi_coef_cfg_out) use cfunctions use bitmasks diff --git a/src/dav_general_mat/dav_general.irp.f b/src/dav_general_mat/dav_general.irp.f index 96775c50..e899e0fb 100644 --- a/src/dav_general_mat/dav_general.irp.f +++ b/src/dav_general_mat/dav_general.irp.f @@ -228,7 +228,7 @@ subroutine davidson_general(u_in,H_jj,energies,dim_in,sze,N_st,N_st_diag_in,conv shift = N_st_diag*(iter-1) shift2 = N_st_diag*iter - if ((iter > 1).or.(itertot == 1)) then +! if ((iter > 1).or.(itertot == 1)) then ! Compute |W_k> = \sum_i |i> ! ----------------------------------- @@ -238,10 +238,10 @@ subroutine davidson_general(u_in,H_jj,energies,dim_in,sze,N_st,N_st_diag_in,conv ! call H_S2_u_0_nstates_openmp(W(:,shift+1),U(:,shift+1),N_st_diag,sze) call hpsi(W(:,shift+1),U(:,shift+1),N_st_diag,sze,h_mat) - else - ! Already computed in update below - continue - endif +! else +! ! Already computed in update below +! continue +! endif ! Compute h_kl = = ! ------------------------------------------- diff --git a/src/davidson/davidson_parallel.irp.f b/src/davidson/davidson_parallel.irp.f index 9e212094..4f88506a 100644 --- a/src/davidson/davidson_parallel.irp.f +++ b/src/davidson/davidson_parallel.irp.f @@ -508,7 +508,7 @@ subroutine H_S2_u_0_nstates_zmq(v_0,s_0,u_0,N_st,sze) endif - call omp_set_max_active_levels(5) + call set_multiple_levels_omp(.True.) !$OMP PARALLEL DEFAULT(shared) NUM_THREADS(2) PRIVATE(ithread) ithread = omp_get_thread_num() @@ -546,21 +546,6 @@ end -!BEGIN_PROVIDER [ integer, nthreads_davidson ] -! implicit none -! BEGIN_DOC -! ! Number of threads for Davidson -! END_DOC -! nthreads_davidson = nproc -! character*(32) :: env -! call getenv('QP_NTHREADS_DAVIDSON',env) -! if (trim(env) /= '') then -! read(env,*) nthreads_davidson -! call write_int(6,nthreads_davidson,'Target number of threads for ') -! endif -!END_PROVIDER - - integer function zmq_put_N_states_diag(zmq_to_qp_run_socket,worker_id) use f77_zmq implicit none diff --git a/src/davidson/davidson_parallel_csf.irp.f b/src/davidson/davidson_parallel_csf.irp.f index fe651b1d..d8e9bffa 100644 --- a/src/davidson/davidson_parallel_csf.irp.f +++ b/src/davidson/davidson_parallel_csf.irp.f @@ -464,7 +464,8 @@ subroutine H_u_0_nstates_zmq(v_0,u_0,N_st,sze) print *, irp_here, ': Failed in zmq_set_running' endif - call omp_set_max_active_levels(4) + call set_multiple_levels_omp(.True.) + !$OMP PARALLEL DEFAULT(shared) NUM_THREADS(2) PRIVATE(ithread) ithread = omp_get_thread_num() if (ithread == 0 ) then diff --git a/src/davidson/davidson_parallel_nos2.irp.f b/src/davidson/davidson_parallel_nos2.irp.f index 84cbe3af..597b001f 100644 --- a/src/davidson/davidson_parallel_nos2.irp.f +++ b/src/davidson/davidson_parallel_nos2.irp.f @@ -464,7 +464,8 @@ subroutine H_u_0_nstates_zmq(v_0,u_0,N_st,sze) print *, irp_here, ': Failed in zmq_set_running' endif - call omp_set_max_active_levels(4) + call set_multiple_levels_omp(.True.) + !$OMP PARALLEL DEFAULT(shared) NUM_THREADS(2) PRIVATE(ithread) ithread = omp_get_thread_num() if (ithread == 0 ) then diff --git a/src/davidson/diagonalization_hcsf_dressed.irp.f b/src/davidson/diagonalization_hcsf_dressed.irp.f index 3020ecd8..f531fb3a 100644 --- a/src/davidson/diagonalization_hcsf_dressed.irp.f +++ b/src/davidson/diagonalization_hcsf_dressed.irp.f @@ -300,7 +300,7 @@ subroutine davidson_diag_csf_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N shift = N_st_diag*(iter-1) shift2 = N_st_diag*iter - if ((iter > 1).or.(itertot == 1)) then +! if ((iter > 1).or.(itertot == 1)) then ! Compute |W_k> = \sum_i |i> ! ----------------------------------- @@ -310,10 +310,10 @@ subroutine davidson_diag_csf_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N else call H_u_0_nstates_openmp(W,U,N_st_diag,sze) endif - else - ! Already computed in update below - continue - endif +! else +! ! Already computed in update below +! continue +! endif if (dressing_state > 0) then diff --git a/src/davidson/diagonalization_hs2_dressed.irp.f b/src/davidson/diagonalization_hs2_dressed.irp.f index 0fd9f091..f4c05595 100644 --- a/src/davidson/diagonalization_hs2_dressed.irp.f +++ b/src/davidson/diagonalization_hs2_dressed.irp.f @@ -14,15 +14,6 @@ BEGIN_PROVIDER [ character*(64), diag_algorithm ] endif END_PROVIDER -!BEGIN_PROVIDER [ double precision, threshold_davidson_pt2 ] -! implicit none -! BEGIN_DOC -! ! Threshold of Davidson's algorithm, using PT2 as a guide -! END_DOC -! threshold_davidson_pt2 = threshold_davidson -! -!END_PROVIDER - BEGIN_PROVIDER [ integer, dressed_column_idx, (N_states) ] @@ -86,9 +77,9 @@ subroutine davidson_diag_hs2(dets_in,u_in,s2_out,dim_in,energies,sze,N_st,N_st_d !$OMP END PARALLEL if (dressing_state > 0) then - do k = 1, N_st + do k=1,N_st - do i = 1, sze + do i=1,sze H_jj(i) += u_in(i,k) * dressing_column_h(i,k) enddo @@ -154,7 +145,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_ character*(16384) :: write_buffer double precision :: to_print(3,N_st) double precision :: cpu, wall - integer :: shift, shift2, itermax, istate, ii + integer :: shift, shift2, itermax, istate double precision :: r1, r2, alpha logical :: state_ok(N_st_diag_in*davidson_sze_max) integer :: nproc_target @@ -354,27 +345,20 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_ shift = N_st_diag*(iter-1) shift2 = N_st_diag*iter - if ((iter > 1).or.(itertot == 1)) then +! if ((iter > 1).or.(itertot == 1)) then ! Compute |W_k> = \sum_i |i> ! ----------------------------------- if ((sze > 100000).and.distributed_davidson) then call H_S2_u_0_nstates_zmq (W(1,shift+1),S_d,U(1,shift+1),N_st_diag,sze) else - double precision :: irp_rdtsc - double precision :: ticks_0, ticks_1 - integer*8 :: irp_imax - irp_imax = 1 - !ticks_0 = irp_rdtsc() call H_S2_u_0_nstates_openmp(W(1,shift+1),S_d,U(1,shift+1),N_st_diag,sze) - !ticks_1 = irp_rdtsc() - !print *,' ----Cycles:',(ticks_1-ticks_0)/dble(irp_imax)," ----" endif S(1:sze,shift+1:shift+N_st_diag) = real(S_d(1:sze,1:N_st_diag)) - else - ! Already computed in update below - continue - endif +! else +! ! Already computed in update below +! continue +! endif if (dressing_state > 0) then diff --git a/src/davidson/diagonalization_nonsym_h_dressed.irp.f b/src/davidson/diagonalization_nonsym_h_dressed.irp.f index cd576b02..3ff060a6 100644 --- a/src/davidson/diagonalization_nonsym_h_dressed.irp.f +++ b/src/davidson/diagonalization_nonsym_h_dressed.irp.f @@ -317,7 +317,7 @@ subroutine davidson_diag_nonsym_hjj(dets_in, u_in, H_jj, energies, dim_in, sze, shift = N_st_diag*(iter-1) shift2 = N_st_diag*iter - if( (iter > 1) .or. (itertot == 1) ) then +! if( (iter > 1) .or. (itertot == 1) ) then ! Gram-Schmidt to orthogonalize all new guess with the previous vectors call ortho_qr(U, size(U, 1), sze, shift2) @@ -331,10 +331,10 @@ subroutine davidson_diag_nonsym_hjj(dets_in, u_in, H_jj, energies, dim_in, sze, else call H_u_0_nstates_openmp(W(1,shift+1), U(1,shift+1), N_st_diag, sze) endif - else - ! Already computed in update below - continue - endif +! else +! ! Already computed in update below +! continue +! endif if(dressing_state > 0) then diff --git a/src/davidson/diagonalize_ci.irp.f b/src/davidson/diagonalize_ci.irp.f index a34637a0..76d8b65f 100644 --- a/src/davidson/diagonalize_ci.irp.f +++ b/src/davidson/diagonalize_ci.irp.f @@ -299,6 +299,7 @@ subroutine diagonalize_CI ! eigenstates of the |CI| matrix. END_DOC integer :: i,j + PROVIDE distributed_davidson do j=1,N_states do i=1,N_det psi_coef(i,j) = CI_eigenvectors(i,j) diff --git a/src/davidson/input.irp.f b/src/davidson/input.irp.f deleted file mode 100644 index 83f5c09e..00000000 --- a/src/davidson/input.irp.f +++ /dev/null @@ -1,39 +0,0 @@ -!BEGIN_PROVIDER [ integer, n_states_diag ] -! implicit none -! BEGIN_DOC -!! Number of states to consider during the Davdison diagonalization -! END_DOC -! -! logical :: has -! PROVIDE ezfio_filename -! if (mpi_master) then -! -! call ezfio_has_davidson_n_states_diag(has) -! if (has) then -! call ezfio_get_davidson_n_states_diag(n_states_diag) -! else -! print *, 'davidson/n_states_diag not found in EZFIO file' -! stop 1 -! endif -! n_states_diag = max(2,N_states * N_states_diag) -! endif -! IRP_IF MPI_DEBUG -! print *, irp_here, mpi_rank -! call MPI_BARRIER(MPI_COMM_WORLD, ierr) -! IRP_ENDIF -! IRP_IF MPI -! include 'mpif.h' -! integer :: ierr -! call MPI_BCAST( n_states_diag, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) -! if (ierr /= MPI_SUCCESS) then -! stop 'Unable to read n_states_diag with MPI' -! endif -! IRP_ENDIF -! -! call write_time(6) -! if (mpi_master) then -! write(6, *) 'Read n_states_diag' -! endif -! -!END_PROVIDER -! diff --git a/src/davidson_keywords/README.rst b/src/davidson_keywords/README.rst index cdc7cc1f..9567cdb1 100644 --- a/src/davidson_keywords/README.rst +++ b/src/davidson_keywords/README.rst @@ -2,3 +2,4 @@ davidson_keywords ================= +Keywords used for Davidson algorithms. diff --git a/src/determinants/EZFIO.cfg b/src/determinants/EZFIO.cfg index c6323cd0..9eefa66c 100644 --- a/src/determinants/EZFIO.cfg +++ b/src/determinants/EZFIO.cfg @@ -42,13 +42,13 @@ default: 2 [weight_selection] type: integer -doc: Weight used in the selection. 0: input state-average weight, 1: 1./(c_0^2), 2: rPT2 matching, 3: variance matching, 4: variance and rPT2 matching, 5: variance minimization and matching, 6: CI coefficients 7: input state-average multiplied by variance and rPT2 matching 8: input state-average multiplied by rPT2 matching 9: input state-average multiplied by variance matching +doc: Weight used in the selection. 0: input state-average weight, 1: 1./(c_0^2), 2: PT2 matching, 3: variance matching, 4: variance and PT2 matching, 5: variance minimization and matching, 6: CI coefficients 7: input state-average multiplied by variance and PT2 matching 8: input state-average multiplied by PT2 matching 9: input state-average multiplied by variance matching interface: ezfio,provider,ocaml default: 1 [threshold_generators] type: Threshold -doc: Thresholds on generators (fraction of the square of the norm) +doc: Thresholds on generators (fraction of the square of the norm) interface: ezfio,provider,ocaml default: 0.999 @@ -80,7 +80,7 @@ type: integer [psi_coef] interface: ezfio doc: Coefficients of the wave function -type: double precision +type: double precision size: (determinants.n_det,determinants.n_states) [psi_det] @@ -92,7 +92,7 @@ size: (determinants.n_int*determinants.bit_kind/8,2,determinants.n_det) [psi_coef_qp_edit] interface: ezfio doc: Coefficients of the wave function -type: double precision +type: double precision size: (determinants.n_det_qp_edit,determinants.n_states) [psi_det_qp_edit] @@ -126,13 +126,13 @@ default: 1. [thresh_sym] type: Threshold -doc: Thresholds to check if a determinant is connected with HF +doc: Thresholds to check if a determinant is connected with HF interface: ezfio,provider,ocaml default: 1.e-15 [pseudo_sym] type: logical -doc: If |true|, discard any Slater determinants with an interaction smaller than thresh_sym with HF. +doc: If |true|, discard any Slater determinants with an interaction smaller than thresh_sym with HF. interface: ezfio,provider,ocaml default: False diff --git a/src/determinants/density_matrix.irp.f b/src/determinants/density_matrix.irp.f index 7c4a7fec..fb461848 100644 --- a/src/determinants/density_matrix.irp.f +++ b/src/determinants/density_matrix.irp.f @@ -304,7 +304,7 @@ BEGIN_PROVIDER [ double precision, c0_weight, (N_states) ] c0_weight(i) = c0_weight(i) * c enddo else - c0_weight = 1.d0 + c0_weight(:) = 1.d0 endif END_PROVIDER @@ -321,7 +321,7 @@ BEGIN_PROVIDER [ double precision, state_average_weight, (N_states) ] if (weight_one_e_dm == 0) then state_average_weight(:) = c0_weight(:) else if (weight_one_e_dm == 1) then - state_average_weight(:) = 1./N_states + state_average_weight(:) = 1.d0/N_states else call ezfio_has_determinants_state_average_weight(exists) if (exists) then @@ -384,6 +384,14 @@ END_PROVIDER END_PROVIDER +BEGIN_PROVIDER [ double precision, one_e_dm_ao, (ao_num, ao_num)] + implicit none + BEGIN_DOC + ! one_e_dm_ao = one_e_dm_ao_alpha + one_e_dm_ao_beta + END_DOC + one_e_dm_ao = one_e_dm_ao_alpha + one_e_dm_ao_beta +END_PROVIDER + subroutine get_occupation_from_dets(istate,occupation) implicit none diff --git a/src/determinants/s2.irp.f b/src/determinants/s2.irp.f index d73b2dbf..a80ba923 100644 --- a/src/determinants/s2.irp.f +++ b/src/determinants/s2.irp.f @@ -103,13 +103,17 @@ BEGIN_PROVIDER [ double precision, expected_s2] END_PROVIDER -BEGIN_PROVIDER [ double precision, s2_values, (N_states) ] + BEGIN_PROVIDER [ double precision, s2_values, (N_states) ] +&BEGIN_PROVIDER [ double precision, s_values, (N_states) ] implicit none BEGIN_DOC ! array of the averaged values of the S^2 operator on the various states END_DOC integer :: i call u_0_S2_u_0(s2_values,psi_coef,n_det,psi_det,N_int,N_states,psi_det_size) + do i = 1, N_states + s_values(i) = 0.5d0 *(-1.d0 + dsqrt(1.d0 + 4.d0 * s2_values(i))) + enddo END_PROVIDER diff --git a/src/determinants/slater_rules.irp.f b/src/determinants/slater_rules.irp.f index 04cf861f..a34608f9 100644 --- a/src/determinants/slater_rules.irp.f +++ b/src/determinants/slater_rules.irp.f @@ -438,7 +438,7 @@ subroutine bitstring_to_list_ab( string, list, n_elements, Nint) use bitmasks implicit none BEGIN_DOC - ! Gives the inidices(+1) of the bits set to 1 in the bit string + ! Gives the indices(+1) of the bits set to 1 in the bit string ! For alpha/beta determinants. END_DOC integer, intent(in) :: Nint @@ -652,7 +652,6 @@ subroutine i_H_j(key_i,key_j,Nint,hij) case (1) call get_single_excitation(key_i,key_j,exc,phase,Nint) !DIR$ FORCEINLINE - call bitstring_to_list_ab(key_i, occ, n_occ_ab, Nint) if (exc(0,1,1) == 1) then ! Single alpha m = exc(1,1,1) @@ -671,10 +670,6 @@ subroutine i_H_j(key_i,key_j,Nint,hij) end select end - - - - subroutine i_H_j_verbose(key_i,key_j,Nint,hij,hmono,hdouble,phase) use bitmasks implicit none @@ -1009,7 +1004,6 @@ subroutine i_H_psi(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array) end - subroutine i_H_psi_minilist(key,keys,idx_key,N_minilist,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array) use bitmasks implicit none diff --git a/src/determinants/slater_rules_wee_mono.irp.f b/src/determinants/slater_rules_wee_mono.irp.f index 4c1c9330..7c2ad148 100644 --- a/src/determinants/slater_rules_wee_mono.irp.f +++ b/src/determinants/slater_rules_wee_mono.irp.f @@ -282,9 +282,7 @@ subroutine i_H_j_two_e(key_i,key_j,Nint,hij) double precision :: get_two_e_integral integer :: m,n,p,q integer :: i,j,k - integer :: occ(Nint*bit_kind_size,2) double precision :: diag_H_mat_elem, phase,phase_2 - integer :: n_occ_ab(2) PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals ref_bitmask_two_e_energy ASSERT (Nint > 0) @@ -342,7 +340,6 @@ subroutine i_H_j_two_e(key_i,key_j,Nint,hij) case (1) call get_single_excitation(key_i,key_j,exc,phase,Nint) !DIR$ FORCEINLINE - call bitstring_to_list_ab(key_i, occ, n_occ_ab, Nint) if (exc(0,1,1) == 1) then ! Mono alpha m = exc(1,1,1) diff --git a/src/determinants/spindeterminants.irp.f b/src/determinants/spindeterminants.irp.f index dea4a566..dd55e112 100644 --- a/src/determinants/spindeterminants.irp.f +++ b/src/determinants/spindeterminants.irp.f @@ -585,7 +585,7 @@ END_PROVIDER enddo !$OMP ENDDO !$OMP END PARALLEL - call i8radix_sort(to_sort, psi_bilinear_matrix_transp_order, N_det,-1) + call i8sort(to_sort, psi_bilinear_matrix_transp_order, N_det) call iset_order(psi_bilinear_matrix_transp_rows,psi_bilinear_matrix_transp_order,N_det) call iset_order(psi_bilinear_matrix_transp_columns,psi_bilinear_matrix_transp_order,N_det) !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(l) diff --git a/src/dft_utils_in_r/ao_prod_mlti_pl.irp.f b/src/dft_utils_in_r/ao_prod_mlti_pl.irp.f new file mode 100644 index 00000000..1af34d74 --- /dev/null +++ b/src/dft_utils_in_r/ao_prod_mlti_pl.irp.f @@ -0,0 +1,116 @@ +BEGIN_PROVIDER [ double precision, ao_overlap_abs_grid, (ao_num, ao_num)] + implicit none + integer :: i,j,ipoint + double precision :: contrib, weight,r(3) + ao_overlap_abs_grid = 0.D0 + do ipoint = 1,n_points_final_grid + r(:) = final_grid_points(:,ipoint) + weight = final_weight_at_r_vector(ipoint) + do i = 1, ao_num + do j = 1, ao_num + contrib = dabs(aos_in_r_array(j,ipoint) * aos_in_r_array(i,ipoint)) * weight + ao_overlap_abs_grid(j,i) += contrib + enddo + enddo + enddo + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, ao_prod_center, (3, ao_num, ao_num)] + implicit none + BEGIN_DOC +! ao_prod_center(1:3,j,i) = \int dr |phi_i(r) phi_j(r)| x/y/z / \int |phi_i(r) phi_j(r)| +! +! if \int |phi_i(r) phi_j(r)| < 1.d-15 then ao_prod_center = 0. + END_DOC + integer :: i,j,m,ipoint + double precision :: contrib, weight,r(3) + ao_prod_center = 0.D0 + do ipoint = 1,n_points_final_grid + r(:) = final_grid_points(:,ipoint) + weight = final_weight_at_r_vector(ipoint) + do i = 1, ao_num + do j = 1, ao_num + contrib = dabs(aos_in_r_array(j,ipoint) * aos_in_r_array(i,ipoint)) * weight + do m = 1, 3 + ao_prod_center(m,j,i) += contrib * r(m) + enddo + enddo + enddo + enddo + do i = 1, ao_num + do j = 1, ao_num + if(dabs(ao_overlap_abs_grid(j,i)).gt.1.d-10)then + do m = 1, 3 + ao_prod_center(m,j,i) *= 1.d0/ao_overlap_abs_grid(j,i) + enddo + endif + enddo + enddo + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, ao_prod_sigma, (ao_num, ao_num)] + implicit none + BEGIN_DOC +! ao_prod_sigma(i,j) = \int |phi_i(r) phi_j(r)| dsqrt((x - <|i|x|j|>)^2 + (y - <|i|y|j|>)^2 +(z - <|i|z|j|>)^2) / \int |phi_i(r) phi_j(r)| +! +! gives you a precise idea of the spatial extension of the distribution phi_i(r) phi_j(r) + END_DOC + ao_prod_sigma = 0.d0 + integer :: i,j,m,ipoint + double precision :: contrib, weight,r(3),contrib_x2 + do ipoint = 1,n_points_final_grid + r(:) = final_grid_points(:,ipoint) + weight = final_weight_at_r_vector(ipoint) + do i = 1, ao_num + do j = 1, ao_num + contrib = dabs(aos_in_r_array(j,ipoint) * aos_in_r_array(i,ipoint)) * weight + contrib_x2 = 0.d0 + do m = 1, 3 + contrib_x2 += (r(m) - ao_prod_center(m,j,i)) * (r(m) - ao_prod_center(m,j,i)) + enddo + contrib_x2 = dsqrt(contrib_x2) + ao_prod_sigma(j,i) += contrib * contrib_x2 + enddo + enddo + enddo + + do i = 1, ao_num + do j = 1, ao_num + if(dabs(ao_overlap_abs_grid(j,i)).gt.1.d-10)then + ao_prod_sigma(j,i) *= 1.d0/ao_overlap_abs_grid(j,i) + endif + enddo + enddo + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, ao_prod_dist_grid, (ao_num, ao_num, n_points_final_grid)] + implicit none + BEGIN_DOC + ! ao_prod_dist_grid(j,i,ipoint) = distance between the center of |phi_i(r) phi_j(r)| and the grid point r(ipoint) + END_DOC + integer :: i,j,m,ipoint + double precision :: distance,r(3) + do ipoint = 1, n_points_final_grid + r(:) = final_grid_points(:,ipoint) + do i = 1, ao_num + do j = 1, ao_num + distance = 0.d0 + do m = 1, 3 + distance += (ao_prod_center(m,j,i) - r(m))*(ao_prod_center(m,j,i) - r(m)) + enddo + distance = dsqrt(distance) + ao_prod_dist_grid(j,i,ipoint) = distance + enddo + enddo + enddo + +END_PROVIDER + + +!BEGIN_PROVIDER [ double precision, ao_abs_prod_j1b, (ao_num, ao_num)] +! implicit none +! +!END_PROVIDER diff --git a/src/dressing/alpha_factory.irp.f b/src/dressing/alpha_factory.irp.f index 5eeeb1a6..c7adffe3 100644 --- a/src/dressing/alpha_factory.irp.f +++ b/src/dressing/alpha_factory.irp.f @@ -1179,7 +1179,7 @@ subroutine bitstring_to_list_in_selection( string, list, n_elements, Nint) use bitmasks implicit none BEGIN_DOC - ! Gives the inidices(+1) of the bits set to 1 in the bit string + ! Gives the indices(+1) of the bits set to 1 in the bit string END_DOC integer, intent(in) :: Nint integer(bit_kind), intent(in) :: string(Nint) diff --git a/src/dressing/run_dress_slave.irp.f b/src/dressing/run_dress_slave.irp.f index a33fb1dd..08b654c9 100644 --- a/src/dressing/run_dress_slave.irp.f +++ b/src/dressing/run_dress_slave.irp.f @@ -72,7 +72,7 @@ subroutine run_dress_slave(thread,iproce,energy) provide psi_energy ending = dress_N_cp+1 ntask_tbd = 0 - call omp_set_max_active_levels(8) + call set_multiple_levels_omp(.True.) !$OMP PARALLEL DEFAULT(SHARED) & !$OMP PRIVATE(interesting, breve_delta_m, task_id) & @@ -84,7 +84,7 @@ subroutine run_dress_slave(thread,iproce,energy) zmq_socket_push = new_zmq_push_socket(thread) integer, external :: connect_to_taskserver !$OMP CRITICAL - call omp_set_max_active_levels(1) + call set_multiple_levels_omp(.False.) if (connect_to_taskserver(zmq_to_qp_run_socket,worker_id,thread) == -1) then print *, irp_here, ': Unable to connect to task server' stop -1 @@ -296,7 +296,7 @@ subroutine run_dress_slave(thread,iproce,energy) !$OMP END CRITICAL !$OMP END PARALLEL - call omp_set_max_active_levels(1) + call set_multiple_levels_omp(.False.) ! do i=0,dress_N_cp+1 ! call omp_destroy_lock(lck_sto(i)) ! end do diff --git a/src/ezfio_files/output.irp.f b/src/ezfio_files/output.irp.f index 48512f92..7b2663a0 100644 --- a/src/ezfio_files/output.irp.f +++ b/src/ezfio_files/output.irp.f @@ -25,7 +25,7 @@ subroutine write_time(iunit) ct = ct - output_cpu_time_0 call wall_time(wt) wt = wt - output_wall_time_0 - write(6,'(A,F14.6,A,F14.6,A)') & + write(6,'(A,F14.2,A,F14.2,A)') & '.. >>>>> [ WALL TIME: ', wt, ' s ] [ CPU TIME: ', ct, ' s ] <<<<< ..' write(6,*) end diff --git a/src/iterations/print_extrapolation.irp.f b/src/iterations/print_extrapolation.irp.f index cb46fb67..7c6dbb9b 100644 --- a/src/iterations/print_extrapolation.irp.f +++ b/src/iterations/print_extrapolation.irp.f @@ -35,12 +35,13 @@ subroutine print_extrapolated_energy do k=2,min(N_iter,8) write(*,'(F11.4,X,3(X,F18.8))') pt2_iterations(i,N_iter+1-k), extrapolated_energy(k,i), & extrapolated_energy(k,i) - extrapolated_energy(k,1), & - (extrapolated_energy(k,i) - extrapolated_energy(k,1) ) * 27.211396641308d0 + (extrapolated_energy(k,i) - extrapolated_energy(k,1) ) * ha_to_ev enddo write(*,*) '=========== ', '=================== ', '=================== ', '===================' enddo print *, '' + call ezfio_set_fci_energy_extrapolated(extrapolated_energy(min(N_iter,3),1:N_states)) end subroutine diff --git a/src/iterations/print_summary.irp.f b/src/iterations/print_summary.irp.f index 641ee209..a0db3534 100644 --- a/src/iterations/print_summary.irp.f +++ b/src/iterations/print_summary.irp.f @@ -36,7 +36,7 @@ subroutine print_summary(e_,pt2_data,pt2_data_err,n_det_,n_configuration_,n_st,s write(*,fmt) '# E ', e_(1:N_states_p) if (N_states_p > 1) then write(*,fmt) '# Excit. (au)', e_(1:N_states_p)-e_(1) - write(*,fmt) '# Excit. (eV)', (e_(1:N_states_p)-e_(1))*27.211396641308d0 + write(*,fmt) '# Excit. (eV)', (e_(1:N_states_p)-e_(1))*ha_to_ev endif write(fmt,*) '(A13,', 2*N_states_p, '(1X,F14.8))' write(*,fmt) '# PT2 '//pt2_string, (pt2_data % pt2(k), pt2_data_err % pt2(k), k=1,N_states_p) @@ -47,8 +47,8 @@ subroutine print_summary(e_,pt2_data,pt2_data_err,n_det_,n_configuration_,n_st,s if (N_states_p > 1) then write(*,fmt) '# Excit. (au)', ( (e_(k)+pt2_data % pt2(k)-e_(1)-pt2_data % pt2(1)), & dsqrt(pt2_data_err % pt2(k)*pt2_data_err % pt2(k)+pt2_data_err % pt2(1)*pt2_data_err % pt2(1)), k=1,N_states_p) - write(*,fmt) '# Excit. (eV)', ( (e_(k)+pt2_data % pt2(k)-e_(1)-pt2_data % pt2(1))*27.211396641308d0, & - dsqrt(pt2_data_err % pt2(k)*pt2_data_err % pt2(k)+pt2_data_err % pt2(1)*pt2_data_err % pt2(1))*27.211396641308d0, k=1,N_states_p) + write(*,fmt) '# Excit. (eV)', ( (e_(k)+pt2_data % pt2(k)-e_(1)-pt2_data % pt2(1))*ha_to_ev, & + dsqrt(pt2_data_err % pt2(k)*pt2_data_err % pt2(k)+pt2_data_err % pt2(1)*pt2_data_err % pt2(1))*ha_to_ev, k=1,N_states_p) endif write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))' write(*,fmt) @@ -82,23 +82,23 @@ subroutine print_summary(e_,pt2_data,pt2_data_err,n_det_,n_configuration_,n_st,s print *, 'Variational Energy difference (au | eV)' do i=2, N_states_p print*,'Delta E = ', (e_(i) - e_(1)), & - (e_(i) - e_(1)) * 27.211396641308d0 + (e_(i) - e_(1)) * ha_to_ev enddo print *, '-----' print*, 'Variational + perturbative Energy difference (au | eV)' do i=2, N_states_p print*,'Delta E = ', (e_(i)+ pt2_data % pt2(i) - (e_(1) + pt2_data % pt2(1))), & - (e_(i)+ pt2_data % pt2(i) - (e_(1) + pt2_data % pt2(1))) * 27.211396641308d0 + (e_(i)+ pt2_data % pt2(i) - (e_(1) + pt2_data % pt2(1))) * ha_to_ev enddo print *, '-----' print*, 'Variational + renormalized perturbative Energy difference (au | eV)' do i=2, N_states_p print*,'Delta E = ', (e_(i)+ pt2_data % rpt2(i) - (e_(1) + pt2_data % rpt2(1))), & - (e_(i)+ pt2_data % rpt2(i) - (e_(1) + pt2_data % rpt2(1))) * 27.211396641308d0 + (e_(i)+ pt2_data % rpt2(i) - (e_(1) + pt2_data % rpt2(1))) * ha_to_ev enddo endif - call print_energy_components() +! call print_energy_components() end subroutine diff --git a/src/mo_basis/mos_in_r.irp.f b/src/mo_basis/mos_in_r.irp.f index ee2795d0..e5d3b243 100644 --- a/src/mo_basis/mos_in_r.irp.f +++ b/src/mo_basis/mos_in_r.irp.f @@ -1,6 +1,9 @@ subroutine give_all_mos_at_r(r,mos_array) implicit none + BEGIN_DOC +! mos_array(i) = ith MO function evaluated at "r" + END_DOC double precision, intent(in) :: r(3) double precision, intent(out) :: mos_array(mo_num) double precision :: aos_array(ao_num) diff --git a/src/mo_two_e_erf_ints/map_integrals_erf.irp.f b/src/mo_two_e_erf_ints/map_integrals_erf.irp.f index 73050ec5..3405ec2b 100644 --- a/src/mo_two_e_erf_ints/map_integrals_erf.irp.f +++ b/src/mo_two_e_erf_ints/map_integrals_erf.irp.f @@ -235,11 +235,11 @@ subroutine get_mo_two_e_integrals_erf_ij(k,l,sze,out_array,map) logical :: integral_is_in_map if (key_kind == 8) then - call i8radix_sort(hash,iorder,kk,-1) + call i8sort(hash,iorder,kk) else if (key_kind == 4) then - call iradix_sort(hash,iorder,kk,-1) + call isort(hash,iorder,kk) else if (key_kind == 2) then - call i2radix_sort(hash,iorder,kk,-1) + call i2sort(hash,iorder,kk) endif call map_get_many(mo_integrals_erf_map, hash, tmp_val, kk) @@ -290,11 +290,11 @@ subroutine get_mo_two_e_integrals_erf_i1j1(k,l,sze,out_array,map) logical :: integral_is_in_map if (key_kind == 8) then - call i8radix_sort(hash,iorder,kk,-1) + call i8sort(hash,iorder,kk) else if (key_kind == 4) then - call iradix_sort(hash,iorder,kk,-1) + call isort(hash,iorder,kk) else if (key_kind == 2) then - call i2radix_sort(hash,iorder,kk,-1) + call i2sort(hash,iorder,kk) endif call map_get_many(mo_integrals_erf_map, hash, tmp_val, kk) diff --git a/src/mo_two_e_ints/mo_bi_integrals.irp.f b/src/mo_two_e_ints/mo_bi_integrals.irp.f index d58932ce..56d8cf28 100644 --- a/src/mo_two_e_ints/mo_bi_integrals.irp.f +++ b/src/mo_two_e_ints/mo_bi_integrals.irp.f @@ -38,7 +38,7 @@ BEGIN_PROVIDER [ logical, mo_two_e_integrals_in_map ] print*, 'MO integrals provided' return else - PROVIDE ao_two_e_integrals_in_map + PROVIDE ao_two_e_integrals_in_map endif print *, '' @@ -53,7 +53,11 @@ BEGIN_PROVIDER [ logical, mo_two_e_integrals_in_map ] ! call four_idx_novvvv call four_idx_novvvv_old else - call add_integrals_to_map(full_ijkl_bitmask_4) + if (dble(ao_num)**4 * 32.d-9 < dble(qp_max_mem)) then + call four_idx_dgemm + else + call add_integrals_to_map(full_ijkl_bitmask_4) + endif endif call wall_time(wall_2) @@ -77,6 +81,94 @@ BEGIN_PROVIDER [ logical, mo_two_e_integrals_in_map ] END_PROVIDER +subroutine four_idx_dgemm + implicit none + integer :: p,q,r,s,i,j,k,l + double precision, allocatable :: a1(:,:,:,:) + double precision, allocatable :: a2(:,:,:,:) + + allocate (a1(ao_num,ao_num,ao_num,ao_num)) + + print *, 'Getting AOs' + !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(q,r,s) + do s=1,ao_num + do r=1,ao_num + do q=1,ao_num + call get_ao_two_e_integrals(q,r,s,ao_num,a1(1,q,r,s)) + enddo + enddo + enddo + !$OMP END PARALLEL DO + + print *, '1st transformation' + ! 1st transformation + allocate (a2(ao_num,ao_num,ao_num,mo_num)) + call dgemm('T','N', (ao_num*ao_num*ao_num), mo_num, ao_num, 1.d0, a1, ao_num, mo_coef, ao_num, 0.d0, a2, (ao_num*ao_num*ao_num)) + + ! 2nd transformation + print *, '2nd transformation' + deallocate (a1) + allocate (a1(ao_num,ao_num,mo_num,mo_num)) + call dgemm('T','N', (ao_num*ao_num*mo_num), mo_num, ao_num, 1.d0, a2, ao_num, mo_coef, ao_num, 0.d0, a1, (ao_num*ao_num*mo_num)) + + ! 3rd transformation + print *, '3rd transformation' + deallocate (a2) + allocate (a2(ao_num,mo_num,mo_num,mo_num)) + call dgemm('T','N', (ao_num*mo_num*mo_num), mo_num, ao_num, 1.d0, a1, ao_num, mo_coef, ao_num, 0.d0, a2, (ao_num*mo_num*mo_num)) + + ! 4th transformation + print *, '4th transformation' + deallocate (a1) + allocate (a1(mo_num,mo_num,mo_num,mo_num)) + call dgemm('T','N', (mo_num*mo_num*mo_num), mo_num, ao_num, 1.d0, a2, ao_num, mo_coef, ao_num, 0.d0, a1, (mo_num*mo_num*mo_num)) + + deallocate (a2) + + integer :: n_integrals, size_buffer + integer(key_kind) , allocatable :: buffer_i(:) + real(integral_kind), allocatable :: buffer_value(:) + size_buffer = min(ao_num*ao_num*ao_num,16000000) + + !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,k,l,buffer_value,buffer_i,n_integrals) + allocate ( buffer_i(size_buffer), buffer_value(size_buffer) ) + + n_integrals = 0 + !$OMP DO + do l=1,mo_num + do k=1,mo_num + do j=1,l + do i=1,k + if (abs(a1(i,j,k,l)) < mo_integrals_threshold) then + cycle + endif + n_integrals += 1 + buffer_value(n_integrals) = a1(i,j,k,l) + !DIR$ FORCEINLINE + call mo_two_e_integrals_index(i,j,k,l,buffer_i(n_integrals)) + if (n_integrals == size_buffer) then + call map_append(mo_integrals_map, buffer_i, buffer_value, n_integrals) + n_integrals = 0 + endif + enddo + enddo + enddo + enddo + !$OMP END DO + + call map_append(mo_integrals_map, buffer_i, buffer_value, n_integrals) + + deallocate(buffer_i, buffer_value) + !$OMP END PARALLEL + + deallocate (a1) + + call map_unique(mo_integrals_map) + + integer*8 :: get_mo_map_size, mo_map_size + mo_map_size = get_mo_map_size() + +end subroutine subroutine add_integrals_to_map(mask_ijkl) use bitmasks @@ -153,14 +245,14 @@ subroutine add_integrals_to_map(mask_ijkl) return endif - size_buffer = min(ao_num*ao_num*ao_num,16000000) - print*, 'Buffers : ', 8.*(mo_num*(n_j)*(n_k+1) + mo_num+& - ao_num+ao_num*ao_num+ size_buffer*3)/(1024*1024), 'MB / core' - double precision :: accu_bis accu_bis = 0.d0 call wall_time(wall_1) + size_buffer = min( (qp_max_mem/(nproc*5)),mo_num*mo_num*mo_num) + print*, 'Buffers : ', 8.*(mo_num*(n_j)*(n_k+1) + mo_num+& + ao_num+ao_num*ao_num+ size_buffer*3)/(1024*1024), 'MB / core' + !$OMP PARALLEL PRIVATE(l1,k1,j1,i1,i2,i3,i4,i,j,k,l,c, ii1,kmax, & !$OMP two_e_tmp_0_idx, two_e_tmp_0, two_e_tmp_1,two_e_tmp_2,two_e_tmp_3,& !$OMP buffer_i,buffer_value,n_integrals,wall_2,i0,j0,k0,l0, & @@ -171,6 +263,10 @@ subroutine add_integrals_to_map(mask_ijkl) !$OMP mo_coef_transp_is_built, list_ijkl, & !$OMP mo_coef_is_built, wall_1, & !$OMP mo_coef,mo_integrals_threshold,mo_integrals_map) + + thread_num = 0 + !$ thread_num = omp_get_thread_num() + n_integrals = 0 wall_0 = wall_1 allocate(two_e_tmp_3(mo_num, n_j, n_k), & @@ -181,8 +277,6 @@ subroutine add_integrals_to_map(mask_ijkl) buffer_i(size_buffer), & buffer_value(size_buffer) ) - thread_num = 0 - !$ thread_num = omp_get_thread_num() !$OMP DO SCHEDULE(guided) do l1 = 1,ao_num two_e_tmp_3 = 0.d0 diff --git a/src/non_h_ints_mu/j12_nucl_utils.irp.f b/src/non_h_ints_mu/j12_nucl_utils.irp.f index f3b68f43..a515e0b8 100644 --- a/src/non_h_ints_mu/j12_nucl_utils.irp.f +++ b/src/non_h_ints_mu/j12_nucl_utils.irp.f @@ -237,6 +237,23 @@ end function j12_mu ! --- +double precision function j12_mu_r12(r12) + + include 'constants.include.F' + + implicit none + double precision, intent(in) :: r12 + double precision :: mu_r12 + + mu_r12 = mu_erf * r12 + + j12_mu_r12 = 0.5d0 * r12 * (1.d0 - derf(mu_r12)) - inv_sq_pi_2 * dexp(-mu_r12*mu_r12) / mu_erf + + return +end function j12_mu_r12 + +! --- + double precision function j12_mu_gauss(r1, r2) implicit none diff --git a/src/non_hermit_dav/biorthog.irp.f b/src/non_hermit_dav/biorthog.irp.f index 89e5b4f4..b621206a 100644 --- a/src/non_hermit_dav/biorthog.irp.f +++ b/src/non_hermit_dav/biorthog.irp.f @@ -315,7 +315,7 @@ subroutine non_hrmt_bieig(n, A, thr_d, thr_nd, leigvec, reigvec, n_real_eigv, ei ! write(*, '(1000(F16.10,X))') VL(:,i) !enddo - thr_diag = 1d-07 + thr_diag = 1d-06 thr_norm = 1d+10 call check_EIGVEC(n, n, A, WR, VL, VR, thr_diag, thr_norm, .false.) @@ -337,19 +337,21 @@ subroutine non_hrmt_bieig(n, A, thr_d, thr_nd, leigvec, reigvec, n_real_eigv, ei else print*, 'Found an imaginary component to eigenvalue on i = ', i print*, 'Re(i) + Im(i)', WR(i), WI(i) - stop endif enddo + if(n_good.ne.n)then + print*,'there are some imaginary eigenvalues ' + thr_diag = 1d-03 + n_good = n + endif allocate(list_good(n_good), iorder(n_good)) n_good = 0 do i = 1, n - if( dabs(WI(i)).lt.thr ) then - n_good += 1 - list_good(n_good) = i - eigval(n_good) = WR(i) - endif + n_good += 1 + list_good(n_good) = i + eigval(n_good) = WR(i) enddo deallocate( WR, WI ) diff --git a/src/scf_utils/diagonalize_fock.irp.f b/src/scf_utils/diagonalize_fock.irp.f index da1d44a7..a6f19c05 100644 --- a/src/scf_utils/diagonalize_fock.irp.f +++ b/src/scf_utils/diagonalize_fock.irp.f @@ -72,11 +72,6 @@ BEGIN_PROVIDER [ double precision, eigenvectors_Fock_matrix_mo, (ao_num,mo_num) liwork = -1 F_save = F - !print *, ' Fock matrix' - !do i = 1, mo_num - ! write(*, '(1000(F16.10,X))') F_save(:,i) - !enddo - call dsyevd( 'V', 'U', mo_num, F, & size(F,1), diag, work, lwork, iwork, liwork, info) @@ -107,16 +102,6 @@ BEGIN_PROVIDER [ double precision, eigenvectors_Fock_matrix_mo, (ao_num,mo_num) endif endif - !print *, ' eigenvalues' - !do i = 1, mo_num - ! write(*, '(1000(F16.10,X))') diag(i) - !enddo - !print *, ' eigenvectors' - !do i = 1, mo_num - ! write(*, '(1000(F16.10,X))') F(:,i) - !enddo - - call dgemm('N','N',ao_num,mo_num,mo_num, 1.d0, & mo_coef, size(mo_coef,1), F, size(F,1), & 0.d0, eigenvectors_Fock_matrix_mo, size(eigenvectors_Fock_matrix_mo,1)) diff --git a/src/scf_utils/fock_matrix.irp.f b/src/scf_utils/fock_matrix.irp.f index 61633d3b..9a95caa1 100644 --- a/src/scf_utils/fock_matrix.irp.f +++ b/src/scf_utils/fock_matrix.irp.f @@ -5,6 +5,90 @@ ! Fock matrix on the MO basis. ! For open shells, the ROHF Fock Matrix is :: ! + ! | Rcc | F^b | Fcv | + ! |-----------------------| + ! | F^b | Roo | F^a | + ! |-----------------------| + ! | Fcv | F^a | Rvv | + ! + ! C: Core, O: Open, V: Virtual + ! + ! Rcc = Acc Fcc^a + Bcc Fcc^b + ! Roo = Aoo Foo^a + Boo Foo^b + ! Rcc = Avv Fvv^a + Bvv Fvv^b + ! Fcv = (F^a + F^b)/2 + ! + ! F^a: Fock matrix alpha (MO), F^b: Fock matrix beta (MO) + ! A,B: Coupling parameters + ! + ! J. Chem. Phys. 133, 141102 (2010), https://doi.org/10.1063/1.3503173 + ! Coupling parameters from J. Chem. Phys. 125, 204110 (2006); https://doi.org/10.1063/1.2393223. + ! cc oo vv + ! A -0.5 0.5 1.5 + ! B 1.5 0.5 -0.5 + ! + END_DOC + integer :: i,j,n + if (elec_alpha_num == elec_beta_num) then + Fock_matrix_mo = Fock_matrix_mo_alpha + else + ! Core + do j = 1, elec_beta_num + ! Core + do i = 1, elec_beta_num + fock_matrix_mo(i,j) = - 0.5d0 * fock_matrix_mo_alpha(i,j) & + + 1.5d0 * fock_matrix_mo_beta(i,j) + enddo + ! Open + do i = elec_beta_num+1, elec_alpha_num + fock_matrix_mo(i,j) = fock_matrix_mo_beta(i,j) + enddo + ! Virtual + do i = elec_alpha_num+1, mo_num + fock_matrix_mo(i,j) = 0.5d0 * fock_matrix_mo_alpha(i,j) & + + 0.5d0 * fock_matrix_mo_beta(i,j) + enddo + enddo + ! Open + do j = elec_beta_num+1, elec_alpha_num + ! Core + do i = 1, elec_beta_num + fock_matrix_mo(i,j) = fock_matrix_mo_beta(i,j) + enddo + ! Open + do i = elec_beta_num+1, elec_alpha_num + fock_matrix_mo(i,j) = 0.5d0 * fock_matrix_mo_alpha(i,j) & + + 0.5d0 * fock_matrix_mo_beta(i,j) + enddo + ! Virtual + do i = elec_alpha_num+1, mo_num + fock_matrix_mo(i,j) = fock_matrix_mo_alpha(i,j) + enddo + enddo + ! Virtual + do j = elec_alpha_num+1, mo_num + ! Core + do i = 1, elec_beta_num + fock_matrix_mo(i,j) = 0.5d0 * fock_matrix_mo_alpha(i,j) & + + 0.5d0 * fock_matrix_mo_beta(i,j) + enddo + ! Open + do i = elec_beta_num+1, elec_alpha_num + fock_matrix_mo(i,j) = fock_matrix_mo_alpha(i,j) + enddo + ! Virtual + do i = elec_alpha_num+1, mo_num + fock_matrix_mo(i,j) = 1.5d0 * fock_matrix_mo_alpha(i,j) & + - 0.5d0 * fock_matrix_mo_beta(i,j) + enddo + enddo + endif + + ! Old + ! BEGIN_DOC + ! Fock matrix on the MO basis. + ! For open shells, the ROHF Fock Matrix is :: + ! ! | F-K | F + K/2 | F | ! |---------------------------------| ! | F + K/2 | F | F - K/2 | @@ -16,64 +100,64 @@ ! ! K = Fb - Fa ! - END_DOC - integer :: i,j,n - if (elec_alpha_num == elec_beta_num) then - Fock_matrix_mo = Fock_matrix_mo_alpha - else + ! END_DOC + !integer :: i,j,n + !if (elec_alpha_num == elec_beta_num) then + ! Fock_matrix_mo = Fock_matrix_mo_alpha + !else - do j=1,elec_beta_num - ! F-K - do i=1,elec_beta_num !CC - Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))& - - (Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j)) - enddo - ! F+K/2 - do i=elec_beta_num+1,elec_alpha_num !CA - Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))& - + 0.5d0*(Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j)) - enddo - ! F - do i=elec_alpha_num+1, mo_num !CV - Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j)) - enddo - enddo + ! do j=1,elec_beta_num + ! ! F-K + ! do i=1,elec_beta_num !CC + ! Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))& + ! - (Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j)) + ! enddo + ! ! F+K/2 + ! do i=elec_beta_num+1,elec_alpha_num !CA + ! Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))& + ! + 0.5d0*(Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j)) + ! enddo + ! ! F + ! do i=elec_alpha_num+1, mo_num !CV + ! Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j)) + ! enddo + ! enddo - do j=elec_beta_num+1,elec_alpha_num - ! F+K/2 - do i=1,elec_beta_num !AC - Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))& - + 0.5d0*(Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j)) - enddo - ! F - do i=elec_beta_num+1,elec_alpha_num !AA - Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j)) - enddo - ! F-K/2 - do i=elec_alpha_num+1, mo_num !AV - Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))& - - 0.5d0*(Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j)) - enddo - enddo + ! do j=elec_beta_num+1,elec_alpha_num + ! ! F+K/2 + ! do i=1,elec_beta_num !AC + ! Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))& + ! + 0.5d0*(Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j)) + ! enddo + ! ! F + ! do i=elec_beta_num+1,elec_alpha_num !AA + ! Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j)) + ! enddo + ! ! F-K/2 + ! do i=elec_alpha_num+1, mo_num !AV + ! Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))& + ! - 0.5d0*(Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j)) + ! enddo + ! enddo - do j=elec_alpha_num+1, mo_num - ! F - do i=1,elec_beta_num !VC - Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j)) - enddo - ! F-K/2 - do i=elec_beta_num+1,elec_alpha_num !VA - Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))& - - 0.5d0*(Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j)) - enddo - ! F+K - do i=elec_alpha_num+1,mo_num !VV - Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j)) & - + (Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j)) - enddo - enddo + ! do j=elec_alpha_num+1, mo_num + ! ! F + ! do i=1,elec_beta_num !VC + ! Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j)) + ! enddo + ! ! F-K/2 + ! do i=elec_beta_num+1,elec_alpha_num !VA + ! Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))& + ! - 0.5d0*(Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j)) + ! enddo + ! ! F+K + ! do i=elec_alpha_num+1,mo_num !VV + ! Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j)) & + ! + (Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j)) + ! enddo + ! enddo - endif + !endif do i = 1, mo_num Fock_matrix_diag_mo(i) = Fock_matrix_mo(i,i) diff --git a/src/scf_utils/scf_density_matrix_ao.irp.f b/src/scf_utils/scf_density_matrix_ao.irp.f index d3e9f196..7759431d 100644 --- a/src/scf_utils/scf_density_matrix_ao.irp.f +++ b/src/scf_utils/scf_density_matrix_ao.irp.f @@ -3,24 +3,15 @@ BEGIN_PROVIDER [double precision, SCF_density_matrix_ao_alpha, (ao_num,ao_num) ] BEGIN_DOC ! $C.C^t$ over $\alpha$ MOs END_DOC - SCF_density_matrix_ao_alpha = 0.d0 - if(elec_alpha_num.gt.0)then - call dgemm('N','T',ao_num,ao_num,elec_alpha_num,1.d0, & - mo_coef, size(mo_coef,1), & - mo_coef, size(mo_coef,1), 0.d0, & - SCF_density_matrix_ao_alpha, size(SCF_density_matrix_ao_alpha,1)) + if(elec_alpha_num > 0)then + call dgemm('N','T',ao_num,ao_num,elec_alpha_num,1.d0, & + mo_coef, size(mo_coef,1), & + mo_coef, size(mo_coef,1), 0.d0, & + SCF_density_matrix_ao_alpha, size(SCF_density_matrix_ao_alpha,1)) + else + SCF_density_matrix_ao_alpha = 0.d0 endif -! integer :: i, j -! double precision :: trace_density -! trace_density = 0.d0 -! do i = 1, ao_num !elec_alpha_num -! do j = 1, ao_num !elec_alpha_num -! trace_density = trace_density & -! + SCF_density_matrix_ao_alpha(j,i) * ao_overlap(j,i) -! enddo -! enddo -! print *, ' trace of SCF_density_matrix_ao_alpha =', trace_density END_PROVIDER @@ -29,12 +20,13 @@ BEGIN_PROVIDER [ double precision, SCF_density_matrix_ao_beta, (ao_num,ao_num) BEGIN_DOC ! $C.C^t$ over $\beta$ MOs END_DOC - SCF_density_matrix_ao_beta = 0.d0 - if(elec_beta_num.gt.0)then - call dgemm('N','T',ao_num,ao_num,elec_beta_num,1.d0, & - mo_coef, size(mo_coef,1), & - mo_coef, size(mo_coef,1), 0.d0, & - SCF_density_matrix_ao_beta, size(SCF_density_matrix_ao_beta,1)) + if(elec_beta_num > 0)then + call dgemm('N','T',ao_num,ao_num,elec_beta_num,1.d0, & + mo_coef, size(mo_coef,1), & + mo_coef, size(mo_coef,1), 0.d0, & + SCF_density_matrix_ao_beta, size(SCF_density_matrix_ao_beta,1)) + else + SCF_density_matrix_ao_beta = 0.d0 endif END_PROVIDER diff --git a/src/tc_bi_ortho/test_tc_fock.irp.f b/src/tc_bi_ortho/test_tc_fock.irp.f index a49a5958..26446daf 100644 --- a/src/tc_bi_ortho/test_tc_fock.irp.f +++ b/src/tc_bi_ortho/test_tc_fock.irp.f @@ -32,7 +32,7 @@ subroutine test_3e print*,'htot = ',htot print*,'' print*,'' - print*,'TC_one= ',TC_HF_one_electron_energy + print*,'TC_one= ',tc_hf_one_e_energy print*,'TC_two= ',TC_HF_two_e_energy print*,'TC_3e = ',diag_three_elem_hf print*,'TC_tot= ',TC_HF_energy diff --git a/src/tc_keywords/EZFIO.cfg b/src/tc_keywords/EZFIO.cfg index 9df73154..a206dfa9 100644 --- a/src/tc_keywords/EZFIO.cfg +++ b/src/tc_keywords/EZFIO.cfg @@ -78,7 +78,7 @@ default: True [symetric_fock_tc] type: logical -doc: If |true|, using F+F^\dagger as Fock TC +doc: If |true|, using F+F^t as Fock TC interface: ezfio,provider,ocaml default: False diff --git a/src/tc_scf/test_int.irp.f b/src/tc_scf/test_int.irp.f new file mode 100644 index 00000000..6961d2f0 --- /dev/null +++ b/src/tc_scf/test_int.irp.f @@ -0,0 +1,346 @@ +program test_ints + + BEGIN_DOC +! TODO : Put the documentation of the program here + END_DOC + + implicit none + + print *, 'starting ...' + + my_grid_becke = .True. +! my_n_pt_r_grid = 30 +! my_n_pt_a_grid = 50 + my_n_pt_r_grid = 10 ! small grid for quick debug + my_n_pt_a_grid = 26 ! small grid for quick debug + touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid + call routine_int2_u_grad1u_j1b2 + call routine_v_ij_erf_rk_cst_mu_j1b + call routine_x_v_ij_erf_rk_cst_mu_tmp_j1b + call routine_v_ij_u_cst_mu_j1b + +! +! call routine_test_j1b + +! call routine_int2_grad1u2_grad2u2_j1b2 +end + +subroutine routine_test_j1b + implicit none + integer :: i,icount,j + icount = 0 + do i = 1, List_all_comb_b3_size + if(dabs(List_all_comb_b3_coef(i)).gt.1.d-10)then + print*,'' + print*,List_all_comb_b3_expo(i),List_all_comb_b3_coef(i) + print*,List_all_comb_b3_cent(1:3,i) + print*,'' + icount += 1 + endif + + enddo + print*,'List_all_comb_b3_coef,icount = ',List_all_comb_b3_size,icount + do i = 1, ao_num + do j = 1, ao_num + do icount = 1, List_comb_b3_size_thr(j,i) + print*,'',j,i + print*,List_comb_thr_b3_expo(icount,j,i),List_comb_thr_b3_coef(icount,j,i) + print*,List_comb_thr_b3_cent(1:3,icount,j,i) + print*,'' + enddo +! enddo + enddo + enddo + print*,'max_List_comb_b3_size_thr = ',max_List_comb_b3_size_thr,List_all_comb_b3_size + +end + +subroutine routine_int2_u_grad1u_j1b2 + implicit none + integer :: i,j,ipoint,k,l + double precision :: weight,accu_relat, accu_abs, contrib + double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:) + + allocate(array(ao_num, ao_num, ao_num, ao_num)) + array = 0.d0 + allocate(array_ref(ao_num, ao_num, ao_num, ao_num)) + array_ref = 0.d0 + do ipoint = 1, n_points_final_grid + weight = final_weight_at_r_vector(ipoint) + do k = 1, ao_num + do l = 1, ao_num + do i = 1, ao_num + do j = 1, ao_num + array(j,i,l,k) += int2_u_grad1u_j1b2_test(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight + array_ref(j,i,l,k) += int2_u_grad1u_j1b2(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight + enddo + enddo + enddo + enddo + enddo + accu_relat = 0.d0 + accu_abs = 0.d0 + do k = 1, ao_num + do l = 1, ao_num + do i = 1, ao_num + do j = 1, ao_num + contrib = dabs(array(j,i,l,k) - array_ref(j,i,l,k)) + accu_abs += contrib + if(dabs(array_ref(j,i,l,k)).gt.1.d-10)then + accu_relat += contrib/dabs(array_ref(j,i,l,k)) + endif + enddo + enddo + enddo + enddo + print*,'accu_abs = ',accu_abs/dble(ao_num)**4 + print*,'accu_relat = ',accu_relat/dble(ao_num)**4 + + + +end + +subroutine routine_v_ij_erf_rk_cst_mu_j1b + implicit none + integer :: i,j,ipoint,k,l + double precision :: weight,accu_relat, accu_abs, contrib + double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:) +! print*,'ao_overlap_abs = ' +! do i = 1, ao_num +! write(*,'(100(F10.5,X))')ao_overlap_abs(i,:) +! enddo +! print*,'center = ' +! do i = 1, ao_num +! write(*,'(100(F10.5,X))')ao_prod_center(2,i,:) +! enddo +! print*,'sigma = ' +! do i = 1, ao_num +! write(*,'(100(F10.5,X))')ao_prod_sigma(i,:) +! enddo + + + allocate(array(ao_num, ao_num, ao_num, ao_num)) + array = 0.d0 + allocate(array_ref(ao_num, ao_num, ao_num, ao_num)) + array_ref = 0.d0 + do ipoint = 1, n_points_final_grid + weight = final_weight_at_r_vector(ipoint) + do k = 1, ao_num + do l = 1, ao_num + do i = 1, ao_num + do j = 1, ao_num + array(j,i,l,k) += v_ij_erf_rk_cst_mu_j1b_test(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight + array_ref(j,i,l,k) += v_ij_erf_rk_cst_mu_j1b(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight + enddo + enddo + enddo + enddo + enddo + accu_relat = 0.d0 + accu_abs = 0.d0 + do k = 1, ao_num + do l = 1, ao_num + do i = 1, ao_num + do j = 1, ao_num + contrib = dabs(array(j,i,l,k) - array_ref(j,i,l,k)) + accu_abs += contrib + if(dabs(array_ref(j,i,l,k)).gt.1.d-10)then + accu_relat += contrib/dabs(array_ref(j,i,l,k)) + endif + enddo + enddo + enddo + enddo + print*,'accu_abs = ',accu_abs/dble(ao_num)**4 + print*,'accu_relat = ',accu_relat/dble(ao_num)**4 + + + +end + + +subroutine routine_x_v_ij_erf_rk_cst_mu_tmp_j1b + implicit none + integer :: i,j,ipoint,k,l,m + double precision :: weight,accu_relat, accu_abs, contrib + double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:) +! print*,'ao_overlap_abs = ' +! do i = 1, ao_num +! write(*,'(100(F10.5,X))')ao_overlap_abs(i,:) +! enddo +! print*,'center = ' +! do i = 1, ao_num +! write(*,'(100(F10.5,X))')ao_prod_center(2,i,:) +! enddo +! print*,'sigma = ' +! do i = 1, ao_num +! write(*,'(100(F10.5,X))')ao_prod_sigma(i,:) +! enddo + + + allocate(array(ao_num, ao_num, ao_num, ao_num)) + array = 0.d0 + allocate(array_ref(ao_num, ao_num, ao_num, ao_num)) + array_ref = 0.d0 + do ipoint = 1, n_points_final_grid + weight = final_weight_at_r_vector(ipoint) + do k = 1, ao_num + do l = 1, ao_num + do i = 1, ao_num + do j = 1, ao_num + do m = 1, 3 + array(j,i,l,k) += x_v_ij_erf_rk_cst_mu_tmp_j1b_test(m,j,i,ipoint) * aos_grad_in_r_array_transp(m,k,ipoint) * aos_in_r_array(l,ipoint) * weight + array_ref(j,i,l,k) += x_v_ij_erf_rk_cst_mu_tmp_j1b(m,j,i,ipoint) * aos_grad_in_r_array_transp(m,k,ipoint) * aos_in_r_array(l,ipoint) * weight + enddo + enddo + enddo + enddo + enddo + enddo + accu_relat = 0.d0 + accu_abs = 0.d0 + do k = 1, ao_num + do l = 1, ao_num + do i = 1, ao_num + do j = 1, ao_num + contrib = dabs(array(j,i,l,k) - array_ref(j,i,l,k)) + accu_abs += contrib + if(dabs(array_ref(j,i,l,k)).gt.1.d-10)then + accu_relat += contrib/dabs(array_ref(j,i,l,k)) + endif + enddo + enddo + enddo + enddo + print*,'accu_abs = ',accu_abs/dble(ao_num)**4 + print*,'accu_relat = ',accu_relat/dble(ao_num)**4 + + + +end + + + +subroutine routine_v_ij_u_cst_mu_j1b + implicit none + integer :: i,j,ipoint,k,l + double precision :: weight,accu_relat, accu_abs, contrib + double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:) +! print*,'ao_overlap_abs = ' +! do i = 1, ao_num +! write(*,'(100(F10.5,X))')ao_overlap_abs(i,:) +! enddo +! print*,'center = ' +! do i = 1, ao_num +! write(*,'(100(F10.5,X))')ao_prod_center(2,i,:) +! enddo +! print*,'sigma = ' +! do i = 1, ao_num +! write(*,'(100(F10.5,X))')ao_prod_sigma(i,:) +! enddo + + + allocate(array(ao_num, ao_num, ao_num, ao_num)) + array = 0.d0 + allocate(array_ref(ao_num, ao_num, ao_num, ao_num)) + array_ref = 0.d0 + do ipoint = 1, n_points_final_grid + weight = final_weight_at_r_vector(ipoint) + do k = 1, ao_num + do l = 1, ao_num + do i = 1, ao_num + do j = 1, ao_num + array(j,i,l,k) += v_ij_u_cst_mu_j1b_test(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight + array_ref(j,i,l,k) += v_ij_u_cst_mu_j1b(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight + enddo + enddo + enddo + enddo + enddo + accu_relat = 0.d0 + accu_abs = 0.d0 + do k = 1, ao_num + do l = 1, ao_num + do i = 1, ao_num + do j = 1, ao_num + contrib = dabs(array(j,i,l,k) - array_ref(j,i,l,k)) + accu_abs += contrib + if(dabs(array_ref(j,i,l,k)).gt.1.d-10)then + accu_relat += contrib/dabs(array_ref(j,i,l,k)) + endif + enddo + enddo + enddo + enddo + print*,'accu_abs = ',accu_abs/dble(ao_num)**4 + print*,'accu_relat = ',accu_relat/dble(ao_num)**4 + + + +end + +subroutine routine_int2_grad1u2_grad2u2_j1b2 + implicit none + integer :: i,j,ipoint,k,l + double precision :: weight,accu_relat, accu_abs, contrib + double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:) + double precision, allocatable :: ints(:,:,:) + allocate(ints(ao_num, ao_num, n_points_final_grid)) + do ipoint = 1, n_points_final_grid + do i = 1, ao_num + do j = 1, ao_num + read(33,*)ints(j,i,ipoint) + enddo + enddo + enddo + + allocate(array(ao_num, ao_num, ao_num, ao_num)) + array = 0.d0 + allocate(array_ref(ao_num, ao_num, ao_num, ao_num)) + array_ref = 0.d0 + do ipoint = 1, n_points_final_grid + weight = final_weight_at_r_vector(ipoint) + do k = 1, ao_num + do l = 1, ao_num + do i = 1, ao_num + do j = 1, ao_num + array(j,i,l,k) += int2_grad1u2_grad2u2_j1b2_test_no_v(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight +! !array(j,i,l,k) += int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight + array_ref(j,i,l,k) += int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight +! !array(j,i,l,k) += ints(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight +! !array_ref(j,i,l,k) += int2_grad1u2_grad2u2_j1b2(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight +! !array_ref(j,i,l,k) += ints(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight + if(dabs(int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint)).gt.1.d-6)then + if(dabs(int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint) - int2_grad1u2_grad2u2_j1b2_test_no_v(j,i,ipoint)).gt.1.d-6)then + print*,j,i,ipoint + print*,int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint) , int2_grad1u2_grad2u2_j1b2_test_no_v(j,i,ipoint), dabs(int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint) - int2_grad1u2_grad2u2_j1b2_test_no_v(j,i,ipoint)) +! print*,int2_grad1u2_grad2u2_j1b2_test(i,j,ipoint) , int2_grad1u2_grad2u2_j1b2_test_no_v(i,j,ipoint), dabs(int2_grad1u2_grad2u2_j1b2_test(i,j,ipoint) - int2_grad1u2_grad2u2_j1b2_test_no_v(i,j,ipoint)) + stop + endif + endif + enddo + enddo + enddo + enddo + enddo + accu_relat = 0.d0 + accu_abs = 0.d0 + do k = 1, ao_num + do l = 1, ao_num + do i = 1, ao_num + do j = 1, ao_num + contrib = dabs(array(j,i,l,k) - array_ref(j,i,l,k)) + accu_abs += contrib + if(dabs(array_ref(j,i,l,k)).gt.1.d-10)then + accu_relat += contrib/dabs(array_ref(j,i,l,k)) + endif + enddo + enddo + enddo + enddo + print*,'accu_abs = ',accu_abs/dble(ao_num)**4 + print*,'accu_relat = ',accu_relat/dble(ao_num)**4 + + + +end diff --git a/src/tools/NEED b/src/tools/NEED index 46507542..f9a74966 100644 --- a/src/tools/NEED +++ b/src/tools/NEED @@ -3,3 +3,4 @@ mo_two_e_erf_ints aux_quantities hartree_fock dft_utils_in_r +two_body_rdm diff --git a/src/tools/molden.irp.f b/src/tools/molden.irp.f index 417b25ad..830a141e 100644 --- a/src/tools/molden.irp.f +++ b/src/tools/molden.irp.f @@ -52,8 +52,8 @@ program molden l += 1 if (l > ao_num) exit enddo - write(i_unit_output,*)'' enddo + write(i_unit_output,*)'' enddo diff --git a/src/tools/print_cfg.irp.f b/src/tools/print_cfg.irp.f new file mode 100644 index 00000000..9e9be7cf --- /dev/null +++ b/src/tools/print_cfg.irp.f @@ -0,0 +1,27 @@ +program print_energy + implicit none + BEGIN_DOC + ! Prints the energy of the wave function stored in the |EZFIO| directory. + END_DOC + + ! this has to be done in order to be sure that N_det, psi_det and + ! psi_coef_sorted are the wave function stored in the |EZFIO| directory. + read_wf = .True. + touch read_wf + PROVIDE N_states + call run +end + +subroutine run + implicit none + integer :: i,j + + do i=1,N_configuration + print *, i, sum(popcnt(psi_configuration(:,1,i))) + enddo + + print *, '' + do i=0,elec_num + print *, i, cfg_seniority_index(i) + enddo +end diff --git a/src/tools/print_dipole.irp.f b/src/tools/print_dipole.irp.f index 8351308e..92f2db36 100644 --- a/src/tools/print_dipole.irp.f +++ b/src/tools/print_dipole.irp.f @@ -1,5 +1,7 @@ program print_dipole implicit none + read_wf = .True. + SOFT_TOUCH read_wf call print_z_dipole_moment_only end diff --git a/src/tools/print_wf.irp.f b/src/tools/print_wf.irp.f index 3b9f6978..46cae2dc 100644 --- a/src/tools/print_wf.irp.f +++ b/src/tools/print_wf.irp.f @@ -33,8 +33,9 @@ subroutine routine double precision :: norm_mono_a,norm_mono_b double precision :: norm_mono_a_2,norm_mono_b_2 double precision :: norm_mono_a_pert_2,norm_mono_b_pert_2 - double precision :: norm_mono_a_pert,norm_mono_b_pert + double precision :: norm_mono_a_pert,norm_mono_b_pert,norm_double_1 double precision :: delta_e,coef_2_2 + norm_mono_a = 0.d0 norm_mono_b = 0.d0 norm_mono_a_2 = 0.d0 @@ -43,6 +44,7 @@ subroutine routine norm_mono_b_pert = 0.d0 norm_mono_a_pert_2 = 0.d0 norm_mono_b_pert_2 = 0.d0 + norm_double_1 = 0.d0 do i = 1, min(N_det_print_wf,N_det) print*,'' print*,'i = ',i @@ -94,6 +96,7 @@ subroutine routine print*,'h1,p1 = ',h1,p1 print*,'s2',s2 print*,'h2,p2 = ',h2,p2 + norm_double_1 += dabs(psi_coef_sorted(i,1)/psi_coef_sorted(1,1)) endif print*,' = ',hij @@ -110,6 +113,7 @@ subroutine routine print*,'' print*,'L1 norm of mono alpha = ',norm_mono_a print*,'L1 norm of mono beta = ',norm_mono_b + print*,'L1 norm of double exc = ',norm_double_1 print*, '---' print*,'L2 norm of mono alpha = ',norm_mono_a_2 print*,'L2 norm of mono beta = ',norm_mono_b_2 diff --git a/src/tools/truncate_wf.irp.f b/src/tools/truncate_wf.irp.f new file mode 100644 index 00000000..8ed97391 --- /dev/null +++ b/src/tools/truncate_wf.irp.f @@ -0,0 +1,111 @@ +program truncate_wf + implicit none + BEGIN_DOC +! Truncate the wave function + END_DOC + read_wf = .True. + if (s2_eig) then + call routine_s2 + else + call routine + endif +end + +subroutine routine + implicit none + integer :: ndet_max + print*, 'Max number of determinants ?' + read(5,*) ndet_max + integer(bit_kind), allocatable :: psi_det_tmp(:,:,:) + double precision, allocatable :: psi_coef_tmp(:,:) + allocate(psi_det_tmp(N_int,2,ndet_max),psi_coef_tmp(ndet_max, N_states)) + + integer :: i,j + double precision :: accu(N_states) + accu = 0.d0 + do i = 1, ndet_max + do j = 1, N_int + psi_det_tmp(j,1,i) = psi_det_sorted(j,1,i) + psi_det_tmp(j,2,i) = psi_det_sorted(j,2,i) + enddo + do j = 1, N_states + psi_coef_tmp(i,j) = psi_coef_sorted(i,j) + accu(j) += psi_coef_tmp(i,j) **2 + enddo + enddo + do j = 1, N_states + accu(j) = 1.d0/dsqrt(accu(j)) + enddo + do j = 1, N_states + do i = 1, ndet_max + psi_coef_tmp(i,j) = psi_coef_tmp(i,j) * accu(j) + enddo + enddo + + call save_wavefunction_general(ndet_max,N_states,psi_det_tmp,size(psi_coef_tmp,1),psi_coef_tmp) + +end + +subroutine routine_s2 + implicit none + integer :: ndet_max + double precision :: wmin + integer(bit_kind), allocatable :: psi_det_tmp(:,:,:) + double precision, allocatable :: psi_coef_tmp(:,:) + integer :: i,j,k + double precision :: accu(N_states) + + integer :: weights(0:16), ix + double precision :: x + + weights(:) = 0 + do i=1,N_det + x = -dlog(1.d-32+sum(weight_configuration(det_to_configuration(i),:)))/dlog(10.d0) + ix = min(int(x), 16) + weights(ix) += 1 + enddo + + print *, 'Histogram of the weights of the CFG' + do i=0,15 + print *, ' 10^{-', i, '} ', weights(i) + end do + print *, '< 10^{-', 15, '} ', weights(16) + + + print*, 'Min weight of the configuration?' + read(5,*) wmin + + ndet_max = 0 + do i=1,N_det + if (maxval(weight_configuration( det_to_configuration(i),:)) < wmin) cycle + ndet_max = ndet_max+1 + enddo + + allocate(psi_det_tmp(N_int,2,ndet_max),psi_coef_tmp(ndet_max, N_states)) + + accu = 0.d0 + k=0 + do i = 1, N_det + if (maxval(weight_configuration( det_to_configuration(i),:)) < wmin) cycle + k = k+1 + do j = 1, N_int + psi_det_tmp(j,1,k) = psi_det(j,1,i) + psi_det_tmp(j,2,k) = psi_det(j,2,i) + enddo + do j = 1, N_states + psi_coef_tmp(k,j) = psi_coef(i,j) + accu(j) += psi_coef_tmp(k,j)**2 + enddo + enddo + do j = 1, N_states + accu(j) = 1.d0/dsqrt(accu(j)) + enddo + do j = 1, N_states + do i = 1, ndet_max + psi_coef_tmp(i,j) = psi_coef_tmp(i,j) * accu(j) + enddo + enddo + + call save_wavefunction_general(ndet_max,N_states,psi_det_tmp,size(psi_coef_tmp,1),psi_coef_tmp) + +end diff --git a/src/two_rdm_routines/davidson_like_state_av_2rdm.irp.f b/src/two_rdm_routines/davidson_like_state_av_2rdm.irp.f index eb247dea..26ed5ae6 100644 --- a/src/two_rdm_routines/davidson_like_state_av_2rdm.irp.f +++ b/src/two_rdm_routines/davidson_like_state_av_2rdm.irp.f @@ -529,10 +529,14 @@ subroutine orb_range_2_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,list_ c_average += c_1(l) * c_1(l) * state_weights(l) enddo - call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm) + if (nkeys > 0) then + call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm) + endif nkeys = 0 call orb_range_diag_to_all_2_rdm_dm_buffer(tmp_det,c_average,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) - call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm) + if (nkeys > 0) then + call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm) + endif nkeys = 0 end do diff --git a/src/utils/constants.include.F b/src/utils/constants.include.F index a96fabe6..d1727701 100644 --- a/src/utils/constants.include.F +++ b/src/utils/constants.include.F @@ -14,7 +14,7 @@ double precision, parameter :: thresh = 1.d-15 double precision, parameter :: cx_lda = -0.73855876638202234d0 double precision, parameter :: c_2_4_3 = 2.5198420997897464d0 double precision, parameter :: cst_lda = -0.93052573634909996d0 -double precision, parameter :: c_4_3 = 1.3333333333333333d0 -double precision, parameter :: c_1_3 = 0.3333333333333333d0 +double precision, parameter :: c_4_3 = 4.d0/3.d0 +double precision, parameter :: c_1_3 = 1.d0/3.d0 double precision, parameter :: sq_op5 = dsqrt(0.5d0) double precision, parameter :: dlog_2pi = dlog(2.d0*dacos(-1.d0)) diff --git a/src/utils/format_w_error.irp.f b/src/utils/format_w_error.irp.f new file mode 100644 index 00000000..1378d367 --- /dev/null +++ b/src/utils/format_w_error.irp.f @@ -0,0 +1,71 @@ +subroutine format_w_error(value,error,size_nb,max_nb_digits,format_value,str_error) + + implicit none + + BEGIN_DOC + ! Format for double precision, value(error) + END_DOC + + ! in + ! | value | double precision | value... | + ! | error | double precision | error... | + ! | size_nb | integer | X in FX.Y | + ! | max_nb_digits | integer | Max Y in FX.Y | + + ! out + ! | format_value | character | string FX.Y for the format | + ! | str_error | character | string of the error | + + ! internal + ! | str_size | character | size in string format | + ! | nb_digits | integer | number of digits Y in FX.Y depending of the error | + ! | str_nb_digits | character | nb_digits in string format | + ! | str_exp | character | string of the value in exponential format | + + ! in + double precision, intent(in) :: error, value + integer, intent(in) :: size_nb, max_nb_digits + + ! out + character(len=20), intent(out) :: str_error, format_value + + ! internal + character(len=20) :: str_size, str_nb_digits, str_exp + integer :: nb_digits + + ! max_nb_digit: Y max + ! size_nb = Size of the double: X (FX.Y) + write(str_size,'(I3)') size_nb + + ! Error + write(str_exp,'(1pE20.0)') error + str_error = trim(adjustl(str_exp)) + + ! Number of digit: Y (FX.Y) from the exponent + str_nb_digits = str_exp(19:20) + read(str_nb_digits,*) nb_digits + + ! If the error is 0d0 + if (error <= 1d-16) then + write(str_nb_digits,*) max_nb_digits + endif + + ! If the error is too small + if (nb_digits > max_nb_digits) then + write(str_nb_digits,*) max_nb_digits + str_error(1:1) = '0' + endif + + ! If the error is too big (>= 0.5) + if (error >= 0.5d0) then + str_nb_digits = '1' + str_error(1:1) = '*' + endif + + ! FX.Y,A1,A1,A1 for value(str_error) + !string = 'F'//trim(adjustl(str_size))//'.'//trim(adjustl(str_nb_digits))//',A1,A1,A1' + + ! FX.Y just for the value + format_value = 'F'//trim(adjustl(str_size))//'.'//trim(adjustl(str_nb_digits)) + +end diff --git a/src/utils/linear_algebra.irp.f b/src/utils/linear_algebra.irp.f index 2f84e753..38a8cad2 100644 --- a/src/utils/linear_algebra.irp.f +++ b/src/utils/linear_algebra.irp.f @@ -1981,3 +1981,54 @@ subroutine diag_nonsym_right(n, A, A_ldim, V, V_ldim, energy, E_ldim) end subroutine diag_nonsym_right ! --- + + +subroutine pivoted_cholesky( A, rank, tol, ndim, U) +! +! A = U**T * U +! +! matrix A is destroyed inside this subroutine +! Cholesky vectors are stored in U +! dimension of U: U(1:rank, 1:n) +! U is allocated inside this subroutine +! rank is the number of Cholesky vectors depending on tol +! +integer :: ndim +integer, intent(inout) :: rank +double precision, dimension(ndim, ndim), intent(inout) :: A +double precision, dimension(ndim, rank), intent(out) :: U +double precision, intent(in) :: tol + +integer, dimension(:), allocatable :: piv +double precision, dimension(:), allocatable :: work +character, parameter :: uplo = "U" +integer :: N, LDA +integer :: info +integer :: k, l, rank0 +external :: dpstrf + +rank0 = rank +N = size(A, dim=1) +LDA = N +allocate(piv(N)) +allocate(work(2*N)) +call dpstrf(uplo, N, A, LDA, piv, rank, tol, work, info) + +if (rank > rank0) then + print *, 'Bug: rank > rank0 in pivoted cholesky. Increase rank before calling' + stop +end if + +do k = 1, N + A(k+1:, k) = 0.00D+0 +end do +! TODO: It should be possible to use only one vector of size (1:rank) as a buffer +! to do the swapping in-place +U = 0.00D+0 +do k = 1, N + l = piv(k) + U(l, :) = A(1:rank, k) +end do + +end subroutine pivoted_cholesky + diff --git a/src/utils/memory.irp.f b/src/utils/memory.irp.f index 3ea242b0..d5a066a1 100644 --- a/src/utils/memory.irp.f +++ b/src/utils/memory.irp.f @@ -114,7 +114,7 @@ subroutine print_memory_usage() call resident_memory(rss) call total_memory(mem) - write(*,'(A,F14.6,A,F14.6,A)') & + write(*,'(A,F14.3,A,F14.3,A)') & '.. >>>>> [ RES MEM : ', rss , & ' GB ] [ VIRT MEM : ', mem, ' GB ] <<<<< ..' end diff --git a/src/utils/set_multiple_levels_omp.irp.f b/src/utils/set_multiple_levels_omp.irp.f new file mode 100644 index 00000000..572a13f4 --- /dev/null +++ b/src/utils/set_multiple_levels_omp.irp.f @@ -0,0 +1,26 @@ +subroutine set_multiple_levels_omp(activate) + + BEGIN_DOC +! If true, activate OpenMP nested parallelism. If false, deactivate. + END_DOC + + implicit none + logical, intent(in) :: activate + + if (activate) then + call omp_set_max_active_levels(3) + + IRP_IF SET_NESTED + call omp_set_nested(.True.) + IRP_ENDIF + + else + + call omp_set_max_active_levels(1) + + IRP_IF SET_NESTED + call omp_set_nested(.False.) + IRP_ENDIF + end if + +end diff --git a/src/utils/units.irp.f b/src/utils/units.irp.f new file mode 100644 index 00000000..1850b28b --- /dev/null +++ b/src/utils/units.irp.f @@ -0,0 +1,22 @@ +BEGIN_PROVIDER [double precision, ha_to_ev] + + implicit none + BEGIN_DOC + ! Converstion from Hartree to eV + END_DOC + + ha_to_ev = 27.211396641308d0 + +END_PROVIDER + +BEGIN_PROVIDER [double precision, au_to_D] + + implicit none + BEGIN_DOC + ! Converstion from au to Debye + END_DOC + + au_to_D = 2.5415802529d0 + +END_PROVIDER + diff --git a/src/utils/util.irp.f b/src/utils/util.irp.f index c9e35ee1..e7f00ce2 100644 --- a/src/utils/util.irp.f +++ b/src/utils/util.irp.f @@ -37,6 +37,10 @@ double precision function binom_func(i,j) else binom_func = dexp( logfact(i)-logfact(j)-logfact(i-j) ) endif + + ! To avoid .999999 numbers + binom_func = floor(binom_func + 0.5d0) + end @@ -328,7 +332,7 @@ BEGIN_PROVIDER [ integer, nproc ] ! Number of current OpenMP threads END_DOC - integer :: omp_get_num_threads + integer, external :: omp_get_num_threads nproc = 1 !$OMP PARALLEL !$OMP MASTER diff --git a/src/zmq/f77_zmq_free.h b/src/zmq/f77_zmq_free.h deleted file mode 120000 index ac5e33cd..00000000 --- a/src/zmq/f77_zmq_free.h +++ /dev/null @@ -1 +0,0 @@ -../../include/f77_zmq_free.h \ No newline at end of file diff --git a/src/zmq/f77_zmq_module.f90 b/src/zmq/f77_zmq_module.F90 similarity index 50% rename from src/zmq/f77_zmq_module.f90 rename to src/zmq/f77_zmq_module.F90 index 1e4a5af3..9171f11c 100644 --- a/src/zmq/f77_zmq_module.f90 +++ b/src/zmq/f77_zmq_module.F90 @@ -1,4 +1,4 @@ module f77_zmq - include 'f77_zmq_free.h' +#include "f77_zmq_free.h" end module