Coder Social home page Coder Social logo

jjgoings / mcmurchie-davidson Goto Github PK

View Code? Open in Web Editor NEW
75.0 3.0 17.0 1.79 MB

do a simple closed shell Hartree-Fock using McMurchie-Davidson to compute integrals

License: BSD 3-Clause "New" or "Revised" License

Python 84.95% Cython 15.05%
mcmurchie-davidson hartree-fock scf integrals quantum-chemistry electronic-structure chemistry

mcmurchie-davidson's Introduction

Build Status codecov Total alerts Language grade: Python

McMurchie-Davidson

This contains some simple routines to compute one and two electron integrals necessary for Hartree Fock calculations using the McMurchie-Davidson algorithm. The Hartree-Fock code can only handle closed shell molecules at the moment. It's not fast (though getting faster with Cython interface), but should be somewhat readable.

I'm slowly porting the integrals over to Cython and reoptimizing. I'm also thinking about reorganizing so as to make it easy to add functionality in the future.

Installation

Installation should be simple. For a local install, just

python setup.py build_ext --inplace install --user

Dependencies

You'll need numpy, scipy, and cython (for the integrals). Some post-SCF functionality requires bitstring, in order to handle the bitstring/integer representations for Slater determinants and to allow for the evaluation of Slater-Condon rules between arbitrary determinants. The install script should yell at you if you don't have the requisite dependencies. You can install them all at once if you have pip:

pip install numpy scipy cython bitstring

Testing

You can test the install with nosetests. In the head directory, just do

nosetests tests

it should take a few seconds, but everything should pass.

Running

Once you've installed, you can try running the input script sample-input.py:

python sample-input.py

which should do an SCF on water with an STO-3G basis and dump out to your terminal:

E(SCF)    =  -74.942079928060 in 10 iterations
  Convergence:
    FPS-SPF  =  3.377372360032819e-13
    RMS(P)   =  2.35e-12
    dE(SCF)  =  -9.95e-10
  Dipole X =  0.00000000
  Dipole Y =  1.53400931
  Dipole Z =  -0.00000000
E(MP2) =  -74.99122956422062

Input file specification

The input is fairly straightforward. Here is a minimal example using water.

from mmd.molecule import *
from mmd.scf import *

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='sto-3g')

# do the SCF
mol.RHF()

The first lines import the molecule and scf modules, which we need to specify our molecule and the SCF routines. The molecular geometry follows afterward and is specified by the stuff in triple quotes. The first line is charge and multiplicity, followed by each atom and its Cartesian coordinates (in Angstrom).

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

Then we generate create the molecule (Molecule object) and build the integrals, and finish by running the SCF.

At any point you can inspect the molecule. For example, you can dump out the (full) integral arrays:

print(mol.S)

which dumps out the overlap matrix:

[[ 1.     0.237  0.     0.     0.     0.038  0.038]
 [ 0.237  1.     0.     0.     0.     0.386  0.386]
 [ 0.     0.     1.     0.     0.     0.268 -0.268]
 [ 0.     0.     0.     1.     0.     0.21   0.21 ]
 [ 0.     0.     0.     0.     1.     0.     0.   ]
 [ 0.038  0.386  0.268  0.21   0.     1.     0.182]
 [ 0.038  0.386 -0.268  0.21   0.     0.182  1.   ]]

There is also some limited post-SCF functionality, hopefully more useful as a template for adding later post-SCF methods.

# do MP2
mp2 = PostSCF(mol)
mp2.MP2()

Examples

In the examples folder you can find some different scripts for setting up more complex jobs. For example, there is a simple script that does Born-Oppenheimer molecular dynamics on minimal basis hydrogen, aptly titled bomd.py. There are some real-time electronic dynamics in real-time.py. It's a good idea to check out the tests folder. In particular, tests/README.md contains an index to the different tests used, which gives some insight into the current functionality.

mcmurchie-davidson's People

Contributors

jjgoings avatar pwborthwick avatar

Stargazers

 avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar

Watchers

 avatar  avatar  avatar

mcmurchie-davidson's Issues

Tests fail on Apple M1

Compiling as usual, no apparent errors, but the tests fail miserably on Apple M1 with Python 3.9.

SCF does not converge or converges to very wrong answers, e.g. the sample-input.py case of STO-3G water:

E(SCF)    =  -92.762209945728 in 36 iterations
  Convergence:
    FPS-SPF  =  5.6769413580857e-09
    RMS(P)   =  5.74e-09
    dE(SCF)  =  1.82e-10
  Dipole X =  0.00000000
  Dipole Y =  -10.13611489
  Dipole Z =  -0.00000000
