Add resize MC move. #58

g-bauer · 2016-11-28T23:29:17Z

This adds a MC move for changing the size of the simulation cell.

Todo issue #59 :
I added a check to compare the smallest inscribed spheres' radius vs the cut off radius rc_max. When the cell gets too small (so that the cut off is too large) it should cancel the simulation and report an error. I'd like to get the (maximum) cut off out of the PairInteractions when invoking new() but it is not clear how to get it at the moment.

@Luthaf we wrote the main part of this together. I created a new branch since the history was messed up in the old branch. Now, you are not reported as author. Can we change this somehow?

Luthaf · 2016-11-29T10:04:09Z

@Luthaf we wrote the main part of this together. I created a new branch since the history was messed up in the old branch. Now, you are not reported as author. Can we change this somehow?

This should be possible, I've never done it. Maybe something using rebase, and then amend the commit to add an author. I don't even know if you can have more than one author for a commit. Anyway, you can leave it as it is, I can live without being cited as author on this commit 😉.

I'll give this PR an eye later today.

g-bauer · 2016-11-29T10:23:02Z

I can live without being cited as author on this commit 😉

So do I. If we find a way to add/change author, I'd be happy.

I'll give this PR an eye later today.

Great! I will write a test case.

Luthaf · 2016-11-29T11:06:15Z

So do I. If we find a way to add/change author, I'd be happy.

I checked that, and git can not associate two authors with the a given commit. You can add me in the commit message (`Implemented with ...); or I can do this before merging.

Luthaf · 2016-11-29T20:22:42Z

src/core/src/sim/mc/moves/resize.rs

+            range: Range::new(-delta, delta),
+            old_system: System::new(),
+            pressure: pressure,
+            rc_max: 100.0, // arbitrarily large value


You could use std::f64::INFINITY here.

I recognized that there has to be a zero. std::f64::INFINITY would lead to the error, since any side length is smaller than infinity.

Luthaf · 2016-11-29T20:24:11Z

src/core/src/sim/mc/moves/resize.rs

+        // cut off distance.
+        // abort simulation when box gets too small
+        if new_cell.lengths_as_vec().iter().any(|&d| 0.5 * d < self.rc_max) {
+            fatal_error!(


I'd argue this should not be a fatal error. It should send an error message saying that the energy will be wrong for now on, but not abort the process, as fatal_error! does.

Well, wrong energies lead to wrong configurations. In the end, the ensemble averages computed are incorrect. We could adjust the cutoff so that it fits into the new box. That is done in some codes. Personally, I think this solution is also not a very good idea, since it essentially changes the force field.

I don't like the idea of changing the cutoff too.

I prefer to leave control to the user, instead of doing things behind him. We can leave a fatal_error, and reconsider if this is being very annoying.

Luthaf · 2016-11-29T20:24:49Z

src/core/src/sim/mc/moves/resize.rs

+        cache.unused();
+        // get system energy before change: this is stored in `self.old_system`
+        let old_energy = self.old_system.potential_energy();
+        // get


Unfinished comment

Thanks for the catch.

Luthaf · 2016-11-29T20:31:13Z

src/core/src/sys/cells.rs

@@ -174,6 +174,23 @@ impl UnitCell {
        return (x, y, z);
    }

+    /// Get a vector of the distances between faces of the unit cell
+    pub fn lengths_as_vec(&self) -> Vec<f64> {


We do not need another function here, the UnitCell::lengths can do this too! You seems to need iteration on the lengths, so why not change the signature to be pub fn lengths(&self) -> [f64; 3]? This way you get iteration, and do not pay for a dynamic allocation.

The current function is only used once, for the exact same reason as here.

Good catch! I will change that. Actually, I tried to get the maximal value out of the tuple but it was kind of messy. I am not used to the fact that there is no max function for arrays and tuples.

g-bauer · 2016-11-30T09:06:36Z

You can add me in the commit message (`Implemented with ...); or I can do this before merging.

Please, be so kind and do this when merging 😄

Luthaf · 2016-11-30T09:12:18Z

This is missing some integration tests before merging.

I would at least have a check that for a Lennard-Jones fluid, the pressure does converge to the requested value, adding a new case here.

Please, be so kind and do this when merging 😄

Yes, will do so!

g-bauer · 2016-11-30T09:14:28Z

I will add test case(s) tonight.

Luthaf · 2016-12-01T13:45:17Z

