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