Scripts/uff.rb

   1 class Molby::Molecule
   2
   3 #
   4 # UFF force field parameters
   5 # Taken from RDKit source code
   6 # http://rdkit.svn.sourceforge.net/viewvc/rdkit/trunk/Code/ForceField/UFF/Params.cpp?revision=2044&view=markup
   7 # Accessed on 2012.5.8.
   8 #
   9 # Atom    an   r1     theta0  x1     D1     zeta    Z1     Vi     Uj    Xi      Hard    Radius Description
  10 UFFParams = [
  11 ["H_",    1,   0.354, 180,    2.886, 0.044, 12,     0.712, 0,     0,    4.528,  6.9452, 0.371, "H generic"],
  12 ["H_b",   1,   0.46,  83.5,   2.886, 0.044, 12,     0.712, 0,     0,    4.528,  6.9452, 0.371, "H bridging B-H-B"],
  13 ["He4+4", 2,   0.849, 90,     2.362, 0.056, 15.24,  0.098, 0,     0,    9.66,   14.92,  1.3,   "He"],
  14 ["Li",    3,   1.336, 180,    2.451, 0.025, 12,     1.026, 0,     2,    3.006,  2.386,  1.557, "Li"],
  15 ["Be3+2", 4,   1.074, 109.47, 2.745, 0.085, 12,     1.565, 0,     2,    4.877,  4.443,  1.24,  "Be2+ tetrahedral"],
  16 ["B_3",   5,   0.838, 109.47, 4.083, 0.18,  12.052, 1.755, 0,     2,    5.11,   4.75,   0.822, "B tetrahedral"],
  17 ["B_2",   5,   0.828, 120,    4.083, 0.18,  12.052, 1.755, 0,     2,    5.11,   4.75,   0.822, "B trigonal"],
  18 ["C_3",   6,   0.757, 109.47, 3.851, 0.105, 12.73,  1.912, 2.119, 2,    5.343,  5.063,  0.759, "C sp3"],
  19 ["C_R",   6,   0.729, 120,    3.851, 0.105, 12.73,  1.912, 0,     2,    5.343,  5.063,  0.759, "C aromatic"],
  20 ["C_2",   6,   0.732, 120,    3.851, 0.105, 12.73,  1.912, 0,     2,    5.343,  5.063,  0.759, "C sp2"],
  21 ["C_1",   6,   0.706, 180,    3.851, 0.105, 12.73,  1.912, 0,     2,    5.343,  5.063,  0.759, "C sp"],
  22 ["N_3",   7,   0.7,   106.7,  3.66,  0.069, 13.407, 2.544, 0.45,  2,    6.899,  5.88,   0.715, "N sp3"],
  23 ["N_R",   7,   0.699, 120,    3.66,  0.069, 13.407, 2.544, 0,     2,    6.899,  5.88,   0.715, "N aromatic"],
  24 ["N_2",   7,   0.685, 111.2,  3.66,  0.069, 13.407, 2.544, 0,     2,    6.899,  5.88,   0.715, "N sp2"],
  25 ["N_1",   7,   0.656, 180,    3.66,  0.069, 13.407, 2.544, 0,     2,    6.899,  5.88,   0.715, "N sp"],
  26 ["O_3",   8,   0.658, 104.51, 3.5,   0.06,  14.085, 2.3,   0.018, 2,    8.741,  6.682,  0.669, "O tetrahedral"],
  27 ["O_3_z", 8,   0.528, 146,    3.5,   0.06,  14.085, 2.3,   0.018, 2,    8.741,  6.682,  0.669, "O tetrahedral (zeolite)"],
  28 ["O_R",   8,   0.68,  110,    3.5,   0.06,  14.085, 2.3,   0,     2,    8.741,  6.682,  0.669, "O aromatic"],
  29 ["O_2",   8,   0.634, 120,    3.5,   0.06,  14.085, 2.3,   0,     2,    8.741,  6.682,  0.669, "O trigonal"],
  30 ["O_1",   8,   0.639, 180,    3.5,   0.06,  14.085, 2.3,   0,     2,    8.741,  6.682,  0.669, "O linear"],
  31 ["F_",    9,   0.668, 180,    3.364, 0.05,  14.762, 1.735, 0,     2,    10.874, 7.474,  0.706, "F"],
  32 ["Ne4+4", 10,  0.92,  90,     3.243, 0.042, 15.44,  0.194, 0,     2,    11.04,  10.55,  1.768, "Ne"],
  33 ["Na",    11,  1.539, 180,    2.983, 0.03,  12,     1.081, 0,     1.25, 2.843,  2.296,  2.085, "Na"],
  34 ["Mg3+2", 12,  1.421, 109.47, 3.021, 0.111, 12,     1.787, 0,     1.25, 3.951,  3.693,  1.5,   "Mg2+ tetrahedral"],
  35 ["Al3",   13,  1.244, 109.47, 4.499, 0.505, 11.278, 1.792, 0,     1.25, 4.06,   3.59,   1.201, "Al tetrahedral"],
  36 ["Si3",   14,  1.117, 109.47, 4.295, 0.402, 12.175, 2.323, 1.225, 1.25, 4.168,  3.487,  1.176, "Si tetrahedral"],
  37 ["P_3+3", 15,  1.101, 93.8,   4.147, 0.305, 13.072, 2.863, 2.4,   1.25, 5.463,  4,      1.102, "P3+ tetrahedral"],
  38 ["P_3+5", 15,  1.056, 109.47, 4.147, 0.305, 13.072, 2.863, 2.4,   1.25, 5.463,  4,      1.102, "P5+ tetrahedral"],
  39 ["P_3+q", 15,  1.056, 109.47, 4.147, 0.305, 13.072, 2.863, 2.4,   1.25, 5.463,  4,      1.102, "P organometallic"],
  40 ["S_3+2", 16,  1.064, 92.1,   4.035, 0.274, 13.969, 2.703, 0.484, 1.25, 6.928,  4.486,  1.047, "S2+ tetrahedral"],
  41 ["S_3+4", 16,  1.049, 103.2,  4.035, 0.274, 13.969, 2.703, 0.484, 1.25, 6.928,  4.486,  1.047, "S4+ tetrahedral"],
  42 ["S_3+6", 16,  1.027, 109.47, 4.035, 0.274, 13.969, 2.703, 0.484, 1.25, 6.928,  4.486,  1.047, "S6+ tetrahedral"],
  43 ["S_R",   16,  1.077, 92.2,   4.035, 0.274, 13.969, 2.703, 0,     1.25, 6.928,  4.486,  1.047, "S aromatic"],
  44 ["S_2",   16,  0.854, 120,    4.035, 0.274, 13.969, 2.703, 0,     1.25, 6.928,  4.486,  1.047, "S trigonal"],
  45 ["Cl",    17,  1.044, 180,    3.947, 0.227, 14.866, 2.348, 0,     1.25, 8.564,  4.946,  0.994, "Cl"],
  46 ["Ar4+4", 18,  1.032, 90,     3.868, 0.185, 15.763, 0.3,   0,     1.25, 9.465,  6.355,  2.108, "Ar"],
  47 ["K_",    19,  1.953, 180,    3.812, 0.035, 12,     1.165, 0,     0.7,  2.421,  1.92,   2.586, "K"],
  48 ["Ca6+2", 20,  1.761, 90,     3.399, 0.238, 12,     2.141, 0,     0.7,  3.231,  2.88,   2,     "Ca2+ octahedral"],
  49 ["Sc3+3", 21,  1.513, 109.47, 3.295, 0.019, 12,     2.592, 0,     0.7,  3.395,  3.08,   1.75,  "Sc3+ tetrahedral"],
  50 ["Ti3+4", 22,  1.412, 109.47, 3.175, 0.017, 12,     2.659, 0,     0.7,  3.47,   3.38,   1.607, "Ti4+ tetrahedral"],
  51 ["Ti6+4", 22,  1.412, 90,     3.175, 0.017, 12,     2.659, 0,     0.7,  3.47,   3.38,   1.607, "Ti6+ octahedral"],
  52 ["V_3+5", 23,  1.402, 109.47, 3.144, 0.016, 12,     2.679, 0,     0.7,  3.65,   3.41,   1.47,  "V5+ tetrahedral"],
  53 ["Cr6+3", 24,  1.345, 90,     3.023, 0.015, 12,     2.463, 0,     0.7,  3.415,  3.865,  1.402, "Cr3+ octahedral"],
  54 ["Mn6+2", 25,  1.382, 90,     2.961, 0.013, 12,     2.43,  0,     0.7,  3.325,  4.105,  1.533, "Mn2+ octahedral"],
  55 ["Fe3+2", 26,  1.27,  109.47, 2.912, 0.013, 12,     2.43,  0,     0.7,  3.76,   4.14,   1.393, "Fe2+ tetrahedral"],
  56 ["Fe6+2", 26,  1.335, 90,     2.912, 0.013, 12,     2.43,  0,     0.7,  3.76,   4.14,   1.393, "Fe2+ octahedral"],
  57 ["Co6+3", 27,  1.241, 90,     2.872, 0.014, 12,     2.43,  0,     0.7,  4.105,  4.175,  1.406, "Co3+ octahedral"],
  58 ["Ni4+2", 28,  1.164, 90,     2.834, 0.015, 12,     2.43,  0,     0.7,  4.465,  4.205,  1.398, "Ni2+ square planar"],
  59 ["Cu3+1", 29,  1.302, 109.47, 3.495, 0.005, 12,     1.756, 0,     0.7,  4.2,    4.22,   1.434, "Cu+ tetrahedral"],
  60 ["Zn3+2", 30,  1.193, 109.47, 2.763, 0.124, 12,     1.308, 0,     0.7,  5.106,  4.285,  1.4,   "Zn2+ tetrahedral"],
  61 ["Ga3+3", 31,  1.26,  109.47, 4.383, 0.415, 11,     1.821, 0,     0.7,  3.641,  3.16,   1.211, "Ga3+ tetrahedral"],
  62 ["Ge3",   32,  1.197, 109.47, 4.28,  0.379, 12,     2.789, 0.701, 0.7,  4.051,  3.438,  1.189, "Ge tetrahedral"],
  63 ["As3+3", 33,  1.211, 92.1,   4.23,  0.309, 13,     2.864, 1.5,   0.7,  5.188,  3.809,  1.204, "As3+ tetrahedral"],
  64 ["Se3+2", 34,  1.19,  90.6,   4.205, 0.291, 14,     2.764, 0.335, 0.7,  6.428,  4.131,  1.224, "Se2+ tetrahedral"],
  65 ["Br",    35,  1.192, 180,    4.189, 0.251, 15,     2.519, 0,     0.7,  7.79,   4.425,  1.141, "Br"],
  66 ["Kr4+4", 36,  1.147, 90,     4.141, 0.22,  16,     0.452, 0,     0.7,  8.505,  5.715,  2.27,  "Kr"],
  67 ["Rb",    37,  2.26,  180,    4.114, 0.04,  12,     1.592, 0,     0.2,  2.331,  1.846,  2.77,  "Rb"],
  68 ["Sr6+2", 38,  2.052, 90,     3.641, 0.235, 12,     2.449, 0,     0.2,  3.024,  2.44,   2.415, "Sr2+ octahedral"],
  69 ["Y_3+3", 39,  1.698, 109.47, 3.345, 0.072, 12,     3.257, 0,     0.2,  3.83,   2.81,   1.998, "Y3+ tetrahedral"],
  70 ["Zr3+4", 40,  1.564, 109.47, 3.124, 0.069, 12,     3.667, 0,     0.2,  3.4,    3.55,   1.758, "Zr4+ tetrahedral"],
  71 ["Nb3+5", 41,  1.473, 109.47, 3.165, 0.059, 12,     3.618, 0,     0.2,  3.55,   3.38,   1.603, "Nb5+ tetrahedral"],
  72 ["Mo6+6", 42,  1.467, 90,     3.052, 0.056, 12,     3.4,   0,     0.2,  3.465,  3.755,  1.53,  "Mo6+ octahedral"],
  73 ["Mo3+6", 42,  1.484, 109.47, 3.052, 0.056, 12,     3.4,   0,     0.2,  3.465,  3.755,  1.53,  "Mo6+ tetrahedral"],
  74 ["Tc6+5", 43,  1.322, 90,     2.998, 0.048, 12,     3.4,   0,     0.2,  3.29,   3.99,   1.5,   "Tc5+ octahedral"],
  75 ["Ru6+2", 44,  1.478, 90,     2.963, 0.056, 12,     3.4,   0,     0.2,  3.575,  4.015,  1.5,   "Ru2+ octahedral"],
  76 ["Rh6+3", 45,  1.332, 90,     2.929, 0.053, 12,     3.5,   0,     0.2,  3.975,  4.005,  1.509, "Rh3+ octahedral"],
  77 ["Pd4+2", 46,  1.338, 90,     2.899, 0.048, 12,     3.21,  0,     0.2,  4.32,   4,      1.544, "Pd2+ square planar"],
  78 ["Ag1+1", 47,  1.386, 180,    3.148, 0.036, 12,     1.956, 0,     0.2,  4.436,  3.134,  1.622, "Ag+ linear"],
  79 ["Cd3+2", 48,  1.403, 109.47, 2.848, 0.228, 12,     1.65,  0,     0.2,  5.034,  3.957,  1.6,   "Cd2+ tetrahedral"],
  80 ["In3+3", 49,  1.459, 109.47, 4.463, 0.599, 11,     2.07,  0,     0.2,  3.506,  2.896,  1.404, "In3+ tetrahedral"],
  81 ["Sn3",   50,  1.398, 109.47, 4.392, 0.567, 12,     2.961, 0.199, 0.2,  3.987,  3.124,  1.354, "Sn tetrahedral"],
  82 ["Sb3+3", 51,  1.407, 91.6,   4.42,  0.449, 13,     2.704, 1.1,   0.2,  4.899,  3.342,  1.404, "Sb3+ tetrahedral"],
  83 ["Te3+2", 52,  1.386, 90.25,  4.47,  0.398, 14,     2.882, 0.3,   0.2,  5.816,  3.526,  1.38,  "Te2+ tetrahedral"],
  84 ["I_",    53,  1.382, 180,    4.5,   0.339, 15,     2.65,  0,     0.2,  6.822,  3.762,  1.333, "I"],
  85 ["Xe4+4", 54,  1.267, 90,     4.404, 0.332, 12,     0.556, 0,     0.2,  7.595,  4.975,  2.459, "Xe"],
  86 ["Cs",    55,  2.57,  180,    4.517, 0.045, 12,     1.573, 0,     0.1,  2.183,  1.711,  2.984, "Cs"],
  87 ["Ba6+2", 56,  2.277, 90,     3.703, 0.364, 12,     2.727, 0,     0.1,  2.814,  2.396,  2.442, "Ba2+ octahedral"],
  88 ["La3+3", 57,  1.943, 109.47, 3.522, 0.017, 12,     3.3,   0,     0.1,  2.8355, 2.7415, 2.071, "La3+ tetrahedral"],
  89 ["Ce6+3", 58,  1.841, 90,     3.556, 0.013, 12,     3.3,   0,     0.1,  2.774,  2.692,  1.925, "Ce3+ octahedral"],
  90 ["Pr6+3", 59,  1.823, 90,     3.606, 0.01,  12,     3.3,   0,     0.1,  2.858,  2.564,  2.007, "Pr3+ octahedral"],
  91 ["Nd6+3", 60,  1.816, 90,     3.575, 0.01,  12,     3.3,   0,     0.1,  2.8685, 2.6205, 2.007, "Nd3+ octahedral"],
  92 ["Pm6+3", 61,  1.801, 90,     3.547, 0.009, 12,     3.3,   0,     0.1,  2.881,  2.673,  2,     "Pm3+ octahedral"],
  93 ["Sm6+3", 62,  1.78,  90,     3.52,  0.008, 12,     3.3,   0,     0.1,  2.9115, 2.7195, 1.978, "Sm3+ octahedral"],
  94 ["Eu6+3", 63,  1.771, 90,     3.493, 0.008, 12,     3.3,   0,     0.1,  2.8785, 2.7875, 2.227, "Eu3+ octahedral"],
  95 ["Gd6+3", 64,  1.735, 90,     3.368, 0.009, 12,     3.3,   0,     0.1,  3.1665, 2.9745, 1.968, "Gd3+ octahedral"],
  96 ["Tb6+3", 65,  1.732, 90,     3.451, 0.007, 12,     3.3,   0,     0.1,  3.018,  2.834,  1.954, "Tb3+ octahedral"],
  97 ["Dy6+3", 66,  1.71,  90,     3.428, 0.007, 12,     3.3,   0,     0.1,  3.0555, 2.8715, 1.934, "Dy3+ octahedral"],
  98 ["Ho6+3", 67,  1.696, 90,     3.409, 0.007, 12,     3.416, 0,     0.1,  3.127,  2.891,  1.925, "Ho3+ octahedral"],
  99 ["Er6+3", 68,  1.673, 90,     3.391, 0.007, 12,     3.3,   0,     0.1,  3.1865, 2.9145, 1.915, "Er3+ octahedral"],
 100 ["Tm6+3", 69,  1.66,  90,     3.374, 0.006, 12,     3.3,   0,     0.1,  3.2514, 2.9329, 2,     "Tm3+ octahedral"],
 101 ["Yb6+3", 70,  1.637, 90,     3.355, 0.228, 12,     2.618, 0,     0.1,  3.2889, 2.965,  2.158, "Yb3+ octahedral"],
 102 ["Lu6+3", 71,  1.671, 90,     3.64,  0.041, 12,     3.271, 0,     0.1,  2.9629, 2.4629, 1.896, "Lu3+ octahedral"],
 103 ["Hf3+4", 72,  1.611, 109.47, 3.141, 0.072, 12,     3.921, 0,     0.1,  3.7,    3.4,    1.759, "Hf4+ tetrahedral"],
 104 ["Ta3+5", 73,  1.511, 109.47, 3.17,  0.081, 12,     4.075, 0,     0.1,  5.1,    2.85,   1.605, "Ta5+ tetrahedral"],
 105 ["W_6+6", 74,  1.392, 90,     3.069, 0.067, 12,     3.7,   0,     0.1,  4.63,   3.31,   1.538, "W6+ octahedral"],
 106 ["W_3+4", 74,  1.526, 109.47, 3.069, 0.067, 12,     3.7,   0,     0.1,  4.63,   3.31,   1.538, "W4+ tetrahedral"],
 107 ["W_3+6", 74,  1.38,  109.47, 3.069, 0.067, 12,     3.7,   0,     0.1,  4.63,   3.31,   1.538, "W6+ tetrahedral"],
 108 ["Re6+5", 75,  1.372, 90,     2.954, 0.066, 12,     3.7,   0,     0.1,  3.96,   3.92,   1.6,   "Re5+ octahedral"],
 109 ["Re3+7", 75,  1.314, 109.47, 2.954, 0.066, 12,     3.7,   0,     0.1,  3.96,   3.92,   1.6,   "Re7+ tetrahedral"],
 110 ["Os6+6", 76,  1.372, 90,     3.12,  0.037, 12,     3.7,   0,     0.1,  5.14,   3.63,   1.7,   "Os6+ octahedral"],
 111 ["Ir6+3", 77,  1.371, 90,     2.84,  0.073, 12,     3.731, 0,     0.1,  5,      4,      1.866, "Ir3+ octahedral"],
 112 ["Pt4+2", 78,  1.364, 90,     2.754, 0.08,  12,     3.382, 0,     0.1,  4.79,   4.43,   1.557, "Pt2+ square planar"],
 113 ["Au4+3", 79,  1.262, 90,     3.293, 0.039, 12,     2.625, 0,     0.1,  4.894,  2.586,  1.618, "Au3+ tetrahedral"],
 114 ["Hg1+2", 80,  1.34,  180,    2.705, 0.385, 12,     1.75,  0,     0.1,  6.27,   4.16,   1.6,   "Hg2+ linear"],
 115 ["Tl3+3", 81,  1.518, 120,    4.347, 0.68,  11,     2.068, 0,     0.1,  3.2,    2.9,    1.53,  "Tl3+ tetrahedral"],
 116 ["Pb3",   82,  1.459, 109.47, 4.297, 0.663, 12,     2.846, 0.1,   0.1,  3.9,    3.53,   1.444, "Pb tetrahedral"],
 117 ["Bi3+3", 83,  1.512, 90,     4.37,  0.518, 13,     2.47,  1,     0.1,  4.69,   3.74,   1.514, "Bi3+ tetrahedral"],
 118 ["Po3+2", 84,  1.5,   90,     4.709, 0.325, 14,     2.33,  0.3,   0.1,  4.21,   4.21,   1.48,  "Po2+ tetrahedral"],
 119 ["At",    85,  1.545, 180,    4.75,  0.284, 15,     2.24,  0,     0.1,  4.75,   4.75,   1.47,  "At"],
 120 ["Rn4+4", 86,  1.42,  90,     4.765, 0.248, 16,     0.583, 0,     0.1,  5.37,   5.37,   2.2,   "Rn"],
 121 ["Fr",    87,  2.88,  180,    4.9,   0.05,  12,     1.847, 0,     0,    2,      2,      2.3,   "Fr"],
 122 ["Ra6+2", 88,  2.512, 90,     3.677, 0.404, 12,     2.92,  0,     0,    2.843,  2.434,  2.2,   "Ra2+ octahedral"],
 123 ["Ac6+3", 89,  1.983, 90,     3.478, 0.033, 12,     3.9,   0,     0,    2.835,  2.835,  2.108, "Ac3+ octahedral"],
 124 ["Th6+4", 90,  1.721, 90,     3.396, 0.026, 12,     4.202, 0,     0,    3.175,  2.905,  2.018, "Th4+ octahedral"],
 125 ["Pa6+4", 91,  1.711, 90,     3.424, 0.022, 12,     3.9,   0,     0,    2.985,  2.905,  1.8,   "Pa4+ octahedral"],
 126 ["U_6+4", 92,  1.684, 90,     3.395, 0.022, 12,     3.9,   0,     0,    3.341,  2.853,  1.713, "U4+ octahedral"],
 127 ["Np6+4", 93,  1.666, 90,     3.424, 0.019, 12,     3.9,   0,     0,    3.549,  2.717,  1.8,   "Np4+ octahedral"],
 128 ["Pu6+4", 94,  1.657, 90,     3.424, 0.016, 12,     3.9,   0,     0,    3.243,  2.819,  1.84,  "Pu4+ octahedral"],
 129 ["Am6+4", 95,  1.66,  90,     3.381, 0.014, 12,     3.9,   0,     0,    2.9895, 3.0035, 1.942, "Am4+ octahedral"],
 130 ["Cm6+3", 96,  1.801, 90,     3.326, 0.013, 12,     3.9,   0,     0,    2.8315, 3.1895, 1.9,   "Cm3+ octahedral"],
 131 ["Bk6+3", 97,  1.761, 90,     3.339, 0.013, 12,     3.9,   0,     0,    3.1935, 3.0355, 1.9,   "Bk3+ octahedral"],
 132 ["Cf6+3", 98,  1.75,  90,     3.313, 0.013, 12,     3.9,   0,     0,    3.197,  3.101,  1.9,   "Cf3+ octahedral"],
 133 ["Es6+3", 99,  1.724, 90,     3.299, 0.012, 12,     3.9,   0,     0,    3.333,  3.089,  1.9,   "Es3+ octahedral"],
 134 ["Fm6+3", 100, 1.712, 90,     3.286, 0.012, 12,     3.9,   0,     0,    3.4,    3.1,    1.9,   "Fm3+ octahedral"],
 135 ["Md6+3", 101, 1.689, 90,     3.274, 0.011, 12,     3.9,   0,     0,    3.47,   3.11,   1.9,   "Md3+ octahedral"],
 136 ["No6+3", 102, 1.679, 90,     3.248, 0.011, 12,     3.9,   0,     0,    3.475,  3.175,  1.9,   "No3+ octahedral"],
 137 ["Lw6+3", 103, 1.698, 90,     3.236, 0.011, 12,     3.9,   0,     0,    3.5,    3.2,    1.9,   "Lr3+ octahedral"]
 138 ]
 139
 140 #  Calculate UFF bond length
 141 def uff_bond_length(r1, r2, x1, x2, bond_order)
 142   bond_order_correction = -0.1332 * (r1 + r2) * log(bond_order)
 143   sq = sqrt(x1) - sqrt(x2)
 144   electronegativity_correction = r1 * r2 * (sqrt(x1) - sqrt(x2)) ** 2 / (x1 * r1 + x2 * r2)
 145   return r1 + r2 + bond_order_correction + electronegativity_correction
 146 end
 147
 148 #  UFF bond force constant
 149 def uff_bond_force(idx1, idx2, bond_order)
 150   r1 = UFFParams[idx1][2]
 151   r2 = UFFParams[idx2][2]
 152   r12 = uff_bond_length(r1, r2, UFFParams[idx1][10], UFFParams[idx2][10], bond_order)
 153   return 332.06 * UFFParams[idx1][7] * UFFParams[idx2][7] / (r12 ** 3)
 154 end
 155
 156 #  UFF angle force constant (equilibrium angle should be given in degree)
 157 def uff_angle_force(idx1, idx2, idx3, bond_order_12, bond_order_23, angle)
 158   r1 = UFFParams[idx1][2]
 159   r2 = UFFParams[idx2][2]
 160   r3 = UFFParams[idx3][2]
 161   cost = cos(angle * 3.1415927 / 180.0)
 162   r12 = uff_bond_length(r1, r2, UFFParams[idx1][10], UFFParams[idx2][10], bond_order_12)
 163   r23 = uff_bond_length(r2, r3, UFFParams[idx2][10], UFFParams[idx3][10], bond_order_23)
 164   r13 = sqrt(r12 * r12 + r23 * r23 - 2 * r12 * r23 * cost)
 165   return 332.06 * UFFParams[idx1][7] * UFFParams[idx3][7] / (r13 ** 5) * (r12 * r23 * (1.0 - cost * cost) - r13 * r13 * cost)
 166 end
 167
 168 def guess_uff_parameter_dialog(current_value, indices)
 169
 170   #  TODO: this dialog is soon to be made obsolete
 171
 172   indices = indices.split(/-/)
 173
 174   if indices.length == 2
 175     par_type = "bond"
 176   elsif indices.length == 3
 177     par_type = "angle"
 178   else
 179     message_box("UFF parameter guess of this type is not implemented.", "Guess UFF Force", "OK")
 180         return nil
 181   end
 182
 183   #  Look up the UFFParams entry
 184   uff_types = indices.map { |i|
 185     ap = atoms[i]
 186     sel = []
 187     UFFParams.each_with_index { |p, j|
 188       if p[1] == ap.atomic_number
 189         sel.push(j)
 190       end
 191     }
 192     if sel.length == 0
 193       message_box("No UFF parameter is available for atom #{ap.name}", "Guess UFF Force", "OK")
 194       return nil
 195     end
 196     sel
 197   }
 198   names = indices.map { |i| sprintf("%d:%s", atoms[i].res_seq, atoms[i].name) }
 199   types = indices.map { |i| atoms[i].atom_type }
 200   utypes = indices.map { |i| atoms[i].uff_type }
 201   names_str = names.join("-")
 202   types_str = types.join("-")
 203   recalc = lambda { |i1, i2, i3, b1, b2|
 204     i1 = uff_types[0][i1] rescue 0
 205         i2 = uff_types[1][i2] rescue 0
 206         i3 = uff_types[2][i3] rescue 0
 207     if par_type == "bond"
 208       k = uff_bond_force(i1, i2, b1)
 209     else
 210       k = uff_angle_force(i1, i2, i3, b1, b2, current_value.to_f)
 211     end
 212     sprintf("%.5f", k)
 213   }
 214   hash = Dialog.run("Guess UFF Force", "Accept", "Cancel") {
 215     action_proc = lambda { |it|
 216           t1 = value("uff_type_0").to_i rescue 0
 217           t2 = value("uff_type_1").to_i rescue 0
 218           t3 = value("uff_type_2").to_i rescue 0
 219           b1 = value("bond_order_0").to_f rescue 1.0
 220           b2 = value("bond_order_1").to_f rescue 1.0
 221           current_value = value("current_value") rescue nil
 222       set_value("guessed", recalc.call(t1, t2, t3, b1, b2)) rescue nil
 223     }
 224     type_selects = []
 225     uff_types.each_with_index { |p, i|
 226       type_selects.push(item(:text, :title=>names[i] + "->"))
 227       if p.length == 1
 228         type_selects.push(item(:text, :title=>UFFParams[p[0]][13]))
 229       else
 230             subitems = p.map { |pp| UFFParams[pp][0] }
 231                 uff_idx = subitems.index(utypes[i]) || 0
 232         type_selects.push(item(:popup, :subitems=>p.map { |pp| UFFParams[pp][13] }, :tag=>"uff_type_#{i}", :action=>action_proc, :value=>uff_idx))
 233       end
 234     }
 235     bond_orders = []
 236     (0..indices.length - 2).each { |i|
 237       bond_orders.push(item(:text, :title=>names[i] + "-" + names[i + 1]))
 238       bond_orders.push(item(:textfield, :width=>60, :value=>"1.0", :tag=>"bond_order_#{i}", :action=>action_proc))
 239     }
 240     layout(1,
 241       item(:text, :title=>"Guess UFF force parameter from atom types and bond order"),
 242       layout(2,
 243         item(:text, :title=>"UFF Atom Types:"),
 244         item(:text, :title=>"Bond order:"),
 245         layout(2, *type_selects),
 246         layout(2, *bond_orders)
 247       ),
 248           (par_type == "bond" ? nil :
 249                 layout(2,
 250                   item(:text, :title=>"Equilibrium angle = "),
 251                   item(:textfield, :width=>100, :value=>current_value, :tag=>"current_value", :action=>action_proc))
 252           ),
 253       layout(2,
 254         item(:text, :title=>"Guessed force = "),
 255         item(:textfield, :editable=>false, :width=>100, :value=>recalc.call(0, 0, 0, 1.0, 1.0), :tag=>"guessed")
 256       )
 257     )
 258   }
 259   if hash[:status] == 0
 260     3.times { |i|
 261           idx = indices[i]
 262           next unless idx
 263           ii = uff_types[i][hash["uff_type_#{i}"].to_i]
 264           if ii
 265         atoms[idx].uff_type = UFFParams[ii][0]
 266       end
 267         }
 268     return hash["guessed"], current_value
 269   else
 270     return nil
 271   end
 272 end
 273
 274 def guess_uff_parameters
 275   mol = self
 276   #  Look up the non-conventional atoms
 277   exclude = [1, 6, 7, 8, 9, 15, 16, 17, 35, 53]
 278   arena = mol.md_arena
 279   arena.prepare(true)
 280   xatoms = IntGroup[]
 281   xbonds = xangles = xdihedrals = xfragments = []
 282   h = Dialog.new("Guess UFF Parameters: #{mol.name}", nil, nil, :resizable=>true) {
 283     update_xatoms = lambda {
 284       xatoms = mol.atom_group { |ap| !exclude.member?(ap.atomic_number) }
 285       xfragments = mol.fragments(xatoms)
 286       xbonds = (0...mol.nbonds).select { |i| b = mol.bonds[i]; xatoms.include?(b[0]) || xatoms.include?(b[1]) }
 287       xangles = (0...mol.nangles).select { |i| a = mol.angles[i]; xatoms.include?(a[0]) || xatoms.include?(a[1]) || xatoms.include?(a[2]) }
 288       xdihedrals = (0...mol.ndihedrals).select { |i| d = mol.dihedrals[i]; xatoms.include?(d[0]) || xatoms.include?(d[1]) || xatoms.include?(d[2]) || xatoms.include?(d[3]) }
 289       xbonds.each { |i| xatoms.add(mol.bonds[i]) }
 290       xangles.each { |i|
 291             ang = mol.angles[i]
 292                 xatoms.add(ang)
 293                 2.times { |j|
 294                   b0 = ang[j]
 295                   b1 = ang[j + 1]
 296                   k = 0
 297                   mol.bonds.each { |b|
 298                     break if (b[0] == b0 && b[1] && b1) || (b[0] == b1 && b[1] && b0)
 299                         k += 1
 300                   }
 301                   if k < mol.nbonds && !xbonds.include?(k)
 302                     xbonds.push(k)
 303                   end
 304                 }
 305           }
 306           xbonds.sort!
 307       item_with_tag("table")[:refresh] = true
 308     }
 309     columns = {
 310       "atoms"=>[["no", 40], ["name", 60], ["type", 40], ["element", 40], ["int_charge", 40], ["uff_type", 120], ["weight", 70], ["eps", 70], ["r", 70], ["eps14", 70], ["r14", 70]],
 311       "bonds"=>[["no", 40], ["bond", 110], ["nums", 80], ["types", 90], ["par_types", 90], ["order", 80], ["k", 80], ["r0", 80], ["real_r", 80]],
 312       "angles"=>[["no", 40], ["angle", 140], ["nums", 100], ["types", 100], ["par_types", 100], ["k", 80], ["a0", 80], ["real_a", 80]],
 313       "dihedrals"=>[["no", 40], ["dihedral", 180], ["nums", 100], ["types", 110], ["par_types", 110], ["k", 50], ["period", 20], ["phi0", 50], ["real_phi", 60]],
 314       "fragments"=>[["no", 40], ["fragment", 360]]
 315     }
 316     tab = "atoms"
 317     update_selection = lambda {
 318       g = mol.selection
 319       sel = IntGroup[]
 320       case tab
 321       when "atoms"
 322         xatoms.each_with_index { |n, i|
 323           sel.add(i) if g.include?(n)
 324         }
 325       when "bonds"
 326         xbonds.each_with_index { |n, i|
 327           b = mol.bonds[n]
 328           sel.add(i) if g.include?(b[0]) && g.include?(b[1])
 329         }
 330       when "angles"
 331         xangles.each_with_index { |n, i|
 332           an = mol.angles[n]
 333           sel.add(i) if g.include?(an[0]) && g.include?(an[1]) && g.include?(an[2])
 334         }
 335       when "dihedrals"
 336         xdihedrals.each_with_index { |n, i|
 337           di = mol.dihedrals[n]
 338           sel.add(i) if g.include?(di[0]) && g.include?(di[1]) && g.include?(di[2]) && g.include?(di[3])
 339         }
 340       when "fragments"
 341         xfragments.each_with_index { |f, i|
 342           sel.add(i) if (f - g).count == 0
 343         }
 344       end
 345       it = item_with_tag("table")
 346       if sel != it[:selection]
 347         @dont_change_mol_selection = true
 348         it[:selection] = sel
 349       end
 350     }
 351     selection_changed = lambda { |it|
 352       if @dont_change_mol_selection
 353         @dont_change_mol_selection = false
 354         return
 355       end
 356       sel = it[:selection]
 357       g = IntGroup[]
 358       case tab
 359       when "atoms"
 360         sel.each { |idx|
 361           g.add(xatoms[idx])
 362         }
 363       when "bonds"
 364         sel.each { |idx|
 365           g.add(mol.bonds[xbonds[idx]])
 366         }
 367       when "angles"
 368         sel.each { |idx|
 369           g.add(mol.angles[xangles[idx]])
 370         }
 371       when "dihedrals"
 372         sel.each { |idx|
 373           g.add(mol.dihedrals[xdihedrals[idx]])
 374         }
 375       when "fragments"
 376         sel.each { |idx|
 377           g.add(xfragments[idx])
 378         }
 379       end
 380       mol.selection = g
 381     }
 382     select_tab = lambda { |tag|
 383       table = item_with_tag("table")
 384       table[:columns] = columns[tag]
 385       tab = tag
 386       table[:refresh] = true
 387       update_selection.call
 388     }
 389     tab_button_pressed = lambda { |it|
 390       ["atoms", "bonds", "angles", "dihedrals", "fragments"].each { |tag|
 391         next if tag == it[:tag]
 392         item_with_tag(tag)[:value] = 0
 393       }
 394       select_tab.call(it[:tag])
 395     }
 396     uff_popup_titles = []
 397     uff_type_for_title = Hash.new
 398     uff_title_for_type = Hash.new
 399     uff_title_for_type[""] = "-- select --"
 400     uff_popup = lambda { |an|
 401       if uff_popup_titles[an] == nil
 402             if an == 0
 403                   titles = ["(no type)"]
 404                 else
 405           titles = []
 406           Molby::Molecule::UFFParams.each { |u|
 407             if u[1] == an
 408               titles.push(u[13])
 409               uff_type_for_title[u[13]] = u[0]
 410               uff_title_for_type[u[0]] = u[13]
 411             end
 412           }
 413                 end
 414         uff_popup_titles[an] = titles
 415       end
 416       uff_popup_titles[an]
 417     }
 418     get_count = lambda { |it|
 419       case tab
 420       when "atoms"
 421         return xatoms.count
 422       when "bonds"
 423         return xbonds.count
 424       when "angles"
 425         return xangles.count
 426       when "dihedrals"
 427         return xdihedrals.count
 428       when "fragments"
 429         return xfragments.count
 430       end
 431       return 0
 432     }
 433     get_value = lambda { |it, row, col|
 434       case tab
 435       when "atoms"
 436         idx = xatoms[row]
 437         ap0 = mol.atoms[idx]
 438         case col
 439         when 0
 440           return idx.to_s
 441         when 1
 442           return "#{ap0.res_seq}:#{ap0.name}"
 443         when 2
 444           return ap0.atom_type
 445         when 3
 446           return ap0.element
 447         when 4
 448           return ap0.int_charge.to_s
 449         when 5
 450           return (ap0.atomic_number == 0 ? "(no type)" : (uff_title_for_type[ap0.uff_type] || "..."))
 451         when 6
 452           return sprintf("%.3f", ap0.weight)
 453         when 7
 454           return sprintf("%.3f", arena.vdw_par(idx).eps)
 455         when 8
 456           return sprintf("%.3f", arena.vdw_par(idx).r_eq)
 457         when 9
 458           return sprintf("%.3f", arena.vdw_par(idx).eps14)
 459         when 10
 460           return sprintf("%.3f", arena.vdw_par(idx).r_eq14)
 461         end
 462       when "bonds"
 463         idx = xbonds[row]
 464         b = mol.bonds[idx]
 465         ap0 = mol.atoms[b[0]]
 466         ap1 = mol.atoms[b[1]]
 467         case col
 468         when 0
 469           return idx.to_s
 470         when 1
 471           return "#{ap0.res_seq}:#{ap0.name}-#{ap1.res_seq}:#{ap1.name}"
 472         when 2
 473           return "#{b[0]}-#{b[1]}"
 474         when 3
 475           return "#{ap0.atom_type}-#{ap1.atom_type}"
 476         when 4
 477           atom_types = arena.bond_par(idx).atom_types
 478           return "#{atom_types[0]}-#{atom_types[1]}"
 479         when 5
 480           return (o = mol.get_bond_order(idx)) && sprintf("%.3f", o)
 481         when 6
 482           return sprintf("%.3f", arena.bond_par(idx).k)
 483         when 7
 484           return sprintf("%.3f", arena.bond_par(idx).r0)
 485         when 8
 486           return sprintf("%.3f", mol.calc_bond(b[0], b[1]))
 487         end
 488       when "angles"
 489         idx = xangles[row]
 490         an = mol.angles[idx]
 491         ap0 = mol.atoms[an[0]]
 492         ap1 = mol.atoms[an[1]]
 493         ap2 = mol.atoms[an[2]]
 494         case col
 495         when 0
 496           return idx.to_s
 497         when 1
 498           return "#{ap0.res_seq}:#{ap0.name}-#{ap1.res_seq}:#{ap1.name}-#{ap2.res_seq}:#{ap2.name}"
 499         when 2
 500           return "#{an[0]}-#{an[1]}-#{an[2]}"
 501         when 3
 502           return "#{ap0.atom_type}-#{ap1.atom_type}-#{ap2.atom_type}"
 503         when 4
 504           atom_types = arena.angle_par(idx).atom_types
 505           return "#{atom_types[0]}-#{atom_types[1]}-#{atom_types[2]}"
 506         when 5
 507           return sprintf("%.3f", arena.angle_par(idx).k)
 508         when 6
 509           return sprintf("%.2f", arena.angle_par(idx).a0)
 510         when 7
 511           return sprintf("%.2f", mol.calc_angle(an[0], an[1], an[2]))
 512         end
 513       when "dihedrals"
 514         idx = xdihedrals[row]
 515         di = mol.dihedrals[idx]
 516         ap0 = mol.atoms[di[0]]
 517         ap1 = mol.atoms[di[1]]
 518         ap2 = mol.atoms[di[2]]
 519         ap3 = mol.atoms[di[3]]
 520         case col
 521         when 0
 522           return idx.to_s
 523         when 1
 524           return "#{ap0.res_seq}:#{ap0.name}-#{ap1.res_seq}:#{ap1.name}-#{ap2.res_seq}:#{ap2.name}-#{ap2.res_seq}:#{ap2.name}"
 525         when 2
 526           return "#{di[0]}-#{di[1]}-#{di[2]}-#{di[3]}"
 527         when 3
 528           return "#{ap0.atom_type}-#{ap1.atom_type}-#{ap2.atom_type}-#{ap3.atom_type}"
 529         when 4
 530           atom_types = arena.dihedral_par(idx).atom_types
 531           return "#{atom_types[0]}-#{atom_types[1]}-#{atom_types[2]}-#{atom_types[3]}"
 532         when 5
 533           return sprintf("%.3f", arena.dihedral_par(idx).k)
 534         when 6
 535           return arena.dihedral_par(idx).period.to_s
 536         when 7
 537           return sprintf("%.2f", arena.dihedral_par(idx).phi0)
 538         when 8
 539           return sprintf("%.2f", mol.calc_dihedral(di[0], di[1], di[2], di[3]))
 540         end
 541       when "fragments"
 542         case col
 543         when 0
 544           return row.to_s
 545         when 1
 546           return xfragments[row].to_s[9..-2]  #  Remove "IntGroup[" and "]"
 547         end
 548       end
 549       return "..."
 550     }
 551     is_item_editable = lambda { |it, row, col|
 552       case tab
 553       when "atoms"
 554         case col
 555         when 1..4, 6..10
 556           return true
 557         end
 558       when "bonds"
 559         case col
 560         when 5..7
 561           return true
 562         end
 563       when "angles"
 564         case col
 565         when 5..6
 566           return true
 567         end
 568       when "dihedrals"
 569         case col
 570         when 5..7
 571           return true
 572         end
 573       end
 574       return false
 575     }
 576     modify_parameter = lambda { |mol, partype, idx, attr, val|
 577       arena = mol.md_arena
 578       case partype
 579       when "vdw"
 580         pref = arena.vdw_par(idx)
 581         pen = mol.parameter.vdws
 582       when "bond"
 583         pref = arena.bond_par(idx)
 584         pen = mol.parameter.bonds
 585       when "angle"
 586         pref = arena.angle_par(idx)
 587         pen = mol.parameter.angles
 588       when "dihedral"
 589         pref = arena.dihedral_par(idx)
 590         pen = mol.parameter.dihedrals
 591       when "improper"
 592         pref = arena.improper_par(idx)
 593         pen = mol.parameter.impropers
 594       end
 595       pref_new = pen.lookup(pref.atom_types, :create, :local, :missing, :nobasetype, :nowildcard)
 596       pref.keys.each { |k|
 597         next if k == :source || k == :index || k == :par_type
 598         pref_new.set_attr(k, pref.get_attr(k))
 599       }
 600       if attr != nil
 601         pref_new.set_attr(attr, val)
 602         arena.prepare(true)
 603       end
 604       return pref_new
 605     }
 606     set_value = lambda { |it, row, col, val|
 607       case tab
 608       when "atoms"
 609         idx = xatoms[row]
 610         ap0 = mol.atoms[idx]
 611         case col
 612         when 1
 613           val = val.to_s
 614           if val =~ /\d:/
 615             val = Regexp.last_match.post_match
 616           end
 617           ap0.name = val
 618         when 2
 619           ap0.atom_type = val
 620         when 3
 621           ap0.element = val
 622         when 4
 623           ap0.int_charge = val.to_i
 624         when 6
 625           ap0.weight = val.to_f
 626         when 7
 627           pp = modify_parameter.call(mol, "vdw", idx, :eps, val)
 628           if pp.eps14 == 0.0
 629             pp.eps14 = val.to_f
 630           end
 631         when 8
 632           pp = modify_parameter.call(mol, "vdw", idx, :r_eq, val)
 633           if pp.r_eq14 == 0.0
 634             pp.r_eq14 = val.to_f
 635           end
 636         when 9
 637           modify_parameter.call(mol, "vdw", idx, :eps14, val)
 638         when 10
 639           modify_parameter.call(mol, "vdw", idx, :r_eq14, val)
 640         end
 641       when "bonds"
 642         idx = xbonds[row]
 643         case col
 644         when 5
 645           mol.assign_bond_order(idx, val.to_f)
 646         when 6
 647           modify_parameter.call(mol, "bond", idx, :k, val)
 648         when 7
 649           modify_parameter.call(mol, "bond", idx, :r0, val)
 650         end
 651       when "angles"
 652         idx = xangles[row]
 653         case col
 654         when 5
 655           modify_parameter.call(mol, "angle", idx, :k, val)
 656         when 6
 657           modify_parameter.call(mol, "angle", idx, :a0, val)
 658         end
 659       when "dihedrals"
 660         idx = xdihedrals[row]
 661         case col
 662         when 5
 663           modify_parameter.call(mol, "dihedral", idx, :k, val)
 664         when 6
 665           modify_parameter.call(mol, "dihedral", idx, :period, val)
 666         when 7
 667           modify_parameter.call(mol, "dihedral", idx, :phi0, val)
 668         end
 669       end
 670     }
 671     has_popup_menu = lambda { |it, row, col|
 672       if tab == "atoms" && col == 5
 673         #  UFF type popup
 674         ap = mol.atoms[xatoms[row]]
 675         val = uff_popup.call(ap.atomic_number)
 676         return val
 677       else
 678         return nil
 679       end
 680     }
 681     popup_menu_selected = lambda { |it, row, col, sel|
 682       return if tab != "atoms" || col != 5
 683       ap = mol.atoms[xatoms[row]]
 684       title = uff_popup.call(ap.atomic_number)[sel]
 685       ap.uff_type = (uff_type_for_title[title] || "")
 686     }
 687     guess_uff_types = lambda { |recalc_all|
 688       xatoms.each { |idx|
 689         ap = mol.atoms[idx]
 690         next if !recalc_all && ap.uff_type != ""
 691         u = uff_popup.call(ap.atomic_number)
 692         if u.length == 1
 693           ap.uff_type = (uff_type_for_title[u[0]] || "")
 694           next
 695         end
 696         case ap.atom_type
 697         when "c3", "cx", "cy"
 698           ap.uff_type = "C_3"
 699         when "ca", "cc", "cd", "cp", "cq"
 700           ap.uff_type = "C_R"
 701         when "c2", "ce", "cf", "cu", "cv", "c"
 702           ap.uff_type = "C_2"
 703         when "c1", "cg", "ch"
 704           ap.uff_type = "C_1"
 705         when "n3", "n4", "nh"
 706           ap.uff_type = "N_3"
 707         when "nb", "nc", "nd"
 708           ap.uff_type = "N_R"
 709         when "n", "n2", "na", "ne", "nf"
 710           ap.uff_type = "N_2"
 711         when "n1"
 712           ap.uff_type = "N_1"
 713         when "oh", "os", "ow"
 714           ap.uff_type = "O_2"
 715         else
 716           ap.uff_type = ""
 717         end
 718       }
 719     }
 720     set_color = lambda { |it, row, col|
 721       @red_color ||= [1.0, 0.2, 0.2]
 722       @yellow_color ||= [1.0, 1.0, 0.6]
 723       arena = mol.md_arena
 724       case tab
 725       when "atoms"
 726         pp = arena.vdw_par(xatoms[row])
 727       when "bonds"
 728         pp = arena.bond_par(xbonds[row])
 729       when "angles"
 730         pp = arena.angle_par(xangles[row])
 731       when "dihedrals"
 732         pp = arena.dihedral_par(xdihedrals[row])
 733       end
 734       src = pp.source
 735       if src == nil
 736         return [nil, @yellow_color]
 737       elsif src == false
 738         return [nil, @red_color]
 739       else
 740         return nil
 741       end
 742     }
 743     guess_parameters_for_fragments = lambda { |*d|
 744       name = mol.name
 745       xfragments.each_with_index { |frag, i|
 746         fmol = mol.extract(frag)
 747                 mol.selection = frag
 748                 frag_str = frag.to_s[9..-2]  #  Remove "IntGroup[" and "]"
 749                 mes = "Guess MM/MD Parameters for #{mol.name}.fragment #{i}"
 750                 n = fmol.ambertools_dialog("antechamber", mes, frag_str)
 751                 break if n == 0
 752                 next if n == -1
 753         n = fmol.invoke_antechamber(false, mes)
 754         break if n == 1
 755         calc_charge = get_global_settings("antechamber.calc_charge").to_i
 756         guess_atom_types = get_global_settings("antechamber.guess_atom_types").to_i
 757         if calc_charge
 758           #  Copy partial charges
 759           frag.each_with_index { |n, i|
 760             mol.atoms[n].charge = fmol.atoms[i].charge
 761           }
 762         end
 763         if guess_atom_types
 764           #  Copy atom types and local parameters
 765           frag.each_with_index { |n, i|
 766             mol.atoms[n].atom_type = fmol.atoms[i].atom_type
 767           }
 768           [:bond, :angle, :dihedral, :improper, :vdw].each { |ptype|
 769             case ptype
 770             when :bond
 771               pen = fmol.parameter.bonds
 772             when :angle
 773               pen = fmol.parameter.angles
 774             when :dihedral
 775               pen = fmol.parameter.dihedrals
 776             when :improper
 777               pen = fmol.parameter.impropers
 778             when :vdw
 779               pen = fmol.parameter.vdws
 780             end
 781             pen.each { |pref|
 782               next if pref.source != nil
 783               pref_new = mol.parameter.lookup(ptype, pref.atom_types, :local, :missing, :create, :nowildcard, :nobasetype)
 784               if pref_new != nil
 785                 pref.keys.each { |k|
 786                   next if k == :index || k == :par_type || k == :source
 787                   pref_new.set_attr(k, pref.get_attr(k))
 788                 }
 789               end
 790             }
 791           }
 792                   guess_uff_types.call(false)
 793         end
 794       }
 795     }
 796     guess_parameters_for_metals = lambda { |*d|
 797       catch(:exit) {
 798         #  Atoms
 799         xatoms.each { |idx|
 800                   pref = arena.vdw_par(idx)
 801                   next if pref.source != false   #  Already defined
 802           ap0 = mol.atoms[idx]
 803           next if exclude.member?(ap0.atomic_number)
 804                   next if ap0.anchor_list != nil
 805           uff_type = ap0.uff_type
 806           u = UFFParams.find { |u| u[0] == uff_type }
 807           if u == nil
 808             error_message_box("The UFF type for atom #{idx} (#{ap0.name}) is not defined.")
 809             throw(:exit)
 810           end
 811           pref = mol.parameter.lookup(:vdw, ap0.index, :local, :missing, :create, :nowildcard, :nobasetype)
 812           pref.atom_type = idx
 813           pref.eps = pref.eps14 = u[5]
 814           pref.r_eq = pref.r_eq14 = u[4] * 0.5
 815                   pref.atomic_number = ap0.atomic_number
 816                   pref.weight = ap0.weight
 817         }
 818                 #  Bonds
 819                 pars = []
 820                 xbonds.each { |idx|
 821                   pref = arena.bond_par(idx)
 822                   next if pref.source != false && pref.k > 0.0   #  Already defined
 823                   b = mol.bonds[idx]
 824                   is = []
 825                   aps = [mol.atoms[b[0]], mol.atoms[b[1]]]
 826                   2.times { |i|
 827                     if aps[i].anchor_list != nil
 828                           is[i] = -1
 829                           next
 830                         end
 831                     uff_type = aps[i].uff_type
 832                     UFFParams.each_with_index { |u, j|
 833                           if u[0] == uff_type
 834                             is[i] = j
 835                                 break
 836                           end
 837                         }
 838                     if is[i] == nil
 839                           error_message_box("The UFF type for atom #{b[i]} (#{aps[i].name}) is not defined.")
 840                           throw(:exit)
 841                         end
 842                   }
 843                   bo = mol.get_bond_order(idx)
 844                   if bo == nil || bo == 0.0
 845                     bo = 1.0
 846                   end
 847                   if pref.source != false && pref.r0 > 0.0
 848                     len = pref.r0
 849               else
 850                     len = mol.calc_bond(b[0], b[1])
 851               end
 852                   if is[0] == -1 && is[1] == -1
 853                     #  Bond between anchors: no force
 854                         force = 0.0
 855               elsif is[0] == -1 || is[1] == -1
 856                     i = (is[0] == -1 ? 1 : 0)
 857                         case aps[i].atomic_number
 858                         when 0..23
 859                           force = 135.0
 860                         when 24
 861                           force = 150.0
 862                         when 25..27
 863                           force = 205.0
 864                         when 28..36
 865                           force = 140.0
 866                         when 37..54
 867                           force = 205.0
 868                         else
 869                           force = 260.0
 870                         end
 871                   else
 872                     force = mol.uff_bond_force(is[0], is[1], bo)
 873                   end
 874                   pars[idx] = [b, force, len]
 875                 }
 876                 pars.each { |pp|
 877                   next if pp == nil
 878                   pref = mol.parameter.lookup(:bond, pp[0], :local, :missing, :create, :nowildcard, :nobasetype)
 879                   if pref.source == false
 880                     pref.atom_types = pp[0]
 881                   end
 882                   pref.k = pp[1]
 883                   pref.r0 = pp[2]
 884             }
 885                 #  Angles
 886                 pars.clear
 887                 xangles.each { |idx|
 888                   pref = arena.angle_par(idx)
 889                   next if pref.source != false && pref.k > 0.0   #  Already defined
 890                   a = mol.angles[idx]
 891                   is = []
 892                   aps = [mol.atoms[a[0]], mol.atoms[a[1]], mol.atoms[a[2]]]
 893                   3.times { |i|
 894                     if aps[i].anchor_list != nil
 895                           is[i] = -1
 896                           next
 897                         end
 898                     uff_type = aps[i].uff_type
 899                     UFFParams.each_with_index { |u, j|
 900                           if u[0] == uff_type
 901                             is[i] = j
 902                                 break
 903                           end
 904                         }
 905                     if is[i] == nil
 906                           error_message_box("The UFF type for atom #{a[i]} (#{aps[i].name}) is not defined.")
 907                           throw(:exit)
 908                         end
 909                   }
 910                   bos = [1.0, 1.0]
 911                   2.times { |i|
 912                     #  Look up the bond and get the bond order
 913                         mol.bonds.each_with_index { |b, j|
 914                           if (a[i] == b[0] && a[i + 1] == b[1]) || (a[i] == b[1] && a[i + 1] == b[0])
 915                             bos[i] = mol.get_bond_order(j)
 916                                 if bos[i] == nil || bos[i] == 0.0
 917                                   bos[i] = 1.0
 918                                 end
 919                           end
 920                         }
 921                   }
 922                   if pref.source != false && pref.a0 > 0.0
 923                     ang = pref.a0
 924                   else
 925                     ang = mol.calc_angle(a[0], a[1], a[2])
 926                   end
 927                   met = [aps[0].atomic_number, aps[1].atomic_number, aps[2].atomic_number].max
 928                   if is[1] == -1 && is[0] != -1 && is[2] != 0
 929                     #  Metal-Dummy-C angle
 930                         case met
 931                         when 0..23
 932                           force = 85.0
 933                         when 24
 934                           force = 93.0
 935                         when 25..27
 936                           force = 100.0
 937                         when 28..36
 938                           force = 28.0
 939                         when 37..54
 940                           force = 125.0
 941                         else
 942                           force = 130.0
 943                         end
 944                   elsif is[1] != -1 && is[0] == -1 && is[2] == -1
 945                     #  Dummy-Metal-Dummy angle
 946                         case met
 947                         when 0..23
 948                           force = 50.0
 949                         when 24
 950                           force = 40.0
 951                         when 25..27
 952                           force = 40.0
 953                         when 28..36
 954                           force = 40.0
 955                         when 37..54
 956                           force = 52.0
 957                         else
 958                           force = 56.0
 959                         end
 960                   elsif is[0] >= 0 && is[1] >= 0 && is[2] >= 0
 961                     force = mol.uff_angle_force(is[0], is[1], is[2], bos[0], bos[1], ang)
 962                   else
 963                     force = 0.0
 964                   end
 965                   pars[idx] = [a, force, ang]
 966                 }
 967                 pars.each { |pp|
 968                   next if pp == nil
 969                   pref = mol.parameter.lookup(:angle, pp[0], :local, :missing, :create, :nowildcard, :nobasetype)
 970                   if pref.source == false
 971                     pref.atom_types = pp[0]
 972                   end
 973                   pref.k = pp[1]
 974                   pref.a0 = pp[2]
 975             }
 976                 xdihedrals.each { |idx|
 977                   pref = arena.dihedral_par(idx)
 978                   next if pref.source != false   #  Already defined
 979                   d = mol.dihedrals[idx]
 980                   is = []
 981                   aps = [mol.atoms[d[0]], mol.atoms[d[1]], mol.atoms[d[2]], mol.atoms[d[3]]]
 982                   4.times { |i|
 983                     if aps[i].anchor_list != nil
 984                           is[i] = -1
 985                           next
 986                         end
 987                     uff_type = aps[i].uff_type
 988                     UFFParams.each_with_index { |u, j|
 989                           if u[0] == uff_type
 990                             is[i] = j
 991                                 break
 992                           end
 993                         }
 994                   }
 995                   if is[1] == -1 && is[2] == -1 && is[0] != -1 && is[3] != -1
 996                     #  X-##-##-X dihedral
 997                         met = (aps[1].connects & aps[2].connects)[0]
 998                         if met != nil
 999                           case mol.atoms[met].atomic_number
1000                           when 0..36
1001                             force = 0.36
1002                       when 37..54
1003                             force = 3.4
1004                           else
1005                             force = 3.4
1006                           end
1007                           pref = mol.parameter.lookup(:dihedral, ["X", d[1], d[2], "X"], :local, :missing, :create)
1008                           pref.atom_types = ["X", d[1], d[2], "X"]
1009                           pref.period = 5
1010                           pref.phi0 = 180.0
1011                           pref.k = force
1012                           arena.prepare(true)
1013                         end
1014                   else
1015                     pref = mol.parameter.lookup(:dihedral, d, :local, :missing, :create, :nowildcard, :nobasetype)
1016                         pref.atom_types = d
1017                         pref.period = 2
1018                         pref.phi0 = 180.0
1019                         pref.k = 0.0
1020                   end
1021             }
1022       }
1023           arena.prepare(true)
1024           item_with_tag("table")[:refresh] = true
1025     }
1026     layout(1,
1027 #      layout(2,
1028 #        item(:text, :title=>"Total charge: "),
1029 #        item(:textfield, :width=>"80", :tag=>"total_charge")),
1030       layout(5,
1031         item(:togglebutton, :width=>80, :height=>24, :title=>"Atoms", :tag=>"atoms",
1032           :value=>1,
1033           :action=>tab_button_pressed),
1034         item(:togglebutton, :width=>80, :height=>24, :title=>"Bonds", :tag=>"bonds",
1035           :action=>tab_button_pressed),
1036         item(:togglebutton, :width=>80, :height=>24, :title=>"Angles", :tag=>"angles",
1037           :action=>tab_button_pressed),
1038         item(:togglebutton, :width=>80, :height=>24, :title=>"Dihedrals", :tag=>"dihedrals",
1039           :action=>tab_button_pressed),
1040         item(:togglebutton, :width=>80, :height=>24, :title=>"Fragments", :tag=>"fragments",
1041           :action=>tab_button_pressed),
1042         :padding=>0, :margin=>0),
1043       item(:table, :width=>640, :height=>240, :flex=>[0,0,0,0,1,1], :tag=>"table",
1044         :columns=>columns["atoms"],
1045         :on_count=>get_count,
1046         :on_get_value=>get_value,
1047         :on_set_value=>set_value,
1048         :is_item_editable=>is_item_editable,
1049         :on_selection_changed=>selection_changed,
1050         :has_popup_menu=>has_popup_menu,
1051         :on_popup_menu_selected=>popup_menu_selected,
1052         :on_set_color=>set_color),
1053       item(:button, :title=>"Run Antechamber for Non-Metal Fragments",
1054         :action=>guess_parameters_for_fragments, :flex=>[1,1,1,0,0,0], :align=>:center),
1055       item(:button, :title=>"Guess UFF Parameters for Bonds and Angles Including Metal Atoms",
1056         :action=>guess_parameters_for_metals, :flex=>[1,1,1,0,0,0], :align=>:center),
1057       item(:button, :title=>"Close", :action=>lambda { |it| hide }, :flex=>[1,1,1,0,0,0], :align=>:center),
1058       :flex=>[0,0,0,0,1,1]
1059     )
1060     size = self.size
1061     set_min_size(size)
1062     set_size(size[0] + 100, size[1] + 50);
1063     listen(mol, "documentModified", lambda { |d| update_xatoms.call; update_selection.call })
1064     listen(mol, "documentWillClose", lambda { |d| hide })
1065     update_xatoms.call
1066     guess_uff_types.call(true)
1067     show
1068   }
1069 end
1070
1071 end
1072