1
#!/usr/bin/env python
2

3
# =============================================================================================
4
# MODULE DOCSTRING
5
# =============================================================================================
6 8
"""
7
Parameter assignment tools for the SMIRNOFF (SMIRKS Native Open Force Field) format.
8

9
.. codeauthor:: John D. Chodera <john.chodera@choderalab.org>
10
.. codeauthor:: David L. Mobley <dmobley@mobleylab.org>
11
.. codeauthor:: Peter K. Eastman <peastman@stanford.edu>
12

13
.. todo ::
14

15
   * Speed up overall import time by putting non-global imports only where they are needed
16

17
"""
18

19 8
__all__ = [
20
    "get_available_force_fields",
21
    "MAX_SUPPORTED_VERSION",
22
    "ParameterHandlerRegistrationError",
23
    "SMIRNOFFVersionError",
24
    "SMIRNOFFAromaticityError",
25
    "ParseError",
26
    "ForceField",
27
]
28

29

30
# =============================================================================================
31
# GLOBAL IMPORTS
32
# =============================================================================================
33

34 8
import copy
35 8
import logging
36 8
import os
37 8
import pathlib
38 8
from collections import OrderedDict
39

40 8
from simtk import openmm, unit
41

42 8
from openforcefield.topology.molecule import DEFAULT_AROMATICITY_MODEL
43 8
from openforcefield.typing.engines.smirnoff.io import ParameterIOHandler
44 8
from openforcefield.typing.engines.smirnoff.parameters import ParameterHandler
45 8
from openforcefield.typing.engines.smirnoff.plugins import load_handler_plugins
46 8
from openforcefield.utils import (
47
    MessageException,
48
    all_subclasses,
49
    convert_0_1_smirnoff_to_0_2,
50
    convert_0_2_smirnoff_to_0_3,
51
    convert_all_quantities_to_string,
52
    convert_all_strings_to_quantity,
53
)
54

55
# =============================================================================================
56
# CONFIGURE LOGGER
57
# =============================================================================================
58

59 8
logger = logging.getLogger(__name__)
60

61
# =============================================================================================
62
# PRIVATE METHODS
63
# =============================================================================================
64

65
# Directory paths used by ForceField to discover offxml files.
66 8
_installed_offxml_dir_paths = []
67

68

69 8
def _get_installed_offxml_dir_paths():
70
    """Return the list of directory paths where to search for offxml files.
71

72
    This function load the information by calling all the entry points
73
    registered in the "openff.forcefielddirs" group. Each entry point
74
    (i.e. a function) should return a list of paths to directories
75
    containing the offxml files.
76

77
    Returns
78
    -------
79
    installed_offxml_dir_paths : List[str]
80
        All the installed directory paths where ``ForceField`` will
81
        look for offxml files.
82

83
    """
84
    global _installed_offxml_dir_paths
85 8
    if len(_installed_offxml_dir_paths) == 0:
86 8
        from pkg_resources import iter_entry_points
87

88
        # Find all registered entry points that should return a list of
89
        # paths to directories where to search for offxml files.
90 8
        for entry_point in iter_entry_points(
91
            group="openforcefield.smirnoff_forcefield_directory"
92
        ):
93 8
            _installed_offxml_dir_paths.extend(entry_point.load()())
94 8
    return _installed_offxml_dir_paths
95

96

97 8
def get_available_force_fields(full_paths=False):
98
    """
99
     Get the filenames of all available .offxml force field files.
100

101
     Availability is determined by what is discovered through the
102
    `openforcefield.smirnoff_forcefield_directory` entry point. If the
103
    `openforcefields` package is installed, this should include several
104
    .offxml files such as `openff-1.0.0.offxml`.
105

106
     Parameters
107
     ----------
108
     full_paths : bool, default=False
109
         If False, return the name of each available *.offxml file.
110
         If True, return the full path to each available .offxml file.
111

112
     Returns
113
     -------
114
     available_force_fields : List[str]
115
         List of available force field files
116

117
    """
118 8
    installed_paths = _get_installed_offxml_dir_paths()
119 8
    available_force_fields = []
120 8
    for installed_path in installed_paths:
121 8
        for globbed in pathlib.Path(installed_path).rglob("*.offxml"):
122 8
            if full_paths:
123 0
                available_force_fields.append(globbed.as_posix())
124
            else:
125 8
                available_force_fields.append(globbed.name)
126 8
    return available_force_fields
127

128

129
# TODO: Instead of having a global version number, alow each Force to have a separate version number
130 8
MAX_SUPPORTED_VERSION = (
131
    "1.0"  # maximum version of the SMIRNOFF spec supported by this SMIRNOFF forcefield
132
)
133

134

135 8
class ParameterHandlerRegistrationError(MessageException):
136
    """
137
    Exception for errors in ParameterHandler registration
138
    """
139

140 8
    pass
141

142

143 8
class SMIRNOFFVersionError(MessageException):
144
    """
145
    Exception thrown when an incompatible SMIRNOFF version data structure in attempted to be read.
146
    """
147

148 8
    pass
149

150

151 8
class SMIRNOFFAromaticityError(MessageException):
152
    """
153
    Exception thrown when an incompatible SMIRNOFF aromaticity model is checked for compatibility.
154
    """
155

156 8
    pass
157

158

159 8
class ParseError(MessageException):
160
    """
161
    Error for when a SMIRNOFF data structure is not parseable by a ForceField
162
    """
163

164 8
    pass
165

166

167
# =============================================================================================
168
# FORCEFIELD
169
# =============================================================================================
170

171
# QUESTION: How should we document private object fields?
172

173
# TODO: How do we serialize/deserialize `ForceField`'s object model? Can we rely on pickle?
174

175

176 8
class ForceField:
177
    """A factory that assigns SMIRNOFF parameters to a molecular system
178

179
    :class:`ForceField` is a factory that constructs an OpenMM :class:`simtk.openmm.System` object from a
180
    :class:`openforcefield.topology.Topology` object defining a (bio)molecular system containing one or more molecules.
181

182
    When a :class:`ForceField` object is created from one or more specified SMIRNOFF serialized representations,
183
    all :class:`ParameterHandler` subclasses currently imported are identified and registered to handle different
184
    sections of the SMIRNOFF force field definition file(s).
185

186
    All :class:`ParameterIOHandler` subclasses currently imported are identified and registered to handle different
187
    serialization formats (such as XML).
188

189
    The force field definition is processed by these handlers to populate the ``ForceField`` object model data
190
    structures that can easily be manipulated via the API:
191

192
    Processing a :class:`Topology` object defining a chemical system will then call all :class`ParameterHandler`
193
    objects in an order guaranteed to satisfy the declared processing order constraints of each
194
    :class`ParameterHandler`.
195

196
    Attributes
197
    ----------
198
    parameters : dict of str : list of ParameterType
199
        ``parameters[tagname]`` is the instantiated :class:`ParameterHandler` class that handles parameters associated
200
        with the force ``tagname``.
201
        This is the primary means of retrieving and modifying parameters, such as
202
        ``parameters['vdW'][0].sigma *= 1.1``
203
    parameter_object_handlers : dict of str : ParameterHandler class
204
        Registered list of :class:`ParameterHandler` classes that will handle different forcefield tags to create the parameter object model.
205
        ``parameter_object_handlers[tagname]`` is the :class:`ParameterHandler` that will be instantiated to process the force field definition section ``tagname``.
206
        :class:`ParameterHandler` classes are registered when the ForceField object is created, but can be manipulated afterwards.
207
    parameter_io_handlers : dict of str : ParameterIOHandler class
208
        Registered list of :class:`ParameterIOHandler` classes that will handle serializing/deserializing the parameter object model to string or file representations, such as XML.
209
        ``parameter_io_handlers[iotype]`` is the :class:`ParameterHandler` that will be instantiated to process the serialization scheme ``iotype``.
210
        :class:`ParameterIOHandler` classes are registered when the ForceField object is created, but can be manipulated afterwards.
211

212
    Examples
213
    --------
214

215
    Create a new ForceField containing the smirnoff99Frosst parameter set:
216

217
    >>> from openforcefield.typing.engines.smirnoff import ForceField
218
    >>> forcefield = ForceField('test_forcefields/smirnoff99Frosst.offxml')
219

220
    Create an OpenMM system from a :class:`openforcefield.topology.Topology` object:
221

222
    >>> from openforcefield.topology import Molecule, Topology
223
    >>> ethanol = Molecule.from_smiles('CCO')
224
    >>> topology = Topology.from_molecules(molecules=[ethanol])
225
    >>> system = forcefield.create_openmm_system(topology)
226

227
    Modify the long-range electrostatics method:
228

229
    >>> forcefield.get_parameter_handler('Electrostatics').method = 'PME'
230

231
    Inspect the first few vdW parameters:
232

233
    >>> low_precedence_parameters = forcefield.get_parameter_handler('vdW').parameters[0:3]
234

235
    Retrieve the vdW parameters by SMIRKS string and manipulate it:
236

237
    >>> parameter = forcefield.get_parameter_handler('vdW').parameters['[#1:1]-[#7]']
238
    >>> parameter.rmin_half += 0.1 * unit.angstroms
239
    >>> parameter.epsilon *= 1.02
240

241
    Make a child vdW type more specific (checking modified SMIRKS for validity):
242

243
    >>> forcefield.get_parameter_handler('vdW').parameters[-1].smirks += '~[#53]'
244

245
    .. warning ::
246

247
       While we check whether the modified SMIRKS is still valid and has the appropriate valence type,
248
       we currently don't check whether the typing remains hierarchical, which could result in some types
249
       no longer being assignable because more general types now come *below* them and preferentially match.
250

251
    Delete a parameter:
252

253
    >>> del forcefield.get_parameter_handler('vdW').parameters['[#1:1]-[#6X4]']
254

255
    Insert a parameter at a specific point in the parameter tree:
256

257
    >>> from openforcefield.typing.engines.smirnoff import vdWHandler
258
    >>> new_parameter = vdWHandler.vdWType(smirks='[*:1]', epsilon=0.0157*unit.kilocalories_per_mole, rmin_half=0.6000*unit.angstroms)
259
    >>> forcefield.get_parameter_handler('vdW').parameters.insert(0, new_parameter)
260

261
    .. warning ::
262

263
       We currently don't check whether removing a parameter could accidentally remove the root type, so it's possible to no longer type all molecules this way.
264

265
    """
266

267 8
    def __init__(
268
        self,
269
        *sources,
270
        aromaticity_model=DEFAULT_AROMATICITY_MODEL,
271
        parameter_handler_classes=None,
272
        parameter_io_handler_classes=None,
273
        disable_version_check=False,
274
        allow_cosmetic_attributes=False,
275
        load_plugins=False,
276
    ):
277
        """Create a new :class:`ForceField` object from one or more SMIRNOFF parameter definition files.
278

279
        Parameters
280
        ----------
281
        sources : string or file-like object or open file handle or URL (or iterable of these)
282
            A list of files defining the SMIRNOFF force field to be loaded.
283
            Currently, only `the SMIRNOFF XML format <https://github.com/openforcefield/openforcefield/blob/master/The-SMIRNOFF-force-field-format.md>`_ is supported.
284
            Each entry may be an absolute file path, a path relative to the current working directory, a path relative to this module's data subdirectory
285
            (for built in force fields), or an open file-like object with a ``read()`` method from which the forcefield XML data can be loaded.
286
            If multiple files are specified, any top-level tags that are repeated will be merged if they are compatible,
287
            with files appearing later in the sequence resulting in parameters that have higher precedence.
288
            Support for multiple files is primarily intended to allow solvent parameters to be specified by listing them last in the sequence.
289
        aromaticity_model : string, default='OEAroModel_MDL'
290
            The aromaticity model used by the force field. Currently, only 'OEAroModel_MDL' is supported
291
        parameter_handler_classes : iterable of ParameterHandler classes, optional, default=None
292
            If not None, the specified set of ParameterHandler classes will be instantiated to create the parameter object model.
293
            By default, all imported subclasses of ParameterHandler are automatically registered.
294
        parameter_io_handler_classes : iterable of ParameterIOHandler classes
295
            If not None, the specified set of ParameterIOHandler classes will be used to parse/generate serialized parameter sets.
296
            By default, all imported subclasses of ParameterIOHandler are automatically registered.
297
        disable_version_check : bool, optional, default=False
298
            If True, will disable checks against the current highest supported forcefield version.
299
            This option is primarily intended for forcefield development.
300
        allow_cosmetic_attributes : bool, optional. Default = False
301
            Whether to retain non-spec kwargs from data sources.
302
        load_plugins: bool, optional. Default = False
303
            Whether to load ``ParameterHandler`` classes which have been registered
304
            by installed plugins.
305

306
        Examples
307
        --------
308

309
        Load one SMIRNOFF parameter set in XML format (searching the package data directory by default, which includes some standard parameter sets):
310

311
        >>> forcefield = ForceField('test_forcefields/smirnoff99Frosst.offxml')
312

313
        Load multiple SMIRNOFF parameter sets:
314

315
        forcefield = ForceField('test_forcefields/smirnoff99Frosst.offxml', 'test_forcefields/tip3p.offxml')
316

317
        Load a parameter set from a string:
318

319
        >>> offxml = '<SMIRNOFF version="0.2" aromaticity_model="OEAroModel_MDL"/>'
320
        >>> forcefield = ForceField(offxml)
321

322
        """
323
        # Clear all object fields
324 8
        self._initialize()
325

326 8
        self.aromaticity_model = aromaticity_model
327
        # Store initialization options
328 8
        self.disable_version_check = disable_version_check
329
        # if True, we won't check which SMIRNOFF version number we're parsing
330

331
        # Register all ParameterHandler objects that will process SMIRNOFF force definitions
332
        # TODO: We need to change this to just find all ParameterHandler objects in this file;
333
        # otherwise, we can't define two different ParameterHandler subclasses to compare for a new type of energy term
334
        # since both will try to register themselves for the same XML tag and an Exception will be raised.
335 8
        if parameter_handler_classes is None:
336 8
            parameter_handler_classes = all_subclasses(ParameterHandler)
337 8
        if load_plugins:
338

339 8
            registered_handlers = load_handler_plugins()
340

341
            # Make sure the same handlers aren't added twice.
342 8
            parameter_handler_classes += [
343
                handler
344
                for handler in registered_handlers
345
                if handler not in parameter_handler_classes
346
            ]
347

348 8
        self._register_parameter_handler_classes(parameter_handler_classes)
349

350
        # Register all ParameterIOHandler objects that will process serialized parameter representations
351 8
        if parameter_io_handler_classes is None:
352 8
            parameter_io_handler_classes = all_subclasses(ParameterIOHandler)
353

354 8
        self._register_parameter_io_handler_classes(parameter_io_handler_classes)
355

356
        # Parse all sources containing SMIRNOFF parameter definitions
357 8
        self.parse_sources(sources, allow_cosmetic_attributes=allow_cosmetic_attributes)
358

359 8
    def _initialize(self):
360
        """
361
        Initialize all object fields.
362
        """
363 8
        self._MIN_SUPPORTED_SMIRNOFF_VERSION = 0.1
364 8
        self._MAX_SUPPORTED_SMIRNOFF_VERSION = 0.3
365 8
        self._disable_version_check = (
366
            False  # if True, will disable checking compatibility version
367
        )
368 8
        self._aromaticity_model = None
