

function B = matrix_B_mod_Ener_Aug(Aux,C,C1,C2,phi,d_phi)

global Invariant ini_Val Material Constraints dt

if isfield(ini_Val,'iopen')
    iopen = ini_Val.iopen;
else
    iopen = 0;
end

% For the (alpha/2)*log(beta*phi) term added to the energy
if isfield(Material,'alpha')
    alpha = Material.alpha;
else
    alpha = 0;
end

% To assign the bounds in material properties due to proteins
Chi = Material.Chi;
RT_a0 = Material.RT_a0;
C0 = Constraints.C0;

if (alpha >0)
    tmp4 = alpha./phi./phi + alpha./(Material.phi_max-phi).^2;
    tmp5 = -alpha./(phi).^2 + alpha./(Material.phi_max-phi).^2 + 2*alpha*phi./(Material.phi_max-phi).^3;
else
    tmp4=0;
    tmp5=0;
end


% tmp4 = alpha*beta^2*sech(beta*(Material.phi_max-phi)).^2.*tanh(beta*(Material.phi_max-phi));
% tmp5 = alpha*beta^3.*phi.*(2*sech(beta*(Material.phi_max-phi)).^2.*tanh(beta*(Material.phi_max-phi)).^2 ...
%     - sech(beta*(Material.phi_max-phi)).^4);
% 
% tmp6 = -gamma*delta^2.*sech(delta*(phi+0.05)).^2.*tanh(delta*(phi+0.05));
% tmp7 = gamma*delta^3.*phi.*(2*sech(delta*(phi+0.05)).^2.*tanh(delta*(phi+0.05)).^2 + sech(delta*(phi+0.05)).^4); 

B=zeros(ini_Val.ncontrol);
connect = Invariant.connect;
aux =  (C(:,1)./ Aux.au) .* (RT_a0 * Material.phi_max./(Material.phi_max - phi).^2 ...
    + Material.k0*C0.^2 + Chi + tmp4 + tmp5).*d_phi .* Invariant.gweight*dt;

for iel=1:size(connect,1)
    ig=(1:ini_Val.gaussno);
    igglob=(iel-1)*ini_Val.gaussno+ig;
    B_loc=Invariant.Nshape(igglob,:)'*diag(aux(igglob))*Invariant.dNsh(igglob,:);
    B(connect(iel,:),connect(iel,:)) = ...
        B(connect(iel,:),connect(iel,:)) + B_loc;
end

% if (iopen == 3)
%     B(:,end) = 0;
%     B(end,:) = 0;
% end
