diff --git a/Minimization/FitTools.py b/Minimization/FitTools.py index d34cc5b..be325ae 100755 --- a/Minimization/FitTools.py +++ b/Minimization/FitTools.py @@ -9,7 +9,7 @@ import math sys.path.append("/mntdirect/_users/semeraro/python_tools") import FormFactor #import PLUV_POPC_RecBuf -from PLUV import SDP_base_POPC_RecBuf, SDP_POPC_RecBuf +from PLUV import SDP_base_POPC_RecBuf, SDP_POPC_RecBuf, SDP_POPC_RecBuf_LogNormal ########################################################################################### ########################################################################################### @@ -168,7 +168,9 @@ def ChooseFunction (function ): elif function=="SDP_POPC_RecBuf": #intensity = PLUV.SDP_POPC_RecBuf intensity = SDP_POPC_RecBuf - + elif function=="SDP_POPC_RecBuf_LogNormal": + #intensity = PLUV.SDP_POPC_RecBuf + intensity = SDP_POPC_RecBuf_LogNormal elif function=="SDP_base_POPC_RecBuf": intensity = SDP_base_POPC_RecBuf diff --git a/Minimization/PLUV.py b/Minimization/PLUV.py index 64bc7bf..b290fa4 100644 --- a/Minimization/PLUV.py +++ b/Minimization/PLUV.py @@ -108,6 +108,14 @@ CONST_V_EDTA = 0.56430 # (nm^3) def PDF_normal(x, mu, sig) : return np.exp(-(x-mu)**2 / (2*sig**2) ) /( sig*np.sqrt(2*np.pi)) +######################################################### +#@njit(parallel=True) +def PDF_lognormal(x, mu_x, sig_x) : + # https://en.wikipedia.org/wiki/Log-normal_distribution + mu = np.log(mu_x**2 / np.sqrt(mu_x**2 + sig_x**2)) + sig = np.sqrt(np.log(1 + sig_x**2 / mu_x**2)) + return np.exp(-(np.log(x)-mu)**2 / (2*sig**2) ) /( x*sig*np.sqrt(2*np.pi)) + ######################################################### def lipid_volume(T) : return (1-CONST_x_PG) * (CONST_a0_V_POPC +T * CONST_a1_V_POPC) + CONST_x_PG * (CONST_a0_V_POPG + T * CONST_a1_V_POPG) @@ -397,4 +405,139 @@ class SDP_POPC_RecBuf: ################################################################################################################## ################################################################################################################## +class SDP_POPC_RecBuf_LogNormal: + +################## + def __init__(self, q, PAR) : + self.q = q + [self.Norm, self.nv, + self.Rm, self.Z, + self.n_TR, self.d_TR, self.s_TR, + self.d_Chol, self.s_Chol, self.d_PCN, self.s_PCN, self.d_CG, self.s_CG, + self.A_L, self.s_D_C, + self.s_CH2, self.d_CH, self.s_CH, self.s_CH3, + self.r_PCN, self.r_CG, self.r12, self.r32, + self.T, self.V_BW, + self.Con] = PAR + + # [Example 1], fixed parameters: + # Norm 1e5 # Normalization + # n_TR 0.0 # Tris fraction + # d_TR 1.0 # Tris width (nm) + # s_TR 0.29 # Tris position (nm) + # d_CH 0.90 # CH position (nm) + # s_CH 0.305 # CH width (nm) + # r12 0.81 # V_CH/V_CH2 + # T 37 # Temperature (°C) + + # [example 1] fixed + Cw = CONST_p0_Cw + CONST_p1_Cw*self.T + CONST_p2_Cw*self.T**2 + CONST_p3_Cw*self.T**3 + xtris = CONST_ctris / Cw # mole fraction of free TRIS in bulk + xEDTA = CONST_cEDTA / Cw # mole fraction of free EDTA in bulk + + # Volumes + # [example 1] fixed + self.V_L = lipid_volume(self.T) + V_HW = water_volume(self.T) + V_HC = self.V_L - ( (1-CONST_x_PG) * CONST_V_HL_PC + CONST_x_PG * CONST_V_HL_PG ) + + # Calculation of mean D_C + self.D_C = V_HC / self.A_L + + # Quasi-molecular volumes + # [example 1] r12 fixed + V_CH2 = V_HC / ( CONST_n_CH2 + CONST_n_CH*self.r12 + CONST_n_CH3*self.r32 ) # Volume of CH2 groups + V_CH = V_CH2 * self.r12 # Volume of CH groups + V_CH3 = V_CH2 * self.r32 # Volume of CH3 groups + + self.V_CG = CONST_V_HL_PC * self.r_CG # Volume of CG group + self.V_PCN = CONST_V_HL_PC * self.r_PCN # Volume of PCN group + self.V_Chol = CONST_V_HL_PC * (1-self.r_PCN-self.r_CG) # Volume of CholCH3 group + + # CONST + CONST_V_PG1 = CONST_V_HL_PG * 0.16 # Kucerka 2012 + CONST_V_PG2 = CONST_V_HL_PG * ( 1 - 0.51 - 0.16) # Kucerka 2012 + + ############### X-ray scattering lengths (nm) + # [example 1] rho_sol fixed + rho_sol = ( CONST_b_HW + xtris*CONST_b_tris + xEDTA*CONST_b_EDTA ) / V_HW + drho_Chol = ( (1-CONST_x_PG)*CONST_b_Chol/self.V_Chol + CONST_x_PG*CONST_b_PG2/CONST_V_PG2 ) - rho_sol + drho_PCN = ( (1-CONST_x_PG)*CONST_b_PCN/self.V_PCN + CONST_x_PG*CONST_b_PG1/CONST_V_PG1 ) - rho_sol + drho_CG = CONST_b_CG / self.V_CG - rho_sol + drho_TR = CONST_b_tris/ CONST_V_tris - rho_sol + drho_CH = CONST_b_CH / V_CH - rho_sol + drho_CH2 = CONST_b_CH2 / V_CH2 - rho_sol + drho_CH3 = CONST_b_CH3 / V_CH3 - rho_sol + drho_HW = CONST_b_HW / self.V_BW - rho_sol + + ############### D_C polydispersity + N = 21 + HC_array = np.linspace(self.D_C-3*self.s_D_C, self.D_C+3*self.s_D_C, N) + LogNormal = PDF_lognormal(HC_array, self.D_C, self.s_D_C) + + ############### calculating scattering amplitude ----------------------------------------------- + self.Am = np.zeros([HC_array.shape[0],self.q.shape[0]],dtype=float) + c_CH = np.zeros(HC_array.shape[0],dtype=float) + c_CH3 = np.zeros(HC_array.shape[0],dtype=float) + + ############### c-prefactors + c_Chol = ( (1-CONST_x_PG)*self.V_Chol + CONST_x_PG*CONST_V_PG2 ) / self.A_L + c_PCN = ( (1-CONST_x_PG)*self.V_PCN + CONST_x_PG*CONST_V_PG1 ) / self.A_L + c_CG = self.V_CG / self.A_L + c_TR = CONST_V_tris*self.n_TR / self.A_L + + for hc in range(HC_array.shape[0]): + c_CH[hc] = V_CH * CONST_n_CH / (V_HC / HC_array[hc] ) + c_CH3[hc] = V_CH3 * CONST_n_CH3 / (V_HC / HC_array[hc] ) + + for hc in range(HC_array.shape[0]): + + # Adding hydrocarbon-chain envelope + self.Am[hc] += 2 * drho_CH2 *HC_array[hc] * FTreal_erf(self.q, 0, 2*HC_array[hc], self.s_CH2) + # Adding CH and CH3 groups + self.Am[hc] += 2 * (drho_CH - drho_CH2) * c_CH[hc] * FTreal_gauss(self.q, self.d_CH, self.s_CH) + self.Am[hc] += 2 * (drho_CH3 - drho_CH2) * c_CH3[hc] * FTreal_gauss(self.q, 0, self.s_CH3) + + # Adding hydration-water envelope + self.Am[hc] += 4 * drho_HW * ( self.d_CG + self.d_PCN + self.d_Chol + CONST_d_shl) * FTreal_erf(self.q, (HC_array[hc]+(self.d_CG+self.d_PCN+self.d_Chol+CONST_d_shl)/2.), (self.d_CG+self.d_PCN+self.d_Chol+CONST_d_shl), self.s_CH2) + # Adding CG, PCN and CholCH3 groups + self.Am[hc] += 2 * (drho_TR - drho_HW) * c_TR * FTreal_gauss(self.q, (HC_array[hc]+self.d_TR/2.), self.s_TR) + self.Am[hc] += 2 * (drho_CG - drho_HW) * c_CG * FTreal_gauss(self.q, (HC_array[hc]+self.d_CG/2.), self.s_CG) + self.Am[hc] += 2 * (drho_PCN - drho_HW) * c_PCN * FTreal_gauss(self.q, (HC_array[hc]+self.d_CG+self.d_PCN/2.), self.s_PCN) + self.Am[hc] += 2 * (drho_Chol - drho_HW) * c_Chol * FTreal_gauss(self.q, (HC_array[hc]+self.d_CG+self.d_PCN+self.d_Chol/2.), self.s_Chol) + + ############### Ensemble average + self.I = np.zeros(self.q.shape[0], dtype=float) + + for hc in range(HC_array.shape[0]): + if hc==0 or hc==N-1 : self.I+= self.Am[hc]**2 * LogNormal[hc] / 2 + else : self.I+= self.Am[hc]**2 * LogNormal[hc] + + self.I*= 6*self.s_D_C/(N-1) + +################## + def intensity(self): + alp = self.Rm/(self.Z+1) + return ( self.Norm * self.nv*1e-6 ) * self.I * ( 16*np.pi**2*mu4(self.q,self.Z,alp) ) + self.Con*( 0.99*(1./(1+np.exp(-8*(self.q-1.)))) + 0.01 ) + +################## + def negative_water(self): + + self.check = 0 + z_array = np.linspace(0.,4.,81) + + CG = Gauss(z_array, self.V_CG, self.D_C+self.d_CG/2., self.s_CG, self.A_L) + PCN = Gauss(z_array, self.V_PCN, self.D_C+self.d_CG+self.d_PCN/2., self.s_PCN, self.A_L) + Chol = Gauss(z_array, self.V_Chol, self.D_C+self.d_CG+self.d_PCN+self.d_Chol/2., self.s_Chol, self.A_L) + TRIS = Gauss(z_array, self.n_TR*CONST_V_tris, self.D_C+self.d_TR/2., self.s_TR, self.A_L) + BW = Slab(z_array, self.D_C+(self.d_CG+self.d_PCN+self.d_Chol+CONST_d_shl)/2., self.d_CG+self.d_PCN+self.d_Chol+CONST_d_shl, self.s_CH2) - CG - PCN - Chol - TRIS + + for i in(BW) : + if i <-0.001 : self.check+= 1 + + return self.check + +################################################################################################################## +################################################################################################################## + # vim ts=4,sts=4,sw=4 diff --git a/plotting/plot_error.py b/plotting/plot_error.py new file mode 100644 index 0000000..1510f13 --- /dev/null +++ b/plotting/plot_error.py @@ -0,0 +1,23 @@ +import matplotlib.pyplot as plt +import numpy as np +import sys +import os + +def plot_error(dir): + # q | I_data | err | I_fit + data = np.genfromtxt(os.path.join(dir, "POPC-test.fit")) + + plt.loglog(data[:,0], data[:, 1], ".", label="I_data") + plt.loglog(data[:,0], data[:, 3], label="I_fit") + + plt.legend() + + plt.show() + +if __name__ == "__main__": + if len(sys.argv) > 1 and os.path.isdir(sys.argv[1]): + dir = sys.argv[1] + else: + dir = "POPC-test" + + plot_error(dir) diff --git a/plotting/stats.py b/plotting/stats.py index 48d88e2..16762c2 100644 --- a/plotting/stats.py +++ b/plotting/stats.py @@ -121,3 +121,4 @@ if __name__ == "__main__": dir = "POPC-test" print_x2_stats(dir) + std_vs_threshold(dir, 13)