diff --git a/control.in b/control.in index 150585f..28820ac 100644 --- a/control.in +++ b/control.in @@ -1,6 +1,6 @@ # diceplayer maxcyc = 3 -opt = YES +opt = no nprocs = 4 qmprog = g16 lps = no @@ -18,4 +18,5 @@ outname = phb randominit = first # Gaussian -level = MP2/aug-cc-pVTZ +level = MP2/aug-cc-pVDZ +keywords = freq=intmodes \ No newline at end of file diff --git a/diceplayer/DPpack/.Dice.py b/diceplayer/DPpack/.Dice.py deleted file mode 100644 index b2aac92..0000000 --- a/diceplayer/DPpack/.Dice.py +++ /dev/null @@ -1,527 +0,0 @@ -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/diceplayer/DPpack/.MolHandling.py b/diceplayer/DPpack/.MolHandling.py deleted file mode 100644 index 87cc0e4..0000000 --- a/diceplayer/DPpack/.MolHandling.py +++ /dev/null @@ -1,600 +0,0 @@ -import sys, math -import textwrap -from copy import deepcopy - -import numpy as np -from numpy import linalg - -from DPpack.PTable import * -from DPpack.SetGlobals import * - - -####################################### functions ###################################### - -# def center_of_mass(molecule): - -# com = np.zeros(3) -# total_mass = 0.0 -# for atom in molecule: -# total_mass += atom['mass'] -# position = np.array([atom['rx'], atom['ry'], atom['rz']]) -# com += atom['mass'] * position - -# com = com / total_mass - -# return com - - - -# def center_of_mass_distance(molecule1, molecule2): - -# com1 = center_of_mass(molecule1) -# com2 = center_of_mass(molecule2) -# dx = com1[0] - com2[0] -# dy = com1[1] - com2[1] -# dz = com1[2] - com2[2] -# distance = math.sqrt(dx**2 + dy**2 + dz**2) - -# return distance - - - -# def center_of_mass_to_origin(molecule): - -# com = center_of_mass(molecule) -# for atom in molecule: -# atom['rx'] -= com[0] -# atom['ry'] -= com[1] -# atom['rz'] -= com[2] - -# return - - - -# def charges_and_dipole(molecule): - -# eA_to_Debye = 1/0.20819434 -# charge = 0 -# dipole = np.zeros(3) -# for atom in molecule: -# position = np.array([ atom['rx'], atom['ry'], atom['rz'] ]) -# dipole += atom['chg'] * position -# charge += atom['chg'] - -# dipole *= eA_to_Debye -# total_dipole = math.sqrt(dipole[0]**2 + dipole[1]**2 + dipole[2]**2) - -# return [charge, dipole[0], dipole[1], dipole[2], total_dipole] - - - -# def distances_between_atoms(molecule): - -# distances = [] -# dim = len(molecule) -# for atom1 in molecule: -# if atom1['na'] != ghost_number: -# for atom2 in molecule: -# 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 np.array(distances).reshape(dim, dim) - - - -# def eixos(molecule): - -# eixos = np.zeros(3) -# if len(molecule) == 2: -# position1 = np.array([ molecule[0]['rx'], molecule[0]['ry'], molecule[0]['rz'] ]) -# position2 = np.array([ molecule[1]['rx'], molecule[1]['ry'], molecule[1]['rz'] ]) -# eixos = position2 - position1 -# eixos /= linalg.norm(eixos) -# elif len(molecule) > 2: -# position1 = np.array([ molecule[0]['rx'], molecule[0]['ry'], molecule[0]['rz'] ]) -# position2 = np.array([ molecule[1]['rx'], molecule[1]['ry'], molecule[1]['rz'] ]) -# position3 = np.array([ molecule[2]['rx'], molecule[2]['ry'], molecule[2]['rz'] ]) -# v1 = position2 - position1 -# v2 = position3 - position1 -# v3 = np.cross(v1, v2) -# v2 = np.cross(v1, v3) -# v1 /= linalg.norm(v1) -# v2 /= linalg.norm(v2) -# v3 /= linalg.norm(v3) -# eixos = np.array([[v1[0], v1[1], v1[2]], -# [v2[0], v2[1], v2[2]], -# [v3[0], v3[1], v3[2]]]) - -# return eixos - - - -# def inertia_tensor(molecule): - -# com = center_of_mass(molecule) -# Ixx = Ixy = Ixz = Iyy = Iyz = Izz = 0.0 -# for atom in molecule: -# #### Obtain the displacement from the center of mass -# dx = atom['rx'] - com[0] -# dy = atom['ry'] - com[1] -# dz = atom['rz'] - com[2] -# #### Update the diagonal components of the tensor -# Ixx += atom['mass'] * (dy**2 + dz**2) -# Iyy += atom['mass'] * (dz**2 + dx**2) -# Izz += atom['mass'] * (dx**2 + dy**2) -# #### Update the off-diagonal components of the tensor -# Ixy += atom['mass'] * dx * dy * -1 -# Ixz += atom['mass'] * dx * dz * -1 -# Iyz += atom['mass'] * dy * dz * -1 - -# return np.array([[Ixx, Ixy, Ixz], -# [Ixy, Iyy, Iyz], -# [Ixz, Iyz, Izz]]) - - - -# def minimum_distance(molecule1, molecule2): - -# distances = [] -# for atom1 in molecule1: -# if atom1['na'] != ghost_number: -# for atom2 in molecule2: -# 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 nearest_image(refmol, molecule, lx, ly, lz, criterium="com"): - -# if criterium != "com" and criterium != "min": -# sys.exit("Error in value passed to function nearest_image") -# min_dist = 1e20 -# for i in range(-1, 2): -# for j in range(-1, 2): -# for k in range(-1, 2): - -# tr_vector = [i * lx, j * ly, k * lz] -# new_molecule = translate(molecule, tr_vector) -# if criterium == "com": -# dist = center_of_mass_distance(refmol, new_molecule) -# else: -# dist = minimum_distance(refmol, new_molecule) - -# if dist < min_dist: -# min_dist = dist -# nearestmol = deepcopy(new_molecule) - -# return min_dist, nearestmol - - - -# def calculate_step(gradient, hessian, fh): - -# invhessian = linalg.inv(hessian) -# pre_step = -1 * np.matmul(invhessian, gradient.T).T -# maxstep = np.amax(np.absolute(pre_step)) -# factor = min(1, player['maxstep']/maxstep) -# step = factor * pre_step - -# fh.write("\nCalculated step:\n") -# pre_step_list = pre_step.tolist() - -# fh.write("-----------------------------------------------------------------------\n" -# "Center Atomic Step (Bohr)\n" -# "Number Number X Y Z\n" -# "-----------------------------------------------------------------------\n") -# for i in range(len(molecules[0])): -# fh.write(" {:>5d} {:>3d} {:>14.9f} {:>14.9f} {:>14.9f}\n".format( -# i + 1, molecules[0][i]['na'], -# pre_step_list.pop(0), pre_step_list.pop(0), pre_step_list.pop(0))) - -# fh.write("-----------------------------------------------------------------------\n") - -# fh.write("Maximum step is {:>11.6}\n".format(maxstep)) -# fh.write("Scaling factor = {:>6.4f}\n".format(factor)) -# fh.write("\nFinal step (Bohr):\n") -# step_list = step.tolist() - -# fh.write("-----------------------------------------------------------------------\n" -# "Center Atomic Step (Bohr)\n" -# "Number Number X Y Z\n" -# "-----------------------------------------------------------------------\n") -# for i in range(len(molecules[0])): -# fh.write(" {:>5d} {:>3d} {:>14.9f} {:>14.9f} {:>14.9f}\n".format( -# i + 1, molecules[0][i]['na'], -# step_list.pop(0), step_list.pop(0), step_list.pop(0))) - -# fh.write("-----------------------------------------------------------------------\n") - -# step_max = np.amax(np.absolute(step)) -# step_rms = np.sqrt(np.mean(np.square(step))) - -# fh.write(" Max Step = {:>14.9f} RMS Step = {:>14.9f}\n\n".format( -# step_max, step_rms)) - -# return step - - - -# def read_position(molecule): - -# position_list = [] -# for atom in molecule: -# position_list.extend([ atom['rx'], atom['ry'], atom['rz'] ]) -# position = np.array(position_list) -# position *= ang2bohr - -# return position - - - -# def update_molecule(position, fh): - -# position_in_ang = (position * bohr2ang).tolist() -# new_molecule = deepcopy(molecules[0]) -# for atom in new_molecule: -# atom['rx'] = position_in_ang.pop(0) -# atom['ry'] = position_in_ang.pop(0) -# atom['rz'] = position_in_ang.pop(0) - -# rmsd, molecules[0] = rmsd_fit(new_molecule, molecules[0]) - -# fh.write("\nProjected new conformation of reference molecule with RMSD fit\n") -# fh.write("RMSD = {:>8.5f} Angstrom\n".format(rmsd)) - -# return - - - -# def update_hessian(step, cur_gradient, old_gradient, hessian): ## According to the BFGS - -# dif_gradient = cur_gradient - old_gradient - -# mat1 = 1/np.dot(dif_gradient, step) * np.matmul(dif_gradient.T, dif_gradient) -# mat2 = 1/np.dot(step, np.matmul(hessian, step.T).T) -# mat2 *= np.matmul( np.matmul(hessian, step.T), np.matmul(step, hessian) ) - -# hessian += mat1 - mat2 - -# return hessian - -### Antes de aplicar devemos aplicar a aplicação de Classes para Dice, Player e etc - -def populate_asec_vdw(cycle, fh): - - asec_charges = [] # (rx, ry, rz, chg) - vdw_meanfield = [] # (rx, ry, rz, eps, sig) - - if dice['nstep'][-1] % dice['isave'] == 0: - nconfigs = round(dice['nstep'][-1] / dice['isave']) - else: - nconfigs = int(dice['nstep'][-1] / dice['isave']) - - norm_factor = nconfigs * player['nprocs'] - - nsitesref = len(molecules[0]) + len(ghost_atoms) + len(lp_atoms) - - nsites_total = dice['nmol'][0] * nsitesref - for i in range(1, len(dice['nmol'])): - nsites_total += dice['nmol'][i] * len(molecules[i]) - - thickness = [] - picked_mols = [] - - for proc in range(1, player['nprocs'] + 1): ## Run over folders - - path = "step{:02d}".format(cycle) + os.sep + "p{:02d}".format(proc) - file = path + os.sep + dice['outname'] + ".xyz" - if not os.path.isfile(file): - sys.exit("Error: cannot find file {}".format(file)) - try: - with open(file) as xyzfh: - xyzfile = xyzfh.readlines() - except: - sys.exit("Error: cannot open file {}".format(file)) - - for config in range(nconfigs): ## Run over configs in a folder - - if int( xyzfile.pop(0).split()[0] ) != nsites_total: - sys.exit("Error: wrong number of sites in file {}".format(file)) - - box = xyzfile.pop(0).split()[-3:] - box = [ float(box[0]), float(box[1]), float(box[2]) ] - sizes = sizes_of_molecule(molecules[0]) - thickness.append( min([ (box[0] - sizes[0])/2, (box[1] - sizes[1])/2, - (box[2] - sizes[2])/2 ]) ) - - xyzfile = xyzfile[nsitesref:] ## Skip the first (reference) molecule - mol_count = 0 - for type in range(len(dice['nmol'])): ## Run over types of molecules - - if type == 0: - nmols = dice['nmol'][0] - 1 - else: - nmols = dice['nmol'][type] - - for mol in range(nmols): ## Run over molecules of each type - - new_molecule = [] - for site in range(len(molecules[type])): ## Run over sites of each molecule - - new_molecule.append({}) - line = xyzfile.pop(0).split() - - if line[0].title() != atomsymb[molecules[type][site]['na']].strip(): - sys.exit("Error reading file {}".format(file)) - - new_molecule[site]['na'] = molecules[type][site]['na'] - new_molecule[site]['rx'] = float(line[1]) - new_molecule[site]['ry'] = float(line[2]) - new_molecule[site]['rz'] = float(line[3]) - new_molecule[site]['chg'] = molecules[type][site]['chg'] - new_molecule[site]['eps'] = molecules[type][site]['eps'] - new_molecule[site]['sig'] = molecules[type][site]['sig'] - - dist = minimum_distance(molecules[0], new_molecule) - if dist < thickness[-1]: - mol_count += 1 - for atom in new_molecule: - asec_charges.append({}) - vdw_meanfield.append({}) - - asec_charges[-1]['rx'] = atom['rx'] - asec_charges[-1]['ry'] = atom['ry'] - asec_charges[-1]['rz'] = atom['rz'] - asec_charges[-1]['chg'] = atom['chg'] / norm_factor - - if player['vdwforces'] == "yes": - vdw_meanfield[-1]['rx'] = atom['rx'] - vdw_meanfield[-1]['ry'] = atom['ry'] - vdw_meanfield[-1]['rz'] = atom['rz'] - vdw_meanfield[-1]['eps'] = atom['eps'] - vdw_meanfield[-1]['sig'] = atom['sig'] - - #### Read lines with ghosts or lps in molecules of type 0 (reference) - #### and, if dist < thickness, appends to asec - if type == 0: - for ghost in ghost_atoms: - line = xyzfile.pop(0).split() - if line[0] != dice_ghost_label: - sys.exit("Error reading file {}".format(file)) - if dist < thickness[-1]: - asec_charges.append({}) - asec_charges[-1]['rx'] = float(line[1]) - asec_charges[-1]['ry'] = float(line[2]) - asec_charges[-1]['rz'] = float(line[3]) - asec_charges[-1]['chg'] = ghost['chg'] / norm_factor - - for lp in lp_atoms: - line = xyzfile.pop(0).split() - if line[0] != dice_ghost_label: - sys.exit("Error reading file {}".format(file)) - if dist < thickness[-1]: - asec_charges.append({}) - asec_charges[-1]['rx'] = float(line[1]) - asec_charges[-1]['ry'] = float(line[2]) - asec_charges[-1]['rz'] = float(line[3]) - asec_charges[-1]['chg'] = lp['chg'] / norm_factor - - picked_mols.append(mol_count) - - fh.write("Done\n") - - string = "In average, {:^7.2f} molecules ".format(sum(picked_mols)/norm_factor) - string += "were selected from each of the {} configurations ".format(len(picked_mols)) - string += "of the production simulations to form the ASEC, comprising a shell with " - string += "minimum thickness of {:>6.2f} Angstrom\n".format(sum(thickness)/norm_factor) - - fh.write(textwrap.fill(string, 86)) - fh.write("\n") - - otherfh = open("ASEC.dat", "w") - for charge in asec_charges: - otherfh.write("{:>10.5f} {:>10.5f} {:>10.5f} {:>11.8f}\n".format( - charge['rx'], charge['ry'], charge['rz'], charge['chg'])) - otherfh.close() - - return asec_charges - - - -# def principal_axes(inertia_tensor): - -# try: -# evals, evecs = linalg.eigh(inertia_tensor) -# except: -# sys.exit("Error: diagonalization of inertia tensor did not converge") - -# return evals, evecs - - - -# def print_geom(cycle, fh): - -# fh.write("{}\n".format(len(molecules[0]))) -# fh.write("Cycle # {}\n".format(cycle)) -# for atom in molecules[0]: -# symbol = atomsymb[atom['na']] -# fh.write("{:<2s} {:>10.6f} {:>10.6f} {:>10.6f}\n".format(symbol, -# atom['rx'], atom['ry'], atom['rz'])) - -# return - - - -# def print_mol_info(molecule, fh): - -# com = center_of_mass(molecule) -# fh.write(" Center of mass = ( {:>10.4f} , {:>10.4f} , {:>10.4f} )\n".format(com[0], -# com[1], com[2])) -# inertia = inertia_tensor(molecule) -# evals, evecs = principal_axes(inertia) - -# fh.write(" Moments of inertia = {:>9E} {:>9E} {:>9E}\n".format(evals[0], -# evals[1], evals[2])) - -# fh.write(" Major principal axis = ( {:>10.6f} , {:>10.6f} , {:>10.6f} )\n".format( -# evecs[0,0], evecs[1,0], evecs[2,0])) -# fh.write(" Inter principal axis = ( {:>10.6f} , {:>10.6f} , {:>10.6f} )\n".format( -# evecs[0,1], evecs[1,1], evecs[2,1])) -# fh.write(" Minor principal axis = ( {:>10.6f} , {:>10.6f} , {:>10.6f} )\n".format( -# evecs[0,2], evecs[1,2], evecs[2,2])) - -# sizes = sizes_of_molecule(molecule) -# fh.write(" Characteristic lengths = ( {:>6.2f} , {:>6.2f} , {:>6.2f} )\n".format( -# sizes[0], sizes[1], sizes[2])) -# mol_mass = total_mass(molecule) -# fh.write(" Total mass = {:>8.2f} au\n".format(mol_mass)) - -# chg_dip = charges_and_dipole(molecule) -# fh.write(" Total charge = {:>8.4f} e\n".format(chg_dip[0])) -# 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])) - -# return - - - -# def sizes_of_molecule(molecule): - -# x_list = [] -# y_list = [] -# z_list = [] -# for atom in molecule: -# if atom['na'] != ghost_number: -# x_list.append(atom['rx']) -# y_list.append(atom['ry']) -# z_list.append(atom['rz']) - -# x_max = max(x_list) -# x_min = min(x_list) -# y_max = max(y_list) -# y_min = min(y_list) -# z_max = max(z_list) -# z_min = min(z_list) - -# sizes = [x_max - x_min, y_max - y_min, z_max - z_min] - -# return sizes - - - -# def standard_orientation(molecule): - -# center_of_mass_to_origin(molecule) -# tensor = inertia_tensor(molecule) -# evals, evecs = principal_axes(tensor) -# if round(linalg.det(evecs)) == -1: -# evecs[0,2] *= -1 -# evecs[1,2] *= -1 -# evecs[2,2] *= -1 -# if round(linalg.det(evecs)) != 1: -# sys.exit("Error: could not make a rotation matrix while adopting the standard orientation") - -# rot_matrix = evecs.T -# for atom in molecule: -# position = np.array([ atom['rx'], atom['ry'], atom['rz'] ]) -# new_position = np.matmul(rot_matrix, position.T).T -# atom['rx'] = new_position[0] -# atom['ry'] = new_position[1] -# atom['rz'] = new_position[2] - -# return - - - -# def total_mass(molecule): - -# mass = 0 -# for atom in molecule: -# mass += atom['mass'] - -# return mass - - - -# def translate(molecule, vector): - -# new_molecule = deepcopy(molecule) -# for atom in new_molecule: -# atom['rx'] += vector[0] -# atom['ry'] += vector[1] -# atom['rz'] += vector[2] - -# return new_molecule - - - -# def rmsd_fit(projecting_mol, reference_mol): - -# if len(projecting_mol) != len(reference_mol): -# sys.exit("Error in RMSD fit procedure: molecules have different number of atoms") -# dim = len(projecting_mol) - -# new_projecting_mol = deepcopy(projecting_mol) -# new_reference_mol = deepcopy(reference_mol) - -# center_of_mass_to_origin(new_projecting_mol) -# center_of_mass_to_origin(new_reference_mol) - -# x = [] -# y = [] - -# for atom in new_projecting_mol: -# x.extend([ atom['rx'], atom['ry'], atom['rz'] ]) - -# for atom in new_reference_mol: -# y.extend([ atom['rx'], atom['ry'], atom['rz'] ]) - -# x = np.array(x).reshape(dim, 3) -# y = np.array(y).reshape(dim, 3) - -# r = np.matmul(y.T, x) -# rr = np.matmul(r.T, r) - -# try: -# evals, evecs = linalg.eigh(rr) -# except: -# sys.exit("Error: diagonalization of RR matrix did not converge") - -# a1 = evecs[:,2].T -# a2 = evecs[:,1].T -# a3 = np.cross(a1, a2) - -# A = np.array([ a1[0], a1[1], a1[2], a2[0], a2[1], a2[2], a3[0], a3[1], a3[2] ]) -# A = A.reshape(3,3) - -# b1 = np.matmul(r, a1.T).T # or np.dot(r, a1) -# b1 /= linalg.norm(b1) -# b2 = np.matmul(r, a2.T).T # or np.dot(r, a2) -# b2 /= linalg.norm(b2) -# b3 = np.cross(b1, b2) - -# B = np.array([ b1[0], b1[1], b1[2], b2[0], b2[1], b2[2], b3[0], b3[1], b3[2] ]) -# B = B.reshape(3,3).T - -# rot_matrix = np.matmul(B, A) -# x = np.matmul(rot_matrix, x.T).T - -# rmsd = 0 -# for i in range(dim): -# rmsd += (x[i,0] - y[i,0])**2 + (x[i,1] - y[i,1])**2 + (x[i,2] - y[i,2])**2 -# rmsd = math.sqrt(rmsd/dim) - -# for i in range(dim): -# new_projecting_mol[i]['rx'] = x[i,0] -# new_projecting_mol[i]['ry'] = x[i,1] -# new_projecting_mol[i]['rz'] = x[i,2] - -# tr_vector = center_of_mass(reference_mol) -# projected_mol = translate(new_projecting_mol, tr_vector) - -# return rmsd, projected_mol \ No newline at end of file diff --git a/diceplayer/DPpack/.SetGlobals.py b/diceplayer/DPpack/.SetGlobals.py deleted file mode 100644 index 6bd1317..0000000 --- a/diceplayer/DPpack/.SetGlobals.py +++ /dev/null @@ -1,813 +0,0 @@ -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/diceplayer/DPpack/Gaussian.py b/diceplayer/DPpack/Gaussian.py deleted file mode 100644 index 177b683..0000000 --- a/diceplayer/DPpack/Gaussian.py +++ /dev/null @@ -1,373 +0,0 @@ -import os, sys -import textwrap -import subprocess - -from DPpack.PTable import * -from DPpack.SetGlobals import * -from DPpack.MolHandling import * -from DPpack.Misc import * - -####################################### functions ###################################### - -def read_forces_fchk(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(molecules[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(molecules[0])): - fh.write(" {:>5d} {:>3d} {:>14.9f} {:>14.9f} {:>14.9f}\n".format( - i + 1, molecules[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(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(molecules[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(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(molecules[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(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(molecules[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_force_input(cycle, asec_charges): - - path = "step{:02d}".format(cycle) + 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 gaussian['mem'] != None: - fh.write("%Mem={}MB\n".format(gaussian['mem'])) - fh.write("%Nprocs={}\n".format(player['nprocs'] * dice['ncores'])) - - kword_line = "#P " + gaussian['level'] + " " + gaussian['keywords'] - kword_line += " Force Charge NoSymm" - - if cycle >= player['switchcyc']: - kword_line += " Pop={} Density=Current".format(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(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'])) - - ## If also performing charge fit in the same calculation - 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'])) - - 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'])) - - 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 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) - - fh.write("\n") - - fh.close() - - return - - - -# 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)) - -# 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" - -# if player['opt'] != "no" or cycle > 1: -# kword_line += " Guess=Read" - -# kword_line += " Pop={} Density=Current\n".format(gaussian['pop']) - -# 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])) - -# 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'])) - -# 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 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) - -# 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 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) - -# fh.write("\n") - -# fh.close() - -# return - - - -def read_charges(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 molecules[0]: - 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 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") - - return \ No newline at end of file diff --git a/diceplayer/DPpack/MolHandling.py b/diceplayer/DPpack/MolHandling.py index d0f8f50..6c4509f 100644 --- a/diceplayer/DPpack/MolHandling.py +++ b/diceplayer/DPpack/MolHandling.py @@ -1,31 +1,317 @@ -import os, sys -import math -import shutil -import textwrap -import sys, math -from copy import deepcopy - -import numpy as np -from numpy import linalg - -from diceplayer.DPpack.Misc import * from diceplayer.DPpack.PTable import * +from diceplayer.DPpack.Misc import * + +from typing import IO, Tuple, List, TextIO, Union + +from numpy import linalg +import numpy as np + +from copy import deepcopy +import sys, math +import textwrap +import os, sys +import shutil +import math env = ["OMP_STACKSIZE"] bohr2ang = 0.52917721092 ang2bohr = 1/bohr2ang +class Atom: + + def __init__(self, lbl: int, na: int, + rx: float,ry: float, rz: float, + chg: float, eps: float, sig: float + ) -> None: + + self.lbl = lbl # Integer + self.na = na # Integer + self.rx = rx # Double + self.ry = ry # Double + self.rz = rz # Double + self.chg = chg # Double + self.eps = eps # Double + self.sig = sig # Double + self.mass = atommass[self.na] # Double + +# Classe que conterá toda informação e funções relacionadas a uma unica molecula + +class Molecule: + + def __init__(self, molname: str ) -> None: + + self.molname = molname + + self.atom: List[Atom] = [] # Lista de instancias de Atom + self.position = None # Array Numpy + self.energy = None # Array Numpy + self.gradient = None # Array Numpy + self.hessian = None # Array Numpy + + self.total_mass = 0 + self.com = None + self.ghost_atoms: List[int] = [] # Stores the index of the ghost atoms in the atoms array + self.lp_atoms: List[int] = [] + + def add_atom(self, a: Atom) -> None: + + self.atom.append(a) # Inserção de um novo atomo + self.total_mass += a.mass + + if (a.na == ghost_number): + + self.ghost_atoms.append(self.atom.index(a)) + + self.center_of_mass() + + def center_of_mass(self) -> None: + + self.com = np.zeros(3) + + for atom in self.atom: + + self.com += atom.mass * np.array([atom.rx, atom.ry, atom.rz]) + + self.com = self.com / self.total_mass + + def center_of_mass_to_origin(self)-> None: + + self.center_of_mass() + + for atom in self.atom: + + atom.rx -= self.com[0] + atom.ry -= self.com[1] + atom.rz -= self.com[2] + + def charges_and_dipole(self) -> List[float]: + + eA_to_Debye = 1/0.20819434 + charge = 0 + dipole = np.zeros(3) + for atom in self.atom: + position = np.array([ atom.rx, atom.ry, atom.rz ]) + dipole += atom.chg * position + charge += atom.chg + + dipole *= eA_to_Debye + total_dipole = math.sqrt(dipole[0]**2 + dipole[1]**2 + dipole[2]**2) + + return [charge, dipole[0], dipole[1], dipole[2], total_dipole] + + def distances_between_atoms(self) -> np.ndarray: + + distances = [] + dim = len(self.atom) + for atom1 in self.atom: + if atom1.na != ghost_number: + for atom2 in self.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 np.array(distances).reshape(dim, dim) + + def inertia_tensor(self) -> np.ndarray: + + self.center_of_mass() + Ixx = Ixy = Ixz = Iyy = Iyz = Izz = 0.0 + + for atom in self.atom: + + #### Obtain the displacement from the center of mass + dx = atom.rx - self.com[0] + dy = atom.ry - self.com[1] + dz = atom.rz - self.com[2] + #### Update the diagonal components of the tensor + Ixx += atom.mass * (dy**2 + dz**2) + Iyy += atom.mass * (dz**2 + dx**2) + Izz += atom.mass * (dx**2 + dy**2) + #### Update the off-diagonal components of the tensor + Ixy += atom.mass * dx * dy * -1 + Ixz += atom.mass * dx * dz * -1 + Iyz += atom.mass * dy * dz * -1 + + return np.array([[Ixx, Ixy, Ixz], + [Ixy, Iyy, Iyz], + [Ixz, Iyz, Izz]]) + + def eixos(self) -> np.ndarray: + + eixos = np.zeros(3) + if len(self.atom) == 2: + + position1 = np.array([ self.atom[0].rx, self.atom[0].ry, self.atom[0].rz ]) + position2 = np.array([ self.atom[1].rx, self.atom[1].ry, self.atom[1].rz ]) + eixos = position2 - position1 + eixos /= linalg.norm(eixos) + + elif len(self.atom) > 2: + + position1 = np.array([ self.atom[0].rx, self.atom[0].ry, self.atom[0].rz ]) + position2 = np.array([ self.atom[1].rx, self.atom[1].ry, self.atom[1].rz ]) + position3 = np.array([ self.atom[2].rx, self.atom[2].ry, self.atom[2].rz ]) + v1 = position2 - position1 + v2 = position3 - position1 + v3 = np.cross(v1, v2) + v2 = np.cross(v1, v3) + v1 /= linalg.norm(v1) + v2 /= linalg.norm(v2) + v3 /= linalg.norm(v3) + eixos = np.array([[v1[0], v1[1], v1[2]], + [v2[0], v2[1], v2[2]], + [v3[0], v3[1], v3[2]]]) + + return eixos + + def principal_axes(self) -> Tuple[np.ndarray, np.ndarray]: + + try: + evals, evecs = linalg.eigh(self.inertia_tensor()) + except: + sys.exit("Error: diagonalization of inertia tensor did not converge") + + return evals, evecs + + def read_position(self) -> np.ndarray: + + position_list = [] + for atom in self.atom: + position_list.extend([ atom.rx, atom.ry, atom.rz ]) + position = np.array(position_list) + position *= ang2bohr + + return position + + def update_hessian(self, step: np.ndarray, + cur_gradient: np.ndarray, + old_gradient: np.ndarray, + hessian: np.ndarray + ) -> None: ## According to the BFGS + + dif_gradient = cur_gradient - old_gradient + + mat1 = 1/np.dot(dif_gradient, step) * np.matmul(dif_gradient.T, dif_gradient) + mat2 = 1/np.dot(step, np.matmul(hessian, step.T).T) + mat2 *= np.matmul( np.matmul(hessian, step.T), np.matmul(step, hessian) ) + + return hessian + mat1 - mat2 + + def sizes_of_molecule(self) -> List[float]: + + x_list = [] + y_list = [] + z_list = [] + + for atom in self.atom: + if atom.na != ghost_number: + x_list.append(atom.rx) + y_list.append(atom.ry) + z_list.append(atom.rz) + + x_max = max(x_list) + x_min = min(x_list) + y_max = max(y_list) + y_min = min(y_list) + z_max = max(z_list) + z_min = min(z_list) + + sizes = [x_max - x_min, y_max - y_min, z_max - z_min] + + return sizes + + def standard_orientation(self) -> None: + + self.center_of_mass_to_origin() + evals, evecs = self.principal_axes() + + if round(linalg.det(evecs)) == -1: + + evecs[0,2] *= -1 + evecs[1,2] *= -1 + evecs[2,2] *= -1 + + if round(linalg.det(evecs)) != 1: + + sys.exit("Error: could not make a rotation matrix while adopting the standard orientation") + + rot_matrix = evecs.T + + for atom in self.atom: + + position = np.array([ atom.rx, atom.ry, atom.rz ]) + new_position = np.matmul(rot_matrix, position.T).T + + atom.rx = new_position[0] + atom.ry = new_position[1] + atom.rz = new_position[2] + + def translate(self, vector: np.ndarray) -> 'Molecule': + + new_molecule = deepcopy(self) + + for atom in new_molecule.atom: + + atom.rx += vector[0] + atom.ry += vector[1] + atom.rz += vector[2] + + return new_molecule + + def print_mol_info(self, fh: TextIO) -> None: + + fh.write(" Center of mass = ( {:>10.4f} , {:>10.4f} , {:>10.4f} )\n".format(self.com[0], + self.com[1], self.com[2])) + inertia = self.inertia_tensor() + evals, evecs = self.principal_axes() + + fh.write(" Moments of inertia = {:>9E} {:>9E} {:>9E}\n".format(evals[0], + evals[1], evals[2])) + + fh.write(" Major principal axis = ( {:>10.6f} , {:>10.6f} , {:>10.6f} )\n".format( + evecs[0,0], evecs[1,0], evecs[2,0])) + fh.write(" Inter principal axis = ( {:>10.6f} , {:>10.6f} , {:>10.6f} )\n".format( + evecs[0,1], evecs[1,1], evecs[2,1])) + fh.write(" Minor principal axis = ( {:>10.6f} , {:>10.6f} , {:>10.6f} )\n".format( + evecs[0,2], evecs[1,2], evecs[2,2])) + + sizes = self.sizes_of_molecule() + fh.write(" Characteristic lengths = ( {:>6.2f} , {:>6.2f} , {:>6.2f} )\n".format( + sizes[0], sizes[1], sizes[2])) + fh.write(" Total mass = {:>8.2f} au\n".format(self.total_mass)) + + chg_dip = self.charges_and_dipole() + fh.write(" Total charge = {:>8.4f} e\n".format(chg_dip[0])) + 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: 'Molecule') -> List[float]: + + 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) + # Usaremos uma nova classe que ira conter toda interação entre moleculas class System: - def __init__(self): + def __init__(self) -> None: - self.molecule = [] - self.nmols = [] + self.molecule: List[Molecule] = [] + self.nmols: List[int] = [] - def add_type(self, nmols, m): + def add_type(self, nmols: int, m: Molecule) -> None: self.molecule.append(m) self.nmols.append(nmols) @@ -33,7 +319,7 @@ class System: # Função que calcula a distância entre dois centros de massa # e por se tratar de uma função de dois atomos não deve ser # inserida dentro de Molecule - def center_of_mass_distance(self, a, b): + def center_of_mass_distance(self, a: Molecule, b: Molecule) -> float: com1 = self.molecule[a].center_of_mass() com2 = self.molecule[b].center_of_mass() @@ -44,10 +330,10 @@ class System: return distance - def rmsd_fit(self, index_p, index_r): + def rmsd_fit(self, p_index: int, r_index: int) -> Tuple[float, Molecule]: - projecting_mol = self.molecule[index_p] - reference_mol = self.molecule[index_r] + projecting_mol = self.molecule[p_index] + reference_mol = self.molecule[r_index] if len(projecting_mol.atom) != len(reference_mol.atom): sys.exit("Error in RMSD fit procedure: molecules have different number of atoms") @@ -62,10 +348,10 @@ class System: x = [] y = [] - for atom in new_projecting_mol: + for atom in new_projecting_mol.atom: x.extend([ atom.rx, atom.ry, atom.rz ]) - for atom in new_reference_mol: + for atom in new_reference_mol.atom: y.extend([ atom.rx, atom.ry, atom.rz ]) x = np.array(x).reshape(dim, 3) @@ -108,15 +394,16 @@ class System: new_projecting_mol.atom[i].ry = x[i,1] new_projecting_mol.atom[i].rz = x[i,2] - tr_vector = reference_mol.center_of_mass() - projected_mol = new_projecting_mol.translate(tr_vector) + reference_mol.center_of_mass() + + projected_mol = new_projecting_mol.translate(reference_mol.com) return rmsd, projected_mol - def update_molecule(self, position, fh): + def update_molecule(self, position: np.ndarray, fh: TextIO) -> None: position_in_ang = (position * bohr2ang).tolist() - self.add_molecule(deepcopy(self.molecule[0])) + self.add_type(self.nmols[0],deepcopy(self.molecule[0])) for atom in self.molecule[-1].atom: @@ -130,7 +417,7 @@ class System: fh.write("\nProjected new conformation of reference molecule with RMSD fit\n") fh.write("RMSD = {:>8.5f} Angstrom\n".format(rmsd)) - def nearest_image(self, index_r, index_m, lx, ly, lz, criterium=None): + def nearest_image(self, index_r: int, index_m: int, lx: float, ly: float, lz: float, criterium=None) -> Tuple[float, Molecule]: if criterium in None: criterium = "com" @@ -160,10 +447,10 @@ class System: return min_dist, nearestmol - def print_geom(self, cycle, fh): + def print_geom(self, cycle: int, fh: TextIO) -> None: - fh.write("{}\n".format(len(self.molecule[0].atom))) fh.write("Cycle # {}\n".format(cycle)) + fh.write("Number of site: {}\n".format(len(self.molecule[0].atom))) for atom in self.molecule[0].atom: symbol = atomsymb[atom.na] fh.write("{:<2s} {:>10.6f} {:>10.6f} {:>10.6f}\n".format(symbol, @@ -171,279 +458,3 @@ class System: -# Classe que conterá toda informação e funções relacionadas a uma unica molecula - -class Molecule: - - def __init__(self, molname): - - self.molname = molname - - self.atom = [] # Lista de instancias de Atom - self.position = None # Array Numpy - self.energy = None # Array Numpy - self.gradient = None # Array Numpy - self.hessian = None # Array Numpy - - self.total_mass = 0 - self.com = None - self.ghost_atoms = [] # Stores the index of the ghost atoms in the atoms array - self.lp_atoms = [] - - def add_atom(self, a): - - self.atom.append(a) # Inserção de um novo atomo - self.total_mass += a.mass - - if (a.na == ghost_number): - - self.ghost_atoms.append(self.atom.index(a)) - - self.center_of_mass() - - def center_of_mass(self): - - self.com = np.zeros(3) - - for atom in self.atom: - - self.com += atom.mass * np.array([atom.rx, atom.ry, atom.rz]) - - self.com = self.com / self.total_mass - - def center_of_mass_to_origin(self): - - self.center_of_mass() - - for atom in self.atom: - - atom.rx -= self.com[0] - atom.ry -= self.com[1] - atom.rz -= self.com[2] - - def charges_and_dipole(self): - - eA_to_Debye = 1/0.20819434 - charge = 0 - dipole = np.zeros(3) - for atom in self.atom: - position = np.array([ atom.rx, atom.ry, atom.rz ]) - dipole += atom.chg * position - charge += atom.chg - - dipole *= eA_to_Debye - total_dipole = math.sqrt(dipole[0]**2 + dipole[1]**2 + dipole[2]**2) - - return [charge, dipole[0], dipole[1], dipole[2], total_dipole] - - def distances_between_atoms(self): - - distances = [] - dim = len(self.atom) - for atom1 in self.atom: - if atom1.na != ghost_number: - for atom2 in self.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 np.array(distances).reshape(dim, dim) - - def inertia_tensor(self): - - self.center_of_mass() - Ixx = Ixy = Ixz = Iyy = Iyz = Izz = 0.0 - - for atom in self.atom: - - #### Obtain the displacement from the center of mass - dx = atom.rx - self.com[0] - dy = atom.ry - self.com[1] - dz = atom.rz - self.com[2] - #### Update the diagonal components of the tensor - Ixx += atom.mass * (dy**2 + dz**2) - Iyy += atom.mass * (dz**2 + dx**2) - Izz += atom.mass * (dx**2 + dy**2) - #### Update the off-diagonal components of the tensor - Ixy += atom.mass * dx * dy * -1 - Ixz += atom.mass * dx * dz * -1 - Iyz += atom.mass * dy * dz * -1 - - return np.array([[Ixx, Ixy, Ixz], - [Ixy, Iyy, Iyz], - [Ixz, Iyz, Izz]]) - - def eixos(self): - - eixos = np.zeros(3) - if len(self.atom) == 2: - - position1 = np.array([ self.atom[0].rx, self.atom[0].ry, self.atom[0].rz ]) - position2 = np.array([ self.atom[1].rx, self.atom[1].ry, self.atom[1].rz ]) - eixos = position2 - position1 - eixos /= linalg.norm(eixos) - - elif len(self.atom) > 2: - - position1 = np.array([ self.atom[0].rx, self.atom[0].ry, self.atom[0].rz ]) - position2 = np.array([ self.atom[1].rx, self.atom[1].ry, self.atom[1].rz ]) - position3 = np.array([ self.atom[2].rx, self.atom[2].ry, self.atom[2].rz ]) - v1 = position2 - position1 - v2 = position3 - position1 - v3 = np.cross(v1, v2) - v2 = np.cross(v1, v3) - v1 /= linalg.norm(v1) - v2 /= linalg.norm(v2) - v3 /= linalg.norm(v3) - eixos = np.array([[v1[0], v1[1], v1[2]], - [v2[0], v2[1], v2[2]], - [v3[0], v3[1], v3[2]]]) - - return eixos - - def principal_axes(self): - - try: - evals, evecs = linalg.eigh(self.inertia_tensor()) - except: - sys.exit("Error: diagonalization of inertia tensor did not converge") - - return evals, evecs - - def read_position(self): - - position_list = [] - for atom in self.atom: - position_list.extend([ atom.rx, atom.ry, atom.rz ]) - position = np.array(position_list) - position *= ang2bohr - - return position - - def update_hessian(self, step, cur_gradient): ## According to the BFGS - - dif_gradient = cur_gradient - self.gradient - - mat1 = 1/np.dot(dif_gradient, step) * np.matmul(dif_gradient.T, dif_gradient) - mat2 = 1/np.dot(step, np.matmul(self.hessian, step.T).T) - mat2 *= np.matmul( np.matmul(self.hessian, step.T), np.matmul(step, self.hessian) ) - - self.hessian += mat1 - mat2 - - def sizes_of_molecule(self): - - x_list = [] - y_list = [] - z_list = [] - - for atom in self.atom: - if atom.na != ghost_number: - x_list.append(atom.rx) - y_list.append(atom.ry) - z_list.append(atom.rz) - - x_max = max(x_list) - x_min = min(x_list) - y_max = max(y_list) - y_min = min(y_list) - z_max = max(z_list) - z_min = min(z_list) - - sizes = [x_max - x_min, y_max - y_min, z_max - z_min] - - return sizes - - def standard_orientation(self): - - self.center_of_mass_to_origin() - evals, evecs = self.principal_axes() - - if round(linalg.det(evecs)) == -1: - - evecs[0,2] *= -1 - evecs[1,2] *= -1 - evecs[2,2] *= -1 - - if round(linalg.det(evecs)) != 1: - - sys.exit("Error: could not make a rotation matrix while adopting the standard orientation") - - rot_matrix = evecs.T - - for atom in self.atom: - - position = np.array([ atom.rx, atom.ry, atom.rz ]) - new_position = np.matmul(rot_matrix, position.T).T - - atom.rx = new_position[0] - atom.ry = new_position[1] - atom.rz = new_position[2] - - def translate(self, vector): - - new_molecule = deepcopy(self) - - for atom in new_molecule.atom: - - atom.rx += vector[0] - atom.ry += vector[1] - atom.rz += vector[2] - - return new_molecule - - def print_mol_info(self, fh): - - fh.write(" Center of mass = ( {:>10.4f} , {:>10.4f} , {:>10.4f} )\n".format(self.com[0], - self.com[1], self.com[2])) - inertia = self.inertia_tensor() - evals, evecs = self.principal_axes() - - fh.write(" Moments of inertia = {:>9E} {:>9E} {:>9E}\n".format(evals[0], - evals[1], evals[2])) - - fh.write(" Major principal axis = ( {:>10.6f} , {:>10.6f} , {:>10.6f} )\n".format( - evecs[0,0], evecs[1,0], evecs[2,0])) - fh.write(" Inter principal axis = ( {:>10.6f} , {:>10.6f} , {:>10.6f} )\n".format( - evecs[0,1], evecs[1,1], evecs[2,1])) - fh.write(" Minor principal axis = ( {:>10.6f} , {:>10.6f} , {:>10.6f} )\n".format( - evecs[0,2], evecs[1,2], evecs[2,2])) - - sizes = self.sizes_of_molecule() - fh.write(" Characteristic lengths = ( {:>6.2f} , {:>6.2f} , {:>6.2f} )\n".format( - sizes[0], sizes[1], sizes[2])) - fh.write(" Total mass = {:>8.2f} au\n".format(self.total_mass)) - - chg_dip = self.charges_and_dipole() - fh.write(" Total charge = {:>8.4f} e\n".format(chg_dip[0])) - 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): - - self.lbl = lbl # Integer - self.na = na # Integer - self.rx = rx # Double - self.ry = ry # Double - self.rz = rz # Double - self.chg = chg # Double - self.eps = eps # Double - self.sig = sig # Double - self.mass = atommass[self.na] # Double \ No newline at end of file diff --git a/diceplayer/DPpack/SetGlobals.py b/diceplayer/DPpack/SetGlobals.py index 4bd5381..25c6baa 100644 --- a/diceplayer/DPpack/SetGlobals.py +++ b/diceplayer/DPpack/SetGlobals.py @@ -1,17 +1,19 @@ -import setproctitle -import os, sys -import shutil -import textwrap -import types - -from numpy.core.fromnumeric import partition - from diceplayer.DPpack.MolHandling import * from diceplayer.DPpack.PTable import * from diceplayer.DPpack.Misc import * + +from typing import IO, Tuple, List, TextIO, Union + +from numpy.core.numeric import partition from numpy import random + +import setproctitle import subprocess +import os, sys +import shutil +import textwrap +import types dice_end_flag = "End of simulation" ## The normal end flag @@ -22,7 +24,7 @@ max_seed = 4294967295 ## Maximum allowed value for a seed (numpy) class Internal: - def __init__(self, infile, outfile): + def __init__(self, infile: TextIO, outfile: TextIO) -> None: self.cyc = 1 self.infile = infile @@ -53,7 +55,7 @@ class Internal: ## Dice: self.combrule = None - def read_keywords(self): + def read_keywords(self) -> None: try: controlfile = self.infile.readlines() @@ -234,7 +236,7 @@ class Internal: # #### End - def check_keywords(self): + def check_keywords(self) -> None: min_steps = 20000 @@ -309,7 +311,7 @@ class Internal: # isave value is always the nearest multiple of 100 self.dice.isave = round(self.dice.isave / 100) * 100 - def print_keywords(self): + def print_keywords(self) -> None: self.outfile.write("##########################################################################################\n" "############# Welcome to DICEPLAYER version 1.0 #############\n" @@ -388,7 +390,7 @@ class Internal: # self.outfile.write("\n") - def read_potential(self): # Deve ser atualizado para o uso de + def read_potential(self) -> None: # Deve ser atualizado para o uso de try: with open(self.dice.ljname) as file: @@ -496,7 +498,7 @@ class Internal: if _var in locals() or _var in globals(): exec(f'del {_var}') - def print_potential(self): + def print_potential(self) -> None: formatstr = "{:<3d} {:>3d} {:>10.5f} {:>10.5f} {:>10.5f} {:>10.6f} {:>9.5f} {:>7.4f} {:>9.4f}\n" self.outfile.write("\n" @@ -577,7 +579,7 @@ class Internal: self.outfile.write("\n" "==========================================================================================\n") - def check_executables(self): + def check_executables(self) -> None: self.outfile.write("\n") self.outfile.write(90 * "=") @@ -604,64 +606,53 @@ class Internal: else: sys.exit("Error: cannot find formchk executable") - def calculate_step(self): - - invhessian = linalg.inv(self.system.molecule[0].hessian) - pre_step = -1 * np.matmul(invhessian, self.system.molecule[0].gradient.T).T + def calculate_step(self, cycle: int, gradient: np.ndarray, + hessian: np.ndarray) -> np.ndarray: + + invhessian = np.linalg.inv(hessian) + pre_step = -1 * np.matmul(invhessian, gradient.T).T maxstep = np.amax(np.absolute(pre_step)) factor = min(1, self.player.maxstep/maxstep) step = factor * pre_step - self.player.outfile.write("\nCalculated step:\n") + self.outfile.write("\nCalculated step-{}:\n".format(cycle)) pre_step_list = pre_step.tolist() - self.player.outfile.write("-----------------------------------------------------------------------\n" + self.outfile.write("-----------------------------------------------------------------------\n" "Center Atomic Step (Bohr)\n" "Number Number X Y Z\n" "-----------------------------------------------------------------------\n") for i in range(len(self.system.molecule[0].atom)): - self.player.outfile.write(" {:>5d} {:>3d} {:>14.9f} {:>14.9f} {:>14.9f}\n".format( + self.outfile.write(" {:>5d} {:>3d} {:>14.9f} {:>14.9f} {:>14.9f}\n".format( i + 1, self.system.molecule[0].atom[i].na, pre_step_list.pop(0), pre_step_list.pop(0), pre_step_list.pop(0))) - self.player.outfile.write("-----------------------------------------------------------------------\n") + self.outfile.write("-----------------------------------------------------------------------\n") - self.player.outfile.write("Maximum step is {:>11.6}\n".format(maxstep)) - self.player.outfile.write("Scaling factor = {:>6.4f}\n".format(factor)) - self.player.outfile.write("\nFinal step (Bohr):\n") + self.outfile.write("Maximum step is {:>11.6}\n".format(maxstep)) + self.outfile.write("Scaling factor = {:>6.4f}\n".format(factor)) + self.outfile.write("\nFinal step (Bohr):\n") step_list = step.tolist() - self.player.outfile.write("-----------------------------------------------------------------------\n" + self.outfile.write("-----------------------------------------------------------------------\n" "Center Atomic Step (Bohr)\n" "Number Number X Y Z\n" "-----------------------------------------------------------------------\n") for i in range(len(self.system.molecule[0].atom)): - self.player.outfile.write(" {:>5d} {:>3d} {:>14.9f} {:>14.9f} {:>14.9f}\n".format( + self.outfile.write(" {:>5d} {:>3d} {:>14.9f} {:>14.9f} {:>14.9f}\n".format( i + 1, self.system.molecule[0].atom[i].na, step_list.pop(0), step_list.pop(0), step_list.pop(0))) - self.player.outfile.write("-----------------------------------------------------------------------\n") + self.outfile.write("-----------------------------------------------------------------------\n") step_max = np.amax(np.absolute(step)) step_rms = np.sqrt(np.mean(np.square(step))) - self.player.outfile.write(" Max Step = {:>14.9f} RMS Step = {:>14.9f}\n\n".format( + self.outfile.write(" Max Step = {:>14.9f} RMS Step = {:>14.9f}\n\n".format( step_max, step_rms)) return step - def read_initial_cicle(self): - - try: - with open(self.infile) as self.outfile: - controlfile = self.outfile.readlines() - except EnvironmentError: - sys.exit("Error: cannot open file {}".format(self.infile)) - - for line in controlfile: - - pass - ### I still have to talk with Herbet about this function def populate_asec_vdw(self, cycle): @@ -785,15 +776,15 @@ class Internal: picked_mols.append(mol_count) - self.player.outfile.write("Done\n") + self.outfile.write("Done\n") string = "In average, {:^7.2f} molecules ".format(sum(picked_mols)/norm_factor) string += "were selected from each of the {} configurations ".format(len(picked_mols)) string += "of the production simulations to form the ASEC, comprising a shell with " string += "minimum thickness of {:>6.2f} Angstrom\n".format(sum(thickness)/norm_factor) - self.player.outfile.write(textwrap.fill(string, 86)) - self.player.outfile.write("\n") + self.outfile.write(textwrap.fill(string, 86)) + self.outfile.write("\n") otherfh = open("ASEC.dat", "w") for charge in asec_charges: @@ -805,7 +796,7 @@ class Internal: ## Dice related Upper fuctions - def print_last_config(self, cycle, proc): + def print_last_config(self, cycle: int, proc: int) -> None: sim_dir = "simfiles" step_dir = "step{:02d}".format(cycle) @@ -835,7 +826,7 @@ class Internal: for line in xyzfile: fh.write(line) - def new_density(self, cycle, proc): + def new_density(self, cycle: int, proc: int) -> float: sim_dir = "simfiles" step_dir = "step{:02d}".format(cycle-1) @@ -856,13 +847,13 @@ class Internal: total_mass = 0 for i in range(len(self.system.molecule)): - total_mass += self.system.molecule[i].total_mass * self.dice.nmol[i] + total_mass += self.system.molecule[i].total_mass * self.system.nmols[i] density = (total_mass / volume) * umaAng3_to_gcm3 return density - def simulation_process(self, cycle, proc): + def simulation_process(self, cycle: int, proc: int) -> None: setproctitle.setproctitle("diceplayer-step{:0d}-p{:0d}".format(cycle,proc)) @@ -873,7 +864,7 @@ class Internal: except Exception as err: sys.exit(err) - def make_dice_inputs(self, cycle, proc): + def make_dice_inputs(self, cycle: int, proc: int) -> None: sim_dir = "simfiles" step_dir = "step{:02d}".format(cycle) @@ -917,7 +908,7 @@ class Internal: # shutil.copyfile(last_path + os.sep + "phb.dat", path + os.sep + "phb.dat") - def make_nvt_ter(self,cycle, path): + def make_nvt_ter(self,cycle: int, path: str) -> None: file = path + os.sep + "NVT.ter" try: @@ -957,7 +948,7 @@ class Internal: fh.close() - def make_nvt_eq(self, path): + def make_nvt_eq(self, path: str) -> None: file = path + os.sep + "NVT.eq" try: @@ -989,7 +980,7 @@ class Internal: fh.close() - def make_npt_ter(self, cycle, path): + def make_npt_ter(self, cycle: int, path: str) -> None: file = path + os.sep + "NPT.ter" try: @@ -1029,7 +1020,7 @@ class Internal: fh.close() - def make_npt_eq(self, path): + def make_npt_eq(self, path: str) -> None: file = path + os.sep + "NPT.eq" try: @@ -1063,7 +1054,7 @@ class Internal: fh.close() - def make_init_file(self, path, file): + def make_init_file(self, path: str, file: TextIO) -> None: if not os.path.isfile(file): sys.exit("Error: cannot find the xyz file {} in main directory".format(file)) @@ -1109,7 +1100,7 @@ class Internal: fh.close() - def make_potential(self, path): + def make_potential(self, path: str) -> None: fstr = "{:<3d} {:>3d} {:>10.5f} {:>10.5f} {:>10.5f} {:>10.6f} {:>9.5f} {:>7.4f}\n" @@ -1146,10 +1137,11 @@ class Internal: for atom in mol.atom: 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): + #---------------------------------------------------------------------------------------------------------------------------------------- + def read_forces_fchk(self, file: str, fh: TextIO) -> np.ndarray: forces = [] try: @@ -1162,9 +1154,9 @@ class Internal: while start.find("Cartesian Gradient") != 0: ## expression in begining of line start = fchkfile.pop(0).strip() - degrees = 3 * len(self.system.molecule[0]) + degrees = 3 * len(self.system.molecule[0].atom) count = 0 - while True: + while len(forces) < degrees: values = fchkfile.pop(0).split() forces.extend([ float(x) for x in values ]) count += len(values) @@ -1179,9 +1171,9 @@ class Internal: "Center Atomic Forces (Hartree/Bohr)\n" "Number Number X Y Z\n" "-----------------------------------------------------------------------\n") - for i in range(len(self.system.molecule[0])): + for i in range(len(self.system.molecule[0].atom)): 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))) + i + 1, self.system.molecule[0].atom[i].na, forces.pop(0), forces.pop(0), forces.pop(0))) fh.write("-----------------------------------------------------------------------\n") @@ -1193,9 +1185,7 @@ class Internal: return gradient - - - def read_hessian_fchk(self, file): + def read_hessian_fchk(self, file: str) -> np.ndarray: force_const = [] try: @@ -1208,16 +1198,23 @@ class Internal: while start.find("Cartesian Force Constants") != 0: start = fchkfile.pop(0).strip() - degrees = 3 * len(self.system.molecule[0]) + degrees = 3 * len(self.system.molecule[0].atom) 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 + + while len(force_const) < last: + + value = fchkfile.pop(0).split() + force_const.extend([ float(x) for x in value ]) + + # while len(force_const) < last: + + # 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): @@ -1227,9 +1224,7 @@ class Internal: return hessian - - - def read_hessian_log(self, file): + def read_hessian_log(self, file: str) -> np.ndarray: try: with open(file) as tmpfh: @@ -1241,7 +1236,7 @@ class Internal: while start.find("The second derivative matrix:") != 0: start = logfile.pop(0).strip() - degrees = 3 * len(self.system.molecule[0]) + degrees = 3 * len(self.system.molecule[0].atom) hessian = np.zeros((degrees, degrees)) k = 0 @@ -1256,9 +1251,9 @@ class Internal: return hessian - - - def print_grad_hessian(self, cycle, cur_gradient, hessian): + def print_grad_hessian(self, cycle: int, + cur_gradient: np.ndarray, hessian: np.ndarray + ) -> None: try: fh = open("grad_hessian.dat", "w") @@ -1267,16 +1262,17 @@ class Internal: fh.write("Optimization cycle: {}\n".format(cycle)) fh.write("Cartesian Gradient\n") - degrees = 3 * len(self.system.molecule[0]) + degrees = 3 * len(self.system.molecule[0].atom) 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") + n = int(np.sqrt(2*degrees)) last = degrees * (degrees + 1) / 2 count = 0 - for i in range(degrees): + for i in range(n): for j in range(i + 1): count += 1 fh.write(" {:>11.8g}".format(hessian[i,j])) @@ -1285,11 +1281,8 @@ class Internal: fh.close() - return - - ## Change the name to make_gaussian_input - def make_gaussian_input(self, cycle, asec_charges=None): + def make_gaussian_input(self, cycle: int, asec_charges=None) -> None: simdir="simfiles" stepdir="step{:02d}".format(cycle) @@ -1333,7 +1326,7 @@ class Internal: 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: @@ -1387,7 +1380,7 @@ class Internal: # fh.close() - def read_charges(self, file, fh): + def read_charges(self, file: str, fh: TextIO) -> None: try: with open(file) as tmpfh: @@ -1429,7 +1422,7 @@ class Internal: class Player: - def __init__(self): + def __init__(self) -> None: self.maxcyc = None self.nprocs = 1 @@ -1449,7 +1442,7 @@ class Internal: class Dice: - def __init__(self): + def __init__(self) -> None: self.title = "Diceplayer run" self.progname = "dice" @@ -1468,13 +1461,13 @@ class Internal: 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 + self.nmol: List[int] = [] # Up to 4 integer values related to up to 4 molecule types + self.nstep: List[int] = [] # 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 self.upbuf = 360 - def make_proc_dir(self, cycle, proc): + def make_proc_dir(self, cycle: int, proc: int) -> None: sim_dir = "simfiles" step_dir = "step{:02d}".format(cycle) @@ -1485,7 +1478,7 @@ class Internal: except: sys.exit("Error: cannot make directory {}".format(path)) - def run_dice(self, cycle, proc, fh): + def run_dice(self, cycle: int, proc: int, fh: str) -> None: sim_dir = "simfiles" step_dir = "step{:02d}".format(cycle) @@ -1528,13 +1521,13 @@ class Internal: sys.exit() if exit_status != 0: - sys.exit("Dice process p{:02d} did not exit properly".format(proc)) + sys.exit("Dice process step{:02d}-p{:02d} did not exit properly".format(cycle,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)) + sys.exit("Dice process step{:02d}-p{:02d} did not exit properly".format(cycle,proc)) ## NVT production fh.write("p{:02d}> NVT production initiated on {}\n".format(proc, date_time())) @@ -1554,13 +1547,13 @@ class Internal: sys.exit() if exit_status != 0: - sys.exit("Dice process p{:02d} did not exit properly".format(proc)) + sys.exit("Dice process step{:02d}-p{:02d} did not exit properly".format(cycle,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)) + sys.exit("Dice process step{:02d}-p{:02d} did not exit properly".format(cycle,proc)) fh.write("p{:02d}> ----- NVT production finished on {}\n".format(proc, date_time())) @@ -1587,13 +1580,13 @@ class Internal: sys.exit() if exit_status != 0: - sys.exit("Dice process p{:02d} did not exit properly".format(proc)) + sys.exit("Dice process step{:02d}-p{:02d} did not exit properly".format(cycle,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)) + sys.exit("Dice process step{:02d}-p{:02d} did not exit properly".format(cycle,proc)) ## NPT thermalization if not self.randominit == 'always' or (self.randominit == 'first' and cycle == 1): @@ -1618,13 +1611,13 @@ class Internal: sys.exit() if exit_status != 0: - sys.exit("Dice process p{:02d} did not exit properly".format(proc)) + sys.exit("Dice process step{:02d}-p{:02d} did not exit properly".format(cycle,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)) + sys.exit("Dice process step{:02d}-p{:02d} did not exit properly".format(cycle,proc)) ## NPT production fh.write("p{:02d}> NPT production initiated on {}\n".format(proc, date_time())) @@ -1644,13 +1637,13 @@ class Internal: sys.exit() if exit_status != 0: - sys.exit("Dice process p{:02d} did not exit properly".format(proc)) + sys.exit("Dice process step{:02d}-p{:02d} did not exit properly".format(cycle,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)) + sys.exit("Dice process step{:02d}-p{:02d} did not exit properly".format(cycle,proc)) fh.write("p{:02d}> ----- NPT production finished on {}\n".format(proc, date_time())) @@ -1659,7 +1652,7 @@ class Internal: class Gaussian: - def __init__(self): + def __init__(self) -> None: self.qmprog = "g09" self.path = None @@ -1672,7 +1665,7 @@ class Internal: self.pop = "chelpg" self.level = None - def run_gaussian(self, cycle, type, fh): + def run_gaussian(self, cycle: int, type: str, fh: TextIO) -> None: simdir="simfiles" stepdir="step{:02d}".format(cycle) @@ -1701,7 +1694,7 @@ class Internal: os.chdir(work_dir) - def run_formchk(self, cycle, fh): + def run_formchk(self, cycle: int, fh: TextIO): simdir="simfiles" stepdir="step{:02d}".format(cycle) @@ -1710,9 +1703,9 @@ class Internal: work_dir = os.getcwd() os.chdir(path) - fh.write("Formatting the checkpoint file... ") + fh.write("Formatting the checkpoint file... \n") - exit_status = subprocess.call(["formchk", "asec.chk"]) + exit_status = subprocess.call(["formchk", "asec.chk"], stdout=fh) fh.write("Done\n") diff --git a/diceplayer/__main__.py b/diceplayer/__main__.py index 363985c..eab4667 100644 --- a/diceplayer/__main__.py +++ b/diceplayer/__main__.py @@ -2,6 +2,7 @@ import os, sys, time, signal import setproctitle +import numpy as np import argparse import shutil from multiprocessing import Process, connection @@ -134,6 +135,7 @@ if __name__ == '__main__': except EnvironmentError as err: sys.exit(err) internal.system.print_geom(0, geomsfh) + geomsfh.write(40 * "-"+"\n") else: try: path = "geoms.xyz" @@ -196,8 +198,6 @@ if __name__ == '__main__': for proc in range(1, internal.player.nprocs + 1): internal.print_last_config(cycle, proc) - - internal.outfile.write("\n+" + 88 * "-" + "+\n") ### ### End of parallel simulations block @@ -221,47 +221,62 @@ if __name__ == '__main__': if internal.player.qmprog in ("g03", "g09", "g16"): if cycle > 1: + src = "simfiles" + os.sep + "step{:02d}".format(cycle - 1) + os.sep + "qm" + os.sep + "asec.chk" dst = "simfiles" + os.sep + "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec.chk" shutil.copyfile(src, dst) internal.make_gaussian_input(cycle) internal.gaussian.run_gaussian(cycle, "force", internal.outfile) -# internal.gaussian.run_formchk(cycle, internal.outfile) + internal.gaussian.run_formchk(cycle, internal.outfile) -# ## Read the gradient -# 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 + ## Read the gradient + file = "simfiles" + os.sep + "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec.fchk" + + try: + gradient + old_gradient = gradient + except: + pass + + gradient = internal.read_forces_fchk(file, internal.outfile) + + print(type(gradient),"\n",gradient.shape,"\n",gradient) -# ## If 1st step, read the hessian -# if cycle == 1: -# if internal.player.readhessian == "yes": -# file = "grad_hessian.dat" -# internal.outfile.write("\nReading the hessian matrix from file {}\n".format(file)) -# hessian = internal.read_hessian_fchk(file) -# else: -# 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) + # If 1st step, read the hessian + if cycle == 1: + + if internal.player.readhessian == "yes": + + file = "grad_hessian.dat" + internal.outfile.write("\nReading the hessian matrix from file {}\n".format(file)) + hessian = internal.read_hessian_log(file) + + else: + + 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.read_hessian_fchk(file) + + print(type(hessian), "\n", hessian.shape, "\n", hessian) -# ## From 2nd step on, update the hessian -# else: -# internal.outfile.write("\nUpdating the hessian matrix using the BFGS method... ") -# hessian = internal.system.molecule[0].update_hessian(step, cur_gradient, old_gradient, hessian) -# internal.outfile.write("Done\n") + # From 2nd step on, update the hessian + else: + internal.outfile.write("\nUpdating the hessian matrix using the BFGS method... ") + hessian = internal.system.molecule[0].update_hessian(step, gradient, old_gradient, hessian) + internal.outfile.write("Done\n") -# ## Save gradient and hessian -# internal.gaussian.print_grad_hessian(cycle, cur_gradient, hessian) + # Save gradient and hessian + internal.print_grad_hessian(cycle, gradient, hessian) + + # Calculate the step and update the position + step = internal.calculate_step(cycle, gradient, hessian) + + print(type(step), "\n", step.shape, "\n", step) + position += step -# ## Calculate the step and update the position -# step = internal.calculate_step(cur_gradient, hessian, internal.outfile) -# position += step - -# ## Update the geometry of the reference molecule -# internal.system.molecule[0].update_molecule(position, internal.outfile) + # ## Update the geometry of the reference molecule + internal.system.update_molecule(position, internal.outfile) # ## If needed, calculate the charges # if cycle < internal.player.switchcyc: @@ -326,17 +341,19 @@ if __name__ == '__main__': # #if player['qmprog'] == "molcas": - - -# #### -# #### End of the iterative process -# #### + internal.system.print_geom(cycle, geomsfh) + geomsfh.write(40 * "-"+"\n") -# # ## imprimir ultimas mensagens, criar um arquivo de potencial para ser usado em eventual -# # ## continuacao, fechar arquivos (geoms.xyz, run.log, ...) + internal.outfile.write("\n+" + 88 * "-" + "+\n") + #### + #### End of the iterative process + #### -# # internal.outfile.write("\nDiceplayer finished normally!\n") -# # internal.outfile.close() -# # #### -# # #### End of the program -# # #### +## 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 +#### diff --git a/geoms.xyz b/geoms.xyz deleted file mode 100644 index 74643c1..0000000 --- a/geoms.xyz +++ /dev/null @@ -1,33 +0,0 @@ -31 -Cycle # 0 -C -4.394497 0.698769 0.421354 -C -4.504613 -0.640849 -0.141092 -C -3.266831 -1.267916 -0.601745 -C -2.064301 -0.654971 -0.479838 -C -1.966148 0.666069 0.106967 -C -3.195469 1.321761 0.499266 -O -5.604026 -1.213715 -0.270977 -N -0.863607 1.365979 0.273247 -C 0.394047 0.812820 0.140269 -C 1.418955 1.646567 -0.341744 -C 0.757396 -0.471466 0.587991 -C 2.717900 1.193894 -0.471504 -C 2.068121 -0.911914 0.507136 -C 3.079906 -0.111825 -0.067434 -N 4.367937 -0.572075 -0.209272 -C 5.421862 0.397502 -0.451662 -C 4.752115 -1.767634 0.518351 -H 6.365079 -0.133224 -0.533034 -H 5.255789 0.923150 -1.389081 -H 5.504293 1.135755 0.351310 -H 4.140103 -2.613567 0.214988 -H 5.783262 -2.003292 0.275116 -H 4.665889 -1.644881 1.602460 -H 0.022604 -1.091753 1.084030 -H 2.307262 -1.879075 0.920974 -H 3.462774 1.870724 -0.861110 -H 1.165469 2.658467 -0.627261 -H -3.100673 2.322526 0.898873 -H -5.305838 1.177165 0.752078 -H -3.356658 -2.230234 -1.086826 -H -1.173708 -1.112751 -0.887177 diff --git a/phb.ljc b/phb.ljc new file mode 100644 index 0000000..597371e --- /dev/null +++ b/phb.ljc @@ -0,0 +1,33 @@ +* +2 +12 BENZENO + 1 6 0.000000 1.400000 0.000000 0.00000000 0.1100 3.7500 + 1 6 1.212436 0.700000 0.000000 0.00000000 0.1100 3.7500 + 1 6 1.212436 -0.700000 0.000000 0.00000000 0.1100 3.7500 + 1 6 0.000000 -1.400000 0.000000 0.00000000 0.1100 3.7500 + 1 6 -1.212436 -0.700000 0.000000 0.00000000 0.1100 3.7500 + 1 6 -1.212436 0.700000 0.000000 0.00000000 0.1100 3.7500 + 2 1 0.000000 2.488100 0.000000 0.00000000 0.0000 0.0000 + 2 1 2.154671 1.244000 0.000000 0.00000000 0.0000 0.0000 + 2 1 2.154671 -1.244000 0.000000 0.00000000 0.0000 0.0000 + 2 1 0.000000 -2.488100 0.000000 0.00000000 0.0000 0.0000 + 2 1 -2.154671 -1.244000 0.000000 0.00000000 0.0000 0.0000 + 2 1 -2.154671 1.244000 0.000000 0.00000000 0.0000 0.0000 +16 PLACEHOLDER + 1 6 0.672026 -2.823446 0.002631 -0.11500000 0.0700 3.5500 + 1 6 2.072026 -2.823446 0.002631 -0.11500000 0.0700 3.5500 + 1 6 2.768235 -1.617645 0.002633 -0.11500000 0.0700 3.5500 + 1 6 2.068235 -0.405210 0.002635 -0.11500000 0.0700 3.5500 + 1 14 0.675894 -0.405219 0.002635 0.15000000 0.0700 3.5500 + 1 10 -0.024198 -1.617602 0.002633 -0.11500000 0.0700 3.5500 + 2 1 0.132026 -3.758753 0.002629 0.11500000 0.0300 2.4200 + 2 1 2.612026 -3.758753 0.002629 0.11500000 0.0300 2.4200 + 2 1 2.608235 0.530098 0.002636 0.11500000 0.0300 2.4200 + 2 1 -1.104198 -1.617602 0.002633 0.11500000 0.0300 2.4200 + 3 8 -0.004114 0.772571 0.002637 -0.58500000 0.1700 3.0700 + 3 1 0.619778 1.502200 0.002638 0.43500000 0.0000 0.0000 + 4 6 4.278235 -1.617645 0.002633 0.11500000 0.1700 3.8000 + 4 1 4.634902 -0.743994 0.507036 0.00000000 0.0000 0.0000 + 4 1 4.634902 -2.491297 0.507036 0.00000000 0.0000 0.0000 + 4 1 4.634902 -1.617645 -1.006173 0.00000000 0.0000 0.0000 + diff --git a/run.log b/run.log deleted file mode 100644 index ade01f2..0000000 --- a/run.log +++ /dev/null @@ -1,212 +0,0 @@ -########################################################################################## -############# Welcome to DICEPLAYER version 1.0 ############# -########################################################################################## - -Your python version is 3.6.9 (default, Jan 26 2021, 15:33:00) -[GCC 8.4.0] - -Program started on Friday, 03 Dec 2021 at 21:20:24 - -Environment variables: -OMP_STACKSIZE = Not set - -========================================================================================== - CONTROL variables being used in this run: ------------------------------------------------------------------------------------------- - -altsteps = 20000 -freq = no -ghosts = no -initcyc = 1 -lps = no -maxcyc = 3 -maxstep = 0.3 -nprocs = 4 -opt = yes -qmprog = g16 -readhessian = no -switchcyc = 3 -tol_factor = 1.2 -vdwforces = no - ------------------------------------------------------------------------------------------- - DICE variables being used in this run: ------------------------------------------------------------------------------------------- - -combrule = * -dens = 0.75 -isave = 1000 -ljname = phb.ljc -ncores = 3 -nmol = 1 50 -nstep = 20000 20000 -outname = phb -press = 1.0 -progname = dice -randominit = first -temp = 300.0 -title = Diceplayer run -upbuf = 360 - ------------------------------------------------------------------------------------------- - GAUSSIAN variables being used in this run: ------------------------------------------------------------------------------------------- - -chgmult = 0 1 -level = MP2/aug-cc-pVTZ -pop = chelpg -qmprog = g16 - - -========================================================================================== - -Program dice found at /usr/local/bin/dice -Program g16 found at /usr/local/g16/g16 -Program formchk found at /usr/local/g16/formchk - -========================================================================================== - Potential parameters from file phb.ljc: ------------------------------------------------------------------------------------------- - -Combination rule: * -Types of molecules: 2 - -31 atoms in molecule type 1: ---------------------------------------------------------------------------------- -Lbl AN X Y Z Charge Epsilon Sigma Mass ---------------------------------------------------------------------------------- -1 6 -4.40034 0.66551 0.41431 -0.387593 0.07000 3.5500 12.0110 -1 6 -4.49957 -0.67279 -0.15327 0.804001 0.07000 3.5500 12.0110 -1 6 -3.25630 -1.28922 -0.61351 -0.338491 0.07000 3.5500 12.0110 -1 6 -2.05849 -0.66807 -0.48674 -0.220783 0.07000 3.5500 12.0110 -1 6 -1.97114 0.65148 0.10511 0.563527 0.07000 3.5500 12.0110 -1 6 -3.20601 1.29684 0.49712 -0.145372 0.07000 3.5500 12.0110 -2 8 -5.59453 -1.25309 -0.28765 -0.714461 0.21000 2.9600 15.9990 -3 7 -0.87404 1.35870 0.27636 -0.721792 0.17000 3.2500 14.0070 -4 6 0.38785 0.81512 0.14410 0.607123 0.07000 3.5500 12.0110 -4 6 1.40776 1.65800 -0.33261 -0.295476 0.07000 3.5500 12.0110 -4 6 0.75950 -0.46815 0.58791 -0.346661 0.07000 3.5500 12.0110 -4 6 2.71021 1.21519 -0.46119 -0.298546 0.07000 3.5500 12.0110 -4 6 2.07354 -0.89884 0.50831 -0.234284 0.07000 3.5500 12.0110 -4 6 3.08076 -0.08936 -0.06112 0.294651 0.07000 3.5500 12.0110 -5 7 4.37238 -0.53979 -0.20183 -0.190497 0.17000 3.2500 14.0070 -6 6 5.41980 0.43824 -0.43836 -0.213183 0.06600 3.5000 12.0110 -6 6 4.76361 -1.73522 0.52225 -0.225420 0.06600 3.5000 12.0110 -7 1 6.36700 -0.08536 -0.51962 0.109840 0.03000 2.5000 1.0079 -7 1 5.25196 0.96612 -1.37421 0.077140 0.03000 2.5000 1.0079 -7 1 5.49517 1.17412 0.36749 0.125758 0.03000 2.5000 1.0079 -7 1 4.15838 -2.58442 0.21446 0.083229 0.03000 2.5000 1.0079 -7 1 5.79695 -1.96254 0.28041 0.115361 0.03000 2.5000 1.0079 -7 1 4.67416 -1.61707 1.60661 0.121440 0.03000 2.5000 1.0079 -7 1 0.02813 -1.09554 1.08007 0.179434 0.03000 2.4200 1.0079 -7 1 2.31876 -1.86576 0.91913 0.177250 0.03000 2.4200 1.0079 -7 1 3.45102 1.89880 -0.84669 0.198844 0.03000 2.4200 1.0079 -7 1 1.14759 2.66909 -0.61498 0.173388 0.03000 2.4200 1.0079 -7 1 -3.11930 2.29679 0.90060 0.151534 0.03000 2.4200 1.0079 -7 1 -5.31582 1.13611 0.74479 0.205289 0.03000 2.4200 1.0079 -7 1 -3.33813 -2.25037 -1.10231 0.184367 0.03000 2.4200 1.0079 -7 1 -1.16374 -1.11791 -0.89380 0.160382 0.03000 2.4200 1.0079 - -16 atoms in molecule type 2: ---------------------------------------------------------------------------------- -Lbl AN X Y Z Charge Epsilon Sigma Mass ---------------------------------------------------------------------------------- -1 6 0.67203 -2.82345 0.00263 -0.115000 0.07000 3.5500 12.0110 -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 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 -2 1 -1.10420 -1.61760 0.00263 0.115000 0.03000 2.4200 1.0079 -3 8 -0.00411 0.77257 0.00264 -0.585000 0.17000 3.0700 15.9990 -3 1 0.61978 1.50220 0.00264 0.435000 0.00000 0.0000 1.0079 -4 6 4.27823 -1.61764 0.00263 0.115000 0.17000 3.8000 12.0110 -4 1 4.63490 -0.74399 0.50704 0.000000 0.00000 0.0000 1.0079 -4 1 4.63490 -2.49130 0.50704 0.000000 0.00000 0.0000 1.0079 -4 1 4.63490 -1.61764 -1.00617 0.000000 0.00000 0.0000 1.0079 - - -========================================================================================== - -Molecule type 1 - PHENOL_BLUE: - - Center of mass = ( -0.0000 , 0.0000 , 0.0000 ) - Moments of inertia = 3.144054E+02 2.801666E+03 3.027366E+03 - Major principal axis = ( 0.999972 , 0.007210 , 0.002184 ) - Inter principal axis = ( -0.007218 , 0.999967 , 0.003660 ) - Minor principal axis = ( -0.002157 , -0.003676 , 0.999991 ) - Characteristic lengths = ( 11.96 , 5.25 , 2.98 ) - Total mass = 226.28 au - Total charge = -0.0000 e - Dipole moment = ( 9.8367 , 0.6848 , 0.8358 ) Total = 9.8959 Debye - - Translating and rotating molecule to standard orientation... Done - - New values: - Center of mass = ( -0.0000 , -0.0000 , -0.0000 ) - Moments of inertia = 3.144054E+02 2.801666E+03 3.027366E+03 - 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 = ( 11.97 , 5.27 , 2.99 ) - Total mass = 226.28 au - Total charge = -0.0000 e - Dipole moment = ( 9.8432 , 0.6168 , 0.8120 ) Total = 9.8959 Debye - - -Molecule type 2 - PLACEHOLDER: - - 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 = 108.14 au - Total charge = -0.0000 e - Dipole moment = ( 2.3293 , 0.1593 , -0.0000 ) Total = 2.3347 Debye - - Translating and rotating molecule to standard orientation... Done - - New values: - Center of mass = ( 0.0000 , -0.0000 , 0.0000 ) - 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.65 , 4.95 , 1.51 ) - Total mass = 108.14 au - Total charge = -0.0000 e - Dipole moment = ( 1.9373 , 1.3031 , -0.0000 ) Total = 2.3347 Debye - -========================================================================================== - -Starting the iterative process. - ------------------------------------------------------------------------------------------- - Step # 1 ------------------------------------------------------------------------------------------- - -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 - -+----------------------------------------------------------------------------------------+ - -Calculation of forces initiated with Gaussian on 03 Dec 2021 at 21:22:39 diff --git a/run.log.backup b/run.log.backup deleted file mode 100644 index 384d5c1..0000000 --- a/run.log.backup +++ /dev/null @@ -1,200 +0,0 @@ -########################################################################################## -############# Welcome to DICEPLAYER version 1.0 ############# -########################################################################################## - -Your python version is 3.6.9 (default, Jan 26 2021, 15:33:00) -[GCC 8.4.0] - -Program started on Friday, 03 Dec 2021 at 21:18:49 - -Environment variables: -OMP_STACKSIZE = Not set - -========================================================================================== - CONTROL variables being used in this run: ------------------------------------------------------------------------------------------- - -altsteps = 20000 -freq = no -ghosts = no -initcyc = 1 -lps = no -maxcyc = 3 -maxstep = 0.3 -nprocs = 4 -opt = yes -qmprog = g16 -readhessian = no -switchcyc = 3 -tol_factor = 1.2 -vdwforces = no - ------------------------------------------------------------------------------------------- - DICE variables being used in this run: ------------------------------------------------------------------------------------------- - -combrule = * -dens = 0.75 -isave = 1000 -ljname = phb.ljc -ncores = 3 -nmol = 1 50 -nstep = 20000 20000 -outname = phb -press = 1.0 -progname = dice -randominit = first -temp = 300.0 -title = Diceplayer run -upbuf = 360 - ------------------------------------------------------------------------------------------- - GAUSSIAN variables being used in this run: ------------------------------------------------------------------------------------------- - -chgmult = 0 1 -level = MP2/aug-cc-pVTZ -pop = chelpg -qmprog = g16 - - -========================================================================================== - -Program dice found at /usr/local/bin/dice -Program g16 found at /usr/local/g16/g16 -Program formchk found at /usr/local/g16/formchk - -========================================================================================== - Potential parameters from file phb.ljc: ------------------------------------------------------------------------------------------- - -Combination rule: * -Types of molecules: 2 - -31 atoms in molecule type 1: ---------------------------------------------------------------------------------- -Lbl AN X Y Z Charge Epsilon Sigma Mass ---------------------------------------------------------------------------------- -1 6 -4.40034 0.66551 0.41431 -0.387593 0.07000 3.5500 12.0110 -1 6 -4.49957 -0.67279 -0.15327 0.804001 0.07000 3.5500 12.0110 -1 6 -3.25630 -1.28922 -0.61351 -0.338491 0.07000 3.5500 12.0110 -1 6 -2.05849 -0.66807 -0.48674 -0.220783 0.07000 3.5500 12.0110 -1 6 -1.97114 0.65148 0.10511 0.563527 0.07000 3.5500 12.0110 -1 6 -3.20601 1.29684 0.49712 -0.145372 0.07000 3.5500 12.0110 -2 8 -5.59453 -1.25309 -0.28765 -0.714461 0.21000 2.9600 15.9990 -3 7 -0.87404 1.35870 0.27636 -0.721792 0.17000 3.2500 14.0070 -4 6 0.38785 0.81512 0.14410 0.607123 0.07000 3.5500 12.0110 -4 6 1.40776 1.65800 -0.33261 -0.295476 0.07000 3.5500 12.0110 -4 6 0.75950 -0.46815 0.58791 -0.346661 0.07000 3.5500 12.0110 -4 6 2.71021 1.21519 -0.46119 -0.298546 0.07000 3.5500 12.0110 -4 6 2.07354 -0.89884 0.50831 -0.234284 0.07000 3.5500 12.0110 -4 6 3.08076 -0.08936 -0.06112 0.294651 0.07000 3.5500 12.0110 -5 7 4.37238 -0.53979 -0.20183 -0.190497 0.17000 3.2500 14.0070 -6 6 5.41980 0.43824 -0.43836 -0.213183 0.06600 3.5000 12.0110 -6 6 4.76361 -1.73522 0.52225 -0.225420 0.06600 3.5000 12.0110 -7 1 6.36700 -0.08536 -0.51962 0.109840 0.03000 2.5000 1.0079 -7 1 5.25196 0.96612 -1.37421 0.077140 0.03000 2.5000 1.0079 -7 1 5.49517 1.17412 0.36749 0.125758 0.03000 2.5000 1.0079 -7 1 4.15838 -2.58442 0.21446 0.083229 0.03000 2.5000 1.0079 -7 1 5.79695 -1.96254 0.28041 0.115361 0.03000 2.5000 1.0079 -7 1 4.67416 -1.61707 1.60661 0.121440 0.03000 2.5000 1.0079 -7 1 0.02813 -1.09554 1.08007 0.179434 0.03000 2.4200 1.0079 -7 1 2.31876 -1.86576 0.91913 0.177250 0.03000 2.4200 1.0079 -7 1 3.45102 1.89880 -0.84669 0.198844 0.03000 2.4200 1.0079 -7 1 1.14759 2.66909 -0.61498 0.173388 0.03000 2.4200 1.0079 -7 1 -3.11930 2.29679 0.90060 0.151534 0.03000 2.4200 1.0079 -7 1 -5.31582 1.13611 0.74479 0.205289 0.03000 2.4200 1.0079 -7 1 -3.33813 -2.25037 -1.10231 0.184367 0.03000 2.4200 1.0079 -7 1 -1.16374 -1.11791 -0.89380 0.160382 0.03000 2.4200 1.0079 - -16 atoms in molecule type 2: ---------------------------------------------------------------------------------- -Lbl AN X Y Z Charge Epsilon Sigma Mass ---------------------------------------------------------------------------------- -1 6 0.67203 -2.82345 0.00263 -0.115000 0.07000 3.5500 12.0110 -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 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 -2 1 -1.10420 -1.61760 0.00263 0.115000 0.03000 2.4200 1.0079 -3 8 -0.00411 0.77257 0.00264 -0.585000 0.17000 3.0700 15.9990 -3 1 0.61978 1.50220 0.00264 0.435000 0.00000 0.0000 1.0079 -4 6 4.27823 -1.61764 0.00263 0.115000 0.17000 3.8000 12.0110 -4 1 4.63490 -0.74399 0.50704 0.000000 0.00000 0.0000 1.0079 -4 1 4.63490 -2.49130 0.50704 0.000000 0.00000 0.0000 1.0079 -4 1 4.63490 -1.61764 -1.00617 0.000000 0.00000 0.0000 1.0079 - - -========================================================================================== - -Molecule type 1 - PHENOL_BLUE: - - Center of mass = ( -0.0000 , 0.0000 , 0.0000 ) - Moments of inertia = 3.144054E+02 2.801666E+03 3.027366E+03 - Major principal axis = ( 0.999972 , 0.007210 , 0.002184 ) - Inter principal axis = ( -0.007218 , 0.999967 , 0.003660 ) - Minor principal axis = ( -0.002157 , -0.003676 , 0.999991 ) - Characteristic lengths = ( 11.96 , 5.25 , 2.98 ) - Total mass = 226.28 au - Total charge = -0.0000 e - Dipole moment = ( 9.8367 , 0.6848 , 0.8358 ) Total = 9.8959 Debye - - Translating and rotating molecule to standard orientation... Done - - New values: - Center of mass = ( -0.0000 , -0.0000 , -0.0000 ) - Moments of inertia = 3.144054E+02 2.801666E+03 3.027366E+03 - 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 = ( 11.97 , 5.27 , 2.99 ) - Total mass = 226.28 au - Total charge = -0.0000 e - Dipole moment = ( 9.8432 , 0.6168 , 0.8120 ) Total = 9.8959 Debye - - -Molecule type 2 - PLACEHOLDER: - - 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 = 108.14 au - Total charge = -0.0000 e - Dipole moment = ( 2.3293 , 0.1593 , -0.0000 ) Total = 2.3347 Debye - - Translating and rotating molecule to standard orientation... Done - - New values: - Center of mass = ( 0.0000 , -0.0000 , 0.0000 ) - 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.65 , 4.95 , 1.51 ) - Total mass = 108.14 au - Total charge = -0.0000 e - Dipole moment = ( 1.9373 , 1.3031 , -0.0000 ) Total = 2.3347 Debye - -========================================================================================== - -Starting the iterative process. - ------------------------------------------------------------------------------------------- - Step # 1 ------------------------------------------------------------------------------------------- - -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