Skip to main content

Enhanced sampling in the Schrödinger suite

 Introduction

In the past I've been using the Schrödinger Molecular Modeling suite quite a bit and I had the chance of working with some talented guys that have taught me a great deal about how that's working. 

One of the hidden gems of it is the ability of doing enhanced sampling using Desmond molecular dynamics code on GPU.

Desmond is a great piece of software but it is closed source and the users base consists mostly of people who are using it via the Maestro GUI to run FEP+, MD or some pre-canned metadynamics workflows. Therefore is no surprise there's very little out there on how to access many capabilities which sits there. 

In this blog entry I'd like to explain a little of how it works and give you an initial glimpse of the capabilities of it.  The case below is worked up using the version 23-1 of the suite but this is a part of the code which is rarely touched so it should be fairly retrocompatible (within a limit).

This blog entry is for you if you have ever found yourself doing some metadynamics but not really pleased with the control you have on walls, or perhaps the CVs and you need more flexibility. Similarly, If you have thought that "Hey, how's it possible that I cannot run a simple umbrella sampling with this state-of-the-art piece of kit?" this is for you! Steered-MD, Random Accelerated MD and quite a number of other algorithms also fall in this category.

Workflow

What I want to show you is how you can start from a "pre-canned" metadynamics example produced by the Maestro and transform in a way to expose the so-called m-expression, which is a cool scripting language that allows you to script up a bunch of things. So these are the steps

  1. Produce a metadynamics input

  2. Process it with a script (link provided) that will expose the m-expression

  3. Run your modified metadynamics 

1- Produce a metadynamics input

I assume you already know how to do this. So you just need to prepare a system in Maestro and use the tutorial to setup a metadynamics run.  Then, instead of clicking on "Run" you click on the cog-wheel symbol and you click on "Write". This should be able to write for you the files without running them, leaving this to you. 

At this stage you may have a folder structure like this one below

desmond_metadynamics_job_1
|-desmond_metadynamics_job_1.cfg
|-desmond_metadynamics_job_1.cms
|-desmond_metadynamics_job_1.msj
|-desmond_metadynamics_job_1.sh

You are now ready to explore the magics! :) 

2- Process with script to get the m-expression

Now you need to download the retrieve_m_expression_from_msj.py  script from here. 

Save it in the folder that contains the desmond_metadynamics_job_1 folder and execute the command below. I assume you work under linux in the following and you have the $SCHRODINGER environment variable already set correctly.


$SCHRODINGER/run retrieve_m_expression_from_msj.py desmond_metadynamics_job_1/desmond_metadynamics_job_1.msj desmond_metadynamics_job_1/desmond_metadynamics_job_1.cms newmsj.msj mexpr.m

The script prints out something like this

Input msj  desmond_metadynamics_job_1/desmond_metadynamics_job_1.msj
Input cms  desmond_metadynamics_job_1/desmond_metadynamics_job_1.cms
Output msj  newmsj.msj
Output mexpression  mexpr.m
Stage  0 task
Stage  1 simulate
Stage  2 simulate
Stage  3 simulate
Stage  4 simulate
Stage  5 simulate
Stage  6 simulate
This stage has metadynamics!
Stage  7 analysis
File written  /tmp_mnt/filer1/unix/dbranduardi/Projects/schrod_tests/blog_material/1_producing_mexpression/desmond_metadynamics_job_1/newmsj.msj
File written  /tmp_mnt/filer1/unix/dbranduardi/Projects/schrod_tests/blog_material/1_producing_mexpression/desmond_metadynamics_job_1/mexpr.m

What happened is that you now have a new MSJ file (newmsj.msj) and a M-expression file (mexpr.m) : these are the files we're interested in and we want to instruct Desmond to use them instead of the original desmond_metadynamics_job_1.msj.

But first... what has happened?

We can have a quick look in the msj section regarding the metadynamics:

simulate {
   cfg_file = "desmond_metadynamics_job_1.cfg"
   jobname  = "$MAINJOBNAME"
   dir      = "."
   compress = ""
   meta = {
      cv = [
         {atom = ["atom.n 3295,3296,3297,3298,3299,3300,3301,3302,3303,3304,3305,3306,3307,3308,3309,3310,3311,3312,3313,3314,3315,3316,3317,3318,3319,3320,3321,3322,3323,3324,3325,3326,3327,3328,3329,3330,3331,3332,333
3,3334,3335,3336,3337,3338,3339,3340,3341,3342,3343" "atom.n 8,23,38,54,69,83,103,113,133,150,160,175,194,204,221,240,257,268,287,306,325,339,353,373,394,405,419,441,456,475,495,514,538,553,572,591,602,616,627,638,650,6
60,679,691,713,732,756,777,792,803,822,836,848,862,873,895,914,926,937,944,966,981,1000,1017,1036,1050,1069,1088,1102,1116,1138,1155,1167,1191,1205,1224,1238,1257,1273,1285,1299,1306,1325,1332,1349,1363,1385,1395,1407,1
426,1445,1459,1473,1492,1499,1513,1532,1542,1564,1575,1582,1596,1618,1628,1648,1665,1680,1690,1709,1726,1736,1743,1753,1765,1784,1795,1812,1831,1838,1855,1875,1882,1898,1905,1925,1946,1957,1967,1988,2007,2023,2033,2048,
2070,2086,2100,2116,2135,2149,2171,2187,2201,2213,2225,2240,2257,2278,2288,2312,2327,2338,2349,2359,2366,2373,2384,2404,2418,2434,2458,2472,2484,2498,2505,2520,2534,2551,2558,2582,2589,2603,2625,2641,2660,2679,2695,2714
,2736,2751,2763,2780,2794,2809,2830,2849,2864,2879,2903,2927,2946,2968,2983,3002,3018,3040,3062,3078,3089,3106,3126,3145,3152,3173,3187,3206,3220,3239,3259,3275,3291" ]
          type = dist
          width = 0.05
         }
      ]
      cv_name = "$JOBNAME$[_replica$REPLICA$].cvseq"
      first = 0.0
      height = 0.03
      interval = 0.09
      name = "$JOBNAME$[_replica$REPLICA$].kerseq"
   } 
           backend ={
               # set cvseq interval to trajectory output
               force.term.ES.interval = '@*.*.*.*.trajectory.interval'
           }
           checkpt.write_last_step = yes
        
}

Now the equivalent section in stage 7 looks pretty much like this:

simulate {
  backend = {
     force = {
        term = {
           ES = {
              interval = "@*.*.*.*.trajectory.interval"
           }
        }
     }
  }
  cfg_file = "desmond_metadynamics_job_1.cfg"
  checkpt = {
     write_last_step = true
  }
  compress = ""
  dir = "."
  jobname = "$MAINJOBNAME"
  meta = FILE
  meta_file = "/tmp_mnt/filer1/unix/dbranduardi/Projects/schrod_tests/blog_material/1_producing_mexpression/desmond_metadynamics_job_1/mexpr.m"
}


Basically, instead of specifying the metadynamics section in the msj, I instructed Desmond (or better, the desmond driver... but this is another story) to read the metadynamics from an external file (meta=FILE) and I provide the path of the file. 

Interestingly that file looks very different:

declare_meta(
    dimension = 1,
    cutoff    = 9.000000,
    first     = 0.000000,
    interval  = 0.090000,
    name      = "$JOBNAME.kerseq",
    initial   = "");

declare_output(
    name = "$JOBNAME.cvseq",
    first = 0.000000,
    interval= 0.090000);


# distance definition for group of atoms
cv_00_g0 = center_of_mass ( atomsel("atom. 3295, 3296, 3297, 3298, 3299, 3300, 3301, 3302, 3303, 3304, 3305, 3306, 3307, 3308, 3309, 3310, 3311, 3312, 3313, 3314, 3315, 3316, 3317, 3318, 3319, 3320, 3321, 3322, 3323, 33
24, 3325, 3326, 3327, 3328, 3329, 3330, 3331, 3332, 3333, 3334, 3335, 3336, 3337, 3338, 3339, 3340, 3341, 3342, 3343") );
cv_00_g1 = center_of_mass ( atomsel("atom. 8, 23, 38, 54, 69, 83, 103, 113, 133, 150, 160, 175, 194, 204, 221, 240, 257, 268, 287, 306, 325, 339, 353, 373, 394, 405, 419, 441, 456, 475, 495, 514, 538, 553, 572, 591, 602
, 616, 627, 638, 650, 660, 679, 691, 713, 732, 756, 777, 792, 803, 822, 836, 848, 862, 873, 895, 914, 926, 937, 944, 966, 981, 1000, 1017, 1036, 1050, 1069, 1088, 1102, 1116, 1138, 1155, 1167, 1191, 1205, 1224, 1238, 12
57, 1273, 1285, 1299, 1306, 1325, 1332, 1349, 1363, 1385, 1395, 1407, 1426, 1445, 1459, 1473, 1492, 1499, 1513, 1532, 1542, 1564, 1575, 1582, 1596, 1618, 1628, 1648, 1665, 1680, 1690, 1709, 1726, 1736, 1743, 1753, 1765,
 1784, 1795, 1812, 1831, 1838, 1855, 1875, 1882, 1898, 1905, 1925, 1946, 1957, 1967, 1988, 2007, 2023, 2033, 2048, 2070, 2086, 2100, 2116, 2135, 2149, 2171, 2187, 2201, 2213, 2225, 2240, 2257, 2278, 2288, 2312, 2327, 23
38, 2349, 2359, 2366, 2373, 2384, 2404, 2418, 2434, 2458, 2472, 2484, 2498, 2505, 2520, 2534, 2551, 2558, 2582, 2589, 2603, 2625, 2641, 2660, 2679, 2695, 2714, 2736, 2751, 2763, 2780, 2794, 2809, 2830, 2849, 2864, 2879,
 2903, 2927, 2946, 2968, 2983, 3002, 3018, 3040, 3062, 3078, 3089, 3106, 3126, 3145, 3152, 3173, 3187, 3206, 3220, 3239, 3259, 3275, 3291") );
cv_00 = norm(min_image(cv_00_g0 - cv_00_g1));
print ("cv_00", cv_00);
# the width for cv_00 will be set to: 0.05


# height used for this run is: 0.030000
meta(0, array(0.03, 0.05),
        array(cv_00));

and here is where the magics happen. This is a scripting language much more powerful that the one in the msj. Here you can define the CVs yourself: for example the distance between two center of mass (ligand and protein) is here defined by first defining the two COMs and then calculating the difference vector, applying minimal image convention and then calculating the norm to that vector. So you're in the belly of the metadynamics. 

Additionally here you can add and mix all the restraints you want. But beware: you need to know what you're doing! So do your homework before trying to reinvent the next greatest CV or enhanced sampling algo. 

For the time being you can try to play around with some parameters (like hill height 0.03 or width 0.05).

In chapter 11 of the Desmond manual you can find a more details about the syntax of the m-expression. Not a reading for the faint hearted.

3-Run your modified metadynamics 

In order to do that we need to hack the shell script which is produced when clicking on the "Write" in Maestro. 

The content of the original file is something like this below


"${SCHRODINGER}/utilities/multisim" -JOBNAME desmond_metadynamics_job_1 -HOST HOSTNAME -maxjob 1 -cpu 1 -m desmond_metadynamics_job_1.msj -c desmond_metadynamics_job_1.cfg -description Metadynamics desmond_metadynamics_jo
b_1.cms -mode umbrella -o desmond_metadynamics_job_1-out.cms -lic DESMOND_GPGPU:16

We just need to change so that the newmsj.msj file will be used


"${SCHRODINGER}/utilities/multisim" -JOBNAME desmond_metadynamics_job_1 -HOST HOSTNAME -maxjob 1 -cpu 1 -m newmsj.msj -c desmond_metadynamics_job_1.cfg -description Metadynamics desmond_metadynamics_jo
b_1.cms -mode umbrella -o desmond_metadynamics_job_1-out.cms -lic DESMOND_GPGPU:16

Then you can run as usual and see for yourself that the changed parameters are in effect.

Labels:

Comments

Post a Comment