An Original FORTRAN Code Combined with an XMGRACE^© Graphical Routine Aimed at the Parameter Calculations and Drawing of Newton Diagrams for Crossed Molecular Beams (CMB) Experiments Planning

This work presents a detailed computational method for analyzing and visually representing molecular collision dynamics in crossed molecular beam (CMB) experiments. An original UNIX-based FORTRAN code was created to calculate important kinematic parameters and build Newton diagrams for reactive systems where two monoenergetic beams intersect at any angle between 0° and 180°. The program finds center-of-mass velocities, collisional energies, product scattering velocities, angular ranges, and coordinate changes between laboratory and center-of-mass reference frames, using the principles of classical mechanics. A related XMGRACE-compatible routine creates high-resolution Newton diagrams that are ready for publication. This allows users to visualize reaction kinematics, scattering circles, and time-of-flight (TOF) detection angles accurately. An extra DOS/GW-Basic executable is included for quick on-screen visualization. The method has been thoroughly tested with different reactive channels and various chemical systems. By combining analytical derivations with automated graphical output, the tool provides a flexible platform for planning CMB experiments, interpreting dynamic results, and comparing theoretical predictions with observed scattering behavior. The codes are free for scientific use with proper credit.

Newton Diagrams are vector representations of the collisional events between two particles treated based on classical physics. Such types of diagrams are extensively used for the planning and the cinematic analysis of collisional events that typically occur in Crossed Molecular Beams (CMB) experiments performed in a high vacuum environment, which are at the basis of the Reaction Dynamics investigations [1-7].

The present work describes in detail a UNIX-based code aimed at the collisional parameter calculation and at the drawing of Newton Diagrams for two reacting molecular beams intersected at any angle g between 0 and 180 degrees. The global code, written in FORTRAN, was implemented by the author during his Ph.D. course in Chemistry and has been used and tested for a great variety of reactions.

It consists of two main parts: The first one calculates the reaction parameters of two monoenergetic beams of particle interacting at a given angle g, for a given elementary reaction; the second one generates and save an XMGRACE^©-readable file which after the launch/opening of an UNIX^© executable, draws an high-resolution graphical scheme of the reactions knematics, which in turn can be modified, printed or exported in the desired format for its final inclusion in a publication.

Moreover, it is also alleged to be a simple GW-Basic^©/DOS^© executable which can generate a Newton Diagram for a quick on-screen visualization for a given reaction.

The codes cited in this article can be freely used or copied by the reader(s), provided that due credit is given by the user(s) to the author.

II. Theoretical treatment of the collision of two particles, according to the classical physics

II.a. Vector derivation of the expressions for Center-Mass velocity (V_c.m.) and angle (Q_c.m.): Considering the generic reaction:

A + BC → AB + C (I)

Between two particle reagent beams A and BC crossed at whatever angle g (with 0⁰ g< 1800) e with velocities V_A e V_BC, for the conservation of the linear momentum in the C.M. reference system:

m_AU_A + m_BC<gBC = 0 (A.I.1)

And the vectorial differences:

UA = V_A - V_c.m. (A.I.2a)

UBC = V_BC - V_c.m. (A.I.2b)

(Figure II.1.) We obtain from (A.I.1):

mAV_A - mAV_c.m. + m_BCV_BC - m_BCV_c.m. = 0

Namely:

V_c.m. = (mAV_A + m_BCV_BC) / (mA + m_BC) (A.I.3)

Using the unit vectors i, j for representing V_A and V_BC (Figure A.I.1), for g = 900:

V_BC = V_BC i (A.I.4a)

V_A = V_A j (A.I.4b)

By substituting the (A.I.4a, b) into (A.I.3):

V_c.m. = [(m_BCV_BC)/(m_A+m_BC)]( i ) + [(mAV_A)/(m_A+m_BC)]( j ) (A.I.5)

By expressing the vector V_c.m. as a function of its components parallel to the X and Y axes of the Laboratory reference system:

V_c.m. = [a]( i ) + [b]( j ) (A.I.6)

The modules of these component vectors a e b correspond to the Cartesian coordinates of the point P (the tip of the V_c.m. vector) in the reference system XOY.

By confronting (A.I.5) with (A.I.6) we get:

a = (m_BCV_BC) / (m_A+m_BC) e b = (mAV_A) / (m_A+m_BC) (A.I.7)

Moreover, by letting p_BC = m_BCV_BC and pA = mAV_A,:

V_c.m. = (pA² + p_BC²)^1/2 / (mA + m_BC) (A.I.8a)

and

tan (Q_c.m.) = (a / b) = (p_BC / pA) (A.I.8b)

For 00 < g <900, we have:

V_BC = (V_BC)Xi + (V_BC)Yj (A.I.9a)

V_A = V_A j (A.I.9b)

The modules of the component vectors of V_BC are given by:

(V_BC)X = vBCsin g; (V_BC)Y = V_BCcos g (A.I.10a,b)

By substituting the (A.I.10a, b) into (A.I.9a) and then the (A.I.9a, b) into (A.I.3), we obtain:

V_c.m. = [(m_BCV_BCsing)(i)+(mAvA+m_BCV_BCcosg)(j)]/(m_A+m_BC), (A.I.11)

from which:

a = (m_BCV_BCsin g) / (m_A+m_BC), (A.I.12a)

b = (mAV_A+m_BCV_BCcos g) / (m_A+m_BC) (A.I.12b)

And:

V_c.m.=(pA²+p_BC²+2p_Ap_BCcosg)^1/2/(m_A+m_BC) (A.I.13)

For 900 < g <1800:

V_BC = (V_BC)X i + (V_BC)Y (-j) (A.I.14a)

V_A = V_Aj (A.I.14b)

The modules of the component vectors of vBC are given by:

(V_BC)X = V_BCsin g (A.I.15a)

(V_BC)Y = - V_BCcos g (A.I.15b)

By substituting the (A.I.15 a,b) into (A.I.14a) and then the (A.I.14a,b) into (A.I.3), we get the (A.I.11) again:

V_c.m.=[(m_BC·V_BCsing)(i)+(mA·vA+m_BC·V_BC·cosg)(j)]/(m_A+m_BC), (A.I.11)

the (A.I.12a,b):

a = (m_BC·V_BCsing) / (m_A+m_BC),

b = (mA·V_A+m_BC·V_BCcos g) / (m_A+m_BC)

And the:

V_c.m.= (pA² + p_BC² + 2p_A·p_BCcos g)^1/2 /(m_A+m_BC) (A.I.13)

Generally, it is more convenient to express the Q_c.m. angle as a function of V_A, V_c.m., and u_A rather than directly as the arctangent of the ratio a / b, because in this fashion it can be obtained a single general formula effective for each value of g between 00 and 1800.

Based on a generic reactive system geometry (Figure I.2.; Figure II.2.), we have:

UA² = V_A² + V_c.m.2 - 2 V_A · V_c.m.cos(Q_c.m.) (A.I.14)

hence:

Q_c.m.= arccos[(V_A² + V_c.m.2 - UA²) / (2 V_A· V_c.m.)]. (A.I.15)

But:

V_c.m.2 = a² + b2 and UA² = a² + (V_A - b)²

Then:

Q_c.m.= arccos[(b) / (a² + b2)]. (A.I.16)

By substituting into (A.I.16) the (A.I.12a, b) we obtain finally:

Q_c.m.= arccos[(pA + p_BCcos g) /(pA² + p_BC² +

+ 2p_A·p_BC·cos g)^1/2], (A.I.17)

Valid for each g between 00 and 1800.

II.b. Relations between the Cartesian coordinate system of the Center of Mass and the Laboratory system for a generic point P in the velocity plane.

The Center of Mass reference system can be seen as the result of a translation of the origin O(0,0) of the Laboratory reference system axes V_xLAB, V_yLAB into the point P(a, b), which is the Center of Mass followed by an anti-clockwise rotation of f degrees (Figure A.II.1), where f is given by :

f = arctan(V_BC sing / [V_A - V_BC cosg]) (A.II.1)

and 00 < g < 1800.

It follows that the relations which allow the passage from Center of Mass coordinates to the Laboratory ones for a generic point P are given by:

V_xLAB = V_xc.m. cos f - V_yc.m. sin f + a (A.II.2 a,b)

V_yLAB = V_xc.m. sin f + V_yc.m. cos f + b

II.b.1. Relations between Cartesian and Polar Coordinates in the two reference systems

