mirror of
https://github.com/LCPQ/quantum_package
synced 2024-07-05 19:05:59 +02:00
commit
6a99f5528d
|
@ -10,7 +10,7 @@
|
||||||
#
|
#
|
||||||
#
|
#
|
||||||
[COMMON]
|
[COMMON]
|
||||||
FC : gfortran -g -ffree-line-length-none -mavx -I .
|
FC : gfortran -g -ffree-line-length-none -I . -static-libgcc
|
||||||
LAPACK_LIB : -llapack -lblas
|
LAPACK_LIB : -llapack -lblas
|
||||||
IRPF90 : irpf90
|
IRPF90 : irpf90
|
||||||
IRPF90_FLAGS : --ninja --align=32
|
IRPF90_FLAGS : --ninja --align=32
|
||||||
|
|
9
configure
vendored
9
configure
vendored
|
@ -438,11 +438,12 @@ def create_ninja_and_rc(l_installed):
|
||||||
print str_info("qp_root"),
|
print str_info("qp_root"),
|
||||||
python_path = [join(QP_ROOT, "scripts"), join(QP_ROOT, "install")]
|
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 dir_ in python_path:
|
||||||
for folder in os.listdir(dir_):
|
for folder in os.listdir(dir_):
|
||||||
path = join(dir_, folder)
|
path = join(dir_, folder)
|
||||||
if os.path.isdir(path):
|
if os.path.isdir(path):
|
||||||
|
path = path.replace(QP_ROOT,"${QP_ROOT}")
|
||||||
l_python.append(path)
|
l_python.append(path)
|
||||||
|
|
||||||
path_ezfio = find_path('ezfio', l_installed, var_for_qp_root=True)
|
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 = [
|
l_rc = [
|
||||||
'export QP_ROOT={0}'.format(QP_ROOT),
|
'export QP_ROOT={0}'.format(QP_ROOT),
|
||||||
'export QP_EZFIO={0}'.format(path_ezfio),
|
'export QP_EZFIO={0}'.format(path_ezfio.replace(QP_ROOT,"${QP_ROOT}")),
|
||||||
'export IRPF90={0}'.format(path_irpf90),
|
'export IRPF90={0}'.format(path_irpf90.replace(QP_ROOT,"${QP_ROOT}")),
|
||||||
'export NINJA={0}'.format(path_ninja),
|
'export NINJA={0}'.format(path_ninja.replace(QP_ROOT,"${QP_ROOT}")),
|
||||||
'export QP_PYTHON={0}'.format(":".join(l_python)), "",
|
'export QP_PYTHON={0}'.format(":".join(l_python)), "",
|
||||||
'export PYTHONPATH="${QP_EZFIO}":"${QP_PYTHON}":"${PYTHONPATH}"',
|
'export PYTHONPATH="${QP_EZFIO}":"${QP_PYTHON}":"${PYTHONPATH}"',
|
||||||
'export PATH="${QP_PYTHON}":"${QP_ROOT}"/bin:"${QP_ROOT}"/ocaml:"${PATH}"',
|
'export PATH="${QP_PYTHON}":"${QP_ROOT}"/bin:"${QP_ROOT}"/ocaml:"${PATH}"',
|
||||||
|
|
|
@ -8,8 +8,8 @@ type t =
|
||||||
coord : Point3d.t ;
|
coord : Point3d.t ;
|
||||||
} with sexp
|
} with sexp
|
||||||
|
|
||||||
(** Read xyz coordinates of the atom with unit u *)
|
(** Read xyz coordinates of the atom *)
|
||||||
let of_string u s =
|
let of_string ~units s =
|
||||||
let buffer = s
|
let buffer = s
|
||||||
|> String.split ~on:' '
|
|> String.split ~on:' '
|
||||||
|> List.filter ~f:(fun x -> x <> "")
|
|> List.filter ~f:(fun x -> x <> "")
|
||||||
|
@ -18,21 +18,21 @@ let of_string u s =
|
||||||
| [ name; charge; x; y; z ] ->
|
| [ name; charge; x; y; z ] ->
|
||||||
{ element = Element.of_string name ;
|
{ element = Element.of_string name ;
|
||||||
charge = Charge.of_string charge ;
|
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 ] ->
|
| [ name; x; y; z ] ->
|
||||||
let e = Element.of_string name in
|
let e = Element.of_string name in
|
||||||
{ element = e ;
|
{ element = e ;
|
||||||
charge = Element.to_charge 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)
|
| _ -> raise (AtomError s)
|
||||||
;;
|
;;
|
||||||
|
|
||||||
let to_string u a =
|
let to_string ~units a =
|
||||||
[ Element.to_string a.element ;
|
[ Element.to_string a.element ;
|
||||||
Charge.to_string a.charge ;
|
Charge.to_string a.charge ;
|
||||||
Point3d.to_string u a.coord ]
|
Point3d.to_string ~units a.coord ]
|
||||||
|> String.concat ~sep:" "
|
|> String.concat ~sep:" "
|
||||||
;;
|
;;
|
||||||
|
|
||||||
|
|
|
@ -5,5 +5,5 @@ type t = { element : Element.t; charge : Charge.t; coord : Point3d.t; }
|
||||||
val t_of_sexp : Sexplib.Sexp.t -> t
|
val t_of_sexp : Sexplib.Sexp.t -> t
|
||||||
val sexp_of_t : t -> Sexplib.Sexp.t
|
val sexp_of_t : t -> Sexplib.Sexp.t
|
||||||
|
|
||||||
val of_string : Units.units -> string -> t
|
val of_string : units:Units.units -> string -> t
|
||||||
val to_string : Units.units -> t -> string
|
val to_string : units:Units.units -> t -> string
|
||||||
|
|
218
ocaml/Element.ml
218
ocaml/Element.ml
|
@ -1,4 +1,5 @@
|
||||||
open Core.Std;;
|
open Core.Std
|
||||||
|
open Qptypes
|
||||||
|
|
||||||
exception ElementError of string
|
exception ElementError of string
|
||||||
|
|
||||||
|
@ -8,49 +9,49 @@ type t =
|
||||||
|Li|Be |B |C |N |O |F |Ne
|
|Li|Be |B |C |N |O |F |Ne
|
||||||
|Na|Mg |Al|Si|P |S |Cl|Ar
|
|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
|
|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 =
|
let of_string x =
|
||||||
match (String.capitalize (String.lowercase x)) with
|
match (String.capitalize (String.lowercase x)) with
|
||||||
| "X" | "Dummy" -> X
|
| "X" | "Dummy" -> X
|
||||||
| "H" | "Hydrogen" -> H
|
| "H" | "Hydrogen" -> H
|
||||||
| "He" | "Helium" -> He
|
| "He" | "Helium" -> He
|
||||||
| "Li" | "Lithium" -> Li
|
| "Li" | "Lithium" -> Li
|
||||||
| "Be" | "Beryllium" -> Be
|
| "Be" | "Beryllium" -> Be
|
||||||
| "B" | "Boron" -> B
|
| "B" | "Boron" -> B
|
||||||
| "C" | "Carbon" -> C
|
| "C" | "Carbon" -> C
|
||||||
| "N" | "Nitrogen" -> N
|
| "N" | "Nitrogen" -> N
|
||||||
| "O" | "Oxygen" -> O
|
| "O" | "Oxygen" -> O
|
||||||
| "F" | "Fluorine" -> F
|
| "F" | "Fluorine" -> F
|
||||||
| "Ne" | "Neon" -> Ne
|
| "Ne" | "Neon" -> Ne
|
||||||
| "Na" | "Sodium" -> Na
|
| "Na" | "Sodium" -> Na
|
||||||
| "Mg" | "Magnesium" -> Mg
|
| "Mg" | "Magnesium" -> Mg
|
||||||
| "Al" | "Aluminum" -> Al
|
| "Al" | "Aluminum" -> Al
|
||||||
| "Si" | "Silicon" -> Si
|
| "Si" | "Silicon" -> Si
|
||||||
| "P" | "Phosphorus" -> P
|
| "P" | "Phosphorus" -> P
|
||||||
| "S" | "Sulfur" -> S
|
| "S" | "Sulfur" -> S
|
||||||
| "Cl" | "Chlorine" -> Cl
|
| "Cl" | "Chlorine" -> Cl
|
||||||
| "Ar" | "Argon" -> Ar
|
| "Ar" | "Argon" -> Ar
|
||||||
| "K" | "Potassium" -> K
|
| "K" | "Potassium" -> K
|
||||||
| "Ca" | "Calcium" -> Ca
|
| "Ca" | "Calcium" -> Ca
|
||||||
| "Sc" | "Scandium" -> Sc
|
| "Sc" | "Scandium" -> Sc
|
||||||
| "Ti" | "Titanium" -> Ti
|
| "Ti" | "Titanium" -> Ti
|
||||||
| "V" | "Vanadium" -> V
|
| "V" | "Vanadium" -> V
|
||||||
| "Cr" | "Chromium" -> Cr
|
| "Cr" | "Chromium" -> Cr
|
||||||
| "Mn" | "Manganese" -> Mn
|
| "Mn" | "Manganese" -> Mn
|
||||||
| "Fe" | "Iron" -> Fe
|
| "Fe" | "Iron" -> Fe
|
||||||
| "Co" | "Cobalt" -> Co
|
| "Co" | "Cobalt" -> Co
|
||||||
| "Ni" | "Nickel" -> Ni
|
| "Ni" | "Nickel" -> Ni
|
||||||
| "Cu" | "Copper" -> Cu
|
| "Cu" | "Copper" -> Cu
|
||||||
| "Zn" | "Zinc" -> Zn
|
| "Zn" | "Zinc" -> Zn
|
||||||
| "Ga" | "Gallium" -> Ga
|
| "Ga" | "Gallium" -> Ga
|
||||||
| "Ge" | "Germanium" -> Ge
|
| "Ge" | "Germanium" -> Ge
|
||||||
| "As" | "Arsenic" -> As
|
| "As" | "Arsenic" -> As
|
||||||
| "Se" | "Selenium" -> Se
|
| "Se" | "Selenium" -> Se
|
||||||
| "Br" | "Bromine" -> Br
|
| "Br" | "Bromine" -> Br
|
||||||
| "Kr" | "Krypton" -> Kr
|
| "Kr" | "Krypton" -> Kr
|
||||||
| x -> raise (ElementError ("Element "^x^" unknown"))
|
| x -> raise (ElementError ("Element "^x^" unknown"))
|
||||||
;;
|
|
||||||
|
|
||||||
let to_string = function
|
let to_string = function
|
||||||
| X -> "X"
|
| X -> "X"
|
||||||
|
@ -90,7 +91,7 @@ let to_string = function
|
||||||
| Se -> "Se"
|
| Se -> "Se"
|
||||||
| Br -> "Br"
|
| Br -> "Br"
|
||||||
| Kr -> "Kr"
|
| Kr -> "Kr"
|
||||||
;;
|
|
||||||
|
|
||||||
let to_long_string = function
|
let to_long_string = function
|
||||||
| X -> "Dummy"
|
| X -> "Dummy"
|
||||||
|
@ -130,7 +131,7 @@ let to_long_string = function
|
||||||
| Se -> "Selenium"
|
| Se -> "Selenium"
|
||||||
| Br -> "Bromine"
|
| Br -> "Bromine"
|
||||||
| Kr -> "Krypton"
|
| Kr -> "Krypton"
|
||||||
;;
|
|
||||||
|
|
||||||
let to_charge c =
|
let to_charge c =
|
||||||
let result = match c with
|
let result = match c with
|
||||||
|
@ -172,7 +173,7 @@ let to_charge c =
|
||||||
| Br -> 35
|
| Br -> 35
|
||||||
| Kr -> 36
|
| Kr -> 36
|
||||||
in Charge.of_int result
|
in Charge.of_int result
|
||||||
;;
|
|
||||||
|
|
||||||
let of_charge c = match (Charge.to_int c) with
|
let of_charge c = match (Charge.to_int c) with
|
||||||
| 0 -> X
|
| 0 -> X
|
||||||
|
@ -213,5 +214,134 @@ let of_charge c = match (Charge.to_int c) with
|
||||||
| 35 -> Br
|
| 35 -> Br
|
||||||
| 36 -> Kr
|
| 36 -> Kr
|
||||||
| x -> raise (ElementError ("Element of charge "^(string_of_int x)^" unknown"))
|
| 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
|
||||||
|
|
||||||
|
|
|
@ -13,6 +13,8 @@ val of_string : string -> t
|
||||||
val to_string : t -> string
|
val to_string : t -> string
|
||||||
val to_long_string : t -> string
|
val to_long_string : t -> string
|
||||||
|
|
||||||
(** get the positive charge *)
|
(** Properties *)
|
||||||
val to_charge : t -> Charge.t
|
val to_charge : t -> Charge.t
|
||||||
val of_charge : Charge.t -> t
|
val of_charge : Charge.t -> t
|
||||||
|
val covalent_radius : t -> Qptypes.Positive_float.t
|
||||||
|
val vdw_radius : t -> Qptypes.Positive_float.t
|
||||||
|
|
|
@ -147,7 +147,7 @@ nucl_coord = %s
|
||||||
(b.nucl_charge |> Array.to_list |> List.map
|
(b.nucl_charge |> Array.to_list |> List.map
|
||||||
~f:(Charge.to_string) |> String.concat ~sep:", " )
|
~f:(Charge.to_string) |> String.concat ~sep:", " )
|
||||||
(b.nucl_coord |> Array.to_list |> List.map
|
(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"
|
Printf.sprintf " %-3s %d %s"
|
||||||
(b.nucl_label.(i) |> Element.to_string)
|
(b.nucl_label.(i) |> Element.to_string)
|
||||||
(b.nucl_charge.(i) |> Charge.to_int )
|
(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"
|
) |> String.concat ~sep:"\n"
|
||||||
in
|
in
|
||||||
Printf.sprintf "
|
Printf.sprintf "
|
||||||
|
|
|
@ -10,27 +10,36 @@ type t = {
|
||||||
} with sexp
|
} with sexp
|
||||||
|
|
||||||
let get_charge { nuclei ; elec_alpha ; elec_beta } =
|
let get_charge { nuclei ; elec_alpha ; elec_beta } =
|
||||||
let result = (Elec_alpha_number.to_int elec_alpha) +
|
let result =
|
||||||
(Elec_beta_number.to_int elec_beta) in
|
(Elec_alpha_number.to_int elec_alpha) +
|
||||||
|
(Elec_beta_number.to_int elec_beta)
|
||||||
|
in
|
||||||
let rec nucl_charge = function
|
let rec nucl_charge = function
|
||||||
| a::rest -> (Charge.to_float a.Atom.charge) +. nucl_charge rest
|
| a::rest -> (Charge.to_float a.Atom.charge) +. nucl_charge rest
|
||||||
| [] -> 0.
|
| [] -> 0.
|
||||||
in
|
in
|
||||||
Charge.of_float (nucl_charge nuclei -. (Float.of_int result))
|
Charge.of_float (nucl_charge nuclei -. (Float.of_int result))
|
||||||
;;
|
|
||||||
|
|
||||||
let get_multiplicity m =
|
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
|
Multiplicity.of_alpha_beta elec_alpha m.elec_beta
|
||||||
;;
|
|
||||||
|
|
||||||
let get_nucl_num m =
|
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
|
Nucl_number.of_int nmax ~max:nmax
|
||||||
;;
|
|
||||||
|
|
||||||
let name m =
|
let name m =
|
||||||
let cm = Charge.to_int (get_charge m) in
|
let cm =
|
||||||
|
get_charge m
|
||||||
|
|> Charge.to_int
|
||||||
|
in
|
||||||
let c =
|
let c =
|
||||||
match cm with
|
match cm with
|
||||||
| 0 -> ""
|
| 0 -> ""
|
||||||
|
@ -39,8 +48,12 @@ let name m =
|
||||||
| i when i>1 -> Printf.sprintf " (%d+)" i
|
| i when i>1 -> Printf.sprintf " (%d+)" i
|
||||||
| i -> Printf.sprintf " (%d-)" (-i)
|
| i -> Printf.sprintf " (%d-)" (-i)
|
||||||
in
|
in
|
||||||
let mult = Multiplicity.to_string (get_multiplicity m) in
|
let mult =
|
||||||
let { nuclei ; elec_alpha ; elec_beta } = m in
|
get_multiplicity m
|
||||||
|
|> Multiplicity.to_string
|
||||||
|
in
|
||||||
|
let { nuclei ; elec_alpha ; elec_beta } = m
|
||||||
|
in
|
||||||
let rec build_list accu = function
|
let rec build_list accu = function
|
||||||
| a::rest ->
|
| a::rest ->
|
||||||
begin
|
begin
|
||||||
|
@ -53,7 +66,9 @@ let name m =
|
||||||
in
|
in
|
||||||
let rec build_name accu = function
|
let rec build_name accu = function
|
||||||
| (a, n)::rest ->
|
| (a, n)::rest ->
|
||||||
let a = Element.to_string a in
|
let a =
|
||||||
|
Element.to_string a
|
||||||
|
in
|
||||||
begin
|
begin
|
||||||
match n with
|
match n with
|
||||||
| 1 -> build_name (a::accu) rest
|
| 1 -> build_name (a::accu) rest
|
||||||
|
@ -64,19 +79,25 @@ let name m =
|
||||||
end
|
end
|
||||||
| [] -> accu
|
| [] -> accu
|
||||||
in
|
in
|
||||||
let result = build_list [] nuclei |> build_name [c ; ", " ; mult]
|
let result =
|
||||||
|
build_list [] nuclei |> build_name [c ; ", " ; mult]
|
||||||
in
|
in
|
||||||
String.concat (result)
|
String.concat (result)
|
||||||
;;
|
|
||||||
|
|
||||||
let to_string m =
|
let to_string m =
|
||||||
let { nuclei ; elec_alpha ; elec_beta } = m in
|
let { nuclei ; elec_alpha ; elec_beta } = m
|
||||||
let n = List.length nuclei in
|
in
|
||||||
let title = name m in
|
let n =
|
||||||
[ Int.to_string n ; title ] @ (List.map ~f:(fun x -> Atom.to_string
|
List.length nuclei
|
||||||
Units.Angstrom x) 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"
|
|> String.concat ~sep:"\n"
|
||||||
;;
|
|
||||||
|
|
||||||
let of_xyz_string
|
let of_xyz_string
|
||||||
?(charge=(Charge.of_int 0)) ?(multiplicity=(Multiplicity.of_int 1))
|
?(charge=(Charge.of_int 0)) ?(multiplicity=(Multiplicity.of_int 1))
|
||||||
|
@ -94,7 +115,9 @@ let of_xyz_string
|
||||||
) + 1 - (Charge.to_int charge)
|
) + 1 - (Charge.to_int charge)
|
||||||
|> Elec_number.of_int
|
|> Elec_number.of_int
|
||||||
in
|
in
|
||||||
let (na,nb) = Multiplicity.to_alpha_beta ne multiplicity in
|
let (na,nb) =
|
||||||
|
Multiplicity.to_alpha_beta ne multiplicity
|
||||||
|
in
|
||||||
let result =
|
let result =
|
||||||
{ nuclei = l ;
|
{ nuclei = l ;
|
||||||
elec_alpha = na ;
|
elec_alpha = na ;
|
||||||
|
@ -109,7 +132,7 @@ let of_xyz_string
|
||||||
raise (MultiplicityError msg);
|
raise (MultiplicityError msg);
|
||||||
else () ;
|
else () ;
|
||||||
result
|
result
|
||||||
;;
|
|
||||||
|
|
||||||
|
|
||||||
let of_xyz_file
|
let of_xyz_file
|
||||||
|
@ -121,8 +144,33 @@ let of_xyz_file
|
||||||
let (_,buffer) = String.lsplit2_exn buffer ~on:'\n' in
|
let (_,buffer) = String.lsplit2_exn buffer ~on:'\n' in
|
||||||
of_xyz_string ~charge:charge ~multiplicity:multiplicity
|
of_xyz_string ~charge:charge ~multiplicity:multiplicity
|
||||||
~units:units buffer
|
~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
|
let to_md5 = to_md5 sexp_of_t
|
||||||
;;
|
|
||||||
|
|
|
@ -34,5 +34,9 @@ val of_xyz_string :
|
||||||
?multiplicity:Multiplicity.t ->
|
?multiplicity:Multiplicity.t ->
|
||||||
?units:Units.units -> string -> 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 *)
|
(** Computes the MD5 hash *)
|
||||||
val to_md5 : t -> Qptypes.MD5.t
|
val to_md5 : t -> Qptypes.MD5.t
|
||||||
|
|
|
@ -7,9 +7,16 @@ type t = {
|
||||||
z : float ;
|
z : float ;
|
||||||
} with sexp
|
} 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 *)
|
(** Read x y z coordinates in string s with units u *)
|
||||||
let of_string u s =
|
let of_string ~units s =
|
||||||
let f = match u with
|
let f = match units with
|
||||||
| Units.Bohr -> 1.
|
| Units.Bohr -> 1.
|
||||||
| Units.Angstrom -> Units.angstrom_to_bohr
|
| Units.Angstrom -> Units.angstrom_to_bohr
|
||||||
in
|
in
|
||||||
|
@ -22,7 +29,6 @@ let of_string u s =
|
||||||
{ x = l.(0) *. f ;
|
{ x = l.(0) *. f ;
|
||||||
y = l.(1) *. f ;
|
y = l.(1) *. f ;
|
||||||
z = l.(2) *. f }
|
z = l.(2) *. f }
|
||||||
;;
|
|
||||||
|
|
||||||
|
|
||||||
let distance2 p1 p2 =
|
let distance2 p1 p2 =
|
||||||
|
@ -30,17 +36,18 @@ let distance2 p1 p2 =
|
||||||
and { x=x2 ; y=y2 ; z=z2 } = p2 in
|
and { x=x2 ; y=y2 ; z=z2 } = p2 in
|
||||||
(x2-.x1)*.(x2-.x1) +. (y2-.y1)*.(y2-.y1) +. (z2-.z1)*.(z2-.z1)
|
(x2-.x1)*.(x2-.x1) +. (y2-.y1)*.(y2-.y1) +. (z2-.z1)*.(z2-.z1)
|
||||||
|> Positive_float.of_float
|
|> Positive_float.of_float
|
||||||
;;
|
|
||||||
|
|
||||||
let distance p1 p2 = sqrt (Positive_float.to_float (distance2 p1 p2))
|
|
||||||
;;
|
|
||||||
|
|
||||||
let to_string u p =
|
let distance p1 p2 =
|
||||||
let f = match u with
|
sqrt (Positive_float.to_float (distance2 p1 p2))
|
||||||
|
|
||||||
|
|
||||||
|
let to_string ~units p =
|
||||||
|
let f = match units with
|
||||||
| Units.Bohr -> 1.
|
| Units.Bohr -> 1.
|
||||||
| Units.Angstrom -> Units.bohr_to_angstrom
|
| Units.Angstrom -> Units.bohr_to_angstrom
|
||||||
in
|
in
|
||||||
let { x=x ; y=y ; z=z } = p in
|
let { x=x ; y=y ; z=z } = p in
|
||||||
Printf.sprintf "%16.8f %16.8f %16.8f" (x*.f) (y*.f) (z*.f)
|
Printf.sprintf "%16.8f %16.8f %16.8f" (x*.f) (y*.f) (z*.f)
|
||||||
;;
|
|
||||||
|
|
||||||
|
|
|
@ -4,11 +4,14 @@ type t =
|
||||||
z : float;
|
z : float;
|
||||||
} with sexp
|
} with sexp
|
||||||
|
|
||||||
|
(** Create from a tuple of floats *)
|
||||||
|
val of_tuple : units:Units.units -> float*float*float -> t
|
||||||
|
|
||||||
(** Create from an xyz string *)
|
(** 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 *)
|
(** 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 *)
|
(** Computes the squared distance between 2 points *)
|
||||||
val distance2 : t -> t -> Qptypes.Positive_float.t
|
val distance2 : t -> t -> Qptypes.Positive_float.t
|
||||||
|
|
|
@ -11,26 +11,92 @@ let spec =
|
||||||
~doc:"string Name of basis set."
|
~doc:"string Name of basis set."
|
||||||
+> flag "c" (optional_with_default 0 int)
|
+> flag "c" (optional_with_default 0 int)
|
||||||
~doc:"int Total charge of the molecule. Default is 0."
|
~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)
|
+> flag "m" (optional_with_default 1 int)
|
||||||
~doc:"int Spin multiplicity (2S+1) of the molecule. Default is 1."
|
~doc:"int Spin multiplicity (2S+1) of the molecule. Default is 1."
|
||||||
+> flag "p" no_arg
|
+> flag "p" no_arg
|
||||||
~doc:" Using pseudopotentials"
|
~doc:" Using pseudopotentials"
|
||||||
+> anon ("xyz_file" %: string)
|
+> 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 *)
|
(* Read molecule *)
|
||||||
let molecule =
|
let molecule =
|
||||||
(Molecule.of_xyz_file xyz_file ~charge:(Charge.of_int c)
|
(Molecule.of_xyz_file xyz_file ~charge:(Charge.of_int c)
|
||||||
~multiplicity:(Multiplicity.of_int m) )
|
~multiplicity:(Multiplicity.of_int m) )
|
||||||
in
|
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
|
let basis_table = Hashtbl.Poly.create () in
|
||||||
(* Open basis set channels *)
|
(* Open basis set channels *)
|
||||||
let basis_channel element =
|
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
|
match Hashtbl.find basis_table key with
|
||||||
| Some in_channel ->
|
| Some in_channel ->
|
||||||
in_channel
|
in_channel
|
||||||
|
@ -80,7 +146,8 @@ let run ?o b c m p xyz_file =
|
||||||
in
|
in
|
||||||
Unix.unlink filename;
|
Unix.unlink filename;
|
||||||
List.iter nuclei ~f:(fun elem->
|
List.iter nuclei ~f:(fun elem->
|
||||||
let key = Element.to_string elem.Atom.element
|
let key =
|
||||||
|
Element.to_string elem.Atom.element
|
||||||
in
|
in
|
||||||
match Hashtbl.add basis_table ~key:key ~data:new_channel with
|
match Hashtbl.add basis_table ~key:key ~data:new_channel with
|
||||||
| `Ok -> ()
|
| `Ok -> ()
|
||||||
|
@ -92,7 +159,8 @@ let run ?o b c m p xyz_file =
|
||||||
let elem = Element.of_string key
|
let elem = Element.of_string key
|
||||||
and basis = String.lowercase basis
|
and basis = String.lowercase basis
|
||||||
in
|
in
|
||||||
let key = Element.to_string elem
|
let key =
|
||||||
|
Element.to_string elem
|
||||||
in
|
in
|
||||||
let command =
|
let command =
|
||||||
Qpackage.root ^ "/scripts/get_basis.sh \"" ^ temp_filename ^
|
Qpackage.root ^ "/scripts/get_basis.sh \"" ^ temp_filename ^
|
||||||
|
@ -175,7 +243,12 @@ let run ?o b c m p xyz_file =
|
||||||
|> List.rev
|
|> List.rev
|
||||||
|> List.map ~f:(fun (x,i) ->
|
|> List.map ~f:(fun (x,i) ->
|
||||||
try
|
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
|
with
|
||||||
| End_of_file ->
|
| End_of_file ->
|
||||||
begin
|
begin
|
||||||
|
@ -264,7 +337,7 @@ let run ?o b c m p xyz_file =
|
||||||
| None -> failwith "Error in basis"
|
| None -> failwith "Error in basis"
|
||||||
| Some x -> Input.Ao_basis.write x
|
| Some x -> Input.Ao_basis.write x
|
||||||
|
|
||||||
;;
|
|
||||||
|
|
||||||
let command =
|
let command =
|
||||||
Command.basic
|
Command.basic
|
||||||
|
@ -279,13 +352,13 @@ elements can be defined as follows:
|
||||||
|
|
||||||
")
|
")
|
||||||
spec
|
spec
|
||||||
(fun o b c m p xyz_file () ->
|
(fun o b c d m p xyz_file () ->
|
||||||
run ?o b c m p xyz_file )
|
run ?o b c d m p xyz_file )
|
||||||
;;
|
|
||||||
|
|
||||||
let () =
|
let () =
|
||||||
Command.run command
|
Command.run command
|
||||||
;;
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
|
@ -39,6 +39,15 @@ H 0.54386314 0.00000000 -0.92559535
|
||||||
let m = Molecule.of_xyz_file "c2h6.xyz" in
|
let m = Molecule.of_xyz_file "c2h6.xyz" in
|
||||||
print_string (Molecule.to_string m);
|
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 ();;
|
test_molecule ();;
|
||||||
|
|
|
@ -10,6 +10,12 @@ doc: Maximum number of SCF iterations
|
||||||
interface: ezfio,provider,ocaml
|
interface: ezfio,provider,ocaml
|
||||||
default: 200
|
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]
|
[mo_guess_type]
|
||||||
type: MO_guess
|
type: MO_guess
|
||||||
doc: Initial MO guess. Can be [ Huckel | HCore ]
|
doc: Initial MO guess. Can be [ Huckel | HCore ]
|
||||||
|
|
|
@ -73,8 +73,13 @@
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
endif
|
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
|
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
|
enddo
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
|
@ -108,9 +113,10 @@ END_PROVIDER
|
||||||
END_DOC
|
END_DOC
|
||||||
|
|
||||||
integer :: i,j,k,l,k1,r,s
|
integer :: i,j,k,l,k1,r,s
|
||||||
|
integer :: i0,j0,k0,l0
|
||||||
integer*8 :: p,q
|
integer*8 :: p,q
|
||||||
double precision :: integral
|
double precision :: integral, c0, c1, c2
|
||||||
double precision :: ao_bielec_integral
|
double precision :: ao_bielec_integral, local_threshold
|
||||||
double precision, allocatable :: ao_bi_elec_integral_alpha_tmp(:,:)
|
double precision, allocatable :: ao_bi_elec_integral_alpha_tmp(:,:)
|
||||||
double precision, allocatable :: ao_bi_elec_integral_beta_tmp(:,:)
|
double precision, allocatable :: ao_bi_elec_integral_beta_tmp(:,:)
|
||||||
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: 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
|
if (do_direct_integrals) then
|
||||||
|
|
||||||
!$OMP PARALLEL DEFAULT(NONE) &
|
!$OMP PARALLEL DEFAULT(NONE) &
|
||||||
!$OMP PRIVATE(i,j,l,k1,k,integral,ii,jj,kk,ll,i8,keys,values,p,q,r,s, &
|
!$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)&
|
!$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 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_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(keys(1), values(1))
|
||||||
allocate(ao_bi_elec_integral_alpha_tmp(ao_num_align,ao_num), &
|
allocate(ao_bi_elec_integral_alpha_tmp(ao_num_align,ao_num), &
|
||||||
|
@ -152,14 +159,16 @@ END_PROVIDER
|
||||||
< ao_integrals_threshold) then
|
< ao_integrals_threshold) then
|
||||||
cycle
|
cycle
|
||||||
endif
|
endif
|
||||||
if (ao_bielec_integral_schwartz(k,l)*ao_bielec_integral_schwartz(i,j) &
|
local_threshold = ao_bielec_integral_schwartz(k,l)*ao_bielec_integral_schwartz(i,j)
|
||||||
< ao_integrals_threshold) then
|
if (local_threshold < ao_integrals_threshold) then
|
||||||
cycle
|
|
||||||
endif
|
|
||||||
values(1) = ao_bielec_integral(k,l,i,j)
|
|
||||||
if (abs(values(1)) < ao_integrals_threshold) then
|
|
||||||
cycle
|
cycle
|
||||||
endif
|
endif
|
||||||
|
i0 = i
|
||||||
|
j0 = j
|
||||||
|
k0 = k
|
||||||
|
l0 = l
|
||||||
|
values(1) = 0.d0
|
||||||
|
local_threshold = ao_integrals_threshold/local_threshold
|
||||||
do k2=1,8
|
do k2=1,8
|
||||||
if (kk(k2)==0) then
|
if (kk(k2)==0) then
|
||||||
cycle
|
cycle
|
||||||
|
@ -168,12 +177,21 @@ END_PROVIDER
|
||||||
j = jj(k2)
|
j = jj(k2)
|
||||||
k = kk(k2)
|
k = kk(k2)
|
||||||
l = ll(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_alpha_tmp(i,j) += integral
|
||||||
ao_bi_elec_integral_beta_tmp (i,j) += integral
|
ao_bi_elec_integral_beta_tmp (i,j) += integral
|
||||||
integral = values(1)
|
integral = values(1)
|
||||||
ao_bi_elec_integral_alpha_tmp(l,j) -= HF_density_matrix_ao_alpha(k,i) * integral
|
ao_bi_elec_integral_alpha_tmp(l,j) -= c1 * integral
|
||||||
ao_bi_elec_integral_beta_tmp (l,j) -= HF_density_matrix_ao_beta (k,i) * integral
|
ao_bi_elec_integral_beta_tmp (l,j) -= c2 * integral
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
!$OMP END DO NOWAIT
|
!$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
|
! Fock matrix in AO basis set
|
||||||
END_DOC
|
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
|
integer :: i,j
|
||||||
do j=1,ao_num
|
do j=1,ao_num
|
||||||
!DIR$ VECTOR ALIGNED
|
!DIR$ VECTOR ALIGNED
|
||||||
|
|
|
@ -42,16 +42,13 @@ subroutine run
|
||||||
BEGIN_DOC
|
BEGIN_DOC
|
||||||
! Run SCF calculation
|
! Run SCF calculation
|
||||||
END_DOC
|
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
|
double precision :: E0
|
||||||
integer :: i_it, i, j, k
|
integer :: i_it, i, j, k
|
||||||
|
|
||||||
E0 = HF_energy
|
E0 = HF_energy
|
||||||
|
|
||||||
thresh_SCF = 1.d-10
|
|
||||||
call damping_SCF
|
|
||||||
mo_label = "Canonical"
|
mo_label = "Canonical"
|
||||||
TOUCH mo_label mo_coef
|
call damping_SCF
|
||||||
call save_mos
|
|
||||||
|
|
||||||
end
|
end
|
||||||
|
|
|
@ -86,7 +86,7 @@ subroutine damping_SCF
|
||||||
if ((E_half > E).and.(E_new < E)) then
|
if ((E_half > E).and.(E_new < E)) then
|
||||||
lambda = 1.d0
|
lambda = 1.d0
|
||||||
exit
|
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
|
lambda = 0.5d0 * lambda
|
||||||
E_new = E_half
|
E_new = E_half
|
||||||
else
|
else
|
||||||
|
|
|
@ -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))
|
(Fock_matrix_diag_mo(p1) + Fock_matrix_diag_mo(p2))
|
||||||
delta_e = 1.d0/delta_e
|
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(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)
|
|
||||||
h = diag_H_mat_elem(det_pert,Nint)
|
h = diag_H_mat_elem(det_pert,Nint)
|
||||||
do i =1,n_st
|
do i =1,n_st
|
||||||
H_pert_diag(i) = h
|
H_pert_diag(i) = h
|
||||||
|
|
|
@ -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_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,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)
|
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)
|
oii = diag_O1_mat_elem_alpha_beta(mo_integrated_delta_rho_one_point,det_pert,N_int)
|
||||||
|
|
|
@ -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_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,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)
|
h = diag_H_mat_elem(det_pert,Nint)
|
||||||
oii = diag_O1_mat_elem(mo_dipole_z,det_pert,N_int)
|
oii = diag_O1_mat_elem(mo_dipole_z,det_pert,N_int)
|
||||||
|
|
|
@ -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 == N_int)
|
||||||
ASSERT (Nint > 0)
|
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,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)
|
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
|
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,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)
|
h = diag_H_mat_elem(det_pert,Nint)
|
||||||
do i =1,N_st
|
do i =1,N_st
|
||||||
|
|
|
@ -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 == N_int)
|
||||||
ASSERT (Nint > 0)
|
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,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)
|
h = diag_H_mat_elem(det_pert,Nint)
|
||||||
|
|
|
@ -62,7 +62,7 @@ program e_curve
|
||||||
endif
|
endif
|
||||||
enddo
|
enddo
|
||||||
call compute_energy(psi_bilinear_matrix_values_save,E,m,norm)
|
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 * m ) / &
|
||||||
dble( elec_alpha_num**3 + elec_alpha_num**2 * n ), norm, E
|
dble( elec_alpha_num**3 + elec_alpha_num**2 * n ), norm, E
|
||||||
if (E < target_energy) then
|
if (E < target_energy) then
|
||||||
|
@ -93,7 +93,7 @@ subroutine compute_energy(psi_bilinear_matrix_values_save, E, m, norm)
|
||||||
m = 0
|
m = 0
|
||||||
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(k,l,det_i,det_j,ci,cj,hij) REDUCTION(+:norm,m,num)
|
!$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))
|
allocate( det_i(N_int,2), det_j(N_int,2))
|
||||||
!$OMP DO
|
!$OMP DO schedule(guided)
|
||||||
do k=1,n_det
|
do k=1,n_det
|
||||||
if (psi_bilinear_matrix_values_save(k) == 0.d0) then
|
if (psi_bilinear_matrix_values_save(k) == 0.d0) then
|
||||||
cycle
|
cycle
|
||||||
|
|
|
@ -35,7 +35,8 @@ except ImportError:
|
||||||
|
|
||||||
from qp_path import QP_ROOT, QP_SRC, QP_EZFIO
|
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")
|
ROOT_BUILD_NINJA = join(QP_ROOT, "config", "build.ninja")
|
||||||
|
|
||||||
header = r"""#
|
header = r"""#
|
||||||
|
@ -94,7 +95,7 @@ def ninja_create_env_variable(pwd_config_file):
|
||||||
l_string.append(str_)
|
l_string.append(str_)
|
||||||
|
|
||||||
lib_lapack = get_compilation_option(pwd_config_file, "LAPACK_LIB")
|
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("")
|
l_string.append("")
|
||||||
|
|
||||||
|
|
|
@ -9,7 +9,6 @@
|
||||||
#
|
#
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
if [[ -z ${QP_ROOT} ]]
|
if [[ -z ${QP_ROOT} ]]
|
||||||
then
|
then
|
||||||
print "The QP_ROOT environment variable is not set."
|
print "The QP_ROOT environment variable is not set."
|
||||||
|
|
|
@ -1,4 +1,5 @@
|
||||||
name_to_elec = {"H": 1,
|
name_to_elec = {"X": 0,
|
||||||
|
"H": 1,
|
||||||
"He": 2,
|
"He": 2,
|
||||||
"Li": 3,
|
"Li": 3,
|
||||||
"Be": 4,
|
"Be": 4,
|
||||||
|
|
|
@ -58,17 +58,35 @@ def get_pseudo_str(l_atom):
|
||||||
str_ = ""
|
str_ = ""
|
||||||
|
|
||||||
for a in l_atom:
|
for a in l_atom:
|
||||||
l_cmd_atom = ["--atom", a]
|
|
||||||
|
|
||||||
l_cmd_head = [EMSL_path, "get_basis_data",
|
if a is not 'X':
|
||||||
"--db_path", db_path,
|
l_cmd_atom = ["--atom", a]
|
||||||
"--basis", "BFD-Pseudo"]
|
|
||||||
|
|
||||||
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()
|
process = Popen(l_cmd_head + l_cmd_atom, stdout=PIPE, stderr=PIPE)
|
||||||
str_ += stdout.strip() + "\n"
|
|
||||||
|
|
||||||
|
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_
|
return str_
|
||||||
|
|
||||||
|
|
||||||
|
|
|
@ -1,7 +1,7 @@
|
||||||
#!/usr/bin/env python
|
#!/usr/bin/env python
|
||||||
|
|
||||||
import subprocess
|
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)
|
shell=True, stdout=subprocess.PIPE)
|
||||||
result = pipe.stdout.read()
|
result = pipe.stdout.read()
|
||||||
is_master_repository = "LCPQ/quantum_package" in result
|
is_master_repository = "LCPQ/quantum_package" in result
|
||||||
|
|
|
@ -40,6 +40,12 @@ doc: Force the wave function to be an eigenfunction of S^2
|
||||||
interface: ezfio,provider,ocaml
|
interface: ezfio,provider,ocaml
|
||||||
default: False
|
default: False
|
||||||
|
|
||||||
|
[threshold_davidson]
|
||||||
|
type: Threshold
|
||||||
|
doc: Thresholds of Davidson's algorithm
|
||||||
|
interface: ezfio,provider,ocaml
|
||||||
|
default: 1.e-8
|
||||||
|
|
||||||
[threshold_generators]
|
[threshold_generators]
|
||||||
type: Threshold
|
type: Threshold
|
||||||
doc: Thresholds on generators (fraction of the norm)
|
doc: Thresholds on generators (fraction of the norm)
|
||||||
|
|
|
@ -2,11 +2,11 @@ use bitmasks
|
||||||
use omp_lib
|
use omp_lib
|
||||||
|
|
||||||
type H_apply_buffer_type
|
type H_apply_buffer_type
|
||||||
integer :: N_det
|
integer :: N_det
|
||||||
integer :: sze
|
integer :: sze
|
||||||
integer(bit_kind), pointer :: det(:,:,:)
|
integer(bit_kind), pointer :: det(:,:,:)
|
||||||
double precision , pointer :: coef(:,:)
|
double precision , pointer :: coef(:,:)
|
||||||
double precision , pointer :: e2(:,:)
|
double precision , pointer :: e2(:,:)
|
||||||
end type H_apply_buffer_type
|
end type H_apply_buffer_type
|
||||||
|
|
||||||
type(H_apply_buffer_type), pointer :: H_apply_buffer(:)
|
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))
|
call omp_init_lock(H_apply_buffer_lock(1,iproc))
|
||||||
!$OMP END PARALLEL
|
!$OMP END PARALLEL
|
||||||
endif
|
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
|
END_PROVIDER
|
||||||
|
|
||||||
|
@ -111,7 +120,6 @@ subroutine copy_H_apply_buffer_to_wf
|
||||||
double precision, allocatable :: buffer_coef(:,:)
|
double precision, allocatable :: buffer_coef(:,:)
|
||||||
integer :: i,j,k
|
integer :: i,j,k
|
||||||
integer :: N_det_old
|
integer :: N_det_old
|
||||||
integer :: iproc
|
|
||||||
|
|
||||||
PROVIDE H_apply_buffer_allocated
|
PROVIDE H_apply_buffer_allocated
|
||||||
|
|
||||||
|
@ -158,7 +166,7 @@ subroutine copy_H_apply_buffer_to_wf
|
||||||
enddo
|
enddo
|
||||||
!$OMP PARALLEL DEFAULT(SHARED) &
|
!$OMP PARALLEL DEFAULT(SHARED) &
|
||||||
!$OMP PRIVATE(j,k,i) FIRSTPRIVATE(N_det_old) &
|
!$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=0
|
||||||
!$ j=omp_get_thread_num()
|
!$ j=omp_get_thread_num()
|
||||||
do k=0,j-1
|
do k=0,j-1
|
||||||
|
|
|
@ -90,54 +90,73 @@ end function
|
||||||
|
|
||||||
subroutine tamiser(key, idx, no, n, Nint, N_key)
|
subroutine tamiser(key, idx, no, n, Nint, N_key)
|
||||||
use bitmasks
|
use bitmasks
|
||||||
|
|
||||||
implicit none
|
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,intent(in) :: no, n, Nint, N_key
|
||||||
|
integer(bit_kind),intent(inout) :: key(Nint, 2, N_key)
|
||||||
integer,intent(inout) :: idx(N_key)
|
integer,intent(inout) :: idx(N_key)
|
||||||
integer :: k,j,tmpidx
|
integer :: k,j,tmpidx
|
||||||
integer(bit_kind) :: tmp(Nint, 2)
|
integer(bit_kind) :: tmp(Nint, 2)
|
||||||
logical :: det_inf
|
logical :: det_inf
|
||||||
|
integer :: ni
|
||||||
|
|
||||||
k = no
|
k = no
|
||||||
j = 2*k
|
j = 2*k
|
||||||
do while(j <= n)
|
do while(j <= n)
|
||||||
if(j < n .and. det_inf(key(:,:,j), key(:,:,j+1), Nint)) then
|
if(j < n) then
|
||||||
j = j+1
|
if (det_inf(key(1,1,j), key(1,1,j+1), Nint)) then
|
||||||
end if
|
j = j+1
|
||||||
if(det_inf(key(:,:,k), key(:,:,j), Nint)) then
|
endif
|
||||||
tmp(:,:) = key(:,:,k)
|
endif
|
||||||
key(:,:,k) = key(:,:,j)
|
if(det_inf(key(1,1,k), key(1,1,j), Nint)) then
|
||||||
key(:,:,j) = tmp(:,:)
|
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)
|
tmpidx = idx(k)
|
||||||
idx(k) = idx(j)
|
idx(k) = idx(j)
|
||||||
idx(j) = tmpidx
|
idx(j) = tmpidx
|
||||||
k = j
|
k = j
|
||||||
j = 2*k
|
j = k+k
|
||||||
else
|
else
|
||||||
return
|
return
|
||||||
end if
|
endif
|
||||||
end do
|
enddo
|
||||||
end subroutine
|
end subroutine
|
||||||
|
|
||||||
|
|
||||||
subroutine sort_dets_ba_v(key_in, key_out, idx, shortcut, version, N_key, Nint)
|
subroutine sort_dets_ba_v(key_in, key_out, idx, shortcut, version, N_key, Nint)
|
||||||
use bitmasks
|
use bitmasks
|
||||||
implicit none
|
implicit none
|
||||||
integer(bit_kind),intent(in) :: key_in(Nint,2,N_key)
|
BEGIN_DOC
|
||||||
integer(bit_kind) :: key(Nint,2,N_key)
|
! Uncodumented : TODO
|
||||||
integer(bit_kind),intent(out) :: key_out(Nint,N_key)
|
END_DOC
|
||||||
integer,intent(out) :: idx(N_key)
|
integer, intent(in) :: Nint, N_key
|
||||||
integer,intent(out) :: shortcut(0:N_key+1)
|
integer(bit_kind),intent(in) :: key_in(Nint,2,N_key)
|
||||||
integer(bit_kind),intent(out) :: version(Nint,N_key+1)
|
integer(bit_kind),intent(out) :: key_out(Nint,N_key)
|
||||||
integer, intent(in) :: Nint, N_key
|
integer,intent(out) :: idx(N_key)
|
||||||
integer(bit_kind) :: tmp(Nint, 2,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)
|
allocate ( key(Nint,2,N_key) )
|
||||||
key(:,2,:N_key) = key_in(:,1,: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)
|
call sort_dets_ab_v(key, key_out, idx, shortcut, version, N_key, Nint)
|
||||||
|
deallocate ( key )
|
||||||
end subroutine
|
end subroutine
|
||||||
|
|
||||||
|
|
||||||
|
@ -146,18 +165,25 @@ subroutine sort_dets_ab_v(key_in, key_out, idx, shortcut, version, N_key, Nint)
|
||||||
use bitmasks
|
use bitmasks
|
||||||
implicit none
|
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),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(bit_kind),intent(out) :: key_out(Nint,N_key)
|
||||||
integer,intent(out) :: idx(N_key)
|
integer,intent(out) :: idx(N_key)
|
||||||
integer,intent(out) :: shortcut(0:N_key+1)
|
integer,intent(out) :: shortcut(0:N_key+1)
|
||||||
integer(bit_kind),intent(out) :: version(Nint,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(bit_kind) :: tmp(Nint, 2)
|
||||||
integer :: tmpidx,i,ni
|
integer :: tmpidx,i,ni
|
||||||
|
|
||||||
key(:,:,:) = key_in(:,:,:)
|
allocate (key(Nint,2,N_key))
|
||||||
do i=1,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
|
idx(i) = i
|
||||||
end do
|
end do
|
||||||
|
|
||||||
|
@ -166,9 +192,14 @@ subroutine sort_dets_ab_v(key_in, key_out, idx, shortcut, version, N_key, Nint)
|
||||||
end do
|
end do
|
||||||
|
|
||||||
do i=N_key,2,-1
|
do i=N_key,2,-1
|
||||||
tmp(:,:) = key(:,:,i)
|
do ni=1,Nint
|
||||||
key(:,:,i) = key(:,:,1)
|
tmp(ni,1) = key(ni,1,i)
|
||||||
key(:,:,1) = tmp(:,:)
|
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)
|
tmpidx = idx(i)
|
||||||
idx(i) = idx(1)
|
idx(i) = idx(1)
|
||||||
idx(1) = tmpidx
|
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(0) = 1
|
||||||
shortcut(1) = 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 i=2,N_key
|
||||||
do ni=1,nint
|
do ni=1,nint
|
||||||
if(key(ni,1,i) /= key(ni,1,i-1)) then
|
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
|
||||||
end do
|
end do
|
||||||
shortcut(shortcut(0)+1) = N_key+1
|
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
|
end subroutine
|
||||||
|
|
||||||
c
|
|
||||||
|
|
||||||
subroutine sort_dets_ab(key, idx, shortcut, N_key, Nint)
|
subroutine sort_dets_ab(key, idx, shortcut, N_key, Nint)
|
||||||
use bitmasks
|
use bitmasks
|
||||||
implicit none
|
implicit none
|
||||||
|
|
||||||
|
|
||||||
|
BEGIN_DOC
|
||||||
|
! Uncodumented : TODO
|
||||||
|
END_DOC
|
||||||
integer(bit_kind),intent(inout) :: key(Nint,2,N_key)
|
integer(bit_kind),intent(inout) :: key(Nint,2,N_key)
|
||||||
integer,intent(out) :: idx(N_key)
|
integer,intent(out) :: idx(N_key)
|
||||||
integer,intent(out) :: shortcut(0:N_key+1)
|
integer,intent(out) :: shortcut(0:N_key+1)
|
||||||
|
@ -214,9 +255,15 @@ subroutine sort_dets_ab(key, idx, shortcut, N_key, Nint)
|
||||||
end do
|
end do
|
||||||
|
|
||||||
do i=N_key,2,-1
|
do i=N_key,2,-1
|
||||||
tmp(:,:) = key(:,:,i)
|
do ni=1,Nint
|
||||||
key(:,:,i) = key(:,:,1)
|
tmp(ni,1) = key(ni,1,i)
|
||||||
key(:,:,1) = tmp(:,:)
|
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)
|
tmpidx = idx(i)
|
||||||
idx(i) = idx(1)
|
idx(i) = idx(1)
|
||||||
idx(1) = tmpidx
|
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
|
abort_here = abort_all
|
||||||
end
|
end
|
||||||
|
|
||||||
BEGIN_PROVIDER [ character(64), davidson_criterion ]
|
BEGIN_PROVIDER [ character(64), davidson_criterion ]
|
||||||
&BEGIN_PROVIDER [ double precision, davidson_threshold ]
|
|
||||||
implicit none
|
implicit none
|
||||||
BEGIN_DOC
|
BEGIN_DOC
|
||||||
! Can be : [ energy | residual | both | wall_time | cpu_time | iterations ]
|
! Can be : [ energy | residual | both | wall_time | cpu_time | iterations ]
|
||||||
END_DOC
|
END_DOC
|
||||||
davidson_criterion = 'residual'
|
davidson_criterion = 'residual'
|
||||||
davidson_threshold = 1.d-10
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
subroutine davidson_converged(energy,residual,wall,iterations,cpu,N_st,converged)
|
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
|
E = energy - energy_old
|
||||||
energy_old = energy
|
energy_old = energy
|
||||||
if (davidson_criterion == 'energy') then
|
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
|
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
|
else if (davidson_criterion == 'both') then
|
||||||
converged = dabs(maxval(residual(1:N_st))) + dabs(maxval(E(1:N_st)) ) &
|
converged = dabs(maxval(residual(1:N_st))) + dabs(maxval(E(1:N_st)) ) &
|
||||||
< davidson_threshold
|
< threshold_davidson
|
||||||
else if (davidson_criterion == 'wall_time') then
|
else if (davidson_criterion == 'wall_time') then
|
||||||
call wall_time(time)
|
call wall_time(time)
|
||||||
converged = time - wall > davidson_threshold
|
converged = time - wall > threshold_davidson
|
||||||
else if (davidson_criterion == 'cpu_time') then
|
else if (davidson_criterion == 'cpu_time') then
|
||||||
call cpu_time(time)
|
call cpu_time(time)
|
||||||
converged = time - cpu > davidson_threshold
|
converged = time - cpu > threshold_davidson
|
||||||
else if (davidson_criterion == 'iterations') then
|
else if (davidson_criterion == 'iterations') then
|
||||||
converged = iterations >= int(davidson_threshold)
|
converged = iterations >= int(threshold_davidson)
|
||||||
endif
|
endif
|
||||||
converged = converged.or.abort_here
|
converged = converged.or.abort_here
|
||||||
end
|
end
|
||||||
|
|
|
@ -8,6 +8,7 @@ BEGIN_PROVIDER [ integer, N_det ]
|
||||||
logical :: exists
|
logical :: exists
|
||||||
character*64 :: label
|
character*64 :: label
|
||||||
PROVIDE ezfio_filename
|
PROVIDE ezfio_filename
|
||||||
|
PROVIDE nproc
|
||||||
if (read_wf) then
|
if (read_wf) then
|
||||||
call ezfio_has_determinants_n_det(exists)
|
call ezfio_has_determinants_n_det(exists)
|
||||||
if (exists) then
|
if (exists) then
|
||||||
|
|
162
src/Determinants/guess_lowest_state.irp.f
Normal file
162
src/Determinants/guess_lowest_state.irp.f
Normal file
|
@ -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
|
|
@ -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*,'<D_i|H|D_j> = ',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*,'<D_i|h(1)|D_j> = ',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*,'<ij|kl> = ',hij
|
|
||||||
|
|
||||||
|
|
||||||
end
|
|
|
@ -109,81 +109,81 @@ end
|
||||||
subroutine get_s2_u0(psi_keys_tmp,psi_coefs_tmp,n,nmax,s2)
|
subroutine get_s2_u0(psi_keys_tmp,psi_coefs_tmp,n,nmax,s2)
|
||||||
implicit none
|
implicit none
|
||||||
use bitmasks
|
use bitmasks
|
||||||
integer(bit_kind), intent(in) :: psi_keys_tmp(N_int,2,nmax)
|
integer(bit_kind), intent(in) :: psi_keys_tmp(N_int,2,nmax)
|
||||||
integer, intent(in) :: n,nmax
|
integer, intent(in) :: n,nmax
|
||||||
double precision, intent(in) :: psi_coefs_tmp(nmax)
|
double precision, intent(in) :: psi_coefs_tmp(nmax)
|
||||||
double precision, intent(out) :: s2
|
double precision, intent(out) :: s2
|
||||||
double precision :: s2_tmp
|
double precision :: s2_tmp
|
||||||
integer :: i,j,l,jj,ii
|
integer :: i,j,l,jj,ii
|
||||||
integer, allocatable :: idx(:)
|
integer, allocatable :: idx(:)
|
||||||
|
|
||||||
integer :: shortcut(0:n+1), sort_idx(n)
|
integer, allocatable :: shortcut(:), sort_idx(:)
|
||||||
integer(bit_kind) :: sorted(N_int,n), version(N_int,n)
|
integer(bit_kind), allocatable :: sorted(:,:), version(:,:)
|
||||||
integer :: sh, sh2, ni, exa, ext, org_i, org_j, endi, pass
|
integer :: sh, sh2, ni, exa, ext, org_i, org_j, endi, pass
|
||||||
double precision :: davidson_threshold_bis
|
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
|
s2 = 0.d0
|
||||||
davidson_threshold_bis = davidson_threshold
|
davidson_threshold_bis = threshold_davidson
|
||||||
!$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
|
|
||||||
call sort_dets_ab_v(psi_keys_tmp, sorted, sort_idx, shortcut, version, n, N_int)
|
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)
|
!$OMP DO SCHEDULE(dynamic)
|
||||||
do sh=1,shortcut(0)
|
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
|
do sh2=1,sh
|
||||||
if(sh==sh2) then
|
exa = 0
|
||||||
endi = i-1
|
do ni=1,N_int
|
||||||
else
|
exa += popcnt(xor(version(ni,sh), version(ni,sh2)))
|
||||||
endi = shortcut(sh2+1)-1
|
end do
|
||||||
|
if(exa > 2) then
|
||||||
|
cycle
|
||||||
end if
|
end if
|
||||||
|
|
||||||
do j=shortcut(sh2),endi
|
do i=shortcut(sh),shortcut(sh+1)-1
|
||||||
ext = exa
|
if(sh==sh2) then
|
||||||
do ni=1,N_int
|
endi = i-1
|
||||||
ext += popcnt(xor(sorted(ni,i), sorted(ni,j)))
|
else
|
||||||
end do
|
endi = shortcut(sh2+1)-1
|
||||||
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
|
|
||||||
end if
|
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
|
end do
|
||||||
end do
|
|
||||||
enddo
|
enddo
|
||||||
!$OMP END DO
|
!$OMP END DO
|
||||||
|
|
||||||
!$OMP SINGLE
|
!$OMP END PARALLEL
|
||||||
|
|
||||||
call sort_dets_ba_v(psi_keys_tmp, sorted, sort_idx, shortcut, version, n, N_int)
|
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)
|
!$OMP DO SCHEDULE(dynamic)
|
||||||
do sh=1,shortcut(0)
|
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
|
do j=shortcut(sh),i-1
|
||||||
ext = 0
|
ext = 0
|
||||||
do ni=1,N_int
|
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_i = sort_idx(i)
|
||||||
org_j = sort_idx(j)
|
org_j = sort_idx(j)
|
||||||
|
|
||||||
if ( dabs(psi_coefs_tmp(org_j)) + dabs(psi_coefs_tmp(org_i)) &
|
if ( dabs(psi_coefs_tmp(org_j)) + dabs(psi_coefs_tmp(org_i))&
|
||||||
> davidson_threshold ) then
|
> threshold_davidson ) then
|
||||||
call get_s2(psi_keys_tmp(1,1,org_i),psi_keys_tmp(1,1,org_j),s2_tmp,N_int)
|
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
|
s2 = s2 + psi_coefs_tmp(org_i)*psi_coefs_tmp(org_j)*s2_tmp
|
||||||
endif
|
endif
|
||||||
end if
|
end if
|
||||||
end do
|
end do
|
||||||
end do
|
end do
|
||||||
enddo
|
enddo
|
||||||
!$OMP END DO
|
!$OMP END DO
|
||||||
|
|
||||||
deallocate(idx)
|
!$OMP END PARALLEL
|
||||||
!$OMP END PARALLEL
|
s2 = s2+s2
|
||||||
s2 = s2+s2
|
do i=1,n
|
||||||
do i=1,n
|
call get_s2(psi_keys_tmp(1,1,i),psi_keys_tmp(1,1,i),s2_tmp,N_int)
|
||||||
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
|
||||||
s2 = s2 + psi_coefs_tmp(i)*psi_coefs_tmp(i)*s2_tmp
|
enddo
|
||||||
enddo
|
s2 = s2 + S_z2_Sz
|
||||||
s2 = s2 + S_z2_Sz
|
deallocate (shortcut, sort_idx, sorted, version)
|
||||||
end
|
end
|
||||||
|
|
||||||
|
|
||||||
|
|
|
@ -15,7 +15,7 @@ subroutine get_excitation_degree(key1,key2,degree,Nint)
|
||||||
|
|
||||||
degree = popcnt(xor( key1(1,1), key2(1,1))) + &
|
degree = popcnt(xor( key1(1,1), key2(1,1))) + &
|
||||||
popcnt(xor( key1(1,2), key2(1,2)))
|
popcnt(xor( key1(1,2), key2(1,2)))
|
||||||
!DEC$ NOUNROLL
|
!DIR$ NOUNROLL
|
||||||
do l=2,Nint
|
do l=2,Nint
|
||||||
degree = degree+ popcnt(xor( key1(l,1), key2(l,1))) + &
|
degree = degree+ popcnt(xor( key1(l,1), key2(l,1))) + &
|
||||||
popcnt(xor( key1(l,2), key2(l,2)))
|
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 :: exc(0:2,2,2)
|
||||||
integer :: degree
|
integer :: degree
|
||||||
double precision :: get_mo_bielec_integral
|
double precision :: get_mo_bielec_integral_schwartz
|
||||||
integer :: m,n,p,q
|
integer :: m,n,p,q
|
||||||
integer :: i,j,k
|
integer :: i,j,k
|
||||||
integer :: occ(Nint*bit_kind_size,2)
|
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)
|
ASSERT (sum(popcnt(key_j(:,2))) == elec_beta_num)
|
||||||
|
|
||||||
hij = 0.d0
|
hij = 0.d0
|
||||||
!DEC$ FORCEINLINE
|
!DIR$ FORCEINLINE
|
||||||
call get_excitation_degree(key_i,key_j,degree,Nint)
|
call get_excitation_degree(key_i,key_j,degree,Nint)
|
||||||
select case (degree)
|
select case (degree)
|
||||||
case (2)
|
case (2)
|
||||||
call get_double_excitation(key_i,key_j,exc,phase,Nint)
|
call get_double_excitation(key_i,key_j,exc,phase,Nint)
|
||||||
if (exc(0,1,1) == 1) then
|
if (exc(0,1,1) == 1) then
|
||||||
! Mono alpha, mono beta
|
! Mono alpha, mono beta
|
||||||
hij = phase*get_mo_bielec_integral( &
|
hij = phase*get_mo_bielec_integral_schwartz( &
|
||||||
exc(1,1,1), &
|
exc(1,1,1), &
|
||||||
exc(1,1,2), &
|
exc(1,1,2), &
|
||||||
exc(1,2,1), &
|
exc(1,2,1), &
|
||||||
exc(1,2,2) ,mo_integrals_map)
|
exc(1,2,2) ,mo_integrals_map)
|
||||||
else if (exc(0,1,1) == 2) then
|
else if (exc(0,1,1) == 2) then
|
||||||
! Double alpha
|
! Double alpha
|
||||||
hij = phase*(get_mo_bielec_integral( &
|
hij = phase*(get_mo_bielec_integral_schwartz( &
|
||||||
exc(1,1,1), &
|
exc(1,1,1), &
|
||||||
exc(2,1,1), &
|
exc(2,1,1), &
|
||||||
exc(1,2,1), &
|
exc(1,2,1), &
|
||||||
exc(2,2,1) ,mo_integrals_map) - &
|
exc(2,2,1) ,mo_integrals_map) - &
|
||||||
get_mo_bielec_integral( &
|
get_mo_bielec_integral_schwartz( &
|
||||||
exc(1,1,1), &
|
exc(1,1,1), &
|
||||||
exc(2,1,1), &
|
exc(2,1,1), &
|
||||||
exc(2,2,1), &
|
exc(2,2,1), &
|
||||||
exc(1,2,1) ,mo_integrals_map) )
|
exc(1,2,1) ,mo_integrals_map) )
|
||||||
else if (exc(0,1,2) == 2) then
|
else if (exc(0,1,2) == 2) then
|
||||||
! Double beta
|
! Double beta
|
||||||
hij = phase*(get_mo_bielec_integral( &
|
hij = phase*(get_mo_bielec_integral_schwartz( &
|
||||||
exc(1,1,2), &
|
exc(1,1,2), &
|
||||||
exc(2,1,2), &
|
exc(2,1,2), &
|
||||||
exc(1,2,2), &
|
exc(1,2,2), &
|
||||||
exc(2,2,2) ,mo_integrals_map) - &
|
exc(2,2,2) ,mo_integrals_map) - &
|
||||||
get_mo_bielec_integral( &
|
get_mo_bielec_integral_schwartz( &
|
||||||
exc(1,1,2), &
|
exc(1,1,2), &
|
||||||
exc(2,1,2), &
|
exc(2,1,2), &
|
||||||
exc(2,2,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
|
do k = 1, elec_alpha_num
|
||||||
i = occ(k,1)
|
i = occ(k,1)
|
||||||
if (.not.has_mipi(i)) then
|
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)
|
||||||
miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map)
|
miip(i) = get_mo_bielec_integral_schwartz(m,i,i,p,mo_integrals_map)
|
||||||
has_mipi(i) = .True.
|
has_mipi(i) = .True.
|
||||||
endif
|
endif
|
||||||
enddo
|
enddo
|
||||||
do k = 1, elec_beta_num
|
do k = 1, elec_beta_num
|
||||||
i = occ(k,2)
|
i = occ(k,2)
|
||||||
if (.not.has_mipi(i)) then
|
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.
|
has_mipi(i) = .True.
|
||||||
endif
|
endif
|
||||||
enddo
|
enddo
|
||||||
|
@ -459,15 +459,15 @@ subroutine i_H_j(key_i,key_j,Nint,hij)
|
||||||
do k = 1, elec_beta_num
|
do k = 1, elec_beta_num
|
||||||
i = occ(k,2)
|
i = occ(k,2)
|
||||||
if (.not.has_mipi(i)) then
|
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)
|
||||||
miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map)
|
miip(i) = get_mo_bielec_integral_schwartz(m,i,i,p,mo_integrals_map)
|
||||||
has_mipi(i) = .True.
|
has_mipi(i) = .True.
|
||||||
endif
|
endif
|
||||||
enddo
|
enddo
|
||||||
do k = 1, elec_alpha_num
|
do k = 1, elec_alpha_num
|
||||||
i = occ(k,1)
|
i = occ(k,1)
|
||||||
if (.not.has_mipi(i)) then
|
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.
|
has_mipi(i) = .True.
|
||||||
endif
|
endif
|
||||||
enddo
|
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) :: exc(0:2,2,2)
|
||||||
integer,intent(out) :: degree
|
integer,intent(out) :: degree
|
||||||
double precision :: get_mo_bielec_integral
|
double precision :: get_mo_bielec_integral_schwartz
|
||||||
integer :: m,n,p,q
|
integer :: m,n,p,q
|
||||||
integer :: i,j,k
|
integer :: i,j,k
|
||||||
integer :: occ(Nint*bit_kind_size,2)
|
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)
|
ASSERT (sum(popcnt(key_j(:,2))) == elec_beta_num)
|
||||||
|
|
||||||
hij = 0.d0
|
hij = 0.d0
|
||||||
!DEC$ FORCEINLINE
|
!DIR$ FORCEINLINE
|
||||||
call get_excitation_degree(key_i,key_j,degree,Nint)
|
call get_excitation_degree(key_i,key_j,degree,Nint)
|
||||||
select case (degree)
|
select case (degree)
|
||||||
case (2)
|
case (2)
|
||||||
call get_double_excitation(key_i,key_j,exc,phase,Nint)
|
call get_double_excitation(key_i,key_j,exc,phase,Nint)
|
||||||
if (exc(0,1,1) == 1) then
|
if (exc(0,1,1) == 1) then
|
||||||
! Mono alpha, mono beta
|
! Mono alpha, mono beta
|
||||||
hij = phase*get_mo_bielec_integral( &
|
hij = phase*get_mo_bielec_integral_schwartz( &
|
||||||
exc(1,1,1), &
|
exc(1,1,1), &
|
||||||
exc(1,1,2), &
|
exc(1,1,2), &
|
||||||
exc(1,2,1), &
|
exc(1,2,1), &
|
||||||
exc(1,2,2) ,mo_integrals_map)
|
exc(1,2,2) ,mo_integrals_map)
|
||||||
else if (exc(0,1,1) == 2) then
|
else if (exc(0,1,1) == 2) then
|
||||||
! Double alpha
|
! Double alpha
|
||||||
hij = phase*(get_mo_bielec_integral( &
|
hij = phase*(get_mo_bielec_integral_schwartz( &
|
||||||
exc(1,1,1), &
|
exc(1,1,1), &
|
||||||
exc(2,1,1), &
|
exc(2,1,1), &
|
||||||
exc(1,2,1), &
|
exc(1,2,1), &
|
||||||
exc(2,2,1) ,mo_integrals_map) - &
|
exc(2,2,1) ,mo_integrals_map) - &
|
||||||
get_mo_bielec_integral( &
|
get_mo_bielec_integral_schwartz( &
|
||||||
exc(1,1,1), &
|
exc(1,1,1), &
|
||||||
exc(2,1,1), &
|
exc(2,1,1), &
|
||||||
exc(2,2,1), &
|
exc(2,2,1), &
|
||||||
exc(1,2,1) ,mo_integrals_map) )
|
exc(1,2,1) ,mo_integrals_map) )
|
||||||
else if (exc(0,1,2) == 2) then
|
else if (exc(0,1,2) == 2) then
|
||||||
! Double beta
|
! Double beta
|
||||||
hij = phase*(get_mo_bielec_integral( &
|
hij = phase*(get_mo_bielec_integral_schwartz( &
|
||||||
exc(1,1,2), &
|
exc(1,1,2), &
|
||||||
exc(2,1,2), &
|
exc(2,1,2), &
|
||||||
exc(1,2,2), &
|
exc(1,2,2), &
|
||||||
exc(2,2,2) ,mo_integrals_map) - &
|
exc(2,2,2) ,mo_integrals_map) - &
|
||||||
get_mo_bielec_integral( &
|
get_mo_bielec_integral_schwartz( &
|
||||||
exc(1,1,2), &
|
exc(1,1,2), &
|
||||||
exc(2,1,2), &
|
exc(2,1,2), &
|
||||||
exc(2,2,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
|
do k = 1, elec_alpha_num
|
||||||
i = occ(k,1)
|
i = occ(k,1)
|
||||||
if (.not.has_mipi(i)) then
|
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)
|
||||||
miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map)
|
miip(i) = get_mo_bielec_integral_schwartz(m,i,i,p,mo_integrals_map)
|
||||||
has_mipi(i) = .True.
|
has_mipi(i) = .True.
|
||||||
endif
|
endif
|
||||||
enddo
|
enddo
|
||||||
do k = 1, elec_beta_num
|
do k = 1, elec_beta_num
|
||||||
i = occ(k,2)
|
i = occ(k,2)
|
||||||
if (.not.has_mipi(i)) then
|
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.
|
has_mipi(i) = .True.
|
||||||
endif
|
endif
|
||||||
enddo
|
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
|
do k = 1, elec_beta_num
|
||||||
i = occ(k,2)
|
i = occ(k,2)
|
||||||
if (.not.has_mipi(i)) then
|
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)
|
||||||
miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map)
|
miip(i) = get_mo_bielec_integral_schwartz(m,i,i,p,mo_integrals_map)
|
||||||
has_mipi(i) = .True.
|
has_mipi(i) = .True.
|
||||||
endif
|
endif
|
||||||
enddo
|
enddo
|
||||||
do k = 1, elec_alpha_num
|
do k = 1, elec_alpha_num
|
||||||
i = occ(k,1)
|
i = occ(k,1)
|
||||||
if (.not.has_mipi(i)) then
|
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.
|
has_mipi(i) = .True.
|
||||||
endif
|
endif
|
||||||
enddo
|
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 :: exc(0:2,2,2)
|
||||||
integer :: degree
|
integer :: degree
|
||||||
double precision :: get_mo_bielec_integral
|
double precision :: get_mo_bielec_integral_schwartz
|
||||||
integer :: m,n,p,q
|
integer :: m,n,p,q
|
||||||
integer :: i,j,k
|
integer :: i,j,k
|
||||||
integer :: occ(Nint*bit_kind_size,2)
|
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
|
hij = 0.d0
|
||||||
hmono = 0.d0
|
hmono = 0.d0
|
||||||
hdouble = 0.d0
|
hdouble = 0.d0
|
||||||
!DEC$ FORCEINLINE
|
!DIR$ FORCEINLINE
|
||||||
call get_excitation_degree(key_i,key_j,degree,Nint)
|
call get_excitation_degree(key_i,key_j,degree,Nint)
|
||||||
select case (degree)
|
select case (degree)
|
||||||
case (2)
|
case (2)
|
||||||
call get_double_excitation(key_i,key_j,exc,phase,Nint)
|
call get_double_excitation(key_i,key_j,exc,phase,Nint)
|
||||||
if (exc(0,1,1) == 1) then
|
if (exc(0,1,1) == 1) then
|
||||||
! Mono alpha, mono beta
|
! Mono alpha, mono beta
|
||||||
hij = phase*get_mo_bielec_integral( &
|
hij = phase*get_mo_bielec_integral_schwartz( &
|
||||||
exc(1,1,1), &
|
exc(1,1,1), &
|
||||||
exc(1,1,2), &
|
exc(1,1,2), &
|
||||||
exc(1,2,1), &
|
exc(1,2,1), &
|
||||||
exc(1,2,2) ,mo_integrals_map)
|
exc(1,2,2) ,mo_integrals_map)
|
||||||
else if (exc(0,1,1) == 2) then
|
else if (exc(0,1,1) == 2) then
|
||||||
! Double alpha
|
! Double alpha
|
||||||
hij = phase*(get_mo_bielec_integral( &
|
hij = phase*(get_mo_bielec_integral_schwartz( &
|
||||||
exc(1,1,1), &
|
exc(1,1,1), &
|
||||||
exc(2,1,1), &
|
exc(2,1,1), &
|
||||||
exc(1,2,1), &
|
exc(1,2,1), &
|
||||||
exc(2,2,1) ,mo_integrals_map) - &
|
exc(2,2,1) ,mo_integrals_map) - &
|
||||||
get_mo_bielec_integral( &
|
get_mo_bielec_integral_schwartz( &
|
||||||
exc(1,1,1), &
|
exc(1,1,1), &
|
||||||
exc(2,1,1), &
|
exc(2,1,1), &
|
||||||
exc(2,2,1), &
|
exc(2,2,1), &
|
||||||
exc(1,2,1) ,mo_integrals_map) )
|
exc(1,2,1) ,mo_integrals_map) )
|
||||||
else if (exc(0,1,2) == 2) then
|
else if (exc(0,1,2) == 2) then
|
||||||
! Double beta
|
! Double beta
|
||||||
hij = phase*(get_mo_bielec_integral( &
|
hij = phase*(get_mo_bielec_integral_schwartz( &
|
||||||
exc(1,1,2), &
|
exc(1,1,2), &
|
||||||
exc(2,1,2), &
|
exc(2,1,2), &
|
||||||
exc(1,2,2), &
|
exc(1,2,2), &
|
||||||
exc(2,2,2) ,mo_integrals_map) - &
|
exc(2,2,2) ,mo_integrals_map) - &
|
||||||
get_mo_bielec_integral( &
|
get_mo_bielec_integral_schwartz( &
|
||||||
exc(1,1,2), &
|
exc(1,1,2), &
|
||||||
exc(2,1,2), &
|
exc(2,1,2), &
|
||||||
exc(2,2,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
|
do k = 1, elec_alpha_num
|
||||||
i = occ(k,1)
|
i = occ(k,1)
|
||||||
if (.not.has_mipi(i)) then
|
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)
|
||||||
miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map)
|
miip(i) = get_mo_bielec_integral_schwartz(m,i,i,p,mo_integrals_map)
|
||||||
has_mipi(i) = .True.
|
has_mipi(i) = .True.
|
||||||
endif
|
endif
|
||||||
enddo
|
enddo
|
||||||
do k = 1, elec_beta_num
|
do k = 1, elec_beta_num
|
||||||
i = occ(k,2)
|
i = occ(k,2)
|
||||||
if (.not.has_mipi(i)) then
|
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.
|
has_mipi(i) = .True.
|
||||||
endif
|
endif
|
||||||
enddo
|
enddo
|
||||||
|
@ -733,15 +733,15 @@ subroutine i_H_j_verbose(key_i,key_j,Nint,hij,hmono,hdouble)
|
||||||
do k = 1, elec_beta_num
|
do k = 1, elec_beta_num
|
||||||
i = occ(k,2)
|
i = occ(k,2)
|
||||||
if (.not.has_mipi(i)) then
|
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)
|
||||||
miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map)
|
miip(i) = get_mo_bielec_integral_schwartz(m,i,i,p,mo_integrals_map)
|
||||||
has_mipi(i) = .True.
|
has_mipi(i) = .True.
|
||||||
endif
|
endif
|
||||||
enddo
|
enddo
|
||||||
do k = 1, elec_alpha_num
|
do k = 1, elec_alpha_num
|
||||||
i = occ(k,1)
|
i = occ(k,1)
|
||||||
if (.not.has_mipi(i)) then
|
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.
|
has_mipi(i) = .True.
|
||||||
endif
|
endif
|
||||||
enddo
|
enddo
|
||||||
|
@ -863,9 +863,17 @@ subroutine create_minilist_find_previous(key_mask, fullList, miniList, N_fullLis
|
||||||
end subroutine
|
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
|
use bitmasks
|
||||||
implicit none
|
implicit none
|
||||||
|
BEGIN_DOC
|
||||||
|
! Computes <i|H|Psi> = \sum_J c_J <i|H|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, intent(in) :: Nint, Ndet,Ndet_max,Nstate
|
||||||
integer(bit_kind), intent(in) :: keys(Nint,2,Ndet)
|
integer(bit_kind), intent(in) :: keys(Nint,2,Ndet)
|
||||||
integer(bit_kind), intent(in) :: key(Nint,2)
|
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)
|
integer :: exc(0:2,2,2)
|
||||||
double precision :: hij
|
double precision :: hij
|
||||||
integer :: idx(0:Ndet)
|
integer :: idx(0:Ndet)
|
||||||
BEGIN_DOC
|
|
||||||
! <key|H|psi> for the various Nstates
|
|
||||||
END_DOC
|
|
||||||
|
|
||||||
ASSERT (Nint > 0)
|
ASSERT (Nint > 0)
|
||||||
ASSERT (N_int == Nint)
|
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)
|
call filter_connected_i_H_psi0(keys,key,Nint,Ndet,idx)
|
||||||
do ii=1,idx(0)
|
do ii=1,idx(0)
|
||||||
i = idx(ii)
|
i = idx(ii)
|
||||||
!DEC$ FORCEINLINE
|
!DIR$ FORCEINLINE
|
||||||
call i_H_j(keys(1,1,i),key,Nint,hij)
|
call i_H_j(keys(1,1,i),key,Nint,hij)
|
||||||
do j = 1, Nstate
|
do j = 1, Nstate
|
||||||
i_H_psi_array(j) = i_H_psi_array(j) + coef(i,j)*hij
|
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
|
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
|
use bitmasks
|
||||||
implicit none
|
implicit none
|
||||||
integer, intent(in) :: Nint, Ndet,Ndet_max,Nstate,idx_key(Ndet), N_minilist
|
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
|
double precision :: hij
|
||||||
integer :: idx(0:Ndet)
|
integer :: idx(0:Ndet)
|
||||||
BEGIN_DOC
|
BEGIN_DOC
|
||||||
! <key|H|psi> for the various Nstates
|
! Computes <i|H|Psi> = \sum_J c_J <i|H|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
|
END_DOC
|
||||||
|
|
||||||
ASSERT (Nint > 0)
|
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)
|
ASSERT (Ndet_max >= Ndet)
|
||||||
i_H_psi_array = 0.d0
|
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)
|
call filter_connected_i_H_psi0(keys,key,Nint,N_minilist,idx)
|
||||||
do ii=1,idx(0)
|
do ii=1,idx(0)
|
||||||
!i = idx_key(idx(ii))
|
|
||||||
i_in_key = idx(ii)
|
i_in_key = idx(ii)
|
||||||
i_in_coef = idx_key(idx(ii))
|
i_in_coef = idx_key(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
|
|
||||||
! ! enddo
|
|
||||||
call i_H_j(keys(1,1,i_in_key),key,Nint,hij)
|
call i_H_j(keys(1,1,i_in_key),key,Nint,hij)
|
||||||
do j = 1, Nstate
|
do j = 1, Nstate
|
||||||
i_H_psi_array(j) = i_H_psi_array(j) + coef(i_in_coef,j)*hij
|
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
|
n_interact = 0
|
||||||
do ii=1,idx(0)
|
do ii=1,idx(0)
|
||||||
i = idx(ii)
|
i = idx(ii)
|
||||||
!DEC$ FORCEINLINE
|
!DIR$ FORCEINLINE
|
||||||
call i_H_j(keys(1,1,i),key,Nint,hij)
|
call i_H_j(keys(1,1,i),key,Nint,hij)
|
||||||
if(dabs(hij).ge.1.d-8)then
|
if(dabs(hij).ge.1.d-8)then
|
||||||
if(i.ne.1)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)
|
call filter_connected_i_H_psi0_SC2(keys,key,Nint,Ndet,idx,idx_repeat)
|
||||||
do ii=1,idx(0)
|
do ii=1,idx(0)
|
||||||
i = idx(ii)
|
i = idx(ii)
|
||||||
!DEC$ FORCEINLINE
|
!DIR$ FORCEINLINE
|
||||||
call i_H_j(keys(1,1,i),key,Nint,hij)
|
call i_H_j(keys(1,1,i),key,Nint,hij)
|
||||||
do j = 1, Nstate
|
do j = 1, Nstate
|
||||||
i_H_psi_array(j) = i_H_psi_array(j) + coef(i,j)*hij
|
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)
|
do ii=1,idx(0)
|
||||||
print*,'--'
|
print*,'--'
|
||||||
i = idx(ii)
|
i = idx(ii)
|
||||||
!DEC$ FORCEINLINE
|
!DIR$ FORCEINLINE
|
||||||
call i_H_j(keys(1,1,i),key,Nint,hij)
|
call i_H_j(keys(1,1,i),key,Nint,hij)
|
||||||
if (i==1)then
|
if (i==1)then
|
||||||
print*,'i==1 !!'
|
print*,'i==1 !!'
|
||||||
|
@ -1167,7 +1169,7 @@ subroutine get_excitation_degree_vector(key1,key2,degree,Nint,sze,idx)
|
||||||
!DIR$ LOOP COUNT (1000)
|
!DIR$ LOOP COUNT (1000)
|
||||||
do i=1,sze
|
do i=1,sze
|
||||||
d = 0
|
d = 0
|
||||||
!DEC$ LOOP COUNT MIN(4)
|
!DIR$ LOOP COUNT MIN(4)
|
||||||
do m=1,Nint
|
do m=1,Nint
|
||||||
d = d + popcnt(xor( key1(m,1,i), key2(m,1))) &
|
d = d + popcnt(xor( key1(m,1,i), key2(m,1))) &
|
||||||
+ popcnt(xor( key1(m,2,i), key2(m,2)))
|
+ 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(1) = 0
|
||||||
nexc(2) = 0
|
nexc(2) = 0
|
||||||
do i=1,Nint
|
do i=1,Nint
|
||||||
hole(i,1) = xor(det_in(i,1),ref_bitmask(i,1))
|
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,2) = xor(det_in(i,2),ref_bitmask(i,2))
|
||||||
particle(i,1) = iand(hole(i,1),det_in(i,1))
|
particle(i,1) = iand(hole(i,1),det_in(i,1))
|
||||||
particle(i,2) = iand(hole(i,2),det_in(i,2))
|
particle(i,2) = iand(hole(i,2),det_in(i,2))
|
||||||
hole(i,1) = iand(hole(i,1),ref_bitmask(i,1))
|
hole(i,1) = iand(hole(i,1),ref_bitmask(i,1))
|
||||||
hole(i,2) = iand(hole(i,2),ref_bitmask(i,2))
|
hole(i,2) = iand(hole(i,2),ref_bitmask(i,2))
|
||||||
nexc(1) += popcnt(hole(i,1))
|
nexc(1) = nexc(1) + popcnt(hole(i,1))
|
||||||
nexc(2) += popcnt(hole(i,2))
|
nexc(2) = nexc(2) + popcnt(hole(i,2))
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
diag_H_mat_elem = ref_bitmask_energy
|
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 :: i,j,k,l, jj,ii
|
||||||
integer :: i0, j0
|
integer :: i0, j0
|
||||||
|
|
||||||
integer :: shortcut(0:n+1), sort_idx(n)
|
integer, allocatable :: shortcut(:), sort_idx(:)
|
||||||
integer(bit_kind) :: sorted(Nint,n), version(Nint,n)
|
integer(bit_kind), allocatable :: sorted(:,:), version(:,:)
|
||||||
|
integer(bit_kind) :: sorted_i(Nint)
|
||||||
|
|
||||||
integer :: sh, sh2, ni, exa, ext, org_i, org_j, endi
|
integer :: sh, sh2, ni, exa, ext, org_i, org_j, endi
|
||||||
!
|
double precision :: local_threshold
|
||||||
|
|
||||||
|
|
||||||
ASSERT (Nint > 0)
|
ASSERT (Nint > 0)
|
||||||
ASSERT (Nint == N_int)
|
ASSERT (Nint == N_int)
|
||||||
ASSERT (n>0)
|
ASSERT (n>0)
|
||||||
PROVIDE ref_bitmask_energy
|
PROVIDE ref_bitmask_energy davidson_criterion
|
||||||
!$OMP PARALLEL DEFAULT(NONE) &
|
|
||||||
!$OMP PRIVATE(i,hij,j,k,idx,jj,vt,ii,sh,sh2,ni,exa,ext,org_i,org_j,endi) &
|
allocate (shortcut(0:n+1), sort_idx(n), sorted(Nint,n), version(Nint,n))
|
||||||
!$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
|
|
||||||
v_0 = 0.d0
|
v_0 = 0.d0
|
||||||
|
|
||||||
!$OMP SINGLE
|
|
||||||
call sort_dets_ab_v(keys_tmp, sorted, sort_idx, shortcut, version, n, Nint)
|
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)
|
!$OMP DO SCHEDULE(dynamic)
|
||||||
do sh=1,shortcut(0)
|
do sh=1,shortcut(0)
|
||||||
do sh2=1,sh
|
do sh2=1,sh
|
||||||
exa = 0
|
exa = 0
|
||||||
do ni=1,Nint
|
do ni=1,Nint
|
||||||
exa += popcnt(xor(version(ni,sh), version(ni,sh2)))
|
exa = exa + popcnt(xor(version(ni,sh), version(ni,sh2)))
|
||||||
end do
|
end do
|
||||||
if(exa > 2) then
|
if(exa > 2) then
|
||||||
cycle
|
cycle
|
||||||
end if
|
|
||||||
|
|
||||||
do i=shortcut(sh),shortcut(sh+1)-1
|
|
||||||
if(sh==sh2) then
|
|
||||||
endi = i-1
|
|
||||||
else
|
|
||||||
endi = shortcut(sh2+1)-1
|
|
||||||
end if
|
end if
|
||||||
|
|
||||||
do j=shortcut(sh2),endi
|
do i=shortcut(sh),shortcut(sh+1)-1
|
||||||
ext = exa
|
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
|
do ni=1,Nint
|
||||||
ext += popcnt(xor(sorted(ni,i), sorted(ni,j)))
|
sorted_i(ni) = sorted(ni,i)
|
||||||
end do
|
enddo
|
||||||
if(ext <= 4) then
|
|
||||||
org_i = sort_idx(i)
|
do j=shortcut(sh2),endi
|
||||||
org_j = sort_idx(j)
|
org_j = sort_idx(j)
|
||||||
if ( dabs(u_0(org_j)) + dabs(u_0(org_i)) > davidson_threshold ) then
|
if ( dabs(u_0(org_j)) > local_threshold ) then
|
||||||
call i_H_j(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),Nint,hij)
|
ext = exa
|
||||||
vt (org_i) = vt (org_i) + hij*u_0(org_j)
|
do ni=1,Nint
|
||||||
vt (org_j) = vt (org_j) + hij*u_0(org_i)
|
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
|
||||||
endif
|
enddo
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
enddo
|
|
||||||
!$OMP END DO
|
!$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)
|
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)
|
!$OMP DO SCHEDULE(dynamic)
|
||||||
do sh=1,shortcut(0)
|
do sh=1,shortcut(0)
|
||||||
do i=shortcut(sh),shortcut(sh+1)-1
|
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
|
do j=shortcut(sh),i-1
|
||||||
ext = 0
|
org_j = sort_idx(j)
|
||||||
do ni=1,Nint
|
if ( dabs(u_0(org_j)) > local_threshold ) then
|
||||||
ext += popcnt(xor(sorted(ni,i), sorted(ni,j)))
|
ext = 0
|
||||||
end do
|
do ni=1,Nint
|
||||||
if(ext == 4) then
|
ext = ext + popcnt(xor(sorted(ni,i), sorted(ni,j)))
|
||||||
org_i = sort_idx(i)
|
end do
|
||||||
org_j = sort_idx(j)
|
if(ext == 4) then
|
||||||
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)
|
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_i) = vt (org_i) + hij*u_0(org_j)
|
||||||
vt (org_j) = vt (org_j) + hij*u_0(org_i)
|
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)
|
v_0(i) = v_0(i) + vt(i)
|
||||||
enddo
|
enddo
|
||||||
!$OMP END CRITICAL
|
!$OMP END CRITICAL
|
||||||
deallocate(idx,vt)
|
deallocate(vt)
|
||||||
!$OMP END PARALLEL
|
!$OMP END PARALLEL
|
||||||
|
|
||||||
do i=1,n
|
do i=1,n
|
||||||
v_0(i) += H_jj(i) * u_0(i)
|
v_0(i) += H_jj(i) * u_0(i)
|
||||||
enddo
|
enddo
|
||||||
|
deallocate (shortcut, sort_idx, sorted, version)
|
||||||
end
|
end
|
||||||
|
|
||||||
|
|
|
@ -442,13 +442,14 @@ BEGIN_PROVIDER [ double precision, psi_bilinear_matrix, (N_det_alpha_unique,N_de
|
||||||
enddo
|
enddo
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
subroutine create_wf_of_psi_bilinear_matrix
|
subroutine create_wf_of_psi_bilinear_matrix(truncate)
|
||||||
use bitmasks
|
use bitmasks
|
||||||
implicit none
|
implicit none
|
||||||
BEGIN_DOC
|
BEGIN_DOC
|
||||||
! Generate a wave function containing all possible products
|
! Generate a wave function containing all possible products
|
||||||
! of alpha and beta determinants
|
! of alpha and beta determinants
|
||||||
END_DOC
|
END_DOC
|
||||||
|
logical, intent(in) :: truncate
|
||||||
integer :: i,j,k
|
integer :: i,j,k
|
||||||
integer(bit_kind) :: tmp_det(N_int,2)
|
integer(bit_kind) :: tmp_det(N_int,2)
|
||||||
integer :: idx
|
integer :: idx
|
||||||
|
@ -488,8 +489,10 @@ subroutine create_wf_of_psi_bilinear_matrix
|
||||||
norm(1) = 0.d0
|
norm(1) = 0.d0
|
||||||
do i=1,N_det
|
do i=1,N_det
|
||||||
norm(1) += psi_average_norm_contrib_sorted(i)
|
norm(1) += psi_average_norm_contrib_sorted(i)
|
||||||
if (norm(1) >= 0.999999d0) then
|
if (truncate) then
|
||||||
exit
|
if (norm(1) >= 0.999999d0) then
|
||||||
|
exit
|
||||||
|
endif
|
||||||
endif
|
endif
|
||||||
enddo
|
enddo
|
||||||
N_det = min(i,N_det)
|
N_det = min(i,N_det)
|
||||||
|
@ -532,7 +535,6 @@ subroutine generate_all_alpha_beta_det_products
|
||||||
!$OMP END DO NOWAIT
|
!$OMP END DO NOWAIT
|
||||||
deallocate(tmp_det)
|
deallocate(tmp_det)
|
||||||
!$OMP END PARALLEL
|
!$OMP END PARALLEL
|
||||||
deallocate (tmp_det)
|
|
||||||
call copy_H_apply_buffer_to_wf
|
call copy_H_apply_buffer_to_wf
|
||||||
SOFT_TOUCH psi_det psi_coef N_det
|
SOFT_TOUCH psi_det psi_coef N_det
|
||||||
end
|
end
|
||||||
|
|
|
@ -8,10 +8,10 @@ program cisd
|
||||||
N_det=10000
|
N_det=10000
|
||||||
do i=1,N_det
|
do i=1,N_det
|
||||||
do k=1,N_int
|
do k=1,N_int
|
||||||
psi_det(k,1,i) = psi_det_sorted(k,1,i)
|
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,2,i) = psi_det_sorted(k,2,i)
|
||||||
enddo
|
enddo
|
||||||
psi_coef(k,:) = psi_coef_sorted(k,:)
|
psi_coef(i,:) = psi_coef_sorted(i,:)
|
||||||
enddo
|
enddo
|
||||||
TOUCH psi_det psi_coef psi_det_sorted psi_coef_sorted psi_average_norm_contrib_sorted N_det
|
TOUCH psi_det psi_coef psi_det_sorted psi_coef_sorted psi_average_norm_contrib_sorted N_det
|
||||||
call save_wavefunction
|
call save_wavefunction
|
||||||
|
|
|
@ -204,7 +204,7 @@ double precision function ao_bielec_integral_schwartz_accel(i,j,k,l)
|
||||||
integral = general_primitive_integral(dim1, &
|
integral = general_primitive_integral(dim1, &
|
||||||
P_new,P_center,fact_p,pp,p_inv,iorder_p, &
|
P_new,P_center,fact_p,pp,p_inv,iorder_p, &
|
||||||
Q_new,Q_center,fact_q,qq,q_inv,iorder_q)
|
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 ! s
|
||||||
enddo ! r
|
enddo ! r
|
||||||
enddo ! q
|
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(1),J_power(1),K_power(1),L_power(1), &
|
||||||
I_power(2),J_power(2),K_power(2),L_power(2), &
|
I_power(2),J_power(2),K_power(2),L_power(2), &
|
||||||
I_power(3),J_power(3),K_power(3),L_power(3))
|
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 ! s
|
||||||
enddo ! r
|
enddo ! r
|
||||||
enddo ! q
|
enddo ! q
|
||||||
|
@ -307,12 +307,20 @@ subroutine compute_ao_bielec_integrals(j,k,l,sze,buffer_value)
|
||||||
buffer_value = 0._integral_kind
|
buffer_value = 0._integral_kind
|
||||||
return
|
return
|
||||||
endif
|
endif
|
||||||
|
if (ao_bielec_integral_schwartz(j,l) < thresh ) then
|
||||||
|
buffer_value = 0._integral_kind
|
||||||
|
return
|
||||||
|
endif
|
||||||
|
|
||||||
do i = 1, ao_num
|
do i = 1, ao_num
|
||||||
if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < thresh) then
|
if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < thresh) then
|
||||||
buffer_value(i) = 0._integral_kind
|
buffer_value(i) = 0._integral_kind
|
||||||
cycle
|
cycle
|
||||||
endif
|
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
|
!DIR$ FORCEINLINE
|
||||||
buffer_value(i) = ao_bielec_integral(i,k,j,l)
|
buffer_value(i) = ao_bielec_integral(i,k,j,l)
|
||||||
enddo
|
enddo
|
||||||
|
@ -378,8 +386,9 @@ BEGIN_PROVIDER [ logical, ao_bielec_integrals_in_map ]
|
||||||
!$OMP DEFAULT(NONE) &
|
!$OMP DEFAULT(NONE) &
|
||||||
!$OMP SHARED (ao_num, jl_pairs, ao_integrals_map,thresh, &
|
!$OMP SHARED (ao_num, jl_pairs, ao_integrals_map,thresh, &
|
||||||
!$OMP cpu_1,wall_1,lock, lmax,n_centers,ao_nucl, &
|
!$OMP cpu_1,wall_1,lock, lmax,n_centers,ao_nucl, &
|
||||||
!$OMP ao_overlap_abs,ao_overlap,abort_here, &
|
!$OMP ao_overlap_abs,ao_overlap,abort_here, &
|
||||||
!$OMP wall_0,progress_bar,progress_value)
|
!$OMP wall_0,progress_bar,progress_value, &
|
||||||
|
!$OMP ao_bielec_integral_schwartz)
|
||||||
|
|
||||||
allocate(buffer_i(size_buffer))
|
allocate(buffer_i(size_buffer))
|
||||||
allocate(buffer_value(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
|
if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < thresh) then
|
||||||
cycle
|
cycle
|
||||||
endif
|
endif
|
||||||
|
if (ao_bielec_integral_schwartz(i,k)*ao_bielec_integral_schwartz(j,l) < thresh ) then
|
||||||
|
cycle
|
||||||
|
endif
|
||||||
!DIR$ FORCEINLINE
|
!DIR$ FORCEINLINE
|
||||||
integral = ao_bielec_integral(i,k,j,l)
|
integral = ao_bielec_integral(i,k,j,l)
|
||||||
if (abs(integral) < thresh) then
|
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
|
integer :: n_pt_sup
|
||||||
double precision :: p,q,denom,coeff
|
double precision :: p,q,denom,coeff
|
||||||
double precision :: I_f
|
double precision :: I_f
|
||||||
|
integer :: nx,ny,nz
|
||||||
include 'Utils/constants.include.F'
|
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
|
ERI = 0.d0
|
||||||
return
|
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
|
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
|
end
|
||||||
|
|
||||||
|
|
||||||
|
|
|
@ -291,19 +291,42 @@ double precision function get_mo_bielec_integral(i,j,k,l,map)
|
||||||
PROVIDE mo_bielec_integrals_in_map
|
PROVIDE mo_bielec_integrals_in_map
|
||||||
!DIR$ FORCEINLINE
|
!DIR$ FORCEINLINE
|
||||||
call bielec_integrals_index(i,j,k,l,idx)
|
call bielec_integrals_index(i,j,k,l,idx)
|
||||||
|
!DIR$ FORCEINLINE
|
||||||
call map_get(map,idx,tmp)
|
call map_get(map,idx,tmp)
|
||||||
get_mo_bielec_integral = dble(tmp)
|
get_mo_bielec_integral = dble(tmp)
|
||||||
end
|
end
|
||||||
|
|
||||||
|
double precision function get_mo_bielec_integral_schwartz(i,j,k,l,map)
|
||||||
|
use map_module
|
||||||
|
implicit none
|
||||||
|
BEGIN_DOC
|
||||||
|
! Returns one integral <ij|kl> 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)
|
double precision function mo_bielec_integral(i,j,k,l)
|
||||||
implicit none
|
implicit none
|
||||||
BEGIN_DOC
|
BEGIN_DOC
|
||||||
! Returns one integral <ij|kl> in the MO basis
|
! Returns one integral <ij|kl> in the MO basis
|
||||||
END_DOC
|
END_DOC
|
||||||
integer, intent(in) :: i,j,k,l
|
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
|
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
|
return
|
||||||
end
|
end
|
||||||
|
|
||||||
|
@ -502,7 +525,8 @@ integer function load_$ao_integrals(filename)
|
||||||
integer*8 :: i
|
integer*8 :: i
|
||||||
integer(cache_key_kind), pointer :: key(:)
|
integer(cache_key_kind), pointer :: key(:)
|
||||||
real(integral_kind), pointer :: val(:)
|
real(integral_kind), pointer :: val(:)
|
||||||
integer :: iknd, kknd, n, j
|
integer :: iknd, kknd
|
||||||
|
integer*8 :: n, j
|
||||||
load_$ao_integrals = 1
|
load_$ao_integrals = 1
|
||||||
open(unit=66,file=filename,FORM='unformatted',STATUS='UNKNOWN')
|
open(unit=66,file=filename,FORM='unformatted',STATUS='UNKNOWN')
|
||||||
read(66,err=98,end=98) iknd, kknd
|
read(66,err=98,end=98) iknd, kknd
|
||||||
|
@ -532,12 +556,8 @@ integer function load_$ao_integrals(filename)
|
||||||
return
|
return
|
||||||
99 continue
|
99 continue
|
||||||
call map_deinit($ao_integrals_map)
|
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
|
98 continue
|
||||||
|
stop 'Problem reading $ao_integrals_map file in work/'
|
||||||
|
|
||||||
end
|
end
|
||||||
|
|
||||||
|
|
|
@ -488,3 +488,19 @@ END_PROVIDER
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
END_PROVIDER
|
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
|
||||||
|
|
136
src/Integrals_Bielec/test_integrals.irp.f
Normal file
136
src/Integrals_Bielec/test_integrals.irp.f
Normal file
|
@ -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
|
|
@ -5,106 +5,105 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integral, (ao_num_align,ao_num)]
|
||||||
END_DOC
|
END_DOC
|
||||||
if (do_pseudo) then
|
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_non_local
|
||||||
! ao_pseudo_integral = ao_pseudo_integral_local
|
|
||||||
! ao_pseudo_integral = ao_pseudo_integral_non_local
|
|
||||||
else
|
else
|
||||||
ao_pseudo_integral = 0.d0
|
ao_pseudo_integral = 0.d0
|
||||||
endif
|
endif
|
||||||
END_PROVIDER
|
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
|
implicit none
|
||||||
BEGIN_DOC
|
BEGIN_DOC
|
||||||
! Local pseudo-potential
|
! Local pseudo-potential
|
||||||
END_DOC
|
END_DOC
|
||||||
double precision :: alpha, beta, gama, delta
|
double precision :: alpha, beta, gama, delta
|
||||||
integer :: num_A,num_B
|
integer :: num_A,num_B
|
||||||
double precision :: A_center(3),B_center(3),C_center(3)
|
double precision :: A_center(3),B_center(3),C_center(3)
|
||||||
integer :: power_A(3),power_B(3)
|
integer :: power_A(3),power_B(3)
|
||||||
integer :: i,j,k,l,n_pt_in,m
|
integer :: i,j,k,l,n_pt_in,m
|
||||||
double precision :: Vloc, Vpseudo
|
double precision :: Vloc, Vpseudo
|
||||||
|
|
||||||
double precision :: cpu_1, cpu_2, wall_1, wall_2, wall_0
|
double precision :: cpu_1, cpu_2, wall_1, wall_2, wall_0
|
||||||
integer :: thread_num
|
integer :: thread_num
|
||||||
|
!$ integer :: omp_get_thread_num
|
||||||
|
|
||||||
ao_pseudo_integral_local = 0.d0
|
ao_pseudo_integral_local = 0.d0
|
||||||
|
|
||||||
!! Dump array
|
!! Dump array
|
||||||
integer, allocatable :: n_k_dump(:)
|
integer, allocatable :: n_k_dump(:)
|
||||||
double precision, allocatable :: v_k_dump(:), dz_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))
|
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)'
|
||||||
! _
|
|
||||||
! / _. | _ |
|
|
||||||
! \_ (_| | (_ |_| |
|
|
||||||
!
|
|
||||||
|
|
||||||
print*, 'Providing the nuclear electron pseudo integrals '
|
|
||||||
|
|
||||||
call wall_time(wall_1)
|
call wall_time(wall_1)
|
||||||
call cpu_time(cpu_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)
|
!$OMP DO SCHEDULE (guided)
|
||||||
|
|
||||||
do j = 1, ao_num
|
do j = 1, ao_num
|
||||||
|
|
||||||
num_A = ao_nucl(j)
|
num_A = ao_nucl(j)
|
||||||
power_A(1:3)= ao_power(j,1:3)
|
power_A(1:3)= ao_power(j,1:3)
|
||||||
A_center(1:3) = nucl_coord(num_A,1:3)
|
A_center(1:3) = nucl_coord(num_A,1:3)
|
||||||
|
|
||||||
do i = 1, ao_num
|
do i = 1, ao_num
|
||||||
|
|
||||||
num_B = ao_nucl(i)
|
num_B = ao_nucl(i)
|
||||||
power_B(1:3)= ao_power(i,1:3)
|
power_B(1:3)= ao_power(i,1:3)
|
||||||
B_center(1:3) = nucl_coord(num_B,1:3)
|
B_center(1:3) = nucl_coord(num_B,1:3)
|
||||||
|
|
||||||
do l=1,ao_prim_num(j)
|
do l=1,ao_prim_num(j)
|
||||||
alpha = ao_expo_ordered_transp(l,j)
|
alpha = ao_expo_ordered_transp(l,j)
|
||||||
|
|
||||||
do m=1,ao_prim_num(i)
|
do m=1,ao_prim_num(i)
|
||||||
beta = ao_expo_ordered_transp(m,i)
|
beta = ao_expo_ordered_transp(m,i)
|
||||||
double precision :: c
|
double precision :: c
|
||||||
c = 0.d0
|
c = 0.d0
|
||||||
|
|
||||||
do k = 1, nucl_num
|
if (dabs(ao_coef_normalized_ordered_transp(l,j)*ao_coef_normalized_ordered_transp(m,i))&
|
||||||
double precision :: Z
|
< 1.d-10) then
|
||||||
Z = nucl_charge(k)
|
cycle
|
||||||
|
endif
|
||||||
C_center(1:3) = nucl_coord(k,1:3)
|
do k = 1, nucl_num
|
||||||
|
double precision :: Z
|
||||||
v_k_dump = pseudo_v_k(k,1:pseudo_klocmax)
|
Z = nucl_charge(k)
|
||||||
n_k_dump = pseudo_n_k(k,1:pseudo_klocmax)
|
|
||||||
dz_k_dump = pseudo_dz_k(k,1:pseudo_klocmax)
|
C_center(1:3) = nucl_coord(k,1:3)
|
||||||
|
|
||||||
c = c + Vloc(pseudo_klocmax, v_k_dump,n_k_dump, dz_k_dump, &
|
v_k_dump = pseudo_v_k(k,1:pseudo_klocmax)
|
||||||
A_center,power_A,alpha,B_center,power_B,beta,C_center)
|
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
|
enddo
|
||||||
ao_pseudo_integral_local(i,j) = ao_pseudo_integral_local(i,j) + &
|
enddo
|
||||||
ao_coef_normalized_ordered_transp(l,j)*ao_coef_normalized_ordered_transp(m,i)*c
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
call wall_time(wall_2)
|
call wall_time(wall_2)
|
||||||
if (thread_num == 0) then
|
if (thread_num == 0) then
|
||||||
if (wall_2 - wall_0 > 1.d0) then
|
if (wall_2 - wall_0 > 1.d0) then
|
||||||
wall_0 = wall_2
|
wall_0 = wall_2
|
||||||
print*, 100.*float(j)/float(ao_num), '% in ', &
|
print*, 100.*float(j)/float(ao_num), '% in ', &
|
||||||
wall_2-wall_1, 's'
|
wall_2-wall_1, 's'
|
||||||
endif
|
endif
|
||||||
endif
|
endif
|
||||||
enddo
|
enddo
|
||||||
|
@ -121,106 +120,108 @@ END_PROVIDER
|
||||||
BEGIN_PROVIDER [ double precision, ao_pseudo_integral_non_local, (ao_num_align,ao_num)]
|
BEGIN_PROVIDER [ double precision, ao_pseudo_integral_non_local, (ao_num_align,ao_num)]
|
||||||
implicit none
|
implicit none
|
||||||
BEGIN_DOC
|
BEGIN_DOC
|
||||||
! Local pseudo-potential
|
! Local pseudo-potential
|
||||||
END_DOC
|
END_DOC
|
||||||
double precision :: alpha, beta, gama, delta
|
double precision :: alpha, beta, gama, delta
|
||||||
integer :: num_A,num_B
|
integer :: num_A,num_B
|
||||||
double precision :: A_center(3),B_center(3),C_center(3)
|
double precision :: A_center(3),B_center(3),C_center(3)
|
||||||
integer :: power_A(3),power_B(3)
|
integer :: power_A(3),power_B(3)
|
||||||
integer :: i,j,k,l,n_pt_in,m
|
integer :: i,j,k,l,n_pt_in,m
|
||||||
double precision :: Vloc, Vpseudo
|
double precision :: Vloc, Vpseudo
|
||||||
|
!$ integer :: omp_get_thread_num
|
||||||
double precision :: cpu_1, cpu_2, wall_1, wall_2, wall_0
|
|
||||||
integer :: thread_num
|
double precision :: cpu_1, cpu_2, wall_1, wall_2, wall_0
|
||||||
|
integer :: thread_num
|
||||||
|
|
||||||
ao_pseudo_integral_non_local = 0.d0
|
ao_pseudo_integral_non_local = 0.d0
|
||||||
|
|
||||||
!! Dump array
|
!! Dump array
|
||||||
integer, allocatable :: n_kl_dump(:,:)
|
integer, allocatable :: n_kl_dump(:,:)
|
||||||
double precision, allocatable :: v_kl_dump(:,:), dz_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))
|
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)'
|
||||||
! / _. | _ |
|
|
||||||
! \_ (_| | (_ |_| |
|
|
||||||
!
|
|
||||||
|
|
||||||
print*, 'Providing the nuclear electron pseudo integrals '
|
|
||||||
|
|
||||||
call wall_time(wall_1)
|
call wall_time(wall_1)
|
||||||
call cpu_time(cpu_1)
|
call cpu_time(cpu_1)
|
||||||
|
thread_num = 0
|
||||||
!$OMP PARALLEL &
|
!$OMP PARALLEL &
|
||||||
!$OMP DEFAULT (NONE) &
|
!$OMP DEFAULT (NONE) &
|
||||||
!$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,C_center,power_A,power_B, &
|
!$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 num_A,num_B,Z,c,n_pt_in, &
|
||||||
!$OMP n_kl_dump, v_kl_dump, dz_kl_dump, &
|
!$OMP n_kl_dump, v_kl_dump, dz_kl_dump, &
|
||||||
!$OMP wall_0,wall_2,thread_num) &
|
!$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 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 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 pseudo_klocmax,pseudo_lmax,pseudo_kmax,pseudo_n_kl, pseudo_v_kl, pseudo_dz_kl,&
|
||||||
!$OMP wall_1)
|
!$OMP wall_1)
|
||||||
|
|
||||||
|
!$ thread_num = omp_get_thread_num()
|
||||||
!$OMP DO SCHEDULE (guided)
|
!$OMP DO SCHEDULE (guided)
|
||||||
|
|
||||||
do j = 1, ao_num
|
do j = 1, ao_num
|
||||||
|
|
||||||
num_A = ao_nucl(j)
|
num_A = ao_nucl(j)
|
||||||
power_A(1:3)= ao_power(j,1:3)
|
power_A(1:3)= ao_power(j,1:3)
|
||||||
A_center(1:3) = nucl_coord(num_A,1:3)
|
A_center(1:3) = nucl_coord(num_A,1:3)
|
||||||
|
|
||||||
do i = 1, ao_num
|
do i = 1, ao_num
|
||||||
|
|
||||||
num_B = ao_nucl(i)
|
num_B = ao_nucl(i)
|
||||||
power_B(1:3)= ao_power(i,1:3)
|
power_B(1:3)= ao_power(i,1:3)
|
||||||
B_center(1:3) = nucl_coord(num_B,1:3)
|
B_center(1:3) = nucl_coord(num_B,1:3)
|
||||||
|
|
||||||
do l=1,ao_prim_num(j)
|
do l=1,ao_prim_num(j)
|
||||||
alpha = ao_expo_ordered_transp(l,j)
|
alpha = ao_expo_ordered_transp(l,j)
|
||||||
|
|
||||||
do m=1,ao_prim_num(i)
|
do m=1,ao_prim_num(i)
|
||||||
beta = ao_expo_ordered_transp(m,i)
|
beta = ao_expo_ordered_transp(m,i)
|
||||||
double precision :: c
|
double precision :: c
|
||||||
c = 0.d0
|
c = 0.d0
|
||||||
|
|
||||||
do k = 1, nucl_num
|
if (dabs(ao_coef_normalized_ordered_transp(l,j)*ao_coef_normalized_ordered_transp(m,i))&
|
||||||
double precision :: Z
|
< 1.d-10) then
|
||||||
Z = nucl_charge(k)
|
cycle
|
||||||
|
endif
|
||||||
C_center(1:3) = nucl_coord(k,1:3)
|
|
||||||
|
do k = 1, nucl_num
|
||||||
n_kl_dump = pseudo_n_kl(k,1:pseudo_kmax,0:pseudo_lmax)
|
double precision :: Z
|
||||||
v_kl_dump = pseudo_v_kl(k,1:pseudo_kmax,0:pseudo_lmax)
|
Z = nucl_charge(k)
|
||||||
dz_kl_dump = pseudo_dz_kl(k,1:pseudo_kmax,0:pseudo_lmax)
|
|
||||||
|
C_center(1:3) = nucl_coord(k,1:3)
|
||||||
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)
|
|
||||||
|
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
|
enddo
|
||||||
ao_pseudo_integral_non_local(i,j) = ao_pseudo_integral_non_local(i,j) + &
|
enddo
|
||||||
ao_coef_normalized_ordered_transp(l,j)*ao_coef_normalized_ordered_transp(m,i)*c
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
call wall_time(wall_2)
|
call wall_time(wall_2)
|
||||||
if (thread_num == 0) then
|
if (thread_num == 0) then
|
||||||
if (wall_2 - wall_0 > 1.d0) then
|
if (wall_2 - wall_0 > 1.d0) then
|
||||||
wall_0 = wall_2
|
wall_0 = wall_2
|
||||||
print*, 100.*float(j)/float(ao_num), '% in ', &
|
print*, 100.*float(j)/float(ao_num), '% in ', &
|
||||||
wall_2-wall_1, 's'
|
wall_2-wall_1, 's'
|
||||||
endif
|
endif
|
||||||
endif
|
endif
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
!$OMP END DO
|
!$OMP END DO
|
||||||
!$OMP END PARALLEL
|
!$OMP END PARALLEL
|
||||||
|
|
||||||
|
|
||||||
deallocate(n_kl_dump,v_kl_dump, dz_kl_dump)
|
deallocate(n_kl_dump,v_kl_dump, dz_kl_dump)
|
||||||
|
|
||||||
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
|
@ -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, intent(in) :: v_kl(kmax,0:lmax),dz_kl(kmax,0:lmax)
|
||||||
|
|
||||||
!
|
!
|
||||||
! | _ _ _. | _
|
! | _ _ _. |
|
||||||
! |_ (_) (_ (_| | (/_
|
! |_ (_) (_ (_| |
|
||||||
!
|
!
|
||||||
|
|
||||||
double precision :: fourpi,f,prod,prodp,binom_func,accu,bigR,bigI,ylm
|
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
|
double precision :: f1, f2, f3
|
||||||
|
|
||||||
! _
|
|
||||||
! / _. | _ |
|
|
||||||
! \_ (_| | (_ |_| |
|
|
||||||
!
|
|
||||||
|
|
||||||
if (kmax.eq.1.and.lmax.eq.0.and.v_kl(1,0).eq.0.d0) then
|
if (kmax.eq.1.and.lmax.eq.0.and.v_kl(1,0).eq.0.d0) then
|
||||||
Vpseudo=0.d0
|
Vpseudo=0.d0
|
||||||
return
|
return
|
||||||
|
@ -236,7 +231,7 @@ end if
|
||||||
fourpi=4.d0*dacos(-1.d0)
|
fourpi=4.d0*dacos(-1.d0)
|
||||||
ac=dsqrt((a(1)-c(1))**2+(a(2)-c(2))**2+(a(3)-c(3))**2)
|
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)
|
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
|
if(arg.gt.-dlog(1.d-20))then
|
||||||
Vpseudo=0.d0
|
Vpseudo=0.d0
|
||||||
|
@ -290,6 +285,21 @@ if(ac.eq.0.d0.and.bc.eq.0.d0)then
|
||||||
enddo
|
enddo
|
||||||
enddo
|
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 !
|
! 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
|
double precision, intent(in) :: rmax
|
||||||
|
|
||||||
!
|
!
|
||||||
! | _ _ _. | _
|
! | _ _ _. |
|
||||||
! |_ (_) (_ (_| | (/_
|
! |_ (_) (_ (_| |
|
||||||
!
|
!
|
||||||
|
|
||||||
integer :: l,m,k,kk
|
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 :: s_q_0, s_q_k, s_0_0, a_over_b_square
|
||||||
double precision :: int_prod_bessel_loc
|
double precision :: int_prod_bessel_loc
|
||||||
double precision :: inverses(0:300)
|
double precision :: inverses(0:300)
|
||||||
|
double precision :: two_qkmp1, qk
|
||||||
|
|
||||||
logical done
|
logical done
|
||||||
|
|
||||||
|
@ -2008,19 +2019,18 @@ double precision function int_prod_bessel(l,gam,n,m,a,b,arg)
|
||||||
stop 'pseudopot.f90 : q > 300'
|
stop 'pseudopot.f90 : q > 300'
|
||||||
endif
|
endif
|
||||||
|
|
||||||
|
two_qkmp1 = dble(2*(q+m)+1)
|
||||||
|
qk = dble(q)
|
||||||
do k=0,q-1
|
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
|
sum=sum+s_q_k
|
||||||
|
two_qkmp1 = two_qkmp1-2.d0
|
||||||
|
qk = qk-1.d0
|
||||||
enddo
|
enddo
|
||||||
inverses(q) = a_over_b_square/(dble(2*(q+n)+3) * dble(q+1))
|
inverses(q) = a_over_b_square/(dble(2*(q+n)+3) * dble(q+1))
|
||||||
! do k=0,q
|
! do k=0,q
|
||||||
! sum=sum+s_q_k
|
! 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
|
! 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
|
! enddo
|
||||||
|
|
||||||
int=int+sum
|
int=int+sum
|
||||||
|
|
Loading…
Reference in New Issue
Block a user