Saturday, August 29, 2015

AFNMR: automated fragmentation quantum mechanical calculation of NMR chemical shifts for biomolecules


Contributed by Jan Jensen

Xiao He and co-workers have in the last few years been working on a fragmentation methodology for computing NMR chemical shifts and this methodology has now been automated an implemented in the David Case's SHIFTS program.

In this method a H-capped model system is created for each residue, while the remaining protein is treated as point charges.  Bulk solvent effects are represented by surface charges computed using a Poisson-Boltzmann calculation based on a point charge representation of the protein while explicit water molecules also can be added using the PLACEVENT program.  Much of this has been automated in the new AFNMR module in SHIFTS program:
We have attempted to automate the process as much as possible, so that default calculations require only a PDB file as input. The preliminary processing creates fragment input files for the Gaussian, ORCA, Q-Chem or deMon3k programs; analysis programs parse the quantum chemistry output files to create tables of computed shifts and to make comparisons with experimental data if it is available. Optional parameters control the level of calculation and basis set, and the type of explicit or implicit solvent model that is used.
Furthermore, Xiao He tells me that if both SHIFTS and MEAD are installed then AFNMR will generate the surface charges automatically (alternatively they can be generated by Delphi or DIVCON in a non-automated fashion).  While the PLACEVENT program is available in AMBER it has not been fully integrated with AFNMR yet.

The authors report that it takes 1-3 hrs on a single 8-core node per residue so calculations on a 100-residue are feasible with relatively modest computational resources and with access to supercomputers one could imagine computing chemical shifts for comparatively large proteins relatively routinely.

The RMSD relatively to conventional full-DFT calculations on very small proteins are 0.2, 0.8, and 1.2 ppm for H, C, and N atoms respectively.  However, the agreement with experimental values is considerably worse and empirical predictors such as SHIFTX2 give rise to lower RMSD values.  This is a well known "problem" and may, in part, reflect small errors in the protein structures and the neglect of conformational averaging that has been "parameterized away" in the empirical methods. In fact, SHIFTX2-results in particular is rather insensitive to the quality of the structure (see Figure 4 in this paper), which is a plus if the prime objective is to get chemical shifts that agree well with experiment, but may be a minus if the information is to be used for protein structure validation and refinement.  This obviously needs further study and the AFNMR method, which promises to make QM-prediction of chemical shifts routine for proteins, is an important step in that direction.

I thank Kresten Lindorff-Larsen for first alerting me to this paper


This work is licensed under a Creative Commons Attribution 4.0

Friday, August 28, 2015

Four Decades of the Chemistry of Planar Hypercoordinate Compounds

Yang, L.-M.; Ganz, E.; Chen, Z.; Wang, Z.-X.; Schleyer, P. v. R. Angew. Chem. Int. Ed. 2015, 54, 9468-9501
Contributed by Steven Bachrach.
Reposted from Computational Organic Chemistry with permission

Yang, Ganz, Chen, Wang, and Schleyer have published a very interesting and comprehensive review of planar hypercoordinate compounds, with a particular emphasis on planar tetrahedral carbon compounds.1 A good deal of this review covers computational results.

There are two major motifs for constructing planar tetrahedral carbon compounds. The first involves some structural constraints that hold (or force) the carbon into planarity. A fascinating example is 1 computed by Rasmussen and Radom in 1999.2 This molecule taxed their computational resources, and as was probably quite typical for that time, there is no supplementary materials. But since this compound has high symmetry (D2h) I reoptimized its structure at ω-B97X-D/6-311+G(d) and computed its frequencies in just a few hours. This structure is shown in Figure 1. However, it should be noted that at this computational level, 1 possesses a single imaginary frequency corresponding to breaking the planarity of the central carbon atom. Rasmussen and Radom computed the structure of 1 at MP2/6-31G(d) with numerical frequencies all being positive. They also note that the B3LYP/6-311+G(3df,2p) structure also has a single imaginary frequency.

A second approach toward planar tetrahedral carbon compounds is electronic: having π-acceptor ligands to stabilize the p-lone pair on carbon and σ-donating ligands to help supply sufficient electrons to cover the four bonds. Perhaps the premier simple example of this is the dication 2¸ whose ω-B97X-D/6-311+G(d,p) structure is also shown in Figure 1.

The review covers heteroatom planar hypercoordinate species as well. It also provides brief coverage of some synthesized examples.

1

2
Figure 1. Optimized structures of 1 and 2.


References

