Implemented LISP for a generic ring

SlipperySlope.py

@@ -2,6 +2,9 @@
 import sys
 import numpy as np
 
+__version__ = "0.1.1"
+
+
 class Molecule:
     """Molecule loading"""
     def __init__(self):
@@ -13,15 +16,16 @@ def __init__(self):
             self.xyz = []
             for line in lines[2:]:
                 line = line.split()
-                self.element.append(line[0])  
+                self.element.append(line[0])
                 self.xyz.append(
-                    np.array([float(line[i]) for i in range(1,4)])
-                    )
+                    np.array([float(line[i]) for i in range(1, 4)])
+                )
 
         except:
             print('Error during importing molecule. Check XMol format!')
             quit()
 
+
 def dist(atom_a, atom_b):
     """Calulate the distance between two atoms"""
     return np.linalg.norm(atom_a-atom_b)
@@ -34,55 +38,87 @@ def ang(atom_a, atom_b, atom_c):
     rad_angle = np.arccos(
         np.dot(ba, bc) / (dist(atom_a, atom_b) * dist(atom_b, atom_c))
     )
-    return np.degrees(rad_angle)
+    return rad_angle
+
+
+def lisp(ring_coor, M_coor):
+    """Calculate the Label Independednt Slippage Parameter"""
+    ring_ctd = np.mean(ring_coor, axis=0)
+    d_ring_ctd_M = dist(ring_ctd, M_coor)
+    d_list_norm = []
+    for idx, val in enumerate(ring_coor):
+        if idx != len(ring_coor)-1:
+            d_ring_pair_mid = np.mean(
+                [ring_coor[idx], ring_coor[idx+1]], axis=0)
+            print("<{}-{}>...\t\t{:.6f}  {:.6f}  {:.6f}".format(
+                ring_idx[idx] + 1,
+                ring_idx[idx+1] + 1,
+                d_ring_pair_mid[0],
+                d_ring_pair_mid[1],
+                d_ring_pair_mid[2]
+                ))
+        else:
+            d_ring_pair_mid = np.mean([ring_coor[idx], ring_coor[0]], axis=0)
+            print("<{}-{}>...\t\t{:.6f}  {:.6f}  {:.6f}".format(
+                ring_idx[idx] + 1,
+                ring_idx[0] + 1,
+                d_ring_pair_mid[0],
+                d_ring_pair_mid[1],
+                d_ring_pair_mid[2]
+                ))
 
+        theta = ang(d_ring_pair_mid, ring_ctd, M_coor)
+        d_list_norm.append(abs(d_ring_ctd_M*np.sin(theta-np.pi/2)))
+
+    return np.mean(d_list_norm), ring_ctd, d_ring_ctd_M
 
 print("""
     ================================
-                 LISP   
-    ================================
-""")
+             SlipperySlope
+    ====================== Mk. {}
+""".format(__version__))
 
-mol=Molecule()
+MOL = Molecule()
 print("Load molecule...\tDONE!\n")
 
 print("--> INPUT\n")
-ring_idx=[int(i)-1 for i in input("Index for the ring...\t").split()]
-M_idx=int(input("Index for the M atom...\t"))
+try:
+    ring_idx = [int(i)-1 for i in input("Index for the ring...\t").split()]
+    M_idx = int(input("Index for the M atom...\t"))-1
+except ValueError as val_err:
+    print("The input index is not valid: ", val_err)
+    quit()
 
-##sanity check
 print("\n--> SANITY CHECK\n")
 
 try:
-    print("Ring definition:\t\t", end="")
-    if len(ring_idx)>=3:
+    print("Ring definition...\t\t", end="")
+    if len(ring_idx) >= 3:
         print("OK!")
     else:
         print("ERROR!")
         raise ValueError()
-    d_list=[]
 
+    d_list = []
     for idx, val in enumerate(ring_idx):
-        if val != ring_idx[-1]:
+        if idx != len(ring_idx)-1:
             d_list.append(
-                dist(mol.xyz[val], mol.xyz[ring_idx[idx+1]])
-                )
+                dist(MOL.xyz[val], MOL.xyz[ring_idx[idx+1]])
+            )
         else:
             d_list.append(
-                dist(mol.xyz[val], mol.xyz[ring_idx[0]])
-                )
-    d_list_var=np.var(d_list)
-    d_list_std=np.std(d_list)
-    print("Variance:\t\t{:.3f}".format(d_list_var))
-
-    print("Std. deviation:\t\t", end="")
+                dist(MOL.xyz[val], MOL.xyz[ring_idx[0]])
+            )
+    d_list_var = np.var(d_list)
+    d_list_std = np.std(d_list)
+    print("Variance...\t\t{:.3f}".format(d_list_var))
+    print("Std. deviation...\t", end="")
     if d_list_std < 0.1:
-        print("{:.3f} => OK!".format(d_list_std))
+        print("{:.3f}\tOK!".format(d_list_std))
     else:
-        print("{:.3f} => TOO HIGH!".format(d_list_std))
+        print("{:.3f}\tTOO HIGH!".format(d_list_std))
         raise ValueError()
 
-
 except ValueError:
     print("\nWARNING: CALCULATION ABORTED DUE TO UNRELIABLE RESULTS")
     print("""
@@ -92,4 +128,16 @@ def ang(atom_a, atom_b, atom_c):
     """)
     quit()
 
-###
+print("\n--> LISP CALCULATION\n")
+ring_coor = [MOL.xyz[i] for i in ring_idx]
+LISP, ring_ctd, d_ring_ctd_M = lisp(ring_coor, MOL.xyz[M_idx])
+
+print("\nRing centroid...\t{:.6f}  {:.6f}  {:.6f}".format(ring_ctd[0], ring_ctd[1], ring_ctd[2]))
+print("Dist. centroid-M...\t{:.3f}".format(d_ring_ctd_M))
+print("LISP...\t\t\t{:.3f}".format(LISP))
+
+print("""
+    ================================
+           NORMAL TERMINATION
+    ================================
+""")

example.xyz

@@ -0,0 +1,27 @@
+25
+CoC11NH12
+C        0.750298000     -1.296668000     -0.481423000
+C        0.901534000     -0.475935000     -1.643779000
+C        2.172643000      0.205795000     -1.571490000
+C        2.787859000     -0.146664000     -0.350605000
+C        1.897636000     -1.040418000      0.348138000
+Co       0.793021000      0.711520000      0.089459000
+C        0.117023000      1.122610000      2.027283000
+C        0.794435000      2.210374000      1.353286000
+C       -1.872001000      1.549291000      0.799484000
+C       -1.354374000      1.114183000      1.967798000
+H        3.733106000      0.225487000      0.025111000
+H        2.084335000     -1.480787000      1.320770000
+H       -0.070025000     -1.967917000     -0.261778000
+H       -1.953659000      0.821287000      2.827680000
+H       -2.938621000      1.652487000      0.605915000
+H        1.697202000      2.693944000      1.718266000
+H        0.615909000      0.674598000      2.884193000
+H       -1.982391000      1.435733000     -2.034179000
+N        0.285850000      2.636936000      0.144528000
+C       -0.883841000      1.941356000     -0.233703000
+H       -2.050990000      3.148900000     -1.588712000
+C       -1.394148000      2.265155000     -1.623765000
+H       -0.566933000      2.490576000     -2.302756000
+H        2.553088000      0.908147000     -2.302967000
+H        0.195016000     -0.400597000     -2.462019000