lattice_qcd_rs/
field.rs

1//! Represent the fields on the lattice.
2
3use std::iter::{FromIterator, FusedIterator};
4use std::{
5    ops::{Add, AddAssign, Div, DivAssign, Index, IndexMut, Mul, MulAssign, Neg, Sub, SubAssign},
6    vec::Vec,
7};
8
9use na::{ComplexField, Matrix3, SVector};
10use rayon::iter::FromParallelIterator;
11use rayon::prelude::*;
12#[cfg(feature = "serde-serialize")]
13use serde::{Deserialize, Serialize};
14
15use super::{
16    lattice::{Direction, LatticeCyclic, LatticeElementToIndex, LatticeLink, LatticePoint},
17    su3,
18    su3::GENERATORS,
19    thread::{run_pool_parallel_vec_with_initializations_mutable, ThreadError},
20    utils::levi_civita,
21    CMatrix3, Complex, Real, Vector8, I,
22};
23
24/// Adjoint representation of SU(3), it is su(3) (i.e. the lie algebra).
25/// See [`su3::GENERATORS`] to view the order of generators.
26/// Note that the generators are normalize such that `Tr[T^a T^b] = \delta^{ab} / 2`
27#[derive(Debug, Copy, Clone, PartialEq)]
28#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
29pub struct Su3Adjoint {
30    data: Vector8<Real>,
31}
32
33#[allow(clippy::len_without_is_empty)]
34impl Su3Adjoint {
35    /// create a new Su3Adjoint representation where `M = M^a T^a`, where `T` are generators given in [`su3::GENERATORS`].
36    /// # Example
37    /// ```
38    /// use lattice_qcd_rs::field::Su3Adjoint;
39    /// use nalgebra::SVector;
40    ///
41    /// let su3 = Su3Adjoint::new(SVector::<f64, 8>::from_element(1_f64));
42    /// ```
43    pub const fn new(data: Vector8<Real>) -> Self {
44        Self { data }
45    }
46
47    /// create a new Su3Adjoint representation where `M = M^a T^a`, where `T` are generators given in [`su3::GENERATORS`].
48    /// # Example
49    /// ```
50    /// # use lattice_qcd_rs::field::Su3Adjoint;
51    /// let su3 = Su3Adjoint::new_from_array([0_f64; 8]);
52    /// ```
53    pub fn new_from_array(data: [Real; 8]) -> Self {
54        Su3Adjoint::new(Vector8::from(data))
55    }
56
57    /// get the data inside the Su3Adjoint.
58    pub const fn data(&self) -> &Vector8<Real> {
59        &self.data
60    }
61
62    /// Get the su3 adjoint as a [`Vector8`].
63    ///
64    /// # Example
65    /// ```
66    /// # use lattice_qcd_rs::field::Su3Adjoint;
67    /// #
68    /// # let adj = Su3Adjoint::default();
69    /// let max = adj.as_vector().max();
70    /// let norm = adj.as_ref().norm();
71    /// ```
72    pub const fn as_vector(&self) -> &Vector8<Real> {
73        self.data()
74    }
75
76    /// Get the su3 adjoint as mut ref to a [`Vector8`].
77    ///
78    /// # Example
79    /// ```
80    /// # use lattice_qcd_rs::{field::Su3Adjoint, Vector8};
81    /// #
82    /// let mut adj = Su3Adjoint::default(); // filled with 0.
83    /// adj.as_vector_mut().apply(|el| *el += 1_f64);
84    /// assert_eq!(adj.as_vector(), &Vector8::from_element(1_f64));
85    ///
86    /// adj.as_mut().set_magnitude(1_f64);
87    ///
88    /// let mut v = Vector8::from_element(1_f64);
89    /// v.set_magnitude(1_f64);
90    ///
91    /// assert_eq!(adj.as_vector(), &v);
92    /// ```
93    pub fn as_vector_mut(&mut self) -> &mut Vector8<Real> {
94        self.data_mut()
95    }
96
97    /// Returns the su(3) (Lie algebra) matrix.
98    ///
99    /// # Example
100    /// ```
101    /// # use lattice_qcd_rs::{field::Su3Adjoint};
102    /// let su3 = Su3Adjoint::new_from_array([1_f64, 0_f64, 0_f64, 0_f64, 0_f64, 0_f64, 0_f64, 0_f64]);
103    /// assert_eq!(su3.to_matrix(), *lattice_qcd_rs::su3::GENERATORS[0]);
104    /// ```
105    // TODO self non consumed ?? passé en &self ? TODO bench
106    pub fn to_matrix(self) -> Matrix3<na::Complex<Real>> {
107        self.data
108            .into_iter()
109            .enumerate()
110            .map(|(pos, el)| *GENERATORS[pos] * na::Complex::<Real>::from(el))
111            .sum()
112    }
113
114    /// Return the SU(3) matrix associated with this generator.
115    /// Note that the function consume self.
116    ///
117    /// # Example
118    /// ```
119    /// # use lattice_qcd_rs::{field::Su3Adjoint};
120    /// let su3 = Su3Adjoint::new_from_array([1_f64, 0_f64, 0_f64, 0_f64, 0_f64, 0_f64, 0_f64, 0_f64]);
121    /// assert_eq!(su3.to_su3().determinant(), nalgebra::Complex::from(1_f64));
122    /// ```
123    pub fn to_su3(self) -> Matrix3<na::Complex<Real>> {
124        // NOTE: should it consume ? the user can manually clone and there is use because
125        // where the value is not necessary anymore.
126        su3::su3_exp_i(self)
127    }
128
129    /// Return exp( T^a v^a) where v is self.
130    /// Note that the function consume self.
131    pub fn exp(self) -> Matrix3<na::Complex<Real>> {
132        su3::su3_exp_r(self)
133    }
134
135    /// Create a new random SU3 adjoint.
136    ///
137    /// # Example
138    /// ```
139    /// use lattice_qcd_rs::field::Su3Adjoint;
140    ///
141    /// let mut rng = rand::thread_rng();
142    /// let distribution = rand::distributions::Uniform::from(-1_f64..1_f64);
143    /// let su3 = Su3Adjoint::random(&mut rng, &distribution);
144    /// ```
145    pub fn random<Rng>(rng: &mut Rng, d: &impl rand_distr::Distribution<Real>) -> Self
146    where
147        Rng: rand::Rng + ?Sized,
148    {
149        Self {
150            data: Vector8::<Real>::from_fn(|_, _| d.sample(rng)),
151        }
152    }
153
154    /// Returns the trace squared `Tr(X^2)`.
155    ///
156    /// It is more accurate and faster than computing
157    /// ```
158    /// # use lattice_qcd_rs::field::Su3Adjoint;
159    /// # use lattice_qcd_rs::ComplexField;
160    /// # let m = Su3Adjoint::default();
161    /// (m.to_matrix() * m.to_matrix()).trace().real();
162    /// ```
163    #[inline]
164    pub fn trace_squared(&self) -> Real {
165        // TODO investigate.
166        self.data.iter().map(|el| el * el).sum::<Real>() / 2_f64
167    }
168
169    /// Return the t coeff `t = - 1/2 * Tr(X^2)`.
170    /// If you are looking for the trace square use [Self::trace_squared] instead.
171    ///
172    /// It is used for [`su3::su3_exp_i`].
173    ///
174    /// # Example
175    /// ```
176    /// # use lattice_qcd_rs::field::Su3Adjoint;
177    /// let su3 = Su3Adjoint::from([1_f64; 8]);
178    /// let m = su3.to_matrix();
179    /// assert_eq!(
180    ///     nalgebra::Complex::new(su3.t(), 0_f64),
181    ///     -nalgebra::Complex::from(0.5_f64) * (m * m).trace()
182    /// );
183    /// ```
184    #[inline]
185    pub fn t(&self) -> Real {
186        -0.5_f64 * self.trace_squared()
187    }
188
189    /// Return the t coeff `d = i * det(X)`.
190    /// Used for [`su3::su3_exp_i`].
191    ///
192    /// # Example
193    /// ```
194    /// # use lattice_qcd_rs::field::Su3Adjoint;
195    /// let su3 = Su3Adjoint::from([1_f64; 8]);
196    /// let m = su3.to_matrix();
197    /// assert_eq!(
198    ///     su3.d(),
199    ///     nalgebra::Complex::new(0_f64, 1_f64) * m.determinant()
200    /// );
201    /// ```
202    #[inline]
203    pub fn d(&self) -> na::Complex<Real> {
204        self.to_matrix().determinant() * I
205    }
206
207    /// Return the number of data. This number is 8
208    /// ```
209    /// # use lattice_qcd_rs::field::Su3Adjoint;
210    /// # let su3 = Su3Adjoint::new(nalgebra::SVector::<f64, 8>::zeros());
211    /// assert_eq!(su3.len(), 8);
212    /// ```
213    #[allow(clippy::unused_self)]
214    pub fn len(&self) -> usize {
215        self.data.len()
216        //8
217    }
218
219    /// Return the number of data. This number is 8
220    /// ```
221    /// # use lattice_qcd_rs::field::Su3Adjoint;
222    /// let su3 = Su3Adjoint::new(nalgebra::SVector::<f64, 8>::zeros());
223    /// assert_eq!(Su3Adjoint::len_const(), su3.len());
224    /// ```
225    pub const fn len_const() -> usize {
226        8
227    }
228
229    /// Get an iterator over the elements.
230    ///
231    /// # Example
232    /// ```
233    /// # use lattice_qcd_rs::field::Su3Adjoint;
234    /// # let adj = Su3Adjoint::default();
235    /// let sum_abs = adj.iter().map(|el| el.abs()).sum::<f64>();
236    /// ```
237    pub fn iter(&self) -> impl Iterator<Item = &Real> + ExactSizeIterator + FusedIterator {
238        self.data.iter()
239    }
240
241    /// Get an iterator over the mutable ref of elements.
242    ///
243    /// # Example
244    /// ```
245    /// # use lattice_qcd_rs::field::Su3Adjoint;
246    /// # let mut adj = Su3Adjoint::default();
247    /// adj.iter_mut().for_each(|el| *el = *el / 2_f64);
248    /// ```
249    pub fn iter_mut(
250        &mut self,
251    ) -> impl Iterator<Item = &mut Real> + ExactSizeIterator + FusedIterator {
252        self.data.iter_mut()
253    }
254
255    /// Get a mutable reference over the data.
256    pub fn data_mut(&mut self) -> &mut Vector8<Real> {
257        &mut self.data
258    }
259}
260
261impl AsRef<Vector8<f64>> for Su3Adjoint {
262    fn as_ref(&self) -> &Vector8<f64> {
263        self.as_vector()
264    }
265}
266
267impl AsMut<Vector8<f64>> for Su3Adjoint {
268    fn as_mut(&mut self) -> &mut Vector8<f64> {
269        self.as_vector_mut()
270    }
271}
272
273impl<'a> IntoIterator for &'a Su3Adjoint {
274    type IntoIter = <&'a Vector8<Real> as IntoIterator>::IntoIter;
275    type Item = &'a Real;
276
277    fn into_iter(self) -> Self::IntoIter {
278        self.data.iter()
279    }
280}
281
282impl<'a> IntoIterator for &'a mut Su3Adjoint {
283    type IntoIter = <&'a mut Vector8<Real> as IntoIterator>::IntoIter;
284    type Item = &'a mut Real;
285
286    fn into_iter(self) -> Self::IntoIter {
287        self.data.iter_mut()
288    }
289}
290
291impl AddAssign for Su3Adjoint {
292    fn add_assign(&mut self, other: Self) {
293        self.data += other.data();
294    }
295}
296
297impl Add<Su3Adjoint> for Su3Adjoint {
298    type Output = Self;
299
300    fn add(mut self, rhs: Self) -> Self::Output {
301        self += rhs;
302        self
303    }
304}
305
306impl Add<&Su3Adjoint> for Su3Adjoint {
307    type Output = Self;
308
309    fn add(self, rhs: &Self) -> Self::Output {
310        self + *rhs
311    }
312}
313
314impl Add<Su3Adjoint> for &Su3Adjoint {
315    type Output = Su3Adjoint;
316
317    fn add(self, rhs: Su3Adjoint) -> Self::Output {
318        rhs + self
319    }
320}
321
322impl Add<&Su3Adjoint> for &Su3Adjoint {
323    type Output = Su3Adjoint;
324
325    fn add(self, rhs: &Su3Adjoint) -> Self::Output {
326        self + *rhs
327    }
328}
329
330impl MulAssign<f64> for Su3Adjoint {
331    fn mul_assign(&mut self, rhs: f64) {
332        self.data *= rhs;
333    }
334}
335
336impl Mul<Real> for Su3Adjoint {
337    type Output = Self;
338
339    fn mul(mut self, rhs: Real) -> Self::Output {
340        self *= rhs;
341        self
342    }
343}
344
345impl Mul<&Real> for Su3Adjoint {
346    type Output = Self;
347
348    fn mul(self, rhs: &Real) -> Self::Output {
349        self * (*rhs)
350    }
351}
352
353impl Mul<Real> for &Su3Adjoint {
354    type Output = Su3Adjoint;
355
356    fn mul(self, rhs: Real) -> Self::Output {
357        *self * rhs
358    }
359}
360
361impl Mul<&Real> for &Su3Adjoint {
362    type Output = Su3Adjoint;
363
364    fn mul(self, rhs: &Real) -> Self::Output {
365        *self * rhs
366    }
367}
368
369impl Mul<Su3Adjoint> for Real {
370    type Output = Su3Adjoint;
371
372    fn mul(self, rhs: Su3Adjoint) -> Self::Output {
373        rhs * self
374    }
375}
376
377impl Mul<&Su3Adjoint> for Real {
378    type Output = Su3Adjoint;
379
380    fn mul(self, rhs: &Su3Adjoint) -> Self::Output {
381        rhs * self
382    }
383}
384
385impl Mul<Su3Adjoint> for &Real {
386    type Output = Su3Adjoint;
387
388    fn mul(self, rhs: Su3Adjoint) -> Self::Output {
389        rhs * self
390    }
391}
392
393impl Mul<&Su3Adjoint> for &Real {
394    type Output = Su3Adjoint;
395
396    fn mul(self, rhs: &Su3Adjoint) -> Self::Output {
397        rhs * self
398    }
399}
400
401impl DivAssign<f64> for Su3Adjoint {
402    fn div_assign(&mut self, rhs: f64) {
403        self.data /= rhs;
404    }
405}
406
407impl DivAssign<&f64> for Su3Adjoint {
408    fn div_assign(&mut self, rhs: &f64) {
409        self.data /= *rhs;
410    }
411}
412
413impl Div<Real> for Su3Adjoint {
414    type Output = Self;
415
416    fn div(mut self, rhs: Real) -> Self::Output {
417        self /= rhs;
418        self
419    }
420}
421
422impl Div<&Real> for Su3Adjoint {
423    type Output = Self;
424
425    fn div(self, rhs: &Real) -> Self::Output {
426        self / (*rhs)
427    }
428}
429
430impl Div<Real> for &Su3Adjoint {
431    type Output = Su3Adjoint;
432
433    fn div(self, rhs: Real) -> Self::Output {
434        *self / rhs
435    }
436}
437
438impl Div<&Real> for &Su3Adjoint {
439    type Output = Su3Adjoint;
440
441    fn div(self, rhs: &Real) -> Self::Output {
442        *self / rhs
443    }
444}
445
446impl SubAssign for Su3Adjoint {
447    fn sub_assign(&mut self, other: Self) {
448        self.data -= other.data();
449    }
450}
451
452impl Sub<Su3Adjoint> for Su3Adjoint {
453    type Output = Self;
454
455    fn sub(mut self, rhs: Self) -> Self::Output {
456        self -= rhs;
457        self
458    }
459}
460
461impl Sub<&Su3Adjoint> for Su3Adjoint {
462    type Output = Self;
463
464    fn sub(self, rhs: &Self) -> Self::Output {
465        self - *rhs
466    }
467}
468
469impl Sub<Su3Adjoint> for &Su3Adjoint {
470    type Output = Su3Adjoint;
471
472    fn sub(self, rhs: Su3Adjoint) -> Self::Output {
473        rhs - self
474    }
475}
476
477impl Sub<&Su3Adjoint> for &Su3Adjoint {
478    type Output = Su3Adjoint;
479
480    fn sub(self, rhs: &Su3Adjoint) -> Self::Output {
481        *self - rhs
482    }
483}
484
485impl Neg for Su3Adjoint {
486    type Output = Self;
487
488    fn neg(self) -> Self::Output {
489        Su3Adjoint::new(-self.data)
490    }
491}
492
493impl Neg for &Su3Adjoint {
494    type Output = Su3Adjoint;
495
496    fn neg(self) -> Self::Output {
497        Su3Adjoint::new(-self.data())
498    }
499}
500
501/// Return the representation for the zero matrix.
502impl Default for Su3Adjoint {
503    /// Return the representation for the zero matrix.
504    ///
505    /// # Example
506    /// ```
507    /// # use lattice_qcd_rs::field::Su3Adjoint;
508    /// assert_eq!(
509    ///     Su3Adjoint::default(),
510    ///     Su3Adjoint::new_from_array([0_f64; 8])
511    /// );
512    /// ```
513    fn default() -> Self {
514        Su3Adjoint::new(Vector8::from_element(0_f64))
515    }
516}
517
518impl std::fmt::Display for Su3Adjoint {
519    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
520        write!(f, "{}", self.to_matrix())
521    }
522}
523
524impl Index<usize> for Su3Adjoint {
525    type Output = Real;
526
527    /// Get the element at position `pos`.
528    ///
529    /// # Panic
530    /// Panics if the position is out of bound (greater or equal to 8).
531    /// ```should_panic
532    /// # use lattice_qcd_rs::field::Su3Adjoint;
533    /// let su3 = Su3Adjoint::new_from_array([0_f64; 8]);
534    /// let _ = su3[8];
535    /// ```
536    fn index(&self, pos: usize) -> &Self::Output {
537        &self.data[pos]
538    }
539}
540
541impl IndexMut<usize> for Su3Adjoint {
542    /// Get the element at position `pos`.
543    ///
544    /// # Panic
545    /// Panics if the position is out of bound (greater or equal to 8).
546    /// ```should_panic
547    /// # use lattice_qcd_rs::field::Su3Adjoint;
548    /// let mut su3 = Su3Adjoint::new_from_array([0_f64; 8]);
549    /// su3[8] += 1_f64;
550    /// ```
551    fn index_mut(&mut self, pos: usize) -> &mut Self::Output {
552        &mut self.data[pos]
553    }
554}
555
556impl From<Vector8<Real>> for Su3Adjoint {
557    fn from(v: Vector8<Real>) -> Self {
558        Su3Adjoint::new(v)
559    }
560}
561
562impl From<Su3Adjoint> for Vector8<Real> {
563    fn from(v: Su3Adjoint) -> Self {
564        v.data
565    }
566}
567
568impl From<&Su3Adjoint> for Vector8<Real> {
569    fn from(v: &Su3Adjoint) -> Self {
570        v.data
571    }
572}
573
574impl From<[Real; 8]> for Su3Adjoint {
575    fn from(v: [Real; 8]) -> Self {
576        Su3Adjoint::new_from_array(v)
577    }
578}
579
580/// Represents the link matrices
581// TODO more doc
582#[derive(Debug, PartialEq, Clone)]
583#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
584pub struct LinkMatrix {
585    data: Vec<Matrix3<na::Complex<Real>>>,
586}
587
588impl LinkMatrix {
589    /// Create a new link matrix field from preexisting data.
590    pub const fn new(data: Vec<Matrix3<na::Complex<Real>>>) -> Self {
591        Self { data }
592    }
593
594    /// Get the raw data.
595    pub const fn data(&self) -> &Vec<Matrix3<na::Complex<Real>>> {
596        &self.data
597    }
598
599    /// Get a mutable reference to the data
600    pub fn data_mut(&mut self) -> &mut Vec<Matrix3<na::Complex<Real>>> {
601        &mut self.data
602    }
603
604    /// Get the link_matrix as a [`Vec`].
605    pub const fn as_vec(&self) -> &Vec<Matrix3<na::Complex<Real>>> {
606        self.data()
607    }
608
609    /// Get the link_matrix as a mutable [`Vec`].
610    pub fn as_vec_mut(&mut self) -> &mut Vec<Matrix3<na::Complex<Real>>> {
611        self.data_mut()
612    }
613
614    /// Get the link_matrix as a slice.
615    pub fn as_slice(&self) -> &[Matrix3<na::Complex<Real>>] {
616        self.data()
617    }
618
619    /// Get the link_matrix as a mut ref to a slice
620    pub fn as_slice_mut(&mut self) -> &mut [Matrix3<na::Complex<Real>>] {
621        &mut self.data
622    }
623
624    /// Single threaded generation with a given random number generator.
625    /// useful to produce a deterministic set of data but slower than
626    /// [`LinkMatrix::new_random_threaded`].
627    ///
628    /// # Example
629    /// ```
630    /// use lattice_qcd_rs::{field::LinkMatrix, lattice::LatticeCyclic};
631    /// use rand::{rngs::StdRng, SeedableRng};
632    /// # use std::error::Error;
633    ///
634    /// # fn main() -> Result<(), Box<dyn Error>> {
635    /// let mut rng_1 = StdRng::seed_from_u64(0);
636    /// let mut rng_2 = StdRng::seed_from_u64(0);
637    /// // They have the same seed and should generate the same numbers
638    /// let lattice = LatticeCyclic::<4>::new(1_f64, 4)?;
639    /// assert_eq!(
640    ///     LinkMatrix::new_determinist(&lattice, &mut rng_1),
641    ///     LinkMatrix::new_determinist(&lattice, &mut rng_2)
642    /// );
643    /// # Ok(())
644    /// # }
645    /// ```
646    pub fn new_determinist<Rng: rand::Rng + ?Sized, const D: usize>(
647        l: &LatticeCyclic<D>,
648        rng: &mut Rng,
649    ) -> Self {
650        // l.get_links_space().map(|_| Su3Adjoint::random(rng, d).to_su3()).collect()
651        // using a for loop improves performance. ( probably because the vector is pre allocated).
652        let mut data = Vec::with_capacity(l.number_of_canonical_links_space());
653        for _ in l.get_links() {
654            // the iterator *should* be in order
655            let matrix = su3::random_su3(rng);
656            data.push(matrix);
657        }
658        Self { data }
659    }
660
661    /// Multi threaded generation of random data. Due to the non deterministic way threads
662    /// operate a set cannot be reduced easily. If you want deterministic
663    /// generation you can use [`LinkMatrix::new_random_threaded`].
664    ///
665    /// # Example
666    /// ```
667    /// use lattice_qcd_rs::{field::LinkMatrix, lattice::LatticeCyclic};
668    /// # use std::error::Error;
669    ///
670    /// # fn main() -> Result<(), Box<dyn Error>> {
671    /// let lattice = LatticeCyclic::<3>::new(1_f64, 4)?;
672    /// let links = LinkMatrix::new_random_threaded(&lattice, 4)?;
673    /// assert!(!links.is_empty());
674    /// # Ok(())
675    /// # }
676    /// ```
677    ///
678    /// # Errors
679    /// Returns [`ThreadError::ThreadNumberIncorrect`] if `number_of_thread` is 0.
680    pub fn new_random_threaded<const D: usize>(
681        l: &LatticeCyclic<D>,
682        number_of_thread: usize,
683    ) -> Result<Self, ThreadError> {
684        if number_of_thread == 0 {
685            return Err(ThreadError::ThreadNumberIncorrect);
686        }
687        else if number_of_thread == 1 {
688            let mut rng = rand::thread_rng();
689            return Ok(LinkMatrix::new_determinist(l, &mut rng));
690        }
691        let data = run_pool_parallel_vec_with_initializations_mutable(
692            l.get_links(),
693            &(),
694            &|rng, _, _| su3::random_su3(rng),
695            rand::thread_rng,
696            number_of_thread,
697            l.number_of_canonical_links_space(),
698            l,
699            &CMatrix3::zeros(),
700        )?;
701        Ok(Self { data })
702    }
703
704    /// Create a cold configuration ( where the link matrices is set to the identity).
705    ///
706    /// # Example
707    /// ```
708    /// use lattice_qcd_rs::{field::LinkMatrix, lattice::LatticeCyclic};
709    /// # use std::error::Error;
710    ///
711    /// # fn main() -> Result<(), Box<dyn Error>> {
712    /// let lattice = LatticeCyclic::<3>::new(1_f64, 4)?;
713    /// let links = LinkMatrix::new_cold(&lattice);
714    /// assert!(!links.is_empty());
715    /// # Ok(())
716    /// # }
717    /// ```
718    pub fn new_cold<const D: usize>(l: &LatticeCyclic<D>) -> Self {
719        Self {
720            data: vec![CMatrix3::identity(); l.number_of_canonical_links_space()],
721        }
722    }
723
724    /// get the link matrix associated to given link using the notation
725    /// $`U_{-i}(x) = U^\dagger_{i}(x-i)`$
726    pub fn matrix<const D: usize>(
727        &self,
728        link: &LatticeLink<D>,
729        l: &LatticeCyclic<D>,
730    ) -> Option<Matrix3<na::Complex<Real>>> {
731        let link_c = l.into_canonical(*link);
732        let matrix = self.data.get(link_c.to_index(l))?;
733        if link.is_dir_negative() {
734            // that means the the link was in the negative direction
735            Some(matrix.adjoint())
736        }
737        else {
738            Some(*matrix)
739        }
740    }
741
742    /// Get $`S_ij(x) = U_j(x) U_i(x+j) U^\dagger_j(x+i)`$.
743    pub fn sij<const D: usize>(
744        &self,
745        point: &LatticePoint<D>,
746        dir_i: &Direction<D>,
747        dir_j: &Direction<D>,
748        lattice: &LatticeCyclic<D>,
749    ) -> Option<Matrix3<na::Complex<Real>>> {
750        let u_j = self.matrix(&LatticeLink::new(*point, *dir_j), lattice)?;
751        let point_pj = lattice.add_point_direction(*point, dir_j);
752        let u_i_p_j = self.matrix(&LatticeLink::new(point_pj, *dir_i), lattice)?;
753        let point_pi = lattice.add_point_direction(*point, dir_i);
754        let u_j_pi_d = self
755            .matrix(&LatticeLink::new(point_pi, *dir_j), lattice)?
756            .adjoint();
757        Some(u_j * u_i_p_j * u_j_pi_d)
758    }
759
760    /// Get the plaquette $`P_{ij}(x) = U_i(x) S^\dagger_ij(x)`$.
761    pub fn pij<const D: usize>(
762        &self,
763        point: &LatticePoint<D>,
764        dir_i: &Direction<D>,
765        dir_j: &Direction<D>,
766        lattice: &LatticeCyclic<D>,
767    ) -> Option<Matrix3<na::Complex<Real>>> {
768        let s_ij = self.sij(point, dir_i, dir_j, lattice)?;
769        let u_i = self.matrix(&LatticeLink::new(*point, *dir_i), lattice)?;
770        Some(u_i * s_ij.adjoint())
771    }
772
773    /// Take the average of the trace of all plaquettes
774    #[allow(clippy::as_conversions)] // no try into for f64
775    pub fn average_trace_plaquette<const D: usize>(
776        &self,
777        lattice: &LatticeCyclic<D>,
778    ) -> Option<na::Complex<Real>> {
779        if lattice.number_of_canonical_links_space() != self.len() {
780            return None;
781        }
782        // the order does not matter as we sum
783        let sum = lattice
784            .get_points()
785            .par_bridge()
786            .map(|point| {
787                Direction::positive_directions()
788                    .iter()
789                    .map(|dir_i| {
790                        Direction::positive_directions()
791                            .iter()
792                            .filter(|dir_j| dir_i.index() < dir_j.index())
793                            .map(|dir_j| {
794                                self.pij(&point, dir_i, dir_j, lattice).map(|el| el.trace())
795                            })
796                            .sum::<Option<na::Complex<Real>>>()
797                    })
798                    .sum::<Option<na::Complex<Real>>>()
799            })
800            .sum::<Option<na::Complex<Real>>>()?;
801        let number_of_directions = (D * (D - 1)) / 2;
802        let number_of_plaquette = (lattice.number_of_points() * number_of_directions) as f64;
803        Some(sum / number_of_plaquette)
804    }
805
806    /// Get the clover, used for F_mu_nu tensor
807    pub fn clover<const D: usize>(
808        &self,
809        point: &LatticePoint<D>,
810        dir_i: &Direction<D>,
811        dir_j: &Direction<D>,
812        lattice: &LatticeCyclic<D>,
813    ) -> Option<CMatrix3> {
814        Some(
815            self.pij(point, dir_i, dir_j, lattice)?
816                + self.pij(point, dir_j, &-dir_i, lattice)?
817                + self.pij(point, &-dir_i, &-dir_j, lattice)?
818                + self.pij(point, &-dir_j, dir_i, lattice)?,
819        )
820    }
821
822    /// Get the `F^{ij}` tensor using the clover appropriation. The direction should be set to positive.
823    /// See <https://arxiv.org/abs/1512.02374>.
824    // TODO negative directions
825    pub fn f_mu_nu<const D: usize>(
826        &self,
827        point: &LatticePoint<D>,
828        dir_i: &Direction<D>,
829        dir_j: &Direction<D>,
830        lattice: &LatticeCyclic<D>,
831    ) -> Option<CMatrix3> {
832        let m = self.clover(point, dir_i, dir_j, lattice)?
833            - self.clover(point, dir_j, dir_i, lattice)?;
834        Some(m / Complex::from(8_f64 * lattice.size() * lattice.size()))
835    }
836
837    /// Get the chromomagnetic field at a given point.
838    pub fn magnetic_field_vec<const D: usize>(
839        &self,
840        point: &LatticePoint<D>,
841        lattice: &LatticeCyclic<D>,
842    ) -> Option<SVector<CMatrix3, D>> {
843        let mut vec = SVector::<CMatrix3, D>::zeros();
844        for dir in &Direction::<D>::positive_directions() {
845            vec[dir.index()] = self.magnetic_field(point, dir, lattice)?;
846        }
847        Some(vec)
848    }
849
850    /// Get the chromomagnetic field at a given point alongside a given direction.
851    pub fn magnetic_field<const D: usize>(
852        &self,
853        point: &LatticePoint<D>,
854        dir: &Direction<D>,
855        lattice: &LatticeCyclic<D>,
856    ) -> Option<CMatrix3> {
857        let sum = Direction::<D>::positive_directions()
858            .iter()
859            .map(|dir_i| {
860                Direction::<D>::positive_directions()
861                    .iter()
862                    .map(|dir_j| {
863                        let f_mn = self.f_mu_nu(point, dir_i, dir_j, lattice)?;
864                        let lc = Complex::from(
865                            levi_civita(&[dir.index(), dir_i.index(), dir_j.index()]).to_f64(),
866                        );
867                        Some(f_mn * lc)
868                    })
869                    .sum::<Option<CMatrix3>>()
870            })
871            .sum::<Option<CMatrix3>>()?;
872        Some(sum / Complex::new(0_f64, 2_f64))
873    }
874
875    /// Get the chromomagnetic field at a given point alongside a given direction given by lattice link.
876    pub fn magnetic_field_link<const D: usize>(
877        &self,
878        link: &LatticeLink<D>,
879        lattice: &LatticeCyclic<D>,
880    ) -> Option<Matrix3<na::Complex<Real>>> {
881        self.magnetic_field(link.pos(), link.dir(), lattice)
882    }
883
884    /// Return the number of elements.
885    pub fn len(&self) -> usize {
886        self.data.len()
887    }
888
889    /// Returns wether the there is no data.
890    pub fn is_empty(&self) -> bool {
891        self.data.is_empty()
892    }
893
894    /// Correct the numerical drift, reprojecting all the matrices to SU(3).
895    ///
896    /// You can look at the example of [`super::simulation::LatticeStateDefault::normalize_link_matrices`]
897    pub fn normalize(&mut self) {
898        self.data.par_iter_mut().for_each(|el| {
899            su3::orthonormalize_matrix_mut(el);
900        });
901    }
902
903    /// Iter on the data.
904    pub fn iter(&self) -> impl Iterator<Item = &CMatrix3> + ExactSizeIterator + FusedIterator {
905        self.data.iter()
906    }
907
908    /// Iter mutably on the data.
909    pub fn iter_mut(
910        &mut self,
911    ) -> impl Iterator<Item = &mut CMatrix3> + ExactSizeIterator + FusedIterator {
912        self.data.iter_mut()
913    }
914}
915
916impl AsRef<Vec<CMatrix3>> for LinkMatrix {
917    fn as_ref(&self) -> &Vec<CMatrix3> {
918        self.as_vec()
919    }
920}
921
922impl AsMut<Vec<CMatrix3>> for LinkMatrix {
923    fn as_mut(&mut self) -> &mut Vec<CMatrix3> {
924        self.as_vec_mut()
925    }
926}
927
928impl AsRef<[CMatrix3]> for LinkMatrix {
929    fn as_ref(&self) -> &[CMatrix3] {
930        self.as_slice()
931    }
932}
933
934impl AsMut<[CMatrix3]> for LinkMatrix {
935    fn as_mut(&mut self) -> &mut [CMatrix3] {
936        self.as_slice_mut()
937    }
938}
939
940impl<'a> IntoIterator for &'a LinkMatrix {
941    type IntoIter = <&'a Vec<CMatrix3> as IntoIterator>::IntoIter;
942    type Item = &'a CMatrix3;
943
944    fn into_iter(self) -> Self::IntoIter {
945        self.data.iter()
946    }
947}
948
949impl<'a> IntoIterator for &'a mut LinkMatrix {
950    type IntoIter = <&'a mut Vec<CMatrix3> as IntoIterator>::IntoIter;
951    type Item = &'a mut CMatrix3;
952
953    fn into_iter(self) -> Self::IntoIter {
954        self.data.iter_mut()
955    }
956}
957
958impl Index<usize> for LinkMatrix {
959    type Output = CMatrix3;
960
961    fn index(&self, pos: usize) -> &Self::Output {
962        &self.data[pos]
963    }
964}
965
966impl IndexMut<usize> for LinkMatrix {
967    fn index_mut(&mut self, pos: usize) -> &mut Self::Output {
968        &mut self.data[pos]
969    }
970}
971
972impl<A> FromIterator<A> for LinkMatrix
973where
974    Vec<CMatrix3>: FromIterator<A>,
975{
976    fn from_iter<T>(iter: T) -> Self
977    where
978        T: IntoIterator<Item = A>,
979    {
980        Self::new(Vec::from_iter(iter))
981    }
982}
983
984impl<A> FromParallelIterator<A> for LinkMatrix
985where
986    Vec<CMatrix3>: FromParallelIterator<A>,
987    A: Send,
988{
989    fn from_par_iter<I>(par_iter: I) -> Self
990    where
991        I: IntoParallelIterator<Item = A>,
992    {
993        Self::new(Vec::from_par_iter(par_iter))
994    }
995}
996
997impl<T> ParallelExtend<T> for LinkMatrix
998where
999    Vec<CMatrix3>: ParallelExtend<T>,
1000    T: Send,
1001{
1002    fn par_extend<I>(&mut self, par_iter: I)
1003    where
1004        I: IntoParallelIterator<Item = T>,
1005    {
1006        self.data.par_extend(par_iter);
1007    }
1008}
1009
1010impl<A> Extend<A> for LinkMatrix
1011where
1012    Vec<CMatrix3>: Extend<A>,
1013{
1014    fn extend<T>(&mut self, iter: T)
1015    where
1016        T: IntoIterator<Item = A>,
1017    {
1018        self.data.extend(iter);
1019    }
1020}
1021
1022/// Represent an electric field.
1023#[derive(Debug, PartialEq, Clone)]
1024#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
1025pub struct EField<const D: usize> {
1026    data: Vec<SVector<Su3Adjoint, D>>, // use a Vec<[Su3Adjoint; D]> instead ?
1027}
1028
1029impl<const D: usize> EField<D> {
1030    /// Create a new "Electrical" field.
1031    pub fn new(data: Vec<SVector<Su3Adjoint, D>>) -> Self {
1032        Self { data }
1033    }
1034
1035    /// Get the raw data.
1036    pub const fn data(&self) -> &Vec<SVector<Su3Adjoint, D>> {
1037        &self.data
1038    }
1039
1040    /// Get a mut ref to the data data.
1041    pub fn data_mut(&mut self) -> &mut Vec<SVector<Su3Adjoint, D>> {
1042        &mut self.data
1043    }
1044
1045    /// Get the e_field as a Vec of Vector of [`Su3Adjoint`]
1046    pub const fn as_vec(&self) -> &Vec<SVector<Su3Adjoint, D>> {
1047        self.data()
1048    }
1049
1050    /// Get the e_field as a slice of Vector of [`Su3Adjoint`]
1051    pub fn as_slice(&self) -> &[SVector<Su3Adjoint, D>] {
1052        &self.data
1053    }
1054
1055    /// Get the e_field as mut ref to slice of Vector of [`Su3Adjoint`]
1056    pub fn as_slice_mut(&mut self) -> &mut [SVector<Su3Adjoint, D>] {
1057        &mut self.data
1058    }
1059
1060    /// Return the number of elements.
1061    pub fn len(&self) -> usize {
1062        self.data.len()
1063    }
1064
1065    /// Single threaded generation with a given random number generator.
1066    /// useful to reproduce a set of data.
1067    ///
1068    /// # Example
1069    /// ```
1070    /// # use lattice_qcd_rs::{field::EField, lattice::LatticeCyclic};
1071    /// # fn main () -> Result<(), Box<dyn std::error::Error>> {
1072    /// use rand::{rngs::StdRng, SeedableRng};
1073    ///
1074    /// let mut rng_1 = StdRng::seed_from_u64(0);
1075    /// let mut rng_2 = StdRng::seed_from_u64(0);
1076    /// // They have the same seed and should generate the same numbers
1077    /// let distribution = rand::distributions::Uniform::from(-1_f64..1_f64);
1078    /// let lattice = LatticeCyclic::<4>::new(1_f64, 4)?;
1079    /// assert_eq!(
1080    ///     EField::new_determinist(&lattice, &mut rng_1, &distribution),
1081    ///     EField::new_determinist(&lattice, &mut rng_2, &distribution)
1082    /// );
1083    /// # Ok(())
1084    /// # }
1085    /// ```
1086    pub fn new_determinist<Rng: rand::Rng + ?Sized>(
1087        l: &LatticeCyclic<D>,
1088        rng: &mut Rng,
1089        d: &impl rand_distr::Distribution<Real>,
1090    ) -> Self {
1091        let mut data = Vec::with_capacity(l.number_of_points());
1092        for _ in l.get_points() {
1093            // iterator *should* be ordered
1094            data.push(SVector::<Su3Adjoint, D>::from_fn(|_, _| {
1095                Su3Adjoint::random(rng, d)
1096            }));
1097        }
1098        Self { data }
1099    }
1100
1101    /// Single thread generation by seeding a new rng number.
1102    /// To create a seedable and reproducible set use [`EField::new_determinist`].
1103    ///
1104    /// # Example
1105    /// ```
1106    /// use lattice_qcd_rs::{field::EField, lattice::LatticeCyclic};
1107    /// # use std::error::Error;
1108    ///
1109    /// # fn main() -> Result<(), Box<dyn Error>> {
1110    /// let distribution = rand::distributions::Uniform::from(-1_f64..1_f64);
1111    /// let lattice = LatticeCyclic::<3>::new(1_f64, 4)?;
1112    /// let e_field = EField::new_random(&lattice, &distribution);
1113    /// assert!(!e_field.is_empty());
1114    /// # Ok(())
1115    /// # }
1116    /// ```
1117    pub fn new_random(l: &LatticeCyclic<D>, d: &impl rand_distr::Distribution<Real>) -> Self {
1118        let mut rng = rand::thread_rng();
1119        EField::new_determinist(l, &mut rng, d)
1120    }
1121
1122    /// Create a new cold configuration for the electrical field, i.e. all E ar set to 0.
1123    ///
1124    /// # Example
1125    /// ```
1126    /// use lattice_qcd_rs::{field::EField, lattice::LatticeCyclic};
1127    /// # use std::error::Error;
1128    ///
1129    /// # fn main() -> Result<(), Box<dyn Error>> {
1130    /// let lattice = LatticeCyclic::<3>::new(1_f64, 4)?;
1131    /// let e_field = EField::new_cold(&lattice);
1132    /// assert!(!e_field.is_empty());
1133    /// # Ok(())
1134    /// # }
1135    /// ```
1136    pub fn new_cold(l: &LatticeCyclic<D>) -> Self {
1137        let p1 = Su3Adjoint::new_from_array([0_f64; 8]);
1138        Self {
1139            data: vec![SVector::<Su3Adjoint, D>::from_element(p1); l.number_of_points()],
1140        }
1141    }
1142
1143    /// Get `E(point) = [E_x(point), E_y(point), E_z(point)]`.
1144    pub fn e_vec(
1145        &self,
1146        point: &LatticePoint<D>,
1147        l: &LatticeCyclic<D>,
1148    ) -> Option<&SVector<Su3Adjoint, D>> {
1149        self.data.get(point.to_index(l))
1150    }
1151
1152    /// Get `E_{dir}(point)`. The sign of the direction does not change the output. i.e.
1153    /// `E_{-dir}(point) = E_{dir}(point)`.
1154    pub fn e_field(
1155        &self,
1156        point: &LatticePoint<D>,
1157        dir: &Direction<D>,
1158        l: &LatticeCyclic<D>,
1159    ) -> Option<&Su3Adjoint> {
1160        let value = self.e_vec(point, l);
1161        match value {
1162            Some(vec) => vec.get(dir.index()),
1163            None => None,
1164        }
1165    }
1166
1167    /// Returns wether there is no data
1168    pub fn is_empty(&self) -> bool {
1169        self.data.is_empty()
1170    }
1171
1172    /// Return the Gauss parameter `G(x) = \sum_i E_i(x) - U_{-i}(x) E_i(x - i) U^\dagger_{-i}(x)`.
1173    #[inline]
1174    pub fn gauss(
1175        &self,
1176        link_matrix: &LinkMatrix,
1177        point: &LatticePoint<D>,
1178        lattice: &LatticeCyclic<D>,
1179    ) -> Option<CMatrix3> {
1180        if lattice.number_of_points() != self.len()
1181            || lattice.number_of_canonical_links_space() != link_matrix.len()
1182        {
1183            return None;
1184        }
1185        Direction::positive_directions()
1186            .iter()
1187            .map(|dir| {
1188                let e_i = self.e_field(point, dir, lattice)?;
1189                let u_mi = link_matrix.matrix(&LatticeLink::new(*point, -*dir), lattice)?;
1190                let p_mi = lattice.add_point_direction(*point, &-dir);
1191                let e_m_i = self.e_field(&p_mi, dir, lattice)?;
1192                Some(e_i.to_matrix() - u_mi * e_m_i.to_matrix() * u_mi.adjoint())
1193            })
1194            .sum::<Option<CMatrix3>>()
1195    }
1196
1197    /// Get the deviation from the Gauss law
1198    #[inline]
1199    pub fn gauss_sum_div(
1200        &self,
1201        link_matrix: &LinkMatrix,
1202        lattice: &LatticeCyclic<D>,
1203    ) -> Option<Real> {
1204        if lattice.number_of_points() != self.len()
1205            || lattice.number_of_canonical_links_space() != link_matrix.len()
1206        {
1207            return None;
1208        }
1209        lattice
1210            .get_points()
1211            .par_bridge()
1212            .map(|point| {
1213                self.gauss(link_matrix, &point, lattice).map(|el| {
1214                    (su3::GENERATORS.iter().copied().sum::<CMatrix3>() * el)
1215                        .trace()
1216                        .abs()
1217                })
1218            })
1219            .sum::<Option<Real>>()
1220    }
1221
1222    /// project to that the gauss law is approximately respected ( up to `f64::EPSILON * 10` per point).
1223    ///
1224    /// It is mainly use internally but can be use to correct numerical drift in simulations.
1225    ///
1226    /// # Example
1227    /// ```
1228    /// use lattice_qcd_rs::error::ImplementationError;
1229    /// use lattice_qcd_rs::integrator::SymplecticEulerRayon;
1230    /// use lattice_qcd_rs::simulation::{
1231    ///     LatticeState, LatticeStateDefault, LatticeStateEFSyncDefault, LatticeStateWithEField,
1232    ///     SimulationStateSynchronous,
1233    /// };
1234    /// use rand::SeedableRng;
1235    ///
1236    /// # use std::error::Error;
1237    /// #
1238    /// # fn main() -> Result<(), Box<dyn Error>> {
1239    /// let mut rng = rand::rngs::StdRng::seed_from_u64(0); // change with your seed
1240    /// let distribution =
1241    ///     rand::distributions::Uniform::from(-std::f64::consts::PI..std::f64::consts::PI);
1242    /// let mut state = LatticeStateEFSyncDefault::new_random_e_state(
1243    ///     LatticeStateDefault::<3>::new_determinist(1_f64, 6_f64, 4, &mut rng)?,
1244    ///     &mut rng,
1245    /// ); // <- here internally when choosing randomly the EField it is projected.
1246    ///
1247    /// let integrator = SymplecticEulerRayon::default();
1248    /// for _ in 0..2 {
1249    ///     for _ in 0..10 {
1250    ///         state = state.simulate_sync(&integrator, 0.0001_f64)?;
1251    ///     }
1252    ///     // we correct the numerical drift of the EField.
1253    ///     let new_e_field = state
1254    ///         .e_field()
1255    ///         .project_to_gauss(state.link_matrix(), state.lattice())
1256    ///         .ok_or(ImplementationError::OptionWithUnexpectedNone)?;
1257    ///     state.set_e_field(new_e_field);
1258    /// }
1259    /// #
1260    /// #     Ok(())
1261    /// # }
1262    /// ```
1263    #[allow(clippy::as_conversions)] // no try into for f64
1264    #[inline]
1265    pub fn project_to_gauss(
1266        &self,
1267        link_matrix: &LinkMatrix,
1268        lattice: &LatticeCyclic<D>,
1269    ) -> Option<Self> {
1270        // TODO improve
1271        const NUMBER_FOR_LOOP: usize = 4;
1272
1273        if lattice.number_of_points() != self.len()
1274            || lattice.number_of_canonical_links_space() != link_matrix.len()
1275        {
1276            return None;
1277        }
1278        let mut return_val = self.project_to_gauss_step(link_matrix, lattice);
1279        loop {
1280            let val_dif = return_val.gauss_sum_div(link_matrix, lattice)?;
1281            //println!("diff : {}", val_dif);
1282            if val_dif.is_nan() {
1283                return None;
1284            }
1285            if val_dif <= f64::EPSILON * (lattice.number_of_points() * 4 * 8 * 10) as f64 {
1286                break;
1287            }
1288            for _ in 0_usize..NUMBER_FOR_LOOP {
1289                return_val = return_val.project_to_gauss_step(link_matrix, lattice);
1290                //println!("{}", return_val[0][0][0]);
1291            }
1292        }
1293        Some(return_val)
1294    }
1295
1296    /// Done one step to project to gauss law.
1297    ///
1298    /// # Panic
1299    /// panics if the link matrix and lattice is not of the correct size.
1300    #[inline]
1301    fn project_to_gauss_step(&self, link_matrix: &LinkMatrix, lattice: &LatticeCyclic<D>) -> Self {
1302        /// see <https://arxiv.org/pdf/1512.02374.pdf>
1303        // TODO verify
1304        const K: na::Complex<f64> = na::Complex::new(0.12_f64, 0_f64);
1305        let data = lattice
1306            .get_points()
1307            .collect::<Vec<LatticePoint<D>>>()
1308            .par_iter()
1309            .map(|point| {
1310                let e = self.e_vec(point, lattice).unwrap();
1311                SVector::<_, D>::from_fn(|index_dir, _| {
1312                    let dir = Direction::<D>::positive_directions()[index_dir];
1313                    let u = link_matrix
1314                        .matrix(&LatticeLink::new(*point, dir), lattice)
1315                        .unwrap();
1316                    let gauss = self.gauss(link_matrix, point, lattice).unwrap();
1317                    let gauss_p = self
1318                        .gauss(
1319                            link_matrix,
1320                            &lattice.add_point_direction(*point, &dir),
1321                            lattice,
1322                        )
1323                        .unwrap();
1324                    Su3Adjoint::new(Vector8::from_fn(|index, _| {
1325                        2_f64
1326                            * (su3::GENERATORS[index]
1327                                * ((u * gauss * u.adjoint() * gauss_p - gauss) * K
1328                                    + su3::GENERATORS[index]
1329                                        * na::Complex::from(e[dir.index()][index])))
1330                            .trace()
1331                            .real()
1332                    }))
1333                })
1334            })
1335            .collect();
1336        Self::new(data)
1337    }
1338}
1339
1340impl<const D: usize> AsRef<Vec<SVector<Su3Adjoint, D>>> for EField<D> {
1341    fn as_ref(&self) -> &Vec<SVector<Su3Adjoint, D>> {
1342        self.as_vec()
1343    }
1344}
1345
1346impl<const D: usize> AsMut<Vec<SVector<Su3Adjoint, D>>> for EField<D> {
1347    fn as_mut(&mut self) -> &mut Vec<SVector<Su3Adjoint, D>> {
1348        self.data_mut()
1349    }
1350}
1351
1352impl<const D: usize> AsRef<[SVector<Su3Adjoint, D>]> for EField<D> {
1353    fn as_ref(&self) -> &[SVector<Su3Adjoint, D>] {
1354        self.as_slice()
1355    }
1356}
1357
1358impl<const D: usize> AsMut<[SVector<Su3Adjoint, D>]> for EField<D> {
1359    fn as_mut(&mut self) -> &mut [SVector<Su3Adjoint, D>] {
1360        self.as_slice_mut()
1361    }
1362}
1363
1364impl<'a, const D: usize> IntoIterator for &'a EField<D> {
1365    type IntoIter = <&'a Vec<SVector<Su3Adjoint, D>> as IntoIterator>::IntoIter;
1366    type Item = &'a SVector<Su3Adjoint, D>;
1367
1368    fn into_iter(self) -> Self::IntoIter {
1369        self.data.iter()
1370    }
1371}
1372
1373impl<'a, const D: usize> IntoIterator for &'a mut EField<D> {
1374    type IntoIter = <&'a mut Vec<SVector<Su3Adjoint, D>> as IntoIterator>::IntoIter;
1375    type Item = &'a mut SVector<Su3Adjoint, D>;
1376
1377    fn into_iter(self) -> Self::IntoIter {
1378        self.data.iter_mut()
1379    }
1380}
1381
1382impl<const D: usize> Index<usize> for EField<D> {
1383    type Output = SVector<Su3Adjoint, D>;
1384
1385    fn index(&self, pos: usize) -> &Self::Output {
1386        &self.data[pos]
1387    }
1388}
1389
1390impl<const D: usize> IndexMut<usize> for EField<D> {
1391    fn index_mut(&mut self, pos: usize) -> &mut Self::Output {
1392        &mut self.data[pos]
1393    }
1394}
1395
1396impl<A, const D: usize> FromIterator<A> for EField<D>
1397where
1398    Vec<SVector<Su3Adjoint, D>>: FromIterator<A>,
1399{
1400    fn from_iter<T>(iter: T) -> Self
1401    where
1402        T: IntoIterator<Item = A>,
1403    {
1404        Self::new(Vec::from_iter(iter))
1405    }
1406}
1407
1408impl<A, const D: usize> FromParallelIterator<A> for EField<D>
1409where
1410    Vec<SVector<Su3Adjoint, D>>: FromParallelIterator<A>,
1411    A: Send,
1412{
1413    fn from_par_iter<I>(par_iter: I) -> Self
1414    where
1415        I: IntoParallelIterator<Item = A>,
1416    {
1417        Self::new(Vec::from_par_iter(par_iter))
1418    }
1419}
1420
1421impl<T, const D: usize> ParallelExtend<T> for EField<D>
1422where
1423    Vec<SVector<Su3Adjoint, D>>: ParallelExtend<T>,
1424    T: Send,
1425{
1426    fn par_extend<I>(&mut self, par_iter: I)
1427    where
1428        I: IntoParallelIterator<Item = T>,
1429    {
1430        self.data.par_extend(par_iter);
1431    }
1432}
1433
1434impl<A, const D: usize> Extend<A> for EField<D>
1435where
1436    Vec<SVector<Su3Adjoint, D>>: Extend<A>,
1437{
1438    fn extend<T>(&mut self, iter: T)
1439    where
1440        T: IntoIterator<Item = A>,
1441    {
1442        self.data.extend(iter);
1443    }
1444}
1445
1446#[cfg(test)]
1447mod test {
1448    use approx::*;
1449    use rand::SeedableRng;
1450
1451    use super::super::{lattice::*, Complex};
1452    use super::*;
1453
1454    const EPSILON: f64 = 0.000_000_001_f64;
1455    const SEED_RNG: u64 = 0x45_78_93_f4_4a_b0_67_f0;
1456
1457    #[test]
1458    fn test_get_e_field_pos_neg() {
1459        use super::super::lattice;
1460
1461        let l = LatticeCyclic::new(1_f64, 4).unwrap();
1462        let e = EField::new(vec![SVector::<_, 4>::from([
1463            Su3Adjoint::from([1_f64; 8]),
1464            Su3Adjoint::from([2_f64; 8]),
1465            Su3Adjoint::from([3_f64; 8]),
1466            Su3Adjoint::from([2_f64; 8]),
1467        ])]);
1468        assert_eq!(
1469            e.e_field(
1470                &LatticePoint::new([0, 0, 0, 0].into()),
1471                &lattice::DirectionEnum::XPos.into(),
1472                &l
1473            ),
1474            e.e_field(
1475                &LatticePoint::new([0, 0, 0, 0].into()),
1476                &lattice::DirectionEnum::XNeg.into(),
1477                &l
1478            )
1479        );
1480    }
1481
1482    #[test]
1483    #[allow(clippy::eq_op)]
1484    #[allow(clippy::op_ref)]
1485    fn test_su3_adj() {
1486        let mut rng = rand::rngs::StdRng::seed_from_u64(SEED_RNG);
1487        let d = rand::distributions::Uniform::from(-1_f64..1_f64);
1488        for _ in 0_u32..100_u32 {
1489            let v = Su3Adjoint::random(&mut rng, &d);
1490            let m = v.to_matrix();
1491            assert_abs_diff_eq!(
1492                v.trace_squared(),
1493                (m * m).trace().modulus(),
1494                epsilon = EPSILON
1495            );
1496            assert_eq_complex!(
1497                v.d(),
1498                nalgebra::Complex::new(0_f64, 1_f64) * m.determinant(),
1499                EPSILON
1500            );
1501            assert_eq_complex!(v.t(), -(m * m).trace() / Complex::from(2_f64), EPSILON);
1502
1503            // ----
1504            let adj_1 = Su3Adjoint::default();
1505            let adj_2 = Su3Adjoint::new_from_array([1_f64; 8]);
1506            assert_eq!(adj_2, adj_2 + adj_1);
1507            assert_eq!(adj_2, &adj_2 + &adj_1);
1508            assert_eq!(adj_2, &adj_2 - &adj_1);
1509            assert_eq!(adj_1, &adj_2 - &adj_2);
1510            assert_eq!(adj_1, &adj_2 - adj_2);
1511            assert_eq!(adj_1, adj_2 - &adj_2);
1512            assert_eq!(adj_1, -&adj_1);
1513            let adj_3 = Su3Adjoint::new_from_array([2_f64; 8]);
1514            assert_eq!(adj_3, &adj_2 + &adj_2);
1515            assert_eq!(adj_3, &adj_2 * &2_f64);
1516            assert_eq!(adj_3, &2_f64 * &adj_2);
1517            assert_eq!(adj_3, 2_f64 * adj_2);
1518            assert_eq!(adj_3, &2_f64 * adj_2);
1519            assert_eq!(adj_3, 2_f64 * &adj_2);
1520            assert_eq!(adj_2, &adj_3 / &2_f64);
1521            assert_eq!(adj_2, &adj_3 / 2_f64);
1522            let mut adj_5 = Su3Adjoint::new_from_array([2_f64; 8]);
1523            adj_5 /= &2_f64;
1524            assert_eq!(adj_2, adj_5);
1525            let adj_4 = Su3Adjoint::new_from_array([-1_f64; 8]);
1526            assert_eq!(adj_2, -adj_4);
1527        }
1528
1529        use crate::su3::su3_exp_r;
1530        let mut rng = rand::rngs::StdRng::seed_from_u64(SEED_RNG);
1531        let d = rand::distributions::Uniform::from(-1_f64..1_f64);
1532        for _ in 0_u32..10_u32 {
1533            let v = Su3Adjoint::random(&mut rng, &d);
1534            assert_eq!(su3_exp_r(v), v.exp());
1535        }
1536    }
1537
1538    #[test]
1539    fn link_matrix() {
1540        let lattice = LatticeCyclic::<3>::new(1_f64, 4).unwrap();
1541        match LinkMatrix::new_random_threaded(&lattice, 0) {
1542            Err(ThreadError::ThreadNumberIncorrect) => {}
1543            _ => panic!("unexpected output"),
1544        }
1545        let link_s = LinkMatrix::new_random_threaded(&lattice, 2);
1546        assert!(link_s.is_ok());
1547        let mut link = link_s.unwrap();
1548        assert!(!link.is_empty());
1549        let l2 = LinkMatrix::new(vec![]);
1550        assert!(l2.is_empty());
1551
1552        let _: &[_] = link.as_ref();
1553        let _: &Vec<_> = link.as_ref();
1554        let _: &mut [_] = link.as_mut();
1555        let _: &mut Vec<_> = link.as_mut();
1556        let _ = link.iter();
1557        let _ = link.iter_mut();
1558        let _ = (&link).into_iter();
1559        let _ = (&mut link).into_iter();
1560    }
1561
1562    #[test]
1563    fn e_field() {
1564        let lattice = LatticeCyclic::<3>::new(1_f64, 4).unwrap();
1565        let e_field_s = LinkMatrix::new_random_threaded(&lattice, 2);
1566        assert!(e_field_s.is_ok());
1567        let mut e_field = e_field_s.unwrap();
1568
1569        let _: &[_] = e_field.as_ref();
1570        let _: &Vec<_> = e_field.as_ref();
1571        let _: &mut [_] = e_field.as_mut();
1572        let _: &mut Vec<_> = e_field.as_mut();
1573        let _ = e_field.iter();
1574        let _ = e_field.iter_mut();
1575        let _ = (&e_field).into_iter();
1576        let _ = (&mut e_field).into_iter();
1577    }
1578
1579    #[test]
1580    fn magnetic_field() {
1581        let lattice = LatticeCyclic::<3>::new(1_f64, 4).unwrap();
1582        let mut link_matrix = LinkMatrix::new_cold(&lattice);
1583        let point = LatticePoint::from([0, 0, 0]);
1584        let dir_x = Direction::new(0, true).unwrap();
1585        let dir_y = Direction::new(1, true).unwrap();
1586        let dir_z = Direction::new(2, true).unwrap();
1587        let clover = link_matrix
1588            .clover(&point, &dir_x, &dir_y, &lattice)
1589            .unwrap();
1590        assert_eq_matrix!(CMatrix3::identity() * Complex::from(4_f64), clover, EPSILON);
1591        let f = link_matrix
1592            .f_mu_nu(&point, &dir_x, &dir_y, &lattice)
1593            .unwrap();
1594        assert_eq_matrix!(CMatrix3::zeros(), f, EPSILON);
1595        let b = link_matrix
1596            .magnetic_field(&point, &dir_x, &lattice)
1597            .unwrap();
1598        assert_eq_matrix!(CMatrix3::zeros(), b, EPSILON);
1599        let b_vec = link_matrix.magnetic_field_vec(&point, &lattice).unwrap();
1600        for i in &b_vec {
1601            assert_eq_matrix!(CMatrix3::zeros(), i, EPSILON);
1602        }
1603        // ---
1604        link_matrix[0] = CMatrix3::identity() * Complex::new(0_f64, 1_f64);
1605        let clover = link_matrix
1606            .clover(&point, &dir_x, &dir_y, &lattice)
1607            .unwrap();
1608        assert_eq_matrix!(
1609            CMatrix3::identity() * Complex::new(2_f64, 0_f64),
1610            clover,
1611            EPSILON
1612        );
1613        let clover = link_matrix
1614            .clover(&point, &dir_y, &dir_x, &lattice)
1615            .unwrap();
1616        assert_eq_matrix!(
1617            CMatrix3::identity() * Complex::new(2_f64, 0_f64),
1618            clover,
1619            EPSILON
1620        );
1621        let f = link_matrix
1622            .f_mu_nu(&point, &dir_x, &dir_y, &lattice)
1623            .unwrap();
1624        assert_eq_matrix!(
1625            CMatrix3::identity() * Complex::new(0_f64, 0_f64),
1626            f,
1627            EPSILON
1628        );
1629        let b = link_matrix
1630            .magnetic_field(&point, &dir_x, &lattice)
1631            .unwrap();
1632        assert_eq_matrix!(CMatrix3::zeros(), b, EPSILON);
1633        let b_vec = link_matrix.magnetic_field_vec(&point, &lattice).unwrap();
1634        for i in &b_vec {
1635            assert_eq_matrix!(CMatrix3::zeros(), i, EPSILON);
1636        }
1637        assert_eq_matrix!(
1638            link_matrix
1639                .magnetic_field_link(&LatticeLink::new(point, dir_x), &lattice)
1640                .unwrap(),
1641            b,
1642            EPSILON
1643        );
1644        //--
1645        let mut link_matrix = LinkMatrix::new_cold(&lattice);
1646        let link = LatticeLinkCanonical::new([1, 0, 0].into(), dir_y).unwrap();
1647        link_matrix[link.to_index(&lattice)] = CMatrix3::identity() * Complex::new(0_f64, 1_f64);
1648        let clover = link_matrix
1649            .clover(&point, &dir_x, &dir_y, &lattice)
1650            .unwrap();
1651        assert_eq_matrix!(
1652            CMatrix3::identity() * Complex::new(3_f64, 1_f64),
1653            clover,
1654            EPSILON
1655        );
1656        let clover = link_matrix
1657            .clover(&point, &dir_y, &dir_x, &lattice)
1658            .unwrap();
1659        assert_eq_matrix!(
1660            CMatrix3::identity() * Complex::new(3_f64, -1_f64),
1661            clover,
1662            EPSILON
1663        );
1664        let f = link_matrix
1665            .f_mu_nu(&point, &dir_x, &dir_y, &lattice)
1666            .unwrap();
1667        assert_eq_matrix!(
1668            CMatrix3::identity() * Complex::new(0_f64, 0.25_f64),
1669            f,
1670            EPSILON
1671        );
1672        let b = link_matrix
1673            .magnetic_field(&point, &dir_x, &lattice)
1674            .unwrap();
1675        assert_eq_matrix!(CMatrix3::zeros(), b, EPSILON);
1676        assert_eq_matrix!(
1677            link_matrix
1678                .magnetic_field_link(&LatticeLink::new(point, dir_x), &lattice)
1679                .unwrap(),
1680            b,
1681            EPSILON
1682        );
1683        let b = link_matrix
1684            .magnetic_field(&point, &dir_z, &lattice)
1685            .unwrap();
1686        assert_eq_matrix!(
1687            link_matrix
1688                .magnetic_field_link(&LatticeLink::new(point, dir_z), &lattice)
1689                .unwrap(),
1690            b,
1691            EPSILON
1692        );
1693        assert_eq_matrix!(
1694            CMatrix3::identity() * Complex::new(0.25_f64, 0_f64),
1695            b,
1696            EPSILON
1697        );
1698        let b_2 = link_matrix
1699            .magnetic_field(&[4, 0, 0].into(), &dir_z, &lattice)
1700            .unwrap();
1701        assert_eq_matrix!(b, b_2, EPSILON);
1702        let b_vec = link_matrix.magnetic_field_vec(&point, &lattice).unwrap();
1703        for (index, m) in b_vec.iter().enumerate() {
1704            if index == 2 {
1705                assert_eq_matrix!(m, b, EPSILON);
1706            }
1707            else {
1708                assert_eq_matrix!(CMatrix3::zeros(), m, EPSILON);
1709            }
1710        }
1711    }
1712
1713    #[test]
1714    fn test_len() {
1715        let su3 = Su3Adjoint::new(nalgebra::SVector::<f64, 8>::zeros());
1716        assert_eq!(Su3Adjoint::len_const(), su3.len());
1717    }
1718}
lattice_qcd_rs/field.rs

lattice_qcd_rs/
field.rs
