diff --git a/src/apps/molresponse/FrequencyResponse.cpp b/src/apps/molresponse/FrequencyResponse.cpp
index 6090a66ebc2..34779ea6538 100644
--- a/src/apps/molresponse/FrequencyResponse.cpp
+++ b/src/apps/molresponse/FrequencyResponse.cpp
@@ -545,15 +545,7 @@ X_space QuadraticResponse::compute_UPSILON(World &world, const X_space &XB, cons
     return UPSILON_BC;
 }
 
-X_space QuadraticResponse::compute_VBC(World &world, const X_space &XB, const X_space &XC) {
-
-    auto conj_XB = to_conjugate_X_space(to_response_matrix(XB));
-
-    auto conj_XC = to_conjugate_X_space(to_response_matrix(XC));
-
-    auto UPSILON_BC = XB * conj_XC + XC * conj_XB;
-    return UPSILON_BC;
-}
+X_space QuadraticResponse::compute_VBC(World &world, const X_space &VB, const X_space &XC) { return VB * XC; }
 
 
 Tensor<double> QuadraticResponse::compute_beta(World &world) {
@@ -574,10 +566,19 @@ Tensor<double> QuadraticResponse::compute_beta(World &world) {
     // now compute UPSILON BC
     auto PQ_BC = setup_rhs_BC(world, MUA);
 
+
+    auto PBC = PQ_BC.first * XC + XB * PQ_BC.second;
+
+
+    auto gamma_BC = compute_gamma_quadratic(world, XB, XC, ground_orbitals);
+    auto V_BC = gamma_BC + PBC;
+
     auto UPSILON_BC = compute_UPSILON(world, XB, XC);
+
     // create PQ vector for A
 
-    auto beta_one = inner(MUA, XA);
+    auto beta_one = inner(MUA, UPSILON_BC);
+    auto beta_two = inner(XA, V_BC);
 
 
     // now compute
@@ -585,8 +586,102 @@ Tensor<double> QuadraticResponse::compute_beta(World &world) {
     // B takes x[1] and x[2] and C takes x[0] and x[1]
 
 
-    return beta;
+    return beta_one + beta_two;
 }
+X_space QuadraticResponse::compute_JBC(World &world, const X_space &B, const X_space &C, const poperatorT &coulomb_ops,
+                                       const vector_real_function_3d &phi0) const {
+
+    X_space J = X_space::zero_functions(world, B.num_states(), B.num_orbitals());
+    auto c_xc = xcf.hf_exchange_coefficient();
+    QProjector<double, 3> projector(world, phi0);
+    vector_real_function_3d rhoX(B.num_states());
+    vector_real_function_3d temp_J(B.num_states());
+
+    vector_real_function_3d x_phi, y_phi;
 
+    for (const auto &b: B.active) {
 
+        x_phi = mul(world, B.x[b], phi0, false);
+        y_phi = mul(world, B.y[b], phi0, false);
+
+        world.gop.fence();
+        rhoX[b] = sum(world, x_phi, true);
+        rhoX[b] += sum(world, y_phi, true);
+        world.gop.fence();
+    }
+    truncate(world, rhoX);
+    for (int k = 0; k < B.num_states(); k++) {
+        temp_J[k] = apply(*coulomb_ops, rhoX[k]);
+        J.x[k] = mul(world, temp_J[k], C.x[k], false);
+        J.y[k] = mul(world, temp_J[k], C.y[k], false);
+    }
+    world.gop.fence();
+
+    return J;
+}
+
+auto Koperator(const vecfuncT &ket, const vecfuncT &bra) {
+    const double lo = 1.e-10;
+    auto &world = ket[0].world();
+    Exchange<double, 3> k{world, lo};
+    k.set_bra_and_ket(bra, ket);
+    k.set_algorithm(k.multiworld_efficient);
+    return k;
+};
+
+X_space QuadraticResponse::compute_KBC(World &world, const X_space &B, const X_space &C,
+                                       const vector_real_function_3d &phi0) const {
+
+    // if the frequecy of B is 0 we run the static case
+    // else we run the dynamic case
+    auto K = X_space::zero_functions(world, B.num_states(), B.num_orbitals());
+
+    vector_real_function_3d k1x, k1y, k2x, k2y;
+    vector_real_function_3d xb;
+    vector_real_function_3d yb;
+    if (world.rank() == 0) { print("K1StrategyFull"); }
+
+    auto k1x_temp = create_response_matrix(B.num_states(), B.num_orbitals());
+    auto k1y_temp = create_response_matrix(B.num_states(), B.num_orbitals());
+    auto k2x_temp = create_response_matrix(B.num_states(), B.num_orbitals());
+    auto k2y_temp = create_response_matrix(B.num_states(), B.num_orbitals());
+
+
+    for (int k = 0; k < B.num_states(); k++) {
+
+        auto K1X = Koperator(B.x[k], phi0);
+        auto K1Y = Koperator(phi0, B.y[k]);
+
+        auto K2X = Koperator(B.y[k], phi0);
+        auto K2Y = Koperator(phi0, B.x[k]);
+        world.gop.fence();
+
+        k1x_temp[k] = K1X(C.x[k]);
+        k1y_temp[k] = K1Y(C.x[k]);
+        k2x_temp[k] = K2X(C.y[k]);
+        k2y_temp[k] = K2Y(C.y[k]);
+        world.gop.fence();
+        K.x[k] = gaxpy_oop(1.0, k1x_temp[k], 1.0, k1y_temp[k], false);
+        K.y[k] = gaxpy_oop(1.0, k2x_temp[k], 1.0, k2y_temp[k], false);
+    }
+
+    return K;
+}
+
+auto QuadraticResponse::compute_gamma_quadratic(World &world, const X_space &XB, const X_space &XC,
+                                                const vector_real_function_3d &phi0) const -> X_space {
+
+    auto JBC = compute_JBC(world, XB, XC, shared_coulomb_operator, phi0);
+    auto JCB = compute_JBC(world, XC, XB, shared_coulomb_operator, phi0);
+
+    auto KBC = compute_KBC(world, XB, XC, phi0);
+    auto KCB = compute_KBC(world, XC, XB, phi0);
+
+    auto gamma = X_space::zero_functions(world, XB.num_states(), XB.num_orbitals());
+    gamma = JBC + JCB + KBC + KCB;
+
+
+    return gamma;
+    // Get sizes
+}
 //
diff --git a/src/apps/molresponse/FrequencyResponse.hpp b/src/apps/molresponse/FrequencyResponse.hpp
index 09f33453320..59260bc05e0 100644
--- a/src/apps/molresponse/FrequencyResponse.hpp
+++ b/src/apps/molresponse/FrequencyResponse.hpp
@@ -199,6 +199,11 @@ class QuadraticResponse : public ResponseBase {
     RHS_Generator generator;
     X_space compute_VBC(World &world, const X_space &XB, const X_space &XC);
     std::pair<X_space, X_space> setup_rhs_BC(World &world, const X_space &PQ);
+    X_space compute_gamma_quadratic(World &world, const X_space &XB, const X_space &XC,
+                                    const vector_real_function_3d &phi0) const;
+    X_space compute_JBC(World &world, const X_space &B, const X_space &C, const poperatorT &coulomb_ops,
+                        const vector_real_function_3d &phi0) const;
+    X_space compute_KBC(World &world, const X_space &B, const X_space &C, const vector_real_function_3d &phi0) const;
 };
 
 class FrequencyResponse : public ResponseBase {
diff --git a/src/apps/molresponse/x_space.cc b/src/apps/molresponse/x_space.cc
index fa52a49dc35..d1bed0ceb80 100644
--- a/src/apps/molresponse/x_space.cc
+++ b/src/apps/molresponse/x_space.cc
@@ -168,7 +168,7 @@ namespace madness {
     auto inner(const X_space &A, const X_space &B) -> Tensor<double> {
         MADNESS_ASSERT(size_states(A) > 0);
         MADNESS_ASSERT(size_orbitals(A) > 0);
-        MADNESS_ASSERT(same_size(A, B));
+        MADNESS_ASSERT(A.num_orbitals()==B.num_orbitals());
 
         // return response_space_inner(A.x, B.x) + response_space_inner(A.y, B.y);
 
@@ -181,7 +181,7 @@ namespace madness {
         auto b_transpose = transposeResponseMatrix(b);
 
         world.gop.fence();
-        Tensor<double> result(a.size(), a.size());
+        Tensor<double> result(A.num_states(), B.num_states());
         int p = 0;
         std::for_each(a_transpose.begin(), a_transpose.end(),
                       [&](const auto &ati) {