Generation of hybridised orbitals from s- and p-atomic orbitals Matlab code 

% Program to plot different hybridized orbitals sp, sp2, sp3, sp + 2

% p-orbitals (as in ethyne), and sp2 + 1 p-orbital (as in ethene).

% The lobes are purely schematic and do not follow true quantum-mechanical

% wavefunctions, as these are anyway approximate for multi-electron atoms


clear; close all


disp('Choose one of the following orbitals: sp1, sp2, sp3,');

prompt = 'sp1pp (sp + 2 p-orbitals), sp2p (sp2 + 1 p-orbital) ';

str1 = input(prompt,'s');

str2 = '-orbital.mp4';

fileName = [str1,str2]; % Name of file depends on type of crystal sample

vid = VideoWriter(fileName,'MPEG-4');

vid.Quality = 100;

vid.FrameRate = 30;

open(vid);


figure('units','pixels','position',[0 0 1920 1080],'ToolBar','none');

set(0,'defaultfigurecolor',[1 1 1]);


orng = [1 0.5 0]; % Orange

lp = [-2 -1 1]; % Lighting angle

view(0,16); % Viewing angle

L = 0.7; % Long lobe

S = 0.5; % Short lobe

F = 1.4; % Fat lobe

T = 1.0; % Thin lobe

scale = 3; % Amount to reduce opposite lobe of hybridized orbitals

opac1 = 0.8; % Opacity of lobe

opac2 = 0.64; % Alternative opacity of lobe

sz = [5 5 3.5 5 5]; % Size of transparent lines to keep fixed xyz limits



% Set up parameters for drawing orbital lobes

% First index defines which lobe (1 to 4); the second index defines the

% five parameters associated with that lobe (phi, theta axis-rotation functions,

% length of lobe, width of lobe, opacity); the third index defines which of the five

% orbitals to plot (sp1, sp2, sp3, sp+pp, sp2+p)

%

% sp1

pars(1,:,1) = [90 0 L T opac2]; % Lobe-1 parameters, sp-orbital

pars(2,:,1) = [-90 0 L T opac2]; % Lobe-2 parameters, sp-orbital

pars(3,:,1) = [0 0 L T 0]; % Lobe-3 parameters, sp-orbital. Doesn't exist, so opacity = 0

pars(4,:,1) = [0 0 L T 0]; % Lobe-4 parameters, sp-orbital. Doesn't exist, so opacity = 0

% sp2

pars(1,:,2) = [90 0 L T opac1]; % Lobe-1 parameters, sp2-orbital

pars(2,:,2) = [90 120 L T opac1]; % Lobe-2 parameters, sp2-orbital

pars(3,:,2) = [90 240 L T opac1]; % Lobe-3 parameters, sp2-orbital

pars(4,:,2) = [0 0 L T 0]; % Lobe-4 parameters, sp2-orbital. Doesn't exist, so opacity = 0

% sp3

pars(1,:,3) = [0 0 S T opac1]; % Lobe-1 parameters, sp3-orbital

pars(2,:,3) = [109.5 0 S T opac1]; % Lobe-2 parameters, sp3-orbital

pars(3,:,3) = [109.5 120 S T opac1]; % Lobe-3 parameters, sp3-orbital

pars(4,:,3) = [109.5 240 S T opac1]; % Lobe-4 parameters, sp3-orbital

% sp1pp

pars(1,:,4) = [90 0 L T opac2]; % Lobe-1 parameters, sp-orbital

pars(2,:,4) = [-90 0 L T opac2]; % Lobe-2 parameters, sp-orbital

pars(3,:,4) = [0 0 S F opac1]; % Lobe-3 parameters, p-orbital

pars(4,:,4) = [90 90 S F opac1]; % Lobe-4 parameters, p-orbital

% sp2p

pars(1,:,5) = [90 0 L T opac1]; % Lobe-1 parameters, sp2-orbital

pars(2,:,5) = [90 120 L T opac1]; % Lobe-2 parameters, sp2-orbital

pars(3,:,5) = [90 240 L T opac1]; % Lobe-3 parameters, sp2-orbital

pars(4,:,5) = [0 0 S F opac1]; % Lobe-4 parameters, p-orbital


% Lobe equations

conecoord1 = 0:pi/320:pi/2; % coordinates of light cone schematic

conecoord2 = 0:pi/150:2*pi; % coordinates of light cone schematic

[conecoord1,conecoord2]=meshgrid(conecoord1,conecoord2);

xcoordcone1 = 4*(1 - sin(conecoord1)).*sin(conecoord1).*cos(conecoord2);

ycoordcone1 = 4*(1 - sin(conecoord1)).*sin(conecoord1).*sin(conecoord2);

zcoordcone1 = 7*(abs(cos(conecoord1))).^3;


if strcmp('sp1',str1)

    orb = 1;

elseif strcmp('sp2',str1)

    orb = 2;

elseif strcmp('sp3',str1)

    orb = 3;

elseif strcmp('sp1pp',str1)

    orb = 4;

elseif strcmp('sp2p',str1)

    orb = 5;

else

    disp('Unknown orbital');

end


%_________________________

% Plotting part of program

