# Rigid and Relaxed Potential Energy Surface Scans (PES Scan) in Gaussian 03 and Gaussian 09

Another common question on CCL…

The use of Internal Redundant coordinates (through the Opt=ModRedundant option) must not be overlooked! This option performes a geometry optimization at each step while maintaining the scanned variable constant, which is referred to as a Relaxed Potential Energy Surface (PES) Scan. Using internal coordinates becomes compulsory and a well-defined Z-Matrix is preferable. A common mistake is to define the variable to scan without taking into consideration all the other variables which depend from the first. This means that functional groups can be destroyed on each step and in the best case scenario the optimization algorithm will put it back together again at, of course, some unnecessary computational cost. The worst case scenario would consist of the molecule’s complete distortion from any physically achievable structure, making the calculation end with an error or by providing meaningless results. Therefore the use of wildcards is compulsory and it is illustrated below.

Lets say we want to calculate the energy barrier associated with the rotation of a dihedral angle on ethane. If we define the dihedral angle 2-1-5-8 (fig.1) to be scanned from 180.0 deg to 0.0 by 45.0 degrees increments then the second step of the scan would look distorted like the molecule on figure 2. Since the value of D[2-1-5-8] is kept constant until convergence then the algorithm (provided a decent level of theory is indicated when it comes to molecules more interesting than ethane) will increase the values of dihedral angles D[2-1-5-6] and D[2-1-5-7] until reaching the convergence criteria at a regular structure.

This would result in an unefficient method for computing the desired barrier. Hence, before defining the variable to be scanned we must select carefully, through the use of wildcards all other variables associated with it.

* 1 5 *

2 1 5 8 S *#steps #increment*

blank line

The first line will make all other dihedral angles around the C2-C3 bond to behave in the same way as D[2-1-5-8], *steps* is the number of increments in the variable to be performed from the starting geometry and *increment* is the step size to be taken in the variable (Note this values are always in integer format so a decimal point is always to be used! e.g. for our hypothetical calculation we could write the second line as 2 1 5 8 S 4 45.0 with which the dihedral angle would be scanned from 0.0 degrees to 180.0 degrees in 45.0 degrees steps) In this way, the second step of our hypothetical calculation would look like the molecule on figure 3, and the convergence criteria would be reached much faster.

This same procedure holds for scanning any other kind of geometrical variable. Common mistakes and error messages during geometry scans are associated to lack of memory, maximum number of steps reached or energy divergence due to inconsistent geometries; when analyzing for an error is always best to take a look at the optimization progress in order to find out at which step the molecule started to break apart.

Rigid PES are much easier to perform although they can be less inforative in terms of a real dynamic process. In this case the “scan” option is the one we want to use. The molecule must be defined in Z-Matrix format below of which a blank line must be placed and after that the following information

*Variable as defined in Z-Matrix initial value Number of steps *and *Increment value*

So a rigid PES scan for the previous example would look like

*Z-Matrix*

D1 180.0 3 60.0

If another variable is placed with the same info below this line then the program will perform all possible combinations while printing a summary with the energy at each conformation. Sometimes this summary doesn’t have enough space for the numbers to be printed and a set of stars ****** are shown where the energy is supposed to be. No panic the energy was calculated it just cannot be printed; you either have to use a visualization program such as GaussView in order to read the energy at each conformation or just browse directly -and patiently- through the output file

The intention of this blog was not to become a Gaussian 0x support forum but I’m glad that so many people -specially grad students- have found it helpful in their research.

As usual all comments and ratings are welcome

Posted on April 19, 2010, in Computational Chemistry, Gaussian, White papers and tagged Computational Chemistry, Gaussian. Bookmark the permalink. 251 Comments.

Interesting article. However, I’m still looking for figure 2 …

Hi Moritz

Thanks for pointing the missing figures out. Lately I’ve been a bit careless with the scheduling of my posts. Anyway, the mistakes were taken care of.

Cheers!

Dear Sir

I am trying to converge my oligomer with carbon nanotube. but now i did not get any yes. now it gives the energy oscillation. RMS do no goes to minimum level. Can you help me about it.

You probably need to start from a better geometry. Try looking for pre-optimized CNT’s on some database.

I hope this helps!

Dear Joaquin,

Can you help me for CASSCF//CASPT2 type of calculations with Gaussian 09 ?

All the papers say, “CASSCF//CASPT2” and “CASSCF MP2” methods were used to study excited states. They dont explain in detail or may be I dont understand.

The main points are

i) How to choose the “Active Space” to generate the reference wave function for computing different type of excited states(pi-pi*, n-pi*, n-sigma*)

ii) How to deal with Rydberg States, means how to identify them by looking at molecular orbital coefficients and orbitals visually so that I can delete them from active space when calculating valance states ?

iii) when to use State-averaged CASSCF and what does it mean physically?

iv)How to compute oscillator strength of different transitions from a CASSCF calculations ?

Can you please give an example of any molecule (say Indole) step by step, so that I can learn it in general for any molecule.

Please refer me some review articles or book chapter to read about these Ab initio methods.

Hello Sayan

Your question is a bit complex to be answered as a comment reply, therefore you’ve inspired me to write a post on multiconfigurational methods!

I will try, however, to answer quickly to your main concerns:

Selecting the active space is a compromise between accuracy and computing resources. You could use a Complete Active Space (CASSCF) for best accuracy but then you would need to have a lot of computer resources! Probably those guys whose articles you are reading had them and thats why they went for it. Usually you can choose as many virtual orbitals as valence orbitals you have, this way you can assure that all low -energy- lying transitions will be taken into consideration.

For an introduction to the methods -specially the math behind them- check “Reviews in Computational Chemistry Vol.14 T. Daniel Crawford and Henry F. Schaefer III, An Introduction to Coupled Cluster Theory for Computational Chemists, pp. 33-136. (2000)”

For explicit and step-by-step examples with G09 check:

http://gaussian.com/g_tech/g_ur/k_casscf.htm

also try a book called “Exploring Chemistry with Electronic Structure” by James B. Foresman and AEleen Frisch. If you have any trouble finding it just let me know and maybe I can send it to you in pdf format.

I hope this helps you get started with your calculations. Stay tuned for a larger post on CASSCF and thanks for the idea!

Thanks for reading

Joaquin

Hi Joaquin,

Thanks a lot.

I an reading those articles.

It will be great if you start a post on Multiconfigurational methods.

sayan

Hello again Sayan!

Like a month later but that post on multiconfigurational methods is finally out there! Please let me know if you find it useful, also if you feel it should have something else.

Have a great day!

Very interesting blog. Big up!

I’m glad you find it interesting! Thanks for reading!

Dear Joaquin,

congratulations for the interesting and useful article! I perform relaxed scans rather often. Maybe you have the answers to two old important questions of mine (plus a third stupid one):

1) when the number of total scans to perform is small (less then 50), the CPU is active all the time, as it should be. But going above, the CPU seems to be active only cyclically: the CPU usage hystory on the task manager looks like a chainsaw. Increasing the number of scans more and more, the “space” between the teeth of the chainsaw becames larger and larger, making the “time the PC is calculating”/”time the PC is not calculating” ratio too unconvenient above let’s say 200-300 scans total. And int he periods when the CPU is not working, the HD is active. Any idea if it is possible to do anything about it?

2) I have seen (with gaussview 4) that when you scan 2 variables (say values A, B, C, D for one variable and a, b, c, d for the other), Gaussian calculates the conformations sometimes in the order Aa, Ab, Ac, Ad, Bd, Bc, Bb, Ba, Ca, Cb….. (more typically) and sometimes in the order Aa, Ab, Ac, Ad, Ba, Bb, Bc, Bd, Ca, Cb…… Any idea if it is possible to control it rather then making the software decide for you??

3) About saving a structure witn GaussView in Z-matrix format, any idea how you can save it in the format

O

H 1 r2

H 1 r3 2 a3

Variables:

r2 1.10000

r3 1.10000

a3 109.50000

rather then

O

H 1 1.10000000

H 1 1.10000000 2 109.50000024

I would really appreciate any suggestion!!!

Gius.

Hello Gius!

You have very interesting questions here (that usually means that the answers aren’t quite available).

1) I assume you mean a single scan with more than 50 steps or combinations, right? I’m affraid I don’t have a straight answer at this moment but I can look into it for you. I guess there is some issue regarding writing and reading from disk when you have a large calculation. I know this is not the optimum procedure but have you tried splitting it into smaller scans which you can later put together to get the full profile? It could also be related to the processor you are using, try running it in another computer with different processor architechture to see what happens. I will look into it, probably you can share an input for testing (my email address is joaco_barroso’at’yahoo’dot’com)

2) This has to do with standard and input orientations, the program uses whatever it suits best for the computation. There should be an IOp to control whether to always use one or the other but I don’t have it at hand. I suggest you check gaussian’s official website and browse through the IOp’s in the apropriate section.

3) This is a weird one. Right now I just did a test on GaussView 3.0 (tonight I will try on 4.0) and I get it just the way you want it! What version are you using? I’m trying to find a control for it on the preferences menu but there is nothing like that. Are you saving in *.gjf(com) file-format? This HAS to be controlable, thats for sure! I will get back at you when I get the answer.

Sorry for not sheding too much light on your questions but I will work a bit on them, I promise. Bare with me, I just came back a few days ago from a long trip :D

I will stay in touch. Cheers!

Just a quick solution to 3) mentioned in this comment – Im using GaussView 5 and can force the z-matrix to be as you want it by using the Connection Editor and selecting Opt All from the Z-matrix Tools drop down list. Hope this helps

Thanks Richard!

I will include your suggestion in another post if you don’t mind :) This way others will be able to find it more easily.

Thank you all for reading!

Stupid question, but is it normal that for quite big molecule, after performing full 360 degrees scan, I found very close but not exactly equal energy for the conformer 0 and 360?

If the molecule is very big (and by big I mean it has many degrees of freedom and not necessarily having too many atoms) then it is quite reasonable not having the same energy after the full scan (I am assuming a RELAXED scan). Some of the other variables may have shifted from their original values. You also have to consider just how close is very close when you talk about energy in atomic units since 1 a.u. is 627.5 kCal/mol

Have a nice day and thanks for reading!

Dear Joaquin,

I have a small question for you. I’m interested in performing the relaxed scan but also in calculating the thermochemistry parameters in each of the optimized points, in order to evaluate the Gibbs Free Energy profile of a reaction.

Do you know if that is possible in G09?

Thanks a lot,

Saludos Cordiales,

Juan Manuel.

In order to evaluate the thermochemistry parameters you need the Hessian matrix at each point, in other words, you would need to have the vibrational analysis of each one of those points along your scan. I don’t know if you can set it to be performed. I also don’t know if the approach is valid since you wont have minima on the PES along the scan. I would look into the G09 manual but I think you will have to extract each step manually (there is a tool within Gaussian to do so, if I can remember how to I will let you know) perform the vibrational analysis and with that get the thermochemistry associated to it.

I hope this halps a bit.

¡Saludos!

Dear Joaquin,

I did scan for torsion angle in biphenyl rings form zero to 180 degree. My question is there a keyword in Gaussian to generate the energies of molecular orbitals at each step?

Thanks,

Moro

Hi Moro

I’m not aware of any keyword that can do what you want. Have you tried using pop=full? you can always do it manually but if you are taking too many steps then it would be quite troublesome. I’m pretty sure you can’t do it with the rigid scheme but with the opt=modredundant option you should be able to do it. Let me know what is it that you have tried so I can be of more help, ok?

Have a nice day!

I think I found the option you need: Try pop=always along with other options you want (e.g. pop=(full,always) or something like that). This option calculates the Mulliken population at every step of an optimization and not just at the beginning and end. I think you should try it.

Have a nice day!

Hi,

You had mentioned about 2 1 5 8 in relaxed scan. could you please send me an input file where in * 1 5* is shown.

In z-matrix, its hard to find no. written like that. and in xyz coordinates, i don’t know how to give input for scan range.

thanks.

Hi ashutosh!

I’m not sure I understand. The first line (* 1 5 *) indicates that ALL dihedral angles will behave the same (whether they are 2 1 5 8 or 4 1 5 6) this is done so ALL three H atoms in each C atom rotate at the same time or eventually you would put H2 right on top of H4 and the calculation would crash (not to mention it would be meaningless)

The procedure is the same both for cartesian coordinates and z-mat formats.

Here is an example of the syntax

“Link0 chk and stuff”

“charge multiplicity”

“molecule specification, i.e. coordinates”

blank line

* 1 5 *

2 1 5 8 S 3 20.0

blank line

This means: Make all dihedral along the 1-5 bond behave the same way. Then take dihedral 2 1 5 8 and SCAN it (that is what the S stands for); do three additional steps (besides the starting point) increasing the dihedral 20.0 degrees in each step (include .0 or it will crash). Here you have scanned the rotation of the C-C bond a total of 60 degrees.

I hope this helps but if it doesn’t get in touch again and I’ll be glad to help.

Have a nice day!

Hi,

I do have small doubt.I want to perform geometry optimization with fixed bond length two atoms from different molecules. how can i do that? My input file in cartesian coordinates form. how can i change it to Z-matrix format.

please help me

Thanks in advance,

satesh

Hi Satesh!

You have to define the distance between those two atoms and then use the normal procedure for a scan (as described in this post). You don’t need to have it in Z-matrix format but if you insist, you can change it with GaussView:

Open the file and then save it again but before you hit SAVE notice there is a checkbox in that window that says “write cartesians” uncheck it if it’s checked and that will save your input file in z-matrix format. A similar procedure works with many more visualization programs; I think even OpenBabel can do this (Ask it to convert your gaussian-cartesian file into a gaussian-z-matrix one)

I hope this was helpful.

Have a nice day and please excuse me for my delay.

Dear Joaquin,

I am trying to perform a relaxed scan with “Opt=ModRed” option, where few atoms are fixed in my system. The syntax used to start scan of bond length is “B atom1 atom2 bond length S #steps stepsize”, which is placed right after the Z-matrix is defined and follows other frozen atom details.

Each time I submit the job, it does specify in the output ”

Number of optimizations in scan= 21

Number of steps in this run= 732 maximum allowed number of steps= 732.”

but after finding the optimization geometry coordinates the job fails with the following error:

” Iteration 1 RMS(Cart)= 0.00294869 RMS(Int)= 0.00319231

New curvilinear step failed, DQL= 4.79D-03 SP=-2.72D-02.

RedCar failed.

Error termination via Lnk1e in /opt/gaussian/gaussian09/g09/l103.exe”

I have tried different options and systems, but all are failing and I am begining to think about the right way of defining the system.

My complete input looks like this:

%chk=SiC13-O_scan1.chk

%mem=22GB

%nproc=8

#p Opt=modred b3lyp/6-31G nosymm

SiC 13X4 run with Oxygen

0 1

H

H 1 r2

H 1 r3 2 a3

H 3 r4 1 a4 2 d4

H 1 r5 2 a5 3 d5

H 5 r6 1 a6 2 d6

–

–

–

Si 68 r119 62 a119 61 d119

Si 70 r120 8 a120 5 d120

