## Orbital Contributions to Excited States

This is a guest post by our very own Gustavo “Gus” Mondragón whose work centers around the study of excited states chemistry of photosynthetic pigments.

When you’re calculating excited states (no matter the method you’re using, TD-DFT, CI-S(D), EOM-CCS(D)) the analysis of the orbital contributions to electronic transitions poses a challenge. In this post, I’m gonna guide you through the CI-singles excited states calculation and the analysis of the electronic transitions.

I’ll use adenine molecule for this post. After doing the corresponding geometry optimization by the method of your choice, you can do the excited states calculation. For this, I’ll use two methods: CI-Singles and TD-DFT.

The route section for the CI-Singles calculation looks as follows:

%nprocshared=8
%mem=1Gb

adenine excited states with CI-Singles method

0 1
--blank line--

I use the same geometry from the optimization step, and I request only for 10 singlet excited states. The CPCP implicit solvation model (solvent=water) is requested. If you want to do TD-DFT, the route section should look as follows:

%nprocshared=8
%mem=1Gb

#p FUNCTIONAL/6-31G(d,p) TD(NStates=10,singlets) geom=check guess=read scrf=(cpcm,solvent=water)

adenine excited states with CI-Singles method

0 1
--blank line--

Where FUNCTIONAL is the DFT exchange-correlation functional of your choice. Here I strictly not recommend using B3LYP, but CAM-B3LYP is a noble choice to start.

Both calculations give to us the excited states information: excitation energy, oscillator strength (as f value), excitation wavelength and multiplicity:

Excitation energies and oscillator strengths:

Excited State   1:      Singlet-A      6.3258 eV  196.00 nm  f=0.4830  <S**2>=0.000
11 -> 39        -0.00130
11 -> 42        -0.00129
11 -> 43         0.00104
11 -> 44        -0.00256
11 -> 48         0.00129
11 -> 49         0.00307
11 -> 52        -0.00181
11 -> 53         0.00100
11 -> 57        -0.00167
11 -> 59         0.00152
11 -> 65         0.00177

The data below corresponds to all the electron transitions involved in this excited state. I have to cut all the electron transitions because there are a lot of them for all excited states. If you have done excited states calculations before, you realize that the HOMO-LUMO transition is always an important one, but not the only one to be considered. Here is when we calculate the Natural Transition Orbitals (NTO), by these orbitals we can analyze the electron transitions.

For the example, I’ll show you first the HOMO-LUMO transition in the first excited state of adenine. It appears in the long list as follows:

35 -> 36         0.65024

The 0.65024 value corresponds to the transition amplitude, but it doesn’t mean anything for excited state analysis. We must calculate the NTOs of an excited state from a new Gaussian input file, requesting from the checkpoint file we used to calculate excited states. The file looks as follows:

%nproc=8
%mem=1Gb

#p SP geom=allcheck guess=(read,only) density=(Check,Transition=1) pop=(minimal,NTO,SaveNTO)

I want to say some important things right here for this last file. See that no level of theory is needed, all the calculation data is requested from the checkpoint file “adenine.chk”, and saved into the new checkpoint file “adNTO1.chk”, we must use the previous calculated density and specify the transition of interest, it means the excited state we want to analyze. As we don’t need to specify charge, multiplicity or even the comment line, this file finishes really fast.

After doing this last calculation, we use the new checkpoint file “adNTO1.chk” and we format it:

If we open this formatted checkpoint file with GaussView, chemcraft or the visualizer you want, we will see something interesting by watching he MOs diagram, as follows:

We can realize that frontier orbitals shows the same value of 0.88135, which means the real transition contribution to the first excited state. As these orbitals are contributing the most, we can plot them by using the cubegen routine:

This last command line is for plotting the equivalent as the HOMO orbital. If we want to plot he LUMO, just change the “homo” keyword for “lumo”, it doesn’t matter if it is written with capital letters or not.

You must realize that the Natural Transition Orbitals are quite different from Molecular Orbitals. For visual comparisson, I’ve printed also the molecular orbitals, given from the optimization and from excited states calculations, without calculating NTOs:

These are the molecular frontier orbitals, plotted with Chimera with 0.02 as the isovalue for both phase spaces:

The frontier NTOs look qualitatively the same, but that’s not necessarily always the case:

If we analyze these NTOs on a hole-electron model, the HOMO refers to the hole space and the LUMO refers to the electron space.

Maybe both orbitals look the same, but both frontier orbitals are quite different between them, and these last orbitals are the ones implied on first excited state of adenine. The electron transition will be reported as follows:

If I can do a graphic summary for this topic, it will be the next one:

NTOs analysis is useful no matter if you calculate excited states by using CIS(D), EOM-CCS(D), TD-DFT, CASSCF, or any of the excited states method of your election. These NTOs are useful for population analysis in excited states, but these calculations require another software, MultiWFN is an open-source code that allows you to do this analysis, and another one is called TheoDORE, which we’ll cover in a later post.

## NIST CCCBDB – Vibrational Scaling Factors & ThermoChem Data

The Computational Chemistry Comparison and Benchmark DataBase (CCCBDB) from the National Institute of Standards and Technology (NIST) collects experimental and calculated thermochemistry—related values for 1968 common molecules, constituting a vast source of benchmarks for various kinds of calculations.

In particular, scaling factors for vibrational frequencies are very useful when calculating vibrational spectra. These scaling factors are arranged by levels of theory ranging from HF to MP2, DFT, and multireference methods. These scaling factors are obtained by least squares regression between experimental and calculated frequencies for a set of molecules at a given level of theory.

Aside from vibrational spectroscopy, a large number of structural and energetic properties can be found and estimated for small molecules. A quick formation enthalpy can be calculated from experimental data and then compared to the reported theoretical values at a large number of levels of theory. Moments of inertia, enthalpies, entropies, charges, frontier orbital gaps, and even some odd values or even calculations gone awry are pointed out for you to know if you’re dealing with a particularly problematic system. The CCCB Database includes tutorials and input/output files for performing these kinds of calculations around thermochemistry, making it also a valuable learning resource.

Every computational chemist should be aware of this site, particularly when collaborating with experimentalists or when carrying calculations trying to replicate experimental data. The vastness of the site calls for a long dive to explore their possibilities and capabilities for more accurate calculations.

## On Putting up with Rejections Letters

Having a paper rejected is one of the certainties of academic life. While there are some strategies to decrease the probability of facing a rejection, today I want to focus on my tips to deal with them—particularly for the benefit of younger scientists.

There are two broad kinds of rejections: Desk Rejections and Rejections from reviewers. In any case, the best advice is never to take action after receiving the dreaded rejection letter. Take a day or two, then react accordingly with a cooler head. Remember, this isn’t about you it’s hard not to make it personal but trust me it isn’t.

The first kind, desk rejections, are provided directly from the chief or associated editors of the journal to which you submitted your work. They tend to be quick and rather uninformative except for maybe the incompatibility—to put it nicely—of your work with the scope of the journal. These are also sometimes the hardest to face since they make you feel your work is simply not good enough to be published; but they’re also the quickest and in the publish-or-perish scheme of things, time is key. After getting a desk rejection, if no other input is given, just try again; one tip—though not infallible—to chose a proper journal is to look at which journals are you citing in your own work and chose one with the highest frequency. Sometimes, editors might offer a transfer to another journal from the same publishing house; my advice is always say yes to transfers: the submission is made for you by the editorial staff, it sort of becomes recommended between the involved editors, and expedites the start-again process. Of course, a transfer does not mean you’re manuscript will get accepted but whenever offered there is a good chance the first editor thinks your work should be kept inside their editorial instead of risking you going to another publishing house. Appealing to a desk rejection is highly discouraged since it practically never works. Sure, you may think the editor will kick himself in the rear once you get the Nobel prize but telling them so, particularly in a colorful language, will not make them change their minds.

The peer review system is far from perfect, but if done properly it is still the best thing we’ve got. Some other alternatives are being tested nowadays to reduce biases like open reviews signed and published by reviewers themselves; double and even triple blind peer review (in the latter not even the editor knows the identities of authors or reviewers) but until proven useful we have to largely cope and adapt to single blind peer review (just play nice, people). In some instances the dreaded third reviewer appears, and even a fourth and a fifth. Since there are no written laws and I’m not aware of any journal specifying the number of referees to be involved in the handling of a manuscript there may be varied opinions among reviewers, so different as from ranging from accept to reject. This may be due to the editor thinking one or more of the reviewers didn’t do their job properly (in either direction) and then brings another one to sort of break the tie or outweigh the opinion of a clearly biased reviewer. If you think there are bias, consult with the editor if a new set of reviewers may be included to complete the process, more often than not they will say no but if you raise a good point they might feel compelled to do so.

Science is a process that starts at the library and ends at the library

Dr. Jesús Gracia-Mora, School of Chemistry UNAM ca. the nineteen nineties

These are truths we must learn from a young age. Any science project does not end at the lab but at the library, therefore I let my students—even the undergrads—do the submission process of their manuscripts along with me, and involve them in the peer review process (sometimes and to some limited extent even when I’m the reviewer) just so they now that getting a rejection letter is part of the process and should never be equated with the relative quality or self-worth of a scientist since that is hardly what the publication process looks at.

So, in a nutshell, if you got a rejection letter, get back on the proverbial saddle and try again. And again. And once again.

## The #LatinxChem Twitter Poster Contest

For the past few weeks, some chemists of the worldwide Latinx community have been cooking an online project devoted to showcase the important contributions to chemistry made by workers, students, and researchers from Latinamerican origin.

The result is the #LatinXChem Twitter Poster Contest which will take place 7th September during a 24 hour span and the corresponding Twitter account @latinxchem (go follow it now! I’ll wait right here.)

