Klauda Lab Simulation Wiki

This page is read only. You can view the source, but not change it. Ask your administrator if you think this is wrong.
[[bilayer simulations]]


======PMF of Permeation for Small Molecules======

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''

=====Building Systems=====

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. 


=====Umbrella Sampling=====

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
=====WHAM for PMF=====

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.