11
$\begingroup$

I want to make a python script that will load an xyz file. From the xyz parameters, I need to find the distance between atoms, angle and dihedral between atoms.

The file xyz have this structure:

 60  Buckminsterfullerene  (C60 Bucky Ball)
 1  C      0.56182991    1.03720708   -3.34153745     2     2     6    16
 2  C     -0.02867841   -0.20922332   -3.53733326     2     1     3    15
 3  C      0.68868170   -1.41687735   -3.17419275     2     2     4     7
 4  C      1.96798784   -1.33001731   -2.62971513     2     3     5     8
 5  C      2.58298268   -0.03190141   -2.42580025     2     4     6    12
 6  C      1.89418484    1.12766893   -2.77448204     2     1     5    11
 7  C     -0.26911535   -2.36054907   -2.62920302     2     3     9    21
 8  C      2.34254586   -2.18322716   -1.51766993     2     4    10    19
 9  C      0.09052931   -3.17978756   -1.56143494     2     7    10    43
10  C      1.42288437   -3.08932569   -0.99437947     2     8     9    36
11  C      1.93147043    2.28439304   -1.89955106     2     6    14    20
12  C      3.33762840   -0.08283139   -1.18772874     2     5    13    19
13  C      3.37342937    1.02783673   -0.34763393     2    12    14    32
14  C      2.65606922    2.23549084   -0.71077448     2    11    13    26
15  C     -1.42982841   -0.40652370   -3.21677674     2     2    17    21
16  C     -0.22432522    2.13802273   -2.81706602     2     1    18    20
17  C     -2.18468203    0.65046195   -2.71318765     2    15    18    27
18  C     -1.56968716    1.94857800   -2.50927274     2    16    17    22
19  C      3.18903027   -1.41242386   -0.62647337     2     8    12    34
20  C      0.62215914    2.90882602   -1.92586946     2    11    16    23
21  C     -1.57842655   -1.73611615   -2.65552140     2     7    15    29
22  C     -2.12435266    2.52208091   -1.29751971     2    18    24    28
23  C      0.08957811    3.45949446   -0.76236333     2    20    24    25
24  C     -1.31157175    3.26219407   -0.44180686     2    22    23    37
25  C      0.84422386    3.40856444    0.47570821     2    23    26    42
26  C      2.10140370    2.80899380    0.50097868     2    14    25    31
27  C     -3.11943529    0.42168510   -1.62746090     2    17    28    30
28  C     -3.08214969    1.57840920   -0.75252997     2    22    27    39
29  C     -2.47596183   -1.95578408   -1.61302381     2    21    30    44
30  C     -3.26211692   -0.85496848   -1.08855240     2    27    29    40
31  C      2.47596171    1.95578403    1.61302382     2    26    32    46
32  C      3.26211681    0.85496841    1.08855240     2    13    31    33
33  C      3.11943527   -0.42168513    1.62746090     2    32    34    57
34  C      3.08214970   -1.57840928    0.75252996     2    19    33    35
35  C      2.12435273   -2.52208091    1.29751964     2    34    36    56
36  C      1.31157182   -3.26219401    0.44180679     2    10    35    47
37  C     -1.42288432    3.08932572    0.99437953     2    24    38    42
38  C     -2.34254574    2.18322716    1.51766996     2    37    39    52
39  C     -3.18903022    1.41242378    0.62647331     2    28    38    41
40  C     -3.37342947   -1.02783679    0.34763388     2    30    41    48
41  C     -3.33762837    0.08283136    1.18772868     2    39    40    49
42  C     -0.09052935    3.17978761    1.56143497     2    25    37    45
43  C     -0.84422391   -3.40856443   -0.47570819     2     9    44    47
44  C     -2.10140382   -2.80899386   -0.50097868     2    29    43    48
45  C      0.26911532    2.36054910    2.62920311     2    42    46    54
46  C      1.57842653    1.73611617    2.65552148     2    31    45    60
47  C     -0.08957818   -3.45949441    0.76236329     2    36    43    53
48  C     -2.65606930   -2.23549092    0.71077442     2    40    44    50
49  C     -2.58298263    0.03190137    2.42580019     2    41    51    52
50  C     -1.93147042   -2.28439306    1.89955101     2    48    51    53
51  C     -1.89418482   -1.12766893    2.77448197     2    49    50    58
52  C     -1.96798776    1.33001737    2.62971509     2    38    49    54
53  C     -0.62215921   -2.90882598    1.92586939     2    47    50    55
54  C     -0.68868164    1.41687744    3.17419279     2    45    52    59
55  C      0.22432525   -2.13802261    2.81706605     2    53    56    58
56  C      1.56968719   -1.94857793    2.50927279     2    35    55    57
57  C      2.18468202   -0.65046198    2.71318771     2    33    56    60
58  C     -0.56182981   -1.03720706    3.34153744     2    51    55    59
59  C      0.02867849    0.20922332    3.53733334     2    54    58    60
60  C      1.42982839    0.40652370    3.21677684     2    46    57    59

