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

Note that interaction energy is an estimate of the enthalpy, and would only correlate strongly to binding if the entropy term of the free energy is small or constant. Also, the "rule of thumb" has been that it's difficult to distinguish between cases where the binding constant difference is significantly smaller than 3 orders of magnitude.

Dear all,

I would like to examine the interaction energy between peptide residue and lipid bilayer, but as a function of distance from membrane center instead of time. Could you please guide me how to get that?

Thank you very much.

Extract the distance you are interested in as well, and put that instead of the time (or as an extra column) in the output table. See eg COOR AXIS command (corman.doc).

Hi there.
I'm trying to determine the interaction energy between two lipids in a lipid bilayer. How do I determine the starting point, i.e. how do I tell CHARMM where (in the TRJ) to start the reading? I tried with BEGIN and START, but it's not working. The trajectory was save with SKIP=1000, whereas I'm reading it with SKIP=10000. The trajectory covers 30 ns of simulation, but I would like to analyse only the last 20 ns.

My input is as follows:
* measure interaction energy

bomlev -2

! parameters needed
! only parts before PSF/COR/TRJ needed
! c - composition of membrane system
! r - residue CHL1/SM
! 1 - residue number, must be a CHL1
! 2 - name of segi 2

! name of main lipid segment
set 2 main

! Load parameters
open read unit 1 card name top_all27_lipid.rtf
read rtf card unit 1
clos unit 10

open read unit 1 card name par_all27_lipid.prm
read param card unit 1

stream toppar_all27_lipid_cholesterol.str

! Load coordinates
open read unit 1 card name sys8fixed.psf
read psf card unit 1
open read unit 1 card name dyn.cor
read coor card unit 1

!defi test sele resi @1 .and. segi main show end

!stop

!coor stat

!stop

crystal define tetragonal 95.0 95.0 59.5 90.0 90.0 90.0
crystal build cutoff 14.0

imag byres sele resnam TIP3 .or. resnam CHL1 .or. resn SM end

nbond atom vatom cutnb 14.0 ctofnb 10. cdie eps 1. -
ctonnb 7. vswitch cutim 14.0 imgfrq -1 wmin 1.0 -
ewald pmew fftx 128 ffty 128 fftz 72 kappa .33 spline order 6

ener

!stop

! Load trajectory
open read unit 11 file name dyn_c-05.0-35.0ns.trj
traj firstu 11 nunit 1 skip 10000 !begin 15001000

!stop

! Open output data file
! asterisks denote which interaction we are measuring
open write unit 99 form name inter-ener-mem@c-@r-@s-@1-@2.dat
write title unit 99
* time ep evdw el
*

set t 0.0 ! keep track of time, assume traj is from 0.0 till the end

label loop

traj read
! we have to update lists every time
! use the same parameters as during simulation!
upda atom vatom cutnb 14.0 ctofnb 10. cdie eps 1. -
ctonnb 7. vswitch cutim 14.0 imgfrq -1 wmin 1.0 -
ewald pmew fftx 128 ffty 128 fftz 72 kappa .33 spline order 6

inte sele (resi @1 .and. segi main) end sele .not. (resi @1 .and. segi main) end

! ?ener is total energy. other terms may also be extracted,
! see energy.doc, section Substitution, for more details
set ep ?ener
set evdw ?vdw
set el ?elec

write title unit 99
* @t @ep @evdw @el
*

incr t by 100.0
if t le 3500.0 goto loop

stop
stop

------
Thanks in advance,
zee

Why do you need bomlev -2?
BEGIn is the correct keyword, so what was the problem you had ("not working" is a bit brief if you want help).
Interaction energies with pme are not particularly well-defined (pme is not really pairwise decomposable), and I don't think you need the images either for this kind of calculation. Your interaction_energy command does not calculate the interaction energy between two lipids.
Note that your OPEN READ and TRAJ FIRST commands do NOT load the trajectory (they do not transfer any data).

