Scripts/gamess.rb

   1 # coding: utf-8
   2 #
   3 #  gamess.rb
   4 #
   5 #  Created by Toshi Nagata on 2009/11/22.
   6 #  Copyright 2009 Toshi Nagata. All rights reserved.
   7 #
   8 # This program is free software; you can redistribute it and/or modify
   9 # it under the terms of the GNU General Public License as published by
  10 # the Free Software Foundation version 2 of the License.
  11 #
  12 # This program is distributed in the hope that it will be useful,
  13 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15 # GNU General Public License for more details.
  16
  17 class MemoryIO < IO
  18   def initialize
  19     @buffer = ""
  20   end
  21   def write(s)
  22         @buffer << s
  23   end
  24   def buffer
  25     @buffer
  26   end
  27   def empty
  28     @buffer = ""
  29   end
  30 end
  31
  32 class Molecule
  33
  34   #  Import nbo log.  The argument is either an array of Strings
  35   #  or an IO object.
  36   def import_nbo_log(lines)
  37     if lines.is_a?(IO)
  38           getline = lambda { lines.gets }
  39         else
  40           getline = lambda { lines.shift }
  41     end
  42     nbo = Hash.new
  43         while (ln = getline.call) != nil
  44           if ln =~ /done with NBO analysis/
  45             break
  46           end
  47           if ln =~ /NATURAL POPULATIONS:  Natural atomic orbital occupancies/
  48             #  List of Natural AOs
  49                 getline.call
  50                 getline.call
  51                 getline.call
  52                 nbo["nao"] = []
  53                 while (ln = getline.call) != nil
  54                   if ln =~ /^\s*$/
  55                     #  Skip a blank line
  56                     ln = getline.call
  57                         if ln =~ /^\s*$/
  58                           #  Double blank lines indicates the end of the table
  59                           break
  60                         end
  61                   end
  62                   ln.chomp!
  63                   next unless ln =~ /^ *(\d+) *([A-Za-z]+) *(\d+) *([a-z]+) +([A-Za-z]+)\( *(\d+)([a-z]+)\)/
  64                   i = $1.to_i - 1
  65                   an = $3.to_i - 1
  66                   atom_label = $2 + $3.to_s
  67                   label = $4
  68                   type = $5
  69                   pn = $6
  70                   vals = $~.post_match.split.map { |e| e.to_f }
  71                   #  atom_index, type, pn+label, occupancy[, energy]
  72                   nbo["nao"].push([atom_label, type, pn + label] + vals)
  73                 end
  74           elsif ln =~ / NATURAL BOND ORBITALS \(Summary\):/
  75             nbo["nbo"] = []
  76                 getline.call
  77             while (ln = getline.call) != nil
  78                   break if ln =~ /Total Lewis/
  79                   if ln =~ /^\s*$/
  80                         #  Skip a blank line
  81                         ln = getline.call
  82                         if ln =~ /^\s*$/
  83                           #  Double blank lines indicates the end of the table
  84                           break
  85                         end
  86                   end
  87                   if ln =~ /^\s*(\d+)\. *([A-Za-z]+\*?) *\( *(\d+)\) *([A-Za-z]+) *(\d+)(- *([A-Za-z]+) *(\d+))? *(\d\.\d+)/
  88                         idx = $1.to_i
  89                         orb_kind = $2
  90                         orb_idx = $3
  91                         atom_label = $4 + ($5.to_i - 1).to_s
  92                         if $6 != nil
  93                           atom_label += "-" + $7 + ($8.to_i - 1).to_s
  94                         end
  95                         occ = $9.to_f
  96                         nbo["nbo"].push([orb_kind + "_" + orb_idx, atom_label, occ])
  97                   end
  98                 end
  99           elsif ln =~ /([A-Z]+)s in the ([A-Z]+) basis:/ || ln =~ /([A-Z]+) Fock matrix:/
 100             #  Read matrix
 101             dst = $1
 102                 src = $2
 103                 if src == nil
 104                   src = "F"   #  Fock
 105                 end
 106                 key = src + "/" + dst
 107                 getline.call
 108                 getline.call
 109                 getline.call
 110                 elements = []
 111                 labels = []
 112                 idx = 0
 113                 block = 0
 114                 while (ln = getline.call) != nil
 115                   if ln =~ /^\s*$/
 116                     #  Blank line: end of one block
 117                     ln = getline.call
 118                         if ln =~ /^\s*$/
 119                           #  Double blank lines indicates the end of the table
 120                           break
 121                         else
 122                           #  Begin next section
 123                           ln = getline.call
 124                           idx = 0
 125                           block += 1
 126                           next
 127                         end
 128                   end
 129                   ln.chomp!
 130                   ln =~ /-?\d+\.\d+/
 131                   lab = $~.pre_match.strip
 132                   a = ([$~[0]] + $~.post_match.split).map { |e| e.to_f }
 133                   if block == 0
 134                     lab =~ /^[0-9]+\. *(.*)$/
 135                     labels.push($1)
 136                   end
 137                   (elements[idx] ||= []).concat(a)
 138                   idx += 1
 139                 end
 140                 nbo[key] = LAMatrix.new(elements).transpose!
 141                 key2 = (src == "F" ? dst : src) + "_L"
 142                 if nbo[key2] == nil
 143                   nbo[key2] = labels
 144                 end
 145           end
 146         end
 147         @nbo = nbo
 148         if @nbo["AO/NAO"]
 149           #  Convert NAO based matrix into AO based matrix
 150           @nbo.keys.each { |key|
 151             if key[0..3] == "NAO/"
 152                   key2 = "AO/" + key[4..-1]
 153                   @nbo[key2] = @nbo["AO/NAO"] * @nbo[key]
 154                 end
 155           }
 156         end
 157         # puts @nbo.inspect
 158   end
 159
 160   def Molecule.read_gamess_basis_sets(fname)
 161     # $gamess_basis = Hash.new unless $gamess_basis
 162         $gamess_ecp = Hash.new unless $gamess_ecp
 163         # $gamess_basis_desc = Hash.new unless $gamess_basis_desc
 164         # $gamess_basis_keys = [] unless $gamess_basis_keys
 165         basename = File.basename(fname, ".*")
 166         keys = []
 167         descname = nil
 168     File.open(fname, "r") { |fp|
 169           while (s = fp.gets)
 170             if s =~ /^\s*!\s*/ && descname == nil
 171               #  Get the descriptive name from the first comment line
 172                   s = Regexp.last_match.post_match
 173                   if s =~ /  EMSL/
 174                     descname = Regexp.last_match.pre_match
 175                   else
 176                     descname = (s.split)[0]
 177                   end
 178                   $gamess_basis_desc[basename] = descname
 179                   next
 180                 end
 181                 ss, bas = (s.split)[0..1]     #  Tokens delimited by whitespaces
 182                 next if ss == nil || ss == ""
 183                 if ss == "$ECP"
 184                 #  Read the ECP section
 185                   ecp = $gamess_ecp[basename]
 186                   if !ecp
 187                     ecp = []
 188                         $gamess_ecp[basename] = ecp
 189                   end
 190                   keys << basename unless keys.include?(basename)
 191                   while (s = fp.gets)
 192                     break if s=~ /\$END/
 193                         ary = s.split  #  (PNAME, PTYPE, IZCORE, LMAX+1)
 194                         elem = ary[0].match(/\w{1,2}/).to_a[0]  #  ary[0] should begin with an element name
 195                         ap = Parameter.builtin.elements.find { |p| p.name.casecmp(elem) == 0 }
 196                         raise MolbyError, "the ECP definition does not begin with an element name: #{s}" if !ap
 197                         ecpdef = s
 198                         ln = 1
 199                         (0..Integer(ary[3])).each {
 200                           s = fp.gets
 201                           raise MolbyError, "the ECP definition ends unexpectedly at line #{ln} for: #{s}" if !s
 202                           ln += 1
 203                           ary = s.split  #  (NGPOT, comment)
 204                           ecpdef += s
 205                           (1..Integer(ary[0])).each {
 206                             s = fp.gets
 207                                 raise MolbyError, "the ECP definition ends unexpectedly at line #{ln} for: #{s}" if !s
 208                                 ln += 1
 209                                 ecpdef += s
 210                           }
 211                         }
 212                     ecp[ap.index] = ecpdef
 213                   end
 214                 elsif ss =~ /\W/
 215                   #  Comments or other unrecognizable lines
 216                   next
 217                 elsif (ap = Parameter.builtin.elements.find { |p| p.name.casecmp(ss) == 0 || p.fullname.casecmp(ss) == 0 })
 218                   #  Valid basis definition
 219                   if bas == nil || bas =~ /\W/
 220                     bas = basename
 221                   end
 222                   basis = $gamess_basis[bas]
 223                   if !basis
 224                     basis = []
 225                         $gamess_basis[bas] = basis
 226                   end
 227                   keys << bas unless keys.include?(bas)
 228                   basdef = ""
 229                   while (s = fp.gets) && s =~ /\S/
 230                     basdef += s
 231                   end
 232                   basis[ap.index] = basdef
 233                 else
 234                   raise MolbyError, "ss is not a valid element symbol or name: #{s}"
 235                 end
 236           end
 237     }
 238         unless $gamess_basis_keys.include?(basename)
 239           $gamess_basis_keys.push(basename)
 240         end
 241   end
 242
 243   #  Execute GAMESS (inferior copy of rungms script)
 244   #  inpname is the input file
 245   #  mol (optional) is the molecule from which the GAMESS input was built.
 246   #  If mol is specified and RUNTYP=OPTIMIZE, then the intermediate structures are
 247   #  displayed real-time.
 248   def Molecule.execute_gamess(inpname, mol = nil)
 249     gmsname = get_global_settings("gamess.executable_path")
 250     gmsdir = nil
 251     gmsvers = nil
 252
 253         cygwin_version = false  #  For windows: old cygwin version
 254
 255     while 1
 256       if gmsname == nil || !File.exist?(gmsname)
 257         gmsname = Dialog.open_panel("Please locate the GAMESS executable")
 258         exit if gmsname == nil
 259       end
 260       gmsbase = File.basename(gmsname)
 261       gmsdir = File.dirname(gmsname)
 262       if gmsbase =~ /gamess\.(.*)\.(exe|x)$/i
 263         gmsvers = $1
 264         break
 265       else
 266         gmsname = nil
 267         error_message_box(gmsbase + " does not look like a GAMESS executable!")
 268       end
 269     end
 270
 271 #       if mol == nil
 272 #               mol = Molecule.open(inpname)
 273 #               if mol == nil
 274 #                       error_message_box("Cannot open #{inpname} as GAMESS input")
 275 #                       return
 276 #               end
 277 #       end
 278
 279     inpbase = File.basename(inpname)
 280     inpdir = File.dirname(inpname)
 281     inpbody = inpbase.sub(/\.inp$/, "")
 282     logbase = inpbody + ".log"
 283
 284     set_global_settings("gamess.executable_path", gmsname)
 285
 286     ncpus = get_global_settings("gamess.ncpus").to_i
 287         if ncpus == 0
 288           ncpus = 1
 289         end
 290
 291     #  Prepare the scratch directory
 292         begin
 293           scrdir = create_temp_dir("gamess", inpbody)
 294         rescue
 295           error_message_box($!.to_s)
 296           return
 297         end
 298 #    scrdir = $home_directory + "/Molby/gamess"
 299 #       begin
 300 #         mkdir_recursive(scrdir)
 301 #       rescue
 302 #         error_message_box("Cannot create directory #{scrdir}: " + $!.to_s)
 303 #         return
 304 #       end
 305
 306 #    scrdir = scrdir + "/" + inpbody + "." + $$.to_s + ".0"
 307 #    n = 0
 308 #    while File.exist?(scrdir)
 309 #      scrdir = scrdir.sub(".#{n}", ".#{n + 1}")
 310 #      n += 1
 311 #    end
 312 #    Dir.mkdir(scrdir)
 313
 314     if $platform == "win"
 315       sep = "\\"
 316       scrdir.gsub!("/", sep)
 317       gmsdir.gsub!("/", sep)
 318     else
 319       sep = "/"
 320     end
 321
 322         #  Old (obsolete) cygwin version, using ddikick
 323     if $platform == "win"
 324           if gmsvers == "11"
 325             cygwin_version = true
 326           end
 327         end
 328
 329     #  Get the host name etc.
 330     hostname = backquote("hostname").chomp
 331     if $platform == "win"
 332           s = backquote("cmd.exe /c dir \"#{scrdir}\"")
 333       freebytes = s.split("\n").pop.match(/([0-9,]+)[^0-9]*$/).to_a[1]
 334       if freebytes
 335         freebytes = (freebytes.gsub(",","").to_i / 1024).to_s + " Kbytes"
 336       else
 337         freebytes = "(unknown)"
 338       end
 339       uname = backquote("cmd.exe /c ver").to_s.gsub("\n", "")
 340     else
 341       freebytes = `df -k "#{scrdir}"`
 342       uname = `uname`.chomp
 343     end
 344
 345     #  Redirect standard output to the log file
 346     logname = scrdir + sep + logbase
 347     fpout = File.open(logname, "w")
 348         if cygwin_version
 349           #  The cygwin version uses LF as the eol character
 350           fpout.binmode
 351         end
 352     fpout.print "----- GAMESS execution script -----\n"
 353     fpout.print "This job is running on host #{hostname}\n"
 354     fpout.print "under operating system #{uname} at #{Time.now.to_s}\n"
 355     fpout.print "Available scratch disk space (Kbyte units) at beginning of the job is\n"
 356     fpout.print freebytes + "\n"
 357
 358     #  Copy the input file
 359         scrprefix = scrdir + sep + inpbody
 360         if $platform == "win"
 361           scrbody = inpbody
 362         else
 363       scrbody = scrprefix
 364         end
 365     filecopy(inpname, scrprefix + ".F05")
 366     File.open("#{scrdir}/.in_use", "w") { |fp| }
 367
 368     #  Prepare environmental variables
 369     auxdir = "#{gmsdir}#{sep}auxdata"
 370     ENV["ERICFMT"] = "#{auxdir}#{sep}ericfmt.dat"
 371     ENV["MCPPATH"] = "#{auxdir}#{sep}MCP"
 372     ENV["BASPATH"] = "#{auxdir}#{sep}BASES"
 373     ENV["QUANPOL"] = "#{auxdir}#{sep}QUANPOL"
 374     ENV["EXTBAS"] = "/dev/null"
 375     ENV["IRCDATA"] = "#{scrbody}.irc"
 376         ENV["RESTART"] = "#{scrbody}.rst"
 377         ENV["TRAJECT"] = "#{scrbody}.trj"
 378     ENV["PUNCH"] = "#{scrbody}.dat"
 379     ENV["INPUT"] = "#{scrbody}.F05"
 380     ENV["AOINTS"] = "#{scrbody}.F08"
 381     ENV["MOINTS"] = "#{scrbody}.F09"
 382     ENV["DICTNRY"] = "#{scrbody}.F10"
 383     ENV["DRTFILE"] = "#{scrbody}.F11"
 384     ENV["CIVECTR"] = "#{scrbody}.F12"
 385     ENV["CASINTS"] = "#{scrbody}.F13"
 386     ENV["CIINTS"] = "#{scrbody}.F14"
 387     ENV["WORK15"] = "#{scrbody}.F15"
 388     ENV["WORK16"] = "#{scrbody}.F16"
 389     ENV["CSFSAVE"] = "#{scrbody}.F17"
 390     ENV["FOCKDER"] = "#{scrbody}.F18"
 391     ENV["WORK19"] = "#{scrbody}.F19"
 392     ENV["DASORT"] = "#{scrbody}.F20"
 393     ENV["DFTINTS"] = "#{scrbody}.F21"
 394     ENV["DFTGRID"] = "#{scrbody}.F22"
 395     ENV["JKFILE"] = "#{scrbody}.F23"
 396     ENV["ORDINT"] = "#{scrbody}.F24"
 397     ENV["EFPIND"] = "#{scrbody}.F25"
 398     ENV["PCMDATA"] = "#{scrbody}.F26"
 399     ENV["PCMINTS"] = "#{scrbody}.F27"
 400     ENV["MLTPL"] = "#{scrbody}.F28"
 401     ENV["MLTPLT"] = "#{scrbody}.F29"
 402     ENV["DAFL30"] = "#{scrbody}.F30"
 403     ENV["SOINTX"] = "#{scrbody}.F31"
 404     ENV["SOINTY"] = "#{scrbody}.F32"
 405     ENV["SOINTZ"] = "#{scrbody}.F33"
 406     ENV["SORESC"] = "#{scrbody}.F34"
 407     ENV["SIMEN"] = "#{scrbody}.simen"
 408     ENV["SIMCOR"] = "#{scrbody}.simcor"
 409     ENV["GCILIST"] = "#{scrbody}.F37"
 410     ENV["HESSIAN"] = "#{scrbody}.F38"
 411     ENV["SOCCDAT"] = "#{scrbody}.F40"
 412     ENV["AABB41"] = "#{scrbody}.F41"
 413     ENV["BBAA42"] = "#{scrbody}.F42"
 414     ENV["BBBB43"] = "#{scrbody}.F43"
 415     ENV["MCQD50"] = "#{scrbody}.F50"
 416     ENV["MCQD51"] = "#{scrbody}.F51"
 417     ENV["MCQD52"] = "#{scrbody}.F52"
 418     ENV["MCQD53"] = "#{scrbody}.F53"
 419     ENV["MCQD54"] = "#{scrbody}.F54"
 420     ENV["MCQD55"] = "#{scrbody}.F55"
 421     ENV["MCQD56"] = "#{scrbody}.F56"
 422     ENV["MCQD57"] = "#{scrbody}.F57"
 423     ENV["MCQD58"] = "#{scrbody}.F58"
 424     ENV["MCQD59"] = "#{scrbody}.F59"
 425     ENV["MCQD60"] = "#{scrbody}.F60"
 426     ENV["MCQD61"] = "#{scrbody}.F61"
 427     ENV["MCQD62"] = "#{scrbody}.F62"
 428     ENV["MCQD63"] = "#{scrbody}.F63"
 429     ENV["MCQD64"] = "#{scrbody}.F64"
 430     ENV["NMRINT1"] = "#{scrbody}.F61"
 431     ENV["NMRINT2"] = "#{scrbody}.F62"
 432     ENV["NMRINT3"] = "#{scrbody}.F63"
 433     ENV["NMRINT4"] = "#{scrbody}.F64"
 434     ENV["NMRINT5"] = "#{scrbody}.F65"
 435     ENV["NMRINT6"] = "#{scrbody}.F66"
 436     ENV["DCPHFH2"] = "#{scrbody}.F67"
 437     ENV["DCPHF21"] = "#{scrbody}.F68"
 438     ENV["GVVPT"] = "#{scrbody}.F69"
 439
 440     #    next files are used only during coupled cluster runs, so let's
 441     #    display the numerous definitions only if they are to be used.
 442         #    Don't set these variables on windows, where the memory for environmental variables
 443         #    is limited.
 444         if $platform != "win"
 445     ENV["CCREST"] = "#{scrbody}.F70"
 446     ENV["CCDIIS"] = "#{scrbody}.F71"
 447     ENV["CCINTS"] = "#{scrbody}.F72"
 448     ENV["CCT1AMP"] = "#{scrbody}.F73"
 449     ENV["CCT2AMP"] = "#{scrbody}.F74"
 450     ENV["CCT3AMP"] = "#{scrbody}.F75"
 451     ENV["CCVM"] = "#{scrbody}.F76"
 452     ENV["CCVE"] = "#{scrbody}.F77"
 453     ENV["EOMSTAR"] = "#{scrbody}.F80"
 454     ENV["EOMVEC1"] = "#{scrbody}.F81"
 455     ENV["EOMVEC2"] = "#{scrbody}.F82"
 456     ENV["EOMHC1"] = "#{scrbody}.F83"
 457     ENV["EOMHC2"] = "#{scrbody}.F84"
 458     ENV["EOMHHHH"] = "#{scrbody}.F85"
 459     ENV["EOMPPPP"] = "#{scrbody}.F86"
 460     ENV["EOMRAMP"] = "#{scrbody}.F87"
 461     ENV["EOMRTMP"] = "#{scrbody}.F88"
 462     ENV["EOMDG12"] = "#{scrbody}.F89"
 463     ENV["MMPP"] = "#{scrbody}.F90"
 464     ENV["MMHPP"] = "#{scrbody}.F91"
 465     ENV["MMCIVEC"] = "#{scrbody}.F92"
 466     ENV["MMCIVC1"] = "#{scrbody}.F93"
 467     ENV["MMCIITR"] = "#{scrbody}.F94"
 468     ENV["MMNEXM"] = "#{scrbody}.F95"
 469     ENV["MMNEXE"] = "#{scrbody}.F96"
 470     ENV["MMNREXM"] = "#{scrbody}.F97"
 471     ENV["MMNREXE"] = "#{scrbody}.F98"
 472         end
 473     #
 474     #     next are for TDHFX code, not used by current GAMESS
 475     #
 476         if false
 477     ENV["OLI201"] = "#{scrbody}.F201"
 478     ENV["OLI202"] = "#{scrbody}.F202"
 479     ENV["OLI203"] = "#{scrbody}.F203"
 480     ENV["OLI204"] = "#{scrbody}.F204"
 481     ENV["OLI205"] = "#{scrbody}.F205"
 482     ENV["OLI206"] = "#{scrbody}.F206"
 483     ENV["OLI207"] = "#{scrbody}.F207"
 484     ENV["OLI208"] = "#{scrbody}.F208"
 485     ENV["OLI209"] = "#{scrbody}.F209"
 486     ENV["OLI210"] = "#{scrbody}.F210"
 487     ENV["OLI211"] = "#{scrbody}.F211"
 488     ENV["OLI212"] = "#{scrbody}.F212"
 489     ENV["OLI213"] = "#{scrbody}.F213"
 490     ENV["OLI214"] = "#{scrbody}.F214"
 491     ENV["OLI215"] = "#{scrbody}.F215"
 492     ENV["OLI216"] = "#{scrbody}.F216"
 493     ENV["OLI217"] = "#{scrbody}.F217"
 494     ENV["OLI218"] = "#{scrbody}.F218"
 495     ENV["OLI219"] = "#{scrbody}.F219"
 496     ENV["OLI220"] = "#{scrbody}.F220"
 497     ENV["OLI221"] = "#{scrbody}.F221"
 498     ENV["OLI222"] = "#{scrbody}.F222"
 499     ENV["OLI223"] = "#{scrbody}.F223"
 500     ENV["OLI224"] = "#{scrbody}.F224"
 501     ENV["OLI225"] = "#{scrbody}.F225"
 502     ENV["OLI226"] = "#{scrbody}.F226"
 503     ENV["OLI227"] = "#{scrbody}.F227"
 504     ENV["OLI228"] = "#{scrbody}.F228"
 505     ENV["OLI229"] = "#{scrbody}.F229"
 506     ENV["OLI230"] = "#{scrbody}.F230"
 507     ENV["OLI231"] = "#{scrbody}.F231"
 508     ENV["OLI232"] = "#{scrbody}.F232"
 509     ENV["OLI233"] = "#{scrbody}.F233"
 510     ENV["OLI234"] = "#{scrbody}.F234"
 511     ENV["OLI235"] = "#{scrbody}.F235"
 512     ENV["OLI236"] = "#{scrbody}.F236"
 513     ENV["OLI237"] = "#{scrbody}.F237"
 514     ENV["OLI238"] = "#{scrbody}.F238"
 515     ENV["OLI239"] = "#{scrbody}.F239"
 516     end
 517
 518     if $platform == "win"
 519         #  if ncpus < 2
 520         #    ncpus = 2
 521         #  end
 522           if cygwin_version
 523             #  Old (obsolete) cygwin version, using ddikick
 524         fpout.print "ddikick will run #{ncpus} compute process\n"
 525                 ENV["CYGWIN"] = "nodosfilewarning"
 526           else
 527         fpout.print "Microsoft MPI will be running GAMESS on 1 node.\n"
 528         fpout.print "The binary kicked off by 'mpiexec' is gamess.#{gmsvers}.exe\n"
 529         fpout.print "MS-MPI will run #{ncpus*2} compute process\n"
 530         #  File containing environmental variables
 531         envfil = "#{scrprefix}.GMS.ENV"
 532         fp = File.open(envfil, "w")
 533         ENV.each { |k, v| fp.print "#{k}=#{v}\n" }
 534         fp.close
 535         #  File containing arguments to mpiexec
 536         procfil = "#{scrprefix}.processes.mpd"
 537         fp = File.open(procfil, "w")
 538         fp.print "-env ENVFIL \"#{envfil}\" -wdir \"#{scrdir}\" -n #{ncpus*2} \"#{gmsdir}#{sep}gamess.#{gmsvers}.exe\"\n"
 539         fp.close
 540           end
 541     end
 542
 543     fpout.close
 544     fplog = File.open(logname, "r")
 545     size = 0
 546     lines = []
 547         lineno = 0
 548     last_line = ""
 549     search_mode = 0
 550         nserch = -1
 551         rflag = 0
 552         ne_alpha = ne_beta = 0
 553         mo_count = 0
 554         ncomps = 0
 555
 556     #  Callback procs
 557         term_callback = lambda { |m, n|
 558           msg = "GAMESS execution on #{inpbase} "
 559           hmsg = "GAMESS "
 560           if n == 0
 561             msg += "succeeded."
 562                 hmsg += "Completed"
 563                 icon = :info
 564           else
 565             msg += "failed with status #{n}."
 566                 hmsg += "Failed"
 567                 icon = :error
 568           end
 569           msg += "\n(In directory #{inpdir})"
 570
 571           File.delete("#{scrdir}/.in_use")
 572
 573           ext_to_keep = [".dat", ".rst", ".trj", ".efp", ".gamma", ".log"]
 574           ext_to_keep.each { |ex|
 575             if File.exists?("#{scrprefix}#{ex}")
 576                   filecopy("#{scrprefix}#{ex}", "#{inpdir}#{sep}#{inpbody}#{ex}")
 577             end
 578           }
 579           Dir.foreach(scrdir) { |file|
 580             if file != "." && file != ".." && !ext_to_keep.include?(File.extname(file))
 581                   File.delete("#{scrdir}#{sep}#{file}")
 582             end
 583           }
 584
 585           begin
 586             erase_old_logs(scrdir, "latest", 5)
 587           rescue
 588             error_message_box($!.to_s)
 589                 return
 590           end
 591
 592           if (script = get_global_settings("gamess.postfix_script")) != nil && script != ""
 593             eval(script)
 594           end
 595
 596           if mol != nil
 597             message_box(msg, hmsg, :ok, icon)
 598       end
 599     }
 600
 601     timer_callback = lambda { |m, n|
 602       fplog.seek(0, IO::SEEK_END)
 603       sizec = fplog.tell
 604       if sizec > size
 605         #  Read new lines
 606         fplog.seek(size, IO::SEEK_SET)
 607         fplog.each_line { |line|
 608           if line[-1, 1] == "\n"
 609             lines.push(last_line + line)
 610             last_line = ""
 611           else
 612             last_line += line
 613             break
 614           end
 615         }
 616         size = fplog.tell
 617         last_i = nil
 618         i = 0
 619         while i < lines.count
 620           line = lines[i]
 621           if line =~ /GEOMETRY SEARCH POINT NSERCH= *(\d+)/
 622             nserch = $1.to_i
 623             last_i = i
 624                         search_mode = 1
 625                         energy = nil
 626                   elsif line =~ /START OF DRC CALCULATION/
 627                     search_mode = 3
 628                         nserch = 1
 629                         energy = nil
 630                   elsif line =~ /CONSTRAINED OPTIMIZATION POINT/
 631                     search_mode = 2
 632                         energy = nil
 633                   elsif line =~ /POINT *([0-9]+) *ON THE REACTION PATH/
 634                     nserch = $1.to_i
 635                         last_i = i
 636                   elsif line =~ /ATOMIC BASIS SET/
 637                         j = i
 638                         while j < lines.count
 639                           line = lines[j]
 640                           break if line =~ /TOTAL NUMBER OF BASIS SET/
 641                           j += 1
 642                         end
 643                         if j < lines.count
 644                           #  Found
 645                           bs_lines = []
 646                           ii = i
 647                           while ii <= j
 648                             bs_lines.push(lines[ii].chomp)
 649                                 ii += 1
 650                           end
 651                           begin
 652                             if mol
 653                                   ncomps = mol.sub_load_gamess_log_basis_set(bs_lines, lineno + i)
 654                                 end
 655                           rescue
 656           $stderr.write($!.to_s + "\n")
 657           $stderr.write($!.backtrace.inspect + "\n")
 658                           end
 659                           last_i = j
 660                         else
 661                           break  #  Wait until all basis set lines are read
 662                         end
 663                   elsif line =~ /NUMBER OF OCCUPIED ORBITALS/
 664                         line =~ /=\s*(\d+)/
 665                         n = Integer($1)
 666                         if line =~ /ALPHA/
 667                           ne_alpha = n
 668                         else
 669                           ne_beta = n
 670                         end
 671                   elsif line =~ /SCFTYP=(\w+)/
 672                         scftyp = $1
 673                         if ne_alpha > 0 || ne_beta > 0
 674                           rflag = 0
 675                           case scftyp
 676                           when "RHF"
 677                                 rflag = 1
 678                           when "ROHF"
 679                                 rflag = 2
 680                           end
 681                         end
 682                   elsif line =~ /^\s*(EIGENVECTORS|MOLECULAR ORBITALS)\s*$/
 683                         if mo_count == 0 && mol
 684                           mol.clear_mo_coefficients
 685                           mol.set_mo_info(:type=>["UHF", "RHF", "ROHF"][rflag], :alpha=>ne_alpha, :beta=>ne_beta)
 686                         end
 687                         i += 2
 688                         j = i
 689                         mo_count += 1
 690                         while j < lines.count
 691                           line = lines[j]
 692                           break if line =~ /\.\.\.\.\.\./ || line =~ /----------------/
 693                           j += 1
 694                         end
 695                         if j == lines.count
 696                           break  #  Wait until complete MO info are read
 697                         end
 698                         ii = i
 699                         mo_lines = []
 700                         while ii < j
 701                           mo_lines.push(lines[ii].chomp)
 702                           ii += 1
 703                         end
 704                         begin
 705                           if mol
 706                             mol.sub_load_gamess_log_mo_coefficients(mo_lines, lineno + i, ncomps)
 707                           end
 708                         rescue
 709                           $stderr.write($!.to_s + "\n")
 710                           $stderr.write($!.backtrace.inspect + "\n")
 711                         end
 712                         last_i = j
 713           elsif line =~ /NSERCH: *([0-9]+)/
 714           #  print line
 715                         if mol
 716                           n = $1
 717                           line =~ /E= +([-.0-9]+)/
 718                           energy = $1.to_f
 719                           line =~ /GRAD\. MAX= +([-.0-9]+)/
 720                           grad = $1
 721                           mol.show_text("Search: #{n}\nGradient: #{grad}")
 722                           # mol.set_property("energy", energy)
 723                         end
 724                         last_i = i
 725                   elsif line =~ /TOTAL ENERGY += *([-.0-9]+)/
 726                     energy = $1
 727                         if mol && search_mode == 2
 728                           mol.show_text("Point: #{nserch}\nEnergy: #{energy}")
 729                         end
 730                         energy = energy.to_f
 731                   elsif line =~ /FINAL .* ENERGY IS +([-.0-9]+) *AFTER/
 732                     energy = $1.to_f
 733                         if mol && search_mode == 0
 734                           mol.set_property("energy", energy)
 735                         end
 736           elsif nserch > 0 && line =~ /ddikick\.x/
 737             last_i = -1
 738             break
 739           elsif mol && nserch > 0 && (line =~ /COORDINATES OF ALL ATOMS/ || line =~ /COORDINATES \(IN ANGSTROM\) FOR \$DATA GROUP ARE/)
 740             #  There should be (natoms) lines
 741                         if line =~ /COORDINATES OF ALL ATOMS/
 742                           #  Skip header lines
 743                           i += 2
 744                     end
 745             if i + mol.natoms + 1 <= lines.count
 746               coords = []
 747               (i + 1...i + 1 + mol.natoms).each { |j|
 748                 name, charge, x, y, z = lines[j].split
 749                 coords.push(Vector3D[x.to_f, y.to_f, z.to_f])
 750               }
 751               mol.create_frame([coords])
 752                           if search_mode == 1 && energy
 753                             mol.set_property("energy", energy)
 754                           end
 755               mol.display
 756               last_i = i + mol.natoms
 757               i = last_i   #  Skip the processed lines
 758             end
 759                   elsif line =~ /N A T U R A L   B O N D   O R B I T A L   A N A L Y S I S/
 760                     #  NBO output
 761                     j = i + 1
 762                         while j < lines.count
 763                           break if lines[j] =~ /done with NBO analysis/
 764                           j += 1
 765                         end
 766                         if j < lines.count
 767                           nbo_lines = lines[i..j]
 768                           mol.import_nbo_log(nbo_lines) rescue puts "Error: #{$!}, #{$!.backtrace.inspect}"
 769                           last_i = j
 770                           i = last_i  #  Skip the processed lines
 771                         else
 772                           break  #  Wait until NBO is done
 773                         end
 774           end
 775           i += 1
 776         end
 777         if last_i == -1
 778           lines.clear
 779           break
 780         elsif last_i
 781           lines[0..last_i] = nil
 782                   lineno += last_i + 1
 783         end
 784       end
 785       true
 786     }
 787
 788     if (script = get_global_settings("gamess.prefix_script")) != nil && script != ""
 789           eval(script)
 790         end
 791
 792     if $platform == "win"
 793           if gmsvers == "11"
 794             hosts = "localhost " * ncpus
 795             cmdline = "cmd.exe /c \"#{gmsdir}/ddikick.exe #{gmsdir}/gamess.#{gmsvers}.exe #{inpbody} -ddi #{ncpus} #{ncpus} #{hosts} -scr #{scrdir} < NUL >>#{logname}"
 796           else
 797             cmdline = "cmd.exe /c \"mpiexec -configfile #{procfil} >>#{logname}"
 798           end
 799         else
 800           hosts = "localhost " * ncpus
 801           cmdline = "/bin/sh -c '#{gmsdir}/ddikick.x #{gmsdir}/gamess.#{gmsvers}.x #{inpbody} -ddi #{ncpus} #{ncpus} #{hosts} -scr #{scrdir} < /dev/null >>#{logname}'"
 802         end
 803
 804         if mol
 805           pid = mol.call_subprocess_async(cmdline, term_callback, timer_callback)
 806           if pid < 0
 807             error_message_box("GAMESS failed to start. Please examine GAMESS installation.")
 808             return
 809           end
 810         else
 811           pid = call_subprocess(cmdline, "Running GAMESS")
 812           term_callback.call(nil, pid)
 813           return pid
 814         end
 815
 816   end
 817
 818   def export_gamess(fname, hash)
 819
 820     def reorder_array(ary, ordered_sub_ary)
 821           return (ordered_sub_ary & ary) + (ary - ordered_sub_ary)
 822         end
 823
 824         now = Time.now.to_s
 825         if fname
 826           basename = File.basename(fname, ".*")
 827         else
 828           basename = File.basename(self.name, ".*")
 829         end
 830
 831         #  Various settings
 832         icharg = hash["charge"]
 833         mult = hash["mult"]
 834         runtyp = ["ENERGY", "PROP", "OPTIMIZE", "SADPOINT", "IRC"][hash["runtype"].to_i]
 835         scftyp = ["RHF", "ROHF", "UHF"][hash["scftype"].to_i]
 836         bssname = hash["basis"]
 837         bssname2 = hash["secondary_basis"]
 838         if hash["use_secondary_basis"].to_i != 0 && bssname2 != bssname
 839           use_2nd = true
 840           element2 = hash["secondary_elements"].split(/[\s,]+/).map { |name| name.capitalize }
 841         else
 842           use_2nd = false
 843           bssname2 = bssname
 844         end
 845         basis = $gamess_basis[bssname]
 846         basis2 = $gamess_basis[bssname2]
 847
 848         #  Use effective core potentials?
 849         ecp = $gamess_ecp[bssname]
 850         ecp2 = $gamess_ecp[bssname2]
 851         ecp_read = (ecp || ecp2 ? "ECP=READ " : "")
 852
 853         #  Use only one built-in basis set?
 854         gbasis = nil
 855         if !use_2nd
 856           case bssname
 857           when "PM3"
 858                 gbasis = "PM3 NGAUSS=3"
 859           when "STO3G"
 860                 gbasis = "STO NGAUSS=3"
 861           when "321G"
 862                 gbasis = "N21 NGAUSS=3"
 863           when "631G"
 864                 gbasis = "N31 NGAUSS=6"
 865           when "631Gd"
 866                 gbasis = "N31 NGAUSS=6 NDFUNC=1"
 867           when "631Gdp"
 868                 gbasis = "N31 NGAUSS=6 NDFUNC=1 NPFUNC=1"
 869           when "6311G"
 870                 gbasis = "N311 NGAUSS=6"
 871           when "6311Gdp"
 872                 gbasis = "N311 NGAUSS=6 NDFUNC=1 NPFUNC=1"
 873           end
 874         end
 875
 876         #  Count non-dummy atoms
 877         natoms = 0
 878         each_atom { |ap|
 879           natoms += 1 if ap.atomic_number != 0
 880         }
 881
 882         #  Fill hash with default values
 883         h = (hash["CONTRL"] ||= Hash.new)
 884         h["COORD"] ||= "UNIQUE"
 885         h["EXETYP"] ||= "RUN"
 886         h["ICHARG"] ||= (icharg || 0).to_s
 887         h["ICUT"] ||= "20"
 888         h["INTTYP"] ||= "HONDO"
 889         h["ITOL"] ||= "30"
 890         h["MAXIT"] ||= "200"
 891         h["MOLPLT"] ||= ".T."
 892         h["MPLEVL"] ||= "0"
 893         h["MULT"] ||= mult.to_s
 894         h["QMTTOL"] ||= "1e-08"
 895         h["RUNTYP"] ||= runtyp
 896         if (hash["dft"] || 0) != 0 && hash["dfttype"]
 897           h["DFTTYP"] ||= hash["dfttype"]
 898         end
 899         if (hash["use_internal"] || 0) != 0 && (hash["runtype"] >= 2 || h["RUNTYP"] == "OPTIMIZE" || h["RUNTYP"] == "SADPOINT" || h["RUNTYP"] == "IRC")
 900           nzvar = natoms * 3 - 6  #  TODO: 3N-5 for linear molecules
 901           h["NZVAR"] ||= nzvar.to_s
 902         else
 903           nzvar = 0
 904         end
 905         h["SCFTYP"] ||= scftyp
 906         if ecp_read != ""
 907           h["ECP"] ||= "READ"
 908         end
 909         h["UNITS"] ||= "ANGS"
 910
 911         h = (hash["SCF"] ||= Hash.new)
 912         h["CONV"] ||= "1.0E-06"
 913         h["DIRSCF"] ||= ".T."
 914         h["FDIFF"] ||= ".T."
 915         h["DAMP"] ||= ".T."
 916
 917         h = (hash["STATPT"] ||= Hash.new)
 918         h["NSTEP"] ||= "400"
 919     if runtyp == "SADPOINT"
 920         h["OPTTOL"] ||= "1.0E-05"
 921         h["HESS"] ||= "CALC"
 922         h["HSSEND"] ||= ".T."
 923     elsif runtyp == "IRC"
 924         h["OPTTOL"] ||= "1.0E-05"
 925         h["HESS"] ||= "READ"
 926         h["HSSEND"] ||= ".T."
 927     else
 928         h["OPTTOL"] ||= "1.0E-06"
 929     end
 930     if hash["eliminate_freedom"] == 0
 931       h["PROJCT"] ||= ".F."
 932     end
 933
 934         h = (hash["SYSTEM"] ||= Hash.new)
 935         h["MEMDDI"] ||= "0"
 936     if runtyp == "SADPOINT" || runtyp == "IRC"
 937         h["MWORDS"] ||= "32"
 938     else
 939         h["MWORDS"] ||= "16"
 940     end
 941         h["TIMLIM"] ||= "50000"
 942
 943         h = (hash["GUESS"] ||= Hash.new)
 944         h["GUESS"] ||= "HUCKEL"
 945
 946         if gbasis
 947           h = (hash["BASIS"] ||= Hash.new)
 948           h["GBASIS"] ||= gbasis
 949         end
 950
 951         if nzvar > 0
 952           h = (hash["ZMAT"] ||= Hash.new)
 953           h["DLC"] ||= ".T."
 954           h["AUTO"] ||= ".T."
 955         end
 956
 957     if runtyp == "IRC"
 958         h = (hash["IRC"] ||= Hash.new)
 959         h["SADDLE"] = ".T."
 960         h["TSENGY"] = ".T."
 961         h["FORWRD"] = ".T."
 962         h["NPOINT"] = "100"
 963         h["STRIDE"] = "0.2"
 964         h["OPTTOL"] = "1.0E-05"
 965     end
 966
 967         if (hash["esp"] || 0) != 0
 968           h = (hash["ELPOT"] ||= Hash.new)
 969           h["IEPOT"] ||= "1"
 970           h["OUTPUT"] ||= "PUNCH"
 971           h["WHERE"] ||= "PDC"
 972           h = (hash["PDC"] ||= Hash.new)
 973           h["CONSTR"] ||= "NONE"
 974           h["PTSEL"] ||= "CONNOLLY"
 975         end
 976
 977     if (hash["include_nbo"] || 0) != 0
 978       h = (hash["NBO"] ||= Hash.new)
 979       s = ""
 980       ["nao", "nbo", "nho", "nlmo", "pnao", "pnbo", "pnho", "pnlmo"].each { |nao|
 981         if (hash[nao] || 0) != 0
 982           s += " f" + nao + " ao" + nao
 983         end
 984       }
 985       h[s] = ""
 986     end
 987
 988         if fname
 989           fp = File.open(fname, "wb")
 990         else
 991           fp = MemoryIO.new
 992         end
 993     if fp
 994           fp.print "!  GAMESS input\n"
 995           fp.print "!  Generated by Molby at #{now}\n"
 996           fp.print "!  Basis set: " + ($gamess_basis_desc[bssname] || "(not specified)") + "\n"
 997           if use_2nd
 998             fp.print "!  [" + element2.join(", ") + "]: " + ($gamess_basis_desc[bssname2] || "(not specified)") + "\n"
 999           end
