Hi all,
For a while now I have been building a small application to generate AMBER forcefield compatible force constants (stretch and angle). I use the full Force Constant matrix for a Gaussian OPT->FREQ calculation. I am not working on normal modes, I am using the vibrations projected on the bond through a displaced of one atom from the other:
# B3LYP/ OPT=(calcfc,verytight,noraman) int=ultrafine FREQ=noraman geom=connectivity
Testing on Methane (CH4) over a number levels of theory I sometime get unequal force constants for the C-H bonds, but only in pairs, e.g., C-H1 and C-H3 are always equal, and C-H2 and C-H4 are always equal, e.g.,
B3LYP/6-311++g**
101.12 101.12 101.12 101.12
B3LYP/6-311++g(2d,2p)
98.76 98.76 207.12 207.12
After looking at the two different frequency results I see that one is very slightly different to the other - see attached. Observing the actual vibrations in motion and I see whilst one calculation (where the force constants are all equal) demonstrates a simultaneous stretch vibration, the other (where the force constants are not equal) pins two of the hydrogens in place, keeping them still whilst the remaining two hydrogens vibrate.
How can I address this obvious error in the frequency calculations? I have checked that the wave function is stable, it is.
See attached for a better description.
First of all You can not use harmonic approximation, Anthony, because of two well nature of the potential of N's movement perpendicular to H's plane. If N is in the plane of H, then symmetry of system is other, then in ground state, therefore electronic state is changing. This all means non-adiabaticity. No software can account for this automatically.
Dear Anthony,
B3LYP/6-311++g** bases set has polarization and diffuse functions, whereas B3LYP/6-311++g(2d,2p) has only diffusion functions; this results in false values calculated by the latter basis.
Hoping this will be helpful,
Rafik
Dear Rafik,
I completely see your point. Thanks for the message, I will apply polarisation.
Many thanks
Dear Rafik,
I may be completely wrong (please correct me if I am), but according to the description in "Exploring Chemistry with Electronic Structure Methods" (the Gaussian book):
"...the 6-31G(2d) basis set adds two functions per heavy atom instead of just one, while the 6-311++G(3df,3dp) basis set contains three sets of valence region functions, diffuse function on both heavy atoms and hydrogens, and multiple polarization functions: 3 d functions and 1 f function on heavy atoms and 3 p functions and 1 d function of hydrogen atoms."
Therefore the theory and basis set that I used, B3LYP/6-311++g(2d,2p), does have polarization. Unless I am very mistaken.
I an OPT + FREQ calculation of B3LYP/6-311++g** on ammonia, expecting symmetrical vibrations.
Unfortunately the vibrations were not symmetrical, with one of the hydrogens remaining fixed (hydrogen label 3) whilst the other two vibrate (this is not then repeated around the molecule). Two frequencies showed symmetrical vibrations. This has biased the projected force constants onto the bonds, making one of the bond stretch force constants 'weaker' than the other two (which are identical).
Anthony: did you start the optimization with an exactly symmetrical (Td) molecule for methane? Remember, that geometry convergence in gaussian (and any QM program) is iterative and may finish at slightly different geometries depending on the starting point. If you don't start with a proper symmetry it's unlikely, that gaussian optimization will converge to a fully symmetrical system.
Try also Symmetry=Follow keyword.
Many thanks Bartosz, I will give this a try. And my apologies, I said methane when I meant to say ammonia (see the linked images above). However, I have done similar tests on both with similar results.
On your recommendation I have just ran calculations on ammonia with symmetry=follow. Before running the simulation I performed a 'clean' in GView, all three H-N-H angles were 109.47122. After performing B3LYP/6-311++g** opt+freq I am still not getting symmetrical vibrations i.e, one of the three sets of projected force constants on the N-H bond is not equal to the other two, and this is always for the case of the hydrogen that is essentially 'static' during one of the vibrations, which is then not repeated across the other two hydrogens for any of the other computed frequency vibrations.
Interestingly, the optimised angles H-N-H are now:
107.89147
107.89147
107.89150
Also if it makes any difference, Gaussian identifies the C3v symmetry ammonia should have
Hello Anthony,
Three thinks that you could try are:
1) Change the Pople basis to Dunning ones (Sometimes I observed noticeable differences in some results).
2)Try other functional. For example B3PW91 is one of the functional famous to give very good frequencies but I think that the type of movements will not change.
3)Also, you could try with a MP2 optimization. The results will not be too much good but the different types of frequencies should be located correctly.
I hope it helps you,
Joaquim.
Hello Joaquim, thanks for the post. I have tried with MP2 but had the same results. I will give your points, 1) and 2) a try.
Dear Eugene,
This is where I'm going to have to show my ignorance. Your comment led me on a unsuccessful six hour reading mission to figure out how non-adiabatic/adiabatic systems could quite relate to the variation in bond stretch vibrations of that system.
I would be most grateful if you would elaborate on your earlier comment.
Many thanks.
First of all You can not use harmonic approximation, Anthony, because of two well nature of the potential of N's movement perpendicular to H's plane. If N is in the plane of H, then symmetry of system is other, then in ground state, therefore electronic state is changing. This all means non-adiabaticity. No software can account for this automatically.
Many thanks Eugene. I am encroaching in an area of science I am very unfamiliar with. But I shall persevere.
It is non-rigid molecules, Anthony. Banker (or Bunker, don't remember) has written the book about this.
- Badges
- Science topic