diff --git a/prody/dynamics/functions.py b/prody/dynamics/functions.py
index 92f69266a..fc867836e 100644
--- a/prody/dynamics/functions.py
+++ b/prody/dynamics/functions.py
@@ -396,22 +396,37 @@ def parseScipionModes(metadata_file, title=None, pdb=None):
     vectors = np.zeros((dof, n_modes))
     vectors[:, 0] = mode1
 
-    eigvals = np.zeros(n_modes)
+    eigvals = {}
+    indices = []
 
     try:
-        eigvals[0] = float(row1['_nmaEigenval'])
+        indices.append(int(row1['_order_'])-1)
+        found_indices = True
+    except KeyError:
+        found_indices = False
+
+    try:
+        if found_indices:
+            eigvals[indices[0]] = float(row1['_nmaEigenval'])
+        else:
+            eigvals[0] = float(row1['_nmaEigenval'])
         found_eigvals = True
-    except:
+    except KeyError:
         found_eigvals = False
 
     for i, row in enumerate(star_loop[1:]):
         vectors[:, i+1] = parseArray(top_dirs + row['_nmaModefile']).reshape(-1)
-        if found_eigvals:
+        if found_eigvals and not found_indices:
             eigvals[i+1] = float(row['_nmaEigenval'])
+        if found_indices:
+            index = int(row['_order_'])-1
+            indices.append(index)
+            if found_eigvals:
+                eigvals[index] = float(row['_nmaEigenval'])
     
     if not found_eigvals:
-        log_fname = run_path + '/logs/run.stdout'
-        fi = open(log_fname, 'r')
+        logFileName = run_path + '/logs/run.stdout'
+        fi = open(logFileName, 'r')
         lines = fi.readlines()
         fi.close()
 
@@ -435,7 +450,22 @@ def parseScipionModes(metadata_file, title=None, pdb=None):
     else:
         nma = GNM(title)
 
-    nma.setEigens(vectors, eigvals)
+    nma.setEigens(vectors, np.array(list(eigvals.values())))
+
+    if found_indices:
+        n_modes = np.max(indices) + 1
+        vectors2 = np.ones((dof, n_modes))
+        eigvals2 = np.ones((n_modes))
+
+        for i, index in enumerate(indices):
+            vectors2[:, index] = vectors[:, i]
+            eigvals2[index] = eigvals[index]
+
+        nma2 = type(nma)(title)
+        nma2.setEigens(vectors2, eigvals2)
+
+        nma = nma2[indices]
+
     return nma
 
 
@@ -495,7 +525,7 @@ def writeScipionModes(output_path, modes, write_star=False, scores=None,
         raise TypeError('collectivityThreshold should be float, not {0}'
                         .format(type(collectivityThreshold)))
 
-    if modes.numModes() == 1 and not isinstance(modes, NMA):
+    if modes.numModes() == 1 and not isinstance(modes, (NMA, ModeSet)):
         old_modes = modes
         modes = NMA(old_modes)
         modes.setEigens(old_modes.getArray().reshape(-1, 1))