New decomposition code

@Decompose/Decompose.m

@@ -0,0 +1,27 @@
+classdef Decompose < handle
+
+   properties
+      fullIn % Gaussian object of full molecule (calculation already run)
+      full % full calc 
+      frags % cell array of the fragments
+      fragList % cell array, with lists of atoms in each fragment
+      links % link atoms (currently only 2 supported)
+      rlinks % array with bond lengths for the link atoms
+      keywords % string: should eventually pull from full, but user for now
+
+      maps % full.orbs(maps{ifrag),:) = full MOs in the AO basis of 
+           % ifrag (but without the link atom present)
+      nonLink % nonLink{ifrag} = all AO's on ifrag except those on link atom
+      overlap % overlap{ifrag}(fragMO,fullMO)
+   end
+
+   methods
+      function obj = Decompose(fullIn,fragList,links,rlinks,keywords)
+         obj.fullIn = fullIn;
+         obj.fragList = fragList;
+         obj.links = links;
+         obj.rlinks = rlinks;
+         obj.keywords = keywords;
+      end
+   end
+end

@Decompose/draw.m

@@ -0,0 +1,65 @@
+function draw(obj,piOnly,figNum)
+if (nargin < 2)
+   figNum = 1;
+end
+
+xp = {[0    0.5];
+   [2.25 2.75];
+   [4.5 5]};
+% draw levels
+ft = {obj.frags{1}, obj.full, obj.frags{2}};
+nbasis = [size(obj.frags{1}.orb,1), size(obj.full.orb,1), ...
+   size(obj.frags{2}.orb,1)];
+
+figure(figNum);
+clf;
+hold on;
+for ifrag = 1:3
+   for iorb = 1:nbasis(ifrag)
+      e = ft{ifrag}.Eorb(iorb);
+      if (iorb <= ft{ifrag}.Nelectrons/2)
+         format = 'b'; % filled orbs are blue
+      else
+         format = 'g'; % empty orbs are green
+      end
+      if (piOnly && (piCharacter(ft{ifrag},iorb) < 0.1))
+         format = 'y'; % non-pi will be yellow
+      end
+      plot(xp{ifrag}, [e e], format);
+   end
+end
+
+% draw overlap lines
+dx = 0.0;
+xc = {[xp{1}(2)+dx xp{2}(1)-dx] [xp{3}(1)-dx xp{2}(2)+dx]};
+
+threshold = .25;
+for ifrag = 1:2
+   ol = obj.overlap{ifrag}.^2;
+   Efrag = obj.frags{ifrag}.Eorb;
+   Efull = obj.full.Eorb;
+   for j = 1:length(Efrag)
+      for k = 1:length(Efull)
+         if (~piOnly || (piCharacter(obj.full,k)>0.1))
+            if ol(j,k)>threshold
+               e1 = Efrag(j);
+               e2 = Efull(k);
+               if (k <= obj.full.Nelectrons/2)
+                  format = 'b';
+               else
+                  format = 'g';
+               end
+               plot(xc{ifrag},[e1 e2],[format,':']);
+            end
+         end
+      end
+   end
+end
+
+end
+
+function res = piCharacter(m, iorb)
+  % m is a guassian calc for a molecule lying in x,y plane
+  a1 = m.orb((m.type == 1) & (m.subtype == 3) , iorb);
+  res = sum(a1.^2);
+end

@Decompose/initialize.m

