Add molecular aligning external field

examples/ortho-terphenyl/5-single-molecule-field-test/.gitignore

@@ -0,0 +1 @@
+output/

examples/ortho-terphenyl/5-single-molecule-field-test/README.md

@@ -0,0 +1,55 @@
+# Single-molecule aligning-field check
+
+This is a minimal validation of the aligning term
+
+\[
+U_h=-\mathbf h\cdot\mathbf n
+\]
+
+for the three-site ortho-terphenyl model. The system contains one isolated
+molecule and uses displacement moves only. Its orientation is the ordered plane
+normal from `PlaneNormalOrientation(1, 2, 3)`. The field points along `+z`, while
+the initial molecular normal is perpendicular to it.
+
+Run the scan from the repository root:
+
+```bash
+julia --project=. examples/ortho-terphenyl/5-single-molecule-field-test/run-field-scan.jl
+```
+
+The script runs independent simulations for
+`|h| = 0, 0.5, 1, 3, 10` at `T = 1`, samples after burn-in, and writes
+`alignment-vs-field.csv`. The measured quantity is
+
+\[
+\langle \hat{\mathbf h}\cdot\mathbf n\rangle.
+\]
+
+It should be close to zero without a field and should approach one as the field
+magnitude increases. Finite sampling means individual estimates need not be
+strictly monotonic. This test intentionally uses a one-molecule system; applying
+a field to a selected molecule inside the full liquid would require a separate
+molecule-selection option.
+
+One run with the default seed produced:
+
+| field magnitude | mean alignment | standard error |
+|---:|---:|---:|
+| 0.0 | 0.0310 | 0.0129 |
+| 0.5 | 0.3134 | 0.0117 |
+| 1.0 | 0.3415 | 0.0113 |
+| 3.0 | 0.6584 | 0.0075 |
+| 10.0 | 0.8995 | 0.0022 |
+
+The detailed values are stored in `alignment-vs-field.csv`.
+
+An XYZ trajectory is also saved for each field magnitude at:
+
+```text
+output/h_FIELD/chains/1/trajectory.xyz
+```
+
+For example, the strongest-field trajectory is
+`output/h_10.0/chains/1/trajectory.xyz`. Frames are written every 100 MC steps,
+including the initial configuration, and can be opened with OVITO, VMD, or
+another XYZ-compatible molecular viewer.

examples/ortho-terphenyl/5-single-molecule-field-test/alignment-vs-field.csv

@@ -0,0 +1,6 @@
+h,mean_alignment,standard_error,n_samples
+0.0,0.031035990242472784,0.012856472952436386,2002
+0.5,0.31335135785002893,0.011689660988449533,2002
+1.0,0.34146412397251336,0.011250510226195084,2002
+3.0,0.658383296722763,0.007500344387261211,2002
+10.0,0.8994559221111872,0.0022077828788680583,2002

examples/ortho-terphenyl/5-single-molecule-field-test/run-field-scan.jl

