Tutorial

“What I cannot create, I do not understand.” — Feynman

DNA from sequence

This section allows you to quickly simulate a piece of DNA.

First import the open3SPN2 module.

import open3SPN2

If you obtain a ModuleNotFoundError error then you may need to check that open3SPN2 is in the installation path as detailed in Installation.

After importing the module, follow the next steps to create a DNA system and add the 3SPN2 forces. Make sure you have installed X3DNA before this step.

# Initialize the DNA from a sequence.
# DNA type can be changed to 'A' or 'B'

seq='ATACAAAGGTGCGAGGTTTCTATGCTCCCACG'
dna=open3SPN2.DNA.fromSequence(seq,dna_type='B_curved')

# Compute the topology for the DNA structure.
# Since the dna was generated from the sequence using X3DNA,
# it is not necesary to recompute the geometry.

dna.computeTopology(template_from_X3DNA=False)

# Create the system.
# To set periodic boundary conditions (periodicBox=[50,50,50]).
# The periodic box size is in nanometers.
dna.periodic=False
s=open3SPN2.System(dna, periodicBox=None)

#Add 3SPN2 forces
s.add3SPN2forces(verbose=True)

Bond
Angle
Stacking
Dihedral
BasePair
CrossStacking
Exclusion
Electrostatics

To set up the simulation you will need the openmm package as detailed in Installation.

import simtk.openmm
import simtk.openmm.app
import simtk.unit
import sys
import numpy as np

#Initialize Molecular Dynamics simulations
s.initializeMD(temperature=300 * simtk.unit.kelvin,platform_name='OpenCL')
simulation=s.simulation

#Set initial positions
simulation.context.setPositions(s.coord.getPositions())

energy_unit=simtk.openmm.unit.kilojoule_per_mole
#Total energy
state = simulation.context.getState(getEnergy=True)
energy = state.getPotentialEnergy().value_in_unit(energy_unit)
print('TotalEnergy',round(energy,6),energy_unit.get_symbol())

#Detailed energy
energies = {}
for force_name, force in s.forces.items():
    group=force.getForceGroup()
    state = simulation.context.getState(getEnergy=True, groups=2**group)
    energies[force_name] =state.getPotentialEnergy().value_in_unit(energy_unit)

for force_name in s.forces.keys():
    print(force_name, round(energies[force_name],6),energy_unit.get_symbol())

TotalEnergy -2046.579102 kJ/mol
Bond 2.1e-05 kJ/mol
Angle 0.001745 kJ/mol
Stacking -596.408569 kJ/mol
Dihedral -839.999756 kJ/mol
BasePair -518.252075 kJ/mol
CrossStacking -135.335464 kJ/mol
Exclusion 0.11901 kJ/mol
Electrostatics 43.296059 kJ/mol

Please make sure that the energies obtained coincide with the energies shown here. Also you can check the energy obtained using other platforms by changing OpenCL to Reference, CUDA or CPU.

#Add simulation reporters
dcd_reporter=simtk.openmm.app.DCDReporter(f'output.dcd', 1000)
energy_reporter=simtk.openmm.app.StateDataReporter(sys.stdout, 1000, step=True,time=True,
                                                   potentialEnergy=True, temperature=True)
simulation.reporters.append(dcd_reporter)
simulation.reporters.append(energy_reporter)

#Run simulation
simulation.step(10000)

#"Step","Time (ps)","Potential Energy (kJ/mole)","Temperature (K)"
1000,2.0000000000000013,-1651.121826171875,304.6066812070446
2000,3.999999999999781,-1646.61328125,309.5945230237376
3000,5.999999999999561,-1661.6788330078125,318.46432160647703
4000,7.999999999999341,-1676.956298828125,268.04874840144447
5000,10.000000000000009,-1629.8892822265625,271.8654648104738
6000,12.000000000000677,-1622.474853515625,312.1083112301662
7000,14.000000000001345,-1704.033203125,283.5259033832464
8000,16.00000000000201,-1608.751708984375,281.82371603990293
9000,18.000000000000902,-1623.486572265625,283.86225823944585
10000,19.999999999999794,-1671.9105224609375,296.18167366285144

