1. Limited time only! Sign up for a free 30min personal tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Molecular Vibrations - Numerical

  1. Jul 23, 2009 #1

    CNX

    User Avatar

    1. The problem statement, all variables and given/known data

    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?
     
  2. jcsd
  3. Jul 23, 2009 #2

    CNX

    User Avatar

    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:

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

    [itex]K_{i,k,j}[/itex] can be obtained from this table:
    Code (Text):

    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:

    Code (Text):


    % 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.,

    [tex](GF-\omega^2 1)\cdot a = 0[/tex]

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

    I swear I had good numbers coming out of this and now it just doesn't make sense. Any help is appreciated.
     
    Last edited: Jul 23, 2009
  4. Aug 1, 2009 #3

    Redbelly98

    User Avatar
    Staff Emeritus
    Science Advisor
    Homework Helper

    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/T2

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

    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.
     
    Last edited: Aug 1, 2009
  5. Aug 1, 2009 #4

    CNX

    User Avatar

    For the unit of time calculation I get

    Code (Text):

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

    ans =

         1.28840987267251e-013
     
    Now, my intention was to convert the [itex]\omega[/itex] 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 [itex][cm^{-1}][/itex] 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):

    Code (Text):

    %% 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 [[itex]cm^{-1}[/itex]]

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

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

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

    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.
     
    Last edited: Aug 1, 2009
  6. Aug 2, 2009 #5

    Redbelly98

    User Avatar
    Staff Emeritus
    Science Advisor
    Homework Helper

    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 :smile:
     
    Last edited: Aug 2, 2009
  7. Aug 2, 2009 #6

    CNX

    User Avatar

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

    Code (Text):

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

    ans =

         4.07430975749267e-015
     
    and the resulting wavenumbers,

    Code (Text):

               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.
     
  8. Aug 2, 2009 #7

    Redbelly98

    User Avatar
    Staff Emeritus
    Science Advisor
    Homework Helper

    Well, at least you're in the ballpark now. Can't help anymore than I have, sorry.
     
  9. Aug 2, 2009 #8

    CNX

    User Avatar

    Thanks for your help -- I'll keep working at it
     
    Last edited: Aug 3, 2009
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Molecular Vibrations - Numerical
  1. Sound vibrations (Replies: 1)

  2. Vibration - Dynamics (Replies: 3)

  3. Period of vibration (Replies: 3)

  4. Vibrating string (Replies: 4)

Loading...