From e2a5ad12b8f101ce427c6f37cb19c4a49990d3bb Mon Sep 17 00:00:00 2001 From: juleaduriz Date: Tue, 17 Dec 2024 10:20:38 +0100 Subject: [PATCH] Deployed 8905d23 with MkDocs version: 1.6.1 --- index.html | 2 +- metadynamics-with-Q/index.html | 148 +++++++++++---------------------- search/search_index.json | 2 +- 3 files changed, 51 insertions(+), 101 deletions(-) diff --git a/index.html b/index.html index 21dc255..9b38377 100644 --- a/index.html +++ b/index.html @@ -164,5 +164,5 @@

Running simulat diff --git a/metadynamics-with-Q/index.html b/metadynamics-with-Q/index.html index 7ddec27..be77b2a 100644 --- a/metadynamics-with-Q/index.html +++ b/metadynamics-with-Q/index.html @@ -95,31 +95,21 @@

Defining native contacts (Q) as collective variable

In order to understand how to implement Q as a collective variable in PLUMED one needs to understand how Q is defined:

-

- -

- -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_, -is the distance between _i_ and _j_ in the native state, +

+

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, +is the distance between i and j in the native state, is a smoothing parameter taken to be 5 and the factor 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.

+

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.

pdb = md.load("Heliat-1/Heliat-1.pdb")
- - -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): - - +

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):

ca = {}
 
 for i in range(1,17):
@@ -128,14 +118,10 @@ 
Defining native conta
     else:
         ca[top.select("residue %i and name CA"%i)[0]] = i
- - -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: - - +

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:

NATIVE_CUTOFF = 1.0  # nanometers
 
 alpha_carbons = [x for x in ca.keys()]
@@ -149,37 +135,25 @@ 
Defining native conta
 
 native_contacts = Ca_pairs[Ca_pairs_distances <= NATIVE_CUTOFF]
- - -In this case, I am considering as a native contact any interaction between to C$`\alpha`$ atoms that are +

In this case, I am considering as a native contact any interaction between to C$\alpha$ 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 native_contacts array, which should look like this:

array([[  3,  38],
        [  3,  59],
        [  3,  69],
        [  7,  59],
 ...
- - -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: - - +

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:

pairs = []
 
 for i in range(len(native_contacts)):
     pairs.append([ca[native_contacts[i,0]], ca[native_contacts[i,1]]])
- - -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: - - +

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:

file = open("Heliat-1/CA-list.dat", "w")
 
 file.write("MOLINFO STRUCTURE=Heliat-1.pdb \n")
@@ -194,54 +168,41 @@ 
Defining native conta
 
 file.close()
- - -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 +

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 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: - - +$\beta$, $\lambda$, r $^{0}_{ij}$ and weight (which should be 1 divided by the number of native contacts). +The synthax is the following:

__ARG1__: CONTACTMAP ...
    ATOMS=atom1,atom2 SWITCH1={Q R_0=1.0 BETA=50.0 LAMBDA=1.2 REF=__ARG2__} WEIGHT=__ARG3__
    ...
    SUM
 ...
- -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: - - +

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:

metad: METAD ARG=cmap ...
    PACE=500 HEIGHT=1.2
    SIGMA=__ARG4__
    FILE=HILLS GRID_MIN=0 GRID_MAX=1
 ...
- - -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: - - +

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:

file = open("Heliat-1/plumed.dat", "w")
 
 w = 1 / len(pairs)
@@ -269,30 +230,19 @@ 
Defining native conta
 
 file.close()
- - -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: - - +

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:

gmx_mpi mdrun -s *tpr -plumed plumed.dat -ntomp 1
- - -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: - - +

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:

plumed sum_hills --hills HILLS
- - -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 +

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