LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-11-12 16:33:37 +01:00
Code Releases Activity

Merge pull request #4 from eginer/dev

Browse Source 
Dev
This commit is contained in:
AbdAmmar 2022-12-10 14:20:46 +01:00 committed by GitHub
commit 441ed8ee6b
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
102 changed files with 2479 additions and 1993 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
.gitignore vendored
View File

@ -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

170
RELEASE_NOTES.org
View File

@ -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

21
configure vendored
View File

@ -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

2
etc/qp.rc
View File

@ -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
;;

2
external/.gitignore vendored
View File

@ -1,2 +1,2 @@
#*
*

2
external/qp2-dependencies vendored

@ -1 +1 @@
Subproject commit 242151e03d1d6bf042387226431d82d35845686a
Subproject commit 90ee61f5041c7c94a0c605625a264860292813a0

9
include/.gitignore vendored
View File

@ -1,8 +1 @@
zmq.h
gmp.h
zconf.h
zconf.h
zlib.h
zmq_utils.h
f77_zmq_free.h
f77_zmq.h
*

19
ocaml/Command_line.ml
View File

@ -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.")
)
;;

2
ocaml/Command_line.mli
View File

@ -59,6 +59,8 @@ let () =
*)
exception Error of string
type short_opt = char
type long_opt = string

113
ocaml/Input_ao_two_e_eff_pot.ml
View File

@ -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

87
ocaml/Input_bi_ortho_mos.ml
View File

@ -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

113
ocaml/Input_cassd.ml
View File

@ -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

243
ocaml/Input_cipsi_deb.ml
View File

@ -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

351
ocaml/Input_tc_h_clean.ml
View File

@ -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

143
ocaml/Input_tc_scf.ml
View File

@ -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

2
ocaml/Makefile
View File

@ -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 ; \

4
ocaml/Molecule.ml
View File

@ -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 ;

2
ocaml/Progress_bar.ml
View File

@ -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 () }

4
ocaml/Qputils.ml
View File

@ -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 )

2
ocaml/TaskServer.ml
View File

@ -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 =

13
ocaml/qp_create_ezfio.ml
View File

@ -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

139
ocaml/qp_tunnel.ml
View File

