Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
else:
# G0 and H0 are used for determining the expected function change.
G0 = G.copy()
H0 = H.copy()
G = np.delete(G, self.excision)
H = np.delete(H, self.excision, axis=0)
H = np.delete(H, self.excision, axis=1)
logger.debug("Inverting Hessian:\n") ###
logger.debug(" G:\n") ###
pvec1d(G,precision=5, loglevel=DEBUG) ###
logger.debug(" H:\n") ###
pmat2d(H,precision=5, loglevel=DEBUG) ###
Hi = invert_svd(np.mat(H))
dx = flat(-1 * Hi * col(G))
logger.debug(" dx:\n") ###
pvec1d(dx,precision=5, loglevel=DEBUG) ###
# dxa = -solve(H, G) # Take Newton Raphson Step ; use -1*G if want steepest descent.
# dxa = flat(dxa)
# print " dxa:" ###
# pvec1d(dxa,precision=5) ###
logger.info('\n') ###
for i in self.excision: # Reinsert deleted coordinates - don't take a step in those directions
dx = np.insert(dx, i, 0)
def para_solver(L):
# Levenberg-Marquardt
# HT = H + (L-1)**2*np.diag(np.diag(H))
# Attempt to use plain Levenberg
HT = H + (L-1)**2*np.eye(len(H))
def para_solver(L):
# Levenberg-Marquardt
# HT = H + (L-1)**2*np.diag(np.diag(H))
# Attempt to use plain Levenberg
HT = H + (L-1)**2*np.eye(len(H))
logger.debug("Inverting Scaled Hessian:\n")
logger.debug(" G:\n")
pvec1d(G,precision=5, loglevel=DEBUG)
logger.debug(" HT: (Scal = %.4f)\n" % (1+(L-1)**2))
pmat2d(HT,precision=5, loglevel=DEBUG)
Hi = invert_svd(HT)
dx = flat(-1 * np.dot(Hi, col(G)))
logger.debug(" dx:\n")
pvec1d(dx,precision=5, loglevel=DEBUG)
sol = flat(0.5*multi_dot([row(dx), H, col(dx)]))[0] + np.dot(dx,G)
for i in self.excision: # Reinsert deleted coordinates - don't take a step in those directions
dx = np.insert(dx, i, 0)
return dx, sol
# Build the first Hvap derivative.
# We don't pass it back, but nice for printing.
GHvap = np.mean(G,axis=1)
GHvap += mBeta * (flat(np.mat(G) * col(Energies)) / N - Pot_avg * np.mean(G, axis=1))
GHvap /= 216
GHvap -= np.mean(mG,axis=1)
GHvap -= mBeta * (flat(np.mat(mG) * col(mEnergies)) / N - mPot_avg * np.mean(mG, axis=1))
GHvap *= -1
GHvap -= mBeta * (flat(np.mat(G) * col(pV)) / N - np.mean(pV) * np.mean(G, axis=1))
print "pV_avg = ", pV_avg, "pV_err = ", pV_err
print "kT, -pV", kT, -1 * np.mean(pV)
print "Derivative of pV term (note it is SUBTRACTED from enthalpy of vaporization, so the contribution to Hvap is minus one times this):"
FF.print_map(vals = mBeta * (flat(np.mat(G) * col(pV)) / N - np.mean(pV) * np.mean(G, axis=1)))
bar = printcool("Density: % .4f +- % .4f kg/m^3, Derivatives below" % (Rho_avg, Rho_err))
FF.print_map(vals=GRho)
print bar
print "Box energy:", np.mean(Energies)
print "Monomer energy:", np.mean(mEnergies)
bar = printcool("Enthalpy of Vaporization: % .4f +- %.4f kJ/mol, Derivatives below" % (Hvap_avg, Hvap_err))
FF.print_map(vals=GHvap)
print bar
# Print the final force field.
pvals = FF.make(mvals)
with open(os.path.join('npt_result.p'),'w') as f: lp_dump((Rhos, pV, Energies, G, mEnergies, mG, Rho_err, Hvap_err),f)
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)
Sep = printcool("Difference (Absolute, Fractional):")
absfrac = ["% .4e % .4e" % (i-j, (i-j)/j) for i,j in zip(GEps0,GEps0_fd)]
FF.print_map(vals=absfrac)
logger.info("Writing final force field.\n")
pvals = FF.make(mvals)
logger.info("Writing all simulation data to disk.\n")
lp_dump((Rhos, Volumes, Potentials, Energies, Dips, G, [GDx, GDy, GDz], mPotentials, mEnergies, mG, Rho_err, Hvap_err, Alpha_err, Kappa_err, Cp_err, Eps0_err, NMol),'npt_result.p')
dEg = results['Potential_Derivatives']
dEaq = results['Potential_Derivatives'] + results['Hydration_Derivatives']
data[p]['dHyd'] = (flat(np.dot(dEaq,col(expmbH))/L)-np.mean(dEg,axis=1)*np.mean(expmbH)) / np.mean(expmbH)
elif p == "liq":
Eg = results['Potentials'] - results['Hydration']
Eaq = results['Potentials']
# Mean and standard error of the exponentiated hydration energy.
exppbH = np.exp(+1.0*beta*results['Hydration'])
data[p]['Hyd'] = +kT*np.log(np.mean(exppbH))
# Estimate standard error by bootstrap method. We also multiply by the
# square root of the statistical inefficiency of the hydration energy time series.
data[p]['HydErr'] = np.std([+kT*np.log(np.mean(exppbH[np.random.randint(L,size=L)])) for i in range(100)]) * np.sqrt(statisticalInefficiency(results['Hydration']))
if AGrad:
dEg = results['Potential_Derivatives'] - results['Hydration_Derivatives']
dEaq = results['Potential_Derivatives']
data[p]['dHyd'] = -(flat(np.dot(dEg, col(exppbH))/L)-np.mean(dEaq,axis=1)*np.mean(exppbH)) / np.mean(exppbH)
os.chdir('..')
# Calculate the hydration free energy using gas phase, liquid phase or the average of both.
# Note that the molecular dynamics methods return energies in kJ/mol.
if self.hfemode == 'exp_gas':
self.hfe_dict[label] = data['gas']['Hyd'] / 4.184
self.hfe_err[label] = data['gas']['HydErr'] / 4.184
elif self.hfemode == 'exp_liq':
self.hfe_dict[label] = data['liq']['Hyd'] / 4.184
self.hfe_err[label] = data['liq']['HydErr'] / 4.184
elif self.hfemode == 'exp_both':
self.hfe_dict[label] = 0.5*(data['liq']['Hyd']+data['gas']['Hyd']) / 4.184
self.hfe_err[label] = 0.5*(data['liq']['HydErr']+data['gas']['HydErr']) / 4.184
if AGrad:
# Calculate the derivative of the hydration free energy.
if self.hfemode == 'exp_gas':
dD[:, ilabel] = self.whfe[ilabel]*data['gas']['dHyd'] / 4.184
G = zeros(np, dtype=float)
H = zeros((np, np), dtype=float)
for i in range(self.ns):
P = self.whamboltz_wts[i]
Z += P
dVdqP = mat(self.invdists[i])
espqval = espqvals[i]
espmval = dVdqP * col(getqatoms(mvals))
desp = flat(espmval) - espqval
X += P * dot(desp, desp) / self.nesp
D += P * (dot(espqval, espqval) / self.nesp - (sum(espqval) / self.nesp)**2)
if AGrad:
dVdqM = (dVdqP * dqPdqM).T
for p, vsd in ddVdqPdVS.items():
dVdqM[p,:] += flat(vsd[i] * col(getqatoms(mvals)))
G += flat(P * 2 * dVdqM * col(desp)) / self.nesp
if AHess:
d2VdqM2 = zeros(dVdqM.shape, dtype=float)
for p, vsd in dddVdqPdVS2.items():
d2VdqM2[p,:] += flat(vsd[i] * col(getqatoms(mvals)))
H += array(P * 2 * (dVdqM * dVdqM.T + d2VdqM2 * col(desp))) / self.nesp
# Redundant but we keep it anyway
D /= Z
X /= Z
X /= D
G /= Z
G /= D
H /= Z
H /= D
if not in_fd():
self.esp_err = sqrt(X)
# Following is the restraint part
dddVdqPdVS2 = {}
if AGrad:
for p in range(np):
if 'VSITE' in self.FF.plist[p]:
ddVdqPdVS[p], dddVdqPdVS2[p] = f12d3p(fdwrap(self.build_invdist,mvals,p), h = self.h, f0 = self.invdists)
X = 0
D = 0
G = zeros(np, dtype=float)
H = zeros((np, np), dtype=float)
for i in range(self.ns):
P = self.whamboltz_wts[i]
Z += P
dVdqP = mat(self.invdists[i])
espqval = espqvals[i]
espmval = dVdqP * col(getqatoms(mvals))
desp = flat(espmval) - espqval
X += P * dot(desp, desp) / self.nesp
D += P * (dot(espqval, espqval) / self.nesp - (sum(espqval) / self.nesp)**2)
if AGrad:
dVdqM = (dVdqP * dqPdqM).T
for p, vsd in ddVdqPdVS.items():
dVdqM[p,:] += flat(vsd[i] * col(getqatoms(mvals)))
G += flat(P * 2 * dVdqM * col(desp)) / self.nesp
if AHess:
d2VdqM2 = zeros(dVdqM.shape, dtype=float)
for p, vsd in dddVdqPdVS2.items():
d2VdqM2[p,:] += flat(vsd[i] * col(getqatoms(mvals)))
H += array(P * 2 * (dVdqM * dVdqM.T + d2VdqM2 * col(desp))) / self.nesp
# Redundant but we keep it anyway
D /= Z
X /= Z
X /= D
def para_solver(L):
# Levenberg-Marquardt
# HT = H + (L-1)**2*np.diag(np.diag(H))
# Attempt to use plain Levenberg
HT = H + (L-1)**2*np.eye(len(H))
logger.debug("Inverting Scaled Hessian:\n")
logger.debug(" G:\n")
pvec1d(G,precision=5, loglevel=DEBUG)
logger.debug(" HT: (Scal = %.4f)\n" % (1+(L-1)**2))
pmat2d(HT,precision=5, loglevel=DEBUG)
Hi = invert_svd(HT)
dx = flat(-1 * np.dot(Hi, col(G)))
logger.debug(" dx:\n")
pvec1d(dx,precision=5, loglevel=DEBUG)
sol = flat(0.5*multi_dot([row(dx), H, col(dx)]))[0] + np.dot(dx,G)
for i in self.excision: # Reinsert deleted coordinates - don't take a step in those directions
dx = np.insert(dx, i, 0)
return dx, sol
def _compute(self, dx):
self.dx = dx.copy()
Tmp = np.mat(self.H)*col(dx)
Reg_Term = self.Penalty.compute(xkd+flat(dx), Obj0)
self.Val = (X + np.dot(dx, G) + 0.5*row(dx)*Tmp + Reg_Term[0] - data['X'])[0,0]
self.Grad = flat(col(G) + Tmp) + Reg_Term[1]
def compute_val(self, dx):
HT = H + (L-1)**2*np.eye(len(H))
logger.debug("Inverting Scaled Hessian:\n") ###
logger.debug(" G:\n") ###
pvec1d(G,precision=5, loglevel=DEBUG) ###
logger.debug(" HT: (Scal = %.4f)\n" % (1+(L-1)**2)) ###
pmat2d(HT,precision=5, loglevel=DEBUG) ###
Hi = invert_svd(np.mat(HT))
dx = flat(-1 * Hi * col(G))
logger.debug(" dx:\n") ###
pvec1d(dx,precision=5, loglevel=DEBUG) ###
# dxa = -solve(HT, G)
# dxa = flat(dxa)
# print " dxa:" ###
# pvec1d(dxa,precision=5) ###
# print ###
sol = flat(0.5*row(dx)*np.mat(H)*col(dx))[0] + np.dot(dx,G)
for i in self.excision: # Reinsert deleted coordinates - don't take a step in those directions
dx = np.insert(dx, i, 0)
return dx, sol