From 6017f272b335e9e987e71e0258bc7dd2cabcbd11 Mon Sep 17 00:00:00 2001 From: Vitor Hideyoshi Nakazone Batista Date: Sun, 24 Oct 2021 00:55:11 -0300 Subject: [PATCH] SetGlobas/MolHandling(fixes)/DicePlayer Translation Fixed Numerus functions in the SetGlobals file and MolHandling, translated DicePlayer and tested till, such as inertia_tensor, add_type, center_of_mass_to_origin. standard_orientation; and removed chglevel keyword Signed-off-by: Vitor Hideyoshi --- DPpack/MolHandling.py | 68 ++++++++++++++++++++----------------------- DPpack/SetGlobals.py | 32 +++++++++----------- diceplayer.py | 7 +++-- run.log | 25 ++++++++-------- run.log.backup | 9 ++---- 5 files changed, 64 insertions(+), 77 deletions(-) diff --git a/DPpack/MolHandling.py b/DPpack/MolHandling.py index 6f5d78e..1951233 100644 --- a/DPpack/MolHandling.py +++ b/DPpack/MolHandling.py @@ -28,7 +28,7 @@ class System: def add_type(self,nmols, m): self.molecule.append(m) - self.nmols = nmols + self.nmols.append(nmols) # Função que calcula a distância entre dois centros de massa # e por se tratar de uma função de dois atomos não deve ser @@ -208,21 +208,17 @@ class Molecule: def center_of_mass(self): - com = np.zeros(3) - total_mass = 0.0 + self.com = np.zeros(3) for atom in self.atom: - total_mass += atom.mass - com += atom.mass * np.array([atom.rx, atom.ry, atom.rz]) + self.com += atom.mass * np.array([atom.rx, atom.ry, atom.rz]) - com = com / total_mass - - self.com = com + self.com = self.com / self.total_mass def center_of_mass_to_origin(self): - com = self.center_of_mass() + self.center_of_mass() for atom in self.atom: @@ -260,6 +256,30 @@ class Molecule: return np.array(distances).reshape(dim, dim) + def inertia_tensor(self): + + self.center_of_mass() + Ixx = Ixy = Ixz = Iyy = Iyz = Izz = 0.0 + + for atom in self.atom: + + #### Obtain the displacement from the center of mass + dx = atom.rx - self.com[0] + dy = atom.ry - self.com[1] + dz = atom.rz - self.com[2] + #### Update the diagonal components of the tensor + Ixx += atom.mass * (dy**2 + dz**2) + Iyy += atom.mass * (dz**2 + dx**2) + Izz += atom.mass * (dx**2 + dy**2) + #### Update the off-diagonal components of the tensor + Ixy += atom.mass * dx * dy * -1 + Ixz += atom.mass * dx * dz * -1 + Iyz += atom.mass * dy * dz * -1 + + return np.array([[Ixx, Ixy, Ixz], + [Ixy, Iyy, Iyz], + [Ixz, Iyz, Izz]]) + def eixos(self): eixos = np.zeros(3) @@ -288,36 +308,13 @@ class Molecule: return eixos - def inertia_tensor(self): - - Ixx = Ixy = Ixz = Iyy = Iyz = Izz = 0.0 - - for atom in self.atom: - - #### Obtain the displacement from the center of mass - dx = atom.rx - self.com[0] - dy = atom.ry - self.com[1] - dz = atom.rz - self.com[2] - #### Update the diagonal components of the tensor - Ixx += atom.mass * (dy**2 + dz**2) - Iyy += atom.mass * (dz**2 + dx**2) - Izz += atom.mass * (dx**2 + dy**2) - #### Update the off-diagonal components of the tensor - Ixy += atom.mass * dx * dy * -1 - Ixz += atom.mass * dx * dz * -1 - Iyz += atom.mass * dy * dz * -1 - - return np.array([ [Ixx, Ixy, Ixz], - [Ixy, Iyy, Iyz], - [Ixz, Iyz, Izz] ]) - def principal_axes(self): try: evals, evecs = linalg.eigh(self.inertia_tensor()) except: sys.exit("Error: diagonalization of inertia tensor did not converge") - + return evals, evecs def read_position(self): @@ -366,7 +363,6 @@ class Molecule: def standard_orientation(self): self.center_of_mass_to_origin() - tensor = self.inertia_tensor() evals, evecs = self.principal_axes() if round(linalg.det(evecs)) == -1: @@ -435,6 +431,4 @@ class Atom: self.chg = chg # Double self.eps = eps # Double self.sig = sig # Double - self.mass = atommass[self.na] # Double - - + self.mass = atommass[self.na] # Double \ No newline at end of file diff --git a/DPpack/SetGlobals.py b/DPpack/SetGlobals.py index 6fd3c29..5fb54b4 100644 --- a/DPpack/SetGlobals.py +++ b/DPpack/SetGlobals.py @@ -1,6 +1,7 @@ import os, sys import shutil import textwrap +import types from DPpack.MolHandling import * from DPpack.PTable import * @@ -179,8 +180,8 @@ class Internal: except ValueError: sys.exit(err) - elif key in ('level', 'chglevel'): - setattr(self.gaussian, key, value[0]) + elif key == 'level': + setattr(self.gaussian, key, value) elif key in ('gmiddle', 'gbottom'): setattr(self.gaussian, key, value[0]) @@ -243,9 +244,6 @@ class Internal: if self.gaussian.pop != "chelpg" and (self.player.ghosts == "yes" or self.player.lps == "yes"): sys.exit("Error: ghost atoms or lone pairs only available with 'pop = chelpg')") - - if self.gaussian.chglevel == None: - self.gaussian.chglevel = self.gaussian.level ## Check only if QM program is Molcas: # if self.player.qmprog == "molcas": @@ -636,7 +634,7 @@ class Internal: pass - + ### I still have to talk with Herbet about def populate_asec_vdw(self, cycle): asec_charges = [] # (rx, ry, rz, chg) @@ -651,17 +649,17 @@ class Internal: nsitesref = len(self.system.molecule[0]) + len(ghost_atoms) + len(lp_atoms) - nsites_total = dice['nmol'][0] * nsitesref - for i in range(1, len(dice['nmol'])): - nsites_total += dice['nmol'][i] * len(molecules[i]) + nsites_total = self.dice.nmol[0] * nsitesref + for i in range(1, len(self.dice.nmol)): + nsites_total += self.dice.nmol[i] * len(self.system.molecule[i]) thickness = [] picked_mols = [] - for proc in range(1, player['nprocs'] + 1): ## Run over folders + for proc in range(1, self.player.nprocs + 1): ## Run over folders path = "step{:02d}".format(cycle) + os.sep + "p{:02d}".format(proc) - file = path + os.sep + dice['outname'] + ".xyz" + file = path + os.sep + self.dice.outname + ".xyz" if not os.path.isfile(file): sys.exit("Error: cannot find file {}".format(file)) try: @@ -677,23 +675,23 @@ class Internal: box = xyzfile.pop(0).split()[-3:] box = [ float(box[0]), float(box[1]), float(box[2]) ] - sizes = sizes_of_molecule(molecules[0]) + sizes = self.system.molecule[0].sizes_of_molecule() thickness.append( min([ (box[0] - sizes[0])/2, (box[1] - sizes[1])/2, (box[2] - sizes[2])/2 ]) ) xyzfile = xyzfile[nsitesref:] ## Skip the first (reference) molecule mol_count = 0 - for type in range(len(dice['nmol'])): ## Run over types of molecules + for type in range(len(self.dice.nmol)): ## Run over types of molecules if type == 0: - nmols = dice['nmol'][0] - 1 + nmols = self.dice.nmol[0] - 1 else: - nmols = dice['nmol'][type] + nmols = self.dice.nmol[type] for mol in range(nmols): ## Run over molecules of each type new_molecule = [] - for site in range(len(molecules[type])): ## Run over sites of each molecule + for site in range(len(self.system.molecule[types].atoms)): ## Run over sites of each molecule new_molecule.append({}) line = xyzfile.pop(0).split() @@ -831,8 +829,6 @@ class Internal: self.gmiddle = None # In each case, if a filename is given, its content will be self.gbottom = None # inserted in the gaussian input self.pop = "chelpg" - self.chglevel = None - self.level = None # class Molcas: diff --git a/diceplayer.py b/diceplayer.py index 821fdc9..4c64081 100644 --- a/diceplayer.py +++ b/diceplayer.py @@ -95,6 +95,7 @@ if __name__ == '__main__': #### Bring the molecules to standard orientation and prints info about them for i in range(len(internal.system.molecule)): + internal.outfile.write("\nMolecule type {}:\n\n".format(i + 1)) internal.system.molecule[i].print_mol_info(internal.outfile) internal.outfile.write(" Translating and rotating molecule to standard orientation...") @@ -134,9 +135,9 @@ if __name__ == '__main__': # #Molcas.read_forces("grad_hessian.dat") # #Molcas.read_hessian("grad_hessian.dat") - # #### - # #### Start the iterative process - # #### + #### + #### Start the iterative process + #### # for cycle in range(internal.player.cyc, internal.player.cyc + internal.player.maxcyc): diff --git a/run.log b/run.log index 45c1227..08780f9 100644 --- a/run.log +++ b/run.log @@ -2,10 +2,10 @@ ############# Welcome to DICEPLAYER version 1.0 ############# ########################################################################################## -Your python version is 3.8.8 (default, Apr 13 2021, 19:58:26) -[GCC 7.3.0] +Your python version is 3.8.12 (default, Oct 7 2021, 05:40:48) +[GCC 11.1.0] -Program started on Saturday, 25 Sep 2021 at 15:24:31 +Program started on Sunday, 24 Oct 2021 at 00:46:13 Environment variables: OMP_STACKSIZE = Not set @@ -50,7 +50,6 @@ title = Diceplayer run GAUSSIAN variables being used in this run: ------------------------------------------------------------------------------------------ -chglevel = MP2/aug-cc-pVTZ chgmult = 0 1 level = MP2/aug-cc-pVTZ pop = chelpg @@ -139,11 +138,11 @@ Molecule type 1: Translating and rotating molecule to standard orientation... Done New values: - Center of mass = ( -0.0000 , 0.0000 , 0.0000 ) + Center of mass = ( -0.0000 , -0.0000 , -0.0000 ) Moments of inertia = 3.144054E+02 2.801666E+03 3.027366E+03 Major principal axis = ( 1.000000 , 0.000000 , 0.000000 ) Inter principal axis = ( -0.000000 , 1.000000 , -0.000000 ) - Minor principal axis = ( -0.000000 , 0.000000 , 1.000000 ) + Minor principal axis = ( 0.000000 , 0.000000 , 1.000000 ) Characteristic lengths = ( 11.97 , 5.27 , 2.99 ) Total mass = 226.28 au Total charge = -0.0000 e @@ -165,14 +164,14 @@ Molecule type 2: Translating and rotating molecule to standard orientation... Done New values: - Center of mass = ( 1.6067 , -1.0929 , 0.0026 ) - Moments of inertia = 1.561638E+02 6.739391E+02 8.270247E+02 - Major principal axis = ( 0.899747 , -0.436411 , 0.000541 ) - Inter principal axis = ( 0.436411 , 0.899747 , -0.001219 ) - Minor principal axis = ( 0.000045 , 0.001333 , 0.999999 ) - Characteristic lengths = ( 5.67 , 5.05 , 1.51 ) + Center of mass = ( 0.0000 , -0.0000 , 0.0000 ) + Moments of inertia = 1.205279E+02 2.859254E+02 4.033761E+02 + Major principal axis = ( 1.000000 , -0.000000 , 0.000000 ) + Inter principal axis = ( 0.000000 , 1.000000 , -0.000000 ) + Minor principal axis = ( -0.000000 , 0.000000 , 1.000000 ) + Characteristic lengths = ( 5.67 , 5.13 , 1.51 ) Total mass = 112.20 au Total charge = -0.0000 e - Dipole moment = ( 1.8148 , 1.4689 , -0.0016 ) Total = 2.3347 Debye + Dipole moment = ( 1.7575 , 1.5369 , -0.0000 ) Total = 2.3347 Debye ========================================================================================== diff --git a/run.log.backup b/run.log.backup index 1f9ac41..cd56ae8 100644 --- a/run.log.backup +++ b/run.log.backup @@ -2,10 +2,10 @@ ############# Welcome to DICEPLAYER version 1.0 ############# ########################################################################################## -Your python version is 3.8.8 (default, Apr 13 2021, 19:58:26) -[GCC 7.3.0] +Your python version is 3.8.12 (default, Oct 7 2021, 05:40:48) +[GCC 11.1.0] -Program started on Saturday, 25 Sep 2021 at 15:23:44 +Program started on Sunday, 24 Oct 2021 at 00:44:02 Environment variables: OMP_STACKSIZE = Not set @@ -50,7 +50,6 @@ title = Diceplayer run GAUSSIAN variables being used in this run: ------------------------------------------------------------------------------------------ -chglevel = MP2/aug-cc-pVTZ chgmult = 0 1 level = MP2/aug-cc-pVTZ pop = chelpg @@ -176,5 +175,3 @@ Molecule type 2: Dipole moment = ( 1.8148 , 1.4689 , -0.0016 ) Total = 2.3347 Debye ========================================================================================== - -Starting the iterative process.