Dear Tahereh,

With the scan keyword you will carry out what is called a rigid scan, this means that the program performs single-point calculations for different structures that vary ONLY the geometrical parameter you selected. This yields a very rough estimation of the potential energy surface, given that you are neglecting all possible relaxations that can take place upon modifying that coordinate. It is usually better (and, for your system, it is perfectly feasible) to perform a relaxed scan, in which the rest of the molecule is optimized at each step of the scan. This can be achieved with the opt=modredundant keyword.

To use this keyword, it is easier to work with Cartesian coordinates instead of z-matrix, just to simplify the input, which will look like this:

-----------------------------------------------

# opt=modredundant wb97xd/aug-cc-pVTZ

\blank line

title

\blank line

2 1

C coord_x coord_y coord_z

H coord_x coord_y coord_z

H coord_x coord_y coord_z

H coord_x coord_y coord_z

Cl coord_x coord_y coord_z

\blank line

B 1 5 S 20 0.10

\blank line

-----------------------------------------------

If you take a look to the end of the input, the last written line says that for the bond (B) between atoms 1 and 5 (which correspond to C and Cl according to the Cartesian coordinates) you want to perform a scan (S) of 20 steps, with a step size of +0.10 angstrom (yes, you can do scans that reduce a distance/angle/dihedral just changing the sign of the last number).

Now, you can see that I have removed the td keyword. Given that TD-DFT excited states do not interact with the ground state (this is not CASSCF or any multireference method), the ground state PES will not be modified by them. And forcing the program to calculate the excited states at each step of the intermediate optimizations will increase the computational cost a lot. I think it is better to perform the scan just in the S0 PES, and once it is finished, carry out TD-DFT single-points at each point of the scan. While it involves submitting more calculations, the computational cost will be reduced.

You could have performed the relaxed scan optimizing the geometry not in the S0 surface, but on an excited state, of course. In that case, you need the td keyword, specifying the excited state you are interested. For instance, if you want to do the scan in the S1, the keyword should read td=(nstates=1,root=1). Again, you do not need the upper excited states, unless there are degeneracies.

Hope you find it helpful

Dear Enrique Manuel Arpa thanks for your explanation, I have one more question in case we want to do a scan job over a bond coordinate, I'm not sure how we can extract the energy at each steps from the output file?

Tahereh Alavi ,

If you are working with Gaussian, then you should use its graphical interface, GaussView. There you can visualize the scan output, with all the steps, and in Results → Scan you have all the information you need. For instance, you can visualize a 2D plot (or 3D plot in case you scanned two coordinates simultaneously) with the total energy with respect to the coordinates, and then save the data to a txt file.

If you want to look for them in the output file directly, it is much more complicated. You need to perform the calculation with #p in the keyword line (for additional print), and then search Summary of Optimized Potential Surface Scan in the output (it is at the very end). The values shown as Eigenvalues are the energies of the optimized geometry at each step. You can also identify which coordinate you have scanned, as all internal coordinates are provided. But this is way less friendly that using GaussView.

is it possible to scan the energy over two coordinates simultaneously?

Tahereh Alavi ,

Yes, there is no problem with that. You can add the second coordinate even with the GaussView menu (Tools → Redundant coordinates → Add). Once the calculation succesfully finishes, visualizing it with GaussView will generate a 3D plot, where the x and y axis are the scan coordinates and z is the energy

Thanks very much for your answer Enrique Manuel Arpa I tried scanning the PES over two coordinates (I'll paste the input file here) but I got a segmentation violation error I was wondering if you can help me to find the issue

%chk=CCl.chk

%mem=16GB

%NProcShared=4

#scan wb97xd/aug-cc-pvdz

MCSC neutral GS syn

2 1

C

H 1 B1

H 1 B2 2 A1

H 1 B3 3 A2 2 D1

O 1 B4 4 A3 3 D2

C 5 B5 1 A4 4 D3

O 6 B6 5 A5 1 D4

S 6 B7 7 A6 5 D5

Cl 8 B8 6 A7 7 D6

B1=1.0959

B2=1.09589

B3=1.09221

B4=1.43872

B5=1.33443

B6=1.19571

B7=1.80113 20 0.2

B8=2.04024 20 0.2

A1=35.09978

A2=110.87010

A3=105.16903

A4=115.28502

A5=127.32526

A6=127.09023

A7=99.98365

D1=-122.85308

D2=-118.88444

D3=179.99363

D4=0.00219

D5=-179.99539

D6=0.00249

B 6 8 S 20 0.2

B 8 9 S 20 0.2

Tahereh Alavi ,

I don't see anything intrinsically wrong in the input. When I copy-pasted it, the protons in the methyl group are placed in a very weird way (they are not distributed as a regular, tetrahedral -CH3). That is one of the reasons I am not a fan of rigid scans with z-matrices, relaxed scans with Cartesian coordinates are less risky in that way.

But it shouldn't affect the outcome. Could you please provide the output, or indicate exactly at which point of the calculation it crashes (if it is at the very beginning, after several steps...)?

Enrique Manuel Arpa thanks for your reply. It was a conversion failure issue, I increased the maxcycle (scf) and it terminated normally, I agree that a simple scan doesn't give us an accurate PES but it gave me a rough idea of how it'll change as these two bonds break (I guess!)

Tahereh Alavi ,

In those situations, I like to use the scf=xqc keyword. Basically, it initially tries the standard algorithm for SCF convergence, and if (and only if) it fails, it switches to the quadratically convergent procedure. So that keyword is completely free to use, because it does nothing unless the calculation fails