mirror of https://gitlab.com/scemama/QCaml.git
Compare commits
5 Commits
d6ef25d55a
...
7eba47e0a4
| Author | SHA1 | Date |
|---|---|---|
 Anthony Scemama
|
7eba47e0a4 | |
|
Anthony Scemama
|
50698ce0ff | |
|
Anthony Scemama
|
ac44492ffe | |
|
Anthony Scemama
|
7e6ff1e8e7 | |
|
Anthony Scemama
|
feebc5a0e8 |
|
|
@ -1,3 +0,0 @@
|
|||
[submodule "docs/org-html-themes"]
|
||||
path = docs/org-html-themes
|
||||
url = https://github.com/fniessen/org-html-themes.git
|
||||
|
|
@ -18,7 +18,9 @@ cd $DIR
|
|||
|
||||
for i in *.org
|
||||
do
|
||||
if [[ $(basename $i) != "README.org" ]] ; then
|
||||
echo "--- $i ----"
|
||||
emacs --batch ./$i --load=../docs/config_tangle.el -f org-babel-tangle
|
||||
fi
|
||||
done
|
||||
|
||||
|
|
|
|||
|
|
@ -1 +0,0 @@
|
|||
Subproject commit 885ddd8451a3b1952628a52c4eeebf0a3d86d33a
|
||||
|
|
@ -1,7 +1,5 @@
|
|||
open Linear_algebra
|
||||
open Common
|
||||
|
||||
(** Types *)
|
||||
open Linear_algebra
|
||||
|
||||
type localization_kind =
|
||||
| Edmiston
|
||||
|
|
@ -29,6 +27,8 @@ type t =
|
|||
selected_mos : int list ;
|
||||
}
|
||||
|
||||
open Common
|
||||
|
||||
|
||||
(** Edmiston-Rudenberg *)
|
||||
|
||||
|
|
@ -197,8 +197,6 @@ let kappa_boys in_basis =
|
|||
r2 (phi_x_phi%:(i,i)) (phi_y_phi%:(i,i)) (phi_z_phi%:(i,i)))
|
||||
|> Vector.sum
|
||||
)
|
||||
(* Boys:1 ends here *)
|
||||
|
||||
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -1,7 +1,4 @@
|
|||
(** Orbital localization *)
|
||||
|
||||
|
||||
(** Types *)
|
||||
(** Molecular orbital localization *)
|
||||
|
||||
open Linear_algebra
|
||||
|
||||
|
|
@ -20,15 +17,23 @@ type t
|
|||
(** Access *)
|
||||
|
||||
val kind : t -> localization_kind
|
||||
(** Returns the kind of localized MOs *)
|
||||
|
||||
val simulation : t -> Simulation.t
|
||||
(** Returns the simulation environment in which the MOs are localized *)
|
||||
|
||||
val selected_mos : t -> int list
|
||||
(** List of indices of the orbitals involved in the localization *)
|
||||
|
||||
|
||||
val kappa :
|
||||
kind:localization_kind ->
|
||||
Basis.t ->
|
||||
( ao,loc) Matrix.t ->
|
||||
(loc,loc) Matrix.t * float
|
||||
(** Returns the $\kappa$ antisymmetric matrix used for the rotation matrix and the value of the localization function *)
|
||||
(** Returns the $\kappa$ antisymmetric matrix used for the rotation matrix and
|
||||
* the value of the localization function
|
||||
*)
|
||||
|
||||
val make :
|
||||
kind:localization_kind ->
|
||||
|
|
@ -44,4 +49,4 @@ val to_basis : t -> Basis.t
|
|||
|
||||
(** Printers *)
|
||||
|
||||
val pp : Format.formatter -> t -> unit
|
||||
val pp : Format.formatter -> t -> unit
|
||||
|
|
|
|||
|
|
@ -1,520 +0,0 @@
|
|||
#+begin_src elisp tangle: no :results none :exports none
|
||||
(setq pwd (file-name-directory buffer-file-name))
|
||||
(setq name (file-name-nondirectory (substring buffer-file-name 0 -4)))
|
||||
(setq lib (concat pwd "lib/"))
|
||||
(setq testdir (concat pwd "test/"))
|
||||
(setq mli (concat lib name ".mli"))
|
||||
(setq ml (concat lib name ".ml"))
|
||||
(setq test-ml (concat testdir name ".ml"))
|
||||
(org-babel-tangle)
|
||||
#+end_src
|
||||
|
||||
* Orbital localization
|
||||
:PROPERTIES:
|
||||
:header-args: :noweb yes :comments both
|
||||
:END:
|
||||
|
||||
Molecular orbital localization function.
|
||||
|
||||
Boys:
|
||||
|
||||
Edmiston-Rudenberg:
|
||||
|
||||
|
||||
** Type
|
||||
|
||||
#+NAME: types
|
||||
#+begin_src ocaml :tangle (eval mli)
|
||||
open Linear_algebra
|
||||
|
||||
type localization_kind =
|
||||
| Edmiston
|
||||
| Boys
|
||||
|
||||
type mo = Mo_dim.t
|
||||
type ao = Ao.Ao_dim.t
|
||||
type loc
|
||||
#+end_src
|
||||
|
||||
#+begin_src ocaml :tangle (eval mli)
|
||||
type localization_data
|
||||
type t
|
||||
#+end_src
|
||||
|
||||
#+begin_src ocaml :tangle (eval ml) :exports none
|
||||
<<types>>
|
||||
|
||||
type localization_data =
|
||||
{
|
||||
coefficients : (ao, loc) Matrix.t ;
|
||||
kappa : (loc, loc) Matrix.t ;
|
||||
scaling : float ;
|
||||
loc_value : float ;
|
||||
convergence : float ;
|
||||
iteration : int ;
|
||||
}
|
||||
|
||||
type t =
|
||||
{
|
||||
kind : localization_kind ;
|
||||
mo_basis : Basis.t ;
|
||||
data : localization_data option lazy_t array ;
|
||||
selected_mos : int list ;
|
||||
}
|
||||
|
||||
open Common
|
||||
#+end_src
|
||||
|
||||
** Edmiston-Rudenberg
|
||||
|
||||
#+begin_src ocaml :tangle (eval ml) :exports none
|
||||
let kappa_edmiston in_basis m_C =
|
||||
let ao_basis =
|
||||
Basis.simulation in_basis
|
||||
|> Simulation.ao_basis
|
||||
in
|
||||
|
||||
let ee_ints = Ao.Basis.ee_ints ao_basis in
|
||||
let n_ao = Ao.Basis.size ao_basis in
|
||||
|
||||
let n_mo = Matrix.dim2 m_C in
|
||||
|
||||
(* Temp arrays for integral transformation *)
|
||||
let m_pqr =
|
||||
Bigarray.(Array3.create Float64 fortran_layout n_ao n_ao n_ao)
|
||||
in
|
||||
let m_qr_i = Matrix.create (n_ao*n_ao) n_mo in
|
||||
let m_ri_j = Matrix.create (n_ao*n_mo) n_mo in
|
||||
let m_ij_k = Matrix.create (n_mo*n_mo) n_mo in
|
||||
|
||||
|
||||
let m_a12 = Bigarray.(Array2.create Float64 fortran_layout n_mo n_mo) in
|
||||
let m_b12 = Bigarray.(Array2.create Float64 fortran_layout n_mo n_mo) in
|
||||
Bigarray.Array2.fill m_b12 0.;
|
||||
Bigarray.Array2.fill m_a12 0.;
|
||||
let v_d =
|
||||
Vector.init n_mo (fun _ -> 0.)
|
||||
|> Vector.to_bigarray_inplace
|
||||
in
|
||||
|
||||
Array.iter (fun s ->
|
||||
|
||||
Array.iter (fun r ->
|
||||
Array.iter (fun q ->
|
||||
Array.iter (fun p ->
|
||||
m_pqr.{p,q,r} <- Four_idx_storage.get_phys ee_ints p q r s
|
||||
) (Util.array_range 1 n_ao)
|
||||
) (Util.array_range 1 n_ao)
|
||||
) (Util.array_range 1 n_ao);
|
||||
|
||||
let m_p_qr =
|
||||
Bigarray.reshape (Bigarray.genarray_of_array3 m_pqr) [| n_ao ; n_ao*n_ao |]
|
||||
|> Bigarray.array2_of_genarray
|
||||
|> Matrix.of_bigarray_inplace
|
||||
in
|
||||
|
||||
(* (qr,i) = <i r|q s> = \sum_p <p r | q s> C_{pi} *)
|
||||
Matrix.gemm_tn_inplace ~c:m_qr_i m_p_qr m_C;
|
||||
|
||||
let m_q_ri =
|
||||
let x = Matrix.to_bigarray_inplace m_qr_i |> Bigarray.genarray_of_array2 in
|
||||
Bigarray.reshape_2 x n_ao (n_ao*n_mo) |> Matrix.of_bigarray_inplace
|
||||
in
|
||||
|
||||
(* (ri,j) = <i r | j s> = \sum_q <i r | q s> C_{qj} *)
|
||||
Matrix.gemm_tn_inplace ~c:m_ri_j m_q_ri m_C;
|
||||
|
||||
let m_r_ij =
|
||||
let x = Matrix.to_bigarray_inplace m_ri_j |> Bigarray.genarray_of_array2 in
|
||||
Bigarray.reshape_2 x n_ao (n_mo*n_mo) |> Matrix.of_bigarray_inplace
|
||||
in
|
||||
|
||||
(* (ij,k) = <i k | j s> = \sum_r <i r | j s> C_{rk} *)
|
||||
Matrix.gemm_tn_inplace ~c:m_ij_k m_r_ij m_C;
|
||||
|
||||
let m_ijk =
|
||||
let x = Matrix.to_bigarray_inplace m_ij_k |> Bigarray.genarray_of_array2 in
|
||||
Bigarray.reshape x [| n_mo ; n_mo ; n_mo |]
|
||||
|> Bigarray.array3_of_genarray
|
||||
in
|
||||
|
||||
let m_Cx = Matrix.to_bigarray_inplace m_C in
|
||||
Array.iter (fun j ->
|
||||
Array.iter (fun i ->
|
||||
m_b12.{i,j} <- m_b12.{i,j} +. m_Cx.{s,j} *. (m_ijk.{i,i,i} -. m_ijk.{j,i,j});
|
||||
m_a12.{i,j} <- m_a12.{i,j} +. m_ijk.{i,i,j} *. m_Cx.{s,j} -.
|
||||
0.25 *. ( (m_ijk.{i,i,i} -. m_ijk.{j,i,j}) *. m_Cx.{s,i} +.
|
||||
(m_ijk.{j,j,j} -. m_ijk.{i,j,i}) *. m_Cx.{s,j})
|
||||
) (Util.array_range 1 n_mo);
|
||||
v_d.{j} <- v_d.{j} +. m_ijk.{j,j,j} *. m_Cx.{s,j}
|
||||
) (Util.array_range 1 n_mo)
|
||||
|
||||
) (Util.array_range 1 n_ao);
|
||||
|
||||
let f i j =
|
||||
if i=j then
|
||||
0.
|
||||
else
|
||||
begin
|
||||
let x = 1./. sqrt (m_b12.{i,j} *. m_b12.{i,j} +. m_a12.{i,j} *. m_a12.{i,j} ) in
|
||||
if asin (m_b12.{i,j} /. x) > 0. then
|
||||
0.25 *. acos( -. m_a12.{i,j} *. x)
|
||||
else
|
||||
-. 0.25 *. acos( -. m_a12.{i,j} *. x )
|
||||
end
|
||||
in
|
||||
(
|
||||
Matrix.init_cols n_mo n_mo ( fun i j -> if i<=j then f i j else -. (f j i) ),
|
||||
Vector.sum (Vector.of_bigarray_inplace v_d)
|
||||
)
|
||||
|
||||
|
||||
|
||||
#+end_src
|
||||
|
||||
** Boys
|
||||
|
||||
#+begin_src ocaml :tangle (eval ml) :exports none
|
||||
let kappa_boys in_basis =
|
||||
let ao_basis =
|
||||
Basis.simulation in_basis
|
||||
|> Simulation.ao_basis
|
||||
in
|
||||
let multipole = Ao.Basis.multipole ao_basis in
|
||||
fun m_C ->
|
||||
let n_mo = Matrix.dim2 m_C in
|
||||
|
||||
let phi_x_phi = Matrix.xt_o_x ~x:m_C ~o:(multipole "x") in
|
||||
let phi_y_phi = Matrix.xt_o_x ~x:m_C ~o:(multipole "y") in
|
||||
let phi_z_phi = Matrix.xt_o_x ~x:m_C ~o:(multipole "z") in
|
||||
|
||||
let m_b12 =
|
||||
let g x i j =
|
||||
let x_ii = x%:(i,i) in
|
||||
let x_jj = x%:(j,j) in
|
||||
let x_ij = x%:(i,j) in
|
||||
(x_ii -. x_jj) *. x_ij
|
||||
in
|
||||
Matrix.init_cols n_mo n_mo (fun i j ->
|
||||
let x =
|
||||
(g phi_x_phi i j) +. (g phi_y_phi i j) +. (g phi_z_phi i j)
|
||||
in
|
||||
if (abs_float x > 1.e-15) then x else 0.
|
||||
)
|
||||
in
|
||||
|
||||
let m_a12 =
|
||||
let g x i j =
|
||||
let x_ii = x%:(i,i) in
|
||||
let x_jj = x%:(j,j) in
|
||||
let x_ij = x%:(i,j) in
|
||||
let y = x_ii -. x_jj in
|
||||
(x_ij *. x_ij) -. 0.25 *. y *. y
|
||||
in
|
||||
Matrix.init_cols n_mo n_mo (fun i j ->
|
||||
let x =
|
||||
(g phi_x_phi i j) +. (g phi_y_phi i j) +. (g phi_z_phi i j)
|
||||
in
|
||||
if (abs_float x > 1.e-15) then x else 0.
|
||||
)
|
||||
in
|
||||
|
||||
let f i j =
|
||||
let pi = Constants.pi in
|
||||
if i=j then
|
||||
0.
|
||||
else
|
||||
let x = atan2 (m_b12%:(i,j)) (m_a12%:(i,j)) in
|
||||
if x >= 0. then
|
||||
0.25 *. (pi -. x)
|
||||
else
|
||||
-. 0.25 *. ( pi +. x)
|
||||
in
|
||||
(
|
||||
Matrix.init_cols n_mo n_mo ( fun i j -> if i>j then f i j else -. (f j i) ),
|
||||
let r2 x y z = x*.x +. y*.y +. z*.z in
|
||||
Vector.init n_mo ( fun i ->
|
||||
r2 (phi_x_phi%:(i,i)) (phi_y_phi%:(i,i)) (phi_z_phi%:(i,i)))
|
||||
|> Vector.sum
|
||||
)
|
||||
|
||||
#+end_src
|
||||
|
||||
|
||||
** Access
|
||||
|
||||
#+begin_src ocaml :tangle (eval mli)
|
||||
val kind : t -> localization_kind
|
||||
val simulation : t -> Simulation.t
|
||||
val selected_mos : t -> int list
|
||||
|
||||
val kappa :
|
||||
kind:localization_kind ->
|
||||
Basis.t ->
|
||||
( ao,loc) Matrix.t ->
|
||||
(loc,loc) Matrix.t * float
|
||||
|
||||
val make :
|
||||
kind:localization_kind ->
|
||||
?max_iter:int ->
|
||||
?convergence:float ->
|
||||
Basis.t ->
|
||||
int list ->
|
||||
t
|
||||
|
||||
val to_basis : t -> Basis.t
|
||||
|
||||
#+end_src
|
||||
|
||||
| ~kappa~ | Returns the $\kappa$ antisymmetric matrix used for the rotation matrix and the value of the localization function |
|
||||
| ~make~ | Performs the orbital localization |
|
||||
|
||||
#+begin_src ocaml :tangle (eval ml) :exports none
|
||||
let kind t = t.kind
|
||||
let simulation t = Basis.simulation t.mo_basis
|
||||
let selected_mos t = t.selected_mos
|
||||
|
||||
let kappa ~kind =
|
||||
match kind with
|
||||
| Edmiston -> kappa_edmiston
|
||||
| Boys -> kappa_boys
|
||||
|
||||
|
||||
let n_iterations t =
|
||||
Array.fold_left (fun accu x ->
|
||||
match Lazy.force x with
|
||||
| Some _ -> accu + 1
|
||||
| None -> accu
|
||||
) 0 t.data
|
||||
|
||||
let last_iteration t =
|
||||
Util.of_some @@ Lazy.force t.data.(n_iterations t - 1)
|
||||
|
||||
(*
|
||||
let ao_basis t = Simulation.ao_basis (simulation t)
|
||||
,*)
|
||||
|
||||
|
||||
let make ~kind ?(max_iter=500) ?(convergence=1.e-6) mo_basis selected_mos =
|
||||
|
||||
let kappa_loc = kappa ~kind mo_basis in
|
||||
|
||||
let mo_array = Matrix.to_col_vecs (Basis.mo_coef mo_basis) in
|
||||
let mos_vec_list = List.map (fun i -> mo_array.(i-1)) selected_mos in
|
||||
let x: (ao,loc) Matrix.t =
|
||||
Matrix.of_col_vecs_list mos_vec_list
|
||||
in
|
||||
|
||||
let next_coef kappa x =
|
||||
let r = Matrix.exponential_antisymmetric kappa in
|
||||
let m_C = Matrix.gemm_nt x r in
|
||||
m_C
|
||||
in
|
||||
|
||||
let data_initial =
|
||||
let iteration = 0
|
||||
and scaling = 0.1
|
||||
in
|
||||
let kappa, loc_value = kappa_loc x in
|
||||
let convergence = abs_float (Matrix.amax kappa) in
|
||||
let kappa = Matrix.scale scaling kappa in
|
||||
let coefficients = next_coef kappa x in
|
||||
{ coefficients ; kappa ; scaling ; convergence ; loc_value ; iteration }
|
||||
in
|
||||
|
||||
let iteration data =
|
||||
let iteration = data.iteration+1 in
|
||||
let new_kappa, new_loc_value = kappa_loc data.coefficients in
|
||||
let new_convergence = abs_float (Matrix.amax new_kappa) in
|
||||
let accept =
|
||||
new_loc_value >= data.loc_value*.0.98
|
||||
in
|
||||
if accept then
|
||||
let coefficients = next_coef new_kappa data.coefficients in
|
||||
let scaling = min 1. (data.scaling *. 1.1) in
|
||||
let kappa = Matrix.scale scaling new_kappa in
|
||||
let convergence = new_convergence in
|
||||
let loc_value = new_loc_value in
|
||||
{ coefficients ; kappa ; scaling ; convergence ; loc_value ; iteration }
|
||||
else
|
||||
let scaling =
|
||||
data.scaling *. 0.5
|
||||
in
|
||||
{ data with scaling ; iteration }
|
||||
in
|
||||
|
||||
let array_data =
|
||||
|
||||
let storage =
|
||||
Array.make max_iter None
|
||||
in
|
||||
|
||||
let rec loop = function
|
||||
| 0 -> Some (iteration data_initial)
|
||||
| n -> begin
|
||||
match storage.(n) with
|
||||
| Some result -> Some result
|
||||
| None -> begin
|
||||
let data = loop (n-1) in
|
||||
match data with
|
||||
| None -> None
|
||||
| Some data -> begin
|
||||
(* Check convergence *)
|
||||
let converged =
|
||||
data.convergence < convergence
|
||||
in
|
||||
if converged then
|
||||
None
|
||||
else
|
||||
begin
|
||||
storage.(n-1) <- Some data ;
|
||||
storage.(n) <- Some (iteration data);
|
||||
storage.(n)
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
in
|
||||
Array.init max_iter (fun i -> lazy (loop i))
|
||||
in
|
||||
{ kind ; mo_basis ; data = array_data ; selected_mos }
|
||||
|
||||
|
||||
|
||||
let to_basis t =
|
||||
let mo_basis = t.mo_basis in
|
||||
let simulation = Basis.simulation mo_basis in
|
||||
let mo_occupation = Basis.mo_occupation mo_basis in
|
||||
|
||||
let data = last_iteration t in
|
||||
|
||||
let mo_coef_array = Matrix.to_col_vecs (Basis.mo_coef mo_basis) in
|
||||
let new_mos =
|
||||
Matrix.to_col_vecs data.coefficients
|
||||
in
|
||||
selected_mos t
|
||||
|> List.iteri (fun i j -> mo_coef_array.(j-1) <- new_mos.(i)) ;
|
||||
let mo_coef = Matrix.of_col_vecs mo_coef_array in
|
||||
Basis.make ~simulation ~mo_type:(Localized "Boys") ~mo_occupation ~mo_coef ()
|
||||
|
||||
#+end_src
|
||||
|
||||
** Printers
|
||||
|
||||
#+begin_src ocaml :tangle (eval mli)
|
||||
val pp : Format.formatter -> t -> unit
|
||||
#+end_src
|
||||
|
||||
#+begin_src ocaml :tangle (eval ml) :exports none
|
||||
let linewidth = 60
|
||||
|
||||
let pp_iterations ppf t =
|
||||
let line = (String.make linewidth '-') in
|
||||
Format.fprintf ppf "@[%4s%s@]@." "" line;
|
||||
Format.fprintf ppf "@[%4s@[%5s@]@,@[%16s@]@,@[%16s@]@,@[%11s@]@]@."
|
||||
"" "#" "Localization " "Convergence" "Scaling";
|
||||
Format.fprintf ppf "@[%4s%s@]@." "" line;
|
||||
Array.iter (fun data ->
|
||||
let data = Lazy.force data in
|
||||
match data with
|
||||
| None -> ()
|
||||
| Some data ->
|
||||
let loc = data.loc_value in
|
||||
let conv = data.convergence in
|
||||
let scaling = data.scaling in
|
||||
let iteration = data.iteration in
|
||||
begin
|
||||
Format.fprintf ppf "@[%4s@[%5d@]@,@[%16.8f@]@,@[%16.4e@]@,@[%11.4f@]@]@." ""
|
||||
iteration loc conv scaling ;
|
||||
end
|
||||
) t.data;
|
||||
Format.fprintf ppf "@[%4s%s@]@." "" line
|
||||
|
||||
|
||||
let pp ppf t =
|
||||
Format.fprintf ppf "@.@[%s@]@." (String.make 70 '=');
|
||||
Format.fprintf ppf "@[%34s %-34s@]@." (match t.kind with
|
||||
| Boys -> "Boys"
|
||||
| Edmiston -> "Edmiston-Ruedenberg"
|
||||
) "MO Localization";
|
||||
Format.fprintf ppf "@[%s@]@.@." (String.make 70 '=');
|
||||
Format.fprintf ppf "@[%a@]@." pp_iterations t;
|
||||
|
||||
#+end_src
|
||||
|
||||
** Tests
|
||||
|
||||
#+begin_src ocaml :tangle (eval test-ml) :exports none
|
||||
|
||||
let test_localization =
|
||||
|
||||
let nuclei =
|
||||
Particles.Nuclei.of_xyz_string
|
||||
" 10
|
||||
Hydrogen chain, d=1.8 Angstrom
|
||||
H -4.286335 0.000000 0.000000
|
||||
H -3.333816 0.000000 0.000000
|
||||
H -2.381297 0.000000 0.000000
|
||||
H -1.428778 0.000000 0.000000
|
||||
H -0.476259 0.000000 0.000000
|
||||
H 0.476259 0.000000 0.000000
|
||||
H 1.428778 0.000000 0.000000
|
||||
H 2.381297 0.000000 0.000000
|
||||
H 3.333816 0.000000 0.000000
|
||||
H 4.286335 0.000000 0.000000
|
||||
" in
|
||||
|
||||
let basis_file = "/home/scemama/qp2/data/basis/sto-6g" in
|
||||
|
||||
let ao_basis =
|
||||
Ao.Basis.of_nuclei_and_basis_filename ~nuclei basis_file
|
||||
in
|
||||
|
||||
|
||||
let charge = 0 in
|
||||
|
||||
let multiplicity = 1 in
|
||||
|
||||
let simulation =
|
||||
Simulation.make ~charge ~multiplicity ~nuclei ao_basis
|
||||
in
|
||||
|
||||
let hf =
|
||||
Mo.Hartree_fock.make ~guess:`Hcore simulation
|
||||
in
|
||||
|
||||
let mo_basis =
|
||||
Mo.Basis.of_hartree_fock hf
|
||||
in
|
||||
|
||||
let localized_mo_basis =
|
||||
Mo.Localization.make
|
||||
~kind:Mo.Localization.Boys
|
||||
mo_basis
|
||||
[4;5;6;7;8]
|
||||
|> Mo.Localization.to_basis
|
||||
in
|
||||
|
||||
Format.printf "%a" (Mo.Basis.pp ~start:1 ~finish:10) localized_mo_basis
|
||||
|
||||
|
||||
(*
|
||||
open Common
|
||||
open Alcotest
|
||||
let wd = Qcaml.root ^ Filename.dir_sep ^ "test" in
|
||||
|
||||
let test_xyz molecule length repulsion charge core =
|
||||
let xyz = Nuclei.of_xyz_file (wd^Filename.dir_sep^molecule^".xyz") in
|
||||
check int "length" length (Array.length xyz);
|
||||
check (float 1.e-4) "repulsion" repulsion (Nuclei.repulsion xyz);
|
||||
check int "charge" charge (Charge.to_int @@ Nuclei.charge xyz);
|
||||
check int "small_core" core (Nuclei.small_core xyz);
|
||||
()
|
||||
|
||||
let tests = [
|
||||
"caffeine", `Quick, (fun () -> test_xyz "caffeine" 24 917.0684 102 28);
|
||||
"water", `Quick, (fun () -> test_xyz "water" 3 9.19497 10 2);
|
||||
]
|
||||
,*)
|
||||
|
||||
#+end_src
|
||||
|
|
@ -1,5 +1,7 @@
|
|||
(* Tests *)
|
||||
|
||||
|
||||
(* [[file:~/QCaml/mo/localization.org::*Tests][Tests:1]] *)
|
||||
let test_localization =
|
||||
|
||||
let nuclei =
|
||||
|
|
@ -70,3 +72,4 @@ let tests = [
|
|||
"water", `Quick, (fun () -> test_xyz "water" 3 9.19497 10 2);
|
||||
]
|
||||
*)
|
||||
(* Tests:1 ends here *)
|
||||
|
|
|
|||
|
|
@ -1,160 +0,0 @@
|
|||
#+begin_src elisp tangle: no :results none :exports none
|
||||
(setq pwd (file-name-directory buffer-file-name))
|
||||
(setq name (file-name-nondirectory (substring buffer-file-name 0 -4)))
|
||||
(setq lib (concat pwd "lib/"))
|
||||
(setq testdir (concat pwd "test/"))
|
||||
(setq mli (concat lib name ".mli"))
|
||||
(setq ml (concat lib name ".ml"))
|
||||
(setq test-ml (concat testdir name ".ml"))
|
||||
(org-babel-tangle)
|
||||
#+end_src
|
||||
|
||||
* Electrons
|
||||
:PROPERTIES:
|
||||
:header-args: :noweb yes :comments both
|
||||
:END:
|
||||
|
||||
Data structure which contains the number of \alpha and \beta electrons.
|
||||
|
||||
** Type
|
||||
|
||||
#+NAME: types
|
||||
#+begin_src ocaml :tangle (eval mli)
|
||||
type t
|
||||
#+end_src
|
||||
|
||||
#+begin_src ocaml :tangle (eval ml) :exports none
|
||||
type t = {
|
||||
n_alfa : int ;
|
||||
n_beta : int ;
|
||||
}
|
||||
#+end_src
|
||||
|
||||
#+begin_src ocaml :tangle (eval test-ml) :exports none
|
||||
open Common
|
||||
open Particles
|
||||
open Alcotest
|
||||
|
||||
let test_all () =
|
||||
let nuclei =
|
||||
"3
|
||||
Water
|
||||
O 0. 0. 0.
|
||||
H -0.756950272703377558 0. -0.585882234512562827
|
||||
H 0.756950272703377558 0. -0.585882234512562827
|
||||
"
|
||||
|> Nuclei.of_xyz_string
|
||||
in
|
||||
let e = Electrons.of_atoms nuclei in
|
||||
#+end_src
|
||||
|
||||
** Creation
|
||||
|
||||
#+begin_src ocaml :tangle (eval mli)
|
||||
open Common
|
||||
|
||||
val make : int -> int -> t
|
||||
|
||||
val of_atoms : ?multiplicity:int -> ?charge:int -> Nuclei.t -> t
|
||||
(* @param multiplicity default is 1
|
||||
@param charge default is 0
|
||||
@raise Invalid_argument if the spin multiplicity is not compatible with
|
||||
the molecule and the total charge.
|
||||
*)
|
||||
#+end_src
|
||||
|
||||
| ~make~ | ~make n_alfa n_beta~ |
|
||||
| ~of_atoms~ | Creates the data relative to electrons for a molecular system described by ~Nuclei.t~ for a given total charge and spin multiplicity. |
|
||||
|
||||
#+begin_src ocaml :tangle (eval ml) :exports none
|
||||
open Common
|
||||
|
||||
let make n_alfa n_beta =
|
||||
{ n_alfa ; n_beta }
|
||||
|
||||
|
||||
let of_atoms ?multiplicity:(multiplicity=1) ?charge:(charge=0) nuclei =
|
||||
let positive_charges =
|
||||
Array.fold_left (fun accu (e, _) -> accu + Charge.to_int (Element.to_charge e) )
|
||||
0 nuclei
|
||||
in
|
||||
let negative_charges = charge - positive_charges in
|
||||
let n_elec = - negative_charges in
|
||||
let n_beta = ((n_elec - multiplicity)+1)/2 in
|
||||
let n_alfa = n_elec - n_beta in
|
||||
let result = { n_alfa ; n_beta } in
|
||||
if multiplicity <> (n_alfa - n_beta)+1 then
|
||||
invalid_arg (__FILE__^": make");
|
||||
result
|
||||
|
||||
#+end_src
|
||||
|
||||
#+begin_src ocaml :tangle (eval test-ml) :exports none
|
||||
check int "of_atoms alfa" 5 (Electrons.n_alfa e);
|
||||
check int "of_atoms beta" 5 (Electrons.n_beta e);
|
||||
#+end_src
|
||||
|
||||
** Access
|
||||
|
||||
#+begin_src ocaml :tangle (eval mli)
|
||||
val charge : t -> Charge.t
|
||||
val n_elec : t -> int
|
||||
val n_alfa : t -> int
|
||||
val n_beta : t -> int
|
||||
val multiplicity : t -> int
|
||||
#+end_src
|
||||
|
||||
| ~charge~ | Sum of the charges of the electrons |
|
||||
| ~n_elec~ | Number of electrons |
|
||||
| ~n_alfa~ | Number of alpha electrons |
|
||||
| ~n_beta~ | Number of beta electrons |
|
||||
| ~multiplicity~ | Spin multiplicity: $2S+1$ |
|
||||
|
||||
#+begin_src ocaml :tangle (eval ml) :exports none
|
||||
let charge e =
|
||||
- (e.n_alfa + e.n_beta)
|
||||
|> Charge.of_int
|
||||
|
||||
let n_alfa t = t.n_alfa
|
||||
|
||||
let n_beta t = t.n_beta
|
||||
|
||||
let n_elec t = t.n_alfa + t.n_beta
|
||||
|
||||
let multiplicity t = t.n_alfa - t.n_beta + 1
|
||||
#+end_src
|
||||
|
||||
#+begin_src ocaml :tangle (eval test-ml) :exports none
|
||||
check int "charge " (-10) (Charge.to_int @@ Electrons.charge e);
|
||||
check int "n_elec" 10 (Electrons.n_elec e);
|
||||
check int "multiplicity" 1 (Electrons.multiplicity e);
|
||||
check int "of_atoms alfa m3" 6 (Electrons.(of_atoms ~multiplicity:3 nuclei |> n_alfa));
|
||||
check int "of_atoms beta m3" 4 (Electrons.(of_atoms ~multiplicity:3 nuclei |> n_beta));
|
||||
check int "of_atoms n_elec m3" 10 (Electrons.(of_atoms ~multiplicity:3 nuclei |> n_elec));
|
||||
check int "of_atoms alfa m2 c1" 5 (Electrons.(of_atoms ~multiplicity:2 ~charge:1 nuclei |> n_alfa));
|
||||
check int "of_atoms beta m2 c1" 4 (Electrons.(of_atoms ~multiplicity:2 ~charge:1 nuclei |> n_beta));
|
||||
check int "of_atoms beta m2 c1" 9 (Electrons.(of_atoms ~multiplicity:2 ~charge:1 nuclei |> n_elec));
|
||||
check int "of_atoms mult m2 c1" 2 (Electrons.(of_atoms ~multiplicity:2 ~charge:1 nuclei |> multiplicity));
|
||||
check bool "make" true (Electrons.make 6 4 = Electrons.(of_atoms ~multiplicity:3 nuclei));
|
||||
#+end_src
|
||||
|
||||
** Printers
|
||||
|
||||
#+begin_src ocaml :tangle (eval mli)
|
||||
val pp : Format.formatter -> t -> unit
|
||||
#+end_src
|
||||
|
||||
#+begin_src ocaml :tangle (eval ml) :exports none
|
||||
let pp ppf t =
|
||||
Format.fprintf ppf "@[n_alfa=%d, n_beta=%d@]" t.n_alfa t.n_beta
|
||||
#+end_src
|
||||
|
||||
** Tests
|
||||
|
||||
#+begin_src ocaml :tangle (eval test-ml) :exports none
|
||||
()
|
||||
|
||||
let tests = [
|
||||
"all", `Quick, test_all
|
||||
]
|
||||
#+end_src
|
||||
|
|
@ -1,328 +0,0 @@
|
|||
#+begin_src elisp tangle: no :results none :exports none
|
||||
(setq pwd (file-name-directory buffer-file-name))
|
||||
(setq name (file-name-nondirectory (substring buffer-file-name 0 -4)))
|
||||
(setq lib (concat pwd "lib/"))
|
||||
(setq testdir (concat pwd "test/"))
|
||||
(setq mli (concat lib name ".mli"))
|
||||
(setq ml (concat lib name ".ml"))
|
||||
(setq test-ml (concat testdir name ".ml"))
|
||||
(org-babel-tangle)
|
||||
#+end_src
|
||||
|
||||
* Element
|
||||
:PROPERTIES:
|
||||
:header-args: :noweb yes :comments both
|
||||
:END:
|
||||
|
||||
Chemical elements.
|
||||
|
||||
** Type
|
||||
|
||||
#+NAME: types
|
||||
#+begin_src ocaml :tangle (eval mli)
|
||||
type t =
|
||||
|X
|
||||
|H |He
|
||||
|Li|Be |B |C |N |O |F |Ne
|
||||
|Na|Mg |Al|Si|P |S |Cl|Ar
|
||||
|K |Ca|Sc|Ti|V |Cr|Mn|Fe|Co|Ni|Cu|Zn|Ga|Ge|As|Se|Br|Kr
|
||||
|Rb|Sr|Y |Zr|Nb|Mo|Tc|Ru|Rh|Pd|Ag|Cd|In|Sn|Sb|Te|I |Xe
|
||||
|Pt
|
||||
|
||||
exception ElementError of string
|
||||
|
||||
open Common
|
||||
#+end_src
|
||||
|
||||
#+begin_src ocaml :tangle (eval ml) :exports none
|
||||
<<types>>
|
||||
#+end_src
|
||||
|
||||
** Conversion
|
||||
|
||||
#+begin_src ocaml :tangle (eval mli)
|
||||
val of_string : string -> t
|
||||
val to_string : t -> string
|
||||
val to_long_string : t -> string
|
||||
|
||||
val to_int : t -> int
|
||||
val of_int : int -> t
|
||||
|
||||
val to_charge : t -> Charge.t
|
||||
val of_charge : Charge.t -> t
|
||||
#+end_src
|
||||
|
||||
| ~of_string~ | Creates an ~Element.t~ from a chemical symbol or from the full name of the element (case insensitive) |
|
||||
| ~to_string~ | Gets the chemical symbol of the ~Element.t~ in a string |
|
||||
| ~to_long_string~ | Gets the full name of the ~Element.t~ in a string |
|
||||
| ~to_int~ | Convert to the atomic charge, with ~int~ type |
|
||||
| ~of_int~ | Create from the atomic charge, with ~int~ type |
|
||||
| ~to_charge~ | Convert to the atomic charge, with ~Charge.t~ type |
|
||||
| ~of_charge~ | Create from the atomic charge, with ~Charge.t~ type |
|
||||
|
||||
#+begin_example
|
||||
Element.of_string "Fe" ;;
|
||||
- : Element.t = Particles.Element.Fe
|
||||
|
||||
Element.of_string "hydrogen" ;;
|
||||
- : Element.t = Particles.Element.H
|
||||
|
||||
Element.of_string "Kryptonite" ;;
|
||||
Exception: Particles.Element.ElementError "Element Kryptonite unknown".
|
||||
|
||||
Element.(to_long_string Fe) ;;
|
||||
- : string = "Iron"
|
||||
|
||||
Element.(to_string Fe);;
|
||||
- : string = "Fe"
|
||||
#+end_example
|
||||
|
||||
#+begin_src ocaml :tangle (eval ml) :exports none
|
||||
let of_string x =
|
||||
match (String.capitalize_ascii (String.lowercase_ascii x)) with
|
||||
| "X" | "Dummy" -> X | "H" | "Hydrogen" -> H
|
||||
| "He" | "Helium" -> He | "Li" | "Lithium" -> Li
|
||||
| "Be" | "Beryllium" -> Be | "B" | "Boron" -> B
|
||||
| "C" | "Carbon" -> C | "N" | "Nitrogen" -> N
|
||||
| "O" | "Oxygen" -> O | "F" | "Fluorine" -> F
|
||||
| "Ne" | "Neon" -> Ne | "Na" | "Sodium" -> Na
|
||||
| "Mg" | "Magnesium" -> Mg | "Al" | "Aluminum" -> Al
|
||||
| "Si" | "Silicon" -> Si | "P" | "Phosphorus" -> P
|
||||
| "S" | "Sulfur" -> S | "Cl" | "Chlorine" -> Cl
|
||||
| "Ar" | "Argon" -> Ar | "K" | "Potassium" -> K
|
||||
| "Ca" | "Calcium" -> Ca | "Sc" | "Scandium" -> Sc
|
||||
| "Ti" | "Titanium" -> Ti | "V" | "Vanadium" -> V
|
||||
| "Cr" | "Chromium" -> Cr | "Mn" | "Manganese" -> Mn
|
||||
| "Fe" | "Iron" -> Fe | "Co" | "Cobalt" -> Co
|
||||
| "Ni" | "Nickel" -> Ni | "Cu" | "Copper" -> Cu
|
||||
| "Zn" | "Zinc" -> Zn | "Ga" | "Gallium" -> Ga
|
||||
| "Ge" | "Germanium" -> Ge | "As" | "Arsenic" -> As
|
||||
| "Se" | "Selenium" -> Se | "Br" | "Bromine" -> Br
|
||||
| "Kr" | "Krypton" -> Kr | "Rb" | "Rubidium" -> Rb
|
||||
| "Sr" | "Strontium" -> Sr | "Y" | "Yttrium" -> Y
|
||||
| "Zr" | "Zirconium" -> Zr | "Nb" | "Niobium" -> Nb
|
||||
| "Mo" | "Molybdenum" -> Mo | "Tc" | "Technetium" -> Tc
|
||||
| "Ru" | "Ruthenium" -> Ru | "Rh" | "Rhodium" -> Rh
|
||||
| "Pd" | "Palladium" -> Pd | "Ag" | "Silver" -> Ag
|
||||
| "Cd" | "Cadmium" -> Cd | "In" | "Indium" -> In
|
||||
| "Sn" | "Tin" -> Sn | "Sb" | "Antimony" -> Sb
|
||||
| "Te" | "Tellurium" -> Te | "I" | "Iodine" -> I
|
||||
| "Xe" | "Xenon" -> Xe | "Pt" | "Platinum" -> Pt
|
||||
| x -> raise (ElementError ("Element "^x^" unknown"))
|
||||
|
||||
|
||||
let to_string = function
|
||||
| X -> "X" | H -> "H" | He -> "He" | Li -> "Li"
|
||||
| Be -> "Be" | B -> "B" | C -> "C" | N -> "N"
|
||||
| O -> "O" | F -> "F" | Ne -> "Ne" | Na -> "Na"
|
||||
| Mg -> "Mg" | Al -> "Al" | Si -> "Si" | P -> "P"
|
||||
| S -> "S" | Cl -> "Cl" | Ar -> "Ar" | K -> "K"
|
||||
| Ca -> "Ca" | Sc -> "Sc" | Ti -> "Ti" | V -> "V"
|
||||
| Cr -> "Cr" | Mn -> "Mn" | Fe -> "Fe" | Co -> "Co"
|
||||
| Ni -> "Ni" | Cu -> "Cu" | Zn -> "Zn" | Ga -> "Ga"
|
||||
| Ge -> "Ge" | As -> "As" | Se -> "Se" | Br -> "Br"
|
||||
| Kr -> "Kr" | Rb -> "Rb" | Sr -> "Sr" | Y -> "Y"
|
||||
| Zr -> "Zr" | Nb -> "Nb" | Mo -> "Mo" | Tc -> "Tc"
|
||||
| Ru -> "Ru" | Rh -> "Rh" | Pd -> "Pd" | Ag -> "Ag"
|
||||
| Cd -> "Cd" | In -> "In" | Sn -> "Sn" | Sb -> "Sb"
|
||||
| Te -> "Te" | I -> "I" | Xe -> "Xe" | Pt -> "Pt"
|
||||
|
||||
let to_long_string = function
|
||||
| X -> "Dummy" | H -> "Hydrogen" | He -> "Helium"
|
||||
| Li -> "Lithium" | Be -> "Beryllium" | B -> "Boron"
|
||||
| C -> "Carbon" | N -> "Nitrogen" | O -> "Oxygen"
|
||||
| F -> "Fluorine" | Ne -> "Neon" | Na -> "Sodium"
|
||||
| Mg -> "Magnesium" | Al -> "Aluminum" | Si -> "Silicon"
|
||||
| P -> "Phosphorus" | S -> "Sulfur" | Cl -> "Chlorine"
|
||||
| Ar -> "Argon" | K -> "Potassium" | Ca -> "Calcium"
|
||||
| Sc -> "Scandium" | Ti -> "Titanium" | V -> "Vanadium"
|
||||
| Cr -> "Chromium" | Mn -> "Manganese" | Fe -> "Iron"
|
||||
| Co -> "Cobalt" | Ni -> "Nickel" | Cu -> "Copper"
|
||||
| Zn -> "Zinc" | Ga -> "Gallium" | Ge -> "Germanium"
|
||||
| As -> "Arsenic" | Se -> "Selenium" | Br -> "Bromine"
|
||||
| Kr -> "Krypton" | Rb -> "Rubidium" | Sr -> "Strontium"
|
||||
| Y -> "Yttrium" | Zr -> "Zirconium" | Nb -> "Niobium"
|
||||
| Mo -> "Molybdenum" | Tc -> "Technetium" | Ru -> "Ruthenium"
|
||||
| Rh -> "Rhodium" | Pd -> "Palladium" | Ag -> "Silver"
|
||||
| Cd -> "Cadmium" | In -> "Indium" | Sn -> "Tin"
|
||||
| Sb -> "Antimony" | Te -> "Tellurium" | I -> "Iodine"
|
||||
| Xe -> "Xenon" | Pt -> "Platinum"
|
||||
|
||||
|
||||
let to_int = function
|
||||
| X -> 0 | H -> 1 | He -> 2 | Li -> 3
|
||||
| Be -> 4 | B -> 5 | C -> 6 | N -> 7
|
||||
| O -> 8 | F -> 9 | Ne -> 10 | Na -> 11
|
||||
| Mg -> 12 | Al -> 13 | Si -> 14 | P -> 15
|
||||
| S -> 16 | Cl -> 17 | Ar -> 18 | K -> 19
|
||||
| Ca -> 20 | Sc -> 21 | Ti -> 22 | V -> 23
|
||||
| Cr -> 24 | Mn -> 25 | Fe -> 26 | Co -> 27
|
||||
| Ni -> 28 | Cu -> 29 | Zn -> 30 | Ga -> 31
|
||||
| Ge -> 32 | As -> 33 | Se -> 34 | Br -> 35
|
||||
| Kr -> 36 | Rb -> 37 | Sr -> 38 | Y -> 39
|
||||
| Zr -> 40 | Nb -> 41 | Mo -> 42 | Tc -> 43
|
||||
| Ru -> 44 | Rh -> 45 | Pd -> 46 | Ag -> 47
|
||||
| Cd -> 48 | In -> 49 | Sn -> 50 | Sb -> 51
|
||||
| Te -> 52 | I -> 53 | Xe -> 54 | Pt -> 78
|
||||
|
||||
|
||||
let to_charge c =
|
||||
to_int c |> Charge.of_int
|
||||
|
||||
|
||||
let of_int = function
|
||||
| 0 -> X | 1 -> H | 2 -> He | 3 -> Li
|
||||
| 4 -> Be | 5 -> B | 6 -> C | 7 -> N
|
||||
| 8 -> O | 9 -> F | 10 -> Ne | 11 -> Na
|
||||
| 12 -> Mg | 13 -> Al | 14 -> Si | 15 -> P
|
||||
| 16 -> S | 17 -> Cl | 18 -> Ar | 19 -> K
|
||||
| 20 -> Ca | 21 -> Sc | 22 -> Ti | 23 -> V
|
||||
| 24 -> Cr | 25 -> Mn | 26 -> Fe | 27 -> Co
|
||||
| 28 -> Ni | 29 -> Cu | 30 -> Zn | 31 -> Ga
|
||||
| 32 -> Ge | 33 -> As | 34 -> Se | 35 -> Br
|
||||
| 36 -> Kr | 37 -> Rb | 38 -> Sr | 39 -> Y
|
||||
| 40 -> Zr | 41 -> Nb | 42 -> Mo | 43 -> Tc
|
||||
| 44 -> Ru | 45 -> Rh | 46 -> Pd | 47 -> Ag
|
||||
| 48 -> Cd | 49 -> In | 50 -> Sn | 51 -> Sb
|
||||
| 52 -> Te | 53 -> I | 54 -> Xe | 78 -> Pt
|
||||
| x -> raise (ElementError ("Element of charge "^(string_of_int x)^" unknown"))
|
||||
|
||||
|
||||
let of_charge c =
|
||||
Charge.to_int c |> of_int
|
||||
|
||||
#+end_src
|
||||
|
||||
** Database information
|
||||
|
||||
#+begin_src ocaml :tangle (eval mli)
|
||||
val covalent_radius : t -> Non_negative_float.t
|
||||
val vdw_radius : t -> Non_negative_float.t
|
||||
val mass : t -> Mass.t
|
||||
val small_core : t -> int
|
||||
val large_core : t -> int
|
||||
#+end_src
|
||||
|
||||
| ~covalent_radius~ | Covalent radii of the elements, in atomic units |
|
||||
| ~vdw_radius~ | Van der Waals radii of the elements, in atomic units |
|
||||
| ~mass~ | Atomic mass of the elements, in atomic units) |
|
||||
| ~small_core~ | Number of electrons in the small core model (all except the outermost two shells) |
|
||||
| ~large_core~ | Number of electrons in the large core model (all except the outermost shell) |
|
||||
|
||||
#+begin_src ocaml :tangle (eval ml) :exports none
|
||||
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 | Rb -> 2.20 | Sr -> 1.95 | Y -> 1.90
|
||||
| Zr -> 1.75 | Nb -> 1.64 | Mo -> 1.54 | Tc -> 1.47
|
||||
| Ru -> 1.46 | Rh -> 1.42 | Pd -> 1.39 | Ag -> 1.45
|
||||
| Cd -> 1.44 | In -> 1.42 | Sn -> 1.39 | Sb -> 1.39
|
||||
| Te -> 1.38 | I -> 1.39 | Xe -> 1.40 | Pt -> 1.30
|
||||
in
|
||||
Constants.a0 *. (result x)
|
||||
|> Non_negative_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 | Rb -> 3.03 | Sr -> 2.49 | Y -> 0.
|
||||
| Zr -> 0. | Nb -> 0. | Mo -> 0. | Tc -> 0.
|
||||
| Ru -> 0. | Rh -> 0. | Pd -> 1.63 | Ag -> 1.72
|
||||
| Cd -> 1.58 | In -> 1.93 | Sn -> 2.17 | Sb -> 2.06
|
||||
| Te -> 2.06 | I -> 1.98 | Xe -> 2.16 | Pt -> 1.75
|
||||
in
|
||||
Constants.a0 *. (result x)
|
||||
|> Non_negative_float.of_float
|
||||
|
||||
|
||||
let mass c =
|
||||
begin
|
||||
match c with
|
||||
| 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 | Rb -> 85.4678 | Sr -> 87.62 | Y -> 88.90584
|
||||
| Zr -> 91.224 | Nb -> 92.90637 | Mo -> 95.95 | Tc -> 98.
|
||||
| Ru -> 101.07 | Rh -> 102.90550 | Pd -> 106.42 | Ag -> 107.8682
|
||||
| Cd -> 112.414 | In -> 114.818 | Sn -> 118.710 | Sb -> 121.760
|
||||
| Te -> 127.60 | I -> 126.90447 | Xe -> 131.293 | Pt -> 195.084
|
||||
end
|
||||
|> Mass.of_float
|
||||
|
||||
|
||||
let noble_gas =
|
||||
[ He ; Ne ; Ar ; Kr ; Xe ]
|
||||
|
||||
|
||||
let large_core t =
|
||||
let num = to_int t in
|
||||
let rec loop = function
|
||||
| gas :: rest ->
|
||||
if gas < num then
|
||||
gas
|
||||
else
|
||||
loop rest
|
||||
| [] -> 0
|
||||
in
|
||||
List.rev_map to_int noble_gas
|
||||
|> loop
|
||||
|
||||
|
||||
let small_core t =
|
||||
let num = to_int t in
|
||||
let rec loop = function
|
||||
| large :: small :: rest ->
|
||||
if large < num then
|
||||
small
|
||||
else
|
||||
loop (small :: rest)
|
||||
| small :: [] ->
|
||||
if small < num then
|
||||
small
|
||||
else
|
||||
0
|
||||
| [] -> 0
|
||||
in
|
||||
List.rev_map to_int noble_gas
|
||||
|> loop
|
||||
|
||||
#+end_src
|
||||
|
||||
** Printers
|
||||
|
||||
#+begin_src ocaml :tangle (eval mli)
|
||||
val pp : Format.formatter -> t -> unit
|
||||
val pp_long : Format.formatter -> t -> unit
|
||||
#+end_src
|
||||
|
||||
#+begin_src ocaml :tangle (eval ml) :exports none
|
||||
let pp ppf t =
|
||||
Format.fprintf ppf "@[%s@]" (to_string t)
|
||||
|
||||
let pp_long ppf t =
|
||||
Format.fprintf ppf "@[%s@]" (to_long_string t)
|
||||
#+end_src
|
||||
|
||||
|
|
@ -1,20 +1,16 @@
|
|||
(* [[file:~/QCaml/particles/electrons.org::*Type][Type:2]] *)
|
||||
(** Type *)
|
||||
type t = {
|
||||
n_alfa : int ;
|
||||
n_beta : int ;
|
||||
}
|
||||
(* Type:2 ends here *)
|
||||
|
||||
|
||||
|
||||
(* | ~make~ | ~make n_alfa n_beta~ |
|
||||
* | ~of_atoms~ | Creates the data relative to electrons for a molecular system described by ~Nuclei.t~ for a given total charge and spin multiplicity. | *)
|
||||
|
||||
|
||||
(* [[file:~/QCaml/particles/electrons.org::*Creation][Creation:2]] *)
|
||||
open Common
|
||||
|
||||
let make n_alfa n_beta =
|
||||
(** Creation *)
|
||||
|
||||
let make n_alfa n_beta =
|
||||
{ n_alfa ; n_beta }
|
||||
|
||||
|
||||
|
|
@ -31,18 +27,9 @@ let of_atoms ?multiplicity:(multiplicity=1) ?charge:(charge=0) nuclei =
|
|||
if multiplicity <> (n_alfa - n_beta)+1 then
|
||||
invalid_arg (__FILE__^": make");
|
||||
result
|
||||
(* Creation:2 ends here *)
|
||||
|
||||
|
||||
|
||||
(* | ~charge~ | Sum of the charges of the electrons |
|
||||
* | ~n_elec~ | Number of electrons |
|
||||
* | ~n_alfa~ | Number of alpha electrons |
|
||||
* | ~n_beta~ | Number of beta electrons |
|
||||
* | ~multiplicity~ | Spin multiplicity: $2S+1$ | *)
|
||||
|
||||
|
||||
(* [[file:~/QCaml/particles/electrons.org::*Access][Access:2]] *)
|
||||
(** Access *)
|
||||
let charge e =
|
||||
- (e.n_alfa + e.n_beta)
|
||||
|> Charge.of_int
|
||||
|
|
@ -54,9 +41,10 @@ let n_beta t = t.n_beta
|
|||
let n_elec t = t.n_alfa + t.n_beta
|
||||
|
||||
let multiplicity t = t.n_alfa - t.n_beta + 1
|
||||
(* Access:2 ends here *)
|
||||
|
||||
(* [[file:~/QCaml/particles/electrons.org::*Printers][Printers:2]] *)
|
||||
|
||||
(** Printers *)
|
||||
|
||||
let pp ppf t =
|
||||
Format.fprintf ppf "@[n_alfa=%d, n_beta=%d@]" t.n_alfa t.n_beta
|
||||
(* Printers:2 ends here *)
|
||||
|
||||
|
|
|
|||
|
|
@ -1,41 +1,43 @@
|
|||
(* Type
|
||||
*
|
||||
* #+NAME: types *)
|
||||
(** Data structure containing the number of alpha and beta electrons.
|
||||
*)
|
||||
|
||||
(* [[file:~/QCaml/particles/electrons.org::types][types]] *)
|
||||
type t
|
||||
(* types ends here *)
|
||||
|
||||
(* Creation *)
|
||||
|
||||
|
||||
(* [[file:~/QCaml/particles/electrons.org::*Creation][Creation:1]] *)
|
||||
open Common
|
||||
|
||||
(** Creation *)
|
||||
|
||||
val make : int -> int -> t
|
||||
(* make n_alfa n_beta *)
|
||||
|
||||
val of_atoms : ?multiplicity:int -> ?charge:int -> Nuclei.t -> t
|
||||
(* @param multiplicity default is 1
|
||||
(** Creates the data relative to electrons for a molecular system described by
|
||||
Nuclei.t for a given total charge and spin multiplicity.
|
||||
|
||||
@param multiplicity default is 1
|
||||
@param charge default is 0
|
||||
@raise Invalid_argument if the spin multiplicity is not compatible with
|
||||
the molecule and the total charge.
|
||||
*)
|
||||
(* Creation:1 ends here *)
|
||||
|
||||
(* Access *)
|
||||
(** Access *)
|
||||
|
||||
|
||||
(* [[file:~/QCaml/particles/electrons.org::*Access][Access:1]] *)
|
||||
val charge : t -> Charge.t
|
||||
(** Sum of the charges of the electrons *)
|
||||
|
||||
val n_elec : t -> int
|
||||
(* Number of electrons *)
|
||||
|
||||
val n_alfa : t -> int
|
||||
(* Number of alpha electrons *)
|
||||
|
||||
val n_beta : t -> int
|
||||
(* Number of beta electrons *)
|
||||
|
||||
val multiplicity : t -> int
|
||||
(* Access:1 ends here *)
|
||||
|
||||
(* Printers *)
|
||||
(* Spin multiplicity: 2S+1 *)
|
||||
|
||||
|
||||
(* [[file:~/QCaml/particles/electrons.org::*Printers][Printers:1]] *)
|
||||
(** Printers *)
|
||||
|
||||
val pp : Format.formatter -> t -> unit
|
||||
(* Printers:1 ends here *)
|
||||
|
|
|
|||
|
|
@ -1,4 +1,5 @@
|
|||
(* [[file:~/QCaml/particles/element.org::*Type][Type:2]] *)
|
||||
(** Type *)
|
||||
|
||||
type t =
|
||||
|X
|
||||
|H |He
|
||||
|
|
@ -6,43 +7,16 @@ type t =
|
|||
|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
|
||||
|Rb|Sr|Y |Zr|Nb|Mo|Tc|Ru|Rh|Pd|Ag|Cd|In|Sn|Sb|Te|I |Xe
|
||||
|Pt
|
||||
|Pt
|
||||
|
||||
exception ElementError of string
|
||||
|
||||
open Common
|
||||
(* Type:2 ends here *)
|
||||
|
||||
|
||||
(** Conversion *)
|
||||
|
||||
(* | ~of_string~ | Creates an ~Element.t~ from a chemical symbol or from the full name of the element (case insensitive) |
|
||||
* | ~to_string~ | Gets the chemical symbol of the ~Element.t~ in a string |
|
||||
* | ~to_long_string~ | Gets the full name of the ~Element.t~ in a string |
|
||||
* | ~to_int~ | Convert to the atomic charge, with ~int~ type |
|
||||
* | ~of_int~ | Create from the atomic charge, with ~int~ type |
|
||||
* | ~to_charge~ | Convert to the atomic charge, with ~Charge.t~ type |
|
||||
* | ~of_charge~ | Create from the atomic charge, with ~Charge.t~ type |
|
||||
*
|
||||
* #+begin_example
|
||||
* Element.of_string "Fe" ;;
|
||||
* - : Element.t = Particles.Element.Fe
|
||||
*
|
||||
* Element.of_string "hydrogen" ;;
|
||||
* - : Element.t = Particles.Element.H
|
||||
*
|
||||
* Element.of_string "Kryptonite" ;;
|
||||
* Exception: Particles.Element.ElementError "Element Kryptonite unknown".
|
||||
*
|
||||
* Element.(to_long_string Fe) ;;
|
||||
* - : string = "Iron"
|
||||
*
|
||||
* Element.(to_string Fe);;
|
||||
* - : string = "Fe"
|
||||
* #+end_example *)
|
||||
|
||||
|
||||
(* [[file:~/QCaml/particles/element.org::*Conversion][Conversion:2]] *)
|
||||
let of_string x =
|
||||
let of_string x =
|
||||
match (String.capitalize_ascii (String.lowercase_ascii x)) with
|
||||
| "X" | "Dummy" -> X | "H" | "Hydrogen" -> H
|
||||
| "He" | "Helium" -> He | "Li" | "Lithium" -> Li
|
||||
|
|
@ -130,7 +104,7 @@ let to_int = function
|
|||
| Te -> 52 | I -> 53 | Xe -> 54 | Pt -> 78
|
||||
|
||||
|
||||
let to_charge c =
|
||||
let to_charge c =
|
||||
to_int c |> Charge.of_int
|
||||
|
||||
|
||||
|
|
@ -154,19 +128,11 @@ let of_int = function
|
|||
|
||||
let of_charge c =
|
||||
Charge.to_int c |> of_int
|
||||
(* Conversion:2 ends here *)
|
||||
|
||||
|
||||
(** Database information *)
|
||||
|
||||
(* | ~covalent_radius~ | Covalent radii of the elements, in atomic units |
|
||||
* | ~vdw_radius~ | Van der Waals radii of the elements, in atomic units |
|
||||
* | ~mass~ | Atomic mass of the elements, in atomic units) |
|
||||
* | ~small_core~ | Number of electrons in the small core model (all except the outermost two shells) |
|
||||
* | ~large_core~ | Number of electrons in the large core model (all except the outermost shell) | *)
|
||||
|
||||
|
||||
(* [[file:~/QCaml/particles/element.org::*Database information][Database information:2]] *)
|
||||
let covalent_radius x =
|
||||
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
|
||||
|
|
@ -187,7 +153,7 @@ let covalent_radius x =
|
|||
|> Non_negative_float.of_float
|
||||
|
||||
|
||||
let vdw_radius x =
|
||||
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
|
||||
|
|
@ -233,7 +199,7 @@ let noble_gas =
|
|||
[ He ; Ne ; Ar ; Kr ; Xe ]
|
||||
|
||||
|
||||
let large_core t =
|
||||
let large_core t =
|
||||
let num = to_int t in
|
||||
let rec loop = function
|
||||
| gas :: rest ->
|
||||
|
|
@ -247,7 +213,7 @@ let large_core t =
|
|||
|> loop
|
||||
|
||||
|
||||
let small_core t =
|
||||
let small_core t =
|
||||
let num = to_int t in
|
||||
let rec loop = function
|
||||
| large :: small :: rest ->
|
||||
|
|
@ -264,12 +230,13 @@ let small_core t =
|
|||
in
|
||||
List.rev_map to_int noble_gas
|
||||
|> loop
|
||||
(* Database information:2 ends here *)
|
||||
|
||||
(* [[file:~/QCaml/particles/element.org::*Printers][Printers:2]] *)
|
||||
|
||||
(** Printers *)
|
||||
|
||||
let pp ppf t =
|
||||
Format.fprintf ppf "@[%s@]" (to_string t)
|
||||
|
||||
let pp_long ppf t =
|
||||
Format.fprintf ppf "@[%s@]" (to_long_string t)
|
||||
(* Printers:2 ends here *)
|
||||
|
||||
|
|
|
|||
|
|
@ -1,8 +1,5 @@
|
|||
(* Type
|
||||
*
|
||||
* #+NAME: types *)
|
||||
(** Chemical elements *)
|
||||
|
||||
(* [[file:~/QCaml/particles/element.org::types][types]] *)
|
||||
type t =
|
||||
|X
|
||||
|H |He
|
||||
|
|
@ -10,43 +7,73 @@ type t =
|
|||
|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
|
||||
|Rb|Sr|Y |Zr|Nb|Mo|Tc|Ru|Rh|Pd|Ag|Cd|In|Sn|Sb|Te|I |Xe
|
||||
|Pt
|
||||
|Pt
|
||||
|
||||
exception ElementError of string
|
||||
|
||||
open Common
|
||||
(* types ends here *)
|
||||
|
||||
(* Conversion *)
|
||||
|
||||
|
||||
(* [[file:~/QCaml/particles/element.org::*Conversion][Conversion:1]] *)
|
||||
(** Conversion *)
|
||||
|
||||
(**
|
||||
* Element.of_string "Fe" ;;
|
||||
* - : Element.t = Particles.Element.Fe
|
||||
*
|
||||
* Element.of_string "hydrogen" ;;
|
||||
* - : Element.t = Particles.Element.H
|
||||
*
|
||||
* Element.of_string "Kryptonite" ;;
|
||||
* Exception: Particles.Element.ElementError "Element Kryptonite unknown".
|
||||
*
|
||||
* Element.(to_long_string Fe) ;;
|
||||
* - : string = "Iron"
|
||||
*
|
||||
* Element.(to_string Fe);;
|
||||
* - : string = "Fe"
|
||||
*)
|
||||
|
||||
val of_string : string -> t
|
||||
val to_string : t -> string
|
||||
val to_long_string : t -> string
|
||||
(** Creates an ~Element.t~ from a chemical symbol or from the full name of the element (case insensitive) *)
|
||||
|
||||
val to_string : t -> string
|
||||
(** Gets the chemical symbol of the ~Element.t~ in a string *)
|
||||
|
||||
val to_long_string : t -> string
|
||||
(** Gets the full name of the ~Element.t~ in a string *)
|
||||
|
||||
val to_int : t -> int
|
||||
(** Convert to the atomic charge, with ~int~ type *)
|
||||
|
||||
val to_int : t -> int
|
||||
val of_int : int -> t
|
||||
(** Create from the atomic charge, with ~int~ type *)
|
||||
|
||||
val to_charge : t -> Charge.t
|
||||
(** Convert to the atomic charge, with ~Charge.t~ type *)
|
||||
|
||||
val of_charge : Charge.t -> t
|
||||
(* Conversion:1 ends here *)
|
||||
|
||||
(* Database information *)
|
||||
(** Create from the atomic charge, with ~Charge.t~ type *)
|
||||
|
||||
|
||||
(* [[file:~/QCaml/particles/element.org::*Database information][Database information:1]] *)
|
||||
(** Database information *)
|
||||
|
||||
val covalent_radius : t -> Non_negative_float.t
|
||||
(** Covalent radii of the elements, in atomic units *)
|
||||
|
||||
val vdw_radius : t -> Non_negative_float.t
|
||||
(** Van der Waals radii of the elements, in atomic units *)
|
||||
|
||||
val mass : t -> Mass.t
|
||||
(** Atomic mass of the elements, in atomic units) *)
|
||||
|
||||
val small_core : t -> int
|
||||
(** Number of electrons in the small core model (all except the outermost two shells) *)
|
||||
|
||||
val large_core : t -> int
|
||||
(* Database information:1 ends here *)
|
||||
|
||||
(* Printers *)
|
||||
(** Number of electrons in the large core model (all except the outermost shell) *)
|
||||
|
||||
|
||||
(* [[file:~/QCaml/particles/element.org::*Printers][Printers:1]] *)
|
||||
(** Printers *)
|
||||
|
||||
val pp : Format.formatter -> t -> unit
|
||||
val pp_long : Format.formatter -> t -> unit
|
||||
(* Printers:1 ends here *)
|
||||
|
|
|
|||
|
|
@ -1,3 +1,3 @@
|
|||
(* [[file:~/QCaml/particles/mass.org::*Atomic mass][Atomic mass:2]] *)
|
||||
(** Atomic mass *)
|
||||
|
||||
include Common.Non_negative_float
|
||||
(* Atomic mass:2 ends here *)
|
||||
|
|
|
|||
|
|
@ -1,12 +1,3 @@
|
|||
(* Atomic mass
|
||||
* :PROPERTIES:
|
||||
* :header-args: :noweb yes :comments both
|
||||
* :END:
|
||||
*
|
||||
* Atomic mass, a non-negative float.
|
||||
*
|
||||
* #+NAME: types *)
|
||||
(** Atomic mass *)
|
||||
|
||||
(* [[file:~/QCaml/particles/mass.org::types][types]] *)
|
||||
include module type of Common.Non_negative_float
|
||||
(* types ends here *)
|
||||
|
|
|
|||
|
|
@ -1,21 +1,12 @@
|
|||
(* [[file:~/QCaml/particles/nuclei.org::*Type][Type:2]] *)
|
||||
(** Type *)
|
||||
open Common
|
||||
|
||||
type t = (Element.t * Coordinate.t) array
|
||||
open Xyz_ast
|
||||
(* Type:2 ends here *)
|
||||
|
||||
|
||||
(** Conversion *)
|
||||
|
||||
(* | ~of_xyz_string~ | Create from a string, in xyz format |
|
||||
* | ~of_xyz_file~ | Create from a file, in xyz format |
|
||||
* | ~of_zmt_string~ | Create from a string, in z-matrix format |
|
||||
* | ~of_zmt_file~ | Create from a file, in z-matrix format |
|
||||
* | ~to_string~ | Transform to a string, for printing |
|
||||
* | ~of_filename~ | Detects the type of file (xyz, z-matrix) and reads the file | *)
|
||||
|
||||
|
||||
(* [[file:~/QCaml/particles/nuclei.org::*Conversion][Conversion:2]] *)
|
||||
let of_xyz_lexbuf lexbuf =
|
||||
let data =
|
||||
Xyz_parser.input Nuclei_lexer.read_all lexbuf
|
||||
|
|
@ -115,18 +106,10 @@ let to_xyz_string t =
|
|||
) t
|
||||
|> Array.to_list )
|
||||
|> String.concat "\n"
|
||||
(* Conversion:2 ends here *)
|
||||
|
||||
|
||||
(** Query *)
|
||||
|
||||
(* | ~formula~ | Returns the chemical formula |
|
||||
* | ~repulsion~ | Nuclear repulsion energy, in atomic units |
|
||||
* | ~charge~ | Sum of the charges of the nuclei |
|
||||
* | ~small_core~ | Number of core electrons in the small core model |
|
||||
* | ~large_core~ | Number of core electrons in the large core model | *)
|
||||
|
||||
|
||||
(* [[file:~/QCaml/particles/nuclei.org::*Query][Query:2]] *)
|
||||
let formula t =
|
||||
let dict = Hashtbl.create 67 in
|
||||
Array.iter (fun (e,_) ->
|
||||
|
|
@ -174,9 +157,10 @@ let small_core a =
|
|||
|
||||
let large_core a =
|
||||
Array.fold_left (fun accu (e,_) -> accu + (Element.large_core e)) 0 a
|
||||
(* Query:2 ends here *)
|
||||
|
||||
(* [[file:~/QCaml/particles/nuclei.org::*Read][Read:2]] *)
|
||||
|
||||
(** Read *)
|
||||
|
||||
let of_trexio f =
|
||||
let num = Trexio.read_nucleus_num f in
|
||||
let charge = Trexio.read_nucleus_charge f
|
||||
|
|
@ -188,9 +172,10 @@ let of_trexio f =
|
|||
z = coord.(3*i+2) } in
|
||||
(Element.of_charge charge.(i), Coordinate.make coord)
|
||||
)
|
||||
(* Read:2 ends here *)
|
||||
|
||||
(* [[file:~/QCaml/particles/nuclei.org::*Write][Write:2]] *)
|
||||
|
||||
(** Write *)
|
||||
|
||||
let to_trexio f t =
|
||||
let num = Array.length t in
|
||||
Trexio.write_nucleus_num f num;
|
||||
|
|
@ -210,9 +195,10 @@ let to_trexio f t =
|
|||
|
||||
repulsion t
|
||||
|> Trexio.write_nucleus_repulsion f
|
||||
(* Write:2 ends here *)
|
||||
|
||||
(* [[file:~/QCaml/particles/nuclei.org::*Printers][Printers:2]] *)
|
||||
|
||||
(** Printers *)
|
||||
|
||||
let pp ppf t =
|
||||
Format.fprintf ppf "@[%s@]" (to_string t)
|
||||
(* Printers:2 ends here *)
|
||||
|
||||
|
|
|
|||
|
|
@ -1,58 +1,64 @@
|
|||
(* Type
|
||||
* <<<~Nuclei.t~>>>
|
||||
*
|
||||
* #+NAME: types *)
|
||||
(** Type *)
|
||||
|
||||
(* [[file:~/QCaml/particles/nuclei.org::types][types]] *)
|
||||
open Common
|
||||
|
||||
type t = (Element.t * Coordinate.t) array
|
||||
(* types ends here *)
|
||||
|
||||
(* Conversion *)
|
||||
(** Conversion *)
|
||||
|
||||
|
||||
(* [[file:~/QCaml/particles/nuclei.org::*Conversion][Conversion:1]] *)
|
||||
val of_xyz_string : string -> t
|
||||
(** Create from a string, in xyz format *)
|
||||
|
||||
val to_xyz_string : t -> string
|
||||
(** Transform into a string, in xyz format *)
|
||||
|
||||
val of_xyz_file : string -> t
|
||||
(** Create from a file, in xyz format *)
|
||||
|
||||
val of_zmt_string : string -> t
|
||||
(** Create from a string, in z-matrix format *)
|
||||
|
||||
val of_zmt_file : string -> t
|
||||
(** Create from a file, in z-matrix format *)
|
||||
|
||||
val to_string : t -> string
|
||||
(** Transform to a string, for printing *)
|
||||
|
||||
val of_filename : string -> t
|
||||
(* Conversion:1 ends here *)
|
||||
|
||||
(* Query *)
|
||||
(** Detects the type of file (xyz, z-matrix) and reads the file *)
|
||||
|
||||
|
||||
(* [[file:~/QCaml/particles/nuclei.org::*Query][Query:1]] *)
|
||||
(** Query *)
|
||||
|
||||
val formula : t -> string
|
||||
(** Returns the chemical formula *)
|
||||
|
||||
val repulsion : t -> float
|
||||
(** Nuclear repulsion energy, in atomic units *)
|
||||
|
||||
val charge : t -> Charge.t
|
||||
(** Sum of the charges of the nuclei *)
|
||||
|
||||
val small_core : t -> int
|
||||
(** Number of core electrons in the small core model *)
|
||||
|
||||
val large_core : t -> int
|
||||
(* Query:1 ends here *)
|
||||
|
||||
(* Read *)
|
||||
(** Number of core electrons in the large core model | *)
|
||||
|
||||
|
||||
(* [[file:~/QCaml/particles/nuclei.org::*Read][Read:1]] *)
|
||||
(** Read *)
|
||||
|
||||
val of_trexio : Trexio.trexio_file -> t
|
||||
(* Read:1 ends here *)
|
||||
(** Read from a file in TREXIO format *)
|
||||
|
||||
(* Write *)
|
||||
(** Write *)
|
||||
|
||||
|
||||
(* [[file:~/QCaml/particles/nuclei.org::*Write][Write:1]] *)
|
||||
val to_trexio : Trexio.trexio_file -> t -> unit
|
||||
(* Write:1 ends here *)
|
||||
|
||||
(* Printers *)
|
||||
(** Write to a file in TREXIO format *)
|
||||
|
||||
|
||||
(* [[file:~/QCaml/particles/nuclei.org::*Printers][Printers:1]] *)
|
||||
(** Printers *)
|
||||
|
||||
val pp : Format.formatter -> t -> unit
|
||||
(* Printers:1 ends here *)
|
||||
|
||||
|
|
|
|||
|
|
@ -1,10 +1,5 @@
|
|||
(* Lexer
|
||||
*
|
||||
* =nuclei_lexer.mll= contains the description of the lexemes used in
|
||||
* an xyz file. *)
|
||||
(** Lexer to read xyz files *)
|
||||
|
||||
|
||||
(* [[file:~/QCaml/particles/nuclei.org::*Lexer][Lexer:1]] *)
|
||||
{
|
||||
open Xyz_parser
|
||||
}
|
||||
|
|
@ -45,4 +40,3 @@ rule read_all = parse
|
|||
done;
|
||||
*)
|
||||
}
|
||||
(* Lexer:1 ends here *)
|
||||
|
|
|
|||
|
|
@ -1,10 +1,9 @@
|
|||
(** When an xyz file is read by =xyz_parser.mly=, it is converted into
|
||||
* an ~xyz_file~ data structure. *)
|
||||
|
||||
|
||||
(* When an xyz file is read by =xyz_parser.mly=, it is converted into
|
||||
* an ~xyz_file~ data structure. *)
|
||||
(** Parser *)
|
||||
|
||||
|
||||
(* [[file:~/QCaml/particles/nuclei.org::*Parser][Parser:2]] *)
|
||||
open Common
|
||||
|
||||
type nucleus =
|
||||
|
|
@ -19,4 +18,4 @@ type xyz_file =
|
|||
file_title : string ;
|
||||
nuclei : nucleus list ;
|
||||
}
|
||||
(* Parser:2 ends here *)
|
||||
|
||||
|
|
|
|||
|
|
@ -1,4 +1,5 @@
|
|||
(* [[file:~/QCaml/particles/zmatrix.org::*Type][Type:2]] *)
|
||||
(** Type *)
|
||||
|
||||
module StringMap = Map.Make(String)
|
||||
|
||||
type atom_id = int
|
||||
|
|
@ -6,7 +7,7 @@ type angle = Label of string | Value of float
|
|||
type distance = Label of string | Value of float
|
||||
type dihedral = Label of string | Value of float
|
||||
|
||||
type line =
|
||||
type line =
|
||||
| First of Element.t
|
||||
| Second of (Element.t * distance)
|
||||
| Third of (Element.t * atom_id * distance * atom_id * angle)
|
||||
|
|
@ -14,59 +15,10 @@ type line =
|
|||
| Coord of (string * float)
|
||||
|
||||
type t = (line array * float StringMap.t)
|
||||
(* Type:2 ends here *)
|
||||
|
||||
|
||||
(** Conversion *)
|
||||
|
||||
(* | ~of_string~ | Reads a z-matrix from a string |
|
||||
* | ~to_xyz~ | Converts to xyz format, as in the ~Nuclei~ module |
|
||||
* | ~to_xyz_string~ | Converts to xyz format, as a string |
|
||||
*
|
||||
* #+begin_example
|
||||
* let zmt = Zmatrix.of_string "
|
||||
* n
|
||||
* n 1 nn
|
||||
* h 1 hn 2 hnn
|
||||
* h 2 hn 1 hnn 3 dih4
|
||||
* h 1 hn 2 hnn 4 dih5
|
||||
* h 2 hn 1 hnn 3 dih5
|
||||
*
|
||||
* nn 1.446
|
||||
* hn 1.016
|
||||
* hnn 106.0
|
||||
* dih4 -54.38
|
||||
* dih5 54.38
|
||||
* " ;;
|
||||
* - : Zmatrix.t = N
|
||||
* N 1 1.446000
|
||||
* H 1 1.016000 2 106.000000
|
||||
* H 2 1.016000 1 106.000000 3 -54.380000
|
||||
* H 1 1.016000 2 106.000000 4 54.380000
|
||||
* H 2 1.016000 1 106.000000 3 54.380000
|
||||
*
|
||||
*
|
||||
* Zmatrix.to_xyz zmt ;;
|
||||
* - : (Element.t * float * float * float) array =
|
||||
* [|(N, 0., 0., 0.); (N, 0., 0., 1.446);
|
||||
* (H, -0.976641883073332107, 0., -0.280047553510071046);
|
||||
* (H, -0.568802835186988709, 0.793909757123734683, 1.726047553510071);
|
||||
* (H, 0.314092649983635563, 0.924756819385119, -0.280047553510071101);
|
||||
* (H, -0.568802835186988709, -0.793909757123734683, 1.726047553510071)|]
|
||||
*
|
||||
*
|
||||
* Zmatrix.to_xyz_string zmt ;;
|
||||
* - : string =
|
||||
* "N 0.000000 0.000000 0.000000
|
||||
* N 0.000000 0.000000 1.446000
|
||||
* H -0.976642 0.000000 -0.280048
|
||||
* H -0.568803 0.793910 1.726048
|
||||
* H 0.314093 0.924757 -0.280048
|
||||
* H -0.568803 -0.793910 1.726048"
|
||||
* #+end_example *)
|
||||
|
||||
|
||||
|
||||
(* [[file:~/QCaml/particles/zmatrix.org::*Conversion][Conversion:2]] *)
|
||||
let pi = Common.Constants.pi
|
||||
let to_radian = pi /. 180.
|
||||
|
||||
|
|
@ -78,43 +30,43 @@ let rec in_range (xmin, xmax) x =
|
|||
else
|
||||
x
|
||||
|
||||
let atom_id_of_int : int -> atom_id =
|
||||
let atom_id_of_int : int -> atom_id =
|
||||
fun x -> ( assert (x>0) ; x)
|
||||
|
||||
let distance_of_float : float -> distance =
|
||||
let distance_of_float : float -> distance =
|
||||
fun x -> ( assert (x>=0.) ; Value x)
|
||||
|
||||
let angle_of_float : float -> angle =
|
||||
let angle_of_float : float -> angle =
|
||||
fun x -> Value (in_range (-180., 180.) x)
|
||||
|
||||
let dihedral_of_float : float -> dihedral =
|
||||
let dihedral_of_float : float -> dihedral =
|
||||
fun x -> Value (in_range (-360., 360.) x)
|
||||
|
||||
|
||||
let atom_id_of_string : string -> atom_id =
|
||||
let atom_id_of_string : string -> atom_id =
|
||||
fun i -> atom_id_of_int @@ int_of_string i
|
||||
|
||||
let distance_of_string : string -> distance =
|
||||
fun s ->
|
||||
fun s ->
|
||||
try
|
||||
distance_of_float @@ float_of_string s
|
||||
with _ -> Label s
|
||||
|
||||
let angle_of_string : string -> angle =
|
||||
fun s ->
|
||||
fun s ->
|
||||
try
|
||||
angle_of_float @@ float_of_string s
|
||||
with _ -> Label s
|
||||
|
||||
let dihedral_of_string : string -> dihedral =
|
||||
fun s ->
|
||||
fun s ->
|
||||
try
|
||||
dihedral_of_float @@ float_of_string s
|
||||
with _ -> Label s
|
||||
|
||||
let int_of_atom_id : atom_id -> int = fun x -> x
|
||||
|
||||
let float_of_distance : float StringMap.t -> distance -> float =
|
||||
let float_of_distance : float StringMap.t -> distance -> float =
|
||||
fun map -> function
|
||||
| Value x -> x
|
||||
| Label s -> StringMap.find s map
|
||||
|
|
@ -130,7 +82,7 @@ let float_of_dihedral : float StringMap.t -> dihedral -> float =
|
|||
| Label s -> StringMap.find s map
|
||||
|
||||
|
||||
let string_of_line map =
|
||||
let string_of_line map =
|
||||
let f_r = float_of_distance map
|
||||
and f_a = float_of_angle map
|
||||
and f_d = float_of_dihedral map
|
||||
|
|
@ -153,7 +105,7 @@ let line_of_string l =
|
|||
| e :: _ :: r :: [] -> Second
|
||||
(Element.of_string e,
|
||||
distance_of_string r)
|
||||
| e :: i :: r :: j :: a :: [] -> Third
|
||||
| e :: i :: r :: j :: a :: [] -> Third
|
||||
(Element.of_string e,
|
||||
atom_id_of_string i,
|
||||
distance_of_string r,
|
||||
|
|
@ -179,7 +131,7 @@ let of_string t =
|
|||
|> List.map line_of_string
|
||||
in
|
||||
|
||||
let l =
|
||||
let l =
|
||||
match l with
|
||||
| First _ :: Second _ :: Third _ :: _
|
||||
| First _ :: Second _ :: Coord _ :: []
|
||||
|
|
@ -200,7 +152,7 @@ let of_string t =
|
|||
work [] (StringMap.empty) l
|
||||
in
|
||||
(Array.of_list l, m)
|
||||
|
||||
|
||||
|
||||
(** Linear algebra *)
|
||||
|
||||
|
|
@ -225,17 +177,17 @@ let normalized u =
|
|||
let cross (x,y,z) (x',y',z') =
|
||||
((y *. z' -. z *. y'), -. (x *. z' -. z *. x'), (x *. y' -. y *. x'))
|
||||
|
||||
let rotation_matrix axis angle =
|
||||
let rotation_matrix axis angle =
|
||||
(* Euler-Rodrigues formula for rotation matrix, taken from
|
||||
https://github.com/jevandezande/zmatrix/blob/master/converter.py
|
||||
*)
|
||||
let a =
|
||||
let a =
|
||||
(cos (angle *. to_radian *. 0.5))
|
||||
in
|
||||
let (b, c, d) =
|
||||
let (b, c, d) =
|
||||
(-. sin (angle *. to_radian *. 0.5)) |. (normalized axis)
|
||||
in
|
||||
Array.of_list @@
|
||||
Array.of_list @@
|
||||
[(a *. a +. b *. b -. c *. c -. d *. d,
|
||||
2. *. (b *. c -. a *. d),
|
||||
2. *. (b *. d +. a *. c));
|
||||
|
|
@ -245,16 +197,16 @@ let rotation_matrix axis angle =
|
|||
(2. *. (b *. d -. a *. c),
|
||||
2. *. (c *. d +. a *. b),
|
||||
a *. a +. d *. d -. b *. b -. c *. c)]
|
||||
|
||||
|
||||
|
||||
let apply_rotation_matrix rot u =
|
||||
(dot rot.(0) u, dot rot.(1) u, dot rot.(2) u)
|
||||
|
||||
|
||||
|
||||
let to_xyz (z,map) =
|
||||
let result =
|
||||
Array.make (Array.length z) None
|
||||
in
|
||||
in
|
||||
|
||||
let get_cartesian_coord i =
|
||||
match result.(i-1) with
|
||||
|
|
@ -265,14 +217,14 @@ let to_xyz (z,map) =
|
|||
|
||||
let append_line i' =
|
||||
match z.(i') with
|
||||
| First e ->
|
||||
| First e ->
|
||||
result.(i') <- Some (e, 0., 0., 0.)
|
||||
| Second (e, r) ->
|
||||
| Second (e, r) ->
|
||||
let r =
|
||||
float_of_distance map r
|
||||
in
|
||||
result.(i') <- Some (e, 0., 0., r)
|
||||
| Third (e, i, r, j, a) ->
|
||||
| Third (e, i, r, j, a) ->
|
||||
begin
|
||||
let i, r, j, a =
|
||||
int_of_atom_id i,
|
||||
|
|
@ -281,13 +233,13 @@ let to_xyz (z,map) =
|
|||
float_of_angle map a
|
||||
in
|
||||
let ui, uj =
|
||||
get_cartesian_coord i,
|
||||
get_cartesian_coord i,
|
||||
get_cartesian_coord j
|
||||
in
|
||||
let u_ij =
|
||||
let u_ij =
|
||||
(uj |- ui)
|
||||
in
|
||||
let rot =
|
||||
let rot =
|
||||
rotation_matrix (0., 1., 0.) a
|
||||
in
|
||||
let new_vec =
|
||||
|
|
@ -298,7 +250,7 @@ let to_xyz (z,map) =
|
|||
in
|
||||
result.(i') <- Some (e, x, y, z)
|
||||
end
|
||||
| Other (e, i, r, j, a, k, d) ->
|
||||
| Other (e, i, r, j, a, k, d) ->
|
||||
begin
|
||||
let i, r, j, a, k, d =
|
||||
int_of_atom_id i,
|
||||
|
|
@ -309,7 +261,7 @@ let to_xyz (z,map) =
|
|||
float_of_dihedral map d
|
||||
in
|
||||
let ui, uj, uk =
|
||||
get_cartesian_coord i,
|
||||
get_cartesian_coord i,
|
||||
get_cartesian_coord j,
|
||||
get_cartesian_coord k
|
||||
in
|
||||
|
|
@ -319,7 +271,7 @@ let to_xyz (z,map) =
|
|||
let normal =
|
||||
cross u_ij u_kj
|
||||
in
|
||||
let new_vec =
|
||||
let new_vec =
|
||||
r |. (normalized u_ij)
|
||||
|> apply_rotation_matrix (rotation_matrix normal a)
|
||||
|> apply_rotation_matrix (rotation_matrix u_ij d)
|
||||
|
|
@ -332,7 +284,7 @@ let to_xyz (z,map) =
|
|||
| Coord _ -> ()
|
||||
in
|
||||
Array.iteri (fun i _ -> append_line i) z;
|
||||
let result =
|
||||
let result =
|
||||
Array.map (function
|
||||
| Some x -> x
|
||||
| None -> failwith "Some atoms were not defined" ) result
|
||||
|
|
@ -341,15 +293,16 @@ let to_xyz (z,map) =
|
|||
|
||||
|
||||
let to_xyz_string (l,map) =
|
||||
String.concat "\n"
|
||||
( to_xyz (l,map)
|
||||
|> Array.map (fun (e,x,y,z) ->
|
||||
String.concat "\n"
|
||||
( to_xyz (l,map)
|
||||
|> Array.map (fun (e,x,y,z) ->
|
||||
Printf.sprintf "%s %f %f %f" (Element.to_string e) x y z)
|
||||
|> Array.to_list
|
||||
)
|
||||
(* Conversion:2 ends here *)
|
||||
|
||||
(* [[file:~/QCaml/particles/zmatrix.org::*Printers][Printers:2]] *)
|
||||
|
||||
(** Printers *)
|
||||
|
||||
let pp ppf (a, map) =
|
||||
let f = string_of_line map in
|
||||
Format.fprintf ppf "@[";
|
||||
|
|
@ -357,4 +310,4 @@ let pp ppf (a, map) =
|
|||
Format.fprintf ppf "%s@." (f line)
|
||||
) a;
|
||||
Format.fprintf ppf "@]"
|
||||
(* Printers:2 ends here *)
|
||||
|
||||
|
|
|
|||
|
|
@ -1,23 +1,61 @@
|
|||
(* Type
|
||||
*
|
||||
* #+NAME: types *)
|
||||
(** Z-matrix representation of nuclear coordinates *)
|
||||
|
||||
(** let zmt = Zmatrix.of_string "
|
||||
* n
|
||||
* n 1 nn
|
||||
* h 1 hn 2 hnn
|
||||
* h 2 hn 1 hnn 3 dih4
|
||||
* h 1 hn 2 hnn 4 dih5
|
||||
* h 2 hn 1 hnn 3 dih5
|
||||
*
|
||||
* nn 1.446
|
||||
* hn 1.016
|
||||
* hnn 106.0
|
||||
* dih4 -54.38
|
||||
* dih5 54.38
|
||||
* " ;;
|
||||
* - : Zmatrix.t = N
|
||||
* N 1 1.446000
|
||||
* H 1 1.016000 2 106.000000
|
||||
* H 2 1.016000 1 106.000000 3 -54.380000
|
||||
* H 1 1.016000 2 106.000000 4 54.380000
|
||||
* H 2 1.016000 1 106.000000 3 54.380000
|
||||
*
|
||||
*
|
||||
* Zmatrix.to_xyz zmt ;;
|
||||
* - : (Element.t * float * float * float) array =
|
||||
* [|(N, 0., 0., 0.); (N, 0., 0., 1.446);
|
||||
* (H, -0.976641883073332107, 0., -0.280047553510071046);
|
||||
* (H, -0.568802835186988709, 0.793909757123734683, 1.726047553510071);
|
||||
* (H, 0.314092649983635563, 0.924756819385119, -0.280047553510071101);
|
||||
* (H, -0.568802835186988709, -0.793909757123734683, 1.726047553510071)|]
|
||||
*
|
||||
*
|
||||
* Zmatrix.to_xyz_string zmt ;;
|
||||
* - : string =
|
||||
* "N 0.000000 0.000000 0.000000
|
||||
* N 0.000000 0.000000 1.446000
|
||||
* H -0.976642 0.000000 -0.280048
|
||||
* H -0.568803 0.793910 1.726048
|
||||
* H 0.314093 0.924757 -0.280048
|
||||
* H -0.568803 -0.793910 1.726048"
|
||||
*)
|
||||
|
||||
(* [[file:~/QCaml/particles/zmatrix.org::types][types]] *)
|
||||
type t
|
||||
(* types ends here *)
|
||||
|
||||
(* Conversion *)
|
||||
|
||||
|
||||
(* [[file:~/QCaml/particles/zmatrix.org::*Conversion][Conversion:1]] *)
|
||||
(** Conversion *)
|
||||
|
||||
val of_string : string -> t
|
||||
(** Reads a z-matrix from a string *)
|
||||
|
||||
val to_xyz : t -> (Element.t * float * float * float) array
|
||||
(** Converts to xyz format, as in the ~Nuclei~ module *)
|
||||
|
||||
val to_xyz_string : t -> string
|
||||
(* Conversion:1 ends here *)
|
||||
|
||||
(* Printers *)
|
||||
(** Converts to xyz format, as a string *)
|
||||
|
||||
|
||||
(* [[file:~/QCaml/particles/zmatrix.org::*Printers][Printers:1]] *)
|
||||
(** Printers *)
|
||||
|
||||
val pp : Format.formatter -> t -> unit
|
||||
(* Printers:1 ends here *)
|
||||
|
|
|
|||
|
|
@ -1,27 +0,0 @@
|
|||
#+begin_src elisp tangle: no :results none :exports none
|
||||
(setq pwd (file-name-directory buffer-file-name))
|
||||
(setq name (file-name-nondirectory (substring buffer-file-name 0 -4)))
|
||||
(setq lib (concat pwd "lib/"))
|
||||
(setq testdir (concat pwd "test/"))
|
||||
(setq mli (concat lib name ".mli"))
|
||||
(setq ml (concat lib name ".ml"))
|
||||
(setq test-ml (concat testdir name ".ml"))
|
||||
(org-babel-tangle)
|
||||
#+end_src
|
||||
|
||||
* Atomic mass
|
||||
:PROPERTIES:
|
||||
:header-args: :noweb yes :comments both
|
||||
:END:
|
||||
|
||||
Atomic mass, a non-negative float.
|
||||
|
||||
#+NAME: types
|
||||
#+begin_src ocaml :tangle (eval mli)
|
||||
include module type of Common.Non_negative_float
|
||||
#+end_src
|
||||
|
||||
#+begin_src ocaml :tangle (eval ml) :exports none
|
||||
include Common.Non_negative_float
|
||||
#+end_src
|
||||
|
||||
|
|
@ -1,471 +0,0 @@
|
|||
#+begin_src elisp tangle: no :results none :exports none
|
||||
(setq pwd (file-name-directory buffer-file-name))
|
||||
(setq name (file-name-nondirectory (substring buffer-file-name 0 -4)))
|
||||
(setq lib (concat pwd "lib/"))
|
||||
(setq testdir (concat pwd "test/"))
|
||||
(setq mli (concat lib name ".mli"))
|
||||
(setq ml (concat lib name ".ml"))
|
||||
(setq test-ml (concat testdir name ".ml"))
|
||||
(org-babel-tangle)
|
||||
#+end_src
|
||||
|
||||
* Nuclei
|
||||
:PROPERTIES:
|
||||
:header-args: :noweb yes :comments both
|
||||
:END:
|
||||
|
||||
** Type
|
||||
<<<~Nuclei.t~>>>
|
||||
|
||||
#+NAME: types
|
||||
#+begin_src ocaml :tangle (eval mli)
|
||||
open Common
|
||||
|
||||
type t = (Element.t * Coordinate.t) array
|
||||
#+end_src
|
||||
|
||||
#+begin_src ocaml :tangle (eval ml) :exports none
|
||||
open Common
|
||||
|
||||
type t = (Element.t * Coordinate.t) array
|
||||
open Xyz_ast
|
||||
#+end_src
|
||||
|
||||
** xyz file lexer/parser
|
||||
|
||||
*** Lexer
|
||||
|
||||
=nuclei_lexer.mll= contains the description of the lexemes used in
|
||||
an xyz file.
|
||||
|
||||
#+begin_src ocaml :tangle lib/nuclei_lexer.mll :export none
|
||||
{
|
||||
open Xyz_parser
|
||||
}
|
||||
|
||||
let eol = ['\n']
|
||||
let white = [' ' '\t']+
|
||||
let word = [^' ' '\t' '\n']+
|
||||
let letter = ['A'-'Z' 'a'-'z']
|
||||
let integer = ['0'-'9']+
|
||||
let real = '-'? (integer '.' integer | integer '.' | '.' integer) (['e' 'E'] ('+'|'-')? integer)?
|
||||
|
||||
|
||||
rule read_all = parse
|
||||
| eof { EOF }
|
||||
| eol { EOL }
|
||||
| white as w { SPACE w }
|
||||
| integer as i { INTEGER (int_of_string i) }
|
||||
| real as f { FLOAT (float_of_string f) }
|
||||
| word as w { WORD w }
|
||||
|
||||
|
||||
{
|
||||
(* DEBUG
|
||||
let () =
|
||||
let ic = open_in "h2o.xyz" in
|
||||
let lexbuf = Lexing.from_channel ic in
|
||||
while true do
|
||||
let s =
|
||||
match read_all lexbuf with
|
||||
| EOL -> "EOL"
|
||||
| SPACE w -> "SPACE("^w^")"
|
||||
| INTEGER i -> "INTEGER("^(string_of_int i)^")"
|
||||
| FLOAT f -> "FLOAT("^(string_of_float f)^")"
|
||||
| WORD w -> "WORD("^w^")"
|
||||
| EOF -> "EOF"
|
||||
in
|
||||
print_endline s
|
||||
done;
|
||||
*)
|
||||
}
|
||||
#+end_src
|
||||
|
||||
*** Parser
|
||||
|
||||
=xyz_parser.mly= parses nuclear coordinates in xyz format.
|
||||
#+begin_src ocaml :tangle lib/xyz_parser.mly :export none :comments none
|
||||
%{
|
||||
open Common
|
||||
|
||||
let make_angstrom x y z =
|
||||
Coordinate.(make_angstrom {
|
||||
x ; y ; z
|
||||
})
|
||||
|
||||
let output_of f x y z =
|
||||
let a = make_angstrom x y z in
|
||||
fun e ->
|
||||
{
|
||||
Xyz_ast.
|
||||
element = f e;
|
||||
coord = a ;
|
||||
}
|
||||
|
||||
let output_of_string = output_of Element.of_string
|
||||
let output_of_int = output_of Element.of_int
|
||||
|
||||
%}
|
||||
|
||||
%token EOL
|
||||
%token <string> SPACE
|
||||
%token <string> WORD
|
||||
%token <int> INTEGER
|
||||
%token <float> FLOAT
|
||||
%token EOF
|
||||
|
||||
%start input
|
||||
%type <Xyz_ast.xyz_file> input
|
||||
|
||||
%% /* Grammar rules and actions follow */
|
||||
|
||||
input:
|
||||
| integer title atoms_xyz {
|
||||
{
|
||||
number_of_atoms = $1;
|
||||
file_title = $2;
|
||||
nuclei = $3;
|
||||
}
|
||||
}
|
||||
;
|
||||
|
||||
|
||||
integer:
|
||||
| INTEGER EOL { $1 }
|
||||
| INTEGER SPACE EOL { $1 }
|
||||
| SPACE INTEGER EOL { $2 }
|
||||
| SPACE INTEGER SPACE EOL { $2 }
|
||||
;
|
||||
|
||||
title:
|
||||
| title_list EOL { $1 }
|
||||
;
|
||||
|
||||
text:
|
||||
| WORD { $1 }
|
||||
| SPACE { $1 }
|
||||
| FLOAT { (string_of_float $1)}
|
||||
| INTEGER { (string_of_int $1)}
|
||||
;
|
||||
|
||||
title_list:
|
||||
| { "" }
|
||||
| title_list text { ($1 ^ $2) }
|
||||
;
|
||||
|
||||
atoms_xyz:
|
||||
| atoms_list EOL { List.rev $1 }
|
||||
| atoms_list EOF { List.rev $1 }
|
||||
;
|
||||
|
||||
atoms_list:
|
||||
| { [] }
|
||||
| atoms_list WORD SPACE FLOAT SPACE FLOAT SPACE FLOAT EOL { output_of_string $4 $6 $8 $2 :: $1 }
|
||||
| atoms_list WORD SPACE FLOAT SPACE FLOAT SPACE FLOAT SPACE EOL { output_of_string $4 $6 $8 $2 :: $1 }
|
||||
| atoms_list INTEGER SPACE FLOAT SPACE FLOAT SPACE FLOAT EOL { output_of_int $4 $6 $8 $2 :: $1 }
|
||||
| atoms_list INTEGER SPACE FLOAT SPACE FLOAT SPACE FLOAT SPACE EOL { output_of_int $4 $6 $8 $2 :: $1 }
|
||||
| atoms_list SPACE WORD SPACE FLOAT SPACE FLOAT SPACE FLOAT EOL { output_of_string $5 $7 $9 $3 :: $1 }
|
||||
| atoms_list SPACE WORD SPACE FLOAT SPACE FLOAT SPACE FLOAT SPACE EOL { output_of_string $5 $7 $9 $3 :: $1 }
|
||||
| atoms_list SPACE INTEGER SPACE FLOAT SPACE FLOAT SPACE FLOAT EOL { output_of_int $5 $7 $9 $3 :: $1 }
|
||||
| atoms_list SPACE INTEGER SPACE FLOAT SPACE FLOAT SPACE FLOAT SPACE EOL { output_of_int $5 $7 $9 $3 :: $1 }
|
||||
;
|
||||
#+end_src
|
||||
|
||||
When an xyz file is read by =xyz_parser.mly=, it is converted into
|
||||
an ~xyz_file~ data structure.
|
||||
|
||||
#+begin_src ocaml :tangle lib/xyz_ast.mli
|
||||
open Common
|
||||
|
||||
type nucleus =
|
||||
{
|
||||
element: Element.t ;
|
||||
coord : Coordinate.angstrom Coordinate.point;
|
||||
}
|
||||
|
||||
type xyz_file =
|
||||
{
|
||||
number_of_atoms : int ;
|
||||
file_title : string ;
|
||||
nuclei : nucleus list ;
|
||||
}
|
||||
#+end_src
|
||||
|
||||
** Conversion
|
||||
|
||||
#+begin_src ocaml :tangle (eval mli)
|
||||
val of_xyz_string : string -> t
|
||||
val to_xyz_string : t -> string
|
||||
val of_xyz_file : string -> t
|
||||
|
||||
val of_zmt_string : string -> t
|
||||
val of_zmt_file : string -> t
|
||||
|
||||
val to_string : t -> string
|
||||
|
||||
val of_filename : string -> t
|
||||
#+end_src
|
||||
|
||||
| ~of_xyz_string~ | Create from a string, in xyz format |
|
||||
| ~of_xyz_file~ | Create from a file, in xyz format |
|
||||
| ~of_zmt_string~ | Create from a string, in z-matrix format |
|
||||
| ~of_zmt_file~ | Create from a file, in z-matrix format |
|
||||
| ~to_string~ | Transform to a string, for printing |
|
||||
| ~of_filename~ | Detects the type of file (xyz, z-matrix) and reads the file |
|
||||
|
||||
#+begin_src ocaml :tangle (eval ml) :exports none
|
||||
let of_xyz_lexbuf lexbuf =
|
||||
let data =
|
||||
Xyz_parser.input Nuclei_lexer.read_all lexbuf
|
||||
in
|
||||
|
||||
let len = List.length data.nuclei in
|
||||
if len <> data.number_of_atoms then
|
||||
Printf.sprintf "Error: expected %d atoms but %d read"
|
||||
data.number_of_atoms len
|
||||
|> failwith;
|
||||
|
||||
List.map (fun nucleus ->
|
||||
nucleus.element, Coordinate.angstrom_to_bohr nucleus.coord
|
||||
) data.nuclei
|
||||
|> Array.of_list
|
||||
|
||||
|
||||
let of_xyz_string input_string =
|
||||
Lexing.from_string input_string
|
||||
|> of_xyz_lexbuf
|
||||
|
||||
|
||||
let of_xyz_file filename =
|
||||
let ic = open_in filename in
|
||||
let lexbuf =
|
||||
Lexing.from_channel ic
|
||||
in
|
||||
let result =
|
||||
of_xyz_lexbuf lexbuf
|
||||
in
|
||||
close_in ic;
|
||||
result
|
||||
|
||||
|
||||
let of_zmt_string buffer =
|
||||
Zmatrix.of_string buffer
|
||||
|> Zmatrix.to_xyz
|
||||
|> Array.map (fun (e,x,y,z) ->
|
||||
(e, Coordinate.(angstrom_to_bohr @@ make_angstrom { x ; y ; z} ))
|
||||
)
|
||||
|
||||
|
||||
let of_zmt_file filename =
|
||||
let ic = open_in filename in
|
||||
let rec aux accu =
|
||||
try
|
||||
let line = input_line ic in
|
||||
aux (line::accu)
|
||||
with End_of_file ->
|
||||
close_in ic;
|
||||
List.rev accu
|
||||
|> String.concat "\n"
|
||||
in aux []
|
||||
|> of_zmt_string
|
||||
|
||||
|
||||
let to_string atoms =
|
||||
"
|
||||
Nuclear Coordinates (Angstrom)
|
||||
------------------------------
|
||||
|
||||
-----------------------------------------------------------------------
|
||||
Center Atomic Element Coordinates (Angstroms)
|
||||
Number X Y Z
|
||||
-----------------------------------------------------------------------
|
||||
" ^
|
||||
(Array.mapi (fun i (e, coord) ->
|
||||
let open Coordinate in
|
||||
let 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.x coord.y coord.z
|
||||
) atoms
|
||||
|> Array.to_list
|
||||
|> String.concat "\n" ) ^
|
||||
"
|
||||
-----------------------------------------------------------------------
|
||||
|
||||
"
|
||||
|
||||
|
||||
let of_filename filename =
|
||||
of_xyz_file filename
|
||||
|
||||
|
||||
let to_xyz_string t =
|
||||
[ string_of_int (Array.length t) ; "" ] @
|
||||
( Array.map (fun (e, coord) ->
|
||||
let open Coordinate in
|
||||
let coord =
|
||||
bohr_to_angstrom coord
|
||||
in
|
||||
Printf.sprintf " %5s %12.6f %12.6f %12.6f"
|
||||
(Element.to_string e) coord.x coord.y coord.z
|
||||
) t
|
||||
|> Array.to_list )
|
||||
|> String.concat "\n"
|
||||
|
||||
|
||||
#+end_src
|
||||
|
||||
** Query
|
||||
|
||||
#+begin_src ocaml :tangle (eval mli)
|
||||
val formula : t -> string
|
||||
val repulsion : t -> float
|
||||
val charge : t -> Charge.t
|
||||
val small_core : t -> int
|
||||
val large_core : t -> int
|
||||
#+end_src
|
||||
|
||||
| ~formula~ | Returns the chemical formula |
|
||||
| ~repulsion~ | Nuclear repulsion energy, in atomic units |
|
||||
| ~charge~ | Sum of the charges of the nuclei |
|
||||
| ~small_core~ | Number of core electrons in the small core model |
|
||||
| ~large_core~ | Number of core electrons in the large core model |
|
||||
|
||||
#+begin_src ocaml :tangle (eval ml) :exports none
|
||||
let formula t =
|
||||
let dict = Hashtbl.create 67 in
|
||||
Array.iter (fun (e,_) ->
|
||||
let e = Element.to_string e in
|
||||
let value =
|
||||
try (Hashtbl.find dict e) + 1
|
||||
with Not_found -> 1
|
||||
in
|
||||
Hashtbl.replace dict e value
|
||||
) t;
|
||||
Hashtbl.to_seq_keys dict
|
||||
|> List.of_seq
|
||||
|> List.sort String.compare
|
||||
|> List.fold_left (fun accu key ->
|
||||
let x = Hashtbl.find dict key in
|
||||
accu ^ key ^ "_{" ^ (string_of_int x) ^ "}") ""
|
||||
|
||||
|
||||
|
||||
let repulsion nuclei =
|
||||
let get_charge e =
|
||||
Element.to_charge e
|
||||
|> Charge.to_float
|
||||
in
|
||||
Array.fold_left ( fun accu (e1, coord1) ->
|
||||
accu +.
|
||||
Array.fold_left (fun accu (e2, coord2) ->
|
||||
let r = Coordinate.(norm (coord1 |- coord2)) in
|
||||
if r > 0. then
|
||||
accu +. 0.5 *. (get_charge e2) *. (get_charge e1) /. r
|
||||
else accu
|
||||
) 0. nuclei
|
||||
) 0. nuclei
|
||||
|
||||
|
||||
let charge nuclei =
|
||||
Array.fold_left (fun accu (e, _) -> accu + Charge.to_int (Element.to_charge e) )
|
||||
0 nuclei
|
||||
|> Charge.of_int
|
||||
|
||||
|
||||
let small_core a =
|
||||
Array.fold_left (fun accu (e,_) -> accu + (Element.small_core e)) 0 a
|
||||
|
||||
|
||||
let large_core a =
|
||||
Array.fold_left (fun accu (e,_) -> accu + (Element.large_core e)) 0 a
|
||||
|
||||
#+end_src
|
||||
|
||||
** TREXIO
|
||||
|
||||
*** Read
|
||||
|
||||
#+begin_src ocaml :tangle (eval mli)
|
||||
val of_trexio : Trexio.trexio_file -> t
|
||||
#+end_src
|
||||
|
||||
#+begin_src ocaml :tangle (eval ml) :exports none
|
||||
let of_trexio f =
|
||||
let num = Trexio.read_nucleus_num f in
|
||||
let charge = Trexio.read_nucleus_charge f
|
||||
|> Array.map Charge.of_float in
|
||||
let coord = Trexio.read_nucleus_coord f in
|
||||
Array.init num (fun i ->
|
||||
let coord = Coordinate.{ x = coord.(3*i) ;
|
||||
y = coord.(3*i+1) ;
|
||||
z = coord.(3*i+2) } in
|
||||
(Element.of_charge charge.(i), Coordinate.make coord)
|
||||
)
|
||||
#+end_src
|
||||
|
||||
*** Write
|
||||
|
||||
#+begin_src ocaml :tangle (eval mli)
|
||||
val to_trexio : Trexio.trexio_file -> t -> unit
|
||||
#+end_src
|
||||
|
||||
#+begin_src ocaml :tangle (eval ml) :exports none
|
||||
let to_trexio f t =
|
||||
let num = Array.length t in
|
||||
Trexio.write_nucleus_num f num;
|
||||
|
||||
Array.map (fun (e, _) -> Element.to_charge e |> Charge.to_float) t
|
||||
|> Trexio.write_nucleus_charge f;
|
||||
|
||||
Array.map (fun (e, _) -> Element.to_string e) t
|
||||
|> Trexio.write_nucleus_label f;
|
||||
|
||||
let coord = Array.init (num*3) (fun _ -> 0.) in
|
||||
Array.iteri (fun i (_, xyz) ->
|
||||
coord.(3*i) <- Coordinate.(get X xyz) ;
|
||||
coord.(3*i+1) <- Coordinate.(get Y xyz) ;
|
||||
coord.(3*i+2) <- Coordinate.(get Z xyz) ) t;
|
||||
Trexio.write_nucleus_coord f coord;
|
||||
|
||||
repulsion t
|
||||
|> Trexio.write_nucleus_repulsion f
|
||||
|
||||
#+end_src
|
||||
|
||||
** Printers
|
||||
|
||||
#+begin_src ocaml :tangle (eval mli)
|
||||
val pp : Format.formatter -> t -> unit
|
||||
#+end_src
|
||||
|
||||
#+begin_src ocaml :tangle (eval ml) :exports none
|
||||
let pp ppf t =
|
||||
Format.fprintf ppf "@[%s@]" (to_string t)
|
||||
#+end_src
|
||||
|
||||
** Tests
|
||||
|
||||
#+begin_src ocaml :tangle (eval test-ml) :exports none
|
||||
open Common
|
||||
open Particles
|
||||
open Alcotest
|
||||
|
||||
let wd = Common.Qcaml.root ^ Filename.dir_sep ^ "test"
|
||||
|
||||
let test_xyz molecule length repulsion charge core =
|
||||
let xyz = Nuclei.of_xyz_file (wd^Filename.dir_sep^molecule^".xyz") in
|
||||
check int "length" length (Array.length xyz);
|
||||
check (float 1.e-4) "repulsion" repulsion (Nuclei.repulsion xyz);
|
||||
check int "charge" charge (Charge.to_int @@ Nuclei.charge xyz);
|
||||
check int "small_core" core (Nuclei.small_core xyz);
|
||||
()
|
||||
|
||||
let tests = [
|
||||
"caffeine", `Quick, (fun () -> test_xyz "caffeine" 24 917.0684 102 28);
|
||||
"water", `Quick, (fun () -> test_xyz "water" 3 9.19497 10 2);
|
||||
]
|
||||
#+end_src
|
||||
|
||||
|
|
@ -1,397 +0,0 @@
|
|||
#+begin_src elisp tangle: no :results none :exports none
|
||||
(setq pwd (file-name-directory buffer-file-name))
|
||||
(setq name (file-name-nondirectory (substring buffer-file-name 0 -4)))
|
||||
(setq lib (concat pwd "lib/"))
|
||||
(setq testdir (concat pwd "test/"))
|
||||
(setq mli (concat lib name ".mli"))
|
||||
(setq ml (concat lib name ".ml"))
|
||||
(setq test-ml (concat testdir name ".ml"))
|
||||
(org-babel-tangle)
|
||||
#+end_src
|
||||
|
||||
* Z-matrix
|
||||
:PROPERTIES:
|
||||
:header-args: :noweb yes :comments both
|
||||
:END:
|
||||
|
||||
Z-matrix representation of nuclear coordinates.
|
||||
|
||||
** Type
|
||||
|
||||
#+NAME: types
|
||||
#+begin_src ocaml :tangle (eval mli)
|
||||
type t
|
||||
#+end_src
|
||||
|
||||
#+begin_src ocaml :tangle (eval ml) :exports none
|
||||
module StringMap = Map.Make(String)
|
||||
|
||||
type atom_id = int
|
||||
type angle = Label of string | Value of float
|
||||
type distance = Label of string | Value of float
|
||||
type dihedral = Label of string | Value of float
|
||||
|
||||
type line =
|
||||
| First of Element.t
|
||||
| Second of (Element.t * distance)
|
||||
| Third of (Element.t * atom_id * distance * atom_id * angle)
|
||||
| Other of (Element.t * atom_id * distance * atom_id * angle * atom_id * dihedral )
|
||||
| Coord of (string * float)
|
||||
|
||||
type t = (line array * float StringMap.t)
|
||||
#+end_src
|
||||
|
||||
** Conversion
|
||||
|
||||
#+begin_src ocaml :tangle (eval mli)
|
||||
val of_string : string -> t
|
||||
val to_xyz : t -> (Element.t * float * float * float) array
|
||||
val to_xyz_string : t -> string
|
||||
#+end_src
|
||||
|
||||
| ~of_string~ | Reads a z-matrix from a string |
|
||||
| ~to_xyz~ | Converts to xyz format, as in the ~Nuclei~ module |
|
||||
| ~to_xyz_string~ | Converts to xyz format, as a string |
|
||||
|
||||
#+begin_example
|
||||
let zmt = Zmatrix.of_string "
|
||||
n
|
||||
n 1 nn
|
||||
h 1 hn 2 hnn
|
||||
h 2 hn 1 hnn 3 dih4
|
||||
h 1 hn 2 hnn 4 dih5
|
||||
h 2 hn 1 hnn 3 dih5
|
||||
|
||||
nn 1.446
|
||||
hn 1.016
|
||||
hnn 106.0
|
||||
dih4 -54.38
|
||||
dih5 54.38
|
||||
" ;;
|
||||
- : Zmatrix.t = N
|
||||
N 1 1.446000
|
||||
H 1 1.016000 2 106.000000
|
||||
H 2 1.016000 1 106.000000 3 -54.380000
|
||||
H 1 1.016000 2 106.000000 4 54.380000
|
||||
H 2 1.016000 1 106.000000 3 54.380000
|
||||
|
||||
|
||||
Zmatrix.to_xyz zmt ;;
|
||||
- : (Element.t * float * float * float) array =
|
||||
[|(N, 0., 0., 0.); (N, 0., 0., 1.446);
|
||||
(H, -0.976641883073332107, 0., -0.280047553510071046);
|
||||
(H, -0.568802835186988709, 0.793909757123734683, 1.726047553510071);
|
||||
(H, 0.314092649983635563, 0.924756819385119, -0.280047553510071101);
|
||||
(H, -0.568802835186988709, -0.793909757123734683, 1.726047553510071)|]
|
||||
|
||||
|
||||
Zmatrix.to_xyz_string zmt ;;
|
||||
- : string =
|
||||
"N 0.000000 0.000000 0.000000
|
||||
N 0.000000 0.000000 1.446000
|
||||
H -0.976642 0.000000 -0.280048
|
||||
H -0.568803 0.793910 1.726048
|
||||
H 0.314093 0.924757 -0.280048
|
||||
H -0.568803 -0.793910 1.726048"
|
||||
#+end_example
|
||||
|
||||
|
||||
#+begin_src ocaml :tangle (eval ml) :exports none
|
||||
let pi = Common.Constants.pi
|
||||
let to_radian = pi /. 180.
|
||||
|
||||
let rec in_range (xmin, xmax) x =
|
||||
if (x <= xmin) then
|
||||
in_range (xmin, xmax) (x -. xmin +. xmax )
|
||||
else if (x > xmax) then
|
||||
in_range (xmin, xmax) (x -. xmax +. xmin )
|
||||
else
|
||||
x
|
||||
|
||||
let atom_id_of_int : int -> atom_id =
|
||||
fun x -> ( assert (x>0) ; x)
|
||||
|
||||
let distance_of_float : float -> distance =
|
||||
fun x -> ( assert (x>=0.) ; Value x)
|
||||
|
||||
let angle_of_float : float -> angle =
|
||||
fun x -> Value (in_range (-180., 180.) x)
|
||||
|
||||
let dihedral_of_float : float -> dihedral =
|
||||
fun x -> Value (in_range (-360., 360.) x)
|
||||
|
||||
|
||||
let atom_id_of_string : string -> atom_id =
|
||||
fun i -> atom_id_of_int @@ int_of_string i
|
||||
|
||||
let distance_of_string : string -> distance =
|
||||
fun s ->
|
||||
try
|
||||
distance_of_float @@ float_of_string s
|
||||
with _ -> Label s
|
||||
|
||||
let angle_of_string : string -> angle =
|
||||
fun s ->
|
||||
try
|
||||
angle_of_float @@ float_of_string s
|
||||
with _ -> Label s
|
||||
|
||||
let dihedral_of_string : string -> dihedral =
|
||||
fun s ->
|
||||
try
|
||||
dihedral_of_float @@ float_of_string s
|
||||
with _ -> Label s
|
||||
|
||||
let int_of_atom_id : atom_id -> int = fun x -> x
|
||||
|
||||
let float_of_distance : float StringMap.t -> distance -> float =
|
||||
fun map -> function
|
||||
| Value x -> x
|
||||
| Label s -> StringMap.find s map
|
||||
|
||||
let float_of_angle : float StringMap.t -> angle -> float =
|
||||
fun map -> function
|
||||
| Value x -> x
|
||||
| Label s -> StringMap.find s map
|
||||
|
||||
let float_of_dihedral : float StringMap.t -> dihedral -> float =
|
||||
fun map -> function
|
||||
| Value x -> x
|
||||
| Label s -> StringMap.find s map
|
||||
|
||||
|
||||
let string_of_line map =
|
||||
let f_r = float_of_distance map
|
||||
and f_a = float_of_angle map
|
||||
and f_d = float_of_dihedral map
|
||||
and i_i = int_of_atom_id
|
||||
in function
|
||||
| First e -> Printf.sprintf "%-3s" (Element.to_string e)
|
||||
| Second (e, r) -> Printf.sprintf "%-3s %5d %f" (Element.to_string e) 1 (f_r r)
|
||||
| Third (e, i, r, j, a) -> Printf.sprintf "%-3s %5d %f %5d %f" (Element.to_string e) (i_i i) (f_r r) (i_i j) (f_a a)
|
||||
| Other (e, i, r, j, a, k, d) -> Printf.sprintf "%-3s %5d %f %5d %f %5d %f" (Element.to_string e) (i_i i) (f_r r) (i_i j) (f_a a) (i_i k) (f_d d)
|
||||
| Coord (c, f) -> Printf.sprintf "%s %f" c f
|
||||
|
||||
|
||||
let line_of_string l =
|
||||
let line_clean =
|
||||
Str.split (Str.regexp " ") l
|
||||
|> List.filter (fun x -> x <> "")
|
||||
in
|
||||
match line_clean with
|
||||
| e :: [] -> First (Element.of_string e)
|
||||
| e :: _ :: r :: [] -> Second
|
||||
(Element.of_string e,
|
||||
distance_of_string r)
|
||||
| e :: i :: r :: j :: a :: [] -> Third
|
||||
(Element.of_string e,
|
||||
atom_id_of_string i,
|
||||
distance_of_string r,
|
||||
atom_id_of_string j,
|
||||
angle_of_string a)
|
||||
| e :: i :: r :: j :: a :: k :: d :: [] -> Other
|
||||
(Element.of_string e,
|
||||
atom_id_of_string i,
|
||||
distance_of_string r,
|
||||
atom_id_of_string j,
|
||||
angle_of_string a,
|
||||
atom_id_of_string k,
|
||||
dihedral_of_string d)
|
||||
| c :: f :: [] -> Coord (c, float_of_string f)
|
||||
| _ -> failwith ("Syntax error: "^l)
|
||||
|
||||
|
||||
let of_string t =
|
||||
let l =
|
||||
Str.split (Str.regexp "\n") t
|
||||
|> List.map String.trim
|
||||
|> List.filter (fun x -> x <> "")
|
||||
|> List.map line_of_string
|
||||
in
|
||||
|
||||
let l =
|
||||
match l with
|
||||
| First _ :: Second _ :: Third _ :: _
|
||||
| First _ :: Second _ :: Coord _ :: []
|
||||
| First _ :: Second _ :: []
|
||||
| First _ :: [] -> l
|
||||
| _ -> failwith "Syntax error"
|
||||
in
|
||||
|
||||
let (l, m) =
|
||||
let rec work lst map = function
|
||||
| (First _ as x) :: rest
|
||||
| (Second _ as x) :: rest
|
||||
| (Third _ as x) :: rest
|
||||
| (Other _ as x) :: rest -> work (x::lst) map rest
|
||||
| (Coord (c,f)) :: rest -> work lst (StringMap.add c f map) rest
|
||||
| [] -> (List.rev lst, map)
|
||||
in
|
||||
work [] (StringMap.empty) l
|
||||
in
|
||||
(Array.of_list l, m)
|
||||
|
||||
|
||||
(** Linear algebra *)
|
||||
|
||||
let (|-) (x,y,z) (x',y',z') =
|
||||
( x-.x', y-.y', z-.z' )
|
||||
|
||||
let (|+) (x,y,z) (x',y',z') =
|
||||
( x+.x', y+.y', z+.z' )
|
||||
|
||||
let (|.) s (x,y,z) =
|
||||
( s*.x, s*.y, s*.z )
|
||||
|
||||
let dot (x,y,z) (x',y',z') =
|
||||
x*.x' +. y*.y' +. z*.z'
|
||||
|
||||
let norm u =
|
||||
sqrt @@ dot u u
|
||||
|
||||
let normalized u =
|
||||
1. /. (norm u) |. u
|
||||
|
||||
let cross (x,y,z) (x',y',z') =
|
||||
((y *. z' -. z *. y'), -. (x *. z' -. z *. x'), (x *. y' -. y *. x'))
|
||||
|
||||
let rotation_matrix axis angle =
|
||||
(* Euler-Rodrigues formula for rotation matrix, taken from
|
||||
https://github.com/jevandezande/zmatrix/blob/master/converter.py
|
||||
,*)
|
||||
let a =
|
||||
(cos (angle *. to_radian *. 0.5))
|
||||
in
|
||||
let (b, c, d) =
|
||||
(-. sin (angle *. to_radian *. 0.5)) |. (normalized axis)
|
||||
in
|
||||
Array.of_list @@
|
||||
[(a *. a +. b *. b -. c *. c -. d *. d,
|
||||
2. *. (b *. c -. a *. d),
|
||||
2. *. (b *. d +. a *. c));
|
||||
(2. *. (b *. c +. a *. d),
|
||||
a *. a +. c *. c -.b *. b -. d *. d,
|
||||
2. *. (c *. d -. a *. b));
|
||||
(2. *. (b *. d -. a *. c),
|
||||
2. *. (c *. d +. a *. b),
|
||||
a *. a +. d *. d -. b *. b -. c *. c)]
|
||||
|
||||
|
||||
let apply_rotation_matrix rot u =
|
||||
(dot rot.(0) u, dot rot.(1) u, dot rot.(2) u)
|
||||
|
||||
|
||||
let to_xyz (z,map) =
|
||||
let result =
|
||||
Array.make (Array.length z) None
|
||||
in
|
||||
|
||||
let get_cartesian_coord i =
|
||||
match result.(i-1) with
|
||||
| None -> failwith @@ Printf.sprintf "Atom %d is defined in the future" i
|
||||
| Some (_, x, y, z) -> (x, y, z)
|
||||
in
|
||||
|
||||
|
||||
let append_line i' =
|
||||
match z.(i') with
|
||||
| First e ->
|
||||
result.(i') <- Some (e, 0., 0., 0.)
|
||||
| Second (e, r) ->
|
||||
let r =
|
||||
float_of_distance map r
|
||||
in
|
||||
result.(i') <- Some (e, 0., 0., r)
|
||||
| Third (e, i, r, j, a) ->
|
||||
begin
|
||||
let i, r, j, a =
|
||||
int_of_atom_id i,
|
||||
float_of_distance map r,
|
||||
int_of_atom_id j,
|
||||
float_of_angle map a
|
||||
in
|
||||
let ui, uj =
|
||||
get_cartesian_coord i,
|
||||
get_cartesian_coord j
|
||||
in
|
||||
let u_ij =
|
||||
(uj |- ui)
|
||||
in
|
||||
let rot =
|
||||
rotation_matrix (0., 1., 0.) a
|
||||
in
|
||||
let new_vec =
|
||||
apply_rotation_matrix rot ( r |. (normalized u_ij))
|
||||
in
|
||||
let (x, y, z) =
|
||||
new_vec |+ ui
|
||||
in
|
||||
result.(i') <- Some (e, x, y, z)
|
||||
end
|
||||
| Other (e, i, r, j, a, k, d) ->
|
||||
begin
|
||||
let i, r, j, a, k, d =
|
||||
int_of_atom_id i,
|
||||
float_of_distance map r,
|
||||
int_of_atom_id j,
|
||||
float_of_angle map a,
|
||||
int_of_atom_id k,
|
||||
float_of_dihedral map d
|
||||
in
|
||||
let ui, uj, uk =
|
||||
get_cartesian_coord i,
|
||||
get_cartesian_coord j,
|
||||
get_cartesian_coord k
|
||||
in
|
||||
let u_ij, u_kj =
|
||||
(uj |- ui) , (uj |- uk)
|
||||
in
|
||||
let normal =
|
||||
cross u_ij u_kj
|
||||
in
|
||||
let new_vec =
|
||||
r |. (normalized u_ij)
|
||||
|> apply_rotation_matrix (rotation_matrix normal a)
|
||||
|> apply_rotation_matrix (rotation_matrix u_ij d)
|
||||
in
|
||||
let (x, y, z) =
|
||||
new_vec |+ ui
|
||||
in
|
||||
result.(i') <- Some (e, x, y, z)
|
||||
end
|
||||
| Coord _ -> ()
|
||||
in
|
||||
Array.iteri (fun i _ -> append_line i) z;
|
||||
let result =
|
||||
Array.map (function
|
||||
| Some x -> x
|
||||
| None -> failwith "Some atoms were not defined" ) result
|
||||
in
|
||||
result
|
||||
|
||||
|
||||
let to_xyz_string (l,map) =
|
||||
String.concat "\n"
|
||||
( to_xyz (l,map)
|
||||
|> Array.map (fun (e,x,y,z) ->
|
||||
Printf.sprintf "%s %f %f %f" (Element.to_string e) x y z)
|
||||
|> Array.to_list
|
||||
)
|
||||
#+end_src
|
||||
|
||||
** Printers
|
||||
|
||||
#+begin_src ocaml :tangle (eval mli)
|
||||
val pp : Format.formatter -> t -> unit
|
||||
#+end_src
|
||||
|
||||
#+begin_src ocaml :tangle (eval ml) :exports none
|
||||
let pp ppf (a, map) =
|
||||
let f = string_of_line map in
|
||||
Format.fprintf ppf "@[";
|
||||
Array.iter (fun line ->
|
||||
Format.fprintf ppf "%s@." (f line)
|
||||
) a;
|
||||
Format.fprintf ppf "@]"
|
||||
#+end_src
|
||||
|
||||
|
|
@ -1,24 +1,13 @@
|
|||
(* [[file:~/QCaml/perturbation/mp2.org::*Type][Type:2]] *)
|
||||
(** Type *)
|
||||
|
||||
type t = {
|
||||
energy : float ;
|
||||
mo_basis : Mo.Basis.t ;
|
||||
frozen_core : Mo.Frozen_core.t ;
|
||||
}
|
||||
(* Type:2 ends here *)
|
||||
|
||||
(** Creation *)
|
||||
|
||||
|
||||
(* | ~make~ | Creates an MP2 data structure |
|
||||
*
|
||||
* #+begin_example
|
||||
* let mp2 =
|
||||
* Perturbation.Mp2.make ~frozen_core:(Mo.Frozen_core.(make Small nuclei)) mo_basis
|
||||
* ;;
|
||||
* val mp2 : Perturbation.Mp2.t = E(MP2)=-0.185523
|
||||
* #+end_example *)
|
||||
|
||||
|
||||
(* [[file:~/QCaml/perturbation/mp2.org::*Creation][Creation:2]] *)
|
||||
open Linear_algebra
|
||||
|
||||
let make_rmp2 mo_basis mo_class =
|
||||
|
|
@ -37,22 +26,22 @@ let make_rmp2 mo_basis mo_class =
|
|||
List.fold_left (fun accu b ->
|
||||
match b with Mo.Class.Virtual b ->
|
||||
let eps = -. (epsilon%.(b)) in
|
||||
accu +.
|
||||
accu +.
|
||||
List.fold_left (fun accu a ->
|
||||
match a with Mo.Class.Virtual a ->
|
||||
let eps = eps -. (epsilon%.(a)) in
|
||||
accu +.
|
||||
accu +.
|
||||
List.fold_left (fun accu j ->
|
||||
match j with Mo.Class.Inactive j ->
|
||||
let eps = eps +. epsilon%.(j) in
|
||||
accu +.
|
||||
accu +.
|
||||
List.fold_left (fun accu i ->
|
||||
match i with Mo.Class.Inactive i ->
|
||||
let eps = eps +. epsilon%.(i) in
|
||||
let ijab = Four_idx_storage.get_phys eri i j a b
|
||||
and abji = Four_idx_storage.get_phys eri a b j i in
|
||||
let abij = ijab in
|
||||
accu +. ijab *. ( abij +. abij -. abji) /. eps
|
||||
accu +. ijab *. ( abij +. abij -. abji) /. eps
|
||||
| _ -> accu
|
||||
) 0. inactives
|
||||
| _ -> accu
|
||||
|
|
@ -70,7 +59,7 @@ let make ~frozen_core mo_basis =
|
|||
|> Mo.Class.to_list
|
||||
in
|
||||
|
||||
let energy =
|
||||
let energy =
|
||||
match Mo.Basis.mo_type mo_basis with
|
||||
| RHF -> make_rmp2 mo_basis mo_class
|
||||
| ROHF -> Common.Util.not_implemented "ROHF MP2"
|
||||
|
|
@ -78,26 +67,17 @@ let make ~frozen_core mo_basis =
|
|||
| _ -> invalid_arg "MP2 needs RHF or ROHF MOs"
|
||||
in
|
||||
{ energy ; mo_basis ; frozen_core }
|
||||
(* Creation:2 ends here *)
|
||||
|
||||
|
||||
(** Access *)
|
||||
|
||||
(* | ~energy~ | Returns the MP2 energy |
|
||||
* | ~mo_basis~ | Returns the MO basis on which the MP2 energy was computed |
|
||||
* | ~frozen_core~ | Returns the frozen_core scheme used to compute the MP2 energy |
|
||||
*
|
||||
* #+begin_example
|
||||
*
|
||||
* #+end_example *)
|
||||
|
||||
|
||||
(* [[file:~/QCaml/perturbation/mp2.org::*Access][Access:2]] *)
|
||||
let energy t = t.energy
|
||||
let mo_basis t = t.mo_basis
|
||||
let frozen_core t = t.frozen_core
|
||||
(* Access:2 ends here *)
|
||||
|
||||
(* [[file:~/QCaml/perturbation/mp2.org::*Printers][Printers:2]] *)
|
||||
|
||||
(** Printers *)
|
||||
|
||||
let pp ppf t =
|
||||
Format.fprintf ppf "@[E(MP2)=%f@]" t.energy
|
||||
(* Printers:2 ends here *)
|
||||
|
||||
|
|
|
|||
|
|
@ -1,29 +1,36 @@
|
|||
(* Type *)
|
||||
(** Data structure for an MP2 calculation *)
|
||||
|
||||
(**
|
||||
* let mp2 =
|
||||
* Perturbation.Mp2.make ~frozen_core:(Mo.Frozen_core.(make Small nuclei)) mo_basis
|
||||
* ;;
|
||||
* val mp2 : Perturbation.Mp2.t = E(MP2)=-0.185523
|
||||
*
|
||||
*)
|
||||
|
||||
|
||||
(* [[file:~/QCaml/perturbation/mp2.org::*Type][Type:1]] *)
|
||||
(** Type *)
|
||||
|
||||
type t
|
||||
(* Type:1 ends here *)
|
||||
|
||||
(* Creation *)
|
||||
|
||||
|
||||
(* [[file:~/QCaml/perturbation/mp2.org::*Creation][Creation:1]] *)
|
||||
(** Creation *)
|
||||
|
||||
val make : frozen_core:Mo.Frozen_core.t -> Mo.Basis.t -> t
|
||||
(* Creation:1 ends here *)
|
||||
|
||||
(* Access *)
|
||||
|
||||
|
||||
(* [[file:~/QCaml/perturbation/mp2.org::*Access][Access:1]] *)
|
||||
(** Access *)
|
||||
|
||||
val energy : t -> float
|
||||
(* Returns the MP2 energy *)
|
||||
|
||||
val mo_basis : t -> Mo.Basis.t
|
||||
(* Returns the MO basis on which the MP2 energy was computed *)
|
||||
|
||||
val frozen_core : t -> Mo.Frozen_core.t
|
||||
(* Access:1 ends here *)
|
||||
|
||||
(* Printers *)
|
||||
(* Returns the frozen_core scheme used to compute the MP2 energy *)
|
||||
|
||||
|
||||
(* [[file:~/QCaml/perturbation/mp2.org::*Printers][Printers:1]] *)
|
||||
(** Printers *)
|
||||
|
||||
val pp : Format.formatter -> t -> unit
|
||||
(* Printers:1 ends here *)
|
||||
|
|
|
|||
|
|
@ -1,180 +0,0 @@
|
|||
#+begin_src elisp tangle: no :results none :exports none
|
||||
(setq pwd (file-name-directory buffer-file-name))
|
||||
(setq name (file-name-nondirectory (substring buffer-file-name 0 -4)))
|
||||
(setq lib (concat pwd "lib/"))
|
||||
(setq testdir (concat pwd "test/"))
|
||||
(setq mli (concat lib name ".mli"))
|
||||
(setq ml (concat lib name ".ml"))
|
||||
(setq test-ml (concat testdir name ".ml"))
|
||||
(org-babel-tangle)
|
||||
#+end_src
|
||||
|
||||
* MP2
|
||||
:PROPERTIES:
|
||||
:header-args: :noweb yes :comments both
|
||||
:END:
|
||||
|
||||
** Type
|
||||
|
||||
#+begin_src ocaml :tangle (eval mli)
|
||||
type t
|
||||
#+end_src
|
||||
|
||||
#+begin_src ocaml :tangle (eval ml) :exports none
|
||||
type t = {
|
||||
energy : float ;
|
||||
mo_basis : Mo.Basis.t ;
|
||||
frozen_core : Mo.Frozen_core.t ;
|
||||
}
|
||||
#+end_src
|
||||
|
||||
** Creation
|
||||
|
||||
#+begin_src ocaml :tangle (eval mli)
|
||||
val make : frozen_core:Mo.Frozen_core.t -> Mo.Basis.t -> t
|
||||
#+end_src
|
||||
|
||||
| ~make~ | Creates an MP2 data structure |
|
||||
|
||||
#+begin_example
|
||||
let mp2 =
|
||||
Perturbation.Mp2.make ~frozen_core:(Mo.Frozen_core.(make Small nuclei)) mo_basis
|
||||
;;
|
||||
val mp2 : Perturbation.Mp2.t = E(MP2)=-0.185523
|
||||
#+end_example
|
||||
|
||||
#+begin_src ocaml :tangle (eval ml) :exports none
|
||||
open Linear_algebra
|
||||
|
||||
let make_rmp2 mo_basis mo_class =
|
||||
|
||||
let epsilon = Mo.Basis.mo_energies mo_basis in
|
||||
let eri = Mo.Basis.ee_ints mo_basis in
|
||||
|
||||
let inactives =
|
||||
List.filter (fun i ->
|
||||
match i with Mo.Class.Inactive _ -> true | _ -> false) mo_class
|
||||
and virtuals =
|
||||
List.filter (fun i ->
|
||||
match i with Mo.Class.Virtual _ -> true | _ -> false) mo_class
|
||||
in
|
||||
|
||||
List.fold_left (fun accu b ->
|
||||
match b with Mo.Class.Virtual b ->
|
||||
let eps = -. (epsilon%.(b)) in
|
||||
accu +.
|
||||
List.fold_left (fun accu a ->
|
||||
match a with Mo.Class.Virtual a ->
|
||||
let eps = eps -. (epsilon%.(a)) in
|
||||
accu +.
|
||||
List.fold_left (fun accu j ->
|
||||
match j with Mo.Class.Inactive j ->
|
||||
let eps = eps +. epsilon%.(j) in
|
||||
accu +.
|
||||
List.fold_left (fun accu i ->
|
||||
match i with Mo.Class.Inactive i ->
|
||||
let eps = eps +. epsilon%.(i) in
|
||||
let ijab = Four_idx_storage.get_phys eri i j a b
|
||||
and abji = Four_idx_storage.get_phys eri a b j i in
|
||||
let abij = ijab in
|
||||
accu +. ijab *. ( abij +. abij -. abji) /. eps
|
||||
| _ -> accu
|
||||
) 0. inactives
|
||||
| _ -> accu
|
||||
) 0. inactives
|
||||
| _ -> accu
|
||||
) 0. virtuals
|
||||
| _ -> accu
|
||||
) 0. virtuals
|
||||
|
||||
|
||||
let make ~frozen_core mo_basis =
|
||||
|
||||
let mo_class =
|
||||
Mo.Class.cas_sd mo_basis ~frozen_core 0 0
|
||||
|> Mo.Class.to_list
|
||||
in
|
||||
|
||||
let energy =
|
||||
match Mo.Basis.mo_type mo_basis with
|
||||
| RHF -> make_rmp2 mo_basis mo_class
|
||||
| ROHF -> Common.Util.not_implemented "ROHF MP2"
|
||||
| UHF -> Common.Util.not_implemented "UHF MP2"
|
||||
| _ -> invalid_arg "MP2 needs RHF or ROHF MOs"
|
||||
in
|
||||
{ energy ; mo_basis ; frozen_core }
|
||||
#+end_src
|
||||
|
||||
** Access
|
||||
|
||||
#+begin_src ocaml :tangle (eval mli)
|
||||
val energy : t -> float
|
||||
val mo_basis : t -> Mo.Basis.t
|
||||
val frozen_core : t -> Mo.Frozen_core.t
|
||||
#+end_src
|
||||
|
||||
| ~energy~ | Returns the MP2 energy |
|
||||
| ~mo_basis~ | Returns the MO basis on which the MP2 energy was computed |
|
||||
| ~frozen_core~ | Returns the frozen_core scheme used to compute the MP2 energy |
|
||||
|
||||
#+begin_example
|
||||
|
||||
#+end_example
|
||||
|
||||
#+begin_src ocaml :tangle (eval ml) :exports none
|
||||
let energy t = t.energy
|
||||
let mo_basis t = t.mo_basis
|
||||
let frozen_core t = t.frozen_core
|
||||
#+end_src
|
||||
|
||||
** Printers
|
||||
|
||||
#+begin_src ocaml :tangle (eval mli)
|
||||
val pp : Format.formatter -> t -> unit
|
||||
#+end_src
|
||||
|
||||
#+begin_src ocaml :tangle (eval ml) :exports none
|
||||
let pp ppf t =
|
||||
Format.fprintf ppf "@[E(MP2)=%f@]" t.energy
|
||||
#+end_src
|
||||
|
||||
** Tests
|
||||
|
||||
#+begin_src ocaml :tangle (eval test-ml) :exports none
|
||||
open Alcotest
|
||||
open Particles
|
||||
|
||||
let wd = Common.Qcaml.root ^ Filename.dir_sep ^ "test"
|
||||
|
||||
let tests =
|
||||
[ "HF Water", `Quick, fun () ->
|
||||
let nuclei =
|
||||
wd ^ Filename.dir_sep ^ "water.xyz"
|
||||
|> Nuclei.of_xyz_file
|
||||
in
|
||||
let basis_filename =
|
||||
wd ^ Filename.dir_sep ^ "cc-pvdz"
|
||||
in
|
||||
let ao_basis =
|
||||
Ao.Basis.of_nuclei_and_basis_filename ~kind:`Gaussian
|
||||
~cartesian:false ~nuclei basis_filename
|
||||
in
|
||||
|
||||
let simulation = Simulation.make ~nuclei ao_basis in
|
||||
|
||||
let hf = Mo.Hartree_fock.make ~guess:`Huckel simulation in
|
||||
Format.printf "%a" (Mo.Hartree_fock.pp) hf;
|
||||
check (float 2.e-10) "Energy" (-76.0267987005) (Mo.Hartree_fock.energy hf);
|
||||
|
||||
let mo_basis = Mo.Basis.of_hartree_fock hf in
|
||||
|
||||
let frozen_core = Mo.Frozen_core.(make Small nuclei) in
|
||||
|
||||
let mp2 = Perturbation.Mp2.make ~frozen_core mo_basis in
|
||||
|
||||
let e_mp2 = Perturbation.Mp2.energy mp2 in
|
||||
|
||||
check (float 1.e-9) "MP2" (-0.2016211415) (e_mp2)
|
||||
]
|
||||
#+end_src
|
||||
|
||||
|
|
@ -1,10 +1,12 @@
|
|||
(* [[file:../simulation.org::*Simulation][Simulation:2]] *)
|
||||
(** Simulation *)
|
||||
|
||||
open Common
|
||||
open Particles
|
||||
open Operators
|
||||
(* Simulation:2 ends here *)
|
||||
|
||||
(* [[file:../simulation.org::*Type][Type:2]] *)
|
||||
|
||||
(** Type *)
|
||||
|
||||
type t = {
|
||||
charge : Charge.t;
|
||||
electrons : Electrons.t;
|
||||
|
|
@ -12,36 +14,20 @@ type t = {
|
|||
ao_basis : Ao.Basis.t;
|
||||
operators : Operator.t list;
|
||||
}
|
||||
(* Type:2 ends here *)
|
||||
|
||||
|
||||
(** Access *)
|
||||
|
||||
(* | ~nuclei~ | Nuclear coordinates used in the smiulation |
|
||||
* | ~charge~ | Total charge (electrons + nuclei) |
|
||||
* | ~electrons~ | Electrons used in the simulation |
|
||||
* | ~ao_basis~ | Atomic basis set |
|
||||
* | ~nuclear_repulsion~ | Nuclear repulsion energy |
|
||||
* | ~operators~ | List of extra operators (range-separation, f12, etc) | *)
|
||||
|
||||
|
||||
(* [[file:../simulation.org::*Access][Access:2]] *)
|
||||
let nuclei t = t.nuclei
|
||||
let charge t = t.charge
|
||||
let electrons t = t.electrons
|
||||
let ao_basis t = t.ao_basis
|
||||
let nuclear_repulsion t = Nuclei.repulsion @@ nuclei t
|
||||
let operators t = t.operators
|
||||
(* Access:2 ends here *)
|
||||
|
||||
|
||||
(** Creation *)
|
||||
|
||||
(* Defaults:
|
||||
* - multiplicity : ~1~
|
||||
* - charge : ~0~
|
||||
* - operators : ~[]~ *)
|
||||
|
||||
|
||||
(* [[file:../simulation.org::*Creation][Creation:2]] *)
|
||||
let make
|
||||
?(multiplicity=1)
|
||||
?(charge=0)
|
||||
|
|
@ -58,18 +44,19 @@ let make
|
|||
Electrons.of_atoms ~multiplicity ~charge nuclei
|
||||
in
|
||||
|
||||
let charge =
|
||||
let charge =
|
||||
Charge.(Nuclei.charge nuclei + Electrons.charge electrons)
|
||||
in
|
||||
|
||||
{ charge ; nuclei ; electrons ; ao_basis ; operators}
|
||||
(* Creation:2 ends here *)
|
||||
|
||||
(* [[file:../simulation.org::*Printers][Printers:2]] *)
|
||||
|
||||
(** Printers *)
|
||||
|
||||
let pp ppf t =
|
||||
let formula = Nuclei.formula t.nuclei in
|
||||
let n_aos = Ao.Basis.size t.ao_basis in
|
||||
let n_ops = List.length t.operators in
|
||||
Format.fprintf ppf "@[@[%s@], @[%a@], @[%d AOs@], @[%d operators@]@]"
|
||||
formula Electrons.pp t.electrons n_aos n_ops
|
||||
(* Printers:2 ends here *)
|
||||
|
||||
|
|
|
|||
|
|
@ -1,50 +1,45 @@
|
|||
(* Simulation
|
||||
* :PROPERTIES:
|
||||
* :header-args: :noweb yes :comments both
|
||||
* :END:
|
||||
*
|
||||
* Contains the state of the simulation.
|
||||
*
|
||||
* #+NAME: open *)
|
||||
(** Contains the state of the simulation. *)
|
||||
|
||||
(* [[file:../simulation.org::open][open]] *)
|
||||
open Common
|
||||
open Particles
|
||||
open Operators
|
||||
(* open ends here *)
|
||||
|
||||
(* Type
|
||||
*
|
||||
* #+NAME: types *)
|
||||
|
||||
(* [[file:../simulation.org::types][types]] *)
|
||||
type t
|
||||
(* types ends here *)
|
||||
|
||||
(* Access *)
|
||||
|
||||
|
||||
(* [[file:../simulation.org::*Access][Access:1]] *)
|
||||
(** Access *)
|
||||
|
||||
val nuclei : t -> Nuclei.t
|
||||
(** Nuclear coordinates used in the smiulation *)
|
||||
|
||||
val charge : t -> Charge.t
|
||||
(** Total charge (electrons + nuclei) *)
|
||||
|
||||
val electrons : t -> Electrons.t
|
||||
(** Electrons used in the simulation *)
|
||||
|
||||
val ao_basis : t -> Ao.Basis.t
|
||||
(** Atomic basis set *)
|
||||
|
||||
val nuclear_repulsion : t -> float
|
||||
(** Nuclear repulsion energy *)
|
||||
|
||||
val operators : t -> Operator.t list
|
||||
(* Access:1 ends here *)
|
||||
|
||||
(* Creation *)
|
||||
(** List of extra operators (range-separation, f12, etc) *)
|
||||
|
||||
|
||||
(* [[file:../simulation.org::*Creation][Creation:1]] *)
|
||||
val make : ?multiplicity:int -> ?charge:int ->
|
||||
(** Creation *)
|
||||
|
||||
val make : ?multiplicity:int -> ?charge:int ->
|
||||
?operators:Operator.t list-> nuclei:Nuclei.t ->
|
||||
Ao.Basis.t -> t
|
||||
(* Creation:1 ends here *)
|
||||
|
||||
(* Printers *)
|
||||
(** Defaults:
|
||||
* - multiplicity : ~1~
|
||||
* - charge : ~0~
|
||||
* - operators : ~[]~
|
||||
*)
|
||||
|
||||
|
||||
(* [[file:../simulation.org::*Printers][Printers:1]] *)
|
||||
(** Printers *)
|
||||
|
||||
val pp : Format.formatter -> t -> unit
|
||||
(* Printers:1 ends here *)
|
||||
|
|
|
|||
|
|
@ -1,131 +0,0 @@
|
|||
#+begin_src elisp tangle: no :results none :exports none
|
||||
(setq pwd (file-name-directory buffer-file-name))
|
||||
(setq name (file-name-nondirectory (substring buffer-file-name 0 -4)))
|
||||
(setq lib (concat pwd "lib/"))
|
||||
(setq testdir (concat pwd "test/"))
|
||||
(setq mli (concat lib name ".mli"))
|
||||
(setq ml (concat lib name ".ml"))
|
||||
(setq test-ml (concat testdir name ".ml"))
|
||||
(org-babel-tangle)
|
||||
#+end_src
|
||||
|
||||
* Simulation
|
||||
:PROPERTIES:
|
||||
:header-args: :noweb yes :comments both
|
||||
:END:
|
||||
|
||||
Contains the state of the simulation.
|
||||
|
||||
#+NAME: open
|
||||
#+begin_src ocaml :tangle (eval mli)
|
||||
open Common
|
||||
open Particles
|
||||
open Operators
|
||||
#+end_src
|
||||
|
||||
#+begin_src ocaml :tangle (eval ml) :exports none
|
||||
<<open>>
|
||||
#+end_src
|
||||
|
||||
** Type
|
||||
|
||||
#+NAME: types
|
||||
#+begin_src ocaml :tangle (eval mli)
|
||||
type t
|
||||
#+end_src
|
||||
|
||||
#+begin_src ocaml :tangle (eval ml) :exports none
|
||||
type t = {
|
||||
charge : Charge.t;
|
||||
electrons : Electrons.t;
|
||||
nuclei : Nuclei.t;
|
||||
ao_basis : Ao.Basis.t;
|
||||
operators : Operator.t list;
|
||||
}
|
||||
#+end_src
|
||||
|
||||
** Access
|
||||
|
||||
#+begin_src ocaml :tangle (eval mli)
|
||||
val nuclei : t -> Nuclei.t
|
||||
val charge : t -> Charge.t
|
||||
val electrons : t -> Electrons.t
|
||||
val ao_basis : t -> Ao.Basis.t
|
||||
val nuclear_repulsion : t -> float
|
||||
val operators : t -> Operator.t list
|
||||
#+end_src
|
||||
|
||||
| ~nuclei~ | Nuclear coordinates used in the smiulation |
|
||||
| ~charge~ | Total charge (electrons + nuclei) |
|
||||
| ~electrons~ | Electrons used in the simulation |
|
||||
| ~ao_basis~ | Atomic basis set |
|
||||
| ~nuclear_repulsion~ | Nuclear repulsion energy |
|
||||
| ~operators~ | List of extra operators (range-separation, f12, etc) |
|
||||
|
||||
#+begin_src ocaml :tangle (eval ml) :exports none
|
||||
let nuclei t = t.nuclei
|
||||
let charge t = t.charge
|
||||
let electrons t = t.electrons
|
||||
let ao_basis t = t.ao_basis
|
||||
let nuclear_repulsion t = Nuclei.repulsion @@ nuclei t
|
||||
let operators t = t.operators
|
||||
#+end_src
|
||||
|
||||
** Creation
|
||||
|
||||
#+begin_src ocaml :tangle (eval mli)
|
||||
val make : ?multiplicity:int -> ?charge:int ->
|
||||
?operators:Operator.t list-> nuclei:Nuclei.t ->
|
||||
Ao.Basis.t -> t
|
||||
#+end_src
|
||||
|
||||
Defaults:
|
||||
- multiplicity : ~1~
|
||||
- charge : ~0~
|
||||
- operators : ~[]~
|
||||
|
||||
#+begin_src ocaml :tangle (eval ml) :exports none
|
||||
let make
|
||||
?(multiplicity=1)
|
||||
?(charge=0)
|
||||
?(operators=[])
|
||||
~nuclei
|
||||
ao_basis
|
||||
=
|
||||
|
||||
(* Tune Garbage Collector *)
|
||||
let gc = Gc.get () in
|
||||
Gc.set { gc with space_overhead = 1000 };
|
||||
|
||||
let electrons =
|
||||
Electrons.of_atoms ~multiplicity ~charge nuclei
|
||||
in
|
||||
|
||||
let charge =
|
||||
Charge.(Nuclei.charge nuclei + Electrons.charge electrons)
|
||||
in
|
||||
|
||||
{ charge ; nuclei ; electrons ; ao_basis ; operators}
|
||||
#+end_src
|
||||
|
||||
** Printers
|
||||
|
||||
#+begin_src ocaml :tangle (eval mli)
|
||||
val pp : Format.formatter -> t -> unit
|
||||
#+end_src
|
||||
|
||||
#+begin_src ocaml :tangle (eval ml) :exports none
|
||||
let pp ppf t =
|
||||
let formula = Nuclei.formula t.nuclei in
|
||||
let n_aos = Ao.Basis.size t.ao_basis in
|
||||
let n_ops = List.length t.operators in
|
||||
Format.fprintf ppf "@[@[%s@], @[%a@], @[%d AOs@], @[%d operators@]@]"
|
||||
formula Electrons.pp t.electrons n_aos n_ops
|
||||
#+end_src
|
||||
|
||||
#+RESULTS:
|
||||
: Line 2, characters 16-30:
|
||||
: 2 | let formula = Nuclei.formula t.nuclei in
|
||||
: ^^^^^^^^^^^^^^
|
||||
: Error: Unbound module Nuclei
|
||||
|
||||
|
|
@ -21,17 +21,17 @@
|
|||
|
||||
#+NAME: raw_printers
|
||||
#+begin_src sh :results output raw :exports none
|
||||
grep "val pp " ../*/*.org \
|
||||
| grep -v docs \
|
||||
grep "val pp " ../*/lib/*.mli \
|
||||
| cut -d ":" -f 1 \
|
||||
| sed 's|^../||' \
|
||||
| sed 's|.org||' \
|
||||
| sed 's|/| |' \
|
||||
| sed 's|.mli||' \
|
||||
| sed 's|/lib/| |' \
|
||||
| sed 's/[^ ]*/\u&/g' \
|
||||
| sed 's| |.|' \
|
||||
| grep -v Top.Install_printers
|
||||
#+end_src
|
||||
|
||||
|
||||
#+NAME: printers
|
||||
#+begin_src python :results output raw :var p=raw_printers :exports none
|
||||
result = []
|
||||
|
|
|
|||
|
|
@ -1,7 +1,28 @@
|
|||
(** Intall printers printers:3]] *)
|
||||
(* [[file:~/QCaml/top/install_printers.org::*Intall printers][Intall printers:3]] *)
|
||||
let printers =
|
||||
[
|
||||
"Ao.Basis.pp" ;
|
||||
"Ao.Basis_gaussian.pp" ;
|
||||
"Common.Angular_momentum.pp" ;
|
||||
"Common.Bitstring.pp" ;
|
||||
"Common.Charge.pp" ;
|
||||
"Common.Coordinate.pp" ;
|
||||
"Common.Powers.pp" ;
|
||||
"Common.Range.pp" ;
|
||||
"Common.Spin.pp" ;
|
||||
"Common.Zkey.pp" ;
|
||||
"Gaussian.Atomic_shell.pp" ;
|
||||
"Gaussian.Atomic_shell_pair.pp" ;
|
||||
"Gaussian.Atomic_shell_pair_couple.pp" ;
|
||||
"Gaussian.Basis.pp" ;
|
||||
"Gaussian.Contracted_shell.pp" ;
|
||||
"Gaussian.General_basis.pp" ;
|
||||
"Linear_algebra.Matrix.pp" ;
|
||||
"Mo.Basis.pp" ;
|
||||
"Mo.Class.pp" ;
|
||||
"Mo.Fock.pp" ;
|
||||
"Mo.Frozen_core.pp" ;
|
||||
"Mo.Hartree_fock.pp" ;
|
||||
"Mo.Localization.pp" ;
|
||||
"Particles.Electrons.pp" ;
|
||||
"Particles.Element.pp" ;
|
||||
|
|
@ -9,7 +30,7 @@ let printers =
|
|||
"Particles.Zmatrix.pp" ;
|
||||
"Perturbation.Mp2.pp" ;
|
||||
"Simulation.pp" ;
|
||||
|
||||
|
||||
]
|
||||
|
||||
let eval_exn str =
|
||||
|
|
@ -27,3 +48,4 @@ let rec install_printers = function
|
|||
let () =
|
||||
if not (install_printers printers) then
|
||||
Format.eprintf "Problem installing QCaml-printers@."
|
||||
(* Intall printers:3 ends here *)
|
||||
|
|
|
|||
Loading…
Reference in New Issue