In [1]:

    

import sys
import warnings
import simtk.openmm as mm
from simtk.openmm import app
from simtk import unit


    

In [4]:

    

from IPython.display import Image
Image('04-frozen-nanotube-water.png')


    

    
Out[4]:

In [2]:

    

# Parameters from Walther, J. H., et al. "Molecular dynamics simulations of carbon
# nanotubes in water." Center for turbulence research: proceedings of the summer
# program, (2000): 5-20.
CARBON_LJ_SIGMA = 3.191 * unit.angstrom
CARBON_LJ_EPSILON = 0.4396 * unit.kilojoules_per_mole

with warnings.catch_warnings():
    warnings.simplefilter('ignore')
    pdb = app.PDBFile('04-nanotube-input.pdb')

# Use modeller to remove the non-water atoms from the PDB
modeller = app.Modeller(pdb.topology, pdb.positions)
modeller.delete([r for r in modeller.topology.residues() if r.name != 'HOH'])
# and now make a system for just the water using tip3p
system = app.ForceField('tip3p.xml').createSystem(modeller.topology, nonbondedMethod=app.PME, nonbondedCutoff=1*unit.nanometer)
initialPositions = modeller.positions

# get the nonbonded force from the newly created system
nbforce = [system.getForce(i) for i in range(system.getNumForces()) if isinstance(system.getForce(i), mm.NonbondedForce)][0]

# add the carbon nanotube back to the topology as a new residue named "CNT"
topology = modeller.topology
cntReside = topology.addResidue('CNT', topology.addChain())

# now iterate though all the carbon atoms from the original pdb topology
for i, atom in enumerate(pdb.topology.atoms()):
    if atom.element.symbol == 'C':
        # and add them back to the system with zero mass. This freezes
        # the carbon atoms in place. They'll just stand still in the
        # simulation.
        index = system.addParticle(0)
        # we'll give them no charge, but a little LJ 
        nbforce.addParticle(0.0, CARBON_LJ_SIGMA, CARBON_LJ_EPSILON)
        initialPositions.append(pdb.positions[i])
        topology.addAtom('C', atom.element, cntReside)


    

In [3]:

    

totalSteps = 100000    # 200 picoseconds
reportInterval = 1000  #   2 picoseconds

# translate the center of the box (looks better in VMD)
for i in range(len(initialPositions)):
    initialPositions[i] += topology.getUnitCellDimensions() / 2

# set up our simulation
simulation = app.Simulation(topology, system, mm.LangevinIntegrator(300*unit.kelvin, 1/unit.picosecond, 2*unit.femtosecond))
simulation.context.setPositions(initialPositions)
simulation.context.setVelocitiesToTemperature(300*unit.kelvin)

# write the first frame of the simulation to pdb
with open('nanotube-simulation.pdb', 'w') as f:
    app.PDBFile.writeFile(simulation.topology, simulation.context.getState(getPositions=True).getPositions(), f)
# and then attach reporters
simulation.reporters.append(app.DCDReporter('nanotube-simulation.dcd', reportInterval))
simulation.reporters.append(app.StateDataReporter(sys.stdout, reportInterval, step=True, time=True,
    remainingTime=True, separator='\t', temperature=True, speed=True, totalSteps=totalSteps))

# and roll.
simulation.step(totalSteps)


    

    




