19
$\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 this?

$\endgroup$
7
  • 4
    $\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, 2021 at 19:51
  • 3
    $\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, 2021 at 8:49
  • 3
    $\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$ Apr 21, 2021 at 19:28
  • 1
    $\begingroup$ @Jose do you know what the last 4 columns in your file represent? I mean, the columns that come after the XYZ coordinates? I'm curious to know what file format this is. $\endgroup$ Apr 11, 2023 at 1:28
  • 1
    $\begingroup$ @uhoh your theory makes a lot of sense, because these numbers go up to 60, and there's 60 atoms. I'm curious what programs generate geometry files in this format, and which ones read such a format, or what such a format is called. I have a database of about 62,000 molecular geometries in standard XYZ format, and providing these in other formats (or a least tools to convert from one format to another) seems like it might be worthwhile. $\endgroup$ Apr 11, 2023 at 3:38

7 Answers 7

17
$\begingroup$

ASE

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

Here's an example script inspired by rsdel's answer which uses ASE:

import ase.io
import io      

xyz_file = """60
            Buckminsterfullerene  (C60 Bucky Ball)
            C      0.56182991    1.03720708   -3.34153745   
            C     -0.02867841   -0.20922332   -3.53733326   
            C      0.68868170   -1.41687735   -3.17419275   
            C      1.96798784   -1.33001731   -2.62971513   
            C      2.58298268   -0.03190141   -2.42580025   
            C      1.89418484    1.12766893   -2.77448204   
            C     -0.26911535   -2.36054907   -2.62920302   
            C      2.34254586   -2.18322716   -1.51766993   
            C      0.09052931   -3.17978756   -1.56143494   
            C      1.42288437   -3.08932569   -0.99437947   
            C      1.93147043    2.28439304   -1.89955106   
            C      3.33762840   -0.08283139   -1.18772874   
            C      3.37342937    1.02783673   -0.34763393   
            C      2.65606922    2.23549084   -0.71077448   
            C     -1.42982841   -0.40652370   -3.21677674   
            C     -0.22432522    2.13802273   -2.81706602   
            C     -2.18468203    0.65046195   -2.71318765   
            C     -1.56968716    1.94857800   -2.50927274   
            C      3.18903027   -1.41242386   -0.62647337   
            C      0.62215914    2.90882602   -1.92586946   
            C     -1.57842655   -1.73611615   -2.65552140   
            C     -2.12435266    2.52208091   -1.29751971   
            C      0.08957811    3.45949446   -0.76236333   
            C     -1.31157175    3.26219407   -0.44180686   
            C      0.84422386    3.40856444    0.47570821   
            C      2.10140370    2.80899380    0.50097868   
            C     -3.11943529    0.42168510   -1.62746090   
            C     -3.08214969    1.57840920   -0.75252997   
            C     -2.47596183   -1.95578408   -1.61302381   
            C     -3.26211692   -0.85496848   -1.08855240   
            C      2.47596171    1.95578403    1.61302382   
            C      3.26211681    0.85496841    1.08855240   
            C      3.11943527   -0.42168513    1.62746090   
            C      3.08214970   -1.57840928    0.75252996   
            C      2.12435273   -2.52208091    1.29751964   
            C      1.31157182   -3.26219401    0.44180679   
            C     -1.42288432    3.08932572    0.99437953   
            C     -2.34254574    2.18322716    1.51766996   
            C     -3.18903022    1.41242378    0.62647331   
            C     -3.37342947   -1.02783679    0.34763388   
            C     -3.33762837    0.08283136    1.18772868   
            C     -0.09052935    3.17978761    1.56143497   
            C     -0.84422391   -3.40856443   -0.47570819   
            C     -2.10140382   -2.80899386   -0.50097868   
            C      0.26911532    2.36054910    2.62920311   
            C      1.57842653    1.73611617    2.65552148   
            C     -0.08957818   -3.45949441    0.76236329   
            C     -2.65606930   -2.23549092    0.71077442   
            C     -2.58298263    0.03190137    2.42580019   
            C     -1.93147042   -2.28439306    1.89955101   
            C     -1.89418482   -1.12766893    2.77448197   
            C     -1.96798776    1.33001737    2.62971509   
            C     -0.62215921   -2.90882598    1.92586939   
            C     -0.68868164    1.41687744    3.17419279   
            C      0.22432525   -2.13802261    2.81706605   
            C      1.56968719   -1.94857793    2.50927279   
            C      2.18468202   -0.65046198    2.71318771   
            C     -0.56182981   -1.03720706    3.34153744   
            C      0.02867849    0.20922332    3.53733334   
            C      1.42982839    0.40652370    3.21677684  
            """

#read xyz file
atom = ase.io.read(io.StringIO(xyz_file), format = 'xyz')

# Take some example indices
i, j, k, l = 0, 1, 2, 3

# Distances
# between two indices 
atom.get_distance(i, j)
# between all atoms (distance matrix)
atom.get_all_distances()

# Angles
atom.get_angle(i, j, k)

# Dihedrals
atom.get_dihedral(i, j, k, l)
$\endgroup$
4
  • 4
    $\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$ Apr 20, 2021 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$ Apr 20, 2021 at 17:43
  • $\begingroup$ On a similar note, RDKit can also do the same thing $\endgroup$
    – user430
    Apr 22, 2021 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$ Apr 24, 2021 at 14:14
16
$\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
  • 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$ Apr 20, 2021 at 19:02
  • 1
    $\begingroup$ @JosePinto It works. I've added the result to my answer. $\endgroup$ Apr 20, 2021 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$ Apr 22, 2021 at 20:54
  • 2
    $\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$ Apr 23, 2021 at 0:45
  • 2
    $\begingroup$ The online tool you link is just running an old version of Open Babel with a web interface. $\endgroup$ Apr 24, 2021 at 4:27
16
$\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
  • 5
    $\begingroup$ And of course, if you want the z-matrix form, you just request it: obabel c60.xyz -ofh $\endgroup$ Apr 24, 2021 at 4:40
9
$\begingroup$

Maybe this comes a bit late. There are some tools which can do this. But I was also interested in doing such a calculation using python. I first started with an iterative process which was terribly slow. Calculating the distance matrix in scipy is much faster.

The distance matrix alone is not very helpful. You have to check every distance with the sum of the covalent radii. Please check my solution for the bond lengths. Finding angles and dihedral angles is more difficult, since you have to check every connection from atom A to atom B and C (and D). This information is in the distance matrix, but more difficult to obtain.

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

import pandas as pd                                       
from scipy.spatial.distance import pdist, squareform     
import numpy as np                                      
import io      

pd.set_option("display.max_rows", None, "display.max_columns", None)

cov_radii = {
      "C" : 0.78,
      "H" : 0.41,
      "O" : 0.66,
}

xyz_file = """60
            Buckminsterfullerene  (C60 Bucky Ball)
            C      0.56182991    1.03720708   -3.34153745   
            C     -0.02867841   -0.20922332   -3.53733326   
            C      0.68868170   -1.41687735   -3.17419275   
            C      1.96798784   -1.33001731   -2.62971513   
            C      2.58298268   -0.03190141   -2.42580025   
            C      1.89418484    1.12766893   -2.77448204   
            C     -0.26911535   -2.36054907   -2.62920302   
            C      2.34254586   -2.18322716   -1.51766993   
            C      0.09052931   -3.17978756   -1.56143494   
            C      1.42288437   -3.08932569   -0.99437947   
            C      1.93147043    2.28439304   -1.89955106   
            C      3.33762840   -0.08283139   -1.18772874   
            C      3.37342937    1.02783673   -0.34763393   
            C      2.65606922    2.23549084   -0.71077448   
            C     -1.42982841   -0.40652370   -3.21677674   
            C     -0.22432522    2.13802273   -2.81706602   
            C     -2.18468203    0.65046195   -2.71318765   
            C     -1.56968716    1.94857800   -2.50927274   
            C      3.18903027   -1.41242386   -0.62647337   
            C      0.62215914    2.90882602   -1.92586946   
            C     -1.57842655   -1.73611615   -2.65552140   
            C     -2.12435266    2.52208091   -1.29751971   
            C      0.08957811    3.45949446   -0.76236333   
            C     -1.31157175    3.26219407   -0.44180686   
            C      0.84422386    3.40856444    0.47570821   
            C      2.10140370    2.80899380    0.50097868   
            C     -3.11943529    0.42168510   -1.62746090   
            C     -3.08214969    1.57840920   -0.75252997   
            C     -2.47596183   -1.95578408   -1.61302381   
            C     -3.26211692   -0.85496848   -1.08855240   
            C      2.47596171    1.95578403    1.61302382   
            C      3.26211681    0.85496841    1.08855240   
            C      3.11943527   -0.42168513    1.62746090   
            C      3.08214970   -1.57840928    0.75252996   
            C      2.12435273   -2.52208091    1.29751964   
            C      1.31157182   -3.26219401    0.44180679   
            C     -1.42288432    3.08932572    0.99437953   
            C     -2.34254574    2.18322716    1.51766996   
            C     -3.18903022    1.41242378    0.62647331   
            C     -3.37342947   -1.02783679    0.34763388   
            C     -3.33762837    0.08283136    1.18772868   
            C     -0.09052935    3.17978761    1.56143497   
            C     -0.84422391   -3.40856443   -0.47570819   
            C     -2.10140382   -2.80899386   -0.50097868   
            C      0.26911532    2.36054910    2.62920311   
            C      1.57842653    1.73611617    2.65552148   
            C     -0.08957818   -3.45949441    0.76236329   
            C     -2.65606930   -2.23549092    0.71077442   
            C     -2.58298263    0.03190137    2.42580019   
            C     -1.93147042   -2.28439306    1.89955101   
            C     -1.89418482   -1.12766893    2.77448197   
            C     -1.96798776    1.33001737    2.62971509   
            C     -0.62215921   -2.90882598    1.92586939   
            C     -0.68868164    1.41687744    3.17419279   
            C      0.22432525   -2.13802261    2.81706605   
            C      1.56968719   -1.94857793    2.50927279   
            C      2.18468202   -0.65046198    2.71318771   
            C     -0.56182981   -1.03720706    3.34153744   
            C      0.02867849    0.20922332    3.53733334   
            C      1.42982839    0.40652370    3.21677684  
            """

#read xyz file
xyz_df = pd.read_csv(io.StringIO(xyz_file), 
      delim_whitespace=True, 
      skiprows=2, 
      names=["element", "x", "y", "z"])

#append a column with the covalent radii of the atoms
xyz_df['cov_rad'] = xyz_df.apply(lambda x: cov_radii[x.element], axis=1)

#calculate the distance matrix
dist_mat=pd.DataFrame(squareform(pdist(xyz_df.iloc[:,1:4],'euclid')))

#calculate the radii matrix
radii_sum_mat = pd.DataFrame([[x + y for x in xyz_df['cov_rad']] for y in xyz_df['cov_rad']])

#fill diagonals with zero
#distance of atom 0 to atom 0 is zero
np.fill_diagonal(radii_sum_mat.values, 0)

#print('\nDistance matrix:\n', dist_mat,'\n')
#print('Sum of the covalent radii:\n',radii_sum_mat,'\n')

#becomes 'True' if the calculated disctance is lower than the sum of atomic radii
bond_mat_true = dist_mat < radii_sum_mat
#print('Is this atom bonded to the other atom:\n',bond_mat_true,'\n')

#combine the distance matrix with the 'True' or 'False' bond table
bond_mat = dist_mat.combine(bond_mat_true,np.multiply)

