import Integrals
from Integrals import erivector,overlapmatrix,nuclearmatrix,kineticmatrix,normalization
import Tools
from Tools import euclidean_norm,printmatrix,eri_index
import algopy
from algopy import UTPM, zeros,transpose
import numpy as np
import Eigen
from Eigen import eigensolver
'''
This module contains functions to obtain energy
'''
def mo_naturalorbital(D):
eigsys = eigensolver(D)
return [eigsys[0]]
def nuclearrepulsion(atoms,charges,natoms):
tmp = 0.0
for i, xyza in enumerate(atoms):
for j in range(i+1,natoms):
xyzb = atoms[j]
tmp = tmp + charges[i]*charges[j]/euclidean_norm(np.subtract(xyzb,xyza))
return tmp
def fockmatrix(Hcore,Eri,D,nbasis,alpha,dtype):
F = algopy.zeros((nbasis,nbasis),dtype=dtype)
for i in range(nbasis):
for j in range(i,nbasis):
tmp = Hcore[i,j]
for k in range(nbasis):
for l in range(nbasis):
tmp = D[k,l]*(2.0*Eri[eri_index(i,j,k,l,nbasis)]-Eri[eri_index(i,k,j,l,nbasis)])+ tmp
F[i,j] = F[j,i] = tmp
return F
def cannonicalputication(D,S):
D2 = np.dot(np.dot(D,S),D)
D3 = np.dot(np.dot(D2,S),D)
c = np.trace(np.subtract(D2,D3))/np.trace(np.subtract(D,D2)) ## Eq 17 Manolopoulos
if c <= 0.5:
tmp = np.subtract(np.add(np.multiply((1.0-2.0*c),D),np.multiply((1.0+c),D2)),D3)/(1.0-c)
else:
tmp = np.subtract(np.multiply((1.0+c),D2),D3)/(c)
return np.divide(tmp,np.trace(np.dot(D,S))) ## I am not dividing by one here!!
def newdensity(F,Sinv,nbasis,ne):
F_offdiag = []
for i in range(nbasis):
F_offdiag.append(0.0)
for i in range(nbasis):
for j in range(nbasis):
if i != j:
F_offdiag[i] = np.absolute(F[i,j]) + F_offdiag[i]
listFmin = []
listFmax = []
for i in range(nbasis):
listFmin.append(F[i][i] - F_offdiag[i])
listFmax.append(F[i][i] + F_offdiag[i])
Fmax = max(listFmax)
Fmin = min(listFmin)
mubar = np.trace(F)/nbasis
lamb = min([(nbasis-ne)/(mubar-Fmin),ne/(Fmax-mubar)])
return np.add(np.multiply(lamb/nbasis,np.subtract(mubar*Sinv,np.dot(np.dot(Sinv,F),Sinv))),np.multiply(np.float64(1.0*ne/nbasis),Sinv))
[docs]def rhfenergy(alpha_old,coef2,xyz,l,charges,xyz_atom,natoms,nbasis,contr_list,ne,max_scf,max_d,log,eigen,printguess,readguess,name,write,dtype):
'''
This function returns the rhf function
Parameters:
alpha_old : array
Gaussian exponents
coef2 : array
Contraction coeffients
xyz : array 3N
Gaussian centers
l : array 3N
Angular momentum each entry is a vector
eg. s orbital (0,0,0) or pz (1,0,0)
charges : array
Atom charges
nbasis : int
Number of basis
contr_list: list of integers
Specify the number of orbitals in each atom
ne : int
Number of electrons
max_scf : int
maximum number of scf cycles
log : bool
The exponents are given in log
printguess: str or None
File to print coeff matrix initial guess
readguess : str or None
File that contains coeff matrix initial guess
name : str
Output file name
write : bool
True if printing
dtype : type of output
This is the directive to know if algopy will be used or not
np.float64(1.0) if it is a single point calculation
otherwise, it specify the size of the UTMP, autodifferentiation
Returns:
energy : float
RHF energy
'''
tool_D = 1e-8
tool = 1e-8
if log:
alpha = algopy.exp(alpha_old)
else:
alpha = alpha_old
if type(xyz_atom) != np.ndarray: ## Cover the case of diff xyz atom
coef = normalization(alpha,coef2,l,contr_list,dtype= np.float64(1.0))
V = nuclearmatrix(alpha,coef,xyz,l,nbasis,charges,xyz_atom,natoms,contr_list,dtype=dtype)
S = overlapmatrix(alpha,coef,xyz,l,nbasis,contr_list,dtype=np.float64(1.0))
T = kineticmatrix(alpha,coef,xyz,l,nbasis,contr_list,dtype=np.float64(1.0))
Eri = erivector(alpha,coef,xyz,l,nbasis,contr_list,dtype=np.float(1.0))
else:
coef = normalization(alpha,coef2,l,contr_list,dtype=dtype)
S = overlapmatrix(alpha,coef,xyz,l,nbasis,contr_list,dtype=dtype)
V = nuclearmatrix(alpha,coef,xyz,l,nbasis,charges,xyz_atom,natoms,contr_list,dtype=dtype)
T = kineticmatrix(alpha,coef,xyz,l,nbasis,contr_list,dtype=dtype)
Eri = erivector(alpha,coef,xyz,l,nbasis,contr_list,dtype=dtype)
Hcore = T + V
if eigen:
eigsys = eigensolver(S)
SqrtLambda = algopy.diag(1./algopy.sqrt(eigsys[0]))
L = eigsys[1]
LT = algopy.transpose(L)
SqrtS = algopy.dot(algopy.dot(L, SqrtLambda), LT)
SqrtST = algopy.transpose(SqrtS)
else:
Sinv = np.linalg.inv(S)
if readguess != None:
C = np.load(readguess)
D = np.zeros((nbasis,nbasis))
for i in range(nbasis):
for j in range(nbasis):
tmp = 0.0
for k in range(ne):
tmp = tmp + C[i,k]*C[j,k]
D[i,j] = tmp
F = fockmatrix(Hcore,Eri,D,nbasis,alpha,dtype)
else:
F = Hcore
OldE = 1e8
status = False
E_step = []
for scf_iter in range(max_scf):
if eigen:
Fprime = algopy.dot(algopy.dot(SqrtST,F), SqrtS)
eigsysFockOp = eigensolver(Fprime)
Cprime = eigsysFockOp[1]
C = algopy.dot(SqrtS, Cprime)
Fprime = algopy.dot(algopy.dot(SqrtST,F), SqrtS)
eigsysFockOp = eigensolver(Fprime)
Cprime = eigsysFockOp[1]
C = algopy.dot(SqrtS, Cprime)
D = algopy.zeros((nbasis,nbasis),dtype=dtype)
for i in range(nbasis):
for j in range(nbasis):
tmp = 0.0
for k in range(ne):
tmp = tmp + C[i,k]*C[j,k]
D[i,j] = tmp
else:
D = newdensity(F,Sinv,nbasis,ne)
for i in range(max_d):
D = cannonicalputication(D,S)
err = np.linalg.norm(D - np.dot(np.dot(D,S),D))
if err < tool_D:
break
F = fockmatrix(Hcore,Eri,D,nbasis,alpha,dtype)
E_elec = algopy.sum(np.multiply(D,Hcore+F))
E_step.append(E_elec)
E_nuc = nuclearrepulsion(xyz_atom,charges,natoms)
if np.absolute(E_elec - OldE) < tool:
status = True
break
OldE = E_elec
E_nuc = nuclearrepulsion(xyz_atom,charges,natoms)
if printguess !=None:
np.save(printguess, C)
def update_system():
mol.energy = E_elec + E_nuc
mol.erepulsion = Eri
mol.hcore = Hcore
mol.mo_coeff = C
return
def write_molden():
import Data
from Data import select_atom
## Details of calculation
tape.write('[Energy] \n')
tape.write('E_elec: '+str(E_elec)+'\n')
tape.write('E_nuc: '+str(E_nuc)+'\n')
tape.write('E_tot: '+str(E_nuc+E_elec)+'\n')
tape.write('SCF Details\n')
line = 'Eigen: '
if eigen:
tape.write(line+'True')
else:
tape.write(line+'False')
tape.write('\n')
for i, step in enumerate(E_step):
line = 'Step: '+str(i)+' '+str(step)
tape.write(line+'\n')
### C Matrix
tape.write('[CM] \n')
tape.write('C AO times MO\n')
printmatrix(C,tape)
### D Matrix
tape.write('[DM] \n')
tape.write('D \n')
printmatrix(D,tape)
### D Matrix
tape.write('[NMO] \n')
tape.write('NMO \n')
printmatrix(mo_naturalorbital(D),tape)
### MO energies
tape.write('[MOE] \n')
tape.write('MOE \n')
for i,energ in enumerate(eigsysFockOp[0]):
tape.write(str(i)+' '+str(energ)+'\n')
### MO energies
tape.write('[INPUT] \n')
line = 'mol = ['
for i,coord in enumerate(xyz_atom):
line += '('+str(charges[i])+','
line +='('+ str(coord[0])+','+str(coord[1])+','+str(coord[2]) + ')),\n'
tape.write(line)
cont = 0
line = 'basis = ['
for i,ci in enumerate(contr_list):
line +='['
line += '('+str(l[i][0])+','+str(l[i][1])+','+str(l[i][2])+'),'
for ii in range(ci):
line += str(alpha[cont])+','+str(coef[i])
line +=',('+ str(xyz[i,0])+','+ str(xyz[i,1])+','+str(xyz[i,2]) + ')],\n'
cont +=1
line +=']\n'
tape.write(line)
### Atom coordinates
tape.write('[Atoms]\n')
for i,coord in enumerate(xyz_atom):
line = select_atom.get(charges[i])
line +=' '+str(i+1)+' '+str(charges[i])
line +=' '+ str(coord[0])+' '+str(coord[1])+' '+str(coord[2]) + '\n'
tape.write(line)
### Basis coordinates
for i,coord in enumerate(xyz):
line = 'XX'
line +=' '+str(i+natoms+1)+' '+str(0)
line +=' '+ str(coord[0])+' '+ str(coord[1])+' '+str(coord[2]) + '\n'
tape.write(line)
### Basis set
cont = 0
tape.write('[GTO]\n')
for i,ci in enumerate(contr_list):
tape.write(' '+str(i+1+natoms)+' 0\n')
if np.sum(l[i]) == 0 :
tape.write(' s '+str(ci)+' 1.0 '+ str(l[i][0])+' '+str(l[i][1])+' '+str(l[i][2])+'\n')
else:
tape.write(' p '+str(ci)+' 1.0 '+ str(l[i][0])+' '+str(l[i][1])+' '+str(l[i][2])+'\n')
#tape.write(' p '+str(1)+' 1.0 '+ str(l[i])+'\n')
for ii in range(ci):
line =' '+str(alpha[cont])+' '+str(coef[cont])+'\n'
tape.write(line)
cont +=1
line =' \n'
tape.write(line)
### MOs
tape.write('[MO]\n')
for j in range(nbasis):
tape.write(' Sym= None\n')
tape.write(' Ene= '+str(eigsysFockOp[0][j])+'\n')
tape.write(' Spin= Alpha\n')
if j > ne:
tape.write(' Occup= 0.0\n')
else:
tape.write(' Occup= 2.0\n')
for i in range(nbasis):
tape.write(str(i+1+natoms)+' '+str(C[i,j])+'\n')
if status:
if write:
tape = open(name+'.molden',"w")
write_molden()
tape.close()
return E_elec+E_nuc
else:
print('E_elec: '+str(E_elec)+'\n')
print('E_nuc: '+str(E_nuc)+'\n')
print('E_tot: '+str(E_nuc+E_elec)+'\n')
print('SCF DID NOT CONVERGED')
return 99999
return E_elec+E_nuc
def penalty_inverse(alpha,coef3,x,y,z,l,charges,x_atom,y_atom,z_atom,natoms,nbasis,ne,max_scf,max_d,lbda,eigen,name,write):
energy = rhfenergy(alpha,coef3,x,y,z,l,charges,x_atom,y_atom,z_atom,natoms,nbasis,ne,max_scf,max_d,eigen,name,write)
## Rebuild xyz for basis
x = np.reshape(x,(1,nbasis))
y = np.reshape(y,(1,nbasis))
z = np.reshape(z,(1,nbasis))
xyz = np.concatenate((np.concatenate((x,y),axis=0),z),axis=0).T
coef = np.ones(nbasis,dtype=dtype)
coef = normalization(alpha,coef,l,nbasis)
S = overlapmatrix2(alpha,coef,l,nbasis)
penalty = lbda*1.0e-3/np.sqrt(np.linalg.det(S))
return energy + penalty