A practical question: in #58, I added an example case. Somehow, I couldn't manage to make this a #[test]. I compiles, but is not able to find the configuration file. I placed the config in lumol/tests/data/ and the test in lumol/tests. Using cargo test --test my_test --release fails.

Are you reading the file like this ?

let data_dir = Path::new(file!()).parent().unwrap();
let configuration = data_dir.join("data").join("argon.xyz");
let mut system = Trajectory::open(configuration)
                                            .and_then(|mut traj| traj.read())
                                            .unwrap();

Because your error looks like a problem with the current directory of the test process.

g-bauer · 2016-12-01T13:57:25Z

That's exactly what I am doing. I copied this part from the existing test helium.rs (which throws the same error).

Edit: Oh, i see. I have to be in the examples directory to actually run cargo run --example mc-npt_argon --release otherwise it will use my current path.

Luthaf · 2016-12-01T14:14:12Z

Weird. I'll clone this branch and have a look. Can you commit your changes and push them ?

g-bauer · 2016-12-01T15:38:52Z

Changes are commited. Both, example and test compile. Example just works when I call it from inside lumol/examples/. Test does not run.

I did:

$ cd /lumol/tests
$ cargo test --test mc_npt_argon

Luthaf · 2016-12-01T15:47:31Z

Ok, I got it ! You need to run cargo from the lumol folder, not the lumol/tests folder.

Path::new(file!()).parent().unwrap(); will give you the directory containing the current file, so configuration will be tests/data/argon.xyz.

At the same time, the examples do not use the Path::new(file!()).parent().unwrap(); trick, so they must be run from the example directory. I realize this is actually confusing, we might want to change examples to use paths from the root directory too. Or we could change everything to use absolute paths only.

g-bauer · 2016-12-01T15:52:32Z

I would at least have a check that for a Lennard-Jones fluid, the pressure does converge to the requested value, adding a new case here.

Energies and volumes look good. I need to perform some long runs over night to actually check the statistics. Internal pressure seems (no good statistics) to work fine for the case without tail corrections. With tail corrections, there is an issue. I'll check that.

Luthaf · 2016-12-17T13:25:10Z

Can you rebase this on master now? All the fixes should be here !

g-bauer · 2016-12-17T13:26:27Z

I did a lot more work on this. I will rebase and push the changes later today.

g-bauer · 2016-12-19T23:19:12Z

This is still WIP

I changed a lot in the latest iteration. Ran some short simulations, while ~~energy~~ and density seem to be correct, internal pressure is not correct yet.

I changed the struct to contain the new_system instead of the old/current system. This is more in line with the translation and rotation moves, only the apply function will modify the system.
Also, positions of molecules ares scaled such that the intramolecular distances are unchanged.
At the moment, energy computation is inefficient. We can fetch the current (old) energy from cache. The new energy is computed from scratch. Since there is no computation via cache yet, when applying the move, we recompute the whole energy again (cache.unused() -> cache.init).

Edit: I realized I left some stuff in the PR that I included for testing, i.e. the whole particle... access in system is not used. Just ignore it for now, I will delete that.

Luthaf · 2016-12-20T23:31:50Z

src/core/src/sim/mc/moves/resize.rs

+        // system.
+        // TODO: include elecetrostatic interactions
+        self.rc_max = system.interactions()
+            .all_pairs()


This would look better if aligned with the .interactions() =)

Luthaf · 2016-12-20T23:34:40Z

Also, positions of molecules ares scaled such that the intramolecular distances are unchanged.

Nice !

At the moment, energy computation is inefficient. We can fetch the current (old) energy from cache. The new energy is computed from scratch. Since there is noch computation via cache yet, when applying the move, we recompute the whole energy again (cache.unused() -> cache.init).

We can add a EnergyCache::recompute_all(system) -> f64 (bad name) function that recompute all the terms, and call it here. The EnergyCache::recompute_all function can set the right updater closure to skip computing everything twice in the case we accept the move.

g-bauer · 2016-12-21T11:09:22Z

Another take would be to make this more like the translation move, i.e. not to clone the whole system, but to get only the components that actually change. For example

/// Monte-Carlo move that changes the size of the simulation cell
pub struct Resize {
    /// Delta for translation of the box length
    delta: f64,
    /// Sampling range for volume scaling
    range: Range<f64>,
    /// Positions after applying changes to the simulation cell
    new_pos: Vec<Vector3D>,
    /// New simulation cell
    new_cell: UnitCell,
    /// target pressure
    pressure: f64,
    /// largest cut off diameter of `PairPotentials`
    rc_max: f64,
}

We could then create a caching function like

