added the "-report_duplicates" which warns the user when multiple bon…

…ds, angles, dihedrals, or impropers were automatically generated for the same set of atoms. This is important for the upcomming 2023 version of the OPLSAA force-field examples, which contain many duplicate or ambiguous bonded interaction definitions. This feature makes it possible for users to find out when other bonded interactions exist, and override the default choice manually.
jewettaij · Dec 2, 2024 · 361f9cc · 361f9cc
1 parent 9b29953
commit 361f9cc
Show file tree

Hide file tree

Showing 4 changed files with 167 additions and 36 deletions.
diff --git a/doc/moltemplate_manual_src/moltemplate_manual.tex b/doc/moltemplate_manual_src/moltemplate_manual.tex
@@ -1292,7 +1292,6 @@ \subsection{moltemplate.sh command line arguments:}
  See appendix \ref{sec:manual_assignment}.)
 \\
 \hline
-
 \begin{tabular}[t]{l}
 -overlay-bonds
 \\
@@ -1405,6 +1404,15 @@ \subsection{moltemplate.sh command line arguments:}
                 script. The commands are commented out and need manual editing
 		to get a meaningful simulation.
 \\
+\hline
+-report-duplicates filter\_id filter\_type &
+Generate a warning file to report whenever multiple bonds, angles, dihedrals,
+or impropers were found between the same atoms (regardless of whether ``-overlay-'' was used).
+The filter\_id and filter\_type arguments are strings which will ignore interactions
+whose id and type strings unless they contain filter\_id and filter\_type, respectively
+\textit{(If they are set to the empty string, "", then all duplicates are reported.
+When using the OPLSAA force-field, use ``-report-duplicates bytype \_\_''.)}
+\\
 \hline\hline
 \end{longtable}
 