Si 69 r121 7 a121 6 d121

O 111 r122 88 a122 85 d122

r2= 1.7384

r3= 3.2335

a3= 151.42

r4= 1.7469

a4= 24.48

d4= 16.82

r5= 12.4080

–

–

–

r120= 1.8639

a120= 107.34

d120= 162.43

r121= 1.8818

a121= 106.03

d121= 117.49

r122= 1.7084

a122= 111.29

d122= 50.92

X 1 F

X 2 F

X 3 F

X 4 F

X 5 F

X 6 F

X 7 F

–

–

–

X 56 F

X 57 F

X 58 F

X 59 F

X 60 F

X 61 F

X 62 F

B 122 89 3.7 S 20 -0.1

I have put – – – to cut down on the input file. I would appreciate if you can suggest me the way around it Or pointout any mistake in the input.

Thanks for your time.

Regards,

Satyender Goel

Hi Satyender!

Your input seems fine! well, it seems a bit odd to request for 22GB in memory, which it wont take (depending on the architecture I believe the most you can take is 2GB). Anyway, you have a convergence failure, the RMS is NOT below the threshold. I assume you have tried this optimization with the Quadratic Convergence option (SCF=QC), haven’t you?

Are you reaching the first optimized state? I mean have you seen the famous 4 YES’s?

Another possibility could be the definition of your system. If you have a surface (which seems to be the case) and then you are trying to get something close to it by decreasing a bond length then you may be deforming the surface by not properly freezing every atom in it. Could this be the case?

What you should definitely look into is more ways of controlling the convergence criteria or methods. I hope this helps but if it doesn’t then send me more info about your output and about those options you’ve tried already.

Have a nice day!

PS Sorry for the delay in my response

Dear joaquinbarroso,

I have the same problem during the first optimized state, means i have seen the 4 YES’s…. could you please suggest me the solution for this kind of error ?? i would appreciate your help.

Neeraj Kumar

Hi Neeraj

If you also have a convergence problem then I suggest you use a different convergence criteria in the SCF. Try SCF=QC in your route section.

I hope this helps but if it doesn’t then send me more info about your system.

Have a nice day

Thank you for this article, it was very informative, I hadn’t used the wildcards like that previously.

I also constantly get the “RedCar failed” error however and I can’t figure out what to do. In my system I am dehydrating a water molecule from a carbon backbone. The first geometry converges (4 YESs), but I get the RedCar error almost immediately after. My .log file shows 4 YESs, then Optimized Parameters in an ascii box with their values, then straight to the iterations, curvilinear step failure, and RedCar. My system gets slowed drastically by scf=qc, do you have any other suggestions? Could the issue be that my dehydration coordinate is moving the oxygen and leaving the hydrogens behind? Is there anyway to rigidly move the hydrogens with the oxygen in between steps and then release them to be optimized?

Thank you so much for your time,

Tim Courtney

Hi Tim

Probably you are right about the hydrogens being left behind, you could check if that is the case with the help of any visualization tool such as GaussView. A way to prevent it is with the use of wildcards: Let’s say your water molecule has the following numbering scheme H1-O2-H3 and the rest of the molecule is the C backbone you talk about; what you want to do is type the following

* 2

2 *

* 2 *

and then the scan parameters (I imagine you are stretching the C# – O2 bond). The first two lines indicate all bonds to O2 in the connectivity list (make sure it is present with geom=connectivity, also make sure to generate it with GaussView by clicking on show connectivity at the calculation panel), the third line defines all bond angles with O2 in the middle.

This is only from the top of my head, I HOPE it works but as usual feel free to ask again if it doesn’t.

Now, maybe what you want is to look into an IRC (Internal Reaction Coordinate) calculation in order to assess the bond breaking process of a water molecule from an organic fragment. This will tell you about the energetics of the process. Hmmm I have an idea for another post :)

Have a nice day, let me know if my suggestion helped.

Dear sir,

I am new to this post session.Few days ago only i saw this post.i have read some of your answer.Very useful to clear the doubts.Now only i have started to the PES scan.But i am not getting results.My molecule is 2-amino-5-nitropyrimidine.I dont know which dihedral angle sholud be scaned.can i do scan of two dihedral angles at atime.Kindly help me .

Thanks a lot in advance.

Hello Meenakshi

You say you don’t know which dihedral angle to scan, well the answer is you can scan whatever dihedral angle you want! although in your molecule you only have two bonds that can be rotated (which is what you want to do with a dihedral angle scan) the H3N-C bond or the O2N-C bond. In any case your defined angle should look like

