matrix.cpp

#include "matrix.h"

MatrixXcd create_first_boundary_matrix(int M, int N_tot)
{
	// function creates the first boundary condition matrix U

	MatrixXcd U = MatrixXcd::Zero(M, N_tot * M);

	for (int i = 0; i < M; ++i)
	{
		for (int j = 0; j < 4 * M; ++j)
		{
			if (j % M == i)
				U(i, j) = complex<double> (1.0);
		}
	}

	return U;
}

/****************************************************	
 * DNA ends z_0-axis collinear boundary condition: 
 * 		ref: Eq. E53;
 * 		(without exp(b_0*f))	
 * **************************************************/

MatrixXcd create_matrix_V(int max_mode)
{
	MatrixXcd V = MatrixXcd::Zero(max_mode, max_mode);

	for (int p = 0; p < max_mode; ++p)
		for (int q = 0; q < max_mode; ++q)
			V(p, q) = complex<double> (std::sqrt((2 * p + 1) * (2 * q + 1)) / 4.0 / pi);

	return V;
}

MatrixXcd create_second_boundary_matrix(Params& p, double mu_L, double mu_P, double mu_S, double f, int M, int N_tot, MatrixXcd& V)
{
	MatrixXcd O = MatrixXcd::Zero(N_tot * M, M);

	O.block(0,     0, M, M) = std::exp(p.b_S * f) * std::exp(-p.q * (mu_S - 2.0 * pi * p.tau * p.lk_S_0)) * V;
	O.block(M,     0, M, M) = std::exp(p.b_P * f) * std::exp(-p.q * (mu_P - 2.0 * pi * p.tau * p.lk_P_0)) * V;
	O.block(2 * M, 0, M, M) = std::exp(p.b_L * f) * std::exp(-p.q * (mu_L - 2.0 * pi * p.tau * p.lk_L_0)) * V;
	O.block(3 * M, 0, M, M) = std::exp(p.b_B * f) * V;

	return O;
}

/******************************************************************
 *	Creating small matrix blocks,
 *	index [max_mode - 1] represents [0] in equations:
*******************************************************************/

