I have written a short program which calculates the force constant of bond stretching and bond angle using a weighted cartesian coordinate hessian matrix from the output of Gaussian (opt and freq).
I am using a water molecule as my first test case. With my code I can extract the bond stretch force constant for O-->H1, H1-->O, O-->H2, and H2-->O. In theory, everything should be identical (and the same for a molecule like CH4).
Unfortunately, I am at a stage where I am not sure whether my maths in my code is incorrect, or am I observing an actually effect of the vibrational frequencies. I have test HF, MP2 and B3LYP, with a varying degree in basis sets. I have attached a .pdf of the force constants assigned to each bond according to the type of QM experiment used.
At a small basis set, I see a discrepancy between O-->H1 and H1--> (which should be the same). I also see a discrepancy between both sets of force constants shared across O-->H1 with those shared across O-->H2, when those should be the same to. Any thoughts why I get large inconsistencies?
the values in your attached file depend on your program, you are the only one who knows what is going on there. but, you need to make sure that the calculated frequencies are accurate.
vibrational frequencies depend on your choice of functional/basis-set combinations (2nd derivative of the potential energy surface), the basis sets dependency becomes less significant for larger basis sets, you need to test these combinations until you get something accurate experimentally (water and methan are well established). you need to correct your calculated frequencies, you can use a scaling factor (depends on your combination) from NIST website, they have a large set of scaling factors for several functional/basis combinations. you can use the SCALE keyword in your input file. once you get accurate frequencies, your force constants depend completely on your program.
Many thanks for all of that information. I wasn't aware of the SCALE keyword. I'm now quite confident that my maths is correct, as you said, I just need to work on the functional/basis-set combination. Thanks again.
Mr. Nash,
Some of your data really look strange, because of for example you have 4 different force constants associated with two bond types, accounting for the fact that there is only one tensor component of the dipole moment according my computations. In paralle, only changing the method you have obtained four equal values for the same system, in fact. Most probably there is a problem with (i) the optimization method of the geometry, because of you have mentioned QM, but which one you have reflected in mind? Is those results in gas - phase of an isolated molecule H2O or ensemble of two- or more molecules? Or/and (ii) the path reading the geometry parameters!?
Despite those questions/comments, please find as attachment few computations using your theoretical levels and only two basis sets (6-31G and 6-31+G) (attachment). There are shown force constants and frequencies allowing with an experimental IR of H2O, but not matrix-isolated spectrum). I hope that those data shall help you to improve your algorithm/s. I shall try to compute the frequencies, soon, using other shown basis sets also, and shall post the data here. Most probably you shall find any correlation on which base to determine the scaling factor of the data using your program.
Hi Bojidarka,
These are in gas phase, with just one H2O.
The reason for two force constants over a single bond is because these have been calculated according to the second order tensor of atom A 'feeling' the force from the displacement of atom B. Which in turn, you could measure the other way i.e., atom B being displaced by atom A. However, in theory once you have converted the hessian into a mass weighted hessian matrix the force constant for A-B and B-A should be identical.
In some cases they are, but in other cases they are not. For example, using MP2/6-311++g**, O->H1 is 90 kcal/(mol.ang^2) whilst measuring the displacement of O relative to H1 i.e., H1->O the weighted hessian returns 93 kcal/(mol.ang^2). However, the opposite bond in the water (O-H2) returns 90 kcal/(mol.ang^2) if either the oxygen is displaced relative to the hydrogen or the hydrogren is displaced relative to the oxygen.
Thanks
Mr. Nash: Before you go much further, you may want to read the Gaussian whitepaper on "Vibrational Analysis" available from Gaussian.com website.
Hi Mr Ramachandra, I've read the document. It goes into detail about constructed a mass weighted hessian, that's easy enough. Now to reconstruct individual force constants for bond distances and bond angles (as per mechanical force fields e.g., Amber), you can use the eigenvectors from the subset of two atoms, A and B, from the hessian. This has been performed a number of times for parameterising Zinc, in zinc binding centres of proteins. However, when testing this process with a number of molecules, I tend to find that I have a bond stretch force constant of 70 for two of the C-H bonds in CH4, and a bond stretch force constant of 140 kcal/mol.ang^2 for the other two C-H bonds. This inconsistency must come from the values stored in the Hessian matrix as produced by Gaussian. Therefore, it must be an artefact from the level of theory and basis set used.
Thanks
Thank you for this additional information. So, you need local mode force constants rather than normal modes. Yes, I agree with your deduction of what must be going on.
Sergey, this is no longer a 'total mess' as you put it, just efforts in trying to understanding something from first principles. And yes, I understand that a force constant does not have a direction, however you should still be able to calculate a given bond stretch force constant between either A-B or B-A and yield the same result. Either way, this problem was fixed using the implementation of Jorge M. Seminario (Calculation of Intermolecular Force Fields from Second-Derivative Tensors).
The only issue that remains is what level of theory and basis set type should be used. I am working my way through several test cases to ensure appropriate results, but I find inconsistencies depending on the QM chemistry used. For example, in CH4. If I use B3LYP/6-311++g**, the force constant on the C-H bonds are not the same. Yet when I move to MP2/aug-cc-pvtz, all four calculated force constants are identical. Unfortunately MP2/aug-cc-pvtz for systems of around 70 to 80 atoms is proving too computationall expensive.
Many thanks Sergey. I will give a change of the integral a try.
Thanks again.
I have tried B3LYP/6-31+g(d,p) with Integral(Ultrafinegrid) on a large molecule (a crosslinked lysine dimer). Unfortunately, each C-H force constant on the methly branches (on the crosslink) yield unique force constants e.g., 101 kcal/(mol A^2), 184 (mol A^2) and 330 kcal/(mol A^2). I have identical force constants on smaller molecules (some test alkanes) using MP2/aug-cc-PVQZ, but I fear that that level of chemistry is too much for my molecules of interest.
- Badges
- Science topic
- Similar topics
- Filing