Does anyone know how to do?

$\endgroup$
3
  • 3
    $\begingroup$ The file that you have posted is not in xyz format. Take a look here en.wikipedia.org/wiki/XYZ_file_format This is what most people and programs do use as xyz format. $\endgroup$ – Hans Wurst Apr 20 at 19:51
  • 2
    $\begingroup$ Start with a smaller molecule in first place. This gives you the possibility to check and compute manually (or, calculator) the distances within reasonable effort to check if the algorithm works well. Start with methane, proceed with formic acid, acteone, DMF, ethyl acetate; because given all the symmetry relationships, C60 does not offer an easy discern of the atoms considered for the distance measurements. On the other hand, though, given the symmetry $I_h$ , you may reduce the number of measurements needed by a lot. $\endgroup$ – Buttonwood Apr 21 at 8:49
  • 2
    $\begingroup$ Your question about PyGSM is a separate one, please ask it in a new post! Otherwise, you can keep editing your questions and the existing answers will have to keep changing, which is unfair for the volunteers answering your questions. At least, asking another question in a new post, tends to also give you more points, so it's not so bad after all. $\endgroup$ – Nike Dattani Apr 21 at 19:28
12
$\begingroup$

Easiest free tool to do this (no download necessary):

What you are describing, is the conversion from XYZ format to ZMAT format. How to do this, has been asked here on the Chemistry Stack Exchange, and my answer was to use this free tool in which you just copy and paste the 3 spatial coordinates along with the name of the element, for each row of the XYZ file. I have done it for you and provided the ZMAT file at the bottom of this answer.

A python code to do this (requires downloading the code):

There's many such tools available for free online, including this Python program.

If for some reason you want to code this in Python yourself:

If you don't want to use one of the above free tools, and don't even want to look at the how it was done in the free Python code I mentioned, then still the procedure is doable (though a bit tedious).

Distances: Calculate the Euclidean distances between each pair of coordinates in your molecule. Your molecule is $\ce{C60}$ so you have ${{60}\choose{2}}= 1770$ distances between atoms, which will be calculated as:

$$ D_{ij} = \sqrt{(x_i-x_j)^2 + (y_i -y_j)^2 - (z_i - z_j)^2}\tag{1} $$

for each of the 1770 pairs $(i,j)$ in your XYZ file. Because of symmetry (i.e. you know in $\ce{C60}$ that some bond lengths are expected to be exactly the same) the number of calculations that your Python program has to do, can be greatly reduced, but calculating 1770 distances will take so little time in Python that you don't need to worry about this if computation time is your concern.

Angles within planes: For every three atoms, there's three angles, which you can obtain using the cosine law.

Dihedral angles: For every four atoms, there's dihedral angles between the planes formed by each possible set of three atoms, and there's again a formula to obtain these angles.

