Tip

All input files can be downloaded: Files.

charmm

This keyword defines the implementation details of CHARMM force field.

Options 

parameters

Value

One or more file names

Default

None

Give the CHARMM force field parameter file names (paths may be included). There is no default one so at least one file name must be given. There can be more than one files:

 charmm
   parameters par_all35_ethers.prm /home/zhang/userdef/par_all35_proteins.prm
 end

Note that if the same term, say bond CT-CT, appears in more than one parameter files, then the latter one will overwrite the previous ones.

topology

Value	A file name
Default	None

Give the topology file name of the molecule in PSF format (path may be included).

scaling14

Value	A real number
Default	`1.0`

Define the scaling factor for electrostatic interactions of 1-4 term. For most cases, it is set to 1.0.

rcutoff

Value	A real number
Default	`15.0`

Define the distance in Angstrom to cutoff non-bonded interactions.

When use_PBC is turned on, rcutoff should be smaller than the half of periodic box size to avoid self interactions.

If rswitch is larger than rcutoff, switching function will not be used and no cutoff is applied on long range forces.

rswitch

Value	A real number
Default	`10.0`

Define the distance in Angstrom to activate switching function.

If rswitch is larger than rcutoff, switching function will not be used and no cutoff is applied on long range forces.

use_PBC: Turn on periodic boundary condition for the system. If this option is not present, the system will be a non-periodic, gas phase one.

Lbox

Value	Three real numbers
Default	`0. 0. 0.`

When use_PBC is present, this option gives the lattice lengths along X, Y, and Z direction in Angstrom.

charmm
  use_PBC
  Lbox 86 86 86
end

This creates a simulation box of 86 Angstrom x 86 Angstrom x 86 Angstrom.

electrostatic

Value	`Cutoff` Cutoff scheme for gas phase or periodic systems.
	`PME` Particle mesh Ewald (PME) method. Only available for periodic systems.
Default	`Cutoff`

Assign the scheme to calculate electrostatic interaction.

For gas phase, Cutoff is the only possible way.

For MM calculations in periodic systems, both are available and it is recommended to use PME.

For QM/MM calculations in gas phase or periodic systems, Cutoff is currently the only possible way.

PMEk

Value	Three integers
Default	`64 64 64`

Assign the number of grids for PME long range force calculations. It is recommended that the number should be chosen as an integer close to the lattice length while is a multiply of only 2,3 and/or 5. For example, if Lbox is 109. 109. 77, then PMEk can be set to 108 108 80.

Theoretical Background 

CHARMM Force Field 

To do a CHARMM force field calculation, besides system coordinates, one needs to provide the system topology (by topology) and force field parameter files (by parameters). The CHARMM force field is decomposed into bonded interactions and nonbonded interactions. The bonded interactions include bond stretching, angle bending, dihedral angle rotation, and improper torsion. The nonbonded interactions include Lennard-Jones interactions and electrostatic interactions. They are all defined in topology and parameters.

For nonbonded interactions, one can calculate all pair interactions (by setting rswitch larger than rcutoff) for non-periodic systems. However,sometimes it is computationally expensive to calculate all pair interactions. In this case (for both non-periodic and periodic systems), one can use a cutoff scheme to truncate the long-range interactions. In the cutoff scheme (electrostatic cutoff), the interactions are truncated at a certain distance rcutoff and a switching function is used to smoothly turn off the interactions at a distance of rswitch. Such switching functions are different for Lennard-Jones and electrostatic interactions. The expressions are shown below:

\[\begin{split}S_{\text{LJ}}(r) = \begin{cases} 1 & \text{if } r \le r_\text{s} \\ \frac{\left(r_\text{c}^{2}-r^{2}\right)^2\left(r_\text{c}^{2}-r^{2}-3\left(r_\text{s}^{2}-r^{2}\right)\right)}{\left(r_\text{c}^{2}-r_\text{s}^{2}\right)^3} & \text{if } r_\text{s} \le r \le r_\text{c} \\ 0 & \text{if } r \ge r_\text{c} \\ \end{cases}\end{split}\]

\[\begin{split}S_{\text{el}}(r) = \begin{cases} \left(1-\frac{r^2}{r_\text{c}^{2}}\right)^2 & \text{if } r \le r_\text{c}\\ 0 & \text{if } r \le r_\text{c} \\ \end{cases}\end{split}\]

The plots of these functions are shown below:

For periodic systems, one can use the particle mesh Ewald (PME) (electrostatic pme) method to calculate long-range electrostatic interactions. Also, rcutoff should be smaller than the half of periodic box size to avoid self interactions.

PSF File Format 

A protein structure file (PSF) is a file format utilized to describe the topological information of biological macromolecules. PSF files primarily store information regarding atoms, bonds, angles, dihedral angles, and improper force terms. However, they do not contain atomic coordinates intrinsically.

The PSF file consists of the following sections related to Qbics:

Atoms: Lists the ID, residue information, atom type, charge, mass, etc. for each atom.
Bonds: Defines the covalent bonding connections. Each line contains the IDs of pairs of atoms.
Angles: Records bond angles formed by three atoms. Each row contains three sets of IDs.
Dihedrals: Defines dihedral (torsion) angles formed by four atoms. Each row contains four atomic IDs.
Impropers: Mechanical constraints used to maintain molecular planes (e.g., the coplanar structure of the benzene ring), four IDs per row.

      6 !NTITLE
