From ecf61058f996ad3fa31d2ecf61bad0a91d436602 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 22 Feb 2019 19:19:11 +0100 Subject: [PATCH] S2 is programmed in CI --- CI/CI.ml | 157 ++++++++++------------------------------- CI/CIMatrixElement.ml | 159 ++++++++++++++++++++++++++++++++++++++++++ CI/Determinant.ml | 11 ++- CI/Determinant.mli | 6 ++ MOBasis/MOBasis.ml | 13 ++-- Utils/Spin.ml | 43 +++++++++++- Utils/Spin.mli | 1 + run_cis.ml | 2 +- run_hartree_fock.ml | 2 +- run_integrals.ml | 2 +- 10 files changed, 265 insertions(+), 131 deletions(-) create mode 100644 CI/CIMatrixElement.ml diff --git a/CI/CI.ml b/CI/CI.ml index 13972bb..c4e76f7 100644 --- a/CI/CI.ml +++ b/CI/CI.ml @@ -1,13 +1,12 @@ open Lacaml.D -module De = Determinant -module Ex = Excitation -module Sp = Spindeterminant +module Ds = Determinant_space type t = { - det_space : Determinant_space.t ; + det_space : Ds.t ; h_matrix : Mat.t lazy_t ; + s2_matrix : Mat.t lazy_t ; eigensystem : (Mat.t * Vec.t) lazy_t; } @@ -15,6 +14,8 @@ let det_space t = t.det_space let h_matrix t = Lazy.force t.h_matrix +let s2_matrix t = Lazy.force t.s2_matrix + let eigensystem t = Lazy.force t.eigensystem let eigenvectors t = @@ -23,131 +24,49 @@ let eigenvectors t = let eigenvalues t = let (_,x) = eigensystem t in x - -let h_ij mo_basis ki kj = +let h_integrals mo_basis = let one_e_ints = MOBasis.one_e_ints mo_basis and two_e_ints = MOBasis.two_e_ints mo_basis - and degree_a, degree_b = De.degrees ki kj in - if degree_a+degree_b > 2 then - 0. - else - begin - let kia = De.alfa ki and kib = De.beta ki - and kja = De.alfa kj and kjb = De.beta kj - in - let integral = - ERI.get_phys two_e_ints - in - let anti_integral i j k l = - integral i j k l -. integral i j l k - in - let single h p same opposite = - let same_spin = - Sp.to_list same - |> List.fold_left (fun accu i -> accu +. anti_integral h i p i) 0. - and opposite_spin = - Sp.to_list opposite - |> List.fold_left (fun accu i -> accu +. integral h i p i) 0. - and one_e = - one_e_ints.{h,p} - in one_e +. same_spin +. opposite_spin - in - let diag_element () = - let mo_a = Sp.to_list kia - and mo_b = Sp.to_list kib - in - let one_e = - List.concat [mo_a ; mo_b] - |> List.fold_left (fun accu i -> accu +. one_e_ints.{i,i}) 0. - and two_e = - let two_index i j = integral i j i j - and anti_two_index i j = anti_integral i j i j - in - let rec aux_same accu = function - | [] -> accu - | i :: rest -> - let new_accu = - List.fold_left (fun accu j -> accu +. anti_two_index i j) accu rest - in - aux_same new_accu rest - in - let rec aux_opposite accu other = function - | [] -> accu - | i :: rest -> - let new_accu = - List.fold_left (fun accu j -> accu +. two_index i j) accu other - in - aux_opposite new_accu other rest - in - (aux_same 0. mo_a) +. (aux_same 0. mo_b) +. (aux_opposite 0. mo_a mo_b) - in - one_e +. two_e - in - match degree_a, degree_b with - | 1, 1 -> (* alpha-beta double *) - begin - let ha, pa, phase_a = Ex.single_of_spindet kia kja in - let hb, pb, phase_b = Ex.single_of_spindet kib kjb in - match phase_a, phase_b with - | Phase.Pos, Phase.Pos - | Phase.Neg, Phase.Neg -> +. integral ha hb pa pb - | Phase.Neg, Phase.Pos - | Phase.Pos, Phase.Neg -> -. integral ha hb pa pb - end + ( (fun i j _ -> one_e_ints.{i,j}), + (fun i j k l s s' -> + if s' = Spin.other s then + ERI.get_phys two_e_ints i j k l + else + (ERI.get_phys two_e_ints i j k l) -. + (ERI.get_phys two_e_ints i j l k) + ) ) - | 2, 0 -> (* alpha double *) - begin - let h1, p1, h2, p2, phase = Ex.double_of_spindet kia kja in - match phase with - | Phase.Pos -> +. anti_integral h1 h2 p1 p2 - | Phase.Neg -> -. anti_integral h1 h2 p1 p2 - end - | 0, 2 -> (* beta double *) - begin - let h1, p1, h2, p2, phase = Ex.double_of_spindet kib kjb in - match phase with - | Phase.Pos -> +. anti_integral h1 h2 p1 p2 - | Phase.Neg -> -. anti_integral h1 h2 p1 p2 - end - | 1, 0 -> (* alpha single *) - begin - let h, p, phase = Ex.single_of_spindet kia kja in - match phase with - | Phase.Pos -> +. single h p kia kib - | Phase.Neg -> -. single h p kia kib - end +let h_ij mo_basis ki kj = + let integrals = + List.map (fun f -> f mo_basis) + [ h_integrals ] + in + CIMatrixElement.make integrals ki kj + |> List.hd - | 0, 1 -> (* beta single *) - begin - let h, p, phase = Ex.single_of_spindet kib kjb in - match phase with - | Phase.Pos -> +. single h p kib kia - | Phase.Neg -> -. single h p kib kia - end - - | 0, 0 -> (* diagonal element *) - diag_element () - - | _ -> assert false - end let make det_space = - let ndet = Determinant_space.size det_space in - let det = Determinant_space.determinants det_space in - let mo_basis = Determinant_space.mo_basis det_space in + let ndet = Ds.size det_space in + let det = Ds.determinants det_space in + let mo_basis = Ds.mo_basis det_space in let h_matrix = lazy ( - Array.init ndet (fun i -> - let ki = det.(i) in - Array.init ndet (fun j -> - let kj = det.(j) in - h_ij mo_basis ki kj - ) - ) - |> Mat.of_array + Array.init ndet (fun i -> let ki = det.(i) in + Array.init ndet (fun j -> let kj = det.(j) in + h_ij mo_basis ki kj + )) + |> Mat.of_array + ) + in + let s2_matrix = lazy ( + Array.init ndet (fun i -> let ki = det.(i) in + Array.init ndet (fun j -> let kj = det.(j) in + CIMatrixElement.make_s2 ki kj + )) + |> Mat.of_array ) in let eigensystem = lazy ( @@ -155,5 +74,5 @@ let make det_space = |> Util.diagonalize_symm ) in - { det_space ; h_matrix ; eigensystem } + { det_space ; h_matrix ; s2_matrix ; eigensystem } diff --git a/CI/CIMatrixElement.ml b/CI/CIMatrixElement.ml new file mode 100644 index 0000000..217dcd5 --- /dev/null +++ b/CI/CIMatrixElement.ml @@ -0,0 +1,159 @@ +open Lacaml.D + +module De = Determinant +module Ex = Excitation +module Sp = Spindeterminant + +type t = float list + + +let non_zero integrals degree_a degree_b ki kj = + let kia = De.alfa ki and kib = De.beta ki + and kja = De.alfa kj and kjb = De.beta kj + in + let single h p spin same opposite = + let same_spin_mo_list = + Sp.to_list same + and opposite_spin_mo_list = + Sp.to_list opposite + in + fun one_e two_e -> + let same_spin = + List.fold_left (fun accu i -> accu +. two_e h i p i spin spin) 0. same_spin_mo_list + and opposite_spin = + List.fold_left (fun accu i -> accu +. two_e h i p i spin (Spin.other spin) ) 0. opposite_spin_mo_list + in (one_e h p spin) +. same_spin +. opposite_spin + in + let diag_element = + let mo_a = Sp.to_list kia + and mo_b = Sp.to_list kib + in + fun one_e two_e -> + let one = + (List.fold_left (fun accu i -> accu +. one_e i i Spin.Alfa) 0. mo_a) + +. + (List.fold_left (fun accu i -> accu +. one_e i i Spin.Beta) 0. mo_b) + in + let two = + let rec aux_same spin accu = function + | [] -> accu + | i :: rest -> + let new_accu = + List.fold_left (fun accu j -> accu +. two_e i j i j spin spin) accu rest + in + aux_same spin new_accu rest + in + let rec aux_opposite accu other = function + | [] -> accu + | i :: rest -> + let new_accu = + List.fold_left (fun accu j -> accu +. two_e i j i j Spin.Alfa Spin.Beta) accu other + in + aux_opposite new_accu other rest + in + (aux_same Spin.Alfa 0. mo_a) +. (aux_same Spin.Beta 0. mo_b) +. + (aux_opposite 0. mo_a mo_b) + in + one +. two + in + let result = + match degree_a, degree_b with + | 1, 1 -> (* alpha-beta double *) + begin + let ha, pa, phase_a = Ex.single_of_spindet kia kja in + let hb, pb, phase_b = Ex.single_of_spindet kib kjb in + match phase_a, phase_b with + | Phase.Pos, Phase.Pos + | Phase.Neg, Phase.Neg -> fun _ two_e -> two_e ha hb pa pb Spin.Alfa Spin.Beta + | Phase.Neg, Phase.Pos + | Phase.Pos, Phase.Neg -> fun _ two_e -> -. two_e ha hb pa pb Spin.Alfa Spin.Beta + end + + | 2, 0 -> (* alpha double *) + begin + let h1, p1, h2, p2, phase = Ex.double_of_spindet kia kja in + match phase with + | Phase.Pos -> fun _ two_e -> two_e h1 h2 p1 p2 Spin.Alfa Spin.Alfa + | Phase.Neg -> fun _ two_e -> -. two_e h1 h2 p1 p2 Spin.Alfa Spin.Alfa + end + + | 0, 2 -> (* beta double *) + begin + let h1, p1, h2, p2, phase = Ex.double_of_spindet kib kjb in + match phase with + | Phase.Pos -> fun _ two_e -> two_e h1 h2 p1 p2 Spin.Beta Spin.Beta + | Phase.Neg -> fun _ two_e -> -. two_e h1 h2 p1 p2 Spin.Beta Spin.Beta + end + + | 1, 0 -> (* alpha single *) + begin + let h, p, phase = Ex.single_of_spindet kia kja in + match phase with + | Phase.Pos -> fun one_e two_e -> single h p Spin.Alfa kia kib one_e two_e + | Phase.Neg -> fun one_e two_e -> -. single h p Spin.Alfa kia kib one_e two_e + end + + | 0, 1 -> (* beta single *) + begin + let h, p, phase = Ex.single_of_spindet kib kjb in + match phase with + | Phase.Pos -> fun one_e two_e -> single h p Spin.Beta kib kia one_e two_e + | Phase.Neg -> fun one_e two_e -> -. single h p Spin.Beta kib kia one_e two_e + end + + | 0, 0 -> (* diagonal element *) + diag_element + + | _ -> assert false + in + List.map (fun (one_e, two_e) -> result one_e two_e) integrals + + +let make integrals ki kj = + let degree_a, degree_b = + De.degrees ki kj + in + if degree_a+degree_b > 2 then + List.map (fun _ -> 0.) integrals + else + non_zero integrals degree_a degree_b ki kj + + + + +let make_s2 ki kj = + let degree_a = De.degree_alfa ki kj in + let kia = De.alfa ki in + let kja = De.alfa kj in + if degree_a > 1 then 0. + else + let degree_b = De.degree_beta ki kj in + let kib = De.beta ki in + let kjb = De.beta kj in + match degree_a, degree_b with + | 1, 1 -> (* alpha-beta double *) + let ha, pa, phase_a = Ex.single_of_spindet kia kja in + let hb, pb, phase_b = Ex.single_of_spindet kib kjb in + if ha = pb && hb = pa then + begin + match phase_a, phase_b with + | Phase.Pos, Phase.Pos + | Phase.Neg, Phase.Neg -> -1. + | Phase.Neg, Phase.Pos + | Phase.Pos, Phase.Neg -> 1. + end + else 0. + | 0, 0 -> + let ba = Sp.bitstring kia and bb = Sp.bitstring kib in + let tmp = Z.(logxor ba bb) in + let popcount x = + if x = Z.zero then 0 else + Z.popcount x + in + let n_a = Z.(logand ba tmp) |> popcount in + let n_b = Z.(logand bb tmp) |> popcount in + let s_z = 0.5 *. float_of_int (n_a - n_b) in + float_of_int n_a +. s_z *. (s_z -. 1.) + | _ -> 0. + + diff --git a/CI/Determinant.ml b/CI/Determinant.ml index ae5b918..8afbf9e 100644 --- a/CI/Determinant.ml +++ b/CI/Determinant.ml @@ -41,13 +41,18 @@ let negate_phase t = { t with alfa = Spindeterminant.negate_phase t.alfa } +let degree_alfa t t' = + Spindeterminant.degree t.alfa t'.alfa + +let degree_beta t t' = + Spindeterminant.degree t.beta t'.beta + let degrees t t' = - Spindeterminant.degree t.alfa t'.alfa, - Spindeterminant.degree t.beta t'.beta + degree_alfa t t', degree_beta t t' let degree t t' = - let a, b = degrees t t' in a+b + (degree_alfa t t') + (degree_beta t t') let of_lists a b = diff --git a/CI/Determinant.mli b/CI/Determinant.mli index a7b2d6a..839f846 100644 --- a/CI/Determinant.mli +++ b/CI/Determinant.mli @@ -44,6 +44,12 @@ val single_excitation : Spin.t -> hole -> particle -> t -> t val double_excitation : Spin.t -> hole -> particle -> Spin.t -> hole -> particle -> t -> t (** Double excitation operator {% $T_{hh'}^{pp'} = a^\dagger_p a^\dagger_{p'} a_{h'} a_h$ %}. *) +val degree_alfa : t -> t -> int +(** Returns degree of excitation between two determinants in the {% $\alpha$ %} spin. *) + +val degree_beta : t -> t -> int +(** Returns degree of excitation between two determinants in the {% $\beta$ %} spin. *) + val degrees : t -> t -> int*int (** Returns degrees of excitation between two determinants in {% $\alpha$ %} and {% $\beta$ %} spins as a pair. *) diff --git a/MOBasis/MOBasis.ml b/MOBasis/MOBasis.ml index f437743..7222b4c 100644 --- a/MOBasis/MOBasis.ml +++ b/MOBasis/MOBasis.ml @@ -21,6 +21,7 @@ type t = eN_ints : NucInt.t lazy_t; (* Electron-nucleus potential integrals *) ee_ints : ERI.t lazy_t; (* Electron-electron potential integrals *) kin_ints : KinInt.t lazy_t; (* Kinetic energy integrals *) + one_e_ints : Mat.t lazy_t; (* Kinetic energy integrals *) } @@ -36,9 +37,7 @@ let eN_ints t = Lazy.force t.eN_ints let ee_ints t = Lazy.force t.ee_ints let kin_ints t = Lazy.force t.kin_ints let two_e_ints t = Lazy.force t.ee_ints -let one_e_ints t = - Mat.add (NucInt.matrix @@ Lazy.force t.eN_ints) - (KinInt.matrix @@ Lazy.force t.kin_ints) +let one_e_ints t = Lazy.force t.one_e_ints let mo_matrix_of_ao_matrix ~mo_coef ao_matrix = xt_o_x ~x:mo_coef ~o:ao_matrix @@ -150,8 +149,12 @@ let make ~simulation ~mo_type ~mo_occupation ~mo_coef () = |> four_index_transform ~mo_coef ) in - { simulation ; mo_type ; mo_occupation ; mo_coef ; - eN_ints ; ee_ints ; kin_ints } + let one_e_ints = lazy ( + Mat.add (NucInt.matrix @@ Lazy.force eN_ints) + (KinInt.matrix @@ Lazy.force kin_ints) ) + in + { simulation ; mo_type ; mo_occupation ; mo_coef ; + eN_ints ; ee_ints ; kin_ints ; one_e_ints } let of_rhf hf = diff --git a/Utils/Spin.ml b/Utils/Spin.ml index dffb792..63ffa8b 100644 --- a/Utils/Spin.ml +++ b/Utils/Spin.ml @@ -1,6 +1,47 @@ (** Electron spin *) -type t = Alfa | Beta +type t = (* m_s *) +| Alfa (* {% $m_s = +1/2$ %} *) +| Beta (* {% $m_s = -1/2$ %} *) + +let other = function +| Alfa -> Beta +| Beta -> Alfa + +(* +let half = 1. /. 2. +(** {% $\alpha(m_s)$ %} *) +let alfa = function +| n, Alfa -> n +| _, Beta -> 0. + +(** {% $\beta(m_s)$ %} *) +let beta = function +| _, Alfa -> 0. +| n, Beta -> n + + +(** {% $S_z(m_s)$ %} *) +let s_z = function +| n, Alfa -> half *. n, Alfa +| n, Beta -> -. half *. n, Beta + +let s_plus = function +| n, Beta -> n , Alfa +| _, Alfa -> 0., Alfa + +let s_minus = function +| n, Alfa -> n , Beta +| _, Beta -> 0., Beta + +let ( ++ ) (n, t) (m,t) = + (m. +n.) + +let s2 s = + s_minus @@ s_plus s +. + s_z s +. + s_z @@ s_z s + *) diff --git a/Utils/Spin.mli b/Utils/Spin.mli index c09b56a..89e3a82 100644 --- a/Utils/Spin.mli +++ b/Utils/Spin.mli @@ -8,5 +8,6 @@ letters as [Beta], so the alignment of the code is nicer. type t = Alfa | Beta +val other : t -> t diff --git a/run_cis.ml b/run_cis.ml index 11e6271..d4bef7a 100644 --- a/run_cis.ml +++ b/run_cis.ml @@ -29,7 +29,7 @@ let run ~out = | Some x -> x and nuclei_file = match !nuclei_file with - | None -> raise (Invalid_argument "Basis set file should be specified with -c") + | None -> raise (Invalid_argument "Coordinate file should be specified with -c") | Some x -> x and charge = !charge and multiplicity = !multiplicity diff --git a/run_hartree_fock.ml b/run_hartree_fock.ml index ec8aee3..d378bc7 100644 --- a/run_hartree_fock.ml +++ b/run_hartree_fock.ml @@ -29,7 +29,7 @@ let run ~out = | Some x -> x and nuclei_file = match !nuclei_file with - | None -> raise (Invalid_argument "Basis set file should be specified with -c") + | None -> raise (Invalid_argument "Coordinate file should be specified with -c") | Some x -> x and charge = !charge and multiplicity = !multiplicity diff --git a/run_integrals.ml b/run_integrals.ml index 20098db..738e19f 100644 --- a/run_integrals.ml +++ b/run_integrals.ml @@ -23,7 +23,7 @@ let run ~out = | Some x -> x and nuclei_file = match !nuclei_file with - | None -> raise (Invalid_argument "Basis set file should be specified with -c") + | None -> raise (Invalid_argument "Coordinate file should be specified with -c") | Some x -> x in