All chemists from Latinx origin are called to participate by registering their posters in our website latinxchem.org before August 25th. Upon registration, each poster should be classified into one of the eleven categories available and use the corresponding hashtag during the event (e.g. #LatinxchemTheo for the readers of this blog), in which prominent Latinx chemist will serve as reviewers and cast their votes for the best one in each category. Some prizes will be available, thanks to our kind sponsors (RSC, Chemical Science, ACS, Carbomex, The Brazilian Chemical Society, and more to come), but just for those registered works; if anyone wishes to present a poster without being registered at the website they can do so but eligibility for prizes remain for those who complete the register. Official languages for the poster are Spanish, Portuguese, and English.

Each category is organized by young prominent Latinx chemists; for the particular case of Computational Chemistry –the recurring theme of this blog– Prof. Fernanda Duarte (@fjduarteg) from Chile now working at Oxford University in the UK and yours truly (@joaquinbarroso) will be in charge of the #LatinXChemTheo section. Please check the website to learn about the other sections and the wonderful people working hard in the organizing committee (see below for the full list of the organizers and their Twitter handles).

The main goal of the event is to celebrate and showcase the espectacular research, education, and innovation brought to chemistry by a large and vibrant community dispersed throughout the globe of Latinx identification. We want to celebrate diversity by showcasing our contributions in the context of a global science interconnected with people from other groups.

So please visit our website, help us spread the word and get those posters ready, we’re eager to read, comment, Tweet and Retweet your work and show the world the drive and passion of Latinxs for chemistry, knowledge, and the betterment of the world through science.

Go follow us all and of course @LatinXchem too!

Gabriel Merino Cinvestav Mérida, México @theochemmerida
Miguel A. Méndez-Rojas UDLAP, México @nanoprofe
Joaquín Barroso UNAM, México @joaquinbarroso
Javier Vela Iowa State University, USA @vela_group
Diego Solís-Ibarra UNAM, México @piketin
Braulio Rodríguez-Molina UNAM, México @MolinaGroup
Paula X. García-Reynaldos Science Communicator, México @paux_gr
Liliana Quintanar Cinvestav Zacatenco, México @lilquintanar
María Gallardo-Williams North Carolina State University, USA @Teachforaliving
Fernanda Duarte University of Oxford, UK @fjduarteg
Yadira Vega Tec de Monterrey, México @yivega
Gabriel Gomes University of Toronto, Canadá @gpassosgomes
Luciana Oliveira UNICAMP, Brasil @LuBruGonzaga
Cesar A. Urbina-Blanco Ghent University, Belgium @cesapo
Ariane Nunes HITS, Germany @anunesalves
Walter Waldman Brazil, @waldmanlab

## I Want You!

This call applies to Mexican and foreign students alike but foreigners should be aware that the scholarships/stipends are awarded –mostly– by the National Council for Science and Technology (CONACYT), not by PIs themselves. Also, should be aware that the admission test will take place until February 2021 (but the registration should take place during the upcoming dates) at your corresponding Mexican Embassy or Consulate.

## Natural Transition Orbitals (NTOs) Gaussian

The canonical molecular orbital depiction of an electronic transition is often a messy business in terms of a ‘chemical‘ interpretation of ‘which electrons‘ go from ‘which occupied orbitals‘ to ‘which virtual orbitals‘.

Natural Transition Orbitals provide a more intuitive picture of the orbitals, whether mixed or not, involved in any hole-particle excitation. This transformation is particularly useful when working with the excited states of molecules with extensively delocalized chromophores or multiple chromophoric sites. The elegance of the NTO method relies on its simplicity: separate unitary transformations are performed on the occupied and on the virtual set of orbitals in order to get a localized picture of the transition density matrix.

[1] R. L. Martin, J. Chem. Phys., 2003, DOI:10.1063/1.1558471.

In Gaussian09:
After running a TD-DFT calculation with the keyword TD(Nstates=n) (where n = number of states to be requested) we need to take that result and launch a new calculation for the NTOs but lets take it one step at a time. As an example here’s phenylalanine which was already optimized to a minimum at the B3LYP/6-31G(d,p) level of theory. If we take that geometry and launch a new calculation with the TD(Nstates=40) in the route section we obtain the UV-Vis spectra and the output looks like this (only the first three states are shown):

Excitation energies and oscillator strengths:

Excited State 1: Singlet-A 5.3875 eV 230.13 nm f=0.0015 <S**2>=0.000
42 -> 46 0.17123
42 -> 47 0.12277
43 -> 46 -0.40383
44 -> 45 0.50838
44 -> 47 0.11008
This state for optimization and/or second-order correction.
Total Energy, E(TD-HF/TD-KS) = -554.614073682
Copying the excited state density for this state as the 1-particle RhoCI density.

Excited State 2: Singlet-A 5.5137 eV 224.86 nm f=0.0138 <S**2>=0.000
41 -> 45 -0.20800
41 -> 47 0.24015
42 -> 45 0.32656
42 -> 46 0.10906
42 -> 47 -0.24401
43 -> 45 0.20598
43 -> 47 -0.14839
44 -> 45 -0.15344
44 -> 47 0.34182

Excited State 3: Singlet-A 5.9254 eV 209.24 nm f=0.0042 <S**2>=0.000
41 -> 45 0.11844
41 -> 47 -0.12539
42 -> 45 -0.10401
42 -> 47 0.16068
43 -> 45 -0.27532
43 -> 46 -0.11640
43 -> 47 0.16780
44 -> 45 -0.18555
44 -> 46 -0.29184
44 -> 47 0.43124

The oscillator strength is listed on each Excited State as “f” and it is a measure of the probability of that excitation to occur. If we look at the third one for this phenylalanine we see f=0.0042, a very low probability, but aside from that the following list shows what orbital transitions compose that excitation and with what energy, so the first line indicates a transition from orbital 41 (HOMO-3) to orbital 45 (LUMO); there are 10 such transitions composing that excitation, visualizing them all with canonical orbitals is not an intuitive picture, so lets try the NTO approach, we’re going to take excitation #10 for phenylalanine as an example just because it has a higher oscillation strength:

%chk=Excited State 10: Singlet-A 7.1048 eV 174.51 nm f=0.3651 <S**2>=0.000
41 -> 45 0.35347
41 -> 47 0.34685
42 -> 45 0.10215
42 -> 46 0.17248
42 -> 47 0.13523
43 -> 45 -0.26596
43 -> 47 -0.22995
44 -> 46 0.23277

Each set of NTOs for each transition must be calculated separately. First, copy you filename.chk file from the TD-DFT result to a new one and name it after the Nth state of interest as shown below (state 10 in this case). NOTE: In the route section, replace N with the number of the excitation of interest according to the results in filename.log. Run separately for each transition your interested in:

#chk=state10.chk

#p B3LYP/6-31G(d,p) Geom=AllCheck Guess=(Read,Only) Density=(Check,Transition=N) Pop=(Minimal,NTO,SaveNTO)

0 1
--blank line--

By requesting SaveNTO, the canonical orbitals in the state10.chk file are replaced with the NTOs for the 10th excitation, this makes it easier to plot since most visualizers just plot whatever set of orbitals they read in the chk file but if they find the canonical MOs then one would need to do some re-processing of them. This is much more straightforward.

Now we format our chk files into fchk with the formchk utility:

formchk -3 filename.chk filename.fchk
formchk -3 state10.chk state10.fchk

If we open filename.fchk (the file where the original TD-DFT calculation is located) with GaussView we can plot all orbitals involved in excited state number ten, those would be seven orbitals from 41 (HOMO-3) to 47 (LUMO+2) as shown in figure 1.

If we now open state10.fchk we see that the numbers at the side of the orbitals are not their energy but their occupation number particular to this state of interest, so we only need to plot those with highest occupations, in our example those are orbitals 44 and 45 (HOMO and LUMO) which have occupations = 0.81186; you may include 43 and 46 (HOMO-1 and LUMO+1, respectively) for a much more complete description (occupations = 0.18223) but we’re still dealing with 4 orbitals instead of 7.

The NTO transition 44 -> 45 is far easier to conceptualize than all the 10 combinations given in the canonical basis from the direct TD-DFT calculation. TD-DFT provides us with the correct transitions, NTOs just paint us a picture more readily available to the chemist mindset.

NOTE: for G09 revC and above, the %OldChk option is available, I haven’t personally tried it but using it to specify where the excitations are located and then write the NTOs of interest into a new chk file in the following way, thus eliminating the need of copying the original chk file for each state:

%OldChk=filename.chk
%chk=stateN.chk

NTOs are based on the Natural Hybrid orbitals vision by Löwdin and others, and it is said to be so straightforward that it has been re-discovered from time to time. Be that as it may, the NTO visualization provides a much clearer vision of the excitations occurring during a TD calculation.

Thanks for reading, stay home and stay safe during these harsh days everyone. Please share, rate and comment this and other posts.

## DFT Estimation of pKb Values – New Paper in JCIM

As a continuation of our previous work on estimating pKa values from DFT calculations for carboxylic acids, we now present the complementary pKb values for amino groups by the same method, and the coupling of both methodologies for predicting the isoelectric point -pI- values of amino acids as a proof of concept.

Analogously to our work on pKa, we now used the Minimum Surface Electrostatic Potentia, VS,min, as a descriptor of the availability of Nitrogen’s lone pair and correlated it with the experimental basicity of a large number of amines, separated into three groups: primary, secondary and tertiary amines.

Interestingly, the correlation coefficient between experimental and calculated pKb values decreases in the following order: primary (R2 = 0.9519) > secondary (R2 = 0.9112) > tertiary (R2 = 0.8172). This could be due to steric effects, the change in s-character of the lone pair or just plain old selection bias. Nevertheless, there is a good correlation between both values and the resulting equations can predict the pKb value of an amino group within less of a unit, which is very good for a statistical method that does not require the calculation of a full thermodynamic cycle.

We then took thirteen amino acids (those without titratable side chains) and calculated simultaneously VS,min and VS,max for the amino and the carboxyl group (this latter with the use of equation 2 from our previous work published in Molecules MDPI) and the arithmetical average of both gave us their corresponding pI values with an agreement of less than one unit.

This work is now available at the Journal of Chemical Information and Modeling (DOI: 10.1021/acs.jcim.9b01173); as always a shoutout is due to the people working on it: Leonardo “Leo” Lugo, Gustavo “Gus” Mondragón and leading the charge Dr. Jacinto Sandoval-Lira.

## A New Gradúate Student. Raúl Márquez

We’re always happy at the lab when a student defends their dissertation thesis and now it was the turn of Raúl Márquez-Avilés to do so with flying colors.

The title of his dissertation is “Molecular Dynamics Simulations of 5 potential entry inhibitors for HIV-1“. He performed 500 ns long molecular dynamics simulations of the CD4 – gp 120 proteins interacting with one or several molecules of various lead compounds with inhibitory properties. The leads were obtained previously in our group (by Durbis Castillo, now at McGill) from a massive docking library of ca. 16 million compounds, all having a central piperazine core (Fig1)

The protein gp120 is a surface glyco-protein located at the surface of the HIV virus which couples to the CD4 protein on lymphocytes-T, being this the first step in the infection process of a healthy cell; generating inhibitors of this coupling could help stop the infection from spreading systemically. Four systems were devised: (SB) The reference state for which only gp-120 and CD4 were considered, (S2) A single ligand molecule was placed in the Phe43 cavity of gp120 to assess their inhibitory capacity, (S3) the ligand was placed right outside the Phe43 cavity to assess their entry capacity, and (S4) five ligand molecules were placed outside the Phe43 cavity of gp120 to force their entry (Fig2). Their binding energies were calculated using MM-PBSA and although all five ligands show statistically similar results as inhibitors all five exhibit a stronger binding energy than the reference proving their efficacy in preventing the coupling of the virus to the healthy cell. As a bonus, his research on system S4 shed light on the existence of an allosteric site on gp120 that will warrant further research in our group.

This work is still pending publication.

Raúl Márquez has always proven to be a hard working person who is also very self-sufficient student, a very cheerful labmate, and, as I just learned yesterday, an avid chess player. I’m sure he has a bright future in whichever endeavor he chooses now. Congratulations Raúl Márquez-Avilés!

## Locked out of your Linux Session

Funny enough I was unable to log into my Linux (Ubuntu) session and I realized this might be a more common problem that it seemed. So, if you keep getting redirected to the login screen after typing your correct password over and over (and over and over), there’s no need to panic.

This usually has to do with the .Xauthority file, so from the login page press Ctrl+Alt+F1 which will bring you to the command line where you can login with your usual credentials. Once logged in, search for the .Xauthority file and check that it is owned by you and not the root

ls -l ~/.Xauthority
-rw------ 1 root root 1 feb 11 13:13 /home/joaquin/.Xauthority

Use the following command to change ownership

in my case both group and username are joaquin. You may need to ‘sudo’ it. If that doesn’t work try deleting the file altogether, upon login it will be created again.

rm -rf ~/.Xauthority

In any case, if any of these solutions worked, press Ctrl+Alt+F7 to go back to the login screen and now you should be able to get in.

These solutions are quite straightfoward but if the problem persist you may need to update the system or downright install it again from the command line we opened at the begining.

## Good Luck to Dr. Jacinto Sandoval

We’re sad to begin this year by saying farewell to Dr. Jacinto Sandoval-Lira who held a postdoc position in our lab for two years with a DGAPA – UNAM scholarship, a very competitive and highly sought-after position here in Mexico. Dr. Sandoval will now relocate to the Technological Institute of San Martín Texmelucan in Puebla, Mexico, whose students will be fortunate to have him as a tutor and a teacher of chemistry in the environmental engineering department.

During the past two years we’ve worked together in various projects, mainly the excitonic transference between photosynthetic pigments but also in calculating reaction mechanisms and solving chemical equilibria problems with various computational approaches, but apart from the research Dr. Sandoval was also a co-organizer of the past Meeting on Physical Chemistry, organized a local course on the use of Dens Tool Kit (DTK), as well as our weekly lab seminars and taught various graduate and undergraduate courses on molecular modeling, chemoinformatics and computational chemistry, not too mention all the collaborations he has brought to our lab in the field of organic chemistry of which I regard him as an expert. He really has been a force of nature!

Aside from a brilliant scientist and a hard working one, Jacinto is an exceptional human being and a great friend. His attention to detail, his drive, and his willingness to help others reach their full potential make him an ideal colleague and an ideal professor. I don’t wish you luck, Jacinto, you don’t need it: I wish you success!