(1) Yang, L.-M.; Ganz, E.; Chen, Z.; Wang, Z.-X.; Schleyer, P. v. R. "Four Decades of the Chemistry of Planar Hypercoordinate Compounds," Angew. Chem. Int. Ed. 201554, 9468-9501, DOI:10.1002/anie.201410407.
(2) Rasmussen, D. R.; Radom, L. "Planar-Tetracoordinate Carbon in a Neutral Saturated Hydrocarbon: Theoretical Design and Characterization," Angew. Chem. Int. Ed. 199938, 2875-2878, DOI:10.1002/(SICI)1521-3773(19991004)38:19<2875::AID-ANIE2875>3.0.CO;2-D.


InChIs

1: InChI=1S/C23H24/c1-7-11-3-15-9-2-10-17-5-13-8(1)14-6-18(10)22-16(9)4-12(7)20(14,22)23(22)19(11,13)21(15,17)23/h7-18H,1-6H2
InChIKey=LMDPKFRIIOUORN-UHFFFAOYSA-N
2: InChI=1S/C5H4/c1-2-5(1)3-4-5/h1-4H/q+2
InChIKey=UGGTXIMRHSZRSQ-UHFFFAOYSA-N




This work is licensed under a Creative Commons Attribution-NoDerivs 3.0 Unported License.

Wednesday, August 19, 2015

Can Density Cumulant Functional Theory Describe Static Correlation Effects?

Mullinax, J. W.; Sokolov, A. Y.; Schaefer, H. F. J. Chem. Theor. Comput. 2015, 11, 2487-2495 
Contributed by Steven Bachrach.
Reposted from Computational Organic Chemistry with permission

I want to update my discussion of m-benzyne, which I present in my book in Chapter 5.5.3. The interesting question concerning m-benzyne concerns its structure: is it a single ring structure 1a or a bicyclic structure 1b? Single configuration methods including closed-shell DFT methods predict the bicylic structure, but multi-configuration methods and unrestricted DFT predict it to be 1a. Experiments support the single ring structure 1a.

The key measurement distinguishing these two structure type is the C1-C3 distance. Table 1 updates Table 5.11 from my book with the computed value of this distance using some new methods. In particular, the state-specific multireference coupled cluster Mk-MRCCSD method1 with the cc-pCVTZ basis set indicates a distance of 2.014 Å.2 The density cumulant functional theory3 ODC-124 with the cc-pCVTZ basis set also predicts the single ring structure with a distance of 2.101 Å.5

Table 1. C1-C3 distance (Å) with different computational methods using the cc-pCVTZ basis set
method
r(C1-C3)
CCSD5
1.556
CCSD(T)5
2.043
OCD-125
2.101
Mk-MRCCSD2
2.014

References

(1) Evangelista, F. A.; Allen, W. D.; Schaefer III, H. F. "Coupling term derivation and general implementation of state-specific multireference coupled cluster theories," J. Chem. Phys 2007127, 024102-024117, DOI:10.1063/1.2743014.
(2) Jagau, T.-C.; Prochnow, E.; Evangelista, F. A.; Gauss, J. "Analytic gradients for Mukherjee’s multireference coupled-cluster method using two-configurational self-consistent-field orbitals," J. Chem. Phys. 2010132, 144110, DOI: 10.1063/1.3370847.
(3) Kutzelnigg, W. "Density-cumulant functional theory," J. Chem. Phys. 2006125, 171101, DOI:10.1063/1.2387955.
(4) Sokolov, A. Y.; Schaefer, H. F. "Orbital-optimized density cumulant functional theory," J. Chem. Phys.2013139, 204110, DOI: 10.1063/1.4833138.
(5) Mullinax, J. W.; Sokolov, A. Y.; Schaefer, H. F. "Can Density Cumulant Functional Theory Describe Static Correlation Effects?," J. Chem. Theor. Comput. 201511, 2487-2495, DOI: 10.1021/acs.jctc.5b00346.

InChIs

1a: InChI=1S/C6H4/c1-2-4-6-5-3-1/h1-3,6H



This work is licensed under a Creative Commons Attribution-NoDerivs 3.0 Unported License.

Friday, August 14, 2015

Domino Tunneling

Schreiner, P. R.; Wagner, J. P.; Reisenauer, H. P.; Gerbig, D.; Ley, D.; Sarka, J.; Császár, A. G.; Vaughn, A.; Allen, W. D. J. Am. Chem. Soc. 2015, 137, 7828-7834
Contributed by Steven Bachrach.
Reposted from Computational Organic Chemistry with permission

A 2013 study of oxalic acid 1 failed to uncover any tunneling between its conformations,1 despite observation of tunneling in other carboxylic acids (see this post). This was rationalized by computations which suggested rather high rearrangement barriers. Schreiner, Csaszar, and Allen have now re-examined oxalic acid using both experiments and computations and find what they call domino tunneling.2


