Fixed open shells CAS

CI/CI.ml

@@ -805,7 +805,7 @@ let pt2_en_reference ci =
 
   let aux_basis = mo_basis in
   let ds =
-    DeterminantSpace.fci_of_mo_basis ~frozen_core:true aux_basis
+    DeterminantSpace.fci_of_mo_basis ~frozen_core:false aux_basis
   in
   let out_dets =
     ds

CI/DeterminantSpace.ml

@@ -292,21 +292,21 @@ let arbitrary_of_mo_basis mo_basis f =
 
 
 
-let cas_of_mo_basis mo_basis n m =
+let cas_of_mo_basis mo_basis ~frozen_core n m =
   let f n_alfa =
-    Ss.cas_of_mo_basis mo_basis n_alfa n m
+    Ss.cas_of_mo_basis mo_basis ~frozen_core n_alfa n m
   in
   spin_of_mo_basis mo_basis f
 
 
-let fci_of_mo_basis ?(frozen_core=true) mo_basis =
+let fci_of_mo_basis mo_basis ~frozen_core =
   let f n_alfa =
-    Ss.fci_of_mo_basis ~frozen_core mo_basis n_alfa
+    Ss.fci_of_mo_basis mo_basis ~frozen_core n_alfa
   in
   spin_of_mo_basis mo_basis f
 
 
-let fci_f12_of_mo_basis ?(frozen_core=true) mo_basis mo_num =
+let fci_f12_of_mo_basis mo_basis ~frozen_core mo_num =
   let s = MOBasis.simulation mo_basis in
   let e = Simulation.electrons s in
   let n_alfa = Electrons.n_alfa e 
@@ -321,7 +321,7 @@ let fci_f12_of_mo_basis ?(frozen_core=true) mo_basis mo_num =
       (mo_num - n_core)
   in
   let f n_alfa =
-    Ss.cas_of_mo_basis mo_basis n_alfa n m
+    Ss.cas_of_mo_basis mo_basis ~frozen_core n_alfa n m
   in
   let r = 
     spin_of_mo_basis mo_basis f
@@ -335,9 +335,9 @@ let fci_f12_of_mo_basis ?(frozen_core=true) mo_basis mo_num =
     |> MOClass.of_list }
     
 
-let cas_f12_of_mo_basis mo_basis n m mo_num =
+let cas_f12_of_mo_basis mo_basis ~frozen_core n m mo_num =
   let f n_alfa =
-    Ss.cas_of_mo_basis mo_basis n_alfa n m
+    Ss.cas_of_mo_basis mo_basis ~frozen_core n_alfa n m
   in
   let r = 
     spin_of_mo_basis mo_basis f

CI/DeterminantSpace.mli

@@ -50,20 +50,20 @@ val fock_diag : t -> Determinant.t -> float array * float array
 *)
 
 
-val fci_of_mo_basis : ?frozen_core:bool -> MOBasis.t -> t
+val fci_of_mo_basis : MOBasis.t -> frozen_core:bool -> t
 (** Creates a space of all possible ways to put [n_alfa] electrons in the {% $\alpha$ %}
     [Active] MOs and [n_beta] electrons in the {% $\beta$ %} [Active] MOs.
     All other MOs are untouched.
 *)
 
-val cas_of_mo_basis : MOBasis.t -> int -> int -> t
+val cas_of_mo_basis : MOBasis.t -> frozen_core:bool -> int -> int -> t
 (** Creates a CAS(n,m) space of determinants.  *)
 
-val fci_f12_of_mo_basis : ?frozen_core:bool -> MOBasis.t -> int -> t
+val fci_f12_of_mo_basis : MOBasis.t -> frozen_core:bool -> int -> t
 (** Creates the active space to perform a FCI-F12 with an
     auxiliary basis set. *)
 
-val cas_f12_of_mo_basis : MOBasis.t -> int -> int -> int -> t
+val cas_f12_of_mo_basis : MOBasis.t -> frozen_core:bool -> int -> int -> int -> t
 (** [cas_of_mo_basis mo_basis m n mo_num] Creates a CAS(n,m) space
   of determinants with an auxiliary basis set defined as the MOs from
   [mo_num+1] to [MOBasis.size mo_basis].

CI/F12CI.ml

@@ -147,7 +147,7 @@ let dressing_vector gamma aux_basis f12_amplitudes ci =
 
 
 
-let make ~simulation ?(threshold=1.e-12) ?(frozen_core=true) ~mo_basis ~aux_basis_filename () =
+let make ~simulation ?(threshold=1.e-12) ~frozen_core ~mo_basis ~aux_basis_filename () =
 
   let gamma = 0.5 in
 

CI/SpindeterminantSpace.ml

@@ -13,7 +13,7 @@ let mo_basis t = t.mo_basis
 let mo_class t = t.mo_class
 let size     t = Array.length t.spin_determinants
 
-let fci_of_mo_basis ?(frozen_core=true) mo_basis elec_num =
+let fci_of_mo_basis ~frozen_core mo_basis elec_num =
   let mo_num = MOBasis.size mo_basis in
   let mo_class = MOClass.fci ~frozen_core mo_basis in
   let m l = 
@@ -35,9 +35,9 @@ let fci_of_mo_basis ?(frozen_core=true) mo_basis elec_num =
   { elec_num ; mo_basis ; mo_class ; spin_determinants }
 
 
-let cas_of_mo_basis mo_basis elec_num n m =
+let cas_of_mo_basis mo_basis ~frozen_core elec_num n m =
   let mo_num = MOBasis.size mo_basis in
-  let mo_class = MOClass.cas_sd mo_basis n m in
+  let mo_class = MOClass.cas_sd ~frozen_core mo_basis n m in
   let m l = 
     List.fold_left (fun accu i -> let j = i-1 in Z.(logor accu (shift_left one j)) 
                    ) Z.zero l

CI/SpindeterminantSpace.mli

@@ -24,12 +24,12 @@ val mo_basis : t -> MOBasis.t
 
 (** {1 Creation} *)
 
-val fci_of_mo_basis : ?frozen_core:bool -> MOBasis.t -> int -> t
+val fci_of_mo_basis : frozen_core:bool -> MOBasis.t -> int -> t
 (** Create a space of all possible ways to put [n_elec-ncore] electrons in the
     [Active] MOs.  All other MOs are untouched.
 *)
 
-val cas_of_mo_basis : MOBasis.t -> int -> int -> int -> t
+val cas_of_mo_basis : MOBasis.t -> frozen_core:bool -> int -> int -> int -> t
 (** [cas_of_mo_basis mo_basis n_elec n m] creates a CAS(n,m) space of 
     [Active] MOs. The unoccupied MOs are [Virtual], and the occupied MOs
     are [Core] and [Inactive].

MOBasis/MOClass.ml

@@ -104,17 +104,19 @@ let fci ?(frozen_core=true) mo_basis =
   )
 
 
-let cas_sd mo_basis n m =
+let cas_sd mo_basis ~frozen_core n m =
   let mo_num = MOBasis.size mo_basis in
   let n_alfa = MOBasis.simulation mo_basis |> Simulation.electrons |> Electrons.n_alfa in
   let n_beta = MOBasis.simulation mo_basis |> Simulation.electrons |> Electrons.n_beta in
   let n_unpaired = n_alfa - n_beta in
-  let n_alfa_in_cas = (n - n_unpaired)/2 in
+  let n_alfa_in_cas = (n - n_unpaired)/2 + n_unpaired in
   let last_inactive = n_alfa - n_alfa_in_cas in
   let last_active  = last_inactive + m in
   let ncore  =
-    (Nuclei.small_core @@ Simulation.nuclei @@ MOBasis.simulation mo_basis) / 2
-    |> min last_inactive
+    if frozen_core then
+      (Nuclei.small_core @@ Simulation.nuclei @@ MOBasis.simulation mo_basis) / 2
+      |> min last_inactive
+    else 0
   in
   of_list (
     List.concat [

MOBasis/MOClass.mli

@@ -20,7 +20,7 @@ val fci : ?frozen_core:bool -> MOBasis.t -> t
     [n] lowest MOs are [Core] if [frozen_core = true].
 *)
 
-val cas_sd: MOBasis.t -> int -> int -> t
+val cas_sd: MOBasis.t -> frozen_core:bool -> int -> int -> t
 (** [cas_sd mo_basis n m ] creates the MO classes for CAS(n,m) + SD
     calculations.  lowest MOs are [Core], then all the next MOs are [Inactive],
     then [Active], then [Virtual].

Perturbation/MP2.ml

@@ -1,6 +1,6 @@
 type t = float
 
-let make ?(frozen_core=true) hf =
+let make ~frozen_core hf =
   let mo_basis = 
     MOBasis.of_hartree_fock hf
   in
@@ -8,7 +8,7 @@ let make ?(frozen_core=true) hf =
     MOBasis.mo_energies mo_basis
   in
   let mo_class =
-    MOClass.cas_sd mo_basis 0 0
+    MOClass.cas_sd mo_basis ~frozen_core 0 0
     |> MOClass.to_list
   in
   let eri = 

run_cas.ml

@@ -67,7 +67,7 @@ let () =
   in
 
   let space =
-    DeterminantSpace.cas_of_mo_basis mos n m
+    DeterminantSpace.cas_of_mo_basis mos ~frozen_core:true n m
   in
   let ci = CI.make space in
   Format.fprintf ppf "CAS-CI energy : %20.16f@." ((CI.eigenvalues ci).{1} +. Simulation.nuclear_repulsion s);

run_mp2.ml

@@ -73,7 +73,7 @@ let () =
 
   let e_hf = HartreeFock.energy hf in
 
-  let mp2 = MP2.make hf in
+  let mp2 = MP2.make ~frozen_core:true hf in
   Format.fprintf ppf "@[MP2   = %15.10f@]@." mp2;
   Format.fprintf ppf "@[E+MP2 = %15.10f@]@." (mp2 +. e_hf)