From 4d4b105380d785d9755efc0ec49dc5a450bd67eb Mon Sep 17 00:00:00 2001 From: Andrew Jewett Date: Sat, 1 Feb 2025 21:50:21 -0500 Subject: [PATCH] Corrected a bug in OPLSAA ("oplsaa.lt") which incorrectly named several angle and dihedral types (@angle:NB_OA_CR, @angle:NB_OA_CR, @dihedral:H~_N~_CT_CT). Those angle and dihedral type names have been restored. Also: the "oplsaa2lt.py" conversion script has been debugged and cleaned up. Most users will not be affecgted by this change. --- .../convert_OPLSAA_to_LT/oplsaa2lt.py | 173 +++++++++++------- moltemplate/force_fields/oplsaa.lt | 18 +- moltemplate/scripts/moltemplate.sh | 4 +- setup.py | 4 +- 4 files changed, 119 insertions(+), 80 deletions(-) diff --git a/moltemplate/force_fields/convert_OPLSAA_to_LT/oplsaa2lt.py b/moltemplate/force_fields/convert_OPLSAA_to_LT/oplsaa2lt.py index c0a35ef..5ad3c0e 100755 --- a/moltemplate/force_fields/convert_OPLSAA_to_LT/oplsaa2lt.py +++ b/moltemplate/force_fields/convert_OPLSAA_to_LT/oplsaa2lt.py @@ -17,8 +17,8 @@ assert version.major > 3 or (version.major == 3 and version.minor >= 7) __author__ = "Domenico Marson and Andrew Jewett" -__version__ = '1.0.0' -__date__ = '2025-1-31' +__version__ = '1.0.1' +__date__ = '2025-2-01' g_program_name = __file__.split('/')[-1] @@ -168,84 +168,126 @@ def get_dihedrals_and_impropers(input_lines) -> tuple[list[Dihedral], list[Impro return loaded_dihedrals, loaded_impropers -# NOTE: PROBLEM WITH SAME-TYPES BONDED INTERACTIONS... -# The same dihedral (from atom types POV) can have != parameters, based on comment... -# e.g.: -# 102 0.000 5.500 0.00 0.0 O -C -OH-HO carboxylic acid - aliphatic -# 210 0.000 5.000 0.00 0.0 O -C -OH-HO benzoic acids -# the same problem is observed also for bonds and angles... -def check_uniqueness( +def delete_redundant_duplicates( interactions: list[Dihedral]|list[Improper]|list[Angle]|list[Bond], ) -> None: + """ OPLSAA files contain many duplicates. Decide which interactions to keep and return the list to the caller.""" + ##### My apologies for this ugly code. ##### - # It wasn't planned that way, but we keep finding - # more special cases that need custom code added. - # PAR and SB files have grown over time and we keep - # finding more redundancy problems that need fixing. - - # Step1: Remove duplicate interactions. - # Create a dictionary of all the duplicate interactions - # which share the same atom types and coefficients - # (excluding comments). - typescoeff_to_interaction = defaultdict(dict) + # It wasn't planned that way. The PAR and SB files + # that store OPLSAA parameters have grown messy over time. + # We keep finding more redundancy problems that need fixing. + # That's why this code has so many edge cases. + + # ----------- Step 1: --------------- + # Organize interactions according to: + # -atom type strings (a tuple which may include wildcards) + # -force-field parameters (This is called "params" in the code.) + # -the "paramstr" (This is the entire string following the atom type list. + # including the parameters AND comments, if present.) + # The result is stored in a dictionary-of-dictionary-of-dictionaries. + # named "types_to_paraminteractions". Lookup interactions this way: + # types_to_paraminteractions[atomtypes][params][paramstr] --> interaction + # + # Later, during Step2, we will use it to figure out which interactions + # are redundant and can be discarded. + types_to_paraminteractions = {} for interaction in interactions: types = tuple(interaction.types) - # Find the text containing the force field parameters ("coeffs") - coeff = interaction.coeff_line - # Typically, this is "angle_coeff @angle:C3_N~_H~ 38. 118.4 # " - # But we only want "38. 118.4 # ". + # Find the text containing the force field parameters ("paramstr") + paramstr = interaction.coeff_line + # Typically, this is "angle_coeff @angle:C3_N~_H~ 38. 118.4 # WJ94" + # But we only want "38. 118.4 # WJ94". # So we ignore the text before the first two spaces: - ispace1 = coeff.find(" ") - ispace2 = coeff.find(" ", ispace1+1) + ispace1 = paramstr.find(" ") + ispace2 = paramstr.find(" ", ispace1+1) ispace2 = ispace2+1 - coeff = coeff[ispace2:].strip() # eg. "38. 118.4 #" - # Finally, strip off the comment - coeff_no_comment = coeff.lstrip("#").split("#")[0].strip() - # Instead of storing the interaction at (types, coeff_no_comment), - # I store an additional dictionary (whose key is (types, coeff)). - # That way I can keep track of all the different comments that - # were added to this (otherwise identical) interaction. - # We want to keep all of them, but discard the ones that are - # trivially similar to each other. - typescoeff_to_interaction[(types, coeff_no_comment)][(types, coeff)] = interaction - # Copy the the non-redudant entries back into "interactions" again. - i = 0 - for coeff_no_comment, coeff_interactions in typescoeff_to_interaction.items(): - for typescoeff, interaction in coeff_interactions.items(): - # Edge case: We want to ignore interactions if - # they lack a comment and are otherwise identical. - # ...so we also check for (coeff != coeff_no_comment). - coeff = typescoeff[1] - coeff_blank_comment = coeff - i_comment = coeff.find("#") - if i_comment > 0: - coeff_blank_comment = coeff[:i_comment+1] - comment = coeff[i_comment+1:] - if len(coeff_interactions) == 1: - interactions[i] = interaction - i += 1 - elif coeff not in (coeff_no_comment, coeff_blank_comment): - # For some reason, a lot of comments only contain ' "' - # We want to skip those too. - if comment.strip() != '"': - interactions[i] = interaction - i += 1 - - del interactions[i:] + paramstr = paramstr[ispace2:].strip() # eg. "38. 118.4 # WJ94" + # Finally, strip off the comment, leaving only the parameters + params = paramstr.lstrip("#").split("#")[0].strip() # eg. "38. 118.4" + if types not in types_to_paraminteractions: + # Multiple interactions can exist for these same atom types + # differing by either the parameters, or comments, or both. + # Store a dictionary that looks up the interaction between these + # atoms according to their params (without comments). + types_to_paraminteractions[types] = {} # We will use this below + + # Create a dictionary that looks up all the interactions + # which share the same types and parameters (.ie params), but may + # have different comments. We might want to delete these duplicates + # later, but for now, we keep track of all of them. + params_to_interactions = types_to_paraminteractions[types] + if params not in params_to_interactions: + params_to_interactions[params] = {} + paramstr_to_interactions = params_to_interactions[params] + + # For these interactions, its possible that multiple interactions + # may exist in the file with identical atom types, parameters (coeffs) + # AND comments. To get rid of these trivial duplicates, we store them in + # a dictionary, indexed by the original paramstr (including comments). + paramstr_to_interactions[paramstr] = interaction + + # ----------- Step 2: --------------- + # Decide which of these interactions should be discarded. + # (The interactions we want to keep will be stored in + # out_interactions, which will be returned to the caller.) + del interactions[:] + for types, params_to_interactions in types_to_paraminteractions.items(): + for params, paramstr_to_interactions in params_to_interactions.items(): + for paramstr, interaction in paramstr_to_interactions.items(): + discarded = True + paramstr_blank_comment = paramstr + i_comment = paramstr.find("#") + if i_comment > 0: + paramstr_blank_comment = paramstr[:i_comment+1] + comment = paramstr[i_comment+1:] + if len(paramstr_to_interactions) == 1: + interactions.append(interaction) + discarded = False + elif ( + paramstr not in (params, paramstr_blank_comment) # case 1 + and (comment.strip() != '"') # case 2 (see below) + ): + # Edge Case 1: + # We want to ignore interactions if + # they lack a comment but they are otherwise identical. + # ...so we also check to make sure that paramstr != params. + # (and also paramstr != paramstr_blank_comment) + # + # Edge Case 2: + # For some reason, a lot of comments only contain '"'. We + # want to skip those too (since they are otherwise identical) + # + # If none of these edge-cases are true, + # then we don't discard the interaction. + interactions.append(interaction) + discarded = False + # If all of the duplicate interactions are identical except for + # the comments following the params, then throw away these + # duplicates (merge them into a single interaction). + if (len(params_to_interactions) == 1) and (not discarded): + break # skip the remaining duplicates for these atoms + +# NOTE: PROBLEM WITH SAME-TYPES BONDED INTERACTIONS... +# The same dihedral (from atom types POV) can have != parameters, based on comment... +# e.g.: +# 102 0.000 5.500 0.00 0.0 O -C -OH-HO carboxylic acid - aliphatic +# 210 0.000 5.000 0.00 0.0 O -C -OH-HO benzoic acids +# the same problem is observed also for bonds and angles... +def count_nonredundant_duplicates( + interactions: list[Dihedral]|list[Improper]|list[Angle]|list[Bond], +) -> None: for idx, it1 in enumerate(interactions): for it2 in interactions[idx+1:]: if it1.types == it2.types: - it1_coeff = it1.coeff_line.lstrip("#").split("#")[0] - it2_coeff = it2.coeff_line.lstrip("#").split("#")[0] - # If we are not skipping this interaction, then - # keep track of the number of duplicate interactions of this type if it1.duplicate_count == 0: it1.duplicate_count = 1 it2.duplicate_count = it1.duplicate_count + 1 + def sort_duplicates( interactions: list[Dihedral]|list[Improper]|list[Angle]|list[Bond], ) -> None: @@ -335,10 +377,13 @@ def main(argv): dihedrals, impropers = get_dihedrals_and_impropers(lines_par) + # Now, let's cleanup all the lists of bonded interactions. + # (Remove redundant entries, and sort by atom-type) for interactions in [bonds, angles, dihedrals, impropers]: interactions.sort(key=lambda x: x.typename) interactions.sort(key=lambda x: x.sort_key) - check_uniqueness(interactions) + delete_redundant_duplicates(interactions) + count_nonredundant_duplicates(interactions) sort_duplicates(interactions) diff --git a/moltemplate/force_fields/oplsaa.lt b/moltemplate/force_fields/oplsaa.lt index dcce222..3db5caa 100644 --- a/moltemplate/force_fields/oplsaa.lt +++ b/moltemplate/force_fields/oplsaa.lt @@ -5731,8 +5731,7 @@ OPLSAA { angle_coeff @angle:CT_CT_CZ 58.35 112.7 # wlj angle_coeff @angle:CT_CT_Cl 69. 109.8 # wlj - from MM2 angle_coeff @angle:CT_CT_F~ 50. 109.5 # PAK F-CT-HC (emd 5-09-94) - angle_coeff @angle:CT_CT_HC__2 37.5 110.7 # CHARMM 22 parameter file - angle_coeff @angle:CT_CT_HC__1 37.5 110.7 # CHARMM 22 + angle_coeff @angle:CT_CT_HC 37.5 110.7 # CHARMM 22 angle_coeff @angle:CT_CT_I~ 75. 112.0 # wlj angle_coeff @angle:CT_CT_N2 80.0 111.2 # ARG JCP 76, 1439 angle_coeff @angle:CT_CT_N3 80.0 111.2 # LYS(OL) JCP 76, 1439 @@ -6296,8 +6295,7 @@ OPLSAA { angle_coeff @angle:NB_NB_OA 70. 110.0 # JT-R oxatriazoles angle_coeff @angle:NB_NB_SA 70. 114.0 # JT-R thiatriazole angle_coeff @angle:NB_NX_CW 56. 113.1 # MKD synonym for NB-NA-CW - angle_coeff @angle:NB_OA_CR__2 70. 101.3 # JT-R oxatriazole - angle_coeff @angle:NB_OA_CR__1 70. 101.3 # JT-R oxatriazoles + angle_coeff @angle:NB_OA_CR 70. 101.3 # JT-R oxatriazoles angle_coeff @angle:NB_OA_NB 70. 103.4 # JT-R oxatriazole angle_coeff @angle:NB_OS_CR 70. 101.3 # JT-R oxatriazoles angle_coeff @angle:NB_OS_DM 10.0 125. # wlj @@ -7166,8 +7164,7 @@ OPLSAA { @angle:CT_CT_CZ @atom:*_b*_aCT*_d*_i* @atom:*_b*_aCT*_d*_i* @atom:*_b*_aCZ*_d*_i* @angle:CT_CT_Cl @atom:*_b*_aCT*_d*_i* @atom:*_b*_aCT*_d*_i* @atom:*_b*_aCl*_d*_i* @angle:CT_CT_F~ @atom:*_b*_aCT*_d*_i* @atom:*_b*_aCT*_d*_i* @atom:*_b*_aF~*_d*_i* - @angle:CT_CT_HC__2 @atom:*_b*_aCT*_d*_i* @atom:*_b*_aCT*_d*_i* @atom:*_b*_aHC*_d*_i* - @angle:CT_CT_HC__1 @atom:*_b*_aCT*_d*_i* @atom:*_b*_aCT*_d*_i* @atom:*_b*_aHC*_d*_i* + @angle:CT_CT_HC @atom:*_b*_aCT*_d*_i* @atom:*_b*_aCT*_d*_i* @atom:*_b*_aHC*_d*_i* @angle:CT_CT_I~ @atom:*_b*_aCT*_d*_i* @atom:*_b*_aCT*_d*_i* @atom:*_b*_aI~*_d*_i* @angle:CT_CT_N2 @atom:*_b*_aCT*_d*_i* @atom:*_b*_aCT*_d*_i* @atom:*_b*_aN2*_d*_i* @angle:CT_CT_N3 @atom:*_b*_aCT*_d*_i* @atom:*_b*_aCT*_d*_i* @atom:*_b*_aN3*_d*_i* @@ -7731,8 +7728,7 @@ OPLSAA { @angle:NB_NB_OA @atom:*_b*_aNB*_d*_i* @atom:*_b*_aNB*_d*_i* @atom:*_b*_aOA*_d*_i* @angle:NB_NB_SA @atom:*_b*_aNB*_d*_i* @atom:*_b*_aNB*_d*_i* @atom:*_b*_aSA*_d*_i* @angle:NB_NX_CW @atom:*_b*_aNB*_d*_i* @atom:*_b*_aNX*_d*_i* @atom:*_b*_aCW*_d*_i* - @angle:NB_OA_CR__2 @atom:*_b*_aNB*_d*_i* @atom:*_b*_aOA*_d*_i* @atom:*_b*_aCR*_d*_i* - @angle:NB_OA_CR__1 @atom:*_b*_aNB*_d*_i* @atom:*_b*_aOA*_d*_i* @atom:*_b*_aCR*_d*_i* + @angle:NB_OA_CR @atom:*_b*_aNB*_d*_i* @atom:*_b*_aOA*_d*_i* @atom:*_b*_aCR*_d*_i* @angle:NB_OA_NB @atom:*_b*_aNB*_d*_i* @atom:*_b*_aOA*_d*_i* @atom:*_b*_aNB*_d*_i* @angle:NB_OS_CR @atom:*_b*_aNB*_d*_i* @atom:*_b*_aOS*_d*_i* @atom:*_b*_aCR*_d*_i* @angle:NB_OS_DM @atom:*_b*_aNB*_d*_i* @atom:*_b*_aOS*_d*_i* @atom:*_b*_aDM*_d*_i* @@ -8915,8 +8911,7 @@ OPLSAA { dihedral_coeff @dihedral:H~_NT_NT_H~ 0.0 0.0 0.3 0.0 # generic dihedral_coeff @dihedral:H~_NT_OH_HO 0.0 0.0 0.3 0.0 # generic dihedral_coeff @dihedral:H~_NT_OS_CT 0.0 0.0 0.3 0.0 # generic - dihedral_coeff @dihedral:H~_N~_CT_CT__2 0.000 0.000 0.000 0.0 # N-ethylformamide - dihedral_coeff @dihedral:H~_N~_CT_CT__1 0.000 0.000 0.000 0.0 # peptides + dihedral_coeff @dihedral:H~_N~_CT_CT 0.000 0.000 0.000 0.0 # peptides dihedral_coeff @dihedral:H~_N~_CT_C~ 0.000 0.000 0.000 0.0 # peptides dihedral_coeff @dihedral:H~_N~_C~_C~ 0.000 4.900 0.000 0.0 # dicarbonyls " dihedral_coeff @dihedral:H~_N~_OH_HO 2.722 -5.154 0.000 0.0 # hydroxamic acids @@ -10114,8 +10109,7 @@ OPLSAA { @dihedral:H~_NT_NT_H~ @atom:*_b*_a*_dH~*_i* @atom:*_b*_a*_dNT*_i* @atom:*_b*_a*_dNT*_i* @atom:*_b*_a*_dH~*_i* @dihedral:H~_NT_OH_HO @atom:*_b*_a*_dH~*_i* @atom:*_b*_a*_dNT*_i* @atom:*_b*_a*_dOH*_i* @atom:*_b*_a*_dHO*_i* @dihedral:H~_NT_OS_CT @atom:*_b*_a*_dH~*_i* @atom:*_b*_a*_dNT*_i* @atom:*_b*_a*_dOS*_i* @atom:*_b*_a*_dCT*_i* - @dihedral:H~_N~_CT_CT__2 @atom:*_b*_a*_dH~*_i* @atom:*_b*_a*_dN~*_i* @atom:*_b*_a*_dCT*_i* @atom:*_b*_a*_dCT*_i* - @dihedral:H~_N~_CT_CT__1 @atom:*_b*_a*_dH~*_i* @atom:*_b*_a*_dN~*_i* @atom:*_b*_a*_dCT*_i* @atom:*_b*_a*_dCT*_i* + @dihedral:H~_N~_CT_CT @atom:*_b*_a*_dH~*_i* @atom:*_b*_a*_dN~*_i* @atom:*_b*_a*_dCT*_i* @atom:*_b*_a*_dCT*_i* @dihedral:H~_N~_CT_C~ @atom:*_b*_a*_dH~*_i* @atom:*_b*_a*_dN~*_i* @atom:*_b*_a*_dCT*_i* @atom:*_b*_a*_dC~*_i* @dihedral:H~_N~_C~_C~ @atom:*_b*_a*_dH~*_i* @atom:*_b*_a*_dN~*_i* @atom:*_b*_a*_dC~*_i* @atom:*_b*_a*_dC~*_i* @dihedral:H~_N~_OH_HO @atom:*_b*_a*_dH~*_i* @atom:*_b*_a*_dN~*_i* @atom:*_b*_a*_dOH*_i* @atom:*_b*_a*_dHO*_i* diff --git a/moltemplate/scripts/moltemplate.sh b/moltemplate/scripts/moltemplate.sh index 61883de..4f57718 100755 --- 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.22.1" -G_DATE="2025-1-31" +G_VERSION="2.22.2" +G_DATE="2025-2-01" echo "${G_PROGRAM_NAME} v${G_VERSION} ${G_DATE}" >&2 echo "" >&2 diff --git a/setup.py b/setup.py index 851674c..efe99fc 100644 --- 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.22.1.zip', + download_url='https://github.com/jewettaij/moltemplate/archive/v2.22.2.zip', - version='2.22.1', + version='2.22.2', keywords=['simulation', 'LAMMPS', 'molecule editor', 'molecule builder', 'ESPResSo'],