Tip

All input files can be downloaded: Files.

Tip

For more information of this section, please refer to these pages:

For generating topology of organic molecules and biomolecules, please refer to pdbtop (http://zhjun-sci.com/pdbtop.html).

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 Geometry Optimization and Standard Molecular Dynamics Simulations.

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:

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:

In the following, we will give 2 examples. Both will be used in the QM/MM tutorial QM/MM Simulations.

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: http://zhjun-sci.com/pdbtop/doc/tutorial-1.html. 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 Model and coordinates 5vbm-sol.pdb in Coordinates:

../_images/a6.png

Now, we need force field parameters. You can go to CHARMM website https://mackerell.umaryland.edu/charmm_ff.shtml 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:

5vbm-sol.inp
 1charmm
 2    # Here:
 3    # toppar_water_ions.str: the parameter file for water and ions.
 4    # par_all36m_prot.prm: the parameter file for proteins.
 5    # Both are obtained from website: https://mackerell.umaryland.edu/charmm_ff.shtml.
 6    parameters toppar_water_ions.str par_all36m_prot.prm
 7    topology   5vbm-sol.psf
 8    scaling14  1.0
 9    rcutoff    12.0
10    rswitch    10.0
11    use_PBC           # Add this line to turn on PBC.
12    Lbox 70 70 70     # Define the box size.
13    electrostatic pme # Use PME for electrostatics.
14    PMEk   72 72 72
15end
16
17mol
18    5vbm-sol.pdb
19end
20
21task
22    energy charmm
23end

Now run the calculation:

$ 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:

5vbm-sol.out
 1CHARMM energies
 2===============
 3Bond energy:                 7661.69410280 kcal/mol
 4Angle energy:                5772.78601511 kcal/mol
 5Dihedral energy:             1557.76413669 kcal/mol
 6Improper energy:               30.51152972 kcal/mol
 7Coulombic energy:           -2120.48409716 kcal/mol
 8Lennard-Jones energy: 11507712434764.78906250 kcal/mol
 9---------------------------------------------------
10CHARMM 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 Geometry Optimization.

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 Density Functional Theory Calculations, 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: http://zhjun-sci.com/pdbtop/doc/tutorial-2.html. 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:

92v.inp
 1charmm
 2    parameters 92v-1-naive.prm
 3    topology   92v-2.psf
 4    scaling14  1.0
 5    rcutoff    12.0
 6    rswitch    100.0 # rswitch>rcutoff means no cutoff is used.
 7end
 8
 9mol
10    92v.pdb
11end
12
13task
14    energy charmm
15end

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:

92v.out
 1CHARMM energies
 2===============
 3Bond energy:                   79.40377426 kcal/mol
 4Angle energy:                  12.30743545 kcal/mol
 5Dihedral energy:               93.12335114 kcal/mol
 6Improper energy:                0.00000000 kcal/mol
 7Coulombic energy:            -116.19015178 kcal/mol
 8Lennard-Jones energy:          31.73702659 kcal/mol
 9---------------------------------------------------
10CHARMM energy:                100.38143566 kcal/mol
11
12Final total energy:         100.38143566 kcal/mol