REMARKS original generated structure x-plor psf file
REMARKS 2 patches were applied to the molecule.
REMARKS topology top_all27_prot_lipid.inp
REMARKS segment U { first NTER; last CTER; auto angles dihedrals }
REMARKS defaultpatch NTER U:1
REMARKS defaultpatch CTER U:76

   1231 !NATOM
      1 U    1    MET  N    NH3   -0.300000       14.0070           0
      2 U    1    MET  HT1  HC     0.330000        1.0080           0
      3 U    1    MET  HT2  HC     0.330000        1.0080           0
      4 U    1    MET  HT3  HC     0.330000        1.0080           0
      5 U    1    MET  CA   CT1    0.210000       12.0110           0
      6 U    1    MET  HA   HB     0.100000        1.0080           0
      7 U    1    MET  CB   CT2   -0.180000       12.0110           0

Key Sections of this file:

Title Section: Contains annotation information of the PSF file, such as generation method, topology file path, and so on.

Atoms Section:

Atom ID: A unique number (e.g. 1, 2).

Segment Name: e.g. U for protein chain.

Residue ID: The number of the residue to which the atom belongs (e.g. 1 for the first residue).

Residue Name: e.g. MET for methionine.

Atom Name: e.g. N, HT1, HT2, CA for nitrogen, hydrogen, carbon, etc.

Atom Type: The type used for force field calculations (e.g. NH3, HC, CT1).

Charge: Partial charge of the atom (e.g. -0.300000).

Mass: Mass of the atom (e.g. 14.0070).

Then, the Bonds, Angles, Dihedrals, and Impropers sections follow, each containing the corresponding information. For example, the Bonds section contains the IDs of atoms that are bonded together:

!NBOND: bonds
       2         1         6         1         7         2         3
      12         3         4         3         8         4         5
       9         5         6         5        10         6        11
      13         7        14        12        13        14        15
      16        17        18        17        19        17        20
      21        19        22        23        24        23        25
      26        23        27        25        28        29        30

Here, we have bonds 1-2, 1-6, 1-7, 2-3, and so on. For others, the format is similar.

CHARMM Parameter File Format 

CHARMM force field files define the rules and parameters of interactions between atoms in a molecular system. It contains several sections related to Qbics:

BONDS: Describe the equilibrium lengths of covalent bonds and their corresponding force constants.
ANGLES: Provides the equilibrium value of the angle within the molecule and its force constant.
DIHEDRALS: Define the potential energy function and the corresponding parameters of the torsion angle.
IMPROPERS: Define the potential energy function and the corresponding parameters of the out-of-plane torsion angle.
NONBONDED: Define the Lennard-Jones potential.

Other sections like NBFIX is NOT used by Qbics.

Hint

CHARMM force field files can be obtained from:

https://mackerell.umaryland.edu/charmm_ff.shtml

Input Examples 

Example: Molecular Dynamics Simulation of Ubiquitin in Water 

Below we provide an example of a Qbics input file for simulations of a protein “ubiquitin” in water. Below shows the system. It contains 24696 atoms, a ubiquitin in a 64x64x64 water box:

This system was built by randomly placing water molecules around the ubiquitin. This Usually leads to a high-energy system, and molecular dyanmics starting with this can be highly unstable. So, we first optimize it:

charmm-1a.inp

charmm
    # Here:
    # toppar_water_ions.str: the parameter file for water and ions.
    # par_all36m_prot.prm: the parameter file for proteins.
    # Both are obtained from website: https://mackerell.umaryland.edu/charmm_ff.shtml.
    parameters toppar_water_ions.str par_all36m_prot.prm
    topology   sys.psf
    scaling14  1.0
    rcutoff    12.0
    rswitch    10.0
    use_PBC         # Add this line to turn on PBC.
    Lbox 64 64 64   # Define the box size.
    electrostatic pme # Use PME for electrostatics.
    PMEk   64 64 64
end

opt
    max_step 500
end

mol
    sys.pdb
end

task
    opt charmm
end

After 500 steps, the energy decreased from 34167803.38107745 kcal/mol to -82402.25839655 kcal/mol.

Using the optimized structure, we can now perform a NVT MD simulation:

charmm-1b.inp

charmm
    # Here:
    # toppar_water_ions.str: the parameter file for water and ions.
    # par_all36m_prot.prm: the parameter file for proteins.
    # Both are obtained from website: https://mackerell.umaryland.edu/charmm_ff.shtml.
    parameters toppar_water_ions.str par_all36m_prot.prm
    topology   sys.psf
    scaling14  1.0
    rcutoff    12.0
    rswitch    10.0
    use_PBC         # Add this line to turn on PBC.
    Lbox 64 64 64   # Define the box size.
    electrostatic pme # Use PME for electrostatics.
    PMEk   64 64 64
end

md
    type        nvt
    dt          0.001 # 1 fs
    num_steps   1000  # 1 ps. Just for demonstration.
    output_freq 100
    temp        298
end

mol
    charmm-1a-opt.xyz # Use the optimized structure from the previous step.
end

task
    md charmm
end

charmm

Options

Theoretical Background

CHARMM Force Field

PSF File Format

CHARMM Parameter File Format

Input Examples

Example: Molecular Dynamics Simulation of Ubiquitin in Water