Klauda Lab Simulation Wiki

Free Energy Perturbation Calculation in CHARMM

Link to FEP files in DT2: /lustre/jbklauda/FEP/FEP_example.zip

To run an FEP calculation, you will need:

1. A coordinate (.crd) file for your molecule in a vacuum, equilibrated. Call it “moleculename_gas.crd,” except with your molecule’s name in the place of “moleculename.”
2. A coordinate file of your molecule in a waterbox that is an appropriate size for the molecule. Call it “moleculename_wbox.crd.” It’s important that this system is NOT equilibrated, but that the waterbox by itself IS. I have provided files for building this system if you need them, and will explain later how to use them. They can be found in the “water” folder in the “example” directory.

Copy the files and follow these steps:

1. Check for proper permission to all .run files. Navigate to the appropriate folder in DT2. Then, you can use the “chmod” command to allow executable permission.
2. Update the path to the CHARMM executable file in the following files:
3. /vac/fep.run, /water/fep.run, /vac/wham_anal.run, /water/wham_anal.run
4. Edit the files /vac/paraset.str and /water/paraset.str to fit your system parameters. Specifically, change molecule or residue names and “nwater,” “temp,” and “boxsize.” The box size should be EXACTLY the size of the waterbox equilibrated before you added the molecule. For example, if my molecule is called DMOE, and I build a 25-Ang. water box for it, I need to then equilibrate that water box and get its size from the restart file. In the “.rst” file produced in a CHARMM equilibration, this will be the line that looks like this: !CRYSTAL PARAMETERS 0.273381676623395D+02 0.000000000000000D+00 0.273381676623395D+02 Where this means your box is 27.3381676623395 Ang. long. Just copy this number and past it after “set boxsize” in the paraset.str file.
5. Move your moleculename_gas.crd file into the /vac folder and the /water folder, and your moleculename_wbox.crd file into the /water folder.
6. Submit the jobs by running “fep.run” in two separate folders, once in the /vac folder and once in the /water folder. Let all the jobs complete. THEN tally up the totals using “xtract.run” separately in each folder. It will create a “job.1” folder where it will store all the outputs. BE SURE TO BACK UP THE FILES before doing this, as they will be deleted. (Alternately, you can deselect this option.) It will also produce a summary of the results.
7. To compute FE, you use the “total” under “final results.” Simply subtract initial from final. That is, subtract the total in the gas phase FROM the total in the water phase. Units of energy are kcal/mol.
8. Lastly, for multiple runs, be sure to change “seed” in the file “fep.run,” otherwise you will get the same exact results.

The basic idea is that you use fep.run to submit a series of simulations of your gas and water systems separately. The file xtract.run does a WHAM analysis (see CHARMM documentation) and a thermodynamic integration analysis of the results. You subtract to find the change in FE from the gas state to the solvated state (water result – vacuum result).

On building an appropriate water box:

1. You can use step4.2_waterbox.inp in the /water folder. Edit parameters “A,” “B,” “C,” and “watboxZ” in the file step3_size.str. Make all of these the same length, the length of your box in angstroms. Make sure your box is big enough to allow for a 7-or-so-angstrom buffer on all sides of your molecule, which should be placed centered at zero.
2. In the command line of DT2, you can just use CHARMM to read step4.2_waterbox.inp. Assuming the path to the CHARMM executable is “charmmlink”: charmmlink < step4.2_waterbox.inp > step4.2_waterbox.out
3. Once you have the resulting files, equilibrate the water box using step6.1_equilibration.inp: charmmlink < step6.1_equilibration.inp > step6.1_equilibration.out Feel free to adjust the dynamics time if you desire. You might want to create a shell script to submit this using sbatch, as it won’t run as fast as the other steps I’m listing here.
4. Place your molecule into the water box using step5_assembly.inp. You’ll need to change the name of your gas molecule. Mine was DMOE, so you can just use edit/replace. NOTE that this will write the file moleculename_wbox.crd, so if you already have such a file and don’t want it deleted, you’ll want to change this.