Molecular Vibrations - Numerical

[CNX]

Homework Statement

I'm trying to do some numerical stuff with vibrations of H20 and I'm working in mdyne, angstroms, atomic mass units, and angles are given in radians. What would the corresponding unit of time be when I calculate my normal mode frequencies? femtosecond, 10e-15?

 

[CNX]

Here is how I'm setting up my eigenvalue problem:

I know how to construct the G matrix elements.

Here's how I construct the F matrix elements:

V(Q_1, Q_2, Q_3) = \sum_{i,k,j} K_{i,k,j} (Q_1)^i (Q_2)^j (Q_2)^k

K_{i,k,j} can be obtained from this table:

i j k Kijk Units
2 0 0 4.227 mdyne °A−1
0 2 0 4.227 mdyne °A−1
0 0 2 0.349 mdyne °A
1 1 0 -0.101 mdyne °A−1
1 0 1 0.219 mdyne
0 1 1 0.219 mdyne

And then I construct the G and F matrices in matlab:

% Speed of light in a vacuum.
c = 2.998e10; % [cm/s]

% AMU from http://www.sisweb.com/referenc/source/exactmaa.htm
m1 = 1.007825; % Mass of hydrogen [amu]
m2 = 15.994915; % Mass of oxygen [amu]
m3 = 1.007825; % Mass of hydrogen [amu]

% Equilibrium values.
phi = 104.54*pi/180; % H20 equilibrium bend angle [rad]
r = 0.9584; % (r12=r32=r) H20 equilibrium stretch [Angstroms]

%% Determining normal mode frequencies and forms.

% Construct G matrix elements.
g11 = 1/m1 + 1/m2;
g22 = 1/m2 + 1/m3;
g12 = cos(phi)/m2;
g13 = -sin(phi)/(m2*r);
g23 = -sin(phi)/(m2*r);
g33 = 1/(m1*r^2) + 1/(m3*r^2) + 1/m2*(1/r^2 + 1/r^2 - 2*cos(phi)/r^2);

% Construct G and F matrices.
G = [[g11 g12 g13]; [g12 g22 g23]; [g13 g23 g33]];
F = [[4.227 -0.101 0.219]; [-0.101 4.227 0.219]; [0.219 0.219 0.349]];

% Solve the eigenvalue/eigenvector problem (GF-omega^2*1).a=0.
[V, omega2] = eig(G*F);

% Divide by 10e-15 to go from femtosec->sec
omega = sqrt(omega2)/1e-15;

% Get wavenumbers
wv = omega/(2*pi*c); % [cm^-1];

i.e.,

(GF-\omega^2 1)\cdot a = 0

My result for the wavenumbers: 4617, 11426, 10994 cm^{-1}. They should be: 3657, 3756, 1595 cm^{-1}

I swear I had good numbers coming out of this and now it just doesn't make sense. Any help is appreciated.

 

[Redbelly98]

Homework Statement

I'm trying to do some numerical stuff with vibrations of H20 and I'm working in mdyne, angstroms, atomic mass units, and angles are given in radians. What would the corresponding unit of time be when I calculate my normal mode frequencies? femtosecond, 10e-15?

Have to think about this one a little.

The "given" units are force, distance, and mass.

We know that Force units, in terms of Mass, Distance, and Time, are

F = M*D/T²​

(Eg., 1 Newton = 1 kg*m/s²)

Solving for Time units, we have

T = (M*D/F)^1/2​

If you convert mdyne, Angstroms, and amu into N, m, and kg, you can calculate the time unit in seconds. The answer is comparable, but not equal, to a femtosecond.

EDIT:
Two things I see wrong with your code.
1. You're using 1 fs for the time units, which is wrong.
2. You're using cm/sec for c. Use Angstroms and your correct time unit instead.

 

[CNX]

For the unit of time calculation I get

>> (1.66e-27*1e-10/1e-11)^(1/2)

ans =

     1.28840987267251e-013

Now, my intention was to convert the \omega values from the eigenvalue problem from [rad/whatever] to [rad/s]. Then I was using the speed of light in [cm/s] to convert these angular frequencies in [rad/s] to wave numbers in [cm^{-1}] for comparison with the accepted values I have (see OP.)

