From b76628020e60cf48332664c8590782c847cf28b2 Mon Sep 17 00:00:00 2001 From: Michael Humbert Date: Thu, 14 Apr 2016 15:06:44 -0400 Subject: [PATCH] removed RadialDistributionPure for properties_abandon and edited Radial Distribution to remove fortran and perform calculations using vectorization. --- rubicon/analysis/lammps/properties_abandon.py | 127 ++++++------------ 1 file changed, 39 insertions(+), 88 deletions(-) diff --git a/rubicon/analysis/lammps/properties_abandon.py b/rubicon/analysis/lammps/properties_abandon.py index d2a88a3..5273761 100644 --- a/rubicon/analysis/lammps/properties_abandon.py +++ b/rubicon/analysis/lammps/properties_abandon.py @@ -692,11 +692,10 @@ def runradial(self, datfilename, mol1, mol2, comx, comy, comz, Lx, Ly, Lz, firststep) (nummol1, nummol2, mol1, mol2) = self.getnummol(moltypel, nummoltype, mol1, mol2) - while count < len(comx): - (count) = self.radialdistribution(g, mol1, mol2, len(comx[1]), - moltype, comx, comy, comz, Lx, - Ly, Lz, binsize, numbins, maxr, - count) + + (count) = self.radialdistribution(g, len(comx[1]), moltype, comx, comy, + comz, Lx, Ly, Lz, binsize, numbins, + maxr, count) (radiuslist) = self.radialnormalization(numbins, binsize, Lx, Ly, Lz, nummol1, nummol2, count, g, firststep) @@ -721,15 +720,41 @@ def getnummol(self, moltypel, nummoltype, mol1, mol2): mol2 = int(moltypel.index(mol2)) return nummol1, nummol2, mol1, mol2 - def radialdistribution(self, g, mol1, mol2, nummol, moltype, comx, comy, + def radialdistribution(self, g, nummol, moltype, comx, comy, comz, Lx, Ly, Lz, binsize, numbins, maxr, count): - # uses FORTRAN code to calculate the number of molecules within each - # shell - g1 = comradial.comradial(mol1, mol2, nummol, moltype, comx[count], - comy[count], comz[count], Lx, Ly, Lz, binsize, - numbins, maxr) - g += g1 - count += 1 + # calculates the number of molecules within each shell + + comxt = np.array(comx[count:]).transpose() + comyt = np.array(comy[count:]).transpose() + comzt = np.array(comz[count:]).transpose() + indexlist = [] + + for i in range(0,len(g)): + indexlist.append(np.array(moltype) == i) + + + for molecule in range(0, nummol - 1): + dx = comxt[molecule+1:] - np.tile(comxt[molecule], + (len(comxt)-molecule-1,1)) + dy = comyt[molecule+1:] - np.tile(comyt[molecule], + (len(comyt)-molecule-1,1)) + dz = comzt[molecule+1:] - np.tile(comzt[molecule], + (len(comzt)-molecule-1,1)) + + dx -= Lx * np.around(dx / Lx) + dy -= Ly * np.around(dy / Ly) + dz -= Lz * np.around(dz / Lz) + + r2 = dx ** 2 + dy ** 2 + dz ** 2 + r = np.sqrt(r2) + for i in range(0,len(indexlist)): + gt,dist = np.histogram(r[indexlist[i][molecule+1:]].ravel(), + bins=numbins, + range=(0.5*binsize,binsize*(numbins+0.5))) + g[moltype[molecule]][i]+= gt + g[i][moltype[molecule]]+= gt + + count = len(comx) return count def radialnormalization(self, numbins, binsize, Lx, Ly, Lz, nummol1, @@ -737,7 +762,7 @@ def radialnormalization(self, numbins, binsize, Lx, Ly, Lz, nummol1, # normalizes g to box density radiuslist = (np.arange(numbins) + 1) * binsize g *= Lx * Ly * Lz / nummol1 / nummol2 / 4 / np.pi / ( - radiuslist) ** 2 / binsize / ( + radiuslist) ** 2 / binsize / ( count - firststep) return radiuslist @@ -975,77 +1000,3 @@ def writetofile(self, numbins, radiuslist, g): str(radiuslist[radius]) + " " + str(g[radius]) + "\n") gfile.close() - -class RadialDistributionPure: - def runradial(self, datfilename, comx, comy, comz, Lx, Ly, Lz, Lx2, Ly2, - Lz2, output, nummoltype, moltypel, moltype, firststep=1): - (maxr, binsize, numbins, count, g) = self.setgparam(Lx2, Ly2, Lz2, - firststep, - moltypel) - while count < len(comx): - (count) = self.radialdistribution(g, len(comx[1]), moltype, comx, - comy, comz, Lx, Ly, Lz, binsize, - numbins, maxr, count) - (radiuslist) = self.radialnormalization(numbins, binsize, Lx, Ly, Lz, - nummoltype, count, g, - firststep) - self.append_dict(radiuslist, g, output, moltypel) - return output - - def setgparam(self, Lx2, Ly2, Lz2, firststep, moltypel): - # uses side lengths to set the maximum radius for box and number of bins - # also sets the first line using data on firststep and number of atoms - maxr = min(Lx2, Ly2, Lz2) - binsize = 0.1 - numbins = int(np.ceil(maxr / binsize)) - count = firststep - g = np.zeros((len(moltypel), len(moltypel), numbins)) - return maxr, binsize, numbins, count, g - - def radialdistribution(self, g, nummol, moltype, comx, comy, comz, Lx, Ly, - Lz, binsize, numbins, maxr, count): - # implements a FORTRAN code to calculate the number of molecules within each shell - - for molecule1 in range(0, nummol - 1): - for molecule2 in range(molecule1 + 1, nummol): - dx = comx[count][molecule1] - comx[count][molecule2] - dy = comy[count][molecule1] - comy[count][molecule2] - dz = comz[count][molecule1] - comz[count][molecule2] - - dx -= Lx * round(dx / Lx) - dy -= Ly * round(dy / Ly) - dz -= Lz * round(dz / Lz) - - r2 = dx ** 2 + dy ** 2 + dz ** 2 - r = np.sqrt(r2) - - if r <= maxr: - g[moltype[molecule1]][moltype[molecule2]][ - int(round(r / binsize) - 1)] += 1 - g[moltype[molecule2]][moltype[molecule1]][ - int(round(r / binsize) - 1)] += 1 - - count += 1 - return count - - def radialnormalization(self, numbins, binsize, Lx, Ly, Lz, nummol, count, - g, firststep): - # normalizes g to box density - radiuslist = (np.arange(numbins) + 1) * binsize - for i in range(0, len(g)): - for j in range(0, len(g)): - g[i][j] *= Lx * Ly * Lz / nummol[i] / nummol[j] / 4 / np.pi / ( - radiuslist) ** 2 / binsize / ( - count - firststep) - return radiuslist - - def append_dict(self, radiuslist, g, output, moltypel): - for i in range(0, len(moltypel)): - for j in range(i, len(moltypel)): - if not all([v == 0 for v in g[i][j]]): - output['RDF']['{0}-{1}'.format(moltypel[i], - moltypel[ - j])] = copy.deepcopy( - g[i][j].tolist()) - if 'distance' not in list(output['RDF'].keys()): - output['RDF']['distance'] = copy.deepcopy(radiuslist.tolist()) \ No newline at end of file