Category Archives: Computational Chemistry

New paper in JACS


Well, I only contributed with the theoretical section by doing electronic structure calculations, so it isn’t really a paper we can ascribe to this particular lab, however it is really nice to see my name in JACS along such a prominent researcher as Prof. Chad Mirkin from Northwestern University, in a work closely related to my area of research interest as macrocyclic recognition agents.

In this manuscript, a calix[4]arene is allosterically opened and closed reversibly by coordinating different kinds of ligands to a platinum center linked to the macrocycle. (This approach has been referred to as the weak link approach.) I recently visited Northwestern and had a great time with José Mendez-Arroyo, the first author, who showed me around and opened the possibility for further work between our research groups.

(Ligands: Green = Chloride; Blue = Cyanide)

Closed, semi-open and fully open conformations; selectivity is modulated through cavity size. (Ligands: Green = Chloride; Blue = Cyanide)

Here at UNAM we calculated the interaction energies for the two guests that were successfully inserted into the cavity: N-methyl-pyridinium (Eint = 57.4 kcal/mol) and Pyridine-N-oxide (Eint = +200.0 kcal/mol). Below you can see the electrostatic potential mapped onto the electron density isosurface for one of the adducts. Relative orientation of the hosts within the cavity follows the expected (anti-) alignment of mutual dipole moments. At this level of theory, we could easily be inclined to assert that the most stable interaction is indeed the one from the semi-open compound and that this in turn is due to the fact that host and guest are packed closer together but there is also an orbital issue: Pyridine Oxide is a better electron acceptor than N-Me-pyridinium and when we take a closer look to the (Natural Bonding) orbitals interacting it becomes evident that a closer location does not necessarily yields a stronger interaction when the electron accepting power of the ligand is weaker (which is, in my opinion, both logic and at the same time a bit counterintuitive, yet fascinating, nonetheless).

Electrostatic potential mapped onto the electron density surface of one of the aducts under study

Electrostatic potential mapped onto the electron density surface of one of the adducts under study

All calculations were performed at the B97D/LANL2DZ level of theory with the use of Gaussian09 and NBO3.1 as provided within the former. Computing time at UNAM’s supercomputer known as ‘Miztli‘ is fully acknowledged.

The full citation follows:

A Multi-State, Allosterically-Regulated Molecular Receptor With Switchable Selectivity
Jose Mendez-Arroyo Joaquín Barroso-Flores §,Alejo M. Lifschitz Amy A. Sarjeant Charlotte L. Stern , and Chad A. Mirkin *

J. Am. Chem. Soc., Article ASAP
DOI: 10.1021/ja503506a
Publication Date (Web): July 9, 2014

 Thanks to José Mendez-Arroyo for contacting me and giving me the opportunity to collaborate with his research; I’m sure this is the first of many joint projects that will mutually benefit our groups. 

 

Some .fchk files wont open in GaussView5.0 (Update)


A couple of weeks ago I posted a solution for a common error regarding .fchk files that will display the error below when opened with GaussView5.0. As I expected, this error has to do with the use of diffuse functions in the basis set and is related to a change of format between Gaussian versions.

CConnectionGFCHK::Parse_GFCHK()
Missing or bad data: Alpha Orbital Energies
Line Number 1234

Although the method described in the previous post works just fine, the following update is a better approach. Due to a change of spelling between G03 and G09 (which has been corrected for G09 but not available for GV versions prior to 5.0.9) one must change “independent” for “independant

To make the change directly from the terminal the following command is needed:

sed -i 's/independent/independant/g' file.fchk

Alternatively you can redirect the output to a new file

sed -e 's/independent/independant/g' file.fchk > newfile.fchk

if you want to keep the old version and work with a new one.

Of course this edition can be performed manually with any text editor available (for example if you work in Windows) but solutions from the terminal always seem easier and a lot more fun to me.

Thanks to Dr. Fernando Cortés for sharing his insight into this issue.

If a .fchk file wont open in GaussView5.0


I’ve found the following error regarding the opening of .fchk files in GaussView5.0.

CConnectionGFCHK::Parse_GFCHK()
Missing or bad data: Alpha Orbital Energies
Line Number 1234

The error is prevented to a first approximation (i.e. it at least will allow GV to open and visualize the file but other issues may arise) by opening the file and modifying the number of basis functions to equal the number of independent functions (which is lower)

FILE HEADER 
FOpt RM062X 6-311++G(d,p) 
Number of atoms I 75
Info1-9 I N= 9
 163 163 0 0 0 110
 2 18 -502
Charge I 0
Multiplicity I 1
Number of electrons I 314
Number of alpha electrons I 157
Number of beta electrons I 157
Number of basis functions I 1199
Number of independent functions I 1199
Number of point charges in /Mol/ I 0
Number of translation vectors I 0
Atomic numbers I N= 75
... ...
... ...

