To begin with, I wrote a script that gets Cartesian coordinates of molecule as input in the below. These are $x,y,z$ coordinates of H2O2 molecule.
1 O -1.7529 -0.5188 1.3324 2 O -0.4737 -0.1091 0.7774 3 H -2.2902 0.1678 0.8933 4 H 0.0636 -0.7957 1.2164
Then, the script constructs a Z-matrix with them, like this:
Z-mat : O O 1 1.45335189476 H 1 0.976176039452 2 96.5694760083 H 2 0.976131061897 1 96.5720363573 3 -179.995395182
Now I need to perform the reverse operation and use this Z-matrix as input and define $x,y,z$ coordinates for each atom. This is converting a Z-matrix to Cartesian coordinates.
My question is after setting first atom as 0,0,0
1 O 0 0 0
and the second one as 0,0,(distance from first) to put it on the z-axis
2 O 0 0 1.45335189476
How should I treat the 3rd and 4th atoms?
The 3rd atom must have coordinates that are something like this if read correctly:
3 H 0 distance*sin(angle) z2+distance*cos(angle)
z2 as the z-coordinate of atom 2. I am not sure if I should calculate this as
z2 + distance.cos(angle) or
z2 - distance.cos(angle) and what it depends on, if both are possible.
For the 4th atom, I use spherical coordinates
r, theta, phi = (0.976, 96.572, -179.995)
to calculate $x,y,z$ values from the formulas
x = r * sin(theta) * cos(phi) y = r * sin(theta) * sin(phi) z = r * cos(theta)
If there aren't any mistakes up to this point, how I will calculate Cartesian coordinates of the 4th atom using these $x,y,z$ values?
While searching I found TMPChem's work on GitHub and it does exactly what I want. However, in his work, there is a mathematical part that I don't understand:
# get local axis system from 3 coordinates def get_local_axes(coords1, coords2, coords3): u21 = get_u12(coords1, coords2) #calculating vector between that points 1-2 u23 = get_u12(coords2, coords3) #calculating vector between that points 2-3 if (abs(get_udp(u21, u23)) >= 1.0): print('\nError: Co-linear atoms in an internal coordinate definition') sys.exit() u23c21 = get_ucp(u23, u21) # unit cross product u21c23c21 = get_ucp(u21, u23c21) # unit cross product z = u21 y = u21c23c21 x = get_ucp(y, z) local_axes = [x, y, z] return local_axes
What is "getting local axis system from 3 coordinates"?
Here is some context of how this function is used:
bond_vector = get_bond_vector(atom.rval, atom.aval, atom.tval) disp_vector = np.array(np.dot(bond_vector, self.atoms[i].local_axes)) for p in range(3): atom.coords[p] = self.atoms[atom.rnum].coords[p] + disp_vector[p]
The bond vector definition is here:
def get_bond_vector(r, a, t): x = r * math.sin(a) * math.sin(t) y = r * math.sin(a) * math.cos(t) z = r * math.cos(a) bond_vector = [x, y, z] return bond_vector
Again, the only part I don't understand is
# get local axis system from 3 coordinates. What is that function doing?