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