Modified iread generator of trajectory_reader to supply lattice vectors.

parser/parser-amber/AMBERParser.py

@@ -113,10 +113,11 @@ class AMBERParser(AmberC.AMBERParserBase):
         self.trajectoryFormat = None
         self.trajectoryFile = None
         self.readChunk = 300
-        self.unitcell = None
+        self.boxlengths = None
+        self.latticevectors = None
         self.atompositions = None
         # start with -1 since zeroth iteration is the initialization
-        self.mdIterNr = -1
+        self.MDiter = -1
         self.singleConfCalcs = []
         self.minConverged = None
         self.parsedLogFile = False
@@ -190,7 +191,10 @@ class AMBERParser(AmberC.AMBERParserBase):
                     self.atompositions = self.trajectory.iread()
                     if self.atompositions is not None:
                         self.topology = self.trajectory.get_topology()
+                        self.MDiter += 1
                         return fileFormat
+            else:
+                self.MDiter += 1
 
     def initializeFileHandlers(self):
         # Files will be loaded using their extensions initially.
@@ -531,9 +535,12 @@ class AMBERParser(AmberC.AMBERParserBase):
 
         if self.atompositions is not None:
             self.trajRefSingleConfigurationCalculation = gIndex
-            SloppyBackend.addArrayValues('atom_positions', np.transpose(
-                np.asarray(self.metaStorage.convertUnits(self.atompositions[0], "Angstrom", self.unitDict)
-                )))
+            unit_cell = np.asarray(self.metaStorage.convertUnits(
+                self.atompositions.unitcell_vectors[0], "Angstrom", self.unitDict))
+            SloppyBackend.addArrayValues('simulation_cell', unit_cell)
+            SloppyBackend.addArrayValues('lattice_vectors', unit_cell)
+            SloppyBackend.addArrayValues('atom_positions', np.transpose(np.asarray(
+                self.metaStorage.convertUnits(self.atompositions.xyz[0], "Angstrom", self.unitDict))))
             if self.topology is not None:
                 atom_labels = self.topologyDict["element"]
                 SloppyBackend.addArrayValues('atom_labels', np.asarray(atom_labels))
@@ -542,37 +549,15 @@ class AMBERParser(AmberC.AMBERParserBase):
             # Read the next step at trajectory in advance
             # If iread returns None, it will be the last step
             self.atompositions = self.trajectory.iread()
+            self.MDiter += 1
 
-#           if atom_vel:
+            #if atom_vel:
             # need to transpose array since its shape is [number_of_atoms,3] in the metadata
             #backend.addArrayValues('atom_velocities', np.transpose(np.asarray(atom_vel)))
 
         if (gIndex is None or gIndex == -1 or gIndex == "-1"):
             SloppyBackend.closeSection("section_system", self.secSystemGIndex)
 
-#        # For MD, the coordinates of the unit cell are not repeated.
-#        # Therefore, we have to store the unit cell, which was read in the beginning, i.e. scfIterNr == -1.
-#        if not self.MD or self.scfIterNr == -1:
-#            # write/store unit cell if present and set flag self.periodicCalc
-#            unit_cell = []
-#            for i in ['x', 'y', 'z']:
-#                uci = section['x_amber_geometry_lattice_vector_' + i]
-#                if uci is not None:
-#                    unit_cell.append(uci)
-#            if unit_cell:
-#                unit_cell = np.transpose(unit_cell)
-#                # from metadata: "The first index is x,y,z and the second index the lattice vector."
-#                # => unit_cell has already the right format
-#                if self.MD:
-#                    self.MDUnitCell = unit_cell
-#                else:
-#                    backend.addArrayValues('simulation_cell', unit_cell)
-#                self.periodicCalc = True
-        # write stored unit cell in case of MD
-#        if self.periodicCalc and self.MD and self.scfIterNr > -1:
-#            backend.addArrayValues('simulation_cell', self.MDUnitCell)
-#            backend.addArrayValues('configuration_periodic_dimensions', np.asarray([True, True, True]))
-
     def onOpen_section_single_configuration_calculation(self, backend, gIndex, section):
         # write the references to section_method and section_system
         backend.addValue('single_configuration_to_calculation_method_ref', self.secMethodGIndex)

parser/parser-amber/trajectory_reader.py

@@ -77,16 +77,29 @@ class TrajectoryReader(object):
                 logger.warning("ASE could not read the file '{}' with format '{}'. The contents might be malformed or wrong format used.".format(filename, file_format))
         return self.trajhandler
     
+    def load_topology(self):
+        """Loads the file handles for topology only
+        """
+        
+        if self.topohandler is None:
+            if self.topofile:
+                self.check_topology_format_support()
+        if self.topohandler is None:
+            logger.error("The topology file format '{}' is not supported by TrajectoryReader.".format(self.topoformat))
+
+        return self.topohandler
+    
     def get_file_format(self, filename, givenformat):
         """Returns extension of file
         """
         file_format = None
         if givenformat:
+            fname = filename.split("/")[-1]
             # MDTraj expects that extensions start with dot
-            if '.' not in filename:
-                file_format = "." + filename.lower()
+            if '.' not in fname:
+                file_format = "." + fname.lower()
             else:
-                file_format = filename.lower()
+                file_format = fname.lower()
         else:
             extension = filename.split(".")[-1]
             file_format = "." + extension.lower()
@@ -185,7 +198,7 @@ class TrajectoryReader(object):
         try:
             while True:
                 self.trajiter = next(iterator_object)
-                return self.trajiter.xyz
+                return self.trajiter
         except StopIteration:
             pass
         finally:
@@ -197,6 +210,12 @@ class TrajectoryReader(object):
         """
         return self.topohandler
 
+    def get_unitcell(self):
+        """Returns an iterator that goes through the given trajectory file one
+        configuration at a time.
+        """
+        return self.trajhandler.cell_lengths
+
     def load_mdtraj_topology(self, file_format, **kwargs):
         """Borrowed from mdtraj.trajectory to remove 
            IOError and TypeError warnings
@@ -382,7 +401,7 @@ class TrajectoryReader(object):
                                 for pos in positions:
                                     yield pos
                 except IOError:
-                    logger.warning("MDTraj could not read the file '{}' with format '{}'. The contents might be malformed or wrong format used.".format(self.trajfile, self.trajformat))
+                    #logger.warning("MDTraj could not read the file '{}' with format '{}'. The contents might be malformed or wrong format used.".format(self.trajfile, self.trajformat))
                     return
 
     def custom_iread(self, file_handle, format):
@@ -420,14 +439,128 @@ if __name__ == "__main__":
     if atom is not None:
         readaccess = True
     
-    if readaccess:
-        while True:
-            if atom is not None:
-                print(atom)
-                print(readdata.natoms())
-                i += 1
-            else:
-                break
-            atom = readdata.iread()
+#    if readaccess:
+#        while True:
+#            if atom is not None:
+#                print(atom)
+#                print(readdata.natoms())
+#                i += 1
+#            else:
+#                break
+#            atom = readdata.iread()
     print(i)
+#    angles = mdtraj.compute_angles(mytopology)
+#    dihedrals = mdtraj.compute_dihedrals(mytopology)
+#    print(dihedrals)
+    #print(' Chains: %s' % [res for res in mytopology.chains])
+    #print(' Residues: %s' % [res for res in mytopology.residues])
+    #print(' Atoms: %s' % [res for res in mytopology.atoms])
+    #print(' Bonds: %s' % [res for res in mytopology.bonds])
+    #print(mytopology.chain(0))
+    table, bonds = mytopology.to_dataframe()
+    print(bonds)
+    print(len(bonds))
+    array = table.get("resSeq")
+    #print([res for res in array])
+    mydict = table.to_dict(orient='list')
+    mylist = mydict["resSeq"]
+    print('mylist')
+    print(mylist)
+    print('mylist')
+    atomindex = np.arange(len(mylist))
+    passarray = np.zeros((len(mylist), 2), dtype=int)
+    passarray[:,0] = atomindex
+    passarray[:,1] = mylist
+    print(table)
+    print(passarray)
+    #print(bonds)
+#    for bond in bonds:
+#        print(bond)
+    mol_to_mol_bond = []
+    mol_to_mol_bond_num = []
+    atom_in_mol_bond = []
+    atom_in_mol_bond_num = []
+    index=0
+    print(len([res for res in mytopology.bonds]))
+    for res in mytopology.bonds:
+        molname1, molatom1 = str(res[0]).split('-')
+        molname2, molatom2 = str(res[1]).split('-')
+        if molname1 in molname2: 
+            atom_in_mol_bond.append(res)
+            atom_in_mol_bond_num.append(bonds[index])
+        else:
+            mol_to_mol_bond.append(res)
+            mol_to_mol_bond_num.append(bonds[index])
+        index += 1
+    print(mol_to_mol_bond)
+    print(np.array(mol_to_mol_bond_num))
+
+    atom_list = list(set(mydict["name"]))
+    atom_type_dict = {}
+    print(atom_list)
+    print(len(atom_list))
+    masses = {}
+    for ielm in range(len(atom_list)):
+        elm = atom_list[ielm]
+        atom_type_dict.update({elm : ielm})
+        for atom in mytopology.atoms:
+            if elm == atom.name:
+                masses.update({atom.name : atom.element})
+    print([atom_list, masses])
+    interDict = {}
+    interTypeDict = {}
+    for elm in atom_list:
+        bondList = []
+        typeList = []
+        bondid = 0
+        for bond in mytopology.bonds:
+            molname1, molatom1 = str(bond[0]).split('-')
+            molname2, molatom2 = str(bond[1]).split('-')
+            if (elm == str(molatom1) or elm == str(molatom2)):
+                bondList.append(list(bonds[bondid]))
+                #typeList.append(list([molatom1,molatom2]))
+                typeList.append(list([
+                    atom_type_dict[str(molatom1)],
+                    atom_type_dict[str(molatom2)]
+                    ]))
+                interDict.update({elm : bondList})
+                interTypeDict.update({elm : typeList})
+            bondid += 1
+
+#    print(interDict)
+    print(interDict["CH3"])
+    print(interTypeDict["CH3"])
+    print(len(interDict["CH3"][0]))
+
+#    print([[atom.element.radius, atom.index, atom.name, atom.residue.name, atom.element] for atom in mytopology.atoms])
+
+# Atom labels:
+    print(table)
+#    for atom in mytopology.atoms:
+#        print(atom)
+    labellist = mydict["element"]
+    #print(labellist)
+#    atomlabels = np.arange(len(labellist))
+#    labelarray = np.zeros((len(labellist), 2), dtype=int)
+#    labelarray[:,0] = atomlabels
+#    labelarray[:,1] = labellist
+#    print(labelarray)
+#    print(readdata.unitcell_lengths)
+    latticevectors = readdata.iread().unitcell_vectors
+    box = readdata.iread().unitcell_lengths
+    step = 0
+    for latvec in latticevectors:
+        print(step, latvec)
+        step += 1
+    for pbc in box:
+        print(step, box)
+        step += 1
+
+
+
+            
+
+
+
+