VMD-L Mailing List

From: Patrick Charchar (patrick.charchar_at_rmit.edu.au)
Date: Tue May 07 2019 - 20:16:32 CDT

Hi Josh,

Thanks for that! I am not very confident with compiling from source, but I guess I will have to do some reading and give this a try if there is no other alternatives.

Best,
Patrick

From: Vermaas, Joshua <Joshua.Vermaas_at_nrel.gov>
Sent: Wednesday, 8 May 2019 10:18 AM
To: vmd-l_at_ks.uiuc.edu; Patrick Charchar <patrick.charchar_at_rmit.edu.au>
Subject: RE: measure rmsd with non-serial atom ordering between two selections

Addendum to my previous post. its just my version of VMD that does that, since I had run into this once myself. This is my edited version of the method in TclMeasure.C that does take the order argument:

//////////////////////////////////////////////
// measure rmsd_qcp $sel1 $sel2 [weight <weights>]
static int vmd_measure_rmsd_qcp(VMDApp *app, int argc, Tcl_Obj * const objv[], Tcl_Interp *interp) {
  int i;
  if (argc !=3 && argc != 5) {
    Tcl_WrongNumArgs(interp, 2, objv-1,
      (char *)"<sel1> <sel2> [weight <weights>] [order <order>]");
    return TCL_ERROR;
  }
  // get the selections
  AtomSel *sel1 = tcl_commands_get_sel(interp, Tcl_GetStringFromObj(objv[1],NULL));
  AtomSel *sel2 = tcl_commands_get_sel(interp, Tcl_GetStringFromObj(objv[2],NULL));
  if (!sel1 || !sel2) {
    Tcl_AppendResult(interp, "measure rmsd: no atom selection", NULL);
    return TCL_ERROR;
  }
  if (sel1->selected != sel2->selected) {
    Tcl_AppendResult(interp, "measure rmsd: selections must have the same number of atoms", NULL);
    return TCL_ERROR;
  }
  if (!sel1->selected) {
    Tcl_AppendResult(interp, "measure rmsd: no atoms selected", NULL);
    return TCL_ERROR;
  }

  Tcl_Obj *weightvar = NULL;
  Tcl_Obj *ordervar = NULL;
  for (i=3; i<argc-1; i+=2) {
    const char *opt = Tcl_GetStringFromObj(objv[i], NULL);
    if (!strcmp(opt, "weight")) {
      weightvar = objv[i+1];
    } else if (!strcmp(opt, "order")) {
      ordervar = objv[i+1];
    } else {
      Tcl_AppendResult(interp, "measure rmsd: unknown option '", opt, "'",
          NULL);
      return TCL_ERROR;
    }
  }
  float *weight = new float[sel1->selected];
  {
    int ret_val;
    if (argc == 3) {
      ret_val = tcl_get_weights(interp, app, sel1, NULL, weight);
    } else {
      ret_val = tcl_get_weights(interp, app, sel1, weightvar, weight);
    }
    if (ret_val < 0) {
      Tcl_AppendResult(interp, "measure rmsd: ", measure_error(ret_val),
      NULL);
      delete [] weight;
      return TCL_ERROR;
    }
  }
  int ret_val;
  // compute the rmsd
  {
    float rmsd = 0;
    const float *x = sel1->coordinates(app->moleculeList);
    const float *y = sel2->coordinates(app->moleculeList);
    if (!x || !y) {
      delete [] weight;
      return TCL_ERROR;
    }
    if (ordervar != NULL) {
      // get the atom order
      int *order = new int[sel1->selected];
      ret_val = tcl_get_orders(interp, sel1->selected, ordervar, order);
      if (ret_val < 0) {
        Tcl_AppendResult(interp, "measure rmsd: ", measure_error(ret_val), NULL);
        delete [] order;
        delete [] weight;
        return TCL_ERROR;
      }
      int i;
      int j = 0;
      float *tmp = new float[3L*sel1->num_atoms];
      for (i=sel1->firstsel; i<sel1->lastsel; i++) {
        if (sel1->on[i]) {
          long ind = 3L * i;
          long idx = 3L * order[j]; // order indices are 0-based
          tmp[ind ] = y[idx ];
          tmp[ind + 1] = y[idx + 1];
          tmp[ind + 2] = y[idx + 2];
          j++;
        }
      }
      ret_val = measure_rmsd_qcp(sel1, sel2, sel1->selected, x, tmp, weight, &rmsd);
      delete[] tmp;
      delete[] order;
    } else {
      ret_val = measure_rmsd_qcp(sel1, sel2, sel1->selected, x, y, weight, &rmsd);
    }
    delete [] weight;
    if (ret_val < 0) {
      Tcl_AppendResult(interp, "measure rmsd: ", measure_error(ret_val),
      NULL);
      return TCL_ERROR;
    }
    Tcl_SetObjResult(interp, Tcl_NewDoubleObj(rmsd));
  }
  return TCL_OK;
}

In principle you'd replace this in the source tree and recompile VMD and you'd be set. However that can at times be non-trivial, so I'm not sure how much that helps you with your current predicament.

-Josh

On 2019-05-07 18:11:09-06:00 Vermaas, Joshua wrote:
Hi Patrick,

In VMD 1.9.3, there is a sparsely documented alternative algorithm that *does* take the order parameter. Try "measure rmsd_qcp $sel $ref order $orderlist"

-Josh

On 2019-05-07 18:00:52-06:00 owner-vmd-l_at_ks.uiuc.edu<mailto:owner-vmd-l_at_ks.uiuc.edu> wrote:
Hi all,

Is anyone aware if an atom order parameter been implemented for the "measure rmsd" command, or if there is a workaround for trajectories that does not involve manually reordering atoms in the file?

Perhaps there is a way to get the RMSD directly from "measure fit"?

Previous related posts:
http://www.ks.uiuc.edu/Research/vmd/mailing_list/vmd-l/10201.htmldbe4417b4d4bdb24b51316b%7C1&sdata=bngrQo%2FrcQCtuExPI3QxWeLb9d2peYGRYLfOPX3%2B95M%3D&reserved=0>
https://www.ks.uiuc.edu/Research/vmd/mailing_list/vmd-l/17002.htmlcdbe4417b4d4bdb24b51316b%7C1&sdata=gwj1ILY1mpIotgDTKk4LRCpbEf%2BQRf4Arb%2FRMonvnwA%3D&reserved=0>

Thanks,
Patrick