This is a tutorial on how to calculate the interaction energy and it is worth noting that there was help from other sources such as a template code on the Charmm source codes website by Lennart Nilsson.
Begin by downloading the source files inte.tar.gz
It is important to re-center the trajectory so that system is always in the center. There is a way to check if everything works by averaging together the total interaction energy, total electrostatic, and total VdW energies and seeing if the sum of the average of VdW+ average ELEC=average(ENER)
To re-center the trajectory examine the recenter_2 folder and look specifically at the recenter1.inp recenter2.inp recenter3.inp and adjust the selection to the specific need. This case is examining a dimer PROA and PROB and both are re-centered.
Run them by typing
./recenter1.scr
./recenter2.scr
./recenter3.scr
Alternatively something similar can be done by adding these lines somewhere between the traj read and the inte command lines
coor stat sele segid PROA .or. segid PROB end
coor tran xdir -?XAVE ydir -?YAVE zdir -?ZAVE
energy
Then adjust the path in the interaction energy script to access the new translated dimer system that is centered.
This script is used to calculate the per/residue interaction energies over several trajectories. Each loop that is incorporated into the code slows everything down a lot and that is why the charmm script has to be invoked individually for each residue and again for each trajectory file and then combined using linux.
The .scr file
sbatch -A energy-bio interaction-energy.scr
However, the analysis runs on each .dcd file individually in separate folders so that I automate the launch process. Adjust for different dcd file numbering by using the copy file to duplicate the folders with the appropriate numbering and desired range
Then adjust the numbering and range in the start file
To officially start the program use the following commands
chmod u+x start
./start
Each line of the program submits a different job that can be specialized for different residue numbers “RES”, and to calculate on different input trajectory (.dcd) files “START” these variable names are passed on and recognized by the input file
I setup the analysis so that all the files in the same folder correspond to a certain .dcd file number. For instance this is the first file and corresponds to dyn30.dcd, then in the folder named 31 I change START to 31 using Replace feature in notepad
$chm RES=1 START=30 < interaction-energy.inp >& inte_1.out
$chm RES=2 START=30 < interaction-energy.inp >& inte_2.out
The .inp file
The input file is what CHARMM reads and tells CHARMM what to do
dimension maxres 150000 needs to be at the top of the input file to override the limits for large simulations
The lines below are standard to load the files of the system and set the number of frames in the .dcd file adjust as needed. It can be found by loading the .psf and then .dcd file into vmd and seeing the number of frames in one single .dcd file (note that it starts at 0 so add one) if it is different from 200 then the time tracker needs to be recalibrated as well at the bottom of the input file
Copy the updated parameter files to the local directory and update the .psf and .crd files by copying replacing them in the local directory or adjust the location in the psfcrd, update the cryst.str file with the updated dimensions of the simulation based off the last simulation in the file dyn???.xsc.gz using the reference positions of the three bold numbers below as the x, y, and z dimensions
5000 101.105358464 0 0 0 101.105358464 0 0 0 102.441023151 0 0 0 1.81888398977e-05 1.81888398977e-05 4.25894846733e-05 0 0 0
stream toppar.str
stream psfcrd.str
stream cryst.str
calc maxt = 200
Update the path to the .dcd files
open unit 1 read file name "/lustre/rallsopp/glucoM4/namd/dyn"@k".dcd"
This section below should be updated from the charmm directory, one directory above the namd directory in the step5_production.inp or the step7_production.inp, with the @fftx … being replaced with the numbers in the checkfft.str file in the same directory
fftx @fftx ffty @ffty fftz @fftz
nbonds atom vatom vfswitch bycb -
ctonnb 10.0 ctofnb 12.0 cutnb 16.0 cutim 16.0 -
inbfrq -1 imgfrq -1 wmin 1.0 cdie eps 1.0 -
ewald pmew fftx 120 ffty 120 fftz 120 kappa .34 spline order 6
energy
Finally it is time to check over the selection of the receptor and the ligand. It is crucial to look over the .pdb, .psf, .crd files and know the different call names to select the different components of the system. For instance as of 2021 a protein/peptide A shows up as “PROA” and a carbohydrate ligand chain is named “CARA”, if there are multiple it would be “CARA”, “CARB”, “PROB”
Start by selecting the first residue of the polysaccharide or peptide using the @id command to substitute in the correct number and then select the receptor, if there are multiple that are interacting it is possible to select them all with the .or. to include multiple. And it is also possible to do more specific analysis of the receptor by specifying individual residues there as well.
inte sele segid CARA .and. resid @id end sele segid PROA .or. segid PROB end
Tabulating the Data
Author: Shyam Patel
-
Next, uncompress and upload the .py file to your directory and uncompress the python file. Move it to one directory outside of the range of directories you wish to tabulate. Replace startDir and endDir with the number component of the directories you wish to run the script over. Update the dir variable with anything else that may be in your directory name. Then update the fileBase variable to the name of the files you wish to tabulate in each directory (i.e inter_NUM.dat).
After, activate your anaconda environment (if you don't have anaconda installed, go to the anaconda page for linux and download it your home directory).
conda activate NAME OF ENVIRONMENT
- If you do not have numpy or pandas installed run pip3 install numpy and pip3 install pandas (assuming you are using python3)
- then run python3 IntEnergyComplexTab.py
pip3 install numpy
pip3 install pandas
python3 IntEnergyComplexTab.py