commit 6242ed2cf1b504ab6d424a5701a739f6aa53f13c Author: hzhang Date: Fri Jun 21 23:42:10 2024 +0800 m4 diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..add57be --- /dev/null +++ b/.gitignore @@ -0,0 +1,5 @@ +bin/ +obj/ +/packages/ +riderModule.iml +/_ReSharper.Caches/ \ No newline at end of file diff --git a/VirtualChemistry.csproj b/VirtualChemistry.csproj new file mode 100644 index 0000000..3b46ad4 --- /dev/null +++ b/VirtualChemistry.csproj @@ -0,0 +1,13 @@ + + + + net6.0 + enable + enable + + + + + + + diff --git a/VirtualChemistry.sln b/VirtualChemistry.sln new file mode 100644 index 0000000..570f21f --- /dev/null +++ b/VirtualChemistry.sln @@ -0,0 +1,27 @@ + +Microsoft Visual Studio Solution File, Format Version 12.00 +Project("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "VirtualChemistry", "VirtualChemistry.csproj", "{36979A00-F434-4001-A80E-ECBF92A0AB5E}" +EndProject +Project("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "Skeleton", "..\RiderProjects\Skeleton\Skeleton.csproj", "{CB0ED85F-CAC0-4B53-96A8-16048A928890}" +EndProject +EndProject +Global + GlobalSection(SolutionConfigurationPlatforms) = preSolution + Debug|Any CPU = Debug|Any CPU + Release|Any CPU = Release|Any CPU + EndGlobalSection + GlobalSection(ProjectConfigurationPlatforms) = postSolution + {36979A00-F434-4001-A80E-ECBF92A0AB5E}.Debug|Any CPU.ActiveCfg = Debug|Any CPU + {36979A00-F434-4001-A80E-ECBF92A0AB5E}.Debug|Any CPU.Build.0 = Debug|Any CPU + {36979A00-F434-4001-A80E-ECBF92A0AB5E}.Release|Any CPU.ActiveCfg = Release|Any CPU + {36979A00-F434-4001-A80E-ECBF92A0AB5E}.Release|Any CPU.Build.0 = Release|Any CPU + {917762E5-D59D-4724-8F1A-73F739573E99}.Debug|Any CPU.ActiveCfg = Debug|Any CPU + {917762E5-D59D-4724-8F1A-73F739573E99}.Debug|Any CPU.Build.0 = Debug|Any CPU + {917762E5-D59D-4724-8F1A-73F739573E99}.Release|Any CPU.ActiveCfg = Release|Any CPU + {917762E5-D59D-4724-8F1A-73F739573E99}.Release|Any CPU.Build.0 = Release|Any CPU + {CB0ED85F-CAC0-4B53-96A8-16048A928890}.Debug|Any CPU.ActiveCfg = Debug|Any CPU + {CB0ED85F-CAC0-4B53-96A8-16048A928890}.Debug|Any CPU.Build.0 = Debug|Any CPU + {CB0ED85F-CAC0-4B53-96A8-16048A928890}.Release|Any CPU.ActiveCfg = Release|Any CPU + {CB0ED85F-CAC0-4B53-96A8-16048A928890}.Release|Any CPU.Build.0 = Release|Any CPU + EndGlobalSection +EndGlobal diff --git a/global.json b/global.json new file mode 100644 index 0000000..1bcf6c0 --- /dev/null +++ b/global.json @@ -0,0 +1,7 @@ +{ + "sdk": { + "version": "6.0.0", + "rollForward": "latestMinor", + "allowPrerelease": false + } +} \ No newline at end of file diff --git a/src/Chemistry/Atoms/Implements/AAtom.cs b/src/Chemistry/Atoms/Implements/AAtom.cs new file mode 100644 index 0000000..4668dc3 --- /dev/null +++ b/src/Chemistry/Atoms/Implements/AAtom.cs @@ -0,0 +1,26 @@ +using VirtualChemistry.Constants; + +namespace VirtualChemistry.Chemistry.Atoms.Implements; +/// +/// atom 12 +/// +public class AAtom : BaseAtom +{ + /// + public override string Element => ChemistryConstant.Elements.A; + + /// + public override int AtomicNumber => 12; + + /// + public override SU3 AdjointState => ChemistryConstant.AtomAdjoints.AtomAdjointStateA; + + /// + public override SU3 KernelState => ChemistryConstant.AtomKernels.AtomKernelStateA; + + /// + public override SU3 EigenStateMatrix => ChemistryConstant.AtomEigenStates.AtomEigenStateA; + + /// + public override double AtomicMass => ChemistryConstant.AtomicMasses.AtomicMassA; +} \ No newline at end of file diff --git a/src/Chemistry/Atoms/Implements/BaseAtom.cs b/src/Chemistry/Atoms/Implements/BaseAtom.cs new file mode 100644 index 0000000..b5cb03d --- /dev/null +++ b/src/Chemistry/Atoms/Implements/BaseAtom.cs @@ -0,0 +1,386 @@ +using System.Text; +using Skeleton.Algebra.AdditionalProperties.Extensions; +using Skeleton.Algebra.DimensionProviders; +using Skeleton.Constants; +using Skeleton.DataStructure; +using Skeleton.Utils.Helpers; +using VirtualChemistry.Chemistry.Bonds.Implements; +using VirtualChemistry.Chemistry.Bonds.Resolver; +using VirtualChemistry.Chemistry.Compounds.Implements; +using VirtualChemistry.Chemistry.Mixtures.Implements; +using VirtualChemistry.Constants; +using VirtualChemistry.DataStructure.Packs; + + +namespace VirtualChemistry.Chemistry.Atoms.Implements; +/// +/// meta atom +/// +public abstract class BaseAtom +{ + /// + /// + /// + protected BaseAtom() + { + LogEigenState = new(_ => EigenStateMatrix.Log()); + BondInformationMatrix = new( + x => Bonds + .Select(bond => bond.LogState.GetFrom(x)) + .DefaultIfEmpty(LieSU3.Zero) + .SpLieSum() + .Exp() + ); + + LogKernelMatrix = new(_ => KernelState.Log()); + StateMatrix = new(x => BondInformationMatrix.GetFrom(x).ConjugateOn(EigenStateMatrix)); + LogStateMatrix = new(x => BondInformationMatrix.GetFrom(x).ConjugateOn(LogEigenState.GetFrom(x))); + StructureMatrix = new(x => + { + var h = Bonds + .Select(b => b.LogConnectionMatrix.GetFrom(x)) + .DefaultIfEmpty(LieSU3.Zero) + .Union(new[] { LogStateMatrix.GetFrom(x), LogKernelMatrix.GetFrom(x) }); + var s = h.SpLieSum().Exp(); + return s; + } + ); + + EnergyDistributionFactor = new CacheItem( + x => + { + C3 v = StructureMatrix.GetFrom(x).Hermitian() * ChemistryConstant.AtomEnergyDistributionVector; + return ChemistryConstant.AtomEnergyDistributionSpectrum.Rayleigh(v); + } + ); + } + + /// + /// singleton to represent no atom + /// + public static readonly BaseAtom Null = new NullAtom(); + /// + /// energy holds in the atom + /// + public double Energy + { + get => Bonds.Sum(bond => bond.Energy); + set + { + double energyFactor = Bonds.Sum(bond => bond.EnergyDistributionFactor.Get); + foreach (BaseBond bond in Bonds) + bond.Energy = value * bond.EnergyDistributionFactor.Get / energyFactor; + } + } + + /// + /// describe how energy in molecular is distributed in this atom + /// + + + public CacheItem EnergyDistributionFactor { get; } + + /// + /// element name/ atom name + /// + public abstract string Element { get; } + /// + /// atomic number, redundant information as element + /// + public abstract int AtomicNumber { get; } + /// + /// atomic mass, wont change as long as atomic number is fixed + /// + public abstract double AtomicMass { get; } + /// + /// compound of this atom + /// + public Compound Compound { get; set; } = Compound.Null; + /// + /// homogeneous mixture of its compound + /// + public HomogeneousMixture HomogeneousMixture => Compound.HomogeneousMixture; + /// + /// bonds + /// + public HashSet Bonds { get; set; } = new HashSet(); + /// + /// bonds that already connected with other bond + /// + public IEnumerable ConnectedBonds => Bonds.Where(x => x.Connected); + + + /// + /// bonds that not connected with any other bond yet + /// + public IEnumerable UnconnectedBonds => Bonds.Where(x => !x.Connected); + /// + /// atoms that connected to this atom directly + /// + public IEnumerable ConnectedAtoms => ConnectedBonds + .Select(bond => bond.ConnectedBond.Atom); + /// + /// mass of all directly connected atoms + /// + public double GroupMass => ConnectedAtoms + .Select(atom => atom.AtomicMass) + .DefaultIfEmpty(0d) + .Sum(); + /// + /// if there is a path to specific atom + /// private use only + /// + /// + /// + /// + private bool HasPathTo(BaseAtom atom, HashSet visited) + { + if (atom == this) + return true; + if (ConnectedAtoms.Contains(atom)) + return true; + visited.Add(this); + return ConnectedAtoms + .Where(connectedAtom => !visited.Contains(connectedAtom)) + .Any(connectedAtom => connectedAtom.HasPathTo(atom, visited)); + } + /// + /// if there is a path to a specific atom + /// + /// + /// + public bool HasPathTo(BaseAtom atom) => HasPathTo(atom, new HashSet()); + /// + /// get the compound by connections + /// + public Compound GrabCompound + { + get + { + Compound res = new (); + res.Atoms.Add(this); + int aCount = 0; + while (aCount != res.Atoms.Count) + { + aCount = res.Atoms.Count; + foreach (BaseAtom a in res.Atoms + .SelectMany(a => a.ConnectedAtoms) + .Where(a => !res.Atoms.Contains(a)) + .ToArray() + ) + res.Atoms.Add(a); + } + + res.HomogeneousMixture = HomogeneousMixture.Null; + res.Amount = 0d; + return res; + } + } + /// + /// how far is an atom from this one(directly connected atoms has distance 1) + /// + /// + /// + /// + private int DistanceTo(BaseAtom atom, HashSet visited) + { + if (Compound != atom.Compound) + return -1; + if (atom == this) + return 0; + if(Compound.CachedDistance(this, atom) > 0) + return Compound.CachedDistance(this, atom); + if (ConnectedAtoms.Contains(atom)) + return 1; + visited.Add(this); + int distance = ConnectedAtoms + .Where(connectedAtom => !visited.Contains(connectedAtom)) + .Select(connectedAtom => connectedAtom.DistanceTo(atom, visited)) + .DefaultIfEmpty(int.MaxValue) + .Min(); + Compound.CacheDistance(this, atom, distance + 1); + return distance + 1; + } + /// + /// how far is an atom from this one(directly connected atoms has distance 1) + /// + /// + /// + public int DistanceTo(BaseAtom atom) => DistanceTo(atom, new HashSet()); + + /// + /// all connections of this atom + /// + public IEnumerable Connections => + Bonds + .Where(bond => bond.Connected) + .Select(bond => new ConnectionInfo(bond, bond.ConnectedBond, bond.ConnectedBond.Atom)); + /// + /// distance info as map + /// + /// + private Dictionary> DistanceMap() + { + Dictionary> res = new Dictionary>(); + foreach (BaseAtom atom in Compound.Atoms) + { + int dist = DistanceTo(atom); + if(!res.ContainsKey(dist)) + res.Add(dist, new HashSet()); + res[dist].Add(atom); + } + return res; + } + + /// + /// order bonds by selectors + /// + public IOrderedEnumerable OrderedBonds => + Bonds + .OrderBy(bond => bond.Connected) + .ThenBy(bond => bond.BondNumber) + .ThenBy(bond => bond.Connected ? bond.ConnectedBond.BondNumber : 0) + .ThenBy(bond => bond.Connected ? bond.ConnectedBond.Atom.AtomicNumber : 0); + /// + /// kernel state matrix of the atom, wont change + /// + public abstract SU3 KernelState { get; } + /// + /// adj matrix wont change + /// + public abstract SU3 AdjointState { get; } + /// + /// eig matrix wont change + /// + public abstract SU3 EigenStateMatrix { get; } + + /// + /// will never expire + /// + public CacheItem LogKernelMatrix { get; } + + /// + /// + /// + public CacheItem BondInformationMatrix { get; } + + /// + /// will never expire + /// + public CacheItem LogEigenState { get; } + + /// + /// state matrix + /// + public CacheItem StateMatrix { get; } + + /// + /// log state matrix + /// + public CacheItem LogStateMatrix { get; } + + /// + /// get an unconnected bond of specific type + /// + /// + /// + public BaseBond NextUnconnected(string bondType) => + Bonds + .Where(bond => bond.BondType == bondType) + .FirstOrDefault(bond => !bond.Connected, BaseBond.Null); + /// + /// get an unconnected bond of specific type, except for bonds in exclude set + /// + /// + /// + /// + public BaseBond NextUnconnected(string bondType, HashSet exclude) => + Bonds + .Where(bond => bond.BondType == bondType) + .Where(bond => !exclude.Contains(bond)) + .FirstOrDefault(bond => !bond.Connected, BaseBond.Null); + /// + /// get an unconnected bond of specific type that is not the excluded bond + /// + /// + /// + /// + public BaseBond NextUnconnected(string bondType, BaseBond exclude) => + Bonds + .Where(bond => bond.BondType == bondType) + .Where(bond => exclude != bond) + .FirstOrDefault(bond => !bond.Connected, BaseBond.Null); + + /// + /// structure matrix + /// + public CacheItem StructureMatrix { get; } + + /// + /// the representation of the compound relative to a center atom + /// + /// + /// + /// + public StringBuilder IsoRepresentation(BaseAtom center, HashSet visited) + { + HashSet toVisit = Bonds + .Where(b => b.Connected) + .Where(b => !visited.Contains(b)) + .Where(b => !visited.Contains(b.ConnectedBond)) + .ToHashSet(); + foreach (BaseBond b in toVisit) + { + visited.Add(b); + visited.Add(b.ConnectedBond); + + } + + Dictionary subRepresentations = new (); + foreach (BaseBond b in toVisit) + subRepresentations[b] = + b.ConnectedBond.Atom.IsoRepresentation(center, new HashSet(visited)); + + StringBuilder res = new StringBuilder($"{Element}#{DistanceTo(center)}-["); + IEnumerable visitOrder = toVisit + .GroupedOrderBy(bond => bond.BondNumber) + .GroupedThenBy(bond => bond.ConnectedBond.BondNumber) + .GroupedThenBy(bond => bond.ConnectedBond.Atom.AtomicNumber) + .GroupedThenBy(bond => bond.ConnectedBond.Atom.ConnectedAtoms.Count()) + .GroupedThenBy(bond => subRepresentations[bond].ToString()) + .SelectMany(g => g); + foreach (BaseBond b in visitOrder) + { + res.Append($"<{b.BondType}={b.ConnectedBond.BondType}*"); + res.Append(subRepresentations[b]); + res.Append(">"); + } + res.Append("]"); + + return res; + } + /// + /// the representation of the compound with this atom as center + /// + public string FullRepresentation => IsoRepresentation(this, new HashSet()).ToString(); + /// + /// initializer + /// + public void C() + { + Bonds = new HashSet(); + for (int bondNumber = 1; bondNumber <= 7; bondNumber++) + { + int xDiff = AtomicNumber - AlgebraConstant.Ludic(bondNumber); + if (AtomicNumber % AlgebraConstant.Ludic(bondNumber) != 0 && + (xDiff < 0 || xDiff % 7 != 0)) + continue; + for (int idx = 0; idx < ChemistryConstant.BondMaxNumbers(bondNumber); idx++) + Bonds.Add(BondResolver.Resolve(bondNumber)); + } + + foreach (BaseBond b in Bonds) + b.Atom = this; + } + +} diff --git a/src/Chemistry/Atoms/Implements/CxAtom.cs b/src/Chemistry/Atoms/Implements/CxAtom.cs new file mode 100644 index 0000000..216b37f --- /dev/null +++ b/src/Chemistry/Atoms/Implements/CxAtom.cs @@ -0,0 +1,27 @@ +using VirtualChemistry.Constants; + +namespace VirtualChemistry.Chemistry.Atoms.Implements; + +/// +/// atom # 13 +/// +public class CxAtom : BaseAtom +{ + /// + public override string Element => ChemistryConstant.Elements.Cx; + + /// + public override int AtomicNumber => 13; + + /// + public override SU3 AdjointState => ChemistryConstant.AtomAdjoints.AtomAdjointStateCx; + + /// + public override SU3 KernelState => ChemistryConstant.AtomKernels.AtomKernelStateCx; + + /// + public override SU3 EigenStateMatrix => ChemistryConstant.AtomEigenStates.AtomEigenStateCx; + + /// + public override double AtomicMass => ChemistryConstant.AtomicMasses.AtomicMassCx; +} \ No newline at end of file diff --git a/src/Chemistry/Atoms/Implements/DAtom.cs b/src/Chemistry/Atoms/Implements/DAtom.cs new file mode 100644 index 0000000..11d69cb --- /dev/null +++ b/src/Chemistry/Atoms/Implements/DAtom.cs @@ -0,0 +1,27 @@ +using VirtualChemistry.Constants; + +namespace VirtualChemistry.Chemistry.Atoms.Implements; + +/// +/// atom # 4 +/// +public class DAtom : BaseAtom +{ + /// + public override string Element => ChemistryConstant.Elements.D; + + /// + public override int AtomicNumber => 4; + + /// + public override SU3 AdjointState => ChemistryConstant.AtomAdjoints.AtomAdjointStateD; + + /// + public override SU3 KernelState => ChemistryConstant.AtomKernels.AtomKernelStateD; + + /// + public override SU3 EigenStateMatrix => ChemistryConstant.AtomEigenStates.AtomEigenStateD; + + /// + public override double AtomicMass => ChemistryConstant.AtomicMasses.AtomicMassD; +} \ No newline at end of file diff --git a/src/Chemistry/Atoms/Implements/EAtom.cs b/src/Chemistry/Atoms/Implements/EAtom.cs new file mode 100644 index 0000000..11db095 --- /dev/null +++ b/src/Chemistry/Atoms/Implements/EAtom.cs @@ -0,0 +1,26 @@ +using VirtualChemistry.Constants; + +namespace VirtualChemistry.Chemistry.Atoms.Implements; +/// +/// atom # 10 +/// +public class EAtom : BaseAtom +{ + /// + public override string Element => ChemistryConstant.Elements.E; + + /// + public override int AtomicNumber => 10; + + /// + public override SU3 AdjointState => ChemistryConstant.AtomAdjoints.AtomAdjointStateE; + + /// + public override SU3 KernelState => ChemistryConstant.AtomKernels.AtomKernelStateE; + + /// + public override SU3 EigenStateMatrix => ChemistryConstant.AtomEigenStates.AtomEigenStateE; + + /// + public override double AtomicMass => ChemistryConstant.AtomicMasses.AtomicMassE; +} \ No newline at end of file diff --git a/src/Chemistry/Atoms/Implements/EsAtom.cs b/src/Chemistry/Atoms/Implements/EsAtom.cs new file mode 100644 index 0000000..980ea7c --- /dev/null +++ b/src/Chemistry/Atoms/Implements/EsAtom.cs @@ -0,0 +1,26 @@ +using VirtualChemistry.Constants; + +namespace VirtualChemistry.Chemistry.Atoms.Implements; +/// +/// atom # 3 +/// +public class EsAtom : BaseAtom +{ + /// + public override string Element => ChemistryConstant.Elements.Es; + + /// + public override int AtomicNumber => 3; + + /// + public override SU3 AdjointState => ChemistryConstant.AtomAdjoints.AtomAdjointStateEs; + + /// + public override SU3 KernelState => ChemistryConstant.AtomKernels.AtomKernelStateEs; + + /// + public override SU3 EigenStateMatrix => ChemistryConstant.AtomEigenStates.AtomEigenStateEs; + + /// + public override double AtomicMass => ChemistryConstant.AtomicMasses.AtomicMassEs; +} \ No newline at end of file diff --git a/src/Chemistry/Atoms/Implements/KAtom.cs b/src/Chemistry/Atoms/Implements/KAtom.cs new file mode 100644 index 0000000..a5c41e9 --- /dev/null +++ b/src/Chemistry/Atoms/Implements/KAtom.cs @@ -0,0 +1,27 @@ +using VirtualChemistry.Constants; + +namespace VirtualChemistry.Chemistry.Atoms.Implements; + +/// +/// atom # 7 +/// +public class KAtom : BaseAtom +{ + /// + public override string Element => ChemistryConstant.Elements.K; + + /// + public override int AtomicNumber => 7; + + /// + public override SU3 AdjointState => ChemistryConstant.AtomAdjoints.AtomAdjointStateK; + + /// + public override SU3 KernelState => ChemistryConstant.AtomKernels.AtomKernelStateK; + + /// + public override SU3 EigenStateMatrix => ChemistryConstant.AtomEigenStates.AtomEigenStateK; + + /// + public override double AtomicMass => ChemistryConstant.AtomicMasses.AtomicMassK; +} \ No newline at end of file diff --git a/src/Chemistry/Atoms/Implements/MAtom.cs b/src/Chemistry/Atoms/Implements/MAtom.cs new file mode 100644 index 0000000..2266d74 --- /dev/null +++ b/src/Chemistry/Atoms/Implements/MAtom.cs @@ -0,0 +1,27 @@ +using VirtualChemistry.Constants; + +namespace VirtualChemistry.Chemistry.Atoms.Implements; + +/// +/// atom # 6 +/// +public class MAtom : BaseAtom +{ + /// + public override string Element => ChemistryConstant.Elements.M; + + /// + public override int AtomicNumber => 6; + + /// + public override SU3 AdjointState => ChemistryConstant.AtomAdjoints.AtomAdjointStateM; + + /// + public override SU3 KernelState => ChemistryConstant.AtomKernels.AtomKernelStateM; + + /// + public override SU3 EigenStateMatrix => ChemistryConstant.AtomEigenStates.AtomEigenStateM; + + /// + public override double AtomicMass => ChemistryConstant.AtomicMasses.AtomicMassM; +} \ No newline at end of file diff --git a/src/Chemistry/Atoms/Implements/NullAtom.cs b/src/Chemistry/Atoms/Implements/NullAtom.cs new file mode 100644 index 0000000..38cae78 --- /dev/null +++ b/src/Chemistry/Atoms/Implements/NullAtom.cs @@ -0,0 +1,30 @@ +namespace VirtualChemistry.Chemistry.Atoms.Implements; + +/// +/// atom # -1, the atom used to represent no atom +/// +public sealed class NullAtom : BaseAtom +{ + /// + public NullAtom() + { + } + + /// + public override string Element => "Null"; + + /// + public override int AtomicNumber => throw new Exception(); + + /// + public override double AtomicMass => throw new Exception(); + + /// + public override SU3 KernelState => throw new Exception(); + + /// + public override SU3 AdjointState => throw new Exception(); + + /// + public override SU3 EigenStateMatrix => throw new Exception(); +} \ No newline at end of file diff --git a/src/Chemistry/Atoms/Implements/PAtom.cs b/src/Chemistry/Atoms/Implements/PAtom.cs new file mode 100644 index 0000000..6b8ccaa --- /dev/null +++ b/src/Chemistry/Atoms/Implements/PAtom.cs @@ -0,0 +1,27 @@ +using VirtualChemistry.Constants; + +namespace VirtualChemistry.Chemistry.Atoms.Implements; + +/// +/// atom # 8 +/// +public class PAtom : BaseAtom +{ + /// + public override string Element => ChemistryConstant.Elements.P; + + /// + public override int AtomicNumber => 8; + + /// + public override SU3 AdjointState => ChemistryConstant.AtomAdjoints.AtomAdjointStateP; + + /// + public override SU3 KernelState => ChemistryConstant.AtomKernels.AtomKernelStateP; + + /// + public override SU3 EigenStateMatrix => ChemistryConstant.AtomEigenStates.AtomEigenStateP; + + /// + public override double AtomicMass => ChemistryConstant.AtomicMasses.AtomicMassP; +} diff --git a/src/Chemistry/Atoms/Implements/QAtom.cs b/src/Chemistry/Atoms/Implements/QAtom.cs new file mode 100644 index 0000000..1487b8c --- /dev/null +++ b/src/Chemistry/Atoms/Implements/QAtom.cs @@ -0,0 +1,26 @@ +using VirtualChemistry.Constants; + +namespace VirtualChemistry.Chemistry.Atoms.Implements; +/// +/// atom # 1 +/// +public class QAtom : BaseAtom +{ + /// + public override string Element => ChemistryConstant.Elements.Q; + + /// + public override int AtomicNumber => 1; + + /// + public override SU3 AdjointState => ChemistryConstant.AtomAdjoints.AtomAdjointStateQ; + + /// + public override SU3 KernelState => ChemistryConstant.AtomKernels.AtomKernelStateQ; + + /// + public override SU3 EigenStateMatrix => ChemistryConstant.AtomEigenStates.AtomEigenStateQ; + + /// + public override double AtomicMass => ChemistryConstant.AtomicMasses.AtomicMassQ; +} \ No newline at end of file diff --git a/src/Chemistry/Atoms/Implements/RAtom.cs b/src/Chemistry/Atoms/Implements/RAtom.cs new file mode 100644 index 0000000..b5c297c --- /dev/null +++ b/src/Chemistry/Atoms/Implements/RAtom.cs @@ -0,0 +1,26 @@ +using VirtualChemistry.Constants; + +namespace VirtualChemistry.Chemistry.Atoms.Implements; +/// +/// atom # 2 +/// +public class RAtom : BaseAtom +{ + /// + public override string Element => ChemistryConstant.Elements.R; + + /// + public override int AtomicNumber => 2; + + /// + public override SU3 AdjointState => ChemistryConstant.AtomAdjoints.AtomAdjointStateR; + + /// + public override SU3 KernelState => ChemistryConstant.AtomKernels.AtomKernelStateR; + + /// + public override SU3 EigenStateMatrix => ChemistryConstant.AtomEigenStates.AtomEigenStateR; + + /// + public override double AtomicMass => ChemistryConstant.AtomicMasses.AtomicMassR; +} \ No newline at end of file diff --git a/src/Chemistry/Atoms/Implements/SoAtom.cs b/src/Chemistry/Atoms/Implements/SoAtom.cs new file mode 100644 index 0000000..a8832e5 --- /dev/null +++ b/src/Chemistry/Atoms/Implements/SoAtom.cs @@ -0,0 +1,26 @@ +using VirtualChemistry.Constants; + +namespace VirtualChemistry.Chemistry.Atoms.Implements; +/// +/// atom # 9 +/// +public class SoAtom : BaseAtom +{ + /// + public override string Element => ChemistryConstant.Elements.So; + + /// + public override int AtomicNumber => 9; + + /// + public override SU3 AdjointState => ChemistryConstant.AtomAdjoints.AtomAdjointStateSo; + + /// + public override SU3 KernelState => ChemistryConstant.AtomKernels.AtomKernelStateSo; + + /// + public override SU3 EigenStateMatrix => ChemistryConstant.AtomEigenStates.AtomEigenStateSo; + + /// + public override double AtomicMass => ChemistryConstant.AtomicMasses.AtomicMassSo; +} \ No newline at end of file diff --git a/src/Chemistry/Atoms/Implements/UAtom.cs b/src/Chemistry/Atoms/Implements/UAtom.cs new file mode 100644 index 0000000..025b874 --- /dev/null +++ b/src/Chemistry/Atoms/Implements/UAtom.cs @@ -0,0 +1,27 @@ +using VirtualChemistry.Constants; + +namespace VirtualChemistry.Chemistry.Atoms.Implements; + +/// +/// atom # 11 +/// +public class UAtom : BaseAtom +{ + /// + public override string Element => ChemistryConstant.Elements.U; + + /// + public override int AtomicNumber => 11; + + /// + public override SU3 AdjointState => ChemistryConstant.AtomAdjoints.AtomAdjointStateU; + + /// + public override SU3 KernelState => ChemistryConstant.AtomKernels.AtomKernelStateU; + + /// + public override SU3 EigenStateMatrix => ChemistryConstant.AtomEigenStates.AtomEigenStateU; + + /// + public override double AtomicMass => ChemistryConstant.AtomicMasses.AtomicMassU; +} \ No newline at end of file diff --git a/src/Chemistry/Atoms/Implements/UeAtom.cs b/src/Chemistry/Atoms/Implements/UeAtom.cs new file mode 100644 index 0000000..befe444 --- /dev/null +++ b/src/Chemistry/Atoms/Implements/UeAtom.cs @@ -0,0 +1,27 @@ +using VirtualChemistry.Constants; + +namespace VirtualChemistry.Chemistry.Atoms.Implements; + +/// +/// atom # 5 +/// +public class UeAtom : BaseAtom +{ + /// + public override string Element => ChemistryConstant.Elements.Ue; + + /// + public override int AtomicNumber => 5; + + /// + public override SU3 AdjointState => ChemistryConstant.AtomAdjoints.AtomAdjointStateUe; + + /// + public override SU3 KernelState => ChemistryConstant.AtomKernels.AtomKernelStateUe; + + /// + public override SU3 EigenStateMatrix => ChemistryConstant.AtomEigenStates.AtomEigenStateUe; + + /// + public override double AtomicMass => ChemistryConstant.AtomicMasses.AtomicMassUe; +} \ No newline at end of file diff --git a/src/Chemistry/Atoms/Implements/VAtom.cs b/src/Chemistry/Atoms/Implements/VAtom.cs new file mode 100644 index 0000000..5df86cd --- /dev/null +++ b/src/Chemistry/Atoms/Implements/VAtom.cs @@ -0,0 +1,26 @@ +using VirtualChemistry.Constants; + +namespace VirtualChemistry.Chemistry.Atoms.Implements; +/// +/// atom # 14 +/// +public class VAtom : BaseAtom +{ + /// + public override string Element => ChemistryConstant.Elements.V; + + /// + public override int AtomicNumber => 14; + + /// + public override SU3 AdjointState => ChemistryConstant.AtomAdjoints.AtomAdjointStateV; + + /// + public override SU3 KernelState => ChemistryConstant.AtomKernels.AtomKernelStateV; + + /// + public override SU3 EigenStateMatrix => ChemistryConstant.AtomEigenStates.AtomEigenStateV; + + /// + public override double AtomicMass => ChemistryConstant.AtomicMasses.AtomicMassV; +} \ No newline at end of file diff --git a/src/Chemistry/Atoms/Implements/ViAtom.cs b/src/Chemistry/Atoms/Implements/ViAtom.cs new file mode 100644 index 0000000..8a5b73a --- /dev/null +++ b/src/Chemistry/Atoms/Implements/ViAtom.cs @@ -0,0 +1,27 @@ +using VirtualChemistry.Constants; + +namespace VirtualChemistry.Chemistry.Atoms.Implements; + +/// +/// atom # 15 +/// +public class ViAtom : BaseAtom +{ + /// + public override string Element => ChemistryConstant.Elements.Vi; + + /// + public override int AtomicNumber => 15; + + /// + public override SU3 AdjointState => ChemistryConstant.AtomAdjoints.AtomAdjointStateVi; + + /// + public override SU3 KernelState => ChemistryConstant.AtomKernels.AtomKernelStateVi; + + /// + public override SU3 EigenStateMatrix => ChemistryConstant.AtomEigenStates.AtomEigenStateVi; + + /// + public override double AtomicMass => ChemistryConstant.AtomicMasses.AtomicMassVi; +} diff --git a/src/Chemistry/Atoms/Resolver/AtomResolver.cs b/src/Chemistry/Atoms/Resolver/AtomResolver.cs new file mode 100644 index 0000000..9c78f34 --- /dev/null +++ b/src/Chemistry/Atoms/Resolver/AtomResolver.cs @@ -0,0 +1,72 @@ +using VirtualChemistry.Chemistry.Atoms.Implements; +using VirtualChemistry.Constants; + +namespace VirtualChemistry.Chemistry.Atoms.Resolver; + +/// +/// get an atom instance from descriptions +/// +public static class AtomResolver +{ + /// + /// create an atom from string + /// + /// + /// + /// + public static BaseAtom Resolve(string element) + { + BaseAtom res = element switch + { + ChemistryConstant.Elements.Q => new QAtom(), + ChemistryConstant.Elements.R => new RAtom(), + ChemistryConstant.Elements.Es => new EsAtom(), + ChemistryConstant.Elements.D => new DAtom(), + ChemistryConstant.Elements.Ue => new UeAtom(), + ChemistryConstant.Elements.M => new MAtom(), + ChemistryConstant.Elements.K => new KAtom(), + ChemistryConstant.Elements.P => new PAtom(), + ChemistryConstant.Elements.So => new SoAtom(), + ChemistryConstant.Elements.E => new EAtom(), + ChemistryConstant.Elements.U => new UAtom(), + ChemistryConstant.Elements.A => new AAtom(), + ChemistryConstant.Elements.Cx => new CxAtom(), + ChemistryConstant.Elements.V => new VAtom(), + ChemistryConstant.Elements.Vi => new ViAtom(), + _ => throw new Exception($"TODO-----{element}") + }; + res.C(); + return res; + } + /// + /// create an atom by number + /// + /// + /// + /// + public static BaseAtom Resolve(int atomicNumber) + { + BaseAtom res = atomicNumber switch + { + 1 => new QAtom(), + 2 => new RAtom(), + 3 => new EsAtom(), + 4 => new DAtom(), + 5 => new UeAtom(), + 6 => new MAtom(), + 7 => new KAtom(), + 8 => new PAtom(), + 9 => new SoAtom(), + 10 => new EAtom(), + 11 => new UAtom(), + 12 => new AAtom(), + 13 => new CxAtom(), + 14 => new VAtom(), + 15 => new ViAtom(), + _ => throw new Exception("TODO") + }; + res.C(); + return res; + } + +} \ No newline at end of file diff --git a/src/Chemistry/Bonds/Implements/BaseBond.cs b/src/Chemistry/Bonds/Implements/BaseBond.cs new file mode 100644 index 0000000..ac8a135 --- /dev/null +++ b/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; +/// +/// meta bond +/// +public abstract class BaseBond +{ + + /// + /// + /// + 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(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(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(s); + LieSU3 f1 = f0 ^ LogState.GetFrom(x); + U3 f2 = f1.CayleyPositive(); + return f2.Det().Imaginary; + } + ); + + } + + /// + /// + /// + public Compound Compound => Atom.Compound; + + /// + /// + /// + public HomogeneousMixture HomogeneousMixture => Compound.HomogeneousMixture; + /// + /// null singleton + /// + public static readonly BaseBond Null = new NullBond(); + /// + /// energy hold in this atom + /// + public double Energy { get; set; } + + /// + /// how energy in the molecular is distributed into this atom + /// + public CacheItem EnergyDistributionFactor { get; } + + /// + /// the connected bond + /// + public BaseBond ConnectedBond { get; set; } = Null; + /// + /// atom that holds this bond + /// + public BaseAtom Atom { get; set; } = BaseAtom.Null; + /// + /// type of the bond + /// + public abstract string BondType { get; } + /// + /// kernel state matrix of the matrix, won't change as long as bond type is fixed + /// + public abstract SU3 KernelStateMatrix { get; } + /// + /// adj state matrix of the bond, won't change as long as bond type fixed + /// + public abstract SU3 AdjointStateMatrix { get; } + /// + /// eigen state matrix of bond, won't change as long as bond type fixed + /// + public abstract SU3 EigenStateMatrix { get; } + + /// + /// log of the kernel state matrix + /// never expire + /// + public CacheItem LogKernel { get; } + + /// + /// skew hermitian matrix that represent the information of bond and its atom + /// + public CacheItem HalfBondInformation { get; } + + /// + /// special unitary matrix that represent the structure of the bond(its atom, connected bond and atom of connected bond) + /// + public CacheItem StructureMatrix { get; } + + /// + /// log of connection matrix + /// + public CacheItem LogConnectionMatrix { get; set; } + + /// + /// a special unitary matrix that represent the state of the bond + /// + + + public CacheItem StateMatrix { get; } + + /// + /// log of state matrix + /// + public CacheItem LogState { get; } + + + + + /// + /// connect with another bond + /// + /// + /// + 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; + } + } + /// + /// disconnect with connected bond + /// + 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); + } + /// + /// if the bond is connected + /// + public bool Connected => ConnectedBond != Null; + /// + /// bond number, redundant information of bond type + /// + public abstract int BondNumber { get; } + private void StructureChanged() + { + //Atom.Compound.ResetDistanceCache(); + Atom.Compound.CacheInit(); + HalfBondInformation.Expire(); + StateMatrix.Expire(); + LogConnectionMatrix.Expire(); + + } + + + /// + /// 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 + /// + /// + 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); + + } + + /// + /// 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 + /// + /// + 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); + } + + /// + /// how antibonding energy affected by mixture + /// + public CacheItem AntibondingCatalyzeFactor { get; set; } + + /// + /// hot bonding energy affected by chemical environment + /// + public CacheItem BondingCatalyzeFactor { get; set; } + + public double BondingEnergy => EigenBondingEnergy() + BondingCatalyzeFactor.Get; + public double AntiBondingEnergy => EigenAntibondingEnergy() + AntibondingCatalyzeFactor.Get; + + /// + /// how to pair bonds + /// + 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 + }; + } + } +} diff --git a/src/Chemistry/Bonds/Implements/DBond.cs b/src/Chemistry/Bonds/Implements/DBond.cs new file mode 100644 index 0000000..3030b9a --- /dev/null +++ b/src/Chemistry/Bonds/Implements/DBond.cs @@ -0,0 +1,23 @@ +using VirtualChemistry.Constants; + +namespace VirtualChemistry.Chemistry.Bonds.Implements; +/// +/// bond # 3 +/// +public class DBond : BaseBond +{ + /// + public override string BondType => ChemistryConstant.BondTypes.D; + + /// + public override int BondNumber => 3; + + /// + public override SU3 KernelStateMatrix => ChemistryConstant.BondKernels.BondKernelD; + + /// + public override SU3 AdjointStateMatrix => ChemistryConstant.BondAdjoints.BondAdjointD; + + /// + public override SU3 EigenStateMatrix => ChemistryConstant.BondEigenStates.BondEigenStateD; +} diff --git a/src/Chemistry/Bonds/Implements/HBond.cs b/src/Chemistry/Bonds/Implements/HBond.cs new file mode 100644 index 0000000..a0214ab --- /dev/null +++ b/src/Chemistry/Bonds/Implements/HBond.cs @@ -0,0 +1,24 @@ +using VirtualChemistry.Constants; + +namespace VirtualChemistry.Chemistry.Bonds.Implements; + +/// +/// bond # 7 +/// +public class HBond : BaseBond +{ + /// + public override string BondType => ChemistryConstant.BondTypes.H; + + /// + public override int BondNumber => 7; + + /// + public override SU3 KernelStateMatrix => ChemistryConstant.BondKernels.BondKernelH; + + /// + public override SU3 AdjointStateMatrix => ChemistryConstant.BondAdjoints.BondAdjointH; + + /// + public override SU3 EigenStateMatrix => ChemistryConstant.BondEigenStates.BondEigenStateH; +} \ No newline at end of file diff --git a/src/Chemistry/Bonds/Implements/MBond.cs b/src/Chemistry/Bonds/Implements/MBond.cs new file mode 100644 index 0000000..4eda6b5 --- /dev/null +++ b/src/Chemistry/Bonds/Implements/MBond.cs @@ -0,0 +1,25 @@ +using VirtualChemistry.Constants; + +namespace VirtualChemistry.Chemistry.Bonds.Implements; + +/// +/// bond # 1 +/// +public class MBond : BaseBond +{ + /// + public override string BondType => ChemistryConstant.BondTypes.M; + + /// + public override int BondNumber => 1; + + /// + public override SU3 KernelStateMatrix => ChemistryConstant.BondKernels.BondKernelM; + + /// + public override SU3 AdjointStateMatrix => ChemistryConstant.BondAdjoints.BondAdjointM; + + /// + public override SU3 EigenStateMatrix => ChemistryConstant.BondEigenStates.BondEigenStateM; + +} \ No newline at end of file diff --git a/src/Chemistry/Bonds/Implements/NullBond.cs b/src/Chemistry/Bonds/Implements/NullBond.cs new file mode 100644 index 0000000..f5135ff --- /dev/null +++ b/src/Chemistry/Bonds/Implements/NullBond.cs @@ -0,0 +1,27 @@ +namespace VirtualChemistry.Chemistry.Bonds.Implements; + +/// +/// bond singleton that represent no bond +/// +public class NullBond : BaseBond +{ + /// + public NullBond() + { + } + + /// + public override string BondType => throw new Exception(); + + /// + public override SU3 KernelStateMatrix => throw new Exception(); + + /// + public override SU3 AdjointStateMatrix => throw new Exception(); + + /// + public override SU3 EigenStateMatrix => throw new Exception(); + + /// + public override int BondNumber => throw new Exception(); +} \ No newline at end of file diff --git a/src/Chemistry/Bonds/Implements/PBond.cs b/src/Chemistry/Bonds/Implements/PBond.cs new file mode 100644 index 0000000..42c3fe6 --- /dev/null +++ b/src/Chemistry/Bonds/Implements/PBond.cs @@ -0,0 +1,24 @@ +using VirtualChemistry.Constants; + +namespace VirtualChemistry.Chemistry.Bonds.Implements; + +/// +/// bond # 6 +/// +public class PBond : BaseBond +{ + /// + public override string BondType => ChemistryConstant.BondTypes.P; + + /// + public override int BondNumber => 6; + + /// + public override SU3 KernelStateMatrix => ChemistryConstant.BondKernels.BondKernelP; + + /// + public override SU3 AdjointStateMatrix => ChemistryConstant.BondAdjoints.BondAdjointP; + + /// + public override SU3 EigenStateMatrix => ChemistryConstant.BondEigenStates.BondEigenStateP; +} \ No newline at end of file diff --git a/src/Chemistry/Bonds/Implements/QBond.cs b/src/Chemistry/Bonds/Implements/QBond.cs new file mode 100644 index 0000000..b0d3594 --- /dev/null +++ b/src/Chemistry/Bonds/Implements/QBond.cs @@ -0,0 +1,25 @@ +using VirtualChemistry.Constants; + +namespace VirtualChemistry.Chemistry.Bonds.Implements; + +/// +/// bond # 5 +/// +public class QBond : BaseBond +{ + /// + public override string BondType => ChemistryConstant.BondTypes.Q; + + /// + public override int BondNumber => 5; + + /// + public override SU3 KernelStateMatrix => ChemistryConstant.BondKernels.BondKernelQ; + + /// + public override SU3 AdjointStateMatrix => ChemistryConstant.BondAdjoints.BondAdjointQ; + + /// + public override SU3 EigenStateMatrix => ChemistryConstant.BondEigenStates.BondEigenStateQ; + +} \ No newline at end of file diff --git a/src/Chemistry/Bonds/Implements/SBond.cs b/src/Chemistry/Bonds/Implements/SBond.cs new file mode 100644 index 0000000..5f07595 --- /dev/null +++ b/src/Chemistry/Bonds/Implements/SBond.cs @@ -0,0 +1,24 @@ +using VirtualChemistry.Constants; + +namespace VirtualChemistry.Chemistry.Bonds.Implements; + +/// +/// bond # 2 +/// +public class SBond : BaseBond +{ + /// + public override string BondType => ChemistryConstant.BondTypes.S; + + /// + public override int BondNumber => 2; + + /// + public override SU3 KernelStateMatrix => ChemistryConstant.BondKernels.BondKernelS; + + /// + public override SU3 AdjointStateMatrix => ChemistryConstant.BondAdjoints.BondAdjointS; + + /// + public override SU3 EigenStateMatrix => ChemistryConstant.BondEigenStates.BondEigenStateS; +} \ No newline at end of file diff --git a/src/Chemistry/Bonds/Implements/YBond.cs b/src/Chemistry/Bonds/Implements/YBond.cs new file mode 100644 index 0000000..7ef29d4 --- /dev/null +++ b/src/Chemistry/Bonds/Implements/YBond.cs @@ -0,0 +1,24 @@ +using VirtualChemistry.Constants; + +namespace VirtualChemistry.Chemistry.Bonds.Implements; + +/// +/// bond # 4 +/// +public class YBond : BaseBond +{ + /// + public override string BondType => ChemistryConstant.BondTypes.Y; + + /// + public override int BondNumber => 4; + + /// + public override SU3 KernelStateMatrix => ChemistryConstant.BondKernels.BondKernelY; + + /// + public override SU3 AdjointStateMatrix => ChemistryConstant.BondAdjoints.BondAdjointY; + + /// + public override SU3 EigenStateMatrix => ChemistryConstant.BondEigenStates.BondEigenStateY; +} \ No newline at end of file diff --git a/src/Chemistry/Bonds/Resolver/BondResolver.cs b/src/Chemistry/Bonds/Resolver/BondResolver.cs new file mode 100644 index 0000000..c59fd83 --- /dev/null +++ b/src/Chemistry/Bonds/Resolver/BondResolver.cs @@ -0,0 +1,53 @@ +using VirtualChemistry.Chemistry.Bonds.Implements; +using VirtualChemistry.Constants; + +namespace VirtualChemistry.Chemistry.Bonds.Resolver; + +/// +/// create a bond instance by descriptions +/// +public static class BondResolver +{ + /// + /// create the bond by type + /// + /// + /// + /// + 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; + } + /// + /// create a bond by number + /// + /// + /// + /// + 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; + } +} \ No newline at end of file diff --git a/src/Chemistry/Compounds/Implements/Compound.cs b/src/Chemistry/Compounds/Implements/Compound.cs new file mode 100644 index 0000000..2f8569b --- /dev/null +++ b/src/Chemistry/Compounds/Implements/Compound.cs @@ -0,0 +1,477 @@ +using Skeleton.Algebra.Extensions; +using Skeleton.DataStructure; +using Skeleton.DataStructure.Packs; +using Skeleton.Utils; +using Skeleton.Utils.Helpers; +using Skeleton.Utils.InverseSampling; +using VirtualChemistry.Chemistry.Atoms.Implements; +using VirtualChemistry.Chemistry.Bonds.Implements; +using VirtualChemistry.Chemistry.Compounds.Resolver; +using VirtualChemistry.Chemistry.Mixtures.Implements; +using VirtualChemistry.DataStructure.Packs; + +namespace VirtualChemistry.Chemistry.Compounds.Implements; + +/// +/// a cluster of molecular of the same kind +/// +public class Compound +{ + /// + /// singleton used to represent no compound + /// + public static readonly Compound Null = new Compound(); + private string? IsoRepresentationCache { get; set; } + /// + /// Constructor + /// + public Compound() + { + CacheInit(); + LogStateMatrix = new(x => Atoms + .Select(atom => atom.LogStateMatrix.GetFrom(x)) + .DefaultIfEmpty(LieSU3.Zero) + .SpLieSum() + ); + StateMatrix = new(x => LogStateMatrix.GetFrom(x).Exp()); + CompressionElasticity = new (x => ( + LogStateMatrix.GetFrom(x).CayleyNegative() * StateMatrix.GetFrom(x)) + .ConjugateOn(new DiagR3(0d, 2d * Math.PI, 4d * Math.PI)) + .Rayleigh(LogStateMatrix.GetFrom(x).ColumnAverage) + ); + CompressionBias = new (x => + { + DiagR3 w1 = new (1d / 11d, 7d / 11d, 10d / 11d); + DiagR3 w2 = new (2d / 17d, 8d / 17d, 16d / 17d); + double p1 = StateMatrix.GetFrom(x).ConjugateOn(w1).Rayleigh(StateMatrix.GetFrom(x).Hermitian().ColumnAverage); + double p2 = StateMatrix.GetFrom(x).Hermitian().ConjugateOn(w2).Rayleigh(StateMatrix.GetFrom(x).ColumnAverageBivariant); + double x1 = InverseSampling.StandardNormal(p1) / 4d - 2d * Math.PI; + double x2 = InverseSampling.StandardNormal(p2) / 3d + 2d * Math.PI; + return Math.Atan(Math.Abs(x1) > Math.Abs(x2) ? x1 : x2) * 2d; + } + ); + FreeDensity = new(x => + { + double avg = (Math.Exp(-4d * Math.PI) + Math.Exp(Math.PI)) / 2d; + DiagR3 spectrum = new DiagR3(Math.Exp(-4d * Math.PI), avg, Math.Exp(Math.PI)); + double tempMod = Math.Exp(-Phase); + return tempMod * LogStateMatrix.GetFrom(x) + .CayleyPositive() + .ConjugateOn(spectrum) + .Rayleigh(StateMatrix.GetFrom(x).ColumnAverageBivariant); + }); + } + /// + /// initialize cache + /// + public void CacheInit() + { + HomogeneousMixture.ResetAllCache(); + ResetDistanceCache(); + IsoRepresentationCache = null; + } + /// + /// + /// + public List TraceLog { get; } = new(); + /// + /// compression elasticity, refer to the homogeneous mixture version + /// + public CacheItem CompressionElasticity { get; } + /// + /// compression bias, refer to the homogeneous mixture version + /// + public CacheItem CompressionBias { get; } + /// + /// all atoms in the molecular + /// + public HashSet Atoms { get; set; } = new (); + /// + /// all connections in this compound + /// + public IOrderedEnumerable Connections + { + get + { + HashSet visited = new (); + List res = new (); + foreach (BaseBond bond in Bonds.Where(bond => bond.Connected)) + { + if(visited.Contains(bond)) + continue; + if(res.All(con => !con.TryAbsorbBond(bond))) + res.Add(new FullConnectionInfo(bond)); + visited.Add(bond); + visited.Add(bond.ConnectedBond); + } + return res.OrderBy(con => con.Representation); + } + } + + + public IEnumerable UnconnectedIndex() => Atoms.Count > 15 ? + Array.Empty(): + Bonds.Where(bond => !bond.Connected && bond.Energy > bond.BondingEnergy); + + + /// + /// string used to save the compound + /// + public string ConnectionInfo => String.Join(':', Connections.Select(con => con.Representation)); + /// + /// all bonds in the compound + /// + public IEnumerable Bonds => Atoms.SelectMany(atom => atom.Bonds); + /// + /// the homogeneous mixture of this compound + /// + public HomogeneousMixture HomogeneousMixture { get; set; } = HomogeneousMixture.Null; + /// + /// state matrix of the compound + /// + public CacheItem StateMatrix { get; } + + /// + /// log state matrix + /// + public CacheItem LogStateMatrix { get; } + + /// + /// a quantity to compare phase among different materials + /// + public double Phase => Math.Tanh(Temperature * CompressionElasticity.Get + CompressionBias.Get); + + /// + /// compression stiffness, refer to the homogeneous mixture version + /// + public double CompressionStiffness => Math.Exp(Math.PI / 2d * (5d * Phase - 3d)); + + /// + /// temperature of the compound + /// + public double Temperature => HomogeneousMixture.Temperature; + + /// + /// unique chemical expression of the compound + /// + public string Expression => + Atoms + .GroupBy(atom => atom.Element) + .OrderBy(group => group.First().AtomicNumber) + .Select(group => $"{group.First().Element}"+ (group.Count() > 1 ? group.Count().ToString() : "") ) + .DefaultIfEmpty("") + .Aggregate((a, b) => a + b); + /// + /// dump the compound into save string + /// + /// + /// + public string Dump() + { + Dictionary atomMap = new Dictionary(); + Dictionary bondMap = new Dictionary(); + HashSet records = new HashSet(); + foreach (BaseAtom atom in Atoms) + atomMap[atom] = $"({atom.Element}.{atomMap.Count})"; + foreach (BaseBond bond in Bonds.Where(bond => bond.Connected)) + bondMap[bond] = $"<{bond.BondType}.{bondMap.Count}>"; + + foreach (BaseAtom atom in Atoms) + { + if (atom.Bonds.All(bond => !bond.Connected)) + { + records.Add($"{atomMap[atom]}---(xxx.0)"); + continue; + } + + foreach (BaseBond bond in atom.Bonds.Where(bond => bond.Connected)) + { + if (bond == bond.ConnectedBond) + throw new Exception("TODO"); + if (bondMap[bond] == bondMap[bond.ConnectedBond]) + throw new Exception("TODO"); + string record = + $"{atomMap[atom]}-{bondMap[bond]}-{bondMap[bond.ConnectedBond]}-{atomMap[bond.ConnectedBond.Atom]}"; + string revRecord = + $"{atomMap[bond.ConnectedBond.Atom]}-{bondMap[bond.ConnectedBond]}-{bondMap[bond]}-{atomMap[atom]}"; + if(records.Contains(revRecord)) + continue; + records.Add(record); + + } + } + return String.Join('|', records) + $"{Amount.ExactDoubleString()}"; + } + /// + /// dump with mapped infomation + /// + /// + public CompoundDumpMap MappedDump() + { + Dictionary atomMap = new(); + Dictionary bondMap = new(); + HashSet records = new(); + foreach (BaseAtom atom in Atoms) + atomMap[atom] = $"{atom.Element}.{atomMap.Count}"; + foreach (BaseBond bond in Bonds.Where(bond => bond.Connected)) + bondMap[bond] = $"{bond.BondType}.{bondMap.Count}"; + + foreach (BaseAtom atom in Atoms) + { + if (atom.Bonds.All(bond => !bond.Connected)) + { + records.Add($"({atomMap[atom]})---(xxx.0)"); + continue; + } + + foreach (BaseBond bond in atom.Bonds.Where(bond => bond.Connected)) + { + string record = + $"({atomMap[atom]})-<{bondMap[bond]}>-<{bondMap[bond.ConnectedBond]}>-({atomMap[bond.ConnectedBond.Atom]})"; + string revRecord = + $"({atomMap[bond.ConnectedBond.Atom]})-<{bondMap[bond.ConnectedBond]}>-<{bondMap[bond]}>-({atomMap[atom]})"; + if(records.Contains(revRecord)) + continue; + records.Add(record); + } + } + string dumpString = String.Join('|', records); + return new CompoundDumpMap(atomMap, bondMap, dumpString); + } + + /// + /// make a copy of the compound + /// + /// + public Compound Copy() => CompoundResolver.Resolve(Dump()); + + /// + /// make a copy with map + /// + /// + /// + public PlainPack2 CopyWithBondMap(BaseBond handle) + { + CompoundDumpMap dm = MappedDump(); + CompoundResolveMap rm = CompoundResolver.MappedResolve(dm.DumpString); + if (dm.BondMap.ContainsKey(handle)) + return new PlainPack2(rm.ResolvedCompound, rm.BondMap[dm.BondMap[handle]]); + return new PlainPack2(rm.ResolvedCompound, rm.AtomMap[dm.AtomMap[handle.Atom]].NextUnconnected(handle.BondType)); + } + /// + /// distance between atoms + /// + public Dictionary> DistanceCache { get; set; } = + new Dictionary>(); + /// + /// get the cached distance between two atoms in the compound + /// + /// + /// + /// + public int CachedDistance(BaseAtom a1, BaseAtom a2) + { + if(DistanceCache.Keys.Contains(a1)) + if (DistanceCache[a1].Keys.Contains(a2)) + return DistanceCache[a1][a2]; + if(DistanceCache.Keys.Contains(a2)) + if (DistanceCache[a2].Keys.Contains(a1)) + return DistanceCache[a2][a1]; + return -1; + } + /// + /// save distance of two atoms in this compound into cache + /// + /// + /// + /// + public void CacheDistance(BaseAtom a1, BaseAtom a2, int distance) + { + if (!DistanceCache.ContainsKey(a1)) + DistanceCache[a1] = new Dictionary(); + DistanceCache[a1][a2] = distance; + } + /// + /// reset the cache + /// + public void ResetDistanceCache() => DistanceCache = new Dictionary>(); + /// + /// if two compound has the same structure + /// + /// + /// + public bool IsometricTo(Compound compound) + { + if (!Expression.Equals(compound.Expression)) + return false; + if(!ConnectionInfo.Equals(compound.ConnectionInfo)) + return false; + if (!StateMatrix.Get.IsEqualApprox(compound.StateMatrix.Get)) + return false; + return IsoRepresentation.Equals(compound.IsoRepresentation, StringComparison.OrdinalIgnoreCase); + } + /// + /// amount of molecular in this compound + /// + public double Amount { get; set; } + /// + /// energy per amount + /// + public double UnitEnergy + { + get => Energy / Amount; + set => Energy = value * Amount; + } + /// + /// energy + /// + public double Energy + { + get => Atoms.Sum(atom => atom.Energy) * Amount; + set + { + double energyPerAmount = value / Amount; + double energyFactor = Atoms.Sum(atom => atom.EnergyDistributionFactor.Get); + foreach (BaseAtom atom in Atoms) + atom.Energy = energyPerAmount * atom.EnergyDistributionFactor.Get / energyFactor; + } + } + /// + /// amount : total mixture amount + /// + public double Concentration => + Amount / (HomogeneousMixture == HomogeneousMixture.Null ? Amount : HomogeneousMixture.Amount); + /// + /// split into two compounds + /// + /// + /// + /// + public Compound Split(double amount, bool byRatio = false) + { + if (amount.IsEqualApprox(0)) + return Null; + if ((!byRatio && amount.IsEqualApprox(Amount)) || (byRatio && amount.IsEqualApprox(1))) + { + HomogeneousMixture.Compounds.Remove(this); + HomogeneousMixture = HomogeneousMixture.Null; + return this; + } + double ebs = Energy; + Compound res = Copy(); + if(!byRatio) + { + double ratio = amount / Amount; + Amount -= amount; + res.Amount = amount; + res.Energy = ebs * ratio; + Energy = ebs * (1d - ratio); + } + else + { + res.Amount = Amount * amount; + Amount *= (1d - amount); + res.Energy = ebs * amount; + Energy = ebs * (1d - amount); + } + res.HomogeneousMixture = HomogeneousMixture.Null; + return res; + } + /// + /// split with bond information + /// + /// + /// + /// + /// + public PlainPack2 SplitWithBondMap(double amount, BaseBond handle, bool byRatio = false) + { + /*if (Amount.IsEqualApprox(0)) + throw new Exception("DD"); + if ((Amount * 0.5).IsEqualApprox(0)) + throw new Exception("FS");*/ + PlainPack2 res = CopyWithBondMap(handle); + if (!byRatio) + { + double ratio = amount / Amount; + double splitEnergy = Energy * ratio; + double energy = Energy - splitEnergy; + Amount -= amount; + res.Item1.Amount = amount; + Energy = energy; + res.Item1.Energy = splitEnergy; + } + else + { + double splitEnergy = Energy * amount; + double energy = Energy * (1 - amount); + res.Item1.Amount = Amount * amount; + Amount *= 1 - amount; + Energy = energy; + res.Item1.Energy = splitEnergy; + } + HomogeneousMixture.AddCompound(res.Item1, true, true); + return res; + } + /// + /// mass per molecular + /// + public double MolecularMass => Atoms + .Select(atom => atom.AtomicMass) + .DefaultIfEmpty(0d) + .Sum(); + /// + /// unique representation of a compound + /// + public string IsoRepresentation => IsoRepresentationCache ??= Atoms + .GroupedMinBy(a => -a.AtomicNumber) + .GroupedMinBy(a => a.ConnectedAtoms.Count()) + .GroupedApproxMinBy(a => a.GroupMass) + .GroupedApproxMinBy(a => a.StateMatrix.Get.Det().Real) + .GroupedApproxMinBy(a => a.StateMatrix.Get.Trace().Real) + .GroupedApproxMinBy(a => a.StructureMatrix.Get.Det().Real) + .Select(a => a.FullRepresentation) + .Min() ?? ""; + + + /// + /// a real number from 0 to 1 + /// depend on temperature + /// Resolve Level + /// when A is next to B(above or below) + /// the one with higher resolve level will try to absorb the other one + /// + /// Higher Level => behaviour like solvent + /// Lower Level => behaviour like solute + /// + public double ResolveLevel + { + get + { + double Normalizer(double x) => 2d * (Math.Tanh(-x * 3d) + 1d) / 5d + 1d / 5d; + double C(double x) => 1 - Math.Exp(-x); + double z = C(C(Phase)); + C3 divisor = (StateMatrix.Get.ColumnAverage + LogStateMatrix.Get.ColumnAverageBivariant.Conj()) / 2d; + double factor = StateMatrix.Get.ConjugateOn(new DiagR3(1d, 1d / 9d, 1d / 3d)).Rayleigh(divisor); + return factor * Normalizer(Phase) / Math.Cosh(z); + } + } + + + /* + * Phase of compound is determined by compression stiffness + * higher stiffness => more like gas + * lower stiffness => more like solid + * + * resolvability of one solute A in solvent B(resolve level of A is less than B) + * if both A and B are solid => 0 + * + * otherwise determined by many factors + * + */ + + /// + /// + /// + public CacheItem FreeDensity { get; } +} \ No newline at end of file diff --git a/src/Chemistry/Compounds/Resolver/CompoundResolver.cs b/src/Chemistry/Compounds/Resolver/CompoundResolver.cs new file mode 100644 index 0000000..3e69d0e --- /dev/null +++ b/src/Chemistry/Compounds/Resolver/CompoundResolver.cs @@ -0,0 +1,111 @@ +using System.Text.RegularExpressions; +using Skeleton.Utils; +using VirtualChemistry.Chemistry.Atoms.Implements; +using VirtualChemistry.Chemistry.Atoms.Resolver; +using VirtualChemistry.Chemistry.Bonds.Implements; +using VirtualChemistry.Chemistry.Compounds.Implements; +using VirtualChemistry.DataStructure.Packs; + +namespace VirtualChemistry.Chemistry.Compounds.Resolver; + +/// +/// get compound by save string +/// +public static class CompoundResolver +{ + private const string AtomRegex = @"\(([a-zA-Z]+\.[0-9]+)\)"; + private const string BondRegex = @"<([a-zA-Z]+\.[0-9]+)>"; + + + /// + /// resolve into map + /// + /// + /// + public static CompoundResolveMap MappedResolve(string expression) + { + Compound res = new Compound(); + IEnumerable records = expression.Split('|'); + Dictionary atomMap = new(); + Dictionary bondMap = new(); + foreach (string record in records) + { + List connection = record.Split('-').ToList(); + string atom1 = new Regex(AtomRegex).Match(connection[0]).Result("$1"); + string bond1 = new Regex(BondRegex).Match(connection[1]).Result("$1"); + string bond2 = new Regex(BondRegex).Match(connection[2]).Result("$1"); + string atom2 = new Regex(AtomRegex).Match(connection[3]).Result("$1"); + if (bond1.Equals("xxx.0", StringComparison.OrdinalIgnoreCase)) + { + BaseAtom single = AtomResolver.Resolve(atom1.Split('.')[0]); + atomMap[atom1] = single; + single.Compound = res; + res.Atoms.Add(single); + res.CacheInit(); + return new CompoundResolveMap(atomMap, bondMap, res); + } + + if (!atomMap.Keys.Contains(atom1)) + atomMap[atom1] = AtomResolver.Resolve(atom1.Split('.')[0]); + if (!atomMap.Keys.Contains(atom2)) + atomMap[atom2] = AtomResolver.Resolve(atom2.Split('.')[0]); + atomMap[atom1].Compound = res; + atomMap[atom2].Compound = res; + res.Atoms.Add(atomMap[atom1]); + res.Atoms.Add(atomMap[atom2]); + BaseBond ba = atomMap[atom1].NextUnconnected(bond1.Split('.')[0]); + BaseBond bb = atomMap[atom2].NextUnconnected(bond2.Split('.')[0], ba); + + ba.Connect(bb); + bondMap[bond1] = ba; + bondMap[bond2] = bb; + } + res.CacheInit(); + return new CompoundResolveMap(atomMap, bondMap, res); + } + /// + /// get the compound by save string + /// + /// + /// + public static Compound Resolve(string expression) + { + string[] parts = expression.Split(""); + expression = parts[0]; + Compound res = new Compound(); + res.Amount = parts[1].ByteStringToDouble(); + IEnumerable records = expression.Split('|'); + Dictionary atomMap = new Dictionary(); + foreach (string record in records) + { + List connection = record.Split('-').ToList(); + string atom1 = new Regex(AtomRegex).Match(connection[0]).Result("$1"); + string bond1 = new Regex(BondRegex).Match(connection[1]).Result("$1"); + string bond2 = new Regex(BondRegex).Match(connection[2]).Result("$1"); + string atom2 = new Regex(AtomRegex).Match(connection[3]).Result("$1"); + if (bond1.Equals("xxx.0", StringComparison.OrdinalIgnoreCase)) + { + BaseAtom single = AtomResolver.Resolve(atom1.Split('.')[0]); + single.Compound = res; + res.Atoms.Add(single); + res.CacheInit(); + return res; + } + + if (!atomMap.Keys.Contains(atom1)) + atomMap[atom1] = AtomResolver.Resolve(atom1.Split('.')[0]); + if (!atomMap.Keys.Contains(atom2)) + atomMap[atom2] = AtomResolver.Resolve(atom2.Split('.')[0]); + atomMap[atom1].Compound = res; + atomMap[atom2].Compound = res; + res.Atoms.Add(atomMap[atom1]); + res.Atoms.Add(atomMap[atom2]); + BaseBond b1 = atomMap[atom1].NextUnconnected(bond1.Split('.')[0]); + BaseBond b2 = atomMap[atom2].NextUnconnected(bond2.Split('.')[0], b1); + b1.Connect(b2); + } + res.CacheInit(); + return res; + } + +} \ No newline at end of file diff --git a/src/Chemistry/Containers/IChemicalContainer.cs b/src/Chemistry/Containers/IChemicalContainer.cs new file mode 100644 index 0000000..866607e --- /dev/null +++ b/src/Chemistry/Containers/IChemicalContainer.cs @@ -0,0 +1,27 @@ +using VirtualChemistry.Chemistry.Mixtures.Implements; + +namespace VirtualChemistry.Chemistry.Containers; + +/// +/// +/// +public interface IChemicalContainer +{ + /// + /// + /// + /// + double Volume(); + /// + /// + /// + HeterogeneousMixture Content { get; set; } + /// + /// + /// + public double EnvironmentPressure { get; set; } + /// + /// + /// + public double EnvironmentTemperature { get; set; } +} \ No newline at end of file diff --git a/src/Chemistry/Containers/NullContainer.cs b/src/Chemistry/Containers/NullContainer.cs new file mode 100644 index 0000000..d744eda --- /dev/null +++ b/src/Chemistry/Containers/NullContainer.cs @@ -0,0 +1,29 @@ +using VirtualChemistry.Chemistry.Mixtures.Implements; + +namespace VirtualChemistry.Chemistry.Containers; +/// +/// container singleton used to represent no container +/// +public sealed class NullContainer : IChemicalContainer +{ + private NullContainer() + { + } + + /// + /// instance + /// + public static readonly NullContainer Null = new(); + + /// + public double Volume() => Double.PositiveInfinity; + + /// + public HeterogeneousMixture Content { get; set; } = HeterogeneousMixture.Null; + + /// + public double EnvironmentPressure { get; set; } + + /// + public double EnvironmentTemperature { get; set; } +} \ No newline at end of file diff --git a/src/Chemistry/Environments/Implements/Environment.cs b/src/Chemistry/Environments/Implements/Environment.cs new file mode 100644 index 0000000..4901c06 --- /dev/null +++ b/src/Chemistry/Environments/Implements/Environment.cs @@ -0,0 +1,12 @@ +namespace VirtualChemistry.Chemistry.Environments.Implements; + +/// +/// the environment +/// +public class Environment +{ + /// + /// environment temperature + /// + public double Temperature { get; set; } +} \ No newline at end of file diff --git a/src/Chemistry/Mixtures/Implements/HeterogeneousMixture.cs b/src/Chemistry/Mixtures/Implements/HeterogeneousMixture.cs new file mode 100644 index 0000000..17d6630 --- /dev/null +++ b/src/Chemistry/Mixtures/Implements/HeterogeneousMixture.cs @@ -0,0 +1,158 @@ +using Skeleton.DataStructure.Link; +using VirtualChemistry.Chemistry.Containers; +using VirtualChemistry.Chemistry.Mixtures.Resolver; + +namespace VirtualChemistry.Chemistry.Mixtures.Implements; + +/// +/// a collection of homogeneous mixtures +/// +public class HeterogeneousMixture +{ + public HeterogeneousMixture() + { + } + + public HeterogeneousMixture(IChemicalContainer c) + { + Container = c; + } + + /// + /// singleton used to represent no mixture + /// + public static readonly HeterogeneousMixture Null = new HeterogeneousMixture(); + /// + /// container of the mixture + /// + public IChemicalContainer Container { get; set; } = NullContainer.Null; + /// + /// volume of the container + /// + public double ContainerVolume => Container.Volume(); + /// + /// all homogeneous mixtures within + /// + public HashSet Layers { get; set; } = new(); + /// + /// homogeneous mixtures in order from bottom to top + /// + public Link LayerOrder { get; set; } = new(); + /// + /// add an empty homo mixture to the layers + /// + /// + public HomogeneousMixture AddLayer() + { + HomogeneousMixture res = new HomogeneousMixture(); + AddLayer(res); + return res; + } + + /// + /// remove a homo from the layers + /// + /// + public void RemoveLayer(HomogeneousMixture m) + { + if (this == Null) + return; + Layers.Remove(m); + LayerOrder.Remove(m.Layer); + m.HeterogeneousMixture = Null; + m.Layer = LinkNode.Null; + } + + /// + /// add a homo mixture into this as a new layer + /// + /// + /// + public void AddLayer(HomogeneousMixture m, bool atHead=true) + { + m.HeterogeneousMixture.RemoveLayer(m); + m.HeterogeneousMixture = this; + Layers.Add(m); + if (atHead) + LayerOrder.AddFirst(m); + else + LayerOrder.AddLast(m); + } + /// + /// total amount of molecular + /// + public double Amount + { + get => Layers + .Select(l => l.Amount) + .DefaultIfEmpty(0d) + .Sum(); + set + { + Dictionary res = new Dictionary(); + foreach (HomogeneousMixture l in Layers) + res[l] = value * l.Ratio; + foreach (HomogeneousMixture l in Layers) + l.Amount = res[l]; + } + } + /// + /// if this degenerate into a single homo mixture + /// + public bool IsHomogeneous => Layers.Count == 1; + + /// + /// dump the mixture into save string + /// + /// + public string Dump() => + $"`LAYERS{String.Join(HeterogeneousMixtureResolver.LayerSplit, LayerOrder.Select(layer => layer.Dump()))}`%LAYERS"; + + /// + /// volume of contents + /// + public double Volume => Layers.Sum(layer => layer.Volume); + + /// + /// the total stiffness of all contents + /// + public double EquivalentStiffness => 1d / Layers.Sum(h => 1d / h.CompressionStiffness); + + /// + /// volume of contents under 0 pressure + /// + public double FreeVolume => Layers.Sum(h => h.FreeVolume); + + private double ContentVolume => Layers.Sum(h => h.Volume); + /// + /// the force from container + /// + public double ForceFromContainer => (Math.Min(FreeVolume, Volume) - ContentVolume) * EquivalentStiffness; + + /// + /// process one step of all its homogeneous mixture
+ /// including reaction, resolving/precipitating compounds,
+ /// heat exchange and volume update + /// + public void OneStep() + { + foreach (HomogeneousMixture h in Layers) + { + if(h.Layer.Next.IsEnding) + continue; + if(h.Density > h.Layer.Next.Value.Density) + h.Layer.MoveForward(); + } + foreach (HomogeneousMixture h in Layers) + h.React(); + foreach (HomogeneousMixture h in Layers.ToArray()) + h.PrecipitateAndResolve(); + foreach (HomogeneousMixture h in Layers) + { + h.Degenerate(); + h.VolumeUpdate(); + h.HeatExchange(); + } + } +} + diff --git a/src/Chemistry/Mixtures/Implements/HomogeneousMixture.cs b/src/Chemistry/Mixtures/Implements/HomogeneousMixture.cs new file mode 100644 index 0000000..b6cabc5 --- /dev/null +++ b/src/Chemistry/Mixtures/Implements/HomogeneousMixture.cs @@ -0,0 +1,772 @@ +using System.Numerics; +using Skeleton.Algebra.Extensions; +using Skeleton.Algebra.ScalarFieldStructure; +using Skeleton.DataStructure; +using Skeleton.DataStructure.Link; +using Skeleton.DataStructure.Packs; +using Skeleton.Utils; +using Skeleton.Utils.Helpers; +using Skeleton.Utils.InverseSampling; +using VirtualChemistry.Chemistry.Bonds.Implements; +using VirtualChemistry.Chemistry.Compounds.Implements; +using VirtualChemistry.Chemistry.Mixtures.Resolver; +using VirtualChemistry.Constants; +using VirtualChemistry.DataStructure.Packs; + +namespace VirtualChemistry.Chemistry.Mixtures.Implements; +/// +/// mixture of many well mixed compounds that can not be seperated by psy methods +/// +public class HomogeneousMixture +{ + /// + /// Constructor by default + /// + public HomogeneousMixture() + { + LogStateMatrix = new CacheItem(x => + Compounds.Select(compound => compound.LogStateMatrix.GetFrom(x) * compound.Concentration) + .DefaultIfEmpty(LieSU3.Zero) + .SpLieSum() + ); + StateMatrix = new(x => LogStateMatrix.GetFrom(x).Exp()); + CompressionElasticity = new(x => + { + DiagR3 spectrum = new (0d, 2d * Math.PI, 4d * Math.PI); + return (LogStateMatrix.GetFrom(x).CayleyNegative() * StateMatrix.GetFrom(x)) + .ConjugateOn(spectrum) + .Rayleigh(LogStateMatrix.GetFrom(x).ColumnAverage); + } + ); + CompressionBias = new(x => + { + DiagR3 w1 = new (1d / 11d, 7d / 11d, 10d / 11d); + DiagR3 w2 = new (2d / 17d, 8d / 17d, 16d / 17d); + double p1 = StateMatrix.GetFrom(x) + .ConjugateOn(w1) + .Rayleigh(StateMatrix.GetFrom(x).Hermitian().ColumnAverage); + double p2 = StateMatrix.GetFrom(x) + .Hermitian() + .ConjugateOn(w2) + .Rayleigh(StateMatrix.GetFrom(x).ColumnAverageBivariant); + double x1 = InverseSampling.StandardNormal(p1) / 4d - 2d * Math.PI; + double x2 = InverseSampling.StandardNormal(p2) / 3d + 2d * Math.PI; + return Math.Atan(Math.Abs(x1) > Math.Abs(x2) ? x1 : x2) * 2d; + } + ); + FreeDensity = new(x => + { + double avg = (Math.Exp(-4d * Math.PI) + Math.Exp(Math.PI)) / 2d; + DiagR3 spectrum = new DiagR3(Math.Exp(-4d * Math.PI), avg, Math.Exp(Math.PI)); + double tempMod = Math.Exp(-Phase); + return tempMod * LogStateMatrix.GetFrom(x) + .CayleyPositive() + .ConjugateOn(spectrum) + .Rayleigh(StateMatrix.GetFrom(x).ColumnAverageBivariant); + } + ); + HeatConductivity = new(x => + StateMatrix.GetFrom(x).ConjugateOn(ChemistryConstant.EnergyConductivitySpectrum) + .Rayleigh(LogStateMatrix.GetFrom(x) * StateMatrix.GetFrom(x).ColumnAverageBivariant) + ); + ResetAllCache(); + } + + /// + /// singleton used to represent no mixture + /// + public static readonly HomogeneousMixture Null = new(); + + private IEnumerable ConnectedIndex() + { + HashSet bs = new(); + foreach (Compound c in Compounds) + { + if(c.Amount.IsEqualApprox(0)) + continue; + + + foreach (BaseBond b in c.Bonds) + { + if (b.Compound.Amount.IsEqualApprox(0)) + throw new Exception("P??"); + + if (b.Compound != c) + throw new Exception("??!!?"); + if ( + b.Connected && + b.Energy > b.AntiBondingEnergy && + b.ConnectedBond.Energy > b.ConnectedBond.AntiBondingEnergy && + !bs.Contains(b.ConnectedBond) + ) + { + bs.Add(b); + } + } + } + + return bs; + //Compounds.Where(c => !(0.5 * c.Amount).IsEqualApprox(0)).SelectMany(x => x.ConnectedIndex()); + } + public void RPC() + { + + int s1 = 0; + int s2 = 0; + int s3 = 0; + int s4 = 0; + foreach (Compound c in Compounds) + { + foreach (BaseBond b in c.Bonds) + { + switch (b.BondingGroup) + { + case 0: + s1++; + break; + case 1: + s2++; + break; + case 2: + s3++; + break; + default: + s4++; + break; + } + } + } + Console.WriteLine($"{s1} ---- {s2} ---- {s3} --- {s4}"); + } + + private BondGroups ReactionGroup() + { + HashSet group1 = new(); + HashSet group2 = new(); + HashSet group3 = new(); + HashSet group4 = new(); + foreach (Compound c in Compounds) + { + if ((c.Amount * 0.125).IsEqualApprox(0) || c.Atoms.Count > 15) + continue; + if(c.Expression.Equals("P2So")) + Console.WriteLine("?"); + foreach(BaseBond b in c.Bonds) + { + if(b.Connected || !(b.Energy > b.BondingEnergy)) + continue; + switch (b.BondingGroup) + { + case 0: + group1.Add(b); + break; + case 1: + group2.Add(b); + break; + case 2: + group3.Add(b); + break; + case 3: + group4.Add(b); + break; + } + } + } + return new(group1, group2, group3, group4); + } + + /// + /// elasticity range: 0, 4pi + /// elasticity controls the speed k drops from e^pi to e^-4pi + /// elasticity -> 0; k -> constant + /// b -> 4pi; k -> e^pi + /// b -> -4pi; k -> e^-4pi + /// + public CacheItem CompressionElasticity { get; } + + /// + /// bias range: -4pi, 4pi + /// bias controls the temp of triple point + /// higher b -> higher melting point and higher boiling point + /// + public CacheItem CompressionBias { get; } + + /// + /// phase -> 1 implies material is more like gas + /// phase -> -1 implies material is more like solid + /// + /// + public double Phase => Compounds.Sum(c => c.Concentration * c.Phase); + + /// + /// stiffness range: e^-4pi, e^pi + /// higher stiffness -> harder to compress -> solid like + /// lower stiffness -> easier to compress -> gas like + /// Stiffness K, Density D, Free Density D0, Temperature T + /// P = K(1/D - 1/D0) + /// + /// Temperature range: -pi, pi + /// + public double CompressionStiffness => Math.Exp(Math.PI / 2d * (5d * Phase - 3d)); + + /// + /// e^-4pi ... e^pi + /// + public CacheItem FreeDensity { get; } + + /// + /// volume that under 0 compression + /// + public double FreeVolume => Amount / FreeDensity.Get; + + /// + /// P = K (1/D - 1/D0) + /// + public double Pressure => (Volume - FreeVolume) * CompressionStiffness; + + /// + /// Volume of the mixture + /// + public double Volume { get; set; } + + /// + /// Heat transfer ability + /// + public CacheItem HeatConductivity { get; } + + /// + /// energy / amount = temperature + /// + public double Temperature + { + get => Energy / Amount; + set => Energy = value * Amount; + } + /// + /// total Energy the mixture holds + /// + public double Energy + { + get => Compounds.Sum(c => c.Energy); + set + { + foreach (Compound c in Compounds) + c.Energy = value * c.Concentration; + } + } + /// + /// reset all cache + /// + public void ResetAllCache() + { + LayerCache = LinkNode.Null; + } + /// + /// the heterogeneous mixture that holds this homogeneous mixture + /// + public HeterogeneousMixture HeterogeneousMixture { get; set; } = HeterogeneousMixture.Null; + + /// + /// Compounds of this mixture + /// + public HashSet Compounds { get; set; } = new HashSet(); + /// + /// Compounds whose ratio is greater than 1% + /// + public IEnumerable MainCompounds => Compounds.Where(c => c.Concentration > 0.001d); + private LinkNode LayerCache { get; set; } = LinkNode.Null; + /// + /// + /// + public LinkNode Layer + { + get => LayerCache == LinkNode.Null + ? HeterogeneousMixture.LayerOrder.Find(this) + : LayerCache; + set => LayerCache = value; + } + + /// + /// State matrix of this mixture + /// + public CacheItem StateMatrix { get; } + + /// + /// Log of the state matrix + /// + public CacheItem LogStateMatrix { get; } + + /// + /// Split the homogeneous mixture into another heterogeneous mixture by ratio or amount; + /// + /// The amount or ratio + /// another heterogeneous mixture + /// determine if the amount variable is ratio or amount + /// + public HomogeneousMixture Split(double amount, HeterogeneousMixture into, bool byRatio = false) + { + double ebs = Energy; + double vbs = Volume; + double ratio = byRatio ? amount : amount / Amount; + HomogeneousMixture res = new HomogeneousMixture(); + foreach (Compound c in Compounds) + { + Compound split = c.Split(amount, byRatio); + if(split != Compound.Null) + res.AddCompound(split); + } + res.Energy = ebs*ratio; + res.Volume = vbs*ratio; + Energy = ebs * (1 - ratio); + Volume = vbs * (1 - ratio); + into.AddLayer(res); + return res; + } + + /// + /// Add compound into the homogeneous mixture + /// + /// compound to add + /// + /// add without combine to isomorphic compound + public void AddCompound(Compound compound, bool resetCache = true, bool skipCombine=false) + { + if (compound.HomogeneousMixture == this) + return; + bool added = false; + if(!skipCombine) + { + foreach (Compound sCompound in Compounds.Where(x => x.Expression == compound.Expression)) + { + if (sCompound.IsometricTo(compound)) + { + sCompound.Amount += compound.Amount; + compound.Amount = 0d; + sCompound.Energy += compound.Energy; + added = true; + break; + } + } + } + + compound.HomogeneousMixture = added ? Null : this; + if (!added) + Compounds.Add(compound); + if(resetCache) + ResetAllCache(); + } + /// + /// Amount of this mixture + /// + public double Amount + { + get => Compounds + .Select(c => c.Amount) + .DefaultIfEmpty(0) + .Sum(); + set + { + Dictionary res = new Dictionary(); + foreach (Compound c in Compounds) + res[c] = value * c.Concentration; + foreach (Compound c in Compounds) + c.Amount = res[c]; + } + } + + /// + /// ratio of this mixture in the heterogeneous mixture + /// + public double Ratio => + Amount / (HeterogeneousMixture == HeterogeneousMixture.Null ? Amount : HeterogeneousMixture.Amount); + + + /// + /// + /// + /// + public void Resolve(HomogeneousMixture homo) + { + foreach (Compound c in homo.Compounds) + Resolve(c, false); + ResetAllCache(); + + } + + /// + /// + /// + /// + /// + public void Resolve(Compound compound, bool resetCache = true) + { + AddCompound(compound, resetCache); + } + + /// + /// density of the mixture, amount of molecular per volume + /// + public double Density => Amount / Volume; + + /// + /// save mixture into string + /// + /// + public string Dump() => + $"`ENV?ENG{Energy.ExactDoubleString()}?V{Volume.ExactDoubleString()}`%ENV`COMPS{String.Join(HomogeneousMixtureResolver.CompoundsSplit, Compounds.Select(c => c.Dump()))}`%COMPS"; + + /// + /// update volume + /// + public void VolumeUpdate() + { + double sVolume = Volume + HeterogeneousMixture.ForceFromContainer / CompressionStiffness; + Volume = Math.Max(sVolume, 1E-6); + } + + /// + /// exchange energy to neighbor mixtures + /// + public void HeatExchange() + { + if (Layer != Layer.Parent.First) + { + HomogeneousMixture mPrevious = Layer.Previous.Value; + double balanceTemp = (Energy + mPrevious.Energy) / (Amount + mPrevious.Amount); + double dTemp = balanceTemp - Temperature; + double conductivity = Math.Sqrt(HeatConductivity.Get * mPrevious.HeatConductivity.Get); + double transferTemp = dTemp * conductivity; + Temperature += transferTemp; + mPrevious.Temperature -= transferTemp; + } + if (Layer != Layer.Parent.Last) + { + HomogeneousMixture mNext = Layer.Next.Value; + double balanceTemp = (Energy + mNext.Energy) / (Amount + mNext.Amount); + double dTemp = balanceTemp - Temperature; + double conductivity = Math.Sqrt(HeatConductivity.Get * mNext.HeatConductivity.Get); + double transferTemp = dTemp * conductivity; + Temperature += transferTemp; + mNext.Temperature -= transferTemp; + } + + double envBalanceTemp = HeterogeneousMixture.Container.EnvironmentTemperature; + double envDTemp = envBalanceTemp - Temperature; + double envConductivity = HeatConductivity.Get; + Temperature += envDTemp * envConductivity; + } + private void ReactConnecting(HashSet g1, HashSet g2) + { + double maxCd = -1; + BaseBond xb1 = BaseBond.Null; + BaseBond xb2 = BaseBond.Null; + foreach (BaseBond b1 in g1) + { + foreach (BaseBond b2 in g2) + { + H3 cd = new(Complex.ImaginaryOne / 2 * (b1.LogState.Get + b2.LogState.Get)); + double cdX = cd.Rayleigh((b1.LogState.Get + b2.LogState.Get).ColumnAverage); + if (cdX > maxCd) + { + maxCd = cdX; + xb1 = b1; + xb2 = b2; + } + } + } + if (xb1 == BaseBond.Null) + return; + if (xb1 == xb2) + { + PlainPack2 s = xb1.Compound.SplitWithBondMap(0.5, xb2, true); + xb2 = s.Item2; + } + + if (xb1.Compound.Amount.IsSignificantlyGreaterThen(xb2.Compound.Amount)) + { + double amountDiff = xb1.Compound.Amount - xb2.Compound.Amount; + Compound oth = xb1.Compound.Split(amountDiff); + xb1.HomogeneousMixture.AddCompound(oth); + } + + if (xb1.Compound.Amount.IsSignificantlyLessThen(xb2.Compound.Amount)) + { + double amountDiff = xb2.Compound.Amount - xb1.Compound.Amount; + Compound oth = xb2.Compound.Split(amountDiff); + xb2.HomogeneousMixture.AddCompound(oth); + } + + string log = $"From Composite React of {xb1.Compound.Amount}*{xb1.Compound.Expression} and {xb2.Compound.Amount}*{xb2.Compound.Expression}"; + + xb1.Connect(xb2); + xb1.Compound.TraceLog.Add(log); + } + + private int ReactionType() + { + Complex d = LogStateMatrix.Get.CayleyPositive().Det(); + return d.Phase 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 + }; + } + + /// + /// react within the homogeneous mixture + /// + public void React() + { + BaseBond[] connected = ConnectedIndex().ToArray(); + foreach (BaseBond bond in connected) + { + if (bond.Compound.Amount.IsEqualApprox(0)) + continue; + if ( + bond.Energy > bond.AntiBondingEnergy && + bond.ConnectedBond.Energy > bond.ConnectedBond.AntiBondingEnergy + ) + { + bond.Disconnect(); + } + + } + for (int i = 1; i <= 3; i++) + { + BondGroups bg = ReactionGroup(); + int type = ReactionType(); + if(type == 0) + { + ReactConnecting(bg.Group1, bg.Group3); + ReactConnecting(bg.Group2, bg.Group4); + } + else if (type == 1) + { + ReactConnecting(bg.Group1, bg.Group1); + ReactConnecting(bg.Group3, bg.Group3); + } + else if (type == 2) + { + ReactConnecting(bg.Group2, bg.Group2); + ReactConnecting(bg.Group3, bg.Group3); + } + else + { + ReactConnecting(bg.Group1, bg.Group4); + ReactConnecting(bg.Group2, bg.Group3); + } + } + + } + /// + /// higher level mixture can resolve lower level mixture + /// + public double ResolveLevel => Compounds.Sum(c => c.Concentration * c.ResolveLevel); + /// + /// ability to resolve other compound + /// + /// + /// + public double Resolvability(Compound c) + { + LieSU3 s = c.LogStateMatrix.Get ^ LogStateMatrix.Get; + C3 a = new(1, 0.5, 1); + C3 b = new(0.5, 1, 0.5); + Complex w = a * s * b; + w = ComplexFieldStructure.Structure.Fix(w); + double res = -Math.Tan((Math.Sin(w.Phase) * Math.PI - Math.PI) / 2); + if (Math.Abs(res) < 1E-6) + return 0; + return res; + } + /// + /// a tensor of order 3 abstracted from this mixture + /// + public C3 Ita => Compounds + .Select(c => c.Concentration * c.StateMatrix.Get.ColumnAverageBivariant) + .Aggregate((a, b) => a + b); + /// + /// Red component + /// + public Byte ColorRed => (Byte) + ( + StateMatrix.Get + .ConjugateOn(ChemistryConstant.OpticalFilter.RedDefinite) + .Rayleigh(Ita) * 255 + ); + /// + /// Green component + /// + public Byte ColorGreen => (Byte) + ( + StateMatrix.Get + .ConjugateOn(ChemistryConstant.OpticalFilter.GreenDefinite) + .Rayleigh(Ita) * 255 + ); + /// + /// Blue component + /// + public Byte ColorBlue => (Byte) + ( + StateMatrix.Get + .ConjugateOn(ChemistryConstant.OpticalFilter.BlueDefinite) + .Rayleigh(Ita) * 255 + ); + /// + /// Transparent component + /// + public Byte ColorTransparent => (Byte) + ( + StateMatrix.Get + .ConjugateOn(ChemistryConstant.OpticalFilter.OpacityDefinite) + .Rayleigh(Ita) * 255 + ); + + /// + /// Should not be used by any methods except ResolveNext and ResolvePrev + /// + private void PureResolve(HomogeneousMixture hOther) + { + HashSet toAbsorb = new (); + foreach (Compound c in hOther.Compounds) + { + double maxRes = Resolvability(c); + double oRes = hOther.Resolvability(c); + if(oRes < 0 || maxRes < 0 || maxRes < oRes) + continue; + double maxResAmount = Amount; + foreach (Compound ct in Compounds) + if (ct.IsometricTo(c)) + maxResAmount -= ct.Amount; + maxResAmount *= maxRes; + if(maxResAmount < 0) + continue; + double splitAmount = Math.Min(maxResAmount, c.Amount); + if(splitAmount.IsEqualApprox(0)) + continue; + string log = $"{splitAmount} split from {c.Amount} * {c.Expression}"; + Compound cSplit = c.Split(splitAmount); + cSplit.TraceLog.Add(log); + toAbsorb.Add(cSplit); + } + foreach (Compound c in toAbsorb) + AddCompound(c); + if(hOther.Amount.IsEqualApprox(0)) + HeterogeneousMixture.RemoveLayer(hOther); + } + private void ResolveNext() + { + if (HeterogeneousMixture == HeterogeneousMixture.Null || Layer.Next.IsTail) + return; + PureResolve(Layer.Next.Value); + } + private void ResolvePrevious() + { + if (HeterogeneousMixture == HeterogeneousMixture.Null || Layer.Previous.IsHead) + return; + PureResolve(Layer.Previous.Value); + } + + private bool Resolve() + { + ResolveNext(); + ResolvePrevious(); + return false; + } + + private bool Precipitate() + { + foreach (Compound c in Compounds) + { + double maxRes = Resolvability(c); + if (maxRes < 0) + continue; + double maxResAmount = maxRes * (Amount - c.Amount); + if(maxResAmount.IsEqualApprox(0)) + continue; + + double aDiff = c.Amount - maxResAmount; + if (aDiff.IsSignificantlyGreaterThen(0) ) + { + Compound oth = c.Split(c.Amount - maxResAmount); + if (oth.FreeDensity.Get > FreeDensity.Get) + { + if (Layer.Next.IsEnding) + { + HomogeneousMixture m = HeterogeneousMixture.AddLayer(); + m.Layer.MoveAfter(Layer); + m.AddCompound(oth); + m.Volume = m.FreeVolume; + return true; + } + Layer.Next.Value.AddCompound(oth); + } + else + { + if (Layer.Previous.IsEnding) + { + HomogeneousMixture m = HeterogeneousMixture.AddLayer(); + m.Layer.MoveBefore(Layer); + m.AddCompound(oth); + m.Volume = m.FreeVolume; + return true; + } + Layer.Previous.Value.AddCompound(oth); + } + } + } + return false; + } + + /// + /// + /// + /// if any new layer created or deleted + public bool PrecipitateAndResolve() + { + bool a = Resolve(); + bool b = Precipitate(); + return a || b; + } + + /// + /// combine identical compounds in the miture + /// + public void Degenerate() + { + HashSet toRemove = new(); + foreach (IGrouping g in Compounds.GroupBy(c => c.Expression)) + { + HashSet classified = new(); + foreach (Compound c in g) + { + classified.Add(c); + foreach (Compound cx in g.Where(x => !classified.Contains(x))) + { + if (!c.IsometricTo(cx)) + continue; + classified.Add(cx); + c.Amount += cx.Amount; + cx.Amount = 0; + toRemove.Add(cx); + } + } + } + + foreach (Compound c in toRemove) + { + c.HomogeneousMixture = Null; + Compounds.Remove(c); + } + } + + public static void Test() + { + Func a = x => -Math.Tan((Math.PI * x - Math.PI) / 2); + Console.WriteLine(a(0)); + Console.WriteLine(a(-1)); + Console.WriteLine(a(1)); + } +} diff --git a/src/Chemistry/Mixtures/Resolver/HeterogeneousMixtureResolver.cs b/src/Chemistry/Mixtures/Resolver/HeterogeneousMixtureResolver.cs new file mode 100644 index 0000000..70aced4 --- /dev/null +++ b/src/Chemistry/Mixtures/Resolver/HeterogeneousMixtureResolver.cs @@ -0,0 +1,32 @@ +using System.Text.RegularExpressions; +using VirtualChemistry.Chemistry.Mixtures.Implements; + +namespace VirtualChemistry.Chemistry.Mixtures.Resolver; + +/// +/// get heterogeneous mixture from descriptions +/// +public static class HeterogeneousMixtureResolver +{ + private const string LayersRegex = @"`LAYERS(.*)`%LAYERS"; + private static readonly Regex Matcher = new Regex(LayersRegex); + /// + /// split different homogeneous mixtures + /// + public const string LayerSplit = ""; + + /// + /// + /// + /// + /// + public static HeterogeneousMixture Resolve(string dumpString) + { + HeterogeneousMixture res = new HeterogeneousMixture(); + string layers = Matcher.Match(dumpString).Result("$1"); + foreach (string layerString in layers.Split(LayerSplit)) + res.AddLayer(HomogeneousMixtureResolver.Resolve(layerString), false); + return res; + } + +} diff --git a/src/Chemistry/Mixtures/Resolver/HomogeneousMixtureResolver.cs b/src/Chemistry/Mixtures/Resolver/HomogeneousMixtureResolver.cs new file mode 100644 index 0000000..298333c --- /dev/null +++ b/src/Chemistry/Mixtures/Resolver/HomogeneousMixtureResolver.cs @@ -0,0 +1,45 @@ +using System.Text.RegularExpressions; +using Skeleton.Utils; +using VirtualChemistry.Chemistry.Compounds.Resolver; +using VirtualChemistry.Chemistry.Mixtures.Implements; + +namespace VirtualChemistry.Chemistry.Mixtures.Resolver; + +/// +/// get homogeneous mixture from description +/// +public static class HomogeneousMixtureResolver +{ + private const string EnvironmentRegex = @"`ENV(.*)`%ENV"; + private const string CompoundsRegex = @"`COMPS(.*)`%COMPS"; + private const string MatchString = EnvironmentRegex + CompoundsRegex; + private static readonly Regex Matcher = new Regex(MatchString); + + private const string VolumeRegex = @"\?V\(([0-9,]+)\)"; + private const string EnergyRegex = @"\?E\(([0-9,]+)\)"; + private const string EnvMatchString = VolumeRegex + EnergyRegex; + private static readonly Regex EnvMatcher = new Regex(EnvMatchString); + /// + /// split different compounds + /// + public const string CompoundsSplit = ""; + /// + /// + /// + /// + /// + public static HomogeneousMixture Resolve(string dumpString) + { + Match result = Matcher.Match(dumpString); + HomogeneousMixture res = new HomogeneousMixture(); + string envString = result.Result("$1"); + Match envRes = EnvMatcher.Match(envString); + double volume = envRes.Result("$1").ByteStringToDouble(); + double energy = envRes.Result("$2").ByteStringToDouble(); + foreach (string compoundString in result.Result("$2").Split(CompoundsSplit)) + res.AddCompound(CompoundResolver.Resolve(compoundString)); + res.Volume = volume; + res.Energy = energy; + return res; + } +} \ No newline at end of file diff --git a/src/Constants/CachedKeys.cs b/src/Constants/CachedKeys.cs new file mode 100644 index 0000000..6e25998 --- /dev/null +++ b/src/Constants/CachedKeys.cs @@ -0,0 +1,103 @@ +namespace VirtualChemistry.Constants; + +/// +/// keys for tensors +/// +public static class CachedKeys +{ + /// + /// + /// + public const string StateMatrix = $"{TypePrefix.SU3}_StateMatrix"; + /// + /// + /// + public const string LogStateMatrix = $"{TypePrefix.LieSU3}_LogStateMatrix"; + /// + /// + /// + public const string StructureMatrix = $"{TypePrefix.SU3}_StructureMatrix"; + /// + /// + /// + public const string LogConnectionMatrix = $"{TypePrefix.LieSU3}_LogConnectionMatrix"; + /// + /// + /// + public const string EnergyDistributionFactor = $"{TypePrefix.Complex}_EnergyDistributionFactor"; + /// + /// + /// + public const string HeatConductivity = $"{TypePrefix.Double}_HeatConductivity"; + /// + /// + /// + public const string FreeDensity = $"{TypePrefix.Double}_FreeDensity"; + /// + /// + /// + public const string CompressionStiffness = $"{TypePrefix.Double}_CompressionStiffness"; + /// + /// + /// + public const string CompressionElasticity = $"{TypePrefix.Double}_CompressionElasticity"; + /// + /// + /// + public const string CompressionBias = $"{TypePrefix.Double}_CompressionBias"; + /// + /// + /// + public const string BondInformationMatrix = $"{TypePrefix.SU3}_BondInformationMatrix"; + /// + /// + /// + public const string LogEigenStateMatrix = $"{TypePrefix.LieSU3}_LogEigenStateMatrix"; + /// + /// + /// + public const string LogKernelMatrix = $"{TypePrefix.LieSU3}_LogKernelMatrix"; + /// + /// + /// + public const string HalfBondInformation = $"{TypePrefix.LieSU3}_HalfBondInformation"; +} +/// +/// tensor type info +/// +public static class TypePrefix +{ + /// + /// real + /// + public const string Double = "D"; + /// + /// complex + /// + public const string Complex = "C"; + /// + /// c3 + /// + public const string C3 = "C3"; + /// + /// c33 + /// + public const string C33 = "C33"; + /// + /// u3 + /// + public const string U3 = "U3"; + /// + /// su3 + /// + public const string SU3 = "SU3"; + /// + /// lie u3 + /// + public const string LieU3 = "LieU3"; + /// + /// lie su3 + /// + public const string LieSU3 = "LieSU3"; + +} diff --git a/src/Constants/ChemistryConstant.cs b/src/Constants/ChemistryConstant.cs new file mode 100644 index 0000000..c584bcb --- /dev/null +++ b/src/Constants/ChemistryConstant.cs @@ -0,0 +1,653 @@ +using System.Numerics; +using Skeleton.Constants; +using VirtualChemistry.Utils.Helpers; + +namespace VirtualChemistry.Constants; + +/// +/// chemistry related constants +/// +public static class ChemistryConstant +{ + private static readonly Complex I = Complex.ImaginaryOne; + /// + /// volume of a single atom + /// + public static readonly double AtomVolume = Math.Exp(-2 * Math.PI); + + /// + /// to calculate edv for atom + /// + public static readonly C3 AtomEnergyDistributionVector = new C3 + ( + Complex.Exp(Math.PI * I / 4d), + Complex.Exp(Math.PI * I * 11d / 12d), + Complex.Exp(Math.PI * I * 19d / 12d) + ); + /// + /// to calculate edv for bond + /// + public static readonly C3 BondEnergyDistributionVector = new C3 + ( + Complex.Exp(Math.PI * I * 5d / 4d), + Complex.Exp(Math.PI * 23d / 12d), + Complex.Exp(Math.PI * I * 7d / 12d) + ); + + /// + /// to calculate edv from atom + /// + public static readonly DiagR3 AtomEnergyDistributionSpectrum = new DiagR3(1d / 4d, 11d / 12d, 19d / 12d); + /// + /// to calculate edv for bond + /// + public static readonly DiagR3 BondEnergyDistributionSpectrum = new DiagR3(5d / 4d, 23d / 12d, 7d / 12d); + /// + /// to calculate energy conductivity + /// + public static readonly DiagR3 EnergyConductivitySpectrum = new DiagR3(2d/15d, 7d/15d, 11d/15d); + /// + /// all bond types + /// + public static class BondTypes + { + /// + /// type m + /// + public const string M = "m"; + /// + /// type s + /// + public const string S = "s"; + /// + /// type d + /// + public const string D = "d"; + /// + /// type y + /// + public const string Y = "y"; + /// + /// type q + /// + public const string Q = "q"; + /// + /// type p + /// + public const string P = "p"; + /// + /// type h + /// + public const string H = "h"; + /// + /// all types + /// + public static readonly string[] All = new[] { M, S, D, Y, Q, P, H }; + } + + /// + /// all atom types + /// + public static class Elements + { + /// + /// atom q + /// + public const string Q = "Q"; + /// + /// atom r + /// + public const string R = "R"; + /// + /// atom es + /// + public const string Es = "Es"; + /// + /// atom d + /// + public const string D = "D"; + /// + /// atom ue + /// + public const string Ue = "Ue"; + /// + /// atom m + /// + public const string M = "M"; + /// + /// atom k + /// + public const string K = "K"; + /// + /// atom p + /// + public const string P = "P"; + /// + /// atom so + /// + public const string So = "So"; + /// + /// atom e + /// + public const string E = "E"; + /// + /// atom u + /// + public const string U = "U"; + /// + /// atom a + /// + public const string A = "A"; + /// + /// atom cx + /// + public const string Cx = "Cx"; + /// + /// atom v + /// + public const string V = "V"; + /// + /// atom vi + /// + public const string Vi = "Vi"; + } + + private static double AtomicMass(int atomicNumber) => + Math.Tan(ChemistryHelper.AtomHeightOrder(atomicNumber) * Math.PI / 2d); + + private static SU3 AtomEigenState(int atomicNumber) + { + SU3 dam = ChemistryHelper.AtomAdjointStateMatrix(atomicNumber); + return dam.ConjugateOn(ChemistryHelper.AtomKernel(atomicNumber)); + } + + /// + /// adjoint matrices for atom + /// + public static class AtomAdjoints + { + /// + /// for q + /// + public static readonly SU3 AtomAdjointStateQ = ChemistryHelper.AtomAdjointStateMatrix(1); + /// + /// for r + /// + public static readonly SU3 AtomAdjointStateR = ChemistryHelper.AtomAdjointStateMatrix(2); + /// + /// for es + /// + public static readonly SU3 AtomAdjointStateEs = ChemistryHelper.AtomAdjointStateMatrix(3); + /// + /// + /// + public static readonly SU3 AtomAdjointStateD = ChemistryHelper.AtomAdjointStateMatrix(4); + /// + /// + /// + public static readonly SU3 AtomAdjointStateUe = ChemistryHelper.AtomAdjointStateMatrix(5); + /// + /// + /// + public static readonly SU3 AtomAdjointStateM = ChemistryHelper.AtomAdjointStateMatrix(6); + /// + /// + /// + public static readonly SU3 AtomAdjointStateK = ChemistryHelper.AtomAdjointStateMatrix(7); + /// + /// + /// + public static readonly SU3 AtomAdjointStateP = ChemistryHelper.AtomAdjointStateMatrix(8); + /// + /// + /// + public static readonly SU3 AtomAdjointStateSo = ChemistryHelper.AtomAdjointStateMatrix(9); + /// + /// + /// + public static readonly SU3 AtomAdjointStateE = ChemistryHelper.AtomAdjointStateMatrix(10); + /// + /// + /// + public static readonly SU3 AtomAdjointStateU = ChemistryHelper.AtomAdjointStateMatrix(11); + /// + /// + /// + public static readonly SU3 AtomAdjointStateA = ChemistryHelper.AtomAdjointStateMatrix(12); + /// + /// + /// + public static readonly SU3 AtomAdjointStateCx = ChemistryHelper.AtomAdjointStateMatrix(13); + /// + /// + /// + public static readonly SU3 AtomAdjointStateV = ChemistryHelper.AtomAdjointStateMatrix(14); + /// + /// + /// + public static readonly SU3 AtomAdjointStateVi = ChemistryHelper.AtomAdjointStateMatrix(15); + } + + /// + /// + /// + public static class AtomKernels + { + /// + /// + /// + public static readonly SU3 AtomKernelStateQ = ChemistryHelper.AtomKernel(1); + /// + /// + /// + public static readonly SU3 AtomKernelStateR = ChemistryHelper.AtomKernel(2); + /// + /// + /// + public static readonly SU3 AtomKernelStateEs = ChemistryHelper.AtomKernel(3); + /// + /// + /// + public static readonly SU3 AtomKernelStateD = ChemistryHelper.AtomKernel(4); + /// + /// + /// + public static readonly SU3 AtomKernelStateUe = ChemistryHelper.AtomKernel(5); + /// + /// + /// + public static readonly SU3 AtomKernelStateM = ChemistryHelper.AtomKernel(6); + /// + /// + /// + public static readonly SU3 AtomKernelStateK = ChemistryHelper.AtomKernel(7); + /// + /// + /// + public static readonly SU3 AtomKernelStateP = ChemistryHelper.AtomKernel(8); + /// + /// + /// + public static readonly SU3 AtomKernelStateSo = ChemistryHelper.AtomKernel(9); + /// + /// + /// + public static readonly SU3 AtomKernelStateE = ChemistryHelper.AtomKernel(10); + /// + /// + /// + public static readonly SU3 AtomKernelStateU = ChemistryHelper.AtomKernel(11); + /// + /// + /// + public static readonly SU3 AtomKernelStateA = ChemistryHelper.AtomKernel(12); + /// + /// + /// + public static readonly SU3 AtomKernelStateCx = ChemistryHelper.AtomKernel(13); + /// + /// + /// + public static readonly SU3 AtomKernelStateV = ChemistryHelper.AtomKernel(14); + /// + /// + /// + public static readonly SU3 AtomKernelStateVi = ChemistryHelper.AtomKernel(15); + + } + + /// + /// atom eigen state matrices + /// + public static class AtomEigenStates + { + /// + /// + /// + public static readonly SU3 AtomEigenStateQ = AtomEigenState(1); + /// + /// + /// + public static readonly SU3 AtomEigenStateR = AtomEigenState(2); + /// + /// + /// + public static readonly SU3 AtomEigenStateEs = AtomEigenState(3); + /// + /// + /// + public static readonly SU3 AtomEigenStateD = AtomEigenState(4); + /// + /// + /// + public static readonly SU3 AtomEigenStateUe = AtomEigenState(5); + /// + /// + /// + public static readonly SU3 AtomEigenStateM = AtomEigenState(6); + /// + /// + /// + public static readonly SU3 AtomEigenStateK = AtomEigenState(7); + /// + /// + /// + public static readonly SU3 AtomEigenStateP = AtomEigenState(8); + /// + /// + /// + public static readonly SU3 AtomEigenStateSo = AtomEigenState(9); + /// + /// + /// + public static readonly SU3 AtomEigenStateE = AtomEigenState(10); + /// + /// + /// + public static readonly SU3 AtomEigenStateU = AtomEigenState(11); + /// + /// + /// + public static readonly SU3 AtomEigenStateA = AtomEigenState(12); + /// + /// + /// + public static readonly SU3 AtomEigenStateCx = AtomEigenState(13); + /// + /// + /// + public static readonly SU3 AtomEigenStateV = AtomEigenState(14); + /// + /// + /// + public static readonly SU3 AtomEigenStateVi = AtomEigenState(15); + } + + /// + /// atomic masses for atoms + /// + public static class AtomicMasses + { + /// + /// + /// + public static readonly double AtomicMassQ = AtomicMass(1); + /// + /// + /// + public static readonly double AtomicMassR = AtomicMass(2); + /// + /// + /// + public static readonly double AtomicMassEs = AtomicMass(3); + /// + /// + /// + public static readonly double AtomicMassD = AtomicMass(4); + /// + /// + /// + public static readonly double AtomicMassUe = AtomicMass(5); + /// + /// + /// + public static readonly double AtomicMassM = AtomicMass(6); + /// + /// + /// + public static readonly double AtomicMassK = AtomicMass(7); + /// + /// + /// + public static readonly double AtomicMassP = AtomicMass(8); + /// + /// + /// + public static readonly double AtomicMassSo = AtomicMass(9); + /// + /// + /// + public static readonly double AtomicMassE = AtomicMass(10); + /// + /// + /// + public static readonly double AtomicMassU = AtomicMass(11); + /// + /// + /// + public static readonly double AtomicMassA = AtomicMass(12); + /// + /// + /// + public static readonly double AtomicMassCx = AtomicMass(13); + /// + /// + /// + public static readonly double AtomicMassV = AtomicMass(14); + /// + /// + /// + public static readonly double AtomicMassVi = AtomicMass(15); + + } + + private static SU3 BondEigenState(int bondNumber) + { + SU3 dam = ChemistryHelper.BondAdjointStateMatrix(bondNumber); + return dam.ConjugateOn(ChemistryHelper.BondKernel(bondNumber)); + } + /// + /// kernel matrices for bonds + /// + public static class BondKernels + { + /// + /// + /// + public static readonly SU3 BondKernelM = ChemistryHelper.BondKernel(1); + /// + /// + /// + public static readonly SU3 BondKernelS = ChemistryHelper.BondKernel(2); + /// + /// + /// + public static readonly SU3 BondKernelD = ChemistryHelper.BondKernel(3); + /// + /// + /// + public static readonly SU3 BondKernelY = ChemistryHelper.BondKernel(4); + /// + /// + /// + public static readonly SU3 BondKernelQ = ChemistryHelper.BondKernel(5); + /// + /// + /// + public static readonly SU3 BondKernelP = ChemistryHelper.BondKernel(6); + /// + /// + /// + public static readonly SU3 BondKernelH = ChemistryHelper.BondKernel(7); + } + /// + /// adjoint matrices for bond + /// + public static class BondAdjoints + { + /// + /// + /// + public static readonly SU3 BondAdjointM = ChemistryHelper.BondAdjointStateMatrix(1); + /// + /// + /// + public static readonly SU3 BondAdjointS = ChemistryHelper.BondAdjointStateMatrix(2); + /// + /// + /// + public static readonly SU3 BondAdjointD = ChemistryHelper.BondAdjointStateMatrix(3); + /// + /// + /// + public static readonly SU3 BondAdjointY = ChemistryHelper.BondAdjointStateMatrix(4); + /// + /// + /// + public static readonly SU3 BondAdjointQ = ChemistryHelper.BondAdjointStateMatrix(5); + /// + /// + /// + public static readonly SU3 BondAdjointP = ChemistryHelper.BondAdjointStateMatrix(6); + /// + /// + /// + public static readonly SU3 BondAdjointH = ChemistryHelper.BondAdjointStateMatrix(7); + } + + /// + /// eigen state matrices for bond + /// + public static class BondEigenStates + { + /// + /// + /// + public static readonly SU3 BondEigenStateM = BondEigenState(1); + /// + /// + /// + public static readonly SU3 BondEigenStateS = BondEigenState(2); + /// + /// + /// + public static readonly SU3 BondEigenStateD = BondEigenState(3); + /// + /// + /// + public static readonly SU3 BondEigenStateY = BondEigenState(4); + /// + /// + /// + public static readonly SU3 BondEigenStateQ = BondEigenState(5); + /// + /// + /// + public static readonly SU3 BondEigenStateP = BondEigenState(6); + /// + /// + /// + public static readonly SU3 BondEigenStateH = BondEigenState(7); + + } + + /// + /// max bond numbers for specific bond number + /// + /// + /// + public static int BondMaxNumbers(int bondNumber) => bondNumber switch + { + 1 => 1, + 2 => 2, + 3 => 1, + 4 => 2, + 5 => 2, + 6 => 1, + 7 => 2, + _ => -1 + }; + + /// + /// filters to calculate the color of mixtures + /// + public static class OpticalFilter + { + /// + /// to calculate redness + /// + public static readonly U3 RedFilter = new LieU3 + ( + I,0, -1, + 0, I, 0, + 1, 0, -I + ).Exp(); + + /// + /// to calculate greenness + /// + public static readonly U3 GreenFilter = new LieU3 + ( + 0, -I, 0, + -I, 0, -1, + 0, 1, 0 + ).Exp(); + /// + /// to calculate blueness + /// + public static readonly U3 BlueFilter = new LieU3 + ( + 0, 0, -I, + 0, 0, -1, + -I, 1, 0 + ).Exp(); + /// + /// to calculate alpha channel + /// + public static readonly U3 OpacityFilter = new LieU3( + I, -I, -I, + -I, I, -1, + -I, 1, -I + ).Exp(); + + /// + /// red + /// + public static readonly H3 RedDefinite = RedFilter.ConjugateOn(AlgebraConstant.UniformSpectrum3); + /// + /// green + /// + public static readonly H3 GreenDefinite = GreenFilter.ConjugateOn(AlgebraConstant.UniformSpectrum3); + /// + /// blue + /// + public static readonly H3 BlueDefinite = BlueFilter.ConjugateOn(AlgebraConstant.UniformSpectrum3); + /// + /// alpha + /// + public static readonly H3 OpacityDefinite = OpacityFilter.ConjugateOn(new DiagR3(0.05d, 0.74d, 1d)); + } + + + /// + /// all elements + /// + public static readonly string[] ElementSymbols = + { + "Q", "R", "Es", "D", "Ue", + "M", "K", "P", "So", "E", + "U", "A", "Cx", "V", "Vi" + }; + + + /// + /// x + /// + public static class CommonResources + { + /// + /// x + /// + public const string Letroline = "res://Data/Chemicals/HighEnergyMixtureLetroline.btx"; + } + + /// + /// Used to measure Cayley value of homogeneous mixture + /// + public static readonly C3 CayleyMeasure = new (1 - I, I, -1 + I); + /// + /// Used to measure Euclid value of homogeneous mixture + /// + public static readonly C3 EuclidMeasure = new(1, I, -I); + +} \ No newline at end of file diff --git a/src/DataStructure/Caches/Cache.cs b/src/DataStructure/Caches/Cache.cs new file mode 100644 index 0000000..5ee00be --- /dev/null +++ b/src/DataStructure/Caches/Cache.cs @@ -0,0 +1,54 @@ +using Skeleton.Utils.Helpers; + +namespace VirtualChemistry.DataStructure.Caches; + +/// +/// cache to store a specific type of data +/// +/// +public class Cache : Dictionary, ICache +{ + /// + public Cache() + { + UpdateRequest = new HashSet(); + } + + /// + /// indexer + /// + /// + /// re calculate by getter if key is expired of no value stored in the specified key + public TValue this[string key, Func getter] + { + get + { + if(ValueOutdated(key)) + { + UpdateRequest.Remove(key); + return this[key] = getter(); + } + return this[key]; + } + } + + /// + public void DeleteCache() + { + Keys.ForEach(x => Remove(x)); + UpdateRequest = new HashSet(); + } + + /// + public void DeleteCache(string key) => Remove(key); + private HashSet UpdateRequest { get; set; } + private bool ValueOutdated(string key) => !ContainsKey(key) || UpdateRequest.Contains(key); + + /// + public void Expire(string key) => UpdateRequest.Add(key); + + /// + public void Expire(IEnumerable keys) => keys.ForEach(Expire); + + +} diff --git a/src/DataStructure/Caches/GeneralCache.cs b/src/DataStructure/Caches/GeneralCache.cs new file mode 100644 index 0000000..8dd0af4 --- /dev/null +++ b/src/DataStructure/Caches/GeneralCache.cs @@ -0,0 +1,156 @@ +using System.Numerics; +using Skeleton.Utils.Helpers; + +namespace VirtualChemistry.DataStructure.Caches; + +/// +/// a general propose cache to store many type of tensors +/// +public class GeneralCache +{ + private Cache C3Cache { get; set; } = new(); + private Cache C33Cache { get; set; } = new(); + private Cache ComplexCache { get; set; } = new(); + private Cache DoubleCache { get; set; } = new(); + private Cache U3Cache { get; set; } = new(); + private Cache SU3Cache { get; set; } = new(); + private Cache LieU3Cache { get; set; } = new(); + private Cache LieSU3Cache { get; set; } = new(); + + private static T Identity(T x) => x; + private Func RegisterOnExpire(string key, Action? e) + { + + if (e != null && !OnExpire.ContainsKey(key)) + OnExpire[key] = e; + return Identity; + } + + /// + /// C3 tensor getter + /// + /// + /// + /// a function to expire other keys if this key is expired + public C3 this[string key, Func getter, Action? onExpire = null] => + RegisterOnExpire(key, onExpire)(C3Cache[key, getter]); + + /// + /// C33 tensor getter + /// + /// + /// + /// + public C33 this[string key, Func getter, Action? onExpire = null] => + RegisterOnExpire(key, onExpire)(C33Cache[key, getter]); + /// + /// complex getter + /// + /// + /// + /// + public Complex this[string key, Func getter, Action? onExpire = null] => + RegisterOnExpire(key, onExpire)(ComplexCache[key, getter]); + /// + /// real getter + /// + /// + /// + /// + public double this[string key, Func getter, Action? onExpire = null] => + RegisterOnExpire(key, onExpire)(DoubleCache[key, getter]); + + /// + /// U3 getter + /// + /// + /// + /// + public U3 this[string key, Func getter, Action? onExpire = null] => + RegisterOnExpire(key, onExpire)(U3Cache[key, getter]); + + /// + /// SU3 getter + /// + /// + /// + /// + public SU3 this[string key, Func getter, Action? onExpire = null] => + RegisterOnExpire(key, onExpire)(SU3Cache[key, getter]); + + /// + /// LieU3 getter + /// + /// + /// + /// + public LieU3 this[string key, Func getter, Action? onExpire = null] => + RegisterOnExpire(key, onExpire)(LieU3Cache[key, getter]); + + /// + /// LieSU3 getter + /// + /// + /// + /// + public LieSU3 this[string key, Func getter, Action? onExpire = null] => + RegisterOnExpire(key, onExpire)(LieSU3Cache[key, getter]); + + /// + /// general propose getter + /// + /// + /// + public ICache this[string key] => key.Split('_')[0] switch + { + "D" => DoubleCache, + "C" => ComplexCache, + "C3" => C3Cache, + "C33" => C33Cache, + "U3" => U3Cache, + "SU3" => SU3Cache, + "LieU3" => LieU3Cache, + "LieSU3" => LieSU3Cache, + _ => throw new Exception() + }; + + /// + /// general expire + /// + /// + public void Expire(string key) + { + if (OnExpire.ContainsKey(key)) + OnExpire[key](); + this[key].Expire(key); + } + + /// + /// general expire + /// + /// + public void Expire(IEnumerable keys) => keys.ForEach(Expire); + /// + /// general delete cache + /// + /// + public void DeleteCache(string key) => this[key].DeleteCache(key); + + /// + /// trigger the action when key is expired + /// + public Dictionary OnExpire { get; set; } = new(); + + /// + /// get raw real(constant) + /// + /// + /// + public double GetRaw(string key) => DoubleCache.GetValueOrDefault(key); + /// + /// init a constant real + /// + /// + /// + public void SetRaw(string key, double value) => DoubleCache[key] = value; +} \ No newline at end of file diff --git a/src/DataStructure/Caches/ICache.cs b/src/DataStructure/Caches/ICache.cs new file mode 100644 index 0000000..c0ce4a1 --- /dev/null +++ b/src/DataStructure/Caches/ICache.cs @@ -0,0 +1,27 @@ +namespace VirtualChemistry.DataStructure.Caches; + +/// +/// interface for cache +/// +public interface ICache +{ + /// + /// delete all cache + /// + public void DeleteCache(); + /// + /// delete for a specific key + /// + /// + public void DeleteCache(string key); + /// + /// expire a key, the value need to be re calculated from getter + /// + /// + public void Expire(string key); + /// + /// expire a set of keys + /// + /// + public void Expire(IEnumerable keys); +} \ No newline at end of file diff --git a/src/DataStructure/Packs/BondGroups.cs b/src/DataStructure/Packs/BondGroups.cs new file mode 100644 index 0000000..8bdc2c8 --- /dev/null +++ b/src/DataStructure/Packs/BondGroups.cs @@ -0,0 +1,33 @@ +using Skeleton.DataStructure.Packs; +using VirtualChemistry.Chemistry.Bonds.Implements; + +namespace VirtualChemistry.DataStructure.Packs; + +using BondGroup = HashSet; + +/// +public class BondGroups : Pack4 +{ + /// + /// + /// + public BondGroup Group1 => Item1; + /// + /// + /// + public BondGroup Group2 => Item2; + /// + /// + /// + public BondGroup Group3 => Item3; + /// + /// + /// + public BondGroup Group4 => Item4; + + /// + public BondGroups(BondGroup g1, BondGroup g2, BondGroup g3, BondGroup g4) : base(g1, g2, g3, g4) + { + + } +} \ No newline at end of file diff --git a/src/DataStructure/Packs/CompoundMap.cs b/src/DataStructure/Packs/CompoundMap.cs new file mode 100644 index 0000000..0679b5d --- /dev/null +++ b/src/DataStructure/Packs/CompoundMap.cs @@ -0,0 +1,67 @@ +using Skeleton.DataStructure.Packs; +using VirtualChemistry.Chemistry.Atoms.Implements; +using VirtualChemistry.Chemistry.Bonds.Implements; +using VirtualChemistry.Chemistry.Compounds.Implements; + +namespace VirtualChemistry.DataStructure.Packs; + +/// +/// compound dump map +/// +public class CompoundDumpMap : Pack3, Dictionary, string> +{ + /// + /// atom map + /// + public Dictionary AtomMap => Item1; + /// + /// bond map + /// + public Dictionary BondMap => Item2; + /// + /// save string + /// + public string DumpString => Item3; + + /// + /// constructor + /// + /// + /// + /// + public CompoundDumpMap(Dictionary atomMap, Dictionary bondMap, string dumpString) : + base(atomMap, bondMap, dumpString) + { + } +} + +/// +/// resolve map for compound +/// +public class CompoundResolveMap : Pack3, Dictionary, Compound> +{ + /// + /// atom map + /// + public Dictionary AtomMap => Item1; + /// + /// bond map + /// + public Dictionary BondMap => Item2; + /// + /// created compound + /// + public Compound ResolvedCompound => Item3; + + /// + /// constructor + /// + /// + /// + /// + public CompoundResolveMap(Dictionary atomMap, Dictionary bondMap, + Compound resolvedCompound) : + base(atomMap, bondMap, resolvedCompound) + { + } +} \ No newline at end of file diff --git a/src/DataStructure/Packs/ConnectionInfo.cs b/src/DataStructure/Packs/ConnectionInfo.cs new file mode 100644 index 0000000..4eab789 --- /dev/null +++ b/src/DataStructure/Packs/ConnectionInfo.cs @@ -0,0 +1,37 @@ +using Skeleton.DataStructure.Packs; +using VirtualChemistry.Chemistry.Atoms.Implements; +using VirtualChemistry.Chemistry.Bonds.Implements; + +namespace VirtualChemistry.DataStructure.Packs; + +/// +/// connection info of an atom and a bond +/// +public class ConnectionInfo : Pack3 +{ + /// + /// constructor + /// + /// + /// + /// + public ConnectionInfo(BaseBond bondFrom, BaseBond bondTo, BaseAtom connectedAtom) + { + Item1 = bondFrom; + Item2 = bondTo; + Item3 = connectedAtom; + } + + /// + /// from bond + /// + public BaseBond BondFrom => Item1; + /// + /// to bond + /// + public BaseBond BondTo => Item2; + /// + /// to atom + /// + public BaseAtom ConnectedAtom => Item3; +} \ No newline at end of file diff --git a/src/DataStructure/Packs/FullConnectionInfo.cs b/src/DataStructure/Packs/FullConnectionInfo.cs new file mode 100644 index 0000000..77076ca --- /dev/null +++ b/src/DataStructure/Packs/FullConnectionInfo.cs @@ -0,0 +1,64 @@ +using Skeleton.DataStructure.Packs; +using VirtualChemistry.Chemistry.Atoms.Implements; +using VirtualChemistry.Chemistry.Bonds.Implements; + +namespace VirtualChemistry.DataStructure.Packs; + +/// +/// all connections from one atom to another +/// +public class FullConnectionInfo : Pack3, BaseAtom> +{ + /// + /// constructor + /// + /// + public FullConnectionInfo(BaseBond bond) : base(bond.Atom, new HashSet() { bond }, bond.ConnectedBond.Atom) + { + } + + /// + /// from atom + /// + public BaseAtom Atom1 => Item1; + /// + /// connected bonds + /// + public HashSet Bonds => Item2; + /// + /// bonds of atom 1 + /// + public BaseAtom Atom2 => Item3; + private string OnwRepresentation => $"{Atom1.Element}-[{String.Join('|', Bonds.Select(bond => $"{bond.BondType}-{bond.ConnectedBond.BondType}"))}]-{Atom2.Element}"; + private string RevRepresentation => $"{Atom2.Element}-[{String.Join('|', Bonds.Select(bond => $"{bond.ConnectedBond.BondType}-{bond.BondType}"))}]-{Atom1.Element}"; + + /// + /// save string + /// + public string Representation => Atom1.AtomicNumber + Bonds.Sum(bond => bond.BondNumber) > + Atom2.AtomicNumber + Bonds.Sum(bond => bond.ConnectedBond.BondNumber) + ? OnwRepresentation + : RevRepresentation; + + /// + /// try to add a bond to bonds + /// + /// + /// + public bool TryAbsorbBond(BaseBond bond) + { + if(bond.Atom == Atom1 && bond.ConnectedBond.Atom == Atom2) + { + Bonds.Add(bond); + return true; + } + + if (bond.Atom == Atom2 && bond.ConnectedBond.Atom == Atom1) + { + Bonds.Add(bond.ConnectedBond); + return true; + } + + return false; + } +} \ No newline at end of file diff --git a/src/Using.cs b/src/Using.cs new file mode 100644 index 0000000..1edca34 --- /dev/null +++ b/src/Using.cs @@ -0,0 +1,42 @@ +global using C2 = Skeleton.Algebra.CategoryOf.OnField.FVector; +global using C3 = Skeleton.Algebra.CategoryOf.OnField.FVector; +global using C4 = Skeleton.Algebra.CategoryOf.OnField.FVector; + +global using R2 = Skeleton.Algebra.CategoryOf.OnField.FVector; +global using R3 = Skeleton.Algebra.CategoryOf.OnField.FVector; +global using R4 = Skeleton.Algebra.CategoryOf.OnField.FVector; + +global using C22 = Skeleton.Algebra.CategoryOf.OnField.FMatrix; +global using C33 = Skeleton.Algebra.CategoryOf.OnField.FMatrix; +global using C44 = Skeleton.Algebra.CategoryOf.OnField.FMatrix; + +global using R22 = Skeleton.Algebra.CategoryOf.OnField.FMatrix; +global using R33 = Skeleton.Algebra.CategoryOf.OnField.FMatrix; +global using R44 = Skeleton.Algebra.CategoryOf.OnField.FMatrix; + +global using U2 = Skeleton.Algebra.CategoryOf.FUnitaryMatrix; +global using U3 = Skeleton.Algebra.CategoryOf.FUnitaryMatrix; +global using U4 = Skeleton.Algebra.CategoryOf.FUnitaryMatrix; + +global using SU2 = Skeleton.Algebra.CategoryOf.FSpecialUnitaryMatrix; +global using SU3 = Skeleton.Algebra.CategoryOf.FSpecialUnitaryMatrix; +global using SU4 = Skeleton.Algebra.CategoryOf.FSpecialUnitaryMatrix; + +global using LieU2 = Skeleton.Algebra.CategoryOf.FLieUnitaryMatrix; +global using LieU3 = Skeleton.Algebra.CategoryOf.FLieUnitaryMatrix; +global using LieU4 = Skeleton.Algebra.CategoryOf.FLieUnitaryMatrix; + +global using LieSU2 = Skeleton.Algebra.CategoryOf.FSpecialLieUnitaryMatrix; +global using LieSU3 = Skeleton.Algebra.CategoryOf.FSpecialLieUnitaryMatrix; +global using LieSU4 = Skeleton.Algebra.CategoryOf.FSpecialLieUnitaryMatrix; + +global using C2Space = Skeleton.Algebra.CategoryOf.OnField.FVectorSpace; +global using C3Space = Skeleton.Algebra.CategoryOf.OnField.FVectorSpace; +global using C4Space = Skeleton.Algebra.CategoryOf.OnField.FVectorSpace; + +global using R2Space = Skeleton.Algebra.CategoryOf.OnField.FVectorSpace; +global using R3Space = Skeleton.Algebra.CategoryOf.OnField.FVectorSpace; +global using R4Space = Skeleton.Algebra.CategoryOf.OnField.FVectorSpace; + +global using DiagR3 = Skeleton.Algebra.CategoryOf.OnField.FDiagonalMatrix; +global using H3 = Skeleton.Algebra.CategoryOf.FHermitianMatrix; \ No newline at end of file diff --git a/src/Utils/Helpers/ChemistryHelper.cs b/src/Utils/Helpers/ChemistryHelper.cs new file mode 100644 index 0000000..244de38 --- /dev/null +++ b/src/Utils/Helpers/ChemistryHelper.cs @@ -0,0 +1,282 @@ +using System.Numerics; + +namespace VirtualChemistry.Utils.Helpers; + +/// +/// helper functions for chemistry +/// +public static class ChemistryHelper +{ + private static readonly Complex I = Complex.ImaginaryOne; + /// + /// + /// + /// + /// + public static double BondHeightOrder(int n) => Math.Pow(4d / 5d, n) - 1; + /// + /// + /// + /// + /// + public static double BondWeightOrder(int n) + { + double ra = n / 11d; + int det = (int)Math.Floor(ra / 0.5d); + if (det % 2 == 0) + return ra % 0.5d; + return 0.5 - (ra % 0.5d); + } + + /// + /// + /// + /// + /// + public static Complex BondHeight(int n) => Complex.Exp(I * Math.PI * BondHeightOrder(n)); + /// + /// + /// + /// + /// + public static Complex BondWeight(int n) => Complex.Exp(I * Math.PI * BondWeightOrder(n)); + + /// + /// + /// + /// + /// + public static Complex BondE1(int n) => BondHeight(n); + /// + /// + /// + /// + /// + public static Complex BondE2(int n) => BondWeight(n); + /// + /// + /// + /// + /// + public static Complex BondE3(int n) => 1 / (BondE1(n)*BondE2(n)); + /// + /// + /// + /// + /// + public static SU3 BondKernel(int p) => new + ( + BondE1(p), 0, 0, + 0, BondE2(p), 0, + 0, 0, BondE3(p) + ); + + + + private static readonly double Sqrt13 = Math.Sqrt(13d); + private static readonly double Sqrt7 = Math.Sqrt(7d); + private static readonly double Sqrt3 = Math.Sqrt(3d); + private static readonly double Sqrt2 = Math.Sqrt(2d); + private static readonly double Sqrt21 = Sqrt3 * Sqrt7; + + + /// + /// + /// + public static readonly Complex BondAdjointEigen1 = I * Math.Sqrt(0.5d * (5d + Sqrt21)); + /// + /// + /// + public static readonly Complex BondAdjointEigen2 = I * Math.Sqrt(0.5d * (5d - Sqrt21)); + + + /// + /// + /// + public static readonly C3 BondAdjointEigenVector1 = new C3 + ( + 1d/Sqrt2, + -I/2d, + 1/2 + ); + + /// + /// + /// + public static readonly C3 BondAdjointEigenVector2 = new C3 + ( + 0, + I/Sqrt2, + 1d/Sqrt2 + ); + + /// + /// + /// + public static readonly C3 BondAdjointEigenVector3 = new C3 + ( + -1d/Sqrt2, + I/Sqrt2, + 1d/2d + ); + + + /// + /// + /// + public static readonly C33 BondAdjointCoMatrix = + new C33(new[] { BondAdjointEigenVector1, BondAdjointEigenVector2, BondAdjointEigenVector3 }, false); + + /// + /// + /// + /// + /// + public static SU3 BondAdjointStateMatrix(int bondNumber) + { + double t = Math.Cos(bondNumber * Math.E) + Sqrt3; + Complex a1 = Complex.Exp(t * BondAdjointEigen1); + Complex a2 = Complex.Exp(t * BondAdjointEigen2); + Complex a3 = 1 / (a1 * a2); + C33 expBondAdjointEigenDiagonal = new + ( + a1, 0d, 0d, + 0d, a2, 0d, + 0d, 0d, a3 + ); + return new SU3(BondAdjointCoMatrix * expBondAdjointEigenDiagonal * BondAdjointCoMatrix.Inv()); + } + + + /// + /// + /// + /// + /// + public static double AtomHeightOrder(int n) => 1 - Math.Pow(6d / 7d, n); + /// + /// + /// + /// + /// + public static double AtomWeightOrder(int n) + { + double ra = n / 7d; + + if ((int)Math.Floor(ra / 0.5) % 2 == 0) + return ra % 0.5d; + return 0.5d - (ra % 0.5d); + } + + /// + /// + /// + /// + /// + public static Complex AtomHeight(int n) => Complex.Exp(I * Math.PI * AtomHeightOrder(n)); + /// + /// + /// + /// + /// + public static Complex AtomWeight(int n) => Complex.Exp(I * Math.PI * AtomWeightOrder(n)); + + /// + /// + /// + /// + /// + public static Complex AtomE1(int n) => AtomHeight(n); + /// + /// + /// + /// + /// + public static Complex AtomE2(int n) => AtomWeight(n); + /// + /// + /// + /// + /// + public static Complex AtomE3(int n) => 1 / (AtomE1(n) * AtomE2(n)); + + /// + /// + /// + /// + /// + public static SU3 AtomKernel(int atomicNumber) => new + ( + AtomE1(atomicNumber), 0, 0, + 0, AtomE2(atomicNumber), 0, + 0, 0, AtomE3(atomicNumber) + ); + + /// + /// + /// + public static readonly Complex AtomAdjointEigen1 = I * Complex.Sqrt(0.5 * (5 + Sqrt21)); + /// + /// + /// + public static readonly Complex AtomAdjointEigen2 = I * Complex.Sqrt(0.5 * (5 - Sqrt21)); + + /// + /// + /// + public static readonly C3 AtomAdjointEigenVector1 = new C3 + ( + 0d, + -I * (Sqrt2 - 1d) / (Complex.Sqrt(4d - 2d * Sqrt2)), + 1d / (Complex.Sqrt(4 - 2 * Sqrt2)) + ); + + /// + /// + /// + public static readonly C3 AtomAdjointEigenVector2 = new C3 + ( + 1d, + 0d, + 0d + ); + + /// + /// + /// + public static readonly C3 AtomAdjointEigenVector3 = new C3 + ( + 0d, + I * (1d + Sqrt2) / Complex.Sqrt(2d * (2d + Sqrt2)), + 1d / Complex.Sqrt(2d * (2d + Sqrt2)) + ); + + + /// + /// + /// + public static readonly SU3 AtomAdjointCoMatrix = + new ( new []{ AtomAdjointEigenVector1, AtomAdjointEigenVector2, AtomAdjointEigenVector3 }, false); + + /// + /// + /// + /// + /// + public static SU3 AtomAdjointStateMatrix(int atomicNumber) + { + double t = Math.Cos(atomicNumber * Math.E) + Sqrt3; + Complex a1 = Complex.Exp(t * AtomAdjointEigen1); + Complex a2 = Complex.Exp(t * AtomAdjointEigen2); + Complex a3 = 1 / (a1 * a2); + SU3 expAtomAdjointEigenDiagonal = new SU3 + ( + a1, 0, 0, + 0, a2 , 0, + 0, 0, a3 + ); + return AtomAdjointCoMatrix.ConjugateOn(expAtomAdjointEigenDiagonal); + } + + +}