diff --git a/README.md b/README.md index e1273a7..665a502 100644 --- a/README.md +++ b/README.md @@ -1,12 +1,29 @@ # xyz2pov -This Python script converts xyz coordinate files into povray scenes. +This implementation in Python and `numpy` converts xyz coordinate files into +povray scenes by + +``` shell +python ./xyz2povray.py input.xyz +``` + +to yield `input.pov` and `input.ini`. A subsequent run of `povray ./input.pov` +writes a single 800x600 px `input.png`. With `povray ./input.ini` a sequence +of 36 frames (640x420 px) yields a full rotation of the structure around +POVRay's x-axis, which can serve as template for further exploration, and to +write an animated .gif e.g. on [ezgif](https://ezgif.com/) + +![Example for glucose](glucose.gif) + +All atoms starting from 1 (hydrogen, H) up to 96 (curium, Cm) are supported +by a color scheme close to the one used by Jmol. If the interatomic +distance is less or equal the sum of the corresponding covalent radii +compiled by [Wikipedia](https://en.wikipedia.org/wiki/Covalent_radius), two +atoms are considered to be bound together; the model depicts them as +connected by a tube. An atom "too far out" to bind with an other atom is +going to be displayed as sphere. -So far, only hydrogen and carbon atoms are supported, but this is extendable for -other atoms. The style is similar to Avogadro's tube model. The camera is automatically directed to the molecule's center of mass and placed at an appropriate distance. Two lights are placed relative to the camera automatically. -TODO: automatic rotation of the camera to get a horizontal view of the -molecule. diff --git a/glucose.gif b/glucose.gif new file mode 100644 index 0000000..3f9e6b6 Binary files /dev/null and b/glucose.gif differ diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..a05b67c --- /dev/null +++ b/requirements.txt @@ -0,0 +1,3 @@ +# Written for Linux Debian 12/bookworm with Python 3.11.2; equally +# known to work in Linux Xubuntu 22.04.2 Jammy with Python 3.10.6. +numpy==1.24.3 diff --git a/benzene.xyz b/test_molecules/benzene.xyz similarity index 100% rename from benzene.xyz rename to test_molecules/benzene.xyz diff --git a/test_molecules/diethylether.xyz b/test_molecules/diethylether.xyz new file mode 100644 index 0000000..4d0f712 --- /dev/null +++ b/test_molecules/diethylether.xyz @@ -0,0 +1,17 @@ +15 + +C 1.05226 -0.05118 -0.07983 +C 2.56736 -0.00919 -0.07150 +O 3.06571 -1.31944 -0.33155 +C 4.49062 -1.34380 -0.33639 +C 4.95531 -2.75883 -0.61723 +H 0.63467 0.94053 0.11699 +H 0.68037 -0.74731 0.67921 +H 0.68037 -0.40467 -1.04720 +H 2.92318 0.33323 0.90628 +H 2.92318 0.68069 -0.84440 +H 4.86916 -1.01676 0.63832 +H 4.86916 -0.66931 -1.11230 +H 4.57512 -3.44928 0.14301 +H 6.04819 -2.81488 -0.62835 +H 4.57512 -3.10661 -1.58352 diff --git a/test_molecules/ethylpyridine.xyz b/test_molecules/ethylpyridine.xyz new file mode 100644 index 0000000..a955652 --- /dev/null +++ b/test_molecules/ethylpyridine.xyz @@ -0,0 +1,19 @@ +17 + +C 1.41296 0.10302 0.00141 +C 0.64113 1.26139 0.02360 +C -0.73734 1.13219 0.03626 +N -1.37948 -0.05512 0.03497 +C -0.60408 -1.16294 0.01072 +C 0.78738 -1.14561 -0.01734 +C 1.58000 -2.42169 -0.05489 +C 1.88600 -2.95600 1.33548 +H 2.49742 0.17960 -0.00851 +H 1.10384 2.24198 0.02822 +H -1.38148 2.00661 0.04696 +H -1.15587 -2.09946 0.00969 +H 2.51732 -2.25594 -0.59970 +H 1.02859 -3.17773 -0.62711 +H 2.48501 -2.24360 1.91252 +H 2.44981 -3.89166 1.26320 +H 0.96737 -3.15694 1.89715 diff --git a/test_molecules/glucose.xyz b/test_molecules/glucose.xyz new file mode 100644 index 0000000..6507784 --- /dev/null +++ b/test_molecules/glucose.xyz @@ -0,0 +1,26 @@ +24 + +O 4.59506 0.44106 -3.18607 +C 4.05938 0.72532 -1.89257 +C 2.98776 -0.30438 -1.50340 +O 1.91330 -0.23175 -2.44604 +C 0.85139 -1.11865 -2.16127 +O -0.18177 -0.90552 -3.11301 +C 1.33510 -2.58932 -2.18448 +O 0.27087 -3.51063 -1.92486 +C 2.46110 -2.74119 -1.15473 +O 3.00118 -4.06417 -1.13001 +C 3.58139 -1.71941 -1.41795 +O 4.55971 -1.75973 -0.36979 +H 5.16453 1.19425 -3.42069 +H 3.61335 1.72785 -1.93053 +H 4.88306 0.73137 -1.17545 +H 2.60456 -0.01108 -0.51419 +H 0.41460 -0.88077 -1.18320 +H 0.04931 -1.42082 -3.90159 +H 1.71521 -2.84259 -3.18064 +H 0.08274 -3.51136 -0.97050 +H 2.05939 -2.58506 -0.14388 +H 3.30065 -4.28731 -2.03033 +H 4.08740 -1.99062 -2.35084 +H 4.30606 -2.46202 0.25763 diff --git a/test_molecules/iodo.xyz b/test_molecules/iodo.xyz new file mode 100644 index 0000000..f7eb6fa --- /dev/null +++ b/test_molecules/iodo.xyz @@ -0,0 +1,26 @@ +24 + +C -3.05938 -3.21987 0.72725 +N -4.16473 -3.37773 1.50103 +O -4.09462 -2.71557 2.76828 +P -2.89848 -3.15224 3.79443 +S -1.38288 -2.34994 2.64644 +C -1.80899 -2.78556 1.00485 +C -0.73339 -2.67949 -0.05568 +C 0.42871 -3.68938 0.06469 +O 1.21780 -3.42398 1.23177 +C 0.03967 -5.17720 0.08063 +C -0.71431 -5.69635 -1.15397 +I 0.52239 -5.56703 -2.95916 +C -1.09584 -7.16870 -0.97145 +H -3.20027 -3.51007 -0.31360 +H -0.32054 -1.66251 -0.02114 +H -1.18046 -2.78232 -1.05277 +H 1.10831 -3.50409 -0.77306 +H 0.59060 -3.35347 1.97817 +H -0.56338 -5.37234 0.97735 +H 0.95274 -5.76943 0.23428 +H -1.62604 -5.11610 -1.32373 +H -0.21026 -7.80096 -0.83274 +H -1.64183 -7.54723 -1.84341 +H -1.74008 -7.29709 -0.09483 diff --git a/test_molecules/pentylbenzene.xyz b/test_molecules/pentylbenzene.xyz new file mode 100644 index 0000000..f2a3912 --- /dev/null +++ b/test_molecules/pentylbenzene.xyz @@ -0,0 +1,29 @@ +27 + +C -0.18369 0.78488 -0.95826 +C -1.06696 -0.27641 -0.75656 +C -0.90247 -1.12526 0.33602 +C 0.14900 -0.91472 1.22608 +C 1.03081 0.14805 1.02734 +C 0.87632 1.00904 -0.06840 +C 1.81056 2.17724 -0.27248 +C 3.29089 1.79735 -0.39488 +C 3.59115 0.92840 -1.61860 +C 5.08149 0.59670 -1.69771 +C 5.39771 -0.27792 -2.89994 +H -0.32884 1.43729 -1.81648 +H -1.88811 -0.43862 -1.45070 +H -1.59375 -1.94708 0.49639 +H 0.28098 -1.57690 2.07761 +H 1.83998 0.30093 1.73818 +H 1.68607 2.85763 0.57974 +H 1.51423 2.74785 -1.16117 +H 3.62859 1.28375 0.51354 +H 3.87610 2.72310 -0.46124 +H 3.27954 1.45419 -2.52869 +H 3.01784 -0.00437 -1.57039 +H 5.39494 0.07698 -0.78485 +H 5.66594 1.52203 -1.76125 +H 4.85596 -1.22749 -2.84735 +H 6.46886 -0.49940 -2.93416 +H 5.12350 0.22420 -3.83345 diff --git a/test_molecules/thiophene.xyz b/test_molecules/thiophene.xyz new file mode 100644 index 0000000..02fe7d0 --- /dev/null +++ b/test_molecules/thiophene.xyz @@ -0,0 +1,11 @@ +9 + +C 1.38296 0.71961 0.00140 +C 0.10409 1.23203 -0.00494 +S -1.08316 -0.00090 -0.00999 +C 0.11198 -1.22614 -0.00521 +C 1.38752 -0.70558 0.00147 +H 2.27214 1.33793 0.00541 +H -0.19362 2.27159 -0.00668 +H -0.17913 -2.26757 -0.00725 +H 2.28065 -1.31820 0.00587 diff --git a/xyz2povray.py b/xyz2povray.py index a313bc7..b264b9c 100755 --- a/xyz2povray.py +++ b/xyz2povray.py @@ -9,26 +9,155 @@ class Atom: """define appearance of atoms as colored spheres""" - def __init__(self, species, tag, position=[0, 0, 0]): + # define the properties for different species + # + RGB values of the Jmol scheme as referenced by Jmol, see + # https://jmol.sourceforge.net/jscolors/ (accessed [2023-04-25 Tue]) + # + average covalent radii as referenced in pm on Wikipedia, see + # https://en.wikipedia.org/wiki/Covalent_radius (accessed [2023-04-26 Wed]) + properties = { + 'H': {'mass': 1, 'rad': 0.25, 'rgb': [1.00, 1.00, 1.00], 'r': 31, 'sd': 5}, + 'He': {'mass': 4, 'rad': 0.25, 'rgb': [0.85, 1.00, 1.00], 'r': 28, 'sd': 0}, + 'Li': {'mass': 6, 'rad': 0.25, 'rgb': [0.80, 0.51, 1.00], 'r':128, 'sd': 7}, + 'Be': {'mass': 9, 'rad': 0.25, 'rgb': [0.76, 1.00, 0.00], 'r': 96, 'sd': 3}, + 'B': {'mass': 11, 'rad': 0.25, 'rgb': [1.00, 0.71, 0.71], 'r': 84, 'sd': 3}, + 'C': {'mass': 12, 'rad': 0.25, 'rgb': [0.56, 0.56, 0.56], 'r': 76, 'sd': 1}, + 'N': {'mass': 14, 'rad': 0.25, 'rgb': [0.19, 0.31, 0.97], 'r': 71, 'sd': 1}, + 'O': {'mass': 16, 'rad': 0.25, 'rgb': [1.00, 0.05, 0.05], 'r': 66, 'sd': 2}, + 'F': {'mass': 19, 'rad': 0.25, 'rgb': [0.56, 0.88, 0.31], 'r': 57, 'sd': 3}, + 'Ne': {'mass': 20, 'rad': 0.25, 'rgb': [0.70, 0.89, 0.96], 'r': 58, 'sd': 0}, + 'Na': {'mass': 22, 'rad': 0.25, 'rgb': [0.67, 0.36, 0.95], 'r':166, 'sd': 9}, + 'Mg': {'mass': 24, 'rad': 0.25, 'rgb': [0.54, 1.00, 0.00], 'r':141, 'sd': 7}, + 'Al': {'mass': 27, 'rad': 0.25, 'rgb': [0.75, 0.65, 0.65], 'r':121, 'sd': 4}, + 'Si': {'mass': 29, 'rad': 0.25, 'rgb': [0.94, 0.78, 0.63], 'r':111, 'sd': 2}, + 'P': {'mass': 31, 'rad': 0.25, 'rgb': [1.00, 0.50, 0.00], 'r':107, 'sd': 3}, + 'S': {'mass': 32, 'rad': 0.25, 'rgb': [1.00, 1.00, 0.19], 'r':105, 'sd': 3}, + 'Cl': {'mass': 35, 'rad': 0.25, 'rgb': [0.12, 0.94, 0.12], 'r':102, 'sd': 4}, + 'Ar': {'mass': 40, 'rad': 0.25, 'rgb': [0.50, 0.82, 0.89], 'r':106, 'sd': 10}, + 'K': {'mass': 39, 'rad': 0.25, 'rgb': [0.56, 0.25, 0.83], 'r':203, 'sd': 12}, + 'Ca': {'mass': 40, 'rad': 0.25, 'rgb': [0.24, 1.00, 0.00], 'r':176, 'sd': 10}, + 'Sc': {'mass': 45, 'rad': 0.25, 'rgb': [0.90, 0.90, 0.90], 'r':170, 'sd': 7}, + 'Ti': {'mass': 48, 'rad': 0.25, 'rgb': [0.75, 0.76, 0.78], 'r':160, 'sd': 8}, + 'V': {'mass': 51, 'rad': 0.25, 'rgb': [0.65, 0.65, 0.78], 'r':153, 'sd': 8}, + 'Cr': {'mass': 52, 'rad': 0.25, 'rgb': [0.54, 0.60, 0.78], 'r':139, 'sd': 5}, + 'Mn': {'mass': 55, 'rad': 0.25, 'rgb': [0.61, 0.48, 0.78], 'r':161, 'sd': 8}, + 'Fe': {'mass': 56, 'rad': 0.25, 'rgb': [0.88, 0.40, 0.20], 'r':152, 'sd': 6}, + 'Co': {'mass': 59, 'rad': 0.25, 'rgb': [0.94, 0.56, 0.62], 'r':150, 'sd': 7}, + 'Ni': {'mass': 59, 'rad': 0.25, 'rgb': [0.31, 0.82, 0.31], 'r':124, 'sd': 4}, + 'Cu': {'mass': 63, 'rad': 0.25, 'rgb': [0.78, 0.50, 0.20], 'r':132, 'sd': 4}, + 'Zn': {'mass': 65, 'rad': 0.25, 'rgb': [0.49, 0.50, 0.69], 'r':122, 'sd': 4}, + 'Ga': {'mass': 70, 'rad': 0.25, 'rgb': [0.76, 0.56, 0.56], 'r':122, 'sd': 3}, + 'Ge': {'mass': 73, 'rad': 0.25, 'rgb': [0.40, 0.56, 0.56], 'r':120, 'sd': 4}, + 'As': {'mass': 75, 'rad': 0.25, 'rgb': [0.74, 0.50, 0.89], 'r':119, 'sd': 4}, + 'Se': {'mass': 79, 'rad': 0.25, 'rgb': [1.00, 0.63, 0.00], 'r':120, 'sd': 4}, + 'Br': {'mass': 80, 'rad': 0.25, 'rgb': [0.65, 0.16, 0.16], 'r':120, 'sd': 3}, + 'Kr': {'mass': 84, 'rad': 0.25, 'rgb': [0.36, 0.72, 0.82], 'r':160, 'sd': 4}, + 'Rb': {'mass': 86, 'rad': 0.25, 'rgb': [0.43, 0.18, 0.69], 'r':220, 'sd': 9}, + 'Sr': {'mass': 88, 'rad': 0.25, 'rgb': [0.00, 1.00, 0.00], 'r':195, 'sd': 10}, + 'Y': {'mass': 89, 'rad': 0.25, 'rgb': [0.58, 1.00, 1.00], 'r':190, 'sd': 7}, + 'Zr': {'mass': 91, 'rad': 0.25, 'rgb': [0.58, 0.88, 0.88], 'r':175, 'sd': 7}, + 'Nb': {'mass': 93, 'rad': 0.25, 'rgb': [0.45, 0.76, 0.79], 'r':164, 'sd': 6}, + 'Mo': {'mass': 96, 'rad': 0.25, 'rgb': [0.32, 0.71, 0.71], 'r':154, 'sd': 5}, + 'Tc': {'mass': 97, 'rad': 0.25, 'rgb': [0.23, 0.62, 0.62], 'r':147, 'sd': 7}, + 'Ru': {'mass': 101, 'rad': 0.25, 'rgb': [0.14, 0.56, 0.56], 'r':146, 'sd': 7}, + 'Rh': {'mass': 103, 'rad': 0.25, 'rgb': [0.04, 0.49, 0.55], 'r':142, 'sd': 7}, + 'Pd': {'mass': 106, 'rad': 0.25, 'rgb': [0.00, 0.41, 0.52], 'r':139, 'sd': 6}, + 'Ag': {'mass': 108, 'rad': 0.25, 'rgb': [0.75, 0.75, 0.75], 'r':145, 'sd': 5}, + 'Cd': {'mass': 112, 'rad': 0.25, 'rgb': [1.00, 0.85, 0.56], 'r':144, 'sd': 9}, + 'In': {'mass': 115, 'rad': 0.25, 'rgb': [0.65, 0.46, 0.45], 'r':142, 'sd': 5}, + 'Sn': {'mass': 117, 'rad': 0.25, 'rgb': [0.40, 0.50, 0.50], 'r':139, 'sd': 4}, + 'Sb': {'mass': 122, 'rad': 0.25, 'rgb': [0.62, 0.39, 0.71], 'r':139, 'sd': 5}, + 'Te': {'mass': 128, 'rad': 0.25, 'rgb': [0.83, 0.48, 0.00], 'r':138, 'sd': 4}, + 'I': {'mass': 127, 'rad': 0.25, 'rgb': [0.58, 0.00, 0.58], 'r':139, 'sd': 3}, + 'Xe': {'mass': 131, 'rad': 0.25, 'rgb': [0.26, 0.62, 0.69], 'r':140, 'sd': 9}, + 'Cs': {'mass': 133, 'rad': 0.25, 'rgb': [0.34, 0.09, 0.56], 'r':244, 'sd': 11}, + 'Ba': {'mass': 137, 'rad': 0.25, 'rgb': [0.00, 0.79, 0.00], 'r':215, 'sd': 11}, + 'La': {'mass': 139, 'rad': 0.25, 'rgb': [0.44, 0.83, 1.00], 'r':207, 'sd': 8}, + 'Ce': {'mass': 140, 'rad': 0.25, 'rgb': [1.00, 1.00, 0.78], 'r':204, 'sd': 9}, + 'Pr': {'mass': 141, 'rad': 0.25, 'rgb': [0.85, 1.00, 0.78], 'r':203, 'sd': 7}, + 'Nd': {'mass': 144, 'rad': 0.25, 'rgb': [0.78, 1.00, 0.78], 'r':201, 'sd': 6}, + 'Pm': {'mass': 145, 'rad': 0.25, 'rgb': [0.64, 1.00, 0.78], 'r':199, 'sd': 0}, + 'Sm': {'mass': 150, 'rad': 0.25, 'rgb': [0.56, 1.00, 0.78], 'r':198, 'sd': 8}, + 'Eu': {'mass': 152, 'rad': 0.25, 'rgb': [0.38, 1.00, 0.78], 'r':198, 'sd': 6}, + 'Gd': {'mass': 157, 'rad': 0.25, 'rgb': [0.27, 1.00, 0.78], 'r':196, 'sd': 6}, + 'Tb': {'mass': 159, 'rad': 0.25, 'rgb': [0.19, 1.00, 0.78], 'r':194, 'sd': 5}, + 'Dy': {'mass': 163, 'rad': 0.25, 'rgb': [0.12, 1.00, 0.78], 'r':192, 'sd': 7}, + 'Ho': {'mass': 165, 'rad': 0.25, 'rgb': [0.00, 1.00, 0.61], 'r':192, 'sd': 7}, + 'Er': {'mass': 167, 'rad': 0.25, 'rgb': [0.00, 0.90, 0.46], 'r':189, 'sd': 6}, + 'Tm': {'mass': 169, 'rad': 0.25, 'rgb': [0.00, 0.83, 0.32], 'r':190, 'sd': 10}, + 'Yb': {'mass': 173, 'rad': 0.25, 'rgb': [0.00, 0.75, 0.22], 'r':187, 'sd': 8}, + 'Lu': {'mass': 175, 'rad': 0.25, 'rgb': [0.00, 0.67, 0.14], 'r':175, 'sd': 10}, + 'Hf': {'mass': 179, 'rad': 0.25, 'rgb': [0.30, 0.76, 1.00], 'r':187, 'sd': 8}, + 'Ta': {'mass': 181, 'rad': 0.25, 'rgb': [0.30, 0.65, 1.00], 'r':170, 'sd': 8}, + 'W': {'mass': 184, 'rad': 0.25, 'rgb': [0.13, 0.58, 0.84], 'r':162, 'sd': 7}, + 'Re': {'mass': 186, 'rad': 0.25, 'rgb': [0.15, 0.49, 0.67], 'r':151, 'sd': 7}, + 'Os': {'mass': 190, 'rad': 0.25, 'rgb': [0.15, 0.40, 0.59], 'r':144, 'sd': 4}, + 'Ir': {'mass': 192, 'rad': 0.25, 'rgb': [0.09, 0.33, 0.53], 'r':141, 'sd': 6}, + 'Pt': {'mass': 195, 'rad': 0.25, 'rgb': [0.82, 0.82, 0.88], 'r':136, 'sd': 5}, + 'Au': {'mass': 197, 'rad': 0.25, 'rgb': [1.00, 0.82, 0.14], 'r':136, 'sd': 6}, + 'Hg': {'mass': 201, 'rad': 0.25, 'rgb': [0.72, 0.72, 0.82], 'r':132, 'sd': 5}, + 'Tl': {'mass': 204, 'rad': 0.25, 'rgb': [0.65, 0.33, 0.30], 'r':145, 'sd': 7}, + 'Pb': {'mass': 207, 'rad': 0.25, 'rgb': [0.34, 0.35, 0.38], 'r':146, 'sd': 5}, + 'Bi': {'mass': 209, 'rad': 0.25, 'rgb': [0.62, 0.31, 0.71], 'r':148, 'sd': 4}, + 'Po': {'mass': 209, 'rad': 0.25, 'rgb': [0.67, 0.36, 0.00], 'r':140, 'sd': 4}, + 'At': {'mass': 210, 'rad': 0.25, 'rgb': [0.46, 0.31, 0.27], 'r':150, 'sd': 0}, + 'Rn': {'mass': 222, 'rad': 0.25, 'rgb': [0.26, 0.51, 0.59], 'r':150, 'sd': 0}, + 'Fr': {'mass': 223, 'rad': 0.25, 'rgb': [0.26, 0.00, 0.40], 'r':260, 'sd': 0}, + 'Ra': {'mass': 226, 'rad': 0.25, 'rgb': [0.00, 0.49, 0.00], 'r':221, 'sd': 2}, + 'Ac': {'mass': 227, 'rad': 0.25, 'rgb': [0.44, 0.67, 1.00], 'r':215, 'sd': 0}, + 'Th': {'mass': 232, 'rad': 0.25, 'rgb': [0.00, 0.73, 1.00], 'r':206, 'sd': 6}, + 'Pa': {'mass': 231, 'rad': 0.25, 'rgb': [0.00, 0.63, 1.00], 'r':200, 'sd': 0}, + 'U': {'mass': 238, 'rad': 0.25, 'rgb': [0.00, 0.56, 1.00], 'r':197, 'sd': 7}, + 'Np': {'mass': 237, 'rad': 0.25, 'rgb': [0.00, 0.50, 1.00], 'r':190, 'sd': 1}, + 'Pu': {'mass': 244, 'rad': 0.25, 'rgb': [0.00, 0.42, 1.00], 'r':187, 'sd': 1}, + 'Am': {'mass': 243, 'rad': 0.25, 'rgb': [0.33, 0.36, 0.95], 'r':180, 'sd': 6}, + 'Cm': {'mass': 247, 'rad': 0.25, 'rgb': [0.47, 0.36, 0.89], 'r':169, 'sd': 3}, +# 'Bk': {'mass': 247, 'rad': 0.25, 'rgb': [0.54, 0.31, 0.89]}, +# 'Cf': {'mass': 251, 'rad': 0.25, 'rgb': [0.63, 0.21, 0.83]}, +# 'Es': {'mass': 252, 'rad': 0.25, 'rgb': [0.70, 0.12, 0.83]}, +# 'Fm': {'mass': 257, 'rad': 0.25, 'rgb': [0.70, 0.12, 0.73]}, +# 'Md': {'mass': 258, 'rad': 0.25, 'rgb': [0.70, 0.05, 0.65]}, +# 'No': {'mass': 259, 'rad': 0.25, 'rgb': [0.74, 0.05, 0.53]}, +# 'Lr': {'mass': 260, 'rad': 0.25, 'rgb': [0.78, 0.00, 0.40]}, +# 'Rf': {'mass': 261, 'rad': 0.25, 'rgb': [0.80, 0.00, 0.35]}, +# 'Db': {'mass': 262, 'rad': 0.25, 'rgb': [0.82, 0.00, 0.31]}, +# 'Sg': {'mass': 269, 'rad': 0.25, 'rgb': [0.85, 0.00, 0.27]}, +# 'Bh': {'mass': 270, 'rad': 0.25, 'rgb': [0.88, 0.00, 0.22]}, +# 'Hs': {'mass': 269, 'rad': 0.25, 'rgb': [0.90, 0.00, 0.18]}, +# 'Mt': {'mass': 278, 'rad': 0.25, 'rgb': [0.92, 0.00, 0.15]}, + } + + def __init__(self, species, tag, position=[0, 0, 0], covalent_radius=0.0, + sd_covalent_radius=0.0): self.species = species self.tag = tag self.position = np.array(position) - - if self.species == 'C': - self.mass = 12 - self.rad = 0.25 - self.rgb = [0.4, 0.4, 0.4] - else: + self.covalent_radius = covalent_radius + self.sd_covalent_radius = sd_covalent_radius + + # set the properties based on the species + if self.species in Atom.properties: + prop = Atom.properties[self.species] + self.mass = prop['mass'] + self.rad = prop['rad'] + self.rgb = prop['rgb'] + self.covalent_radius = prop['r'] + self.sd_covalent_radius = prop['sd'] + else: # default to hydrogen self.mass = 1 self.rad = 0.25 self.rgb = [0.75, 0.75, 0.75] + self.covalent_radius = 31 + self.sd_covalent_radius = 5 + print(f"WARNING: structure contains {self.species}, an atom unknown to the program.") + def __repr__(self): return "%r %r tag: %r pos: %r, %r, %r" % ( self.species, self.mass, self.tag, self.position[0], self.position[1], self.position[2]) - def toPOV(self): + def to_pov(self): + """assemble atom coordinates and RGB for the .pov file""" return "Atom(<{:6.3f},{:6.3f},{:6.3f}>, <{:4.3f}, {:4.3f}, {:4.3f}>, {:4.2f})\n".format( self.position[0], self.position[1], self.position[2], self.rgb[0], self.rgb[1], self.rgb[2], self.rad) @@ -45,7 +174,7 @@ def __init__(self, atom_a, atom_b): self.atom_b = atom_b self.ID = str(atom_a.tag) + str(atom_b.tag) - def toPOV(self): + def to_pov(self): halfway_point = (self.atom_a.position - self.atom_b.position) / 2 + self.atom_b.position @@ -64,8 +193,9 @@ def toPOV(self): def get_structure(data): + """access atomic coordinates""" atoms = np.array([]) - with open(data) as xyz: + with open(data, mode="r", encoding="utf8") as xyz: for ii, line in enumerate(xyz): line = line.split() if len(line) == 4: #and line[0] == 'C': #no hydrogen @@ -78,7 +208,7 @@ def get_structure(data): return atoms -def get_center_of_mass(Molecule): +def get_center_of_mass(molecule): """determine the molecule's center of gravity Note: this is literally by the atoms' masses, and not only by mere dimension @@ -86,7 +216,7 @@ def get_center_of_mass(Molecule): center_of_mass = np.array([0.0, 0.0, 0.0]) total_mass = 0.0 - for atom in Molecule: + for atom in molecule: center_of_mass += atom.mass * atom.position total_mass += atom.mass center_of_mass /= total_mass @@ -94,9 +224,9 @@ def get_center_of_mass(Molecule): return np.around(center_of_mass, decimals=2) -def move2origin(Molecule, center_of_mass): +def move2origin(molecule, center_of_mass): """align molecule's centre of gravity and origin of the coordinate system""" - for atom in Molecule: + for atom in molecule: atom.translate(center_of_mass) @@ -131,10 +261,66 @@ def get_args(): parser.add_argument("source_file", metavar="FILE", help="Input .xyz file about the structure.") + parser.add_argument("-v", "--verbose", + default=False, + action='store_true', + help="Write an optional detailed report to the CLI.") return parser.parse_args() +def plausible_bond(atom1, atom2,report_level=False): + """ check if two atoms could form a bond + + By comparison of the sum of the corresponding covalent radii with + the distance derived from data in the .xyz file, this procedure + shall check if the two are close enough to form any bond (i.e., + bond order is irrelevant). For testing purpose, the output now is + a print to the screen; the intent for later is the provision of a + return value.""" + + check_value = False + # computing the sum of the radii: + observed_distance = abs(npl.norm(atom1.position - atom2.position)) + theoretical_threshold = (atom1.covalent_radius + atom2.covalent_radius) / 100 + + # add the sd into the picture + # one sigma: mean value +/- 68% of the Gaussian distribution + # two sigma: mean value +/- 95% of the Gaussian distribution + # three sigma: +/- 99.7% of the Gaussian distribution + sum_sd = (atom1.sd_covalent_radius + atom2.sd_covalent_radius) / 100 + ubound_threshold_with_sd = theoretical_threshold + (3* sum_sd) + + # a covalent bound now is set .true. below the limit of ubound + # (different to xyz2mol, bond order is not of interest here) + observed_distance = abs(npl.norm(atom1.position - atom2.position)) + if observed_distance <= ubound_threshold_with_sd: + check_value = True + + report_level = args.verbose + if report_level: + print(f"index: {atom1.tag:4} type: {atom1.species:>2}\ + radius: {atom1.covalent_radius:>3} pm sd(radius): {atom1.sd_covalent_radius:3} pm") + + print(f"index: {atom2.tag:4} type: {atom2.species:>2}\ + radius: {atom2.covalent_radius:>3} pm sd(radius): {atom2.sd_covalent_radius:3} pm") + + print(f"{'tabulated threshold (sum of both radii):':47} {theoretical_threshold:7.2f} A") + print(f"{'sum of both sd_covalent_radii:':47} {sum_sd:7.2f} A") + print(57*"-") + print(f"{'upper bound (sum radii + 3 sd of radii):':47} {ubound_threshold_with_sd:7.2f} A") + print(57*"-") + print(f"{'Interatomic distance calculated from .xyz file:':47} {observed_distance:7.2f} A") + + if check_value: + print("This distance is less than the sum of the covalent radii.") + else: + print("Warning: This distance is too large for a covalent bond.") + print("") + + return check_value + + if __name__ == '__main__': args = get_args() @@ -186,7 +372,7 @@ def get_args(): global_settings {ambient_light rgb <0.200, 0.200, 0.200> max_trace_level 15} - background {color rgb <1,1,1>} + background {color rgb <0.8,0.8,0.8>} camera { orthographic @@ -233,7 +419,7 @@ def get_args(): povfile.write(default_settings) for atom in molecule: - povfile.write(atom.toPOV()) + povfile.write(atom.to_pov()) bond_list = [ ] #keeps track of the bond halfwaypoints just to avoid double counting @@ -241,15 +427,15 @@ def get_args(): for atom2 in molecule: bond = Bond(atom1, atom2) if (atom1 != atom2 and - abs(npl.norm(atom1.position - atom2.position)) <= 1.6 + # abs(npl.norm(atom1.position - atom2.position)) <= 1.6 + plausible_bond(atom1, atom2) and bond.ID not in bond_list and bond.ID[::-1] not in bond_list): bond_list.append(bond.ID) - povfile.write(bond.toPOV()) - + povfile.write(bond.to_pov()) povfile.write('\n}') - # declare possibilty for a rotation around x in the .pov of scene + # declare possibility for a rotation around x in the .pov of scene rotation_block_a = """ union{ molecule @@ -262,7 +448,7 @@ def get_args(): # To eventually render a sequence of 36 frames, this requires a run # `povray benzene.ini` instead of `povray benzene.pov`. Ascertain the # simultaneous presence of `benzene.pov and `benzene.ini` in the very - # same writeable folder.) + # same folder with permission to write records.) rotation_block_b = "" rotation_block_b += "".join(['Input_File_Name="', output_pov, '"']) rotation_block_b +="""