369 8
        self._parameter_handler_classes = (
370
            OrderedDict()
371
        )  # Parameter handler classes that _can_ be initialized if needed
372 8
        self._parameter_handlers = (
373
            OrderedDict()
374
        )  # ParameterHandler classes to be instantiated for each parameter type
375 8
        self._parameter_io_handler_classes = (
376
            OrderedDict()
377
        )  # ParameterIOHandler classes that _can_ be initialiazed if needed
378 8
        self._parameter_io_handlers = (
379
            OrderedDict()
380
        )  # ParameterIO classes to be used for each file type
381 8
        self._author = None
382 8
        self._date = None
383

384 8
    def _check_smirnoff_version_compatibility(self, version):
385
        """
386
        Raise a parsing exception if the given file version is incompatible with this ForceField class.
387

388
        Parameters
389
        ----------
390
        version : str
391
            The SMIRNOFF version being read.
392

393
        Raises
394
        ------
395
        SMIRNOFFVersionError if an incompatible version is passed in.
396

397
        """
398 8
        import packaging.version
399

400
        # Use PEP-440 compliant version number comparison, if requested
401 8
        if (not self.disable_version_check) and (
402
            (
403
                packaging.version.parse(str(version))
404
                > packaging.version.parse(str(self._MAX_SUPPORTED_SMIRNOFF_VERSION))
405
            )
406
            or (
407
                packaging.version.parse(str(version))
408
                < packaging.version.parse(str(self._MIN_SUPPORTED_SMIRNOFF_VERSION))
409
            )
410
        ):
411 0
            raise SMIRNOFFVersionError(
412
                "SMIRNOFF offxml file was written with version {}, but this version of ForceField only supports "
413
                "version {} to version {}".format(
414
                    version,
415
                    self._MIN_SUPPORTED_SMIRNOFF_VERSION,
416
                    self._MAX_SUPPORTED_SMIRNOFF_VERSION,
417
                )
418
            )
419

420 8
    @property
421
    def aromaticity_model(self):
422
        """Returns the aromaticity model for this ForceField object.
423

424
        Returns
425
        -------
426
        aromaticity_model
427
            The aromaticity model for this force field.
428
        """
429 8
        return self._aromaticity_model
430

431 8
    @aromaticity_model.setter
432
    def aromaticity_model(self, aromaticity_model):
433
        """
434
        Register that this forcefield is using an aromaticity model. Will check for
435
        compatibility with other aromaticity model(s) already in use.
436

437
        Parameters
438
        ----------
439
        aromaticity_model : str
440
            The aromaticity model to register.
441

442
        Raises
443
        ------
444
        SMIRNOFFAromaticityError if an incompatible aromaticity model is passed in.
445

446
        .. notes ::
447
           * Currently, the only supported aromaticity model is 'OEAroModel_MDL'.
448

449
        """
450
        # Implement better logic here if we ever support another aromaticity model
451 8
        if aromaticity_model != "OEAroModel_MDL":
452 8
            raise SMIRNOFFAromaticityError(
453
                "Read aromaticity model {}. Currently only "
454
                "OEAroModel_MDL is supported.".format(aromaticity_model)
455
            )
456

457 8
        self._aromaticity_model = aromaticity_model
458

459 8
    def _add_author(self, author):
460
        """
461
        Add an author to this forcefield. If this functional is called multiple times, all provided authors
462
        will be concatenated with the string " AND ". No redundancy checking is performed by this function.
463

464
        Parameters
465
        ----------
466
        author : str
467
            The author to add to this ForceField object
468
        """
469 8
        if self._author is None:
470 8
            self._author = author
471
        else:
472 8
            self._author += " AND " + author
473

474 8
    def _add_date(self, date):
475
        """
476
        Add an date to this forcefield. If this functional is called multiple times, all provided dates
477
        will be concatenated with the string " AND ". No redundancy checking is performed by this function.
478

479
        Parameters
480
        ----------
481
        date : str
482
            The author to add to this ForceField object
483
        """
484 8
        if self._date is None:
485 8
            self._date = date
486
        else:
487 8
            self._date += " AND " + date
488

489 8
    @property
490
    def author(self):
491
        """Returns the author data for this ForceField object. If not defined in any loaded files, this will be None.
492

493
        Returns
494
        -------
495
        author : str
496
            The author data for this forcefield.
497
        """
498 8
        return self._author
499

500 8
    @author.setter
501
    def author(self, author):
502
        """Set the author data for this ForceField object. If not defined in any loaded files, this will be None.
503

504
        Parameters
505
        ----------
506
        author : str
507
            The author data to set for this forcefield.
508
        """
509 8
        self._author = author
510

511 8
    @property
512
    def date(self):
513
        """Returns the date data for this ForceField object. If not defined in any loaded files, this will be None.
514

515
        Returns
516
        -------
517
        date : str
518
            The date data for this forcefield.
519
        """
520 8
        return self._date
521

522 8
    @date.setter
523
    def date(self, date):
524
        """Set the author data for this ForceField object. If not defined in any loaded files, this will be None.
525

526
        Parameters
527
        ----------
528
        date : str
529
            The date data to set for this forcefield.
530
        """
531 8
        self._date = date
532

533 8
    def _register_parameter_handler_classes(self, parameter_handler_classes):
534
        """
535
        Register multiple ParameterHandler classes, ensuring they specify unique tags to process
536

537
        .. warning :: This API is experimental and subject to change.
538

539
        Parameters
540
        ----------
541
        parameter_handler_classes : iterable of ParameterHandler subclasses
542
            List of ParameterHandler classes to register for this ForceField.
543
        """
544 8
        for parameter_handler_class in parameter_handler_classes:
545 8
            tagname = parameter_handler_class._TAGNAME
546 8
            if tagname is not None:
547 8
                if tagname in self._parameter_handler_classes:
548 0
                    raise Exception(
549
                        "Attempting to register ParameterHandler {}, which provides a parser for tag"
550
                        " '{}', but ParameterHandler {} has already been registered to handle that tag.".format(
551
                            parameter_handler_class,
552
                            tagname,
553
                            self._parameter_handler_classes[tagname],
554
                        )
555
                    )
556 8
                self._parameter_handler_classes[tagname] = parameter_handler_class
557

558 8
    def _register_parameter_io_handler_classes(self, parameter_io_handler_classes):
559
        """
560
        Register multiple ParameterIOHandler classes, ensuring they specify unique suffixes
561

562
        .. warning :: This API is experimental and subject to change.
563

564
        Parameters
565
        ----------
566
        parameter_io_handler_classes : iterable of ParameterIOHandler subclasses
567
            All specified ParameterIOHandler classes will be registered as ways to translate to/from the object model
568
            to serialized parameter sets.
569

570
        Raises
571
        ------
572
        Exception if two ParameterIOHandlers are attempted to be registered for the same file format.
573

574
        """
575 8
        for parameter_io_handler_class in parameter_io_handler_classes:
576 8
            serialization_format = parameter_io_handler_class._FORMAT
577 8
            if serialization_format is not None:
578 8
                if serialization_format in self._parameter_io_handler_classes.keys():
579 0
                    raise Exception(
580
                        "Attempting to register ParameterIOHandler {}, which provides a IO parser for format "
581
                        "'{}', but ParameterIOHandler {} has already been registered to handle that tag.".format(
582
                            parameter_io_handler_class,
583
                            serialization_format,
584
                            self._parameter_io_handler_classes[serialization_format],
585
                        )
586
                    )
587 8
                self._parameter_io_handler_classes[
588
                    serialization_format
589
                ] = parameter_io_handler_class
590

591 8
    def register_parameter_handler(self, parameter_handler):
592
        """
593
        Register a new ParameterHandler for a specific tag, making it
594
        available for lookup in the ForceField.
595

596
        .. warning :: This API is experimental and subject to change.
597

598
        Parameters
599
        ----------
600
        parameter_handler : A ParameterHandler object
601
            The ParameterHandler to register. The TAGNAME attribute of this object will be used as the key for
602
            registration.
603

604
        """
