[molpro-user] DFT calculations produced different energy with a different guesses
Kirk Peterson
kipeters at wsu.edu
Thu Nov 4 18:00:34 GMT 2010
Dear David,
this turned out to be a very interesting problem. As I suspected there were issues over which orbitals were occupied in the HF and DFT calculations. The RHF seems to converge to a state where all three metals atoms have a singly occupied valence s orbital. This state has a dipole moment of nearly zero and is strongly dominated by this HF configuration. However after experimenting with some CASSCF calculations (which are actually quite small if the active space only includes the valence s and d orbitals of the metals, which are nearly all doubly occupied), I find that there are two states much lower in energy with relatively large dipole moments. These represent occupations whereby there is significant d character in the open-shell orbital on the central Cu atom. So I think in this case there is significant charge transfer from the Cu to the terminal Au atoms. These two ionic states are also quite multireference in character and are within just 3 mEh of each other. A UHF calculation actually converges to one of these states, but the RHF always goes to the strongly single reference case. Presumably this is also what causes you problems in your ROHF-CCSD(T) calculations. When I can get RHF-B3LYP to converge (an initial guess from a preceeding UHF seems to work well), it converges to an ionic state but it's probably a wash on which of the two you actually get, which explains the varying energies I believe. So I think a CASPT2 calculation would seem to be the minimum approach for this system, but a MRCI would be quite doable as well.
regards,
Kirk
On Nov 4, 2010, at 8:31 AM, David Danovich wrote:
> Dear Kirk,
>
> I attached log, out and com files for two cases (1 and 3 in my first post). You may see that for HF calculations orders of energies is Ok but when I used it as a guess for DFT I got not correct order.
> (I renamed com files as com.txt because mail agent does not allow to send com files).
> Let me know if you want me send files with ccsd(t) calculations.
>
> Thank you very much. David
>
> On Thu, Nov 4, 2010 at 5:05 PM, Kirk Peterson <kipeters at wsu.edu> wrote:
> Dear David,
>
> yes, please send me a representative output file.
>
> regards,
>
> Kirk
>
> On Nov 4, 2010, at 3:01 AM, David Danovich wrote:
>
>> Dear Kirk,
>>
>> At least it looks like it was converged correctly (no error or warning messages). The same order of energies also for the optimized structures. Do you want me to send input and output files to you?
>> I found this quite often (the order of the energies is not correct, even convergence is Ok) for the systems I am calculating (like 2CuAu, 2AuCu, 2Cu2Au in high spin state). May be for this reason I have troubles with CCSD(T) calculations? It took many cycles and often not converges.
>>
>> David
>>
>> On Wed, Nov 3, 2010 at 11:26 PM, Kirk Peterson <kipeters at wsu.edu> wrote:
>> David,
>>
>> you could try as well to do the HF calculation with the multi program to see if this has better convergence characteristics. I often do this for difficult cases. Just specify the same occ and closed as your final occupations below. I agree that it's strange your orbitals are not ordered by energy. Are you sure it converged correctly?
>>
>> -Kirk
>>
>> On Nov 3, 2010, at 1:45 PM, David Danovich wrote:
>>
>>> Dear Kirk,
>>>
>>> In general, you are right and it is always possible to converge to different electronic state, especially for the systems with transition metals. For DFT methods I always do TDDFT calculations in order to be sure that the state I get is a lower one (for such systems with 3 transition metal atoms I have no ideas how to do CASSCF with 33 valence electrons and I do not know a priory which d electrons are important). But for this particular case I think I have the same state with the same occupancy as you can see below and also with almost the same eigenvectors. But as I said before difference in total energy is 3.8 kcal/mol. What is also strange that for irreducible representation A1 Aufbau principle is not satisfy. You may see that orbital 11.1 has much lower energy than orbitals 6.1, 7.1, 8.1, 9.1, 10.1 for both calculations. My feeling is that something goes wrong with DFT calculations in Molpro. I will try to carry out this type of calculations with other programs.
>>>
>>> Thank you David
>>>
>>> Case 1)
>>>
>>> Final alpha occupancy: 12 5 9 4
>>> Final beta occupancy: 10 5 8 4
>>>
>>> Orbital energies:
>>>
>>> 1.1 2.1 3.1 4.1 5.1 6.1 7.1 8.1 9.1 10.1
>>> -4.671414 -4.118564 -3.065611 -2.328939 -2.323338 -0.438461 -0.379638 -0.343993 -0.313299 -0.288287
>>>
>>> 11.1 12.1 13.1 14.1
>>> -0.509678 -0.218071 -0.055982 -0.024797
>>>
>>> 1.2 2.2 3.2 4.2 5.2 6.2 7.2
>>> -3.028087 -2.319030 -0.478500 -0.324411 -0.306541 -0.066386 0.004508
>>>
>>> 1.3 2.3 3.3 4.3 5.3 6.3 7.3 8.3 9.3 10.3
>>> -4.118395 -3.028437 -2.326183 -2.321552 -0.481322 -0.301593 -0.288912 -0.282174 -0.196122 -0.041215
>>>
>>> 11.3
>>> 0.035134
>>>
>>> 1.4 2.4 3.4 4.4 5.4 6.4
>>> -2.318707 -0.440376 -0.298955 -0.284142 0.018513 0.139495
>>>
>>> Case 3)
>>>
>>> Final alpha occupancy: 12 5 9 4
>>> Final beta occupancy: 10 5 8 4
>>> 1.1 2.1 3.1 4.1 5.1 6.1 7.1 8.1 9.1 10.1
>>> -4.673062 -4.117760 -3.014878 -2.327371 -2.323105 -0.441349 -0.360882 -0.345206 -0.319350 -0.288610
>>>
>>> 11.1 12.1 13.1 14.1
>>> -0.514110 -0.216705 -0.051875 -0.024126
>>>
>>> 1.2 2.2 3.2 4.2 5.2 6.2 7.2
>>> -3.061378 -2.318440 -0.460149 -0.323223 -0.306348 -0.070494 0.001935
>>>
>>> 1.3 2.3 3.3 4.3 5.3 6.3 7.3 8.3 9.3 10.3
>>> -4.117591 -3.050251 -2.324891 -2.321034 -0.451631 -0.300020 -0.288446 -0.281090 -0.195381 -0.042282
>>>
>>> 11.3
>>> 0.035544
>>>
>>> 1.4 2.4 3.4 4.4 5.4 6.4
>>> -2.318118 -0.492930 -0.298691 -0.283211 0.018322 0.139423
>>>
>>>
>>>
>>>
>>> On Wed, Nov 3, 2010 at 9:40 PM, Kirk Peterson <kipeters at wsu.edu> wrote:
>>> David,
>>>
>>> it is certainly not always the case that the resulting HF or DFT solution is independent of the initial guess. In a perfect world that would be nice, but often a poorly constructed initial guess (or just bad luck) can result in convergence to a local minimum, i.e., an upper electronic state. This happens quite often for transition metal complexes because of the plethora of possible d-orbital occupations. The only way to get a feeling for what is going on is to print the orbitals and inspect the higher lying occupied ones. Is the bonding as you might expect? Another option is to carry out a CASSCF calculation as I suggested before in order to determine what the character of the higher lying occupied orbitals should really be. One can always then use the rotate command to swap in the correct orbitals in your HF or DFT calculation.
>>>
>>> -Kirk
>>>
>>> On Nov 3, 2010, at 9:03 AM, David Danovich wrote:
>>>
>>>> Dear Kirk,
>>>>
>>>> You are right, the calculations converge to different structures. But the question is why. In all calculations as you can see in the input I sent in my original post I used the same starting geometry, just guess wave functions for B3LYP 4B2 state were different. I got different energy already in the first step (at the point where geometry was the same for all three cases). It means that even for single pint calculations solution is depending on starting guess. In general, I can use any guess wave function and should get the same converged result. But as you see I got different results. Occupancy and eigenvalues are more or less similar for the cases 1) and 3) but in any case the difference in energy is around 3.8 kcal/mol. How one can be sure that he really get the lower solution (for particular geometry)? Do you have any suggestion.
>>>>
>>>> Thank you David
>>>>
>>>>
>>>> On Wed, Nov 3, 2010 at 5:01 PM, Kirk Peterson <kipeters at wsu.edu> wrote:
>>>> David,
>>>>
>>>> it looks to me that you're getting stuck into a few local minima. Do these different calculations converge to different structures? How different are the energies at the starting geometries? I would strongly recommend printing the orbitals and taking a careful look. You could also try doing a CASSCF calculation (print the orbitals and CI vector) to see how things should really be occupied.
>>>>
>>>> regards,
>>>>
>>>> Kirk
>>>>
>>>>
>>>> On Nov 3, 2010, at 6:19 AM, David Danovich wrote:
>>>>
>>>>> Hello,
>>>>>
>>>>> I am calculating CuAu2 molecule in high spin state (4B2 state) using B3LYP method. Below you can find three different possibilities I have used in the calculations. In the first one 1) I have calculated first HF wave function and used it as a guess for DFT calculation. Energy I got was -468.86467917 au. In the second calculation 2) I directly have calculated B3LYP without HF wave function. The energy I got was -468.83995914 au. In the third calculation 3) I first have done B3LYP calculation for state (4A2) and then used it as a guess for the calculations of 4B2. The energy I got was
>>>>> -468.85824088 au. As you can see there is quite large difference in energy for the same state. I was trying to use different grids but it does not solved the problem. In my opinion result should not depend on the guess so drastically. What can be a solution for this problem?
>>>>>
>>>>> Thank you in advance David
>>>>>
>>>>> 1)
>>>>>
>>>>> ***, Peterson PP tz basis set
>>>>> memory,250,m
>>>>> r = 2.72684681 ang;
>>>>> a = 60.65125012 degree;
>>>>> geometry={Cu1; !z-matrix geometry input
>>>>> Au2,Cu1,r;
>>>>> Au3,Cu1,r,Au2,a;
>>>>> }
>>>>> basis=cc-pVTZ-PP
>>>>>
>>>>> {hf,maxit=500;
>>>>> wf,57,3,3}
>>>>>
>>>>> {rks,b3lyp3,maxit=500;
>>>>> wf,57,3,3}
>>>>>
>>>>> optg
>>>>>
>>>>> Total energy -468.86467917 au
>>>>> ____________________________________________________
>>>>>
>>>>> 2)
>>>>> ***, Peterson PP tz basis set
>>>>> memory,250,m
>>>>> r = 2.72684681 ang;
>>>>> a = 60.65125012 degree;
>>>>> geometry={Cu1; !z-matrix geometry input
>>>>> Au2,Cu1,r;
>>>>> Au3,Cu1,r,Au2,a;
>>>>> }
>>>>> basis=cc-pVTZ-PP
>>>>>
>>>>> {ks,b3lyp3,maxit=500;
>>>>> shift,-0.1,0.0;
>>>>> wf,57,3,3}
>>>>>
>>>>> optg
>>>>>
>>>>> Total energy -468.83995914 au
>>>>>
>>>>> ______________________________________________________
>>>>>
>>>>> 3)
>>>>> ***, Peterson PP tz basis set
>>>>> memory,250,m
>>>>> r = 2.72684681 ang;
>>>>> a = 60.65125012 degree;
>>>>> geometry={Cu1; !z-matrix geometry input
>>>>> Au2,Cu1,r;
>>>>> Au3,Cu1,r,Au2,a;
>>>>> }
>>>>>
>>>>> ! cc-pVTZ-PP
>>>>>
>>>>> basis=cc-pVTZ-PP
>>>>>
>>>>> {ks,b3lyp3,maxit=500;
>>>>> wf,57,4,3}
>>>>>
>>>>> {ks,b3lyp3,maxit=500;
>>>>> shift,-1.0,0.0;
>>>>> wf,57,3,3}
>>>>>
>>>>> optg
>>>>>
>>>>> Total energy -468.85824088 au
>>>>>
>>>>> _______________________________________________
>>>>> Molpro-user mailing list
>>>>> Molpro-user at molpro.net
>>>>> http://www.molpro.net/mailman/listinfo/molpro-user
>>>>
>>>>
>>>
>>>
>>
>>
>
>
> <2aucu-b3lyp-3-com.txt><2aucu-b3lyp-3.log><2aucu-b3lyp-3.out><2aucu-b3lyp-3-2.log><2aucu-b3lyp-3-2.out><2aucu-b3lyp-3-2-com.txt>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.molpro.net/pipermail/molpro-user/attachments/20101104/093d6dd7/attachment.html>
More information about the Molpro-user
mailing list