VMD-L Mailing List

From: Norbert Straeter (strater_at_bbz.uni-leipzig.de)
Date: Mon Aug 29 2022 - 11:00:01 CDT

Hi Josh,

thanks for your comments, this was quite helpful.

Indeed the simulations were run using Amber and the periodic boundary
box appears to be specified as a CRYST card in the pdb file (as I did
not run these simulations myself, I do not know for sure right now):

CRYST1   80.078   79.692   82.368  90.00  90.00 90.00               1
ATOM      1  N   GLU     1     -32.373   6.361 -10.893  1.00
0.00           N
ATOM      2  H1  GLU     1     -33.300   6.763 -10.894  1.00
0.00           H

When importing this structure as a pdb file and the dcd trajectory, the
periodic box ended up as:

(VMD) 17 % pbc get
{79.934036 79.504913 82.228027 -3.577219 2.599454 -3.765918}

After setting the box with the pbc set command, the problems with the
measure bond and measure angle commands disappeared.

There may be a better way to read the pdb and dcd files to avoid this
problem. Changing to the newest VMD alpha release did not make any
difference.

Best regards,

Norbert

Am 29.08.2022 um 15:49 schrieb Vermaas, Josh:

> 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
>
>