E(MP2) =  -92.81887498561838

Bug in integral code

Hi,

I tried to use you code (used integrals.py) to run a calculation on He with the basisset 6-31G, and got the following overlap matrix:

3.51E-01 | 3.75E-01
3.75E-01 | 1.00E+00

I think the diagonal, i.e. the selfoverlap should always be equal to 1.0 ?

I found this problem because I have in my own code (largely based on your code). I have not been able to figure it out. But have a suspicion that it have to do with the normalization of contracted basissets.

Please if you find a solution please send me an email: [email protected]

Changing source code

Hi,
I rewrote your integral engine in C++ for using it in my own code. However, I have some bugs regarding the nuclear attraction integrals in my code and I would like to add some print()-statements in your source code to see e.g. the values of nuclear attraction integrals between primitive gaussians.
If I change the code in mmd/integrals/reference.py and rebuild with python setup.py build_ext --inplace install --user the changes are not working. I suggest that this has to do with the already existing files in /cython. I am very inexperienced with Python, could you please tell me 1) whether I am allowed to change the source code at all and 2) how to solve this cython issue?
Thanks a lot!

FCI Eigenvectors

Hi Josh,
The eigenvectors produced by the FCI class do not seem to give correct results when used in eg transition dipole computations. For water in STO-3G using Crawfords geometry - I'm using this geometry as in his projects he gives reference Hamiltonians. Using CIS the first few printed energies and oscillator strengths are

Configuration Interaction Singles (CIS)
------------------------------
# Determinants:  40
nOcc * nVirt:    40
CIS state  1 (eV):       7.8166 (f=0.0000)
CIS state  2 (eV):       7.8166 (f=0.0000)
CIS state  3 (eV):       7.8166 (f=0.0000)
CIS state  4 (eV):       9.3723 (f=0.0000)
CIS state  5 (eV):       9.3723 (f=0.0000)
CIS state  6 (eV):       9.3723 (f=0.0000)
CIS state  7 (eV):       9.6998 (f=0.0023)

Using bitstrings via the commented out code (and making small modifications to the transition dipole computation line) you get

Configuration Interaction Singles (CIS)
------------------------------
# Determinants:  40
nOcc * nVirt:    40
CIS state  1 (eV):       7.8166 (f=0.0000)
CIS state  2 (eV):       7.8166 (f=0.0019)
CIS state  3 (eV):       7.8166 (f=0.0000)
CIS state  4 (eV):       9.3723 (f=0.2856)
CIS state  5 (eV):       9.3723 (f=0.0000)
CIS state  6 (eV):       9.3723 (f=0.0000)
CIS state  7 (eV):       9.6998 (f=0.0000)

The Hamiltonian in the bitstring case is substantially different from CIS because in CIS the 1st row is ordered (if $^a_i$ represents excitation from occupied orbital i to virtual orbital a and we have n occupied and u virtual spin-orbitals)
$^1_1$ $^2_1$ ... $^u_1$, $^1_2$ $^2_2$... $^u_2$,..., $^1_n$ $^2_n$ ... $^u_n$
whereas itertools gives an order of
$^1_n$ $^2_n$ ... $^u_n$, $^{n-1}_1$ $^{n-1}_2$ ... $^{n-1}_2$ ,..., $^1_1$ $^2_1$ ... $^u_1$
To get the bitstring Hamiltonian to the same form as the CIS one (which is reproduced in Crawford here you can reverse the whole determinant list and then reverse each 'number of virtual orbitals' sub-block. This leads to a Hamiltonian identical in absolute values to the CIS one but with some sign changes (the dominant diagonal is all positive - so this is not just arbitary phase of wavefunction). For example on the first row sixth column the bitstring Hamiltonian has -2.61966629e-02 whereas CIS has 2.61966629e-02 (and Crawford 0.0261967).
I'm wondering now if the FCI energies are correct? Given the strong diagonal dominance of the sparse Hamiltonian some sign errors in off-diagonal terms may not show up except in higher precision than we're considering here.
It would be good to get this resolved because not having the eigenvectors greatly reduces the usefulness of the bitstring FCI method.

I've used

det_list = list(zip(*(iter(det_list[::-1]),) * 4))
det_list = sum([list(i)[::-1] for i in det_list], [])
for i in det_list:
     print('{:6}     {:20}'.format(i, bin(i)[2:].zfill(nOrb)))

to make to order of the det_list same as CIS ordering. You might have to do some extra processing because your det_list is a list of arrays. What is reason for the array as everything seems to work fine without it?
Peter

row col        det A            det B        excitation     degree    phase               value
37   34   00011111111110    00101111111101    [[ 2  2]        2         1         -0.02619666286190461
                                               [ 1  0]
                                               [10 11]] 

Clearly the there are 2 new holes and 2 new particles [2, 2], occuring at 1 and 10 and 0 and 11. The level of excitation is clearly 2 and there is 1 permutation needed to move 1->0 and one needed to move 10->11 with no new holes between the excitations => $-1^{2}= +1$ . Can't see anything wrong there.

Handling of complex arithmetic

It would be good to toggle on or off real vs complex arithmetic.

Right now we (for the most part) treat everything as complex-valued, because the goal at the time was to allow for real-time electronic dynamics. For most users this is just additional computational overhead, so it would be good to implement the ability to switch between real and complex, and include safeguards to force complex if doing RT-TDSCF.

BOMD

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?

H2O/6-31++G** SCF not convering

Hi,

I tried run your gode with H2O/6-31++G**, but running:

from mmd.molecule import * 
from mmd.postscf import * 

water = """
0 1
H	0.866811829	0.601435778	0.0
O	0.000000000	-0.075791844	0.0
H	-0.866811829	0.601435778	0.0

"""

# init molecule and build integrals
mol = Molecule(geometry=water,basis='6-31ppgss')
mol.build()

# do the SCF
mol.RHF()

I took the basisset from the PSI4 github, since it uses the same format as your program. As an output I got:

NOT CONVERGED

I think it might be because your SCF solver dosent have DIIS? Atleast I couldnt see it when I looked in the script. I tried the same basis with my program with and without DIIS, see:

http://slowquant.readthedocs.io/en/latest/illustrativecalc.html

and it seems to yield no convergence too, when I turn off DIIS.

Question: Does the ERI routine require that E and R routines both be computed using Hermite Integrals?

I am trying to write a small HF code that can perform energy calculations on molecules (inspired by your excellent blog and this repository). However, I am presently calculating the overlaps via iteration (as explained in https://www.mathematica-journal.com/2012/02/16/evaluation-of-gaussian-molecular-integrals/). I found the other integral schemes (HGP and OS) to be pretty confusing, while MD algorithm is relatively straightforward. My query is regarding the calculation of ERI integrals, where both overlaps and repulsion integrals are calculated. Can I use the unnormalized overlaps calculated from iteration in conjunction with the recursive calculation of RI?

Thanks
Hemanth

Is the setup working?

Hi,

Please let me know if the current code setup is working.

I tried to install mmd and tried to run a simple script from the examples, but it throws me an error,

python sample-input.py
/home/varadarajan/softwares/mmd_overlap/McMurchie-Davidson/mmd/molecule.py:71: RuntimeWarning: divide by zero encountered in scalar divide
self.bfs.append(Basis(np.asarray(atom.origin),
/home/varadarajan/softwares/mmd_overlap/McMurchie-Davidson/mmd/molecule.py:71: RuntimeWarning: invalid value encountered in scalar multiply
self.bfs.append(Basis(np.asarray(atom.origin),
/home/varadarajan/softwares/mmd_overlap/McMurchie-Davidson/mmd/molecule.py:71: RuntimeWarning: invalid value encountered in scalar add
self.bfs.append(Basis(np.asarray(atom.origin),
Traceback (most recent call last):
File "/home/varadarajan/softwares/mmd_overlap/McMurchie-Davidson/sample-input.py", line 15, in
mol.RHF()
File "/home/varadarajan/softwares/mmd_overlap/McMurchie-Davidson/mmd/scf.py", line 21, in RHF
self.build(self.direct) # build integrals
File "/home/varadarajan/softwares/mmd_overlap/McMurchie-Davidson/mmd/molecule.py", line 92, in build
self.one_electron_integrals()
File "/home/varadarajan/softwares/mmd_overlap/McMurchie-Davidson/mmd/molecule.py", line 285, in one_electron_integrals
self.X = mat_pow(self.S,-0.5)
File "/home/varadarajan/softwares/my-jupyter-env/lib/python3.10/site-packages/scipy/linalg/_matfuncs.py", line 143, in fractional_matrix_power
return scipy.linalg._matfuncs_inv_ssq._fractional_matrix_power(A, t)
File "/home/varadarajan/softwares/my-jupyter-env/lib/python3.10/site-packages/scipy/linalg/_matfuncs_inv_ssq.py", line 684, in _fractional_matrix_power
s = svdvals(A)
File "/home/varadarajan/softwares/my-jupyter-env/lib/python3.10/site-packages/scipy/linalg/_decomp_svd.py", line 225, in svdvals
a = _asarray_validated(a, check_finite=check_finite)
File "/home/varadarajan/softwares/my-jupyter-env/lib/python3.10/site-packages/scipy/_lib/_util.py", line 306, in _asarray_validated
a = toarray(a)
File "/home/varadarajan/softwares/my-jupyter-env/lib/python3.10/site-packages/numpy/lib/function_base.py", line 630, in asarray_chkfinite
raise ValueError(
ValueError: array must not contain infs or NaNs

Scf

from mmd.integrals.fock import formPT

No such file exist

Wrong orbital energy of He with cc-pVTZ

I have run a calculation of the orbital energies of He by running:

from mmd.molecule import * 
from mmd.scf import * 
from mmd.postscf import * 
import numpy as np

# read in geometry
geometry = './geoms/he.dat'

# init molecule and build integrals
mol = Molecule(filename=geometry,basis='cc-pvtz')
mol.build()

# do the SCF
scf = SCF(mol)
scf.RHF()

A = np.dot(np.transpose(mol.C), np.dot(mol.F, mol.C))
for i in range(len(A)):
    print(A[i,i])

This gives:
(-0.91763022293+0j)
(0.590153664669+0j)
(1.50050938826+0j)
(1.50050938826+0j)
(1.50050938826+0j)
(3.44176528807+0j)
(6.35759720056+0j)
(6.35759720056+0j)
(6.35759720056+0j)
(6.35759720056+0j)
(6.35759720056+0j)
(7.84834677038+0j)
(7.84834677038+0j)
(7.84834677038+0j)
(8.33659446818+0j)

When the same calculation is performed with the Dalton program the following energies are found:

Symmetry 1A . Hartree-Fock orbital energies:
-0.91762549 0.63664267 1.50051327 1.50051327 1.50051327
4.60856244 6.35760376 6.35760376 6.35760376 6.35760376
6.35760376 7.84835826 7.84835826 7.84835826

It can be seen that the second orbital energy is quite different (0.590153664669 and 0.63664267). I have only been able to observe differences larger than the SCF convergence when D orbitals are included in the basisset. The orbital energies might be correct when only S and P orbitals are included, maybe. I have no idea right now what the problem is.
Here is data.py with cc-pVTZ for He included (rename it from .pdf to .py, github cannot attach .py in issues).
data.pdf

Speed up by detecting which two-e integrals are zeros?

Hi,

Thanks for the great repo! I'm learning a lot from this.
It seems like the two-electron integrals produces a lot of zeros, even if the central positions are close to each other (or even exactly the same).
Is it possible to speed it up by detecting which two-e integrals are zeros?

Angular Momentum

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

Real Time Magnus

Hi Josh,
Sorry me again! In your /examples/real-time.py you have

      plt.plot(rt.time,np.asarray(rt.shape)*rt.field,label='Applied field')

however, at some point you've removed the

      import numpy as np

I have some questions, for example the pulse is set to None so the shape is zero always. This means that the dipole is never added to the field and the energy does not change which you can see by plotting rt.Energy. I'm wondering if the, rich in detail, dipole curves are an artifact of the calculation? The dipole is ~1e-15 and given that the matrix exponent is a Pade approximation...An estimate of the local error can be calculated as the norn of the difference between third and fourth order Magnus Numerical integrators based on the Magnus expansion for nonlinear dynamical systems, Hajiketabi & Casas . I notice that if you apply a pulse (your Gaussian one) then you get a smoothly oscillating dipole where both Magnus2 and Magnus4 coincide, I'm still trying to puzzle out why? I think this code is the only one on github using the Magnus expansion in this way and so is important and worth expanding (time permitting). Best Peter.

Test suite organization

Unit tests need to be labeled and indexed for ease of reference. Currently the names are meaningless.

TODO: Describe the tests in the file McMurchie-Davidson/tests/index for reference, maybe add to README or Wiki(?)

Combinations

Hi Josh,
just a small typo...

    if comb(nEle,nOrb) > 5000:
        print("Number determinants: ",comb(nEle,nOrb))
        sys.exit("FCI too expensive. Quitting.")

This is the wrong way round (eg for water 'taking 10 things 14 at a time'), it will always return 0 so 'if' is always False. You've got it the right way round a couple of lines down! Hope all well.

SCF-Direct

Hi Josh,
As promised I've raised this as an issue. The problem is with the direct=True option. With H2O 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 H2 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 H2 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

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo D3

    Bring data to life with SVG, Canvas and HTML. 📊📈🎉

Recommend Topics

  • javascript

    JavaScript (JS) is a lightweight interpreted programming language with first-class functions.

  • web

    Some thing interesting about web. New door for the world.

  • server

    A server is a program made to process requests and deliver data to clients.

  • Machine learning

    Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google ❤️ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.