If you have Molpro ab initio quantum chemistry package, you can use it instead of Gaussian. You need to set some environment variables that will be used by GRRM17 program as follows:
- bash) alias submol='molpro -W./ -d./'
- csh/tcsh) alias submol 'molpro -W./ -d./'
molpro is the Molpro execution command. -W./ and -d./ specify the scratch directory (see Molpro users manual for details). If the different Molpro execution command is used in your system, you need to set the alias submol appropriately.
The general input file format for calculations with Molpro 2006 is:
%link=molpro2006
# [Job type]
(blank line)
[Charge] [Spin multiplicity]
[Element and Cartesian coordinates of atom 1]
[Element and Cartesian coordinates of atom 2]
[Element and Cartesian coordinates of atom 3]
[Element and Cartesian coordinates of atom 4]
┆
Options
[Option 1]
[Option 2]
[Option 3]
[Option 4]
[Option 5]
Molpro Input
[A template of Molpro input, line 1]
[A template of Molpro input, line 2]
[A template of Molpro input, line 3]
┆
FORCE
---
XXXXX STATE n.1
END
[Option 6]
[Option 7]
┆
(blank line)
If you use Molpro 2006, you must declare "%link=molpro2006" in the first line of the input file. When a different version of Molpro is used, the first line should be modified appropriately, e.g. "%link=molpro2010" for Molpro 2010. Our tests were made only with versions 2006, 2008, 2010, and 2012.
Lines between the keywords "Molpro Input" and "---" will be copied to each Molpro input. The line between "---" and "END" specifies a keyword shown in the line presenting the total energy value in the Molpro punch file (e.g., the line should be MCSCF STATE 3.1 for the CASSCF energy for state 3.1). All Molpro calculations are performed with the nosymm option, so that the state symmetry should be n.1.
The following input is for a geometry optimization of formaldehyde in the first excited singlet (S1) state by the two-state CASPT2/6-311+G(d,p) method with an 8-electrons in 8-orbitals (8e, 8o) active space:
%link=molpro2006
# MIN
0 1
C -0.079213112255 0.000038936045 -0.592031587948
O 0.018800714554 -0.000018308130 0.717081381477
H 0.322739833191 0.909186505065 -0.947006485140
H 0.322739833191 -0.909359858691 -0.946723086445
Options
MaxStepSize=0.1
MO GUESS = mmm.wfu
MolMEM = 50
Molpro Input
BASIS=6-311+G(d,p)
{MCSCF,maxit=39,GRADIENT=1.d-8
{ITERATIONS;DO,UNCOUPLE,11,TO,39;}
START,*
ORBITAL,*
Occ,12;
Closed,4;
WF,16,1,0;
STATE,2;}
{RS2,Mix=2,Root=2,SHIFT=0.3,MaxIt=99,MaxIti=999;
ORBITAL,IGNORE_ERROR;
WF,16,1,0;
STATE,2;}
FORCE;
---
RSPT2 STATE 2.1
END
Use of MaxStepSize=x (x ≤ 0.1) is recommended, because large geometrical displacements often causes sudden changes in MOs. The dynamical memory (memory, 50, M) in Molpro calculations is set by the MolMEM=50 option. When the lines START,* and ORBITAL,* are copied to a Molpro input file, * are replaced by a proper MO file name. If these lines do not exist, MOs are re-calculated from a semi-empirical guess in every optimization steps.
It is recommended to prepare initial MOs beforehand, and the MO file name mmm.wfu should be given in the input file by the MO GUESS option. MOs can be prepared automatically at the initial geometry by using the following option:
Stable=Opt
BASIS=6-31G
{HF;
WF,16,1,2;}
{MCSCF,maxit=39,GRADIENT=1.0d-6
{ITERATIONS;DO,UNCOUPLE,11,TO,39;}
START,2100.2;
ORBITAL,*
Occ,12;
Closed,2;
WF,16,1,0;
STATE,4;}
END
With this option, initial MOs are prepared by an HF calculation for the lowest triplet state, and then obtained MOs are further optimized by a four-state-averaged SA4-(12e, 10o)-CASSCF/6-31G calculation. CASSCF with a large active space and with many averaged states usually gives a reasonable set of MOs.
To speed up the optimization, initial Hessian should be estimated by a low level method. Usage of the HESLOW option is little different from the case of G09. An example is:
HESLOW
BASIS=6-31G
{MCSCF,maxit=39,GRADIENT=1.0d-8
{ITERATIONS;DO,UNCOUPLE,11,TO,39;}
START,*
ORBITAL,*
Occ,12;
Closed,4;
WF,16,1,0;
STATE,2;
CPMCSCF,GRAD,2.1,SPIN=0.0,ACC=1.0d-8}
FORCE
---
MCSCF STATE 2.1
END
With this option, the initial Hessian is computed by the SA2-(8e, 8o)-CASSCF/6-31G method. Then, the Hessian will be updated by the 2S-(8e, 8o)-CASPT2/6-311+G(d,p) gradients in the above formaldehyde example. The SPHIGH option can also be used likewise.
Multi-surface calculations
The multistate options, i.e., OptX and ModelF, work also with Molpro. The format of the second state input is different from the one for G09. The following example is for an optimization of S1/S0 minimum energy conical intersection (MECI) point with the 2S-(8e, 8o)-CASPT2/6-311+G(d,p) method.
%link=molpro2006
# MIN
0 1
C -0.079213112255 0.000038936045 -0.592031587948
O 0.018800714554 -0.000018308130 0.717081381477
H 0.322739833191 0.909186505065 -0.947006485140
H 0.322739833191 -0.909359858691 -0.946723086445
Options
OptX(Conical,Guess=Same)
MaxStepSize=0.1
MO GUESS=mmm.wfu
MolMEM=50
Molpro Input
BASIS=6-311+G(d,p)
{MCSCF,maxit=39,GRADIENT=1.0d-8
{ITERATIONS;DO,UNCOUPLE,11,TO,39;}
START,*
ORBITAL,*
Occ,12;
Closed,4;
WF,16,1,0;
STATE,2;}
{RS2,Mix=2,Root=2,SHIFT=0.3,MaxIt=99,MaxIti=999;
ORBITAL,IGNORE_ERROR;
WF,16,1,0;
STATE,2;}
FORCE
---
RSPT2 STATE 2.1
@@@
BASIS=6-311+G(d,p)
{MCSCF,maxit=39,GRADIENT=1.0d-8
{ITERATIONS;DO,UNCOUPLE,11,TO,39;}
START,*
ORBITAL,*
Occ,12;
Closed,4;
WF,16,1,0;
STATE,2;}
{RS2,Mix=2,Root=1,SHIFT=0.3,MaxIt=99,MaxIti=999;
ORBITAL,IGNORE_ERROR;
WF,16,1,0;
STATE,2;}
FORCE
---
RSPT2 STATE 1.1
END
As seen in this example, the Molpro-input templates for the two states are given with a partition keyword @@@. In the Stable=Opt---END, HESLOW---END, and SPHIGH---END options, a template for the second state should be given likewise. When the MO GUESS option is used without the Guess=Same option (the default is Guess=Different), two MO files, i.e., mmm.wfu and mmm.wfu_1, are required for the first and second states, respectively.
Hints to avoid SCF convergence errors in the ADDF or AFIR searches
SCF convergence errors occur frequently during a global reaction route mapping. In the ADDF or AFIR searches, use of very small active space such as (2e, 2o) is strongly recommended. If the Stable=Opt---END option is used, a set of MOs is reconstructed in certain intervals of the search procedures by the method specified between Stable=Opt---END to obtain a correct set of active orbitals. When an orbital exchange (between inside and outside the small active space) happens, SCF does not converge unless the orbitals are not altered manually. This can be done automatically by modifying a Molpro input section (Molpro Input---END) as follows.
%link=molpro2006
# ADDF
0 1
C -0.079213112255 0.000038936045 -0.592031587948
O 0.018800714554 -0.000018308130 0.717081381477
H 0.322739833191 0.909186505065 -0.947006485140
H 0.322739833191 -0.909359858691 -0.946723086445
Options
NRUN=5
MaxStepSize=0.1
UpSize=5
DownSize=5
MolMEM = 50
Molpro Input
BASIS=6-31G(d)
{MCSCF,maxit=39,GRADIENT=1.0d-8
{ITERATIONS;DO,UNCOUPLE,11,TO,39;}
START,*
ORBITAL,*
Occ,9;
Closed,6;
WF,16,1,0;
STATE,1;
CPMCSCF,GRAD,1.1,ACC=1.0d-8;}
FORCE
---
MCSCF STATE 1.1
---
BASIS=6-31G
{HF;
WF,16,1,0;}
{MCSCF,maxit=39,GRADIENT=1.0d-8
{ITERATIONS;DO,UNCOUPLE,11,TO,39;}
START,2100.2;
ORBITAL,*
Occ,11;
Closed,5;
WF,16,1,0;
STATE,2;}
---
BASIS=6-31G
{HF;
WF,16,1,0;}
{MCSCF,maxit=39,GRADIENT=1.0d-8
{ITERATIONS;DO,UNCOUPLE,11,TO,39;}
START,2100.2;
ORBITAL,*
Occ,12;
Closed,4;
WF,16,1,0;
STATE,2;}
---
BASIS=6-31G
{HF;
WF,16,1,0;}
{MCSCF,maxit=39,GRADIENT=1.0d-8
{ITERATIONS;DO,UNCOUPLE,11,TO,39;}
START,2100.2;
ORBITAL,*
Occ,12;
Closed,2;
WF,16,1,0;
STATE,2;}
END
Stable=Opt
BASIS=6-31G
{HF;
WF,16,1,0;}
{MCSCF,maxit=39,GRADIENT=1.0d-6
{ITERATIONS;DO,UNCOUPLE,11,TO,39;}
START,2100.2;
ORBITAL,*
Occ,12;
Closed,3;
WF,16,1,0;
STATE,2;}
END
The above calculation uses SA1-(4e, 3o)-CASSCF/6-31G(d) method in the ADDF search. When a SCF convergence error happens, the program tries to obtain a set of correct active orbitals by a SA2-(6e, 6o)-CASSCF/6-31G single point calculation. Then, force is computed by SA1-(4e, 3o)-CASSCF/6-31G(d) using a set of MOs by SA2-(6e, 6o)-CASSCF/6-31G as an initial guess to continue the ADDF search. If this does not remove the convergence error, single point calculation(s) with computation level(s) specified below are further considered to obtain correct active orbitals. There is no limitation of the number of the trial computation levels. "---" is the keyword to partition each computation level. One can use this option also for the second state in multi-surface calculations as follows.
Molpro Input
BASIS=6-31G(d)
{MCSCF,maxit=39,GRADIENT=1.0d-8
{ITERATIONS;DO,UNCOUPLE,11,TO,39;}
START,*
ORBITAL,*
Occ,9;
Closed,6;
WF,16,1,0;
STATE,2;
CPMCSCF,GRAD,2.1,ACC=1.0d-8;}
FORCE
---
MCSCF STATE 2.1
---
BASIS=6-31G
{HF;
WF,16,1,0;}
{MCSCF,maxit=39,GRADIENT=1.0d-8
{ITERATIONS;DO,UNCOUPLE,11,TO,39;}
START,2100.2;
ORBITAL,*
Occ,12;
Closed,3;
WF,16,1,0;
STATE,3;}
@@@
BASIS=6-31G(d)
{MCSCF,maxit=39,GRADIENT=1.0d-8
{ITERATIONS;DO,UNCOUPLE,11,TO,39;}
START,*
ORBITAL,*
Occ,9;
Closed,6;
WF,16,1,0;
STATE,2;
CPMCSCF,GRAD,1.1,ACC=1.0d-8;}
FORCE
---
MCSCF STATE 1.1
---
BASIS=6-31G
{HF;
WF,16,1,0;}
{MCSCF,maxit=39,GRADIENT=1.0d-8
{ITERATIONS;DO,UNCOUPLE,11,TO,39;}
START,2100.2;
ORBITAL,*
Occ,12;
Closed,3;
WF,16,1,0;
STATE,3;}
END
Use of options limiting step sizes are also recommended; MaxStepSize=x (x ≤ 0.1), UpSize=n (n ≤ 5), and DownSize=n (n ≤ 5).