H(#) N(#) C(#) N(#)

where the second N atom is one of the N atoms in the heterocyclic ring. Go with:

* N(#) C(#) *

H(#) N(#) C(#) N(#) S nsteps size

and try it again

Remember “size” is an integer in degrees so asking for 60 degrees wont work, you need to specify 60.0 degrees, ok?

I hope this helps but if it doesn’t send me your input file and I will take a look into it.

Have a nice day!

Dear Joaquin,

Could you please talk on two coordinate scan ( as two angle at same time) How one can do that ? How to view result in gaussviev?

Thanks

Vishal

Dear Vishal,

Performing a two coordinate scan is possible in the same way as you perform a single scan, just add the second variable to be scanned below the first, e.g., if you are scanning an angle you get to write something like:

1 2 3 S 1.0 10

which will scan the 1-2-3 angle by 10 1.0 degrees increment. A second angle would be:

1 2 3 S 1.0 10

4 5 6 S 1.0 10

You will see in your gaussian output a matrix covering all combinations between those two coordinate changes.

I hope this helps!

Dear sir,

I am new to this post session.Few days ago only i saw this post.i have read some of your answer.Very useful to clear the doubts.Now only i have started to do the PES scan.But i am not getting results.My molecule is 2-amino-5-nitropyrimidine.I dont know which dihedral angle sholud be scaned.can i do scan of two dihedral angles at atime.Kindly help me .

Thanks a lot in advance.

Hello Joaquin Barroso,

I am trying to locate minima occurring when two molecules containing a certain functional group are close to each other. I performed a rigid PES at a given distance, scanning valence angles and dihedral angles, and got nice surfaces showing some minima. This was performed with structures optimized at MP2/6-311++G** which were transformed to a X-matrix and I used the Scan keyword putting the steps and lengths behind the variable definitions of the two parameters.

If I try to repeat this procedure using Opt=ModRedundant and

* #A #A *

#A #A #A #A S x y

#A #A #A *

#A #A #A S x y

(where #A are numbered atoms and x,y are the steps and steplengths), Gaussian lets the molecules relax to distances and angles that are all over the place, and not comparable at all.

Can I still trust the that my results with the rigid PES, which indicate the interaction is enthalpically favoured, are somewhat energetically accurate?

(given the basis set of course – I could do a scan with ccVTZ)

Thanks for your help!

Hello Kaupang!

Wow! Very interesting question you have here! I would like to understand more about it but I would have to take a look at your files on the relaxed scans. What exactly do you mean by “all over the place”? are the molecules breaking/bonding/distorting/pushing-each-other-away?

I am positive your results with the rigid PES are to be TRUSTED, unless a chemical reaction (i.e. bond breaking/forming) could occur. Do you know if that is the case? About the energetic accuracy, try using a benchmark to prove that the approach is valid: take a small bimolecular system of which you know the interaction enthalpy and run the RIGID scan for it, see how well it compares with the experiment.

THe employed level of theory for the optimization seems quite robust (although I know nothing about the molecular identity of your system); what level of theory are you employing for the rigid scan? I would try using the same one. Using ccVTZ will yield a more accurate energy but if your molecule contains only light atoms then I don’t think you will gain much more accuracy than you already have achieved.

Anyway I would correct the interaction energy by taking the basis set superposition error into consideration (COUNTERPOISE keyword in Gaussian)

Thanks for reading and have a nice day!

Hi again!

Thanks a lot for your answers, I now put some more faith in my energies!

I was looking at the ModRedundant keyword, and it seems I need to revise my input file. I think it may be missing some constraints, because the molecules move away from each other (downhill in energy) along a third parameter that I hadn’t specified, while the angles and distances I have specified are adhered to.

Would locking say a dihedral angle, while scanning a distance and and angle, cause the calculated energies to be less valid (i.e. too high)?

As for the rigid PES, it was performed at MP2/6-311++G** and the formula of each molecule is C3H3N3O, so I think you’re right about the accuracy gain using ccVTZ.

Using a small reference system for which there are thermodynamic data available is a great idea! I will look into this! Also I will have a look at the counterpoise keyword.

Thanks again for a great blog!

Hi Kaupang!

Thanks for your kind words! I’m very glad to know you find this blog useful.

As for your question, You can’t know a priori whether the energy will be higher or lower given a particular constrain, BUT you don’t have to worry about them being less “accurate”, you are only taking some restrictions on a complex hypersurface, in other words you are not taking every variable into consideration but whether that variable is important or not depends from one system to another. I guess what I’m trying to say is that you can report your findings always stating that this constraints you imposed (maybe for technical reasons) only represent the results of a PART of the complex potential energy hypersurface.

I think you are on the right path, congratulations and keep up the good work.

Have a nice day!

Hello sir

My molecule name is 2-amino-5-nitropyrimidine.I would like to perform 3D PES SCAN .sir help me do and explain me that how to visualize that using gaussview.thanks in advance.

Hello

Sorry for the delay on my response. If you want to visualize the hypersurface generated by the scan I am affraid you cannot do it with GaussView, you need other visualization software or generate enough points during your PES scan so you may get a nice graph/plot with a data analysis program such as Origin, GNUplot or in the worst case scenario with MS Excel.

Good luck!

Dear Dr. Joaquin,

I am experimental physicist and trying to explain one of my result on the basis of PES using Gaussian 03. In brief, I have found multiple peak in Kinetic energy distribution of H3+ ion, ejected from singly ionized CH3SOCH3 (DMSO). H3+ ejection is due to association. Now, I want to calculate transition states and PES. Please guide me how I should proceed? There are many excited states are responsible for formation of H3+!..

Regards, Rajesh

Hi Rajesh

The PES is calculated the way described above; that is you have to select a variable (or many) and ask for a scan or an optimization in the opt=modredundant way.

I am no expert in Transition State theory but I suggest you to search for the QST1 and QST2 keywords in the Gaussian website.

I hope this helps. Thanks for reading!

hi dr.

how i can scan calculation simple

HI,

I want to know how to calculate free energy change (delta G) from the thermochemistry data in G09. Can i calculate for an isolated molecule or it should be a reaction.

thanks

Hi Marawan

Wow, your question has a rather long answer. I think it is time to write a new ‘white paper’ on the topic :) thanks for the inspiration :D

In the mean time take a look at this:

thermochemistry G09

it is very helpful and comprehensive.

Thanks for reading and I hope you keep on visiting!

Dear Joaquin,

I am performing a rigid scan calculation with G03. I am changing 3 variables (torsions). By using gaussview, can I select the visualization of couples of scans? In other words, lets name A, B, and C the three torsions that I wish to vary. Can I visualize three separate graphics where the energy is function of scan(A)-scan(B), scan(A)-scan(C), and scan(B)-scan(C)? From the “partial” log file I can only visualize a 2D graphic with Energy as function of “unorganized” torsions.

Thank you,

luciano

Hi Luciano!

I’m afraid the short answer to your question is NO. Even the latest GaussView5 can only generate surfaces when you scan two and only two variables. You will have to use Origin, gnuplot, MSexcel or other program to plot your surfaces.

What do you mean by partial? I assume your calculation ended correctly. As for “unorganized”, you should get three columns (at the “Summary of the potential surface scan” section) with your torsions, where the first is being scanned and the remaining two retain their original values, then the second moves to the next value and so on. Are you seeing something different?

I hope this helps. Thanks for reading!

Dear Joaquin,

Thank you for the answer. By “partial” I mean that the calculation is not finished yet. Overall I am performing more than 2000 scan, so that I am routinely checking the calculation.

I thought that with GV I could choose which scan visualize, but I was wrong.

I will then need to wait the end of it to use the results and paste them into origin, excell etc.

Thank you again,

luciano

Dear Joaquin,

I am new in this blog. I am a Ph.D student in physical organic chemistry. First of all, I should confess that I don’t know much about how to perform fluorescence calculation of a compound in gaussian 03. One of my friend have asked me to perform fluorescence and UV calculation of a ruthenium based complex. I have already calculated TD calculation of UV which is very straight-forward in g03 program. But, I am not able to calculate fluorescence for him. Actually I am confused about that how can I get first excited singlet structure which is necessary (I think) from where I can relax or…..I am really confused. I have only gaussian 03 program and gauss view 3. Is it possible in g03. If possible please help me. Please provide me an example of a very small molecule. I am waiting for your kind reply.

With best regards,

Thanking you,

Debabrata.

HI Debabrata!

I have to confess I’m no expert in excited states calculations either. Try looking for Dr. Tomas Rocha Rinza, he is also mexican and an expert on the topic. Besides, doing calculations with transition Metals is difficult enough, let alone calculating their excited states where you have to take into consideration the number of available terms depending on the wavelength of the excitation.

It seems to me you have to enroll yourself in a Theoretical Chemistry class, your confusions will then disappear. I mean this with the best intention possible

Have a nice day!

HI Dr,

I’m trying to find the global minimal structure of Cytosine using G09, i found only one imaginary frequency in the enolic form but none in the Keto form, what does this mean..?, second, could you supply me with sample inputs to perform QM/MM in Gaussian.

Thanks

It means you reached a TS with the enolic form; this is a good thing! remember there is a phenomenon called keto-enol tautomerism, this means the keto form is more stable than the enolic one, hence, the last one is a Transition State, which you were able to find! if you try to distort the molecule in the direction of the imaginary frequency you will find it is trying to become the keto form.

As for the second part of your question I think it depends on the kind of calculation you are trying to achieve, so let me know a bit more about it so i can help you properly, ok?

Have a nice day!

Dear Dr,

Thank you very much for your extremely useful message. Then, in QM/MM context, i want to do some Docking and MD simulation then refine with the QM/MM methods to extract an accurate free energy of binding.

Thanks

Marawan

Hello again Marawan

I think in Gaussian the only way to do this is by means of an ONIOM calculation. Take a look at this link for a brief introduction to the method and how to work it:

http://gaussian.com/g_tech/g_ur/k_oniom.htm

Best regards!

Hello Dr. Barroso,

I’m trying to do a rigid energy scan to gauge intermolecular forces. As it turns out, this is my first time working with Z-matrices (I’ve been letting GaussView take care of my inputs and have been using cartesian coordinates exclusively). I’m trying to scan the bond distance between two atoms that are not formally bound and I can’t figure out what variable in my Z-matrix to scan. Do you have any suggestions short of re-writing the whole z-matrix?

Thanks,

Scott

Hi Scott!

If you are using GaussView then it is a very simple matter: Whenever you save your files (i.e. molecule specification along maybe with the calculation settings) there is a check box that says something like “write cartesians” which is usually checked by default; just uncheck it and it will save your input file as a Z-matrix. Then you just need to find which atoms you want to scan, say x and y, you will probably see something like:

x y B1 z A1 w D1

This means: x is bonded to y with a bond length of B1 Angstoms, and they both form an angle of A1 degrees with z, the three of which form a plane which has a dihedral angle of D1 degrees respect to the plane defined by y,z,w. (Sorry if I’m being too descriptive) Then you find B1 in the variables section and scan it with the procedure described in the post

B1 S #steps incrementsize

Hope this helps but if it doesn’t I’m always here!

Have a nice day

dear sir

can i do a pes scan using CBS-QB3 method?i did a 3-D scan(by changing 2 bond lengths)but the scan was not visualized from gaussview3(also the scaned points were not in the output file).sir,Do i able to do this kind of a scan using Gaussian98?

Hi Hasara

You must need a lot of accuracy on your scan! I think you could use it, but I’ve never tried it. I guess the reason why you can’t see it is because of the versions: G98 is too old and Gview3 also is, so probably visualization of 3D hypersurfaces is not available. The problem is that Gaussian no longer gives any support to these old versions. I suggest you take your question to the ccl.net forum since I don’t have access to either versions as to make a test.

In short, my suggestion is: Try to get an update :S Sorry.

On the other hand I think it is very interesting what you try to do with a complete basis set method! Run the calculations (even if you have to run them separately, I mean, not in a single scan file but doing the scan manually) and then use some other program to visualize the surface (origin, gnuplot, even MSExcel)

I hope this helps! Thanks for reading!

Dear Dr Barroso,

I’ve been trying to performe a PES Scan for the adsorption of Xe on Benzene, but i cant define the correct input. I need that Xe gets closer to the hollow site on benzene. Could you help me?

Thanks a lot!

Hello K!

This is a bit of a tricky one, if you try to define the scan in terms of the internal coordinates of the benzene + Xe atom only. If you are using gaussview, use it to place a Dummy (X) atom (not a ghost (Bq) one! there is a difference) in the center of the benzene ring, then place the Xe atom at any given length of that X atom making sure the angle with any other C atom is 90.0 degrees, as well as remaining dihedrals, i.e., Xe-X-C = 90.0 Xe-X-C-C = 90.0

Scanning the Xe-X length then becomes straightforward.

I hope this helps

Have a nice day!

Joaquin,

Excellent site that is accesible and respectful to newcomers! I appreciate that. When I was new to all of this I encountered a number of really simple issues that were “too simple” for anyone on CCL to take seriously.

My question is this: What is your take on the relation between a relaxed potential energy surface scan and calculated zero-point energies for the (global) minimum structure? For instance, I am scanning the C-Cl bond length on a medium sized system and I am interested in the potential well for S0, S1, and T1 Surfaces. The C-Cl is bound in the S0 state, dissociated in the S1 state and barely bound in the T1 State. I fully expected the S0 and T1 States to be bound, and the excited state to be marginally bound or unbound (as this unbound intermediate is thought to be seen in experimental photochemical studies).

So, this all seems well and good. However, I got to be concerned about the depth of the potential well for the T1 state, and I have gone back and calculated it with many small steps to more closely follow the C-Cl surface. I end up with a small barrier to dissociation of about 7 kcal/mol. So, I revisited the S0 potential energy surface to evaluate the ZPE versus the derived C-Cl barrier. The fully optimized geometry is calculated to have a ~90 kcal/mol zero-point energy, but the optimized scan shows a barrier of only ~75 kcal/mol.

So, what gives? I know the C-Cl is bound in the ground state, and not so loosely that any vibrational energy at RT would cause dissociation. Yet, it seems that I have found an inconsistency here. I suspect that I am missing something trivial and comparing apples and oranges. Then the question is – what constitutes a bound state when scanning a bond cleavage?

Regards,

Soren

Dear Barroso,

i’m trying to perform a scan over the FCOH molecule, so that i fix the value for the dihedral angle and leave free the other parameters to get an optimization…theproblem is that i don’t know what do i have exactly to write down in the input line.

Would you please help me?

Hi Maite

If I understand correctly you want to FREEZE a dihedral angle and then optimize the rest, right? Well, this is not scan, in fact it would be the opposite to a scan. What you need to do is to write Opt=ReadFreeze in the route section, then in the end of the file, after a blank line following the molecular coordinates write:

noatoms atoms=i,j,k,..

where i j k are the atoms to be frozen.

I hope this helps!

Dear Joaquin,

Thank you very much for your interesting and informative homepage. I have also a question concerning the ModRedundant scans. Is it possible to scan a coordinate, optimize the structure and perform another calculation where I calculate the energy on a different level of theory? The result would be two PES on two different levels of theory? I want to do a scan with a cheap method and switch after that to an expensive method to get better “energies”, without optimization.

Is this possible?

Thank you very much in advance,

patrick

Hi Patrick!

Sorry for the delay in my response and thank you very much for your kind words.

I’m not sure if it is possible to get it done in one go but I’m sure it is not very informative, it might not even be valid! you see, the potential energy surface changes with the level of theory employed (there is one real PES but in order to fully describe it you would need an infinitely large basis set). This means, among other things, that a minimum on a PES at a certain level fo theory might not be so at another. You might try to achieve this two PES scans by hand but you should then somehow validate the process, which might be done by calculating frequencies at both levels of theory on each point and showing they are both minimum, although this would be more time consuming than just generating the PES scan at the expensive level of theory.

I hope I understand your question correctly and if I did then I hope this was helpful.

Cheers!

Dear All;

Thank you very much for sharing such a vital informations. even though i am new to Gaussian; i have been learning and practicing the package so that i can effectively apply it in my research work. however, i found it difficult to prepare an input file for a one electron reduction reactions. could any one help me on this issue? i mean is it only changing the spin and charge? or is there any other way? just for your information; i have tried to do just that and in some cases the job ends with error.

Hi Ermi

It would be very helpful if you gave us more information about your problem so we can help you solve it.

Take care!

Good morning sir,

I am Ph.D student in Physics Department

I want to do TDDFT caluculaton(UV)and PES scan, So how can i run Gaussian 03 programme through z matrix.

Thanking you Sir

Hi Venkata

Have you actually tried to do this? I know that TDDFT is only available for energies and gradients so I don’t think it is possible. Try running the calculation and if it fails then I think you should do it by hand, calculating the UV spectra at each point previously defined by yourself.

Thanks for reading!

Hi Dr.Jacquin,

I want to perform a relaxed scan and i saved the coordinates file as a Z matrix with Molden, but, the dihedral angle of interest was not in the obtained Z matrix. I viewed it using the Z mat editor in molden and captured its numbering in the dihedral sequence and manually added it in the Z matrix coordinates file but, i got the error

Symbol “dih3” not found in Z-matrix.

so, what i have to do, the same happened when i tried with saving the coordinates as mol2 file using gaussview and tried to convert to the Z-matrix using openBabel. Would you plz tell me what to do.?

Thanks

Marawan

Hello Marawan!

With Gaussian you don’t need to find the specific variable and that is a relief! if you are trying to scan the dihedral angle defined by atoms 1 2 3 4 type the following at the end of your input file:

* 2 3 *

1 2 3 4 S n deg.deg

Don’t forget to use the option Opt=ModRedundant in your route section. In this same post there are some instructions to do what you are asking about; give it a try.

Take care and thanks for reading!

Dear Dr Barroso,

I am trying to adsorb a molecule over a substrate using Gaussian 09. In this regard I want to calculate ENERGY of the system of substrate and adsorbed molecule as a function of distance. I am aware of SCAN utility in Gaussian. However, I am not able to use SCAN for my calculations. Though I am using Z-matrix, it is difficult to move the whole adsorb molecule or set variables for it. I hereby request you to kindly help me in this regard.

Thanks and regards

Deepak

Hi Deepak!

This can be a bit tricky, specially if the adsorbate is made of many atoms. Try using the ONIOM scheme applying a high layer to the adsorbate and a medium layer to the substrate. That way you can move whatever is defined as a single layer without distorting it.

Another possibility, if you don’t want to try ONIOM, is taking very small steps and letting the adsorbate optimize at each of them so if you only scan the distance between the surface and say the central atom, the rest of the molecule will follow.

I hope this helps

Thanks for reading!

Sir,

Good Evening.

Yes, your advice has helped me a lot in making final decision. I will try out things suggested by you.

Thanks and warm regards

Deepak

dear sir,

im uma.i want to use the external potential to read the gaussian software. is it possiple?and one more thing how can i use the keyword pseudo potential.

looking forward for your reply. thans in advance.

with regards,

uma

Hi Uma!

I don’t understand what do you want to do :( You want Gaussian to read an external potential? what kind of potential? Electrostatic? Now, when you ask about “pseudo potential” this means to use an effective core in which the inner electrons of a given atom(s) are substituted with a parametrized potential. To use an external pseudopotential you need the instruction pseudo=read and then you can just copy it from https://bse.pnl.gov/bse/portal

I hope this helps Uma but feel free to ask again anything you need, ok?

Take care and thanks for reading!

dear sir,

thank you very much for your help. i want to use the pseudopotential in my job.that is my query. now im trying for that. again thanks a lot.

Great post; very helpful.

I’ve gotten this to work for RHF and AM1, but gaussian doesn’t seem to scan if I select UFF — it just runs an energy calculation. Can gaussian perform dihedral scans using molecular mechanics (particularly for a method in which i don’t have to specify atom types or parameters)?

Thanks,

chris

Hi Chris!

Thanks for your comment. I’m not aware about the capabilities of G09 to run scans with MM methods. Now, have you tried the opt=ModRedundant option? It might work. Hope this helps (even a little)

Thanks for reading! Greetings from Mexico!

I have opt=modredundant and dihedral scan work with RHF…

When I switch to Amber=… It says that no variable in z-matrix so cartesian optimization will be done.

I tried to bring the redundant definition before the FF definition but gaussian complains about that.

any comment on this?

Thank in advance

Regards

I found this! :D

http://www.ccl.net/chemistry/resources/messages/2007/09/27.002-dir/

Cheers

Dear Prof.Dr.Jaquine,

Thanks for the help the job is running now. one more question, if i want to restart, shall i write:

opt=ModRedundant=restart or opt=restart=ModRedundant

and what about the BERNY Algorithm provoked by the opt=Z-matrix option, do you think it’s more efficient or the opt=ModRedundant is better.

Full of respect

Marawan

Hello Marawan! I’m glad to know my advice worked!

The correct way to restart your calculation would be typing opt=(ModRedundant,Restart)

Just be careful of providing the same checkpoint filename in the link0, that is:

%chk=filename.chk

where filename is the name of your file. Keep it at the location where it was originally generated or it wont find it. You should also include ALL the keywords of your original route section (the one that begins with the symbol ‘#’)

I hope this helps!

Ok, so right now I am having a little trouble with restarting my scan. I didn’t think to check the .chk file name. But it seems that HPC i submit it to appends the full path to the name of my file. So when I place opt=restart in the route section I only get the optimization of the last structure as my output file.

Maybe this is obvious for some people but I’ve just started using a HPC and it’s killing me until I saw this post and put 2+2 together to get 4.

And further reading you suggest to add:

example1.chk

%nosave

So I don’t accidentally overwrite the last good file from my scan and you recommend to always work with a copy too. All good points and I hope to pass them along to my labmates.

Thanks for the great resource!

Best Regards,

Agapito

Dear Dr.Jacquine,

Thanks for your reply. another question, what does it mean if after opt freq output normally with no imaginary frequencies an the 4 yes , at the end of the file i found as if the job is not yet converged (Max displacement NO) just before the end of the output..??!!

Thanks

Hello again Marawan!

I have never seen such error!!! I wouldn’t worry if the rest of your output looks fine but let me take a look into it so we know for sure this is not a problem but just a glitch.

Have a nice day! Thanks for reading!

Dear Dr.Jacquine,

I am a beginner in Gaussian calculation. Your blog is quite helpful for me. I have tried the dihedral scan for 4-nitroimidazole molecule, using “opt=modredundant” and it gives me successful results. But I have two questions:

1. Is it possible to use “opt=z-matrix” to do the same dihedral scan. Of course I mean that all the connected dihedrals will act the same way with the defined scan dihedral, just as Fig3 in your blog.

2. I try the rigid scan. the input file is as below:

#p scan mp2/6-31g(d) nosymm

scan dihedral

0 1

N

C 1 B1

N 2 B2 1 A1

C 3 B3 2 A2 1 D1

C 1 B4 2 A3 3 D2

H 1 B5 2 A4 3 D3

H 2 B6 1 A5 5 D4

H 5 B7 1 A6 2 D5

N 4 B8 3 A7 2 D6

O 9 B9 4 A8 3 D7

O 9 B10 4 A9 3 D8

B1 1.36750604

B2 1.32573787

B3 1.35954392

B4 1.37008763

B5 1.00835493

B6 1.07676154

B7 1.07391005

B8 1.44435416

B9 1.24534003

B10 1.23906260

A1 111.84868839

A2 103.61776449

A3 108.06000051

A4 126.24622045

A5 122.49515194

A6 124.18323846

A7 122.03317257

A8 115.77832107

A9 118.27564307

D1 -0.02965322

D2 0.03161998

D3 -179.98523495

D4 -179.99122908

D5 179.98804274

D6 -179.98502187

D7 -179.98527891

D8 0.0 6 30.000000

But the result shows that only the dihedral D8 change in the scan job and the connected dihedral D7 doesn’t change at all. How can I get the similar dihedral scan as using “opt=modredundant”. I want D7 dihedral to behave as the same way as D8 in the rigid scan. I will appreciate for your help. Thanks!

Hi Julian!

Question 1 is simple to answer: Yes you can do it with the opt=modredundant option as discussed in the post.

About question 2: If you scan a variable only that variable will change; if you want two variables to be scanned then you have to ask for it, BUT if you do then you will get ALL possible combinations, most of which wont make sense (actually, only the diagonal elements of the combinatorial matrix will make sense). I would suggest to perform the scan manually and not with the option scan.

I hope this helps and thanks for reading

Have a nice day!

Hello

I will try the dihedral scan for a kind of proline molecule, using “opt=modredundant” . I will use and imaginary dihedral angle C(a-1)-O(-1)-Cd-Ca (19 2 4 7)

How I should consider which other coordinates have to be relaxed?. I know that using standard dihedral angles such as 1 2 3 4, the distance 2-3 have to be consider * 2 3*R.

But in this case seems more complicated . I will really appreciate your advice in this issue.

Thanks

Here is my structure. I was using z-matrix input, but in order to save space here I saved it in Cartesian coord.

# b3lyp/6-311g(d,p) opt=modredundant

proline cis

0 1

C 1.80714800 -0.98263300 0.05629900

O 2.83895700 -1.06984500 -0.58928800

N 1.01564000 0.13907900 -0.02673600

C 1.42383600 1.25017500 -0.91017800

H 1.46420900 0.92183000 -1.95089500

H 2.43209300 1.57860000 -0.63975100

C -0.09976400 0.50377400 0.85490000

H 0.06771100 0.15242400 1.87439800

C -0.09154600 2.04577300 0.79247800

H -1.06723700 2.45967200 1.04373500

H 0.64053500 2.42986100 1.50829500

C 0.36738300 2.33212900 -0.64592300

H -0.47321900 2.22291000 -1.33698300

H 0.76998400 3.33819300 -0.77170100

C -1.48821500 -0.04145300 0.45333700

O -2.44941100 0.17413200 1.17321900

N -1.54588900 -0.73272300 -0.71413800

H -0.68797600 -0.84544100 -1.23140000

C 1.35804700 -2.07396200 1.01412200

H 1.96588100 -2.95601800 0.82295800

H 1.51858900 -1.75945600 2.05001500

H 0.30069300 -2.32221200 0.90201100

C -2.78121000 -1.29796800 -1.23117800

H -3.04740100 -0.85057700 -2.19302700

H -2.69502300 -2.38070600 -1.35679900

H -3.56753300 -1.08288000 -0.50932300

Hola Fernanda!

I’m sorry for the delay in my response.

I just built your molecule from your cartesian coordinates and the atoms 19 2 4 7 seem like a very odd selection for scanning. If you scan this dihedral then you will break the C-N bond within the first steps. Since atoms 2 and 4 (Oxygen and Carbon, respectively) are not bound, this selection doesn’t make too much chemical sense. Please send me more information about what you are trying to achieve with this calculation so I can help you, ok?

Now, if you insist on trying this scan then be warned about the lost of the N-C bond.

I would think you are trying to scan the N-C bond, sorry but I dont understand what you are trying to do. Please send me more info.

Have a nice day!

Dear Dr Joaquin,

About 1 month ago a CCL user asked a question on the listserv (http://www.ccl.net/cgi-bin/ccl/message-new?2011+06+11+001 ) but this question has never revived a informative response (as much as I can see on the web). Since I have exactly the same problem, I copy and paste that question here and I hope you can help me. Thank you. Somaye

Dear CCL users:

I was trying to run transition state search with qst3. I wanted to freeze a few atoms in the geometry. However, it seems that Gaussian03 cannot read in the third modredundant input section for the initial guess of transition state. I have also modified a Gaussian test job and see the same error message. It always complains that ” WANTED AN INTEGER AS INPUT. FOUND A STRING AS INPUT.

4 F

?

”

Here are the test job and the error message in the output file. I have also tried qst2 with only the reactant and the product sections which works well. I have tried to remove the third modredundant section but it then asked for “AN INTEGER AS INPUT”.

Could you please help me to fix this?

Test job input:

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

%nproc=2

%mem=1GB

%chk=test4131.chk

#p uhf/3-21g opt=(modRedundant,qst3)

Gaussian Test Job 413:

Path optimization, Reactant

0,2

C

O 1 R1

H 1 R2 2 A2

H 1 R3 2 A3 3 D3

H 1 R3 2 A3 3 -D3

R1=1.48

R2=1.085

r3=1.08

A2=110.0

A3=110.0

D3=120.0

4 F

5 F

Gaussian Test Job 413:

Path optimization, Product

0,2

C

O 1 R1

H 1 R2 2 A2

H 1 R3 2 A3 3 D3

H 1 R3 2 A3 3 -D3

R1=1.48

R2=1.90

r3=1.08

A2=30.0

A3=110.0

D3=120.0

4 F

5 F

Gaussian Test Job 413:

Path optimization, TS

0,2

C

O 1 R1

H 1 R2 2 A2

H 1 R3 2 A3 3 D3

H 1 R3 2 A3 3 -D3

R1=1.48

R2=1.40

r3=1.08

A2=70.0

A3=110.0

D3=120.0

4 F

5 F

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Test job output with the error message: the first two sections are read in very well…

————————————————–

Gaussian Test Job 413: Path optimization, Reactant

————————————————–

Symbolic Z-matrix:

Charge = 0 Multiplicity = 2

C

O 1 R1

H 1 R2 2 A2

H 1 R3 2 A3 3 D3 0

H 1 R3 2 A3 3 -D3 0

Variables:

R1 1.48

R2 1.085

r3 1.08

A2 110.

A3 110.

D3 120.

The following ModRedundant input section has been read:

X 4 F

X 5 F

————————————————-

Gaussian Test Job 413: Path optimization, Product

————————————————-

Symbolic Z-matrix:

Charge = 0 Multiplicity = 2

C

O 1 R1

H 1 R2 2 A2

H 1 R3 2 A3 3 D3 0

H 1 R3 2 A3 3 -D3 0

Variables:

R1 1.48

R2 1.9

r3 1.08

A2 30.

A3 110.

D3 120.

The following ModRedundant input section has been read:

X 4 F

X 5 F

——————————————–

Gaussian Test Job 413: Path optimization, TS

——————————————–

Symbolic Z-matrix:

Charge = 0 Multiplicity = 2

C

O 1 R1

H 1 R2 2 A2

H 1 R3 2 A3 3 D3 0

H 1 R3 2 A3 3 -D3 0

Variables:

R1 1.48

R2 1.4

r3 1.08

A2 70.

A3 110.

D3 120.

WANTED AN INTEGER AS INPUT.

FOUND A STRING AS INPUT.

4 F

?

Error termination via Lnk1e in /share/apps/gaussian/g03/l101.exe at Sat Jun 11 03:03:56 2011.

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Here is some information from Gaussian 03 manual which is not so helpful…

This method is requested with the QST2 and QST3 options. QST2 requires two molecule specifications, for the reactants and products, as its input, while QST3 requires three molecule specifications: the reactants, the products, and an initial structure for the transition state, in that order. The order of the atoms must be identical within all molecule specifications.

The QST3 option allows you to specify a better initial structure for the transition state. It requires the two title and molecule specification sections for the reactants and products as for QST2 and also additional, third title and molecule specification sections for the initial transition state geometry (along with the usual blank line separators), as well as three corresponding modifications to the redundant internal coordinates if the ModRedundant option is specified. The program will then locate the transition structure connecting the reactants and products closest to the specified initial geometry.

ModRedundant

Add, delete or modify redundant internal coordinate definitions (including scan and constraint information). This option requires a separate input section following the geometry specification. When used in conjunction with QST2 or QST3, a ModRedundant input section must follow each geometry specification.

Thank you very much for your time and help!

Hello Somaye!

Wow! I have naver taken so long to reply to a comment, I’m awfully sorry about that. Even more sorry I’m about the fact that I don’t know exactly what the problem might be sine the syntax of the input is correct. Try the following:

a) Try not freezing atoms 4 and 5 in the proposed TS (third part of your input file). At first I thought there could be a geometrical incompatibility but the Error Message doesn’t really correspond to this issue.

b) Freeze them but use the Opt=ReadFreeze option instead of the ModRedundat after molecule specification leave a blank line then type atoms=4,5 and then leave another blank line. You can do this for only the suggested TS structure or for the three molecule specs in your file.

I’m sorry I don’t have any more information available at the moment but if I know of something I will let you know, ok? Also please let us know through this blog about another solution you found to circumvent this problem.

I hope this helps

Have a nice day! and once again I’m sorry for such a long delay

Thanks Dr Joaquin for your reply. Actually the solution to this problem is that the freezed atoms (modredundant part) have not be written for TS in QST3, just you need to indicate those freezed atoms for the reactant and product at the end of their coordinates. Obviously, the freezed atoms for reactant and product need to be the same. Somaye

dear sir,

this is uma. i want to perform a QM/MM calculation for protein molecules. i searched a lot, but couldnt understand properly. so please help me regarding this.

thank you in advance

Hello again Uma!

Take a look at the following links. I think what you need is to perform an ONIOM calculation. It’s so exciting to work with such large molecules! I hope it helps but if it doesn’t let me know and I will help you from closer

http://www.gaussian.com/g_tech/g_ur/k_oniom.htm

http://www.gaussian.com/g_tech/g_ur/k_mm.htm

http://www.gaussian.com/g_tech/gv5ref/tut0.htm#O1

http://www.gaussian.com/g_tech/gv5ref/tut0.htm#O2

Thanks for reading!

dear sir,

again uma. thanks for the link you provide. i tried with gauss view. i got it. thanks a lot

Dear Dr.Jacquine,

Thanks for the reply, and plz don’t leave us for long time again. Another question, I’m doing optimization at the B3LYP/6-311++G** for a nuclesoides and i got 1 imaginary frequency at 90cm-1, when i redo a single point calculation at the same level of theory i got the same , but, once i use B3LYP/TZVP i don’t get any imaginary frequency, is it possible, and if it’s possible, doesn’t mean that i’m on a global minimal structure or not. I want to say that i have done a PES scan at B3LYP/6-31g then i did the optimization from the deepest minima.

Thanks

Marawan

Dear Marawan

Sorry for the long absence, hopefully this wont happen again.

A potential energy surface is constructed by an infinite number of single point calculations (at least in theory) each point would then be dependent of the level of theory used so it is not weird for a point to be a minimum at one level while the same geometry corresponds to a TS at a second level. If you are finding the global -or near global- minimum with a scan at a level1 and then you reoptimize the structure with level2 then you are ok!

Sometimes its a good idea to help the calculation by performing an optimization at a low level of theory and then improve on it by performing a second opt with a higher level; thus disregarding the imaginary frequencies previously found.

I hope I got your question right. Thanks for being a frequent reader! I appreciate it.

Joaquin

Dear Dr.Jacquin,

I’m planning to run QM/MM on a drug design project and i would like to use ADMP, do you agree or you think that CPMD or CP2k is a better option to choose. CPMD is extremely time consuming and difficult to setup and run. But i think ADMP is more user friendly. Does it have an accuracy concern, does it have any other advantage-or disadvantages- over CPMD..?

Marawan

Hi Marawan,

I’m not quite familiar with any of them but from what I hear you are taking the right choice by going with ADMP.

I wish you success with your project, I personally find drug design a very exciting area of work

Thanks for stopping by

Hi sir,

Can you give me some ideas on spin-orbit coupling leading to the mixing of singlet and triplet states. Also can you mention the best program to calculate SOC. I couldn’t find any input example for the SOC calculation using Gaussian.If you have already written about this please share me the link.

Thanking

Jo

Dear Dr.Joacquin,

I’m asking if Gaussian can use external FF for QM/MM like Charmm or i have to stuck with the available ones, also, using electronic embedding in QM/MM does require any manual intervention after setting up the input using gaussview bcz i have read that i have to treat the MM atom near the model system in a way that to reduce the charges of the MM atoms in order not to cause overpolarization of the near by QM region.?

Thanks

Hi Marawan

If you want to specify a new force field as an external file, you have to use the following input

# UFF=SoftOnly

TITLECARDRead-in force field example

molecule specification

@$g09root/g09/file.prm {specify the location of the FF file)

where file.prm is the external FF you want to load

About the border this is a very tricky issue. I suggest trying a couple of options in a trial and error fashion.

I hope this helps!

Dear Dr.Joacquin,

How can simulate the vibrational spectra at different temperature, i used the Temperature=x key word but it seems that i got the same spectra, do i need to specify something else in the input.

Thanks

Hi Marawan!

You can’t calculate the vibrations at different temperature. The vibration analysis is performed by analysis of the Hessian Matrix which includes all the net forces on atoms, which in turn are calculated from the Fock Matrix (in the HF formalism or from the density matrix in the DFT formalism). All this calculations are performed on a single molecule in vacuum; the concept of temperature lacks any meaning in these conditions.

The Temperature=x keyword will only use a different temperature for the Thermochemical analysis performed after the vibrations are calculated. The software uses these vibrations as part of the particion function (statistichal mechanics) to calculate the thermodynamics of a homogeneous molecular ensemble.

Sorry about that. I hope this helps you in some way. Please excuse the lateness of my response but I’ve been very busy as you can possibly imagine.

Have a nice day!

Dear sir,

we are the utilizer of this blog, Can we restart the incomplete any gaussian job due to power failure.Sir we read that any incomplete job can be restarted from checkpoint file but we don’t know how to restart the work from checkpoint file.

Another query is how to submit a structure of dimer of 6-anilinopurine.

Hello Meenakshi

If you ran a calculation examplefile-1.com or examplefile-1.gjf and the output is examplefile.log then you surely have a file called examplefile.chk modify the input file with the following instructions:

Create a new file called examplefile-2.com and edit it the following way

%chk=examplefile-1.chk

%NoSave

# HF/6-31… (all the keywords from the original file) opt=restart (if you were doing an optimization) or guess=check (if this was only a single point calculation) geom=check (this will recover the geometry from the checkpoint file)

No further specifications are needed.

The %NoSave line is very important! this will make the program read the checkpoint file but it wont re-write it! if you forget this line you will loose the first checkpoint file and the program will stop very quickly. It is always a safer option to work with a copy of the checkpoint file.

Hope this helps!

Dear sir ,

Thank you so much sir. I have got clear idea about the restart process.Already i am doing my PES SCAN job by your help only.Regarding that can I copy the SCAN FILE AS A MOVIE FILE TO SUMMIT in International journal.

I’m very glad I could help!

I’ve never heard of a movie file sent as supplementary information to a journal! Anyway, if you have GaussView 4.x or 5.x you can save a movie of the optimization by opening your log file (click on the “Read Intermediate Geometries” checkbox) and then go File -> Save Movie.

If you are using another software such as Molekel then you have to go to Menu -> Render -> Animation -> Save this will generate many (and I do mean many) image files that you can then put together as a movie with Microsoft Movie Maker or something.

Have a nice day and thanks for clicking!

Dear Dr.Joacquin,

I have a doubt about the method you mentioned to calculate the exact barrier (dihedral barrier using scan run in Gaussian-09). If we use wildcard as suggested by you, we will no longer be scanning only one dihedral angle, but a collective dihedrals based on the specific structure at hand. So how reliable would the calculated barrier be by this (wildcard) method and also if we want to report this barrier in a publication would you still give it as a specific dihedral angle or mention that it is a collective barrier?

Thank you.

Chidambar Kulkarni

Hi Chidambar!

I think this collective method, as you call it, is the ONLY reliable one since its the one that makes sense chemically, i.e. in the ethane example, rotation of the C-C bond involves the motion of the entire CH3 group and that is achieved by the use of wildcards. If you omit the wildcards then one of the CH bonds will start getting closer to the other on the sam C atom and the energy will rise until diverging by collaps of both H atoms when the calculation will surely crash.

I don’t think its necessary to clarify in publications the use of a collective set of internal coordinates or not since, as I already said, this is what makes sense chemically. Only if you were to scan a dihedral in which you only want a single atom to move then you would have to specify it.

I hope this helps, but as usual if I can help with something else, just let me know.

Have a nice day!

Dear Dr.Joacquin,

Thanks a lot for the above suggestions. I want to know; how to model (classical mainly using MD) of organic molecules. I am asking this question because most force fields exits for biomolecules and very few for organic molecules in general. Can you suggest how do I go about the simulation of organic molecules. I have tried AMBER’S GAFF, but it has not been that helpful. Always we end up doing more parametrisation.

Thanks

Chidambar

Dear Chidambar

I’m not familiar with force fields for small organic molecules, I’m sorry. It would seem to me that being so small, so many changes would occur in them from one kind of system to the next, therefore the external/tedious parametrization would be almost mandatory, right?

Sorry I can’t help you with this one, still if I find out something I will let you know.

Thanks for reading. Best wishes!

Hello Joaquin,

I really enjoyed reading your blog and finding out new ways to use Gaussian. I have been trying to use the opt=z-matrix Geom=Modify option in order to check the potential energy surface of a molecule. A problem I encountered is that the bond I want to twist looks fine when I use the “Modify Dihedral” option, but when I save it, more then one dihedral is changed, so there is no use in turning just one dihedral. If possible, I would be happy to know how can I re-build the molecule so that it’s Z matrix will be “diagonal” with respect to the dihedrals I am interested in.

Hi Amir,

I’m sorry but I need more information about your molecular system to fully understand your problem and help you.

I wait for your info.

Hi Joaquin, I think the best way for me to explain my goal is with this short video:

I would be happy to send you one of the .gjf files if you have the time :)

Amir

Dear Dr.Jacquin,

Is there’s a way to extract the 2D PES relaxed scan data from G09, Gaussview can only enable the extraction of 1D PES but not the 2D PES….??????

Thanks

Marawan

Hello Marawan!

If gaussview can’t display it then you need to use another program like origin to plot the data. Copy the summary into origin (or gnuplot if you know how to use it) and generate a surface (x and y axis relate to the variables you were scanning while the z axis is the potential energy)

Gaussview has its shortcomings I know.

Have a nice day!

Many thanks! so useful!

Dear Dr.Jacquin,

I am a PhD candidate of mechanical engineering, and I am not so familiar

with computational chemistry. But I need to calculate barrier height

(in eV) for passing a gas molecule through a carbon structure. So I will

appreciate if you let me know if it is possible to calculate with

Guassian and how? I am so sorry if my question is elementary!!

Hi Yali,

It is possible to calculate such a barrier with Gaussian if you are talking about a nanotube or something similar. If you are talking about a mesoscopic or macroscopic wall then you have to reduce it somehow so it can be actually calculated. I would suggest a hybrid approach such as ONIOM for calculating the barrier with the methodology described in this post.

Send me an email with more details so I can help you properly, ok?

Best wishes!

Thank you very much for your answer.

I sent the email to jbarraso@unam.mx

Is it right?

Thank you again.

Sorry, Yali. You have only one little mistake in that email: It is barroso with two ‘o’s

I’m waiting for your files

Thank you very much for letting me know my mistake in the email address.

I sent again an email to jbarroso@unam.mx

I hope this time you recieve my email!

Thank you again.

Can you provide a brief example of how to do this using Gamess?

Sorry, I have no experience working with Gamess.

Best wishes!

Hi sir,

Can you give me some ideas on spin-orbit coupling leading to the mixing of singlet and triplet states. Also can you mention the best program to calculate SOC. I couldn’t find any input example for the SOC calculation using Gaussian. Can you get me a sample Gaussian input for the SOC calculation? If you have already written about this please share me the link.

thanking

Jo

hello sir

i am doing phd.how can calculated pes scan in gausview 5.0 give me idea with example

hi,

This is thirupathi, I am trying to do NBO analysis for the conformations that are possible around the scanning a dihedral angle(360 degrees with 72 steps ). But i dont want to optimize the structure. I used the z matrix formulation as u mentioned in the articles. But How can I make the all the dependent angles and distances changes as dihedral angles scan?

one more question is

is it possible to do a opt=modredundant with pm3 and the NBO with b3lyp for all the structures in a single calculation?

Hi Thirupathi,

You have to make the whole group move around as you scan a dihedral. This is described in the post. Lets say you want to scan the dihedral defined by 1-2-3-4, so the bond that is twisting is the 2-3 bond. You have to make all other atoms bonded to atom 2 to move accordingly, you accomplish that by typing the following lines:

* 2 3 *

1 2 3 4 “scan variables”

this means: Make all dihedrals around the 2-3 bond behave the same way. Check the post carefuly and you find the answer.

Thanks for reading!

how to give input in gaussview

Hi Joaquin

i want to optimize a FeC structure with scan, but it doesnt work.

i´m working with z-matrix. What i want to do? I´ve been trying to use the opt=Modredundant

Help me, please

sir

electron density diagram for my sample not get red region,potential diagram not clear. why ?

I don’t understand what you mean but if you are trying to plot the electron density of a molecule you shouldn’t get a red region because you don’t have two phases. This only happens when you plot orbitals, there you have a positive and a negative phase which are usually depicted with two different colors.

I hope this helps you finding an answer.

Thanks for reading and have a nice day!

hi

Im a new user of Gaussian03

And I want to know what is the difference between chk file and the log file?

thanks

A chk file is a binary file in which a lot of the information about the calculation is stored in a way that may be used for further calculations. A log file is a text file which contains all the results of the calculation in a (pseudo) human readable fashion.

Most of the time you can get the results from just the log file but if you require some post-processing of the data (say, plotting molecular orbitals) you might need to use the chk file or its text transformed version, the fchk file.

Please go to http://www.gaussian.com/g_tech/1.htm for further information.

Have a nice day!

hi,

thanks for your explanations

Have a nice day 2

Dear Dr.Joacquin,

I am devi research schlolar in computational chemistry.I am running the guassian 03 for NMR calculation,giao method with B3LYP with 6-311++G (d,p) basis set with DMSO solvent(IEFPCM) of my compound 3-(2-oxo-2-phenyl-ethylidene)-1,3-dihydro-indol-2-one, I got the error message like this

requested convergence on Max density matrix = 1.00D-06

requested convergence on energy = 1.00D-06

No special action if energy rises.

Out of memory error in routine PCMQM-DMIVC[1 End =5971753 Max core=5431964)

use %mem=7MW to provide minimum energy to complete this step.

I request you to help me to solve this problem.

Thanks in advance,

with regards,

Devi.

Dear Joaquim

I want to optimize an FeC structure with scan, but it doesn´t work.

I need freeze some atoms in the z-matrix, scan the dislocation of carbide and optimize each point.

I am running the gaussian 09 for scan calculation. I´ve been trying to use the opt=loose in Z-matrix

I´ve been reading this forum to finding an example, but until today i didn´t find.

Have you got windows messenger? I´ve some doubts about optimization with PES.

Help me please

———————-

See the error

———————–

CompJT: unrecognized IType=10.

Error termination via Lnk1e in C:\G09W\l1.exe at Thu Mar 29 16:12:40 2012.

Job cpu time: 0 days 0 hours 0 minutes 0.0 seconds.

File lengths (MBytes): RWF= 1 Int= 0 D2E= 0 Chk= 4 Scr= 1

————————————

my input

—————————————-

%chk=scanteste.chk

%mem=96MW

%nprocshared=4

#p bpw91 lanl2dz popt=(loose) scf=(maxcycles=400,xqc) scan test

Fe8multiplicidade 23 (2S+1)

0 23

Fe(Fragment=52)

Fe(Fragment=58) 1 2.48321456

Fe(Fragment=59) 2 2.48229999 1 70.52682490

Fe(Fragment=60) 3 2.48266151 2 70.53348065 1 60.00159306 0

Fe(Fragment=80) 3 2.48193698 2 70.54565842 1 -59.99973565 0

Fe(Fragment=81) 3 2.48269605 2 109.48955843 1 0.01169723 0

Fe(Fragment=54) 6 2.48167099 3 109.47510844 2 59.99181548 0

Fe(Fragment=61) 4 2.48279578 3 70.53467144 2 179.97879635 0

C(Fragment=1) 8 R 3 42.98654691 2 -44.81871601 0

R 2.01 9 0.01

What’s wrong?

I need you help

thanks in advance

Dear Joaquim,

Thx for your explanations which give us another look of how to perform a quick calculation.

I Have a little question. Can a scan be performed on a cartesian system ? Because I would like to perform a bonding scan on this (A,B,C, three atoms non bounded, for example)

A 0 0 0

B 0 0 3

C 0 0 -3

What I want to do : Scan between A-B 0.05 10 (so 10 points with a step of 0.05) and I want to freeze B (understand, only A is moving)

Is it possible ? If yes, how can I perform such type of calculation ?

Thanks in advance and continue to provide us all your experiences

Dear Camille,

Thank you very much for your kind words, I’m glad you have found useful information in this site.

Yes, you can perform a scan on a cartesian system. I’m just not sure of what you are trying to do. Are the numbers real or just an example? Bear in mind that you need to specify real numbers, so 0 is unacceptable, 0.0 is the correct way to do it.

I would use the opt=modredundant option with a Z matrix as described in the post but if you insist on doing it with cartesian coordinates you can always use the following input for such a simple molecule

X atom S x-steps x-size y-steps y-size z-steps z-size

This line should be included after the coordinates and a blank line. Example. Scan the atom number 5, two steps on the x axes with 0.02 increments, zero steps on the y axes and three steps on the z axes with 0.1 increments

X 5 2 0.02 0 0.0 3 0.1

I hope this helps! Have a nice day and thanks for reading!

Dear joaquinbarroso,

Thanks for the time you gave to me. Unfortunately it does not work (No variables in Z-matrix.). Here is my example:

%Nproc=16

%mem=20GB

%chk=AgI

#P PBE1PBE/GEN pseudo=read SCAN

Ag11Se12_I3_C3h

1 1

53 0.000000 0.000000 0.000000

…..

47 0.000000 0.000000 4.078987

47 0.000000 0.000000 -4.078987

blank line

basis

blank line

Pseudo

blank line

I want to do 10 steps on z (10 0.1) on the first atom (53=I), so what and where may I place the key words ?

Thx in advance again.

Dr. Barroso:

Estoy tratando de evaluar la distancia entre un yodobenceno unido a un dummy atom (yodo-X) y éste a su vez con el oxígeno de un formaldehido, es decir, Ph-I-X…..O=CH2. Necesito, primero hacer que el dummy atom tenga carga positiva (no se como hacerlo) y segundo hacer un scan de la distancia entre el yodo y el dummy atom.

Espero que sepas como ayudarme…

Muchas gracias!

Nice blog!

Desafortunadamente ambas cosas son imposibles de hacer. Definir un scan con un dummy atom y asignarle una carga específica a un átomo en particular no pueden ser definidos en Gaussian. El scan lo tendrías que hacer manualmente (cálculo a cálculo para cada punto que desees evaluar) mientras que asignar una carga al dummy atom no se puede. Podrías asignar una carga total al sistema pero como el dummy no tiene orbitales ni electrones entonces la carga no correspondería con la multiplicidad de las capas cerradas que tienes en las otras dos moléculas (el cálculo finalizaría de inmediato) Intenta con un átomo ghost (Bq) y que utilice la misma base que el resto, pero no creo que funcione por la misma razón que en el caso anterior.

Espero que te sea de ayuda

Muchas gracias, pense lo mismo pero como no soy usuario habitual de gaussian necesitaba de alguien que si supiera.

Buscaré otro medio.

Saludos!

Dr Joaquin,

i’m trying to search for the transtion state ot the clch3cl- and I get this error message

#2070

Is something wrong with the scan coordiantes?

thanks

Dr.joaquin

I have a problem . I’m trying to optimize Fe(2+) complex in M062X/6-311++G(d,p) level but PES is oscillate. How I can solve this problem?

thanks.

Hi ,

How to calculate excited state vibrational frequencies and optimized parameters using G03 or 09

Dear Sir,

Is it better to scan a molecule before optimization or after optimization.

HI

I think its better after opt

sir,

Geometry optimization and Scan both give minimum energy values. If we get minimum through geometry optimization, then what is the need for scanning after optimization. However if we scan the molecule and get a minimum, I think we should optimize the geometry on that configuration then. please suggest.

Hi Ahmed,

I guess the Dr. means that you need to start your scan from a reasonable minima, local or global..Local to find the global, global to make sure you R in the global minima…If you start your scan from a very messy geometry, most probably you will end with wrong results…

Cheers

hello sir

I want to know that how can we get information of transition from ground state to excited state is pi-pi* or sigma->pi* with G09.please help me.

Estimado Joaquin,

tengo una duda relacionada con un scan rigido en g03, sé que para que sea rígido debo poner la sentencia scan metodo y nosymm, pero mi duda está relacionada con que quiero definir un nuevo diedro, en el caso de modredun puedo ponerlo abajo despues de los valores de distancias, angulos y diedros, pero acá nosé como hacerlo.

de antemano gracias

Dear Dr.Jacquin,

I want to know what is the best method to describe the difference in reactivity between two electrophilic sites “will be attacked by a nucleophile”..Do you think NBO is a good candidate, or, Fukui functions, or any other similar method or population analyses…

Full of respect

Marawan

sir i run Opt=ModRedundant key words in my compound.what are factor taken from output.please tellme

Dear Marawan,

You need to use Fukui functions. The trick is what population analysis paradigm will you use as input for them! I suggest you run an NBO population analysis and use the population numbers you obtain as input for the Fukui functions the same way it is described in this blog. Find the corresponding post under the name “How to calculate Fukui functions”

I hope this helps. Have a nice day!

Dear Dr.Joaquin,

is it possible that a conformer which has a higher energy in calculations (by around 5 kj/mol) to be more populated in a gas phase experiment than another one which has a lower energy in calculations..???

Regards

Marawan

Hi joaq,

No words to thank you for your help, we so much appreciate your time and effort to help us.I am a research student working on gaussian03 and 09 for the first time. My aim is to model an iron crystal surface of around 24 Fe to examine the mechanism of synthetic gases on the surface; But the issue is i am struggling to calculate the SCF for months now but yet no success talk less of the optimisation and the rest. Please Please i will highly appreciate it if you can give me a helping hand on this.

Many thanks

Mala

Hello Sainna

Thank you very much for such nice words about this little blog.

24 Fe atoms sound like a lot! Have you tried performing a molecular mechanics optimization first? That is if you haven’t started from the actual X-Ray diffraction results.

Also I’d like to know what level of theory are you using, because if you use an all-electron basis set (3-21, 6-31, etc) you are calculating the structure for EVERY atom in there. Try using a pseudopotential or ECP like SDD or LANL2DZ. These will replace the entire core of electrons for Fe and approximate their effect on the valence ones (they even include relativistic corrections but since you are dealing with Fe they aren’t that significant) This will dramatically decrease the amount of time needed but it will still be pretty demanding if you don’t have a proper calculation platform. You have to place either one of those keywords (or any other if you want to use a different pseudopotential) instead of the basis set in your input file.

I hope this helps. Have a nice day!

First, great blog!

I would like to perform a relaxed PES scan in which two bond lengths change in concert. My problem is that if I request two variables to be scanned, g09 will calculate all permutations of possibilities but I only want the energies for when these two bonds are equal! Do you know how to request such a job? Thanks!

Hi Brad,

Thanks for your kind words about this blog. This feature you are requesting is unavailable in Gaussian. Therefore, you either do it manually or decrease the number of steps so its easier for you to find/read the diagonal elements of the permutation matrix you are obtaining.

One way to do it manually is to add several jobs to a single input file (–Link1–) and use geom=check addredundant and then change the value for the variable(s) of interest.

I hope this helps! Have a nice day

Hello,

I am a Ph’d student and I find this to be a very useful blog.

I trying to do a Relaxed PES scan over two dihedrals. The first step optimizes but the second fails. I am providing the input file for inspection. It completes the first optimization quickly but then fails on the second with the following error.

Atom 9 needs variable 9= 1.0207179493 but is 1.0202339455

Input z-matrix variables are not compatible with final structure.

%mem=50MB

%chk=alaninePESscan.chk

#b3lyp/6-311+g* opt=modredundant scrf( Read, pcm, solvent=CH3OH) maxdisk=16gb

alaninescan

0 1

1

6 1 R12

6 2 R23 1 106.889233

8 3 R34 2 A234 1 D1234

8 3 R35 2 111.843906 1 -42.755509

1 5 R56 3 107.419341 4 1.094446

7 2 R27 1 108.358101 4 103.974630

1 7 R78 2 110.929900 1 -57.519541

1 7 R78 2 110.164565 1 -175.092061

6 2 R23 1 108.342631 5 -96.556293

1 10 R12 1 A1_10_11 5 -54.048838

1 10 R12 1 A1_10_11 5 -D51_10_12

1 10 R12 1 A1_10_13 5 D51_10_13

Variables:

R12 = 1.09262603

R23 = 1.53445181

R34 = 1.20559737

A234 = 125.58993489

D1234 = 137.99568825

R35 = 1.35639477

R56 = 0.96944157

R27 = 1.45478305

R78 = 1.01555476

A1_10_11 = 94.44086879

D51_10_12 = 163.77452496

A1_10_13 = 139.36712081

D51_10_13 = 72.16917436

D * 2 3 * s 15 24.0

D * 2 7 * s 15 24.0

rmin=0.5

Ofac=0.8

Hi Jason,

Thanks for your kind words.

I’ll take a look at your input later on but it seems to me that the problem is the PCM keyword. Have you tried to do this scan in vacuo? Here is the thing, PCM will create a cavity by placing spheres on atoms (or groups of atoms) and then ‘gluing’ the outermost surfaces. But it doesn’t do it at every step. So probably after the first optimization (the starting point from which the cavity is generated) your molecule is trying to get out of the cavity because of the constraints imposed by the second step and this causes the PCM part to fail.

I’m not sure this is the case but if it is I would switch to other solvation model such as SCIPCM which re-calculates the cavity at every step of an optimization. Takes a lot of time to do it, though, so make sure you have a proper platform to perform it.

Hope this helps. Have a nice day!

Hi Jason,

Thanks for your kind words.

I’ll take a look at your input later on but it seems to me that the problem is the PCM keyword. Have you tried to do this scan in vacuo? Here is the thing, PCM will create a cavity by placing spheres on atoms (or groups of atoms) and then ‘gluing’ the outermost surfaces. But it doesn’t do it at every step. So probably after the first optimization (the starting point from which the cavity is generated) your molecule is trying to get out of the cavity because of the constraints imposed by the second step and this causes the PCM part to fail.

I’m not sure this is the case but if it is I would switch to other solvation model such as SCIPCM which re-calculates the cavity at every step of an optimization. Takes a lot of time to do it, though, so make sure you have a proper platform to perform it.

Hope this helps. Have a nice day!

The calculation runs without the use of wildcards. The wildcards cause the calculation to fail because the frozen dihedrals are not independent.

Dear Prof. Joaquin,

I am fairly new to using gaussian. I want to know 3 things

1. How to determine whether a uv-vis transition from TD-DFT is n -> pi*/ pi -> pi*/ d-d transition ?

2. How to study reaction mechanism steps using gaussian ?

3. Is there any open source software to analyse all the results of gaussian outputs ?

Sincerely,

Arnab

Student, Dept of Chemistry

Tripura University, India

Dear sir,

I am trying to do PES scan in guassain 03, I gave input file like this

%chk=OXINDOLE PES SCAN 1.chk

%mem=1MW

# p opt=modred b3lyp/6-311++g(d,p) nosymm geom=connectivity

Symbolic Z-matrix:

Charge = 0 Multiplicity = 1

N

C 1 B1

C 2 B2 1 A1

C 3 B3 2 A2 1 D1 0

C 4 B4 3 A3 2 D2 0

C 5 B5 4 A4 3 D3 0

C 6 B6 5 A5 4 D4 0

C 7 B7 6 A6 5 D5 0

C 8 B8 7 A7 6 D6 0

O 2 B9 1 A8 9 D7 0

C 3 B10 2 A9 1 D8 0

C 11 B11 3 A10 2 D9 0

C 12 B12 11 A11 3 D10 0

C 13 B13 12 A12 11 D11 0

C 14 B14 13 A13 12 D12 0

C 15 B15 14 A14 13 D13 0

C 16 B16 15 A15 14 D14 0

C 13 B17 12 A16 11 D15 0

O 12 B18 11 A17 3 D16 0

H 1 B19 9 A18 8 D17 0

H 5 B20 4 A19 3 D18 0

H 6 B21 5 A20 4 D19 0

H 7 B22 6 A21 5 D20 0

H 8 B23 7 A22 6 D21 0

H 11 B24 3 A23 2 D22 0

H 14 B25 13 A24 12 D23 0

H 15 B26 14 A25 13 D24 0

H 16 B27 15 A26 14 D25 0

H 17 B28 16 A27 15 D26 0

H 18 B29 13 A28 12 D27 0

Variables:

B1 1.78503

B2 1.51698

B3 1.50298

B4 1.41998

B5 1.41998

B6 1.54909

B7 1.41998

B8 1.41998

B9 1.20799

B10 1.33699

B11 1.51698

B12 1.51698

B13 1.41998

B14 1.41998

B15 1.41998

B16 1.41993

B17 1.41998

B18 1.20799

B19 1.04999

B20 1.09999

B21 1.09999

B22 1.09999

B23 1.09999

B24 1.09999

B25 1.09999

B26 1.09999

B27 1.09999

B28 1.09999

B29 1.09999

A1 99.41261

A2 110.99998

A3 120.

A4 120.00003

A5 114.5018

A6 114.50314

A7 119.99998

A8 129.728

A9 123.29839

A10 117.59999

A11 114.66671

A12 119.99995

A13 119.99996

A14 120.00003

A15 119.99999

A16 119.99879

A17 122.66512

A18 126.57188

A19 119.99879

A20 129.87568

A21 129.87515

A22 119.9988

A23 121.19869

A24 119.99876

A25 119.99872

A26 119.99877

A27 119.99816

A28 119.99874

D1 -1.91618

D2 159.645

D3 -144.20917

D4 -25.83842

D5 38.47666

D6 -25.83866

D7 -179.2926

D8 178.72377

D9 0.

D10 -180.

D11 0.

D12 179.42711

D13 0.

D14 0.

D15 179.42705

D16 0.61844

D17 -0.82464

D18 36.36379

D19 154.84087

D20 -140.84398

D21 154.73417

D22 179.40642

D23 0.

D24 -179.42702

D25 -179.42701

D26 -179.42707

D27 0.

D 3 11 12 13 0.0 S18 +10.0

I got the error message,WANTED AN INTEGER AS INPUT.

FOUND A STRING AS INPUT.

D 3 11 12 13 0.0 S18 +10.0

please help me to solve this problem.

Thanks in advance,

Devi.

As it seems in the input you have to leave a space between S and 18. Also try eliminating the + sign before 10.0.

Hope this helps!

Have a nice day

dear Joacquine, sorry for replying to such an old post.

I am performing relaxed hindered rotors calculations using g09, and I was encountering some problems, that’s how I ended up here.

It seems to me there is an error in your first post (at least using g09)

I didn’t read all the previous comments and answers (so, if my suggestion is already in there, just ignore it).

If I understood it right, the procedure you suggested in the very first post is the following

#HF/6-31G(d) opt=modredundant

[molecule specification]

blank

* 1 5 *

D 2 1 5 6 S 4 90.0

blank

I tried with gaussian 09 @HF/6-31G(d), just for curiosity (on the same machine)

well, the line * 1 5 * does not do anything.

It takes 5 minutes and (about) 10 seconds either including or not the * 1 5 *.

I checked the geometry convergence using molden (that unlike gaussview) shows the energy also for the intermediate steps and in both cases they look absolutely the same, with sharp peaks at each new step, with the geometry being comparable to the one you show in figure 2. In both cases, it takes about 40 steps to reach convergence.

If I use the

* 1 5 * R

D 2 1 5 6 S 4 90.0

(note the R in the first line), then I have exactly the behaviour you are looking for, and the molecule converges in about 20 steps.

It is not clear to me if this is a different behaviour between 03 and 09, but in 09 the procedure you showed seems not to work as expected.

Best regards,

Samuele Sommariva

Hello,

Thanks for the wonderful post. I have a question regarding something that seems to be contradictory with the g09 manual. There they write this regarding a relaxed PES:

——–

Wildcards in the ModRedundant input may also be useful in setting up relaxed PES scans. For example, the following input is appropriate for a potential energy surface scan involving the N1-N2-N3-N4 dihedral angle. Note that all other dihedrals around the bond should be removed:

* N2 N3 * R Remove all dihedrals involving the N2-N3 bond

N1 N2 N3 N4 S 20 2.0 Specify a relaxed PES scan of 20 steps in 2° increments

———

So instead of adding the angles they suggest removing them. Why the difference between your approach and theirs? Is it equivalent?

Dear Joaquin,

I read all the post and I would like to thank you for your comments. I have a question concerning G03. Indeed, I wanted to perform a PES scan but after few hours G03 crashed and the program closes immediately. I try by different way but I saw always the same problem. I’m beginner… Below you can find the “script”:

%chk=test.chk

%mem=64MW

%nproc=2

#p opt=Z-Matrix HF/6-31g nosymm

Molecule specification in Z-Matrix

* 23 25 *

21 23 25 30 0.0 S 18 +10.0

I tried different modifications such as SCF=QC, SCF=DM, with and without nosymm, … I saw always this problem when the program start the link502.

Thank for your answer.

Yours,

Eric.

Could you provide more info about the output file? What is the error message? I guess your scan works at first since it takes ‘a few hours’, right?

I need a little more info in order to help you.

Dear Joaquin,

I’ve been trying to perform a Scan with Gaussian 09 for tetrahydrogenfuran,

# b3lyp/6-31g Opt=ModRedundant

tetrahydrofuran

0,1

O 0.77079654 0.79787235 0.00000000

C -0.65957446 0.79787235 0.00000000

C -1.15511846 2.24007135 -0.0000010

C 0.09623254 3.10742935 -0.00000100

C 1.27278754 2.13726235 0.00000000

H -0.98373546 0.23889335 -0.91513200

H -0.98373546 0.23889335 0.91513200

H 1.90962654 2.24476135 0.91536000

H 1.90962654 2.24476135 -0.91536000

H 0.12761954 3.76662135 0.90190200

H 0.12762154 3.76662035 -0.90190500

H -1.78339746 2.44202535 0.90190200

H -1.78339746 2.44202435 -0.90190500

1 2 3 4 s 10 10.0

2 3 4 5 s 10 10.0

then I got this error

Iteration 99 RMS(Cart)= 0.00004405 RMS(Int)= 0.06810661

Iteration100 RMS(Cart)= 0.00004221 RMS(Int)= 0.06815077

New curvilinear step not converged.

Error imposing constraints

Error termination via Lnk1e in /share/apps/g09/l103.exe

Please could you let me know what I’m doing wrong?

Thanks

Monica

Dear Monica,

This kind of error usually relates to geometric constraints that increase the energy (or more accurately make it diverge) of the system. In this case, it seems you are trying to perform a Dihedral Angle scan on a cyclic molecule! The first 10.0° step is for sure trying to break the molecule apart. I don’t know what you are trying to calculate, maybe the energy of the bond breakage? if so, I would do it inversely: I’d take the open linear molecule as a diradical and then scan it so the final steps would point to the formation of the cyclic compound.

I hope this helps. Have a nice day and thanks for reading!

Thanks so much Dr. Joaquin for your quick answer! Well, I’ve tried to find different pseudorotation conformation of tetrahydrofuran ring. Please if you have a better idea how to do that let me know.

Thanks,

Monica

Send me an email with more details about what you want to do and I’ll be happy to help :)

Gracias a ti, Monica!

how can get pes of my molecule using gaus view 5.hear molecule contain 33 atoms.please give idea.

Quick question– I performed a relaxed scan to find a TS and found a structure for the TS and the product (with the TS looking about 5-10 kJ mol-1 higher in energy than the product, as expected). When I actually ran the thermochem on it though, they are the opposite way around. Is there any reason why gaussian would calculate this? (i.e. the energies from scanning not lining up those calculated by a freq calc?)

I have never come across this before and can’t figure out what this could be so any help would be appreciated

As an update to this, it looks like the H and G correction terms are the problem– They are higher in the TS which is causing the problem. It seems to me that they should be lower? Any thoughts?

Saludos cordiales estimado doctor, quisiera saber si es posible compilar varias salidas (outputs) de gaussian en una sola.

Cuando estoy realizando un scan relajado (PES) alrededor de un ángulo diedro y por ejemplo se detiene el cálculo por un corte de energía eléctrica, al restaurar el cálculo tendré un nuevo archivo de salida (output2) y si se volviere a cortar el cálculo por alguna otra razón generaría otro archivo de salida (output3). Esto lo realizo usando la palabra clave restart, y el cálculo reiniciará tomando el check hasta ese momento desarrollado y continuará desde el punto en que se quedó. Entonces es posible compilar todos los archivos de salida como uno solo o esto no se puede realizar. Disculpe doctor si la pregunta es un poco tonta y sin sentido, de antemano muchas gracias por su ayuda y felicitaciones por su blog, que es muy bueno y de mucha ayuda.

Hola Christian,

Gracias por leer el blog y por tus palabras sobre él. Acerca de tu pregunta, no creo que Gaussian pueda ‘pegar’ los archivos de salida de manera automática. Yo lo que haría sería pegarlos manualmente; asumo que usas linux y en la linea de comando escribes cat file1 file2 file3 > filecomplete

de esta manera se unirán todos los archivos en el orden que tu le indiques. Para que gaussview, por ejemplo, supiera que es un solo archivo quizá habría que editar un poco cada archivo, pero no lo sé.

Espero que te haya sido de ayuda. Saludos y que tengas un excelente día.

Gracias estimado doctor por su oportuna respuesta, como usted supuso uso linux cualquier resultado se lo comento. Y la verdad, su blog es de mucha ayuda.

Saludos estimado doctor, le comento que realicé la compilación de los archivos como me sugirió y gaussview abrió el archivo sin inconvenientes, así que muchas gracias por la ayuda y a futuro seguro molesto con otras preguntas, ya que Gaussian a veces resulta impredecible.

Saludos estimado doctor,

El motivo del presente, es preguntar acerca de un error recurrente que tengo en un output de Gaussian y albergo la posibilidad de que talvez me pueda brindar su ayuda.

Estoy realizando un PES scan de una molécula y en uno de los ángulos diedros en lo que desarrolla el cálculo, arroja el siguiente error:

Small interatomic distances encountered: 24 4

Problem with the distance matrix.

Error termination via Lnk1e in /home/calculos/Documentos/Gaussian09W/tar/g09/l202.exe at Sat Mar 9 11:16:39 2013.

Los átomos involucrados son los siguientes C12_C10_C15_C2 y la entrada que estoy usando es la siguiente:

%NProcShared=6

%mem=12000MB

%chk=Scan7.chk

# opt=addredundant b3lyp/6-31+G* Nosymm

Metoxi-Cromona Ale Structura B Scan7 C12_C10_C15_C2

0 1

C 4.74828600 1.19411500 -0.01431300

C 2.57701000 0.06499200 0.00280300

C 3.26956700 -1.15916600 0.10633700

C 4.65070400 -1.21613000 0.14840600

C 5.38964300 -0.02687500 0.08707500

H 5.30153200 2.12414500 -0.06309200

H 2.71402100 -2.08612400 0.15591000

H 5.15419700 -2.17172200 0.22859100

H 6.47316200 -0.06058300 0.11905700

C 0.29001700 -1.04206100 0.00418000

H 0.75908200 -2.00670900 0.07868600

C -1.08449000 -1.01355100 -0.04387900

C 3.34723500 1.26047500 -0.05755100

O 2.79113200 2.47646800 -0.15612300

C 1.10064100 0.14954200 -0.04370600

O 0.55183100 1.28430200 -0.12837700

N -1.83412100 0.09997600 -0.11263800

C -1.82736600 -2.34466400 0.00395800

F -2.62578000 -2.42468500 1.09738000

F -2.62488300 -2.50211700 -1.07890600

F -1.00426700 -3.40195600 0.04313900

H -1.27869600 0.95488800 -0.19159500

C -3.27793700 0.23063300 -0.29078400

H -3.52733900 0.23797000 -1.35795600

H -3.80191200 -0.60273200 0.17463100

C -3.75642400 1.52588300 0.34358200

H -4.84658300 1.60742500 0.22003300

H -3.53086900 1.52760000 1.42019900

O -3.09591000 2.59762100 -0.30097800

H 1.81045400 2.33626800 -0.17310300

C -3.38146900 3.86407800 0.26791300

H -3.07416500 3.90542000 1.32120400

H -2.81631500 4.60347600 -0.29855900

H -4.45228000 4.09891400 0.19964000

12 10 15 2 0. s 12 30.

De antemano muchas gracias por cualquier ayuda.

Hola Christian.

Durante el barrido lis átomos 4 y 24 llegan a quedar demasiado juntos sin que la estructura se pueda relajar, modificala para que no choquen.

Saludos

Gracias Doctor, por su ayuda logré solucionar el problema con su ayuda. Saludos cordiales.

Dear Sir,

what is the difference between CIS and td

how to do a excited state optimization using b3lyp/6-31g(d,p)

how to calculate the UV and fluorescence spectrum using Gaussian 03.

please help me in the route section for the above items.

thank you sir.

Hello Dr.Barroso,

I did a rigid potential energy surface scan between one molecule of carbonic acid and a water molecule, varying the distance between O(h2co3) and H(h2o) from 1.8 to 1.9 in steps of 0.02 using mp2 method and cc-aug-pvqz basis set. But, I found that instead of the change in the intermolecular distance

there is a change in the intra-molecular distances.

Later,I did the scans manually for each point to get the corresponding energy.

Can you please tell me why such changes in intra-molecular distances happen and how I can avoid that?

Thanks in advance,

ph

Are you feeling like you do not have the car you

deserve? This is because the salesman is trying to get a great

commission out of you. Don’t treat a salesperson like a friend. Continue reading for some valuable tips and information.

Hola Joaquín. Mi nombre es Andrea y quiero hacer un scan alejando dos moléculas. El problema es que en la segunda molécula debo alejarla desde su centro. Pensé en usar un dummy atom, pero gaussview no me permite utilizar el dummy atom en el editor de coordenadas redundantes. Me aparece dummy atom not allowed. Sabes de alguna solución a este problema?

Hola Andrea. No hay. Lamentablemente el uso de átomos dummy no esta permitido en el scan. Sugiero que lo definas mediante alguna distancia ínter molecular o de plano a mano haciendo varios cálculos a varias distancias

Saludos y que tengas un buen día!

hello Dr.barroso,

i would like to convey my greetings…ur blog is really very informative and usefull to understand computational chemistry..

Thank you very much Anitha! Your kind words encourage me to continue blogging :-) you made my day.

I’m glad to know you’ve found my site useful!

Have a wonderful day :-)

