Scripts/md.rb

   1 #
   2 #  md.rb
   3 #
   4 #  Created by Toshi Nagata.
   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 cmd_md_sub
  19     arena = self.md_arena
  20     keys = arena.to_hash.keys
  21         #  Sort the keys so that they are (a little more) readable
  22         [:timestep, :temperature, :cutoff, :electro_cutoff, :pairlist_distance,
  23          :scale14_vdw, :scale14_elect, :use_xplor_shift, :dielectric,
  24          :andersen_freq, :andersen_coupling, :random_seed, :relocate_center,
  25          :use_graphite, :surface_probe_radius, :surface_tension, :surface_potential_freq,
  26          :gradient_convergence, :coordinate_convergence,
  27          :log_file, :coord_file, :vel_file, :force_file, :debug_file, :debug_output_level,
  28          :coord_output_freq, :energy_output_freq].each_with_index { |k, i|
  29           keys.delete(k)
  30           keys.insert(i, k)
  31         }
  32         #  Arrange the items in vertical direction
  33         n = (keys.count + 1) / 2
  34         i = 0
  35         keys = keys.sort_by { |k| i += 1; (i > n ? i - n + 0.5 : i) }
  36         #  Do dialog
  37     hash = RubyDialog.run("Molecular Dynamics Advanced Settings") {
  38           items = []
  39           keys.each { |k|
  40             enabled = (k == :step || k == :coord_frame ? false : true)
  41             items.push(item(:text, :title=>k.to_s))
  42                 it = item(:textfield, :width=>120, :value=>arena[k].to_s)
  43                 if enabled
  44                   set_attr(it, :tag=>k)
  45                 else
  46                   set_attr(it, :enabled=>false)
  47                 end
  48                 items.push(it)
  49           }
  50           layout(4, *items)
  51         }
  52         if hash
  53           hash.keys.each { |k|
  54             arena[k] = hash[k]
  55           }
  56 #         arena.prepare
  57           return hash
  58         else
  59           return nil
  60         end
  61   end
  62
  63   def prepare_arena
  64
  65         #  Get the MDArena
  66         arena = self.md_arena
  67         if !arena.prepare(true)   #  Parameter check only
  68           RubyDialog.run("MD Error") {
  69                 layout(1,
  70                   item(:text, :title=>"Some parameters are missing. Please open\nthe 'parameters' table and examine."))
  71                   set_attr(1, :hidden=>true)  #  Hide the cancel button
  72           }
  73           return nil
  74     end
  75
  76     #  Initialize some fields at first invocation
  77     if !@arena_created
  78           if arena.log_file == nil && self.dir != nil
  79             arena.log_file = self.dir + "/" + self.name.gsub(/\.\w+$/, ".log")
  80           end
  81           @md_spf = 500
  82           @md_nf = 200
  83           @arena_created = true
  84         end
  85
  86     return arena
  87
  88   end
  89
  90   def cmd_md(minimize)
  91
  92         #  Create arena
  93         arena = self.prepare_arena
  94     if !arena
  95           return -1
  96         end
  97
  98     #  Basic parameters
  99         spf = @md_spf
 100         nf = @md_nf
 101         mol = self
 102         files_save = Hash.new
 103         [:log_file, :coord_file, :vel_file, :force_file, :debug_file].each { |k|
 104           s = arena[k]
 105           files_save[k] = s
 106           if s != nil
 107             arena[k] = File.basename(s)
 108           end
 109         }
 110         title = (minimize ? "Minimize" : "Molecular Dynamics")
 111         hash = RubyDialog.run(title) {
 112           items = []
 113           if !minimize
 114             items.push item(:text, :title=>"Timestep (fs)")
 115             items.push item(:textfield, :width=>120, :value=>arena.timestep.to_s, :tag=>"timestep")
 116             items.push item(:text, :title=>"Target temperature (K)")
 117             items.push item(:textfield, :width=>120, :value=>arena.temperature.to_s, :tag=>"temperature")
 118           end
 119           items.push item(:text, :title=>"Steps per frame")
 120           items.push item(:textfield, :width=>120, :value=>arena.coord_output_freq.to_s, :tag=>"steps_per_frame")
 121           items.push item(:text, :title=>"Number of frames")
 122           items.push item(:textfield, :width=>120, :value=>nf.to_s, :tag=>"number_of_frames")
 123           items.push item(:text, :title=>"Log file")
 124           items.push item(:textfield, :width=>120, :value=>(arena.log_file || ""), :tag=>"log_file")
 125           items.push item(:button, :title=>"Advanced...",
 126                      :action=>proc {
 127                            if mol.cmd_md_sub
 128                              if !minimize
 129                                set_value("timestep", arena[:timestep].to_s)
 130                                    set_value("temperature", arena[:temperature].to_s)
 131                                  end
 132                                  set_value("log_file", arena[:log_file])
 133                            end
 134                          }
 135                         )
 136           layout(2, *items)
 137         }
 138         dirstr = (self.dir || document_home)
 139         if hash
 140           arena[:log_file] = hash[:log_file]
 141         end
 142         [:log_file, :coord_file, :vel_file, :force_file, :debug_file].each { |k|
 143           if hash
 144             s = arena[k]
 145             if s != nil && s != ""
 146               arena[k] = dirstr + "/" + s
 147                 else
 148                   arena[k] = nil
 149                 end
 150           else
 151             arena[k] = files_save[k]
 152       end
 153         }
 154         if !hash
 155           return -1
 156         end
 157         if !minimize
 158           arena.temperature = Float(hash["temperature"])
 159           arena.timestep = Float(hash["timestep"])
 160         end
 161         arena.coord_output_freq = Integer(hash["steps_per_frame"])
 162         arena.energy_output_freq = arena.coord_output_freq
 163         return Integer(hash["number_of_frames"])
 164   end
 165
 166   def cmd_define_unit_cell
 167     mol = self
 168     hash = RubyDialog.run("Define Unit Cell") {
 169           @mol = mol
 170           def set_box_value(n)
 171             h = Hash.new
 172                 ["o0", "o1", "o2", "a0", "a1", "a2", "b0", "b1", "b2", "c0", "c1", "c2"].each { |k|
 173                   begin
 174                     s = value(k)
 175                         if s == nil || s == ""
 176                           h[k] = 0.0
 177                         else
 178                       h[k] = Float(eval(s))
 179                           set_value(k, h[k].to_s)
 180                         end
 181                   rescue
 182                     mes = "Cannot evaluate #{value(k)}: " + $!.to_s
 183                         RubyDialog.run("Value Error") {
 184                           layout(1, item(:text, :title=>mes))
 185                           set_attr(1, :hidden=>true)
 186                         }
 187                         return nil
 188                   end
 189                 }
 190                 ax = Vector3D[h["a0"], h["a1"], h["a2"]]
 191                 bx = Vector3D[h["b0"], h["b1"], h["b2"]]
 192                 cx = Vector3D[h["c0"], h["c1"], h["c2"]]
 193                 ox = Vector3D[h["o0"], h["o1"], h["o2"]]
 194                 @mol.set_box(ax, bx, cx, ox)
 195                 return @mol
 196           end
 197           def action(n)
 198             if n == 0
 199                   if !set_box_value(n)
 200                     return  #  Cannot set box: dialog is not dismissed
 201                   end
 202                 end
 203                 super
 204           end
 205           box = @mol.box
 206           layout(4,
 207             item(:text, :title=>"Unit cell:"),
 208                 -1, -1, -1,
 209             item(:text, :title=>"origin"),
 210                 item(:textfield, :width=>140, :tag=>"o0", :value=>(box ? box[3].x.to_s : "")),
 211                 item(:textfield, :width=>140, :tag=>"o1", :value=>(box ? box[3].y.to_s : "")),
 212                 item(:textfield, :width=>140, :tag=>"o2", :value=>(box ? box[3].z.to_s : "")),
 213             item(:text, :title=>"a-axis"),
 214                 item(:textfield, :width=>140, :tag=>"a0", :value=>(box ? box[0].x.to_s : "")),
 215                 item(:textfield, :width=>140, :tag=>"a1", :value=>(box ? box[0].y.to_s : "")),
 216                 item(:textfield, :width=>140, :tag=>"a2", :value=>(box ? box[0].z.to_s : "")),
 217             item(:text, :title=>"b-axis"),
 218                 item(:textfield, :width=>140, :tag=>"b0", :value=>(box ? box[1].x.to_s : "")),
 219                 item(:textfield, :width=>140, :tag=>"b1", :value=>(box ? box[1].y.to_s : "")),
 220                 item(:textfield, :width=>140, :tag=>"b2", :value=>(box ? box[1].z.to_s : "")),
 221             item(:text, :title=>"c-axis"),
 222                 item(:textfield, :width=>140, :tag=>"c0", :value=>(box ? box[2].x.to_s : "")),
 223                 item(:textfield, :width=>140, :tag=>"c1", :value=>(box ? box[2].y.to_s : "")),
 224                 item(:textfield, :width=>140, :tag=>"c2", :value=>(box ? box[2].z.to_s : "")),
 225                 item(:button, :title=>"Set", :action=>:set_box_value),
 226                 item(:text, :title=>"(Ruby expressions are allowed as the values)"),
 227                 -1, -1)
 228         }
 229   end
 230
 231   def cmd_show_periodic_image
 232     mol = self
 233     hash = RubyDialog.run("Show Periodic Image") {
 234           @mol = mol
 235           def set_periodic_image(n)
 236             a = []
 237             ["amin", "amax", "bmin", "bmax", "cmin", "cmax"].each_with_index { |k, i|
 238                   s = value(k)
 239                   if s == nil || s == ""
 240                     a[i] = 0
 241                   else
 242                     a[i] = Integer(s)
 243                   end
 244                 }
 245             @mol.show_periodic_image(a)
 246           end
 247           def action(n)
 248             if n == 0
 249                   set_periodic_image(n)
 250                 end
 251                 super
 252           end
 253           pimage = @mol.show_periodic_image
 254           layout(4,
 255             item(:text, :title=>"Show Periodic Image:"),
 256                 -1, -1, -1,
 257                 item(:text, :title=>"a-axis"),
 258                 item(:textfield, :width=>80, :tag=>"amin", :value=>pimage[0].to_s),
 259                 item(:text, :title=>"to"),
 260                 item(:textfield, :width=>80, :tag=>"amax", :value=>pimage[1].to_s),
 261                 item(:text, :title=>"b-axis"),
 262                 item(:textfield, :width=>80, :tag=>"bmin", :value=>pimage[2].to_s),
 263                 item(:text, :title=>"to"),
 264                 item(:textfield, :width=>80, :tag=>"bmax", :value=>pimage[3].to_s),
 265                 item(:text, :title=>"c-axis"),
 266                 item(:textfield, :width=>80, :tag=>"cmin", :value=>pimage[4].to_s),
 267                 item(:text, :title=>"to"),
 268                 item(:textfield, :width=>80, :tag=>"cmax", :value=>pimage[5].to_s),
 269                 item(:button, :title=>"Set", :action=>:set_periodic_image))
 270         }
 271   end
 272
 273   def cmd_pressure_control
 274
 275     if box == nil
 276           RubyDialog.run("Pressure Control Error") {
 277                 layout(1,
 278                   item(:text, :title=>"Unit cell is not defined. Please open 'Unit Cell...' \nfrom the MM/MD menu and set up the unit cell."))
 279                 set_attr(1, :hidden=>true)  #  Hide the cancel button
 280           }
 281           return nil
 282         end
 283
 284         #  Create arena
 285         arena = self.prepare_arena
 286     if !arena
 287           return nil
 288         end
 289     mol = self
 290
 291         #  Create pressure arena
 292         if arena.pressure_freq == nil
 293           arena.pressure_freq = 20   #  This assignment will automatically create pressure arena
 294         end
 295
 296         hash = RubyDialog.run("Pressure Control") {
 297           @mol = mol
 298           layout(2,
 299             item(:text, :title=>"Frequency"),
 300                 item(:textfield, :width=>100, :tag=>"freq", :value=>arena.pressure_freq.to_s),
 301                 item(:text, :title=>"Coupling"),
 302                 item(:textfield, :width=>100, :tag=>"coupling", :value=>arena.pressure_coupling.to_s),
 303             item(:text, :title=>"Pressure a"),
 304                 item(:textfield, :width=>100, :tag=>"pa", :value=>arena.pressure[0].to_s),
 305             item(:text, :title=>"Pressure b"),
 306                 item(:textfield, :width=>100, :tag=>"pb", :value=>arena.pressure[1].to_s),
 307             item(:text, :title=>"Pressure c"),
 308                 item(:textfield, :width=>100, :tag=>"pc", :value=>arena.pressure[2].to_s),
 309             item(:text, :title=>"Cell flexibility a"),
 310                 item(:textfield, :width=>100, :tag=>"fa", :value=>arena.pressure_cell_flexibility[0].to_s),
 311             item(:text, :title=>"Cell flexibility b"),
 312                 item(:textfield, :width=>100, :tag=>"fb", :value=>arena.pressure_cell_flexibility[1].to_s),
 313             item(:text, :title=>"Cell flexibility c"),
 314                 item(:textfield, :width=>100, :tag=>"fc", :value=>arena.pressure_cell_flexibility[2].to_s),
 315             item(:text, :title=>"Cell flexibility bc"),
 316                 item(:textfield, :width=>100, :tag=>"fbc", :value=>arena.pressure_cell_flexibility[3].to_s),
 317             item(:text, :title=>"Cell flexibility ca"),
 318                 item(:textfield, :width=>100, :tag=>"fca", :value=>arena.pressure_cell_flexibility[4].to_s),
 319             item(:text, :title=>"Cell flexibility ab"),
 320                 item(:textfield, :width=>100, :tag=>"fab", :value=>arena.pressure_cell_flexibility[5].to_s),
 321             item(:text, :title=>"Cell flexibility origin"),
 322                 item(:textfield, :width=>100, :tag=>"forig", :value=>arena.pressure_cell_flexibility[6].to_s),
 323             item(:text, :title=>"Cell flexibility orientation"),
 324                 item(:textfield, :width=>100, :tag=>"forient", :value=>arena.pressure_cell_flexibility[7].to_s),
 325             item(:text, :title=>"Fluctuate cell origin"),
 326                 item(:textfield, :width=>100, :tag=>"orig", :value=>arena.pressure_fluctuate_cell_origin.to_s),
 327             item(:text, :title=>"Fluctuate cell orientation"),
 328                 item(:textfield, :width=>100, :tag=>"orient", :value=>arena.pressure_fluctuate_cell_orientation.to_s))
 329     }
 330         if hash
 331           arena.pressure_freq = hash["freq"]
 332           arena.pressure_coupling = hash["coupling"]
 333           arena.pressure = [hash["pa"], hash["pb"], hash["pc"]]
 334           arena.pressure_cell_flexibility = [hash["fa"], hash["fb"], hash["fc"], hash["fbc"], hash["fca"], hash["fab"], hash["forig"], hash["forient"]]
 335           arena.pressure_fluctuate_cell_origin = hash["orig"]
 336           arena.pressure_fluctuate_cell_orientation = hash["orient"]
 337         end
 338   end
 339
 340   def ambertools_dialog(tool)
 341         ante_dir = get_global_settings("antechamber.ante_dir")
 342         log_dir = get_global_settings("antechamber.log_dir")
 343         if !ante_dir
 344           ante_dir = ($platform == "mac" ? "/Applications/amber10" : ($platform == "win" ? "c:/opt/amber10" : ""))
 345         end
 346         if !log_dir
 347           log_dir = document_home + "/amber10"
 348         end
 349         if $platform == "win"
 350           suffix = ".exe"
 351         else
 352           suffix = ""
 353         end
 354         log_level = (get_global_settings("antechamber.log_level") || "none")
 355     hash = RubyDialog.run("Run " + tool.capitalize) {
 356           @toolname = tool + suffix
 357       def valid_antechamber_dir(s)
 358             FileTest.exist?(s + "/" + @toolname)
 359           end
 360           def action(n)
 361             if n == 4
 362                   if valid_antechamber_dir(attr(n, :value))
 363                     set_attr(0, :enabled=>true)
 364                   end
 365                 elsif n == 0 || n == 1
 366                   #  Settings are stored in the global settings (even if Cancel is pressed)
 367                   self.each_item { |i|
 368                     if (tag = attr(i, :tag))
 369                           value = attr(i, :value)
 370                   v = [["log_none", "none"], ["log_error_only", "error_only"], ["log_keep_latest", "latest"], ["log_all", "all"]].assoc(tag)
 371                           if v
 372                             next if value != 1
 373                             tag = "log_level"
 374                                 value = v[1]
 375                           end
 376                           set_global_settings("antechamber.#{tag}", value)
 377                         end
 378                   }
 379                 end
 380                 super
 381           end
 382           layout(2,
 383                 item(:text, :title=>"Ambertools directory:"),
 384                 [ item(:button, :title=>"Choose...",
 385                         :action=>proc { |n|
 386                         #  print "action method called\n"
 387                           dir = RubyDialog.open_panel(nil, nil, nil, true)
 388                           if dir
 389                                 if valid_antechamber_dir(dir)
 390                                   set_value("ante_dir", dir)
 391                           set_attr(0, :enabled=>true)
 392                                 else
 393                                   RubyDialog.run {
 394                                     layout(1,
 395                                           item(:text, :title=>"Cannot find #{tool} in #{dir}."))
 396                                   }
 397                                   set_attr(0, :enabled=>false)
 398                                 end
 399                           end
 400                         }
 401                   ), {:align=>:right} ],
 402                 item(:textfield, :width=>360, :height=>40, :tag=>"ante_dir", :value=>ante_dir),
 403                 -1,
 404                 item(:text, :title=>"Log directory:"),
 405                 [ item(:button, :title=>"Choose...",
 406                         :action=>proc { |n|
 407                           dir = RubyDialog.open_panel(nil, nil, nil, true)
 408                           if dir
 409                                 set_value("log_dir", dir)
 410                           end
 411                         }
 412                   ), {:align=>:right} ],
 413                 item(:textfield, :width=>360, :height=>40, :tag=>"log_dir", :value=>log_dir),
 414                 -1,
 415                 item(:text, :title=>"Log handling"),
 416                 -1,
 417                 layout(1,
 418                   item(:radio, :title=>"Do not keep logs", :tag=>"log_none", :value=>(log_level == "none" ? 1 : 0)),
 419                   item(:radio, :title=>"Keep only when error occurred", :tag=>"log_error_only", :value=>(log_level == "error_only" ? 1 : 0)),
 420                   layout(3,
 421                     item(:radio, :title=>"Keep latest ", :tag=>"log_keep_latest", :value=>(log_level == "latest" ? 1 : 0)),
 422                     item(:textfield, :width=>80, :tag=>"log_keep_number", :value=>(get_global_settings("antechamber.log_keep_number") || "10")),
 423                     item(:text, :title=>"logs"),
 424                     {:margin=>0, :padding=>0}),
 425                   item(:radio, :title=>"Keep All", :tag=>"log_all", :value=>(log_level == "all" ? 1 : 0)),
 426                   {:margin=>0, :padding=>0}
 427                 ),
 428                 -1,
 429                 item(:text, :title=>"Net Molecular Charge:"),
 430                 item(:textfield, :width=>"80", :tag=>"nc", :value=>(get_global_settings("antechamber.nc") || "0"))
 431           )
 432           set_attr(0, :enabled=>valid_antechamber_dir(ante_dir))
 433     }
 434         #  The hash values are set in the action() method
 435 #       print "retval = #{hash ? 1 : 0}\n"
 436         return (hash ? 1 : 0)
 437   end
 438
 439   def create_ante_log_dir(name, key)
 440     log_dir = get_global_settings("antechamber.log_dir")
 441         msg = ""
 442         begin
 443           if !FileTest.directory?(log_dir)
 444             Dir.mkdir(log_dir) rescue ((msg = "Cannot create directory #{log_dir}") && raise)
 445           end
 446           n = 1
 447           while FileTest.exist?(dname = log_dir + "/#{name}_#{key}.#{n}")
 448             n += 1
 449           end
 450           Dir.mkdir(dname) rescue ((msg = "Cannot create directory #{dname}") && raise)
 451           return dname
 452         rescue
 453           error_message_box(msg + ": " + $!.to_s)
 454           return
 455         end
 456   end
 457
 458   def clean_ante_log_dir(nkeep)
 459     def rm_recursive(f)
 460           if FileTest.directory?(f)
 461             Dir.entries(f).each { |file|
 462                   next if file == "." || file == ".."
 463                   rm_recursive(f + "/" + file)
 464                 }
 465                 Dir.rmdir(f)
 466           else
 467             File.delete(f)
 468           end
 469         end
 470     log_dir = get_global_settings("antechamber.log_dir")
 471         cwd = Dir.pwd
 472         count = 0
 473         begin
 474           Dir.chdir(log_dir)
 475           #  Get subdirectories and sort by last modified time
 476           dirs = Dir.entries(log_dir).reject! { |x|
 477             !FileTest.directory?(x) || x == "." || x == ".."
 478           }.sort_by { |x|
 479             File.mtime(x)
 480           }
 481           dirs[0..-(nkeep + 1)].each { |d|
 482             rm_recursive(d)
 483                 count += 1
 484           }
 485         rescue
 486           error_message_box $!.to_s
 487           Dir.chdir(cwd)
 488           return
 489         end
 490         Dir.chdir(cwd)
 491         return count
 492   end
 493
 494   def cmd_antechamber
 495     return ambertools_dialog("antechamber")
 496   end
 497
 498   def import_ac(acfile)
 499     open(acfile, "r") { |fp|
 500           while (s = fp.gets)
 501             next if s !~ /^ATOM/
 502                 s.chomp!
 503                 idx = Integer(s[4..11]) - 1
 504                 charge = Float(s[54..63])
 505                 type = s[72..-1]
 506                 type.gsub!(/ /, "")
 507                 ap = atoms[idx]
 508                 ap.charge = charge
 509                 ap.atom_type = type
 510           end
 511         }
 512   end
 513
 514   def import_frcmod(file)
 515     self.md_arena.prepare(true)  #  Clean up existing parameters
 516         par = self.parameter
 517     open(file, "r") { |fp|
 518           wtable = Hash.new
 519           state = 0
 520           while (s = fp.gets)
 521             s.chomp!
 522                 case s
 523                 when /^MASS/
 524                   state = 1
 525                   next
 526                 when /^BOND/
 527                   state = 2
 528                   next
 529                 when /^ANGLE/
 530                   state = 3
 531                   next
 532                 when /^DIHE/
 533                   state = 4
 534                   next
 535                 when /^IMPR/
 536                   state = 5
 537                   next
 538                 when /^NONB/
 539                   state = 6
 540                   next
 541                 when ""
 542                   state = 0
 543                   next
 544                 else
 545                   case state
 546                   when 1
 547                         name, weight = s.split
 548                         wtable[name] = Float(weight)
 549                   when 2
 550                     types, k, r0, com = s.split(nil, 4)
 551                     pp = par.bonds.lookup(types, :local, :missing) || par.bonds.insert
 552                         pp.atom_types = types
 553                         pp.k = k
 554                         pp.r0 = r0
 555                         pp.comment = com
 556                   when 3
 557                     types, k, a0, com = s.split(nil, 4)
 558                         pp = par.angles.lookup(types, :local, :missing) || par.angles.insert
 559                         pp.atom_types = types
 560                         pp.k = k
 561                         pp.a0 = a0
 562                         pp.comment = com
 563                   when 4
 564                     types, n, k, phi0, period, com = s.split(nil, 6)
 565                         pp = par.dihedrals.lookup(types, :local, :missing) || par.dihedrals.insert
 566                         pp.atom_types = types
 567                         pp.mult = 1
 568                         pp.k = k
 569                         pp.phi0 = phi0
 570                         pp.period = Float(period).round
 571                         pp.comment = com
 572                   when 5
 573                     types, k, phi0, period, com = s.split(nil, 5)
 574                         pp = par.impropers.lookup(types, :local, :missing) || par.impropers.insert
 575                         pp.atom_types = types
 576                         pp.mult = 1
 577                         pp.k = k
 578                         pp.phi0 = phi0
 579                         pp.period = Float(period).round
 580                         pp.comment = com
 581                   when 6
 582                     name, r_eq, eps, com = s.split(nil, 4)
 583                         pp = par.vdws.lookup(name, :local, :missing) || par.vdws.insert
 584                         pp.atom_type = name
 585                         pp.r_eq = r_eq
 586                         pp.r_eq14 = r_eq
 587                         pp.eps = eps
 588                         pp.eps14 = eps
 589                         if wtable[name]
 590                           pp.weight = wtable[name]
 591                         end
 592                         pp.comment = com
 593                   end
 594                 end
 595           end
 596     }
 597   end
 598
 599   def count_elements
 600     elements = []
 601         each_atom { |ap|
 602           if (p = elements.assoc(ap.atomic_number))
 603             p[1] += 1
 604           else
 605             elements.push [ap.atomic_number, 1]
 606           end
 607         }
 608         elements.sort_by { |p| p[0] }
 609   end
 610
 611   def cmd_run_resp
 612     if natoms == 0
 613           error_message_box "Molecule is empty"
 614           return
 615         elsif nelpots == 0
 616           error_message_box "No ESP information is loaded"
 617           return
 618         end
 619         return unless ambertools_dialog("resp")
 620         nc = get_global_settings("antechamber.nc")
 621         ante_dir = get_global_settings("antechamber.ante_dir")
 622
 623         #  Create the temporary directory
 624         dname = create_ante_log_dir((self.path ? File.basename(self.path, ".*") : self.name), "rs")
 625         return unless dname
 626         cwd = Dir.pwd
 627         Dir.chdir(dname)
 628
 629         if true
 630           eq = search_equivalent_atoms
 631           #  search for methyl (and methylene??) hydrogens
 632           methyl = []
 633           each_atom { |ap|
 634             if ap.element == "C" && ap.connects.length == 4
 635                   hs = ap.connects.select { |n| atoms[n].element == "H" }
 636                   if hs.length == 3
 637                     mc = hs.min
 638                     hs.each { |n| methyl[n] = mc }  #  methyl hydrogens
 639                         methyl[ap.index] = -1           #  methyl carbon
 640                   end
 641                 end
 642           }
 643           for i in [1,2]
 644             open("resp.input#{i}", "w") { |fp|
 645                   fp.print " #{self.name} \n"
 646                   fp.print " &cntrl\n"
 647                   fp.print " ioutput=1, IQOPT=#{i}, nmol=1, ihfree=1, irstrnt=1, qwt=#{i*0.0005},\n"
 648                   fp.print " &end\n"
 649                   fp.print "  1.0\n"
 650                   fp.print " #{self.name} \n"
 651                   fp.printf "%4d%5d\n", nc, self.natoms
 652                   each_atom { |ap|
 653                     idx = ap.index
 654                         if i == 1
 655                           n = (methyl[idx] ? 0 : eq[idx] + 1)
 656                           n = 0 if n == idx + 1
 657                         else
 658                           n = (methyl[idx] ? methyl[idx] + 1 : -1)
 659                           n = 0 if n == idx + 1
 660                         end
 661                         fp.printf "%4d%5d\n", ap.atomic_number, n
 662                   }
 663                   fp.print "\n\n\n\n\n\n"
 664                 }
 665           end
 666         else
 667       #  Export as an antechamber format
 668           c = count_elements
 669           formula = c.map { |p| Parameter.atoms[p[0]].name + p[1].to_s }.join(" ")
 670           open("respgen_in.ac", "w") { |fp|
 671             fp.printf("CHARGE %9.2f ( %d )\n", nc, nc)
 672             fp.printf("Formula: #{formula}\n")
 673             each_atom { |ap|
 674               fp.printf("ATOM %6d  %-4s%-3s%6d%12.3f%8.3f%8.3f%10.6f%8s\n", ap.index + 1, ap.name, ap.res_name, ap.res_seq, ap.r.x, ap.r.y, ap.r.z, ap.charge, ap.atom_type)
 675             }
 676             bonds.each_with_index { |b, i|
 677               fp.printf("BOND%5d%5d%5d%5d  %5s%5s\n", i + 1, b[0] + 1, b[1] + 1, 0, atoms[b[0]].name, atoms[b[1]].name)
 678             }
 679           }
 680
 681           #  Create resp input by respgen
 682           if !system("#{ante_dir}/respgen -i respgen_in.ac -o resp.input1 -f resp1") \
 683           || !system("#{ante_dir}/respgen -i respgen_in.ac -o resp.input2 -f resp2")
 684             error_message_box("Cannot run respgen.")
 685             Dir.chdir(cwd)
 686             return
 687           end
 688         end
 689
 690         #  Create ESP file
 691         a2b = 1.8897259885   #  angstrom to bohr
 692         open("resp.esp", "w") { |fp|
 693           fp.printf("%5d%5d%5d\n", natoms, nelpots, nc)
 694           each_atom { |ap|
 695             fp.printf("                %16.7E%16.7E%16.7E\n", ap.r.x * a2b, ap.r.y * a2b, ap.r.z * a2b)
 696           }
 697           nelpots.times { |i|
 698             pos, esp = elpot(i)
 699             fp.printf("%16.7E%16.7E%16.7E%16.7E\n", esp, pos.x, pos.y, pos.z)
 700           }
 701           fp.print "\n\n"
 702         }
 703
 704     #  Run resp
 705         if !system("#{ante_dir}/resp -O -i resp.input1 -o resp.output1 -e resp.esp -t qout_stage1") \
 706         || !system("#{ante_dir}/resp -O -i resp.input2 -o resp.output2 -e resp.esp -q qout_stage1 -t qout_stage2")
 707           error_message_box("Cannot run resp.")
 708           Dir.chdir(cwd)
 709           return
 710         end
 711
 712         #  Import resp output
 713         open("punch", "r") { |fp|
 714           while (s = fp.gets)
 715             next unless s =~ /Point charges/ && s =~ /after optimization/
 716                 s = fp.gets
 717                 i = 0
 718                 while (s = fp.gets)
 719                   ary = s.split
 720                   break if ary.count == 0
 721                   if Integer(ary[1]) != atoms[i].atomic_number
 722                     error_message_box(sprintf("The atom %d has inconsistent atomic number (%d in mol, %d in resp output)", i, atoms[i].atomic_number, ary[1]))
 723                         Dir.chdir(cwd)
 724                         return
 725                   end
 726                   atoms[i].charge = Float(ary[3])
 727                   i += 1
 728                 end
 729                 break
 730           end
 731     }
 732         Dir.chdir(cwd)
 733         return true
 734   end
 735
 736   def cmd_gamess_resp
 737         if natoms == 0
 738           error_message_box "Molecule is empty"
 739           return
 740         end
 741         mol = self;
 742         RubyDialog.run("#{name}:GAMESS/RESP", nil, nil) {
 743           layout(1,
 744             item(:text, :title=>"Step 1:\nCreate GAMESS input for ESP calculation"),
 745                 [item(:button, :title=>"Create GAMESS Input...",
 746                   :action=>proc {
 747                     esp_save = get_global_settings("gamess.esp")
 748                         set_global_settings("gamess.esp", 1)
 749                         mol.cmd_create_gamess_input
 750                         set_global_settings("gamess.esp", esp_save)
 751                   }),
 752                   {:align=>:right}],
 753                 item(:line),
 754                 item(:text, :title=>"Step 2:\nImport GAMESS .dat file"),
 755                 [item(:button, :title=>"Import GAMESS dat...",
 756                   :action=>proc {
 757                     fname = RubyDialog.open_panel("Select GAMESS .dat file", nil, "*.dat")
 758                         if fname
 759                           errmsg = nil
 760                           begin
 761                             mol.loaddat(fname)
 762                                 errmsg = "Cannot find ESP results in the dat file." if mol.nelpots == 0
 763                           rescue
 764                             errmsg = "Error reading the dat file."
 765                           end
 766                           if errmsg
 767                             message_box(errmsg + "\nYou may want to try another .dat file.", "GAMESS Import Error", :ok, :warning)
 768                           else
 769                             set_attr("resp", :enabled=>true)
 770                           end
 771                         end
 772                   }),
 773                   {:align=>:right}],
 774                 item(:line),
 775                 item(:text, :title=>"Step 3:\nRun RESP for charge fitting"),
 776                 [item(:button, :title=>"Run RESP...", :tag=>"resp",
 777                   :action=>proc {
 778                     if mol.cmd_run_resp
 779                           if get_global_settings("antechamber.log_level") == "latest"
 780                             mol.clean_ante_log_dir(Integer(get_global_settings("antechamber.log_keep_number")))
 781                           end
 782                           action(0)
 783                         end
 784                   }),
 785                   {:align=>:right}],
 786                 item(:line),
 787                 [item(:button, :title=>"Close", :action=>proc { action(1) }),
 788                   {:align=>:center}]
 789           )
 790           if mol.nelpots == 0
 791             set_attr("resp", :enabled=>false)
 792           end
 793         }
 794   end
 795
 796   def cmd_edit_local_parameter_in_mainview(ptype, names, types, value, params)
 797     #  The parameters are given as space separated strings
 798         #  e.g. "1.862 200.000"
 799         case ptype
 800         when "bond"
 801           k = ["r0", "k"]
 802           pen = self.parameter.bonds
 803         when "angle"
 804           k = ["a0", "k"]
 805           pen = self.parameter.angles
 806         when "dihedral"
 807           k = ["k", "period", "phi0"]
 808           pen = self.parameter.dihedrals
 809         when "improper"
 810           k = ["k", "period", "phi0"]
 811           pen = self.parameter.impropers
 812         else
 813           return
 814         end
 815         p = params.split
 816         hash = RubyDialog.run("Edit local parameter") {
 817       layout(1,
 818             item(:text, :title=>"Edit #{ptype} parameter for #{names} (#{types})"),
 819                 item(:text, :title=>"(Current value = #{value})"),
 820             layout(4,
 821               [item(:text, :title=>"types"), {:align=>:center}],
 822               [item(:text, :title=>k[0]), {:align=>:center}],
 823               [item(:text, :title=>k[1]), {:align=>:center}],
 824               (k[2] ? [item(:text, :title=>k[2]), {:align=>:center}] : -1),
 825
 826                   item(:textfield, :width=>100, :value=>types, :tag=>"types"),
 827                   item(:textfield, :width=>100, :value=>p[0], :tag=>k[0]),
 828                   item(:textfield, :width=>100, :value=>p[1], :tag=>k[1]),
 829                   (k[2] ? item(:textfield, :width=>100, :value=>p[2], :tag=>k[2]) : -1)
 830                 )
 831           )
 832         }
 833         if hash
 834           pref = pen.lookup(hash["types"], :create, :local, :missing, :nobasetype, :nowildcard)
 835           raise "Cannot create new #{ptype} parameter" if !pref
 836           k.each { |key| pref.set_attr(key, hash[key]) }
 837           self.md_arena.prepare(true)  #  Check parameter only
 838           return 1
 839         else
 840           return 0
 841         end
 842   end
 843
 844 end