Commit 323b7402 authored by Samuel Tan's avatar Samuel Tan
Browse files

major changes; FMO cluster almost done

parent 3834a1a0
......@@ -22,7 +22,7 @@ class Atom(object):
self.atom_id = atom_id
# fill in other info using imported module periodictable
if not atomic_num:
self.atomic_num = int(molDB.getAtomicNumber(self.elem_sym))
self.atomic_num = molDB.getAtomicNumber(self.elem_sym)
if not charge:
self.charge = 0
# self.ptabHandle = getattr(ptab, self.elem_sym)
......@@ -109,12 +109,16 @@ class Molecule(object):
class System(object):
""" usually from xyz file """
def __init__(self, read_file = None, v = False):
if read_file:
self.f = read_xyz(read_file)
def __init__(self, read_xyz = None, v = False):
# expect an xyz file, other formats not supported
if read_xyz:
self.f = read_file(read_xyz)
self.lineno = len(self.f) - 2 # no of lines in xyz file
self.f = [line.split() for line in self.f]
self.size = int(self.f[0][0]) # number of atoms
if self.size != self.lineno:
print("WARNING: number of lines (", self.lineno, ") and number of atoms (", self.size, ") do not match in xyz file")
if self.f[1]:
self.comment = str(self.f[1]) # comment line in xyz files
......@@ -123,6 +127,8 @@ class System(object):
self.molDict = {}
for i in range(self.size):
a = self.f[i+2]
if len(a) != 4:
print("WARNING: number of columns in xyz file is incorrect!")
# create Atom
self.molDict[i] = Atom(atom_id = i, elem_sym = a[0], coord = a[1:4])
# create molecule with key "i"; note the lists
......@@ -151,6 +157,9 @@ class System(object):
# line as assignment
minPair.sort(key = itemgetter(2))
mol1, mol2 = minPair[0][:2]
self.molDict[mol1].printMol()
self.molDict[mol2].printMol()
for key, atom in self.molDict[mol2].atomDict.items():
self.molDict[mol1].addAtom(key, atom)
......@@ -189,12 +198,26 @@ def euclid_Dist(a, b):
# separate function might be useful, e.g. to purely read file
# mainly used by the class System though
def read_xyz(xyz_f):
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()]
def scriptPath(f = None):
p = os.path.dirname(os.path.realpath(sys.argv[0]))
"""
determines the location of the script, as opposed
to the directory it is being run from, so that
paths to template files can be resolved
relative to ppqc's location
"""
if sys.argv[0]:
p = os.path.dirname(os.path.realpath(sys.argv[0]))
else:
p = os.path.dirname(os.path.realpath("."))
if not f:
return p
elif f: # f a string
......@@ -220,14 +243,11 @@ def insert_List(inList, tag, insertThis):
effectively replaces "tag" with "insertThis"
and returns the modified list
"""
if not isinstance(insertThis, list):
warning("TypeError: insertList(), arg insertThis must be a list")
if len(inList) > 1:
tag_occurence = pos_uniq(inList, tag)
if tag_occurence:
if isinstance(tag_occurence, int):
l_top = inList[:tag_occurence]
l_bot = inList[tag_occurence + 1:] # note that the tag position is gone
return l_top + insertThis + l_bot
......
task=gen-geod-inp
mwords=200
# no need memddi it seems
title=geodesic file
using=xyz_files/c1mim-br-p1.xyz
......@@ -19,9 +19,9 @@ def genPsi4(curr_sys, tsk_d, template = None):
return 1
if "template" in tsk_d:
out_f = read_xyz(template)
out_f = read_file(template)
else:
out_f = read_xyz(scriptPath("templates/psi4sapt.template"))
out_f = read_file(scriptPath("templates/psi4sapt.template"))
# replace anything mentioned in task file with
......@@ -83,14 +83,17 @@ def genPsi4(curr_sys, tsk_d, template = None):
### GAMESS formats ###
def optGamess(curr_sys, tsk_d):
# note also used for geodesic GAMESS
# only diff is default template file
# Fri 12 Feb 201/
if not curr_sys:
print("missing system arg in optGamess()")
return 1
if "template" in tsk_d:
out_f = read_xyz(tsk_d["template"])
out_f = read_file(tsk_d["template"])
else:
out_f = read_xyz(scriptPath("templates/optGamess.template"))
out_f = read_file(scriptPath("templates/optGamess.template"))
# replace anything mentioned in task file with
# "variable = value"
......@@ -117,9 +120,9 @@ def efpGamess(curr_sys, tsk_d):
return 1
if "template" in tsk_d:
tmp = read_xyz(tsk_d["template"])
tmp = read_file(tsk_d["template"])
else:
tmp = read_xyz(scriptPath("templates/efp.template"))
tmp = read_file(scriptPath("templates/efp.template"))
# replace anything mentioned in task file with
# "variable = value"
......@@ -154,16 +157,16 @@ def efpGamess(curr_sys, tsk_d):
return out_f
def ie_Gamess(curr_sys, tsk_d):
def FMOintEnGamess(curr_sys, tsk_d):
if not curr_sys:
print("missing system arg in ie_Gamess()")
return 1
if "template" in tsk_d:
tmp = read_xyz(tsk_d["template"])
tmp = read_file(tsk_d["template"])
else:
tmp = read_xyz(scriptPath("templates/fmo.template"))
tmp = read_file("templates/fmo.template")
# replacements
for key,val in tsk_d.items():
......@@ -182,9 +185,9 @@ def ie_Gamess(curr_sys, tsk_d):
# cluster part
if "cluster_template" in tsk_d:
frag_d["cluster"] = read_xyz(tsk_d["cluster_template"])
frag_d["cluster"] = read_file(tsk_d["cluster_template"])
else:
frag_d["cluster"] = read_xyz(scriptPath("templates/fmo_cluster.template"))
frag_d["cluster"] = read_file(scriptPath("templates/fmo_cluster.template"))
cluster_xyz = [] # xyz_data in one big list
charge_list = [] # to join later with ","
......@@ -203,7 +206,7 @@ def ie_Gamess(curr_sys, tsk_d):
if mDB.isAnion(mol):
frag_d[i] = mySed(frag_d[i], tag = "ICHARG", repl = "ICHARG=-1")
frag_d["anions"].append(i)
frag_d[i][comment_line] += " anion " + getAnion(mol)
frag_d[i][comment_line] += " anion " + mDB.getAnion(mol)
charge_list.append("-1")
line_no.append("0," + str(n) + ",-" + str(n + len(xyz)))
......@@ -213,7 +216,7 @@ def ie_Gamess(curr_sys, tsk_d):
elif mDB.isCation(mol):
frag_d[i] = mySed(frag_d[i], tag = "ICHARG", repl = "ICHARG=1")
frag_d["cations"].append(i)
frag_d[i][comment_line] += " cation " + getCation(mol)
frag_d[i][comment_line] += " cation " + mDB.getCation(mol)
charge_list.append("1")
line_no.append("0," + str(n) + ",-" + str(n + len(xyz)))
......@@ -241,13 +244,28 @@ def ie_Gamess(curr_sys, tsk_d):
# now only do we put everything together for cluster
charge_list = ",".join(charge_list)
frag_d["cluster"] = mySed(frag_d["cluster"], tag = "INDAT\\(1\\)", repl = "INDAT\\(1\\)=" + charge_list)
frag_d["cluster"] = mySed(frag_d["cluster"], tag = "INDAT(1)", repl = "INDAT(1)=" + charge_list)
# flatten first, [[ ... ]]
cluster_xyz = [i for j in cluster_xyz for i in j]
# add list of elements in $DATA
# get list of elements in cluster for each atom [C, H, H, F, H]
element_set = [a for m in curr_sys.molDict.values() for a in m.atomListAsElem_Sym()]
# use Counter() to create non-repeated list of elements [C,H,F]
element_set = [a for a in col.Counter(element_set).keys()]
# get atomic numbers
for i,n in enumerate(element_set):
element_set[i] = ' ' + str(n) + ' ' + str(mDB.getAtomicNumber(n))
# insert this in the cluster input file
frag_d["cluster"] = insert_List(frag_d["cluster"], tag = "element_set", insertThis = element_set)
# supply string for repl to mySed
frag_d["cluster"] = insert_List(frag_d["cluster"], tag = "xyz_data", insertThis = cluster_xyz)
frag_d["cluster"] = [" " + line if line.startswith("$") else line for line in frag_d["cluster"]]
# todo
# charge_list (which we randomly put in INDAT??)
# and INDAT, which I don't know the meaning of
return frag_d
task=gen_int_e
task=gen_cluster_intE
nfrag=4
sam=hello
mwords=200
guess_charges=true
using=xyz_files/c1mim-br-p1.xyz
using=xyz_files/rand.xyz
title=My_Awesome_Input_File
RUNTYP=ENERGY
template=templates/fmo.template
#!/usr/bin/env python3
#from chemF import warning
import chemF
import chemF as cf
import collections as col
# functions for accessing databases
......@@ -12,17 +12,28 @@ def isAnion(mol):
# checking done via Counter() from collections, no sorting required, duplicates included
#return next((key for key, anion in AnionDB.items()
# if col.Counter(a) == col.Counter(anion)), None)
isAni = False
for key, anion in AnionDB.items():
if col.Counter(a) == col.Counter(anion):
warning("Anion detected ", key)
return True
cf.warning("Anion detected ", key)
isAni = True
break
#return True
#else:
# print("no match ", col.Counter(a), col.Counter(anion))
# return False
return isAni
def isCation(mol):
a = mol.atomListAsElem_Sym()
isCat = False
for key, cation in CationDB.items():
if col.Counter(a) == col.Counter(cation):
warning("Cation detected ", key)
return True
cf.warning("Cation detected ", key)
#return True
isCat = True
break
return isCat
def getCation(mol):
a = mol.atomListAsElem_Sym()
......@@ -41,7 +52,7 @@ def getAtomicNumber(elem_sym):
""" function to return atomic number given chemical symbol """
for num, sym in periodic_table.items():
if elem_sym in sym:
return num
return float(num)
# molecule databases
AnionDB = {"br" : ["Br"]}
......
......@@ -14,14 +14,15 @@ def readTaskFile(tsk_f):
in a global dictionary, for all to see
ignores comment lines, which start with '#'
"""
f = read_xyz(tsk_f)
f = read_file(tsk_f)
global taskDict
for line in f:
if not line.strip().startswith("#"):
l = line.split("=")
taskDict[l[0]] = l[1]
# ensure correct types
taskDict["nfrag"] = int(taskDict["nfrag"])
if "nfrag" in taskDict:
taskDict["nfrag"] = int(taskDict["nfrag"])
if "guess_charges" in taskDict:
taskDict["guess_charges"] = taskDict["guess_charges"].lower()
......@@ -33,8 +34,14 @@ def prettyPrint(contentsList, out_filename = None):
if out_filename:
f = open(out_filename, 'w+')
f.write("\n".join(contentsList))
else:
elif not contentsList:
print(" empty list")
elif isinstance(contentsList[0], int):
print("\n".join(str(contentsList)))
elif isinstance(contentsList[0], str):
print("\n".join(contentsList))
else:
print("TypeError: prettyPrint() arg must be a list")
......@@ -84,9 +91,10 @@ def main():
and must be prefixed with a marker""")
usage()
else:
readTaskFile("intEgamess.gen")
print("Task file cannot be found: did you forget to specify it?")
usage()
sys.exit(2)
#sys.exit(2)
# override task file defaults with optional arguments
if "-i" in optDict:
......@@ -108,11 +116,12 @@ def main():
#sys.exit(5)
# System class from chemF
# System class from chemF: read in xyz coordinates
try:
curr_sys = System(read_file = taskDict["using"])
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
......@@ -161,23 +170,51 @@ def main():
return 0
# GAMESS interaction energy (FMO version)
elif taskDict["task"] == "int_e_gam":
warning("running ie_Gamess()")
a = ie_Gamess(curr_sys = curr_sys, tsk_d = taskDict)
elif taskDict["task"] == "gen_cluster_intE":
# generates multiple input files, one for the whole cluster
# and individual files for each cation and anion
a = FMOintEnGamess(curr_sys = curr_sys, tsk_d = taskDict)
#for j,f in a.items():
# print("this is ", j)
# prettyPrint(f)
#print(a)
# these print nicely now
# ironing out the last bits
# of the cluster input file
for cation in a["cations"]:
print("Cation ", cation)
prettyPrint(a[cation])
print('end cation')
for anion in a["anions"]:
print("Anion ", anion)
prettyPrint(a[anion])
# print cluster
prettyPrint(a["cluster"])
return 0
# geodesic GAMESS
elif taskDict["task"] == "gen-geod-inp":
key_list = [str(a).lower() for a in taskDict.keys()]
if [b for b in key_list if b in ["basis", "gbasis"]]:
print("No support for Basis group in GAMESS at the moment")
print("Please remove line from task file and use a template file instead")
return 0
if "template" in taskDict:
pass
else:
taskDict["template"] = scriptPath("templates/geodesic.template")
# note using same function as optimise, since
# basically just replacing xyz_data
a = optGamess(curr_sys = curr_sys, tsk_d = taskDict)
a = prettyPrint(a)
return 0
else:
print("Task not recognized")
print("Task not recognized", taskDict["task"])
return 1
#return a
......
$SYSTEM MWORDS=125 $END
$CONTRL SCFTYP=RHF RUNTYP=ENERGY ICHARG=0 MULT=1 ISPHER=1 $END
$BASIS GBASIS=ACCT $END
$SCF DIRSCF=.TRUE. DIIS=.TRUE. FDIFF=.FALSE. $END
$ELPOT IEPOT=1 WHERE=PDC $END
$ELMOM IEMOM=2 $END
$PDC PTSEL=GEODESIC $END
$DATA
title
C1
xyz_data
$END
42
C 0.032724238 0.098304735 0.334389491
H -0.840117910 0.250178044 -0.290535965
C 1.367710233 0.012249626 0.010236470
H 1.869831682 0.075730104 -0.948557842
N -0.052063379 -0.031577762 1.698957674
C 1.178348797 -0.176471785 2.202916287
H 1.424986752 -0.298531463 3.250403505
N 2.050018627 -0.166941754 1.188537483
C 3.497792332 -0.205008581 1.369009272
H 3.891198183 0.811465464 1.265604884
H 3.692692060 -0.576553135 2.380270448
H 3.932863121 -0.881375394 0.623757235
C -1.250230936 0.100591854 2.521778367
H -1.521490658 1.159715822 2.579418927
H -1.004911853 -0.273766167 3.520974874
H -2.055551565 -0.496668469 2.078041556
F 1.183610621 2.681418232 1.449040402
B 1.520063647 2.862627555 2.796653816
F 1.689891622 4.189110194 3.143234908
F 2.729317581 2.116983596 3.075608894
F 0.487416095 2.263212652 3.615529732
C 1.994171901 1.245126497 8.150396029
H 1.489676907 1.182008733 9.107966323
C 3.329899266 1.158547929 7.829451946
H 4.201146259 1.006445358 8.456531436
N 1.314776422 1.424527757 6.970419862
C 2.188914611 1.433686980 5.958148888
H 1.945122099 1.555614407 4.909976537
N 3.418031031 1.288342408 6.465106416
C 4.618028407 1.155181862 5.645126586
H 4.888664804 0.095860130 5.588304705
H 4.375362055 1.529640592 4.645325775
H 5.422824849 1.751786407 6.090708181
C -0.132555286 1.463624380 6.786545524
H -0.526953720 0.447482736 6.889454218
H -0.324795106 1.834828972 5.774641018
H -0.568714310 2.140891145 7.530334560
F 2.189380427 -1.424768457 6.710546157
B 1.851311545 -1.605759197 5.363319145
F 1.682165670 -2.932207363 5.016269838
F 2.882054652 -1.004750473 4.543264729
F 0.640788225 -0.861355171 5.086370708
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