Wednesday, March 14, 2018

DeePCG: A Deep Neutral Network Molecular Force Field

DeePCG: constructing coarse-grained models via deep neural networks. L Zhang, J Han, H Wang, R Car, Weinan E. arXiv:1802.08549v2 [physics.chem-ph]
Contributed by Jesper Madsen

The idea of “learning” a molecular force field (FF) using neutral networks can be traced back to Blank et al. in 1995.[1] Modern variations (reviewed recently by Behler[2]), such as the DeePCG scheme[3] that I highlight here, seem to have two key innovations to set them apart from earlier work: network depth and atomic environment descriptors. The latter was the topic of my recent highlight and Zhang et al.[3] take advantage of similar ideas.
Figure 1: “Schematic plot of the neural network input for the environment of CG particle i, using water as an example. Red and white balls represent the oxygen and the hydrogen atoms of the microscopic system, respectively. Purple balls denote CG particles, which, in our example, are centered at the positions of the oxygens.)” from ref. [3]    
Zhang et al. simulate liquid water using ab initio molecular dynamics (AIMD) on the DFT/PBE0 level of theory in order to train a coarse-grained (CG) molecular water model. The training is done by a standard protocol used in CGing where mean forces are fitted by minimizing a loss-function (the natural choice is the residual sum of squares) over the sampled configurations. CGing liquid water is difficult because of the necessity of many-body contributions to interactions, especially so upon integrating out degrees-of-freedom. One would therefore expect that a FF capable of capturing such many-body effects to perform well, just as DeePCG does, and I think this is a very nice example of exactly how much can be gained by using faithful representations of atomic neighborhoods instead of radially symmetric pair potentials. Recall that traditional force-matching, while provably exact in the limit of the complete many-body expansion,[4] still shows non-negligible deviations from the target distributions for most simple liquids when standard approximations are used.

FF transferability, however, is likely where the current grand challenge is to be found. Zhang et al. remark that it would be convenient to have an accurate yet cheap (e.g., CG) model for describing phase transitions in water. They do not attempt this in the current preprint paper, but I suspect that it is not *that* easy to make a decent CG model that can correctly get subtle long-range correlations right at various densities, let alone different phases of water and ice, coexistences, interfaces, impurities (non-water moieties), etc. Machine-learnt potentials continuously demonstrate excellent accuracy over the parameterization space of states or configurations, but for transferability and extrapolations, we are still waiting to see how far they can get.


[1] Neural network models of potential energy surfaces. TB Blank, SD Brown, AW Calhoun, DJ Doren. J Chem Phys 103, 4129 (1995)
[2] Perspective: Machine learning potentials for atomistic simulations. J Behler. J Chem Phys 145, 170901 (2016)
[3] DeePCG: constructing coarse-grained models via deep neural networks. L Zhang, J Han, H Wang, R Car, Weinan E. arXiv:1802.08549v2 [physics.chem-ph]
[4] The multiscale coarse-graining method. I. A rigorous bridge between atomistic and coarse-grained models. WG Noid, J-W Chu, GS Ayton, V Krishna, S Izvekov, GA Voth, A Das, HC Andersen. J Chem Phys 128, 244114 (2008)

Monday, March 12, 2018

Comprehensive theoretical study of all 1812 C60 isomers

Sure, R.; Hansen, A.; Schwerdtfeger, P.; Grimme, S., Phys. Chem. Chem. Phys. 2017, 19, 14296
Contributed by Steven Bacharach
Reposted from Computational Organic Chemistry with permission

The Grimme group has examined all 1812 C60 isomers, in part to benchmark some computational methods.1 They computed all of these structures at PW6B95-D3/def2-QZVP//PBE-D3/def2-TZVP. The lowest energy structure is the expected fullerene 1 and the highest energy structure is the nanorod 2 (see Figure 1).


Figure 1. Optimized structures of the lowest (1) and highest (2) energy C60 isomers.

About 70% of the isomers like in the range of 150-250 kcal mol-1 above the fullerene 1, and the highest energy isomer 2 lies 549.1 kcal mol-1 above 1. To benchmark some computational methods, they selected the five lowest energy isomers and five other isomers with higher energy to serve as a new database (C60ISO), with energies computed at DLPNO-CCSD(T)/CBS*. The mean absolute deviation of the PBE-D3/def2-TZVP relative energies with the DLPNO-CCSD(T)/CBS* energies is relative large 10.7 kcal mol-1. However, the PW6B95-D3/def2-QZVP//PBE-D3/def2-TZVP method is considerably better, with a MAD of only 1.7 kcal mol-1. This is clearly a reasonable compromise method for fullerene-like systems, balancing accuracy with computational time.

