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('');

177. Jerry-rigging g09 UV/VIS spectra in gnuplot and/or octave

EDIT: I had a nicer post with lots of figures before. Because I realised that the data is good enough to be included in a future paper we're working on, I had to take everything down again. All the data in the plots now is made up (hence 'fakeuv.dat'), and I haven't made the plots look nice.

I don't like proprietary formats for anything. They never, ever benefit anyone other than the software vendor.

Almost as bad as using binary proprietary formats is not providing export facilities to ascii formats.

I may have missed it, but I was using gaussview to look at td-dft calculated uv/vis spectra -- and couldn't find a way of exporting the data. Sure, I could export the graph as a png, svg etc. file. But not double column tab-separated ascii file.

There's a bit of fudging in what I'm doing  -- I'll be the first one to admit that.

So here's single line to export the wavelengths and intensities:
cat g03.g03out|grep Excited|grep -v singles|sed 's/=/\t/g'|gawk '{print $7,$10}'>uvvis.dat

You can plot them in gnuplot using
plot 'uvvis.dat' u 1:2 w impulse

The problem is that these are just spikes -- not the smooth uv/vis like spectra we're used to. On the other hand, if I understand things correctly, this is the REAL data, while the smoothed uv/vis spectrum above is more for presentation purposes. I might obviously be wrong, and I am by no stretch a computational or theoretical chemist - I just like their tools.

We've got an immensely powerful tool at our hands: Octave!

data=load('fakeuv.dat');
gauss= @(x,c,r,s) r.*1./(s.*sqrt(2*pi)).*exp(-0.5*((x-c)./s).^2)
x=linspace(250,850,600);
plot(x,cumsum(gauss(x,data(:,1),data(:,2),20)))

where 20 is an abitrary value. Anyway, this is how it looks:

We can try s=30 instead:

We export it

outdata=cumsum(gauss(x,data(:,1),data(:,2),30));
exportdata=[x' outdata'];
save 'uvvis2.sim' exportdata

and plot it in gnuplot
plot 'uvvis2.sim' u 1:48 w lines

It might not look like the UV/VIS spectrum you're used to, but as I said in the beginning, the data's all made up -- using 'real' calculated data I got a beautiful spectrum.