CHARMM Development Project Forums CHARMM Community Script Archive interaction-energy.inp

*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,

Why you use a cut-offs while computing interaction energies from simulation. I realize that it is impossible to use a long-range corrections for an interactions, but should it be a cutnb of 999, for example? and why shift but not switch? Any specific reason?

I normally use the same options for the interaction energy as I used during the simulation, so that the result is the same as would have been obtained if these terms had been saved during the simulation. It is of course possible to use something different, but why complicate things - remember that you will have to write an understandable methods section at some point...
One can turn the argument around:
If you find an important interaction by using an infinite cutoff in the interaction energy calculation, then why was this not important enough to have been included in the first place?

I strongly agree on using the same interaction potential as was used during the simulation , at least for the non-"covalent " force-field terms while , for instance, the popular approach to post-process the solvent effects with a continuum approximation, may be better of using an "infinite" cutoff.

It make sence if you want to estimate site-site or let's say amino-acid - amino-acid interactions, but... for post-processing such as computation of receptor-ligand(molecule-molecule) interaction the inclusion of all atoms in the system may even change the sign of total non-bonded interactions for higly charged species.(it is the case jb007 mentioned above as an example, post-processing with some continuum electrostatic methods).

Thanks a lot for your reply.

Dear CHARMMERs,

I was using the script lennart provided here to calculate the interaction energy between two portions of the system. I wonder how "750 ps" used here comes out. Say my trajectory includes trj24..trj30, and each trajectory contains 50 ps. In this case, should I make changes to this script as follows:

========================
set t 1200.0
...
incr t by 5.0
if t lt 1500.0 goto loop
=======================

However, when I used something else such as 3000.00 instead of 1500.0 here in the IF..loop command, the calculation still continued after the 1500.0 ps, and gave non-zero interaction energy values.

I'm little confused here, could someone provide some helpful comments?

Thank you very much.

Shirley

You have 7 files, so the time interval is 350 ps, not 300.

Thank Rick very much for the helpful comment.

So, for traj 24..30, which of the following settings should I use?

========= 1 ===============
set t 1200.0 !! 1200 = 24*50
...
incr t by 5.0
if t lt 1500.0 goto loop !! 1500 = 30 * 50
=======================

========== 2 ==============
set t 0.0
...
incr t by 5.0
if t lt 350.0 goto loop
=======================

Thank you again for your kind help.

Shirley

Neither-- dyn24.trj ends with the 1200 ps point, but the first time point in that file should be 1151 ps.

Thank you for your help, Rich.

I wanted to see what would happen if I assign the last maximum time point a value exceeding the real maximum, so I did two test runs for traj24-25 which includes a total interval of 100 ps covering the time points of 1151-1250 ps. The following settings were used for the two runs respectively:

======= param for run1 =================
set tmax 1250.0 !! This value will be diff. in run2
set t 1151.0
...
incr t by 5.0
if t lt @tmax goto loop
...
=== end of param for run1 ===============

======= param for run2 =================
set tmax 1500.0 !! NOTE the diff. from the value in run1

set t 1151.0
...
incr t by 5.0
if t lt @tmax goto loop
...
=== end of param for run2 ===============

What I expect from these tests is that test 2 should give an error, since the assigned value for @tmax exceeds the maximum available trajectory time points, or at most it gives the same amount of time-energy output. However as a matter of fact, test 2 gives the output for time points from 1151.0 ps all the way to 1446.0 ps with total time interval of 300 ps. I wonder how the output for the last 200 ps was generated.

I guess there must be something wrong in my input, or I didn't understand the scheme correctly. Could you give some suggestions?

Thank you very much.

Regards,

Shirley