Thank you very much, Anitha! I hope this site continues providing you with help for your research.

Hi Joaquin,

Just found your nice blog and am surprised to find you still answering questions about PES scans three and a half years later. Thanks for all your help. How about one more question?

I have trying to do a scan with g09 but i’m running into a very frustrating problem. My input is as follows:

%NProc=4

# MP2/aug-cc-pVDZ Opt=ModRedundant

title

1 1

C -1.914 3.074 -1.879

S -0.773 1.939 -1.044

C 0.064 3.082 0.089

C -1.908 1.108 0.100

H -2.610 2.516 -2.505

H -2.488 3.652 -1.154

H -1.361 3.768 -2.512

H 0.696 2.528 0.783

H 0.691 3.775 -0.471

H -0.661 3.658 0.664

H -1.349 0.480 0.795

H -2.482 1.834 0.675

H -2.603 0.476 -0.452

5 1 2 3 -180.0 S 12 15.0

If i run this using gaussin03 on an older computer, it works no problem. However, to speed things up, i’d like to use a different cluster that is running g09. The same input using g09 on the newer cluster doesn’t perform the scan. Instead, the log file mentions:

The following ModRedundant input section has been read:

D 5 1 2 3 -180.0 B

Therefore, only a geometry optimization is performed for the molecule.

