Command disabled: backlink

Calculation of interaction energies

You can do this nicely and easily using the NAMDenergy plugin of VMD.
Or, you can do it directly with NAMD as follows:

  • Make a clean directory.
  • Copy to it the solute-only DCD & PSF files from your run (my.dcd and my.psf).
  • Copy to it the parameter file (par_all27_prot_na.inp).
  • Copy to it the solute-only PSF file from PSFGEN (my.pdb).
  • Edit the my.pdb file and make sure the occupancies the B-factor columns are both 1.0. Your file should look like this
REMARK original generated coordinate pdb file
ATOM      1  N   ARG     1      -9.063   0.480  -0.011  1.00  1.00      A   
ATOM      2  HT1 ARG     1      -9.839  -0.024   0.368  1.00  1.00      A   
ATOM      3  HT2 ARG     1      -9.096   0.458  -1.010  1.00  1.00      A   
ATOM      4  HT3 ARG     1      -9.082   1.427   0.310  1.00  1.00      A   
ATOM      5  CA  ARG     1      -7.825  -0.148   0.437  1.00  1.00      A   
....
  • Generate a script file for NAMD:
structure                my.psf
paraTypeCharmm           on
parameters               par_all27_prot_na.inp
outputname               junk
numsteps                 1
exclude                  scaled1-4
temperature              0
COMmotion                yes
cutoff                   16         # <====
switchdist               14         # <====
dielectric               1.0        # <====
pairInteraction          on
pairInteractionGroup1    1
pairInteractionFile      my.pdb
pairInteractionSelf      on
coordinates              my.pdb
set ts 0
coorfile open dcd my.dcd
while { ![coorfile read] } {
   firstTimestep $ts
   run 0
   incr ts 1
}
coorfile close
  • Run it with something like /usr/local/namd/namd2 script | grep '^ENERGY:' > LOG &
  • awk '{print $12}' LOG | plot to plot TOTAL
  • awk '{print $7}' LOG | plot to plot ELECT
  • awk '{print $8}' LOG | plot to plot VDW
  • awk '{print $7+$8}' LOG | plot to plot non-bonded.
  • etc (the columns in LOG are the same as always from NAMD, ie. ETITLE: TS BOND ANGLE DIHED IMPRP ELECT VDW BOUNDARY MISC KINETIC TOTAL TEMP TOTAL2 TOTAL3 TEMPAVG).

If you want to calculate interaction energies between two species, you need to modify the B-factor column of my.pdb so that the first species will have the number 1.00, the second the number 2.00. In addition you need to tell NAMD what you have done by adding a couple of lines in the script. On the whole, I think it will be safer (and surely easier) to avoid all this by using VMD's NAMDenergy plugin for this type of calculation.

Take a good look at the dielectric constant used in the NAMD script. The value of 1.0 (corresponding to vacuum) immediately tells you that you should consider the magnitude of the electrostatics contribution to the total energy a very rough (if not very wrong) approximation.

research/howto/calculation_of_interaction_energies_using_namd.txt · Last modified: 2009/02/24 19:41 (external edit)