QM/MM Boundary Potentials

Two boundary potential methods are available: SMBP and GSBP.

QM/MM/SMBP calculations

The solvated macromolecule boundary potential (SMBP) extends QM/MM to a three layer method. Details of this method are explained in J. Chem. Theory Comput. 2009, 5, 3114-3128, a brief description is given below:

In the SMBP, the macromolecule is separated into an inner region which is modeled explicitly and an outer region that is mimicked by the boundary potential. Since only a small subset of all atoms are modeled explicitly, the SMBP can reduce computational costs. Preliminary evaluations estimate savings on the order of 30-50% for a single energy and gradient evaluation. In the SMBP, the outer solvent region is described by a polarisable dielectric continuum (PDC) and the electrostatic potential induced by the outer macromolecule and solvent region is computed as a finite-difference solution to the Poisson-Boltzmann (PB) equation. Hence, the SMBP needs the following information from the input:

• Definition of the inner region, provided as a Tcl List of the atom numbers of the atoms in the inner region via the inner_region= keyword
• Definition of the outer region, provided as a Tcl List of the atom numbers of the atoms in the outer region via the outer_region= keyword
• Position of the dielectric boundary which is defined by the envelope of the atomic van der Waals radii. Therefore, this information is passed to the programm as a Tcl List of the atomic van der Waals radii using the vdw_radii= keyword. A list of van der Waals radii with values from the CHARMM22 force field can be generated from a list of atom types using the pb_types2vdwradii procedure.
  set rad [ pb_types2vdwradii atom_types= $types ]
• Centre of the inner region via the centre= or central_atom= keyword
• The inner region must have a spherical shape and is comprised by a virtual sphere. Since all atoms inside the sphere are modeled explicitly the dielectric constant inside the sphere is set to 1.0. Therefore, the sphere contributes to the definition of the dielectric boundary. The radius of the sphere must be such that no active atom (surrounded by its atomic van der Waals sphere) touches the dielectric boundary. Therefore, the radius of the sphere must be 2-3 Å larger than necessary (for details see JCTC paper). E.g., if all atoms within 17 Å of the centre belong to the inner region, the radius of the sphere should be set to 19 or 20 Å. The value is provided via the inner_radius keyword.
• The charge of the QM subsystem via the qm_charge= keyword. This information is necessary for the computation of the ESP charges.

The input to these keywords defines the model system and is therefore specific to each system. The boundary separating inner and outer region has to be chosen sensibly.

Moreover, the SMBP takes a few input values that are not specific to the system. Therefore, reliable default values have been determined so that it is not necessary (but possible) to change these values:

• Mesh size of the grid used for finite-difference solution of the PB equation via the grid_size keyword. The default value is 2.362 au (1.25 Å).
• Usually a focusing approach is used for solving the PB equation. Then, the mesh size of the inner grid focusing on the inner region has to be determined via the inner_grid_size keyword. The default value of 0.472 au (0.25 Å) is a careful choice and yields very high accuracy. Larger values up to 1.134 au (0.6 Å) can be chosen with acceptable loss of accuracy to reduce computational costs. The focusing approach can only be switched off by chosing inner_grid_size= undefined.
• Number of virtual surface charges used to represent the boundary potential in the QM calculations via the nrpc keyword. The default value is 90. Values should typically be in the range from 60 to 200.
• The side length of the grids in the PB solution are determined automatically such that the full systems or the inner region fit into the grid. The radius_extension= keyword can be used to increase the size of the side lengths. The default value of 5 bohr should be increased if the program shows lack of stability.
• A distribution function is used to interpolate between the grid points. The spline_order= keyword determines the shape of the function. For geometry optimisations, 3rd order B-splines (spline_order= 3) are recommended and are the default for all calculations. The trilinear distribution function (spline_order= 2) cannot be recommended.
• The dielectric constant of the outer macromolecule region is by default 1.0 (eps_mol= 1.0).
• The dielectric constant of the outer solvent region is by default 80.0 (eps_solv= 80.0) corresponding to a calculation in bulk water. Under vacuum conditions (eps_mol= 1.0 and eps_solv= 1.0), QM/MM/SMBP reproduces standard QM/MM results with reduced computational costs.

An example of SMBP input can be found in the SMBP tutorial.

Summary of SMBP keywords

QM/MM/GSBP calculations

The theory behind the generalized solvent boundary potential (GSBP) is very similar to the SMBP and is explained in detail in J. Chem. Theory Comput. 2008, 4, 1600-1609. In contrast to the SMBP, application of the GSBP is only sensible for MD simulations. Therefore, it has only been implemented in ChemShell for those QM methods that can be used for MD simulations. The GSBP can be used with mndo and pointcharge as qm_theory=. However, only the most recent versions of MNDO can be used which have been modified for application with the GSBP.

Compared to the SMBP, setup of QM/MM/GSBP calculations require more manual preparation. In QM/MM/GSBP, it is necessary to cut down the system size to achieve all computational savings that the GSBP offers. The gsbp.prepare procedure helps to remove all outer region atoms from the coordinate fragment. Moreover, it also adapts the lists of charges, atom types, etc. that the hybrid module requires.

Summary of gsbp.prepare keywords

After preparation, the GSBP can be invoked inside the hybrid module by the gsbp= keyword. Since the GSBP is conceptually very similar to the SMBP, most of the keywords are identical. However, the GSBP uses a Green's function approach in combination with a basis set expansion to circumvent repeated solution of the PB equation during the MD (see original GSBP paper in J. Chem. Phys. 2001, 114, 2914-2937). Therefore, the GSBP requires some additional information from the input:

• Size of the basis set used to represent the inner region charge distribution. It is provided as the highest multipole moment that is included via the nr_multipoles= keyword. A default value of 19 requires application of 400 basis functions.
• The coarsening of the inner region (CIR) factor determines the mesh size of the grid that is used to set the boundary values in computation of the reaction field matrix via the cir_factor= keyword. A default value of 2.5 is recommended.
• Different boundary conditions can be used during computation of the reaction field matrix and are invoked by the bc= keyword: Debye-Hückel (dh), Debye-Hückel with linear interpolation (dhli, default), and zero boundary conditions (zero).
• Basis set projection transforms the Green's function into a matrix. This reaction field matrix is computed at the beginning of each QM/MM/GSBP calculation and stored in the matrix object M_rf. If a QM/MM/GSBP calculation is restarted, computation of the reaction field matrix can be skipped by setting restart= yes. Then the reaction field matrix is read from M_rf. An existing reaction field matrix may only be used if the radius and the centre of the inner region has not changed and if the same GSBP basis set is used (same number of multipole moments).

The QM/MM/GSBP calculation is performed on the reduced system that was prepared with gsbp.prepare. However, for computation of the reaction field matrix and the electrostatic potential that is induced by the outer macromolecule region, the GSBP needs information about the full system:

• The coordinates of the full system, provided by the tot_coords= keyword.
• The atomic charges of the full sytem, provided via the tot_charges= keyword.
• The atomic van der Waals radii of all atoms of the full system, provided via the vdw_radii= keyword.

Summary of GSBP keywords

Example for QM/MM/GSBP input

The full system is first cut down to the inner region that is simulated explicitly. The coordinates of the full system are stored in tcoords.

gsbp.prepare coords= ${sys_name_id}.c  \
        groups = $groups \
        atom_charges= $charges \
        qm_region= $qmatoms \
        inner_region= $inner_mm_numbers \
        outer_region= $outer_mm_numbers \
        frozen= $frozen_atom_numbers \
        atom_types= $types \
        active= $new_active_atoms \
        tot_coords= tcoords

If Charmm is used as MM method, one has to create a .psf file of the reduced system. There is no procedure for this task in ChemShell but it can be realized with Charmm. In this example, the .psf file is saved as ${sys_name_id}_in.psf.

Finally, the dynamics module can be initialised with the QM/MM/GSBP theory:

set rad [ pb_types2vdwradii atom_types= $types ]
...
dynamics dyn1 coords= ${sys_name_id}.c  \
        frozen= $prep_frozen \
	...
	theory=hybrid : [ list \
		...
                groups= $prep_groups \
                atom_charges= $prep_charges \
                qm_region= $prep_qm_region \
                gsbp= yes : [ list \
                        tot_coords= tcoords \
                        tot_charges= $tot_charges \
                        inner_region= $inner_atom_numbers \
                        outer_region= $outer_atom_numbers \
                        centre= $centre \
                        inner_radius= 39.683 \
                        vdw_radii= $rad \
                        grid_size= 2.362 \
                        inner_grid_size= 0.472 \
                        radius_extension= 10.0 \
                        nr_multipoles= 1 \
                        spline_order= 3 \
                        eps_mol= 1.0 \
                        eps_solv= 80.0 ] \
                qm_theory= $qm_theory : [ list \
			...
			\ ]
                mm_theory= dl_poly : [ list \
			frozen= $prep_frozen \
			...
			charmm_psf_file= ${sys_name_id}_in.psf \ ] ]

For further discussion of the GSBP method please see the GSBP tutorial.