tooba_gmx.py

import os
import subprocess
import numpy as np
from scipy.signal import find_peaks
from scipy.signal import savgol_filter
import hashlib
#cmd="gmx covar -s test/eq_run.tpr -f test/eq_traj.trr -n test/gmx_ndx/test_\(1\,\ 0\,\ 3\).ndx"
#os.system(cmd)

#cmd="gmx anaeig -v test/eigenvec.trr -s test/eq_run.tpr -f test/eq_traj.trr -n test/gmx_ndx/test_\(1\,\ 0\,\ 3\).ndx -proj -first 1 -last 20"
#os.system(cmd)

#cmd="gmx_clusterByFeatures cluster -s test/eq_run.tpr -f test/eq_traj.trr -n test/gmx_ndx/test_\(1\,\ 0\,\ 3\).ndx \
#                                             -feat proj.xvg -method kmeans -nfeature 20 -cmetric ssr-sst -ncluster 5 \
#                                             -fit2central -sort features -cpdb test/central.pdb \
#                                             -fout test/cluster.xtc -plot test/pca_cluster.png"


#os.system(cmd)
def ndx_index(SYSTEM_NAME):
    uniq_id={}
    path='./'+SYSTEM_NAME+'/gmx_ndx/'
    for f in os.listdir(path):
        text=SYSTEM_NAME+'_'+f
        iidd=hashlib.md5(text.encode('utf-8')).hexdigest()
        dom=f.split('_')[2].split('.')[0]
        uniq_id[iidd]={}
        uniq_id[iidd]={'system':SYSTEM_NAME,'ndx_file':f,'domain':dom}
    return uniq_id

def read_xvg(XVG):
    x=[]
    y=[]
    with open(XVG) as f:
        for line in f:
            if line.startswith(("#","@")):
                continue
            else:
                cols = line.split()
                x.append(cols[0])
                y.append(cols[1])
    x = np.asarray(x)
    y = np.asarray(y)
    return x,y

def density_picks(TRR,TPR,IND,SLC,ST,EN,normal,fld,arg,dist_pk=1):

    cmd=str(arg)+'\n'+str(arg)+'\n'
    if EN==-1:
        p = subprocess.Popen(['gmx', 'density', '-symm', '-b', str(ST), '-sl',str(SLC), '-d', normal,
                            '-f', TRR,
                            '-s',TPR,
                            '-n',IND,
                            '-o', fld+'/tmp.xvg'],
                        stdin=subprocess.PIPE)
    else:
        p = subprocess.Popen(['gmx', 'density', '-symm', '-b', str(ST), '-e', str(EN),'-sl',str(SLC), '-d', normal,
                            '-f', TRR,
                            '-s',TPR,
                            '-n',IND,
                            '-o', fld+'/tmp.xvg'],
                        stdin=subprocess.PIPE)
    p.communicate(cmd.encode('UTF-8'))
    p.wait()

    x,y=read_xvg(XVG=fld+'/tmp.xvg')
    yhat = savgol_filter(y, 15, 4) # window size 15, polynomial order 4
    peaks, _ = find_peaks(yhat, distance=dist_pk)

    pathname = os.path.abspath(os.path.join(fld, 'tmp.xvg'))
    os.remove(pathname)

    return peaks #res