VMD-L Mailing List

From: Vermaas, Josh (vermaasj_at_msu.edu)
Date: Mon Aug 29 2022 - 08:49:37 CDT

Hi Norbert,

I don't have concrete advice, other than to perhaps check that the PBC dimensions or periodic image transformations are doing what you want them to do. At the root of measure bonds is this little function:

int calculate_bond(MoleculeList *mlist, int *molid, int *atmid, float *value) {

  // get coords to calculate distance
  int ret_val;
  float pos1[3], pos2[3];
  if ((ret_val=normal_atom_coord(mlist->mol_from_id(molid[0]), atmid[0], pos1))<0)
    return ret_val;
  if ((ret_val=normal_atom_coord(mlist->mol_from_id(molid[1]), atmid[1], pos2))<0)
    return ret_val;

  vec_sub(pos2, pos2, pos1);
  *value = norm(pos2);

  return MEASURE_NOERR;
}

The most likely way for this to go wrong is if somehow normal_atom_coord gives you a NAN or an INF or something. This is the definition for normal_atom_coord:

// for the given Molecule, find the UNTRANSFORMED coords for the given atom
// return Molecule pointer if successful, NULL otherwise.
int normal_atom_coord(Molecule *mol, int a, float *pos) {
  Timestep *now;

  int cell[3];
  memset(cell, 0, 3L*sizeof(int));

  // get the molecule pointer, and get the coords for the current timestep
  int ret_val = check_mol(mol, a);
  if (ret_val<0)
    return ret_val;

  if ((now = mol->current())) {
    memcpy((void *)pos, (void *)(now->pos + 3L*a), 3L*sizeof(float));

    // Apply periodic image transformation before returning
    Matrix4 mat;
    now->get_transform_from_cell(cell, mat);
    mat.multpoint3d(pos, pos);

    return MEASURE_NOERR;
  }

  // if here, error (i.e. molecule contains no frames)
  return MEASURE_ERR_NOFRAMES;
}

The place where this is dangerous is if the transformation matrix is undefined somehow, since you might be multiplying with uninitialized memory, which can be all sorts of things. Again, this is all pure speculation, since at least for .dcd files, the PBC dimensions are always included in the file, so the transformation should totally work. But if you are working with a different file type, it is at least conceivable that the PBC info is wrong in 1.9.3. I recall there being some issues with some Amber file formats not recognizing PBC info correctly back then. Perhaps one of the 1.9.4 alphas works?

-Josh

On 8/29/22, 9:33 AM, "owner-vmd-l_at_ks.uiuc.edu on behalf of Norbert Straeter" <owner-vmd-l_at_ks.uiuc.edu on behalf of strater_at_bbz.uni-leipzig.de> wrote:

    Dear all,

    I am analyzing a trajectory with VMD to retrieve an interatomic distance
    and angle to analyze a hydrogen bonding interaction. For some frames the
    value "-NaN" is returned for both measure commands. The problem appears
    using VMD 1.9.3. on a Windows notebook and on a Linux computer. There is
    nothing unusual that I can spot in the corresponding snapshots of the
    frames and I can calculate the distance and angle in vmd using the
    coordinates. So in principle, there is a workaround, but I wonder if
    anybody has an idea concerning this problem, as it would be preferable
    to use the measure bond and measure angle commands directly.

    This is the output, for simplicity only using the measure bond command:

>Main< (with_protonated_his) 20 % set frame 1
    1
>Main< (with_protonated_his) 21 % set atom1 [[atomselect top "name OH
    and resid 267"] get index]
    4340
>Main< (with_protonated_his) 22 % set atom3 [[atomselect top "name O4
    and resid 272"] get index]
    4442
>Main< (with_protonated_his) 23 % set at1xyz [lindex [[atomselect top
    "name OH and resid 267" frame $frame] get {x y z}] 0]
    6.038471221923828 -8.148454666137695 -1.7956215143203735
>Main< (with_protonated_his) 24 % set at3xyz [lindex [[atomselect top
    "name O4 and resid 272" frame $frame] get {x y z}] 0]
    5.896842956542969 -6.978078842163086 1.9930213689804077
>Main< (with_protonated_his) 25 % set distance [measure bond [list
    $atom1 $atom3] frame $frame]
    3.967827320098877
>Main< (with_protonated_his) 26 % set frame 2
    2
>Main< (with_protonated_his) 27 % set at1xyz [lindex [[atomselect top
    "name OH and resid 267" frame $frame] get {x y z}] 0]
    6.15404748916626 -8.125819206237793 -1.7144218683242798
>Main< (with_protonated_his) 28 % set at3xyz [lindex [[atomselect top
    "name O4 and resid 272" frame $frame] get {x y z}] 0]
    5.544281482696533 -6.590084075927734 1.5558415651321411
>Main< (with_protonated_his) 29 % set distance [measure bond [list
    $atom1 $atom3] frame $frame]
    -NaN
>Main< (with_protonated_his) 30 %

    The error occurs in about 10 % of the 20000 frames, but I cannot see any
    pattern,

    Has anybody encountered this problem or has suggestions?

    Best regards,

    Norbert