How to make code that can read an .xyz file and calculate distances between atoms?

Question

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.I have a xyz coordinate in this matrix form:

 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

I have tried with pyGSM but it give a error 'NoneType' object has no attribute 'group' and I don't know how to correct this because I never worked with groups in python.

from pygsm.utilities.manage_xyz import read_xyz,xyz_to_np
from pygsm.coordinate_systems.slots import Distance,Angle,Dihedral


bond = Distance(0,1)
angle = Angle(0,1,2)
dihedral=Dihedral(0,1,2,3)

geom = read_xyz('bucky.xyz')
xyz = xyz_to_np(geom)

bond_value = bond.value(xyz)
angle_value = angle.value(xyz)
dihedral_value = dihedral.value(xyz)

print(bond_value,angle_value,dihedral_value)

    import numpy as np
try:
    from . import units 
except:
    import units

import re
#import openbabel as ob

# => XYZ File Utility <= #

def read_xyz(
    filename, 
    scale=1.):

    lines = open(filename).readlines()
    lines = lines[2:]
    geom = []
    for line in lines:
        mobj = re.match(r'^\s*(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s*$', line)
        geom.append((
            scale*float(mobj.group(1),
            scale*float(mobj.group(2)),
            scale*float(mobj.group(3)),
            scale*float(mobj.group(4)),
            )))
    return geom

Can anyone help me do this?

Answer 1

My knowledge of python isn't big, but i think you can do something like this:

filename = "filename.xyz"
xyz = open(filename, "r")
line = xyz.readline()
atom, x, y, z = line.split()
# Calculate what you need here

xyz.close()

Or with a for loop:

xyz = open(filename, "r")
for line in xyz:
    print(line)
    # Calculate what you need here
xyz.close()

Answer 2

Manual parsing of the file as @MateusJunges suggested, may be a good idea because there is no generally accepted format for these xyz-files. A common issue, however, could be that many programs expect rather something like this:

60
Buckminsterfullerene  (C60 Bucky Ball)
 C      0.56182991    1.03720708   -3.34153745     2     2     6    16
 C     -0.02867841   -0.20922332   -3.53733326     2     1     3    15
 C      0.68868170   -1.41687735   -3.17419275     2     2     4     7
 ...

So the number of atoms in the first line, followed by a comment line, followed by the atom entries (without indices in the first column).

I am not familiar with PyGSM but there are many packages out there to read xyz-files. Just as an example, take MDAnalysis.

import MDAnalysis

u = MDAnalysis.Universe("structure.xyz")
atom_positions = u.coord.positions

This will give you the atom_positions as a NumPy-array of shape (n_atoms, 3) which you can use for your calculations.

Answer 3

I believe no one has answered about the distances here yet.

There are 1000 libraries for reading in xyz files. My favorite is ase https://databases.fysik.dtu.dk/ase/index.html then build a neighbor matrix in numpy like this:

import numpy as np
from ase.io import read
atoms = read(<your file>)
natoms = len(atoms)

# First pull out the atom positions arranged like an N x 3 array
positions = atoms.positions
positions = np.array(positions, ndmin=3) # This makes it a 3d array

# Here we make the actual neighbor matrix
distances = positions.repeat(natoms, axis=0) \
           - positions.repeat(natoms, axis=1) \
                        .reshape((natoms, natoms, 3))
neigh = np.linalg.norm(distances, axis=2)

# Now if you know a charecteristic bond length, you can get an 
# average coordination number
bond_length = <some number>
ave_coordination = np.count_nonzero(neigh < bond_length) / natoms