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

Merge pull request #1 from HideyoshiNakazone/devel

Browse Source 
Translation of Gaussian Processes and Step Calculations Fixes
This commit is contained in:
Vitor Hideyoshi Nakazone Batista
2022-01-08 13:22:59 -03:00
committed by GitHub
parent 2a4e9eff0c 926ffc5c6b
commit 4d122a43aa
12 changed files with 515 additions and 3218 deletions
Download Patch File Download Diff File Expand all files Collapse all files

5
control.in
View File

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

527
diceplayer/DPpack/.Dice.py
View File

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

600
diceplayer/DPpack/.MolHandling.py
View File

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

813
diceplayer/DPpack/.SetGlobals.py
View File

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

373
diceplayer/DPpack/Gaussian.py
View File

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

619
diceplayer/DPpack/MolHandling.py
View File

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

207
diceplayer/DPpack/SetGlobals.py
View File

@@ -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):
def calculate_step(self, cycle: int, gradient: np.ndarray,
hessian: np.ndarray) -> np.ndarray:
invhessian = linalg.inv(self.system.molecule[0].hessian)
pre_step = -1 * np.matmul(invhessian, self.system.molecule[0].gradient.T).T
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"
@@ -1147,9 +1138,10 @@ class Internal:
fh.write(fstr.format(atom.lbl, atom.na, atom.rx, atom.ry,
atom.rz, atom.chg, atom.eps, atom.sig))
#----------------------------------------------------------------------------------------------------------------------------------------
# Gaussian related methods
def read_forces_fchk(self, file, fh):
#----------------------------------------------------------------------------------------------------------------------------------------
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)
@@ -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")

99
diceplayer/__main__.py
View File

@@ -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"
@@ -197,8 +199,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"
# ## 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)
try:
gradient
old_gradient = gradient
except:
pass
# ## 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")
gradient = internal.read_forces_fchk(file, internal.outfile)
# ## Save gradient and hessian
# internal.gaussian.print_grad_hessian(cycle, cur_gradient, hessian)
print(type(gradient),"\n",gradient.shape,"\n",gradient)
# ## Calculate the step and update the position
# step = internal.calculate_step(cur_gradient, hessian, internal.outfile)
# position += step
# 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, gradient, old_gradient, hessian)
internal.outfile.write("Done\n")
# 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
# ## Update the geometry of the reference molecule
# internal.system.molecule[0].update_molecule(position, internal.outfile)
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":
internal.system.print_geom(cycle, geomsfh)
geomsfh.write(40 * "-"+"\n")
internal.outfile.write("\n+" + 88 * "-" + "+\n")
####
#### End of the iterative process
####
# ####
# #### End of the iterative process
# ####
## imprimir ultimas mensagens, criar um arquivo de potencial para ser usado em eventual
## continuacao, fechar arquivos (geoms.xyz, run.log, ...)
# # ## imprimir ultimas mensagens, criar um arquivo de potencial para ser usado em eventual
# # ## continuacao, fechar arquivos (geoms.xyz, run.log, ...)
# # internal.outfile.write("\nDiceplayer finished normally!\n")
# # internal.outfile.close()
# # ####
# # #### End of the program
# # ####
internal.outfile.write("\nDiceplayer finished normally!\n")
internal.outfile.close()
####
#### End of the program
####

33
geoms.xyz
View File

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

33
phb.ljc Normal file
View File

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

212
run.log
View File

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

200
run.log.backup
View File

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