Demo entry 6169900

Redox Fitting Script

   

Submitted by anonymous on Oct 09, 2016 at 22:05
Language: Python. Code size: 7.0 kB.

##############################################################################
# Script for fitting of redox based experimental
# data of biochemecial relevant molecules
#
#
#
# written by Dominik Lemm, Bsc. 2016


import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
from scipy.optimize import minimize


##############################################################################
#           Global Variables, Conversion Factors, Data Paths


pH_range = np.arange(2, 14, 0.1)  # pH range for the later plots
cpH_conv = -1.364936  # Conversion factor from pH in chemical potential
#mu_h = pH_range * -1.364936  # delete if not needed
F = 23.061        # kcal/mol
R = 1.985877e-3   # kcal/mol
T = 298.5         # Kelvin
beta = 1./(R*T)


ox = "data\\ox.dat"        # Experimental data measured by Anderson or
ox_dp = "data\\ox_dp.dat"  # Draper and Ingraham
red = "data\\red.dat"      # as described in the thesis
red_dp = "data\\red_dp.dat"


#############################################################################
#                       Redox Systems, 6-States and 9-States
#
#
#    States in the form [Proton, Electron]
#    Will later be extended into
#    [Proton, Electron, G]
#    with G being the intrisic free energy
#    during the routine of the cost_function
#
#          6 state model
states = np.array(      [[0, 0],  #  0
                         [1, 0],  #  1
                         [1, 1],  #  2
                         [2, 1],  #  3
                         [2, 2],  #  4
                         [3, 2]]) #  5


#           9 state model
states_9 = np.array(      [[0, 0],  #  0
                         [1, 0],  #  1
                         [2, 0],  #  2
                         [1, 1],  #  3
                         [2, 1],  #  4
                         [3, 1],  #  5
                         [2, 2],  #  6
                         [3, 2],  #  7
                         [4, 2]]) #  8


#############################################################################
#                   Experimental Data Processing
#
#
# Reads the experimental data either of Anderson or Draper and Ingraham
# The data file has to consist of two columns
def read_exp_data(redox_curve):
    try:
        red = np.loadtxt(redox_curve, dtype = float, usecols=([1]))
        ph = np.loadtxt(redox_curve, dtype = float, usecols=([0]))
    except:
        print ("File %s not readable!" % str(redox_curve))
        exit()
    return ph, red


# Converts the experimental data points of the read_exp_data function
# from V in kcal/mol
def calc_delta_G(curve):
    f = read_exp_data(curve)
    delta_G = -F * f[1]
    return f[0], delta_G # in kcal/mol


############################################################################
# Theoretical Model for the Calculation of the Energy G in Depency of the pH
#
#
# Calculates the energy in kcal/mol of the given state in dependency of the pH
# Converts the pH into a chemical potential
# retuern an array of all energies at the given pH values
def calc_one_state(pH,P, E, G):
    G_i = G - P*cpH_conv*pH
    return G_i


# Calculates the partition function of all states with the
# same amount of electrons
# Returns an array of the partition function in dependency of the pH
def calc_state_total(electrons,state_energy, pH):
    if electrons == 0:
        Z = np.sum(np.exp(- beta * np.array([calc_one_state(pH,*state) \
        for state in state_energy[0:2]])), axis=0)
    elif electrons == 1:
        Z = np.sum(np.exp(- beta * np.array([calc_one_state(pH,*state) \
        for state in state_energy[2:4]])), axis=0)
    elif electrons == 2:
        Z = np.sum(np.exp(- beta * np.array([calc_one_state(pH,*state) \
        for state in state_energy[4:6]])), axis=0)
    return Z


# Calculates the transition curve between two redox zones by dividing two
# partition function by each other
# Converts the resulting probability into an energy
# Return the energy in kcal/mol in dependency of the pH
def calc_curve(trans, state_energy, pH):
    if trans == "oxsem":
        G = - 1/beta * np.log(calc_state_total(1,state_energy, pH) / \
        calc_state_total(0,state_energy, pH))
    elif trans == "midpoint":
        G = ((- 1/beta * np.log((calc_state_total(2,state_energy, pH)/ \
        (calc_state_total(0,state_energy, pH)) )))/ 2)
    elif trans == "semred":
        G = - 1/beta * np.log(calc_state_total(2,state_energy, pH) / \
        calc_state_total(1,state_energy, pH))
    return G #in kcal/mol


#############################################################################
#                  Fitting Functions
#
#
# Set up cost functions that calculates the delta between every experimental
# data point and the calculated model
# Converts the states array into the form [proton, electron, G]
# Returns the sum of all deltas in dependency of the pH
# THIS FUNCTION HAS TO BE MINIMIZED IN ORDER TO ALIGN THE POINTS!
def cost_function(vars):
    vars = np.array(vars)
    vars = np.append(0, vars)
    state_energy = np.append(states, vars.reshape((len(vars),1)), axis=1)
    f_ox =calc_delta_G(ox)
    f_red = calc_delta_G(red)
    summe = np.sum((f_ox[1] - calc_curve("oxsem",state_energy, f_ox[0]))**2 +\
    (f_red[1] - calc_curve("semred", state_energy, f_red[0]))**2)
    return summe/float(len(f_ox[0]))


# Minimizing function that tries to find the global minimum
# of the cost function
# Needs slices for each Parameter that has to be fitted in the form of
# slice(-min, max, amount of slices)
# optimize.brute calculates every slice and tries to find minima
# Calls the optimize.minimize function with the lowest calculated parameters
def minimizer():
    ranges = (slice(-20,20,40), slice(-20,20,40),slice(-20,20,40),\
    slice(-20,20,40),slice(-20,20,40))
    x0, fval, _, _ = optimize.brute(cost_function,ranges,full_output=1,\
    disp=True)
    print(x0)
    print(fval)
    res = minimize(cost_function, x0 , method='Powell' ,\
     options={'disp':True , 'maxiter':1000000, 'maxfev':1000000})
    print(res)
    print(res.x)
    result_E = np.append(0, res.x)
    states_E = np.append(states,result_E.reshape((len(result_E),1)), axis=1)
    print(states_E)
    f_ox = calc_curve("oxsem",states_E, pH_range)
    f_red = calc_curve("semred",states_E, pH_range)
    E_1 = - calc_curve("semred",states_E, 7.)/F
    E_2 = - calc_curve("oxsem",states_E, 7.)/F
    print("E_1-Value at pH 7: " + str(E_1))
    print("E_2-Value at pH 7: " + str(E_2))
    ox_may = calc_delta_G(ox)
    red_may = calc_delta_G(red)
    plt.figure(figsize=(10,10))
    plt.plot(red_may[0], red_may[1],"ko", label="Exp. data")
    plt.plot(ox_may[0], ox_may[1],"ko")
    plt.plot(pH_range, f_red, label="Red-Curve-Model")
    plt.plot(pH_range, f_ox, label="Ox-Curve-Model")
    plt.xlabel('pH', fontsize=20)
    plt.ylabel('G [kcal/mol]', fontsize=20)
    plt.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc=3,
           ncol=2, mode="expand", borderaxespad=0.)
    plt.show()

This snippet took 0.01 seconds to highlight.

Back to the Entry List or Home.

Delete this entry (admin only).