search

jjgoings / McMurchie-Davidson

do a simple closed shell Hartree-Fock using McMurchie-Davidson to compute integrals
BSD 3-Clause "New" or "Revised" License
78 stars 17 forks source link

Angular Momentum #12

Closed pwborthwick closed 3 years ago

pwborthwick commented 3 years ago

Hi Josh, As requested this is all I know at the moment. ANGULAR MOMENTUM MMD

  water = """
  0 1
  O    0.000000000   -0.075791844    0.000000
  H    0.866811829    0.601435779    0.000000
  H   -0.866811829    0.601435779    0.000000
  """   
  mol = Molecule(geometry=water,basis='sto-3g')
  mol.RHF(direct=False)
  print(mol.L[0])

E(SCF) = -74.942079897739 in 9 iterations Convergence: FPS-SPF = 8.609961786824215e-11 RMS(P) = 5.71e-10 dE(SCF) = 1.57e-07 Dipole X = 0.00000000 Dipole Y = 1.53400920 Dipole Z = -0.00000000
0. 0. 0. 0. 0.25648831 0. 0.
0. 0. 0. 0. 0.16698098 0. 0.
0. 0. 0. 0. 0. 0. 0.
0. 0. 0. 0. -1. 0. -0.
-0.25648831 -0.16698098 0. 1. . 0.16778154 0.16778154
0. 0. 0. 0. -0.16778154 0. 0.
0. 0. 0. 0. -0.16778154 0. 0.

psi4 I'm new to psi4 so excuse the example, I've tried to stop psi4 reorienting the molecule to make comparisons valid but not sure how successful I've been

 h2o = psi4.geometry("""
 O    0.000000000   -0.075791844    0.000000
 H    0.866811829    0.601435779    0.000000
 H   -0.866811829    0.601435779    0.000000

 units angstrom
 symmetry c1
 noreorient
 nocom
 """)

 psi4.set_options({'basis': 'sto-3g',
              'scf_type': 'direct',
              'e_convergence': 1e-8,
               'print':5})

 e,wfn = psi4.energy('scf/sto-3g',puream=False,return_wfn = True)
 mints = psi4.core.MintsHelper(wfn.basisset())
 L = mints.ao_angular_momentum()
 print(e)
 print(np.asarray(L[0]))
-74.9420798970327
0. 0. 0.14352429 0. 0. 0. 0.
0. 0. 0.09343828 0. 0. 0. 0.
-0.14352429 -0.09343828 0. 0. 1. 0.18625536 0.18625536
0. 0. 0. 0. 0. 0. 0.
0. 0. -1. 0. 0. 0. 0.
0. 0. -0.18625536 0. 0. 0. 0.
0. 0. -0.18625536 0. 0. 0. 0.

I've checked L[1] and L[2] just in case the axes have been swapped. Psi4 swaps px and pz although I thought I told it not to.

I've looked at the code in angular and can't see why it shouldn't be right based on the original McMurchie-Davidson paper equations 2.34 which gives Ψ∇Ψ and 3.6 which gives [NLM|xc].

I've also tried the commutators in the molecular basis in both programs - and they still don't commute. (In theory H, L2, Lz, S2 and Sz should all commute).

NABLA As a check I coded a del (nabla) integral based the angular code as 

 def nabla(a, lmn1, A, b, lmn2, B, direction):
     l1,m1,n1 = lmn1
     l2,m2,n2 = lmn2
     P = gaussian_product_center(a,A,b,B)

     S0x =    E(l1,l2,0,A[0]-B[0],a,b) 
     S0y =    E(m1,m2,0,A[1]-B[1],a,b) 
     S0z =    E(n1,n2,0,A[2]-B[2],a,b)     

     D1x = l2*E(l1,l2-1,0,A[0]-B[0],a,b) - 2*b*E(l1,l2+1,0,A[0]-B[0],a,b)
     D1y = m2*E(m1,m2-1,0,A[1]-B[1],a,b) - 2*b*E(m1,m2+1,0,A[1]-B[1],a,b)
     D1z = n2*E(n1,n2-1,0,A[2]-B[2],a,b) - 2*b*E(n1,n2+1,0,A[2]-B[2],a,b)

     if direction.lower() == 'x':
         return D1x*S0y*S0z*np.power(pi/(a+b),1.5) 

     elif direction.lower() == 'y':
         return S0x*D1y*S0z*np.power(pi/(a+b),1.5) 

     elif direction.lower() == 'z':
         return S0x*S0y*D1z*np.power(pi/(a+b),1.5)

This agrees with psi4 ao_nabla() in all directions.

DIPOLE psi4 and mmd both give the same dipole magnitude

   mmd   Dipole Y =  1.53400920
   psi4  Y:     1.5340

The x-direction integrals (M[0]) agree, but M[1] and M[2] do not agree with ao_dipole() y and z components. I'm in the middle of programming a quadrupole and that again seems to agree in xx direction only at the moment.

Summary The dipole is a derived quantity from the integral matices and the dipole vectors agree in both programs ie both magnitude and orientation. However we know that the precursor matrices do not agree in all components in both programs. From this we might conclude that the matrices can differ and still be correct (in that programs internal scheme). mmd's code for the angular momentum looks correct (you have used Helgaker and I have used McMurchie-Davidson), but there is no way to derive a quantity by which to check (which is why the commutators would have been nice). The nabla agrees in both programs (so both are right or wrong) - why does it agree in all components while dipole doesn't? My feeling is that mmd is correct but then why does psi4 output things not in the original orientation. My 'feeling' is that the angular momentum is gauge invariant (in absence of external field) but I also need to check that. There is more work to do here!

postscript - quadrupole This is my quadrupole code (my parameter order differ a bit from mmd)

  def quadrupole(ia, ja, ie, je, ir, jr, kr, direction):
     # quadrupole moment
     p = ie + je
     q = ((ie*ir + je*jr)/p) - kr
     ijr = ir - jr

     sx = e(ia[0], ja[0], 0, ijr[0], ie, je)
     sy = e(ia[1], ja[1], 0, ijr[1], ie, je) 
     sz = e(ia[2], ja[2], 0, ijr[2], ie, je) 

     tx = e(ia[0], ja[0], 1, ijr[0], ie, je) + q[0]* e(ia[0], ja[0], 0, ijr[0], ie, je) 
     ty = e(ia[1], ja[1], 1, ijr[1], ie, je) + q[1]* e(ia[1], ja[1], 0, ijr[1], ie, je)
     tz = e(ia[2], ja[2], 1, ijr[2], ie, je) + q[2]* e(ia[2], ja[2], 0, ijr[2], ie, je)

     if direction == 'xx':
         u = 2.0 * e(ia[0], ja[0], 2, ijr[0], ie, je) + 2.0 * q[0]* e(ia[0], ja[0], 1, ijr[0], ie, je)  \
               + (q[0]*q[0] + (0.5 / p)) * e(ia[0], ja[0], 0, ijr[0], ie, je)
         return u * sy * sz * pow(pi/p, 1.5)
     if direction == 'yy':
         u = 2.0 * e(ia[1], ja[1], 2, ijr[1], ie, je) + 2.0 * q[1]* e(ia[1], ja[1], 1, ijr[1], ie, je)  \
               + (q[1]*q[1] + (0.5 / p)) * e(ia[1], ja[1], 0, ijr[1], ie, je) 
         return sx * u * sz * pow(pi/p, 1.5)
     if direction == 'zz':
         u = 2.0 * e(ia[2], ja[2], 2, ijr[2], ie, je) + 2.0 * q[2]* e(ia[2], ja[2], 1, ijr[2], ie, je)  \
               + (q[2]*q[2] + (0.5 / p)) * e(ia[2], ja[2], 0, ijr[2], ie, je) 
         return sx * sy * u * pow(pi/p, 1.5)
     if direction == 'xy':
         return  tx * ty * sz * pow(pi/p, 1.5)
     if direction == 'yz':
         return  sx * ty * tz * pow(pi/p, 1.5)
     if direction == 'zx':
         return  tx * sy * tz * pow(pi/p, 1.5)

The xx component agrees with psi4 (similar to dipole) but other components differ.

Finally, thanks for the attachment in your e-mail I'm always grateful for papers etc. Peter

jjgoings commented 3 years ago

Awesome, thank you for this work! I'll get to it when I have time, but please share any progress you make! I'm curious about the commutators, but assuming we can get it to work (or understand why it doesn't), it would be a great addition to the test suite.

jjgoings commented 3 years ago

@pwborthwick The discrepancy seems to arise because in mmd the dipole and angular routines allow for a gauge origin C different from the origin (0,0,0). I had set C to be equivalent to the center of positive charge throughout the program, but by setting self.center_of_charge to zero, we obtain the same integrals as psi4 and Gaussian16. I made this change years ago when I was very interested in gauge origin issues in molecular property calculations, but I'm not sure if this is the desired behavior, although for gauge origin independent properties it makes no difference (which is why computed dipole did not change). I like retaining the ability to choose the origin of the multipole moment integrals, but I am not sure the center of positive charge is the best idea. Thoughts?

For the time being though, here is how to reproduce the issue. Running:

import numpy as np
from mmd.molecule import Molecule

np.set_printoptions(precision=5)

water = """
  0 1
  O    0.000000000   -0.075791844    0.000000
  H    0.866811829    0.601435779    0.000000
  H   -0.866811829    0.601435779    0.000000
  """

mol = Molecule(geometry=water,basis='sto-3g')
mol.RHF()

print("\nCenter of charge before: ",mol.center_of_charge)

print("Y-dipole (length gauge)")
print(mol.M[1])

print("R x Del X")
print(mol.L[0])

# rebuild integrals with center of charge zeroed out
mol.center_of_charge *= 0
mol.build()

print()
print("Center of charge after: ",mol.center_of_charge)

print("Y-dipole (length gauge)")
print(mol.M[1])

print("R x Del X")
print(mol.L[0])

yields

E(SCF)    =  -74.942079897838 in 9 iterations
  Convergence:
    FPS-SPF  =  8.805776323613215e-11
    RMS(P)   =  5.20e-10
    dE(SCF)  =  1.57e-07
  Dipole X =  0.00000000
  Dipole Y =  1.53400920
  Dipole Z =  -0.00000000

Center of charge before:  [0.      0.11273 0.     ]
Y-dipole (length gauge)
[[-0.25595 -0.06059  0.       0.05079  0.      -0.00809 -0.00809]
 [-0.06059 -0.25595  0.       0.64117  0.       0.10644  0.10644]
 [ 0.       0.      -0.25595  0.       0.       0.07869 -0.07869]
 [ 0.05079  0.64117  0.      -0.25595  0.       0.31045  0.31045]
 [ 0.       0.       0.       0.      -0.25595  0.       0.     ]
 [-0.00809  0.10644  0.07869  0.31045  0.       1.02382  0.18609]
 [-0.00809  0.10644 -0.07869  0.31045  0.       0.18609  1.02382]]
R x Del X
[[-0.      -0.      -0.      -0.       0.25649 -0.      -0.     ]
 [ 0.      -0.      -0.      -0.       0.16698 -0.      -0.     ]
 [ 0.       0.      -0.      -0.      -0.      -0.      -0.     ]
 [ 0.       0.       0.      -0.      -1.      -0.      -0.     ]
 [-0.25649 -0.16698  0.       1.      -0.       0.16778  0.16778]
 [ 0.       0.       0.       0.      -0.16778 -0.      -0.     ]
 [ 0.       0.       0.       0.      -0.16778  0.      -0.     ]]

Center of charge after:  [0. 0. 0.]
Y-dipole (length gauge)
[[-0.14323 -0.0339   0.       0.05079  0.      -0.00376 -0.00376]
 [-0.0339  -0.14323  0.       0.64117  0.       0.14997  0.14997]
 [ 0.       0.      -0.14323  0.       0.       0.10895 -0.10895]
 [ 0.05079  0.64117  0.      -0.14323  0.       0.33409  0.33409]
 [ 0.       0.       0.       0.      -0.14323  0.       0.     ]
 [-0.00376  0.14997  0.10895  0.33409  0.       1.13655  0.20658]
 [-0.00376  0.14997 -0.10895  0.33409  0.       0.20658  1.13655]]
R x Del X
[[-0.      -0.      -0.      -0.       0.14352 -0.      -0.     ]
 [ 0.      -0.      -0.      -0.       0.09344 -0.      -0.     ]
 [ 0.       0.      -0.      -0.      -0.      -0.      -0.     ]
 [ 0.       0.       0.      -0.      -1.      -0.      -0.     ]
 [-0.14352 -0.09344  0.       1.      -0.       0.18626  0.18626]
 [ 0.       0.       0.       0.      -0.18626 -0.      -0.     ]
 [ 0.       0.       0.       0.      -0.18626  0.      -0.     ]]
jjgoings commented 3 years ago

Also, as I think of it, I think the molecular wave function will not be an eigenstate of the angular momentum operator because it breaks spherical symmetry. So the angular momentum commutation relations should not be expected to hold for most systems of interest except for the case of atoms. Kind of like how once we add in SO coupling we break the spin symmetry. See these notes, for example.

jjgoings commented 3 years ago

Looking more closely, I believe behavior in mmd is correct, but the handling is a bit more subtle. By default Gaussian reorients the molecules to be centered at the center of nuclear charge, thereby setting the multipole origin to (0,0,0). Turning off symmetry removes this behavior. I am not sure Psi4 reorients properly. mmd takes into consideration the center of nuclear charge explicitly at the integral level, without requiring a reorientation. It might be worth exposing the C multipole origin argument so that the user can set a custom origin, if desired. Apparently MOLCAS and MOLPRO allow for this type of functionality.

So whereas mmd gives for the above example:

Center of charge before:  [0.      0.11273 0.     ]
Y-dipole (length gauge)
[[-0.25595 -0.06059  0.       0.05079  0.      -0.00809 -0.00809]
 [-0.06059 -0.25595  0.       0.64117  0.       0.10644  0.10644]
 [ 0.       0.      -0.25595  0.       0.       0.07869 -0.07869]
 [ 0.05079  0.64117  0.      -0.25595  0.       0.31045  0.31045]
 [ 0.       0.       0.       0.      -0.25595  0.       0.     ]
 [-0.00809  0.10644  0.07869  0.31045  0.       1.02382  0.18609]
 [-0.00809  0.10644 -0.07869  0.31045  0.       0.18609  1.02382]]
R x Del X
[[-0.      -0.      -0.      -0.       0.25649 -0.      -0.     ]
 [ 0.      -0.      -0.      -0.       0.16698 -0.      -0.     ]
 [ 0.       0.      -0.      -0.      -0.      -0.      -0.     ]
 [ 0.       0.       0.      -0.      -1.      -0.      -0.     ]
 [-0.25649 -0.16698  0.       1.      -0.       0.16778  0.16778]
 [ 0.       0.       0.       0.      -0.16778 -0.      -0.     ]
 [ 0.       0.       0.       0.      -0.16778  0.      -0.     ]]

Gaussian16 yields (note it aligns to z-axis):

 Multipole matrices IBuc=  518 IX=    3 IJ=           1:
                1             2             3             4             5
      1  0.255955D+00
      2  0.605855D-01  0.255955D+00
      3  0.000000D+00  0.000000D+00  0.255955D+00
      4  0.000000D+00  0.000000D+00  0.000000D+00  0.255955D+00
      5  0.507919D-01  0.641173D+00  0.000000D+00  0.000000D+00  0.255955D+00
      6  0.808810D-02 -0.106443D+00  0.000000D+00 -0.786914D-01  0.310448D+00
      7  0.808810D-02 -0.106443D+00  0.000000D+00  0.786914D-01  0.310448D+00
                6             7
      6 -0.102382D+01
      7 -0.186089D+00 -0.102382D+01

 R x Del for IMat=    2:
                1             2             3             4             5
      1  0.000000D+00
      2  0.000000D+00  0.000000D+00
      3 -0.256488D+00 -0.166981D+00  0.000000D+00
      4  0.000000D+00  0.000000D+00  0.000000D+00  0.000000D+00
      5  0.000000D+00  0.000000D+00  0.100000D+01  0.000000D+00  0.000000D+00
      6  0.000000D+00  0.000000D+00 -0.167782D+00  0.000000D+00  0.000000D+00
      7  0.000000D+00  0.000000D+00 -0.167782D+00  0.000000D+00  0.000000D+00
                6             7
      6  0.000000D+00
      7  0.000000D+00  0.000000D+00

Which are equal up to a sign or permutation, due to the extra "reorienting" Gaussian does.

pwborthwick commented 3 years ago

Hi @jjgoings , OK firstly, thanks for putting your time into these issues. As I wrote I had suspected that both programs were (in their way) correct. I couldn't see anything wrong with the mmd code and psi4 should be well tested. Ironically I looked at psi4 code and in the quadrupole code found 

    Vector3 geom = mol->xyz(i) - origin;
    ret[0] += mol->Z(i) * geom[0] * geom[0];  // xx

I immediately thought 'ah they are subtracting the charge center' just like mmd. It's only now you've pointed it out that I find..

        std::string str = options["PROPERTIES_ORIGIN"][0].to_string();
        if (str == "COM") {
            for (int atom = 0; atom < natoms; ++atom) property[atom] = mol->mass(atom);
        } else if (str == "NUCLEAR_CHARGE") {
            for (int atom = 0; atom < natoms; ++atom) property[atom] = mol->charge(atom);
        } else {
            throw PSIEXCEPTION("Invalid specification of PROPERTIES_ORIGIN.  Please consult the manual.");
        }

So psi4 allows for these 3 gauge origins as you said. The cartesian origin is the one most easily conceptualised, if you've set the geometry you'll know where your origin is. However, it lacks any physical/chemical significance. The COM and nuclear charge centres (sorry I'm English) have a physical significance but perhaps lack an immediately identified location. Personally I like the charge centre, it seems more 'chemical' :) I certainly think there should be an option to choose the gauge center, but what default behaviour? From what you said in your earlier post Gaussian and psi4 adopt the cartesian origin so it might be as well to adopt it if it's an 'industry standard'? I will check the integrals in the next few days and try to get quadrupole working. This also explains why nabla was OK.

I read the article re angular momentum, but will go back to my notes to convince myself as I feel there's more to be learned here.

Anyway thanks as always for your efforts and time. Best Peter

pwborthwick commented 3 years ago

Hi @jjgoings, Just some observations on the angular momentum. Checking the spherically symmetric case...

o = """
0 1
O    0.0 0.0 0.0
"""

mol = Molecule(geometry=o, basis = 'sto-3g')
mol.RHF(direct=False , DIIS = False)

lox = np.dot(mol.X.T, np.dot(mol.L[0], mol.X))
loy = np.dot(mol.X.T, np.dot(mol.L[1], mol.X))
loz = np.dot(mol.X.T, np.dot(mol.L[2], mol.X))
lo2 = np.dot(lox,lox) + np.dot(loy,loy) + np.dot(loz,loz)

#check [Ly,Lz] = Lz
comm = np.dot(loy,loz) - np.dot(loz,loy)
print(np.allclose(comm, lox))

#check [LL, Lz] = 0
zero = np.zeros((5,5))
comm = np.dot(lo2, loz) - np.dot(loz,lo2)
print(np.allclose(comm, zero))

results in...

True
True

(orthogonalisation makes no difference in this case). In the spherically symmetric case L's and L2 are gauge invariant. Lx2+Ly2+Lz2 = L2 in this case acts like x2+y2+z2 = ds2 in Euclidean space ie it's a metric in angular momentum space and as such has the same sort of invariance properties as distance between points in Euclidean space. If you deform the flatness of the space then ds2 > 0 or < 0, destroying the wavefunction symmetry is analogous to deforming the space and L2 invariance is destroyed. One could conjecture that [L2,Lz] is a measure of the divergence of the system from symmetry. For example, for water || [L2, Li || is 0.20, but for methane (a fairly symmetric molecule) the value is 1e-15! In fact, to all intents and purposes, the (L2) commutator relation is satisfied. Note L2 matrix is Hermitian.

jjgoings commented 3 years ago

Came across this in Barron's "Molecular Light Scattering and Optical Activity" book, and seemed worth mentioning: Screen Shot 2020-12-29 at 2 18 39 PM Testing this out in Gaussian16, for neutral systems disabling symmetry/reorientation makes no difference for the computed dipole moment. For charged systems disabling reorientation can make a significant difference (I put the +2 charged water molecules about 1000A away from Cartesian origin --- BIG changes). Seems like center of nuclear charge is a good way to enforce consistency in results. But I'm not sure if center of nuclear charge will be identical to center of charge overall?