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.
Idea
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.
Prerequisites
- You have a running (and hopefully equilibrated) simulation of the larger molecule.
- 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 } {
cleardrops
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
break
}
}
if { $thisaffected } {
set affected [expr $affected + 1]
} else {
dropatom
}
}
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.