DivCon Command Line Tutorials
Note: These tutorials have been updated for DivCon Discovery Suite DEV.1013 or later. If you do not have a current version, please email sales@quantumbioinc.com for download instructions.
The DivCon Command Line arguments cover various standard functions and features available within the toolbox. If you are looking for X-ray refinement tools and options, please see the X-ray Refinement Manual. Otherwise, the documentation below will provide you with a starting point.
- Structure Preparation
- Single Point Calculation
- Interaction Energy Decomposition
- Structure Optimization
- Active Site Structure Optimization
- ONIOM (mixed-QM/MM) Simulations
- DivCon-based Protonation
(If you would like to be notified when these other functions are documented, email support@quantumbioinc.com to register your interest).
Structure Preparation
Structure preparation for QM-based methods isn’t all that different from preparation for MM-based methods. You are free to use various 3rd party graphical user interface (GUI) platforms such as MOE, Maestro, Sybyl, and so on. If you have any trouble, Contact Us for help and recommendations. It should be noted that QM, being a physics-based, all-atom method, is very sensitive to protonation, structural hot spots, poor structure, and so on. This is an underlying benefit of the method. Generally, any sort of preparation will follow the steps below:
Protonation
When a PDB is first downloaded from PDB.org or from an internal repository, the file is usually devoid of protons and only consists of heavy atom positions. Since QM is an all-atom method, the structure must absolutely have all of the protonation states of the biomolecule prior to the simulation. Further, since QM is exquisitely sensitive to protonation states, this step in the preparation process is an opportunity to explore the chemistry of the target:ligand complex and double-check whether the states you have set are correct. If you have access to the Phenix package along with the original structure factors, you can even use the Phenix/DivCon plugin to explore various states in order to get a more definitive understanding of questionable states.
Beta! DivCon-based Protonation & Structure Optimization
With the most recent release, we have announced the availability of the beta version of a DivCon-based protonation and structure optimization tool! This is the culmination of a year of development to include various protonation algorithms, methods for perception and atom typing, H-bond network optimization, database-driven pKa prediction, and topology determination methods. As such, it should still be considered Beta and you are encouraged to report any problems to support@quantumbioinc.com. Generally, this method is much faster than other methods and will soon be used not only for protonation but also for tautomer/protomer scoring. The following two tutorials are provided which use this protonation method:
- ONIOM (mixed-QM/MM) Simulations – with optimization
- DivCon-based Protonation (beta) – with no optimization
Structure Optimization
Once the structure has been protonated, depending upon the quality of the structure, you may want to perform a structure optimization using a standard force field – AMBER, OPLS, and so on – in order to clean up any structural “hot spots” or other problems. These hot spots can lead to longer QM calculations and they will impact results. Often this step may include a tether or restraint in order to limit the movement of the heavy atoms if these atom positions are fairly well defined. MOE, Maestro, or several other 3rd party software packages can be used to perform this optimization.
Docking/Scoring
Finally, if you are working with a novel ligand, you may need to dock this ligand into the target structure. This sort of preparation is beyond the scope of this document as there is a myriad of possible docking software packages. In our experience, more often than not, it seems that the atomic positions of any docked structures should be optimized using a force field prior to further analysis. Even when the docking software is configured to “refine” or “optimize” or “minimize” the ligand pose, it is not unusual to see the ligand atoms placed too close to the heavy atoms within the target.
From there you can either run additional preparation within your platform of choice, or you can proceed to one of the next stages with DivCon. When performing any QM calculations using the DivCon command line, you will use the qmechanic software available within $QBHOME/bin directory. This software can be invoked directly or through plugins to MOE and Phenix. As will be detailed in a subsequent version, it can also be used as a WebService and be called from Pipeline Pilot, KNIME, or other software.
Single Point Calculation
The default calculation is the single-point characterization of the input structure as provided. This type of calculation can be used to generate energies, atomic charges, various analyses, and so on. The following command line can be used to invoke this functionality:
% wget http://downloads.quantumbioinc.com/media/tutorials/cli/1TOW-H.pdb
% mv 1TOW-H.pdb 1TOW.pdb
% /path/to/DivConSuite/bin/qmechanic 1TOW.pdb --np 8 -h pm6 -v 2
Where the command line switches are defined:
| Switch | Default | Description |
| Input File or PDBid | Required | The qmechanic command line can accept most standard PDB, MOL2, and CML/XML files. In this case, this is the fully prepared PDB file. |
| –np | 2 | The qmechanic tool – along with most QuantumBio technologies – is optimized for multi-threaded parallelism and can therefore be used on multi-core machines in order to maximize performance. Generally, the larger the QM region of the calculation, the better (e.g. more linear) parallelism you will experience. |
| -h | pm6 | qmechanic supports several methods or Hamiltonians for both semiempirical quantum mechanics and soon molecular mechanics simulations. The -h switch allows you to pick which method should be used within the calculation. The system will automatically perform any atom typing, formal charge determination, and so on required for a particular method. |
| -v | 1 | Verbosity can be set to provide more information on the screen in order to keep track of the simulation. If the default verbosity is chosen, then the screen output will only include a minimal amount of information and the log file will need to be used to measure progress. A setting of 2 will provide information corresponding to each cycle of the simulation. |
| -p | NONE | If PDB is chosen as the option, upon completion of the calculation, output the pdb formatted version of the final coordinates. |
Unless a verbosity level of 0 is chosen, qmechanic will provide information on any errors or warnings. If warnings such as these are seen, this indicates a possible problem in protonation. The –fix-open-shell switch will automatically “fix” these problems by adding protons to fill any open valences, but depending upon your goals, you may not want these valences addressed.
Perception Warning: 1 Hydrogen atom(s) are missing at atom CB in SER 51 A Perception Warning: 3 Hydrogen atom(s) are missing at atom CB in VAL 79 A
Other output will include the following general system information. These characteristics are provided in order to aide you in debugging any problems.
System Information Atoms: 8893 Total Charge: 5 Electrons: 24526 Molecular Orbitals: 22225 Requested Memory: 106071 mb Total Memory Requirements: 1898 mb
Interaction Energy Decomposition
As discussed in J. Phys. Chem. A 1999, 103, 3321-3329, the total interaction energy between a target and its ligand can be decomposed into polarization, electrostatics, and charge transfer using a multi-step quantum mechanics workflow. This workflow has been simplified and captured within the qmechanic application allowing you to calculate the total interaction energy along with its decomposition using a single set of command line arguments.
The 1LRI-addH.pdb used in the following example is available for download here.
% wget http://downloads.quantumbioinc.com/media/tutorials/cli/1LRI-addH.pdb
% /path/to/DivConSuite/bin/qmechanic 1LRI-addH.pdb -h pm6 --np 8 -v 2 --decompose "/A/CLR/99//" --dc
Where the command line switches are defined:
| Switch | Default | Description |
| Input File | Required | The qmechanic command line can accept most standard PDB, MOL2, and CML/XML files. In this case, this is the fully prepared PDB file (the –input-file switch itself can be implied with the inclusion of a recognized input file). |
| –np | 2 | The qmechanic tool – along with most QuantumBio technologies – is optimized for multi-threaded parallelism and can therefore be used on multi-core machines in order to maximize performance. Generally, the larger the QM region of the calculation, the better (e.g. more linear) parallelism you will experience. |
| -h | pm6 | qmechanic supports several methods or Hamiltonians for both semiempirical quantum mechanics and soon molecular mechanics simulations. The -h switch allows you to pick which method should be used within the calculation. The system will automatically perform any atom typing, formal charge determination, and so on required for a particular method. |
| -v | 1 | Verbosity can be set to provide more information on the screen in order to keep track of the simulation. If the default verbosity is chosen, then the screen output will only include a minimal amount of information and the log file will need to be used to measure progress. A setting of 2 will provide information corresponding to each cycle of the simulation. |
| –decompose | Required | The decomposition workflow is invoked using this command line switch where the (required) selection designates the ligand or analogous species to which the interaction energy should be calculated. In this case “/L/CLR/99//” is defined where the “L” denotes the chain, the “CLR” references the residue name, and the “99” corresponds to the residue ID. A selection is required otherwise the calculation will throw an error. |
The output from the simulation will follow very closely with that noted in the Single Point Calculated listed above where in this case four individual calculations will be performed instead of one. Whenever possible, the density matrix from a previous step is used in subsequent steps in order to speed up the calculation.
Processors: 8
Warning Single alternate A A LYS 39 CG
Warning Single alternate A A LYS 39 CD
Warning Single alternate A A LYS 39 CE
Warning Single alternate A A LYS 39 NZ
Warning Single alternate A A THR 65 CB
Warning Single alternate A A THR 65 OG1
Warning Single alternate A A THR 65 CG2
Warning Single alternate A A SER 92 C
Warning Single alternate A A SER 92 O
Warning Single alternate A A SER 92 CB
Warning Single alternate A A SER 92 OG
Warning Single alternate A A ASN 93 CB
Warning Single alternate A A ASN 93 CG
Warning Single alternate A A ASN 93 ND2
Warning Single alternate A A ASN 93 OD1
System Information
Atoms: 1804
Total Charge: 2
Electrons: 4964
Molecular Orbitals: 4386
Running Selection
Requested Memory: 103246 mb
Total Memory Requirements: 2 mb
--------------------------------------------------------------------------- -----
Single Point
Cycle Energy (eV) Delta E D RMSD Max Delta D Damp Time
--------------------------------------------------------------------------- -----
1 -4.28428072e+04 -1.0492e+01 7.4846e-03 1.2211e-01 1.00 0.01
2 -4.28437318e+04 -9.2463e-01 2.2559e-03 3.5269e-02 1.00 0.01
3 -4.28438995e+04 -1.6766e-01 9.7325e-04 1.7827e-02 1.00 0.01
...
Total SCF Time (Seconds): 0.33699
Electronic Energy (eV): -42843.96982
Total Energy (eV): -4225.95743
Heat of Formation (kcal/mol): -102.18634
Ionization Potential (eV): 9.34891
Electron Affinity (eV): -1.29458
Running Compliment
Requested Memory: 103246 mb
Total Memory Requirements: 475 mb
--------------------------------------------------------------------------- -----
Single Point
Cycle Energy (eV) Delta E D RMSD Max Delta D Damp Time
--------------------------------------------------------------------------- -----
1 -1.20302473e+07 -5.0345e+02 9.7810e-03 1.9166e+00 1.00 1.74
2 -1.20302573e+07 -1.0002e+01 9.5584e-03 1.9639e+00 0.06 1.73
3 -1.20302871e+07 -2.9795e+01 5.8163e-03 1.8498e+00 0.27 1.92
....
Total SCF Time (Seconds): 51.10683
Electronic Energy (eV): -12030429.13432
Total Energy (eV): -161196.45642
Heat of Formation (kcal/mol): -11806.30630
Ionization Potential (eV): 0.00000
Electron Affinity (eV): 0.00000
Running Complex with No Charge Transfer
--------------------------------------------------------------------------- -----
Single Point
Cycle Energy (eV) Delta E D RMSD Max Delta D Damp Time
--------------------------------------------------------------------------- -----
1 -1.29683270e+07 -6.3867e-02 9.9470e-05 3.2056e-02 1.00 1.91
2 -1.29683270e+07 -9.4217e-03 3.6061e-05 9.9193e-03 1.00 1.77
3 -1.29683270e+07 -2.1501e-03 1.7306e-05 4.0261e-03 1.00 1.79
....
Total SCF Time (Seconds): 20.06500
Electronic Energy (eV): -12968327.01374
Total Energy (eV): -165421.55448
Heat of Formation (kcal/mol): -11888.67504
Ionization Potential (eV): 0.00000
Electron Affinity (eV): 0.00000
Running Complex
Requested Memory: 103246 mb
Total Memory Requirements: 712 mb
--------------------------------------------------------------------------- -----
Single Point
Cycle Energy (eV) Delta E D RMSD Max Delta D Damp Time
--------------------------------------------------------------------------- -----
1 -1.29683286e+07 -1.6010e+00 3.0341e-04 1.0391e-01 1.00 2.48
2 -1.29683288e+07 -1.3640e-01 9.3024e-05 3.6148e-02 1.00 2.36
3 -1.29683288e+07 -1.5096e-02 3.2478e-05 1.3109e-02 1.00 2.37
....
Total SCF Time (Seconds): 23.63456
Electronic Energy (eV): -12968328.76866
Total Energy (eV): -165423.30940
Heat of Formation (kcal/mol): -11929.14442
Ionization Potential (eV): 0.00000
Electron Affinity (eV): 0.00000
Total Interaction Energy (kcal/mol): -20.65
Electrostatic Energy (kcal/mol): 21.58 (33.82%)
Polarization Energy (kcal/mol): -1.76 ( 2.76%)
Charge Transfer Formation (kcal/mol): -40.47 (63.42%)
e- (Compliment --> Selection): -0.02
For this example simulation, the alternate atom warning is noted and corresponds to the fact that there are alternate atom locations within the PDB file. Since qmechanic supports X-ray refinement methods through its integration with Phenix, it will process alternate atom conformations (and use the “A” conformation for all simulations). This warning gives you the opportunity to double check to make sure that you have provided qmechanic with the correct conformation.
Once the simulation has completed, the final interaction energies are presented including the Total, Electrostatic, Polarization, and Charge Transfer. The number of electrons transferred one way or the other is provided as well.
Structure Optimization
There are times when QM-based structure optimization would be preferred prior to interaction energy calculation or other simulations. Just as with the other applications, it is strongly suggested that you still perform a full preparation prior to performing the optimization. Once the structure is fully prepared, the optimization is invoked using the command line options detailed below.
The 1LRI-addH.pdb used in the following example is available for download here.
% wget http://downloads.quantumbioinc.com/media/tutorials/cli/1LRI-addH.pdb
% /path/to/DivConSuite/bin/qmechanic 1LRI-addH.pdb -h pm6 --np 8 -v 2 --opt all 3 0.01 --symmetry -p pdb
Where the command line switches are defined:
| Switch | Default | Description |
| Input File or PDBid | Required | The qmechanic command line can accept most standard PDB, MOL2, and CML/XML files. In this case, this is the fully prepared PDB file. |
| –prepare | This command line (coupled with a PDBid) will download the PDB file and automatically protonate it. | |
| –np | 1 | The qmechanic tool – along with most QuantumBio technologies – is optimized for multi-threaded parallelism and can therefore be used on multi-core machines in order to maximize performance. Generally, the larger the QM region of the calculation, the better (e.g. more linear) parallelism you will experience. |
| -h | pm6 | qmechanic supports several methods or Hamiltonians for both semiempirical quantum mechanics and soon molecular mechanics simulations. The -h switch allows you to pick which method should be used within the calculation. The system will automatically perform any atom typing, formal charge determination, and so on required for a particular method. |
| -v | 1 | Verbosity can be set to provide more information on the screen in order to keep track of the simulation. If the default verbosity is chosen, then the screen output will only include a minimal amount of information and the log file will need to be used to measure progress. A setting of 2 will provide information corresponding to each cycle of the simulation. |
| -p | NONE | If PDB is chosen as the option, upon completion of the calculation, output the pdb formatted version of the final coordinates. |
The output from the simulation will follow very closely with that noted in the Single Point Calculated listed above:
Processors: 8
Warning Single alternate A A LYS 39 CG
Warning Single alternate A A LYS 39 CD
Warning Single alternate A A LYS 39 CE
Warning Single alternate A A LYS 39 NZ
Warning Single alternate A A THR 65 CB
Warning Single alternate A A THR 65 OG1
Warning Single alternate A A THR 65 CG2
Warning Single alternate A A SER 92 C
Warning Single alternate A A SER 92 O
Warning Single alternate A A SER 92 CB
Warning Single alternate A A SER 92 OG
Warning Single alternate A A ASN 93 CB
Warning Single alternate A A ASN 93 CG
Warning Single alternate A A ASN 93 ND2
Warning Single alternate A A ASN 93 OD1
System Information
Atoms: 1804
Total Charge: 2
Electrons: 4964
Molecular Orbitals: 4386
Requested Memory: 103971 mb
Total Memory Requirements: 882 mb
--------------------------------------------------------------------------- -----
Single Point
Cycle Energy (eV) Delta E D RMSD Max Delta D Damp Time
--------------------------------------------------------------------------- -----
1 -1.29681355e+07 -4.7933e+02 7.5518e-03 1.9236e+00 1.00 3.16
2 -1.29681411e+07 -5.5259e+00 1.1341e-02 1.9998e+00 0.02 2.92
3 -1.29681448e+07 -3.7057e+00 8.4044e-03 1.9684e+00 0.02 3.08
4 -1.29681577e+07 -1.2962e+01 6.2165e-03 1.9171e+00 0.09 3.13
....
Total SCF Time (Seconds): 105.28504
Electronic Energy (eV): -12968356.73306
Total Energy (eV): -165451.27381
Heat of Formation (kcal/mol): -12574.01925
Ionization Potential (eV): 0.00000
Electron Affinity (eV): 0.00000
Total Population Analysis Time (Seconds): 0.67444
--------------------------------------------------------------------------- -----
Optimization
Cycle Total Energy (eV) Delta E Grad. Norm Grad. Max Coord. Norm
--------------------------------------------------------------------------- -----
0 -1.65451274e+05 -1.6545e+05 5.6550e+01 1.2487e+01 0.0000e+00
1 -1.65464363e+05 -1.3089e+01 6.0254e+01 1.5780e+01 0.0000e+00
2 -1.65476360e+05 -1.1998e+01 2.5711e+01 5.7909e+00 9.8262e-01
3 -1.65481805e+05 -5.4445e+00 1.7337e+01 3.0524e+00 3.2795e-01
4 -1.65488170e+05 -6.3653e+00 1.6233e+01 5.3432e+00 2.6463e-01
5 -1.65492458e+05 -4.2879e+00 1.8380e+01 4.0560e+00 5.3172e-01
6 -1.65495287e+05 -2.8290e+00 1.0116e+01 1.0943e+00 6.0731e-01
Key = Chain Name/Residue Name/Residue Seq./Atom Index/Atom Type/Transform Index
Collisions between symmetry related residues:
--------------------------------------------------------------------------- -----
A/CL/100/1432/CL/0 A/CL/100/1432/CL/12288 0.162261
A/CL/100/1432/CL/12288 A/CL/100/1432/CL/0 0.162261
A/CL/100/1432/CL/0 A/CL/100/1432/CL/12288 0.162261
A/CL/100/1432/CL/12288 A/CL/100/1432/CL/0 0.162261
A/CL/100/1432/CL/0 A/CL/100/1432/CL/12288 0.162261
A/CL/100/1432/CL/12288 A/CL/100/1432/CL/0 0.162261
A/CL/100/1432/CL/0 A/CL/100/1432/CL/12288 0.162261
A/CL/100/1432/CL/12288 A/CL/100/1432/CL/0 0.162261
A/CL/100/1432/CL/0 A/CL/100/1432/CL/12288 0.162261
A/CL/100/1432/CL/12288 A/CL/100/1432/CL/0 0.162261
A/CL/100/1432/CL/0 A/CL/100/1432/CL/12288 0.162261
A/CL/100/1432/CL/12288 A/CL/100/1432/CL/0 0.162261
A/CL/100/1432/CL/0 A/CL/100/1432/CL/12288 0.162261
A/CL/100/1432/CL/12288 A/CL/100/1432/CL/0 0.162261
A/CL/100/1432/CL/0 A/CL/100/1432/CL/12288 0.162261
A/CL/100/1432/CL/12288 A/CL/100/1432/CL/0 0.162261
A/CL/100/1432/CL/0 A/CL/100/1432/CL/12288 0.162261
A/CL/100/1432/CL/12288 A/CL/100/1432/CL/0 0.162261
....
13 -1.65719453e+05 9.3967e+01 5.7710e+01 1.9696e+01 5.8442e-02
Total Optimization Time (Seconds): 4830.65934
Electronic Energy (eV): -12958414.29832
Total Energy (eV): -165719.45305
Heat of Formation (kcal/mol): -18758.38257
Ionization Potential (eV): 0.00000
Electron Affinity (eV): 0.00000
Total Population Analysis Time (Seconds): 0.37609
As in other cases, the software will provide any structural warnings, information about the input, and so on. Then it will perform the first self consistent field (SCF) calculation and run until SCF convergence. This will provide starting energy data on the initial structure of the optimization and create an initial density matrix that will be used in the first optimization cycle.
For each optimization cycle, another SCF convergence step (not shown) is performed, and the atomic gradients for the structure are generated and provided to the L-BFGS optimization algorithm. Each cycle repeats the same set of calculations until the optimization converges. In this case, due to changes in coordinates of the special position atoms within the structure, a number of collisions were detected during the simulation. These are addressed and the optimization will continue. Finally, upon completion of the optimization, the final energy values are provided. In this case, there is a drop of over 6000 kcal/mol in heat of formation.
The qmechanic tool uses the non-proprietary, HDF5 binary file format in order to store the data from the simulation. The qbreporter tool is used to process this h5 file and output the PDB files for the initial and the final calculations. If you wish to view the contents of the file for other pieces of data, use the hdfview tool distributed within the package. You are also free to use the various tools available on the HDF5 website to process the file in Python or other languages.
Active Site Structure Optimization
There are times when QM-based structure optimization would be preferred prior to interaction energy calculation or other simulations. Just as with the other applications, it is strongly suggested that you still perform a full preparation prior to performing the optimization. Once the structure is fully prepared, the optimization is invoked using the command line options detailed below.
The 1P2Y-addH.pdb used in the following example is available for download here.
% wget http://downloads.quantumbioinc.com/media/tutorials/cli/1P2Y-addH.pdb
% /path/to/DivConSuite/bin/qmechanic 1P2Y-addH.pdb --opt /*/HEM/*//,/*/NCT/*// 5.0 0.0 --np 8 -h pm6 -v 2 -p pdb
Where the command line switches are defined:
| Switch | Default | Description |
| Input file | Required | The qmechanic command line can accept most standard PDB, MOL2, and CML/XML files. In this case, this is the fully prepared PDB file (the –input-file switch itself can be implied with the inclusion of a recognized input file). |
| –np | 1 | The qmechanic tool – along with most QuantumBio technologies – is optimized for multi-threaded parallelism and can therefore be used on multi-core machines in order to maximize performance. Generally, the larger the QM region of the calculation, the better (e.g. more linear) parallelism you will experience. |
| -h | pm6 | qmechanic supports several methods or Hamiltonians for both semiempirical quantum mechanics and soon molecular mechanics simulations. The -h switch allows you to pick which method should be used within the calculation. The system will automatically perform any atom typing, formal charge determination, and so on required for a particular method. |
| -v | 1 | Verbosity can be set to provide more information on the screen in order to keep track of the simulation. If the default verbosity is chosen, then the screen output will only include a minimal amount of information and the log file will need to be use to measure progress. A setting of 2 will provide information corresponding to each cycle of the simulation. |
| –opt | all |
This switch invokes the optimization routines within qmechanic in order to perform an energy minimization on the structure. Currently, you can choose to optimize all of the atoms in the structure or just the hydrogens (by using the “hydrogen” argument). In this example, we use the same selection syntax as with other switches. Using the qmechanic-standard format for selections, this switch can be used (optionally) to choose the ligand, a particular amino acid residue, or some other species to center the QM region(s) of the QM/MM calculation. In this case, all atoms will more as per the ONIOM Hamiltonian, but only these atoms will be treated quantum mechanically. The format for selections is: “/Chain/Residue/UID//” (in this case /*/HEM/*//,/*/NCT/*//). The second and third values for this command line switch are the selection radius and buffer respectively. The buffer is not used for the current context and should be set to 0.0. The selection radius will expand the selection to be treated quantum mechanically within the ONIOM calculation to include all residues within that distance (in this case 3.0 Å). As always, selections are built in a residue-extended manner where if a single atom of any residue is within 3.0 Å, then the entire active site residue is included in the calculation. |
| -p | NONE | If PDB is chosen as the option, upon completion of the calculation, output the pdb formatted version of the final coordinates. Mol2 may also be chosen. |
The output from the simulation will follow very closely with that noted in the Structure Optimization listed above:
Processors: 15
System Information
Atoms: 6876
Total Charge: -10
Electrons: 19158
Molecular Orbitals: 17183
Requested Memory: 114316 mb
Total Memory Requirements: 2717 mb
--------------------------------------------------------------------------------
Single Point
Cycle Energy (eV) Delta E D RMSD Max Delta D Damp Time
--------------------------------------------------------------------------------
1 -1.19916772e+08 -2.4574e+01 1.4276e-02 1.9988e+00 0.01 9.80
2 -1.19916814e+08 -4.2653e+01 1.0323e-02 1.9744e+00 0.03 9.86
3 -1.19917377e+08 -5.6298e+02 7.4994e-03 1.9055e+00 0.56 9.59
4 -1.19917402e+08 -2.4549e+01 1.3371e-02 1.9212e+00 0.01 9.71
5 -1.19917450e+08 -4.8184e+01 8.4422e-03 1.6556e+00 0.03 9.54
6 -1.19917564e+08 -1.1453e+02 5.6907e-03 1.5794e+00 0.13 9.63
7 -1.19917955e+08 -3.9032e+02 4.6832e-03 1.3368e+00 0.66 9.66
8 -1.19917964e+08 -8.9591e+00 1.1833e-02 1.9678e+00 0.01 9.68
9 -1.19917976e+08 -1.2615e+01 8.3181e-03 1.7839e+00 0.02 9.50
10 -1.19918007e+08 -3.0408e+01 5.3553e-03 1.4487e+00 0.06 9.70
11 -1.19918146e+08 -1.3936e+02 3.4657e-03 1.3648e+00 0.41 9.66
12 -1.19918215e+08 -6.8978e+01 2.7581e-03 9.9511e-01 0.29 9.58
....
70 -1.19919300e+08 0.0000e+00 1.4015e-04 8.6507e-02 -0.00 13.89
71 -1.19919300e+08 0.0000e+00 1.4022e-04 8.6549e-02 -0.00 14.82
72 -1.19919300e+08 0.0000e+00 1.4025e-04 8.6569e-02 -0.00 13.43
73 -1.19919300e+08 0.0000e+00 1.4026e-04 8.6579e-02 -0.00 13.12
Total SCF Time (Seconds): 838.60212
Electronic Energy (eV): -119919299.52772
Total Energy (eV): -620647.21206
Heat of Formation (kcal/mol): -32682.16011
Ionization Potential (eV): 0.00000
Electron Affinity (eV): 0.00000
Total Population Analysis Time (Seconds): 1.55600
--------------------------------------------------------------------------------
Optimization
Cycle Total Energy (eV) Delta E Grad. Norm Grad. Max Coord. Norm
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
Selected Optimization
Cycle Total Energy (eV) Delta E Grad. Norm Grad. Max Coord. Norm
--------------------------------------------------------------------------------
0 -6.20647212e+05 -1.4312e+07 1.1926e+02 1.1592e+01 0.0000e+00
1 -6.20652037e+05 -1.1127e+02 1.2237e+02 1.1582e+01 0.0000e+00
2 -6.20662616e+05 -2.4395e+02 1.1298e+02 1.1586e+01 9.7508e-01
3 -6.20664537e+05 -4.4301e+01 1.1206e+02 1.1591e+01 3.4767e-01
....
28 -6.20677830e+05 -2.0422e+00 1.1169e+02 1.1579e+01 3.0179e-01
29 -6.20677946e+05 -2.6686e+00 1.1163e+02 1.1579e+01 1.5510e-01
30 -6.20677849e+05 2.2456e+00 1.1164e+02 1.1579e+01 1.4718e-01
Total Optimization Time (Seconds): 25188.92760
Electronic Energy (eV): -119923161.39025
Total Energy (eV): -620677.84870
Heat of Formation (kcal/mol): -33388.65811
Ionization Potential (eV): 0.00000
Electron Affinity (eV): 0.00000
Total Population Analysis Time (Seconds): 1.97603
Job Complete Total Computation Time (Seconds): 26041.11343
As in other cases, the software will provide any structural warnings, information about the input, and so on. Then it will perform the first self-consistent field (SCF) calculation and run until SCF convergence. This will provide starting energy data on the initial structure of the optimization and create an initial density matrix that will be used in the first optimization cycle.
For each optimization cycle, another SCF convergence step (not shown) is performed, and the atomic gradients for the structure are generated and provided to the L-BFGS optimization algorithm. Each cycle repeats the same set of calculations until the optimization converges.
The qmechanic tool uses the non-proprietary, HDF5 binary file format in order to store the data from the simulation. The qbreporter tool is used to process this h5 file and output the PDB files for the initial and the final calculations. If you wish to view the contents of the file for other pieces of data, use the hdfview tool distributed within the package. You are also free to use the various tools available on the HDF5 website to process the file in Python or other languages
ONIOM (mixed-QM/MM) Simulations
The mixed-QM/MM (ONIOM) method as implemented in the DivCon Discovery Suite can be used with most energy-based applications (including single point, optimization, and X-ray refinement). As implied by the description, the method will apply quantum mechanics to the selections (in this case the ligand and the surrounding active site) and then use a molecular mechanics force field for the remaining structure. Often, this is a faster calculation than pure QM calculations, and it is going through continued development to make it even faster over time. Unlike other implementations, all-atom typing, perception, topology determination, selection expansion, and capping are done completely automatically. Just as with the other applications, it is strongly suggested that you still perform a full preparation prior to performing the optimization.
% wget http://downloads.quantumbioinc.com/media/tutorials/cli/1LRI-addH.pdb
% /path/to/DivConSuite/bin/qmechanic 1LRI-addH.pdb -h pm6 amberff14sb --opt 100 0.01 --qm-region /A/CLR/99// 3.0 0.0 --np 8 -v 1 -p pdb
Where the command line switches are defined:
| Switch | Default | Description |
| Input file or PDBid | Required | In this case, the 4-character PDBid is used for the input. qmechanic will automatically download the data from the RCSB.org website and begin the calculation. You may download the PDB file directly from the RCSB and place it on the command line at this same location. |
| –np | 1 | The qmechanic tool – along with most QuantumBio technologies – is optimized for multi-threaded parallelism and can therefore be used on multi-core machines in order to maximize performance. Generally, the larger the QM region of the calculation, the better (e.g. more linear) parallelism you will experience. |
| -h | Required | For ONIOM (QM/MM), this option is required. By default the value for the -h command line option is PM6; however, to run QM/MM two indicators must be supplied: the QM Hamiltonian (PM6, PM3, or AM1) followed by the designation for the AMBER force field (available options include amberff14, amberff12, amberff10, amberff99, and amberff94). |
| -v | 2 | Verbosity can be set to provide more information on the screen in order to keep track of the simulation. If the default verbosity is chosen, then the screen output will only include a minimal amount of information and the log file will need to be use to measure progress. A setting of 2 will provide information corresponding to each cycle of the simulation. |
| –opt | all |
This switch invokes the optimization routines within qmechanic in order to perform an energy minimization on the structure. Currently, you can choose to optimize all of the atoms in the structure or just the hydrogens (by using the “hydrogen” argument). The second two options are convergence limits in the number of cycles and step-to-step change in energy (in kcal/mol). Note: if your structure is a crystal structure and if you have access to the Phenix crystallographic platform, it is recommended that you perform an X-ray refinement using Phenix/DivCon instead of performing a “free form” QM optimization. |
| –qm-region | NONE |
Using the qmechanic-standard format for selections, this switch can be used (optionally) to choose the ligand, a particular amino acid residue, or some other species to center the QM region(s) of the QM/MM calculation. In this case, all atoms will move as per the ONIOM Hamiltonian, but only these atoms will be treated quantum mechanically. The format for selections is: “/Chain/Residue/UID//” (in this case /A/CLR/99//). The second and third values for this command line switch are the selection radius and buffer respectively. The buffer is not used for the current context and should be set to 0.0. The selection radius will expand the selection to be treated quantum mechanically within the ONIOM calculation to include all residues within that distance (in this case 3.0 Å). As always, selections are built in a residue-extended manner where if a single atom of any residue is within 3.0 Å, then the entire active site residue is included in the calculation. |
| -p | NONE | If PDB is chosen as the option, upon completion of the calculation, output the pdb formatted version of the final coordinates. |
The output from the simulation will follow very closely with that noted in the Structure Optimization listed above:
System Information
Atoms: 541
Total Charge: 0
Electrons: 1340
Molecular Orbitals: 1261
Nonstandard Residues: A/99/CLR/
1 -4.28361063e+04 -1. 0949e+01 7.7903e-03 1.5400e-01 1 .00 0.55
2 -4.28371071e+04 -1. 0008e+00 2.3832e-03 6.5336e-02 1 .00 0.56
3 -4.28373108e+04 -2. 0369e-01 1.0931e-03 3.7302e-02 1 .00 0.56
....
----------------------------------------------------------------------------- ---
Single Point
Cycle Energy (eV) & nbsp; Delta E D RMSD Max Delta D Damp ; Time
----------------------------------------------------------------------------- ---
1 -1.37559539e+06 -1. 3712e+02 3.6266e-03 3.3384e-01 1 .00 0.55
2 -1.37560524e+06 -9. 8545e+00 9.2582e-04 6.8333e-02 1 .00 0.55
....
; Total SCF Time (Seconds): ; & nbsp; 12.98250
; Electronic Energy (eV): -1375607.71945
; & nbsp;Total Energy (eV): & nbsp; -39871.67084
; Heat of Formation (kcal/mol): & nbsp; -1106.71048
; Ionization Potential (eV): & nbsp; ; 8.09942
; Electron Affinity (eV): ; 0.77856
qm total energy (selection eV): -39871.67084
qm total energy (selection kcal/mol) & nbsp; : -919463.01362
mm total energy (entire system kcal/mol) : -317.75037
mm total energy (selection system kcal/mol ):& nbsp;-8.43205
oniom energy kcal/mol = -919772.33193
----------------------------------------------------------------------------- ---
Optimization
Cycle Total Energy (eV) & nbsp; Delta E Grad. Norm & nbsp; Grad. Max Coord. Norm
----------------------------------------------------------------------------- ---
0 -9.19772332e+05 & nbsp;-9.1977e+05 4.1462e+01 6. 5795e+00 0.0000e+00
1 -9.19806611e+05 & nbsp;-3.4279e+01 6.4851e+01 1. 0231e+01 0.0000e+00
2 -9.20100949e+05 & nbsp;-2.9434e+02 2.0436e+01 3. 1061e+00 9.7645e-01
.....
100 -9.20872961e+05 -1. 5095e+00 1.8788e+00 3.6345e-01& nbsp; 1.3096e-01
As in other cases, the software will provide any structural warnings, information about the input, and so on. Then it will perform the first self consistent field (SCF) calculation and run until SCF convergence. This will provide starting energy data on the initial structure of the optimization and create an initial density matrix that will be used in the first optimization cycle.
Since the calculation also includes an “unknown” residue (CLR in this case) it will be treated quantum mechanically in order to determine charges and the like.
For each optimization cycle, another SCF convergence step (not shown) is performed, and the atomic gradients for the structure are generated and provided to the L-BFGS optimization algorithm. Each cycle repeats the same set of calculations until the optimization converges or until the criteria are reached.
Upon completion, the new, optimized and prepared structure is written out to a PDB file with the name, in this case of 1LRI_out.pdb (download this file here). The image blow depicts the downward trend of the optimization (as measured by total energy) for the 100 cycle steps chosen.
Protonation
Protonation functionality has been added to the DivCon Discovery Suite and it is currently available as a beta feature. With current packages, the method is not available, but it will be available in the coming weeks.
N/A - check back soon!
| Switch | Default | Description |
| Input file or PDBid | Required | In this case, the 4-character PDBid is used for the input. qmechanic will automatically download the data from the RCSB.org website and begin the calculation. You may download the PDB file directly from the RCSB and place it on the command line at this same location. |
| –prepare | NONE | Then the –prepare switch is provided, protonation and H-bond network optimization will be performed on the entire structure using a combination of dead end elimination and graph theory along with crystallographic symmetry and density – as available. Note: by default, only preparation will be performed and no calculations (QM, MM, or QM/MM) will execute. If you wish to optimize the added protons using amberff14, add the following command line arguments: -h amberff14 –opt hydrogen |
| –np | 1 | The qmechanic tool – along with most QuantumBio technologies – is optimized for multi-threaded parallelism and can therefore be used on multi-core machines in order to maximize performance. |
| -v | 2 | Verbosity can be set to provide more information on the screen in order to keep track of the simulation. If the default verbosity is chosen, then the screen output will only include a minimal amount of information and the log file will need to be use to measure progress. A setting of 2 will provide information corresponding to each cycle of the simulation. |
| -p | NONE | If PDB is chosen as the option, upon completion of the protonation, a PDB-formatted version of the file will be written. If MOL2 is entered as the command line option, the mol2 formatted version will be written. |