mirror of
https://gitlab.com/scemama/QCaml.git
synced 2024-09-01 01:03:41 +02:00
Compare commits
4 Commits
701ac05e23
...
1e778594ec
Author | SHA1 | Date | |
---|---|---|---|
1e778594ec | |||
16a0fbacaa | |||
58f51dbecf | |||
91a2ec1288 |
16
.merlin
16
.merlin
@ -1,12 +1,8 @@
|
||||
PKG str unix bigarray lacaml zarith mpi alcotest
|
||||
EXCLUDE_QUERY_DIR
|
||||
B _build/default/.run_hartree_fock.eobjs/byte
|
||||
B _build/default/Basis/.Basis.objs/byte
|
||||
B _build/default/Nuclei/.Nuclei.objs/byte
|
||||
S .
|
||||
S Test
|
||||
S Nuclei
|
||||
S Parallel
|
||||
S CI
|
||||
S Utils
|
||||
S Basis
|
||||
S SCF
|
||||
S MOBasis
|
||||
S CI
|
||||
B _build/**
|
||||
S Nuclei
|
||||
FLG -w @a-4-29-40-41-42-44-45-48-58-59-60-40 -strict-sequence -strict-formats -short-paths -keep-locs
|
||||
|
@ -756,14 +756,14 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
|
||||
in
|
||||
let empty = Array.make (maxm+1) 0. in
|
||||
let center_qc_tmp = Array.init nq (fun cd ->
|
||||
Bohr.make { Point.
|
||||
Coordinate.make { Coordinate.
|
||||
x = (center_qc Co.X).(cd) ;
|
||||
y = (center_qc Co.Y).(cd) ;
|
||||
z = (center_qc Co.Z).(cd) ;
|
||||
})
|
||||
in
|
||||
Array.iteri (fun ab shell_ab ->
|
||||
let center_pa = Bohr.make { Point.
|
||||
let center_pa = Coordinate.make { Coordinate.
|
||||
x = (center_pa Co.X).(ab) ;
|
||||
y = (center_pa Co.Y).(ab) ;
|
||||
z = (center_pa Co.Z).(ab) ;
|
||||
@ -789,7 +789,7 @@ let contracted_class_shell_pairs ~zero_m ?schwartz_p ?schwartz_q shell_p shell_q
|
||||
in
|
||||
zero_m Zp.{
|
||||
maxm ; expo_p_inv ; expo_q_inv ; norm_pq_sq ;
|
||||
center_pq = Bohr.make Point.{x ; y ; z} ;
|
||||
center_pq = Coordinate.make Coordinate.{x ; y ; z} ;
|
||||
center_pa ; center_qc = center_qc_tmp.(cd) ;
|
||||
normalization ;
|
||||
}
|
||||
|
@ -28,7 +28,9 @@ let of_xyz_file filename =
|
||||
let of_zmt_string buffer =
|
||||
Zmatrix.of_string buffer
|
||||
|> Zmatrix.to_xyz
|
||||
|> Array.map (fun (e,x,y,z) -> (e, Coordinate.angstrom_to_bohr (Angstrom.make {Point.x ; y ; z} )))
|
||||
|> Array.map (fun (e,x,y,z) ->
|
||||
(e, Coordinate.(angstrom_to_bohr @@ make_angstrom { x ; y ; z} ))
|
||||
)
|
||||
|
||||
|
||||
let of_zmt_file filename =
|
||||
@ -56,12 +58,13 @@ let to_string atoms =
|
||||
-----------------------------------------------------------------------
|
||||
" ^
|
||||
(Array.mapi (fun i (e, coord) ->
|
||||
let open Coordinate in
|
||||
let coord =
|
||||
Coordinate.bohr_to_angstrom coord
|
||||
bohr_to_angstrom coord
|
||||
in
|
||||
Printf.sprintf " %5d %5d %5s %12.6f %12.6f %12.6f"
|
||||
(i+1) (Element.to_int e) (Element.to_string e)
|
||||
coord.Angstrom.x coord.Angstrom.y coord.Angstrom.z
|
||||
coord.x coord.y coord.z
|
||||
) atoms
|
||||
|> Array.to_list
|
||||
|> String.concat "\n" ) ^
|
||||
@ -100,11 +103,12 @@ let charge nuclei =
|
||||
let to_xyz_string t =
|
||||
[ string_of_int (Array.length t) ; "" ] @
|
||||
( Array.mapi (fun i (e, coord) ->
|
||||
let open Coordinate in
|
||||
let coord =
|
||||
Coordinate.bohr_to_angstrom coord
|
||||
bohr_to_angstrom coord
|
||||
in
|
||||
Printf.sprintf " %5s %12.6f %12.6f %12.6f"
|
||||
(Element.to_string e) coord.Angstrom.x coord.Angstrom.y coord.Angstrom.z
|
||||
(Element.to_string e) coord.x coord.y coord.z
|
||||
) t
|
||||
|> Array.to_list )
|
||||
|> String.concat "\n"
|
||||
@ -117,8 +121,9 @@ let to_t2_string t =
|
||||
(2 * Array.length t);
|
||||
"# Znuc x y z" ]
|
||||
@ (Array.mapi (fun i (e, coord) ->
|
||||
let open Coordinate in
|
||||
Printf.sprintf " %5f %12.6f %12.6f %12.6f"
|
||||
(Element.to_int e |> float_of_int) coord.Bohr.x coord.Bohr.y coord.Bohr.z
|
||||
(Element.to_int e |> float_of_int) coord.x coord.y coord.z
|
||||
) t
|
||||
|> Array.to_list )
|
||||
|> String.concat "\n"
|
||||
|
@ -2,7 +2,7 @@
|
||||
of tuples ({!Element.t}, {!Coordinate.t}).
|
||||
*)
|
||||
|
||||
type t = (Element.t * Bohr.t) array
|
||||
type t = (Element.t * Coordinate.t) array
|
||||
|
||||
|
||||
val of_xyz_file : string -> t
|
||||
|
@ -4,7 +4,7 @@
|
||||
type nucleus =
|
||||
{
|
||||
element: Element.t ;
|
||||
coord : Angstrom.t;
|
||||
coord : Coordinate.angstrom Coordinate.point;
|
||||
}
|
||||
|
||||
type xyz_file =
|
||||
|
@ -3,10 +3,9 @@
|
||||
%{
|
||||
exception InputError of string
|
||||
let make_angstrom x y z =
|
||||
Angstrom.make {
|
||||
Point.
|
||||
Coordinate.(make_angstrom {
|
||||
x ; y ; z
|
||||
}
|
||||
})
|
||||
|
||||
let output_of f x y z =
|
||||
let a = make_angstrom x y z in
|
||||
|
@ -1 +0,0 @@
|
||||
include Bohr
|
@ -1,4 +0,0 @@
|
||||
(** 3D point in cartesian coordinates, where units are Angstroms. *)
|
||||
|
||||
include module type of Bohr
|
||||
|
@ -1,11 +0,0 @@
|
||||
type t = {
|
||||
x : float ;
|
||||
y : float ;
|
||||
z : float ;
|
||||
}
|
||||
|
||||
|
||||
external make : Point.t -> t = "%identity"
|
||||
|
||||
|
||||
|
@ -1,12 +0,0 @@
|
||||
(** 3D point in cartesian coordinates, where units are atomic units (Bohr). *)
|
||||
|
||||
type t = private {
|
||||
x : float ;
|
||||
y : float ;
|
||||
z : float ;
|
||||
}
|
||||
|
||||
|
||||
val make : Point.t -> t
|
||||
|
||||
|
@ -3,13 +3,16 @@ open Lacaml.D
|
||||
let full_ldl m_A =
|
||||
|
||||
let n = Mat.dim1 m_A in
|
||||
|
||||
assert (Mat.dim2 m_A = n);
|
||||
|
||||
let v_D = Vec.make0 n in
|
||||
|
||||
let m_Lt = Mat.identity n in
|
||||
|
||||
let v_D_inv = Vec.make0 n in
|
||||
|
||||
|
||||
let compute_d j =
|
||||
let l_jk =
|
||||
Mat.col m_Lt j
|
||||
@ -20,6 +23,7 @@ let full_ldl m_A =
|
||||
m_A.{j,j} -. dot ~n:(j-1) l_jk l_jk__d_k
|
||||
in
|
||||
|
||||
|
||||
let compute_l i =
|
||||
let l_ik__d_k =
|
||||
Mat.col m_Lt i
|
||||
@ -33,6 +37,7 @@ let full_ldl m_A =
|
||||
v_D_inv.{j} *. ( m_A.{j,i} -. dot ~n:(j-1) l_ik__d_k l_jk )
|
||||
in
|
||||
|
||||
|
||||
for i=1 to n do
|
||||
for j=1 to (i-1) do
|
||||
m_Lt.{j,i} <- compute_l i j;
|
||||
@ -85,6 +90,7 @@ let pivoted_ldl threshold m_A =
|
||||
v_D_inv.{j} *. ( m_A.{pi.(j),pi.(i)} -. dot ~n:(j-1) l_ik__d_k l_jk )
|
||||
in
|
||||
|
||||
|
||||
let maxloc v i =
|
||||
let rec aux pos value i =
|
||||
if i > n then
|
||||
@ -103,7 +109,9 @@ let pivoted_ldl threshold m_A =
|
||||
let p_i, p_j = pi.(i), pi.(j) in
|
||||
pi.(i) <- p_j;
|
||||
pi.(j) <- p_i;
|
||||
in
|
||||
in
|
||||
|
||||
|
||||
|
||||
|
||||
let () =
|
||||
@ -124,7 +132,11 @@ let pivoted_ldl threshold m_A =
|
||||
m_Lt, v_D, pi
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
let make_ldl ?(threshold=Constants.epsilon) m_A =
|
||||
pivoted_ldl m_A
|
||||
|
||||
|
||||
|
||||
|
||||
|
@ -1,7 +1,19 @@
|
||||
type bohr = Bohr.t
|
||||
type angstrom = Angstrom.t
|
||||
type bohr
|
||||
type angstrom
|
||||
|
||||
type t = bohr
|
||||
type xyz = {
|
||||
x : float ;
|
||||
y : float ;
|
||||
z : float ;
|
||||
}
|
||||
|
||||
type 'a point = xyz
|
||||
|
||||
type t = bohr point
|
||||
|
||||
external make : 'a point -> t = "%identity"
|
||||
external make_angstrom : 'a point -> angstrom point = "%identity"
|
||||
|
||||
|
||||
type axis = X | Y | Z
|
||||
|
||||
@ -9,19 +21,17 @@ type axis = X | Y | Z
|
||||
let a_to_b a = Constants.a0_inv *. a
|
||||
let b_to_a b = Constants.a0 *. b
|
||||
|
||||
let bohr_to_angstrom { Bohr.x ; y ; z } =
|
||||
Angstrom.make
|
||||
let bohr_to_angstrom { x ; y ; z } =
|
||||
make
|
||||
{
|
||||
Point.
|
||||
x = b_to_a x ;
|
||||
y = b_to_a y ;
|
||||
z = b_to_a z ;
|
||||
}
|
||||
|
||||
let angstrom_to_bohr { Angstrom.x ; y ; z } =
|
||||
Bohr.make
|
||||
let angstrom_to_bohr { x ; y ; z } =
|
||||
make
|
||||
{
|
||||
Point.
|
||||
x = a_to_b x ;
|
||||
y = a_to_b y ;
|
||||
z = a_to_b z ;
|
||||
@ -29,30 +39,27 @@ let angstrom_to_bohr { Angstrom.x ; y ; z } =
|
||||
|
||||
|
||||
let zero =
|
||||
Bohr.make { Point.x = 0. ; y = 0. ; z = 0. }
|
||||
make { x = 0. ; y = 0. ; z = 0. }
|
||||
|
||||
|
||||
(** Linear algebra *)
|
||||
let ( |. ) s { Bohr.x ; y ; z } =
|
||||
Bohr.make {
|
||||
Point.
|
||||
let ( |. ) s { x ; y ; z } =
|
||||
make {
|
||||
x = s *. x ;
|
||||
y = s *. y ;
|
||||
z = s *. z ;
|
||||
}
|
||||
|
||||
let ( |+ ) { Bohr.x = x1 ; y = y1 ; z = z1 } { Bohr.x = x2 ; y = y2 ; z = z2 } =
|
||||
Bohr.make {
|
||||
Point.
|
||||
let ( |+ ) { x = x1 ; y = y1 ; z = z1 } { x = x2 ; y = y2 ; z = z2 } =
|
||||
make {
|
||||
x = x1 +. x2 ;
|
||||
y = y1 +. y2 ;
|
||||
z = z1 +. z2 ;
|
||||
}
|
||||
|
||||
|
||||
let ( |- ) { Bohr.x = x1 ; y = y1 ; z = z1 } { Bohr.x = x2 ; y = y2 ; z = z2 } =
|
||||
Bohr.make {
|
||||
Point.
|
||||
let ( |- ) { x = x1 ; y = y1 ; z = z1 } { x = x2 ; y = y2 ; z = z2 } =
|
||||
make {
|
||||
x = x1 -. x2 ;
|
||||
y = y1 -. y2 ;
|
||||
z = z1 -. z2 ;
|
||||
@ -62,7 +69,7 @@ let ( |- ) { Bohr.x = x1 ; y = y1 ; z = z1 } { Bohr.x = x2 ; y = y2 ; z = z2 } =
|
||||
let neg a = -1. |. a
|
||||
|
||||
|
||||
let dot { Bohr.x = x1 ; y = y1 ; z = z1 } { Bohr.x = x2 ; y = y2 ; z = z2 } =
|
||||
let dot { x = x1 ; y = y1 ; z = z1 } { x = x2 ; y = y2 ; z = z2 } =
|
||||
x1 *. x2 +. y1 *. y2 +. z1 *. z2
|
||||
|
||||
|
||||
@ -70,7 +77,7 @@ let norm u =
|
||||
sqrt ( dot u u )
|
||||
|
||||
|
||||
let get axis { Bohr.x ; y ; z } =
|
||||
let get axis { x ; y ; z } =
|
||||
match axis with
|
||||
| X -> x
|
||||
| Y -> y
|
||||
@ -79,14 +86,12 @@ let get axis { Bohr.x ; y ; z } =
|
||||
|
||||
open Format
|
||||
let pp ppf c =
|
||||
let open Bohr in
|
||||
fprintf ppf "@[@[%8.4f@] @[%8.4f@] @[%8.4f@]@]" c.x c.y c.z
|
||||
|
||||
let pp_bohr ppf c =
|
||||
let open Bohr in
|
||||
fprintf ppf "@[(@[%10f@], @[%10f@], @[%10f@] Bohr)@]" c.x c.y c.z
|
||||
|
||||
let pp_angstrom ppf c =
|
||||
let c = bohr_to_angstrom c in
|
||||
let open Angstrom in
|
||||
fprintf ppf "@[(@[%10f@], @[%10f@], @[%10f@] Angs)@]" c.x c.y c.z
|
||||
|
||||
|
@ -1,92 +1,108 @@
|
||||
(** Type for coordinates in 3D space.
|
||||
|
||||
All operations on points are done in atomic units. Therefore,
|
||||
all coordinates are given in {!Bohr.t} and manipulated with this
|
||||
all coordinates are given in bohr and manipulated with this
|
||||
module.
|
||||
*)
|
||||
|
||||
type t = Bohr.t
|
||||
type bohr
|
||||
type angstrom
|
||||
|
||||
type xyz = {
|
||||
x : float ;
|
||||
y : float ;
|
||||
z : float ;
|
||||
}
|
||||
|
||||
type 'a point = xyz
|
||||
|
||||
type t = bohr point
|
||||
|
||||
val make : 'a point -> t
|
||||
(** Creates a point in atomic units *)
|
||||
|
||||
val make_angstrom : 'a point -> angstrom point
|
||||
(** Creates a point in Angstrom *)
|
||||
|
||||
type axis = X | Y | Z
|
||||
|
||||
val bohr_to_angstrom : Bohr.t -> Angstrom.t
|
||||
(** Converts a point in Bohr to Angstrom. *)
|
||||
val bohr_to_angstrom : bohr point -> angstrom point
|
||||
(** Converts a point in bohr to angstrom. *)
|
||||
|
||||
|
||||
val angstrom_to_bohr : Angstrom.t -> Bohr.t
|
||||
(** Converts a point in Angstrom to Bohr. *)
|
||||
val angstrom_to_bohr : angstrom point -> bohr point
|
||||
(** Converts a point in Angstrom to bohr. *)
|
||||
|
||||
|
||||
val zero : Bohr.t
|
||||
(** [zero = { Bohr.x = 0. ; y=0. ; z=0. }] *)
|
||||
val zero : bohr point
|
||||
(** [zero = { x = 0. ; y=0. ; z=0. }] *)
|
||||
|
||||
|
||||
val get : axis -> Bohr.t -> float
|
||||
val get : axis -> bohr point -> float
|
||||
(** Extracts the projection of the coordinate on an axis.
|
||||
|
||||
Example:
|
||||
|
||||
[Coordinate.(get Y) { Bohr.x=1. ; y=2. ; z=3. } -> 2.]
|
||||
[Coordinate.(get Y) { x=1. ; y=2. ; z=3. } -> 2.]
|
||||
|
||||
*)
|
||||
|
||||
|
||||
(** {1 Vector operations} *)
|
||||
|
||||
|
||||
val ( |. ) : float -> Bohr.t -> Bohr.t
|
||||
val ( |. ) : float -> t -> t
|
||||
(** Scale by a float.
|
||||
|
||||
Example:
|
||||
|
||||
[ 2. |. { Bohr.x=1. ; y=2. ; z=3. } -> { Bohr.x=2. ; y=4. ; z=6. } ]
|
||||
[ 2. |. { x=1. ; y=2. ; z=3. } -> { x=2. ; y=4. ; z=6. } ]
|
||||
|
||||
*)
|
||||
|
||||
|
||||
val ( |+ ) : Bohr.t -> Bohr.t -> Bohr.t
|
||||
val ( |+ ) : t -> t -> t
|
||||
(** Add two vectors.
|
||||
|
||||
Example:
|
||||
|
||||
{[{ Bohr.x=1. ; y=2. ; z=3. } |+ { Bohr.x=2. ; y=3. ; z=1. } ->
|
||||
{ Bohr.x=3. ; y=5. ; z=4. }]}
|
||||
{[{ x=1. ; y=2. ; z=3. } |+ { x=2. ; y=3. ; z=1. } ->
|
||||
{ x=3. ; y=5. ; z=4. }]}
|
||||
|
||||
*)
|
||||
|
||||
|
||||
val ( |- ) : Bohr.t -> Bohr.t -> Bohr.t
|
||||
val ( |- ) : t -> t -> t
|
||||
(** Subtract two vectors.
|
||||
|
||||
Example:
|
||||
|
||||
{[{ Bohr.x=1. ; y=2. ; z=3. } |- { Bohr.x=2. ; y=3. ; z=1. } ->
|
||||
{ Bohr.x=-1. ; y=-1. ; z=2. }]}
|
||||
{[{ x=1. ; y=2. ; z=3. } |- { x=2. ; y=3. ; z=1. } ->
|
||||
{ x=-1. ; y=-1. ; z=2. }]}
|
||||
|
||||
*)
|
||||
|
||||
|
||||
val neg : Bohr.t -> Bohr.t
|
||||
val neg : t -> t
|
||||
(** Example:
|
||||
|
||||
{[Coordinate.neg { Bohr.x=1. ; y=2. ; z=-3. } ->
|
||||
{ Bohr.x=-1. ; y=-2. ; z=3. }]}
|
||||
{[Coordinate.neg { x=1. ; y=2. ; z=-3. } ->
|
||||
{ x=-1. ; y=-2. ; z=3. }]}
|
||||
|
||||
*)
|
||||
|
||||
val dot : Bohr.t -> Bohr.t -> float
|
||||
val dot : t -> t -> float
|
||||
(** Dot product. *)
|
||||
|
||||
|
||||
val norm : Bohr.t -> float
|
||||
val norm : t -> float
|
||||
(** L{^2} norm of the vector. *)
|
||||
|
||||
|
||||
(** {2 Printers} *)
|
||||
|
||||
val pp: Format.formatter -> Bohr.t -> unit
|
||||
val pp: Format.formatter -> t -> unit
|
||||
|
||||
val pp_bohr: Format.formatter -> Bohr.t -> unit
|
||||
val pp_bohr: Format.formatter -> t -> unit
|
||||
|
||||
val pp_angstrom : Format.formatter -> Bohr.t -> unit
|
||||
val pp_angstrom : Format.formatter -> t -> unit
|
||||
|
||||
|
Loading…
Reference in New Issue
Block a user