The forces and system in open3SPN can be treated as forces and systems from openmm. Please refer to openmm documentation to learn more about running the simulations or adding forces.

Protein DNA system

You can find this example on the examples/Protein_DNA folder. For this example you need to download the structure of the Lambda repressor-operator complex 1mlb.pdb. You will also need to have installed the openAWSEM library.

# If you want to specify the package address
# you can add them to the PYTHONPATH environment variable.
# Also you can add them on the run time uncommenting the lines below
# import sys
# open3SPN2_HOME = '/Users/weilu/open3spn2/'
# openAWSEM_HOME = '/Users/weilu/openmmawsem/'
# sys.path.insert(0,open3SPN2_HOME)
# sys.path.insert(0,openAWSEM_HOME)

#Import openAWSEM, open3SPN2 and other libraries
import open3SPN2
import ffAWSEM

import pandas
import numpy as np
import simtk.openmm

from functools import partial
import sys

#Fix the system (adds missing atoms)
fix=open3SPN2.fixPDB("1lmb.pdb")

#Create a table containing both the proteins and the DNA
complex_table=open3SPN2.pdb2table(fix)

# Create a single memory file
ffAWSEM.create_single_memory(fix)

#Generate a coarse-grained model of the DNA molecules
dna_atoms=open3SPN2.DNA.CoarseGrain(complex_table)

#Generate a coarse-grained model of the Protein molecules
protein_atoms=ffAWSEM.Protein.CoarseGrain(complex_table)

#Merge the models
Coarse=pandas.concat([protein_atoms,dna_atoms],sort=False)
Coarse.index=range(len(Coarse))
Coarse['serial']=list(Coarse.index)

#Save the protein_sequence
ffAWSEM.save_protein_sequence(Coarse,sequence_file='protein.seq')

# Create a merged PDB
ffAWSEM.writePDB(Coarse,'clean.pdb')

#Create the merged system

pdb=simtk.openmm.app.PDBFile('clean.pdb')
top=pdb.topology
coord=pdb.positions
forcefield=simtk.openmm.app.ForceField(ffAWSEM.xml,open3SPN2.xml)
s=forcefield.createSystem(top)

dna=open3SPN2.DNA.fromCoarsePDB('clean.pdb')
with open('protein.seq') as ps:
    protein_sequence_one=ps.readlines()[0]
protein=ffAWSEM.Protein.fromCoarsePDB('clean.pdb',sequence=protein_sequence_one)
dna.periodic=False
protein.periodic=False

#Create the DNA and Protein Objects
dna=open3SPN2.DNA.fromCoarsePDB('clean.pdb')
with open('protein.seq') as ps:
    protein_sequence_one=ps.readlines()[0]
protein=ffAWSEM.Protein.fromCoarsePDB('clean.pdb',sequence=protein_sequence_one)
dna.periodic=False
protein.periodic=False

#Copy the AWSEM parameter files
ffAWSEM.copy_parameter_files()

#Clear Forces from the system (optional)
keepCMMotionRemover=True
j=0
for i, f in enumerate(s.getForces()):
    if keepCMMotionRemover and i == 0 and f.__class__ == simtk.openmm.CMMotionRemover:
        # print('Kept ', f.__class__)
        j += 1
        continue
    else:
        # print('Removed ', f.__class__)
        s.removeForce(j)
if keepCMMotionRemover == False:
    assert len(s.getForces()) == 0, 'Not all the forces were removed'
else:
    assert len(s.getForces()) <= 1, 'Not all the forces were removed'

#Initialize the force dictionary
forces={}
for i in range(s.getNumForces()):
    force = s.getForce(i)
    force_name="CMMotionRemover"

#Add 3SPN2 forces
for force_name in open3SPN2.forces:
    print(force_name)
    force = open3SPN2.forces[force_name](dna)
    if force_name in ['BasePair','CrossStacking']:
        force.addForce(s)
    else:
        s.addForce(force)
    forces.update({force_name:force})

