Frenkel-Ladd

[joaander]

Implement the springs necessary to perform Frenkel-Ladd (or Einstein molecule) free energy integration. Users would greatly appreciate a full working example.

[tokiro16]

Hi, I am a heavy HOOMD-blue user. Recently I’m trying to give a possible working example on calculating hard-sphere free energies using the Einstein Molecule method. I've managed to integrate this into hoomd-workflow by customizing the siteproperties/bodyproperties, implementing a custom " External Hamiltonian" and a custom “localtrial”.

It’s working, numerical result looks good. But I’m curious if there are better Implementation which is more elegent. Would love to get some eyes on this and hear your thoughts on it.

---

For tested system ($\rho^*=1.04086$), extrapolated result given by Vega and Noya(2007) is 4.9590(2).
Polyfit: F/NkT = 248.192438 * (1/N)^2 + -8.834959 * (1/N) + 4.959137

---

But it seems difficult to upload the code directly to a GitHub issue comment. How can I contribute?

[janbridley]

You should be able to place the code in a code block, which can be started and ended by sets of 3 backticks (```). Hope this helps, and I'm excited the code is working well for you!

[tokiro16]

Thank you. I hope this lengthy comment won't cause any trouble.

Here is the project directory.

• src/ - source code
  • bin/actions.rs - No change from hoomd-workflow-trunk
  • model/lattice_point.rs - lattice definition (just FCC for this example)
  • lib.rs - just change model name
  • model.rs - model definition
  • state_point.rs - statepoint definition
• cmd.ipynb - command notebook
• post_analysis.py - data post analysis
• populate_workspace.py - create signac jobs
• other things from hoomd-workflow-trunk

Here are the main pieces of code.

src/model/lattice_point.rs

// lattice_point.rs

pub fn lattice_point_fcc_cube(
    replicate: [u32; 3], 
    rho_star: f64
) -> (Vec<[f64; 3]>, [f64; 3]) { // return (lattice points, box length)
    for &dim in &replicate {
        assert!(dim > 0, "Replicate must be positive. ");
    }

    let a = (4.0 / rho_star).powf(1.0 / 3.0);

    let box_length = [
        replicate[0] as f64 * a,
        replicate[1] as f64 * a,
        replicate[2] as f64 * a,
    ];

    let basis = [
        [0.0, 0.0, 0.0],
        [0.5, 0.5, 0.0],
        [0.5, 0.0, 0.5],
        [0.0, 0.5, 0.5],
    ];

    let mut points = Vec::with_capacity((replicate[0] * replicate[1] * replicate[2] * 4) as usize);

    for ix in 0..replicate[0] {
        for iy in 0..replicate[1] {
            for iz in 0..replicate[2] {
                for b in &basis {
                    let x = (ix as f64 + b[0]) * a;
                    let y = (iy as f64 + b[1]) * a;
                    let z = (iz as f64 + b[2]) * a;
                    points.push([x, y, z]);
                }
            }
        }
    }

    if !points.is_empty() {
        let num_points = points.len() as f64;
        
        let mut com = [0.0, 0.0, 0.0];
        for p in &points {
            com[0] += p[0];
            com[1] += p[1];
            com[2] += p[2];
        }
        com[0] /= num_points;
        com[1] /= num_points;
        com[2] /= num_points;

        for p in &mut points {
            p[0] -= com[0];
            p[1] -= com[1];
            p[2] -= com[2];
        }
    }
    (points, box_length)
}

src/model.rs

// #![allow(unused_imports, nonstandard_style)]

use rand::{Rng, random};
use hoomd_geometry::shape::Cuboid;
use hoomd_interaction::{
    MaximumInteractionRange, PairwiseCutoff, External, SiteEnergy, TotalEnergy, DeltaEnergyOne,
    pairwise::HardSphere ,
};
use hoomd_gsd::hoomd::{Dimensions, HoomdGsdFile};
use hoomd_mc::{LocalTrial, Sweep, Trial, Adjust, tune_by_scaling, Count};
use hoomd_microstate::{Microstate, Body, Transform, AppendMicrostate,  SiteKey, boundary::Periodic, property::Position};
use hoomd_simulation::{Simulation, macrostate::Isothermal};
use hoomd_spatial::VecCell;
use hoomd_vector::{Cartesian, Metric};
use hoomd_utility::valid::{OpenUnitIntervalNumber, PositiveReal};
use log::debug;
use serde::{Deserialize, Serialize};

use super::StatePoint;
mod lattice_point;

const RELAX_STEPS: u64 = 10_000;
const TARGET_ACCEPT_RATIO: f64 = 0.3;

// Definitions of types
type PositionVector = Cartesian<3>;
#[derive(Clone, Copy, Default, Position, Serialize, Deserialize)]
pub struct BodyProperties {position: PositionVector, index: usize}
#[derive(Clone, Copy, Default, Position, Serialize, Deserialize)]
pub struct SiteProperties {position: PositionVector, index: usize}
type Boundary = Periodic<Cuboid>;

impl Transform<SiteProperties> for BodyProperties {
    fn transform(&self, site_properties: &SiteProperties) -> SiteProperties {
        SiteProperties {
            position: self.position + site_properties.position,
            ..*site_properties
        }
    }
}

// State Machine
#[derive(Serialize, Deserialize)]
pub enum Phase {
    Relax,
    Equilibrate
}

// Definition of trial move and tuning (fix one particle). Using tune_by_scaling in Relax. 
#[derive(Serialize, Deserialize)]
pub struct FixOneParticleTranslate { 
    d: f64 
}
impl LocalTrial<BodyProperties> for FixOneParticleTranslate {
    fn propose<R:Rng>(
        &self,
        _rng: &mut R,
        body_properties: BodyProperties,
    ) -> BodyProperties {
        
        let mut trial = body_properties;
        if trial.index != 0 {
            let r: f64 = random();
            trial.position.coordinates[0] += self.d * (r * 2.0 - 1.0);
            let r: f64 = random();
            trial.position.coordinates[1] += self.d * (r * 2.0 - 1.0);
            let r: f64 = random();
            trial.position.coordinates[2] += self.d * (r * 2.0 - 1.0);
        };
        trial
    }
}
impl Adjust for FixOneParticleTranslate {
    fn adjust(
        &mut self,
        factor: PositiveReal
    ) {
        self.d *= factor.get();
    }
}

// Definition of 3d EinsteinCrystal potential field.
#[derive(Serialize, Deserialize)] 
pub struct EinsteinCrystalField {
    lambda: f64,
    reference_position: Vec<[f64; 3]>
}
impl SiteEnergy<SiteProperties> for EinsteinCrystalField {
    fn site_energy(
        &self,
        site_properties: &SiteProperties
    ) -> f64 {
        let index = &site_properties.index;
        let particle_pos = &site_properties.position;
        let r2 = PositionVector::from(self.reference_position[*index]).distance_squared(particle_pos);
        self.lambda * r2
    }
}

// Definition of Total Hamiltonian.
#[derive(TotalEnergy, DeltaEnergyOne, MaximumInteractionRange, Serialize, Deserialize)]
pub struct Hamiltonian {
    external: External<EinsteinCrystalField>,
    pairwise_cutoff: PairwiseCutoff<HardSphere>,
}


// The Frenkel-Ladd simulation model for HardSphere. Using Einstein Molecule method. 
// struct definition
#[derive(Serialize, Deserialize)]
pub struct FrenkelLaddHardSphereModel {
    pub macrostate: Isothermal,
    pub hamiltonian: Hamiltonian ,
    pub sweep: Sweep<FixOneParticleTranslate>,
    pub microstate: Microstate<BodyProperties, SiteProperties, VecCell<SiteKey, 3>, Boundary>,
    pub phase: Phase,
    pub translate_count: Count,
    pub relax_step: u64,
}

// Simulation trait.
impl Simulation for FrenkelLaddHardSphereModel {
    #[inline]
    fn advance(&mut self) -> anyhow::Result<()> {
        if self.microstate.step().is_multiple_of(300) {
            self.microstate.sort_sites();
        }

        match self.phase {
            Phase::Relax => self.relax(),
            Phase::Equilibrate => self.equilibrate(),
        }

        self.microstate.increment_step();

        Ok(())
    }

    #[inline]
    fn step(&self) -> u64 {
        self.microstate.step()
    }
}

// Model definition
impl FrenkelLaddHardSphereModel {
    pub fn new(state_point: StatePoint) -> anyhow::Result<Self> {
        // macrostate
        let macrostate = Isothermal { temperature: 1.0 }; // HPMC set kT = 1

        // reference_position
        let reference_position: Vec<[f64;3]>;
        let box_length:[f64;3];
        println!(
            "Generating {} lattice ({}x{}x{}) with rho_star={}",
            state_point.lattice,
            state_point.replicate[0],
            state_point.replicate[1],
            state_point.replicate[2],
            state_point.rho_star
        );
        match state_point.lattice {
            0 => (reference_position, box_length) = lattice_point::lattice_point_fcc_cube(state_point.replicate, state_point.rho_star),
            _ => panic!("未实现的晶格类型标识: {}", state_point.lattice),
        };

        // hamiltonian
        let external = External(EinsteinCrystalField {lambda: state_point.lambda, reference_position: reference_position.clone()});
        let pairwise = PairwiseCutoff( HardSphere { diameter: 1.0 } );
        let hamiltonian = Hamiltonian {external, pairwise_cutoff: pairwise};

        // sweep
        let trial = FixOneParticleTranslate {d: 0.07};
        let sweep = Sweep(trial);

        // microstate
        let cuboid = Cuboid { edge_lengths: [box_length[0].try_into()?, box_length[1].try_into()?, box_length[2].try_into()?] };
        let periodic_cuboid = Periodic::new(hamiltonian.maximum_interaction_range(), cuboid)?;
        let vec_cell = VecCell::builder()
            .nominal_search_radius(hamiltonian.maximum_interaction_range().try_into()?)
            .build();
        let mut microstate = Microstate::builder()
            .seed(state_point.seed)
            .boundary(periodic_cuboid)
            .spatial_data(vec_cell)
            .try_build()?;

        // place particles on lattice
        for (index, pos) in reference_position.iter().enumerate() {
            let body_properties = BodyProperties {
                position: PositionVector::from(*pos),
                index: index,
            };
            let site_properties = SiteProperties {
                position: PositionVector::from([0.0,0.0,0.0]),
                index: index
            };
            let body = Body {
                properties: body_properties,
                sites: vec![site_properties],
            };
            microstate.add_body(body)?;
        };


        Ok(FrenkelLaddHardSphereModel {
            macrostate,
            hamiltonian,
            sweep,
            microstate,
            phase: Phase::Relax,
            translate_count: Count::default(),
            relax_step: 0,
        })
    }

    fn relax(&mut self) {
        self.translate_count +=
            self.sweep
                .apply(&mut self.microstate, &self.hamiltonian, &self.macrostate);

        self.relax_step += 1;
        if self.relax_step.is_multiple_of(300) {
            let v = OpenUnitIntervalNumber::try_from(TARGET_ACCEPT_RATIO).expect("TARGET_ACCEPT_RATIO must be between 0 and 1");
            tune_by_scaling(&mut self.sweep.0, v , &self.translate_count);
        };
        
        if self.relax_step == RELAX_STEPS {
            debug!(
                "Relax phase complete at step {}, retuning trial moves.",
                self.microstate.step()
            );
            self.phase = Phase::Equilibrate;
        };
    }

    fn equilibrate(&mut self) {
        self.translate_count +=
            self.sweep
                .apply(&mut self.microstate, &self.hamiltonian, &self.macrostate);
    }

    pub fn clear_move_counts(&mut self) {
        self.translate_count = Count::default()
    }
}

// implement AppendMicrostate for HoomdGsdFile
impl<X> AppendMicrostate<BodyProperties, SiteProperties, X, Boundary>
    for HoomdGsdFile
{
    #[inline]
    fn append_microstate(
        &mut self,
        microstate: &Microstate<
            BodyProperties,
            SiteProperties,
            X,
            Boundary,
        >,
    ) -> Result<hoomd_gsd::hoomd::Frame<'_>, hoomd_gsd::hoomd::AppendError>
    {
        self.append_frame(microstate.step())?
            .configuration_box(microstate.boundary().shape().to_gsd_box())?
            .configuration_dimensions(Dimensions::Three)?
            .particles_position(
                microstate
                    .iter_sites_tag_order()
                    .map(|s| s.properties.position),
            )
    }
}

src/state_points.rs

use serde::{Deserialize, Serialize};

/// Model parameters that describe a single simulation.
#[derive(Serialize, Deserialize)]
pub struct StatePoint {
    pub lattice: u32, // lattice code. 0 for FCC.
    pub replicate: [u32;3], // replicate for system. [nx, ny, nz]
    pub rho_star: f64, // sigma**3 * N / V
    pub lambda:f64, // spring constant
    pub seed: u32, // simulate random seed
}

populate_workspace.py

"""Populate the workspace."""

import signac
import numpy as np

project = signac.get_project()

def lattice_code(lattice: str):
    # convert lattice_name to lattice_index for statepoint
    if lattice == 'fcc_cube':
        return 0
replicates = [
    [3,3,3],
    [4,4,4],
    [5,5,5],
    [6,6,6],
    [7,7,7]
    ]
lamdas = np.logspace(-2.5,np.log10(3000),60)
for replicate in replicates:
    for lamda in lamdas:
        job = project.open_job(
            {
                "lattice": lattice_code("fcc_cube"),
                "replicate": replicate,
                "rho_star": 1.04086,
                "lambda": lamda,
                "seed": 12254
            }
        ).init()

I have question about rng:&mut Rng in LocalTrial. How should I get a random float number from it? Current implementation used rand::random() but not rng. If you see any ways to improve this code, please let me know!

[janbridley]

In rust, Rng implementors are methods of generating random bits (e.g., via next_u64). To sample a particular distribution - in this case uniform floats - you first construct that distribution object, then call the distribution.sample(rng). I've linked a relevant piece of code from the internal translation trial moves. Hope this helps!

use rand::{
    Rng,
    distr::{Distribution, Uniform},
};
use hoomd_vector::{Cartesian, distribution::Ball};


// hoomd-mc/src/translate/cartesian.rs::L40
fn propose<R: Rng>(&self, rng: &mut R, body_properties: B) -> B {
   let mut trial = body_properties;

   let ball = Ball {
       radius: self.maximum_distance,
   };

   let delta_r = ball.sample(rng);
   *trial.position_mut() += delta_r;

   trial
}

I'll ping Josh about this PR when he's back in office and will see if he has any further thoughts on the code!

[joaander]

It’s working, numerical result looks good. But I’m curious if there are better Implementation which is more elegent.

Your implementation looks great! I don't see much room to make it more elegant with hoomd-rs 1.0. You could save a couple lines of code by reusing Translate via composition. In a future hoomd-rs release, I could add #derive(Adjust) to save a couple more lines of code:

#[derive(Adjust, Serialize, Deserialize)]
pub struct FixOneParticleTranslate(pub LocalTrial<PositionVector>);

impl LocalTrial<BodyProperties> for FixOneParticleTranslate {
    fn propose<R:Rng>(
        &self,
        rng: &mut R,
        body_properties: BodyProperties,
    ) -> BodyProperties {
        if trial.index != 0 {
            self.0.propose(rng, body_properties);
        } else {
            body_properties
        }
    }
}

I could instead add an apply_filtered method to Sweep so that there is no need for a custom impl LocalTrial at all:

let translate = Translate::with_maximum_distance(...);
let sweep = Sweep(translate);

sweep.apply_filtered(&mut self.microstate, &self.hamiltonian, &self.macrostate, |b| b.index != 0);

I haven't decided which is the better option: apply_filtered or a custom local trial. I don't want to provide a method that is easy to use incorrectly. I also need to figure out a solution so that move sizes are tuned appropriately when moves are filtered in general.

The form of the potential drives the need for a custom site type which by extension requires the corresponding impl AppendMicrostate and impl Transform<SiteProperties>. There really isn't anything that can be simplified there.

[tokiro16]

Thank you very much for your reply! It would be great if adjust supported derive. I'm really looking forward to the new version of hoomd-rs!

By the way, I'm very interested in what you mentioned about "incorrect use of filter" and the "appropriate tune". Could you give an example of what the improper use looks like?

[joaander]

An incorrect use of filter would be to apply trial moves with apply_filtered but tune those trial moves with the unfiltered tune_default. The user would obtain an unexpected acceptance ratio.

tune_default also (currently) tunes inappropriately when the LocalTrial skips trials on some bodies. For example, if a LocalTrial filtered out half the system, then tune_default would always read a minimum acceptance ratio of 50% and it would always scale the move size down.

[joaander]

In hoomd-rs 1.1, you can now use Sweep::apply_with_filter to skip trial moves for some particles. In the sample code discussed here, that removes the need for a custom trial move and for custom body properties. The new Seeded Self-Assembly tutorial demonstrates apply_with_filter.