MolLib/MolAction.c

   1 /*
   2  *  MolAction.c
   3  *
   4  *  Created by Toshi Nagata on 07/06/23.
   5  *  Copyright 2007-2008 Toshi Nagata. All rights reserved.
   6  *
   7  This program is free software; you can redistribute it and/or modify
   8  it under the terms of the GNU General Public License as published by
   9  the Free Software Foundation version 2 of the License.
  10
  11  This program is distributed in the hope that it will be useful,
  12  but WITHOUT ANY WARRANTY; without even the implied warranty of
  13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14  GNU General Public License for more details.
  15  */
  16
  17 #include "Ruby_bind/Molby.h"  /*  Support Ruby binding  */
  18
  19 #include "IntGroup.h"
  20 #include "Molecule.h"
  21 #include "MolAction.h"
  22 #include "MD/MDCore.h"
  23
  24 #include <stdarg.h>
  25 #include <ctype.h>
  26
  27 const char *gMolActionNone            = "none";
  28 const char *gMolActionAddAnAtom       = "addAtom:Ai;i";
  29 const char *gMolActionDeleteAnAtom    = "deleteAtom:i";
  30 const char *gMolActionMergeMolecule   = "mergeMol:MG";
  31 const char *gMolActionUnmergeMolecule = "unmergeMol:G";
  32 const char *gMolActionAddBonds        = "addBonds:I";
  33 const char *gMolActionDeleteBonds     = "deleteBonds:I";
  34 const char *gMolActionAddAngles       = "addAngles:IG";
  35 const char *gMolActionDeleteAngles    = "deleteAngles:G";
  36 const char *gMolActionAddDihedrals    = "addDihedrals:IG";
  37 const char *gMolActionDeleteDihedrals = "deleteDihedrals:G";
  38 const char *gMolActionAddImpropers    = "addImpropers:IG";
  39 const char *gMolActionDeleteImpropers = "deleteImpropers:G";
  40 const char *gMolActionTranslateAtoms  = "translateAtoms:vG";
  41 const char *gMolActionRotateAtoms     = "rotate:vdvG";
  42 const char *gMolActionTransformAtoms  = "transform:tG";
  43 const char *gMolActionSetAtomPositions = "atomPosition:GV";
  44 const char *gMolActionSetAtomVelocities = "atomVelocity:GV";
  45 const char *gMolActionSetAtomForces   = "atomForce:GV";
  46 const char *gMolActionInsertFrames    = "insertFrames:GVV";
  47 const char *gMolActionRemoveFrames    = "removeFrames:G";
  48 const char *gMolActionSetSelection    = "selection:G";
  49 const char *gMolActionChangeResidueNumber = "changeResSeq:Gi";
  50 const char *gMolActionChangeResidueNumberForUndo = "changeResSeqForUndo:GIi";
  51 const char *gMolActionChangeResidueNames = "changeResNames:IC";
  52 const char *gMolActionOffsetResidueNumbers = "offsetResSeq:Gii";
  53 const char *gMolActionChangeNumberOfResidues = "changeNres:i";
  54 const char *gMolActionRenumberAtoms    = "renumberAtoms:I";
  55 const char *gMolActionExpandBySymmetry = "expandSym:Giiii;I";
  56 const char *gMolActionDeleteSymmetryOperation = "deleteSymop";
  57 const char *gMolActionAddSymmetryOperation = "addSymop:t";
  58 const char *gMolActionSetCell         = "setCell:Di";
  59 const char *gMolActionSetBox          = "setBox:vvvvi";
  60 const char *gMolActionClearBox        = "clearBox";
  61 const char *gMolActionSetBoxForFrames = "setBoxForFrames:V";
  62 const char *gMolActionSetCellFlexibility = "setCellFlexibility:i";
  63 const char *gMolActionAddParameters   = "addParameters:iGU";
  64 const char *gMolActionDeleteParameters = "deleteParameters:iG";
  65 const char *gMolActionCartesianToXtal  = "cartesianToXtal";
  66 const char *gMolActionXtalToCartesian  = "xtalToCartesian";
  67 const char *gMolActionAmendBySymmetry = "amendBySymmetry:G;G";
  68
  69 /*  A Ruby array to retain objects used in MolActionArg  */
  70 static VALUE sMolActionArgValues = Qfalse;
  71
  72 /*  Action arguments  */
  73 /*  (Simple types)  i: Int, d: double, s: string, v: Vector, t: Transform, u: UnionPar
  74  (Array types)   I: array of Int, D: array of double, V: array of Vector, C: array of char, T: array of Transform, U: array of UnionPars
  75  (Complex types) M: Molecule, G: IntGroup, A: Atom
  76  (Ruby value)    b: Ruby boolean, r: Ruby object, R: an array of Ruby object (a Ruby array)
  77  (Return value)  i, d, s, v, t, r, G + 0x80 */
  78 typedef struct MolActionArg {
  79         Byte type;
  80         union {
  81                 Int ival;
  82                 double dval;
  83                 struct {  /*  Array type: the size of the element is defined by the argument type  */
  84                         Int nitems; /* Number of items */
  85                         void *ptr;
  86                 } arval;
  87                 struct IntGroup *igval; /* The record is retained but not duplicated */
  88                 struct Molecule *mval; /* The record is retained but not duplicated */
  89                 struct Atom *aval; /* The value is duplicated, so any pointer can be passed */
  90                 VALUE vval;        /* The value is retained in sMolActionArgValues  */
  91                 struct {
  92                         void *ptr;    /*  Return value pointer  */
  93                         Int *nptr;    /*  If the return value is an array, then *ptr contains a malloc'ed pointer, and *nptr contains the number of items.  */
  94                 } retval;          /*  Store address for the return value; usually used in performing Ruby script */
  95         } u;
  96 } MolActionArg;
  97
  98 /*  For debug  */
  99 /* static int sMolActionCount = 0; */
 100
 101 /*  Create a new MolAction with the given arguments.
 102  *  Note: If the argument is an array type, it is represented by an integer (number
 103  *  of items) followed by the pointer.  */
 104 MolAction *
 105 MolActionNewArgv(const char *name, va_list ap)
 106 {
 107         MolAction *action;
 108         const char *p;
 109
 110         action = (MolAction *)malloc(sizeof(MolAction));
 111         if (action == NULL)
 112                 goto low_memory;
 113         memset(action, 0, sizeof(MolAction));
 114         action->refCount = 1;
 115         action->frame = -1;
 116         action->name = strdup(name);
 117
 118         /*  For debug  */
 119 /*      fprintf(stderr, "MolAction %p created, count = %d\n", action, ++sMolActionCount); */
 120
 121         /*  Handle arguments  */
 122         p = strchr(name, ':');
 123         if (p == NULL)
 124                 return action;
 125         p++;
 126
 127         while (*p) {
 128                 MolActionArg arg;
 129                 memset(&arg, 0, sizeof(MolActionArg));
 130                 arg.type = *p++;
 131                 switch (arg.type) {
 132                         case 'i':
 133                                 arg.u.ival = va_arg(ap, Int);
 134                                 break;
 135                         case 'd':
 136                                 arg.u.dval = va_arg(ap, double);
 137                                 break;
 138                         case 'I':
 139                         case 'D':
 140                         case 'C':
 141                         case 's':
 142                         case 'V':
 143                         case 'v':
 144                         case 'T':
 145                         case 't':
 146                         case 'U':
 147                         case 'u': {
 148                                 /*  Array type of some kind  */
 149                                 void *ptr;
 150                                 Int nitems, itemsize, allocsize;
 151                                 if (arg.type != 's' && arg.type != 'v' && arg.type != 't' && arg.type != 'u')
 152                                         nitems = va_arg(ap, Int);
 153                                 ptr = va_arg(ap, void *);
 154                                 switch (arg.type) {
 155                                         case 'I': itemsize = sizeof(Int); break;
 156                                         case 'D': itemsize = sizeof(double); break;
 157                                         case 'C':
 158                                         case 's': itemsize = sizeof(char); break;
 159                                         case 'V':
 160                                         case 'v': itemsize = sizeof(Vector); break;
 161                                         case 'T':
 162                                         case 't': itemsize = sizeof(Transform); break;
 163                                         case 'U':
 164                                         case 'u': itemsize = sizeof(UnionPar); break;
 165                                 }
 166                                 if (arg.type == 's')
 167                                         nitems = (ptr == NULL ? 0 : strlen((char *)ptr));
 168                                 else if (arg.type == 'v' || arg.type == 't' || arg.type == 'u')
 169                                         nitems = 1;
 170                                 arg.u.arval.nitems = nitems;
 171                                 allocsize = itemsize * nitems + (arg.type == 's'); /* +1 for string type */
 172                                 if (allocsize == 0)
 173                                         arg.u.arval.ptr == NULL;
 174                                 else {
 175                                         arg.u.arval.ptr = calloc(1, allocsize);
 176                                         if (arg.u.arval.ptr == NULL)
 177                                                 goto low_memory;
 178                                         memmove(arg.u.arval.ptr, ptr, itemsize * nitems);
 179                                 }
 180                                 break;
 181                         }
 182                         case 'G':
 183                                 arg.u.igval = va_arg(ap, IntGroup *);
 184                                 if (arg.u.igval != NULL)
 185                                         IntGroupRetain(arg.u.igval);
 186                                 break;
 187                         case 'M':
 188                                 arg.u.mval = va_arg(ap, Molecule *);
 189                                 MoleculeRetain(arg.u.mval);
 190                                 break;
 191                         case 'A': {
 192                                 Atom *aval = va_arg(ap, Atom *);
 193                                 arg.u.aval = (Atom *)malloc(gSizeOfAtomRecord);
 194                                 if (arg.u.aval == NULL)
 195                                         goto low_memory;
 196                                 AtomDuplicate(arg.u.aval, aval);
 197                                 break;
 198                         }
 199                         case 'r':
 200                         case 'R': {
 201                                 int rtype;
 202                                 arg.u.vval = va_arg(ap, VALUE);
 203                                 rtype = TYPE(arg.u.vval);
 204                                 if (rtype != T_NIL && rtype != T_FIXNUM && rtype != T_SYMBOL && rtype != T_TRUE && rtype != T_FALSE) {
 205                                         if (sMolActionArgValues == Qfalse) {
 206                                                 sMolActionArgValues = rb_ary_new();
 207                                                 rb_global_variable(&sMolActionArgValues);
 208                                         }
 209                                         rb_ary_push(sMolActionArgValues, arg.u.vval);
 210                                 }
 211                                 break;
 212                         }
 213                         case 'b':
 214                                 arg.u.vval = (va_arg(ap, Int) ? Qtrue : Qfalse);
 215                                 break;
 216                         case ';':
 217                                 if (*p == 'i' || *p == 'd' || *p == 's' || *p == 'v' || *p == 't' || *p == 'r' || *p == 'G') {
 218                                         arg.type = (*p | 0x80);
 219                                         p++;
 220                                         arg.u.retval.ptr = va_arg(ap, void *);
 221                                 } else if (*p == 'I' || *p == 'D' || *p == 'V') {
 222                                         arg.type = (*p | 0x80);
 223                                         p++;
 224                                         arg.u.retval.nptr = va_arg(ap, Int *);
 225                                         arg.u.retval.ptr = va_arg(ap, void *);
 226                                 } else {
 227                                         fprintf(stderr, "Internal error: unknown return-type \'%c\' in NewMolAction()", *p);
 228                                 }
 229                                 break;
 230                         default:
 231                                 fprintf(stderr, "Internal error: unknown argument type \'%c\' in NewMolAction()", arg.type);
 232                                 break;
 233                 }
 234                 if (AssignArray(&action->args, &action->nargs, sizeof(arg), action->nargs, &arg) == NULL)
 235                         goto low_memory;
 236         }
 237         return action;
 238
 239   low_memory:
 240         if (action != NULL)
 241                 free(action);
 242         Panic("Low memory during creation of action record");
 243         return NULL;  /* Not reached */
 244 }
 245
 246 MolAction *
 247 MolActionNew(const char *name, ...)
 248 {
 249         va_list ap;
 250         MolAction *retval;
 251         va_start(ap, name);
 252         retval = MolActionNewArgv(name, ap);
 253         va_end(ap);
 254         return retval;
 255 }
 256
 257 MolAction *
 258 MolActionRetain(MolAction *action)
 259 {
 260         if (action != NULL)
 261                 action->refCount++;
 262         return action;
 263 }
 264
 265 void
 266 MolActionRelease(MolAction *action)
 267 {
 268         int i;
 269         if (action == NULL)
 270                 return;
 271         if (--action->refCount > 0)
 272                 return;
 273         if (action->name != NULL)
 274                 free((void *)action->name);
 275         for (i = 0; i < action->nargs; i++) {
 276                 MolActionArg *argp = action->args + i;
 277                 switch (argp->type) {
 278                         case 'i':
 279                         case 'd':
 280                         case 'b':
 281                                 break;
 282                         case 'I':
 283                         case 'D':
 284                         case 'C':
 285                         case 's':
 286                         case 'V':
 287                         case 'v':
 288                         case 'T':
 289                         case 't':
 290                         case 'U':
 291                         case 'u':
 292                                 if (argp->u.arval.ptr != NULL)
 293                                         free(argp->u.arval.ptr);
 294                                 break;
 295                         case 'G':
 296                                 if (argp->u.igval != NULL)
 297                                         IntGroupRelease(argp->u.igval);
 298                                 break;
 299                         case 'M':
 300                                 MoleculeRelease(argp->u.mval);
 301                                 break;
 302                         case 'A':
 303                                 if (argp->u.aval != NULL) {
 304                                         AtomClean(argp->u.aval);
 305                                         free(argp->u.aval);
 306                                 }
 307                                 break;
 308                         case 'r':
 309                         case 'R': {
 310                                 int rtype = TYPE(argp->u.vval);
 311                                 if (rtype != T_NIL && rtype != T_FIXNUM && rtype != T_SYMBOL && rtype != T_TRUE && rtype != T_FALSE) {
 312                                         /*  Look for the same value in sMolActionArgValues and remove it  */
 313                                         if (sMolActionArgValues != Qfalse) {
 314                                                 VALUE *ptr = RARRAY_PTR(sMolActionArgValues);
 315                                                 int n = RARRAY_LEN(sMolActionArgValues);
 316                                                 while (--n >= 0) {
 317                                                         if (ptr[n] == argp->u.vval) {
 318                                                                 rb_ary_delete_at(sMolActionArgValues, n);
 319                                                                 break;
 320                                                         }
 321                                                 }
 322                                         }
 323                                 }
 324                                 break;
 325                         }
 326                         default:
 327                                 break;
 328                 }
 329         }
 330         if (action->args != NULL)
 331                 free(action->args);
 332         free(action);
 333
 334         /*  For debug  */
 335 /*      fprintf(stderr, "MolAction %p disposed, count = %d\n", action, --sMolActionCount);      */
 336 }
 337
 338 void
 339 MolActionSetFrame(MolAction *action, int frame)
 340 {
 341         if (action != NULL)
 342                 action->frame = frame;
 343 }
 344
 345 typedef struct MolRubyActionInfo {
 346         Molecule *mol;
 347         MolAction *action;
 348         VALUE receiver;
 349         ID method_id;
 350         VALUE ary;
 351         char enable_interrupt;
 352         char return_type;
 353         void *return_ptr;
 354         Int *return_nptr;
 355 } MolRubyActionInfo;
 356
 357 static VALUE
 358 s_MolActionToRubyArguments(VALUE vinfo)
 359 {
 360         MolRubyActionInfo *info = (MolRubyActionInfo *)vinfo;
 361         VALUE val;
 362         int i, argc;
 363         char *s, *p;
 364
 365         argc = info->action->nargs - 1;
 366         info->return_type = 0;
 367         info->return_ptr = NULL;
 368         info->ary = rb_ary_new2(argc);
 369         for (i = 1; i <= argc; i++) {
 370                 MolActionArg *argp = &(info->action->args[i]);
 371                 if (argp->type & 0x80) {
 372                         /*  Return value specified  */
 373                         info->return_type = (argp->type & 0x7f);
 374                         info->return_ptr = argp->u.retval.ptr;
 375                         info->return_nptr = argp->u.retval.nptr;  /*  May be unused  */
 376                         break;  /*  This should be the last argument  */
 377                 }
 378                 switch (argp->type) {
 379                         case 'i':
 380                                 val = INT2NUM(argp->u.ival);
 381                                 break;
 382                         case 'd':
 383                                 val = rb_float_new(argp->u.dval);
 384                                 break;
 385                         case 's':
 386                         case 'C':
 387                                 val = rb_str_new((char *)argp->u.arval.ptr, argp->u.arval.nitems);
 388                                 break;
 389                         case 'v':
 390                                 val = ValueFromVector((Vector *)argp->u.arval.ptr);
 391                                 break;
 392                         case 't':
 393                                 val = ValueFromTransform((Transform *)argp->u.arval.ptr);
 394                                 break;
 395                         case 'I':
 396                         case 'D':
 397                         case 'V':
 398                         case 'T': {
 399                                 int j;
 400                                 val = rb_ary_new();
 401                                 for (j = 0; j < argp->u.arval.nitems; j++) {
 402                                         VALUE val1;
 403                                         void *ptr = argp->u.arval.ptr;
 404                                         switch (argp->type) {
 405                                                 case 'I': val1 = INT2NUM(((Int *)ptr)[j]); break;
 406                                                 case 'D': val1 = rb_float_new(((double *)ptr)[j]); break;
 407                                                 case 'V': val1 = ValueFromVector(&((Vector *)ptr)[j]); break;
 408                                                 case 'T': val1 = ValueFromTransform(&((Transform *)ptr)[j]); break;
 409                                         }
 410                                         rb_ary_push(val, val1);
 411                                 }
 412                                 break;
 413                         }
 414                         case 'G':
 415                                 val = ValueFromIntGroup(info->action->args[i].u.igval);
 416                                 break;
 417                         case 'M':
 418                                 val = ValueFromMolecule(info->action->args[i].u.mval);
 419                                 break;
 420                         case 'r':
 421                         case 'b':
 422                                 val = (VALUE)info->action->args[i].u.vval;
 423                                 break;
 424                         case 'R':
 425                                 rb_ary_concat(info->ary, (VALUE)info->action->args[i].u.vval);
 426                                 continue;  /*  Skip rb_ary_push()  */
 427                         default:
 428                                 rb_raise(rb_eArgError, "Unsupported argument type '%c'", info->action->args[i].type);
 429                 }
 430                 rb_ary_push(info->ary, val);
 431         }
 432
 433         if (info->mol != NULL) {
 434                 info->receiver = ValueFromMolecule(info->mol);
 435         } else {
 436                 info->receiver = rb_cMolecule;  /*  Assume class method  */
 437         }
 438
 439         s = info->action->args[0].u.arval.ptr;
 440         for (p = s; *p != 0; p++) {
 441                 if (isalpha(*p) || *p == '_' || (p > s && isdigit(*p)) || (p[1] == 0 && (*p == '?' || *p == '!')))
 442                         continue;
 443                 break;
 444         }
 445
 446         if (*p == 0) {
 447                 /*  Can be a method name  */
 448                 info->method_id = rb_intern(s);
 449         } else {
 450                 /*  Cannot be a method name: try to create a proc object and call it  */
 451                 VALUE procstr = rb_str_new2(s);
 452                 VALUE proc = rb_obj_instance_eval(1, &procstr, info->receiver);
 453                 info->receiver = proc;
 454                 info->method_id = rb_intern("call");
 455         }
 456
 457         return Qnil;
 458 }
 459
 460 static void
 461 s_MolActionStoreReturnValue(MolRubyActionInfo *info, VALUE val)
 462 {
 463         if (info->return_type != 0 && info->return_ptr != NULL) {
 464                 switch (info->return_type) {
 465                         case 'i': *((Int *)(info->return_ptr)) = NUM2INT(rb_Integer(val)); break;
 466                         case 'd': *((Double *)(info->return_ptr)) = NUM2DBL(rb_Float(val)); break;
 467                         case 's': *((char **)(info->return_ptr)) = strdup(StringValuePtr(val)); break;
 468                         case 'v': VectorFromValue(val, (Vector *)(info->return_ptr)); break;
 469                         case 't': TransformFromValue(val, (Transform *)(info->return_ptr)); break;
 470                         case 'r': *((RubyValue *)(info->return_ptr)) = (RubyValue)val; break;
 471                         case 'G': *((IntGroup **)(info->return_ptr)) = (val == Qnil ? NULL : IntGroupFromValue(val)); break;
 472                         case 'I':
 473                         case 'D':
 474                         case 'V': {
 475                                 int i, n, size;
 476                                 void *p;
 477                                 val = rb_ary_to_ary(val);
 478                                 n = RARRAY_LEN(val);
 479                                 *(info->return_nptr) = n;
 480                                 if (n == 0) {
 481                                         *((void **)(info->return_ptr)) = NULL;
 482                                         break;
 483                                 }
 484                                 if (info->return_type == 'I')
 485                                         size = sizeof(Int);
 486                                 else if (info->return_type == 'D')
 487                                         size = sizeof(Double);
 488                                 else if (info->return_type == 'V')
 489                                         size = sizeof(Vector);
 490                                 else break;
 491                                 *((void **)(info->return_ptr)) = p = calloc(size, n);
 492                                 for (i = 0; i < n; i++) {
 493                                         VALUE rval = (RARRAY_PTR(val))[i];
 494                                         if (info->return_type == 'I') {
 495                                                 ((Int *)p)[i] = NUM2INT(rb_Integer(rval));
 496                                         } else if (info->return_type == 'D') {
 497                                                 ((Double *)p)[i] = NUM2DBL(rb_Float(rval));
 498                                         } else if (info->return_type == 'V') {
 499                                                 VectorFromValue(rval, (Vector *)p + i);
 500                                         }
 501                                 }
 502                                 break;
 503                         }
 504                         default: break;
 505                 }
 506         }
 507 }
 508
 509 static VALUE
 510 s_MolActionExecuteRubyScript(VALUE vinfo)
 511 {
 512         VALUE val;
 513         MolRubyActionInfo *info = (MolRubyActionInfo *)vinfo;
 514         val = rb_funcall2(info->receiver, info->method_id, RARRAY_LEN(info->ary), RARRAY_PTR(info->ary));
 515         s_MolActionStoreReturnValue(info, val);
 516         return val;
 517 }
 518
 519 static int
 520 s_MolActionPerformRubyScript(Molecule *mol, MolAction *action)
 521 {
 522         int result;
 523         MolRubyActionInfo info;
 524         memset(&info, 0, sizeof(info));
 525         if (gMolbyIsCheckingInterrupt) {
 526                 MolActionAlertRubyIsRunning();
 527                 return -1;
 528         }
 529         info.mol = mol; /*  May be NULL  */
 530         info.action = action;
 531         gMolbyRunLevel++;
 532         rb_protect(s_MolActionToRubyArguments, (VALUE)&info, &result);
 533         gMolbyRunLevel--;
 534         if (result == 0) {
 535                 VALUE save_interrupt;
 536                 MyAppCallback_beginUndoGrouping();
 537                 save_interrupt = Ruby_SetInterruptFlag(Qtrue);
 538                 gMolbyRunLevel++;
 539                 rb_protect(s_MolActionExecuteRubyScript, (VALUE)&info, &result);
 540                 gMolbyRunLevel--;
 541                 Ruby_SetInterruptFlag(save_interrupt);
 542                 MyAppCallback_endUndoGrouping();
 543         }
 544         if (result != 0) {
 545                 Molby_showError(result);
 546         }
 547         MyAppCallback_hideProgressPanel();  /*  In case when the progress panel is still onscreen */
 548         return (result == 0 ? 0 : -1);
 549 }
 550
 551 static void
 552 s_UpdateSelection(Molecule *mol, const IntGroup *ig, int is_insert)
 553 {
 554         int n, old_natoms;
 555         IntGroup *sel, *orig_atoms, *ig3;
 556         MolAction *act;
 557
 558         if (ig == NULL || (n = IntGroupGetCount(ig)) == 0)
 559                 return;
 560
 561         /*  Insert: orig_atoms = set of original atoms in new indices  */
 562         /*  Delete: orig_atoms = set of remaining atoms in old indices  */
 563         orig_atoms = IntGroupNew();
 564         old_natoms = (is_insert ? mol->natoms - n : mol->natoms + n);
 565         ig3 = IntGroupNewWithPoints(0, (is_insert ? mol->natoms : old_natoms), -1);
 566         IntGroupXor(ig, ig3, orig_atoms);
 567         IntGroupRelease(ig3);
 568
 569         /*  Update selection  */
 570         sel = MoleculeGetSelection(mol);
 571         if (sel != NULL && IntGroupGetCount(sel) > 0) {
 572                 IntGroup *sel2 = IntGroupNew();
 573                 if (is_insert)
 574                         IntGroupConvolute(sel, orig_atoms, sel2);  /*  Renumber  */
 575                 else {
 576                         IntGroup *selx = IntGroupNew();
 577                         IntGroupIntersect(sel, orig_atoms, selx);
 578                         IntGroupDeconvolute(selx, orig_atoms, sel2);
 579                         IntGroupRelease(selx);
 580                 }
 581                 IntGroupRetain(sel);  /*  To avoid being released during MoleculeSetSelection()  */
 582                 MoleculeSetSelection(mol, sel2);
 583                 if (IntGroupGetCount(sel) != IntGroupGetCount(sel2)) {
 584                         /*  Register undo action  */
 585                         act = MolActionNew(gMolActionSetSelection, sel);
 586                         MolActionCallback_registerUndo(mol, act);
 587                         MolActionRelease(act);
 588                 }
 589                 IntGroupRelease(sel2);
 590                 IntGroupRelease(sel);
 591         }
 592
 593         IntGroupRelease(orig_atoms);
 594 }
 595
 596 static int
 597 s_MolActionAddAnAtom(Molecule *mol, MolAction *action, MolAction **actp)
 598 {
 599         Int *ip, n1;
 600         IntGroup *ig;
 601         n1 = MoleculeCreateAnAtom(mol, action->args[0].u.aval, action->args[1].u.ival);
 602         if ((ip = action->args[2].u.retval.ptr) != NULL)
 603                 *ip = n1;
 604         if (n1 < 0)
 605                 return -1;
 606         ig = IntGroupNewWithPoints(n1, 1, -1);
 607         s_UpdateSelection(mol, ig, 1);
 608         IntGroupRelease(ig);
 609         *actp = MolActionNew(gMolActionDeleteAnAtom, n1);
 610         return 0;
 611 }
 612
 613 static int
 614 s_MolActionMergeMolecule(Molecule *mol, MolAction *action, MolAction **actp)
 615 {
 616         int n1, result, minreg, maxreg;
 617         Atom *ap;
 618         IntGroup *ig, *ig2;
 619         Molecule *mol2;
 620         ig = action->args[1].u.igval;
 621         mol2 = action->args[0].u.mval;
 622         if (ig != NULL) {
 623                 ig2 = ig;
 624                 IntGroupRetain(ig2);
 625         } else {
 626                 ig2 = IntGroupNew();
 627                 IntGroupAdd(ig2, mol->natoms, mol2->natoms);
 628         }
 629
 630         /*  Calculate the offset for the residue number  */
 631         maxreg = -1;
 632         for (n1 = 0, ap = mol->atoms; n1 < mol->natoms; n1++, ap = ATOM_NEXT(ap)) {
 633                 if (ap->resSeq > maxreg)
 634                         maxreg = ap->resSeq;
 635         }
 636         minreg = ATOMS_MAX_NUMBER;
 637         for (n1 = 0, ap = mol2->atoms; n1 < mol2->natoms; n1++, ap = ATOM_NEXT(ap)) {
 638                 if (ap->resSeq < minreg)
 639                         minreg = ap->resSeq;
 640         }
 641         if (maxreg < 0)
 642                 maxreg = 0;
 643         if (minreg == ATOMS_MAX_NUMBER)
 644                 minreg = 0;
 645         if (maxreg >= minreg)
 646                 n1 = maxreg - minreg + 1;
 647         else n1 = 0;
 648
 649         if ((result = MoleculeMerge(mol, mol2, ig, n1)) != 0)
 650                 return result;
 651
 652         s_UpdateSelection(mol, ig, 1);
 653
 654         *actp = MolActionNew(gMolActionUnmergeMolecule, ig2);
 655         IntGroupRelease(ig2);
 656         return 0;
 657 }
 658
 659 static int
 660 s_MolActionDeleteAtoms(Molecule *mol, MolAction *action, MolAction **actp)
 661 {
 662         Int n1, result, *ip;
 663         IntGroup *bg, *ig, *pg;
 664         UnionPar *up;
 665         Molecule *mol2;
 666         if (strcmp(action->name, gMolActionDeleteAnAtom) == 0) {
 667                 ig = IntGroupNew();
 668                 if (ig == NULL || IntGroupAdd(ig, action->args[0].u.ival, 1) != 0)
 669                         return -1;
 670         } else {
 671                 ig = action->args[0].u.igval;
 672                 IntGroupRetain(ig);
 673         }
 674         /*  Search bonds crossing the molecule border  */
 675         bg = MoleculeSearchBondsAcrossAtomGroup(mol, ig);
 676         ip = NULL;
 677         if (bg != NULL) {
 678                 n1 = IntGroupGetCount(bg);
 679                 if (n1 > 0) {
 680                         IntGroupIterator iter;
 681                         Int i, idx;
 682                         ip = (Int *)calloc(sizeof(Int), n1 * 2 + 1);
 683                         if (ip == NULL) {
 684                                 IntGroupRelease(bg);
 685                                 IntGroupRelease(ig);
 686                                 return -1;
 687                         }
 688                         IntGroupIteratorInit(bg, &iter);
 689                         idx = 0;
 690                         while ((i = IntGroupIteratorNext(&iter)) >= 0) {
 691                                 /*  The atoms at the border  */
 692                                 ip[idx++] = mol->bonds[i * 2];
 693                                 ip[idx++] = mol->bonds[i * 2 + 1];
 694                         }
 695                         IntGroupIteratorRelease(&iter);
 696                 }
 697                 IntGroupRelease(bg);
 698         }
 699         /*  Search parameters that may disappear after unmerging  */
 700         pg = NULL;
 701         up = NULL;
 702         if (mol->par != NULL) {
 703                 Int i, j, n;
 704                 IntGroup *rg = IntGroupNewWithPoints(0, mol->natoms, -1);
 705                 IntGroupRemoveIntGroup(rg, ig);  /*  Remaining atoms  */
 706                 pg = IntGroupNew();
 707                 for (j = kFirstParType; j <= kLastParType; j++) {
 708                         n = ParameterGetCountForType(mol->par, j);
 709                         for (i = 0; i < n; i++) {
 710                                 UnionPar *up1 = ParameterGetUnionParFromTypeAndIndex(mol->par, j, i);
 711                                 if (!ParameterIsRelevantToAtomGroup(j, up1, mol->atoms, ig) && !ParameterIsRelevantToAtomGroup(j, up1, mol->atoms, rg)) {
 712                                         IntGroupAdd(pg, i + (j - kFirstParType) * kParameterIndexOffset, 1);
 713                                 }
 714                         }
 715                 }
 716                 n = IntGroupGetCount(pg);
 717                 if (n > 0) {
 718                         up = (UnionPar *)calloc(sizeof(UnionPar), n);
 719                         ParameterCopy(mol->par, kFirstParType, up, pg);
 720                 } else {
 721                         IntGroupRelease(pg);
 722                         pg = NULL;
 723                 }
 724         }
 725         /*  Unmerge molecule  */
 726         if ((result = MoleculeUnmerge(mol, &mol2, ig, 0)) != 0) {
 727                 if (ip != NULL)
 728                         free(ip);
 729                 return result;
 730         }
 731
 732         s_UpdateSelection(mol, ig, 0);
 733
 734         if (mol2 == NULL)
 735                 *actp = NULL;
 736         else {
 737                 MolAction *act2;
 738                 /*  If there exist bonds crossing the molecule border, then register
 739                  an undo action to restore them  */
 740                 if (ip != NULL) {
 741                         act2 = MolActionNew(gMolActionAddBonds, n1 * 2, ip);
 742                         MolActionCallback_registerUndo(mol, act2);
 743                         MolActionRelease(act2);
 744                         free(ip);
 745                 }
 746                 /*  Register an undo action to restore the disappearing parameters  */
 747                 if (up != NULL) {
 748                         act2 = MolActionNew(gMolActionAddParameters, kFirstParType, pg, IntGroupGetCount(pg), up);
 749                         MolActionCallback_registerUndo(mol, act2);
 750                         MolActionRelease(act2);
 751                         IntGroupRelease(pg);
 752                         free(up);
 753                 }
 754                 *actp = MolActionNew(gMolActionMergeMolecule, mol2, ig);
 755                 MoleculeRelease(mol2);
 756         }
 757         IntGroupRelease(ig);
 758         return 0;
 759 }
 760
 761 static int
 762 s_MolActionAddStructuralElements(Molecule *mol, MolAction *action, MolAction **actp, int type)
 763 {
 764         Int *ip, n1, result;
 765         IntGroup *ig;
 766
 767         ip = (Int *)action->args[0].u.arval.ptr;
 768
 769         if (type == 0) {  /*  bond  */
 770                 n1 = action->args[0].u.arval.nitems / 2;
 771                 if ((result = MoleculeAddBonds(mol, n1, ip)) <= 0)
 772                         return result;
 773                 ip = (Int *)malloc(sizeof(Int) * 2 * result);
 774                 if (ip == NULL)
 775                         return -4;
 776                 memmove(ip, mol->bonds + (mol->nbonds - result) * 2, sizeof(Int) * 2 * result);
 777                 *actp = MolActionNew(gMolActionDeleteBonds, result * 2, ip);
 778                 free(ip);
 779         } else if (type == 1) {  /*  angle  */
 780                 n1 = action->args[0].u.arval.nitems / 3;
 781                 ig = action->args[1].u.igval;
 782                 if (ig == NULL)
 783                         ig = IntGroupNewWithPoints(mol->nangles, n1, -1);
 784                 else
 785                         IntGroupRetain(ig);
 786                 if ((result = MoleculeAddAngles(mol, ip, ig)) < 0) {
 787                         IntGroupRelease(ig);
 788                         return result;
 789                 }
 790                 *actp = MolActionNew(gMolActionDeleteAngles, ig);
 791                 IntGroupRelease(ig);
 792         } else if (type == 2) {  /*  dihedral  */
 793                 n1 = action->args[0].u.arval.nitems / 4;
 794                 ig = action->args[1].u.igval;
 795                 if (ig == NULL)
 796                         ig = IntGroupNewWithPoints(mol->ndihedrals, n1, -1);
 797                 else
 798                         IntGroupRetain(ig);
 799                 if ((result = MoleculeAddDihedrals(mol, ip, ig)) < 0) {
 800                         IntGroupRelease(ig);
 801                         return result;
 802                 }
 803                 *actp = MolActionNew(gMolActionDeleteDihedrals, ig);
 804                 IntGroupRelease(ig);
 805         } else if (type == 3) {  /*  improper  */
 806                 n1 = action->args[0].u.arval.nitems / 4;
 807                 ig = action->args[1].u.igval;
 808                 if (ig == NULL)
 809                         ig = IntGroupNewWithPoints(mol->nimpropers, n1, -1);
 810                 else
 811                         IntGroupRetain(ig);
 812                 if ((result = MoleculeAddImpropers(mol, ip, ig)) < 0) {
 813                         IntGroupRelease(ig);
 814                         return result;
 815                 }
 816                 *actp = MolActionNew(gMolActionDeleteImpropers, ig);
 817                 IntGroupRelease(ig);
 818         }
 819         return 0;
 820 }
 821
 822 static int
 823 s_MolActionDeleteStructuralElements(Molecule *mol, MolAction *action, MolAction **actp, int type)
 824 {
 825         Int *ip, n1, result;
 826         IntGroup *ig;
 827         if (type == 0) {  /*  bond  */
 828                 ip = (Int *)action->args[0].u.arval.ptr;
 829                 n1 = action->args[0].u.arval.nitems / 2;
 830                 if ((result = MoleculeDeleteBonds(mol, n1, ip)) < 0)
 831                         return result;
 832                 *actp = MolActionNew(gMolActionAddBonds, n1 * 2, ip);
 833         } else if (type == 1) {  /*  angle  */
 834                 ig = action->args[0].u.igval;
 835                 n1 = IntGroupGetCount(ig) * 3;
 836                 ip = (Int *)malloc(sizeof(Int) * n1);
 837                 if ((result = MoleculeDeleteAngles(mol, ip, ig)) < 0) {
 838                         free(ip);
 839                         return result;
 840                 }
 841                 *actp = MolActionNew(gMolActionAddAngles, n1, ip, ig);
 842                 free(ip);
 843         } else if (type == 2) {  /*  dihedral  */
 844                 ig = action->args[0].u.igval;
 845                 n1 = IntGroupGetCount(ig) * 4;
 846                 ip = (Int *)malloc(sizeof(Int) * n1);
 847                 if ((result = MoleculeDeleteDihedrals(mol, ip, ig)) < 0) {
 848                         free(ip);
 849                         return result;
 850                 }
 851                 *actp = MolActionNew(gMolActionAddDihedrals, n1, ip, ig);
 852                 free(ip);
 853         } else if (type == 3) {  /*  improper  */
 854                 ig = action->args[0].u.igval;
 855                 n1 = IntGroupGetCount(ig) * 4;
 856                 ip = (Int *)malloc(sizeof(Int) * n1);
 857                 if ((result = MoleculeDeleteImpropers(mol, ip, ig)) < 0) {
 858                         free(ip);
 859                         return result;
 860                 }
 861                 *actp = MolActionNew(gMolActionAddImpropers, n1, ip, ig);
 862                 free(ip);
 863         }
 864         return 0;
 865 }
 866
 867 static int
 868 s_MolActionTransformAtoms(Molecule *mol, MolAction *action, MolAction **actp, IntGroup **igp, int type)
 869 {
 870         Vector *vp, v;
 871         if (type == 0) {  /*  translate  */
 872                 vp = (Vector *)action->args[0].u.arval.ptr;
 873                 if (vp == NULL)
 874                         return -1;
 875                 *igp = action->args[1].u.igval;
 876                 MoleculeTranslate(mol, vp, *igp);
 877                 VecScale(v, *vp, -1);
 878                 *actp = MolActionNew(gMolActionTranslateAtoms, &v, *igp);
 879                 (*actp)->frame = mol->cframe;
 880         } else if (type == 1) {  /*  rotate  */
 881                 Double ang;
 882                 Vector *vp2;
 883                 vp = (Vector *)action->args[0].u.arval.ptr;
 884                 ang = action->args[1].u.dval;
 885                 vp2 = (Vector *)action->args[2].u.arval.ptr;
 886                 *igp = action->args[3].u.igval;
 887                 MoleculeRotate(mol, vp, ang, vp2, *igp);
 888                 *actp = MolActionNew(gMolActionRotateAtoms, vp, -ang, vp2, *igp);
 889                 (*actp)->frame = mol->cframe;
 890         } else if (type == 2) {  /*  general transform  */
 891                 Transform *trp, tr_inv;
 892                 int atomPositions = 0; /* Save atom positions for undo? */
 893                 trp = (Transform *)action->args[0].u.arval.ptr;
 894                 *igp = action->args[1].u.igval;
 895                 if (TransformInvert(tr_inv, *trp) != 0)
 896                         atomPositions = 1;
 897                 else {
 898                         /*  Check if inverse transform is reliable enough  */
 899                         Transform temp;
 900                         int j;
 901                         TransformMul(temp, *trp, tr_inv);
 902                         temp[0] -= 1.0;
 903                         temp[4] -= 1.0;
 904                         temp[8] -= 1.0;
 905                         for (j = 0; j < 12; j++) {
 906                                 if (temp[j] < -1e-4 || temp[j] > 1e-4)
 907                                         break;
 908                         }
 909                         if (j != 12)
 910                                 atomPositions = 1;
 911                 }
 912                 if (atomPositions) {
 913                         IntGroupIterator iter;
 914                         int j, k, n1;
 915                         n1 = IntGroupGetCount(*igp);
 916                         vp = (Vector *)malloc(sizeof(Vector) * n1);
 917                         if (vp == NULL) {
 918                                 *igp = NULL;
 919                                 return -1;
 920                         }
 921                         IntGroupIteratorInit(*igp, &iter);
 922                         k = 0;
 923                         while ((j = IntGroupIteratorNext(&iter)) >= 0) {
 924                                 vp[k++] = (ATOM_AT_INDEX(mol->atoms, j))->r;
 925                         }
 926                         *actp = MolActionNew(gMolActionSetAtomPositions, *igp, n1, vp);
 927                         free(vp);
 928                         IntGroupIteratorRelease(&iter);
 929                 } else {
 930                         *actp = MolActionNew(gMolActionTransformAtoms, &tr_inv, *igp);
 931                 }
 932                 MoleculeTransform(mol, *trp, *igp);
 933                 (*actp)->frame = mol->cframe;
 934         }
 935         if (*igp != NULL)
 936                 IntGroupRetain(*igp);
 937         return 0;
 938 }
 939
 940 static int
 941 s_MolActionSetAtomGeometry(Molecule *mol, MolAction *action, MolAction **actp, IntGroup **igp, int type)
 942 {
 943         IntGroupIterator iter;
 944         int j, k, n2;
 945         Vector *vp;
 946         Atom *ap;
 947         *igp = action->args[0].u.igval;
 948         n2 = IntGroupGetCount(*igp);
 949         vp = (Vector *)malloc(sizeof(Vector) * n2);
 950         if (vp == NULL) {
 951                 *igp = NULL;
 952                 return -1;
 953         }
 954         IntGroupIteratorInit(*igp, &iter);
 955         k = 0;
 956         while ((j = IntGroupIteratorNext(&iter)) >= 0) {
 957                 ap = ATOM_AT_INDEX(mol->atoms, j);
 958                 vp[k++] = (type == 0 ? ap->r : (type == 1 ? ap->v : ap->f));
 959         }
 960         *actp = MolActionNew(gMolActionSetAtomPositions, *igp, n2, vp);
 961         free(vp);
 962         vp = (Vector *)action->args[1].u.arval.ptr;
 963         IntGroupIteratorReset(&iter);
 964         k = 0;
 965         while ((j = IntGroupIteratorNext(&iter)) >= 0) {
 966                 Vector w = vp[k++];
 967                 ap = ATOM_AT_INDEX(mol->atoms, j);
 968                 if (type == 0)
 969                         ap->r = w;
 970                 else if (type == 1)
 971                         ap->v = w;
 972                 else ap->f = w;
 973         }
 974         IntGroupIteratorRelease(&iter);
 975         (*actp)->frame = mol->cframe;
 976         if (*igp != NULL)
 977                 IntGroupRetain(*igp);
 978         return 0;
 979 }
 980
 981 static int
 982 s_MolActionInsertFrames(Molecule *mol, MolAction *action, MolAction **actp)
 983 {
 984         Int n1, old_nframes, new_nframes, old_cframe;
 985         Vector *vp, *vp2;
 986         IntGroup *ig;
 987         MolAction *act2;
 988
 989         ig = action->args[0].u.igval;
 990         vp = (Vector *)action->args[1].u.arval.ptr;
 991         vp2 = (Vector *)action->args[2].u.arval.ptr;
 992         old_cframe = mol->cframe;
 993         n1 = IntGroupGetCount(ig);
 994         if (n1 == 0)
 995                 return 0;  /*  Do nothing  */
 996         if (vp != NULL && action->args[1].u.arval.nitems != n1 * mol->natoms)
 997                 return -1;  /*  Internal inconsistency  */
 998         if (vp2 != NULL && action->args[2].u.arval.nitems != n1 * 4)
 999                 return -1;  /*  Internal inconsistency  */
1000         old_nframes = MoleculeGetNumberOfFrames(mol);
1001         if (MoleculeInsertFrames(mol, ig, vp, vp2) < 0)
1002                 return -1;  /*  Error  */
1003
1004         /*  Undo action for restoring old cframe  */
1005         act2 = MolActionNew(gMolActionNone);
1006         act2->frame = old_cframe;
1007         MolActionCallback_registerUndo(mol, act2);
1008         MolActionRelease(act2);
1009
1010         new_nframes = MoleculeGetNumberOfFrames(mol);
1011         if (old_nframes + n1 < new_nframes) {
1012                 /*  "Extra" frames were automatically inserted because large frame indices were specified  */
1013                 /*  Register undo operation to remove these extra frames  */
1014                 IntGroup *ig2 = IntGroupNewWithPoints(old_nframes, new_nframes - (old_nframes + n1), -1);
1015                 act2 = MolActionNew(gMolActionRemoveFrames, ig2);
1016                 IntGroupRelease(ig2);
1017                 MolActionCallback_registerUndo(mol, act2);
1018                 MolActionRelease(act2);
1019         }
1020         *actp = MolActionNew(gMolActionRemoveFrames, ig);
1021         (*actp)->frame = mol->cframe;
1022         return 0;
1023 }
1024
1025 static int
1026 s_MolActionRemoveFrames(Molecule *mol, MolAction *action, MolAction **actp)
1027 {
1028         Vector *vp, *vp2;
1029         IntGroup *ig, *ig2;
1030         Int n1, n2, nframes, old_cframe;
1031         MolAction *act2;
1032
1033         ig = ig2 = action->args[0].u.igval;
1034         old_cframe = mol->cframe;
1035         n1 = IntGroupGetCount(ig);
1036         if (n1 == 0)
1037                 return 0;  /*  Do nothing  */
1038         nframes = MoleculeGetNumberOfFrames(mol);
1039         n2 = IntGroupGetEndPoint(ig, IntGroupGetIntervalCount(ig) - 1);  /*  Max point + 1  */
1040         if (n2 > nframes) {
1041                 /*  Remove extra points  */
1042                 ig2 = IntGroupNewFromIntGroup(ig);
1043                 IntGroupRemove(ig2, nframes, n2 - nframes);
1044                 n1 = IntGroupGetCount(ig2);
1045         }
1046         if (nframes == n1 && nframes >= 2) {
1047                 /*  Remove all frames: keep the current frame  */
1048                 if (ig2 == ig)
1049                         ig2 = IntGroupNewFromIntGroup(ig);
1050                 IntGroupRemove(ig2, mol->cframe, 1);
1051                 n1--;
1052         }
1053         if (n1 == 0) {
1054                 if (ig2 != ig)
1055                         IntGroupRelease(ig2);
1056                 return 0;  /*  Do nothing  */
1057         }
1058         vp = (Vector *)calloc(sizeof(Vector), n1 * mol->natoms);
1059         if (mol->cell != NULL && mol->frame_cells != NULL)
1060                 vp2 = (Vector *)calloc(sizeof(Vector) * 4, n1);
1061         else vp2 = NULL;
1062         if (MoleculeRemoveFrames(mol, ig2, vp, vp2) < 0) {
1063                 if (ig2 != ig)
1064                         IntGroupRelease(ig2);
1065                 return -1;  /*  Error  */
1066         }
1067         /*  Undo action for restoring old cframe  */
1068         act2 = MolActionNew(gMolActionNone);
1069         act2->frame = old_cframe;
1070         MolActionCallback_registerUndo(mol, act2);
1071         MolActionRelease(act2);
1072
1073         *actp = MolActionNew(gMolActionInsertFrames, ig2, n1 * mol->natoms, vp, (vp2 != NULL ? n1 * 4 : 0), vp2);
1074         (*actp)->frame = mol->cframe;
1075
1076         free(vp);
1077         if (vp2 != NULL)
1078                 free(vp2);
1079         if (ig2 != ig)
1080                 IntGroupRelease(ig2);
1081         return 0;
1082 }
1083
1084 static int
1085 s_MolActionSetSelection(Molecule *mol, MolAction *action, MolAction **actp)
1086 {
1087         IntGroup *ig, *ig2;
1088         ig2 = MoleculeGetSelection(mol);
1089         if (ig2 != NULL)
1090                 IntGroupRetain(ig2);  /*  To avoid releasing during MoleculeSetSelection() */
1091         ig = action->args[0].u.igval;
1092         MoleculeSetSelection(mol, ig);
1093         *actp = MolActionNew(gMolActionSetSelection, ig2);
1094         if (ig2 != NULL)
1095                 IntGroupRelease(ig2);
1096         return 0;
1097 }
1098
1099 static int
1100 s_MolActionRenumberAtoms(Molecule *mol, MolAction *action, MolAction **actp)
1101 {
1102         Int *ip, n1, result;
1103         Int *ip2 = (Int *)malloc(sizeof(Int) * (mol->natoms + 1));
1104         if (ip2 == NULL)
1105                 return -1;
1106         ip = (Int *)action->args[0].u.arval.ptr;
1107         n1 = action->args[0].u.arval.nitems;
1108         result = MoleculeRenumberAtoms(mol, ip, ip2, n1);
1109         if (result != 0) {
1110                 free(ip2);
1111                 return result;
1112         }
1113         *actp = MolActionNew(gMolActionRenumberAtoms, mol->natoms, ip2);
1114         return 0;
1115 }
1116
1117 static int
1118 s_MolActionChangeResidueNumber(Molecule *mol, MolAction *action, MolAction **actp, int forUndo)
1119 {
1120         IntGroup *ig;
1121         IntGroupIterator iter;
1122         Int i, n1, *ip, nresidues;
1123
1124         ig = action->args[0].u.igval;
1125         n1 = IntGroupGetCount(ig);
1126         ip = (Int *)calloc(sizeof(Int), n1 + 1);
1127         IntGroupIteratorInit(ig, &iter);
1128         i = 0;
1129         while ((n1 = IntGroupIteratorNext(&iter)) >= 0) {
1130                 ip[i++] = ATOM_AT_INDEX(mol->atoms, n1)->resSeq;
1131         }
1132         n1 = i;
1133         ip[i] = kInvalidIndex;
1134         nresidues = mol->nresidues;
1135         if (forUndo) {
1136                 MoleculeChangeResidueNumberWithArray(mol, ig, (Int *)action->args[1].u.arval.ptr);
1137                 MoleculeChangeNumberOfResidues(mol, action->args[2].u.ival);
1138         } else {
1139                 MoleculeChangeResidueNumber(mol, ig, action->args[1].u.ival);
1140         }
1141         *actp = MolActionNew(gMolActionChangeResidueNumberForUndo, ig, n1, ip, nresidues);
1142         return 0;
1143 }
1144
1145 static int
1146 s_MolActionOffsetResidueNumbers(Molecule *mol, MolAction *action, MolAction **actp)
1147 {
1148         IntGroup *ig;
1149         Int n1, result;
1150         ig = action->args[0].u.igval;
1151         n1 = mol->nresidues;
1152         result = MoleculeOffsetResidueNumbers(mol, ig, action->args[1].u.ival, action->args[2].u.ival);
1153         if (result != 0)
1154                 return result;  /*  The molecule is not modified  */
1155         *actp = MolActionNew(gMolActionOffsetResidueNumbers, ig, -(action->args[1].u.ival), n1);
1156         return 0;
1157 }
1158
1159 static int
1160 s_MolActionChangeResidueNames(Molecule *mol, MolAction *action, MolAction **actp)
1161 {
1162         char *new_names, *old_names;
1163         Int *ip, i, argc;
1164         ip = (Int *)action->args[0].u.arval.ptr;
1165         argc = action->args[0].u.arval.nitems;
1166         new_names = action->args[1].u.arval.ptr;
1167         old_names = (char *)malloc(argc * 4 + 1);
1168         for (i = 0; i < argc; i++) {
1169                 if (ip[i] >= mol->nresidues)
1170                         old_names[i * 4] = 0;
1171                 else
1172                         strncpy(old_names + i * 4, mol->residues[ip[i]], 4);
1173         }
1174         MoleculeChangeResidueNames(mol, argc, ip, new_names);
1175         *actp = MolActionNew(gMolActionChangeResidueNames, argc, ip, argc * 4, old_names);
1176         free(old_names);
1177         return 0;
1178 }
1179
1180 static int
1181 s_MolActionChangeNumberOfResidues(Molecule *mol, MolAction *action, MolAction **actp)
1182 {
1183         Int n1, nresidues = mol->nresidues;
1184         n1 = action->args[0].u.ival;
1185         if (n1 < nresidues) {
1186                 /*  The residue names will be lost, so undo must be registered to restore the names  */
1187                 int argc = nresidues - n1;
1188                 char *names = (char *)malloc(4 * argc + 1);
1189                 Int *ip = (Int *)malloc(sizeof(Int) * (argc + 1));
1190                 int i;
1191                 MolAction *act2;
1192                 for (i = 0; i < argc; i++) {
1193                         strncpy(names + i * 4, mol->residues[i + n1], 4);
1194                         ip[i] = i + n1;
1195                 }
1196                 ip[i] = kInvalidIndex;
1197                 act2 = MolActionNew(gMolActionChangeResidueNames, argc, ip, argc * 4, names);
1198                 MolActionCallback_registerUndo(mol, act2);
1199                 free(ip);
1200                 free(names);
1201         }
1202         MoleculeChangeNumberOfResidues(mol, n1);
1203         *actp = MolActionNew(gMolActionChangeNumberOfResidues, nresidues);
1204         return 0;
1205 }
1206
1207 static int
1208 s_MolActionExpandBySymmetry(Molecule *mol, MolAction *action, MolAction **actp)
1209 {
1210         Int n1, *ip, count;
1211         Symop symop;
1212         IntGroup *ig = action->args[0].u.igval;
1213         count = IntGroupGetCount(ig);
1214         if (count == 0)
1215                 return 0;  /*  No operation  */
1216         symop.dx = action->args[1].u.ival;
1217         symop.dy = action->args[2].u.ival;
1218         symop.dz = action->args[3].u.ival;
1219         symop.sym = action->args[4].u.ival;
1220         symop.alive = (symop.dx != 0 || symop.dy != 0 || symop.dz != 0 || symop.sym != 0);
1221         if (action->args[5].u.retval.ptr != NULL) {
1222                 /*  Request the indices of the atoms  */
1223                 ip = (Int *)calloc(sizeof(Int), count);
1224         } else ip = NULL;
1225         n1 = MoleculeAddExpandedAtoms(mol, symop, ig, ip);
1226         if (n1 > 0) {
1227                 ig = IntGroupNew();
1228                 IntGroupAdd(ig, mol->natoms - n1, n1);
1229                 (*actp) = MolActionNew(gMolActionUnmergeMolecule, ig);
1230                 IntGroupRelease(ig);
1231         }
1232         if (ip != NULL) {
1233                 if (n1 < 0) {
1234                         free(ip);
1235                         ip = NULL;
1236                         count = 0;
1237                 }
1238                 *((Int **)(action->args[5].u.retval.ptr)) = ip;
1239                 *(action->args[5].u.retval.nptr) = count;
1240         }
1241         return n1;
1242 }
1243
1244 static int
1245 s_MolActionAmendBySymmetry(Molecule *mol, MolAction *action, MolAction **actp)
1246 {
1247         IntGroup *ig1, *ig2;
1248         Vector *vp2;
1249         int n1;
1250         ig1 = action->args[0].u.igval;
1251         n1 = MoleculeAmendBySymmetry(mol, ig1, &ig2, &vp2);
1252         if (action->args[1].u.retval.ptr != NULL) {
1253                 *((IntGroup **)(action->args[1].u.retval.ptr)) = ig2;
1254                 IntGroupRetain(ig2);
1255         }
1256         if (n1 > 0) {
1257                 *actp = MolActionNew(gMolActionSetAtomPositions, ig2, n1, vp2);
1258                 free(vp2);
1259                 IntGroupRelease(ig2);
1260         }
1261         mol->needsMDCopyCoordinates = 1;
1262         return 0;
1263 }
1264
1265 static int
1266 s_MolActionAddSymmetryOperation(Molecule *mol, MolAction *action, MolAction **actp)
1267 {
1268         Int n1;
1269         Transform *trp;
1270         Transform itr = {1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0};
1271         trp = (Transform *)action->args[0].u.arval.ptr;
1272         if (mol->nsyms == 0) {
1273                 for (n1 = 0; n1 < 12; n1++) {
1274                         if (fabs((*trp)[n1] - itr[n1]) > 1e-8)
1275                                 break;
1276                 }
1277                 if (n1 < 12) {
1278                         MolAction *act2;
1279                         if (AssignArray(&mol->syms, &mol->nsyms, sizeof(Transform), mol->nsyms, &itr) == 0)
1280                                 return -1;
1281                         act2 = MolActionNew(gMolActionDeleteSymmetryOperation);
1282                         MolActionCallback_registerUndo(mol, act2);
1283                         MolActionRelease(act2);
1284                 }
1285         }
1286         if (AssignArray(&mol->syms, &mol->nsyms, sizeof(Transform), mol->nsyms, trp) == 0)
1287                 return -1;
1288         *actp = MolActionNew(gMolActionDeleteSymmetryOperation);
1289         return 0;
1290 }
1291
1292 static int
1293 s_MolActionDeleteSymmetryOperation(Molecule *mol, MolAction *action, MolAction **actp)
1294 {
1295         if (mol->nsyms == 0)
1296                 return -1;
1297         *actp = MolActionNew(gMolActionAddSymmetryOperation, &(mol->syms[mol->nsyms - 1]));
1298         mol->nsyms--;
1299         if (mol->nsyms == 0) {
1300                 free(mol->syms);
1301                 mol->syms = NULL;
1302         }
1303         return 0;
1304 }
1305
1306 static int
1307 s_MolActionSetCell(Molecule *mol, MolAction *action, MolAction **actp)
1308 {
1309         double *dp, d[12];
1310         Int convertCoord, n1, n2;
1311         if (mol->cell == NULL) {
1312                 d[0] = 0.0;
1313                 n1 = 0;
1314         } else {
1315                 for (n1 = 0; n1 < 6; n1++)
1316                         d[n1] = mol->cell->cell[n1];
1317                 if (mol->cell->has_sigma) {
1318                         for (n1 = 6; n1 < 12; n1++)
1319                                 d[n1] = mol->cell->cellsigma[n1 - 6];
1320                 }
1321         }
1322         convertCoord = action->args[1].u.ival;
1323         dp = action->args[0].u.arval.ptr;
1324         n2 = action->args[0].u.arval.nitems;
1325         if (n2 == 0)
1326                 MoleculeSetCell(mol, 0, 0, 0, 0, 0, 0, convertCoord);
1327         else {
1328                 MoleculeSetCell(mol, dp[0], dp[1], dp[2], dp[3], dp[4], dp[5], convertCoord);
1329                 if (n2 == 12) {
1330                         mol->cell->has_sigma = 1;
1331                         for (n2 = 6; n2 < 12; n2++)
1332                                 mol->cell->cellsigma[n2 - 6] = dp[n2];
1333                 } else mol->cell->has_sigma = 0;
1334         }
1335         *actp = MolActionNew(gMolActionSetCell, n1, (n1 == 0 ? NULL : d), convertCoord);
1336         return 0;
1337 }
1338
1339 static int
1340 s_MolActionSetBox(Molecule *mol, MolAction *action, MolAction **actp)
1341 {
1342         Int n1, n2;
1343         if (action == NULL) {
1344                 /*  Clear box  */
1345                 if (mol->cell == NULL)
1346                         return 0;  /*  Do nothing  */
1347                 n1 = ((mol->cell->flags[0] != 0) * 4 + (mol->cell->flags[1] != 0) * 2 + (mol->cell->flags[2] != 0));
1348                 *actp = MolActionNew(gMolActionSetBox, &(mol->cell->axes[0]), &(mol->cell->axes[1]), &(mol->cell->axes[2]), &(mol->cell->origin), n1);
1349                 MoleculeSetPeriodicBox(mol, NULL, NULL, NULL, NULL, NULL);
1350         } else {
1351                 /*  Set box  */
1352                 Vector v[4];
1353                 char flags[3];
1354                 if (mol->cell == NULL)
1355                         *actp = MolActionNew(gMolActionClearBox);
1356                 else {
1357                         n1 = ((mol->cell->flags[0] != 0) * 4 + (mol->cell->flags[1] != 0) * 2 + (mol->cell->flags[2] != 0));
1358                         *actp = MolActionNew(gMolActionSetBox, &(mol->cell->axes[0]), &(mol->cell->axes[1]), &(mol->cell->axes[2]), &(mol->cell->origin), n1);
1359                 }
1360                 for (n1 = 0; n1 < 4; n1++)
1361                         v[n1] = *((Vector *)(action->args[n1].u.arval.ptr));
1362                 n2 = action->args[4].u.ival;
1363                 if (n2 < 0) {
1364                         /*  Keep existing flags; if not present, set all flags to 1.  */
1365                         if (mol->cell == NULL)
1366                                 flags[0] = flags[1] = flags[2] = 1;
1367                         else {
1368                                 flags[0] = mol->cell->flags[0];
1369                                 flags[1] = mol->cell->flags[1];
1370                                 flags[2] = mol->cell->flags[2];
1371                         }
1372                 } else {
1373                         for (n1 = 0; n1 < 3; n1++)
1374                                 flags[n1] = ((n2 >> (2 - n1)) & 1);
1375                 }
1376                 MoleculeSetPeriodicBox(mol, &v[0], &v[1], &v[2], &v[3], flags);
1377         }
1378         return 0;
1379 }
1380
1381 /*  This action is used for undoing "cell_flexibility = false"  */
1382 static int
1383 s_MolActionSetBoxForFrames(Molecule *mol, MolAction *action, MolAction **actp)
1384 {
1385         Int i, n1, n2;
1386         Vector *vp1, *vp2;
1387         n2 = MoleculeGetNumberOfFrames(mol);
1388         if (n2 == 0 || mol->cell == NULL)
1389                 return 0;  /*  Do nothing  */
1390         n1 = action->args[0].u.arval.nitems / 4;
1391         vp1 = (Vector *)(action->args[0].u.arval.ptr);
1392         if (mol->nframe_cells < n2) {
1393                 /*  Expand the array before processing  */
1394                 i = mol->nframe_cells * 4;
1395                 AssignArray(&(mol->frame_cells), &(mol->nframe_cells), sizeof(Vector) * 4, n2 - 1, NULL);
1396                 while (i < n2 * 4) {
1397                         /*  Copy the current cell  */
1398                         mol->frame_cells[i++] = mol->cell->axes[0];
1399                         mol->frame_cells[i++] = mol->cell->axes[1];
1400                         mol->frame_cells[i++] = mol->cell->axes[2];
1401                         mol->frame_cells[i++] = mol->cell->origin;
1402                 }
1403         }
1404
1405         vp2 = (Vector *)malloc(sizeof(Vector) * n2 * 4);
1406         memmove(vp2, mol->frame_cells, sizeof(Vector) * n2 * 4);
1407         memmove(mol->frame_cells, vp1, sizeof(Vector) * 4 * (n1 < n2 ? n1 : n2));
1408         *actp = MolActionNew(gMolActionSetBoxForFrames, n2 * 4, vp2);
1409         free(vp2);
1410
1411         /*  Set the current cell (no change on the periodic flags)  */
1412         vp2 = mol->frame_cells + mol->cframe * 4;
1413         MoleculeSetPeriodicBox(mol, vp2, vp2 + 1, vp2 + 2, vp2 + 3, mol->cell->flags);
1414
1415         return 0;
1416 }
1417
1418 static int
1419 s_MolActionSetCellFlexibility(Molecule *mol, MolAction *action, MolAction **actp)
1420 {
1421         Int n1;
1422         n1 = action->args[0].u.ival;
1423         if ((n1 != 0) == (mol->useFlexibleCell != 0))
1424                 return 0;  /*  Do nothing  */
1425         mol->useFlexibleCell = (n1 != 0);
1426         if (n1 == 0) {
1427                 /*  Clear the existing cells, and register undo  */
1428                 if (mol->nframe_cells > 0) {
1429                         MolAction *act2 = MolActionNew(gMolActionSetBoxForFrames, mol->nframe_cells * 4, mol->frame_cells);
1430                         MolActionSetFrame(act2, mol->cframe);
1431                         MolActionCallback_registerUndo(mol, act2);
1432                         MolActionRelease(act2);
1433                 }
1434                 free(mol->frame_cells);
1435                 mol->frame_cells = NULL;
1436                 mol->nframe_cells = 0;
1437         } else {
1438                 /*  Allocate cells for all frames and copy the current cell  */
1439                 Int i, nframes = MoleculeGetNumberOfFrames(mol);
1440                 if (nframes != 0 && mol->cell != NULL) {
1441                         if (mol->nframe_cells < nframes) {
1442                                 /*  Expand the array  */
1443                                 AssignArray(&(mol->frame_cells), &(mol->nframe_cells), sizeof(Vector) * 4, nframes - 1, NULL);
1444                         }
1445                         /*  Copy the current cell  */
1446                         /*  (No undo action is registered; actually, the frame_cells array should be empty)  */
1447                         for (i = 0; i < nframes; i++) {
1448                                 mol->frame_cells[i * 4] = mol->cell->axes[0];
1449                                 mol->frame_cells[i * 4 + 1] = mol->cell->axes[1];
1450                                 mol->frame_cells[i * 4 + 2] = mol->cell->axes[2];
1451                                 mol->frame_cells[i * 4 + 3] = mol->cell->origin;
1452                         }
1453                 }
1454         }
1455         *actp = MolActionNew(gMolActionSetCellFlexibility, (n1 == 0));
1456         return 0;
1457 }
1458
1459 static int
1460 s_MolActionAddParameters(Molecule *mol, MolAction *action, MolAction **actp)
1461 {
1462         UnionPar *up;
1463         IntGroup *ig;
1464         Int parType;
1465         if (mol->par == NULL)
1466                 return -1;
1467         parType = action->args[0].u.ival;
1468         ig = action->args[1].u.igval;
1469         up = action->args[2].u.arval.ptr;
1470         ParameterInsert(mol->par, parType, up, ig);
1471         *actp = MolActionNew(gMolActionDeleteParameters, parType, ig);
1472         mol->parameterTableSelectionNeedsClear = 1;
1473         return 0;
1474 }
1475
1476 static int
1477 s_MolActionDeleteParameters(Molecule *mol, MolAction *action, MolAction **actp)
1478 {
1479         UnionPar *up;
1480         IntGroup *ig;
1481         Int parType, n1;
1482         if (mol->par == NULL)
1483                 return -1;
1484         parType = action->args[0].u.ival;
1485         ig = action->args[1].u.igval;
1486         n1 = IntGroupGetCount(ig);
1487         up = (UnionPar *)calloc(sizeof(UnionPar), n1);
1488         ParameterDelete(mol->par, parType, up, ig);
1489         *actp = MolActionNew(gMolActionAddParameters, parType, ig, n1, up);
1490         free(up);
1491         mol->parameterTableSelectionNeedsClear = 1;
1492         return 0;
1493 }
1494
1495 int
1496 MolActionPerform(Molecule *mol, MolAction *action)
1497 {
1498         int result;
1499         IntGroup *ig = NULL;
1500         MolAction *act2 = NULL;
1501         int needsSymmetryAmendment = 0;
1502         int needsRebuildMDArena = 0;
1503
1504         if (action == NULL)
1505                 return -1;
1506
1507         if (action->frame >= 0 && mol->cframe != action->frame)
1508                 MoleculeSelectFrame(mol, action->frame, 1);
1509
1510         /*  Ruby script execution  */
1511         if (strncmp(action->name, kMolActionPerformScript, strlen(kMolActionPerformScript)) == 0) {
1512                 return s_MolActionPerformRubyScript(mol, action);
1513         }
1514
1515         if (mol == NULL)
1516                 return -1;
1517
1518         if (strcmp(action->name, gMolActionAddAnAtom) == 0) {
1519                 if ((result = s_MolActionAddAnAtom(mol, action, &act2)) != 0)
1520                         return result;
1521                 needsRebuildMDArena = 1;
1522         } else if (strcmp(action->name, gMolActionMergeMolecule) == 0) {
1523                 if ((result = s_MolActionMergeMolecule(mol, action, &act2)) != 0)
1524                         return result;
1525                 needsRebuildMDArena = 1;
1526         } else if (strcmp(action->name, gMolActionDeleteAnAtom) == 0 || strcmp(action->name, gMolActionUnmergeMolecule) == 0) {
1527                 if ((result = s_MolActionDeleteAtoms(mol, action, &act2)) != 0)
1528                         return result;
1529                 needsRebuildMDArena = 1;
1530         } else if (strcmp(action->name, gMolActionAddBonds) == 0) {
1531                 if ((result = s_MolActionAddStructuralElements(mol, action, &act2, 0)) != 0)
1532                         return result;
1533                 needsRebuildMDArena = 1;
1534         } else if (strcmp(action->name, gMolActionDeleteBonds) == 0) {
1535                 if ((result = s_MolActionDeleteStructuralElements(mol, action, &act2, 0)) != 0)
1536                         return result;
1537                 needsRebuildMDArena = 1;
1538         } else if (strcmp(action->name, gMolActionAddAngles) == 0) {
1539                 if ((result = s_MolActionAddStructuralElements(mol, action, &act2, 1)) != 0)
1540                         return result;
1541                 needsRebuildMDArena = 1;
1542         } else if (strcmp(action->name, gMolActionDeleteAngles) == 0) {
1543                 if ((result = s_MolActionDeleteStructuralElements(mol, action, &act2, 1)) != 0)
1544                         return result;
1545                 needsRebuildMDArena = 1;
1546         } else if (strcmp(action->name, gMolActionAddDihedrals) == 0) {
1547                 if ((result = s_MolActionAddStructuralElements(mol, action, &act2, 2)) != 0)
1548                         return result;
1549                 needsRebuildMDArena = 1;
1550         } else if (strcmp(action->name, gMolActionDeleteDihedrals) == 0) {
1551                 if ((result = s_MolActionDeleteStructuralElements(mol, action, &act2, 2)) != 0)
1552                         return result;
1553                 needsRebuildMDArena = 1;
1554         } else if (strcmp(action->name, gMolActionAddImpropers) == 0) {
1555                 if ((result = s_MolActionAddStructuralElements(mol, action, &act2, 3)) != 0)
1556                         return result;
1557                 needsRebuildMDArena = 1;
1558         } else if (strcmp(action->name, gMolActionDeleteImpropers) == 0) {
1559                 if ((result = s_MolActionDeleteStructuralElements(mol, action, &act2, 3)) != 0)
1560                         return result;
1561                 needsRebuildMDArena = 1;
1562         } else if (strcmp(action->name, gMolActionTranslateAtoms) == 0) {
1563                 if ((result = s_MolActionTransformAtoms(mol, action, &act2, &ig, 0)) != 0)
1564                         return result;
1565                 needsSymmetryAmendment = 1;
1566         } else if (strcmp(action->name, gMolActionRotateAtoms) == 0) {
1567                 if ((result = s_MolActionTransformAtoms(mol, action, &act2, &ig, 1)) != 0)
1568                         return result;
1569                 needsSymmetryAmendment = 1;
1570         } else if (strcmp(action->name, gMolActionTransformAtoms) == 0) {
1571                 if ((result = s_MolActionTransformAtoms(mol, action, &act2, &ig, 2)) != 0)
1572                         return result;
1573                 needsSymmetryAmendment = 1;
1574         } else if (strcmp(action->name, gMolActionSetAtomPositions) == 0) {
1575                 if ((result = s_MolActionSetAtomGeometry(mol, action, &act2, &ig, 0)) != 0)
1576                         return result;
1577                 needsSymmetryAmendment = 1;
1578         } else if (strcmp(action->name, gMolActionSetAtomVelocities) == 0) {
1579                 if ((result = s_MolActionSetAtomGeometry(mol, action, &act2, &ig, 1)) != 0)
1580                         return result;
1581                 needsSymmetryAmendment = 1;
1582         } else if (strcmp(action->name, gMolActionSetAtomForces) == 0) {
1583                 if ((result = s_MolActionSetAtomGeometry(mol, action, &act2, &ig, 2)) != 0)
1584                         return result;
1585                 needsSymmetryAmendment = 1;
1586         } else if (strcmp(action->name, gMolActionInsertFrames) == 0) {
1587                 if ((result = s_MolActionInsertFrames(mol, action, &act2)) != 0)
1588                         return result;
1589                 needsSymmetryAmendment = 1;
1590         } else if (strcmp(action->name, gMolActionRemoveFrames) == 0) {
1591                 if ((result = s_MolActionRemoveFrames(mol, action, &act2)) != 0)
1592                         return result;
1593                 needsSymmetryAmendment = 1;
1594         } else if (strcmp(action->name, gMolActionSetSelection) == 0) {
1595                 if ((result = s_MolActionSetSelection(mol, action, &act2)) != 0)
1596                         return result;
1597         } else if (strcmp(action->name, gMolActionRenumberAtoms) == 0) {
1598                 if ((result = s_MolActionRenumberAtoms(mol, action, &act2)) != 0)
1599                         return result;
1600                 needsRebuildMDArena = 1;
1601         } else if (strcmp(action->name, gMolActionChangeResidueNumber) == 0) {
1602                 if ((result = s_MolActionChangeResidueNumber(mol, action, &act2, 0)) != 0)
1603                         return result;
1604         } else if (strcmp(action->name, gMolActionChangeResidueNumberForUndo) == 0) {
1605                 if ((result = s_MolActionChangeResidueNumber(mol, action, &act2, 1)) != 0)
1606                         return result;
1607         } else if (strcmp(action->name, gMolActionOffsetResidueNumbers) == 0) {
1608                 if ((result = s_MolActionOffsetResidueNumbers(mol, action, &act2)) != 0)
1609                         return result;
1610         } else if (strcmp(action->name, gMolActionChangeResidueNames) == 0) {
1611                 if ((result = s_MolActionChangeResidueNames(mol, action, &act2)) != 0)
1612                         return result;
1613         } else if (strcmp(action->name, gMolActionChangeNumberOfResidues) == 0) {
1614                 if ((result = s_MolActionChangeNumberOfResidues(mol, action, &act2)) != 0)
1615                         return result;
1616         } else if (strcmp(action->name, gMolActionAmendBySymmetry) == 0) {
1617                 if ((result = s_MolActionAmendBySymmetry(mol, action, &act2)) != 0)
1618                         return result;
1619                 needsRebuildMDArena = 1;
1620         } else if (strcmp(action->name, gMolActionExpandBySymmetry) == 0) {
1621                 if ((result = s_MolActionExpandBySymmetry(mol, action, &act2)) != 0)
1622                         return result;
1623         } else if (strcmp(action->name, gMolActionAddSymmetryOperation) == 0) {
1624                 if ((result = s_MolActionAddSymmetryOperation(mol, action, &act2)) != 0)
1625                         return result;
1626                 needsRebuildMDArena = 1;
1627         } else if (strcmp(action->name, gMolActionDeleteSymmetryOperation) == 0) {
1628                 if ((result = s_MolActionDeleteSymmetryOperation(mol, action, &act2)) != 0)
1629                         return result;
1630                 needsRebuildMDArena = 1;
1631         } else if (strcmp(action->name, gMolActionSetCell) == 0) {
1632                 if ((result = s_MolActionSetCell(mol, action, &act2)) != 0)
1633                         return result;
1634                 if (mol->arena != NULL)
1635                         md_set_cell(mol->arena);
1636                 needsRebuildMDArena = 1;
1637         } else if (strcmp(action->name, gMolActionSetBox) == 0) {
1638                 if ((result = s_MolActionSetBox(mol, action, &act2)) != 0)
1639                         return result;
1640                 if (mol->arena != NULL)
1641                         md_set_cell(mol->arena);
1642                 needsSymmetryAmendment = 1;
1643         } else if (strcmp(action->name, gMolActionClearBox) == 0) {
1644                 if ((result = s_MolActionSetBox(mol, NULL, &act2)) != 0)
1645                         return result;
1646                 if (mol->arena != NULL)
1647                         md_set_cell(mol->arena);
1648                 needsSymmetryAmendment = 1;
1649         } else if (strcmp(action->name, gMolActionSetBoxForFrames) == 0) {
1650                 if ((result = s_MolActionSetBoxForFrames(mol, action, &act2)) != 0)
1651                         return result;
1652         } else if (strcmp(action->name, gMolActionSetCellFlexibility) == 0) {
1653                 if ((result = s_MolActionSetCellFlexibility(mol, action, &act2)) != 0)
1654                         return result;
1655         } else if (strcmp(action->name, gMolActionAddParameters) == 0) {
1656                 if ((result = s_MolActionAddParameters(mol, action, &act2)) != 0)
1657                         return result;
1658                 needsRebuildMDArena = 1;
1659         } else if (strcmp(action->name, gMolActionDeleteParameters) == 0) {
1660                 if ((result = s_MolActionDeleteParameters(mol, action, &act2)) != 0)
1661                         return result;
1662                 needsRebuildMDArena = 1;
1663         } else if (strcmp(action->name, gMolActionNone) == 0) {
1664                 /*  Do nothing  */
1665         } else {
1666                 fprintf(stderr, "Internal error: unknown action name %s\n", action->name);
1667                 return -1;
1668         }
1669         if (act2 != NULL) {
1670                 MolActionCallback_registerUndo(mol, act2);
1671                 MolActionRelease(act2);
1672         }
1673
1674         if (needsSymmetryAmendment) {
1675                 MolAction *act3;
1676                 act3 = MolActionNew(gMolActionAmendBySymmetry, ig, NULL);
1677                 if (ig != NULL)
1678                         IntGroupRelease(ig);
1679                 act2 = NULL;
1680                 result = s_MolActionAmendBySymmetry(mol, act3, &act2);
1681                 MolActionRelease(act3);
1682                 if (result != 0)
1683                         return result;
1684                 if (act2 != NULL) {
1685                         MolActionCallback_registerUndo(mol, act2);
1686                         MolActionRelease(act2);
1687                 }
1688         }
1689
1690         if (needsRebuildMDArena) {
1691                 mol->needsMDRebuild = 1;
1692         }
1693
1694         MoleculeCallback_notifyModification(mol, 0);
1695         return 0;
1696 }
1697
1698 int
1699 MolActionCreateAndPerform(Molecule *mol, const char *name, ...)
1700 {
1701         va_list ap;
1702         MolAction *action;
1703         va_start(ap, name);
1704         action = MolActionNewArgv(name, ap);
1705         va_end(ap);
1706         if (action != NULL) {
1707                 int result = MolActionPerform(mol, action);
1708                 MolActionRelease(action);
1709                 return result;
1710         } else return -1;
1711 }
1712
1713 void
1714 MolActionAlertRubyIsRunning(void)
1715 {
1716         MyAppCallback_errorMessageBox("Cannot perform operation (Ruby script is running background)");
1717 }