path: root/sim/src/parameters/minimize.py

#! /usr/bin/python

from Numeric import *
from VQT import *
from LinearAlgebra import *
from string import *
import re
import os
import sys
from struct import pack
from stat import *

datapat = re.compile("\ \$DATA")
gradpat = re.compile("\ \$GRAD")
hesspat = re.compile(" \$HESS")
vecpat = re.compile(" \$VEC")
endpat = re.compile("\ \$END")
failpat = re.compile("-ABNORMALLY-")
blankpat = re.compile("\s*$")
equilpat = re.compile("1     \*\*\*\*\* EQUILIBRIUM GEOMETRY LOCATED \*\*\*\*\*")
coordpat = re.compile(" COORDINATES")
erecpat = re.compile("E\=\ +([\d\.-]+)\ +GMAX\=\ +([\d\.-]+)\ +GRMS\=\ +([\d\.-]+)")
frecpat = re.compile(".+\ +\d+\.\ +([\d\.E+-]+)\ +([\d\.E+-]+)\ +([\d\.E+-]+)")
irecpat = re.compile(" (\w+) +\d+\.\d* +([\d\.E+-]+) +([\d\.E+-]+) +([\d\.E+-]+)")
xyzpat = re.compile("(\w+) +([\d\.E+-]+) +([\d\.E+-]+) +([\d\.E+-]+)")
pdbname= re.compile("([^\.]+)\.pdb")
xyzname= re.compile("([^\.]+)\.xyz")
inpname= re.compile("([^\.]+)\.inp")

am1p = re.compile("GBASIS=AM1")

preface=["""! Gamess control file
 $CONTRL SCFTYP=RHF MAXIT=200 RUNTYP=optimize MULT=1 ICHARG=0 
ICUT=9 ITOL=20 INTTYP=best NPRINT=-5 $END     
 $SCF NCONV=8 dirscf=.t. DAMP=.t. SHIFT=.t. DIIS=.T. SOSCF=.T. npunch=0 $END 
 $STATPT NSTEP=50 OPTTOL=0.0001 $END     
 $FORCE VIBANL=.f. PRTIFC=.f. METHOD=analytic $END 
 $SYSTEM TIMLIM=10000 MWORDS=250 $END                                       
""",
        """! Gamess control file
 $CONTRL SCFTYP=RHF MAXIT=500 RUNTYP=energy MULT=1 ICHARG=0 
ICUT=9 ITOL=20 INTTYP=best NPRINT=-5 $END     
 $SCF NCONV=8 dirscf=.t. DAMP=.t. SHIFT=.t. DIIS=.T. SOSCF=.T. npunch=0 $END 
 $FORCE VIBANL=.f. PRTIFC=.f. METHOD=analytic $END 
 $SYSTEM TIMLIM=10000 MWORDS=250 $END                                       
""",
        """! Gamess control file
 $CONTRL SCFTYP=RHF MAXIT=500 RUNTYP=optimize MULT=1 ICHARG=0 
ICUT=9 ITOL=20 INTTYP=best NPRINT=-5 $END     
 $SCF NCONV=8 dirscf=.t. DAMP=.t. SHIFT=.t. DIIS=.T. SOSCF=.T. npunch=0 $END 
 $STATPT NSTEP=50 OPTTOL=0.0001 HSSEND=.true. $END     
 $FORCE VIBANL=.f. PRTIFC=.f. METHOD=seminum $END 
 $SYSTEM TIMLIM=10000 MWORDS=250 $END                                       
"""]
preface2=" $DATA\nGamess theory ladder: "
preface3="\nC1\n"
postface=""" $END
"""
dftface=""" $DFT DFTTYP=B3LYP $END
"""
levels=[" $BASIS GBASIS=AM1 $END\n",
        " $BASIS GBASIS=N21 NGAUSS=3 $END\n",
        " $BASIS GBASIS=N31 NGAUSS=6 $END\n",
        " $BASIS GBASIS=N31 NGAUSS=6 NDFUNC=1 NPFUNC=1 $END\n",
        " $BASIS GBASIS=N31 NGAUSS=6 NDFUNC=1 NPFUNC=1 DIFFSP=.t. $END\n",
        " $BASIS GBASIS=N311 NGAUSS=6 NDFUNC=2 NPFUNC=2 DIFFSP=.t. $END\n"
        ]
schedule=[#[0, 0, 0, " am1, no dft"], 
          [2, 0, 0, " am1, hessian"], 
          [0, 1, 1, " 3-21G, b3lyp"]]
##           [0, 2, 1, " 6-31G, b3lyp"], 
##           [0, 3, 1, " 6-31G(d,p), b3lyp"], 
##           [0, 3, 1, " 6-31G(d,p), b3lyp"], 
##           [0, 4, 1, " 6-31+G(d,p), b3lyp"], 
##           [1, 5, 1, " 6-31+G(2d,2p), b3lyp, single point energy"]] 
          
    
# distances are in angstroms, and delta T is 1e-16 seconds

# gradients from Gamess are given in Hartrees/Bohr. To convert to N:
gradU = 82.38729477e-9

# given in 1e-27 kg
elmnts=[("H",   1,   1.6737),
        ("He",  2,   6.646),
        ("Li",  3,  11.525),
        ("Be",  4,  14.964),
        ("B",   5,  17.949),
        ("C",   6,  19.925),
        ("N",   7,  23.257),
        ("O",   8,  26.565),
        ("F",   9,  31.545),
        ("Ne", 10,  33.49),

        ("Na", 11,  38.1726),
        ("Mg", 12,  40.356),
        ("Al", 13,  44.7997),
        ("Si", 14,  46.6245),
        ("P",  15,  51.429),
        ("S",  16,  53.233),
        ("Cl", 17,  58.867),
        ("Ar", 18,  66.33),
        
        ("K",  19,  64.9256),
        ("Ca", 20,  66.5495),
        ("Sc", 21,  74.646),
        ("Ti", 22,  79.534),
        ("V",  23,  84.584),
        ("Cr", 24,  86.335),
        ("Mn", 25,  91.22),
        ("Fe", 26,  92.729),
        ("Co", 27,  97.854),
        ("Ni", 28,  97.483),
        ("Cu", 29, 105.513),
        ("Zn", 30, 108.541),
        ("Ga", 31, 115.764),
        ("Ge", 32, 120.53),
        ("As", 33, 124.401),
        ("Se", 34, 131.106),
        ("Br", 35, 132.674),
        ("Kr", 36, 134.429)]


sym2num={}
num2mass=[0]
for (sym, num, mass) in elmnts:
    sym2num[sym] = num
    num2mass += [mass]
    
def datread(fname):
    f=open(fname,"r")

    grads=zeros((0,3),Float)

    while 1:
        card = f.readline()
        if gradpat.match(card): break

    f.readline() # R= ...

    while 1:
        card = f.readline()
        if endpat.match(card): break
        m=frecpat.match(card)
        v=A([map(float,(m.group(1),m.group(2), m.group(3)))])
        grads = concatenate((grads,v),axis=0)

    f.close()
    return grads

def inpread(fname):

    f=open(fname+'.inp',"r")

    atoms = zeros((0,3),Float)
    elem = []
    preface=" $DATA\n"
    postface=""

    while 1:
        card = f.readline()
        if datapat.match(card): break

    preface += f.readline() # the comment line
    preface += f.readline() # the symmetry line

    while 1:
        card = f.readline()
        if blankpat.match(card):
            preface += card # if non-C1 sym
            continue
        if endpat.match(card): break
        m=irecpat.match(card)
        elem += [m.group(1)]

        v=A([map(float,(m.group(2),m.group(3), m.group(4)))])
        atoms = concatenate((atoms,v),axis=0)

    postface = card
    f.close()
    return elem, atoms, preface
    
def inpwrite(fname, elem, pos, pref=" $DATA\ncomment\nC1\n"):
    f=open(fname+'.inp',"w")
    f.write(pref)

    for i in range(len(elem)):
        f.write(" %-10s %2d." % (elem[i], sym2num[elem[i]]) + 
                " %12.7f %12.7f %12.7f" % tuple(pos[i]) + "\n")
        
    f.write(" $END\n")
    f.close()

def xyzread(f):
    n=int(f.readline())
    elems=[]
    atoms = zeros((n,3),Float)
    f.readline() # skip comment line
    for i in range(n):
        m=xyzpat.match(f.readline())
        elems += [capitalize(m.group(1))]
        atoms[i]=A(map(float,(m.group(2),m.group(3), m.group(4))))
    return elems, atoms

def xyzwrite(f, elem, pos):
    f.write(str(len(elem))+"\n--comment line--\n")
    for i in range(len(elem)):
        f.write(elem[i] + " %12.7f %12.7f %12.7f" % tuple(pos[i]) + "\n")

def logread(name):
    f=open(name+'.log',"r")
    while 1:
        card = f.readline()
        if failpat.search(card):
            print 'GAMESS bombs in',name
            return None, None
        if card == "1     ***** EQUILIBRIUM GEOMETRY LOCATED *****\n": break
        if card == " **** THE GEOMETRY SEARCH IS NOT CONVERGED! ****":
            print 'GAMESS bombs in',name
            return None, None

    f.readline() # COORDINATES ...
    f.readline() #    ATOM   CHARGE       X              Y              Z
    f.readline() #  ------------------------------------------------------------

    atoms = zeros((0,3),Float)
    elem = []
    
    while 1:
        card = f.readline()
        if len(card)<10: break
        if coordpat.match(card): break
        m=irecpat.match(card)
        elem += [capitalize(m.group(1))]

        v=A([map(float,(m.group(2),m.group(3), m.group(4)))])
        atoms = concatenate((atoms,v),axis=0)

    return elem, atoms


def dpbwrite(f, pos):
    global ipos
    npos=floor(pos*100)
    delta=npos-ipos
    ipos = npos
    for line in delta:
        f.write(pack("bbb",int(line[0]),int(line[1]),int(line[2])))

def rungms(name, elem, pos, pref):
    print "*****  running",name,"   *****"
    inpwrite(name, elem, pos, pref)
    if am1p.search(pref):
        os.system("rungms "+name+" 1 > "+name+".log")
    else:
        os.system("rcp "+name+".inp cluster:/home/josh/working.inp")
        os.system("rsh cluster /home/gamess32/rungms working Nov222004R1 4 >"+name+".log")
    return name

def makexyz():
    print "making xyz files:"
    xyzs = {}
    files = os.listdir('level0')
    for file in files:
        x = xyzname.match(file)
        if x: xyzs[x.group(1)] = True
    for file in files:
        x = pdbname.match(file)
        if x and not x.group(1) in xyzs:
            print 'Making',x.group(1)+'.xyz'
            os.system("pdb2xyz " + 'level0/'+x.group(1))

def fexist(fname):
    try: os.stat(fname)
    except OSError: return False
    return True
            
    
def dirsetup():
    print "making input files:"
    files = os.listdir('level0')
    for file in files:
        x = xyzname.match(file)
        if x:
            name = x.group(1)
            if fexist('level0/'+name+'.inp'): continue
            elems, atoms = xyzread(open('level0/'+name+'.xyz'))
            inpwrite('level0/'+name, elems, atoms)

def dorunrun(n):
    print "\nrunning Gamess at level",n
    fd = 'level'+str(n-1)
    td = 'level'+str(n)
    header = open(td+'/theory').read()
    files = os.listdir(fd)
    for file in files:
        x = inpname.match(file)
        if not x: continue
        name = x.group(1)
        fname = fd+'/'+name
        tname = td+'/'+name
        if fexist(tname+'.log'): continue
        elem, atoms, pref = inpread(fname)
        if fexist(fname+'.log'): elem, atoms = logread(fname)
        if elem: rungms(tname, elem, atoms, header+pref)

            

if __name__ == "__main__":
    makexyz()
    dirsetup()
    n=1
    while 1:
        dorunrun(n)
        n += 1
        if not fexist('level'+str(n)): break