Scripts/gamess.rb

   1 #
   2 #  gamess.rb
   3 #
   4 #  Created by Toshi Nagata on 2009/11/22.
   5 #  Copyright 2009 Toshi Nagata. All rights reserved.
   6 #
   7 # This program is free software; you can redistribute it and/or modify
   8 # it under the terms of the GNU General Public License as published by
   9 # the Free Software Foundation version 2 of the License.
  10 #
  11 # This program is distributed in the hope that it will be useful,
  12 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14 # GNU General Public License for more details.
  15
  16 class Molecule
  17
  18   def Molecule.read_gamess_basis_sets(fname)
  19     # $gamess_basis = Hash.new unless $gamess_basis
  20         $gamess_ecp = Hash.new unless $gamess_ecp
  21         # $gamess_basis_desc = Hash.new unless $gamess_basis_desc
  22         # $gamess_basis_keys = [] unless $gamess_basis_keys
  23         basename = File.basename(fname, ".*")
  24         keys = []
  25         descname = nil
  26     File.open(fname, "r") { |fp|
  27           while (s = fp.gets)
  28             if s =~ /^\s*!\s*/ && descname == nil
  29               #  Get the descriptive name from the first comment line
  30                   s = Regexp.last_match.post_match
  31                   if s =~ /  EMSL/
  32                     descname = Regexp.last_match.pre_match
  33                   else
  34                     descname = (s.split)[0]
  35                   end
  36                   $gamess_basis_desc[basename] = descname
  37                   next
  38                 end
  39                 ss, bas = (s.split)[0..1]     #  Tokens delimited by whitespaces
  40                 next if ss == nil || ss == ""
  41                 if ss == "$ECP"
  42                 #  Read the ECP section
  43                   ecp = $gamess_ecp[basename]
  44                   if !ecp
  45                     ecp = []
  46                         $gamess_ecp[basename] = ecp
  47                   end
  48                   keys << basename unless keys.include?(basename)
  49                   while (s = fp.gets)
  50                     break if s=~ /\$END/
  51                         ary = s.split  #  (PNAME, PTYPE, IZCORE, LMAX+1)
  52                         elem = ary[0].match(/\w{1,2}/).to_a[0]  #  ary[0] should begin with an element name
  53                         ap = Parameter.builtin.elements.find { |p| p.name.casecmp(elem) == 0 }
  54                         raise MolbyError, "the ECP definition does not begin with an element name: #{s}" if !ap
  55                         ecpdef = s
  56                         ln = 1
  57                         (0..Integer(ary[3])).each {
  58                           s = fp.gets
  59                           raise MolbyError, "the ECP definition ends unexpectedly at line #{ln} for: #{s}" if !s
  60                           ln += 1
  61                           ary = s.split  #  (NGPOT, comment)
  62                           ecpdef += s
  63                           (1..Integer(ary[0])).each {
  64                             s = fp.gets
  65                                 raise MolbyError, "the ECP definition ends unexpectedly at line #{ln} for: #{s}" if !s
  66                                 ln += 1
  67                                 ecpdef += s
  68                           }
  69                         }
  70                     ecp[ap.index] = ecpdef
  71                   end
  72                 elsif ss =~ /\W/
  73                   #  Comments or other unrecognizable lines
  74                   next
  75                 elsif (ap = Parameter.builtin.elements.find { |p| p.name.casecmp(ss) == 0 || p.fullname.casecmp(ss) == 0 })
  76                   #  Valid basis definition
  77                   if bas == nil || bas =~ /\W/
  78                     bas = basename
  79                   end
  80                   basis = $gamess_basis[bas]
  81                   if !basis
  82                     basis = []
  83                         $gamess_basis[bas] = basis
  84                   end
  85                   keys << bas unless keys.include?(bas)
  86                   basdef = ""
  87                   while (s = fp.gets) && s =~ /\S/
  88                     basdef += s
  89                   end
  90                   basis[ap.index] = basdef
  91                 else
  92                   raise MolbyError, "ss is not a valid element symbol or name: #{s}"
  93                 end
  94           end
  95     }
  96         unless $gamess_basis_keys.include?(basename)
  97           $gamess_basis_keys.push(basename)
  98         end
  99   end
 100
 101   def export_gamess(fname, hash)
 102
 103         now = Time.now.to_s
 104         basename = File.basename(fname, ".*")
 105
 106         #  Various settings
 107         icharg = hash["charge"]
 108         mult = hash["mult"]
 109         runtyp = ["ENERGY", "PROP", "OPTIMIZE"][Integer(hash["runtype"])]
 110         scftyp = ["RHF", "ROHF", "UHF"][Integer(hash["scftype"])]
 111         bssname = hash["basis"]
 112         bssname2 = hash["secondary_basis"]
 113         if hash["use_secondary_basis"] != 0 && bssname2 != bssname
 114           use_2nd = true
 115           element2 = hash["secondary_elements"].split(/[\s,]+/).map { |name| name.capitalize }
 116         else
 117           use_2nd = false
 118           bssname2 = bssname
 119         end
 120         basis = $gamess_basis[bssname]
 121         basis2 = $gamess_basis[bssname2]
 122         if !basis || !basis2
 123           raise MolbyError, "Unknown basis set name??? \"#{bssname}\" or \"#{bssname2}\""
 124         end
 125
 126         #  Use effective core potentials?
 127         ecp = $gamess_ecp[bssname]
 128         ecp2 = $gamess_ecp[bssname2]
 129         ecp_read = (ecp || ecp2 ? "ECP=READ " : "")
 130
 131         #  Use only one built-in basis set?
 132         gbasis = nil
 133         if !use_2nd
 134           case bssname
 135           when "PM3"
 136                 gbasis = "PM3"
 137           when "STO3G"
 138                 gbasis = "STO NGAUSS=3"
 139           when "321G"
 140                 gbasis = "N21 NGAUSS=3"
 141           when "631G"
 142                 gbasis = "N31 NGAUSS=6"
 143           when "631Gd"
 144                 gbasis = "N31 NGAUSS=6 NDFUNC=1"
 145           when "631Gdp"
 146                 gbasis = "N31 NGAUSS=6 NDFUNC=1 NPFUNC=1"
 147           when "6311G"
 148                 gbasis = "N311 NGAUSS=6"
 149           when "6311Gdp"
 150                 gbasis = "N311 NGAUSS=6 NDFUNC=1 NPFUNC=1"
 151           end
 152         end
 153
 154         #  Count non-dummy atoms
 155         natoms = 0
 156         each_atom { |ap|
 157           natoms += 1 if ap.atomic_number != 0
 158         }
 159
 160     File.open(fname, "wb") { |fp|
 161           fp.print "!  GAMESS input\n"
 162           fp.print "!  Generated by Molby at #{now}\n"
 163           fp.print " $CONTRL COORD=UNIQUE EXETYP=RUN ICHARG=#{icharg}\n"
 164           fp.print "         ICUT=20 INTTYP=HONDO ITOL=30\n"
 165           fp.print "         MAXIT=200 MOLPLT=.T. MPLEVL=0\n"
 166           fp.print "         MULT=#{mult} QMTTOL=1e-08 RUNTYP=#{runtyp}\n"
 167           if hash["dft"] != 0
 168             dfttyp = hash["dfttype"]
 169             fp.print "         DFTTYP=#{dfttyp}\n"
 170           end
 171           if hash["use_internal"] != 0 && hash["runtype"] == 2
 172             nzvar = natoms * 3 - 6  #  TODO: 3N-5 for linear molecules
 173             fp.print "         NZVAR=#{nzvar}\n"
 174           else
 175             nzvar = 0
 176       end
 177           fp.print "         SCFTYP=#{scftyp} #{ecp_read}UNITS=ANGS $END\n"
 178           fp.print " $SCF    CONV=1.0E-06 DIRSCF=.T. FDIFF=.T. DAMP=.T. $END\n"
 179           fp.print " $STATPT NSTEP=400 OPTTOL=1.0E-06               $END\n"
 180           fp.print " $SYSTEM MEMDDI=0 MWORDS=16 TIMLIM=50000        $END\n"
 181           fp.print " $GUESS  GUESS=HUCKEL                           $END\n"
 182       if gbasis
 183             fp.print " $BASIS  GBASIS=#{gbasis} $END\n"
 184           end
 185           if nzvar > 0
 186                 fp.print " $ZMAT   DLC=.T. AUTO=.T. $END\n"
 187           end
 188           if hash["esp"] != 0
 189             fp.print " $ELPOT  IEPOT=1 OUTPUT=PUNCH WHERE=PDC         $END\n"
 190                 fp.print " $PDC    CONSTR=NONE PTSEL=CONNOLLY             $END\n"
 191           end
 192           fp.print " $DATA\n#{basename}\nC1 0\n"
 193           secondary = []
 194           each_atom { |ap|
 195             next if ap.atomic_number == 0
 196                 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
 197                 if use_2nd && element2.include?(ap.element)
 198                   secondary[ap.index] = true
 199                 end
 200             if !gbasis
 201                   #  Basis specification followed by a blank line
 202                   bas = (secondary[ap.index] ? basis2 : basis)
 203                   if bas.is_a?(Array)
 204                     bas = bas[ap.atomic_number]
 205                   end
 206                   if !bas
 207                     raise MolbyError, "Basis set is not defined for atom #{ap.index}, element #{ap.element}"
 208                   end
 209                   fp.print bas
 210                   fp.print "\n"
 211                 end
 212           }
 213           fp.print " $END\n"
 214           ecp_ary = []
 215           if ecp || ecp2
 216             fp.print " $ECP\n"
 217                 each_atom { |ap|
 218                   an = ap.atomic_number
 219                   next if an == 0
 220                   ecpp = (secondary[ap.index] ? ecp2 : ecp)
 221                   e = ecp_ary[an] || (ecpp && ecpp[an])
 222                   if e
 223                     #  Cache the PNAME of the $ECP entry and re-use it
 224                     ecp_ary[an] ||= (e.split)[0] + "\n"
 225                   else
 226                     e = ap.element.upcase + "-ECP NONE\n"
 227                   end
 228                   fp.print e
 229                 }
 230             fp.print " $END\n"
 231           end
 232         }
 233         fname
 234   end
 235
 236   def cmd_create_gamess_input
 237     if natoms == 0
 238       raise MolbyError, "cannot create GAMESS input; the molecule is empty"
 239     end
 240
 241         #  Descriptive text and internal string for popup menus
 242 #    bset_desc = ["PM3", "STO-3G", "3-21G", "6-31G", "6-31G(d)", "6-31G(d,p)", "6-311G", "6-311G(d,p)", "LanL2DZ"]
 243 #       bset_internal = ["PM3", "STO3G", "321G", "631G", "631Gd", "631Gdp", "6311G", "6311Gdp", "LanL2DZ"]
 244     bset_desc = $gamess_basis_keys.map { |key| $gamess_basis_desc[key] }
 245         dft_desc = ["B3LYP"]
 246         dft_internal = ["B3LYP"]
 247
 248         defaults = {"scftype"=>0, "runtype"=>0, "charge"=>"0", "mult"=>"1",
 249           "basis"=>4, "use_secondary_basis"=>0, "secondary_elements"=>"",
 250           "secondary_basis"=>8, "esp"=>0}
 251
 252     hash = Dialog.run("GAMESS Export") {
 253       def load_basis_set_sub(item)
 254             fname = Dialog.open_panel("Select a file containing GAMESS basis set:")
 255                 if fname
 256                   Molecule.read_gamess_basis_sets(fname)
 257                   bset_desc_new = $gamess_basis_keys.map { |key| $gamess_basis_desc[key] }
 258                   sel1 = attr("basis", :value)
 259                   sel2 = attr("secondary_basis", :value)
 260                   set_attr("basis", :subitems=>bset_desc_new)
 261                   set_attr("basis", :value=>sel1)
 262                   set_attr("secondary_basis", :subitems=>bset_desc_new)
 263                   set_attr("secondary_basis", :value=>sel2)
 264                 end
 265       end
 266           layout(4,
 267                 item(:text, :title=>"SCF type"),
 268                 item(:popup, :subitems=>["RHF", "ROHF", "UHF"], :tag=>"scftype"),
 269             item(:text, :title=>"Run type"),
 270                 item(:popup, :subitems=>["Energy", "Property", "Optimize"], :tag=>"runtype",
 271                   :action=>proc { |it| set_attr("use_internal", :enabled=>(it[:value] == 2)) } ),
 272
 273                 item(:checkbox, :title=>"Use internal coordinates for structure optimization", :tag=>"use_internal"),
 274                 -1, -1, -1,
 275
 276                 item(:text, :title=>"Charge"),
 277                 item(:textfield, :width=>80, :tag=>"charge"),
 278                 item(:text, :title=>"Multiplicity"),
 279                 item(:textfield, :width=>80, :tag=>"mult"),
 280
 281                 item(:checkbox, :title=>"Use DFT", :tag=>"dft",
 282                   :action=>proc { |it| set_attr("dfttype", :enabled=>(it[:value] != 0)) } ),
 283                 -1,
 284                 item(:text, :title=>"DFT type"),
 285                 item(:popup, :subitems=>dft_desc, :tag=>"dfttype"),
 286
 287                 item(:line),
 288                 -1, -1, -1,
 289
 290                 item(:text, :title=>"Basis set"),
 291                 item(:popup, :subitems=>bset_desc, :tag=>"basis"),
 292                 -1,
 293                 -1,
 294
 295                 item(:button, :title=>"Load Basis Set...", :action=>:load_basis_set_sub),
 296                 -1, -1, -1,
 297
 298                 item(:checkbox, :title=>"Use secondary basis set", :tag=>"use_secondary_basis",
 299                   :action=>proc { |it|
 300                     flag = (it[:value] != 0)
 301                         set_attr("secondary_elements", :enabled=>flag)
 302                         set_attr("secondary_basis", :enabled=>flag)
 303                   }),
 304                 -1, -1, -1,
 305
 306                 item(:text, :title=>"   Elements"),
 307                 item(:textfield, :width=>80, :tag=>"secondary_elements"),
 308                 item(:text, :title=>"Basis set"),
 309                 item(:popup, :subitems=>bset_desc, :tag=>"secondary_basis"),
 310
 311                 item(:line),
 312                 -1, -1, -1,
 313
 314                 item(:checkbox, :title=>"Calculate electrostatic potential (ESP)", :tag=>"esp"),
 315                 -1, -1, -1,
 316
 317                 item(:line),
 318                 -1, -1, -1
 319           )
 320           values = Hash.new
 321           each_item { |it|
 322             tag = it[:tag]
 323                 if tag
 324                   values[tag] = (get_global_settings("gamess.#{tag}") || defaults[tag])
 325                   it[:value] = values[tag]
 326                 end
 327           }
 328           set_attr("secondary_elements", :enabled=>(values["use_secondary_basis"] == 1))
 329           set_attr("secondary_basis", :enabled=>(values["use_secondary_basis"] == 1))
 330           set_attr("dfttype", :enabled=>(values["dft"] == 1))
 331           set_attr("use_internal", :enabled=>(values["runtype"] == 2))
 332         }
 333         hash.each_pair { |key, value|
 334           next if key == :status
 335           set_global_settings("gamess.#{key}", value)
 336         }
 337         if hash[:status] == 0
 338           #  Specify basis by internal keys
 339           hash["basis"] = $gamess_basis_keys[hash["basis"]]
 340           hash["secondary_basis"] = $gamess_basis_keys[hash["secondary_basis"]]
 341           hash["dfttype"] = dft_internal[hash["dfttype"]]
 342           basename = (self.path ? File.basename(self.path, ".*") : self.name)
 343       fname = Dialog.save_panel("GAMESS input file name", self.dir, basename + ".inp", "GAMESS input file (*.inp)|*.inp|All files|*.*")
 344           return nil if !fname
 345           export_gamess(fname, hash)
 346         else
 347           nil
 348         end
 349   end
 350
 351 end
 352
 353 $gamess_basis = {
 354   "PM3"   => " PM3 0\n",
 355   "STO3G" => " STO 3\n",
 356   "321G"  => " N21 3\n",
 357   "631G"  => " N31 6\n" }
 358 $gamess_basis_desc = {
 359   "PM3"   => "PM3",
 360   "STO3G" => "STO-3G",
 361   "321G"  => "3-21G",
 362   "631G"  => "6-31G" }
 363 $gamess_basis_keys = ["PM3", "STO3G", "321G", "631G"]
 364
 365 ["631Gd", "631Gdp", "631+Gd", "631++Gdp", "6311Gdp", "6311+Gd", "6311++Gdp", "6311++G2d2p", "6311++G3df3pd", "LanL2DZ"].each { |n|
 366   Molecule.read_gamess_basis_sets("basis_sets/#{n}.txt")
 367 }