• ## Outreach

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