File: solvation.py.html

package info (click to toggle)
mmtk 2.7.9-1
  • links: PTS, VCS
  • area: main
  • in suites: jessie, jessie-kfreebsd
  • size: 11,788 kB
  • ctags: 6,600
  • sloc: python: 18,050; ansic: 12,400; makefile: 129; csh: 3
file content (289 lines) | stat: -rw-r--r-- 16,157 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289


<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN"
  "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">


<html xmlns="http://www.w3.org/1999/xhtml">
  <head>
    <meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
    
    <title>Solvation of a protein in water &mdash; MMTK User Guide 2.7.7 documentation</title>
    
    <link rel="stylesheet" href="../../_static/default.css" type="text/css" />
    <link rel="stylesheet" href="../../_static/pygments.css" type="text/css" />
    
    <script type="text/javascript">
      var DOCUMENTATION_OPTIONS = {
        URL_ROOT:    '../../',
        VERSION:     '2.7.7',
        COLLAPSE_INDEX: false,
        FILE_SUFFIX: '.html',
        HAS_SOURCE:  true
      };
    </script>
    <script type="text/javascript" src="../../_static/jquery.js"></script>
    <script type="text/javascript" src="../../_static/underscore.js"></script>
    <script type="text/javascript" src="../../_static/doctools.js"></script>
    <link rel="top" title="MMTK User Guide 2.7.7 documentation" href="../../index.html" /> 
  </head>
  <body>
    <div class="related">
      <h3>Navigation</h3>
      <ul>
        <li class="right" style="margin-right: 10px">
          <a href="../../genindex.html" title="General Index"
             accesskey="I">index</a></li>
        <li class="right" >
          <a href="../../py-modindex.html" title="Python Module Index"
             >modules</a> |</li>
        <li><a href="../../index.html">MMTK User Guide 2.7.7 documentation</a> &raquo;</li> 
      </ul>
    </div>  

    <div class="document">
      <div class="documentwrapper">
        <div class="bodywrapper">
          <div class="body">
            
  <div class="section" id="solvation-of-a-protein-in-water">
<h1>Solvation of a protein in water<a class="headerlink" href="#solvation-of-a-protein-in-water" title="Permalink to this headline">ΒΆ</a></h1>
<div class="highlight-python"><table class="highlighttable"><tr><td class="linenos"><div class="linenodiv"><pre> 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95</pre></div></td><td class="code"><div class="highlight"><pre><span class="c"># Solvation of a protein with water.</span>
<span class="c">#</span>
<span class="c"># The solvation procedure consists of three steps:</span>
<span class="c">#</span>
<span class="c"># - The universe containing the protein is scaled up, and the required</span>
<span class="c">#   number of water molecules is added at random positions, but without</span>
<span class="c">#   any overlap between molecules. The system is so dilute that random</span>
<span class="c">#   placements are easily possible.</span>
<span class="c">#</span>
<span class="c"># - The universe is slowly scaled down do its original size, with</span>
<span class="c">#   each scaling step followed by some energy minimization and</span>
<span class="c">#   molecular dynamics steps.</span>
<span class="c">#</span>
<span class="c"># - A molecular dynamics run at constant pressure and temperature is</span>
<span class="c">#   used to put the system into a well-defined thermodynamic state.</span>
<span class="c">#</span>

<span class="kn">from</span> <span class="nn">MMTK</span> <span class="kn">import</span> <span class="o">*</span>
<span class="kn">from</span> <span class="nn">MMTK.Proteins</span> <span class="kn">import</span> <span class="n">Protein</span>
<span class="kn">from</span> <span class="nn">MMTK.ForceFields</span> <span class="kn">import</span> <span class="n">Amber94ForceField</span>
<span class="kn">from</span> <span class="nn">MMTK.Environment</span> <span class="kn">import</span> <span class="n">NoseThermostat</span><span class="p">,</span> <span class="n">AndersenBarostat</span>
<span class="kn">from</span> <span class="nn">MMTK.Trajectory</span> <span class="kn">import</span> <span class="n">Trajectory</span><span class="p">,</span> <span class="n">TrajectoryOutput</span><span class="p">,</span> <span class="n">LogOutput</span>
<span class="kn">from</span> <span class="nn">MMTK.Dynamics</span> <span class="kn">import</span> <span class="n">VelocityVerletIntegrator</span><span class="p">,</span> <span class="n">VelocityScaler</span><span class="p">,</span> \
                          <span class="n">BarostatReset</span><span class="p">,</span> <span class="n">TranslationRemover</span>
<span class="kn">import</span> <span class="nn">MMTK.Solvation</span>