diff --git a/moltemplate/remove_duplicates_nbody.py b/moltemplate/remove_duplicates_nbody.py
@@ -15,7 +15,9 @@
 
 """
 
+import os
 import sys
+from collections import defaultdict
 
 try:
     from .ttree_lex import SplitQuotedString
@@ -27,14 +29,27 @@
 def main():
     in_stream = sys.stdin
 
-    if len(sys.argv) == 2:
+    interaction_style = ""
+    log_warning_filename = ""
+    filter_id_str = ""
+    filter_type_str = ""
+    if len(sys.argv) > 1:
         n = int(sys.argv[1])
-    if (len(sys.argv) != 2) or (n < 1):
+    if len(sys.argv) > 3:
+        interaction_style = sys.argv[2]  # eg. "Dihedral"
+        log_warning_filename = sys.argv[3]  # eg. "warning_duplicate_dihedrals.txt"
+    if len(sys.argv) > 4:
+        filter_id_str = sys.argv[4]  # eg. "bytype"
+    if len(sys.argv) > 5:
+        filter_type_str = sys.argv[5]  # eg. "__"
+    if (len(sys.argv) not in (2,4,5,6)) or (n < 1):
         sys.stderr.write(
-            'Error (remove_duplicates_nbody.py): expected a positive integer argument.\n')
+            'Error (remove_duplicates_nbody.py): expected 1,3,4, or 5 arguments.\n')
         sys.exit(-1)
 
-    atom_ids_in_use = set([])
+    # Keep track of what interaction(s) were assigned to each set of atoms
+    atomids2interactions = defaultdict(list)  # Dict[Tuple[str], List[str])
+    atomids2interaction = {}  # of type Dict[Tuple[str], str]
 
     lines = in_stream.readlines()
 
@@ -54,18 +69,69 @@ def main():
                                    quotes='{',
                                    endquote='}')
 
-        if len(tokens) == 2 + n:
+        if len(tokens) == 0:
+            del lines[i]  # skip blank lines
+        elif len(tokens) == 2 + n:
             atom_ids = tuple(tokens[2:2 + n])
-            if atom_ids in atom_ids_in_use:
+            if atom_ids in atomids2interactions:
+                # If an interaction already exists between these atoms
+                # then delete this one (this line)
                 del lines[i]
             else:
-                atom_ids_in_use.add(atom_ids)
-        elif len(tokens) == 0:
-            del lines[i]
+                # Then this is the remaining interaction assigned to these atoms
+                atomids2interaction[atom_ids] = line
+            # Also keep track of all the interactions we did not keep.
+            atomids2interactions[atom_ids].append(line)
 
     for line in lines:
         sys.stdout.write(line)
 
+    # This portion of the code generates warning messages when duplicate
+    # interactions were deleted.  But the code is a little confusing because
+    # in order to be reported the interactions might need to satisfy a name
+    # filter.  Specifically the id_names ($) or type_names (@) might need
+    # to contain some string we are looking for.
+    if log_warning_filename != "":
+        with open(log_warning_filename, 'w') as log_warning:
+            # Did the user ask us to report duplicates?
+            for atom_ids, interactions in atomids2interactions.items():
+                if len(interactions) == 1:
+                    continue
+                chosen_interaction = atomids2interaction[atom_ids]
+                chosen_interaction_tokens = chosen_interaction.split()
+                interaction_id = chosen_interaction_tokens[0]  # eg. "$/dihedral:bytype6910"
+                interaction_type = chosen_interaction_tokens[1]  # eg. "@/dihedral:OPLSAA/HC_CM_CM_HC"
+                # Look for id_filter_string somewhere inside the interaction_id string.
+                if filter_id_str != "" and interaction_id.find(filter_id_str) == -1:
+                    # If that string is not found, then skip over this interaction.  This is used
+                    # to avoid complaining about interactions which lack a certain string.  (For
+                    # example, we only complain about interactions which were generated automatically.
+                    # Automatically generated interactions have id strings beginning with "bytype".)
+                    continue
+                # Look for type_filter_string somewhere inside the interaction_type string.
+                if filter_type_str != "" and interaction_type.find(filter_type_str) == -1:
+                    # If that string is not found, then skip over this interaction.  This is used
+                    # to avoid complaining about interactions which lack a certain string.  (For
+                    # example, in the OPLSAA force field, we only complain about "__" interactions.)
+                    continue
+                log_warning.write(
+                    f"Duplicate {interaction_style.lower()}s detected involving atoms:\n"
+                    "    "
+                    + "\n    ".join(atom_ids) + "\n"
+                    "  ...of types:\n"
+                    f"    "
+                    + "\n    ".join([int_str.split()[1] for int_str in interactions]) + "\n"
+                    f"  ...but only this {interaction_style.lower()} type was kept:\n"
+                    f"    {atomids2interaction[atom_ids].split()[1]}\n"
+                    f"  Was this the correct {interaction_style.lower()} type for these atoms?\n"
+                    f"  If not, create an explicit {interaction_style.lower()} interaction between those atoms\n"
+                    f"  in the \"Data {interaction_style}\" section of your molecule to override this choice.\n"
+                    "\n"
+                )
+            # If the warnings file is empty, there were no warnings.
+            # In that case, delete the file:
+            if os.stat(log_warning_filename).st_size == 0:
+                os.remove(log_warning_filename)
     return
 
 

diff --git a/moltemplate/scripts/moltemplate.sh b/moltemplate/scripts/moltemplate.sh
@@ -6,8 +6,8 @@
 # Copyright (c) 2013
 
 G_PROGRAM_NAME="moltemplate.sh"
-G_VERSION="2.20.22"
-G_DATE="2024-11-26"
+G_VERSION="2.20.23"
+G_DATE="2024-12-01"
 
 echo "${G_PROGRAM_NAME} v${G_VERSION} ${G_DATE}" >&2
 echo "" >&2
@@ -363,6 +363,12 @@ Optional arguments:
                     reordering of atom indexes, e.g. by allowing to define a bond
                     between atom 3 and 2 instead of 2 and 3.
 
+-report-duplicates filter_id filter_type
+              Warn when multiple bonds, angles, dihedrals, or impropers were found
+              between the same atoms (regardless of whether ``-overlay-'' was used).
+              The filter_id and filter_type arguments will ignore interactions whose
+              ids or types lack those strings. (Normally set them to "".  See manual.)
+
 -angle-symmetry file.py     Normally moltemplate.sh reorders the atoms in each
 -dihedral-symmetry file.py  angle, dihedral, improper, and bond interaction.
 -improper-symmetry file.py  TO TURN OFF ATOM REORDERING, SET file.py to "NONE"
@@ -438,6 +444,9 @@ REMOVE_DUPLICATE_BONDS="true"
 REMOVE_DUPLICATE_ANGLES="true"
 REMOVE_DUPLICATE_DIHEDRALS="true"
 REMOVE_DUPLICATE_IMPROPERS="true"
+REPORT_DUPLICATES=""
+RD_ID_FILTER=""
+RD_TYPE_FILTER=""
 SETTINGS_MOLC=""
 CHECKFF=""
 RUN_VMD_AT_END=""
@@ -500,6 +509,18 @@ while [ "$i" -lt "$ARGC" ]; do
         unset REMOVE_DUPLICATE_ANGLES
         unset REMOVE_DUPLICATE_DIHEDRALS
         unset REMOVE_DUPLICATE_IMPROPERS
+    elif [ "$A" = "-report-duplicates" ]; then
+        REPORT_DUPLICATES="true"
+        if [ "$i+1" -eq "$ARGC" ]; then
+            echo "ERROR: Expected 2 string arguments following the -report-duplicates argument" >&2
+            exit 7
+        fi
+        i=$((i+1))
+        eval A=\${ARGV${i}}
+        RD_ID_FILTER="$A"
+        i=$((i+1))
+        eval A=\${ARGV${i}}
+        RD_TYPE_FILTER="$A"
     elif [ "$A" = "-vmd" ]; then
         RUN_VMD_AT_END="true"
     elif [ "$A" = "-molc" ]; then
@@ -674,7 +695,7 @@ while [ "$i" -lt "$ARGC" ]; do
             # the transformation formula from the LAMMPS documentation 
             # https://lammps.sandia.gov/doc/Howto_triclinic.html  (which matches
             # http://www.ccl.net/cca/documents/molecular-modeling/node4.html)
-            TRICLINIC="True"
+            TRICLINIC="true"
             PI=3.1415926535897931
             BOXSIZE_X=$BOXSIZE_A
             BOXSIZE_Y=`awk -v PI="$PI" -v BOXSIZE_B="$BOXSIZE_B" -v GAMMA="$GAMMA" 'BEGIN{print BOXSIZE_B*sin(GAMMA*PI/180.0)}'`
@@ -767,7 +788,7 @@ while [ "$i" -lt "$ARGC" ]; do
           BOXSIZE_XY=${box[2]}
           BOXSIZE_XZ=${box[5]}
           BOXSIZE_YZ=${box[8]}
-          TRICLINIC="True"
+          TRICLINIC="true"
         fi
 
         # Save the coordinates.
@@ -1588,11 +1609,19 @@ if [ -s "${data_bonds}" ]; then
             ERR_INTERNAL
         fi
         mv "${data_bonds}.tmp" "${data_bonds}"
-        if ! $PYTHON_COMMAND "${PY_SCR_DIR}/remove_duplicates_nbody.py" 2 \
-                             < "${data_bonds}.template" \
-                             > "${data_bonds}.tmp"; then
-            ERR_INTERNAL
-        fi
+	if [ ! -z "$REPORT_DUPLICATES" ]; then
+          if ! $PYTHON_COMMAND "${PY_SCR_DIR}/remove_duplicates_nbody.py" 2 Bond warning_duplicate_bonds.txt "$RD_ID_FILTER" "$RD_TYPE_FILTER" \
+                               < "${data_bonds}.template" \
+                               > "${data_bonds}.tmp"; then
+              ERR_INTERNAL
+          fi
+	else
+          if ! $PYTHON_COMMAND "${PY_SCR_DIR}/remove_duplicates_nbody.py" 2 \
+                               < "${data_bonds}.template" \
+                               > "${data_bonds}.tmp"; then
+              ERR_INTERNAL
+          fi
+	fi
         mv "${data_bonds}.tmp" "${data_bonds}.template"
     fi
     if ! $PYTHON_COMMAND "${PY_SCR_DIR}/renumber_DATA_first_column.py" \
@@ -1626,10 +1655,18 @@ if [ -s "${data_angles}" ]; then
             ERR_INTERNAL
         fi
         mv "${data_angles}.tmp" "${data_angles}"
-        if ! $PYTHON_COMMAND "${PY_SCR_DIR}/remove_duplicates_nbody.py" 3 \
-                             < "${data_angles}.template" \
-                             > "${data_angles}.tmp"; then
-            ERR_INTERNAL
+	if [ ! -z "$REPORT_DUPLICATES" ]; then
+          if ! $PYTHON_COMMAND "${PY_SCR_DIR}/remove_duplicates_nbody.py" 3 Angle warning_duplicate_angles.txt "$RD_ID_FILTER" "$RD_TYPE_FILTER" \
+                               < "${data_angles}.template" \
+                               > "${data_angles}.tmp"; then
+              ERR_INTERNAL
+          fi
+	else
+	  if ! $PYTHON_COMMAND "${PY_SCR_DIR}/remove_duplicates_nbody.py" 3 \
+                               < "${data_angles}.template" \
+                               > "${data_angles}.tmp"; then
+              ERR_INTERNAL
+          fi
         fi
         mv "${data_angles}.tmp" "${data_angles}".template
     fi
@@ -1667,6 +1704,7 @@ fi
 
 
 
+
 if [ -s "${data_dihedrals}" ]; then
     SUBGRAPH_SCRIPT="nbody_Dihedrals.py"
     if [ -n "$SUBGRAPH_SCRIPT_DIHEDRALS" ]; then
@@ -1687,11 +1725,19 @@ if [ -s "${data_dihedrals}" ]; then
             ERR_INTERNAL
         fi
         mv "${data_dihedrals}.tmp" "${data_dihedrals}"
-        if ! $PYTHON_COMMAND "${PY_SCR_DIR}/remove_duplicates_nbody.py" 4 \
-                             < "${data_dihedrals}.template" \
-                             > "${data_dihedrals}.tmp"; then
-            ERR_INTERNAL
-        fi
+	if [ ! -z "$REPORT_DUPLICATES" ]; then
+          if ! $PYTHON_COMMAND "${PY_SCR_DIR}/remove_duplicates_nbody.py" 4 Dihedral warning_duplicate_dihedrals.txt "$RD_ID_FILTER" "$RD_TYPE_FILTER" \
+                               < "${data_dihedrals}.template" \
+                               > "${data_dihedrals}.tmp"; then
+              ERR_INTERNAL
+          fi
+	else
+          if ! $PYTHON_COMMAND "${PY_SCR_DIR}/remove_duplicates_nbody.py" 4 \
+                               < "${data_dihedrals}.template" \
+                               > "${data_dihedrals}.tmp"; then
+              ERR_INTERNAL
+          fi
+	fi
         mv "${data_dihedrals}.tmp" "${data_dihedrals}.template"
     fi
     if ! $PYTHON_COMMAND "${PY_SCR_DIR}/renumber_DATA_first_column.py" \
@@ -1726,6 +1772,9 @@ EOF
 fi
 
 
+
+
+
 if [ -s "${data_impropers}" ]; then
     SUBGRAPH_SCRIPT="nbody_Impropers.py"
     if [ -n "$SUBGRAPH_SCRIPT_IMPROPERS" ]; then
@@ -1746,11 +1795,19 @@ if [ -s "${data_impropers}" ]; then
             ERR_INTERNAL
         fi
         mv "${data_impropers}.tmp" "${data_impropers}"
-        if ! $PYTHON_COMMAND "${PY_SCR_DIR}/remove_duplicates_nbody.py" 4 \
-                             < "${data_impropers}.template" \
-                             > "${data_impropers}.tmp"; then
-            ERR_INTERNAL
-        fi
+	if [ ! -z "$REPORT_DUPLICATES" ]; then
+          if ! $PYTHON_COMMAND "${PY_SCR_DIR}/remove_duplicates_nbody.py" 4 Improper warning_duplicate_impropers.txt "$RD_ID_FILTER" "$RD_TYPE_FILTER" \
+                               < "${data_impropers}.template" \
+                               > "${data_impropers}.tmp"; then
+              ERR_INTERNAL
+          fi
+	else
+          if ! $PYTHON_COMMAND "${PY_SCR_DIR}/remove_duplicates_nbody.py" 4 \
+                               < "${data_impropers}.template" \
+                               > "${data_impropers}.tmp"; then
+              ERR_INTERNAL
+          fi
+	fi
         mv "${data_impropers}.tmp" "${data_impropers}.template"
     fi
     if ! $PYTHON_COMMAND "${PY_SCR_DIR}/renumber_DATA_first_column.py" \
@@ -1983,7 +2040,7 @@ if [ -s "$data_boundary" ]; then
     if [ -n "$BOXSIZE_XY" ] || [ -n "$BOXSIZE_XZ" ] || [ -n "$BOXSIZE_YZ" ]; then
         if [ -n "$BOXSIZE_XY" ] && [ -n "$BOXSIZE_XZ" ] && [ -n "$BOXSIZE_YZ" ]; then
             #echo "triclinic_parameters: XY XZ YZ = $BOXSIZE_XY $BOXSIZE_XZ $BOXSIZE_YZ" >&2
-            TRICLINIC="True"
+            TRICLINIC="true"
         else
             echo "Error: Problem with triclinic format (\"xy xz yz\") in \"$data_boundary\"" >&2
             exit 13

diff --git a/setup.py b/setup.py
@@ -45,9 +45,9 @@
 
   url='https://github.com/jewettaij/moltemplate',
 
-  download_url='https://github.com/jewettaij/moltemplate/archive/v2.20.22.zip',
+  download_url='https://github.com/jewettaij/moltemplate/archive/v2.20.23.zip',
 
-  version='2.20.22',
+  version='2.20.23',
 
   keywords=['simulation', 'LAMMPS', 'molecule editor', 'molecule builder',
             'ESPResSo'],