Here is the modified code (the matrix construction is more drawn-out but should be similar--though I still don't know about the factors of 1/2 for off-diagonal elements of matrix F):

%% Define Constants

% Define the accepted values for normal mode wave numbers of H20.
theory = [1594 3657 3756] % [cm^-1];

% Speed of light in a vacuum.
c = 2.998e10; % [cm/s]

% AMU from http://www.sisweb.com/referenc/source/exactmaa.htm
m1 = 1.007825; % Mass of hydrogen [amu]
m2 = 15.994915; % Mass of oxygen [amu]
m3 = 1.007825; % Mass of hydrogen [amu]

% So m2 is oxygen, with two hydrogens attached to it to form H20.

% Equilibrium values.
phi = 104.54*pi/180; % H20 equilibrium bend angle [rad]
r12 = 0.9584; % H20 equilibrium stretch (between m1 and m2) [Angstroms]
r32 = 0.9584; % H20 equilibrium stretch (between m3 and m2) [Angstroms]

%% Matrix Construction

% Construct G matrix elements.
g11 = 1/m1 + 1/m2;
g22 = 1/m2 + 1/m3;
g33 = 1/(m1*r12^2) + 1/(m3*r32^2) + 1/m2*( 1/r32^2 + 1/r12^2 - (2*cos(phi))/(r12*r32) );

g12 = cos(phi)/m2;
g21 = g12; % G matrix is symmetric.
g13 = -sin(phi)/(m2*r32);
g31 = g13;
g23 = -sin(phi)/(m2*r12);
g32 = g23;

% Construct G matrix
G = [[g11 g12 g13]; [g21 g22 g23]; [g31 g32 g33]];

% Construct F matrix elements.
f11 = 4.227;
f22 = 4.227;
f33 = 0.349;

f12 = -0.101/2;
f21 = f12;
f13 = 0.219/2;
f31 = f13;
f23 = 0.219/2;
f32 = f23;

% Construct F matrix.
F = [[f11 f12 f13]; [f21 f22 f23]; [f31 f32 f33]];

%% Solving the eigenvalue/eigenvector problem (GF-omega^2*1).a=0.

% Get eigenvectors and eigenvalues.
[V, omega2] = eig(G*F);

% Get angular frequencies and convert to rad/s.
omega = sqrt(omega2)/1.28840987267251e-013; % [rad/s]

% Calculate wavenumbers to compare with theory.
wn = omega/(2*pi*c); % [cm^-1];

Unfortunately the wave number results are 88, 86, 37 [cm^{-1}]

The odd thing is when I was playing around I noticed that replacing the 1.288e-013s unit of time with 2 \times 1.66e-15s when trying to convert to [rad/s], I get (after converting to wave number),

1430, 3325, 3421 [cm^{-1}]

Again, theory: 1594, 3657, 3756 [cm^{-1}]

i.e. I know it's random but they're the closest values I've seen yet and I'm not sure why. If all else is OK then my unit of time is two orders of magnitude off to match the accepted values... I feel like I'm soo close.

 

[Redbelly98]

For the unit of time calculation I get

>> (1.66e-27*1e-10/1e-11)^(1/2)

ans =

     1.28840987267251e-013

1 millidyne is 1e-8 N, not 1e-11. Or did you mean microdynes?

Looks like you still need to convert the speed of light into Angstroms/(your time unit) in your code.
[EDIT: Never mind, I just realized c is only used for a final units conversion]

p.s. Welcome to Physics Forums

 

[CNX]

Oops I misread the conversion for mdyne. The unit of time is,

>> sqrt(1.66e-27*1e-10/1e-8)

ans =

     4.07430975749267e-015

and the resulting wavenumbers,

           2788.1297896707                         0                         0
                         0          2709.85261638262                         0
                         0                         0          1165.45560695961

I'm not sure how close I can expect these to get but right now they're sitting at about 25% difference from accepted values.

 

[Redbelly98]

Well, at least you're in the ballpark now. Can't help anymore than I have, sorry.

 

[CNX]

Thanks for your help -- I'll keep working at it