MolLib/Ruby_bind/ruby_md.c

   1 /*
   2  *  ruby_md.c
   3  *  Ruby binding
   4  *
   5  *  Created by Toshi Nagata on 09/08/11.
   6  *  Copyright 2007-2009 Toshi Nagata. All rights reserved.
   7  *
   8  This program is free software; you can redistribute it and/or modify
   9  it under the terms of the GNU General Public License as published by
  10  the Free Software Foundation version 2 of the License.
  11
  12  This program is distributed in the hope that it will be useful,
  13  but WITHOUT ANY WARRANTY; without even the implied warranty of
  14  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15  GNU General Public License for more details.
  16 */
  17
  18 #include "Molby.h"
  19 #include "../MD/MDCore.h"
  20 #include "../MD/MDSurface.h"
  21 #include "../MD/MDPressure.h"
  22
  23 //#include "env.h"  /*  For ruby_frame  */
  24
  25 #include <errno.h>
  26 #include <string.h>
  27 #include <ctype.h>
  28 #include <limits.h>
  29
  30 #pragma mark ====== Global Values ======
  31
  32 VALUE rb_cMDArena;
  33
  34 #pragma mark ====== MDArena class ======
  35
  36 struct MDArena *
  37 MDArenaFromValue(VALUE val)
  38 {
  39         if (rb_obj_is_kind_of(val, rb_cMDArena)) {
  40                 MDArena *arena;
  41                 Data_Get_Struct(val, MDArena, arena);
  42                 return arena;
  43         } else return NULL;
  44 }
  45
  46 VALUE
  47 ValueFromMDArena(struct MDArena *arena)
  48 {
  49         if (arena != NULL)
  50                 md_arena_retain(arena);
  51         return Data_Wrap_Struct(rb_cMDArena, 0, (void (*)(void *))md_arena_release, arena);
  52 }
  53
  54 static void
  55 s_MDArena_panic_func(MDArena *arena, const char *msg)
  56 {
  57         rb_raise(rb_eMolbyError, "MD Error: %s", msg);
  58 }
  59
  60 static VALUE
  61 s_MDArena_Run_or_minimize(VALUE self, VALUE arg, int minimize)
  62 {
  63         MDArena *arena;
  64         int nsteps, retval, start_step;
  65         Data_Get_Struct(self, MDArena, arena);
  66         nsteps = NUM2INT(rb_Integer(arg));
  67         if (arena->is_running)
  68                 rb_raise(rb_eMolbyError, "the simulation is already running. You cannot do simulation recursively.");
  69         if (nsteps < 0)
  70                 rb_raise(rb_eMolbyError, "the number of steps should be non-negative integer");
  71
  72         if (arena->is_initialized < 2 || arena->xmol->needsMDRebuild) {
  73                 const char *msg = md_prepare(arena, 0);
  74                 if (msg != NULL)
  75                         rb_raise(rb_eMolbyError, "cannot initialize for MD: %s", msg);
  76         }
  77         arena->end_step = arena->start_step + nsteps;
  78         md_log(arena, "%s for %d steps\n", (minimize ? "Minimizing" : "Running"), nsteps);
  79         start_step = arena->start_step;
  80
  81         /*  Create a new frame with current coordinates  */
  82 /*      if (nsteps > 0 && MoleculeGetNumberOfFrames(arena->xmol) == 0) {
  83                 ig = IntGroupNewWithPoints(0, 1, -1);
  84                 MolActionCreateAndPerform(arena->xmol, gMolActionInsertFrames, ig, 0, NULL);
  85                 IntGroupRelease(ig);
  86         } */
  87
  88         /*  Run simulation  */
  89         retval = md_main(arena, minimize);
  90
  91         if (retval == 0 || arena->step > start_step) {
  92                 int i, natoms = arena->xmol->natoms;
  93                 Atom *ap;
  94                 Vector *vp = (Vector *)malloc(sizeof(Vector) * (natoms > 4 ? natoms : 4));
  95                 IntGroup *ig = IntGroupNewWithPoints(0, natoms, -1);
  96                 int copy_cell = 0;
  97                 if (arena->pressure != NULL && arena->pressure->disabled == 0)
  98                         copy_cell = 1;
  99 #if MINIMIZE_CELL
 100                 if (arena->minimize_cell && minimize)
 101                         copy_cell = 1;
 102 #endif
 103                 if (arena->step > start_step) {
 104                         /*  Copy coordinates and velocities from mol to xmol */
 105                         for (i = 0, ap = arena->mol->atoms; i < natoms; i++, ap = ATOM_NEXT(ap))
 106                                 vp[i] = ap->r;
 107                         MolActionCreateAndPerform(arena->xmol, gMolActionSetAtomPositions, ig, natoms, vp);
 108                         for (i = 0, ap = arena->mol->atoms; i < natoms; i++, ap = ATOM_NEXT(ap))
 109                                 vp[i] = ap->v;
 110                         MolActionCreateAndPerform(arena->xmol, gMolActionSetAtomVelocities, ig, natoms, vp);
 111                         if (copy_cell && arena->mol->cell != NULL) {
 112                                 vp[0] = arena->mol->cell->axes[0];
 113                                 vp[1] = arena->mol->cell->axes[1];
 114                                 vp[2] = arena->mol->cell->axes[2];
 115                                 vp[3] = arena->mol->cell->origin;
 116                                 MolActionCreateAndPerform(arena->xmol, gMolActionSetBox, vp, vp + 1, vp + 2, vp + 3, -1, 0);
 117                         }
 118                 }
 119                 /*  Copy forces (this is valid even for "zero-step" run)  */
 120                 for (i = 0, ap = arena->mol->atoms; i < natoms; i++, ap = ATOM_NEXT(ap))
 121                         vp[i] = ap->f;
 122                 MolActionCreateAndPerform(arena->xmol, gMolActionSetAtomForces, ig, natoms, vp);
 123                 free(vp);
 124                 IntGroupRelease(ig);
 125         }
 126
 127         if (retval != 0)
 128                 rb_raise(rb_eMolbyError, "MD Error: %s", arena->errmsg);
 129
 130         if (minimize && arena->minimize_complete && rb_block_given_p())
 131                 rb_yield(self);
 132
 133         return self;
 134 }
 135
 136 /*
 137  *  call-seq:
 138  *     run(n)       -> self
 139  *
 140  *  Run the simulation for n steps.
 141  */
 142 static VALUE
 143 s_MDArena_Run(VALUE self, VALUE arg)
 144 {
 145         return s_MDArena_Run_or_minimize(self, arg, 0);
 146 }
 147
 148 /*
 149  *  call-seq:
 150  *     minimize(n)       -> self
 151  *     minimize(n) { ... } -> self
 152  *
 153  *  Minimize the energy for n steps. If a block is given, it is executed when minimization is complete.
 154  */
 155 static VALUE
 156 s_MDArena_Minimize(VALUE self, VALUE arg)
 157 {
 158         return s_MDArena_Run_or_minimize(self, arg, 1);
 159 }
 160
 161 /*
 162  *  call-seq:
 163  *     prepare(check_only = false)         -> self or nil
 164  *
 165  *  Prepare for the MD calculation; refresh the internal information even if initialization is
 166  *  already done.
 167  *  If check_only is true, then only the parameter checking is done.
 168  *  Returns self when initialization is complete; returns nil if some MM parameters are missing;
 169  *  otherwise, an exception is raised.
 170  */
 171 static VALUE
 172 s_MDArena_Prepare(int argc, VALUE *argv, VALUE self)
 173 {
 174         MDArena *arena;
 175         Int check_only = 0, status;
 176         char *msg;
 177         VALUE fval;
 178         Data_Get_Struct(self, MDArena, arena);
 179         rb_scan_args(argc, argv, "01", &fval);
 180         if (RTEST(fval))
 181                 check_only = 1;
 182         status = MoleculePrepareMDArena(arena->xmol, check_only, &msg);
 183         if (status < 0) {
 184                 /*  Exception object is created first to have a chance to do free(msg)  */
 185                 VALUE exval = rb_exc_new2(rb_eMolbyError, msg);
 186                 free(msg);
 187                 rb_exc_raise(exval);
 188         } else if (status > 0) {
 189                 free(msg);
 190                 return Qnil;
 191         } else return self;
 192 }
 193
 194 /*
 195  *  call-seq:
 196  *     energies -> [total, bond, angle, dihedral, improper, vdw, electrostatic, auxiliary, surface, kinetic, net [, ewald]
 197  *
 198  *  Get the current energies.
 199  */
 200 static VALUE
 201 s_MDArena_Energies(VALUE self)
 202 {
 203         MDArena *arena;
 204         VALUE val;
 205         int i;
 206         static Int s_indices[] = {kBondIndex, kAngleIndex, kDihedralIndex, kImproperIndex, kVDWIndex, kElectrostaticIndex, kAuxiliaryIndex, kSurfaceIndex, kKineticIndex };
 207         Data_Get_Struct(self, MDArena, arena);
 208         val = rb_ary_new();
 209         rb_ary_push(val, rb_float_new(arena->total_energy * INTERNAL2KCAL));
 210         for (i = 0; i < sizeof(s_indices) / sizeof(s_indices[0]); i++)
 211                 rb_ary_push(val, rb_float_new(arena->energies[s_indices[i]] * INTERNAL2KCAL));
 212         rb_ary_push(val, rb_float_new((arena->energies[kKineticIndex] + arena->total_energy) * INTERNAL2KCAL));
 213         if (arena->use_ewald)
 214                 rb_ary_push(val, rb_float_new((arena->energies[kESCorrectionIndex] + arena->energies[kPMEIndex]) * INTERNAL2KCAL));
 215         return val;
 216 }
 217
 218 static VALUE s_LogFileSym, s_CoordFileSym, s_VelFileSym, s_ForceFileSym,
 219 s_DebugFileSym, s_DebugOutputLevelSym, s_StepSym, s_CoordOutputFreqSym,
 220 s_EnergyOutputFreqSym, s_CoordFrameSym, s_TimestepSym, s_CutoffSym,
 221 s_ElectroCutoffSym, s_PairlistDistanceSym, s_SwitchDistanceSym, s_TemperatureSym, s_TransientTempSym,
 222 s_AverageTempSym, s_AndersenFreqSym, s_AndersenCouplingSym, s_RandomSeedSym,
 223 s_DielectricSym, s_GradientConvergenceSym, s_CoordinateConvergenceSym, s_UseXplorShiftSym,
 224 s_Scale14VdwSym, s_Scale14ElectSym, s_RelocateCenterSym,
 225 s_SurfaceProbeRadiusSym, s_SurfaceTensionSym, s_SurfacePotentialFreqSym, s_UseGraphiteSym,
 226 s_AlchemicalLambdaSym, s_AlchemicalDeltaLambdaSym, s_AlchemicalEnergySym, s_MinimizeCellSym,
 227 s_UseEwaldSym, s_EwaldBetaSym, s_EwaldGridSym, s_EwaldFreqSym, s_EwaldOrderSym;
 228
 229 struct s_MDArenaAttrDef {
 230         char *name;
 231         VALUE *symref;  /*  Address of s_LogFileSym etc. */
 232         ID id;                  /*  Ruby ID of the symbol; will be set within Init_MolbyMDTypes()  */
 233         ID sid;         /*  Ruby ID of the symbol plus '='; will be set within Init_MolbyMDTypes()  */
 234         char type;      /*  s: string (const char *), i: Int, f: Double, e: Double in energy dimension (unit conversion is necessary). */
 235                                         /*  Uppercase: read-only.  */
 236         int  offset;    /*  Offset in the MDArena structure.  */
 237 };
 238 static struct s_MDArenaAttrDef s_MDArenaAttrDefTable[] = {
 239         {"log_file",          &s_LogFileSym,          0, 0, 's', offsetof(MDArena, log_result_name)},
 240         {"coord_file",        &s_CoordFileSym,        0, 0, 's', offsetof(MDArena, coord_result_name)},
 241         {"vel_file",          &s_VelFileSym,          0, 0, 's', offsetof(MDArena, vel_result_name)},
 242         {"force_file",        &s_ForceFileSym,        0, 0, 's', offsetof(MDArena, force_result_name)},
 243         {"debug_file",        &s_DebugFileSym,        0, 0, 's', offsetof(MDArena, debug_result_name)},
 244         {"debug_output_level", &s_DebugOutputLevelSym, 0, 0, 'i', offsetof(MDArena, debug_output_level)},
 245         {"step",              &s_StepSym,             0, 0, 'I', offsetof(MDArena, step)},
 246         {"coord_output_freq", &s_CoordOutputFreqSym,  0, 0, 'i', offsetof(MDArena, coord_output_freq)},
 247         {"energy_output_freq", &s_EnergyOutputFreqSym, 0, 0, 'i', offsetof(MDArena, energy_output_freq)},
 248         {"coord_frame",       &s_CoordFrameSym,       0, 0, 'I', offsetof(MDArena, coord_result_frame)},
 249         {"timestep",          &s_TimestepSym,         0, 0, 'f', offsetof(MDArena, timestep)},
 250         {"cutoff",            &s_CutoffSym,           0, 0, 'f', offsetof(MDArena, cutoff)},
 251         {"electro_cutoff",    &s_ElectroCutoffSym,    0, 0, 'f', offsetof(MDArena, electro_cutoff)},
 252         {"pairlist_distance", &s_PairlistDistanceSym, 0, 0, 'f', offsetof(MDArena, pairlist_distance)},
 253         {"switch_distance",   &s_SwitchDistanceSym,   0, 0, 'f', offsetof(MDArena, switch_distance)},
 254         {"temperature",       &s_TemperatureSym,      0, 0, 'f', offsetof(MDArena, temperature)},
 255         {"transient_temperature", &s_TransientTempSym, 0, 0, 'F', offsetof(MDArena, transient_temperature)},
 256         {"average_temperature", &s_AverageTempSym,    0, 0, 'F', offsetof(MDArena, average_temperature)},
 257         {"andersen_freq",     &s_AndersenFreqSym,     0, 0, 'i', offsetof(MDArena, andersen_thermo_freq)},
 258         {"andersen_coupling", &s_AndersenCouplingSym, 0, 0, 'f', offsetof(MDArena, andersen_thermo_coupling)},
 259         {"random_seed",       &s_RandomSeedSym,       0, 0, 'i', offsetof(MDArena, random_seed)},
 260         {"dielectric",        &s_DielectricSym,       0, 0, 'f', offsetof(MDArena, dielectric)},
 261         {"gradient_convergence", &s_GradientConvergenceSym, 0, 0, 'f', offsetof(MDArena, gradient_convergence)},
 262         {"coordinate_convergence", &s_CoordinateConvergenceSym, 0, 0, 'f', offsetof(MDArena, coordinate_convergence)},
 263         {"use_xplor_shift",   &s_UseXplorShiftSym,    0, 0, 'i', offsetof(MDArena, use_xplor_shift)},
 264         {"scale14_vdw",       &s_Scale14VdwSym,       0, 0, 'f', offsetof(MDArena, scale14_vdw)},
 265         {"scale14_elect",     &s_Scale14ElectSym,     0, 0, 'f', offsetof(MDArena, scale14_elect)},
 266         {"relocate_center",   &s_RelocateCenterSym,   0, 0, 'i', offsetof(MDArena, relocate_center)},
 267         {"surface_probe_radius", &s_SurfaceProbeRadiusSym, 0, 0, 'f', offsetof(MDArena, probe_radius)},
 268         {"surface_tension",   &s_SurfaceTensionSym,   0, 0, 'f', offsetof(MDArena, surface_tension)},
 269         {"surface_potential_freq", &s_SurfacePotentialFreqSym, 0, 0, 'i', offsetof(MDArena, surface_potential_freq)},
 270         {"use_graphite",      &s_UseGraphiteSym,      0, 0, 'i', offsetof(MDArena, use_graphite)},
 271         {"alchemical_lambda", &s_AlchemicalLambdaSym, 0, 0, 'f', offsetof(MDArena, alchem_lambda)},
 272         {"alchemical_delta_lambda", &s_AlchemicalDeltaLambdaSym, 0, 0, 'f', offsetof(MDArena, alchem_dlambda)},
 273         {"alchemical_energy", &s_AlchemicalEnergySym, 0, 0, 'E', offsetof(MDArena, alchem_energy)},
 274         {"minimize_cell",     &s_MinimizeCellSym,     0, 0, 'b', offsetof(MDArena, minimize_cell)},
 275         {"use_ewald",         &s_UseEwaldSym,         0, 0, 'i', offsetof(MDArena, use_ewald)},
 276         {"ewald_beta",        &s_EwaldBetaSym,        0, 0, 'f', offsetof(MDArena, ewald_beta)},
 277         {"ewald_grid",        &s_EwaldGridSym,        0, 0, 0,   offsetof(MDArena, ewald_grid_x)},
 278         {"ewald_freq",        &s_EwaldFreqSym,        0, 0, 'i', offsetof(MDArena, ewald_freq)},
 279         {"ewald_order",       &s_EwaldOrderSym,       0, 0, 'i', offsetof(MDArena, ewald_order)},
 280         {NULL} /* Sentinel */
 281 };
 282
 283 static VALUE s_PresFreqSym, s_PresCouplingSym, s_PresSym, s_PresCellFlexSym, s_PresFluctCellOriginSym, s_PresFluctCellOrientSym;
 284 static struct s_MDArenaAttrDef s_MDPressureAttrDefTable[] = {
 285         {"pressure_freq",     &s_PresFreqSym,         0, 0, 'i', offsetof(MDPressureArena, freq)},
 286         {"pressure_coupling", &s_PresCouplingSym,     0, 0, 'f', offsetof(MDPressureArena, coupling)},
 287         {"pressure",          &s_PresSym,             0, 0, 'X', offsetof(MDPressureArena, apply)},
 288         {"pressure_cell_flexibility", &s_PresCellFlexSym, 0, 0, 'Y', offsetof(MDPressureArena, cell_flexibility)},
 289         {"pressure_fluctuate_cell_origin", &s_PresFluctCellOriginSym, 0, 0, 'f', offsetof(MDPressureArena, fluctuate_cell_origin)},
 290         {"pressure_fluctuate_cell_orientation", &s_PresFluctCellOrientSym, 0, 0, 'f', offsetof(MDPressureArena, fluctuate_cell_orientation)},
 291         {NULL} /* Sentinel */
 292 };
 293
 294 /*
 295  *  call-seq:
 296  *     self[attr]
 297  *
 298  *  Get the attribute value. The attribute values also can be accessed by self.attribute_name.
 299  */
 300 static VALUE
 301 s_MDArena_Get(VALUE self, VALUE attr)
 302 {
 303         MDArena *arena;
 304         int i;
 305         struct s_MDArenaAttrDef *dp;
 306         ID aid = rb_to_id(attr);
 307         Data_Get_Struct(self, MDArena, arena);
 308         if (aid == SYM2ID(s_EwaldGridSym)) {
 309                 /*  Array of three grid values  */
 310                 return rb_ary_new3(3, INT2NUM(arena->ewald_grid_x), INT2NUM(arena->ewald_grid_y), INT2NUM(arena->ewald_grid_z));
 311         }
 312         for (i = 0, dp = s_MDArenaAttrDefTable; dp->name != NULL; i++, dp++) {
 313                 if (dp->id == aid) {
 314                         char *p = (char *)arena + dp->offset;
 315                         switch (dp->type) {
 316                                 case 's':
 317                                 case 'S': {
 318                                         const char *cp = *((const char **)p);
 319                                         if (cp == NULL)
 320                                                 return Qnil;
 321                                         else
 322                                                 return Ruby_NewFileStringValue(cp);
 323                                 }
 324                                 case 'b':
 325                                 case 'B':
 326                                         return INT2NUM((Int)(*((Byte *)p)));
 327                                 case 'i':
 328                                 case 'I':
 329                                         return INT2NUM(*((Int *)p));
 330                                 case 'f':
 331                                 case 'F':
 332                                         return rb_float_new(*((Double *)p));
 333                                 case 'e':
 334                                 case 'E':
 335                                         return rb_float_new(*((Double *)p) * INTERNAL2KCAL);
 336                                 default:
 337                                         rb_raise(rb_eMolbyError, "Internal inconsistency: unknown type field");
 338                         }
 339                 }
 340         }
 341         for (i = 0, dp = s_MDPressureAttrDefTable; dp->name != NULL; i++, dp++) {
 342                 if (dp->id == aid) {
 343                         char *pp;
 344                         MDPressureArena *pres = arena->pressure;
 345                         if (pres == NULL)
 346                                 return Qnil;
 347                         pp = (char *)pres + dp->offset;
 348                         switch (dp->type) {
 349                                 case 'i':
 350                                 case 'I':
 351                                         return INT2NUM(*((Int *)pp));
 352                                 case 'b':
 353                                 case 'B':
 354                                         return INT2NUM((Int)(*((Byte *)pp)));
 355                                 case 'f':
 356                                 case 'F':
 357                                         return rb_float_new(*((Double *)pp));
 358                                 case 'e':
 359                                 case 'E':
 360                                         return rb_float_new(*((Double *)pp) * INTERNAL2KCAL);
 361                                 case 'X':
 362                                         /*  Isotropic pressure only  */
 363                                         return rb_ary_new3(3, rb_float_new(pres->apply[0]), rb_float_new(pres->apply[4]), rb_float_new(pres->apply[8]));
 364                                 case 'Y': {
 365                                         VALUE aval = rb_ary_new();
 366                                         int j;
 367                                         for (j = 0; j < 8; j++)
 368                                                 rb_ary_push(aval, rb_float_new(pres->cell_flexibility[j]));
 369                                         return aval;
 370                                 }
 371                                 default:
 372                                         rb_raise(rb_eMolbyError, "Internal inconsistency: unknown type field");
 373                         }
 374                 }
 375         }
 376         rb_raise(rb_eMolbyError, "unknown attribute name (%s)", rb_id2name(aid));
 377         return Qnil;  /*  Not reached  */
 378 }
 379
 380 static VALUE
 381 s_MDArena_GetAttr(VALUE self)
 382 {
 383         ID aid = rb_frame_this_func();
 384         return s_MDArena_Get(self, ID2SYM(aid));
 385 }
 386
 387 /*
 388  *  call-seq:
 389  *     self[attr] = value
 390  *
 391  *  Set the attribute value. The attributes can also be modified by self.attribute_name = value.
 392  */
 393 static VALUE
 394 s_MDArena_Set(VALUE self, VALUE attr, VALUE val)
 395 {
 396         MDArena *arena;
 397         int i, j;
 398         struct s_MDArenaAttrDef *dp;
 399         ID aid = rb_to_id(attr);
 400         attr = ID2SYM(aid);  /*  May be used later  */
 401         Data_Get_Struct(self, MDArena, arena);
 402         if (aid == SYM2ID(s_EwaldGridSym)) {
 403                 if (rb_obj_is_kind_of(val, rb_cNumeric)) {
 404                         i = NUM2INT(rb_Integer(val));
 405                         if (i <= 0)
 406                                 rb_raise(rb_eMolbyError, "The ewald grid must be positive integer");
 407                         arena->ewald_grid_x = arena->ewald_grid_y = arena->ewald_grid_z = i;
 408                 } else {
 409                         int ival[3];
 410                         val = rb_ary_to_ary(val);
 411                         j = RARRAY_LEN(val);
 412                         if (j < 3)
 413                                 rb_raise(rb_eMolbyError, "The ewald grid must be an integer or an array of three integers");
 414                         for (i = 0; i < 3; i++) {
 415                                 ival[i] = NUM2INT(rb_Integer(RARRAY_PTR(val)[i]));
 416                                 if (ival[i] <= 0)
 417                                         rb_raise(rb_eMolbyError, "The ewald grid must be positive integer");
 418                         }
 419                         arena->ewald_grid_x = ival[0];
 420                         arena->ewald_grid_y = ival[1];
 421                         arena->ewald_grid_z = ival[2];
 422                 }
 423                 return val;
 424         }
 425         for (i = 0, dp = s_MDArenaAttrDefTable; dp->name != NULL; i++, dp++) {
 426                 if (dp->id == aid) {
 427                         char *p = (char *)arena + dp->offset;
 428                         switch (dp->type) {
 429                                 case 's': {
 430                                         const char *cp = (val == Qnil ? NULL : (const char *)strdup(FileStringValuePtr(val)));
 431                                         const char **cpp = (const char **)p;
 432                                         FILE **fpp;
 433                                         if (*cpp == cp || (*cpp != NULL && cp != NULL && strcmp(*cpp, cp) == 0))
 434                                                 return val;  /*  No need to change  */
 435                                         if (*cpp != NULL)
 436                                                 free((void *)*cpp);
 437                                         if (cp != NULL && cp[0] == 0) {
 438                                                 free((void *)cp);
 439                                                 cp = NULL;
 440                                         }
 441                                         /*  Close the corresponding FILE if necessary  */
 442                                         if (attr == s_LogFileSym)
 443                                                 fpp = &(arena->log_result);
 444                                         else if (attr == s_CoordFileSym)
 445                                                 fpp = &(arena->coord_result);
 446                                         else if (attr == s_VelFileSym)
 447                                                 fpp = &(arena->vel_result);
 448                                         else if (attr == s_ForceFileSym)
 449                                                 fpp = &(arena->force_result);
 450                                         else if (attr == s_DebugFileSym)
 451                                                 fpp = &(arena->debug_result);
 452                                         else fpp = NULL;
 453                                         if (fpp != NULL && *fpp != NULL) {
 454                                                 fclose(*fpp);
 455                                                 *fpp = NULL;
 456                                         }
 457                                         *cpp = cp;
 458                                         return val;
 459                                 }
 460                                 case 'i':
 461                                         *((Int *)p) = NUM2INT(rb_Integer(val));
 462                                         return val;
 463                                 case 'b':
 464                                         *((Byte *)p) = NUM2INT(rb_Integer(val));
 465                                         return val;
 466                                 case 'f':
 467                                         *((Double *)p) = NUM2DBL(rb_Float(val));
 468                                         return val;
 469                                 case 'e':
 470                                         *((Double *)p) = NUM2DBL(rb_Float(val) * KCAL2INTERNAL);
 471                                         return val;
 472                                 case 'S': case 'I': case 'B': case 'F': case 'E':
 473                                         rb_raise(rb_eMolbyError, "The attribute '%s' is read-only", rb_id2name(aid));
 474                                 default:
 475                                         rb_raise(rb_eMolbyError, "Internal inconsistency: unknown type field");
 476                         }
 477                 }
 478         }
 479         for (i = 0, dp = s_MDPressureAttrDefTable; dp->name != NULL; i++, dp++) {
 480                 if (dp->id == aid) {
 481                         char *pp;
 482                         MDPressureArena *pres = arena->pressure;
 483                         if (pres == NULL)
 484                                 arena->pressure = pres = pressure_new();
 485                         pp = (char *)pres + dp->offset;
 486                         switch (dp->type) {
 487                                 case 'i':
 488                                         *((Int *)pp) = NUM2INT(rb_Integer(val));
 489                                         return val;
 490                                 case 'b':
 491                                         *((Byte *)pp) = NUM2INT(rb_Integer(val));
 492                                         return val;
 493                                 case 'f':
 494                                         *((Double *)pp) = NUM2DBL(rb_Float(val));
 495                                         return val;
 496                                 case 'e':
 497                                         *((Double *)pp) = NUM2DBL(rb_Float(val) * KCAL2INTERNAL);
 498                                         return val;
 499                                 case 'X':
 500                                         /*  Isotropic pressure only  */
 501                                         val = rb_ary_to_ary(val);
 502                                         memset(pres->apply, 0, sizeof(Mat33));
 503                                         for (j = 0; j < 3 && j < RARRAY_LEN(val); j++)
 504                                                 pres->apply[j * 4] = NUM2DBL(rb_Float(RARRAY_PTR(val)[j]));
 505                                         return val;
 506                                 case 'Y':
 507                                         val = rb_ary_to_ary(val);
 508                                         for (j = 0; j < 8; j++) {
 509                                                 if (j < RARRAY_LEN(val))
 510                                                         pres->cell_flexibility[j] = NUM2DBL(rb_Float(RARRAY_PTR(val)[j]));
 511                                                 else pres->cell_flexibility[j] = 0.0;
 512                                         }
 513                                         return val;
 514                                 case 'S': case 'I': case 'B': case 'F': case 'E':
 515                                         rb_raise(rb_eMolbyError, "The attribute '%s' is read-only", rb_id2name(aid));
 516                                 default:
 517                                         rb_raise(rb_eMolbyError, "Internal inconsistency: unknown type field");
 518                         }
 519                 }
 520         }
 521         rb_raise(rb_eMolbyError, "unknown attribute name (%s)", rb_id2name(aid));
 522         return Qnil;  /*  Not reached  */
 523 }
 524
 525 static VALUE
 526 s_MDArena_SetAttr(VALUE self, VALUE val)
 527 {
 528         int i;
 529         struct s_MDArenaAttrDef *dp;
 530 //      ID aid = ruby_frame->last_func;
 531         ID aid = rb_frame_this_func();
 532         for (i = 0, dp = s_MDArenaAttrDefTable; dp->name != NULL; i++, dp++) {
 533                 if (dp->sid == aid)
 534                         return s_MDArena_Set(self, *(dp->symref), val);
 535         }
 536         for (i = 0, dp = s_MDPressureAttrDefTable; dp->name != NULL; i++, dp++) {
 537                 if (dp->sid == aid)
 538                         return s_MDArena_Set(self, *(dp->symref), val);
 539         }
 540         rb_raise(rb_eMolbyError, "unknown attribute name (%s)", rb_id2name(aid));
 541         return Qnil;  /*  Not reached  */
 542 }
 543
 544 /*
 545  *  call-seq:
 546  *     to_hash -> Hash
 547  *
 548  *  Returns a (frozen) hash that contains the current value for all attribute keys.
 549  */
 550 static VALUE
 551 s_MDArena_ToHash(VALUE self)
 552 {
 553         int i;
 554         VALUE hash;
 555         struct s_MDArenaAttrDef *dp;
 556         hash = rb_hash_new();
 557         for (i = 0, dp = s_MDArenaAttrDefTable; dp->name != NULL; i++, dp++) {
 558                 VALUE attr = ID2SYM(dp->id);
 559                 rb_hash_aset(hash, attr, s_MDArena_Get(self, attr));
 560         }
 561         rb_obj_freeze(hash);
 562         return hash;
 563 }
 564
 565 /*
 566  *  call-seq:
 567  *     set_alchemical_perturbation(group1, group2) -> [group1, group2]
 568  *
 569  *  Set vanishing and appearing atom groups for alchemical perturbation.
 570  */
 571 static VALUE
 572 s_MDArena_SetAlchemicalPerturbation(VALUE self, VALUE gval1, VALUE gval2)
 573 {
 574         IntGroup *ig1, *ig2;
 575         MDArena *arena;
 576         char *flags;
 577         int i, n;
 578         Data_Get_Struct(self, MDArena, arena);
 579         if (gval1 == Qnil && gval2 == Qnil) {
 580                 md_set_alchemical_flags(arena, 0, NULL);
 581                 return Qnil;
 582         }
 583         if (arena->mol == NULL)
 584                 rb_raise(rb_eMolbyError, "Molecule is not set");
 585         n = arena->xmol->natoms;
 586         flags = (char *)calloc(1, n);
 587         ig1 = (gval1 == Qnil ? NULL : IntGroupFromValue(gval1));
 588         ig2 = (gval2 == Qnil ? NULL : IntGroupFromValue(gval2));
 589         for (i = 0; i < n; i++) {
 590                 if (ig1 != NULL && IntGroupLookupPoint(ig1, i) >= 0)
 591                         flags[i] = 1;
 592                 if (ig2 != NULL && IntGroupLookupPoint(ig2, i) >= 0) {
 593                         if (flags[i] == 1)
 594                                 rb_raise(rb_eMolbyError, "duplicate atom (%d) in vanishing and appearing groups", i);
 595                         flags[i] = 2;
 596                 }
 597         }
 598         if (md_set_alchemical_flags(arena, n, flags) != 0)
 599                 rb_raise(rb_eMolbyError, "cannot set alchemical flags");
 600         free(flags);
 601         if (ig1 != NULL) {
 602                 gval1 = ValueFromIntGroup(ig1);
 603                 IntGroupRelease(ig1);
 604         }
 605         if (ig2 != NULL) {
 606                 gval2 = ValueFromIntGroup(ig2);
 607                 IntGroupRelease(ig2);
 608         }
 609         return rb_ary_new3(2, gval1, gval2);
 610 }
 611
 612 /*
 613  *  call-seq:
 614  *     get_alchemical_perturbation -> [group1, group2] or nil
 615  *
 616  *  If alchemical perturbation is enabled, get the current vanishing and appearing atom groups.
 617  *  Otherwise, return nil.
 618  */
 619 static VALUE
 620 s_MDArena_GetAlchemicalPerturbation(VALUE self)
 621 {
 622         IntGroup *ig1, *ig2;
 623         VALUE gval1, gval2;
 624         MDArena *arena;
 625         int i;
 626         Data_Get_Struct(self, MDArena, arena);
 627         if (arena->nalchem_flags == 0)
 628                 return Qnil;
 629         if (arena->mol == NULL)
 630                 rb_raise(rb_eMolbyError, "Molecule is not set");
 631         ig1 = IntGroupNew();
 632         ig2 = IntGroupNew();
 633         for (i = 0; i < arena->nalchem_flags; i++) {
 634                 if (arena->alchem_flags[i] == 1)
 635                         IntGroupAdd(ig1, i, 1);
 636                 else if (arena->alchem_flags[i] == 2)
 637                         IntGroupAdd(ig2, i, 1);
 638         }
 639         gval1 = ValueFromIntGroup(ig1);
 640         gval2 = ValueFromIntGroup(ig2);
 641         IntGroupRelease(ig1);
 642         IntGroupRelease(ig2);
 643         return rb_ary_new3(2, gval1, gval2);
 644 }
 645
 646 /*
 647  *  call-seq:
 648  *    set_external_forces(ary) -> self
 649  *
 650  *  Set external forces. Ary should be an array of objects that can be converted to Vector3D.
 651  */
 652 static VALUE
 653 s_MDArena_SetExternalForces(VALUE self, VALUE aval)
 654 {
 655         Vector *vp;
 656         MDArena *arena;
 657         int i, n;
 658         Data_Get_Struct(self, MDArena, arena);
 659         if (arena->mol == NULL)
 660                 rb_raise(rb_eMolbyError, "Molecule is not set");
 661         if (aval == Qnil) {
 662                 md_set_external_forces(arena, 0, NULL);
 663                 return self;
 664         }
 665         aval = rb_ary_to_ary(aval);
 666         n = RARRAY_LEN(aval);
 667         if (n == 0) {
 668                 md_set_external_forces(arena, 0, NULL);
 669                 return self;
 670         }
 671         vp = (Vector *)calloc(sizeof(Vector), n);
 672         for (i = 0; i < n; i++) {
 673                 VectorFromValue(RARRAY_PTR(aval)[i], vp + i);
 674                 VecScaleSelf(vp[i], KCAL2INTERNAL);
 675         }
 676         md_set_external_forces(arena, n, vp);
 677         free(vp);
 678         return self;
 679 }
 680
 681 /*
 682  *  call-seq:
 683  *    get_external_force(index) -> Vector3D or nil
 684  *
 685  *  Get the current external force for the atom. If the external force is not set, nil is returned.
 686  */
 687 static VALUE
 688 s_MDArena_GetExternalForce(VALUE self, VALUE ival)
 689 {
 690         int i;
 691         VALUE vval;
 692         Vector v;
 693         MDArena *arena;
 694         Data_Get_Struct(self, MDArena, arena);
 695         if (arena->mol == NULL)
 696                 rb_raise(rb_eMolbyError, "Molecule is not set");
 697         i = NUM2INT(rb_Integer(ival));
 698         if (i < 0 || i >= arena->nexforces)
 699                 return Qnil;
 700         v = arena->exforces[i];
 701         VecScaleSelf(v, INTERNAL2KCAL);
 702         vval = ValueFromVector(&v);
 703         return vval;
 704 }
 705
 706 /*
 707  *  call-seq:
 708  *     init_velocities([temperature]) -> self
 709  *
 710  *  Give random (Boltzmann-weighted) velocities to all atoms. If temperature is given,
 711  *  it is set before giving velocities.
 712  */
 713 static VALUE
 714 s_MDArena_InitVelocities(int argc, VALUE *argv, VALUE self)
 715 {
 716         MDArena *arena;
 717         VALUE tval;
 718         Data_Get_Struct(self, MDArena, arena);
 719         rb_scan_args(argc, argv, "01", &tval);
 720         if (tval != Qnil)
 721                 s_MDArena_Set(self, s_TemperatureSym, tval);
 722         md_init_velocities(arena);
 723         return self;
 724 }
 725
 726 /*
 727  *  call-seq:
 728  *     scale_velocities([temperature]) -> self
 729  *
 730  *  Scale atom velocities to match the current temperature. If temperature is given,
 731  *  it is set before giving velocities.
 732  */
 733 static VALUE
 734 s_MDArena_ScaleVelocities(int argc, VALUE *argv, VALUE self)
 735 {
 736         MDArena *arena;
 737         VALUE tval;
 738         Data_Get_Struct(self, MDArena, arena);
 739         rb_scan_args(argc, argv, "01", &tval);
 740         if (tval != Qnil)
 741                 s_MDArena_Set(self, s_TemperatureSym, tval);
 742         md_scale_velocities(arena);
 743         return self;
 744 }
 745
 746 /*
 747  *  call-seq:
 748  *     keys -> Array
 749  *
 750  *  Returns an array of valid attributes.
 751  */
 752 static VALUE
 753 s_MDArena_Keys(VALUE self)
 754 {
 755         int i;
 756         VALUE ary;
 757         struct s_MDArenaAttrDef *dp;
 758         ary = rb_ary_new();
 759         for (i = 0, dp = s_MDArenaAttrDefTable; dp->name != NULL; i++, dp++) {
 760                 VALUE attr = ID2SYM(dp->id);
 761                 rb_ary_push(ary, attr);
 762         }
 763         return ary;
 764 }
 765
 766 /*
 767  *  call-seq:
 768  *     print_surface_area
 769  *
 770  *  Print the surface area information to standard output. (for debug)
 771  */
 772 static VALUE
 773 s_MDArena_PrintSurfaceArea(VALUE self)
 774 {
 775         MDArena *arena;
 776         Int i, natoms;
 777         VALUE retval, outval;
 778         Data_Get_Struct(self, MDArena, arena);
 779         if (arena->sp_arena == NULL)
 780                 rb_raise(rb_eMolbyError, "surface potential is not available");
 781         natoms = arena->mol->natoms;
 782         retval = Ruby_NewEncodedStringValue2("Atom     area    energy         forcex1000\n");
 783         for (i = 0; i < natoms; i++) {
 784                 char buf[256];
 785                 Vector f = arena->forces[kSurfaceIndex * natoms + i];
 786                 Double area = arena->sp_arena->atom_area[i];
 787                 Double energy = area * arena->sp_arena->atom_pot[i];
 788                 VecScaleSelf(f, INTERNAL2KCAL * 1000);
 789                 energy *= INTERNAL2KCAL;
 790                 snprintf(buf, sizeof buf, "%5d %11.5f %11.8f %11.5f %11.5f %11.5f\n",
 791                            i+1, area, energy, f.x, f.y, f.z);
 792                 rb_str_cat(retval, buf, strlen(buf));
 793         }
 794         outval = rb_gv_get("$stdout");
 795         rb_funcall(outval, rb_intern("write"), 1, retval);
 796         return self;
 797 }
 798
 799 /*
 800  *  call-seq:
 801  *     bond_par(idx) -> ParameterRef
 802  *
 803  *  Returns a parameter that is used for the idx-th bond.
 804  */
 805 static VALUE
 806 s_MDArena_BondPar(VALUE self, VALUE val)
 807 {
 808         MDArena *arena;
 809         Int i, j;
 810         Data_Get_Struct(self, MDArena, arena);
 811         i = NUM2INT(rb_Integer(val));
 812         if (arena->xmol->needsMDRebuild || arena->par == NULL)
 813                 md_prepare(arena, 1);
 814         if (i < 0 || i >= arena->xmol->nbonds) {
 815                 rb_raise(rb_eMolbyError, "bond index (%d) out of range (0...%d)", i, arena->xmol->nbonds);
 816         }
 817         j = arena->bond_par_i[i];
 818         if (j < 0 || j >= arena->par->nbondPars)
 819                 return Qnil;  /*  No parameter assigned  */
 820         return ValueFromMoleculeWithParameterTypeAndIndex(arena->xmol, kBondParType, -j - 1);
 821 }
 822
 823 /*
 824  *  call-seq:
 825  *     angle_par(idx) -> ParameterRef
 826  *
 827  *  Returns a parameter that is used for the idx-th angle.
 828  */
 829 static VALUE
 830 s_MDArena_AnglePar(VALUE self, VALUE val)
 831 {
 832         MDArena *arena;
 833         Int i, j;
 834         Data_Get_Struct(self, MDArena, arena);
 835         i = NUM2INT(rb_Integer(val));
 836         if (arena->xmol->needsMDRebuild || arena->par == NULL)
 837                 md_prepare(arena, 1);
 838         if (i < 0 || i >= arena->xmol->nangles) {
 839                 rb_raise(rb_eMolbyError, "angle index (%d) out of range (0...%d)", i, arena->xmol->nangles);
 840         }
 841         j = arena->angle_par_i[i];
 842         if (j < 0 || j >= arena->par->nanglePars)
 843                 return Qnil;  /*  No parameter assigned  */
 844         return ValueFromMoleculeWithParameterTypeAndIndex(arena->xmol, kAngleParType, -j - 1);
 845 }
 846
 847 /*
 848  *  call-seq:
 849  *     dihedral_par(idx) -> ParameterRef
 850  *
 851  *  Returns a parameter that is used for the idx-th dihedral.
 852  */
 853 static VALUE
 854 s_MDArena_DihedralPar(VALUE self, VALUE val)
 855 {
 856         MDArena *arena;
 857         Int i, j;
 858         Data_Get_Struct(self, MDArena, arena);
 859         i = NUM2INT(rb_Integer(val));
 860         if (arena->xmol->needsMDRebuild || arena->par == NULL)
 861                 md_prepare(arena, 1);
 862         if (i < 0 || i >= arena->xmol->ndihedrals) {
 863                 rb_raise(rb_eMolbyError, "dihedral index (%d) out of range (0...%d)", i, arena->xmol->ndihedrals);
 864         }
 865         j = arena->dihedral_par_i[i];
 866         if (j < 0 || j >= arena->par->ndihedralPars)
 867                 return Qnil;  /*  No parameter assigned  */
 868         return ValueFromMoleculeWithParameterTypeAndIndex(arena->xmol, kDihedralParType, -j - 1);
 869 }
 870
 871 /*
 872  *  call-seq:
 873  *     improper_par(idx) -> ParameterRef
 874  *
 875  *  Returns a parameter that is used for the idx-th improper.
 876  */
 877 static VALUE
 878 s_MDArena_ImproperPar(VALUE self, VALUE val)
 879 {
 880         MDArena *arena;
 881         Int i, j;
 882         Data_Get_Struct(self, MDArena, arena);
 883         i = NUM2INT(rb_Integer(val));
 884         if (arena->xmol->needsMDRebuild || arena->par == NULL)
 885                 md_prepare(arena, 1);
 886         if (i < 0 || i >= arena->xmol->nimpropers) {
 887                 rb_raise(rb_eMolbyError, "improper index (%d) out of range (0...%d)", i, arena->xmol->nimpropers);
 888         }
 889         j = arena->improper_par_i[i];
 890         if (j < 0 || j >= arena->par->nimproperPars)
 891                 return Qnil;  /*  No parameter assigned  */
 892         return ValueFromMoleculeWithParameterTypeAndIndex(arena->xmol, kImproperParType, -j - 1);
 893 }
 894
 895 /*
 896  *  call-seq:
 897  *     vdw_par(idx) -> ParameterRef
 898  *
 899  *  Returns a vdw parameter that is used for the idx-th atom.
 900  */
 901 static VALUE
 902 s_MDArena_VdwPar(VALUE self, VALUE val)
 903 {
 904         MDArena *arena;
 905         Int i, j;
 906         Data_Get_Struct(self, MDArena, arena);
 907         i = NUM2INT(rb_Integer(val));
 908         if (arena->xmol->needsMDRebuild || arena->par == NULL)
 909                 md_prepare(arena, 1);
 910         if (i < 0 || i >= arena->xmol->natoms) {
 911                 rb_raise(rb_eMolbyError, "atom index (%d) out of range (0...%d)", i, arena->xmol->natoms);
 912         }
 913         j = arena->vdw_par_i[i];
 914         if (j < 0 || j >= arena->par->nvdwPars)
 915                 return Qnil;  /*  No parameter assigned  */
 916         return ValueFromMoleculeWithParameterTypeAndIndex(arena->xmol, kVdwParType, -j - 1);
 917 }
 918
 919 static VALUE
 920 s_MDArena_testPME(int argc, VALUE *argv, VALUE self)
 921 {
 922         extern int pme_test(MDArena *);
 923         MDArena *arena;
 924         if (rb_obj_is_kind_of(self, rb_cMDArena)) {
 925                 Data_Get_Struct(self, MDArena, arena);
 926                 if (argc == 0 || RTEST(argv[0])) {
 927                         arena->use_ewald = 1;
 928                         arena->xmol->needsMDRebuild = 1;
 929                         if (argc > 0 && rb_obj_is_kind_of(argv[0], rb_cNumeric)) {
 930                                 Int n = NUM2INT(rb_Integer(argv[0]));
 931                                 if (n > 1) {
 932                                         arena->ewald_grid_x = n;
 933                                         arena->ewald_grid_y = n;
 934                                         arena->ewald_grid_z = n;
 935                                 } else {
 936                                         arena->use_ewald = 2;  /*  Direct Ewald  */
 937                                 }
 938                         }
 939                 } else {
 940                         arena->use_ewald = 0;
 941                         arena->xmol->needsMDRebuild = 1;
 942                 }
 943                 s_MDArena_Prepare(0, NULL, self);
 944         } else {
 945                 rb_raise(rb_eMolbyError, "class method MDArena.test_pme is not supported now");
 946         }
 947         return self;
 948 }
 949
 950 void
 951 Init_MolbyMDTypes(void)
 952 {
 953         int i;
 954         struct s_MDArenaAttrDef *dp;
 955         char name[40];
 956
 957         /*  class MDArena  */
 958         rb_cMDArena = rb_define_class("MDArena", rb_cObject);
 959 /*      rb_define_alloc_func(rb_cMDArena, s_MDArena_Alloc);
 960     rb_define_private_method(rb_cMDArena, "initialize", s_MDArena_Initialize, 1); */
 961         rb_define_method(rb_cMDArena, "run", s_MDArena_Run, 1);
 962     rb_define_method(rb_cMDArena, "minimize", s_MDArena_Minimize, 1);
 963     rb_define_method(rb_cMDArena, "prepare", s_MDArena_Prepare, -1);
 964     rb_define_method(rb_cMDArena, "energies", s_MDArena_Energies, 0);
 965         rb_define_method(rb_cMDArena, "[]", s_MDArena_Get, 1);
 966         rb_define_method(rb_cMDArena, "[]=", s_MDArena_Set, 2);
 967         rb_define_method(rb_cMDArena, "to_hash", s_MDArena_ToHash, 0);
 968         rb_define_method(rb_cMDArena, "keys", s_MDArena_Keys, 0);
 969         rb_define_method(rb_cMDArena, "set_alchemical_perturbation", s_MDArena_SetAlchemicalPerturbation, 2);
 970         rb_define_method(rb_cMDArena, "get_alchemical_perturbation", s_MDArena_GetAlchemicalPerturbation, 0);
 971         rb_define_method(rb_cMDArena, "set_external_forces", s_MDArena_SetExternalForces, 1);
 972         rb_define_method(rb_cMDArena, "get_external_force", s_MDArena_GetExternalForce, 1);
 973         rb_define_method(rb_cMDArena, "init_velocities", s_MDArena_InitVelocities, -1);
 974         rb_define_method(rb_cMDArena, "scale_velocities", s_MDArena_ScaleVelocities, -1);
 975         rb_define_method(rb_cMDArena, "print_surface_area", s_MDArena_PrintSurfaceArea, 0);
 976
 977         rb_define_method(rb_cMDArena, "bond_par", s_MDArena_BondPar, 1);
 978         rb_define_method(rb_cMDArena, "angle_par", s_MDArena_AnglePar, 1);
 979         rb_define_method(rb_cMDArena, "dihedral_par", s_MDArena_DihedralPar, 1);
 980         rb_define_method(rb_cMDArena, "improper_par", s_MDArena_ImproperPar, 1);
 981         rb_define_method(rb_cMDArena, "vdw_par", s_MDArena_VdwPar, 1);
 982
 983         rb_define_method(rb_cMDArena, "pme_test", s_MDArena_testPME, -1);
 984 //      rb_define_singleton_method(rb_cMDArena, "pme_test", s_MDArena_testPME, -1);
 985
 986         /*  All setter and getter are handled with the same C function (attribute name is taken
 987             from ruby_frame)  */
 988         for (i = 0, dp = s_MDArenaAttrDefTable; dp->name != NULL; i++, dp++) {
 989                 rb_define_method(rb_cMDArena, dp->name, s_MDArena_GetAttr, 0);
 990                 strncpy(name, dp->name, 38);
 991                 name[38] = 0;
 992                 strcat(name, "=");
 993                 rb_define_method(rb_cMDArena, name, s_MDArena_SetAttr, 1);
 994                 dp->id = rb_intern(dp->name);
 995                 dp->sid = rb_intern(name);
 996                 *(dp->symref) = ID2SYM(dp->id);
 997         }
 998         for (i = 0, dp = s_MDPressureAttrDefTable; dp->name != NULL; i++, dp++) {
 999                 rb_define_method(rb_cMDArena, dp->name, s_MDArena_GetAttr, 0);
1000                 strncpy(name, dp->name, 38);
1001                 name[38] = 0;
1002                 strcat(name, "=");
1003                 rb_define_method(rb_cMDArena, name, s_MDArena_SetAttr, 1);
1004                 dp->id = rb_intern(dp->name);
1005                 dp->sid = rb_intern(name);
1006                 *(dp->symref) = ID2SYM(dp->id);
1007         }
1008 }