Mr. Nairat,

The mentioned in your posting values of states weights are "a priori" probabilities, which are different from "a posteriori" probabilities. Both Gaussian and Gamess use state-average CAS. This method provides the best description of a state, which however represents sum of states when a priori probabilities are 0.5 and 0.5 or when both the orbitals are treated equivalently during the optimization. In addition the number of adiabatic states are preliminary fixed as well.

These a priori probabilities and are input parameters running SA-CASSCF, and your computations are fixed within the framework of these parameters.

You have not optimization towards these values in order to obtain a posteriori probabilities. The convergence criteria determine this ensemble of states in a stationary/equilibrium conditions. So when you have treated in a different way the orbitals, there is a high probability of a convergence failure. In this context would be better if you have pay attention to ref. 1.

I can deduce that you would like, in fact, to search the conical intersection, introducing optimization. If yes, then you again operate with probabilities 0.5/0.5 by default. There is treated an equilibrium system, where the distribution of all the states is equal; in terms Boltzmann's distribution and the free Gibbs energy is a scalar at equilibrium. Presumably you should obtain equal energies for both the states of interest at the conical intersection point.

There is possible to operate not at equilibrium states by CASSCF, using RASSCF. Details about CASSCF, RASSCF and MCSCF can be found in ref. 1, as well.

By Gaussian you should add:

iop(5/13=1) (L502: Action on convergence failure, http://www.gaussian.com/g_tech/g09iop.htm)

In this case you have received a set energy values to given state, which are not converged along with the converged states. These are conditions out of equilibrium.

An example as output:

Davidson Disk Diagonalization is being used:

MVecIn= 3 MSekIn= 3 NVConv= 3 MConv= 3 WhenSc= 0 IFact= 1 FinIt=F

NRoot= 3 Conv= 0.00D+00 MaxIt= 4665 Dim= 4665 Max sub= 200

Iteration 1 internal iteration 1 dimension is 3

Eigenvalues at current iteration:

Root 1 : -153.158624110448200

Root 2 : -152.862406291481700

Root 3 : -152.788076202762000

DEMax= 0.00D+00

State 1 2nd order lowering -0.176092888214E-03 norm 1.25D-03

State 2 2nd order lowering -0.260172055668E-03 norm 1.52D-03

State 3 2nd order lowering -0.372762479031E-03 norm 1.82D-03

Iteration 2 internal iteration 2 dimension is 6

Root 1 not converged, maximum delta is 0.012779504037237

Root 2 not converged, maximum delta is 0.018518217598507

Root 3 not converged, maximum delta is 0.026171638967021

Eigenvalues at current iteration:

Root 1 : -153.170077193306800 Change is -0.011453082858651

...............................................

State 3 2nd order lowering -0.334135431834E-10 norm 5.44D-07

Iteration 15 internal iteration 15 dimension is 44

Root 1 has converged.

Root 2 not converged, maximum delta is 0.000007262106266

Root 3 not converged, maximum delta is 0.000027577835841

Eigenvalues at current iteration:

Root 1 : -153.170711905692900 Change is -0.000000000002473

Root 2 : -152.882972740971700 Change is -0.000000000619679

Root 3 : -152.817770599536500 Change is -0.000000002530072

DEMax= 2.53D-09

State 1 2nd order lowering -0.949671893843E-13 norm 2.89D-08

State 2 2nd order lowering -0.211847714938E-11 norm 1.37D-07

State 3 2nd order lowering -0.134747754909E-10 norm 3.45D-07

Iteration 16 internal iteration 16 dimension is 46

Root 1 has converged.

Root 2 not converged, maximum delta is 0.000003194008866

Root 3 not converged, maximum delta is 0.000019783691404

Eigenvalues at current iteration:

...........................

But if you fix only the probabilities in your input parameters within the frame of CASSCF computations, depending on the number of orbitals/electrons you can obtain set a posteriori probabilities, different from your starting probabilities (0.5, 0.5 for example). Please pay attention to a computation:

H2O, CASSCF, root = 2, Starting w.: 0.5, 0.5:

IBUJAK length= 840

Integral file not found: evaluate integrals

Symmetry not used in FoFDir.

MinBra= 0 MaxBra= 1 Meth= 1.

IRaf= 0 NMat= 1 IRICut= 1 DoRegI=T DoRafI=F ISym2E= 0 JSym2E=0.

Defining IBUGAM

State Average Calculation. The weights are:

St.: 1 w.=******** # St.: 2 w.=0.000000 # St.:

2ND ORD PT ENERGY CV 0.000000 CU -12.177303 UV -0.041710

TOTAL -87.940139

WARNING! : large rotation I J = 13 2

Step scaled by 0.5757176152963087

ITN= 1 MaxIt= 64 E= -75.7211263136 DE=-7.57D+01 Acc= 1.00D-08 Lan= 0

ITN= 2 MaxIt= 64 E= -68.8183467109 DE= 6.90D+00 Acc= 1.00D-08 Lan= 0

WARNING! : large rotation I J = 13 2

Step scaled by 0.4172629602417476

ITN= 3 MaxIt= 64 E= -67.4823592914 DE= 1.34D+00 Acc= 1.00D-08 Lan= 0

ITN= 4 MaxIt= 64 E= -70.5169032556 DE=-3.03D+00 Acc= 1.00D-08 Lan= 0

..............................................

ITN= 8 MaxIt= 64 E= -72.0108600441 DE= 4.04D-02 Acc= 1.00D-08 Lan= 0

ITN= 9 MaxIt= 64 E= -71.9699469402 DE= 4.09D-02 Acc= 1.00D-08 Lan= 0

ITN= 10 MaxIt= 64 E= -71.9271655222 DE= 4.28D-02 Acc= 1.00D-08 Lan=

...........

XXXZ= 0.0000 YYYX= 0.0000 YYYZ= 0.0000 ZZZX= 0.0000

ZZZY= 0.0000 XXYY= 0.0318 XXZZ= -0.8341 YYZZ= -0.8356

XXYZ= 0.0000 YYXZ= 0.0000 ZZXY= 0.0000

N-N= 9.343865761660D+00 E-N=-1.012543056861D+02 KE= 3.864258625820D+01

Orbital energies and kinetic energies (alpha):

1 2

1 O 0.000000 21.097084

2 O 0.000000 4.214473

3 O 0.000000 4.617222

......................................................

10 V 0.000000 2.196784

11 V 0.000000 0.937641

.....................................................

13 V 0.000000 2.559239

14 V 0.000000 117.539358

Total kinetic energy from orbitals= 7.027999075653D+01

..........................................

As you could expect the shown a posteriori probabilities have not physical meaning (w.=******** (state 1), while w.=0.000000 (state 2)). And the obtained energy (via iop(5/13=1)) does not correspond to the equilibrium state. It corresponds to an unknown state.

Ref. 1. C. Cramer, Essential of Computational Chemistry, 2nd Ed., Wiley, 2004, West Sussex, pp. 1-596.

In my opinion, your system can't control the weights(during the optimization) because automatically it is done from program's default (there is no a scientific reason for that). As you know by default, the weights of all states are identical, which is normally the most sensible choice. The symmetric situation during the optimization is also considerable .

In Molpro  it is  possible to optimize states of different symmetries together. In this case several wf / state directives can follow each other. I think in Gaussian09 it is possible too.Other thing averaging of states with different spin-multiplicity might be important which is possible only for CASSCF(not for restricted RASSCF or MCSCF, wavefunctions).

Although I don't know what is the condition of your system in the Gaussian calculation, I am sure in some conditions optimization of the system follows the default .for example : In the Gaussian  State Average is not allowed in combination with Opt=Conical or CASSCF = SpinOrbit, both of which perform state-averaged calculations by default.

Dear Dr. Ivanova,

Thanks a lot for your insights regarding states averaging and for providing such detailed and fruitful information, so much appreciated

Dear Dr. Monajjemi,

This was my first impression, system size (~107000 configurations), I just wanted to make sure about so from experts, thanks a lot

An overview about my problem, am trying to find a minima on the S2 surface for a system with a multi-reference character

I was planing to start with the most sensible situation, which is equal states averaging for a restricted closed shell system, hence spin is conserved, with 12 electrons in 11 orbitals for the first three states, the active space is carefully selected

single point cassscf converges in less than 20 steps with good diagonal density matrix values.

input:

#P CASSCF(12,11,NRoot=3,stateaverage)/6-31g(d) Guess=read Geom=allcheck

SCF(Conver=6,MaxCycle=512)

0.33 0.33 0.34

part of the output:

....

State Average Calculation. The weights are:

St.: 1 w.=0.330000 # St.: 2 w.=0.330000 # St.: 3 w.=0.340000 # St.:

....

1 2 3 4 5

1 0.191198D+01

2 0.691703D-02 0.155991D+01

3 -0.417758D-02 -0.335454D-01 0.187891D+01

4 -0.664441D-02 -0.255851D-01 0.129494D-01 0.185117D+01

5 -0.227274D-01 -0.557250D-01 0.885035D-01 0.102661D+00 0.104519D+01

6 0.100090D-01 0.215347D-01 -0.272887D-01 -0.434592D-01 -0.177645D+00

7 0.397191D-02 -0.849646D-01 -0.284044D-01 -0.557946D-01 0.476164D-01

8 0.226621D-02 -0.717392D-02 0.398236D-01 -0.163246D-01 -0.650156D-01

9 -0.100221D-01 0.241527D+00 0.407830D-01 -0.953562D-01 0.331086D+00

10 0.133762D+00 0.699834D-01 0.802615D-02 0.376610D-01 0.419056D-01

11 0.468519D-01 -0.323698D+00 -0.261217D-02 -0.109366D+00 -0.204965D+00

6 7 8 9 10

6 0.195374D+01

7 -0.183872D-01 0.149375D+00

8 0.936240D-02 -0.193215D-02 0.114312D+00

9 -0.466928D-01 0.210082D-01 0.269241D-02 0.398234D+00

10 -0.577987D-01 0.924515D-03 0.135055D-02 -0.427468D-01 0.170649D+00

11 0.968450D-01 -0.464819D-01 0.130495D-01 0.153042D-01 -0.495661D+00

11

11 0.966522D+00

MCSCF converged.

While when doing geometry optimization, starting from the converged casscf orbitals.,with same state averaging, output is:

State Average Calculation. The weights are:

St.: 1 w.=0.000000 # St.: 2 w.=0.500000 # St.: 3 w.=0.500000 # St.:

then it takes a lot of iterations to converge before changing geometry, and goes to structures that I know are not what am looking for.

Dear Gadzhiev,

Thats what I have been trying to say , I will clarify again for you

This input:

#P CASSCF(12,11,NRoot=3,stateaverage)/6-31g(d) Guess=read Geom=allcheck

SCF(Conver=6,MaxCycle=512)

0.33 0.33 0.34

gives this, which is right:

St.: 1 w.=0.330000 # St.: 2 w.=0.330000 # St.: 3 w.=0.340000 # St.:

while this input, which is geometry optimization for the previous run:

#P OPT CASSCF(12,11,NRoot=3,stateaverage)/6-31g(d) Guess=read Geom=allcheck

SCF(Conver=6,MaxCycle=512)

0.33 0.33 0.34

gives this:

St.: 1 w.=0.000000 # St.: 2 w.=0.500000 # St.: 3 w.=0.500000 # St.:

And that is the reason I am asking this question

Dear Muath

Please look at to the  rout section(Used the excited state of allyl in the C2v

symmetry ) as follows:

#P cas(3,3,nroot=2)/6-31g  test opt=(conical,z-matrix) nosymm

Optimization with state averaged gradients, and

appriximated state averaged second derivatives.

No direct CP-MCSCF.

It works completely

0   2

C

H,1,ch1

C,1,cc,2,cch

C,1,cc,2,cch,3,180.,0

H,3,ch2,1,hcc1,2,180.,0

H,3,ch3,1,hcc2,2,0.,0

H,4,ch2,1,hcc1,2,180.,0

H,4,ch3,1,hcc2,2,0.,0

Variables:

cc=1.46752553

ch1=1.07337584

ch2=1.07255071

ch3=1.07143914

cch=118.19703547

hcc1=120.92526605

hcc2=120.54512373

What do you think?