#Add AWSEM forces
openAWSEMforces = dict(Connectivity=ffAWSEM.functionTerms.basicTerms.con_term,
                       Chain=ffAWSEM.functionTerms.basicTerms.chain_term,
                       Chi=ffAWSEM.functionTerms.basicTerms.chi_term,
                       Excl=ffAWSEM.functionTerms.basicTerms.excl_term,
                       rama=ffAWSEM.functionTerms.basicTerms.rama_term,
                       rama_pro=ffAWSEM.functionTerms.basicTerms.rama_proline_term,
                       contact=ffAWSEM.functionTerms.contactTerms.contact_term,
                       frag  = partial(ffAWSEM.functionTerms.templateTerms.fragment_memory_term,
                                       frag_file_list_file="./single_frags.mem",
                                       npy_frag_table="./single_frags.npy",
                                       UseSavedFragTable=False,
                                       k_fm=0.04184/3),
                       beta1 = ffAWSEM.functionTerms.hydrogenBondTerms.beta_term_1,
                       beta2 = ffAWSEM.functionTerms.hydrogenBondTerms.beta_term_2,
                       beta3 = ffAWSEM.functionTerms.hydrogenBondTerms.beta_term_3,
                       pap1 = ffAWSEM.functionTerms.hydrogenBondTerms.pap_term_1,
                       pap2 = ffAWSEM.functionTerms.hydrogenBondTerms.pap_term_2,
                      )
protein.setup_virtual_sites(s)

#Add DNA-protein interaction forces
for force_name in open3SPN2.protein_dna_forces:
    print(force_name)
    force = open3SPN2.protein_dna_forces[force_name](dna,protein)
    s.addForce(force)
    forces.update({force_name: force})

#Fix exclussions
for force_name in openAWSEMforces:
    print(force_name)
    if force_name in ['contact']:
        force = openAWSEMforces[force_name](protein, withExclusion=False,periodic=False)
        print(force.getNumExclusions())
        open3SPN2.addNonBondedExclusions(dna,force)
        print(force.getNumExclusions())
    elif force_name in ['Excl']:
        force = openAWSEMforces[force_name](protein)
        print(force.getNumExclusions())
        open3SPN2.addNonBondedExclusions(dna,force)
        print(force.getNumExclusions())
    else:
        force = openAWSEMforces[force_name](protein)
    s.addForce(force)
    forces.update({force_name: force})

Bond
Angle
Stacking
Dihedral
BasePair
CrossStacking
Exclusion
Electrostatics
ExclusionProteinDNA
ElectrostaticsProteinDNA
Connectivity
Chain
Chi
Excl
1205
1844
rama
rama_pro
contact
Number of atom:  1171 Number of residue:  179
Contact cutoff  1.0 nm
NonbondedMethod:  1
0
639
frag
Loading Fragment files(Gro files)
Saving fragment table as npy file to speed up future calculation.
All gro files information have been stored in the ./single_frags.npy.
You might want to set the 'UseSavedFragTable'=True to speed up the loading next time.
But be sure to remove the .npy file if you modify the .mem file. otherwise it will keep using the old frag memeory.
beta1
beta_1 term ON
beta2
beta_2 term ON
beta3
beta_3 term ON
pap1
pap_1 term ON
No ssweight given, assume all zero
pap2
pap_2 term ON
No ssweight given, assume all zero

# Set up the simulation
temperature=300 * simtk.openmm.unit.kelvin
platform_name='OpenCL' #'Reference','CPU','CUDA', 'OpenCL'

integrator = simtk.openmm.LangevinIntegrator(temperature, 1 / simtk.openmm.unit.picosecond, 2 * simtk.openmm.unit.femtoseconds)
platform = simtk.openmm.Platform.getPlatformByName(platform_name)
simulation = simtk.openmm.app.Simulation(top,s, integrator, platform)
simulation.context.setPositions(coord)
energy_unit=simtk.openmm.unit.kilojoule_per_mole
state = simulation.context.getState(getEnergy=True)
energy = state.getPotentialEnergy().value_in_unit(energy_unit)
print(energy)

-899.144763339983

#Obtain total energy

energy_unit=simtk.openmm.unit.kilojoule_per_mole
state = simulation.context.getState(getEnergy=True)
energy = state.getPotentialEnergy().value_in_unit(energy_unit)
print('TotalEnergy',round(energy,6),energy_unit.get_symbol())

