*FILENAME: interaction-energy.inp
*PURPOSE: compute interaction energies from trajectory
*AUTHOR: Lennart Nilsson, Karolinska Institutet, October 2003
!unix environment variable CHM_HOME points to CHARMM installation directory
read rtf card name $CHM_HOME/toppar/top_all22_prot.inp
read para card name $CHM_HOME/toppar/par_all22_prot.inp
read psf card name my_psf.psf

open unit 51 read unform name my_traj1.cor
open unit 52 read unform name my_traj2.cor
! specify how we are going to read the trajectory
traj firstu 51 nunit 2 skip 2500 ! use whole trajectory, pick frames every 5 ps

open write unit 21 form name inter_1_2.dat
write title unit 21
* time res1-res2 prot-ligand
set t 0.0 ! keep track of time, assume traj is from 0 to 750 ps

label loop
! get next coordinate set according to specifications above
traj read
! we have to update lists every time, things can move a lot in 5ps
update cutnb 12.0 ctofnb 12.0 fshift vshift cdie
! first get interaction energy between residues 1 and 2 in the protein
inte sele segid prot .and. resi 1 end sele segid prot .and. resi 2 end
! ?ener is total energy. other terms may also be extracted,
! see energy.doc, section Substitution, for more details
set e1 ?ener
! and protein ligand interaction
inte sele segid prot end sele segid lig end
set e2 ?ener
write title unit 21
* @t @e1 @e2

incr t by 5.0
if t lt 750.0 goto loop

Lennart Nilsson
Karolinska Institutet
Stockholm, Sweden