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

Initial Translation of Gaussian Processes and Packaging of DicePlayer python module

Browse Source 
This commit adds the methods that were present in the Gaussian.py file into the SetGlobals.py file and packages the program into a diceplayer module so it can be ran using 'python3 -m diceplayer'

Signed-off-by: Vitor Hideyoshi <vitor.h.n.batista@gmail.com>
This commit is contained in:
Vitor Hideyoshi
2021-12-03 21:28:27 +00:00
parent 49d509029f
commit 2a4e9eff0c
15 changed files with 1088 additions and 337 deletions
Download Patch File Download Diff File Expand all files Collapse all files

527
diceplayer/DPpack/.Dice.py Normal file
View File

@@ -0,0 +1,527 @@
import sys, os, time
import subprocess
from copy import deepcopy
from numpy import random
from DPpack.PTable import *
from DPpack.SetGlobals import *
from DPpack.Misc import *
####################################### functions ######################################
def make_inputs(cycle, proc):
step_dir = "step{:02d}".format(cycle)
proc_dir = "p{:02d}".format(proc)
path = step_dir + os.sep + proc_dir
num = time.time() ## Take the decimal places 7 to 12 of the
num = (num - int(num)) * 1e6 ## time in seconds as a floating point
num = int((num - int(num)) * 1e6) ## to make an integer in the range 1-1e6
random.seed( (os.getpid() * num) % (max_seed + 1) )
if not dice['randominit']:
xyzfile = dice['outname'] + ".xyz.last-" + "p{:02d}".format(proc)
make_init_file(path, xyzfile)
if len(dice['nstep']) == 2: ## Means NVT simulation
make_nvt_ter(path)
make_nvt_eq(path)
elif len(dice['nstep']) == 3: ## Means NPT simulation
if dice['randominit']:
make_nvt_ter(path)
else:
dice['dens'] = new_density(proc)
make_npt_ter(path)
make_npt_eq(path)
else:
sys.exit("Error: bad number of entries for 'nstep'")
make_potential(path)
return
def make_nvt_ter(path):
file = path + os.sep + "NVT.ter"
try:
fh = open(file, "w")
except:
sys.exit("Error: cannot open file {}".format(file))
fh.write("title = {} - NVT Thermalization\n".format(dice['title']))
fh.write("ncores = {}\n".format(dice['ncores']))
fh.write("ljname = {}\n".format(dice['ljname']))
fh.write("outname = {}\n".format(dice['outname']))
string = " ".join(str(x) for x in dice['nmol'])
fh.write("nmol = {}\n".format(string))
fh.write("dens = {}\n".format(dice['dens']))
fh.write("temp = {}\n".format(dice['temp']))
if dice['randominit']:
fh.write("init = yes\n")
fh.write("nstep = {}\n".format(dice['nstep'][0]))
else:
fh.write("init = yesreadxyz\n")
fh.write("nstep = {}\n".format(player['altsteps']))
fh.write("vstep = 0\n")
fh.write("mstop = 1\n")
fh.write("accum = no\n")
fh.write("iprint = 1\n")
fh.write("isave = 0\n")
fh.write("irdf = 0\n")
seed = int(1e6 * random.random())
fh.write("seed = {}\n".format(seed))
fh.close()
return
def make_nvt_eq(path):
file = path + os.sep + "NVT.eq"
try:
fh = open(file, "w")
except:
sys.exit("Error: cannot open file {}".format(file))
fh.write("title = {} - NVT Production\n".format(dice['title']))
fh.write("ncores = {}\n".format(dice['ncores']))
fh.write("ljname = {}\n".format(dice['ljname']))
fh.write("outname = {}\n".format(dice['outname']))
string = " ".join(str(x) for x in dice['nmol'])
fh.write("nmol = {}\n".format(string))
fh.write("dens = {}\n".format(dice['dens']))
fh.write("temp = {}\n".format(dice['temp']))
fh.write("init = no\n")
fh.write("nstep = {}\n".format(dice['nstep'][1]))
fh.write("vstep = 0\n")
fh.write("mstop = 1\n")
fh.write("accum = no\n")
fh.write("iprint = 1\n")
fh.write("isave = {}\n".format(dice['isave']))
fh.write("irdf = {}\n".format(10 * player['nprocs']))
seed = int(1e6 * random.random())
fh.write("seed = {}\n".format(seed))
fh.close()
return
def make_npt_ter(path):
file = path + os.sep + "NPT.ter"
try:
fh = open(file, "w")
except:
sys.exit("Error: cannot open file {}".format(file))
fh.write("title = {} - NPT Thermalization\n".format(dice['title']))
fh.write("ncores = {}\n".format(dice['ncores']))
fh.write("ljname = {}\n".format(dice['ljname']))
fh.write("outname = {}\n".format(dice['outname']))
string = " ".join(str(x) for x in dice['nmol'])
fh.write("nmol = {}\n".format(string))
fh.write("press = {}\n".format(dice['press']))
fh.write("temp = {}\n".format(dice['temp']))
if dice['randominit']:
fh.write("init = no\n") ## Because there will be a previous NVT simulation
fh.write("vstep = {}\n".format(int(dice['nstep'][1] / 5)))
else:
fh.write("init = yesreadxyz\n")
fh.write("dens = {:<8.4f}\n".format(dice['dens']))
fh.write("vstep = {}\n".format(int(player['altsteps'] / 5)))
fh.write("nstep = 5\n")
fh.write("mstop = 1\n")
fh.write("accum = no\n")
fh.write("iprint = 1\n")
fh.write("isave = 0\n")
fh.write("irdf = 0\n")
seed = int(1e6 * random.random())
fh.write("seed = {}\n".format(seed))
fh.close()
return
def make_npt_eq(path):
file = path + os.sep + "NPT.eq"
try:
fh = open(file, "w")
except:
sys.exit("Error: cannot open file {}".format(file))
fh.write("title = {} - NPT Production\n".format(dice['title']))
fh.write("ncores = {}\n".format(dice['ncores']))
fh.write("ljname = {}\n".format(dice['ljname']))
fh.write("outname = {}\n".format(dice['outname']))
string = " ".join(str(x) for x in dice['nmol'])
fh.write("nmol = {}\n".format(string))
fh.write("press = {}\n".format(dice['press']))
fh.write("temp = {}\n".format(dice['temp']))
fh.write("nstep = 5\n")
fh.write("vstep = {}\n".format(int(dice['nstep'][2] / 5)))
fh.write("init = no\n")
fh.write("mstop = 1\n")
fh.write("accum = no\n")
fh.write("iprint = 1\n")
fh.write("isave = {}\n".format(dice['isave']))
fh.write("irdf = {}\n".format(10 * player['nprocs']))
seed = int(1e6 * random.random())
fh.write("seed = {}\n".format(seed))
fh.close()
return
def make_init_file(path, file):
if not os.path.isfile(file):
sys.exit("Error: cannot find the xyz file {} in main directory".format(file))
try:
with open(file) as fh:
xyzfile = fh.readlines()
except:
sys.exit("Error: cannot open file {}".format(file))
nsites_mm = 0
for i in range(1, len(dice['nmol'])):
nsites_mm += dice['nmol'][i] * len(molecules[i])
nsites_mm *= -1 ## Become an index to count from the end of xyzfile (list)
xyzfile = xyzfile[nsites_mm :] ## Only the MM atoms of the last configuration remains
file = path + os.sep + dice['outname'] + ".xy"
try:
fh = open(file, "w")
except:
sys.exit("Error: cannot open file {}".format(file))
for atom in molecules[0]:
fh.write("{:>10.6f} {:>10.6f} {:>10.6f}\n".format(atom['rx'], atom['ry'],
atom['rz']))
for ghost in ghost_atoms:
fh.write("{:>10.6f} {:>10.6f} {:>10.6f}\n".format(ghost['rx'], ghost['ry'],
ghost['rz']))
for lps in lp_atoms:
fh.write("{:>10.6f} {:>10.6f} {:>10.6f}\n".format(lps['rx'], lps['ry'],
lps['rz']))
for line in xyzfile:
atom = line.split()
rx = float(atom[1])
ry = float(atom[2])
rz = float(atom[3])
fh.write("{:>10.5f} {:>10.5f} {:>10.5f}\n".format(rx, ry, rz))
fh.write("$end")
fh.close()
return
def make_potential(path):
fstr = "{:<3d} {:>3d} {:>10.5f} {:>10.5f} {:>10.5f} {:>10.6f} {:>9.5f} {:>7.4f}\n"
file = path + os.sep + dice['ljname']
try:
fh = open(file, "w")
except:
sys.exit("Error: cannot open file {}".format(file))
fh.write("{}\n".format(dice['combrule']))
fh.write("{}\n".format(len(dice['nmol'])))
nsites_qm = len(molecules[0]) + len(ghost_atoms) + len(lp_atoms)
## Print the sites of the QM molecule
fh.write("{}\n".format(nsites_qm))
for atom in molecules[0]:
fh.write(fstr.format(atom['lbl'], atom['na'], atom['rx'], atom['ry'], atom['rz'],
atom['chg'], atom['eps'], atom['sig']))
ghost_label = molecules[0][-1]['lbl'] + 1
for ghost in ghost_atoms:
fh.write(fstr.format(ghost_label, ghost_number, ghost['rx'], ghost['ry'],
ghost['rz'], ghost['chg'], 0, 0))
ghost_label += 1
for lp in lp_atoms:
fh.write(fstr.format(ghost_label, ghost_number, lp['rx'], lp['ry'], lp['rz'],
lp['chg'], 0, 0))
## Print the sites of the other molecules
for mol in molecules[1:]:
fh.write("{}\n".format(len(mol)))
for atom in mol:
fh.write(fstr.format(atom['lbl'], atom['na'], atom['rx'], atom['ry'],
atom['rz'], atom['chg'], atom['eps'], atom['sig']))
return
def make_proc_dir(cycle, proc):
step_dir = "step{:02d}".format(cycle)
proc_dir = "p{:02d}".format(proc)
path = step_dir + os.sep + proc_dir
try:
os.makedirs(path)
except:
sys.exit("Error: cannot make directory {}".format(path))
return
def run_dice(cycle, proc, fh):
step_dir = "step{:02d}".format(cycle)
proc_dir = "p{:02d}".format(proc)
path = step_dir + os.sep + proc_dir
working_dir = os.getcwd()
os.chdir(path)
fh.write("Simulation process {} initiated with pid {}\n".format(proc_dir, os.getpid()))
if len(dice['nstep']) == 2: ## Means NVT simulation
## NVT thermalization
string = "(from " + ("random" if dice['randominit'] else "previous") + " configuration)"
fh.write("p{:02d}> NVT thermalization initiated {} on {}\n".format(proc, string,
date_time()))
infh = open("NVT.ter")
outfh = open("NVT.ter.out", "w")
exit_status = subprocess.call(dice['progname'], stdin=infh, stdout=outfh)
infh.close()
outfh.close()
if os.getppid() == 1: ## Parent process is dead
sys.exit()
if exit_status != 0:
sys.exit("Dice process p{:02d} did not exit properly".format(proc))
else:
outfh = open("NVT.ter.out") ## Open again to seek the normal end flag
flag = outfh.readlines()[dice_flag_line].strip()
outfh.close()
if flag != dice_end_flag:
sys.exit("Dice process p{:02d} did not exit properly".format(proc))
## NVT production
fh.write("p{:02d}> NVT production initiated on {}\n".format(proc, date_time()))
infh = open("NVT.eq")
outfh = open("NVT.eq.out", "w")
exit_status = subprocess.call(dice['progname'], stdin=infh, stdout=outfh)
infh.close()
outfh.close()
if os.getppid() == 1: ## Parent process is dead
sys.exit()
if exit_status != 0:
sys.exit("Dice process p{:02d} did not exit properly".format(proc))
else:
outfh = open("NVT.eq.out") ## Open again to seek the normal end flag
flag = outfh.readlines()[dice_flag_line].strip()
outfh.close()
if flag != dice_end_flag:
sys.exit("Dice process p{:02d} did not exit properly".format(proc))
fh.write("p{:02d}> ----- NVT production finished on {}\n".format(proc,
date_time()))
elif len(dice['nstep']) == 3: ## Means NPT simulation
## NVT thermalization if randominit
if dice['randominit']:
string = "(from random configuration)"
fh.write("p{:02d}> NVT thermalization initiated {} on {}\n".format(proc,
string, date_time()))
infh = open("NVT.ter")
outfh = open("NVT.ter.out", "w")
exit_status = subprocess.call(dice['progname'], stdin=infh, stdout=outfh)
infh.close()
outfh.close()
if os.getppid() == 1: ## Parent process is dead
sys.exit()
if exit_status != 0:
sys.exit("Dice process p{:02d} did not exit properly".format(proc))
else:
outfh = open("NVT.ter.out") ## Open again to seek the normal end flag
flag = outfh.readlines()[dice_flag_line].strip()
outfh.close()
if flag != dice_end_flag:
sys.exit("Dice process p{:02d} did not exit properly".format(proc))
## NPT thermalization
string = (" (from previous configuration) " if not dice['randominit'] else " ")
fh.write("p{:02d}> NPT thermalization initiated{}on {}\n".format(proc, string,
date_time()))
infh = open("NPT.ter")
outfh = open("NPT.ter.out", "w")
exit_status = subprocess.call(dice['progname'], stdin=infh, stdout=outfh)
infh.close()
outfh.close()
if os.getppid() == 1: ## Parent process is dead
sys.exit()
if exit_status != 0:
sys.exit("Dice process p{:02d} did not exit properly".format(proc))
else:
outfh = open("NPT.ter.out") ## Open again to seek the normal end flag
flag = outfh.readlines()[dice_flag_line].strip()
outfh.close()
if flag != dice_end_flag:
sys.exit("Dice process p{:02d} did not exit properly".format(proc))
## NPT production
fh.write("p{:02d}> NPT production initiated on {}\n".format(proc, date_time()))
infh = open("NPT.eq")
outfh = open("NPT.eq.out", "w")
exit_status = subprocess.call(dice['progname'], stdin=infh, stdout=outfh)
infh.close()
outfh.close()
if os.getppid() == 1: ## Parent process is dead
sys.exit()
if exit_status != 0:
sys.exit("Dice process p{:02d} did not exit properly".format(proc))
else:
outfh = open("NPT.eq.out") ## Open again to seek the normal end flag
flag = outfh.readlines()[dice_flag_line].strip()
outfh.close()
if flag != dice_end_flag:
sys.exit("Dice process p{:02d} did not exit properly".format(proc))
fh.write("p{:02d}> ----- NPT production finished on {}\n".format(proc,
date_time()))
os.chdir(working_dir)
return
def print_last_config(cycle, proc):
step_dir = "step{:02d}".format(cycle)
proc_dir = "p{:02d}".format(proc)
path = step_dir + os.sep + proc_dir
file = path + os.sep + dice['outname'] + ".xyz"
if not os.path.isfile(file):
sys.exit("Error: cannot find the xyz file {}".format(file))
try:
with open(file) as fh:
xyzfile = fh.readlines()
except:
sys.exit("Error: cannot open file {}".format(file))
nsites = ( len(molecules[0]) + len(ghost_atoms) + len(lp_atoms) ) * dice['nmol'][0]
for i in range(1, len(dice['nmol'])):
nsites += dice['nmol'][i] * len(molecules[i])
nsites += 2 ## To include the comment line and the number of atoms (xyz file format)
nsites *= -1 ## Become an index to count from the end of xyzfile (list)
xyzfile = xyzfile[nsites :] ## Take the last configuration
file = dice['outname'] + ".xyz.last-" + proc_dir
fh = open(file, "w")
for line in xyzfile:
fh.write(line)
fh.close()
return
def new_density(proc):
file = dice['outname'] + ".xyz.last-" + "p{:02d}".format(proc)
if not os.path.isfile(file):
sys.exit("Error: cannot find the xyz file {} in main directory".format(file))
try:
with open(file) as fh:
xyzfile = fh.readlines()
except:
sys.exit("Error: cannot open file {}".format(file))
box = xyzfile[1].split()
volume = float(box[-3]) * float(box[-2]) * float(box[-1])
total_mass = 0
for i in range(len(molecules)):
mol_mass = 0
for atom in molecules[i]:
mol_mass += atom['mass']
total_mass += mol_mass * dice['nmol'][i]
density = (total_mass / volume) * umaAng3_to_gcm3
return density
def simulation_process(cycle, proc, logfh):
try:
make_proc_dir(cycle, proc)
make_inputs(cycle, proc)
run_dice(cycle, proc, logfh)
except Exception as err:
sys.exit(err)
return

600
diceplayer/DPpack/.MolHandling.py Normal file
View File

@@ -0,0 +1,600 @@
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 Normal file
View File

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

373
diceplayer/DPpack/Gaussian.py Normal file
View File

@@ -0,0 +1,373 @@
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

67
diceplayer/DPpack/Misc.py Normal file
View File

@@ -0,0 +1,67 @@
import os, sys, time
from posixpath import sep
import shutil, gzip
####################################### functions ######################################
def weekday_date_time():
return time.strftime("%A, %d %b %Y at %H:%M:%S")
def date_time():
return time.strftime("%d %b %Y at %H:%M:%S")
def compress_files_1mb(path):
working_dir = os.getcwd()
os.chdir(path)
files = filter(os.path.isfile, os.listdir(os.curdir))
for file in files:
if os.path.getsize(file) > 1024 * 1024: ## If bigger than 1MB
filegz = file + ".gz"
try:
with open(file, 'rb') as f_in:
with gzip.open(filegz, 'wb') as f_out:
shutil.copyfileobj(f_in, f_out)
except:
sys.exit("Error: cannot compress file {}".format(file))
os.chdir(working_dir)
return
def make_simulation_dir():
sim_dir = "simfiles"
if os.path.exists(sim_dir):
sys.exit("Error: a file or a directory {} already exists, move or delete de simfiles directory to continue.".format(sim_dir))
try:
os.makedirs(sim_dir)
except:
sys.exit("Error: cannot make directory {}".format(sim_dir))
def make_step_dir(cycle):
sim_dir = "simfiles"
step_dir = "step{:02d}".format(cycle)
path = sim_dir + os.sep + step_dir
if os.path.exists(path):
sys.exit("Error: a file or directory {} already exists".format(step_dir))
try:
os.makedirs(path)
except:
sys.exit("Error: cannot make directory {}".format(step_dir))
def make_qm_dir(cycle):
sim_dir = "simfiles"
step_dir = "step{:02d}".format(cycle)
path = sim_dir + os.sep + step_dir + os.sep + "qm"
try:
os.makedirs(path)
except:
sys.exit("Error: cannot make directory {}".format(path))

449
diceplayer/DPpack/MolHandling.py Normal file
View File

@@ -0,0 +1,449 @@
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 *
env = ["OMP_STACKSIZE"]
bohr2ang = 0.52917721092
ang2bohr = 1/bohr2ang
# Usaremos uma nova classe que ira conter toda interação entre moleculas
class System:
def __init__(self):
self.molecule = []
self.nmols = []
def add_type(self, nmols, m):
self.molecule.append(m)
self.nmols.append(nmols)
# 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):
com1 = self.molecule[a].center_of_mass()
com2 = self.molecule[b].center_of_mass()
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 rmsd_fit(self, index_p, index_r):
projecting_mol = self.molecule[index_p]
reference_mol = self.molecule[index_r]
if len(projecting_mol.atom) != len(reference_mol.atom):
sys.exit("Error in RMSD fit procedure: molecules have different number of atoms")
dim = len(projecting_mol.atom)
new_projecting_mol = deepcopy(projecting_mol)
new_reference_mol = deepcopy(reference_mol)
new_projecting_mol.center_of_mass_to_origin()
new_reference_mol.center_of_mass_to_origin()
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.atom[i].rx = x[i,0]
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)
return rmsd, projected_mol
def update_molecule(self, position, fh):
position_in_ang = (position * bohr2ang).tolist()
self.add_molecule(deepcopy(self.molecule[0]))
for atom in self.molecule[-1].atom:
atom.rx = position_in_ang.pop(0)
atom.ry = position_in_ang.pop(0)
atom.rz = position_in_ang.pop(0)
rmsd, self.molecule[0] = self.rmsd_fit(-1, 0)
self.molecule.pop(-1)
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):
if criterium in None:
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]
self.add_molecule(self.molecule[index_m].translate(tr_vector))
if criterium == "com":
dist = self.center_of_mass_distance(index_r, -1)
else:
dist = self.minimum_distance(index_r, -1)
if dist < min_dist:
min_dist = dist
nearestmol = deepcopy(self.molecule[-1])
self.molecule.pop(-1)
return min_dist, nearestmol
def print_geom(self, cycle, fh):
fh.write("{}\n".format(len(self.molecule[0].atom)))
fh.write("Cycle # {}\n".format(cycle))
for atom in self.molecule[0].atom:
symbol = atomsymb[atom.na]
fh.write("{:<2s} {:>10.6f} {:>10.6f} {:>10.6f}\n".format(symbol,
atom.rx, atom.ry, atom.rz))
# 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

348
diceplayer/DPpack/Molcas.py Normal file
View File

@@ -0,0 +1,348 @@
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_log(file, fh):
forces = []
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("Molecular gradients") < 0: ## expression not found
start = logfile.pop(0).strip()
logfile = logfile[7:] ## skip next 7 lines
for i in range(len(molecules[0])):
values = logfile.pop(0).split()
values = values[1:]
forces.extend([ float(x) for x in values ])
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_log(file):
force_const = []
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("Force constant matrix") < 0:
start = logfile.pop(0).strip()
logfile = logfile[1:] ## skip next 1 line
degrees = 3 * len(molecules[0])
dim = degrees
last = round(dim * dim)
count = 0
while True:
values = logfile.pop(0).rstrip()
while len(values) != 0:
new_value = values[:16]
values = values[16:]
force_const.append(float(new_value))
count += 1
if count >= last:
break
hessian = np.array(force_const).reshape(dim, dim)
hessian = hessian[:degrees, :degrees] ## remove degrees related to ghost atoms
for i in range(degrees):
for j in range(i + 1):
hessian[j,i] = hessian[i,j] ## force the hessian to be symmetric
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
def make_asec_file(cycle, asec_charges):
path = "step{:02d}".format(cycle) + os.sep + "qm"
file = path + os.sep + "asec.xfield"
try:
fh = open(file, "w")
except:
sys.exit("Error: cannot open file {}".format(file))
fh.write("{} Angstrom\n".format(len(asec_charges)))
## Write the ASEC:
for charge in asec_charges:
fh.write("{:>10.5f} {:>10.5f} {:>10.5f} {:>11.8f} 0.0 0.0 0.0\n".format(
charge['rx'], charge['ry'], charge['rz'], charge['chg']))
fh.write("End of input\n")
fh.close()
return
def make_force_input(cycle, asec_charges):
path = "step{:02d}".format(cycle) + os.sep + "qm"
file = path + os.sep + "asec.input"
try:
fh = open(file, "w")
except:
sys.exit("Error: cannot open file {}".format(file))
fh.write(" &Gateway\n")
fh.write(" Coord\n")
nsites = len(molecules[0])
if cycle >= player['switchcyc']:
nsites += len(ghost_atoms) + len(lp_atoms)
fh.write("{}\n".format(nsites))
fh.write("\nForce calculation - Cycle number {}\n".format(cycle))
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("basis = {}\n".format(molcas['basis']))
fh.write("group= nosym\n")
fh.write(" XFIELD\n")
fh.write(">>> Include asec.xfield\n")
if not os.path.isfile(molcas['mbottom']):
sys.exit("Error: cannot find file {} in main directory".format(molcas['mbottom']))
try:
with open(molcas['mbottom']) as mbottomfile:
mbottom = mbottomfile.readlines()
except:
sys.exit("Error: cannot open file {}".format(molcas['mbottom']))
for line in mbottom:
fh.write(line)
fh.write(" &Alaska\nPNEW\n &SLAPAF\nCartesian\n")
fh.close()
return
def make_charge_input(cycle, asec_charges):
path = "step{:02d}".format(cycle) + os.sep + "qm"
file = path + os.sep + "asec.input"
try:
fh = open(file, "w")
except:
sys.exit("Error: cannot open file {}".format(file))
fh.write(" &Gateway\n")
fh.write(" Coord\n")
nsites = len(molecules[0])
if cycle >= player['switchcyc']:
nsites += len(ghost_atoms) + len(lp_atoms)
fh.write("{}\n".format(nsites))
fh.write("\nForce calculation - Cycle number {}\n".format(cycle))
for atom in molecules[0]:
symbol = atomsymb[atom['na']]
fh.write("{:<2s} {:>10.5f} {:>10.5f} {:>10.5f}\n".format(symbol,
atom['rx'], atom['ry'], atom['rz']))
for 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("basis = {}\n".format(molcas['basis']))
fh.write("group= nosym\n")
fh.write(" XFIELD\n")
fh.write(">>> Include asec.xfield\n")
if not os.path.isfile(molcas['mbottom']):
sys.exit("Error: cannot find file {} in main directory".format(molcas['mbottom']))
try:
with open(molcas['mbottom']) as mbottomfile:
mbottom = mbottomfile.readlines()
except:
sys.exit("Error: cannot open file {}".format(molcas['mbottom']))
for line in mbottom:
fh.write(line)
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
def run_gaussian(cycle, type, fh):
path = "step{:02d}".format(cycle) + os.sep + "qm"
work_dir = os.getcwd()
os.chdir(path)
if type == "force":
infile = "asec.gjf"
elif type == "charge":
infile = "asec2.gjf"
fh.write("\nCalculation of {}s initiated with Gaussian on {}\n".format(type, date_time()))
exit_status = subprocess.call([player['qmprog'], infile])
if exit_status != 0:
sys.exit("Gaussian process did not exit properly")
fh.write("Calculation of {}s finished on {}\n".format(type, date_time()))
os.chdir(work_dir)
return
def run_formchk(cycle, fh):
path = "step{:02d}".format(cycle) + os.sep + "qm"
work_dir = os.getcwd()
os.chdir(path)
fh.write("Formatting the checkpoint file... ")
exit_status = subprocess.call(["formchk", "asec.chk"])
fh.write("Done\n")
os.chdir(work_dir)
return

266
diceplayer/DPpack/Optimization.py Normal file
View File

@@ -0,0 +1,266 @@
import sys, math
from copy import deepcopy
import numpy as np
from numpy import linalg
from DPpack.SetGlobals import *
epsilon = 1e-8
####################################### functions ######################################
def best_previous_point():
min_energy = 0
idx = 0
for energy in internal['energy'][:-1]:
if energy < min_energy or abs(energy - min_energy) < 1e-10:
min_energy = energy
min_idx = idx
idx += 1
return min_idx
def best_point():
min_energy = 0
idx = 0
for energy in internal['energy']:
if energy < min_energy or abs(energy - min_energy) < 1e-10:
min_energy = energy
min_idx = idx
idx += 1
return min_idx
def line_search(fh):
X1 = internal['position'][-1] # numpy array
e1 = internal['energy'][-1]
G1 = internal['gradient'][-1] # numpy array
idx = best_previous_point()
X0 = internal['position'][idx] # numpy array
e0 = internal['energy'][idx]
G0 = internal['gradient'][idx] # numpy array
# First try a quartic fit
fh.write("Attempting a quartic fit.\n")
success, y0 = quartic_fit(X0, X1, e0, e1, G0, G1, fh)
if success and y0 > 0:
if y0 < 1:
new_point = X0 + y0*(X1 - X0)
new_gradient = interpolate_gradient(G0, G1, y0)
new_gradient = perpendicular_projection(new_gradient, X1 - X0)
fh.write("Line search succeded.\n")
return True, new_point, new_gradient
else:
idx = best_point()
if idx == len(internal['energy']) - 1:
new_point = X0 + y0*(X1 - X0)
new_gradient = interpolate_gradient(G0, G1, y0)
new_gradient = perpendicular_projection(new_gradient, X1 - X0)
fh.write("Line search succeded.\n")
return True, new_point, new_gradient
else:
fh.write("Quartic step is not acceptable. ")
elif success:
fh.write("Quartic step is not acceptable. ")
# If no condition is met, then y0 is unacceptable. Try the cubic fit next
fh.write("Attempting a cubic fit.\n")
success, y0 = cubic_fit(X0, X1, e0, e1, G0, G1, fh)
if success and y0 > 0:
if y0 < 1:
new_point = X0 + y0*(X1 - X0)
new_gradient = interpolate_gradient(G0, G1, y0)
new_gradient = perpendicular_projection(new_gradient, X1 - X0)
fh.write("Line search succeded.\n")
return True, new_point, new_gradient
else:
previous_step = X1 - internal['position'][-2]
previous_step_size = linalg.norm(previous_step)
new_point = X0 + y0*(X1 - X0)
step = new_point - X1
step_size = linalg.norm(step)
if step_size < previous_step_size:
new_gradient = interpolate_gradient(G0, G1, y0)
new_gradient = perpendicular_projection(new_gradient, X1 - X0)
fh.write("Line search succeded.\n")
return True, new_point, new_gradient
else:
fh.write("Cubic step is not acceptable. ")
elif success:
fh.write("Cubic step is not acceptable. ")
# If no condition is met again, then all fits fail.
fh.write("All fits fail. ")
# Then, if the latest point is not the best, use y0 = 0.5 (step to the midpoint)
idx = best_point()
if idx < len(internal['energy']) - 1:
y0 = 0.5
new_point = X0 + y0*(X1 - X0)
new_gradient = interpolate_gradient(G0, G1, y0)
new_gradient = perpendicular_projection(new_gradient, X1 - X0)
fh.write("Moving to the midpoint.\n")
return True, new_point, new_gradient
# If the latest point is the best point, no linear search is done
fh.write("No linear search will be used in this step.\n")
return False, None, None
## For cubic and quartic fits, G0 and G1 are the gradient vectors
def cubic_fit(X0, X1, e0, e1, G0, G1, fh):
line = X1 - X0
line /= linalg.norm(line)
g0 = np.dot(G0, line)
g1 = np.dot(G1, line)
De = e1 - e0
fh.write("De = {:<18.15e} g0 = {:<12.8f} g1 = {:<12.8f}\n".format(De, g0, g1))
alpha = g1 + g0 - 2*De
if abs(alpha) < epsilon:
fh.write("Cubic fit failed: alpha too small\n")
return False, None
beta = 3*De - 2*g0 - g1
discriminant = 4 * (beta**2 - 3*alpha*g0)
if discriminant < 0:
fh.write("Cubic fit failed: no minimum found (negative Delta)\n")
return False, None
if abs(discriminant) < epsilon:
fh.write("Cubic fit failed: no minimum found (null Delta)\n")
return False, None
y0 = (-beta + math.sqrt(discriminant/4)) / (3*alpha)
fh.write("Minimum found with y0 = {:<8.4f}\n".format(y0))
return True, y0
def quartic_fit(X0, X1, e0, e1, G0, G1, fh):
line = X1 - X0
line /= linalg.norm(line)
g0 = np.dot(G0, line)
g1 = np.dot(G1, line)
De = e1 - e0
Dg = g1 - g0
fh.write("De = {:<18.15e} g0 = {:<12.8f} g1 = {:<12.8f}\n".format(De, g0, g1))
if Dg < 0 or De - g0 < 0:
fh.write("Quartic fit failed: negative alpha\n")
return False, None
if abs(Dg) < epsilon or abs(De - g0) < epsilon:
fh.write("Quartic fit failed: alpha too small\n")
return False, None
discriminant = 16 * (Dg**2 - 3*(g1 + g0 - 2*De)**2)
if discriminant < 0:
fh.write("Quartic fit failed: no minimum found (negative Delta)\n")
return False, None
alpha1 = (Dg + math.sqrt(discriminant/16)) / 2
alpha2 = (Dg - math.sqrt(discriminant/16)) / 2
fh.write("alpha1 = {:<7.4e} alpha2 = {:<7.4e}\n".format(alpha1, alpha2))
alpha = alpha1
beta = g1 + g0 - 2*De - 2*alpha
gamma = De - g0 - alpha - beta
y0 = (-1/(2*alpha)) * ((beta**3 - 4*alpha*beta*gamma + 8*g0*alpha**2)/4)**(1/3)
fh.write("Minimum found with y0 = {:<8.4f}\n".format(y0))
return True, y0
def rfo_step(gradient, hessian, type):
dim = len(gradient)
aug_hessian = []
for i in range(dim):
aug_hessian.extend(hessian[i,:].tolist())
aug_hessian.append(gradient[i])
aug_hessian.extend(gradient.tolist())
aug_hessian.append(0)
aug_hessian = np.array(aug_hessian).reshape(dim + 1, dim + 1)
evals, evecs = linalg.eigh(aug_hessian)
if type == "min":
step = np.array(evecs[:-1,0])
elif type == "ts":
step = np.array(evecs[:-1,1])
return step
def update_trust_radius():
if internal['trust_radius'] == None:
internal['trust_radius'] = player['maxstep']
elif len(internal['energy']) > 1:
X1 = internal['position'][-1]
X0 = internal['position'][-2]
Dx = X1 - X0
displace = linalg.norm(Dx)
e1 = internal['energy'][-1]
e0 = internal['energy'][-2]
De = e1 - e0
g0 = internal['gradient'][-2]
h0 = internal['hessian'][-2]
rho = De / (np.dot(g0, Dx) + 0.5*np.dot(Dx, np.matmul(h0, Dx.T).T))
if rho > 0.75 and displace > 0.8*internal['trust_radius']:
internal['trust_radius'] = 2*internal['trust_radius']
elif rho < 0.25:
internal['trust_radius'] = 0.25*displace
return
def interpolate_gradient(G0, G1, y0):
DG = G1 - G0
gradient = G0 + y0*DG
return gradient
def perpendicular_projection(vector, line):
direction = line / linalg.norm(line)
projection = np.dot(vector, direction) * direction
return vector - projection

35
diceplayer/DPpack/PTable.py Normal file
View File

@@ -0,0 +1,35 @@
#### Label used in Dice for a ghost atom
dice_ghost_label = "Xx"
#### Tuple of atom symbols
atomsymb = ( "00",
"H ", "He",
"Li","Be", "B ","C ","N ","O ","F ","Ne",
"Na","Mg", "Al","Si","P ","S ","Cl","Ar",
"K ","Ca","Sc","Ti","V ","Cr","Mn","Fe","Co","Ni","Cu","Zn","Ga","Ge","As","Se","Br","Kr",
"Rb","Sr","Y ","Zr","Nb","Mo","Tc","Ru","Rh","Pd","Ag","Cd","In","Sn","Sb","Te","I ","Xe",
"Cs","Ba",
"La","Ce","Pr","Nd","Pm","Sm","Eu","Gd","Tb","Dy","Ho","Er","Tm","Yb","Lu",
"Hf","Ta","W ","Re","Os","Ir","Pt","Au","Hg","Ti","Pb","Bi","Po","At","Rn",
"Fr","Ra",
"Ac","Th","Pa","U ","Np","Pu","Am","Cm","Bk","Cf","Es","Fm","Md","No","Lr",
dice_ghost_label )
#### Tuple of atom masses
atommass = ( 0.0,
1.0079, 4.0026,
6.9410,9.0122, 10.811,12.011,14.007,15.999,18.998,20.180,
22.990,24.305, 26.982,28.086,30.974,32.065,35.453,39.948,
39.098,40.078,44.956,47.867,50.942,51.996,54.938,55.845,58.933,58.693,63.546,65.409,69.723,72.640,74.922,78.960,79.904,83.798,
85.468,87.620,88.906,91.224,92.906,95.940,98.000,101.07,102.91,106.42,107.87,112.41,114.82,118.71,121.76,127.60,126.90,131.29,
132.91,137.33,
138.91,140.12,140.91,144.24,145.00,150.36,151.96,157.25,158.93,162.50,164.93,167.26,168.93,173.04,174.97,
178.49,180.95,183.84,186.21,190.23,192.22,195.08,196.97,200.59,204.38,207.20,208.98,209.00,210.00,222.00,
223.00,226.00,
227.00,232.04,231.04,238.03,237.00,244.00,243.00,247.00,247.00,251.00,252.00,257.00,258.00,259.00,262.00,
0.000 )
#### Number of the ghost atom
ghost_number = len(atomsymb) - 1

1729
diceplayer/DPpack/SetGlobals.py Normal file
View File

File diff suppressed because it is too large Load Diff

0
diceplayer/DPpack/__init__.py Normal file
View File

342
diceplayer/__main__.py Normal file
View File

@@ -0,0 +1,342 @@
#!/usr/bin/python3
import os, sys, time, signal
import setproctitle
import argparse
import shutil
from multiprocessing import Process, connection
from diceplayer.DPpack.PTable import *
from diceplayer.DPpack.SetGlobals import *
from diceplayer.DPpack.MolHandling import *
from diceplayer.DPpack.Misc import *
_version = 'dev'
setproctitle.setproctitle("diceplayer-{}".format(_version))
if __name__ == '__main__':
#### Read and store the arguments passed to the program ####
#### and set the usage and help messages ####
parser = argparse.ArgumentParser(prog='Diceplayer')
parser.add_argument('--continue', dest='opt_continue' , default=False, action='store_true')
parser.add_argument('--version', action='version', version='%(prog)s 1.0')
parser.add_argument('-i', dest='infile', default='control.in', metavar='INFILE',
help='input file of diceplayer [default = control.in]')
parser.add_argument('-o', dest='outfile', default='run.log', metavar='OUTFILE',
help='output file of diceplayer [default = run.log]')
## Study the option of a parameter for continuing the last process via data from control.in and run.log files
args = parser.parse_args()
#### Open OUTFILE for writing and print keywords and initial info
try:
if args.opt_continue and os.path.exists(args.outfile):
outfile = open(args.outfile,'r')
run_file = outfile.readlines()
control_sequence = ' Step # '
sucessfull_sequence = '+----------------------------------------------------------------------------------------+'
for line in run_file:
if control_sequence in line:
cyc = int(line[-2])
if sucessfull_sequence in line:
cyc += 1
outfile.close()
if os.path.isfile(args.outfile+".backup"):
os.remove(args.outfile+".backup")
os.rename(args.outfile,args.outfile+".backup")
outfile = open(args.outfile,'w',1)
elif os.path.exists(args.outfile):
os.rename(args.outfile, args.outfile+".backup")
outfile = open(args.outfile,'w',1)
else:
outfile = open(args.outfile,"w",1)
except EnvironmentError as err:
sys.exit(err)
try:
if os.path.exists(args.infile):
infile = open(args.infile,"r")
except EnvironmentError as err:
sys.exit(err)
#### Read and check the keywords in INFILE
internal = Internal(infile, outfile)
internal.read_keywords()
if args.opt_continue:
try:
internal.player.initcyc = cyc
except:
sys.exit("Error: There is no sutable run.log file to continue the previous process")
internal.check_keywords()
internal.print_keywords()
# #### Check whether the executables are in the path
internal.check_executables()
# #### Read the potential, store the info in 'molecules' and prints the info in OUTFILE
internal.read_potential()
# if internal.player.lps == "yes":
# read_lps()
# if internal.player.ghosts == "yes":
# read_ghosts()
internal.print_potential()
#### Bring the molecules to standard orientation and prints info about them
for i in range(len(internal.system.molecule)):
internal.outfile.write("\nMolecule type {} - {}:\n\n".format(i + 1, internal.system.molecule[i].molname))
internal.system.molecule[i].print_mol_info(internal.outfile)
internal.outfile.write(" Translating and rotating molecule to standard orientation...")
internal.system.molecule[i].standard_orientation()
internal.outfile.write(" Done\n\n New values:\n")
internal.system.molecule[i].print_mol_info(internal.outfile)
internal.outfile.write(90 * "=")
internal.outfile.write("\n")
if not args.opt_continue:
make_simulation_dir()
else:
simdir = "simfiles"
stepdir = "step{:02d}".format(internal.player.initcyc)
if os.path.exists(simdir+os.sep+stepdir):
shutil.rmtree(simdir+os.sep+stepdir)
#### Open the geoms.xyz file and prints the initial geometry if starting from zero
if internal.player.initcyc == 1:
try:
path = "geoms.xyz"
geomsfh = open(path, "w", 1)
except EnvironmentError as err:
sys.exit(err)
internal.system.print_geom(0, geomsfh)
else:
try:
path = "geoms.xyz"
geomsfh = open(path, "a", 1)
except EnvironmentError as err:
sys.exit(err)
internal.outfile.write("\nStarting the iterative process.\n")
## Initial position (in Bohr)
position = internal.system.molecule[0].read_position()
## If restarting, read the last gradient and hessian
# if internal.player.initcyc > 1:
# if internal.player.qmprog in ("g03", "g09", "g16"):
# Gaussian.read_forces("grad_hessian.dat")
# Gaussian.read_hessian_fchk("grad_hessian.dat")
#if player['qmprog'] == "molcas":
#Molcas.read_forces("grad_hessian.dat")
#Molcas.read_hessian("grad_hessian.dat")
###
### Start the iterative process
###
internal.outfile.write("\n" + 90 * "-" + "\n")
for cycle in range(internal.player.initcyc, internal.player.initcyc + internal.player.maxcyc):
internal.outfile.write("{} Step # {}\n".format(40 * " ", cycle))
internal.outfile.write(90 * "-" + "\n\n")
make_step_dir(cycle)
####
#### Start block of parallel simulations
####
procs = []
sentinels = []
for proc in range(1, internal.player.nprocs + 1):
p = Process(target=internal.simulation_process, args=(cycle, proc))
p.start()
procs.append(p)
sentinels.append(p.sentinel)
while procs:
finished = connection.wait(sentinels)
for proc_sentinel in finished:
i = sentinels.index(proc_sentinel)
status = procs[i].exitcode
procs.pop(i)
sentinels.pop(i)
if status != 0:
for p in procs:
p.terminate()
sys.exit(status)
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
###
## Make ASEC
# internal.outfile.write("\nBuilding the ASEC and vdW meanfields... ")
# asec_charges = internal.populate_asec_vdw(cycle)
# ## After ASEC is built, compress files bigger than 1MB
# for proc in range(1, internal.player.nprocs + 1):
# path = "step{:02d}".format(cycle) + os.sep + "p{:02d}".format(proc)
# compress_files_1mb(path)
###
### Start QM calculation
###
make_qm_dir(cycle)
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)
# ## 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
# ## 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)
# ## 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")
# ## Save gradient and hessian
# internal.gaussian.print_grad_hessian(cycle, cur_gradient, hessian)
# ## Calculate the step and update the position
# step = internal.calculate_step(cur_gradient, hessian, internal.outfile)
# position += step
# ## Update the geometry of the reference molecule
# internal.system.molecule[0].update_molecule(position, internal.outfile)
# ## If needed, calculate the charges
# if cycle < internal.player.switchcyc:
# # internal.gaussian.make_charge_input(cycle, asec_charges)
# internal.gaussian.run_gaussian(cycle, "charge", internal.outfile)
# ## Read the new charges and update molecules[0]
# if cycle < internal.player.switchcyc:
# file = "simfiles" + os.sep + "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec2.log"
# internal.gaussian.read_charges(file, internal.outfile)
# else:
# file = "simfiles" + os.sep + "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec.log"
# internal.gaussian.read_charges(file, internal.outfile)
# ## Print new info for molecule[0]
# internal.outfile.write("\nNew values for molecule type 1:\n\n")
# internal.system.molecule[0].print_mol_info()
# ## Print new geometry in geoms.xyz
# internal.system.molecule[0].print_geom(cycle, geomsfh)
# ##
# ## Molcas block
# ##
# #if player['qmprog'] == "molcas":
# #elif player['opt'] == "ts":
# ##
# ## Gaussian block
# ##
# #if player['qmprog'] in ("g03", "g09", "g16"):
# ##
# ## Molcas block
# ##
# #if player['qmprog'] == "molcas":
# else: ## Only relax the charge distribution
# 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.gaussian.make_charge_input(cycle, asec_charges)
# internal.gaussian.run_gaussian(cycle, "charge", internal.outfile)
# file = "simfiles" + os.sep + "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec2.log"
# internal.read_charges(file)
# ## Print new info for molecule[0]
# internal.outfile.write("\nNew values for molecule type 1:\n\n")
# internal.system.molecule[0].print_mol_info()
# #if player['qmprog'] == "molcas":
# ####
# #### 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, ...)
# # internal.outfile.write("\nDiceplayer finished normally!\n")
# # internal.outfile.close()
# # ####
# # #### End of the program
# # ####