Based on Figure A.II.1, we can also write:

V_xLAB = V_P sin Q (A.II.3)

V_yLAB = V_P cos Q

And:

V_xc.m. = U_P sin q (A.II.4)

V_yc.m. = U_P cos q

Where V_P, U_P are respectively the P product velocity in the Laboratory reference system and in the Center of Mass one and Q e q the corresponding angles.

By substituting the (A.II.3) and (A.II.4) into (A.II.2), we have:

V_P sin Q = U_P sin q cos f - U_P cos q sin f + a (A.II.5)

V_P cos Q = U_P sin q sin f + U_P cos q cos f + b

Namely:

V_P sin Q = U_P sin( q - f ) + a (A.II.6 a, b)

V_P cos Q = U_P cos(q - f ) + b

Dividing member by member the (A.II.6):

tan Q = (U_P sin( q - f ) + a) / (U_P cos( q - f ) + b) (A.II.7)

By combining the (A.II.6 a, b) as a function of the angle Q:

tan (q - f) = (V_P sin(Q) - a) / (V_P cos(Q) - b) (A.II.8)

A.II.3. Jacobian operator and relations for the conversion of the Laboratory reference frame angles Q into the Center of Mass reference system ones.

Between the product flux in the two reference systems, the following relation holds:

I_LAB(Q,V_P) = (V_P²/U_P²)· I_c.m.(q,U_P) (A.II.9)

The product's density measured in the Laboratory reference system is:

N(Q,V_P) = (V_P/U_P²)· I_c.m.(q,U_P) (A.II.10)

Where (V²/U²) is the Jacobian operator, which converts the effective Center of Mass flux into the one observed in the Laboratory reference frame.

By expressing the velocity of a generic product’s particle P as V_P in the Laboratory frame as a function of the corresponding velocity in the Center of Mass reference system U_P and of the Center of Mass velocity V_c.m. :

V_P² = U_P² + V_c.m.2 - 2U_P ·V_c.m. cos w. (A.II.11)

Where w is the {OPK} triangle angle defined by the half-lines where are the vectors U_P and V_c.m. (Figure A.II.3.). The expressions for w as a function of the angle Q formed by the vectors V_P and V_A are given by:

w = p - f - Q_c.m. + q for Q < Q_c.m.

w = p for Q = Q_c.m.

w = p + f + Q_c.m. - q for Q > Q_c.m. (A.II.12 a, b, c)

Dividing each part of (A.II.11) by U_P², it can be obtained the expression for the "Jacobian" as a function of w:

(V_P²/U_P²) = 1 + (V_c.m. / U_P)² - 2(V_c.m. / U_P) cos w. (A.II.13)

or:

(V_P/U_P²) = {[1 + (V_c.m. / U_P)² - 2(V_c.m. / U_P) cos w]½}/U_P (A.II.14)

Applying the sine’s theorem to the same triangle {OPK} above, by varying the angle, Q it is possible to obtain expressions for the conversion of the Q angles of the Laboratory system to the angles q of the center of Mass reference system (Figure A.II.4):

for Q < Q_c.m. q = Q + f - arcsin {(V_c.m. / U_AB)sin(Q_c.m. - Q)}

for Q = Q_c.m. q = Q_c.m. + f

for Q > Q_c.m. q = Q + f + arcsin {(V_c.m. / U_AB)sin(Q - Q_c.m.)}. (A.II.14 a, b, c)

A.II.4. Collisional energy (Ec), angular range amplitude (DA), velocity (vCM), and angle (QCM) of the center of mass dependence on g (intersection angle of the reagents beams).

Considering two reagent beams A and BC crossed at an angle g between 00 and 1800, the collisional energy is given by:

E_c = (½) m V_r² (A.II.15)

Where m is the reduced mass of the reagents and their relative velocity:

V_r = (V_A 2+V_BC²- 2V_A·V_BC cos g)^1/2 (A.II.16)

E_c as a function of g in the range (0, p) is in general always growing when the angle g increases, with a flection at the point (p/2, ½ m(V_A² + V_BC²) of the plane {Ec, g}, corresponding to a beam configuration of 900. The velocity vector of a generic product AB in the Center of mass reference system is given by:

U_AB = (2E_T / [{m_AB² / m^C} + m_AB])^1/2 (A.II.17)

Where the total energy E_T is given by the sum of the collisional energy and the exothermicity of the considered reaction:

E_T = (½)m V_r² - DH⁰₀ (A.II.18)

It is always growing in the interval (0,p) of the angle g (Figure A.II.3b ). If the quantity:

(2V_A·V_BC) / (V_A 2+ V_BC² - 2DH⁰) (A.II.19)

Is less then 1, the U_AB = f(g) a flex point at g = (2V_A V_BC) / (V_A 2+V_BC² - 2DH⁰₀).

The V_c dependence g is given by:

V_c.m. =(pA² + p_BC² + 2p_A·p_BCcos g)^1/2 / (m_A+m_BC) (A.II.20)

It is always decreasing when the value of g increases, with an inflexion point at g_F = arccos(-p_y / p_x), if p_y/p_x < 1, and g_F = arccos(-p_x/p_y) if p_y/p_x> 1, with 900 < g F £ 1800.

The angular range amplitude of a scattering circle with respect to the Laboratory frame is:

DA(g) = Q_max - Q_min = 2 arcsin (U_P / V_c.m..) (A.II.23)

Where:

Q_Max = Q_c.m. + arcsin (U_P /V_c.m..) (A.II.21)

QMin_Min = Q_c.m. - arcsin (U_P / V_c.m..). (A.II.22)

And Q_c.m. is the center of mass vector angle, U_P is the scattering velocity of the product P(= AB, C) with respect to the center of mass frame, and V_c.m. is the center of mass vector module.

The function DA(g) = 2 arcsin (U_P / V_c.m..) (with U_P £ V_c.m..) is always growing in the interval (0, gI), in cui gI is the limit angle (lesser than 1800) between the beams when uP = vc.m.. (isotachic point: iso - tachis =same velocity).

The general expression for the center of mass angle is given by:

Q_c.m. = arccos[(pA + p_BC·cos g) / (pA² + p_BC² + 2p_Ap_BCcos g)^1/2] (A.II.24)

By the substitution: R = (pA / p_BC ), it can be written as:

Q_c.m. = arccos[(R + cos g) / (R² + 1 + 2Rcos g)^1/2] (A.II.25)

It can be demonstrated that for reactions with 0BC), the corresponding curves Q_c.m. = f(g) are always growing when g increasing and are comprised in the cartesian plane {Q_c.m., g} region between the straight lines :

Q_c.m. = g and Q_c.m. = ½(g) (Figure A.II.5).

Q_c.m. = ½(g) is the (A.II.25) for R = 1, namely for those reactions with pA = p_BC. If R > 1 (namely for ractions with pA > p_BC), the corresponding curves Q_c.m. = f(g) show a relative maximum for gMAX = arccos(-1 / R) = arccos(- p_BC / pA), and are comprised between the two straight lines:

Q_c.m. = ½(g), Q_c.m. = 0.

For R tending to +¥, the curves Q_c.m. = f(g) approximate to the line Q_c.m. = 0, while for R tending to 0, they approach asymptotically the Q_c.m. = g.

(Figure A.II.5)

III. Algorithm for the drawing of Newton Diagrams

Considering a generic reactive channel relative to a given multichannel reaction:

A + BC ® AB(i) + C(i)

Between two beams of reagents A and BC crossed at a whatever angle g (between 00 and 1800) with velocities V_A and V_BC, respectively. If N is the number of the global reaction, it is possible to define the following algorithm for the Newton Diagrams drawing:

Input data subroutine

Read:

masses velocities of the reagents: mA, V_A, m_BC, V_BC

Angle between the two beams: g

number of reactive channels: N

TOF spectra angles number: NTOF

TOF spectra angles:

for i = 1 to NTOF { Q(i) } end

Produces masses and reactive channels’ Enthalpies:

For i = 1 to N { m_ABc⁰₀ (i) } end

Interval extremes for the drawing: V_Xmax=V_Ymax

SUBROUTINE for the cinematic parameters calculations

Calculate:

Reagents total mass: m_T = m_ABC

Reagents reduced mass: m = m_A · m_BC / m_T

Reagents linear momenta:

p_BC = m_BCV_BC;

p_A = m_AV_A

Scattering circles coordinates:

a = (p_BC sin g) / m_T ; b = ( p_A + p_BC cos g) / m_T

Center of Mass velocity modulus:

V_c.m. = ( a² + b2 )

Center of Mass velocity vector angle:

Q_c.m.= arccos [(b) / (a² + b2)].

Relative velocity modulus:V_r² = V_A 2+V_BC² -2V_AV_BCcos g

Collisional energy:

E_c = (1/2)mV_r²

Total energies for the reactive channels: For i = 1 to N

{E_T(i) = Ec - DH⁰₀ (i)}

end

Products scattering circles rays:

For = i = 1 a N

U_AB (i) ={2E_T/([ m_AB 2 (i) / mC (i)] + m_AB (i))}1/2 ;

U_C (i) = ( m_AB (i) / mC (i)) · U_AB (i)

End

Reactive circles range EXTREMES angles SCATTERING calculation SUBROUTINE

Calculate:

For i = 1 to N

if U_{A B} (i) £ V_c then

i-th Circle Angular range amplitude:

DA(i) = 2 · arcsin (U_{A B} (i) / V_c.m. )

Minimum angle: Amin(i) = Q_c.m. - DA(i) / 2

Maximum angle: Amax(i) = Q_c.m. + DA(i) / 2

if U_{A B} (i) > V_c.m. then i = i+1; continue

end

Coordinates of the scattering circles points calculations subroutine

Calculate:

For i = 1 to N

for q = 0 to 2p step (50 / 1800)·p

x(i) = a + U_{A B} (i)· cos q

y(i) = b + U_{A B} (i)· sin q

end

end

SUBROUTINE for the calculations of the limiting angles of the straight lines which define the RANGES of SCATTERING and of the TOF spectra angles

Calculate:

For i = 1 to N

for e = Amin(i), Amax(i), Q(i)

y_R = Xmax·cos |e|

x_R = y_R / {tan([p/2]-e)}

end

SUBROUTINE for the T.O.F. spectra angles calculation in the Center of Mass reference system

Calcolate:

f = arctan (V_BC sin g / [V_A - V_BC cos g])

for i = 1 to NTOF

if Q(i) < Q_c.m. then

q(i) = Q(i) + f - arcsin {(V_c.m./U_AB (j) )sin(Q_c.m. - Q(i))}

if Q(i) = Q_c.m. then

q(i) = Q_c.m. + f

if Q(i) > Q_c.m. then

q (i) = Q(i) + f + arcsin {(V_c.m./U_AB (j))sin(Q(i) - Q_c.m.)}

end

Evaluated results OUTPUT SUBROUTINE

Print:

Center of Mass angle of the C.M. velocity vector: Q_c.m.

C. M. velocity vector modulus: V_c.m.

Collisional energy: Ec

Total energies of the N reactive channels E_T(i)

Reagents' elastic scattering velocities:

Ve_A, Ve_BC

Products' maximum reactive scattering velocities:

U_{A B} (i), U_C (i)

Angular Ranges extremes of the scattering circles:

Amin(i), Amax(i)

(Figure A.III.1)

Graphical Output:

Newton Diagram drawn in the velocity plane

between: (Vxmin, Vxmax) - (Vymin, Vymax)

(Figure A.III.2)

1. Fisk GA, McDonald JD, Herschbach DR. Discuss Faraday Soc. 1967;44:228.
2. Polanyi JC. Some Concepts in Reaction Dynamics. Chemica Scripta. 1986;27:229. Available from: https://www.nobelprize.org/uploads/2018/06/polanyi-lecture.pdf
3. Casavecchia P. Chemical reaction dynamics with molecular beams. Rep Prog Phys. 2000;63:355. Available from: https://iopscience.iop.org/article/10.1088/0034-4885/63/3/203/pdf
4. Lee YT. Atomic and Molecular Beam Methods. In: Scoles G, editor. Vol. 1. New York: Oxford University Press; 1987;553–568.
5. Parrish DD, Herschbach DR. J Am Chem Soc. 1973;95:6133.
6. Sibener SJ, Buss RJ, Ng CY, Lee YT. Rev Sci Instrum. 1980;51:167.
7. Levine RD, Bernstein RB. Molecular Reaction Dynamics and Chemical Reactivity. New York: Oxford University Press; 1987.