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