diff --git a/src/abinitio_pairwise.py b/src/abinitio_pairwise.py new file mode 100644 index 000000000..5231e80e8 --- /dev/null +++ b/src/abinitio_pairwise.py @@ -0,0 +1,368 @@ +""" @package forcebalance.abinitio_pairwise Ab-initio fitting module with pairwise energies. + +@author Chapin Cavender +@date 05/2024 +""" +from __future__ import division +from __future__ import print_function + +from builtins import zip +from builtins import range +import os +import shutil +from forcebalance.abinitio import AbInitio, plot_qm_vs_mm +from forcebalance.nifty import col, eqcgmx, flat, floatornan, fqcgmx, invert_svd, kb, printcool, bohr2ang, warn_press_key, warn_once, pvec1d, commadash, uncommadash, isint +import numpy as np +from forcebalance.target import Target +from forcebalance.molecule import Molecule, format_xyz_coord +from re import match, sub +import subprocess +from subprocess import PIPE +from forcebalance.finite_difference import fdwrap, f1d2p, f12d3p, in_fd +from collections import defaultdict, OrderedDict +import itertools + +from forcebalance.output import getLogger +logger = getLogger(__name__) + +class AbInitioPairwise(AbInitio): + + """ Subclass of AbInitio with objective function composed of pairwise energy differences.""" + + def __init__(self,options,tgt_opts,forcefield): + """ + Initialization; set up pairwise references energies and weights. + """ + + # Raise an error before base class init if force/nft fitting was requested, + # because AbInitio.__init__ resets self.force=0 when no force data exists, + # which would silently swallow the user's (invalid) request. + if tgt_opts.get('force', False) or tgt_opts.get('w_netforce', 0.0) > 0 or tgt_opts.get('w_torque', 0.0) > 0: + raise RuntimeError( + "Fitting using forces, net forces, or net torques is not " + "implemented for this target." + ) + + # Initialize the base class + super(AbInitioPairwise,self).__init__(options,tgt_opts,forcefield) + + # Also guard against use_nft set by the base class (e.g. via w_netforce/w_torque defaults) + if self.use_nft: + raise RuntimeError( + "Fitting using forces, net forces, or net torques is not " + "implemented for this target." + ) + + # Set up pairwise reference energies and weights + N_pairs = int(self.ns * (self.ns - 1) / 2) + self.eqm_pairs = np.zeros(N_pairs) + self.boltz_wt_pairs = np.ones(N_pairs) + + for i, pair in enumerate(itertools.combinations(range(self.ns), 2)): + self.eqm_pairs[i] = self.eqm[pair[0]] - self.eqm[pair[1]] + self.boltz_wt_pairs[i] = np.sqrt( + self.boltz_wts[pair[0]] * self.boltz_wts[pair[1]] + ) + + self.boltz_wt_pairs /= self.boltz_wt_pairs.sum() + + + def get_energy_force(self, mvals, AGrad=False, AHess=False): + """ + This code computes the least squares objective function for pairwise + energy differences. The numerator is a weighted sum of squared + differences between (E_MM,i - E_MM,j) and (E_QM,i - E_QM,j) for all + pairs of grid points (i, j). The weighted variance of the QM energies + -^2 is the denominator. The indicators are set to the + square roots of the numerator and denominator above. + + Fitting using forces, net forces, or net torques is not implemented for + this target. + + In equation form, the objective function is given by: + + \[ = {\frac{{\left( {\sum\limits_{i,j} {{w_i}{{\left( + {\left( E_i^{MM} - E_j^{MM} \right) - \left( E_i^{QM} - E_j^{QM} \right)} + \right)}^2}} } \right)}} {{\sum\limits_{i,j} {{w_i}{{\left( + {E_i^{QM} - {{\bar E}^{QM}}} \right)}^2}} }}} \] + + @param[in] mvals Mathematical parameter values + @param[in] AGrad Switch to turn on analytic gradient + @param[in] AHess Switch to turn on analytic Hessian + @return Answer Contribution to the objective function + """ + Answer = {} + # Create the new force field!! + pvals = self.FF.make(mvals) + + # Number of atoms containing forces being fitted + nat = len(self.fitatoms) + # Number of net forces on molecules + nnf = self.nnf + # Number of torques on molecules + ntq = self.ntq + # Basically the size of the matrix + NCP1 = 3*nat+1 + NParts = 1 + NP = self.FF.np + NS = self.ns + #==============================================================# + # STEP 1: Form all of the arrays. # + #==============================================================# + if (self.w_energy == 0.0 and self.w_force == 0.0 and self.w_netforce == 0.0 and self.w_torque == 0.0): + AGrad = False + AHess = False + # Sum of all the weights + Z = 0.0 + # All vectors with NCP1 elements are ordered as + # [E F_1x F_1y F_1z F_2x ... NF_1x NF_1y ... TQ_1x TQ_1y ... ] + # Vector of QM-quantities + Q = np.zeros(NCP1) + # Mean quantities over the trajectory + M0 = np.zeros(NCP1) + Q0 = np.zeros(NCP1) + X0 = np.zeros(NCP1) + # The mean squared QM-quantities + QQ0 = np.zeros(NCP1) + # Derivatives of the MM-quantity + M_p = np.zeros((NP,NCP1)) + M_pp = np.zeros((NP,NCP1)) + # Means of gradients + M0_p = np.zeros((NP,NCP1)) + M0_pp = np.zeros((NP,NCP1)) + # Vector of objective function terms + SPX = np.zeros(NCP1) + if AGrad: + SPX_p = np.zeros((NP,NCP1)) + # Derivatives of "parts" of objective functions - i.e. + # the sum is taken over the components of force, net force, torque + # but the components haven't been weighted and summed. + X2_Parts_p = np.zeros((NP,NParts)) + if AHess: + SPX_pq = np.zeros((NP,NP,NCP1)) + X2_Parts_pq = np.zeros((NP,NP,NParts)) + # Storage of the MM-quantities and derivatives for all snapshots. + # This saves time because we don't need to execute the external program + # once per snapshot, but requires memory. + M_all = np.zeros((NS,NCP1)) + if AGrad and self.all_at_once: + dM_all = np.zeros((NS,NP,NCP1)) + ddM_all = np.zeros((NS,NP,NCP1)) + #==============================================================# + # STEP 2: Loop through the snapshots. # + #==============================================================# + if self.all_at_once: + logger.debug("\rExecuting\r") + M_all = self.energy_force_transform() + if self.energy_mode == 'qm_minimum': + M_all[:, 0] -= M_all[self.smin, 0] + if AGrad or AHess: + def callM(mvals_): + logger.debug("\r") + pvals = self.FF.make(mvals_) + return self.energy_force_transform() + for p in self.pgrad: + dM_all[:,p,:], ddM_all[:,p,:] = f12d3p(fdwrap(callM, mvals, p), h = self.h, f0 = M_all) + if self.energy_mode == 'qm_minimum': + dM_all[:, p, 0] -= dM_all[self.smin, p, 0] + ddM_all[:, p, 0] -= ddM_all[self.smin, p, 0] + else: + for i in range(NS): + if i % 100 == 0: + logger.debug("Shot %i\r" % i) + M = self.energy_force_transform_one(i) + M_all[i,:] = M.copy() + + if not AGrad: continue + for p in self.pgrad: + def callM(mvals_): + if i % 100 == 0: + logger.debug("\r") + pvals = self.FF.make(mvals_) + return self.energy_force_transform_one(i) + dM_all[i, p, :], ddM_all[i, p, :] = f12d3p(fdwrap(callM, mvals, p), h = self.h, f0 = M) + for i, pair in enumerate(itertools.combinations(range(NS), 2)): + if i % 100 == 0: + logger.debug("\rIncrementing quantities for pair %i, %i\r" % pair) + # Build Boltzmann weights and increment partition function. + P = self.boltz_wt_pairs[i] + Z += P + # Load reference (QM) data + Q[0] = self.eqm_pairs[i] + QQ = Q*Q + # Load MM quantities from M_all array + M = M_all[pair[0]] - M_all[pair[1]] + # MM pair - QM pair difference + X = M-Q + # For asymmetric fit, MM pair energy differences with opposite sign + # as the corresponding QM pair energy difference are given a boost + # factor. + boost = self.energy_asymmetry if ((M[0] < 0.0) != (Q[0] < 0.0)) else 1.0 + # Increment the average quantities + # The [0] indicates that we are fitting the RMS force and not the RMSD + # (without the covariance, subtracting a mean force doesn't make sense.) + # The rest of the array is empty. + M0[0] += P*M[0] + Q0[0] += P*Q[0] + X0[0] += P*X[0] + # We store all elements of the mean-squared QM quantities. + QQ0 += P*QQ + # Increment the objective function. + Xi = X**2 + Xi[0] *= boost + # SPX contains the sum over snapshots + SPX += P * Xi + #==============================================================# + # STEP 2a: Increment gradients and mean quantities. # + #==============================================================# + for p in self.pgrad: + if not AGrad: continue + M_p[p] = dM_all[pair[0], p] - dM_all[pair[1], p] + M_pp[p] = ddM_all[pair[0], p] - ddM_all[pair[1], p] + if all(M_p[p] == 0): continue + M0_p[p][0] += P * M_p[p][0] + Xi_p = 2 * X * M_p[p] + Xi_p[0] *= boost + SPX_p[p] += P * Xi_p + if not AHess: continue + # This formula is more correct, but perhapsively convergence is slower. + #Xi_pq = 2 * (M_p[p] * M_p[p] + X * M_pp[p]) + # Gauss-Newton formula for approximate Hessian + Xi_pq = 2 * (M_p[p] * M_p[p]) + Xi_pq[0] *= boost + SPX_pq[p,p] += P * Xi_pq + for q in range(p): + if all(M_p[q] == 0): continue + if q not in self.pgrad: continue + Xi_pq = 2 * M_p[p] * M_p[q] + Xi_pq[0] *= boost + SPX_pq[p,q] += P * Xi_pq + + #==============================================================# + # STEP 2b: Write energies and forces to disk. # + #==============================================================# + M_all_print = M_all.copy() + Q_all_print = col(self.eqm) + if self.energy_mode == 'average': + M_all_print[:,0] -= np.average(M_all_print[:,0], weights=self.boltz_wts) + Q_all_print[:,0] -= np.average(Q_all_print[:,0], weights=self.boltz_wts) + elif self.energy_mode == 'qm_minimum': + M_all_print[:,0] -= M_all_print[self.smin,0] + Q_all_print[:,0] -= Q_all_print[self.smin,0] + if self.writelevel > 1: + np.savetxt('M.txt',M_all_print) + np.savetxt('Q.txt',Q_all_print) + if self.writelevel > 0: + EnergyComparison = np.hstack((col(np.arange(NS)), + col(Q_all_print[:,0]), + col(M_all_print[:,0]), + col(M_all_print[:,0])-col(Q_all_print[:,0]), + col(self.boltz_wts))) + np.savetxt("EnergyCompare.txt", EnergyComparison, fmt=" %12i % 12.5f % 12.5f % 13.5f % 12.5e", header="%11s %12s %12s %13s %12s" % ("Num", "QMEnergy", "MMEnergy", "DeltaE(MM-QM)", "Weight")) + if self.writelevel > 1: + plot_qm_vs_mm(Q_all_print[:,0], M_all_print[:,0], + M_orig=self.M_orig[:,0] if self.M_orig is not None else None, + title='Abinitio '+self.name) + if self.M_orig is None: + self.M_orig = M_all_print.copy() + + #==============================================================# + # STEP 3: Build the objective function. # + #==============================================================# + + logger.debug("Done with snapshots, building objective function now\r") + W_Components = np.zeros(NParts) + W_Components[0] = self.w_energy + if np.sum(W_Components) > 0 and self.w_normalize: + W_Components /= np.sum(W_Components) + + if self.energy_rms_override != 0.0: + QQ0[0] = self.energy_rms_override ** 2 + Q0[0] = 0.0 + + def compute_objective(SPX_like,divide=1,L=None,R=None,L2=None,R2=None): + a = 0 + n = 1 + X2E = compute_objective_part(SPX_like,QQ0,Q0,Z,a,n,energy=False,subtract_mean=(self.energy_mode=='average'), + divide=divide,L=L,R=R,L2=L2,R2=R2) + objs = [X2E] + return np.array(objs) + + # The objective function components (i.e. energy, force, net force, torque) + X2_Components = compute_objective(SPX,L=X0,R=X0) + # The squared residuals before they are normalized + X2_Physical = compute_objective(SPX,divide=0,L=X0,R=X0) + # The normalization factors + X2_Normalize = compute_objective(SPX,divide=2,L=X0,R=X0) + # The derivatives of the objective function components + for p in self.pgrad: + if not AGrad: continue + X2_Parts_p[p,:] = compute_objective(SPX_p[p],L=2*X0,R=M0_p[p]) + if not AHess: continue + X2_Parts_pq[p,p,:] = compute_objective(SPX_pq[p,p],L=2*M0_p[p],R=M0_p[p],L2=2*X0,R2=M0_pp[p]) + for q in range(p): + if q not in self.pgrad: continue + X2_Parts_pq[p,q,:] = compute_objective(SPX_pq[p,q],L=2*M0_p[p],R=M0_p[q]) + # Get the other half of the Hessian matrix. + X2_Parts_pq[q,p,:] = X2_Parts_pq[p,q,:] + # The objective function as a weighted sum of the components + X2 = np.dot(X2_Components, W_Components) + # Derivatives of the objective function + G = np.zeros(NP) + H = np.zeros((NP,NP)) + for p in self.pgrad: + if not AGrad: continue + G[p] = np.dot(X2_Parts_p[p], W_Components) + if not AHess: continue + for q in self.pgrad: + H[p,q] = np.dot(X2_Parts_pq[p,q], W_Components) + + #==============================================================# + # STEP 3a: Build the indicators. # + #==============================================================# + # Save values to qualitative indicator if not inside finite difference code. + if not in_fd(): + # Contribution from energy and force parts. + tw = self.w_energy + self.w_force + self.w_netforce + self.w_torque + self.w_resp + self.e_trm = X2_Components[0] + self.e_ctr = X2_Components[0]*W_Components[0] + if self.w_normalize: self.e_ctr /= tw + self.e_ref = np.sqrt(X2_Normalize[0]) + self.e_err = np.sqrt(X2_Physical[0]) + self.e_err_pct = self.e_err/self.e_ref + + pvals = self.FF.make(mvals) # Write a force field that isn't perturbed by finite differences. + Answer = {'X':X2, 'G':G, 'H':H} + return Answer + + +def compute_objective_part(SPX,QQ0,Q0,Z,a,n,energy=False,subtract_mean=False,divide=1,L=None,R=None,L2=None,R2=None): + # Divide by Z to normalize + XiZ = SPX[a:a+n]/Z + QQ0iZ = QQ0[a:a+n]/Z + Q0iZ = Q0[a:a+n]/Z + # Subtract out the product of L and R if provided. + if subtract_mean: + if L is not None and R is not None: + LiZ = L[a:a+n]/Z + RiZ = R[a:a+n]/Z + XiZ -= LiZ*RiZ + elif L2 is not None and R2 is not None: + L2iZ = L2[a:a+n]/Z + R2iZ = R2[a:a+n]/Z + XiZ -= L2iZ*R2iZ + else: + raise RuntimeError("subtract_mean is set to True, must provide L/R or L2/R2") + if energy: + QQ0iZ -= Q0iZ*Q0iZ + + # Return the answer. + if divide == 1: + X2 = np.sum(XiZ)/np.sum(QQ0iZ) + elif divide == 0: + X2 = np.sum(XiZ) + elif divide == 2: + X2 = np.sum(QQ0iZ) + else: + raise RuntimeError('Please pass either 0, 1, 2 to divide') + return X2 diff --git a/src/objective.py b/src/objective.py index ca9da81b0..6c1cda7a9 100644 --- a/src/objective.py +++ b/src/objective.py @@ -37,7 +37,7 @@ logger.warning("OpenMM module import failed; check OpenMM package\n") try: - from forcebalance.smirnoffio import AbInitio_SMIRNOFF, Liquid_SMIRNOFF, Vibration_SMIRNOFF, Hessian_SMIRNOFF, OptGeoTarget_SMIRNOFF, TorsionProfileTarget_SMIRNOFF, smirnoff_analyze_parameter_coverage + from forcebalance.smirnoffio import AbInitio_SMIRNOFF, AbInitioPairwise_SMIRNOFF, Liquid_SMIRNOFF, Vibration_SMIRNOFF, Hessian_SMIRNOFF, OptGeoTarget_SMIRNOFF, TorsionProfileTarget_SMIRNOFF, smirnoff_analyze_parameter_coverage except: logger.warning(traceback.format_exc()) logger.warning("SMIRNOFF module import failed; check SMIRNOFF package\n") @@ -90,6 +90,7 @@ 'ABINITIO_TINKER':AbInitio_TINKER, 'ABINITIO_OPENMM':AbInitio_OpenMM, 'ABINITIO_SMIRNOFF':AbInitio_SMIRNOFF, + 'ABINITIOPAIRWISE_SMIRNOFF':AbInitioPairwise_SMIRNOFF, 'ABINITIO_AMBER':AbInitio_AMBER, 'ABINITIO_INTERNAL':AbInitio_Internal, 'VIBRATION_TINKER':Vibration_TINKER, diff --git a/src/smirnoffio.py b/src/smirnoffio.py index f167c4a7e..012775e8f 100644 --- a/src/smirnoffio.py +++ b/src/smirnoffio.py @@ -10,6 +10,7 @@ import os from forcebalance import BaseReader from forcebalance.abinitio import AbInitio +from forcebalance.abinitio_pairwise import AbInitioPairwise from forcebalance.binding import BindingEnergy from forcebalance.liquid import Liquid from forcebalance.interaction import Interaction @@ -833,6 +834,24 @@ def submit_jobs(self, mvals, AGrad=False, AHess=False): # we update the self.pgrads here so it's not overwritten in rtarget.py smirnoff_update_pgrads(self) +class AbInitioPairwise_SMIRNOFF(AbInitioPairwise): + """ Pairwise energy matching using OpenMM. """ + def __init__(self,options,tgt_opts,forcefield): + ## Default file names for coordinates and key file. + self.set_option(tgt_opts,'pdb',default="conf.pdb") + # List of .mol2 files for SMIRNOFF to set up the system + self.set_option(tgt_opts,'mol2',forceprint=True) + self.set_option(tgt_opts,'coords',default="all.gro") + self.set_option(tgt_opts,'openmm_precision','precision',default="double", forceprint=True) + self.set_option(tgt_opts,'openmm_platform','platname',default="Reference", forceprint=True) + self.engine_ = SMIRNOFF + ## Initialize base class. + super(AbInitioPairwise_SMIRNOFF,self).__init__(options,tgt_opts,forcefield) + + def submit_jobs(self, mvals, AGrad=False, AHess=False): + # we update the self.pgrads here so it's not overwritten in rtarget.py + smirnoff_update_pgrads(self) + class Vibration_SMIRNOFF(Vibration): """ Vibrational frequency matching using using SMIRNOFF format powered by OpenMM. """ def __init__(self,options,tgt_opts,forcefield): diff --git a/src/tests/files/forcefield/ethanol-smirnoff.offxml b/src/tests/files/forcefield/ethanol-smirnoff.offxml new file mode 100644 index 000000000..49132259a --- /dev/null +++ b/src/tests/files/forcefield/ethanol-smirnoff.offxml @@ -0,0 +1,40 @@ + + + + The Open Force Field Initiative + 2026-01-02 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/src/tests/files/targets/ethanol-pairwise/all.gro b/src/tests/files/targets/ethanol-pairwise/all.gro new file mode 100644 index 000000000..57075d55b --- /dev/null +++ b/src/tests/files/targets/ethanol-pairwise/all.gro @@ -0,0 +1,60 @@ +Ethanol conformation 1 + 9 + 1ETH C1 1 0.000 0.000 0.000 + 1ETH C2 2 0.152 0.000 0.000 + 1ETH O 3 0.206 0.124 0.000 + 1ETH HO 4 0.258 0.124 -0.089 + 1ETH H1 5 -0.039 0.103 0.000 + 1ETH H2 6 -0.039 -0.051 0.089 + 1ETH H3 7 -0.039 -0.051 -0.089 + 1ETH H4 8 0.191 -0.051 0.089 + 1ETH H5 9 0.191 -0.051 -0.089 + 2.00000 2.00000 2.00000 +Ethanol conformation 2 + 9 + 1ETH C1 1 0.000 0.000 0.000 + 1ETH C2 2 0.160 0.000 0.000 + 1ETH O 3 0.214 0.124 0.000 + 1ETH HO 4 0.266 0.124 -0.089 + 1ETH H1 5 -0.039 0.103 0.000 + 1ETH H2 6 -0.039 -0.051 0.089 + 1ETH H3 7 -0.039 -0.051 -0.089 + 1ETH H4 8 0.199 -0.051 0.089 + 1ETH H5 9 0.199 -0.051 -0.089 + 2.00000 2.00000 2.00000 +Ethanol conformation 3 + 9 + 1ETH C1 1 0.000 0.000 0.000 + 1ETH C2 2 0.145 0.000 0.000 + 1ETH O 3 0.199 0.124 0.000 + 1ETH HO 4 0.251 0.124 -0.089 + 1ETH H1 5 -0.039 0.103 0.000 + 1ETH H2 6 -0.039 -0.051 0.089 + 1ETH H3 7 -0.039 -0.051 -0.089 + 1ETH H4 8 0.184 -0.051 0.089 + 1ETH H5 9 0.184 -0.051 -0.089 + 2.00000 2.00000 2.00000 +Ethanol conformation 4 + 9 + 1ETH C1 1 0.000 0.000 0.000 + 1ETH C2 2 0.156 0.000 0.000 + 1ETH O 3 0.210 0.124 0.000 + 1ETH HO 4 0.262 0.124 -0.089 + 1ETH H1 5 -0.039 0.103 0.000 + 1ETH H2 6 -0.039 -0.051 0.089 + 1ETH H3 7 -0.039 -0.051 -0.089 + 1ETH H4 8 0.195 -0.051 0.089 + 1ETH H5 9 0.195 -0.051 -0.089 + 2.00000 2.00000 2.00000 +Ethanol conformation 5 + 9 + 1ETH C1 1 0.000 0.000 0.000 + 1ETH C2 2 0.148 0.000 0.000 + 1ETH O 3 0.202 0.124 0.000 + 1ETH HO 4 0.254 0.124 -0.089 + 1ETH H1 5 -0.039 0.103 0.000 + 1ETH H2 6 -0.039 -0.051 0.089 + 1ETH H3 7 -0.039 -0.051 -0.089 + 1ETH H4 8 0.187 -0.051 0.089 + 1ETH H5 9 0.187 -0.051 -0.089 + 2.00000 2.00000 2.00000 diff --git a/src/tests/files/targets/ethanol-pairwise/conf.pdb b/src/tests/files/targets/ethanol-pairwise/conf.pdb new file mode 100644 index 000000000..a99f68869 --- /dev/null +++ b/src/tests/files/targets/ethanol-pairwise/conf.pdb @@ -0,0 +1,21 @@ +REMARK Ethanol for pairwise SMIRNOFF test +CRYST1 20.000 20.000 20.000 90.00 90.00 90.00 P 1 1 +ATOM 1 C1 ETH A 1 0.000 0.000 0.000 1.00 0.00 C +ATOM 2 C2 ETH A 1 1.520 0.000 0.000 1.00 0.00 C +ATOM 3 O ETH A 1 2.061 1.244 0.000 1.00 0.00 O +ATOM 4 HO ETH A 1 2.584 1.244 -0.891 1.00 0.00 H +ATOM 5 H1 ETH A 1 -0.388 1.028 0.000 1.00 0.00 H +ATOM 6 H2 ETH A 1 -0.389 -0.514 0.891 1.00 0.00 H +ATOM 7 H3 ETH A 1 -0.389 -0.514 -0.891 1.00 0.00 H +ATOM 8 H4 ETH A 1 1.909 -0.508 0.893 1.00 0.00 H +ATOM 9 H5 ETH A 1 1.909 -0.508 -0.893 1.00 0.00 H +CONECT 1 2 5 6 7 +CONECT 2 1 3 8 9 +CONECT 3 2 4 +CONECT 4 3 +CONECT 5 1 +CONECT 6 1 +CONECT 7 1 +CONECT 8 2 +CONECT 9 2 +END diff --git a/src/tests/files/targets/ethanol-pairwise/ethanol.mol2 b/src/tests/files/targets/ethanol-pairwise/ethanol.mol2 new file mode 100644 index 000000000..9e7894280 --- /dev/null +++ b/src/tests/files/targets/ethanol-pairwise/ethanol.mol2 @@ -0,0 +1,25 @@ +@MOLECULE +ethanol +9 8 0 0 0 +SMALL +GASTEIGER + +@ATOM + 1 C1 0.0000 0.0000 0.0000 C.3 1 ETH1 0.0000 + 2 C2 1.5200 0.0000 0.0000 C.3 1 ETH1 0.0000 + 3 O 2.0610 1.2440 0.0000 O.3 1 ETH1 0.0000 + 4 HO 2.5840 1.2440 -0.8910 H 1 ETH1 0.0000 + 5 H1 -0.3880 1.0280 0.0000 H 1 ETH1 0.0000 + 6 H2 -0.3890 -0.5140 0.8910 H 1 ETH1 0.0000 + 7 H3 -0.3890 -0.5140 -0.8910 H 1 ETH1 0.0000 + 8 H4 1.9090 -0.5080 0.8930 H 1 ETH1 0.0000 + 9 H5 1.9090 -0.5080 -0.8930 H 1 ETH1 0.0000 +@BOND + 1 1 2 1 + 2 1 5 1 + 3 1 6 1 + 4 1 7 1 + 5 2 3 1 + 6 2 8 1 + 7 2 9 1 + 8 3 4 1 diff --git a/src/tests/files/targets/ethanol-pairwise/ethanol.sdf b/src/tests/files/targets/ethanol-pairwise/ethanol.sdf new file mode 100644 index 000000000..21c86740f --- /dev/null +++ b/src/tests/files/targets/ethanol-pairwise/ethanol.sdf @@ -0,0 +1,23 @@ + + RDKit 3D + + 9 8 0 0 0 0 0 0 0 0999 V2000 + 0.8817 -0.0448 -0.0147 C 0 0 0 0 0 0 0 0 0 0 0 0 + -0.5817 -0.3757 0.0510 C 0 0 0 0 0 0 0 0 0 0 0 0 + -1.3500 0.7581 0.1762 O 0 0 0 0 0 0 0 0 0 0 0 0 + 1.2650 0.1742 1.0122 H 0 0 0 0 0 0 0 0 0 0 0 0 + 1.0165 0.8705 -0.6090 H 0 0 0 0 0 0 0 0 0 0 0 0 + 1.4764 -0.8945 -0.3919 H 0 0 0 0 0 0 0 0 0 0 0 0 + -0.7854 -0.9968 0.9683 H 0 0 0 0 0 0 0 0 0 0 0 0 + -0.8355 -1.0035 -0.8159 H 0 0 0 0 0 0 0 0 0 0 0 0 + -1.0869 1.5126 -0.3762 H 0 0 0 0 0 0 0 0 0 0 0 0 + 1 2 1 0 + 2 3 1 0 + 1 4 1 0 + 1 5 1 0 + 1 6 1 0 + 2 7 1 0 + 2 8 1 0 + 3 9 1 0 +M END +$$$$ diff --git a/src/tests/files/targets/ethanol-pairwise/qdata.txt b/src/tests/files/targets/ethanol-pairwise/qdata.txt new file mode 100644 index 000000000..a17ad484c --- /dev/null +++ b/src/tests/files/targets/ethanol-pairwise/qdata.txt @@ -0,0 +1,20 @@ +JOB 0 +COORDS 0.0000000000e+00 0.0000000000e+00 0.0000000000e+00 2.8723848542e+00 0.0000000000e+00 0.0000000000e+00 3.8947270951e+00 2.3508202360e+00 0.0000000000e+00 4.8830542522e+00 2.3508202360e+00 -1.6837466481e+00 -7.3321402858e-01 1.9426392304e+00 0.0000000000e+00 -7.3510375545e-01 -9.7131961518e-01 1.6837466481e+00 -7.3510375545e-01 -9.7131961518e-01 -1.6837466481e+00 3.6074886097e+00 -9.5998125391e-01 1.6875261019e+00 3.6074886097e+00 -9.5998125391e-01 -1.6875261019e+00 +ENERGY -1.549012000000e+02 + +JOB 1 +COORDS 0.0000000000e+00 0.0000000000e+00 0.0000000000e+00 2.8723848542e+00 0.0000000000e+00 0.0000000000e+00 3.8947270951e+00 2.3508202360e+00 0.0000000000e+00 3.1302955502e+00 3.1130730155e+00 -1.6267785521e+00 -7.3321402858e-01 1.9426392304e+00 0.0000000000e+00 -7.3510375545e-01 -9.7131961518e-01 1.6837466481e+00 -7.3510375545e-01 -9.7131961518e-01 -1.6837466481e+00 3.6074886097e+00 -9.5998125391e-01 1.6875261019e+00 3.6074886097e+00 -9.5998125391e-01 -1.6875261019e+00 +ENERGY -1.548978000000e+02 + +JOB 2 +COORDS 0.0000000000e+00 0.0000000000e+00 0.0000000000e+00 2.8723848542e+00 0.0000000000e+00 0.0000000000e+00 3.8947270951e+00 2.3508202360e+00 0.0000000000e+00 2.2991588662e+00 3.4745239368e+00 5.6968095959e-02 -7.3321402858e-01 1.9426392304e+00 0.0000000000e+00 -7.3510375545e-01 -9.7131961518e-01 1.6837466481e+00 -7.3510375545e-01 -9.7131961518e-01 -1.6837466481e+00 3.6074886097e+00 -9.5998125391e-01 1.6875261019e+00 3.6074886097e+00 -9.5998125391e-01 -1.6875261019e+00 +ENERGY -1.548934000000e+02 + +JOB 3 +COORDS 0.0000000000e+00 0.0000000000e+00 0.0000000000e+00 2.8723848542e+00 0.0000000000e+00 0.0000000000e+00 3.8947270951e+00 2.3508202360e+00 0.0000000000e+00 3.2207808840e+00 3.0737220785e+00 1.6837466481e+00 -7.3321402858e-01 1.9426392304e+00 0.0000000000e+00 -7.3510375545e-01 -9.7131961518e-01 1.6837466481e+00 -7.3510375545e-01 -9.7131961518e-01 -1.6837466481e+00 3.6074886097e+00 -9.5998125391e-01 1.6875261019e+00 3.6074886097e+00 -9.5998125391e-01 -1.6875261019e+00 +ENERGY -1.549002000000e+02 + +JOB 4 +COORDS 0.0000000000e+00 0.0000000000e+00 0.0000000000e+00 2.8723848542e+00 0.0000000000e+00 0.0000000000e+00 3.8947270951e+00 2.3508202360e+00 0.0000000000e+00 4.9735395860e+00 2.3114692990e+00 1.6267785521e+00 -7.3321402858e-01 1.9426392304e+00 0.0000000000e+00 -7.3510375545e-01 -9.7131961518e-01 1.6837466481e+00 -7.3510375545e-01 -9.7131961518e-01 -1.6837466481e+00 3.6074886097e+00 -9.5998125391e-01 1.6875261019e+00 3.6074886097e+00 -9.5998125391e-01 -1.6875261019e+00 +ENERGY -1.548956000000e+02 + diff --git a/src/tests/test_abinitio_pairwise.py b/src/tests/test_abinitio_pairwise.py new file mode 100644 index 000000000..a711d8d3e --- /dev/null +++ b/src/tests/test_abinitio_pairwise.py @@ -0,0 +1,100 @@ +"""Tests for AbInitioPairwise and AbInitioPairwise_SMIRNOFF targets.""" +from __future__ import absolute_import + +import os +import sys +import shutil + +import numpy as np +import pytest + +import forcebalance +import forcebalance.smirnoffio +from .__init__ import ForceBalanceTestCase +from .test_target import TargetTests +from .test_system import skip_openff_py39 + +has_openff_toolkit = True +try: + import openff.toolkit +except ModuleNotFoundError: + has_openff_toolkit = False + +try: + try: + from openmm.app import * + from openmm import * + from openmm.unit import * + except ImportError: + from simtk.openmm.app import * + from simtk.openmm import * + from simtk.unit import * + no_openmm = False +except ImportError: + no_openmm = True + + +@skip_openff_py39 +@pytest.mark.skipif( + not has_openff_toolkit, reason="openff.toolkit not found" +) +@pytest.mark.skipif(no_openmm, reason="OpenMM not found") +class TestAbInitioPairwise_SMIRNOFF(TargetTests): + """Test AbInitioPairwise_SMIRNOFF using a 5-conformation ethanol dataset.""" + + def setup_method(self, method): + super(TestAbInitioPairwise_SMIRNOFF, self).setup_method(method) + test_files_root = os.path.join(os.path.dirname(__file__), 'files') + self.options['root'] = test_files_root + os.chdir(test_files_root) + self.options.update({ + 'jobtype': 'NEWTON', + 'forcefield': ['ethanol-smirnoff.offxml'], + }) + self.tgt_opt.update({ + 'type': 'ABINITIOPAIRWISE_SMIRNOFF', + 'name': 'ethanol-pairwise', + 'mol2': ['ethanol.sdf'], + 'energy': True, + 'force': False, + 'w_energy': 1.0, + 'w_force': 0.0, + }) + self.ff = forcebalance.forcefield.FF(self.options) + self.mvals = np.array([0.0] * self.ff.np) + + self.target = forcebalance.smirnoffio.AbInitioPairwise_SMIRNOFF( + self.options, self.tgt_opt, self.ff + ) + + def teardown_method(self): + shutil.rmtree('temp', ignore_errors=True) + super(TestAbInitioPairwise_SMIRNOFF, self).teardown_method() + + def test_force_raises(self): + """Enabling force fitting should raise RuntimeError.""" + bad_opt = self.tgt_opt.copy() + bad_opt['force'] = True + with pytest.raises(RuntimeError): + forcebalance.smirnoffio.AbInitioPairwise_SMIRNOFF( + self.options, bad_opt, self.ff + ) + + def test_pairwise_pairs_count(self): + """Number of pairs should be n*(n-1)/2 for n snapshots.""" + ns = self.target.ns + expected_pairs = ns * (ns - 1) // 2 + assert len(self.target.eqm_pairs) == expected_pairs + assert len(self.target.boltz_wt_pairs) == expected_pairs + + def test_pairwise_weights_sum_to_one(self): + """Boltzmann weights for pairs should be normalized.""" + assert abs(self.target.boltz_wt_pairs.sum() - 1.0) < 1e-10 + + def test_pairwise_energy_differences(self): + """Pairwise QM energy differences should match manual computation.""" + import itertools + eqm = self.target.eqm + for i, (a, b) in enumerate(itertools.combinations(range(self.target.ns), 2)): + expected = eqm[a] - eqm[b] + assert abs(self.target.eqm_pairs[i] - expected) < 1e-10 diff --git a/src/tests/test_objective.py b/src/tests/test_objective.py index 9c9de12e5..8ba8bdec1 100644 --- a/src/tests/test_objective.py +++ b/src/tests/test_objective.py @@ -40,6 +40,7 @@ def test_no_unlisted_classes_derived_from_Target(self): # Basically, platform-independent targets are excluded. exclude = ['Target', 'AbInitio', + 'AbInitioPairwise', 'Interaction', 'Interaction_GMX', 'Liquid',