Stupid-Simple (*nix-Specific) Sed Scripts To Get (All Current) Gaussian09 Output Files Working With aClimax

The following three snippets of Gaussian output are for an optimization and normal mode analysis of simple olde methane (CH4).

 Gaussian 03:  EM64L-G03RevE.01 11-Sep-2007
 incident light, reduced masses (AMU), force constants (mDyne/A),
 and normal coordinates:
                     1                      2                      3
                     T                      T                      T
 Frequencies --  1356.0070              1356.0070              1356.0070
 Red. masses --     1.1789                 1.1789                 1.1789
 Frc consts  --     1.2771                 1.2771                 1.2771
 IR Inten    --    14.1122                14.1122                14.1122
 Atom AN      X      Y      Z        X      Y      Z        X      Y      Z
   1   1     0.02  -0.42   0.43    -0.34  -0.13  -0.08    -0.36  -0.23  -0.23
   2   6     0.00   0.08  -0.09     0.00   0.09   0.08     0.12   0.00   0.00
 - Thermochemistry -
 Temperature   298.150 Kelvin.  Pressure   1.00000 Atm.
 Atom  1 has atomic number  1 and mass   1.00783

Continue reading “Stupid-Simple (*nix-Specific) Sed Scripts To Get (All Current) Gaussian09 Output Files Working With aClimax”

Commensurate Urea Inclusion Crystals With The Guest (E,E)‐1,4-Diiodo-1,3-Butadiene

Published in Crystal Growth & Design (Cryst. Growth Des., 2013, 13 (9), pp. 3852–3855) earlier this year. The theory work is less impressive than the successful crystal growth, with initial solid-state efforts in Crystal09 only very recently now producing good results (leaving the molecular calculations to Gaussian09 in this paper). The procedure leading to the observed crystal structure of this inclusion complex is a significant step in the direction of testing the theory proposed in Bond Alternation In Infinite Periodic Polyacetylene: Dynamical Treatment Of The Anharmonic Potential published earlier this year in J. Mol. Struct.


Caption: Two views along the ba and ca crystal axes of the (E,E)‐1,4-Diiodo-1,3-Butadiene : Urea Inclusion Complex.

Amanda F. Lashua, Tiffany M. Smith, Hegui Hu, Lihui Wei, Damian G. Allis, Michael B. Sponsler, and Bruce S. Hudson

Abstract: The urea inclusion compound (UIC) with (E,E)-1,4-diiodo-1,3-butadiene (DIBD) as a guest (DIBD:UIC) has been prepared and crystallographically characterized at 90 and 298 K as a rare example of a commensurate, fully ordered UIC. The crystal shows nearly hexagonal channels in the monoclinic space group P21/n. The DIBD guest molecules are arranged end-to-end with the nonbonding iodine atoms in the van der Waals contact. The guest structure is compared with that for DIBD at 90 K and with computations for the periodic UIC and isolated DIBD molecule.

Dipole Derivative, Polarizability Derivative, And Vibrational Polarizability Contribution Output From Gaussian09 With IOp(7/33)

For those itching for polarizability derivative orientation information and wondering where it is when you ask for it… what’s included below is a combination of a few points in one, specifically pointing out that the IOp options are not just “another part” of the Gaussian input file (with the IOp Overlays currently linked HERE).

The problem I realized after an email from Gaussian HQ was that, as was the case for the KMLYP density functional call discussed in previous posts about [18]-annulene, “opt” and “freq” keyword combinations are seen as two distinct runs in Gaussian that don’t pass the IOp information along (and, admittedly, I should have remembered that). Specifically, the additional print-out for the polarizability info is called by IOp(7/33=3).

Continue reading “Dipole Derivative, Polarizability Derivative, And Vibrational Polarizability Contribution Output From Gaussian09 With IOp(7/33)”

Bond Alternation In Infinite Periodic Polyacetylene: Dynamical Treatment Of The Anharmonic Potential

In press (DOI:10.1016/j.molstruc.2012.07.051) in the Journal Of Molecular Structure. May go down in history as a hardest-fought paper acceptance. In a similar line of research as the [18]-annulene study, but exploring the infinite limit of geometry and bond length alternation energy barrier for this infinite case. If the numbers are correct, the infinite polyene chains (polyacetylene) do not exhibit bond length alternation because the Peierls’ barrier between the single-double and double-single bond alternate minima is below the vibrational zero-point level. Plenty of ramifications.

Bruce S. Hudson and Damian G. Allis

Abstract. The potential energy of the infinite periodic chain model of polyacetylene (pPA) is symmetric with two equivalent minima separated by the Peierls’ stabilization barrier. In this work it is shown how an energy scale and vibrational energy levels for this highly anharmonic Peierls’ degree of freedom can be estimated. Attention is given to the potential energy increase for large deformations. The Born-Kármán treatment of translational symmetry is applied. Two empirical methods and a direct periodic boundary condition (PBC) density functional theory (DFT) calculations are in semi-quantitative agreement, each leading to the conclusion that pPA has a zero-point level that is above the Peierls’ barrier. The argument does not depend critically on the barrier height or the other parameters of the model or the computation method. It is concluded that pPA will not exhibit bond alternation and that the zero-point average geometry does not preclude possible conductivity.

The Structure Of [18]-Annulene: Computed Raman Spectra, Zero-Point Level And Proton NMR Chemical Shifts

In press (DOI:10.1016/j.molstruc.2012.05.016) in the Journal Of Molecular Structure (Volume 1023, 12 September 2012, Pages 212–215) in the special issue: MOLECULAR VIBRATIONS AND STRUCTURES: THEORY AND EXPERIMENT — A collection of papers dedicated to Professor Jaan Laane on the occasion of his 70th birthday.

This paper on the “actual” geometry of [18]-annulene is part of several larger stories addressing a larger polyene (or larger-polyene) issue. First among these is the meaning of experimental results obtained by various spectroscopic methods (in this case, using previous X-ray, Raman (with the C2 (blue) and D6h (red) simulated spectra shown in the image above), IR, and NMR data that produce different results within the limitations of the methods to study the single molecule). Second is the quality of the theoretical method for reproducing certain types of spectroscopic data. In the case of the [N]-annulene series, the ever-present B3LYP density functional is found to produce the time-average geometry of [18]-annulene found in X-ray data, but another density functional (in this case, KMLYP), finds that bond-alternate minima exist. Third is the importance of the zero-point level in the treatment of systems for which bond-alternate geometries exist with transition-state barriers calculated to be below the zero-point level in the classical approximation of nuclear positions (the Born-Oppenheimer Approximation).

NOTE 1: The KMYLP density functional is called in Gaussian with the following keyword set:

BLYP iop(3/76=1000005570) iop(3/77=0000004430) iop(3/78=0448010000)

NOTE 2: Optimization and Frequency calculations must be performed as TWO SEPARATE CALCULATIONS. The iop-called density functional does not carry itself over between opt + freq (or other properties) in the same input file. If you opt + freq in the same input file, you will Opt with KMLYP but freq with BLYP. This will be obvious by the number of imaginary modes.

Bruce S. Hudson and Damian G. Allis

Abstract. [18]-annulene has been of great interest from the structural point of view of its bond alternation. High-level calculations based on structures selected for agreement with NMR spectra lead to a bond-alternate C2 form over a non-alternating planar D6h structure deduced from diffraction, infrared (IR) and electronic spectral studies. Here it is shown that computed Raman spectra for the D6h and C2 forms are expected to be very different. However, two equivalent non-D6h bond-alternate minima of D3h or C2 geometries are separated by only a small barrier along a motion that involves CC stretching and compression. It is shown here that the zero-point level is above the barrier for this species. In light of that fact, the NMR calculations are reconsidered with inclusion of zero-point level averaging.

CCSD(T) and MP4 Z-Matrix Coordinate Oddity In Gaussian: When In Doubt, Change Your Format

This post was instigated by Syracuse University Professor of Chemistry and well-known non-blogger Tim Korter concerning efforts to, I believe, generate proper Møller–Plesset Perturbation Theory Of The 4th Order (MP4, and also testing coupled cluster CCSD(T) calculations) intermolecular potentials for improving terms for Grimme dispersion-corrected density functional theory (DFT) calculations with the Gaussian09 package (a program for which many people grumble about various issues but which is, by nearly all metrics, a fantastic set of quantum chemical programs). The examples below, using water only, are just for ease-of-testing, which produce the following results based on the form of the input of the molecular coordinates. For those wondering why, z-matrices are the preferred format for performing SCAN or other automated trajectory calculations (an absolutely useless format, in my opinion, now that we have computers that can handle more than five atoms).

Continue reading “CCSD(T) and MP4 Z-Matrix Coordinate Oddity In Gaussian: When In Doubt, Change Your Format”

/var/lib/ureadahead/debugfs And Disk Space “Recovery” In Ubuntu 10.04

And happy new year.

I had thought this was something involving Gaussian09 memory usage until I restarted a machine and found the same problem occurring in Ubuntu. Below is a quick fix (and reminder for the future).

Checking disk space with df -all :

user@machine:~$ df -all
Filesystem           1K-blocks      Used Available Use% Mounted on
/dev/sda1            147550696   5863896 134191632   5% /
proc                         0         0         0   -  /proc
none                         0         0         0   -  /sys
none                         0         0         0   -  /sys/fs/fuse/connections
none                         0         0         0   -  /sys/kernel/debug
none                         0         0         0   -  /sys/kernel/security
none                   3053752       260   3053492   1% /dev
none                         0         0         0   -  /dev/pts
none                   3058264         0   3058264   0% /dev/shm
none                   3058264        84   3058180   1% /var/run
none                   3058264         0   3058264   0% /var/lock
none                   3058264         0   3058264   0% /lib/init/rw
none                         0         0         0   -  /var/lib/ureadahead/debugfs
nfsd                         0         0         0   -  /proc/fs/nfsd
binfmt_misc                  0         0         0   -  /proc/sys/fs/binfmt_misc

Continue reading “/var/lib/ureadahead/debugfs And Disk Space “Recovery” In Ubuntu 10.04″

The Vibrational Spectrum Of Parabanic Acid By Inelastic Neutron Scattering Spectroscopy And Simulation By Solid-State DFT

Available as an ASAP in The Journal of Physical Chemistry A. As a general rule in computational chemistry, the smaller the molecule, the harder it is to get right. As a brief summary, parabanic acid has several interesting properties of significance to computational chemists as both a model for other systems containing similar sub-structures and as a complicated little molecule in its own right.

1. The solid-state spectrum requires solid-state modeling. This should be of no surprise (see the figure below for the difference in solid-state (top) and isolated-molecule (bottom)). This task was undertaken with both DMol3 and Crystal06, with DMol3 calculations responsible for the majority of the analysis of this system (as has always been the case in the neutron studies reported on this site).

2. The agreement in the hydrogen-bonded N-H…O vibrations is, starting from the crystal structure, in poor agreement with experiment. You’ll note the region between 750 and 900 cm-1 is a little too high (and for clarification, the simulated spectrum is in red below). According to the kitchen sink that Matt threw at the structure, the problem is not the same anharmonicity one would acknowledge by Dr. Walnut’s “catalytic handwaving” approach to spectrum assignment (Dr. Walnut does not engage in this behavior, rather endeavors to find it in others where it should not be).

3. The local geometry of the hydrogen-bonding network in this molecular solid leads to notable changes in parabanic acid structure that, in turn, leads to the different behavior of the N-H…O vibrational motions. There is one potentially inflammatory comment in the Conclusions section that results from this identification. The parabanic acid molecule is, at its sub-structure, a set of three constrained peptide linkages that under go subtle but vibrationally-observable changes to their geometry because of crystal packing and intermolecular hydrogen bond formation. This means that the isolated molecule and solid-state forms are different and that peptide groups are influenced by neighboring interactions.

So, why should one care? Suppose one is parameterizing a biomolecular force field (CHARMM, AMBER, GROMOS, etc.) using bond lengths, bond angles, etc., for the amino acid geometry and vibrational data for some aspect of the force constant analysis. The structural data for these force fields often originates with solid-state studies (diffraction results). This means, to those very concerned with structural accuracy, that a geometry we know to be influenced by solid-state interactions is being used as the basis for molecular dynamics calculations that will NOT be used in their solid-state forms. Coupled with the different spectral properties due to intermolecular interactions, the description being used as the basis for the biomolecular force field likely being used in solution (solvent box approaches) is based on data in a phase where the structure and dynamics are altered from their less conformationally-restricted counterpart (in this case, solid-state).

A subtle point, but that’s where applied theoreticians do some of their best work.

Matthew R. Hudson, Damian G. Allis, and Bruce S. Hudson

Department of Chemistry, 1-014 Center for Science and Technology, Syracuse University, Syracuse, New York 13244-4100

