m4

src/Chemistry/Bonds/Implements/BaseBond.cs

@@ -0,0 +1,318 @@
+using Skeleton.Algebra.AdditionalProperties.Extensions;
+using Skeleton.Algebra.DimensionProviders;
+using Skeleton.Algebra.Extensions;
+using Skeleton.DataStructure;
+using VirtualChemistry.Chemistry.Atoms.Implements;
+using VirtualChemistry.Chemistry.Compounds.Implements;
+using VirtualChemistry.Chemistry.Mixtures.Implements;
+using VirtualChemistry.Constants;
+
+namespace VirtualChemistry.Chemistry.Bonds.Implements;
+/// <summary>
+/// meta bond
+/// </summary>
+public abstract class BaseBond
+{
+
+    /// <summary>
+    /// 
+    /// </summary>
+    protected BaseBond()
+    {
+        StateMatrix = new (_ =>
+            Connected ? EigenStateMatrix * ConnectedBond.EigenStateMatrix : EigenStateMatrix
+        );
+        LogKernel = new (_ => KernelStateMatrix.Log());
+        HalfBondInformation = new(x =>
+            LogKernel.GetFrom(x) + Atom.StructureMatrix.GetFrom(x).Log()
+        );
+        StructureMatrix = new(x =>
+            (HalfBondInformation.GetFrom(x) + (Connected ? ConnectedBond.HalfBondInformation.GetFrom(x) : LieSU3.Zero))
+            .Exp()
+        );
+        EnergyDistributionFactor = new(x =>
+            {
+                C3 v = StructureMatrix.GetFrom(x).Hermitian() * ChemistryConstant.BondEnergyDistributionVector;
+                return ChemistryConstant.BondEnergyDistributionSpectrum.Rayleigh<IDim3>(v);
+            }
+        );
+        LogConnectionMatrix = new(x => Connected
+            ? (AdjointStateMatrix * ConnectedBond.AdjointStateMatrix)
+            .ConjugateOn(Atom.LogKernelMatrix.GetFrom(x))
+            : LieSU3.Zero
+        );
+        LogState = new (x => StateMatrix.GetFrom(x).Log());
+        AntibondingCatalyzeFactor = new(x =>
+            {
+                if (Connected || HomogeneousMixture == HomogeneousMixture.Null)
+                    return 0;
+                U3 s = LogState.GetFrom(x).CayleyPositive();
+                LieSU3 f0 = HomogeneousMixture.LogStateMatrix.GetFrom(x).ConjugatedBy<IDim3, LieSU3>(s);
+                LieSU3 f1 = f0 ^ LogState.GetFrom(x);
+                U3 f2 = f1.CayleyNegative();
+                return f2.Det().Imaginary;
+            }
+        );
+        BondingCatalyzeFactor = new(x =>
+            {
+                if (!Connected || HomogeneousMixture == HomogeneousMixture.Null)
+                    return 0;
+                U3 s = LogState.GetFrom(x).CayleyNegative();
+                LieSU3 f0 = HomogeneousMixture.LogStateMatrix.GetFrom(x).ConjugatedBy<IDim3, LieSU3>(s);
+                LieSU3 f1 = f0 ^ LogState.GetFrom(x);
+                U3 f2 = f1.CayleyPositive();
+                return f2.Det().Imaginary;
+            }
+        );
+
+    }
+
+    /// <summary>
+    /// 
+    /// </summary>
+    public Compound Compound => Atom.Compound;
+
+    /// <summary>
+    /// 
+    /// </summary>
+    public HomogeneousMixture HomogeneousMixture => Compound.HomogeneousMixture; 
+    /// <summary>
+    /// null singleton
+    /// </summary>
+    public static readonly BaseBond Null = new NullBond();
+    /// <summary>
+    /// energy hold in this atom
+    /// </summary>
+    public double Energy { get; set; }
+
+    /// <summary>
+    /// how energy in the molecular is distributed into this atom
+    /// </summary>
+    public CacheItem<double> EnergyDistributionFactor { get; }
+
+    /// <summary>
+    /// the connected bond
+    /// </summary>
+    public BaseBond ConnectedBond { get; set; } = Null;
+    /// <summary>
+    /// atom that holds this bond
+    /// </summary>
+    public BaseAtom Atom { get; set; } = BaseAtom.Null;
+    /// <summary>
+    /// type of the bond
+    /// </summary>
+    public abstract string BondType { get; }
+    /// <summary>
+    /// kernel state matrix of the matrix, won't change as long as bond type is fixed
+    /// </summary>
+    public abstract SU3 KernelStateMatrix { get; }
+    /// <summary>
+    /// adj state matrix of the bond, won't change as long as bond type fixed
+    /// </summary>
+    public abstract SU3 AdjointStateMatrix { get; }
+    /// <summary>
+    /// eigen state matrix of bond, won't change as long as bond type fixed
+    /// </summary>
+    public abstract SU3 EigenStateMatrix { get; }
+
+    /// <summary>
+    /// log of the kernel state matrix
+    /// never expire
+    /// </summary>
+    public CacheItem<LieSU3> LogKernel { get; }
+
+    /// <summary>
+    /// skew hermitian matrix that represent the information of bond and its atom
+    /// </summary>
+    public CacheItem<LieSU3> HalfBondInformation { get; }
+
+    /// <summary>
+    /// special unitary matrix that represent the structure of the bond(its atom, connected bond and atom of connected bond)
+    /// </summary>
+    public CacheItem<SU3> StructureMatrix { get; }
+
+    /// <summary>
+    /// log of connection matrix
+    /// </summary>
+    public CacheItem<LieSU3> LogConnectionMatrix { get; set; }
+
+    /// <summary>
+    /// a special unitary matrix that represent the state of the bond
+    /// </summary>
+
+
+    public CacheItem<SU3> StateMatrix { get; }
+
+    /// <summary>
+    /// log of state matrix
+    /// </summary>
+    public CacheItem<LieSU3> LogState { get; }
+
+
+
+
+    /// <summary>
+    /// connect with another bond
+    /// </summary>
+    /// <param name="bond"></param>
+    /// <exception cref="Exception"></exception>
+    public void Connect(BaseBond bond)
+    {
+        if(Connected || bond.Connected)
+            return;
+        if(Atom.HomogeneousMixture != HomogeneousMixture.Null)
+            if (!Atom.Compound.Amount.IsEqualApprox(bond.Atom.Compound.Amount))
+                return;
+        if (bond == this)
+            throw new Exception("TODO");
+        ConnectedBond = bond;
+        bond.ConnectedBond = this;
+        StructureChanged();
+        bond.StructureChanged();
+        if (Atom.Compound.Atoms.Contains(bond.Atom))
+            return;
+        bond.Atom.HomogeneousMixture.Compounds.Remove(bond.Atom.Compound);
+        bond.Atom.Compound.Amount = 0d;
+        
+        foreach (BaseAtom atom in bond.Atom.Compound.Atoms.ToHashSet())
+        {
+            Atom.Compound.Atoms.Add(atom);
+            atom.Compound = Atom.Compound;
+        }
+    }
+    /// <summary>
+    /// disconnect with connected bond
+    /// </summary>
+    public void Disconnect()
+    {
+        if (!Connected)
+            return;
+        BaseBond disconnectedBond = ConnectedBond;
+        ConnectedBond = BaseBond.Null;
+        disconnectedBond.ConnectedBond = BaseBond.Null;
+        StructureChanged();
+        disconnectedBond.StructureChanged();
+        if (Atom.HasPathTo(disconnectedBond.Atom))
+            return;
+        Compound newCompound = disconnectedBond.Atom.GrabCompound;
+        newCompound.Amount = Atom.Compound.Amount;
+        
+        /*
+        newCompound.HomogeneousMixture = Atom.HomogeneousMixture;*/
+        foreach (BaseAtom atom in newCompound.Atoms)
+        {
+            atom.Compound = newCompound;
+            if (Atom.Compound.Atoms.Contains(atom))
+                Atom.Compound.Atoms.Remove(atom);
+        }
+        Atom.HomogeneousMixture.AddCompound(newCompound);
+    }
+    /// <summary>
+    /// if the bond is connected
+    /// </summary>
+    public bool Connected => ConnectedBond != Null;
+    /// <summary>
+    /// bond number, redundant information of bond type
+    /// </summary>
+    public abstract int BondNumber { get; }
+    private void StructureChanged()
+    {
+        //Atom.Compound.ResetDistanceCache();
+        Atom.Compound.CacheInit();
+        HalfBondInformation.Expire();
+        StateMatrix.Expire();
+        LogConnectionMatrix.Expire();
+        
+    }
+    
+    
+    /// <summary>
+    /// Conditions for a pair of unconnected bonds tp connect:
+    /// 1. Energy difference of the two bond is less than a threshold
+    /// 2. Energy in both bonds is higher than bonding energy
+    ///
+    /// if bond is connected, the antibonding energy is NaN
+    /// </summary>
+    /// <returns></returns>
+    public double EigenBondingEnergy()
+    {
+        if (Connected)
+            return double.NaN;
+        SU3 s1 = StructureMatrix.Get * Atom.StructureMatrix.Get;
+        C33 s2 = StructureMatrix.Get + Atom.StructureMatrix.Get;
+        double lowerLimit = -Math.PI;
+        double upperLimit = Math.PI;
+        C3 coVector = s2.Hermitian().ColumnAverage;
+        return s1
+            .ConjugateOn(new DiagR3(lowerLimit, 0d, upperLimit))
+            .Rayleigh(coVector);
+        
+    }
+    
+    /// <summary>
+    /// Conditions for a connected bond to break:
+    /// 1. energy difference of two connected bond is greater than a threshold
+    /// 2. energy in either bond is higher than antibonding energy
+    ///
+    /// if bond is unconnected, the antibonding energy is 
+    /// </summary>
+    /// <returns></returns>
+    public double EigenAntibondingEnergy()
+    {
+        if(!Connected) 
+            return double.NaN;
+        SU3 s1 = Atom.StructureMatrix.Get.Hermitian() * StructureMatrix.Get;
+        SU3 s2 = ConnectedBond.Atom.StructureMatrix.Get.Hermitian() * ConnectedBond.StructureMatrix.Get;
+        C33 sw = s1 + s2;
+        SU3 w = (s1.Log() + s2.Log()).Exp();
+        C3 coVector = sw.ColumnAverageBivariant;
+        double upperLimit = -Math.PI;
+        double lowerLimit = Math.PI;
+        return w
+            .ConjugateOn(new DiagR3(lowerLimit, 0d, upperLimit))
+            .Rayleigh(coVector);
+    }
+
+    /// <summary>
+    /// how antibonding energy affected by mixture
+    /// </summary>
+    public CacheItem<double> AntibondingCatalyzeFactor { get; set; }
+
+    /// <summary>
+    /// hot bonding energy affected by chemical environment
+    /// </summary>
+    public CacheItem<double> BondingCatalyzeFactor { get; set; }
+
+    public double BondingEnergy => EigenBondingEnergy() + BondingCatalyzeFactor.Get;
+    public double AntiBondingEnergy => EigenAntibondingEnergy() + AntibondingCatalyzeFactor.Get;
+
+    /// <summary>
+    /// how to pair bonds
+    /// </summary>
+    public int BondingGroup
+    {
+        get
+        {
+            U3 sP = HalfBondInformation.Get.CayleyPositive();
+            U3 eN = HomogeneousMixture.LogStateMatrix.Get.CayleyNegative();
+            DiagR3 h = new(0, Math.PI, Math.PI * 2);
+            U3 w = sP * eN;
+            double t = w.ConjugateOn(h).Rayleigh(sP.ColumnAverageBivariant);
+            /*
+            //Complex t = (sP * eN).Trace();// + eN.Trace() * eN.Trace();
+            C33 a1 = sP * eN;
+            C33 a2 = eN * sP;
+            C33 w = a1 - a2;
+            Complex t = (w).Det() + a1.Trace() + a2.Trace();// * (eN*sP ).Trace();// * EigenStateMatrix.Trace();*/
+            //Console.WriteLine($"TRACE  {t} -");
+
+            return t switch
+            {
+                > 7 * Math.PI / 4 or < Math.PI / 4 => 0,
+                > Math.PI / 4 and < 3 * Math.PI / 4 => 1,
+                > 3 * Math.PI / 4 and < 5 * Math.PI / 4 => 2,
+                _ => 3
+            };
+        }
+    }
+}

src/Chemistry/Bonds/Implements/DBond.cs

@@ -0,0 +1,23 @@
+using VirtualChemistry.Constants;
+
+namespace VirtualChemistry.Chemistry.Bonds.Implements;
+/// <summary>
+/// bond # 3
+/// </summary>
+public class DBond : BaseBond
+{
+    /// <inheritdoc />
+    public override string BondType => ChemistryConstant.BondTypes.D;
+
+    /// <inheritdoc />
+    public override int BondNumber => 3;
+
+    /// <inheritdoc />
+    public override SU3 KernelStateMatrix => ChemistryConstant.BondKernels.BondKernelD;
+
+    /// <inheritdoc />
+    public override SU3 AdjointStateMatrix => ChemistryConstant.BondAdjoints.BondAdjointD;
+
+    /// <inheritdoc />
+    public override SU3 EigenStateMatrix => ChemistryConstant.BondEigenStates.BondEigenStateD;
+}

src/Chemistry/Bonds/Implements/HBond.cs

@@ -0,0 +1,24 @@
+using VirtualChemistry.Constants;
+
+namespace VirtualChemistry.Chemistry.Bonds.Implements;
+
+/// <summary>
+/// bond # 7
+/// </summary>
+public class HBond : BaseBond
+{
+    /// <inheritdoc />
+    public override string BondType => ChemistryConstant.BondTypes.H;
+
+    /// <inheritdoc />
+    public override int BondNumber => 7;
+
+    /// <inheritdoc />
+    public override SU3 KernelStateMatrix => ChemistryConstant.BondKernels.BondKernelH;
+
+    /// <inheritdoc />
+    public override SU3 AdjointStateMatrix => ChemistryConstant.BondAdjoints.BondAdjointH;
+
+    /// <inheritdoc />
+    public override SU3 EigenStateMatrix => ChemistryConstant.BondEigenStates.BondEigenStateH;
+}

src/Chemistry/Bonds/Implements/MBond.cs

@@ -0,0 +1,25 @@
+using VirtualChemistry.Constants;
+
+namespace VirtualChemistry.Chemistry.Bonds.Implements;
+
+/// <summary>
+/// bond # 1
+/// </summary>
+public class MBond : BaseBond
+{
+    /// <inheritdoc />
+    public override string BondType => ChemistryConstant.BondTypes.M;
+
+    /// <inheritdoc />
+    public override int BondNumber => 1;
+
+    /// <inheritdoc />
+    public override SU3 KernelStateMatrix => ChemistryConstant.BondKernels.BondKernelM;
+
+    /// <inheritdoc />
+    public override SU3 AdjointStateMatrix => ChemistryConstant.BondAdjoints.BondAdjointM;
+
+    /// <inheritdoc />
+    public override SU3 EigenStateMatrix => ChemistryConstant.BondEigenStates.BondEigenStateM;
+
+}

src/Chemistry/Bonds/Implements/NullBond.cs

@@ -0,0 +1,27 @@
+namespace VirtualChemistry.Chemistry.Bonds.Implements;
+
+/// <summary>
+/// bond singleton that represent no bond
+/// </summary>
+public class NullBond : BaseBond
+{
+    /// <inheritdoc />
+    public NullBond()
+    {
+    }
+
+    /// <inheritdoc />
+    public override string BondType => throw new Exception();
+
+    /// <inheritdoc />
+    public override SU3 KernelStateMatrix => throw new Exception();
+
+    /// <inheritdoc />
+    public override SU3 AdjointStateMatrix => throw new Exception();
+
+    /// <inheritdoc />
+    public override SU3 EigenStateMatrix => throw new Exception();
+
+    /// <inheritdoc />
+    public override int BondNumber => throw new Exception();
+}