<span class="c"># Create the solute.</span>
<span class="n">protein</span> <span class="o">=</span> <span class="n">Protein</span><span class="p">(</span><span class="s">&#39;bala1&#39;</span><span class="p">)</span>

<span class="c"># Put the solvent in a standard configuration: center of mass at the</span>
<span class="c"># coordinate origin, principal axes of inertia parallel to the coordinate axes.</span>
<span class="n">protein</span><span class="o">.</span><span class="n">normalizeConfiguration</span><span class="p">()</span>

<span class="c"># Define density, pressure,  and temperature of the solvent.</span>
<span class="n">water_density</span> <span class="o">=</span> <span class="mf">1.</span><span class="o">*</span><span class="n">Units</span><span class="o">.</span><span class="n">g</span><span class="o">/</span><span class="n">Units</span><span class="o">.</span><span class="n">cm</span><span class="o">**</span><span class="mi">3</span>
<span class="n">temperature</span> <span class="o">=</span> <span class="mf">300.</span><span class="o">*</span><span class="n">Units</span><span class="o">.</span><span class="n">K</span>
<span class="n">pressure</span> <span class="o">=</span> <span class="mf">1.</span><span class="o">*</span><span class="n">Units</span><span class="o">.</span><span class="n">atm</span>

<span class="c"># Calculate the box size as the boundary box of the protein plus an</span>
<span class="c"># offset. Note: a much larger offset should be used in real applications.</span>
<span class="n">box</span> <span class="o">=</span> <span class="n">protein</span><span class="o">.</span><span class="n">boundingBox</span><span class="p">()</span>
<span class="n">box</span> <span class="o">=</span> <span class="n">box</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span><span class="o">-</span><span class="n">box</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">+</span><span class="n">Vector</span><span class="p">(</span><span class="mf">0.5</span><span class="p">,</span> <span class="mf">0.5</span><span class="p">,</span> <span class="mf">0.5</span><span class="p">)</span>

<span class="c"># Create a periodic universe. The force field is intentionally created with</span>
<span class="c"># a rather small cutoff to speed up the solvation process.</span>
<span class="n">universe</span> <span class="o">=</span> <span class="n">OrthorhombicPeriodicUniverse</span><span class="p">(</span><span class="nb">tuple</span><span class="p">(</span><span class="n">box</span><span class="p">),</span>
                                        <span class="n">Amber94ForceField</span><span class="p">(</span><span class="mf">1.</span><span class="p">,</span> <span class="mf">1.</span><span class="p">))</span>
<span class="n">universe</span><span class="o">.</span><span class="n">protein</span> <span class="o">=</span> <span class="n">protein</span>

<span class="c"># Find the number of solvent molecules.</span>
<span class="k">print</span> <span class="n">MMTK</span><span class="o">.</span><span class="n">Solvation</span><span class="o">.</span><span class="n">numberOfSolventMolecules</span><span class="p">(</span><span class="n">universe</span><span class="p">,</span><span class="s">&#39;water&#39;</span><span class="p">,</span><span class="n">water_density</span><span class="p">),</span>\
      <span class="s">&quot;water molecules will be added&quot;</span>

<span class="c"># Scale up the universe and add the solvent molecules.</span>
<span class="n">MMTK</span><span class="o">.</span><span class="n">Solvation</span><span class="o">.</span><span class="n">addSolvent</span><span class="p">(</span><span class="n">universe</span><span class="p">,</span> <span class="s">&#39;water&#39;</span><span class="p">,</span> <span class="n">water_density</span><span class="p">)</span>
<span class="k">print</span> <span class="s">&quot;Solvent molecules have been added, now shrinking universe...&quot;</span>

<span class="c"># Shrink the universe back to its original size, thereby compressing</span>
<span class="c"># the solvent to its real density.</span>
<span class="n">MMTK</span><span class="o">.</span><span class="n">Solvation</span><span class="o">.</span><span class="n">shrinkUniverse</span><span class="p">(</span><span class="n">universe</span><span class="p">,</span> <span class="n">temperature</span><span class="p">,</span> <span class="s">&#39;solvation.nc&#39;</span><span class="p">)</span>
<span class="k">print</span> <span class="s">&quot;Universe has been compressed, now equilibrating...&quot;</span>

<span class="c"># Set a better force field and add thermostat and barostat.</span>
<span class="c">#</span>
<span class="c"># Note: For efficiency, optimized Ewald parameters should be used</span>
<span class="c"># in a real application. The barostat relaxation time must be</span>
<span class="c"># adjusted to the system size; it should be chosen smaller than</span>
<span class="c"># for a realistic simulation in order to reach the final</span>
<span class="c"># volume faster.</span>
<span class="n">universe</span><span class="o">.</span><span class="n">setForceField</span><span class="p">(</span><span class="n">Amber94ForceField</span><span class="p">(</span><span class="mf">1.4</span><span class="p">,</span> <span class="p">{</span><span class="s">&#39;method&#39;</span><span class="p">:</span> <span class="s">&#39;ewald&#39;</span><span class="p">}))</span>
<span class="n">universe</span><span class="o">.</span><span class="n">thermostat</span> <span class="o">=</span> <span class="n">NoseThermostat</span><span class="p">(</span><span class="n">temperature</span><span class="p">)</span>
<span class="n">universe</span><span class="o">.</span><span class="n">barostat</span> <span class="o">=</span> <span class="n">AndersenBarostat</span><span class="p">(</span><span class="n">pressure</span><span class="p">,</span> <span class="mf">0.1</span><span class="o">*</span><span class="n">Units</span><span class="o">.</span><span class="n">ps</span><span class="p">)</span>

<span class="c"># Create an integrator and a trajectory.</span>
<span class="n">integrator</span> <span class="o">=</span> <span class="n">VelocityVerletIntegrator</span><span class="p">(</span><span class="n">universe</span><span class="p">,</span> <span class="n">delta_t</span><span class="o">=</span><span class="mf">1.</span><span class="o">*</span><span class="n">Units</span><span class="o">.</span><span class="n">fs</span><span class="p">)</span>

<span class="n">trajectory</span> <span class="o">=</span> <span class="n">Trajectory</span><span class="p">(</span><span class="n">universe</span><span class="p">,</span> <span class="s">&quot;equilibration.nc&quot;</span><span class="p">,</span> <span class="s">&quot;w&quot;</span><span class="p">,</span>
                        <span class="s">&quot;Equilibration (NPT ensemble)&quot;</span><span class="p">)</span>

<span class="c"># Start an NPT integration with periodic rescaling of velocities</span>
<span class="c"># and resetting of the barostat. The number of steps required to</span>
<span class="c"># reach a stable volume depends strongly on the system!</span>
<span class="n">output_actions</span> <span class="o">=</span> <span class="p">[</span><span class="n">TrajectoryOutput</span><span class="p">(</span><span class="n">trajectory</span><span class="p">,</span>
                                   <span class="p">(</span><span class="s">&#39;configuration&#39;</span><span class="p">,</span> <span class="s">&#39;energy&#39;</span><span class="p">,</span> <span class="s">&#39;thermodynamic&#39;</span><span class="p">,</span>
                                    <span class="s">&#39;time&#39;</span><span class="p">,</span> <span class="s">&#39;auxiliary&#39;</span><span class="p">),</span> <span class="mi">0</span><span class="p">,</span> <span class="bp">None</span><span class="p">,</span> <span class="mi">100</span><span class="p">),</span>
                  <span class="n">LogOutput</span><span class="p">(</span><span class="s">&quot;equilibration.log&quot;</span><span class="p">,</span> <span class="p">(</span><span class="s">&#39;time&#39;</span><span class="p">,</span> <span class="s">&#39;energy&#39;</span><span class="p">),</span>
                            <span class="mi">0</span><span class="p">,</span> <span class="bp">None</span><span class="p">,</span> <span class="mi">100</span><span class="p">)]</span>
<span class="n">integrator</span><span class="p">(</span><span class="n">steps</span> <span class="o">=</span> <span class="mi">1000</span><span class="p">,</span>
           <span class="n">actions</span> <span class="o">=</span> <span class="p">[</span><span class="n">TranslationRemover</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span> <span class="bp">None</span><span class="p">,</span> <span class="mi">200</span><span class="p">),</span>
                      <span class="n">BarostatReset</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span> <span class="bp">None</span><span class="p">,</span> <span class="mi">20</span><span class="p">),</span>
                      <span class="n">VelocityScaler</span><span class="p">(</span><span class="n">temperature</span><span class="p">,</span> <span class="mf">0.</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="bp">None</span><span class="p">,</span> <span class="mi">20</span><span class="p">)]</span>
           <span class="o">+</span> <span class="n">output_actions</span><span class="p">)</span>

<span class="c"># Close the equilibration trajectory</span>
<span class="n">trajectory</span><span class="o">.</span><span class="n">close</span><span class="p">()</span>
</pre></div>
</td></tr></table></div>
</div>


          </div>
        </div>
      </div>
      <div class="sphinxsidebar">
        <div class="sphinxsidebarwrapper">
  <h3>This Page</h3>
  <ul class="this-page-menu">
    <li><a href="../../_sources/Examples/MolecularDynamics/solvation.py.txt"
           rel="nofollow">Show Source</a></li>
  </ul>
<div id="searchbox" style="display: none">
  <h3>Quick search</h3>
    <form class="search" action="../../search.html" method="get">
      <input type="text" name="q" />
      <input type="submit" value="Go" />
      <input type="hidden" name="check_keywords" value="yes" />
      <input type="hidden" name="area" value="default" />
    </form>
    <p class="searchtip" style="font-size: 90%">
    Enter search terms or a module, class or function name.
    </p>
</div>
<script type="text/javascript">$('#searchbox').show(0);</script>
        </div>
      </div>
      <div class="clearer"></div>
    </div>
    <div class="related">
      <h3>Navigation</h3>
      <ul>
        <li class="right" style="margin-right: 10px">
          <a href="../../genindex.html" title="General Index"
             >index</a></li>
        <li class="right" >
          <a href="../../py-modindex.html" title="Python Module Index"
             >modules</a> |</li>
        <li><a href="../../index.html">MMTK User Guide 2.7.7 documentation</a> &raquo;</li> 
      </ul>
    </div>
    <div class="footer">
        &copy; Copyright 2010, Konrad Hinsen.
      Created using <a href="http://sphinx.pocoo.org/">Sphinx</a> 1.1.3.
    </div>
  </body>
</html>