Merge pull request #283 from leeping/v195_final_fixes

V195 final fixes

bin/ForceBalance.py

@@ -60,7 +60,7 @@ def process(word, color):
     return Answer
 
 def main():
-    printcool("Welcome to ForceBalance version 1.9.4! =D\nForce Field Optimization System\n\nAuthors:\nLee-Ping Wang\nYudong Qiu, Keri A. McKiernan\nJeffrey R. Wagner, Hyesu Jang, Simon Boothroyd\nArthur Vigil, Erik G. Brandt, John Stoppelman\nJohnny Israeli, Matt Thompson", ansi="1", bold=True, minwidth=64)
+    printcool("Welcome to ForceBalance version 1.9.5! =D\nForce Field Optimization System\n\nAuthors:\nLee-Ping Wang\nYudong Qiu, Keri A. McKiernan\nJeffrey R. Wagner, Hyesu Jang, Simon Boothroyd\nArthur Vigil, Erik G. Brandt, John Stoppelman\nJohnny Israeli, Matt Thompson", ansi="1", bold=True, minwidth=64)
     logostr = """
                           ,'+++                                        
                        ,++++++.      .:,,.                              

devtools/conda-envs/test_env.yaml

@@ -19,12 +19,12 @@ dependencies:
   - swig
   - future
   - pymbar =3
-  - openmm
+  - openmm >= 8
   - ambertools
   - ndcctools
   - geometric
-  - gromacs =2019.1
+  # - gromacs =2019.1
   - openff-toolkit >=0.11.3
-  - openff-evaluator-base
-  - openff-recharge
+  - openff-evaluator-base >= 0.4.1
+  # - openff-recharge
   - openeye-toolkits

devtools/conda-recipe/meta.yaml

@@ -1,6 +1,6 @@
 package:
   name: forcebalance-dev
-  version: !!str 1.9.4
+  version: !!str 1.9.5
 
 source:
   path: ../..

doc/api_header.tex

@@ -57,7 +57,7 @@
 \vspace*{1cm}
 \begin{center}
 
-{\Large ForceBalance Developer API Guide version 1.9.4}\\
+{\Large ForceBalance Developer API Guide version 1.9.5}\\
 \vspace*{2cm}
 {\large Generated by Doxygen 1.8.11}\\
 \vspace*{2.5 cm}

doc/doxygen.cfg

@@ -32,7 +32,7 @@ PROJECT_NAME           = ForceBalance
 # This could be handy for archiving the generated documentation or
 # if some version control system is used.
 
-PROJECT_NUMBER         = 1.9.4
+PROJECT_NUMBER         = 1.9.5
 
 # Using the PROJECT_BRIEF tag one can provide an optional one line description
 # for a project that appears at the top of each page and should give viewer

doc/header.tex

@@ -57,7 +57,7 @@
 \vspace*{1cm}
 \begin{center}
 
-{\Large ForceBalance version 1.9.4}\\
+{\Large ForceBalance version 1.9.5}\\
 \vspace*{2cm}
 {\large Generated by Doxygen 1.8.11}\\
 \vspace*{2.5 cm}

setup.py

@@ -21,10 +21,10 @@
 #| doc/api_header.tex              |#
 #| bin/ForceBalance.py             |#
 #===================================#
-__version__ = "v1.9.4"
+__version__ = "v1.9.5"
 try:
     # use git to find current version
-    git_describe = subprocess.check_output(["git", "describe"]).strip()
+    git_describe = subprocess.check_output(["git", "describe", "--tags"]).strip()
     __version__ = re.sub('-g[0-9a-f]*$','',git_describe)
 except: pass
 

src/evaluator_io.py

@@ -308,6 +308,12 @@ def _parameter_value_from_gradient_key(self, gradient_key):
         bool
             Returns True if the parameter is a cosmetic one.
         """
+        # try:
+        #     import openmm.unit as simtk_unit
+        # except ImportError:
+        #     import simtk.unit as simtk_unit
+        from openff.units import unit as openff_unit
+
 
         parameter_handler = self.FF.openff_forcefield.get_parameter_handler(
             gradient_key.tag
@@ -346,6 +352,10 @@ def _parameter_value_from_gradient_key(self, gradient_key):
         ):
             is_cosmetic = True
 
+        if not isinstance(parameter_value, openff_unit.Quantity):
+            parameter_value = parameter_value * openff_unit.dimensionless
+
+        #return openmm_quantity_to_pint(parameter_value), is_cosmetic
         return parameter_value, is_cosmetic
 
     def _extract_physical_parameter_values(self):

src/gmxio.py

@@ -544,7 +544,7 @@ def setopts(self, **kwargs):
         if 'gmxsuffix' in kwargs:
             self.gmxsuffix = kwargs['gmxsuffix']
         else:
-            warn_once("The 'gmxsuffix' option were not provided; using default.")
+            warn_once("The 'gmxsuffix' option were not provided; using default (blank).")
             self.gmxsuffix = ''
 
         ## Barostat keyword for equilibration
@@ -893,11 +893,12 @@ def optimize(self, shot, crit=1e-4, align=True, **kwargs):
 
         self.warngmx("grompp -c %s.gro -p %s.top -f %s-min.mdp -o %s-min.tpr" % (self.name, self.name, self.name, self.name))
         self.callgmx("mdrun -deffnm %s-min -nt 1" % self.name)
-        self.callgmx("trjconv -f %s-min.trr -s %s-min.tpr -o %s-min.gro -ndec 9" % (self.name, self.name, self.name), stdin="System")
+        # self.callgmx("trjconv -f %s-min.trr -s %s-min.tpr -o %s-min.gro -ndec 9" % (self.name, self.name, self.name), stdin="System")
+        self.callgmx("trjconv -f %s-min.trr -s %s-min.tpr -o %s-min.g96" % (self.name, self.name, self.name), stdin="System")
         self.callgmx("g_energy -xvg no -f %s-min.edr -o %s-min-e.xvg" % (self.name, self.name), stdin='Potential')
 
         E = float(open("%s-min-e.xvg" % self.name).readlines()[-1].split()[1])
-        M = Molecule("%s.gro" % self.name, build_topology=False) + Molecule("%s-min.gro" % self.name, build_topology=False)
+        M = Molecule("%s.gro" % self.name, build_topology=False) + Molecule("%s-min.g96" % self.name)
         if not self.pbc:
             M.align(center=False)
         rmsd = M.ref_rmsd(0)[1]
@@ -968,15 +969,6 @@ def evaluate_trajectory(self, force=False, dipole=False, traj=None):
         self.mol[0].write("%s.gro" % self.name)
         return self.evaluate_(force, dipole, traj)
 
-    def make_gro_trajectory(self, fout=None):
-        """ Return the MD trajectory as a Molecule object. """
-        if fout is None:
-            fout = "%s-mdtraj.gro" % self.name
-        if not hasattr(self, 'mdtraj'):
-            raise RuntimeError('Engine does not have an associated trajectory.')
-        self.callgmx("trjconv -f %s -o %s -ndec 9 -pbc mol -novel -noforce" % (self.mdtraj, fout), stdin='System')
-        return fout
-
     def energy_one(self, shot):
 
         """ Compute the energy using GROMACS for a snapshot. """
@@ -1126,8 +1118,8 @@ def normal_modes(self, shot=0, optimize=True):
         self.callgmx("mdrun -deffnm %s-nm -nt 1 -mtx %s-nm.mtx -v" % (self.name, self.name))
         self.callgmx("g_nmeig -s %s-nm.tpr -f %s-nm.mtx -of %s-nm.xvg -v %s-nm.trr -last 10000 -xvg no" % \
                          (self.name, self.name, self.name, self.name))
-        self.callgmx("trjconv -s %s-nm.tpr -f %s-nm.trr -o %s-nm.gro -ndec 9" % (self.name, self.name, self.name), stdin="System")
-        NM = Molecule("%s-nm.gro" % self.name, build_topology=False)
+        self.callgmx("trjconv -s %s-nm.tpr -f %s-nm.trr -o %s-nm.g96" % (self.name, self.name, self.name), stdin="System")
+        NM = Molecule("%s-nm.g96" % self.name, build_topology=False)
 
         calc_eigvals = np.array([float(line.split()[1]) for line in open("%s-nm.xvg" % self.name).readlines()])
         calc_eigvecs = NM.xyzs[1:]
