diff --git a/config/gfortran.cfg b/config/gfortran.cfg index 291ea7fa..b713aaf0 100644 --- a/config/gfortran.cfg +++ b/config/gfortran.cfg @@ -10,7 +10,7 @@ # # [COMMON] -FC : gfortran -g -ffree-line-length-none -mavx -I . +FC : gfortran -g -ffree-line-length-none -I . -static-libgcc LAPACK_LIB : -llapack -lblas IRPF90 : irpf90 IRPF90_FLAGS : --ninja --align=32 diff --git a/configure b/configure index 1bd08ee2..7370d540 100755 --- a/configure +++ b/configure @@ -438,11 +438,12 @@ def create_ninja_and_rc(l_installed): print str_info("qp_root"), python_path = [join(QP_ROOT, "scripts"), join(QP_ROOT, "install")] - l_python = [join(QP_ROOT, "scripts")] + l_python = [join("${QP_ROOT}", "scripts")] for dir_ in python_path: for folder in os.listdir(dir_): path = join(dir_, folder) if os.path.isdir(path): + path = path.replace(QP_ROOT,"${QP_ROOT}") l_python.append(path) path_ezfio = find_path('ezfio', l_installed, var_for_qp_root=True) @@ -451,9 +452,9 @@ def create_ninja_and_rc(l_installed): l_rc = [ 'export QP_ROOT={0}'.format(QP_ROOT), - 'export QP_EZFIO={0}'.format(path_ezfio), - 'export IRPF90={0}'.format(path_irpf90), - 'export NINJA={0}'.format(path_ninja), + 'export QP_EZFIO={0}'.format(path_ezfio.replace(QP_ROOT,"${QP_ROOT}")), + 'export IRPF90={0}'.format(path_irpf90.replace(QP_ROOT,"${QP_ROOT}")), + 'export NINJA={0}'.format(path_ninja.replace(QP_ROOT,"${QP_ROOT}")), 'export QP_PYTHON={0}'.format(":".join(l_python)), "", 'export PYTHONPATH="${QP_EZFIO}":"${QP_PYTHON}":"${PYTHONPATH}"', 'export PATH="${QP_PYTHON}":"${QP_ROOT}"/bin:"${QP_ROOT}"/ocaml:"${PATH}"', diff --git a/ocaml/Atom.ml b/ocaml/Atom.ml index 2aafe644..832cfa5b 100644 --- a/ocaml/Atom.ml +++ b/ocaml/Atom.ml @@ -8,8 +8,8 @@ type t = coord : Point3d.t ; } with sexp -(** Read xyz coordinates of the atom with unit u *) -let of_string u s = +(** Read xyz coordinates of the atom *) +let of_string ~units s = let buffer = s |> String.split ~on:' ' |> List.filter ~f:(fun x -> x <> "") @@ -18,21 +18,21 @@ let of_string u s = | [ name; charge; x; y; z ] -> { element = Element.of_string name ; charge = Charge.of_string charge ; - coord = Point3d.of_string u (String.concat [x; y; z] ~sep:" ") + coord = Point3d.of_string ~units (String.concat [x; y; z] ~sep:" ") } | [ name; x; y; z ] -> let e = Element.of_string name in { element = e ; charge = Element.to_charge e; - coord = Point3d.of_string u (String.concat [x; y; z] ~sep:" ") + coord = Point3d.of_string ~units (String.concat [x; y; z] ~sep:" ") } | _ -> raise (AtomError s) ;; -let to_string u a = +let to_string ~units a = [ Element.to_string a.element ; Charge.to_string a.charge ; - Point3d.to_string u a.coord ] + Point3d.to_string ~units a.coord ] |> String.concat ~sep:" " ;; diff --git a/ocaml/Atom.mli b/ocaml/Atom.mli index 93f7d885..28915993 100644 --- a/ocaml/Atom.mli +++ b/ocaml/Atom.mli @@ -5,5 +5,5 @@ type t = { element : Element.t; charge : Charge.t; coord : Point3d.t; } val t_of_sexp : Sexplib.Sexp.t -> t val sexp_of_t : t -> Sexplib.Sexp.t -val of_string : Units.units -> string -> t -val to_string : Units.units -> t -> string +val of_string : units:Units.units -> string -> t +val to_string : units:Units.units -> t -> string diff --git a/ocaml/Element.ml b/ocaml/Element.ml index d282340b..6bc2de4e 100644 --- a/ocaml/Element.ml +++ b/ocaml/Element.ml @@ -1,4 +1,5 @@ -open Core.Std;; +open Core.Std +open Qptypes exception ElementError of string @@ -8,49 +9,49 @@ type t = |Li|Be |B |C |N |O |F |Ne |Na|Mg |Al|Si|P |S |Cl|Ar |K |Ca|Sc|Ti|V |Cr|Mn|Fe|Co|Ni|Cu|Zn|Ga|Ge|As|Se|Br|Kr -with sexp;; +with sexp let of_string x = match (String.capitalize (String.lowercase x)) with -| "X" | "Dummy" -> X -| "H" | "Hydrogen" -> H -| "He" | "Helium" -> He -| "Li" | "Lithium" -> Li -| "Be" | "Beryllium" -> Be -| "B" | "Boron" -> B -| "C" | "Carbon" -> C -| "N" | "Nitrogen" -> N -| "O" | "Oxygen" -> O -| "F" | "Fluorine" -> F -| "Ne" | "Neon" -> Ne -| "Na" | "Sodium" -> Na -| "Mg" | "Magnesium" -> Mg -| "Al" | "Aluminum" -> Al -| "Si" | "Silicon" -> Si -| "P" | "Phosphorus" -> P -| "S" | "Sulfur" -> S -| "Cl" | "Chlorine" -> Cl -| "Ar" | "Argon" -> Ar -| "K" | "Potassium" -> K -| "Ca" | "Calcium" -> Ca -| "Sc" | "Scandium" -> Sc -| "Ti" | "Titanium" -> Ti -| "V" | "Vanadium" -> V -| "Cr" | "Chromium" -> Cr -| "Mn" | "Manganese" -> Mn -| "Fe" | "Iron" -> Fe -| "Co" | "Cobalt" -> Co -| "Ni" | "Nickel" -> Ni -| "Cu" | "Copper" -> Cu -| "Zn" | "Zinc" -> Zn -| "Ga" | "Gallium" -> Ga -| "Ge" | "Germanium" -> Ge -| "As" | "Arsenic" -> As -| "Se" | "Selenium" -> Se -| "Br" | "Bromine" -> Br -| "Kr" | "Krypton" -> Kr +| "X" | "Dummy" -> X +| "H" | "Hydrogen" -> H +| "He" | "Helium" -> He +| "Li" | "Lithium" -> Li +| "Be" | "Beryllium" -> Be +| "B" | "Boron" -> B +| "C" | "Carbon" -> C +| "N" | "Nitrogen" -> N +| "O" | "Oxygen" -> O +| "F" | "Fluorine" -> F +| "Ne" | "Neon" -> Ne +| "Na" | "Sodium" -> Na +| "Mg" | "Magnesium" -> Mg +| "Al" | "Aluminum" -> Al +| "Si" | "Silicon" -> Si +| "P" | "Phosphorus" -> P +| "S" | "Sulfur" -> S +| "Cl" | "Chlorine" -> Cl +| "Ar" | "Argon" -> Ar +| "K" | "Potassium" -> K +| "Ca" | "Calcium" -> Ca +| "Sc" | "Scandium" -> Sc +| "Ti" | "Titanium" -> Ti +| "V" | "Vanadium" -> V +| "Cr" | "Chromium" -> Cr +| "Mn" | "Manganese" -> Mn +| "Fe" | "Iron" -> Fe +| "Co" | "Cobalt" -> Co +| "Ni" | "Nickel" -> Ni +| "Cu" | "Copper" -> Cu +| "Zn" | "Zinc" -> Zn +| "Ga" | "Gallium" -> Ga +| "Ge" | "Germanium" -> Ge +| "As" | "Arsenic" -> As +| "Se" | "Selenium" -> Se +| "Br" | "Bromine" -> Br +| "Kr" | "Krypton" -> Kr | x -> raise (ElementError ("Element "^x^" unknown")) -;; + let to_string = function | X -> "X" @@ -90,7 +91,7 @@ let to_string = function | Se -> "Se" | Br -> "Br" | Kr -> "Kr" -;; + let to_long_string = function | X -> "Dummy" @@ -130,7 +131,7 @@ let to_long_string = function | Se -> "Selenium" | Br -> "Bromine" | Kr -> "Krypton" -;; + let to_charge c = let result = match c with @@ -172,7 +173,7 @@ let to_charge c = | Br -> 35 | Kr -> 36 in Charge.of_int result -;; + let of_charge c = match (Charge.to_int c) with | 0 -> X @@ -213,5 +214,134 @@ let of_charge c = match (Charge.to_int c) with | 35 -> Br | 36 -> Kr | x -> raise (ElementError ("Element of charge "^(string_of_int x)^" unknown")) -;; + + +let covalent_radius x = + let result = function + | X -> 0. + | H -> 0.37 + | He -> 0.70 + | Li -> 1.23 + | Be -> 0.89 + | B -> 0.90 + | C -> 0.85 + | N -> 0.74 + | O -> 0.74 + | F -> 0.72 + | Ne -> 0.70 + | Na -> 1.00 + | Mg -> 1.36 + | Al -> 1.25 + | Si -> 1.17 + | P -> 1.10 + | S -> 1.10 + | Cl -> 0.99 + | Ar -> 0.70 + | K -> 2.03 + | Ca -> 1.74 + | Sc -> 1.44 + | Ti -> 1.32 + | V -> 1.22 + | Cr -> 0.00 + | Mn -> 1.16 + | Fe -> 0.00 + | Co -> 1.15 + | Ni -> 1.17 + | Cu -> 1.25 + | Zn -> 1.25 + | Ga -> 1.20 + | Ge -> 1.21 + | As -> 1.16 + | Se -> 0.70 + | Br -> 1.24 + | Kr -> 1.91 + in + Units.angstrom_to_bohr *. (result x) + |> Positive_float.of_float + +let vdw_radius x = + let result = function + | X -> 0. + | H -> 1.20 + | He -> 1.70 + | Li -> 1.70 + | Be -> 1.70 + | B -> 1.70 + | C -> 1.70 + | N -> 1.55 + | O -> 1.52 + | F -> 1.47 + | Ne -> 1.70 + | Na -> 1.70 + | Mg -> 1.70 + | Al -> 1.94 + | Si -> 2.10 + | P -> 1.80 + | S -> 1.80 + | Cl -> 1.75 + | Ar -> 1.70 + | K -> 1.70 + | Ca -> 1.70 + | Sc -> 1.70 + | Ti -> 1.70 + | V -> 1.98 + | Cr -> 1.94 + | Mn -> 1.93 + | Fe -> 1.93 + | Co -> 1.92 + | Ni -> 1.70 + | Cu -> 1.70 + | Zn -> 1.70 + | Ga -> 2.02 + | Ge -> 1.70 + | As -> 1.96 + | Se -> 1.70 + | Br -> 2.10 + | Kr -> 1.70 + in + Units.angstrom_to_bohr *. (result x) + |> Positive_float.of_float + +let mass x = + let result = function + | X -> 0. + | H -> 1.0079 + | He -> 4.00260 + | Li -> 6.941 + | Be -> 9.01218 + | B -> 10.81 + | C -> 12.011 + | N -> 14.0067 + | O -> 15.9994 + | F -> 18.998403 + | Ne -> 20.179 + | Na -> 22.98977 + | Mg -> 24.305 + | Al -> 26.98154 + | Si -> 28.0855 + | P -> 30.97376 + | S -> 32.06 + | Cl -> 35.453 + | Ar -> 39.948 + | K -> 39.0983 + | Ca -> 40.08 + | Sc -> 44.9559 + | Ti -> 47.90 + | V -> 50.9415 + | Cr -> 51.996 + | Mn -> 54.9380 + | Fe -> 55.9332 + | Co -> 58.9332 + | Ni -> 58.70 + | Cu -> 63.546 + | Zn -> 65.38 + | Ga -> 69.72 + | Ge -> 72.59 + | As -> 74.9216 + | Se -> 78.96 + | Br -> 79.904 + | Kr -> 83.80 + in + result x + |> Positive_float.of_float diff --git a/ocaml/Element.mli b/ocaml/Element.mli index 48bd3c61..8d9862c9 100644 --- a/ocaml/Element.mli +++ b/ocaml/Element.mli @@ -13,6 +13,8 @@ val of_string : string -> t val to_string : t -> string val to_long_string : t -> string -(** get the positive charge *) +(** Properties *) val to_charge : t -> Charge.t val of_charge : Charge.t -> t +val covalent_radius : t -> Qptypes.Positive_float.t +val vdw_radius : t -> Qptypes.Positive_float.t diff --git a/ocaml/Input_nuclei.ml b/ocaml/Input_nuclei.ml index cbcb9f46..d050ded9 100644 --- a/ocaml/Input_nuclei.ml +++ b/ocaml/Input_nuclei.ml @@ -147,7 +147,7 @@ nucl_coord = %s (b.nucl_charge |> Array.to_list |> List.map ~f:(Charge.to_string) |> String.concat ~sep:", " ) (b.nucl_coord |> Array.to_list |> List.map - ~f:(Point3d.to_string Units.Bohr) |> String.concat ~sep:"\n" ) + ~f:(Point3d.to_string ~units:Units.Bohr) |> String.concat ~sep:"\n" ) ;; @@ -161,7 +161,7 @@ nucl_coord = %s Printf.sprintf " %-3s %d %s" (b.nucl_label.(i) |> Element.to_string) (b.nucl_charge.(i) |> Charge.to_int ) - (b.nucl_coord.(i) |> Point3d.to_string Units.Angstrom) ) + (b.nucl_coord.(i) |> Point3d.to_string ~units:Units.Angstrom) ) ) |> String.concat ~sep:"\n" in Printf.sprintf " diff --git a/ocaml/Molecule.ml b/ocaml/Molecule.ml index f295420b..f0800f7f 100644 --- a/ocaml/Molecule.ml +++ b/ocaml/Molecule.ml @@ -10,27 +10,36 @@ type t = { } with sexp let get_charge { nuclei ; elec_alpha ; elec_beta } = - let result = (Elec_alpha_number.to_int elec_alpha) + - (Elec_beta_number.to_int elec_beta) in + let result = + (Elec_alpha_number.to_int elec_alpha) + + (Elec_beta_number.to_int elec_beta) + in let rec nucl_charge = function | a::rest -> (Charge.to_float a.Atom.charge) +. nucl_charge rest | [] -> 0. in Charge.of_float (nucl_charge nuclei -. (Float.of_int result)) -;; + let get_multiplicity m = - let elec_alpha = m.elec_alpha in + let elec_alpha = + m.elec_alpha + in Multiplicity.of_alpha_beta elec_alpha m.elec_beta -;; + let get_nucl_num m = - let nmax = (List.length m.nuclei) in + let nmax = + List.length m.nuclei + in Nucl_number.of_int nmax ~max:nmax -;; + let name m = - let cm = Charge.to_int (get_charge m) in + let cm = + get_charge m + |> Charge.to_int + in let c = match cm with | 0 -> "" @@ -39,8 +48,12 @@ let name m = | i when i>1 -> Printf.sprintf " (%d+)" i | i -> Printf.sprintf " (%d-)" (-i) in - let mult = Multiplicity.to_string (get_multiplicity m) in - let { nuclei ; elec_alpha ; elec_beta } = m in + let mult = + get_multiplicity m + |> Multiplicity.to_string + in + let { nuclei ; elec_alpha ; elec_beta } = m + in let rec build_list accu = function | a::rest -> begin @@ -53,7 +66,9 @@ let name m = in let rec build_name accu = function | (a, n)::rest -> - let a = Element.to_string a in + let a = + Element.to_string a + in begin match n with | 1 -> build_name (a::accu) rest @@ -64,19 +79,25 @@ let name m = end | [] -> accu in - let result = build_list [] nuclei |> build_name [c ; ", " ; mult] + let result = + build_list [] nuclei |> build_name [c ; ", " ; mult] in String.concat (result) -;; + let to_string m = - let { nuclei ; elec_alpha ; elec_beta } = m in - let n = List.length nuclei in - let title = name m in - [ Int.to_string n ; title ] @ (List.map ~f:(fun x -> Atom.to_string - Units.Angstrom x) nuclei) + let { nuclei ; elec_alpha ; elec_beta } = m + in + let n = + List.length nuclei + in + let title = + name m + in + [ Int.to_string n ; title ] @ + (List.map ~f:(fun x -> Atom.to_string Units.Angstrom x) nuclei) |> String.concat ~sep:"\n" -;; + let of_xyz_string ?(charge=(Charge.of_int 0)) ?(multiplicity=(Multiplicity.of_int 1)) @@ -94,7 +115,9 @@ let of_xyz_string ) + 1 - (Charge.to_int charge) |> Elec_number.of_int in - let (na,nb) = Multiplicity.to_alpha_beta ne multiplicity in + let (na,nb) = + Multiplicity.to_alpha_beta ne multiplicity + in let result = { nuclei = l ; elec_alpha = na ; @@ -109,7 +132,7 @@ let of_xyz_string raise (MultiplicityError msg); else () ; result -;; + let of_xyz_file @@ -121,8 +144,33 @@ let of_xyz_file let (_,buffer) = String.lsplit2_exn buffer ~on:'\n' in of_xyz_string ~charge:charge ~multiplicity:multiplicity ~units:units buffer -;; -include To_md5;; + + +let distance_matrix molecule = + let coord = + molecule.nuclei + |> List.map ~f:(fun x -> x.Atom.coord) + |> Array.of_list + in + let n = + Array.length coord + in + let result = + Array.make_matrix ~dimx:n ~dimy:n 0. + in + for i = 0 to (n-1) + do + for j = 0 to (n-1) + do + result.(i).(j) <- Point3d.distance coord.(i) coord.(j) + done; + done; + result + + + + +include To_md5 let to_md5 = to_md5 sexp_of_t -;; + diff --git a/ocaml/Molecule.mli b/ocaml/Molecule.mli index 176441d4..1a3d9715 100644 --- a/ocaml/Molecule.mli +++ b/ocaml/Molecule.mli @@ -34,5 +34,9 @@ val of_xyz_string : ?multiplicity:Multiplicity.t -> ?units:Units.units -> string -> t +(** Creates the distance matrix between all the atoms *) +val distance_matrix : + t -> (float array) array + (** Computes the MD5 hash *) val to_md5 : t -> Qptypes.MD5.t diff --git a/ocaml/Point3d.ml b/ocaml/Point3d.ml index 18f091f1..5717ed39 100644 --- a/ocaml/Point3d.ml +++ b/ocaml/Point3d.ml @@ -7,9 +7,16 @@ type t = { z : float ; } with sexp +let of_tuple ~units (x,y,z) = + let f = match units with + | Units.Bohr -> 1. + | Units.Angstrom -> Units.angstrom_to_bohr + in + { x = x *. f ; y = y *. f ; z = z *. f } + (** Read x y z coordinates in string s with units u *) -let of_string u s = - let f = match u with +let of_string ~units s = + let f = match units with | Units.Bohr -> 1. | Units.Angstrom -> Units.angstrom_to_bohr in @@ -22,7 +29,6 @@ let of_string u s = { x = l.(0) *. f ; y = l.(1) *. f ; z = l.(2) *. f } -;; let distance2 p1 p2 = @@ -30,17 +36,18 @@ let distance2 p1 p2 = and { x=x2 ; y=y2 ; z=z2 } = p2 in (x2-.x1)*.(x2-.x1) +. (y2-.y1)*.(y2-.y1) +. (z2-.z1)*.(z2-.z1) |> Positive_float.of_float -;; -let distance p1 p2 = sqrt (Positive_float.to_float (distance2 p1 p2)) -;; -let to_string u p = - let f = match u with +let distance p1 p2 = + sqrt (Positive_float.to_float (distance2 p1 p2)) + + +let to_string ~units p = + let f = match units with | Units.Bohr -> 1. | Units.Angstrom -> Units.bohr_to_angstrom in let { x=x ; y=y ; z=z } = p in Printf.sprintf "%16.8f %16.8f %16.8f" (x*.f) (y*.f) (z*.f) -;; + diff --git a/ocaml/Point3d.mli b/ocaml/Point3d.mli index 066a4514..6d7428ec 100644 --- a/ocaml/Point3d.mli +++ b/ocaml/Point3d.mli @@ -4,11 +4,14 @@ type t = z : float; } with sexp +(** Create from a tuple of floats *) +val of_tuple : units:Units.units -> float*float*float -> t + (** Create from an xyz string *) -val of_string : Units.units -> string -> t +val of_string : units:Units.units -> string -> t (** Convert to a string for printing *) -val to_string : Units.units -> t -> string +val to_string : units:Units.units -> t -> string (** Computes the squared distance between 2 points *) val distance2 : t -> t -> Qptypes.Positive_float.t diff --git a/ocaml/qp_create_ezfio_from_xyz.ml b/ocaml/qp_create_ezfio_from_xyz.ml index 544e6e09..fb10ff2f 100644 --- a/ocaml/qp_create_ezfio_from_xyz.ml +++ b/ocaml/qp_create_ezfio_from_xyz.ml @@ -11,26 +11,92 @@ let spec = ~doc:"string Name of basis set." +> flag "c" (optional_with_default 0 int) ~doc:"int Total charge of the molecule. Default is 0." + +> flag "d" (optional_with_default 0. float) + ~doc:"float Add dummy atoms. x * (covalent radii of the atoms)" +> flag "m" (optional_with_default 1 int) ~doc:"int Spin multiplicity (2S+1) of the molecule. Default is 1." +> flag "p" no_arg ~doc:" Using pseudopotentials" +> anon ("xyz_file" %: string) -;; -let run ?o b c m p xyz_file = + +let dummy_centers ~threshold ~molecule ~nuclei = + let d = + Molecule.distance_matrix molecule + in + let n = + Array.length d + in + let nuclei = + Array.of_list nuclei + in + let rec aux accu = function + | (-1,_) -> accu + | (i,-1) -> aux accu (i-1,i-1) + | (i,j) when (i>j) -> + let new_accu = + let x,y = + Element.covalent_radius (nuclei.(i)).Atom.element |> Positive_float.to_float, + Element.covalent_radius (nuclei.(j)).Atom.element |> Positive_float.to_float + in + let r = + ( x +. y ) *. threshold + in + if d.(i).(j) < r then + (i,x,j,y,d.(i).(j)) :: accu + else + accu + in aux new_accu (i,j-1) + | (i,j) when (i=j) -> aux accu (i,j-1) + | _ -> assert false + in + aux [] (n-1,n-1) + |> List.map ~f:(fun (i,x,j,y,r) -> + let f = + x /. (x +. y) + in + let u = + Point3d.of_tuple ~units:Units.Bohr + ( nuclei.(i).Atom.coord.Point3d.x +. + (nuclei.(j).Atom.coord.Point3d.x -. nuclei.(i).Atom.coord.Point3d.x) *. f, + nuclei.(i).Atom.coord.Point3d.y +. + (nuclei.(j).Atom.coord.Point3d.y -. nuclei.(i).Atom.coord.Point3d.y) *. f, + nuclei.(i).Atom.coord.Point3d.z +. + (nuclei.(j).Atom.coord.Point3d.z -. nuclei.(i).Atom.coord.Point3d.z) *. f) + in + Atom.{ element = Element.X ; charge = Charge.of_int 0 ; coord = u } + ) + + + + + +let run ?o b c d m p xyz_file = (* Read molecule *) let molecule = (Molecule.of_xyz_file xyz_file ~charge:(Charge.of_int c) ~multiplicity:(Multiplicity.of_int m) ) in - let nuclei = molecule.Molecule.nuclei in + let dummy = + dummy_centers ~threshold:d ~molecule ~nuclei:molecule.Molecule.nuclei + in +(* + List.iter dummy ~f:(fun x -> + Printf.printf "%s\n" (Atom.to_string ~units:Units.Angstrom x) + ); +*) + let nuclei = + molecule.Molecule.nuclei @ dummy + in + let basis_table = Hashtbl.Poly.create () in (* Open basis set channels *) let basis_channel element = - let key = Element.to_string element in + let key = + Element.to_string element + in match Hashtbl.find basis_table key with | Some in_channel -> in_channel @@ -80,7 +146,8 @@ let run ?o b c m p xyz_file = in Unix.unlink filename; List.iter nuclei ~f:(fun elem-> - let key = Element.to_string elem.Atom.element + let key = + Element.to_string elem.Atom.element in match Hashtbl.add basis_table ~key:key ~data:new_channel with | `Ok -> () @@ -92,7 +159,8 @@ let run ?o b c m p xyz_file = let elem = Element.of_string key and basis = String.lowercase basis in - let key = Element.to_string elem + let key = + Element.to_string elem in let command = Qpackage.root ^ "/scripts/get_basis.sh \"" ^ temp_filename ^ @@ -175,7 +243,12 @@ let run ?o b c m p xyz_file = |> List.rev |> List.map ~f:(fun (x,i) -> try - Basis.read_element (basis_channel x.Atom.element) i x.Atom.element + let e = + match x.Atom.element with + | Element.X -> Element.H + | e -> e + in + Basis.read_element (basis_channel x.Atom.element) i e with | End_of_file -> begin @@ -264,7 +337,7 @@ let run ?o b c m p xyz_file = | None -> failwith "Error in basis" | Some x -> Input.Ao_basis.write x -;; + let command = Command.basic @@ -279,13 +352,13 @@ elements can be defined as follows: ") spec - (fun o b c m p xyz_file () -> - run ?o b c m p xyz_file ) -;; + (fun o b c d m p xyz_file () -> + run ?o b c d m p xyz_file ) + let () = Command.run command -;; + diff --git a/ocaml/test_molecule.ml b/ocaml/test_molecule.ml index 3abd7e9a..5d2935ed 100644 --- a/ocaml/test_molecule.ml +++ b/ocaml/test_molecule.ml @@ -39,6 +39,15 @@ H 0.54386314 0.00000000 -0.92559535 let m = Molecule.of_xyz_file "c2h6.xyz" in print_string (Molecule.to_string m); + print_string "\nDistance matrix\n"; + print_string "---------------\n"; + let d = + Molecule.distance_matrix m + in + Array.iter d ~f:(fun x -> + Array.iter x ~f:(fun y -> Printf.printf "%12.8f " y); + print_newline (); + ) ;; test_molecule ();; diff --git a/plugins/Hartree_Fock/EZFIO.cfg b/plugins/Hartree_Fock/EZFIO.cfg index c39c3483..1b8486ee 100644 --- a/plugins/Hartree_Fock/EZFIO.cfg +++ b/plugins/Hartree_Fock/EZFIO.cfg @@ -10,6 +10,12 @@ doc: Maximum number of SCF iterations interface: ezfio,provider,ocaml default: 200 +[level_shift] +type: Positive_float +doc: Energy shift on the virtual MOs to improve SCF convergence +interface: ezfio,provider,ocaml +default: 0.0 + [mo_guess_type] type: MO_guess doc: Initial MO guess. Can be [ Huckel | HCore ] diff --git a/plugins/Hartree_Fock/Fock_matrix.irp.f b/plugins/Hartree_Fock/Fock_matrix.irp.f index 2561ad03..82540de3 100644 --- a/plugins/Hartree_Fock/Fock_matrix.irp.f +++ b/plugins/Hartree_Fock/Fock_matrix.irp.f @@ -73,8 +73,13 @@ enddo endif + ! Introduce level shift here + do i = elec_alpha_num+1, mo_tot_num + Fock_matrix_mo(i,i) += level_shift + enddo + do i = 1, mo_tot_num - Fock_matrix_diag_mo(i) = Fock_matrix_mo(i,i) + Fock_matrix_diag_mo(i) = Fock_matrix_mo(i,i) enddo END_PROVIDER @@ -108,9 +113,10 @@ END_PROVIDER END_DOC integer :: i,j,k,l,k1,r,s + integer :: i0,j0,k0,l0 integer*8 :: p,q - double precision :: integral - double precision :: ao_bielec_integral + double precision :: integral, c0, c1, c2 + double precision :: ao_bielec_integral, local_threshold double precision, allocatable :: ao_bi_elec_integral_alpha_tmp(:,:) double precision, allocatable :: ao_bi_elec_integral_beta_tmp(:,:) !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: ao_bi_elec_integral_beta_tmp @@ -121,11 +127,12 @@ END_PROVIDER if (do_direct_integrals) then !$OMP PARALLEL DEFAULT(NONE) & - !$OMP PRIVATE(i,j,l,k1,k,integral,ii,jj,kk,ll,i8,keys,values,p,q,r,s, & - !$OMP ao_bi_elec_integral_alpha_tmp,ao_bi_elec_integral_beta_tmp)& + !$OMP PRIVATE(i,j,l,k1,k,integral,ii,jj,kk,ll,i8,keys,values,p,q,r,s,i0,j0,k0,l0, & + !$OMP ao_bi_elec_integral_alpha_tmp,ao_bi_elec_integral_beta_tmp, c0, c1, c2, & + !$OMP local_threshold)& !$OMP SHARED(ao_num,ao_num_align,HF_density_matrix_ao_alpha,HF_density_matrix_ao_beta,& !$OMP ao_integrals_map,ao_integrals_threshold, ao_bielec_integral_schwartz, & - !$OMP ao_overlap_abs, ao_bi_elec_integral_alpha, ao_bi_elec_integral_beta) + !$OMP ao_overlap_abs, ao_bi_elec_integral_alpha, ao_bi_elec_integral_beta) allocate(keys(1), values(1)) allocate(ao_bi_elec_integral_alpha_tmp(ao_num_align,ao_num), & @@ -152,14 +159,16 @@ END_PROVIDER < ao_integrals_threshold) then cycle endif - if (ao_bielec_integral_schwartz(k,l)*ao_bielec_integral_schwartz(i,j) & - < ao_integrals_threshold) then - cycle - endif - values(1) = ao_bielec_integral(k,l,i,j) - if (abs(values(1)) < ao_integrals_threshold) then + local_threshold = ao_bielec_integral_schwartz(k,l)*ao_bielec_integral_schwartz(i,j) + if (local_threshold < ao_integrals_threshold) then cycle endif + i0 = i + j0 = j + k0 = k + l0 = l + values(1) = 0.d0 + local_threshold = ao_integrals_threshold/local_threshold do k2=1,8 if (kk(k2)==0) then cycle @@ -168,12 +177,21 @@ END_PROVIDER j = jj(k2) k = kk(k2) l = ll(k2) - integral = (HF_density_matrix_ao_alpha(k,l)+HF_density_matrix_ao_beta(k,l)) * values(1) + c0 = HF_density_matrix_ao_alpha(k,l)+HF_density_matrix_ao_beta(k,l) + c1 = HF_density_matrix_ao_alpha(k,i) + c2 = HF_density_matrix_ao_beta(k,i) + if ( dabs(c0)+dabs(c1)+dabs(c2) < local_threshold) then + cycle + endif + if (values(1) == 0.d0) then + values(1) = ao_bielec_integral(k0,l0,i0,j0) + endif + integral = c0 * values(1) ao_bi_elec_integral_alpha_tmp(i,j) += integral ao_bi_elec_integral_beta_tmp (i,j) += integral integral = values(1) - ao_bi_elec_integral_alpha_tmp(l,j) -= HF_density_matrix_ao_alpha(k,i) * integral - ao_bi_elec_integral_beta_tmp (l,j) -= HF_density_matrix_ao_beta (k,i) * integral + ao_bi_elec_integral_alpha_tmp(l,j) -= c1 * integral + ao_bi_elec_integral_beta_tmp (l,j) -= c2 * integral enddo enddo !$OMP END DO NOWAIT @@ -315,7 +333,9 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_ao, (ao_num_align, ao_num) ] ! Fock matrix in AO basis set END_DOC - if (elec_alpha_num == elec_beta_num) then + if ( (elec_alpha_num == elec_beta_num).and. & + (level_shift == 0.) ) & + then integer :: i,j do j=1,ao_num !DIR$ VECTOR ALIGNED diff --git a/plugins/Hartree_Fock/SCF.irp.f b/plugins/Hartree_Fock/SCF.irp.f index 864e9f3f..dead61ee 100644 --- a/plugins/Hartree_Fock/SCF.irp.f +++ b/plugins/Hartree_Fock/SCF.irp.f @@ -42,16 +42,13 @@ subroutine run BEGIN_DOC ! Run SCF calculation END_DOC - double precision :: SCF_energy_before,SCF_energy_after,diag_H_mat_elem,get_mo_bielec_integral + double precision :: SCF_energy_before,SCF_energy_after,diag_H_mat_elem double precision :: E0 integer :: i_it, i, j, k E0 = HF_energy - thresh_SCF = 1.d-10 - call damping_SCF mo_label = "Canonical" - TOUCH mo_label mo_coef - call save_mos + call damping_SCF end diff --git a/plugins/Hartree_Fock/damping_SCF.irp.f b/plugins/Hartree_Fock/damping_SCF.irp.f index d7d9c2bf..c7c005c9 100644 --- a/plugins/Hartree_Fock/damping_SCF.irp.f +++ b/plugins/Hartree_Fock/damping_SCF.irp.f @@ -86,7 +86,7 @@ subroutine damping_SCF if ((E_half > E).and.(E_new < E)) then lambda = 1.d0 exit - else if ((E_half > E).and.(lambda > 5.d-2)) then + else if ((E_half > E).and.(lambda > 5.d-4)) then lambda = 0.5d0 * lambda E_new = E_half else diff --git a/plugins/Perturbation/Moller_plesset.irp.f b/plugins/Perturbation/Moller_plesset.irp.f index 2e8ba8e1..7435f70c 100644 --- a/plugins/Perturbation/Moller_plesset.irp.f +++ b/plugins/Perturbation/Moller_plesset.irp.f @@ -31,8 +31,7 @@ subroutine pt2_moller_plesset(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_s (Fock_matrix_diag_mo(p1) + Fock_matrix_diag_mo(p2)) delta_e = 1.d0/delta_e - !call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det,psi_selectors_size,n_st,i_H_psi_array) - call i_H_psi_nominilist(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det,psi_selectors_size,n_st,i_H_psi_array) + call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det,psi_selectors_size,n_st,i_H_psi_array) h = diag_H_mat_elem(det_pert,Nint) do i =1,n_st H_pert_diag(i) = h diff --git a/plugins/Perturbation/delta_rho_perturbation.irp.f b/plugins/Perturbation/delta_rho_perturbation.irp.f index 77972c88..c95972a6 100644 --- a/plugins/Perturbation/delta_rho_perturbation.irp.f +++ b/plugins/Perturbation/delta_rho_perturbation.irp.f @@ -51,7 +51,7 @@ subroutine pt2_delta_rho_one_point(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,nde call i_O1_psi_alpha_beta(mo_integrated_delta_rho_one_point,det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_O1_psi_array) !call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array) - call i_H_psi(det_pert,minilist,idx_minilist,N_minilist,psi_selectors_coef,Nint,N_minilist,psi_selectors_size,N_st,i_H_psi_array) + call i_H_psi_minilist(det_pert,minilist,idx_minilist,N_minilist,psi_selectors_coef,Nint,N_minilist,psi_selectors_size,N_st,i_H_psi_array) h = diag_H_mat_elem(det_pert,Nint) oii = diag_O1_mat_elem_alpha_beta(mo_integrated_delta_rho_one_point,det_pert,N_int) diff --git a/plugins/Perturbation/dipole_moment.irp.f b/plugins/Perturbation/dipole_moment.irp.f index 4af9ea6b..53beb081 100644 --- a/plugins/Perturbation/dipole_moment.irp.f +++ b/plugins/Perturbation/dipole_moment.irp.f @@ -51,7 +51,7 @@ subroutine pt2_dipole_moment_z(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_ call i_O1_psi(mo_dipole_z,det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_O1_psi_array) !call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array) - call i_H_psi(det_pert,minilist,idx_minilist,N_minilist,psi_selectors_coef,Nint,N_minilist,psi_selectors_size,N_st,i_H_psi_array) + call i_H_psi_minilist(det_pert,minilist,idx_minilist,N_minilist,psi_selectors_coef,Nint,N_minilist,psi_selectors_size,N_st,i_H_psi_array) h = diag_H_mat_elem(det_pert,Nint) oii = diag_O1_mat_elem(mo_dipole_z,det_pert,N_int) diff --git a/plugins/Perturbation/epstein_nesbet.irp.f b/plugins/Perturbation/epstein_nesbet.irp.f index b4ad9bfe..ede75df9 100644 --- a/plugins/Perturbation/epstein_nesbet.irp.f +++ b/plugins/Perturbation/epstein_nesbet.irp.f @@ -28,7 +28,7 @@ subroutine pt2_epstein_nesbet(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_s ASSERT (Nint == N_int) ASSERT (Nint > 0) !call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array) - call i_H_psi(det_pert,minilist,idx_minilist,N_minilist,psi_selectors_coef,Nint,N_minilist,psi_selectors_size,N_st,i_H_psi_array) + call i_H_psi_minilist(det_pert,minilist,idx_minilist,N_minilist,psi_selectors_coef,Nint,N_minilist,psi_selectors_size,N_st,i_H_psi_array) h = diag_H_mat_elem(det_pert,Nint) @@ -79,7 +79,7 @@ subroutine pt2_epstein_nesbet_2x2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet PROVIDE CI_electronic_energy !call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array) - call i_H_psi(det_pert,minilist,idx_minilist,N_minilist,psi_selectors_coef,Nint,N_minilist,psi_selectors_size,N_st,i_H_psi_array) + call i_H_psi_minilist(det_pert,minilist,idx_minilist,N_minilist,psi_selectors_coef,Nint,N_minilist,psi_selectors_size,N_st,i_H_psi_array) h = diag_H_mat_elem(det_pert,Nint) do i =1,N_st diff --git a/plugins/Perturbation/pert_sc2.irp.f b/plugins/Perturbation/pert_sc2.irp.f index 15399e4e..1ae6ce73 100644 --- a/plugins/Perturbation/pert_sc2.irp.f +++ b/plugins/Perturbation/pert_sc2.irp.f @@ -221,7 +221,7 @@ subroutine pt2_epstein_nesbet_sc2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet ASSERT (Nint == N_int) ASSERT (Nint > 0) !call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array) - call i_H_psi(det_pert,minilist,idx_minilist,N_minilist,psi_selectors_coef,Nint,N_minilist,psi_selectors_size,N_st,i_H_psi_array) + call i_H_psi_minilist(det_pert,minilist,idx_minilist,N_minilist,psi_selectors_coef,Nint,N_minilist,psi_selectors_size,N_st,i_H_psi_array) h = diag_H_mat_elem(det_pert,Nint) diff --git a/plugins/QmcChem/target_pt2_qmc.irp.f b/plugins/QmcChem/target_pt2_qmc.irp.f index e5142a73..228bcb5e 100644 --- a/plugins/QmcChem/target_pt2_qmc.irp.f +++ b/plugins/QmcChem/target_pt2_qmc.irp.f @@ -62,7 +62,7 @@ program e_curve endif enddo call compute_energy(psi_bilinear_matrix_values_save,E,m,norm) - print '(E9.1,2X,I8,2X,F10.2,2X,F10.8,2X,F12.6)', norm_sort(n), m, & + print '(E9.1,2X,I8,2X,F10.2,2X,F10.6,2X,F12.6)', norm_sort(n), m, & dble( elec_alpha_num**3 + elec_alpha_num**2 * m ) / & dble( elec_alpha_num**3 + elec_alpha_num**2 * n ), norm, E if (E < target_energy) then @@ -93,7 +93,7 @@ subroutine compute_energy(psi_bilinear_matrix_values_save, E, m, norm) m = 0 !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(k,l,det_i,det_j,ci,cj,hij) REDUCTION(+:norm,m,num) allocate( det_i(N_int,2), det_j(N_int,2)) - !$OMP DO + !$OMP DO schedule(guided) do k=1,n_det if (psi_bilinear_matrix_values_save(k) == 0.d0) then cycle diff --git a/scripts/compilation/qp_create_ninja.py b/scripts/compilation/qp_create_ninja.py index 094ec7a6..b0fe2f8c 100755 --- a/scripts/compilation/qp_create_ninja.py +++ b/scripts/compilation/qp_create_ninja.py @@ -35,7 +35,8 @@ except ImportError: from qp_path import QP_ROOT, QP_SRC, QP_EZFIO -EZFIO_LIB = join(QP_ROOT, "lib", "libezfio.a") +LIB = "" # join(QP_ROOT, "lib", "rdtsc.o") +EZFIO_LIB = join(QP_ROOT, "lib", "libezfio.a") ROOT_BUILD_NINJA = join(QP_ROOT, "config", "build.ninja") header = r"""# @@ -94,7 +95,7 @@ def ninja_create_env_variable(pwd_config_file): l_string.append(str_) lib_lapack = get_compilation_option(pwd_config_file, "LAPACK_LIB") - l_string.append("{0} = {1} {2}".format("LIB", lib_lapack, EZFIO_LIB)) + l_string.append("LIB = {0} {1} {2}".format(LIB, lib_lapack, EZFIO_LIB)) l_string.append("") diff --git a/scripts/get_basis.sh b/scripts/get_basis.sh index 1e969d5c..4d1fc61d 100755 --- a/scripts/get_basis.sh +++ b/scripts/get_basis.sh @@ -9,7 +9,6 @@ # - if [[ -z ${QP_ROOT} ]] then print "The QP_ROOT environment variable is not set." diff --git a/scripts/pseudo/elts_num_ele.py b/scripts/pseudo/elts_num_ele.py index 3c4ad09f..8f31f4f7 100644 --- a/scripts/pseudo/elts_num_ele.py +++ b/scripts/pseudo/elts_num_ele.py @@ -1,4 +1,5 @@ -name_to_elec = {"H": 1, +name_to_elec = {"X": 0, + "H": 1, "He": 2, "Li": 3, "Be": 4, diff --git a/scripts/pseudo/put_pseudo_in_ezfio.py b/scripts/pseudo/put_pseudo_in_ezfio.py index f44fb097..6ad69f10 100755 --- a/scripts/pseudo/put_pseudo_in_ezfio.py +++ b/scripts/pseudo/put_pseudo_in_ezfio.py @@ -58,17 +58,35 @@ def get_pseudo_str(l_atom): str_ = "" for a in l_atom: - l_cmd_atom = ["--atom", a] - l_cmd_head = [EMSL_path, "get_basis_data", - "--db_path", db_path, - "--basis", "BFD-Pseudo"] + if a is not 'X': + l_cmd_atom = ["--atom", a] - process = Popen(l_cmd_head + l_cmd_atom, stdout=PIPE, stderr=PIPE) + l_cmd_head = [EMSL_path, "get_basis_data", + "--db_path", db_path, + "--basis", "BFD-Pseudo"] - stdout, _ = process.communicate() - str_ += stdout.strip() + "\n" + process = Popen(l_cmd_head + l_cmd_atom, stdout=PIPE, stderr=PIPE) + stdout, _ = process.communicate() + str_ += stdout.strip() + "\n" + + else: # Dummy atoms + str_ += """Element Symbol: X +Number of replaced protons: 0 +Number of projectors: 0 + +Pseudopotential data: + +Local component: +Coeff. r^n Exp. +0.0 -1 0. +0.0 1 0. +0.0 0 0. + +Non-local component: +Coeff. r^n Exp. Proj. +""" return str_ diff --git a/scripts/utility/is_master_repository.py b/scripts/utility/is_master_repository.py index 9ead14a2..da5fb56f 100755 --- a/scripts/utility/is_master_repository.py +++ b/scripts/utility/is_master_repository.py @@ -1,7 +1,7 @@ #!/usr/bin/env python import subprocess -pipe = subprocess.Popen("git config --local --get remote.origin.url", \ +pipe = subprocess.Popen("git config --get remote.origin.url", \ shell=True, stdout=subprocess.PIPE) result = pipe.stdout.read() is_master_repository = "LCPQ/quantum_package" in result diff --git a/src/Determinants/EZFIO.cfg b/src/Determinants/EZFIO.cfg index 9fe2f23c..9613c6c1 100644 --- a/src/Determinants/EZFIO.cfg +++ b/src/Determinants/EZFIO.cfg @@ -40,6 +40,12 @@ doc: Force the wave function to be an eigenfunction of S^2 interface: ezfio,provider,ocaml default: False +[threshold_davidson] +type: Threshold +doc: Thresholds of Davidson's algorithm +interface: ezfio,provider,ocaml +default: 1.e-8 + [threshold_generators] type: Threshold doc: Thresholds on generators (fraction of the norm) diff --git a/src/Determinants/H_apply.irp.f b/src/Determinants/H_apply.irp.f index 6e2b3a5a..7e9861fe 100644 --- a/src/Determinants/H_apply.irp.f +++ b/src/Determinants/H_apply.irp.f @@ -2,11 +2,11 @@ use bitmasks use omp_lib type H_apply_buffer_type -integer :: N_det -integer :: sze -integer(bit_kind), pointer :: det(:,:,:) -double precision , pointer :: coef(:,:) -double precision , pointer :: e2(:,:) + integer :: N_det + integer :: sze + integer(bit_kind), pointer :: det(:,:,:) + double precision , pointer :: coef(:,:) + double precision , pointer :: e2(:,:) end type H_apply_buffer_type type(H_apply_buffer_type), pointer :: H_apply_buffer(:) @@ -41,6 +41,15 @@ type(H_apply_buffer_type), pointer :: H_apply_buffer(:) call omp_init_lock(H_apply_buffer_lock(1,iproc)) !$OMP END PARALLEL endif + do iproc=2,nproc-1 + if (.not.associated(H_apply_buffer(iproc)%det)) then + print *, ' ===================== Error =================== ' + print *, 'H_apply_buffer_allocated should be provided outside' + print *, 'of an OpenMP section' + print *, ' =============================================== ' + stop + endif + enddo END_PROVIDER @@ -111,7 +120,6 @@ subroutine copy_H_apply_buffer_to_wf double precision, allocatable :: buffer_coef(:,:) integer :: i,j,k integer :: N_det_old - integer :: iproc PROVIDE H_apply_buffer_allocated @@ -158,7 +166,7 @@ subroutine copy_H_apply_buffer_to_wf enddo !$OMP PARALLEL DEFAULT(SHARED) & !$OMP PRIVATE(j,k,i) FIRSTPRIVATE(N_det_old) & - !$OMP SHARED(N_int,H_apply_buffer,psi_det,psi_coef,N_states) + !$OMP SHARED(N_int,H_apply_buffer,psi_det,psi_coef,N_states,psi_det_size) j=0 !$ j=omp_get_thread_num() do k=0,j-1 diff --git a/src/Determinants/davidson.irp.f b/src/Determinants/davidson.irp.f index bceb145a..3432ab2e 100644 --- a/src/Determinants/davidson.irp.f +++ b/src/Determinants/davidson.irp.f @@ -90,54 +90,73 @@ end function subroutine tamiser(key, idx, no, n, Nint, N_key) use bitmasks - implicit none - integer(bit_kind),intent(inout) :: key(Nint, 2, N_key) + + BEGIN_DOC +! Uncodumented : TODO + END_DOC integer,intent(in) :: no, n, Nint, N_key + integer(bit_kind),intent(inout) :: key(Nint, 2, N_key) integer,intent(inout) :: idx(N_key) integer :: k,j,tmpidx integer(bit_kind) :: tmp(Nint, 2) logical :: det_inf + integer :: ni k = no j = 2*k do while(j <= n) - if(j < n .and. det_inf(key(:,:,j), key(:,:,j+1), Nint)) then - j = j+1 - end if - if(det_inf(key(:,:,k), key(:,:,j), Nint)) then - tmp(:,:) = key(:,:,k) - key(:,:,k) = key(:,:,j) - key(:,:,j) = tmp(:,:) + if(j < n) then + if (det_inf(key(1,1,j), key(1,1,j+1), Nint)) then + j = j+1 + endif + endif + if(det_inf(key(1,1,k), key(1,1,j), Nint)) then + do ni=1,Nint + tmp(ni,1) = key(ni,1,k) + tmp(ni,2) = key(ni,2,k) + key(ni,1,k) = key(ni,1,j) + key(ni,2,k) = key(ni,2,j) + key(ni,1,j) = tmp(ni,1) + key(ni,2,j) = tmp(ni,2) + enddo tmpidx = idx(k) idx(k) = idx(j) idx(j) = tmpidx k = j - j = 2*k + j = k+k else return - end if - end do + endif + enddo end subroutine subroutine sort_dets_ba_v(key_in, key_out, idx, shortcut, version, N_key, Nint) use bitmasks implicit none - integer(bit_kind),intent(in) :: key_in(Nint,2,N_key) - integer(bit_kind) :: key(Nint,2,N_key) - integer(bit_kind),intent(out) :: key_out(Nint,N_key) - integer,intent(out) :: idx(N_key) - integer,intent(out) :: shortcut(0:N_key+1) - integer(bit_kind),intent(out) :: version(Nint,N_key+1) - integer, intent(in) :: Nint, N_key - integer(bit_kind) :: tmp(Nint, 2,N_key) + BEGIN_DOC +! Uncodumented : TODO + END_DOC + integer, intent(in) :: Nint, N_key + integer(bit_kind),intent(in) :: key_in(Nint,2,N_key) + integer(bit_kind),intent(out) :: key_out(Nint,N_key) + integer,intent(out) :: idx(N_key) + integer,intent(out) :: shortcut(0:N_key+1) + integer(bit_kind),intent(out) :: version(Nint,N_key+1) + integer(bit_kind), allocatable :: key(:,:,:) + integer :: i,ni - key(:,1,:N_key) = key_in(:,2,:N_key) - key(:,2,:N_key) = key_in(:,1,:N_key) + allocate ( key(Nint,2,N_key) ) + do i=1,N_key + do ni=1,Nint + key(ni,1,i) = key_in(ni,2,i) + key(ni,2,i) = key_in(ni,1,i) + enddo + enddo - call sort_dets_ab_v(key, key_out, idx, shortcut, version, N_key, Nint) + deallocate ( key ) end subroutine @@ -146,18 +165,25 @@ subroutine sort_dets_ab_v(key_in, key_out, idx, shortcut, version, N_key, Nint) use bitmasks implicit none + BEGIN_DOC +! Uncodumented : TODO + END_DOC + integer, intent(in) :: Nint, N_key integer(bit_kind),intent(in) :: key_in(Nint,2,N_key) - integer(bit_kind) :: key(Nint,2,N_key) integer(bit_kind),intent(out) :: key_out(Nint,N_key) integer,intent(out) :: idx(N_key) integer,intent(out) :: shortcut(0:N_key+1) integer(bit_kind),intent(out) :: version(Nint,N_key+1) - integer, intent(in) :: Nint, N_key + integer(bit_kind), allocatable :: key(:,:,:) integer(bit_kind) :: tmp(Nint, 2) integer :: tmpidx,i,ni - key(:,:,:) = key_in(:,:,:) + allocate (key(Nint,2,N_key)) do i=1,N_key + do ni=1,Nint + key(ni,1,i) = key_in(ni,1,i) + key(ni,2,i) = key_in(ni,2,i) + enddo idx(i) = i end do @@ -166,9 +192,14 @@ subroutine sort_dets_ab_v(key_in, key_out, idx, shortcut, version, N_key, Nint) end do do i=N_key,2,-1 - tmp(:,:) = key(:,:,i) - key(:,:,i) = key(:,:,1) - key(:,:,1) = tmp(:,:) + do ni=1,Nint + tmp(ni,1) = key(ni,1,i) + tmp(ni,2) = key(ni,2,i) + key(ni,1,i) = key(ni,1,1) + key(ni,2,i) = key(ni,2,1) + key(ni,1,1) = tmp(ni,1) + key(ni,2,1) = tmp(ni,2) + enddo tmpidx = idx(i) idx(i) = idx(1) idx(1) = tmpidx @@ -177,7 +208,9 @@ subroutine sort_dets_ab_v(key_in, key_out, idx, shortcut, version, N_key, Nint) shortcut(0) = 1 shortcut(1) = 1 - version(:,1) = key(:,1,1) + do ni=1,Nint + version(ni,1) = key(ni,1,1) + enddo do i=2,N_key do ni=1,nint if(key(ni,1,i) /= key(ni,1,i-1)) then @@ -189,15 +222,23 @@ subroutine sort_dets_ab_v(key_in, key_out, idx, shortcut, version, N_key, Nint) end do end do shortcut(shortcut(0)+1) = N_key+1 - key_out(:,:) = key(:,2,:) + do i=1,N_key + do ni=1,Nint + key_out(ni,i) = key(ni,2,i) + enddo + enddo + deallocate (key) end subroutine -c subroutine sort_dets_ab(key, idx, shortcut, N_key, Nint) use bitmasks implicit none + + BEGIN_DOC +! Uncodumented : TODO + END_DOC integer(bit_kind),intent(inout) :: key(Nint,2,N_key) integer,intent(out) :: idx(N_key) integer,intent(out) :: shortcut(0:N_key+1) @@ -214,9 +255,15 @@ subroutine sort_dets_ab(key, idx, shortcut, N_key, Nint) end do do i=N_key,2,-1 - tmp(:,:) = key(:,:,i) - key(:,:,i) = key(:,:,1) - key(:,:,1) = tmp(:,:) + do ni=1,Nint + tmp(ni,1) = key(ni,1,i) + tmp(ni,2) = key(ni,2,i) + key(ni,1,i) = key(ni,1,1) + key(ni,2,i) = key(ni,2,1) + key(ni,1,1) = tmp(ni,1) + key(ni,2,1) = tmp(ni,2) + enddo + tmpidx = idx(i) idx(i) = idx(1) idx(1) = tmpidx @@ -546,14 +593,12 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun abort_here = abort_all end - BEGIN_PROVIDER [ character(64), davidson_criterion ] -&BEGIN_PROVIDER [ double precision, davidson_threshold ] +BEGIN_PROVIDER [ character(64), davidson_criterion ] implicit none BEGIN_DOC ! Can be : [ energy | residual | both | wall_time | cpu_time | iterations ] END_DOC davidson_criterion = 'residual' - davidson_threshold = 1.d-10 END_PROVIDER subroutine davidson_converged(energy,residual,wall,iterations,cpu,N_st,converged) @@ -576,20 +621,20 @@ subroutine davidson_converged(energy,residual,wall,iterations,cpu,N_st,converged E = energy - energy_old energy_old = energy if (davidson_criterion == 'energy') then - converged = dabs(maxval(E(1:N_st))) < davidson_threshold + converged = dabs(maxval(E(1:N_st))) < threshold_davidson else if (davidson_criterion == 'residual') then - converged = dabs(maxval(residual(1:N_st))) < davidson_threshold + converged = dabs(maxval(residual(1:N_st))) < threshold_davidson else if (davidson_criterion == 'both') then converged = dabs(maxval(residual(1:N_st))) + dabs(maxval(E(1:N_st)) ) & - < davidson_threshold + < threshold_davidson else if (davidson_criterion == 'wall_time') then call wall_time(time) - converged = time - wall > davidson_threshold + converged = time - wall > threshold_davidson else if (davidson_criterion == 'cpu_time') then call cpu_time(time) - converged = time - cpu > davidson_threshold + converged = time - cpu > threshold_davidson else if (davidson_criterion == 'iterations') then - converged = iterations >= int(davidson_threshold) + converged = iterations >= int(threshold_davidson) endif converged = converged.or.abort_here end diff --git a/src/Determinants/determinants.irp.f b/src/Determinants/determinants.irp.f index d1c36163..5fe18c49 100644 --- a/src/Determinants/determinants.irp.f +++ b/src/Determinants/determinants.irp.f @@ -8,6 +8,7 @@ BEGIN_PROVIDER [ integer, N_det ] logical :: exists character*64 :: label PROVIDE ezfio_filename + PROVIDE nproc if (read_wf) then call ezfio_has_determinants_n_det(exists) if (exists) then diff --git a/src/Determinants/guess_lowest_state.irp.f b/src/Determinants/guess_lowest_state.irp.f new file mode 100644 index 00000000..f6d0a004 --- /dev/null +++ b/src/Determinants/guess_lowest_state.irp.f @@ -0,0 +1,162 @@ +program first_guess + use bitmasks + implicit none + BEGIN_DOC + ! Select all the determinants with the lowest energy as a starting point. + END_DOC + integer :: i,j + double precision, allocatable :: orb_energy(:) + double precision :: E + integer, allocatable :: kept(:) + integer :: nelec_kept(2) + character :: occ_char, keep_char + + PROVIDE H_apply_buffer_allocated psi_det + allocate (orb_energy(mo_tot_num), kept(0:mo_tot_num)) + nelec_kept(1:2) = 0 + kept(0) = 0 + + print *, 'Orbital energies' + print *, '================' + print *, '' + do i=1,mo_tot_num + keep_char = ' ' + occ_char = '-' + orb_energy(i) = mo_mono_elec_integral(i,i) + do j=1,elec_beta_num + if (i==j) cycle + orb_energy(i) += mo_bielec_integral_jj_anti(i,j) + enddo + do j=1,elec_alpha_num + orb_energy(i) += mo_bielec_integral_jj(i,j) + enddo + if ( (orb_energy(i) > -.5d0).and.(orb_energy(i) < .1d0) ) then + kept(0) += 1 + keep_char = 'X' + kept( kept(0) ) = i + if (i <= elec_beta_num) then + nelec_kept(2) += 1 + endif + if (i <= elec_alpha_num) then + nelec_kept(1) += 1 + endif + endif + if (i <= elec_alpha_num) then + if (i <= elec_beta_num) then + occ_char = '#' + else + occ_char = '+' + endif + endif + print '(I4, 3X, A, 3X, F10.6, 3X, A)', i, occ_char, orb_energy(i), keep_char + enddo + + + integer, allocatable :: list (:,:) + integer(bit_kind), allocatable :: string(:,:) + allocate ( list(N_int*bit_kind_size,2), string(N_int,2) ) + + string = ref_bitmask + call bitstring_to_list( string(1,1), list(1,1), elec_alpha_num, N_int) + call bitstring_to_list( string(1,2), list(1,2), elec_beta_num , N_int) + + psi_det_alpha_unique(:,1) = string(:,1) + psi_det_beta_unique (:,1) = string(:,2) + N_det_alpha_unique = 1 + N_det_beta_unique = 1 + + integer :: i1,i2,i3,i4,i5,i6,i7,i8,i9 + + psi_det_size = kept(0)**(nelec_kept(1)+nelec_kept(2)) + print *, kept(0), nelec_kept(:) + call write_int(6,psi_det_size,'psi_det_size') + TOUCH psi_det_size + +BEGIN_SHELL [ /usr/bin/python ] + +template_alpha_ext = """ +do %(i2)s = %(i1)s-1,1,-1 + list(elec_alpha_num-%(i)d,1) = kept(%(i2)s) + call list_to_bitstring( string(1,1), list(1,1), elec_alpha_num, N_int) +""" + +template_alpha = """ +do %(i2)s = %(i1)s-1,1,-1 + list(elec_alpha_num-%(i)d,1) = kept(%(i2)s) + call list_to_bitstring( string(1,1), list(1,1), elec_alpha_num, N_int) + N_det_alpha_unique += 1 + psi_det_alpha_unique(:,N_det_alpha_unique) = string(:,1) +""" + +template_beta_ext = """ +do %(i2)s = %(i1)s-1,1,-1 + list(elec_beta_num-%(i)d,2) = kept(%(i2)s) + call list_to_bitstring( string(1,2), list(1,2), elec_beta_num, N_int) +""" +template_beta = """ +do %(i2)s = %(i1)s-1,1,-1 + list(elec_beta_num-%(i)d,2) = kept(%(i2)s) + call list_to_bitstring( string(1,2), list(1,2), elec_beta_num, N_int) + N_det_beta_unique += 1 + psi_det_beta_unique(:,N_det_beta_unique) = string(:,2) +""" + +def write(template_ext,template,imax): + print "case(%d)"%(imax) + def aux(i2,i1,i,j): + if (i==imax-1): + print template%locals() + else: + print template_ext%locals() + i += 1 + j -= 1 + if (i != imax): + i1 = "i%d"%(i) + i2 = "i%d"%(i+1) + aux(i2,i1,i,j) + print "enddo" + + i2 = "i1" + i1 = "kept(0)+1" + i = 0 + aux(i2,i1,i,imax) + +def main(): + print """ + select case (nelec_kept(1)) + case(0) + continue + """ + for imax in range(1,10): + write(template_alpha_ext,template_alpha,imax) + + print """ + end select + + select case (nelec_kept(2)) + case(0) + continue + """ + for imax in range(1,10): + write(template_beta_ext,template_beta,imax) + print "end select" + +main() + +END_SHELL + + TOUCH N_det_alpha_unique N_det_beta_unique psi_det_alpha_unique psi_det_beta_unique + call create_wf_of_psi_bilinear_matrix(.False.) + call diagonalize_ci + j= N_det + do i=1,N_det + if (psi_average_norm_contrib_sorted(i) < 1.d-6) then + j = i-1 + exit + endif +! call debug_det(psi_det_sorted(1,1,i),N_int) + enddo + call save_wavefunction_general(j,N_states,psi_det_sorted,size(psi_coef_sorted,1),psi_coef_sorted) + + deallocate(orb_energy, kept, list, string) +end diff --git a/src/Determinants/program_initial_determinants.irp.f b/src/Determinants/program_initial_determinants.irp.f deleted file mode 100644 index 6375af22..00000000 --- a/src/Determinants/program_initial_determinants.irp.f +++ /dev/null @@ -1,138 +0,0 @@ -program pouet - implicit none - print*,'HF energy = ',ref_bitmask_energy + nuclear_repulsion - call routine - -end -subroutine routine - use bitmasks - implicit none - integer :: i,j,k,l - double precision :: hij,get_mo_bielec_integral - double precision :: hmono,h_bi_ispin,h_bi_other_spin - integer(bit_kind),allocatable :: key_tmp(:,:) - integer, allocatable :: occ(:,:) - integer :: n_occ_alpha, n_occ_beta - ! First checks - print*,'N_int = ',N_int - print*,'mo_tot_num = ',mo_tot_num - print*,'mo_tot_num / 64+1= ',mo_tot_num/64+1 - ! We print the HF determinant - do i = 1, N_int - print*,'ref_bitmask(i,1) = ',ref_bitmask(i,1) - print*,'ref_bitmask(i,2) = ',ref_bitmask(i,2) - enddo - print*,'' - print*,'Hartree Fock determinant ...' - call debug_det(ref_bitmask,N_int) - allocate(key_tmp(N_int,2)) - ! We initialize key_tmp to the Hartree Fock one - key_tmp = ref_bitmask - integer :: i_hole,i_particle,ispin,i_ok,other_spin - ! We do a mono excitation on the top of the HF determinant - write(*,*)'Enter the (hole, particle) couple for the mono excitation ...' - read(5,*)i_hole,i_particle -!!i_hole = 4 -!!i_particle = 20 - write(*,*)'Enter the ispin variable ...' - write(*,*)'ispin = 1 ==> alpha ' - write(*,*)'ispin = 2 ==> beta ' - read(5,*)ispin - if(ispin == 1)then - other_spin = 2 - else if(ispin == 2)then - other_spin = 1 - else - print*,'PB !! ' - print*,'ispin must be 1 or 2 !' - stop - endif -!!ispin = 1 - call do_mono_excitation(key_tmp,i_hole,i_particle,ispin,i_ok) - ! We check if it the excitation was possible with "i_ok" - if(i_ok == -1)then - print*,'i_ok = ',i_ok - print*,'You can not do this excitation because of Pauli principle ...' - print*,'check your hole particle couple, there must be something wrong ...' - stop - - endif - print*,'New det = ' - call debug_det(key_tmp,N_int) - call i_H_j(key_tmp,ref_bitmask,N_int,hij) - ! We calculate the H matrix element between the new determinant and HF - print*,' = ',hij - print*,'' - print*,'' - print*,'Recalculating it old school style ....' - print*,'' - print*,'' - ! We recalculate this old school style !!! - ! Mono electronic part - hmono = mo_mono_elec_integral(i_hole,i_particle) - print*,'' - print*,'Mono electronic part ' - print*,'' - print*,' = ',hmono - h_bi_ispin = 0.d0 - h_bi_other_spin = 0.d0 - print*,'' - print*,'Getting all the info for the calculation of the bi electronic part ...' - print*,'' - allocate (occ(N_int*bit_kind_size,2)) - ! We get the occupation of the alpha electrons in occ(:,1) - call bitstring_to_list(key_tmp(1,1), occ(1,1), n_occ_alpha, N_int) - print*,'n_occ_alpha = ',n_occ_alpha - print*,'elec_alpha_num = ',elec_alpha_num - ! We get the occupation of the beta electrons in occ(:,2) - call bitstring_to_list(key_tmp(1,2), occ(1,2), n_occ_beta, N_int) - print*,'n_occ_beta = ',n_occ_beta - print*,'elec_beta_num = ',elec_beta_num - ! We print the occupation of the alpha electrons - print*,'Alpha electrons !' - do i = 1, n_occ_alpha - print*,'i = ',i - print*,'occ(i,1) = ',occ(i,1) - enddo - ! We print the occupation of the beta electrons - print*,'Alpha electrons !' - do i = 1, n_occ_beta - print*,'i = ',i - print*,'occ(i,2) = ',occ(i,2) - enddo - integer :: exc(0:2,2,2),degree,h1,p1,h2,p2,s1,s2 - double precision :: phase - - call get_excitation_degree(key_tmp,ref_bitmask,degree,N_int) - print*,'degree = ',degree - call get_mono_excitation(ref_bitmask,key_tmp,exc,phase,N_int) - call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) - print*,'h1 = ',h1 - print*,'p1 = ',p1 - print*,'s1 = ',s1 - print*,'phase = ',phase - do i = 1, elec_num_tab(ispin) - integer :: orb_occupied - orb_occupied = occ(i,ispin) - h_bi_ispin += get_mo_bielec_integral(i_hole,orb_occupied,i_particle,orb_occupied,mo_integrals_map) & - -get_mo_bielec_integral(i_hole,i_particle,orb_occupied,orb_occupied,mo_integrals_map) - enddo - print*,'h_bi_ispin = ',h_bi_ispin - - do i = 1, elec_num_tab(other_spin) - orb_occupied = occ(i,other_spin) - h_bi_other_spin += get_mo_bielec_integral(i_hole,orb_occupied,i_particle,orb_occupied,mo_integrals_map) - enddo - print*,'h_bi_other_spin = ',h_bi_other_spin - print*,'h_bi_ispin + h_bi_other_spin = ',h_bi_ispin + h_bi_other_spin - - print*,'Total matrix element = ',phase*(h_bi_ispin + h_bi_other_spin + hmono) -!i = 1 -!j = 1 -!k = 1 -!l = 1 -!hij = get_mo_bielec_integral(i,j,k,l,mo_integrals_map) -!print*,' = ',hij - - -end diff --git a/src/Determinants/s2.irp.f b/src/Determinants/s2.irp.f index e836d25d..6da7b8ec 100644 --- a/src/Determinants/s2.irp.f +++ b/src/Determinants/s2.irp.f @@ -109,81 +109,81 @@ end subroutine get_s2_u0(psi_keys_tmp,psi_coefs_tmp,n,nmax,s2) implicit none use bitmasks - integer(bit_kind), intent(in) :: psi_keys_tmp(N_int,2,nmax) - integer, intent(in) :: n,nmax - double precision, intent(in) :: psi_coefs_tmp(nmax) - double precision, intent(out) :: s2 - double precision :: s2_tmp - integer :: i,j,l,jj,ii + integer(bit_kind), intent(in) :: psi_keys_tmp(N_int,2,nmax) + integer, intent(in) :: n,nmax + double precision, intent(in) :: psi_coefs_tmp(nmax) + double precision, intent(out) :: s2 + double precision :: s2_tmp + integer :: i,j,l,jj,ii integer, allocatable :: idx(:) - - integer :: shortcut(0:n+1), sort_idx(n) - integer(bit_kind) :: sorted(N_int,n), version(N_int,n) + + integer, allocatable :: shortcut(:), sort_idx(:) + integer(bit_kind), allocatable :: sorted(:,:), version(:,:) integer :: sh, sh2, ni, exa, ext, org_i, org_j, endi, pass double precision :: davidson_threshold_bis - - !PROVIDE davidson_threshold + allocate (shortcut(0:n+1), sort_idx(n), sorted(N_int,n), version(N_int,n)) s2 = 0.d0 - davidson_threshold_bis = davidson_threshold - !$OMP PARALLEL DEFAULT(NONE) & - !$OMP PRIVATE(i,j,s2_tmp,idx,sh, sh2, ni, exa, ext, org_i, org_j, endi, pass) & - !$OMP SHARED(n,psi_coefs_tmp,psi_keys_tmp,N_int,davidson_threshold,shortcut,sorted,sort_idx,version)& - !$OMP REDUCTION(+:s2) - allocate(idx(0:n)) - - - !$OMP SINGLE + davidson_threshold_bis = threshold_davidson call sort_dets_ab_v(psi_keys_tmp, sorted, sort_idx, shortcut, version, n, N_int) - !$OMP END SINGLE + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(i,j,s2_tmp,sh, sh2, ni, exa, ext, org_i, org_j, endi, pass)& + !$OMP SHARED(n,psi_coefs_tmp,psi_keys_tmp,N_int,threshold_davidson,shortcut,sorted,sort_idx,version)& + !$OMP REDUCTION(+:s2) !$OMP DO SCHEDULE(dynamic) do sh=1,shortcut(0) - - do sh2=1,sh - exa = 0 - do ni=1,N_int - exa += popcnt(xor(version(ni,sh), version(ni,sh2))) - end do - if(exa > 2) then - cycle - end if - do i=shortcut(sh),shortcut(sh+1)-1 - if(sh==sh2) then - endi = i-1 - else - endi = shortcut(sh2+1)-1 + do sh2=1,sh + exa = 0 + do ni=1,N_int + exa += popcnt(xor(version(ni,sh), version(ni,sh2))) + end do + if(exa > 2) then + cycle end if - do j=shortcut(sh2),endi - ext = exa - do ni=1,N_int - ext += popcnt(xor(sorted(ni,i), sorted(ni,j))) - end do - if(ext <= 4) then - org_i = sort_idx(i) - org_j = sort_idx(j) - - if ( dabs(psi_coefs_tmp(org_j)) + dabs(psi_coefs_tmp(org_i)) & - > davidson_threshold ) then - call get_s2(psi_keys_tmp(1,1,org_i),psi_keys_tmp(1,1,org_j),s2_tmp,N_int) - s2 = s2 + psi_coefs_tmp(org_i)*psi_coefs_tmp(org_j)*s2_tmp - endif + do i=shortcut(sh),shortcut(sh+1)-1 + if(sh==sh2) then + endi = i-1 + else + endi = shortcut(sh2+1)-1 end if + + do j=shortcut(sh2),endi + ext = exa + do ni=1,N_int + ext += popcnt(xor(sorted(ni,i), sorted(ni,j))) + end do + if(ext <= 4) then + org_i = sort_idx(i) + org_j = sort_idx(j) + + if ( dabs(psi_coefs_tmp(org_j)) + dabs(psi_coefs_tmp(org_i))& + > threshold_davidson ) then + call get_s2(psi_keys_tmp(1,1,org_i),psi_keys_tmp(1,1,org_j),s2_tmp,N_int) + s2 = s2 + psi_coefs_tmp(org_i)*psi_coefs_tmp(org_j)*s2_tmp + endif + end if + end do end do end do - end do enddo !$OMP END DO - !$OMP SINGLE + !$OMP END PARALLEL + call sort_dets_ba_v(psi_keys_tmp, sorted, sort_idx, shortcut, version, n, N_int) - !$OMP END SINGLE + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(i,j,s2_tmp,sh, sh2, ni, exa, ext, org_i, org_j, endi, pass)& + !$OMP SHARED(n,psi_coefs_tmp,psi_keys_tmp,N_int,threshold_davidson,shortcut,sorted,sort_idx,version)& + !$OMP REDUCTION(+:s2) !$OMP DO SCHEDULE(dynamic) do sh=1,shortcut(0) - do i=shortcut(sh),shortcut(sh+1)-1 + do i=shortcut(sh),shortcut(sh+1)-1 do j=shortcut(sh),i-1 ext = 0 do ni=1,N_int @@ -193,25 +193,25 @@ subroutine get_s2_u0(psi_keys_tmp,psi_coefs_tmp,n,nmax,s2) org_i = sort_idx(i) org_j = sort_idx(j) - if ( dabs(psi_coefs_tmp(org_j)) + dabs(psi_coefs_tmp(org_i)) & - > davidson_threshold ) then + if ( dabs(psi_coefs_tmp(org_j)) + dabs(psi_coefs_tmp(org_i))& + > threshold_davidson ) then call get_s2(psi_keys_tmp(1,1,org_i),psi_keys_tmp(1,1,org_j),s2_tmp,N_int) s2 = s2 + psi_coefs_tmp(org_i)*psi_coefs_tmp(org_j)*s2_tmp endif end if end do - end do + end do enddo !$OMP END DO - deallocate(idx) - !$OMP END PARALLEL - s2 = s2+s2 - do i=1,n - call get_s2(psi_keys_tmp(1,1,i),psi_keys_tmp(1,1,i),s2_tmp,N_int) - s2 = s2 + psi_coefs_tmp(i)*psi_coefs_tmp(i)*s2_tmp - enddo - s2 = s2 + S_z2_Sz + !$OMP END PARALLEL + s2 = s2+s2 + do i=1,n + call get_s2(psi_keys_tmp(1,1,i),psi_keys_tmp(1,1,i),s2_tmp,N_int) + s2 = s2 + psi_coefs_tmp(i)*psi_coefs_tmp(i)*s2_tmp + enddo + s2 = s2 + S_z2_Sz + deallocate (shortcut, sort_idx, sorted, version) end diff --git a/src/Determinants/slater_rules.irp.f b/src/Determinants/slater_rules.irp.f index ea3e9f3c..b60bf0f7 100644 --- a/src/Determinants/slater_rules.irp.f +++ b/src/Determinants/slater_rules.irp.f @@ -15,7 +15,7 @@ subroutine get_excitation_degree(key1,key2,degree,Nint) degree = popcnt(xor( key1(1,1), key2(1,1))) + & popcnt(xor( key1(1,2), key2(1,2))) - !DEC$ NOUNROLL + !DIR$ NOUNROLL do l=2,Nint degree = degree+ popcnt(xor( key1(l,1), key2(l,1))) + & popcnt(xor( key1(l,2), key2(l,2))) @@ -365,7 +365,7 @@ subroutine i_H_j(key_i,key_j,Nint,hij) integer :: exc(0:2,2,2) integer :: degree - double precision :: get_mo_bielec_integral + double precision :: get_mo_bielec_integral_schwartz integer :: m,n,p,q integer :: i,j,k integer :: occ(Nint*bit_kind_size,2) @@ -383,38 +383,38 @@ subroutine i_H_j(key_i,key_j,Nint,hij) ASSERT (sum(popcnt(key_j(:,2))) == elec_beta_num) hij = 0.d0 - !DEC$ FORCEINLINE + !DIR$ FORCEINLINE call get_excitation_degree(key_i,key_j,degree,Nint) select case (degree) case (2) call get_double_excitation(key_i,key_j,exc,phase,Nint) if (exc(0,1,1) == 1) then ! Mono alpha, mono beta - hij = phase*get_mo_bielec_integral( & + hij = phase*get_mo_bielec_integral_schwartz( & exc(1,1,1), & exc(1,1,2), & exc(1,2,1), & exc(1,2,2) ,mo_integrals_map) else if (exc(0,1,1) == 2) then ! Double alpha - hij = phase*(get_mo_bielec_integral( & + hij = phase*(get_mo_bielec_integral_schwartz( & exc(1,1,1), & exc(2,1,1), & exc(1,2,1), & exc(2,2,1) ,mo_integrals_map) - & - get_mo_bielec_integral( & + get_mo_bielec_integral_schwartz( & exc(1,1,1), & exc(2,1,1), & exc(2,2,1), & exc(1,2,1) ,mo_integrals_map) ) else if (exc(0,1,2) == 2) then ! Double beta - hij = phase*(get_mo_bielec_integral( & + hij = phase*(get_mo_bielec_integral_schwartz( & exc(1,1,2), & exc(2,1,2), & exc(1,2,2), & exc(2,2,2) ,mo_integrals_map) - & - get_mo_bielec_integral( & + get_mo_bielec_integral_schwartz( & exc(1,1,2), & exc(2,1,2), & exc(2,2,2), & @@ -432,15 +432,15 @@ subroutine i_H_j(key_i,key_j,Nint,hij) do k = 1, elec_alpha_num i = occ(k,1) if (.not.has_mipi(i)) then - mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) - miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map) + mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) + miip(i) = get_mo_bielec_integral_schwartz(m,i,i,p,mo_integrals_map) has_mipi(i) = .True. endif enddo do k = 1, elec_beta_num i = occ(k,2) if (.not.has_mipi(i)) then - mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) + mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) has_mipi(i) = .True. endif enddo @@ -459,15 +459,15 @@ subroutine i_H_j(key_i,key_j,Nint,hij) do k = 1, elec_beta_num i = occ(k,2) if (.not.has_mipi(i)) then - mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) - miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map) + mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) + miip(i) = get_mo_bielec_integral_schwartz(m,i,i,p,mo_integrals_map) has_mipi(i) = .True. endif enddo do k = 1, elec_alpha_num i = occ(k,1) if (.not.has_mipi(i)) then - mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) + mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) has_mipi(i) = .True. endif enddo @@ -501,7 +501,7 @@ subroutine i_H_j_phase_out(key_i,key_j,Nint,hij,phase,exc,degree) integer,intent(out) :: exc(0:2,2,2) integer,intent(out) :: degree - double precision :: get_mo_bielec_integral + double precision :: get_mo_bielec_integral_schwartz integer :: m,n,p,q integer :: i,j,k integer :: occ(Nint*bit_kind_size,2) @@ -519,38 +519,38 @@ subroutine i_H_j_phase_out(key_i,key_j,Nint,hij,phase,exc,degree) ASSERT (sum(popcnt(key_j(:,2))) == elec_beta_num) hij = 0.d0 - !DEC$ FORCEINLINE + !DIR$ FORCEINLINE call get_excitation_degree(key_i,key_j,degree,Nint) select case (degree) case (2) call get_double_excitation(key_i,key_j,exc,phase,Nint) if (exc(0,1,1) == 1) then ! Mono alpha, mono beta - hij = phase*get_mo_bielec_integral( & + hij = phase*get_mo_bielec_integral_schwartz( & exc(1,1,1), & exc(1,1,2), & exc(1,2,1), & exc(1,2,2) ,mo_integrals_map) else if (exc(0,1,1) == 2) then ! Double alpha - hij = phase*(get_mo_bielec_integral( & + hij = phase*(get_mo_bielec_integral_schwartz( & exc(1,1,1), & exc(2,1,1), & exc(1,2,1), & exc(2,2,1) ,mo_integrals_map) - & - get_mo_bielec_integral( & + get_mo_bielec_integral_schwartz( & exc(1,1,1), & exc(2,1,1), & exc(2,2,1), & exc(1,2,1) ,mo_integrals_map) ) else if (exc(0,1,2) == 2) then ! Double beta - hij = phase*(get_mo_bielec_integral( & + hij = phase*(get_mo_bielec_integral_schwartz( & exc(1,1,2), & exc(2,1,2), & exc(1,2,2), & exc(2,2,2) ,mo_integrals_map) - & - get_mo_bielec_integral( & + get_mo_bielec_integral_schwartz( & exc(1,1,2), & exc(2,1,2), & exc(2,2,2), & @@ -568,15 +568,15 @@ subroutine i_H_j_phase_out(key_i,key_j,Nint,hij,phase,exc,degree) do k = 1, elec_alpha_num i = occ(k,1) if (.not.has_mipi(i)) then - mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) - miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map) + mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) + miip(i) = get_mo_bielec_integral_schwartz(m,i,i,p,mo_integrals_map) has_mipi(i) = .True. endif enddo do k = 1, elec_beta_num i = occ(k,2) if (.not.has_mipi(i)) then - mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) + mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) has_mipi(i) = .True. endif enddo @@ -595,15 +595,15 @@ subroutine i_H_j_phase_out(key_i,key_j,Nint,hij,phase,exc,degree) do k = 1, elec_beta_num i = occ(k,2) if (.not.has_mipi(i)) then - mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) - miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map) + mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) + miip(i) = get_mo_bielec_integral_schwartz(m,i,i,p,mo_integrals_map) has_mipi(i) = .True. endif enddo do k = 1, elec_alpha_num i = occ(k,1) if (.not.has_mipi(i)) then - mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) + mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) has_mipi(i) = .True. endif enddo @@ -637,7 +637,7 @@ subroutine i_H_j_verbose(key_i,key_j,Nint,hij,hmono,hdouble) integer :: exc(0:2,2,2) integer :: degree - double precision :: get_mo_bielec_integral + double precision :: get_mo_bielec_integral_schwartz integer :: m,n,p,q integer :: i,j,k integer :: occ(Nint*bit_kind_size,2) @@ -657,38 +657,38 @@ subroutine i_H_j_verbose(key_i,key_j,Nint,hij,hmono,hdouble) hij = 0.d0 hmono = 0.d0 hdouble = 0.d0 - !DEC$ FORCEINLINE + !DIR$ FORCEINLINE call get_excitation_degree(key_i,key_j,degree,Nint) select case (degree) case (2) call get_double_excitation(key_i,key_j,exc,phase,Nint) if (exc(0,1,1) == 1) then ! Mono alpha, mono beta - hij = phase*get_mo_bielec_integral( & + hij = phase*get_mo_bielec_integral_schwartz( & exc(1,1,1), & exc(1,1,2), & exc(1,2,1), & exc(1,2,2) ,mo_integrals_map) else if (exc(0,1,1) == 2) then ! Double alpha - hij = phase*(get_mo_bielec_integral( & + hij = phase*(get_mo_bielec_integral_schwartz( & exc(1,1,1), & exc(2,1,1), & exc(1,2,1), & exc(2,2,1) ,mo_integrals_map) - & - get_mo_bielec_integral( & + get_mo_bielec_integral_schwartz( & exc(1,1,1), & exc(2,1,1), & exc(2,2,1), & exc(1,2,1) ,mo_integrals_map) ) else if (exc(0,1,2) == 2) then ! Double beta - hij = phase*(get_mo_bielec_integral( & + hij = phase*(get_mo_bielec_integral_schwartz( & exc(1,1,2), & exc(2,1,2), & exc(1,2,2), & exc(2,2,2) ,mo_integrals_map) - & - get_mo_bielec_integral( & + get_mo_bielec_integral_schwartz( & exc(1,1,2), & exc(2,1,2), & exc(2,2,2), & @@ -706,15 +706,15 @@ subroutine i_H_j_verbose(key_i,key_j,Nint,hij,hmono,hdouble) do k = 1, elec_alpha_num i = occ(k,1) if (.not.has_mipi(i)) then - mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) - miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map) + mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) + miip(i) = get_mo_bielec_integral_schwartz(m,i,i,p,mo_integrals_map) has_mipi(i) = .True. endif enddo do k = 1, elec_beta_num i = occ(k,2) if (.not.has_mipi(i)) then - mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) + mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) has_mipi(i) = .True. endif enddo @@ -733,15 +733,15 @@ subroutine i_H_j_verbose(key_i,key_j,Nint,hij,hmono,hdouble) do k = 1, elec_beta_num i = occ(k,2) if (.not.has_mipi(i)) then - mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) - miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map) + mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) + miip(i) = get_mo_bielec_integral_schwartz(m,i,i,p,mo_integrals_map) has_mipi(i) = .True. endif enddo do k = 1, elec_alpha_num i = occ(k,1) if (.not.has_mipi(i)) then - mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) + mipi(i) = get_mo_bielec_integral_schwartz(m,i,p,i,mo_integrals_map) has_mipi(i) = .True. endif enddo @@ -863,9 +863,17 @@ subroutine create_minilist_find_previous(key_mask, fullList, miniList, N_fullLis end subroutine -subroutine i_H_psi_nominilist(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array) +subroutine i_H_psi(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array) use bitmasks implicit none + BEGIN_DOC +! Computes = \sum_J c_J . +! +! Uses filter_connected_i_H_psi0 to get all the |J> to which |i> +! is connected. +! The i_H_psi_minilist is much faster but requires to build the +! minilists + END_DOC integer, intent(in) :: Nint, Ndet,Ndet_max,Nstate integer(bit_kind), intent(in) :: keys(Nint,2,Ndet) integer(bit_kind), intent(in) :: key(Nint,2) @@ -877,9 +885,6 @@ subroutine i_H_psi_nominilist(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_ar integer :: exc(0:2,2,2) double precision :: hij integer :: idx(0:Ndet) - BEGIN_DOC - ! for the various Nstates - END_DOC ASSERT (Nint > 0) ASSERT (N_int == Nint) @@ -891,7 +896,7 @@ subroutine i_H_psi_nominilist(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_ar call filter_connected_i_H_psi0(keys,key,Nint,Ndet,idx) do ii=1,idx(0) i = idx(ii) - !DEC$ FORCEINLINE + !DIR$ FORCEINLINE call i_H_j(keys(1,1,i),key,Nint,hij) do j = 1, Nstate i_H_psi_array(j) = i_H_psi_array(j) + coef(i,j)*hij @@ -900,7 +905,7 @@ subroutine i_H_psi_nominilist(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_ar end -subroutine i_H_psi(key,keys,idx_key,N_minilist,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array) +subroutine i_H_psi_minilist(key,keys,idx_key,N_minilist,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array) use bitmasks implicit none integer, intent(in) :: Nint, Ndet,Ndet_max,Nstate,idx_key(Ndet), N_minilist @@ -915,7 +920,10 @@ subroutine i_H_psi(key,keys,idx_key,N_minilist,coef,Nint,Ndet,Ndet_max,Nstate,i_ double precision :: hij integer :: idx(0:Ndet) BEGIN_DOC - ! for the various Nstates +! Computes = \sum_J c_J . +! +! Uses filter_connected_i_H_psi0 to get all the |J> to which |i> +! is connected. The |J> are searched in short pre-computed lists. END_DOC ASSERT (Nint > 0) @@ -925,17 +933,11 @@ subroutine i_H_psi(key,keys,idx_key,N_minilist,coef,Nint,Ndet,Ndet_max,Nstate,i_ ASSERT (Ndet_max >= Ndet) i_H_psi_array = 0.d0 - !call filter_connected_i_H_psi0(keys,key,Nint,Ndet,idx) call filter_connected_i_H_psi0(keys,key,Nint,N_minilist,idx) do ii=1,idx(0) - !i = idx_key(idx(ii)) i_in_key = idx(ii) i_in_coef = idx_key(idx(ii)) - !DEC$ FORCEINLINE -! ! call i_H_j(keys(1,1,i),key,Nint,hij) -! ! do j = 1, Nstate -! ! i_H_psi_array(j) = i_H_psi_array(j) + coef(i,j)*hij -! ! enddo + !DIR$ FORCEINLINE call i_H_j(keys(1,1,i_in_key),key,Nint,hij) do j = 1, Nstate i_H_psi_array(j) = i_H_psi_array(j) + coef(i_in_coef,j)*hij @@ -973,7 +975,7 @@ subroutine i_H_psi_sec_ord(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array n_interact = 0 do ii=1,idx(0) i = idx(ii) - !DEC$ FORCEINLINE + !DIR$ FORCEINLINE call i_H_j(keys(1,1,i),key,Nint,hij) if(dabs(hij).ge.1.d-8)then if(i.ne.1)then @@ -1028,7 +1030,7 @@ subroutine i_H_psi_SC2(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array,idx call filter_connected_i_H_psi0_SC2(keys,key,Nint,Ndet,idx,idx_repeat) do ii=1,idx(0) i = idx(ii) - !DEC$ FORCEINLINE + !DIR$ FORCEINLINE call i_H_j(keys(1,1,i),key,Nint,hij) do j = 1, Nstate i_H_psi_array(j) = i_H_psi_array(j) + coef(i,j)*hij @@ -1077,7 +1079,7 @@ subroutine i_H_psi_SC2_verbose(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_a do ii=1,idx(0) print*,'--' i = idx(ii) - !DEC$ FORCEINLINE + !DIR$ FORCEINLINE call i_H_j(keys(1,1,i),key,Nint,hij) if (i==1)then print*,'i==1 !!' @@ -1167,7 +1169,7 @@ subroutine get_excitation_degree_vector(key1,key2,degree,Nint,sze,idx) !DIR$ LOOP COUNT (1000) do i=1,sze d = 0 - !DEC$ LOOP COUNT MIN(4) + !DIR$ LOOP COUNT MIN(4) do m=1,Nint d = d + popcnt(xor( key1(m,1,i), key2(m,1))) & + popcnt(xor( key1(m,2,i), key2(m,2))) @@ -1211,14 +1213,14 @@ double precision function diag_H_mat_elem(det_in,Nint) nexc(1) = 0 nexc(2) = 0 do i=1,Nint - hole(i,1) = xor(det_in(i,1),ref_bitmask(i,1)) - hole(i,2) = xor(det_in(i,2),ref_bitmask(i,2)) + hole(i,1) = xor(det_in(i,1),ref_bitmask(i,1)) + hole(i,2) = xor(det_in(i,2),ref_bitmask(i,2)) particle(i,1) = iand(hole(i,1),det_in(i,1)) particle(i,2) = iand(hole(i,2),det_in(i,2)) hole(i,1) = iand(hole(i,1),ref_bitmask(i,1)) hole(i,2) = iand(hole(i,2),ref_bitmask(i,2)) - nexc(1) += popcnt(hole(i,1)) - nexc(2) += popcnt(hole(i,2)) + nexc(1) = nexc(1) + popcnt(hole(i,1)) + nexc(2) = nexc(2) + popcnt(hole(i,2)) enddo diag_H_mat_elem = ref_bitmask_energy @@ -1380,83 +1382,102 @@ subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint) integer :: i,j,k,l, jj,ii integer :: i0, j0 - integer :: shortcut(0:n+1), sort_idx(n) - integer(bit_kind) :: sorted(Nint,n), version(Nint,n) - + integer, allocatable :: shortcut(:), sort_idx(:) + integer(bit_kind), allocatable :: sorted(:,:), version(:,:) + integer(bit_kind) :: sorted_i(Nint) integer :: sh, sh2, ni, exa, ext, org_i, org_j, endi -! + double precision :: local_threshold ASSERT (Nint > 0) ASSERT (Nint == N_int) ASSERT (n>0) - PROVIDE ref_bitmask_energy - !$OMP PARALLEL DEFAULT(NONE) & - !$OMP PRIVATE(i,hij,j,k,idx,jj,vt,ii,sh,sh2,ni,exa,ext,org_i,org_j,endi) & - !$OMP SHARED(n,H_jj,u_0,keys_tmp,Nint,v_0,davidson_threshold,sorted,shortcut,sort_idx,version,davidson_criterion_is_built) - allocate(idx(0:n), vt(n)) - Vt = 0.d0 + PROVIDE ref_bitmask_energy davidson_criterion + + allocate (shortcut(0:n+1), sort_idx(n), sorted(Nint,n), version(Nint,n)) v_0 = 0.d0 - - !$OMP SINGLE + call sort_dets_ab_v(keys_tmp, sorted, sort_idx, shortcut, version, n, Nint) - !$OMP END SINGLE + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(i,hij,j,k,jj,vt,ii,sh,sh2,ni,exa,ext,org_i,org_j,endi,local_threshold,sorted_i)& + !$OMP SHARED(n,H_jj,u_0,keys_tmp,Nint,v_0,threshold_davidson,sorted,shortcut,sort_idx,version) + allocate(vt(n)) + Vt = 0.d0 !$OMP DO SCHEDULE(dynamic) do sh=1,shortcut(0) - do sh2=1,sh - exa = 0 - do ni=1,Nint - exa += popcnt(xor(version(ni,sh), version(ni,sh2))) - end do - if(exa > 2) then - cycle - end if - - do i=shortcut(sh),shortcut(sh+1)-1 - if(sh==sh2) then - endi = i-1 - else - endi = shortcut(sh2+1)-1 + do sh2=1,sh + exa = 0 + do ni=1,Nint + exa = exa + popcnt(xor(version(ni,sh), version(ni,sh2))) + end do + if(exa > 2) then + cycle end if - do j=shortcut(sh2),endi - ext = exa + do i=shortcut(sh),shortcut(sh+1)-1 + org_i = sort_idx(i) + local_threshold = threshold_davidson - dabs(u_0(org_i)) + if(sh==sh2) then + endi = i-1 + else + endi = shortcut(sh2+1)-1 + end if do ni=1,Nint - ext += popcnt(xor(sorted(ni,i), sorted(ni,j))) - end do - if(ext <= 4) then - org_i = sort_idx(i) + sorted_i(ni) = sorted(ni,i) + enddo + + do j=shortcut(sh2),endi org_j = sort_idx(j) - if ( dabs(u_0(org_j)) + dabs(u_0(org_i)) > davidson_threshold ) then - call i_H_j(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),Nint,hij) - vt (org_i) = vt (org_i) + hij*u_0(org_j) - vt (org_j) = vt (org_j) + hij*u_0(org_i) + if ( dabs(u_0(org_j)) > local_threshold ) then + ext = exa + do ni=1,Nint + ext = ext + popcnt(xor(sorted_i(ni), sorted(ni,j))) + end do + if(ext <= 4) then + call i_H_j(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),Nint,hij) + vt (org_i) = vt (org_i) + hij*u_0(org_j) + vt (org_j) = vt (org_j) + hij*u_0(org_i) + endif endif - endif + enddo enddo enddo enddo - enddo !$OMP END DO - !$OMP SINGLE + !$OMP CRITICAL + do i=1,n + v_0(i) = v_0(i) + vt(i) + enddo + !$OMP END CRITICAL + + deallocate(vt) + !$OMP END PARALLEL + call sort_dets_ba_v(keys_tmp, sorted, sort_idx, shortcut, version, n, Nint) - !$OMP END SINGLE - + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(i,hij,j,k,jj,vt,ii,sh,sh2,ni,exa,ext,org_i,org_j,endi,local_threshold)& + !$OMP SHARED(n,H_jj,u_0,keys_tmp,Nint,v_0,threshold_davidson,sorted,shortcut,sort_idx,version) + allocate(vt(n)) + Vt = 0.d0 + !$OMP DO SCHEDULE(dynamic) do sh=1,shortcut(0) do i=shortcut(sh),shortcut(sh+1)-1 + local_threshold = threshold_davidson - dabs(u_0(org_i)) + org_i = sort_idx(i) do j=shortcut(sh),i-1 - ext = 0 - do ni=1,Nint - ext += popcnt(xor(sorted(ni,i), sorted(ni,j))) - end do - if(ext == 4) then - org_i = sort_idx(i) - org_j = sort_idx(j) - if ( dabs(u_0(org_j)) + dabs(u_0(org_i)) > davidson_threshold ) then + org_j = sort_idx(j) + if ( dabs(u_0(org_j)) > local_threshold ) then + ext = 0 + do ni=1,Nint + ext = ext + popcnt(xor(sorted(ni,i), sorted(ni,j))) + end do + if(ext == 4) then call i_H_j(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),Nint,hij) vt (org_i) = vt (org_i) + hij*u_0(org_j) vt (org_j) = vt (org_j) + hij*u_0(org_i) @@ -1472,10 +1493,12 @@ subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint) v_0(i) = v_0(i) + vt(i) enddo !$OMP END CRITICAL - deallocate(idx,vt) + deallocate(vt) !$OMP END PARALLEL + do i=1,n v_0(i) += H_jj(i) * u_0(i) enddo + deallocate (shortcut, sort_idx, sorted, version) end diff --git a/src/Determinants/spindeterminants.irp.f b/src/Determinants/spindeterminants.irp.f index b8f1b68c..0ca6301a 100644 --- a/src/Determinants/spindeterminants.irp.f +++ b/src/Determinants/spindeterminants.irp.f @@ -442,13 +442,14 @@ BEGIN_PROVIDER [ double precision, psi_bilinear_matrix, (N_det_alpha_unique,N_de enddo END_PROVIDER -subroutine create_wf_of_psi_bilinear_matrix +subroutine create_wf_of_psi_bilinear_matrix(truncate) use bitmasks implicit none BEGIN_DOC ! Generate a wave function containing all possible products ! of alpha and beta determinants END_DOC + logical, intent(in) :: truncate integer :: i,j,k integer(bit_kind) :: tmp_det(N_int,2) integer :: idx @@ -488,8 +489,10 @@ subroutine create_wf_of_psi_bilinear_matrix norm(1) = 0.d0 do i=1,N_det norm(1) += psi_average_norm_contrib_sorted(i) - if (norm(1) >= 0.999999d0) then - exit + if (truncate) then + if (norm(1) >= 0.999999d0) then + exit + endif endif enddo N_det = min(i,N_det) @@ -532,7 +535,6 @@ subroutine generate_all_alpha_beta_det_products !$OMP END DO NOWAIT deallocate(tmp_det) !$OMP END PARALLEL - deallocate (tmp_det) call copy_H_apply_buffer_to_wf SOFT_TOUCH psi_det psi_coef N_det end diff --git a/src/Determinants/truncate_wf.irp.f b/src/Determinants/truncate_wf.irp.f index f867ad7e..42340c71 100644 --- a/src/Determinants/truncate_wf.irp.f +++ b/src/Determinants/truncate_wf.irp.f @@ -8,10 +8,10 @@ program cisd N_det=10000 do i=1,N_det do k=1,N_int - psi_det(k,1,i) = psi_det_sorted(k,1,i) - psi_det(k,2,i) = psi_det_sorted(k,2,i) + psi_det(k,1,i) = psi_det_sorted(k,1,i) + psi_det(k,2,i) = psi_det_sorted(k,2,i) enddo - psi_coef(k,:) = psi_coef_sorted(k,:) + psi_coef(i,:) = psi_coef_sorted(i,:) enddo TOUCH psi_det psi_coef psi_det_sorted psi_coef_sorted psi_average_norm_contrib_sorted N_det call save_wavefunction diff --git a/src/Integrals_Bielec/ao_bi_integrals.irp.f b/src/Integrals_Bielec/ao_bi_integrals.irp.f index aa186454..19f4d2a7 100644 --- a/src/Integrals_Bielec/ao_bi_integrals.irp.f +++ b/src/Integrals_Bielec/ao_bi_integrals.irp.f @@ -204,7 +204,7 @@ double precision function ao_bielec_integral_schwartz_accel(i,j,k,l) integral = general_primitive_integral(dim1, & P_new,P_center,fact_p,pp,p_inv,iorder_p, & Q_new,Q_center,fact_q,qq,q_inv,iorder_q) - ao_bielec_integral_schwartz_accel += coef4 * integral + ao_bielec_integral_schwartz_accel = ao_bielec_integral_schwartz_accel + coef4 * integral enddo ! s enddo ! r enddo ! q @@ -264,7 +264,7 @@ double precision function ao_bielec_integral_schwartz_accel(i,j,k,l) I_power(1),J_power(1),K_power(1),L_power(1), & I_power(2),J_power(2),K_power(2),L_power(2), & I_power(3),J_power(3),K_power(3),L_power(3)) - ao_bielec_integral_schwartz_accel += coef4 * integral + ao_bielec_integral_schwartz_accel = ao_bielec_integral_schwartz_accel + coef4 * integral enddo ! s enddo ! r enddo ! q @@ -307,12 +307,20 @@ subroutine compute_ao_bielec_integrals(j,k,l,sze,buffer_value) buffer_value = 0._integral_kind return endif + if (ao_bielec_integral_schwartz(j,l) < thresh ) then + buffer_value = 0._integral_kind + return + endif do i = 1, ao_num if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < thresh) then buffer_value(i) = 0._integral_kind cycle endif + if (ao_bielec_integral_schwartz(i,k)*ao_bielec_integral_schwartz(j,l) < thresh ) then + buffer_value(i) = 0._integral_kind + cycle + endif !DIR$ FORCEINLINE buffer_value(i) = ao_bielec_integral(i,k,j,l) enddo @@ -378,8 +386,9 @@ BEGIN_PROVIDER [ logical, ao_bielec_integrals_in_map ] !$OMP DEFAULT(NONE) & !$OMP SHARED (ao_num, jl_pairs, ao_integrals_map,thresh, & !$OMP cpu_1,wall_1,lock, lmax,n_centers,ao_nucl, & - !$OMP ao_overlap_abs,ao_overlap,abort_here, & - !$OMP wall_0,progress_bar,progress_value) + !$OMP ao_overlap_abs,ao_overlap,abort_here, & + !$OMP wall_0,progress_bar,progress_value, & + !$OMP ao_bielec_integral_schwartz) allocate(buffer_i(size_buffer)) allocate(buffer_value(size_buffer)) @@ -418,6 +427,9 @@ IRP_ENDIF if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < thresh) then cycle endif + if (ao_bielec_integral_schwartz(i,k)*ao_bielec_integral_schwartz(j,l) < thresh ) then + cycle + endif !DIR$ FORCEINLINE integral = ao_bielec_integral(i,k,j,l) if (abs(integral) < thresh) then @@ -665,32 +677,44 @@ double precision function ERI(alpha,beta,delta,gama,a_x,b_x,c_x,d_x,a_y,b_y,c_y, integer :: n_pt_sup double precision :: p,q,denom,coeff double precision :: I_f + integer :: nx,ny,nz include 'Utils/constants.include.F' - if(iand(a_x+b_x+c_x+d_x,1).eq.1.or.iand(a_y+b_y+c_y+d_y,1).eq.1.or.iand(a_z+b_z+c_z+d_z,1).eq.1)then + nx = a_x+b_x+c_x+d_x + if(iand(nx,1) == 1) then ERI = 0.d0 return - else - - ASSERT (alpha >= 0.d0) - ASSERT (beta >= 0.d0) - ASSERT (delta >= 0.d0) - ASSERT (gama >= 0.d0) - p = alpha + beta - q = delta + gama - ASSERT (p+q >= 0.d0) - coeff = pi_5_2 / (p * q * dsqrt(p+q)) - !DIR$ FORCEINLINE - n_pt = n_pt_sup(a_x,b_x,c_x,d_x,a_y,b_y,c_y,d_y,a_z,b_z,c_z,d_z) - - if (n_pt == 0) then - ERI = coeff - return - endif - - call integrale_new(I_f,a_x,b_x,c_x,d_x,a_y,b_y,c_y,d_y,a_z,b_z,c_z,d_z,p,q,n_pt) - - ERI = I_f * coeff endif + + ny = a_y+b_y+c_y+d_y + if(iand(ny,1) == 1) then + ERI = 0.d0 + return + endif + + nz = a_z+b_z+c_z+d_z + if(iand(nz,1) == 1) then + ERI = 0.d0 + return + endif + + ASSERT (alpha >= 0.d0) + ASSERT (beta >= 0.d0) + ASSERT (delta >= 0.d0) + ASSERT (gama >= 0.d0) + p = alpha + beta + q = delta + gama + ASSERT (p+q >= 0.d0) + n_pt = ishft( nx+ny+nz,1 ) + + coeff = pi_5_2 / (p * q * dsqrt(p+q)) + if (n_pt == 0) then + ERI = coeff + return + endif + + call integrale_new(I_f,a_x,b_x,c_x,d_x,a_y,b_y,c_y,d_y,a_z,b_z,c_z,d_z,p,q,n_pt) + + ERI = I_f * coeff end diff --git a/src/Integrals_Bielec/map_integrals.irp.f b/src/Integrals_Bielec/map_integrals.irp.f index 83950c5c..23deec02 100644 --- a/src/Integrals_Bielec/map_integrals.irp.f +++ b/src/Integrals_Bielec/map_integrals.irp.f @@ -291,19 +291,42 @@ double precision function get_mo_bielec_integral(i,j,k,l,map) PROVIDE mo_bielec_integrals_in_map !DIR$ FORCEINLINE call bielec_integrals_index(i,j,k,l,idx) + !DIR$ FORCEINLINE call map_get(map,idx,tmp) get_mo_bielec_integral = dble(tmp) end +double precision function get_mo_bielec_integral_schwartz(i,j,k,l,map) + use map_module + implicit none + BEGIN_DOC + ! Returns one integral in the MO basis + END_DOC + integer, intent(in) :: i,j,k,l + integer(key_kind) :: idx + type(map_type), intent(inout) :: map + real(integral_kind) :: tmp + PROVIDE mo_bielec_integrals_in_map + if (mo_bielec_integral_schwartz(i,k)*mo_bielec_integral_schwartz(j,l) > mo_integrals_threshold) then + !DIR$ FORCEINLINE + call bielec_integrals_index(i,j,k,l,idx) + !DIR$ FORCEINLINE + call map_get(map,idx,tmp) + else + tmp = 0.d0 + endif + get_mo_bielec_integral_schwartz = dble(tmp) +end + double precision function mo_bielec_integral(i,j,k,l) implicit none BEGIN_DOC ! Returns one integral in the MO basis END_DOC integer, intent(in) :: i,j,k,l - double precision :: get_mo_bielec_integral + double precision :: get_mo_bielec_integral_schwartz PROVIDE mo_bielec_integrals_in_map - mo_bielec_integral = get_mo_bielec_integral(i,j,k,l,mo_integrals_map) + mo_bielec_integral = get_mo_bielec_integral_schwartz(i,j,k,l,mo_integrals_map) return end @@ -502,7 +525,8 @@ integer function load_$ao_integrals(filename) integer*8 :: i integer(cache_key_kind), pointer :: key(:) real(integral_kind), pointer :: val(:) - integer :: iknd, kknd, n, j + integer :: iknd, kknd + integer*8 :: n, j load_$ao_integrals = 1 open(unit=66,file=filename,FORM='unformatted',STATUS='UNKNOWN') read(66,err=98,end=98) iknd, kknd @@ -532,12 +556,8 @@ integer function load_$ao_integrals(filename) return 99 continue call map_deinit($ao_integrals_map) - FREE $ao_integrals_map - if (.True.) then - PROVIDE $ao_integrals_map - endif - stop 'Problem reading $ao_integrals_map file in work/' 98 continue + stop 'Problem reading $ao_integrals_map file in work/' end diff --git a/src/Integrals_Bielec/mo_bi_integrals.irp.f b/src/Integrals_Bielec/mo_bi_integrals.irp.f index 7d772ac3..98737e2a 100644 --- a/src/Integrals_Bielec/mo_bi_integrals.irp.f +++ b/src/Integrals_Bielec/mo_bi_integrals.irp.f @@ -488,3 +488,19 @@ END_PROVIDER enddo END_PROVIDER + +BEGIN_PROVIDER [ double precision, mo_bielec_integral_schwartz,(mo_tot_num,mo_tot_num) ] + implicit none + BEGIN_DOC + ! Needed to compute Schwartz inequalities + END_DOC + + integer :: i,k + + do i=1,mo_tot_num + do k=1,mo_tot_num + mo_bielec_integral_schwartz(k,i) = dsqrt(mo_bielec_integral_jj(k,i)) + enddo + enddo + +END_PROVIDER diff --git a/src/Integrals_Bielec/test_integrals.irp.f b/src/Integrals_Bielec/test_integrals.irp.f new file mode 100644 index 00000000..8d590f27 --- /dev/null +++ b/src/Integrals_Bielec/test_integrals.irp.f @@ -0,0 +1,136 @@ +program bench_maps + implicit none + BEGIN_DOC +! Performs timing benchmarks on integral access + END_DOC + integer :: i,j,k,l + integer*8 :: ii,jj + double precision :: r, cpu + integer*8 :: cpu0, cpu1, count_rate, count_max + + PROVIDE mo_bielec_integrals_in_map + print *, mo_tot_num, 'MOs' + + ! Time random function + call system_clock(cpu0, count_rate, count_max) + do ii=1,100000000_8 + call random_number(r) + i = int(r*mo_tot_num)+1 + call random_number(r) + j = int(r*mo_tot_num)+1 + call random_number(r) + k = int(r*mo_tot_num)+1 + call random_number(r) + l = int(r*mo_tot_num)+1 + enddo + call system_clock(cpu1, count_rate, count_max) + cpu = (cpu1-cpu0)/count_rate + print *, 'Random function : ', cpu/dble(ii) + + call system_clock(cpu0, count_rate, count_max) + do ii=1,100000000_8 + call random_number(r) + i = int(r*mo_tot_num)+1 + call random_number(r) + j = int(r*mo_tot_num)+1 + call random_number(r) + k = int(r*mo_tot_num)+1 + call random_number(r) + l = int(r*mo_tot_num)+1 + call get_mo_bielec_integral(i,j,k,l,mo_integrals_map) + enddo + call system_clock(cpu1, count_rate, count_max) + cpu = -cpu + (cpu1 - cpu0)/count_rate + print *, 'Random access : ', cpu/dble(ii) + + ii=0_8 + call system_clock(cpu0, count_rate, count_max) + do jj=1,10 + do l=1,mo_tot_num + do k=1,mo_tot_num + do j=1,mo_tot_num + do i=1,mo_tot_num + ii += 1 + call get_mo_bielec_integral(i,j,k,l,mo_integrals_map) + enddo + enddo + enddo + enddo + enddo + call system_clock(cpu1, count_rate, count_max) + cpu = (cpu1 - cpu0)/count_rate + print *, 'loop ijkl : ', cpu/dble(ii) + + ii=0 + call system_clock(cpu0, count_rate, count_max) + do jj=1,10 + do l=1,mo_tot_num + do j=1,mo_tot_num + do k=1,mo_tot_num + do i=1,mo_tot_num + ii += 1 + call get_mo_bielec_integral(i,j,k,l,mo_integrals_map) + enddo + enddo + enddo + enddo + enddo + call system_clock(cpu1, count_rate, count_max) + cpu = (cpu1 - cpu0)/count_rate + print *, 'loop ikjl : ', cpu/dble(ii) + + ii=0 + call system_clock(cpu0, count_rate, count_max) + do jj=1,10 + do k=1,mo_tot_num + do l=1,mo_tot_num + do j=1,mo_tot_num + do i=1,mo_tot_num + ii += 1 + call get_mo_bielec_integral(i,j,k,l,mo_integrals_map) + enddo + enddo + enddo + enddo + enddo + call system_clock(cpu1, count_rate, count_max) + cpu = (cpu1 - cpu0)/count_rate + print *, 'loop ijlk : ', cpu/dble(ii) + + ii=0 + call system_clock(cpu0, count_rate, count_max) + do jj=1,10 + do i=1,mo_tot_num + do j=1,mo_tot_num + do k=1,mo_tot_num + do l=1,mo_tot_num + ii += 1 + call get_mo_bielec_integral(i,j,k,l,mo_integrals_map) + enddo + enddo + enddo + enddo + enddo + call system_clock(cpu1, count_rate, count_max) + cpu = (cpu1 - cpu0)/count_rate + print *, 'loop lkji : ', cpu/dble(ii) + + ii=0 + call system_clock(cpu0, count_rate, count_max) + do jj=1,10 + do j=1,mo_tot_num + do i=1,mo_tot_num + do k=1,mo_tot_num + do l=1,mo_tot_num + ii += 1 + call get_mo_bielec_integral(i,j,k,l,mo_integrals_map) + enddo + enddo + enddo + enddo + enddo + call system_clock(cpu1, count_rate, count_max) + cpu = (cpu1 - cpu0)/count_rate + print *, 'loop lkij : ', cpu/dble(ii) + +end diff --git a/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f b/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f index 95023177..e18bc006 100644 --- a/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f +++ b/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f @@ -5,106 +5,105 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integral, (ao_num_align,ao_num)] END_DOC if (do_pseudo) then ao_pseudo_integral = ao_pseudo_integral_local + ao_pseudo_integral_non_local -! ao_pseudo_integral = ao_pseudo_integral_local -! ao_pseudo_integral = ao_pseudo_integral_non_local else ao_pseudo_integral = 0.d0 endif END_PROVIDER - BEGIN_PROVIDER [ double precision, ao_pseudo_integral_local, (ao_num_align,ao_num)] +BEGIN_PROVIDER [ double precision, ao_pseudo_integral_local, (ao_num_align,ao_num)] implicit none BEGIN_DOC -! Local pseudo-potential + ! Local pseudo-potential END_DOC - double precision :: alpha, beta, gama, delta - integer :: num_A,num_B - double precision :: A_center(3),B_center(3),C_center(3) - integer :: power_A(3),power_B(3) - integer :: i,j,k,l,n_pt_in,m - double precision :: Vloc, Vpseudo - - double precision :: cpu_1, cpu_2, wall_1, wall_2, wall_0 - integer :: thread_num - + double precision :: alpha, beta, gama, delta + integer :: num_A,num_B + double precision :: A_center(3),B_center(3),C_center(3) + integer :: power_A(3),power_B(3) + integer :: i,j,k,l,n_pt_in,m + double precision :: Vloc, Vpseudo + + double precision :: cpu_1, cpu_2, wall_1, wall_2, wall_0 + integer :: thread_num + !$ integer :: omp_get_thread_num + ao_pseudo_integral_local = 0.d0 - - !! Dump array - integer, allocatable :: n_k_dump(:) - double precision, allocatable :: v_k_dump(:), dz_k_dump(:) - - allocate(n_k_dump(1:pseudo_klocmax), v_k_dump(1:pseudo_klocmax), dz_k_dump(1:pseudo_klocmax)) - - - ! _ - ! / _. | _ | - ! \_ (_| | (_ |_| | - ! - - print*, 'Providing the nuclear electron pseudo integrals ' - + + !! Dump array + integer, allocatable :: n_k_dump(:) + double precision, allocatable :: v_k_dump(:), dz_k_dump(:) + + allocate(n_k_dump(1:pseudo_klocmax), v_k_dump(1:pseudo_klocmax), dz_k_dump(1:pseudo_klocmax)) + + print*, 'Providing the nuclear electron pseudo integrals (local)' + call wall_time(wall_1) call cpu_time(cpu_1) - - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,C_center,power_A,power_B, & - !$OMP num_A,num_B,Z,c,n_pt_in, & - !$OMP v_k_dump,n_k_dump, dz_k_dump, & - !$OMP wall_0,wall_2,thread_num) & - !$OMP SHARED (ao_num,ao_prim_num,ao_expo_ordered_transp,ao_power,ao_nucl,nucl_coord,ao_coef_normalized_ordered_transp, & - !$OMP ao_pseudo_integral_local,nucl_num,nucl_charge, & - !$OMP pseudo_klocmax,pseudo_lmax,pseudo_kmax,pseudo_v_k,pseudo_n_k, pseudo_dz_k, & - !$OMP wall_1) + thread_num = 0 + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,C_center,power_A,power_B,& + !$OMP num_A,num_B,Z,c,n_pt_in, & + !$OMP v_k_dump,n_k_dump, dz_k_dump, & + !$OMP wall_0,wall_2,thread_num) & + !$OMP SHARED (ao_num,ao_prim_num,ao_expo_ordered_transp,ao_power,ao_nucl,nucl_coord,ao_coef_normalized_ordered_transp,& + !$OMP ao_pseudo_integral_local,nucl_num,nucl_charge, & + !$OMP pseudo_klocmax,pseudo_lmax,pseudo_kmax,pseudo_v_k,pseudo_n_k, pseudo_dz_k,& + !$OMP wall_1) + + !$ thread_num = omp_get_thread_num() !$OMP DO SCHEDULE (guided) - + do j = 1, ao_num - - num_A = ao_nucl(j) - power_A(1:3)= ao_power(j,1:3) - A_center(1:3) = nucl_coord(num_A,1:3) - - do i = 1, ao_num - - num_B = ao_nucl(i) - power_B(1:3)= ao_power(i,1:3) - B_center(1:3) = nucl_coord(num_B,1:3) - - do l=1,ao_prim_num(j) - alpha = ao_expo_ordered_transp(l,j) - - do m=1,ao_prim_num(i) - beta = ao_expo_ordered_transp(m,i) - double precision :: c - c = 0.d0 - - do k = 1, nucl_num - double precision :: Z - Z = nucl_charge(k) - - C_center(1:3) = nucl_coord(k,1:3) - - v_k_dump = pseudo_v_k(k,1:pseudo_klocmax) - n_k_dump = pseudo_n_k(k,1:pseudo_klocmax) - dz_k_dump = pseudo_dz_k(k,1:pseudo_klocmax) - - c = c + Vloc(pseudo_klocmax, v_k_dump,n_k_dump, dz_k_dump, & - A_center,power_A,alpha,B_center,power_B,beta,C_center) - + + num_A = ao_nucl(j) + power_A(1:3)= ao_power(j,1:3) + A_center(1:3) = nucl_coord(num_A,1:3) + + do i = 1, ao_num + + num_B = ao_nucl(i) + power_B(1:3)= ao_power(i,1:3) + B_center(1:3) = nucl_coord(num_B,1:3) + + do l=1,ao_prim_num(j) + alpha = ao_expo_ordered_transp(l,j) + + do m=1,ao_prim_num(i) + beta = ao_expo_ordered_transp(m,i) + double precision :: c + c = 0.d0 + + if (dabs(ao_coef_normalized_ordered_transp(l,j)*ao_coef_normalized_ordered_transp(m,i))& + < 1.d-10) then + cycle + endif + do k = 1, nucl_num + double precision :: Z + Z = nucl_charge(k) + + C_center(1:3) = nucl_coord(k,1:3) + + v_k_dump = pseudo_v_k(k,1:pseudo_klocmax) + n_k_dump = pseudo_n_k(k,1:pseudo_klocmax) + dz_k_dump = pseudo_dz_k(k,1:pseudo_klocmax) + + c = c + Vloc(pseudo_klocmax, v_k_dump,n_k_dump, dz_k_dump,& + A_center,power_A,alpha,B_center,power_B,beta,C_center) + + enddo + ao_pseudo_integral_local(i,j) = ao_pseudo_integral_local(i,j) +& + ao_coef_normalized_ordered_transp(l,j)*ao_coef_normalized_ordered_transp(m,i)*c + enddo enddo - ao_pseudo_integral_local(i,j) = ao_pseudo_integral_local(i,j) + & - ao_coef_normalized_ordered_transp(l,j)*ao_coef_normalized_ordered_transp(m,i)*c - enddo - enddo - enddo - + enddo + call wall_time(wall_2) if (thread_num == 0) then if (wall_2 - wall_0 > 1.d0) then wall_0 = wall_2 - print*, 100.*float(j)/float(ao_num), '% in ', & - wall_2-wall_1, 's' + print*, 100.*float(j)/float(ao_num), '% in ', & + wall_2-wall_1, 's' endif endif enddo @@ -121,106 +120,108 @@ END_PROVIDER BEGIN_PROVIDER [ double precision, ao_pseudo_integral_non_local, (ao_num_align,ao_num)] implicit none BEGIN_DOC -! Local pseudo-potential + ! Local pseudo-potential END_DOC - double precision :: alpha, beta, gama, delta - integer :: num_A,num_B - double precision :: A_center(3),B_center(3),C_center(3) - integer :: power_A(3),power_B(3) - integer :: i,j,k,l,n_pt_in,m - double precision :: Vloc, Vpseudo - - double precision :: cpu_1, cpu_2, wall_1, wall_2, wall_0 - integer :: thread_num - + double precision :: alpha, beta, gama, delta + integer :: num_A,num_B + double precision :: A_center(3),B_center(3),C_center(3) + integer :: power_A(3),power_B(3) + integer :: i,j,k,l,n_pt_in,m + double precision :: Vloc, Vpseudo + !$ integer :: omp_get_thread_num + + double precision :: cpu_1, cpu_2, wall_1, wall_2, wall_0 + integer :: thread_num + ao_pseudo_integral_non_local = 0.d0 - - !! Dump array - integer, allocatable :: n_kl_dump(:,:) - double precision, allocatable :: v_kl_dump(:,:), dz_kl_dump(:,:) - - allocate(n_kl_dump(pseudo_kmax,0:pseudo_lmax), v_kl_dump(pseudo_kmax,0:pseudo_lmax), dz_kl_dump(pseudo_kmax,0:pseudo_lmax)) - - ! _ - ! / _. | _ | - ! \_ (_| | (_ |_| | - ! - - print*, 'Providing the nuclear electron pseudo integrals ' - + + !! Dump array + integer, allocatable :: n_kl_dump(:,:) + double precision, allocatable :: v_kl_dump(:,:), dz_kl_dump(:,:) + + allocate(n_kl_dump(pseudo_kmax,0:pseudo_lmax), v_kl_dump(pseudo_kmax,0:pseudo_lmax), dz_kl_dump(pseudo_kmax,0:pseudo_lmax)) + + print*, 'Providing the nuclear electron pseudo integrals (non-local)' + call wall_time(wall_1) call cpu_time(cpu_1) - - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,C_center,power_A,power_B, & - !$OMP num_A,num_B,Z,c,n_pt_in, & - !$OMP n_kl_dump, v_kl_dump, dz_kl_dump, & - !$OMP wall_0,wall_2,thread_num) & - !$OMP SHARED (ao_num,ao_prim_num,ao_expo_ordered_transp,ao_power,ao_nucl,nucl_coord,ao_coef_normalized_ordered_transp, & - !$OMP ao_pseudo_integral_non_local,nucl_num,nucl_charge, & - !$OMP pseudo_klocmax,pseudo_lmax,pseudo_kmax,pseudo_n_kl, pseudo_v_kl, pseudo_dz_kl, & - !$OMP wall_1) + thread_num = 0 + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,C_center,power_A,power_B,& + !$OMP num_A,num_B,Z,c,n_pt_in, & + !$OMP n_kl_dump, v_kl_dump, dz_kl_dump, & + !$OMP wall_0,wall_2,thread_num) & + !$OMP SHARED (ao_num,ao_prim_num,ao_expo_ordered_transp,ao_power,ao_nucl,nucl_coord,ao_coef_normalized_ordered_transp,& + !$OMP ao_pseudo_integral_non_local,nucl_num,nucl_charge,& + !$OMP pseudo_klocmax,pseudo_lmax,pseudo_kmax,pseudo_n_kl, pseudo_v_kl, pseudo_dz_kl,& + !$OMP wall_1) + !$ thread_num = omp_get_thread_num() !$OMP DO SCHEDULE (guided) - + do j = 1, ao_num - - num_A = ao_nucl(j) - power_A(1:3)= ao_power(j,1:3) - A_center(1:3) = nucl_coord(num_A,1:3) - - do i = 1, ao_num - - num_B = ao_nucl(i) - power_B(1:3)= ao_power(i,1:3) - B_center(1:3) = nucl_coord(num_B,1:3) - - do l=1,ao_prim_num(j) - alpha = ao_expo_ordered_transp(l,j) - - do m=1,ao_prim_num(i) - beta = ao_expo_ordered_transp(m,i) - double precision :: c - c = 0.d0 - - do k = 1, nucl_num - double precision :: Z - Z = nucl_charge(k) - - C_center(1:3) = nucl_coord(k,1:3) - - n_kl_dump = pseudo_n_kl(k,1:pseudo_kmax,0:pseudo_lmax) - v_kl_dump = pseudo_v_kl(k,1:pseudo_kmax,0:pseudo_lmax) - dz_kl_dump = pseudo_dz_kl(k,1:pseudo_kmax,0:pseudo_lmax) - - c = c + Vpseudo(pseudo_lmax,pseudo_kmax,v_kl_dump,n_kl_dump,dz_kl_dump,A_center,power_A,alpha,B_center,power_B,beta,C_center) - + + num_A = ao_nucl(j) + power_A(1:3)= ao_power(j,1:3) + A_center(1:3) = nucl_coord(num_A,1:3) + + do i = 1, ao_num + + num_B = ao_nucl(i) + power_B(1:3)= ao_power(i,1:3) + B_center(1:3) = nucl_coord(num_B,1:3) + + do l=1,ao_prim_num(j) + alpha = ao_expo_ordered_transp(l,j) + + do m=1,ao_prim_num(i) + beta = ao_expo_ordered_transp(m,i) + double precision :: c + c = 0.d0 + + if (dabs(ao_coef_normalized_ordered_transp(l,j)*ao_coef_normalized_ordered_transp(m,i))& + < 1.d-10) then + cycle + endif + + do k = 1, nucl_num + double precision :: Z + Z = nucl_charge(k) + + C_center(1:3) = nucl_coord(k,1:3) + + n_kl_dump = pseudo_n_kl(k,1:pseudo_kmax,0:pseudo_lmax) + v_kl_dump = pseudo_v_kl(k,1:pseudo_kmax,0:pseudo_lmax) + dz_kl_dump = pseudo_dz_kl(k,1:pseudo_kmax,0:pseudo_lmax) + + c = c + Vpseudo(pseudo_lmax,pseudo_kmax,v_kl_dump,n_kl_dump,dz_kl_dump,A_center,power_A,alpha,B_center,power_B,beta,C_center) + + enddo + ao_pseudo_integral_non_local(i,j) = ao_pseudo_integral_non_local(i,j) +& + ao_coef_normalized_ordered_transp(l,j)*ao_coef_normalized_ordered_transp(m,i)*c + enddo enddo - ao_pseudo_integral_non_local(i,j) = ao_pseudo_integral_non_local(i,j) + & - ao_coef_normalized_ordered_transp(l,j)*ao_coef_normalized_ordered_transp(m,i)*c - enddo - enddo - enddo - + enddo + call wall_time(wall_2) if (thread_num == 0) then if (wall_2 - wall_0 > 1.d0) then wall_0 = wall_2 - print*, 100.*float(j)/float(ao_num), '% in ', & - wall_2-wall_1, 's' + print*, 100.*float(j)/float(ao_num), '% in ', & + wall_2-wall_1, 's' endif endif enddo - - !$OMP END DO - !$OMP END PARALLEL - - + + !$OMP END DO + !$OMP END PARALLEL + + deallocate(n_kl_dump,v_kl_dump, dz_kl_dump) - - - END_PROVIDER + + +END_PROVIDER diff --git a/src/Integrals_Monoelec/pseudopot.f90 b/src/Integrals_Monoelec/pseudopot.f90 index e05b883b..c262105f 100644 --- a/src/Integrals_Monoelec/pseudopot.f90 +++ b/src/Integrals_Monoelec/pseudopot.f90 @@ -197,8 +197,8 @@ integer, intent(in) :: n_a(3),n_b(3) double precision, intent(in) :: v_kl(kmax,0:lmax),dz_kl(kmax,0:lmax) ! -! | _ _ _. | _ -! |_ (_) (_ (_| | (/_ +! | _ _ _. | +! |_ (_) (_ (_| | ! double precision :: fourpi,f,prod,prodp,binom_func,accu,bigR,bigI,ylm @@ -223,11 +223,6 @@ double precision, allocatable :: array_I_B(:,:,:,:,:) double precision :: f1, f2, f3 -! _ -! / _. | _ | -! \_ (_| | (_ |_| | -! - if (kmax.eq.1.and.lmax.eq.0.and.v_kl(1,0).eq.0.d0) then Vpseudo=0.d0 return @@ -236,7 +231,7 @@ end if fourpi=4.d0*dacos(-1.d0) ac=dsqrt((a(1)-c(1))**2+(a(2)-c(2))**2+(a(3)-c(3))**2) bc=dsqrt((b(1)-c(1))**2+(b(2)-c(2))**2+(b(3)-c(3))**2) -arg=g_a*ac**2+g_b*bc**2 +arg= g_a*ac*ac + g_b*bc*bc if(arg.gt.-dlog(1.d-20))then Vpseudo=0.d0 @@ -290,6 +285,21 @@ if(ac.eq.0.d0.and.bc.eq.0.d0)then enddo enddo enddo +! do k=1,kmax +! do l=0,lmax +! ktot=ntot+n_kl(k,l) +! do m=-l,l +! prod =bigI(0,0,l,m,n_a(1),n_a(2),n_a(3))*v_kl(k,l) +! prodp=bigI(0,0,l,m,n_b(1),n_b(2),n_b(3))*prod +! if (dabs (prodp) < 1.d-15) then +! cycle +! endif +! +! accu=accu+prodp*int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),0,0,areal,breal,arg) +! +! enddo +! enddo +! enddo !=!=!=!=! ! E n d ! @@ -625,8 +635,8 @@ double precision, intent(in) :: v_kl(kmax,0:lmax),dz_kl(kmax,0:lmax) double precision, intent(in) :: rmax ! -! | _ _ _. | _ -! |_ (_) (_ (_| | (/_ +! | _ _ _. | +! |_ (_) (_ (_| | ! integer :: l,m,k,kk @@ -1950,6 +1960,7 @@ double precision function int_prod_bessel(l,gam,n,m,a,b,arg) double precision :: s_q_0, s_q_k, s_0_0, a_over_b_square double precision :: int_prod_bessel_loc double precision :: inverses(0:300) + double precision :: two_qkmp1, qk logical done @@ -2008,19 +2019,18 @@ double precision function int_prod_bessel(l,gam,n,m,a,b,arg) stop 'pseudopot.f90 : q > 300' endif + two_qkmp1 = dble(2*(q+m)+1) + qk = dble(q) do k=0,q-1 - s_q_k = ( dble(2*(q-k+m)+1)*dble(q-k)*inverses(k) ) * s_q_k + s_q_k = ( two_qkmp1*qk*inverses(k) ) * s_q_k sum=sum+s_q_k + two_qkmp1 = two_qkmp1-2.d0 + qk = qk-1.d0 enddo inverses(q) = a_over_b_square/(dble(2*(q+n)+3) * dble(q+1)) ! do k=0,q ! sum=sum+s_q_k ! s_q_k = a_over_b_square * ( dble(2*(q-k+m)+1)*dble(q-k)/(dble(2*(k+n)+3) * dble(k+1)) ) * s_q_k -! enddo - ! Iteration of k -! do k=0,q -! sum=sum+s_q_k -! s_q_k = a_over_b_square * ( dble(2*(q-k+m)+1)*dble(q-k)/(dble(2*(k+n)+3) * dble(k+1)) ) * s_q_k ! enddo int=int+sum