diff --git a/examples/ortho-terphenyl/5-single-molecule-field-test/.gitignore b/examples/ortho-terphenyl/5-single-molecule-field-test/.gitignore new file mode 100644 index 0000000..ea1472e --- /dev/null +++ b/examples/ortho-terphenyl/5-single-molecule-field-test/.gitignore @@ -0,0 +1 @@ +output/ diff --git a/examples/ortho-terphenyl/5-single-molecule-field-test/README.md b/examples/ortho-terphenyl/5-single-molecule-field-test/README.md new file mode 100644 index 0000000..910bcb0 --- /dev/null +++ b/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. diff --git a/examples/ortho-terphenyl/5-single-molecule-field-test/alignment-vs-field.csv b/examples/ortho-terphenyl/5-single-molecule-field-test/alignment-vs-field.csv new file mode 100644 index 0000000..34ff7b6 --- /dev/null +++ b/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 diff --git a/examples/ortho-terphenyl/5-single-molecule-field-test/run-field-scan.jl b/examples/ortho-terphenyl/5-single-molecule-field-test/run-field-scan.jl new file mode 100644 index 0000000..941cd50 --- /dev/null +++ b/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: = $mean_alignment +/- $standard_error") + end +end +println("Wrote $summary_path") diff --git a/notes/particlesmc_architecture_moves_energy.tex b/notes/particlesmc_architecture_moves_energy.tex index 3df1407..290ff1b 100644 --- a/notes/particlesmc_architecture_moves_energy.tex +++ b/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} diff --git a/src/IO/IO.jl b/src/IO/IO.jl index 519663c..ed78bbd 100644 --- a/src/IO/IO.jl +++ b/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 diff --git a/src/ParticlesMC.jl b/src/ParticlesMC.jl index bf021b9..0507acd 100644 --- a/src/ParticlesMC.jl +++ b/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, ) diff --git a/src/external_fields.jl b/src/external_fields.jl new file mode 100644 index 0000000..7ed3893 --- /dev/null +++ b/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 diff --git a/src/molecules.jl b/src/molecules.jl index 5eea846..1375be4 100644 --- a/src/molecules.jl +++ b/src/molecules.jl @@ -2,7 +2,7 @@ Molecules - type managing a collection of molecules (chains) of particles. Description -`Molecules{D,VS,C,T,SM} <: Particles` represents a particle system organized into +`Molecules{D,VS,C,T,SM,F} <: Particles` represents a particle system organized into molecules (or chains). It stores particle positions, molecule membership, molecule lengths and start indices, thermodynamic properties, interaction models, and the neighbour list structure used for energy and neighbour queries. @@ -10,6 +10,7 @@ models, and the neighbour list structure used for energy and neighbour queries. Fields - `position::Vector{SVector{D,T}}`: positions of particles. - `mass::Vector{T}`: mass of each particle. +- `external_field::F`: external one-body field acting on molecular orientation. - `Φ::Vector{Vector{SVector{3,T}}}`: rotation vector of molecules per threshold index k. - `species::VS`: species/type of each particle. - `molecule::VS`: molecule identifier for each particle. @@ -23,9 +24,10 @@ Fields - `neighbour_list::C`: neighbour list structure (CellList/LinkedList/EmptyList). - `bonds::Vector{Vector{Int}}`: lists of bonded neighbours for each particle. """ -struct Molecules{D, VS<:AbstractVector, C<:NeighbourList, T<:AbstractFloat, SM<:AbstractArray} <: Particles +struct Molecules{D, VS<:AbstractVector, C<:NeighbourList, T<:AbstractFloat, SM<:AbstractArray, F<:ExternalField} <: Particles position::Vector{SVector{D,T}} mass::Vector{T} + external_field::F Φ::Vector{Vector{SVector{3,T}}} # Φ[k][m] = rotation vector for molecule m, threshold k species::VS molecule::VS @@ -61,7 +63,7 @@ struct NonBonded end """ Create a `Molecules` system and initialize neighbour list and energy. -`System(position, species, molecule, density, temperature, model_matrix, bonds; masses=nothing, molecule_species=nothing, list_type=EmptyList)` builds +`System(position, species, molecule, density, temperature, model_matrix, bonds; masses=nothing, external_field=NoExternalField(), molecule_species=nothing, list_type=EmptyList)` builds and returns a `Molecules` instance. It computes `start_mol`/`length_mol`, constructs the neighbour list of type `list_type`, builds it, and computes the initial total energy (with a check for Inf/NaN). @@ -77,7 +79,7 @@ Arguments Returns - `Molecules` instance with neighbour list built and `energy[1]` set. """ -function System(position, species, molecule, density::T, temperature::T, model_matrix, bonds; masses=nothing, molecule_species=nothing, list_type=EmptyList, list_parameters=nothing) where {T<:AbstractFloat} +function System(position, species, molecule, density::T, temperature::T, model_matrix, bonds; masses=nothing, external_field::ExternalField=NoExternalField(), molecule_species=nothing, list_type=EmptyList, list_parameters=nothing) where {T<:AbstractFloat} @assert length(position) == length(species) N = length(position) mass = isnothing(masses) ? ones(T, N) : T.(collect(masses)) @@ -89,14 +91,17 @@ function System(position, species, molecule, density::T, temperature::T, model_m start_mol, length_mol = get_first_and_counts(molecule) molecule_species = something(molecule_species, ones(Int, N)) d = length(Array(position)[1]) + field_d = field_dimension(external_field) + isnothing(field_d) || field_d == d || + throw(DimensionMismatch("external field has dimension $field_d, but the system has dimension $d")) box = @SVector fill(T((N / density)^(1 / d)), d) energy = zeros(T, 1) maxcut = maximum([model.rcut for model in model_matrix]) neighbour_list = list_type(box, maxcut, N; list_parameters=list_parameters) - system = Molecules(position, mass, Φ, species, molecule, molecule_species, start_mol, length_mol, density, temperature, energy, model_matrix, d, N, Nmol, box, neighbour_list, bonds) + system = Molecules(position, mass, external_field, Φ, species, molecule, molecule_species, start_mol, length_mol, density, temperature, energy, model_matrix, d, N, Nmol, box, neighbour_list, bonds) build_neighbour_list!(system) local_energy = [compute_energy_particle(system, i, neighbour_list) for i in eachindex(position)] - energy = sum(local_energy) / 2 + energy = sum(local_energy) / 2 + total_field_energy(external_field, system) if isinf(energy) || isnan(energy) error("Initial configuration has infinite or NaN energy.") end diff --git a/src/moves.jl b/src/moves.jl index 71d35d0..2fdfbbe 100644 --- a/src/moves.jl +++ b/src/moves.jl @@ -66,6 +66,34 @@ function perform_action!(system::Particles, action::Displacement) return e₁, e₂ end +""" +Return the external-field energy associated with the molecule affected by a +particle displacement. Future move types, including a geometric molecule flip, +can provide their own action-specific method without changing the field model. +""" +function action_field_energy(system::Molecules, action::Displacement) + molecule_id = system.molecule[action.i] + return field_energy(system.external_field, system, molecule_id) +end + +""" +Perform a molecular particle displacement, including the field-energy change of +the displaced particle's molecule exactly once. +""" +function perform_action!(system::Molecules, action::Displacement) + neighbour_list = get_neighbour_list(system) + e₁ = compute_energy_particle(system, action.i, neighbour_list) + + action_field_energy(system, action) + + update_position!(system, action) + update_neighbour_list!(system, action.i, neighbour_list) + + e₂ = compute_energy_particle(system, action.i, neighbour_list) + + action_field_energy(system, action) + action.δe = e₂ - e₁ + return e₁, e₂ +end + """ Revert a previously applied `Displacement` action. diff --git a/test/runtests.jl b/test/runtests.jl index c5f4124..40f321a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -5,8 +5,90 @@ using StaticArrays using Distributions using ComponentArrays using DelimitedFiles +using LinearAlgebra using Pkg +@testset "Molecular orientation" begin + box = SVector(10.0, 10.0, 10.0) + triangle = [ + SVector(1.0, 1.0, 1.0), + SVector(2.0, 1.0, 1.0), + SVector(1.0, 2.0, 1.0), + ] + equal_masses = ones(3) + unequal_masses = [1.0, 2.0, 1.0] + + @test molecule_center_of_mass(triangle, equal_masses, box) ≈ + SVector(4 / 3, 4 / 3, 1.0) + @test molecule_center_of_mass(triangle, unequal_masses, box) ≈ + SVector(1.5, 1.25, 1.0) + @test orientation(CenterToAtomOrientation(1), triangle, equal_masses, box) ≈ + normalize(SVector(-1.0, -1.0, 0.0)) + @test orientation(CenterToAtomOrientation(1), triangle, unequal_masses, box) ≈ + normalize(SVector(-0.5, -0.25, 0.0)) + @test orientation(PlaneNormalOrientation(), triangle, box) ≈ + SVector(0.0, 0.0, 1.0) + @test orientation(PlaneNormalOrientation(1, 3, 2), triangle, box) ≈ + SVector(0.0, 0.0, -1.0) + + wrapped_triangle = [ + SVector(9.5, 1.0, 1.0), + SVector(0.5, 1.0, 1.0), + SVector(9.5, 2.0, 1.0), + ] + @test orientation(CenterToAtomOrientation(1), wrapped_triangle, equal_masses, box) ≈ + normalize(SVector(-1.0, -1.0, 0.0)) + @test orientation(PlaneNormalOrientation(), wrapped_triangle, box) ≈ + SVector(0.0, 0.0, 1.0) + + @test_throws ArgumentError PlaneNormalOrientation(1, 1, 3) + @test_throws ArgumentError molecule_center_of_mass(triangle, [1.0, 2.0], box) + @test_throws ArgumentError molecule_center_of_mass(triangle, [1.0, 0.0, 1.0], box) + @test_throws ArgumentError orientation( + PlaneNormalOrientation(), + fill(SVector(1.0, 1.0, 1.0), 3), + box, + ) +end + +@testset "Aligning field with displacement" begin + common_args = Dict( + "temperature" => [2.0], + "model" => ["Trimer"], + "list_type" => "LinkedList", + ) + field = AligningField(SVector(0.0, 0.0, 0.7), PlaneNormalOrientation()) + no_field_system = load_chains("molecule.exyz", args=common_args)[1] + field_system = load_chains( + "molecule.exyz", + args=merge(common_args, Dict("external_field" => field)), + )[1] + + @test field_system.energy[1] - no_field_system.energy[1] ≈ + total_field_energy(field, field_system) + + molecule_id = field_system.molecule[1] + old_field_energy = field_energy(field, field_system, molecule_id) + no_field_action = Displacement(1, SVector(0.0, 0.0, 0.02), 0.0) + field_action = Displacement(1, SVector(0.0, 0.0, 0.02), 0.0) + ParticlesMC.perform_action!(no_field_system, no_field_action) + ParticlesMC.perform_action!(field_system, field_action) + new_field_energy = field_energy(field, field_system, molecule_id) + + @test field_action.δe - no_field_action.δe ≈ + new_field_energy - old_field_energy atol=1e-10 + @test field_system.external_field === field + + config_field = ParticlesMC.external_field_from_config(Dict( + "type" => "align", + "h" => [0.0, 0.0, 0.7], + "orientation" => "center_to_atom", + "atom" => 2, + )) + @test config_field isa AligningField + @test config_field.orientation_definition == CenterToAtomOrientation(2) +end + @testset "Running from CLI" begin Pkg.build("ParticlesMC") path_sep = Sys.iswindows() ? ";" : ":" @@ -235,4 +317,3 @@ end run!(simulation) println("OK") end -