The truth is rarely pure and never simple

embedding molecules or: g_membed for namd

Embedding one molecule (e.g., a protein) into another one (e.g., a membrane) has to be done carefully in order to avoid any conformational clashes. For the GROMACS world, there is g_membed, but for namd there is no tool that I know of. So here is a TCL script to help you out. The main idea of which is to grow forbidden areas around the atoms of the molecule to be inserted while doing a normal MD simulation.


the idea of the script is to grow spheres around the atom positions of the molecule to be inserted as to model the van der Waals radii. In general, this works in the following steps:
1. Equilibrate everything.
2. Simulation of the larger system with the tcl script.
3. Manually combine the two systems.
4. Equilibrate again.

Lets suppose you have a large system (A) and a small system (B) to be embedded into A. Now you have an equilibration run for A in namd. Below, there are two separate listings of code. Copy the second one and save it as “push.tcl” in the directory of the namd.conf file, then add the content of the first listing to your namd.conf. In case you already have a tclBCScript block in the file, you have to merge them instead of adding the whole block. Extract a snapshot of B as PDB file and save it in the same directory as namd.conf (adjust the filename in the tclBCScript block as needed). If the center of mass position of B is different from the insertion location in the coordinate system from A, you will have to set the “pushvector” setting accordingly. The “radius” setting should match the typical van der Waals radius of the atoms of B as to minimise both the distortion to A and the equilibration time afterwards.

The settings for “increment” and “force” are to be found empirically. If your simulation is unstable, then probably the energy added due to the arbitrary forces is too high to dissipate within reasonable time. Decreasing the “increment” parameter is safe – it just takes longer to form the holes. The force parameter depends mostly on the cohesive forces in system A. If they are too high, then “force” has to be increased such that the force is strong enough to locally push the atoms of system A apart. On the other hand, if this parameter is too high, the additional forces imposed on the atoms will rip A apart or create unreasonable high local distortions. The default values on the website are found to work well for DPPC membranes. If you find other values to work better for different systems, I would be interested in this, as well.

In your simulation, you must not include the molecule that is to be embedded. If you do this, the spheres are grown at the same coordinates as the ones of the atoms in the combined system that belong to the small (to be embedded) molecules. This gives very small distances, and, due to the nature of the Lennard-Jones potential that is used for van der Waals repulsion, very high forces which – in turn – result in high velocities.


  1. You have a running (and hopefully equilibrated) simulation of the larger molecule.
  2. You have a PDB file of the smaller molecule.

What to do

Append the following section to your namd.conf

tclBC 			on
tclBCScript		{ 
	set place		test.pdb	;# PDB file format
	set radius		5.0 		;# maximum radius for shells in Å
	set increment	0.01 		;# growth rate of shells in Å
	set force		10.0 		;# force to apply per atom in kcal/(mol*Å)
	set refresh		50			;# number of time steps a exclusion list is valid
	set pushvector	{0 0 0}		;# all atoms in the $place file are shifted by this vector in Å

	source push.tcl

Adjust the settings as described and put the following file into the same directory as namd.conf:

set pushlist ""
set valid_keywords {"ATOM  " "HETATM"}

print "Reading $place as molecule to make room for."
set pdbfile [open $place r]
foreach line [split [read $pdbfile] \n] {
  # grep for atom definitions
  set keyword [string range $line 0 5]
  if { [lsearch $valid_keywords $keyword] >= 0 } {
    set atomvector [vecadd $pushvector [string range $line 30 53]]
    lappend pushlist $atomvector
print "Found [llength $pushlist] atoms."

# caps the growing spheres at a certain $radius
proc calc_radius { step increment radius } {
  set R [expr $increment * $step]
  if { $R > $radius } { set R $radius }
  return $R

proc calcforces { step unique } {
  # import settings
  global increment radius refresh pushlist n_atom force
  # calculate current radius of the spheres (identical for all atoms)
  set R [calc_radius $step $increment $radius]

  # refresh exclusion: only consider those atoms that currently are
  # within the sphere that will active in $refresh steps from now
  # caveat: all $refresh steps are not applying any forces as there
  # is no way to reset the nextatom iterator state from the TCL interface

  if { $step % $refresh == 0 } {

    set refresh_R [calc_radius [expr $step + $refresh] $increment $radius]  
    print "Refreshing exclusion list at step $step with prospective radius $refresh_R."

    set affected 0
    while {[nextatom]} {
      set rvec [getcoord]
      set thisaffected 0
      foreach pushatom $pushlist {
        set distance [veclength [vecsub $pushatom $rvec]]
        if { $distance < $refresh_R } {
          set thisaffected 1
      if { $thisaffected } {
        set affected [expr $affected + 1]
      } else {

    print "Up to $affected atoms will be affected by additional forces."
  } else {    
    while {[nextatom]} {
      set rvec [getcoord]
      set shiftvector {0 0 0}
      foreach pushatom $pushlist {
        set diffvector [vecsub $rvec $pushatom]
        set distance [veclength $diffvector]

        if { $distance < $R } {
          set diffvector [vecscale $diffvector [expr 1 / $distance]]
          set shiftvector [vecadd $shiftvector $diffvector]

      set length [veclength $shiftvector]
      if { $length > 0 } {
        set shiftvector [vecscale $shiftvector [expr $force / $length]]
        addforce $shiftvector
        addenergy 1

Done. Now start your simulation and monitor the movements of the atoms.

Leave a comment

Your email address will not be published.