Skip to content

Commit

Permalink
preparations for new mass weighting
Browse files Browse the repository at this point in the history
  • Loading branch information
@degiacom
degiacom committed Aug 24, 2016
1 parent 77ab640 commit 0506071
Show file tree
Hide file tree
Showing 2 changed files with 116 additions and 83 deletions.
91 changes: 49 additions & 42 deletions Polymer.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@
import numpy as np
from copy import deepcopy
import logging
from ForceField import ForceField

class Polymer(object):

Expand All @@ -27,61 +28,37 @@ def __init__(self,db,ff,molname,mode):
self.chain=""
self.search_grid=self._make_search_grid()

self.mass=0

self.logger=logging.getLogger('assemble')


'''
#in gromacs mode, calculate mass of chain
def get_mass(self):

if len(self.chain)==0:
raise Exception("no chain provided!")

mass=0

try:
for c in self.chain:
mol=self.db.molecules[c]
mapping=mol.topology.mapping
for l in mol.data[:,1]:
atomtype=mapping[np.where(mapping==l)[0],1]
mass+=self.ff.nonbonded[atomtype][1]
for c in self.chain:
mol=self.db.molecules[c]
mapping=mol.topology.mapping
for l in mol.data[:,0]:
atomname=mol.atom.keys()[np.where(mol.atom.values()==l)[0]]

return mass
except:
raise Exception("could not calculate polymer mass!")
try:
atomtype=mapping[np.where(mapping==atomname)[0],1][0]
thismass=self.ff.nonbonded[atomtype][1]
mass+=float(thismass)

except Exception, e:
raise Exception("Could not find mass of atom %s"%atomname)

self.mass=mass
return mass
'''

def get_mass_2(self):

masslist=[]
for j in xrange(0,len(self.poly),1):

data_list=self.poly[j].mapping(self.poly[j].data)

#get topology information of current molecule
top=self.poly[j].topology

#if molecule is terminal, modify its topology accordingly
if j==0:
top.make_terminal("nterminal")
if j==len(self.poly)-1:
top.make_terminal("cterminal")

for i in xrange(0,len(data_list),1):

atomtype=top.mapping[top.mapping[:,0]==data_list[i][1],1][0]
mass=self.ff.nonbonded[atomtype][1]
masslist.append(float(mass))

polymass=sum(masslist)
return(polymass)





def make(self,chain):

#debug=0 #debug mode prints hooks positions in a separate file
Expand All @@ -91,6 +68,15 @@ def make(self,chain):
self.logger.info("\n> generating polymer %s..."%self.molname)
self.logger.info(">> sequence: %s"%self.chain)

#if gromacs mode: calculate mass
if self.mode=="gromacs":
try:
self.get_mass()
except Exception, e:
raise Exception("%s for polymer %s"%(e, self.molname))
self.logger.info(">> mass: %s Da"%self.mass)


#add first element in the chain
m=deepcopy(self.db.molecules[self.chain[0]])
self.poly.append(m)
Expand Down Expand Up @@ -849,3 +835,24 @@ def _random_orthonormal(self,normal):
v /= np.linalg.norm(v)
alpha = np.random.uniform(0.0, np.pi*2)
return np.cos(alpha)*u + np.sin(alpha)*v



if __name__=="__main__":

import os,sys
cwd=os.getcwd()
assembled=os.path.abspath(os.path.dirname(str(sys.argv[0])))
os.environ["ASSEMBLEPATH"]="%s;%s"%(cwd,assembled)

#test mass estimation
from Database import Database
D=Database()
D.load("database\\database.txt", "gromacs")
F=ForceField()
F.load("database\\forcefield\\trappe.ff.txt")
P=Polymer(D,F,"test","gromacs")
P.chain="ccctttccCCt"

print P.get_mass()

108 changes: 67 additions & 41 deletions System.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,28 +24,41 @@ def __init__(self,polymers,ff,params):
self.logger=logging.getLogger('assemble')

#generate a random distribution of molecules in a box, given their percentages
def make_box(self, dim, percentage):
def make_box(self, dim, data_in, use_fractional_mass):

#prepare data
names=[]
v=[]
for d in data_in:
names.append(d[0])
v.append(d[1])

#if reference percentage is a fractional mass, tweak it using monomers mass
if use_fractional_mass:
percentage=self.convert_concentration(np.array(v), get_fractional=True)
else:
percentage=np.array(v)

#normalize to 100 percent
test=0
for l in percentage:
test+=float(l[1])
test+=float(l)

for i in xrange(len(percentage)):
v=float(percentage[i][1])
percentage[i][1]=(v/test)*100.0
#v=float(percentage[i])
percentage[i]=(percentage[i]/test)*100.0