605 8
        tagname = parameter_handler._TAGNAME
606 8
        if tagname in self._parameter_handlers.keys():
607 0
            raise ParameterHandlerRegistrationError(
608
                "Tried to register parameter handler '{}' for tag '{}', but "
609
                "tag is already registered to {}".format(
610
                    parameter_handler, tagname, self._parameter_handlers[tagname]
611
                )
612
            )
613

614 8
        self._parameter_handlers[parameter_handler._TAGNAME] = parameter_handler
615

616 8
    def register_parameter_io_handler(self, parameter_io_handler):
617
        """
618
        Register a new ParameterIOHandler, making it available for lookup in the ForceField.
619

620
        .. warning :: This API is experimental and subject to change.
621

622
        Parameters
623
        ----------
624
        parameter_io_handler :  A ParameterIOHandler object
625
            The ParameterIOHandler to register. The FORMAT attribute of this object will be used
626
            to associate it to a file format/suffix.
627

628
        """
629 8
        io_format = parameter_io_handler._FORMAT
630 8
        if io_format in self._parameter_io_handlers.keys():
631 0
            raise ParameterHandlerRegistrationError(
632
                "Tried to register parameter IO handler '{}' for tag '{}', but "
633
                "tag is already registered to {}".format(
634
                    parameter_io_handler,
635
                    io_format,
636
                    self._parameter_io_handlers[io_format],
637
                )
638
            )
639 8
        self._parameter_io_handlers[io_format] = parameter_io_handler
640

641 8
    @property
642
    def registered_parameter_handlers(self):
643
        """
644
        Return the list of registered parameter handlers by name
645

646
        .. warning :: This API is experimental and subject to change.
647

648
        Returns
649
        -------
650
            registered_parameter_handlers: iterable of names of ParameterHandler objects in this ForceField
651

652
        """
653 8
        return [*self._parameter_handlers.keys()]
654

655
    # TODO: Do we want to make this optional?
656

657 8
    @staticmethod
658
    def _check_for_missing_valence_terms(
659
        name, topology, assigned_terms, topological_terms
660
    ):
661
        """
662
        Check to ensure there are no missing valence terms in the given topology, identifying potential gaps in
663
        parameter coverage.
664

665
        .. warning :: This API is experimental and subject to change.
666

667
        Parameters
668
        ----------
669
        name : str
670
            Name of the calling force Handler
671
        topology : openforcefield.topology.Topology
672
            The Topology object
673
        assigned_terms : iterable of ints or int tuples
674
            Atom index tuples defining added valence terms
675
        topological_terms : iterable of atoms or atom tuples
676
            Atom tuples defining topological valence atomsets to which forces should be assigned
677

678
        """
679
        # Convert assigned terms and topological terms to lists
680 0
        assigned_terms = [item for item in assigned_terms]
681 0
        topological_terms = [item for item in topological_terms]
682

683 0
        def ordered_tuple(atoms):
684 0
            atoms = list(atoms)
685 0
            if atoms[0] < atoms[-1]:
686 0
                return tuple(atoms)
687
            else:
688 0
                return tuple(reversed(atoms))
689

690 0
        try:
691 0
            topology_set = set(
692
                [
693
                    ordered_tuple(atom.index for atom in atomset)
694
                    for atomset in topological_terms
695
                ]
696
            )
697 0
            assigned_set = set(
698
                [
699
                    ordered_tuple(index for index in atomset)
700
                    for atomset in assigned_terms
701
                ]
702
            )
703 0
        except TypeError as te:
704 0
            topology_set = set([atom.index for atom in topological_terms])
705 0
            assigned_set = set([atomset[0] for atomset in assigned_terms])
706

707 0
        def render_atoms(atomsets):
708 0
            msg = ""
709 0
            for atomset in atomsets:
710 0
                msg += f"{atomset:30} :"
711 0
                try:
712 0
                    for atom_index in atomset:
713 0
                        atom = atoms[atom_index]
714 0
                        msg += f" {atom.residue.index:5} {atom.residue.name:3} {atom.name:3}"
715 0
                except TypeError as te:
716 0
                    atom = atoms[atomset]
717 0
                    msg += (
718
                        f" {atom.residue.index:5} {atom.residue.name:3} {atom.name:3}"
719
                    )
720

721 0
                msg += "\n"
722 0
            return msg
723

724 0
        if set(assigned_set) != set(topology_set):
725
            # Form informative error message
726 0
            msg = f"{name}: Mismatch between valence terms added and topological terms expected.\n"
727 0
            atoms = [atom for atom in topology.topology_atoms]
728 0
            if len(assigned_set.difference(topology_set)) > 0:
729 0
                msg += "Valence terms created that are not present in Topology:\n"
730 0
                msg += render_atoms(assigned_set.difference(topology_set))
731 0
            if len(topology_set.difference(assigned_set)) > 0:
732 0
                msg += "Topological atom sets not assigned parameters:\n"
733 0
                msg += render_atoms(topology_set.difference(assigned_set))
734 0
            msg += "topology_set:\n"
735 0
            msg += str(topology_set) + "\n"
736 0
            msg += "assigned_set:\n"
737 0
            msg += str(assigned_set) + "\n"
738 0
            raise Exception(
739
                msg
740
            )  # TODO: Should we raise a more specific exception here?
741

742 8
    def get_parameter_handler(
743
        self, tagname, handler_kwargs=None, allow_cosmetic_attributes=False
744
    ):
745
        """Retrieve the parameter handlers associated with the provided tagname.
746

747
        If the parameter handler has not yet been instantiated, it will be created and returned.
748
        If a parameter handler object already exists, it will be checked for compatibility
749
        and an Exception raised if it is incompatible with the provided kwargs. If compatible, the
750
        existing ParameterHandler will be returned.
751

752
        Parameters
753
        ----------
754
        tagname : str
755
            The name of the parameter to be handled.
756
        handler_kwargs : dict, optional. Default = None
757
            Dict to be passed to the handler for construction or checking compatibility. If this is None and no
758
            existing ParameterHandler exists for the desired tag, a handler will be initialized with all default
759
            values. If this is None and a handler for the desired tag exists, the existing ParameterHandler will
760
            be returned.
761
        allow_cosmetic_attributes : bool, optional. Default = False
762
            Whether to permit non-spec kwargs in smirnoff_data.
763

764
        Returns
765
        -------
766
        handler : An openforcefield.engines.typing.smirnoff.ParameterHandler
767

768
        Raises
769
        ------
770
        KeyError if there is no ParameterHandler for the given tagname
771
        """
772
        # If there are no kwargs for the handler, initialize handler_kwargs as an empty dict
773 8
        skip_version_check = False
774 8
        if handler_kwargs is None:
775 8
            handler_kwargs = dict()
776 8
            skip_version_check = True
777

778
        # Ensure that the ForceField has a ParameterHandler class registered that can handle this tag
779 8
        ph_class = self._get_parameter_handler_class(tagname)
780

781 8
        if tagname in self._parameter_handlers:
782
            # If a handler of this class already exists, ensure that the two handlers encode compatible science
783 8
            old_handler = self._parameter_handlers[tagname]
784
            # If no handler kwargs were provided, skip the compatibility check
785 8
            if handler_kwargs != {}:
786
                # Initialize a new instance of this parameter handler class with the given kwargs
787 8
                new_handler = ph_class(
788
                    **handler_kwargs,
789
                    allow_cosmetic_attributes=allow_cosmetic_attributes,
790
                    skip_version_check=skip_version_check,
791
                )
792 8
                old_handler.check_handler_compatibility(new_handler)
793 8
            return_handler = old_handler
794 8
        elif tagname in self._parameter_handler_classes:
795
            # Otherwise, register this handler in the forcefield
796
            # Initialize a new instance of this parameter handler class with the given kwargs
797 8
            new_handler = ph_class(
798
                **handler_kwargs,
799
                allow_cosmetic_attributes=allow_cosmetic_attributes,
800
                skip_version_check=skip_version_check,
801
            )
802 8
            self.register_parameter_handler(new_handler)
803 8
            return_handler = new_handler
804

805 8
        return return_handler
806

807 8
    def get_parameter_io_handler(self, io_format):
808
        """Retrieve the parameter handlers associated with the provided tagname.
809
        If the parameter IO handler has not yet been instantiated, it will be created.
810

811
        Parameters
812
        ----------
813
        io_format : str
814
            The name of the io format to be handled.
815

816
        Returns
817
        -------
818
        io_handler : An openforcefield.engines.typing.smirnoff.ParameterIOHandler
819

820
        Raises
821
        ------
822
        KeyError if there is no ParameterIOHandler for the given tagname
823
        """
824
        # Uppercase the format string to avoid case mismatches
825 8
        io_format = io_format.upper()
826

827
        # Remove "." (ParameterIOHandler tags do not include it)
828 8
        io_format = io_format.strip(".")
829

830
        # Find or initialize ParameterIOHandler for this format
831 8
        io_handler = None
832 8
        if io_format in self._parameter_io_handlers.keys():
833 8
            io_handler = self._parameter_io_handlers[io_format]
834 8
        elif io_format in self._parameter_io_handler_classes.keys():
835 8
            new_handler_class = self._parameter_io_handler_classes[io_format]
836 8
            io_handler = new_handler_class()
837 8
            self.register_parameter_io_handler(io_handler)
838 8
        if io_handler is None:
839 0
            msg = "Cannot find a registered parameter IO handler for format '{}'\n".format(
840
                io_format
841
            )
842 0
            msg += "Registered parameter IO handlers: {}\n".format(
843
                self._parameter_io_handlers.keys()
844
            )
845 0
            raise KeyError(msg)
846

847 8
        return io_handler
848

849 8
    def parse_sources(self, sources, allow_cosmetic_attributes=True):
850
        """Parse a SMIRNOFF force field definition.
851

852
        Parameters
853
        ----------
854
        sources : string or file-like object or open file handle or URL (or iterable of these)
855
            A list of files defining the SMIRNOFF force field to be loaded.
856
            Currently, only `the SMIRNOFF XML format <https://github.com/openforcefield/openforcefield/blob/master/The-SMIRNOFF-force-field-format.md>`_ is supported.
857
            Each entry may be an absolute file path, a path relative to the current working directory, a path relative to this module's data subdirectory
858
            (for built in force fields), or an open file-like object with a ``read()`` method from which the forcefield XML data can be loaded.
859
            If multiple files are specified, any top-level tags that are repeated will be merged if they are compatible,
860
            with files appearing later in the sequence resulting in parameters that have higher precedence.
861
            Support for multiple files is primarily intended to allow solvent parameters to be specified by listing them last in the sequence.
862
        allow_cosmetic_attributes : bool, optional. Default = False
863
            Whether to permit non-spec kwargs present in the source.
864

865
        .. notes ::
866
           * New SMIRNOFF sections are handled independently, as if they were specified in the same file.
867
           * If a SMIRNOFF section that has already been read appears again, its definitions are appended to the end of the previously-read
868
             definitions if the sections are configured with compatible attributes; otherwise, an ``IncompatibleTagException`` is raised.
869

870
        """
871
        # Ensure that we are working with an iterable
872 8
        try:
873 8
            sources = iter(sources)
874 0
        except TypeError as te:
875
            # Make iterable object
876 0
            sources = [sources]
877

878
        # TODO: If a non-first source fails here, the forcefield might be partially modified
879 8
        for source in sources:
880 8
            smirnoff_data = self.parse_smirnoff_from_source(source)
881 8
            self._load_smirnoff_data(
882
                smirnoff_data, allow_cosmetic_attributes=allow_cosmetic_attributes
883
            )
884

885 8
    def _to_smirnoff_data(self, discard_cosmetic_attributes=False):
886
        """
887
        Convert this ForceField and all related ParameterHandlers to an OrderedDict representing a SMIRNOFF
888
        data object.
889

890
        Returns
891
        -------
892
        smirnoff_dict : OrderedDict
893
            A nested OrderedDict representing this ForceField as a SMIRNOFF data object.
894
        discard_cosmetic_attributes : bool, optional. Default=False
895
            Whether to discard any non-spec attributes stored in the ForceField.
896

897
        """
898 8
        l1_dict = OrderedDict()
899

900
        # Assume we will write out SMIRNOFF data in compliance with the max supported spec version
901 8
        l1_dict["version"] = self._MAX_SUPPORTED_SMIRNOFF_VERSION
902

903
        # Write out the aromaticity model used
904 8
        l1_dict["aromaticity_model"] = self._aromaticity_model
905

906
        # Write out author and date (if they have been set)
907 8
        if not (self._author is None):
908 8
            l1_dict["Author"] = self._author
909

910
        # Write out author and date (if they have been set)
911 8
        if not (self._date is None):
912 8
            l1_dict["Date"] = self._date
913

914 8
        for handler_format, parameter_handler in self._parameter_handlers.items():
915 8
            handler_tag = parameter_handler._TAGNAME
916 8
            l1_dict[handler_tag] = parameter_handler.to_dict(
917
                discard_cosmetic_attributes=discard_cosmetic_attributes
918
            )
919

920 8
        smirnoff_dict = OrderedDict()
921 8
        smirnoff_dict["SMIRNOFF"] = l1_dict
922 8
        smirnoff_dict = convert_all_quantities_to_string(smirnoff_dict)
923 8
        return smirnoff_dict
924

925
    # TODO: Should we call this "from_dict"?
926 8
    def _load_smirnoff_data(self, smirnoff_data, allow_cosmetic_attributes=False):
927
        """
928
        Add parameters from a SMIRNOFF-format data structure to this ForceField.
929

930
        Parameters
931
        ----------
932
        smirnoff_data : OrderedDict
933
            A representation of a SMIRNOFF-format data structure. Begins at top-level 'SMIRNOFF' key.
934
        allow_cosmetic_attributes : bool, optional. Default = False
935
            Whether to permit non-spec kwargs in smirnoff_data.
936
        """
937 8
        import packaging.version
938

939
        # Check that the SMIRNOFF version of this data structure is supported by this ForceField implementation
940

941 8
        if "SMIRNOFF" in smirnoff_data:
942 8
            version = smirnoff_data["SMIRNOFF"]["version"]
943 8
        elif "SMIRFF" in smirnoff_data:
944 8
            version = smirnoff_data["SMIRFF"]["version"]
945
        else:
946 0
            raise ParseError("'version' attribute must be specified in SMIRNOFF tag")
947

948 8
        self._check_smirnoff_version_compatibility(str(version))
949
        # Convert 0.1 spec files to 0.3 SMIRNOFF data format by converting
950
        # from 0.1 spec to 0.2, then 0.2 to 0.3
951 8
        if packaging.version.parse(str(version)) == packaging.version.parse("0.1"):
952
            # NOTE: This will convert the top-level "SMIRFF" tag to "SMIRNOFF"
953 8
            smirnoff_data = convert_0_1_smirnoff_to_0_2(smirnoff_data)
954 8
            smirnoff_data = convert_0_2_smirnoff_to_0_3(smirnoff_data)
955

956
        # Convert 0.2 spec files to 0.3 SMIRNOFF data format by removing units
957
        # from section headers and adding them to quantity strings at all levels.
958 8
        elif packaging.version.parse(str(version)) == packaging.version.parse("0.2"):
959 8
            smirnoff_data = convert_0_2_smirnoff_to_0_3(smirnoff_data)
