Subject: Question about the CPHF

From: Ed Brothers <[email-masked]>
Date: Mon, 13 Jan 2003 11:44:25 -0500

I sent something similar to this question before Christmas, but I
thought I would try the CCL again as I am still a bit stuck.
    I am trying to implement the CPHF in order to calculate the second
derivative of the energy with respect to nuclear motion.  In addition to
the famous Pople article from 1979, I have been using Frisch,
Head-Gordon, and Pople, Chem. Phys. 141, 1990, p 189-196.  My difficulty
lies in forming the first derivative of the density matrix.
    I am currently conducting test calculations on an H2 molecule at the
equilibrium geometry.  Each atom has one s type basis function on it.
It would appear to me that the first derivative of P(1,1) and P(2,2)
(where P is the denisty matrix) should be equal, and finite difference
calculations seem to bear this out.  But I get a unphysical answer using
the CPHF, which means I must be using it incorrectly.  Here's where my
problem arises:

As there is only one occupied orbital, the elements of the density
matrix are simply:

(1)    P(mu,nu) = c(mu,1) c(nu,1)

Where c is the MO orbital coefficients, and mu and nu refer to atomic
basis functions.  Taking the derivative:

(2)    dP(mu,nu)/dx = dc(mu,1)/dx c(nu,1) + dc(nu,1)/dx c(mu,1)

Now according to the Pople paper from 1979:

(3)     dc(nu,1)/dx = (Sum over q=1, N) c(nu,q) u(q,1)

Where u are the values derived from the CPHF.  (This eq. 29 in the
paper).  Since q can equal either 1 or 2:

(4)    dc(nu,1)/dx =  c(nu,1) u(1,1)  + c(nu,2) u(2,1)

Substituting (4) into (2):

(5)    dP(mu,nu)/dx = [c(mu,1) u(1,1)  + c(mu,2) u(2,1)]c(nu,1) +
[c(nu,1) u(1,1)  + c(nu,2) u(2,1)] c(mu,1)

Now lets look at two specific derivative density matrix elements,
dP(1,1)/dx and dP(2,2)/dx, which should be equal.

(6a)   dP(1,1)/dx = [c(1,1) u(1,1)  + c(1,2) u(2,1)]c(1,1) +
[c(1,1)u(1,1)  + c(1,2) u(2,1)] c(1,1)
(6b)   dP(2,2)/dx = [c(2,1) u(1,1)  + c(2,2) u(2,1)]c(2,1) +
[c(2,1)u(1,1)  + c(2,2) u(2,1)] c(2,1)

Rearranging:

(7a)    dP(1,1)/dx = 2 c(1,1)^2 u(1,1) + 2 c(1,1) c(1,2) u(2,1)
(7b)    dP(2,2)/dx = 2 c(2,1)^2 u(1,1) + 2 c(2,2) c(2,1) u(2,1)

> From the the coefficients of the first MO, we know c(1,1) = c(2,1). From
the second MO, we know that c(1,2) = -c(2,2).  If we substitute this
into (7) we get:

(8a)   dP(1,1)/dx = 2 c(1,1)^2 u(1,1) + 2 c(1,1) c(1,2) u(2,1)
(8b)   dP(2,2)/dx = 2 c(1,1)^2 u(1,1) -  2 c(1,2) c(1,2) u(2,1)

Thus it would appear that dP(1,1)/dx .ne. dP(2,2)/dx.  This is
incorrect.

I am guessing my error is in the way I use (3), but I am unsure as to
what I am doing wrong.  Any help you can give me would be enormously
appreciated.

Thanks for your time!
Ed Brothers.
Merz Group.
PSU.