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.atoms.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.atoms.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     File.open(fname, "wb") { |fp|
 138           fp.print "!  GAMESS input\n"
 139           fp.print "!  Generated by Molby at #{now}\n"
 140           fp.print " $CONTRL COORD=UNIQUE EXETYP=RUN ICHARG=#{icharg}\n"
 141           fp.print "         ICUT=20 INTTYP=HONDO ITOL=30\n"
 142           fp.print "         MAXIT=200 MOLPLT=.T. MPLEVL=0\n"
 143           fp.print "         MULT=#{mult} QMTTOL=1e-08 RUNTYP=#{runtyp}\n"
 144           if hash["dft"] != 0
 145             dfttyp = hash["dfttype"]
 146             fp.print "         DFTTYP=#{dfttyp}\n"
 147           end
 148           fp.print "         SCFTYP=#{scftyp} #{ecp_read}UNITS=ANGS $END\n"
 149           fp.print " $SCF    CONV=1.0E-06 DIRSCF=.T. FDIFF=.T. DAMP=.T. $END\n"
 150           fp.print " $STATPT NSTEP=400 OPTTOL=1.0E-06               $END\n"
 151           fp.print " $SYSTEM MEMDDI=0 MWORDS=16 TIMLIM=50000        $END\n"
 152           fp.print " $GUESS  GUESS=HUCKEL                           $END\n"
 153       if gbasis
 154             fp.print " $BASIS  GBASIS=#{gbasis} $END\n"
 155           end
 156           if hash["esp"] != 0
 157             fp.print " $ELPOT  IEPOT=1 OUTPUT=PUNCH WHERE=PDC         $END\n"
 158                 fp.print " $PDC    CONSTR=NONE PTSEL=CONNOLLY             $END\n"
 159           end
 160           fp.print " $DATA\n#{basename}\nC1 0\n"
 161           secondary = []
 162           each_atom { |ap|
 163                 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
 164                 if use_2nd && element2.include?(ap.element)
 165                   secondary[ap.index] = true
 166                 end
 167             if !gbasis
 168                   #  Basis specification followed by a blank line
 169                   bas = (secondary[ap.index] ? basis2 : basis)
 170                   if bas.is_a?(Array)
 171                     bas = bas[ap.atomic_number]
 172                   end
 173                   if !bas
 174                     raise MolbyError, "Basis set is not defined for atom #{ap.index}, element #{ap.element}"
 175                   end
 176                   fp.print bas
 177                   fp.print "\n"
 178                 end
 179           }
 180           fp.print " $END\n"
 181           ecp_ary = []
 182           if ecp || ecp2
 183             fp.print " $ECP\n"
 184                 each_atom { |ap|
 185                   an = ap.atomic_number
 186                   ecpp = (secondary[ap.index] ? ecp2 : ecp)
 187                   e = ecp_ary[an] || (ecpp && ecpp[an])
 188                   if e
 189                     #  Cache the PNAME of the $ECP entry and re-use it
 190                     ecp_ary[an] ||= (e.split)[0] + "\n"
 191                   else
 192                     e = ap.element.upcase + "-ECP NONE\n"
 193                   end
 194                   fp.print e
 195                 }
 196             fp.print " $END\n"
 197           end
 198         }
 199         fname
 200   end
 201
 202   def cmd_create_gamess_input
 203     if natoms == 0
 204       raise MolbyError, "cannot create GAMESS input; the molecule is empty"
 205     end
 206
 207         #  Descriptive text and internal string for popup menus
 208     bset_desc = ["PM3", "STO-3G", "3-21G", "6-31G", "6-31G(d)", "6-31G(d,p)", "6-311G", "6-311G(d,p)", "LanL2DZ"]
 209         bset_internal = ["PM3", "STO3G", "321G", "631G", "631Gd", "631Gdp", "6311G", "6311Gdp", "LanL2DZ"]
 210         dft_desc = ["B3LYP"]
 211         dft_internal = ["B3LYP"]
 212
 213         defaults = {"scftype"=>0, "runtype"=>0, "charge"=>"0", "mult"=>"1",
 214           "basis"=>4, "use_secondary_basis"=>0, "secondary_elements"=>"",
 215           "secondary_basis"=>8, "esp"=>0}
 216
 217     hash = RubyDialog.run("GAMESS Export") {
 218           def action(n)
 219             if n == 0 || n == 1
 220               each_item { |i|
 221                 tag = attr(i, :tag)
 222                         if tag
 223                           set_global_settings("gamess.#{tag}", attr(i, :value))
 224                         end
 225                   }
 226                 end
 227                 super
 228           end
 229           layout(4,
 230                 item(:text, :title=>"SCF type"),
 231                 item(:popup, :subitems=>["RHF", "ROHF", "UHF"], :tag=>"scftype"),
 232             item(:text, :title=>"Run type"),
 233                 item(:popup, :subitems=>["Energy", "Property", "Optimize"], :tag=>"runtype"),
 234
 235                 item(:text, :title=>"Charge"),
 236                 item(:textfield, :width=>80, :tag=>"charge"),
 237                 item(:text, :title=>"Multiplicity"),
 238                 item(:textfield, :width=>80, :tag=>"mult"),
 239
 240                 item(:checkbox, :title=>"Use DFT", :tag=>"dft",
 241                   :action=>proc { |n| set_attr("dfttype", :enabled=>(attr(n, :value) != 0)) } ),
 242                 -1,
 243                 item(:text, :title=>"DFT type"),
 244                 item(:popup, :subitems=>dft_desc, :tag=>"dfttype"),
 245
 246                 item(:line),
 247                 -1, -1, -1,
 248
 249                 item(:text, :title=>"Basis set"),
 250                 item(:popup, :subitems=>bset_desc, :tag=>"basis"),
 251                 -1, -1,
 252
 253                 item(:checkbox, :title=>"Use secondary basis set", :tag=>"use_secondary_basis",
 254                   :action=>proc { |n|
 255                     flag = (attr(n, :value) != 0)
 256                         set_attr("secondary_elements", :enabled=>flag)
 257                         set_attr("secondary_basis", :enabled=>flag)
 258                   }),
 259                 -1, -1, -1,
 260
 261                 item(:text, :title=>"   Elements"),
 262                 item(:textfield, :width=>80, :tag=>"secondary_elements"),
 263                 item(:text, :title=>"Basis set"),
 264                 item(:popup, :subitems=>bset_desc, :tag=>"secondary_basis"),
 265
 266                 item(:line),
 267                 -1, -1, -1,
 268
 269                 item(:checkbox, :title=>"Calculate electrostatic potential (ESP)", :tag=>"esp"),
 270                 -1, -1, -1,
 271
 272                 item(:line),
 273                 -1, -1, -1
 274           )
 275           values = Hash.new
 276           each_item { |i|
 277             tag = attr(i, :tag)
 278                 if tag
 279                   values[tag] = (get_global_settings("gamess.#{tag}") || defaults[tag])
 280                   set_attr(i, :value=>values[tag])
 281                 end
 282           }
 283           set_attr("secondary_elements", :enabled=>(values["use_secondary_basis"] == 1))
 284           set_attr("secondary_basis", :enabled=>(values["use_secondary_basis"] == 1))
 285           set_attr("dfttype", :enabled=>(values["dft"] == 1))
 286         }
 287         if hash
 288           #  Specify basis by internal keys
 289           hash["basis"] = bset_internal[hash["basis"]]
 290           hash["secondary_basis"] = bset_internal[hash["secondary_basis"]]
 291           hash["dfttype"] = dft_internal[hash["dfttype"]]
 292           basename = (self.path ? File.basename(self.path, ".*") : self.name)
 293       fname = RubyDialog.save_panel("GAMESS input file name", self.dir, basename + ".inp", "GAMESS input file (*.inp)|*.inp|All files|*.*")
 294           return nil if !fname
 295           export_gamess(fname, hash)
 296         else
 297           nil
 298         end
 299   end
 300
 301 end
 302
 303 Molecule.read_gamess_basis_sets("LanL2DZ.txt")
 304 Molecule.read_gamess_basis_sets("631Gd.txt")
 305 Molecule.read_gamess_basis_sets("631Gdp.txt")
 306 Molecule.read_gamess_basis_sets("6311Gdp.txt")
 307 $gamess_basis["PM3"]   = " PM3 0\n"
 308 $gamess_basis["STO3G"] = " STO 3\n"
 309 $gamess_basis["321G"]  = " N21 3\n"
 310 $gamess_basis["631G"]  = " N31 6\n"