@ -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 <host> &
%!"
%!"
(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"

3
scripts/.gitignore vendored
View File

@ -2,3 +2,6 @@
*.pyo
docopt.py
resultsFile/
verif_omp/a.out
src/*/Makefile
src/*/*/

29
scripts/compilation/qp_create_ninja
View File

@ -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("")

11
src/.gitignore vendored Normal file
View File

@ -0,0 +1,11 @@
*
!README.rst
!*/
*/*
!*/*.*
*/*.o
*/build.ninja
*/ezfio_interface.irp.f
*/.gitignore
*/*.swp

1
src/ao_basis/EZFIO.cfg
View File

@ -67,4 +67,3 @@ doc: Use normalized primitive functions
interface: ezfio, provider
default: true

2
src/ao_basis/spherical_to_cartesian.irp.f
View File

@ -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) ]

284
src/ao_many_one_e_ints/grad2_jmu_manu.irp.f Normal file
View File

@ -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

7
src/ao_many_one_e_ints/grad2_jmu_modif.irp.f
View File

@ -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

287
src/ao_many_one_e_ints/grad_lapl_jmu_manu.irp.f Normal file
View File

@ -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
! ---

1
src/ao_many_one_e_ints/listj1b.irp.f
View File

@ -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)

191
src/ao_many_one_e_ints/listj1b_sorted.irp.f Normal file
View File

@ -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

9
src/ao_one_e_ints/pot_ao_ints.irp.f
View File

@ -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

4
src/ao_one_e_ints/pot_ao_pseudo_ints.irp.f
View File

@ -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

4
src/ao_one_e_ints/pseudopot.f90
View File

@ -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

18
src/ao_two_e_ints/map_integrals.irp.f
View File

@ -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 : <ij|kl>
! physicist convention : <ij|kl>
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 : <ij|kl>
! physicist convention : <ij|kl>
END_DOC
implicit none
integer, intent(in) :: j,k,l, sze

191
src/ao_two_e_ints/test_cosgtos_1e.irp.f
View File

@ -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
! ---

165
src/ao_two_e_ints/test_cosgtos_2e.irp.f
View File

@ -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
! ---

12
src/ao_two_e_ints/two_e_integrals.irp.f
View File

@ -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

15
src/bitmask/bitmasks_routines.irp.f
View File

@ -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

12
src/cipsi/EZFIO.cfg
View File

@ -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

1
src/cipsi/NEED
View File

@ -2,5 +2,4 @@ perturbation
zmq
mpi
iterations
two_body_rdm
csf

116
src/cipsi/pt2_stoch_routines.irp.f
View File

@ -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)

96
src/cipsi/selection.irp.f
View File

@ -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)

92
src/cipsi/selection_buffer.irp.f
View File

@ -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)

24
src/cipsi/selection_weight.irp.f
View File

@ -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

6
src/cipsi/slave_cipsi.irp.f
View File

@ -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'
!---

2
src/cipsi/zmq_selection.irp.f
View File

@ -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')

9
src/cis/cis.irp.f
View File

@ -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)

2
src/cisd/30.cisd.bats
View File

@ -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

2
src/cisd/cisd.irp.f
View File

@ -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

24
src/csf/configurations.irp.f
View File

@ -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

1
src/csf/conversion.irp.f
View File

@ -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

10
src/dav_general_mat/dav_general.irp.f
View File

@ -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><i|H|u_k>
! -----------------------------------
@ -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 = <u_k | W_l> = <u_k| H |u_l>
! -------------------------------------------

17
src/davidson/davidson_parallel.irp.f
View File

@ -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 <Psi|H|Psi>')
! endif
!END_PROVIDER
integer function zmq_put_N_states_diag(zmq_to_qp_run_socket,worker_id)
use f77_zmq
implicit none

3
src/davidson/davidson_parallel_csf.irp.f
View File

@ -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

3
src/davidson/davidson_parallel_nos2.irp.f
View File

@ -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

10
src/davidson/diagonalization_hcsf_dressed.irp.f
View File

@ -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><i|H|u_k>
! -----------------------------------
@ -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

32
src/davidson/diagonalization_hs2_dressed.irp.f
View File

@ -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><i|H|u_k>
! -----------------------------------
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

10
src/davidson/diagonalization_nonsym_h_dressed.irp.f
View File

@ -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

1
src/davidson/diagonalize_ci.irp.f
View File

@ -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)

39
src/davidson/input.irp.f
View File

@ -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
!

1
src/davidson_keywords/README.rst
View File

@ -2,3 +2,4 @@
davidson_keywords
=================
Keywords used for Davidson algorithms.

12
src/determinants/EZFIO.cfg
View File

@ -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

12
src/determinants/density_matrix.irp.f
View File

@ -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

6
src/determinants/s2.irp.f
View File

@ -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

8
src/determinants/slater_rules.irp.f
View File

@ -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

3
src/determinants/slater_rules_wee_mono.irp.f
View File

@ -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)

2
src/determinants/spindeterminants.irp.f
View File

@ -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)

116
src/dft_utils_in_r/ao_prod_mlti_pl.irp.f Normal file
View File

@ -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

2
src/dressing/alpha_factory.irp.f
View File

@ -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)

6
src/dressing/run_dress_slave.irp.f
View File

@ -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

2
src/ezfio_files/output.irp.f
View File

@ -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

3
src/iterations/print_extrapolation.irp.f
View File

@ -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

14
src/iterations/print_summary.irp.f
View File

@ -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

3
src/mo_basis/mos_in_r.irp.f
View File

@ -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)

12
src/mo_two_e_erf_ints/map_integrals_erf.irp.f
View File

@ -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)

110
src/mo_two_e_ints/mo_bi_integrals.irp.f
View File

@ -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

17
src/non_h_ints_mu/j12_nucl_utils.irp.f
View File

@ -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

16
src/non_hermit_dav/biorthog.irp.f
View File

@ -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 )

15
src/scf_utils/diagonalize_fock.irp.f
View File

@ -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))

192
src/scf_utils/fock_matrix.irp.f
View File

@ -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)

36
src/scf_utils/scf_density_matrix_ao.irp.f
View File

@ -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

2
src/tc_bi_ortho/test_tc_fock.irp.f
View File

@ -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

2
src/tc_keywords/EZFIO.cfg
View File

@ -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

346
src/tc_scf/test_int.irp.f Normal file
View File

@ -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

1
src/tools/NEED
View File

@ -3,3 +3,4 @@ mo_two_e_erf_ints
aux_quantities
hartree_fock
dft_utils_in_r
two_body_rdm

2
src/tools/molden.irp.f
View File

@ -52,8 +52,8 @@ program molden
l += 1
if (l > ao_num) exit
enddo
write(i_unit_output,*)''
enddo
write(i_unit_output,*)''
enddo

27
src/tools/print_cfg.irp.f Normal file
View File

@ -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

2
src/tools/print_dipole.irp.f
View File

@ -1,5 +1,7 @@
program print_dipole
implicit none
read_wf = .True.
SOFT_TOUCH read_wf
call print_z_dipole_moment_only
end

6
src/tools/print_wf.irp.f
View File

@ -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*,'<Ref| H |D_I> = ',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

111
src/tools/truncate_wf.irp.f Normal file
View File

@ -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

8
src/two_rdm_routines/davidson_like_state_av_2rdm.irp.f
View File

@ -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

4
src/utils/constants.include.F
View File

@ -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))

71
src/utils/format_w_error.irp.f Normal file
View File

@ -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

51
src/utils/linear_algebra.irp.f
View File

@ -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

2
src/utils/memory.irp.f
View File

@ -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

26
src/utils/set_multiple_levels_omp.irp.f Normal file
View File

@ -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

22
src/utils/units.irp.f Normal file
View File

@ -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

6
src/utils/util.irp.f
View File

@ -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

Some files were not shown because too many files have changed in this diff Show More