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 *gMolActionReplaceTables   = "replaceTables:III"; */
  41 const char *gMolActionTranslateAtoms  = "translateAtoms:vG";
  42 const char *gMolActionRotateAtoms     = "rotate:vdvG";
  43 const char *gMolActionTransformAtoms  = "transform:tG";
  44 const char *gMolActionSetAtomPositions = "atomPosition:GV";
  45 const char *gMolActionInsertFrames    = "insertFrames:GVV";
  46 const char *gMolActionRemoveFrames    = "removeFrames:G";
  47 const char *gMolActionSetSelection    = "selection:G";
  48 const char *gMolActionChangeResidueNumber = "changeResSeq:Gi";
  49 const char *gMolActionChangeResidueNumberForUndo = "changeResSeqForUndo:GIi";
  50 const char *gMolActionChangeResidueNames = "changeResNames:IC";
  51 const char *gMolActionOffsetResidueNumbers = "offsetResSeq:Gii";
  52 const char *gMolActionChangeNumberOfResidues = "changeNres:i";
  53 const char *gMolActionRenumberAtoms    = "renumberAtoms:I";
  54 const char *gMolActionExpandBySymmetry = "expandSym:Giiii";
  55 const char *gMolActionDeleteSymmetryOperation = "deleteSymop";
  56 const char *gMolActionAddSymmetryOperation = "addSymop:t";
  57 const char *gMolActionSetCell         = "setCell:Di";
  58 const char *gMolActionSetBox          = "setBox:vvvvi";
  59 const char *gMolActionClearBox        = "clearBox";
  60 /*const char *gMolActionSetParameterAttribute = "setParameterAttr:iirri"; */
  61 const char *gMolActionAddParameters   = "addParameters:iGU";
  62 const char *gMolActionDeleteParameters = "deleteParameters:iG";
  63 const char *gMolActionCartesianToXtal  = "cartesianToXtal";
  64 const char *gMolActionXtalToCartesian  = "xtalToCartesian";
  65
  66 /*  A Ruby array to retain objects used in MolActionArg  */
  67 static VALUE sMolActionArgValues = Qfalse;
  68
  69 /*  Action arguments  */
  70 /*  (Simple types)  i: Int, d: double, s: string, v: Vector, t: Transform, u: UnionPar
  71  (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
  72  (Complex types) M: Molecule, G: IntGroup, A: Atom
  73  (Ruby value)    b: Ruby boolean, r: Ruby object, R: an array of Ruby object (a Ruby array)
  74  (Return value from Ruby script) i, d, s, v, t + 0x80 */
  75 typedef struct MolActionArg {
  76         Byte type;
  77         union {
  78                 Int ival;
  79                 double dval;
  80                 struct {  /*  Array type: the size of the element is defined by the argument type  */
  81                         Int nitems; /* Number of items */
  82                         void *ptr;
  83                 } arval;
  84                 struct IntGroup *igval; /* The record is retained but not duplicated */
  85                 struct Molecule *mval; /* The record is retained but not duplicated */
  86                 struct Atom *aval; /* The value is duplicated, so any pointer can be passed */
  87                 VALUE vval;        /* The value is retained in sMolActionArgValues  */
  88                 void *pval;        /* Store address for the return value; only meaningful in performing Ruby script  */
  89         } u;
  90 } MolActionArg;
  91
  92 /*  For debug  */
  93 /* static int sMolActionCount = 0; */
  94
  95 /*  Create a new MolAction with the given arguments.
  96  *  Note: If the argument is an array type, it is represented by an integer (number
  97  *  of items) followed by the pointer.  */
  98 MolAction *
  99 MolActionNewArgv(const char *name, va_list ap)
 100 {
 101         MolAction *action;
 102         const char *p;
 103
 104         action = (MolAction *)malloc(sizeof(MolAction));
 105         if (action == NULL)
 106                 goto low_memory;
 107         memset(action, 0, sizeof(MolAction));
 108         action->refCount = 1;
 109         action->frame = -1;
 110         action->name = strdup(name);
 111
 112         /*  For debug  */
 113 /*      fprintf(stderr, "MolAction %p created, count = %d\n", action, ++sMolActionCount); */
 114
 115         /*  Handle arguments  */
 116         p = strchr(name, ':');
 117         if (p == NULL)
 118                 return action;
 119         p++;
 120
 121         while (*p) {
 122                 MolActionArg arg;
 123                 memset(&arg, 0, sizeof(MolActionArg));
 124                 arg.type = *p++;
 125                 switch (arg.type) {
 126                         case 'i':
 127                                 arg.u.ival = va_arg(ap, Int);
 128                                 break;
 129                         case 'd':
 130                                 arg.u.dval = va_arg(ap, double);
 131                                 break;
 132                         case 'I':
 133                         case 'D':
 134                         case 'C':
 135                         case 's':
 136                         case 'V':
 137                         case 'v':
 138                         case 'T':
 139                         case 't':
 140                         case 'U':
 141                         case 'u': {
 142                                 /*  Array type of some kind  */
 143                                 void *ptr;
 144                                 Int nitems, itemsize, allocsize;
 145                                 if (arg.type != 's' && arg.type != 'v' && arg.type != 't' && arg.type != 'u')
 146                                         nitems = va_arg(ap, Int);
 147                                 ptr = va_arg(ap, void *);
 148                                 switch (arg.type) {
 149                                         case 'I': itemsize = sizeof(Int); break;
 150                                         case 'D': itemsize = sizeof(double); break;
 151                                         case 'C':
 152                                         case 's': itemsize = sizeof(char); break;
 153                                         case 'V':
 154                                         case 'v': itemsize = sizeof(Vector); break;
 155                                         case 'T':
 156                                         case 't': itemsize = sizeof(Transform); break;
 157                                         case 'U':
 158                                         case 'u': itemsize = sizeof(UnionPar); break;
 159                                 }
 160                                 if (arg.type == 's')
 161                                         nitems = (ptr == NULL ? 0 : strlen((char *)ptr));
 162                                 else if (arg.type == 'v' || arg.type == 't' || arg.type == 'u')
 163                                         nitems = 1;
 164                                 arg.u.arval.nitems = nitems;
 165                                 allocsize = itemsize * nitems + (arg.type == 's'); /* +1 for string type */
 166                                 if (allocsize == 0)
 167                                         arg.u.arval.ptr == NULL;
 168                                 else {
 169                                         arg.u.arval.ptr = calloc(1, allocsize);
 170                                         if (arg.u.arval.ptr == NULL)
 171                                                 goto low_memory;
 172                                         memmove(arg.u.arval.ptr, ptr, itemsize * nitems);
 173                                 }
 174                                 break;
 175                         }
 176 /*                      case 's':
 177                                 arg.u.sval = va_arg(ap, char *);
 178                                 arg.u.sval = strdup(arg.u.sval);
 179                                 break; */
 180 /*                      case 'v': {
 181                                 Vector *vp = va_arg(ap, Vector *);
 182                                 arg.u.vval = (Vector *)malloc(sizeof(Vector));
 183                                 if (arg.u.vval == NULL)
 184                                         goto low_memory;
 185                                 *(arg.u.vval) = *vp;
 186                                 break;
 187                         } */
 188                         case 'G':
 189                                 arg.u.igval = va_arg(ap, IntGroup *);
 190                                 if (arg.u.igval != NULL)
 191                                         IntGroupRetain(arg.u.igval);
 192                                 break;
 193                         case 'M':
 194                                 arg.u.mval = va_arg(ap, Molecule *);
 195                                 MoleculeRetain(arg.u.mval);
 196                                 break;
 197                         case 'A': {
 198                                 Atom *aval = va_arg(ap, Atom *);
 199                                 arg.u.aval = (Atom *)malloc(gSizeOfAtomRecord);
 200                                 if (arg.u.aval == NULL)
 201                                         goto low_memory;
 202                                 AtomDuplicate(arg.u.aval, aval);
 203                                 break;
 204                         }
 205                         case 'r':
 206                         case 'R': {
 207                                 int rtype;
 208                                 arg.u.vval = va_arg(ap, VALUE);
 209                                 rtype = TYPE(arg.u.vval);
 210                                 if (rtype != T_NIL && rtype != T_FIXNUM && rtype != T_SYMBOL && rtype != T_TRUE && rtype != T_FALSE) {
 211                                         if (sMolActionArgValues == Qfalse) {
 212                                                 sMolActionArgValues = rb_ary_new();
 213                                                 rb_global_variable(&sMolActionArgValues);
 214                                         }
 215                                         rb_ary_push(sMolActionArgValues, arg.u.vval);
 216                                 }
 217                                 break;
 218                         }
 219                         case 'b':
 220                                 arg.u.vval = (va_arg(ap, Int) ? Qtrue : Qfalse);
 221                                 break;
 222                         case ';':
 223                                 if (*p == 'i' || *p == 'd' || *p == 's' || *p == 'v' || *p == 't' || *p == 'r') {
 224                                         arg.type = (*p | 0x80);
 225                                         p++;
 226                                         arg.u.pval = va_arg(ap, void *);
 227                                 } else {
 228                                         fprintf(stderr, "Internal error: unknown return-type \'%c\' in NewMolAction()", *p);
 229                                 }
 230                                 break;
 231                         default:
 232                                 fprintf(stderr, "Internal error: unknown argument type \'%c\' in NewMolAction()", arg.type);
 233                                 break;
 234                 }
 235                 if (AssignArray(&action->args, &action->nargs, sizeof(arg), action->nargs, &arg) == NULL)
 236                         goto low_memory;
 237         }
 238         return action;
 239
 240   low_memory:
 241         if (action != NULL)
 242                 free(action);
 243         Panic("Low memory during creation of action record");
 244         return NULL;  /* Not reached */
 245 }
 246
 247 MolAction *
 248 MolActionNew(const char *name, ...)
 249 {
 250         va_list ap;
 251         MolAction *retval;
 252         va_start(ap, name);
 253         retval = MolActionNewArgv(name, ap);
 254         va_end(ap);
 255         return retval;
 256 }
 257
 258 MolAction *
 259 MolActionRetain(MolAction *action)
 260 {
 261         if (action != NULL)
 262                 action->refCount++;
 263         return action;
 264 }
 265
 266 void
 267 MolActionRelease(MolAction *action)
 268 {
 269         int i;
 270         if (action == NULL)
 271                 return;
 272         if (--action->refCount > 0)
 273                 return;
 274         if (action->name != NULL)
 275                 free((void *)action->name);
 276         for (i = 0; i < action->nargs; i++) {
 277                 MolActionArg *argp = action->args + i;
 278                 switch (argp->type) {
 279                         case 'i':
 280                         case 'd':
 281                         case 'b':
 282                                 break;
 283                         case 'I':
 284                         case 'D':
 285                         case 'C':
 286                         case 's':
 287                         case 'V':
 288                         case 'v':
 289                         case 'T':
 290                         case 't':
 291                         case 'U':
 292                         case 'u':
 293                                 if (argp->u.arval.ptr != NULL)
 294                                         free(argp->u.arval.ptr);
 295                                 break;
 296
 297                 /*      case 'a':
 298                                 if (argp->u.iaval != NULL)
 299                                         free(argp->u.iaval);
 300                                 break;
 301                         case 's':
 302                                 if (argp->u.sval != NULL)
 303                                         free(argp->u.sval);
 304                                 break;
 305                         case 'v':
 306                                 if (argp->u.vval != NULL)
 307                                         free(argp->u.vval);
 308                                 break; */
 309                         case 'G':
 310                                 if (argp->u.igval != NULL)
 311                                         IntGroupRelease(argp->u.igval);
 312                                 break;
 313                         case 'M':
 314                                 MoleculeRelease(argp->u.mval);
 315                                 break;
 316                         case 'A':
 317                                 if (argp->u.aval != NULL) {
 318                                         AtomClean(argp->u.aval);
 319                                         free(argp->u.aval);
 320                                 }
 321                                 break;
 322                         case 'r':
 323                         case 'R': {
 324                                 int rtype = TYPE(argp->u.vval);
 325                                 if (rtype != T_NIL && rtype != T_FIXNUM && rtype != T_SYMBOL && rtype != T_TRUE && rtype != T_FALSE) {
 326                                         /*  Look for the same value in sMolActionArgValues and remove it  */
 327                                         if (sMolActionArgValues != Qfalse) {
 328                                                 VALUE *ptr = RARRAY_PTR(sMolActionArgValues);
 329                                                 int n = RARRAY_LEN(sMolActionArgValues);
 330                                                 while (--n >= 0) {
 331                                                         if (ptr[n] == argp->u.vval) {
 332                                                                 rb_ary_delete_at(sMolActionArgValues, n);
 333                                                                 break;
 334                                                         }
 335                                                 }
 336                                         }
 337                                 }
 338                                 break;
 339                         }
 340                         default:
 341                                 break;
 342                 }
 343         }
 344         if (action->args != NULL)
 345                 free(action->args);
 346         free(action);
 347
 348         /*  For debug  */
 349 /*      fprintf(stderr, "MolAction %p disposed, count = %d\n", action, --sMolActionCount);      */
 350 }
 351
 352 void
 353 MolActionSetFrame(MolAction *action, int frame)
 354 {
 355         if (action != NULL)
 356                 action->frame = frame;
 357 }
 358
 359 typedef struct MolRubyActionInfo {
 360         Molecule *mol;
 361         MolAction *action;
 362         VALUE receiver;
 363         ID method_id;
 364         VALUE ary;
 365         char enable_interrupt;
 366         char return_type;
 367         void *return_ptr;
 368 } MolRubyActionInfo;
 369
 370 static VALUE
 371 s_MolActionToRubyArguments(VALUE vinfo)
 372 {
 373         MolRubyActionInfo *info = (MolRubyActionInfo *)vinfo;
 374         VALUE val;
 375         int i, argc;
 376         char *s, *p;
 377
 378         argc = info->action->nargs - 1;
 379         info->return_type = 0;
 380         info->return_ptr = NULL;
 381         info->ary = rb_ary_new2(argc);
 382         for (i = 1; i <= argc; i++) {
 383                 MolActionArg *argp = &(info->action->args[i]);
 384                 if (argp->type & 0x80) {
 385                         /*  Return value specified  */
 386                         info->return_type = (argp->type & 0x7f);
 387                         info->return_ptr = argp->u.pval;
 388                         break;  /*  This should be the last argument  */
 389                 }
 390                 switch (argp->type) {
 391                         case 'i':
 392                                 val = INT2NUM(argp->u.ival);
 393                                 break;
 394                         case 'd':
 395                                 val = rb_float_new(argp->u.dval);
 396                                 break;
 397                         case 's':
 398                         case 'C':
 399                                 val = rb_str_new((char *)argp->u.arval.ptr, argp->u.arval.nitems);
 400                                 break;
 401                         case 'v':
 402                                 val = ValueFromVector((Vector *)argp->u.arval.ptr);
 403                                 break;
 404                         case 't':
 405                                 val = ValueFromTransform((Transform *)argp->u.arval.ptr);
 406                                 break;
 407                         case 'I':
 408                         case 'D':
 409                         case 'V':
 410                         case 'T': {
 411                                 int j;
 412                                 val = rb_ary_new();
 413                                 for (j = 0; j < argp->u.arval.nitems; j++) {
 414                                         VALUE val1;
 415                                         void *ptr = argp->u.arval.ptr;
 416                                         switch (argp->type) {
 417                                                 case 'I': val1 = INT2NUM(((Int *)ptr)[j]); break;
 418                                                 case 'D': val1 = rb_float_new(((double *)ptr)[j]); break;
 419                                                 case 'V': val1 = ValueFromVector(&((Vector *)ptr)[j]); break;
 420                                                 case 'T': val1 = ValueFromTransform(&((Transform *)ptr)[j]); break;
 421                                         }
 422                                         rb_ary_push(val, val1);
 423                                 }
 424                                 break;
 425                         }
 426                         case 'G':
 427                                 val = ValueFromIntGroup(info->action->args[i].u.igval);
 428                                 break;
 429                         case 'M':
 430                                 val = ValueFromMolecule(info->action->args[i].u.mval);
 431                                 break;
 432                         case 'r':
 433                         case 'b':
 434                                 val = (VALUE)info->action->args[i].u.vval;
 435                                 break;
 436                         case 'R':
 437                                 rb_ary_concat(info->ary, (VALUE)info->action->args[i].u.vval);
 438                                 continue;  /*  Skip rb_ary_push()  */
 439                         default:
 440                                 rb_raise(rb_eArgError, "Unsupported argument type '%c'", info->action->args[i].type);
 441                 }
 442                 rb_ary_push(info->ary, val);
 443         }
 444
 445         if (info->mol != NULL) {
 446                 info->receiver = ValueFromMolecule(info->mol);
 447         } else {
 448                 info->receiver = rb_cMolecule;  /*  Assume class method  */
 449         }
 450
 451         s = info->action->args[0].u.arval.ptr;
 452         for (p = s; *p != 0; p++) {
 453                 if (isalpha(*p) || *p == '_' || (p > s && isdigit(*p)) || (p[1] == 0 && (*p == '?' || *p == '!')))
 454                         continue;
 455                 break;
 456         }
 457
 458         if (*p == 0) {
 459                 /*  Can be a method name  */
 460                 info->method_id = rb_intern(s);
 461         } else {
 462                 /*  Cannot be a method name: try to create a proc object and call it  */
 463                 VALUE procstr = rb_str_new2(s);
 464                 VALUE proc = rb_obj_instance_eval(1, &procstr, info->receiver);
 465                 info->receiver = proc;
 466                 info->method_id = rb_intern("call");
 467         }
 468
 469         return Qnil;
 470 }
 471
 472 static void
 473 s_MolActionStoreReturnValue(MolRubyActionInfo *info, VALUE val)
 474 {
 475         if (info->return_type != 0 && info->return_ptr != NULL) {
 476                 switch (info->return_type) {
 477                         case 'i': *((Int *)(info->return_ptr)) = NUM2INT(rb_Integer(val)); break;
 478                         case 'd': *((Double *)(info->return_ptr)) = NUM2DBL(rb_Float(val)); break;
 479                         case 's': *((char **)(info->return_ptr)) = strdup(StringValuePtr(val)); break;
 480                         case 'v': VectorFromValue(val, (Vector *)(info->return_ptr)); break;
 481                         case 't': TransformFromValue(val, (Transform *)(info->return_ptr)); break;
 482                         case 'r': *((RubyValue *)(info->return_ptr)) = (RubyValue)val; break;
 483                 }
 484         }
 485 }
 486
 487 static VALUE
 488 s_MolActionPerformRubyScript(VALUE vinfo)
 489 {
 490 /*      VALUE save_interrupt_flag; */
 491         VALUE val;
 492         MolRubyActionInfo *info = (MolRubyActionInfo *)vinfo;
 493 /*      if (info->enable_interrupt)
 494                 save_interrupt_flag = Ruby_SetInterruptFlag(Qtrue); */
 495         val = rb_funcall2(info->receiver, info->method_id, RARRAY_LEN(info->ary), RARRAY_PTR(info->ary));
 496         s_MolActionStoreReturnValue(info, val);
 497 /*      if (info->enable_interrupt)
 498                 Ruby_SetInterruptFlag(save_interrupt_flag); */
 499         return val;
 500 }
 501
 502 static void
 503 sMolActionUpdateSelectionAndParameterNumbering(Molecule *mol, const IntGroup *ig, int is_insert)
 504 {
 505         int i, j, n, old_natoms;
 506         IntGroup *sel, *orig_atoms, *ig3;
 507         UnionPar *up;
 508         MolAction *act;
 509
 510         if (ig == NULL || (n = IntGroupGetCount(ig)) == 0)
 511                 return;
 512
 513         /*  Insert: orig_atoms = set of original atoms in new indices  */
 514         /*  Delete: orig_atoms = set of remaining atoms in old indices  */
 515         orig_atoms = IntGroupNew();
 516         old_natoms = (is_insert ? mol->natoms - n : mol->natoms + n);
 517         ig3 = IntGroupNewWithPoints(0, (is_insert ? mol->natoms : old_natoms), -1);
 518         IntGroupXor(ig, ig3, orig_atoms);
 519         IntGroupRelease(ig3);
 520
 521         /*  Update selection  */
 522         sel = MoleculeGetSelection(mol);
 523         if (sel != NULL && IntGroupGetCount(sel) > 0) {
 524                 IntGroup *sel2 = IntGroupNew();
 525                 if (is_insert)
 526                         IntGroupConvolute(sel, orig_atoms, sel2);  /*  Renumber  */
 527                 else {
 528                         IntGroup *selx = IntGroupNew();
 529                         IntGroupIntersect(sel, orig_atoms, selx);
 530                         IntGroupDeconvolute(selx, orig_atoms, sel2);
 531                         IntGroupRelease(selx);
 532                 }
 533                 IntGroupRetain(sel);  /*  To avoid being released during MoleculeSetSelection()  */
 534                 MoleculeSetSelection(mol, sel2);
 535                 if (IntGroupGetCount(sel) != IntGroupGetCount(sel2)) {
 536                         /*  Register undo action  */
 537                         act = MolActionNew(gMolActionSetSelection, sel);
 538                         MolActionCallback_registerUndo(mol, act);
 539                         MolActionRelease(act);
 540                 }
 541                 IntGroupRelease(sel2);
 542                 IntGroupRelease(sel);
 543         }
 544
 545         /*  Update parameters  */
 546         if (mol->par != NULL) {
 547                 Int *ip = (Int *)malloc(sizeof(Int) * old_natoms);
 548                 if (is_insert) {
 549                         for (i = 0; i < old_natoms; i++)
 550                                 ip[i] = IntGroupGetNthPoint(orig_atoms, i);
 551                 } else {
 552                         for (i = 0; i < old_natoms; i++)
 553                                 ip[i] = IntGroupLookupPoint(orig_atoms, i);
 554                 }
 555                 for (i = kFirstParType; i <= kLastParType; i++) {
 556                         UnionPar usave;
 557                         UnionPar *upary = NULL;
 558                         Int count_upary = 0;
 559                         if (!is_insert)
 560                                 ig3 = IntGroupNew();
 561                         for (j = 0; (up = ParameterGetUnionParFromTypeAndIndex(mol->par, i, j)) != NULL; j++) {
 562                                 usave = *up;
 563                                 if (ParameterRenumberAtoms(i, up, old_natoms, ip) && !is_insert) {
 564                                         IntGroupAdd(ig3, j, 1);  /*  This parameter is to be restored on undo  */
 565                                         AssignArray(&upary, &count_upary, sizeof(UnionPar), count_upary, &usave);
 566                                 }
 567                         }
 568                         if (count_upary > 0) {
 569                                 /*  Register undo for modifying parameters  */
 570                                 act = MolActionNew(gMolActionAddParameters, i, ig3, count_upary, upary);
 571                                 MolActionCallback_registerUndo(mol, act);
 572                                 MolActionRelease(act);
 573                                 act = MolActionNew(gMolActionDeleteParameters, i, ig3);
 574                                 MolActionCallback_registerUndo(mol, act);
 575                                 MolActionRelease(act);
 576                                 free(upary);
 577                         }
 578                         if (!is_insert)
 579                                 IntGroupRelease(ig3);
 580                 }
 581                 free(ip);
 582         }
 583         IntGroupRelease(orig_atoms);
 584 }
 585
 586 int
 587 MolActionPerform(Molecule *mol, MolAction *action)
 588 {
 589         int n1, result, natoms;
 590         Molecule *mol2;
 591         IntGroup *ig;
 592         MolAction *act2 = NULL;
 593         int needsSymmetryAmendment = 0;
 594         int needsRebuildMDArena = 0;
 595         Int *ip;
 596         Vector *vp, v;
 597         n1 = 0;
 598
 599         if (action == NULL)
 600                 return -1;
 601
 602         if (action->frame >= 0 && mol->cframe != action->frame)
 603                 MoleculeSelectFrame(mol, action->frame, 1);
 604
 605         /*  Ruby script execution  */
 606         if (strncmp(action->name, kMolActionPerformScript, strlen(kMolActionPerformScript)) == 0) {
 607                 MolRubyActionInfo info;
 608                 memset(&info, 0, sizeof(info));
 609                 if (gMolbyIsCheckingInterrupt) {
 610                         MolActionAlertRubyIsRunning();
 611                         return -1;
 612                 }
 613         /*      gMolbyBacktrace = Qnil; */
 614                 info.mol = mol; /*  May be NULL  */
 615                 info.action = action;
 616                 gMolbyRunLevel++;
 617                 rb_protect(s_MolActionToRubyArguments, (VALUE)&info, &result);
 618                 gMolbyRunLevel--;
 619                 if (result == 0) {
 620                         VALUE save_interrupt;
 621                         MyAppCallback_beginUndoGrouping();
 622                 /*      info.enable_interrupt = 1; */
 623                         save_interrupt = Ruby_SetInterruptFlag(Qtrue);
 624                         gMolbyRunLevel++;
 625                         rb_protect(s_MolActionPerformRubyScript, (VALUE)&info, &result);
 626                         gMolbyRunLevel--;
 627                         Ruby_SetInterruptFlag(save_interrupt);
 628                         MyAppCallback_endUndoGrouping();
 629                 }
 630                 if (result != 0) {
 631                         Molby_showError(result);
 632                 }
 633                 MyAppCallback_hideProgressPanel();  /*  In case when the progress panel is still onscreen */\v
 634                 return (result == 0 ? 0 : -1);
 635         }
 636
 637         if (mol == NULL)
 638                 return -1;
 639         natoms = mol->natoms;
 640
 641         if (strcmp(action->name, gMolActionAddAnAtom) == 0) {
 642                 n1 = MoleculeCreateAnAtom(mol, action->args[0].u.aval, action->args[1].u.ival);
 643                 if ((ip = action->args[2].u.pval) != NULL)
 644                         *ip = n1;
 645                 if (n1 < 0)
 646                         return -1;
 647                 ig = IntGroupNewWithPoints(n1, 1, -1);
 648
 649                 /*  Update selection  */
 650                 sMolActionUpdateSelectionAndParameterNumbering(mol, ig, 1);
 651
 652                 act2 = MolActionNew(gMolActionDeleteAnAtom, n1);
 653                 needsRebuildMDArena = 1;
 654         } else if (strcmp(action->name, gMolActionMergeMolecule) == 0) {
 655                 IntGroup *ig2;
 656                 Molecule *mol2;
 657                 ig = action->args[1].u.igval;
 658                 mol2 = action->args[0].u.mval;
 659                 if (ig != NULL) {
 660                         ig2 = ig;
 661                         IntGroupRetain(ig2);
 662                 } else {
 663                         ig2 = IntGroupNew();
 664                         IntGroupAdd(ig2, mol->natoms, mol2->natoms);
 665                 }
 666                 /*  Calculate the offset for the residue number  */
 667                 {
 668                         int minreg, maxreg;
 669                         Atom *ap;
 670                         maxreg = -1;
 671                         for (n1 = 0, ap = mol->atoms; n1 < mol->natoms; n1++, ap = ATOM_NEXT(ap)) {
 672                                 if (ap->resSeq > maxreg)
 673                                         maxreg = ap->resSeq;
 674                         }
 675                         minreg = ATOMS_MAX_NUMBER;
 676                         for (n1 = 0, ap = mol2->atoms; n1 < mol2->natoms; n1++, ap = ATOM_NEXT(ap)) {
 677                                 if (ap->resSeq < minreg)
 678                                         minreg = ap->resSeq;
 679                         }
 680                         if (maxreg < 0)
 681                                 maxreg = 0;
 682                         if (minreg == ATOMS_MAX_NUMBER)
 683                                 minreg = 0;
 684                         if (maxreg >= minreg)
 685                                 n1 = maxreg - minreg + 1;
 686                         else n1 = 0;
 687                 }
 688                 if ((result = MoleculeMerge(mol, mol2, ig, n1)) != 0)
 689                         return result;
 690
 691                 sMolActionUpdateSelectionAndParameterNumbering(mol, ig, 1);
 692
 693                 act2 = MolActionNew(gMolActionUnmergeMolecule, ig2);
 694                 IntGroupRelease(ig2);
 695                 needsRebuildMDArena = 1;
 696         } else if (strcmp(action->name, gMolActionDeleteAnAtom) == 0 || (strcmp(action->name, gMolActionUnmergeMolecule) == 0 && (n1 = 1))) {
 697                 IntGroup *bg;
 698                 if (n1 == 0) {
 699                         ig = IntGroupNew();
 700                         if (ig == NULL || IntGroupAdd(ig, action->args[0].u.ival, 1) != 0)
 701                                 return -1;
 702                 } else {
 703                         ig = action->args[0].u.igval;
 704                         IntGroupRetain(ig);
 705                 }
 706                 /*  Search bonds crossing the molecule border  */
 707                 bg = MoleculeSearchBondsAcrossAtomGroup(mol, ig);
 708                 ip = NULL;
 709                 if (bg != NULL) {
 710                         n1 = IntGroupGetCount(bg);
 711                         if (n1 > 0) {
 712                                 IntGroupIterator iter;
 713                                 Int i, idx;
 714                                 ip = (Int *)calloc(sizeof(Int), n1 * 2 + 1);
 715                                 if (ip == NULL)
 716                                         return -1;
 717                                 IntGroupIteratorInit(bg, &iter);
 718                                 idx = 0;
 719                                 while ((i = IntGroupIteratorNext(&iter)) >= 0) {
 720                                         /*  The atoms at the border  */
 721                                         ip[idx++] = mol->bonds[i * 2];
 722                                         ip[idx++] = mol->bonds[i * 2 + 1];
 723                                 }
 724                                 IntGroupIteratorRelease(&iter);
 725                         }
 726                         IntGroupRelease(bg);
 727                 }
 728                 /*  Unmerge molecule  */
 729                 if ((result = MoleculeUnmerge(mol, &mol2, ig, 0)) != 0) {
 730                         if (ip != NULL)
 731                                 free(ip);
 732                         return result;
 733                 }
 734
 735                 sMolActionUpdateSelectionAndParameterNumbering(mol, ig, 0);
 736
 737 #if 0
 738                 /*  Update selection  */
 739                 {
 740                         IntGroup *sel = MoleculeGetSelection(mol);
 741                         if (sel != NULL && IntGroupGetCount(sel) > 0) {
 742                                 IntGroup *remain_atoms = IntGroupNew();  /*  Remaining atoms (with original indices)  */
 743                                 IntGroup *sel2 = IntGroupNew();
 744                                 IntGroup *sel3 = IntGroupNew();
 745                                 IntGroupXor(ig, IntGroupNewWithPoints(0, mol->natoms + IntGroupGetCount(ig), -1), remain_atoms);
 746                                 IntGroupIntersect(sel, remain_atoms, sel2);  /*  Atoms to select after deletion  */
 747                                 IntGroupDeconvolute(sel2, remain_atoms, sel3);  /*  Renumber  */
 748                                 if (!IntGroupIsEqual(sel, sel3)) {
 749                                         MoleculeSetSelection(mol, sel3);
 750                                         act2 = MolActionNew(gMolActionSetSelection, sel);
 751                                         MolActionCallback_registerUndo(mol, act2);
 752                                         act2 = NULL;
 753                                 }
 754                                 IntGroupRelease(sel3);
 755                                 IntGroupRelease(sel2);
 756                                 IntGroupRelease(remain_atoms);
 757                         }
 758                 }
 759 #endif
 760                 if (mol2 == NULL)
 761                         act2 = NULL;
 762                 else {
 763                         /*  If there exist bonds crossing the molecule border, then register
 764                             an action to restore them  */
 765                         if (ip != NULL) {
 766                                 act2 = MolActionNew(gMolActionAddBonds, n1 * 2, ip);
 767                                 MolActionCallback_registerUndo(mol, act2);
 768                                 free(ip);
 769                         }
 770                         act2 = MolActionNew(gMolActionMergeMolecule, mol2, ig);
 771                         MoleculeRelease(mol2);
 772                 }
 773                 IntGroupRelease(ig);
 774                 needsRebuildMDArena = 1;
 775         } else if (strcmp(action->name, gMolActionAddBonds) == 0) {
 776                 ip = (Int *)action->args[0].u.arval.ptr;
 777                 n1 = action->args[0].u.arval.nitems / 2;
 778                 if ((result = MoleculeAddBonds(mol, n1, ip)) <= 0)
 779                         return result;
 780                 ip = (Int *)malloc(sizeof(Int) * 2 * result);
 781                 if (ip == NULL)
 782                         return -4;
 783                 memmove(ip, mol->bonds - result * 2, sizeof(Int) * 2 * result);
 784                 act2 = MolActionNew(gMolActionDeleteBonds, result * 2, ip);
 785                 free(ip);
 786                 needsRebuildMDArena = 1;
 787         } else if (strcmp(action->name, gMolActionDeleteBonds) == 0) {
 788                 ip = (Int *)action->args[0].u.arval.ptr;
 789                 n1 = action->args[0].u.arval.nitems / 2;
 790                 if ((result = MoleculeDeleteBonds(mol, n1, ip)) < 0)
 791                         return result;
 792                 act2 = MolActionNew(gMolActionAddBonds, n1 * 2, ip);
 793                 needsRebuildMDArena = 1;
 794         } else if (strcmp(action->name, gMolActionAddAngles) == 0) {
 795                 ip = (Int *)action->args[0].u.arval.ptr;
 796                 n1 = action->args[0].u.arval.nitems / 3;
 797                 ig = action->args[1].u.igval;
 798                 if (ig == NULL)
 799                         ig = IntGroupNewWithPoints(mol->nangles, n1, -1);
 800                 else
 801                         IntGroupRetain(ig);
 802                 if ((result = MoleculeAddAngles(mol, ip, ig)) < 0) {
 803                         IntGroupRelease(ig);
 804                         return result;
 805                 }
 806                 act2 = MolActionNew(gMolActionDeleteAngles, ig);
 807                 IntGroupRelease(ig);
 808                 needsRebuildMDArena = 1;
 809         } else if (strcmp(action->name, gMolActionDeleteAngles) == 0) {
 810                 ig = action->args[0].u.igval;
 811                 n1 = IntGroupGetCount(ig) * 3;
 812                 ip = (Int *)malloc(sizeof(Int) * n1);
 813                 if ((result = MoleculeDeleteAngles(mol, ip, ig)) < 0) {
 814                         free(ip);
 815                         return result;
 816                 }
 817                 act2 = MolActionNew(gMolActionAddAngles, n1, ip, ig);
 818                 free(ip);
 819                 needsRebuildMDArena = 1;
 820         } else if (strcmp(action->name, gMolActionAddDihedrals) == 0) {
 821                 ip = (Int *)action->args[0].u.arval.ptr;
 822                 n1 = action->args[0].u.arval.nitems / 4;
 823                 ig = action->args[1].u.igval;
 824                 if (ig == NULL)
 825                         ig = IntGroupNewWithPoints(mol->ndihedrals, n1, -1);
 826                 else
 827                         IntGroupRetain(ig);
 828                 if ((result = MoleculeAddDihedrals(mol, ip, ig)) < 0) {
 829                         IntGroupRelease(ig);
 830                         return result;
 831                 }
 832                 act2 = MolActionNew(gMolActionDeleteDihedrals, ig);
 833                 IntGroupRelease(ig);
 834                 needsRebuildMDArena = 1;
 835         } else if (strcmp(action->name, gMolActionDeleteDihedrals) == 0) {
 836                 ig = action->args[0].u.igval;
 837                 n1 = IntGroupGetCount(ig) * 4;
 838                 ip = (Int *)malloc(sizeof(Int) * n1);
 839                 if ((result = MoleculeDeleteDihedrals(mol, ip, ig)) < 0) {
 840                         free(ip);
 841                         return result;
 842                 }
 843                 act2 = MolActionNew(gMolActionAddDihedrals, n1, ip, ig);
 844                 free(ip);
 845                 needsRebuildMDArena = 1;
 846         } else if (strcmp(action->name, gMolActionAddImpropers) == 0) {
 847                 ip = (Int *)action->args[0].u.arval.ptr;
 848                 n1 = action->args[0].u.arval.nitems / 4;
 849                 ig = action->args[1].u.igval;
 850                 if (ig == NULL)
 851                         ig = IntGroupNewWithPoints(mol->nimpropers, n1, -1);
 852                 else
 853                         IntGroupRetain(ig);
 854                 if ((result = MoleculeAddImpropers(mol, ip, ig)) < 0) {
 855                         IntGroupRelease(ig);
 856                         return result;
 857                 }
 858                 act2 = MolActionNew(gMolActionDeleteImpropers, ig);
 859                 IntGroupRelease(ig);
 860                 needsRebuildMDArena = 1;
 861         } else if (strcmp(action->name, gMolActionDeleteImpropers) == 0) {
 862                 ig = action->args[0].u.igval;
 863                 n1 = IntGroupGetCount(ig) * 4;
 864                 ip = (Int *)malloc(sizeof(Int) * n1);
 865                 if ((result = MoleculeDeleteImpropers(mol, ip, ig)) < 0) {
 866                         free(ip);
 867                         return result;
 868                 }
 869                 act2 = MolActionNew(gMolActionAddImpropers, n1, ip, ig);
 870                 free(ip);
 871                 needsRebuildMDArena = 1;
 872 /*      } else if (strcmp(action->name, gMolActionReplaceTables) == 0) {
 873                 Int i1, i2, i3;
 874                 Int *ip1, *ip2, *ip3;
 875                 if ((i1 = MoleculeReplaceAllAngles(mol, action->args[0].u.arval.nitems / 3, (Int *)action->args[0].u.arval.ptr, &ip1)) < 0)
 876                         return -1;
 877                 if ((i2 = MoleculeReplaceAllDihedrals(mol, action->args[1].u.arval.nitems / 4, (Int *)action->args[1].u.arval.ptr, &ip2)) < 0)
 878                         return -1;
 879                 if ((i3 = MoleculeReplaceAllImpropers(mol, action->args[2].u.arval.nitems / 4, (Int *)action->args[2].u.arval.ptr, &ip3)) < 0)
 880                         return -1;
 881                 act2 = MolActionNew(gMolActionReplaceTables, i1, ip1, i2, ip2, i3, ip3);
 882                 free(ip1);
 883                 free(ip2);
 884                 free(ip3);
 885                 needsRebuildMDArena = 1; */
 886         } else if (strcmp(action->name, gMolActionTranslateAtoms) == 0) {
 887                 vp = (Vector *)action->args[0].u.arval.ptr;
 888                 ig = action->args[1].u.igval;
 889                 if (vp == NULL)
 890                         return -1;
 891                 MoleculeTranslate(mol, vp, ig);
 892                 VecScale(v, *vp, -1);
 893                 act2 = MolActionNew(gMolActionTranslateAtoms, &v, ig);
 894                 act2->frame = mol->cframe;
 895                 needsSymmetryAmendment = 1;
 896         } else if (strcmp(action->name, gMolActionRotateAtoms) == 0) {
 897                 Vector *vp2;
 898                 Double ang;
 899                 vp = (Vector *)action->args[0].u.arval.ptr;
 900                 ang = action->args[1].u.dval;
 901                 vp2 = (Vector *)action->args[2].u.arval.ptr;
 902                 ig = action->args[3].u.igval;
 903                 MoleculeRotate(mol, vp, ang, vp2, ig);
 904                 act2 = MolActionNew(gMolActionRotateAtoms, vp, -ang, vp2, ig);
 905                 act2->frame = mol->cframe;
 906                 needsSymmetryAmendment = 1;
 907         } else if (strcmp(action->name, gMolActionTransformAtoms) == 0) {
 908                 Transform *trp, tr_inv;
 909                 int atomPositions = 0; /* Save atom positions for undo? */
 910                 trp = (Transform *)action->args[0].u.arval.ptr;
 911                 ig = action->args[1].u.igval;
 912                 if (TransformInvert(tr_inv, *trp) != 0)
 913                         atomPositions = 1;
 914                 else {
 915                         /*  Check if inverse transform is reliable enough  */
 916                         Transform temp;
 917                         int j;
 918                         TransformMul(temp, *trp, tr_inv);
 919                         temp[0] -= 1.0;
 920                         temp[4] -= 1.0;
 921                         temp[8] -= 1.0;
 922                         for (j = 0; j < 12; j++) {
 923                                 if (temp[j] < -1e-4 || temp[j] > 1e-4)
 924                                         break;
 925                         }
 926                         if (j != 12)
 927                                 atomPositions = 1;
 928                 }
 929                 if (atomPositions) {
 930                         IntGroupIterator iter;
 931                         int j, k;
 932                         n1 = IntGroupGetCount(ig);
 933                         vp = (Vector *)malloc(sizeof(Vector) * n1);
 934                         if (vp == NULL)
 935                                 return -1;
 936                         IntGroupIteratorInit(ig, &iter);
 937                         k = 0;
 938                         while ((j = IntGroupIteratorNext(&iter)) >= 0) {
 939                                 vp[k++] = (ATOM_AT_INDEX(mol->atoms, j))->r;
 940                         }
 941                         act2 = MolActionNew(gMolActionSetAtomPositions, ig, n1, vp);
 942                         free(vp);
 943                         IntGroupIteratorRelease(&iter);
 944                 } else {
 945                         act2 = MolActionNew(gMolActionTransformAtoms, &tr_inv, ig);
 946                 }
 947                 MoleculeTransform(mol, *trp, ig);
 948                 act2->frame = mol->cframe;
 949                 needsSymmetryAmendment = 1;
 950         } else if (strcmp(action->name, gMolActionSetAtomPositions) == 0) {
 951                 IntGroupIterator iter;
 952                 int j, k;
 953                 ig = action->args[0].u.igval;
 954                 n1 = IntGroupGetCount(ig);
 955                 vp = (Vector *)malloc(sizeof(Vector) * n1);
 956                 if (vp == NULL)
 957                         return -1;
 958                 IntGroupIteratorInit(ig, &iter);
 959                 k = 0;
 960                 while ((j = IntGroupIteratorNext(&iter)) >= 0) {
 961                         vp[k++] = (ATOM_AT_INDEX(mol->atoms, j))->r;
 962                 }
 963                 act2 = MolActionNew(gMolActionSetAtomPositions, ig, n1, vp);
 964                 free(vp);
 965                 vp = (Vector *)action->args[1].u.arval.ptr;
 966                 IntGroupIteratorReset(&iter);
 967                 k = 0;
 968                 while ((j = IntGroupIteratorNext(&iter)) >= 0) {
 969                         (ATOM_AT_INDEX(mol->atoms, j))->r = vp[k++];
 970                 }
 971                 IntGroupIteratorRelease(&iter);
 972                 act2->frame = mol->cframe;
 973                 needsSymmetryAmendment = 1;
 974         } else if (strcmp(action->name, gMolActionInsertFrames) == 0) {
 975                 int old_nframes, new_nframes;
 976                 Vector *vp2;
 977                 ig = action->args[0].u.igval;
 978                 vp = (Vector *)action->args[1].u.arval.ptr;
 979                 vp2 = (Vector *)action->args[2].u.arval.ptr;
 980                 n1 = IntGroupGetCount(ig);
 981                 if (n1 == 0)
 982                         return 0;  /*  Do nothing  */
 983                 if (vp != NULL && action->args[1].u.arval.nitems != n1 * mol->natoms)
 984                         return -1;  /*  Internal inconsistency  */
 985                 if (vp2 != NULL && action->args[2].u.arval.nitems != n1 * 4)
 986                         return -1;  /*  Internal inconsistency  */
 987                 old_nframes = MoleculeGetNumberOfFrames(mol);
 988                 if (MoleculeInsertFrames(mol, ig, vp, vp2) < 0)
 989                         return -1;  /*  Error  */
 990                 new_nframes = MoleculeGetNumberOfFrames(mol);
 991                 if (old_nframes + n1 < new_nframes) {
 992                         /*  "Extra" frames were automatically inserted because large frame indices were specified  */
 993                         /*  Register undo operation to remove these extra frames  */
 994                         IntGroup *ig2 = IntGroupNewWithPoints(old_nframes, new_nframes - (old_nframes + n1), -1);
 995                         act2 = MolActionNew(gMolActionRemoveFrames, ig2);
 996                         IntGroupRelease(ig2);
 997                         MolActionCallback_registerUndo(mol, act2);
 998                         MolActionRelease(act2);
 999                 }
1000                 act2 = MolActionNew(gMolActionRemoveFrames, ig);
1001                 act2->frame = mol->cframe;
1002         } else if (strcmp(action->name, gMolActionRemoveFrames) == 0) {
1003                 Vector *vp2;
1004                 ig = action->args[0].u.igval;
1005                 n1 = IntGroupGetCount(ig);
1006                 if (n1 == 0)
1007                         return 0;  /*  Do nothing  */
1008                 vp = (Vector *)calloc(sizeof(Vector), n1 * mol->natoms);
1009                 if (mol->cell != NULL && mol->frame_cells != NULL)
1010                         vp2 = (Vector *)calloc(sizeof(Vector) * 4, n1);
1011                 else vp2 = NULL;
1012                 if (MoleculeRemoveFrames(mol, ig, vp, vp2) < 0)
1013                         return -1;  /*  Error  */
1014                 act2 = MolActionNew(gMolActionInsertFrames, ig, n1 * mol->natoms, vp, (vp2 != NULL ? n1 * 4 : 0), vp2);
1015                 act2->frame = mol->cframe;
1016                 free(vp);
1017                 if (vp2 != NULL)
1018                         free(vp2);
1019         } else if (strcmp(action->name, gMolActionSetSelection) == 0) {
1020                 IntGroup *ig2;
1021                 ig2 = MoleculeGetSelection(mol);
1022                 if (ig2 != NULL)
1023                         IntGroupRetain(ig2);  /*  To avoid releasing during MoleculeSetSelection() */
1024                 ig = action->args[0].u.igval;
1025                 MoleculeSetSelection(mol, ig);
1026                 act2 = MolActionNew(gMolActionSetSelection, ig2);
1027                 if (ig2 != NULL)
1028                         IntGroupRelease(ig2);
1029         } else if (strcmp(action->name, gMolActionRenumberAtoms) == 0) {
1030                 Int *ip2 = (Int *)malloc(sizeof(Int) * (mol->natoms + 1));
1031                 if (ip2 == NULL)
1032                         return -1;
1033                 ip = (Int *)action->args[0].u.arval.ptr;
1034                 n1 = action->args[0].u.arval.nitems;
1035         /*      for (n1 = 0; n1 < mol->natoms; n1++) {
1036                         if (ip[n1] < 0)
1037                                 break;
1038                 } */
1039                 result = MoleculeRenumberAtoms(mol, ip, ip2, n1);
1040                 if (result != 0) {
1041                         free(ip2);
1042                         return result;
1043                 }
1044                 act2 = MolActionNew(gMolActionRenumberAtoms, mol->natoms, ip2);
1045                 needsRebuildMDArena = 1;
1046         } else if (strcmp(action->name, gMolActionChangeResidueNumber) == 0 || (strcmp(action->name, gMolActionChangeResidueNumberForUndo) == 0 && (n1 = 1))) {
1047                 IntGroupIterator iter;
1048                 int i, nresidues;
1049                 int forUndo;
1050                 forUndo = n1;
1051                 ig = action->args[0].u.igval;
1052                 n1 = IntGroupGetCount(ig);
1053                 ip = (Int *)calloc(sizeof(Int), n1 + 1);
1054                 IntGroupIteratorInit(ig, &iter);
1055                 i = 0;
1056                 while ((n1 = IntGroupIteratorNext(&iter)) >= 0) {
1057                         ip[i++] = ATOM_AT_INDEX(mol->atoms, n1)->resSeq;
1058                 }
1059                 n1 = i;
1060                 ip[i] = kInvalidIndex;
1061                 nresidues = mol->nresidues;
1062                 if (forUndo) {
1063                         MoleculeChangeResidueNumberWithArray(mol, ig, (Int *)action->args[1].u.arval.ptr);
1064                         MoleculeChangeNumberOfResidues(mol, action->args[2].u.ival);
1065                 } else {
1066                         MoleculeChangeResidueNumber(mol, ig, action->args[1].u.ival);
1067                 }
1068                 act2 = MolActionNew(gMolActionChangeResidueNumberForUndo, ig, n1, ip, nresidues);
1069         } else if (strcmp(action->name, gMolActionChangeResidueNames) == 0) {
1070                 char *new_names, *old_names;
1071                 int i, argc;
1072                 ip = (Int *)action->args[0].u.arval.ptr;
1073                 argc = action->args[0].u.arval.nitems;
1074                 new_names = action->args[1].u.arval.ptr;
1075                 old_names = (char *)malloc(argc * 4 + 1);
1076                 for (i = 0; i < argc; i++) {
1077                         if (ip[i] >= mol->nresidues)
1078                                 old_names[i * 4] = 0;
1079                         else
1080                                 strncpy(old_names + i * 4, mol->residues[ip[i]], 4);
1081                 }
1082                 MoleculeChangeResidueNames(mol, argc, ip, new_names);
1083                 act2 = MolActionNew(gMolActionChangeResidueNames, argc, ip, argc * 4, old_names);
1084                 free(old_names);
1085         } else if (strcmp(action->name, gMolActionOffsetResidueNumbers) == 0) {
1086                 ig = action->args[0].u.igval;
1087                 n1 = mol->nresidues;
1088                 result = MoleculeOffsetResidueNumbers(mol, ig, action->args[1].u.ival, action->args[2].u.ival);
1089                 if (result != 0)
1090                         return result;  /*  The molecule is not modified  */
1091                 act2 = MolActionNew(gMolActionOffsetResidueNumbers, ig, -(action->args[1].u.ival), n1);
1092         } else if (strcmp(action->name, gMolActionChangeNumberOfResidues) == 0) {
1093                 int nresidues = mol->nresidues;
1094                 n1 = action->args[0].u.ival;
1095                 if (n1 < nresidues) {
1096                         /*  The residue names will be lost, so undo must be registered to restore the names  */
1097                         int argc = nresidues - n1;
1098                         char *names = (char *)malloc(4 * argc + 1);
1099                         Int *ip = (Int *)malloc(sizeof(Int) * (argc + 1));
1100                         int i;
1101                         for (i = 0; i < argc; i++) {
1102                                 strncpy(names + i * 4, mol->residues[i + n1], 4);
1103                                 ip[i] = i + n1;
1104                         }
1105                         ip[i] = kInvalidIndex;
1106                         act2 = MolActionNew(gMolActionChangeResidueNames, argc, ip, argc * 4, names);
1107                         MolActionCallback_registerUndo(mol, act2);
1108                         free(ip);
1109                         free(names);
1110                 }
1111                 MoleculeChangeNumberOfResidues(mol, n1);
1112                 act2 = MolActionNew(gMolActionChangeNumberOfResidues, nresidues);
1113         } else if (strcmp(action->name, gMolActionExpandBySymmetry) == 0) {
1114                 Symop symop;
1115                 ig = action->args[0].u.igval;
1116                 symop.dx = action->args[1].u.ival;
1117                 symop.dy = action->args[2].u.ival;
1118                 symop.dz = action->args[3].u.ival;
1119                 symop.sym = action->args[4].u.ival;
1120                 symop.alive = (symop.dx != 0 || symop.dy != 0 || symop.dz != 0 || symop.sym != 0);
1121                 n1 = MoleculeAddExpandedAtoms(mol, symop, ig);
1122                 if (n1 > 0) {
1123                         ig = IntGroupNew();
1124                         IntGroupAdd(ig, mol->natoms - n1, n1);
1125                         act2 = MolActionNew(gMolActionUnmergeMolecule, ig);
1126                         IntGroupRelease(ig);
1127                 } else return n1;
1128                 needsRebuildMDArena = 1;
1129         } else if (strcmp(action->name, gMolActionAddSymmetryOperation) == 0) {
1130                 Transform *trp;
1131                 Transform itr = {1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0};
1132                 trp = (Transform *)action->args[0].u.arval.ptr;
1133                 if (mol->nsyms == 0) {
1134                         for (n1 = 0; n1 < 12; n1++) {
1135                                 if (fabs((*trp)[n1] - itr[n1]) > 1e-8)
1136                                         break;
1137                         }
1138                         if (n1 < 12) {
1139                                 if (AssignArray(&mol->syms, &mol->nsyms, sizeof(Transform), mol->nsyms, &itr) == 0)
1140                                         return -1;
1141                                 act2 = MolActionNew(gMolActionDeleteSymmetryOperation);
1142                                 MolActionCallback_registerUndo(mol, act2);
1143                                 MolActionRelease(act2);
1144                         }
1145                 }
1146                 if (AssignArray(&mol->syms, &mol->nsyms, sizeof(Transform), mol->nsyms, trp) == 0)
1147                                 return -1;
1148                 act2 = MolActionNew(gMolActionDeleteSymmetryOperation);
1149         } else if (strcmp(action->name, gMolActionDeleteSymmetryOperation) == 0) {
1150                 if (mol->nsyms == 0)
1151                         return -1;
1152                 act2 = MolActionNew(gMolActionAddSymmetryOperation, &(mol->syms[mol->nsyms - 1]));
1153                 mol->nsyms--;
1154         /*      if (mol->nsyms == 1)
1155                         mol->nsyms--;  *//*  Remove the identity operation  */
1156                 if (mol->nsyms == 0) {
1157                         free(mol->syms);
1158                         mol->syms = NULL;
1159                 }
1160         } else if (strcmp(action->name, gMolActionSetCell) == 0) {
1161                 double *dp, d[6];
1162                 if (mol->cell == NULL)
1163                         d[0] = 0.0;
1164                 else {
1165                         for (n1 = 0; n1 < 6; n1++)
1166                                 d[n1] = mol->cell->cell[n1];
1167                 }
1168                 n1 = action->args[1].u.ival;
1169                 dp = action->args[0].u.arval.ptr;
1170                 if (action->args[0].u.arval.nitems == 0)
1171                         MoleculeSetCell(mol, 0, 0, 0, 0, 0, 0, n1);
1172                 else
1173                         MoleculeSetCell(mol, dp[0], dp[1], dp[2], dp[3], dp[4], dp[5], n1);
1174                 act2 = MolActionNew(gMolActionSetCell, (d[0] == 0.0 ? 0 : 6), d, n1);
1175         } else if (strcmp(action->name, gMolActionSetBox) == 0) {
1176                 Vector v[4];
1177                 char flags[3];
1178                 if (mol->cell == NULL)
1179                         act2 = MolActionNew(gMolActionClearBox);
1180                 else {
1181                         n1 = ((mol->cell->flags[0] != 0) * 4 + (mol->cell->flags[1] != 0) * 2 + (mol->cell->flags[2] != 0));
1182                         act2 = MolActionNew(gMolActionSetBox, &(mol->cell->axes[0]), &(mol->cell->axes[1]), &(mol->cell->axes[2]), &(mol->cell->origin), n1);
1183                 }
1184                 for (n1 = 0; n1 < 4; n1++)
1185                         v[n1] = *((Vector *)(action->args[n1].u.arval.ptr));
1186                 for (n1 = 0; n1 < 3; n1++)
1187                         flags[n1] = ((action->args[4].u.ival >> (2 - n1)) & 1);
1188                 MoleculeSetPeriodicBox(mol, &v[0], &v[1], &v[2], &v[3], flags);
1189                 needsRebuildMDArena = 1;
1190         } else if (strcmp(action->name, gMolActionClearBox) == 0) {
1191                 if (mol->cell == NULL)
1192                         return 0;  /*  Do nothing  */
1193                 n1 = ((mol->cell->flags[0] != 0) * 4 + (mol->cell->flags[1] != 0) * 2 + (mol->cell->flags[2] != 0));
1194                 act2 = MolActionNew(gMolActionSetBox, &(mol->cell->axes[0]), &(mol->cell->axes[1]), &(mol->cell->axes[2]), &(mol->cell->origin), n1);
1195                 MoleculeSetPeriodicBox(mol, NULL, NULL, NULL, NULL, NULL);
1196                 needsRebuildMDArena = 1;
1197         } else if (strcmp(action->name, gMolActionAddParameters) == 0) {
1198                 UnionPar *up;
1199                 Int parType;
1200                 if (mol->par == NULL)
1201                         return -1;
1202                 parType = action->args[0].u.ival;
1203                 ig = action->args[1].u.igval;
1204                 up = action->args[2].u.arval.ptr;
1205                 ParameterInsert(mol->par, parType, up, ig);
1206                 act2 = MolActionNew(gMolActionDeleteParameters, parType, ig);
1207         } else if (strcmp(action->name, gMolActionDeleteParameters) == 0) {
1208                 UnionPar *up;
1209                 Int parType;
1210                 if (mol->par == NULL)
1211                         return -1;
1212                 parType = action->args[0].u.ival;
1213                 ig = action->args[1].u.igval;
1214                 n1 = IntGroupGetCount(ig);
1215                 up = (UnionPar *)calloc(sizeof(UnionPar), n1);
1216                 ParameterDelete(mol->par, parType, up, ig);
1217                 act2 = MolActionNew(gMolActionAddParameters, parType, ig, n1, up);
1218                 free(up);
1219 /*      } else if (strcmp(action->name, gMolActionCartesianToXtal) == 0 || (strcmp(action->name, gMolActionXtalToCartesian) == 0 && (n1 = 1))) {
1220                 Int i, j;
1221                 Atom *ap;
1222                 if (mol->cell == NULL || (n1 == 0 && mol->is_xtal_coord) || (n1 == 1 && !mol->is_xtal_coord))
1223                         return 0;
1224                 if (n1 == 0) {
1225                         for (i = 0, ap = mol->atoms; i < mol->natoms; i++, ap = ATOM_NEXT(ap)) {
1226                                 TransformVec(&(ap->r), mol->cell->tr, &(ap->r));
1227                                 if (ap->nframes > 0) {
1228                                         for (j = 0; j < ap->nframes; j++)
1229                                                 TransformVec(&(ap->frames[j]), mol->cell->tr, &(ap->frames[j]));
1230                                 }
1231                         }
1232                         mol->is_xtal_coord = 1;
1233                         act2 = MolActionNew(gMolActionXtalToCartesian);
1234                 } else {
1235                         for (i = 0, ap = mol->atoms; i < mol->natoms; i++, ap = ATOM_NEXT(ap))
1236                                 TransformVec(&(ap->r), mol->cell->rtr, &(ap->r));
1237                         if (ap->nframes > 0) {
1238                                 for (j = 0; j < ap->nframes; j++)
1239                                         TransformVec(&(ap->frames[j]), mol->cell->rtr, &(ap->frames[j]));
1240                         }
1241                         mol->is_xtal_coord = 0;
1242                         act2 = MolActionNew(gMolActionCartesianToXtal);
1243                 } */
1244         } else {
1245                 fprintf(stderr, "Internal error: unknown action name %s\n", action->name);
1246                 return -1;
1247         }
1248         if (act2 != NULL) {
1249                 MolActionCallback_registerUndo(mol, act2);
1250                 MolActionRelease(act2);
1251         }
1252         if (needsSymmetryAmendment) {
1253                 //  ig should be valid
1254                 IntGroup *ig2;
1255                 Vector *vp2;
1256                 n1 = MoleculeAmendBySymmetry(mol, ig, &ig2, &vp2);
1257                 if (n1 > 0) {
1258                         act2 = MolActionNew(gMolActionSetAtomPositions, ig2, n1, vp2);
1259                         act2->frame = mol->cframe;
1260                         MolActionCallback_registerUndo(mol, act2);
1261                         MolActionRelease(act2);
1262                         free(vp2);
1263                         IntGroupRelease(ig2);
1264                 }
1265         }
1266
1267         if (needsRebuildMDArena) {
1268                 mol->needsMDRebuild = 1;
1269         }
1270
1271         MoleculeCallback_notifyModification(mol, 0);
1272         return 0;
1273 }
1274
1275 int
1276 MolActionCreateAndPerform(Molecule *mol, const char *name, ...)
1277 {
1278         va_list ap;
1279         MolAction *action;
1280         va_start(ap, name);
1281         action = MolActionNewArgv(name, ap);
1282         va_end(ap);
1283         if (action != NULL) {
1284                 int result = MolActionPerform(mol, action);
1285                 MolActionRelease(action);
1286                 return result;
1287         } else return -1;
1288 }
1289
1290 void
1291 MolActionAlertRubyIsRunning(void)
1292 {
1293         MyAppCallback_errorMessageBox("Cannot perform operation (Ruby script is running background)");
1294 }