if len(dim)!=3:
raise IOError("ERROR: expected 3 values for dimensions")

length=dim[0]*dim[1]*dim[2]

#make roulette array
roulette=[percentage[0][1]]
t=[percentage[0][1]]
roulette=[percentage[0]]
t=[percentage[0]]
for r in xrange(1,len(percentage),1):
t.append(percentage[r][1])
roulette.append(roulette[r-1]+percentage[r][1])
t.append(percentage[r])
roulette.append(roulette[r-1]+percentage[r])

target=np.array(t)
error=100000.0
Expand All @@ -66,7 +79,7 @@ def make_box(self, dim, percentage):
break
else:
index+=1
c.append(percentage[index][0])
c.append(names[index])
counter[index]+=1

#chain completed, report produced percentages:
Expand All @@ -79,18 +92,58 @@ def make_box(self, dim, percentage):
error=tmp_err
cbest=np.array(c)
counterbest=counter.copy()

print counterbest

#if reference percentage is a fractional mass, tweak it using monomers mass
if use_fractional_mass:
perc=self.convert_concentration(counterbest, get_fractional=False) #convert into percentage
perc/=np.sum(perc)
perc*=100.0
mw=counterbest.copy() #this is already a fractional mass
else:
perc=counterbest.copy()
mw=self.convert_concentration(counterbest, get_fractional=True)
mw/=np.sum(mw)
mw*=100.0

self.logger.info("> polymer concentrations in box:")
self.logger.info("> polymer units in box:")
for x in xrange(0,len(counterbest),1):
self.logger.info(">> %s: %s percent"%(percentage[x][0],counterbest[x]))

self.logger.info(">> %s: %s percent"%(names[x],perc[x]))


self.logger.info("\n> polymers fractional weight:")
for x in xrange(0,len(counterbest),1):
self.logger.info(">> %s: %s percent"%(names[x],mw[x]))


return np.reshape(cbest,dim)


#if get_fractional=True, percentage is edited using monomers mass
#if get_fractional=False, fractional mass is converted into percentage
def convert_concentration(self, c, get_fractional=True):
m=[]
for x in xrange(0,len(self.polymers),1):
m.append(self.polymers[x].mass)
allmasses=np.array(m)
allmasses/=np.sum(allmasses)

for x in xrange(0,len(allmasses),1):
if get_fractional:
rescaled=c[x]/allmasses[x]
else:
rescaled=c[x]*allmasses[x]

c[x]=rescaled

return c


def create_system(self,mypath="."):

### CREATE BOX ###
self.systembox=self.make_box(self.params.box_grid_shape, self.params.concentration)
### CREATE BOX ###
self.systembox=self.make_box(self.params.box_grid_shape, self.params.concentration, use_fractional_mass=False)

#get maximal box between existing polymers, to define voxel size
voxel_size=np.array([0.,0.,0.])
Expand Down Expand Up @@ -169,30 +222,6 @@ def create_system(self,mypath="."):
cnt=len(self.systembox[self.systembox==name])
contmols.append([name,cnt])


# print weight concentration
masspoly=[]
for x in xrange(0,len(self.polymers),1):
masspoly.append(self.polymers[x].get_mass_2())

weighted=[]
for x in xrange(0,len(masspoly),1):
weighted.append(float(masspoly[x])*float(contmols[x][1]))

weightfrac=[]
for x in xrange(0,len(masspoly),1):
weightfrac.append(weighted[x]/sum(weighted))

self.logger.info("\n> polymer weight concentrations in box:")
for x in xrange(0,len(contmols),1):
self.logger.info(">> %s: %s weight percent"%(contmols[x][0],weightfrac[x]))

# print number of molecules in box
self.logger.info("\n> number of molecules in box:")
for x in xrange(0,len(contmols),1):
self.logger.info(">> %s: %s "%(contmols[x][0],contmols[x][1]))


#counters needed for index file creation
index_full=1 #atom counter without wrapping
index_per_mol=1 #atom counter per molecule type (for pretty formatting)
Expand Down Expand Up @@ -257,9 +286,6 @@ def create_system(self,mypath="."):
maxbox=np.max(np.array(maxpos),axis=0)
box=maxbox-minbox

# print total number of beads in box
self.logger.info("\n> total number of beads in box: %s", index_full-1)

self.logger.info("\n> box size: %10.5f x %10.5f x %10.5f nm^3"%(box[0],box[1],box[2]))

### BOX INFORMATION TO ADD! ###
Expand Down

0 comments on commit 0506071

Please sign in to comment.