Hi,

I have been inspecting some of the output of psi4 and found what may be bug. In my test I had a H2O trimer with O in (0, 0, 0) and the two H at (0.7cos(t), 0.7sin(t), -0.7) and (0.7cos(t+PI), 0.7sin(t+PI), -0.7). H-O-H angle is not the correct one for water most likely but it really does not matter. I was checking how the density matrix elements change when the trimer is rotated around the z axis, and computed the same molecule with t=0 -> Pi in steps of Pi/6.

All CCSD (cc-pvdz, symmetry C1) calculations returned the same energy. Then I looked at some DM elements and noticed a non-smooth behviour. I checked the much more intuitive overlap matrix elements and saw the same.

Here are the values of the overlap between the first p of O and the first s of one of the H:

0.189045, 0.189045, 0, 0, 0, -0.189045, -0.189045

corresponding to the 7 rotations. These really do not look right. All overlap elements between O and H that I inspected have the same step-wise trend, either resembling a sine or cosine.

I manually repeated the integrals calculation with my own libcint hooks and got the followings:

0.18905, 0.16372, 0.09452, 0, -0.09452, -0.16372, -0.18905

When I plotted these, they look like a nice cosine function that is expected.

psi4 matrices were obtained from the wavefunction object with:

wf.S().to_array(copy=False, dense=True)

from the RHF initial calculation.

Is this a bug or am I missing something?

For completeness here are the positions of the H atoms in the 7 runs (O is always in 0,0,0):

H1 0.7 0.0 -0.7

H2 -0.7 0.0 -0.7

H1 0.606217782649 0.35 -0.7

H2 -0.606217782649 -0.35 -0.7

H1 0.35 0.606217782649 -0.7

H2 -0.35 -0.606217782649 -0.7

H1 0.0 0.7 -0.7

H2 0.0 -0.7 -0.7

H1 -0.35 0.606217782649 -0.7

H2 0.35 -0.606217782649 -0.7

H1 -0.606217782649 0.35 -0.7

H2 0.606217782649 -0.35 -0.7

H1 -0.7 0.0 -0.7

H2 0.7 0.0 -0.7