Resgrp:comp-photo-CHD
Overview
In the following sections it will be described how one can calculate the conical intersection seam of CHD (cyclo-hexa-1,3-diene) respectively cZc-HT (CZc-hexa-1,3,5-triene) with the new algorithm. Therefore the
- Optimization of a minimum point on the S0/S1 crossing seam
- Computation of Intersection-Space Frequencies
- Optimization of a Saddle point on the S0/S1 CHD crossing seam
- Computation of the Intersection Space Minimum Energy Path (MEP)
will be demonstrated.
Development version and Substitutions
For the calculations you have to use the Gaussian Development version, for the use of the current stable version your jobscript must contain following lines:
module load gaussian/devel-modules
module load gdvg01_08
Also the beginning of your Input-file (.com) has to contain the following lines:
%subst l103 /home/gaussian-devel/gaussiandvg01_pgi_lindafix/gdv_fabrizio/l103
%subst l115 /home/gaussian-devel/gaussiandvg01_pgi_lindafix/gdv_fabrizio/l115
%subst l510 /home/gaussian-devel/gaussiandvg01_pgi_lindafix/gdv_fabrizio/l510
%subst l716 /home/gaussian-devel/gaussiandvg01_pgi_lindafix/gdv_fabrizio/l716
%subst l1003 /home/gaussian-devel/gaussiandvg01_pgi_lindafix/gdv_fabrizio/l1003
Optimization of a minimum point on the S0/S1 crossing seam
For the optimization of a minimum point with the new algorithm you can use a standard input file. When you want to optimize a conical intersection geometry between the q-th and the (q-1)-th electronic state with an active space of n electrons in m orbitals the input file must contain the following lines:
#P CAS(m,n,nr = q)/[basis set]
Opt = Conical
IOp(1/5=-1,10/10=500005)
[Input Title]
[Charge] [Spin Multiplicity]
[Geometry]
The IOp(1/5) set to minus one specifies that we are attempting to optimize a minimum conical intersection point along the crossing seam. The IOp(10/10) is set to 500005 in order to compute the average gradient, which is used in the new algorithm, by solving the SA-CPMCSCF equations. The following files are the input and output files for the optimization of the 2A1/1A1 conical intersection minimum of CHD:
Here also the IOp(1/16) is set to 5 in order to set the maximum allowable magnitude of the the eigenvalues of the second derative matrix. In order to check if this point is a saddle point or a minimum point a frequency calculation is necessary, which is described in the next section.
Computation of Intersection-Space Frequencies
For the calculation of the Intersection-Space Frequencies no keywords exist to run such calculation yet. Therefore, a non-standard root has to be composed. A non-standard root is a complete sequence of links with associated options specified in an input file and read in by Gaussian. It is important to note that the program executes each link in the specified order, therefore extreme caution has to be taken when composing a non-standard route. Thus, a non-standard route is generally created by modifying the sequence of links produced by the standard route most similar to the one desired. This is accomplished most easily with the testrt utility, for more information see here. In computing intersection-space frequencies, testrt is used to produce the route for a standard state-average CASSCF frequency calculation:
#P CAS(n,m,stateaveraged,nroot = q)/[basis set]
Nosymm
Freq
The route obtained with testrt gives the correct sequence of links in a SA-CASSCF frequency analysis with the corresponding default options. However, in an intersection-space frequency analysis one needs the Hessian of both the intersecting electronic states and the Hessian of the rotated optimized wavefunctions. Create the route using testrt and repeat the sequence of links 5, 8, 11, 10, 6, 7 another two times within the new route. Then add/edit the specific options as below. Thus, the non-standard route used for an intersection-space frequency analysis should have a form similar to:
1/.../1,3;
2/.../2;
3/.../1,2,3;
4/.../1,5;
5/..,17=41000200,.../10;
8/.../1;
11/.../1;
10/...,10=600007,.../3;
6/.../1;
7/...,45=52/1,2,3,16;
5/...,17=41000200,...,97=100/10;
8/.../1;
11/.../1;
10/...,10=600007,...,97=100/3;
6/.../1;
7/...,45=51/1,2,3,16;
5/...,17=41000200,...,79=1/10;
8/.../1;
11/.../1;
10/...,10=600007,...,79=1/3;
6/.../1;
7/...,45=53/1,2,3,16;
1/.../3;
99//99;
The following files are the input and output files for the frequency calculation of the optimized conical intersection minimum of CHD, computed in the previous section:
With the IOp(1/29) and the IOp(4/5) set to one, the geometry and the initial guess is read from the checkpoint file, which you got from the optimization. The IOp(4/17) and IOp(4/18) are set to six, which specifies the number of electrons and orbitals in the active space, one might need to add more keywords depending on the specific study and system.
So how you can see in the output file there is no negative frequency, therefore the structure of section 2 should be a conical intersection minimum.
Viewing the normal modes and frequencies of the seam (or the normal modes and frequencies of the energy difference or of the bilinear coupling)
The scripts dffx_modes_g01.sh, seam_modes_g01.sh and eta_modes_g01.sh can be used to create log files that gaussview can read to visualise the normal modes and frequencies of the energy difference, the modes of the seam and/or the modes of the bilinear coupling respectively.
First, you need a directory (on your mac or on the cluster) that we will call /dir/. In here, you put the the scripts (dffx_modes_g01.sh, seam_modes_g01.sh, eta_modes_g01.sh) and a file called template.log, and you create a subdirectory /dir/gaussian_output/ in which you put your log file (say it's called jobname.log). Then you create a subdirectory /dir/data/. Then you execute the script this way:
./dffx_modes_g01.sh jobname
and you will get a new file in data that is called jobname_DFFX.log (for energy difference frequencies) or jobname_GV.log (for seam modes) and that can be read by GaussView.
N.B.: ./ might not be needed on the cluster (also this may depend on the type of shell), but it is on my Mac.
Now, there used to be a small printing error in Fabrizio's code that may not have been fixed. It's about the mass-weighting of the branching-space vectors, gradient difference and derivative coupling (the last two modes, given with zero frequencies). The consequence is that the motion of H atoms is exaggerated when you visualise these last two vectors with Gaussview. The other ones are fine.
NB the convention for the energy difference is S0 - S1, so if you want the modes that reduce the energy difference, take those that have a large positive frequency (with respect to the energy difference).
Optimization of a Saddle point on the S0/S1 CHD crossing seam
Generally, in order to optimize a saddle point structure, one must start from a geometry that is close enough to the quadratic region of such a point. In addition, the Hessian at this initial structure must have a negative eigenvalue, and the corresponding eigevector must be a suitable approximation for the direction connecting the two valleys. When one tries to optimize a saddle point within the intersection space, the initial geometry should also be a point of degeneracy (or near-degeneracy), since the intersection-space Hessian must be computed. Thus it is sometimes useful first to make a single point calculation and then a energy gap optimization between the two states along the gradient difference vector. Input and output files of the single point calculation:
Input and output files of the energy gap optimization:
Then it is also quite useful first to make a frequency calculation of this point:
And then one can do with this checkpoint file the optimization of the first-order saddle point:
The IOp(1/8) set to one specifies the maximum step size allowed during the optimization and the IOp(5/7) set to 250 specifies the maximum number of iterations. Because this time we want to get a first-order saddle point instead of a minimum also the IOp(1/5) is set minus 2. To prove that this is a first-order saddle point a frequency calculation must be done, which is also necessary for the Computation of the Intersection Space Minimum Energy Path (section 5):
Media:CI_trans_opt_freq_com.ogg
Media:CI_trans_opt_freq_log.ogg
The one negative frequency indicates that the point is really a first-order saddle point.
Computation of the Intersection Space Minimum Energy Path (MEP)
Currently no standard route exists to compute a minimum energy path limited to the intersection-space. Therefore, as discussed in section 3, a non-standard route has to be specified in the input file in order to perform these calculations. Using testrt and carefully modifying the output obtained for a standard intrinsic reaction path computation, which is specified by the IRC keyword, one should obtain something similar to:
1/...,44=-3,.../1,15;
2/.../2;
3/.../1,2,3;
4/.../1,5;
5/.../10;
8/.../1;
11/.../1;
10/10=500005,.../3;
6/.../1;
7/.../16;
1/...,44=-3,../15(2);
2/.../2;
99//99;
2/.../2;
3/.../1,2,3;
4/.../1,5;
5/.../10;
8/.../1;
11/.../1;
10/10=500005,.../3;
7/.../16;
1/...,44=-3,.../15(-8);
2/.../2;
6/.../1;
99/9=1/99;
In order to evaluate the intersection-space minimum energy path, there are added three more options to the existing IOp(1/44):
1/44 = -1 Compute the intersection-space minimum energy path in mass-weighted Cartesian coordinates.
1/44 = -2 Compute the intersection-space minimum energy path in internal coordinates.
1/44 = -3 Compute the intersection-space minimum energy path in mass-weighted internal coordinates.
For the following example the intersection-space minimum energy path (MEP) is computed in mass-weighted internal coordinates. For the computation of the MEP you need the transition point structure and the search vector from the optimization of the Saddle point from section 4. The transition point structure is the optimized Saddle point structure, while the search vector is the normal mode of the one negative frequency. But because in the output the normal mode is given in Cartesian coordinates, you first have to transform them in internal coordinates, by adding the normal mode to the optimized Saddle point structure, write out this structure in internal coordinates and calculate e.g with Exel the differences of the internal coordinates in comparison to the optimized Saddle point structure. These differences of the internal coordinates (search vector) have to be add to the transition point structure (in internal coordinates) in a 8F10.6 format in the input file:
So the calculated MEP leads from the transition point (CI_trans) to the minimum point (CI_min) on the S0/S1 crossing seam.