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])