pub fn resize_cell_cost(&mut self, system: &System, new_pos: &[Vector3D], new_cell: UnitCell) -> f64

that only recomputes intermolecular changes. This would work nicely with the update closure.
There are currently two things preventing such an implementation: 1) We'd need to also implement this for electrostatic interactions, 2) applying the move is not as elegant; we'd have to reiterate over all particles in the system and change each position individually.

I will debug the current approach, as it is the most comfortable to work with for now.

g-bauer · 2016-12-26T17:12:36Z

I added a function move_rigid_molecules_cost (maybe change to move_resize_cell_cost?) that takes a (new) system and recomputes only the intermolecular interactions + tail corrections + global interactions from scratch which makes resize more efficient to use. I didn't dive into the caching of coulomb and global - these are recomputed from scratch at the moment, I added a TODO.

I also changed the sign of the tail correction to the virial. I did two quick runs for Lennard Jones particles and compared to NIST. Energies and densities are in very good (as far as I can tell from the short runs) agreement, internal pressure fits external pressure.

Edit: Added a check in Resize::setup which aborts the simulation if the system cell is infinite.

Todo:

I will add tests/examples for TraPPE Ethane and Water to check if molecules/electrostatics also work as intended.

Luthaf · 2016-12-29T22:18:03Z

src/core/src/sys/cache.rs

+
+        let mut new_system = system.clone();
+        // translate the center of mass
+        for pi in &system.molecules()[0] {


This can be system.molecule(0) until we get better API for molecules.

Luthaf · 2016-12-29T22:19:44Z

src/core/src/sys/compute.rs

@@ -201,7 +200,7 @@ impl Compute for Virial {
                let nj = composition[&j] as f64;
                let potentials = system.interactions().pairs(i, j);
                for potential in potentials {
-                    virial -= 2.0 * PI * ni * nj * potential.tail_virial() / volume;
+                    virial += 2.0 * PI * ni * nj * potential.tail_virial() / volume;


This should be in a separated commit (commits should do only one thing), and maybe with an explanation message (the commit message) and associated tests.

It is fine for this one, can you try to do it for your next commits?

You are right. I'll try to keep separate issues in separate commits in the future.

Luthaf · 2016-12-29T22:21:55Z

Energies and densities are in very good (as far as I can tell from the short runs) agreement, internal pressure fits external pressure.

Yay ! \o/

I didn't dive into the caching of coulomb and global - these are recomputed from scratch at the moment, I added a TODO.

This can be a separated PR.

g-bauer · 2017-01-03T22:32:38Z

I added functionalities to wrap molecules into the simulation cell. I use these to shift molecules as a whole, as soon as the center of mass of the molecule leaves the cell (some particles may still be outside the cell, but this way the com is always inside the cell). I used this mainly to debug some weird stuff happening to the bond lengths. In the end, the molecule wrapping is not necessary but it might come in handy in the future, for example when writing the configuration to a file.

I added a property dr_max to Translate which is set via Translate::setup that is used to limit the amplitude for translations. I found that for the ethane system, i.e. a dilute gas phase, the amplitude explodes since virtually all translations are accepted. That leads to numerical problems.

I added an example case for ethane. I found very good agreement for energy and density w.r.t other codes, and internal pressure agrees very well with external pressure. Going to add a case for SPC/E water.

Luthaf · 2017-01-04T14:29:32Z

examples/mc_npt_ethane.rs

+    system.set_cell(UnitCell::cubic(100.0));
+
+    // Add intermolecular interactions
+    // TraPPE parameters for CH3 - note that masses are not correct.


You can iterate over the particles and set the right mass like this:

for particle in &mut system { assert_eq!(particle.name(), "C"); particle.mass = ...; }

Ah I see. Here, it would not make a difference, since we don't output the density.

small changes added largest cutoff computation to setup routine added TODO to skip com computation

…ite.

…ation cell.

…ng for com functions.

…n the system.

…make functions more transparent.

…_amplitude.

Luthaf · 2017-01-18T13:10:18Z

examples/mc_npt_argon.rs

+
+    let temperature = units::from(102.231042155, "K").unwrap();
+    let pressure = units::from(20.1586274874, "bar").unwrap();
+    //let mut mc = MonteCarlo::from_rng(temperature, rng); // This is not working, why?


Cleanup? And maybe open an issue with the error if this is not working =)

Luthaf

Ok, I performed a last review on the whole code, can you do a last bit of cleanup? Then we'll get that merged!

Luthaf · 2017-01-18T13:13:18Z

examples/mc_npt_spce.rs

+        match p.name() {
+            "H" => p.charge = h,
+            "O" => p.charge = o,
+            _ => panic!("Particle {}", p.name()),


Maybe use a better message here ("Unknown particle in charge setting")? And replace p by particle?

Thanks for the comment, I added a better description.

Luthaf · 2017-01-18T13:18:01Z

src/core/src/sys/systems.rs

+    ///
+    /// This function does not check for the particles' positions'
+    /// nearest images. To use this function properly, make sure
+    /// that all particles of the molecule are adjacent.  


I don't understand this warning. Is it in the case where we have some particles of the molecule on a side of the cell and some other particles on the other side? Why not check particles nearest image here?

Is it in the case where we have some particles of the molecule on a side of the cell and some other particles on the other side?

Exactly. To avoid that, we can use wrap_molecule. It get's messy when an input file is used and parts of a molecule are already wrapped.

I would like to address this in a separate PR, if that's ok for you.

Luthaf · 2017-01-18T13:19:33Z

src/core/src/sys/systems.rs

+        // cells, a wrap will return just the same vector and
+        // hence `delta` will be zero; but it is faster.
+        match self.cell().shape() {
+            CellShape::Infinite => (),


This is not needed, as the wrap_vector function already manage the Infinite cell case. The same code ca be used everywhere.

That's true. Removed the match statement.

Luthaf · 2017-01-18T13:19:50Z

src/core/src/sys/systems.rs

+    type Output = [Particle];
+    #[inline]
+    fn index(&self, index: Range<usize>) -> &[Particle] {
+        // Index::index(&self.particles, index)


Should we keep indexing function via Range? We currently don't use it.

If we do not use it, let's remove it for now. It is easy to add back if needed.

Yep, I removed it.

Luthaf · 2017-01-18T13:20:39Z

tests/mc_npt_argon.rs

+
+    let temperature = units::from(120.0, "K").unwrap();
+    let pressure = units::from(12.6, "bar").unwrap();
+    //let mut mc = MonteCarlo::from_rng(temperature, rng); // This is not working, why?


Luthaf · 2017-01-18T13:21:04Z

tests/mc_npt_argon.rs

+
+    // rng: set this here in case we want to change the seed
+    //let mut rng = Box::new(rand::XorShiftRng::new_unseeded());
+    //rng.reseed([2015u32, 42u32, 3u32, 12u32]);


cleanup? This is not adding a lot of information to this test.

Luthaf · 2017-01-18T13:22:16Z

tests/mc_npt_argon.rs

+    // run simulation for at least 500 cycles (so that a volume update is performed)
+    simulation.run(&mut system, 500000);
+
+    // [1] Lotfi et al.: 


Why are we not checking that in assert_approx_eq!

What do you mean? Checking versus energy and density from the paper?

Yes, or move this from the test folder to the examples folder. It does not seems to test anything yet =)

We already have a NPT example for argon. I remove this one. The simulation has to run very long to get comparable results.

Luthaf · 2017-01-18T13:24:20Z

I rebased the whole branch on top of master, and have added our names to the first commit.

g-bauer · 2017-01-18T21:19:25Z

Thanks for the review. Can I make my changes locally and then simply push as usual or do I have to fetch your changes beforehand?

Luthaf · 2017-01-18T21:43:27Z

You have to get the changes first. Use git fetch origin and then git reset --hard origin/resize-cell (assuming github.com/g-bauer/lumol is your origin)

g-bauer · 2017-01-21T16:38:32Z

Changed according to review. @Luthaf Again, thanks for your patience and your time investment.

The examples for SPCE water and ethane are still in the example folder although they currently do not produce correct samples. For ethane, pressure is wrong due to using the number of atoms for the ideal gas pressure which we discussed above and which will be resolved in a separate PR. The same is true for SPCE which is also affected by the bugs that were found and worked on in #75.

I left the current implementation of the center of mass computation - without checking if a molecule is partly wrapped. This would result in a wrong center of mass. I'd address changes of the com in a separate PR.

Luthaf · 2017-01-21T19:01:46Z

Thank you @g-bauer, that was an epic PR! More than 6 month between the initial implementation and the merge, that is the best I have done.

g-bauer · 2017-01-21T19:20:16Z

Thank you @g-bauer, that was an epic PR! More than 6 month between the initial implementation and the merge, that is the best I have done.

Yeah! 🎉 Thanks @Luthaf for mentoring!

g-bauer mentioned this pull request Nov 29, 2016

Implement NPT Monte-Carlo #15

Closed

Luthaf requested changes Nov 29, 2016

View reviewed changes

Luthaf approved these changes Nov 30, 2016

View reviewed changes

g-bauer force-pushed the resize-cell branch from 4d9d068 to e4708bb Compare November 30, 2016 22:53

g-bauer mentioned this pull request Dec 1, 2016

Benchmarking Lumol #62

Open

22 tasks

g-bauer force-pushed the resize-cell branch from 5ee1ddb to 8b9c9a2 Compare December 19, 2016 23:08

Luthaf reviewed Dec 20, 2016

View reviewed changes

Luthaf reviewed Dec 29, 2016

View reviewed changes

Luthaf reviewed Jan 4, 2017

View reviewed changes

g-bauer added 16 commits January 18, 2017 14:07

Added setup to MC propagator

97922c7

small changes added largest cutoff computation to setup routine added TODO to skip com computation

Changed struct to contain new_system.

fe85d06

Some changes to example

72dcb7e

Added caching function for resize move. Changed sign in p_tail.

f12ebfe

Add check in setup to abort simulation when system CellShape is Infin…

c85f9b1

…ite.

Add function that returns a new vector that is wrapped into the simul…

6a198ce

…ation cell.

Add function to wrap a molecule using its center-of-mass. Added warni…

255978a

…ng for com functions.

Fixed bug in virial and stress computation to account for molecules i…

a5ce69d

…n the system.

Molecules are now wrapped when a move is accepted. Added comments to …

efc04cd

…make functions more transparent.

Added comments.

16e9b1a

Add restriction to the amplitude for translations in setup and update…

e4ee605

…_amplitude.

Added example for TraPPE ethane in the NPT ensemble using MC.

8a928b2

Small change to cache due to revision.

c730288

Added example for SPCE water.

6234cb6

Changes due to revision.

f11a3e9

Cleaned up remaining code of removed functions.

fbd15b9

Luthaf force-pushed the resize-cell branch from dbb5cce to fbd15b9 Compare January 18, 2017 13:08

Luthaf reviewed Jan 18, 2017

View reviewed changes

Luthaf requested changes Jan 18, 2017

View reviewed changes

Changes due to review.

f21e2a0

Luthaf merged commit 398b60d into lumol-org:master Jan 21, 2017

g-bauer deleted the resize-cell branch January 21, 2017 19:20

g-bauer mentioned this pull request Jan 21, 2017

Add cut off distance to resize MC move to catch too small volumes #59

Closed

Luthaf mentioned this pull request Sep 18, 2018

Use the molecular definition of virial for rigid molecules #253

Merged

Add resize MC move. #58

Add resize MC move. #58

Conversation

g-bauer commented Nov 28, 2016 • edited Loading

Luthaf commented Nov 29, 2016

g-bauer commented Nov 29, 2016

Luthaf commented Nov 29, 2016

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

g-bauer commented Nov 30, 2016

Luthaf commented Nov 30, 2016

g-bauer commented Nov 30, 2016

Luthaf commented Dec 1, 2016

g-bauer commented Dec 1, 2016 • edited Loading

Luthaf commented Dec 1, 2016

g-bauer commented Dec 1, 2016

Luthaf commented Dec 1, 2016

g-bauer commented Dec 1, 2016

Luthaf commented Dec 17, 2016

g-bauer commented Dec 17, 2016

g-bauer commented Dec 19, 2016 • edited Loading

This is still WIP

Choose a reason for hiding this comment

Luthaf commented Dec 20, 2016

g-bauer commented Dec 21, 2016 • edited Loading

g-bauer commented Dec 26, 2016 • edited Loading

Todo:

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Luthaf commented Dec 29, 2016

g-bauer commented Jan 3, 2017

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Luthaf left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Luthaf commented Jan 18, 2017

g-bauer commented Jan 18, 2017

Luthaf commented Jan 18, 2017

g-bauer commented Jan 21, 2017 • edited Loading

Luthaf commented Jan 21, 2017

g-bauer commented Jan 21, 2017

g-bauer commented Nov 28, 2016 •

edited

Loading

g-bauer commented Dec 1, 2016 •

edited

Loading

g-bauer commented Dec 19, 2016 •

edited

Loading

g-bauer commented Dec 21, 2016 •

edited

Loading

g-bauer commented Dec 26, 2016 •

edited

Loading

g-bauer commented Jan 21, 2017 •

edited

Loading