3

I have a set of points $A,B,C,D$ in 3-D space: $$A = (x_a, y_a, z_a)$$ $$B = (x_b, y_b, z_b)$$ $$C = (x_c, y_c, z_c)$$ $$D = (x_d, y_d, z_d)$$

They belong to a 3-D figure, e.g.:

enter image description here

I'm trying to define the entire plane $ABCD$ by a set of functions, based on the points I have:

$ x(y,z)$ defines the x-coordinates on the plane, depending on y and z

$ y(x,z)$ ...

$ z(x,y)$ ...

I'm finding it a difficult problem. This, I think, is the best way to solve the problem, but maybe there is a better way based on basis vectors (which I also have).

The reason: I'm trying to develop a molecular dynamics software that can handle any periodic system, so I'm trying to find the way to universalize periodic boundary conditions. For example, in a simple cubic system, I would just do (this pseudocode):

move_the_particle();    
if (particle[x] < -x_length/2)
        particle[x] = particle[x] + x_length

...to check if my particle escaped the box on the left side, and move it over to the right side.

If this is unclear please let me know. I can try to clarify further.

khaverim
  • 643
  • Are you trying to work in any periodic boundary conditions, or just a parralelepiped (6-sided) periodic space? If it's the latter, transforming your coordinate system might be a simpler approach. – Kajelad Mar 08 '17 at 03:47
  • @Kajelad any -- including triclinic, the most complicated (I think). i.e. the basis angles could all be different. – khaverim Mar 08 '17 at 03:48
  • 1
    A plane is defined with 3 points. You have an overconstrained system here. What if the plane defined by ABC is different from the plane defined by BCD? – John Alexiou Mar 09 '17 at 03:17

3 Answers3

3

A plane is defined from three points ABC using the following algorithm. How you handle the 4th point is up to you.

  1. Plane normal direction $$\mathbf{n} = (B-A) \times (C-B)$$
  2. Scalar Component $$d=-\mathbf{n} \cdot A$$
  3. Equation of plane $$ \mathbf{n}_x x+\mathbf{n}_y y+\mathbf{n}_z z = d $$

See answer to related question for how to interpret the plane equation, in terms of the properties of the plane.

Community
  • 1
John Alexiou
  • 13,816
  • I realized that pretty early on and answered thus... But you made me realized that the 4th point may indeed be in a different plane when considered, e.g. as B,C,D. I'm not sure if my systems will ever have that problem, but do you have a suggestion? (I may have to just re-calculate another plane with the 4th point and make sure they match?) – khaverim Mar 09 '17 at 03:31
  • 1
    Yes, check for a match first. Take the coordinates of D and plug them into the equation above for the plane to see how close you are. The error in distance units would be $$h = \frac{ \mathbf{n}_x x_D +\mathbf{n}_y y_D + \mathbf{n}_z z_D -d}{| \mathbf{n} |}$$ – John Alexiou Mar 09 '17 at 03:37
1

This approach may not be general enough to account for all possible periodic boundary conditions, but it might be useful to consider nonetheless. If your periodic box is a parallelepiped (the volume enclosed by 3 sets of parallel planes), it can be defined in terms of 3 vectors. Let's call them $\vec a$, $\vec b$, and $\vec c$.

enter image description here

(image courtesy of Wikipedia)

If you define your coordinates such that $\vec a$, $\vec b$, and $\vec c$ are all drawn from the origin in the figure above, then there exists an invertible linear transform defined by a matrix $A$ that transforms your periodic box into the unit cube $x,y,z\in[0,1)$ We can find this matrix from the fact that associated inverse transform goes back to your original box, mapping the coordinate vectors $\hat x$, $\hat y$, and $\hat z$ (or $\hat i$, $\hat j$, $\hat k$ if you prefer) to $\vec a$, $\vec b$, and $\vec c$ respectively. That is:

$$A^{-1}\begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix} = \vec a,\ \ \ A^{-1}\begin{bmatrix} 0 \\ 1 \\ 0 \end{bmatrix} = \vec b,\ \ \ A^{-1}\begin{bmatrix} 0 \\ 0 \\ 1 \end{bmatrix} = \vec c$$

Since these unit vectors just take the 1st, 2nd, or 3rd column of $A^{-1}$, it follows that $A^{-1}$ is just a matrix whose columns are $\vec a$, $\vec b$, and $\vec c$. We can calculate the inverse of this matrix and obtain are original transform $A$. Left multiplying our coordinate vector (which I'll call $\vec x$) by $A$ allows us to get a new vector in our normalized coordinates $\vec x' = A\vec x$. Imposing periodic boundary conditions on $\vec x'$ amounts to taking each coordinate modulo 1, since our periodic box in $\vec x'$ is the unit cube. Depending on what you're doing, it may be useful to do your computations in $\vec x'$ coordinates or simply compute $A\vec x$ each time to check if you have moved out of the periodic box. If you have, you can take the coordinates of $\vec x'$ modulo 1 and transform back to $x$.

I can't comment on the computational efficiency of this approach vs others, but it seems like a conceptually simple way to deal with the types of boundary conditions you mention. Let me know if there's anything I can add.

Kajelad
  • 14,951
1

The problem was a bit easier than I thought at first. The convention for writing an equation of a plane in 3D is simply $Ax + By + Cz + D = 0$. The steps I followed are below:

1) A plane can be defined by 3 points, so I select 3 points from the desired plane (the left side of the box), here points $\vec{A},\vec{B},$ and $\vec{C}$. Let's say $\vec{A} = \{A_x, A_y, A_z\}$.

2) We need two vectors in the plane, from which to calculate a normal vector to the plane later. The two vectors I select are $\vec{AB} = \vec{B} - \vec{A}$ and $\vec{AC} = \vec{C} - \vec{A}$.

3) The normal vector is always perpendicular to two the vectors in the plane, thus $\vec{n} = \vec{AB} \times \vec{AC}$ (cross product). I'll write $\vec{n}$ as $\{n_x,n_y,n_z\}$.

4) So the equation of the plane can be computed as

$$ n_x(x-A_x) + n_y(y-A_y) + n_z(z-A_z) = 0 $$

or

$$ n_xx - n_xA_x + n_yy - n_yA_y + n_zz - n_zA_z = 0 $$

I recognized that $ -n_xA_x - n_yA_y - n_zA_z = -(\vec{n} \cdot \vec{A})$ (dot product).

That expression is simply a scalar value which is $D$ in the final form. Thus the equation for the plane of interest is:

$$n_xx + n_yy + n_zz -(\vec{n} \cdot \vec{A}) = 0$$

which means, for the basic form $Ax + By + Cz + D = 0$,

$ A = n_x$

$ B = n_y$

$ C = n_z$

$ D = -(\vec{n} \cdot \vec{A})$

khaverim
  • 643