24 KiB
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
Liens utiles¶
- Documentation des modules : https://sanette.github.io/ocaml-api/4.10/index.html
- Lacaml : http://mmottl.github.io/lacaml/api/lacaml/
- Gnuplot : https://github.com/c-cube/ocaml-gnuplot
Initialisation¶
Bloc a executer avant de pouvoir utiliser QCaml dans le Notebook
let png_image = print_endline ;;
(* --------- *)
(*Mettre le bon chemin ici *)
#cd "/home/scemama/QCaml";;
#use "topfind";;
#require "jupyter.notebook";;
#require "gnuplot";;
let png_image name =
Jupyter_notebook.display_file ~base64:true "image/png" ("Notebooks/images/"^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$_2$O¶
let xyz_string =
"3
Water
O 0. 0. 0.
H -0.756950272703377558 0. -0.585882234512562827
H 0.756950272703377558 0. -0.585882234512562827
"
Base atomique¶
Les bases atomiques sont telechargeables ici: https://www.basissetexchange.org/
On telecharge la base cc-pvdz depuis le site:
let basis_string = "
HYDROGEN
S 4
1 1.301000E+01 1.968500E-02
2 1.962000E+00 1.379770E-01
3 4.446000E-01 4.781480E-01
4 1.220000E-01 5.012400E-01
S 1
1 1.220000E-01 1.000000E+00
P 1
1 7.270000E-01 1.0000000
OXYGEN
S 9
1 1.172000E+04 7.100000E-04
2 1.759000E+03 5.470000E-03
3 4.008000E+02 2.783700E-02
4 1.137000E+02 1.048000E-01
5 3.703000E+01 2.830620E-01
6 1.327000E+01 4.487190E-01
7 5.025000E+00 2.709520E-01
8 1.013000E+00 1.545800E-02
9 3.023000E-01 -2.585000E-03
S 9
1 1.172000E+04 -1.600000E-04
2 1.759000E+03 -1.263000E-03
3 4.008000E+02 -6.267000E-03
4 1.137000E+02 -2.571600E-02
5 3.703000E+01 -7.092400E-02
6 1.327000E+01 -1.654110E-01
7 5.025000E+00 -1.169550E-01
8 1.013000E+00 5.573680E-01
9 3.023000E-01 5.727590E-01
S 1
1 3.023000E-01 1.000000E+00
P 4
1 1.770000E+01 4.301800E-02
2 3.854000E+00 2.289130E-01
3 1.046000E+00 5.087280E-01
4 2.753000E-01 4.605310E-01
P 1
1 2.753000E-01 1.000000E+00
D 1
1 1.185000E+00 1.0000000
"
Une orbitale atomique centree sur l'atome A est composee d'une contraction de m Gaussiennes: $$ \chi(r) = \sum_{i=1}^m 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$.
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¶
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¶
let hf = HartreeFock.make simulation
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:
let mo_coef = MOBasis.mo_coef mo_basis
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) ; result.(2) ]
in
plot data "test_data.png"
Calcul Hartree-Fock a la main¶
Methode¶
$$ \hat{H} = \hat{T} + \hat{V}^{\text{NN}} + \hat{V}^{\text{eN}} + \hat{V}^{\text{ee}} $$On exprime les differentes quantites dans la base des orbitales atomiques pour chacun des termes:
\begin{eqnarray} T_{ij} & = & \langle \chi_i | \hat{T} | \chi_j \rangle = \iiint \chi_i(\mathbf{r}) \left( -\frac{1}{2} \Delta \chi_j(\mathbf{r}) \right) d\mathbf{r} \\ V^{NN} & = & \text{constante} \\ V^{eN}_{ij} & = & \langle \chi_i | \hat{V}^{\text{eN}} | \chi_j \rangle = \sum_A \iiint \chi_i(\mathbf{r}) \frac{-Z_A}{|\mathbf{r} - \mathbf{R}_A|} \chi_j(\mathbf{r}) d\mathbf{r} \\ V^{ee}_{ijkl} & = & \langle \chi_i \chi_j | \hat{V}^{\text{ee}} | \chi_k \chi_l \rangle = \iiint \iiint \chi_i(\mathbf{r}_1) \chi_j(\mathbf{r}_2) \frac{1}{|\mathbf{r}_1 - \mathbf{r}_2|} \chi_k(\mathbf{r}) \chi_l(\mathbf{r}) d\mathbf{r}_1 d\mathbf{r}_2 \\ \end{eqnarray}ou $\mathbf{R}_A$ est la position du noyau $A$ et $Z_A$ sa charge.
Dans la methode Hartree-Fock, l'electron 1 ne "voit" pas directement l'electron 2, mais il voit l'electron 2 comme une densite de charge.
On definit donc 2 operateurs pour chaque orbitale:
- Coulomb : $$ \hat{J}_j \chi_i(\mathbf{r}_1) = \chi_i(\mathbf{r}_1) \iiint \frac{1}{|\mathbf{r}_1 - \mathbf{r}_2|} |\chi_j(\mathbf{r}_2)|^2 d \mathbf{r}_2 $$
- Echange : $$ \hat{K}_j \chi_i(\mathbf{r}_1) = \chi_j(\mathbf{r}_1) \iiint \frac{1}{|\mathbf{r}_1 - \mathbf{r}_2|} \chi_i(\mathbf{r}_2) \chi_j(\mathbf{r}_2) d \mathbf{r}_2 $$ et on n'a plus que des operateurs a 1 electron.
ou $P$ est la matrice densite definie comme $$ P_{ij} = \sum_{k=1}^{N_{\text{occ}}} 2 C_{ik} C_{kj} $$ et $C$ est la matrice des coefficients des orbitales moleculaires exprimees dans la base des orbitales atomiques.
Une orbitale moleculaire est une combinaison lineaire d'orbitale atomique: $$ \phi_k(\mathbf{r}) = \sum_i C_{ik} \chi_i(\mathbf{r}) $$
Les orbitales moleculaires sont toutes orthogonales entre elles, et normalisees a 1.
La methode Hartree-Fock permet d'obtenir les orbitales moleculaires qui minimisent l'energie quand la fonction d'onde est exprimee comme un determinant de Slater: $$ \Psi(\mathbf{r}_1,\dots,\mathbf{r}_N) = \left| \begin{array}{ccc} \phi_1(\mathbf{r}_1) & \dots & \phi_N(\mathbf{r}_1) \\ \vdots & \ddots & \vdots \\ \phi_1(\mathbf{r}_N) & \dots & \phi_N(\mathbf{r}_N) \end{array}\right| $$
Application¶
On commence par creer une base orthogonale $X$ a partir des OA
let m_X =
AOBasis.ortho ao_basis
Pour pouvoir calculer les matrices $J$ et $K$, il faut avoir une matrice densite. On part donc avec une matrice densite d'essai. Une facon de faire est de diagonaliser l'Hamiltonien sans $V_{ee}$.
let nocc = 5 (* On a 10 electrons, donc 5 orbitales occupees. *)
let c_of_h m_H =
(* On exprime H dans la base orthonormale *)
let m_Hmo =
Util.xt_o_x m_H m_X (* H_mo = X^t H X *)
in
(* On diagonalise cet Hamiltonien *)
let m_C', _ =
Util.diagonalize_symm m_Hmo
in
(* On re-exprime les MOs dans la base des AOs (non-orthonormales) *)
gemm m_X m_C' (* C = X.C' *)
let m_C =
match Guess.make ~nocc ~guess:`Hcore ao_basis with
| Hcore m_H -> c_of_h m_H
| _ -> assert false
On construit la matrice densite
let m_P =
(* P = 2 C.C^t *)
gemm ~alpha:2. ~transb:`T ~k:nocc m_C m_C
On construit les matrices $H_c = T+V_{\text{eN}}$, et $J$ et $K$:
let m_Hc, m_J, m_K =
let f =
Fock.make_rhf ~density:m_P ao_basis
in
Fock.(core f, coulomb f, exchange f)
On construit l'operateur de Fock: $ F = H_c + J - K $
let m_F =
Mat.add m_Hc (Mat.sub m_J m_K)
On l'exprime dans la base orthonormale, on le diagonalise, et on le re-exprime ses vecteurs propres dans la base d'OA.
let m_C =
let m_C', _ =
Util.xt_o_x m_F m_X
|> Util.diagonalize_symm
in
gemm m_X m_C'
On calcule l'energie avec ces nouvelles orbitales:
let energy =
(Simulation.nuclear_repulsion simulation) +. 0.5 *.
Mat.gemm_trace m_P (Mat.add m_Hc m_F)
On itere:
let rec iteration m_C n =
let m_P =
(* P = 2 C.C^t *)
gemm ~alpha:2. ~transb:`T ~k:nocc m_C m_C
in
let m_Hc, m_J, m_K =
let f =
Fock.make_rhf ~density:m_P ao_basis
in
Fock.(core f, coulomb f, exchange f)
in
let m_F =
Mat.add m_Hc (Mat.sub m_J m_K)
in
let m_C =
let m_C', _ =
Util.xt_o_x m_F m_X
|> Util.diagonalize_symm
in
gemm m_X m_C'
in
let energy =
(Simulation.nuclear_repulsion simulation) +. 0.5 *.
Mat.gemm_trace m_P (Mat.add m_Hc m_F)
in
Printf.printf "%f\n%!" energy;
if n > 0 then
iteration m_C (n-1)
let () = iteration m_C 20
let ee_ints = AOBasis.ee_ints ao_basis;;
let m_C = MOBasis.mo_coef mo_basis;;
let n_ao = Mat.dim1 m_C ;;
let sum a =
Array.fold_left (fun accu x -> accu +. x) 0. a
;;
let i = 1;;
let j = 2;;
Array.map (fun a ->
let v =
Array.map (fun b ->
let u =
Array.map (fun e ->
let t = Array.map (fun f ->
m_C.{f,j} *. ERI.get_phys ee_ints a b e f
) (Util.array_range 1 n_ao)
in m_C.{e,i} *. sum t
) (Util.array_range 1 n_ao)
in m_C.{b,j} *. sum u
) (Util.array_range 1 n_ao)
in m_C.{a,i} *. sum v
) (Util.array_range 1 n_ao)
|> sum
$\langle \phi_i \phi_j | \phi_i \phi_j \rangle = \sum_a c_{ai} \langle \chi_a \phi_j | \phi_i \phi_j \rangle = \sum_a \sum_b c_{ai} c_{bj} \langle \chi_a \chi_b | \phi_i \phi_j \rangle$ $\sum_a \sum_b \sum_e \sum_f \left(c_{ai} \left(c_{bj} \left(c_{ei} \left(c_{fj} \langle \chi_a \chi_b | \chi_e \chi_f \rangle \right)\right)\right)\right) $
$\phi_i ^2 = \left( \sum_a c_{ai} \chi_a \right)^2$ $= \sum_a c_{ai} \sum_b c_{bi} \chi_a \chi_b$
$\phi_i ^2 - \phi_j ^2 = \left( \sum_a c_{ai} \chi_a \right)^2 - \left( \sum_a c_{aj} \chi_a \right)^2$ $= \sum_a c_{ai} \sum_b c_{bi} \chi_a \chi_b - \sum_a c_{aj} \sum_b c_{bj} \chi_a \chi_b$
$[\phi_i ^2 - \phi_j ^2 | \phi_i \phi_j ] = [\left( \sum_a c_{ai} \chi_i \right)^2 - \left( \sum_a c_{aj} \chi_j \right)^2| \phi_i \phi_j ] $ $= \sum_a c_{ai} \sum_b c_{bi} [\chi_a \chi_b| \phi_i \phi_j ] - \sum_a c_{aj} \sum_b c_{bj} [ \chi_a \chi_b| \phi_i \phi_j ] $ $= \left(\sum_a c_{ai} \sum_b c_{bi} - \sum_a c_{aj} \sum_b c_{bj} \right) [ \chi_a \chi_b| \phi_i \phi_j ] $ $= \left(\sum_a c_{ai} \sum_b c_{bi} - \sum_a c_{aj} \sum_b c_{bj} \right) \sum_e \sum_f c_{ei} c_{fj} [ \chi_a \chi_b| \chi_e \chi_f ] $
$\langle \phi_i | x | \phi_j \rangle$ $=\langle \sum_a c_{ia} \chi_a | x | \sum_b c_{jb} \chi_b\rangle$ $=\sum_a \sum_b c_{jb} c_{ia} \langle \chi_a | x | \chi_b\rangle$
let ints = MOBasis.ee_ints mo_basis in
ERI.get_phys ints 1 2 1 2;;