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