[molpro-user] MRCI symmetry breaking

[Jacky LIEVIN]

Dear Toru and interested users ,

I've tested two ways to work around the MRCI symmetry problem.

1) the first one is to adopt at all points the Cs configuration list (17310 CSFs with the VDZ basis set) by
 using the mrci directive  :  thresh,THRNRM=1.0d-10  
The number of pair functions is 49
Note that it is not necessary to change parameter THRDLP as in done previous message.
This works for all points but the exact D3h geometry (alpha=90.0). The solution is to calculate this point at a slightly displaced geometry (aplha=90.001). This is not elegant but it works...
90.001	-39.71754708   (Davidson corrected MRCI energies)
90.02	-39.71754708
90.04	-39.71754707
90.06	-39.71754705
90.08	-39.71754702
90.1	-39.71754699

2) the second one is to enforce mrci to use the C2v configuration list (17183 CSFs) using thresh,THRNRM=1.0d-6
The number of pair functions is now lowered to 48
90		-39.71754590   (Davidson corrected MRCI energies)
90.02	-39.71754589
90.04	-39.71754588
90.06	-39.71754586
90.08	-39.71754584
90.1	-39.71754580

These potential curves are parallel with an energy shift of about 1.E-6 hartrees.

Application of both solutions to the frequency calculation now provides correct out of plane frequencies: 494.40 and 495.10 cm-1 for solution 1 and 2 respectively with AVTZ.
These numbers are to be compared to 208 cm-1 obtained using the standard procedure. The latter result is meaningless because it uses numerical differences between points of both series.
Changing the basis set or the step size make the value change drastically. The frequency is imaginary in some cases.

In conclusion:
1) I do prefer solution 1 because it takes additional pair functions into account and they may gain importance at larger displacements from the higher symmetry point.
2) I still don't understand the real origin of the discontinuity which suggests that something is missing at the D3h(C2v) symmetry point. I'm wondering if it would be possible to enforce molpro to generate the configuration list of the lowest symmetry group (here Cs) even at points of higher symmetry (here C2v). It would help calculating smooth potential surfaces
3) I guess that similar problems should occur for many other molecules, probably those belonging to non abelian groups. I thus suggest users to carefully check what happens at the vicinity of a symmetry change before calculating frequencies. In the present case (CH3) the concerned frequency is low and errors in the numerical differentiation are important enough to realize that the frequency calculation fails. I'm afraid that the diagnostic can be less evident in other cases.
Any comment or suggestion are welcomed

Jacky

Le 25 mars 2011 à 11:26, Jacky LIEVIN a écrit :

> Dear Toru,
> 
> many thanks for your help.
> The additional thresh directive actually eliminates the jump at 90.06, but it still occurs at 90.0:
> 
>    ANGLE       ECASSCF        EMRCI          EMRCIDAV
>    90.00   -39.61901934   -39.71401806   -39.71754590
>    90.02   -39.61901934   -39.71401906   -39.71754708
>    90.04   -39.61901933   -39.71401905   -39.71754707
>    90.06   -39.61901932   -39.71401903   -39.71754705
>    90.08   -39.61901930   -39.71401900   -39.71754702
>    90.10   -39.61901928   -39.71401897   -39.71754699
> 
> Any other suggestion?
> According to your comment, I've now calculated all points with Cs symmetry, using the following input:
> 
> *** CH3 fondamental
> gthresh,energy=1.d-8,step=1.d-5
> gprint,orbitals
> basis=vdz
> rx=1.0 ang
> r=1.09503    ang
> hx=90.0
> dhh=0.02
> hh=hx
> d=120
> symmetry,y
> geometry={C;H2 C r;X1 C rx H2 hh;H3 C r X1 hh H2 d;H4 C r X1 hh H2 -d}
> {hf
> wf,9,1,1}
> hh=hh-dhh
> do i=1,6
> hh=hh+dhh
> symmetry,y
> geometry={C;H2 C r;X1 C rx H2 hh;H3 C r X1 hh H2 d;H4 C r X1 hh H2 -d}
> {multi
> closed,1
> wf,9,1,1
> print,civector}
> angle(i)=hh
> ecasscf(i)=energy(1)
> {ci
> wf,9,1,1
> thresh,THRDLP=1.0d-10,THRNRM=1.0d-10}
> emrcidav(i)=energd(1)
> emrci(i)=energy(1)
> enddo
> table,angle,ecasscf,emrci,emrcidav
> ---
> 
> 
> Le 25 mars 2011 à 09:50, Toru Shiozaki a écrit :
> 
>> Dear Jacky,
>> the jump in 90.06 is related to the pair elimination when orthogonal pair functions are generated.
>> Could you try the tighter thresholds for thrdlp and thrnrm as follows (or tighter value if you want):
>> {mrci; thresh,THRDLP=1.0d-10,THRNRM=1.0d-10}
>> 
>> In addition, it is not a good practice in general to use MRCI calculations with different symmetry (C2v and Cs in your case)
>> to construct the PES, since MRCI may not be invariant (but only slightly) with respect to the choice of the symmetry.
>> This is due to how we generate the CI basis functions.
>> 
>> Best wishes,
>> Toru
>> 
>> On Mar 24, 2011, at 9:30 PM, Jacky LIEVIN wrote:
>> 
>>> dear users,
>>> 
>>> The input below calculates few points of the potential curve for the rigid umbrella motion of CH3+ in its ground electronic state (2 A2" (D3h) --> 2 B2 (C2v) --> 2 A' (Cs)).
>>> The umbrella angle is changed from 90.0 to 90.1 degrees by steps of 0.02 degrees and the results are:
>>> 
>>> ANGLE         ECASSCF          EMRCI             EMRCIDAV
>>> 90.00      -39.61901934   -39.71401806   -39.71754590
>>> 90.02      -39.61901934   -39.71401805   -39.71754589
>>> 90.04      -39.61901933   -39.71401804   -39.71754588
>>> 90.06      -39.61901932   -39.71401903   -39.71754705
>>> 90.08      -39.61901930   -39.71401900   -39.71754702
>>> 90.10      -39.61901928   -39.71401897   -39.71754699
>>> 
>>> The first point (D3h geometry) is calculated in C2v symmetry while the others (C3v) are calculated in Cs.
>>> A sudden stabilization (below the planar equilibrium energy) occurs in the MRCI and MRCI+Davidson curves at 90.06 degrees . Note that the CASSCF curve doesn't exhibit such a feature.
>>> Comparison of MRCI calculations below and above 90.06 shows that the number of N-2 electron functions increase from 48 to 49 and that the number of contracted configurations moves accordingly from 17183 to 17310.
>>> Would it be possible to control these numbers in order to obtain smooth curves correlating C3v point to planar D3h?
>>> 
>>> Note that this problem occurs in the frequency calculation with as a consequence wrong frequencies (too low) for the out of plane mode. 
>>> 
>>> Any ideas about the origin of this symmetry breaking and how to solve it are welcomed?
>>> 
>>> Many thanks in advance
>>> 
>>> Jacky
>>> 
>>> 
>>> 
>>> *** CH3 fondamental
>>> gthresh,energy=1.d-8,step=1.d-5
>>> gprint,orbitals,civector
>>> basis=vdz
>>> rx=1.0 ang
>>> r=1.09503    ang
>>> hx=90.0
>>> dhh=0.02
>>> hh=hx
>>> d=120
>>> geometry={C;H2 C r;X1 C rx H2 hh;H3 C r X1 hh H2 d;H4 C r X1 hh H2 -d}
>>> angle(1)=hh
>>> {hf
>>> wf,9,3,1}
>>> {multi
>>> wf,9,3,1}
>>> ecasscf(1)=energy(1)
>>> {ci
>>> maxiter,50
>>> wf,9,3,1}
>>> emrci(1)=energy(1)
>>> emrcidav(1)=energd(1)
>>> do i=2,6
>>> hh=hh+dhh
>>> symmetry,nosym
>>> geometry={C;H2 C r;X1 C rx H2 hh;H3 C r X1 hh H2 d;H4 C r X1 hh H2 -d}
>>> hf
>>> {multi
>>> closed,1
>>> start,2140.2
>>> wf,9,1,1}
>>> angle(i)=hh
>>> ecasscf(i)=energy(1)
>>> {ci
>>> wf,9,1,1}
>>> emrcidav(i)=energd(1)
>>> emrci(i)=energy(1)
>>> enddo
>>> table,angle,ecasscf,emrci,emrcidav
>>> ---
>>> 
>>> 
>>> _____________________________________
>>> Jacky Liévin
>>> Service de Chimie Quantique et Photophysique
>>> Université Libre de Bruxelles, CPi 160/09
>>> 50 av F.D. Roosevelt
>>> B-1050 Bruxelles
>>> Belgium
>>> Tel: +32-2-650 4089
>>> Fax: +32-2-650 4232
>>> _____________________________________
>>> 
>>> _______________________________________________
>>> Molpro-user mailing list
>>> Molpro-user at molpro.net
>>> http://www.molpro.net/mailman/listinfo/molpro-user
>> 
> 
> 
> _____________________________________
> Jacky Liévin
> Service de Chimie Quantique et Photophysique
> Université Libre de Bruxelles, CPi 160/09
> 50 av F.D. Roosevelt
> B-1050 Bruxelles
> Belgium
> Tel: +32-2-650 4089
> Fax: +32-2-650 4232
> _____________________________________
> 
> _______________________________________________
> Molpro-user mailing list
> Molpro-user at molpro.net
> http://www.molpro.net/mailman/listinfo/molpro-user


_____________________________________
Jacky Liévin
Service de Chimie Quantique et Photophysique
Université Libre de Bruxelles, CPi 160/09
50 av F.D. Roosevelt
B-1050 Bruxelles
Belgium
Tel: +32-2-650 4089
Fax: +32-2-650 4232
_____________________________________