Have you ever seen anything like this or can you spot any errors in the input? It’s driving me crazy that this is working on an older version of gaussian but not g09. Thanks in advance,

ryan

Hello Ryan,

I’ve never seen something like this but I’d suggest you to try adding the following line before the scan definition

* 1 2 *

so it looks like

* 1 2 *

5 1 2 3 -180.0 S 12 15.0

leaving the necessary blank lines, of course.

Maybe your system is hindered and needs to have all related dihedrals defined as one.

Try it and let me know how it went, ok? if it doesn’t work we’ll try something different.

Have a nice day!

Hi!

I am quite new to Gaussian, so I am sorry for such a simple question.

Anyway, the problem is, how can I do relaxed scan sequentially with different basis sets? For example I would like to get PES for ethane, as it is described. But when I change angle, I would like to do optimization with RHF/3-21g and after that with B3LYP/6-31g, than change angle with next step and again do optimization with RHF/3-21g and after that with B3LYP/6-31g. The idea is that I have a lot of steps and can save some computational time doing it that way. And even if it is not true, I would like to know how to do this.

If I use –Link1– it does PES using RHF/3-21g for all angles, and after that optimisation with B3LYP/6-31g only for last point of PES scan.

Thanks!

Hi Denis!

If I understand correctly, what you’d like to do is having the same PES scan done twice, with 2 different levels of theory, in one go. Am I right?

