From 316acf47613052b7cdb5d4e702866b433f1c7bdf Mon Sep 17 00:00:00 2001 From: Vitor Hideyoshi Nakazone Batista Date: Sat, 25 Sep 2021 15:30:40 -0300 Subject: [PATCH] SetGlobas/MolHandling(fixes)/DicePlayer Translation Fixed Numerus functions in the SetGlobals file and MolHandling, translated DicePlayer and tested till printing of initial geometry for cyc==1 Signed-off-by: Vitor Hideyoshi --- DPpack/MolHandling.py | 84 +--- DPpack/SetGlobals.py | 380 ++++++++++++----- DPpack/__pycache__/Dice.cpython-38.pyc | Bin 0 -> 11245 bytes DPpack/__pycache__/Gaussian.cpython-38.pyc | Bin 0 -> 9011 bytes DPpack/__pycache__/Misc.cpython-38.pyc | Bin 0 -> 1627 bytes DPpack/__pycache__/MolHandling.cpython-38.pyc | Bin 0 -> 11684 bytes DPpack/__pycache__/PTable.cpython-38.pyc | Bin 0 -> 2017 bytes DPpack/__pycache__/SetGlobals.cpython-38.pyc | Bin 0 -> 19799 bytes DPpack/__pycache__/__init__.cpython-38.pyc | Bin 0 -> 155 bytes control.in | 8 +- diceplayer.py | 399 ++++++++++-------- run.log | 97 +---- run.log.backup | 180 ++++++++ 13 files changed, 725 insertions(+), 423 deletions(-) create mode 100644 DPpack/__pycache__/Dice.cpython-38.pyc create mode 100644 DPpack/__pycache__/Gaussian.cpython-38.pyc create mode 100644 DPpack/__pycache__/Misc.cpython-38.pyc create mode 100644 DPpack/__pycache__/MolHandling.cpython-38.pyc create mode 100644 DPpack/__pycache__/PTable.cpython-38.pyc create mode 100644 DPpack/__pycache__/SetGlobals.cpython-38.pyc create mode 100644 DPpack/__pycache__/__init__.cpython-38.pyc create mode 100644 run.log.backup diff --git a/DPpack/MolHandling.py b/DPpack/MolHandling.py index d47d817..6f5d78e 100644 --- a/DPpack/MolHandling.py +++ b/DPpack/MolHandling.py @@ -1,4 +1,3 @@ -from DPpack.MolHandling import total_mass import os, sys import math import shutil @@ -11,7 +10,11 @@ from numpy import linalg from DPpack.Misc import * from DPpack.PTable import * -from DPpack.SetGlobals import * + +env = ["OMP_STACKSIZE"] + +bohr2ang = 0.52917721092 +ang2bohr = 1/bohr2ang # Usaremos uma nova classe que ira conter toda interação entre moleculas @@ -173,9 +176,9 @@ class System: def print_geom(self, cycle, fh): - fh.write("{}\n".format(len(self.molecule[0]))) + fh.write("{}\n".format(len(self.molecule[0].atom))) fh.write("Cycle # {}\n".format(cycle)) - for atom in self.molecule[0].atoms: + for atom in self.molecule[0].atom: symbol = atomsymb[atom.na] fh.write("{:<2s} {:>10.6f} {:>10.6f} {:>10.6f}\n".format(symbol, atom.rx, atom.ry, atom.rz)) @@ -201,6 +204,8 @@ class Molecule: self.atom.append(a) # Inserção de um novo atomo self.total_mass += a.mass + self.center_of_mass() + def center_of_mass(self): com = np.zeros(3) @@ -221,9 +226,9 @@ class Molecule: for atom in self.atom: - atom.rx -= com[0] - atom.ry -= com[1] - atom.rz -= com[2] + atom.rx -= self.com[0] + atom.ry -= self.com[1] + atom.rz -= self.com[2] def charges_and_dipole(self): @@ -285,15 +290,14 @@ class Molecule: def inertia_tensor(self): - com = 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 - com[0] - dy = atom.ry - com[1] - dz = atom.rz - com[2] + 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) @@ -332,7 +336,7 @@ class Molecule: mat1 = 1/np.dot(dif_gradient, step) * np.matmul(dif_gradient.T, dif_gradient) mat2 = 1/np.dot(step, np.matmul(self.hessian, step.T).T) - mat2 *= np.matmul( np.matmul(self.hessian, step.T), np.matmul(step, hessian) ) + mat2 *= np.matmul( np.matmul(self.hessian, step.T), np.matmul(step, self.hessian) ) self.hessian += mat1 - mat2 @@ -393,10 +397,9 @@ class 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])) + + fh.write(" Center of mass = ( {:>10.4f} , {:>10.4f} , {:>10.4f} )\n".format(self.com[0], + self.com[1], self.com[2])) inertia = self.inertia_tensor() evals, evecs = self.principal_axes() @@ -413,60 +416,13 @@ class Molecule: 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)) + fh.write(" Total mass = {:>8.2f} au\n".format(self.total_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: def __init__(self, lbl,na,rx,ry,rz,chg,eps,sig): diff --git a/DPpack/SetGlobals.py b/DPpack/SetGlobals.py index e582316..6fd3c29 100644 --- a/DPpack/SetGlobals.py +++ b/DPpack/SetGlobals.py @@ -6,16 +6,13 @@ from DPpack.MolHandling import * from DPpack.PTable import * from DPpack.Misc import * -env = ["OMP_STACKSIZE"] - -bohr2ang = 0.52917721092 -ang2bohr = 1/bohr2ang - class Internal: - def __init__(self, infile): + def __init__(self, infile, outfile): self.infile = infile + self.outfile = outfile + self.system = System() self.player = self.Player() @@ -45,10 +42,9 @@ class Internal: def read_keywords(self): try: - with open(self.infile) as fh: - controlfile = fh.readlines() + controlfile = self.infile.readlines() except EnvironmentError: - sys.exit("Error: cannot open file {}".format(self.infile)) + sys.exit("Error: cannot read file {}".format(self.infile)) for line in controlfile: @@ -77,7 +73,7 @@ class Internal: elif key in ('readhessian', 'vdwforces') and value[0].lower() in ("yes", "no"): setattr(self.player, key, value[0].lower()) - elif key in ('maxcyc', 'initcyc', 'nprocs', 'altsteps', 'switchcyc'): + elif key in ('maxcyc', 'nprocs', 'altsteps', 'switchcyc'): err = "Error: expected a positive integer for keyword {} in file {}".format(key, self.infile) try: new_value = int(value[0]) @@ -95,7 +91,7 @@ class Internal: if new_value < 0.01: sys.exit(err) else: - setattr(self.player, key).append(new_value) + setattr(self.player, key, new_value) except ValueError: sys.exit(err) @@ -136,7 +132,7 @@ class Internal: if new_value < 1: sys.exit(err) else: - setattr(self.dice, key, new_value) + getattr(self.dice, key).append(new_value) elif i == 0: sys.exit(err) else: @@ -153,7 +149,7 @@ class Internal: if new_value < 1: sys.exit(err) else: - setattr(self.dice, key, new_value) + getattr(self.dice, key).append(new_value) elif i < 2: sys.exit(err) else: @@ -179,7 +175,7 @@ class Internal: sys.exit(err) for i in range (2): try: - setattr(self.gaussian, key)[i] = int(value[i]) + getattr(self.gaussian, key)[i] = int(value[i]) except ValueError: sys.exit(err) @@ -192,24 +188,24 @@ class Internal: elif key == 'pop' and value[0].lower() in ("chelpg", "mk", "nbo"): setattr(self.gaussian, key, value[0].lower()) - #### Read the Molcas related keywords - elif key in self.molcas_keywords and len(value) != 0: ## 'value' is not empty! + # #### Read the Molcas related keywords + # elif key in self.molcas_keywords and len(value) != 0: ## 'value' is not empty! - if key == 'root': # If defined, must be well defined (only positive integer values) - err = "Error: expected a positive integer for keyword {} in file {}".format(key, self.infile) - if not value[0].isdigit(): - sys.exit(err) - new_value = int(value[0]) - if new_value >= 1: - setattr(self.molcas, key, new_value) + # if key == 'root': # If defined, must be well defined (only positive integer values) + # err = "Error: expected a positive integer for keyword {} in file {}".format(key, self.infile) + # if not value[0].isdigit(): + # sys.exit(err) + # new_value = int(value[0]) + # if new_value >= 1: + # setattr(self.molcas, key, new_value) - elif key in ('mbottom', 'orbfile'): - setattr(self.molcas, key, value[0]) + # elif key in ('mbottom', 'orbfile'): + # setattr(self.molcas, key, value[0]) - elif key == 'basis': - setattr(self.molcas ,key, value[0]) + # elif key == 'basis': + # setattr(self.molcas ,key, value[0]) - #### End + # #### End def check_keywords(self): @@ -225,10 +221,10 @@ class Internal: if self.dice.dens == None: sys.exit("Error: 'dens' keyword not specified in file {}".format(self.infile)) - if len(self.dice.nmol) == 0: + if self.dice.nmol == 0: sys.exit("Error: 'nmol' keyword not defined appropriately in file {}".format(self.infile)) - if len(self.dice.nstep) == 0: + if self.dice.nstep == 0: sys.exit("Error: 'nstep' keyword not defined appropriately in file {}".format(self.infile)) ## Check only if QM program is Gaussian: @@ -290,22 +286,22 @@ class Internal: # isave value is always the nearest multiple of 100 self.dice.isave = round(self.dice.isave / 100) * 100 - def print_keywords(self, fh): + def print_keywords(self): - fh.write("##########################################################################################\n" + self.outfile.write("##########################################################################################\n" "############# Welcome to DICEPLAYER version 1.0 #############\n" "##########################################################################################\n" "\n") - fh.write("Your python version is {}\n".format(sys.version)) - fh.write("\n") - fh.write("Program started on {}\n".format(weekday_date_time())) - fh.write("\n") - fh.write("Environment variables:\n") + self.outfile.write("Your python version is {}\n".format(sys.version)) + self.outfile.write("\n") + self.outfile.write("Program started on {}\n".format(weekday_date_time())) + self.outfile.write("\n") + self.outfile.write("Environment variables:\n") for var in env: - fh.write("{} = {}\n".format(var, + self.outfile.write("{} = {}\n".format(var, (os.environ[var] if var in os.environ else "Not set"))) - fh.write("\n==========================================================================================\n" + self.outfile.write("\n==========================================================================================\n" " CONTROL variables being used in this run:\n" "------------------------------------------------------------------------------------------\n" "\n") @@ -314,13 +310,13 @@ class Internal: if getattr(self.player,key) != None: if isinstance(getattr(self.player,key), list): string = " ".join(str(x) for x in getattr(self.player,key)) - fh.write("{} = {}\n".format(key, string)) + self.outfile.write("{} = {}\n".format(key, string)) else: - fh.write("{} = {}\n".format(key, getattr(self.player,key))) + self.outfile.write("{} = {}\n".format(key, getattr(self.player,key))) - fh.write("\n") + self.outfile.write("\n") - fh.write("------------------------------------------------------------------------------------------\n" + self.outfile.write("------------------------------------------------------------------------------------------\n" " DICE variables being used in this run:\n" "------------------------------------------------------------------------------------------\n" "\n") @@ -329,15 +325,15 @@ class Internal: if getattr(self.dice,key) != None: if isinstance(getattr(self.dice,key), list): string = " ".join(str(x) for x in getattr(self.dice,key)) - fh.write("{} = {}\n".format(key, string)) + self.outfile.write("{} = {}\n".format(key, string)) else: - fh.write("{} = {}\n".format(key, getattr(self.dice,key))) + self.outfile.write("{} = {}\n".format(key, getattr(self.dice,key))) - fh.write("\n") + self.outfile.write("\n") if self.player.qmprog in ("g03", "g09", "g16"): - fh.write("------------------------------------------------------------------------------------------\n" + self.outfile.write("------------------------------------------------------------------------------------------\n" " GAUSSIAN variables being used in this run:\n" "------------------------------------------------------------------------------------------\n" "\n") @@ -346,15 +342,15 @@ class Internal: if getattr(self.gaussian,key) != None: if isinstance(getattr(self.gaussian,key), list): string = " ".join(str(x) for x in getattr(self.gaussian,key)) - fh.write("{} = {}\n".format(key, string)) + self.outfile.write("{} = {}\n".format(key, string)) else: - fh.write("{} = {}\n".format(key, getattr(self.gaussian,key))) + self.outfile.write("{} = {}\n".format(key, getattr(self.gaussian,key))) - fh.write("\n") + self.outfile.write("\n") # elif self.player.qmprog == "molcas": - # fh.write("------------------------------------------------------------------------------------------\n" + # self.outfile.write("------------------------------------------------------------------------------------------\n" # " MOLCAS variables being used in this run:\n" # "------------------------------------------------------------------------------------------\n" # "\n") @@ -363,11 +359,11 @@ class Internal: # if molcas[key] != None: # if isinstance(molcas[key], list): # string = " ".join(str(x) for x in molcas[key]) - # fh.write("{} = {}\n".format(key, string)) + # self.outfile.write("{} = {}\n".format(key, string)) # else: - # fh.write("{} = {}\n".format(key, molcas[key])) + # self.outfile.write("{} = {}\n".format(key, molcas[key])) - # fh.write("\n") + # self.outfile.write("\n") def read_potential(self): # Deve ser atualizado para o uso de @@ -390,7 +386,6 @@ class Internal: if ntypes != len(self.dice.nmol): sys.exit("Error: number of molecule types in file {} must match that of 'nmol' keyword in file {}".format( self.dice.ljname, self.infile)) - line = 2 for i in range(ntypes): @@ -400,7 +395,7 @@ class Internal: sys.exit("Error: expected an integer in line {} of file {}".format(line, self.dice.ljname)) nsites = int(nsites) - self.system.add_type(Molecule()) + self.system.add_type(nsites,Molecule()) for j in range(nsites): @@ -410,8 +405,6 @@ class Internal: if len(new_atom) < 8: sys.exit("Error: expected at least 8 fields in line {} of file {}".format(line, self.dice.ljname)) - self.system.molecule[i].add_atom() - if not new_atom[0].isdigit(): sys.exit("Error: expected an integer in field 1, line {} of file {}".format(line, self.dice.ljname)) lbl = int(new_atom[0]) @@ -468,7 +461,7 @@ class Internal: "Error: expected a positive float after 'mass=' in field 9, line {} of file {}".format( line, self.dice.ljname)) - self.system.molecule[i].add_atom(Atom(lbl,na,rx,ry,rz,chg,eps,sig,mass)) + self.system.molecule[i].add_atom(Atom(lbl,na,rx,ry,rz,chg,eps,sig)) to_delete = ['lbl','na','rx','ry','rz','chg','eps','sig','mass'] for _var in to_delete: @@ -476,44 +469,44 @@ class Internal: exec(f'del {_var}') - def print_potential(self, fh): + def print_potential(self): formatstr = "{:<3d} {:>3d} {:>10.5f} {:>10.5f} {:>10.5f} {:>10.6f} {:>9.5f} {:>7.4f} {:>9.4f}\n" - fh.write("\n" + self.outfile.write("\n" "==========================================================================================\n") - fh.write(" Potential parameters from file {}:\n".format(self.dice.ljname)) - fh.write("------------------------------------------------------------------------------------------\n" + self.outfile.write(" Potential parameters from file {}:\n".format(self.dice.ljname)) + self.outfile.write("------------------------------------------------------------------------------------------\n" "\n") - fh.write("Combination rule: {}\n".format(self.dice.combrule)) - fh.write("Types of molecules: {}\n\n".format(len(self.system.molecule))) + self.outfile.write("Combination rule: {}\n".format(self.dice.combrule)) + self.outfile.write("Types of molecules: {}\n\n".format(len(self.system.molecule))) i = 0 for mol in self.system.molecule: i += 1 - fh.write("{} atoms in molecule type {}:\n".format(len(mol), i)) - fh.write("---------------------------------------------------------------------------------\n" + self.outfile.write("{} atoms in molecule type {}:\n".format(len(mol.atom), i)) + self.outfile.write("---------------------------------------------------------------------------------\n" "Lbl AN X Y Z Charge Epsilon Sigma Mass\n") - fh.write("---------------------------------------------------------------------------------\n") + self.outfile.write("---------------------------------------------------------------------------------\n") for atom in mol.atom: - fh.write(formatstr.format(atom.lbl, atom.na, atom.rx, atom.ry, atom.rz, + self.outfile.write(formatstr.format(atom.lbl, atom.na, atom.rx, atom.ry, atom.rz, atom.chg, atom.eps, atom.sig, atom.mass)) - fh.write("\n") + self.outfile.write("\n") if self.player.ghosts == "yes" or self.player.lps == "yes": - fh.write("\n" + self.outfile.write("\n" "------------------------------------------------------------------------------------------\n" " Aditional potential parameters:\n" "------------------------------------------------------------------------------------------\n") # if player['ghosts'] == "yes": - # fh.write("\n") - # fh.write("{} ghost atoms appended to molecule type 1 at:\n".format(len(ghost_types))) - # fh.write("---------------------------------------------------------------------------------\n") + # self.outfile.write("\n") + # self.outfile.write("{} ghost atoms appended to molecule type 1 at:\n".format(len(ghost_types))) + # self.outfile.write("---------------------------------------------------------------------------------\n") # atoms_string = "" # for ghost in ghost_types: @@ -522,19 +515,19 @@ class Internal: # atoms_string += "{}{} ".format(atom_sym,atom) # if ghost['type'] == "g": - # fh.write(textwrap.fill("* Geometric center of atoms {}".format(atoms_string), 80)) + # self.outfile.write(textwrap.fill("* Geometric center of atoms {}".format(atoms_string), 80)) # elif ghost['type'] == "m": - # fh.write(textwrap.fill("* Center of mass of atoms {}".format(atoms_string), 80)) + # self.outfile.write(textwrap.fill("* Center of mass of atoms {}".format(atoms_string), 80)) # elif ghost['type'] == "z": - # fh.write(textwrap.fill("* Center of atomic number of atoms {}".format(atoms_string), 80)) + # self.outfile.write(textwrap.fill("* Center of atomic number of atoms {}".format(atoms_string), 80)) - # fh.write("\n") + # self.outfile.write("\n") # if player['lps'] == 'yes': - # fh.write("\n") - # fh.write("{} lone pairs appended to molecule type 1:\n".format(len(lp_types))) - # fh.write("---------------------------------------------------------------------------------\n") + # self.outfile.write("\n") + # self.outfile.write("{} lone pairs appended to molecule type 1:\n".format(len(lp_types))) + # self.outfile.write("---------------------------------------------------------------------------------\n") # for lp in lp_types: # # LP type 1 or 2 @@ -546,36 +539,33 @@ class Internal: # atom3_num = lp['numbers'][2] # atom3_sym = atomsymb[ molecules[0][atom3_num - 1]['na'] ].strip() - # fh.write(textwrap.fill( + # self.outfile.write(textwrap.fill( # "* Type {} on atom {}{} with {}{} {}{}. Alpha = {:<5.1f} Deg and D = {:<4.2f} Angs".format( # lp['type'], atom1_sym, atom1_num, atom2_sym, atom2_num, atom3_sym, atom3_num, lp['alpha'], # lp['dist']), 86)) - # fh.write("\n") + # self.outfile.write("\n") # # Other LP types - fh.write("\n" + self.outfile.write("\n" "==========================================================================================\n") - - return - - def check_executables(self, fh): + def check_executables(self): - fh.write("\n") - fh.write(90 * "=") - fh.write("\n\n") + self.outfile.write("\n") + self.outfile.write(90 * "=") + self.outfile.write("\n\n") dice_path = shutil.which(self.dice.progname) if dice_path != None: - fh.write("Program {} found at {}\n".format(self.dice.progname, dice_path)) + self.outfile.write("Program {} found at {}\n".format(self.dice.progname, dice_path)) self.dice.path = dice_path else: sys.exit("Error: cannot find dice executable") qmprog_path = shutil.which(self.gaussian.qmprog) if qmprog_path != None: - fh.write("Program {} found at {}\n".format(self.gaussian.qmprog, qmprog_path)) + self.outfile.write("Program {} found at {}\n".format(self.gaussian.qmprog, qmprog_path)) self.gaussian.path = qmprog_path else: sys.exit("Error: cannot find {} executable".format(self.gaussian.qmprog)) @@ -583,16 +573,211 @@ class Internal: if self.gaussian.qmprog in ("g03", "g09", "g16"): formchk_path = shutil.which("formchk") if formchk_path != None: - fh.write("Program formchk found at {}\n".format(formchk_path)) + self.outfile.write("Program formchk found at {}\n".format(formchk_path)) else: sys.exit("Error: cannot find formchk executable") + + def calculate_step(self): + + invhessian = linalg.inv(self.system.molecule[0].hessian) + pre_step = -1 * np.matmul(invhessian, self.system.molecule[0].gradient.T).T + maxstep = np.amax(np.absolute(pre_step)) + factor = min(1, self.player.maxstep/maxstep) + step = factor * pre_step + + self.player.outfile.write("\nCalculated step:\n") + pre_step_list = pre_step.tolist() + + self.player.outfile.write("-----------------------------------------------------------------------\n" + "Center Atomic Step (Bohr)\n" + "Number Number X Y Z\n" + "-----------------------------------------------------------------------\n") + for i in range(len(self.system.molecule[0].atom)): + self.player.outfile.write(" {:>5d} {:>3d} {:>14.9f} {:>14.9f} {:>14.9f}\n".format( + i + 1, self.system.molecule[0].atom[i].na, + pre_step_list.pop(0), pre_step_list.pop(0), pre_step_list.pop(0))) + + self.player.outfile.write("-----------------------------------------------------------------------\n") + + self.player.outfile.write("Maximum step is {:>11.6}\n".format(maxstep)) + self.player.outfile.write("Scaling factor = {:>6.4f}\n".format(factor)) + self.player.outfile.write("\nFinal step (Bohr):\n") + step_list = step.tolist() + + self.player.outfile.write("-----------------------------------------------------------------------\n" + "Center Atomic Step (Bohr)\n" + "Number Number X Y Z\n" + "-----------------------------------------------------------------------\n") + for i in range(len(self.system.molecule[0].atom)): + self.player.outfile.write(" {:>5d} {:>3d} {:>14.9f} {:>14.9f} {:>14.9f}\n".format( + i + 1, self.system.molecule[0].atom[i].na, + step_list.pop(0), step_list.pop(0), step_list.pop(0))) + + self.player.outfile.write("-----------------------------------------------------------------------\n") + + step_max = np.amax(np.absolute(step)) + step_rms = np.sqrt(np.mean(np.square(step))) + + self.player.outfile.write(" Max Step = {:>14.9f} RMS Step = {:>14.9f}\n\n".format( + step_max, step_rms)) + + return step + + def read_initial_cicle(self): + + try: + with open(self.infile) as self.outfile: + controlfile = self.outfile.readlines() + except EnvironmentError: + sys.exit("Error: cannot open file {}".format(self.infile)) + + for line in controlfile: + + pass + + + def populate_asec_vdw(self, cycle): + + asec_charges = [] # (rx, ry, rz, chg) + vdw_meanfield = [] # (rx, ry, rz, eps, sig) + + if self.dice.nstep[-1] % self.dice.isave == 0: + nconfigs = round(self.dice.nstep[-1] / self.dice.isave) + else: + nconfigs = int(self.dice.nstep[-1] / self.dice.isave) + + norm_factor = nconfigs * self.player.nprocs + + 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]) + + thickness = [] + picked_mols = [] + + for proc in range(1, player['nprocs'] + 1): ## Run over folders + + path = "step{:02d}".format(cycle) + os.sep + "p{:02d}".format(proc) + file = path + os.sep + dice['outname'] + ".xyz" + if not os.path.isfile(file): + sys.exit("Error: cannot find file {}".format(file)) + try: + with open(file) as xyzfh: + xyzfile = xyzfh.readlines() + except: + sys.exit("Error: cannot open file {}".format(file)) + + for config in range(nconfigs): ## Run over configs in a folder + + if int( xyzfile.pop(0).split()[0] ) != nsites_total: + sys.exit("Error: wrong number of sites in file {}".format(file)) + + box = xyzfile.pop(0).split()[-3:] + box = [ float(box[0]), float(box[1]), float(box[2]) ] + sizes = sizes_of_molecule(molecules[0]) + 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 + + if type == 0: + nmols = dice['nmol'][0] - 1 + else: + nmols = 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 + + new_molecule.append({}) + line = xyzfile.pop(0).split() + + if line[0].title() != atomsymb[molecules[type][site]['na']].strip(): + sys.exit("Error reading file {}".format(file)) + + new_molecule[site]['na'] = molecules[type][site]['na'] + new_molecule[site]['rx'] = float(line[1]) + new_molecule[site]['ry'] = float(line[2]) + new_molecule[site]['rz'] = float(line[3]) + new_molecule[site]['chg'] = molecules[type][site]['chg'] + new_molecule[site]['eps'] = molecules[type][site]['eps'] + new_molecule[site]['sig'] = molecules[type][site]['sig'] + + dist = minimum_distance(molecules[0], new_molecule) + if dist < thickness[-1]: + mol_count += 1 + for atom in new_molecule: + asec_charges.append({}) + vdw_meanfield.append({}) + + asec_charges[-1]['rx'] = atom['rx'] + asec_charges[-1]['ry'] = atom['ry'] + asec_charges[-1]['rz'] = atom['rz'] + asec_charges[-1]['chg'] = atom['chg'] / norm_factor + + if player['vdwforces'] == "yes": + vdw_meanfield[-1]['rx'] = atom['rx'] + vdw_meanfield[-1]['ry'] = atom['ry'] + vdw_meanfield[-1]['rz'] = atom['rz'] + vdw_meanfield[-1]['eps'] = atom['eps'] + vdw_meanfield[-1]['sig'] = atom['sig'] + + #### Read lines with ghosts or lps in molecules of type 0 (reference) + #### and, if dist < thickness, appends to asec + if type == 0: + for ghost in ghost_atoms: + line = xyzfile.pop(0).split() + if line[0] != dice_ghost_label: + sys.exit("Error reading file {}".format(file)) + if dist < thickness[-1]: + asec_charges.append({}) + asec_charges[-1]['rx'] = float(line[1]) + asec_charges[-1]['ry'] = float(line[2]) + asec_charges[-1]['rz'] = float(line[3]) + asec_charges[-1]['chg'] = ghost['chg'] / norm_factor + + for lp in lp_atoms: + line = xyzfile.pop(0).split() + if line[0] != dice_ghost_label: + sys.exit("Error reading file {}".format(file)) + if dist < thickness[-1]: + asec_charges.append({}) + asec_charges[-1]['rx'] = float(line[1]) + asec_charges[-1]['ry'] = float(line[2]) + asec_charges[-1]['rz'] = float(line[3]) + asec_charges[-1]['chg'] = lp['chg'] / norm_factor + + picked_mols.append(mol_count) + + self.player.outfile.write("Done\n") + + string = "In average, {:^7.2f} molecules ".format(sum(picked_mols)/norm_factor) + string += "were selected from each of the {} configurations ".format(len(picked_mols)) + string += "of the production simulations to form the ASEC, comprising a shell with " + string += "minimum thickness of {:>6.2f} Angstrom\n".format(sum(thickness)/norm_factor) + + self.player.outfile.write(textwrap.fill(string, 86)) + self.player.outfile.write("\n") + + otherfh = open("ASEC.dat", "w") + for charge in asec_charges: + otherfh.write("{:>10.5f} {:>10.5f} {:>10.5f} {:>11.8f}\n".format( + charge['rx'], charge['ry'], charge['rz'], charge['chg'])) + otherfh.close() + + return asec_charges + class Player: def __init__(self): self.maxcyc = None - # self.initcyc = 1 # Eliminated self.nprocs = 1 self.switchcyc = 3 self.altsteps = 20000 @@ -604,9 +789,10 @@ class Internal: self.ghosts = "no" self.vdwforces = "no" self.tol_factor = 1.2 - - + self.qmprog = "g16" + self.cyc = 1 + class Dice: def __init__(self): diff --git a/DPpack/__pycache__/Dice.cpython-38.pyc b/DPpack/__pycache__/Dice.cpython-38.pyc new file mode 100644 index 0000000000000000000000000000000000000000..b171a9653d3cb944313dd05c3ebf7330980999d4 GIT binary patch literal 11245 zcmds7%aa^OTCZ1ES6BD+JT#i8G}C%m));Ff%dfQ%^4iE=FzldRX>D(VcebaqW~MdW zRXtfXlBQjaHyFc01P>9h2M!#12rjzq?jJyKfSZHhu)(2%;KG3r2Ou~Q8y&%l-}hDZ zqgo?bm>?iL6a8gYzA7{8`+cv>Z;g%F3jU1fv+w@s4Mq7!dKmuG@bFdq{J#PbicsrH zIZoBG%FkL^<5Vx}NVU4r$dofGwbAQlBU{dL-l$u*sLZJ6lGH|;Dd&Y*E{JS-OjzZj z$O-#Xr93Y3qJVTljEN%BNii-akWPt7F@^M)I3}i%PK)DW2I+AzE9Q{Sh&JV{w^K zu<2uKUE2?w=EqkqT^8+bHp%=G?cDc1U482zHcN@@_#@QZS`cNuy}+$BoUXlkZ@;x# zulT`*DC5!~4{6UI(l>|<;kbV5{5w*5@=9r=B1-OFW5bd3^&B^lHODV)db0F9cYJg7v_WQNI7=>_&Olxw}r2tOfhMcFj~fJfWTL@B5MA+{4P_S2Vo< z3#e7Sh-LK+ZMpJs#$f7M+7d07K^vpSa@CmbVjwc38Xz*puBEtmL@NteU-j$?~edt;m08PL)3&y zpxjo$viP;tRx!{*RcHrVTl+w9PpAt0115&hI~jryG1s>hSNl*2wT>x_PA13_`^Pn0 zW^!3JRc3KnE>&i8S-w{$3LPt!IT&$#UbzmS#E`{%Ch&HB=pB>XcdFNy-+PS9rh8== z<@mgEUzHZu&d|&`H_Pqj0M+?kTX%_T7E*0;VzF1xwI}+uC;Pdjer~y+EA?}yC>Pk` zbdc}l50z=GsBQ3AXRsclzn;P&cA24icCDLTd;iudv9fx%yO3b3T5(-3D0xlCEfIf} zK5n;4K`lV;7mzPqD4|p7*0uw5)oU%5ns$`o7Ttnd^(3H3E$zsRpSp$mZ9;gjG=A!i zCBTo1Wge{3RZA%rBK#9pefrISgvv+7wi-af9;okQA-hmXHOICpOo)@C+|@0|b?!CgwO^vK3=l)nl-6<|!`Oz^t97qhsr%Pf`>i9xZFu#rO*Ecp z*Cv!Flb*mhT|00Z&17)~dC z5cK`7eTPzXx@33lh97v8e+h$Cs@1&)@~(%6T2n%%#Wi-#uiT-HMvq)gicQW+f+6HG zrWRFGozbQtVeqfbsS}iE$+G|& z=SU9WUzZ|by23cB- z=TMs{aF;f1XtR4OaC-h#b-CHy+dK|Hr8O0kQF@rjBjur^2$EZ^8!Y86g>d7yy)&v-~#3 z?&f>?6+j!n^dWJ)d zVtoT=`9M9;t}7pv5P34dcqj?+Lv5bY19e(?9nq=XW@rfW4dg-{Ipw#tdAyY!_QEW= zm2k1)xLb+ij`@k)5sv#~vkCG-C-?fJvuzV=853~CLsP&J2Q4ri13=(z!w+bO}R z>p9Y;np**0chP_l_gw&g>FARn>K+tLveV~yP^eSGEGOP5w(+iaJP z%-WGA@8P$PUrRDJ#xYi}->ZXu>8#x{&N{AGITm|_?8=!npPac-6-8<*%8{R1uepvN zWe~;Na^xb4V@ELuC)ug))g7NbOV(*oerwzFgY`<_HT=k`Hqsdy%gI;b8(&TSN2<0b-#8L;!a``h^%5;#`1AEDDn7m)=3#3n5{Rg1(85 z*DJG73UE^VikLF6U>I`X7diU?kv61xN(tKWIK604+!D}u8^Qr~p^7(W1R`gkHd6v> zh9Q+9NFfiXfC40-6PS&~(EYP`EHMcK7yZQ)-iL5B#=Z_(S z&j84&*QtLV{a7PPLLJ+?L3|GfxaiXCVjR%+R*dW{VU zh{awg_)7r{HtJ=~t&D)6tXH?U$~po`2#wXYE+aA?6 ziIAUMF2sP2J}9|)aREjY6Dq$wIB)>4=7Oi#8!acfI0#jVl4s7o3@ zrJPd^F1_^g4|9NCwF5;-yg!0o9mg>!rna!X`(G(^t-LzpOmB$Ky=8XCo`V; zA~VK&em6a5l>GVVn|anV&ot|Q812imo=wmCzUJ14zSSJQl_59cFLhs&eanFxSz(y;b;sFBZUYi^)n(9KmZRlP|aQg}$6F-yT_sF26H68(zuM zGkiH+{>f-x{&LIdF9+@|C9^*D1!w)!ky#sZdi2da>r>A#YeTM$_T^ch9+@@C4U;4U zYt%Eta`VNT{bTwPNi&M^;6RGvt+h&kvyHp8V7r9FTsqER#hDh{ z>+}k5i%z{(v?${Tf=*D7P(2&^q`6x!(1C1!Ey-(+qsfnRD ziJcEkDw)xTG_)ND;Aujw3;nJa=#N)6Gypu6n!kQ^83#i@|IFa9m&WxjTvv|qC+bQP+m~Nde6>vZ=vi-dc%#qkKGH{GoXRjv8(gg!B zD-35MmCOrD=^u! zf=?9N4hZ2`bjb%H95XaKIam$ZBSSdU&x8a_`o$JwU?Xh8J{@wfawq^3z!KpZz_+71 zTph9k+@nH`9{YPqm?sxnaY+b2j7Fj~3G-OHpTS$Q;agcq^9c7?VHUYK+>?#hClK)A zD@z$#11um7YskV^nt(lJcE&@Sq%(Y;^LSzliY0Y%NK(hK3xn?Sv`(ju-J_x+LopBv zf5h^UBN*(mkQK@@b!Bn{rLptECU2b@og#9Y2s<w4@SCF!Tt@!)mQRE7)!%z z;Lyae2s4gF$bbg8OC%j=h$i_gew{8ct;dSxGlCWAh#stFz}&;(hx2jx0nA>2YhwdY zxNxJ%@1sX5u#jThpQ1JZ>w+1;m@uyqVmKrxUk$X|xI)8`#mfPL1_&ppJ(Ct}2`9fF zS+u~G`jKo&{W61oy)XoYADEdDNko(|VNv8b;DDQIRGP3`8zBU8;e!IP?tuyhdaa#B znT0!O=*bKPW`(JaAOo@F?q7x`m{*CTeZ=A@HFtlA(#)a~6v6U#X!9Mj_&}H$gZvY@ zMdT6@g7(Azr1%-Rb+ONjF95OFsQ5lR%`4+E-%L^ic7&hKL6X|}M3`^-E^#!z*QmVZ zZoRx7c-$N_Qpm6So(%{uWa7&>U#~RYzJOYsju^7h>OY=T;TYT{p4(Ng=~be zFjo<`DGD0*NO^Pc-xP{YvE(CN(mIQaG<21QE)0T! zlAK1lP&?4sw|&uo+}p9`4He9S;B+w{smB3ZVjY3f-_$=G@^$&R9WE?+lYA+RWpI>x z=U&xm(qaDqrYpt(dF4#q+rs4$-H~!+@zK=ae&%|jVQme2K{7%?%#O~oa5DraTk`i& zva%Szg1_xFP2&>~=w#dgmFTCGGYDQFPqtt*7T<(hy>Y9ufxCB;$>U8Y_*&iDz-^!? zpA>z)=2zvXG~RC!`8^_YL@p4aJ1KnnE|)1+B0|wY4miZSLCT1M%^!#7IDp&5pC&RVg~ST9)5S(lM3TAwIIJTLwivx&Q> literal 0 HcmV?d00001 diff --git a/DPpack/__pycache__/Gaussian.cpython-38.pyc b/DPpack/__pycache__/Gaussian.cpython-38.pyc new file mode 100644 index 0000000000000000000000000000000000000000..e482efe8e402e5c0933d37d3b45ee1b5182ce506 GIT binary patch literal 9011 zcmdT}O>i97ah{o-ot>Ti!4Dx2ph$U46BOZ407X-ZjgSmY(4vxvFf7s*g-C%G`v$wq zo!!}eGa!NGtX&};9a5@_Qm%5jQaLP+%rS=?bIKtHmk&8krBbQn#7hpYoKUG!zV4a* z13{Wh>5|==p7+=9XI{Vly8HdOP_PvI^yt^0{qgT9%70R3^rsTE!Yi>Gu8_7F^cS`6#S2j$M76#r~kr5`|mdFYV@0`eqJl=Ux z5JkKTVnR&fT@+Jd8t(~lOw8atDb(A_;%iSaBV}z-jnr?KX`20C0f{`GRXqM9RZ%KR zSL@h;5-N|?zbbZ9VXCL@D_#8~t+UREbyD5*2~48 zAT)2F6^#>S0y8uptG}B$wlTxvGG!1h{B)84MR zj{UGJu}Mm2dbK12$FG%|_Dxw5HK!RoO`)Ppwf!G$w2tRC9)AX;f^GdY-&L`_>SOh< zwYGj+d7?Z~qx81xm4YW)q%B>D)LofIVUcjrCTkTwlIf+VN>}ORc@gy6o2cpzJkTS5MJIXhw#jDV;OY zDT&Q(PnI2jbGuxvFI{%Ma>@0tBGG5*oU-~0w5h$lExktnTdY{62`u=nujgm2;ww7R zrr-AiuTd*=@-TUiCv2a4A3Go&XZgBUm5bKePGc+K9Z3#&?+#K);^V^w8`jseinpxJ z72AHea`he27R6Hr}l@9Y0F>dw!%lyR{%PAUTav5T#mPD@yx;thFM2yVev@ z$|a_T=Q`yb7fflt7dctu0|%wFjrzW5RND#oC+ia4a(k5 z6BD>!a(9q5ae(;9Mrv&%%bjduA(sBCL(=YnNfuCI-nqdi#}hD;X>`@ zW2Af`T>LTw9w%_2w*R&L%b?IJQvcfw8No!D?gIkS7LCgJ^W7Yjb_@d9&7Fe-~J&st#-F=PKaENjX8?O%j+{0G z$y4cL4Yth|Y^DBqg?$A!mBQANNo|5+m?0q#4PijkHRf3AW+D7$Ebfn09FI(xWr0sp z5co7iK0{?VM1M!zu*r<}7hxK9OlDN@Pu;0Hw(nr0g)JOeyI%@w_Z=GwT-J71MueZl z_WvStEMoF)aI=81>_sV4W}{Z3I*dQJXjWbEIx{}DNJA3leGH46*5Gf5b_nQX<9 zfE}Pd=$w0=*?hQh=lBm>L9J2iKnlI4UEV9Z&I*ict8;9qp!;^G6-y_XMIDQ6m$#R` zTj{c{mhYl8QmxL3A>SdzZpme|S?EGKyR|-!dh@W-&xEeLGz08;@+(Z9^?Sa8R4?=QVgy}#=g9W{phuLsa&1;m_ zSreL&+gY`!nVS4Vv>)>x2@^;svWXxd)$j)4V(hDddQXE);pxImNOCFUQjl(FIsu(T zKIOoZB%d)2UP$Nzz*mFU25)U9kuH%3 zcVS+^e2VbsGGc<#|NM$-CuBmwiUgW=^Y?^yMMzR^#ikZ4D{{uCzaULVF?oEr? z#Ala+w54|f8aK}9?br9ydPg`@)7oTvNsi%?Wgy$^IN=Lw2Swsb0$c| zZ8)-Eg;sAGc$pi+OFPkai9GcCq;)3Ig&gKP3EvrOgV#*!*V@~|nT6!XpYIm9f!14p9^>F3FMbVA4cH_V-0a}+T_qHWjJ8W+r6c<*j z^($Nfzr+F-I*kKs;g-|5^00gBx^<9SSZhhI>?3PMMw8Qcmp$qD2dQtZ+fmwe?mO;* zYDZ?>*?Ztg;df>dS8cUgk`>2Z^KS1o8d27NPz%acD6GzLd);fHPy2?`^lQQ1mDL?7 z;aSJlA{n>`Ir}F3ApgpT@avv#(2V~8q_befsuM8)IBv0Oe#>4ZC1^Ke{~A+0FdlZ# z)08viYNfH`25ec&B6WQJ;mYNUe%Hp2oZAaa?+hH>XVO-r$z42qcsd8?uD4&3D^qFI zgm9gM^W$hku4g|mVb_{4Z9qg#I_1EVdj}VO=Y(kdKD=yWeJQLj%2c+z0KteSGmI6H z=8+TZw2nd(IF?D=B{EtZkBM}F*i2*waD5-h5g`S;KbJN9qMF2}j$Q-}7pUQ3Kc>w2Va*Bd2}nu!2z^Z{!sm<92h{de#!#SDJq-+v#xN8}1`H)mMo&=U>|^NXV;D+pN6-`L(YAj- zTn}b{W`LnT9iR=$7;WYlZ5DXFTsDF>$iEVX4)|Yyp|9Art8 zT^xr`Nj9y#LF5NS-Xn682u0oGAAv+B6f`5G{1Ymp=t`$(j{{@MYXyinda<{#Mos3P z1CUmXkUylp7!S`{NnbApM7d5JKO#at3L~NXkWv)*mv@N#n8+<6KOsWOR^BD@F%gDf zR}gYwzM?jl^%hhH$1x!9Nl3REPIpdiw+Nc9K|23Na$YSy<0eD)wzmy8nXKx&~jQ z&zps1jED|-J@Bn9q)a*;u<1?c)hvRJ#Bafq&xB^I{n}pyxD>?M%hxq%uzNwC`^>&n zpV{a1`R_piy4~unZHueq2cVU6? zd5mCai8w^IK^AjEGbs*7nU%44RGnISARBCu){`kjN006XW^hXY^JS15reQL{I7T01 zp()alJ_wHS#*99DhoQ zROh0#I^t+}+lc@8Fx+d9u)q;E9Ga@#zsR%^C1N+J-GHA|iG%rP2esro!J&bg)K+Ad!5rqj010>=2tOH+%e+e98v$2rub(dkx_Ai@&kb&P!&#uB*BWZh;uoc z2wXXll6_R&@X0UD_sxvFo^{!CKEj}VVBo8(IX7|CK=XIjwTB! zD+N2x+&ksvPII&GRP2#tokS<3@J*P{Fop5fwHc^+%FmG0{2cXTN|%IH$HP)X=yQbB zAVkDcL(-$uwLM7<_L$`-Freym5laI?mne@F=Rm@E#}RJ6hTa*{JE!4EBDTz0B|zjE zXSiIC)e9|pX=%xhGLVV}T;opcJyR-u!)rR$qIqaDB*F*$rl0 zKc}k4L=GL&9DT}U#9z>toJ6Dkm$WoE5%l#wN7>?(WHI2q+@h{%yr>Z0tuL+LDQ#if zrjqP!C%EZ)Tevolrjw#up8I~ONjIjIDAzBm`DMQ2<68|5LvuWjoe@b$$>*K7k7qG5 ve_jb;XJlS(cmj8-SFr%ar;D6PNPtNW75#N{$vkbIGQVMdVt#C%H&6U8WY_}7 literal 0 HcmV?d00001 diff --git a/DPpack/__pycache__/Misc.cpython-38.pyc b/DPpack/__pycache__/Misc.cpython-38.pyc new file mode 100644 index 0000000000000000000000000000000000000000..b5da37251314c266cf088ca7b17ad0da634db29f GIT binary patch literal 1627 zcmb7EOK;mo5Z+xrL`k-y294tuK|-`qQ58be2YNArphaCEml$r-hb^E`DehXLC5rU! z8V(idl-~0fsAKN^8T|!&D$rB^f=``U#gQDehc3Z+?9R;2d^6;a+uJn)%Qb)e`ln0C zU#M&@2Fh3Psvm%G!s&$cR-E>z)me{O?(`hX-JZ)CcYY&1kGtFh?sK14fLD0HtH1+Z z<8|OwPM?zY*544B9K=XBmjlW-@G4ADlag|}M0<8H+iiW_30pjD4Z_xO80oO}c(3(j zul4j5Gpvat4jmClmfY0m<1YG;C=@Ifm@OV1i=_E_avg<;{h!>*dzmoBK;f z$d{( zF__)`R?1B7g>f{UW;%?sG?zlD@Fbat@cE+MG%Qob9||2`aO1@zo=D@vsTR_>d89|i zPZV04$|O?brc?G;wG84wlaxpdTaF*o8SEF;COiVq_@&f&`!Kranp^Ph@B^qL9 zGL3Y*B5wfDrhbvhF$_{41A1r=hRAxDgc#0x?nytH_PmpRcCI&?Z2Nojs-C{}?0=RH znh;gp0YX?nYix%GvuNPD_HVhasS+taRo8 zwERfBrNxV#L!sy>_v_ZcpKhrwC!J=w!mE>yoA*N%_d;l06ce?#FJ}iIoai=r+H6kT*1bkS9J1&S=*rUkkvwCc9)BJTH{`}jB< zSy`lg#5?!g-*@l1=bpdIZ!RvD3_M@#exvpJ6~p*e+3X2f!Q}2IyD=n zSpQ(VYu5K z?H&U6`rTox-vMS`ZY+M|_MPV4dvD$R`rY>$?{xlr?WYUwfBxIoo|x#^X77_vPK7O$ zh6Y!TRhxnJp@CJ~frHcu+`vQX28F;!>IKE1gtQPW1ZAXtuo_f?r4P+UF<1^(P*MsO zG1s-!yL%YL;h@b=PaonI-kW%%Wq^*67)=fF^}EI==2qQG{lTaow)gvC>JEWMpO~o= zh5hYM3{^&_|G9Z}XEX?}?sS9ja1`xyuijCkj%p1C-C^hIox^x%G`xDN+YayaTZf^# zdh5<^tNq~W?NR^T)-V9MI#+fNQ@`2l4!d!)d6tuL01V4B$4h$JSEjXXQfc}S=PUDH z2a*^)@tSsjwc7U41A>1$~F7l{D9uRpC zd0X;;$QNSk$WGkAUNaJp@T`#(bb4$~{KUtGIa`XmsyK%y+}paN*Ij9c!#GsU(ROpt zilWpRwBjA?PpmDqhWXS>n!MJqm6kg@qbP0;_Xqbwm0Ak4Q=pw1t0lBlX9&opsg0A75xwDc zhe3F7E!F8}TAY$EO5IkBt+wk2`&)u<@D4UqyXX1u!fJ)Gx_f zHu<|`jn8I#aJttV@yYi5j1d1hd&bK%KAt$SnbpzMY zj>?IDY#uG*q*%w+#FvudiIO5pDoL@ogm%_sL2?VoEg)AOug1$qD@pm-{EjtQoK(JT z4A)J}6KmnBWEQ_;46WZbk_uFb@dGPxwsBs2OOvI%#LY^UG4^V*lq_T96`99U;Ne6s z1+G8`a9K*0rDQo^$%>S$NXZIII+*=rHCa7c3@lmGs-APQ)LS~X@a45+4YLM3HlH^p zYgz4bR$EP~Qd^bUs?2FEsbWrNWKL`Hy*2sXn#`aoB~>Y@O34{1IU^-!@;R+1>(e=5 z7VGle^<)(zpB4OU;%9U418=esR}=TKiT(TTcY8Zj5Xv<|6c+?!JZl^m;z3?EWqO75r`aHo4>|hOdA(5179i{~pMmw$DF!eBaRe@=8zqafn2xtnH)YE(eNUf7Od zCp1CA#-a*^ZtalS)DF6X#_BM9)I8O0g=ORv(f!d>FQ) z)M{PB3$@pJIW^x(t^3H{-)va-k$pQYsZk8wNy4X&8bm>AcDaHk`o<&GSj;962AtSw z0qwI#6D;v2J|BG#z_3bY$*Mv$mn`4(<#kL4>4voMzz!E}AqeFYdjR7)K=tE%VrdICb@iZ_p@m&|BMS*3p6%s|hg& z;@j2wmbP2n)dtOzPw|XzlONB47{@}-t?p^_%HJJDU7q38+TLj_O>;1NoA{&*soq2P)10%kzZ>9$* z+GLs?+ifjl)63ZUb+ALp9It1Z{Zyx|HjUuvL#vF|uxR$50)g3MJn&RGfQF_C1TFp| zQkqpS#$XD)L4iXXlrG3)3rz4M=+;>-nG{dt`qbcwU8QVy085EmH9GcW0SqCGm4KCw zTPe`evDc#F#{JGq&zGmSD5CKa#DC;s5?m+bMS?cphICvEbiUedwcH{VW}c| z()2GRi}w;0QHJ&zq$^5D&7VStABv^ z5bwBAY?ED7cYo0EVUURlr+x#?Q?r+v57gI*IV$&Av!485DVM_-sAx?S_Be3`EUh&uwTko|@0F3PqCM8j56o6@ZdOK`_h$Hah^XDEx}BGACp#h2`&@Jo{1{(Me808?yHNes1fA) zuI}o?>$dq}^-WevFE2umr^S3}S9OPRvlEU6e}}fwDu7`A*gw1AJYU02LaJqq)o6WubR$tdB9$vLVFyl6g9rZ1O-z8`gl;;6mT1BZ( z9W35HRF3QW{(@ie-|<&U7p1n1+PQ~}d;@RvA;2#jmRJlb2u2ixrC?ct5-Y(f(uH6x zs3I)~XM%O4i^17o18F5V7o11B6kG_NL%JM1A6!Jb5_}FSp_cl$vxvmA!x>jFTOG~_ z7}P*GLx%$Z5kdllY6ElyhxEg^aOvnf*9%$}-vi$8!XXSRgeMTt2)eM=(!x#{Mcvjg zEybg_)faDCeH(2RgIp)5E(GUh{R7_25%!?Tr@!&abQQDc28lg=h|`=BAx)WrgU0BF z7MS287_er~{t7}kO&r1_hwjuY-r?$1i0iqH-t_pRNbIhjV+3355hz#ud1Z7LGV(jWW89Hrnyj0+f z(FMAJ^;@!s6pn~yyr%b&WzJNjN*RPC%4j=DI7rIeJvg6;Y*DQaRHslFTw$T1|-P#eriBS*nnKH ziLe390LEe%gK>Hucbx$+;02lr*ne^ga^#e~4b+&gqb&B15P@NUOdLg<4khKP_( zag61W&m^2CZFFtobzJa}x7GRYw?De^!K-UG)3$m2z5eUIPgO^}f>~v=S z41@-W+*kY1k}R0%{UglX#^B*0i0lcWJL0 zmyQ-F(qQ)EKpex9R?lN%^+m0wh5VwIKZnxO=XpM^rgu&DYLi3wn(u#^9luKOBEhc^ z{QO|r**pf{=P`a?DWs%07yM86PY=&@^&}}Svz979WmVwo10kJB1l!)AP^rxBlBZ(SUd4u z=VMd-MLr?|(ElG_&c_o!Qsz%g&ybj&A@qUF5R4_MhkkG%5)W=cA;Sy6y~NclCM!U_ zxtP7Xuto|M7ikun6;QC5($f-V#J1d2g=z$2vw z((-}^g9Kda3nzreWtKima!8CseT7Y41<2SGF^KNLC~6d^e&Mwlx;aB%POXQCN<2hV z;^E7u1^EcW6)>adB_IqhRIE$pCDTVBVOcQZzQlkTk*}B=3>W;&bH)>^<8?a9xKK>EG-9nFg4#7&fxbTsu(`S*c;JhO8)Nu%% zO1=dBUI2d;MW&auJl8e=IDj=$hG@J|>@aGiO!qvC_yx2DCBM|r`K;oePfDQ; z^3qaXihBVmMZ08H>h4{PJ-VLMA>a|{p~OYX*3eYs1<@j^3$kbpLRGIh;x;hssqB7@ zrSsZYiE$$Og^~(vL;5!C_YMvS4hhBphr|wvjmLVhpD(j2T^Ooyx79?@E*hyXkPRp| z4Qsd?}oj38L|Az=snXK@$6SmBJ|2YnzEt6|@^6YpW=LtJ}WE~^KipUO>$EiTM zSgBc{9CQ?S+}awy{_L=3Hi4R6z*DhJcw(G~`Z9qcAeW0PKDBVZ#)2CJr!*>Jh}~hk zy9;x@br43^NgQqLSRe+2YbOXX-*%1=G;(u7ynzh$XU6vDY0!cw7R+)B>~8fOG77GZ zO0mnR28-pKAYdq)vbJ4l3`aUu^6FBWS%dlNz^N6PlF}@49@PS=N z^T0l~s*ns>E!$86xbRpUTv_0SU60#0me1EojzmQe=-2|WVEKyui!-=^^9^qYP8-f2 z+%=v$n7;uM?FHDuC{`~s`wBs*>v_U()ftDojdFWmHK)ED{Lb$7EXONR+Mp2sbvo!z>-oxU=CNz*9buq&sZ5hdPIZ zxv|b6!Nr3cb;n_JdIwgdRIaZx-~Zd(1kQ=v%)mLVF&zt~mdk$Dx-~wZbNEw(%Ihmgu z!G*Csrq&5$PjQBP@PW}Qf+OC=&eBYFdx(wqaIx8|h&|$=bP2l2>)|33?V)#r0xTD| ziD0Q*)w#_@V7`=9_I@bmrnn^m%}*??xz?~w;JE?=&Z3z+H2y{=v)iNnzP@i7v>t@D zRt*;-I`pn@MQV?Bx(Lp-g3&H-YijXMSmR#?f|d$uxK%>*M4HdOTVEFEQ;YKmU&TGq zXdjnOl=QgT3FDcN*C$%kYAzU6jWhi`!3zX015E8HPwNpiQws00i41TLSfebLMwuKv z&2H7{#CaVyo%nFZapbNtu4MGT8decr^}v1W$de1l&(An6Bf!aV+9$(26bB@}h=eU5 z^v7^WsYU3GU8lB?Kmpk7oECr(w$aG5jlazVT=o>&&J7MPZvF>K0bR%1R(3fdJpxEoCrtGfD;tY z87JyXQon$nPEZ!OI9Feh`X#AZ7Pykvw zdcq(;&njxf>&WU@c1Fs$o~c)|F8En|+o8OdQ@pV?{%Q@+O)zJucrDEEUhU1=r7Qye z%J!u8(lc;GI6!^vxX3{tU*lM}@i(qv9C5IzG4`ee)L(x`2i8xdrSV(RrPafYcWyG( zS_j>elVC)CI*HSGy@WWVOoITph@cD?5ky>n<~>>D|M4>%fBKmWm(LK-R@230^)I>_ zm0%n-A^3fcc!Tu0iGWP2&D)(Q?zU^Vi0{NZIfY)mvN-*Oh1S_C`>KqmTS!E!R1(+s ziG%e)&c4sJ-X#zpM~w1_*$)V|2?hXJY)mdscmpPW$_cqVMw>f@qxd-TWPiRcbG7`W5&4gXUQA+R!Twmif8)U$o2#8CPahT2)_Dq2*hdGt8t0L;F) zRbN%#!}nB+;64F^`f7)uO2GK3EL07N(F>nxG_}v{0fA73#MoXk6LOZS0zP_4TxP5quL-2kr^`f2q zxZ6+t7=NRAq*}Wv{#UtssQ+~1;vYZ2#f-7pdO`mu_R8&U)J_+(%$<9!`}ohT9S$F= zEe;$L&_0!2m+zifbw0t=e?#9K{u(A}L@xl~sK{?wBF2mw`<4f**>@^`U|hyK_+Lle B+ZO-; literal 0 HcmV?d00001 diff --git a/DPpack/__pycache__/PTable.cpython-38.pyc b/DPpack/__pycache__/PTable.cpython-38.pyc new file mode 100644 index 0000000000000000000000000000000000000000..f0daee8b92264daf814094d614d34445fd544a78 GIT binary patch literal 2017 zcmaiyYiv_x7{|MJ9UB8=&WYTa2)ZHAiGWke%I?xwbQ+Nu5eeY!f6mTH7C*>o-qYW6`#;a~ zw(YZK6}$Mj$3O4CQ0Q`9H;t)JoQCAd-sB4|-Gw{|bMP4EVjkvW0TyBr7Gnv@@i>;^2`s~MJc*}Jfl92v(^!dBSPc)V@C?@A zS**o6tjBZMfQ{IM&Desis74JWyzs$~TGXK)4QNCYo<}oU(24-spkW)k1`$FS?dZUE zbfOD7cAy(Q=tUp;5kV9K01OO58Lq*&o`;w&jC><+UmoIfUAh}PL$1nQy6D`ki>j)e zXe=KWjfO5-;=1q|IpHhM2|rJOr;Vq+To(bx0})-cDqVOJx(LN{LSrn%n4h*clM~G* zGbc2li!kFIrk!C;N*DF4_ou1B9{v$s==3!(w}F0dmNoqLF}Il~7^O|UF4pU`1GHON z8>ViMwL$Kqfpwh4k8#mNt&n+`7v}s;tm`!UQP;;FEj)hOnyJM)Ey-M-I%+l04zsRy zR2PlxQOBGxV_o#MFy1&wJCYOIsTrh38#SAlqYbma8RrZdbA8-PW4)LDpt&z|g6!)x z@4}py8tv2w43SyqWe-1(kH3STJ?f}iPrH@#1gO)JsfAtbgju%weIr%H4%En&%WzRnz z2i8b+vf)~;Z1q~HwibW=&d^`$ZBJl>R7IxH{Y_LJK0MbwXNy!{wMI)i?^a87dh!0B zw-rgLY7ZnP?w|2GA-`0d#(&@-!3s^UmFni6nsvRoda1tuV(E8xZ#GEv@=i+bZgR#T z{PrcKcQ;GbedyfniW>w|Ch=mcRBXBNR43uvOLs!oFSJQD&)xCRDqWMx+^;URO{z+( zEGSj)*`}G|BZ9TgCn2f&t*!*^)i+Y5x583gJYM(Vt)6zN8oBz(34%51#tx}QxM;^E z!lY-rJ;Ct4N>yfFS#6h8=7nuNOL)n29no#UI|M5iri9OH84^ozQy&QvE@_bC(I$IvQ-xQG)fG{^)oX8TLHK3Fe~3R|ro&ZAqpvk}#?TC9~h~<<1vExI`Y$)X#(LsksjSx5n{L{Wf#{KkvcO M2W9;3p6#CfFJDqv#{d8T literal 0 HcmV?d00001 diff --git a/DPpack/__pycache__/SetGlobals.cpython-38.pyc b/DPpack/__pycache__/SetGlobals.cpython-38.pyc new file mode 100644 index 0000000000000000000000000000000000000000..00317eadcc3d357938641670a339e1df658a3a67 GIT binary patch literal 19799 zcmdsfdyE`sR&Tw!ySk>Q=hgG{*seHF)We#`aXI0ck>5@iYtCKlbWE58R80bhUq6arckC54h3_U2#_z5J*6=5+EXkK|*Lj z3(5lFcg|P+oUzvnv5@Fief7Q0`Ofz~-#O>2&khe;3V!=mUM+t6rlR~Q1$uuv0vGUj zpN8Woj#^dn@~!4o4r_Uh-+ErhTdNwcQ@U17#P8(m$_t94JH~s8V{B{rl%rlz7819C zPgz<}BlT#ByuH6pq3H8N5wr z*vaCZc1E00ye((U8OM9bnQ(G=XF$}XGX zBlD%2?}}QnS|a+|Kk~O9@5=~;>gSb$5^8W&U)xsj)<3Vbj2-eNeneSTxhAHc{JgT7 z*j6to+XRt99+UH=IggIo$~@^#9t(MfI8TQ27@a(sPM%>W0c=?Z*b5^Luogxg6YiKp z<&U#F!EVk`P_Oh-Np05!Zw*fp&jcRtD4enqDpwCR)lhArUw2f`@KG(0u7~QjwvBM8 z{y=vh*ySyMrCvLBsZw$;GlpZA zE;ov$wPROY|AlJ(TCwUaHnt+`Y_;O~rTTi~X|9CHHPRGI8C8OZ+kLoLC{$_{zff2_ zTdjje?`e42aGz9*2m@s=7=jv9q{ty+BEvxy1kF1!OwI^7qvVXinZo1s5`CYkG0VkG z&#M$`ccf@`kfL(&oHBF9Ey&E2|D7mgPh)ZJp7G7M9*byLh2+H_MlM`P`3 zbaLc}{d%<^*1f`VU6foITrb{?gPC?JNNVIcM2T@oYcGysSZOsy=RE7S&-ZC7um7qL3>W-=^Zz(nH zO)OBClv|0bA2qd5Yb70Rht_ed&_ArtDsO8|y=jDcE49PRQBx1G+^!|X_0XWWrcpED zOOr3uf!SJ1-qOT7{!nOat4?A~rx0M0Ynu3t5R3FI<8e~J@f&^y5W|4$#Bb}Wax1xU z8W1KR{-K-Vqyg~{cZHy`RdKQhW=O`|6)MAdmUlrd?+sP^7HAu0n7?weXiciMvZFN< z!OMPDV%b)k$uQ}RtR+CtFLm%FsdlwPfFE7ckn(?#Cn0_9JWhi0pzmz--N>4STz}E+ zy)e=3wOS2$#sHlZr=9UWI!^hceRNE$Dd>fx)c+DgE!!?HWSV-)(?ovn@;1(JO~+7k zL;QMA{~`oRa-c*LKg)Q>sYiQyauP6__?2!Nje?%^DWvNn)8CqLe@yIvcJIS3?R&Jb zA->kZ&Gm^#&@|)y0cYl}tuu?((Z!4165`gOS`4<1E^hAKI&q0}Xaz$|Is1A`^mER< zY_Dytner#Y6m8B|mz$V5p^5oGb=1Y*Qq3sb89fyH$31k7`%@Teple%eO@~Pu3$^*s z`H|glyua}x;FyCz(?lY9Qt%c^<-kF^c zU%SItN^HExJhA-XL{DON4Py{y`|dst+yLznDJX! z3g(^@qAt$ZrDCmC_if=89h)|4``zY8>bcv*3od9;!ms$)`^@&9EQnu>k{jy{FlS)* z^z!k?$axCRiPQPydc9gIdXZjlV3UmWE!T^*T0PQy56JXt!;6yTmAXek@hPAZkIHg} zsQeX|d66i6-MPWrwucgu*qBROr6^fLktJ-%#i~#Hv=^m`%Pg%RF*qAn#=Y5aOTO#a zMY~b=Dt_guC70HLjN!Chrf#yn7yyA5i zxeahjjVRGTBv>KGt$D$j-Bfj|DapHICCmTgfBoJc{89La=OUwqP70nKM0LXU>-MQ# z&F0;?ez%FkjS@9(uwD30+Gy{`zAJoml#29qcm1{lpLm#psnSY$eY5HZrv?c=+3wl9 zkVqrA>RxxN;Nw=x>_b$lyk2n}K9*do`+j{r(i`;#bx>)=tv1S$w!W6vYuD-vbIjsJ zX&Rl?N)6pJ`dsaLMbvBSZp~+kA{`8EWVkmgew0M>uNVC&-M|?Lqpx0z5}q$A=*DXO zhAR+iR4cxiOZqIAu@QIKJy zD9T=>RwPT?Ii5tr)CL{?noft_H{d9mre^V<(JXC7&8VhoX{J60KR=0iyyqzuc{HS3 z>ZqDi$Kj6m#AP&&0i4D0Q=FDXZi3jIGZ&|4fR%D&v>dyc4jjb~r0M-wN0FBvqMpjp z)Ok=vc?f$BH3+6UJb7Q``xtW?vy6S)}UXTmH!EXlsn~XCVam7U#A&H}Z7qlv#y@0fzt;>}MBZd`G>ddPntks;8!A7iluXAWphH%87nhtkhrikip9SBc>*d=RqnvWr>_&P$$ zD%xpQak6U$^8Nzj7$F?~Oh`S1+}Nn`E=>EgG4$v@uw>3?j%cI!Ie#B~$|~eAyy0JH zYB4n8LJbTYVPOARg|Wo(PA{f7{g%4%4S&X=mQ+^fC1hdV-|vhM#1J*e;fyggV@!=w znCGo~V(GTMgyY}wcs)EKHY?tnkg#Su(I@d!B!aDroP zY`pUKW81TdW@3Pt5wfW*m+PBYR$nHh#IXdI;yld9*hOp~9&L5idd;;P#fm_vhDsH$ z7b`6MvPl+ZAHfb|pR?r#b7VoBCR8uRsM;)z{j_~tLStL^Dv+E>;1l*|Ua`@bZV7@2 z`w?na;5uQ~PTIA{?Aob?+kcK_NIVCemEVG_Y2g6x7UD@F`Z776A?G!6SZ47Q`OcGb zft+W_Aqj&xOb+cH;stVEf)i~D=C1bWGeM=c{hO=-$54dc{qw@Xc=q?1CO4K zW3atwS!^&l?I^Y)Yyl(L#JFndf%L?7EQE>Sd(-)t?>T5~vO6>C(yI(n}Yg zyZrJqZ$0;#eccrvSX=wV;_?3Mw%v!UA9Y{Szcyf5!Q5N*O<^}!FVG>W;^Am#MXD7{ zU6$GdoAnXcazHSWg6SUq-M(JL7KBySJ7WbV&V=VUOKORBP}dK>ZJoPsKWnh*I^m12 zF1`NRt1ownv#+^Of!Ld#+|K$_taMIf@cDLrhi^ zZzzu|pHW^@-Xvk{qfoOLP?23$m8Pyjf#xS#N$@x-DPbOm8my`Rfc|D%QKq$Etgxe! zUfu!k1(g>x?rSM7V!%xFXsE`{H_)(DAoF*5@CY0_TE<*3d*xyD2Q`m1JFbs zf0*O5+bTG=KjXZKc3z~6NDjeCHj_~Iu-(8d{pw6J1*KBz2?YwM!E(*eq;heRKFQcj zU`xhjGY5&XO_XgSZmcyPn*IcMR1Li3D8)lbs9pU^6EM5+H!V;LZWBs6;1^$D{MI1; z(4tn0X^^!a8q5P4I2Y&`f-lv{FIyn9XMW^D{!Ezhr5b8JM2kZjpow{=WOxuN%ru9C zQK(`pC&{o`jDE?7e&1bBcaK@kCR&W9^yi4%!mwxex zf53StocWK+f5IBVDz!(L1Almlr5eyz5DWB?`l*h5qet0rUv8@!g2&?%e;1arN(BW#{1w+eS#x%Qy%B@eHrKB~`2d6+X@I}gaJ?Q8 zioFi&1shW6(h4ale9Al^iRsFO#Dl1lltr}VNq3-(D*_g|TbJ7;b`GA{RUuT`ErNt6 zQE9j8Fm1b`gQs=@_O%4V+9w{BP0?5Iy@8(G6SOT6VxbQDt1HBzCftL})4RyDrKW`& z5%sel+p8c!TT}3zfluuT{LH-pKfV|6Z4yS`L&>N21cqVBy#PP4Com*W_X7OnUcllI zfFonQ=y|~_yW4y~%~f0`9o!MheeQ@H{L~ds?S&a?$7@w+A&ZeFZsOU(6GS@nL13q# zM$|E11L)VP_QJu)s6&a;BRqqoRS~a}(-ot!;Olb~$4uKn^1V&YTja2?>jL?5AJpKxDT8r;PHt07YQKX8AD+K>pj%&|$iG*_Gjb})4m6Qf- zP-#gcI8sEprG;6^%o!xD%R}wRQae_v5w*}Y@~R<dw!dMp2uG!%@&oY8u)8U@?-5|%Q&h{Ar6|Nt~jy3 zCVgxx%}hFc(c0jhL>=3zA`&AK8;UKC#fK2S!NzhxsA%-$NW?jM6iQNT=o6C%aUluw zD8G|}r5DFw^ZQYV_cX@@!$@&dQX0ZZu)w6=2sO!M#wk`$X+!KsF5q|)j4;E3#nwV& z)xrbY5(pd(Sa%t#uxnWz!2_$L)iFHdJ86R6_-`>@77T9G5SoCpXPA{eq`H|{os=PQ zD9#DRQl~|@b$B}-SqiSad*=dNCr_jhnA}oP^q2QjK*zNhaGr74 z=Ac-G)kVCL_6)82Dlc4^l6J5G)t6;c3yv60y(&PfN3b~Q#7>lAsxoD-8&M17aE@q)(r)7Sh_87vXVE@Qq9na`fuvLH` zEBYK_U`wDKPz}f?@rNiQ3p<9)$zWzcO4&i!x73qRlXxj!>c0t3yr}V_9Q*s0d|K>C zLya_p>;y=L!TOuFrpJJ+g404ZUZ_cjZ4^>9%g5>7k$NPF0a|`wfCh&#{9Ig~#DKKU zQcv4}1q?aHkmM2=HePzkT7Zmq|B7|zIAO6#Flm)OjB%!=1cma)bnwQ(OZCtaPlh0W z--PJr0u7pg>9(jJOH$0&`@D`FiuuUc#0{!voCsk@ZWE8R>67F?u=mIVJlz~D4b7)$ zn3Y!6g86n~aYTQK5AT7kow18@A=8~K5VNgp`jsl=fGd^Kiex>{V8F;(v`)S*>&cqg z!-&Jhf-GLR;X-!tq@2wUNuxX&$i$^`XmBcwNtI`pAXwvJGCTx?3K%8)!SoPYV_^!k zeL18M1C-N(xy}^Wg{W^5ko-S`hbO_GVZQm`IWTavRh@)G#NRZ*`HWEsR?iadm353PH95jhx#!V-+n7h)-z$ zwJze^fCC}@=vY)Kv3F1B3QD&hdA7bH7OW*HPN94~-VRs~iy`BB7!4Vb|D?8fD>%oD z?Bm4H;@9SFWiiCzo?3hg46EdJcloX0#ZRg!gV|S#H!JI#>)cN^X;e}AiN(`klIO1A zwiM|imy0DI@-FC0&pdsaIB6@$Scg6vHs zu!BNH`xoSV;`<84KCMNG`gHST*z~@ zR=*513XC2Puy%zGsdx;qEjSY6pse8kPo|MYjvUMzcc)oMBV$LHJu2GANxkO)+h1DW zg07>`jy(*CJc7eg?*gq07S;ms6W~Cw@Zw5RkLHCgCz0mmaLyhWnSz}UMx(1q+WmZLOmj8ukVzs8hkPh#rWL%0 z=nb_I?s&lz5~(5OPx%?gKq|*#7u0aFx1kJ*W9ViIQitUYm7D^JL29D~>eF~EhE1?( zr7?{?7|J38H_C=!wv+{Ps*>^NG$aoQfel6Jq^#jMln2a6OCxp|ain^JStry~$Nh;` z&anW4Fzt(h4Q0VFtMjs9wl#%81tm<+pqh^d6TuC2;iU2gswLb>W6IuXj-mxepsoN9 z!O974qcFv0r@zeQX1efwIbk0A!^qtngK6hj+_vO1xRuG5C0|)bjyZo==4g(G>DE4f zUP6T96gvYm)P0JCL=50NF~EQrxH6`Un#ip7Gd`$ih<@gR{<=96PPg_|V5$cG zsi(b~v*9e&p|yWU@ehzn2Ncg>G(k0n*n>M{^s9+)a!ZbOTe3fAI4xn=Omr9V z$8lT6apQ3utHt8D9OAmAL-`fQPxi%6b>gS`;-UD8^G~0JdUK=Y+p@<(m;pEE!nrU_ zvyyK}-3%uj1EYR|l}7)wxexqgyJhr_Z=VB29}Zc7a!`kGpEKSGgMNgmrDQs2N0?|% zU>wLeqC4SGOU|;o?5oH@7@0fn!6?IBj4~`SQf*}zwFNxYQiiFP`o>o~BjkPG!Yuz% zbAGSB^&dbV8saO1ec0D3?{`}P^P$~`o%(Z}>_ZvX--q3HnBUzFxDEEcKg9IyzaxDK z@wvOumk@_PA^P_3L0@9yFviULcrJ%lYhO4)vjn{J^^I@vDtxE~H6z90K48{=h!o`1 z)8L;Ega@2?#HA>1N5!b3nEmnl0Ouz8E9d4|TGfHif4FrRIS+Q^vK;H!REj@Is!fdB zgCVWmk2dY_ASiecv(M%+Wb;_Gc?8;2pB;nXaP#hM*)tul~@I_$6#pQnGh} zd}Ns{9)lC$F2CD&_ssE=PBTcw9zmwVfw6dVD;T^&#^>wUa$5X%%623^4&Q*GO4*ib z1kxa~uAxU~dHWyfr0^gWOw3m>$#8rpjk!&EEv-lQK)dRm3WyVw@F~KPxCDc2D>(R4 z4W`(xD3;wv!L7dY#NtV?(QSdHX9p*6FUN(sAB?N%aKiFew+Q7a=p=P4gaNofvs@`} z3YL=sAcT9QuR4E@ZN=Po`91h_9p203NBXnVM0RF!t! z!4bMnPP{gZ`bulKUj*avWyw?%na|Y9;8E-A)@@cd396S#C>LQ86sb3Y|rT zyx%WFCl;Rso4)XflrTI??ZhV84f3((iVeW2k}f;O92@g##9s@p6=}WlV_qe2y+VDt(8ebUl!F;~5==dSv?UP5*`{HX_oCLE=)X);)LO z0cPkC+5+E~&lEkkR4B3V*vk*&CSrlue^y6Bp(TfcmmjL3hYOO?&!?sGOt{PWjPw@# zx?ilaL=w8McHicQ8t7%$p*oUDJ?=BgR{9dvD1j-pyb`78#cDfQR-|96-^?etePA<> zuNstKS&vf*s}d7*7r+xcgdjtTyLR*C4Vq-k?$S3WFpOHr6D7EPBCXmG+XR!UV~mUC zl{~J#k-U76H6FvY1wOVHxLwd3KSZ4I?*f+OriQ%d1zD+0#|ia-jML`eC&Mp|baK#{ zK@$h*=_sVEhj2tWs*aOsm}(%NZg-Ixm;u2O?qX?&k!uEr7W&WNKaZmgiJ(aPW}yuF zlYzhV!EcaIit>!&!M!aBNgBX`I8cakVEm8#Q=lY^@(6Ddf?E7P#=W&((bletm#EIY zKVtVU;PG;Bc7Le`i>ncgbj}O%a~FNbDF{ILn=7=Md2hmTa3&}~D+m{IF*;C@w5x!x zFu=tp@z2N=1%r5aqU{)OQY`xw#R!Kox)2CqvAL57QaQ# zZ<9kaTC9+>LC$;R{5(0oK+b<5=U2%23OQdVhvu80E);~aXB2fi9L2KO$s=oJQ!@{m zKd2lvm&{kq74tLZyt!x&BPN~AWoHHDg=mP_5XpD^BKiIgoM>ErGHUUadi6zo%7+A* o(Qq7h`E}9;MMLc@6|cmo_M#^;&&oySY3f5S1BYx;Kp4CK45Ivg`k0+w565YetdjpUS>&r Yyk0@&Ee@O9{FKt1R6CH#pMjVG00ZYGc>n+a literal 0 HcmV?d00001 diff --git a/control.in b/control.in index 98a38a1..36ed7f6 100644 --- a/control.in +++ b/control.in @@ -1,5 +1,5 @@ # diceplayer -initcyc = 1 +initcyc = 1 maxcyc = 3 opt = NO nprocs = 2 @@ -10,12 +10,12 @@ altsteps = 20000 # dice ncores = 1 -nmol = 1 100 -dens = 0.75 +nmol = 1 100 +dens = 0.75 nstep = 40000 60000 50000 isave = 1000 ljname = phb.pot outname = phb # Gaussian -level = MP2/aug-cc-pVTZ +level = MP2/aug-cc-pVTZ \ No newline at end of file diff --git a/diceplayer.py b/diceplayer.py index fb79f77..821fdc9 100644 --- a/diceplayer.py +++ b/diceplayer.py @@ -18,6 +18,7 @@ if __name__ == '__main__': #### and set the usage and help messages #### parser = argparse.ArgumentParser(prog='Diceplayer') + parser.add_argument('--continue', dest='opt_continue' , default=False, action='store_true') parser.add_argument('--version', action='version', version='%(prog)s 1.0') parser.add_argument('-i', dest='infile', default='control.in', metavar='INFILE', help='input file of diceplayer [default = control.in]') @@ -27,270 +28,300 @@ if __name__ == '__main__': args = parser.parse_args() -#### Read and check the keywords in INFILE - - read_keywords(args.infile) - check_keywords(args.infile) - #### Open OUTFILE for writing and print keywords and initial info try: - if player['initcyc'] > 1 and os.path.exists(args.outfile): - oldname = args.outfile + ".old" - os.replace(args.outfile, oldname) - logfh = open(args.outfile, 'w', 1) + + if args.opt_continue and os.path.exists(args.outfile): + + outfile = open(args.outfile,'r') + run_file = outfile.readlines() + control_sequence = ' Step # ' + + for line in run_file: + if control_sequence in line: + cyc = int(line[-2]) + 1 + + outfile.close() + os.rename(os.path.abspath(args.outfile),os.path.abspath(args.outfile)+".backup") + outfile = open(args.outfile,'w') + + + if os.path.exists(args.outfile): + os.rename(os.path.abspath(args.outfile),os.path.abspath(args.outfile)+".backup") + outfile = open(args.outfile,'w') + else: + outfile = open(args.outfile,"w") + except EnvironmentError as err: sys.exit(err) - print_keywords(logfh) + try: -#### Check whether the executables are in the path + if os.path.exists(args.infile): + infile = open(args.infile,"r") - check_executables(logfh) + except EnvironmentError as err: + sys.exit(err) -#### Read the potential, store the info in 'molecules' and prints the info in OUTFILE +#### Read and check the keywords in INFILE - read_potential(args.infile) + internal = Internal(infile, outfile) - if player['lps'] == "yes": - read_lps() + internal.read_keywords() - if player['ghosts'] == "yes": - read_ghosts() + if args.opt_continue: + internal.player.cyc = cyc - print_potential(logfh) + internal.check_keywords() + internal.print_keywords() + +# #### Check whether the executables are in the path + + # internal.check_executables() + +# #### Read the potential, store the info in 'molecules' and prints the info in OUTFILE + + internal.read_potential() + + # if internal.player.lps == "yes": + # read_lps() + + # if internal.player.ghosts == "yes": + # read_ghosts() + + internal.print_potential() #### Bring the molecules to standard orientation and prints info about them - for i in range(len(molecules)): - logfh.write("\nMolecule type {}:\n\n".format(i + 1)) - print_mol_info(molecules[i], logfh) - logfh.write(" Translating and rotating molecule to standard orientation...") - standard_orientation(molecules[i]) - logfh.write(" Done\n\n New values:\n") - print_mol_info(molecules[i], logfh) + 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...") + internal.system.molecule[i].standard_orientation() + internal.outfile.write(" Done\n\n New values:\n") + internal.system.molecule[i].print_mol_info(internal.outfile) - logfh.write(90 * "=") - logfh.write("\n") + internal.outfile.write(90 * "=") + internal.outfile.write("\n") #### Open the geoms.xyz file and prints the initial geometry if starting from zero - if player['initcyc'] == 1: + if internal.player.cyc == 1: try: geomsfh = open("geoms.xyz", "w", 1) except EnvironmentError as err: sys.exit(err) - print_geom(0, geomsfh) + internal.system.print_geom(0, geomsfh) else: try: geomsfh = open("geoms.xyz", "A", 1) except EnvironmentError as err: sys.exit(err) - - logfh.write("\nStarting the iterative process.\n") + # internal.outfile.write("\nStarting the iterative process.\n") - ## Initial position (in Bohr) - position = read_position(molecules[0]) + # ## Initial position (in Bohr) + # position = internal.system.molecule[0].read_position() - ## If restarting, read the last gradient and hessian - if player['initcyc'] > 1: - if player['qmprog'] in ("g03", "g09", "g16"): - Gaussian.read_forces("grad_hessian.dat") - Gaussian.read_hessian_fchk("grad_hessian.dat") + # ## If restarting, read the last gradient and hessian + # if internal.player.cyc > 1: + # if internal.player.qmprog in ("g03", "g09", "g16"): + # Gaussian.read_forces("grad_hessian.dat") + # Gaussian.read_hessian_fchk("grad_hessian.dat") - #if player['qmprog'] == "molcas": - #Molcas.read_forces("grad_hessian.dat") - #Molcas.read_hessian("grad_hessian.dat") + # #if player['qmprog'] == "molcas": + # #Molcas.read_forces("grad_hessian.dat") + # #Molcas.read_hessian("grad_hessian.dat") - #### - #### Start the iterative process - #### + # #### + # #### Start the iterative process + # #### - for cycle in range(player['initcyc'], player['initcyc'] + player['maxcyc']): +# for cycle in range(internal.player.cyc, internal.player.cyc + internal.player.maxcyc): - logfh.write("\n" + 90 * "-" + "\n") - logfh.write("{} Step # {}\n".format(40 * " ", cycle)) - logfh.write(90 * "-" + "\n\n") +# internal.outfile.write("\n" + 90 * "-" + "\n") +# internal.outfile.write("{} Step # {}\n".format(40 * " ", cycle)) +# internal.outfile.write(90 * "-" + "\n\n") - make_step_dir(cycle) +# make_step_dir(cycle) - if player['altsteps'] == 0 or cycle == 1: - dice['randominit'] = True - else: - dice['randominit'] = False +# if internal.player.altsteps == 0 or cycle == 1: +# internal.dice.randominit = True +# else: +# internal.dice.randominit = False - #### - #### Start block of parallel simulations - #### +# #### +# #### Start block of parallel simulations +# #### - procs = [] - sentinels = [] - for proc in range(1, player['nprocs'] + 1): +# procs = [] +# sentinels = [] +# for proc in range(1, internal.player.nprocs + 1): - p = Process(target=Dice.simulation_process, args=(cycle, proc, logfh)) - p.start() - procs.append(p) - sentinels.append(p.sentinel) +# p = Process(target=Dice.simulation_process, args=(cycle, proc, internal.outfile)) +# p.start() +# procs.append(p) +# sentinels.append(p.sentinel) - while procs: - finished = connection.wait(sentinels) - for proc_sentinel in finished: - i = sentinels.index(proc_sentinel) - status = procs[i].exitcode - procs.pop(i) - sentinels.pop(i) - if status != 0: - for p in procs: - p.terminate() - sys.exit(status) +# while procs: +# finished = connection.wait(sentinels) +# for proc_sentinel in finished: +# i = sentinels.index(proc_sentinel) +# status = procs[i].exitcode +# procs.pop(i) +# sentinels.pop(i) +# if status != 0: +# for p in procs: +# p.terminate() +# sys.exit(status) - for proc in range(1, player['nprocs'] + 1): - Dice.print_last_config(cycle, proc) +# for proc in range(1, internal.player.nprocs + 1): +# Dice.print_last_config(cycle, proc) - #### - #### End of parallel simulations block - #### +# #### +# #### End of parallel simulations block +# #### - ## Make ASEC - logfh.write("\nBuilding the ASEC and vdW meanfields... ") - asec_charges = populate_asec_vdw(cycle, logfh) +# ## Make ASEC +# internal.outfile.write("\nBuilding the ASEC and vdW meanfields... ") +# asec_charges = internal.populate_asec_vdw(cycle) - ## After ASEC is built, compress files bigger than 1MB - for proc in range(1, player['nprocs'] + 1): - path = "step{:02d}".format(cycle) + os.sep + "p{:02d}".format(proc) - compress_files_1mb(path) +# ## After ASEC is built, compress files bigger than 1MB +# for proc in range(1, internal.player.nprocs + 1): +# path = "step{:02d}".format(cycle) + os.sep + "p{:02d}".format(proc) +# compress_files_1mb(path) - #### - #### Start QM calculation - #### +# #### +# #### Start QM calculation +# #### - make_qm_dir(cycle) +# make_qm_dir(cycle) - if player['opt'] == "yes": +# if internal.player.opt == "yes": - ## - ## Gaussian block - ## - if player['qmprog'] in ("g03", "g09", "g16"): +# ## +# ## Gaussian block +# ## +# if internal.player.qmprog in ("g03", "g09", "g16"): - if cycle > 1: - src = "step{:02d}".format(cycle - 1) + os.sep + "qm" + os.sep + "asec.chk" - dst = "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec.chk" - shutil.copyfile(src, dst) +# if cycle > 1: +# src = "step{:02d}".format(cycle - 1) + os.sep + "qm" + os.sep + "asec.chk" +# dst = "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec.chk" +# shutil.copyfile(src, dst) - Gaussian.make_force_input(cycle, asec_charges) - Gaussian.run_gaussian(cycle, "force", logfh) - Gaussian.run_formchk(cycle, logfh) +# Gaussian.make_force_input(cycle, asec_charges) +# Gaussian.run_gaussian(cycle, "force", internal.outfile) +# Gaussian.run_formchk(cycle, internal.outfile) - ## Read the gradient - file = "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec.fchk" - gradient = Gaussian.read_forces(file, logfh) - if len(cur_gradient) > 0: - old_gradient = cur_gradient - cur_gradient = gradient +# ## Read the gradient +# file = "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec.fchk" +# gradient = Gaussian.read_forces(file, internal.outfile) +# if len(cur_gradient) > 0: +# old_gradient = cur_gradient +# cur_gradient = gradient - ## If 1st step, read the hessian - if cycle == 1: - if player['readhessian'] == "yes": - file = "grad_hessian.dat" - logfh.write("\nReading the hessian matrix from file {}\n".format(file)) - hessian = Gaussian.read_hessian_fchk(file) - else: - file = "step01" + os.sep + "qm" + os.sep + "asec.fchk" - logfh.write("\nReading the hessian matrix from file {}\n".format(file)) - hessian = Gaussian.read_hessian(file) +# ## If 1st step, read the hessian +# if cycle == 1: +# if internal.player.readhessian == "yes": +# file = "grad_hessian.dat" +# internal.outfile.write("\nReading the hessian matrix from file {}\n".format(file)) +# hessian = Gaussian.read_hessian_fchk(file) +# else: +# file = "step01" + os.sep + "qm" + os.sep + "asec.fchk" +# internal.outfile.write("\nReading the hessian matrix from file {}\n".format(file)) +# hessian = internal.gaussian.read_hessian(file) - ## From 2nd step on, update the hessian - else: - logfh.write("\nUpdating the hessian matrix using the BFGS method... ") - hessian = update_hessian(step, cur_gradient, old_gradient, hessian) - logfh.write("Done\n") +# ## From 2nd step on, update the hessian +# else: +# internal.outfile.write("\nUpdating the hessian matrix using the BFGS method... ") +# hessian = internal.system.molecule[0].update_hessian(step, cur_gradient, old_gradient, hessian) +# internal.outfile.write("Done\n") - ## Save gradient and hessian - Gaussian.print_grad_hessian(cycle, cur_gradient, hessian) +# ## Save gradient and hessian +# internal.gaussian.print_grad_hessian(cycle, cur_gradient, hessian) - ## Calculate the step and update the position - step = calculate_step(cur_gradient, hessian, logfh) - position += step +# ## Calculate the step and update the position +# step = internal.calculate_step(cur_gradient, hessian, internal.outfile) +# position += step - ## Update the geometry of the reference molecule - update_molecule(position, logfh) +# ## Update the geometry of the reference molecule +# internal.system.molecule[0].update_molecule(position, internal.outfile) - ## If needed, calculate the charges - if cycle < player['switchcyc']: +# ## If needed, calculate the charges +# if cycle < internal.player.switchcyc: - Gaussian.make_charge_input(cycle, asec_charges) - Gaussian.run_gaussian(cycle, "charge", logfh) +# internal.gaussian.make_charge_input(cycle, asec_charges) +# internal.gaussian.run_gaussian(cycle, "charge", internal.outfile) - ## Read the new charges and update molecules[0] - if cycle < player['switchcyc']: - file = "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec2.log" - Gaussian.read_charges(file, logfh) - else: - file = "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec.log" - Gaussian.read_charges(file, logfh) +# ## Read the new charges and update molecules[0] +# if cycle < internal.player.switchcyc: +# file = "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec2.log" +# internal.gaussian.read_charges(file, internal.outfile) +# else: +# file = "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec.log" +# internal.gaussian.read_charges(file, internal.outfile) - ## Print new info for molecule[0] - logfh.write("\nNew values for molecule type 1:\n\n") - print_mol_info(molecules[0], logfh) +# ## Print new info for molecule[0] +# internal.outfile.write("\nNew values for molecule type 1:\n\n") +# internal.system.molecule[0].print_mol_info() - ## Print new geometry in geoms.xyz - print_geom(cycle, geomsfh) +# ## Print new geometry in geoms.xyz +# internal.system.molecule[0].print_geom(cycle, geomsfh) - ## - ## Molcas block - ## - #if player['qmprog'] == "molcas": +# ## +# ## Molcas block +# ## +# #if player['qmprog'] == "molcas": - #elif player['opt'] == "ts": +# #elif player['opt'] == "ts": - ## - ## Gaussian block - ## - #if player['qmprog'] in ("g03", "g09", "g16"): +# ## +# ## Gaussian block +# ## +# #if player['qmprog'] in ("g03", "g09", "g16"): - ## - ## Molcas block - ## - #if player['qmprog'] == "molcas": +# ## +# ## Molcas block +# ## +# #if player['qmprog'] == "molcas": - else: ## Only relax the charge distribution +# else: ## Only relax the charge distribution - if player['qmprog'] in ("g03", "g09", "g16"): +# if internal.player.qmprog in ("g03", "g09", "g16"): - if cycle > 1: - src = "step{:02d}".format(cycle - 1) + os.sep + "qm" + os.sep + "asec.chk" - dst = "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec.chk" - shutil.copyfile(src, dst) +# if cycle > 1: +# src = "step{:02d}".format(cycle - 1) + os.sep + "qm" + os.sep + "asec.chk" +# dst = "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec.chk" +# shutil.copyfile(src, dst) - Gaussian.make_charge_input(cycle, asec_charges) - Gaussian.run_gaussian(cycle, "charge", logfh) +# Gaussian.make_charge_input(cycle, asec_charges) +# Gaussian.run_gaussian(cycle, "charge", internal.outfile) - file = "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec2.log" - Gaussian.read_charges(file) +# file = "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec2.log" +# Gaussian.read_charges(file) - ## Print new info for molecule[0] - logfh.write("\nNew values for molecule type 1:\n\n") - print_mol_info(molecules[0], logfh) +# ## Print new info for molecule[0] +# internal.outfile.write("\nNew values for molecule type 1:\n\n") +# internal.system.molecule[0].print_mol_info() - #if player['qmprog'] == "molcas": +# #if player['qmprog'] == "molcas": - #### - #### End of the iterative process - #### +# #### +# #### End of the iterative process +# #### -## imprimir ultimas mensagens, criar um arquivo de potencial para ser usado em eventual -## continuacao, fechar arquivos (geoms.xyz, run.log, ...) +# ## imprimir ultimas mensagens, criar um arquivo de potencial para ser usado em eventual +# ## continuacao, fechar arquivos (geoms.xyz, run.log, ...) - logfh.write("\nDiceplayer finished normally!\n") - logfh.close() -#### -#### End of the program -#### \ No newline at end of file +# internal.outfile.write("\nDiceplayer finished normally!\n") +# internal.outfile.close() +# #### +# #### End of the program +# #### \ No newline at end of file diff --git a/run.log b/run.log index 0be58a8..45c1227 100644 --- a/run.log +++ b/run.log @@ -2,34 +2,38 @@ ############# Welcome to DICEPLAYER version 1.0 ############# ########################################################################################## -Your python version is 3.6.1 (default, Apr 24 2017, 06:18:27) -[GCC 4.2.1 Compatible Apple LLVM 8.0.0 (clang-800.0.42.1)] +Your python version is 3.8.8 (default, Apr 13 2021, 19:58:26) +[GCC 7.3.0] -Program started on Wednesday, 14 Jun 2017 at 12:30:02 +Program started on Saturday, 25 Sep 2021 at 15:24:31 Environment variables: -OMP_STACKSIZE = 32M +OMP_STACKSIZE = Not set ========================================================================================== CONTROL variables being used in this run: ------------------------------------------------------------------------------------------ altsteps = 20000 +cyc = 1 +freq = no ghosts = no -initcyc = 1 lps = no maxcyc = 3 maxstep = 0.3 nprocs = 2 opt = no qmprog = g09 -switch_step = 3 -zipprog = gzip +readhessian = no +switchcyc = 3 +tol_factor = 1.2 +vdwforces = no ------------------------------------------------------------------------------------------ DICE variables being used in this run: ------------------------------------------------------------------------------------------ +combrule = * dens = 0.75 isave = 1000 ljname = phb.pot @@ -38,6 +42,7 @@ nmol = 1 100 nstep = 40000 60000 50000 outname = phb press = 1.0 +progname = dice temp = 300.0 title = Diceplayer run @@ -45,15 +50,13 @@ 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 +qmprog = g09 -========================================================================================== - -Program dice found at /usr/local/bin/dice - ========================================================================================== Potential parameters from file phb.pot: ------------------------------------------------------------------------------------------ @@ -136,11 +139,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 @@ -162,68 +165,14 @@ Molecule type 2: Translating and rotating molecule to standard orientation... Done New values: - 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 ) + 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 ) Total mass = 112.20 au Total charge = -0.0000 e - Dipole moment = ( 1.7575 , 1.5369 , -0.0000 ) Total = 2.3347 Debye + Dipole moment = ( 1.8148 , 1.4689 , -0.0016 ) Total = 2.3347 Debye ========================================================================================== - -Starting the iterative process. - ------------------------------------------------------------------------------------------- - Step # 1 ------------------------------------------------------------------------------------------- - -Simulation process p01 initiated with pid 6822 -Simulation process p02 initiated with pid 6823 -p01> NVT thermalization initiated (from random configuration) on 14 Jun 2017 at 12:30:02 -p02> NVT thermalization initiated (from random configuration) on 14 Jun 2017 at 12:30:02 -p02> NPT thermalization initiated on 14 Jun 2017 at 12:41:16 -p01> NPT thermalization initiated on 14 Jun 2017 at 12:41:17 -p02> NPT production initiated on 14 Jun 2017 at 13:01:25 -p01> NPT production initiated on 14 Jun 2017 at 13:01:27 -p02> ----- NPT production finished on 14 Jun 2017 at 13:57:41 -p01> ----- NPT production finished on 14 Jun 2017 at 13:57:53 - -Building the ASEC and vdW meanfields... In average, 99.97 molecules were selected from the production simulations to form the -ASEC comprising a shell with minimum thickness of 15.352263400006278 AngstromDone. - ------------------------------------------------------------------------------------------- - Step # 2 ------------------------------------------------------------------------------------------- - -Simulation process p01 initiated with pid 6903 -Simulation process p02 initiated with pid 6904 -p02> NPT thermalization initiated (from previous configuration) on 14 Jun 2017 at 13:58:52 -p01> NPT thermalization initiated (from previous configuration) on 14 Jun 2017 at 13:58:52 -p02> NPT production initiated on 14 Jun 2017 at 14:05:14 -p01> NPT production initiated on 14 Jun 2017 at 14:05:20 -p02> ----- NPT production finished on 14 Jun 2017 at 14:23:58 -p01> ----- NPT production finished on 14 Jun 2017 at 14:24:05 - -Building the ASEC and vdW meanfields... In average, 99.66 molecules were selected from the production simulations to form the -ASEC comprising a shell with minimum thickness of 14.942583400006276 AngstromDone. - ------------------------------------------------------------------------------------------- - Step # 3 ------------------------------------------------------------------------------------------- - -Simulation process p01 initiated with pid 6945 -p01> NPT thermalization initiated (from previous configuration) on 14 Jun 2017 at 14:24:45 -Simulation process p02 initiated with pid 6946 -p02> NPT thermalization initiated (from previous configuration) on 14 Jun 2017 at 14:24:45 -p02> NPT production initiated on 14 Jun 2017 at 14:31:05 -p01> NPT production initiated on 14 Jun 2017 at 14:31:06 -p02> ----- NPT production finished on 14 Jun 2017 at 14:48:15 -p01> ----- NPT production finished on 14 Jun 2017 at 14:48:15 - -Building the ASEC and vdW meanfields... In average, 99.0 molecules were selected from the production simulations to form the -ASEC comprising a shell with minimum thickness of 14.79446440000628 AngstromDone. - -Diceplayer finished normally! \ No newline at end of file diff --git a/run.log.backup b/run.log.backup new file mode 100644 index 0000000..1f9ac41 --- /dev/null +++ b/run.log.backup @@ -0,0 +1,180 @@ +########################################################################################## +############# Welcome to DICEPLAYER version 1.0 ############# +########################################################################################## + +Your python version is 3.8.8 (default, Apr 13 2021, 19:58:26) +[GCC 7.3.0] + +Program started on Saturday, 25 Sep 2021 at 15:23:44 + +Environment variables: +OMP_STACKSIZE = Not set + +========================================================================================== + CONTROL variables being used in this run: +------------------------------------------------------------------------------------------ + +altsteps = 20000 +cyc = 1 +freq = no +ghosts = no +lps = no +maxcyc = 3 +maxstep = 0.3 +nprocs = 2 +opt = no +qmprog = g09 +readhessian = no +switchcyc = 3 +tol_factor = 1.2 +vdwforces = no + +------------------------------------------------------------------------------------------ + DICE variables being used in this run: +------------------------------------------------------------------------------------------ + +combrule = * +dens = 0.75 +isave = 1000 +ljname = phb.pot +ncores = 1 +nmol = 1 100 +nstep = 40000 60000 50000 +outname = phb +press = 1.0 +progname = dice +temp = 300.0 +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 +qmprog = g09 + + +========================================================================================== + Potential parameters from file phb.pot: +------------------------------------------------------------------------------------------ + +Combination rule: * +Types of molecules: 2 + +31 atoms in molecule type 1: +--------------------------------------------------------------------------------- +Lbl AN X Y Z Charge Epsilon Sigma Mass +--------------------------------------------------------------------------------- +1 6 -4.40034 0.66551 0.41431 -0.387593 0.07000 3.5500 12.0110 +1 6 -4.49957 -0.67279 -0.15327 0.804001 0.07000 3.5500 12.0110 +1 6 -3.25630 -1.28922 -0.61351 -0.338491 0.07000 3.5500 12.0110 +1 6 -2.05849 -0.66807 -0.48674 -0.220783 0.07000 3.5500 12.0110 +1 6 -1.97114 0.65148 0.10511 0.563527 0.07000 3.5500 12.0110 +1 6 -3.20601 1.29684 0.49712 -0.145372 0.07000 3.5500 12.0110 +2 8 -5.59453 -1.25309 -0.28765 -0.714461 0.21000 2.9600 15.9990 +3 7 -0.87404 1.35870 0.27636 -0.721792 0.17000 3.2500 14.0070 +4 6 0.38785 0.81512 0.14410 0.607123 0.07000 3.5500 12.0110 +4 6 1.40776 1.65800 -0.33261 -0.295476 0.07000 3.5500 12.0110 +4 6 0.75950 -0.46815 0.58791 -0.346661 0.07000 3.5500 12.0110 +4 6 2.71021 1.21519 -0.46119 -0.298546 0.07000 3.5500 12.0110 +4 6 2.07354 -0.89884 0.50831 -0.234284 0.07000 3.5500 12.0110 +4 6 3.08076 -0.08936 -0.06112 0.294651 0.07000 3.5500 12.0110 +5 7 4.37238 -0.53979 -0.20183 -0.190497 0.17000 3.2500 14.0070 +6 6 5.41980 0.43824 -0.43836 -0.213183 0.06600 3.5000 12.0110 +6 6 4.76361 -1.73522 0.52225 -0.225420 0.06600 3.5000 12.0110 +7 1 6.36700 -0.08536 -0.51962 0.109840 0.03000 2.5000 1.0079 +7 1 5.25196 0.96612 -1.37421 0.077140 0.03000 2.5000 1.0079 +7 1 5.49517 1.17412 0.36749 0.125758 0.03000 2.5000 1.0079 +7 1 4.15838 -2.58442 0.21446 0.083229 0.03000 2.5000 1.0079 +7 1 5.79695 -1.96254 0.28041 0.115361 0.03000 2.5000 1.0079 +7 1 4.67416 -1.61707 1.60661 0.121440 0.03000 2.5000 1.0079 +7 1 0.02813 -1.09554 1.08007 0.179434 0.03000 2.4200 1.0079 +7 1 2.31876 -1.86576 0.91913 0.177250 0.03000 2.4200 1.0079 +7 1 3.45102 1.89880 -0.84669 0.198844 0.03000 2.4200 1.0079 +7 1 1.14759 2.66909 -0.61498 0.173388 0.03000 2.4200 1.0079 +7 1 -3.11930 2.29679 0.90060 0.151534 0.03000 2.4200 1.0079 +7 1 -5.31582 1.13611 0.74479 0.205289 0.03000 2.4200 1.0079 +7 1 -3.33813 -2.25037 -1.10231 0.184367 0.03000 2.4200 1.0079 +7 1 -1.16374 -1.11791 -0.89380 0.160382 0.03000 2.4200 1.0079 + +16 atoms in molecule type 2: +--------------------------------------------------------------------------------- +Lbl AN X Y Z Charge Epsilon Sigma Mass +--------------------------------------------------------------------------------- +1 6 0.67203 -2.82345 0.00263 -0.115000 0.07000 3.5500 12.0110 +1 6 2.07203 -2.82345 0.00263 -0.115000 0.07000 3.5500 12.0110 +1 6 2.76823 -1.61764 0.00263 -0.115000 0.07000 3.5500 12.0110 +1 6 2.06824 -0.40521 0.00264 -0.115000 0.07000 3.5500 12.0110 +1 14 0.67589 -0.40522 0.00264 0.150000 0.07000 3.5500 28.0860 +1 104 -0.02420 -1.61760 0.00263 -0.115000 0.07000 3.5500 0.0000 +2 1 0.13203 -3.75875 0.00263 0.115000 0.03000 2.4200 1.0079 +2 1 2.61203 -3.75875 0.00263 0.115000 0.03000 2.4200 1.0079 +2 1 2.60824 0.53010 0.00264 0.115000 0.03000 2.4200 1.0079 +2 1 -1.10420 -1.61760 0.00263 0.115000 0.03000 2.4200 1.0079 +3 8 -0.00411 0.77257 0.00264 -0.585000 0.17000 3.0700 15.9990 +3 1 0.61978 1.50220 0.00264 0.435000 0.00000 0.0000 1.0079 +4 6 4.27823 -1.61764 0.00263 0.115000 0.17000 3.8000 12.0110 +4 1 4.63490 -0.74399 0.50704 0.000000 0.00000 0.0000 1.0079 +4 1 4.63490 -2.49130 0.50704 0.000000 0.00000 0.0000 1.0079 +4 1 4.63490 -1.61764 -1.00617 0.000000 0.00000 0.0000 1.0079 + + +========================================================================================== + +Molecule type 1: + + Center of mass = ( -0.0000 , 0.0000 , 0.0000 ) + Moments of inertia = 3.144054E+02 2.801666E+03 3.027366E+03 + Major principal axis = ( 0.999972 , 0.007210 , 0.002184 ) + Inter principal axis = ( -0.007218 , 0.999967 , 0.003660 ) + Minor principal axis = ( -0.002157 , -0.003676 , 0.999991 ) + Characteristic lengths = ( 11.96 , 5.25 , 2.98 ) + Total mass = 226.28 au + Total charge = -0.0000 e + Dipole moment = ( 9.8367 , 0.6848 , 0.8358 ) Total = 9.8959 Debye + + Translating and rotating molecule to standard orientation... Done + + New values: + 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 ) + Characteristic lengths = ( 11.97 , 5.27 , 2.99 ) + Total mass = 226.28 au + Total charge = -0.0000 e + Dipole moment = ( 9.8432 , 0.6168 , 0.8120 ) Total = 9.8959 Debye + + +Molecule type 2: + + Center of mass = ( 1.6067 , -1.0929 , 0.0026 ) + Moments of inertia = 1.205279E+02 2.859254E+02 4.033761E+02 + Major principal axis = ( 0.795913 , -0.605411 , -0.000001 ) + Inter principal axis = ( 0.605411 , 0.795913 , 0.000001 ) + Minor principal axis = ( 0.000000 , -0.000002 , 1.000000 ) + Characteristic lengths = ( 5.74 , 5.26 , 1.51 ) + Total mass = 112.20 au + Total charge = -0.0000 e + Dipole moment = ( 2.3293 , 0.1593 , -0.0000 ) Total = 2.3347 Debye + + 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 ) + Total mass = 112.20 au + Total charge = -0.0000 e + Dipole moment = ( 1.8148 , 1.4689 , -0.0016 ) Total = 2.3347 Debye + +========================================================================================== + +Starting the iterative process.