DivCon Command Line Tutorials

The DivCon Command Line arguments cover various standard functions and features available within the tool box. 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.

(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) platform 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 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:

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.

 

Additional Examples



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:

$ qmechanic 4EK4 --prepare --np 8 -v 2 -p pdb
$ mv 4EK4_out.pdb 4EK4.pdb
$ qmechanic 4EK4_out.pdb --np 8 -h pm6 -v 2 -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 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 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

 

Additional Examples



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
$ qmechanic 1LRI-addH.pdb -h pm6 --np 8 -v 2 --decompose "/A/CLR/99//" --dc

Where the command line switches are defined:

Switch Defaut 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.
--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 transfered one way or the other are provided as well.

 

Additional Examples



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
$ qmechanic 1LRI-addH.pdb -h pm6 --np 8 -v 2 --opt all --symmetry -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). 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.
-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.

Additional Examples



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
$ 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 Defaut 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 swtich 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 specie(s) to center the QM region 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 Optmization 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.

Additional Examples



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 that 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. In the example below, we have switched on the use of our Beta --prepare command line function to protonate the structure using our built in protonation tool.

$ qmechanic 1LRI --prepare -h pm6 amberff14 --opt 100 0.01 --region-selection /A/CLR/99// 3.0 0.0 --np 8 -v 2 -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.
--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.
--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.

--region-selection
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 specie(s) to center the QM region 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 /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.

Additional Examples



Protonation

Protonation functionality has been added to the DivCon Discovery Suite and it is currently available as a beta feature. The ONIOM tutorial above shows how to use protonation within a calculation while the present tutorial shows how to just protonate the structure and write out a PDB file.

$ qmechanic 4EK4 --prepare --np 8 -v 2 -p pdb

Switch Defaut 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.

Additional Examples