HideyoshiNakazone/DicePlayer
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Implements More Tests and Begins Dice Refactor Implementation

Browse Source
This commit is contained in:
Vitor Hideyoshi
2023-04-22 11:27:20 -03:00
parent 71744641ff
commit b869ee74fb
19 changed files with 719 additions and 297 deletions
Download Patch File Download Diff File Expand all files Collapse all files

121
diceplayer/shared/environment/molecule.py
View File

@@ -79,14 +79,13 @@ class Molecule:
"""
Updated positions based on the center of mass of the molecule
"""
self.center_of_mass()
for atom in self.atom:
atom.rx -= self.com[0]
atom.ry -= self.com[1]
atom.rz -= self.com[2]
self.center_of_mass()
def charges_and_dipole(self) -> List[float]:
"""
Calculates the charges and dipole of the molecule atoms
@@ -119,16 +118,15 @@ class Molecule:
distances = []
dim = len(self.atom)
for atom1 in self.atom:
if atom1.na != ghost_number:
for atom2 in self.atom:
if atom2.na != ghost_number:
dx = atom1.rx - atom2.rx
dy = atom1.ry - atom2.ry
dz = atom1.rz - atom2.rz
distances.append(math.sqrt(dx ** 2 + dy ** 2 + dz ** 2))
for index1, atom1 in enumerate(self.atom):
for index2, atom2 in enumerate(self.atom):
if index1 != index2:
dx = atom1.rx - atom2.rx
dy = atom1.ry - atom2.ry
dz = atom1.rz - atom2.rz
distances.append(math.sqrt(dx ** 2 + dy ** 2 + dz ** 2))
return np.array(distances).reshape(dim, dim)
return np.array(distances).reshape(dim, dim-1)
def inertia_tensor(self) -> NDArray[Shape["3, 3"], Float]:
"""
@@ -156,40 +154,6 @@ class Molecule:
return np.array([[Ixx, Ixy, Ixz], [Ixy, Iyy, Iyz], [Ixz, Iyz, Izz]])
def axes(self) -> NDArray[Shape["3, 3"], Float]:
"""
Calculates the axes of the molecule
Returns:
NDArray[Shape["3, 3"], Float]: Returns the axes of molecule
"""
eixos = np.zeros(3)
if len(self.atom) == 2:
position1 = np.array([self.atom[0].rx, self.atom[0].ry, self.atom[0].rz])
position2 = np.array([self.atom[1].rx, self.atom[1].ry, self.atom[1].rz])
eixos = position2 - position1
eixos /= linalg.norm(eixos)
elif len(self.atom) > 2:
position1 = np.array([self.atom[0].rx, self.atom[0].ry, self.atom[0].rz])
position2 = np.array([self.atom[1].rx, self.atom[1].ry, self.atom[1].rz])
position3 = np.array([self.atom[2].rx, self.atom[2].ry, self.atom[2].rz])
v1 = position2 - position1
v2 = position3 - position1
v3 = np.cross(v1, v2)
v2 = np.cross(v1, v3)
v1 /= linalg.norm(v1)
v2 /= linalg.norm(v2)
v3 /= linalg.norm(v3)
eixos = np.array(
[[v1[0], v1[1], v1[2]], [v2[0], v2[1], v2[2]], [v3[0], v3[1], v3[2]]]
)
return eixos
def principal_axes(self) -> Tuple[np.ndarray, np.ndarray]:
"""
Calculates the principal axes of the molecule
@@ -201,7 +165,7 @@ class Molecule:
try:
evals, evecs = linalg.eigh(self.inertia_tensor())
except:
except ValueError:
raise RuntimeError("Error: diagonalization of inertia tensor did not converge")
return evals, evecs
@@ -221,38 +185,38 @@ class Molecule:
return position
def updateCharges(self, charges: List[float]) -> None:
def update_charges(self, charges: List[float]) -> None:
for i, atom in enumerate(self.atom):
atom.chg = charges[i]
def update_hessian(
self,
step: np.ndarray,
cur_gradient: np.ndarray,
old_gradient: np.ndarray,
hessian: np.ndarray,
) -> np.ndarray:
"""
Updates the Hessian of the molecule based on the current hessian, the current gradient and the previous gradient
Args:
step (np.ndarray): step value of the iteration
cur_gradient (np.ndarray): current gradient
old_gradient (np.ndarray): previous gradient
hessian (np.ndarray): current hessian
Returns:
np.ndarray: updated hessian of the molecule
"""
dif_gradient = cur_gradient - old_gradient
mat1 = 1 / np.dot(dif_gradient, step) * np.matmul(dif_gradient.T, dif_gradient)
mat2 = 1 / np.dot(step, np.matmul(hessian, step.T).T)
mat2 *= np.matmul(np.matmul(hessian, step.T), np.matmul(step, hessian))
return hessian + mat1 - mat2
# @staticmethod
# def update_hessian(
# step: np.ndarray,
# cur_gradient: np.ndarray,
# old_gradient: np.ndarray,
# hessian: np.ndarray,
# ) -> np.ndarray:
# """
# Updates the Hessian of the molecule based on the current hessian, the current gradient and the previous gradient
#
# Args:
# step (np.ndarray): step value of the iteration
# cur_gradient (np.ndarray): current gradient
# old_gradient (np.ndarray): previous gradient
# hessian (np.ndarray): current hessian
#
# Returns:
# np.ndarray: updated hessian of the molecule
# """
#
# dif_gradient = cur_gradient - old_gradient
#
# mat1 = 1 / np.dot(dif_gradient, step) * np.matmul(dif_gradient.T, dif_gradient)
# mat2 = 1 / np.dot(step, np.matmul(hessian, step.T).T)
# mat2 *= np.matmul(np.matmul(hessian, step.T), np.matmul(step, hessian))
#
# return hessian + mat1 - mat2
def sizes_of_molecule(self) -> List[float]:
"""
@@ -267,10 +231,9 @@ class Molecule:
z_list = []
for atom in self.atom:
if atom.na != ghost_number:
x_list.append(atom.rx)
y_list.append(atom.ry)
z_list.append(atom.rz)
x_list.append(atom.rx)
y_list.append(atom.ry)
z_list.append(atom.rz)
x_max = max(x_list)
x_min = min(x_list)

194
diceplayer/shared/environment/system.py
View File

@@ -1,64 +1,47 @@
import math
from copy import deepcopy
from typing import List, Tuple, TextIO
import numpy as np
from numpy import linalg
from diceplayer.shared.environment.molecule import Molecule
from diceplayer.shared.utils.misc import BOHR2ANG
from diceplayer.shared.utils.ptable import atomsymb
from diceplayer.shared.utils.misc import BOHR2ANG
from typing import List, Tuple, TextIO
from copy import deepcopy
from numpy import linalg
import numpy as np
import math
class System:
"""
System class declaration. This class is used throughout the DicePlayer program to represent the system containing the molecules.
System class declaration. This class is used throughout the DicePlayer program to represent the system containing the molecules.
Atributes:
molecule (List[Molecule]): List of molecules of the system
nmols (List[int]): List of number of molecules in the system
"""
Atributes:
molecule (List[Molecule]): List of molecules of the system
nmols (List[int]): List of number of molecules in the system
"""
def __init__(self) -> None:
"""
Initializes a empty system object that will be populated afterwards
"""
Initializes a empty system object that will be populated afterwards
"""
self.molecule: List[Molecule] = []
self.nmols: List[int] = []
def add_type(self, nmols: int, m: Molecule) -> None:
"""
Adds a new molecule type to the system
Adds a new molecule type to the system
Args:
nmols (int): Number of molecules of the new type in the system
m (Molecule): The instance of the new type of molecule
"""
Args:
nmols (int): Number of molecules of the new type in the system
m (Molecule): The instance of the new type of molecule
"""
if isinstance(m, Molecule) is False:
raise TypeError("Error: molecule is not a Molecule instance")
self.molecule.append(m)
if isinstance(nmols, int) is False:
raise TypeError("Error: nmols is not an integer")
self.nmols.append(nmols)
def center_of_mass_distance(self, a: int, b: int) -> float:
"""
Calculates the distance between the center of mass of two molecules
Args:
a (Molecule): First Molecule Instance
b (Molecule): Second Molecule Instance
Returns:
float: module of the distance between the two center of masses
"""
com1 = self.molecule[a].center_of_mass()
com2 = self.molecule[b].center_of_mass()
dx = com1[0] - com2[0]
dy = com1[1] - com2[1]
dz = com1[2] - com2[2]
distance = math.sqrt(dx**2 + dy**2 + dz**2)
return distance
def rmsd_fit(self, p_index: int, r_index: int) -> Tuple[float, Molecule]:
projecting_mol = self.molecule[p_index]
@@ -118,9 +101,9 @@ class System:
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
(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)
@@ -135,13 +118,71 @@ class System:
return rmsd, projected_mol
# def center_of_mass_distance(self, a: int, b: int) -> float:
# """
# Calculates the distance between the center of mass of two molecules
#
# Args:
# a (Molecule): First Molecule Instance
# b (Molecule): Second Molecule Instance
#
# Returns:
# float: module of the distance between the two center of masses
# """
#
# com1 = self.molecule[a].center_of_mass()
# com2 = self.molecule[b].center_of_mass()
# dx = com1[0] - com2[0]
# dy = com1[1] - com2[1]
# dz = com1[2] - com2[2]
# distance = math.sqrt(dx**2 + dy**2 + dz**2)
#
# return distance
# def nearest_image(
# self,
# index_r: int,
# index_m: int,
# lx: float,
# ly: float,
# lz: float,
# criterium=None,
# ) -> Tuple[float, Molecule]:
#
# if criterium in None:
# criterium = "com"
#
# if criterium != "com" and criterium != "min":
# raise RuntimeError("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 update_molecule(self, position: np.ndarray, fh: TextIO) -> None:
"""Updates the position of the molecule in the Output file
Args:
position (np.ndarray): numpy position vector
fh (TextIO): Output file
"""
Args:
position (np.ndarray): numpy position vector
fh (TextIO): Output file
"""
position_in_ang = (position * BOHR2ANG).tolist()
self.add_type(self.nmols[0], deepcopy(self.molecule[0]))
@@ -158,52 +199,15 @@ class System:
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: int,
index_m: int,
lx: float,
ly: float,
lz: float,
criterium=None,
) -> Tuple[float, Molecule]:
if criterium in None:
criterium = "com"
if criterium != "com" and criterium != "min":
raise RuntimeError("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: int, fh: TextIO) -> None:
"""
Print the geometry of the molecule in the Output file
Print the geometry of the molecule in the Output file
Args:
cycle (int): Number of the cycle
fh (TextIO): Output file
"""
Args:
cycle (int): Number of the cycle
fh (TextIO): Output file
"""
fh.write("Cycle # {}\n".format(cycle))
fh.write("Number of site: {}\n".format(len(self.molecule[0].atom)))
@@ -217,12 +221,12 @@ class System:
def printChargesAndDipole(self, cycle: int, fh: TextIO) -> None:
"""
Print the charges and dipole of the molecule in the Output file
Print the charges and dipole of the molecule in the Output file
Args:
cycle (int): Number of the cycle
fh (TextIO): Output file
"""
Args:
cycle (int): Number of the cycle
fh (TextIO): Output file
"""
fh.write("Cycle # {}\n".format(cycle))
fh.write("Number of site: {}\n".format(len(self.molecule[0].atom)))