import numpy as np
from cpolymer.polymer import Polymer
from cpolymer.lsimu import LSimu
from cpolymer.constrain import Box,Sphere,Point
from cpolymer.halley.constrain import Spherical,Nowhere
from cpolymer.halley.vectors import V
len_chrom = [46, 162, 63, 306, 115, 54, 218, 112, 87, 149, 133, 365, 184, 156, 218, 189]
dist_centro = [30, 47, 22, 89, 30, 29, 99, 21, 71, 87, 88, 30, 53, 125, 65, 111]
Radius = 16.6
Mt=0.3*16.6
Nchromosomes = 16
dnuc =3
def create_nucleus(simple=True,oSim=None,dnuc=3,boxf=1.1,Radius=16.6,quad=1):
#If simple is True,
#create one copy of the chromosomes,
#else create two copies
Sim = LSimu()
nucleus = Sphere(position=[0,0,0],radius=Radius)
bead_type = 1 # All the beads are going to be the same type
liaison = {"1-1":[1,1],"1-2":[1,2],"1-3":[1,3],"1-4":[(dnuc+1)/2.,4],"1-5":[0,5],
"2-2":[1,6],"2-3":[1,7],"2-4":[(dnuc+1)/2.,8],"2-5":[0,9],
"3-3":[1,10],"3-4":[(dnuc+1)/2.,11],"3-5":[Mt,12],
"4-4":[dnuc,13],"4-5":[0,14],
"5-5":[0,15]}
fact = 1
if not simple:
fact = 2
for X in range(fact*Nchromosomes):
#We need to define geometrical constrain:
#The centromere must be at a distance mt of the spb positioned at (-Radius,0,0)
#We use the module halley for that
S_spb = Spherical(V(-Radius,0,0),radius=Mt) #we define a sphere centered on the spb
S_centered = Spherical(V(0,0,0),radius=Radius-Mt*0.9) # a sphere centered that intersect the spb sphere
circle = S_spb * S_centered
centromere = circle.get_random()
#We must then construct a sphere centered on the centromere with a radius sqrt(Nbead)
#and look at its intersection with the nucleus
Nucleus = Spherical(V(0,0,0),radius=Radius*0.95) # a sphere centered that intersect the spb sphere
d1 = dist_centro[X % Nchromosomes]
Telo1_possible = Spherical(centromere,radius=np.sqrt(d1)) * Nucleus
telo1 = Telo1_possible.get_random()
d2 = len_chrom[X % Nchromosomes]-dist_centro[X % Nchromosomes]
Telo2_possible = Spherical(centromere,radius=np.sqrt(d2)) * Nucleus
telo2 = Telo1_possible.get_random()
if X != 11 and X != 11+Nchromosomes:
Sim.add(Polymer(N=len_chrom[X % Nchromosomes],type_bead=[2]+[1]*(d1-2)+[3]+[1]*(d2-1) + [2],
liaison=liaison,
gconstrain=[nucleus],
lconstrain=[Point(index=0,position=telo1._v),
Point(index=d1,position=centromere._v),
Point(index=len_chrom[X % Nchromosomes]-1,position=telo2._v)]))
else:
# This chromosome is different because it has a nucleole
delta = 0
if X != 11:
delta = 1 # to avoid the superposition of the center of the nucleole
Sim.add(Polymer(N=len_chrom[X % Nchromosomes],type_bead=[2]+[1]*(d1-2)+[3]+[1]*(90-d1)+[4]*150+\
[1]*(len_chrom[X % Nchromosomes]-150-90-1) + [2],
liaison=liaison,
gconstrain=[nucleus],
lconstrain=[Point(index=0,position=telo1._v),
Point(index=d1,position=centromere._v),
Point(index=90+75,position=(0.66*Radius,delta,0)), #We add a new constrain: the center
#of the nucleole must be at 2/3 of
#the radius at the
#opposite of the spb:
Point(index=len_chrom[X % Nchromosomes]-1,position=telo2._v)]))
#Then Add the spb
Sim.add(Polymer(N=1,type_bead=5,liaison=liaison))
Sim.molecules[-1].coords = numpy.array([[-Radius,0,0]])
if not simple:
Sim.add(Polymer(N=1,type_bead=5,liaison=liaison))
Sim.molecules[-1].coords = numpy.array([[-Radius,0,0]])
for i,c in enumerate(dist_centro,1):
if c is not None:
Sim.add_extra_bond(mol1=[len(Sim.molecules),1],mol2=[i,c],typeb=liaison["3-5"][1])
if not simple:
Sim.add_extra_bond(mol1=[len(Sim.molecules)-1,1],mol2=[i+Nchromosomes,c],typeb=liaison["3-5"][1])
#This create the extra bond between the centromeres.
#Should be 3-3
Sim.add_extra_bond(mol1=[i,c],mol2=[i+Nchromosomes,c],typeb=liaison["2-2"][1]) #Should be 3-3
#Now create the interactions
simsoft = LSimu()
print liaison
for k,(bond_size,bond_type) in liaison.items():
simsoft.add_bond(typeb="harmonic",idbond=bond_type,K=350,R0=bond_size)
if bond_size == 0:
Sim.add_bond(typeb="harmonic",idbond=bond_type,K=0,R0=0)
elif bond_type == 6:
Sim.add_bond(typeb="harmonic",idbond=bond_type,K=350,R0=bond_size)
elif bond_type == 12:
Sim.add_bond(typeb="harmonic",idbond=bond_type,K=350,R0=bond_size)
else:
Sim.add_bond(typeb="fene",idbond=bond_type,K=30./(bond_size*bond_size),
R0=1.5*bond_size,epsilon=1,sigma=bond_size)
bead1,bead2=map(int,k.split("-"))
if bead1 == 5 or bead2 == 5:
simsoft.add_pair(typep="soft",idpair1=bead1,idpair2=bead2,A=0,rc=bond_size)
Sim.add_pair(typep="lj/cut",idpair1=bead1,idpair2=bead2,epsilon=0,sigma=bond_size,cutoff1=bond_size*1.15)
else:
simsoft.add_pair(typep="soft",idpair1=bead1,idpair2=bead2,A=5,rc=bond_size)
Sim.add_pair(typep="lj/cut",idpair1=bead1,idpair2=bead2,epsilon=1,sigma=bond_size,cutoff1=bond_size*1.15)
Sim.add_box(Box([-1*boxf*Radius,-1*boxf*Radius,-1*boxf*Radius],[quad*boxf*Radius,boxf*Radius,boxf*Radius]))
if not simple and oSim != None:
for i in range(Nchromosomes):
Sim.molecules[i].coords = oSim.molecules[i].coords
#Create the copy with shifted chromosomes
Sim.molecules[i+Nchromosomes].coords = oSim.molecules[i].coords + 1
return Sim,simsoft
#Now we have to create a template for the simulation.
Template="""
################################
#Template
#Must contain the variables
# typecell
# outtraj
# outfile
# radius
# interaction
# run_length
# samplingrate
# particle
# VARIABLES
variable tcell index $typecell # Define the variable from the tempalte
variable fname index ${tcell}conf2.txt # configuration initiale
# Initialization
#correspond to x=y=z=1
lattice fcc 4
units lj
boundary f f f
atom_style molecular
log log.txt
read_data ${fname}
neighbor 2.0 multi
comm_modify cutoff 30.0
include $softinteractions
group particle type 1 2 3 4
group normal type 1 3
group telo type 2
group ribo type 4
compute hic particle pair/local dist
compute hicp particle property/local patom1 patom2
dump init all dcd $samplingrate $outtraj.${tcell}.comp.dcd
###########################################################
#Definiton of nucleus and its interaction
#the telomere part is added when the nuceus has the right size
variable rad equal $radius
region mySphere sphere 0.0 0.0 0.0 v_rad side in
fix wall1 normal wall/region mySphere lj126 1 0.5 0.56
fix wall2 ribo wall/region mySphere lj126 1 1.62 1.81
fix wall telo wall/region mySphere lj93 4 2 6
#####################################################
# Equilibration (Langevin dynamics )
velocity particle create 1.0 1231
fix 1 particle nve/limit 0.0005
fix lang particle langevin 1.0 1.0 1.0 904297
run 30000
unfix 1
fix 1 particle nve/limit 0.005
run 30000
unfix 1
fix 1 particle nve/limit 0.05
run 30000
include $interaction
thermo_style custom step temp
thermo 10000
fix 1 particle nve/limit 0.05
timestep 0.005
run $run_length
write_data $outfile
"""
with open("tscript","w") as f:
f.writelines(Template)
Sim,simsoft = create_nucleus(dnuc=1)
{'2-2': [1, 6], '2-3': [1, 7], '2-4': [1.0, 8], '2-5': [0, 9], '1-1': [1, 1], '1-3': [1, 3], '1-2': [1, 2], '1-5': [0, 5], '1-4': [1.0, 4], '5-5': [0, 15], '4-4': [1, 13], '4-5': [0, 14], '3-5': [4.98, 12], '3-4': [1.0, 11], '3-3': [1, 10]}
#Then let's generate all the files needed to run the simulation:
REP="./"
cell="nucleus_yeast"
Sim.generate_xyz(REP+"/%sconf2.txt"%cell,Mass={"1":1,"2":1,"3":1,"4":1,"5":1})
Sim.generate_pdb(REP+"/%snoyau2.pdb"%cell,shift=1,
traduction={"1":"bead","2":"telo","3":"cent","4":"ribo","5":"spbb","6":"rcut","7":"scut"})
Sim.generate_interactions(REP+"/interactions",info_bond=["special_bonds fene"])
simsoft.generate_interactions(REP+"/softinteractions")
Sim.generate_script(REP+"/nucleus_init.txt",template_name="./tscript",outfile="final.xyz",
outtraj="dump_init",samplingrate=100,run_length=100000,
interaction="interactions",
softinteractions="softinteractions",typecell=cell,
radius=Radius)
Sim.cmd="mpirun -np 4 lammps"
r = Sim.run("nucleus_init.txt")
#cat log.txt
mpirun -np 4 lammps < nucleus_init.txt
#Now we have to create a template for the simulation.
Template="""
################################
#Template
#Must contain the variables
# typecell
# outtraj
# outfile
# radius
# interaction
# run_length
# samplingrate
# particle
# VARIABLES
variable tcell index $typecell # Define the variable from the tempalte
variable fname index ${tcell}conf2.txt # configuration initiale
# Initialization
#correspond to x=y=z=1
lattice fcc 4
units lj
boundary f f f
atom_style molecular
log log.txt
read_data ${fname}
neighbor 2.0 multi
comm_modify cutoff 30.0
include $softinteractions
group particle type 1 2 3 4
group normal type 1 3
group telo type 2
group ribo type 4
group spb id $spb
compute hic particle pair/local dist
compute hicp particle property/local patom1 patom2
dump init all dcd $samplingrate $outtraj.${tcell}.comp.dcd
###########################################################
#Definiton of nucleus and its interaction
#the telomere part is added when the nuceus has the right size
variable rad equal $radius
region mySphere sphere 0.0 0.0 0.0 v_rad side in
fix wall1 normal wall/region mySphere lj126 1 0.5 0.56
fix wall2 ribo wall/region mySphere lj126 1 1.62 1.81
fix wall telo wall/region mySphere lj93 4 2 6
#####################################################
# Equilibration (Langevin dynamics )
velocity particle create 1.0 1231
fix 1 particle nve/limit 0.0005
fix lang particle langevin 1.0 1.0 1.0 904297
run 30000
unfix 1
fix 1 particle nve/limit 0.005
run 30000
unfix 1
fix 1 particle nve/limit 0.05
run 200000
#Make the spb rotate
fix 10 spb move rotate 0.0 0.0 0.0 0.0 1.0 0.0 $period
#Increase the size of the spring between spb and centromere
#from Mt ot radius size
label loop
variable a loop 10
include $interaction
variable cd equal $($Mt + v_a *(1.1*$radius-$Mt)/10.0)
bond_coeff 12 harmonic 100 ${cd}
thermo_style custom step temp
thermo 1000
fix 1 particle nve/limit 0.05
timestep 0.005
run $run_length
next a
jump SELF loop
label break
variable a delete
write_data $outfile
"""
with open("tscript","w") as f:
f.writelines(Template)
#Take the result from the first simulation,
#Then duplicate the chromosomes
Sim.get_atoms_bonds_angles("final.xyz")
Sim2,simsoft2 = create_nucleus(simple=False,oSim=Sim,dnuc=1,boxf=1.3)
{'2-2': [1, 6], '2-3': [1, 7], '2-4': [1.0, 8], '2-5': [0, 9], '1-1': [1, 1], '1-3': [1, 3], '1-2': [1, 2], '1-5': [0, 5], '1-4': [1.0, 4], '5-5': [0, 15], '4-4': [1, 13], '4-5': [0, 14], '3-5': [4.98, 12], '3-4': [1.0, 11], '3-3': [1, 10]}
#Then let's generate all the files needed to run the simulation:
REP="./"
cell="rnucleus_yeast"
Sim2.generate_xyz(REP+"/%sconf2.txt"%cell,Mass={"1":1,"2":1,"3":1,"4":1,"5":1})
Sim2.generate_pdb(REP+"/%snoyau2.pdb"%cell,shift=1,traduction={"1":"bead","2":"telo","3":"cent","4":"ribo","5":"spbb","6":"rcut","7":"scut"}) # Not necellary to run the simulation but usefull for the analysis
Sim2.generate_interactions(REP+"/interactions",info_bond=["special_bonds fene"])
simsoft.generate_interactions(REP+"/softinteractions")
Sim2.generate_script(REP+"/nucleus_init.txt",template_name="./tscript",outfile="rfinal.xyz",
outtraj="dump_init",samplingrate=100,run_length=100000,period=2*0.005*1000000,
interaction="interactions",Mt = Mt,
softinteractions="softinteractions",typecell=cell,spb=len(Sim2.Atom),
radius=Radius+1)
Sim2.cmd ="mpirun -np 4 lammps"
Sim2.run("nucleus_init.txt")
#!cat log.txt
mpirun -np 4 lammps < nucleus_init.txt
--------------------------------------------------------------------------- CalledProcessError Traceback (most recent call last) <ipython-input-176-1bb0bb4bd811> in <module>() 1 Sim2.cmd ="mpirun -np 4 lammps" ----> 2 Sim2.run("nucleus_init.txt") 3 #!cat log.txt /home/jarbona/cpolymer/cpolymer/lsimu.pyc in run(self, script) 718 if self.script != None: 719 print "{0} < {1}".format(self.cmd,self.script) --> 720 output = subprocess.check_output("{0} < {1}".format(self.cmd,self.script), shell=True) 721 return output 722 else: /usr/lib/python2.7/subprocess.pyc in check_output(*popenargs, **kwargs) 571 if cmd is None: 572 cmd = popenargs[0] --> 573 raise CalledProcessError(retcode, cmd, output=output) 574 return output 575 CalledProcessError: Command 'mpirun -np 4 lammps < nucleus_init.txt' returned non-zero exit status 1
#Now we have to create a template for the simulation.
Template="""
################################
#Template
#Must contain the variables
# typecell
# outtraj
# outfile
# radius
# interaction
# run_length
# samplingrate
# particle
# VARIABLES
variable tcell index $typecell # Define the variable from the tempalte
variable fname index ${tcell}conf2.txt # configuration initiale
# Initialization
#correspond to x=y=z=1
lattice fcc 4
units lj
boundary f f f
atom_style molecular
log log.txt
read_data ${fname}
neighbor 2.0 multi
comm_modify cutoff 30.0
include $softinteractions
group particle type 1 2 3 4
group normal type 1 3
group telo type 2
group ribo type 4
group spb id $spb
compute hic particle pair/local dist
compute hicp particle property/local patom1 patom2
dump init all dcd $samplingrate $outtraj.${tcell}.comp.dcd
###########################################################
#Definiton of nucleus and its interaction
#the telomere part is added when the nuceus has the right size
variable rad equal $radius
region mySphere sphere 0.0 0.0 0.0 v_rad side in
fix wall1 normal wall/region mySphere lj126 1 0.5 0.56
fix wall2 ribo wall/region mySphere lj126 1 1.62 1.81
fix wall telo wall/region mySphere lj93 4 2 6
#####################################################
# Equilibration (Langevin dynamics )
include $interaction
variable mtl equal $(1.1*$radius)
bond_coeff 12 harmonic 100 ${mtl}
thermo_style custom step temp
thermo 10000
fix 1 particle nve/limit 0.05
timestep 0.005
run $run_length
write_data $outfile
"""
with open("tscript","w") as f:
f.writelines(Template)
a = Sim2.get_atoms_bonds_angles("rfinal.xyz")
#Sim2,simsoft2 = create_nucleus(simple=False,oSim=Sim,dnuc=1,boxf=1.3)
#Then let's generate all the files needed to run the simulation:
REP="./"
cell="equilibriate_yeast"
Sim2.generate_xyz(REP+"/%sconf2.txt"%cell,Mass={"1":1,"2":1,"3":1,"4":1,"5":1})
Sim2.generate_pdb(REP+"/%snoyau2.pdb"%cell,shift=1,traduction={"1":"bead","2":"telo","3":"cent","4":"ribo","5":"spbb","6":"rcut","7":"scut"}) # Not necellary to run the simulation but usefull for the analysis
Sim2.generate_interactions(REP+"/interactions",info_bond=["special_bonds fene"])
simsoft.generate_interactions(REP+"/softinteractions")
Sim2.generate_script(REP+"/nucleus_equi.txt",template_name="./tscript",outfile="refinal.xyz",
outtraj="dump_init",samplingrate=1000,run_length=1000000,
interaction="interactions",
softinteractions="softinteractions",typecell=cell,spb=len(Sim2.Atom),
radius=Radius+1)
Sim2.cmd ="mpirun -np 4 lammps"
i = Sim2.run("nucleus_equi.txt")
#Now we have to create a template for the simulation.
Template="""
################################
#Template
#Must contain the variables
# typecell
# outtraj
# outfile
# radius
# interaction
# run_length
# samplingrate
# particle
# VARIABLES
variable tcell index $typecell # Define the variable from the tempalte
variable fname index ${tcell}conf2.txt # configuration initiale
# Initialization
#correspond to x=y=z=1
lattice fcc 4
units lj
boundary f f f
atom_style molecular
log log.txt
read_data ${fname}
neighbor 2.0 multi
comm_modify cutoff 60
group particle type 1 2 3 4
group normal type 1 3
group telo type 2
group ribo type 4
group spb id $spb
compute hic particle pair/local dist
compute hicp particle property/local patom1 patom2
dump init all dcd $samplingrate $outtraj.${tcell}.comp.dcd
###########################################################
#Definiton of nucleus and its interaction
variable rad equal $radius
include $interaction
bond_coeff 12 harmonic 100 $Startm
#Make one of the spb move
fix 10 spb move linear $speed NULL NULL
#Make the two regions and the first one will move
label loop
variable a universe $rampf
region nuc1 sphere $a 0.0 0.0 v_rad side in
region nuc2 sphere 0.0 0.0 0.0 v_rad side in
region mySphere union 2 nuc1 nuc2
fix wall1 normal wall/region mySphere lj126 1 2 2.24
fix wall2 ribo wall/region mySphere lj126 1 2 2.24
fix wall telo wall/region mySphere lj93 4 3 8
variable ce equal $(350 - v_a*350/($stop))
variable cd equal $(1 + v_a *(2*$radius-2*$Mt-1)/($stop))
variable mts equal $($Startm + v_a *($Mt-$Startm)/($stop))
# at some point separate all the centromere from each other
if "$a < $stop" then "bond_coeff 6 harmonic ${ce} ${cd}"
if "$a == $stop" then "delete_bonds all bond 6 remove special"
# and set the spb bond to the right size
if "$a < $stop" then "bond_coeff 12 harmonic 100 ${mts} "
if "$a == $stop" then "bond_coeff 12 harmonic 100 $Mt "
#variable A0 equal $A
#variable pspb equal $($radius+ v_a)
#print $A
#
print "$a ${ce} ${cd}"
thermo_style custom step temp
thermo 1000
fix 1 particle nve/limit 0.05
timestep 0.005
run $run_length
unfix wall1
unfix wall2
unfix wall
region mySphere delete
region nuc1 delete
region nuc2 delete
#unfix 10
next a
jump SELF loop
label break
variable a delete
#variable Spos equal $dradius
#run 1000000
#write_data $outfile
"""
with open("tscript","w") as f:
f.writelines(Template)
a = Sim2.get_atoms_bonds_angles("refinal.xyz")
Sim2.box.tr[0] = 4*1.1*Radius
#Then let's generate all the files needed to run the simulation:
increase=0.5
stepi = 10000
REP="./"
cell="divide_yeast"
nramp = ["%.1f"%n for n in np.arange(increase,36.5+2,increase)]
ramp = " ".join(nramp)
Sim2.generate_xyz(REP+"/%sconf2.txt"%cell,Mass={"1":1,"2":1,"3":1,"4":1,"5":1})
Sim2.generate_pdb(REP+"/%snoyau2.pdb"%cell,shift=1,traduction={"1":"bead","2":"telo","3":"cent","4":"ribo","5":"spbb","6":"rcut","7":"scut"}) # Not necellary to run the simulation but usefull for the analysis
Sim2.generate_interactions(REP+"/interactions",info_bond=["special_bonds fene"])
simsoft.generate_interactions(REP+"/softinteractions")
Sim2.generate_script(REP+"/nucleus_divide.txt",template_name="./tscript",outfile="Ffinal.xyz",
outtraj="dump_init",samplingrate=100,run_length=stepi,
rampf=ramp,increase=increase,stop=nramp[10],
interaction="interactions",Startm=1.1*Radius,Mt=0.30*Radius,
softinteractions="softinteractions",typecell=cell,spb=len(Sim2.Atom),
speed = increase/(stepi*0.005),
radius=Radius+2.5)
#create trajectory of nucleus:
t = []
for i in range(13590):
t.append("2\nFix\na 0 0 0 \nb 0 0 0\n")
for i in range(13590,20604):
p =(i-13590) *(2*17.6)/(20604-13590)
t.append("2\nFix\na 0 0 0 \nb %i 0 0\n"%p)
with open("test.xyz","w") as f:
f.write("".join(t))