Author: Eric Wang

Please email me at ericzwang@gmail.com if something isn't clear

This procedure is used to study permeation of a small molecule through a bilayer. Conventional MD may not be satisfactory due to large energy barriers which hinder ergodicity. Therefore, we restrain the molecule at a certain plane and calculate the free energy at that plane. Many independent replicates are run to sample the full system. I use CHARMM instead of NAMD because restraints are easier to use in CHARMM.

This is the general procedure
1. Run in MD until equilibration
2. Take final structure and insert small molecule
3. Run initial restrained MD to get molecule to proper position
4. Run restrained MD
5. Get PMF with WHAM

At this point, I will assume you have equilibrated the system in MD, the instructions for which can be found in this wiki.

Before beginning, you need to plan out how you will restrain the molecules in each window. Umbrella sampling is expensive, so we insert multiple molecules into a system and hold them in different planes. This way, we can obtain multiple samples from one replicate. Look at the following excel sheet for an example

umb-zvalues.xlsx

Here, I have 4 ETOH molecules. ccycle is the window number, which I named after looking through some CGUI scripts. zmax is the maximum distance I sample at. This will vary depending on the thickness of your system. However, you need to have points in the water because these become reference values to 0 the PMF. You may notice that there are no negative values because I assume symmetry between leaflets. Make sure to properly space apart the molecules in each window. Big molecules will need more space, so talk to Dr. Klauda if necessary.

Refer to the following path in DT1 as a reference: /lustre/ewang125/sc/part1/tri-cer16-32-umb

Recenter your system using the instructions found in the EDP section.

Use the following scripts in recenter directory to extract the membrane and water. You need to change dyn# and frame#. The cryst.str file should be the same as for recentering.
getframewat.inp.gz getframelip.inp.gz

After extracting, move the getframe*.crd and getframe*.psf files to an umbrella sampling directory. Create a new file getframe.str like in the example and fill in the proper info. The etoh.crd file is also new and is used to insert ethanol into the system.

step5_assembly.inp has been modified. It takes the getframe* and etoh.crd files and combines them. You should adjust the numbers as necessary and look at the system in VMD to make sure it looks OK (small molecules are in system and are properly spaced). Adjust cryst.str in umbstrt directory.

Adjust job name and CHARMM path in step5.csh. Type ./step5.scr to run. If all goes well, the step5* files will be created.

Now, each window will go through an initial run to place the molecules in the right plane. This occurs in the umbstrt directory.

1. Adjust files in psfcrd.str, toppar.str
2. Adjust umbstrt.inp. Change constraint calculation if necessary. Change temperature. Change force constants if necessary (probably depends on molecule size). Adjust nstep if needed.
3. Adjust window range in sbatch.csh
4. Submit sbatch.csh as ./sbatch.csh 1

Look at a few trajectories and verify that molecules have moved into proper planes. Then copy win* directories into umb directory. Remove files in win* directories except for .rea file.

Next we run the actual umbrella sampling in the umb directory
1. Adjust umb.inp similarly to how umbstrt.inp was adjusted. Talk to Dr. Klauda about what force constant is needed. I used 0.75 because 1.5 is in the literature, and CHARMM doesn't multiply by 1/2.
2. Adjust final run in sdyn.csh in the line $nxt > 10.
3. Adjust cryst.str
4. Adjust window range in sbatch.csh.
5. Submit by typing ./sbatch.csh 1

Note: The *.csh files are written to be run on Bridges. sdyn.csh checks if all windows have a .txt file and submits if they do. This method helps identify jobs that randomly fail which is somewhat common on Bridges.
Note: The sbatchnew.csh file is for when a single job fails. I adjust the window # and type ./sbatchnew.csh # where # is the run number.
Note: dirslurm is just a directory for containing slurm output files

Now use WHAM to extract a PMF. I use the code written at http://membrane.urmc.rochester.edu/sites/default/files/wham/doc.html. This method takes in a position time series and uses it to calculate the PMF.
Why don't I use WHAM written in CHARMM? Because there's very little documentation. The documentation in the above link is good and you could figure out how to use it yourself.

Each molecule gets its own directory. For example, ETOH1 is in the e1 directory

Position time series
1. Adjust pos*.inp. Adjust path to psf,crd,dcd files. Adjust molecule segid if necessary. Increase mxa,mxs,mxt if necessary.
2. Adjust pos*.csh. Adjust run numbers and window range.
3. Submit by typing ./pos.scr
The pos*.dat files are the position time series used to calculate PMF.

WHAM
1. To do WHAM, we need a meta file. It's called meta. Adjust the path name to pos*.dat. Adjust the z position. Adjust the force constant (multiply by 2 because this code uses 1/2 factor).
2. Adjust wham.scr. First number is minimum window (0 in mine). Second is maximum window (40 in mine). Third is # window (41). Fourth is tolerance (use .0001). Fifth is temperature. Sixth should be 0.
3. Execute by typing ./wham.scr. It will take only a few seconds. Output is in pmf
4. Repeat in each directory

Excel
1. The bottom section of the pmf file is the energy against window number. Take this section. Using Excel, reorder the energies based on z position.
2. Zero your PMF using the bulk water values.
3. Calculate average and standard error using replicates. If the standard error is too large, you may need to run another replicate. Make sure your PMF makes intuitive sense. Polar molecules should have high energy in the hydrophobic region, for example.