Adding Custom Torsion

This adds the custom torsion force reading the expression, parameters, and atoms from the XML from SMOG2

Author: contessoto

.gitignore

@@ -21,3 +21,4 @@ OpenSMOG/Test_Custom_Dihedral/testODdih/testODdih.gro
 OpenSMOG/Test_Custom_Dihedral/testODdih/testODdih.ndx
 OpenSMOG/Test_Custom_Dihedral/testODdih/testODdih.top
 OpenSMOG/Test_Custom_Dihedral/testODdih/testODdih.xml
+*.ipynb

OpenSMOG/OpenSMOG.py

@@ -463,6 +463,15 @@ def _splitForces_nb(self):
             forces[nb_data[0][n]] = [nb_data[1][n], nb_data[2][n], nb_data[3][n]]
         self.nonbond = forces
 
+    def _splitForces_dihedrals(self):
+        #Custom Dihedrals
+        dihedrals_data=self.data['dihedrals']
+        n_forces =  len(dihedrals_data[3])
+        forces = {}
+        for n in range(n_forces):
+            forces[dihedrals_data[3][n]] = [dihedrals_data[0][n], dihedrals_data[1][n], dihedrals_data[2][n]]
+        self.dihedrals = forces
+
     def _customSmogForce(self, name, data):
         #first set the equation
         contacts_ff = CustomBondForce(data[0])
@@ -488,6 +497,31 @@ def _customSmogForce(self, name, data):
         contacts_ff.setForceGroup(self.forceCount)
         self.forceCount +=1
 
+    def _customSmogForce_cd(self, name, data):
+        #first set the equation
+        dihedrals_ff = CustomTorsionForce(data[0])
+
+        #second set the number of variable
+        for pars in data[1]:
+            dihedrals_ff.addPerTorsionParameter(pars)
+
+        #third, apply the bonds from each pair of atoms and the related variables.
+        pars = [pars for pars in data[1]]
+
+        for iteraction in data[2]:
+            atom_index_i = int(iteraction['i'])-1 
+            atom_index_j = int(iteraction['j'])-1
+            atom_index_k = int(iteraction['k'])-1 
+            atom_index_l = int(iteraction['l'])-1
+            parameters = [float(iteraction[k]) for k in pars]
+
+            dihedrals_ff.addTorsion(atom_index_i, atom_index_j, atom_index_k, atom_index_l, parameters)
+
+        self.forcesDict[name] =  dihedrals_ff
+        dihedrals_ff.setForceGroup(self.forceCount)
+        self.forceCount +=1
+
+
     def _customSmogForce_nb(self, name, data):
         #first set the equation
         nonbond_ff = CustomNonbondedForce(data[0])
@@ -710,6 +744,41 @@ def import_xml2OpenSMOG(file_xml):
                 xml_data['nonbond']=[NonBond_Num,NBExpression,NBExpressionParameters,NBParameters]
             else:
                 self.nonbond_present=False
+
+
+            ## Custom Dihedrals 
+            CDForce_Names=[]
+            CDExpression=[]
+            CDParameters=[]
+            ijkl=[]
+            
+            self.dihedrals_present=False
+            if root.find('dihedrals') == None:
+                print('''
+        Only dihedrals contribution in the top file will be added
+        ''')
+            else:
+                self.dihedrals_present=True
+                dihedrals_xml=root.find('dihedrals')
+                for i in range(len(dihedrals_xml)):
+                    for name in dihedrals_xml[i].iter('dihedrals_type'):
+                        CDForce_Names.append(name.attrib['name'])
+
+                    for expr in dihedrals_xml[i].iter('expression'):
+                        CDExpression.append(expr.attrib['expr'])
+                        
+                    CDinternal_Param=[]
+                    for par in dihedrals_xml[i].iter('parameter'):
+                        CDinternal_Param.append(par.text)
+                    CDParameters.append(CDinternal_Param)
+
+                    internal_ijkl=[]
+                    for atoms_ijkl in dihedrals_xml[i].iter('interaction'):
+                            internal_ijkl.append(atoms_ijkl.attrib)
+                    ijkl.append(internal_ijkl)
+
+                xml_data['dihedrals']=[CDExpression,CDParameters,ijkl,CDForce_Names]
+
             return xml_data
 
         if not (self.forceApplied):
@@ -725,6 +794,13 @@ def import_xml2OpenSMOG(file_xml):
                     print("Creating Contacts force {:} from xml file".format(force))
                     self._customSmogForce(force, self.contacts[force])
                     self.system.addForce(self.forcesDict[force])
+
+            if self.dihedrals_present==True: 
+                self._splitForces_dihedrals()
+                for force in self.dihedrals:
+                    print("Creating Dihedrals force {:} from xml file".format(force))
+                    self._customSmogForce_cd(force, self.dihedrals[force])
+                    self.system.addForce(self.forcesDict[force])
 
             if self.nonbond_present==True: 
                 self._splitForces_nb()

OpenSMOG/share/OpenSMOG.xsd

@@ -1,11 +1,19 @@
 <?xml version="1.0" encoding="ISO-8859-1"?>
 <xs:schema xmlns:xs="http://www.w3.org/2001/XMLSchema">
   <!-- INTERACTION TYPE -->
-  <xs:complexType name="interactionType">
+  <xs:complexType name="contactInteractionType">
     <xs:attribute name="i" type="xs:integer" use="required"/>
     <xs:attribute name="j" type="xs:integer" use="required"/>
     <xs:anyAttribute processContents="skip"/>
   </xs:complexType>
+  <xs:complexType name="dihedralInteractionType">
+    <xs:attribute name="i" type="xs:integer" use="required"/>
+    <xs:attribute name="j" type="xs:integer" use="required"/>
+    <xs:attribute name="k" type="xs:integer" use="required"/>
+    <xs:attribute name="l" type="xs:integer" use="required"/>
+    <xs:anyAttribute processContents="skip"/>
+  </xs:complexType>
+ 
   <!-- nonbond_param TYPE -->
   <xs:complexType name="nonbondparamType">
     <xs:attribute name="type1" type="xs:string" use="required"/>
@@ -34,10 +42,20 @@
     <xs:sequence>
       <xs:element name="expression" type="expressionType" minOccurs="1" maxOccurs="unbounded"/>
       <xs:element name="parameter" type="parameterType" minOccurs="1" maxOccurs="unbounded"/>
-      <xs:element name="interaction" type="interactionType" minOccurs="1" maxOccurs="unbounded"/>
+      <xs:element name="interaction" type="contactInteractionType" minOccurs="1" maxOccurs="unbounded"/>
     </xs:sequence>
     <xs:attribute name="name" type="xs:string" use="required"/>
   </xs:complexType>
+  <!-- DIHEDRAL POTENTIAL TYPE -->
+  <xs:complexType name="dihedralpotType">
+    <xs:sequence>
+      <xs:element name="expression" type="expressionType" minOccurs="1" maxOccurs="unbounded"/>
+      <xs:element name="parameter" type="parameterType" minOccurs="1" maxOccurs="unbounded"/>
+      <xs:element name="interaction" type="dihedralInteractionType" minOccurs="1" maxOccurs="unbounded"/>
+    </xs:sequence>
+    <xs:attribute name="name" type="xs:string" use="required"/>
+  </xs:complexType>
+ 
   <!-- CONTACTS TYPE -->
   <xs:complexType name="contactsType">
     <xs:sequence>
@@ -55,6 +73,12 @@
       <xs:element name="constant" type="constantType" maxOccurs="unbounded"/>
     </xs:sequence>
   </xs:complexType>
+  <!-- DIHEDRALS TYPE -->
+  <xs:complexType name="dihedralsType">
+    <xs:sequence>
+      <xs:element name="dihedrals_type" type="dihedralpotType" maxOccurs="unbounded"/>
+    </xs:sequence>
+  </xs:complexType>
   <!-- NB TYPE -->
   <xs:complexType name="nbType">
     <xs:sequence>
@@ -74,6 +98,7 @@
   <xs:complexType name="OpenSMOGType">
     <xs:all>
       <xs:element name="contacts" type="contactsType" minOccurs="0" maxOccurs="1"/>
+      <xs:element name="dihedrals" type="dihedralsType" minOccurs="0" maxOccurs="1"/>
       <xs:element name="constants" type="constantsType" minOccurs="0" maxOccurs="1"/>
       <xs:element name="nonbond" type="nonbondType" minOccurs="0" maxOccurs="1"/>
     </xs:all>