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.

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

  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

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

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

Leave a comment

Your email address will not be published.