From 7fb36a99dde391fcfa1a79f76a47483fe1de6cf2 Mon Sep 17 00:00:00 2001 From: Gaspard Jankowiak Date: Wed, 13 Mar 2024 12:07:46 +0100 Subject: [PATCH] fix integration bounds for the polydispersity --- Minimization/FitTools.py | 3 +- Minimization/PLUV.py | 140 ++++++--------------------------------- 2 files changed, 24 insertions(+), 119 deletions(-) diff --git a/Minimization/FitTools.py b/Minimization/FitTools.py index be325ae..7015a80 100755 --- a/Minimization/FitTools.py +++ b/Minimization/FitTools.py @@ -9,7 +9,8 @@ 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, SDP_POPC_RecBuf_LogNormal +from PLUV import SDP_POPC_RecBuf +#from PLUV import SDP_base_POPC_RecBuf, SDP_POPC_RecBuf, SDP_POPC_RecBuf_LogNormal ########################################################################################### ########################################################################################### diff --git a/Minimization/PLUV.py b/Minimization/PLUV.py index b290fa4..f289e7e 100644 --- a/Minimization/PLUV.py +++ b/Minimization/PLUV.py @@ -167,109 +167,6 @@ def mu4(q, Z, a) : ######### liposomes and proteoliposomes ################# ######################################################### -class SDP_base_POPC_RecBuf: - -################## - 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_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 - - 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 - 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 ) - - # Quasi-molecular volumes - 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 - - V_PG1 = CONST_V_HL_PG * 0.16 # Kucerka 2012 - V_PG2 = CONST_V_HL_PG * ( 1 - 0.51 - 0.16) # Kucerka 2012 - - # Calculation of mean D_C - self.D_C = V_HC / self.A_L - - # X-ray scattering lengths (nm) - 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/V_PG2 ) - rho_sol - drho_PCN = ( (1-CONST_x_PG)*CONST_b_PCN/self.V_PCN + CONST_x_PG*CONST_b_PG1/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 - - # c-prefactors - c_Chol = ( (1-CONST_x_PG)*self.V_Chol + CONST_x_PG*V_PG2 ) / self.A_L - c_PCN = ( (1-CONST_x_PG)*self.V_PCN + CONST_x_PG*V_PG1 ) / self.A_L - c_CG = self.V_CG / self.A_L - c_TR = CONST_V_tris*self.n_TR / self.A_L - c_CH = V_CH * CONST_n_CH / self.A_L - c_CH3 = V_CH3 * CONST_n_CH3 / self.A_L - - # calculating scattering amplitude - self.Am=np.zeros(q.shape[0],dtype=float) - - # Adding hydrocarbon-chain envelope - self.Am += 2 * drho_CH2 *self.D_C * FTreal_erf(self.q, 0, 2*self.D_C, self.s_CH2) - # Adding CH and CH3 groups - self.Am += 2 * (drho_CH - drho_CH2) * c_CH * FTreal_gauss(self.q, self.d_CH, self.s_CH) - self.Am += 2 * (drho_CH3 - drho_CH2) * c_CH3 * FTreal_gauss(self.q, 0, self.s_CH3) - - # Adding hydration-water envelope - self.Am += 4 * drho_HW * ( self.d_CG + self.d_PCN + self.d_Chol + CONST_d_shl) * FTreal_erf(self.q, (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) - # Adding CG, PCN and CholCH3 groups - self.Am += 2 * (drho_TR - drho_HW) * c_TR * FTreal_gauss(self.q, (self.D_C+self.d_TR/2.), self.s_TR) - self.Am += 2 * (drho_CG - drho_HW) * c_CG * FTreal_gauss(self.q, (self.D_C+self.d_CG/2.), self.s_CG) - self.Am += 2 * (drho_PCN - drho_HW) * c_PCN * FTreal_gauss(self.q, (self.D_C+self.d_CG+self.d_PCN/2.), self.s_PCN) - self.Am += 2 * (drho_Chol - drho_HW) * c_Chol * FTreal_gauss(self.q, (self.D_C+self.d_CG+self.d_PCN+self.d_Chol/2.), self.s_Chol) - -################## - def amplitude(self): - return self.Am - -################## - def intensity(self): - alp = self.Rm/(self.Z+1) - return ( self.Norm * self.nv*1e-6 ) * self.Am**2 * ( 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 -################################################################################################################## -################################################################################################################## - class SDP_POPC_RecBuf: ################## @@ -336,14 +233,25 @@ class SDP_POPC_RecBuf: 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) + # Number of integration points + N = 201 + + # Parameter dependent integration bounds + # HC_min, HC_max = self.D_C-3*self.s_D_C, self.D_C+3*self.s_D_C + # + # Hardcoded integration bounds for POPC, fixed T = 37, A_L in [0.598, 0.719] + HC_min, HC_max = V_HC/0.719-0.5, V_HC/0.589+0.5 + # Minimum and maximum samples (offset by 1/2 of the integration step w.r.t integration bounds for the midpoint rule) + integration_step = (HC_max - HC_min) / N + HC_first, HC_last = HC_min + 0.5*integration_step, HC_max - 0.5*integration_step + + HC_array = np.linspace(HC_first, HC_last, N) Normal = PDF_normal(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) + self.Am = np.zeros([N,self.q.shape[0]],dtype=float) + c_CH = np.zeros(N,dtype=float) + c_CH3 = np.zeros(N,dtype=float) ############### c-prefactors c_Chol = ( (1-CONST_x_PG)*self.V_Chol + CONST_x_PG*CONST_V_PG2 ) / self.A_L @@ -351,11 +259,11 @@ class SDP_POPC_RecBuf: 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]): + for hc in range(N): 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]): + for hc in range(N): # 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) @@ -372,13 +280,9 @@ class SDP_POPC_RecBuf: 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 * Normal[hc] / 2 - else : self.I+= self.Am[hc]**2 * Normal[hc] - - self.I*= 6*self.s_D_C/(N-1) + # multiply each columns of Am by Normal and sum along the columns, + # then multiply by integration step + self.I = np.einsum('ij,i->j', self.Am**2, Normal) * integration_step ################## def intensity(self): @@ -471,7 +375,7 @@ class SDP_POPC_RecBuf_LogNormal: drho_HW = CONST_b_HW / self.V_BW - rho_sol ############### D_C polydispersity - N = 21 + N = 200 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)