<gh_stars>0
import numpy as np
#from matplotlib import pyplot as plt

'''
This is a code that takes some of your  atomistic simulation snapshots (frames/*.gro), 
looks  up  in  polyIndex  which  atomistic  molecules  correspond  to  a  CG bead, and
calculates the CG  centre  of mass of the CG bead and outputs it to frames2/*.gro.

This can be used to calculate angular distributions or bond distance distributions.

You have to specify how many CG segments exist in each molecule in nMono, and also how 
many atomistic frames you want to generate.
'''

numframe  = 40
numMono   = 2
listPairs = []
f = open("polyIndex", 'r')
for line in f:
    words = line.split()
    for i in range(len(words)):
        words[i] = int(words[i])

    listPairs.append(words)
f.close()


##############################################

class Pairs:
    def __init__(self, filename, pair, data):
        self.pair = pair
        self.posit = []
        framefile = open(filename, 'r')
        lineList = framefile.readlines()
        self.boxDim = lineList[len(lineList)-1].split()
        self.boxDim = float(self.boxDim[0])

        for i in pair:
            self.posit.append(data[i])

        for j in range(0, len(self.posit)):
            for i in range(0, len(self.posit[j])):
                if self.posit[j][i] - self.posit[0][i] < -1*self.boxDim/2:
                    self.posit[j][i] = self.posit[j][i] + self.boxDim
                elif self.posit[j][i] - self.posit[0][i] > self.boxDim/2:
                    self.posit[j][i] = self.posit[j][i] - self.boxDim

        cascade = np.array([])
        for i in range(0, len(self.posit)):
            if i == 0:
                cascade = np.array(self.posit[i])
            else:
                cascade = cascade + np.array(self.posit[i])
        self.positAv = cascade/len(pair)
        self.dist = dist(self.posit[0], self.posit[1])


def dist(vec1, vec2):
    dist = 0
    for i in range(0, len(vec1)):
        dist = dist + (vec1[i]-vec2[i])**2
    return dist**.5


def extractor(filename):
    f=open(filename, 'r')
    pos=[]
    f.readline()
    nPart=int(f.readline())
    for i in range(nPart):
        words = f.readline().split()
        if len(words) > 7:
            x = float(words[-6])
            y = float(words[-5])
            z = float(words[-4])
        else:
            x = float(words[-3])
            y = float(words[-2])
            z = float(words[-1])
        pos.append([x,y,z])
    f.close()
    return pos



###########################################


for i in range(0,int(numframe)):
        filename = "frames/%d.gro" % i
        filename2 = "frames2/%d.gro" % i
        data = extractor(filename)
        outFile  =  open("frames2/%d.gro" % i, 'w')
        pairList = np.ndarray.tolist(np.copy(listPairs))

        outFile.write('Generated by trjconv : MAZCHAIN Pv at 617K t= %11.5f\n %i\n' % (i, len(listPairs)))
        molNum = 0
        for j in range(0, len(pairList)):
            pairList[j] = Pairs(filename, pairList[j], data)
            if j%numMono == 0:
                molNum = molNum + 1 

            outFile.write('%5d%-5s%5s%5d%8.3f%8.3f%8.3f' %(molNum,'PLE','PL',(j+1),pairList[j].positAv[0],pairList[j].positAv[1],pairList[j].positAv[2]))
            outFile.write('\n')

        inFile = open("frames/%d.gro" % i, 'r')
        linez = inFile.readlines()
        linez = linez[len(linez)-1]

        outFile.write(linez)
        print(i+1, 'done')
        outFile.close()


import numpy as np
#from matplotlib import pyplot as plt

'''
This is a code that takes some of your  atomistic simulation snapshots (frames/*.gro), 
looks  up  in  polyIndex  which  atomistic  molecules  correspond  to  a  CG bead, and
calculates the CG  centre  of mass of the CG bead and outputs it to frames2/*.gro.

This can be used to calculate angular distributions or bond distance distributions.

You have to specify how many CG segments exist in each molecule in nMono, and also how 
many atomistic frames you want to generate.
'''

numframe  = 40
numMono   = 2
listPairs = []
f = open("polyIndex", 'r')
for line in f:
    words = line.split()
    for i in range(len(words)):
        words[i] = int(words[i])

    listPairs.append(words)
f.close()


##############################################

class Pairs:
    def __init__(self, filename, pair, data):
        self.pair = pair
        self.posit = []
        framefile = open(filename, 'r')
        lineList = framefile.readlines()
        self.boxDim = lineList[len(lineList)-1].split()
        self.boxDim = float(self.boxDim[0])

        for i in pair:
            self.posit.append(data[i])

        for j in range(0, len(self.posit)):
            for i in range(0, len(self.posit[j])):
                if self.posit[j][i] - self.posit[0][i] < -1*self.boxDim/2:
                    self.posit[j][i] = self.posit[j][i] + self.boxDim
                elif self.posit[j][i] - self.posit[0][i] > self.boxDim/2:
                    self.posit[j][i] = self.posit[j][i] - self.boxDim

        cascade = np.array([])
        for i in range(0, len(self.posit)):
            if i == 0:
                cascade = np.array(self.posit[i])
            else:
                cascade = cascade + np.array(self.posit[i])
        self.positAv = cascade/len(pair)
        self.dist = dist(self.posit[0], self.posit[1])


def dist(vec1, vec2):
    dist = 0
    for i in range(0, len(vec1)):
        dist = dist + (vec1[i]-vec2[i])**2
    return dist**.5


def extractor(filename):
    f=open(filename, 'r')
    pos=[]
    f.readline()
    nPart=int(f.readline())
    for i in range(nPart):
        words = f.readline().split()
        if len(words) > 7:
            x = float(words[-6])
            y = float(words[-5])
            z = float(words[-4])
        else:
            x = float(words[-3])
            y = float(words[-2])
            z = float(words[-1])
        pos.append([x,y,z])
    f.close()
    return pos



###########################################


for i in range(0,int(numframe)):
        filename = "frames/%d.gro" % i
        filename2 = "frames2/%d.gro" % i
        data = extractor(filename)
        outFile  =  open("frames2/%d.gro" % i, 'w')
        pairList = np.ndarray.tolist(np.copy(listPairs))

        outFile.write('Generated by trjconv : MAZCHAIN Pv at 617K t= %11.5f\n %i\n' % (i, len(listPairs)))
        molNum = 0
        for j in range(0, len(pairList)):
            pairList[j] = Pairs(filename, pairList[j], data)
            if j%numMono == 0:
                molNum = molNum + 1 

            outFile.write('%5d%-5s%5s%5d%8.3f%8.3f%8.3f' %(molNum,'PLE','PL',(j+1),pairList[j].positAv[0],pairList[j].positAv[1],pairList[j].positAv[2]))
            outFile.write('\n')

        inFile = open("frames/%d.gro" % i, 'r')
        linez = inFile.readlines()
        linez = linez[len(linez)-1]

        outFile.write(linez)
        print(i+1, 'done')
        outFile.close()
#from matplotlib import pyplot as plt This is a code that takes some of your  atomistic simulation snapshots (frames/*.gro),  looks  up  in  polyIndex  which  atomistic  molecules  correspond  to  a  CG bead, and calculates the CG  centre  of mass of the CG bead and outputs it to frames2/*.gro.  This can be used to calculate angular distributions or bond distance distributions.  You have to specify how many CG segments exist in each molecule in nMono, and also how  many atomistic frames you want to generate. ############################################## ###########################################