Tutorials: Configuration Relaxation¶
Import the necessary packages¶
[1]:
import parmed as pmd
import json
import sys
from sys import platform
import mdtraj as md
import simtk.openmm.app as mmapp
import mdtraj.reporters
import simtk.unit as u
import random
from openNucleome import OpenNucleome
from openNucleome.utils import coor_transformation, final_frame
import warnings
import numpy as np
Important parameters¶
We defined the important parameters in the next block. First, we set the transition probability between dP particles and P particles as 0.2, and the transition frequency as 4000. When creating the system, we set type 6, 7 as the dP and P particles, respectively, so we kept using 6, 7 here. In this example, we ran a simulation with total length of 100,000 steps, and output one configuration and the energy every 2000 steps.
[2]:
prob_P_dP = 0.2 # Transition probability from P to dP
prob_dP_P = 0.2 # Transition probability from dP to P
transition_freq = 4000
sampling_freq = 2000
dP_type = 6
P_type = 7
total_steps = 100000
Initialize the system¶
We first set up an example “model” of class “OpenNucleome” with the conserved temperature, damping coefficient, timestep, and the mass scale. In this folder, we also included the initial configuration “human.pdb” used for the simulation and created a system according to the initial configuration.
In this example, we freezed all the lamina beads, and would not consider the dynamics of the membrane, and that is why we set “False” (off) for membrane dynamics. Consequently, there was also no need to set the bond between lamina beads, so we set “None” for the variable “lam_bond”.
[3]:
model = OpenNucleome(1.0, 0.1, 0.005, 1.0) # 1.0: temperature (LJ reduced unit);
# 0.1: damping coefficient (LJ reduced unit);
# 0.005: timestep (LJ reduced unit);
# 1.0: mass_scale
PDB_file = "../configs_generation/init_config_pool/human_1.pdb" #The initial configuration
# Generate new elements and construct topology as well
# flag_membrane: True for including lamina dynamics, False for excluding lamina dynamics;
# lam_bond: A file contains the lamina bond when membrane_dynamics is on.
warnings.filterwarnings("ignore")
model.create_system(PDB_file, flag_membrane = False, lam_bond = None)
Add the force field¶
In this example, we loaded default settings for the interactions between chromosomes and chromosomes, and chromosomes and nuclear landmarks. All the types of potential can be found in “Section: Energy function of the whole nucleus model” in Supporting Information. According to the order of added potential, the index of speckle-speckle potential is 6, and this index is needed when transiting the speckle particle between type 6 and type 7.
Because during the simulations, we would transit the dP speckle and P speckle particles, here, we logged the start index and end index of speckle particle, and computed the number of speckle particles.
[4]:
# Add the default force field
model.load_default_settings()
index_spec_spec_potential = 6 #The 6-th potential is speckle-speckle interaction in the default settings
start_spec_index = model.N_chr_nuc+1
end_spec_index = model.N_chr_nuc_spec+1
N_spec = end_spec_index-start_spec_index #Currently, the number of speckle is 1600.
Perform the simulation¶
We first created the simulation with a specific Platform, and here, we used “CUDA” but users can also use “CPU”, “Reference”, and “OpenCL” according to their hardware. Before the simulation, we minimized the energy to make the system much more reasonable and stable. After randomly setting velocity, we started our simulation with a total length of 100,000 steps and output the configuration and energy every 2000 steps, and change the speckle types every 4000 steps as we mentioned in the previous blocks.
[5]:
#model.save_system("model_before_simulation_0.xml")
simulation = model.create_simulation(platform_type = "CUDA") # Users can also use CPU, Reference, OpenCL.
simulation.context.setPositions(model.chr_positions)
simulation.minimizeEnergy()
simulation.reporters.append(mdtraj.reporters.DCDReporter('results/step_100000.dcd', sampling_freq))
def setVelocity(context):
sigma = u.sqrt(1.0*u.kilojoule_per_mole / model.chr_system.getParticleMass(1))
velocs = u.Quantity(1.0 * np.random.normal(size=(model.chr_system.getNumParticles(), 3)),
u.meter) * (sigma / u.meter)
context.setVelocities(velocs)
setVelocity(simulation.context)
simulation.reporters.append(mmapp.statedatareporter.StateDataReporter(sys.stdout, sampling_freq, step=True,
potentialEnergy=True, kineticEnergy=True, temperature=True, progress=True,
remainingTime=True, separator='\t', totalSteps = total_steps))
for i in range(total_steps//transition_freq):
simulation.step(transition_freq)
# Change the type of speckles every 4000 steps, non-equilibrium scheme.
# Do the chemical modification, and change the spec-spec potential on the fly.
for j in np.random.randint(start_spec_index-1, end_spec_index-1, N_spec):
if model.compart_type[j] == dP_type-1:
model.compart_type[j] = P_type-1 if random.random() < prob_dP_P else dP_type-1
else:
model.compart_type[j] = dP_type-1 if random.random() < prob_P_dP else P_type-1
# Update the context after changing the type of speckles.
for m in range(model.chr_system.getNumParticles()):
model.chr_system.getForce(index_spec_spec_potential).setParticleParameters(m, [model.compart_type[m]])
model.chr_system.getForce(index_spec_spec_potential).updateParametersInContext(simulation.context)
# Keep the final result of bead types in case constructing the configuration for the continuous simulation.
np.savetxt('results/compt_final_frame.txt', (np.array(model.compart_type)+1).reshape((-1,1)), fmt='%d')
#"Progress (%)" "Step" "Potential Energy (kJ/mole)" "Kinetic Energy (kJ/mole)" "Temperature (K)" "Time Remaining"
2.0% 2000 305517.57944424567 239729.9228943532 307.3442102983061 --
4.0% 4000 207185.90673694137 151506.21446981165 194.23757067669288 14:35
6.0% 6000 158288.3406172468 121372.09834880882 155.60432034887066 14:01
8.0% 8000 132847.48986336985 107250.5440399466 137.49987220638914 13:29
10.0% 10000 121596.96969540659 100929.74502448138 129.39633236285206 13:05
12.0% 12000 113962.34265134309 97423.7777640623 124.90152951982073 12:41
14.0% 14000 112257.73056958467 96318.60754455603 123.4846531271768 12:20
16.0% 16000 109864.41145441416 94804.20224058055 121.54311951884611 11:59
18.0% 18000 109899.56832139855 94648.97944137751 121.34411712454603 11:39
20.0% 20000 108446.11521933178 94627.91730774056 121.31711455118477 11:19
22.0% 22000 109445.79017165932 94263.81045584117 120.85031369666426 11:00
24.0% 24000 108907.84712871362 93942.05751382257 120.43781240069778 10:40
26.0% 26000 109715.24762966722 94501.93586251422 121.15560084723053 10:22
28.0% 28000 109012.5285738545 94112.84046071759 120.65676358256171 10:03
30.0% 30000 109207.79834143931 93955.38224187843 120.4548952828525 9:45
32.0% 32000 108795.9361688504 93830.04515171243 120.29420767017125 9:27
34.0% 34000 108592.09765447844 94863.32459109789 121.618916949185 9:09
36.0% 36000 109532.47760620397 93814.96253464038 120.27487109767635 8:51
38.0% 38000 108626.15925753233 95301.87573854289 122.1811586354108 8:34
40.0% 40000 108300.33539516928 94183.20974602466 120.74698007352303 8:17
42.0% 42000 109516.66493095538 93640.43991764485 120.05112549574196 8:00
44.0% 44000 108996.12333174223 93769.1005199753 120.21607399588052 7:42
46.0% 46000 109280.6361911603 93836.14837655456 120.30203227044113 7:25
48.0% 48000 108364.16466056774 93853.77001240787 120.32462397568223 7:08
50.0% 50000 109210.26718781143 94154.07247884736 120.70962482699086 6:51
52.0% 52000 108715.65516172761 94486.79727268191 121.1361925152173 6:35
54.0% 54000 108919.86091815002 94448.96149973289 121.08768540514258 6:18
56.0% 56000 108705.23455818793 94537.86794517544 121.20166734337985 6:01
58.0% 58000 109063.11636587992 93659.14819804284 120.07511032665309 5:44
60.0% 60000 108135.8393093122 94429.40193950695 121.06260919638713 5:28
62.0% 62000 109449.988285381 94249.21790594471 120.83160541164911 5:11
64.0% 64000 108265.71772550428 94931.98987010222 121.7069488297934 4:54
66.0% 66000 108873.47916040457 93522.16191102509 119.89948793585242 4:38
68.0% 68000 108727.68894074738 93940.03245892211 120.435216192049 4:21
70.0% 70000 109450.27407430168 94213.81964361423 120.78622329643153 4:05
72.0% 72000 108723.09269815465 94072.10699834587 120.60454150834501 3:48
74.0% 74000 108844.80312130705 94326.39749392454 120.9305530074787 3:32
76.0% 76000 108184.15437296525 93909.8220192326 120.39648509157962 3:15
78.0% 78000 108687.59792681698 93698.60431995375 120.12569479472339 2:59
80.0% 80000 107990.74859291193 94525.97634042453 121.18642178776784 2:42
82.0% 82000 108185.52120780546 94680.5950308404 121.38464968827138 2:26
84.0% 84000 108176.67632892226 94440.20228596928 121.07645571130833 2:10
86.0% 86000 108495.17923021212 94492.26184807981 121.14319833906276 1:53
88.0% 88000 107678.69027082807 94359.94717841972 120.97356516533716 1:37
90.0% 90000 107893.99383572253 94416.27027058577 121.04577382445795 1:21
92.0% 92000 106406.35815377338 93848.41319502125 120.31775630230409 1:04
94.0% 94000 108033.8972673082 93972.11037766634 120.47634148206463 0:48
96.0% 96000 107101.25632630658 93948.42374930475 120.4459741920721 0:32
98.0% 98000 108023.6646569059 94298.59749964751 120.8949122030926 0:16
100.0% 100000 106862.39324199257 93993.76087071122 120.5040983578318 0:00
Post-process the output files¶
We converted the trajectory with the default OpenMM unit to the nex trajectory with the reduced unit, and saved the final configuration according to the last frame.
[6]:
coor_transformation('results/step_100000.dcd', 'results/reduced_traj.dcd')
final_frame('results/step_100000.dcd', 'results/compt_final_frame.txt', 'results/final_config.pdb')