16 KiB
16 KiB
None
<html>
<head>
</head>
</html>
Installation de QCaml¶
- Clonage de QCaml:
git clone https://gitlab.com/scemama/QCaml.git
- Installation des dependances:
opam install ocamlbuild ocamlfind lacaml gnuplot getopt alcotest zarith cd QCaml ./configure
- Compilation
make
Initialisation¶
Bloc a executer avant de pouvoir utiliser QCaml dans le Notebook
In [ ]:
let png_image = print_endline ;;
(* --------- *)
(*Mettre le bon chemin ici *)
#cd "/home/ydamour/QCaml";;
#use "topfind";;
#require "jupyter.notebook";;
#require "gnuplot";;
let png_image name =
Jupyter_notebook.display_file ~base64:true "image/png" (name)
;;
#require "lacaml.top";;
#require "alcotest";;
#require "str";;
#require "bigarray";;
#require "zarith";;
#require "getopt";;
#directory "_build";;
#directory "_build/Basis";;
#directory "_build/CI";;
#directory "_build/MOBasis";;
#directory "_build/Nuclei";;
#directory "_build/Parallel";;
#directory "_build/Perturbation";;
#directory "_build/SCF";;
#directory "_build/Utils";;
#load "Constants.cmo";;
#load_rec "Util.cma";;
#load_rec "Matrix.cmo";;
#load_rec "Simulation.cmo";;
#load_rec "Spindeterminant.cmo";;
#load_rec "Determinant.cmo";;
#load_rec "HartreeFock.cmo";;
#load_rec "MOBasis.cmo";;
#load_rec "F12CI.cmo";;
#install_printer AngularMomentum.pp_string ;;
#install_printer Basis.pp ;;
#install_printer Charge.pp ;;
#install_printer Coordinate.pp ;;
#install_printer Vector.pp;;
#install_printer Matrix.pp;;
#install_printer Util.pp_float_2darray;;
#install_printer Util.pp_float_array;;
#install_printer Util.pp_matrix;;
#install_printer HartreeFock.pp ;;
#install_printer Fock.pp ;;
#install_printer MOClass.pp ;;
#install_printer DeterminantSpace.pp;;
#install_printer SpindeterminantSpace.pp;;
#install_printer Bitstring.pp;;
let pp_mo ppf t = MOBasis.pp ~start:1 ~finish:0 ppf t ;;
#install_printer pp_mo;;
(* --------- *)
open Lacaml.D
Calculs¶
H_{n}¶
Generation du fichier xyz¶
Pour faire plusieurs chaines d'hydrogene, il faut d'abord
construire une fonction qui genere un fichier xyz pour
une chaine de n hydrogenes avec une distance H-H de d bohr.
Le fichier xyz est en Angstrom, donc il faut utiliser
Constants.a0 pour convertir les bohr en Angstroms.
Voila un exemple de ce qu'on doit obtenir pour 20 hydrogenes:
20
Hydrogen chain, d=1.8 Bohr
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.000000In [ ]:
let string_of_xyz d n =
" 10
Hydrogen chain, d=1.8 Bohr
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
"
let xyz_string =
string_of_xyz 1.8 10;;
Base atomique¶
Les bases atomiques sont telechargeables ici: https://www.basissetexchange.org/
Voila la base STO-6G pour l'hydrogene:
In [ ]:
let basis_string =
"
HYDROGEN
S 6
1 0.3552322122E+02 0.9163596281E-02
2 0.6513143725E+01 0.4936149294E-01
3 0.1822142904E+01 0.1685383049E+00
4 0.6259552659E+00 0.3705627997E+00
5 0.2430767471E+00 0.4164915298E+00
6 0.1001124280E+00 0.1303340841E+00
"
Une orbitale atomique STO-6G centree sur l'atome A est composee d'une contraction de 6 Gaussiennes: $$ \chi(r) = \sum_{i=1}^6 a_i \exp \left( -\alpha_i |r-r_A|^2 \right) $$ Dans le fichier de base, la 2ieme colonne represente l'exposant $\alpha_i$ et la 3ieme colonne le coefficient de contraction $a_i$.
In [ ]:
let nuclei =
Nuclei.of_xyz_string xyz_string
let basis =
Basis.of_nuclei_and_basis_string nuclei basis_string
let simulation =
Simulation.make ~charge:0 ~multiplicity:1 ~nuclei basis
let ao_basis =
Simulation.ao_basis simulation
Plot des orbitales atomiques¶
In [ ]:
let plot data filename =
let output = Gnuplot.Output.create (`Png filename) in
let gp = Gnuplot.create () in
Gnuplot.set gp ~use_grid:true;
List.map Gnuplot.Series.lines_xy data
|> Gnuplot.plot_many gp ~output;
Gnuplot.close gp ;
Jupyter_notebook.display_file ~base64:true "image/png" filename
;;
let x_values =
let n = 1000 in
let xmin, xmax =
let coord =
Array.map snd nuclei
|> Array.map (fun a -> Coordinate.(get X) a)
in
Array.sort compare coord;
coord.(0) -. 4. ,
coord.(Array.length coord -1) +. 4.
in
let dx =
(xmax -. xmin) /. (float_of_int n -. 1.)
in
Array.init n (fun i -> xmin +. (float_of_int i)*.dx)
in
let data =
Array.map (fun x ->
let point = Coordinate.make_angstrom
{ Coordinate.
x ; y = 0. ; z = 0.
} in
AOBasis.values ao_basis point
) x_values
|> Mat.of_col_vecs
|> Mat.transpose_copy
|> Mat.to_col_vecs
|> Array.map Vec.to_list
|> Array.map (fun l -> List.mapi (fun i y -> (x_values.(i),y)) l)
|> Array.to_list
in
plot data "test_data.png"
Calcul Hartree-Fock pour generer des orbitales moleculaires¶
In [ ]:
let hf = HartreeFock.make simulation
In [ ]:
let mo_basis = MOBasis.of_hartree_fock hf
Orbitales moleculaires : $$ \phi_j(r) = \sum_{i=1}^{N_b} C_{ij} \chi_i(r) $$
- $i$: lignes
- $j$: colonnes
Extraction des OM de la structure de donnees mo_basis comme une matrice $C$ utilisable avec Lacaml:
In [ ]:
let mo_coef = MOBasis.mo_coef mo_basis
In [ ]:
let plot data filename =
let output = Gnuplot.Output.create (`Png filename) in
let gp = Gnuplot.create () in
Gnuplot.set gp ~use_grid:true;
List.map Gnuplot.Series.lines_xy data
|> Gnuplot.plot_many gp ~output;
Gnuplot.close gp ;
Jupyter_notebook.display_file ~base64:true "image/png" filename
;;
let x_values =
let n = 1000 in
let xmin, xmax =
let coord =
Array.map snd nuclei
|> Array.map (fun a -> Coordinate.(get X) a)
in
Array.sort compare coord;
coord.(0) -. 4. ,
coord.(Array.length coord -1) +. 4.
in
let dx =
(xmax -. xmin) /. (float_of_int n -. 1.)
in
Array.init n (fun i -> xmin +. (float_of_int i)*.dx)
in
let data =
let result =
Array.map (fun x ->
let point = Coordinate.make_angstrom
{ Coordinate.
x ; y = 0. ; z = 0.
} in
MOBasis.values mo_basis point
) x_values
|> Mat.of_col_vecs
|> Mat.transpose_copy
|> Mat.to_col_vecs
|> Array.map Vec.to_list
|> Array.map (fun l -> List.mapi (fun i y -> (x_values.(i),y)) l)
in
[ result.(0) ; result.(1) ]
in
plot data "test_data.png"
Calcul des integrales necessaires¶
In [ ]:
let dx, x_values =
let n = 101 in
let xmin, xmax =
let coord =
Array.map snd nuclei
|> Array.map (fun a -> Coordinate.(get X) a)
in
Array.sort compare coord;
coord.(0) -. 10. ,
coord.(Array.length coord -1) +. 10.
in
let dx =
(xmax -. xmin) /. (float_of_int n -. 1.)
in
dx, Array.init n (fun i -> xmin +. (float_of_int i)*.dx)
in
let dy, y_values =
let n = 51 in
let ymin, ymax =
let coord =
Array.map snd nuclei
|> Array.map (fun a -> Coordinate.(get Y) a)
in
Array.sort compare coord;
coord.(0) -. 10. ,
coord.(Array.length coord -1) +. 10.
in
let dy =
(ymax -. ymin) /. (float_of_int n -. 1.)
in
dy, Array.init n (fun i -> ymin +. (float_of_int i)*.dy)
in
let dz, z_values =
let n = 51 in
let zmin, zmax =
let coord =
Array.map snd nuclei
|> Array.map (fun a -> Coordinate.(get Z) a)
in
Array.sort compare coord;
coord.(0) -. 10. ,
coord.(Array.length coord -1) +. 10.
in
let dz =
(zmax -. zmin) /. (float_of_int n -. 1.)
in
dz, Array.init n (fun i -> zmin +. (float_of_int i)*.dz)
in
let points =
Array.map (fun x ->
Array.map (fun y ->
Array.map (fun z ->
Coordinate.make_angstrom
{ Coordinate.
x ; y ; z
}
) x_values
) y_values
) z_values
|> Array.map Array.to_list
|> Array.to_list
|> List.concat
|> Array.concat
in
let dv = dx *. dy *. dz in
Printf.printf "%f %f %f %f\n%!" dx dy dz dv;
let ao_val =
Array.map (fun point ->
AOBasis.values ao_basis point
) points
|> Mat.of_col_vecs
in
(*
let r =
gemm ~transb:`T ao_val ao_val
in
Mat.scal dv r;
r;;
AOBasis.overlap ao_basis |> Overlap.matrix ;;
*)
let ao_x_val =
Array.map (fun point ->
let r = AOBasis.values ao_basis point in
scal (Coordinate.(get X) point) r;
r
) points
|> Mat.of_col_vecs
in
let r =
gemm ~transb:`T ao_val ao_x_val
in
Mat.scal dv r;
r;;
AOBasis.multipole ao_basis |> Multipole.matrix_x ;;
Localisation des orbitales¶
In [ ]:
let local;;