First, they determined the structures of the three conformations of 1 along with the two transition states interconnecting them using the focal point method. These geometries and relative energies are shown in Figure 1. The barrier for the two rearrangement steps are smaller than previous computations suggest, which suggests that tunneling may be possible.

1tTt
(0.0)

TS1
(9.7)

1cTt
(-1.4)

TS2
(9.0)

1cTc
(-4.0)
Figure 1. Geometries of the conformers of 1 and the TS for rearrangement and relative energies (kcal mol-1)

Placing oxalic acid in a neon matrix at 3 K and then exposing it to IR radiation populates the excited 1tTtconformation. This state then decays to both 1cTt and 1cTc, which can only happen through a tunneling process at this very cold temperature. Kinetic analysis indicates that there are two different rates for decay from both 1tTt and 1cTc, with the two rates associated with different types of sites within the matrix.

The intrinsic reaction paths for the two rearrangements: 1tTt → 1cTt and → 1cTc were obtained at MP2/aug-cc-pVTZ. Numerical integration and the WKB method yield similar half-lives: about 42 h for the first reaction and 23 h for the second reaction. These match up very well with the experimental half-lives from the fast matrix sites of 43 ± 4 h and 30 ± 20 h, respectively. Thus, the two steps take place sequentially via tunneling, like dominos falling over.


References

(1) Olbert-Majkut, A.; Ahokas, J.; Pettersson, M.; Lundell, J. "Visible Light-Driven Chemistry of Oxalic Acid in Solid Argon, Probed by Raman Spectroscopy," J. Phys. Chem. A 2013117, 1492-1502, DOI:10.1021/jp311749z.
(2) Schreiner, P. R.; Wagner, J. P.; Reisenauer, H. P.; Gerbig, D.; Ley, D.; Sarka, J.; Császár, A. G.; Vaughn, A.; Allen, W. D. "Domino Tunneling," J. Am. Chem. Soc. 2015137, 7828-7834, DOI:10.1021/jacs.5b03322.


InChIs

1: InChI=1S/C2H2O4/c3-1(4)2(5)6/h(H,3,4)(H,5,6)
InChIKey=MUBZPKHOEPUJKR-UHFFFAOYSA-N




This work is licensed under a Creative Commons Attribution-NoDerivs 3.0 Unported License.

Friday, August 7, 2015

A Molecular Nanotube with Three-Dimensional π-Conjugation

Neuhaus, P.; Cnossen, A.; Gong, J. Q.; Herz, L. M.; Anderson, H. L. Angew. Chem. Int. Ed. 2015, 54, 7344-7348
Contributed by Steven Bachrach.
Reposted from Computational Organic Chemistry with permission

I don’t really have anything to say about this recent paper by Anderson, et al.1 They have simply prepared a very beautiful structure, an aryllated analogue of 1. They even optimized the structure of 1 at BLYP/6-31G(d) and it’s shown in Figure 1. That must have taken some time!

Figure 1. BLYP/6-31G(d) optimized structure of 1.
(Remember that you can manipulate this structure by simply clicking on in, which will launch the JMol app.)


References

(1) Neuhaus, P.; Cnossen, A.; Gong, J. Q.; Herz, L. M.; Anderson, H. L. "A Molecular Nanotube with Three-Dimensional π-Conjugation," Angew. Chem. Int. Ed. 201554, 7344-7348, DOI:10.1002/anie.201502735.



This work is licensed under a Creative Commons Attribution-NoDerivs 3.0 Unported License.

Benzene-Fused Azacorannulene Bearing an Internal Nitrogen Atom

Ito, S.; Tokimaru, Y.; Nozaki, K. Angew. Chem. Int. Ed. 2015, 54, 7256-7260
Contributed by Steven Bachrach.
Reposted from Computational Organic Chemistry with permission

A heterosubstituted corranulene analogue has now been prepared. Ito, Tokimaru, and Nozaki report the synthesis of 1 and compare it with corranulene.1 The x-ray structure of 1 shows it to be a deeper bowl than corranulene, and the bond distances suggest the Kekule structure with a central pyrrole and five Clar-type phenyl rings.
The B3LYP/6-311+G(2d,p) optimized structure of 2, then analogue of 1 missing the t-butyl group, is shown in Figure 1. Its geometry is very similar to that of 1 observed in the crystal structure. The NICS(0) values are shown in Scheme 1. These values support the notion of a central (aromatic) pyrrole surrounded by a periphery of five aromatic phenyl rings.