Once both numbers match you can open the file normally and work with it. My guess is this will continue to happen with highly polarized basis sets but I need to run some tests.

Virtual Conference in Computational Chemistry VCCC-2014 (1st call)


So many things have happened since I last updated this blog but I will come to write on them when appropriate. Right now I’d like to share an invitation by Prof. Ponnadurai Ramasami from the University of Mauritius to the upcoming Virtual Conference on Computational Chemistry from the 1st to the 31st of August. Deadlines can be consulted here and the most important is the abstract submission on June 30th. This conference is part of the official celebrations of the International Year of Crystallography so talks involving experimental determination of electron densities will be well suited.

I participated in the latest edition and I must say it was a very enriching opportunity to learn from so many other researchers from across the world without leaving my desk. I already know what my talk will be about, now that we are so close to finish and submit a paper on the absence of reactivity for an anti-aromatic set of molecules. (I think I’ll call it “The reactivity of molecules that never existed [but that maybe should have.].) All talks are sent either as pdf, powerpoint presentations, youtube videos, etc. and Q&A are done over e-mail.

So this is a calling to other computational chemists out there who want to participate in this virtual conference. Kudos to Prof. Ponnadurai Ramasami and lets hope we can crystallize his visit to Mexico during 2014 the International Year of Crystallography (pun intended) and here’s to me going sometime to Mauritius!

New paper in Computational and Theoretical Chemistry


I always get very happy to have a new paper out there! I find it exciting but most of all liberating since it makes you feel like your work is going somewhere but most of all that it is making its way ‘out there'; there is a strong feeling of validation, I guess.

Two very different families of calix[n]arenes (Fig 1) were tested as drug carriers for a very small molecule with a huge potential as a chemotherapeutic agent against Leukemia, namely, 3-phenyl-1H-[1]benzofuro[3,2-c]pyrazole a.k.a. GTP which has proven to be an effective in vitro Tyrosine Kinase III inhibitor. Having such a low molecular weight it is expected to have a very high excretion rate therefore the use of a carrier could increase its retention time and hence its activity. This time we considered n = 4, 5, 6 and 8 for the size of the cavities and R = -SO3H and -OEt as functional groups on the upper rim as to evaluate only a polar coordinating group and a non-polar non-coordinating one since GTP offers two H-bond acceptor sites and one H-bond donor one along the π electron density that could form π – π stacking interactions between the aromatic groups on GTP and the walls of the calixarene.

Fig 1. Calixarenes under study and their complexes with GTP

Fig 1. Calixarenes under study and their complexes with GTP

Once again calculations were carried out at the B97D/6-31G(d,p) level of theory along with molecular dynamics simulations for over 100 ns of production runs. NBO Deletion interaction energies were computed in order to discern which hosts formed the most stable complexes.

NBO Del interaction energies B97D/6-31G(d,p)

NBO Del interaction energies B97D/6-31G(d,p)

You may find a link to the ScienceDirect website for downloading the paper from this link. Last, but certainly not least, I’d like to thank all coauthors for their contributions and patience in getting this study published: Dr. Rodrigo Galindo-Murillo; Alberto Olmedo-Romero; Eduardo Cruz-Flores; Dr. Petronela M. Petrar and Prof. Dr. Kunsági-Máté Sándor. Thanks a lot for everything!

fig8

Donor and acceptor H-bond sites increases the probability of keeping the drug in place for a higher retention rate

Donor and acceptor H-bond sites increases the probability of keeping the drug in place for a higher retention rate

Internal Symposium at CCIQS – 2014 edition


Once again as every year we celebrate our internal symposium here at CCIQS, and like every year, my students presented some of their progress with their research projects. This time, three students, from three different levels, present posters regarding some of the data they’ve obtained.

20140220-221553.jpg

María Eugenia ‘Maru’ Sandoval presented a poster regarding the molecular dynamics simulations performed for the drug Imatinb and a family of calix- and thia-calix[n]arenes as published here and reported in this blog here. ‘Maru’ is now a first year grad student at the National University, UNAM, after spending a year working for a pharmaceutical company. Her research in the realm of photosynthesis has only begun recently, that is why we had to rely on some other data.

20140220-221627.jpg

Luis Enrique Aguilar is researching cation-π interactions within the aromatic cavities of calix[n]arenes in order to find suitable leads among these, our favorite macrocyles, for designing extraction agents of heavy (toxic) metals. Luis Enrique is an undergrad student here at the State University who should finish this year and has shown some interest (threatened us) in writing his dissertation thesis in our research group.

20140220-221638.jpg

