Section 2: Build the windows and equilibrate the system
The APR scripts will build a series of umbrella sampling windows that belong to at least two phases: the attach phase and the pulling phase. A release phase may be necessary if conformational restraints are applied. The number of windows in the attach phase is equal to the number of elements in the attach_list in apr.in. In this tutorial, we have 14 numbers in the attach_list, and therefore, the windows in the attach phase range from a00 to a13. In the attach phase, the guest always stays inside the host cavity, and the restraints for the guest molecule are gradually increased from a00 to a13. For instance, the 11th element in the attach_list is 24.40, so only a weight of 24.40% of the restraints will be imposed in window a10. The reason why only 74.80% of the restraint is imposed in a13, instead of 100% is that the first window of the pulling phase: p00, is essentially the last window of the attach phase, and there is no need to simulate the same window twice. The full strengths of the restraints are listed as the distance_force and angle_force entries in apr.in.
Following the attach phase, the ligand will be pulled out gradually during the pulling phase. Similarly, the number of windows in the pulling phase is determined by the length of the translate_list in apr.in. In this case there are 46 pulling windows ranging from p00 to p45, and each number in translate_list indicates the distance between the guest and its original position. Umbrella sampling windows are spaced at 0.4 Å intervals, and extended to a maximum distance of 18 Å. We assume that at this distance, the interaction between the host and the guest is negligible. When conformational restraints are used, an additional series of "release" windows will be appended.
To build the windows and equilibrate, issue the command
(in ./APR)
python2 apr.py eq -i apr.in -s continue
A folder called "windows" and all 60 subfolders (a00-a13, p00-p45) will be created immediately. Files included in "param_files" and "input_files" folder in addition to align_z.pdb are copied over to each window subfolder.
The apr.py script simulates the equilibration and production phases separately. Writing "eq" to the command line tells the program to run the equilibration phase. The flag "-i" indicates the APR user input file and the flag "-s" provides the option of two modes: "overwrite" and "continue". When "overwrite" is used, all files and subfolders in the preexisting windows folder will be removed, so that the building process can start from scratch. This is useful when you make changes to your settings or input files, and you want to start over the equilibration. The "continue" mode will keep the previous files, and pick up from where the equilibration stops. This can be useful when the simulations are disrupted abruptly yet no changes are made to the settings or parameters. Both "continue" and "overwrite" can be used for the first time run.
Once you issue the command, the screen will show the script conducting the following actions:
Setting up the APR framework ...
Continue mode enabled. Checking previous simulations...
pmemd.cuda is on.
Translational and rotational restraints will be imposed.
No conformational restraints applied.
Preparing folder windows/a00 2016-11-20 02:18:51 PM
Solvating...
Adding restraints...
Equilibrating...
In brief, apr.py imports functions from a group of "child" scripts, each carrying out a unique task as follow:
It takes about 3 minutes to finish equilibration in one window using one single GPU (GeForce GTX 980 or alike). If you would like to continue the tutorial without waiting for the equilibration and the following production phase to finish, a completed windows folder is available as apr_preRun.tar.gz. You can download and extract it to continue with this tutorial without running these equilibrations or productions yourself (Delete or rename the existing windows folder in your directory before you extract the tar.gz file):
(in ./APR)
tar zxvf apr_preRun.tar.gz
The disang.rest file in each window provides a full view of the values of all restraints, as well as which subsets of atoms they are imposed on. We strongly recommend you to take a close look at this file in windows at different phases so that you can gain a better understanding of the method. For instance, the disang.rest file in window a00 looks like this:
(in ./APR/windows/a00)
#AnchorAtoms :OCT@C11 :OCT@C23 :OCT@C51 :MOL@N1 :MOL@C6 #Type a #Weight 0.00 #TransDist 0.00 #Scale 0.002
&rst iat= 200,13, r1=0.0000, r2=4.0097, r3=4.0097, r4=999.0000, rk2=5.0000000, rk3=5.0000000, &end #HostTR
...
&rst iat= 200,192, r1=0.0000, r2=6.0000, r3=6.0000, r4=999.0000, rk2=0.0000000, rk3=0.0000000, &end #GuestTR
...
The first line of the disang.rest file prints out the five atoms we picked in the last section. The letter "a" following the "Type" keyword indicates that the current window belongs to the attach phrase. The value following the "Weight" keyword controls the force constant of guest restraints in the current window. There are no restraints imposed on the guest molecule in the a00 window, which is also called the "bound state", and that is why the weight value is zero. In addition, the "TransDist" keyword defines the pulling distance of the guest, and is always zero in all windows of the attach phase (i.e. in windows a00-a13).
Line 2 to line 7 ended with a keyword "HostTR" provide the values for six host restraints. Since the host restraints are kept constant during the entire APR process, those lines should look exactly the same in all 60 windows. The values of r1, r2, r3 and r4 together specify the range while rk2 and rk3 define the strength of the harmonic restraints.
Values of three guest restraints are provided in the last three lines of the disang.rest file, all ended with a keyword "GuestTR". The values of rk2 and rk3 are both zero because guest restraints are not applied in the bound state.
You can also load solvated.pdb in VMD in each window to see how the guest is gradually pulled away.
The equilibration phase consists of the following:
(1) Running 1 ps NVT at 10 K, using the input file therm1.in;
(2) Heating the system from 10 K to 298.15 K within 100 ps, using therm2.in;
(3) Equilibrating the system under constant pressure, using eqnpt.in. This step contains 50 NPT cycles; The simulation time of
each cycle can be controlled via the eq_dt and eq_nstlim options in apr.in.
Make sure you have the output file eqnpt50.rst7 in all 60 windows before we get the production phase in the next section. This will be the starting input file (same file as rst.00) for the following production runs in each window.