.. tip:: All input files can be downloaded: :download:`Files `. charmm ====== .. contents:: :local: This keyword defines the implementation details of CHARMM force field. Options ------------ .. contents:: :local: .. option:: parameters .. list-table:: :stub-columns: 1 :widths: 5 20 * - 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: .. code-block:: bash :linenos: 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. .. option:: topology .. list-table:: :stub-columns: 1 :widths: 5 20 * - Value - A file name * - Default - None Give the topology file name of the molecule in PSF format (path may be included). .. option:: scaling14 .. list-table:: :stub-columns: 1 :widths: 5 20 * - 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``. .. option:: rcutoff .. list-table:: :stub-columns: 1 :widths: 5 20 * - 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. .. option:: rswitch .. list-table:: :stub-columns: 1 :widths: 5 20 * - 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. .. option:: 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. .. option:: Lbox .. list-table:: :stub-columns: 1 :widths: 5 20 * - 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. .. code-block:: bash :linenos: charmm use_PBC Lbox 86 86 86 end This creates a simulation box of 86 Angstrom x 86 Angstrom x 86 Angstrom. .. option:: electrostatic .. list-table:: :stub-columns: 1 :widths: 5 20 * - 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. .. option:: PMEk .. list-table:: :stub-columns: 1 :widths: 5 20 * - 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: .. math:: 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} .. math:: 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} The plots of these functions are shown below: .. figure:: figs/charmm-2.png 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. .. code-block:: bash :linenos: 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: .. code-block:: bash :linenos: 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: #. ``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: .. figure:: figs/charmm-1.jpg 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: .. code-block:: bash :linenos: :caption: 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: .. code-block:: bash :linenos: :caption: 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