Monserrat Enriquez is a PhD student at CINVESTAV under the joint supervision of Dr. Eddie López-Honorato and myself (Dr. Eddie is her principal advisor), her research project involves both theoretical calculations and synthesis of the leads for extraction agents for several Arsenic species. For the time being, Monserrat is here with us, far from her home on the north side of the country, for this semester in which we have to finish with the theoretical section of her work. Besides her research concerning calixarenes she is also running calculations on the interactions between graphene oxide and the aforementioned As species. We are very excited about working with such a complex yet simple material that has such an exciting electronic structure.

This symposium is always interesting and important in bringing our research projects closer to all the comunity of this center. And since symposium comes from the Greek meaning ‘drinking together‘, then lets raise our glasses and toast for the data to come!

Cheers!

20140220-221653.jpg

New paper in Journal of Chemical Theory and Computation


Happy new year to all my readers!

Having a new paper published is always a matter of happiness for this computational chemist but this time I’m excedingly excited about anouncing the publishing of a paper in the Journal of Chemical Theory and Computation, which is my highest ranked publication so far! It also establishes the consolidation of our research group at CCIQS as a solid and competitive group within the field of theoretical and computational chemistry. The title of our paper is “In Silico design of monomolecular drug carriers for the tyrosine kinase inhibitor drug Imatinib based on calix- and thiacalix[n]arene host molecules. A DFT and Molecular Dynamics study“.

In this article we aimed towards finding a suitable (thia-) calix[n]arene based drug delivery agent for the drug Imatinib (Gleevec by Novartis), which is a broadly used powerful Tyrosine Kinase III inhibitor used in the treatment of Chronic Myeloid Leukaemia and, to a lesser extent, Gastrointestinal Stromal Tumors; although Imatinib (IMB) exhibits a bioavailability close to 90% most of it is excreted, becomes bound to serum proteins or gets accumulated in other tissues such as the heart causing several undesired side effects which ultimately limit its use. By using a molecular capsule we can increase the molecular weight of the drug thus increasing its retention, and at the same time we can prevent Imatinib to bind, in its active form, to undesired proteins.

We suggested 36 different calix and thia-calix[n]arenes (CX) as possible candidates; IMB-CX complexes were manually docked and then optimized at the B97D/6-31G(d,p) level of theory; Stephan Grimme’s B97D functional was selected for its inclusion of dispersion terms, so important in describing π-π interactions. Intermolecular interaction energies were calculated under the Natural Bond Order approximation; a stable complex was needed but a too stable complex would never deliver its drug payload! This brings us to the next part of the study. A monomolecular drug delivery agent must be able to form a stable complex with the drug but it must also be able to release it. Molecular Dynamics simulations (+100 ns) and umbrella sampling methods were used to analyse the release of the drug into the aqueous media.

Optimized geometries for all complexes under study (B97D/6-31G*)

Optimized geometries for the 20 most stable complexes under study (B97D/6-31G*)

Potential Mean Force profiles for the four most stable complexes for position N1 and  N2 from the QM simulations are shown below (Red, complexes in the N1 position, blue, N2 position). These plots, derived from the MD simulations  give us an idea of the final destination of the drug respect of the calixarene carrier. In the next image, the three preferred structures (rotaxane-like; inside; released) for the final outcome of the delivery process are shown. The stability of the complexes was also assessed by calculating the values of ΔG binding through the use of the Poisson equations.

PMF for the most stable compounds

PMF for the most stable compounds

General MD simulation final structures

General MD simulation final structures

Thanks to my co-authors Maria Eugenia Sandoval-Salinas and Dr. Rodrigo Galindo-Murillo for their enormous contributions to this work; without their hard work and commitment to the project this paper wouldn’t have been possible.

Transition State Search (QST2 & QST3) and IRC with Gaussian09


Theoretical evaluation of a reaction mechanism is all about finding the right transition states (TS) but there are no guarantees within the available methods to actually find the one we need. Chemical intuition in the proposal of a mechanism is paramount. Let’s remember that a TS is a critical point on a Potential Energy Surface (PES) that is a minimum in every dimension but one. For a PES with more than two degrees of freedom, a hyper-surface, envisioning the location of a TS is a bit tricky, in the case of a three dimensional PES (two degrees of freedom) the saddle point constitutes the location of the TS as depicted in figure 1 by a section of a revolution hyperboloid.

400px-Saddle_point

Fig1. Saddle point on a surface (min in one direction; max in the other)

Fig 1a Pringles chips -Yuck-. They exhibit a maximum on the direction parallel to the screen and a minimum on the direction perpendicular to the screen at the same point.

Fig 1a Pringles chips -Yuck-. They exhibit a maximum on the direction parallel to the screen and a minimum on the direction perpendicular to the screen at the same point.

The following procedure considers gas phase calculations. Nevertheless, the use of the SCRF keyword activates the implicit solvent calculation of choice in order to evaluate to some degree the solvent influence on the reaction energetics at different temperatures with the use of the temperature keyword.