@@ -0,0 +1,92 @@
+using Arianna
+using ComponentArrays
+using LinearAlgebra
+using ParticlesMC
+using Statistics
+using StaticArrays
+
+const FIELD_DIRECTION = SVector(0.0, 0.0, 1.0)
+
+# This callback deliberately measures molecule 1: the example contains exactly
+# one molecule, so no molecule-selection feature is needed for this first test.
+Arianna.@callback function field_alignment(system::Molecules)
+    n = orientation(PlaneNormalOrientation(), system, 1)
+    return dot(FIELD_DIRECTION, n)
+end
+
+function isolated_otp(field_magnitude)
+    # The initial triangle lies in the yz plane, so its normal starts
+    # perpendicular to the field direction (+z).
+    positions = [
+        SVector(2.0, 2.0, 2.0),
+        SVector(2.0, 2.8, 2.0),
+        SVector(2.0, 2.4, 2.7),
+    ]
+    species = [1, 2, 3]
+    molecule = ones(Int, 3)
+    bonds = [[2, 3], [1, 3], [1, 2]]
+    box_length = 5.0
+    density = length(positions) / box_length^3
+    temperature = 1.0
+    field = AligningField(
+        field_magnitude * FIELD_DIRECTION,
+        PlaneNormalOrientation(),
+    )
+
+    return System(
+        positions,
+        species,
+        molecule,
+        density,
+        temperature,
+        Trimer(),
+        bonds;
+        masses=ones(3),
+        external_field=field,
+        list_type=EmptyList,
+    )
+end
+
+function read_callback_values(path)
+    lines = filter(!isempty, readlines(path))
+    return [parse(Float64, split(line)[2]) for line in lines]
+end
+
+function run_one_field(field_magnitude;
+                       steps=50_000,
+                       burn=10_000,
+                       interval=20,
+                       trajectory_interval=100)
+    system = isolated_otp(field_magnitude)
+    policy = SimpleGaussian()
+    parameters = ComponentArray(σ=0.05)
+    pool = (
+        Move(Displacement(0, zero(system.box), 0.0), policy, parameters, 1.0),
+    )
+    scheduler = build_schedule(steps, burn, interval)
+    trajectory_scheduler = build_schedule(steps, 0, trajectory_interval)
+    output = joinpath(@__DIR__, "output", "h_$(field_magnitude)")
+    algorithms = (
+        (algorithm=Metropolis, pool=pool, seed=1234, parallel=false, sweepstep=system.N),
+        (algorithm=StoreCallbacks, callbacks=(field_alignment,), scheduler=scheduler),
+        (algorithm=StoreTrajectories, scheduler=trajectory_scheduler, fmt=XYZ()),
+    )
+    simulation = Simulation([system], algorithms, steps; path=output, verbose=false)
+    run!(simulation)
+
+    values = read_callback_values(joinpath(output, "chains", "1", "field_alignment.dat"))
+    return mean(values), std(values) / sqrt(length(values)), length(values)
+end
+
+field_magnitudes = [0.0, 0.5, 1.0, 3.0, 10.0]
+results = [run_one_field(h) for h in field_magnitudes]
+
+summary_path = joinpath(@__DIR__, "alignment-vs-field.csv")
+open(summary_path, "w") do io
+    println(io, "h,mean_alignment,standard_error,n_samples")
+    for (h, (mean_alignment, standard_error, n_samples)) in zip(field_magnitudes, results)
+        println(io, "$h,$mean_alignment,$standard_error,$n_samples")
+        println("h=$h: <hhat dot n> = $mean_alignment +/- $standard_error")
+    end
+end
+println("Wrote $summary_path")

notes/particlesmc_architecture_moves_energy.tex

@@ -249,9 +249,50 @@ \subsection{Triangular Plane-Normal Convention}
 or coincident defining coordinates are rejected because they do not define a
 unique normal.
 
-The external-field Hamiltonian and its contribution to trial-move energy
-changes are not yet implemented; the current code establishes only mass storage
-and the orientation geometry on which that later term will depend.
+\subsection{Aligning External Field}
+
+The aligning field is represented by \texttt{AligningField(h, definition)} and
+adds the one-body molecular energy
+
+\begin{equation}
+H_{\mathrm{field}}
+=
+-\sum_{m=1}^{N_{\mathrm{mol}}}\mathbf h\cdot\mathbf n_m.
+\end{equation}
+
+It is intentionally separate from pair-interaction models. During system
+initialization, the pair contribution is still divided by two, while each field
+term is added exactly once. A system without a field uses
+\texttt{NoExternalField()}, preserving previous behavior.
+
+An aligning field can be specified in TOML as:
+
+\begin{lstlisting}
+[system.external_field]
+type = "align"
+h = [0.0, 0.0, 0.7]
+orientation = "plane_normal"
+atoms = [1, 2, 3]
+\end{lstlisting}
+
+For the center-to-atom convention, use
+\texttt{orientation = "center\_to\_atom"} and \texttt{atom = 1}.
+
+For an ordinary particle displacement, only the orientation of the displaced
+particle's molecule is recomputed. The additional trial energy is
+
+\begin{equation}
+\Delta U_h
+=
+-\mathbf h\cdot
+\left(\mathbf n_m^{\mathrm{new}}-\mathbf n_m^{\mathrm{old}}\right).
+\end{equation}
+
+The implementation dispatches field scoring on the action type. Only
+\texttt{Displacement} is field-aware at present. A future geometric flip move
+can provide its own action-specific method because its affected orientation and
+energy accounting need not be identical to a one-particle displacement. The
+current \texttt{MoleculeFlip}, which swaps species labels, is unchanged.
 
 \section{Interaction Models and Hamiltonians}
 

