First w need to define all the physical values of the simulation:
We are going to create a simple nucleus with 16 chromosome of size 200 kilo base. The physical values that we need to define are the compaction that we will take to 50 base / nm. So every chromosome will we of length 4 micron. We next need to define the width of the fiber. We can choose 30 nm. This will give 4000/30 = 133 beads, where each bead will represent 1500 base. Finally we can also define a persistence length.
Then to scale the simulation in Lennard-Jones units.
To scale all the simulation, the size of the most common beads is set to one and this is the reference. for example if we choose beads of 60 nm in Diameter, an nucleus of Radius one micron, the diameter of the beads in the simulation is one, and the Radius = 1000/60 = 16.6
import numpy as np
from cpolymer.polymer import Polymer
from cpolymer.lsimu import LSimu
from cpolymer.constrain import Box,Sphere
diameter =0.030
Radius = 5/diameter
Nchromosomes = 20
Schromosome = int(3000000/5/20)
nucleus = Sphere(position=[0,0,0],radius=Radius)
Sim = LSimu()
bead_type = 1 # ALl the beads are going to be the same type
for X in range(Nchromosomes):
bead_size=1
bond_type=1
Sim.add(Polymer(N=Schromosome,type_bead=1,
liaison={"{0}-{0}".format(bead_type):[bead_size,bond_type]},
gconstrain=[nucleus]))
Sim.add_bond(typeb="harmonic",idbond=bond_type,K=350,R0=bead_size)
Sim.add_pair(typep="lj/cut",idpair1=bead_type,idpair2=bead_type,epsilon=1,sigma=bead_size,cutoff1=1.15)
Sim.add_box(Box([-1.1*Radius,-1.1*Radius,-1.1*Radius],[1.1*Radius,1.1*Radius,1.1*Radius]))
3000000/5/20
30000.0
#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
include $interaction
group particle type 1
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 particle wall/region mySphere lj126 1 1 1.12
#####################################################
# 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 20000
unfix 1
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)
#Then let's generate all the files needed to run the simulation:
!mkdir human
REP="./human"
cell="simple_yeast"
Sim.generate_xyz(REP+"/%sconf2.txt"%cell,Mass={str(bead_type):1})
Sim.generate_pdb(REP+"/%snoyau2.pdb"%cell,shift=1) # Not necellary to run the simulation but usefull for the analysis
Sim.generate_interactions(REP+"/interactions")
Sim.generate_script(REP+"/nucleus_init.txt",template_name="./tscript",outfile="final.xyz",
outtraj="dump_init",samplingrate=1000,run_length=20000,
interaction="interactions",typecell=cell,
radius=Radius)
mkdir: cannot create directory ‘human’: File exists
r = Sim.run("nucleus_init.txt")
lammps < nucleus_init.txt
#Requier chemview to be installed:
from chemview import MolecularViewer, RepresentationViewer
from chemview import enable_notebook
enable_notebook()
#r = lambda: random.randint(0,255)
#print ['0x%02X%02X%02X' % (r(),r(),r()) for i in range(16)]
color = [0x2CDC48,0x40C5E2,0xDA88A6,0x1C9AB8,0x7AEAA7,0xA66EEB,
0x840CD7,0x439505,0x9DA7D6,0x8C7D9B,0x130D08,0xCFEDC9,
0x494FFC,0xA380A8,0xF55109,0xD8F950]
def draw(Atoms,lengts,rv):
start=0
for i,n1 in enumerate(lengths):
import random
#print color
#print np.array(Atoms,dtype=float)[start:start+n1,2:5].shape
coordinates = np.array(Atoms,dtype=float)[start:start+n1,2:5]
rv.add_representation('smoothline',{"coordinates":coordinates,
'color':color[i],
"resolution":10})
start += n1
rv = RepresentationViewer()
Atoms,Bonds,Angles,lengths = Sim.get_atoms_bonds_angles("final.xyz")
draw(Atoms,lengths,rv)
rv