From 7fa1527a3de79015fc70aac4053ed84cc1ab50fc Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sun, 16 May 2021 01:32:55 +0200 Subject: [PATCH 01/29] Move restore_symm keyword --- RELEASE_NOTES.org | 2 ++ bin/qp_export_as_tgz | 6 ++++-- src/mo_one_e_ints/EZFIO.cfg | 7 ------- 3 files changed, 6 insertions(+), 9 deletions(-) diff --git a/RELEASE_NOTES.org b/RELEASE_NOTES.org index f9dbcbb0..02176ffa 100644 --- a/RELEASE_NOTES.org +++ b/RELEASE_NOTES.org @@ -52,6 +52,8 @@ - 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 *** Code diff --git a/bin/qp_export_as_tgz b/bin/qp_export_as_tgz index 7332d1ed..20624dba 100755 --- a/bin/qp_export_as_tgz +++ b/bin/qp_export_as_tgz @@ -99,7 +99,9 @@ function find_libs () { } function find_exec () { - find ${QP_ROOT}/$1 -perm /u+x -type f + for i in $@ ; do + find ${QP_ROOT}/$i -perm /u+x -type f + done } @@ -119,7 +121,7 @@ fi echo "Copying binary files" # -------------------- -FORTRAN_EXEC=$(find_exec src) +FORTRAN_EXEC=$(find_exec src/*/) if [[ -z $FORTRAN_EXEC ]] ; then error 'No Fortran binaries found.' exit 1 diff --git a/src/mo_one_e_ints/EZFIO.cfg b/src/mo_one_e_ints/EZFIO.cfg index ca4dabf4..d58b3da1 100644 --- a/src/mo_one_e_ints/EZFIO.cfg +++ b/src/mo_one_e_ints/EZFIO.cfg @@ -47,10 +47,3 @@ type: Disk_access doc: Read/Write |MO| one-electron integrals from/to disk [ Write | Read | None ] interface: ezfio,provider,ocaml default: None - - -[restore_symm] -type: logical -doc: If true, try to find symmetry in the MO coefficient matrices -interface: ezfio,provider,ocaml -default: True From c226ffd9dfdbd05f822884f932956cd356191808 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sun, 16 May 2021 11:24:05 +0200 Subject: [PATCH 02/29] forgot file --- src/utils/EZFIO.cfg | 5 +++++ 1 file changed, 5 insertions(+) create mode 100644 src/utils/EZFIO.cfg diff --git a/src/utils/EZFIO.cfg b/src/utils/EZFIO.cfg new file mode 100644 index 00000000..7d367f0c --- /dev/null +++ b/src/utils/EZFIO.cfg @@ -0,0 +1,5 @@ +[restore_symm] +type: logical +doc: If true, try to find symmetry in the MO coefficient matrices +interface: ezfio,provider,ocaml +default: False From f543d386a9aa6ecc97abe06be45b3cf295a6a903 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 17 May 2021 15:08:44 +0200 Subject: [PATCH 03/29] Fix minor bug in multi-state PT2 --- src/cipsi/selection.irp.f | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f index e688a7bb..62c5dc27 100644 --- a/src/cipsi/selection.irp.f +++ b/src/cipsi/selection.irp.f @@ -735,7 +735,7 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d double precision :: eigvalues(N_states+1) double precision :: work(1+6*(N_states+1)+2*(N_states+1)**2) - integer :: iwork(3+5*(N_states+1)), info, k + integer :: info, k , iwork(N_states+1) if (do_diag) then double precision :: pt2_matrix(N_states+1,N_states+1) @@ -747,8 +747,8 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d pt2_matrix(N_states+1,istate) = mat(istate,p1,p2) enddo - call DSYEVD( 'V', 'U', N_states+1, pt2_matrix, N_states+1, eigvalues, & - work, size(work), iwork, size(iwork), info ) + call DSYEV( 'V', 'U', N_states+1, pt2_matrix, N_states+1, eigvalues, & + work, size(work), info ) if (info /= 0) then print *, 'error in '//irp_here stop -1 @@ -756,7 +756,7 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d pt2_matrix = dabs(pt2_matrix) iwork(1:N_states+1) = maxloc(pt2_matrix,DIM=1) do k=1,N_states - e_pert(iwork(k)) = eigvalues(k) - E0(iwork(k)) + e_pert(k) = eigvalues(iwork(k)) - E0(k) enddo endif From bb0743ac1ee439f387809e83f6927883aa320275 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 19 May 2021 18:38:11 +0200 Subject: [PATCH 04/29] Updated IRPF90 --- etc/local.rc | 2 +- external/irpf90 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/etc/local.rc b/etc/local.rc index 954b315b..0e212c24 100644 --- a/etc/local.rc +++ b/etc/local.rc @@ -19,4 +19,4 @@ # export QP_NIC=lo # export QP_NIC=ib0 - +export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/home/scemama/TREX/trexio/_install/lib diff --git a/external/irpf90 b/external/irpf90 index 132a4a16..33ca5e10 160000 --- a/external/irpf90 +++ b/external/irpf90 @@ -1 +1 @@ -Subproject commit 132a4a1661c9878d21dcbf0ac14f7fe9a3b110d0 +Subproject commit 33ca5e1018f3bbb5e695e6ee558f5dac0753b271 From 1dc2350b328decf2b2b35f4fe5cdcf61e0d85e9d Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 19 May 2021 18:43:48 +0200 Subject: [PATCH 05/29] revert etc/local.rc --- etc/local.rc | 1 - 1 file changed, 1 deletion(-) diff --git a/etc/local.rc b/etc/local.rc index 0e212c24..182599c6 100644 --- a/etc/local.rc +++ b/etc/local.rc @@ -19,4 +19,3 @@ # export QP_NIC=lo # export QP_NIC=ib0 -export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/home/scemama/TREX/trexio/_install/lib From 3ceea5e8b27b1ce82286ac8254964fd8e784e6df Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 20 May 2021 18:19:44 +0200 Subject: [PATCH 06/29] Comment in 2rdm --- src/two_body_rdm/two_e_dm_mo.irp.f | 14 +++++--------- 1 file changed, 5 insertions(+), 9 deletions(-) diff --git a/src/two_body_rdm/two_e_dm_mo.irp.f b/src/two_body_rdm/two_e_dm_mo.irp.f index 19e8d75e..4dadd2e6 100644 --- a/src/two_body_rdm/two_e_dm_mo.irp.f +++ b/src/two_body_rdm/two_e_dm_mo.irp.f @@ -5,18 +5,14 @@ BEGIN_PROVIDER [double precision, two_e_dm_mo, (mo_num,mo_num,mo_num,mo_num)] ! ! ! - ! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO ALL OCCUPIED ORBITALS : core, inactive and active + ! where the indices (i,j,k,l) belong to all MOs. ! - ! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{elec} * (N_{elec} - 1)/2 + ! The normalization (i.e. sum of diagonal elements) is set to $N_{elec} * (N_{elec} - 1)/2$ ! - ! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act" + ! !!!!! WARNING !!!!! IF "no_core_density" then all elements involving at least one CORE MO are set to zero + ! The state-averaged two-electron energy : ! - ! !!!!! WARNING !!!!! IF "no_core_density" then all elements involving at least one CORE MO is set to zero - ! The two-electron energy of each state can be computed as: - ! - ! \sum_{i,j,k,l = 1, n_core_inact_act_orb} two_e_dm_mo(i,j,k,l,istate) * < ii jj | kk ll > - ! - ! with ii = list_core_inact_act(i), jj = list_core_inact_act(j), kk = list_core_inact_act(k), ll = list_core_inact_act(l) + ! \sum_{i,j,k,l = 1, mo_num} two_e_dm_mo(i,j,k,l) * < ii jj | kk ll > END_DOC two_e_dm_mo = 0.d0 integer :: i,j,k,l,iorb,jorb,korb,lorb,istate From 68b33e4b356a26723123e65ec90763627babb946 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 21 May 2021 16:42:48 +0200 Subject: [PATCH 07/29] Changed Symmetry into Angmom in OCaml --- ocaml/{Symmetry.ml => Angmom.ml} | 6 ++-- ocaml/{Symmetry.mli => Angmom.mli} | 0 ocaml/GaussianPrimitive.ml | 4 +-- ocaml/Gto.ml | 8 ++--- ocaml/Gto.mli | 2 +- ocaml/Input_ao_basis.ml | 18 ++++++------ ocaml/Long_basis.ml | 8 ++--- ocaml/Long_basis.mli | 8 ++--- ocaml/Primitive.mli | 4 +-- ocaml/qp_create_ezfio.ml | 7 +++-- src/basis/EZFIO.cfg | 47 ++++++++++++++++++++++++++++++ src/basis/NEED | 1 + src/basis/README.rst | 8 +++++ src/basis/basis.irp.f | 32 ++++++++++++++++++++ 14 files changed, 121 insertions(+), 32 deletions(-) rename ocaml/{Symmetry.ml => Angmom.ml} (95%) rename ocaml/{Symmetry.mli => Angmom.mli} (100%) create mode 100644 src/basis/EZFIO.cfg create mode 100644 src/basis/NEED create mode 100644 src/basis/README.rst create mode 100644 src/basis/basis.irp.f diff --git a/ocaml/Symmetry.ml b/ocaml/Angmom.ml similarity index 95% rename from ocaml/Symmetry.ml rename to ocaml/Angmom.ml index eb4b637b..ed13e8dc 100644 --- a/ocaml/Symmetry.ml +++ b/ocaml/Angmom.ml @@ -26,7 +26,7 @@ let of_string = function | "J" | "j" -> J | "K" | "k" -> K | "L" | "l" -> L - | x -> raise (Failure ("Symmetry should be S|P|D|F|G|H|I|J|K|L, + | x -> raise (Failure ("Angmom should be S|P|D|F|G|H|I|J|K|L, not "^x^".")) let of_char = function @@ -40,7 +40,7 @@ let of_char = function | 'J' | 'j' -> J | 'K' | 'k' -> K | 'L' | 'l' -> L - | x -> raise (Failure ("Symmetry should be S|P|D|F|G|H|I|J|K|L")) + | x -> raise (Failure ("Angmom should be S|P|D|F|G|H|I|J|K|L")) let to_l = function | S -> Positive_int.of_int 0 @@ -68,7 +68,7 @@ let of_l i = | 7 -> J | 8 -> K | 9 -> L - | x -> raise (Failure ("Symmetry should be S|P|D|F|G|H|I|J|K|L")) + | x -> raise (Failure ("Angmom should be S|P|D|F|G|H|I|J|K|L")) type st = t diff --git a/ocaml/Symmetry.mli b/ocaml/Angmom.mli similarity index 100% rename from ocaml/Symmetry.mli rename to ocaml/Angmom.mli diff --git a/ocaml/GaussianPrimitive.ml b/ocaml/GaussianPrimitive.ml index d2144dec..c59a68de 100644 --- a/ocaml/GaussianPrimitive.ml +++ b/ocaml/GaussianPrimitive.ml @@ -2,14 +2,14 @@ open Qptypes open Sexplib.Std type t = - { sym : Symmetry.t ; + { sym : Angmom.t ; expo : AO_expo.t ; } [@@deriving sexp] let to_string p = let { sym = s ; expo = e } = p in Printf.sprintf "(%s, %22e)" - (Symmetry.to_string s) + (Angmom.to_string s) (AO_expo.to_float e) diff --git a/ocaml/Gto.ml b/ocaml/Gto.ml index 465e5627..200b4fa6 100644 --- a/ocaml/Gto.ml +++ b/ocaml/Gto.ml @@ -9,7 +9,7 @@ type fmt = | Gaussian type t = -{ sym : Symmetry.t ; +{ sym : Angmom.t ; lc : ((GaussianPrimitive.t * AO_coef.t) list) } [@@deriving sexp] @@ -47,7 +47,7 @@ let read_one in_channel = in let sym_str = String.sub buffer 0 2 in let n_str = String.sub buffer 2 ((String.length buffer)-2) in - let sym = Symmetry.of_string (String_ext.strip sym_str) in + let sym = Angmom.of_string (String_ext.strip sym_str) in let n = int_of_string (String_ext.strip n_str) in (* Read all the primitives *) let rec read_lines result = function @@ -82,7 +82,7 @@ let read_one in_channel = (** Write the GTO in Gamess format *) let to_string_gamess { sym = sym ; lc = lc } = let result = - Printf.sprintf "%s %3d" (Symmetry.to_string sym) (List.length lc) + Printf.sprintf "%s %3d" (Angmom.to_string sym) (List.length lc) in let rec do_work accu i = function | [] -> List.rev accu @@ -102,7 +102,7 @@ let to_string_gamess { sym = sym ; lc = lc } = (** Write the GTO in Gaussian format *) let to_string_gaussian { sym = sym ; lc = lc } = let result = - Printf.sprintf "%s %3d 1.00" (Symmetry.to_string sym) (List.length lc) + Printf.sprintf "%s %3d 1.00" (Angmom.to_string sym) (List.length lc) in let rec do_work accu i = function | [] -> List.rev accu diff --git a/ocaml/Gto.mli b/ocaml/Gto.mli index 91534ebe..ccbac7fe 100644 --- a/ocaml/Gto.mli +++ b/ocaml/Gto.mli @@ -5,7 +5,7 @@ type fmt = | Gaussian type t = - { sym : Symmetry.t ; + { sym : Angmom.t ; lc : (GaussianPrimitive.t * Qptypes.AO_coef.t) list; } [@@deriving sexp] diff --git a/ocaml/Input_ao_basis.ml b/ocaml/Input_ao_basis.ml index a73d641c..95d37a7a 100644 --- a/ocaml/Input_ao_basis.ml +++ b/ocaml/Input_ao_basis.ml @@ -9,7 +9,7 @@ module Ao_basis : sig ao_prim_num : AO_prim_number.t array; ao_prim_num_max : AO_prim_number.t; ao_nucl : Nucl_number.t array; - ao_power : Symmetry.Xyz.t array; + ao_power : Angmom.Xyz.t array; ao_coef : AO_coef.t array; ao_expo : AO_expo.t array; ao_cartesian : bool; @@ -32,7 +32,7 @@ end = struct ao_prim_num : AO_prim_number.t array; ao_prim_num_max : AO_prim_number.t; ao_nucl : Nucl_number.t array; - ao_power : Symmetry.Xyz.t array; + ao_power : Angmom.Xyz.t array; ao_coef : AO_coef.t array; ao_expo : AO_expo.t array; ao_cartesian : bool; @@ -87,7 +87,7 @@ end = struct if (data.(2*dim+i-1) > 0) then result.(i-1) <- result.(i-1)^"z"^(string_of_int data.(2*dim+i-1)); done; - Array.map Symmetry.Xyz.of_string result + Array.map Angmom.Xyz.of_string result ;; let read_ao_coef () = @@ -133,7 +133,7 @@ end = struct let ao_num = AO_number.to_int b.ao_num in let gto_array = Array.init (AO_number.to_int b.ao_num) (fun i -> - let s = Symmetry.Xyz.to_symmetry b.ao_power.(i) in + let s = Angmom.Xyz.to_symmetry b.ao_power.(i) in let ao_prim_num = AO_prim_number.to_int b.ao_prim_num.(i) in let prims = List.init ao_prim_num (fun j -> let prim = { GaussianPrimitive.sym = s ; @@ -217,9 +217,9 @@ end = struct let ao_power = let l = Array.to_list ao_power in List.concat [ - (list_map (fun a -> Positive_int.to_int a.Symmetry.Xyz.x) l) ; - (list_map (fun a -> Positive_int.to_int a.Symmetry.Xyz.y) l) ; - (list_map (fun a -> Positive_int.to_int a.Symmetry.Xyz.z) l) ] + (list_map (fun a -> Positive_int.to_int a.Angmom.Xyz.x) l) ; + (list_map (fun a -> Positive_int.to_int a.Angmom.Xyz.y) l) ; + (list_map (fun a -> Positive_int.to_int a.Angmom.Xyz.z) l) ] in Ezfio.set_ao_basis_ao_power(Ezfio.ezfio_array_of_list ~rank:2 ~dim:[| ao_num ; 3 |] ~data:ao_power) ; @@ -409,7 +409,7 @@ end = struct | [] -> [] | (i,n,x)::tail -> (Printf.sprintf " %5d %6d %-8s\n" i (Nucl_number.to_int n) - (Symmetry.Xyz.to_string x) + (Angmom.Xyz.to_string x) )::(do_work tail) in do_work l |> String.concat "" @@ -496,7 +496,7 @@ md5 = %s (b.ao_nucl |> Array.to_list |> list_map Nucl_number.to_string |> String.concat ", ") (b.ao_power |> Array.to_list |> list_map (fun x-> - "("^(Symmetry.Xyz.to_string x)^")" )|> String.concat ", ") + "("^(Angmom.Xyz.to_string x)^")" )|> String.concat ", ") (b.ao_coef |> Array.to_list |> list_map AO_coef.to_string |> String.concat ", ") (b.ao_expo |> Array.to_list |> list_map AO_expo.to_string diff --git a/ocaml/Long_basis.ml b/ocaml/Long_basis.ml index 6c88954a..a8ea3c66 100644 --- a/ocaml/Long_basis.ml +++ b/ocaml/Long_basis.ml @@ -2,7 +2,7 @@ open Qptypes open Qputils open Sexplib.Std -type t = (Symmetry.Xyz.t * Gto.t * Nucl_number.t ) list [@@deriving sexp] +type t = (Angmom.Xyz.t * Gto.t * Nucl_number.t ) list [@@deriving sexp] let of_basis b = let rec do_work accu = function @@ -10,7 +10,7 @@ let of_basis b = | (g,n)::tail -> begin let new_accu = - Symmetry.Xyz.of_symmetry g.Gto.sym + Angmom.Xyz.of_symmetry g.Gto.sym |> List.rev_map (fun x-> (x,g,n)) in do_work (new_accu@accu) tail @@ -25,7 +25,7 @@ let to_basis b = | [] -> List.rev accu | (s,g,n)::tail -> let first_sym = - Symmetry.Xyz.of_symmetry g.Gto.sym + Angmom.Xyz.of_symmetry g.Gto.sym |> List.hd in let new_accu = @@ -42,7 +42,7 @@ let to_basis b = let to_string b = let middle = list_map (fun (x,y,z) -> "( "^((string_of_int (Nucl_number.to_int z)))^", "^ - (Symmetry.Xyz.to_string x)^", "^(Gto.to_string y) + (Angmom.Xyz.to_string x)^", "^(Gto.to_string y) ^" )" ) b |> String.concat ",\n" diff --git a/ocaml/Long_basis.mli b/ocaml/Long_basis.mli index 26009c10..2c41ad74 100644 --- a/ocaml/Long_basis.mli +++ b/ocaml/Long_basis.mli @@ -5,16 +5,16 @@ open Qptypes;; * all the D orbitals are converted to xx, xy, xz, yy, yx * etc *) -type t = (Symmetry.Xyz.t * Gto.t * Nucl_number.t) list [@@deriving sexp] +type t = (Angmom.Xyz.t * Gto.t * Nucl_number.t) list [@@deriving sexp] (** Transform a basis to a long basis *) val of_basis : - (Gto.t * Nucl_number.t) list -> (Symmetry.Xyz.t * Gto.t * Nucl_number.t) list + (Gto.t * Nucl_number.t) list -> (Angmom.Xyz.t * Gto.t * Nucl_number.t) list (** Transform a long basis to a basis *) val to_basis : - (Symmetry.Xyz.t * Gto.t * Nucl_number.t) list -> (Gto.t * Nucl_number.t) list + (Angmom.Xyz.t * Gto.t * Nucl_number.t) list -> (Gto.t * Nucl_number.t) list (** Convert the basis into its string representation *) val to_string : - (Symmetry.Xyz.t * Gto.t * Nucl_number.t) list -> string + (Angmom.Xyz.t * Gto.t * Nucl_number.t) list -> string diff --git a/ocaml/Primitive.mli b/ocaml/Primitive.mli index f7d8809d..f63ecd8c 100644 --- a/ocaml/Primitive.mli +++ b/ocaml/Primitive.mli @@ -1,5 +1,5 @@ type t = -{ sym : Symmetry.t; +{ sym : Angmom.t; expo : Qptypes.AO_expo.t; } [@@deriving sexp] @@ -7,5 +7,5 @@ type t = val to_string : t -> string (** Creation *) -val of_sym_expo : Symmetry.t -> Qptypes.AO_expo.t -> t +val of_sym_expo : Angmom.t -> Qptypes.AO_expo.t -> t diff --git a/ocaml/qp_create_ezfio.ml b/ocaml/qp_create_ezfio.ml index f4c7f3b6..6a4a36a1 100644 --- a/ocaml/qp_create_ezfio.ml +++ b/ocaml/qp_create_ezfio.ml @@ -512,13 +512,14 @@ let run ?o b au c d m p cart xyz_file = let ao_num = List.length long_basis in Ezfio.set_ao_basis_ao_num ao_num; Ezfio.set_ao_basis_ao_basis b; + Ezfio.set_basis_basis b; let ao_prim_num = list_map (fun (_,g,_) -> List.length g.Gto.lc) long_basis and ao_nucl = list_map (fun (_,_,n) -> Nucl_number.to_int n) long_basis and ao_power= let l = list_map (fun (x,_,_) -> x) long_basis in - (list_map (fun t -> Positive_int.to_int Symmetry.Xyz.(t.x)) l)@ - (list_map (fun t -> Positive_int.to_int Symmetry.Xyz.(t.y)) l)@ - (list_map (fun t -> Positive_int.to_int Symmetry.Xyz.(t.z)) l) + (list_map (fun t -> Positive_int.to_int Angmom.Xyz.(t.x)) l)@ + (list_map (fun t -> Positive_int.to_int Angmom.Xyz.(t.y)) l)@ + (list_map (fun t -> Positive_int.to_int Angmom.Xyz.(t.z)) l) in let ao_prim_num_max = List.fold_left (fun s x -> if x > s then x diff --git a/src/basis/EZFIO.cfg b/src/basis/EZFIO.cfg new file mode 100644 index 00000000..d560a291 --- /dev/null +++ b/src/basis/EZFIO.cfg @@ -0,0 +1,47 @@ +[basis] +type: character*(256) +doc: Name of the |AO| basis set +interface: ezfio + +[shell_num] +type: integer +doc: Number of shells +interface: ezfio, provider + +[shell_normalization_factor] +type: double precision +doc: Number of primitives per |AO| +size: (basis.shell_num) +interface: ezfio, provider + +[shell_prim_num] +type: integer +doc: Max number of primitives in a shell +size: (basis.shell_num) +interface: ezfio, provider + +[shell_prim_index] +type: integer +doc: Max number of primitives in a shell +size: (basis.shell_num) +interface: ezfio, provider + +[shell_nucl] +type: integer +doc: Index of the nucleus on which the shell is centered +size: (basis.shell_num) +interface: ezfio, provider + +[shell_coef] +type: double precision +doc: Primitive coefficients +size: (basis.shell_prim_num) +interface: ezfio, provider + +[shell_expo] +type: double precision +doc: Exponents in the shell +size: (basis.shell_prim_num) +interface: ezfio, provider + + diff --git a/src/basis/NEED b/src/basis/NEED new file mode 100644 index 00000000..d2066b18 --- /dev/null +++ b/src/basis/NEED @@ -0,0 +1 @@ +nuclei diff --git a/src/basis/README.rst b/src/basis/README.rst new file mode 100644 index 00000000..17493b4f --- /dev/null +++ b/src/basis/README.rst @@ -0,0 +1,8 @@ +====== +basis +====== + +This module contains the basis set information, which will then be used to build the atomic orbitals. + + + diff --git a/src/basis/basis.irp.f b/src/basis/basis.irp.f new file mode 100644 index 00000000..2674b9f5 --- /dev/null +++ b/src/basis/basis.irp.f @@ -0,0 +1,32 @@ +BEGIN_PROVIDER [ double precision, basis_normalization_factor, (shell_num) ] + implicit none + BEGIN_DOC + ! Normalization factors of the shells + END_DOC + double precision :: norm,overlap_x,overlap_y,overlap_z,C_A(3), c + integer :: l, powA(3), nz + integer :: i,j,k + nz=100 + C_A(1) = 0.d0 + C_A(2) = 0.d0 + C_A(3) = 0.d0 + + do i=1,shell_num + powA(1) = shell_ang_mom(i) + powA(2) = 0 + powA(3) = 0 + + ! Normalization of the contracted basis functions + norm = 0.d0 + do j=shell_prim_index(i), shell_prim_index(i)+shell_prim_num(i)-1 + do k=shell_prim_index(i), shell_prim_index(i)+shell_prim_num(i)-1 + call overlap_gaussian_xyz(C_A,C_A,shell_prim_expo(j),shell_prim_expo(k), & + powA,powA,overlap_x,overlap_y,overlap_z,c,nz) + norm = norm+c*shell_prim_coef(j)*shell_prim_coef(k) + enddo + enddo + basis_normalization_factor(i) = basis_normalization_factor(i) * sqrt(norm) + + enddo + +END_PROVIDER From 1685e412584b4ee521cfb33571547b5fc52c037a Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 21 May 2021 18:26:50 +0200 Subject: [PATCH 08/29] Added basis module for easy export --- data/basis/cc-pvdz | 2 +- ocaml/qp_create_ezfio.ml | 51 ++++++++++++++++++++++++++++++++++++++++ src/ao_basis/NEED | 1 + src/ao_basis/aos.irp.f | 33 -------------------------- src/basis/EZFIO.cfg | 28 +++++++++++++++++----- 5 files changed, 75 insertions(+), 40 deletions(-) diff --git a/data/basis/cc-pvdz b/data/basis/cc-pvdz index 3d1a77e6..156e46ac 100644 --- a/data/basis/cc-pvdz +++ b/data/basis/cc-pvdz @@ -3606,4 +3606,4 @@ D 5 5 1.507524E+00 2.667560E-01 D 1 1 5.030000E-01 1.000000E+00 -$END \ No newline at end of file +$END diff --git a/ocaml/qp_create_ezfio.ml b/ocaml/qp_create_ezfio.ml index 6a4a36a1..bf3018b3 100644 --- a/ocaml/qp_create_ezfio.ml +++ b/ocaml/qp_create_ezfio.ml @@ -562,6 +562,57 @@ let run ?o b au c d m p cart xyz_file = and ao_expo = create_expo_coef `Expos in let () = + let shell_num = List.length basis in + let lc : (GaussianPrimitive.t * Qptypes.AO_coef.t) list list = + list_map ( fun (g,_) -> g.Gto.lc ) basis + in + let ang_mom = + list_map (fun (l : (GaussianPrimitive.t * Qptypes.AO_coef.t) list) -> + let x, _ = List.hd l in + Angmom.to_l x.GaussianPrimitive.sym |> Qptypes.Positive_int.to_int + ) lc + in + let expo = + list_map (fun l -> list_map (fun (x,_) -> Qptypes.AO_expo.to_float x.GaussianPrimitive.expo) l ) lc + |> List.concat + in + let coef = + list_map (fun l -> + list_map (fun (_,x) -> Qptypes.AO_coef.to_float x) l + ) lc + |> List.concat + in + let shell_prim_num = + list_map List.length lc + in + let shell_prim_idx = + let rec aux count accu = function + | [] -> List.rev accu + | l::rest -> + let newcount = count+(List.length l) in + aux newcount (count::accu) rest + in + aux 1 [] lc + in + let prim_num = List.length coef in + Ezfio.set_basis_typ "Gaussian"; + Ezfio.set_basis_shell_num shell_num; + Ezfio.set_basis_prim_num prim_num ; + Ezfio.set_basis_shell_prim_num (Ezfio.ezfio_array_of_list + ~rank:1 ~dim:[| shell_num |] ~data:shell_prim_num); + Ezfio.set_basis_shell_ang_mom (Ezfio.ezfio_array_of_list + ~rank:1 ~dim:[| shell_num |] ~data:ang_mom ) ; + Ezfio.set_basis_shell_prim_index (Ezfio.ezfio_array_of_list + ~rank:1 ~dim:[| shell_num |] ~data:shell_prim_idx) ; + Ezfio.set_basis_shell_nucl (Ezfio.ezfio_array_of_list + ~rank:1 ~dim:[| shell_num |] + ~data:(list_map (fun (_,n) -> Nucl_number.to_int n) basis) ) ; + Ezfio.set_basis_shell_prim_coef (Ezfio.ezfio_array_of_list + ~rank:1 ~dim:[| prim_num |] ~data:coef) ; + Ezfio.set_basis_shell_prim_expo (Ezfio.ezfio_array_of_list + ~rank:1 ~dim:[| prim_num |] ~data:expo) ; + + Ezfio.set_ao_basis_ao_prim_num (Ezfio.ezfio_array_of_list ~rank:1 ~dim:[| ao_num |] ~data:ao_prim_num) ; Ezfio.set_ao_basis_ao_nucl(Ezfio.ezfio_array_of_list diff --git a/src/ao_basis/NEED b/src/ao_basis/NEED index d2066b18..a69fe768 100644 --- a/src/ao_basis/NEED +++ b/src/ao_basis/NEED @@ -1 +1,2 @@ nuclei +basis diff --git a/src/ao_basis/aos.irp.f b/src/ao_basis/aos.irp.f index e95ac711..a8384ad4 100644 --- a/src/ao_basis/aos.irp.f +++ b/src/ao_basis/aos.irp.f @@ -56,39 +56,6 @@ END_PROVIDER END_PROVIDER -BEGIN_PROVIDER [ double precision, ao_coef_normalization_libint_factor, (ao_num) ] - implicit none - BEGIN_DOC - ! |AO| normalization for interfacing with libint - END_DOC - double precision :: norm,overlap_x,overlap_y,overlap_z,C_A(3), c - integer :: l, powA(3), nz - integer :: i,j,k - nz=100 - C_A(1) = 0.d0 - C_A(2) = 0.d0 - C_A(3) = 0.d0 - - do i=1,ao_num - powA(1) = ao_l(i) - powA(2) = 0 - powA(3) = 0 - - ! Normalization of the contracted basis functions - norm = 0.d0 - do j=1,ao_prim_num(i) - do k=1,ao_prim_num(i) - call overlap_gaussian_xyz(C_A,C_A,ao_expo(i,j),ao_expo(i,k),powA,powA,overlap_x,overlap_y,overlap_z,c,nz) - norm = norm+c*ao_coef_normalized(i,j)*ao_coef_normalized(i,k) - enddo - enddo - ao_coef_normalization_libint_factor(i) = ao_coef_normalization_factor(i) * sqrt(norm) - - enddo - -END_PROVIDER - - BEGIN_PROVIDER [ double precision, ao_coef_normalized_ordered, (ao_num,ao_prim_num_max) ] &BEGIN_PROVIDER [ double precision, ao_expo_ordered, (ao_num,ao_prim_num_max) ] implicit none diff --git a/src/basis/EZFIO.cfg b/src/basis/EZFIO.cfg index d560a291..57bd427b 100644 --- a/src/basis/EZFIO.cfg +++ b/src/basis/EZFIO.cfg @@ -1,7 +1,12 @@ [basis] type: character*(256) doc: Name of the |AO| basis set -interface: ezfio +interface: ezfio, provider + +[typ] +type: character*(32) +doc: Type of basis set. Only 'Gaussian' is supported +interface: ezfio, provider [shell_num] type: integer @@ -14,9 +19,20 @@ doc: Number of primitives per |AO| size: (basis.shell_num) interface: ezfio, provider +[prim_num] +type: integer +doc: Total number of primitives +interface: ezfio, provider + +[shell_ang_mom] +type: integer +doc: Angular momentum of each shell +size: (basis.shell_num) +interface: ezfio, provider + [shell_prim_num] type: integer -doc: Max number of primitives in a shell +doc: Number of primitives in a shell size: (basis.shell_num) interface: ezfio, provider @@ -32,16 +48,16 @@ doc: Index of the nucleus on which the shell is centered size: (basis.shell_num) interface: ezfio, provider -[shell_coef] +[shell_prim_coef] type: double precision doc: Primitive coefficients -size: (basis.shell_prim_num) +size: (basis.prim_num) interface: ezfio, provider -[shell_expo] +[shell_prim_expo] type: double precision doc: Exponents in the shell -size: (basis.shell_prim_num) +size: (basis.prim_num) interface: ezfio, provider From e7dd374ab9bbdb045537272f4126a51d62bdeaa0 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 24 May 2021 21:55:14 +0200 Subject: [PATCH 09/29] Removed omp_set_nested --- src/ao_basis/aos.irp.f | 49 ++++++++----------- src/basis/EZFIO.cfg | 2 +- src/basis/basis.irp.f | 78 +++++++++++++++++++++--------- src/cipsi/pt2_stoch_routines.irp.f | 3 +- src/dressing/run_dress_slave.irp.f | 6 +-- 5 files changed, 80 insertions(+), 58 deletions(-) diff --git a/src/ao_basis/aos.irp.f b/src/ao_basis/aos.irp.f index a8384ad4..2ff72898 100644 --- a/src/ao_basis/aos.irp.f +++ b/src/ao_basis/aos.irp.f @@ -6,6 +6,23 @@ BEGIN_PROVIDER [ integer, ao_prim_num_max ] ao_prim_num_max = maxval(ao_prim_num) END_PROVIDER +BEGIN_PROVIDER [ integer, ao_shell, (ao_num) ] + implicit none + BEGIN_DOC + ! Index of the shell to which the AO corresponds + END_DOC + integer :: i, j, k, n + k=0 + do i=1,shell_num + n = shell_ang_mom(i)+1 + do j=1,(n*(n+1))/2 + k = k+1 + ao_shell(k) = i + enddo + enddo + +END_PROVIDER + BEGIN_PROVIDER [ double precision, ao_coef_normalized, (ao_num,ao_prim_num_max) ] &BEGIN_PROVIDER [ double precision, ao_coef_normalization_factor, (ao_num) ] implicit none @@ -30,7 +47,8 @@ END_PROVIDER ! Normalization of the primitives if (primitives_normalized) then do j=1,ao_prim_num(i) - call overlap_gaussian_xyz(C_A,C_A,ao_expo(i,j),ao_expo(i,j),powA,powA,overlap_x,overlap_y,overlap_z,norm,nz) + call overlap_gaussian_xyz(C_A,C_A,ao_expo(i,j),ao_expo(i,j), & + powA,powA,overlap_x,overlap_y,overlap_z,norm,nz) ao_coef_normalized(i,j) = ao_coef(i,j)/sqrt(norm) enddo else @@ -198,38 +216,11 @@ END_PROVIDER END_DOC do i = 1, nucl_num Nucl_num_shell_Aos(i) = 0 - do j = 1, Nucl_N_Aos(i) - if(ao_l(Nucl_Aos(i,j))==0)then - ! S type function - Nucl_num_shell_Aos(i)+=1 - Nucl_list_shell_Aos(i,Nucl_num_shell_Aos(i))=Nucl_Aos(i,j) - elseif(ao_l(Nucl_Aos(i,j))==1)then - ! P type function - if(ao_power(Nucl_Aos(i,j),1)==1)then + if (ao_power(Nucl_Aos(i,j),1) == ao_l(Nucl_Aos(i,j))) then Nucl_num_shell_Aos(i)+=1 Nucl_list_shell_Aos(i,Nucl_num_shell_Aos(i))=Nucl_Aos(i,j) endif - elseif(ao_l(Nucl_Aos(i,j))==2)then - ! D type function - if(ao_power(Nucl_Aos(i,j),1)==2)then - Nucl_num_shell_Aos(i)+=1 - Nucl_list_shell_Aos(i,Nucl_num_shell_Aos(i))=Nucl_Aos(i,j) - endif - elseif(ao_l(Nucl_Aos(i,j))==3)then - ! F type function - if(ao_power(Nucl_Aos(i,j),1)==3)then - Nucl_num_shell_Aos(i)+=1 - Nucl_list_shell_Aos(i,Nucl_num_shell_Aos(i))=Nucl_Aos(i,j) - endif - elseif(ao_l(Nucl_Aos(i,j))==4)then - ! G type function - if(ao_power(Nucl_Aos(i,j),1)==4)then - Nucl_num_shell_Aos(i)+=1 - Nucl_list_shell_Aos(i,Nucl_num_shell_Aos(i))=Nucl_Aos(i,j) - endif - endif - enddo enddo diff --git a/src/basis/EZFIO.cfg b/src/basis/EZFIO.cfg index 57bd427b..1c66e758 100644 --- a/src/basis/EZFIO.cfg +++ b/src/basis/EZFIO.cfg @@ -17,7 +17,7 @@ interface: ezfio, provider type: double precision doc: Number of primitives per |AO| size: (basis.shell_num) -interface: ezfio, provider +interface: ezfio [prim_num] type: integer diff --git a/src/basis/basis.irp.f b/src/basis/basis.irp.f index 2674b9f5..e25e6a46 100644 --- a/src/basis/basis.irp.f +++ b/src/basis/basis.irp.f @@ -1,32 +1,62 @@ -BEGIN_PROVIDER [ double precision, basis_normalization_factor, (shell_num) ] +BEGIN_PROVIDER [ double precision, shell_normalization_factor , (shell_num) ] implicit none BEGIN_DOC - ! Normalization factors of the shells + ! Number of primitives per |AO| END_DOC - double precision :: norm,overlap_x,overlap_y,overlap_z,C_A(3), c - integer :: l, powA(3), nz - integer :: i,j,k - nz=100 - C_A(1) = 0.d0 - C_A(2) = 0.d0 - C_A(3) = 0.d0 - do i=1,shell_num - powA(1) = shell_ang_mom(i) - powA(2) = 0 - powA(3) = 0 + logical :: has + PROVIDE ezfio_filename + if (mpi_master) then + if (size(shell_normalization_factor) == 0) return - ! Normalization of the contracted basis functions - norm = 0.d0 - do j=shell_prim_index(i), shell_prim_index(i)+shell_prim_num(i)-1 - do k=shell_prim_index(i), shell_prim_index(i)+shell_prim_num(i)-1 - call overlap_gaussian_xyz(C_A,C_A,shell_prim_expo(j),shell_prim_expo(k), & - powA,powA,overlap_x,overlap_y,overlap_z,c,nz) - norm = norm+c*shell_prim_coef(j)*shell_prim_coef(k) - enddo - enddo - basis_normalization_factor(i) = basis_normalization_factor(i) * sqrt(norm) + call ezfio_has_basis_shell_normalization_factor(has) + if (has) then + write(6,'(A)') '.. >>>>> [ IO READ: shell_normalization_factor ] <<<<< ..' + call ezfio_get_basis_shell_normalization_factor(shell_normalization_factor) + else - enddo + + double precision :: norm,overlap_x,overlap_y,overlap_z,C_A(3), c + integer :: l, powA(3), nz + integer :: i,j,k + nz=100 + C_A(1) = 0.d0 + C_A(2) = 0.d0 + C_A(3) = 0.d0 + + do i=1,shell_num + powA(1) = shell_ang_mom(i) + powA(2) = 0 + powA(3) = 0 + + ! Normalization of the contracted basis functions + norm = 0.d0 + do j=shell_prim_index(i), shell_prim_index(i)+shell_prim_num(i)-1 + do k=shell_prim_index(i), shell_prim_index(i)+shell_prim_num(i)-1 + call overlap_gaussian_xyz(C_A,C_A,shell_prim_expo(j),shell_prim_expo(k),& + powA,powA,overlap_x,overlap_y,overlap_z,c,nz) + norm = norm+c*shell_prim_coef(j)*shell_prim_coef(k) + enddo + enddo + shell_normalization_factor(i) = dsqrt(norm) + enddo + + + endif + 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( shell_normalization_factor, (shell_num), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + if (ierr /= MPI_SUCCESS) then + stop 'Unable to read shell_normalization_factor with MPI' + endif + IRP_ENDIF + + call write_time(6) END_PROVIDER diff --git a/src/cipsi/pt2_stoch_routines.irp.f b/src/cipsi/pt2_stoch_routines.irp.f index 7554c39e..6f0d6683 100644 --- a/src/cipsi/pt2_stoch_routines.irp.f +++ b/src/cipsi/pt2_stoch_routines.irp.f @@ -286,7 +286,7 @@ 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_nested(.false.) + call omp_set_max_active_levels(1) print '(A)', '========== ======================= ===================== ===================== ===========' @@ -313,6 +313,7 @@ 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) print '(A)', '========== ======================= ===================== ===================== ===========' diff --git a/src/dressing/run_dress_slave.irp.f b/src/dressing/run_dress_slave.irp.f index 8a92962c..a33fb1dd 100644 --- a/src/dressing/run_dress_slave.irp.f +++ b/src/dressing/run_dress_slave.irp.f @@ -72,7 +72,7 @@ subroutine run_dress_slave(thread,iproce,energy) provide psi_energy ending = dress_N_cp+1 ntask_tbd = 0 - call omp_set_nested(.true.) + call omp_set_max_active_levels(8) !$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_nested(.false.) + call omp_set_max_active_levels(1) 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_nested(.false.) + call omp_set_max_active_levels(1) ! do i=0,dress_N_cp+1 ! call omp_destroy_lock(lck_sto(i)) ! end do From 94a069c2d5be74f2af87f49105a22a6460dc3e3c Mon Sep 17 00:00:00 2001 From: vijay gopal chilkuri Date: Wed, 26 May 2021 20:11:39 +0530 Subject: [PATCH 10/29] Simplified calculation of n_CSF. #158 --- src/csf/sigma_vector.irp.f | 33 +++++++++++++++------------------ 1 file changed, 15 insertions(+), 18 deletions(-) diff --git a/src/csf/sigma_vector.irp.f b/src/csf/sigma_vector.irp.f index 265d2384..1d6e663a 100644 --- a/src/csf/sigma_vector.irp.f +++ b/src/csf/sigma_vector.irp.f @@ -28,25 +28,22 @@ detDimperBF = 0 MS = elec_alpha_num-elec_beta_num ! number of cfgs = number of dets for 0 somos - n_CSF = cfg_seniority_index(0)-1 + n_CSF = 0 ncfgprev = cfg_seniority_index(0) - do i = 0-iand(MS,1)+2, NSOMOMax,2 - if(cfg_seniority_index(i) .EQ. -1)then - ncfgpersomo = N_configuration + 1 - else - ncfgpersomo = cfg_seniority_index(i) - endif - ncfg = ncfgpersomo - ncfgprev - !detDimperBF = max(1,nint((binom(i,(i+1)/2)))) - if (i > 2) then - dimcsfpercfg = max(1,nint((binom(i-2,(i-2+1)/2)-binom(i-2,((i-2+1)/2)+1)))) - else - dimcsfpercfg = 1 - endif - n_CSF += ncfg * dimcsfpercfg - !if(cfg_seniority_index(i+2) == -1) EXIT - !if(detDimperBF > maxDetDimPerBF) maxDetDimPerBF = detDimperBF - ncfgprev = cfg_seniority_index(i) + do i = iand(MS,1), NSOMOMax-2,2 + if(cfg_seniority_index(i+2) .EQ. -1) then + ncfgpersomo = N_configuration + 1 + else + ncfgpersomo = cfg_seniority_index(i+2) + endif + ncfg = ncfgpersomo - ncfgprev + if(iand(MS,1) .EQ. 0) then + dimcsfpercfg = max(1,nint((binom(i,i/2)-binom(i,i/2+1)))) + else + dimcsfpercfg = max(1,nint((binom(i,(i+1)/2)-binom(i,(i+3)/2)))) + endif + n_CSF += ncfg * dimcsfpercfg + ncfgprev = cfg_seniority_index(i+2) enddo END_PROVIDER From 2557c768450b9bcabef021a8e1b8c2fb6db4fee6 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 27 May 2021 12:28:36 +0200 Subject: [PATCH 11/29] Travis.com --- README.md | 4 ++-- external/irpf90 | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/README.md b/README.md index d5496f5e..b8141b88 100644 --- a/README.md +++ b/README.md @@ -36,8 +36,8 @@ https://arxiv.org/abs/1902.08154 # Build status -* Master [![master build status](https://travis-ci.org/QuantumPackage/qp2.svg?branch=master)](https://travis-ci.org/QuantumPackage/qp2) -* Development [![dev build status](https://travis-ci.org/QuantumPackage/qp2.svg?branch=dev)](https://travis-ci.org/QuantumPackage/qp2) +* Master [![master build status](https://travis-ci.com/QuantumPackage/qp2.svg?branch=master)](https://travis-ci.org/QuantumPackage/qp2) +* Development [![dev build status](https://travis-ci.com/QuantumPackage/qp2.svg?branch=dev)](https://travis-ci.org/QuantumPackage/qp2) * Documentation [![Documentation Status](https://readthedocs.org/projects/quantum-package/badge/?version=master)](https://quantum-package.readthedocs.io/en/master/?badge=master) diff --git a/external/irpf90 b/external/irpf90 index 132a4a16..33ca5e10 160000 --- a/external/irpf90 +++ b/external/irpf90 @@ -1 +1 @@ -Subproject commit 132a4a1661c9878d21dcbf0ac14f7fe9a3b110d0 +Subproject commit 33ca5e1018f3bbb5e695e6ee558f5dac0753b271 From 32e2afca90dac6d12a6614fdbec39a03ad4c8f8a Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 31 May 2021 01:48:34 +0200 Subject: [PATCH 12/29] Using Intel IPP for sorting --- config/ifort.cfg | 4 +- config/ifort_avx.cfg | 4 +- config/ifort_avx_mpi.cfg | 4 +- config/ifort_debug.cfg | 4 +- config/ifort_mpi.cfg | 4 +- config/ifort_rome.cfg | 4 +- config/ifort_xHost.cfg | 6 +- src/ao_one_e_ints/pseudopot.f90 | 19 ++- src/cipsi/selection.irp.f | 7 +- src/csf/sigma_vector.irp.f | 106 ++++++------ .../diagonalization_hcsf_dressed.irp.f | 1 + src/utils/intel.f90 | 70 ++++++++ src/utils/sort.irp.f | 151 ++++++++++++++++-- 13 files changed, 297 insertions(+), 87 deletions(-) create mode 100644 src/utils/intel.f90 diff --git a/config/ifort.cfg b/config/ifort.cfg index 866aae3d..714c4b10 100644 --- a/config/ifort.cfg +++ b/config/ifort.cfg @@ -7,9 +7,9 @@ # [COMMON] FC : ifort -fpic -LAPACK_LIB : -mkl=parallel +LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps IRPF90 : irpf90 -IRPF90_FLAGS : --ninja --align=32 +IRPF90_FLAGS : --ninja --align=32 -DINTEL # Global options ################ diff --git a/config/ifort_avx.cfg b/config/ifort_avx.cfg index d689050e..a2cb4c8a 100644 --- a/config/ifort_avx.cfg +++ b/config/ifort_avx.cfg @@ -7,9 +7,9 @@ # [COMMON] FC : ifort -fpic -LAPACK_LIB : -mkl=parallel +LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps IRPF90 : irpf90 -IRPF90_FLAGS : --ninja --align=32 +IRPF90_FLAGS : --ninja --align=32 -DINTEL # Global options ################ diff --git a/config/ifort_avx_mpi.cfg b/config/ifort_avx_mpi.cfg index 9a839e66..f2bb8889 100644 --- a/config/ifort_avx_mpi.cfg +++ b/config/ifort_avx_mpi.cfg @@ -7,9 +7,9 @@ # [COMMON] FC : mpiifort -fpic -LAPACK_LIB : -mkl=parallel +LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps IRPF90 : irpf90 -IRPF90_FLAGS : --ninja --align=32 -DMPI +IRPF90_FLAGS : --ninja --align=32 -DMPI -DINTEL # Global options ################ diff --git a/config/ifort_debug.cfg b/config/ifort_debug.cfg index f2fbf8ce..9b718380 100644 --- a/config/ifort_debug.cfg +++ b/config/ifort_debug.cfg @@ -7,9 +7,9 @@ # [COMMON] FC : ifort -fpic -LAPACK_LIB : -mkl=parallel +LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps IRPF90 : irpf90 -IRPF90_FLAGS : --ninja --align=32 --assert +IRPF90_FLAGS : --ninja --align=32 --assert -DINTEL # Global options ################ diff --git a/config/ifort_mpi.cfg b/config/ifort_mpi.cfg index 57087847..e0d489a0 100644 --- a/config/ifort_mpi.cfg +++ b/config/ifort_mpi.cfg @@ -7,9 +7,9 @@ # [COMMON] FC : mpiifort -fpic -LAPACK_LIB : -mkl=parallel +LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps IRPF90 : irpf90 -IRPF90_FLAGS : --ninja --align=32 -DMPI +IRPF90_FLAGS : --ninja --align=32 -DMPI -DINTEL # Global options ################ diff --git a/config/ifort_rome.cfg b/config/ifort_rome.cfg index a4bce680..3b70f98b 100644 --- a/config/ifort_rome.cfg +++ b/config/ifort_rome.cfg @@ -7,9 +7,9 @@ # [COMMON] FC : ifort -fpic -LAPACK_LIB : -mkl=parallel +LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps IRPF90 : irpf90 -IRPF90_FLAGS : --ninja --align=32 +IRPF90_FLAGS : --ninja --align=32 -DINTEL # Global options ################ diff --git a/config/ifort_xHost.cfg b/config/ifort_xHost.cfg index 5d952e54..ddb4aa2d 100644 --- a/config/ifort_xHost.cfg +++ b/config/ifort_xHost.cfg @@ -6,10 +6,10 @@ # --align=32 : Align all provided arrays on a 32-byte boundary # [COMMON] -FC : ifort -fpic -LAPACK_LIB : -mkl=parallel +FC : ifort -fpic +LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps IRPF90 : irpf90 -IRPF90_FLAGS : --ninja --align=64 +IRPF90_FLAGS : --ninja --align=64 -DINTEL # Global options ################ diff --git a/src/ao_one_e_ints/pseudopot.f90 b/src/ao_one_e_ints/pseudopot.f90 index 6c8e5c83..48e3803e 100644 --- a/src/ao_one_e_ints/pseudopot.f90 +++ b/src/ao_one_e_ints/pseudopot.f90 @@ -96,8 +96,12 @@ end ! x=cos(theta) double precision function ylm_real(l,m,x,phi) - implicit double precision (a-h,o-z) - DIMENSION PM(0:100,0:100) + implicit none + integer :: MM, iabs_m, m, l + double precision :: pi, fourpi, factor, x, phi, coef + double precision :: xchap, ychap, zchap + double precision, external :: fact + double precision :: PM(0:100,0:100), plm MM=100 pi=dacos(-1.d0) fourpi=4.d0*pi @@ -1150,8 +1154,10 @@ end ! Output: PM(m,n) --- Pmn(x) ! ===================================================== ! - IMPLICIT DOUBLE PRECISION (P,X) - DIMENSION PM(0:MM,0:(N+1)) + implicit none +! IMPLICIT DOUBLE PRECISION (P,X) + integer :: MM, N, I, J, M + double precision :: PM(0:MM,0:(N+1)), X, XQ, XS DOUBLE PRECISION, SAVE :: INVERSE(100) = 0.D0 DOUBLE PRECISION :: LS, II, JJ IF (INVERSE(1) == 0.d0) THEN @@ -1202,8 +1208,9 @@ end ! P_l^|m|(cos(theta)) exp(i m phi) subroutine erreur(x,n,rmoy,error) - implicit double precision(a-h,o-z) - dimension x(n) + implicit none + integer :: i, n + double precision :: x(n), rn, rn1, error, rmoy ! calcul de la moyenne rmoy=0.d0 do i=1,n diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f index 11492652..eda9642c 100644 --- a/src/cipsi/selection.irp.f +++ b/src/cipsi/selection.irp.f @@ -253,12 +253,7 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d deallocate(exc_degree) nmax=k-1 - allocate(iorder(nmax)) - do i=1,nmax - iorder(i) = i - enddo - call isort(indices,iorder,nmax) - deallocate(iorder) + call isort_noidx(indices,nmax) ! Start with 32 elements. Size will double along with the filtering. allocate(preinteresting(0:32), prefullinteresting(0:32), & diff --git a/src/csf/sigma_vector.irp.f b/src/csf/sigma_vector.irp.f index 265d2384..0f3f093c 100644 --- a/src/csf/sigma_vector.irp.f +++ b/src/csf/sigma_vector.irp.f @@ -1,54 +1,62 @@ - BEGIN_PROVIDER [ integer, NSOMOMax] - &BEGIN_PROVIDER [ integer, NCSFMax] - &BEGIN_PROVIDER [ integer*8, NMO] - &BEGIN_PROVIDER [ integer, NBFMax] - &BEGIN_PROVIDER [ integer, n_CSF] - &BEGIN_PROVIDER [ integer, maxDetDimPerBF] - implicit none - BEGIN_DOC - ! Documentation for NSOMOMax - ! The maximum number of SOMOs for the current calculation. - ! required for the calculation of prototype arrays. - END_DOC - NSOMOMax = min(elec_num, cfg_nsomo_max + 2) - ! Note that here we need NSOMOMax + 2 sizes - NCSFMax = max(1,nint((binom(NSOMOMax,(NSOMOMax+1)/2)-binom(NSOMOMax,((NSOMOMax+1)/2)+1)))) ! TODO: NCSFs for MS=0 - NBFMax = NCSFMax - maxDetDimPerBF = max(1,nint((binom(NSOMOMax,(NSOMOMax+1)/2)))) - NMO = n_act_orb - integer i,j,k,l - integer startdet,enddet - integer ncfg,ncfgprev - integer NSOMO - integer dimcsfpercfg - integer detDimperBF - real*8 :: coeff - integer MS - integer ncfgpersomo - detDimperBF = 0 - MS = elec_alpha_num-elec_beta_num - ! number of cfgs = number of dets for 0 somos - n_CSF = cfg_seniority_index(0)-1 - ncfgprev = cfg_seniority_index(0) - do i = 0-iand(MS,1)+2, NSOMOMax,2 - if(cfg_seniority_index(i) .EQ. -1)then - ncfgpersomo = N_configuration + 1 + BEGIN_PROVIDER [ integer, NSOMOMax] +&BEGIN_PROVIDER [ integer, NCSFMax] +&BEGIN_PROVIDER [ integer*8, NMO] +&BEGIN_PROVIDER [ integer, NBFMax] +&BEGIN_PROVIDER [ integer, n_CSF] +&BEGIN_PROVIDER [ integer, maxDetDimPerBF] + implicit none + + BEGIN_DOC + ! The maximum number of SOMOs for the current calculation. + ! required for the calculation of prototype arrays. + END_DOC + + integer :: i,j,k,l + integer :: startdet,enddet + integer :: ncfg,ncfgprev + integer :: NSOMO + integer :: dimcsfpercfg + integer :: detDimperBF + real*8 :: coeff + integer :: MS + integer :: ncfgpersomo + double precision :: bin_1, bin_2 + + NSOMOMax = min(elec_num, cfg_nsomo_max + 2) + + bin_1 = binom(NSOMOMax, (NSOMOMax+1)/2) + bin_2 = binom(NSOMOMax,((NSOMOMax+1)/2)+1) + + ! Note that here we need NSOMOMax + 2 sizes + NCSFMax = max(1,int(bin_1-bin_2)) + NBFMax = NCSFMax + maxDetDimPerBF = max(1,nint((binom(NSOMOMax,(NSOMOMax+1)/2)))) + NMO = n_act_orb + + detDimperBF = 0 + MS = elec_alpha_num-elec_beta_num + + ! number of cfgs = number of dets for 0 somos + n_CSF = cfg_seniority_index(0)-1 + ncfgprev = cfg_seniority_index(0) + do i = 0-iand(MS,1)+2, cfg_nsomo_max,2 + if(cfg_seniority_index(i) == -1)then + ncfgpersomo = N_configuration + 1 else - ncfgpersomo = cfg_seniority_index(i) + ncfgpersomo = cfg_seniority_index(i) endif - ncfg = ncfgpersomo - ncfgprev - !detDimperBF = max(1,nint((binom(i,(i+1)/2)))) - if (i > 2) then - dimcsfpercfg = max(1,nint((binom(i-2,(i-2+1)/2)-binom(i-2,((i-2+1)/2)+1)))) - else - dimcsfpercfg = 1 - endif - n_CSF += ncfg * dimcsfpercfg - !if(cfg_seniority_index(i+2) == -1) EXIT - !if(detDimperBF > maxDetDimPerBF) maxDetDimPerBF = detDimperBF - ncfgprev = cfg_seniority_index(i) - enddo - END_PROVIDER + ncfg = ncfgpersomo - ncfgprev + if (i > 2) then + bin_1 = binom(i-2, (i-2+1)/2) + bin_2 = binom(i-2,((i-2+1)/2)+1) + dimcsfpercfg = max(1,int(bin_1-bin_2)) + else + dimcsfpercfg = 1 + endif + n_CSF += ncfg * dimcsfpercfg + ncfgprev = cfg_seniority_index(i) + enddo +END_PROVIDER subroutine get_phase_qp_to_cfg(Ialpha, Ibeta, phaseout) diff --git a/src/davidson/diagonalization_hcsf_dressed.irp.f b/src/davidson/diagonalization_hcsf_dressed.irp.f index 2a83cc28..da23b919 100644 --- a/src/davidson/diagonalization_hcsf_dressed.irp.f +++ b/src/davidson/diagonalization_hcsf_dressed.irp.f @@ -197,6 +197,7 @@ subroutine davidson_diag_csf_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N call write_int(6,N_st,'Number of states') call write_int(6,N_st_diag,'Number of states in diagonalization') call write_int(6,sze,'Number of determinants') + call write_int(6,sze_csf,'Number of CSFs') call write_int(6,nproc_target,'Number of threads for diagonalization') call write_double(6, r1, 'Memory(Gb)') if (disk_based) then diff --git a/src/utils/intel.f90 b/src/utils/intel.f90 new file mode 100644 index 00000000..681adde9 --- /dev/null +++ b/src/utils/intel.f90 @@ -0,0 +1,70 @@ +module intel + use, intrinsic :: iso_c_binding + interface + subroutine ippsSortRadixIndexGetBufferSize(len, dataType, pBufSize) bind(C, name='ippsSortRadixIndexGetBufferSize') + use iso_c_binding + integer, intent(in), value :: len + integer, intent(in), value :: dataType + integer, intent(out) :: pBufSize + end + end interface + interface + subroutine ippsSortAscend_32s_I(pSrc, len) bind(C, name='ippsSortAscend_32s_I') + use iso_c_binding + integer, intent(in), value :: len + integer, intent(inout) :: pSrc(len) + end + end interface + interface + subroutine ippsSortRadixAscend_32s_I(pSrc, len, pTmp) bind(C, name='ippsSortRadixAscend_32s_I') + use iso_c_binding + integer, intent(in), value :: len + integer, intent(inout) :: pSrc(len) + integer, intent(inout) :: pTmp(len) + end + end interface + interface + subroutine ippsSortRadixIndexAscend_32s(pSrc, srcStrideBytes, pDstIndx, len, pTmpIndx) bind(C, name='ippsSortRadixIndexAscend_32s') + use iso_c_binding + integer, intent(in), value :: len + integer, intent(inout) :: pSrc(len) + integer, intent(in), value :: srcStrideBytes + integer, intent(inout) :: pDstIndx(len) + integer, intent(inout) :: pTmpIndx(len) + end + end interface + interface + subroutine ippsSortRadixIndexAscend_32f(pSrc, srcStrideBytes, pDstIndx, len, pTmpIndx) bind(C,name='ippsSortRadixIndexAscend_32f') + use iso_c_binding + integer, intent(in), value :: len + real , intent(inout) :: pSrc(len) + integer, intent(in), value :: srcStrideBytes + integer, intent(inout) :: pDstIndx(len) + integer, intent(inout) :: pTmpIndx(len) + end + end interface + interface + subroutine ippsSortIndexAscend_32f_I(pSrcDst, pDstIndx, len) bind(C,name='ippsSortIndexAscend_32f_I') + use iso_c_binding + real(4), intent(in) :: pSrcDst(*) + integer(4), intent(inout) :: pDstIndx(*) + integer(4), intent(in), value :: len + end + end interface + interface + subroutine ippsSortIndexAscend_32s_I(pSrcDst, pDstIndx, len) bind(C,name='ippsSortIndexAscend_32s_I') + use iso_c_binding + integer(4), intent(in) :: pSrcDst(*) + integer(4), intent(inout) :: pDstIndx(*) + integer(4), intent(in), value :: len + end + end interface + interface + subroutine ippsSortIndexAscend_64f_I(pSrcDst, pDstIndx, len) bind(C,name='ippsSortIndexAscend_64f_I') + use iso_c_binding + real(8), intent(in) :: pSrcDst(*) + integer(4), intent(inout) :: pDstIndx(*) + integer(4), intent(in), value :: len + end + end interface +end module diff --git a/src/utils/sort.irp.f b/src/utils/sort.irp.f index 2a655eed..be58e7a5 100644 --- a/src/utils/sort.irp.f +++ b/src/utils/sort.irp.f @@ -57,7 +57,7 @@ BEGIN_TEMPLATE $type :: c, tmp integer :: itmp integer :: i, j - + if(isize<2)return c = x( shiftr(first+last,1) ) @@ -262,6 +262,104 @@ SUBST [ X, type ] i2 ; integer*2 ;; END_TEMPLATE +IRP_IF INTEL + + subroutine sort(x,iorder,isize) + use intel + implicit none + BEGIN_DOC + ! Sort array x(isize). + ! iorder in input should be (1,2,3,...,isize), and in output + ! contains the new order of the elements. + END_DOC + integer,intent(in) :: isize + real ,intent(inout) :: x(isize) + integer,intent(inout) :: iorder(isize) + integer :: n + call ippsSortIndexAscend_32f_I(x, iorder, isize) + iorder(:) = iorder(:)+1 + end subroutine sort + + subroutine dsort(x,iorder,isize) + use intel + implicit none + BEGIN_DOC + ! Sort array x(isize). + ! iorder in input should be (1,2,3,...,isize), and in output + ! contains the new order of the elements. + END_DOC + integer,intent(in) :: isize + real(8) ,intent(inout) :: x(isize) + integer,intent(inout) :: iorder(isize) + integer :: n + call ippsSortIndexAscend_64f_I(x, iorder, isize) + iorder(:) = iorder(:)+1 + end subroutine dsort + + subroutine isort(x,iorder,isize) + use intel + implicit none + BEGIN_DOC + ! Sort array x(isize). + ! iorder in input should be (1,2,3,...,isize), and in output + ! contains the new order of the elements. + END_DOC + integer,intent(in) :: isize + integer ,intent(inout) :: x(isize) + integer,intent(inout) :: iorder(isize) + integer :: n + integer, allocatable :: iorder1(:) + allocate(iorder1(isize*2)) + n=4 + call ippsSortRadixIndexAscend_32s(x, n, iorder, isize, iorder1) + iorder(1:isize) = iorder(1:isize)+1 + deallocate(iorder1) + call iset_order(x,iorder,isize) + end subroutine isort + + subroutine isort_noidx(x,isize) + use intel + implicit none + BEGIN_DOC + ! Sort array x(isize). + ! iorder in input should be (1,2,3,...,isize), and in output + ! contains the new order of the elements. + END_DOC + integer,intent(in) :: isize + integer ,intent(inout) :: x(isize) + integer, allocatable :: iorder1(:) + integer :: n + call ippsSortRadixIndexGetBufferSize(isize, 11, n) + n = n/4 + allocate(iorder1(n)) + call ippsSortRadixAscend_32s_I(x, isize, iorder1) + deallocate(iorder1) + end subroutine isort_noidx + + +BEGIN_TEMPLATE + subroutine $Xsort(x,iorder,isize) + implicit none + BEGIN_DOC + ! Sort array x(isize). + ! iorder in input should be (1,2,3,...,isize), and in output + ! contains the new order of the elements. + END_DOC + integer,intent(in) :: isize + $type,intent(inout) :: x(isize) + integer,intent(inout) :: iorder(isize) + integer :: n +! call $Xradix_sort(x,iorder,isize,-1) + call quick_$Xsort(x,iorder,isize) + end subroutine $Xsort + +SUBST [ X, type ] + i8 ; integer*8 ;; + i2 ; integer*2 ;; +END_TEMPLATE + +IRP_ELSE + BEGIN_TEMPLATE subroutine $Xsort(x,iorder,isize) implicit none @@ -289,9 +387,9 @@ BEGIN_TEMPLATE endif end subroutine $Xsort -SUBST [ X, type, Y ] - ; real ; i ;; - d ; double precision ; i8 ;; +SUBST [ X, type ] + ; real ;; + d ; double precision ;; END_TEMPLATE BEGIN_TEMPLATE @@ -316,6 +414,22 @@ SUBST [ X, type ] i2 ; integer*2 ;; END_TEMPLATE + subroutine isort_noidx(x,isize) + implicit none + BEGIN_DOC + ! Sort array x(isize). + END_DOC + integer,intent(in) :: isize + $type,intent(inout) :: x(isize) + integer, allocatable :: iorder + allocate(iorder) + iorder=0 + call $Xradix_sort(x,iorder,isize,-1) + deallocate(iorder) + end subroutine $Xsort +IRP_ENDIF + + BEGIN_TEMPLATE subroutine $Xset_order(x,iorder,isize) implicit none @@ -413,10 +527,15 @@ SUBST [ X, type ] i2; integer*2 ;; END_TEMPLATE + BEGIN_TEMPLATE - recursive subroutine $Xradix_sort$big(x,iorder,isize,iradix) +recursive subroutine $Xradix_sort$big(x,iorder,isize,iradix) +IRP_IF INTEL + use intel +IRP_ENDIF implicit none + BEGIN_DOC ! Sort integer array x(isize) using the radix sort algorithm. ! iorder in input should be (1,2,3,...,isize), and in output @@ -448,6 +567,15 @@ BEGIN_TEMPLATE stop endif + IRP_IF INTEL + if ( ($type == 4).and.($integer_size == 32).and.($is_big == .False.) ) then + $intel + iorder(:) = iorder(:)+1 + return + endif + IRP_ENDIF + + i1=1_$int_type i2=1_$int_type do i=1_$int_type,isize @@ -637,12 +765,13 @@ BEGIN_TEMPLATE end -SUBST [ X, type, integer_size, is_big, big, int_type ] - i ; 4 ; 32 ; .False. ; ; 4 ;; - i8 ; 8 ; 64 ; .False. ; ; 4 ;; - i2 ; 2 ; 16 ; .False. ; ; 4 ;; - i ; 4 ; 32 ; .True. ; _big ; 8 ;; - i8 ; 8 ; 64 ; .True. ; _big ; 8 ;; +SUBST [ X, type, integer_size, is_big, big, int_type, intel ] + i ; 4 ; 32 ; .False. ; ; 4 ; call ippsSortRadixIndexAscend_32s(x, 4, iorder, isize, iorder1) ;; + i8 ; 8 ; 64 ; .False. ; ; 4 ; ;; + i2 ; 2 ; 16 ; .False. ; ; 4 ; ;; + i ; 4 ; 32 ; .True. ; _big ; 8 ; ;; + i8 ; 8 ; 64 ; .True. ; _big ; 8 ; ;; END_TEMPLATE + From c4a91c78c9011de18ae5d866b70693981b5485ed Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 31 May 2021 08:28:00 +0200 Subject: [PATCH 13/29] Fix gfortrab compilation --- src/utils/sort.irp.f | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/src/utils/sort.irp.f b/src/utils/sort.irp.f index be58e7a5..b4e30db4 100644 --- a/src/utils/sort.irp.f +++ b/src/utils/sort.irp.f @@ -420,13 +420,16 @@ END_TEMPLATE ! Sort array x(isize). END_DOC integer,intent(in) :: isize - $type,intent(inout) :: x(isize) - integer, allocatable :: iorder - allocate(iorder) - iorder=0 - call $Xradix_sort(x,iorder,isize,-1) + integer,intent(inout) :: x(isize) + integer, allocatable :: iorder(:) + integer :: i + allocate(iorder(isize)) + do i=1,isize + iorder(i)=i + enddo + call iradix_sort(x,iorder,isize,-1) deallocate(iorder) - end subroutine $Xsort + end subroutine isort_noidx IRP_ENDIF From 2c3a25e20a6e924de0ac63ba6a353c1db01455b9 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 31 May 2021 10:45:37 +0200 Subject: [PATCH 14/29] Cleaned sorting --- src/utils/intel.f90 | 125 ++++++++++++++++++++++++++--- src/utils/sort.irp.f | 187 ++++++++++++++----------------------------- 2 files changed, 173 insertions(+), 139 deletions(-) diff --git a/src/utils/intel.f90 b/src/utils/intel.f90 index 681adde9..4c18af8d 100644 --- a/src/utils/intel.f90 +++ b/src/utils/intel.f90 @@ -1,13 +1,5 @@ module intel use, intrinsic :: iso_c_binding - interface - subroutine ippsSortRadixIndexGetBufferSize(len, dataType, pBufSize) bind(C, name='ippsSortRadixIndexGetBufferSize') - use iso_c_binding - integer, intent(in), value :: len - integer, intent(in), value :: dataType - integer, intent(out) :: pBufSize - end - end interface interface subroutine ippsSortAscend_32s_I(pSrc, len) bind(C, name='ippsSortAscend_32s_I') use iso_c_binding @@ -15,12 +7,86 @@ module intel integer, intent(inout) :: pSrc(len) end end interface + interface + subroutine ippsSortAscend_32f_I(pSrc, len) bind(C, name='ippsSortAscend_32f_I') + use iso_c_binding + integer, intent(in), value :: len + real, intent(inout) :: pSrc(len) + end + end interface + interface + subroutine ippsSortAscend_64s_I(pSrc, len) bind(C, name='ippsSortAscend_64s_I') + use iso_c_binding + integer*8, intent(in), value :: len + integer, intent(inout) :: pSrc(len) + end + end interface + interface + subroutine ippsSortAscend_64f_I(pSrc, len) bind(C, name='ippsSortAscend_64f_I') + use iso_c_binding + double precision, intent(in), value :: len + real, intent(inout) :: pSrc(len) + end + end interface + + interface + subroutine ippsSortRadixIndexGetBufferSize(len, dataType, pBufSize) bind(C, name='ippsSortRadixIndexGetBufferSize') + use iso_c_binding + integer, intent(in), value :: len + integer, intent(in), value :: dataType + integer, intent(out) :: pBufSize + end + end interface + + interface + subroutine ippsSortRadixAscend_16s_I(pSrc, len, pTmp) bind(C, name='ippsSortRadixAscend_16s_I') + use iso_c_binding + integer, intent(in), value :: len + integer*2, intent(inout) :: pSrc(len) + character, intent(inout) :: pTmp(len) + end + end interface interface subroutine ippsSortRadixAscend_32s_I(pSrc, len, pTmp) bind(C, name='ippsSortRadixAscend_32s_I') use iso_c_binding integer, intent(in), value :: len integer, intent(inout) :: pSrc(len) - integer, intent(inout) :: pTmp(len) + character, intent(inout) :: pTmp(len) + end + end interface + interface + subroutine ippsSortRadixAscend_32f_I(pSrc, len, pTmp) bind(C, name='ippsSortRadixAscend_32f_I') + use iso_c_binding + integer, intent(in), value :: len + real, intent(inout) :: pSrc(len) + character, intent(inout) :: pTmp(len) + end + end interface + interface + subroutine ippsSortRadixAscend_64s_I(pSrc, len, pTmp) bind(C, name='ippsSortRadixAscend_64s_I') + use iso_c_binding + integer, intent(in), value :: len + integer*8, intent(inout) :: pSrc(len) + character, intent(inout) :: pTmp(len) + end + end interface + interface + subroutine ippsSortRadixAscend_64f_I(pSrc, len, pTmp) bind(C, name='ippsSortRadixAscend_64f_I') + use iso_c_binding + integer, intent(in), value :: len + double precision, intent(inout) :: pSrc(len) + character, intent(inout) :: pTmp(len) + end + end interface + + interface + subroutine ippsSortRadixIndexAscend_16s(pSrc, srcStrideBytes, pDstIndx, len, pTmpIndx) bind(C, name='ippsSortRadixIndexAscend_16s') + use iso_c_binding + integer, intent(in), value :: len + integer*2, intent(inout) :: pSrc(len) + integer, intent(in), value :: srcStrideBytes + integer, intent(inout) :: pDstIndx(len) + character, intent(inout) :: pTmpIndx(len) end end interface interface @@ -30,7 +96,7 @@ module intel integer, intent(inout) :: pSrc(len) integer, intent(in), value :: srcStrideBytes integer, intent(inout) :: pDstIndx(len) - integer, intent(inout) :: pTmpIndx(len) + character, intent(inout) :: pTmpIndx(len) end end interface interface @@ -40,9 +106,30 @@ module intel real , intent(inout) :: pSrc(len) integer, intent(in), value :: srcStrideBytes integer, intent(inout) :: pDstIndx(len) - integer, intent(inout) :: pTmpIndx(len) + character, intent(inout) :: pTmpIndx(len) end end interface + interface + subroutine ippsSortRadixIndexAscend_64s(pSrc, srcStrideBytes, pDstIndx, len, pTmpIndx) bind(C, name='ippsSortRadixIndexAscend_64s') + use iso_c_binding + integer, intent(in), value :: len + integer*8, intent(inout) :: pSrc(len) + integer, intent(in), value :: srcStrideBytes + integer, intent(inout) :: pDstIndx(len) + character, intent(inout) :: pTmpIndx(len) + end + end interface + interface + subroutine ippsSortRadixIndexAscend_64f(pSrc, srcStrideBytes, pDstIndx, len, pTmpIndx) bind(C,name='ippsSortRadixIndexAscend_64f') + use iso_c_binding + integer, intent(in), value :: len + real*8 , intent(inout) :: pSrc(len) + integer, intent(in), value :: srcStrideBytes + integer, intent(inout) :: pDstIndx(len) + character, intent(inout) :: pTmpIndx(len) + end + end interface + interface subroutine ippsSortIndexAscend_32f_I(pSrcDst, pDstIndx, len) bind(C,name='ippsSortIndexAscend_32f_I') use iso_c_binding @@ -67,4 +154,20 @@ module intel integer(4), intent(in), value :: len end end interface + interface + subroutine ippsSortIndexAscend_64s_I(pSrcDst, pDstIndx, len) bind(C,name='ippsSortIndexAscend_64s_I') + use iso_c_binding + integer(8), intent(in) :: pSrcDst(*) + integer(4), intent(inout) :: pDstIndx(*) + integer(4), intent(in), value :: len + end + end interface + interface + subroutine ippsSortIndexAscend_16s_I(pSrcDst, pDstIndx, len) bind(C,name='ippsSortIndexAscend_16s_I') + use iso_c_binding + integer(2), intent(in) :: pSrcDst(*) + integer(4), intent(inout) :: pDstIndx(*) + integer(4), intent(in), value :: len + end + end interface end module diff --git a/src/utils/sort.irp.f b/src/utils/sort.irp.f index b4e30db4..d40b7d1d 100644 --- a/src/utils/sort.irp.f +++ b/src/utils/sort.irp.f @@ -262,83 +262,13 @@ SUBST [ X, type ] i2 ; integer*2 ;; END_TEMPLATE + +!---------------------- INTEL IRP_IF INTEL - subroutine sort(x,iorder,isize) - use intel - implicit none - BEGIN_DOC - ! Sort array x(isize). - ! iorder in input should be (1,2,3,...,isize), and in output - ! contains the new order of the elements. - END_DOC - integer,intent(in) :: isize - real ,intent(inout) :: x(isize) - integer,intent(inout) :: iorder(isize) - integer :: n - call ippsSortIndexAscend_32f_I(x, iorder, isize) - iorder(:) = iorder(:)+1 - end subroutine sort - - subroutine dsort(x,iorder,isize) - use intel - implicit none - BEGIN_DOC - ! Sort array x(isize). - ! iorder in input should be (1,2,3,...,isize), and in output - ! contains the new order of the elements. - END_DOC - integer,intent(in) :: isize - real(8) ,intent(inout) :: x(isize) - integer,intent(inout) :: iorder(isize) - integer :: n - call ippsSortIndexAscend_64f_I(x, iorder, isize) - iorder(:) = iorder(:)+1 - end subroutine dsort - - subroutine isort(x,iorder,isize) - use intel - implicit none - BEGIN_DOC - ! Sort array x(isize). - ! iorder in input should be (1,2,3,...,isize), and in output - ! contains the new order of the elements. - END_DOC - integer,intent(in) :: isize - integer ,intent(inout) :: x(isize) - integer,intent(inout) :: iorder(isize) - integer :: n - integer, allocatable :: iorder1(:) - allocate(iorder1(isize*2)) - n=4 - call ippsSortRadixIndexAscend_32s(x, n, iorder, isize, iorder1) - iorder(1:isize) = iorder(1:isize)+1 - deallocate(iorder1) - call iset_order(x,iorder,isize) - end subroutine isort - - subroutine isort_noidx(x,isize) - use intel - implicit none - BEGIN_DOC - ! Sort array x(isize). - ! iorder in input should be (1,2,3,...,isize), and in output - ! contains the new order of the elements. - END_DOC - integer,intent(in) :: isize - integer ,intent(inout) :: x(isize) - integer, allocatable :: iorder1(:) - integer :: n - call ippsSortRadixIndexGetBufferSize(isize, 11, n) - n = n/4 - allocate(iorder1(n)) - call ippsSortRadixAscend_32s_I(x, isize, iorder1) - deallocate(iorder1) - end subroutine isort_noidx - - BEGIN_TEMPLATE subroutine $Xsort(x,iorder,isize) + use intel implicit none BEGIN_DOC ! Sort array x(isize). @@ -349,17 +279,46 @@ BEGIN_TEMPLATE $type,intent(inout) :: x(isize) integer,intent(inout) :: iorder(isize) integer :: n -! call $Xradix_sort(x,iorder,isize,-1) - call quick_$Xsort(x,iorder,isize) - end subroutine $Xsort + character, allocatable :: tmp(:) + if (isize < 2) return + call ippsSortRadixIndexGetBufferSize(isize, $ippsz, n) + allocate(tmp(n)) + call ippsSortRadixIndexAscend_$ityp(x, $n, iorder, isize, tmp) + deallocate(tmp) + iorder(1:isize) = iorder(1:isize)+1 + call $Xset_order(x,iorder,isize) + end -SUBST [ X, type ] - i8 ; integer*8 ;; - i2 ; integer*2 ;; + subroutine $Xsort_noidx(x,isize) + use intel + implicit none + BEGIN_DOC + ! Sort array x(isize). + ! iorder in input should be (1,2,3,...,isize), and in output + ! contains the new order of the elements. + END_DOC + integer,intent(in) :: isize + $type,intent(inout) :: x(isize) + integer :: n + character, allocatable :: tmp(:) + if (isize < 2) return + call ippsSortRadixIndexGetBufferSize(isize, $ippsz, n) + allocate(tmp(n)) + call ippsSortRadixAscend_$ityp_I(x, isize, tmp) + deallocate(tmp) + end + +SUBST [ X, type, ityp, n, ippsz ] + ; real ; 32f ; 4 ; 13 ;; + d ; double precision ; 64f ; 8 ; 19;; + i ; integer ; 32s ; 4 ; 11 ;; + i8 ; integer*8 ; 64s ; 8 ; 17;; + i2 ; integer*2 ; 16s ; 2 ; 7 ;; END_TEMPLATE +!---------------------- END INTEL IRP_ELSE - +!---------------------- NON-INTEL BEGIN_TEMPLATE subroutine $Xsort(x,iorder,isize) implicit none @@ -387,50 +346,34 @@ BEGIN_TEMPLATE endif end subroutine $Xsort -SUBST [ X, type ] - ; real ;; - d ; double precision ;; -END_TEMPLATE - -BEGIN_TEMPLATE - subroutine $Xsort(x,iorder,isize) + subroutine $Xsort_noidx(x,isize) implicit none BEGIN_DOC ! Sort array x(isize). - ! iorder in input should be (1,2,3,...,isize), and in output - ! contains the new order of the elements. END_DOC integer,intent(in) :: isize $type,intent(inout) :: x(isize) - integer,intent(inout) :: iorder(isize) - integer :: n -! call $Xradix_sort(x,iorder,isize,-1) - call quick_$Xsort(x,iorder,isize) - end subroutine $Xsort + integer, allocatable :: iorder(:) + integer :: i + allocate(iorder(isize)) + do i=1,isize + iorder(i)=i + enddo + call $Xsort(x,iorder,isize) + deallocate(iorder) + end subroutine $Xsort_noidx SUBST [ X, type ] + ; real ;; + d ; double precision ;; i ; integer ;; i8 ; integer*8 ;; i2 ; integer*2 ;; END_TEMPLATE - subroutine isort_noidx(x,isize) - implicit none - BEGIN_DOC - ! Sort array x(isize). - END_DOC - integer,intent(in) :: isize - integer,intent(inout) :: x(isize) - integer, allocatable :: iorder(:) - integer :: i - allocate(iorder(isize)) - do i=1,isize - iorder(i)=i - enddo - call iradix_sort(x,iorder,isize,-1) - deallocate(iorder) - end subroutine isort_noidx IRP_ENDIF +!---------------------- END NON-INTEL + BEGIN_TEMPLATE @@ -534,9 +477,6 @@ END_TEMPLATE BEGIN_TEMPLATE recursive subroutine $Xradix_sort$big(x,iorder,isize,iradix) -IRP_IF INTEL - use intel -IRP_ENDIF implicit none BEGIN_DOC @@ -570,15 +510,6 @@ IRP_ENDIF stop endif - IRP_IF INTEL - if ( ($type == 4).and.($integer_size == 32).and.($is_big == .False.) ) then - $intel - iorder(:) = iorder(:)+1 - return - endif - IRP_ENDIF - - i1=1_$int_type i2=1_$int_type do i=1_$int_type,isize @@ -768,12 +699,12 @@ IRP_ENDIF end -SUBST [ X, type, integer_size, is_big, big, int_type, intel ] - i ; 4 ; 32 ; .False. ; ; 4 ; call ippsSortRadixIndexAscend_32s(x, 4, iorder, isize, iorder1) ;; - i8 ; 8 ; 64 ; .False. ; ; 4 ; ;; - i2 ; 2 ; 16 ; .False. ; ; 4 ; ;; - i ; 4 ; 32 ; .True. ; _big ; 8 ; ;; - i8 ; 8 ; 64 ; .True. ; _big ; 8 ; ;; +SUBST [ X, type, integer_size, is_big, big, int_type ] + i ; 4 ; 32 ; .False. ; ; 4 ;; + i8 ; 8 ; 64 ; .False. ; ; 4 ;; + i2 ; 2 ; 16 ; .False. ; ; 4 ;; + i ; 4 ; 32 ; .True. ; _big ; 8 ;; + i8 ; 8 ; 64 ; .True. ; _big ; 8 ;; END_TEMPLATE From 3837dea58c0a5cc30996b6fec3ecc4824d774db7 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 31 May 2021 11:47:34 +0200 Subject: [PATCH 15/29] Improving sort --- src/utils/sort.irp.f | 75 ++++++++++++++++++++++++++++++++++---------- 1 file changed, 58 insertions(+), 17 deletions(-) diff --git a/src/utils/sort.irp.f b/src/utils/sort.irp.f index d40b7d1d..f3879d4a 100644 --- a/src/utils/sort.irp.f +++ b/src/utils/sort.irp.f @@ -320,6 +320,34 @@ END_TEMPLATE IRP_ELSE !---------------------- NON-INTEL BEGIN_TEMPLATE + + subroutine $Xsort_noidx(x,isize) + implicit none + BEGIN_DOC + ! Sort array x(isize). + END_DOC + integer,intent(in) :: isize + $type,intent(inout) :: x(isize) + integer, allocatable :: iorder(:) + integer :: i + allocate(iorder(isize)) + do i=1,isize + iorder(i)=i + enddo + call $Xsort(x,iorder,isize) + deallocate(iorder) + end subroutine $Xsort_noidx + +SUBST [ X, type ] + ; real ;; + d ; double precision ;; + i ; integer ;; + i8 ; integer*8 ;; + i2 ; integer*2 ;; +END_TEMPLATE + +BEGIN_TEMPLATE + subroutine $Xsort(x,iorder,isize) implicit none BEGIN_DOC @@ -346,26 +374,39 @@ BEGIN_TEMPLATE endif end subroutine $Xsort - subroutine $Xsort_noidx(x,isize) - implicit none - BEGIN_DOC - ! Sort array x(isize). - END_DOC - integer,intent(in) :: isize - $type,intent(inout) :: x(isize) - integer, allocatable :: iorder(:) - integer :: i - allocate(iorder(isize)) - do i=1,isize - iorder(i)=i - enddo - call $Xsort(x,iorder,isize) - deallocate(iorder) - end subroutine $Xsort_noidx - SUBST [ X, type ] ; real ;; d ; double precision ;; +END_TEMPLATE + +BEGIN_TEMPLATE + + subroutine $Xsort(x,iorder,isize) + implicit none + BEGIN_DOC + ! Sort array x(isize). + ! iorder in input should be (1,2,3,...,isize), and in output + ! contains the new order of the elements. + END_DOC + integer,intent(in) :: isize + $type,intent(inout) :: x(isize) + integer,intent(inout) :: iorder(isize) + integer :: n + if (isize < 2) then + return + endif + call sorted_$Xnumber(x,isize,n) + if (isize == n) then + return + endif + if ( isize < 32) then + call insertion_$Xsort(x,iorder,isize) + else + call $Xradix_sort(x,iorder,isize,-1) + endif + end subroutine $Xsort + +SUBST [ X, type ] i ; integer ;; i8 ; integer*8 ;; i2 ; integer*2 ;; From d45f6091daf107a2a18925c9b75d9a8ed4963eee Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 31 May 2021 12:16:00 +0200 Subject: [PATCH 16/29] Sorting compatible with old IPP --- config/ifort_rome.cfg | 4 +-- src/utils/intel.f90 | 8 +++--- src/utils/sort.irp.f | 65 +++++++++++++++++++++++++++++++++++++++++-- 3 files changed, 69 insertions(+), 8 deletions(-) diff --git a/config/ifort_rome.cfg b/config/ifort_rome.cfg index 3b70f98b..5ed01227 100644 --- a/config/ifort_rome.cfg +++ b/config/ifort_rome.cfg @@ -31,8 +31,8 @@ OPENMP : 1 ; Append OpenMP flags # -ftz : Flushes denormal results to zero # [OPT] -FC : -traceback -FCFLAGS : -O2 -ip -g -march=core-avx2 -align array64byte -fma -ftz -fomit-frame-pointer +FC : -traceback -shared-intel +FCFLAGS : -O2 -ip -g -march=core-avx2 -align array64byte -fma -ftz -fomit-frame-pointer # Profiling flags ################# diff --git a/src/utils/intel.f90 b/src/utils/intel.f90 index 4c18af8d..2b61cb38 100644 --- a/src/utils/intel.f90 +++ b/src/utils/intel.f90 @@ -17,15 +17,15 @@ module intel interface subroutine ippsSortAscend_64s_I(pSrc, len) bind(C, name='ippsSortAscend_64s_I') use iso_c_binding - integer*8, intent(in), value :: len - integer, intent(inout) :: pSrc(len) + integer, intent(in), value :: len + integer*8, intent(inout) :: pSrc(len) end end interface interface subroutine ippsSortAscend_64f_I(pSrc, len) bind(C, name='ippsSortAscend_64f_I') use iso_c_binding - double precision, intent(in), value :: len - real, intent(inout) :: pSrc(len) + integer, intent(in), value :: len + double precision, intent(inout) :: pSrc(len) end end interface diff --git a/src/utils/sort.irp.f b/src/utils/sort.irp.f index f3879d4a..21eb8b67 100644 --- a/src/utils/sort.irp.f +++ b/src/utils/sort.irp.f @@ -310,12 +310,73 @@ BEGIN_TEMPLATE SUBST [ X, type, ityp, n, ippsz ] ; real ; 32f ; 4 ; 13 ;; - d ; double precision ; 64f ; 8 ; 19;; i ; integer ; 32s ; 4 ; 11 ;; - i8 ; integer*8 ; 64s ; 8 ; 17;; i2 ; integer*2 ; 16s ; 2 ; 7 ;; END_TEMPLATE +BEGIN_TEMPLATE + + subroutine $Xsort(x,iorder,isize) + implicit none + BEGIN_DOC + ! Sort array x(isize). + ! iorder in input should be (1,2,3,...,isize), and in output + ! contains the new order of the elements. + END_DOC + integer,intent(in) :: isize + $type,intent(inout) :: x(isize) + integer,intent(inout) :: iorder(isize) + integer :: n + if (isize < 2) then + return + endif +! call sorted_$Xnumber(x,isize,n) +! if (isize == n) then +! return +! endif + if ( isize < 32) then + call insertion_$Xsort(x,iorder,isize) + else +! call heap_$Xsort(x,iorder,isize) + call quick_$Xsort(x,iorder,isize) + endif + end subroutine $Xsort + +SUBST [ X, type ] + d ; double precision ;; +END_TEMPLATE + +BEGIN_TEMPLATE + + subroutine $Xsort(x,iorder,isize) + implicit none + BEGIN_DOC + ! Sort array x(isize). + ! iorder in input should be (1,2,3,...,isize), and in output + ! contains the new order of the elements. + END_DOC + integer,intent(in) :: isize + $type,intent(inout) :: x(isize) + integer,intent(inout) :: iorder(isize) + integer :: n + if (isize < 2) then + return + endif + call sorted_$Xnumber(x,isize,n) + if (isize == n) then + return + endif + if ( isize < 32) then + call insertion_$Xsort(x,iorder,isize) + else + call $Xradix_sort(x,iorder,isize,-1) + endif + end subroutine $Xsort + +SUBST [ X, type ] + i8 ; integer*8 ;; +END_TEMPLATE + !---------------------- END INTEL IRP_ELSE !---------------------- NON-INTEL From 00dfa11f4f013f91bf49d99d7300bbe98cbe3da5 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 31 May 2021 19:02:25 +0200 Subject: [PATCH 17/29] Removed OMP in sorting --- src/utils/sort.irp.f | 30 +----------------------------- 1 file changed, 1 insertion(+), 29 deletions(-) diff --git a/src/utils/sort.irp.f b/src/utils/sort.irp.f index 21eb8b67..a63eb4a3 100644 --- a/src/utils/sort.irp.f +++ b/src/utils/sort.irp.f @@ -38,15 +38,7 @@ BEGIN_TEMPLATE $type,intent(inout) :: x(isize) integer,intent(inout) :: iorder(isize) integer, external :: omp_get_num_threads - if (omp_get_num_threads() == 1) then - !$OMP PARALLEL DEFAULT(SHARED) - !$OMP SINGLE - call rec_$X_quicksort(x,iorder,isize,1,isize,nproc) - !$OMP END SINGLE - !$OMP END PARALLEL - else - call rec_$X_quicksort(x,iorder,isize,1,isize,nproc) - endif + call rec_$X_quicksort(x,iorder,isize,1,isize,nproc) end recursive subroutine rec_$X_quicksort(x, iorder, isize, first, last, level) @@ -89,16 +81,11 @@ BEGIN_TEMPLATE endif else if (first < i-1) then - !$OMP TASK DEFAULT(SHARED) FIRSTPRIVATE(isize,first,i,level) call rec_$X_quicksort(x, iorder, isize, first, i-1,level/2) - !$OMP END TASK endif if (j+1 < last) then - !$OMP TASK DEFAULT(SHARED) FIRSTPRIVATE(isize,last,j,level) call rec_$X_quicksort(x, iorder, isize, j+1, last,level/2) - !$OMP END TASK endif - !$OMP TASKWAIT endif end @@ -716,24 +703,14 @@ recursive subroutine $Xradix_sort$big(x,iorder,isize,iradix) endif -! !$OMP PARALLEL DEFAULT(SHARED) if (isize > 1000000) -! !$OMP SINGLE if (i3>1_$int_type) then -! !$OMP TASK FIRSTPRIVATE(iradix_new,i3) SHARED(x,iorder) if(i3 > 1000000) call $Xradix_sort$big(x,iorder,i3,iradix_new-1) -! !$OMP END TASK endif if (isize-i3>1_$int_type) then -! !$OMP TASK FIRSTPRIVATE(iradix_new,i3) SHARED(x,iorder) if(isize-i3 > 1000000) call $Xradix_sort$big(x(i3+1_$int_type),iorder(i3+1_$int_type),isize-i3,iradix_new-1) -! !$OMP END TASK endif -! !$OMP TASKWAIT -! !$OMP END SINGLE -! !$OMP END PARALLEL - return endif @@ -788,16 +765,11 @@ recursive subroutine $Xradix_sort$big(x,iorder,isize,iradix) if (i1>1_$int_type) then - !$OMP TASK FIRSTPRIVATE(i0,iradix,i1) SHARED(x,iorder) if(i1 >1000000) call $Xradix_sort$big(x(i0+1_$int_type),iorder(i0+1_$int_type),i1,iradix-1) - !$OMP END TASK endif if (i0>1) then - !$OMP TASK FIRSTPRIVATE(i0,iradix) SHARED(x,iorder) if(i0 >1000000) call $Xradix_sort$big(x,iorder,i0,iradix-1) - !$OMP END TASK endif - !$OMP TASKWAIT end From d174a1f7cc0a4dda35d7556f9a3f7d5b2b29711a Mon Sep 17 00:00:00 2001 From: vijay gopal chilkuri Date: Mon, 31 May 2021 23:45:05 +0530 Subject: [PATCH 18/29] Fixed and verfied bug in n_CSF for Benzene singlet. #158 --- src/csf/sigma_vector.irp.f | 67 +++++++++++++++++++++++++++++++------- 1 file changed, 56 insertions(+), 11 deletions(-) diff --git a/src/csf/sigma_vector.irp.f b/src/csf/sigma_vector.irp.f index 581a6a6b..75d9f1f7 100644 --- a/src/csf/sigma_vector.irp.f +++ b/src/csf/sigma_vector.irp.f @@ -1,9 +1,15 @@ - BEGIN_PROVIDER [ integer, NSOMOMax] -&BEGIN_PROVIDER [ integer, NCSFMax] -&BEGIN_PROVIDER [ integer*8, NMO] -&BEGIN_PROVIDER [ integer, NBFMax] -&BEGIN_PROVIDER [ integer, n_CSF] -&BEGIN_PROVIDER [ integer, maxDetDimPerBF] + real*8 function lgamma(x) + implicit none + real*8, intent(in) :: x + lgamma = log(abs(gamma(x))) + end function lgamma + + BEGIN_PROVIDER [ integer, NSOMOMax] + &BEGIN_PROVIDER [ integer, NCSFMax] + &BEGIN_PROVIDER [ integer*8, NMO] + &BEGIN_PROVIDER [ integer, NBFMax] + &BEGIN_PROVIDER [ integer, n_CSF] + &BEGIN_PROVIDER [ integer, maxDetDimPerBF] implicit none BEGIN_DOC ! Documentation for NSOMOMax @@ -22,28 +28,67 @@ integer NSOMO integer dimcsfpercfg integer detDimperBF - real*8 :: coeff + real*8 :: coeff, binom1, binom2 integer MS integer ncfgpersomo + real*8, external :: lgamma detDimperBF = 0 MS = elec_alpha_num-elec_beta_num ! number of cfgs = number of dets for 0 somos n_CSF = 0 ncfgprev = cfg_seniority_index(0) + ncfgpersomo = ncfgprev + do i = 1, elec_num + print *,"i=",i," Ncfg= ",cfg_seniority_index(i) + enddo do i = iand(MS,1), NSOMOMax-2,2 + if(cfg_seniority_index(i) .EQ. -1) then + cycle + endif if(cfg_seniority_index(i+2) .EQ. -1) then ncfgpersomo = N_configuration + 1 else - ncfgpersomo = cfg_seniority_index(i+2) + if(cfg_seniority_index(i+2) > ncfgpersomo) then + ncfgpersomo = cfg_seniority_index(i+2) + else + k = 0 + do while(cfg_seniority_index(i+2+k) < ncfgpersomo) + k = k + 2 + ncfgpersomo = cfg_seniority_index(i+2+k) + enddo + endif endif ncfg = ncfgpersomo - ncfgprev if(iand(MS,1) .EQ. 0) then - dimcsfpercfg = max(1,nint((binom(i,i/2)-binom(i,i/2+1)))) + !dimcsfpercfg = max(1,nint((binom(i,i/2)-binom(i,i/2+1)))) + binom1 = dexp(lgamma(1.0d0*(i+1)) & + - lgamma(1.0d0*((i/2)+1)) & + - lgamma(1.0d0*(i-((i/2))+1))); + binom2 = dexp(lgamma(1.0d0*(i+1)) & + - lgamma(1.0d0*(((i/2)+1)+1)) & + - lgamma(1.0d0*(i-((i/2)+1)+1))); + dimcsfpercfg = max(1,nint(binom1 - binom2)) else - dimcsfpercfg = max(1,nint((binom(i,(i+1)/2)-binom(i,(i+3)/2)))) + !dimcsfpercfg = max(1,nint((binom(i,(i+1)/2)-binom(i,(i+3)/2)))) + binom1 = dexp(lgamma(1.0d0*(i+1)) & + - lgamma(1.0d0*(((i+1)/2)+1)) & + - lgamma(1.0d0*(i-(((i+1)/2))+1))); + binom2 = dexp(lgamma(1.0d0*(i+1)) & + - lgamma(1.0d0*((((i+3)/2)+1)+1)) & + - lgamma(1.0d0*(i-(((i+3)/2)+1)+1))); + dimcsfpercfg = max(1,nint(binom1 - binom2)) endif n_CSF += ncfg * dimcsfpercfg - ncfgprev = cfg_seniority_index(i+2) + print *,"i=",i," ncfg= ", ncfg, " dims=", dimcsfpercfg, " n_csf=", n_CSF, ncfgpersomo, ncfgprev + if(cfg_seniority_index(i+2) > ncfgprev) then + ncfgprev = cfg_seniority_index(i+2) + else + k = 0 + do while(cfg_seniority_index(i+2+k) < ncfgprev) + k = k + 2 + ncfgprev = cfg_seniority_index(i+2+k) + enddo + endif enddo END_PROVIDER From d584862f09462058457f01529f20e2906e650752 Mon Sep 17 00:00:00 2001 From: vijay gopal chilkuri Date: Tue, 1 Jun 2021 10:05:18 +0530 Subject: [PATCH 19/29] Remove debug print. --- src/csf/sigma_vector.irp.f | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/csf/sigma_vector.irp.f b/src/csf/sigma_vector.irp.f index 75d9f1f7..c1b0b228 100644 --- a/src/csf/sigma_vector.irp.f +++ b/src/csf/sigma_vector.irp.f @@ -38,9 +38,6 @@ n_CSF = 0 ncfgprev = cfg_seniority_index(0) ncfgpersomo = ncfgprev - do i = 1, elec_num - print *,"i=",i," Ncfg= ",cfg_seniority_index(i) - enddo do i = iand(MS,1), NSOMOMax-2,2 if(cfg_seniority_index(i) .EQ. -1) then cycle @@ -79,7 +76,6 @@ dimcsfpercfg = max(1,nint(binom1 - binom2)) endif n_CSF += ncfg * dimcsfpercfg - print *,"i=",i," ncfg= ", ncfg, " dims=", dimcsfpercfg, " n_csf=", n_CSF, ncfgpersomo, ncfgprev if(cfg_seniority_index(i+2) > ncfgprev) then ncfgprev = cfg_seniority_index(i+2) else From e751901fa82a65ba88da87962782c1540b6d7fe2 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 1 Jun 2021 10:35:33 +0200 Subject: [PATCH 20/29] Renamed lgamma into logabsgamma (lgamma is a Fortran intrinsic) --- src/csf/sigma_vector.irp.f | 32 +++++++++++++++---------------- src/davidson/diagonalize_ci.irp.f | 4 ++-- 2 files changed, 18 insertions(+), 18 deletions(-) diff --git a/src/csf/sigma_vector.irp.f b/src/csf/sigma_vector.irp.f index c1b0b228..85ed5f84 100644 --- a/src/csf/sigma_vector.irp.f +++ b/src/csf/sigma_vector.irp.f @@ -1,8 +1,8 @@ - real*8 function lgamma(x) + real*8 function logabsgamma(x) implicit none real*8, intent(in) :: x - lgamma = log(abs(gamma(x))) - end function lgamma + logabsgamma = log(abs(gamma(x))) + end function logabsgamma BEGIN_PROVIDER [ integer, NSOMOMax] &BEGIN_PROVIDER [ integer, NCSFMax] @@ -31,7 +31,7 @@ real*8 :: coeff, binom1, binom2 integer MS integer ncfgpersomo - real*8, external :: lgamma + real*8, external :: logabsgamma detDimperBF = 0 MS = elec_alpha_num-elec_beta_num ! number of cfgs = number of dets for 0 somos @@ -58,21 +58,21 @@ ncfg = ncfgpersomo - ncfgprev if(iand(MS,1) .EQ. 0) then !dimcsfpercfg = max(1,nint((binom(i,i/2)-binom(i,i/2+1)))) - binom1 = dexp(lgamma(1.0d0*(i+1)) & - - lgamma(1.0d0*((i/2)+1)) & - - lgamma(1.0d0*(i-((i/2))+1))); - binom2 = dexp(lgamma(1.0d0*(i+1)) & - - lgamma(1.0d0*(((i/2)+1)+1)) & - - lgamma(1.0d0*(i-((i/2)+1)+1))); + binom1 = dexp(logabsgamma(1.0d0*(i+1)) & + - logabsgamma(1.0d0*((i/2)+1)) & + - logabsgamma(1.0d0*(i-((i/2))+1))); + binom2 = dexp(logabsgamma(1.0d0*(i+1)) & + - logabsgamma(1.0d0*(((i/2)+1)+1)) & + - logabsgamma(1.0d0*(i-((i/2)+1)+1))); dimcsfpercfg = max(1,nint(binom1 - binom2)) else !dimcsfpercfg = max(1,nint((binom(i,(i+1)/2)-binom(i,(i+3)/2)))) - binom1 = dexp(lgamma(1.0d0*(i+1)) & - - lgamma(1.0d0*(((i+1)/2)+1)) & - - lgamma(1.0d0*(i-(((i+1)/2))+1))); - binom2 = dexp(lgamma(1.0d0*(i+1)) & - - lgamma(1.0d0*((((i+3)/2)+1)+1)) & - - lgamma(1.0d0*(i-(((i+3)/2)+1)+1))); + binom1 = dexp(logabsgamma(1.0d0*(i+1)) & + - logabsgamma(1.0d0*(((i+1)/2)+1)) & + - logabsgamma(1.0d0*(i-(((i+1)/2))+1))); + binom2 = dexp(logabsgamma(1.0d0*(i+1)) & + - logabsgamma(1.0d0*((((i+3)/2)+1)+1)) & + - logabsgamma(1.0d0*(i-(((i+3)/2)+1)+1))); dimcsfpercfg = max(1,nint(binom1 - binom2)) endif n_CSF += ncfg * dimcsfpercfg diff --git a/src/davidson/diagonalize_ci.irp.f b/src/davidson/diagonalize_ci.irp.f index da5fb950..06744934 100644 --- a/src/davidson/diagonalize_ci.irp.f +++ b/src/davidson/diagonalize_ci.irp.f @@ -57,8 +57,8 @@ END_PROVIDER enddo ! Deactivated temporarily: bug in N_csf -! do_csf = s2_eig .and. only_expected_s2 .and. (expected_s2 == 0.d0) - do_csf = .False. + do_csf = s2_eig .and. only_expected_s2 .and. (expected_s2 == 0.d0) +! do_csf = .False. if (diag_algorithm == "Davidson") then From 65e124380af281cbb6b1ad54837ec8f74a12e7e8 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 1 Jun 2021 11:33:39 +0200 Subject: [PATCH 21/29] Fix bug with non-continuous active MOs and deleted --- src/bitmask/bitmasks.irp.f | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/bitmask/bitmasks.irp.f b/src/bitmask/bitmasks.irp.f index 91617397..27c7dc84 100644 --- a/src/bitmask/bitmasks.irp.f +++ b/src/bitmask/bitmasks.irp.f @@ -27,9 +27,7 @@ BEGIN_PROVIDER [ integer(bit_kind), full_ijkl_bitmask, (N_int) ] full_ijkl_bitmask(j) = 0_bit_kind do i=0,bit_kind_size-1 k=k+1 - if (mo_class(k) /= 'Deleted') then - full_ijkl_bitmask(j) = ibset(full_ijkl_bitmask(j),i) - endif + full_ijkl_bitmask(j) = ibset(full_ijkl_bitmask(j),i) if (k == mo_num) exit enddo enddo From 49a43be70b096a70ee0248d45c9b9165aab889ae Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 1 Jun 2021 11:37:34 +0200 Subject: [PATCH 22/29] Release notes --- RELEASE_NOTES.org | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/RELEASE_NOTES.org b/RELEASE_NOTES.org index 02176ffa..758618a9 100644 --- a/RELEASE_NOTES.org +++ b/RELEASE_NOTES.org @@ -29,6 +29,7 @@ - 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 *** User interface @@ -54,6 +55,8 @@ - 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 *** Code @@ -78,6 +81,8 @@ - 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 ao_one_e_integral_zero banned_excitations From d31d8435357d8ec4b415c3340929cdd08023ba8b Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 1 Jun 2021 14:51:50 +0200 Subject: [PATCH 23/29] Update EZFIO --- external/ezfio | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/external/ezfio b/external/ezfio index ccee52d0..ed1df9f3 160000 --- a/external/ezfio +++ b/external/ezfio @@ -1 +1 @@ -Subproject commit ccee52d00c2cde1d628b0d34f4a247143747bf36 +Subproject commit ed1df9f3c1f51752656ca98da5693a4119add05c From e7ed68205878298024d4bf1fcc8336b71a9a2ced Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 1 Jun 2021 17:09:59 +0200 Subject: [PATCH 24/29] Added normalization factors for basis --- external/ezfio | 2 +- src/basis/EZFIO.cfg | 22 ++++++++----- src/basis/basis.irp.f | 75 +++++++++++++++++++++++++++++++++++++------ 3 files changed, 81 insertions(+), 18 deletions(-) diff --git a/external/ezfio b/external/ezfio index ed1df9f3..ccee52d0 160000 --- a/external/ezfio +++ b/external/ezfio @@ -1 +1 @@ -Subproject commit ed1df9f3c1f51752656ca98da5693a4119add05c +Subproject commit ccee52d00c2cde1d628b0d34f4a247143747bf36 diff --git a/src/basis/EZFIO.cfg b/src/basis/EZFIO.cfg index 1c66e758..42123373 100644 --- a/src/basis/EZFIO.cfg +++ b/src/basis/EZFIO.cfg @@ -15,15 +15,10 @@ interface: ezfio, provider [shell_normalization_factor] type: double precision -doc: Number of primitives per |AO| +doc: Normalization factor applied to the whole shell, ex $1/\sqrt{ }$ size: (basis.shell_num) interface: ezfio -[prim_num] -type: integer -doc: Total number of primitives -interface: ezfio, provider - [shell_ang_mom] type: integer doc: Angular momentum of each shell @@ -48,13 +43,24 @@ doc: Index of the nucleus on which the shell is centered size: (basis.shell_num) interface: ezfio, provider -[shell_prim_coef] +[prim_normalization_factor] +type: double precision +doc: Normalization factor applied to each primitive +size: (basis.prim_num) +interface: ezfio + +[prim_num] +type: integer +doc: Total number of primitives +interface: ezfio, provider + +[prim_coef] type: double precision doc: Primitive coefficients size: (basis.prim_num) interface: ezfio, provider -[shell_prim_expo] +[prim_expo] type: double precision doc: Exponents in the shell size: (basis.prim_num) diff --git a/src/basis/basis.irp.f b/src/basis/basis.irp.f index e25e6a46..a247a90e 100644 --- a/src/basis/basis.irp.f +++ b/src/basis/basis.irp.f @@ -15,7 +15,6 @@ BEGIN_PROVIDER [ double precision, shell_normalization_factor , (shell_num) ] call ezfio_get_basis_shell_normalization_factor(shell_normalization_factor) else - double precision :: norm,overlap_x,overlap_y,overlap_z,C_A(3), c integer :: l, powA(3), nz integer :: i,j,k @@ -25,23 +24,22 @@ BEGIN_PROVIDER [ double precision, shell_normalization_factor , (shell_num) ] C_A(3) = 0.d0 do i=1,shell_num + powA(1) = shell_ang_mom(i) powA(2) = 0 powA(3) = 0 - ! Normalization of the contracted basis functions norm = 0.d0 - do j=shell_prim_index(i), shell_prim_index(i)+shell_prim_num(i)-1 - do k=shell_prim_index(i), shell_prim_index(i)+shell_prim_num(i)-1 - call overlap_gaussian_xyz(C_A,C_A,shell_prim_expo(j),shell_prim_expo(k),& - powA,powA,overlap_x,overlap_y,overlap_z,c,nz) - norm = norm+c*shell_prim_coef(j)*shell_prim_coef(k) + do k=shell_prim_index(i),shell_prim_index(i)+shell_prim_num(i)-1 + do j=shell_prim_index(i),shell_prim_index(i)+shell_prim_num(i)-1 + call overlap_gaussian_xyz(C_A,C_A,prim_expo(j),prim_expo(k), & + powA,powA,overlap_x,overlap_y,overlap_z,c,nz) + norm = norm+c*prim_coef(j)*prim_coef(k) enddo enddo - shell_normalization_factor(i) = dsqrt(norm) + shell_normalization_factor(i) = 1.d0/dsqrt(norm) enddo - endif endif IRP_IF MPI_DEBUG @@ -60,3 +58,62 @@ BEGIN_PROVIDER [ double precision, shell_normalization_factor , (shell_num) ] call write_time(6) END_PROVIDER + + +BEGIN_PROVIDER [ double precision, prim_normalization_factor , (prim_num) ] + implicit none + BEGIN_DOC + ! Number of primitives per |AO| + END_DOC + + logical :: has + PROVIDE ezfio_filename + if (mpi_master) then + if (size(prim_normalization_factor) == 0) return + + call ezfio_has_basis_prim_normalization_factor(has) + if (has) then + write(6,'(A)') '.. >>>>> [ IO READ: prim_normalization_factor ] <<<<< ..' + call ezfio_get_basis_prim_normalization_factor(prim_normalization_factor) + else + + double precision :: norm,overlap_x,overlap_y,overlap_z,C_A(3), c + integer :: l, powA(3), nz + integer :: i,j,k + nz=100 + C_A(1) = 0.d0 + C_A(2) = 0.d0 + C_A(3) = 0.d0 + + do i=1,shell_num + + powA(1) = shell_ang_mom(i) + powA(2) = 0 + powA(3) = 0 + + do k=shell_prim_index(i),shell_prim_index(i)+shell_prim_num(i)-1 + call overlap_gaussian_xyz(C_A,C_A,prim_expo(k),prim_expo(k), & + powA,powA,overlap_x,overlap_y,overlap_z,norm,nz) + prim_normalization_factor(k) = 1.d0/dsqrt(norm) + enddo + enddo + + + endif + 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( prim_normalization_factor, (prim_num), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + if (ierr /= MPI_SUCCESS) then + stop 'Unable to read prim_normalization_factor with MPI' + endif + IRP_ENDIF + + call write_time(6) + +END_PROVIDER From 63d407e4d26ebaece5cb3e1ca75846c448f07586 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 1 Jun 2021 19:46:15 +0200 Subject: [PATCH 25/29] Normalization in basis --- ocaml/qp_create_ezfio.ml | 4 ++-- src/basis/basis.irp.f | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/ocaml/qp_create_ezfio.ml b/ocaml/qp_create_ezfio.ml index bf3018b3..e4ead21b 100644 --- a/ocaml/qp_create_ezfio.ml +++ b/ocaml/qp_create_ezfio.ml @@ -607,9 +607,9 @@ let run ?o b au c d m p cart xyz_file = Ezfio.set_basis_shell_nucl (Ezfio.ezfio_array_of_list ~rank:1 ~dim:[| shell_num |] ~data:(list_map (fun (_,n) -> Nucl_number.to_int n) basis) ) ; - Ezfio.set_basis_shell_prim_coef (Ezfio.ezfio_array_of_list + Ezfio.set_basis_prim_coef (Ezfio.ezfio_array_of_list ~rank:1 ~dim:[| prim_num |] ~data:coef) ; - Ezfio.set_basis_shell_prim_expo (Ezfio.ezfio_array_of_list + Ezfio.set_basis_prim_expo (Ezfio.ezfio_array_of_list ~rank:1 ~dim:[| prim_num |] ~data:expo) ; diff --git a/src/basis/basis.irp.f b/src/basis/basis.irp.f index a247a90e..6a406e28 100644 --- a/src/basis/basis.irp.f +++ b/src/basis/basis.irp.f @@ -34,7 +34,7 @@ BEGIN_PROVIDER [ double precision, shell_normalization_factor , (shell_num) ] do j=shell_prim_index(i),shell_prim_index(i)+shell_prim_num(i)-1 call overlap_gaussian_xyz(C_A,C_A,prim_expo(j),prim_expo(k), & powA,powA,overlap_x,overlap_y,overlap_z,c,nz) - norm = norm+c*prim_coef(j)*prim_coef(k) + norm = norm+c*prim_coef(j)*prim_coef(k) * prim_normalization_factor(j) * prim_normalization_factor(k) enddo enddo shell_normalization_factor(i) = 1.d0/dsqrt(norm) From 48652bbe9d2a7979c72ae50db491023e48dce707 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 1 Jun 2021 21:43:13 +0200 Subject: [PATCH 26/29] Update EZFIO --- external/ezfio | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/external/ezfio b/external/ezfio index ccee52d0..ed1df9f3 160000 --- a/external/ezfio +++ b/external/ezfio @@ -1 +1 @@ -Subproject commit ccee52d00c2cde1d628b0d34f4a247143747bf36 +Subproject commit ed1df9f3c1f51752656ca98da5693a4119add05c From 54b01efc231562a99256fc8dbc27356fa786d8ab Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 7 Jun 2021 16:02:42 +0200 Subject: [PATCH 27/29] Removed do_csf : Bug --- src/davidson/diagonalize_ci.irp.f | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/davidson/diagonalize_ci.irp.f b/src/davidson/diagonalize_ci.irp.f index 06744934..da5fb950 100644 --- a/src/davidson/diagonalize_ci.irp.f +++ b/src/davidson/diagonalize_ci.irp.f @@ -57,8 +57,8 @@ END_PROVIDER enddo ! Deactivated temporarily: bug in N_csf - do_csf = s2_eig .and. only_expected_s2 .and. (expected_s2 == 0.d0) -! do_csf = .False. +! do_csf = s2_eig .and. only_expected_s2 .and. (expected_s2 == 0.d0) + do_csf = .False. if (diag_algorithm == "Davidson") then From c2bb6e92f0451b3a08f85ad5b2a62324f5d77b8e Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 7 Jun 2021 16:15:35 +0200 Subject: [PATCH 28/29] Protected internal writes --- src/zmq/utils.irp.f | 26 ++++++++++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/src/zmq/utils.irp.f b/src/zmq/utils.irp.f index 93dbd16a..7cb6c896 100644 --- a/src/zmq/utils.irp.f +++ b/src/zmq/utils.irp.f @@ -127,7 +127,9 @@ function zmq_port(ishift) END_DOC integer, intent(in) :: ishift character*(8) :: zmq_port + !$OMP CRITICAL(write) write(zmq_port,'(I8)') zmq_port_start+ishift + !$OMP END CRITICAL(write) zmq_port = adjustl(trim(zmq_port)) end @@ -518,7 +520,9 @@ subroutine new_parallel_job(zmq_to_qp_run_socket,zmq_socket_pull,name_in) zmq_to_qp_run_socket = new_zmq_to_qp_run_socket() zmq_socket_pull = new_zmq_pull_socket () + !$OMP CRITICAL(write) write(name,'(A,I8.8)') trim(name_in)//'.', icount + !$OMP END CRITICAL(write) sze = len(trim(name)) zmq_state = trim(name) call lowercase(name,sze) @@ -582,7 +586,9 @@ subroutine end_parallel_job(zmq_to_qp_run_socket,zmq_socket_pull,name_in) integer, save :: icount=0 icount = icount+1 + !$OMP CRITICAL(write) write(name,'(A,I8.8)') trim(name_in)//'.', icount + !$OMP END CRITICAL(write) sze = len(trim(name)) call lowercase(name,sze) if (name /= zmq_state) then @@ -704,7 +710,9 @@ integer function disconnect_from_taskserver_state(zmq_to_qp_run_socket, worker_i disconnect_from_taskserver_state = -1 + !$OMP CRITICAL(write) write(message,*) 'disconnect '//trim(state), worker_id + !$OMP END CRITICAL(write) sze = min(510,len(trim(message))) rc = f77_zmq_send(zmq_to_qp_run_socket, trim(message), sze, 0) @@ -781,7 +789,9 @@ integer function zmq_abort(zmq_to_qp_run_socket) character*(512) :: message zmq_abort = 0 + !$OMP CRITICAL(write) write(message,*) 'abort ' + !$OMP END CRITICAL(write) sze = len(trim(message)) @@ -823,7 +833,9 @@ integer function task_done_to_taskserver(zmq_to_qp_run_socket, worker_id, task_i task_done_to_taskserver = 0 + !$OMP CRITICAL(write) write(message,*) 'task_done '//trim(zmq_state), worker_id, task_id + !$OMP END CRITICAL(write) sze = len(trim(message)) rc = f77_zmq_send(zmq_to_qp_run_socket, trim(message), sze, 0) @@ -856,9 +868,11 @@ integer function tasks_done_to_taskserver(zmq_to_qp_run_socket, worker_id, task_ tasks_done_to_taskserver = 0 + !$OMP CRITICAL(write) allocate(character(LEN=64+n_tasks*12) :: message) write(fmt,*) '(A,X,A,I10,X,', n_tasks, '(I11,1X))' write(message,*) 'task_done '//trim(zmq_state), worker_id, (task_id(k), k=1,n_tasks) + !$OMP END CRITICAL(write) sze = len(trim(message)) rc = f77_zmq_send(zmq_to_qp_run_socket, trim(message), sze, 0) @@ -900,7 +914,9 @@ integer function get_task_from_taskserver(zmq_to_qp_run_socket,worker_id,task_id get_task_from_taskserver = 0 + !$OMP CRITICAL(write) write(message,*) 'get_task '//trim(zmq_state), worker_id + !$OMP END CRITICAL(write) sze = len(trim(message)) rc = f77_zmq_send(zmq_to_qp_run_socket, message, sze, 0) @@ -961,7 +977,9 @@ integer function get_tasks_from_taskserver(zmq_to_qp_run_socket,worker_id,task_i get_tasks_from_taskserver = 0 + !$OMP CRITICAL(write) write(message,'(A,A,X,I10,I10)') 'get_tasks ', trim(zmq_state), worker_id, n_tasks + !$OMP END CRITICAL(write) sze = len(trim(message)) rc = f77_zmq_send(zmq_to_qp_run_socket, message, sze, 0) @@ -1061,7 +1079,9 @@ integer function zmq_delete_task(zmq_to_qp_run_socket,zmq_socket_pull,task_id,mo zmq_delete_task = 0 + !$OMP CRITICAL(write) write(message,*) 'del_task ', zmq_state, task_id + !$OMP END CRITICAL(write) rc = f77_zmq_send(zmq_to_qp_run_socket,trim(message),len(trim(message)),0) if (rc /= len(trim(message))) then zmq_delete_task = -1 @@ -1101,7 +1121,9 @@ integer function zmq_delete_task_async_send(zmq_to_qp_run_socket,task_id,sending endif zmq_delete_task_async_send = 0 + !$OMP CRITICAL(write) write(message,*) 'del_task ', zmq_state, task_id + !$OMP END CRITICAL(write) rc = f77_zmq_send(zmq_to_qp_run_socket,trim(message),len(trim(message)),0) if (rc /= len(trim(message))) then zmq_delete_task_async_send = -1 @@ -1159,8 +1181,10 @@ integer function zmq_delete_tasks(zmq_to_qp_run_socket,zmq_socket_pull,task_id,n allocate(character(LEN=64+n_tasks*12) :: message) + !$OMP CRITICAL(write) write(fmt,*) '(A,1X,A,1X,', n_tasks, '(I11,1X))' write(message,*) 'del_task '//trim(zmq_state), (task_id(k), k=1,n_tasks) + !$OMP END CRITICAL(write) rc = f77_zmq_send(zmq_to_qp_run_socket,trim(message),len(trim(message)),0) @@ -1206,8 +1230,10 @@ integer function zmq_delete_tasks_async_send(zmq_to_qp_run_socket,task_id,n_task allocate(character(LEN=64+n_tasks*12) :: message) + !$OMP CRITICAL(write) write(fmt,*) '(A,1X,A,1X,', n_tasks, '(I11,1X))' write(message,*) 'del_task '//trim(zmq_state), (task_id(k), k=1,n_tasks) + !$OMP END CRITICAL(write) rc = f77_zmq_send(zmq_to_qp_run_socket,trim(message),len(trim(message)),0) From 0bf0513fb13623d3a2807c11a8316eea51a72705 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 11 Jun 2021 14:18:23 +0200 Subject: [PATCH 29/29] Added VERSION file --- VERSION | 1 + ocaml/qp_create_ezfio.ml | 27 ++++++++++++++++++++++++--- src/ao_basis/aos.irp.f | 10 +++++++--- src/basis/EZFIO.cfg | 10 ++++++++-- 4 files changed, 40 insertions(+), 8 deletions(-) create mode 100644 VERSION diff --git a/VERSION b/VERSION new file mode 100644 index 00000000..c043eea7 --- /dev/null +++ b/VERSION @@ -0,0 +1 @@ +2.2.1 diff --git a/ocaml/qp_create_ezfio.ml b/ocaml/qp_create_ezfio.ml index e4ead21b..a4865e2b 100644 --- a/ocaml/qp_create_ezfio.ml +++ b/ocaml/qp_create_ezfio.ml @@ -604,9 +604,30 @@ let run ?o b au c d m p cart xyz_file = ~rank:1 ~dim:[| shell_num |] ~data:ang_mom ) ; Ezfio.set_basis_shell_prim_index (Ezfio.ezfio_array_of_list ~rank:1 ~dim:[| shell_num |] ~data:shell_prim_idx) ; - Ezfio.set_basis_shell_nucl (Ezfio.ezfio_array_of_list - ~rank:1 ~dim:[| shell_num |] - ~data:(list_map (fun (_,n) -> Nucl_number.to_int n) basis) ) ; + Ezfio.set_basis_basis_nucleus_index (Ezfio.ezfio_array_of_list + ~rank:1 ~dim:[| nucl_num |] + ~data:( + list_map (fun (_,n) -> Nucl_number.to_int n) basis + |> List.fold_left (fun accu i -> + match accu with + | [] -> [] + | (h,j) :: rest -> if j == i then ((h+1,j)::rest) else ((h+1,i)::(h+1,j)::rest) + ) [(0,0)] + |> List.rev + |> List.map fst + )) ; + Ezfio.set_basis_nucleus_shell_num(Ezfio.ezfio_array_of_list + ~rank:1 ~dim:[| nucl_num |] + ~data:( + list_map (fun (_,n) -> Nucl_number.to_int n) basis + |> List.fold_left (fun accu i -> + match accu with + | [] -> [(1,i)] + | (h,j) :: rest -> if j == i then ((h+1,j)::rest) else ((1,i)::(h,j)::rest) + ) [] + |> List.rev + |> List.map fst + )) ; Ezfio.set_basis_prim_coef (Ezfio.ezfio_array_of_list ~rank:1 ~dim:[| prim_num |] ~data:coef) ; Ezfio.set_basis_prim_expo (Ezfio.ezfio_array_of_list diff --git a/src/ao_basis/aos.irp.f b/src/ao_basis/aos.irp.f index 2ff72898..17aa784d 100644 --- a/src/ao_basis/aos.irp.f +++ b/src/ao_basis/aos.irp.f @@ -40,9 +40,9 @@ END_PROVIDER do i=1,ao_num - powA(1) = ao_power(i,1) - powA(2) = ao_power(i,2) - powA(3) = ao_power(i,3) + powA(1) = ao_power(i,1) + ao_power(i,2) + ao_power(i,3) + powA(2) = 0 + powA(3) = 0 ! Normalization of the primitives if (primitives_normalized) then @@ -57,6 +57,10 @@ END_PROVIDER enddo endif + powA(1) = ao_power(i,1) + powA(2) = ao_power(i,2) + powA(3) = ao_power(i,3) + ! Normalization of the contracted basis functions if (ao_normalized) then norm = 0.d0 diff --git a/src/basis/EZFIO.cfg b/src/basis/EZFIO.cfg index 42123373..7f2ede4c 100644 --- a/src/basis/EZFIO.cfg +++ b/src/basis/EZFIO.cfg @@ -13,6 +13,12 @@ type: integer doc: Number of shells interface: ezfio, provider +[nucleus_shell_num] +type: integer +doc: Number of shells per nucleus +size: (nuclei.nucl_num) +interface: ezfio, provider + [shell_normalization_factor] type: double precision doc: Normalization factor applied to the whole shell, ex $1/\sqrt{ }$ @@ -37,10 +43,10 @@ doc: Max number of primitives in a shell size: (basis.shell_num) interface: ezfio, provider -[shell_nucl] +[basis_nucleus_index] type: integer doc: Index of the nucleus on which the shell is centered -size: (basis.shell_num) +size: (nuclei.nucl_num) interface: ezfio, provider [prim_normalization_factor]