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         IntGroup *ig;
  66         Data_Get_Struct(self, MDArena, arena);
  67         nsteps = NUM2INT(rb_Integer(arg));
  68         if (arena->is_running)
  69                 rb_raise(rb_eMolbyError, "the simulation is already running. You cannot do simulation recursively.");
  70         if (nsteps < 0)
  71                 rb_raise(rb_eMolbyError, "the number of steps should be non-negative integer");
  72         if (!arena->is_initialized || 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 /*      arena->md_panic_func = s_MDArena_panic_func; */
  90         retval = md_main(arena, minimize);
  91         if (retval != 0)
  92                 rb_raise(rb_eMolbyError, "MD Error: %s", arena->errmsg);
  93
  94 /*      arena->md_panic_func = NULL; */
  95
  96         if (arena->step > start_step) {
  97                 /*  Create a new frame and copy new coordinates  */
  98                 ig = IntGroupNewWithPoints(MoleculeGetNumberOfFrames(arena->xmol), 1, -1);
  99                 MolActionCreateAndPerform(arena->xmol, gMolActionInsertFrames, ig, 0, NULL);
 100                 IntGroupRelease(ig);
 101                 md_copy_coordinates_from_internal(arena);
 102         }
 103
 104 /*      if (retval != 0)
 105                 rb_raise(rb_eMolbyError, "Calculation aborted with status %d", retval); */
 106
 107         if (minimize && arena->minimize_complete && rb_block_given_p())
 108                 rb_yield(self);
 109
 110         return self;
 111 }
 112
 113 /*
 114  *  call-seq:
 115  *     run(n)       -> self
 116  *
 117  *  Run the simulation for n steps.
 118  */
 119 static VALUE
 120 s_MDArena_Run(VALUE self, VALUE arg)
 121 {
 122         return s_MDArena_Run_or_minimize(self, arg, 0);
 123 }
 124
 125 /*
 126  *  call-seq:
 127  *     minimize(n)       -> self
 128  *     minimize(n) { ... } -> self
 129  *
 130  *  Minimize the energy for n steps. If a block is given, it is executed when minimization is complete.
 131  */
 132 static VALUE
 133 s_MDArena_Minimize(VALUE self, VALUE arg)
 134 {
 135         return s_MDArena_Run_or_minimize(self, arg, 1);
 136 }
 137
 138 /*
 139  *  call-seq:
 140  *     prepare(check_only = false)         -> self or nil
 141  *
 142  *  Prepare for the MD calculation; refresh the internal information even if initialization is
 143  *  already done.
 144  *  If check_only is true, then only the parameter checking is done.
 145  *  Returns self when initialization is complete; returns nil if some MM parameters are missing;
 146  *  otherwise, an exception is raised.
 147  */
 148 static VALUE
 149 s_MDArena_Prepare(int argc, VALUE *argv, VALUE self)
 150 {
 151         MDArena *arena;
 152         Molecule *mol;
 153         const char *msg;
 154         Int nangles, *angles, ndihedrals, *dihedrals, nimpropers, *impropers;
 155         Int missing = 0;
 156         Int check_only = 0;
 157         IntGroup *ig1, *ig2, *ig3;
 158         VALUE fval;
 159
 160         Data_Get_Struct(self, MDArena, arena);
 161         rb_scan_args(argc, argv, "01", &fval);
 162         if (RTEST(fval))
 163                 check_only = 1;
 164
 165         arena->is_initialized = 0;
 166         mol = arena->xmol;
 167
 168         /*  Rebuild the tables  */
 169         ig1 = ig2 = ig3 = NULL;
 170         nangles = MoleculeFindMissingAngles(mol, &angles);
 171         ndihedrals = MoleculeFindMissingDihedrals(mol, &dihedrals);
 172         nimpropers = MoleculeFindMissingImpropers(mol, &impropers);
 173         if (nangles > 0) {
 174                 ig1 = IntGroupNewWithPoints(mol->nangles, nangles, -1);
 175                 MolActionCreateAndPerform(mol, gMolActionAddAngles, nangles * 3, angles, ig1);
 176                 free(angles);
 177                 IntGroupRelease(ig1);
 178         }
 179         if (ndihedrals > 0) {
 180                 ig2 = IntGroupNewWithPoints(mol->ndihedrals, ndihedrals, -1);
 181                 MolActionCreateAndPerform(mol, gMolActionAddDihedrals, ndihedrals * 4, dihedrals, ig2);
 182                 free(dihedrals);
 183                 IntGroupRelease(ig2);
 184         }
 185         if (nimpropers > 0) {
 186                 ig3 = IntGroupNewWithPoints(mol->nimpropers, nimpropers, -1);
 187                 MolActionCreateAndPerform(mol, gMolActionAddImpropers, nimpropers * 4, impropers, ig3);
 188                 free(impropers);
 189                 IntGroupRelease(ig3);
 190         }
 191
 192         /*  Prepare parameters and internal information  */
 193         msg = md_prepare(arena, check_only);
 194
 195         /*  Some parameters are missing?  */
 196         if (msg != NULL) {
 197                 if (strstr(msg, "parameter") != NULL && strstr(msg, "missing") != NULL)
 198                         missing = 1;
 199                 else
 200                         rb_raise(rb_eMolbyError, "cannot initialize for MD: %s", msg);
 201         }
 202
 203         /*  The local parameter list is updated  */
 204         {
 205                 Int parType, idx;
 206                 if (mol->par == NULL)
 207                         mol->par = ParameterNew();
 208                 for (parType = kFirstParType; parType <= kLastParType; parType++) {
 209                         /*  Delete global and undefined parameters  */
 210                         UnionPar *up, *upbuf;
 211                         Int nparams, count;
 212                         ig1 = IntGroupNew();
 213                         for (idx = 0; (up = ParameterGetUnionParFromTypeAndIndex(mol->par, parType, idx)) != NULL; idx++) {
 214                                 if (up->bond.src != 0)
 215                                         IntGroupAdd(ig1, idx, 1);
 216                         }
 217                         if (IntGroupGetCount(ig1) > 0)
 218                                 MolActionCreateAndPerform(mol, gMolActionDeleteParameters, parType, ig1);
 219                         IntGroupRelease(ig1);
 220                         /*  Copy global and undefined parameters from arena and insert to mol->par  */
 221                         nparams = ParameterGetCountForType(arena->par, parType);
 222                         if (nparams == 0)
 223                                 continue;
 224                         upbuf = (UnionPar *)calloc(sizeof(UnionPar), nparams);
 225                         ig1 = IntGroupNew();
 226                         ig2 = IntGroupNew();
 227                         for (idx = 0; (up = ParameterGetUnionParFromTypeAndIndex(arena->par, parType, idx)) != NULL; idx++) {
 228                                 if (up->bond.src > 0)
 229                                         IntGroupAdd(ig1, idx, 1); /* Global parameter */
 230                                 else if (up->bond.src < 0)
 231                                         IntGroupAdd(ig2, idx, 1); /* Undefined parameter */
 232                         }
 233                         if ((count = IntGroupGetCount(ig1)) > 0) {
 234                                 /*  Insert global parameters (at the top)  */
 235                                 ParameterCopy(arena->par, parType, upbuf, ig1);
 236                                 ig3 = IntGroupNewWithPoints(0, count, -1);
 237                                 MolActionCreateAndPerform(mol, gMolActionAddParameters, parType, ig3, count, upbuf);
 238                                 IntGroupRelease(ig3);
 239                         }
 240                         if ((count = IntGroupGetCount(ig2)) > 0) {
 241                                 /*  Insert undefined parameters (at the bottom)  */
 242                                 ParameterCopy(arena->par, parType, upbuf, ig2);
 243                                 idx = ParameterGetCountForType(mol->par, parType);
 244                                 ig3 = IntGroupNewWithPoints(idx, count, -1);
 245                                 MolActionCreateAndPerform(mol, gMolActionAddParameters, parType, ig3, count, upbuf);
 246                                 IntGroupRelease(ig3);
 247                         }
 248                         IntGroupRelease(ig2);
 249                         IntGroupRelease(ig1);
 250                         free(upbuf);
 251                 }
 252         }
 253
 254         if (missing)
 255                 return Qnil;
 256         else return self;
 257 }
 258
 259 /*
 260  *  call-seq:
 261  *     energies -> [total, bond, angle, dihedral, improper, vdw, electrostatic, auxiliary, surface, kinetic, net]
 262  *
 263  *  Get the current energies.
 264  */
 265 static VALUE
 266 s_MDArena_Energies(VALUE self)
 267 {
 268         MDArena *arena;
 269         VALUE val;
 270         int i;
 271         Data_Get_Struct(self, MDArena, arena);
 272         val = rb_ary_new();
 273         rb_ary_push(val, rb_float_new(arena->total_energy * INTERNAL2KCAL));
 274         for (i = 0; i < kEndIndex; i++)
 275                 rb_ary_push(val, rb_float_new(arena->energies[i] * INTERNAL2KCAL));
 276         rb_ary_push(val, rb_float_new((arena->energies[kKineticIndex] + arena->total_energy) * INTERNAL2KCAL));
 277         return val;
 278 }
 279
 280 static VALUE s_LogFileSym, s_CoordFileSym, s_VelFileSym, s_ForceFileSym,
 281 s_DebugFileSym, s_DebugOutputLevelSym, s_StepSym, s_CoordOutputFreqSym,
 282 s_EnergyOutputFreqSym, s_CoordFrameSym, s_TimestepSym, s_CutoffSym,
 283 s_ElectroCutoffSym, s_PairlistDistanceSym, s_TemperatureSym, s_TransientTempSym,
 284 s_AverageTempSym, s_AndersenFreqSym, s_AndersenCouplingSym, s_RandomSeedSym,
 285 s_DielectricSym, s_GradientConvergenceSym, s_CoordinateConvergenceSym, s_UseXplorShiftSym,
 286 s_Scale14VdwSym, s_Scale14ElectSym, s_RelocateCenterSym, s_SurfaceProbeRadiusSym,
 287 s_SurfaceTensionSym, s_SurfacePotentialFreqSym, s_UseGraphiteSym;
 288
 289 struct s_MDArenaAttrDef {
 290         char *name;
 291         VALUE *symref;  /*  Address of s_LogFileSym etc. */
 292         ID id;                  /*  Ruby ID of the symbol; will be set within Init_MolbyMDTypes()  */
 293         ID sid;         /*  Ruby ID of the symbol plus '='; will be set within Init_MolbyMDTypes()  */
 294         char type;      /*  s: string (const char *), i: Int, f: Double. Uppercase: read-only.  */
 295         int  offset;    /*  Offset in the MDArena structure.  */
 296 };
 297 static struct s_MDArenaAttrDef s_MDArenaAttrDefTable[] = {
 298         {"log_file",          &s_LogFileSym,          0, 0, 's', offsetof(MDArena, log_result_name)},
 299         {"coord_file",        &s_CoordFileSym,        0, 0, 's', offsetof(MDArena, coord_result_name)},
 300         {"vel_file",          &s_VelFileSym,          0, 0, 's', offsetof(MDArena, vel_result_name)},
 301         {"force_file",        &s_ForceFileSym,        0, 0, 's', offsetof(MDArena, force_result_name)},
 302         {"debug_file",        &s_DebugFileSym,        0, 0, 's', offsetof(MDArena, debug_result_name)},
 303         {"debug_output_level", &s_DebugOutputLevelSym, 0, 0, 'i', offsetof(MDArena, debug_output_level)},
 304         {"step",              &s_StepSym,             0, 0, 'I', offsetof(MDArena, step)},
 305         {"coord_output_freq", &s_CoordOutputFreqSym,  0, 0, 'i', offsetof(MDArena, coord_output_freq)},
 306         {"energy_output_freq", &s_EnergyOutputFreqSym, 0, 0, 'i', offsetof(MDArena, energy_output_freq)},
 307         {"coord_frame",       &s_CoordFrameSym,       0, 0, 'I', offsetof(MDArena, coord_result_frame)},
 308         {"timestep",          &s_TimestepSym,         0, 0, 'f', offsetof(MDArena, timestep)},
 309         {"cutoff",            &s_CutoffSym,           0, 0, 'f', offsetof(MDArena, cutoff)},
 310         {"electro_cutoff",    &s_ElectroCutoffSym,    0, 0, 'f', offsetof(MDArena, electro_cutoff)},
 311         {"pairlist_distance", &s_PairlistDistanceSym, 0, 0, 'f', offsetof(MDArena, pairlist_distance)},
 312         {"temperature",       &s_TemperatureSym,      0, 0, 'f', offsetof(MDArena, temperature)},
 313         {"transient_temperature", &s_TransientTempSym, 0, 0, 'F', offsetof(MDArena, transient_temperature)},
 314         {"average_temperature", &s_AverageTempSym,    0, 0, 'F', offsetof(MDArena, average_temperature)},
 315         {"andersen_freq",     &s_AndersenFreqSym,     0, 0, 'i', offsetof(MDArena, andersen_thermo_freq)},
 316         {"andersen_coupling", &s_AndersenCouplingSym, 0, 0, 'f', offsetof(MDArena, andersen_thermo_coupling)},
 317         {"random_seed",       &s_RandomSeedSym,       0, 0, 'i', offsetof(MDArena, random_seed)},
 318         {"dielectric",        &s_DielectricSym,       0, 0, 'f', offsetof(MDArena, dielectric)},
 319         {"gradient_convergence", &s_GradientConvergenceSym, 0, 0, 'f', offsetof(MDArena, gradient_convergence)},
 320         {"coordinate_convergence", &s_CoordinateConvergenceSym, 0, 0, 'f', offsetof(MDArena, coordinate_convergence)},
 321         {"use_xplor_shift",   &s_UseXplorShiftSym,    0, 0, 'i', offsetof(MDArena, use_xplor_shift)},
 322         {"scale14_vdw",       &s_Scale14VdwSym,       0, 0, 'f', offsetof(MDArena, scale14_vdw)},
 323         {"scale14_elect",     &s_Scale14ElectSym,     0, 0, 'f', offsetof(MDArena, scale14_elect)},
 324         {"relocate_center",   &s_RelocateCenterSym,   0, 0, 'i', offsetof(MDArena, relocate_center)},
 325         {"surface_probe_radius", &s_SurfaceProbeRadiusSym, 0, 0, 'f', offsetof(MDArena, probe_radius)},
 326         {"surface_tension",   &s_SurfaceTensionSym,   0, 0, 'f', offsetof(MDArena, surface_tension)},
 327         {"surface_potential_freq", &s_SurfacePotentialFreqSym, 0, 0, 'i', offsetof(MDArena, surface_potential_freq)},
 328         {"use_graphite",      &s_UseGraphiteSym,      0, 0, 'i', offsetof(MDArena, use_graphite)},
 329         {NULL} /* Sentinel */
 330 };
 331
 332 static VALUE s_PresFreqSym, s_PresCouplingSym, s_PresSym, s_PresCellFlexSym, s_PresFluctCellOriginSym, s_PresFluctCellOrientSym;
 333 static struct s_MDArenaAttrDef s_MDPressureAttrDefTable[] = {
 334         {"pressure_freq",     &s_PresFreqSym,         0, 0, 'i', offsetof(MDPressureArena, freq)},
 335         {"pressure_coupling", &s_PresCouplingSym,     0, 0, 'f', offsetof(MDPressureArena, coupling)},
 336         {"pressure",          &s_PresSym,             0, 0, 'X', offsetof(MDPressureArena, apply)},
 337         {"pressure_cell_flexibility", &s_PresCellFlexSym, 0, 0, 'Y', offsetof(MDPressureArena, cell_flexibility)},
 338         {"pressure_fluctuate_cell_origin", &s_PresFluctCellOriginSym, 0, 0, 'f', offsetof(MDPressureArena, fluctuate_cell_origin)},
 339         {"pressure_fluctuate_cell_orientation", &s_PresFluctCellOrientSym, 0, 0, 'f', offsetof(MDPressureArena, fluctuate_cell_orientation)},
 340         {NULL} /* Sentinel */
 341 };
 342
 343 /*
 344  *  call-seq:
 345  *     self[attr]
 346  *
 347  *  Get the attribute value. The attribute values also can be accessed by self.attribute_name.
 348  */
 349 static VALUE
 350 s_MDArena_Get(VALUE self, VALUE attr)
 351 {
 352         MDArena *arena;
 353         int i;
 354         struct s_MDArenaAttrDef *dp;
 355         ID aid = rb_to_id(attr);
 356         Data_Get_Struct(self, MDArena, arena);
 357         for (i = 0, dp = s_MDArenaAttrDefTable; dp->name != NULL; i++, dp++) {
 358                 if (dp->id == aid) {
 359                         char *p = (char *)arena + dp->offset;
 360                         switch (dp->type) {
 361                                 case 's':
 362                                 case 'S': {
 363                                         const char *cp = *((const char **)p);
 364                                         if (cp == NULL)
 365                                                 return Qnil;
 366                                         else
 367                                                 return Ruby_NewFileStringValue(cp);
 368                                 }
 369                                 case 'i':
 370                                 case 'I':
 371                                         return INT2NUM(*((Int *)p));
 372                                 case 'f':
 373                                 case 'F':
 374                                         return rb_float_new(*((Double *)p));
 375                                 default:
 376                                         rb_raise(rb_eMolbyError, "Internal inconsistency: unknown type field");
 377                         }
 378                 }
 379         }
 380         for (i = 0, dp = s_MDPressureAttrDefTable; dp->name != NULL; i++, dp++) {
 381                 if (dp->id == aid) {
 382                         char *pp;
 383                         MDPressureArena *pres = arena->pressure;
 384                         if (pres == NULL)
 385                                 return Qnil;
 386                         pp = (char *)pres + dp->offset;
 387                         switch (dp->type) {
 388                                 case 'i':
 389                                 case 'I':
 390                                         return INT2NUM(*((Int *)pp));
 391                                 case 'f':
 392                                 case 'F':
 393                                         return rb_float_new(*((Double *)pp));
 394                                 case 'X':
 395                                         /*  Isotropic pressure only  */
 396                                         return rb_ary_new3(3, rb_float_new(pres->apply[0]), rb_float_new(pres->apply[4]), rb_float_new(pres->apply[8]));
 397                                 case 'Y': {
 398                                         VALUE aval = rb_ary_new();
 399                                         int j;
 400                                         for (j = 0; j < 8; j++)
 401                                                 rb_ary_push(aval, rb_float_new(pres->cell_flexibility[j]));
 402                                         return aval;
 403                                 }
 404                                 default:
 405                                         rb_raise(rb_eMolbyError, "Internal inconsistency: unknown type field");
 406                         }
 407                 }
 408         }
 409         rb_raise(rb_eMolbyError, "unknown attribute name (%s)", rb_id2name(aid));
 410         return Qnil;  /*  Not reached  */
 411 }
 412
 413 static VALUE
 414 s_MDArena_GetAttr(VALUE self)
 415 {
 416         return s_MDArena_Get(self, ID2SYM(ruby_frame->last_func));
 417 }
 418
 419 /*
 420  *  call-seq:
 421  *     self[attr] = value
 422  *
 423  *  Set the attribute value. The attributes can also be modified by self.attribute_name = value.
 424  */
 425 static VALUE
 426 s_MDArena_Set(VALUE self, VALUE attr, VALUE val)
 427 {
 428         MDArena *arena;
 429         int i, j;
 430         struct s_MDArenaAttrDef *dp;
 431         ID aid = rb_to_id(attr);
 432         attr = ID2SYM(aid);  /*  May be used later  */
 433         Data_Get_Struct(self, MDArena, arena);
 434         for (i = 0, dp = s_MDArenaAttrDefTable; dp->name != NULL; i++, dp++) {
 435                 if (dp->id == aid) {
 436                         char *p = (char *)arena + dp->offset;
 437                         switch (dp->type) {
 438                                 case 's': {
 439                                         const char *cp = (val == Qnil ? NULL : (const char *)strdup(FileStringValuePtr(val)));
 440                                         const char **cpp = (const char **)p;
 441                                         FILE **fpp;
 442                                         if (*cpp == cp || (*cpp != NULL && cp != NULL && strcmp(*cpp, cp) == 0))
 443                                                 return val;  /*  No need to change  */
 444                                         if (*cpp != NULL)
 445                                                 free((void *)*cpp);
 446                                         if (cp != NULL && cp[0] == 0) {
 447                                                 free((void *)cp);
 448                                                 cp = NULL;
 449                                         }
 450                                         /*  Close the corresponding FILE if necessary  */
 451                                         if (attr == s_LogFileSym)
 452                                                 fpp = &(arena->log_result);
 453                                         else if (attr == s_CoordFileSym)
 454                                                 fpp = &(arena->coord_result);
 455                                         else if (attr == s_VelFileSym)
 456                                                 fpp = &(arena->vel_result);
 457                                         else if (attr == s_ForceFileSym)
 458                                                 fpp = &(arena->force_result);
 459                                         else if (attr == s_DebugFileSym)
 460                                                 fpp = &(arena->debug_result);
 461                                         else fpp = NULL;
 462                                         if (fpp != NULL && *fpp != NULL) {
 463                                                 fclose(*fpp);
 464                                                 *fpp = NULL;
 465                                         }
 466                                         *cpp = cp;
 467                                         return val;
 468                                 }
 469                                 case 'i':
 470                                         *((Int *)p) = NUM2INT(rb_Integer(val));
 471                                         return val;
 472                                 case 'f':
 473                                         *((Double *)p) = NUM2DBL(rb_Float(val));
 474                                         return val;
 475                                 case 'S': case 'I': case 'F':
 476                                         rb_raise(rb_eMolbyError, "The attribute '%s' is read-only", rb_id2name(aid));
 477                                 default:
 478                                         rb_raise(rb_eMolbyError, "Internal inconsistency: unknown type field");
 479                         }
 480                 }
 481         }
 482         for (i = 0, dp = s_MDPressureAttrDefTable; dp->name != NULL; i++, dp++) {
 483                 if (dp->id == aid) {
 484                         char *pp;
 485                         MDPressureArena *pres = arena->pressure;
 486                         if (pres == NULL)
 487                                 arena->pressure = pres = pressure_new();
 488                         pp = (char *)pres + dp->offset;
 489                         switch (dp->type) {
 490                                 case 'i':
 491                                         *((Int *)pp) = NUM2INT(rb_Integer(val));
 492                                         return val;
 493                                 case 'f':
 494                                         *((Double *)pp) = NUM2DBL(rb_Float(val));
 495                                         return val;
 496                                 case 'X':
 497                                         /*  Isotropic pressure only  */
 498                                         val = rb_ary_to_ary(val);
 499                                         memset(pres->apply, 0, sizeof(Mat33));
 500                                         for (j = 0; j < 3 && j < RARRAY_LEN(val); j++)
 501                                                 pres->apply[j * 4] = NUM2DBL(rb_Float(RARRAY_PTR(val)[j]));
 502                                         return val;
 503                                 case 'Y':
 504                                         val = rb_ary_to_ary(val);
 505                                         for (j = 0; j < 8; j++) {
 506                                                 if (j < RARRAY_LEN(val))
 507                                                         pres->cell_flexibility[j] = NUM2DBL(rb_Float(RARRAY_PTR(val)[j]));
 508                                                 else pres->cell_flexibility[j] = 0.0;
 509                                         }
 510                                         return val;
 511                                 case 'S': case 'I': case 'F':
 512                                         rb_raise(rb_eMolbyError, "The attribute '%s' is read-only", rb_id2name(aid));
 513                                 default:
 514                                         rb_raise(rb_eMolbyError, "Internal inconsistency: unknown type field");
 515                         }
 516                 }
 517         }
 518         rb_raise(rb_eMolbyError, "unknown attribute name (%s)", rb_id2name(aid));
 519         return Qnil;  /*  Not reached  */
 520 }
 521
 522 static VALUE
 523 s_MDArena_SetAttr(VALUE self, VALUE val)
 524 {
 525         int i;
 526         struct s_MDArenaAttrDef *dp;
 527         ID aid = ruby_frame->last_func;
 528         for (i = 0, dp = s_MDArenaAttrDefTable; dp->name != NULL; i++, dp++) {
 529                 if (dp->sid == aid)
 530                         return s_MDArena_Set(self, *(dp->symref), val);
 531         }
 532         for (i = 0, dp = s_MDPressureAttrDefTable; dp->name != NULL; i++, dp++) {
 533                 if (dp->sid == aid)
 534                         return s_MDArena_Set(self, *(dp->symref), val);
 535         }
 536         rb_raise(rb_eMolbyError, "unknown attribute name (%s)", rb_id2name(aid));
 537         return Qnil;  /*  Not reached  */
 538 }
 539
 540 /*
 541  *  call-seq:
 542  *     to_hash -> Hash
 543  *
 544  *  Returns a (frozen) hash that contains the current value for all attribute keys.
 545  */
 546 static VALUE
 547 s_MDArena_ToHash(VALUE self)
 548 {
 549         int i;
 550         VALUE hash;
 551         struct s_MDArenaAttrDef *dp;
 552         hash = rb_hash_new();
 553         for (i = 0, dp = s_MDArenaAttrDefTable; dp->name != NULL; i++, dp++) {
 554                 VALUE attr = ID2SYM(dp->id);
 555                 rb_hash_aset(hash, attr, s_MDArena_Get(self, attr));
 556         }
 557         rb_obj_freeze(hash);
 558         return hash;
 559 }
 560
 561 /*
 562  *  call-seq:
 563  *     keys -> Array
 564  *
 565  *  Returns an array of valid attributes.
 566  */
 567 static VALUE
 568 s_MDArena_Keys(VALUE self)
 569 {
 570         int i;
 571         VALUE ary;
 572         struct s_MDArenaAttrDef *dp;
 573         ary = rb_ary_new();
 574         for (i = 0, dp = s_MDArenaAttrDefTable; dp->name != NULL; i++, dp++) {
 575                 VALUE attr = ID2SYM(dp->id);
 576                 rb_ary_push(ary, attr);
 577         }
 578         return ary;
 579 }
 580
 581 /*
 582  *  call-seq:
 583  *     print_surface_area
 584  *
 585  *  Print the surface area information to standard output. (for debug)
 586  */
 587 static VALUE
 588 s_MDArena_PrintSurfaceArea(VALUE self)
 589 {
 590         MDArena *arena;
 591         Int i, natoms;
 592         VALUE retval, outval;
 593         Data_Get_Struct(self, MDArena, arena);
 594         if (arena->sp_arena == NULL)
 595                 rb_raise(rb_eMolbyError, "surface potential is not available");
 596         natoms = arena->mol->natoms;
 597         retval = rb_str_new2("Atom     area    energy         forcex1000\n");
 598         for (i = 0; i < natoms; i++) {
 599                 char buf[256];
 600                 Vector f = arena->forces[kSurfaceIndex * natoms + i];
 601                 Double area = arena->sp_arena->atom_area[i];
 602                 Double energy = area * arena->sp_arena->atom_pot[i];
 603                 VecScaleSelf(f, INTERNAL2KCAL * 1000);
 604                 energy *= INTERNAL2KCAL;
 605                 snprintf(buf, sizeof buf, "%5d %11.5f %11.8f %11.5f %11.5f %11.5f\n",
 606                            i+1, area, energy, f.x, f.y, f.z);
 607                 rb_str_cat(retval, buf, strlen(buf));
 608         }
 609         outval = rb_gv_get("$stdout");
 610         rb_funcall(outval, rb_intern("write"), 1, retval);
 611         return self;
 612 }
 613
 614 void
 615 Init_MolbyMDTypes(void)
 616 {
 617         int i;
 618         struct s_MDArenaAttrDef *dp;
 619         char name[40];
 620
 621         /*  class MDArena  */
 622         rb_cMDArena = rb_define_class("MDArena", rb_cObject);
 623 /*      rb_define_alloc_func(rb_cMDArena, s_MDArena_Alloc);
 624     rb_define_private_method(rb_cMDArena, "initialize", s_MDArena_Initialize, 1); */
 625         rb_define_method(rb_cMDArena, "run", s_MDArena_Run, 1);
 626     rb_define_method(rb_cMDArena, "minimize", s_MDArena_Minimize, 1);
 627     rb_define_method(rb_cMDArena, "prepare", s_MDArena_Prepare, -1);
 628     rb_define_method(rb_cMDArena, "energies", s_MDArena_Energies, 0);
 629         rb_define_method(rb_cMDArena, "[]", s_MDArena_Get, 1);
 630         rb_define_method(rb_cMDArena, "[]=", s_MDArena_Set, 2);
 631         rb_define_method(rb_cMDArena, "to_hash", s_MDArena_ToHash, 0);
 632         rb_define_method(rb_cMDArena, "keys", s_MDArena_Keys, 0);
 633         rb_define_method(rb_cMDArena, "print_surface_area", s_MDArena_PrintSurfaceArea, 0);
 634
 635         /*  All setter and getter are handled with the same C function (attribute name is taken
 636             from ruby_frame)  */
 637         for (i = 0, dp = s_MDArenaAttrDefTable; dp->name != NULL; i++, dp++) {
 638                 rb_define_method(rb_cMDArena, dp->name, s_MDArena_GetAttr, 0);
 639                 strncpy(name, dp->name, 38);
 640                 strcat(name, "=");
 641                 rb_define_method(rb_cMDArena, name, s_MDArena_SetAttr, 1);
 642                 dp->id = rb_intern(dp->name);
 643                 dp->sid = rb_intern(name);
 644                 *(dp->symref) = ID2SYM(dp->id);
 645         }
 646         for (i = 0, dp = s_MDPressureAttrDefTable; dp->name != NULL; i++, dp++) {
 647                 rb_define_method(rb_cMDArena, dp->name, s_MDArena_GetAttr, 0);
 648                 strncpy(name, dp->name, 38);
 649                 strcat(name, "=");
 650                 rb_define_method(rb_cMDArena, name, s_MDArena_SetAttr, 1);
 651                 dp->id = rb_intern(dp->name);
 652                 dp->sid = rb_intern(name);
 653                 *(dp->symref) = ID2SYM(dp->id);
 654         }
 655 }