Merge pull request #36 from molssi-seamm/dev

Bugfix: factor of 2 for dihedrals and impropers in ligpargen
molssi-seamm · Jun 29, 2024 · 5b2eaa9 · 5b2eaa9
2 parents d1598eb + f031b50
commit 5b2eaa9
Show file tree

Hide file tree

Showing 14 changed files with 337 additions and 30 deletions.
diff --git a/HISTORY.rst b/HISTORY.rst
@@ -1,6 +1,11 @@
 =======
 History
 =======
+2024.6.29 -- Bugfix: factor of 2 for dihedrals and impropers in ligpargen
+  * The ligpargen tool was missing a factor of 2 in the dihedral and improper parameters.
+  * Corrected the search paths for forcefields.
+  * Improved the documentation.
+
 2024.6.28 -- Added customizable local forcefields and LigParGen interface.
   * Added the machinery to handle local forcefield files in either
     ~/.seamm.d/data/Forcefields (personal) or ~/SEAMM/data/Forcefields (site).

diff --git a/docs/getting_started/images/dialog.png b/docs/getting_started/images/dialog.png
diff --git a/docs/getting_started/images/ff_choices.png b/docs/getting_started/images/ff_choices.png
diff --git a/docs/getting_started/images/initial.png b/docs/getting_started/images/initial.png
diff --git a/docs/getting_started/index.rst b/docs/getting_started/index.rst
@@ -22,8 +22,79 @@ will ensure both that it is installed and up-to-date.
 
 .. _SEAMM Installer: https://molssi-seamm.github.io/installation/index.html
 
-Replace this!
-=============
-Put an example or two here....
+Intoduction to the Forcefield Step
+==================================
+The Forcefield Step is used to read a forcefield into SEAMM and set it as the default
+forcefield until another Forcefield Step changes the default. You can also use the
+Forcefield Step to assign the forcefield to your system, though the assignment will also
+be done automatically when needed so you do not need to explicitly ask the Forcefield
+Step to do the assignment.
 
-That should be enough to get started. For more detail about the functionality in this plug-in, see the :ref:`User Guide <user-guide>`.
+The Forcefield Step resides in the Simulation section of the Flowchart. Typically it is
+one of the firt couple steps in a Flowchart, defining the forcefield to be used by
+subsequent steps:
+
+.. figure:: images/initial.png
+   :align: center
+   :alt: A flowchart with a Forcefield Step
+
+   A flowchart with a single Forcefield Step
+
+Opening the Forcefield Step gives this dialog:
+
+.. figure:: images/dialog.png
+   :align: center
+   :alt: The initial forcefield dialog
+
+   The initial forcefield dialog
+
+The first field, "What to do:", allows you to setup a forcefield or to assign it to the
+current system, as mentioned above. The second, "Forcefield Repository", is used to
+select the type of forcefield to setup:
+
+.. figure:: images/ff_choices.png
+   :align: center
+   :alt: The forcefield choices
+
+   The choice of forcefield to set up.
+
+There are currently four main possibilities:
+
+   #. OpenKIM_ is an online repository that contains many of the "interatomic potentials"
+      used for materials simulations. These are potentials such as the Stillinger-Webber
+      potential for silicon, the potentials for the embedded atom method (EAM), and
+      those for the modified embedded atom method (MEAM).
+
+   #. OPLS-AA_, the Optimized Potentials for Liquid Simulations -- All Atoms forcefield
+      of William Jorgensen and et al. [#]_, has wide coverage of organic and
+      biomolecular systems.
+
+   #. PCFF2018_, the Polymer Consistent Force Field [#]_,also covers organic molecules
+      with an emphasis on (synthetic) polymers.
+
+   #. local:OPLS_AA/ligpargen.frc contains parameters that the user has downloaded from
+      the LigParGen server. See the :ref:`LigParGen` section in the :ref:`User Guide
+      <user-guide>` for more information.
+
+Your installation may have a different set of forcefields than shown above. Once you
+have setup a forcefield you can proceed to e.g. LAMMPS to run the simulation using the
+forcefield.
+
+That should be enough to get started. For more detail about the functionality in this
+plug-in, see the :ref:`User Guide <user-guide>`.
+
+
+.. [#] W.L. Jorgensen, J. Tirado-Rives,
+       Potential energy functions for atomic-level simulations of water and organic and
+       biomolecular systems,
+       Proc. Natl. Acad. Sci. 102 (2005) 6665–6670.
+       https://doi.org/10.1073/pnas.0408037102.
+       
+.. [#] Huai Sun, Stephen J. Mumby, Jon R. Maple, and Arnold T. Hagler,
+       An ab Initio CFF93 All-Atom Force Field for Polycarbonates
+       Journal of the American Chemical Society 1994 116 (7), 2978-2987
+       DOI: 10.1021/ja00086a030
+
+.. _OPLS-AA: https://doi.org/10.1073/pnas.0408037102
+.. _OpenKIM: https://www.openkim.org
+.. _PCFF2018: https://pubs.acs.org/doi/10.1021/ja00086a030
diff --git a/docs/user_guide/index.rst b/docs/user_guide/index.rst
@@ -3,16 +3,28 @@
 **********
 User Guide
 **********
-The XXXX plug-in ...
+The Forcefield plug-in manages one or more forcefields or sets of interatomic potentials
+for your simulations. It reads a forcefield file, selects a particular forcefield from
+the file, and prepares the forcefield parameters for susbsequent simulation steps. This
+step can also be used to assign the forcefield to the molecule or system that you are
+working on, though that happens by default when you start the simulation, so you don't
+need to handle it explicitly.
 
-..
-   The following sections cover accessing and controlling this functionality.
+If you want to understand more about the theory and definition of forcefields, see the
+`Overview of Forcefields`_ in the main SEAMM documentation.
 
-   .. toctree::
-      :maxdepth: 2
-      :titlesonly:
+The following sections cover accessing and controlling forcefields.
+
+.. toctree::
+   :maxdepth: 2
+   :titlesonly:
+
+   ligpargen/index
 
 Index
 =====
 
 * :ref:`genindex`
+
+
+.. _Overview of Forcefields:  https://molssi-seamm.github.io/background/forcefields/overview.html
diff --git a/docs/user_guide/ligpargen/images/ethane.key b/docs/user_guide/ligpargen/images/ethane.key
@@ -0,0 +1,160 @@
+
+
+      ##############################
+      ##                          ##
+      ##  Force Field Definition  ##
+      ##                          ##
+      ##############################
+
+
+forcefield              OPLS-AA
+
+vdwindex                TYPE
+vdwtype                 LENNARD-JONES
+radiusrule              GEOMETRIC
+radiustype              SIGMA
+radiussize              DIAMETER
+epsilonrule             GEOMETRIC
+torsionunit             1.0
+imptorunit              1.0
+vdw-14-scale            2.0
+chg-14-scale            2.0
+electric                332.06
+dielectric              1.0
+
+
+      #############################
+      ##                         ##
+      ##  Atom Type Definitions  ##
+      ##                         ##
+      #############################
+
+
+atom        800  800    CT    "C00"          6     12.011     4 
+atom        801  801    CT    "C01"          6     12.011     4 
+atom        802  802    HC    "H02"          1      1.008     1 
+atom        803  803    HC    "H03"          1      1.008     1 
+atom        804  804    HC    "H04"          1      1.008     1 
+atom        805  805    HC    "H05"          1      1.008     1 
+atom        806  806    HC    "H06"          1      1.008     1 
+atom        807  807    HC    "H07"          1      1.008     1 
+
+
+
+      ################################
+      ##                            ##
+      ##  Van der Waals Parameters  ##
+      ##                            ##
+      ################################
+
+
+vdw         800           3.5000   0.0660 
+vdw         801           3.5000   0.0660 
+vdw         802           2.5000   0.0300 
+vdw         803           2.5000   0.0300 
+vdw         804           2.5000   0.0300 
+vdw         805           2.5000   0.0300 
+vdw         806           2.5000   0.0300 
+vdw         807           2.5000   0.0300 
+
+
+
+      ##################################
+      ##                              ##
+      ##  Bond Stretching Parameters  ##
+      ##                              ##
+      ##################################
+
+
+bond        801  800           268.00   1.5290 
+bond        802  800           340.00   1.0900 
+bond        803  800           340.00   1.0900 
+bond        804  800           340.00   1.0900 
+bond        805  801           340.00   1.0900 
+bond        806  801           340.00   1.0900 
+bond        807  801           340.00   1.0900 
+
+
+
+      ################################
+      ##                            ##
+      ##  Angle Bending Parameters  ##
+      ##                            ##
+      ################################
+
+
+angle       801  800  802    37.50   110.70 
+angle       801  800  803    37.50   110.70 
+angle       801  800  804    37.50   110.70 
+angle       800  801  805    37.50   110.70 
+angle       800  801  806    37.50   110.70 
+angle       800  801  807    37.50   110.70 
+angle       802  800  803    33.00   107.80 
+angle       802  800  804    33.00   107.80 
+angle       805  801  806    33.00   107.80 
+angle       805  801  807    33.00   107.80 
+angle       806  801  807    33.00   107.80 
+angle       803  800  804    33.00   107.80 
+
+
+
+      ################################
+      ##                            ##
+      ##   Urey-Bradley Parameters  ##
+      ##                            ##
+      ################################
+
+
+ureybrad     35   34   35      38.25     1.5537
+
+
+      #####################################
+      ##                                 ##
+      ##  Improper Torsional Parameters  ##
+      ##                                 ##
+      #####################################
+
+
+
+imptors     801  802  800  803        0.000 180.0  2 
+imptors     801  802  800  804        0.000 180.0  2 
+imptors     800  805  801  806        0.000 180.0  2 
+imptors     800  805  801  807        0.000 180.0  2 
+
+
+
+      ############################
+      ##                        ##
+      ##  Torsional Parameters  ##
+      ##                        ##
+      ############################
+
+
+torsion     805  801  800  802        0.000  0.0  1  0.000 180.0  2  0.150  0.0  3 
+torsion     806  801  800  804        0.000  0.0  1  0.000 180.0  2  0.150  0.0  3 
+torsion     805  801  800  804        0.000  0.0  1  0.000 180.0  2  0.150  0.0  3 
+torsion     807  801  800  803        0.000  0.0  1  0.000 180.0  2  0.150  0.0  3 
+torsion     807  801  800  802        0.000  0.0  1  0.000 180.0  2  0.150  0.0  3 
+torsion     807  801  800  804        0.000  0.0  1  0.000 180.0  2  0.150  0.0  3 
+torsion     806  801  800  802        0.000  0.0  1  0.000 180.0  2  0.150  0.0  3 
+torsion     805  801  800  803        0.000  0.0  1  0.000 180.0  2  0.150  0.0  3 
+torsion     806  801  800  803        0.000  0.0  1  0.000 180.0  2  0.150  0.0  3 
+
+torsion       0    0    0    0        0.000  0.0  1  0.000 180.0  2  0.000  0.0  3
+
+      ########################################
+      ##                                    ##
+      ##  Atomic Partial Charge Parameters  ##
+      ##                                    ##
+      ########################################
+
+
+charge         800          -0.2396 
+charge         801          -0.2396 
+charge         802           0.0799 
+charge         803           0.0799 
+charge         804           0.0799 
+charge         805           0.0799 
+charge         806           0.0799 
+charge         807           0.0799 
+
diff --git a/docs/user_guide/ligpargen/images/ligpargen.png b/docs/user_guide/ligpargen/images/ligpargen.png
diff --git a/docs/user_guide/ligpargen/index.rst b/docs/user_guide/ligpargen/index.rst
@@ -0,0 +1,19 @@
+.. _LigParGen:
+
+********************
+The LigParGen Server
+********************
+
+The `LigParGen Server` is a website maintained by the Jorgensen group at Yale University
+for generating OPLS-AA parameters for small molecules, with up to 200 atoms. These
+parameters are compatible with the published OPLS-AA forcefield and are useful for
+molecules that are not covered by OPLS-AA out of the box.
+
+.. toctree::
+   :maxdepth: 2
+   :titlesonly:
+
+   overview
+   utility
+
+.. _LigParGen Server: https://zarbi.chem.yale.edu/ligpargen/index.html
diff --git a/docs/user_guide/ligpargen/overview.rst b/docs/user_guide/ligpargen/overview.rst
@@ -0,0 +1,20 @@
+********
+Overview
+********
+The heart of the main webpage is this section:
+
+.. figure:: images/ligpargen.png
+   :align: center
+   :alt: The LigParGen webpage
+
+   The LigParGen Webpage
+
+At the top you enter your structure as either SMILES or a Mol or PDB file. After
+choosing the charge model and specifying the charge, if any, click on **Submit
+Molecule**. For example, if you enter **CC** for ethane as the SMILES, select the **1.14*
+CMA1-LBCC** charge model and submit it, you will be sent to a page where you can
+download the parameters in formats suitable for various simulation codes. The **KEY**
+file for Tinker is easy to read if  bit verbose:
+
+.. include:: images/ethane.key
+   :literal:
diff --git a/docs/user_guide/ligpargen/utility.rst b/docs/user_guide/ligpargen/utility.rst
@@ -0,0 +1,3 @@
+**********************
+LigParGen utility tool
+**********************
diff --git a/forcefield_step/forcefield.py b/forcefield_step/forcefield.py
@@ -324,19 +324,19 @@ def setup_forcefield(self, P=None):
         else:
             # Find the forcefield file
             ff_file = P["forcefield_file"]
+            printer.important(
+                self.indent + 4 * " " + f"Reading the forcefield file {ff_file}"
+            )
             if ff_file.startswith("local:"):
-                ff_file = ff_file[6:]
+                ff_file = "Forcefields/" + ff_file[6:]
                 path = self.find_data_file(ff_file)
                 ff_file = str(path)
             else:
                 path = pkg_resources.resource_filename(__name__, "data/")
                 ff_file = os.path.join(path, P["forcefield_file"])
 
-            printer.important(self.indent + 4 * " " + "Reading the forcefield file")
             printer.important("")
-            printer.important(
-                self.indent + 8 * " " + P["forcefield_file"] + ": " + str(ff_file)
-            )
+            printer.important(self.indent + 8 * " " + str(ff_file))
             printer.important("")
 
             ff = seamm_ff_util.Forcefield(

diff --git a/forcefield_step/forcefield_parameters.py b/forcefield_step/forcefield_parameters.py
@@ -30,13 +30,22 @@
             if ext == ".frc":
                 forcefields.append(entry.name)
 
-# and local forcefields in the ~/SEAMM/Parameters directory
-local = Path.home() / "SEAMM" / "Parameters"
+# and local forcefields in the ~/SEAMM/data/Forcefields
+local = Path.home() / "SEAMM" / "data" / "Forcefields"
 logger.debug("Looking for forcefields at " + str(local))
 if local.exists():
     for path in local.glob("**/*.frc"):
         logger.debug("   " + str(path))
         forcefields.append("local:" + str(path.relative_to(local)))
+
+# and local forcefields in the ~/.seamm.d/data/Forcefields
+local = Path.home() / ".seamm.d" / "data" / "Forcefields"
+logger.debug("Looking for forcefields at " + str(local))
+if local.exists():
+    for path in local.glob("**/*.frc"):
+        logger.debug("   " + str(path))
+        forcefields.append("local:" + str(path.relative_to(local)))
+
 forcefields = sorted(forcefields)