MolHandling Translation - Most
In this commit most of the work about Molecules and the Quantum Mechanical System has been finished, the only part still missing if a method that populates de ASSEC, but this method requires the Classes: Dice, Player and Gaussian Signed-off-by: Vitor Hideyoshi <vitor.h.n.batista@gmail.com>
This commit is contained in:
@@ -152,75 +152,75 @@ from DPpack.SetGlobals import *
|
|||||||
|
|
||||||
|
|
||||||
|
|
||||||
def nearest_image(refmol, molecule, lx, ly, lz, criterium="com"):
|
# def nearest_image(refmol, molecule, lx, ly, lz, criterium="com"):
|
||||||
|
|
||||||
if criterium != "com" and criterium != "min":
|
# if criterium != "com" and criterium != "min":
|
||||||
sys.exit("Error in value passed to function nearest_image")
|
# sys.exit("Error in value passed to function nearest_image")
|
||||||
min_dist = 1e20
|
# min_dist = 1e20
|
||||||
for i in range(-1, 2):
|
# for i in range(-1, 2):
|
||||||
for j in range(-1, 2):
|
# for j in range(-1, 2):
|
||||||
for k in range(-1, 2):
|
# for k in range(-1, 2):
|
||||||
|
|
||||||
tr_vector = [i * lx, j * ly, k * lz]
|
# tr_vector = [i * lx, j * ly, k * lz]
|
||||||
new_molecule = translate(molecule, tr_vector)
|
# new_molecule = translate(molecule, tr_vector)
|
||||||
if criterium == "com":
|
# if criterium == "com":
|
||||||
dist = center_of_mass_distance(refmol, new_molecule)
|
# dist = center_of_mass_distance(refmol, new_molecule)
|
||||||
else:
|
# else:
|
||||||
dist = minimum_distance(refmol, new_molecule)
|
# dist = minimum_distance(refmol, new_molecule)
|
||||||
|
|
||||||
if dist < min_dist:
|
# if dist < min_dist:
|
||||||
min_dist = dist
|
# min_dist = dist
|
||||||
nearestmol = deepcopy(new_molecule)
|
# nearestmol = deepcopy(new_molecule)
|
||||||
|
|
||||||
return min_dist, nearestmol
|
# return min_dist, nearestmol
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
def calculate_step(gradient, hessian, fh):
|
# def calculate_step(gradient, hessian, fh):
|
||||||
|
|
||||||
invhessian = linalg.inv(hessian)
|
# invhessian = linalg.inv(hessian)
|
||||||
pre_step = -1 * np.matmul(invhessian, gradient.T).T
|
# pre_step = -1 * np.matmul(invhessian, gradient.T).T
|
||||||
maxstep = np.amax(np.absolute(pre_step))
|
# maxstep = np.amax(np.absolute(pre_step))
|
||||||
factor = min(1, player['maxstep']/maxstep)
|
# factor = min(1, player['maxstep']/maxstep)
|
||||||
step = factor * pre_step
|
# step = factor * pre_step
|
||||||
|
|
||||||
fh.write("\nCalculated step:\n")
|
# fh.write("\nCalculated step:\n")
|
||||||
pre_step_list = pre_step.tolist()
|
# pre_step_list = pre_step.tolist()
|
||||||
|
|
||||||
fh.write("-----------------------------------------------------------------------\n"
|
# fh.write("-----------------------------------------------------------------------\n"
|
||||||
"Center Atomic Step (Bohr)\n"
|
# "Center Atomic Step (Bohr)\n"
|
||||||
"Number Number X Y Z\n"
|
# "Number Number X Y Z\n"
|
||||||
"-----------------------------------------------------------------------\n")
|
# "-----------------------------------------------------------------------\n")
|
||||||
for i in range(len(molecules[0])):
|
# for i in range(len(molecules[0])):
|
||||||
fh.write(" {:>5d} {:>3d} {:>14.9f} {:>14.9f} {:>14.9f}\n".format(
|
# fh.write(" {:>5d} {:>3d} {:>14.9f} {:>14.9f} {:>14.9f}\n".format(
|
||||||
i + 1, molecules[0][i]['na'],
|
# i + 1, molecules[0][i]['na'],
|
||||||
pre_step_list.pop(0), pre_step_list.pop(0), pre_step_list.pop(0)))
|
# pre_step_list.pop(0), pre_step_list.pop(0), pre_step_list.pop(0)))
|
||||||
|
|
||||||
fh.write("-----------------------------------------------------------------------\n")
|
# fh.write("-----------------------------------------------------------------------\n")
|
||||||
|
|
||||||
fh.write("Maximum step is {:>11.6}\n".format(maxstep))
|
# fh.write("Maximum step is {:>11.6}\n".format(maxstep))
|
||||||
fh.write("Scaling factor = {:>6.4f}\n".format(factor))
|
# fh.write("Scaling factor = {:>6.4f}\n".format(factor))
|
||||||
fh.write("\nFinal step (Bohr):\n")
|
# fh.write("\nFinal step (Bohr):\n")
|
||||||
step_list = step.tolist()
|
# step_list = step.tolist()
|
||||||
|
|
||||||
fh.write("-----------------------------------------------------------------------\n"
|
# fh.write("-----------------------------------------------------------------------\n"
|
||||||
"Center Atomic Step (Bohr)\n"
|
# "Center Atomic Step (Bohr)\n"
|
||||||
"Number Number X Y Z\n"
|
# "Number Number X Y Z\n"
|
||||||
"-----------------------------------------------------------------------\n")
|
# "-----------------------------------------------------------------------\n")
|
||||||
for i in range(len(molecules[0])):
|
# for i in range(len(molecules[0])):
|
||||||
fh.write(" {:>5d} {:>3d} {:>14.9f} {:>14.9f} {:>14.9f}\n".format(
|
# fh.write(" {:>5d} {:>3d} {:>14.9f} {:>14.9f} {:>14.9f}\n".format(
|
||||||
i + 1, molecules[0][i]['na'],
|
# i + 1, molecules[0][i]['na'],
|
||||||
step_list.pop(0), step_list.pop(0), step_list.pop(0)))
|
# step_list.pop(0), step_list.pop(0), step_list.pop(0)))
|
||||||
|
|
||||||
fh.write("-----------------------------------------------------------------------\n")
|
# fh.write("-----------------------------------------------------------------------\n")
|
||||||
|
|
||||||
step_max = np.amax(np.absolute(step))
|
# step_max = np.amax(np.absolute(step))
|
||||||
step_rms = np.sqrt(np.mean(np.square(step)))
|
# step_rms = np.sqrt(np.mean(np.square(step)))
|
||||||
|
|
||||||
fh.write(" Max Step = {:>14.9f} RMS Step = {:>14.9f}\n\n".format(
|
# fh.write(" Max Step = {:>14.9f} RMS Step = {:>14.9f}\n\n".format(
|
||||||
step_max, step_rms))
|
# step_max, step_rms))
|
||||||
|
|
||||||
return step
|
# return step
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
@@ -236,21 +236,21 @@ def calculate_step(gradient, hessian, fh):
|
|||||||
|
|
||||||
|
|
||||||
|
|
||||||
def update_molecule(position, fh):
|
# def update_molecule(position, fh):
|
||||||
|
|
||||||
position_in_ang = (position * bohr2ang).tolist()
|
# position_in_ang = (position * bohr2ang).tolist()
|
||||||
new_molecule = deepcopy(molecules[0])
|
# new_molecule = deepcopy(molecules[0])
|
||||||
for atom in new_molecule:
|
# for atom in new_molecule:
|
||||||
atom['rx'] = position_in_ang.pop(0)
|
# atom['rx'] = position_in_ang.pop(0)
|
||||||
atom['ry'] = position_in_ang.pop(0)
|
# atom['ry'] = position_in_ang.pop(0)
|
||||||
atom['rz'] = position_in_ang.pop(0)
|
# atom['rz'] = position_in_ang.pop(0)
|
||||||
|
|
||||||
rmsd, molecules[0] = rmsd_fit(new_molecule, molecules[0])
|
# rmsd, molecules[0] = rmsd_fit(new_molecule, molecules[0])
|
||||||
|
|
||||||
fh.write("\nProjected new conformation of reference molecule with RMSD fit\n")
|
# fh.write("\nProjected new conformation of reference molecule with RMSD fit\n")
|
||||||
fh.write("RMSD = {:>8.5f} Angstrom\n".format(rmsd))
|
# fh.write("RMSD = {:>8.5f} Angstrom\n".format(rmsd))
|
||||||
|
|
||||||
return
|
# return
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
@@ -266,7 +266,7 @@ def update_molecule(position, fh):
|
|||||||
|
|
||||||
# return hessian
|
# return hessian
|
||||||
|
|
||||||
|
### Antes de aplicar devemos aplicar a aplicação de Classes para Dice, Player e etc
|
||||||
|
|
||||||
def populate_asec_vdw(cycle, fh):
|
def populate_asec_vdw(cycle, fh):
|
||||||
|
|
||||||
@@ -417,49 +417,49 @@ def populate_asec_vdw(cycle, fh):
|
|||||||
|
|
||||||
|
|
||||||
|
|
||||||
def print_geom(cycle, fh):
|
# def print_geom(cycle, fh):
|
||||||
|
|
||||||
fh.write("{}\n".format(len(molecules[0])))
|
# fh.write("{}\n".format(len(molecules[0])))
|
||||||
fh.write("Cycle # {}\n".format(cycle))
|
# fh.write("Cycle # {}\n".format(cycle))
|
||||||
for atom in molecules[0]:
|
# for atom in molecules[0]:
|
||||||
symbol = atomsymb[atom['na']]
|
# symbol = atomsymb[atom['na']]
|
||||||
fh.write("{:<2s} {:>10.6f} {:>10.6f} {:>10.6f}\n".format(symbol,
|
# fh.write("{:<2s} {:>10.6f} {:>10.6f} {:>10.6f}\n".format(symbol,
|
||||||
atom['rx'], atom['ry'], atom['rz']))
|
# atom['rx'], atom['ry'], atom['rz']))
|
||||||
|
|
||||||
return
|
# return
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
def print_mol_info(molecule, fh):
|
# def print_mol_info(molecule, fh):
|
||||||
|
|
||||||
com = center_of_mass(molecule)
|
# com = center_of_mass(molecule)
|
||||||
fh.write(" Center of mass = ( {:>10.4f} , {:>10.4f} , {:>10.4f} )\n".format(com[0],
|
# fh.write(" Center of mass = ( {:>10.4f} , {:>10.4f} , {:>10.4f} )\n".format(com[0],
|
||||||
com[1], com[2]))
|
# com[1], com[2]))
|
||||||
inertia = inertia_tensor(molecule)
|
# inertia = inertia_tensor(molecule)
|
||||||
evals, evecs = principal_axes(inertia)
|
# evals, evecs = principal_axes(inertia)
|
||||||
|
|
||||||
fh.write(" Moments of inertia = {:>9E} {:>9E} {:>9E}\n".format(evals[0],
|
# fh.write(" Moments of inertia = {:>9E} {:>9E} {:>9E}\n".format(evals[0],
|
||||||
evals[1], evals[2]))
|
# evals[1], evals[2]))
|
||||||
|
|
||||||
fh.write(" Major principal axis = ( {:>10.6f} , {:>10.6f} , {:>10.6f} )\n".format(
|
# fh.write(" Major principal axis = ( {:>10.6f} , {:>10.6f} , {:>10.6f} )\n".format(
|
||||||
evecs[0,0], evecs[1,0], evecs[2,0]))
|
# evecs[0,0], evecs[1,0], evecs[2,0]))
|
||||||
fh.write(" Inter principal axis = ( {:>10.6f} , {:>10.6f} , {:>10.6f} )\n".format(
|
# fh.write(" Inter principal axis = ( {:>10.6f} , {:>10.6f} , {:>10.6f} )\n".format(
|
||||||
evecs[0,1], evecs[1,1], evecs[2,1]))
|
# evecs[0,1], evecs[1,1], evecs[2,1]))
|
||||||
fh.write(" Minor principal axis = ( {:>10.6f} , {:>10.6f} , {:>10.6f} )\n".format(
|
# fh.write(" Minor principal axis = ( {:>10.6f} , {:>10.6f} , {:>10.6f} )\n".format(
|
||||||
evecs[0,2], evecs[1,2], evecs[2,2]))
|
# evecs[0,2], evecs[1,2], evecs[2,2]))
|
||||||
|
|
||||||
sizes = sizes_of_molecule(molecule)
|
# sizes = sizes_of_molecule(molecule)
|
||||||
fh.write(" Characteristic lengths = ( {:>6.2f} , {:>6.2f} , {:>6.2f} )\n".format(
|
# fh.write(" Characteristic lengths = ( {:>6.2f} , {:>6.2f} , {:>6.2f} )\n".format(
|
||||||
sizes[0], sizes[1], sizes[2]))
|
# sizes[0], sizes[1], sizes[2]))
|
||||||
mol_mass = total_mass(molecule)
|
# mol_mass = total_mass(molecule)
|
||||||
fh.write(" Total mass = {:>8.2f} au\n".format(mol_mass))
|
# fh.write(" Total mass = {:>8.2f} au\n".format(mol_mass))
|
||||||
|
|
||||||
chg_dip = charges_and_dipole(molecule)
|
# chg_dip = charges_and_dipole(molecule)
|
||||||
fh.write(" Total charge = {:>8.4f} e\n".format(chg_dip[0]))
|
# fh.write(" Total charge = {:>8.4f} e\n".format(chg_dip[0]))
|
||||||
fh.write(" Dipole moment = ( {:>9.4f} , {:>9.4f} , {:>9.4f} ) Total = {:>9.4f} Debye\n\n".format(
|
# fh.write(" Dipole moment = ( {:>9.4f} , {:>9.4f} , {:>9.4f} ) Total = {:>9.4f} Debye\n\n".format(
|
||||||
chg_dip[1], chg_dip[2], chg_dip[3], chg_dip[4]))
|
# chg_dip[1], chg_dip[2], chg_dip[3], chg_dip[4]))
|
||||||
|
|
||||||
return
|
# return
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
@@ -533,71 +533,68 @@ def print_mol_info(molecule, fh):
|
|||||||
|
|
||||||
|
|
||||||
|
|
||||||
def rmsd_fit(projecting_mol, reference_mol):
|
# def rmsd_fit(projecting_mol, reference_mol):
|
||||||
|
|
||||||
if len(projecting_mol) != len(reference_mol):
|
# if len(projecting_mol) != len(reference_mol):
|
||||||
sys.exit("Error in RMSD fit procedure: molecules have different number of atoms")
|
# sys.exit("Error in RMSD fit procedure: molecules have different number of atoms")
|
||||||
dim = len(projecting_mol)
|
# dim = len(projecting_mol)
|
||||||
|
|
||||||
new_projecting_mol = deepcopy(projecting_mol)
|
# new_projecting_mol = deepcopy(projecting_mol)
|
||||||
new_reference_mol = deepcopy(reference_mol)
|
# new_reference_mol = deepcopy(reference_mol)
|
||||||
|
|
||||||
center_of_mass_to_origin(new_projecting_mol)
|
# center_of_mass_to_origin(new_projecting_mol)
|
||||||
center_of_mass_to_origin(new_reference_mol)
|
# center_of_mass_to_origin(new_reference_mol)
|
||||||
|
|
||||||
x = []
|
# x = []
|
||||||
y = []
|
# y = []
|
||||||
|
|
||||||
for atom in new_projecting_mol:
|
# for atom in new_projecting_mol:
|
||||||
x.extend([ atom['rx'], atom['ry'], atom['rz'] ])
|
# x.extend([ atom['rx'], atom['ry'], atom['rz'] ])
|
||||||
|
|
||||||
for atom in new_reference_mol:
|
# for atom in new_reference_mol:
|
||||||
y.extend([ atom['rx'], atom['ry'], atom['rz'] ])
|
# y.extend([ atom['rx'], atom['ry'], atom['rz'] ])
|
||||||
|
|
||||||
x = np.array(x).reshape(dim, 3)
|
# x = np.array(x).reshape(dim, 3)
|
||||||
y = np.array(y).reshape(dim, 3)
|
# y = np.array(y).reshape(dim, 3)
|
||||||
|
|
||||||
r = np.matmul(y.T, x)
|
# r = np.matmul(y.T, x)
|
||||||
rr = np.matmul(r.T, r)
|
# rr = np.matmul(r.T, r)
|
||||||
|
|
||||||
try:
|
# try:
|
||||||
evals, evecs = linalg.eigh(rr)
|
# evals, evecs = linalg.eigh(rr)
|
||||||
except:
|
# except:
|
||||||
sys.exit("Error: diagonalization of RR matrix did not converge")
|
# sys.exit("Error: diagonalization of RR matrix did not converge")
|
||||||
|
|
||||||
a1 = evecs[:,2].T
|
# a1 = evecs[:,2].T
|
||||||
a2 = evecs[:,1].T
|
# a2 = evecs[:,1].T
|
||||||
a3 = np.cross(a1, a2)
|
# a3 = np.cross(a1, a2)
|
||||||
|
|
||||||
A = np.array([ a1[0], a1[1], a1[2], a2[0], a2[1], a2[2], a3[0], a3[1], a3[2] ])
|
# A = np.array([ a1[0], a1[1], a1[2], a2[0], a2[1], a2[2], a3[0], a3[1], a3[2] ])
|
||||||
A = A.reshape(3,3)
|
# A = A.reshape(3,3)
|
||||||
|
|
||||||
b1 = np.matmul(r, a1.T).T # or np.dot(r, a1)
|
# b1 = np.matmul(r, a1.T).T # or np.dot(r, a1)
|
||||||
b1 /= linalg.norm(b1)
|
# b1 /= linalg.norm(b1)
|
||||||
b2 = np.matmul(r, a2.T).T # or np.dot(r, a2)
|
# b2 = np.matmul(r, a2.T).T # or np.dot(r, a2)
|
||||||
b2 /= linalg.norm(b2)
|
# b2 /= linalg.norm(b2)
|
||||||
b3 = np.cross(b1, b2)
|
# b3 = np.cross(b1, b2)
|
||||||
|
|
||||||
B = np.array([ b1[0], b1[1], b1[2], b2[0], b2[1], b2[2], b3[0], b3[1], b3[2] ])
|
# B = np.array([ b1[0], b1[1], b1[2], b2[0], b2[1], b2[2], b3[0], b3[1], b3[2] ])
|
||||||
B = B.reshape(3,3).T
|
# B = B.reshape(3,3).T
|
||||||
|
|
||||||
rot_matrix = np.matmul(B, A)
|
# rot_matrix = np.matmul(B, A)
|
||||||
x = np.matmul(rot_matrix, x.T).T
|
# x = np.matmul(rot_matrix, x.T).T
|
||||||
|
|
||||||
rmsd = 0
|
# rmsd = 0
|
||||||
for i in range(dim):
|
# for i in range(dim):
|
||||||
rmsd += (x[i,0] - y[i,0])**2 + (x[i,1] - y[i,1])**2 + (x[i,2] - y[i,2])**2
|
# rmsd += (x[i,0] - y[i,0])**2 + (x[i,1] - y[i,1])**2 + (x[i,2] - y[i,2])**2
|
||||||
rmsd = math.sqrt(rmsd/dim)
|
# rmsd = math.sqrt(rmsd/dim)
|
||||||
|
|
||||||
for i in range(dim):
|
|
||||||
new_projecting_mol[i]['rx'] = x[i,0]
|
|
||||||
new_projecting_mol[i]['ry'] = x[i,1]
|
|
||||||
new_projecting_mol[i]['rz'] = x[i,2]
|
|
||||||
|
|
||||||
tr_vector = center_of_mass(reference_mol)
|
|
||||||
projected_mol = translate(new_projecting_mol, tr_vector)
|
|
||||||
|
|
||||||
return rmsd, projected_mol
|
|
||||||
|
|
||||||
|
# for i in range(dim):
|
||||||
|
# new_projecting_mol[i]['rx'] = x[i,0]
|
||||||
|
# new_projecting_mol[i]['ry'] = x[i,1]
|
||||||
|
# new_projecting_mol[i]['rz'] = x[i,2]
|
||||||
|
|
||||||
|
# tr_vector = center_of_mass(reference_mol)
|
||||||
|
# projected_mol = translate(new_projecting_mol, tr_vector)
|
||||||
|
|
||||||
|
# return rmsd, projected_mol
|
||||||
@@ -53,6 +53,133 @@ class System:
|
|||||||
|
|
||||||
return min(distances)
|
return min(distances)
|
||||||
|
|
||||||
|
def rmsd_fit(self, index_p, index_r):
|
||||||
|
|
||||||
|
projecting_mol = self.molecule[index_p]
|
||||||
|
reference_mol = self.molecule[index_r]
|
||||||
|
|
||||||
|
if len(projecting_mol.atom) != len(reference_mol.atom):
|
||||||
|
sys.exit("Error in RMSD fit procedure: molecules have different number of atoms")
|
||||||
|
dim = len(projecting_mol.atom)
|
||||||
|
|
||||||
|
new_projecting_mol = deepcopy(projecting_mol)
|
||||||
|
new_reference_mol = deepcopy(reference_mol)
|
||||||
|
|
||||||
|
new_projecting_mol.center_of_mass_to_origin()
|
||||||
|
new_reference_mol.center_of_mass_to_origin()
|
||||||
|
|
||||||
|
x = []
|
||||||
|
y = []
|
||||||
|
|
||||||
|
for atom in new_projecting_mol:
|
||||||
|
x.extend([ atom.rx, atom.ry, atom.rz ])
|
||||||
|
|
||||||
|
for atom in new_reference_mol:
|
||||||
|
y.extend([ atom.rx, atom.ry, atom.rz ])
|
||||||
|
|
||||||
|
x = np.array(x).reshape(dim, 3)
|
||||||
|
y = np.array(y).reshape(dim, 3)
|
||||||
|
|
||||||
|
r = np.matmul(y.T, x)
|
||||||
|
rr = np.matmul(r.T, r)
|
||||||
|
|
||||||
|
try:
|
||||||
|
evals, evecs = linalg.eigh(rr)
|
||||||
|
except:
|
||||||
|
sys.exit("Error: diagonalization of RR matrix did not converge")
|
||||||
|
|
||||||
|
a1 = evecs[:,2].T
|
||||||
|
a2 = evecs[:,1].T
|
||||||
|
a3 = np.cross(a1, a2)
|
||||||
|
|
||||||
|
A = np.array([ a1[0], a1[1], a1[2], a2[0], a2[1], a2[2], a3[0], a3[1], a3[2] ])
|
||||||
|
A = A.reshape(3,3)
|
||||||
|
|
||||||
|
b1 = np.matmul(r, a1.T).T # or np.dot(r, a1)
|
||||||
|
b1 /= linalg.norm(b1)
|
||||||
|
b2 = np.matmul(r, a2.T).T # or np.dot(r, a2)
|
||||||
|
b2 /= linalg.norm(b2)
|
||||||
|
b3 = np.cross(b1, b2)
|
||||||
|
|
||||||
|
B = np.array([ b1[0], b1[1], b1[2], b2[0], b2[1], b2[2], b3[0], b3[1], b3[2] ])
|
||||||
|
B = B.reshape(3,3).T
|
||||||
|
|
||||||
|
rot_matrix = np.matmul(B, A)
|
||||||
|
x = np.matmul(rot_matrix, x.T).T
|
||||||
|
|
||||||
|
rmsd = 0
|
||||||
|
for i in range(dim):
|
||||||
|
rmsd += (x[i,0] - y[i,0])**2 + (x[i,1] - y[i,1])**2 + (x[i,2] - y[i,2])**2
|
||||||
|
rmsd = math.sqrt(rmsd/dim)
|
||||||
|
|
||||||
|
for i in range(dim):
|
||||||
|
new_projecting_mol.atom[i].rx = x[i,0]
|
||||||
|
new_projecting_mol.atom[i].ry = x[i,1]
|
||||||
|
new_projecting_mol.atom[i].rz = x[i,2]
|
||||||
|
|
||||||
|
tr_vector = reference_mol.center_of_mass()
|
||||||
|
projected_mol = new_projecting_mol.translate(tr_vector)
|
||||||
|
|
||||||
|
return rmsd, projected_mol
|
||||||
|
|
||||||
|
def update_molecule(self, position, fh):
|
||||||
|
|
||||||
|
position_in_ang = (position * bohr2ang).tolist()
|
||||||
|
self.add_molecule(deepcopy(self.molecule[0]))
|
||||||
|
|
||||||
|
for atom in self.molecule[-1].atom:
|
||||||
|
|
||||||
|
atom.rx = position_in_ang.pop(0)
|
||||||
|
atom.ry = position_in_ang.pop(0)
|
||||||
|
atom.rz = position_in_ang.pop(0)
|
||||||
|
|
||||||
|
rmsd, self.molecule[0] = self.rmsd_fit(-1, 0)
|
||||||
|
self.molecule.pop(-1)
|
||||||
|
|
||||||
|
fh.write("\nProjected new conformation of reference molecule with RMSD fit\n")
|
||||||
|
fh.write("RMSD = {:>8.5f} Angstrom\n".format(rmsd))
|
||||||
|
|
||||||
|
def nearest_image(self, index_r, index_m, lx, ly, lz, criterium=None):
|
||||||
|
|
||||||
|
if criterium in None:
|
||||||
|
criterium = "com"
|
||||||
|
|
||||||
|
if criterium != "com" and criterium != "min":
|
||||||
|
sys.exit("Error in value passed to function nearest_image")
|
||||||
|
|
||||||
|
min_dist = 1e20
|
||||||
|
|
||||||
|
for i in range(-1, 2):
|
||||||
|
for j in range(-1, 2):
|
||||||
|
for k in range(-1, 2):
|
||||||
|
|
||||||
|
tr_vector = [i * lx, j * ly, k * lz]
|
||||||
|
self.add_molecule(self.molecule[index_m].translate(tr_vector))
|
||||||
|
|
||||||
|
if criterium == "com":
|
||||||
|
dist = self.center_of_mass_distance(index_r, -1)
|
||||||
|
else:
|
||||||
|
dist = self.minimum_distance(index_r, -1)
|
||||||
|
|
||||||
|
if dist < min_dist:
|
||||||
|
min_dist = dist
|
||||||
|
nearestmol = deepcopy(self.molecule[-1])
|
||||||
|
|
||||||
|
self.molecule.pop(-1)
|
||||||
|
|
||||||
|
return min_dist, nearestmol
|
||||||
|
|
||||||
|
def print_geom(self, cycle, fh):
|
||||||
|
|
||||||
|
fh.write("{}\n".format(len(self.molecule[0])))
|
||||||
|
fh.write("Cycle # {}\n".format(cycle))
|
||||||
|
for atom in self.molecule[0].atoms:
|
||||||
|
symbol = atomsymb[atom.na]
|
||||||
|
fh.write("{:<2s} {:>10.6f} {:>10.6f} {:>10.6f}\n".format(symbol,
|
||||||
|
atom.rx, atom.ry, atom.rz))
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
# Classe que conterá toda informação e funções relacionadas a uma unica molecula
|
# Classe que conterá toda informação e funções relacionadas a uma unica molecula
|
||||||
|
|
||||||
class Molecule:
|
class Molecule:
|
||||||
@@ -262,6 +389,81 @@ class Molecule:
|
|||||||
|
|
||||||
return new_molecule
|
return new_molecule
|
||||||
|
|
||||||
|
def print_mol_info(self, fh):
|
||||||
|
|
||||||
|
com = self.center_of_mass()
|
||||||
|
fh.write(" Center of mass = ( {:>10.4f} , {:>10.4f} , {:>10.4f} )\n".format(com[0],
|
||||||
|
com[1], com[2]))
|
||||||
|
inertia = self.inertia_tensor()
|
||||||
|
evals, evecs = self.principal_axes()
|
||||||
|
|
||||||
|
fh.write(" Moments of inertia = {:>9E} {:>9E} {:>9E}\n".format(evals[0],
|
||||||
|
evals[1], evals[2]))
|
||||||
|
|
||||||
|
fh.write(" Major principal axis = ( {:>10.6f} , {:>10.6f} , {:>10.6f} )\n".format(
|
||||||
|
evecs[0,0], evecs[1,0], evecs[2,0]))
|
||||||
|
fh.write(" Inter principal axis = ( {:>10.6f} , {:>10.6f} , {:>10.6f} )\n".format(
|
||||||
|
evecs[0,1], evecs[1,1], evecs[2,1]))
|
||||||
|
fh.write(" Minor principal axis = ( {:>10.6f} , {:>10.6f} , {:>10.6f} )\n".format(
|
||||||
|
evecs[0,2], evecs[1,2], evecs[2,2]))
|
||||||
|
|
||||||
|
sizes = self.sizes_of_molecule()
|
||||||
|
fh.write(" Characteristic lengths = ( {:>6.2f} , {:>6.2f} , {:>6.2f} )\n".format(
|
||||||
|
sizes[0], sizes[1], sizes[2]))
|
||||||
|
mol_mass = self.total_mass()
|
||||||
|
fh.write(" Total mass = {:>8.2f} au\n".format(mol_mass))
|
||||||
|
|
||||||
|
chg_dip = self.charges_and_dipole()
|
||||||
|
fh.write(" Total charge = {:>8.4f} e\n".format(chg_dip[0]))
|
||||||
|
fh.write(" Dipole moment = ( {:>9.4f} , {:>9.4f} , {:>9.4f} ) Total = {:>9.4f} Debye\n\n".format(
|
||||||
|
chg_dip[1], chg_dip[2], chg_dip[3], chg_dip[4]))
|
||||||
|
|
||||||
|
def calculate_step(self, fh):
|
||||||
|
|
||||||
|
invhessian = linalg.inv(self.hessian)
|
||||||
|
pre_step = -1 * np.matmul(invhessian, self.gradient.T).T
|
||||||
|
maxstep = np.amax(np.absolute(pre_step))
|
||||||
|
factor = min(1, player['maxstep']/maxstep)
|
||||||
|
step = factor * pre_step
|
||||||
|
|
||||||
|
fh.write("\nCalculated step:\n")
|
||||||
|
pre_step_list = pre_step.tolist()
|
||||||
|
|
||||||
|
fh.write("-----------------------------------------------------------------------\n"
|
||||||
|
"Center Atomic Step (Bohr)\n"
|
||||||
|
"Number Number X Y Z\n"
|
||||||
|
"-----------------------------------------------------------------------\n")
|
||||||
|
for i in range(len(molecules[0])):
|
||||||
|
fh.write(" {:>5d} {:>3d} {:>14.9f} {:>14.9f} {:>14.9f}\n".format(
|
||||||
|
i + 1, molecules[0][i]['na'],
|
||||||
|
pre_step_list.pop(0), pre_step_list.pop(0), pre_step_list.pop(0)))
|
||||||
|
|
||||||
|
fh.write("-----------------------------------------------------------------------\n")
|
||||||
|
|
||||||
|
fh.write("Maximum step is {:>11.6}\n".format(maxstep))
|
||||||
|
fh.write("Scaling factor = {:>6.4f}\n".format(factor))
|
||||||
|
fh.write("\nFinal step (Bohr):\n")
|
||||||
|
step_list = step.tolist()
|
||||||
|
|
||||||
|
fh.write("-----------------------------------------------------------------------\n"
|
||||||
|
"Center Atomic Step (Bohr)\n"
|
||||||
|
"Number Number X Y Z\n"
|
||||||
|
"-----------------------------------------------------------------------\n")
|
||||||
|
for i in range(len(molecules[0])):
|
||||||
|
fh.write(" {:>5d} {:>3d} {:>14.9f} {:>14.9f} {:>14.9f}\n".format(
|
||||||
|
i + 1, molecules[0][i]['na'],
|
||||||
|
step_list.pop(0), step_list.pop(0), step_list.pop(0)))
|
||||||
|
|
||||||
|
fh.write("-----------------------------------------------------------------------\n")
|
||||||
|
|
||||||
|
step_max = np.amax(np.absolute(step))
|
||||||
|
step_rms = np.sqrt(np.mean(np.square(step)))
|
||||||
|
|
||||||
|
fh.write(" Max Step = {:>14.9f} RMS Step = {:>14.9f}\n\n".format(
|
||||||
|
step_max, step_rms))
|
||||||
|
|
||||||
|
return step
|
||||||
|
|
||||||
class Atom:
|
class Atom:
|
||||||
|
|
||||||
def __init__(self, lbl,na,rx,ry,rz,chg,eps,sig):
|
def __init__(self, lbl,na,rx,ry,rz,chg,eps,sig):
|
||||||
|
|||||||
Reference in New Issue
Block a user