Commit cbc9190e authored by Samuel Tan's avatar Samuel Tan
Browse files

finished off efpIntGamess()

parent e7deacf4
......@@ -203,8 +203,11 @@ def read_file(xyz_f):
reads a file into list of strings per line
"""
#print(xyz_f)
with open(xyz_f) as f:
return [line.strip() for line in f.readlines()]
try:
with open(xyz_f) as f:
return [line.strip() for line in f.readlines()]
except FileNotFoundError:
print("FileNotFoundError: in read_file(), can't find ", xyz_f)
def scriptPath(f = None):
"""
......
......@@ -5,7 +5,7 @@
from chemF import *
import molDB as mDB
import collections as col
import re
import re, sys
# must ensure last step is atomPrint()
# otherwise the string mangling messes up the tabling space!!
......@@ -156,6 +156,87 @@ def efpGamess(curr_sys, tsk_d):
return out_f
def efpIntGamess(tsk_d):
""" generate GAMESS EFP interaction input file
by combining the .efp files
and using the input files to determine geometry
TODO: right now only 2 fragments, generalise to
multiple fragments in future
probably with inp as a comma separated list of paths
and efp the same
"""
with open(tsk_d["efp_1"]) as f1:
efpL1 = f1.readlines()
# remove ending newlines
efpL1 = [s.rstrip("\n") for s in efpL1]
with open(tsk_d["efp_2"]) as f2:
efpL2 = f2.readlines()
efpL2 = [s.rstrip("\n") for s in efpL2]
# use input files to determine geometry
with open(tsk_d["inp_1"]) as i1:
inpL1 = i1.readlines()
with open(tsk_d["inp_2"]) as i2:
inpL2 = i2.readlines()
if "template" in tsk_d:
tmp = read_file(tsk_d["template"])
else:
tmp = read_file(scriptPath("templates/interactefp.template"))
# get xyz data from input file, after the C1 line
geom1_line = pos_uniq(inList = inpL1, tag = "^C1$")
geom1 = inpL1[geom1_line + 1: geom1_line + 4]
geom2_line = pos_uniq(inList = inpL2, tag = "^C1$")
geom2 = inpL2[geom2_line + 1: geom2_line + 4]
# check that geom from inp matches efp
#efpLX in lines 5 to 8 looks like
#A01N a b c
#A02C a b c
#A03H a b a
# get the element symbol part
# and compare it with element symbol from the inp file
geom1 = [i.split() for i in geom1]
x1 = [i[0] for i in geom1]
y1 = [i.split()[0] for i in efpL1[5:8]]
if all(x1[i] == y1[i][3:] for i in range(3)):
warning("frag1 matches")
else:
print("frag1 doesn't match!")
sys.exit(3)
geom2 = [i.split() for i in geom2]
x2 = [i[0] for i in geom2]
y2 = [i.split()[0] for i in efpL2[5:8]]
if all(x2[i] == y2[i][3:] for i in range(3)):
warning("frag2 matches")
else:
print("frag2 doesn't match!")
sys.exit(3)
# put A01N in front of xyz geom
g1 = [i[2:] for i in geom1] # last 3, i.e. xyz data, discard elem_sym and atomic num
g2 = [i[2:] for i in geom2]
# add strings together
frag1_geom = [ [y1[i]] + g1[i] for i in range(3)]
frag2_geom = [ [y2[i]] + g2[i] for i in range(3)]
# add a space in front and newline at the end
frag1_geom = [" " + " ".join(frag1_geom[i]) for i in range(3)]
frag2_geom = [" " + " ".join(frag2_geom[i]) for i in range(3)]
out_f = tmp + [" COORD=CART"] + \
[" FRAGNAME=EFP1"] + frag1_geom + \
[" FRAGNAME=EFP2"] + frag2_geom + \
["$EFP1"] + efpL1[3:] + \
[" "] + \
["$EFP2"] + efpL2[3:]
# insert space in front of $ groups
out_f = [" " + line if line.startswith("$") else line for line in out_f]
return out_f
def FMOintEnGamess(curr_sys, tsk_d):
......
......@@ -19,8 +19,12 @@ def readTaskFile(tsk_f):
global taskDict
for line in f:
if not line.strip().startswith("#"):
l = line.split("=")
taskDict[l[0]] = l[1]
try:
l = line.split("=")
taskDict[l[0]] = l[1]
except IndexError:
print("Line without equal symbol (=) found in task file--do you have a blank line?")
sys.exit(2)
# ensure correct types
if "nfrag" in taskDict:
taskDict["nfrag"] = int(taskDict["nfrag"])
......@@ -78,7 +82,7 @@ def main():
readTaskFile(args.taskfile)
except FileNotFoundError:
readTaskFile(scriptPath("taskfiles/intEgamess.gen"))
print("task file cannot be found ", args[0])
print("task file cannot be found ", args.taskfile)
usage()
#sys.exit(2)
......@@ -91,34 +95,38 @@ def main():
if args.template:
taskDict["template"] = args.template
essentials = ["task", "using"]
if taskDict["task"] == "interactefp":
essentials = ["task", "inp_1", "inp_2", "efp_1", "efp_2"]
else:
essentials = ["task", "using"]
# check that taskDict contains all the information we need
i = [x for x in essentials if x not in taskDict]
if i: # i not an empty list
warning("Error: not specified in either task file or command line arguments: ", ", ".join(i))
print(taskDict)
usage()
return 1
#sys.exit(5)
# System class from chemF: read in xyz coordinates
curr_sys = System(read_xyz = taskDict["using"])
try:
curr_sys = System(read_xyz = taskDict["using"])
except FileNotFoundError:
print("FileNotFoundError: specified xyz file can't be found")
print("cant find ", taskDict["using"])
sys.exit(2)
# collect atoms into molecules
if "nfrag" in taskDict:
print("nfrag detected: ", taskDict["nfrag"])
curr_sys.aggregateMol(nfrag = taskDict["nfrag"])
for mol in curr_sys.molDict.values():
if not mDB.isAnion(mol, q = True) and not mDB.isCation(mol, q = True):
warning("WARNING! ",
col.Counter(mol.atomListAsElem_Sym()),
" is not found in the database")
if taskDict["task"] != "interactefp":
try:
curr_sys = System(read_xyz = taskDict["using"])
except FileNotFoundError:
print("FileNotFoundError: specified xyz file can't be found")
print("cant find ", taskDict["using"])
sys.exit(2)
# collect atoms into molecules
if "nfrag" in taskDict:
print("nfrag detected: ", taskDict["nfrag"])
curr_sys.aggregateMol(nfrag = taskDict["nfrag"])
for mol in curr_sys.molDict.values():
if not mDB.isAnion(mol, q = True) and not mDB.isCation(mol, q = True):
warning("WARNING! ",
col.Counter(mol.atomListAsElem_Sym()),
" is not found in the database")
# printing using functions in formats.py
......@@ -127,7 +135,6 @@ def main():
# Psi4 SAPT
if taskDict["task"] == "psi4sapt":
a = genPsi4(curr_sys = curr_sys, tsk_d = taskDict)
a = prettyPrint(a, taskDict.get("out_filename"))
try:
prettyPrint(a, taskDict["out_filename"])
except KeyError:
......@@ -169,6 +176,15 @@ def main():
return 0
# GAMESS EFP interaction
elif taskDict["task"] == "interactefp":
a = efpIntGamess(tsk_d = taskDict)
if "out_filename" in taskDict:
prettyPrint(a, out_filename = taskDict["out_filename"])
else:
prettyPrint(a)
return 0
# GAMESS interaction energy (FMO version)
elif taskDict["task"] == "gen_cluster_intE":
# generates multiple input files, one for the whole cluster
......
task=interactefp
inp_1=efp_input_0.inp
inp_2=efp_input_1.inp
efp_1=efp_input_0.efp
efp_2=efp_input_1.efp
mwords=2
memddi=1
template=/short/k96/apps/src/ppqc/templates/interactefp.template
$SYSTEM MWORDS=2 MEMDDI=1 $END
$CONTRL RUNTYP=ENERGY COORD=FRAGONLY
ICHARG=0 MULT=1
$END
$EFRAG
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment