[molpro-user] MRCI symmetry breaking
Jacky LIEVIN
jlievin at ulb.ac.be
Thu Mar 31 07:33:21 BST 2011
Dear Ioannis,
thanks for your answer. I have already tried the "ref" directive, which is indeed made for adjusting the configuration space in case of symmetry lowering, as in the example you quote of a PI state. Unfortunately it doesn't help here because the state of interest is non degenerated (A"2 in D3h --> B2 in C2v).
yours
Jacky
Le 30 mars 2011 à 21:43, Ioannis Kerkines a écrit :
> Dear Jacky,
>
> Concerning your #2 conclusion: I think Molpro has a "ref" directive that you can use in the MRCI part of the input. If I remember correctly, you would use that in order to produce the same number of configurations between different irreps.
>
> I have used that in the past in order to "correctly" describe the cylindrical symmetry in diatomic molecules in MRCI calculations. For example for a Pi state, and using C2v representation (Pi --> B1 + B2), in CASSCF you could state average over B1 and B2, but in MRCI you could only do either B1 or B2 irrep. But with the ref command, e.g.: {ci;wf,9,2,1;ref,3} or {ci;wf,9,3,1;ref,2} you were certain to have exactly the same configuration space size in both cases, and therefore compare directly their energies.
>
> I don't know if you have used that already, or if that would be a solution to what you are asking. You can find the directive in the molpro manual.
>
> Best
> Ioannis
>
> On 30/3/2011 5:39 μμ, Jacky LIEVIN wrote:
>> 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
>> _____________________________________
>>
>> _______________________________________________
>> 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
_____________________________________
More information about the Molpro-user
mailing list