Using geometry to reduce the number of necessary calculations: For each subset of three atoms, you only need to calculate 2 bond lengths and 1 angle (for example), not all 3 bond lengths and all 3 angles. Similar rules exist for larger numbers of atoms, and that is why the ZMAT format for $n$ atoms only requires $(n-1)$ distances, $(n-2)$ angles within planes, and $(n-3)$ dihedral angles.

Conversion of your XYZ to ZMAT using the online tool I recommended:

This is the result I get from my recommended solution, which is once again this online tool:

    60
C   1
C   1 1.393
C   2 1.451  1 120.000
C   3 1.393  2 120.000  1   0.0
C   4 1.451  3 120.000  2 360.0
C   5 1.393  4 120.000  3 360.0
C   3 1.451  2 108.000  1 142.6
C   4 1.451  3 120.000  2 138.2
C   7 1.393  3 120.000  2 217.4
C   8 1.393  4 120.000  3 360.0
C   6 1.451  5 120.000  4 138.2
C   5 1.451  4 108.000  3 142.6
C  12 1.393  5 120.000  4 217.4
C  11 1.393  6 120.000  5 360.0
C   2 1.451  1 120.000  3 138.2
C   1 1.451  2 120.000  3 221.8
C  15 1.393  2 120.000  1 360.0
C  16 1.393  1 120.000  2 360.0
C   8 1.451  4 108.000  3 217.4
C  16 1.451  1 108.000  2 142.6
C   7 1.451  3 108.000  2   0.0
C  18 1.451 16 120.000  1 138.2
C  20 1.393 16 120.000  1 217.4
C  22 1.393 18 120.000 16 360.0
C  23 1.451 20 120.000 16 138.2
C  25 1.393 23 120.000 20 360.0
C  17 1.451 15 120.000  2 221.8
C  27 1.451 17 108.000 15 142.6
C  21 1.393  7 120.000  3 217.4
C  27 1.393 17 120.000 15 360.0
C  26 1.451 25 120.000 23 221.8
C  13 1.451 12 120.000  5 138.2
C  32 1.393 13 120.000 12 360.0
C  19 1.393  8 120.000  4 142.6
C  34 1.451 19 120.000  8   0.0
C  35 1.393 34 120.000 19 360.0
C  24 1.451 22 120.000 18 221.8
C  37 1.393 24 120.000 22   0.0
C  28 1.393 27 120.000 17 217.4
C  30 1.451 27 120.000 17 138.2
C  40 1.393 30 120.000 27 360.0
C  37 1.451 24 108.000 22 142.6
C   9 1.451  7 120.000  3 138.2
C  43 1.393  9 120.000  7   0.0
C  42 1.393 37 120.000 24 217.4
C  31 1.393 26 120.000 25 360.0
C  43 1.451  9 108.000  7 217.4
C  44 1.451 43 120.000  9 221.8
C  41 1.451 40 120.000 30 221.8
C  48 1.393 44 120.000 43 360.0
C  49 1.393 41 120.000 40   0.0
C  38 1.451 37 120.000 24 138.2
C  47 1.393 43 120.000  9 142.6
C  52 1.393 38 120.000 37   0.0
C  53 1.451 47 120.000 43 221.8
C  55 1.393 53 120.000 47 360.0
C  56 1.451 55 120.000 53 138.2
C  55 1.451 53 108.000 47 142.6
C  58 1.393 55 120.000 53 217.4
C  57 1.393 56 120.000 55   0.0

To get this, copy and paste the XYZ file (in the correct format!) and set the input file type to "xyz" and output file type to "fh -- Fenske Hall Z-Matrix format". The tool seems to be a bit flaky for 60 atoms, but with some trial and error by removing lines and adding them back, I managed to get it to work.

$\endgroup$
9
  • $\begingroup$ thank you for your time but neighteir of yours sugestions isn't working for my xyz file $\endgroup$ – Jose Pinto Apr 20 at 18:41
  • 2
    $\begingroup$ That's because you have to remove all those extra numbers and put it in the proper format. Also, can you please copy and paste the file into a code block in your question, rather than doing a screenshot of it? Questions in general can get closed if they don't follow our guidelines, and using text rather than images for something like this, is one of our guidelines. $\endgroup$ – Nike Dattani Apr 20 at 19:02
  • 1
    $\begingroup$ @JosePinto It works. I've added the result to my answer. $\endgroup$ – Nike Dattani Apr 20 at 19:39
  • 2
    $\begingroup$ @Buttonwood you're correct, I should have put "n choose 2" not "n^2". I've fixed the answer now. $\endgroup$ – Nike Dattani Apr 22 at 20:54
  • 1
    $\begingroup$ @TristanMaxson I was the one that added the Python tag. Further more you are incorrect that I have given "extra steps" to calculate the distances. The copying and pasting of the XYZ file into the online tool I suggested, is by far the fastest way to get the distances. Faster than installing ASE and figuring out how to use it (as you said in the comment to the other answer). But I didn't want to just give a "cheap" answer, so I also gave an example of a Python tool that does it, and also told the user how the math is actually done (without any extra steps!). $\endgroup$ – Nike Dattani Apr 23 at 0:45
10
$\begingroup$

ASE

The ASE library has an atom object with built-in get_angle, get_dihedral and get_distance methods that do just that.

$\endgroup$
4
  • 2
    $\begingroup$ Can you by chance add some extra details about how exactly to do this? Maybe provide an example script? I fear the asker may not be able to work it out from this and I do not wish to just answer again as "ASE has it" with a script. $\endgroup$ – Tristan Maxson Apr 20 at 17:22
  • $\begingroup$ I agree with @TristanMaxson. I gave you +1 but this is a bit more of a "comment" than an "answer". $\endgroup$ – Nike Dattani Apr 20 at 17:43
  • $\begingroup$ On a similar note, RDKit can also do the same thing $\endgroup$ – user430 Apr 22 at 7:22
  • $\begingroup$ @Ezze why don't you add RDKit as an answer? There's now been a few suggestions written as answers. $\endgroup$ – Nike Dattani Apr 24 at 14:14
6
$\begingroup$

Open Babel

First, the command-line version: obabel c60.xyz -oreport

Yields:

  • All interatomic distances
  • All bond angles
  • All connected torsions
  • Also, partial charges

e.g.

INTERATOMIC DISTANCES

              C   1      C   2      C   3      C   4      C   5      C   6
              ------------------------------------------------------------------
   C   1    0.0000 
   C   2    1.4534     0.0000 
   C   3    2.3516     1.4534     0.0000 
   C   4    2.3516     2.3516     1.4534     0.0000 
   C   5    1.4534     2.3516     2.3516     1.4534     0.0000 
   C   6    2.4673     3.5759     3.5759     2.4673     1.3955     0.0000 
   C   7    3.6940     4.5105     4.1143     2.8488     2.4673     1.4534 

But if you want to do it in Python, you can easily do this through the openbabel or pybel modules. I posted a script we've used for basic machine learning representations

Key calls:

  • for bond in ob.OBMolBondIter(mol.OBMol): # Loop through bonds
  • for pair in ob.OBMolPairIter(mol.OBMol): # Loop through non-bonded pairs of atoms
  • for angle in ob.OBMolAngleIter(mol.OBMol): # Loop through angles
  • for torsion in ob.OBMolTorsionIter(mol.OBMol): # Loop through torsions

You could use a slightly different loop if you didn't want to separate bonded and non-bonded atoms, e.g.

for idx1 in range(1, mol.OBMol.NumAtoms() + 1): # OB counts atoms from 1
    a = mol.OBMol.GetAtom(idx1)
    for idx2 in range(idx1 + 1, mol.OBMol.NumAtoms() + 1):
        b = mol.OBMol.GetAtom(idx2)
        # calculate distance, whatever
$\endgroup$
1
  • 2
    $\begingroup$ And of course, if you want the z-matrix form, you just request it: obabel c60.xyz -ofh $\endgroup$ – Geoff Hutchison Apr 24 at 4:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.