SCF-Direct

[pwborthwick]

Hi Josh, As promised I've raised this as an issue. The problem is with the direct=True option. With H₂O using aug-cc-pvdz basis fails to converge in 100 iterations (tol = 1e-12). The final iterations are

     96 (-76.00384620636785+0j)        4.465234668405578e-10
     97 (-76.00384620636542+0j)        3.3355716725767773e-10
     98 (-76.00384620632069+0j)        1.5942003869849322e-10
     99 (-76.00384620632447+0j)        2.783263711926929e-10
     NOT CONVERGED

The energy is nearly converged, but your convergence criteria is based on the delta density and that has been around the 1e-10 mark since cycle 20. With direct=False the converged energy is -76.003846206545 in 21 cycles. So it's not that the direct=True is failing it's just converging slowly.

If you look at a non-diffuse basis like dz then with direct=True is converges in 21 cycles compared to 18 without direct, albeit much slower because of the extra computation of the eri's each cycle.

Taking the case of H₂ with aug-cc-pvdz with tol=1e-14 again

      96 (-1.128778275349057+0j)        1.5412132357630283e-12
      97 (-1.1287782753490616+0j)       1.3135489030800225e-12
      98 (-1.1287782753490487+0j)       1.6866085565244321e-12
      99 (-1.1287782753490434+0j)       2.5065950282967338e-12
      NOT CONVERGED

reduce to tol=1e-12 and it converges to -1.128778275349 in 9 iterations. So direct=True again works but converges more slowly especially with diffuse basis sets. Having said that, the diffuse sets happen also to be the biggest so maybe its not a problem with them at all? It would be interesting to see how a full implementation of the screen and density contraction affects the computation time. I notice your sole criterion for convergence (apart from maximum cycles check) is on the delta density and not the energy itself. This is H₂ in sto-3g

       E(SCF)    =  -1.116759307396 in 1 iterations
       Convergence:
       FPS-SPF  =  7.850462293418876e-17
       RMS(P)   =  2.00e-16
       dE(SCF)  =  6.75e-01

The energy is only converged to 7e-01, is this OK? That's it then, don't really think this is that much of issue after all as you can get direct to converge by upping the iterations, but it would be a lot quicker not to use direct at the moment. Still the code is useful as a pointer to the direct techniques. Peter

[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