@@ -0,0 +1,76 @@
+function initialize(obj)
+
+% Calculation on the full molecule
+fullList = [obj.fragList{1}(:);obj.fragList{2}(:)];
+[tempdir,jobname] = writeTPL(obj,'full',fullList);
+obj.full = Gaussian([tempdir,'\'],jobname,{});
+obj.full.run();
+% fragment calcs
+%  vector pointing from link1 to link2
+direction = obj.fullIn.rcart(:,obj.links(2)) - ...
+   obj.fullIn.rcart(:,obj.links(1));
+% hydrogen on frag1 is directed away from link 1 along direction
+rLink{1} = obj.fullIn.rcart(:,obj.links(1)) + ...
+   obj.rlinks(1) * direction/norm(direction);
+% hydrogen on frag2 is directed away from link 2 along -direction
+rLink{2} = obj.fullIn.rcart(:,obj.links(2)) - ...
+   obj.rlinks(2) * direction/norm(direction);
+for ifrag = 1:2
+   [tempdir,jobname] = writeTPL(obj,['frag',int2str(ifrag)], ...
+      obj.fragList{ifrag},rLink{ifrag});
+   obj.frags{ifrag} =  Gaussian([tempdir,'\'],jobname,{});
+   obj.frags{ifrag}.run();
+end
+
+% map fragment AO basis to full AO basis
+icount = 1;
+for ifrag = 1:2
+   ft = obj.frags{ifrag};
+   obj.nonLink{ifrag} = find(ft.atom ~= ft.atom(end));
+   lengthNonLink = length(obj.nonLink{ifrag});
+   % assume same ordering of orbs in full as frag, without link
+   obj.maps{ifrag} = icount:(icount + lengthNonLink - 1);
+   icount = icount + lengthNonLink;
+end
+% calculate overlaps
+for ifrag = 1:2
+   % frag(ao,mo)' * S(ao,ao) * full(ao,mo)
+   noLink = obj.nonLink{ifrag};
+   ofrag = obj.frags{ifrag}.orb(noLink,:);
+   Stemp = obj.frags{ifrag}.overlap(noLink,noLink);
+   ofull = obj.full.orb(obj.maps{ifrag},:);
+   obj.overlap{ifrag} = ofrag' * Stemp * ofull;
+end
+end
+
+
+function [tempDir, jobname] = writeTPL(obj,jobname,atoms,rLink)
+newline = char(10);
+syms{1} = 'H'; syms{6} = 'C'; syms{7} = 'N'; syms{8} = 'O';
+syms{15} = 'P'; syms{16} = 'S';
+gjf_file = [jobname,'.tpl'];
+tempDir = tempname(['c:\G09W','\','Scratch']);
+mkdir(tempDir);
+fid1 = fopen([tempDir,'\',gjf_file],'w');
+fwrite(fid1,['%chk=temp.chk',newline]);
+fwrite(fid1,['# ',obj.keywords,' NoSymmetry iop(3/33=4) pop=regular',newline]);
+fwrite(fid1,newline)
+fwrite(fid1,jobname)
+fwrite(fid1,newline)
+fwrite(fid1,newline)
+fwrite(fid1,'0 1');
+fwrite(fid1,newline)
+for iatom = atoms(:)'
+   fwrite(fid1,[' ',syms{obj.fullIn.Z(iatom)},' ']);
+   for ic = 1:3
+      fwrite(fid1,[num2str(obj.fullIn.rcart(ic,iatom)),' ']);
+   end
+   fwrite(fid1,newline)
+end
+if (nargin > 3)
+   fwrite(fid1,[' H ',num2str(rLink(:)'),newline]);
+end
+fwrite(fid1,newline);
+fclose(fid1);
+end
+

testdat/butadiene.tpl

@@ -0,0 +1,53 @@
+%chk=temp.chk
+# b3lyp/sto-3g geom=connectivity
+
+Title Card Required
+
+0 1
+ C              
+ H                  1            B1
+ H                  1            B2    2            A1
+ C                  1            B3    3            A2    2            D1    0
+ H                  4            B4    1            A3    3            D2    0
+ C                  4            B5    1            A4    3            D3    0
+ H                  6            B6    4            A5    1            D4    0
+ C                  6            B7    4            A6    1            D5    0
+ H                  8            B8    6            A7    4            D6    0
+ H                  8            B9    6            A8    4            D7    0
+
+   B1             1.07000000
+   B2             1.07000000
+   B3             1.35520000
+   B4             1.07000000
+   B5             1.54000000
+   B6             1.07000000
+   B7             1.35520000
+   B8             1.07000000
+   B9             1.07000000
+   A1           120.00000000
+   A2           120.00000000
+   A3           120.00000000
+   A4           120.00000000
+   A5           120.00000000
+   A6           120.00000000
+   A7           120.00000000
+   A8           120.00000000
+   D1          -180.00000000
+   D2          -180.00000000
+   D3             0.00000000
+   D4             0.00000000
+   D5          -180.00000000
+   D6           180.00000000
+   D7             0.00000000
+
+ 1 2 1.0 3 1.0 4 2.0
+ 2
+ 3
+ 4 5 1.0 6 1.0
+ 5
+ 6 7 1.0 8 2.0
+ 7
+ 8 9 1.0 10 1.0
+ 9
+ 10
+

verifyDecompose.m

@@ -0,0 +1,37 @@
+%% start by doing an optimization of butadiene
+clear classes
+qmatlab = pwd;
+
+if (1)
+    gstart = Gaussian([qmatlab, '\testdat\'],'butadiene',{});
+    gstart.run();
+    save([qmatlab, '\testdat\gstart.mat']);
+else
+    load([qmatlab, '\testdat\gstart.mat']);
+end
+%%
+
+if (1)
+    fragList{1} = 1:5;
+    fragList{2} = 6:10;
+    links = [4 6];
+    rlinks = [1.1 1.1];
+    keywords = 'b3lyp/sto-3g';
+    obj = Decompose(gstart,fragList,links,rlinks,keywords);
+    obj.initialize();
+    save([qmatlab, '\testdat\gtemp2.mat']);
+else
+    load([qmatlab, '\testdat\gtemp2.mat']);
+end
+
+disp(['basis size: full ',num2str(length(obj.full.atom)), ...
+   ' frag1 ',num2str(length(obj.frags{1}.atom)), ...   
+   ' frag2 ',num2str(length(obj.frags{1}.atom))]);
+disp(['frag1 ',num2str(obj.maps{1}(:)')]);
+disp(['frag2 ',num2str(obj.maps{2}(:)')]);
+
+s1 = obj.overlap{1};
+s2 = obj.overlap{2};
+max(max(abs(s1.^2 -s2.^2)))
+
+obj.draw(0,10);