Skip to content

Commit

Permalink
Deployed 8905d23 with MkDocs version: 1.6.1
Browse files Browse the repository at this point in the history
  • Loading branch information
@julenaduriz
julenaduriz committed Dec 17, 2024
1 parent b88d4e3 commit e2a5ad1
Show file tree
Hide file tree
Showing 3 changed files with 51 additions and 101 deletions.
2 changes: 1 addition & 1 deletion index.html
View file
Original file line number Diff line number Diff line change
Expand Up @@ -164,5 +164,5 @@ <h3 id="running-simulations-with-plumed-at-dipc"><a href="dipc/">Running simulat

<!--
MkDocs version : 1.6.1
Build Date UTC : 2024-12-17 09:19:08.052268+00:00
Build Date UTC : 2024-12-17 09:20:38.900302+00:00
-->
148 changes: 49 additions & 99 deletions metadynamics-with-Q/index.html
View file
Original file line number Diff line number Diff line change
Expand Up @@ -95,31 +95,21 @@
<h2 id="defining-native-contacts-q-as-collective-variable">Defining native contacts (Q) as collective variable</h2>
<p>In order to understand how to implement Q as a collective variable in PLUMED
one needs to understand how Q is defined:</p>
<p>
<img src="https://latex.codecogs.com/svg.image?Q(X)=\frac{1}{N}\sum_{\left(i,i\right)}\frac{1}{1&plus;exp\left[\beta\left(r_{ij}(X)-\lambda&space;r^{0}_{ij}\right)\right]}">
<p>

where the sum runs over the _N_ pairs of native contacts _(i,j)_, _r_ $`_{ij}`$ _(X)_ is the
distance between _i_ and _j_ in configuration _X_, <img src="https://latex.codecogs.com/svg.image?r_{ij}^{0}">
is the distance between _i_ and _j_ in the native state, <img src="https://latex.codecogs.com/svg.image?\beta">
<p><center><img src="https://latex.codecogs.com/svg.image?Q(X)=\frac{1}{N}\sum_{\left(i,i\right)}\frac{1}{1&plus;exp\left[\beta\left(r_{ij}(X)-\lambda&space;r^{0}_{ij}\right)\right]}"></p>
<p>where the sum runs over the <em>N</em> pairs of native contacts <em>(i,j)</em>, <em>r</em> $<code>_{ij}</code>$ <em>(X)</em> is the
distance between <em>i</em> and <em>j</em> in configuration <em>X</em>, <img src="https://latex.codecogs.com/svg.image?r_{ij}^{0}">
is the distance between <em>i</em> and <em>j</em> in the native state, <img src="https://latex.codecogs.com/svg.image?\beta">
is a smoothing parameter taken to be 5 <img src="https://latex.codecogs.com/svg.image?\AA ^{-1}">
and the factor <img src="https://latex.codecogs.com/svg.image?\lambda"> accounts for fluctuations when the
contact is formed, taken to be 1.2.

Therefore, the first step is to compute the number of native contacts of our system. To do so,
contact is formed, taken to be 1.2.</p>
<p>Therefore, the first step is to compute the number of native contacts of our system. To do so,
I load the native structure. In the case of this example, this native structure has been loaded
as a PDB file, although one could select the file type of their preference.


as a PDB file, although one could select the file type of their preference.</p>
<pre><code>pdb = md.load(&quot;Heliat-1/Heliat-1.pdb&quot;)
</code></pre>


For this particular example, I will consider my native contacts to be exclusively given by
C$`\alpha`$ atoms. Therefore, I make a dictionary to save the atom index of each C$`\alpha`$
atom alongside their residue index number (this will be usefull in a future step):


<p>For this particular example, I will consider my native contacts to be exclusively given by
C$<code>\alpha</code>$ atoms. Therefore, I make a dictionary to save the atom index of each C$<code>\alpha</code>$
atom alongside their residue index number (this will be usefull in a future step):</p>
<pre><code>ca = {}

for i in range(1,17):
Expand All @@ -128,14 +118,10 @@ <h2 id="defining-native-contacts-q-as-collective-variable">Defining native conta
else:
ca[top.select(&quot;residue %i and name CA&quot;%i)[0]] = i
</code></pre>


Notice that I am working with an Heliat, which is cyclated using the sidechain of a CYS and an
ACE. The C atom of this ACE is involved in the formation of $`\alpha`$-helix and therefore I
define it as an C$`\alpha`$ atom. Once the atom selection is done, I compute the number of native
contacts:


<p>Notice that I am working with an Heliat, which is cyclated using the sidechain of a CYS and an
ACE. The C atom of this ACE is involved in the formation of $<code>\alpha</code>$-helix and therefore I
define it as an C$<code>\alpha</code>$ atom. Once the atom selection is done, I compute the number of native
contacts:</p>
<pre><code>NATIVE_CUTOFF = 1.0 # nanometers

alpha_carbons = [x for x in ca.keys()]
Expand All @@ -149,37 +135,25 @@ <h2 id="defining-native-contacts-q-as-collective-variable">Defining native conta

native_contacts = Ca_pairs[Ca_pairs_distances &lt;= NATIVE_CUTOFF]
</code></pre>


In this case, I am considering as a native contact any interaction between to C$`\alpha`$ atoms that are
<p>In this case, I am considering as a native contact any interaction between to C$<code>\alpha</code>$ atoms that are
at least at 2 residues apart and that their distance is lower than 1 nm. These native contacts are
then stored in the _native\_contacts_ array, which should look like this:


then stored in the <em>native_contacts</em> array, which should look like this:</p>
<pre><code>array([[ 3, 38],
[ 3, 59],
[ 3, 69],
[ 7, 59],
...
</code></pre>


where each number belongs to the atom index of each C$`\alpha`$. Since I want to translate this to a synthax
PLUMED can understand, I translate this contacts from atom indexes to residue indexes as follows:


<p>where each number belongs to the atom index of each C$<code>\alpha</code>$. Since I want to translate this to a synthax
PLUMED can understand, I translate this contacts from atom indexes to residue indexes as follows:</p>
<pre><code>pairs = []

for i in range(len(native_contacts)):
pairs.append([ca[native_contacts[i,0]], ca[native_contacts[i,1]]])
</code></pre>


Once this step is performed, I am now ready to start writing the first actual PLUMED input files. To do so,
the first thing I do is to create a file where I select each C$`\alpha`$ internally with PLUMED. I will call
this file `CA-list.dat` and it will contain the following:


<p>Once this step is performed, I am now ready to start writing the first actual PLUMED input files. To do so,
the first thing I do is to create a file where I select each C$<code>\alpha</code>$ internally with PLUMED. I will call
this file <code>CA-list.dat</code> and it will contain the following:</p>
<pre><code>file = open(&quot;Heliat-1/CA-list.dat&quot;, &quot;w&quot;)

file.write(&quot;MOLINFO STRUCTURE=Heliat-1.pdb \n&quot;)
Expand All @@ -194,54 +168,41 @@ <h2 id="defining-native-contacts-q-as-collective-variable">Defining native conta

file.close()
</code></pre>


This file contains the information of where the atom selection is being made from (PDB file in this case) and
makes the atom selection internally making use of the MDTraj synthax thanks to the `@mdt` flag. Once this step
<p>This file contains the information of where the atom selection is being made from (PDB file in this case) and
makes the atom selection internally making use of the MDTraj synthax thanks to the <code>@mdt</code> flag. Once this step
is done, I will write the actual PLUMED input file that contains the information for the metadynamics
simulation. In order to fill this file, one should understand how Q should be defined as a collective variable (CV)
in PLUMED. The first step is to define a contact map where the contacts are your native contacts. Inside the
definition of this contact map, one should define five important parameters: cutoff distance of a native contact,
$`\beta`$, $`\lambda`$, _r_ $`^{0}_{ij}`$ and weight (which should be 1 divided by the number of native contacts).
The synthax is the following:


$<code>\beta</code>$, $<code>\lambda</code>$, <em>r</em> $<code>^{0}_{ij}</code>$ and weight (which should be 1 divided by the number of native contacts).
The synthax is the following:</p>
<pre><code>__ARG1__: CONTACTMAP ...
ATOMS=atom1,atom2 SWITCH1={Q R_0=1.0 BETA=50.0 LAMBDA=1.2 REF=__ARG2__} WEIGHT=__ARG3__
...
SUM
...
</code></pre>

where `__ARG1__` should be substituted by the name with which one wants to define the contact map (for now on,
I will be using `cmap`), `__ARG2__` should be substituted by the reference distance (_r_ $`^{0}_{ij}`$) and
`__ARG3__` should be substituted by the weight of each distance (1 divided by the number of native contacts).
In order to let PLUMED know that it should perform a calculation of Q, we add the Q inside the `SWITCH` function
and add the `SUM` line at the end of the definition.

Once the function is defined, I will add the following lines to define this contact map (`cmap` in my case) as
a CV:


<p>where <code>__ARG1__</code> should be substituted by the name with which one wants to define the contact map (for now on,
I will be using <code>cmap</code>), <code>__ARG2__</code> should be substituted by the reference distance (<em>r</em> $<code>^{0}_{ij}</code>$) and
<code>__ARG3__</code> should be substituted by the weight of each distance (1 divided by the number of native contacts).
In order to let PLUMED know that it should perform a calculation of Q, we add the Q inside the <code>SWITCH</code> function
and add the <code>SUM</code> line at the end of the definition.</p>
<p>Once the function is defined, I will add the following lines to define this contact map (<code>cmap</code> in my case) as
a CV:</p>
<pre><code>metad: METAD ARG=cmap ...
PACE=500 HEIGHT=1.2
SIGMA=__ARG4__
FILE=HILLS GRID_MIN=0 GRID_MAX=1
...
</code></pre>


where `PACE` defines the frequency (in frames) with which a new Gaussian hill will be added,`HEIGHT` defines the hieghts
of the Gaussian hills and `SIGMA` defines the width of the Gaussian hills. Here, `__ARG4__` should be defined depending
on each individual case. PLUMED recommends to set the value of `SIGMA` to be a third of the fluctuations of the CV one
is using. In this particular case, `SIGMA` should be of around 0.005-0.01. Lastly, `GRID_MIN` and `GRID_MAX` define the
lower and upper bound of the grid. The information of the Gaussian hills will then be written to the `HILLS` file. This
will be important in future steps.

In order to create an input file to run a metadynamics calculation with PLUMED which will be named `plumed.dat` I run
the following lines:


<p>where <code>PACE</code> defines the frequency (in frames) with which a new Gaussian hill will be added,<code>HEIGHT</code> defines the hieghts
of the Gaussian hills and <code>SIGMA</code> defines the width of the Gaussian hills. Here, <code>__ARG4__</code> should be defined depending
on each individual case. PLUMED recommends to set the value of <code>SIGMA</code> to be a third of the fluctuations of the CV one
is using. In this particular case, <code>SIGMA</code> should be of around 0.005-0.01. Lastly, <code>GRID_MIN</code> and <code>GRID_MAX</code> define the
lower and upper bound of the grid. The information of the Gaussian hills will then be written to the <code>HILLS</code> file. This
will be important in future steps.</p>
<p>In order to create an input file to run a metadynamics calculation with PLUMED which will be named <code>plumed.dat</code> I run
the following lines:</p>
<pre><code>file = open(&quot;Heliat-1/plumed.dat&quot;, &quot;w&quot;)

w = 1 / len(pairs)
Expand Down Expand Up @@ -269,30 +230,19 @@ <h2 id="defining-native-contacts-q-as-collective-variable">Defining native conta

file.close()
</code></pre>


Once the `plumed.dat` file is created, the metadynamics simulation is ready to be performed. One should run it in an
environment where MDTraj is present (as it is imprescindible to define the C$`\alpha`$ atoms we defined in the `CA-list.dat` file.
The line to run this calculation is the following:


<p>Once the <code>plumed.dat</code> file is created, the metadynamics simulation is ready to be performed. One should run it in an
environment where MDTraj is present (as it is imprescindible to define the C$<code>\alpha</code>$ atoms we defined in the <code>CA-list.dat</code> file.
The line to run this calculation is the following:</p>
<pre><code>gmx_mpi mdrun -s *tpr -plumed plumed.dat -ntomp 1
</code></pre>


Once the calculation is finished, one can extract the evolution of how Q changed during the simulation from the
`COLVAR` file, where the first column gives the time value and the second column gives the value of Q. Regarding
the `HILLS` file, in order to extract information one must run the following line:


<p>Once the calculation is finished, one can extract the evolution of how Q changed during the simulation from the
<code>COLVAR</code> file, where the first column gives the time value and the second column gives the value of Q. Regarding
the <code>HILLS</code> file, in order to extract information one must run the following line:</p>
<pre><code>plumed sum_hills --hills HILLS
</code></pre>


Once this line is run, a file named `fes.dat` will be created, which containts the information of the Free Energy Surface
of your metadynamics calculation.

Lets see
<p>Once this line is run, a file named <code>fes.dat</code> will be created, which containts the information of the Free Energy Surface
of your metadynamics calculation.</p>
<p>Lets see</p>

</div>
</div><footer>
Expand Down
2 changes: 1 addition & 1 deletion search/search_index.json
View file

Large diffs are not rendered by default.

0 comments on commit e2a5ad1

Please sign in to comment.