diff --git a/control.in b/control.in index eec4b66..150585f 100644 --- a/control.in +++ b/control.in @@ -1,6 +1,6 @@ # diceplayer maxcyc = 3 -opt = NO +opt = YES nprocs = 4 qmprog = g16 lps = no @@ -11,7 +11,7 @@ altsteps = 20000 ncores = 3 nmol = 1 50 dens = 0.75 -nstep = 40000 60000 50000 +nstep = 4000 6000 isave = 1000 ljname = phb.ljc outname = phb diff --git a/diceplayer/DPpack/.Dice.py b/diceplayer/DPpack/.Dice.py new file mode 100644 index 0000000..b2aac92 --- /dev/null +++ b/diceplayer/DPpack/.Dice.py @@ -0,0 +1,527 @@ +import sys, os, time +import subprocess + +from copy import deepcopy + +from numpy import random + +from DPpack.PTable import * +from DPpack.SetGlobals import * +from DPpack.Misc import * + +####################################### functions ###################################### + +def make_inputs(cycle, proc): + + step_dir = "step{:02d}".format(cycle) + proc_dir = "p{:02d}".format(proc) + path = step_dir + os.sep + proc_dir + + num = time.time() ## Take the decimal places 7 to 12 of the + num = (num - int(num)) * 1e6 ## time in seconds as a floating point + num = int((num - int(num)) * 1e6) ## to make an integer in the range 1-1e6 + random.seed( (os.getpid() * num) % (max_seed + 1) ) + + if not dice['randominit']: + xyzfile = dice['outname'] + ".xyz.last-" + "p{:02d}".format(proc) + make_init_file(path, xyzfile) + + if len(dice['nstep']) == 2: ## Means NVT simulation + + make_nvt_ter(path) + make_nvt_eq(path) + + elif len(dice['nstep']) == 3: ## Means NPT simulation + + if dice['randominit']: + make_nvt_ter(path) + else: + dice['dens'] = new_density(proc) + + make_npt_ter(path) + make_npt_eq(path) + + else: + sys.exit("Error: bad number of entries for 'nstep'") + + make_potential(path) + + return + + +def make_nvt_ter(path): + + file = path + os.sep + "NVT.ter" + try: + fh = open(file, "w") + except: + sys.exit("Error: cannot open file {}".format(file)) + + fh.write("title = {} - NVT Thermalization\n".format(dice['title'])) + fh.write("ncores = {}\n".format(dice['ncores'])) + fh.write("ljname = {}\n".format(dice['ljname'])) + fh.write("outname = {}\n".format(dice['outname'])) + + string = " ".join(str(x) for x in dice['nmol']) + fh.write("nmol = {}\n".format(string)) + + fh.write("dens = {}\n".format(dice['dens'])) + fh.write("temp = {}\n".format(dice['temp'])) + + if dice['randominit']: + fh.write("init = yes\n") + fh.write("nstep = {}\n".format(dice['nstep'][0])) + else: + fh.write("init = yesreadxyz\n") + fh.write("nstep = {}\n".format(player['altsteps'])) + + fh.write("vstep = 0\n") + fh.write("mstop = 1\n") + fh.write("accum = no\n") + fh.write("iprint = 1\n") + fh.write("isave = 0\n") + fh.write("irdf = 0\n") + + seed = int(1e6 * random.random()) + fh.write("seed = {}\n".format(seed)) + + fh.close() + + return + + +def make_nvt_eq(path): + + file = path + os.sep + "NVT.eq" + try: + fh = open(file, "w") + except: + sys.exit("Error: cannot open file {}".format(file)) + + fh.write("title = {} - NVT Production\n".format(dice['title'])) + fh.write("ncores = {}\n".format(dice['ncores'])) + fh.write("ljname = {}\n".format(dice['ljname'])) + fh.write("outname = {}\n".format(dice['outname'])) + + string = " ".join(str(x) for x in dice['nmol']) + fh.write("nmol = {}\n".format(string)) + + fh.write("dens = {}\n".format(dice['dens'])) + fh.write("temp = {}\n".format(dice['temp'])) + fh.write("init = no\n") + fh.write("nstep = {}\n".format(dice['nstep'][1])) + fh.write("vstep = 0\n") + fh.write("mstop = 1\n") + fh.write("accum = no\n") + fh.write("iprint = 1\n") + fh.write("isave = {}\n".format(dice['isave'])) + fh.write("irdf = {}\n".format(10 * player['nprocs'])) + + seed = int(1e6 * random.random()) + fh.write("seed = {}\n".format(seed)) + + fh.close() + + return + + +def make_npt_ter(path): + + file = path + os.sep + "NPT.ter" + try: + fh = open(file, "w") + except: + sys.exit("Error: cannot open file {}".format(file)) + + fh.write("title = {} - NPT Thermalization\n".format(dice['title'])) + fh.write("ncores = {}\n".format(dice['ncores'])) + fh.write("ljname = {}\n".format(dice['ljname'])) + fh.write("outname = {}\n".format(dice['outname'])) + + string = " ".join(str(x) for x in dice['nmol']) + fh.write("nmol = {}\n".format(string)) + + fh.write("press = {}\n".format(dice['press'])) + fh.write("temp = {}\n".format(dice['temp'])) + + if dice['randominit']: + fh.write("init = no\n") ## Because there will be a previous NVT simulation + fh.write("vstep = {}\n".format(int(dice['nstep'][1] / 5))) + else: + fh.write("init = yesreadxyz\n") + fh.write("dens = {:<8.4f}\n".format(dice['dens'])) + fh.write("vstep = {}\n".format(int(player['altsteps'] / 5))) + + fh.write("nstep = 5\n") + fh.write("mstop = 1\n") + fh.write("accum = no\n") + fh.write("iprint = 1\n") + fh.write("isave = 0\n") + fh.write("irdf = 0\n") + + seed = int(1e6 * random.random()) + fh.write("seed = {}\n".format(seed)) + + fh.close() + + return + + +def make_npt_eq(path): + + file = path + os.sep + "NPT.eq" + try: + fh = open(file, "w") + except: + sys.exit("Error: cannot open file {}".format(file)) + + fh.write("title = {} - NPT Production\n".format(dice['title'])) + fh.write("ncores = {}\n".format(dice['ncores'])) + fh.write("ljname = {}\n".format(dice['ljname'])) + fh.write("outname = {}\n".format(dice['outname'])) + + string = " ".join(str(x) for x in dice['nmol']) + fh.write("nmol = {}\n".format(string)) + + fh.write("press = {}\n".format(dice['press'])) + fh.write("temp = {}\n".format(dice['temp'])) + + fh.write("nstep = 5\n") + + fh.write("vstep = {}\n".format(int(dice['nstep'][2] / 5))) + fh.write("init = no\n") + fh.write("mstop = 1\n") + fh.write("accum = no\n") + fh.write("iprint = 1\n") + fh.write("isave = {}\n".format(dice['isave'])) + fh.write("irdf = {}\n".format(10 * player['nprocs'])) + + seed = int(1e6 * random.random()) + fh.write("seed = {}\n".format(seed)) + + fh.close() + + return + + +def make_init_file(path, file): + + if not os.path.isfile(file): + sys.exit("Error: cannot find the xyz file {} in main directory".format(file)) + try: + with open(file) as fh: + xyzfile = fh.readlines() + except: + sys.exit("Error: cannot open file {}".format(file)) + + nsites_mm = 0 + for i in range(1, len(dice['nmol'])): + nsites_mm += dice['nmol'][i] * len(molecules[i]) + + nsites_mm *= -1 ## Become an index to count from the end of xyzfile (list) + xyzfile = xyzfile[nsites_mm :] ## Only the MM atoms of the last configuration remains + + file = path + os.sep + dice['outname'] + ".xy" + + try: + fh = open(file, "w") + except: + sys.exit("Error: cannot open file {}".format(file)) + + for atom in molecules[0]: + fh.write("{:>10.6f} {:>10.6f} {:>10.6f}\n".format(atom['rx'], atom['ry'], + atom['rz'])) + + for ghost in ghost_atoms: + fh.write("{:>10.6f} {:>10.6f} {:>10.6f}\n".format(ghost['rx'], ghost['ry'], + ghost['rz'])) + + for lps in lp_atoms: + fh.write("{:>10.6f} {:>10.6f} {:>10.6f}\n".format(lps['rx'], lps['ry'], + lps['rz'])) + + for line in xyzfile: + atom = line.split() + rx = float(atom[1]) + ry = float(atom[2]) + rz = float(atom[3]) + fh.write("{:>10.5f} {:>10.5f} {:>10.5f}\n".format(rx, ry, rz)) + + fh.write("$end") + + fh.close() + + return + + +def make_potential(path): + + fstr = "{:<3d} {:>3d} {:>10.5f} {:>10.5f} {:>10.5f} {:>10.6f} {:>9.5f} {:>7.4f}\n" + + file = path + os.sep + dice['ljname'] + try: + fh = open(file, "w") + except: + sys.exit("Error: cannot open file {}".format(file)) + + fh.write("{}\n".format(dice['combrule'])) + fh.write("{}\n".format(len(dice['nmol']))) + + nsites_qm = len(molecules[0]) + len(ghost_atoms) + len(lp_atoms) + + ## Print the sites of the QM molecule + fh.write("{}\n".format(nsites_qm)) + for atom in molecules[0]: + fh.write(fstr.format(atom['lbl'], atom['na'], atom['rx'], atom['ry'], atom['rz'], + atom['chg'], atom['eps'], atom['sig'])) + + ghost_label = molecules[0][-1]['lbl'] + 1 + for ghost in ghost_atoms: + fh.write(fstr.format(ghost_label, ghost_number, ghost['rx'], ghost['ry'], + ghost['rz'], ghost['chg'], 0, 0)) + + ghost_label += 1 + for lp in lp_atoms: + fh.write(fstr.format(ghost_label, ghost_number, lp['rx'], lp['ry'], lp['rz'], + lp['chg'], 0, 0)) + + ## Print the sites of the other molecules + for mol in molecules[1:]: + fh.write("{}\n".format(len(mol))) + for atom in mol: + fh.write(fstr.format(atom['lbl'], atom['na'], atom['rx'], atom['ry'], + atom['rz'], atom['chg'], atom['eps'], atom['sig'])) + + return + + +def make_proc_dir(cycle, proc): + + step_dir = "step{:02d}".format(cycle) + proc_dir = "p{:02d}".format(proc) + path = step_dir + os.sep + proc_dir + try: + os.makedirs(path) + except: + sys.exit("Error: cannot make directory {}".format(path)) + + return + + + +def run_dice(cycle, proc, fh): + + step_dir = "step{:02d}".format(cycle) + proc_dir = "p{:02d}".format(proc) + path = step_dir + os.sep + proc_dir + working_dir = os.getcwd() + os.chdir(path) + + fh.write("Simulation process {} initiated with pid {}\n".format(proc_dir, os.getpid())) + + if len(dice['nstep']) == 2: ## Means NVT simulation + + ## NVT thermalization + string = "(from " + ("random" if dice['randominit'] else "previous") + " configuration)" + fh.write("p{:02d}> NVT thermalization initiated {} on {}\n".format(proc, string, + date_time())) + + infh = open("NVT.ter") + outfh = open("NVT.ter.out", "w") + + exit_status = subprocess.call(dice['progname'], stdin=infh, stdout=outfh) + infh.close() + outfh.close() + + if os.getppid() == 1: ## Parent process is dead + sys.exit() + + if exit_status != 0: + sys.exit("Dice process p{:02d} did not exit properly".format(proc)) + else: + outfh = open("NVT.ter.out") ## Open again to seek the normal end flag + flag = outfh.readlines()[dice_flag_line].strip() + outfh.close() + if flag != dice_end_flag: + sys.exit("Dice process p{:02d} did not exit properly".format(proc)) + + ## NVT production + fh.write("p{:02d}> NVT production initiated on {}\n".format(proc, date_time())) + + infh = open("NVT.eq") + outfh = open("NVT.eq.out", "w") + + exit_status = subprocess.call(dice['progname'], stdin=infh, stdout=outfh) + infh.close() + outfh.close() + + if os.getppid() == 1: ## Parent process is dead + sys.exit() + + if exit_status != 0: + sys.exit("Dice process p{:02d} did not exit properly".format(proc)) + else: + outfh = open("NVT.eq.out") ## Open again to seek the normal end flag + flag = outfh.readlines()[dice_flag_line].strip() + outfh.close() + if flag != dice_end_flag: + sys.exit("Dice process p{:02d} did not exit properly".format(proc)) + + fh.write("p{:02d}> ----- NVT production finished on {}\n".format(proc, + date_time())) + + elif len(dice['nstep']) == 3: ## Means NPT simulation + + ## NVT thermalization if randominit + if dice['randominit']: + string = "(from random configuration)" + fh.write("p{:02d}> NVT thermalization initiated {} on {}\n".format(proc, + string, date_time())) + infh = open("NVT.ter") + outfh = open("NVT.ter.out", "w") + + exit_status = subprocess.call(dice['progname'], stdin=infh, stdout=outfh) + infh.close() + outfh.close() + + if os.getppid() == 1: ## Parent process is dead + sys.exit() + + if exit_status != 0: + sys.exit("Dice process p{:02d} did not exit properly".format(proc)) + else: + outfh = open("NVT.ter.out") ## Open again to seek the normal end flag + flag = outfh.readlines()[dice_flag_line].strip() + outfh.close() + if flag != dice_end_flag: + sys.exit("Dice process p{:02d} did not exit properly".format(proc)) + + ## NPT thermalization + string = (" (from previous configuration) " if not dice['randominit'] else " ") + fh.write("p{:02d}> NPT thermalization initiated{}on {}\n".format(proc, string, + date_time())) + infh = open("NPT.ter") + outfh = open("NPT.ter.out", "w") + + exit_status = subprocess.call(dice['progname'], stdin=infh, stdout=outfh) + infh.close() + outfh.close() + + if os.getppid() == 1: ## Parent process is dead + sys.exit() + + if exit_status != 0: + sys.exit("Dice process p{:02d} did not exit properly".format(proc)) + else: + outfh = open("NPT.ter.out") ## Open again to seek the normal end flag + flag = outfh.readlines()[dice_flag_line].strip() + outfh.close() + if flag != dice_end_flag: + sys.exit("Dice process p{:02d} did not exit properly".format(proc)) + + ## NPT production + fh.write("p{:02d}> NPT production initiated on {}\n".format(proc, date_time())) + + infh = open("NPT.eq") + outfh = open("NPT.eq.out", "w") + + exit_status = subprocess.call(dice['progname'], stdin=infh, stdout=outfh) + infh.close() + outfh.close() + + if os.getppid() == 1: ## Parent process is dead + sys.exit() + + if exit_status != 0: + sys.exit("Dice process p{:02d} did not exit properly".format(proc)) + else: + outfh = open("NPT.eq.out") ## Open again to seek the normal end flag + flag = outfh.readlines()[dice_flag_line].strip() + outfh.close() + if flag != dice_end_flag: + sys.exit("Dice process p{:02d} did not exit properly".format(proc)) + + fh.write("p{:02d}> ----- NPT production finished on {}\n".format(proc, + date_time())) + + os.chdir(working_dir) + + return + + + +def print_last_config(cycle, proc): + + step_dir = "step{:02d}".format(cycle) + proc_dir = "p{:02d}".format(proc) + path = step_dir + os.sep + proc_dir + file = path + os.sep + dice['outname'] + ".xyz" + if not os.path.isfile(file): + sys.exit("Error: cannot find the xyz file {}".format(file)) + try: + with open(file) as fh: + xyzfile = fh.readlines() + except: + sys.exit("Error: cannot open file {}".format(file)) + + nsites = ( len(molecules[0]) + len(ghost_atoms) + len(lp_atoms) ) * dice['nmol'][0] + for i in range(1, len(dice['nmol'])): + nsites += dice['nmol'][i] * len(molecules[i]) + + nsites += 2 ## To include the comment line and the number of atoms (xyz file format) + + nsites *= -1 ## Become an index to count from the end of xyzfile (list) + xyzfile = xyzfile[nsites :] ## Take the last configuration + + + file = dice['outname'] + ".xyz.last-" + proc_dir + fh = open(file, "w") + for line in xyzfile: + fh.write(line) + + fh.close() + + return + + + +def new_density(proc): + + file = dice['outname'] + ".xyz.last-" + "p{:02d}".format(proc) + if not os.path.isfile(file): + sys.exit("Error: cannot find the xyz file {} in main directory".format(file)) + try: + with open(file) as fh: + xyzfile = fh.readlines() + except: + sys.exit("Error: cannot open file {}".format(file)) + + box = xyzfile[1].split() + volume = float(box[-3]) * float(box[-2]) * float(box[-1]) + + total_mass = 0 + for i in range(len(molecules)): + mol_mass = 0 + for atom in molecules[i]: + mol_mass += atom['mass'] + total_mass += mol_mass * dice['nmol'][i] + + density = (total_mass / volume) * umaAng3_to_gcm3 + + return density + + + +def simulation_process(cycle, proc, logfh): + + try: + make_proc_dir(cycle, proc) + make_inputs(cycle, proc) + run_dice(cycle, proc, logfh) + except Exception as err: + sys.exit(err) + + return + + + diff --git a/DPpack/.MolHandling.py b/diceplayer/DPpack/.MolHandling.py similarity index 100% rename from DPpack/.MolHandling.py rename to diceplayer/DPpack/.MolHandling.py diff --git a/DPpack/.SetGlobals.py b/diceplayer/DPpack/.SetGlobals.py similarity index 100% rename from DPpack/.SetGlobals.py rename to diceplayer/DPpack/.SetGlobals.py diff --git a/DPpack/Gaussian.py b/diceplayer/DPpack/Gaussian.py similarity index 69% rename from DPpack/Gaussian.py rename to diceplayer/DPpack/Gaussian.py index 5871a0e..177b683 100644 --- a/DPpack/Gaussian.py +++ b/diceplayer/DPpack/Gaussian.py @@ -148,7 +148,7 @@ def print_grad_hessian(cycle, cur_gradient, hessian): return - +## Change the name to make_gaussian_input def make_force_input(cycle, asec_charges): path = "step{:02d}".format(cycle) + os.sep + "qm" @@ -241,92 +241,92 @@ def make_force_input(cycle, asec_charges): -def make_charge_input(cycle, asec_charges): +# def make_charge_input(cycle, asec_charges): - path = "step{:02d}".format(cycle) + os.sep + "qm" - file = path + os.sep + "asec2.gjf" - try: - fh = open(file, "w") - except: - sys.exit("Error: cannot open file {}".format(file)) +# path = "step{:02d}".format(cycle) + os.sep + "qm" +# file = path + os.sep + "asec2.gjf" +# try: +# fh = open(file, "w") +# except: +# sys.exit("Error: cannot open file {}".format(file)) - fh.write("%Chk=asec.chk\n") - if gaussian['mem'] != None: - fh.write("%Mem={}MB\n".format(gaussian['mem'])) - fh.write("%Nprocs={}\n".format(player['nprocs'] * dice['ncores'])) +# fh.write("%Chk=asec.chk\n") +# if gaussian['mem'] != None: +# fh.write("%Mem={}MB\n".format(gaussian['mem'])) +# fh.write("%Nprocs={}\n".format(player['nprocs'] * dice['ncores'])) - kword_line = "#P " + gaussian['chglevel'] + " " + gaussian['keywords'] + " Charge NoSymm" +# kword_line = "#P " + gaussian['chglevel'] + " " + gaussian['keywords'] + " Charge NoSymm" - if player['opt'] != "no" or cycle > 1: - kword_line += " Guess=Read" +# if player['opt'] != "no" or cycle > 1: +# kword_line += " Guess=Read" - kword_line += " Pop={} Density=Current\n".format(gaussian['pop']) +# kword_line += " Pop={} Density=Current\n".format(gaussian['pop']) - fh.write(textwrap.fill(kword_line, 90)) - fh.write("\n") +# fh.write(textwrap.fill(kword_line, 90)) +# fh.write("\n") - fh.write("\nCharge calculation - Cycle number {}\n".format(cycle)) - fh.write("\n") - fh.write("{},{}\n".format(gaussian['chgmult'][0], gaussian['chgmult'][1])) +# fh.write("\nCharge calculation - Cycle number {}\n".format(cycle)) +# fh.write("\n") +# fh.write("{},{}\n".format(gaussian['chgmult'][0], gaussian['chgmult'][1])) - for atom in molecules[0]: - symbol = atomsymb[atom['na']] - fh.write("{:<2s} {:>10.5f} {:>10.5f} {:>10.5f}\n".format(symbol, - atom['rx'], atom['ry'], atom['rz'])) +# for atom in molecules[0]: +# symbol = atomsymb[atom['na']] +# fh.write("{:<2s} {:>10.5f} {:>10.5f} {:>10.5f}\n".format(symbol, +# atom['rx'], atom['ry'], atom['rz'])) - if cycle >= player['switchcyc']: - for ghost in ghost_atoms: - fh.write("Bq {:>10.5f} {:>10.5f} {:>10.5f}\n".format( - ghost['rx'], ghost['ry'], ghost['rz'])) +# if cycle >= player['switchcyc']: +# for ghost in ghost_atoms: +# fh.write("Bq {:>10.5f} {:>10.5f} {:>10.5f}\n".format( +# ghost['rx'], ghost['ry'], ghost['rz'])) - for lp in lp_atoms: - fh.write("Bq {:>10.5f} {:>10.5f} {:>10.5f}\n".format( - lp['rx'], lp['ry'], lp['rz'])) +# for lp in lp_atoms: +# fh.write("Bq {:>10.5f} {:>10.5f} {:>10.5f}\n".format( +# lp['rx'], lp['ry'], lp['rz'])) - fh.write("\n") +# fh.write("\n") - ## If gmiddle file was informed, write its contents in asec.gjf - if gaussian['gmiddle'] != None: - if not os.path.isfile(gaussian['gmiddle']): - sys.exit("Error: cannot find file {} in main directory".format( - gaussian['gmiddle'])) - try: - with open(gaussian['gmiddle']) as gmiddlefile: - gmiddle = gmiddlefile.readlines() - except: - sys.exit("Error: cannot open file {}".format(gaussian['gmiddle'])) +# ## If gmiddle file was informed, write its contents in asec.gjf +# if gaussian['gmiddle'] != None: +# if not os.path.isfile(gaussian['gmiddle']): +# sys.exit("Error: cannot find file {} in main directory".format( +# gaussian['gmiddle'])) +# try: +# with open(gaussian['gmiddle']) as gmiddlefile: +# gmiddle = gmiddlefile.readlines() +# except: +# sys.exit("Error: cannot open file {}".format(gaussian['gmiddle'])) - for line in gmiddle: - fh.write(line) +# for line in gmiddle: +# fh.write(line) - fh.write("\n") +# fh.write("\n") - ## Write the ASEC: - for charge in asec_charges: - fh.write("{:>10.5f} {:>10.5f} {:>10.5f} {:>11.8f}\n".format( - charge['rx'], charge['ry'], charge['rz'], charge['chg'])) +# ## Write the ASEC: +# for charge in asec_charges: +# fh.write("{:>10.5f} {:>10.5f} {:>10.5f} {:>11.8f}\n".format( +# charge['rx'], charge['ry'], charge['rz'], charge['chg'])) - fh.write("\n") +# fh.write("\n") - ## If gbottom file was informed, write its contents in asec.gjf - if gaussian['gbottom'] != None: - if not os.path.isfile(gaussian['gbottom']): - sys.exit("Error: cannot find file {} in main directory".format( - gaussian['gbottom'])) - try: - with open(gaussian['gbottom']) as gbottomfile: - gbottom = gbottomfile.readlines() - except: - sys.exit("Error: cannot open file {}".format(gaussian['gbottom'])) +# ## If gbottom file was informed, write its contents in asec.gjf +# if gaussian['gbottom'] != None: +# if not os.path.isfile(gaussian['gbottom']): +# sys.exit("Error: cannot find file {} in main directory".format( +# gaussian['gbottom'])) +# try: +# with open(gaussian['gbottom']) as gbottomfile: +# gbottom = gbottomfile.readlines() +# except: +# sys.exit("Error: cannot open file {}".format(gaussian['gbottom'])) - for line in gbottom: - fh.write(line) +# for line in gbottom: +# fh.write(line) - fh.write("\n") +# fh.write("\n") - fh.close() +# fh.close() - return +# return @@ -370,51 +370,4 @@ def read_charges(file, fh): fh.write("------------------------------------\n") - return - - - -def run_gaussian(cycle, type, fh): - - path = "step{:02d}".format(cycle) + os.sep + "qm" - work_dir = os.getcwd() - os.chdir(path) - - if type == "force": - infile = "asec.gjf" - elif type == "charge": - infile = "asec2.gjf" - - fh.write("\nCalculation of {}s initiated with Gaussian on {}\n".format(type, date_time())) - - exit_status = subprocess.call([player['qmprog'], infile]) - - if exit_status != 0: - sys.exit("Gaussian process did not exit properly") - - fh.write("Calculation of {}s finished on {}\n".format(type, date_time())) - - os.chdir(work_dir) - - return - - - -def run_formchk(cycle, fh): - - path = "step{:02d}".format(cycle) + os.sep + "qm" - work_dir = os.getcwd() - os.chdir(path) - - fh.write("Formatting the checkpoint file... ") - - exit_status = subprocess.call(["formchk", "asec.chk"]) - - fh.write("Done\n") - - os.chdir(work_dir) - - return - - - + return \ No newline at end of file diff --git a/DPpack/Misc.py b/diceplayer/DPpack/Misc.py similarity index 100% rename from DPpack/Misc.py rename to diceplayer/DPpack/Misc.py diff --git a/DPpack/MolHandling.py b/diceplayer/DPpack/MolHandling.py similarity index 98% rename from DPpack/MolHandling.py rename to diceplayer/DPpack/MolHandling.py index 5b4584b..d0f8f50 100644 --- a/DPpack/MolHandling.py +++ b/diceplayer/DPpack/MolHandling.py @@ -8,8 +8,8 @@ from copy import deepcopy import numpy as np from numpy import linalg -from DPpack.Misc import * -from DPpack.PTable import * +from diceplayer.DPpack.Misc import * +from diceplayer.DPpack.PTable import * env = ["OMP_STACKSIZE"] @@ -44,20 +44,6 @@ class System: return distance - def minimum_distance(self, index1, index2): - - distances = [] - for atom1 in self.molecule[index1]: - if atom1.na != ghost_number: - for atom2 in self.molecule[index2]: - if atom2.na != ghost_number: - dx = atom1.rx - atom2.rx - dy = atom1.ry - atom2.ry - dz = atom1.rz - atom2.rz - distances.append(math.sqrt(dx**2 + dy**2 + dz**2)) - - return min(distances) - def rmsd_fit(self, index_p, index_r): projecting_mol = self.molecule[index_p] @@ -434,6 +420,20 @@ class Molecule: fh.write(" Dipole moment = ( {:>9.4f} , {:>9.4f} , {:>9.4f} ) Total = {:>9.4f} Debye\n\n".format( chg_dip[1], chg_dip[2], chg_dip[3], chg_dip[4])) + def minimum_distance(self, molec): + + distances = [] + for atom1 in self.atom: + if atom1.na != ghost_number: + for atom2 in molec.atom: + if atom2.na != ghost_number: + dx = atom1.rx - atom2.rx + dy = atom1.ry - atom2.ry + dz = atom1.rz - atom2.rz + distances.append(math.sqrt(dx**2 + dy**2 + dz**2)) + + return min(distances) + class Atom: def __init__(self, lbl,na,rx,ry,rz,chg,eps,sig): diff --git a/DPpack/Molcas.py b/diceplayer/DPpack/Molcas.py similarity index 100% rename from DPpack/Molcas.py rename to diceplayer/DPpack/Molcas.py diff --git a/DPpack/Optimization.py b/diceplayer/DPpack/Optimization.py similarity index 100% rename from DPpack/Optimization.py rename to diceplayer/DPpack/Optimization.py diff --git a/DPpack/PTable.py b/diceplayer/DPpack/PTable.py similarity index 100% rename from DPpack/PTable.py rename to diceplayer/DPpack/PTable.py diff --git a/DPpack/SetGlobals.py b/diceplayer/DPpack/SetGlobals.py similarity index 80% rename from DPpack/SetGlobals.py rename to diceplayer/DPpack/SetGlobals.py index faf337b..4bd5381 100644 --- a/DPpack/SetGlobals.py +++ b/diceplayer/DPpack/SetGlobals.py @@ -6,9 +6,9 @@ import types from numpy.core.fromnumeric import partition -from DPpack.MolHandling import * -from DPpack.PTable import * -from DPpack.Misc import * +from diceplayer.DPpack.MolHandling import * +from diceplayer.DPpack.PTable import * +from diceplayer.DPpack.Misc import * from numpy import random import subprocess @@ -129,7 +129,7 @@ class Internal: elif key == 'randominit': if value in ('always','first'): - setattr(self.dice,key,value) + 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) @@ -207,7 +207,7 @@ class Internal: sys.exit(err) elif key == 'level': - setattr(self.gaussian, key, value) + setattr(self.gaussian, key, value[0]) elif key in ('gmiddle', 'gbottom'): setattr(self.gaussian, key, value[0]) @@ -664,9 +664,11 @@ class Internal: ### I still have to talk with Herbet about this function def populate_asec_vdw(self, cycle): + + ## Both asec_charges and vdw_meanfield will utilize the Molecule() class and Atoms() with some None elements - asec_charges = [] # (rx, ry, rz, chg) - vdw_meanfield = [] # (rx, ry, rz, eps, sig) + asec_charges = Molecule() # (lbl=None, na=None, rx, ry, rz, chg, eps=None, sig=None) + vdw_meanfield = Molecule() # (lbl=None, na=None, rx, ry, rz, chg=None, eps, sig) if self.dice.nstep[-1] % self.dice.isave == 0: nconfigs = round(self.dice.nstep[-1] / self.dice.isave) @@ -719,24 +721,25 @@ class Internal: for mol in range(nmols): ## Run over molecules of each type - new_molecule = [] + new_molecule = Molecule(self.system.molecule[type].molnale) for site in range(len(self.system.molecule[types].atom)): ## Run over sites of each molecule new_molecule.append({}) line = xyzfile.pop(0).split() - if line[0].title() != atomsymb[molecules[type][site]['na']].strip(): + if line[0].title() != atomsymb[self.system.molecule[type].atom[site].na.strip()]: sys.exit("Error reading file {}".format(file)) - new_molecule[site]['na'] = molecules[type][site]['na'] - new_molecule[site]['rx'] = float(line[1]) - new_molecule[site]['ry'] = float(line[2]) - new_molecule[site]['rz'] = float(line[3]) - new_molecule[site]['chg'] = molecules[type][site]['chg'] - new_molecule[site]['eps'] = molecules[type][site]['eps'] - new_molecule[site]['sig'] = molecules[type][site]['sig'] + new_molecule.add_atom(Atom(self.system.molecule[type].atom[site].lbl, + self.system.molecule[type].atom[site].na, + self.system.molecule[type].atom[site].float(line[1]), + self.system.molecule[type].atom[site].float(line[2]), + self.system.molecule[type].atom[site].float(line[3]), + self.system.molecule[type].atom[site].chg, + self.system.molecule[type].atom[site].eps, + self.system.molecule[type].atom[site].sig)) - dist = minimum_distance(molecules[0], new_molecule) + dist = self.system.molecule[0].minimum_distance(new_molecule) if dist < thickness[-1]: mol_count += 1 for atom in new_molecule: @@ -748,37 +751,37 @@ class Internal: asec_charges[-1]['rz'] = atom['rz'] asec_charges[-1]['chg'] = atom['chg'] / norm_factor - if player['vdwforces'] == "yes": + if self.player.vdwforces == "yes": vdw_meanfield[-1]['rx'] = atom['rx'] vdw_meanfield[-1]['ry'] = atom['ry'] vdw_meanfield[-1]['rz'] = atom['rz'] vdw_meanfield[-1]['eps'] = atom['eps'] vdw_meanfield[-1]['sig'] = atom['sig'] - #### Read lines with ghosts or lps in molecules of type 0 (reference) - #### and, if dist < thickness, appends to asec - if type == 0: - for ghost in ghost_atoms: - line = xyzfile.pop(0).split() - if line[0] != dice_ghost_label: - sys.exit("Error reading file {}".format(file)) - if dist < thickness[-1]: - asec_charges.append({}) - asec_charges[-1]['rx'] = float(line[1]) - asec_charges[-1]['ry'] = float(line[2]) - asec_charges[-1]['rz'] = float(line[3]) - asec_charges[-1]['chg'] = ghost['chg'] / norm_factor + # #### Read lines with ghosts or lps in molecules of type 0 (reference) + # #### and, if dist < thickness, appends to asec + # if type == 0: + # for ghost in ghost_atoms: + # line = xyzfile.pop(0).split() + # if line[0] != dice_ghost_label: + # sys.exit("Error reading file {}".format(file)) + # if dist < thickness[-1]: + # asec_charges.append({}) + # asec_charges[-1]['rx'] = float(line[1]) + # asec_charges[-1]['ry'] = float(line[2]) + # asec_charges[-1]['rz'] = float(line[3]) + # asec_charges[-1]['chg'] = ghost['chg'] / norm_factor - for lp in lp_atoms: - line = xyzfile.pop(0).split() - if line[0] != dice_ghost_label: - sys.exit("Error reading file {}".format(file)) - if dist < thickness[-1]: - asec_charges.append({}) - asec_charges[-1]['rx'] = float(line[1]) - asec_charges[-1]['ry'] = float(line[2]) - asec_charges[-1]['rz'] = float(line[3]) - asec_charges[-1]['chg'] = lp['chg'] / norm_factor + # for lp in lp_atoms: + # line = xyzfile.pop(0).split() + # if line[0] != dice_ghost_label: + # sys.exit("Error reading file {}".format(file)) + # if dist < thickness[-1]: + # asec_charges.append({}) + # asec_charges[-1]['rx'] = float(line[1]) + # asec_charges[-1]['ry'] = float(line[2]) + # asec_charges[-1]['rz'] = float(line[3]) + # asec_charges[-1]['chg'] = lp['chg'] / norm_factor picked_mols.append(mol_count) @@ -865,12 +868,12 @@ class Internal: try: self.dice.make_proc_dir(cycle, proc) - self.make_inputs(cycle, proc) + self.make_dice_inputs(cycle, proc) self.dice.run_dice(cycle, proc, self.outfile) except Exception as err: sys.exit(err) - def make_inputs(self, cycle, proc): + def make_dice_inputs(self, cycle, proc): sim_dir = "simfiles" step_dir = "step{:02d}".format(cycle) @@ -978,7 +981,7 @@ class Internal: fh.write("mstop = 1\n") fh.write("accum = no\n") fh.write("iprint = 1\n") - fh.write("isave = {}\n".format(self.isave)) + fh.write("isave = {}\n".format(self.dice.isave)) fh.write("irdf = {}\n".format(10 * self.player.nprocs)) seed = int(1e6 * random.random()) @@ -1144,6 +1147,286 @@ class Internal: fh.write(fstr.format(atom.lbl, atom.na, atom.rx, atom.ry, atom.rz, atom.chg, atom.eps, atom.sig)) + # Gaussian related methods + + def read_forces_fchk(self, file, fh): + + forces = [] + try: + with open(file) as tmpfh: + fchkfile = tmpfh.readlines() + except: + sys.exit("Error: cannot open file {}".format(file)) + + start = fchkfile.pop(0).strip() + while start.find("Cartesian Gradient") != 0: ## expression in begining of line + start = fchkfile.pop(0).strip() + + degrees = 3 * len(self.system.molecule[0]) + count = 0 + while True: + values = fchkfile.pop(0).split() + forces.extend([ float(x) for x in values ]) + count += len(values) + if count >= degrees: + forces = forces[:degrees] + break + + gradient = np.array(forces) + + fh.write("\nGradient read from file {}:\n".format(file)) + fh.write("-----------------------------------------------------------------------\n" + "Center Atomic Forces (Hartree/Bohr)\n" + "Number Number X Y Z\n" + "-----------------------------------------------------------------------\n") + for i in range(len(self.system.molecule[0])): + fh.write(" {:>5d} {:>3d} {:>14.9f} {:>14.9f} {:>14.9f}\n".format( + i + 1, self.system.molecule[0][i]['na'], forces.pop(0), forces.pop(0), forces.pop(0))) + + fh.write("-----------------------------------------------------------------------\n") + + force_max = np.amax(np.absolute(gradient)) + force_rms = np.sqrt(np.mean(np.square(gradient))) + + fh.write(" Max Force = {:>14.9f} RMS Force = {:>14.9f}\n\n".format( + force_max, force_rms)) + + return gradient + + + + def read_hessian_fchk(self, file): + + force_const = [] + try: + with open(file) as tmpfh: + fchkfile = tmpfh.readlines() + except: + sys.exit("Error: cannot open file {}".format(file)) + + start = fchkfile.pop(0).strip() + while start.find("Cartesian Force Constants") != 0: + start = fchkfile.pop(0).strip() + + degrees = 3 * len(self.system.molecule[0]) + last = round(degrees * (degrees + 1) / 2) + count = 0 + while True: + values = fchkfile.pop(0).split() + force_const.extend([ float(x) for x in values ]) + count += len(values) + if count >= last: + force_const = force_const[:last] + break + + hessian = np.zeros((degrees, degrees)) + for i in range(degrees): + for j in range(i + 1): + hessian[i,j] = force_const.pop(0) + hessian[j,i] = hessian[i,j] + + return hessian + + + + def read_hessian_log(self, file): + + try: + with open(file) as tmpfh: + logfile = tmpfh.readlines() + except: + sys.exit("Error: cannot open file {}".format(file)) + + start = logfile.pop(0).strip() + while start.find("The second derivative matrix:") != 0: + start = logfile.pop(0).strip() + + degrees = 3 * len(self.system.molecule[0]) + hessian = np.zeros((degrees, degrees)) + + k = 0 + while k < degrees: + logfile.pop(0) + for i in range(k, degrees): + values = logfile.pop(0).split()[1:] + for j in range(k, min(i + 1, k + 5)): + hessian[i,j] = float(values.pop(0)) + hessian[j,i] = hessian[i,j] + k += 5 + + return hessian + + + + def print_grad_hessian(self, cycle, cur_gradient, hessian): + + try: + fh = open("grad_hessian.dat", "w") + except: + sys.exit("Error: cannot open file grad_hessian.dat") + + fh.write("Optimization cycle: {}\n".format(cycle)) + fh.write("Cartesian Gradient\n") + degrees = 3 * len(self.system.molecule[0]) + for i in range(degrees): + fh.write(" {:>11.8g}".format(cur_gradient[i])) + if (i + 1) % 5 == 0 or i == degrees - 1: + fh.write("\n") + + fh.write("Cartesian Force Constants\n") + last = degrees * (degrees + 1) / 2 + count = 0 + for i in range(degrees): + for j in range(i + 1): + count += 1 + fh.write(" {:>11.8g}".format(hessian[i,j])) + if count % 5 == 0 or count == last: + fh.write("\n") + + fh.close() + + return + + + ## Change the name to make_gaussian_input + def make_gaussian_input(self, cycle, asec_charges=None): + + simdir="simfiles" + stepdir="step{:02d}".format(cycle) + path = simdir + os.sep + stepdir + os.sep + "qm" + + file = path + os.sep + "asec.gjf" + + try: + fh = open(file, "w") + except: + sys.exit("Error: cannot open file {}".format(file)) + + fh.write("%Chk=asec.chk\n") + if self.gaussian.mem != None: + fh.write("%Mem={}MB\n".format(self.gaussian.mem)) + fh.write("%Nprocs={}\n".format(self.player.nprocs * self.dice.ncores)) + + kword_line = "#P " + str(self.gaussian.level) + + if self.gaussian.keywords != None: + kword_line += " " + self.gaussian.keywords + + if self.player.opt == 'yes': + kword_line += " Force" + + # kword_line += " Charge" + kword_line += " NoSymm" + kword_line += " Pop={} Density=Current".format(self.gaussian.pop) + + if cycle > 1: + kword_line += " Guess=Read" + + fh.write(textwrap.fill(kword_line, 90)) + fh.write("\n") + + fh.write("\nForce calculation - Cycle number {}\n".format(cycle)) + fh.write("\n") + fh.write("{},{}\n".format(self.gaussian.chgmult[0], self.gaussian.chgmult[1])) + + for atom in self.system.molecule[0].atom: + symbol = atomsymb[atom.na] + fh.write("{:<2s} {:>10.5f} {:>10.5f} {:>10.5f}\n".format(symbol, + atom.rx, atom.ry, atom.rz)) + + # ## If also performing charge fit in the same calculation + # if cycle >= self.player.switchcyc: + # for ghost in ghost_atoms: + # fh.write("Bq {:>10.5f} {:>10.5f} {:>10.5f}\n".format( + # ghost['rx'], ghost['ry'], ghost['rz'])) + + # for lp in lp_atoms: + # fh.write("Bq {:>10.5f} {:>10.5f} {:>10.5f}\n".format( + # lp['rx'], lp['ry'], lp['rz'])) + + # fh.write("\n") + + ## If gmiddle file was informed, write its contents in asec.gjf + # if self.gaussian.gmiddle != None: + # if not os.path.isfile(self.gaussian.gmiddle): + # sys.exit("Error: cannot find file {} in main directory".format( + # self.gaussian.gmiddle)) + # try: + # with open(self.gaussian.gmiddle) as gmiddlefile: + # gmiddle = gmiddlefile.readlines() + # except: + # sys.exit("Error: cannot open file {}".format(self.gaussian.gmiddle)) + + # for line in gmiddle: + # fh.write(line) + + # fh.write("\n") + + # ## Write the ASEC: + # for charge in asec_charges: + # fh.write("{:>10.5f} {:>10.5f} {:>10.5f} {:>11.8f}\n".format( + # charge['rx'], charge['ry'], charge['rz'], charge['chg'])) + + fh.write("\n") + + # ## If gbottom file was informed, write its contents in asec.gjf + # if self.gaussian.gbottom != None: + # if not os.path.isfile(self.gaussian.gbottom): + # sys.exit("Error: cannot find file {} in main directory".format( + # self.gaussian.gbottom)) + # try: + # with open(self.gaussian.gbottom) as gbottomfile: + # gbottom = gbottomfile.readlines() + # except: + # sys.exit("Error: cannot open file {}".format(self.gaussian.gbottom)) + + # for line in gbottom: + # fh.write(line) + + # fh.write("\n") + + # fh.close() + + def read_charges(self, file, fh): + + try: + with open(file) as tmpfh: + glogfile = tmpfh.readlines() + except: + sys.exit("Error: cannot open file {}".format(file)) + + start = glogfile.pop(0).strip() + while start != "Fitting point charges to electrostatic potential": + start = glogfile.pop(0).strip() + + glogfile = glogfile[3:] ## Consume 3 more lines + + fh.write("\nAtomic charges:\n") + fh.write("------------------------------------\n") + for atom in self.system.molecule[0].atom: + line = glogfile.pop(0).split() + atom_str = line[1] + charge = float(line[2]) + atom.chg = charge + fh.write(" {:<2s} {:>10.6f}\n".format(atom_str, charge)) + + # if self.gaussian.pop == "chelpg": + # for ghost in ghost_atoms: + # line = glogfile.pop(0).split() + # atom_str = line[1] + # charge = float(line[2]) + # ghost['chg'] = charge + # fh.write(" {:<2s} {:>10.6f}\n".format(atom_str, charge)) + + # for lp in lp_atoms: + # line = glogfile.pop(0).split() + # atom_str = line[1] + # charge = float(line[2]) + # lp['chg'] = charge + # fh.write(" {:<2s} {:>10.6f}\n".format(atom_str, charge)) + + fh.write("------------------------------------\n") + class Player: def __init__(self): @@ -1389,6 +1672,52 @@ class Internal: self.pop = "chelpg" self.level = None + def run_gaussian(self, cycle, type, fh): + + simdir="simfiles" + stepdir="step{:02d}".format(cycle) + path = simdir + os.sep + stepdir + os.sep + "qm" + work_dir = os.getcwd() + os.chdir(path) + + # if type == "force": + # infile = "asec.gjf" + # elif type == "charge": + # infile = "asec2.gjf" + + infile = "asec.gjf" + + fh.write("\nCalculation of {}s initiated with Gaussian on {}\n".format(type, date_time())) + + if shutil.which("bash") != None: + exit_status = subprocess.call(["bash","-c","exec -a {}-step{} {} {}".format(self.qmprog, cycle, self.qmprog, infile)]) + else: + exit_status = subprocess.call([self.qmprog, infile]) + + if exit_status != 0: + sys.exit("Gaussian process did not exit properly") + + fh.write("Calculation of {}s finished on {}\n".format(type, date_time())) + + os.chdir(work_dir) + + def run_formchk(self, cycle, fh): + + simdir="simfiles" + stepdir="step{:02d}".format(cycle) + path = simdir + os.sep + stepdir + os.sep + "qm" + + work_dir = os.getcwd() + os.chdir(path) + + fh.write("Formatting the checkpoint file... ") + + exit_status = subprocess.call(["formchk", "asec.chk"]) + + fh.write("Done\n") + + os.chdir(work_dir) + # class Molcas: # def __init(self): diff --git a/DPpack/__init__.py b/diceplayer/DPpack/__init__.py similarity index 100% rename from DPpack/__init__.py rename to diceplayer/DPpack/__init__.py diff --git a/diceplayer.py b/diceplayer/__main__.py similarity index 79% rename from diceplayer.py rename to diceplayer/__main__.py index be773ab..363985c 100644 --- a/diceplayer.py +++ b/diceplayer/__main__.py @@ -1,4 +1,4 @@ -#!/export/apps/python/361/bin/python3 +#!/usr/bin/python3 import os, sys, time, signal import setproctitle @@ -6,19 +6,18 @@ import argparse import shutil from multiprocessing import Process, connection -import DPpack.Gaussian as Gaussian -from DPpack.PTable import * -from DPpack.SetGlobals import * -from DPpack.MolHandling import * -from DPpack.Misc import * +from diceplayer.DPpack.PTable import * +from diceplayer.DPpack.SetGlobals import * +from diceplayer.DPpack.MolHandling import * +from diceplayer.DPpack.Misc import * +_version = 'dev' +setproctitle.setproctitle("diceplayer-{}".format(_version)) if __name__ == '__main__': #### Read and store the arguments passed to the program #### #### and set the usage and help messages #### - setproctitle.setproctitle("diceplayer-current") - parser = argparse.ArgumentParser(prog='Diceplayer') parser.add_argument('--continue', dest='opt_continue' , default=False, action='store_true') parser.add_argument('--version', action='version', version='%(prog)s 1.0') @@ -200,11 +199,11 @@ if __name__ == '__main__': internal.outfile.write("\n+" + 88 * "-" + "+\n") - #### - #### End of parallel simulations block - #### + ### + ### End of parallel simulations block + ### - # ## Make ASEC + ## Make ASEC # internal.outfile.write("\nBuilding the ASEC and vdW meanfields... ") # asec_charges = internal.populate_asec_vdw(cycle) @@ -213,31 +212,26 @@ if __name__ == '__main__': # path = "step{:02d}".format(cycle) + os.sep + "p{:02d}".format(proc) # compress_files_1mb(path) - #### - #### Start QM calculation - #### + ### + ### Start QM calculation + ### -# make_qm_dir(cycle) - -# if internal.player.opt == "yes": + make_qm_dir(cycle) + + if internal.player.qmprog in ("g03", "g09", "g16"): -# ## -# ## Gaussian block -# ## -# if internal.player.qmprog in ("g03", "g09", "g16"): - -# if cycle > 1: -# src = "step{:02d}".format(cycle - 1) + os.sep + "qm" + os.sep + "asec.chk" -# dst = "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec.chk" -# shutil.copyfile(src, dst) - -# Gaussian.make_force_input(cycle, asec_charges) -# Gaussian.run_gaussian(cycle, "force", internal.outfile) -# Gaussian.run_formchk(cycle, internal.outfile) + if cycle > 1: + src = "simfiles" + os.sep + "step{:02d}".format(cycle - 1) + os.sep + "qm" + os.sep + "asec.chk" + dst = "simfiles" + os.sep + "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec.chk" + shutil.copyfile(src, dst) + + internal.make_gaussian_input(cycle) + internal.gaussian.run_gaussian(cycle, "force", internal.outfile) +# internal.gaussian.run_formchk(cycle, internal.outfile) # ## Read the gradient -# file = "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec.fchk" -# gradient = Gaussian.read_forces(file, internal.outfile) +# file = "simfiles" + os.sep + "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec.fchk" +# gradient = internal.read_forces(file, internal.outfile) # if len(cur_gradient) > 0: # old_gradient = cur_gradient # cur_gradient = gradient @@ -247,9 +241,9 @@ if __name__ == '__main__': # if internal.player.readhessian == "yes": # file = "grad_hessian.dat" # internal.outfile.write("\nReading the hessian matrix from file {}\n".format(file)) -# hessian = Gaussian.read_hessian_fchk(file) +# hessian = internal.read_hessian_fchk(file) # else: -# file = "step01" + os.sep + "qm" + os.sep + "asec.fchk" +# file = "simfiles" + os.sep + "step01" + os.sep + "qm" + os.sep + "asec.fchk" # internal.outfile.write("\nReading the hessian matrix from file {}\n".format(file)) # hessian = internal.gaussian.read_hessian(file) @@ -272,15 +266,15 @@ if __name__ == '__main__': # ## If needed, calculate the charges # if cycle < internal.player.switchcyc: -# internal.gaussian.make_charge_input(cycle, asec_charges) +# # internal.gaussian.make_charge_input(cycle, asec_charges) # internal.gaussian.run_gaussian(cycle, "charge", internal.outfile) # ## Read the new charges and update molecules[0] # if cycle < internal.player.switchcyc: -# file = "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec2.log" +# file = "simfiles" + os.sep + "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec2.log" # internal.gaussian.read_charges(file, internal.outfile) # else: -# file = "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec.log" +# file = "simfiles" + os.sep + "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec.log" # internal.gaussian.read_charges(file, internal.outfile) # ## Print new info for molecule[0] @@ -316,15 +310,15 @@ if __name__ == '__main__': # if internal.player.qmprog in ("g03", "g09", "g16"): # if cycle > 1: -# src = "step{:02d}".format(cycle - 1) + os.sep + "qm" + os.sep + "asec.chk" -# dst = "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec.chk" +# src = "simfiles" + os.sep + "step{:02d}".format(cycle - 1) + os.sep + "qm" + os.sep + "asec.chk" +# dst = "simfiles" + os.sep + "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec.chk" # shutil.copyfile(src, dst) -# Gaussian.make_charge_input(cycle, asec_charges) -# Gaussian.run_gaussian(cycle, "charge", internal.outfile) +# # internal.gaussian.make_charge_input(cycle, asec_charges) +# internal.gaussian.run_gaussian(cycle, "charge", internal.outfile) -# file = "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec2.log" -# Gaussian.read_charges(file) +# file = "simfiles" + os.sep + "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec2.log" +# internal.read_charges(file) # ## Print new info for molecule[0] # internal.outfile.write("\nNew values for molecule type 1:\n\n") @@ -338,11 +332,11 @@ if __name__ == '__main__': # #### End of the iterative process # #### -# ## imprimir ultimas mensagens, criar um arquivo de potencial para ser usado em eventual -# ## continuacao, fechar arquivos (geoms.xyz, run.log, ...) +# # ## imprimir ultimas mensagens, criar um arquivo de potencial para ser usado em eventual +# # ## continuacao, fechar arquivos (geoms.xyz, run.log, ...) -# internal.outfile.write("\nDiceplayer finished normally!\n") -# internal.outfile.close() -# #### -# #### End of the program -# #### +# # internal.outfile.write("\nDiceplayer finished normally!\n") +# # internal.outfile.close() +# # #### +# # #### End of the program +# # #### diff --git a/run.log b/run.log index aee2304..ade01f2 100644 --- a/run.log +++ b/run.log @@ -5,7 +5,7 @@ Your python version is 3.6.9 (default, Jan 26 2021, 15:33:00) [GCC 8.4.0] -Program started on Saturday, 20 Nov 2021 at 14:52:20 +Program started on Friday, 03 Dec 2021 at 21:20:24 Environment variables: OMP_STACKSIZE = Not set @@ -22,7 +22,7 @@ lps = no maxcyc = 3 maxstep = 0.3 nprocs = 4 -opt = no +opt = yes qmprog = g16 readhessian = no switchcyc = 3 @@ -39,7 +39,7 @@ isave = 1000 ljname = phb.ljc ncores = 3 nmol = 1 50 -nstep = 40000 60000 50000 +nstep = 20000 20000 outname = phb press = 1.0 progname = dice @@ -115,8 +115,8 @@ Lbl AN X Y Z Charge Epsilon Sigma Mass 1 6 2.07203 -2.82345 0.00263 -0.115000 0.07000 3.5500 12.0110 1 6 2.76823 -1.61764 0.00263 -0.115000 0.07000 3.5500 12.0110 1 6 2.06824 -0.40521 0.00264 -0.115000 0.07000 3.5500 12.0110 -1 14 0.67589 -0.40522 0.00264 0.150000 0.07000 3.5500 28.0860 -1 104 -0.02420 -1.61760 0.00263 -0.115000 0.07000 3.5500 0.0000 +1 6 0.67589 -0.40522 0.00264 0.150000 0.07000 3.5500 12.0110 +1 6 -0.02420 -1.61760 0.00263 -0.115000 0.07000 3.5500 12.0110 2 1 0.13203 -3.75875 0.00263 0.115000 0.03000 2.4200 1.0079 2 1 2.61203 -3.75875 0.00263 0.115000 0.03000 2.4200 1.0079 2 1 2.60824 0.53010 0.00264 0.115000 0.03000 2.4200 1.0079 @@ -159,13 +159,13 @@ Molecule type 1 - PHENOL_BLUE: Molecule type 2 - PLACEHOLDER: - Center of mass = ( 1.6067 , -1.0929 , 0.0026 ) - Moments of inertia = 1.205279E+02 2.859254E+02 4.033761E+02 - Major principal axis = ( 0.795913 , -0.605411 , -0.000001 ) - Inter principal axis = ( 0.605411 , 0.795913 , 0.000001 ) + Center of mass = ( 1.5639 , -1.2534 , 0.0026 ) + Moments of inertia = 1.394830E+02 2.777126E+02 4.141185E+02 + Major principal axis = ( 0.865900 , -0.500217 , -0.000001 ) + Inter principal axis = ( 0.500217 , 0.865900 , 0.000001 ) Minor principal axis = ( 0.000000 , -0.000002 , 1.000000 ) Characteristic lengths = ( 5.74 , 5.26 , 1.51 ) - Total mass = 112.20 au + Total mass = 108.14 au Total charge = -0.0000 e Dipole moment = ( 2.3293 , 0.1593 , -0.0000 ) Total = 2.3347 Debye @@ -173,14 +173,14 @@ Molecule type 2 - PLACEHOLDER: New values: Center of mass = ( 0.0000 , -0.0000 , 0.0000 ) - Moments of inertia = 1.205279E+02 2.859254E+02 4.033761E+02 - Major principal axis = ( -1.000000 , 0.000000 , 0.000000 ) + Moments of inertia = 1.394830E+02 2.777126E+02 4.141185E+02 + Major principal axis = ( 1.000000 , -0.000000 , -0.000000 ) Inter principal axis = ( 0.000000 , 1.000000 , 0.000000 ) - Minor principal axis = ( 0.000000 , 0.000000 , 1.000000 ) - Characteristic lengths = ( 5.67 , 5.13 , 1.51 ) - Total mass = 112.20 au + Minor principal axis = ( 0.000000 , -0.000000 , 1.000000 ) + Characteristic lengths = ( 5.65 , 4.95 , 1.51 ) + Total mass = 108.14 au Total charge = -0.0000 e - Dipole moment = ( 1.7575 , 1.5369 , -0.0000 ) Total = 2.3347 Debye + Dipole moment = ( 1.9373 , 1.3031 , -0.0000 ) Total = 2.3347 Debye ========================================================================================== @@ -190,67 +190,23 @@ Starting the iterative process. Step # 1 ------------------------------------------------------------------------------------------ -Simulation process simfiles/step01/p02 initiated with pid 15516 -p02> NVT thermalization initiated (from random configuration) on 20 Nov 2021 at 14:52:20 -Simulation process simfiles/step01/p01 initiated with pid 15515 -p01> NVT thermalization initiated (from random configuration) on 20 Nov 2021 at 14:52:20 -Simulation process simfiles/step01/p03 initiated with pid 15517 -p03> NVT thermalization initiated (from random configuration) on 20 Nov 2021 at 14:52:20 -Simulation process simfiles/step01/p04 initiated with pid 15518 -p04> NVT thermalization initiated (from random configuration) on 20 Nov 2021 at 14:52:20 -p03> NPT thermalization finished (from previous configuration) on 20 Nov 2021 at 14:54:12 -p04> NPT thermalization finished (from previous configuration) on 20 Nov 2021 at 14:54:13 -p02> NPT thermalization finished (from previous configuration) on 20 Nov 2021 at 14:54:13 -p01> NPT thermalization finished (from previous configuration) on 20 Nov 2021 at 14:54:13 -p03> NPT production initiated on 20 Nov 2021 at 14:57:16 -p02> NPT production initiated on 20 Nov 2021 at 14:57:19 -p04> NPT production initiated on 20 Nov 2021 at 14:57:20 -p01> NPT production initiated on 20 Nov 2021 at 14:57:21 -p03> ----- NPT production finished on 20 Nov 2021 at 15:00:16 -p02> ----- NPT production finished on 20 Nov 2021 at 15:00:21 -p04> ----- NPT production finished on 20 Nov 2021 at 15:00:23 -p01> ----- NPT production finished on 20 Nov 2021 at 15:00:27 +Simulation process simfiles/step01/p01 initiated with pid 8024 +p01> NVT thermalization finished (from random configuration) on 03 Dec 2021 at 21:20:24 +Simulation process simfiles/step01/p03 initiated with pid 8026 +p03> NVT thermalization finished (from random configuration) on 03 Dec 2021 at 21:20:24 +Simulation process simfiles/step01/p02 initiated with pid 8025 +p02> NVT thermalization finished (from random configuration) on 03 Dec 2021 at 21:20:24 +Simulation process simfiles/step01/p04 initiated with pid 8027 +p04> NVT thermalization finished (from random configuration) on 03 Dec 2021 at 21:20:24 +p02> NVT production initiated on 03 Dec 2021 at 21:21:22 +p04> NVT production initiated on 03 Dec 2021 at 21:21:23 +p03> NVT production initiated on 03 Dec 2021 at 21:21:23 +p01> NVT production initiated on 03 Dec 2021 at 21:21:27 +p03> ----- NVT production finished on 03 Dec 2021 at 21:22:35 +p04> ----- NVT production finished on 03 Dec 2021 at 21:22:36 +p02> ----- NVT production finished on 03 Dec 2021 at 21:22:36 +p01> ----- NVT production finished on 03 Dec 2021 at 21:22:39 +----------------------------------------------------------------------------------------+ - Step # 2 ------------------------------------------------------------------------------------------- -Simulation process simfiles/step02/p02 initiated with pid 15559 -p02> NPT thermalization finished (from previous configuration) on 20 Nov 2021 at 15:00:27 -Simulation process simfiles/step02/p01 initiated with pid 15558 -p01> NPT thermalization finished (from previous configuration) on 20 Nov 2021 at 15:00:27 -Simulation process simfiles/step02/p03 initiated with pid 15560 -p03> NPT thermalization finished (from previous configuration) on 20 Nov 2021 at 15:00:27 -Simulation process simfiles/step02/p04 initiated with pid 15561 -p04> NPT thermalization finished (from previous configuration) on 20 Nov 2021 at 15:00:27 -p03> NPT production initiated on 20 Nov 2021 at 15:01:30 -p01> NPT production initiated on 20 Nov 2021 at 15:01:30 -p02> NPT production initiated on 20 Nov 2021 at 15:01:31 -p04> NPT production initiated on 20 Nov 2021 at 15:01:31 -p03> ----- NPT production finished on 20 Nov 2021 at 15:04:29 -p02> ----- NPT production finished on 20 Nov 2021 at 15:04:33 -p01> ----- NPT production finished on 20 Nov 2021 at 15:04:34 -p04> ----- NPT production finished on 20 Nov 2021 at 15:04:37 - -+----------------------------------------------------------------------------------------+ - Step # 3 ------------------------------------------------------------------------------------------- - -Simulation process simfiles/step03/p01 initiated with pid 15661 -p01> NPT thermalization finished (from previous configuration) on 20 Nov 2021 at 15:04:37 -Simulation process simfiles/step03/p02 initiated with pid 15662 -p02> NPT thermalization finished (from previous configuration) on 20 Nov 2021 at 15:04:37 -Simulation process simfiles/step03/p04 initiated with pid 15664 -p04> NPT thermalization finished (from previous configuration) on 20 Nov 2021 at 15:04:37 -Simulation process simfiles/step03/p03 initiated with pid 15663 -p03> NPT thermalization finished (from previous configuration) on 20 Nov 2021 at 15:04:37 -p03> NPT production initiated on 20 Nov 2021 at 15:05:41 -p01> NPT production initiated on 20 Nov 2021 at 15:05:42 -p04> NPT production initiated on 20 Nov 2021 at 15:05:42 -p02> NPT production initiated on 20 Nov 2021 at 15:05:44 -p01> ----- NPT production finished on 20 Nov 2021 at 15:08:44 -p03> ----- NPT production finished on 20 Nov 2021 at 15:08:45 -p04> ----- NPT production finished on 20 Nov 2021 at 15:08:48 -p02> ----- NPT production finished on 20 Nov 2021 at 15:08:51 - -+----------------------------------------------------------------------------------------+ +Calculation of forces initiated with Gaussian on 03 Dec 2021 at 21:22:39 diff --git a/run.log.backup b/run.log.backup index 5d9ccc3..384d5c1 100644 --- a/run.log.backup +++ b/run.log.backup @@ -5,7 +5,7 @@ Your python version is 3.6.9 (default, Jan 26 2021, 15:33:00) [GCC 8.4.0] -Program started on Saturday, 20 Nov 2021 at 13:29:59 +Program started on Friday, 03 Dec 2021 at 21:18:49 Environment variables: OMP_STACKSIZE = Not set @@ -19,10 +19,10 @@ freq = no ghosts = no initcyc = 1 lps = no -maxcyc = 4 +maxcyc = 3 maxstep = 0.3 nprocs = 4 -opt = no +opt = yes qmprog = g16 readhessian = no switchcyc = 3 @@ -39,7 +39,7 @@ isave = 1000 ljname = phb.ljc ncores = 3 nmol = 1 50 -nstep = 40000 60000 50000 +nstep = 20000 20000 outname = phb press = 1.0 progname = dice @@ -115,8 +115,8 @@ Lbl AN X Y Z Charge Epsilon Sigma Mass 1 6 2.07203 -2.82345 0.00263 -0.115000 0.07000 3.5500 12.0110 1 6 2.76823 -1.61764 0.00263 -0.115000 0.07000 3.5500 12.0110 1 6 2.06824 -0.40521 0.00264 -0.115000 0.07000 3.5500 12.0110 -1 14 0.67589 -0.40522 0.00264 0.150000 0.07000 3.5500 28.0860 -1 104 -0.02420 -1.61760 0.00263 -0.115000 0.07000 3.5500 0.0000 +1 6 0.67589 -0.40522 0.00264 0.150000 0.07000 3.5500 12.0110 +1 6 -0.02420 -1.61760 0.00263 -0.115000 0.07000 3.5500 12.0110 2 1 0.13203 -3.75875 0.00263 0.115000 0.03000 2.4200 1.0079 2 1 2.61203 -3.75875 0.00263 0.115000 0.03000 2.4200 1.0079 2 1 2.60824 0.53010 0.00264 0.115000 0.03000 2.4200 1.0079 @@ -159,13 +159,13 @@ Molecule type 1 - PHENOL_BLUE: Molecule type 2 - PLACEHOLDER: - Center of mass = ( 1.6067 , -1.0929 , 0.0026 ) - Moments of inertia = 1.205279E+02 2.859254E+02 4.033761E+02 - Major principal axis = ( 0.795913 , -0.605411 , -0.000001 ) - Inter principal axis = ( 0.605411 , 0.795913 , 0.000001 ) + Center of mass = ( 1.5639 , -1.2534 , 0.0026 ) + Moments of inertia = 1.394830E+02 2.777126E+02 4.141185E+02 + Major principal axis = ( 0.865900 , -0.500217 , -0.000001 ) + Inter principal axis = ( 0.500217 , 0.865900 , 0.000001 ) Minor principal axis = ( 0.000000 , -0.000002 , 1.000000 ) Characteristic lengths = ( 5.74 , 5.26 , 1.51 ) - Total mass = 112.20 au + Total mass = 108.14 au Total charge = -0.0000 e Dipole moment = ( 2.3293 , 0.1593 , -0.0000 ) Total = 2.3347 Debye @@ -173,14 +173,14 @@ Molecule type 2 - PLACEHOLDER: New values: Center of mass = ( 0.0000 , -0.0000 , 0.0000 ) - Moments of inertia = 1.205279E+02 2.859254E+02 4.033761E+02 - Major principal axis = ( -1.000000 , 0.000000 , 0.000000 ) + Moments of inertia = 1.394830E+02 2.777126E+02 4.141185E+02 + Major principal axis = ( 1.000000 , -0.000000 , -0.000000 ) Inter principal axis = ( 0.000000 , 1.000000 , 0.000000 ) - Minor principal axis = ( 0.000000 , 0.000000 , 1.000000 ) - Characteristic lengths = ( 5.67 , 5.13 , 1.51 ) - Total mass = 112.20 au + Minor principal axis = ( 0.000000 , -0.000000 , 1.000000 ) + Characteristic lengths = ( 5.65 , 4.95 , 1.51 ) + Total mass = 108.14 au Total charge = -0.0000 e - Dipole moment = ( 1.7575 , 1.5369 , -0.0000 ) Total = 2.3347 Debye + Dipole moment = ( 1.9373 , 1.3031 , -0.0000 ) Total = 2.3347 Debye ========================================================================================== @@ -190,19 +190,11 @@ Starting the iterative process. Step # 1 ------------------------------------------------------------------------------------------ -Simulation process simfiles/step01/p01 initiated with pid 14014 -p01> NVT thermalization initiated (from random configuration) on 20 Nov 2021 at 13:29:59 -Simulation process simfiles/step01/p02 initiated with pid 14015 -p02> NVT thermalization initiated (from random configuration) on 20 Nov 2021 at 13:29:59 -Simulation process simfiles/step01/p03 initiated with pid 14016 -p03> NVT thermalization initiated (from random configuration) on 20 Nov 2021 at 13:29:59 -Simulation process simfiles/step01/p04 initiated with pid 14017 -p04> NVT thermalization initiated (from random configuration) on 20 Nov 2021 at 13:29:59 -p03> NPT thermalization finished (from previous configuration) on 20 Nov 2021 at 13:31:50 -p01> NPT thermalization finished (from previous configuration) on 20 Nov 2021 at 13:31:51 -p04> NPT thermalization finished (from previous configuration) on 20 Nov 2021 at 13:31:52 -p02> NPT thermalization finished (from previous configuration) on 20 Nov 2021 at 13:31:52 -p03> NPT production initiated on 20 Nov 2021 at 13:34:53 -p02> NPT production initiated on 20 Nov 2021 at 13:34:57 -p01> NPT production initiated on 20 Nov 2021 at 13:34:58 -p04> NPT production initiated on 20 Nov 2021 at 13:35:00 +Simulation process simfiles/step01/p01 initiated with pid 7513 +p01> NVT thermalization finished (from random configuration) on 03 Dec 2021 at 21:18:49 +Simulation process simfiles/step01/p02 initiated with pid 7514 +p02> NVT thermalization finished (from random configuration) on 03 Dec 2021 at 21:18:49 +Simulation process simfiles/step01/p03 initiated with pid 7515 +p03> NVT thermalization finished (from random configuration) on 03 Dec 2021 at 21:18:49 +Simulation process simfiles/step01/p04 initiated with pid 7516 +p04> NVT thermalization finished (from random configuration) on 03 Dec 2021 at 21:18:49