Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
def mean_stderr(ts):
"""Return mean and standard deviation of a time series ts."""
return np.mean(ts), \
np.std(ts)*np.sqrt(statisticalInefficiency(ts, warn=False)/len(ts))
def mean_stderr(ts):
""" Get mean and standard deviation of a time series. """
return np.mean(ts), np.std(ts)*np.sqrt(statisticalInefficiency(ts, warn=False)/len(ts))
v_ = kwargs['v_']
b0 = np.ones(L,dtype=float)
dx = d_[:,0]
dy = d_[:,1]
dz = d_[:,2]
D2 = bzavg(dx**2,b)-bzavg(dx,b)**2
D2 += bzavg(dy**2,b)-bzavg(dy,b)**2
D2 += bzavg(dz**2,b)-bzavg(dz,b)**2
return prefactor*D2/bzavg(v_,b)/T
Eps0 = calc_eps0(None,**{'d_':Dips, 'v_':V})
Eps0boot = []
for i in range(numboots):
boot = np.random.randint(L,size=L)
Eps0boot.append(calc_eps0(None,**{'d_':Dips[boot], 'v_':V[boot]}))
Eps0boot = np.array(Eps0boot)
Eps0_err = np.std(Eps0boot)*np.sqrt(np.mean([statisticalInefficiency(Dips[:,0]),statisticalInefficiency(Dips[:,1]),statisticalInefficiency(Dips[:,2])]))
# Dielectric constant analytic derivative
Dx = Dips[:,0]
Dy = Dips[:,1]
Dz = Dips[:,2]
D2 = avg(Dx**2)+avg(Dy**2)+avg(Dz**2)-avg(Dx)**2-avg(Dy)**2-avg(Dz)**2
GD2 = 2*(flat(np.dot(GDx,col(Dx)))/L - avg(Dx)*(np.mean(GDx,axis=1))) - Beta*(covde(Dx**2) - 2*avg(Dx)*covde(Dx))
GD2 += 2*(flat(np.dot(GDy,col(Dy)))/L - avg(Dy)*(np.mean(GDy,axis=1))) - Beta*(covde(Dy**2) - 2*avg(Dy)*covde(Dy))
GD2 += 2*(flat(np.dot(GDz,col(Dz)))/L - avg(Dz)*(np.mean(GDz,axis=1))) - Beta*(covde(Dz**2) - 2*avg(Dz)*covde(Dz))
GEps0 = prefactor*(GD2/avg(V) - mBeta*covde(V)*D2/avg(V)**2)/T
Sep = printcool("Dielectric constant: % .4e +- %.4e\nAnalytic Derivative:" % (Eps0, Eps0_err))
FF.print_map(vals=GEps0)
if FDCheck:
GEps0_fd = property_derivatives(Liquid, FF, mvals, h, pgrad, kT, calc_eps0, {'d_':Dips,'v_':V})
Sep = printcool("Numerical Derivative:")
FF.print_map(vals=GEps0_fd)
#----
# Isothermal compressibility
#----
def calc_kappa(b=None, **kwargs):
if b is None: b = np.ones(L,dtype=float)
if 'v_' in kwargs:
v_ = kwargs['v_']
return bar_unit / kT * (bzavg(v_**2,b)-bzavg(v_,b)**2)/bzavg(v_,b)
Kappa = calc_kappa(None,**{'v_':V})
Kappaboot = []
for i in range(numboots):
boot = np.random.randint(L,size=L)
Kappaboot.append(calc_kappa(None,**{'v_':V[boot]}))
Kappaboot = np.array(Kappaboot)
Kappa_err = np.std(Kappaboot) * np.sqrt(statisticalInefficiency(V))
# Isothermal compressibility analytic derivative
Sep = printcool("Isothermal compressibility: % .4e +- %.4e bar^-1\nAnalytic Derivative:" % (Kappa, Kappa_err))
GKappa1 = +1 * Beta**2 * avg(V**2) * deprod(V) / avg(V)**2
GKappa2 = -1 * Beta**2 * avg(V) * deprod(V**2) / avg(V)**2
GKappa3 = +1 * Beta**2 * covde(V)
GKappa = bar_unit*(GKappa1 + GKappa2 + GKappa3)
FF.print_map(vals=GKappa)
if FDCheck:
GKappa_fd = property_derivatives(Lipid, FF, mvals, h, pgrad, kT, calc_kappa, {'v_':V})
Sep = printcool("Numerical Derivative:")
FF.print_map(vals=GKappa_fd)
Sep = printcool("Difference (Absolute, Fractional):")
absfrac = ["% .4e % .4e" % (i-j, (i-j)/j) for i,j in zip(GKappa, GKappa_fd)]
FF.print_map(vals=absfrac)
v_ = kwargs['v_']
b0 = np.ones(L,dtype=float)
dx = d_[:,0]
dy = d_[:,1]
dz = d_[:,2]
D2 = bzavg(dx**2,b)-bzavg(dx,b)**2
D2 += bzavg(dy**2,b)-bzavg(dy,b)**2
D2 += bzavg(dz**2,b)-bzavg(dz,b)**2
return prefactor*D2/bzavg(v_,b)/T
Eps0 = calc_eps0(None,**{'d_':Dips, 'v_':V})
Eps0boot = []
for i in range(numboots):
boot = np.random.randint(L,size=L)
Eps0boot.append(calc_eps0(None,**{'d_':Dips[boot], 'v_':V[boot]}))
Eps0boot = np.array(Eps0boot)
Eps0_err = np.std(Eps0boot)*np.sqrt(np.mean([statisticalInefficiency(Dips[:,0]),statisticalInefficiency(Dips[:,1]),statisticalInefficiency(Dips[:,2])]))
# Dielectric constant analytic derivative
Dx = Dips[:,0]
Dy = Dips[:,1]
Dz = Dips[:,2]
D2 = avg(Dx**2)+avg(Dy**2)+avg(Dz**2)-avg(Dx)**2-avg(Dy)**2-avg(Dz)**2
GD2 = 2*(flat(np.dot(GDx,col(Dx)))/L - avg(Dx)*(np.mean(GDx,axis=1))) - Beta*(covde(Dx**2) - 2*avg(Dx)*covde(Dx))
GD2 += 2*(flat(np.dot(GDy,col(Dy)))/L - avg(Dy)*(np.mean(GDy,axis=1))) - Beta*(covde(Dy**2) - 2*avg(Dy)*covde(Dy))
GD2 += 2*(flat(np.dot(GDz,col(Dz)))/L - avg(Dz)*(np.mean(GDz,axis=1))) - Beta*(covde(Dz**2) - 2*avg(Dz)*covde(Dz))
GEps0 = prefactor*(GD2/avg(V) - mBeta*covde(V)*D2/avg(V)**2)/T
Sep = printcool("Dielectric constant: % .4e +- %.4e\nAnalytic Derivative:" % (Eps0, Eps0_err))
FF.print_map(vals=GEps0)
if FDCheck:
GEps0_fd = property_derivatives(Lipid, FF, mvals, h, pgrad, kT, calc_eps0, {'d_':Dips,'v_':V})
Sep = printcool("Numerical Derivative:")
FF.print_map(vals=GEps0_fd)
#----
# Isothermal compressibility
#----
def calc_kappa(b=None, **kwargs):
if b is None: b = np.ones(L,dtype=float)
if 'v_' in kwargs:
v_ = kwargs['v_']
return bar_unit / kT * (bzavg(v_**2,b)-bzavg(v_,b)**2)/bzavg(v_,b)
Kappa = calc_kappa(None,**{'v_':V})
Kappaboot = []
for i in range(numboots):
boot = np.random.randint(L,size=L)
Kappaboot.append(calc_kappa(None,**{'v_':V[boot]}))
Kappaboot = np.array(Kappaboot)
Kappa_err = np.std(Kappaboot) * np.sqrt(statisticalInefficiency(V))
# Isothermal compressibility analytic derivative
Sep = printcool("Isothermal compressibility: % .4e +- %.4e bar^-1\nAnalytic Derivative:" % (Kappa, Kappa_err))
GKappa1 = +1 * Beta**2 * avg(V**2) * deprod(V) / avg(V)**2
GKappa2 = -1 * Beta**2 * avg(V) * deprod(V**2) / avg(V)**2
GKappa3 = +1 * Beta**2 * covde(V)
GKappa = bar_unit*(GKappa1 + GKappa2 + GKappa3)
FF.print_map(vals=GKappa)
if FDCheck:
GKappa_fd = property_derivatives(Liquid, FF, mvals, h, pgrad, kT, calc_kappa, {'v_':V})
Sep = printcool("Numerical Derivative:")
FF.print_map(vals=GKappa_fd)
Sep = printcool("Difference (Absolute, Fractional):")
absfrac = ["% .4e % .4e" % (i-j, (i-j)/j) for i,j in zip(GKappa, GKappa_fd)]
FF.print_map(vals=absfrac)
def mean_stderr(ts):
""" Get mean and standard deviation of a time series. """
if ts.ndim == 1:
return np.mean(ts, axis = 0), np.std(ts, axis = 0)*np.sqrt(statisticalInefficiency(ts, warn=False)/len(ts))
else:
return np.mean(ts, axis = 0), np.std(ts, axis = 0)*np.sqrt(multiD_statisticalInefficiency(ts, warn=False)/len(ts))
def mean_stderr(ts):
""" Get mean and standard deviation of a time series. """
return np.mean(ts), np.std(ts)*np.sqrt(statisticalInefficiency(ts, warn=False)/len(ts))