We develop a formalism to directly evaluate the matrix of force constants within a Quantum Monte Carlo calculation. We utilize the matrix of force constants to accurately relax the positions of atoms in molecules and determine their vibrational modes, using a combination of variational and diffusion Monte Carlo. The computed bond lengths differ by less than 0.007 Å from the experimental results for all four tested molecules. For hydrogen and hydrogen chloride, we obtain fundamental vibrational frequencies within 0.1% of experimental results and ∼10 times more accurate than leading computational methods. For carbon dioxide and methane, the vibrational frequency obtained is on average within 1.1% of the experimental result, which is at least 3 times closer than results using restricted Hartree-Fock and density functional theory with a Perdew-Burke-Ernzerhof functional and comparable or better than density functional theory with a semi-empirical functional.