vector<MatrixXcd> create_matrix_block_SB_SL_SP_SS_Sprotein(Params& ps, int max_mode, vector<double>& w3j_all,
	vector2d& L_0, vector2d& L_a_B, vector2d& L_b_B_f, vector2d& L_b_B_f_plus, vector2d& L_b_B_f_minus, vector2d& L_a_L,
	vector2d& L_b_L_f, vector2d& L_b_L_f_plus, vector2d& L_b_L_f_minus, vector2d& L_a_P, vector2d& L_b_P_f,
	vector2d& L_b_P_f_plus, vector2d& L_b_P_f_minus, vector2d& L_a_S, vector2d& L_b_S_f, vector2d& L_b_S_f_plus,
	vector2d& L_b_S_f_minus, vector2d& L_a_protein, vectorcd& D_n_ml_A_in_octamer, vectorcd& D_n_ml_A_out_octamer,
	vectorcd& D_n_ml_A_in_L_tetramer, vectorcd& D_n_ml_A_out_L_tetramer, vectorcd& D_n_ml_A_in_R_tetramer,
	vectorcd& D_n_ml_A_out_R_tetramer, vectorcd& D_n_ml_A_in_hexamer, vectorcd& D_n_ml_A_out_hexamer)
{
	MatrixXcd m_S_B = MatrixXcd::Zero(max_mode, max_mode);			// S_B is S_00 for B-DNA;
	MatrixXcd m_S_L = MatrixXcd::Zero(max_mode, max_mode);			// S_L is S_00 for L-DNA;
	MatrixXcd m_S_P = MatrixXcd::Zero(max_mode, max_mode);			// S_P is S_00 for P-DNA;
	MatrixXcd m_S_S = MatrixXcd::Zero(max_mode, max_mode);			// S_P is S_00 for P-DNA;

	MatrixXcd m_S_A_in_octamer  = MatrixXcd::Zero(ps.M, ps.M);	// m_S_A_in is S_01 
	MatrixXcd m_S_A_out_octamer = MatrixXcd::Zero(ps.M, ps.M);	// m_S_A_out is S_K0 
	MatrixXcd m_S_octamer       = MatrixXcd::Zero(ps.M, ps.M);	// m_S_octamer is S_12 = S_23 

	MatrixXcd m_S_A_in_L_tetramer  = MatrixXcd::Zero(ps.M, ps.M);	// m_S_A_in is S_01 
	MatrixXcd m_S_A_out_L_tetramer = MatrixXcd::Zero(ps.M, ps.M);	// m_S_A_out is S_K0 
	MatrixXcd m_S_tetramer         = MatrixXcd::Zero(ps.M, ps.M);	// m_S_L_tetramer is S_12 = S_23 

	MatrixXcd m_S_A_in_R_tetramer  = MatrixXcd::Zero(ps.M, ps.M);	// m_S_A_in is S_01 
	MatrixXcd m_S_A_out_R_tetramer = MatrixXcd::Zero(ps.M, ps.M);	// m_S_A_out is S_K0

	MatrixXcd m_S_A_in_hexamer  = MatrixXcd::Zero(ps.M, ps.M);	// m_S_A_in is S_01 
	MatrixXcd m_S_A_out_hexamer = MatrixXcd::Zero(ps.M, ps.M);	// m_S_A_out is S_K0 
	MatrixXcd m_S_hexamer       = MatrixXcd::Zero(ps.M, ps.M);	// m_S_hexamer is S_12 = S_23 

	MatrixXcd m_S_B_plus_f = MatrixXcd::Zero(max_mode, max_mode);
	MatrixXcd m_S_L_plus_f = MatrixXcd::Zero(max_mode, max_mode);
	MatrixXcd m_S_P_plus_f = MatrixXcd::Zero(max_mode, max_mode);
	MatrixXcd m_S_S_plus_f = MatrixXcd::Zero(max_mode, max_mode);

	MatrixXcd m_S_B_plus_t = MatrixXcd::Zero(max_mode, max_mode);
	MatrixXcd m_S_L_plus_t = MatrixXcd::Zero(max_mode, max_mode);
	MatrixXcd m_S_P_plus_t = MatrixXcd::Zero(max_mode, max_mode);
	MatrixXcd m_S_S_plus_t = MatrixXcd::Zero(max_mode, max_mode);

	MatrixXcd m_S_A_in_plus_f_octamer  = MatrixXcd::Zero(ps.M, ps.M);
	MatrixXcd m_S_A_out_plus_f_octamer = MatrixXcd::Zero(ps.M, ps.M);
	MatrixXcd m_S_plus_f_octamer       = MatrixXcd::Zero(ps.M, ps.M);

	MatrixXcd m_S_A_in_plus_f_L_tetramer  = MatrixXcd::Zero(ps.M, ps.M);
	MatrixXcd m_S_A_out_plus_f_L_tetramer = MatrixXcd::Zero(ps.M, ps.M);
	MatrixXcd m_S_plus_f_tetramer         = MatrixXcd::Zero(ps.M, ps.M);

	MatrixXcd m_S_A_in_plus_f_R_tetramer  = MatrixXcd::Zero(ps.M, ps.M);
	MatrixXcd m_S_A_out_plus_f_R_tetramer = MatrixXcd::Zero(ps.M, ps.M);

	MatrixXcd m_S_A_in_plus_f_hexamer  = MatrixXcd::Zero(ps.M, ps.M);
	MatrixXcd m_S_A_out_plus_f_hexamer = MatrixXcd::Zero(ps.M, ps.M);
	MatrixXcd m_S_plus_f_hexamer       = MatrixXcd::Zero(ps.M, ps.M);

	MatrixXcd m_S_B_minus_f = MatrixXcd::Zero(max_mode, max_mode);
	MatrixXcd m_S_L_minus_f = MatrixXcd::Zero(max_mode, max_mode);
	MatrixXcd m_S_P_minus_f = MatrixXcd::Zero(max_mode, max_mode);
	MatrixXcd m_S_S_minus_f = MatrixXcd::Zero(max_mode, max_mode);

	MatrixXcd m_S_B_minus_t = MatrixXcd::Zero(max_mode, max_mode);
	MatrixXcd m_S_L_minus_t = MatrixXcd::Zero(max_mode, max_mode);
	MatrixXcd m_S_P_minus_t = MatrixXcd::Zero(max_mode, max_mode);
	MatrixXcd m_S_S_minus_t = MatrixXcd::Zero(max_mode, max_mode);

	MatrixXcd m_S_A_in_minus_f_octamer  = MatrixXcd::Zero(ps.M, ps.M);
	MatrixXcd m_S_A_out_minus_f_octamer = MatrixXcd::Zero(ps.M, ps.M);
	MatrixXcd m_S_minus_f_octamer       = MatrixXcd::Zero(ps.M, ps.M);

	MatrixXcd m_S_A_in_minus_f_L_tetramer  = MatrixXcd::Zero(ps.M, ps.M);
	MatrixXcd m_S_A_out_minus_f_L_tetramer = MatrixXcd::Zero(ps.M, ps.M);
	MatrixXcd m_S_minus_f_tetramer       = MatrixXcd::Zero(ps.M, ps.M);

	MatrixXcd m_S_A_in_minus_f_R_tetramer  = MatrixXcd::Zero(ps.M, ps.M);
	MatrixXcd m_S_A_out_minus_f_R_tetramer = MatrixXcd::Zero(ps.M, ps.M);

	MatrixXcd m_S_A_in_minus_f_hexamer  = MatrixXcd::Zero(ps.M, ps.M);
	MatrixXcd m_S_A_out_minus_f_hexamer = MatrixXcd::Zero(ps.M, ps.M);
	MatrixXcd m_S_minus_f_hexamer       = MatrixXcd::Zero(ps.M, ps.M);

	for (int n1 = 0; n1 < max_mode; ++n1)
	{// 'n1' is 's' in Eq. E32;
		int index_1 = n1 * (n1 + 1);

		for (int n2 = 0; n2 < max_mode; ++n2)
		{// 'n2' is 's'' in Eq. E32;
			double prefactor_1 = std::sqrt((2.0 * n1 + 1.0) * (2.0 * n2 + 1.0));
			int index_2 = n2 * (n2 + 1);

			/***************** Bare DNA *****************/

			for (int p = 0; p < max_mode; ++p)
			{// 'p' is 't' in Eq. E32;
				double prefactor_2 = prefactor_1 * (2.0 * p + 1.0);

				for (int k1 = std::abs(n1 - p); k1 <= n1 + p; ++k1)
				{// 'k1' is 'k' in Eq. E32;
					double prefactor_3 = prefactor_2 * (2.0 * k1 + 1.0);

					for (int k2 = std::abs(n2 - p); k2 <= n2 + p; ++k2)
					{// 'k2' is 'k'' in Eq. E32;
						double prefactor_4 = prefactor_3 * (2.0 * k2 + 1.0);

						for (int r = -std::min({ p, k1, k2 }); r <= std::min({ p, k1, k2 }); ++r)
						{// 'r' is 'r' in Eq. E32;
							int ind_r = std::abs(r);
							int ind_pr = max_mode + r - 1;
							int ind_mr = max_mode - r - 1;

							double w_3j_sym = prefactor_4 * L_0[k2][ind_pr] * std::pow(w3j_all[index(n1, p, k1, 
								max_mode - 1, ind_mr, max_mode - 1)] * w3j_all[index(n2, p, k2, max_mode - 1, ind_mr, max_mode - 1)], 2.0);

							/***************** B - DNA *****************/

							complex<double> matr_elem_general = ps.coeff_1i_all_B[ind_pr] * ps.coeff_omega_B[ind_pr] *
								ps.besseli_c_B[ind_r] * L_a_B[p][ind_pr] * w_3j_sym;

							complex<double> matr_elem = matr_elem_general * ps.besseli_coeff_B[ind_r];

							m_S_B(n1, n2)         += matr_elem * L_b_B_f[k1][ind_pr];
							m_S_B_plus_f(n1, n2)  += matr_elem * L_b_B_f_plus[k1][ind_pr];
							m_S_B_minus_f(n1, n2) += matr_elem * L_b_B_f_minus[k1][ind_pr];

							matr_elem = matr_elem_general * ps.besseli_coeff_plus_B[ind_r];
							m_S_B_plus_t(n1, n2)  += matr_elem * L_b_B_f[k1][ind_pr];
							
							matr_elem = matr_elem_general * ps.besseli_coeff_minus_B[ind_r];
							m_S_B_minus_t(n1, n2) += matr_elem * L_b_B_f[k1][ind_pr];

							/***************** L - DNA *****************/

							matr_elem_general = ps.coeff_1i_all_L[ind_pr] * ps.coeff_omega_L[ind_pr] *
								ps.besseli_c_L[ind_r] * L_a_L[p][ind_pr] * w_3j_sym;

							matr_elem = matr_elem_general * ps.besseli_coeff_L[ind_r];

							m_S_L(n1, n2)         += matr_elem * L_b_L_f[k1][ind_pr];
							m_S_L_plus_f(n1, n2)  += matr_elem * L_b_L_f_plus[k1][ind_pr];
							m_S_L_minus_f(n1, n2) += matr_elem * L_b_L_f_minus[k1][ind_pr];

							matr_elem = matr_elem_general * ps.besseli_coeff_plus_L[ind_r];
							m_S_L_plus_t(n1, n2)  += matr_elem * L_b_L_f[k1][ind_pr];
							
							matr_elem = matr_elem_general * ps.besseli_coeff_minus_L[ind_r];
							m_S_L_minus_t(n1, n2) += matr_elem * L_b_L_f[k1][ind_pr];

							/***************** P - DNA *****************/

							matr_elem_general = ps.coeff_1i_all_P[ind_pr] * ps.coeff_omega_P[ind_pr] *
								ps.besseli_c_P[ind_r] * L_a_P[p][ind_pr] * w_3j_sym;

							matr_elem = matr_elem_general * ps.besseli_coeff_P[ind_r];

							m_S_P(n1, n2)         += matr_elem * L_b_P_f[k1][ind_pr];
							m_S_P_plus_f(n1, n2)  += matr_elem * L_b_P_f_plus[k1][ind_pr];
							m_S_P_minus_f(n1, n2) += matr_elem * L_b_P_f_minus[k1][ind_pr];

							matr_elem = matr_elem_general * ps.besseli_coeff_plus_P[ind_r];
							m_S_P_plus_t(n1, n2)  += matr_elem * L_b_P_f[k1][ind_pr];

							matr_elem = matr_elem_general * ps.besseli_coeff_minus_P[ind_r];
							m_S_P_minus_t(n1, n2) += matr_elem * L_b_P_f[k1][ind_pr];

							/***************** S - DNA *****************/

							matr_elem_general = ps.coeff_1i_all_S[ind_pr] * ps.coeff_omega_S[ind_pr] *
								ps.besseli_c_S[ind_r] * L_a_S[p][ind_pr] * w_3j_sym;

							matr_elem = matr_elem_general * ps.besseli_coeff_S[ind_r];

							m_S_S(n1, n2)         += matr_elem * L_b_S_f[k1][ind_pr];
							m_S_S_plus_f(n1, n2)  += matr_elem * L_b_S_f_plus[k1][ind_pr];
							m_S_S_minus_f(n1, n2) += matr_elem * L_b_S_f_minus[k1][ind_pr];

							matr_elem = matr_elem_general * ps.besseli_coeff_plus_S[ind_r];
							m_S_S_plus_t(n1, n2)  += matr_elem * L_b_S_f[k1][ind_pr];
							
							matr_elem = matr_elem_general * ps.besseli_coeff_minus_S[ind_r];
							m_S_S_minus_t(n1, n2) += matr_elem * L_b_S_f[k1][ind_pr];
						}
					}
				}
			}

			/***************** S_01 block *****************/

			double multiplier(0.0), multiplier_plus(0.0), multiplier_minus(0.0);

			for (int k = std::abs(n1 - n2); k <= n1 + n2; ++k)
			{
				double matr_elem = (2.0 * k + 1) * std::pow(w3j_all[index(n1, n2, k, max_mode - 1, max_mode - 1, max_mode - 1)], 2.0);

				multiplier       += matr_elem * L_b_B_f[k][max_mode - 1];
				multiplier_plus  += matr_elem * L_b_B_f_plus[k][max_mode - 1];
				multiplier_minus += matr_elem * L_b_B_f_minus[k][max_mode - 1];
			}

			for (int t = -n2; t <= n2; ++t)
			{
				int index_2b = index_2 + t;
				int ind_t = max_mode + t - 1;
				int abs_t = std::abs(t);
				complex<double> prefactor_2 = 4.0 * std::pow(pi, 3.0) * prefactor_1 * ps.besseli_coeff_c_protein[abs_t] *
					L_a_protein[n2][ind_t];
				
				complex<double> prefactor_2_octamer = prefactor_2 * D_n_ml_A_in_octamer[index(n2, max_mode - 1, ind_t, max_mode)];
				m_S_A_in_octamer(index_1, index_2b)         += multiplier       * prefactor_2_octamer;
				m_S_A_in_plus_f_octamer(index_1, index_2b)  += multiplier_plus  * prefactor_2_octamer;
				m_S_A_in_minus_f_octamer(index_1, index_2b) += multiplier_minus * prefactor_2_octamer;

				complex<double> prefactor_2_L_tetramer = prefactor_2 * D_n_ml_A_in_L_tetramer[index(n2, max_mode - 1, ind_t, max_mode)];
				m_S_A_in_L_tetramer(index_1, index_2b)         += multiplier       * prefactor_2_L_tetramer;
				m_S_A_in_plus_f_L_tetramer(index_1, index_2b)  += multiplier_plus  * prefactor_2_L_tetramer;
				m_S_A_in_minus_f_L_tetramer(index_1, index_2b) += multiplier_minus * prefactor_2_L_tetramer;

				complex<double> prefactor_2_R_tetramer = prefactor_2 * D_n_ml_A_in_R_tetramer[index(n2, max_mode - 1, ind_t, max_mode)];
				m_S_A_in_R_tetramer(index_1, index_2b)         += multiplier       * prefactor_2_R_tetramer;
				m_S_A_in_plus_f_R_tetramer(index_1, index_2b)  += multiplier_plus  * prefactor_2_R_tetramer;
				m_S_A_in_minus_f_R_tetramer(index_1, index_2b) += multiplier_minus * prefactor_2_R_tetramer;

				complex<double> prefactor_2_hexamer = prefactor_2 * D_n_ml_A_in_hexamer[index(n2, max_mode - 1, ind_t, max_mode)];
				m_S_A_in_hexamer(index_1, index_2b)         += multiplier       * prefactor_2_hexamer;
				m_S_A_in_plus_f_hexamer(index_1, index_2b)  += multiplier_plus  * prefactor_2_hexamer;
				m_S_A_in_minus_f_hexamer(index_1, index_2b) += multiplier_minus * prefactor_2_hexamer;
			}

			/***************** S_12 block *****************/

			for (int t = -std::min(n1, n2); t <= std::min(n1, n2); ++t)
			{
				int index_1b = index_1 + t;
				int index_2b = index_2 + t;
				int ind_t = max_mode + t - 1;
				int ind_mt = max_mode - t - 1;
				complex<double> prefactor_2 = prefactor_1 * ps.coeff_S12[ind_t];

				for (int k = std::abs(n1 - n2); k <= n1 + n2; ++k)
				{
					complex<double> matr_elem = prefactor_2 * w3j_all[index(n1, n2, k, max_mode - 1, max_mode - 1,
						max_mode - 1)] * w3j_all[index(n1, n2, k, ind_mt, ind_t, max_mode - 1)] * (2.0 * k + 1);

					m_S_octamer(index_1b, index_2b)         += matr_elem * ps.spherical_besseli_coeff_r_pr_f_octamer[k];
					m_S_plus_f_octamer(index_1b, index_2b)  += matr_elem * ps.spherical_besseli_coeff_r_pr_f_plus_octamer[k];
					m_S_minus_f_octamer(index_1b, index_2b) += matr_elem * ps.spherical_besseli_coeff_r_pr_f_minus_octamer[k];

					m_S_tetramer(index_1b, index_2b)         += matr_elem * ps.spherical_besseli_coeff_r_pr_f_tetrasome[k];
					m_S_plus_f_tetramer(index_1b, index_2b)  += matr_elem * ps.spherical_besseli_coeff_r_pr_f_plus_tetrasome[k];
					m_S_minus_f_tetramer(index_1b, index_2b) += matr_elem * ps.spherical_besseli_coeff_r_pr_f_minus_tetrasome[k];

					m_S_hexamer(index_1b, index_2b)         += matr_elem * ps.spherical_besseli_coeff_r_pr_f_hexamer[k];
					m_S_plus_f_hexamer(index_1b, index_2b)  += matr_elem * ps.spherical_besseli_coeff_r_pr_f_plus_hexamer[k];
					m_S_minus_f_hexamer(index_1b, index_2b) += matr_elem * ps.spherical_besseli_coeff_r_pr_f_minus_hexamer[k];
				}
			}

			/***************** S_K0 block *****************/ 

			for (int t = -std::min(n1, n2); t <= std::min(n1, n2); ++t)
			{
				int index_1b = index_1 + t;
				int ind_t = max_mode + t - 1;
				int ind_mt = max_mode - t - 1;
				complex<double> prefactor_2 = 4.0 * std::pow(pi, 2.0) * prefactor_1 * ps.coeff_S12[ind_t] * ps.besseli_coeff_c_protein[0] *
					L_a_protein[n2][max_mode - 1];

				for (int k = std::abs(n1 - n2); k <= n1 + n2; ++k)
				{
					complex<double> matr_elem = prefactor_2 * w3j_all[index(n1, n2, k, max_mode - 1, max_mode - 1,
						max_mode - 1)] * w3j_all[index(n1, n2, k, ind_mt, ind_t, max_mode - 1)] * (2.0 * k + 1);

					complex<double> matr_elem_octamer = matr_elem * D_n_ml_A_out_octamer[index(n2, ind_t, max_mode - 1, max_mode)];
					m_S_A_out_octamer(index_1b, index_2)         += matr_elem_octamer * ps.spherical_besseli_coeff_r_pr_f_octamer[k];
					m_S_A_out_plus_f_octamer(index_1b, index_2)  += matr_elem_octamer * ps.spherical_besseli_coeff_r_pr_f_plus_octamer[k];
					m_S_A_out_minus_f_octamer(index_1b, index_2) += matr_elem_octamer * ps.spherical_besseli_coeff_r_pr_f_minus_octamer[k];

					complex<double> matr_elem_L_tetramer = matr_elem * D_n_ml_A_out_L_tetramer[index(n2, ind_t, max_mode - 1, max_mode)];
					m_S_A_out_L_tetramer(index_1b, index_2)         += matr_elem_L_tetramer * ps.spherical_besseli_coeff_r_pr_f_tetrasome[k];
					m_S_A_out_plus_f_L_tetramer(index_1b, index_2)  += matr_elem_L_tetramer * ps.spherical_besseli_coeff_r_pr_f_plus_tetrasome[k];
					m_S_A_out_minus_f_L_tetramer(index_1b, index_2) += matr_elem_L_tetramer * ps.spherical_besseli_coeff_r_pr_f_minus_tetrasome[k];

					complex<double> matr_elem_R_tetramer = matr_elem * D_n_ml_A_out_R_tetramer[index(n2, ind_t, max_mode - 1, max_mode)];
					m_S_A_out_R_tetramer(index_1b, index_2)         += matr_elem_R_tetramer * ps.spherical_besseli_coeff_r_pr_f_tetrasome[k];
					m_S_A_out_plus_f_R_tetramer(index_1b, index_2)  += matr_elem_R_tetramer * ps.spherical_besseli_coeff_r_pr_f_plus_tetrasome[k];
					m_S_A_out_minus_f_R_tetramer(index_1b, index_2) += matr_elem_R_tetramer * ps.spherical_besseli_coeff_r_pr_f_minus_tetrasome[k];

					complex<double> matr_elem_hexamer = matr_elem * D_n_ml_A_out_hexamer[index(n2, ind_t, max_mode - 1, max_mode)];
					m_S_A_out_hexamer(index_1b, index_2)         += matr_elem_hexamer * ps.spherical_besseli_coeff_r_pr_f_hexamer[k];
					m_S_A_out_plus_f_hexamer(index_1b, index_2)  += matr_elem_hexamer * ps.spherical_besseli_coeff_r_pr_f_plus_hexamer[k];
					m_S_A_out_minus_f_hexamer(index_1b, index_2) += matr_elem_hexamer * ps.spherical_besseli_coeff_r_pr_f_minus_hexamer[k];
				}
			}
		}
	}

	double energy_coeff_octamer = -double(ps.octamer_L) * ps.block_offset_energy + ps.octamer_mu + ps.octamer_mu_correction + 2.0 * pi * ps.tau * ps.octamer_Lk;
	MatrixXcd m_S_oct         = create_protein_matrix(m_S_A_in_octamer,         m_S_A_out_octamer,         m_S_octamer,         energy_coeff_octamer, ps.octamer_size, max_mode);
	MatrixXcd m_S_plus_f_oct  = create_protein_matrix(m_S_A_in_plus_f_octamer,  m_S_A_out_plus_f_octamer,  m_S_plus_f_octamer,  energy_coeff_octamer, ps.octamer_size, max_mode);
	MatrixXcd m_S_minus_f_oct = create_protein_matrix(m_S_A_in_minus_f_octamer, m_S_A_out_minus_f_octamer, m_S_minus_f_octamer, energy_coeff_octamer, ps.octamer_size, max_mode);

	double energy_coeff_L_tetramer = -double(ps.tetrasome_L) * ps.block_offset_energy + ps.L_tetrasome_mu + ps.L_tetrasome_mu_correction + 2.0 * pi * ps.tau * ps.L_tetrasome_Lk;
	MatrixXcd m_S_L_tet         = create_protein_matrix(m_S_A_in_L_tetramer,         m_S_A_out_L_tetramer,         m_S_tetramer,         energy_coeff_L_tetramer, ps.tetrasome_size, max_mode);
	MatrixXcd m_S_plus_f_L_tet  = create_protein_matrix(m_S_A_in_plus_f_L_tetramer,  m_S_A_out_plus_f_L_tetramer,  m_S_plus_f_tetramer,  energy_coeff_L_tetramer, ps.tetrasome_size, max_mode);
	MatrixXcd m_S_minus_f_L_tet = create_protein_matrix(m_S_A_in_minus_f_L_tetramer, m_S_A_out_minus_f_L_tetramer, m_S_minus_f_tetramer, energy_coeff_L_tetramer, ps.tetrasome_size, max_mode);

	double energy_coeff_R_tetramer = -double(ps.tetrasome_L) * ps.block_offset_energy + ps.R_tetrasome_mu + ps.R_tetrasome_mu_correction + 2.0 * pi * ps.tau * ps.R_tetrasome_Lk;
	MatrixXcd m_S_R_tet         = create_protein_matrix(m_S_A_in_R_tetramer,         m_S_A_out_R_tetramer,         m_S_tetramer,         energy_coeff_R_tetramer, ps.tetrasome_size, max_mode);
	MatrixXcd m_S_plus_f_R_tet  = create_protein_matrix(m_S_A_in_plus_f_R_tetramer,  m_S_A_out_plus_f_R_tetramer,  m_S_plus_f_tetramer,  energy_coeff_R_tetramer, ps.tetrasome_size, max_mode);
	MatrixXcd m_S_minus_f_R_tet = create_protein_matrix(m_S_A_in_minus_f_R_tetramer, m_S_A_out_minus_f_R_tetramer, m_S_minus_f_tetramer, energy_coeff_R_tetramer, ps.tetrasome_size, max_mode);

	double energy_coeff_hexamer = -double(ps.hexamer_L) * ps.block_offset_energy + ps.hexamer_mu + ps.hexamer_mu_correction + 2.0 * pi * ps.tau * ps.hexamer_Lk;
	MatrixXcd m_S_hex         = create_protein_matrix(m_S_A_in_hexamer,         m_S_A_out_hexamer,         m_S_hexamer,         energy_coeff_hexamer, ps.hexamer_size, max_mode);
	MatrixXcd m_S_plus_f_hex  = create_protein_matrix(m_S_A_in_plus_f_hexamer,  m_S_A_out_plus_f_hexamer,  m_S_plus_f_hexamer,  energy_coeff_hexamer, ps.hexamer_size, max_mode);
	MatrixXcd m_S_minus_f_hex = create_protein_matrix(m_S_A_in_minus_f_hexamer, m_S_A_out_minus_f_hexamer, m_S_minus_f_hexamer, energy_coeff_hexamer, ps.hexamer_size, max_mode);

	return vector<MatrixXcd> {m_S_B, m_S_L, m_S_P, m_S_S, m_S_oct, m_S_L_tet, m_S_R_tet, m_S_hex,
		m_S_B_plus_f, m_S_L_plus_f, m_S_P_plus_f, m_S_S_plus_f, m_S_plus_f_oct, m_S_plus_f_L_tet, m_S_plus_f_R_tet, m_S_plus_f_hex,
		m_S_B_minus_f, m_S_L_minus_f, m_S_P_minus_f, m_S_S_minus_f, m_S_minus_f_oct, m_S_minus_f_L_tet, m_S_minus_f_R_tet, m_S_minus_f_hex,
		m_S_B_plus_t, m_S_L_plus_t, m_S_P_plus_t, m_S_S_plus_t,
		m_S_B_minus_t, m_S_L_minus_t, m_S_P_minus_t, m_S_S_minus_t};
}


