Created
May 12, 2023 02:29
-
-
Save Eipgen/8af4d5f919a81c26300cb6302223c8df to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
import numpy as np | |
from scipy.special import erf | |
def get_parameters(qeqfile): | |
eV = 3.67493245e-2 | |
Angstrom = 1./0.529177249 | |
parameters = {} | |
with open(qeqfile) as f: | |
for i in range(13): | |
next(f) | |
for line in f: | |
data = line.rstrip().split() | |
element = data[0] | |
elecnegativity = float(data[1])*eV | |
hardness = float(data[2])*eV | |
sradius = float(data[3])*Angstrom | |
basis = 1.0/(sradius*sradius) | |
parameters[element] = [elecnegativity, hardness, basis] | |
return parameters | |
def read_total_charge(infile): | |
with open(infile) as f: | |
for line in f: | |
if line.startswith("TCHARGE"): | |
total_charge = float(line.rstrip().split()[1]) | |
return total_charge | |
def get_elements(file): | |
atomElements = {} | |
F=open(file).readlines()[2:] | |
for line in F: | |
data = line.rstrip().split() | |
atomElements[int(data[4])] = data[0] | |
return atomElements | |
def atom_info(file): | |
F=open(file).readlines()[2:] | |
atoms = [] | |
for line in F: | |
data = line.rstrip().split() | |
atoms.append([float(data[1]),float(data[2]),float(data[3]),int(data[4])]) | |
return atoms | |
def calculate_coulomb_intergral(a,b,R): | |
""" | |
a,b are the basis of the two atoms | |
R is the distance between the two atoms | |
""" | |
p = np.sqrt(a*b/(a+b)) | |
return erf(p*R)/R | |
def fill_J(atoms, J, BasisSet, CoulombMaxDistance): | |
Angstrom = 1./0.529177249 | |
nAtoms = len(atoms) | |
for k in range(0, nAtoms): | |
for l in range(0,nAtoms): | |
if k >l: | |
atom1 = atoms[k] | |
atom2 = atoms[l] | |
xyz_atom1 = np.array(atom1[0:3]) | |
xyz_atom2 = np.array(atom2[0:3]) | |
R = np.linalg.norm(xyz_atom1-xyz_atom2)*Angstrom | |
if R < CoulombMaxDistance: | |
coulomb = calculate_coulomb_intergral(BasisSet[k],BasisSet[l],R) | |
else: | |
coulomb = 1.0/R | |
J[k][l] = coulomb | |
J[l][k] = coulomb | |
for i in range(nAtoms+1): | |
J[nAtoms][i] = 1.0 | |
J[i][nAtoms] = 1.0 | |
J[nAtoms][nAtoms] = 0.0 | |
def compute_Qeq_charges(atoms,atomElments,total_charge,charge_past): | |
CoulombThreshold = 1e-9 | |
nAtoms = len(atoms) | |
ElectroNegativity = np.zeros(nAtoms) | |
J = np.zeros((nAtoms+1,nAtoms+1)) # hardness | |
Voltage = charge_past | |
BasisSet = np.zeros(nAtoms) | |
parameters = get_parameters("qeq_reaxff2016.txt") | |
#atomElements = get_elements("./atom.xyz") | |
for i in range(nAtoms): | |
atomType = atoms[i][3] | |
atomElement = atomElements[atomType] | |
ElectroNegativity[i] = parameters[atomElement][0] | |
J[i][i] = parameters[atomElement][1] | |
BasisSet[i] = parameters[atomElement][2] | |
#total_charge = read_total_charge("./atom.xyz") | |
print("Total charge: ", total_charge) | |
SmallestGaussianExponent = min(BasisSet) | |
CoulombMaxDistance = 2*np.sqrt(-np.log(CoulombThreshold)/SmallestGaussianExponent) | |
fill_J(atoms, J, BasisSet, CoulombMaxDistance) | |
print("Hardness") | |
for j in J: | |
j=[str(round(i,6)) for i in j] | |
print(" ".join(j)) | |
Voltage[:-1] = ElectroNegativity | |
Voltage[-1] = total_charge | |
print("Voltage") | |
for i in Voltage: | |
print(round(i,6)) | |
charges = np.linalg.solve(J,Voltage) | |
return charges | |
def Qeq_charge_equilibration(atoms,atomElements,total_charge): | |
nAtoms = len(atoms) | |
new_charges = np.zeros(nAtoms+1) | |
for i in range(1): | |
charge_past = new_charges | |
new_charges = compute_Qeq_charges(atoms,atomElements,total_charge,charge_past) | |
qeq_charges = new_charges | |
return qeq_charges | |
atoms=atom_info("./atom.xyz") | |
atomElements = get_elements("./atom.xyz") | |
total_charge = read_total_charge("./atom.xyz") | |
qeq_charges = Qeq_charge_equilibration(atoms,atomElements,total_charge) | |
print(qeq_charges) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment