Created
May 14, 2024 22:38
-
-
Save philipturner/38ee141923ab2373afa6b338ad3b76ab to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
import Foundation | |
import HDL | |
import MM4 | |
import Numerics | |
import OpenMM | |
func createGeometry() -> [Entity] { | |
var hexagon = Hexagon() | |
hexagon.center() | |
return hexagon.topology.atoms | |
} | |
struct Hexagon { | |
var topology: Topology | |
init() { | |
let lattice = Self.createLattice() | |
topology = Self.createTopology(lattice: lattice) | |
} | |
static func createLattice() -> Lattice<Hexagonal> { | |
Lattice<Hexagonal> { h, k, l in | |
let h2k = h + 2 * k | |
Bounds { 16 * h + 10 * h2k + 2 * l } | |
Material { .elemental(.carbon) } | |
Volume { | |
Origin { 8 * h + 5 * h2k } | |
var directions: [SIMD3<Float>] = [] | |
directions.append(k) | |
directions.append(h) | |
directions.append(-k) | |
directions.append(-h) | |
directions.append(k + h) | |
directions.append(-k - h) | |
for direction in directions { | |
Convex { | |
Origin { 5 * direction } | |
Plane { direction } | |
} | |
} | |
Replace { .empty } | |
} | |
} | |
} | |
static func createTopology(lattice: Lattice<Hexagonal>) -> Topology { | |
var topology = Topology() | |
topology.insert(atoms: lattice.atoms) | |
guard topology.bonds.count == .zero else { | |
fatalError("The topology should not have any bonds yet.") | |
} | |
// Encapsulate 'insertedBonds' in a scope. This way, you cannot reference | |
// the variable name (by accident) later in the code. | |
do { | |
let matches = topology.match( | |
topology.atoms, algorithm: .covalentBondLength(1.5)) | |
var insertedBonds: [SIMD2<UInt32>] = [] | |
for i in topology.atoms.indices { | |
for j in matches[i] where i < j { | |
let bond = SIMD2(UInt32(i), UInt32(j)) | |
insertedBonds.append(bond) | |
} | |
} | |
topology.insert(bonds: insertedBonds) | |
} | |
guard topology.bonds.count > .zero else { | |
fatalError("No bonds were added.") | |
} | |
// Encapsulate 'nonbondingOrbitals', 'insertedAtoms', and 'insertedBonds' | |
// within their own scope, so they will be deleted as soon as they're not | |
// relevant. | |
do { | |
// 'nonbondingOrbitals' is a list of lists. | |
let nonbondingOrbitals = topology.nonbondingOrbitals(hybridization: .sp3) | |
var insertedAtoms: [Entity] = [] | |
var insertedBonds: [SIMD2<UInt32>] = [] | |
for atomID in topology.atoms.indices { | |
let carbon = topology.atoms[atomID] | |
let carbonPosition = carbon.position | |
// Here, you have fetched a list from the 'list of lists'. | |
let orbitals = nonbondingOrbitals[atomID] | |
for orbital in orbitals { | |
// Get a rough estimate of the bond length. I often substitute this | |
// with precise values from the literature, or MM3/MM4 equilibrium | |
// bond lengths. Usually it doesn't matter too much, because | |
// minimization will relax it to the exact value. | |
let chBondLength = Element.carbon.covalentRadius + | |
Element.hydrogen.covalentRadius | |
// Position + Length * (Direction Vector) | |
// | |
// The direction vector is normalized (unit length), and | |
// 'chBondLength' effectively gives it a length. | |
let hydrogenPosition = carbonPosition + chBondLength * orbital | |
let hydrogen = Entity( | |
position: hydrogenPosition, type: .atom(.hydrogen)) | |
// Predict the passivator's index in the topology, after the current | |
// atoms are merged with the 'inserted atoms'. | |
let hydrogenID = topology.atoms.count + insertedAtoms.count | |
let bond = SIMD2(UInt32(atomID), UInt32(hydrogenID)) | |
insertedAtoms.append(hydrogen) | |
insertedBonds.append(bond) | |
} | |
} | |
topology.insert(atoms: insertedAtoms) | |
topology.insert(bonds: insertedBonds) | |
} | |
return topology | |
} | |
mutating func center() { | |
// Accumulate in double precision to minimize rounding error. | |
var weightedCenterAccumulator: SIMD3<Double> = .zero | |
var massAccumulator: Double = .zero | |
// Iterate over the atoms. | |
for atomID in topology.atoms.indices { | |
let atom = topology.atoms[atomID] | |
// Approximate the atomic mass by using the atomic number (larger atoms | |
// are heavier). We don't have access to atomic masses this early on in | |
// compilation. | |
let mass = Float(atom.atomicNumber) | |
let position: SIMD3<Float> = atom.position | |
weightedCenterAccumulator += SIMD3<Double>(position * mass) | |
massAccumulator += Double(mass) | |
} | |
// Cast the final sum to single precision, because we don't need the extra | |
// precision anymore. | |
let weightedCenter = SIMD3<Float>(weightedCenterAccumulator) | |
let mass = Float(massAccumulator) | |
let centerOfMass = weightedCenter / mass | |
// Iterate over the atoms. | |
for atomID in topology.atoms.indices { | |
var atom = topology.atoms[atomID] | |
atom.position -= centerOfMass | |
topology.atoms[atomID] = atom | |
} | |
} | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment