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

SetGlobas/MolHandling(fixes)/DicePlayer Translation

Browse Source 
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 <vitor.h.n.batista@gmail.com>
This commit is contained in:
Vitor Hideyoshi
2021-09-25 15:30:40 -03:00
committed by Vitor Hideyoshi
parent 2f9cecee34
commit 316acf4761
13 changed files with 725 additions and 423 deletions
Download Patch File Download Diff File Expand all files Collapse all files

380
DPpack/SetGlobals.py
View File

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