09 Mar 2023

Prepare MD simulations in NAMD

We show how to prepare MD simulations for alanine dipeptide and chignolin with NAMD.

To do this, one typically builds a PSF file first starting from a PDB file. While a PDB file contains basic information about the atoms in a molecular system and can often be downloaded from the RCSB Protein Data Bank, a PSF file contains all necessary information (such as force field) to run MD simulations with NAMD. In this post, we show how to use psfgen plugin to generate the PSF file.

Alanine dipeptide.

1.1. PDB file

ad.pdb

ATOM      1 HH31 ACE     1       2.000   1.000  -0.000  1.00  0.00
ATOM      2  CH3 ACE     1       2.000   2.090   0.000  1.00  0.00
ATOM      3 HH32 ACE     1       1.486   2.454   0.890  1.00  0.00
ATOM      4 HH33 ACE     1       1.486   2.454  -0.890  1.00  0.00
ATOM      5  C   ACE     1       3.427   2.641  -0.000  1.00  0.00
ATOM      6  O   ACE     1       4.391   1.877  -0.000  1.00  0.00
ATOM      7  N   ALA     1       3.555   3.970  -0.000  1.00  0.00
ATOM      8  H   ALA     1       2.733   4.556  -0.000  1.00  0.00
ATOM      9  CA  ALA     1       4.853   4.614  -0.000  1.00  0.00
ATOM     10  HA  ALA     1       5.408   4.316   0.890  1.00  0.00
ATOM     11  CB  ALA     1       5.661   4.221  -1.232  1.00  0.00
ATOM     12  HB1 ALA     1       5.123   4.521  -2.131  1.00  0.00
ATOM     13  HB2 ALA     1       6.630   4.719  -1.206  1.00  0.00
ATOM     14  HB3 ALA     1       5.809   3.141  -1.241  1.00  0.00
ATOM     15  C   ALA     1       4.713   6.129   0.000  1.00  0.00
ATOM     16  O   ALA     1       3.601   6.653   0.000  1.00  0.00
ATOM     17  N   NME     1       5.846   6.835   0.000  1.00  0.00
ATOM     18  H   NME     1       6.737   6.359  -0.000  1.00  0.00
ATOM     19  CH3 NME     1       5.846   8.284   0.000  1.00  0.00
ATOM     20 HH31 NME     1       4.819   8.648   0.000  1.00  0.00
ATOM     21 HH32 NME     1       6.360   8.648   0.890  1.00  0.00
ATOM     22 HH33 NME     1       6.360   8.648  -0.890  1.00  0.00
TER
END

1.2. tcl script

See the website of psfgen for detailed explanation of the following tcl script.

One needs to modify and use the actual path to CHARMM36 force field.

build.tcl:

set segname AD

set pdbfile [lindex $argv 0] 

puts "pdb file: $pdbfile"

set name [file rootname [file tail $pdbfile]]

set sys [lindex $argv 1]

puts "Build system in $sys."

package require psfgen

# Read topology file
topology path/to/toppar/top_all36_prot.rtf

# Change names
pdbalias atom ACE CH3 CL
pdbalias atom ACE HH31 HL1
pdbalias atom ACE HH32 HL2
pdbalias atom ACE HH33 HL3
pdbalias atom ACE C CLP
pdbalias atom ACE O OL

pdbalias residue ACE ALAD
pdbalias residue ALA ALAD
pdbalias residue NME ALAD

# Build protein segment
segment $segname {
pdb $pdbfile 
}

# Use coordinates
coordpdb $pdbfile $segname

# Guess unkonwns
guesscoord

# Write structure and coordinate files
writepsf ./${name}_vacuum.psf
writepdb ./${name}_vacuum.pdb

if {$sys == "water"} {

  package require solvate 

  solvate ${name}.psf ${name}.pdb -t 10 -o ${name}_water
}

exit

1.3. To build the system in vacuum, run:

vmd -dispdev text -e build.tcl -args ad.pdb vacuum

or, to build the system in water, run:

vmd -dispdev text -e build.tcl -args ad.pdb water

This will generate the pdb and psf files that can be used to run MD simulation with NAMD.

1.4. Run simulation

1.4.1. Vacuum simulations.

Some discussions can be found here:

[1] https://www.ks.uiuc.edu/Research/namd/mailing_list/namd-l.2011-2012/4142.html

[2] https://www.ks.uiuc.edu/Research/namd/mailing_list/namd-l.2005-2006/1617.html

[3] https://www.ks.uiuc.edu/Research/namd/mailing_list/namd-l.2012-2013/2825.html

1.4.2. Water

To run simulation in water, one needs to provide the information about the simulation box.

Such information can be obtained using vmd:

vmd> source get_cell.tcl 
vmd> get_cell

where the script get_cell.tcl is:

proc get_cell {{molid top}} {
 set all [atomselect $molid all]
 set minmax [measure minmax $all]
 set vec [vecsub [lindex $minmax 1] [lindex $minmax 0]]
 puts "\n==========Cell Info=================="
 puts "cellBasisVector1 [lindex $vec 0] 0 0"
 puts "cellBasisVector2 0 [lindex $vec 1] 0"
 puts "cellBasisVector3 0 0 [lindex $vec 2]"
 set center [measure center $all]
 puts "cellOrigin $center"
 puts "======================================\n"
 $all delete
}

Chignolin

2.1. PDB

We download the PDB file from RCSB Protein Data Bank with PDB code 5awl. The file is named as 5awl.pdb and is saved in the current directory.

2.2. psfgen tcl scipt

The tcl script for building PSF file is similar as above.


set segname CGL

set pdbfile [lindex $argv 0] 

puts "pdb file: $pdbfile"

set name [file rootname [file tail $pdbfile]]

set sys [lindex $argv 1]

puts "Build system in $sys."

package require psfgen

# Read topology file
topology path/to/toppar/top_all36_prot.rtf

# Build protein segment
segment $segname {
pdb 5awl.pdb
}

# Use coordinates
coordpdb 5awl.pdb $segname

# Guess unkonwns
guesscoord

# Write structure and coordinate files
writepsf ./${name}_vacuum.psf
writepdb ./${name}_vacuum.pdb

if {$sys == "water"} {

  package require solvate 
  solvate ${name}_vacuum.psf ${name}_vacuum.pdb -t 12 -o ${name}_water

  package require autoionize
  autoionize -psf ${name}_water.psf -pdb ${name}_water.pdb -neutralize -o ${name}_ionized
}

exit

The ways to generate PDB/PSF files for NAMD and to run MD simulation are the same as the ways described above for alanine dipeptide.