#only one pair is needed, since bond length A-B = bond length B-A are equal
bond_mat.values[np.triu_indices_from(bond_mat, k=1)] = np.nan

#replace zeros with nan
bond_mat = bond_mat.replace(0, np.nan)

#drop all nan
bond_mat = bond_mat.unstack().dropna()

print('Atom_x is bonded to Atom_y & bond length:\n', bond_mat)
$\endgroup$
0
6
$\begingroup$

This is not a xyz file. If you have xyz file then you can use

  1. numpy.genfromtxt() module to read from the data.
  2. create an numpy array having each component as an element of the array.
  3. Using eular's method to find the distance
$\endgroup$
6
$\begingroup$
import math
import csv

filename = input('Nome?')
dic_atomo = {}
dic_atomo_neigh = {}

with open(filename,'r') as f:

    for line in f.readlines()[1:]:

        line = line.split()

        dic_atomo[line[0]] = (float(line[2]),float(line[3]),float(line[4]))

        dic_atomo_neigh[line[0]] =(line[6],line[7],line[8])

       

lista_d = []

print (dic_atomo)



for i in list(dic_atomo):

    x = dic_atomo[i][0]

    y = dic_atomo[i][1]

    z = dic_atomo[i][2]

    for j in list(dic_atomo_neigh[i]):

        if i != j:

            x1 = dic_atomo[j][0]

            y1 = dic_atomo[j][1]

            z1 = dic_atomo[j][2]

            print("{} - ({},{},{}) - {} - ({},{},{}) ".format(i,x,y,z,j,x1,y1,z1))

            d=math.sqrt((x1-x)**2+(y1-y)**2+(z1-z)**2)

            a = (i,j,d)

            lista_d.append(a)

        

print(lista_d)







filename=input('Nome final? ')

with open(filename, 'w') as csvfile:

    csvwriter = csv.writer(csvfile)

    csvwriter.writerows(lista_d)

    
$\endgroup$
6
$\begingroup$

RDKIT

Getting the distance between atoms, angle and dihedral between atoms can be achieved with a python library called rdkit.

For the sake of brevity, I will demonstrate using a simple compound: Ethanamine

Below I illustrate a minimal reproducible example

from rdkit.Chem import MolFromXYZBlock
from rdkit.Chem import rdDetermineBonds
from rdkit.Chem.rdMolTransforms import GetDihedralDeg, GetBondLength, GetAngleDeg

xyz_block = """10
Ethanamine
  C      1.2277      0.4098     -0.0270
  C      0.0679     -0.5616      0.0680
  N     -1.2273      0.1369     -0.0900
  H      1.2237      0.9656     -0.9741
  H      1.2144      1.1497      0.7853
  H      0.1336     -1.1313      1.0206
  H      0.1307     -1.3142     -0.7430
  H     -1.2998      0.8711      0.5835
  H     -1.9719     -0.5117      0.0544
  H      2.1852     -0.1236      0.0348   
"""

mol = MolFromXYZBlock(xyz_block)
rdDetermineBonds.DetermineConnectivity(mol)

# To get an idea of the atom ids, symbols and the bond type
for bond in mol.GetBonds():
  print(bond.GetBeginAtom().GetSymbol(),  bond.GetEndAtom().GetSymbol(), bond.GetBeginAtomIdx(), bond.GetEndAtomIdx(), bond.GetBondTypeAsDouble()) 

conf = mol.GetConformers()[0]

# Here I choose four of the atoms at random, for extracting various required informations
atom_1_id = 0
atom_2_id = 1
atom_3_id = 3
atom_4_id = 4

GetDihedralDeg(conf, atom_1_id, atom_2_id, atom_3_id, atom_4_id) 
# >> -32.793435222825245

GetBondLength(conf, atom_1_id, atom_2_id)
# >> 1.5158426699364285

GetAngleDeg(conf, atom_1_id, atom_2_id, atom_3_id)
# >> 27.824571764428438

Documentation for the three functions used above can be found here

$\endgroup$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .