open Linear_algebra
   
type t = float

let make ~frozen_core hf =
  let mo_basis = 
    Mo.Basis.of_hartree_fock hf
  in
  let epsilon =
    Mo.Basis.mo_energies mo_basis
  in
  let mo_class =
    Mo.Class.cas_sd mo_basis ~frozen_core 0 0
    |> Mo.Class.to_list
  in
  let eri = 
    Mo.Basis.ee_ints mo_basis
  in
  let inactives =
    List.filter (fun i ->
      match i with Mo.Class.Inactive _ -> true | _ -> false) mo_class
  and virtuals =
    List.filter (fun i ->
      match i with Mo.Class.Virtual _ -> true | _ -> false) mo_class
  in

  let rmp2 () = 
    List.fold_left (fun accu b ->
        match b with Mo.Class.Virtual b ->
          let eps = -. (epsilon%.(b)) in
          accu +. 
          List.fold_left (fun accu a ->
              match a with Mo.Class.Virtual a ->
                let eps = eps -. (epsilon%.(a)) in
                accu +. 
                List.fold_left (fun accu j ->
                    match j with Mo.Class.Inactive j ->
                      let eps = eps +. epsilon%.(j) in
                      accu +. 
                      List.fold_left (fun accu i ->
                          match i with Mo.Class.Inactive i ->
                            let eps = eps +. epsilon%.(i) in
                            let ijab = Four_idx_storage.get_phys eri i j a b
                            and abji = Four_idx_storage.get_phys eri a b j i in
                            let abij =  ijab in
                            accu +. ijab *. ( abij +. abij -. abji) /. eps 
                                    | _ -> accu
                        ) 0. inactives
                              | _ -> accu
                  ) 0. inactives
                        | _ -> accu
            ) 0. virtuals
                  | _ -> accu
      ) 0. virtuals
  in


  match Mo.Hartree_fock.kind hf with
  | Mo.Hartree_fock.RHF   -> rmp2 ()
  | _ -> failwith "Not implemented"