@@ -1148,8 +1140,8 @@ def generate_positions(self):
         ## Call grompp followed by mdrun.
         self.warngmx("grompp -c %s.gro -p %s.top -f %s.mdp -o %s.tpr" % (self.name, self.name, self.name, self.name))
         self.callgmx("mdrun -deffnm %s -nt 1 -rerunvsite -rerun %s-all.gro" % (self.name, self.name))
-        self.callgmx("trjconv -f %s.trr -o %s-out.gro -ndec 9 -novel -noforce" % (self.name, self.name), stdin='System')
-        NewMol = Molecule("%s-out.gro" % self.name)
+        self.callgmx("trjconv -f %s.trr -o %s-out.g96 -novel -noforce" % (self.name, self.name), stdin='System')
+        NewMol = Molecule("%s-out.g96" % self.name, build_topology-False)
         return NewMol.xyzs
 
     def n_snaps(self, nsteps, step_interval, timestep):

src/hessian.py

@@ -201,9 +201,9 @@ def compute(mvals_):
                 outfile.writelines('%s\n' % line for line in compute.M_opt.write_xyz([0]))
                 outfile.write('\n')
                 for freq, normal_mode in zip(freqs_rearr, normal_modes_rearr):
-                    outfile.write(f'{freq}\n')
+                    outfile.write("%13.4f\n" % freq)
                     for nx, ny, nz in normal_mode:
-                        outfile.write(f'{nx:13.4f} {ny:13.4f} {nz:13.4f}\n')
+                        outfile.write("%13.4f %13.4f %13.4f\n" % (nx, ny, nz))
                     outfile.write('\n')
             outfile.close()
 