1000           ordered = ["CONTRL", "SCF", "STATPT", "SYSTEM", "GUESS", "BASIS", "ZMAT", "ELPOT", "PDC"]
1001           controls = reorder_array(hash.keys.select { |k| k != "key_order" && hash[k].is_a?(Hash) },
1002                 hash["key_order"] || ordered)
1003           controls.each { |k|
1004             h = hash[k]
1005                 next if h == nil || h.size == 0 || (h["key_order"] && h.size == 1)
1006                 if k == "CONTRL"
1007                   ordered = ["COORD", "EXETYP", "ICHARG", "ICUT", "INTTYP", "ITOL", "MAXIT", "MOLPLT", "MPLEVL",
1008                              "MULT", "QMTTOL", "RUNTYP", "SCFTYP", "ECP", "UNITS", "DFTTYP", "NZVAR"]
1009                 elsif k == "SCF"
1010                   ordered = ["CONV", "DIRSCF", "FDIFF", "DAMP"]
1011                 elsif k == "STATPT"
1012                   ordered = ["NSTEP", "OPTTOL", "HESS", "HSSEND"]
1013                 elsif k == "SYSTEM"
1014                   ordered = ["MEMDDI", "MWORDS", "TIMLIM"]
1015                 elsif k == "GUESS"
1016                   ordered = ["GUESS"]
1017                 elsif k == "BASIS"
1018                   ordered = ["GBASIS"]
1019                 elsif k == "ZMAT"
1020                   ordered = ["DLC", "AUTO"]
1021                 elsif k == "ELPOT"
1022                   ordered = ["IEPOT", "OUTPUT", "WHERE"]
1023                 elsif k == "PDC"
1024                   ordered = ["CONSTR", "PTSEL"]
1025         elsif k == "IRC"
1026           ordered = ["SADDLE", "TSENGY", "FORWRD", "NPOINT", "STRIDE", "OPTTOL"]
1027                 else
1028                   ordered = []
1029                 end
1030             keys = reorder_array(h.keys, h["key_order"] || ordered)
1031                 n = 0
1032                 keys.each_with_index { |kk, i|
1033                   v = h[kk]
1034                   if n == 0
1035                     fp.printf " $%-6s", k
1036                   elsif n % 3 == 0
1037                     fp.print "\n        "
1038                   end
1039                   if v == nil || v == ""
1040                     #  No value
1041                         fp.print " " + kk
1042                   else
1043                     fp.printf " %s=%s", kk, h[kk].to_s
1044                   end
1045                   n += 1
1046                 }
1047                 if n > 0
1048                   fp.print " $END\n"
1049                 end
1050           }
1051           fp.print " $DATA\n#{basename}\nC1 0\n"
1052           secondary = []
1053           each_atom { |ap|
1054             next if ap.atomic_number == 0
1055                 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
1056                 if use_2nd && element2.include?(ap.element)
1057                   secondary[ap.index] = true
1058                 end
1059             if !gbasis
1060                   #  Basis specification followed by a blank line
1061                   bas = (secondary[ap.index] ? basis2 : basis)
1062                   if bas.is_a?(Array)
1063                     bas = bas[ap.atomic_number]
1064                   end
1065                   if bas
1066                     fp.print bas
1067                     fp.print "\n"
1068                   else
1069                     puts "Warning: basis set is not defined for atom #{ap.index}, element #{ap.element}"
1070                   end
1071                 end
1072           }
1073           fp.print " $END\n"
1074           ecp_ary = []
1075           if ecp || ecp2
1076             fp.print " $ECP\n"
1077                 each_atom { |ap|
1078                   an = ap.atomic_number
1079                   next if an == 0
1080                   ecpp = (secondary[ap.index] ? ecp2 : ecp)
1081                   e = ecp_ary[an] || (ecpp && ecpp[an])
1082                   if e
1083                     #  Cache the PNAME of the $ECP entry and re-use it
1084                     ecp_ary[an] ||= (e.split)[0] + "\n"
1085                   else
1086                     e = ap.element.upcase + "-ECP NONE\n"
1087                   end
1088                   fp.print e
1089                 }
1090             fp.print " $END\n"
1091           end
1092         end
1093         if fname == nil
1094           s = fp.buffer
1095       fp.empty
1096           return s
1097         else
1098           fp.close
1099           return fname
1100         end
1101   end
1102
1103   def copy_section_from_gamess_output
1104       if @gamess_output_fname
1105           dir = File.dirname(@gamess_output_fname)
1106       else
1107           dir = self.dir
1108       end
1109       fname = Dialog.open_panel("Select GAMESS Output:", dir, "*.log;*.dat")
1110       if fname
1111           fp = open(fname, "rb")
1112           if fp
1113               mtime = fp.mtime
1114               if @gamess_output_fname == fname && @gamess_output_mtime == mtime
1115                   #  Use the cached value
1116                   k = @gamess_output_sections
1117               else
1118                   pos = 0
1119                   k = Hash.new
1120                   fp.each_line { |ln|
1121                       if ln.start_with?(" $")
1122                           keyword = ln.strip
1123                           if keyword !~ /\$END/
1124                               k[keyword] = pos
1125                           end
1126                       end
1127                       pos += ln.length
1128                   }
1129                   #  Cache the information
1130                   @gamess_output_fname = fname
1131                   @gamess_output_mtime = mtime
1132                   @gamess_output_sections = k
1133               end
1134               keywords = k.keys.sort { |a, b| k[a] <=> k[b] }
1135               h = Dialog.run("Select GAMESS section to copy", "Copy", "Cancel") {
1136                   layout(1,
1137                          item(:popup, :subitems=>keywords, :tag=>"section"))
1138               }
1139               if h[:status] == 0
1140                   fp.seek(k[keywords[h["section"]]])
1141                   s = ""
1142                   fp.each_line { |ln|
1143                       s += ln
1144                       break if ln =~ /\$END/
1145                   }
1146                   export_to_clipboard(s)
1147               end
1148               fp.close
1149           end
1150       end
1151   end
1152
1153   def cmd_edit_gamess_input(s)
1154       mol = self
1155       h = Dialog.run("Edit GAMESS Input", "OK", "Cancel", :resizable=>true) {
1156           layout(1,
1157                  item(:textview, :value=>s, :tag=>"edit", :width=>400, :height=>400, :flex=>[0,0,0,0,1,1]),
1158                  item(:button, :title=>"Copy Section from GAMESS Output...", :action=>lambda { |it| mol.copy_section_from_gamess_output } ),
1159                  :flex=>[0,0,0,0,1,1]
1160                  )
1161                  set_min_size(300, 300)
1162       }
1163       if h[:status] == 0
1164           return h["edit"]
1165       else
1166           return nil
1167       end
1168   end
1169
1170   def cmd_create_gamess_input
1171
1172     mol = self
1173
1174     if natoms == 0
1175       raise MolbyError, "cannot create GAMESS input; the molecule is empty"
1176     end
1177
1178         #  Descriptive text and internal string for popup menus
1179 #    bset_desc = ["PM3", "STO-3G", "3-21G", "6-31G", "6-31G(d)", "6-31G(d,p)", "6-311G", "6-311G(d,p)", "LanL2DZ"]
1180 #       bset_internal = ["PM3", "STO3G", "321G", "631G", "631Gd", "631Gdp", "6311G", "6311Gdp", "LanL2DZ"]
1181     bset_desc = $gamess_basis_keys.map { |key| $gamess_basis_desc[key] }
1182         dft_desc = ["B3LYP"]
1183         dft_internal = ["B3LYP"]
1184
1185         defaults = {"scftype"=>0, "runtype"=>0, "use_internal"=>1, "eliminate_freedom"=>1, "charge"=>"0", "mult"=>"1",
1186           "basis"=>4, "use_secondary_basis"=>0, "secondary_elements"=>"",
1187           "secondary_basis"=>8, "esp"=>0, "ncpus"=>"1"}
1188
1189         gamess_input_direct = nil
1190
1191     user_input = Hash.new
1192         ["CONTRL", "SCF", "STATPT", "SYSTEM", "GUESS", "BASIS"].each { |k|
1193           user_input[k] = Hash.new
1194         }
1195     if ".inp".casecmp(self.name[-4, 4]) == 0
1196           #  Read input and analyze commands
1197           fp = open(self.path, "r")
1198           key = hh = nil
1199           if fp
1200             while ln = fp.gets
1201                   ln.strip!
1202                   break if "$DATA".casecmp(ln[0, 5]) == 0
1203                   ln.split.each { |s|
1204                     if s[0] == "$"
1205                           #  New key
1206                           key = s[1..-1].upcase
1207                           hh = user_input[key] = Hash.new
1208                           (user_input["key_order"] ||= []).push(key)
1209                         else
1210                           k, v = s.split("=")
1211                           if key && hh
1212                             k.upcase!
1213                             hh[k] = v
1214                                 (hh["key_order"] ||= []).push(k)
1215                           end
1216                         end
1217                   }
1218                 end  #  end while
1219                 fp.close
1220                 user_input["charge"] = user_input["CONTRL"]["ICHARG"]
1221                 user_input["mult"] = user_input["CONTRL"]["MULT"]
1222                 user_input["runtype"] = ["ENERGY", "PROP", "OPTIMIZE"].find_index(user_input["CONTRL"]["RUNTYP"])
1223                 user_input["scftype"] = ["RHF", "ROHF", "UHF"].find_index(user_input["CONTRL"]["SCFTYP"])
1224                 dft_type = dft_internal.find_index(user_input["CONTRL"]["DFTTYP"])
1225                 if dft_type
1226                   user_input["dfttype"] = dft_type
1227                   user_input["dft"] = 1
1228                 end
1229                 bssname = nil
1230                 user_input["basis"] = "-1"
1231                 case user_input["BASIS"]["GBASIS"]
1232                 when "PM3"
1233                   bssname = "PM3"
1234                 when "STO"
1235                   bssname = "STO3G"
1236                 when "N21"
1237                   bssname = "321G"
1238                 when "N31"
1239                   if user_input["NDFUNC"] == "1"
1240                     if user_input["NPFUNC"] == "1"
1241                           bssname = "631Gdp"
1242                         else
1243                           bssname = "631Gd"
1244                         end
1245                   else
1246                     bssname = "631G"
1247                   end
1248                 when "N311"
1249                   if user_input["NDFUNC"] == "1" && user_input["NPFUNC"] == "1"
1250                     bssname = "6311Gdp"
1251                   else
1252                     bssname = "6311G"
1253                   end
1254                 end
1255                 if bssname
1256                   user_input["basis"] = $gamess_basis_keys.find_index(bssname).to_s
1257                 end
1258           #  puts user_input.inspect
1259           end
1260         end
1261
1262     hash = Dialog.run("GAMESS Export") {
1263       def load_basis_set_sub(item)
1264             fname = Dialog.open_panel("Select a file containing GAMESS basis set:")
1265                 if fname
1266                   Molecule.read_gamess_basis_sets(fname)
1267                   bset_desc_new = $gamess_basis_keys.map { |key| $gamess_basis_desc[key] }
1268                   sel1 = attr("basis", :value)
1269                   sel2 = attr("secondary_basis", :value)
1270                   set_attr("basis", :subitems=>bset_desc_new)
1271                   set_attr("basis", :value=>sel1)
1272                   set_attr("secondary_basis", :subitems=>bset_desc_new)
1273                   set_attr("secondary_basis", :value=>sel2)
1274                 end
1275       end
1276           def select_gamess_path(item)
1277             while 1
1278               fname = Dialog.open_panel("Locate GAMESS executable:")
1279                   return if fname == nil
1280                   bname = File.basename(fname)
1281                   if bname =~ /gamess\.(.*)\.(exe|x)$/i
1282                     set_attr("executable_path", :value=>fname)
1283                     return
1284                   else
1285                     error_message_box("\"#{bname}\" does not look like a GAMESS executable!  Please try again.")
1286                   end
1287                 end
1288           end
1289           def set_optional_scripts(item)
1290             h = Dialog.run("GAMESS Optional Scripts") {
1291                   s_pre = get_global_settings("gamess.prefix_script")
1292                   s_post = get_global_settings("gamess.postfix_script")
1293                   layout(1,
1294                     item(:text, :title=>"Script to run before GAMESS execution:"),
1295                         item(:textview, :width=>400, :height=>200, :value=>s_pre, :tag=>"prefix"),
1296                     item(:text, :title=>"Script to run after GAMESS execution:"),
1297                         item(:textview, :width=>400, :height=>200, :value=>s_post, :tag=>"postfix"))
1298                 }
1299                 if h[:status] == 0
1300                   set_global_settings("gamess.prefix_script", h["prefix"])
1301                   set_global_settings("gamess.postfix_script", h["postfix"])
1302                 end
1303           end
1304       nbos = ["nao", "nbo", "nho", "nlmo", "pnao", "pnbo", "pnho", "pnlmo"]
1305           layout(4,
1306                 item(:text, :title=>"SCF type"),
1307                 item(:popup, :subitems=>["RHF", "ROHF", "UHF"], :tag=>"scftype"),
1308             item(:text, :title=>"Run type"),
1309                 item(:popup, :subitems=>["Energy", "Property", "Optimize", "Sadpoint", "IRC"], :tag=>"runtype",
1310                   :action=>lambda { |it| set_attr("use_internal", :enabled=>(it[:value] >= 2)) } ),
1311         item(:checkbox, :title=>"Use internal coordinates for structure optimization", :tag=>"use_internal"),
1312         -1, -1, -1,
1313                 item(:checkbox, :title=>"Eliminate translation and rotational degrees of freedom", :tag=>"eliminate_freedom"),
1314                 -1, -1, -1,
1315
1316                 item(:text, :title=>"Charge"),
1317                 item(:textfield, :width=>80, :tag=>"charge"),
1318                 item(:text, :title=>"Multiplicity"),
1319                 item(:textfield, :width=>80, :tag=>"mult"),
1320
1321                 item(:checkbox, :title=>"Use DFT", :tag=>"dft",
1322                   :action=>lambda { |it| set_attr("dfttype", :enabled=>(it[:value] != 0)) } ),
1323                 -1,
1324                 item(:text, :title=>"DFT type"),
1325                 item(:popup, :subitems=>dft_desc, :tag=>"dfttype"),
1326
1327                 item(:line),
1328                 -1, -1, -1,
1329
1330                 item(:text, :title=>"Basis set"),
1331                 item(:popup, :subitems=>bset_desc + ["(no select)"], :tag=>"basis"),
1332                 -1,
1333                 -1,
1334
1335                 item(:button, :title=>"Load Basis Set...", :action=>:load_basis_set_sub),
1336                 -1, -1, -1,
1337
1338                 item(:checkbox, :title=>"Use secondary basis set", :tag=>"use_secondary_basis",
1339                   :action=>lambda { |it|
1340                     flag = (it[:value] != 0)
1341                         set_attr("secondary_elements", :enabled=>flag)
1342                         set_attr("secondary_basis", :enabled=>flag)
1343                   }),
1344                 -1, -1, -1,
1345
1346                 item(:text, :title=>"   Elements"),
1347                 item(:textfield, :width=>80, :tag=>"secondary_elements"),
1348                 item(:text, :title=>"Basis set"),
1349                 item(:popup, :subitems=>bset_desc, :tag=>"secondary_basis"),
1350
1351                 item(:line),
1352                 -1, -1, -1,
1353
1354                 item(:checkbox, :title=>"Calculate electrostatic potential (ESP)", :tag=>"esp"),
1355                 -1, -1, -1,
1356
1357         item(:line),
1358         -1, -1, -1,
1359
1360         item(:checkbox, :title=>"Include NBO instructions", :tag=>"include_nbo",
1361              :action=>lambda { |it|
1362                flag = (it[:value] != 0)
1363                nbos.each { |nbo| set_attr(nbo, :enabled=>flag) }
1364              }),
1365         -1, -1, -1,
1366         item(:checkbox, :title=>"NAO", :tag=>"nao"),
1367         item(:checkbox, :title=>"NBO", :tag=>"nbo"),
1368         item(:checkbox, :title=>"NHO", :tag=>"nho"),
1369         item(:checkbox, :title=>"NLMO", :tag=>"nlmo"),
1370         item(:checkbox, :title=>"PNAO", :tag=>"pnao"),
1371         item(:checkbox, :title=>"PNBO", :tag=>"pnbo"),
1372         item(:checkbox, :title=>"PNHO", :tag=>"pnho"),
1373         item(:checkbox, :title=>"PNLMO", :tag=>"pnlmo"),
1374
1375         item(:line),
1376                 -1, -1, -1,
1377
1378                 item(:checkbox, :title=>"Execute GAMESS on this machine", :tag=>"execute_local",
1379                   :action=>lambda { |it|
1380                     flag = (it[:value] != 0)
1381                         set_attr("executable_path", :enabled=>flag)
1382                         set_attr("select_path", :enabled=>flag)
1383                         set_attr("ncpus", :enabled=>flag)
1384                   }),
1385                 -1, -1, -1,
1386
1387                 item(:text, :title=>"   Path"),
1388                 item(:textfield, :width=>300, :tag=>"executable_path"),
1389                 -1, -1,
1390
1391                 -1,
1392                 item(:button, :title=>"Select Path...", :tag=>"select_path", :action=>:select_gamess_path),
1393                 -1,
1394                 item(:button, :title=>"Optional Scripts...", :action=>:set_optional_scripts),
1395
1396                 item(:text, :title=>"   N of CPUs"),
1397                 item(:textfield, :width=>80, :tag=>"ncpus"),
1398                 -1, -1,
1399
1400                 item(:line),
1401                 -1, -1, -1,
1402
1403                 item(:button, :title=>"Edit GAMESS Input and Go", :action=>lambda { |it|
1404                   h = Hash.new
1405                   each_item { |it2|
1406                     if (tag = it2[:tag]) != nil
1407                           h[tag] = it2[:value]
1408                         end
1409                   }
1410                   h["basis"] = $gamess_basis_keys[h["basis"]]
1411                   h["secondary_basis"] = $gamess_basis_keys[h["secondary_basis"]]
1412                   h["dfttype"] = dft_internal[h["dfttype"]]
1413                   gamess_input_direct = mol.cmd_edit_gamess_input(mol.export_gamess(nil, h))
1414                   if gamess_input_direct
1415                         end_modal(0)
1416                   end
1417                 }),
1418                 -1, -1, -1
1419           )
1420           values = Hash.new
1421           each_item { |it|
1422             tag = it[:tag]
1423                 if tag
1424                   values[tag] = (user_input[tag] || get_global_settings("gamess.#{tag}") || defaults[tag])
1425                   if tag == "basis" && values[tag] == "-1"
1426                     values[tag] = (bset_desc.count).to_s
1427                   end
1428                   it[:value] = values[tag]
1429                 end
1430           }
1431           set_attr("secondary_elements", :enabled=>(values["use_secondary_basis"] == 1))
1432           set_attr("secondary_basis", :enabled=>(values["use_secondary_basis"] == 1))
1433           set_attr("dfttype", :enabled=>(values["dft"] == 1))
1434           set_attr("use_internal", :enabled=>(values["runtype"] >= 2))
1435           set_attr("executable_path", :enabled=>(values["execute_local"] == 1))
1436           set_attr("select_path", :enabled=>(values["execute_local"] == 1))
1437           set_attr("ncpus", :enabled=>(values["execute_local"] == 1))
1438       nbos.each { |nao|
1439         set_attr(nao, :enabled=>(values["include_nbo"] == 1))
1440       }
1441         }
1442         hash.each_pair { |key, value|
1443           next if key == :status
1444           set_global_settings("gamess.#{key}", value)
1445         }
1446         if hash[:status] == 0
1447           #  Specify basis by internal keys
1448           hash["basis"] = $gamess_basis_keys[hash["basis"]]
1449           hash["secondary_basis"] = $gamess_basis_keys[hash["secondary_basis"]]
1450           hash["dfttype"] = dft_internal[hash["dfttype"]]
1451           basename = (self.path ? File.basename(self.path, ".*") : self.name)
1452           fname = Dialog.save_panel("Export GAMESS input file:", self.dir, basename + ".inp", "GAMESS input file (*.inp)|*.inp|All files|*.*")
1453           return nil if !fname
1454           if gamess_input_direct
1455           #  puts "gamess_input_direct = \"#{gamess_input_direct}\""
1456             File.open(fname, "w") { |fp| fp.print(gamess_input_direct) }
1457           else
1458             export_gamess(fname, hash)
1459           end
1460           if hash["execute_local"] == 1
1461             Molecule.execute_gamess(fname, self)
1462           end
1463         else
1464           nil
1465         end
1466
1467   end
1468
1469   def cmd_create_cube
1470     grid = default_MO_grid
1471         if grid == nil
1472           Dialog.run {
1473             layout(1,
1474                   item(:text, :title=>"This molecule does not contain MO information."))
1475           }
1476           return
1477         end
1478
1479     mos = selected_MO
1480         if mos == nil || mos.length == 0
1481       Dialog.run {
1482             layout(1,
1483                   item(:text, :title=>"Please select MO(s) in the MO Info table."))
1484       }
1485           return
1486         end
1487         hash = Dialog.run {
1488           layout(1,
1489             item(:text, :title=>"Please specify cube dimensions (in angstrom units):"),
1490             layout(4,
1491                   item(:text, :title=>"Origin"),
1492                   item(:textfield, :width=>100, :height=>20, :tag=>"originx", :value=>sprintf("%.6f", grid[0].x)),
1493                   item(:textfield, :width=>100, :height=>20, :tag=>"originy", :value=>sprintf("%.6f", grid[0].y)),
1494                   item(:textfield, :width=>100, :height=>20, :tag=>"originz", :value=>sprintf("%.6f", grid[0].z)),
1495                   item(:text, :title=>"Delta"),
1496                   item(:textfield, :width=>100, :height=>20, :tag=>"deltax", :value=>sprintf("%.6f", grid[1])),
1497                   item(:textfield, :width=>100, :height=>20, :tag=>"deltay", :value=>sprintf("%.6f", grid[2])),
1498                   item(:textfield, :width=>100, :height=>20, :tag=>"deltaz", :value=>sprintf("%.6f", grid[3])),
1499                   item(:text, :title=>"Step"),
1500                   item(:textfield, :width=>100, :height=>20, :tag=>"stepx", :value=>grid[4].to_s),
1501                   item(:textfield, :width=>100, :height=>20, :tag=>"stepy", :value=>grid[5].to_s),
1502                   item(:textfield, :width=>100, :height=>20, :tag=>"stepz", :value=>grid[6].to_s)))
1503         }
1504         if hash[:status] == 0
1505           path = self.path || self.name
1506           dir = self.dir || Dir.pwd
1507           origin = Vector3D[hash["originx"], hash["originy"], hash["originz"]]
1508           dx = hash["deltax"]
1509           dy = hash["deltay"]
1510           dz = hash["deltaz"]
1511           nx = hash["stepx"]
1512           ny = hash["stepy"]
1513           nz = hash["stepz"]
1514           basename = File.basename(path, ".*")
1515           filenames = []
1516           mo_type = self.mo_type
1517           mos.each { |n|
1518             fname1 = fname2 = nil
1519             alpha = (mo_type != "UHF" ? "" : "alpha ")
1520                 a = (mo_type != "UHF" ? "" : "a")
1521             fname1 = Dialog.save_panel("Cube file name for #{alpha}MO #{n}", dir, basename + "_#{n}#{a}.cube", "Gaussian cube file (*.cube)|*.cube")
1522                 if (mo_type == "UHF")
1523                   fname2 = Dialog.save_panel("Cube file name for beta MO #{n}", dir, basename + "_#{n}b.cube", "Gaussian cube file (*.cube)|*.cube")
1524                 end
1525                 filenames.push([n, fname1, fname2])
1526           }
1527           filenames.each { |pair|
1528             n = pair[0]
1529                 alpha = (mo_type != "UHF" ? "" : "alpha ")
1530             show_progress_panel("Creating cube file for #{alpha}MO #{n}...")
1531                 if pair[1]
1532                   cubegen(pair[1], n, origin, dx, dy, dz, nx, ny, nz, true)
1533                 end
1534                 if pair[2] && mo_type == "UHF"
1535                   set_progress_message("Creating cube file for beta MO #{n}...")
1536                   cubegen(pair[2], n, origin, dx, dy, dz, nx, ny, nz, true, true)
1537                 end
1538             hide_progress_panel
1539           }
1540         end
1541   end
1542
1543   def export_psi4(fname, hash)
1544     now = Time.now.to_s
1545     if fname
1546       fp = File.open(fname, "wb")
1547     else
1548       fp = MemoryIO.new
1549     end
1550     if fp
1551       fp.print "#  Psi4 input\n"
1552       fp.print "#  Generated by Molby at #{now}\n"
1553       fp.print "import sys\n"
1554       fp.print "import re\n"
1555       fp.print "base = re.sub('\\.\\w*$', '', sys.argv[1])  #  Basename of the input file\n"
1556       fp.print "molecule mol {\n"
1557       charge = Integer(hash["charge"])
1558       mult = Integer(hash["mult"])
1559       if charge != 0 || mult != 1
1560         fp.printf " %d %d\n", charge, mult
1561       end
1562       atoms.each { |ap|
1563         fp.printf "%-8s %16.12f %16.12f %16.12f\n", ap.element, ap.x, ap.y, ap.z
1564       }
1565       use_symmetry = Integer(hash["use_symmetry"])
1566       move_to_com = Integer(hash["move_to_com"])
1567       do_reorient = Integer(hash["do_reorient"])
1568       if move_to_com == 0
1569         fp.print "nocom\n"
1570       end
1571       if do_reorient == 0
1572         fp.print "noreorient\n"
1573       end
1574       if use_symmetry == 0
1575         fp.print "symmetry c1\n"
1576       end
1577       fp.print "units angstrom\n"
1578       fp.print "}\n"
1579       fp.print "set basis #{hash['basis']}\n"
1580       fp.print "set reference #{hash['scftype']}\n"
1581       options = "return_wfn=True"
1582       if hash['dft'] != 0
1583         options += ", dft_functional='#{hash['dfttype']}'"
1584       end
1585       runtype = hash['runtype'].downcase
1586       fp.print "energy, wfn = #{runtype}('scf', #{options})\n"
1587       fp.print "wfn.write_molden(f'{base}.molden')\n"
1588       if hash['run_junpa'] != 0
1589         fp.print "\n"
1590         fp.print "#  Interface for JANPA\n"
1591         fp.print "#  cf. https://sourceforge.net/p/janpa/wiki/psi4Examples/\n"
1592         fp.print "d = wfn.Da().to_array()\n"
1593         fp.print "s = wfn.S().to_array()\n"
1594         fp.print "c = wfn.Ca().to_array()\n"
1595         fp.print "occs = c.T.dot(s.dot(d).dot(s).dot(c))\n"
1596         fp.print "molden(wfn, f'{base}.da.molden', density_a = psi4.core.Matrix.from_array(occs))\n"
1597       end
1598     end
1599     if fname == nil
1600       s = fp.buffer
1601       fp.empty
1602       return s
1603     else
1604       fp.close
1605       return fname
1606     end
1607   end
1608
1609   def Molecule.is_java_available
1610     if $platform == "win"
1611       f = get_global_settings("java_home")
1612       if f
1613         ENV["JAVA_HOME"] = f
1614         if !ENV["PATH"].split(";").find { |e| e == "#{f}\\bin" }
1615           ENV["PATH"] = "#{f}\\bin;" + ENV["PATH"]
1616         end
1617       end
1618     end
1619     return (call_subprocess("java -version", nil) == 0)
1620   end
1621
1622   def Molecule.make_java_available
1623     if $platform == "win"
1624       fname = Dialog.open_panel("Locate JDK Folder (if you have one):", "c:\\", nil, true)
1625       return false if fname == nil
1626       fname.sub!(/\//, "\\")
1627       if File.exists?("#{fname}\\bin\\java.exe")
1628         set_global_settings("java_home", fname)
1629         if Molecule.is_java_available()
1630           return true
1631         end
1632       end
1633       error_message_box("Cannot run Java. Please examine your installation again.")
1634       return false
1635     elsif $platform == "mac"
1636       message_box("Please download OpenJDK, and move it into /Library/Java/JavaVirtualMachines folder.", "Install Java", :ok)
1637       return false
1638     else
1639       message_box("Please install Java virtual machine.", "Install Java", :ok)
1640       return false
1641     end
1642   end
1643
1644   #  Execute JANPA
1645   #  inppath is the input file minus extention
1646   #  mol is the molecule (may be nil)
1647   def Molecule.execute_janpa(inppath, mol, spherical)
1648     #    nbo_desc =   #  JANPA
1649     janpa_dir = "#{ResourcePath}/JANPA"
1650     status = 0
1651     outfile = "#{inppath}.janpa.log"
1652     if spherical
1653       cmd1 = ["java", "-jar", "#{janpa_dir}/molden2molden.jar", "-NormalizeBF", "-i", "#{inppath}.da.molden", "-o", "#{inppath}.in.molden"]
1654       cmd2 = nil
1655     else
1656       cmd1 = ["java", "-jar", "#{janpa_dir}/molden2molden.jar", "-frompsi4v1mo", "-NormalizeBF", "-cart2pure", "-i", "#{inppath}.da.molden", "-o", "#{inppath}.in.molden"]
1657       cmd2 = ["java", "-jar", "#{janpa_dir}/molden2molden.jar", "-frompsi4v1mo", "-NormalizeBF", "-cart2pure", "-i", "#{inppath}.molden", "-o", "#{inppath}.spherical.molden"]
1658     end
1659     cmd3 = ["java", "-jar", "#{janpa_dir}/janpa.jar", "-i", "#{inppath}.in.molden"]
1660     ["nao", "pnao", "aho", "lho", "lpo", "clpo"].each { |type|
1661       generate = (get_global_settings("psi4.#{type}").to_i != 0)
1662       if type == "pnao" || type == "nao" || type == "lho"
1663         #  PLHO is generated within Molby from JANPA NAO/PNAO/LHO
1664         generate ||= (get_global_settings("psi4.plho").to_i != 0)
1665       end
1666       if generate
1667         cmd3.push("-#{type.upcase}_Molden_File", "#{inppath}.#{type.upcase}.molden")
1668       end
1669     }
1670     # show_progress_panel("Executing JANPA...")
1671     procname = "molden2molden"
1672     flag = call_subprocess(cmd1, procname, nil, outfile)
1673     if flag && cmd2 != nil
1674       procname = "molden2molden"
1675       flag = call_subprocess(cmd2, procname, nil, ">>#{outfile}")
1676     end
1677     if flag
1678       procname = "janpa"
1679       flag = call_subprocess(cmd3, procname, nil, ">>#{outfile}")
1680     end
1681     if flag
1682       if mol
1683         #  import JANPA log and molden output
1684         #  Files: inppath.janpa.log, inppath.{NAO,PNAO,AHO,LHO,LPO,CLPO}.molden
1685         mol.sub_load_janpa_log(inppath, spherical)
1686       end
1687       hide_progress_panel
1688     else
1689       status = $?.exitstatus
1690       hide_progress_panel
1691       message_box("Execution of #{procname} failed with status #{status}.", "JANPA Failed")
1692     end
1693     return status
1694   end
1695
1696   #  Execute Psi4
1697   #  inpname is the input file
1698   #  mol (optional) is the molecule from which the Psi4 input was built.
1699   def Molecule.execute_psi4(inpname, mol = nil)
1700
1701     inpbase = File.basename(inpname)
1702     inpbody = File.basename(inpname, ".*")
1703     inpdir = File.dirname(inpname)
1704
1705     #  Set PATH for psi4conda
1706     psi4folder = get_global_settings("psi4.psi4conda_folder")
1707     ncpus = get_global_settings("psi4.ncpus").to_i
1708     orgpath = ENV["PATH"]
1709     orgdir = Dir.pwd
1710     if $platform == "win"
1711       ENV["PATH"] = "#{psi4folder};#{psi4folder}\\Library\\mingw-w64\\bin;#{psi4folder}\\Library\\usr\\bin;#{psi4folder}\\Library\\bin;#{psi4folder}\\Scripts;#{psi4folder}\\bin;#{psi4folder}\\condabin;" + ENV["PATH"]
1712     else
1713       ENV["PATH"] = "#{psi4folder}/bin:#{psi4folder}/condabin:" + ENV["PATH"]
1714     end
1715     Dir.chdir(inpdir)
1716     cmdargv = ["psi4", "#{inpbase}"]
1717     if ncpus > 0
1718       cmdargv.push("-n", "#{ncpus}")
1719     end
1720     hf_type = nil
1721     nalpha = nil
1722     nbeta = nil
1723     spherical = false
1724
1725     outfile = inpdir + "/" + inpbody + ".out"
1726     if File.exists?(outfile)
1727       n = 1
1728       while true
1729         outbackfile = inpdir + "/" + inpbody + "~" + (n == 1 ? "" : "#{n}") + ".out"
1730         break if !File.exists?(outbackfile)
1731         n += 1
1732       end
1733       File.rename(outfile, outbackfile)
1734     else
1735       outbackfile = nil
1736     end
1737
1738     #  Timer callback
1739     timer_count = 0
1740     fplog = nil
1741     last_size = 0
1742     lines = []
1743     last_line = ""
1744     next_index = 0
1745     timer_callback = lambda { |m, n|
1746       begin
1747         timer_count += 1
1748         if timer_count < 10   #  Only handle every 1 seconds
1749           return true
1750         end
1751         timer_count = 0
1752         if fplog == nil
1753           if File.exists?(outfile)
1754             fplog = Kernel.open(outfile, "rt")
1755             if fplog == nil
1756               return true
1757             end
1758             last_size = 0
1759           else
1760             return true  #  Skip until outfile is available
1761           end
1762         end
1763         fplog.seek(0, IO::SEEK_END)
1764         current_size = fplog.tell
1765         if current_size > last_size
1766           #  Read new lines
1767           fplog.seek(last_size, IO::SEEK_SET)
1768           fplog.each_line { |line|
1769             if line[-1, 1] == "\n"
1770               lines.push(last_line + line)
1771               last_line = ""
1772             else
1773               last_line += line
1774               break
1775             end
1776           }
1777           last_size = fplog.tell
1778         end
1779         li = next_index
1780         getline = lambda { if li < lines.length; li += 1; return lines[li - 1]; else return nil; end }
1781         vecs = []
1782         while (line = getline.call)
1783           if line =~ /==> Geometry <==/
1784             #  Skip until line containing "------"
1785             while line = getline.call
1786               break if line =~ /------/
1787             end
1788             vecs.clear
1789             index = 0
1790             #  Read atom positions
1791             while line = getline.call
1792               line.chomp!
1793               break if line =~ /^\s*$/
1794               tokens = line.split(' ')
1795               vecs.push(Vector3D[Float(tokens[1]), Float(tokens[2]), Float(tokens[3])])
1796               index += 1
1797             end
1798             if vecs.length < mol.natoms
1799               break  #  Log file is incomplete
1800             end
1801             #  Does this geometry differ from the last one?
1802             vecs.length.times { |i|
1803               if (mol.atoms[i].r - vecs[i]).length2 > 1.0e-20
1804                 #  Create a new frame and break
1805                 mol.create_frame
1806                 vecs.length.times { |j|
1807                   mol.atoms[j].r = vecs[j]
1808                 }
1809                 break
1810               end
1811             }
1812             next_index = li
1813             #  end geometry
1814           elsif line =~ /Final Energy: +([-.0-9]+)/
1815             #  Energy for this geometry
1816             energy = Float($1)
1817             mol.set_property("energy", energy)
1818             mol.show_text("Frame: #{mol.nframes - 1}\nEnergy: #{energy}")
1819             print("Frame: #{mol.nframes - 1} Energy: #{energy}\n")
1820             if line =~ /RHF/
1821               hf_type = "RHF"
1822             elsif line =~ /UHF/
1823               hf_type = "UHF"
1824             elsif line =~ /ROHF/
1825               hf_type = "ROHF"
1826             end
1827             next_index = li
1828           elsif line =~ /^ *Nalpha *= *(\d+)/
1829             nalpha = Integer($1)
1830           elsif line =~ /^ *Nbeta *= *(\d+)/
1831             nbeta = Integer($1)
1832           elsif line =~ /^ *Spherical Harmonics\?: *(\w+)/
1833             spherical = ($1 == "true")
1834           end
1835         end
1836         if next_index > 0
1837           lines.slice!(0, next_index)
1838           next_index = 0
1839         end
1840         return true
1841       rescue => e
1842         #  Some error occurred during callback. Show the message in the console and abort.
1843         $stderr.write("#{e.message}\n")
1844         $stderr.write("#{e.backtrace.inspect}\n")
1845         return false
1846       end
1847     }
1848
1849     #  Terminate callback
1850     term_callback = lambda { |m, n|
1851       do_janpa = false
1852       begin
1853         msg = "Psi4 execution of #{inpbase} "
1854         hmsg = "Psi4 "
1855         if n == -1
1856           msg += "cannot be started. Please examine Psi4 installation."
1857           hmsg += "Error"
1858           icon = :error
1859         else
1860           if n == 0
1861             msg += "succeeded."
1862             hmsg += "Completed"
1863             icon = :info
1864           else
1865             msg += "failed with status #{n}."
1866             hmsg += "Failed"
1867             icon = :error
1868           end
1869           msg += "\n(In directory #{inpdir})"
1870         end
1871         if n == 0
1872           #  Try to load final lines of the logfile
1873           timer_count = 100
1874           timer_callback.call(m, n)
1875           #  Try to load molden file if available
1876           mol.clear_basis_set
1877           mol.clear_mo_coefficients
1878           mol.set_mo_info(:type => hf_type, :alpha => nalpha, :beta => nbeta)
1879           molden = inpdir + "/" + inpbody +".molden"
1880           if File.exists?(molden)
1881             fp = Kernel.open(molden, "rt")
1882             mol.instance_eval { @lineno = 0 }
1883             begin
1884               #  mol.@hf_type should be set before calling sub_load_molden
1885               mol.sub_load_molden(fp)
1886               fp.close
1887             rescue => e
1888               $stderr.write("#{e.message}\n")
1889               $stderr.write("#{e.backtrace.inspect}\n")
1890             end
1891           end
1892           if (get_global_settings("psi4.run_janpa").to_i == 1)
1893             do_janpa = true
1894             Molecule.execute_janpa(inpdir + "/" + inpbody, mol, spherical)
1895           end
1896         elsif n == -1
1897           #  The child process actually did not start
1898           #  Restore the old file if outbackfile is not nil
1899           if outbackfile && !File.exists?(outfile)
1900             File.rename(outbackfile, outfile)
1901           end
1902         end
1903         if outbackfile && File.exists?(outbackfile)
1904           File.delete(outbackfile)
1905         end
1906         ENV["PATH"] = orgpath
1907         Dir.chdir(orgdir)
1908         if mol != nil
1909           message_box(msg, hmsg, :ok, icon)
1910         end
1911       rescue => e
1912         $stderr.write("#{e.message}\n")
1913         $stderr.write("#{e.backtrace.inspect}\n")
1914       end
1915       print("% ")
1916       true
1917     }
1918
1919     if mol
1920       pid = mol.call_subprocess_async(cmdargv, term_callback, timer_callback)
1921       if pid < 0
1922         #  This may not happen on OSX or Linux (don't know for MSW)
1923         error_message_box("Psi4 failed to start. Please examine Psi4 installation.")
1924         return -1
1925       end
1926     else
1927       status = call_subprocess(cmdargv, "Running Psi4")
1928       term_callback.call(nil, status)
1929       return status
1930     end
1931   end
1932
1933   def cmd_edit_psi4_input(s)
1934     mol = self
1935     h = Dialog.run("Edit Psi4 Input", "OK", "Cancel", :resizable=>true) {
1936       layout(1,
1937         item(:textview, :value=>s, :tag=>"edit", :width=>400, :height=>400, :flex=>[0,0,0,0,1,1]),
1938       )
1939       set_min_size(300, 300)
1940     }
1941     if h[:status] == 0
1942         return h["edit"]
1943     else
1944         return nil
1945     end
1946   end
1947
1948   def cmd_create_psi4_input
1949
1950     mol = self
1951
1952     if natoms == 0
1953       raise MolbyError, "cannot create Psi4 input; the molecule is empty"
1954     end
1955
1956     #  Basis sets Cf. https://psicode.org/psi4manual/master/basissets_tables.html#apdx-basistables
1957     bset_desc = ["STO-3G", "3-21G", "6-31G", "6-31G(d)", "6-31G(d,p)", "6-311G", "6-311G(d)", "6-311G(d,p)", "LanL2DZ"]
1958     dfttype_desc = ["B3LYP"]
1959     runtype_desc = ["Energy", "Optimize"]
1960     scftype_desc = ["RHF", "ROHF", "UHF"]
1961     nbo_desc = ["nao", "pnao", "aho", "lho", "lpo", "clpo"]  #  JANPA
1962     user_input = Hash.new
1963     defaults = {"scftype"=>0, "runtype"=>0, "move_to_com"=>0, "do_reorient"=>0,
1964       "use_symmetry"=>0,  "charge"=>"0", "mult"=>"1",
1965       "basis"=>1, "use_secondary_basis"=>0, "secondary_elements"=>"",
1966       "secondary_basis"=>8, "esp"=>0, "ncpus"=>"1"}
1967     psi4_input_direct = nil
1968
1969     hash = Dialog.run("Psi4 Export") {
1970       def load_basis_set_sub(item)
1971         fname = Dialog.open_panel("Select a file containing GAMESS basis set:")
1972         if fname
1973           Molecule.read_gamess_basis_sets(fname)
1974           bset_desc_new = $gamess_basis_keys.map { |key| $gamess_basis_desc[key] }
1975           sel1 = attr("basis", :value)
1976           sel2 = attr("secondary_basis", :value)
1977           set_attr("basis", :subitems=>bset_desc_new)
1978           set_attr("basis", :value=>sel1)
1979           set_attr("secondary_basis", :subitems=>bset_desc_new)
1980           set_attr("secondary_basis", :value=>sel2)
1981         end
1982       end
1983       def select_psi4_folder(item)
1984         #  By default, psi4conda is installed in the home directory
1985         while 1
1986           fname = Dialog.open_panel("Locate 'psi4conda' Folder:", ENV["HOME"] || ENV["USERPROFILE"], nil, true)
1987           return if fname == nil
1988           bname = File.basename(fname)
1989           if bname == "psi4conda"
1990             break
1991           else
1992             h = Dialog.run("", "OK", "Choose Again") {
1993               layout(1,
1994                 item(:text, :title=>"The folder does not have the name 'psi4conda'."),
1995                 item(:text, :title=>"Do you want to use it anyway?"))
1996             }
1997             if h[:status] == 0
1998               break
1999             end
2000           end
2001         end
2002         set_attr("psi4conda_folder", :value=>fname)
2003       end
2004       layout(4,
2005         #  ------
2006         item(:text, :title=>"SCF type"),
2007         item(:popup, :subitems=>scftype_desc, :tag=>"scftype"),
2008         item(:text, :title=>"Run type"),
2009         item(:popup, :subitems=>runtype_desc, :tag=>"runtype"),
2010         #  ------
2011         item(:checkbox, :title=>"Detect symmetry", :tag=>"use_symmetry"),
2012         -1, -1, -1,
2013         #  ------
2014         item(:checkbox, :title=>"Move the molecule to the center of mass", :tag=>"move_to_com"),
2015         -1, -1, -1,
2016         #  ------
2017         item(:checkbox, :title=>"Rotate the molecule to the symmetry principle axes", :tag=>"do_reorient"),
2018         -1, -1, -1,
2019         #  ------
2020         item(:text, :title=>"Charge"),
2021         item(:textfield, :width=>80, :tag=>"charge"),
2022         item(:text, :title=>"Multiplicity"),
2023         item(:textfield, :width=>80, :tag=>"mult"),
2024         #  ------
2025         item(:checkbox, :title=>"Use DFT", :tag=>"dft",
2026           :action=>lambda { |it| set_attr("dfttype", :enabled=>(it[:value] != 0)) } ),
2027         -1,
2028         item(:text, :title=>"DFT type"),
2029         item(:popup, :subitems=>dfttype_desc, :tag=>"dfttype"),
2030         #  ------
2031         item(:line),
2032         -1, -1, -1,
2033         #  ------
2034         item(:text, :title=>"Basis set"),
2035         item(:popup, :subitems=>bset_desc, :tag=>"basis"),
2036         -1,
2037         -1,
2038         #  ------
2039         #item(:button, :title=>"Load Basis Set...", :action=>:load_basis_set_sub),
2040         #-1, -1, -1,
2041         #  ------
2042         #item(:checkbox, :title=>"Use secondary basis set", :tag=>"use_secondary_basis",
2043         #  :action=>lambda { |it|
2044         #    flag = (it[:value] != 0)
2045         #    set_attr("secondary_elements", :enabled=>flag)
2046         #    set_attr("secondary_basis", :enabled=>flag)
2047         #  }),
2048         #-1, -1, -1,
2049         #  ------
2050         #item(:text, :title=>"   Elements"),
2051         #item(:textfield, :width=>80, :tag=>"secondary_elements"),
2052         #item(:text, :title=>"Basis set"),
2053         #item(:popup, :subitems=>bset_desc, :tag=>"secondary_basis"),
2054         #  ------
2055         item(:line),
2056         -1, -1, -1,
2057         #  ------
2058         #item(:checkbox, :title=>"Calculate electrostatic potential (ESP)", :tag=>"esp"),
2059         #-1, -1, -1,
2060         #  ------
2061         #item(:line),
2062         #-1, -1, -1,
2063         #  ------
2064         item(:checkbox, :title=>"Run JANPA after Psi4", :tag=>"run_janpa",
2065             :action=>lambda { |it|
2066               flag = (it[:value] != 0)
2067               nbo_desc.each { |nbo| set_attr(nbo, :enabled=>flag) }
2068             }),
2069         -1, -1, -1,
2070         #  ------
2071         item(:checkbox, :title=>"NAO", :tag=>"nao"),
2072         item(:checkbox, :title=>"PNAO", :tag=>"pnao"),
2073         item(:checkbox, :title=>"AHO", :tag=>"aho"),
2074         item(:checkbox, :title=>"LHO", :tag=>"lho"),
2075         #  ------
2076         item(:checkbox, :title=>"PLHO*", :tag=>"plho"),
2077         item(:checkbox, :title=>"LPO", :tag=>"lpo"),
2078         item(:checkbox, :title=>"CLPO", :tag=>"clpo"),
2079         -1,
2080         #  ------
2081         item(:text, :title=>"* Not JANPA original; Molby extension", :font=>[9]),
2082         -1, -1, -1,
2083         #  ------
2084         item(:line),
2085         -1, -1, -1,
2086         #  ------
2087         item(:checkbox, :title=>"Execute Psi4 on this machine", :tag=>"execute_local",
2088           :action=>lambda { |it|
2089             flag = (it[:value] != 0)
2090             set_attr("psi4conda_folder", :enabled=>flag)
2091             set_attr("select_folder", :enabled=>flag)
2092             set_attr("ncpus", :enabled=>flag)
2093           }),
2094         -1, -1, -1,
2095         #  ------
2096         item(:text, :title=>" Psi4 Folder"),
2097         item(:textfield, :width=>300, :tag=>"psi4conda_folder"),
2098         -1, -1,
2099         #  ------
2100         -1,
2101         item(:button, :title=>"Select Folder...", :tag=>"select_folder", :action=>:select_psi4_folder),
2102         -1, -1,
2103         # item(:button, :title=>"Optional Scripts...", :action=>:set_optional_scripts),
2104         #  ------
2105         item(:text, :title=>"   N of CPUs"),
2106         item(:textfield, :width=>80, :tag=>"ncpus"),
2107         -1, -1,
2108         #  ------
2109         item(:line),
2110         -1, -1, -1,
2111         #  ------
2112         item(:button, :title=>"Edit Psi4 Input and Go", :action=>lambda { |it|
2113           h = Hash.new
2114           each_item { |it2|
2115             if (tag = it2[:tag]) != nil
2116               h[tag] = it2[:value]
2117             end
2118           }
2119           h["basis"] = bset_desc[h["basis"] || 0]
2120           h["secondary_basis"] = bset_desc[h["secondary_basis"] || 0]
2121           h["dfttype"] = dfttype_desc[h["dfttype"] || 0]
2122           h["runtype"] = runtype_desc[h["runtype"] || 0]
2123           h["scftype"] = scftype_desc[h["scftype"] || 0]
2124           psi4_input_direct = mol.cmd_edit_psi4_input(mol.export_psi4(nil, h))
2125           if psi4_input_direct
2126             end_modal(0)
2127           end
2128         }),
2129         -1, -1, -1
2130       )
2131       values = Hash.new
2132       each_item { |it|
2133         tag = it[:tag]
2134         if tag
2135           values[tag] = (user_input[tag] || get_global_settings("psi4.#{tag}") || defaults[tag])
2136           if values[tag]
2137             it[:value] = values[tag]
2138           end
2139         end
2140       }
2141       #set_attr("secondary_elements", :enabled=>(values["use_secondary_basis"] == 1))
2142       #set_attr("secondary_basis", :enabled=>(values["use_secondary_basis"] == 1))
2143       set_attr("dfttype", :enabled=>(values["dft"] == 1))
2144       set_attr("psi4conda_folder", :enabled=>(values["execute_local"] == 1))
2145       set_attr("select_folder", :enabled=>(values["execute_local"] == 1))
2146       set_attr("ncpus", :enabled=>(values["execute_local"] == 1))
2147       #nbos.each { |nao|
2148       #  set_attr(nao, :enabled=>(values["include_nbo"] == 1))
2149       #}
2150     }  #  end Dialog.run
2151
2152     hash.each_pair { |key, value|
2153       next if key == :status
2154       set_global_settings("psi4.#{key}", value)
2155     }
2156     if hash[:status] == 0
2157       #  Specify basis by internal keys
2158       hash["basis"] = bset_desc[hash["basis"] || 0]
2159       hash["secondary_basis"] = bset_desc[hash["secondary_basis"] || 0]
2160       hash["dfttype"] = dfttype_desc[hash["dfttype"] || 0]
2161       hash["runtype"] = runtype_desc[hash["runtype"] || 0]
2162       hash["scftype"] = scftype_desc[hash["scftype"] || 0]
2163       basename = (self.path ? File.basename(self.path, ".*") : self.name)
2164       fname = Dialog.save_panel("Export Psi4 input file:", self.dir, basename + ".in", "Psi4 input file (*.in)|*.in|All files|*.*")
2165       return nil if !fname
2166       if psi4_input_direct
2167         File.open(fname, "w") { |fp| fp.print(psi4_input_direct) }
2168       else
2169         export_psi4(fname, hash)
2170       end
2171       if hash["execute_local"] == 1
2172         if hash["run_janpa"] == 1
2173           #  Check if Java is available
2174           if !Molecule.is_java_available()
2175             if !Molecule.make_java_available()
2176               return nil
2177             end
2178           end
2179         end
2180         @hf_type = hash["scftype"]
2181         Molecule.execute_psi4(fname, self)
2182       end
2183     else
2184       nil
2185     end
2186
2187   end  #  End def create_psi4_input
2188
2189
2190 end
2191
2192 $gamess_basis = {
2193   "PM3"   => " PM3 0\n",
2194   "STO3G" => " STO 3\n",
2195   "321G"  => " N21 3\n",
2196   "631G"  => " N31 6\n" }
2197 $gamess_basis_desc = {
2198   "PM3"   => "PM3",
2199   "STO3G" => "STO-3G",
2200   "321G"  => "3-21G",
2201   "631G"  => "6-31G" }
2202 $gamess_basis_keys = ["PM3", "STO3G", "321G", "631G"]
2203
2204 ["631Gd", "631Gdp", "631+Gd", "631++Gdp", "6311Gdp", "6311+Gd", "6311++Gdp", "6311++G2d2p", "6311++G3df3pd", "LanL2DZ"].each { |n|
2205   Molecule.read_gamess_basis_sets("basis_sets/#{n}.txt")
2206 }
2207
2208 register_menu("QChem\tCreate GAMESS Input...",
2209   :cmd_create_gamess_input, :non_empty)
2210 register_menu("QChem\tCreate Psi4 Input...",
2211   :cmd_create_psi4_input, :non_empty)
2212 register_menu("QChem\tCreate MOPAC6 Input...",
2213   :cmd_create_mopac_input, :non_empty)    # mopac6.rb
2214 register_menu("QChem\tCreate MO Cube...",
2215   :cmd_create_cube, lambda { |m| m && m.mo_type } )