src/Chemistry/Bonds/Implements/PBond.cs

@@ -0,0 +1,24 @@
+using VirtualChemistry.Constants;
+
+namespace VirtualChemistry.Chemistry.Bonds.Implements;
+
+/// <summary>
+/// bond # 6
+/// </summary>
+public class PBond : BaseBond
+{
+    /// <inheritdoc />
+    public override string BondType => ChemistryConstant.BondTypes.P;
+
+    /// <inheritdoc />
+    public override int BondNumber => 6;
+
+    /// <inheritdoc />
+    public override SU3 KernelStateMatrix => ChemistryConstant.BondKernels.BondKernelP;
+
+    /// <inheritdoc />
+    public override SU3 AdjointStateMatrix => ChemistryConstant.BondAdjoints.BondAdjointP;
+
+    /// <inheritdoc />
+    public override SU3 EigenStateMatrix => ChemistryConstant.BondEigenStates.BondEigenStateP;
+}

src/Chemistry/Bonds/Implements/QBond.cs

@@ -0,0 +1,25 @@
+using VirtualChemistry.Constants;
+
+namespace VirtualChemistry.Chemistry.Bonds.Implements;
+
+/// <summary>
+/// bond # 5
+/// </summary>
+public class QBond : BaseBond
+{
+    /// <inheritdoc />
+    public override string BondType => ChemistryConstant.BondTypes.Q;
+
+    /// <inheritdoc />
+    public override int BondNumber => 5;
+
+    /// <inheritdoc />
+    public override SU3 KernelStateMatrix => ChemistryConstant.BondKernels.BondKernelQ;
+
+    /// <inheritdoc />
+    public override SU3 AdjointStateMatrix => ChemistryConstant.BondAdjoints.BondAdjointQ;
+
+    /// <inheritdoc />
+    public override SU3 EigenStateMatrix => ChemistryConstant.BondEigenStates.BondEigenStateQ;
+    
+}

src/Chemistry/Bonds/Implements/SBond.cs

@@ -0,0 +1,24 @@
+using VirtualChemistry.Constants;
+
+namespace VirtualChemistry.Chemistry.Bonds.Implements;
+
+/// <summary>
+/// bond # 2
+/// </summary>
+public class SBond : BaseBond
+{
+    /// <inheritdoc />
+    public override string BondType => ChemistryConstant.BondTypes.S;
+
+    /// <inheritdoc />
+    public override int BondNumber => 2;
+
+    /// <inheritdoc />
+    public override SU3 KernelStateMatrix => ChemistryConstant.BondKernels.BondKernelS;
+
+    /// <inheritdoc />
+    public override SU3 AdjointStateMatrix => ChemistryConstant.BondAdjoints.BondAdjointS;
+
+    /// <inheritdoc />
+    public override SU3 EigenStateMatrix => ChemistryConstant.BondEigenStates.BondEigenStateS;
+}

src/Chemistry/Bonds/Implements/YBond.cs

@@ -0,0 +1,24 @@
+using VirtualChemistry.Constants;
+
+namespace VirtualChemistry.Chemistry.Bonds.Implements;
+
+/// <summary>
+/// bond # 4
+/// </summary>
+public class YBond : BaseBond
+{
+    /// <inheritdoc />
+    public override string BondType => ChemistryConstant.BondTypes.Y;
+
+    /// <inheritdoc />
+    public override int BondNumber => 4;
+
+    /// <inheritdoc />
+    public override SU3 KernelStateMatrix => ChemistryConstant.BondKernels.BondKernelY;
+
+    /// <inheritdoc />
+    public override SU3 AdjointStateMatrix => ChemistryConstant.BondAdjoints.BondAdjointY;
+
+    /// <inheritdoc />
+    public override SU3 EigenStateMatrix => ChemistryConstant.BondEigenStates.BondEigenStateY;
+}

src/Chemistry/Bonds/Resolver/BondResolver.cs

@@ -0,0 +1,53 @@
+using VirtualChemistry.Chemistry.Bonds.Implements;
+using VirtualChemistry.Constants;
+
+namespace VirtualChemistry.Chemistry.Bonds.Resolver;
+
+/// <summary>
+/// create a bond instance by descriptions 
+/// </summary>
+public static class BondResolver
+{
+    /// <summary>
+    /// create the bond by type
+    /// </summary>
+    /// <param name="type"></param>
+    /// <returns></returns>
+    /// <exception cref="Exception"></exception>
+    public static BaseBond Resolve(string type)
+    {
+        BaseBond res = type switch
+        {
+            ChemistryConstant.BondTypes.M => new MBond(),
+            ChemistryConstant.BondTypes.S => new SBond(),
+            ChemistryConstant.BondTypes.D => new DBond(),
+            ChemistryConstant.BondTypes.Y => new YBond(),
+            ChemistryConstant.BondTypes.Q => new QBond(),
+            ChemistryConstant.BondTypes.P => new PBond(),
+            ChemistryConstant.BondTypes.H => new HBond(),
+            _ => throw new Exception("TODO")
+        };
+        return res;
+    }
+    /// <summary>
+    /// create a bond by number
+    /// </summary>
+    /// <param name="number"></param>
+    /// <returns></returns>
+    /// <exception cref="Exception"></exception>
+    public static BaseBond Resolve(int number)
+    {
+        BaseBond res = number switch
+        {
+            1 => new MBond(),
+            2 => new SBond(),
+            3 => new DBond(),
+            4 => new YBond(),
+            5 => new QBond(),
+            6 => new PBond(),
+            7 => new HBond(),
+            _ => throw new Exception("TODO")
+        };
+        return res;
+    }
+}