for i = 0:2:358 % Rotation around vertical axis

    % Lobes 1 - 4 in blue

    lobe1 = surf(xcoordcone1*pars(1,4,orb),ycoordcone1*pars(1,4,orb),...

        zcoordcone1*pars(1,3,orb),'FaceAlpha',pars(1,5,orb),'FaceColor',...

        'b','LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);

    rotate(lobe1,[1 0 0],pars(1,1,orb),[0 0 0]);

    rotate(lobe1,[0 0 1],pars(1,2,orb)+i,[0 0 0]);

    hold on

    lobe2 = surf(xcoordcone1*pars(2,4,orb),ycoordcone1*pars(2,4,orb),...

        zcoordcone1*pars(2,3,orb),'FaceAlpha',pars(2,5,orb),'FaceColor',...

        'b','LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);

    rotate(lobe2,[1 0 0],pars(2,1,orb),[0 0 0]);

    rotate(lobe2,[0 0 1],pars(2,2,orb)+i,[0 0 0]);


    lobe3 = surf(xcoordcone1*pars(3,4,orb),ycoordcone1*pars(3,4,orb),...

        zcoordcone1*pars(3,3,orb),'FaceAlpha',pars(3,5,orb),'FaceColor',...

        'b','LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);

    rotate(lobe3,[1 0 0],pars(3,1,orb),[0 0 0]);

    rotate(lobe3,[0 0 1],pars(3,2,orb)+i,[0 0 0]);


    lobe4 = surf(xcoordcone1*pars(4,4,orb),ycoordcone1*pars(4,4,orb),...

        zcoordcone1*pars(4,3,orb),'FaceAlpha',pars(4,5,orb),'FaceColor',...

        'b','LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);

    rotate(lobe4,[1 0 0],pars(4,1,orb),[0 0 0]);

    rotate(lobe4,[0 0 1],pars(4,2,orb)+i,[0 0 0]);


    %------------------------------------------------------------------------

    % Opposite lobes to lobes 1 - 4 in orange


    lobe5 = surf(xcoordcone1*pars(1,4,orb)/scale,ycoordcone1*pars(1,4,orb)/scale,...

        -zcoordcone1*pars(1,3,orb)/scale,'FaceAlpha',pars(1,5,orb),'FaceColor',...

        orng,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);

    rotate(lobe5,[1 0 0],pars(1,1,orb),[0 0 0]);

    rotate(lobe5,[0 0 1],pars(1,2,orb)+i,[0 0 0]);


    lobe6 = surf(xcoordcone1*pars(2,4,orb)/scale,ycoordcone1*pars(2,4,orb)/scale,...

        -zcoordcone1*pars(2,3,orb)/scale,'FaceAlpha',pars(2,5,orb),'FaceColor',...

        orng,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);

    rotate(lobe6,[1 0 0],pars(2,1,orb),[0 0 0]);

    rotate(lobe6,[0 0 1],pars(2,2,orb)+i,[0 0 0]);


    if strcmp('sp1pp',str1)

        scale2 = 1; % Unhybridized p-orbital

    else

        scale2 = scale; % Hybridized

    end

    lobe7 = surf(xcoordcone1*pars(3,4,orb)/scale2,ycoordcone1*pars(3,4,orb)/scale2,...

        -zcoordcone1*pars(3,3,orb)/scale2,'FaceAlpha',pars(3,5,orb),'FaceColor',...

        orng,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);

    rotate(lobe7,[1 0 0],pars(3,1,orb),[0 0 0]);

    rotate(lobe7,[0 0 1],pars(3,2,orb)+i,[0 0 0]);


    if strcmp('sp1pp',str1) || strcmp('sp2p',str1)

        scale2 = 1; % Unhybridized p-orbital

    else

        scale2 = scale; % Hybridized

    end

    lobe8 = surf(xcoordcone1*pars(4,4,orb)/scale2,ycoordcone1*pars(4,4,orb)/scale2,...

        -zcoordcone1*pars(4,3,orb)/scale2,'FaceAlpha',pars(4,5,orb),'FaceColor',...

        orng,'LineStyle','none','FaceLighting','gouraud','DiffuseStrength',1);

    rotate(lobe8,[1 0 0],pars(4,1,orb),[0 0 0]);

    rotate(lobe8,[0 0 1],pars(4,2,orb)+i,[0 0 0]);


    % Dummy horizontal and vertical lines to stop image resizing

    % Four-element color definition [rgb alpha] only works with MATLAB 2021

    % onwards

    plot3([sz(orb) -sz(orb)],[0 0],[-sz(orb) -sz(orb)],...

        'color',[1 1 1 0], 'LineWidth', 0.001);

    plot3([0 0],[sz(orb) -sz(orb)],[-sz(orb) -sz(orb)],...

        'color',[1 1 1 0], 'LineWidth', 0.001);

    plot3([-sz(orb) -sz(orb)],[0 0],[sz(orb) -sz(orb)],...

        'color',[1 1 1 0], 'LineWidth', 0.001);


    axis off

    axis equal

    axis tight

    light('Position',lp,'Style','infinite');


    frame = getframe(gcf);

    writeVideo(vid,frame);


    hold off

end

% Output the movie as an mpg file

close(vid);