They also compared the relative energies of all 1812 isomers computed at PW6B95-D3/def2-QZVP//PBE-D3/def2-TZVP with a number of semi-empirical methods. The best results are with the DFTB-D3 method, with an MAD of 5.3 kcal mol-1.


1) Sure, R.; Hansen, A.; Schwerdtfeger, P.; Grimme, S., "Comprehensive theoretical study of all 1812 C60isomers." Phys. Chem. Chem. Phys. 201719, 14296-14305, DOI: 10.1039/C7CP00735C.


1: InChI=1S/C60/c1-2-5-6-3(1)8-12-10-4(1)9-11-7(2)17-21-13(5)23-24-14(6)22-18(8)28-20(12)30-26-16(10)15(9)25-29-19(11)27(17)37-41-31(21)33(23)43-44-34(24)32(22)42-38(28)48-40(30)46-36(26)35(25)45-39(29)47(37)55-49(41)51(43)57-52(44)50(42)56(48)59-54(46)53(45)58(55)60(57)59
2: InChI=1S/C60/c1-11-12-2-21(1)31-41-32-22(1)3-13(11)15-5-24(3)34-43(32)53-55-47-36-26-6-16-17-7(26)28-9-19(17)20-10-29-8(18(16)20)27(6)37-46(36)54(51(41)55)52-42(31)33-23(2)4(14(12)15)25(5)35-44(33)58-56(52)48(37)39(29)50-40(30(9)10)49(38(28)47)57(53)59(45(34)35)60(50)58

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

Wednesday, February 28, 2018

Automated Transition State Theory Calculations for High-Throughput Kinetics

Pierre L. Bhoorasingh, Belinda L. Slakman, Fariba Seyedzadeh Khanshan, Jason Y. Cain, and Richard H. West (2017)
Highlighted by Jan Jensen

Figure 1 from Bhoorasingh et al. J. Phys. Chem. A 2017, 121, 6896. 
Copyright 2017 American Chemical Society

I have written about automated transition state searching before, so I was interested to see how this work differed. Both methods aim at obtaining the best possible guess of the TS structure, which is then used as a starting point for a conventionional TS optimization. In the current work this is done by estimating bond lengths between the reacting atoms using a group contribution method based on known TS structures. These distances are then constrained while a conformational search is performed for the rest of the molecular structure using the UFF force field. The method is described in more detail here.

This approach is thus not too different from the TS template structure approach used in the Schrödinger study, but goes on to perform a conformational search for the TS, which the Schrödinger study did not. So it indeed encouraging to see that the conformational search seems to work and give reasonable results.

Both approaches requires that the atom orders are the same in the reactants and products. In general this is a hard problem and the Schrödinger paper offers one approach to this. However, in the current study the products are automatically generated from the reactants using the Reaction Mechanism Generator (RMG) program in such a way (I believe) that the atom order is preserved.

So if you're interested in a particular TS the current approach is unlikely to be useful since it is rather intimately tied to the RMG program and certain types of chemical reactions. However, if you are interested in these types of chemical reactions then the approach seems quite useful since the entire process is automated and appears quite robust. 

More importantly is an important proof-of-concept of what is possible in terms of automation given a large an carefully constructed training set of chemical reactions.

This work is licensed under a Creative Commons Attribution 4.0 International License.

Wednesday, February 21, 2018

Kinetics of the Strain-Promoted Oxidation-Controlled Cycloalkyne-1,2-quinone Cycloaddition: Experimental and Theoretical Studies

Escorihuela, J.; Das, A.; Looijen, W. J. E.; van Delft, F. L.; Aquino, A. J. A.; Lischka, H.; Zuilhof, H., J. Org. Chem. 2018, 83, 244
Contributed by Steven Bacharach
Reposted from Computational Organic Chemistry with permission

Click chemistry has been used in a broad range of applications. The use of metal catalysts has limited its application to biological system, but the development of strain-promoted cycloaddition to cyclooctyne has opened up click chemistry to bioorthogonal labelling.

An interesting variation on this is the use of 1,2-benzoquinone 1 and substituted analogues as the Diels-Alder diene component. Escorihuela and co-workers have reported on the use of this diene with a number of cyclooctyne derivatives, measuring kinetics and also using computations to assess the mechanism.1

Their computations focused on two reactions using cyclooctyne 2 and the cyclopropane-fused analogue 3:
Reaction 1
Reaction 2
They examined these reactions with a variety of density functionals along with some post-HF methods. The transition states of the two reactions are shown in Figure 1. A variety of different density functionals and MP2 are consistent in finding synchronous or nearly synchronous transition states.


Figure 1. B97D/6-311+G(d,p) transition states for Reactions 1 and 2.