Well, if I’m right then I’m sorry but there is no way to do it, you have to actually launch both jobs separately, or use the –Link1– method and define the second scan from the beginning in the second part of the input but you might as well just send them separately, which I recommend since it would be easier to visualize separate results that way.

Have a nice day!

Dear Joaquin,

Thanks for such fast and professioanl reply!

Maybe it is really hard to explain, but I would like to concern that we understand each other perfectly well.

Imagine this:

1. I want to do PES scan for relaxed system (or constrain some distances/angles – doesn’t matter), and I have a step of 10 degrees for dihedral angle. And I want to get dependence of Energy from this Angle change.

2. As far as I know, when system has done optimisation for some angle it does a step of this angle for 10 degrees and starts geometry optimisation once again. Usually (if I understand correctly) it makes larger moves for molecule conformation if I use some simple levels of theory (like RHF/3-21g) and smaller but a lot of moves for higher levels of theory (like b3lyp/6-311+g). So, what I’m trying to say is that rough RHF will make for example 3 large steps and will be near really good defined local minima, but precise b3lyp will make 10 smaller steps and will be near really good defined local minima (closer of course), so it will take much more time. ( – is it right?)

3. What I want is – to make rough approximation to local minima using RHF and to finish it with more precise b3lyp, but not to do it with b3lyp from the very start. Then, go to next step of 10 degrees and do it once again. (it’s impossible?)

