{"action":"opened","issue":{"url":"https://api.github.com/repos/mdtraj/mdtraj/issues/908","labels_url":"https://api.github.com/repos/mdtraj/mdtraj/issues/908/labels{/name}","comments_url":"https://api.github.com/repos/mdtraj/mdtraj/issues/908/comments","events_url":"https://api.github.com/repos/mdtraj/mdtraj/issues/908/events","html_url":"https://github.com/mdtraj/mdtraj/issues/908","id":99451095,"number":908,"title":"Distance calculation for triclinic boxes in DCD with PBC","user":{"login":"ferchault","id":2736606,"avatar_url":"https://avatars.githubusercontent.com/u/2736606?v=3","gravatar_id":"","url":"https://api.github.com/users/ferchault","html_url":"https://github.com/ferchault","followers_url":"https://api.github.com/users/ferchault/followers","following_url":"https://api.github.com/users/ferchault/following{/other_user}","gists_url":"https://api.github.com/users/ferchault/gists{/gist_id}","starred_url":"https://api.github.com/users/ferchault/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/ferchault/subscriptions","organizations_url":"https://api.github.com/users/ferchault/orgs","repos_url":"https://api.github.com/users/ferchault/repos","events_url":"https://api.github.com/users/ferchault/events{/privacy}","received_events_url":"https://api.github.com/users/ferchault/received_events","type":"User","site_admin":false},"labels":[],"state":"open","locked":false,"assignee":null,"milestone":null,"comments":0,"created_at":"2015-08-06T14:37:29Z","updated_at":"2015-08-06T14:37:29Z","closed_at":null,"body":"I think there may be a problem with the calculation of distances between atoms while including the periodic boundary conditions. I suspect that this comes from a different convention how the H matrix is defined in mdtraj. I found this with a triclinic simulation box but any box where H^T != H should be affected.\r\n\r\nFor easier testing, all files (including the sample DCD file) can be found [here](http://guido.vonrudorff.de/dir/mdtraj-pbc/)\r\n\r\nSample code that generates the DCD file to inspect with vmd (uses MDAnalysis, sorry)\r\n```python\r\n#!/usr/bin/env python\r\n\r\n# generate DCD file, can be verified with vmd\r\nimport MDAnalysis as mda\r\nimport math\r\nimport numpy as np\r\nnumatoms = 2\r\n\r\n# build unit cell and positions\r\nts = mda.coordinates.base.Timestep(numatoms)\r\nts._pos[0] = [0.9, 0.5, 0.1]\r\nts._pos[1] = [0.1, 0.5, 0.9]\r\nts.dimensions = [1, math.sqrt(2), 1, 60, 60, 45]\r\n\r\n# calculate H matrix and its inverse\r\ndef abc_to_hmatrix(a, b, c, alpha, beta, gamma, degrees=True):\r\n        if degrees:\r\n                alpha, beta, gamma = map(math.radians, (alpha, beta, gamma))\r\n        result = np.zeros((3, 3))\r\n\r\n        a = np.array((a, 0, 0))\r\n        b = b * np.array((math.cos(gamma), math.sin(gamma), 0))\r\n        bracket = (math.cos(alpha) - math.cos(beta) * math.cos(gamma)) / math.sin(gamma)\r\n        c = c * np.array((math.cos(beta), bracket, math.sqrt(math.sin(beta) ** 2 - bracket ** 2)))\r\n\r\n        result[:, 0] = a\r\n        result[:, 1] = b\r\n        result[:, 2] = c\r\n\r\n        return result\r\nhmat = abc_to_hmatrix(*ts.dimensions)\r\nhinv = np.linalg.inv(hmat)\r\n\r\n# calculate distance according to Tuckerman\r\ns_1 = np.dot(hinv, ts._pos[0])\r\ns_2 = np.dot(hinv, ts._pos[1])\r\ns_12 = s_1-s_2\r\ns_12 -= np.round(s_12)\r\nprint 'Expected distance by Tuckerman  ', np.linalg.norm(np.dot(hmat, s_12))\r\n\r\n# calculate distance by enumeration\r\ndist = 100000\r\nfor i in range(-1, 2):\r\n\tfor j in range(-1, 2):\r\n\t\tfor k in range(-1, 2):\r\n\t\t\tdist = min(np.linalg.norm(ts._pos[0]+i*hmat[:, 0]+j*hmat[:,1]+k*hmat[:,2]-ts._pos[1]), dist)\r\nprint 'Expected distance by enumeration', dist\r\n\r\n# calculate distance with transposed Hmat\r\ns_1 = np.dot(hinv.T, ts._pos[0])\r\ns_2 = np.dot(hinv.T, ts._pos[1])\r\ns_12 = s_1-s_2\r\ns_12 -= np.round(s_12)\r\nprint 'Wrong distance by Tuckerman     ', np.linalg.norm(np.dot(hmat.T, s_12))\r\n\r\nwriter = mda.coordinates.DCD.DCDWriter('mdtraj-pbc.dcd', numatoms)\r\nwriter.write_next_timestep(ts)\r\nwriter.close()\r\n```\r\n\r\nThis file produces the sample DCD file and this output:\r\n```\r\nExpected distance by Tuckerman   0.366832008509\r\nExpected distance by enumeration 0.366832008509\r\nWrong distance by Tuckerman      0.321717867657\r\n```\r\n\r\nHere, \"wrong distance\" means the distance calculated with a transposed H matrix. Tuckerman's \"Statistical Mechanics\" that is mentioned in the mdtraj source uses the H matrix with the cell vectors in columns (see formula 5.6.2 on page 225), which gives the 0.37 distance that matches the distance measured in VMD or calculated by iterating over the 27 boxes.\r\n\r\nNow using mdtraj:\r\n```python\r\n#!/usr/bin/env python\r\n\r\nimport mdtraj\r\nimport numpy as np\r\n\r\n# open\r\ntraj = mdtraj.load('mdtraj-pbc.dcd', top='mdtraj-pbc.psf')\r\n\r\nprint 'Actual mdtraj distance          ', mdtraj.compute_distances(traj, ((0, 1), ), periodic=True, opt=False)[0][0]\r\n```\r\n\r\nThis prints \r\n```\r\nActual mdtraj distance           0.0321718\r\n```\r\nwhich is identical (up to units) to the \"wrong distance\" mentioned above which comes from following the mdtraj algorithm with an H matrix with the cell vectors in rows. Can anybody confirm this issue? For loading the DCD file generated from the first script, you will need the PSF file:\r\n\r\n```\r\nPSF\r\n\r\n       1 !NTITLE\r\n REMARKS original generated structure x-plor psf file\r\n\r\n       2 !NATOM\r\n       1 UNK  1    LIG  C1   C      0.000000        0.0000           0\r\n       2 UNK  1    LIG  C2   C      0.000000        0.0000          0\r\n\r\n       0 !NBOND: bonds\r\n\r\n\r\n       0 !NTHETA: angles\r\n\r\n\r\n       0 !NPHI: dihedrals\r\n\r\n\r\n       0 !NIMPHI: impropers\r\n\r\n\r\n       0 !NDON: donors\r\n\r\n\r\n       0 !NACC: acceptors\r\n\r\n\r\n       0 !NNB\r\n\r\n\r\n\r\n       1       0 !NGRP\r\n       0       0       0\r\n\r\n```"}}