In terms of activation energies, all of the DFT methods consistently overestimate the barrier by about 5-10 kcal mol-1, with B97D-D3 doing the best. MP2 drastically underestimates the barriers, though the SOS-MP2 or SCS-MP2 improve the estimate. Both CCSD(T) and MR-AQCC provide estimates of about 8.5 kcal mol-1, still 3-4 kcal mol-1 too high. The agreement between CCSD(T), a single reference method, and MR-AQCC, a multireference method, indicate that the transition states have little multireference character. Given the reasonable estimate of the barrier afforded by B97D-D3, and its tremendous performance advantage over SCS-MP2, CCSD(T) and MR-AQCC, this is the preferred method (at least with current technology) for examining Diels-Alder reactions like these, especially with larger molecules.


1) Escorihuela, J.; Das, A.; Looijen, W. J. E.; van Delft, F. L.; Aquino, A. J. A.; Lischka, H.; Zuilhof, H., "Kinetics of the Strain-Promoted Oxidation-Controlled Cycloalkyne-1,2-quinone Cycloaddition: Experimental and Theoretical Studies." J. Org. Chem. 201883, 244-252, DOI: 10.1021/acs.joc.7b02614.


1: InChI=1S/C6H4O2/c7-5-3-1-2-4-6(5)8/h1-4H
2: InChI=1S/C8H12/c1-2-4-6-8-7-5-3-1/h1-6H2
3: InChI=1S/C9H12/c1-2-4-6-9-7-8(9)5-3-1/h8-9H,3-7H2
4: InChI=1S/C14H16O2/c15-13-11-7-8-12(14(13)16)10-6-4-2-1-3-5-9(10)11/h7-8,11-12H,1-6H2
5: InChI=1S/C15H16O2/c16-14-12-5-6-13(15(14)17)11-4-2-9-7-8(9)1-3-10(11)12/h5-6,8-9,12-13H,1-4,7H2/t8-,9+,12?,13?

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

Tuesday, February 6, 2018

Fully Automated Quantum-Chemistry-Based Computation of Spin–Spin-Coupled Nuclear Magnetic Resonance Spectra

Grimme, S.; Bannwarth, C.; Dohm, S.; Hansen, A.; Pisarek, J.; Pracht, P.; Seibert, J.; Neese, F., Angew. Chem. Int. Ed. 2017, 56, 14763-14769
Contributed by Steven Bacharach
Reposted from Computational Organic Chemistry with permission

Computed NMR spectra have become a very useful tool in identifying chemical structures. I have blogged on this multiple times. A recent trend has been the development of computational procedures that lead to computed spectra (again, see that above link). Now, Grimme, Neese and coworkers have offered their approach to computed NMR spectra, including spin-spin splitting.1
Their procedure involves four distinct steps.
  1. Generation of the conformer and rotamer space. This is a critical distinctive element of their method in that they take a number of different tacks for sampling conformational space to insure that they have identified all low-energy structures. This involves a combination of normal mode following, genetic structure crossing (based on genetic algorithms for optimization), and molecular dynamics. Making this all work is their choice of using the computational efficient GFN-xTB2 quantum mechanical method.
  2. The low-energy structures are then subjected to re-optimization at PBEh-3c and then single-point energies obtained at DSD-BLYP-D3/def2-TZVPP including treatment of solvation by COSMO-RS. The low-energy structures that contribute 4% or more of the Boltzmann-weighted population are then carried forward.
  3. Chemical shifts and spin-spin coupling constants are then computed with the PBE0 method and the pcS and pcJ basis sets developed by Jensen for computing NMR shifts.3
  4. Lastly, the chemical shifts and coupling constants are averaged and the spin Hamiltonian is solved.
The paper provides a number of examples of the application of the methodology, all with quite good success. The computer codes to run this method are available for academic use from


1) Grimme, S.; Bannwarth, C.; Dohm, S.; Hansen, A.; Pisarek, J.; Pracht, P.; Seibert, J.; Neese, F., "Fully Automated Quantum-Chemistry-Based Computation of Spin–Spin-Coupled Nuclear Magnetic Resonance Spectra." Angew. Chem. Int. Ed. 201756, 14763-14769, DOI: 10.1002/anie.201708266.
2) Grimme, S.; Bannwarth, C.; Shushkov, P., "A Robust and Accurate Tight-Binding Quantum Chemical Method for Structures, Vibrational Frequencies, and Noncovalent Interactions of Large Molecular Systems Parametrized for All spd-Block Elements (Z = 1–86)." J. Chem. Theory Comput. 201713, 1989-2009, DOI: 10.1021/acs.jctc.7b00118.
3) Jensen, F., "Basis Set Convergence of Nuclear Magnetic Shielding Constants Calculated by Density Functional Methods." J. Chem. Theory Comput. 20084, 719-727, DOI: 10.1021/ct800013z.

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