src/IO/IO.jl

@@ -1,7 +1,7 @@
 module IO
 
 using ..ParticlesMC: Particles, Atoms, Molecules, System
-using ..ParticlesMC: fold_back, volume_sphere
+using ..ParticlesMC: fold_back, volume_sphere, external_field_from_config
 using ..ParticlesMC: EmptyList, LinkedList, CellList, VerletList
 using ..ParticlesMC: Model, GeneralKG, JBB, BHHP, SoftSpheres, KobAndersen, Trimer, LennardJones
 using Arianna
@@ -322,7 +322,8 @@ function load_chains(init_path; args=Dict(), filename="", verbose=false)
         initial_molecule_array = broadcast_dict(config_dict, :molecule)
         initial_bond_array = broadcast_dict(config_dict, :bond)
         masses = get(args, "masses", nothing)
-        chains = [System(initial_position_array[k], initial_species_array[k], initial_molecule_array[k], initial_density_array[k], initial_temperature_array[k], model_matrix, initial_bond_array[k], masses=masses, list_type=list_type, list_parameters=list_parameters) for k in eachindex(initial_position_array)]
+        external_field = external_field_from_config(get(args, "external_field", nothing))
+        chains = [System(initial_position_array[k], initial_species_array[k], initial_molecule_array[k], initial_density_array[k], initial_temperature_array[k], model_matrix, initial_bond_array[k], masses=masses, external_field=external_field, list_type=list_type, list_parameters=list_parameters) for k in eachindex(initial_position_array)]
     else
         chains = [System(initial_position_array[k], initial_species_array[k], initial_density_array[k], initial_temperature_array[k], model_matrix, list_type=list_type, list_parameters=list_parameters) for k in eachindex(initial_position_array)]
     end

src/ParticlesMC.jl

@@ -17,8 +17,9 @@ abstract type Particles <: AriannaSystem end
 include("utils.jl")
 include("neighbours.jl")
 include("models.jl")
-include("molecules.jl")
 include("orientation.jl")
+include("external_fields.jl")
+include("molecules.jl")
 include("atoms.jl")
 include("moves.jl")
 include("rotation.jl")
@@ -139,6 +140,8 @@ export NeighbourList, LinkedList, CellList, EmptyList, VerletList
 export Atoms, Molecules
 export OrientationDefinition, CenterToAtomOrientation, PlaneNormalOrientation
 export orientation, molecule_center_of_mass
+export ExternalField, NoExternalField, AligningField
+export field_energy, total_field_energy
 export Displacement, DiscreteSwap, MoleculeFlip
 export fold_back, System
 export SimpleGaussian, DoubleUniform, EnergyBias
@@ -180,6 +183,7 @@ ParticlesMC implemented in Comonicon.
     list_parameters = get(system, "list_parameters", nothing)  # optional field
     bonds = get(system, "bonds", nothing)
     masses = get(system, "masses", nothing)
+    external_field = external_field_from_config(get(system, "external_field", nothing))
 
     # Extract simulation parameters
     sim = params["simulation"]
@@ -202,6 +206,7 @@ ParticlesMC implemented in Comonicon.
             "list_parameters" => list_parameters,
             "bonds" => bonds,
             "masses" => masses,
+            "external_field" => external_field,
         ),
         filename=filename,
         )
@@ -213,6 +218,7 @@ ParticlesMC implemented in Comonicon.
             "list_type" => list_type,
             "list_parameters" => list_parameters,
             "masses" => masses,
+            "external_field" => external_field,
         ),
         filename=filename,
         )

src/external_fields.jl