Thank you for the prompt reply.
Bomlev is set to -2 for no particular reason, I set it to 0.

Since PME was used during simulation I think it is better to leave it on for the analysis part.

I apologize for not being clear enough. I calculate the interaction energy between a given lipid molecule and the rest of the system.

I uncommented the "begin 15010000" part and re-run the job. It doesn not fail, though it should since I set the value of "t" to 3600 which is outside the trajectory (total number of frames is 35001000). I attach the output file.

What confuses me is CHARMM not giving any indication of the BEGIn being set as it does for skip. Is this normal behaviour?

zee

I'd simply not have a bomlev command at all.
As you can see from your output there is no PME and no image energy from the interaction_energy command.
I must say that I fail to see the relevance of the value of your variable t.
Yes, the TRAJ command for some reason ("historical"?) does not print out any more information than what you see. You can change this if you want by adding code as needed in file source/io/trajio.src.
Your initial energy command gives a very large (0.20334E+19) vdW energy.

Thanks for the advice on the bomlev.

I see there's no PME and no image energy for the interaction_energy but I fail to see why this happens. Though the energy.doc file suggests the values for energy will be zeroed but for the following terms:
Bond - Energy defined by the two atoms involved.
Angles - Energy allocated to the central atom (auto energy only).
Dihedral - Energy defined between central two atoms
Improper - Energy defined by first atom (auto energy only)
van der Waal - ATOM option only. Energy defined by two atoms involved.
Electrostatic - ATOM option only. Energy defined by two atoms involved.
Hbond - Energy defined by heavy atom donor and acceptor atom.
Harmonic cons - Energy allocated to central atom (auto energy only).
Dihedral cons - Energy defined by central two atoms.
User energy - Atom selections may be passed to USERE in the selection
common (DEFIne command). Fill forces and energies as desired.

Thank you for the information on the TRAJ command.

The initial energy command gives a large vdW energy because the initial (step 0) unit cell is not tall (along z axis) enough. This error is corrected in subsequent steps because the trajectory contains the crystal data from the simulation.

The variable "t" is just a measure of time and I use to step out of the "traj read" loop. I don't know else what could I use as a trigger/switch to exit the loop otherwise. Will CHARMM exit the loop when the end of the trajectory is reached if no other conditions are specified?

Thank you again.
zee

Minor addition. Upon adding the "ener" command in the loop I get:
NER ENR: Eval# ENERgy Delta-E GRMS
ENER INTERN: BONDs ANGLes UREY-b DIHEdrals IMPRopers
ENER EXTERN: VDWaals ELEC HBONds ASP USER
ENER IMAGES: IMNBvdw IMELec IMHBnd RXNField EXTElec
ENER EWALD: EWKSum EWSElf EWEXcl EWQCor EWUTil
---------- --------- --------- --------- --------- ---------
ENER> 0 -94460.30516 94407.04284 18.63891
ENER INTERN> 22555.96027 24860.01355 7553.39409 12340.67470 0.00000
ENER EXTERN> -3931.87294-114083.62288 0.00000 0.00000 0.00000
ENER IMAGES> -865.40355 -9797.38733 0.00000 0.00000 0.00000
ENER EWALD> 1773.42130-763869.88948 729004.40712 0.00000 0.00000
---------- --------- --------- --------- --------- ---------

So I think the script is working as expected.

The interaction_energy command simply does not compute the other energy terms, as clearly indicated in the documentation.

You are right in using a loop counter, but in order to interpret t=3600 as meaning 36 ns you have to specify t correctly. If your timestep during the simulation was 0.001 ps and you use SKIP 10000 then your t increment of 100 means the t=1 corresponds to 0.1ps, and t=3600 is 3.6 ns, not 36 ns. The initial t should also match the BEGIN value.
The CHARMM variable ?TIME contains the time into the trajectory (in ps) after each invocation of TRAJ READ.
CHARMM stops reading when it reaches end-of-file.