Plotting a 3d hydrogen atomic orbital plot from a data file

I have a data file of two colums, column1 has data of radial coordinate 'r' and column2 has data of 'wavefunction' of 3d orbital both are 1D arrays of each 45 entries. I need to obtain a 3D plot in cartesian coordinates 'wavefunction(x,y,z)' after multiplying the corresponding spherical harmonic functions of 3dz^2. How to get one. I am very new to MATLAB any help would be appreciated.

The data is as follows:

0.00 0.000000

0.45 0.000627

0.90 0.002467

1.35 0.005404

1.80 0.009263

2.25 0.013828

2.70 0.018859

3.15 0.024107

3.60 0.029329

4.05 0.034296

4.50 0.038798

4.95 0.042653

5.40 0.045707

5.85 0.047850

6.30 0.049018

6.75 0.049205

7.20 0.048467

7.65 0.046926

8.10 0.044759

8.55 0.042183

9.00 0.039422

9.45 0.036660

9.90 0.033972

10.35 0.031248

10.80 0.028078

11.25 0.023511

11.70 0.016026

12.15 0.005061

12.60 -.006867

13.05 -.014665

13.50 -.014647

13.95 -.007211

14.40 0.003441

14.85 0.011842

15.30 0.015687

15.75 0.016839

16.20 0.016349

16.65 0.014982

17.10 0.013238

17.55 0.011416

18.00 0.009680

18.45 0.008110

18.90 0.006736

19.35 0.005560

19.80 0.004567

[Umar]

Hi Rasheed,

To generate a 3D plot of the wavefunction in Cartesian coordinates, you need to perform the following steps:

Start by loading the data from the two columns into MATLAB. You can do this by creating two separate arrays for 'r' and the wavefunction values.

% Define the data

r = [0.00, 0.45, 0.90, 1.35, 1.80, 2.25, 2.70, 3.15, 3.60, 4.05, 4.50, 4.95, 5.40, 5.85, 6.30, 6.75, 7.20, 7.65, 8.10, 8.55, 9.00, 9.45, 9.90, 10.35, 10.80, 11.25, 11.70, 12.15, 12.60, 13.05, 13.50, 13.95, 14.40, 14.85, 15.30, 15.75, 16.20, 16.65, 17.10, 17.55, 18.00, 18.45, 18.90, 19.35, 19.80];

wavefunction = [0.000000, 0.000627, 0.002467, 0.005404, 0.009263, 0.013828, 0.018859, 0.024107, 0.029329, 0.034296, 0.038798, 0.042653, 0.045707, 0.047850, 0.049018, 0.049205, 0.048467, 0.046926, 0.044759, 0.042183, 0.039422, 0.036660, 0.033972, 0.031248, 0.028078, 0.023511, 0.016026, 0.005061, -0.006867, -0.014665, -0.014647, -0.007211, 0.003441, 0.011842, 0.015687, 0.016839, 0.016349, 0.014982, 0.013238, 0.011416, 0.009680, 0.008110, 0.006736, 0.005560, 0.004567];

Next, you need to calculate the spherical harmonic function for the 3dz^2 orbital. The spherical harmonic function for this orbital is given by Y_3^0(theta, phi) = sqrt(7/16pi) * (3cos^2(theta) - 1), where theta is the polar angle and phi is the azimuthal angle.

% Calculate the spherical harmonic function

theta = acos(r./max(r)); % Calculate the polar angle

phi = linspace(0, 2*pi, length(r)); % Generate the azimuthal angle

Y_3_0 = sqrt(7/(16*pi)) * (3*cos(theta).^2 - 1

Multiply the wavefunction values by the corresponding spherical harmonic function values to obtain the modified wavefunction.

% Multiply the wavefunction by the spherical harmonic function

modified_wavefunction = wavefunction .* Y_3_0;

Convert the modified wavefunction from spherical coordinates to Cartesian coordinates using the following equations:

x = r * sin(theta) * cos(phi)

y = r * sin(theta) * sin(phi)

z = r * cos(theta)

% Convert to Cartesian coordinates

x = r .* sin(theta) .* cos(phi);

y = r .* sin(theta) .* sin(phi);

z = r .* cos(theta);

Finally, use the 'plot3' function to generate the 3D plot of the wavefunction in Cartesian coordinates.

% Generate the 3D plot

figure;

plot3(x, y, z, 'b.');

xlabel('x');

ylabel('y');

zlabel('z');

title('3D Plot of the Wavefunction');

grid on;

Feel free to customize code based on your preferences. I hope this answers your question.

[Cris LaPierre]

As a runnable example

% Define the data
r = [0.00, 0.45, 0.90, 1.35, 1.80, 2.25, 2.70, 3.15, 3.60, 4.05, 4.50, 4.95, 5.40, 5.85, 6.30, 6.75, 7.20, 7.65, 8.10, 8.55, 9.00, 9.45, 9.90, 10.35, 10.80, 11.25, 11.70, 12.15, 12.60, 13.05, 13.50, 13.95, 14.40, 14.85, 15.30, 15.75, 16.20, 16.65, 17.10, 17.55, 18.00, 18.45, 18.90, 19.35, 19.80];
wavefunction = [0.000000, 0.000627, 0.002467, 0.005404, 0.009263, 0.013828, 0.018859, 0.024107, 0.029329, 0.034296, 0.038798, 0.042653, 0.045707, 0.047850, 0.049018, 0.049205, 0.048467, 0.046926, 0.044759, 0.042183, 0.039422, 0.036660, 0.033972, 0.031248, 0.028078, 0.023511, 0.016026, 0.005061, -0.006867, -0.014665, -0.014647, -0.007211, 0.003441, 0.011842, 0.015687, 0.016839, 0.016349, 0.014982, 0.013238, 0.011416, 0.009680, 0.008110, 0.006736, 0.005560, 0.004567];

% Calculate the spherical harmonic function
theta = acos(r./max(r)); % Calculate the polar angle
phi = linspace(0, 2*pi, length(r)); % Generate the azimuthal angle
Y_3_0 = sqrt(7/(16*pi)) * (3*cos(theta).^2 - 1);

% Multiply the wavefunction by the spherical harmonic function
modified_wavefunction = wavefunction .* Y_3_0;

% Convert to Cartesian coordinates
x = r .* sin(theta) .* cos(phi);
y = r .* sin(theta) .* sin(phi);
z = r .* cos(theta);

% Generate the 3D plot
figure;
plot3(x, y, z, 'b.');
xlabel('x');
ylabel('y');
zlabel('z');
title('3D Plot of the Wavefunction');
grid on;

[Rasheed Shaik]

Thank you for the help. The plot3 (x,y,z,'b') does not supply the modified_wavefunction data for plotting am i right?

[Umar]

Hi Rashid,

Since the modified_wavefunction data was not provided, hence the 3d plot. But I am pretty confident that you can figure out that part. I put a lot of effort into it, so I hope this above code snippet answers your question.

[Umar]

Hi Rasheed,

I do apologize since I was exhausted on working on a large scale software project. Thanks for noticing, there was a missing line in my code which was multiplying the wavefunction values with the corresponding spherical harmonic functions of 3dz^2. However, in this modified code, I fixed it.

To answer your question in the post, to create a 3D plot in Cartesian coordinates using your given dataset, you need to convert the radial coordinates 'r' to Cartesian coordinates (x, y, z) and then multiply the wavefunction values with the corresponding spherical harmonic functions of 3dz^2. Then, extract the radial coordinates 'r' and wavefunction values from the dataset into separate variables r and wavefunction. Next, convert the radial coordinates to Cartesian coordinates using the formulas x = r * sin(r), y = r * cos(r), and z = r. After that, multiply the wavefunction values with the corresponding spherical harmonic functions of 3dz^2 to obtain the wavefunction_cartesian values. Finally, create a 3D scatter plot using the scatter3 function, where the x, y, and z coordinates represent the Cartesian coordinates, and the wavefunction_cartesian values will determine the color of each point. I added a colorbar to the plot to indicate the color scale. The xlabel, ylabel, and zlabel functions are used to label the axes, and the title function sets the title of the plot. I also adjusted the plot settings using axis equal to ensure equal scaling on all axes and view to set the viewing angle for better visualization. Now, by running this code in MATLAB, you should obtain a 3D plot of the wavefunction in Cartesian coordinates, where the color represents the magnitude of the wavefunction. Here is the snippet code,

 % Load the dataset

data = [0.00 0.000000 0.45 0.000627 0.90 0.002467 1.35 0.005404 1.80 0.009263 2.25 0.013828 2.70 0.018859 3.15 0.024107 3.60 0.029329 4.05 0.034296 4.50 0.038798 4.95 0.042653 5.40 0.045707 5.85 0.047850 6.30 0.049018 6.75 0.049205 7.20 0.048467 7.65 0.046926 8.10 0.044759 8.55 0.042183 9.00 0.039422 9.45 0.036660 9.90 0.033972 10.35 0.031248 10.80 0.028078 11.25 0.023511 11.70 0.016026 12.15 0.005061 12.60 -0.006867 13.05 -0.014665 13.50 -0.014647 13.95 -0.007211 14.40 0.003441 14.85 0.011842 15.30 0.015687 15.75 0.016839 16.20 0.016349 16.65 0.014982 17.10 0.013238 17.55 0.011416 18.00 0.009680 18.45 0.008110 18.90 0.006736 19.35 0.005560 19.80 0.004567];

% Extract the radial coordinates 'r' and wavefunction values

r = data(:, 1);

wavefunction = data(:, 2);

% Convert radial coordinates to Cartesian coordinates

x = r .* sin(r);

y = r .* cos(r);

z = r;

% Multiply wavefunction values with the corresponding spherical harmonic

functions of 3dz^2

modified_wavefunction = wavefunction .* (3 * z.^2);

% Create a 3D plot

figure;

scatter3(x, y, z, 50, modified_wavefunction, 'filled');

colorbar;

xlabel('x');

ylabel('y');

zlabel('z');

title('3D Plot of wavefunction(x, y, z)');

% Adjust the plot settings for better visualization

axis equal;

view(45, 30);

Please see attached snippet code and plot.

I hope this will make you happy now. Please let me know if you have further questions.

Answer 1

[Karan Singh]

Hi @Rasheed Shaik,

Is this what is required?

% Data
r = [0.00, 0.45, 0.90, 1.35, 1.80, 2.25, 2.70, 3.15, 3.60, 4.05, 4.50, 4.95, 5.40, 5.85, 6.30, 6.75, 7.20, 7.65, 8.10, 8.55, 9.00, 9.45, 9.90, 10.35, 10.80, 11.25, 11.70, 12.15, 12.60, 13.05, 13.50, 13.95, 14.40, 14.85, 15.30, 15.75, 16.20, 16.65, 17.10, 17.55, 18.00, 18.45, 18.90, 19.35, 19.80];
wavefunction = [0.000000, 0.000627, 0.002467, 0.005404, 0.009263, 0.013828, 0.018859, 0.024107, 0.029329, 0.034296, 0.038798, 0.042653, 0.045707, 0.047850, 0.049018, 0.049205, 0.048467, 0.046926, 0.044759, 0.042183, 0.039422, 0.036660, 0.033972, 0.031248, 0.028078, 0.023511, 0.016026, 0.005061, -0.006867, -0.014665, -0.014647, -0.007211, 0.003441, 0.011842, 0.015687, 0.016839, 0.016349, 0.014982, 0.013238, 0.011416, 0.009680, 0.008110, 0.006736, 0.005560, 0.004567];

% Define the spherical harmonic function for 3d orbital (3dz^2)
Y_3dz2 = @(theta, phi) sqrt(5/(16*pi)) * (3*cos(theta).^2 - 1);

% Create a meshgrid in spherical coordinates
[theta, phi] = meshgrid(linspace(0, pi, 45), linspace(0, 2*pi, 45));

% Interpolate wavefunction data to the meshgrid
r_interp = interp1(linspace(0, max(r), length(wavefunction)), wavefunction, linspace(0, max(r), 45));

% Compute the wavefunction in spherical coordinates
wavefunction_spherical = r_interp' .* Y_3dz2(theta, phi);

% Convert spherical to Cartesian coordinates
[x, y, z] = sph2cart(phi, pi/2 - theta, wavefunction_spherical);

% Plot the 3D wavefunction
figure;
surf(x, y, z, wavefunction_spherical);
xlabel('x');
ylabel('y');
zlabel('z');
title('3D Plot of Wavefunction in Cartesian Coordinates');
colorbar;
shading interp;

[Rasheed Shaik]

With minor edit in the code in the line 'wavefunction_spherical' I could able to achieve sensible plot. Thank you.

[Umar]

Hi @Rasheed Shaik,

Do you want Karan Singh to help you with minor edit in the code, please let us know. If there are any further questions, please let us know.

[Umar]

Hi @Rasheed,

Hi Rasheed,

I made some modifications to Karen Singh’s code by adding abs() to ensure positive values for plotting the wavefunction in spherical coordinates. This adjustment will help in generating a more visually meaningful and accurate representation of the 3D wavefunction in Cartesian coordinates, the view(0, 0) command, I set the viewing angle to a horizontal position, providing a better perspective from all angles, resembling a 3D hydrogen atomic orbital. This adjustment enhances the visualization for a more comprehensive and visually appealing representation.

>> % Data r = [0.00, 0.45, 0.90, 1.35, 1.80, 2.25, 2.70, 3.15, 3.60, 4.05, 4.50, 4.95, 5.40, 5.85, 6.30, 6.75, 7.20, 7.65, 8.10, 8.55, 9.00, 9.45, 9.90, 10.35, 10.80, 11.25, 11.70, 12.15, 12.60, 13.05, 13.50, 13.95, 14.40, 14.85, 15.30, 15.75, 16.20, 16.65, 17.10, 17.55, 18.00, 18.45, 18.90, 19.35, 19.80]; wavefunction = [0.000000, 0.000627, 0.002467, 0.005404, 0.009263, 0.013828, 0.018859, 0.024107, 0.029329, 0.034296, 0.038798, 0.042653, 0.045707, 0.047850, 0.049018, 0.049205, 0.048467, 0.046926, 0.044759, 0.042183, 0.039422, 0.036660, 0.033972, 0.031248, 0.028078, 0.023511, 0.016026, 0.005061, -0.006867, -0.014665, -0.014647, -0.007211, 0.003441, 0.011842, 0.015687, 0.016839, 0.016349, 0.014982, 0.013238, 0.011416, 0.009680, 0.008110, 0.006736, 0.005560, 0.004567];

% Define the spherical harmonic function for 3d orbital (3dz^2)

Y_3dz2 = @(theta, phi) sqrt(5/(16*pi)) * (3*cos(theta).^2 - 1);

% Create a meshgrid in spherical coordinates

[theta, phi] = meshgrid(linspace(0, pi, 45), linspace(0, 2*pi, 45));

% Interpolate wavefunction data to the meshgrid

r_interp = interp1(linspace(0, max(r), length(wavefunction)), wavefunction, linspace(0, max(r), 45));

% Compute the wavefunction in spherical coordinates

wavefunction_spherical = r_interp' .* Y_3dz2(theta, phi);

% Convert spherical to Cartesian coordinates

[x, y, z] = sph2cart(phi, pi/2 - theta,

abs(wavefunction_spherical)); % Use abs() to ensure positive values for plotting

% Plot the 3D wavefunction with a horizontal view figure;

surf(x, y, z, wavefunction_spherical);

view(0, 0); % Set the view to a horizontal angle

xlabel('x');

ylabel('y');

zlabel('z');

title('3D Plot of Wavefunction in Cartesian Coordinates (Horizontal View)');

colorbar;

shading interp;

Please see attached plot.