#!/usr/bin/env python import math import sys from numpy import * from numpy.linalg import * def read_one_electron_integrals(filename): """Read the one-electron integrals contained in the file 'filename' and return a numpy matrix""" ints_lines = open(filename,"r").readlines() k = len(ints_lines[0].split()) ints = zeros( (k,k) ) for i,line in enumerate(ints_lines): for j,value in enumerate(line.split()): ints[i][j] = float(value) return ints def read_two_electron_integrals(filename,nao): """Read the two-electron integrals contained in the file 'filename' and return a 4-dimensional numpy array that stores the integrals in chemistry notation as ints[m][n][r][s] = (mn|rs). This function needs to know the total number of atomic orbitals 'nao'.""" ints_lines = open(filename,"r").readlines() ints = zeros( (nao,nao,nao,nao) ) for line in ints_lines: line_split = line.split() p,q,r,s = [int(t) for t in line_split[0:4]] value = float(line_split[-1]) ints[p][q][r][s] = ints[q][p][r][s] = value ints[p][q][s][r] = ints[q][p][s][r] = value ints[r][s][p][q] = ints[r][s][q][p] = value ints[s][r][p][q] = ints[s][r][q][p] = value return ints # Define the number of atomic orbitals (nao) nao = 7 # Define the number of electrons N = 10 # Define the number of filled spatial orbitals Nhalf = N / 2 # Read the integrals VNN = float(open("vnn","r").readlines()[0]) S = read_one_electron_integrals("overlap") H = read_one_electron_integrals("one-electron") V = read_two_electron_integrals("two-electron",nao) # Adjust the printing of the numpy matrices set_printoptions(precision=5,linewidth=100) # Print the S and H matrices print "S matrix:\n",S print "H matrix:\n",H # Print the two-electron integrals in chemist notation for m in xrange(nao): for n in xrange(nao): for r in xrange(nao): for s in xrange(nao): print "(%d %d|%d %d) = %20.12f" % V[m][n][r][s]