If it is not possible, what will you reccomend for getting more precise results but saving computing time?

If I use –Link1–, everything that was written before this point will be done first and only after that will go to the part after –link1– string? Will it put everything to check file before starting things that are written after –link1– string?

Sorry for having a lot of words and questions in this one post =) And thanks once again!

Hello again, Denis

Unfortunately what you want to do is not possible with Gaussian, but more than that, it is conceptually wrong.

The ‘step size’ during optimizations (how much will degrees of freedom -lengths, angles, dihedrals- will change in order to search for the minimum) does not depend on the level of theory employed but on the algorithm for finding the minimums; the default on Gaussian is Berny’s algorithm. What actually changes from one level of theory to another is how well are those minimums described, thus, the potential energy surface (PES) is different for the same molecule at various levels. If you take the example of this post -ethane- and rotate the C-C bond and plot the energy vs the angle you will find the classical curve you can find in any organic chemistry textbook; and if you do it at a higher level of theory you’ll find the same qualitative plot but with different numbers. The problem is that for more complicated molecules, some minimums may be completely overlooked by a low level of theory so exploration of the PES is level dependent.

In other words, what you are trying to do is to come closer to a minimum by using two different levels but the reality is that both levels yield different minimum (by virtue of describing two different PES). So, which one is the real PES? usually the one obtained with a higher level since, according to the variational principle, the lower the energy the closer it is to the true value.

I recommend you use the low level for a qualitative description and the higher one for a quantitative one. The question is how large is your molecule? and how many computational resources do you have available. It is always a compromise between both.

Best regards

Thanks a lot! =) Now I feel better, knowing that it’s impossible and will read further about Berny’s opt.

Thx!

Dear joaquinbarroso,

This is a really wonderful blog site, I’ve learnt a lot here.

I am new to Gaussian, and I have some questions about PES scan of a dihedral.

I intend to rotate a methyl group attached to a cytosine. For example dihedral 1-2-3-4

atom 1 and 2 are on the ring, 3 is the carbon and 4 is the hydrogen.

I can used the wildcard * 2 3 * R to rotate the three hydrogens on the methyl group simutaneously. However, i want the position of atom 1 2 to be fixed when gassian try to rotate the dihedral 1-2-3-4, since atom 1 and 2 are in a ring. But I want the possition of atom 1 and 2 could be relaxed in the subsequent process.

How should I specifie the wild cards?

Does

1 2 3 * R (why an R here? The gaussian official website said R should be placed here, but you did not mention this? )

1 2 3 4 S 36 10

will this work ?

Junyan Lu

Hello Juanyan

Thank you very much for your kind words about my blog, I’m glad it has served you well.

I’m not sure about the R, perhaps it means rotate? I have no idea but you can try with the S and with the R, one of them will might end in an error and then you’ll know which is the right way.

I think the 1 2 3 * option is the right one because if you scan 1 2 3 4 it will (at best) rotate the H4 without any consideration of it clashing with the other two H on that methyl group.

I will take a look on the R thing and probably do an update.

Have a nice day!

Estimado Prof. Barroso,

Mi nombre es Leandro, soy un estudiante de doctorado argentino. Antes que nada le agradezco mucho por su blog y por su predisposición para ayudarnos. Estoy dando mis primeros pasos con el modelado y me surgió un problema que no estoy pudiendo solucionar. Quiero realizar un scan con optimización (para el ángulo diedro 13-14-18-19) usando el comando ModRedundant con AM1, pero no funciona.

El comando que uso es de la forma:

# am1 Opt=ModRedundant Test

(Z-Matrix)

13 14 18 19 S (nº pasos) (paso)

El cálculo comienza pero después de unos minutos se detiene y aparece “Test job not archived”

Si se le ocurre qué puedo estar haciendo mal le agradecería mucho que me ayude

Gracias nuevamente

Saludos

Hi, i think that i saw you visited my weblog thus i

came to “return the favor”.I’m attempting to find things to enhance my website!I suppose

its ok to use a few of your ideas!!

Thank you for the good writeup. It in fact was a amusement account it.

Look advanced to more added aagreeable from you! However, how can we communicate?

Dear Joaquin,

i am a korean student. In my chemical class, we are using gaussian 03 & gauss view.

i optimized H3+ ion using MP2/6-311G(d,p) method. but i don`t know how to calculate H3+ ion transition state and h3+ ion transition state structure. please tell me about this problem.

Way cool! Some extremely valid points! I appreciate you penning this article plus

the rest of the site is extremely good.

Hola Joaquin, te felicito por este blog. Lo consulto en ocaciones.

Tengo una pregunta. ¿Es posible realizar scan’s relajados con CCSD(T) en g09?

There are a number ways that you can improve the overall health of the

teeth. The quick change of diet restores the pets

coat to a lush and glossy fur that does not have the shedding problem.

Just keep reading to learn more about how discount veterinary care

can help you care for your pets.

Hey there, You must have done a very good task. Let me surely bing that for my own portion highly recommend in order to my buddies. I believe they’re going to be taken advantage of this amazing site.

Hello, after reading this awesome post i am too delighted to share my experience here with

colleagues.

Have you ever thought about creating an ebook or guest authoring on other websites?

I have a blog based on the same subjects you discuss and would really like

to have you share some stories/information. I know my viewers would appreciate your work.

If you’re even remotely interested, feel free to shoot me an e-mail.

I am getting extremely frustrated by trying to scan a simple bond length. I have tried relaxed and rigid. The molecule often falls apart – a very similar molecule was sucessfully modeled in the literature so I feel there is some syntax issue. I am listing one of my .com files that will optimize but not scan. What am i doing wrong?

I have tried diffent coodinate systems and creating files both with help of Gaussview and Molden.

%chk=test_method.chk

%Mem=32GB

%NProcShared=8

#p opt=modredundant am1 geom=connectivity

***

0 1

N 1.11548900 -0.02669100 1.24033100

C 1.73535800 1.02454700 0.53187700

C 1.51301200 2.39000000 0.61737900

H 0.75065700 2.80773500 1.26331100

C 2.30777000 3.23101200 -0.16717000

H 2.14702300 4.30150900 -0.11955900

C 3.28755100 2.72036000 -1.00904000

H 3.88564600 3.39174900 -1.61252000

C 3.49413000 1.33860600 -1.08758400

H 4.25054300 0.93365500 -1.75194400

C 2.71168500 0.50008900 -0.32032100

C 2.73366300 -1.00913300 -0.16785700

C 1.33053000 -1.23196100 0.46265100

O 0.39906900 -1.27989000 -0.65761300

C -0.86166500 -0.84369600 -0.45829200

C -1.78113100 -1.10005400 -1.48199700

H -1.43961700 -1.63340700 -2.35988500

C -3.09044300 -0.67932000 -1.36166800

H -3.82040400 -0.86609400 -2.13755200

C -3.46923200 0.00118000 -0.20945900

C -2.57142200 0.27398000 0.80880200

H -2.90678000 0.82021300 1.68261400

C -1.25373900 -0.14581900 0.69404600

C -0.22726500 0.12775000 1.76679200

H -0.34383900 1.13567100 2.16623200

H -0.35817100 -0.55961800 2.60874900

C 3.83211900 -1.39494200 0.83619700

H 4.79998500 -1.07340100 0.44662500

H 3.86702600 -2.47752200 0.98536600

H 3.67454300 -0.90780600 1.80099000

C 2.93388700 -1.76352000 -1.47896800

H 3.95202500 -1.59903000 -1.83986800

H 2.23286400 -1.42897000 -2.24165800

H 2.80548300 -2.83968200 -1.33236900

C 1.13459500 -2.50227600 1.26651800

H 0.08100700 -2.62770400 1.52688500

H 1.72477900 -2.47135100 2.18222200

H 1.43477900 -3.36431100 0.66834300

N -4.85907400 0.45578200 -0.07203700

O -5.62521600 0.20553400 -0.97881700

O -5.15996200 1.05359400 0.94074700

1 2 1.0 13 1.0 24 1.0

2 3 1.5 11 1.5

3 4 1.0 5 1.5

4

5 6 1.0 7 1.5

6

7 8 1.0 9 1.5

8

9 10 1.0 11 2.0

10

11 12 1.0

12 13 1.0 27 1.0 31 1.0

13 14 1.0 35 1.0

14 15 1.0

15 16 1.5 23 1.5

16 17 1.0 18 2.0

17

18 19 1.0 20 1.5

19

20 21 2.0 39 1.0

21 22 1.0 23 1.5

22

23 24 1.0

24 25 1.0 26 1.0

25

26

27 28 1.0 29 1.0 30 1.0

28

29

30

31 32 1.0 33 1.0 34 1.0

32

33

34

35 36 1.0 37 1.0 38 1.0

36

37

38

39 40 2.0 41 2.0

40

41

B 13 14 S 10 0.150000

Hello

What is the error you are getting? Your input looks fine but maybe there are some dependencies you havent taken into consideration. Try re building the molecule or add the lines before the scan * 13 and * 14

I will try to check it myself and get back to you.

Have a nice day

thanks a bunch for fast response. I manually scanned but am having another vexing problem. I want to scan a bond-length, and optimize at each step with the bond distance constraint. In Gaussview, I fix the bond length and let it translate. I also fix atom position for two atoms in bond using atom editor by placing “-1” for the two atoms. Retain, save and run.

After optimization, either the molecule relaxes to starting bond length (lowest energy) or the bond breaks. Maybe this is the chemical information I am looking for on bond stability, but the bond is breaking at unreasonable short distances.

I have also tried running as Z-matrix and using appropriate nomenclature or in cartesian coordinates. I am trying to find TS and this bond length tragectory is vexing me. I must be doing something fundamentally wrong.

thanks again. I am having more luck with a conformational analysis and have dropped the project causing me problems for a little while.

Hi Joaquin

I am trying to run a relaxed scan on one of the dihedrals in my molecule, with the ultimate aim of getting a full PES for it. I think there must be something wrong in my input because I don’t get the table of energies that you describe. Any chance you could take a quick look and see where I’m going wrong?

The calculation terminates normally and I get the magic “4 Yeses”, so it’s doing something…

%NProcShared = 12

%chk = sarin_PES.chk

#p B3LYP/6-31G(d) Opt=ModRedundant pop=ESP NoSymm

Sarin-test

0 1

C

P 1 r2

O 2 r3 1 a3

F 2 r4 1 a4 3 d4

O 2 r5 1 a5 3 d5

C 5 r6 2 a6 1 d6

C 6 r7 5 a7 2 d7

C 6 r8 5 a8 2 d8

H 1 r9 2 a9 3 d9

H 1 r10 2 a10 3 d10

H 1 r11 2 a11 3 d11

H 6 r12 5 a12 2 d12

H 7 r13 6 a13 5 d13

H 7 r14 6 a14 5 d14

H 7 r15 6 a15 5 d15

H 8 r16 6 a16 5 d16

H 8 r17 6 a17 5 d17

H 8 r18 6 a18 5 d18

Variables:

r2= 1.8382

r3= 1.4876

a3= 113.92

r4= 1.5846

a4= 109.95

d4= 125.01

r5= 1.6306

a5= 103.26

d5= 233.68

r6= 1.4185

a6= 123.34

d6= 258.03

r7= 1.5374

a7= 111.54

d7= 72.68

r8= 1.5337

a8= 107.21

d8= 195.23

r9= 1.0889

a9= 109.69

d9= 72.28

r10= 1.0891

a10= 109.55

d10= 312.67

r11= 1.0883

a11= 110.30

d11= 192.63

r12= 1.0917

a12= 108.96

d12= 311.97

r13= 1.0898

a13= 111.68

d13= 295.45

r14= 1.0904

a14= 110.83

d14= 176.30

r15= 1.0903

a15= 111.32

d15= 56.44

r16= 1.0903

a16= 111.30

d16= 298.62

r17= 1.0903

a17= 111.30

d17= 178.13

r18= 1.0905

a18= 110.82

d18= 58.41

* 6 7 *

13 7 6 8 S 11 30.0

Hi James,

At a first glance, maybe you could change the order in the last two lines so the numbers in the middle get the same order:

* 6 7 *

8 6 7 13 S 11 30.0

There could also be a problem with atom clashes. Do you get an error or something like that? Does your calculation finish correctly? Try my previous suggestion and let me know how it went.

Have a nice day

Hi Joaquin,

All of my calculations are finishing correctly without an error. I tried the change you suggested and I am getting the same result. I have also tried it on a different dihedral with the same result.

Thanks very much for your reply.

James

Hi there

I think I know what may be happening. The results summary table is only available for rigid scans. For relaxed scans you have to extract the energy profile manually. Try searching for every instance of the ‘Optimized parameters’ string and search for the energy ‘SCF Done’ slightly above the previous string.

Let me know if I’m right please. Have a nice day

This has also been my experience and it is an annoying aspect wrt to making PES. You can use programs like TextCrawler to pull out data but takes some effort.

Hey Joaquin

I *think* I fixed the problem by removing the wildcards. When I was running with the wildcards, I was only getting one “Optimization completed” string, so for whatever reason, it wasn’t doing a scan. I have run quite a few calculations without the wildcards and have had no atom clashes and the results summary table has appeared.

Thanks for your help.

Oh. And William, I am developing a Python script to pull the data from the output files. It’s still very much in Alpha stage and only works for my specific molecule at the moment, but when I have a more complete script, you’re more than welcome to have it if you like.

that would be wonderful. Recently, I wanted to perform a PES but using chemical shift rather than energy for the z axis. It would be nice to have a simple script that would create the output table you get with Scan function.

Well as I say, it’s still in very early development, and I’m very much a beginner at coding, but feel free to email me

Hi Joaquin

If you see this comment, could you do me a favour and delete my comment with my email address. It was foolish of me to put it up and I am getting huge amounts of spam!

Hi James, as per your request, I have deleted your email address from your comment. I’m sorry you got so much spam from posting here. Maybe try to use ‘at’ and ‘dot’ (with letters and not symbols) so no spambots can detect your address.

I hope this helps

Estimado Joaquin,

Saludos cordiales y muchas gracias por el tiempo que se toma para responder a las preguntas.

Me podría dar unas pautas para extraer los datos de un PES scan 3D y realizar el gráfico de la superficie de energía potencial con programas como Origin u otro que me recomiende.

No encuentro la manera de extraer los datos del archivo .log del cálculo.

De antemano muchas gracias.

Saludos

En el caso de los cálculos scan, hay un pequeño sumario al final del cálculo: “scan summary”. Si tu calculo es relajado entonces tienes que buscar la energía de cada punto antes de la leyenda Optimized parameters para cada uno de los pasos. De ahí pasarlo a Origin es a mano.

Saludos

excellent put up, very informative. I wonder why the other specialists

of this sector do not notice this. You should continue your

writing. I am confident, you have a great readers’ base already!

Dear Sir,

I tried to perform PES scan for three selected dihedral angles at a time. It is possible to draw 4d gaussian grid plot in gauss view. please kindly send the reply.

Not with GaussView but probably with Origin by Oracle you could get the corresponding surfaces. That was one bold calculation!

Have a nice day!

Good day

Thank you for your wonderful and informative blog :)

My questions are slightly related to tdcourtney’s question (see bottom).

1) If I rotate the dihedral as follows:

* 1 5 *

2 1 5 8 S #steps #increment

All atoms connected to atoms 1 and 5 will rotate by the #increment for each #step. However, will these atoms be allowed to be optimized w.r.t the dihedral angle or will they effectively become frozen w.r.t the dihedral angle, as well? That is will the dihedral be able to change for the atoms connected to atoms 1 and 5 during the optimization, except for dihedral 2 1 5 8 ?

2) If I want to perform a relaxed PES calculation 2,2-bipyridine (as an example), where the dihedral is set for N(1) C(6) C(7) N(12), would it be possible to allow the entire ring to rotate as well in order to minimize the number of calculations required, but still allow the ring to be fully optimized afterwards? (This is also related to question 1) How would this be done?

2,2-Bipyridine: Numbered from N(1) around the ring and then across to the next ring, finally ending on the other N(12). C(2) being a secondary carbon and C(6) being a tertiary carbon. Hydrogens have 1 at the end of the number.

2 6 7 *

3 6 7 *

4 6 7 *

5 6 7 *

21 6 7 *

31 6 7 *

41 6 7 *

51 6 7 *

* 6 7 *

1 6 7 12 S 36 10

tdcourtney |

“Is there anyway to rigidly move the hydrogens with the oxygen in between steps and then release them to be optimized?”

Kind regards

Cameron

I think that I’ve found a problem with the method I proposed for bipyridine. The atoms on the opposite side of the ring (from Nitrogen) will have a dihedral of the current dihedral (step) +/-180 degrees, i.e. for C(4), C(5), H(41), H51. However, C(3) and H(31) shouldn’t be a problem, since they should be ‘in line’ with C(6)-C(7) ?

How can I solve this?

I have never performed a scan around a bond conecting two rings, however it should be as I said before. Don’t use the wildcards and just scan around 1 6 7 12, the rest of the atoms will follow and will be optimized. I think it would be interesting to plot the C6-C7 bond length as a function of the dihedral angle and try to find a correlation in terms of pi-pi overlap between both rings or pi delocalization from one to another.

Let me know if it worked.

Dear Cameron,

Thank you for your kind words about this little blog of mine and please excuse the lateness of my response.

When you define all dihedrals around a bond with * 1 5 * (as in your example above) all atoms attached to the 1-5 bond will rotate at the same time with each step of the scan, however they will NOT be frozen during the optimization on each step. Only the dihedral 2 1 5 8 (as per your example above again) will be frozen.

In the case of bipyridine you don’t need to use the wildcards because given the connectivity the rest of the ring must rotate as a group. Have you observed the opposite?

I hope this helps, have a nice day!

Dear Joaquin,

Thank you for your response.

As you said the ring is supposed to rotate along the the N-C-C-N dihedral, however I’m not sure if Gaussian considers this or breaks the ring similar to your Fig. 2, i.e. only the atoms 1,6,7,12 and the attached atoms (to 6,7) 5 and 8 will rotate leaving the other atoms 2,3,4 behind. I know that if you set the dihedral using GaussView the ring does rotate as one would expect. I added redundant coordinates for the atoms 2,3,10,11 [in the format :D # 19 6 # (12 or 1) R] (for the hydrogens, as well) (eg. D 2 6 7 12) so that the atoms on the same side of the N rotate as well. There was only a small improvement in performance using this method. The results are the same regardless.

Bipyridine was only an example, but you’re right it would be interesting to look at the C6-C7 bond length as a function of the dihedral.

Kind regards

Cameron

Hello Cameron,

Did it work? I’m curious to know if the modredundant method is hindered for this kind of rings. It shouldn’t be. In GaussView there are various options for moving dihedral angles (rotate groups, rotate atoms, translate atoms, etc.) take a look at the proper option.

Have a nice day!

Yes, it works. However, I went back to using the following modredundant keywords, since there wasn’t much of an improvement:

* 6 7 *

D 1 6 7 12 S 36 10.0

The symmetrical system that I’m currently working with can be described using 4 dihedral angles, cf. Di-(2-picolyl)amine. However, some of the permutations of the aforementioned dihedral angles result in overlapping atoms (or atoms in very close proximity, for example. dihedral: N-C-C-N : 0; C-C-N-C:0;N-C-C-N : 0; C-C-N-C:0) which results in a ‘Problem with the distance matrix’ error when Gaussian reaches this step. I then proceeded to write a small program in order to generate input files for each unique permutation in order to allow the calculation to continue past this error. I’m using a batch method for the submission of jobs.

Have your ever used Optim2 (http://www-wales.ch.cam.ac.uk/~wales/OPTIM2/node1.html) or other software in order to determine/estimate the position/geometries of minima, saddle points, maxima? I wonder if this can be done for a 5D system as well (4 dihedrals + Energy).

I like to know how to perform geometry optimization and potential surface scan at once using GAMESS…?

Thank you,

I’m sorry Iresh but I have little experience with GAMESS.

Have a nice day

Hey there! I’ve been reading your blog for

a long time now and finally got the bravery to go ahead and

give you a shout out from Atascocita Tx! Just wanted to mention keep up the great job!

Hey! This is a great blog.

I have a question. Suppose i want to find the rotational barrier in ethane or similar molecules, should i report the rotational barrier as the highest energy that i obtain on the PES graph or should i report the value that i obtain from ts calculations (QST2 or QST3).

Thank you in advance.

I understand that PES scans, to a small extent change the co-ordinates of other atoms to optimize the structure, so the energy that i obtain at the peak of PES graph is not going to be the same as TS state calculations. So i think it should be TS calculations energy value that should be used for reporting rotational barrier. I hope I am right ( :-|)

Dear Joaquin,

i had ask you about 2D angle scan and you explain it very well.

But i have some more question

1. suppose i have two molecule and i want to do 1D angle scan and i am doing this using opt=modredudant keyword for every 10 degree what should i do if i want just single point energy over every step instead optimized ? since i am not getting the result as i am expecting.

2 Is it any difference between the doing scan using z-matrix and opt=modredudant (redundant coordinate) keyword through Gauss-view ?

3. suppose i want to do scan using z-matrix but angle that i am going to use for scanning is not connected through bond then what should i do? since for this we have to add scan co-ordinate manually and i dont know how to do this.

4.During the scan how many coordinate we can freeze for publication?

Thank you

Vishal

Pingback: Gaussian 09 User Manual Pdf - Queerspectives

Pingback: Manual Gaussian 09 PdfPDF Online | Manual PDF

Pingback: Manual Gaussian 03 PdfPDF Online | Manual PDF