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 #if 0
 194         MDArena *arena;
 195         Molecule *mol;
 196         const char *msg;
 197         Int nangles, *angles, ndihedrals, *dihedrals, nimpropers, *impropers;
 198         Int missing = 0;
 199         Int check_only = 0;
 200         IntGroup *ig1, *ig2, *ig3;
 201         VALUE fval;
 202
 203         Data_Get_Struct(self, MDArena, arena);
 204         rb_scan_args(argc, argv, "01", &fval);
 205         if (RTEST(fval))
 206                 check_only = 1;
 207
 208         arena->is_initialized = 0;
 209         mol = arena->xmol;
 210
 211         /*  Rebuild the tables  */
 212         ig1 = ig2 = ig3 = NULL;
 213         nangles = MoleculeFindMissingAngles(mol, &angles);
 214         ndihedrals = MoleculeFindMissingDihedrals(mol, &dihedrals);
 215         nimpropers = MoleculeFindMissingImpropers(mol, &impropers);
 216         if (nangles > 0) {
 217                 ig1 = IntGroupNewWithPoints(mol->nangles, nangles, -1);
 218                 MolActionCreateAndPerform(mol, gMolActionAddAngles, nangles * 3, angles, ig1);
 219                 free(angles);
 220                 IntGroupRelease(ig1);
 221         }
 222         if (ndihedrals > 0) {
 223                 ig2 = IntGroupNewWithPoints(mol->ndihedrals, ndihedrals, -1);
 224                 MolActionCreateAndPerform(mol, gMolActionAddDihedrals, ndihedrals * 4, dihedrals, ig2);
 225                 free(dihedrals);
 226                 IntGroupRelease(ig2);
 227         }
 228         if (nimpropers > 0) {
 229                 ig3 = IntGroupNewWithPoints(mol->nimpropers, nimpropers, -1);
 230                 MolActionCreateAndPerform(mol, gMolActionAddImpropers, nimpropers * 4, impropers, ig3);
 231                 free(impropers);
 232                 IntGroupRelease(ig3);
 233         }
 234
 235         /*  Prepare parameters and internal information  */
 236         msg = md_prepare(arena, check_only);
 237
 238         /*  Some parameters are missing?  */
 239         if (msg != NULL) {
 240                 if (strstr(msg, "parameter") != NULL && strstr(msg, "missing") != NULL)
 241                         missing = 1;
 242                 else
 243                         rb_raise(rb_eMolbyError, "cannot initialize for MD: %s", msg);
 244         }
 245
 246         /*  The local parameter list is updated  */
 247         {
 248                 Int parType, idx;
 249                 if (mol->par == NULL)
 250                         mol->par = ParameterNew();
 251                 for (parType = kFirstParType; parType <= kLastParType; parType++) {
 252                         /*  Delete global and undefined parameters  */
 253                         UnionPar *up, *upbuf;
 254                         Int nparams, count;
 255                         ig1 = IntGroupNew();
 256                         for (idx = 0; (up = ParameterGetUnionParFromTypeAndIndex(mol->par, parType, idx)) != NULL; idx++) {
 257                                 if (up->bond.src != 0)
 258                                         IntGroupAdd(ig1, idx, 1);
 259                         }
 260                         if (IntGroupGetCount(ig1) > 0)
 261                                 MolActionCreateAndPerform(mol, gMolActionDeleteParameters, parType, ig1);
 262                         IntGroupRelease(ig1);
 263                         /*  Copy global and undefined parameters from arena and insert to mol->par  */
 264                         nparams = ParameterGetCountForType(arena->par, parType);
 265                         if (nparams == 0)
 266                                 continue;
 267                         upbuf = (UnionPar *)calloc(sizeof(UnionPar), nparams);
 268                         ig1 = IntGroupNew();
 269                         ig2 = IntGroupNew();
 270                         for (idx = 0; (up = ParameterGetUnionParFromTypeAndIndex(arena->par, parType, idx)) != NULL; idx++) {
 271                                 if (up->bond.src > 0)
 272                                         IntGroupAdd(ig1, idx, 1); /* Global parameter */
 273                                 else if (up->bond.src < 0)
 274                                         IntGroupAdd(ig2, idx, 1); /* Undefined parameter */
 275                         }
 276                         if ((count = IntGroupGetCount(ig1)) > 0) {
 277                                 /*  Insert global parameters (at the top)  */
 278                                 ParameterCopy(arena->par, parType, upbuf, ig1);
 279                                 ig3 = IntGroupNewWithPoints(0, count, -1);
 280                                 MolActionCreateAndPerform(mol, gMolActionAddParameters, parType, ig3, count, upbuf);
 281                                 IntGroupRelease(ig3);
 282                         }
 283                         if ((count = IntGroupGetCount(ig2)) > 0) {
 284                                 /*  Insert undefined parameters (at the bottom)  */
 285                                 ParameterCopy(arena->par, parType, upbuf, ig2);
 286                                 idx = ParameterGetCountForType(mol->par, parType);
 287                                 ig3 = IntGroupNewWithPoints(idx, count, -1);
 288                                 MolActionCreateAndPerform(mol, gMolActionAddParameters, parType, ig3, count, upbuf);
 289                                 IntGroupRelease(ig3);
 290                         }
 291                         IntGroupRelease(ig2);
 292                         IntGroupRelease(ig1);
 293                         free(upbuf);
 294                 }
 295                 mol->needsMDRebuild = 0;  /*  We know the "modified" parameters are consistent with the MDArena  */
 296         }
 297
 298         if (missing)
 299                 return Qnil;
 300         else return self;
 301 #endif
 302 }
 303
 304 /*
 305  *  call-seq:
 306  *     energies -> [total, bond, angle, dihedral, improper, vdw, electrostatic, auxiliary, surface, kinetic, net [, ewald]
 307  *
 308  *  Get the current energies.
 309  */
 310 static VALUE
 311 s_MDArena_Energies(VALUE self)
 312 {
 313         MDArena *arena;
 314         VALUE val;
 315         int i;
 316         static Int s_indices[] = {kBondIndex, kAngleIndex, kDihedralIndex, kImproperIndex, kVDWIndex, kElectrostaticIndex, kAuxiliaryIndex, kSurfaceIndex, kKineticIndex };
 317         Data_Get_Struct(self, MDArena, arena);
 318         val = rb_ary_new();
 319         rb_ary_push(val, rb_float_new(arena->total_energy * INTERNAL2KCAL));
 320         for (i = 0; i < sizeof(s_indices) / sizeof(s_indices[0]); i++)
 321                 rb_ary_push(val, rb_float_new(arena->energies[s_indices[i]] * INTERNAL2KCAL));
 322         rb_ary_push(val, rb_float_new((arena->energies[kKineticIndex] + arena->total_energy) * INTERNAL2KCAL));
 323         if (arena->use_ewald)
 324                 rb_ary_push(val, rb_float_new((arena->energies[kESCorrectionIndex] + arena->energies[kPMEIndex]) * INTERNAL2KCAL));
 325         return val;
 326 }
 327
 328 static VALUE s_LogFileSym, s_CoordFileSym, s_VelFileSym, s_ForceFileSym,
 329 s_DebugFileSym, s_DebugOutputLevelSym, s_StepSym, s_CoordOutputFreqSym,
 330 s_EnergyOutputFreqSym, s_CoordFrameSym, s_TimestepSym, s_CutoffSym,
 331 s_ElectroCutoffSym, s_PairlistDistanceSym, s_SwitchDistanceSym, s_TemperatureSym, s_TransientTempSym,
 332 s_AverageTempSym, s_AndersenFreqSym, s_AndersenCouplingSym, s_RandomSeedSym,
 333 s_DielectricSym, s_GradientConvergenceSym, s_CoordinateConvergenceSym, s_UseXplorShiftSym,
 334 s_Scale14VdwSym, s_Scale14ElectSym, s_RelocateCenterSym,
 335 s_SurfaceProbeRadiusSym, s_SurfaceTensionSym, s_SurfacePotentialFreqSym, s_UseGraphiteSym,
 336 s_AlchemicalLambdaSym, s_AlchemicalDeltaLambdaSym, s_AlchemicalEnergySym, s_MinimizeCellSym,
 337 s_UseEwaldSym, s_EwaldBetaSym, s_EwaldGridSym, s_EwaldFreqSym, s_EwaldOrderSym;
 338
 339 struct s_MDArenaAttrDef {
 340         char *name;
 341         VALUE *symref;  /*  Address of s_LogFileSym etc. */
 342         ID id;                  /*  Ruby ID of the symbol; will be set within Init_MolbyMDTypes()  */
 343         ID sid;         /*  Ruby ID of the symbol plus '='; will be set within Init_MolbyMDTypes()  */
 344         char type;      /*  s: string (const char *), i: Int, f: Double, e: Double in energy dimension (unit conversion is necessary). */
 345                                         /*  Uppercase: read-only.  */
 346         int  offset;    /*  Offset in the MDArena structure.  */
 347 };
 348 static struct s_MDArenaAttrDef s_MDArenaAttrDefTable[] = {
 349         {"log_file",          &s_LogFileSym,          0, 0, 's', offsetof(MDArena, log_result_name)},
 350         {"coord_file",        &s_CoordFileSym,        0, 0, 's', offsetof(MDArena, coord_result_name)},
 351         {"vel_file",          &s_VelFileSym,          0, 0, 's', offsetof(MDArena, vel_result_name)},
 352         {"force_file",        &s_ForceFileSym,        0, 0, 's', offsetof(MDArena, force_result_name)},
 353         {"debug_file",        &s_DebugFileSym,        0, 0, 's', offsetof(MDArena, debug_result_name)},
 354         {"debug_output_level", &s_DebugOutputLevelSym, 0, 0, 'i', offsetof(MDArena, debug_output_level)},
 355         {"step",              &s_StepSym,             0, 0, 'I', offsetof(MDArena, step)},
 356         {"coord_output_freq", &s_CoordOutputFreqSym,  0, 0, 'i', offsetof(MDArena, coord_output_freq)},
 357         {"energy_output_freq", &s_EnergyOutputFreqSym, 0, 0, 'i', offsetof(MDArena, energy_output_freq)},
 358         {"coord_frame",       &s_CoordFrameSym,       0, 0, 'I', offsetof(MDArena, coord_result_frame)},
 359         {"timestep",          &s_TimestepSym,         0, 0, 'f', offsetof(MDArena, timestep)},
 360         {"cutoff",            &s_CutoffSym,           0, 0, 'f', offsetof(MDArena, cutoff)},
 361         {"electro_cutoff",    &s_ElectroCutoffSym,    0, 0, 'f', offsetof(MDArena, electro_cutoff)},
 362         {"pairlist_distance", &s_PairlistDistanceSym, 0, 0, 'f', offsetof(MDArena, pairlist_distance)},
 363         {"switch_distance",   &s_SwitchDistanceSym,   0, 0, 'f', offsetof(MDArena, switch_distance)},
 364         {"temperature",       &s_TemperatureSym,      0, 0, 'f', offsetof(MDArena, temperature)},
 365         {"transient_temperature", &s_TransientTempSym, 0, 0, 'F', offsetof(MDArena, transient_temperature)},
 366         {"average_temperature", &s_AverageTempSym,    0, 0, 'F', offsetof(MDArena, average_temperature)},
 367         {"andersen_freq",     &s_AndersenFreqSym,     0, 0, 'i', offsetof(MDArena, andersen_thermo_freq)},
 368         {"andersen_coupling", &s_AndersenCouplingSym, 0, 0, 'f', offsetof(MDArena, andersen_thermo_coupling)},
 369         {"random_seed",       &s_RandomSeedSym,       0, 0, 'i', offsetof(MDArena, random_seed)},
 370         {"dielectric",        &s_DielectricSym,       0, 0, 'f', offsetof(MDArena, dielectric)},
 371         {"gradient_convergence", &s_GradientConvergenceSym, 0, 0, 'f', offsetof(MDArena, gradient_convergence)},
 372         {"coordinate_convergence", &s_CoordinateConvergenceSym, 0, 0, 'f', offsetof(MDArena, coordinate_convergence)},
 373         {"use_xplor_shift",   &s_UseXplorShiftSym,    0, 0, 'i', offsetof(MDArena, use_xplor_shift)},
 374         {"scale14_vdw",       &s_Scale14VdwSym,       0, 0, 'f', offsetof(MDArena, scale14_vdw)},
 375         {"scale14_elect",     &s_Scale14ElectSym,     0, 0, 'f', offsetof(MDArena, scale14_elect)},
 376         {"relocate_center",   &s_RelocateCenterSym,   0, 0, 'i', offsetof(MDArena, relocate_center)},
 377         {"surface_probe_radius", &s_SurfaceProbeRadiusSym, 0, 0, 'f', offsetof(MDArena, probe_radius)},
 378         {"surface_tension",   &s_SurfaceTensionSym,   0, 0, 'f', offsetof(MDArena, surface_tension)},
 379         {"surface_potential_freq", &s_SurfacePotentialFreqSym, 0, 0, 'i', offsetof(MDArena, surface_potential_freq)},
 380         {"use_graphite",      &s_UseGraphiteSym,      0, 0, 'i', offsetof(MDArena, use_graphite)},
 381         {"alchemical_lambda", &s_AlchemicalLambdaSym, 0, 0, 'f', offsetof(MDArena, alchem_lambda)},
 382         {"alchemical_delta_lambda", &s_AlchemicalDeltaLambdaSym, 0, 0, 'f', offsetof(MDArena, alchem_dlambda)},
 383         {"alchemical_energy", &s_AlchemicalEnergySym, 0, 0, 'E', offsetof(MDArena, alchem_energy)},
 384         {"minimize_cell",     &s_MinimizeCellSym,     0, 0, 'b', offsetof(MDArena, minimize_cell)},
 385         {"use_ewald",         &s_UseEwaldSym,         0, 0, 'i', offsetof(MDArena, use_ewald)},
 386         {"ewald_beta",        &s_EwaldBetaSym,        0, 0, 'f', offsetof(MDArena, ewald_beta)},
 387         {"ewald_grid",        &s_EwaldGridSym,        0, 0, 0,   offsetof(MDArena, ewald_grid_x)},
 388         {"ewald_freq",        &s_EwaldFreqSym,        0, 0, 'i', offsetof(MDArena, ewald_freq)},
 389         {"ewald_order",       &s_EwaldOrderSym,       0, 0, 'i', offsetof(MDArena, ewald_order)},
 390         {NULL} /* Sentinel */
 391 };
 392
 393 static VALUE s_PresFreqSym, s_PresCouplingSym, s_PresSym, s_PresCellFlexSym, s_PresFluctCellOriginSym, s_PresFluctCellOrientSym;
 394 static struct s_MDArenaAttrDef s_MDPressureAttrDefTable[] = {
 395         {"pressure_freq",     &s_PresFreqSym,         0, 0, 'i', offsetof(MDPressureArena, freq)},
 396         {"pressure_coupling", &s_PresCouplingSym,     0, 0, 'f', offsetof(MDPressureArena, coupling)},
 397         {"pressure",          &s_PresSym,             0, 0, 'X', offsetof(MDPressureArena, apply)},
 398         {"pressure_cell_flexibility", &s_PresCellFlexSym, 0, 0, 'Y', offsetof(MDPressureArena, cell_flexibility)},
 399         {"pressure_fluctuate_cell_origin", &s_PresFluctCellOriginSym, 0, 0, 'f', offsetof(MDPressureArena, fluctuate_cell_origin)},
 400         {"pressure_fluctuate_cell_orientation", &s_PresFluctCellOrientSym, 0, 0, 'f', offsetof(MDPressureArena, fluctuate_cell_orientation)},
 401         {NULL} /* Sentinel */
 402 };
 403
 404 /*
 405  *  call-seq:
 406  *     self[attr]
 407  *
 408  *  Get the attribute value. The attribute values also can be accessed by self.attribute_name.
 409  */
 410 static VALUE
 411 s_MDArena_Get(VALUE self, VALUE attr)
 412 {
 413         MDArena *arena;
 414         int i;
 415         struct s_MDArenaAttrDef *dp;
 416         ID aid = rb_to_id(attr);
 417         Data_Get_Struct(self, MDArena, arena);
 418         if (aid == SYM2ID(s_EwaldGridSym)) {
 419                 /*  Array of three grid values  */
 420                 return rb_ary_new3(3, INT2NUM(arena->ewald_grid_x), INT2NUM(arena->ewald_grid_y), INT2NUM(arena->ewald_grid_z));
 421         }
 422         for (i = 0, dp = s_MDArenaAttrDefTable; dp->name != NULL; i++, dp++) {
 423                 if (dp->id == aid) {
 424                         char *p = (char *)arena + dp->offset;
 425                         switch (dp->type) {
 426                                 case 's':
 427                                 case 'S': {
 428                                         const char *cp = *((const char **)p);
 429                                         if (cp == NULL)
 430                                                 return Qnil;
 431                                         else
 432                                                 return Ruby_NewFileStringValue(cp);
 433                                 }
 434                                 case 'b':
 435                                 case 'B':
 436                                         return INT2NUM((Int)(*((Byte *)p)));
 437                                 case 'i':
 438                                 case 'I':
 439                                         return INT2NUM(*((Int *)p));
 440                                 case 'f':
 441                                 case 'F':
 442                                         return rb_float_new(*((Double *)p));
 443                                 case 'e':
 444                                 case 'E':
 445                                         return rb_float_new(*((Double *)p) * INTERNAL2KCAL);
 446                                 default:
 447                                         rb_raise(rb_eMolbyError, "Internal inconsistency: unknown type field");
 448                         }
 449                 }
 450         }
 451         for (i = 0, dp = s_MDPressureAttrDefTable; dp->name != NULL; i++, dp++) {
 452                 if (dp->id == aid) {
 453                         char *pp;
 454                         MDPressureArena *pres = arena->pressure;
 455                         if (pres == NULL)
 456                                 return Qnil;
 457                         pp = (char *)pres + dp->offset;
 458                         switch (dp->type) {
 459                                 case 'i':
 460                                 case 'I':
 461                                         return INT2NUM(*((Int *)pp));
 462                                 case 'b':
 463                                 case 'B':
 464                                         return INT2NUM((Int)(*((Byte *)pp)));
 465                                 case 'f':
 466                                 case 'F':
 467                                         return rb_float_new(*((Double *)pp));
 468                                 case 'e':
 469                                 case 'E':
 470                                         return rb_float_new(*((Double *)pp) * INTERNAL2KCAL);
 471                                 case 'X':
 472                                         /*  Isotropic pressure only  */
 473                                         return rb_ary_new3(3, rb_float_new(pres->apply[0]), rb_float_new(pres->apply[4]), rb_float_new(pres->apply[8]));
 474                                 case 'Y': {
 475                                         VALUE aval = rb_ary_new();
 476                                         int j;
 477                                         for (j = 0; j < 8; j++)
 478                                                 rb_ary_push(aval, rb_float_new(pres->cell_flexibility[j]));
 479                                         return aval;
 480                                 }
 481                                 default:
 482                                         rb_raise(rb_eMolbyError, "Internal inconsistency: unknown type field");
 483                         }
 484                 }
 485         }
 486         rb_raise(rb_eMolbyError, "unknown attribute name (%s)", rb_id2name(aid));
 487         return Qnil;  /*  Not reached  */
 488 }
 489
 490 static VALUE
 491 s_MDArena_GetAttr(VALUE self)
 492 {
 493         ID aid = rb_frame_this_func();
 494         return s_MDArena_Get(self, ID2SYM(aid));
 495 }
 496
 497 /*
 498  *  call-seq:
 499  *     self[attr] = value
 500  *
 501  *  Set the attribute value. The attributes can also be modified by self.attribute_name = value.
 502  */
 503 static VALUE
 504 s_MDArena_Set(VALUE self, VALUE attr, VALUE val)
 505 {
 506         MDArena *arena;
 507         int i, j;
 508         struct s_MDArenaAttrDef *dp;
 509         ID aid = rb_to_id(attr);
 510         attr = ID2SYM(aid);  /*  May be used later  */
 511         Data_Get_Struct(self, MDArena, arena);
 512         if (aid == SYM2ID(s_EwaldGridSym)) {
 513                 if (rb_obj_is_kind_of(val, rb_cNumeric)) {
 514                         i = NUM2INT(rb_Integer(val));
 515                         if (i <= 0)
 516                                 rb_raise(rb_eMolbyError, "The ewald grid must be positive integer");
 517                         arena->ewald_grid_x = arena->ewald_grid_y = arena->ewald_grid_z = i;
 518                 } else {
 519                         int ival[3];
 520                         val = rb_ary_to_ary(val);
 521                         j = RARRAY_LEN(val);
 522                         if (j < 3)
 523                                 rb_raise(rb_eMolbyError, "The ewald grid must be an integer or an array of three integers");
 524                         for (i = 0; i < 3; i++) {
 525                                 ival[i] = NUM2INT(rb_Integer(RARRAY_PTR(val)[i]));
 526                                 if (ival[i] <= 0)
 527                                         rb_raise(rb_eMolbyError, "The ewald grid must be positive integer");
 528                         }
 529                         arena->ewald_grid_x = ival[0];
 530                         arena->ewald_grid_y = ival[1];
 531                         arena->ewald_grid_z = ival[2];
 532                 }
 533                 return val;
 534         }
 535         for (i = 0, dp = s_MDArenaAttrDefTable; dp->name != NULL; i++, dp++) {
 536                 if (dp->id == aid) {
 537                         char *p = (char *)arena + dp->offset;
 538                         switch (dp->type) {
 539                                 case 's': {
 540                                         const char *cp = (val == Qnil ? NULL : (const char *)strdup(FileStringValuePtr(val)));
 541                                         const char **cpp = (const char **)p;
 542                                         FILE **fpp;
 543                                         if (*cpp == cp || (*cpp != NULL && cp != NULL && strcmp(*cpp, cp) == 0))
 544                                                 return val;  /*  No need to change  */
 545                                         if (*cpp != NULL)
 546                                                 free((void *)*cpp);
 547                                         if (cp != NULL && cp[0] == 0) {
 548                                                 free((void *)cp);
 549                                                 cp = NULL;
 550                                         }
 551                                         /*  Close the corresponding FILE if necessary  */
 552                                         if (attr == s_LogFileSym)
 553                                                 fpp = &(arena->log_result);
 554                                         else if (attr == s_CoordFileSym)
 555                                                 fpp = &(arena->coord_result);
 556                                         else if (attr == s_VelFileSym)
 557                                                 fpp = &(arena->vel_result);
 558                                         else if (attr == s_ForceFileSym)
 559                                                 fpp = &(arena->force_result);
 560                                         else if (attr == s_DebugFileSym)
 561                                                 fpp = &(arena->debug_result);
 562                                         else fpp = NULL;
 563                                         if (fpp != NULL && *fpp != NULL) {
 564                                                 fclose(*fpp);
 565                                                 *fpp = NULL;
 566                                         }
 567                                         *cpp = cp;
 568                                         return val;
 569                                 }
 570                                 case 'i':
 571                                         *((Int *)p) = NUM2INT(rb_Integer(val));
 572                                         return val;
 573                                 case 'b':
 574                                         *((Byte *)p) = NUM2INT(rb_Integer(val));
 575                                         return val;
 576                                 case 'f':
 577                                         *((Double *)p) = NUM2DBL(rb_Float(val));
 578                                         return val;
 579                                 case 'e':
 580                                         *((Double *)p) = NUM2DBL(rb_Float(val) * KCAL2INTERNAL);
 581                                         return val;
 582                                 case 'S': case 'I': case 'B': case 'F': case 'E':
 583                                         rb_raise(rb_eMolbyError, "The attribute '%s' is read-only", rb_id2name(aid));
 584                                 default:
 585                                         rb_raise(rb_eMolbyError, "Internal inconsistency: unknown type field");
 586                         }
 587                 }
 588         }
 589         for (i = 0, dp = s_MDPressureAttrDefTable; dp->name != NULL; i++, dp++) {
 590                 if (dp->id == aid) {
 591                         char *pp;
 592                         MDPressureArena *pres = arena->pressure;
 593                         if (pres == NULL)
 594                                 arena->pressure = pres = pressure_new();
 595                         pp = (char *)pres + dp->offset;
 596                         switch (dp->type) {
 597                                 case 'i':
 598                                         *((Int *)pp) = NUM2INT(rb_Integer(val));
 599                                         return val;
 600                                 case 'b':
 601                                         *((Byte *)pp) = NUM2INT(rb_Integer(val));
 602                                         return val;
 603                                 case 'f':
 604                                         *((Double *)pp) = NUM2DBL(rb_Float(val));
 605                                         return val;
 606                                 case 'e':
 607                                         *((Double *)pp) = NUM2DBL(rb_Float(val) * KCAL2INTERNAL);
 608                                         return val;
 609                                 case 'X':
 610                                         /*  Isotropic pressure only  */
 611                                         val = rb_ary_to_ary(val);
 612                                         memset(pres->apply, 0, sizeof(Mat33));
 613                                         for (j = 0; j < 3 && j < RARRAY_LEN(val); j++)
 614                                                 pres->apply[j * 4] = NUM2DBL(rb_Float(RARRAY_PTR(val)[j]));
 615                                         return val;
 616                                 case 'Y':
 617                                         val = rb_ary_to_ary(val);
 618                                         for (j = 0; j < 8; j++) {
 619                                                 if (j < RARRAY_LEN(val))
 620                                                         pres->cell_flexibility[j] = NUM2DBL(rb_Float(RARRAY_PTR(val)[j]));
 621                                                 else pres->cell_flexibility[j] = 0.0;
 622                                         }
 623                                         return val;
 624                                 case 'S': case 'I': case 'B': case 'F': case 'E':
 625                                         rb_raise(rb_eMolbyError, "The attribute '%s' is read-only", rb_id2name(aid));
 626                                 default:
 627                                         rb_raise(rb_eMolbyError, "Internal inconsistency: unknown type field");
 628                         }
 629                 }
 630         }
 631         rb_raise(rb_eMolbyError, "unknown attribute name (%s)", rb_id2name(aid));
 632         return Qnil;  /*  Not reached  */
 633 }
 634
 635 static VALUE
 636 s_MDArena_SetAttr(VALUE self, VALUE val)
 637 {
 638         int i;
 639         struct s_MDArenaAttrDef *dp;
 640 //      ID aid = ruby_frame->last_func;
 641         ID aid = rb_frame_this_func();
 642         for (i = 0, dp = s_MDArenaAttrDefTable; dp->name != NULL; i++, dp++) {
 643                 if (dp->sid == aid)
 644                         return s_MDArena_Set(self, *(dp->symref), val);
 645         }
 646         for (i = 0, dp = s_MDPressureAttrDefTable; dp->name != NULL; i++, dp++) {
 647                 if (dp->sid == aid)
 648                         return s_MDArena_Set(self, *(dp->symref), val);
 649         }
 650         rb_raise(rb_eMolbyError, "unknown attribute name (%s)", rb_id2name(aid));
 651         return Qnil;  /*  Not reached  */
 652 }
 653
 654 /*
 655  *  call-seq:
 656  *     to_hash -> Hash
 657  *
 658  *  Returns a (frozen) hash that contains the current value for all attribute keys.
 659  */
 660 static VALUE
 661 s_MDArena_ToHash(VALUE self)
 662 {
 663         int i;
 664         VALUE hash;
 665         struct s_MDArenaAttrDef *dp;
 666         hash = rb_hash_new();
 667         for (i = 0, dp = s_MDArenaAttrDefTable; dp->name != NULL; i++, dp++) {
 668                 VALUE attr = ID2SYM(dp->id);
 669                 rb_hash_aset(hash, attr, s_MDArena_Get(self, attr));
 670         }
 671         rb_obj_freeze(hash);
 672         return hash;
 673 }
 674
 675 /*
 676  *  call-seq:
 677  *     set_alchemical_perturbation(group1, group2) -> [group1, group2]
 678  *
 679  *  Set vanishing and appearing atom groups for alchemical perturbation.
 680  */
 681 static VALUE
 682 s_MDArena_SetAlchemicalPerturbation(VALUE self, VALUE gval1, VALUE gval2)
 683 {
 684         IntGroup *ig1, *ig2;
 685         MDArena *arena;
 686         char *flags;
 687         int i, n;
 688         Data_Get_Struct(self, MDArena, arena);
 689         if (gval1 == Qnil && gval2 == Qnil) {
 690                 md_set_alchemical_flags(arena, 0, NULL);
 691                 return Qnil;
 692         }
 693         if (arena->mol == NULL)
 694                 rb_raise(rb_eMolbyError, "Molecule is not set");
 695         n = arena->xmol->natoms;
 696         flags = (char *)calloc(1, n);
 697         ig1 = (gval1 == Qnil ? NULL : IntGroupFromValue(gval1));
 698         ig2 = (gval2 == Qnil ? NULL : IntGroupFromValue(gval2));
 699         for (i = 0; i < n; i++) {
 700                 if (ig1 != NULL && IntGroupLookupPoint(ig1, i) >= 0)
 701                         flags[i] = 1;
 702                 if (ig2 != NULL && IntGroupLookupPoint(ig2, i) >= 0) {
 703                         if (flags[i] == 1)
 704                                 rb_raise(rb_eMolbyError, "duplicate atom (%d) in vanishing and appearing groups", i);
 705                         flags[i] = 2;
 706                 }
 707         }
 708         if (md_set_alchemical_flags(arena, n, flags) != 0)
 709                 rb_raise(rb_eMolbyError, "cannot set alchemical flags");
 710         free(flags);
 711         if (ig1 != NULL) {
 712                 gval1 = ValueFromIntGroup(ig1);
 713                 IntGroupRelease(ig1);
 714         }
 715         if (ig2 != NULL) {
 716                 gval2 = ValueFromIntGroup(ig2);
 717                 IntGroupRelease(ig2);
 718         }
 719         return rb_ary_new3(2, gval1, gval2);
 720 }
 721
 722 /*
 723  *  call-seq:
 724  *     get_alchemical_perturbation -> [group1, group2] or nil
 725  *
 726  *  If alchemical perturbation is enabled, get the current vanishing and appearing atom groups.
 727  *  Otherwise, return nil.
 728  */
 729 static VALUE
 730 s_MDArena_GetAlchemicalPerturbation(VALUE self)
 731 {
 732         IntGroup *ig1, *ig2;
 733         VALUE gval1, gval2;
 734         MDArena *arena;
 735         int i;
 736         Data_Get_Struct(self, MDArena, arena);
 737         if (arena->nalchem_flags == 0)
 738                 return Qnil;
 739         if (arena->mol == NULL)
 740                 rb_raise(rb_eMolbyError, "Molecule is not set");
 741         ig1 = IntGroupNew();
 742         ig2 = IntGroupNew();
 743         for (i = 0; i < arena->nalchem_flags; i++) {
 744                 if (arena->alchem_flags[i] == 1)
 745                         IntGroupAdd(ig1, i, 1);
 746                 else if (arena->alchem_flags[i] == 2)
 747                         IntGroupAdd(ig2, i, 1);
 748         }
 749         gval1 = ValueFromIntGroup(ig1);
 750         gval2 = ValueFromIntGroup(ig2);
 751         IntGroupRelease(ig1);
 752         IntGroupRelease(ig2);
 753         return rb_ary_new3(2, gval1, gval2);
 754 }
 755
 756 /*
 757  *  call-seq:
 758  *    set_external_forces(ary) -> self
 759  *
 760  *  Set external forces. Ary should be an array of objects that can be converted to Vector3D.
 761  */
 762 static VALUE
 763 s_MDArena_SetExternalForces(VALUE self, VALUE aval)
 764 {
 765         Vector *vp;
 766         MDArena *arena;
 767         int i, n;
 768         Data_Get_Struct(self, MDArena, arena);
 769         if (arena->mol == NULL)
 770                 rb_raise(rb_eMolbyError, "Molecule is not set");
 771         if (aval == Qnil) {
 772                 md_set_external_forces(arena, 0, NULL);
 773                 return self;
 774         }
 775         aval = rb_ary_to_ary(aval);
 776         n = RARRAY_LEN(aval);
 777         if (n == 0) {
 778                 md_set_external_forces(arena, 0, NULL);
 779                 return self;
 780         }
 781         vp = (Vector *)calloc(sizeof(Vector), n);
 782         for (i = 0; i < n; i++) {
 783                 VectorFromValue(RARRAY_PTR(aval)[i], vp + i);
 784                 VecScaleSelf(vp[i], KCAL2INTERNAL);
 785         }
 786         md_set_external_forces(arena, n, vp);
 787         free(vp);
 788         return self;
 789 }
 790
 791 /*
 792  *  call-seq:
 793  *    get_external_force(index) -> Vector3D or nil
 794  *
 795  *  Get the current external force for the atom. If the external force is not set, nil is returned.
 796  */
 797 static VALUE
 798 s_MDArena_GetExternalForce(VALUE self, VALUE ival)
 799 {
 800         int i;
 801         VALUE vval;
 802         Vector v;
 803         MDArena *arena;
 804         Data_Get_Struct(self, MDArena, arena);
 805         if (arena->mol == NULL)
 806                 rb_raise(rb_eMolbyError, "Molecule is not set");
 807         i = NUM2INT(rb_Integer(ival));
 808         if (i < 0 || i >= arena->nexforces)
 809                 return Qnil;
 810         v = arena->exforces[i];
 811         VecScaleSelf(v, INTERNAL2KCAL);
 812         vval = ValueFromVector(&v);
 813         return vval;
 814 }
 815
 816 /*
 817  *  call-seq:
 818  *     init_velocities([temperature]) -> self
 819  *
 820  *  Give random (Boltzmann-weighted) velocities to all atoms. If temperature is given,
 821  *  it is set before giving velocities.
 822  */
 823 static VALUE
 824 s_MDArena_InitVelocities(int argc, VALUE *argv, VALUE self)
 825 {
 826         MDArena *arena;
 827         VALUE tval;
 828         Data_Get_Struct(self, MDArena, arena);
 829         rb_scan_args(argc, argv, "01", &tval);
 830         if (tval != Qnil)
 831                 s_MDArena_Set(self, s_TemperatureSym, tval);
 832         md_init_velocities(arena);
 833         return self;
 834 }
 835
 836 /*
 837  *  call-seq:
 838  *     scale_velocities([temperature]) -> self
 839  *
 840  *  Scale atom velocities to match the current temperature. If temperature is given,
 841  *  it is set before giving velocities.
 842  */
 843 static VALUE
 844 s_MDArena_ScaleVelocities(int argc, VALUE *argv, VALUE self)
 845 {
 846         MDArena *arena;
 847         VALUE tval;
 848         Data_Get_Struct(self, MDArena, arena);
 849         rb_scan_args(argc, argv, "01", &tval);
 850         if (tval != Qnil)
 851                 s_MDArena_Set(self, s_TemperatureSym, tval);
 852         md_scale_velocities(arena);
 853         return self;
 854 }
 855
 856 /*
 857  *  call-seq:
 858  *     keys -> Array
 859  *
 860  *  Returns an array of valid attributes.
 861  */
 862 static VALUE
 863 s_MDArena_Keys(VALUE self)
 864 {
 865         int i;
 866         VALUE ary;
 867         struct s_MDArenaAttrDef *dp;
 868         ary = rb_ary_new();
 869         for (i = 0, dp = s_MDArenaAttrDefTable; dp->name != NULL; i++, dp++) {
 870                 VALUE attr = ID2SYM(dp->id);
 871                 rb_ary_push(ary, attr);
 872         }
 873         return ary;
 874 }
 875
 876 /*
 877  *  call-seq:
 878  *     print_surface_area
 879  *
 880  *  Print the surface area information to standard output. (for debug)
 881  */
 882 static VALUE
 883 s_MDArena_PrintSurfaceArea(VALUE self)
 884 {
 885         MDArena *arena;
 886         Int i, natoms;
 887         VALUE retval, outval;
 888         Data_Get_Struct(self, MDArena, arena);
 889         if (arena->sp_arena == NULL)
 890                 rb_raise(rb_eMolbyError, "surface potential is not available");
 891         natoms = arena->mol->natoms;
 892         retval = rb_str_new2("Atom     area    energy         forcex1000\n");
 893         for (i = 0; i < natoms; i++) {
 894                 char buf[256];
 895                 Vector f = arena->forces[kSurfaceIndex * natoms + i];
 896                 Double area = arena->sp_arena->atom_area[i];
 897                 Double energy = area * arena->sp_arena->atom_pot[i];
 898                 VecScaleSelf(f, INTERNAL2KCAL * 1000);
 899                 energy *= INTERNAL2KCAL;
 900                 snprintf(buf, sizeof buf, "%5d %11.5f %11.8f %11.5f %11.5f %11.5f\n",
 901                            i+1, area, energy, f.x, f.y, f.z);
 902                 rb_str_cat(retval, buf, strlen(buf));
 903         }
 904         outval = rb_gv_get("$stdout");
 905         rb_funcall(outval, rb_intern("write"), 1, retval);
 906         return self;
 907 }
 908
 909 /*
 910  *  call-seq:
 911  *     bond_par(idx) -> ParameterRef
 912  *
 913  *  Returns a parameter that is used for the idx-th bond.
 914  */
 915 static VALUE
 916 s_MDArena_BondPar(VALUE self, VALUE val)
 917 {
 918         MDArena *arena;
 919         Int i, j;
 920         Data_Get_Struct(self, MDArena, arena);
 921         i = NUM2INT(rb_Integer(val));
 922         if (arena->xmol->needsMDRebuild || arena->par == NULL)
 923                 md_prepare(arena, 1);
 924         if (i < 0 || i >= arena->xmol->nbonds) {
 925                 rb_raise(rb_eMolbyError, "bond index (%d) out of range (0...%d)", i, arena->xmol->nbonds);
 926         }
 927         j = arena->bond_par_i[i];
 928         if (j < 0 || j >= arena->par->nbondPars)
 929                 return Qnil;  /*  No parameter assigned  */
 930         return ValueFromMoleculeWithParameterTypeAndIndex(arena->xmol, kBondParType, -j - 1);
 931 }
 932
 933 /*
 934  *  call-seq:
 935  *     angle_par(idx) -> ParameterRef
 936  *
 937  *  Returns a parameter that is used for the idx-th angle.
 938  */
 939 static VALUE
 940 s_MDArena_AnglePar(VALUE self, VALUE val)
 941 {
 942         MDArena *arena;
 943         Int i, j;
 944         Data_Get_Struct(self, MDArena, arena);
 945         i = NUM2INT(rb_Integer(val));
 946         if (arena->xmol->needsMDRebuild || arena->par == NULL)
 947                 md_prepare(arena, 1);
 948         if (i < 0 || i >= arena->xmol->nangles) {
 949                 rb_raise(rb_eMolbyError, "angle index (%d) out of range (0...%d)", i, arena->xmol->nangles);
 950         }
 951         j = arena->angle_par_i[i];
 952         if (j < 0 || j >= arena->par->nanglePars)
 953                 return Qnil;  /*  No parameter assigned  */
 954         return ValueFromMoleculeWithParameterTypeAndIndex(arena->xmol, kAngleParType, -j - 1);
 955 }
 956
 957 /*
 958  *  call-seq:
 959  *     dihedral_par(idx) -> ParameterRef
 960  *
 961  *  Returns a parameter that is used for the idx-th dihedral.
 962  */
 963 static VALUE
 964 s_MDArena_DihedralPar(VALUE self, VALUE val)
 965 {
 966         MDArena *arena;
 967         Int i, j;
 968         Data_Get_Struct(self, MDArena, arena);
 969         i = NUM2INT(rb_Integer(val));
 970         if (arena->xmol->needsMDRebuild || arena->par == NULL)
 971                 md_prepare(arena, 1);
 972         if (i < 0 || i >= arena->xmol->ndihedrals) {
 973                 rb_raise(rb_eMolbyError, "dihedral index (%d) out of range (0...%d)", i, arena->xmol->ndihedrals);
 974         }
 975         j = arena->dihedral_par_i[i];
 976         if (j < 0 || j >= arena->par->ndihedralPars)
 977                 return Qnil;  /*  No parameter assigned  */
 978         return ValueFromMoleculeWithParameterTypeAndIndex(arena->xmol, kDihedralParType, -j - 1);
 979 }
 980
 981 /*
 982  *  call-seq:
 983  *     improper_par(idx) -> ParameterRef
 984  *
 985  *  Returns a parameter that is used for the idx-th improper.
 986  */
 987 static VALUE
 988 s_MDArena_ImproperPar(VALUE self, VALUE val)
 989 {
 990         MDArena *arena;
 991         Int i, j;
 992         Data_Get_Struct(self, MDArena, arena);
 993         i = NUM2INT(rb_Integer(val));
 994         if (arena->xmol->needsMDRebuild || arena->par == NULL)
 995                 md_prepare(arena, 1);
 996         if (i < 0 || i >= arena->xmol->nimpropers) {
 997                 rb_raise(rb_eMolbyError, "improper index (%d) out of range (0...%d)", i, arena->xmol->nimpropers);
 998         }
 999         j = arena->improper_par_i[i];
1000         if (j < 0 || j >= arena->par->nimproperPars)
1001                 return Qnil;  /*  No parameter assigned  */
1002         return ValueFromMoleculeWithParameterTypeAndIndex(arena->xmol, kImproperParType, -j - 1);
1003 }
1004
1005 /*
1006  *  call-seq:
1007  *     vdw_par(idx) -> ParameterRef
1008  *
1009  *  Returns a vdw parameter that is used for the idx-th atom.
1010  */
1011 static VALUE
1012 s_MDArena_VdwPar(VALUE self, VALUE val)
1013 {
1014         MDArena *arena;
1015         Int i, j;
1016         Data_Get_Struct(self, MDArena, arena);
1017         i = NUM2INT(rb_Integer(val));
1018         if (arena->xmol->needsMDRebuild || arena->par == NULL)
1019                 md_prepare(arena, 1);
1020         if (i < 0 || i >= arena->xmol->natoms) {
1021                 rb_raise(rb_eMolbyError, "atom index (%d) out of range (0...%d)", i, arena->xmol->natoms);
1022         }
1023         j = arena->vdw_par_i[i];
1024         if (j < 0 || j >= arena->par->nvdwPars)
1025                 return Qnil;  /*  No parameter assigned  */
1026         return ValueFromMoleculeWithParameterTypeAndIndex(arena->xmol, kVdwParType, -j - 1);
1027 }
1028
1029 static VALUE
1030 s_MDArena_testPME(int argc, VALUE *argv, VALUE self)
1031 {
1032         extern int pme_test(MDArena *);
1033         MDArena *arena;
1034         if (rb_obj_is_kind_of(self, rb_cMDArena)) {
1035                 Data_Get_Struct(self, MDArena, arena);
1036                 if (argc == 0 || RTEST(argv[0])) {
1037                         arena->use_ewald = 1;
1038                         arena->xmol->needsMDRebuild = 1;
1039                         if (argc > 0 && rb_obj_is_kind_of(argv[0], rb_cNumeric)) {
1040                                 Int n = NUM2INT(rb_Integer(argv[0]));
1041                                 if (n > 1) {
1042                                         arena->ewald_grid_x = n;
1043                                         arena->ewald_grid_y = n;
1044                                         arena->ewald_grid_z = n;
1045                                 } else {
1046                                         arena->use_ewald = 2;  /*  Direct Ewald  */
1047                                 }
1048                         }
1049                 } else {
1050                         arena->use_ewald = 0;
1051                         arena->xmol->needsMDRebuild = 1;
1052                 }
1053                 s_MDArena_Prepare(0, NULL, self);
1054         } else {
1055                 rb_raise(rb_eMolbyError, "class method MDArena.test_pme is not supported now");
1056         }
1057         return self;
1058 }
1059
1060 void
1061 Init_MolbyMDTypes(void)
1062 {
1063         int i;
1064         struct s_MDArenaAttrDef *dp;
1065         char name[40];
1066
1067         /*  class MDArena  */
1068         rb_cMDArena = rb_define_class("MDArena", rb_cObject);
1069 /*      rb_define_alloc_func(rb_cMDArena, s_MDArena_Alloc);
1070     rb_define_private_method(rb_cMDArena, "initialize", s_MDArena_Initialize, 1); */
1071         rb_define_method(rb_cMDArena, "run", s_MDArena_Run, 1);
1072     rb_define_method(rb_cMDArena, "minimize", s_MDArena_Minimize, 1);
1073     rb_define_method(rb_cMDArena, "prepare", s_MDArena_Prepare, -1);
1074     rb_define_method(rb_cMDArena, "energies", s_MDArena_Energies, 0);
1075         rb_define_method(rb_cMDArena, "[]", s_MDArena_Get, 1);
1076         rb_define_method(rb_cMDArena, "[]=", s_MDArena_Set, 2);
1077         rb_define_method(rb_cMDArena, "to_hash", s_MDArena_ToHash, 0);
1078         rb_define_method(rb_cMDArena, "keys", s_MDArena_Keys, 0);
1079         rb_define_method(rb_cMDArena, "set_alchemical_perturbation", s_MDArena_SetAlchemicalPerturbation, 2);
1080         rb_define_method(rb_cMDArena, "get_alchemical_perturbation", s_MDArena_GetAlchemicalPerturbation, 0);
1081         rb_define_method(rb_cMDArena, "set_external_forces", s_MDArena_SetExternalForces, 1);
1082         rb_define_method(rb_cMDArena, "get_external_force", s_MDArena_GetExternalForce, 1);
1083         rb_define_method(rb_cMDArena, "init_velocities", s_MDArena_InitVelocities, -1);
1084         rb_define_method(rb_cMDArena, "scale_velocities", s_MDArena_ScaleVelocities, -1);
1085         rb_define_method(rb_cMDArena, "print_surface_area", s_MDArena_PrintSurfaceArea, 0);
1086
1087         rb_define_method(rb_cMDArena, "bond_par", s_MDArena_BondPar, 1);
1088         rb_define_method(rb_cMDArena, "angle_par", s_MDArena_AnglePar, 1);
1089         rb_define_method(rb_cMDArena, "dihedral_par", s_MDArena_DihedralPar, 1);
1090         rb_define_method(rb_cMDArena, "improper_par", s_MDArena_ImproperPar, 1);
1091         rb_define_method(rb_cMDArena, "vdw_par", s_MDArena_VdwPar, 1);
1092
1093         rb_define_method(rb_cMDArena, "pme_test", s_MDArena_testPME, -1);
1094 //      rb_define_singleton_method(rb_cMDArena, "pme_test", s_MDArena_testPME, -1);
1095
1096         /*  All setter and getter are handled with the same C function (attribute name is taken
1097             from ruby_frame)  */
1098         for (i = 0, dp = s_MDArenaAttrDefTable; dp->name != NULL; i++, dp++) {
1099                 rb_define_method(rb_cMDArena, dp->name, s_MDArena_GetAttr, 0);
1100                 strncpy(name, dp->name, 38);
1101                 name[38] = 0;
1102                 strcat(name, "=");
1103                 rb_define_method(rb_cMDArena, name, s_MDArena_SetAttr, 1);
1104                 dp->id = rb_intern(dp->name);
1105                 dp->sid = rb_intern(name);
1106                 *(dp->symref) = ID2SYM(dp->id);
1107         }
1108         for (i = 0, dp = s_MDPressureAttrDefTable; dp->name != NULL; i++, dp++) {
1109                 rb_define_method(rb_cMDArena, dp->name, s_MDArena_GetAttr, 0);
1110                 strncpy(name, dp->name, 38);
1111                 name[38] = 0;
1112                 strcat(name, "=");
1113                 rb_define_method(rb_cMDArena, name, s_MDArena_SetAttr, 1);
1114                 dp->id = rb_intern(dp->name);
1115                 dp->sid = rb_intern(name);
1116                 *(dp->symref) = ID2SYM(dp->id);
1117         }
1118 }