Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
print "Creating Double Precision Simulation for parameter derivatives"
mSim, _ = create_simulation_object(mpdb, mSettings, pbc=False, precision="double")
mG = energy_derivatives(mvals, h, mpdb, FF, mXyzs, mSettings, mSim, None, AGrad)
# pV_avg and mean(pV) are exactly the same.
pV = (pressure * Data['volume'] * AVOGADRO_CONSTANT_NA).value_in_unit(kilojoule_per_mole)
kT = (kB * temperature).value_in_unit(kilojoule_per_mole)
# The enthalpy of vaporization in kJ/mol.
Hvap_avg = mPot_avg - Pot_avg / 216 + kT - np.mean(pV) / 216
Hvap_err = np.sqrt(Pot_err**2 / 216**2 + mPot_err**2 + pV_err**2/216**2)
# 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)) / 216
print "The finite difference step size is:",h
Sep = printcool("Density: % .4f +- % .4f kg/m^3, Analytic Derivative" % (Rho_avg, Rho_err))
FF.print_map(vals=GRho)
print Sep
H = Energies + pV
V = np.array(Volumes)
numboots = 1000
L = len(H)
# Run the simulation for the full system and analyze the results. #
#=================================================================#
Data, Xyzs, Boxes, Rhos, Energies = run_simulation(pdb, direct_kwargs)
# Get statistics from our simulation.
Rho_avg, Rho_err, Pot_avg, Pot_err = analyze(Data)
# Now that we have the coordinates, we can compute the energy derivatives.
G, Hd = energy_derivatives(mvals, h, pdb, FF, Xyzs, direct_kwargs, Boxes)
# The density derivative can be computed using the energy derivative.
N = len(Xyzs)
kB = units.BOLTZMANN_CONSTANT_kB * units.AVOGADRO_CONSTANT_NA
mBeta = (-1 / (temperature * kB)).value_in_unit(units.mole / units.kilojoule)
Beta = (1 / (temperature * kB)).value_in_unit(units.mole / units.kilojoule)
# It becomes necessary to make the gradient squared for the second derivative.
GG = G * G
# Build the first density derivative .
GRho = mBeta * (flat(np.mat(G) * col(Rhos)) / N - np.mean(Rhos) * np.mean(G, axis=1))
# The second density derivative is so inaccurate we might as well not compute it. -_-
# HdRho = mBeta * (flat(np.mat(Hd) * col(Rhos)) / N
# - Beta * flat(np.mat(GG) * col(Rhos)) / N
# + Beta * (flat(np.mat(G) * col(Rhos)) / N) * np.mean(G, axis=1)
# - np.mean(Rhos) * (np.mean(Hd, axis=1)
# - Beta * np.mean(GG, axis=1)
# + Beta * np.mean(G, axis=1) * np.mean(G, axis=1)))
#==============================================#
# Now run the simulation for just the monomer. #
#==============================================#
mpdb = PDBFile('mono.pdb')
mData, mXyzs, _trash, _crap, mEnergies = run_simulation(mpdb, mono_kwargs, pbc=False, Trajectory=False)
# Get statistics from our simulation.
_trash, _crap, mPot_avg, mPot_err = analyze(mData)
# Now that we have the coordinates, we can compute the energy derivatives.
Rho_avg, Rho_err, Pot_avg, Pot_err, pV_avg, pV_err = analyze(Data)
# Now that we have the coordinates, we can compute the energy derivatives.
# First create a double-precision simulation object.
DoublePrecisionDerivatives = True
if DoublePrecisionDerivatives and AGrad:
print "Creating Double Precision Simulation for parameter derivatives"
Sim, _ = create_simulation_object(pdb, Settings, pbc=True, precision="double")
G, GDx, GDy, GDz = energy_dipole_derivatives(mvals, h, pdb, FF, Xyzs, Settings, Sim, Boxes, AGrad)
# The density derivative can be computed using the energy derivative.
N = len(Xyzs)
kB = BOLTZMANN_CONSTANT_kB * AVOGADRO_CONSTANT_NA
T = temperature / kelvin
mBeta = (-1 / (temperature * kB)).value_in_unit(mole / kilojoule)
Beta = (1 / (temperature * kB)).value_in_unit(mole / kilojoule)
# Build the first density derivative .
GRho = mBeta * (flat(np.mat(G) * col(Rhos)) / N - np.mean(Rhos) * np.mean(G, axis=1))
#==============================================#
# Now run the simulation for just the monomer. #
#==============================================#
global timestep, nsteps, niterations, nequiliterations
timestep = 1.0 * femtosecond # timestep for integration
nsteps = 100 # number of steps per data record
nequiliterations = 500 # number of equilibration iterations
niterations = 5000 # number of iterations to collect data for
mpdb = PDBFile('mono.pdb')
mData, mXyzs, _trash, _crap, mEnergies, _nah, _dontneed, mSim = run_simulation(mpdb, mSettings, pbc=False, Trajectory=False)
# Get statistics from our simulation.
_trash, _crap, mPot_avg, mPot_err, __trash, __crap = analyze(mData)
# Now that we have the coordinates, we can compute the energy derivatives.
if DoublePrecisionDerivatives and AGrad:
h_ = kwargs['h_']
if 'v_' in kwargs:
v_ = kwargs['v_']
return 1/(kT*T) * (bzavg(h_*v_,b)-bzavg(h_,b)*bzavg(v_,b))/bzavg(v_,b)
Alpha = calc_alpha(None, **{'h_':H, 'v_':V})
Alphaboot = []
for i in range(numboots):
boot = np.random.randint(L,size=L)
Alphaboot.append(calc_alpha(None, **{'h_':H[boot], 'v_':V[boot]}))
Alphaboot = np.array(Alphaboot)
Alpha_err = np.std(Alphaboot) * max([np.sqrt(statisticalInefficiency(V)),np.sqrt(statisticalInefficiency(H))])
## Thermal expansion coefficient analytic derivative
GAlpha1 = mBeta * covde(H*V) / avg(V)
GAlpha2 = Beta * avg(H*V) * covde(V) / avg(V)**2
GAlpha3 = flat(np.mat(G)*col(V))/N/avg(V) - Gbar
GAlpha4 = Beta * covde(H)
GAlpha = (GAlpha1 + GAlpha2 + GAlpha3 + GAlpha4)/(kT*T)
Sep = printcool("Thermal expansion coefficient: % .4e +- %.4e K^-1\nAnalytic Derivative:" % (Alpha, Alpha_err))
FF.print_map(vals=GAlpha)
if FDCheck:
GAlpha_fd = property_derivatives(mvals, h, pdb, FF, Xyzs, Settings, Sim, kT, calc_alpha, {'h_':H,'v_':V}, Boxes)
Sep = printcool("Numerical Derivative:")
FF.print_map(vals=GAlpha_fd)
Sep = printcool("Difference (Absolute, Fractional):")
absfrac = ["% .4e % .4e" % (i-j, (i-j)/j) for i,j in zip(GAlpha, GAlpha_fd)]
FF.print_map(vals=absfrac)
## Isothermal compressibility
bar_unit = 1.0*bar*nanometer**3/kilojoules_per_mole/item
def calc_kappa(b=None, **kwargs):
if b == None: b = np.ones(L,dtype=float)
dVdqP = np.matrix(self.invdists[i])
espqval = espqvals[i]
espmval = dVdqP * col(new_charges(mvals))
desp = flat(espmval) - espqval
X += P * np.dot(desp, desp) / self.nesp
Q += P * np.dot(espqval, espqval) / self.nesp
D += P * (np.dot(espqval, espqval) / self.nesp - (np.sum(espqval) / self.nesp)**2)
if AGrad:
dVdqM = (dVdqP * dqPdqM).T
for p, vsd in ddVdqPdVS.items():
dVdqM[p,:] += flat(vsd[i] * col(new_charges(mvals)))
G += flat(P * 2 * dVdqM * col(desp)) / self.nesp
if AHess:
d2VdqM2 = np.zeros(dVdqM.shape)
for p, vsd in dddVdqPdVS2.items():
d2VdqM2[p,:] += flat(vsd[i] * col(new_charges(mvals)))
H += np.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
Q /= Z
Q /= D
if not in_fd():
self.esp_err = np.sqrt(X)
self.esp_ref = np.sqrt(Q)
self.esp_err_pct = self.esp_err / self.esp_ref
H /= D
if not in_fd():
self.esp_err = sqrt(X)
# Following is the restraint part
# RESP hyperbola "strength" parameter; 0.0005 is weak, 0.001 is strong
# RESP hyperbola "tightness" parameter; don't need to change this
a = self.resp_a
b = self.resp_b
q = getqatoms(mvals)
R = a*sum((q**2 + b**2)**0.5 - b)
dR = a*q*(q**2 + b**2)**-0.5
ddR = a*b**2*(q**2 + b**2)**-1.5
self.respterm = R
X += R
if AGrad:
G += flat(dqPdqM.T * col(dR))
if AHess:
H += diag(flat(dqPdqM.T * col(ddR)))
Answer = {'X':X,'G':G,'H':H}
return Answer
if self.backup:
for fnm in self.FF.fnms:
if os.path.exists(os.path.join(self.resdir, fnm)):
bak(os.path.join(self.resdir, fnm))
self.FF.make(xk,printdir=self.resdir)
# logger.info("The force field has been written to the '%s' directory.\n" % self.resdir)
outfnm = self.save_mvals_to_input(xk)
# logger.info("Input file with optimization parameters saved to %s.\n" % outfnm)
printcool("The force field has been written to the %s directory.\n"
"Input file with optimization parameters saved to %s." % (self.resdir, outfnm), color=0)
#================================#
#| Hessian update for BFGS. |#
#================================#
if b_BFGS:
Hnew = H_stor.copy()
Dx = col(xk - xk_prev)
Dy = col(G - G_prev)
Mat1 = (np.dot(Dy,Dy.T))/(np.dot(Dy.T,Dx))[0,0]
Mat2 = (np.dot(np.dot(Hnew,Dx),np.dot(Hnew,Dx).T))/(multi_dot([Dx.T,Hnew,Dx]))[0,0]
Hnew += Mat1-Mat2
H = Hnew.copy()
data['H'] = H.copy()
# (Experimental): Deleting lines in the parameter file
if len(self.FF.prmdestroy_this) > 0:
self.FF.prmdestroy_save.append(self.FF.prmdestroy_this)
self.FF.linedestroy_save.append(self.FF.linedestroy_this)
# Update objective function history.
X_hist = np.append(X_hist, X)
# Take the stdev over the previous (hist) values.
# Multiply by 2, so when hist=2 this is simply the difference.
stdfront = np.std(X_hist[-self.hist:]) if len(X_hist) > self.hist else np.std(X_hist)
stdfront *= 2
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, Eps0_err = bootstats(calc_eps0,{'d_':Dips, 'v_':V})
Eps0 += 1.0
Eps0_err *= 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.mat(GDx)*col(Dx))/N - avg(Dx)*(np.mean(GDx,axis=1))) - Beta*(covde(Dx**2) - 2*avg(Dx)*covde(Dx))
GD2 += 2*(flat(np.mat(GDy)*col(Dy))/N - avg(Dy)*(np.mean(GDy,axis=1))) - Beta*(covde(Dy**2) - 2*avg(Dy)*covde(Dy))
GD2 += 2*(flat(np.mat(GDz)*col(Dz))/N - 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(mvals, h, FF, args.liquid_xyzfile, args.liquid_keyfile, 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)
## Print the final force field.
pvals = FF.make(mvals)