Scheme 1. NICS(0) values
An interesting feature of bowl compounds is their inversion. The inversion barrier, through the planar TS shown in Figure 2, is computed to be 17.0 kcal mol-1 at B3LYP/6-311+G(2d,p). This is 6-7 kcal mol-1larger than the inversion barrier of corranulene, which is not surprising given the additional phenyl groups about the periphery.

2

bowl inversion TS
Figure 1. B3LYP/6-311+G(2d,p) optimized geometry of 2.

References

(1) Ito, S.; Tokimaru, Y.; Nozaki, K. "Benzene-Fused Azacorannulene Bearing an Internal Nitrogen Atom,"Angew. Chem. Int. Ed. 201554, 7256-7260, DOI: 10.1002/anie.201502599.

InChI

1: InChI=1S/C38H23N/c1-38(2,3)18-16-27-25-14-6-12-23-21-10-4-8-19-20-9-5-11-22-24-13-7-15-26-28(17-18)35(27)39-36(31(23)25)33(29(19)21)34(30(20)22)37(39)32(24)26/h4-17H,1-3H3
InChIKey= JJPREOOFFSVARV-UHFFFAOYSA-N
2: InChI=1S/C34H15N/c1-6-16-17-7-2-9-19-21-11-4-13-23-25-15-5-14-24-22-12-3-10-20-18(8-1)26(16)30-31(27(17)19)34(29(21)23)35(32(24)25)33(30)28(20)22/h1-15H
InChIKey= FFFUKUDWQMOFQQ-UHFFFAOYSA-N




This work is licensed under a Creative Commons Attribution-NoDerivs 3.0 Unported License.

Wednesday, July 22, 2015

Simulations of Chemical Reactions with the Frozen Domain Formulation of the Fragment Molecular Orbital Method

Hiroya Nakata, Dmitri G. Fedorov, Takeshi Nagata, Kazuo Kitaura, and Shinichiro Nakamura. Journal of Chemical Theory and Computation 2015, 11, 3053−3064
Contributed by +Jan Jensen

Reprinted (adapted) with permission from J. Chem. Theory Comput. 2015, 11, 3053−3064. Copyright (2015) American Chemical Society

I am highlighting this paper because it is the first paper I have come across that reports the optimized TS, free energy barrier, and IRC for an enzymatic reaction using an all-QM description of the entire enzyme.

The underlying method is the frozen domain with dimers (FDD) formulation of the fragment molecular orbital method (FMO) implemented in the GAMES program.  The FMO is fragmentation method and in the FDD formulation the system is divided into three domain: A) an active domain where the both the geometry and density are allowed to change, B) a polarizable buffer where only the density is allowed to change and F) a frozen domain where neither is allowed to change.  The latter restraints means that the fragments in the frozen domain only have to be calculated once.

This paper presents the implementation of the fully analytic gradient and Hessian at the RHF level and the interface to the IRC methodology in GAMESS.  (The interface to the TS-optimization code was reported previously using a different variant of FMO.)  In this approach only the Hessian of the active domain is computed.

The method is then applied to the keo-enol tautomerization of phosphoglycolohydroxamic acid (PGH) by triosephosphate isomerase (TIM).  The active domain is PGH and the proton donor Glu410 while the polarizable buffer is any residue or water molecule within 5.2 Å of the active domain, which were both treated at the RHF/6-31G(d) + dispersion level of theory.  The frozen domain, i.e. the rest of the TIM monomer, is treated at the RHF/STO-3G level of theory and the entire system is comprised of 9227 atoms.

The workflow was at follows using 128 2.93 GHz Xeon cores:

1. Constrained optimization to obtain TS guess (12 hrs)
2. Hessian calculation (24 hrs)
3. TS optimization (12 hrs)
4. Hessian calculation (24 hrs)
5. IRC calculations (16 points) in each direction (42 hrs)
6. Optimization of reactants and products starting from last IRC points (12 hrs each)
7. Hessian calculation for reactant and products (24 hrs each)

Of course, at this level of theory and with such a small active domain the resulting free energy barrier (21.4 kcal/mol) is considerably different from the experimental value (14 kcal/mol) but this is only a proof-of-concept study.

An important next step for FMO studies of enzyme catalysis will be the implementation of a mixed DFT/RHF-D model where all or part of the active domain can be treated with DFT.  The computational cost-increase associated with using DFT could be off-set by using, for example, HF-3c for the non-reactive part of of the active region as well as the buffer region since this method gives excellent structures as low computational cost.

I thank +Casper Steinmann for alerting me to this article


This work is licensed under a Creative Commons Attribution 4.0