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

Linear dependencies cutoff

Browse Source
This commit is contained in:
Anthony Scemama 2020-05-25 11:31:28 +02:00
parent 9c52a612dd
commit 75891f14b7
22 changed files with 181 additions and 162 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
ocaml/Generic_input_of_rst.ml
View File

@ -1,6 +1,7 @@
open Sexplib open Sexplib
open Sexplib.Std open Sexplib.Std
open Qptypes open Qptypes
open Qputils
let fail_msg str (ex,range) = let fail_msg str (ex,range) =
@ -25,7 +26,7 @@ let fail_msg str (ex,range) =
in in
let str = String_ext.tr str ~target:'(' ~replacement:' ' let str = String_ext.tr str ~target:'(' ~replacement:' '
|> String_ext.split ~on:')' |> String_ext.split ~on:')'
|> List.map String_ext.strip |> list_map String_ext.strip
|> List.filter (fun x -> |> List.filter (fun x ->
match String_ext.substr_index ~pos:0 ~pattern:"##" x with match String_ext.substr_index ~pos:0 ~pattern:"##" x with
| None -> false | None -> false
@ -48,7 +49,7 @@ let of_rst t_of_sexp s =
Rst_string.to_string s Rst_string.to_string s
|> String_ext.split ~on:'\n' |> String_ext.split ~on:'\n'
|> List.filter (fun line -> String.contains line '=') |> List.filter (fun line -> String.contains line '=')
|> List.map (fun line -> |> list_map (fun line ->
"("^( "("^(
String_ext.tr ~target:'=' ~replacement:' ' line String_ext.tr ~target:'=' ~replacement:' ' line
)^")" ) )^")" )

46
ocaml/Input_ao_basis.ml
View File

@ -202,14 +202,14 @@ end = struct
in in
let ao_prim_num = let ao_prim_num =
Array.to_list ao_prim_num Array.to_list ao_prim_num
|> List.map AO_prim_number.to_int |> list_map AO_prim_number.to_int
in in
Ezfio.set_ao_basis_ao_prim_num (Ezfio.ezfio_array_of_list Ezfio.set_ao_basis_ao_prim_num (Ezfio.ezfio_array_of_list
~rank:1 ~dim:[| ao_num |] ~data:ao_prim_num) ; ~rank:1 ~dim:[| ao_num |] ~data:ao_prim_num) ;
let ao_nucl = let ao_nucl =
Array.to_list ao_nucl Array.to_list ao_nucl
|> List.map Nucl_number.to_int |> list_map Nucl_number.to_int
in in
Ezfio.set_ao_basis_ao_nucl(Ezfio.ezfio_array_of_list Ezfio.set_ao_basis_ao_nucl(Ezfio.ezfio_array_of_list
~rank:1 ~dim:[| ao_num |] ~data:ao_nucl) ; ~rank:1 ~dim:[| ao_num |] ~data:ao_nucl) ;
@ -217,9 +217,9 @@ end = struct
let ao_power = let ao_power =
let l = Array.to_list ao_power in let l = Array.to_list ao_power in
List.concat [ List.concat [
(List.map (fun a -> Positive_int.to_int a.Symmetry.Xyz.x) l) ; (list_map (fun a -> Positive_int.to_int a.Symmetry.Xyz.x) l) ;
(List.map (fun a -> Positive_int.to_int a.Symmetry.Xyz.y) l) ; (list_map (fun a -> Positive_int.to_int a.Symmetry.Xyz.y) l) ;
(List.map (fun a -> Positive_int.to_int a.Symmetry.Xyz.z) l) ] (list_map (fun a -> Positive_int.to_int a.Symmetry.Xyz.z) l) ]
in in
Ezfio.set_ao_basis_ao_power(Ezfio.ezfio_array_of_list Ezfio.set_ao_basis_ao_power(Ezfio.ezfio_array_of_list
~rank:2 ~dim:[| ao_num ; 3 |] ~data:ao_power) ; ~rank:2 ~dim:[| ao_num ; 3 |] ~data:ao_power) ;
@ -230,14 +230,14 @@ end = struct
let ao_coef = let ao_coef =
Array.to_list ao_coef Array.to_list ao_coef
|> List.map AO_coef.to_float |> list_map AO_coef.to_float
in in
Ezfio.set_ao_basis_ao_coef(Ezfio.ezfio_array_of_list Ezfio.set_ao_basis_ao_coef(Ezfio.ezfio_array_of_list
~rank:2 ~dim:[| ao_num ; ao_prim_num_max |] ~data:ao_coef) ; ~rank:2 ~dim:[| ao_num ; ao_prim_num_max |] ~data:ao_coef) ;
let ao_expo = let ao_expo =
Array.to_list ao_expo Array.to_list ao_expo
|> List.map AO_expo.to_float |> list_map AO_expo.to_float
in in
Ezfio.set_ao_basis_ao_expo(Ezfio.ezfio_array_of_list Ezfio.set_ao_basis_ao_expo(Ezfio.ezfio_array_of_list
~rank:2 ~dim:[| ao_num ; ao_prim_num_max |] ~data:ao_expo) ; ~rank:2 ~dim:[| ao_num ; ao_prim_num_max |] ~data:ao_expo) ;
@ -299,8 +299,8 @@ end = struct
| Some (s', g', n') -> | Some (s', g', n') ->
if s <> s' || n <> n' then find2 (s,g,n) a (i+1) if s <> s' || n <> n' then find2 (s,g,n) a (i+1)
else else
let lc = List.map (fun (prim, _) -> prim) g.Gto.lc let lc = list_map (fun (prim, _) -> prim) g.Gto.lc
and lc' = List.map (fun (prim, _) -> prim) g'.Gto.lc and lc' = list_map (fun (prim, _) -> prim) g'.Gto.lc
in in
if lc <> lc' then find2 (s,g,n) a (i+1) else (a.(i) <- None ; i) if lc <> lc' then find2 (s,g,n) a (i+1) else (a.(i) <- None ; i)
in in
@ -314,14 +314,14 @@ end = struct
let of_long_basis long_basis name ao_cartesian = let of_long_basis long_basis name ao_cartesian =
let ao_num = List.length long_basis |> AO_number.of_int in let ao_num = List.length long_basis |> AO_number.of_int in
let ao_prim_num = let ao_prim_num =
List.map (fun (_,g,_) -> List.length g.Gto.lc list_map (fun (_,g,_) -> List.length g.Gto.lc
|> AO_prim_number.of_int ) long_basis |> AO_prim_number.of_int ) long_basis
|> Array.of_list |> Array.of_list
and ao_nucl = and ao_nucl =
List.map (fun (_,_,n) -> n) long_basis list_map (fun (_,_,n) -> n) long_basis
|> Array.of_list |> Array.of_list
and ao_power = and ao_power =
List.map (fun (x,_,_) -> x) long_basis list_map (fun (x,_,_) -> x) long_basis
|> Array.of_list |> Array.of_list
in in
let ao_prim_num_max = Array.fold_left (fun s x -> let ao_prim_num_max = Array.fold_left (fun s x ->
@ -331,15 +331,15 @@ end = struct
in in
let gtos = let gtos =
List.map (fun (_,x,_) -> x) long_basis list_map (fun (_,x,_) -> x) long_basis
in in
let create_expo_coef ec = let create_expo_coef ec =
let coefs = let coefs =
begin match ec with begin match ec with
| `Coefs -> List.map (fun x-> | `Coefs -> list_map (fun x->
List.map (fun (_,coef) -> AO_coef.to_float coef) x.Gto.lc ) gtos list_map (fun (_,coef) -> AO_coef.to_float coef) x.Gto.lc ) gtos
| `Expos -> List.map (fun x-> | `Expos -> list_map (fun x->
List.map (fun (prim,_) -> AO_expo.to_float list_map (fun (prim,_) -> AO_expo.to_float
prim.GaussianPrimitive.expo) x.Gto.lc ) gtos prim.GaussianPrimitive.expo) x.Gto.lc ) gtos
end end
in in
@ -450,7 +450,7 @@ Basis set (read-only) ::
(string_of_bool b.ao_normalized) (string_of_bool b.ao_normalized)
(Basis.to_string short_basis (Basis.to_string short_basis
|> String_ext.split ~on:'\n' |> String_ext.split ~on:'\n'
|> List.map (fun x-> " "^x) |> list_map (fun x-> " "^x)
|> String.concat "\n" |> String.concat "\n"
) print_sym ) print_sym
@ -490,16 +490,16 @@ md5 = %s
" "
(AO_basis_name.to_string b.ao_basis) (AO_basis_name.to_string b.ao_basis)
(AO_number.to_string b.ao_num) (AO_number.to_string b.ao_num)
(b.ao_prim_num |> Array.to_list |> List.map (b.ao_prim_num |> Array.to_list |> list_map
(AO_prim_number.to_string) |> String.concat ", " ) (AO_prim_number.to_string) |> String.concat ", " )
(AO_prim_number.to_string b.ao_prim_num_max) (AO_prim_number.to_string b.ao_prim_num_max)
(b.ao_nucl |> Array.to_list |> List.map Nucl_number.to_string |> (b.ao_nucl |> Array.to_list |> list_map Nucl_number.to_string |>
String.concat ", ") String.concat ", ")
(b.ao_power |> Array.to_list |> List.map (fun x-> (b.ao_power |> Array.to_list |> list_map (fun x->
"("^(Symmetry.Xyz.to_string x)^")" )|> String.concat ", ") "("^(Symmetry.Xyz.to_string x)^")" )|> String.concat ", ")
(b.ao_coef |> Array.to_list |> List.map AO_coef.to_string (b.ao_coef |> Array.to_list |> list_map AO_coef.to_string
|> String.concat ", ") |> String.concat ", ")
(b.ao_expo |> Array.to_list |> List.map AO_expo.to_string (b.ao_expo |> Array.to_list |> list_map AO_expo.to_string
|> String.concat ", ") |> String.concat ", ")
(b.ao_cartesian |> string_of_bool) (b.ao_cartesian |> string_of_bool)
(b.ao_normalized |> string_of_bool) (b.ao_normalized |> string_of_bool)

12
ocaml/Input_determinants_by_hand.ml
View File

@ -377,7 +377,7 @@ end = struct
(coefs_string i) (coefs_string i)
(Determinant.to_string ~mo_num:mo_num b.psi_det.(i) (Determinant.to_string ~mo_num:mo_num b.psi_det.(i)
|> String_ext.split ~on:'\n' |> String_ext.split ~on:'\n'
|> List.map (fun x -> " "^x) |> list_map (fun x -> " "^x)
|> String.concat "\n" |> String.concat "\n"
) )
) )
@ -427,7 +427,7 @@ psi_det = %s
(b.n_det |> Det_number.to_string) (b.n_det |> Det_number.to_string)
(b.n_states |> States_number.to_string) (b.n_states |> States_number.to_string)
(b.expected_s2 |> Positive_float.to_string) (b.expected_s2 |> Positive_float.to_string)
(b.state_average_weight |> Array.to_list |> List.map Positive_float.to_string |> String.concat ",") (b.state_average_weight |> Array.to_list |> list_map Positive_float.to_string |> String.concat ",")
(b.psi_coef |> Array.map Det_coef.to_string |> Array.to_list (b.psi_coef |> Array.map Det_coef.to_string |> Array.to_list
|> String.concat ", ") |> String.concat ", ")
(b.psi_det |> Array.map (Determinant.to_string ~mo_num) |> Array.to_list (b.psi_det |> Array.map (Determinant.to_string ~mo_num) |> Array.to_list
@ -457,7 +457,7 @@ psi_det = %s
else else
( (String.contains line '=') && (line.[0] = ' ') ) ( (String.contains line '=') && (line.[0] = ' ') )
) )
|> List.map (fun line -> |> list_map (fun line ->
"("^( "("^(
String_ext.tr line ~target:'=' ~replacement:' ' String_ext.tr line ~target:'=' ~replacement:' '
|> String.trim |> String.trim
@ -468,7 +468,7 @@ psi_det = %s
(* Handle determinant coefs *) (* Handle determinant coefs *)
let dets = match ( dets let dets = match ( dets
|> String_ext.split ~on:'\n' |> String_ext.split ~on:'\n'
|> List.map String.trim |> list_map String.trim
) with ) with
| _::lines -> lines | _::lines -> lines
| _ -> failwith "Error in determinants" | _ -> failwith "Error in determinants"
@ -481,7 +481,7 @@ psi_det = %s
| ""::c::tail -> | ""::c::tail ->
let c = let c =
String_ext.split ~on:'\t' c String_ext.split ~on:'\t' c
|> List.map (fun x -> Det_coef.of_float (Float.of_string x)) |> list_map (fun x -> Det_coef.of_float (Float.of_string x))
|> Array.of_list |> Array.of_list
in in
read_coefs (c::accu) tail read_coefs (c::accu) tail
@ -499,7 +499,7 @@ psi_det = %s
let i = let i =
i-1 i-1
in in
List.map (fun x -> Det_coef.to_string x.(i)) buffer list_map (fun x -> Det_coef.to_string x.(i)) buffer
|> String.concat " " |> String.concat " "
in in
let rec build_result = function let rec build_result = function

4
ocaml/Input_mo_basis.ml
View File

@ -257,9 +257,9 @@ mo_coef = %s
" "
(MO_label.to_string b.mo_label) (MO_label.to_string b.mo_label)
(MO_number.to_string b.mo_num) (MO_number.to_string b.mo_num)
(b.mo_class |> Array.to_list |> List.map (b.mo_class |> Array.to_list |> list_map
(MO_class.to_string) |> String.concat ", " ) (MO_class.to_string) |> String.concat ", " )
(b.mo_occ |> Array.to_list |> List.map (b.mo_occ |> Array.to_list |> list_map
(MO_occ.to_string) |> String.concat ", " ) (MO_occ.to_string) |> String.concat ", " )
(b.mo_coef |> Array.map (b.mo_coef |> Array.map
(fun x-> Array.map MO_coef.to_string x |> (fun x-> Array.map MO_coef.to_string x |>

22
ocaml/Input_nuclei_by_hand.ml
View File

@ -50,7 +50,7 @@ end = struct
in in
let labels = let labels =
Array.to_list labels Array.to_list labels
|> List.map Element.to_string |> list_map Element.to_string
in in
Ezfio.ezfio_array_of_list ~rank:1 Ezfio.ezfio_array_of_list ~rank:1
~dim:[| nucl_num |] ~data:labels ~dim:[| nucl_num |] ~data:labels
@ -70,7 +70,7 @@ end = struct
in in
let charges = let charges =
Array.to_list charges Array.to_list charges
|> List.map Charge.to_float |> list_map Charge.to_float
in in
Ezfio.ezfio_array_of_list ~rank:1 Ezfio.ezfio_array_of_list ~rank:1
~dim:[| nucl_num |] ~data:charges ~dim:[| nucl_num |] ~data:charges
@ -101,9 +101,9 @@ end = struct
in in
let coord = Array.to_list coord in let coord = Array.to_list coord in
let coord = let coord =
(List.map (fun x-> x.Point3d.x) coord) @ (list_map (fun x-> x.Point3d.x) coord) @
(List.map (fun x-> x.Point3d.y) coord) @ (list_map (fun x-> x.Point3d.y) coord) @
(List.map (fun x-> x.Point3d.z) coord) (list_map (fun x-> x.Point3d.z) coord)
in in
Ezfio.ezfio_array_of_list ~rank:2 Ezfio.ezfio_array_of_list ~rank:2
~dim:[| nucl_num ; 3 |] ~data:coord ~dim:[| nucl_num ; 3 |] ~data:coord
@ -159,11 +159,11 @@ nucl_charge = %s
nucl_coord = %s nucl_coord = %s
" "
(Nucl_number.to_string b.nucl_num) (Nucl_number.to_string b.nucl_num)
(b.nucl_label |> Array.to_list |> List.map (b.nucl_label |> Array.to_list |> list_map
(Element.to_string) |> String.concat ", " ) (Element.to_string) |> String.concat ", " )
(b.nucl_charge |> Array.to_list |> List.map (b.nucl_charge |> Array.to_list |> list_map
(Charge.to_string) |> String.concat ", " ) (Charge.to_string) |> String.concat ", " )
(b.nucl_coord |> Array.to_list |> List.map (b.nucl_coord |> Array.to_list |> list_map
(Point3d.to_string ~units:Units.Bohr) |> String.concat "\n" ) (Point3d.to_string ~units:Units.Bohr) |> String.concat "\n" )
;; ;;
@ -226,11 +226,11 @@ Nuclear coordinates in xyz format (Angstroms) ::
let result = let result =
{ nucl_num = List.length atom_list { nucl_num = List.length atom_list
|> Nucl_number.of_int ~max:nmax; |> Nucl_number.of_int ~max:nmax;
nucl_label = List.map (fun x -> nucl_label = list_map (fun x ->
x.Atom.element) atom_list |> Array.of_list ; x.Atom.element) atom_list |> Array.of_list ;
nucl_charge = List.map (fun x -> nucl_charge = list_map (fun x ->
x.Atom.charge ) atom_list |> Array.of_list ; x.Atom.charge ) atom_list |> Array.of_list ;
nucl_coord = List.map (fun x -> nucl_coord = list_map (fun x ->
x.Atom.coord ) atom_list |> Array.of_list ; x.Atom.coord ) atom_list |> Array.of_list ;
} }
in Some result in Some result

3
ocaml/Long_basis.ml
View File

@ -1,4 +1,5 @@
open Qptypes open Qptypes
open Qputils
open Sexplib.Std open Sexplib.Std
type t = (Symmetry.Xyz.t * Gto.t * Nucl_number.t ) list [@@deriving sexp] type t = (Symmetry.Xyz.t * Gto.t * Nucl_number.t ) list [@@deriving sexp]
@ -39,7 +40,7 @@ let to_basis b =
let to_string b = let to_string b =
let middle = List.map (fun (x,y,z) -> let middle = list_map (fun (x,y,z) ->
"( "^((string_of_int (Nucl_number.to_int z)))^", "^ "( "^((string_of_int (Nucl_number.to_int z)))^", "^
(Symmetry.Xyz.to_string x)^", "^(Gto.to_string y) (Symmetry.Xyz.to_string x)^", "^(Gto.to_string y)
^" )" ^" )"

5
ocaml/MO_class.ml
View File

@ -1,4 +1,5 @@
open Qptypes open Qptypes
open Qputils
open Sexplib.Std open Sexplib.Std
@ -13,7 +14,7 @@ type t =
let to_string x = let to_string x =
let print_list l = let print_list l =
let s = List.map (fun x-> MO_number.to_int x |> string_of_int ) l let s = list_map (fun x-> MO_number.to_int x |> string_of_int ) l
|> (String.concat ", ") |> (String.concat ", ")
in in
"("^s^")" "("^s^")"
@ -43,7 +44,7 @@ let of_string s =
let _mo_number_list_of_range range = let _mo_number_list_of_range range =
Range.of_string range |> List.map MO_number.of_int Range.of_string range |> list_map MO_number.of_int
let create_core range = Core (_mo_number_list_of_range range) let create_core range = Core (_mo_number_list_of_range range)

13
ocaml/Message.ml
View File

@ -1,5 +1,6 @@
open Sexplib.Std open Sexplib.Std
open Qptypes open Qptypes
open Qputils
(** New job : Request to create a new multi-tasked job *) (** New job : Request to create a new multi-tasked job *)
@ -193,12 +194,12 @@ end = struct
} }
let create ~state ~task_ids = let create ~state ~task_ids =
{ state = State.of_string state ; { state = State.of_string state ;
task_ids = List.map Id.Task.of_int task_ids task_ids = list_map Id.Task.of_int task_ids
} }
let to_string x = let to_string x =
Printf.sprintf "del_task %s %s" Printf.sprintf "del_task %s %s"
(State.to_string x.state) (State.to_string x.state)
(String.concat "|" @@ List.map Id.Task.to_string x.task_ids) (String.concat "|" @@ list_map Id.Task.to_string x.task_ids)
end end
@ -219,7 +220,7 @@ end = struct
else "done" else "done"
in in
Printf.sprintf "del_task_reply %s %s" Printf.sprintf "del_task_reply %s %s"
more (String.concat "|" @@ List.map Id.Task.to_string x.task_ids) more (String.concat "|" @@ list_map Id.Task.to_string x.task_ids)
end end
@ -303,7 +304,7 @@ end = struct
"get_tasks_reply ok" "get_tasks_reply ok"
let to_string_list x = let to_string_list x =
"get_tasks_reply ok" :: ( "get_tasks_reply ok" :: (
List.map (fun (task_id, task) -> list_map (fun (task_id, task) ->
match task_id with match task_id with
| Some task_id -> Printf.sprintf "%d %s" (Id.Task.to_int task_id) task | Some task_id -> Printf.sprintf "%d %s" (Id.Task.to_int task_id) task
| None -> Printf.sprintf "0 terminate" | None -> Printf.sprintf "0 terminate"
@ -408,14 +409,14 @@ end = struct
let create ~state ~client_id ~task_ids = let create ~state ~client_id ~task_ids =
{ client_id = Id.Client.of_int client_id ; { client_id = Id.Client.of_int client_id ;
state = State.of_string state ; state = State.of_string state ;
task_ids = List.map Id.Task.of_int task_ids; task_ids = list_map Id.Task.of_int task_ids;
} }
let to_string x = let to_string x =
Printf.sprintf "task_done %s %d %s" Printf.sprintf "task_done %s %d %s"
(State.to_string x.state) (State.to_string x.state)
(Id.Client.to_int x.client_id) (Id.Client.to_int x.client_id)
(String.concat "|" @@ List.map Id.Task.to_string x.task_ids) (String.concat "|" @@ list_map Id.Task.to_string x.task_ids)
end end
(** Terminate *) (** Terminate *)

7
ocaml/Molecule.ml
View File

@ -1,4 +1,5 @@
open Qptypes open Qptypes
open Qputils
open Sexplib.Std open Sexplib.Std
exception MultiplicityError of string exception MultiplicityError of string
@ -96,7 +97,7 @@ let to_string_general ~f m =
let title = let title =
name m name m
in in
[ string_of_int n ; title ] @ (List.map f nuclei) [ string_of_int n ; title ] @ (list_map f nuclei)
|> String.concat "\n" |> String.concat "\n"
let to_string = let to_string =
@ -112,7 +113,7 @@ let of_xyz_string
s = s =
let l = String_ext.split s ~on:'\n' let l = String_ext.split s ~on:'\n'
|> List.filter (fun x -> x <> "") |> List.filter (fun x -> x <> "")
|> List.map (fun x -> Atom.of_string units x) |> list_map (fun x -> Atom.of_string units x)
in in
let ne = ( get_charge { let ne = ( get_charge {
nuclei=l ; nuclei=l ;
@ -186,7 +187,7 @@ let of_file
let distance_matrix molecule = let distance_matrix molecule =
let coord = let coord =
molecule.nuclei molecule.nuclei
|> List.map (fun x -> x.Atom.coord) |> list_map (fun x -> x.Atom.coord)
|> Array.of_list |> Array.of_list
in in
let n = let n =

3
ocaml/Point3d.ml
View File

@ -1,4 +1,5 @@
open Qptypes open Qptypes
open Qputils
open Sexplib.Std open Sexplib.Std
type t = { type t = {
@ -23,7 +24,7 @@ let of_string ~units s =
let l = s let l = s
|> String_ext.split ~on:' ' |> String_ext.split ~on:' '
|> List.filter (fun x -> x <> "") |> List.filter (fun x -> x <> "")
|> List.map float_of_string |> list_map float_of_string
|> Array.of_list |> Array.of_list
in in
{ x = l.(0) *. f ; { x = l.(0) *. f ;

5
ocaml/Pseudo.ml
View File

@ -1,4 +1,5 @@
open Sexplib.Std open Sexplib.Std
open Qputils
open Qptypes open Qptypes
@ -81,7 +82,7 @@ let to_string_local = function
| t -> | t ->
"Local component:" :: "Local component:" ::
( Printf.sprintf "%20s %8s %20s" "Coeff." "r^n" "Exp." ) :: ( Printf.sprintf "%20s %8s %20s" "Coeff." "r^n" "Exp." ) ::
( List.map (fun (l,c) -> Printf.sprintf "%20f %8d %20f" ( list_map (fun (l,c) -> Printf.sprintf "%20f %8d %20f"
(AO_coef.to_float c) (AO_coef.to_float c)
(R_power.to_int l.GaussianPrimitive_local.r_power) (R_power.to_int l.GaussianPrimitive_local.r_power)
(AO_expo.to_float l.GaussianPrimitive_local.expo) (AO_expo.to_float l.GaussianPrimitive_local.expo)
@ -95,7 +96,7 @@ let to_string_non_local = function
| t -> | t ->
"Non-local component:" :: "Non-local component:" ::
( Printf.sprintf "%20s %8s %20s %8s" "Coeff." "r^n" "Exp." "Proj") :: ( Printf.sprintf "%20s %8s %20s %8s" "Coeff." "r^n" "Exp." "Proj") ::
( List.map (fun (l,c) -> ( list_map (fun (l,c) ->
let p = let p =
Positive_int.to_int l.GaussianPrimitive_non_local.proj Positive_int.to_int l.GaussianPrimitive_non_local.proj
in in

4
ocaml/Qpackage.ml
View File

@ -30,7 +30,7 @@ let bit_kind_size = lazy (
in in
begin match (String_ext.rsplit2 ~on:':' line) with begin match (String_ext.rsplit2 ~on:':' line) with
| Some (_ ,buffer) -> | Some (_ ,buffer) ->
begin match (String_ext.split ~on:'=' buffer |> List.map String.trim) with begin match (String_ext.split ~on:'=' buffer |> list_map String.trim) with
| ["bit_kind_size"; x] -> | ["bit_kind_size"; x] ->
int_of_string x |> Bit_kind_size.of_int int_of_string x |> Bit_kind_size.of_int
| _ -> get_data tail | _ -> get_data tail
@ -58,7 +58,7 @@ let executables = lazy (
result result
in in
lines lines
|> List.map (fun x -> |> list_map (fun x ->
let e = String_ext.split ~on:' ' x let e = String_ext.split ~on:' ' x
|> List.filter (fun x -> x <> "") |> List.filter (fun x -> x <> "")
in in

3
ocaml/Qputils.ml
View File

@ -53,3 +53,6 @@ let input_lines ic =
let string_of_string s = s let string_of_string s = s
let list_map f l =
List.rev_map f l
|> List.rev

42
ocaml/qp_create_ezfio.ml
View File

@ -38,7 +38,7 @@ let dummy_centers ~threshold ~molecule ~nuclei =
| _ -> assert false | _ -> assert false
in in
aux [] (n-1,n-1) aux [] (n-1,n-1)
|> List.map (fun (i,x,j,y,r) -> |> list_map (fun (i,x,j,y,r) ->
let f = let f =
x /. (x +. y) x /. (x +. y)
in in
@ -270,7 +270,7 @@ let run ?o b au c d m p cart xyz_file =
(* Write Pseudo *) (* Write Pseudo *)
let pseudo = let pseudo =
List.map (fun x -> list_map (fun x ->
match pseudo_channel x.Atom.element with match pseudo_channel x.Atom.element with
| Some channel -> Pseudo.read_element channel x.Atom.element | Some channel -> Pseudo.read_element channel x.Atom.element
| None -> Pseudo.empty x.Atom.element | None -> Pseudo.empty x.Atom.element
@ -292,7 +292,7 @@ let run ?o b au c d m p cart xyz_file =
|> Elec_beta_number.of_int; |> Elec_beta_number.of_int;
Molecule.nuclei = Molecule.nuclei =
let charges = let charges =
List.map (fun x -> Positive_int.to_int x.Pseudo.n_elec list_map (fun x -> Positive_int.to_int x.Pseudo.n_elec
|> Float.of_int) pseudo |> Float.of_int) pseudo
|> Array.of_list |> Array.of_list
in in
@ -315,13 +315,13 @@ let run ?o b au c d m p cart xyz_file =
(* Write Nuclei *) (* Write Nuclei *)
let labels = let labels =
List.map (fun x->Element.to_string x.Atom.element) nuclei list_map (fun x->Element.to_string x.Atom.element) nuclei
and charges = and charges =
List.map (fun x-> Atom.(Charge.to_float x.charge)) nuclei list_map (fun x-> Atom.(Charge.to_float x.charge)) nuclei
and coords = and coords =
(List.map (fun x-> x.Atom.coord.Point3d.x) nuclei) @ (list_map (fun x-> x.Atom.coord.Point3d.x) nuclei) @
(List.map (fun x-> x.Atom.coord.Point3d.y) nuclei) @ (list_map (fun x-> x.Atom.coord.Point3d.y) nuclei) @
(List.map (fun x-> x.Atom.coord.Point3d.z) nuclei) in (list_map (fun x-> x.Atom.coord.Point3d.z) nuclei) in
let nucl_num = (List.length labels) in let nucl_num = (List.length labels) in
Ezfio.set_nuclei_nucl_num nucl_num ; Ezfio.set_nuclei_nucl_num nucl_num ;
Ezfio.set_nuclei_nucl_label (Ezfio.ezfio_array_of_list Ezfio.set_nuclei_nucl_label (Ezfio.ezfio_array_of_list
@ -365,7 +365,7 @@ let run ?o b au c d m p cart xyz_file =
let kmax = let kmax =
Array.init (lmax+1) (fun i-> Array.init (lmax+1) (fun i->
List.map (fun x -> list_map (fun x ->
List.filter (fun (y,_) -> List.filter (fun (y,_) ->
(Positive_int.to_int y.Pseudo.GaussianPrimitive_non_local.proj) = i) (Positive_int.to_int y.Pseudo.GaussianPrimitive_non_local.proj) = i)
x.Pseudo.non_local x.Pseudo.non_local
@ -478,7 +478,7 @@ let run ?o b au c d m p cart xyz_file =
in in
let result = do_work [] 1 nuclei let result = do_work [] 1 nuclei
|> List.rev |> List.rev
|> List.map (fun (x,i) -> |> list_map (fun (x,i) ->
try try
let e = let e =
match x.Atom.element with match x.Atom.element with
@ -512,30 +512,30 @@ let run ?o b au c d m p cart xyz_file =
let ao_num = List.length long_basis in let ao_num = List.length long_basis in
Ezfio.set_ao_basis_ao_num ao_num; Ezfio.set_ao_basis_ao_num ao_num;
Ezfio.set_ao_basis_ao_basis b; Ezfio.set_ao_basis_ao_basis b;
let ao_prim_num = List.map (fun (_,g,_) -> List.length g.Gto.lc) long_basis let ao_prim_num = list_map (fun (_,g,_) -> List.length g.Gto.lc) long_basis
and ao_nucl = List.map (fun (_,_,n) -> Nucl_number.to_int n) long_basis and ao_nucl = list_map (fun (_,_,n) -> Nucl_number.to_int n) long_basis
and ao_power= and ao_power=
let l = List.map (fun (x,_,_) -> x) long_basis in let l = list_map (fun (x,_,_) -> x) long_basis in
(List.map (fun t -> Positive_int.to_int Symmetry.Xyz.(t.x)) l)@ (list_map (fun t -> Positive_int.to_int Symmetry.Xyz.(t.x)) l)@
(List.map (fun t -> Positive_int.to_int Symmetry.Xyz.(t.y)) l)@ (list_map (fun t -> Positive_int.to_int Symmetry.Xyz.(t.y)) l)@
(List.map (fun t -> Positive_int.to_int Symmetry.Xyz.(t.z)) l) (list_map (fun t -> Positive_int.to_int Symmetry.Xyz.(t.z)) l)
in in
let ao_prim_num_max = List.fold_left (fun s x -> let ao_prim_num_max = List.fold_left (fun s x ->
if x > s then x if x > s then x
else s) 0 ao_prim_num else s) 0 ao_prim_num
in in
let gtos = let gtos =
List.map (fun (_,x,_) -> x) long_basis list_map (fun (_,x,_) -> x) long_basis
in in
let create_expo_coef ec = let create_expo_coef ec =
let coefs = let coefs =
begin match ec with begin match ec with
| `Coefs -> List.map (fun x-> | `Coefs -> list_map (fun x->
List.map (fun (_,coef) -> list_map (fun (_,coef) ->
AO_coef.to_float coef) x.Gto.lc) gtos AO_coef.to_float coef) x.Gto.lc) gtos
| `Expos -> List.map (fun x-> | `Expos -> list_map (fun x->
List.map (fun (prim,_) -> AO_expo.to_float list_map (fun (prim,_) -> AO_expo.to_float
prim.GaussianPrimitive.expo) x.Gto.lc) gtos prim.GaussianPrimitive.expo) x.Gto.lc) gtos
end end
in in

11
src/ao_one_e_ints/ao_ortho_canonical.irp.f
View File

@ -79,7 +79,7 @@ BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_inv, (ao_cart_to_sphe_num,ao_
call get_pseudo_inverse(ao_cart_to_sphe_coef,size(ao_cart_to_sphe_coef,1),& call get_pseudo_inverse(ao_cart_to_sphe_coef,size(ao_cart_to_sphe_coef,1),&
ao_num,ao_cart_to_sphe_num, & ao_num,ao_cart_to_sphe_num, &
ao_cart_to_sphe_inv, size(ao_cart_to_sphe_inv,1)) ao_cart_to_sphe_inv, size(ao_cart_to_sphe_inv,1), lin_dep_cutoff)
END_PROVIDER END_PROVIDER
@ -107,16 +107,13 @@ END_PROVIDER
ao_ortho_canonical_coef(i,i) = 1.d0 ao_ortho_canonical_coef(i,i) = 1.d0
enddo enddo
!call ortho_lowdin(ao_overlap,size(ao_overlap,1),ao_num,ao_ortho_canonical_coef,size(ao_ortho_canonical_coef,1),ao_num) call write_double(6, lin_dep_cutoff, "Linear dependencies cut-off")
!ao_ortho_canonical_num=ao_num
!return
if (ao_cartesian) then if (ao_cartesian) then
ao_ortho_canonical_num = ao_num ao_ortho_canonical_num = ao_num
call ortho_canonical(ao_overlap,size(ao_overlap,1), & call ortho_canonical(ao_overlap,size(ao_overlap,1), &
ao_num,ao_ortho_canonical_coef,size(ao_ortho_canonical_coef,1), & ao_num,ao_ortho_canonical_coef,size(ao_ortho_canonical_coef,1), &
ao_ortho_canonical_num) ao_ortho_canonical_num, lin_dep_cutoff)
else else
@ -131,7 +128,7 @@ END_PROVIDER
ao_ortho_canonical_num = ao_cart_to_sphe_num ao_ortho_canonical_num = ao_cart_to_sphe_num
call ortho_canonical(ao_cart_to_sphe_overlap, size(ao_cart_to_sphe_overlap,1), & call ortho_canonical(ao_cart_to_sphe_overlap, size(ao_cart_to_sphe_overlap,1), &
ao_cart_to_sphe_num, S, size(S,1), ao_ortho_canonical_num) ao_cart_to_sphe_num, S, size(S,1), ao_ortho_canonical_num, lin_dep_cutoff)
call dgemm('N','N', ao_num, ao_ortho_canonical_num, ao_cart_to_sphe_num, 1.d0, & call dgemm('N','N', ao_num, ao_ortho_canonical_num, ao_cart_to_sphe_num, 1.d0, &
ao_cart_to_sphe_coef, size(ao_cart_to_sphe_coef,1), & ao_cart_to_sphe_coef, size(ao_cart_to_sphe_coef,1), &

7
src/ao_one_e_ints/ao_overlap.irp.f
View File

@ -162,7 +162,8 @@ BEGIN_PROVIDER [ double precision, S_inv,(ao_num,ao_num) ]
BEGIN_DOC BEGIN_DOC
! Inverse of the overlap matrix ! Inverse of the overlap matrix
END_DOC END_DOC
call get_pseudo_inverse(ao_overlap,size(ao_overlap,1),ao_num,ao_num,S_inv,size(S_inv,1)) call get_pseudo_inverse(ao_overlap,size(ao_overlap,1),ao_num,ao_num,S_inv, &
size(S_inv,1),lin_dep_cutoff)
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ complex*16, S_inv_complex,(ao_num,ao_num) ] BEGIN_PROVIDER [ complex*16, S_inv_complex,(ao_num,ao_num) ]
@ -170,8 +171,8 @@ BEGIN_PROVIDER [ complex*16, S_inv_complex,(ao_num,ao_num) ]
BEGIN_DOC BEGIN_DOC
! Inverse of the overlap matrix ! Inverse of the overlap matrix
END_DOC END_DOC
call get_pseudo_inverse_complex(ao_overlap_complex, & call get_pseudo_inverse_complex(ao_overlap_complex, size(ao_overlap_complex,1),&
size(ao_overlap_complex,1),ao_num,ao_num,S_inv_complex,size(S_inv_complex,1)) ao_num,ao_num,S_inv_complex,size(S_inv_complex,1),lin_dep_cutoff)
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ double precision, S_half_inv, (AO_num,AO_num) ] BEGIN_PROVIDER [ double precision, S_half_inv, (AO_num,AO_num) ]

2
src/mo_guess/mo_ortho_lowdin.irp.f
View File

@ -13,7 +13,7 @@ BEGIN_PROVIDER [double precision, ao_ortho_lowdin_coef, (ao_num,ao_num)]
do j=1, ao_num do j=1, ao_num
tmp_matrix(j,j) = 1.d0 tmp_matrix(j,j) = 1.d0
enddo enddo
call ortho_lowdin(ao_overlap,ao_num,ao_num,tmp_matrix,ao_num,ao_num) call ortho_lowdin(ao_overlap,ao_num,ao_num,tmp_matrix,ao_num,ao_num,lin_dep_cutoff)
do i=1, ao_num do i=1, ao_num
do j=1, ao_num do j=1, ao_num
ao_ortho_lowdin_coef(j,i) = tmp_matrix(i,j) ao_ortho_lowdin_coef(j,i) = tmp_matrix(i,j)

5
src/mo_one_e_ints/EZFIO.cfg
View File

@ -48,3 +48,8 @@ doc: Read/Write |MO| one-electron integrals from/to disk [ Write | Read | None ]
interface: ezfio,provider,ocaml interface: ezfio,provider,ocaml
default: None default: None
[lin_dep_cutoff]
type: Threshold
doc: Remove linear dependencies when the eigenvalues of the overlap matrix are below this value
interface: ezfio,provider,ocaml
default: 1.e-6

2
src/mo_one_e_ints/orthonormalize.irp.f
View File

@ -3,7 +3,7 @@ subroutine orthonormalize_mos
integer :: m,p,s integer :: m,p,s
m = size(mo_coef,1) m = size(mo_coef,1)
p = size(mo_overlap,1) p = size(mo_overlap,1)
call ortho_lowdin(mo_overlap,p,mo_num,mo_coef,m,ao_num) call ortho_lowdin(mo_overlap,p,mo_num,mo_coef,m,ao_num,lin_dep_cutoff)
mo_label = 'Orthonormalized' mo_label = 'Orthonormalized'
SOFT_TOUCH mo_coef mo_label SOFT_TOUCH mo_coef mo_label
end end

2
src/scf_utils/huckel.irp.f
View File

@ -23,8 +23,6 @@ subroutine huckel_guess
Fock_matrix_ao_alpha(1:ao_num,1:ao_num) = A(1:ao_num,1:ao_num) Fock_matrix_ao_alpha(1:ao_num,1:ao_num) = A(1:ao_num,1:ao_num)
Fock_matrix_ao_beta (1:ao_num,1:ao_num) = A(1:ao_num,1:ao_num) Fock_matrix_ao_beta (1:ao_num,1:ao_num) = A(1:ao_num,1:ao_num)
! TOUCH mo_coef
TOUCH Fock_matrix_ao_alpha Fock_matrix_ao_beta TOUCH Fock_matrix_ao_alpha Fock_matrix_ao_beta
mo_coef = eigenvectors_fock_matrix_mo mo_coef = eigenvectors_fock_matrix_mo
SOFT_TOUCH mo_coef SOFT_TOUCH mo_coef

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

@ -23,7 +23,7 @@ subroutine extrapolate_data(N_data, data, pt2, output)
x(i,2) = pt2_rev(i) x(i,2) = pt2_rev(i)
enddo enddo
do ifit=2,N_data do ifit=2,N_data
call get_pseudo_inverse(x,size(x,1),ifit,2,x_inv,size(x_inv,1)) call get_pseudo_inverse(x,size(x,1),ifit,2,x_inv,size(x_inv,1),1.d-10)
ab = matmul(x_inv(1:2,1:ifit),y(1:ifit)) ab = matmul(x_inv(1:2,1:ifit),y(1:ifit))
output(ifit) = ab(1) output(ifit) = ab(1)
enddo enddo

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

@ -89,7 +89,7 @@ subroutine svd_complex(A,LDA,U,LDU,D,Vt,LDVt,m,n)
end end
subroutine ortho_canonical_complex(overlap,LDA,N,C,LDC,m) subroutine ortho_canonical_complex(overlap,LDA,N,C,LDC,m,cutoff)
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! Compute C_new=C_old.U.s^-1/2 canonical orthogonalization. ! Compute C_new=C_old.U.s^-1/2 canonical orthogonalization.
@ -112,6 +112,7 @@ subroutine ortho_canonical_complex(overlap,LDA,N,C,LDC,m)
integer, intent(in) :: lda, ldc, n integer, intent(in) :: lda, ldc, n
integer, intent(out) :: m integer, intent(out) :: m
complex*16, intent(in) :: overlap(lda,n) complex*16, intent(in) :: overlap(lda,n)
double precision, intent(in) :: cutoff
complex*16, intent(inout) :: C(ldc,n) complex*16, intent(inout) :: C(ldc,n)
complex*16, allocatable :: U(:,:) complex*16, allocatable :: U(:,:)
complex*16, allocatable :: Vt(:,:) complex*16, allocatable :: Vt(:,:)
@ -131,7 +132,7 @@ subroutine ortho_canonical_complex(overlap,LDA,N,C,LDC,m)
D(:) = dsqrt(D(:)) D(:) = dsqrt(D(:))
m=n m=n
do i=1,n do i=1,n
if ( D(i) >= 1.d-6 ) then if ( D(i) >= cutoff ) then
D(i) = 1.d0/D(i) D(i) = 1.d0/D(i)
else else
m = i-1 m = i-1
@ -229,7 +230,7 @@ subroutine ortho_qr_unblocked_complex(A,LDA,m,n)
! deallocate(WORK,jpvt,tau) ! deallocate(WORK,jpvt,tau)
end end
subroutine ortho_lowdin_complex(overlap,LDA,N,C,LDC,m) subroutine ortho_lowdin_complex(overlap,LDA,N,C,LDC,m,cutoff)
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! Compute C_new=C_old.S^-1/2 orthogonalization. ! Compute C_new=C_old.S^-1/2 orthogonalization.
@ -256,6 +257,7 @@ subroutine ortho_lowdin_complex(overlap,LDA,N,C,LDC,m)
complex*16, allocatable :: Vt(:,:) complex*16, allocatable :: Vt(:,:)
double precision, allocatable :: D(:) double precision, allocatable :: D(:)
complex*16, allocatable :: S(:,:) complex*16, allocatable :: S(:,:)
double precision, intent(in) :: cutoff
integer :: info, i, j, k integer :: info, i, j, k
if (n < 2) then if (n < 2) then
@ -267,12 +269,13 @@ subroutine ortho_lowdin_complex(overlap,LDA,N,C,LDC,m)
call svd_complex(overlap,lda,U,ldc,D,Vt,lda,n,n) call svd_complex(overlap,lda,U,ldc,D,Vt,lda,n,n)
!$OMP PARALLEL DEFAULT(NONE) & !$OMP PARALLEL DEFAULT(NONE) &
!$OMP SHARED(S,U,D,Vt,n,C,m) & !$OMP SHARED(S,U,D,Vt,n,C,m,cutoff) &
!$OMP PRIVATE(i,j,k) !$OMP PRIVATE(i,j,k)
!$OMP DO !$OMP DO
do i=1,n do i=1,n
if ( D(i) < 1.d-6 ) then if ( D(i) < cutoff) then
print *, 'Removed Linear dependencies :', 1.d0/D(i)
D(i) = 0.d0 D(i) = 0.d0
else else
D(i) = 1.d0/dsqrt(D(i)) D(i) = 1.d0/dsqrt(D(i))
@ -340,13 +343,14 @@ subroutine get_inverse_complex(A,LDA,m,C,LDC)
end end
subroutine get_pseudo_inverse_complex(A,LDA,m,n,C,LDC) subroutine get_pseudo_inverse_complex(A,LDA,m,n,C,LDC,cutoff)
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! Find C = A^-1 ! Find C = A^-1
END_DOC END_DOC
integer, intent(in) :: m,n, LDA, LDC integer, intent(in) :: m,n, LDA, LDC
complex*16, intent(in) :: A(LDA,n) complex*16, intent(in) :: A(LDA,n)
double precision, intent(in) :: cutoff
complex*16, intent(out) :: C(LDC,m) complex*16, intent(out) :: C(LDC,m)
double precision, allocatable :: D(:), rwork(:) double precision, allocatable :: D(:), rwork(:)
@ -375,7 +379,7 @@ subroutine get_pseudo_inverse_complex(A,LDA,m,n,C,LDC)
endif endif
do i=1,n do i=1,n
if (D(i)/D(1) > 1.d-10) then if (D(i)/D(1) > cutoff) then
D(i) = 1.d0/D(i) D(i) = 1.d0/D(i)
else else
D(i) = 0.d0 D(i) = 0.d0
@ -727,7 +731,7 @@ subroutine matrix_vector_product_complex(u0,u1,matrix,sze,lda)
call zhemv('U', sze, (1.d0,0.d0), matrix, lda, u0, incx, (1.d0,0.d0), u1, incy) call zhemv('U', sze, (1.d0,0.d0), matrix, lda, u0, incx, (1.d0,0.d0), u1, incy)
end end
subroutine ortho_canonical(overlap,LDA,N,C,LDC,m) subroutine ortho_canonical(overlap,LDA,N,C,LDC,m,cutoff)
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! Compute C_new=C_old.U.s^-1/2 canonical orthogonalization. ! Compute C_new=C_old.U.s^-1/2 canonical orthogonalization.
@ -750,6 +754,7 @@ subroutine ortho_canonical(overlap,LDA,N,C,LDC,m)
integer, intent(in) :: lda, ldc, n integer, intent(in) :: lda, ldc, n
integer, intent(out) :: m integer, intent(out) :: m
double precision, intent(in) :: overlap(lda,n) double precision, intent(in) :: overlap(lda,n)
double precision, intent(in) :: cutoff
double precision, intent(inout) :: C(ldc,n) double precision, intent(inout) :: C(ldc,n)
double precision, allocatable :: U(:,:) double precision, allocatable :: U(:,:)
double precision, allocatable :: Vt(:,:) double precision, allocatable :: Vt(:,:)
@ -769,7 +774,7 @@ subroutine ortho_canonical(overlap,LDA,N,C,LDC,m)
D(:) = dsqrt(D(:)) D(:) = dsqrt(D(:))
m=n m=n
do i=1,n do i=1,n
if ( D(i) >= 1.d-6 ) then if ( D(i) >= cutoff ) then
D(i) = 1.d0/D(i) D(i) = 1.d0/D(i)
else else
m = i-1 m = i-1
@ -874,7 +879,7 @@ subroutine ortho_qr_unblocked(A,LDA,m,n)
deallocate(WORK,TAU) deallocate(WORK,TAU)
end end
subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m) subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m,cutoff)
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! Compute C_new=C_old.S^-1/2 orthogonalization. ! Compute C_new=C_old.S^-1/2 orthogonalization.
@ -896,6 +901,7 @@ subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m)
integer, intent(in) :: LDA, ldc, n, m integer, intent(in) :: LDA, ldc, n, m
double precision, intent(in) :: overlap(lda,n) double precision, intent(in) :: overlap(lda,n)
double precision, intent(in) :: cutoff
double precision, intent(inout) :: C(ldc,n) double precision, intent(inout) :: C(ldc,n)
double precision, allocatable :: U(:,:) double precision, allocatable :: U(:,:)
double precision, allocatable :: Vt(:,:) double precision, allocatable :: Vt(:,:)
@ -912,12 +918,13 @@ subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m)
call svd(overlap,lda,U,ldc,D,Vt,lda,n,n) call svd(overlap,lda,U,ldc,D,Vt,lda,n,n)
!$OMP PARALLEL DEFAULT(NONE) & !$OMP PARALLEL DEFAULT(NONE) &
!$OMP SHARED(S,U,D,Vt,n,C,m) & !$OMP SHARED(S,U,D,Vt,n,C,m,cutoff) &
!$OMP PRIVATE(i,j,k) !$OMP PRIVATE(i,j,k)
!$OMP DO !$OMP DO
do i=1,n do i=1,n
if ( D(i) < 1.d-6 ) then if ( D(i) < cutoff ) then
print *, 'Removed Linear dependencies :', 1.d0/D(i)
D(i) = 0.d0 D(i) = 0.d0
else else
D(i) = 1.d0/dsqrt(D(i)) D(i) = 1.d0/dsqrt(D(i))
@ -986,13 +993,14 @@ subroutine get_inverse(A,LDA,m,C,LDC)
deallocate(ipiv,work) deallocate(ipiv,work)
end end
subroutine get_pseudo_inverse(A,LDA,m,n,C,LDC) subroutine get_pseudo_inverse(A,LDA,m,n,C,LDC,cutoff)
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! Find C = A^-1 ! Find C = A^-1
END_DOC END_DOC
integer, intent(in) :: m,n, LDA, LDC integer, intent(in) :: m,n, LDA, LDC
double precision, intent(in) :: A(LDA,n) double precision, intent(in) :: A(LDA,n)
double precision, intent(in) :: cutoff
double precision, intent(out) :: C(LDC,m) double precision, intent(out) :: C(LDC,m)
double precision, allocatable :: U(:,:), D(:), Vt(:,:), work(:), A_tmp(:,:) double precision, allocatable :: U(:,:), D(:), Vt(:,:), work(:), A_tmp(:,:)
@ -1020,7 +1028,7 @@ subroutine get_pseudo_inverse(A,LDA,m,n,C,LDC)
endif endif
do i=1,n do i=1,n
if (D(i)/D(1) > 1.d-10) then if (D(i)/D(1) > cutoff) then
D(i) = 1.d0/D(i) D(i) = 1.d0/D(i)
else else
D(i) = 0.d0 D(i) = 0.d0
@ -1053,7 +1061,7 @@ subroutine find_rotation(A,LDA,B,m,C,n)
double precision, allocatable :: A_inv(:,:) double precision, allocatable :: A_inv(:,:)
allocate(A_inv(LDA,n)) allocate(A_inv(LDA,n))
call get_pseudo_inverse(A,LDA,m,n,A_inv,LDA) call get_pseudo_inverse(A,LDA,m,n,A_inv,LDA,1.d-10)
integer :: i,j,k integer :: i,j,k
call dgemm('N','N',n,n,m,1.d0,A_inv,n,B,LDA,0.d0,C,n) call dgemm('N','N',n,n,m,1.d0,A_inv,n,B,LDA,0.d0,C,n)