#Obtain detailed energy

energies = {}
for force_name, force in forces.items():
    group=force.getForceGroup()
    state = simulation.context.getState(getEnergy=True, groups=2**group)
    energies[force_name] =state.getPotentialEnergy().value_in_unit(energy_unit)

for force_name in forces.keys():
    print(force_name, round(energies[force_name],6),energy_unit.get_symbol())

TotalEnergy -899.147583 kJ/mol
Bond 327.558105 kJ/mol
Angle 973.859009 kJ/mol
Stacking 203.565979 kJ/mol
Dihedral -385.277161 kJ/mol
BasePair -284.232208 kJ/mol
CrossStacking -47.586143 kJ/mol
Exclusion 23.991552 kJ/mol
Electrostatics 23.268274 kJ/mol
ExclusionProteinDNA 296.033478 kJ/mol
ElectrostaticsProteinDNA -10.459805 kJ/mol
Connectivity 1899.296875 kJ/mol
Chain 1899.296875 kJ/mol
Chi 1899.296875 kJ/mol
Excl 1899.296875 kJ/mol
rama -1363.522705 kJ/mol
rama_pro -1363.522705 kJ/mol
contact -1041.547607 kJ/mol
frag -1213.29834 kJ/mol
beta1 -300.796692 kJ/mol
beta2 -300.796692 kJ/mol
beta3 -300.796692 kJ/mol
pap1 0.0 kJ/mol
pap2 0.0 kJ/mol

#Add simulation reporters
dcd_reporter=simtk.openmm.app.DCDReporter(f'output.dcd', 10000)
energy_reporter=simtk.openmm.app.StateDataReporter(sys.stdout, 10000, step=True,time=True,
                                                   potentialEnergy=True, temperature=True)
simulation.reporters.append(dcd_reporter)
simulation.reporters.append(energy_reporter)

#Run simulation
simulation.minimizeEnergy()
simulation.context.setVelocitiesToTemperature(temperature)
simulation.step(100000)

#"Step","Time (ps)","Potential Energy (kJ/mole)","Temperature (K)"
1000,2.0000000000000013,-3507.7216796875,291.6758662031427
2000,3.999999999999781,-3233.395751953125,299.73657773689456
3000,5.999999999999561,-3472.61083984375,314.04882500167656
4000,7.999999999999341,-3199.33251953125,309.8776608226663
5000,10.000000000000009,-3275.126220703125,321.27214200690065
6000,12.000000000000677,-3167.643798828125,295.1803751740272
7000,14.000000000001345,-3187.2998046875,310.9197062404284
8000,16.00000000000201,-3296.966064453125,305.38198987465074
9000,18.000000000000902,-3182.54443359375,316.07187422798313
10000,19.999999999999794,-3229.941650390625,309.6002450725328

#Get the detailed energy after the simulation
energies = {}
for force_name, force in forces.items():
    group=force.getForceGroup()
    state = simulation.context.getState(getEnergy=True, groups=2**group)
    energies[force_name] =state.getPotentialEnergy().value_in_unit(energy_unit)

for force_name in forces.keys():
    print(force_name, round(energies[force_name],6),energy_unit.get_symbol())

Bond 102.787193 kJ/mol
Angle 169.343231 kJ/mol
Stacking -411.438232 kJ/mol
Dihedral -452.744202 kJ/mol
BasePair -242.723328 kJ/mol
CrossStacking -48.275833 kJ/mol
Exclusion 1.823019 kJ/mol
Electrostatics 23.37781 kJ/mol
ExclusionProteinDNA -6.992864 kJ/mol
ElectrostaticsProteinDNA -10.225787 kJ/mol
Connectivity 1487.126465 kJ/mol
Chain 1487.126587 kJ/mol
Chi 1487.126465 kJ/mol
Excl 1487.126587 kJ/mol
rama -1456.633789 kJ/mol
rama_pro -1456.633789 kJ/mol
contact -1326.107178 kJ/mol
frag -922.991699 kJ/mol
beta1 -136.266449 kJ/mol
beta2 -136.266449 kJ/mol
beta3 -136.266449 kJ/mol
pap1 0.0 kJ/mol
pap2 0.0 kJ/mol