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