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

BOMD #7

Closed pwborthwick closed 3 years ago

pwborthwick commented 3 years ago

Hi Josh, I'll open a seperate issue for the BOMD problem as I think it's important to get it sorted as it affects more than just your program. I think we're both agreed that the correct equation should have the force divided by the mass in lines 42 and 53. I think the reason you were not then getting the behaviour you expected is because the mass used is '(unified) atomic mass units' - these are badly named! They are not atomic (Hartree) units being a relative unit ( relative to Carbon 12?). The conversion is 1 amu = 1822.8839 me. Changing the values in BOMD.py should then agree with the psi4numpy/MD-Verlet-Integrator implementation (with the same time step, basis and geometry) - except it doesn't. Psi4 gives a wave time of 175 atu and mmd 125 atu. Psi4 uses the half-step version of velocity-verlet but this should make no difference. Anyway after hours of working through the program and equations I noticed this line (ln74 md_helper.py)...

       vel_new += 0*5*timestep*accel_new

correcting to 

       vel_new += 0.5*timestep*accel_new

makes mmd and psi4numpty agree. As the whether the frequency is correct (or even in the right ball-park) is another question for another day. There are now 3 programs that do BOMD on github that I know about mmd, psi4numpty and slowQuant and all currently giving different results for different reasons - it would be nice to get a concensus?

Best wishes, Peter

PS there is a seperate issue with H2 sometimes causing solve in diis to encounter a singular matrix - should I open an issue for that?

jjgoings commented 3 years ago

@pwborthwick Thanks for bringing this up! I forgot that I wasn't treating mass as atomic units (we don't really use mass anywhere else in the code at the moment). I think that it is best to internally convert everything to atomic units anyway, so I went ahead and did that. Then switching to division by mass yields the correct behavior.

The other problem was obtaining the vibrational frequency from energy vs time instead of bond length vs time. During one vibrational period, the energy reaches its maximum twice due to conservation of energy, so the period appears to be twice as fast. The correct vibrational period can be determined visually from bond length vs time. The example has been updated. Currently we get a vibrational period of ~ 6.08 fs, corresponding to a vibrational frequency of 5486 cm-1. Gaussian16 results for the same system (H2/STO-3G) yield -- in addition to the identical SCF energy -- a frequency of 5481 cm-1. So this verifies the implementation for this simple case.

bomd

Feel free to continue to open cases where the behavior is not as you'd expect. I suspect the DIIS is a fringe issue due to the small size of H2/STO-3G (it's solution is analytic) but it's good to have a record of problems as they arise.

pwborthwick commented 3 years ago

@jjgoings Thanks. I agree that the diis issue is because of the size of molecule/basis and I won't raise an issue for that. I see you've alerted psi4numpy, did you include slowQuant owner as he uses amu rather au's? Thanks Josh this has been an informative exercise for me.