From 63d3d5f40dd5448714b44fccb4d453af22d8e1f0 Mon Sep 17 00:00:00 2001 From: Vitor Hideyoshi Date: Thu, 5 Aug 2021 11:22:06 -0300 Subject: [PATCH] SetGlobas Translation - Final (without ghosts and lp) Internal methods translations: read_potential, print_potential Signed-off-by: Vitor Hideyoshi --- DPpack/.SetGlobals.py | 813 +++++++++++++++++++++ DPpack/MolHandling.py | 2 + DPpack/SetGlobals.py | 1397 +++++++++++++++++-------------------- DPpack/SetGlobalsClass.py | 543 -------------- 4 files changed, 1437 insertions(+), 1318 deletions(-) create mode 100644 DPpack/.SetGlobals.py delete mode 100644 DPpack/SetGlobalsClass.py diff --git a/DPpack/.SetGlobals.py b/DPpack/.SetGlobals.py new file mode 100644 index 0000000..6bd1317 --- /dev/null +++ b/DPpack/.SetGlobals.py @@ -0,0 +1,813 @@ +import os, sys +import shutil +import textwrap + +from DPpack.PTable import * +from DPpack.Misc import * + +#### Global hashes that control the behaviour of Diceplayer + +player = {} +dice = {} +gaussian = {} +molcas = {} +internal = {} + +####################################################################### +#### Global parameters that MAY be given by the user #### +#### (If not given by the user, default values will be used) #### +####################################################################### + +## Diceplayer: +player['maxcyc'] = 1 +player['initcyc'] = 1 # May restart an optimization (append to geoms.xyz from start) +player['nprocs'] = 1 +player['switchcyc'] = 3 # At which step start doing only one QM calculation (geom & chg) +player['altsteps'] = 20000 # Steps for thermalization when starting from previous step +player['maxstep'] = 0.3 # Maxstep for geometry relaxation in Bohr +player['qmprog'] = "g09" +player['opt'] = "yes" +player['freq'] = "no" +player['readhessian'] = "no" +player['lps'] = "no" +player['ghosts'] = "no" +player['vdwforces'] = "no" +player['tol_factor'] = 1.2 # Factor to multiply the default tolerance values + +## Dice: +dice['title'] = "Diceplayer run" +dice['progname'] = "dice" +dice['temp'] = 300.0 +dice['press'] = 1.0 +dice['isave'] = 1000 # ASEC construction will take this into account +dice['ncores'] = 1 + +## Gaussian: +gaussian['mem'] = None +gaussian['keywords'] = None +gaussian['chgmult'] = [0, 1] +gaussian['gmiddle'] = None # In each case, if a filename is given, its content will be +gaussian['gbottom'] = None # inserted in the gaussian input +gaussian['pop'] = "chelpg" +gaussian['chglevel'] = None + +## Molcas: +molcas['orbfile'] = "input.exporb" +molcas['root'] = 1 + + +######################################################################## +#### Global parameters that MUST be given by the user #### +######################################################################## + +## Dice: +dice['dens'] = None # Investigate the possibility of using 'box = Lx Ly Lz' instead. +#dice['box'] = None # So 'geom' would be set by diceplayer and 'cutoff' would be + # switched off. One of them must be given. +dice['ljname'] = None +dice['outname'] = None +dice['nmol'] = [] # Up to 4 integer values related to up to 4 molecule types +dice['nstep'] = [] # 2 or 3 integer values related to 2 or 3 simulations + # (NVT th + NVT eq) or (NVT th + NPT th + NPT eq). + # This will control the 'nstep' keyword of Dice + +## Gaussian: +gaussian['level'] = None + +## Molcas: +molcas['mbottom'] = None +molcas['basis'] = None + + +## The following Dice keywords will be handled automatically by Diceplayer: +## * init ("yes" in the first thermalization and "yesreadxyz" for thermalizations +## starting from a previous step / "no" in subsequent simulations) +## * vstep ("0" for NVT simulations and 'nstep'/5 for NPT simulations) +## * nstep ('nstep' for NVT and "5" for NPT simulations ) +## * irdf ("0" for thermalizations and "10*nprocs" for equilibrium) +## * seed (will be generated randomly each time a Dice input is created) + +## The following Dice keywords will be set constant by Diceplayer for all simulations +## * mstop = 1 (So to guarantee that the ASEC will be correctly built) +## * accum = no (There is never a simulation continuation in Diceplayer) +## * iprint = 1 (Print energy info every step in Dice output) + +## All the other Dice keywords will not be altered from their default values +## and therefore are not mentioned in Diceplayer + + +##################################################################### +#### Global parameters that are not accessible to the user #### +#### (Intended to be used only internally by the program) #### +##################################################################### + +## Diceplayer: +internal['tol_rms_force'] = 3e-4 # Hartree/Bohr +internal['tol_max_force'] = 4.5e-4 # Hartree/Bohr +internal['tol_rms_step'] = 1.2e-3 # Bohr +internal['tol_max_step'] = 1.8e-3 # Bohr +internal['trust_radius'] = None + +## Dice: +internal['combrule'] = None +internal['randominit'] = None + +## Other global variables: +molecules = [] # Armazena todas as informacoes sobre cada tipo de molecula + # (lbl, na, rx, ry, rz, chg, eps, sig, mass) + +internal['ghost_types'] = [] +internal['ghost_atoms'] = [] # Store the ghost atoms (off-atom charge sites) in the QM molecule + # (rx, ry, rz, chg) +internal['lp_types'] = [] +internal['lp_atoms'] = [] # Store the lone pairs (special off-atom charge sites) in the QM molecule + # (rx, ry, rz, chg) + +## Numpy arrays: +step = [] ## Values in Bohr +internal['position'] = [] +internal['energy'] = [] ## Values in Hartree +internal['gradient'] = [] ## Values in Hartree/Bohr +internal['hessian'] = [] ## Values in Hartree/Bohr^2 + +## Conversion factors: +bohr2ang = 0.52917721092 +ang2bohr = 1/bohr2ang + +###################################################################### +#### Environment variables important for executing Diceplayer #### +###################################################################### + +env = ["OMP_STACKSIZE"] + +####################################### functions ###################################### +## Functions to process the input files and store the values in the global variables ## +########################################################################################## + +# def read_keywords(infile): + +# try: +# with open(infile) as fh: +# controlfile = fh.readlines() +# except EnvironmentError: +# sys.exit("Error: cannot open file {}".format(infile)) + +# for line in controlfile: + +# key, value = line.partition("=")[::2] # Discards the '=' +# key = key.strip().lower() +# if key in ('title', 'keywords'): +# value = value.strip() +# else: +# value = value.split() + +# #### Read the Diceplayer related keywords +# if key in player and len(value) != 0: ## 'value' is not empty! + +# if key == 'qmprog' and value[0].lower() in ("g03", "g09", "g16", "molcas"): +# player[key] = value[0].lower() + +# elif key == 'opt' and value[0].lower() in ("yes", "no", "ts"): +# player[key] = value[0].lower() + +# #elif key == 'zipprog' and value[0].lower() in ("zip", "gzip", "bzip"): +# #player[key] = value[0].lower() + +# elif key in ('lps', 'ghosts') and value[0].lower() in ("yes", "no"): +# player[key] = value[0].lower() + +# elif key in ('readhessian', 'vdwforces') and value[0].lower() in ("yes", "no"): +# player[key] = value[0].lower() + +# elif key in ('maxcyc', 'initcyc', 'nprocs', 'altsteps', 'switchcyc'): +# err = "Error: expected a positive integer for keyword {} in file {}".format(key, infile) +# try: +# new_value = int(value[0]) +# if new_value >= 1: +# player[key] = new_value +# elif key == 'altsteps' and new_value == 0: +# player[key] = 0 +# except ValueError: +# sys.exit(err) + +# elif key == 'maxstep': # Cannot be less than 0.01 +# err = "Error: expected a float greater than 0.01 for keyword {} in file {}".format(key, infile) +# try: +# new_value = float(value[0]) +# if new_value < 0.01: +# sys.exit(err) +# else: +# player[key] = new_value +# except ValueError: +# sys.exit(err) + +# #### Read the Dice related keywords +# elif key in dice and len(value) != 0: ## 'value' is not empty! + +# if key == 'title': +# dice[key] = value + +# elif key in ('ljname', 'outname', 'progname'): +# dice[key] = value[0] + +# elif key in ('ncores', 'isave'): +# err = "Error: expected a positive integer for keyword {} in file {}".format(key, infile) +# if not value[0].isdigit(): +# sys.exit(err) +# new_value = int(value[0]) +# if new_value >= 1: +# dice[key] = new_value + +# elif key in ('temp', 'press', 'dens'): # Cannot be less than 1e-10 +# err = "Error: expected a positive float for keyword {} in file {}".format(key, infile) +# try: +# new_value = float(value[0]) +# if new_value < 1e-10: +# sys.exit(err) +# else: +# dice[key] = new_value +# except ValueError: +# sys.exit(err) + +# elif key == 'nmol': # If defined, must be well defined (only positive integer values) +# err = "Error: expected 1 to 4 positive integers for keyword {} in file {}".format(key, infile) +# args = min(4, len(value)) +# for i in range(args): +# if value[i].isdigit(): +# new_value = int(value[i]) +# if new_value < 1: +# sys.exit(err) +# else: +# dice[key].append(new_value) +# elif i == 0: +# sys.exit(err) +# else: +# break + +# elif key == 'nstep': # If defined, must be well defined (only positive integer values) +# err = "Error: expected 2 or 3 positive integers for keyword {} in file {}".format(key, infile) +# if len(value) < 2: +# sys.exit(err) +# args = min(3, len(value)) +# for i in range(args): +# if value[i].isdigit(): +# new_value = int(value[i]) +# if new_value < 1: +# sys.exit(err) +# else: +# dice[key].append(new_value) +# elif i < 2: +# sys.exit(err) +# else: +# break + +# #### Read the Gaussian related keywords +# elif key in gaussian and len(value) != 0: ## 'value' is not empty! + +# if key == 'mem': # Memory in MB (minimum of 100) +# err = "Error: expected a positive integer for keyword {} in file {}".format(key, infile) +# if not value[0].isdigit(): +# sys.exit(err) +# new_value = int(value[0]) +# if new_value >= 100: +# gaussian[key] = new_value + +# elif key == 'keywords': +# gaussian[key] = value + +# elif key == 'chgmult': # If defined, must be well defined (2 integer values) +# err = "Error: expected 2 integers for keyword {} in file {}".format(key, infile) +# if len(value) < 2: +# sys.exit(err) +# for i in range (2): +# try: +# gaussian[key][i] = int(value[i]) +# except ValueError: +# sys.exit(err) + +# elif key in ('level', 'chglevel'): +# gaussian[key] = value[0] + +# elif key in ('gmiddle', 'gbottom'): +# gaussian[key] = value[0] + +# elif key == 'pop' and value[0].lower() in ("chelpg", "mk", "nbo"): +# gaussian[key] = value[0].lower() + +# #### Read the Molcas related keywords +# elif key in molcas 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, infile) +# if not value[0].isdigit(): +# sys.exit(err) +# new_value = int(value[0]) +# if new_value >= 1: +# molcas[key] = new_value + +# elif key in ('mbottom', 'orbfile'): +# molcas[key] = value[0] + +# elif key == 'basis': +# molcas[key] = value[0] + +# #### End +# return + + + +# def check_keywords(infile): + +# min_steps = 20000 + +# if dice['ljname'] == None: +# sys.exit("Error: 'ljname' keyword not specified in file {}".format(infile)) + +# if dice['outname'] == None: +# sys.exit("Error: 'outname' keyword not specified in file {}".format(infile)) + +# if dice['dens'] == None: +# sys.exit("Error: 'dens' keyword not specified in file {}".format(infile)) + +# if len(dice['nmol']) == 0: +# sys.exit("Error: 'nmol' keyword not defined appropriately in file {}".format(infile)) + +# if len(dice['nstep']) == 0: +# sys.exit("Error: 'nstep' keyword not defined appropriately in file {}".format(infile)) + +# ## Check only if QM program is Gaussian: +# if player['qmprog'] in ("g03", "g09", "g16"): +# if gaussian['level'] == None: +# sys.exit("Error: 'level' keyword not specified in file {}".format(infile)) + +# if gaussian['gmiddle'] != None: +# if not os.path.isfile(gaussian['gmiddle']): +# sys.exit("Error: file {} not found".format(gaussian['gmiddle'])) + +# if gaussian['gbottom'] != None: +# if not os.path.isfile(gaussian['gbottom']): +# sys.exit("Error: file {} not found".format(gaussian['gbottom'])) + +# if gaussian['pop'] != "chelpg" and (player['ghosts'] == "yes" or player['lps'] == "yes"): +# sys.exit("Error: ghost atoms or lone pairs only available with 'pop = chelpg')") + +# if gaussian['chglevel'] == None: +# gaussian['chglevel'] = gaussian['level'] + +# ## Check only if QM program is Molcas: +# if player['qmprog'] == "molcas": + +# if molcas['mbottom'] == None: +# sys.exit("Error: 'mbottom' keyword not specified in file {}".format(infile)) +# else: +# if not os.path.isfile(molcas['mbottom']): +# sys.exit("Error: file {} not found".format(molcas['mbottom'])) + +# if molcas['basis'] == None: +# sys.exit("Error: 'basis' keyword not specified in file {}".format(infile)) + + +# if player['altsteps'] != 0: + +# ### Verifica se tem mais de 1 molecula QM +# ### (No futuro usar o RMSD fit para poder substituir todas as moleculas QM +# ### no arquivo outname.xy - Need to change the make_init_file!!) +# if dice['nmol'][0] > 1: +# sys.exit("Error: altsteps > 0 only possible with 1 QM molecule (nmol = 1 n2 n3 n4)") + +# # if not zero, altsteps cannot be less than min_steps +# player['altsteps'] = max(min_steps, player['altsteps']) +# # altsteps value is always the nearest multiple of 1000 +# player['altsteps'] = round(player['altsteps'] / 1000) * 1000 + + +# for i in range(len(dice['nstep'])): +# # nstep can never be less than min_steps +# dice['nstep'][i] = max(min_steps, dice['nstep'][i]) +# # nstep values are always the nearest multiple of 1000 +# dice['nstep'][i] = round(dice['nstep'][i] / 1000) * 1000 + +# # isave must be between 100 and 2000 +# dice['isave'] = max(100, dice['isave']) +# dice['isave'] = min(2000, dice['isave']) +# # isave value is always the nearest multiple of 100 +# dice['isave'] = round(dice['isave'] / 100) * 100 + +# return + + + +# def print_keywords(fh): + +# fh.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") +# for var in env: +# fh.write("{} = {}\n".format(var, +# (os.environ[var] if var in os.environ else "Not set"))) + +# fh.write("\n==========================================================================================\n" +# " CONTROL variables being used in this run:\n" +# "------------------------------------------------------------------------------------------\n" +# "\n") + +# for key in sorted(player): +# if player[key] != None: +# if isinstance(player[key], list): +# string = " ".join(str(x) for x in player[key]) +# fh.write("{} = {}\n".format(key, string)) +# else: +# fh.write("{} = {}\n".format(key, player[key])) + +# fh.write("\n") + +# fh.write("------------------------------------------------------------------------------------------\n" +# " DICE variables being used in this run:\n" +# "------------------------------------------------------------------------------------------\n" +# "\n") + +# for key in sorted(dice): +# if dice[key] != None: +# if isinstance(dice[key], list): +# string = " ".join(str(x) for x in dice[key]) +# fh.write("{} = {}\n".format(key, string)) +# else: +# fh.write("{} = {}\n".format(key, dice[key])) + +# fh.write("\n") + +# if player['qmprog'] in ("g03", "g09", "g16"): + +# fh.write("------------------------------------------------------------------------------------------\n" +# " GAUSSIAN variables being used in this run:\n" +# "------------------------------------------------------------------------------------------\n" +# "\n") + +# for key in sorted(gaussian): +# if gaussian[key] != None: +# if isinstance(gaussian[key], list): +# string = " ".join(str(x) for x in gaussian[key]) +# fh.write("{} = {}\n".format(key, string)) +# else: +# fh.write("{} = {}\n".format(key, gaussian[key])) + +# fh.write("\n") + +# elif player['qmprog'] == "molcas": + +# fh.write("------------------------------------------------------------------------------------------\n" +# " MOLCAS variables being used in this run:\n" +# "------------------------------------------------------------------------------------------\n" +# "\n") + +# for key in sorted(molcas): +# if molcas[key] != None: +# if isinstance(molcas[key], list): +# string = " ".join(str(x) for x in molcas[key]) +# fh.write("{} = {}\n".format(key, string)) +# else: +# fh.write("{} = {}\n".format(key, molcas[key])) + +# fh.write("\n") + +# return + + + +# def read_potential(infile): # Deve ser atualizado para o uso de + +# try: +# with open(dice['ljname']) as file: +# ljfile = file.readlines() +# except EnvironmentError as err: +# sys.exit(err) + +# combrule = ljfile.pop(0).split()[0] +# if combrule not in ("*", "+"): +# sys.exit("Error: expected a '*' or a '+' sign in 1st line of file {}".format(dice['ljname'])) +# dice['combrule'] = combrule + +# ntypes = ljfile.pop(0).split()[0] +# if not ntypes.isdigit(): +# sys.exit("Error: expected an integer in the 2nd line of file {}".format(dice['ljname'])) +# ntypes = int(ntypes) + +# if ntypes != len(dice['nmol']): +# sys.exit("Error: number of molecule types in file {} must match that of 'nmol' keyword in file {}".format( +# dice['ljname'], infile)) + +# line = 2 +# for i in range(ntypes): +# line += 1 +# nsites = ljfile.pop(0).split()[0] +# if not nsites.isdigit(): +# sys.exit("Error: expected an integer in line {} of file {}".format(line, dice['ljname'])) + +# nsites = int(nsites) +# molecules.append([]) + +# for j in range(nsites): +# line += 1 +# new_atom = ljfile.pop(0).split() + +# if len(new_atom) < 8: +# sys.exit("Error: expected at least 8 fields in line {} of file {}".format(line, dice['ljname'])) + +# molecules[i].append({}) + +# if not new_atom[0].isdigit(): +# sys.exit("Error: expected an integer in field 1, line {} of file {}".format(line, dice['ljname'])) +# molecules[i][j]['lbl'] = int(new_atom[0]) + +# if not new_atom[1].isdigit(): +# sys.exit("Error: expected an integer in field 2, line {} of file {}".format(line, dice['ljname'])) + +# atnumber = int(new_atom[1]) +# if atnumber == ghost_number and i == 0: # Ghost atom not allowed in the QM molecule +# sys.exit("Error: found a ghost atom in line {} of file {}".format(line, dice['ljname'])) +# molecules[i][j]['na'] = atnumber + +# try: +# molecules[i][j]['rx'] = float(new_atom[2]) +# except: +# sys.exit("Error: expected a float in field 3, line {} of file {}".format(line, dice['ljname'])) + +# try: +# molecules[i][j]['ry'] = float(new_atom[3]) +# except: +# sys.exit("Error: expected a float in field 4, line {} of file {}".format(line, dice['ljname'])) + +# try: +# molecules[i][j]['rz'] = float(new_atom[4]) +# except: +# sys.exit("Error: expected a float in field 5, line {} of file {}".format(line, dice['ljname'])) + +# try: +# molecules[i][j]['chg'] = float(new_atom[5]) +# except: +# sys.exit("Error: expected a float in field 6, line {} of file {}".format(line, dice['ljname'])) + +# try: +# molecules[i][j]['eps'] = float(new_atom[6]) +# except: +# sys.exit("Error: expected a float in field 7, line {} of file {}".format(line, dice['ljname'])) + +# try: +# molecules[i][j]['sig'] = float(new_atom[7]) +# except: +# sys.exit("Error: expected a float in field 8, line {} of file {}".format(line, dice['ljname'])) + +# molecules[i][j]['mass'] = atommass[molecules[i][j]['na']] + +# if len(new_atom) > 8: +# masskey, mass = new_atom[8].partition("=")[::2] +# if masskey.lower() == 'mass' and len(mass) !=0: +# try: +# new_mass = float(mass) +# if new_mass > 0: +# molecules[i][j]['mass'] = new_mass +# except: +# sys.exit( +# "Error: expected a positive float after 'mass=' in field 9, line {} of file {}".format( +# line, dice['ljname'])) + +# return + + + +def read_ghosts(): + + max_atom_number = len(molecules[0]) + + try: + with open("ghosts.in") as fh: + ghostfile = fh.readlines() + except EnvironmentError: + sys.exit("Error: cannot open file ghosts.in") + + for line in ghostfile: + + if len(line.split()) > 1: # Discard lines with less than 2 fields + + key, *atom_numbers = line.split() + key = key.lower() + + if key in ("g", "m", "z"): # Discard lines that do not start with g|m|z + ghost_types.append({}) + ghost_types[-1]['type'] = key + ghost_types[-1]['numbers'] = [] + + for num in atom_numbers: + if not num.isdigit(): + sys.exit("Error: in file ghosts.in: only positive integers allowed after letter g|m|z") + new_num = int(num) + if new_num > max_atom_number: + sys.exit("Error: in file ghosts.in: there is no atom number {}".format(new_num)) + else: + ghost_types[-1]['numbers'].append(new_num) + + if len(ghost_types[-1]['numbers']) < 2: + sys.exit("Error: in file ghosts.in: at least 2 atoms are necessary to make a ghost") + + if len(ghost_types) == 0: + sys.exit("Error: no ghost atom found in ghosts.in") + + return + + + +def read_lps(): + + lp_alpha = 104.0 # Default values + lp_dist = 0.7 # + max_lp_type = 2 + min_alpha = 90.0 + max_alpha = 150.0 + min_dist = 0.5 + max_dist = 1.5 + max_atom_number = len(molecules[0]) + + try: + with open("lps.in") as fh: + lpfile = fh.readlines() + except EnvironmentError: + sys.exit("Error: cannot open file lps.in") + + for line in lpfile: + + if len(line.split()) > 1: # Discard lines with less than 2 fields + + type, *atom_numbers = line.split() + + if type.isdigit(): # Discard lines that do not start with an integer + new_type = int(type) + if new_type > max_lp_type: + sys.exit("Error: in file lps.in: allowed LP types from 1 to {}".format(max_lp_type)) + lp_types.append({}) + lp_types[-1]['type'] = new_type + lp_types[-1]['numbers'] = [] + + # Read types 1 and 2 + if new_type in (1, 2): + + if len(atom_numbers) < 3: + sys.exit("Error: in file lps.in: at least 3 atoms are necessary to make LPs type 1 and 2") + for i in range(3): + num = atom_numbers.pop(0) + if not num.isdigit(): + sys.exit("Error: in file lps.in: expected 3 atom numbers after LPs type 1 and 2") + new_num = int(num) + if new_num > max_atom_number or new_num < 1: + sys.exit("Error: in file lps.in: there is no atom number {}".format(new_num)) + else: + lp_types[-1]['numbers'].append(new_num) + + lp_types[-1]['alpha'] = lp_alpha + lp_types[-1]['dist'] = lp_dist + + if len(atom_numbers) != 0: + try: + alpha = float(atom_numbers.pop(0)) + if alpha > min_alpha and alpha < max_alpha: + lp_types[-1]['alpha'] = alpha + else: + atom_numbers = [] + except: + atom_numbers = [] + + if len(atom_numbers) != 0: + try: + dist = float(atom_numbers.pop(0)) + if dist > min_dist and dist < max_dist: + lp_types[-1]['dist'] = dist + except: + None + # End of types 1 and 2 + + if len(lp_types) == 0: + sys.exit("Error: no lone pair found in lps.in") + + return + + + +# def print_potential(fh): + +# formatstr = "{:<3d} {:>3d} {:>10.5f} {:>10.5f} {:>10.5f} {:>10.6f} {:>9.5f} {:>7.4f} {:>9.4f}\n" +# fh.write("\n" +# "==========================================================================================\n") +# fh.write(" Potential parameters from file {}:\n".format(dice['ljname'])) +# fh.write("------------------------------------------------------------------------------------------\n" +# "\n") + +# fh.write("Combination rule: {}\n".format(dice['combrule'])) +# fh.write("Types of molecules: {}\n\n".format(len(molecules))) + +# i = 0 +# for mol in molecules: +# i += 1 +# fh.write("{} atoms in molecule type {}:\n".format(len(mol), i)) +# fh.write("---------------------------------------------------------------------------------\n" +# "Lbl AN X Y Z Charge Epsilon Sigma Mass\n") +# fh.write("---------------------------------------------------------------------------------\n") + +# for atom in mol: +# fh.write(formatstr.format(atom['lbl'], atom['na'], atom['rx'], atom['ry'], atom['rz'], +# atom['chg'], atom['eps'], atom['sig'], atom['mass'])) + +# fh.write("\n") + +# if player['ghosts'] == "yes" or player['lps'] == "yes": +# fh.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") + +# atoms_string = "" +# for ghost in ghost_types: +# for atom in ghost['numbers']: +# atom_sym = atomsymb[ molecules[0][atom - 1]['na'] ].strip() +# atoms_string += "{}{} ".format(atom_sym,atom) + +# if ghost['type'] == "g": +# fh.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)) +# elif ghost['type'] == "z": +# fh.write(textwrap.fill("* Center of atomic number of atoms {}".format(atoms_string), 80)) + +# fh.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") + +# for lp in lp_types: +# # LP type 1 or 2 +# if lp['type'] in (1, 2): +# atom1_num = lp['numbers'][0] +# atom1_sym = atomsymb[ molecules[0][atom1_num - 1]['na'] ].strip() +# atom2_num = lp['numbers'][1] +# atom2_sym = atomsymb[ molecules[0][atom2_num - 1]['na'] ].strip() +# atom3_num = lp['numbers'][2] +# atom3_sym = atomsymb[ molecules[0][atom3_num - 1]['na'] ].strip() + +# fh.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") + +# # Other LP types + +# fh.write("\n" +# "==========================================================================================\n") + +# return + +## Creation of continue_function + +# def check_executables(fh): + +# fh.write("\n") +# fh.write(90 * "=") +# fh.write("\n\n") + +# dice_path = shutil.which(dice['progname']) +# if dice_path != None: +# fh.write("Program {} found at {}\n".format(dice['progname'], dice_path)) +# else: +# sys.exit("Error: cannot find dice executable") + +# qmprog_path = shutil.which(player['qmprog']) +# if qmprog_path != None: +# fh.write("Program {} found at {}\n".format(player['qmprog'], qmprog_path)) +# else: +# sys.exit("Error: cannot find {} executable".format(player['qmprog'])) + +# if player['qmprog'] in ("g03", "g09", "g16"): +# formchk_path = shutil.which("formchk") +# if formchk_path != None: +# fh.write("Program formchk found at {}\n".format(formchk_path)) +# else: +# sys.exit("Error: cannot find formchk executable") + + +# return + + + diff --git a/DPpack/MolHandling.py b/DPpack/MolHandling.py index fedc65f..d47d817 100644 --- a/DPpack/MolHandling.py +++ b/DPpack/MolHandling.py @@ -480,3 +480,5 @@ class Atom: self.eps = eps # Double self.sig = sig # Double self.mass = atommass[self.na] # Double + + diff --git a/DPpack/SetGlobals.py b/DPpack/SetGlobals.py index a79979d..a05704a 100644 --- a/DPpack/SetGlobals.py +++ b/DPpack/SetGlobals.py @@ -2,812 +2,659 @@ import os, sys import shutil import textwrap +from DPpack.MolHandling import * from DPpack.PTable import * from DPpack.Misc import * -#### Global hashes that control the behaviour of Diceplayer +env = ["OMP_STACKSIZE"] -player = {} -dice = {} -gaussian = {} -molcas = {} -internal = {} - -####################################################################### -#### Global parameters that MAY be given by the user #### -#### (If not given by the user, default values will be used) #### -####################################################################### - -## Diceplayer: -player['maxcyc'] = 1 -player['initcyc'] = 1 # May restart an optimization (append to geoms.xyz from start) -player['nprocs'] = 1 -player['switchcyc'] = 3 # At which step start doing only one QM calculation (geom & chg) -player['altsteps'] = 20000 # Steps for thermalization when starting from previous step -player['maxstep'] = 0.3 # Maxstep for geometry relaxation in Bohr -player['qmprog'] = "g09" -player['opt'] = "yes" -player['freq'] = "no" -player['readhessian'] = "no" -player['lps'] = "no" -player['ghosts'] = "no" -player['vdwforces'] = "no" -player['tol_factor'] = 1.2 # Factor to multiply the default tolerance values - -## Dice: -dice['title'] = "Diceplayer run" -dice['progname'] = "dice" -dice['temp'] = 300.0 -dice['press'] = 1.0 -dice['isave'] = 1000 # ASEC construction will take this into account -dice['ncores'] = 1 - -## Gaussian: -gaussian['mem'] = None -gaussian['keywords'] = None -gaussian['chgmult'] = [0, 1] -gaussian['gmiddle'] = None # In each case, if a filename is given, its content will be -gaussian['gbottom'] = None # inserted in the gaussian input -gaussian['pop'] = "chelpg" -gaussian['chglevel'] = None - -## Molcas: -molcas['orbfile'] = "input.exporb" -molcas['root'] = 1 - - -######################################################################## -#### Global parameters that MUST be given by the user #### -######################################################################## - -## Dice: -dice['dens'] = None # Investigate the possibility of using 'box = Lx Ly Lz' instead. -#dice['box'] = None # So 'geom' would be set by diceplayer and 'cutoff' would be - # switched off. One of them must be given. -dice['ljname'] = None -dice['outname'] = None -dice['nmol'] = [] # Up to 4 integer values related to up to 4 molecule types -dice['nstep'] = [] # 2 or 3 integer values related to 2 or 3 simulations - # (NVT th + NVT eq) or (NVT th + NPT th + NPT eq). - # This will control the 'nstep' keyword of Dice - -## Gaussian: -gaussian['level'] = None - -## Molcas: -molcas['mbottom'] = None -molcas['basis'] = None - - -## The following Dice keywords will be handled automatically by Diceplayer: -## * init ("yes" in the first thermalization and "yesreadxyz" for thermalizations -## starting from a previous step / "no" in subsequent simulations) -## * vstep ("0" for NVT simulations and 'nstep'/5 for NPT simulations) -## * nstep ('nstep' for NVT and "5" for NPT simulations ) -## * irdf ("0" for thermalizations and "10*nprocs" for equilibrium) -## * seed (will be generated randomly each time a Dice input is created) - -## The following Dice keywords will be set constant by Diceplayer for all simulations -## * mstop = 1 (So to guarantee that the ASEC will be correctly built) -## * accum = no (There is never a simulation continuation in Diceplayer) -## * iprint = 1 (Print energy info every step in Dice output) - -## All the other Dice keywords will not be altered from their default values -## and therefore are not mentioned in Diceplayer - - -##################################################################### -#### Global parameters that are not accessible to the user #### -#### (Intended to be used only internally by the program) #### -##################################################################### - -## Diceplayer: -internal['tol_rms_force'] = 3e-4 # Hartree/Bohr -internal['tol_max_force'] = 4.5e-4 # Hartree/Bohr -internal['tol_rms_step'] = 1.2e-3 # Bohr -internal['tol_max_step'] = 1.8e-3 # Bohr -internal['trust_radius'] = None - -## Dice: -internal['combrule'] = None -internal['randominit'] = None - -## Other global variables: -molecules = [] # Armazena todas as informacoes sobre cada tipo de molecula - # (lbl, na, rx, ry, rz, chg, eps, sig, mass) - -internal['ghost_types'] = [] -internal['ghost_atoms'] = [] # Store the ghost atoms (off-atom charge sites) in the QM molecule - # (rx, ry, rz, chg) -internal['lp_types'] = [] -internal['lp_atoms'] = [] # Store the lone pairs (special off-atom charge sites) in the QM molecule - # (rx, ry, rz, chg) - -## Numpy arrays: -step = [] ## Values in Bohr -internal['position'] = [] -internal['energy'] = [] ## Values in Hartree -internal['gradient'] = [] ## Values in Hartree/Bohr -internal['hessian'] = [] ## Values in Hartree/Bohr^2 - -## Conversion factors: bohr2ang = 0.52917721092 ang2bohr = 1/bohr2ang -###################################################################### -#### Environment variables important for executing Diceplayer #### -###################################################################### +class Internal: -env = ["OMP_STACKSIZE"] + def __init__(self, infile): -####################################### functions ###################################### -## Functions to process the input files and store the values in the global variables ## -########################################################################################## + self.infile = infile + self.system = System() -# def read_keywords(infile): + self.player = self.Player() + self.player_keywords = [a for a in dir(self.player) if not a.startswith('__') and not callable(getattr(self.player, a))] -# try: -# with open(infile) as fh: -# controlfile = fh.readlines() -# except EnvironmentError: -# sys.exit("Error: cannot open file {}".format(infile)) - -# for line in controlfile: + self.dice = self.Dice() + self.dice_keywords = [a for a in dir(self.dice) if not a.startswith('__') and not callable(getattr(self.dice, a))] + + self.gaussian = self.Gaussian() + self.gaussian_keywords = [a for a in dir(self.gaussian) if not a.startswith('__') and not callable(getattr(self.gaussian, a))] + + # self.molcas = self.Molcas() + # self.molcas_keywords = [a for a in dir(self.molcas) if not a.startswith('__') and not callable(getattr(self.molcas, a))] + + ## Constanst that shall be set for global use + + self.tol_rms_force = 3e-4 # Hartree/Bohr + self.tol_max_force = 4.5e-4 # Hartree/Bohr + self.tol_rms_step = 1.2e-3 # Bohr + self.tol_max_step = 1.8e-3 # Bohr + self.trust_radius = None + + ## Dice: + self.combrule = None + self.randominit = None + + def read_keywords(self): + + try: + with open(self.infile) as fh: + controlfile = fh.readlines() + except EnvironmentError: + sys.exit("Error: cannot open file {}".format(self.infile)) -# key, value = line.partition("=")[::2] # Discards the '=' -# key = key.strip().lower() -# if key in ('title', 'keywords'): -# value = value.strip() -# else: -# value = value.split() - -# #### Read the Diceplayer related keywords -# if key in player and len(value) != 0: ## 'value' is not empty! + for line in controlfile: -# if key == 'qmprog' and value[0].lower() in ("g03", "g09", "g16", "molcas"): -# player[key] = value[0].lower() + key, value = line.partition("=")[::2] # Discards the '=' + key = key.strip().lower() + if key in ('title', 'keywords'): + value = value.strip() + else: + value = value.split() -# elif key == 'opt' and value[0].lower() in ("yes", "no", "ts"): -# player[key] = value[0].lower() - -# #elif key == 'zipprog' and value[0].lower() in ("zip", "gzip", "bzip"): -# #player[key] = value[0].lower() - -# elif key in ('lps', 'ghosts') and value[0].lower() in ("yes", "no"): -# player[key] = value[0].lower() - -# elif key in ('readhessian', 'vdwforces') and value[0].lower() in ("yes", "no"): -# player[key] = value[0].lower() - -# elif key in ('maxcyc', 'initcyc', 'nprocs', 'altsteps', 'switchcyc'): -# err = "Error: expected a positive integer for keyword {} in file {}".format(key, infile) -# try: -# new_value = int(value[0]) -# if new_value >= 1: -# player[key] = new_value -# elif key == 'altsteps' and new_value == 0: -# player[key] = 0 -# except ValueError: -# sys.exit(err) - -# elif key == 'maxstep': # Cannot be less than 0.01 -# err = "Error: expected a float greater than 0.01 for keyword {} in file {}".format(key, infile) -# try: -# new_value = float(value[0]) -# if new_value < 0.01: -# sys.exit(err) -# else: -# player[key] = new_value -# except ValueError: -# sys.exit(err) - -# #### Read the Dice related keywords -# elif key in dice and len(value) != 0: ## 'value' is not empty! - -# if key == 'title': -# dice[key] = value - -# elif key in ('ljname', 'outname', 'progname'): -# dice[key] = value[0] - -# elif key in ('ncores', 'isave'): -# err = "Error: expected a positive integer for keyword {} in file {}".format(key, infile) -# if not value[0].isdigit(): -# sys.exit(err) -# new_value = int(value[0]) -# if new_value >= 1: -# dice[key] = new_value - -# elif key in ('temp', 'press', 'dens'): # Cannot be less than 1e-10 -# err = "Error: expected a positive float for keyword {} in file {}".format(key, infile) -# try: -# new_value = float(value[0]) -# if new_value < 1e-10: -# sys.exit(err) -# else: -# dice[key] = new_value -# except ValueError: -# sys.exit(err) - -# elif key == 'nmol': # If defined, must be well defined (only positive integer values) -# err = "Error: expected 1 to 4 positive integers for keyword {} in file {}".format(key, infile) -# args = min(4, len(value)) -# for i in range(args): -# if value[i].isdigit(): -# new_value = int(value[i]) -# if new_value < 1: -# sys.exit(err) -# else: -# dice[key].append(new_value) -# elif i == 0: -# sys.exit(err) -# else: -# break - -# elif key == 'nstep': # If defined, must be well defined (only positive integer values) -# err = "Error: expected 2 or 3 positive integers for keyword {} in file {}".format(key, infile) -# if len(value) < 2: -# sys.exit(err) -# args = min(3, len(value)) -# for i in range(args): -# if value[i].isdigit(): -# new_value = int(value[i]) -# if new_value < 1: -# sys.exit(err) -# else: -# dice[key].append(new_value) -# elif i < 2: -# sys.exit(err) -# else: -# break - -# #### Read the Gaussian related keywords -# elif key in gaussian and len(value) != 0: ## 'value' is not empty! - -# if key == 'mem': # Memory in MB (minimum of 100) -# err = "Error: expected a positive integer for keyword {} in file {}".format(key, infile) -# if not value[0].isdigit(): -# sys.exit(err) -# new_value = int(value[0]) -# if new_value >= 100: -# gaussian[key] = new_value - -# elif key == 'keywords': -# gaussian[key] = value - -# elif key == 'chgmult': # If defined, must be well defined (2 integer values) -# err = "Error: expected 2 integers for keyword {} in file {}".format(key, infile) -# if len(value) < 2: -# sys.exit(err) -# for i in range (2): -# try: -# gaussian[key][i] = int(value[i]) -# except ValueError: -# sys.exit(err) - -# elif key in ('level', 'chglevel'): -# gaussian[key] = value[0] - -# elif key in ('gmiddle', 'gbottom'): -# gaussian[key] = value[0] - -# elif key == 'pop' and value[0].lower() in ("chelpg", "mk", "nbo"): -# gaussian[key] = value[0].lower() - -# #### Read the Molcas related keywords -# elif key in molcas 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, infile) -# if not value[0].isdigit(): -# sys.exit(err) -# new_value = int(value[0]) -# if new_value >= 1: -# molcas[key] = new_value - -# elif key in ('mbottom', 'orbfile'): -# molcas[key] = value[0] - -# elif key == 'basis': -# molcas[key] = value[0] - -# #### End -# return - - - -# def check_keywords(infile): - -# min_steps = 20000 - -# if dice['ljname'] == None: -# sys.exit("Error: 'ljname' keyword not specified in file {}".format(infile)) - -# if dice['outname'] == None: -# sys.exit("Error: 'outname' keyword not specified in file {}".format(infile)) - -# if dice['dens'] == None: -# sys.exit("Error: 'dens' keyword not specified in file {}".format(infile)) - -# if len(dice['nmol']) == 0: -# sys.exit("Error: 'nmol' keyword not defined appropriately in file {}".format(infile)) - -# if len(dice['nstep']) == 0: -# sys.exit("Error: 'nstep' keyword not defined appropriately in file {}".format(infile)) - -# ## Check only if QM program is Gaussian: -# if player['qmprog'] in ("g03", "g09", "g16"): -# if gaussian['level'] == None: -# sys.exit("Error: 'level' keyword not specified in file {}".format(infile)) - -# if gaussian['gmiddle'] != None: -# if not os.path.isfile(gaussian['gmiddle']): -# sys.exit("Error: file {} not found".format(gaussian['gmiddle'])) - -# if gaussian['gbottom'] != None: -# if not os.path.isfile(gaussian['gbottom']): -# sys.exit("Error: file {} not found".format(gaussian['gbottom'])) - -# if gaussian['pop'] != "chelpg" and (player['ghosts'] == "yes" or player['lps'] == "yes"): -# sys.exit("Error: ghost atoms or lone pairs only available with 'pop = chelpg')") - -# if gaussian['chglevel'] == None: -# gaussian['chglevel'] = gaussian['level'] - -# ## Check only if QM program is Molcas: -# if player['qmprog'] == "molcas": - -# if molcas['mbottom'] == None: -# sys.exit("Error: 'mbottom' keyword not specified in file {}".format(infile)) -# else: -# if not os.path.isfile(molcas['mbottom']): -# sys.exit("Error: file {} not found".format(molcas['mbottom'])) - -# if molcas['basis'] == None: -# sys.exit("Error: 'basis' keyword not specified in file {}".format(infile)) - - -# if player['altsteps'] != 0: - -# ### Verifica se tem mais de 1 molecula QM -# ### (No futuro usar o RMSD fit para poder substituir todas as moleculas QM -# ### no arquivo outname.xy - Need to change the make_init_file!!) -# if dice['nmol'][0] > 1: -# sys.exit("Error: altsteps > 0 only possible with 1 QM molecule (nmol = 1 n2 n3 n4)") - -# # if not zero, altsteps cannot be less than min_steps -# player['altsteps'] = max(min_steps, player['altsteps']) -# # altsteps value is always the nearest multiple of 1000 -# player['altsteps'] = round(player['altsteps'] / 1000) * 1000 - - -# for i in range(len(dice['nstep'])): -# # nstep can never be less than min_steps -# dice['nstep'][i] = max(min_steps, dice['nstep'][i]) -# # nstep values are always the nearest multiple of 1000 -# dice['nstep'][i] = round(dice['nstep'][i] / 1000) * 1000 - -# # isave must be between 100 and 2000 -# dice['isave'] = max(100, dice['isave']) -# dice['isave'] = min(2000, dice['isave']) -# # isave value is always the nearest multiple of 100 -# dice['isave'] = round(dice['isave'] / 100) * 100 - -# return - - - -# def print_keywords(fh): - -# fh.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") -# for var in env: -# fh.write("{} = {}\n".format(var, -# (os.environ[var] if var in os.environ else "Not set"))) - -# fh.write("\n==========================================================================================\n" -# " CONTROL variables being used in this run:\n" -# "------------------------------------------------------------------------------------------\n" -# "\n") - -# for key in sorted(player): -# if player[key] != None: -# if isinstance(player[key], list): -# string = " ".join(str(x) for x in player[key]) -# fh.write("{} = {}\n".format(key, string)) -# else: -# fh.write("{} = {}\n".format(key, player[key])) - -# fh.write("\n") - -# fh.write("------------------------------------------------------------------------------------------\n" -# " DICE variables being used in this run:\n" -# "------------------------------------------------------------------------------------------\n" -# "\n") - -# for key in sorted(dice): -# if dice[key] != None: -# if isinstance(dice[key], list): -# string = " ".join(str(x) for x in dice[key]) -# fh.write("{} = {}\n".format(key, string)) -# else: -# fh.write("{} = {}\n".format(key, dice[key])) - -# fh.write("\n") - -# if player['qmprog'] in ("g03", "g09", "g16"): - -# fh.write("------------------------------------------------------------------------------------------\n" -# " GAUSSIAN variables being used in this run:\n" -# "------------------------------------------------------------------------------------------\n" -# "\n") - -# for key in sorted(gaussian): -# if gaussian[key] != None: -# if isinstance(gaussian[key], list): -# string = " ".join(str(x) for x in gaussian[key]) -# fh.write("{} = {}\n".format(key, string)) -# else: -# fh.write("{} = {}\n".format(key, gaussian[key])) - -# fh.write("\n") - -# elif player['qmprog'] == "molcas": - -# fh.write("------------------------------------------------------------------------------------------\n" -# " MOLCAS variables being used in this run:\n" -# "------------------------------------------------------------------------------------------\n" -# "\n") - -# for key in sorted(molcas): -# if molcas[key] != None: -# if isinstance(molcas[key], list): -# string = " ".join(str(x) for x in molcas[key]) -# fh.write("{} = {}\n".format(key, string)) -# else: -# fh.write("{} = {}\n".format(key, molcas[key])) - -# fh.write("\n") - -# return - - - -# def read_potential(infile): # Deve ser atualizado para o uso de - -# try: -# with open(dice['ljname']) as file: -# ljfile = file.readlines() -# except EnvironmentError as err: -# sys.exit(err) - -# combrule = ljfile.pop(0).split()[0] -# if combrule not in ("*", "+"): -# sys.exit("Error: expected a '*' or a '+' sign in 1st line of file {}".format(dice['ljname'])) -# dice['combrule'] = combrule - -# ntypes = ljfile.pop(0).split()[0] -# if not ntypes.isdigit(): -# sys.exit("Error: expected an integer in the 2nd line of file {}".format(dice['ljname'])) -# ntypes = int(ntypes) - -# if ntypes != len(dice['nmol']): -# sys.exit("Error: number of molecule types in file {} must match that of 'nmol' keyword in file {}".format( -# dice['ljname'], infile)) - -# line = 2 -# for i in range(ntypes): -# line += 1 -# nsites = ljfile.pop(0).split()[0] -# if not nsites.isdigit(): -# sys.exit("Error: expected an integer in line {} of file {}".format(line, dice['ljname'])) - -# nsites = int(nsites) -# molecules.append([]) - -# for j in range(nsites): -# line += 1 -# new_atom = ljfile.pop(0).split() - -# if len(new_atom) < 8: -# sys.exit("Error: expected at least 8 fields in line {} of file {}".format(line, dice['ljname'])) - -# molecules[i].append({}) - -# if not new_atom[0].isdigit(): -# sys.exit("Error: expected an integer in field 1, line {} of file {}".format(line, dice['ljname'])) -# molecules[i][j]['lbl'] = int(new_atom[0]) - -# if not new_atom[1].isdigit(): -# sys.exit("Error: expected an integer in field 2, line {} of file {}".format(line, dice['ljname'])) - -# atnumber = int(new_atom[1]) -# if atnumber == ghost_number and i == 0: # Ghost atom not allowed in the QM molecule -# sys.exit("Error: found a ghost atom in line {} of file {}".format(line, dice['ljname'])) -# molecules[i][j]['na'] = atnumber - -# try: -# molecules[i][j]['rx'] = float(new_atom[2]) -# except: -# sys.exit("Error: expected a float in field 3, line {} of file {}".format(line, dice['ljname'])) - -# try: -# molecules[i][j]['ry'] = float(new_atom[3]) -# except: -# sys.exit("Error: expected a float in field 4, line {} of file {}".format(line, dice['ljname'])) - -# try: -# molecules[i][j]['rz'] = float(new_atom[4]) -# except: -# sys.exit("Error: expected a float in field 5, line {} of file {}".format(line, dice['ljname'])) - -# try: -# molecules[i][j]['chg'] = float(new_atom[5]) -# except: -# sys.exit("Error: expected a float in field 6, line {} of file {}".format(line, dice['ljname'])) - -# try: -# molecules[i][j]['eps'] = float(new_atom[6]) -# except: -# sys.exit("Error: expected a float in field 7, line {} of file {}".format(line, dice['ljname'])) - -# try: -# molecules[i][j]['sig'] = float(new_atom[7]) -# except: -# sys.exit("Error: expected a float in field 8, line {} of file {}".format(line, dice['ljname'])) - -# molecules[i][j]['mass'] = atommass[molecules[i][j]['na']] - -# if len(new_atom) > 8: -# masskey, mass = new_atom[8].partition("=")[::2] -# if masskey.lower() == 'mass' and len(mass) !=0: -# try: -# new_mass = float(mass) -# if new_mass > 0: -# molecules[i][j]['mass'] = new_mass -# except: -# sys.exit( -# "Error: expected a positive float after 'mass=' in field 9, line {} of file {}".format( -# line, dice['ljname'])) - -# return - - - -def read_ghosts(): - - max_atom_number = len(molecules[0]) - - try: - with open("ghosts.in") as fh: - ghostfile = fh.readlines() - except EnvironmentError: - sys.exit("Error: cannot open file ghosts.in") - - for line in ghostfile: - - if len(line.split()) > 1: # Discard lines with less than 2 fields - - key, *atom_numbers = line.split() - key = key.lower() - - if key in ("g", "m", "z"): # Discard lines that do not start with g|m|z - ghost_types.append({}) - ghost_types[-1]['type'] = key - ghost_types[-1]['numbers'] = [] + #### Read the Diceplayer related keywords + if key in self.player_keywords and len(value) != 0: ## 'value' is not empty! - for num in atom_numbers: - if not num.isdigit(): - sys.exit("Error: in file ghosts.in: only positive integers allowed after letter g|m|z") - new_num = int(num) - if new_num > max_atom_number: - sys.exit("Error: in file ghosts.in: there is no atom number {}".format(new_num)) - else: - ghost_types[-1]['numbers'].append(new_num) + if key == 'qmprog' and value[0].lower() in ("g03", "g09", "g16", "molcas"): + setattr(self.player, key, value[0].lower()) - if len(ghost_types[-1]['numbers']) < 2: - sys.exit("Error: in file ghosts.in: at least 2 atoms are necessary to make a ghost") - - if len(ghost_types) == 0: - sys.exit("Error: no ghost atom found in ghosts.in") - - return - - - -def read_lps(): - - lp_alpha = 104.0 # Default values - lp_dist = 0.7 # - max_lp_type = 2 - min_alpha = 90.0 - max_alpha = 150.0 - min_dist = 0.5 - max_dist = 1.5 - max_atom_number = len(molecules[0]) - - try: - with open("lps.in") as fh: - lpfile = fh.readlines() - except EnvironmentError: - sys.exit("Error: cannot open file lps.in") - - for line in lpfile: - - if len(line.split()) > 1: # Discard lines with less than 2 fields - - type, *atom_numbers = line.split() - - if type.isdigit(): # Discard lines that do not start with an integer - new_type = int(type) - if new_type > max_lp_type: - sys.exit("Error: in file lps.in: allowed LP types from 1 to {}".format(max_lp_type)) - lp_types.append({}) - lp_types[-1]['type'] = new_type - lp_types[-1]['numbers'] = [] + elif key == 'opt' and value[0].lower() in ("yes", "no", "ts"): + setattr(self.player, key, value[0].lower()) - # Read types 1 and 2 - if new_type in (1, 2): + #elif key == 'zipprog' and value[0].lower() in ("zip", "gzip", "bzip"): + #player[key] = value[0].lower() - if len(atom_numbers) < 3: - sys.exit("Error: in file lps.in: at least 3 atoms are necessary to make LPs type 1 and 2") - for i in range(3): - num = atom_numbers.pop(0) - if not num.isdigit(): - sys.exit("Error: in file lps.in: expected 3 atom numbers after LPs type 1 and 2") - new_num = int(num) - if new_num > max_atom_number or new_num < 1: - sys.exit("Error: in file lps.in: there is no atom number {}".format(new_num)) + elif key in ('lps', 'ghosts') and value[0].lower() in ("yes", "no"): + setattr(self.player, key, value[0].lower()) + + 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'): + err = "Error: expected a positive integer for keyword {} in file {}".format(key, self.infile) + try: + new_value = int(value[0]) + if new_value >= 1: + setattr(self.player, key, new_value) + elif key == 'altsteps' and new_value == 0: + setattr(self.player, key, 0) + except ValueError: + sys.exit(err) + + elif key == 'maxstep': # Cannot be less than 0.01 + err = "Error: expected a float greater than 0.01 for keyword {} in file {}".format(key, self.infile) + try: + new_value = float(value[0]) + if new_value < 0.01: + sys.exit(err) else: - lp_types[-1]['numbers'].append(new_num) - - lp_types[-1]['alpha'] = lp_alpha - lp_types[-1]['dist'] = lp_dist - - if len(atom_numbers) != 0: - try: - alpha = float(atom_numbers.pop(0)) - if alpha > min_alpha and alpha < max_alpha: - lp_types[-1]['alpha'] = alpha + setattr(self.player, key).append(new_value) + except ValueError: + sys.exit(err) + + #### Read the Dice related keywords + elif key in self.dice_keywords and len(value) != 0: ## 'value' is not empty! + + if key == 'title': + setattr(self.dice, key, value) + + elif key in ('ljname', 'outname', 'progname'): + setattr(self.dice, key, value[0]) + + elif key in ('ncores', 'isave'): + 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.dice, key, new_value) + + elif key in ('temp', 'press', 'dens'): # Cannot be less than 1e-10 + err = "Error: expected a positive float for keyword {} in file {}".format(key, self.infile) + try: + new_value = float(value[0]) + if new_value < 1e-10: + sys.exit(err) + else: + setattr(self.dice, key, new_value) + except ValueError: + sys.exit(err) + + elif key == 'nmol': # If defined, must be well defined (only positive integer values) + err = "Error: expected 1 to 4 positive integers for keyword {} in file {}".format(key, self.infile) + args = min(4, len(value)) + for i in range(args): + if value[i].isdigit(): + new_value = int(value[i]) + if new_value < 1: + sys.exit(err) else: - atom_numbers = [] - except: - atom_numbers = [] - - if len(atom_numbers) != 0: - try: - dist = float(atom_numbers.pop(0)) - if dist > min_dist and dist < max_dist: - lp_types[-1]['dist'] = dist - except: - None - # End of types 1 and 2 - - if len(lp_types) == 0: - sys.exit("Error: no lone pair found in lps.in") - - return - - - -def print_potential(fh): - - formatstr = "{:<3d} {:>3d} {:>10.5f} {:>10.5f} {:>10.5f} {:>10.6f} {:>9.5f} {:>7.4f} {:>9.4f}\n" - fh.write("\n" - "==========================================================================================\n") - fh.write(" Potential parameters from file {}:\n".format(dice['ljname'])) - fh.write("------------------------------------------------------------------------------------------\n" - "\n") - - fh.write("Combination rule: {}\n".format(dice['combrule'])) - fh.write("Types of molecules: {}\n\n".format(len(molecules))) - - i = 0 - for mol in molecules: - i += 1 - fh.write("{} atoms in molecule type {}:\n".format(len(mol), i)) - fh.write("---------------------------------------------------------------------------------\n" - "Lbl AN X Y Z Charge Epsilon Sigma Mass\n") - fh.write("---------------------------------------------------------------------------------\n") - - for atom in mol: - fh.write(formatstr.format(atom['lbl'], atom['na'], atom['rx'], atom['ry'], atom['rz'], - atom['chg'], atom['eps'], atom['sig'], atom['mass'])) - - fh.write("\n") - - if player['ghosts'] == "yes" or player['lps'] == "yes": - fh.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") - - atoms_string = "" - for ghost in ghost_types: - for atom in ghost['numbers']: - atom_sym = atomsymb[ molecules[0][atom - 1]['na'] ].strip() - atoms_string += "{}{} ".format(atom_sym,atom) + setattr(self.dice, key, new_value) + elif i == 0: + sys.exit(err) + else: + break + + elif key == 'nstep': # If defined, must be well defined (only positive integer values) + err = "Error: expected 2 or 3 positive integers for keyword {} in file {}".format(key, self.infile) + if len(value) < 2: + sys.exit(err) + args = min(3, len(value)) + for i in range(args): + if value[i].isdigit(): + new_value = int(value[i]) + if new_value < 1: + sys.exit(err) + else: + setattr(self.dice, key, new_value) + elif i < 2: + sys.exit(err) + else: + break - if ghost['type'] == "g": - fh.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)) - elif ghost['type'] == "z": - fh.write(textwrap.fill("* Center of atomic number of atoms {}".format(atoms_string), 80)) + #### Read the Gaussian related keywords + elif key in self.gaussian_keywords and len(value) != 0: ## 'value' is not empty! + + if key == 'mem': # Memory in MB (minimum of 100) + 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 >= 100: + setattr(self.gaussian, key, new_value) + + elif key == 'keywords': + setattr(self.gaussian, key, value) + + elif key == 'chgmult': # If defined, must be well defined (2 integer values) + err = "Error: expected 2 integers for keyword {} in file {}".format(key, self.infile) + if len(value) < 2: + sys.exit(err) + for i in range (2): + try: + setattr(self.gaussian, key)[i] = int(value[i]) + except ValueError: + sys.exit(err) + + elif key in ('level', 'chglevel'): + setattr(self.gaussian, key, value[0]) + + elif key in ('gmiddle', 'gbottom'): + setattr(self.gaussian, key, value[0]) + + 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! + + 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 == 'basis': + setattr(self.molcas ,key, value[0]) + + #### End + + + def check_keywords(self): + + min_steps = 20000 + + if self.dice.ljname == None: + sys.exit("Error: 'ljname' keyword not specified in file {}".format(self.infile)) + + if self.dice.outname == None: + sys.exit("Error: 'outname' keyword not specified in file {}".format(self.infile)) + + if self.dice.dens == None: + sys.exit("Error: 'dens' keyword not specified in file {}".format(self.infile)) + + if len(self.dice.nmol) == 0: + sys.exit("Error: 'nmol' keyword not defined appropriately in file {}".format(self.infile)) + + if len(self.dice.nstep) == 0: + sys.exit("Error: 'nstep' keyword not defined appropriately in file {}".format(self.infile)) + + ## Check only if QM program is Gaussian: + if self.player.qmprog in ("g03", "g09", "g16"): + + if self.gaussian.level == None: + sys.exit("Error: 'level' keyword not specified in file {}".format(self.infile)) + + if self.gaussian.gmiddle != None: + if not os.path.isfile(self.gaussian.gmiddle): + sys.exit("Error: file {} not found".format(self.gaussian.gmiddle)) + + if self.gaussian.gbottom != None: + if not os.path.isfile(self.gaussian.gbottom): + sys.exit("Error: file {} not found".format(self.gaussian.gbottom)) + + if self.gaussian.pop != "chelpg" and (self.player.ghosts == "yes" or self.player.lps == "yes"): + sys.exit("Error: ghost atoms or lone pairs only available with 'pop = chelpg')") + + if self.gaussian.chglevel == None: + self.gaussian.chglevel = self.gaussian.level + + ## Check only if QM program is Molcas: + # if self.player.qmprog == "molcas": + + # if self.molcas.mbottom == None: + # sys.exit("Error: 'mbottom' keyword not specified in file {}".format(self.infile)) + # else: + # if not os.path.isfile(self.molcas.mbottom): + # sys.exit("Error: file {} not found".format(self.molcas.mbottom)) + + # if self.molcas.basis == None: + # sys.exit("Error: 'basis' keyword not specified in file {}".format(self.infile)) + + + if self.player.altsteps != 0: + + ### Verifica se tem mais de 1 molecula QM + ### (No futuro usar o RMSD fit para poder substituir todas as moleculas QM + ### no arquivo outname.xy - Need to change the make_init_file!!) + if self.dice.nmol[0] > 1: + sys.exit("Error: altsteps > 0 only possible with 1 QM molecule (nmol = 1 n2 n3 n4)") + + # if not zero, altsteps cannot be less than min_steps + self.player.altsteps = max(min_steps, self.player.altsteps) + # altsteps value is always the nearest multiple of 1000 + self.player.altsteps = round(self.player.altsteps / 1000) * 1000 + + + for i in range(len(self.dice.nstep)): + # nstep can never be less than min_steps + self.dice.nstep[i] = max(min_steps, self.dice.nstep[i]) + # nstep values are always the nearest multiple of 1000 + self.dice.nstep[i] = round(self.dice.nstep[i] / 1000) * 1000 + + # isave must be between 100 and 2000 + self.dice.isave = max(100, self.dice.isave) + self.dice.isave = min(2000, self.dice.isave) + # isave value is always the nearest multiple of 100 + self.dice.isave = round(self.dice.isave / 100) * 100 + + def print_keywords(self, fh): + + fh.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") + for var in env: + fh.write("{} = {}\n".format(var, + (os.environ[var] if var in os.environ else "Not set"))) + + fh.write("\n==========================================================================================\n" + " CONTROL variables being used in this run:\n" + "------------------------------------------------------------------------------------------\n" + "\n") + + for key in sorted(self.player_keywords): + 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)) + else: + fh.write("{} = {}\n".format(key, getattr(self.player,key))) + + fh.write("\n") + + fh.write("------------------------------------------------------------------------------------------\n" + " DICE variables being used in this run:\n" + "------------------------------------------------------------------------------------------\n" + "\n") + + for key in sorted(self.dice_keywords): + 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)) + else: + fh.write("{} = {}\n".format(key, getattr(self.dice,key))) + + fh.write("\n") + + if self.player.qmprog in ("g03", "g09", "g16"): + + fh.write("------------------------------------------------------------------------------------------\n" + " GAUSSIAN variables being used in this run:\n" + "------------------------------------------------------------------------------------------\n" + "\n") + + for key in sorted(self.gaussian_keywords): + 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)) + else: + fh.write("{} = {}\n".format(key, getattr(self.gaussian,key))) fh.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") - - for lp in lp_types: - # LP type 1 or 2 - if lp['type'] in (1, 2): - atom1_num = lp['numbers'][0] - atom1_sym = atomsymb[ molecules[0][atom1_num - 1]['na'] ].strip() - atom2_num = lp['numbers'][1] - atom2_sym = atomsymb[ molecules[0][atom2_num - 1]['na'] ].strip() - atom3_num = lp['numbers'][2] - atom3_sym = atomsymb[ molecules[0][atom3_num - 1]['na'] ].strip() - - fh.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") - - # Other LP types + # elif self.player.qmprog == "molcas": + + # fh.write("------------------------------------------------------------------------------------------\n" + # " MOLCAS variables being used in this run:\n" + # "------------------------------------------------------------------------------------------\n" + # "\n") - fh.write("\n" - "==========================================================================================\n") - - return + # for key in sorted(molcas): + # if molcas[key] != None: + # if isinstance(molcas[key], list): + # string = " ".join(str(x) for x in molcas[key]) + # fh.write("{} = {}\n".format(key, string)) + # else: + # fh.write("{} = {}\n".format(key, molcas[key])) + + # fh.write("\n") -## Creation of continue_function + def read_potential(self): # Deve ser atualizado para o uso de + + try: + with open(self.dice.ljname) as file: + ljfile = file.readlines() + except EnvironmentError as err: + sys.exit(err) + + combrule = ljfile.pop(0).split()[0] + if combrule not in ("*", "+"): + sys.exit("Error: expected a '*' or a '+' sign in 1st line of file {}".format(self.dice.ljname)) + self.dice.combrule = combrule + + ntypes = ljfile.pop(0).split()[0] + if not ntypes.isdigit(): + sys.exit("Error: expected an integer in the 2nd line of file {}".format(self.dice.ljname)) + ntypes = int(ntypes) + + 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): -def check_executables(fh): - - fh.write("\n") - fh.write(90 * "=") - fh.write("\n\n") - - dice_path = shutil.which(dice['progname']) - if dice_path != None: - fh.write("Program {} found at {}\n".format(dice['progname'], dice_path)) - else: - sys.exit("Error: cannot find dice executable") + line += 1 + nsites = ljfile.pop(0).split()[0] + if not nsites.isdigit(): + sys.exit("Error: expected an integer in line {} of file {}".format(line, self.dice.ljname)) + + nsites = int(nsites) + self.system.add_type(Molecule()) + + for j in range(nsites): + + line += 1 + new_atom = ljfile.pop(0).split() + + if len(new_atom) < 8: + sys.exit("Error: expected at least 8 fields in line {} of file {}".format(line, 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, dice['ljname'])) + lbl = int(new_atom[0]) + + if not new_atom[1].isdigit(): + sys.exit("Error: expected an integer in field 2, line {} of file {}".format(line, dice['ljname'])) + + atnumber = int(new_atom[1]) + if atnumber == ghost_number and i == 0: # Ghost atom not allowed in the QM molecule + sys.exit("Error: found a ghost atom in line {} of file {}".format(line, dice['ljname'])) + na = atnumber + + try: + rx = float(new_atom[2]) + except: + sys.exit("Error: expected a float in field 3, line {} of file {}".format(line, dice['ljname'])) + + try: + ry = float(new_atom[3]) + except: + sys.exit("Error: expected a float in field 4, line {} of file {}".format(line, dice['ljname'])) + + try: + rz = float(new_atom[4]) + except: + sys.exit("Error: expected a float in field 5, line {} of file {}".format(line, dice['ljname'])) + + try: + chg = float(new_atom[5]) + except: + sys.exit("Error: expected a float in field 6, line {} of file {}".format(line, dice['ljname'])) + + try: + eps = float(new_atom[6]) + except: + sys.exit("Error: expected a float in field 7, line {} of file {}".format(line, dice['ljname'])) + + try: + sig = float(new_atom[7]) + except: + sys.exit("Error: expected a float in field 8, line {} of file {}".format(line, dice['ljname'])) + + mass = atommass[na] + + if len(new_atom) > 8: + masskey, mass = new_atom[8].partition("=")[::2] + if masskey.lower() == 'mass' and len(mass) !=0: + try: + new_mass = float(mass) + if new_mass > 0: + mass = new_mass + except: + sys.exit( + "Error: expected a positive float after 'mass=' in field 9, line {} of file {}".format( + line, dice['ljname'])) + + self.system.molecule[i].add_atom(Atom(lbl,na,rx,ry,rz,chg,eps,sig,mass)) + + to_delete = ['lbl','na','rx','ry','rz','chg','eps','sig','mass'] + for _var in to_delete: + if _var in locals() or _var in globals(): + exec(f'del {_var}') - qmprog_path = shutil.which(player['qmprog']) - if qmprog_path != None: - fh.write("Program {} found at {}\n".format(player['qmprog'], qmprog_path)) - else: - sys.exit("Error: cannot find {} executable".format(player['qmprog'])) - if player['qmprog'] in ("g03", "g09", "g16"): - formchk_path = shutil.which("formchk") - if formchk_path != None: - fh.write("Program formchk found at {}\n".format(formchk_path)) + def print_potential(self, fh): + + formatstr = "{:<3d} {:>3d} {:>10.5f} {:>10.5f} {:>10.5f} {:>10.6f} {:>9.5f} {:>7.4f} {:>9.4f}\n" + fh.write("\n" + "==========================================================================================\n") + fh.write(" Potential parameters from file {}:\n".format(dice['ljname'])) + fh.write("------------------------------------------------------------------------------------------\n" + "\n") + + fh.write("Combination rule: {}\n".format(self.dice.combrule)) + fh.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" + "Lbl AN X Y Z Charge Epsilon Sigma Mass\n") + fh.write("---------------------------------------------------------------------------------\n") + + for atom in mol.atom: + + fh.write(formatstr.format(atom.lbl, atom.na, atom.rx, atom.ry, atom.rz, + atom.chg, atom.eps, atom.sig, atom.mass)) + + fh.write("\n") + + if self.player.ghosts == "yes" or self.player.lps == "yes": + fh.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") + + # atoms_string = "" + # for ghost in ghost_types: + # for atom in ghost['numbers']: + # atom_sym = atomsymb[ molecules[0][atom - 1]['na'] ].strip() + # atoms_string += "{}{} ".format(atom_sym,atom) + + # if ghost['type'] == "g": + # fh.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)) + # elif ghost['type'] == "z": + # fh.write(textwrap.fill("* Center of atomic number of atoms {}".format(atoms_string), 80)) + + # fh.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") + + # for lp in lp_types: + # # LP type 1 or 2 + # if lp['type'] in (1, 2): + # atom1_num = lp['numbers'][0] + # atom1_sym = atomsymb[ molecules[0][atom1_num - 1]['na'] ].strip() + # atom2_num = lp['numbers'][1] + # atom2_sym = atomsymb[ molecules[0][atom2_num - 1]['na'] ].strip() + # atom3_num = lp['numbers'][2] + # atom3_sym = atomsymb[ molecules[0][atom3_num - 1]['na'] ].strip() + + # fh.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") + + # # Other LP types + + fh.write("\n" + "==========================================================================================\n") + + return + + + def check_executables(self, fh): + + fh.write("\n") + fh.write(90 * "=") + fh.write("\n\n") + + dice_path = shutil.which(dice['progname']) + if dice_path != None: + fh.write("Program {} found at {}\n".format(self.dice.progname, dice_path)) + self.dice.path = dice_path else: - sys.exit("Error: cannot find formchk executable") + sys.exit("Error: cannot find dice executable") - - return + qmprog_path = shutil.which(self.gaussian.qmprog) + if qmprog_path != None: + fh.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)) + + 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)) + else: + sys.exit("Error: cannot find formchk executable") + + class Player: + + def __init__(self): + + self.maxcyc = None + # self.initcyc = 1 # Eliminated + self.nprocs = 1 + self.switchcyc = 3 + self.altsteps = 20000 + self.maxstep = .3 + self.opt = "yes" + self.freq = "no" + self.readhessian = "no" + self.lps = "no" + self.ghosts = "no" + self.vdwforces = "no" + self.tol_factor = 1.2 + + class Dice: + def __init__(self): + + self.title = "Diceplayer run" + self.progname = "dice" + self.path = None + + self.temp = 300.0 + self.press = 1.0 + self.isave = 1000 # ASEC construction will take this into account + self.ncores = 1 + + self.dens = None # Investigate the possibility of using 'box = Lx Ly Lz' instead. + # self.box = None # So 'geom' would be set by diceplayer and 'cutoff' would be + # switched off. One of them must be given. + self.combrule = "*" + self.ljname = None + self.outname = None + self.nmol = [] # Up to 4 integer values related to up to 4 molecule types + self.nstep = [] # 2 or 3 integer values related to 2 or 3 simulations + # (NVT th + NVT eq) or (NVT th + NPT th + NPT eq). + # This will control the 'nstep' keyword of Dice + + + class Gaussian: + + def __init__(self): + + self.qmprog = "g09" + self.path = None + + self.mem = None + self.keywords = None + self.chgmult = [0, 1] + self.gmiddle = None # In each case, if a filename is given, its content will be + self.gbottom = None # inserted in the gaussian input + self.pop = "chelpg" + self.chglevel = None + + self.level = None + + # class Molcas: + + # def __init(self): + + # self.orbfile = "input.exporb" + # self.root = 1 + + # self.mbottom = None + # self.basis = None \ No newline at end of file diff --git a/DPpack/SetGlobalsClass.py b/DPpack/SetGlobalsClass.py deleted file mode 100644 index e140c9f..0000000 --- a/DPpack/SetGlobalsClass.py +++ /dev/null @@ -1,543 +0,0 @@ -import os, sys -import shutil -import textwrap - -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): - - self.infile = infile - self.system = System() - - self.player = self.Player() - self.player_keywords = [a for a in dir(self.player) if not a.startswith('__') and not callable(getattr(self.player, a))] - - self.dice = self.Dice() - self.dice_keywords = [a for a in dir(self.dice) if not a.startswith('__') and not callable(getattr(self.dice, a))] - - self.gaussian = self.Gaussian() - self.gaussian_keywords = [a for a in dir(self.gaussian) if not a.startswith('__') and not callable(getattr(self.gaussian, a))] - - # self.molcas = self.Molcas() - # self.molcas_keywords = [a for a in dir(self.molcas) if not a.startswith('__') and not callable(getattr(self.molcas, a))] - - ## Constanst that shall be set for global use - - self.tol_rms_force = 3e-4 # Hartree/Bohr - self.tol_max_force = 4.5e-4 # Hartree/Bohr - self.tol_rms_step = 1.2e-3 # Bohr - self.tol_max_step = 1.8e-3 # Bohr - self.trust_radius = None - - ## Dice: - self.combrule = None - self.randominit = None - - def read_keywords(self): - - try: - with open(self.infile) as fh: - controlfile = fh.readlines() - except EnvironmentError: - sys.exit("Error: cannot open file {}".format(self.infile)) - - for line in controlfile: - - key, value = line.partition("=")[::2] # Discards the '=' - key = key.strip().lower() - if key in ('title', 'keywords'): - value = value.strip() - else: - value = value.split() - - #### Read the Diceplayer related keywords - if key in self.player_keywords and len(value) != 0: ## 'value' is not empty! - - if key == 'qmprog' and value[0].lower() in ("g03", "g09", "g16", "molcas"): - setattr(self.player, key, value[0].lower()) - - elif key == 'opt' and value[0].lower() in ("yes", "no", "ts"): - setattr(self.player, key, value[0].lower()) - - #elif key == 'zipprog' and value[0].lower() in ("zip", "gzip", "bzip"): - #player[key] = value[0].lower() - - elif key in ('lps', 'ghosts') and value[0].lower() in ("yes", "no"): - setattr(self.player, key, value[0].lower()) - - 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'): - err = "Error: expected a positive integer for keyword {} in file {}".format(key, self.infile) - try: - new_value = int(value[0]) - if new_value >= 1: - setattr(self.player, key, new_value) - elif key == 'altsteps' and new_value == 0: - setattr(self.player, key, 0) - except ValueError: - sys.exit(err) - - elif key == 'maxstep': # Cannot be less than 0.01 - err = "Error: expected a float greater than 0.01 for keyword {} in file {}".format(key, self.infile) - try: - new_value = float(value[0]) - if new_value < 0.01: - sys.exit(err) - else: - setattr(self.player, key).append(new_value) - except ValueError: - sys.exit(err) - - #### Read the Dice related keywords - elif key in self.dice_keywords and len(value) != 0: ## 'value' is not empty! - - if key == 'title': - setattr(self.dice, key, value) - - elif key in ('ljname', 'outname', 'progname'): - setattr(self.dice, key, value[0]) - - elif key in ('ncores', 'isave'): - 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.dice, key, new_value) - - elif key in ('temp', 'press', 'dens'): # Cannot be less than 1e-10 - err = "Error: expected a positive float for keyword {} in file {}".format(key, self.infile) - try: - new_value = float(value[0]) - if new_value < 1e-10: - sys.exit(err) - else: - setattr(self.dice, key, new_value) - except ValueError: - sys.exit(err) - - elif key == 'nmol': # If defined, must be well defined (only positive integer values) - err = "Error: expected 1 to 4 positive integers for keyword {} in file {}".format(key, self.infile) - args = min(4, len(value)) - for i in range(args): - if value[i].isdigit(): - new_value = int(value[i]) - if new_value < 1: - sys.exit(err) - else: - setattr(self.dice, key, new_value) - elif i == 0: - sys.exit(err) - else: - break - - elif key == 'nstep': # If defined, must be well defined (only positive integer values) - err = "Error: expected 2 or 3 positive integers for keyword {} in file {}".format(key, self.infile) - if len(value) < 2: - sys.exit(err) - args = min(3, len(value)) - for i in range(args): - if value[i].isdigit(): - new_value = int(value[i]) - if new_value < 1: - sys.exit(err) - else: - setattr(self.dice, key, new_value) - elif i < 2: - sys.exit(err) - else: - break - - #### Read the Gaussian related keywords - elif key in self.gaussian_keywords and len(value) != 0: ## 'value' is not empty! - - if key == 'mem': # Memory in MB (minimum of 100) - 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 >= 100: - setattr(self.gaussian, key, new_value) - - elif key == 'keywords': - setattr(self.gaussian, key, value) - - elif key == 'chgmult': # If defined, must be well defined (2 integer values) - err = "Error: expected 2 integers for keyword {} in file {}".format(key, self.infile) - if len(value) < 2: - sys.exit(err) - for i in range (2): - try: - setattr(self.gaussian, key)[i] = int(value[i]) - except ValueError: - sys.exit(err) - - elif key in ('level', 'chglevel'): - setattr(self.gaussian, key, value[0]) - - elif key in ('gmiddle', 'gbottom'): - setattr(self.gaussian, key, value[0]) - - 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! - - 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 == 'basis': - setattr(self.molcas ,key, value[0]) - - #### End - - - def check_keywords(self): - - min_steps = 20000 - - if self.dice.ljname == None: - sys.exit("Error: 'ljname' keyword not specified in file {}".format(self.infile)) - - if self.dice.outname == None: - sys.exit("Error: 'outname' keyword not specified in file {}".format(self.infile)) - - if self.dice.dens == None: - sys.exit("Error: 'dens' keyword not specified in file {}".format(self.infile)) - - if len(self.dice.nmol) == 0: - sys.exit("Error: 'nmol' keyword not defined appropriately in file {}".format(self.infile)) - - if len(self.dice.nstep) == 0: - sys.exit("Error: 'nstep' keyword not defined appropriately in file {}".format(self.infile)) - - ## Check only if QM program is Gaussian: - if self.player.qmprog in ("g03", "g09", "g16"): - - if self.gaussian.level == None: - sys.exit("Error: 'level' keyword not specified in file {}".format(self.infile)) - - if self.gaussian.gmiddle != None: - if not os.path.isfile(self.gaussian.gmiddle): - sys.exit("Error: file {} not found".format(self.gaussian.gmiddle)) - - if self.gaussian.gbottom != None: - if not os.path.isfile(self.gaussian.gbottom): - sys.exit("Error: file {} not found".format(self.gaussian.gbottom)) - - if self.gaussian.pop != "chelpg" and (self.player.ghosts == "yes" or self.player.lps == "yes"): - sys.exit("Error: ghost atoms or lone pairs only available with 'pop = chelpg')") - - if self.gaussian.chglevel == None: - self.gaussian.chglevel = self.gaussian.level - - ## Check only if QM program is Molcas: - # if self.player.qmprog == "molcas": - - # if self.molcas.mbottom == None: - # sys.exit("Error: 'mbottom' keyword not specified in file {}".format(self.infile)) - # else: - # if not os.path.isfile(self.molcas.mbottom): - # sys.exit("Error: file {} not found".format(self.molcas.mbottom)) - - # if self.molcas.basis == None: - # sys.exit("Error: 'basis' keyword not specified in file {}".format(self.infile)) - - - if self.player.altsteps != 0: - - ### Verifica se tem mais de 1 molecula QM - ### (No futuro usar o RMSD fit para poder substituir todas as moleculas QM - ### no arquivo outname.xy - Need to change the make_init_file!!) - if self.dice.nmol[0] > 1: - sys.exit("Error: altsteps > 0 only possible with 1 QM molecule (nmol = 1 n2 n3 n4)") - - # if not zero, altsteps cannot be less than min_steps - self.player.altsteps = max(min_steps, self.player.altsteps) - # altsteps value is always the nearest multiple of 1000 - self.player.altsteps = round(self.player.altsteps / 1000) * 1000 - - - for i in range(len(self.dice.nstep)): - # nstep can never be less than min_steps - self.dice.nstep[i] = max(min_steps, self.dice.nstep[i]) - # nstep values are always the nearest multiple of 1000 - self.dice.nstep[i] = round(self.dice.nstep[i] / 1000) * 1000 - - # isave must be between 100 and 2000 - self.dice.isave = max(100, self.dice.isave) - self.dice.isave = min(2000, self.dice.isave) - # isave value is always the nearest multiple of 100 - self.dice.isave = round(self.dice.isave / 100) * 100 - - def print_keywords(self, fh): - - fh.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") - for var in env: - fh.write("{} = {}\n".format(var, - (os.environ[var] if var in os.environ else "Not set"))) - - fh.write("\n==========================================================================================\n" - " CONTROL variables being used in this run:\n" - "------------------------------------------------------------------------------------------\n" - "\n") - - for key in sorted(self.player_keywords): - 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)) - else: - fh.write("{} = {}\n".format(key, getattr(self.player,key))) - - fh.write("\n") - - fh.write("------------------------------------------------------------------------------------------\n" - " DICE variables being used in this run:\n" - "------------------------------------------------------------------------------------------\n" - "\n") - - for key in sorted(self.dice_keywords): - 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)) - else: - fh.write("{} = {}\n".format(key, getattr(self.dice,key))) - - fh.write("\n") - - if self.player.qmprog in ("g03", "g09", "g16"): - - fh.write("------------------------------------------------------------------------------------------\n" - " GAUSSIAN variables being used in this run:\n" - "------------------------------------------------------------------------------------------\n" - "\n") - - for key in sorted(self.gaussian_keywords): - 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)) - else: - fh.write("{} = {}\n".format(key, getattr(self.gaussian,key))) - - fh.write("\n") - - # elif self.player.qmprog == "molcas": - - # fh.write("------------------------------------------------------------------------------------------\n" - # " MOLCAS variables being used in this run:\n" - # "------------------------------------------------------------------------------------------\n" - # "\n") - - # for key in sorted(molcas): - # if molcas[key] != None: - # if isinstance(molcas[key], list): - # string = " ".join(str(x) for x in molcas[key]) - # fh.write("{} = {}\n".format(key, string)) - # else: - # fh.write("{} = {}\n".format(key, molcas[key])) - - # fh.write("\n") - - def read_potential(self): # Deve ser atualizado para o uso de - - try: - with open(self.dice.ljname) as file: - ljfile = file.readlines() - except EnvironmentError as err: - sys.exit(err) - - combrule = ljfile.pop(0).split()[0] - if combrule not in ("*", "+"): - sys.exit("Error: expected a '*' or a '+' sign in 1st line of file {}".format(self.dice.ljname)) - self.dice.combrule = combrule - - ntypes = ljfile.pop(0).split()[0] - if not ntypes.isdigit(): - sys.exit("Error: expected an integer in the 2nd line of file {}".format(self.dice.ljname)) - ntypes = int(ntypes) - - 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): - - line += 1 - nsites = ljfile.pop(0).split()[0] - if not nsites.isdigit(): - sys.exit("Error: expected an integer in line {} of file {}".format(line, self.dice.ljname)) - - nsites = int(nsites) - self.system.add_type(Molecule()) - - for j in range(nsites): - - line += 1 - new_atom = ljfile.pop(0).split() - - if len(new_atom) < 8: - sys.exit("Error: expected at least 8 fields in line {} of file {}".format(line, 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, dice['ljname'])) - lbl = int(new_atom[0]) - - if not new_atom[1].isdigit(): - sys.exit("Error: expected an integer in field 2, line {} of file {}".format(line, dice['ljname'])) - - atnumber = int(new_atom[1]) - if atnumber == ghost_number and i == 0: # Ghost atom not allowed in the QM molecule - sys.exit("Error: found a ghost atom in line {} of file {}".format(line, dice['ljname'])) - na = atnumber - - try: - rx = float(new_atom[2]) - except: - sys.exit("Error: expected a float in field 3, line {} of file {}".format(line, dice['ljname'])) - - try: - ry = float(new_atom[3]) - except: - sys.exit("Error: expected a float in field 4, line {} of file {}".format(line, dice['ljname'])) - - try: - rz = float(new_atom[4]) - except: - sys.exit("Error: expected a float in field 5, line {} of file {}".format(line, dice['ljname'])) - - try: - chg = float(new_atom[5]) - except: - sys.exit("Error: expected a float in field 6, line {} of file {}".format(line, dice['ljname'])) - - try: - eps = float(new_atom[6]) - except: - sys.exit("Error: expected a float in field 7, line {} of file {}".format(line, dice['ljname'])) - - try: - sig = float(new_atom[7]) - except: - sys.exit("Error: expected a float in field 8, line {} of file {}".format(line, dice['ljname'])) - - mass = atommass[na] - - if len(new_atom) > 8: - masskey, mass = new_atom[8].partition("=")[::2] - if masskey.lower() == 'mass' and len(mass) !=0: - try: - new_mass = float(mass) - if new_mass > 0: - mass = new_mass - except: - sys.exit( - "Error: expected a positive float after 'mass=' in field 9, line {} of file {}".format( - line, dice['ljname'])) - - self.system.molecule[i].add_atom(Atom(lbl,na,rx,ry,rz,chg,eps,sig,mass)) - - to_delete = ['lbl','na','rx','ry','rz','chg','eps','sig','mass'] - for _var in to_delete: - if _var in locals() or _var in globals(): - exec(f'del {_var}') - - class Player: - - def __init__(self): - - self.maxcyc = None - # self.initcyc = 1 Eliminated - self.nprocs = 1 - self.switchcyc = 3 - self.altsteps = 20000 - self.maxstep = .3 - self.qmprog = "g09" - self.opt = "yes" - self.freq = "no" - self.readhessian = "no" - self.lps = "no" - self.ghosts = "no" - self.vdwforces = "no" - self.tol_factor = 1.2 - - - - class Dice: - - def __init__(self): - - self.title = "Diceplayer run" - self.progname = "dice" - self.temp = 300.0 - self.press = 1.0 - self.isave = 1000 # ASEC construction will take this into account - self.ncores = 1 - - self.dens = None # Investigate the possibility of using 'box = Lx Ly Lz' instead. - # self.box = None # So 'geom' would be set by diceplayer and 'cutoff' would be - # switched off. One of them must be given. - self.ljname = None - self.outname = None - self.nmol = [] # Up to 4 integer values related to up to 4 molecule types - self.nstep = [] # 2 or 3 integer values related to 2 or 3 simulations - # (NVT th + NVT eq) or (NVT th + NPT th + NPT eq). - # This will control the 'nstep' keyword of Dice - - - class Gaussian: - - def __init__(self): - - self.mem = None - self.keywords = None - self.chgmult = [0, 1] - self.gmiddle = None # In each case, if a filename is given, its content will be - self.gbottom = None # inserted in the gaussian input - self.pop = "chelpg" - self.chglevel = None - - self.level = None - - # class Molcas: - - # def __init(self): - - # self.orbfile = "input.exporb" - # self.root = 1 - - # self.mbottom = None - # self.basis = None \ No newline at end of file