• Outreach

From: JC Gumbart (gumbart_at_ks.uiuc.edu)
Date: Sat Apr 13 2013 - 08:08:17 CDT

Read the source code in QM tool (see text below). More importantly though,
don¡¯t use Paratool, use fftk. Getting force constants directly from the
hessian is not recommended either.

################################################################

# Unit conversion of internal Hessian matrix. #

# The force constants in the Hessian are not mass weighted and #

# given in Hartree/bohr^2 or Hartree/rad^2. They will be #

# converted to kcal/(mol*Angstrom^2) or kcal/(mol*rad^2) #

# In CHARMM force field k/2 is listed which is why we have to #

# multiply by 0.5: #

# CHARMM term: V = k(q-q0)^2 #

# quadratic expansion term: V = 1/2*F(q-q0)^2 #

# with the second derivative of the potential energy with #

# respect to the coordinates F = d^2V/d(q-q0)^2. #

# Thus we have k=F/2. #

# The coupling constants in the off-diagonals between bonds #

# and angles are converted to kcal/(mol*Angstrom*rad). #

# The force constants are scaled down by a constant factor #

# that cures the systematic overestimation of frequencies by #

# HF calculations and to less extend by B3LYP. #

################################################################

From: owner-vmd-l_at_ks.uiuc.edu [mailto:owner-vmd-l_at_ks.uiuc.edu] On Behalf Of
???
Sent: Friday, April 12, 2013 4:40 PM
To: vmd-l_at_ks.uiuc.edu
Subject: vmd-l: VMD, Paratool

Hello All,

I would like to use paratool plugin for the analysis of force constant for
water molecule.

Firstly, I use Gaussian to do a frequency analysis. The input file was like
this:

%chk=input.chk
%nprocshared=8
%mem=1gb
#p MP2/6-31G** IOP(7/33=1) freq=intmodes

******

0 1
O 0.6735 0.97431 -0.821672
H 1.0774 0.25058 -1.30049
H 0.2246 1.4863 -1.49445

In the output file (.log), I found the the Hessian matrix was like this:

The second derivative matrix:
R1 R2 A1
R1 0.58069
R2 -0.01115 0.58051
A1 0.03420 0.03420 0.16835
As I think, 0.58069 and 0.59051 just represent the force constant for O-H
bonds in water. Considering the unit for each value in the above hessian
matrix is hartree*bohr^-2, the force constant for O-H bond should be about
0.58069*2240.9=1301.3 kcal/mol/A^2.

However, if I import the .log file using paratool in VMD, I found the force
constant for O-H bond is about 515 kcal/mol/A^2.

Why the result calculated from paratool is so different from the result
analyzed directly from hessian matrix?

I appreciate any help!