Subject: First derivative of density matrix using CPHF.

From: Ed Brothers <[email-masked]>
Date: Tue, 03 Dec 2002 14:23:26 -0500

I am trying to construct the first derivative of the density matrix
with respect to an atomic center's perturbation.  I am able to solve the

CPHF and I believe I am getting the right answers, but I am getting
stuck on forming the first derivative of the density matrix.
    For the test case I am using an H2 with one uncontracted Gaussian
with orbital exponent=3 on each atom.  The optimized geometry and MO
coefficients follow:

H          0.116439735      0.116439735      0.116439735
H         -0.116439735     -0.116439735     -0.116439735

c[ 1, 1]= -0.5937423229
c[ 2, 1]= -0.5937423229 = c[ 1, 1]
c[ 1, 2]= 0.9271342977
c[ 2, 2]= -0.9271342977  =  - c[ 1, 2]

Where c[i,j] is AO i’s contribution to MO j.

According to Pople, Krishnan, and Binkley, Int. J. Quant. Chem 13,
225-241 (1979), The first derivative of density matrix element u,v
(Called D(u,v) below) is:

dD(u,v)/dy = (Sum over r=1,2N) (Sum over S=1,2N) P’(r,s) c[u,r] c[v,s]

Where N is the total number of basis functions, and r,s are over MO’s.
If we assume the UHF case and an alpha density matrix, u and v are both
alpha functions and r and s are only over alpha MO’s.  Also, for this
basis set there are only two alpha MO’s.  Thus:

dD(u,v)/dy = (Sum over r=1,2) (Sum over S=1,2) P’(r,s) c[u,r] c[v,s]

dD(u,v)/dy = P’(1,1) c[u,1] c[v,1] + P’(1,2) c[u,1] c[v,2] + P’(2,1)
c[u,2] c[v,1] + P’(2,2) c[u,2] c[v,2]

My problem arises if we look at the derivative of density matrix
elements 1,1 and 2,2.  It would seem to me that for an arbitrary
Cartesian perturbation, dD(1,1)/dy = dD(2,2)/dy.  For instance, if this
molecule was along the x-axis and we took the derivative of Atom 1 in
the x direction, both D(1,1) and D(2,2) would change by the same
amount.  This doesn’t appear to be the case from my understanding of the

mathematics, which indicates I am making an error somewhere.  The
derivative are formally:

dD(1,1)/dy = P’(1,1) c[1,1] c[1,1] + P’(1,2) c[1,1] c[1,2] + P’(2,1)
c[1,2] c[1,1] + P’(2,2) c[1,2] c[1,2]

dD(2,2)/dy = P’(1,1) c[2,1] c[2,1] + P’(1,2) c[2,1] c[2,2] + P’(2,1)
c[2,2] c[2,1] + P’(2,2) c[2,2] c[2,2]

Now P’(1,1) is occupied-occupied, P’(2,1) and P’(1,2) are
virtual-occupied, and P’(2,2) is virtual-virtual.  According to equation

(54) of Pople et. al.,

P’(occ,occ) = - S’(occ,occ)

P’(occ,virt) =  u(occ:virt)  where the occ-virt is a single index and u
is the result of the CPHF equation.

P’(virt,virt) = 0.

Thus the derivatives reduce to:

dD(1,1)/dy = -S’(1,1) c[1,1] c[1,1] + 2u(2:1) c[1,1] c[1,2]
dD(2,2)/dy = -S’(1,1) c[2,1] c[2,1] + 2u(2:1) c[2,1] c[2,2]

(Note:  I assume that P’(1,2) =  P’(2,1).  Is this correct?)

We saw above that c[ 2, 1]= c[ 1, 1] and c[ 2, 2]=  - c[ 1, 2].  If we
substitute this into the relation ship for dD(2,2)/dy we get:

dD(1,1)/dy = -S’(1,1) c[1,1] c[1,1] + 2u(2:1) c[1,1] c[1,2]

dD(2,2)/dy = -S’(1,1) c[1,1] c[1,1] - 2u(2:1) c[1,1] c[1,2]

And these are very obviously not equal.    Does anyone know what I am
doing wrong?  Any assistance would be greatly appreciated, as I think
this is a fairly trivial problem and I am stumped.

 Ed Brothers.
 Merz Group.
 PSU.

PS  Sorry about the length of the post.