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

v0.0.1

Browse Source 
Signed-off-by: Vitor Hideyoshi <vitor.h.n.batista@gmail.com>
This commit is contained in:
Vitor Hideyoshi
2021-06-17 20:01:42 -03:00
committed by Vitor Hideyoshi
commit 2899031b71
31 changed files with 12048 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

BIN
.control.in.swp Normal file
View File

Binary file not shown.

534
DPpack/Dice.py Normal file
View File

@@ -0,0 +1,534 @@
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 *
dice_end_flag = "End of simulation" ## The normal end flag
dice_flag_line = -2 ## must be in the line before the last
umaAng3_to_gcm3 = 1.6605 ## Conversion between uma/Ang3 to g/cm3
max_seed = 4294967295 ## Maximum allowed value for a seed (numpy)
####################################### 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

420
DPpack/Gaussian.py Normal file
View File

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

63
DPpack/Misc.py Normal file
View File

@@ -0,0 +1,63 @@
import os, sys, time
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_step_dir(cycle):
step_dir = "step{:02d}".format(cycle)
if os.path.exists(step_dir):
sys.exit("Error: a file or directory {} already exists".format(step_dir))
try:
os.makedirs(step_dir)
except:
sys.exit("Error: cannot make directory {}".format(step_dir))
return
def make_qm_dir(cycle):
step_dir = "step{:02d}".format(cycle)
path = step_dir + os.sep + "qm"
try:
os.makedirs(path)
except:
sys.exit("Error: cannot make directory {}".format(path))
return

603
DPpack/MolHandling.py Normal file
View File

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

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

813
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'] = 1f
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

0
DPpack/__init__.py Normal file
View File

BIN
DPpack/__pycache__/Dice.cpython-36.pyc Normal file
View File

Binary file not shown.

BIN
DPpack/__pycache__/Misc.cpython-36.pyc Normal file
View File

Binary file not shown.

BIN
DPpack/__pycache__/MolHandling.cpython-36.pyc Normal file
View File

Binary file not shown.

BIN
DPpack/__pycache__/PTable.cpython-36.pyc Normal file
View File

Binary file not shown.

BIN
DPpack/__pycache__/SetGlobals.cpython-36.pyc Normal file
View File

Binary file not shown.

BIN
DPpack/__pycache__/__init__.cpython-36.pyc Normal file
View File

Binary file not shown.

BIN
DPpack/__pycache__/mol_handling.cpython-36.pyc Normal file
View File

Binary file not shown.

21
control.in Normal file
View File

@@ -0,0 +1,21 @@
# diceplayer
initcyc = 1
maxcyc = 3
opt = NO
nprocs = 2
qmprog = g09
lps = no
ghosts = no
altsteps = 20000
# dice
ncores = 1
nmol = 1 100
dens = 0.75
nstep = 40000 60000 50000
isave = 1000
ljname = phb.pot
outname = phb
# Gaussian
level = MP2/aug-cc-pVTZ

295
diceplayer-0.0.1.py Normal file
View File

@@ -0,0 +1,295 @@
#!/export/apps/python/361/bin/python3
import os, sys, time, signal
import argparse
import shutil
from multiprocessing import Process, connection
import DPpack.Dice as Dice
import DPpack.Gaussian as Gaussian
from DPpack.PTable import *
from DPpack.SetGlobals import *
from DPpack.MolHandling import *
from DPpack.Misc import *
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('--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]')
args = parser.parse_args()
#### Read and check the keywords in INFILE
read_keywords(args.infile)
check_keywords(args.infile)
#### Open OUTFILE for writing and print keywords and initial info
try:
if player['initcyc'] > 1 and os.path.exists(args.outfile):
oldname = args.outfile + ".old"
os.replace(args.outfile, oldname)
logfh = open(args.outfile, 'w', 1)
except EnvironmentError as err:
sys.exit(err)
print_keywords(logfh)
#### Check whether the executables are in the path
check_executables(logfh)
#### Read the potential, store the info in 'molecules' and prints the info in OUTFILE
read_potential(args.infile)
if player['lps'] == "yes":
read_lps()
if player['ghosts'] == "yes":
read_ghosts()
print_potential(logfh)
#### Bring the molecules to standard orientation and prints info about them
for i in range(len(molecules)):
logfh.write("\nMolecule type {}:\n\n".format(i + 1))
print_mol_info(molecules[i], logfh)
logfh.write(" Translating and rotating molecule to standard orientation...")
standard_orientation(molecules[i])
logfh.write(" Done\n\n New values:\n")
print_mol_info(molecules[i], logfh)
logfh.write(90 * "=")
logfh.write("\n")
#### Open the geoms.xyz file and prints the initial geometry if starting from zero
if player['initcyc'] == 1:
try:
geomsfh = open("geoms.xyz", "w", 1)
except EnvironmentError as err:
sys.exit(err)
print_geom(0, geomsfh)
else:
try:
geomsfh = open("geoms.xyz", "A", 1)
except EnvironmentError as err:
sys.exit(err)
logfh.write("\nStarting the iterative process.\n")
## Initial position (in Bohr)
position = read_position(molecules[0])
## If restarting, read the last gradient and hessian
if player['initcyc'] > 1:
if 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
####
for cycle in range(player['initcyc'], player['initcyc'] + player['maxcyc']):
logfh.write("\n" + 90 * "-" + "\n")
logfh.write("{} Step # {}\n".format(40 * " ", cycle))
logfh.write(90 * "-" + "\n\n")
make_step_dir(cycle)
if player['altsteps'] == 0 or cycle == 1:
dice['randominit'] = True
else:
dice['randominit'] = False
####
#### Start block of parallel simulations
####
procs = []
sentinels = []
for proc in range(1, player['nprocs'] + 1):
p = Process(target=Dice.simulation_process, args=(cycle, proc, logfh))
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, player['nprocs'] + 1):
Dice.print_last_config(cycle, proc)
####
#### End of parallel simulations block
####
## Make ASEC
logfh.write("\nBuilding the ASEC and vdW meanfields... ")
asec_charges = populate_asec_vdw(cycle, logfh)
## After ASEC is built, compress files bigger than 1MB
for proc in range(1, 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 player['opt'] == "yes":
##
## Gaussian block
##
if player['qmprog'] in ("g03", "g09", "g16"):
if cycle > 1:
src = "step{:02d}".format(cycle - 1) + os.sep + "qm" + os.sep + "asec.chk"
dst = "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec.chk"
shutil.copyfile(src, dst)
Gaussian.make_force_input(cycle, asec_charges)
Gaussian.run_gaussian(cycle, "force", logfh)
Gaussian.run_formchk(cycle, logfh)
## Read the gradient
file = "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec.fchk"
gradient = Gaussian.read_forces(file, logfh)
if len(cur_gradient) > 0:
old_gradient = cur_gradient
cur_gradient = gradient
## If 1st step, read the hessian
if cycle == 1:
if player['readhessian'] == "yes":
file = "grad_hessian.dat"
logfh.write("\nReading the hessian matrix from file {}\n".format(file))
hessian = Gaussian.read_hessian_fchk(file)
else:
file = "step01" + os.sep + "qm" + os.sep + "asec.fchk"
logfh.write("\nReading the hessian matrix from file {}\n".format(file))
hessian = Gaussian.read_hessian(file)
## From 2nd step on, update the hessian
else:
logfh.write("\nUpdating the hessian matrix using the BFGS method... ")
hessian = update_hessian(step, cur_gradient, old_gradient, hessian)
logfh.write("Done\n")
## Save gradient and hessian
Gaussian.print_grad_hessian(cycle, cur_gradient, hessian)
## Calculate the step and update the position
step = calculate_step(cur_gradient, hessian, logfh)
position += step
## Update the geometry of the reference molecule
update_molecule(position, logfh)
## If needed, calculate the charges
if cycle < player['switchcyc']:
Gaussian.make_charge_input(cycle, asec_charges)
Gaussian.run_gaussian(cycle, "charge", logfh)
## Read the new charges and update molecules[0]
if cycle < player['switchcyc']:
file = "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec2.log"
Gaussian.read_charges(file, logfh)
else:
file = "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec.log"
Gaussian.read_charges(file, logfh)
## Print new info for molecule[0]
logfh.write("\nNew values for molecule type 1:\n\n")
print_mol_info(molecules[0], logfh)
## Print new geometry in geoms.xyz
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 player['qmprog'] in ("g03", "g09", "g16"):
if cycle > 1:
src = "step{:02d}".format(cycle - 1) + os.sep + "qm" + os.sep + "asec.chk"
dst = "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec.chk"
shutil.copyfile(src, dst)
Gaussian.make_charge_input(cycle, asec_charges)
Gaussian.run_gaussian(cycle, "charge", logfh)
file = "step{:02d}".format(cycle) + os.sep + "qm" + os.sep + "asec2.log"
Gaussian.read_charges(file)
## Print new info for molecule[0]
logfh.write("\nNew values for molecule type 1:\n\n")
print_mol_info(molecules[0], logfh)
#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, ...)
logfh.write("\nDiceplayer finished normally!\n")
logfh.close()
####
#### End of the program
####

33
geoms.xyz Normal file
View File

@@ -0,0 +1,33 @@
31
Cycle # 0
C -4.394497 0.698769 0.421354
C -4.504613 -0.640849 -0.141092
C -3.266831 -1.267916 -0.601745
C -2.064301 -0.654971 -0.479838
C -1.966148 0.666069 0.106967
C -3.195469 1.321761 0.499266
O -5.604026 -1.213715 -0.270977
N -0.863607 1.365979 0.273247
C 0.394047 0.812820 0.140269
C 1.418955 1.646567 -0.341744
C 0.757396 -0.471466 0.587991
C 2.717900 1.193894 -0.471504
C 2.068121 -0.911914 0.507136
C 3.079906 -0.111825 -0.067434
N 4.367937 -0.572075 -0.209272
C 5.421862 0.397502 -0.451662
C 4.752115 -1.767634 0.518351
H 6.365079 -0.133224 -0.533034
H 5.255789 0.923150 -1.389081
H 5.504293 1.135755 0.351310
H 4.140103 -2.613567 0.214988
H 5.783262 -2.003292 0.275116
H 4.665889 -1.644881 1.602460
H 0.022604 -1.091753 1.084030
H 2.307262 -1.879075 0.920974
H 3.462774 1.870724 -0.861110
H 1.165469 2.658467 -0.627261
H -3.100673 2.322526 0.898873
H -5.305838 1.177165 0.752078
H -3.356658 -2.230234 -1.086826
H -1.173708 -1.112751 -0.887177

6
ghosts.in Normal file
View File

@@ -0,0 +1,6 @@
G 2 4 6 8 10 12 13 15 17 20
M 2 3 4 5 6 7 8 19 20 22
Z 23 25 26 27 28 29

4
lps.in Normal file
View File

@@ -0,0 +1,4 @@
1 12 13 14 109.5 0.95
2 9 10 11 50 1.2
#3 1 5 6
2 17 19 23

1602
new_cyclo.xyz Normal file
View File

File diff suppressed because it is too large Load Diff

1602
new_cyclo.xyz.new Normal file
View File

File diff suppressed because it is too large Load Diff

1782
new_cyclo.xyz.new2 Normal file
View File

File diff suppressed because it is too large Load Diff

11
othercontrol.in Normal file
View File

@@ -0,0 +1,11 @@
initcyc = -4 pode?
maxcyc = 1
forces = yes
nprocs = 2
# dice
#nmol = aslkj 1 20 slkdfj
dens = 0.75
nsteps = 7000 2500 3876
isave = 1340
ljname = phb.pot

51
phb.pot Normal file
View File

@@ -0,0 +1,51 @@
*
2
31 PHENOL BLUE geom/chg MP2/aug-cc-pVTZ Charge MK LJ OPLSAA
1 6 -4.400337 0.665513 0.414312 -0.38759300 0.0700 3.5500
1 6 -4.499567 -0.672788 -0.153273 0.80400100 0.0700 3.5500
1 6 -3.256300 -1.289216 -0.613514 -0.33849100 0.0700 3.5500
1 6 -2.058491 -0.668069 -0.486739 -0.22078300 0.0700 3.5500
1 6 -1.971143 0.651479 0.105111 0.56352700 0.0700 3.5500
1 6 -3.206008 1.296843 0.497122 -0.14537200 0.0700 3.5500
2 8 -5.594533 -1.253085 -0.287655 -0.71446100 0.2100 2.9600
3 7 -0.874044 1.358704 0.276359 -0.72179200 0.1700 3.2500
4 6 0.387854 0.815120 0.144104 0.60712300 0.0700 3.5500
4 6 1.407755 1.658002 -0.332615 -0.29547600 0.0700 3.5500
4 6 0.759497 -0.468150 0.587914 -0.34666100 0.0700 3.5500
4 6 2.710211 1.215186 -0.461194 -0.29854600 0.0700 3.5500
4 6 2.073539 -0.898835 0.508310 -0.23428400 0.0700 3.5500
4 6 3.080760 -0.089365 -0.061116 0.29465100 0.0700 3.5500
5 7 4.372382 -0.539792 -0.201825 -0.19049700 0.1700 3.2500
6 6 5.419802 0.438244 -0.438362 -0.21318300 0.0660 3.5000
6 6 4.763609 -1.735216 0.522254 -0.22542000 0.0660 3.5000
7 1 6.366998 -0.085365 -0.519616 0.10984000 0.0300 2.5000
7 1 5.251961 0.966123 -1.374211 0.07714000 0.0300 2.5000
7 1 5.495169 1.174115 0.367485 0.12575800 0.0300 2.5000
7 1 4.158375 -2.584419 0.214461 0.08322900 0.0300 2.5000
7 1 5.796953 -1.962538 0.280411 0.11536100 0.0300 2.5000
7 1 4.674161 -1.617074 1.606614 0.12144000 0.0300 2.5000
7 1 0.028134 -1.095538 1.080074 0.17943400 0.0300 2.4200
7 1 2.318762 -1.865762 0.919126 0.17725000 0.0300 2.4200
7 1 3.451018 1.898797 -0.846692 0.19884400 0.0300 2.4200
7 1 1.147588 2.669090 -0.614979 0.17338800 0.0300 2.4200
7 1 -3.119300 2.296790 0.900595 0.15153400 0.0300 2.4200
7 1 -5.315819 1.136106 0.744794 0.20528900 0.0300 2.4200
7 1 -3.338131 -2.250367 -1.102310 0.18436700 0.0300 2.4200
7 1 -1.163740 -1.117915 -0.893805 0.16038200 0.0300 2.4200
16
1 6 0.672026 -2.823446 0.002631 -0.11500000 0.0700 3.5500
1 6 2.072026 -2.823446 0.002631 -0.11500000 0.0700 3.5500
1 6 2.768235 -1.617645 0.002633 -0.11500000 0.0700 3.5500
1 6 2.068235 -0.405210 0.002635 -0.11500000 0.0700 3.5500
1 14 0.675894 -0.405219 0.002635 0.15000000 0.0700 3.5500
1 104 -0.024198 -1.617602 0.002633 -0.11500000 0.0700 3.5500
2 1 0.132026 -3.758753 0.002629 0.11500000 0.0300 2.4200
2 1 2.612026 -3.758753 0.002629 0.11500000 0.0300 2.4200
2 1 2.608235 0.530098 0.002636 0.11500000 0.0300 2.4200
2 1 -1.104198 -1.617602 0.002633 0.11500000 0.0300 2.4200
3 8 -0.004114 0.772571 0.002637 -0.58500000 0.1700 3.0700
3 1 0.619778 1.502200 0.002638 0.43500000 0.0000 0.0000
4 6 4.278235 -1.617645 0.002633 0.11500000 0.1700 3.8000
4 1 4.634902 -0.743994 0.507036 0.00000000 0.0000 0.0000
4 1 4.634902 -2.491297 0.507036 0.00000000 0.0000 0.0000
4 1 4.634902 -1.617645 -1.006173 0.00000000 0.0000 0.0000

1633
phb.xyz.last-p01 Normal file
View File

File diff suppressed because it is too large Load Diff

1633
phb.xyz.last-p02 Normal file
View File

File diff suppressed because it is too large Load Diff

229
run.log Normal file
View File

@@ -0,0 +1,229 @@
##########################################################################################
############# Welcome to DICEPLAYER version 1.0 #############
##########################################################################################
Your python version is 3.6.1 (default, Apr 24 2017, 06:18:27)
[GCC 4.2.1 Compatible Apple LLVM 8.0.0 (clang-800.0.42.1)]
Program started on Wednesday, 14 Jun 2017 at 12:30:02
Environment variables:
OMP_STACKSIZE = 32M
==========================================================================================
CONTROL variables being used in this run:
------------------------------------------------------------------------------------------
altsteps = 20000
ghosts = no
initcyc = 1
lps = no
maxcyc = 3
maxstep = 0.3
nprocs = 2
opt = no
qmprog = g09
switch_step = 3
zipprog = gzip
------------------------------------------------------------------------------------------
DICE variables being used in this run:
------------------------------------------------------------------------------------------
dens = 0.75
isave = 1000
ljname = phb.pot
ncores = 1
nmol = 1 100
nstep = 40000 60000 50000
outname = phb
press = 1.0
temp = 300.0
title = Diceplayer run
------------------------------------------------------------------------------------------
GAUSSIAN variables being used in this run:
------------------------------------------------------------------------------------------
chgmult = 0 1
level = MP2/aug-cc-pVTZ
pop = chelpg
==========================================================================================
Program dice found at /usr/local/bin/dice
==========================================================================================
Potential parameters from file phb.pot:
------------------------------------------------------------------------------------------
Combination rule: *
Types of molecules: 2
31 atoms in molecule type 1:
---------------------------------------------------------------------------------
Lbl AN X Y Z Charge Epsilon Sigma Mass
---------------------------------------------------------------------------------
1 6 -4.40034 0.66551 0.41431 -0.387593 0.07000 3.5500 12.0110
1 6 -4.49957 -0.67279 -0.15327 0.804001 0.07000 3.5500 12.0110
1 6 -3.25630 -1.28922 -0.61351 -0.338491 0.07000 3.5500 12.0110
1 6 -2.05849 -0.66807 -0.48674 -0.220783 0.07000 3.5500 12.0110
1 6 -1.97114 0.65148 0.10511 0.563527 0.07000 3.5500 12.0110
1 6 -3.20601 1.29684 0.49712 -0.145372 0.07000 3.5500 12.0110
2 8 -5.59453 -1.25309 -0.28765 -0.714461 0.21000 2.9600 15.9990
3 7 -0.87404 1.35870 0.27636 -0.721792 0.17000 3.2500 14.0070
4 6 0.38785 0.81512 0.14410 0.607123 0.07000 3.5500 12.0110
4 6 1.40776 1.65800 -0.33261 -0.295476 0.07000 3.5500 12.0110
4 6 0.75950 -0.46815 0.58791 -0.346661 0.07000 3.5500 12.0110
4 6 2.71021 1.21519 -0.46119 -0.298546 0.07000 3.5500 12.0110
4 6 2.07354 -0.89884 0.50831 -0.234284 0.07000 3.5500 12.0110
4 6 3.08076 -0.08936 -0.06112 0.294651 0.07000 3.5500 12.0110
5 7 4.37238 -0.53979 -0.20183 -0.190497 0.17000 3.2500 14.0070
6 6 5.41980 0.43824 -0.43836 -0.213183 0.06600 3.5000 12.0110
6 6 4.76361 -1.73522 0.52225 -0.225420 0.06600 3.5000 12.0110
7 1 6.36700 -0.08536 -0.51962 0.109840 0.03000 2.5000 1.0079
7 1 5.25196 0.96612 -1.37421 0.077140 0.03000 2.5000 1.0079
7 1 5.49517 1.17412 0.36749 0.125758 0.03000 2.5000 1.0079
7 1 4.15838 -2.58442 0.21446 0.083229 0.03000 2.5000 1.0079
7 1 5.79695 -1.96254 0.28041 0.115361 0.03000 2.5000 1.0079
7 1 4.67416 -1.61707 1.60661 0.121440 0.03000 2.5000 1.0079
7 1 0.02813 -1.09554 1.08007 0.179434 0.03000 2.4200 1.0079
7 1 2.31876 -1.86576 0.91913 0.177250 0.03000 2.4200 1.0079
7 1 3.45102 1.89880 -0.84669 0.198844 0.03000 2.4200 1.0079
7 1 1.14759 2.66909 -0.61498 0.173388 0.03000 2.4200 1.0079
7 1 -3.11930 2.29679 0.90060 0.151534 0.03000 2.4200 1.0079
7 1 -5.31582 1.13611 0.74479 0.205289 0.03000 2.4200 1.0079
7 1 -3.33813 -2.25037 -1.10231 0.184367 0.03000 2.4200 1.0079
7 1 -1.16374 -1.11791 -0.89380 0.160382 0.03000 2.4200 1.0079
16 atoms in molecule type 2:
---------------------------------------------------------------------------------
Lbl AN X Y Z Charge Epsilon Sigma Mass
---------------------------------------------------------------------------------
1 6 0.67203 -2.82345 0.00263 -0.115000 0.07000 3.5500 12.0110
1 6 2.07203 -2.82345 0.00263 -0.115000 0.07000 3.5500 12.0110
1 6 2.76823 -1.61764 0.00263 -0.115000 0.07000 3.5500 12.0110
1 6 2.06824 -0.40521 0.00264 -0.115000 0.07000 3.5500 12.0110
1 14 0.67589 -0.40522 0.00264 0.150000 0.07000 3.5500 28.0860
1 104 -0.02420 -1.61760 0.00263 -0.115000 0.07000 3.5500 0.0000
2 1 0.13203 -3.75875 0.00263 0.115000 0.03000 2.4200 1.0079
2 1 2.61203 -3.75875 0.00263 0.115000 0.03000 2.4200 1.0079
2 1 2.60824 0.53010 0.00264 0.115000 0.03000 2.4200 1.0079
2 1 -1.10420 -1.61760 0.00263 0.115000 0.03000 2.4200 1.0079
3 8 -0.00411 0.77257 0.00264 -0.585000 0.17000 3.0700 15.9990
3 1 0.61978 1.50220 0.00264 0.435000 0.00000 0.0000 1.0079
4 6 4.27823 -1.61764 0.00263 0.115000 0.17000 3.8000 12.0110
4 1 4.63490 -0.74399 0.50704 0.000000 0.00000 0.0000 1.0079
4 1 4.63490 -2.49130 0.50704 0.000000 0.00000 0.0000 1.0079
4 1 4.63490 -1.61764 -1.00617 0.000000 0.00000 0.0000 1.0079
==========================================================================================
Molecule type 1:
Center of mass = ( -0.0000 , 0.0000 , 0.0000 )
Moments of inertia = 3.144054E+02 2.801666E+03 3.027366E+03
Major principal axis = ( 0.999972 , 0.007210 , 0.002184 )
Inter principal axis = ( -0.007218 , 0.999967 , 0.003660 )
Minor principal axis = ( -0.002157 , -0.003676 , 0.999991 )
Characteristic lengths = ( 11.96 , 5.25 , 2.98 )
Total mass = 226.28 au
Total charge = -0.0000 e
Dipole moment = ( 9.8367 , 0.6848 , 0.8358 ) Total = 9.8959 Debye
Translating and rotating molecule to standard orientation... Done
New values:
Center of mass = ( -0.0000 , 0.0000 , -0.0000 )
Moments of inertia = 3.144054E+02 2.801666E+03 3.027366E+03
Major principal axis = ( 1.000000 , 0.000000 , 0.000000 )
Inter principal axis = ( -0.000000 , 1.000000 , -0.000000 )
Minor principal axis = ( 0.000000 , 0.000000 , 1.000000 )
Characteristic lengths = ( 11.97 , 5.27 , 2.99 )
Total mass = 226.28 au
Total charge = -0.0000 e
Dipole moment = ( 9.8432 , 0.6168 , 0.8120 ) Total = 9.8959 Debye
Molecule type 2:
Center of mass = ( 1.6067 , -1.0929 , 0.0026 )
Moments of inertia = 1.205279E+02 2.859254E+02 4.033761E+02
Major principal axis = ( 0.795913 , -0.605411 , -0.000001 )
Inter principal axis = ( 0.605411 , 0.795913 , 0.000001 )
Minor principal axis = ( 0.000000 , -0.000002 , 1.000000 )
Characteristic lengths = ( 5.74 , 5.26 , 1.51 )
Total mass = 112.20 au
Total charge = -0.0000 e
Dipole moment = ( 2.3293 , 0.1593 , -0.0000 ) Total = 2.3347 Debye
Translating and rotating molecule to standard orientation... Done
New values:
Center of mass = ( 0.0000 , 0.0000 , 0.0000 )
Moments of inertia = 1.205279E+02 2.859254E+02 4.033761E+02
Major principal axis = ( -1.000000 , 0.000000 , 0.000000 )
Inter principal axis = ( 0.000000 , 1.000000 , 0.000000 )
Minor principal axis = ( 0.000000 , 0.000000 , 1.000000 )
Characteristic lengths = ( 5.67 , 5.13 , 1.51 )
Total mass = 112.20 au
Total charge = -0.0000 e
Dipole moment = ( 1.7575 , 1.5369 , -0.0000 ) Total = 2.3347 Debye
==========================================================================================
Starting the iterative process.
------------------------------------------------------------------------------------------
Step # 1
------------------------------------------------------------------------------------------
Simulation process p01 initiated with pid 6822
Simulation process p02 initiated with pid 6823
p01> NVT thermalization initiated (from random configuration) on 14 Jun 2017 at 12:30:02
p02> NVT thermalization initiated (from random configuration) on 14 Jun 2017 at 12:30:02
p02> NPT thermalization initiated on 14 Jun 2017 at 12:41:16
p01> NPT thermalization initiated on 14 Jun 2017 at 12:41:17
p02> NPT production initiated on 14 Jun 2017 at 13:01:25
p01> NPT production initiated on 14 Jun 2017 at 13:01:27
p02> ----- NPT production finished on 14 Jun 2017 at 13:57:41
p01> ----- NPT production finished on 14 Jun 2017 at 13:57:53
Building the ASEC and vdW meanfields... In average, 99.97 molecules were selected from the production simulations to form the
ASEC comprising a shell with minimum thickness of 15.352263400006278 AngstromDone.
------------------------------------------------------------------------------------------
Step # 2
------------------------------------------------------------------------------------------
Simulation process p01 initiated with pid 6903
Simulation process p02 initiated with pid 6904
p02> NPT thermalization initiated (from previous configuration) on 14 Jun 2017 at 13:58:52
p01> NPT thermalization initiated (from previous configuration) on 14 Jun 2017 at 13:58:52
p02> NPT production initiated on 14 Jun 2017 at 14:05:14
p01> NPT production initiated on 14 Jun 2017 at 14:05:20
p02> ----- NPT production finished on 14 Jun 2017 at 14:23:58
p01> ----- NPT production finished on 14 Jun 2017 at 14:24:05
Building the ASEC and vdW meanfields... In average, 99.66 molecules were selected from the production simulations to form the
ASEC comprising a shell with minimum thickness of 14.942583400006276 AngstromDone.
------------------------------------------------------------------------------------------
Step # 3
------------------------------------------------------------------------------------------
Simulation process p01 initiated with pid 6945
p01> NPT thermalization initiated (from previous configuration) on 14 Jun 2017 at 14:24:45
Simulation process p02 initiated with pid 6946
p02> NPT thermalization initiated (from previous configuration) on 14 Jun 2017 at 14:24:45
p02> NPT production initiated on 14 Jun 2017 at 14:31:05
p01> NPT production initiated on 14 Jun 2017 at 14:31:06
p02> ----- NPT production finished on 14 Jun 2017 at 14:48:15
p01> ----- NPT production finished on 14 Jun 2017 at 14:48:15
Building the ASEC and vdW meanfields... In average, 99.0 molecules were selected from the production simulations to form the
ASEC comprising a shell with minimum thickness of 14.79446440000628 AngstromDone.
Diceplayer finished normally!

64
test4.py Normal file
View File

@@ -0,0 +1,64 @@
#!/opt/local/bin/python3
import sys, math
from math import pi
import numpy as np
from numpy import linalg
from numpy.random import rand, random
from copy import deepcopy
from DPpack.MolHandling import *
from DPpack.PTable import *
natoms = input("# of atoms: ")
natoms = int(natoms)
file = input("xyz file: ")
molecules=[]
with open(file) as fh:
xyzfile = fh.readlines()
file += ".new"
file2 = file + "2"
fh = open(file, "w")
fh2 = open(file2, "w")
total_atoms = int(xyzfile.pop(0).split()[0])
fh.write("{}\n".format(total_atoms))
comment = xyzfile.pop(0)
fh.write("{}".format(comment))
nmols = round(total_atoms/natoms)
for i in range(nmols):
molecules.append([])
for j in range(natoms):
molecules[i].append({})
line = xyzfile.pop(0).split()
molecules[i][j]['na'] = int(line[0])
molecules[i][j]['rx'] = float(line[1])
molecules[i][j]['ry'] = float(line[2])
molecules[i][j]['rz'] = float(line[3])
molecules[i][j]['mass'] = atommass[molecules[0][j]['na']]
for atom in molecules[0]:
fh.write("{:>4s} {:>11.5f} {:>11.5f} {:>11.5f}\n".format(atomsymb[atom['na']],
atom['rx'], atom['ry'], atom['rz']))
fh2.write("{:>4s} {:>11.5f} {:>11.5f} {:>11.5f}\n".format(atomsymb[atom['na']],
atom['rx'], atom['ry'], atom['rz']))
for i in range(1, nmols):
fh2.write("{}\n".format(natoms))
fh2.write("{}".format(comment))
rmsd, projected_mol = rmsd_fit(molecules[i], molecules[0])
print("{:>9.5f}".format(rmsd))
for atom in projected_mol:
fh.write("{:>4s} {:>11.5f} {:>11.5f} {:>11.5f}\n".format(atomsymb[atom['na']],
atom['rx'], atom['ry'], atom['rz']))
fh2.write("{:>4s} {:>11.5f} {:>11.5f} {:>11.5f}\n".format(atomsymb[atom['na']],
atom['rx'], atom['ry'], atom['rz']))
fh.close()
fh2.close()