The first step consists of a high level optimization of all minimums involved, such as reagents, products and intermediates, with a subsequent frequency analysis that includes no imaginary eigenvalues.

In order to find the structures of the transition states we use in Gaussian the Synchronous Transit-guided Quasi-Newton method [1] through the keywords QST2 or QST3. In the former case, coordinates for the reagents and products are needed as input; for the latter keyword, coordinates for the TS structure guess is needed also.

QST2)

%chk=file.chk
%nprocshared=n
%mem=nGB

#p opt=(qst2,redundant) m062x/6-31++G(d,p) freq=noraman Temperature=373.15 SCRF=(Solvent=Water)

Title card for reagents

Q M
Cartesian Coordinates for reagents
blank line
Title card for products

Q M
Cartesian Coordinates for products
blank line

QST3)

%chk=file.chk
%nprocshared=n
%mem=nGB

#p opt=(qst3,redundant) m062x/6-31++G(d,p) freq=noraman Temperature=373.15 SCRF=(Solvent=Water)

Title Card for reagents

Q M
Cartesian Coordinates for reagents
blank line
Title card for products

Q M
Cartesian Coordinates for products
blank line–
Title card for TS
Q M
Cartesian Coordinates for TS
blank line

NOTE: It is fundamental that the numbering order is kept constant throughout the molecular specifications of all two, or three, input structures. Hence, I recommend to build a set of molecules, save their structure, and then modified the coordinates on the same file to produce the following structure, that way the number for every atom will remain the same for every step.

As I wrote above, there are no guarantees of finding the right TS so many attempts are probably needed. Once we have the optimized structures for all the species involved in our mechanistic proposal we can plot their energies very simply with MS Excel the way we’ve previously described in this blog (reblogged from eutactic.wordpress.com)

Once we’ve succeeded in finding the structure of our TS we may run an Internal Reaction Coordinate (IRC) calculation. This calculation will connect the TS structure to those of the products and the reagents. Initial constant forces are required and these are commonly retrieved from the TS calculation checkpoint file through the RCFC keyword.

%chk=QST3_2p.chk
%nprocshared=8

#p m062x/6-31++G(d,p) IRC=(Maxpoints=50,RCFC,phase=(2,1))Temperature=373.15 SCRF=(Solvent=Water) geom=allcheck

Title Card

Q M
blank line

Finally, the IRC path can be visualized with GaussView from the Results menu. A successful IRC will link both structures along a single reaction coordinate proving that both reagents and products are linked by the obtained TS.

Hat tip to Howard Diaz who has become quite skillful in calculating these mechanisms as proven by his recent poster at the XII RMFQT a couple of weeks back. And as usual thanks to everyone who reads, comments, likes, recommends, rates and shares my silly little posts.

XIIth Mexican Reunion on Theoretical Physical Chemistry


As every year this month we had the yearly Mexican Reunion on Theoretical Physical Chemistry organized by prominent researchers in the field, such as Dr. Emilio Orgaz (UNAM), Dr. Alberto Vela (CINVESTAV) and many other. Over 150 different works were presented during this edition which took place in Juriquilla, Querétaro at one of the many campuses of the National Autonomous University of Mexico scattered all around the country. Below you can see some pictures from the talks and the first poster session.

20131119-154119.jpg20131119-154102.jpg20131119-154222.jpg

This time we contributed with a small poster on a mechanism proposed by Howard Diaz (an undergrad student from UAEM) on the equilibrium transformation of dihydrocinolines into 1-amino-indoles by an intramolecular rearrangement. May this post also serve as the starting point of a -mini-tutorial on how to evaluate a mechanism theoretically using QST3 and IRC in implicitly solvated environments (PCM)

20131119-154028.jpg

Howard Diaz posing next to his poster

The equilibrium under study and the proposed mechanism  by which it occurs, originally proposed by Frontana-Uribe et al. looks a bit like this:

equilibrium

Dihydrocinolines in equilibrium with 1-aminoindole

mechanism

Mechanistic proposal by Frontana-Uribe et al.

The energy profile, in which all transition states were calculated with the QST3 method, is presented below, calculated at various levels of theory. Also, the Internal Reaction Coordinate (IRC) connecting both states was calculated and is shown further below in the full poster.

Energy Profile

Energy Profile

From this results we believe that a new mechanistic proposal is needed since the energy barrier for the first step is quite high (~60 kcal/mol) and hence a bit unlikely to occur through that transition state. Nevertheless this is a first approach to elucidating a mechanism and the more knowledge about it the higher the control will be on this chemical transformation.

A full version of the poster is shown below for your convenience (Spanish). See you all at the next RMFQT in Morelia 2014!

Full Poster

Full Poster

Follow

Get every new post delivered to your Inbox.

Join 1,223 other followers

%d bloggers like this: