diff --git a/Pipfile b/Pipfile index 1434829..26334cf 100644 --- a/Pipfile +++ b/Pipfile @@ -8,6 +8,8 @@ numpy = "*" pickle5 = "*" argparse = "*" pyinstaller = "*" +setproctitle = "*" +pyyaml = "*" [dev-packages] diff --git a/diceplayer/DPpack/Environment/Molecule.py b/diceplayer/DPpack/Environment/Molecule.py index 4269211..2f41235 100644 --- a/diceplayer/DPpack/Environment/Molecule.py +++ b/diceplayer/DPpack/Environment/Molecule.py @@ -234,6 +234,11 @@ class Molecule: return position + def updateCharges(self, charges: List[float]) -> None: + + for i, atom in enumerate(self.atom): + atom.chg = charges[i] + def update_hessian( self, step: np.ndarray, @@ -418,4 +423,4 @@ class Molecule: dz = atom1.rz - atom2.rz distances.append(math.sqrt(dx**2 + dy**2 + dz**2)) - return min(distances) + return min(distances) \ No newline at end of file diff --git a/diceplayer/DPpack/Environment/System.py b/diceplayer/DPpack/Environment/System.py index 9e5d52a..4783938 100644 --- a/diceplayer/DPpack/Environment/System.py +++ b/diceplayer/DPpack/Environment/System.py @@ -222,3 +222,23 @@ class System: symbol, atom.rx, atom.ry, atom.rz ) ) + + def printChargesAndDipole(self, cycle: int, fh: TextIO) -> None: + """ + Print the charges and dipole of the molecule in the Output file + + Args: + cycle (int): Number of the cycle + fh (TextIO): Output file + """ + + fh.write("Cycle # {}\n".format(cycle)) + fh.write("Number of site: {}\n".format(len(self.molecule[0].atom))) + + chargesAndDipole = self.molecule[0].charges_and_dipole() + + fh.write( + "{:>10.6f} {:>10.6f} {:>10.6f} {:>10.6f} {:>10.6f}\n".format( + chargesAndDipole[0], chargesAndDipole[1], chargesAndDipole[2], chargesAndDipole[3], chargesAndDipole[4] + ) + ) diff --git a/diceplayer/DPpack/External/Dice.py b/diceplayer/DPpack/External/Dice.py index 5b477f8..dbc6b7c 100644 --- a/diceplayer/DPpack/External/Dice.py +++ b/diceplayer/DPpack/External/Dice.py @@ -158,7 +158,7 @@ class Dice: fh.write("ljname = {}\n".format(self.ljname)) fh.write("outname = {}\n".format(self.outname)) - string = " ".join(str(x) for x in self.step.nmol) + string = " ".join(str(x) for x in self.nmol) fh.write("nmol = {}\n".format(string)) fh.write("dens = {}\n".format(self.dens)) @@ -197,7 +197,7 @@ class Dice: fh.write("ljname = {}\n".format(self.ljname)) fh.write("outname = {}\n".format(self.outname)) - string = " ".join(str(x) for x in self.step.nmol) + string = " ".join(str(x) for x in self.nmol) fh.write("nmol = {}\n".format(string)) fh.write("dens = {}\n".format(self.dens)) @@ -229,7 +229,7 @@ class Dice: fh.write("ljname = {}\n".format(self.ljname)) fh.write("outname = {}\n".format(self.outname)) - string = " ".join(str(x) for x in self.step.nmol) + string = " ".join(str(x) for x in self.nmol) fh.write("nmol = {}\n".format(string)) fh.write("press = {}\n".format(self.press)) @@ -268,7 +268,7 @@ class Dice: fh.write("ljname = {}\n".format(self.ljname)) fh.write("outname = {}\n".format(self.outname)) - string = " ".join(str(x) for x in self.step.nmol) + string = " ".join(str(x) for x in self.nmol) fh.write("nmol = {}\n".format(string)) fh.write("press = {}\n".format(self.press)) diff --git a/diceplayer/DPpack/External/Gaussian.py b/diceplayer/DPpack/External/Gaussian.py index 4062342..43b9ee1 100644 --- a/diceplayer/DPpack/External/Gaussian.py +++ b/diceplayer/DPpack/External/Gaussian.py @@ -1,9 +1,11 @@ +from ast import keyword +from asyncore import read import os import shutil import subprocess import sys import textwrap -from typing import TextIO +from typing import Dict, List, TextIO import numpy as np @@ -23,12 +25,20 @@ class Gaussian: gbottom = None # inserted in the gaussian input pop = "chelpg" + keywords = "" + def __init__(self) -> None: pass @NotNull(requiredArgs=["qmprog","level"]) def updateKeywords(self, **data): self.__dict__.update(**data) + self.checkKeywords() + + def checkKeywords(self): + + if self.pop not in ["chelpg", "mk", "nbo"]: + self.pop = "chelpg" def run_formchk(self, cycle: int, fh: TextIO): @@ -47,6 +57,29 @@ class Gaussian: os.chdir(work_dir) + def readChargesFromFchk(self, file: str, fh: TextIO) -> List[float]: + + try: + with open(file) as fchk: + fchkfile = fchk.readlines() + except: + sys.exit("Error: cannot open file {}".format(file)) + + if self.pop in ["chelpg", "mk"]: + CHARGE_FLAG = "ESP Charges" + else: + CHARGE_FLAG = "ESP Charges" + + start = fchkfile.pop(0).strip() + while start.find(CHARGE_FLAG) != 0: # expression in begining of line + start = fchkfile.pop(0).strip() + + charges: List[float] = [] + while len(charges) < len(self.step.molecule[0].atom): + charges.extend([float(x) for x in fchkfile.pop(0).split()]) + + return charges + def read_forces_fchk(self, file: str, fh: TextIO) -> np.ndarray: forces = [] @@ -202,7 +235,7 @@ class Gaussian: fh.close() # Change the name to make_gaussian_input - def make_gaussian_input(self, cycle: int, asec_charges=None) -> None: + def make_gaussian_input(self, cycle: int, asec_charges: List[Dict]) -> None: simdir = "simfiles" stepdir = "step{:02d}".format(cycle) @@ -222,7 +255,7 @@ class Gaussian: kword_line = "#P " + str(self.level) - if self.keywords != None: + if self.keywords != "": kword_line += " " + self.keywords if self.step.opt == "yes": @@ -250,58 +283,18 @@ class Gaussian: ) ) - # ## If also performing charge fit in the same calculation - # if cycle >= self.player.switchcyc: - # for ghost in ghost_atoms: - # fh.write("Bq {:>10.5f} {:>10.5f} {:>10.5f}\n".format( - # ghost['rx'], ghost['ry'], ghost['rz'])) - - # for lp in lp_atoms: - # fh.write("Bq {:>10.5f} {:>10.5f} {:>10.5f}\n".format( - # lp['rx'], lp['ry'], lp['rz'])) - - # fh.write("\n") - - # If gmiddle file was informed, write its contents in asec.gjf - # if self.gmiddle != None: - # if not os.path.isfile(self.gmiddle): - # sys.exit("Error: cannot find file {} in main directory".format( - # self.gmiddle)) - # try: - # with open(self.gmiddle) as gmiddlefile: - # gmiddle = gmiddlefile.readlines() - # except: - # sys.exit("Error: cannot open file {}".format(self.gmiddle)) - - # for line in gmiddle: - # fh.write(line) - - # fh.write("\n") - - # ## Write the ASEC: - # for charge in asec_charges: - # fh.write("{:>10.5f} {:>10.5f} {:>10.5f} {:>11.8f}\n".format( - # charge['rx'], charge['ry'], charge['rz'], charge['chg'])) - fh.write("\n") - # ## If gbottom file was informed, write its contents in asec.gjf - # if self.gbottom != None: - # if not os.path.isfile(self.gbottom): - # sys.exit("Error: cannot find file {} in main directory".format( - # self.gbottom)) - # try: - # with open(self.gbottom) as gbottomfile: - # gbottom = gbottomfile.readlines() - # except: - # sys.exit("Error: cannot open file {}".format(self.gbottom)) + for charge in asec_charges: + fh.write( + "{:>10.5f} {:>10.5f} {:>10.5f} {:>11.8f}\n".format( + charge['rx'], charge['ry'], charge['rz'], charge['chg'] + ) + ) - # for line in gbottom: - # fh.write(line) - - # fh.write("\n") - - # fh.close() + fh.write("\n") + + fh.close() def read_charges(self, file: str, fh: TextIO) -> None: @@ -343,6 +336,90 @@ class Gaussian: fh.write("------------------------------------\n") + def executeOptimizationRoutine(self, cycle: int, outfile: TextIO, readhessian: str): + + try: + gradient + old_gradient = gradient + except: + pass + + gradient = self.read_forces_fchk(file, outfile) + + # If 1st step, read the hessian + if cycle == 1: + + if readhessian == "yes": + + file = "grad_hessian.dat" + outfile.write( + "\nReading the hessian matrix from file {}\n".format(file) + ) + hessian = self.read_hessian_log(file) + + else: + + file = ( + "simfiles" + + os.sep + + "step01" + + os.sep + + "qm" + + os.sep + + "asec.fchk" + ) + outfile.write( + "\nReading the hessian matrix from file {}\n".format(file) + ) + hessian = self.read_hessian_fchk(file) + + # From 2nd step on, update the hessian + else: + outfile.write("\nUpdating the hessian matrix using the BFGS method... ") + hessian = self.step.molecule[0].update_hessian( + step, gradient, old_gradient, hessian + ) + outfile.write("Done\n") + + # Save gradient and hessian + self.print_grad_hessian(cycle, gradient, hessian) + + # Calculate the step and update the position + step = self.calculate_step(cycle, gradient, hessian) + + position += step + + ## If needed, calculate the charges + if cycle < self.step.switchcyc: + + # internal.gaussian.make_charge_input(cycle, asec_charges) + self.run_gaussian(cycle, "charge", outfile) + + file = ( + "simfiles" + + os.sep + + "step{:02d}".format(cycle) + + os.sep + + "qm" + + os.sep + + "asec2.log" + ) + self.read_charges(file, outfile) + else: + file = ( + "simfiles" + + os.sep + + "step{:02d}".format(cycle) + + os.sep + + "qm" + + os.sep + + "asec.log" + ) + self.read_charges(file, outfile) + + self.outfile.write("\nNew values for molecule type 1:\n\n") + self.step.molecule[0].print_mol_info(outfile) + def run_gaussian(self, cycle: int, type: str, fh: TextIO) -> None: simdir = "simfiles" @@ -459,41 +536,39 @@ class Gaussian: self.step = step - def start(self, cycle: int, outfile: TextIO, readhessian: str) -> np.ndarray: + def start(self, cycle: int, outfile: TextIO, asec_charges: List[Dict], readhessian: str) -> StepDTO: make_qm_dir(cycle) - if self.qmprog in ("g03", "g09", "g16"): + if cycle > 1: - if cycle > 1: - - src = ( - "simfiles" - + os.sep - + "step{:02d}".format(cycle - 1) - + os.sep - + "qm" - + os.sep - + "asec.chk" - ) - dst = ( - "simfiles" - + os.sep - + "step{:02d}".format(cycle) - + os.sep - + "qm" - + os.sep - + "asec.chk" - ) - shutil.copyfile(src, dst) - - self.make_gaussian_input(cycle) - self.run_gaussian(cycle, "force", outfile) - self.run_formchk(cycle, outfile) - - ## Read the gradient - file = ( + src = ( "simfiles" + + os.sep + + "step{:02d}".format(cycle - 1) + + os.sep + + "qm" + + os.sep + + "asec.chk" + ) + dst = ( + "simfiles" + + os.sep + + "step{:02d}".format(cycle) + + os.sep + + "qm" + + os.sep + + "asec.chk" + ) + shutil.copyfile(src, dst) + + self.make_gaussian_input(cycle, asec_charges) + self.run_gaussian(cycle, "force", outfile) + self.run_formchk(cycle, outfile) + + ## Read the gradient + file = ( + "simfiles" + os.sep + "step{:02d}".format(cycle) + os.sep @@ -502,152 +577,18 @@ class Gaussian: + "asec.fchk" ) - try: - gradient - old_gradient = gradient - except: - pass + if self.step.opt: - gradient = self.read_forces_fchk(file, outfile) + pass + # position = self.executeOptimizationRoutine(cycle, outfile, readhessian) + # self.step.position = position - # If 1st step, read the hessian - if cycle == 1: + else: - if readhessian == "yes": + charges = self.readChargesFromFchk(file, outfile) + self.step.charges = charges - file = "grad_hessian.dat" - outfile.write( - "\nReading the hessian matrix from file {}\n".format(file) - ) - hessian = self.read_hessian_log(file) - - else: - - file = ( - "simfiles" - + os.sep - + "step01" - + os.sep - + "qm" - + os.sep - + "asec.fchk" - ) - outfile.write( - "\nReading the hessian matrix from file {}\n".format(file) - ) - hessian = self.read_hessian_fchk(file) - - # From 2nd step on, update the hessian - else: - outfile.write("\nUpdating the hessian matrix using the BFGS method... ") - hessian = self.step.molecule[0].update_hessian( - step, gradient, old_gradient, hessian - ) - outfile.write("Done\n") - - # Save gradient and hessian - self.print_grad_hessian(cycle, gradient, hessian) - - # Calculate the step and update the position - step = self.calculate_step(cycle, gradient, hessian) - - position += step - - ## If needed, calculate the charges - if cycle < self.step.switchcyc: - - # internal.gaussian.make_charge_input(cycle, asec_charges) - self.run_gaussian(cycle, "charge", outfile) - - file = ( - "simfiles" - + os.sep - + "step{:02d}".format(cycle) - + os.sep - + "qm" - + os.sep - + "asec2.log" - ) - self.read_charges(file, outfile) - else: - file = ( - "simfiles" - + os.sep - + "step{:02d}".format(cycle) - + os.sep - + "qm" - + os.sep - + "asec.log" - ) - self.read_charges(file, outfile) - - ## Print new info for molecule[0] - self.outfile.write("\nNew values for molecule type 1:\n\n") - self.step.molecule[0].print_mol_info(outfile) - - ## - ## Molcas block - ## - # if player['qmprog'] == "molcas": - - # elif player['opt'] == "ts": - - ## - ## Gaussian block - ## - # if player['qmprog'] in ("g03", "g09", "g16"): - - ## - ## Molcas block - ## - # if player['qmprog'] == "molcas": - - # else: ## Only relax the charge distribution - - # if internal.player.qmprog in ("g03", "g09", "g16"): - - # if cycle > 1: - # src = ( - # "simfiles" - # + os.sep - # + "step{:02d}".format(cycle - 1) - # + os.sep - # + "qm" - # + os.sep - # + "asec.chk" - # ) - # dst = ( - # "simfiles" - # + os.sep - # + "step{:02d}".format(cycle) - # + os.sep - # + "qm" - # + os.sep - # + "asec.chk" - # ) - # shutil.copyfile(src, dst) - - # # internal.gaussian.make_charge_input(cycle, asec_charges) - # internal.gaussian.run_gaussian(cycle, "charge", internal.outfile) - - # file = ( - # "simfiles" - # + os.sep - # + "step{:02d}".format(cycle) - # + os.sep - # + "qm" - # + os.sep - # + "asec2.log" - # ) - # internal.read_charges(file) - - # ## 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" - - return position + return self.step def reset(self): diff --git a/diceplayer/DPpack/Player.py b/diceplayer/DPpack/Player.py index cb81833..664256e 100644 --- a/diceplayer/DPpack/Player.py +++ b/diceplayer/DPpack/Player.py @@ -2,8 +2,7 @@ import os import shutil import sys import textwrap -import types -from typing import TextIO +from typing import Dict, List, TextIO import yaml @@ -29,9 +28,15 @@ class Player: lps = None ghosts = None altsteps = None - combrule = None + switchcyc = 3 + maxstep = 0.3 + freq = "no" + readhessian = "no" + vdwforces = "no" + tol_factor = 1.2 + TOL_RMS_FORCE = 3e-4 TOL_MAX_FORCE = 4.5e-4 TOL_RMS_STEP = 1.2e-3 @@ -62,7 +67,7 @@ class Player: ] @NotNull( - requiredArgs=["maxcyc", "opt", "nprocs", "qmprog", "lps", "ghosts", "altsteps"] + requiredArgs=["maxcyc", "opt", "nprocs", "qmprog", "altsteps"] ) def updateKeywords(self, **data): self.__dict__.update(**data) @@ -590,35 +595,51 @@ class Player: self.dice.reset() def gaussian_start(self, cycle: int, geomsfh: TextIO): - + self.gaussian.configure( - self.initcyc, - self.nprocs, - self.dice.ncores, - self.altsteps, - self.switchcyc, - self.opt, - self.system.nmols, - self.system.molecule, + StepDTO( + initcyc=self.initcyc, + nprocs=self.nprocs, + ncores=self.dice.ncores, + altsteps=self.altsteps, + switchcyc=self.switchcyc, + opt=self.opt, + nmol=self.system.nmols, + molecule=self.system.molecule + ) ) - position = self.gaussian.start(cycle, self.outfile, self.readhessian) + # Make ASEC + self.outfile.write("\nBuilding the ASEC and vdW meanfields... ") + asec_charges = self.populate_asec_vdw(cycle) - ## Update the geometry of the reference molecule - self.system.update_molecule(position, self.outfile) + step = self.gaussian.start(cycle, self.outfile, asec_charges, self.readhessian) - ## Print new geometry in geoms.xyz - self.system.print_geom(cycle, geomsfh) + if self.opt: + + position = step.position + + ## Update the geometry of the reference molecule + self.system.update_molecule(position, self.outfile) + + ## Print new geometry in geoms.xyz + self.system.print_geom(cycle, geomsfh) + + else: + + charges = step.charges + + self.system.molecule[0].updateCharges(charges) + + self.system.printChargesAndDipole(cycle, self.outfile) self.gaussian.reset() - def populate_asec_vdw(self, cycle): + def populate_asec_vdw(self, cycle) -> List[Dict]: # Both asec_charges and vdw_meanfield will utilize the Molecule() class and Atoms() with some None elements - asec_charges = Molecule( - "ASEC_CHARGES" - ) + asec_charges = [] if self.dice.nstep[-1] % self.dice.isave == 0: nconfigs = round(self.dice.nstep[-1] / self.dice.isave) @@ -683,18 +704,14 @@ class Player: for mol in range(nmols): - new_molecule = Molecule(self.system.molecule[type]) - # Run over sites of each molecule - for site in range(len(self.system.molecule[types].atom)): + new_molecule = Molecule("ASEC TMP MOLECULE") + for site in range(len(self.system.molecule[type].atom)): - # new_molecule.append({}) line = xyzfile.pop(0).split() if ( line[0].title() - != atomsymb[ - self.system.molecule[type].atom[site].na.strip() - ] + != atomsymb[self.system.molecule[type].atom[site].na].strip() ): sys.exit("Error reading file {}".format(file)) @@ -702,15 +719,9 @@ class Player: Atom( self.system.molecule[type].atom[site].lbl, self.system.molecule[type].atom[site].na, - self.system.molecule[type] - .atom[site] - .float(line[1]), - self.system.molecule[type] - .atom[site] - .float(line[2]), - self.system.molecule[type] - .atom[site] - .float(line[3]), + float(line[1]), + float(line[2]), + float(line[3]), self.system.molecule[type].atom[site].chg, self.system.molecule[type].atom[site].eps, self.system.molecule[type].atom[site].sig, @@ -721,13 +732,7 @@ class Player: if dist < thickness[-1]: mol_count += 1 for atom in new_molecule.atom: - 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 + asec_charges.append({"lbl": atomsymb[atom.na], "rx": atom.rx, "ry": atom.ry, "rz": atom.rz, "chg": atom.chg}) # if self.vdwforces == "yes": # vdw_meanfield[-1]["rx"] = atom["rx"] @@ -781,13 +786,16 @@ class Player: self.outfile.write(textwrap.fill(string, 86)) self.outfile.write("\n") - otherfh = open("ASEC.dat", "w") + otherfh = open("ASEC.xyz", "w", 1) 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"] + "{} {:>10.5f} {:>10.5f} {:>10.5f}\n".format( + charge['lbl'], charge['rx'], charge['ry'], charge['rz'] ) ) otherfh.close() - return asec_charges \ No newline at end of file + for charge in asec_charges: + charge['chg'] /= norm_factor + + return asec_charges diff --git a/diceplayer/DPpack/Utils/StepDTO.py b/diceplayer/DPpack/Utils/StepDTO.py index 99c01ae..d03d43b 100644 --- a/diceplayer/DPpack/Utils/StepDTO.py +++ b/diceplayer/DPpack/Utils/StepDTO.py @@ -15,4 +15,7 @@ class StepDTO: switchcyc: int = None opt: str = None nmol: List[int] = None - molecule: List[Molecule] = None \ No newline at end of file + molecule: List[Molecule] = None + + charges: List[float] = None + position: List[float] = None \ No newline at end of file diff --git a/diceplayer/__main__.py b/diceplayer/__main__.py index 4f91fb7..247b0fa 100644 --- a/diceplayer/__main__.py +++ b/diceplayer/__main__.py @@ -9,7 +9,8 @@ import setproctitle from diceplayer.DPpack.Player import Player from diceplayer.DPpack.Utils.Misc import * -__VERSION = "dev" +__VERSION = "v0.0.1" +os.nice(+19) setproctitle.setproctitle("diceplayer-{}".format(__VERSION)) if __name__ == "__main__": @@ -178,13 +179,9 @@ if __name__ == "__main__": ### End of parallel simulations block ### - # Make ASEC - player.outfile.write("\nBuilding the ASEC and vdW meanfields... ") - asec_charges = player.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) + path = "simfiles"+os.sep+"step{:02d}".format(cycle) + os.sep + "p{:02d}".format(proc) compress_files_1mb(path) ### diff --git a/setup.py b/setup.py old mode 100644 new mode 100755 index 0c99503..925add5 --- a/setup.py +++ b/setup.py @@ -22,7 +22,7 @@ parser.add_argument( dest="clean", default=False, action="store_true", - help="Cleans the Development Environment", + help="Cleans the Development Environment" ) args = parser.parse_args() @@ -37,9 +37,18 @@ def __build(): def __clean(): - shutil.rmtree("build") - shutil.rmtree("dist") - os.remove("diceplayer.spec") + try: + + shutil.rmtree("build") + shutil.rmtree("dist") + os.remove("diceplayer.spec") + + except: + + print("Workspace clean.") + + + if __name__ == "__main__":