.. tip:: All input files can be downloaded: :download:`Files `. .. tip:: For more information of this section, please refer to these pages: - :doc:`../keywords/charmm` - :doc:`../keywords/mol` For generating topology of organic molecules and biomolecules, please refer to **pdbtop** (``_). Molecular Mechanics Calculations ========================================= This tutorial will lead you step by step to do molecular mechanics (MM) calculations using Qbics. In Qbics, MM means CHARMM potential. For geometry optimization and molecular dynamics, please refer to :doc:`opt` and :doc:`md`. Topology and Parameters ---------------------------------- An MM calculation needs more steps since you must set up **topology** and **parameters** for the molecule. - **Topology** defines the bonds, angles, dihedrals, impropers and atomic charges. Qbics uses **PSF** file to record topology. - **Parameters** are those for bonding and nonbonding energy terms. Qbics uses **CHARMM** format parameter files. .. tip:: For CHARMM files, ``par*`` usually contains only parameters, ``top*`` usually contains only topologies, but ``toppar*`` contains both parameters and topologies. Qbics often fails to parse ``toppar*``, so you have to delete topology information in ``toppar*`` and keep only parameter information. Building topology and parameters of a system for MM calculations is probably the most difficult and error-prone step. There are many different ways to build topology for the same system: - Manually; - (**Strongly recommended**) Using **pdbtop** (``_); - Using **ABCluster** (``_); - Using **VMD** (``_); - Using **CHARMM-GUI** (``_). We strongly recommend to use `pdbtop `_ since it is highly compatible with Qbics. pdbtop cna be used to add hydrogens, add counterions, solvate in box or droplet, and treat covalently bonded ligands. To use pdbtop, please refer: - Download pdbtop: ``_ - pdbtop manual: ``_ In the following, we will give 2 examples. Both will be used in the QM/MM tutorial :doc:`qmmm`. Solvated Protein 5VBM ---------------------------------- Here, we want to do a calculation with a protein, whose PDB ID is 5VBM. To convert a protein from PDB bank to a format that Qbics can read, please follow the procedure given here: ``_. Finally, we get a protein solvated in water, whose coordinates and topology are stored in ``5vbm-sol.pdb`` and ``5vbm-sol.psf``, respectively. You can use `Qbics-MolStar `_ to visualize it, note that load topology ``5vbm-sol.psf`` in :guilabel:`Model` and coordinates ``5vbm-sol.pdb`` in :guilabel:`Coordinates`: .. figure:: figs/a6.png Now, we need force field parameters. You can go to CHARMM website ``_ to download force field parameters, just select ``toppar_c36_jul24.tgz``, download and decompress it. Put ``par_all36m_prot.prm`` (parameters for protein) and ``toppar_water_ions.str`` (parameters for water and ions) to your current path. Now prepare the input file: .. code-block:: :caption: 5vbm-sol.inp :linenos: 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 5vbm-sol.psf scaling14 1.0 rcutoff 12.0 rswitch 10.0 use_PBC # Add this line to turn on PBC. Lbox 70 70 70 # Define the box size. electrostatic pme # Use PME for electrostatics. PMEk 72 72 72 end mol 5vbm-sol.pdb end task energy charmm end Now run the calculation: .. code-block:: bash $ qbics 5vbm-sol.inp -n 4 > 5vbm-sol.out Here, ``-n 4`` means Qbics will use 4 CPU cores for parallization. This calculation is very fast. After calculation, check ``5vbm-sol.out``: .. code-block:: :caption: 5vbm-sol.out :linenos: CHARMM energies =============== Bond energy: 7661.69410280 kcal/mol Angle energy: 5772.78601511 kcal/mol Dihedral energy: 1557.76413669 kcal/mol Improper energy: 30.51152972 kcal/mol Coulombic energy: -2120.48409716 kcal/mol Lennard-Jones energy: 11507712434764.78906250 kcal/mol --------------------------------------------------- CHARMM energy: 11507712447667.06054688 kcal/mol This energy is very large, since the water box are generated randomly. You must optimize structure for any useful calculations. This will be introduced in :doc:`opt`. Now we will explain the input ``5vbm-sol.inp`` of the last example. - Anything after ``#`` is treated as comments. You can write anything anywhere in the input file. - ``charmm ... end`` This block defines the CHARMM force field. * ``parameters`` The CHARMM parameter files. You can provide many ones. * ``topology`` The topology of the system in PSF format. * ``scaling14`` The scaling factor for 1-4 interactions. For CHARMM calculations, ``1.0`` is preferred. * ``rcutoff`` and ``rswitch`` The cutoff and switch distance (in Angstrom) for nonbonded interactions. * For van der Waals interactions, they are the switch and cutoff distances for the switching function. * For electrostatic interaction energies: when ``electrostatic`` is ``cutoff``, they are the switch and cuttoff distances for its switching function; when ``electrostatic`` is ``pme``, ``rcutoff`` is the cutoff distance for short-range part, and ``rswitch`` is not used. * ``use_pbc`` The periodic boundary condition (PBC) is used. If it is not present, a gas phase calculation is applied. * ``Lbox`` The 3 lengths (in Angstrom) of the simulation box. Only valid when ``use_pbc`` is applied. * ``electrostatic`` The calculation method of electrostatic energy, which can be cut off at a certain distance, or particle mesh Ewald (PME) algorithm. For gas phase calculations, only ``cutoff`` is allowed; when PBC is turned on, one can use either ``pme`` (you should almost **ALWAY** use this for pure MM calculations) or ``cutoff`` (only used in QM/MM calculations). * ``PMEk`` The number of K points on 3 lengths of the box for PME calculations. It is recommended to be an integer that is **close** to the box length and is factorized to **only 2, 3, and 5.** - ``mol ... end`` This block gives the molecule coordinates. Unlike in :doc:`dft`, we just give the file name。 - ``task ... end`` This tells Qbics what it should do. ``energy charmm`` means it will use CHARMM force field to calculate energy. You can put several tasks in this block. Oragnic Molecule 92V ------------------------------- In many applications, you need to do MM calculations for small molecule in gas phase. This involves how to generate the topology and parameters for an arbitrary molecule. This can be easily done with `pdbtop `_. Here, we consider a ligand called 92V. To obtain the topology and paramters of this molecule, please follow the proceduer given here: ``_. What we need are: topology ``92v-2.psf``, parameters ``92v-1-naive.prm``, and structure ``92v.pdb``. Qbics can automatically parse PDB format. You can also use XYZ format if you prefer. To calculate its CHARMM energy in gas phase, use the following input: .. code-block:: :caption: 92v.inp :linenos: charmm parameters 92v-1-naive.prm topology 92v-2.psf scaling14 1.0 rcutoff 12.0 rswitch 100.0 # rswitch>rcutoff means no cutoff is used. end mol 92v.pdb end task energy charmm end Here we want to point out: - There is no ``use_pbc`` in ``charmm...end``, suggesting that this is a gas phase calculation. - Here, ``rswitch`` is greater than ``rcutoff``, this means **NO** cutoff is applied: all interactions are calculated! After calculation, we can find energies: .. code-block:: :caption: 92v.out :linenos: CHARMM energies =============== Bond energy: 79.40377426 kcal/mol Angle energy: 12.30743545 kcal/mol Dihedral energy: 93.12335114 kcal/mol Improper energy: 0.00000000 kcal/mol Coulombic energy: -116.19015178 kcal/mol Lennard-Jones energy: 31.73702659 kcal/mol --------------------------------------------------- CHARMM energy: 100.38143566 kcal/mol Final total energy: 100.38143566 kcal/mol