charmm

This keyword defines the implementation details of CHARMM force field.

Options

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.

charmm
  use_PBC
  Lbox 86 86 86
end

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{aligned} 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{aligned}$$

$$\begin{aligned} 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{aligned}$$

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:

1. Atoms: Lists the ID, residue information, atom type, charge, mass, etc. for each atom.
2. Bonds: Defines the covalent bonding connections. Each line contains the IDs of pairs of atoms.
3. Angles: Records bond angles formed by three atoms. Each row contains three sets of IDs.
4. Dihedrals: Defines dihedral (torsion) angles formed by four atoms. Each row contains four atomic IDs.
5. 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:

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

2. 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:

    7517 !NBOND: bonds
 1         2         1         6         1         7         2         3 
 2        12         3         4         3         8         4         5 
 4         9         5         6         5        10         6        11 
 7        13         7        14        12        13        14        15
14        16        17        18        17        19        17        20 
17        21        19        22        23        24        23        25 
23        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:

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

Other sections like NBFIX is NOT used by Qbics.

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
    # 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
    # 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