From c5fc65433aee509f25096ae78113d5b30eccceed Mon Sep 17 00:00:00 2001
From: Oliver Henrich <oliver.henrich@strath.ac.uk>
Date: Thu, 28 Apr 2022 19:00:40 +0100
Subject: [PATCH] Updated and added utility scripts

---
 examples/PACKAGES/cgdna/util/lmp2vis.py |  57 ++++-----
 examples/PACKAGES/cgdna/util/nbps.py    | 161 ++++++++++++++++++++++++
 2 files changed, 190 insertions(+), 28 deletions(-)
 create mode 100755 examples/PACKAGES/cgdna/util/nbps.py

diff --git a/examples/PACKAGES/cgdna/util/lmp2vis.py b/examples/PACKAGES/cgdna/util/lmp2vis.py
index d51db447ce..9ce856b745 100644
--- a/examples/PACKAGES/cgdna/util/lmp2vis.py
+++ b/examples/PACKAGES/cgdna/util/lmp2vis.py
@@ -37,9 +37,9 @@ import sys, math, subprocess
 # converts quaternion DOF into local body reference frame
 def q_to_exyz(q1,q2,q3,q4):
 
-    q=[q1, q2, q3, q4]
-    ex= [0, 0, 0]
-    ey = [0, 0,0]
+    q = [q1, q2, q3, q4]
+    ex = [0, 0, 0]
+    ey = [0, 0, 0]
     ez = [0, 0, 0]
 
     ex[0]=q[0]*q[0]+q[1]*q[1]-q[2]*q[2]-q[3]*q[3]
@@ -57,22 +57,22 @@ def q_to_exyz(q1,q2,q3,q4):
     return ex,ey,ez
 
 # processes line by line of LAMMPS trajectory
-def transform(line): 
+def transform(line, colind): 
 
     list1 = line.split()
-    ident, mol, typ = int(list1[0]), int(list1[1]), int(list1[2])
+    ident, mol, typ = int(list1[colind[0]]), int(list1[colind[1]]), int(list1[colind[2]])
     typb = typ + 10  # defines new backbone types, here offset is 10 from atom (base) type
-    x, y, z = float(list1[3]), float(list1[4]), float(list1[5])
-    vx, vy, vz = float(list1[9]), float(list1[10]), float(list1[11])
+    x, y, z = float(list1[colind[3]]), float(list1[colind[4]]), float(list1[colind[5]])
+    vx, vy, vz = float(list1[colind[6]]), float(list1[colind[7]]), float(list1[colind[8]])
     c_quat1, c_quat2, c_quat3, c_quat4 = \
-            float(list1[12]), float(list1[13]), float(list1[14]), float(list1[15])
+            float(list1[colind[9]]), float(list1[colind[10]]), float(list1[colind[11]]), float(list1[colind[12]])
 
     ex, ey, ez = q_to_exyz(c_quat1, c_quat2, c_quat3, c_quat4)
 
     # position of sugar-phosphate backbone interaction site in oxDNA2
     x1, y1, z1 = x -0.34*ex[0]+0.3408*ey[0],y -0.34*ex[1]+0.3408*ey[1], z-0.34*ex[2]+0.3408*ey[2]
 
-    # position of base interaction site in oxDNA2
+    # position of base interaction site in oxDNA/oxDNA2
     x2, y2, z2 = x +0.4*ex[0], y + 0.4*ex[1], z+0.4*ex[2]
 
     # compose basic output data: id, molecule id, type, position, velocity quaternion
@@ -96,7 +96,7 @@ def transform(line):
     line1 += '\n'
     line2 += '\n'
 
-    line= line1 + line2
+    line = line1 + line2
     return line
 
 ### main part ###
@@ -132,24 +132,25 @@ except:
 r=open(infilename,'r')
 w=open(outfilename,'w+')
 
-line=r.readline() # read first line in file
+pass1 = 0
 
-while line != '':
+for line in r:
 
     sys.stdout.write('# Processed %3d %%\r' % (100*n/nlines))     
 
-    # find number of atoms in timestep and double
-    if line.find('NUMBER OF ATOMS') != -1: 
+    # find number of atoms in timestep
+    if line.find('ITEM: NUMBER OF ATOMS') != -1: 
         w.write(line)
         N=int(r.readline())
         # write to output file and read next line
         w.write('%d'%int(2*N)+'\n')
         line=r.readline()
+        n+=2
 
     # find beginning of atom data section
     if line.find('ITEM: ATOMS') != -1:
-        # first pass: extract column number of ID, molecule ID, type, postion, velocity, quaternion
-        if n==0:
+        # first pass: extract column number of ID, molecule ID, type, postion, velocity, quaternion in header line
+        if pass1 == 0:
             linestring=line.split()
             idindex = linestring.index('id')
             molindex = linestring.index('mol')
@@ -177,24 +178,24 @@ while line != '':
                 header += ' shape[0] shape[1] shape[2]'
             header += '\n'
 
-        ### begin processing atom data
-        i=0
+            # store column number in data line, -2 offset form header line
+            colind = [idindex-2,molindex-2,typeindex-2,xindex-2,yindex-2,zindex-2,\
+                      vxindex-2,vyindex-2,vzindex-2,qwindex-2,qxindex-2,qyindex-2,qzindex-2]
+            pass1 = 1
+
+        # begin processing atom data in timestep
         w.write(header)
-
         # tranform each atom and write to output file
-        while i<N:
+        for i in range(N):
             line=r.readline()
-            w.write(transform(line))
-            i +=1
-            n +=1
-        else:
-            line=r.readline() # read next line
-        ### end processing atom data
+            w.write(transform(line,colind))
+            n+=1
+        # end processing atom data
 
-    # duplicate all lines that are not in atom data section 
+    # duplicate all lines that are not in any other section 
     else:
         w.write(line)
-        line=r.readline()
+        n+=1
 
 print('# Done                                ')     
 
diff --git a/examples/PACKAGES/cgdna/util/nbps.py b/examples/PACKAGES/cgdna/util/nbps.py
new file mode 100755
index 0000000000..16c118cf4f
--- /dev/null
+++ b/examples/PACKAGES/cgdna/util/nbps.py
@@ -0,0 +1,161 @@
+#!/usr/bin/env python
+"""
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/ Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+   Contributing author: Oliver Henrich (University of Strathclyde, Glasgow)
+------------------------------------------------------------------------- */
+
+Program:
+
+Usage: 
+$$ python 
+
+Requirements:
+The LAMMPS trajectory input file needs to contain the following data columns:
+id mol type x y z c_quat[1] c_quat[2] c_quat[3] c_quat[4]
+"""
+
+import sys, math, subprocess
+
+# converts quaternion DOF into local body reference frame
+def q_to_exyz(q1,q2,q3,q4):
+
+    q = [q1, q2, q3, q4]
+    ex = [0, 0, 0]
+    ey = [0, 0, 0]
+    ez = [0, 0, 0]
+
+    ex[0]=q[0]*q[0]+q[1]*q[1]-q[2]*q[2]-q[3]*q[3]
+    ex[1]=2*(q[1]*q[2]+q[0]*q[3])
+    ex[2]=2*(q[1]*q[3]-q[0]*q[2])
+
+    ey[0]=2*(q[1]*q[2]-q[0]*q[3])
+    ey[1]=q[0]*q[0]-q[1]*q[1]+q[2]*q[2]-q[3]*q[3]
+    ey[2]=2*(q[2]*q[3]+q[0]*q[1])
+
+    ez[0]=2*(q[1]*q[3]+q[0]*q[2])
+    ez[1]=2*(q[2]*q[3]-q[0]*q[1])
+    ez[2]=q[0]*q[0]-q[1]*q[1]-q[2]*q[2]+q[3]*q[3]
+
+    return ex,ey,ez
+
+# processes line by line of LAMMPS trajectory
+def transform(line, colind): 
+
+    list1 = line.split()
+    ident, mol, typ = int(list1[colind[0]]), int(list1[colind[1]]), int(list1[colind[2]])
+    x, y, z = float(list1[colind[3]]), float(list1[colind[4]]), float(list1[colind[5]])
+    c_quat1, c_quat2, c_quat3, c_quat4 = \
+            float(list1[colind[6]]), float(list1[colind[7]]), float(list1[colind[8]]), float(list1[colind[9]])
+
+    ex, ey, ez = q_to_exyz(c_quat1, c_quat2, c_quat3, c_quat4)
+
+    # position of base interaction site in oxDNA/oxDNA2
+    x2, y2, z2 = x +0.4*ex[0], y + 0.4*ex[1], z+0.4*ex[2]
+
+    # compose basic output data: id, molecule id, type, position
+    line2 = [ident, mol, typ, x2, y2, z2]
+
+    return line2
+
+### main part ###
+
+# digest command line input
+if len(sys.argv)<3:
+    print("Syntax: $$ python nbps.py input_filename output_filename")
+    sys.exit(1)
+
+if len(sys.argv)==3:
+    infilename  = sys.argv[1]
+    outfilename = sys.argv[2]
+
+print('# Calculating number of base paris from LAMMPS trajectory output')     
+
+# count lines to process for progress report
+n = 0
+try:
+    result = subprocess.run(['wc', '-l', '%s'%infilename], stdout=subprocess.PIPE)
+    reslist=str(result).split()
+    nlines=float(reslist[5])
+except:
+    nlines = 100
+
+r=open(infilename,'r')
+w=open(outfilename,'w+')
+
+pass1 = 0
+
+for line in r:
+
+    sys.stdout.write('# Processed %3d %%\r' % (100*n/nlines))     
+
+   # find timestep number
+    if line.find('ITEM: TIMESTEP') != -1:
+        t=int(r.readline())
+        n+=1
+
+   # find number of atoms in timestep
+    if line.find('ITEM: NUMBER OF ATOMS') != -1:
+        N=int(r.readline())
+        n+=1
+
+    # find beginning of atom data section
+    if line.find('ITEM: ATOMS') != -1:
+        # first pass: extract column number of ID, molecule ID, type, postion, velocity, quaternion in header line
+        if pass1 == 0:
+            linestring=line.split()
+            idindex = linestring.index('id')
+            molindex = linestring.index('mol')
+            typeindex = linestring.index('type')
+            xindex = linestring.index('x')
+            yindex = linestring.index('y')
+            zindex = linestring.index('z')
+            qwindex = linestring.index('c_quat[1]')
+            qxindex = linestring.index('c_quat[2]')
+            qyindex = linestring.index('c_quat[3]')
+            qzindex = linestring.index('c_quat[4]')
+
+            # store column number in data line, -2 offset form header line
+            colind = [idindex-2,molindex-2,typeindex-2,xindex-2,yindex-2,zindex-2,\
+                      qwindex-2,qxindex-2,qyindex-2,qzindex-2]
+            pass1 = 1
+
+        # begin processing atom data
+        atom = []
+        for i in range(N):
+            line=r.readline()
+            atom.append(transform(line,colind))
+            n +=1
+
+        rc2 = 0.2
+        nbp = 0
+        for i in range(N): 
+            for j in range(i+1,N): 
+                if (atom[i][2]+atom[j][2])%4 and atom[i][1]!=atom[j][1]:
+                    r2 = (atom[i][3]-atom[j][3])*(atom[i][3]-atom[j][3]) +\
+                          (atom[i][4]-atom[j][4])*(atom[i][4]-atom[j][4]) +\
+                          (atom[i][5]-atom[j][5])*(atom[i][5]-atom[j][5])
+                    if r2<rc2:
+
+                        nbp += 1
+        w.write('%d %d\n' % (t,nbp))
+        # end processing atom data
+
+    n+=1
+
+print('# Done                                ')     
+
+r.close()
+w.close()