// Just a multiplicatin of all protein-related S matrices? - Yes. 
MatrixXcd create_protein_matrix(const MatrixXcd &m_S_protein_A_in, const MatrixXcd &m_S_protein_A_out,	
	const MatrixXcd &m_S_protein_B, double energy_coeff, double protein_size, int max_mode)
{
	// Function calculates the protein block in the DNA transfer-matrix  

	complex<double> log_coeff;
	MatrixXcd S_protein = MatrixXcd::Zero(max_mode, max_mode);
	MatrixXcd S_protein_temp = m_S_protein_A_in * quick_matrices_product_v3(m_S_protein_B, protein_size - 1, log_coeff) * m_S_protein_A_out;
	S_protein_temp = S_protein_temp * std::exp(energy_coeff + log_coeff);

//	MatrixXcd S_protein_temp = m_S_protein_A_in * matr_pow(m_S_protein_B, protein_size - 1) * m_S_protein_A_out; 

	for (int i = 0; i < max_mode; ++i)
		for (int j = 0; j < max_mode; ++j)
			S_protein(i, j) = S_protein_temp(i * (i + 1), j * (j + 1));

	return S_protein;
}

/*	This function creates the main transfer matrix:  */
MatrixXcd create_transfer_matrix(const MatrixXcd& m_S_B, const MatrixXcd& m_S_L, const MatrixXcd& m_S_P, const MatrixXcd& m_S_S,
	const MatrixXcd& m_S_oct, const MatrixXcd& m_S_L_tet, const MatrixXcd& m_S_R_tet, const MatrixXcd& m_S_hex, Params& p, Input& input)
{
	int M = input.max_mode;
	int size = (p.octamer_L + 4) * M;
	double coeff_coop = std::exp(-p.coop);
	double coeff_offset = std::exp(p.block_offset_energy);
	MatrixXcd L = MatrixXcd::Zero(size, size);
	MatrixXcd I = MatrixXcd::Identity(M, M);

	/***************** Protein part *****************/ 

	MatrixXcd S_B_tilde = I;
	MatrixXcd S_L_tilde = I;
	MatrixXcd S_P_tilde = I;
	MatrixXcd S_S_tilde = I;

	if (p.block_size > 1)
	{
		S_B_tilde = matr_pow(m_S_B, p.block_size - 1);
		S_L_tilde = matr_pow(m_S_L, p.block_size - 1);
		S_P_tilde = matr_pow(m_S_P, p.block_size - 1);
		S_S_tilde = matr_pow(m_S_S, p.block_size - 1);
	}

	L.block(0,     4 * M, M, M) = S_S_tilde * coeff_coop;
	L.block(M,     4 * M, M, M) = S_P_tilde * coeff_coop;
	L.block(2 * M, 4 * M, M, M) = S_L_tilde * coeff_coop;
	L.block(3 * M, 4 * M, M, M) = S_B_tilde;

	for (int i = 4; i <= p.octamer_L + 2; i++) {
		L.block(i * M, (i + 1) * M, M, M) = I * coeff_offset;
	}

	/// Octamer ///  

	double H1_energy_coeff = 1.0 + std::exp(p.H1_mu);

	L.block((p.octamer_L + 3) * M, 0,     M, M) = m_S_oct * coeff_coop * coeff_offset * H1_energy_coeff;
	L.block((p.octamer_L + 3) * M, M,     M, M) = m_S_oct * coeff_coop * coeff_offset * H1_energy_coeff;
	L.block((p.octamer_L + 3) * M, 2 * M, M, M) = m_S_oct * coeff_coop * coeff_offset * H1_energy_coeff;
	L.block((p.octamer_L + 3) * M, 3 * M, M, M) = m_S_oct * coeff_offset * H1_energy_coeff;
 
	/// L- and R-tetrasomes ///  

	MatrixXcd m_S_tet = m_S_L_tet + m_S_R_tet;

	L.block((p.tetrasome_L + 3) * M, 0,     M, M) = m_S_tet * coeff_coop * coeff_offset;
	L.block((p.tetrasome_L + 3) * M, M,     M, M) = m_S_tet * coeff_coop * coeff_offset;
	L.block((p.tetrasome_L + 3) * M, 2 * M, M, M) = m_S_tet * coeff_coop * coeff_offset;
	L.block((p.tetrasome_L + 3) * M, 3 * M, M, M) = m_S_tet * coeff_offset;

	/// Hexamer ///

	L.block((p.hexamer_L + 3) * M, 0,     M, M) = m_S_hex * coeff_coop * coeff_offset;
	L.block((p.hexamer_L + 3) * M, M,     M, M) = m_S_hex * coeff_coop * coeff_offset;
	L.block((p.hexamer_L + 3) * M, 2 * M, M, M) = m_S_hex * coeff_coop * coeff_offset;
	L.block((p.hexamer_L + 3) * M, 3 * M, M, M) = m_S_hex * coeff_offset;

	/***************** Bare DNA part *****************/  

	S_B_tilde *= m_S_B;
	S_L_tilde *= m_S_L;
	S_P_tilde *= m_S_P;
	S_S_tilde *= m_S_S;

	L.block(0, 0,     M, M) = S_S_tilde;
	L.block(0, M,     M, M) = S_S_tilde * coeff_coop;
	L.block(0, 2 * M, M, M) = S_S_tilde * coeff_coop;
	L.block(0, 3 * M, M, M) = S_S_tilde * coeff_coop;

	L.block(M, 0,     M, M) = S_P_tilde * coeff_coop;
	L.block(M, M,     M, M) = S_P_tilde;
	L.block(M, 2 * M, M, M) = S_P_tilde * coeff_coop;
	L.block(M, 3 * M, M, M) = S_P_tilde * coeff_coop;

	L.block(2 * M, 0,     M, M) = S_L_tilde * coeff_coop;
	L.block(2 * M, M,     M, M) = S_L_tilde * coeff_coop;
	L.block(2 * M, 2 * M, M, M) = S_L_tilde;
	L.block(2 * M, 3 * M, M, M) = S_L_tilde * coeff_coop;

	L.block(3 * M, 0,     M, M) = S_B_tilde * coeff_coop;
	L.block(3 * M, M,     M, M) = S_B_tilde * coeff_coop;
	L.block(3 * M, 2 * M, M, M) = S_B_tilde * coeff_coop;
	L.block(3 * M, 3 * M, M, M) = S_B_tilde;

	return L;
}