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

Relaxation of Molecular Charges - v0.0.1

Browse Source
This commit is contained in:
Vitor Hideyoshi
2022-08-12 09:17:37 -03:00
parent ab996196d4
commit 6fcdb0919f
9 changed files with 274 additions and 289 deletions
Download Patch File Download Diff File Expand all files Collapse all files

7
diceplayer/DPpack/Environment/Molecule.py
View File

@@ -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)

20
diceplayer/DPpack/Environment/System.py
View File

@@ -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]
)
)

8
diceplayer/DPpack/External/Dice.py vendored
View File

@@ -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))

389
diceplayer/DPpack/External/Gaussian.py vendored
View File

@@ -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):

106
diceplayer/DPpack/Player.py
View File

@@ -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
for charge in asec_charges:
charge['chg'] /= norm_factor
return asec_charges

5
diceplayer/DPpack/Utils/StepDTO.py
View File

@@ -15,4 +15,7 @@ class StepDTO:
switchcyc: int = None
opt: str = None
nmol: List[int] = None
molecule: List[Molecule] = None
molecule: List[Molecule] = None
charges: List[float] = None
position: List[float] = None