960

961
        # Ensure that SMIRNOFF is a top-level key of the dict
962 8
        if not ("SMIRNOFF" in smirnoff_data):
963 0
            raise ParseError(
964
                "'SMIRNOFF' must be a top-level key in the SMIRNOFF object model"
965
            )
966

967
        # Check that the aromaticity model required by this parameter set is compatible with
968
        # others loaded by this ForceField
969 8
        if "aromaticity_model" in smirnoff_data["SMIRNOFF"]:
970 8
            aromaticity_model = smirnoff_data["SMIRNOFF"]["aromaticity_model"]
971 8
            self.aromaticity_model = aromaticity_model
972

973 0
        elif self._aromaticity_model is None:
974 0
            raise ParseError(
975
                "'aromaticity_model' attribute must be specified in SMIRNOFF "
976
                "tag, or contained in a previously-loaded SMIRNOFF data source"
977
            )
978

979 8
        if "Author" in smirnoff_data["SMIRNOFF"]:
980 8
            self._add_author(smirnoff_data["SMIRNOFF"]["Author"])
981

982 8
        if "Date" in smirnoff_data["SMIRNOFF"]:
983 8
            self._add_date(smirnoff_data["SMIRNOFF"]["Date"])
984

985
        # Go through the whole SMIRNOFF data structure, trying to convert all strings to Quantity
986 8
        smirnoff_data = convert_all_strings_to_quantity(smirnoff_data)
987

988
        # Go through the subsections, delegating each to the proper ParameterHandler
989

990
        # Define keys which are expected from the spec, but are not parameter sections
991 8
        l1_spec_keys = ["Author", "Date", "version", "aromaticity_model"]
992
        # TODO: Throw SMIRNOFFSpecError for unrecognized keywords
993

994 8
        for parameter_name in smirnoff_data["SMIRNOFF"]:
995
            # Skip (for now) cosmetic l1 items. They're handled above
996 8
            if parameter_name in l1_spec_keys:
997 8
                continue
998
            # Handle cases where a parameter name has no info (eg. ToolkitAM1BCC)
999 8
            if smirnoff_data["SMIRNOFF"][parameter_name] is None:
1000
                # "get"ting the parameter handler here will also initialize it.
1001 8
                _ = self.get_parameter_handler(parameter_name, {})
1002 0
                continue
1003

1004
            # Otherwise, we expect this l1_key to correspond to a ParameterHandler
1005 8
            section_dict = smirnoff_data["SMIRNOFF"][parameter_name]
1006

1007
            # TODO: Implement a ParameterHandler.from_dict() that knows how to deserialize itself for extensibility.
1008
            #       We could let it load the ParameterTypes from the dict and append them to the existing handler
1009
            #       after verifying that they are compatible.
1010

1011
            # Get the parameter types serialization that is not passed to the ParameterHandler constructor.
1012 8
            ph_class = self._get_parameter_handler_class(parameter_name)
1013 8
            try:
1014 8
                parameter_list_tagname = ph_class._INFOTYPE._ELEMENT_NAME
1015 8
            except AttributeError:
1016
                # The ParameterHandler doesn't have ParameterTypes (e.g. ToolkitAM1BCCHandler).
1017 8
                parameter_list_dict = {}
1018
            else:
1019 8
                parameter_list_dict = section_dict.pop(parameter_list_tagname, {})
1020

1021
                # If the parameter list isn't empty, it must be transferred into its own tag.
1022
                # This is necessary for deserializing SMIRNOFF force field sections which may or may
1023
                # not have any smirks-based elements (like an empty ChargeIncrementModel section)
1024 8
                if parameter_list_dict != {}:
1025 8
                    parameter_list_dict = {parameter_list_tagname: parameter_list_dict}
1026

1027
            # Retrieve or create parameter handler, passing in section_dict to check for
1028
            # compatibility if a handler for this parameter name already exists
1029 8
            handler = self.get_parameter_handler(
1030
                parameter_name,
1031
                section_dict,
1032
                allow_cosmetic_attributes=allow_cosmetic_attributes,
1033
            )
1034 8
            handler._add_parameters(
1035
                parameter_list_dict, allow_cosmetic_attributes=allow_cosmetic_attributes
1036
            )
1037

1038 8
    def parse_smirnoff_from_source(self, source):
1039
        """
1040
        Reads a SMIRNOFF data structure from a source, which can be one of many types.
1041

1042
        Parameters
1043
        ----------
1044
        source : str or bytes
1045
            sources : string or file-like object or open file handle or URL (or iterable of these)
1046
            A list of files defining the SMIRNOFF force field to be loaded
1047
            Currently, only `the SMIRNOFF XML format <https://github.com/openforcefield/openforcefield/blob/master/The-SMIRNOFF-force-field-format.md>`_ is supported.
1048
            Each entry may be an absolute file path, a path relative to the current working directory, a path relative to this module's data subdirectory
1049
            (for built in force fields), or an open file-like object with a ``read()`` method from which the forcefield XML data can be loaded.
1050

1051
        Returns
1052
        -------
1053
        smirnoff_data : OrderedDict
1054
            A representation of a SMIRNOFF-format data structure. Begins at top-level 'SMIRNOFF' key.
1055

1056
        """
1057 8
        from openforcefield.utils import get_data_file_path
1058

1059
        # Check whether this could be a file path. It could also be a
1060
        # file handler or a simple XML string.
1061 8
        if isinstance(source, str):
1062
            # Try first the simple path.
1063 8
            searched_dirs_paths = [""]
1064
            # Then try a relative file path w.r.t. an installed directory.
1065 8
            searched_dirs_paths.extend(_get_installed_offxml_dir_paths())
1066
            # Finally, search in openforcefield/data/.
1067
            # TODO: Remove this when smirnoff99Frosst 1.0.9 will be released.
1068 8
            searched_dirs_paths.append(get_data_file_path(""))
1069

1070
            # Determine the actual path of the file.
1071
            # TODO: What is desired toolkit behavior if two files with the desired name are available?
1072 8
            for dir_path in searched_dirs_paths:
1073 8
                file_path = os.path.join(dir_path, source)
1074 8
                if os.path.isfile(file_path):
1075 8
                    source = file_path
1076 8
                    break
1077

1078
        # Process all SMIRNOFF definition files or objects
1079
        # QUESTION: Allow users to specify forcefield URLs so they can pull forcefield definitions from the web too?
1080 8
        io_formats_to_try = self._parameter_io_handler_classes.keys()
1081

1082
        # Parse content depending on type
1083 8
        for parameter_io_format in io_formats_to_try:
1084 8
            parameter_io_handler = self.get_parameter_io_handler(parameter_io_format)
1085

1086
            # Try parsing as a forcefield file or file-like object
1087 8
            try:
1088 8
                smirnoff_data = parameter_io_handler.parse_file(source)
1089 8
                return smirnoff_data
1090 8
            except ParseError as e:
1091 0
                exception_msg = e.msg
1092 8
            except (FileNotFoundError, OSError):
1093
                # If this is not a file path or a file handle, attempt parsing as a string.
1094 8
                try:
1095 8
                    smirnoff_data = parameter_io_handler.parse_string(source)
1096 8
                    return smirnoff_data
1097 8
                except ParseError as e:
1098 8
                    exception_msg = e.msg
1099

1100
        # If we haven't returned by now, the parsing was unsuccessful
1101 8
        valid_formats = [
1102
            input_format for input_format in self._parameter_io_handlers.keys()
1103
        ]
1104 8
        msg = f"Source {source} could not be read. If this is a file, ensure that the path is correct.\n"
1105 8
        msg += "If the file is present, ensure it is in a known SMIRNOFF encoding.\n"
1106 8
        msg += f"Valid formats are: {valid_formats}\n"
1107 8
        msg += f"Parsing failed with the following error:\n{exception_msg}\n"
1108 8
        raise IOError(msg)
1109

1110 8
    def to_string(self, io_format="XML", discard_cosmetic_attributes=False):
1111
        """
1112
        Write this Forcefield and all its associated parameters to a string in a given format which
1113
        complies with the SMIRNOFF spec.
1114

1115

1116
        Parameters
1117
        ----------
1118
        io_format : str or ParameterIOHandler, optional. Default='XML'
1119
            The serialization format to write to
1120
        discard_cosmetic_attributes : bool, default=False
1121
            Whether to discard any non-spec attributes stored in the ForceField.
1122

1123
        Returns
1124
        -------
1125
        forcefield_string : str
1126
            The string representation of the serialized forcefield
1127
        """
1128
        # Resolve which IO handler to use
1129 8
        if isinstance(io_format, ParameterIOHandler):
1130 0
            io_handler = io_format
1131
        else:
1132 8
            io_handler = self.get_parameter_io_handler(io_format)
1133

1134 8
        smirnoff_data = self._to_smirnoff_data(
1135
            discard_cosmetic_attributes=discard_cosmetic_attributes
1136
        )
1137 8
        string_data = io_handler.to_string(smirnoff_data)
1138 8
        return string_data
1139

1140 8
    def to_file(self, filename, io_format=None, discard_cosmetic_attributes=False):
1141
        """
1142
        Write this Forcefield and all its associated parameters to a string in a given format which
1143
        complies with the SMIRNOFF spec.
1144

1145

1146
        Parameters
1147
        ----------
1148
        filename : str
1149
            The filename to write to
1150
        io_format : str or ParameterIOHandler, optional. Default=None
1151
            The serialization format to write out. If None, will attempt to be inferred from the filename.
1152
        discard_cosmetic_attributes : bool, default=False
1153
            Whether to discard any non-spec attributes stored in the ForceField.
1154

1155
        Returns
1156
        -------
1157
        forcefield_string : str
1158
            The string representation of the serialized forcefield
1159
        """
1160

1161 8
        if io_format is None:
1162 8
            basename, io_format = os.path.splitext(filename)
1163

1164
        # Resolve which IO handler to use
1165 8
        if isinstance(io_format, ParameterIOHandler):
1166 8
            io_handler = io_format
1167
        else:
1168
            # Handle the fact that .offxml is the same as .xml
1169 8
            if io_format.lower() == "offxml" or io_format.lower() == ".offxml":
1170 8
                io_format = "xml"
1171 8
            io_handler = self.get_parameter_io_handler(io_format)
1172

1173
        # Write out the SMIRNOFF data to the IOHandler
1174 8
        smirnoff_data = self._to_smirnoff_data(
1175
            discard_cosmetic_attributes=discard_cosmetic_attributes
1176
        )
1177 8
        io_handler.to_file(filename, smirnoff_data)
1178

1179 8
    def _resolve_parameter_handler_order(self):
1180
        """Resolve the order in which ParameterHandler objects should execute to satisfy constraints.
1181

1182
        Returns
1183
        -------
1184
        Iterable of ParameterHandlers
1185
            The ParameterHandlers in this ForceField, in the order that they should be called to satisfy constraints.
1186
        """
1187

1188
        # Create a DAG expressing dependencies
1189 8
        import networkx as nx
1190

1191 8
        G = nx.DiGraph()
1192 8
        for tagname, parameter_handler in self._parameter_handlers.items():
1193 8
            G.add_node(tagname)
1194 8
        for tagname, parameter_handler in self._parameter_handlers.items():
1195 8
            if parameter_handler._DEPENDENCIES is not None:
1196 8
                for dependency in parameter_handler._DEPENDENCIES:
1197 8
                    G.add_edge(dependency._TAGNAME, parameter_handler._TAGNAME)
1198

1199
        # Ensure there are no loops in handler order
1200 8
        if not (nx.is_directed_acyclic_graph(G)):
1201 8
            raise RuntimeError(
1202
                "Unable to resolve order in which to run ParameterHandlers. Dependencies do not form "
1203
                "a directed acyclic graph."
1204
            )
1205
        # Resolve order
1206 8
        ordered_parameter_handlers = list()
1207 8
        for tagname in nx.topological_sort(G):
1208 8
            if tagname in self._parameter_handlers:
1209 8
                ordered_parameter_handlers.append(self._parameter_handlers[tagname])
1210 8
        return ordered_parameter_handlers
1211

1212
    # TODO: Should we add convenience methods to parameterize a Topology and export directly to AMBER, gromacs, CHARMM, etc.?
1213
    #       Or should we create an "enhanced" openforcefield System object that knows how to convert to all of these formats?
1214
    #       We could even create a universal applyParameters(format='AMBER') method that allows us to export to whatever system we want.
1215

1216
    # TODO: Should the Topology contain the default box vectors? Or should we require they be specified externally?
1217

1218
    # TODO: How do we know if the system is periodic or not?
1219
    # TODO: Should we also accept a Molecule as an alternative to a Topology?
1220

1221 8
    def create_openmm_system(self, topology, **kwargs):
1222
        """Create an OpenMM System representing the interactions for the specified Topology with the current force field
1223

1224
        Parameters
1225
        ----------
1226
        topology : openforcefield.topology.Topology
1227
            The ``Topology`` corresponding to the system to be parameterized
1228
        charge_from_molecules : List[openforcefield.molecule.Molecule], optional. default =[]
1229
            If specified, partial charges will be taken from the given molecules
1230
            instead of being determined by the force field.
1231
        partial_bond_orders_from_molecules : List[openforcefield.molecule.Molecule], optional. default=[]
1232
            If specified, partial bond orders will be taken from the given molecules
1233
            instead of being determined by the force field.
1234
            **All** bonds on each molecule given must have ``fractional_bond_order`` specified.
1235
            A `ValueError` will be raised if any bonds have ``fractional_bond_order=None``.
1236
            Molecules in the topology not represented in this list will have fractional
1237
            bond orders calculated using underlying toolkits as needed.
1238
        return_topology : bool, optional. default=False
1239
            If ``True``, return tuple of ``(system, topology)``, where
1240
            ``topology`` is the processed topology. Default ``False``. This topology will have the
1241
            final partial charges assigned on its reference_molecules attribute, as well as partial
1242
            bond orders (if they were calculated).
1243

1244
        Returns
1245
        -------
1246
        system : simtk.openmm.System
1247
            The newly created OpenMM System corresponding to the specified ``topology``
1248
        topology : openforcefield.topology.Topology, optional.
1249
            If the `return_topology` keyword argument is used, this method will also return a Topology. This
1250
            can be used to inspect the partial charges and partial bond orders assigned to the molecules
1251
            during parameterization.
1252

1253
        """
1254 8
        return_topology = kwargs.pop("return_topology", False)
1255

1256
        # Make a deep copy of the topology so we don't accidentally modify it
1257 8
        topology = copy.deepcopy(topology)
1258

1259
        # set all fractional_bond_orders in topology to None
1260 8
        for ref_mol in topology.reference_molecules:
1261 8
            for bond in ref_mol.bonds:
1262 8
                bond.fractional_bond_order = None
1263

1264
        # Set the topology aromaticity model to that used by the current forcefield
1265
        # TODO: See openforcefield issue #206 for proposed implementation of aromaticity
1266
        # topology.set_aromaticity_model(self._aromaticity_model)
1267

1268
        # Create an empty OpenMM System
1269 8
        system = openmm.System()
1270

1271
        # Set periodic boundary conditions if specified
1272 8
        if topology.box_vectors is not None:
1273 8
            system.setDefaultPeriodicBoxVectors(*topology.box_vectors)
1274

1275
        # Add particles (both atoms and virtual sites) with appropriate masses
1276 8
        for atom in topology.topology_particles:
1277 8
            system.addParticle(atom.atom.mass)
1278

1279
        # Determine the order in which to process ParameterHandler objects in order to satisfy dependencies
1280 8
        parameter_handlers = self._resolve_parameter_handler_order()
1281

1282
        # Check if any kwargs have been provided that aren't handled by force Handlers
1283
        # TODO: Delete this and kwargs from arguments above?
1284 8
        known_kwargs = set()
1285 8
        for parameter_handler in parameter_handlers:
1286 8
            known_kwargs.update(parameter_handler.known_kwargs)
1287 8
        unknown_kwargs = set(kwargs.keys()).difference(known_kwargs)
1288 8
        if len(unknown_kwargs) > 0:
1289 8
            msg = "The following keyword arguments to create_openmm_system() are not used by any registered force Handler: {}\n".format(
1290
                unknown_kwargs
1291
            )
1292 8
            msg += "Known keyword arguments: {}".format(known_kwargs)
1293 8
            raise ValueError(msg)
1294

1295
        # Add forces and parameters to the System
1296 8
        for parameter_handler in parameter_handlers:
1297 8
            parameter_handler.create_force(system, topology, **kwargs)
1298

1299
        # Let force Handlers do postprocessing
1300 8
        for parameter_handler in parameter_handlers:
1301 8
            parameter_handler.postprocess_system(system, topology, **kwargs)
1302

1303 8
        if return_topology:
1304 8
            return (system, topology)
1305
        else:
1306 8
            return system
1307

1308 8
    def create_parmed_structure(self, topology, positions, **kwargs):
1309
        """Create a ParmEd Structure object representing the interactions for the specified Topology with the current force field
1310

1311
        This method creates a `ParmEd <http://github.com/parmed/parmed>`_ ``Structure`` object containing a topology, positions, and parameters.
1312

1313
        Parameters
1314
        ----------
1315
        topology : openforcefield.topology.Topology
1316
            The ``Topology`` corresponding to the ``System`` object to be created.
1317
        positions : simtk.unit.Quantity of dimension (natoms,3) with units compatible with angstroms
1318
            The positions corresponding to the ``System`` object to be created
1319

1320
        Returns
1321
        -------
1322
        structure : parmed.Structure
1323
            The newly created ``parmed.Structure`` object
1324

1325
        """
1326 0
        raise NotImplementedError
1327
        # import parmed
1328
        # TODO: Automagically handle expansion of virtual sites? Or is Topology supposed to do that?
1329

1330
        # Create OpenMM System
1331
        # system = self.create_openmm_system(
1332
        #    topology, **kwargs)
1333

1334
        # Create a ParmEd Structure object
1335
        # structure = parmed.openmm.topsystem.load_topology(
1336
        #    topology.to_openmm(), system, positions)
1337
        #
1338
        # return structure
1339

1340 8
    def label_molecules(self, topology):
1341
        """Return labels for a list of molecules corresponding to parameters from this force field.
1342
        For each molecule, a dictionary of force types is returned, and for each force type,
1343
        each force term is provided with the atoms involved, the parameter id assigned, and the corresponding SMIRKS.
1344

1345
        Parameters
1346
        ----------
1347
        topology : openforcefield.topology.Topology
1348
            A Topology object containing one or more unique molecules to be labeled
1349

1350
        Returns
1351
        -------
1352
        molecule_labels : list
1353
            List of labels for unique molecules. Each entry in the list corresponds to
1354
            one unique molecule in the Topology and is a dictionary keyed by force type,
1355
            i.e., ``molecule_labels[0]['HarmonicBondForce']`` gives details for the harmonic
1356
            bond parameters for the first molecule. Each element is a list of the form:
1357
            ``[ ( [ atom1, ..., atomN], parameter_id, SMIRKS), ... ]``.
1358

1359
        .. todo ::
1360

1361
           What is the most useful API for this method?
1362
           Should we instead accept :class:`Molecule` objects as input and individually return labels?
1363
           Should we attach the labels to the :class:`Molecule` object?
1364
           Or should we label all interactions in a :class:`Topology` instead of just labeling its ``unique_molecules``?
1365

1366
        """
1367 0
        from openforcefield.topology import Topology
1368

1369
        # Loop over molecules and label
1370 0
        molecule_labels = list()
1371 0
        for molecule_idx, molecule in enumerate(topology.reference_molecules):
1372 0
            top_mol = Topology.from_molecules([molecule])
1373 0
            current_molecule_labels = dict()
1374 0
            for tag, parameter_handler in self._parameter_handlers.items():
1375

1376 0
                matches = parameter_handler.find_matches(top_mol)
1377

1378
                # Remove the chemical environment matches from the
1379
                # matched results.
1380

1381
                # Because we sometimes need to enforce atom ordering,
1382
                # `matches` here is not a normal `dict`, but rather
1383
                # one that transforms keys in arbitrary ways. Thus,
1384
                # we need to make a copy of its specific class here.
1385

1386 0
                parameter_matches = matches.__class__()
1387

1388
                # Now make parameter_matches into a dict mapping
1389
                # match objects to ParameterTypes
1390

1391 0
                for match in matches:
1392 0
                    parameter_matches[match] = matches[match].parameter_type
1393

1394 0
                current_molecule_labels[tag] = parameter_matches
1395

1396 0
            molecule_labels.append(current_molecule_labels)
1397 0
        return molecule_labels
1398

1399 8
    def _get_parameter_handler_class(self, tagname):
1400
        """Retrieve the ParameterHandler class associated to the tagname and throw a custom error if not found."""
1401 8
        try:
1402 8
            ph_class = self._parameter_handler_classes[tagname]
1403 8
        except KeyError:
1404 8
            msg = "Cannot find a registered parameter handler class for tag '{}'\n".format(
1405
                tagname
1406
            )
1407 8
            msg += "Known parameter handler class tags are {}".format(
1408
                self._parameter_handler_classes.keys()
1409
            )
1410 8
            raise KeyError(msg)
1411 8
        return ph_class
1412

1413 8
    def __getitem__(self, val):
1414
        """
1415
        Syntax sugar for lookikng up a ParameterHandler. Note that only
1416
        string-based lookups are currently supported.
1417
        """
1418 8
        if isinstance(val, str):
1419 8
            if val in self._parameter_handlers:
1420 8
                return self.get_parameter_handler(val)
1421
            else:
1422 8
                raise KeyError(f"Parameter handler with name '{val}' not found.")
1423 8
        elif isinstance(val, ParameterHandler) or issubclass(val, ParameterHandler):
1424 8
            raise NotImplementedError
1425

1426 8
    def __hash__(self):
1427
        """Deterministically hash a ForceField object
1428

1429
        Notable behavior:
1430
          * `author` and `date` are stripped from the ForceField
1431
          * `id` and `parent_id` are stripped from each ParameterType"""
1432

1433
        # Completely re-constructing the force field may be overkill
1434
        # compared to deepcopying and modifying, but is not currently slow
1435 8
        ff_copy = ForceField()
1436 8
        ff_copy.date = None
1437 8
        ff_copy.author = None
1438

1439 8
        param_attrs_to_strip = ["_id", "_parent_id"]
1440

1441 8
        for handler_name in self.registered_parameter_handlers:
1442 8
            handler = copy.deepcopy(self.get_parameter_handler(handler_name))
1443

1444 8
            for param in handler._parameters:
1445 8
                for attr in param_attrs_to_strip:
1446
                    # param.__dict__.pop(attr, None) may be faster
1447
                    # https://stackoverflow.com/a/42303681/4248961
1448 8
                    if hasattr(param, attr):
1449 8
                        delattr(param, attr)
1450

1451 8
            ff_copy.register_parameter_handler(handler)
1452

1453 8
        return hash(ff_copy.to_string(discard_cosmetic_attributes=True))

Read our documentation on viewing source code .

Loading