mirror of
https://gitlab.com/scemama/QCaml.git
synced 2024-08-30 00:03:42 +02:00
Added HartreeFock
This commit is contained in:
parent
e1da54cd67
commit
c23821e098
33
Basis/BasisLexer.mll
Normal file
33
Basis/BasisLexer.mll
Normal file
@ -0,0 +1,33 @@
|
|||||||
|
{
|
||||||
|
exception SyntaxError of string
|
||||||
|
|
||||||
|
open GamessParser
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
let eol = ['\n']
|
||||||
|
let white = [' ' '\t']+
|
||||||
|
let element = ['A'-'Z' 'a'-'z']+ white? eol
|
||||||
|
let ang_mom = ['S' 'P' 'D' 'F' 'G' 'H' 'I' 'J' 'K' 'L' 'M' 'N' 'O'
|
||||||
|
's' 'p' 'd' 'f' 'g' 'h' 'i' 'j' 'k' 'l' 'm' 'n' 'o' ]
|
||||||
|
white
|
||||||
|
let integer = ['0'-'9']+
|
||||||
|
let real = '-'? integer '.' integer (['e' 'E'] ('+'|'-')? integer)?
|
||||||
|
|
||||||
|
|
||||||
|
rule read_all_rule = parse
|
||||||
|
| eol { EOL }
|
||||||
|
| white { SPACE }
|
||||||
|
| ang_mom as a { ANG_MOM (a.[0]) }
|
||||||
|
| element as e { ELEMENT (String.trim e) }
|
||||||
|
| integer as i { INTEGER (int_of_string i) }
|
||||||
|
| real as f { FLOAT (float_of_string f) }
|
||||||
|
| eof { EOF }
|
||||||
|
|
||||||
|
|
||||||
|
{
|
||||||
|
let rec read_all lexbuf =
|
||||||
|
match read_all_rule lexbuf with
|
||||||
|
| SPACE -> read_all_rule lexbuf
|
||||||
|
| x -> x
|
||||||
|
}
|
127
Basis/ContractedShell.ml
Normal file
127
Basis/ContractedShell.ml
Normal file
@ -0,0 +1,127 @@
|
|||||||
|
open Util
|
||||||
|
open Constants
|
||||||
|
open Coordinate
|
||||||
|
|
||||||
|
type t = {
|
||||||
|
expo : float array; (* Gaussian exponents *)
|
||||||
|
coef : float array; (* Contraction coefficients *)
|
||||||
|
center : Coordinate.t; (* Center of all the Gaussians *)
|
||||||
|
totAngMom : AngularMomentum.t; (* Total angular momentum *)
|
||||||
|
size : int; (* Number of contracted Gaussians *)
|
||||||
|
norm_coef : float array; (* Normalization coefficient of the class
|
||||||
|
corresponding to the i-th contraction *)
|
||||||
|
norm_coef_scale : float array; (* Inside a class, the norm is the norm
|
||||||
|
of the function with (totAngMom,0,0) *.
|
||||||
|
this scaling factor *)
|
||||||
|
index : int; (* Index in the array of contracted shells *)
|
||||||
|
powers : Zkey.t array; (* Array of Zkeys corresponding to the
|
||||||
|
powers of (x,y,z) in the class *)
|
||||||
|
}
|
||||||
|
|
||||||
|
module Am = AngularMomentum
|
||||||
|
|
||||||
|
(** Normalization coefficient of contracted function i, which depends on the
|
||||||
|
exponent and the angular momentum. Two conventions can be chosen : a single
|
||||||
|
normalisation factor for all functions of the class, or a coefficient which
|
||||||
|
depends on the powers of x,y and z.
|
||||||
|
Returns, for each contracted function, an array of functions taking as
|
||||||
|
argument the [|x;y;z|] powers.
|
||||||
|
*)
|
||||||
|
let compute_norm_coef expo totAngMom =
|
||||||
|
let atot =
|
||||||
|
Am.to_int totAngMom
|
||||||
|
in
|
||||||
|
let factor int_array =
|
||||||
|
let dfa = Array.map (fun j ->
|
||||||
|
( float_of_int (1 lsl j) *. fact j) /. fact (j+j)
|
||||||
|
) int_array
|
||||||
|
in
|
||||||
|
sqrt (dfa.(0) *.dfa.(1) *. dfa.(2))
|
||||||
|
in
|
||||||
|
let expo =
|
||||||
|
if atot mod 2 = 0 then
|
||||||
|
Array.map (fun alpha ->
|
||||||
|
let alpha_2 = alpha +. alpha in
|
||||||
|
(alpha_2 *. pi_inv)**(0.75) *. (pow (alpha_2 +. alpha_2) (atot/2))
|
||||||
|
) expo
|
||||||
|
else
|
||||||
|
Array.map (fun alpha ->
|
||||||
|
let alpha_2 = alpha +. alpha in
|
||||||
|
(alpha_2 *. pi_inv)**(0.75) *. sqrt (pow (alpha_2 +. alpha_2) atot)
|
||||||
|
) expo
|
||||||
|
in
|
||||||
|
Array.map (fun x -> let f a = x *. (factor a) in f) expo
|
||||||
|
|
||||||
|
|
||||||
|
let make ~index ~expo ~coef ~center ~totAngMom =
|
||||||
|
assert (Array.length expo = Array.length coef);
|
||||||
|
assert (Array.length expo > 0);
|
||||||
|
let norm_coef_func =
|
||||||
|
compute_norm_coef expo totAngMom
|
||||||
|
in
|
||||||
|
let powers =
|
||||||
|
Am.zkey_array (Am.Singlet totAngMom)
|
||||||
|
in
|
||||||
|
let norm_coef =
|
||||||
|
Array.map (fun f -> f [| Am.to_int totAngMom ; 0 ; 0 |]) norm_coef_func
|
||||||
|
in
|
||||||
|
let norm_coef_scale =
|
||||||
|
Array.map (fun a ->
|
||||||
|
(norm_coef_func.(0) (Zkey.to_int_array ~kind:Zkey.Kind_3 a)) /. norm_coef.(0)
|
||||||
|
) powers
|
||||||
|
in
|
||||||
|
{ index ; expo ; coef ; center ; totAngMom ; size=Array.length expo ; norm_coef ;
|
||||||
|
norm_coef_scale ; powers }
|
||||||
|
|
||||||
|
|
||||||
|
let with_index a i =
|
||||||
|
{ a with index = i }
|
||||||
|
|
||||||
|
|
||||||
|
let to_string s =
|
||||||
|
let coord = s.center in
|
||||||
|
let open Printf in
|
||||||
|
(match s.totAngMom with
|
||||||
|
| Am.S -> sprintf "%3d " (s.index+1)
|
||||||
|
| _ -> sprintf "%3d-%-3d" (s.index+1) (s.index+(Array.length s.powers))
|
||||||
|
) ^
|
||||||
|
( sprintf "%1s %8.3f %8.3f %8.3f " (Am.to_string s.totAngMom)
|
||||||
|
(get X coord) (get Y coord) (get Z coord) ) ^
|
||||||
|
(Array.map2 (fun e c -> sprintf "%16.8e %16.8e" e c) s.expo s.coef
|
||||||
|
|> Array.to_list |> String.concat (sprintf "\n%36s" " ") )
|
||||||
|
|
||||||
|
|
||||||
|
(** Normalization coefficient of contracted function i, which depends on the
|
||||||
|
exponent and the angular momentum. Two conventions can be chosen : a single
|
||||||
|
normalisation factor for all functions of the class, or a coefficient which
|
||||||
|
depends on the powers of x,y and z.
|
||||||
|
Returns, for each contracted function, an array of functions taking as
|
||||||
|
argument the [|x;y;z|] powers.
|
||||||
|
*)
|
||||||
|
let compute_norm_coef expo totAngMom =
|
||||||
|
let atot =
|
||||||
|
Am.to_int totAngMom
|
||||||
|
in
|
||||||
|
let factor int_array =
|
||||||
|
let dfa = Array.map (fun j ->
|
||||||
|
(float_of_int (1 lsl j) *. fact j) /. fact (j+j)
|
||||||
|
) int_array
|
||||||
|
in
|
||||||
|
sqrt (dfa.(0) *.dfa.(1) *. dfa.(2))
|
||||||
|
in
|
||||||
|
let expo =
|
||||||
|
if atot mod 2 = 0 then
|
||||||
|
Array.map (fun alpha ->
|
||||||
|
let alpha_2 = alpha +. alpha in
|
||||||
|
(alpha_2 *. pi_inv)**(0.75) *. (pow (alpha_2 +. alpha_2) (atot/2))
|
||||||
|
) expo
|
||||||
|
else
|
||||||
|
Array.map (fun alpha ->
|
||||||
|
let alpha_2 = alpha +. alpha in
|
||||||
|
(alpha_2 *. pi_inv)**(0.75) *. sqrt (pow (alpha_2 +. alpha_2) atot)
|
||||||
|
) expo
|
||||||
|
in
|
||||||
|
Array.map (fun x -> let f a = x *. factor a in f) expo
|
||||||
|
|
||||||
|
|
||||||
|
|
26
Basis/ContractedShell.mli
Normal file
26
Basis/ContractedShell.mli
Normal file
@ -0,0 +1,26 @@
|
|||||||
|
type t = private {
|
||||||
|
expo : float array;
|
||||||
|
coef : float array;
|
||||||
|
center : Coordinate.t;
|
||||||
|
totAngMom : AngularMomentum.t;
|
||||||
|
size : int;
|
||||||
|
norm_coef : float array;
|
||||||
|
norm_coef_scale : float array;
|
||||||
|
index : int;
|
||||||
|
powers : Zkey.t array;
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
(** Pretty-printing of the contracted shell in a string *)
|
||||||
|
val to_string : t -> string
|
||||||
|
|
||||||
|
(** Creates a contracted shell *)
|
||||||
|
val make :
|
||||||
|
index:int ->
|
||||||
|
expo:float array ->
|
||||||
|
coef:float array ->
|
||||||
|
center:Coordinate.t -> totAngMom:AngularMomentum.t -> t
|
||||||
|
|
||||||
|
(** Returns a copy of the contracted shell with a modified index *)
|
||||||
|
val with_index : t -> int -> t
|
49
Basis/GamessParser.mly
Normal file
49
Basis/GamessParser.mly
Normal file
@ -0,0 +1,49 @@
|
|||||||
|
/* Parses basis sets GAMESS format */
|
||||||
|
|
||||||
|
%{
|
||||||
|
|
||||||
|
%}
|
||||||
|
|
||||||
|
%token <string> ELEMENT
|
||||||
|
%token <char> ANG_MOM
|
||||||
|
%token <int> INTEGER
|
||||||
|
%token <float> FLOAT
|
||||||
|
%token SPACE
|
||||||
|
%token EOL
|
||||||
|
%token EOF
|
||||||
|
|
||||||
|
%start input
|
||||||
|
%type <GeneralBasis.t> input
|
||||||
|
|
||||||
|
%% /* Grammar rules and actions follow */
|
||||||
|
|
||||||
|
input:
|
||||||
|
| basis { $1 }
|
||||||
|
| EOL input { $2 }
|
||||||
|
|
||||||
|
basis:
|
||||||
|
| element shell_array EOL { ($1, $2) }
|
||||||
|
| element shell_array EOF { ($1, $2) }
|
||||||
|
|
||||||
|
element:
|
||||||
|
| ELEMENT { Element.of_string $1 }
|
||||||
|
|
||||||
|
shell_array:
|
||||||
|
| shell_list { Array.of_list @@ List.rev $1 }
|
||||||
|
|
||||||
|
shell_list:
|
||||||
|
| { [] }
|
||||||
|
| shell_list shell { $2 :: $1 }
|
||||||
|
|
||||||
|
shell:
|
||||||
|
| ANG_MOM INTEGER EOL primitive_list { (AngularMomentum.of_char $1, Array.of_list @@ List.rev $4 ) }
|
||||||
|
|
||||||
|
primitive_list:
|
||||||
|
| { [] }
|
||||||
|
| primitive_list primitive { $2 :: $1 }
|
||||||
|
|
||||||
|
primitive:
|
||||||
|
| INTEGER FLOAT FLOAT EOL { GeneralBasis.{exponent=$2 ; coefficient=$3 } }
|
||||||
|
|
||||||
|
|
||||||
|
|
18
Basis/GamessReader.ml
Normal file
18
Basis/GamessReader.ml
Normal file
@ -0,0 +1,18 @@
|
|||||||
|
(** Read a basis set file in GAMESS format and return an association list where the key is an
|
||||||
|
Element.t and the value is the parsed basis set. *)
|
||||||
|
let read ~filename =
|
||||||
|
let lexbuf =
|
||||||
|
let ic = open_in filename in
|
||||||
|
Lexing.from_channel ic
|
||||||
|
in
|
||||||
|
let rec aux accu =
|
||||||
|
try
|
||||||
|
let key, basis =
|
||||||
|
GamessParser.input BasisLexer.read_all lexbuf
|
||||||
|
in
|
||||||
|
aux ((key, basis)::accu)
|
||||||
|
with
|
||||||
|
| Parsing.Parse_error -> List.rev accu
|
||||||
|
in
|
||||||
|
aux []
|
||||||
|
|
40
Basis/GeneralBasis.ml
Normal file
40
Basis/GeneralBasis.ml
Normal file
@ -0,0 +1,40 @@
|
|||||||
|
(** General basis set read from a file *)
|
||||||
|
type primitive = {
|
||||||
|
exponent: float ;
|
||||||
|
coefficient: float
|
||||||
|
}
|
||||||
|
|
||||||
|
type general_contracted_shell = AngularMomentum.t * (primitive array)
|
||||||
|
|
||||||
|
type t = Element.t * (general_contracted_shell array)
|
||||||
|
|
||||||
|
|
||||||
|
module Am = AngularMomentum
|
||||||
|
|
||||||
|
let string_of_primitive ?id prim =
|
||||||
|
match id with
|
||||||
|
| None -> (string_of_float prim.exponent)^" "^(string_of_float prim.coefficient)
|
||||||
|
| Some i -> (string_of_int i)^" "^(string_of_float prim.exponent)^" "^(string_of_float prim.coefficient)
|
||||||
|
|
||||||
|
|
||||||
|
let string_of_contracted_shell (angular_momentum, prim_array) =
|
||||||
|
let n =
|
||||||
|
Array.length prim_array
|
||||||
|
in
|
||||||
|
Printf.sprintf "%s %d\n%s"
|
||||||
|
(Am.to_string angular_momentum) n
|
||||||
|
(Array.init n (fun i -> string_of_primitive ~id:(i+1) prim_array.(i))
|
||||||
|
|> Array.to_list
|
||||||
|
|> String.concat "\n")
|
||||||
|
|
||||||
|
|
||||||
|
let string_of_contracted_shell_array a =
|
||||||
|
Array.map string_of_contracted_shell a
|
||||||
|
|> Array.to_list
|
||||||
|
|> String.concat "\n"
|
||||||
|
|
||||||
|
|
||||||
|
let to_string (name, contracted_shell_array) =
|
||||||
|
Printf.sprintf "%s\n%s" name (string_of_contracted_shell_array contracted_shell_array)
|
||||||
|
|
||||||
|
|
29
Basis/Orthonormalization.ml
Normal file
29
Basis/Orthonormalization.ml
Normal file
@ -0,0 +1,29 @@
|
|||||||
|
open Util
|
||||||
|
open Lacaml.D
|
||||||
|
|
||||||
|
type t =
|
||||||
|
| Lowdin of Mat.t
|
||||||
|
| Canonical of Mat.t
|
||||||
|
| Svd of Mat.t
|
||||||
|
|
||||||
|
|
||||||
|
let make_lowdin ?(thresh=1.e-12) ~overlap =
|
||||||
|
|
||||||
|
let u_vec, u_val = diagonalize_symm overlap in
|
||||||
|
|
||||||
|
Vec.iter (fun x -> if x < thresh then
|
||||||
|
invalid_arg "Orthonormalization.make_lowdin") u_val;
|
||||||
|
|
||||||
|
let u_val = Vec.reci (Vec.sqrt u_val) in
|
||||||
|
|
||||||
|
let u_vec' =
|
||||||
|
Mat.init_cols (Mat.dim1 u_vec) (Mat.dim2 u_vec) (fun i j -> u_vec.{i,j} *. u_val.{j})
|
||||||
|
in
|
||||||
|
let result =
|
||||||
|
gemm u_vec' ~transb:`T u_vec
|
||||||
|
in
|
||||||
|
|
||||||
|
Lowdin result
|
||||||
|
|
||||||
|
|
||||||
|
let make = make_lowdin
|
33
HartreeFock/Fock.ml
Normal file
33
HartreeFock/Fock.ml
Normal file
@ -0,0 +1,33 @@
|
|||||||
|
open Lacaml.D
|
||||||
|
open Simulation
|
||||||
|
|
||||||
|
type t = Mat.t
|
||||||
|
|
||||||
|
let make ~density simulation =
|
||||||
|
let m_P = density
|
||||||
|
and m_T = Lazy.force simulation.kin_ints
|
||||||
|
and m_V = Lazy.force simulation.eN_ints
|
||||||
|
and m_G = Lazy.force simulation.ee_ints
|
||||||
|
in
|
||||||
|
let nBas = Mat.dim1 m_T
|
||||||
|
in
|
||||||
|
|
||||||
|
let m_F = Mat.add m_T m_V in
|
||||||
|
for sigma = 1 to nBas do
|
||||||
|
for nu = 1 to nBas do
|
||||||
|
for lambda = 1 to nBas do
|
||||||
|
let p = m_P.{lambda,sigma} in
|
||||||
|
for mu = 1 to nu do
|
||||||
|
m_F.{mu,nu} <- m_F.{mu,nu} +. p *.
|
||||||
|
(m_G.{mu,lambda,nu,sigma} -. 0.5 *. m_G.{mu,lambda,sigma,nu})
|
||||||
|
done
|
||||||
|
done
|
||||||
|
done
|
||||||
|
done;
|
||||||
|
for nu = 1 to nBas do
|
||||||
|
for mu = 1 to nu do
|
||||||
|
m_F.{nu,mu} <- m_F.{mu,nu}
|
||||||
|
done
|
||||||
|
done;
|
||||||
|
m_F
|
||||||
|
|
17
HartreeFock/Guess.ml
Normal file
17
HartreeFock/Guess.ml
Normal file
@ -0,0 +1,17 @@
|
|||||||
|
open Lacaml.D
|
||||||
|
|
||||||
|
type guess =
|
||||||
|
| Hcore of Mat.t
|
||||||
|
|
||||||
|
type t = guess
|
||||||
|
|
||||||
|
|
||||||
|
let make ?guess:(guess=`Hcore) simulation =
|
||||||
|
let eN_ints = Lazy.force simulation.Simulation.eN_ints
|
||||||
|
and kin_ints = Lazy.force simulation.Simulation.kin_ints
|
||||||
|
in
|
||||||
|
match guess with
|
||||||
|
| `Hcore -> Hcore (Mat.add eN_ints kin_ints)
|
||||||
|
|
||||||
|
|
||||||
|
|
12
HartreeFock/Guess.mli
Normal file
12
HartreeFock/Guess.mli
Normal file
@ -0,0 +1,12 @@
|
|||||||
|
(** Guess for Hartree-Fock calculations. *)
|
||||||
|
|
||||||
|
type guess =
|
||||||
|
| Hcore of Lacaml.D.Mat.t
|
||||||
|
|
||||||
|
type t = guess
|
||||||
|
|
||||||
|
|
||||||
|
val make : ?guess:[ `Hcore ] -> Simulation.t -> t
|
||||||
|
|
||||||
|
|
||||||
|
|
9
HartreeFock/HartreeFock.ml
Normal file
9
HartreeFock/HartreeFock.ml
Normal file
@ -0,0 +1,9 @@
|
|||||||
|
open Util
|
||||||
|
open Simulation
|
||||||
|
|
||||||
|
let make ?guess simulation =
|
||||||
|
if simulation.electrons.Electrons.multiplicity = 1 then
|
||||||
|
RHF.make ?guess simulation
|
||||||
|
else
|
||||||
|
invalid_arg "UHF or ROHF not implemented"
|
||||||
|
|
11
HartreeFock/HartreeFock_type.ml
Normal file
11
HartreeFock/HartreeFock_type.ml
Normal file
@ -0,0 +1,11 @@
|
|||||||
|
type t =
|
||||||
|
{
|
||||||
|
guess : Guess.t;
|
||||||
|
eigenvectors : Lacaml.D.Mat.t ;
|
||||||
|
eigenvalues : Lacaml.D.Vec.t ;
|
||||||
|
energy : float ;
|
||||||
|
iterations : float array;
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
|
10
HartreeFock/HartreeFock_type.mli
Normal file
10
HartreeFock/HartreeFock_type.mli
Normal file
@ -0,0 +1,10 @@
|
|||||||
|
type t =
|
||||||
|
{
|
||||||
|
guess : Guess.t;
|
||||||
|
eigenvectors : Lacaml.D.Mat.t ;
|
||||||
|
eigenvalues : Lacaml.D.Vec.t ;
|
||||||
|
energy : float ;
|
||||||
|
iterations : float array;
|
||||||
|
}
|
||||||
|
|
||||||
|
|
126
HartreeFock/RHF.ml
Normal file
126
HartreeFock/RHF.ml
Normal file
@ -0,0 +1,126 @@
|
|||||||
|
open Util
|
||||||
|
open Lacaml.D
|
||||||
|
open Simulation
|
||||||
|
|
||||||
|
let make ?guess:(guess=`Hcore) ?max_scf:(max_scf=64)
|
||||||
|
?threshold_SCF:(threshold_SCF=1.e-6) simulation =
|
||||||
|
|
||||||
|
(* Number of occupied MOs *)
|
||||||
|
let nocc =
|
||||||
|
simulation.electrons.Electrons.n_alpha
|
||||||
|
in
|
||||||
|
|
||||||
|
(* Initial guess *)
|
||||||
|
let guess =
|
||||||
|
Guess.make ~guess simulation
|
||||||
|
in
|
||||||
|
|
||||||
|
(* Orthogonalization matrix *)
|
||||||
|
let m_X =
|
||||||
|
match Lazy.force simulation.overlap_ortho with
|
||||||
|
| Lowdin x -> x
|
||||||
|
| Canonical x -> x
|
||||||
|
| Svd x -> x
|
||||||
|
in
|
||||||
|
|
||||||
|
(* Core Hamiltonian *)
|
||||||
|
let m_Hc =
|
||||||
|
let m_T = Lazy.force simulation.kin_ints
|
||||||
|
and m_V = Lazy.force simulation.eN_ints
|
||||||
|
in Mat.add m_T m_V
|
||||||
|
in
|
||||||
|
|
||||||
|
(* Overlap matrix *)
|
||||||
|
let m_S =
|
||||||
|
Lazy.force simulation.overlap
|
||||||
|
in
|
||||||
|
|
||||||
|
|
||||||
|
(* SCF iterations *)
|
||||||
|
let rec loop nSCF iterations m_C =
|
||||||
|
|
||||||
|
(* Density matrix over nocc occupied MOs *)
|
||||||
|
let m_P =
|
||||||
|
gemm ~alpha:2. ~transb:`T ~k:nocc m_C m_C
|
||||||
|
in
|
||||||
|
|
||||||
|
(* Fock matrix in AO basis *)
|
||||||
|
let m_F =
|
||||||
|
Fock.make ~density:m_P simulation
|
||||||
|
in
|
||||||
|
|
||||||
|
(* Fock matrix in MO basis *)
|
||||||
|
let m_Fmo =
|
||||||
|
xt_o_x m_F m_X
|
||||||
|
in
|
||||||
|
|
||||||
|
(* MOs in old MO basis *)
|
||||||
|
let m_C', eigenvalues =
|
||||||
|
diagonalize_symm m_Fmo
|
||||||
|
in
|
||||||
|
|
||||||
|
(* MOs in AO basis *)
|
||||||
|
let m_C =
|
||||||
|
gemm m_X m_C'
|
||||||
|
in
|
||||||
|
|
||||||
|
(* Hartree-Fock energy *)
|
||||||
|
let energy =
|
||||||
|
simulation.nuclear_repulsion +. 0.5 *.
|
||||||
|
Mat.gemm_trace m_P (Mat.add m_Hc m_F)
|
||||||
|
in
|
||||||
|
|
||||||
|
(* Convergence criterion *)
|
||||||
|
let commutator =
|
||||||
|
let fps =
|
||||||
|
gemm m_F (gemm m_P m_S)
|
||||||
|
and spf =
|
||||||
|
gemm m_S (gemm m_P m_F)
|
||||||
|
in
|
||||||
|
Mat.sub fps spf
|
||||||
|
|> Mat.as_vec
|
||||||
|
|> amax
|
||||||
|
in
|
||||||
|
|
||||||
|
let converged =
|
||||||
|
nSCF = max_scf || (abs_float commutator) < threshold_SCF
|
||||||
|
in
|
||||||
|
|
||||||
|
Printf.printf "%d %16.10f %10.4e\n%!" nSCF energy commutator;
|
||||||
|
|
||||||
|
if not converged then
|
||||||
|
loop (nSCF+1) (energy :: iterations) m_C
|
||||||
|
else
|
||||||
|
let iterations =
|
||||||
|
List.rev iterations
|
||||||
|
|> Array.of_list
|
||||||
|
in
|
||||||
|
{ HartreeFock_type.
|
||||||
|
guess ;
|
||||||
|
eigenvectors = m_C ;
|
||||||
|
eigenvalues ;
|
||||||
|
energy = iterations.(Array.length iterations - 1) ;
|
||||||
|
iterations ;
|
||||||
|
}
|
||||||
|
in
|
||||||
|
|
||||||
|
|
||||||
|
(* Guess coefficients *)
|
||||||
|
let m_H =
|
||||||
|
match guess with
|
||||||
|
| Guess.Hcore m_H -> m_H
|
||||||
|
in
|
||||||
|
let m_Hmo =
|
||||||
|
xt_o_x m_H m_X
|
||||||
|
in
|
||||||
|
let m_C', _ =
|
||||||
|
diagonalize_symm m_Hmo
|
||||||
|
in
|
||||||
|
let m_C =
|
||||||
|
gemm m_X m_C'
|
||||||
|
in
|
||||||
|
|
||||||
|
loop 1 [] m_C
|
||||||
|
|
||||||
|
|
||||||
|
|
134
Utils/AngularMomentum.ml
Normal file
134
Utils/AngularMomentum.ml
Normal file
@ -0,0 +1,134 @@
|
|||||||
|
open Powers
|
||||||
|
|
||||||
|
exception AngularMomentumError of string
|
||||||
|
|
||||||
|
type t =
|
||||||
|
| S | P | D | F | G | H | I | J | K | L | M | N | O
|
||||||
|
|
||||||
|
let of_char = function
|
||||||
|
| 's' | 'S' -> S | 'p' | 'P' -> P
|
||||||
|
| 'd' | 'D' -> D | 'f' | 'F' -> F
|
||||||
|
| 'g' | 'G' -> G | 'h' | 'H' -> H
|
||||||
|
| 'i' | 'I' -> I | 'j' | 'J' -> J
|
||||||
|
| 'k' | 'K' -> K | 'l' | 'L' -> L
|
||||||
|
| 'm' | 'M' -> M | 'n' | 'N' -> N
|
||||||
|
| 'o' | 'O' -> O
|
||||||
|
| c -> raise (AngularMomentumError (String.make 1 c))
|
||||||
|
|
||||||
|
let to_string = function
|
||||||
|
| S -> "S" | P -> "P"
|
||||||
|
| D -> "D" | F -> "F"
|
||||||
|
| G -> "G" | H -> "H"
|
||||||
|
| I -> "I" | J -> "J"
|
||||||
|
| K -> "K" | L -> "L"
|
||||||
|
| M -> "M" | N -> "N"
|
||||||
|
| O -> "O"
|
||||||
|
|
||||||
|
let to_char = function
|
||||||
|
| S -> 'S' | P -> 'P'
|
||||||
|
| D -> 'D' | F -> 'F'
|
||||||
|
| G -> 'G' | H -> 'H'
|
||||||
|
| I -> 'I' | J -> 'J'
|
||||||
|
| K -> 'K' | L -> 'L'
|
||||||
|
| M -> 'M' | N -> 'N'
|
||||||
|
| O -> 'O'
|
||||||
|
|
||||||
|
let to_int = function
|
||||||
|
| S -> 0 | P -> 1
|
||||||
|
| D -> 2 | F -> 3
|
||||||
|
| G -> 4 | H -> 5
|
||||||
|
| I -> 6 | J -> 7
|
||||||
|
| K -> 8 | L -> 9
|
||||||
|
| M -> 10 | N -> 11
|
||||||
|
| O -> 12
|
||||||
|
|
||||||
|
let of_int = function
|
||||||
|
| 0 -> S | 1 -> P
|
||||||
|
| 2 -> D | 3 -> F
|
||||||
|
| 4 -> G | 5 -> H
|
||||||
|
| 6 -> I | 7 -> J
|
||||||
|
| 8 -> K | 9 -> L
|
||||||
|
| 10 -> M | 11 -> N
|
||||||
|
| 12 -> O
|
||||||
|
| c -> raise (AngularMomentumError (string_of_int c))
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
type kind =
|
||||||
|
| Singlet of t
|
||||||
|
| Doublet of (t*t)
|
||||||
|
| Triplet of (t*t*t)
|
||||||
|
| Quartet of (t*t*t*t)
|
||||||
|
|
||||||
|
|
||||||
|
let n_functions a =
|
||||||
|
let a =
|
||||||
|
to_int a
|
||||||
|
in
|
||||||
|
(a*a + 3*a + 2)/2
|
||||||
|
|
||||||
|
|
||||||
|
(** Returns an array of Zkeys corresponding to all possible angular momenta *)
|
||||||
|
let zkey_array a =
|
||||||
|
let keys_1d l =
|
||||||
|
let create_z { x ; y ; _ } =
|
||||||
|
Powers.of_int_tuple (x,y,l-(x+y))
|
||||||
|
in
|
||||||
|
let rec create_y accu xyz =
|
||||||
|
let { x ; y ; z } = xyz in
|
||||||
|
match y with
|
||||||
|
| 0 -> (create_z xyz)::accu
|
||||||
|
| i -> let ynew = y-1 in
|
||||||
|
create_y ( (create_z xyz)::accu) (Powers.of_int_tuple (x,ynew,z))
|
||||||
|
in
|
||||||
|
let rec create_x accu xyz =
|
||||||
|
let { x ; y ; z } = xyz in
|
||||||
|
match x with
|
||||||
|
| 0 -> (create_y [] xyz)@accu
|
||||||
|
| i -> let xnew = x-1 in
|
||||||
|
let ynew = l-xnew in
|
||||||
|
create_x ((create_y [] xyz)@accu) (Powers.of_int_tuple (xnew, ynew, z))
|
||||||
|
in
|
||||||
|
create_x [] (Powers.of_int_tuple (l,0,0))
|
||||||
|
|> List.rev
|
||||||
|
in
|
||||||
|
|
||||||
|
begin
|
||||||
|
match a with
|
||||||
|
| Singlet l1 ->
|
||||||
|
List.map (fun x -> Zkey.of_powers (Zkey.Three x)) (keys_1d @@ to_int l1)
|
||||||
|
|
||||||
|
| Doublet (l1, l2) ->
|
||||||
|
List.map (fun a ->
|
||||||
|
List.map (fun b ->
|
||||||
|
Zkey.of_powers (Zkey.Six (a,b))) (keys_1d @@ to_int l2)
|
||||||
|
) (keys_1d @@ to_int l1)
|
||||||
|
|> List.concat
|
||||||
|
|
||||||
|
| Triplet (l1, l2, l3) ->
|
||||||
|
|
||||||
|
List.map (fun a ->
|
||||||
|
List.map (fun b ->
|
||||||
|
List.map (fun c ->
|
||||||
|
Zkey.of_powers (Zkey.Nine (a,b,c))) (keys_1d @@ to_int l3)
|
||||||
|
) (keys_1d @@ to_int l2)
|
||||||
|
|> List.concat
|
||||||
|
) (keys_1d @@ to_int l1)
|
||||||
|
|> List.concat
|
||||||
|
|
||||||
|
| Quartet (l1, l2, l3, l4) ->
|
||||||
|
|
||||||
|
List.map (fun a ->
|
||||||
|
List.map (fun b ->
|
||||||
|
List.map (fun c ->
|
||||||
|
List.map (fun d ->
|
||||||
|
Zkey.of_powers (Zkey.Twelve (a,b,c,d))) (keys_1d @@ to_int l4)
|
||||||
|
) (keys_1d @@ to_int l3)
|
||||||
|
|> List.concat
|
||||||
|
) (keys_1d @@ to_int l2)
|
||||||
|
|> List.concat
|
||||||
|
) (keys_1d @@ to_int l1)
|
||||||
|
|> List.concat
|
||||||
|
end
|
||||||
|
|> Array.of_list
|
||||||
|
|
23
Utils/Electrons.ml
Normal file
23
Utils/Electrons.ml
Normal file
@ -0,0 +1,23 @@
|
|||||||
|
(** Number of alpha and beta electrons *)
|
||||||
|
|
||||||
|
type t = {
|
||||||
|
n_alpha : int ;
|
||||||
|
n_beta : int ;
|
||||||
|
multiplicity : int;
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
let make ?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_alpha = n_elec - n_beta in
|
||||||
|
let result = { n_alpha ; n_beta ; multiplicity } in
|
||||||
|
if multiplicity <> (n_alpha - n_beta)+1 then
|
||||||
|
invalid_arg "Electrons.make";
|
||||||
|
result
|
||||||
|
|
14
Utils/PositiveFloat.ml
Normal file
14
Utils/PositiveFloat.ml
Normal file
@ -0,0 +1,14 @@
|
|||||||
|
type t = float
|
||||||
|
|
||||||
|
let of_float x =
|
||||||
|
assert ( x >= 0. );
|
||||||
|
x
|
||||||
|
|
||||||
|
external to_float : t -> float = "%identity"
|
||||||
|
|
||||||
|
let to_string x =
|
||||||
|
let f = to_float x in string_of_float f
|
||||||
|
|
||||||
|
let of_string x =
|
||||||
|
let f = float_of_string x in of_float f
|
||||||
|
|
5
Utils/PositiveFloat.mli
Normal file
5
Utils/PositiveFloat.mli
Normal file
@ -0,0 +1,5 @@
|
|||||||
|
type t = private float
|
||||||
|
val of_float : float -> t
|
||||||
|
val to_float : t -> float
|
||||||
|
val to_string : t -> string
|
||||||
|
val of_string : string -> t
|
53
run_hartree_fock.ml
Normal file
53
run_hartree_fock.ml
Normal file
@ -0,0 +1,53 @@
|
|||||||
|
let out_file : string option ref = ref None
|
||||||
|
let basis_file : string option ref = ref None
|
||||||
|
let nuclei_file : string option ref = ref None
|
||||||
|
let charge : int ref = ref 0
|
||||||
|
let multiplicity: int ref = ref 1
|
||||||
|
|
||||||
|
|
||||||
|
let speclist = [
|
||||||
|
( "-b" , Arg.String (fun x -> basis_file := Some x),
|
||||||
|
"File containing the atomic basis set") ;
|
||||||
|
( "-c" , Arg.String (fun x -> nuclei_file := Some x),
|
||||||
|
"File containing the nuclear coordinates") ;
|
||||||
|
( "-o" , Arg.String (fun x -> out_file := Some x) ,
|
||||||
|
"Output file") ;
|
||||||
|
( "-charge" , Arg.Int (fun x -> charge := x),
|
||||||
|
"Charge of the system") ;
|
||||||
|
( "-mult" , Arg.Int (fun x -> multiplicity := x),
|
||||||
|
"Spin multiplicity of the system") ;
|
||||||
|
]
|
||||||
|
|
||||||
|
let run ~out =
|
||||||
|
(*
|
||||||
|
let gc = Gc.get () in
|
||||||
|
Gc.set { gc with minor_heap_size=(262144 / 16) };
|
||||||
|
*)
|
||||||
|
let basis_file =
|
||||||
|
match !basis_file with
|
||||||
|
| None -> raise (Invalid_argument "Basis set file should be specified with -b")
|
||||||
|
| Some x -> x
|
||||||
|
and nuclei_file =
|
||||||
|
match !nuclei_file with
|
||||||
|
| None -> raise (Invalid_argument "Basis set file should be specified with -c")
|
||||||
|
| Some x -> x
|
||||||
|
and charge = !charge
|
||||||
|
and multiplicity = !multiplicity
|
||||||
|
in
|
||||||
|
|
||||||
|
let s =
|
||||||
|
Simulation.of_filenames ~charge ~multiplicity ~nuclei:nuclei_file basis_file
|
||||||
|
in
|
||||||
|
|
||||||
|
let _ =
|
||||||
|
HartreeFock.make s
|
||||||
|
in
|
||||||
|
Printf.printf "Done.\n%!";
|
||||||
|
()
|
||||||
|
|
||||||
|
|
||||||
|
let () =
|
||||||
|
let usage_msg = "Available options:" in
|
||||||
|
Arg.parse speclist (fun _ -> ()) usage_msg;
|
||||||
|
run ~out:!out_file
|
||||||
|
|
Loading…
Reference in New Issue
Block a user