@@ -237,13 +237,13 @@ def render_normal_modes(elem, xyz, eigvecs, color, qm=False, ref_eigvals=None, e
             x, y, z = xyz.T
             u, v, w = eigvec.T *5
             origin = np.array([x, y, z])
-            axs[idx].quiver(*origin, u, v, w, color=color)
+            axs[idx].quiver(origin[0], origin[1], origin[2], u, v, w, color=color)
 
             axs[idx].set_xlabel('x')
             axs[idx].set_ylabel('y')
             axs[idx].set_zlabel('z')
             if qm:
-                axs[idx].set_title(f'normal mode #{idx} (blue:QM({ref_eigvals[idx]:.2f}), red:MM({eigvals_rearr[idx]:.2f}))')
+                axs[idx].set_title('normal mode #%i (blue:QM(%.2f), red:MM(%.2f))' % (idx, ref_eigvals[idx], eigvals_rearr[idx]))
                 axs[idx].scatter(x, y, z, color='black', s=30)
                 axs[idx].set_xlim(min(u+x), max(u+x))
                 axs[idx].set_ylim(min(v+y), max(v+y))
@@ -273,7 +273,7 @@ def draw_vibfreq_scatter_plot_n_overlap_matrix(name, engine, ref_eigvals, ref_ei
     axs[0].set_xlabel(r'QM vibrational frequency ($cm^{-1}$)')
     axs[0].set_ylabel(r'MM vibrational frequency ($cm^{-1}$)')
     mae = np.sum(np.abs(ref_eigvals - freqs_rearr))/ len(ref_eigvals)
-    axs[0].set_title(f'QM vs. MM vibrational frequencies\n MAE= {mae:.2f}')
+    axs[0].set_title('QM vs. MM vibrational frequencies\n MAE= %2f' % mae)
     x0,x1 = axs[0].get_xlim()
     y0,y1 = axs[0].get_ylim()
     axs[0].set_aspect((x1-x0)/(y1-y0))
@@ -297,7 +297,7 @@ def draw_vibfreq_scatter_plot_n_overlap_matrix(name, engine, ref_eigvals, ref_ei
     plt.colorbar(im, cax=cax)
     corr_coef = cal_corr_coef(overlap_matrix)
     err = np.linalg.norm(qm_overlap_matrix - overlap_matrix)/np.linalg.norm(qm_overlap_matrix) # measure of error in matrix (Relative error)
-    axs[1].set_title(f'QM vs. MM normal modes\n Correlation coef. ={corr_coef:.4f}, Error={err:.4f}')
+    axs[1].set_title('QM vs. MM normal modes\n Correlation coef. =%.4f, Error=%.4f' % (corr_coef, err))
 
     # # move ax x axis to top
     # axs[2].xaxis.tick_top()
@@ -321,4 +321,4 @@ def draw_vibfreq_scatter_plot_n_overlap_matrix(name, engine, ref_eigvals, ref_ei
     plt.tight_layout()
     plt.subplots_adjust(top=0.85)
     fig.suptitle('Hessian: iteration %i\nSystem: %s' % (Counter(), name))
-    fig.savefig('vibfreq_scatter_plot_n_overlap_matrix.pdf')
+    fig.savefig('vibfreq_scatter_plot_n_overlap_matrix.pdf')

src/molecule.py

@@ -28,7 +28,7 @@
 try:
     import importlib
     package_install_dir = os.path.split(importlib.util.find_spec(__name__.split('.')[0]).origin)[0]
-except ImportError: # importlib is new in version 3.1
+except AttributeError: # importlib is new in version 3.1
     import imp
     package_install_dir = imp.find_module(__name__.split('.')[0])[1]
 
@@ -1209,6 +1209,7 @@ def __init__(self, fnm = None, ftype = None, top = None, ttype = None, **kwargs)
         ## The table of file readers
         self.Read_Tab = {'gaussian' : self.read_com,
                          'gromacs'  : self.read_gro,
+                         'g96'      : self.read_g96,
                          'charmm'   : self.read_charmm,
                          'dcd'      : self.read_dcd,
                          'mdcrd'    : self.read_mdcrd,
@@ -3382,6 +3383,61 @@ def read_gro(self, fnm, **kwargs):
                   'comms'    : comms}
         return Answer
 
+    def read_g96(self, fnm, **kwargs):
+        """ Read a GROMACS .g96 file.
+        This was implemented because as of 2018, GROMACS no longer writes .gro files in variable precision
+        therefore we need to write out .g96 files to get precise coordinates.
+        """
+        xyzs     = []
+        comms    = []
+        boxes    = []
+        xyz      = []
+        read_mode = 'None'
+        for line in open(fnm):
+            sline = line.split()
+            if len(sline) == 1 and sline[0] == 'TITLE':
+                read_mode = 'TITLE'
+            elif len(sline) == 1 and sline[0] in ['POSITION', 'POSITIONRED']:
+                read_mode = 'POSITION'
+            elif len(sline) == 1 and sline[0] == 'BOX':
+                read_mode = 'BOX'
+            elif len(sline) == 1 and sline[0] == 'END':
+                # comms and xyzs should only be incremented when we encounter the END of a TITLE or POSITION section
+                if read_mode == 'POSITION':
+                    xyzs.append(np.array(xyz))
+                    xyz = []
+                elif read_mode == 'TITLE':
+                    comms.append(title)
+                read_mode = 'None'
+            elif read_mode == 'BOX' and len(sline) > 0:
+                # Read box information (should be on a single line)
+                box = [float(i)*10 for i in sline]
+                if len(box) == 3:
+                    a = box[0]
+                    b = box[1]
+                    c = box[2]
+                    alpha = 90.0
+                    beta = 90.0
+                    gamma = 90.0
+                    boxes.append(BuildLatticeFromLengthsAngles(a, b, c, alpha, beta, gamma))
+                elif len(box) == 9:
+                    v1 = np.array([box[0], box[3], box[4]])
+                    v2 = np.array([box[5], box[1], box[6]])
+                    v3 = np.array([box[7], box[8], box[2]])
+                    boxes.append(BuildLatticeFromVectors(v1, v2, v3))
+            elif read_mode == 'POSITION' and len(sline) > 0:
+                xyz.append([float(i)*10 for i in sline])
+            elif read_mode == 'TITLE':
+                title = line.strip()
+        if len(xyzs) != len(boxes):
+            raise IOError('in read_g96: xyzs and boxes should have the same length')
+        if len(xyzs) != len(comms):
+            raise IOError('in read_g96: xyzs and comms should have the same length')
+        Answer = {'xyzs'     : xyzs,
+                  'boxes'    : boxes,
+                  'comms'    : comms}
+        return Answer
+
     def read_charmm(self, fnm, **kwargs):
         """ Read a CHARMM .cor (or .crd) file.
 

src/openmmio.py

@@ -199,7 +199,7 @@ def vsfunc(pos, idx_, wt_):
                     zdir /= np.linalg.norm(zdir)
                     return origin + np.array(np.array([xdir, ydir, zdir]) * vpos[:, None]).sum(axis=0)
             else:
-                raise NotImplementedError(f"The virtual site type {vs.__class__.__name__} is not currently supported.")
+                raise NotImplementedError("The virtual site type %s is not currently supported." % vs.__class__.__name__)
 
         else:
             isvsites.append(0)

src/smirnoff_hack.py

@@ -36,7 +36,9 @@ def hash_molecule(molecule):
 def hash_molecule_args_and_kwargs(molecule, *args, **kwargs):
     arguments = [str(arg) for arg in args]
     keywords = [str(keyword) for keyword in kwargs.values()]
-    return hash((hash_molecule(molecule), *arguments, *keywords))
+    arguments_plus_keywords = arguments + keywords
+    # Deleted * in front of arguments_plus_keywords for Python 2.7 compatibility, convert list to tuple.
+    return hash((hash_molecule(molecule), tuple(arguments_plus_keywords)))
 
 
 if _SHOULD_CACHE:

src/smirnoffio.py

@@ -38,10 +38,18 @@
 
 logger = getLogger(__name__)
 try:
-    from simtk.openmm.app import *
-    from simtk.openmm import *
-    from simtk.unit import *
-    import simtk.openmm._openmm as _openmm
+    try:
+        # Try importing openmm using >=7.6 namespace
+        from openmm.app import *
+        from openmm import *
+        from openmm.unit import *
+        import openmm._openmm as _openmm
+    except ImportError:
+        # Try importing openmm using <7.6 namespace
+        from simtk.openmm.app import *
+        from simtk.openmm import *
+        from simtk.unit import *
+        import simtk.openmm._openmm as _openmm
 except:
     pass
 

src/tests/test_system.py

@@ -34,8 +34,10 @@
 # expected results taken from previous runs. Update this if it changes and seems reasonable (updated 07/23/14)
 EXPECTED_LIPID_RESULTS = array([-6.7553e-03, -2.4070e-02])
 
-# expected results taken from OpenFF torsion profile optimization using OpenFF toolkit 0.4.1 and OpenEye toolkit 2019.10.2. (updated 11/02/19)
-EXPECTED_OPENFF_TORSIONPROFILE_RESULTS = array([-9.4238e-02, 7.3350e-03, -7.9467e-05, 1.7172e-02, -1.3309e-01, 6.0076e-02, 1.7895e-02, 6.5866e-02, -1.4084e-01, -2.2906e-02])
+# expected results taken from OpenFF torsion profile optimization using OpenFF toolkit 0.4.1 and OpenEye toolkit 2019.10.2. (updated 02/06/23)
+# EXPECTED_OPENFF_TORSIONPROFILE_RESULTS = array([-9.4238e-02, 7.3350e-03, -7.9467e-05, 1.7172e-02, -1.3309e-01, 6.0076e-02, 1.7895e-02, 6.5866e-02, -1.4084e-01, -2.2906e-02])
+# 02/06/23: As of toolkit v0.11, the default charge assignment method changed, which caused the following change in the optimization result:
+EXPECTED_OPENFF_TORSIONPROFILE_RESULTS = array([-8.6810e-02, 6.7106e-03, 3.0992e-03, 1.8605e-02, -1.1292e-01, 5.6741e-02, 1.8884e-02, 7.3325e-02, -1.4203e-01, -9.2920e-03])
 
 # expected objective function from 025 recharge methane study. (updated 08/04/20)
 EXPECTED_RECHARGE_METHANE_ESP_OBJECTIVE = array([5.68107e-04])