SCF-Direct about mcmurchie-davidson

[jjgoings]

I think the problem here is the incredibly tight convergence criteria. For some reason I set the SCF convergence tolerance and the 2e integral tolerance to be equal (1E-12). 1E-12 is reasonable for integral screening, but is way too tight for SCF density changes, or even energy changes. RMS(P) changes are more reasonably set as 1E-8. So, the tight SCF convergence criteria, coupled with the integral screening means that in many cases the SCF will never "converge", due to error propagation (in general shouldn't expect SCF to be more accurate than the integral sit depends on).

The solution for now -- and this is not to say the direct SCF could be improved further -- would be to split the integral and SCF convergence criteria, a la:

def RHF(self,doPrint=True,DIIS=True,direct=False,conver=1e-8,acc2e=1e-12):
        """Routine to compute the RHF energy for a closed shell molecule"""

I think these are reasonable values, and are the defaults in Gaussian, for example.

[jjgoings]

Now the input

from mmd.molecule import Molecule
from mmd.postscf import PostSCF

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

# init molecule and build integrals
mol = Molecule(geometry=water,basis='aug-cc-pvdz')

# do the SCF
print("conventional")
mol.RHF(direct=False)
print(" ")
print("direct")
mol.RHF(direct=True)

Results in:

conventional
E(SCF)    =  -76.003846206182 in 15 iterations
  Convergence:
    FPS-SPF  =  5.917280296468369e-09
    RMS(P)   =  6.09e-09
    dE(SCF)  =  1.57e-08
  Dipole X =  -0.00000000
  Dipole Y =  2.10397365
  Dipole Z =  0.00000000

direct

E(SCF)    =  -76.003846206240 in 15 iterations
  Convergence:
    FPS-SPF  =  5.912503216000997e-09
    RMS(P)   =  6.11e-09
    dE(SCF)  =  1.57e-08
  Dipole X =  0.00000000
  Dipole Y =  2.10397365
  Dipole Z =  0.00000000

But direct=True is still very very slow.

[pwborthwick]

That looks a lot better doesn't it? Direct will be a lot slower as the eri's are being recalculated on every cycle at moment. But the density contraction works ok and given a bigger molecule the Schwartz should help too. Thanks, informative as always.

PS in BOMD.py pretty sure you should be dividing by mass as need acceleration in velocity-verlet ie Force/mass.

[jjgoings]

@pwborthwick Thanks! I agree. I'll look at BOMD.py...I think you are right, but naively making the changes doesn't give me the vibrational period I'd expect, so there might be something else missing.

[pwborthwick]

Hi Josh,    Thanks.    It makes little difference for your example as mass of hydrogen is 1 so it doesn't much matter if you multiply or divide. A quick sanity check on the dimensions of the two lines in question suggests the last term in line 42 is force (potential / [L]) ie [M][L]/[TT] multiplied by TT (dt) = [M][L] which should equal left hand side [L] and velocity times dt ([L]/[T] x [T]) ie need to divide last term by [M]. Similarly line 52 last term needs to be a [L]/[T]. I've implemented a version of BOMD and get nice results for water. You can see a symmetric stretch and a scissoring with the energy minima occurring when the resultant forces are smallest which sounds qualitatively right. The frequency of the stretch is about twice that of the bend which again is about right. Further you can see the values of the bond lengths and angle at minimum energy agree very well with those obtained by geometry optimization (using 6-31g).       Peter Ps. As I was writing this I wondered if slowquant does BOMD? I checked and it does and he divides by mass. Sent from my Huawei phone

[pwborthwick]

@pwborthwick Thanks! I agree. I'll look at BOMD.py...I think you are right, but naively making the changes doesn't give me the vibrational period I'd expect, so there might be something else missing.

You are using atomic mass units? Multiply by 1822.8839 to bring to atomic units. Is this why you're not getting the the right quantitative results? See psi4numpty. Psi4numpty still doesn't agree even if you divide by amu2au factor in your program. They use the half-step version of velocity-verlet but that is ok.

I'm getting a problem with BOMD and H2 for dz basis. Diis is failing with singular matrix (iteration 34). I've also had this problem in my program... I've had to trap error switch off diis and return the last Fock matrix for the scf to continue without diis.
Best, Peter