Abstract: The incoherent inelastic neutron scattering spectrum of parabanic acid was measured and simulated using solid-state density functional theory (DFT). This molecule was previously the subject of low-temperature X-ray and neutron diffraction studies. While the simulated spectra from several density functionals account for relative intensities and factor group splitting regardless of functional choice, the hydrogen-bending vibrational energies for the out-of-plane modes are poorly described by all methods. The disagreement between calculated and observed out-of-plane hydrogen bending mode energies is examined along with geometry optimization differences of bond lengths, bond angles, and hydrogen-bonding interactions for different functionals. Neutron diffraction suggests nearly symmetric hydrogen atom positions in the crystalline solid for both heavy-atom and N-H bond distances but different hydrogen-bonding angles. The spectroscopic results suggest a significant factor group splitting for the out-of-plane bending motions associated with the hydrogen atoms (N-H) for both the symmetric and asymmetric bending modes, as is also supported by DFT simulations. The differences between the quality of the crystallographic and spectroscopic simulations by isolated-molecule DFT, cluster-based DFT (that account for only the hydrogen-bonding interactions around a single molecule), and solid-state DFT are considered in detail, with parabanic acid serving as an excellent case study due to its small size and the availability of high-quality structure data. These calculations show that hydrogen bonding results in a change in the bond distances and bond angles of parabanic acid from the free molecule values.

Exploring the Implications of Vitamin B12 Conjugation to Insulin on Insulin Receptor Binding and Cellular Uptake

In press, in the journal ChemMedChem (and, because I think it’s hip, I note that the current “obligatory” image for the wikipedia article for ChemMedChem features the image I made for the review article on the topic addressed in this new study). As with many theory papers (there’s some experiment in there, too), this very brief article summarizes several months of cyanocobalamin (B12) parameterization and molecular dynamics (MD) simulations. The purpose of the theory was to address all of the major structural snapshots in the uptake process associated with the insulin-B12 bioconjugate being developed as part of the much heralded oral insulin project in Robert Doyle’s group here at Syracuse. These structures include:

1. The structure and dynamic properties of the insulin-B12 bioconjugate
2. The binding of B12 to Transcobalamin II (TCII) (for B12 parameterization)
3. The binding of the insulin-B12 bioconjugate to TCII (and the steric demands therein)
4. The interaction of the insulin-B12 bioconjugate, bound to TCII, with the insulin Receptor (IR)

The quantum chemical (for the B12 geometry and missing force constants) and molecular dynamics (GROMACS with the GROMOS96 (53a6)) simulation work is going to serve as the basis for several posts here (eventually) about parameterization, topology generation, and force field development.

As an example of some of the insights modeling provides, the figure above shows the insulin-B12 bioconjugate (the insulin is divided into A and B chains, the A chain in blue and the important division of the insulin B chain in the front half of the rainbow). Insulin is a rather large-scale example of many of the same molecular issues that arise in the analysis of solid-state molecular crystals by either terahertz or inelastic neutron scattering spectroscopy. The packing of molecules in their crystal lattices can lead to significant changes in molecular geometry, be these changes in the stabilization of higher-energy molecular conformations or even deformations in the covalent framework. In the case of insulin, it is found that the crystal geometry (also the geometry of stored insulin in the body) is quite different from the solution-phase form. It’s even worse! The B chain end (B20-B30) in the solid-state geometry covers (protects?) the business-end of the insulin binding region to the Insulin Receptor. One can imagine the difficulty in proposing the original binding model for insulin to its receptor from the original crystal data given that the actual binding region is blocked off in the solid-state form! The “Extended” form in the figure is representative of “multiple other” conformations of the B20-B30 region (which mimics the characterized T-state of insulin), those geometries for which the insulin binding region (blue and green) is completely exposed. This extended geometry is also the one that separates the bulk of the insulin structure from the covalently-linked B12 (at Lys29) and, it is argued from the MD simulations in the paper, enables the B12 to still tightly bind to TCII despite the presence of all this steric bulk.

Amanda K. Petrus1, Damian G. Allis1, Robert P. Smith2, Timothy J. Fairchild3 and Robert P. Doyle1

1. Department of Chemistry, Syracuse University, Syracuse, NY 13244, USA
2. Department of Construction Management and Wood Products Engineering, SUNY, College of Environmental Science and Forestry, Syracuse, NY 13210, USA
3. Department of Exercise Science, Syracuse University, Syracuse, NY 13244, USA

Extract: We recently reported a vitamin B12 (B12) based insulin conjugate that produced significantly decreased blood glucose levels in diabetic STZ-rat models. The results of this study posed a fundamental question, namely what implications does B12 conjugation have on insulin’s interaction with its receptor? To explore this question we used a combination of molecular dynamics (MD) simulations and immuno-electron microscopy (IEM).

The Solid-State Terahertz Spectrum of MDMA (Ecstasy) – A Unique Test for Molecular Modeling Assignments

In press, in Chemical Physics Letters (CPL).  Yes, the blog has taken a bit of a turn, from high explosives to illicit drugs.  I expect my google rating to rise sharply with this post.  The protonated form of the 3,4-methylene-dioxymethamphetamine (Ecstasy, but I’ll keep the post legit, so it is herein referred to as MDMA) molecule (herein referred to as MDMA:H+) and MDMA:H+ in its crystal cell with a chloride ion (Cl, the crystal form herein referred to as MDMA:HCl) is shown below in yet another fantastic NanoEngineer-1 rendering (if I do say so myself).

This CPL article is, to some extent, a response to those in the terahertz community who continue to attempt spectral assignments of crystalline and poly-crystalline samples using isolated-molecule quantum chemical calculations.  The MDMA:HCl sample and MDMA molecule, as a pair, are a very interesting case study of theory and experiment for reasons detailed below.  The spectra, shown below from a previous version of the paper (but the same spectra), show quite a bit of detail that will make sense shortly.

Panel A shows the isolated-molecule calculation for the neutral MDMA molecule at a B3LYP/6-31G(d) level of theory (in red).  You will note that this simulated spectrum is in very good agreement with experiment (in black), reproducing all of the major features and showing a number of smaller features that account for shoulders.  This agreement was the basis for the assignment of the MDMA:HCl spectrum reported in: G Wang, J Shen, Y Jia. “Vibrational spectra of ketamine hydrochloride and 3,4-methylenedioxymethamphetamine in terahertz range.” Journal of Applied Physics 102 (2007) 013106/1-06/4.

The new theoretical analysis reported in the CPL article was instigated by this assignment in this previous publication.  Relevant to the measured sample and the previously reported assignment, two points arise that require address.

1. The previous calculation, as reported, was of the neutral MDMA molecule and is reasonably close to the MDMA spectrum shown above (in red.  This calculation was redone for the CPL article for comparative purposes).  As the experimental THz sample was of solid-state MDMA:HCl, the appropriate form of the molecule to run is not the neutral MDMA molecule, but the protonated form, MDMA:H+.  The protonated form has a different vibrational spectrum (shown in green in Panel A) than the isolated molecule form.  At the very least, the isolated-molecule to consider for the MDMA:HCl sample must be the protonated form.  Interestingly, the re-calculation at B3LYP/6-31G(d) reported for the CPL article predicts a fifth vibrational mode at 48.0 cm-1 that was not reported in the previous study.  We do not know if the previous group missed that peak in the write-up, decided that (since it is a low-intensity mode) it was not worth reporting, or if their starting molecular geometry was somehow different so that the other four modes were predicted to be in the same region and this mode was somehow turned off.

2. The solid-state spectrum shown in Panel B at a BP/DNP level of theory does not agree as well as the isolated-molecule MDMA B3LYP/6-31G(d) calculation. That being said, THAT IS NOT THE POINT.  The goal of a simulated spectrum IS NOT to obtain the closest spectral agreement with experiment.  The goal IS to explain the solid-state spectrum with the best theoretical model possible that, hopefully, is as close to the experimental result as possible.  In this case, the solid-state BP/DNP spectrum contains a finite number of vibrational modes that do group according to features in the THz spectrum, making the assignment reasonably straightforward.  Interestingly, the two most intense modes in the solid-state BP/DNP calculations involve the motions of the Cl…H+-N chains, which CANNOT be accounted for in an isolated-molecule calculation of either the neutral MDMA molecule or the protonated MDMA:H+.

In summary, as taken from the CPL paper:

With all of these considerations taken into account in this re-examination of the MDMA.HCl THz spectrum, it is found that this system serves as a fortuitous example of one whose THz spectrum is predicted quite precisely by two very different approaches, but is only described accurately by one that considers the crystal environment and the actual state of the molecule in its solid-state form.

Damian G. Allis1,2, Patrick M. Hakey1, and Timothy M. Korter1

1. Department of Chemistry, Syracuse University, Syracuse NY 13244-4100 USA
2. Nanorex, Inc., P.O. Box 7188, Bloomfield Hills, MI 48302-7188 USA

Abstract: The terahertz (THz, far-infrared) spectrum of 3,4-methylene-dioxymethamphetamine hydrochloride (Ecstasy) is simulated using solid-state density functional theory.  While a previously reported isolated-molecule calculation is noteworthy for the precision of its solid-state THz reproduction, the solid-state calculation predicts that the isolated-molecule modes account for only half of the spectral features in the THz region, with the remaining structure arising from lattice vibrations that cannot be predicted without solid-state molecular modeling.  The molecular origins of the internal mode contributions to the solid-state THz spectrum, as well as the proper consideration of the protonation state of the molecule, are also considered.