Hacked together in about an hour after I discovered matrix.rb in ruby's stdlib

hf.rb

#!/bin/env ruby
# ------------------------------------------------------------------------
# Copyright (C) 2018 Sebastian Ehlert
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
# ------------------------------------------------------------------------

require 'optparse'
require 'matrix'

options = Hash.new
usage = "Usage: ruby main.rb <input>"

OptionParser.new do |opt|
   opt.banner = usage
   opt.on_tail "-h","--help","Show this message" do
      puts opts
      exit
   end
end.parse!

# ------------------------------------------------------------------------
lines = ARGF.readlines
exit if lines.empty?
nat,nel,nbf = lines.shift.split.map { |x| x.to_i }
xyz  = Array.new(nat) { Hash.new }
at   = Array.new(nat) { 0 }
zeta = Array.new(nbf)
aoc  = Array.new(nat)

cbf = 0
k = 0
(0...nat).each do |i|
   x,y,z,iat,ibf = lines.shift.split.map{ |x| x.to_f }
   at[i] = iat.to_i
   xyz[i] = {:x=>x,:y=>y,:z=>z}
   (0...(ibf.to_i)).each do |j|
      zeta[k] = lines.shift.split.map{ |x| x.to_f }.shift
      k += 1
   end
   aoc[i] = [cbf,k]
   cbf += ibf.to_i
end
raise 'basis function missmatch!' unless cbf==nbf

# ------------------------------------------------------------------------
# nuclear nuclear repulsion
enuc = 0.0
(0...nat).each do |i|
   (0...i).each do |j|
      rx = xyz[j][:x] - xyz[i][:x]
      ry = xyz[j][:y] - xyz[i][:y]
      rz = xyz[j][:z] - xyz[i][:z]
      rab = rx*rx + ry*ry + rz*rz
      enuc += at[i]*at[j]/Math::sqrt(rab)
   end
end

# ------------------------------------------------------------------------
def oneint nat,xyz,chrg,r_a,r_b,alp,bet
   rx = r_b[:x] - r_a[:x]
   ry = r_b[:y] - r_a[:y]
   rz = r_b[:z] - r_a[:z]
   rab = rx*rx + ry*ry + rz*rz
   eab = alp+bet
   oab = 1.0/eab
   xab = alp*bet*oab
   est = rab*xab
   ab = Math::exp -est
   s = ab * (4.0*oab*xab)**0.75
   t = xab*(3.0-2*est)*s
   v = 0.0
   r_p = Hash.new
   r_p[:x] = (alp*r_a[:x]+bet*r_b[:x])*oab
   r_p[:y] = (alp*r_a[:y]+bet*r_b[:y])*oab
   r_p[:z] = (alp*r_a[:z]+bet*r_b[:z])*oab
   (0...nat).each do |k|
      rx = xyz[k][:x] - r_p[:x]
      ry = xyz[k][:y] - r_p[:y]
      rz = xyz[k][:z] - r_p[:z]
      rcp = rx*rx + ry*ry + rz*rz
      if (rcp < 10e-10)
         v -= chrg[k]*(1.0-eab*rcp/3.0)
      else
         rcp = Math::sqrt(rcp)
         v -= chrg[k]*Math::erf(Math::sqrt(eab)*rcp)/rcp
      end
   end
   [s,t,v]
end
# ------------------------------------------------------------------------
# initialize integrals
s = Array.new(nbf) {Array.new(nbf)}
t = Array.new(nbf) {Array.new(nbf)}
v = Array.new(nbf) {Array.new(nbf)}

(0...nat).each do |i|
   (0..i).each do |j|
      ((aoc[i][0])...(aoc[i][1])).each do |ii|
         ((aoc[j][0])...(aoc[j][1])).each do |jj|
            sdum,tdum,vdum = oneint nat,xyz,at,xyz[i],xyz[j],\
                                    zeta[ii],zeta[jj]
            s[ii][jj] = sdum; s[jj][ii] = sdum
            t[ii][jj] = tdum; t[jj][ii] = tdum
            v[ii][jj] = vdum; v[jj][ii] = vdum
         end
      end
   end
end


# ------------------------------------------------------------------------
def twoint r_a,r_b,r_c,r_d,alp,bet,gam,del,sab,scd
   eab = alp+bet
   oab = 1.0/eab
   r_p = Hash.new
   r_p[:x] = (alp*r_a[:x]+bet*r_b[:x])*oab
   r_p[:y] = (alp*r_a[:y]+bet*r_b[:y])*oab
   r_p[:z] = (alp*r_a[:z]+bet*r_b[:z])*oab
   ecd = gam+del
   ocd = 1.0/ecd
   r_q = Hash.new
   r_q[:x] = (gam*r_c[:x]+del*r_d[:x])*ocd
   r_q[:y] = (gam*r_c[:y]+del*r_d[:y])*ocd
   r_q[:z] = (gam*r_c[:z]+del*r_d[:z])*ocd
   rx = r_p[:x] - r_q[:x]
   ry = r_p[:y] - r_q[:y]
   rz = r_p[:z] - r_q[:z]
   rpq = rx*rx + ry*ry + rz*rz
   xpq = eab*ecd/(eab+ecd)
   if (rpq < 10e-10)
      eri = sab*scd*(1.0-xpq*rpq/3.0)
   else
      rcp = Math::sqrt(rpq)
      eri = sab*scd*Math::erf(Math::sqrt(xpq)*rpq)/rpq
   end
   eri
end
# ------------------------------------------------------------------------
eri = Array.new (nbf*(nbf+1)/2)*(nbf*(nbf+1)/2+1)/2

(0...nat).each do |i|
   (0..i).each do |j|
      (0..i).each do |k|
         limit = i==k ? j : k
         (0..limit).each do |l|
            ((aoc[i][0])...(aoc[i][1])).each do |ii|
               ((aoc[j][0])...(aoc[j][1])).each do |jj|
                  ij = ii > jj ? ii*(ii+1)/2 + jj : jj*(jj+1)/2 + ii
                  ((aoc[k][0])...(aoc[k][1])).each do |kk|
                     ((aoc[l][0])...(aoc[l][1])).each do |ll|
                        kl = kk > ll ? kk*(kk+1)/2 + ll : ll*(ll+1)/2 + kk
                        ijkl = ij > kl ? ij*(ij+1)/2 + kl : kl*(kl+1)/2 + ij
                        eri[ijkl] = twoint xyz[i],xyz[j],xyz[k],xyz[l],\
                                           zeta[ii],zeta[jj],zeta[kk],zeta[ll],\
                                           s[ii][jj],s[kk][ll]
                     end
                  end
               end
            end
         end
      end
   end
end

# ------------------------------------------------------------------------
# build orthonormalizer
x = Matrix[*s] ** -0.5

# ------------------------------------------------------------------------
# hcore guess
h = Matrix[*t] + Matrix[*v]
ff = x*h
cc,eigv,cci = ff.eigen
c = x*cc

# ------------------------------------------------------------------------
# initial density
p = c.minor(0...nbf,0...(nel/2))*c.minor(0...nbf,0...(nel/2)).t

# ------------------------------------------------------------------------
def escf h,f,p
   (p*(h+f)).tr
end
# ------------------------------------------------------------------------
# initial energy
eold = escf h,h,p
puts "initial electronic energy: #{eold+enuc}"

# ------------------------------------------------------------------------
def build_fock nbf,h,p,eri
   f = Array.new(nbf) {Array.new(nbf)}
   (0...nbf).each do |i|
      (0...nbf).each do |j|
         f[i][j] = h.element i,j
         ij = i>j ? i*(i+1)/2 + j : j*(j+1)/2 + i
         (0...nbf).each do |k|
            ik = i>k ? i*(i+1)/2 + k : k*(k+1)/2 + i
            (0...nbf).each do |l|
               kl = k>l ? k*(k+1)/2 + l : l*(l+1)/2 + k
               jl = j>l ? j*(j+1)/2 + l : l*(l+1)/2 + j
               ijkl = ij>kl ? ij*(ij+1)/2 + kl : kl*(kl+1)/2 + ij
               ikjl = ik>jl ? ik*(ik+1)/2 + jl : jl*(jl+1)/2 + ik
               f[i][j] += (2*eri[ijkl]-eri[ikjl]) * (p.element k,l)
            end
         end
      end
   end
   Matrix[*f]
end
# ------------------------------------------------------------------------
# SCF procedure
e = enuc
iter = 0
until (e-eold).abs < 10e-9
   eold = e
   f = build_fock nbf,h,p,eri
   ff = x*f
   cc,eigv,cci = ff.eigen # eigenvectors and eigenvalues are *not* sorted!

   cvec = cc.column_vectors
   (0...nbf).each { |i| cvec[i] = [(eigv.element i,i),(cc.column i)] }
   cvec = cvec.sort.map {|x|x.pop}

   c = x*(Matrix.columns cvec)
   p = c.minor(0...nbf,0...(nel/2))*c.minor(0...nbf,0...(nel/2)).t
   e = escf h,f,p
   iter += 1
   puts "%3i %10f %10f" % [iter,enuc+e,e-eold]
end