diff --git a/CI/CI.ml b/CI/CI.ml index 77d15a0..7752661 100644 --- a/CI/CI.ml +++ b/CI/CI.ml @@ -5,16 +5,16 @@ module Ds = Determinant_space type t = { det_space : Ds.t ; - h_matrix : Mat.t lazy_t ; - s2_matrix : Mat.t lazy_t ; + m_H : Mat.t lazy_t ; + m_S2 : Mat.t lazy_t ; eigensystem : (Mat.t * Vec.t) lazy_t; } let det_space t = t.det_space -let h_matrix t = Lazy.force t.h_matrix +let m_H t = Lazy.force t.m_H -let s2_matrix t = Lazy.force t.s2_matrix +let m_S2 t = Lazy.force t.m_S2 let eigensystem t = Lazy.force t.eigensystem @@ -54,7 +54,7 @@ let make det_space = let det = Ds.determinants det_space in let mo_basis = Ds.mo_basis det_space in (* - let h_matrix = lazy ( + let m_H = lazy ( Util.list_range 0 (ndet-1) |> List.map (fun i -> let ki = det.(i) in Array.init ndet (fun j -> let kj = det.(j) in @@ -66,7 +66,7 @@ let make det_space = in *) (* - let h_matrix = lazy ( + let m_H = lazy ( let ntasks = int_of_float @@ sqrt @@ float_of_int ndet in List.init ntasks (fun i -> let m = @@ -93,7 +93,7 @@ let make det_space = |> Mat.of_col_vecs_list ) in *) - let h_matrix = lazy ( + let m_H = lazy ( let h = if Parallel.master then Array.make_matrix ndet ndet 0. @@ -120,7 +120,7 @@ let make det_space = List.iter (fun (i,j,x) -> h.{i+1,j+1} <- x) l); Parallel.broadcast (lazy h) ) in - let s2_matrix = lazy ( + let m_S2 = 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 @@ -129,10 +129,10 @@ let make det_space = ) in let eigensystem = lazy ( - let h_matrix = Lazy.force h_matrix in + let m_H = Lazy.force m_H in Parallel.broadcast @@ - lazy (Util.diagonalize_symm h_matrix) + lazy (Util.diagonalize_symm m_H) ) in - { det_space ; h_matrix ; s2_matrix ; eigensystem } + { det_space ; m_H ; m_S2 ; eigensystem } diff --git a/MOBasis/MOBasis.ml b/MOBasis/MOBasis.ml index 0b1a64b..b976945 100644 --- a/MOBasis/MOBasis.ml +++ b/MOBasis/MOBasis.ml @@ -45,8 +45,7 @@ let mo_matrix_of_ao_matrix ~mo_coef ao_matrix = let ao_matrix_of_mo_matrix ~mo_coef ~ao_overlap mo_matrix = let sc = gemm ao_overlap mo_coef in - gemm sc @@ - gemm ~transb:`T mo_matrix sc + x_o_xt ~x:sc ~o:mo_matrix let four_index_transform ~mo_coef eri_ao = diff --git a/SCF/Fock.ml b/SCF/Fock.ml index a6604cb..194d18f 100644 --- a/SCF/Fock.ml +++ b/SCF/Fock.ml @@ -12,11 +12,16 @@ type t = exchange : Mat.t ; } + +let fock t = t.fock +let core t = t.core +let coulomb t = t.coulomb +let exchange t = t.exchange + + module Ao = AOBasis - - -let make ~density ?(threshold=Constants.epsilon) ao_basis = +let make_rhf ~density ?(threshold=Constants.epsilon) ao_basis = let m_P = density and m_T = Ao.kin_ints ao_basis |> KinInt.matrix and m_V = Ao.eN_ints ao_basis |> NucInt.matrix @@ -41,47 +46,47 @@ let make ~density ?(threshold=Constants.epsilon) ao_basis = match (abs_float pJ > threshold , abs_float pK > threshold, nu < sigma) with | (false, false, _) -> () | (true , true , true) -> - begin - for mu = 1 to nu do - let integral = - ERI.get_phys m_G mu lambda nu sigma - in - if (integral <> 0.) then begin - m_Jnu.(mu-1) <- m_Jnu.(mu-1) +. pJ *. integral; - m_Ksigma.(mu-1) <- m_Ksigma.(mu-1) -. pK *. integral - end - done; - for mu = nu+1 to sigma do - m_Ksigma.(mu-1) <- m_Ksigma.(mu-1) -. pK *. - ERI.get_phys m_G mu lambda nu sigma - done - end - | (true , true , false) -> - begin - for mu = 1 to sigma do - let integral = - ERI.get_phys m_G mu lambda nu sigma - in - if (integral <> 0.) then begin - m_Jnu.(mu-1) <- m_Jnu.(mu-1) +. pJ *. integral; - m_Ksigma.(mu-1) <- m_Ksigma.(mu-1) -. pK *. integral - end - done; - for mu = sigma+1 to nu do - m_Jnu.(mu-1) <- m_Jnu.(mu-1) +. pJ *. - ERI.get_phys m_G mu lambda nu sigma - done - end - | (false, true , _) -> - for mu = 1 to sigma do - m_Ksigma.(mu-1) <- m_Ksigma.(mu-1) -. pK *. - ERI.get_phys m_G mu lambda nu sigma - done - | (true , false, _) -> + begin for mu = 1 to nu do - m_Jnu.(mu-1) <- m_Jnu.(mu-1) +. pJ *. - ERI.get_phys m_G mu lambda nu sigma + let integral = + ERI.get_phys m_G mu lambda nu sigma + in + if (integral <> 0.) then begin + m_Jnu.(mu-1) <- m_Jnu.(mu-1) +. pJ *. integral; + m_Ksigma.(mu-1) <- m_Ksigma.(mu-1) -. pK *. integral + end + done; + for mu = nu+1 to sigma do + m_Ksigma.(mu-1) <- m_Ksigma.(mu-1) -. pK *. + ERI.get_phys m_G mu lambda nu sigma done + end + | (true , true , false) -> + begin + for mu = 1 to sigma do + let integral = + ERI.get_phys m_G mu lambda nu sigma + in + if (integral <> 0.) then begin + m_Jnu.(mu-1) <- m_Jnu.(mu-1) +. pJ *. integral; + m_Ksigma.(mu-1) <- m_Ksigma.(mu-1) -. pK *. integral + end + done; + for mu = sigma+1 to nu do + m_Jnu.(mu-1) <- + m_Jnu.(mu-1) +. pJ *. ERI.get_phys m_G mu lambda nu sigma + done + end + | (false, true , _) -> + for mu = 1 to sigma do + m_Ksigma.(mu-1) <- + m_Ksigma.(mu-1) -. pK *. ERI.get_phys m_G mu lambda nu sigma + done + | (true , false, _) -> + for mu = 1 to nu do + m_Jnu.(mu-1) <- + m_Jnu.(mu-1) +. pJ *. ERI.get_phys m_G mu lambda nu sigma + done done done; for mu = 1 to sigma-1 do @@ -103,7 +108,100 @@ let make ~density ?(threshold=Constants.epsilon) ao_basis = +let make_uhf ~density_same ~density_other ?(threshold=Constants.epsilon) ao_basis = + let m_P_a = density_same + and m_P_b = density_other + and m_T = Ao.kin_ints ao_basis |> KinInt.matrix + and m_V = Ao.eN_ints ao_basis |> NucInt.matrix + and m_G = Ao.ee_ints ao_basis + in + let nBas = Mat.dim1 m_T + in + + let m_Hc = Mat.add m_T m_V + and m_J = Array.make_matrix nBas nBas 0. + and m_K = Array.make_matrix nBas nBas 0. + in + + for sigma = 1 to nBas do + let m_Ksigma = m_K.(sigma-1) in + for nu = 1 to nBas do + let m_Jnu = m_J.(nu-1) in + for lambda = 1 to nBas do + let pJ = m_P_a.{lambda,sigma} +. m_P_b.{lambda,sigma} + and pK = m_P_a.{lambda,nu} + in + match (abs_float pJ > threshold , abs_float pK > threshold, nu < sigma) with + | (false, false, _) -> () + | (true , true , true) -> + begin + for mu = 1 to nu do + let integral = + ERI.get_phys m_G mu lambda nu sigma + in + if (integral <> 0.) then begin + m_Jnu.(mu-1) <- m_Jnu.(mu-1) +. pJ *. integral; + m_Ksigma.(mu-1) <- m_Ksigma.(mu-1) -. pK *. integral + end + done; + for mu = nu+1 to sigma do + m_Ksigma.(mu-1) <- m_Ksigma.(mu-1) -. pK *. + ERI.get_phys m_G mu lambda nu sigma + done + end + | (true , true , false) -> + begin + for mu = 1 to sigma do + let integral = + ERI.get_phys m_G mu lambda nu sigma + in + if (integral <> 0.) then begin + m_Jnu.(mu-1) <- m_Jnu.(mu-1) +. pJ *. integral; + m_Ksigma.(mu-1) <- m_Ksigma.(mu-1) -. pK *. integral + end + done; + for mu = sigma+1 to nu do + m_Jnu.(mu-1) <- + m_Jnu.(mu-1) +. pJ *. ERI.get_phys m_G mu lambda nu sigma + done + end + | (false, true , _) -> + for mu = 1 to sigma do + m_Ksigma.(mu-1) <- + m_Ksigma.(mu-1) -. pK *. ERI.get_phys m_G mu lambda nu sigma + done + | (true , false, _) -> + for mu = 1 to nu do + m_Jnu.(mu-1) <- + m_Jnu.(mu-1) +. pJ *. ERI.get_phys m_G mu lambda nu sigma + done + done + done; + for mu = 1 to sigma-1 do + m_K.(mu-1).(sigma-1) <- m_Ksigma.(mu-1); + done + done; + for nu = 1 to nBas do + let m_Jnu = m_J.(nu-1) in + for mu = 1 to nu-1 do + m_J.(mu-1).(nu-1) <- m_Jnu.(mu-1) + done + done; + + let m_J = Mat.of_array m_J + and m_K = Mat.of_array m_K + in + { fock = Mat.add m_Hc (Mat.add m_J m_K) ; + core = m_Hc ; coulomb = m_J ; exchange = m_K } + + + + + + + let op ~f f1 f2 = + assert (f1.core = f2.core); let m_Hc = f1.core and m_J = f f1.coulomb f2.coulomb and m_K = f f1.exchange f2.exchange @@ -117,8 +215,24 @@ let op ~f f1 f2 = let add = op ~f:(fun a b -> Mat.add a b) + let sub = op ~f:(fun a b -> Mat.sub a b) +let scale alpha f1 = + let m_Hc = f1.core + and m_J = lacpy f1.coulomb + and m_K = lacpy f1.exchange + in + Mat.scal alpha m_J; + Mat.scal alpha m_K; + { + fock = Mat.add m_Hc (Mat.add m_J m_K); + core = m_Hc; + coulomb = m_J; + exchange = m_K; + } + + let pp_fock ppf a = Format.fprintf ppf "@[<2>"; diff --git a/SCF/Fock.mli b/SCF/Fock.mli new file mode 100644 index 0000000..21da4e8 --- /dev/null +++ b/SCF/Fock.mli @@ -0,0 +1,50 @@ +(** Type for the Fock operator in AO basis. *) + +open Lacaml.D + + +type t + +(** {1 Accessors} *) + +val fock : t -> Mat.t +(** Fock matrix in AO basis *) + +val core : t -> Mat.t +(** Core Hamiltonian : {% $\langle i | \hat{h} | j \rangle$ %} *) + +val coulomb : t -> Mat.t +(** Coulomb matrix : {% $\langle i | J | j \rangle$ %} *) + +val exchange : t -> Mat.t +(** Exchange matrix : {% $\langle i | K | j \rangle$ %} *) + + +(** {1 Creators} *) + +val make_rhf : density:Mat.t -> ?threshold:float -> AOBasis.t -> t +(** Create a Fock operator in the RHF formalism. Expected density is + {% $2 \mathbf{C\, C}^\dagger$ %}. [threshold] is a threshold on the + integrals. *) + +val make_uhf : density_same: Mat.t -> density_other:Mat.t -> ?threshold:float -> + AOBasis.t -> t +(** Create a Fock operator in the UHF formalism. Expected density is + {% $\mathbf{C\, C}^\dagger$ %}. When building the {% $\alpha$ %} Fock + operator, [density_same] is the {% $\alpha$ %} density and [density_other] + is the {% $\beta$ %} density. [threshold] is a threshold on the integrals. *) + + +(** {1 Operations} *) + +val add : t -> t -> t +(** Add two Fock operators sharing the same core Hamiltonian. *) + +val sub : t -> t -> t +(** Subtract two Fock operators sharing the same core Hamiltonian. *) + + +(** {1 Printers} *) + +val pp_fock : Format.formatter -> t -> unit + diff --git a/SCF/Guess.ml b/SCF/Guess.ml index 5878dcc..e6d81fa 100644 --- a/SCF/Guess.ml +++ b/SCF/Guess.ml @@ -34,8 +34,8 @@ let huckel_guess ao_basis = | 0 -> invalid_arg "Huckel guess needs a non-zero number of occupied MOs." | nocc -> let density = gemm ~alpha:2. ~transb:`T ~k:nocc m_X m_X in - let fock = Fock.make ~density ao_basis in - let m_F = fock.Fock.fock in + let f = Fock.make_rhf ~density ao_basis in + let m_F = Fock.fock f in for j=1 to ao_num do for i=1 to ao_num do if (i <> j) then diff --git a/SCF/HartreeFock.ml b/SCF/HartreeFock.ml index 5b8305f..9eec1fc 100644 --- a/SCF/HartreeFock.ml +++ b/SCF/HartreeFock.ml @@ -4,7 +4,7 @@ let make ?guess simulation = if Electrons.multiplicity @@ Simulation.electrons simulation = 1 then RHF.make ?guess simulation else - invalid_arg "UHF or ROHF not implemented" + ROHF.make ?guess simulation let to_string = HartreeFock_type.to_string diff --git a/SCF/RHF.ml b/SCF/RHF.ml index 0582ad8..3545a9a 100644 --- a/SCF/RHF.ml +++ b/SCF/RHF.ml @@ -65,14 +65,14 @@ let make ?guess:(guess=`Huckel) ?max_scf:(max_scf=64) ?level_shift:(level_shift= match fock_prev, threshold > 100. *. threshold_SCF with | Some fock_prev, true -> let threshold = 1.e-8 in - Fock.make ~density:(Mat.sub m_P m_P_prev) ~threshold ao_basis + Fock.make_rhf ~density:(Mat.sub m_P m_P_prev) ~threshold ao_basis |> Fock.add fock_prev - | _ -> Fock.make ~density:m_P ao_basis + | _ -> Fock.make_rhf ~density:m_P ao_basis in let m_F, m_Hc, m_J, m_K = let x = fock in - x.Fock.fock, x.Fock.core, x.Fock.coulomb, x.Fock.exchange + Fock.(fock x, core x, coulomb x, exchange x) in (* Add level shift in AO basis *) diff --git a/SCF/ROHF.ml b/SCF/ROHF.ml new file mode 100644 index 0000000..c4aafdd --- /dev/null +++ b/SCF/ROHF.ml @@ -0,0 +1,290 @@ +open Util +open Constants +open Lacaml.D + +module Si = Simulation +module El = Electrons +module Ao = AOBasis +module Ov = Overlap + +let make ?guess:(guess=`Huckel) ?max_scf:(max_scf=64) ?level_shift:(level_shift=0.1) + ?threshold_SCF:(threshold_SCF=1.e-5) simulation = + + (* Number of occupied MOs *) + let n_alfa = + El.n_alfa @@ Si.electrons simulation + in + + let nocc = n_alfa in + + let n_beta = + El.n_beta @@ Si.electrons simulation + in + + let nuclear_repulsion = + Si.nuclear_repulsion simulation + in + + let ao_basis = + Si.ao_basis simulation + in + + (* Initial guess *) + let guess = + Guess.make ~nocc ~guess ao_basis + in + + (* Orthogonalization matrix *) + let m_X = + Ao.ortho ao_basis + in + + + (* Overlap matrix *) + let m_S = + Ao.overlap ao_basis + |> Ov.matrix + in + + let m_T = Ao.kin_ints ao_basis |> KinInt.matrix + and m_V = Ao.eN_ints ao_basis |> NucInt.matrix + in + + (* Level shift in MO basis *) + let m_LSmo = + Array.init (Mat.dim2 m_X) (fun i -> + if i > n_alfa then level_shift else 0.) + |> Vec.of_array + |> Mat.of_diag + in + + (* SCF iterations *) + let rec loop nSCF iterations energy_prev m_C m_P_a_prev m_P_b_prev fock_a_prev fock_b_prev threshold diis = + + + (* Density matrix *) + let m_P_a = + gemm ~alpha:1. ~transb:`T ~k:n_alfa m_C m_C + in + + let m_P_b = + gemm ~alpha:1. ~transb:`T ~k:n_beta m_C m_C + in + + let m_P = + Mat.add m_P_a m_P_b + in + + (* Fock matrix in AO basis *) + let fock_a = + match fock_a_prev, threshold > 100. *. threshold_SCF with + | Some fock_a_prev, true -> + let threshold = 1.e-8 in + Fock.make_uhf ~density_same:(Mat.sub m_P_a m_P_a_prev) ~density_other:(Mat.sub m_P_b m_P_b_prev) ~threshold ao_basis + |> Fock.add fock_a_prev + | _ -> Fock.make_uhf ~density_same:m_P_a ~density_other:m_P_b ao_basis + in + + let fock_b = + match fock_b_prev, threshold > 100. *. threshold_SCF with + | Some fock_b_prev, true -> + let threshold = 1.e-8 in + Fock.make_uhf ~density_same:(Mat.sub m_P_b m_P_b_prev) ~density_other:(Mat.sub m_P_a m_P_a_prev) ~threshold ao_basis + |> Fock.add fock_b_prev + | _ -> Fock.make_uhf ~density_same:m_P_b ~density_other:m_P_a ao_basis + in + + let m_F_a, m_Hc_a, m_J_a, m_K_a = + let x = fock_a in + Fock.(fock x, core x, coulomb x, exchange x) + in + + let m_F_b, m_Hc_b, m_J_b, m_K_b = + let x = fock_b in + Fock.(fock x, core x, coulomb x, exchange x) + in + + + let m_F_mo_a = + xt_o_x ~o:m_F_a ~x:m_C + in + + let m_F_mo_b = + xt_o_x ~o:m_F_b ~x:m_C + in + + let m_F_mo = + let space k = + if k <= n_beta then + `Core + else if k <= n_alfa then + `Active + else + `Virtual + in + Array.init (Mat.dim2 m_F_mo_a) (fun i -> + let i = i+1 in + Array.init (Mat.dim1 m_F_mo_a) (fun j -> + let j = j+1 in + match (space i), (space j) with + | `Core , `Core -> + 0.5 *. (m_F_mo_a.{i,j} +. m_F_mo_b.{i,j}) -. + (m_F_mo_b.{i,j} -. m_F_mo_a.{i,j}) + + | `Active , `Core + | `Core , `Active -> + (* + 0.5 *. (m_F_mo_a.{i,j} +. m_F_mo_b.{i,j}) +. + 0.5 *. (m_F_mo_b.{i,j} -. m_F_mo_a.{i,j}) + *) + m_F_mo_b.{i,j} + + | `Core , `Virtual + | `Virtual , `Core + | `Active , `Active -> + 0.5 *. (m_F_mo_a.{i,j} +. m_F_mo_b.{i,j}) + + | `Virtual , `Active + | `Active , `Virtual -> + (* + 0.5 *. (m_F_mo_a.{i,j} +. m_F_mo_b.{i,j}) -. + 0.5 *. (m_F_mo_b.{i,j} -. m_F_mo_a.{i,j}) + *) + m_F_mo_a.{i,j} + + | `Virtual , `Virtual -> + 0.5 *. (m_F_mo_a.{i,j} +. m_F_mo_b.{i,j}) +. + (m_F_mo_b.{i,j} -. m_F_mo_a.{i,j}) + ) ) + |> Mat.of_array + in + + let m_SC = + gemm m_S m_C + in + + let m_F = + x_o_xt ~x:m_SC ~o:m_F_mo + in + + + (* Add level shift in AO basis *) + let m_F = + x_o_xt ~x:m_SC ~o:m_LSmo + |> Mat.add m_F + in + + (* Fock matrix in orthogonal basis *) + let m_F_ortho = + xt_o_x m_F m_X + in + + let error_fock = + let fps = + gemm m_F (gemm m_P m_S) + and spf = + gemm m_S (gemm m_P m_F) + in + xt_o_x (Mat.sub fps spf) m_X + in + + let diis = + DIIS.append ~p:(Mat.as_vec m_F_ortho) ~e:(Mat.as_vec error_fock) diis + in + + let m_F_diis = + let x = + Bigarray.genarray_of_array1 (DIIS.next diis) + in + Bigarray.reshape_2 x (Mat.dim1 m_F_ortho) (Mat.dim2 m_F_ortho) + in + + + (* MOs in orthogonal MO basis *) + let m_C', eigenvalues = + diagonalize_symm m_F_diis + in + + (* MOs in AO basis *) + let m_C = + gemm m_X m_C' + in + + (* Hartree-Fock energy *) + let energy = + nuclear_repulsion +. 0.5 *. ( Mat.gemm_trace m_P_a (Mat.add m_Hc_a m_F_a) +. + Mat.gemm_trace m_P_b (Mat.add m_Hc_b m_F_b) ) + in + + (* Convergence criterion *) + let error = + error_fock + |> Mat.as_vec + |> amax + |> abs_float + in + + let converged = + nSCF = max_scf || error < threshold_SCF + in + + let gap = + if nocc < Vec.dim eigenvalues then + eigenvalues.{nocc+1} -. eigenvalues.{nocc} + else 0. + in + + let () = + match energy_prev with + | Some energy_prev -> + Printf.eprintf "%3d %16.10f %16.10f %11.4e %10.4f\n%!" nSCF energy (energy -. energy_prev) error gap + | None -> + Printf.eprintf "%3d %16.10f %16s %11.4e %10.4f\n%!" nSCF energy "" error gap + in + + if not converged then + loop (nSCF+1) ( (energy, error, gap) :: iterations) (Some energy) m_C m_P_a m_P_b (Some fock_a) (Some fock_b) error diis + else + let iterations = + List.rev ( (energy, error, gap) :: iterations ) + |> Array.of_list + in + HartreeFock_type.(ROHF + { + simulation; + nocc; + guess ; + eigenvectors = m_C ; + eigenvalues ; + energy ; + nuclear_repulsion; + iterations ; + kin_energy = Mat.gemm_trace m_P m_T; + eN_energy = Mat.gemm_trace m_P m_V; + coulomb_energy = 0.5 *. (Mat.gemm_trace m_P_a m_J_a) +. + 0.5 *. (Mat.gemm_trace m_P_b m_J_b); + exchange_energy = 0.5 *. (Mat.gemm_trace m_P_a m_K_a) +. + 0.5 *. (Mat.gemm_trace m_P_b m_K_b); + }) + in + + + (* Guess coefficients *) + let m_H = + match guess with + | Guess.Hcore m_H -> m_H + | Guess.Huckel 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 + + let diis = DIIS.make () in + loop 1 [] None m_C m_C m_C None None threshold_SCF diis + diff --git a/Utils/Util.ml b/Utils/Util.ml index 52677ad..02f18a9 100644 --- a/Utils/Util.ml +++ b/Utils/Util.ml @@ -227,6 +227,10 @@ let xt_o_x ~o ~x = gemm o x |> gemm ~transa:`T x +let x_o_xt ~o ~x = + gemm o x ~transb:`T + |> gemm x + let canonical_ortho ?thresh:(thresh=1.e-6) ~overlap c = let d, u, _ = gesvd (lacpy overlap) in diff --git a/Utils/Util.mli b/Utils/Util.mli index 915a057..b2ff95e 100644 --- a/Utils/Util.mli +++ b/Utils/Util.mli @@ -84,6 +84,9 @@ val diagonalize_symm : Lacaml.D.mat -> Lacaml.D.mat * Lacaml.D.vec val xt_o_x : o:Lacaml.D.mat -> x:Lacaml.D.mat -> Lacaml.D.mat (** Computes {% $\mathbf{X^\dag\, O\, X}$ %} *) +val x_o_xt : o:Lacaml.D.mat -> x:Lacaml.D.mat -> Lacaml.D.mat +(** Computes {% $\mathbf{X\, O\, X^\dag}$ %} *) + val canonical_ortho: ?thresh:float -> overlap:Lacaml.D.mat -> Lacaml.D.mat -> Lacaml.D.mat (** Canonical orthogonalization. [overlap] is the overlap matrix {% $\mathbf{S}$ %}, and the last argument contains the vectors {% $\mathbf{C}$ %} to orthogonalize. diff --git a/run_fci.ml b/run_fci.ml new file mode 100644 index 0000000..500ce3f --- /dev/null +++ b/run_fci.ml @@ -0,0 +1,68 @@ +let () = + let open Command_line in + begin + set_header_doc (Sys.argv.(0) ^ " - QCaml command"); + set_description_doc "Runs a Hartree-Fock calculation"; + set_specs + [ { short='b' ; long="basis" ; opt=Mandatory; + arg=With_arg ""; + doc="Name of the file containing the basis set"; } ; + + { short='x' ; long="xyz" ; opt=Mandatory; + arg=With_arg ""; + doc="Name of the file containing the nuclear coordinates in xyz format"; } ; + + { short='m' ; long="multiplicity" ; opt=Optional; + arg=With_arg ""; + doc="Spin multiplicity (2S+1). Default is singlet"; } ; + + { short='c' ; long="charge" ; opt=Optional; + arg=With_arg ""; + doc="Total charge of the molecule. Default is 0"; } ; + ] + end; + + (* Handle options *) + let basis_file = + match Command_line.get "basis" with + | Some x -> x + | None -> assert false + in + + let nuclei_file = + match Command_line.get "xyz" with + | Some x -> x + | None -> assert false + in + + let charge = + match Command_line.get "charge" with + | Some x -> int_of_string x + | None -> 0 + in + + let multiplicity = + match Command_line.get "multiplicity" with + | Some x -> int_of_string x + | None -> 1 + in + + let s = + Simulation.of_filenames ~charge ~multiplicity ~nuclei:nuclei_file basis_file + in + + let hf = HartreeFock.make s in + let mos = + MOBasis.of_hartree_fock hf + in + let space = + Determinant_space.fci_of_mo_basis ~frozen_core:false mos + in + let ci = CI.make space in + Format.printf "FCI energy : %20.16f@." ((CI.eigenvalues ci).{1} +. Simulation.nuclear_repulsion s); + (* + let s2 = Util.xt_o_x ~o:(CI.s2_matrix ci) ~x:(CI.eigenvectors ci) in + Util.list_range 1 (Determinant_space.size space) + |> List.iter (fun i -> Format.printf "@[%f@]@;" s2.{i,i}); + *) +