198. Nwchem -- freeze core and tddft on benzene

UPDATE: it's taken a while, but I've tested this on a large polyoxometalate cluster (>10 group 5/6 atoms and >30 oxygens) -- with a total of ca 700 alpha and beta orbitals, respectively, I froze 200 core orbs and used 3-21G/3-21G* w/ ub3lyp. The frozen core calculation took 44% of the time the full calculation took. The spectra look identical to within 3 nm (18 roots -- only 2 intermediate transitions have been shifted. All other transitions are identical). The full calc took 7 days and 4 hours, while the frozen calc took 3 days and 4 hours. 

Original post:
Benzene has 21 occupied α and 21 occupied β orbitals.

How many core orbitals can we freeze when looking at electronic transitions, how does freezing core orbitals affect the energies, and what are the performance gains, if any?

I used the following relevant statements

tddft
    cis
    freeze core X
    nroots 24
end
task tddft energy

Also, xc b3lyp and 6-31++g**.

Freeze core 10 means freeze 10 α and freeze 10 β orbitals.

Results:

Frozen Runtime    Transitions with f>0
0          1111 s     6.7835, 6.8038, 6.8042, 7.3479, 7.3496
5          1107 s     6.7835, 6.8038, 6.8042, 7.3480, 7.3498
10        1063 s     6.7838, 6.8040, 6.8045, 7.3761, 7.3862
15        1063 s     6.5878, 6.7840, 6.8042, 6.8046
20          719 s     6.8038, 6.8355, 7.1692, 7.5334, 8.1866

Evaluation:
We're not really looking at what's right or wrong -- the main goal is to understand how the results are affected (we might accidentally get 'right' answer doing something stupid). Freezing more than 10 α/β orbitals leads to significant differences in predicted transitions.

Taking oscillation strength into account and plotting it, we see that we really don't want to overdo it with the number of frozen cores -- 15 and 20 frozen cores yield results that are about as predictable as a coin toss. We also don't see any overwhelming performance gains, but that may well be due to the size and computational cost (or lack thereof) of the system.

Octave code:

spec1=load('bz631gppdd_cosmo.dat');

spec2=load('bz631gppdd_cosmo_f5.dat');
spec3=load('bz631gppdd_cosmo_f10.dat');
spec4=load('bz631gppdd_cosmo_f15.dat');
spec5=load('bz631gppdd_cosmo_f20.dat');

gau= @(x,c,w,i) i*(1/(w*sqrt(2*pi))).*exp(-0.5.*((x-c)./w).^2);
x=linspace(160,200,200);
y1=cumsum(gau(x,1241.25./spec1(:,2),4,spec1(:,3)));
y2=cumsum(gau(x,1241.25./spec2(:,2),4,spec2(:,3)));
y3=cumsum(gau(x,1241.25./spec3(:,2),4,spec3(:,3)));
y4=cumsum(gau(x,1241.25./spec4(:,2),4,spec4(:,3)));
y5=cumsum(gau(x,1241.25./spec5(:,2),4,spec5(:,3)));


subplot(3,2,1)
 plot(x,y1(rows(y1),:))
 axis([160 200 0 0.2]);
 title('0 frozen');
 subplot(3,2,3)
 plot(x,y2(rows(y2),:))
 axis([160 200 0 0.2]);
 title('5 frozen');
 subplot(3,2,5)
 plot(x,y3(rows(y3),:))
 axis([160 200 0 0.2]);
 title('10 frozen');
 subplot(3,2,2)
 plot(x,y4(rows(y4),:))
 axis([160 200 0 0.2]);
 title('15 frozen');
 subplot(3,2,4)
 plot(x,y5(rows(y5),:))
 axis([160 200 0 0.2]);
 title('20 frozen');
 subplot(3,2,6)
  plot(x,y1(rows(y1),:),x,y2(rows(y2),:), x, y3(rows(y3),:), x,y4(rows(y4),:),x,y5(rows(y5),:))
  axis([160 200 0 0.2]);
 title('');