@@ -0,0 +1,116 @@
+"""
+External fields acting on molecular observables.
+
+Field energies are kept separate from pair-interaction energies because a
+molecular field contribution is counted once per molecule, whereas the sum of
+per-particle pair energies is divided by two.
+"""
+
+using LinearAlgebra
+using StaticArrays
+
+"""Abstract type for an external-field contribution to the Hamiltonian."""
+abstract type ExternalField end
+
+"""
+Marker used when a molecular system has no external field.
+
+Using a concrete no-op field instead of `nothing` keeps field-energy calls
+branch-free and preserves the behavior of systems constructed without a field.
+"""
+struct NoExternalField <: ExternalField end
+
+"""
+    AligningField(h, orientation_definition)
+
+An aligning field whose contribution for molecule `m` is
+
+```math
+U_m = -\\mathbf{h} \\cdot \\mathbf{n}_m,
+```
+
+where `n_m` is the unit vector computed using `orientation_definition`. This linear coupling is polar: reversing `n_m` changes the
+sign of the field energy.
+"""
+struct AligningField{H<:SVector,O<:OrientationDefinition} <: ExternalField
+    h::H
+    orientation_definition::O
+end
+
+function AligningField(h::AbstractVector{<:Real},
+                       orientation_definition::OrientationDefinition)
+    isempty(h) && throw(ArgumentError("the external-field vector cannot be empty"))
+    all(isfinite, h) || throw(ArgumentError("external-field components must be finite"))
+    h_float = float.(h)
+    return AligningField(SVector{length(h_float)}(h_float), orientation_definition)
+end
+
+"""
+    field_energy(field, system, molecule_id)
+
+Return the external-field contribution of one molecule. Move implementations
+use this local quantity to compute a field-energy difference without summing
+over unaffected molecules.
+"""
+field_energy(::NoExternalField, system, molecule_id::Int) = zero(typeof(system.temperature))
+
+function field_energy(field::AligningField, system, molecule_id::Int)
+    n = orientation(field.orientation_definition, system, molecule_id)
+    return -dot(field.h, n)
+end
+
+"""
+    total_field_energy(field, system)
+
+Return the full one-body field contribution
+
+```math
+U_{\\mathrm{field}} = \\sum_{m=1}^{N_{\\mathrm{mol}}} U_m.
+```
+
+Each molecular term is counted exactly once. In particular, this result must
+not be included in the per-particle pair-energy sum that is later halved.
+"""
+function total_field_energy(field::ExternalField, system)
+    return mapreduce(
+        molecule_id -> field_energy(field, system, molecule_id),
+        +,
+        Base.OneTo(system.Nmol);
+        init=zero(typeof(system.temperature)),
+    )
+end
+
+"""Return the field-vector dimension, or `nothing` for a dimensionless no-op field."""
+field_dimension(::NoExternalField) = nothing
+field_dimension(field::AligningField) = length(field.h)
+
+"""Convert an absent or already constructed field specification to a field object."""
+external_field_from_config(::Nothing) = NoExternalField()
+external_field_from_config(field::ExternalField) = field
+
+"""
+Build an external field from a TOML-compatible dictionary. This is the boundary
+between user-facing configuration values and the typed field/orientation
+objects used by the simulation.
+
+Supported orientations are `"plane_normal"` with `atoms = [a,b,c]` and
+`"center_to_atom"` with `atom = a`.
+"""
+function external_field_from_config(config::AbstractDict)
+    field_type = get(config, "type", "align")
+    field_type == "align" || throw(ArgumentError("unsupported external field type: $field_type"))
+    haskey(config, "h") || throw(ArgumentError("an aligning field requires an `h` vector"))
+
+    orientation_name = get(config, "orientation", "plane_normal")
+    definition = if orientation_name == "plane_normal"
+        atoms = get(config, "atoms", [1, 2, 3])
+        length(atoms) == 3 || throw(ArgumentError("plane-normal orientation requires three atoms"))
+        PlaneNormalOrientation(atoms...)
+    elseif orientation_name == "center_to_atom"
+        CenterToAtomOrientation(get(config, "atom", 1))
+    else
+        throw(ArgumentError("unsupported orientation definition: $orientation_name"))
+    end
+
+    return AligningField(config["h"], definition)
+end