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