From aef23bd3121b9c4eaf4ffdcfbf27661df0920af7 Mon Sep 17 00:00:00 2001 From: Christopher Collins Date: Tue, 18 Jun 2013 03:08:55 -0400 Subject: [PATCH] New decomposition code --- @Decompose/Decompose.m | 27 +++++++++++++++ @Decompose/draw.m | 65 +++++++++++++++++++++++++++++++++++ @Decompose/initialize.m | 76 +++++++++++++++++++++++++++++++++++++++++ testdat/butadiene.tpl | 53 ++++++++++++++++++++++++++++ verifyDecompose.m | 37 ++++++++++++++++++++ 5 files changed, 258 insertions(+) create mode 100644 @Decompose/Decompose.m create mode 100644 @Decompose/draw.m create mode 100644 @Decompose/initialize.m create mode 100644 testdat/butadiene.tpl create mode 100644 verifyDecompose.m diff --git a/@Decompose/Decompose.m b/@Decompose/Decompose.m new file mode 100644 index 0000000..65ab9d1 --- /dev/null +++ b/@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 diff --git a/@Decompose/draw.m b/@Decompose/draw.m new file mode 100644 index 0000000..3282b79 --- /dev/null +++ b/@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 diff --git a/@Decompose/initialize.m b/@Decompose/initialize.m new file mode 100644 index 0000000..ceb032b --- /dev/null +++ b/@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 + diff --git a/testdat/butadiene.tpl b/testdat/butadiene.tpl new file mode 100644 index 0000000..fa36d4f --- /dev/null +++ b/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 + diff --git a/verifyDecompose.m b/verifyDecompose.m new file mode 100644 index 0000000..1f006fe --- /dev/null +++ b/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);