#"Step"	"Time (ps)"	"Temperature (K)"	"Speed (ns/day)"	"Time Remaining"
1000	2.0	487.987494077	0	--
2000	4.0	357.420745452	257	1:05
3000	6.0	308.05365382	257	1:05
4000	8.0	319.500745796	257	1:04
5000	10.0	293.894438389	256	1:04
6000	12.0	301.080799064	255	1:03
7000	14.0	310.948865244	255	1:03
8000	16.0	307.765403153	254	1:02
9000	18.0	302.650415263	255	1:01
10000	20.0	304.805679508	255	1:01
11000	22.0	294.838490696	255	1:00
12000	24.0	310.374418597	255	0:59
13000	26.0	300.802588048	255	0:58
14000	28.0	296.700354312	255	0:58
15000	30.0	295.307265853	255	0:57
16000	32.0	296.119434238	255	0:56
17000	34.0	293.489350443	255	0:56
18000	36.0	282.946455886	255	0:55
19000	38.0	296.076877543	255	0:54
20000	40.0	304.831627347	255	0:54
21000	42.0	297.959218194	255	0:53
22000	44.0	321.31210389	256	0:52
23000	46.0	312.158418107	256	0:52
24000	48.0	302.077547245	256	0:51
25000	50.0	305.436995153	256	0:50
26000	52.0	297.323567892	256	0:50
27000	54.0	302.655175566	256	0:49
28000	56.0	308.377683916	256	0:48
29000	58.0	303.088325759	256	0:48
30000	60.0	284.611929039	256	0:47
31000	62.0	288.297989049	256	0:46
32000	64.0	283.478706278	256	0:45
33000	66.0	295.130225102	256	0:45
34000	68.0	286.653824879	256	0:44
35000	70.0	309.918303103	256	0:43
36000	72.0	291.744661914	256	0:43
37000	74.0	312.068939559	255	0:42
38000	76.0	290.574163172	255	0:41
39000	78.0	298.103499821	255	0:41
40000	80.0	305.290808591	255	0:40
41000	82.0	300.688519491	256	0:39
42000	84.0	286.713939475	256	0:39
43000	86.0	324.543581669	256	0:38
44000	88.0	322.00452475	256	0:37
45000	90.0	303.820493159	255	0:37
46000	92.0	292.427581169	255	0:36
47000	94.0	292.914920944	255	0:35
48000	96.0	299.377402469	255	0:35
49000	98.0	296.743999297	256	0:34
50000	99.9999999999	303.203641318	256	0:33
51000	102.0	297.643546736	256	0:33
52000	104.0	293.718019725	256	0:32
53000	106.0	313.950795878	255	0:31
54000	108.0	299.674532108	255	0:31
55000	110.0	306.275319921	255	0:30
56000	112.0	293.709466454	255	0:29
57000	114.0	314.668181283	255	0:29
58000	116.0	291.348303238	255	0:28
59000	118.0	306.627825662	255	0:27
60000	120.0	305.319693548	255	0:27
61000	122.0	311.555212044	255	0:26
62000	124.0	308.499301785	255	0:25
63000	126.0	300.339848451	255	0:25
64000	128.0	317.796954398	255	0:24
65000	130.0	299.349716673	255	0:23
66000	132.0	308.216770798	255	0:23
67000	134.0	296.627442085	255	0:22
68000	136.0	303.083952123	255	0:21
69000	138.0	309.351169251	255	0:21
70000	140.0	289.371084118	255	0:20
71000	142.0	302.075808804	255	0:19
72000	144.0	295.128431173	255	0:18
73000	146.0	306.287846024	255	0:18
74000	148.0	309.278884046	255	0:17
75000	150.0	307.435083718	255	0:16
76000	152.0	302.588928022	255	0:16
77000	154.0	284.26431185	255	0:15
78000	156.0	286.600094832	255	0:14
79000	158.0	298.722242301	255	0:14
80000	160.0	312.524278007	255	0:13
81000	162.0	288.532317408	255	0:12
82000	164.0	311.714008112	255	0:12
83000	166.0	312.352242718	255	0:11
84000	168.0	296.204101945	255	0:10
85000	170.0	313.090484569	255	0:10
86000	172.0	298.193547877	255	0:09
87000	174.0	297.49275318	255	0:08
88000	176.0	307.210156399	255	0:08
89000	178.0	295.392446535	255	0:07
90000	180.0	294.23406572	255	0:06
91000	182.0	317.272963138	255	0:06
92000	184.0	312.412998773	255	0:05
93000	186.0	290.701264865	255	0:04
94000	188.0	308.933957419	254	0:04
95000	190.0	305.583958019	254	0:03
96000	192.0	310.446224968	254	0:02
97000	194.0	297.386872871	254	0:02
98000	196.0	309.416573594	254	0:01
99000	198.0	289.602514752	254	0:00
100000	200.0	266.420975813	254	0:00

In [3]: