bizarre CCSD correlation energy

From NWChem

Viewed 2174 times, With a total of 15 Posts
Jump to: navigation, search

Clicked A Few Times
Threads 3
Posts 15
Dear Forum,

I have been trying to do some CCSD(T) calculations on some open-shell species and I am getting some correlation energies that don't seem to make a lot of sense.

So I have a (OH) substituted phenyl radical and I'm looking at the conformational energy difference between two OH rotomers. If I perform an ROCCSD(T) calculation on the two rotomers I get an energy difference:

c1-c2 = +9.6 kJ mol-1

with energies like the following:

c1:

CCSD correlation energy / hartree =        -1.138729684526594
CCSD(T) correction energy / hartree = -0.054336569695383
CCSD(T) total energy / hartree = -307.359516043510553

c2:

CCSD correlation energy / hartree =        -1.138739493778493
CCSD(T) correction energy / hartree = -0.054526511816299
CCSD(T) total energy / hartree = -307.363175100201033

However if I do this with UCCSD(T) then I get the following:

c1:

CCSD correlation energy / hartree =        -1.111323500626606
CCSD(T) correlation energy / hartree = -1.162727145664610
CCSD(T) total energy / hartree = -307.357498242822828

c2:

CCSD correlation energy / hartree =        -0.720195222066366
CCSD(T) correlation energy / hartree = -0.744704157414741
CCSD(T) total energy / hartree = -306.942326549647646

So you see that the correlation energy for the c2(UCCSD) calculation is about half as much as all the other calculations. How did this happen? Any solution here would be greatly appreciated. I have checked the following and have now run out of possible ideas:

(1) The SCF converges into the in/correct solution. This has been cross-checked against the same calculation with the a different program and was able to get the same SCFs in all 4 calculations.

(2) Ok so then it is something to do only with the correlation calculation. Is it the number of frozen orbitals? Well I have made sure it is the same number in each calculation...can't be this.

(3) Symmetry? This is were I think something is up. conf2 is in cs symmetry whereas conf1 is in c1 symmetry. I have tried: both with and without SYMMETRY CS in the geometry directive and I get the same answer. But because the ROCCSD(T) calculation is correct but the UCCSD(T) is not, is it that all the correlation energy is not being calculated? This is the input file that I used:

memory stack 5000 mb heap 1000 mb global 5000 mb noverify
PERMANENT_DIR /home/gpw501/tryptophan/oxidation/benzene_OH/bz_OH/
SCRATCH_DIR /scratch/gpw501/
GEOMETRY
SYMMETRY CS
C 0.99645638 -0.42648682 0.00000000
C -1.13763952 0.06248425 1.22330296
C -1.13763952 0.06248425 -1.22330296
C -1.83376503 0.20763768 0.00000000
C 0.19113955 -0.24170594 -1.24554896
C 0.19113955 -0.24170594 1.24554896
H 1.44466674 -1.43731034 0.00000000
H 0.73399717 -0.33302689 -2.18585706
H 0.73399717 -0.33302689 2.18585706
H -2.89396286 0.45223898 0.00000000
H -1.67575204 0.20666754 -2.15983510
H -1.67575204 0.20666754 2.15983510
O 2.14601254 0.42733523 0.00000000
H 1.80402493 1.33388138 0.00000000
END
BASIS large
* library 6-311++G(2d,2p)
END
BASIS small
* library 6-31G*
END
set "ao basis" small
SCF
maxiter 500
DOUBLET
UHF
vectors input atomic output small.movecs
THRESH 1e-7
END
task scf energy

set "ao basis" large
 SCF
maxiter 500
DOUBLET
UHF
vectors input project small small.movecs output large.mos
THRESH 1e-7
END
TCE
SCF
tilesize 15
io ga
maxiter 400
CCSD(T)
FREEZE 7
END
TASK TCE ENERGY

Forum Vet
Threads 10
Posts 1623
s^2 and UHF v ROHF energies
Did you check the value of S^2 for the two rotomers and the energy difference between ROHF and UHF for each one of the rotomers?
Edo

Clicked A Few Times
Threads 3
Posts 15
s^2 and UHF v ROHF energies
Hi Edo,

Thanks for the suggestion. I have looked at the SCFs both UHF and ROHF and their associated S^2 values they are all correct and reproducible in other quantum chemistry packages:

i.e.

UHF(c1-c2) = + 7.5 kJ mol-1
ROHF(c1-c2) = +9.1 kJ mol-1
S^2(c1) = 1.16
S^2(c2) = 1.18

but
UCCSD(T)[c1-c2] = -1090.0 kJmol-1
ROCCSD(T)[c1-c2] = +9.7 kJ mol-1

The S^2 values are high but this is why I was careful making sure that I have the correct SCF solution (viz previous posts).

So why is the correlation energy so different for c2 doesn't make much sense to me? I thought it has something to do with the symmetry?

Forum Vet
Threads 10
Posts 1623
I was asking for something different.
For c1 I would like to see E_UHF-E_ROHF
For c2 I would like to see E_UHF-E_ROHF

Clicked A Few Times
Threads 3
Posts 15
Hi Edo,

Sorry my mistake. These numbers are:

c1(E_UHF-E_ROHF) = 72 kJ mol
c2(E_UHF-E_ROHF) = 5 kJ mol

Clicked A Few Times
Threads 3
Posts 15
Opps
Sorry I made a typo in the last post:

The energies are:

c1(E_UHF-E_ROHF) = 72 kJ mol
c2(E_UHF-E_ROHF) = 74 kJ mol

Forum Vet
Threads 10
Posts 1623
What MP2?
Does it exhibit the same behavior?
You can find the MP2 energies just before the beginning of the CCSD cycles.

Clicked A Few Times
Threads 3
Posts 15
For some reason the log file does not have the MP2 energies (or I can't find them)..... this is just before the CCSD iterations (before this is just the SCF):

t1 file size   =             4778
t1 file name = /scratch/gpw501//bz_
t1 file handle = -998
T2-number-of-boxes 1050

t2 file size   =         19139799
t2 file name = /scratch/gpw501//bz_
t2 file handle = -995

CCSD iterations
-----------------------------------------------------------------
Iter Residuum Correlation Cpu Wall V2*C2
-----------------------------------------------------------------
1 0.1921757054645 -0.6598518759381 5711.1 6443.5 1092.3
2 0.0491557688756 -0.7025738833947 1399.0 1527.9 270.6
3 0.0209091708769 -0.7133177166611 2196.1 2296.5 440.9
4 0.0119476548267 -0.7165434872009 1056.3 1104.4 331.7

Clicked A Few Times
Threads 3
Posts 15
What string should I search on to find them?

Forum Vet
Threads 10
Posts 1623
Quote:Gpw501 Sep 25th 7:26 am
What string should I search on to find them?




MP2 Energy (coupled cluster initial guess)
------------------------------------------
Reference energy: -152.139429055398381
MP2 Corr. energy: -0.589012063944303
Total MP2 energy: -152.728441119342676

Clicked A Few Times
Threads 3
Posts 15
So there is no MP2 energy in my log files.

Forum Vet
Threads 10
Posts 1623
Sorry about it.
I was thinking of the old CCSD code.
I would suggest you to try to run the MP2 code on this molecules using the UHF wavefunction to have another point of comparison.

Clicked A Few Times
Threads 3
Posts 15
OK will try this...thanks.

Clicked A Few Times
Threads 3
Posts 15
Symmetry
Hi Edo,

I decided to look at the symmetry aspect of the calculation. The conformation with the seemingly small CCSD and CCSD(T) correlation energy is in CS symmetry. In the new code it automatically checks this symmetry with autosymm. If I use SYMMETRY CS then as I said before I get the same answer as with out this directive present:

CCSD correlation energy / hartree = -0.720195222066366


Just by chance I decided to rerun the calculation with SYMMETRY C1 as a directice and surprisingly I am now getting a correlation energy that matches up with the other 3 calculations. The job hasn't finished but you can see from the iterations that the correlation energy is much larger:

MICROCYCLE DIIS UPDATE:                  480                    5
481 0.0004320511988 -1.1118181107609 130.3 130.4 41.5
482 0.0008659036337 -1.1118394392735 129.3 129.8 41.7
483 0.0057894855964 -1.1118555457920 129.3 129.5 41.9
484 0.0510192426928 -1.1118676172230 130.0 130.3 42.0
485 0.4662714641661 -1.1118759942074 130.3 130.7 43.0
MICROCYCLE DIIS UPDATE: 485 5
486 0.0004166784843 -1.1118190111607 130.7 130.8 42.4
487 0.0004444610069 -1.1118401224515 130.6 132.2 41.5
488 0.0008756038419 -1.1118560832914 128.2 130.3 41.6
489 0.0056253402134 -1.1118681776394 129.4 131.0 42.0
Edited On 8:25:08 AM PDT - Thu, Sep 26th 2013 by Gpw501

  • Karol Forum:Admin, Forum:Mod, NWChemDeveloper, bureaucrat, sysop
    Profile
    Send PM
Clicked A Few Times
Threads 1
Posts 36
Hello,
it looks to me (I may be wrong) that this may be a case of switching between various UHF solutions (or HF instability). Have you tried to "continue" the UHF solution from c1 to c2 and then perform CCSD(T) calculations on the resutling UHF reference?

In order to do so, please move from c1 to c2 geometry using small "geometry increments" and for
next geometry use the converged UHF vectors from the previous geometry. To do so please use the
"vectors input" "vectors output" directives described in the HF section of the NWChem manual. You don't have to perfrom CCSD(T) calculations for each intermediate geometry. The point is to perform one-step CCSD(T) calculations for c2 geometry using UHF reference obtained in such a way.
Best regards,
Karol

Clicked A Few Times
Threads 3
Posts 15
c
c
Edited On 4:23:03 PM PDT - Sat, Sep 28th 2013 by Gpw501


Forum >> NWChem's corner >> Running NWChem



Who's here now Members 0 Guests 1 Bots/Crawler 0


AWC's: 2.5.10 MediaWiki - Stand Alone Forum Extension
Forum theme style by: AWC