Scripts/gamess.rb

   1 #
   2 #  gamess.rb
   3 #
   4 #  Created by Toshi Nagata on 2009/11/22.
   5 #  Copyright 2009 Toshi Nagata. All rights reserved.
   6 #
   7 # This program is free software; you can redistribute it and/or modify
   8 # it under the terms of the GNU General Public License as published by
   9 # the Free Software Foundation version 2 of the License.
  10 #
  11 # This program is distributed in the hope that it will be useful,
  12 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14 # GNU General Public License for more details.
  15
  16 class Molecule
  17
  18   def Molecule.read_gamess_basis_sets(fname)
  19     # $gamess_basis = Hash.new unless $gamess_basis
  20         $gamess_ecp = Hash.new unless $gamess_ecp
  21         # $gamess_basis_desc = Hash.new unless $gamess_basis_desc
  22         # $gamess_basis_keys = [] unless $gamess_basis_keys
  23         basename = File.basename(fname, ".*")
  24         keys = []
  25         descname = nil
  26     File.open(fname, "r") { |fp|
  27           while (s = fp.gets)
  28             if s =~ /^\s*!\s*/ && descname == nil
  29               #  Get the descriptive name from the first comment line
  30                   s = Regexp.last_match.post_match
  31                   if s =~ /  EMSL/
  32                     descname = Regexp.last_match.pre_match
  33                   else
  34                     descname = (s.split)[0]
  35                   end
  36                   $gamess_basis_desc[basename] = descname
  37                   next
  38                 end
  39                 ss, bas = (s.split)[0..1]     #  Tokens delimited by whitespaces
  40                 next if ss == nil || ss == ""
  41                 if ss == "$ECP"
  42                 #  Read the ECP section
  43                   ecp = $gamess_ecp[basename]
  44                   if !ecp
  45                     ecp = []
  46                         $gamess_ecp[basename] = ecp
  47                   end
  48                   keys << basename unless keys.include?(basename)
  49                   while (s = fp.gets)
  50                     break if s=~ /\$END/
  51                         ary = s.split  #  (PNAME, PTYPE, IZCORE, LMAX+1)
  52                         elem = ary[0].match(/\w{1,2}/).to_a[0]  #  ary[0] should begin with an element name
  53                         ap = Parameter.builtin.elements.find { |p| p.name.casecmp(elem) == 0 }
  54                         raise MolbyError, "the ECP definition does not begin with an element name: #{s}" if !ap
  55                         ecpdef = s
  56                         ln = 1
  57                         (0..Integer(ary[3])).each {
  58                           s = fp.gets
  59                           raise MolbyError, "the ECP definition ends unexpectedly at line #{ln} for: #{s}" if !s
  60                           ln += 1
  61                           ary = s.split  #  (NGPOT, comment)
  62                           ecpdef += s
  63                           (1..Integer(ary[0])).each {
  64                             s = fp.gets
  65                                 raise MolbyError, "the ECP definition ends unexpectedly at line #{ln} for: #{s}" if !s
  66                                 ln += 1
  67                                 ecpdef += s
  68                           }
  69                         }
  70                     ecp[ap.index] = ecpdef
  71                   end
  72                 elsif ss =~ /\W/
  73                   #  Comments or other unrecognizable lines
  74                   next
  75                 elsif (ap = Parameter.builtin.elements.find { |p| p.name.casecmp(ss) == 0 || p.fullname.casecmp(ss) == 0 })
  76                   #  Valid basis definition
  77                   if bas == nil || bas =~ /\W/
  78                     bas = basename
  79                   end
  80                   basis = $gamess_basis[bas]
  81                   if !basis
  82                     basis = []
  83                         $gamess_basis[bas] = basis
  84                   end
  85                   keys << bas unless keys.include?(bas)
  86                   basdef = ""
  87                   while (s = fp.gets) && s =~ /\S/
  88                     basdef += s
  89                   end
  90                   basis[ap.index] = basdef
  91                 else
  92                   raise MolbyError, "ss is not a valid element symbol or name: #{s}"
  93                 end
  94           end
  95     }
  96         unless $gamess_basis_keys.include?(basename)
  97           $gamess_basis_keys.push(basename)
  98         end
  99   end
 100
 101   #  Execute GAMESS (inferior copy of rungms script)
 102   #  inpname is the input file
 103   #  mol (optional) is the molecule from which the GAMESS input was built.
 104   #  If mol is specified and RUNTYP=OPTIMIZE, then the intermediate structures are
 105   #  displayed real-time.
 106   def Molecule.execute_gamess(inpname, mol = nil)
 107     gmsname = get_global_settings("gamess.executable_path")
 108     gmsdir = nil
 109     gmsvers = nil
 110     while 1
 111       if gmsname == nil || !File.exist?(gmsname)
 112         gmsname = Dialog.open_panel("Please locate the GAMESS executable")
 113         exit if gmsname == nil
 114       end
 115       gmsbase = File.basename(gmsname)
 116       gmsdir = File.dirname(gmsname)
 117       if gmsbase =~ /gamess\.(.*)\.(exe|x)$/i
 118         gmsvers = $1
 119         break
 120       else
 121         gmsname = nil
 122         error_message_box(gmsbase + " does not look like a GAMESS executable!")
 123       end
 124     end
 125
 126     inpbase = File.basename(inpname)
 127     inpdir = File.dirname(inpname)
 128     inpbody = inpbase.sub(/\.inp$/, "")
 129     logbase = inpbody + ".log"
 130
 131     set_global_settings("gamess.executable_path", gmsname)
 132
 133     ncpus = get_global_settings("gamess.ncpus").to_i
 134         if ncpus == 0
 135           ncpus = 1
 136         end
 137
 138     #  Prepare the scratch directory in the home directory
 139     #  (Not in the document home to avoid space-containing path in Windows)
 140     scrdir = document_home.sub(/\/My Documents/, "") + "/gamess"
 141     n = 0
 142     while File.exist?(scrdir) && !File.directory?(scrdir)
 143       if n == 0
 144         scrdir += ".1"
 145       else
 146         scrdir = scrdir.sub(".#{n}", ".#{n + 1}")
 147       end
 148       n += 1
 149     end
 150     if !File.exist?(scrdir)
 151       Dir.mkdir(scrdir)
 152     end
 153     scrdir = scrdir + "/" + inpbody + "." + $$.to_s + ".0"
 154     n = 0
 155     while File.exist?(scrdir)
 156       scrdir = scrdir.sub(".#{n}", ".#{n + 1}")
 157       n += 1
 158     end
 159     Dir.mkdir(scrdir)
 160
 161     if $platform == "win"
 162       sep = "\\"
 163       scrdir.gsub!("/", sep)
 164       gmsdir.gsub!("/", sep)
 165     else
 166       sep = "/"
 167     end
 168
 169     #  Get the host name etc.
 170     hostname = `hostname`.chomp
 171     if $platform == "win"
 172       freebytes = `dir #{scrdir}`.split("\n").pop.match(/([0-9,]+)[^0-9]*$/).to_a[1]
 173       if freebytes
 174         freebytes = (freebytes.gsub(",","").to_i / 1024).to_s + " Kbytes"
 175       else
 176         freebytes = "(unknown)"
 177       end
 178       uname = `ver`.to_s.gsub("\n", "")
 179     else
 180       freebytes = `df -k #{scrdir}`
 181       uname = `uname`.chomp
 182     end
 183
 184     #  Redirect standard output to the log file
 185     logname = scrdir + sep + logbase
 186     fpout = File.open(logname, "w")
 187     fpout.print "----- GAMESS execution script -----\n"
 188     fpout.print "This job is running on host #{hostname}\n"
 189     fpout.print "under operating system #{uname} at #{Time.now.to_s}\n"
 190     fpout.print "Available scratch disk space (Kbyte units) at beginning of the job is\n"
 191     fpout.print freebytes
 192
 193     #  Copy the input file
 194     scrprefix = scrdir + sep + inpbody
 195     filecopy(inpname, scrprefix + ".F05")
 196
 197     #  Prepare environmental variables
 198     auxdir = "#{gmsdir}#{sep}auxdata"
 199     ENV["ERICFMT"] = "#{auxdir}#{sep}ericfmt.dat"
 200     ENV["MCPPATH"] = "#{auxdir}#{sep}MCP"
 201     ENV["BASPATH"] = "#{auxdir}#{sep}BASES"
 202     ENV["QUANPOL"] = "#{auxdir}#{sep}QUANPOL"
 203     ENV["EXTBAS"] = "/dev/null"
 204     ENV["IRCDATA"] = "#{scrprefix}.irc"
 205     ENV["PUNCH"] = "#{scrprefix}.dat"
 206     ENV["INPUT"] = "#{scrprefix}.F05"
 207     ENV["AOINTS"] = "#{scrprefix}.F08"
 208     ENV["MOINTS"] = "#{scrprefix}.F09"
 209     ENV["DICTNRY"] = "#{scrprefix}.F10"
 210     ENV["DRTFILE"] = "#{scrprefix}.F11"
 211     ENV["CIVECTR"] = "#{scrprefix}.F12"
 212     ENV["CASINTS"] = "#{scrprefix}.F13"
 213     ENV["CIINTS"] = "#{scrprefix}.F14"
 214     ENV["WORK15"] = "#{scrprefix}.F15"
 215     ENV["WORK16"] = "#{scrprefix}.F16"
 216     ENV["CSFSAVE"] = "#{scrprefix}.F17"
 217     ENV["FOCKDER"] = "#{scrprefix}.F18"
 218     ENV["WORK19"] = "#{scrprefix}.F19"
 219     ENV["DASORT"] = "#{scrprefix}.F20"
 220     ENV["DFTINTS"] = "#{scrprefix}.F21"
 221     ENV["DFTGRID"] = "#{scrprefix}.F22"
 222     ENV["JKFILE"] = "#{scrprefix}.F23"
 223     ENV["ORDINT"] = "#{scrprefix}.F24"
 224     ENV["EFPIND"] = "#{scrprefix}.F25"
 225     ENV["PCMDATA"] = "#{scrprefix}.F26"
 226     ENV["PCMINTS"] = "#{scrprefix}.F27"
 227     ENV["MLTPL"] = "#{scrprefix}.F28"
 228     ENV["MLTPLT"] = "#{scrprefix}.F29"
 229     ENV["DAFL30"] = "#{scrprefix}.F30"
 230     ENV["SOINTX"] = "#{scrprefix}.F31"
 231     ENV["SOINTY"] = "#{scrprefix}.F32"
 232     ENV["SOINTZ"] = "#{scrprefix}.F33"
 233     ENV["SORESC"] = "#{scrprefix}.F34"
 234     ENV["SIMEN"] = "#{scrprefix}.simen"
 235     ENV["SIMCOR"] = "#{scrprefix}.simcor"
 236     ENV["GCILIST"] = "#{scrprefix}.F37"
 237     ENV["HESSIAN"] = "#{scrprefix}.F38"
 238     ENV["SOCCDAT"] = "#{scrprefix}.F40"
 239     ENV["AABB41"] = "#{scrprefix}.F41"
 240     ENV["BBAA42"] = "#{scrprefix}.F42"
 241     ENV["BBBB43"] = "#{scrprefix}.F43"
 242     ENV["MCQD50"] = "#{scrprefix}.F50"
 243     ENV["MCQD51"] = "#{scrprefix}.F51"
 244     ENV["MCQD52"] = "#{scrprefix}.F52"
 245     ENV["MCQD53"] = "#{scrprefix}.F53"
 246     ENV["MCQD54"] = "#{scrprefix}.F54"
 247     ENV["MCQD55"] = "#{scrprefix}.F55"
 248     ENV["MCQD56"] = "#{scrprefix}.F56"
 249     ENV["MCQD57"] = "#{scrprefix}.F57"
 250     ENV["MCQD58"] = "#{scrprefix}.F58"
 251     ENV["MCQD59"] = "#{scrprefix}.F59"
 252     ENV["MCQD60"] = "#{scrprefix}.F60"
 253     ENV["MCQD61"] = "#{scrprefix}.F61"
 254     ENV["MCQD62"] = "#{scrprefix}.F62"
 255     ENV["MCQD63"] = "#{scrprefix}.F63"
 256     ENV["MCQD64"] = "#{scrprefix}.F64"
 257     ENV["NMRINT1"] = "#{scrprefix}.F61"
 258     ENV["NMRINT2"] = "#{scrprefix}.F62"
 259     ENV["NMRINT3"] = "#{scrprefix}.F63"
 260     ENV["NMRINT4"] = "#{scrprefix}.F64"
 261     ENV["NMRINT5"] = "#{scrprefix}.F65"
 262     ENV["NMRINT6"] = "#{scrprefix}.F66"
 263     ENV["DCPHFH2"] = "#{scrprefix}.F67"
 264     ENV["DCPHF21"] = "#{scrprefix}.F68"
 265     ENV["GVVPT"] = "#{scrprefix}.F69"
 266
 267     #    next files are used only during coupled cluster runs, so let's
 268     #    display the numerous definitions only if they are to be used.
 269     ENV["CCREST"] = "#{scrprefix}.F70"
 270     ENV["CCDIIS"] = "#{scrprefix}.F71"
 271     ENV["CCINTS"] = "#{scrprefix}.F72"
 272     ENV["CCT1AMP"] = "#{scrprefix}.F73"
 273     ENV["CCT2AMP"] = "#{scrprefix}.F74"
 274     ENV["CCT3AMP"] = "#{scrprefix}.F75"
 275     ENV["CCVM"] = "#{scrprefix}.F76"
 276     ENV["CCVE"] = "#{scrprefix}.F77"
 277     ENV["EOMSTAR"] = "#{scrprefix}.F80"
 278     ENV["EOMVEC1"] = "#{scrprefix}.F81"
 279     ENV["EOMVEC2"] = "#{scrprefix}.F82"
 280     ENV["EOMHC1"] = "#{scrprefix}.F83"
 281     ENV["EOMHC2"] = "#{scrprefix}.F84"
 282     ENV["EOMHHHH"] = "#{scrprefix}.F85"
 283     ENV["EOMPPPP"] = "#{scrprefix}.F86"
 284     ENV["EOMRAMP"] = "#{scrprefix}.F87"
 285     ENV["EOMRTMP"] = "#{scrprefix}.F88"
 286     ENV["EOMDG12"] = "#{scrprefix}.F89"
 287     ENV["MMPP"] = "#{scrprefix}.F90"
 288     ENV["MMHPP"] = "#{scrprefix}.F91"
 289     ENV["MMCIVEC"] = "#{scrprefix}.F92"
 290     ENV["MMCIVC1"] = "#{scrprefix}.F93"
 291     ENV["MMCIITR"] = "#{scrprefix}.F94"
 292     ENV["MMNEXM"] = "#{scrprefix}.F95"
 293     ENV["MMNEXE"] = "#{scrprefix}.F96"
 294     ENV["MMNREXM"] = "#{scrprefix}.F97"
 295     ENV["MMNREXE"] = "#{scrprefix}.F98 "
 296     #
 297     #     next are for TDHFX code, not used by current GAMESS
 298     #
 299     ENV["OLI201"] = "#{scrprefix}.F201"
 300     ENV["OLI202"] = "#{scrprefix}.F202"
 301     ENV["OLI203"] = "#{scrprefix}.F203"
 302     ENV["OLI204"] = "#{scrprefix}.F204"
 303     ENV["OLI205"] = "#{scrprefix}.F205"
 304     ENV["OLI206"] = "#{scrprefix}.F206"
 305     ENV["OLI207"] = "#{scrprefix}.F207"
 306     ENV["OLI208"] = "#{scrprefix}.F208"
 307     ENV["OLI209"] = "#{scrprefix}.F209"
 308     ENV["OLI210"] = "#{scrprefix}.F210"
 309     ENV["OLI211"] = "#{scrprefix}.F211"
 310     ENV["OLI212"] = "#{scrprefix}.F212"
 311     ENV["OLI213"] = "#{scrprefix}.F213"
 312     ENV["OLI214"] = "#{scrprefix}.F214"
 313     ENV["OLI215"] = "#{scrprefix}.F215"
 314     ENV["OLI216"] = "#{scrprefix}.F216"
 315     ENV["OLI217"] = "#{scrprefix}.F217"
 316     ENV["OLI218"] = "#{scrprefix}.F218"
 317     ENV["OLI219"] = "#{scrprefix}.F219"
 318     ENV["OLI220"] = "#{scrprefix}.F220"
 319     ENV["OLI221"] = "#{scrprefix}.F221"
 320     ENV["OLI222"] = "#{scrprefix}.F222"
 321     ENV["OLI223"] = "#{scrprefix}.F223"
 322     ENV["OLI224"] = "#{scrprefix}.F224"
 323     ENV["OLI225"] = "#{scrprefix}.F225"
 324     ENV["OLI226"] = "#{scrprefix}.F226"
 325     ENV["OLI227"] = "#{scrprefix}.F227"
 326     ENV["OLI228"] = "#{scrprefix}.F228"
 327     ENV["OLI229"] = "#{scrprefix}.F229"
 328     ENV["OLI230"] = "#{scrprefix}.F230"
 329     ENV["OLI231"] = "#{scrprefix}.F231"
 330     ENV["OLI232"] = "#{scrprefix}.F232"
 331     ENV["OLI233"] = "#{scrprefix}.F233"
 332     ENV["OLI234"] = "#{scrprefix}.F234"
 333     ENV["OLI235"] = "#{scrprefix}.F235"
 334     ENV["OLI236"] = "#{scrprefix}.F236"
 335     ENV["OLI237"] = "#{scrprefix}.F237"
 336     ENV["OLI238"] = "#{scrprefix}.F238"
 337     ENV["OLI239"] = "#{scrprefix}.F239"
 338
 339     if $platform == "win"
 340           if ncpus < 2
 341             ncpus = 2
 342           end
 343       fpout.print "Microsoft MPI will be running GAMESS on 1 node.\n"
 344       fpout.print "The binary kicked off by 'mpiexec' is gamess.#{gmsvers}.exe\n"
 345       fpout.print "MS-MPI will run #{ncpus} compute process\n"
 346       #  File containing environmental variables
 347       envfil = "#{scrprefix}.GMS.ENV"
 348       fp = File.open(envfil, "w")
 349       ENV.each { |k, v| fp.print "#{k}=#{v}\n" }
 350       fp.close
 351       #  File containing arguments to mpiexec
 352       procfil = "#{scrprefix}.processes.mpd"
 353       fp = File.open(procfil, "w")
 354       fp.print "-env ENVFIL #{envfil} -n #{ncpus} #{gmsdir}#{sep}gamess.#{gmsvers}.exe\n"
 355       fp.close
 356     end
 357
 358     fpout.close
 359     fplog = File.open(logname, "r")
 360     size = 0
 361     lines = []
 362     last_line = ""
 363
 364     #  Callback proc
 365     callback = proc {
 366       fplog.seek(0, IO::SEEK_END)
 367       sizec = fplog.tell
 368       if sizec > size
 369         #  Read new lines
 370         fplog.seek(size, IO::SEEK_SET)
 371         fplog.each_line { |line|
 372           if line[-1, 1] == "\n"
 373             lines.push(last_line + line)
 374             last_line = ""
 375           else
 376             last_line += line
 377             break
 378           end
 379         }
 380         size = fplog.tell
 381         last_i = nil
 382         i = 0
 383         nserch = -1
 384         while i < lines.count
 385           line = lines[i]
 386           if line =~ /GEOMETRY SEARCH POINT NSERCH= *(\d+)/
 387             nserch = $1.to_i
 388             last_i = i
 389           elsif line =~ /NSERCH:/
 390             print line
 391                         if mol
 392                           dummy, n, grad = line.match(/NSERCH:[^0-9]*([0-9]+).*GRAD[^0-9]*([-.0-9]+)/).to_a
 393                           mol.show_text("Search: #{n}\nGradient: #{grad}")
 394                         end
 395                         last_i = i
 396           elsif nserch > 0 && line =~ /ddikick\.x/
 397             last_i = -1
 398             break
 399           elsif mol && nserch > 0 && line =~ /COORDINATES OF ALL ATOMS/
 400             #  There should be (natoms + 2) lines
 401             if i + mol.natoms + 3 <= lines.count
 402               coords = []
 403               (i + 3...i + 3 + mol.natoms).each { |j|
 404                 name, charge, x, y, z = lines[j].split
 405                 coords.push(Vector3D[x.to_f, y.to_f, z.to_f])
 406               }
 407               mol.create_frame([coords])
 408               mol.display
 409               last_i = i + mol.natoms + 2
 410               i = last_i   #  Skip the processed lines
 411             end
 412           end
 413           i += 1
 414         end
 415         if last_i == -1
 416           lines.clear
 417           break
 418         elsif last_i
 419           lines[0..last_i] = nil
 420         end
 421       end
 422       true
 423     }
 424
 425     show_console_window
 426     if mol
 427           mol.make_front
 428         end
 429
 430     if $platform == "win"
 431       status = call_subprocess("cmd.exe /c \"mpiexec -configfile #{procfil} >>#{logname}\"", "GAMESS", callback)
 432     else
 433           hosts = "localhost " * ncpus
 434       status = call_subprocess("/bin/sh -c '#{gmsdir}/ddikick.x #{gmsdir}/gamess.#{gmsvers}.x #{inpbody} -ddi #{ncpus} #{ncpus} #{hosts} -scr #{scrdir} < /dev/null >>#{logname}'", "GAMESS", callback)
 435     end
 436
 437     if status != 0
 438       if status == -1
 439         error_message_box("GAMESS failed to start. Please examine GAMESS installation.")
 440         exit
 441       elsif status == -2
 442         error_message_box("GAMESS execution interrupted.")
 443       else
 444         error_message_box("GAMESS failed with exit status #{status}.")
 445       end
 446         end
 447
 448         ext_to_keep = [".dat", ".rst", ".trj", ".efp", ".gamma", ".log"]
 449         ext_to_keep.each { |ex|
 450           if File.exists?("#{scrprefix}#{ex}")
 451                 filecopy("#{scrprefix}#{ex}", "#{inpdir}#{sep}#{inpbody}#{ex}")
 452           end
 453         }
 454         Dir.foreach(scrdir) { |file|
 455           if file != "." && file != ".." && !ext_to_keep.include?(File.extname(file))
 456                 File.delete("#{scrdir}#{sep}#{file}")
 457           end
 458         }
 459
 460   end
 461
 462   def export_gamess(fname, hash)
 463
 464     def reorder_array(ary, ordered_sub_ary)
 465           return ordered_sub_ary + (ary - ordered_sub_ary)
 466         end
 467
 468         now = Time.now.to_s
 469         basename = File.basename(fname, ".*")
 470
 471         #  Various settings
 472         icharg = hash["charge"]
 473         mult = hash["mult"]
 474         runtyp = ["ENERGY", "PROP", "OPTIMIZE"][Integer(hash["runtype"])]
 475         scftyp = ["RHF", "ROHF", "UHF"][Integer(hash["scftype"])]
 476         bssname = hash["basis"]
 477         bssname2 = hash["secondary_basis"]
 478         if hash["use_secondary_basis"] != 0 && bssname2 != bssname
 479           use_2nd = true
 480           element2 = hash["secondary_elements"].split(/[\s,]+/).map { |name| name.capitalize }
 481         else
 482           use_2nd = false
 483           bssname2 = bssname
 484         end
 485         basis = $gamess_basis[bssname]
 486         basis2 = $gamess_basis[bssname2]
 487         if !basis || !basis2
 488           raise MolbyError, "Unknown basis set name??? \"#{bssname}\" or \"#{bssname2}\""
 489         end
 490
 491         #  Use effective core potentials?
 492         ecp = $gamess_ecp[bssname]
 493         ecp2 = $gamess_ecp[bssname2]
 494         ecp_read = (ecp || ecp2 ? "ECP=READ " : "")
 495
 496         #  Use only one built-in basis set?
 497         gbasis = nil
 498         if !use_2nd
 499           case bssname
 500           when "PM3"
 501                 gbasis = "PM3 NGAUSS=3"
 502           when "STO3G"
 503                 gbasis = "STO NGAUSS=3"
 504           when "321G"
 505                 gbasis = "N21 NGAUSS=3"
 506           when "631G"
 507                 gbasis = "N31 NGAUSS=6"
 508           when "631Gd"
 509                 gbasis = "N31 NGAUSS=6 NDFUNC=1"
 510           when "631Gdp"
 511                 gbasis = "N31 NGAUSS=6 NDFUNC=1 NPFUNC=1"
 512           when "6311G"
 513                 gbasis = "N311 NGAUSS=6"
 514           when "6311Gdp"
 515                 gbasis = "N311 NGAUSS=6 NDFUNC=1 NPFUNC=1"
 516           end
 517         end
 518
 519         #  Count non-dummy atoms
 520         natoms = 0
 521         each_atom { |ap|
 522           natoms += 1 if ap.atomic_number != 0
 523         }
 524
 525         #  Fill hash with default values
 526         h = (hash["CONTRL"] ||= Hash.new)
 527         h["COORD"] ||= "UNIQUE"
 528         h["EXETYP"] ||= "RUN"
 529         h["ICHARG"] ||= icharg.to_s
 530         h["ICUT"] ||= "20"
 531         h["INTTYP"] ||= "HONDO"
 532         h["ITOL"] ||= "30"
 533         h["MAXIT"] ||= "200"
 534         h["MOLPLT"] ||= ".T."
 535         h["MPLEVL"] ||= "0"
 536         h["MULT"] ||= mult.to_s
 537         h["QMTTOL"] ||= "1e-08"
 538         h["RUNTYP"] ||= runtyp
 539         if hash["dft"] != 0
 540           h["DFTTYP"] ||= hash["dfttype"]
 541         end
 542         if hash["use_internal"] != 0 && (hash["runtype"] == 2 || h["RUNTYP"] == "OPTIMIZE")
 543           nzvar = natoms * 3 - 6  #  TODO: 3N-5 for linear molecules
 544           h["NZVAR"] = nzvar.to_s
 545         else
 546           nzvar = 0
 547         end
 548         h["SCFTYP"] ||= scftyp
 549         if ecp_read != ""
 550           h["ECP"] ||= "READ"
 551         end
 552         h["UNITS"] = "ANGS"
 553
 554         h = (hash["SCF"] ||= Hash.new)
 555         h["CONV"] ||= "1.0E-06"
 556         h["DIRSCF"] ||= ".T."
 557         h["FDIFF"] ||= ".T."
 558         h["DAMP"] ||= ".T."
 559
 560         h = (hash["STATPT"] ||= Hash.new)
 561         h["NSTEP"] ||= "400"
 562         h["OPTTOL"] ||= "1.0E-06"
 563
 564         h = (hash["SYSTEM"] ||= Hash.new)
 565         h["MEMDDI"] ||= "0"
 566         h["MWORDS"] ||= "16"
 567         h["TIMLIM"] ||= "50000"
 568
 569         h = (hash["GUESS"] ||= Hash.new)
 570         h["GUESS"] ||= "HUCKEL"
 571
 572         if gbasis
 573           h = (hash["BASIS"] ||= Hash.new)
 574           h["GBASIS"] ||= gbasis
 575         end
 576
 577         if nzvar > 0
 578           h = (hash["ZMAT"] ||= Hash.new)
 579           h["DLC"] ||= ".T."
 580           h["AUTO"] ||= ".T."
 581         end
 582
 583         if hash["esp"] != 0
 584           h = (hash["ELPOT"] ||= Hash.new)
 585           h["IEPOT"] ||= "1"
 586           h["OUTPUT"] ||= "PUNCH"
 587           h["WHERE"] ||= "PDC"
 588           h = (hash["PDC"] ||= Hash.new)
 589           h["CONSTR"] ||= "NONE"
 590           h["PTSEL"] ||= "CONNOLLY"
 591         end
 592
 593     File.open(fname, "wb") { |fp|
 594           fp.print "!  GAMESS input\n"
 595           fp.print "!  Generated by Molby at #{now}\n"
 596           fp.print "!  Basis set: " + $gamess_basis_desc[bssname] + "\n"
 597           if use_2nd
 598             fp.print "!  [" + element2.join(", ") + "]: " + $gamess_basis_desc[bssname2] + "\n"
 599           end
 600           controls = reorder_array(hash.keys.select { |k| hash[k].is_a?(Hash) },
 601                 ["CONTRL", "SCF", "STATPT", "SYSTEM", "GUESS", "BASIS", "ZMAT", "ELPOT", "PDC"])
 602           controls.each { |k|
 603             h = hash[k]
 604                 next if h == nil || h.size == 0
 605                 if k == "CONTRL"
 606                   ordered = ["COORD", "EXETYP", "ICHARG", "ICUT", "INTTYP", "ITOL", "MAXIT", "MOLPLT", "MPLEVL",
 607                              "MULT", "QMTTOL", "RUNTYP", "SCFTYP", "ECP", "UNITS", "DFTTYP", "NZVAR"]
 608                 elsif k == "SCF"
 609                   ordered = ["CONV", "DIRSCF", "FDIFF", "DAMP"]
 610                 elsif k == "STATPT"
 611                   ordered = ["NSTEP", "OPTTOL"]
 612                 elsif k == "SYSTEM"
 613                   ordered = ["MEMDDI", "MWORDS", "TIMLIM"]
 614                 elsif k == "GUESS"
 615                   ordered = ["GUESS"]
 616                 elsif k == "BASIS"
 617                   ordered = ["GBASIS"]
 618                 elsif k == "ZMAT"
 619                   ordered = ["DLC", "AUTO"]
 620                 elsif k == "ELPOT"
 621                   ordered = ["IEPOT", "OUTPUT", "WHERE"]
 622                 elsif k == "PDC"
 623                   ordered = ["CONSTR", "PTSEL"]
 624                 else
 625                   ordered = []
 626                 end
 627             keys = reorder_array(h.keys, ordered)
 628                 n = 0
 629                 keys.each_with_index { |kk, i|
 630                   v = h[kk]
 631                   next if v == nil || v == ""
 632                   if n == 0
 633                     fp.printf " $%-6s", k
 634                   elsif n % 3 == 0
 635                     fp.print "\n        "
 636                   end
 637                   fp.printf " %s=%s", kk, h[kk].to_s
 638                   n += 1
 639                 }
 640                 if n > 0
 641                   fp.print " $END\n"
 642                 end
 643           }
 644           fp.print " $DATA\n#{basename}\nC1 0\n"
 645           secondary = []
 646           each_atom { |ap|
 647             next if ap.atomic_number == 0
 648                 fp.printf "%-6s %4d %10.6f %10.6f %10.6f\n", ap.name, ap.atomic_number, ap.r.x, ap.r.y, ap.r.z
 649                 if use_2nd && element2.include?(ap.element)
 650                   secondary[ap.index] = true
 651                 end
 652             if !gbasis
 653                   #  Basis specification followed by a blank line
 654                   bas = (secondary[ap.index] ? basis2 : basis)
 655                   if bas.is_a?(Array)
 656                     bas = bas[ap.atomic_number]
 657                   end
 658                   if !bas
 659                     raise MolbyError, "Basis set is not defined for atom #{ap.index}, element #{ap.element}"
 660                   end
 661                   fp.print bas
 662                   fp.print "\n"
 663                 end
 664           }
 665           fp.print " $END\n"
 666           ecp_ary = []
 667           if ecp || ecp2
 668             fp.print " $ECP\n"
 669                 each_atom { |ap|
 670                   an = ap.atomic_number
 671                   next if an == 0
 672                   ecpp = (secondary[ap.index] ? ecp2 : ecp)
 673                   e = ecp_ary[an] || (ecpp && ecpp[an])
 674                   if e
 675                     #  Cache the PNAME of the $ECP entry and re-use it
 676                     ecp_ary[an] ||= (e.split)[0] + "\n"
 677                   else
 678                     e = ap.element.upcase + "-ECP NONE\n"
 679                   end
 680                   fp.print e
 681                 }
 682             fp.print " $END\n"
 683           end
 684         }
 685         fname
 686   end
 687
 688   def cmd_create_gamess_input
 689     if natoms == 0
 690       raise MolbyError, "cannot create GAMESS input; the molecule is empty"
 691     end
 692
 693         #  Descriptive text and internal string for popup menus
 694 #    bset_desc = ["PM3", "STO-3G", "3-21G", "6-31G", "6-31G(d)", "6-31G(d,p)", "6-311G", "6-311G(d,p)", "LanL2DZ"]
 695 #       bset_internal = ["PM3", "STO3G", "321G", "631G", "631Gd", "631Gdp", "6311G", "6311Gdp", "LanL2DZ"]
 696     bset_desc = $gamess_basis_keys.map { |key| $gamess_basis_desc[key] }
 697         dft_desc = ["B3LYP"]
 698         dft_internal = ["B3LYP"]
 699
 700         defaults = {"scftype"=>0, "runtype"=>0, "charge"=>"0", "mult"=>"1",
 701           "basis"=>4, "use_secondary_basis"=>0, "secondary_elements"=>"",
 702           "secondary_basis"=>8, "esp"=>0, "ncpus"=>"1"}
 703
 704     hash = Dialog.run("GAMESS Export") {
 705       def load_basis_set_sub(item)
 706             fname = Dialog.open_panel("Select a file containing GAMESS basis set:")
 707                 if fname
 708                   Molecule.read_gamess_basis_sets(fname)
 709                   bset_desc_new = $gamess_basis_keys.map { |key| $gamess_basis_desc[key] }
 710                   sel1 = attr("basis", :value)
 711                   sel2 = attr("secondary_basis", :value)
 712                   set_attr("basis", :subitems=>bset_desc_new)
 713                   set_attr("basis", :value=>sel1)
 714                   set_attr("secondary_basis", :subitems=>bset_desc_new)
 715                   set_attr("secondary_basis", :value=>sel2)
 716                 end
 717       end
 718           def select_gamess_path(item)
 719             while 1
 720               fname = Dialog.open_panel("Locate GAMESS executable:")
 721                   return if fname == nil
 722                   bname = File.basename(fname)
 723                   if bname =~ /gamess\.(.*)\.(exe|x)$/i
 724                     set_attr("executable_path", :value=>fname)
 725                     return
 726                   else
 727                     error_message_box("\"#{bname}\" does not look like a GAMESS executable!  Please try again.")
 728                   end
 729                 end
 730           end
 731           layout(4,
 732                 item(:text, :title=>"SCF type"),
 733                 item(:popup, :subitems=>["RHF", "ROHF", "UHF"], :tag=>"scftype"),
 734             item(:text, :title=>"Run type"),
 735                 item(:popup, :subitems=>["Energy", "Property", "Optimize"], :tag=>"runtype",
 736                   :action=>proc { |it| set_attr("use_internal", :enabled=>(it[:value] == 2)) } ),
 737
 738                 item(:checkbox, :title=>"Use internal coordinates for structure optimization", :tag=>"use_internal"),
 739                 -1, -1, -1,
 740
 741                 item(:text, :title=>"Charge"),
 742                 item(:textfield, :width=>80, :tag=>"charge"),
 743                 item(:text, :title=>"Multiplicity"),
 744                 item(:textfield, :width=>80, :tag=>"mult"),
 745
 746                 item(:checkbox, :title=>"Use DFT", :tag=>"dft",
 747                   :action=>proc { |it| set_attr("dfttype", :enabled=>(it[:value] != 0)) } ),
 748                 -1,
 749                 item(:text, :title=>"DFT type"),
 750                 item(:popup, :subitems=>dft_desc, :tag=>"dfttype"),
 751
 752                 item(:line),
 753                 -1, -1, -1,
 754
 755                 item(:text, :title=>"Basis set"),
 756                 item(:popup, :subitems=>bset_desc, :tag=>"basis"),
 757                 -1,
 758                 -1,
 759
 760                 item(:button, :title=>"Load Basis Set...", :action=>:load_basis_set_sub),
 761                 -1, -1, -1,
 762
 763                 item(:checkbox, :title=>"Use secondary basis set", :tag=>"use_secondary_basis",
 764                   :action=>proc { |it|
 765                     flag = (it[:value] != 0)
 766                         set_attr("secondary_elements", :enabled=>flag)
 767                         set_attr("secondary_basis", :enabled=>flag)
 768                   }),
 769                 -1, -1, -1,
 770
 771                 item(:text, :title=>"   Elements"),
 772                 item(:textfield, :width=>80, :tag=>"secondary_elements"),
 773                 item(:text, :title=>"Basis set"),
 774                 item(:popup, :subitems=>bset_desc, :tag=>"secondary_basis"),
 775
 776                 item(:line),
 777                 -1, -1, -1,
 778
 779                 item(:checkbox, :title=>"Calculate electrostatic potential (ESP)", :tag=>"esp"),
 780                 -1, -1, -1,
 781
 782                 item(:line),
 783                 -1, -1, -1,
 784
 785                 item(:checkbox, :title=>"Execute GAMESS on this machine", :tag=>"execute_local",
 786                   :action=>proc { |it|
 787                     flag = (it[:value] != 0)
 788                         set_attr("executable_path", :enabled=>flag)
 789                         set_attr("select_path", :enabled=>flag)
 790                         set_attr("ncpus", :enabled=>flag)
 791                   }),
 792                 -1, -1, -1,
 793
 794                 item(:text, :title=>"   Path"),
 795                 item(:textfield, :width=>300, :tag=>"executable_path"),
 796                 -1, -1,
 797
 798                 -1,
 799                 item(:button, :title=>"Select Path...", :tag=>"select_path", :action=>:select_gamess_path),
 800                 -1, -1,
 801
 802                 item(:text, :title=>"   N of CPUs"),
 803                 item(:textfield, :width=>80, :tag=>"ncpus"),
 804                 -1, -1,
 805
 806                 item(:line),
 807                 -1, -1, -1
 808           )
 809           values = Hash.new
 810           each_item { |it|
 811             tag = it[:tag]
 812                 if tag
 813                   values[tag] = (get_global_settings("gamess.#{tag}") || defaults[tag])
 814                   it[:value] = values[tag]
 815                 end
 816           }
 817           set_attr("secondary_elements", :enabled=>(values["use_secondary_basis"] == 1))
 818           set_attr("secondary_basis", :enabled=>(values["use_secondary_basis"] == 1))
 819           set_attr("dfttype", :enabled=>(values["dft"] == 1))
 820           set_attr("use_internal", :enabled=>(values["runtype"] == 2))
 821           set_attr("executable_path", :enabled=>(values["execute_local"] == 1))
 822           set_attr("select_path", :enabled=>(values["execute_local"] == 1))
 823           set_attr("ncpus", :enabled=>(values["execute_local"] == 1))
 824         }
 825         hash.each_pair { |key, value|
 826           next if key == :status
 827           set_global_settings("gamess.#{key}", value)
 828         }
 829         if hash[:status] == 0
 830           #  Specify basis by internal keys
 831           hash["basis"] = $gamess_basis_keys[hash["basis"]]
 832           hash["secondary_basis"] = $gamess_basis_keys[hash["secondary_basis"]]
 833           hash["dfttype"] = dft_internal[hash["dfttype"]]
 834           basename = (self.path ? File.basename(self.path, ".*") : self.name)
 835       fname = Dialog.save_panel("Export GAMESS input file:", self.dir, basename + ".inp", "GAMESS input file (*.inp)|*.inp|All files|*.*")
 836           return nil if !fname
 837           export_gamess(fname, hash)
 838           if hash["execute_local"] == 1
 839             Molecule.execute_gamess(fname, self)
 840           end
 841         else
 842           nil
 843         end
 844   end
 845
 846 end
 847
 848 $gamess_basis = {
 849   "PM3"   => " PM3 0\n",
 850   "STO3G" => " STO 3\n",
 851   "321G"  => " N21 3\n",
 852   "631G"  => " N31 6\n" }
 853 $gamess_basis_desc = {
 854   "PM3"   => "PM3",
 855   "STO3G" => "STO-3G",
 856   "321G"  => "3-21G",
 857   "631G"  => "6-31G" }
 858 $gamess_basis_keys = ["PM3", "STO3G", "321G", "631G"]
 859
 860 ["631Gd", "631Gdp", "631+Gd", "631++Gdp", "6311Gdp", "6311+Gd", "6311++Gdp", "6311++G2d2p", "6311++G3df3pd", "LanL2DZ"].each { |n|
 861   Molecule.read_gamess_basis_sets("basis_sets/#{n}.txt")
 862 }