Add computation for truncated potentials with tail corrections about lumol HOT 14 CLOSED

I believe this would be better as two additional functions for PairPotential trait:

trait PairPotential {
   // Default return values is equivalent to today's code
   fn tail_energy(&self, cutoff: f64, density: f64) -> f64 {0.0}
   fn tail_virial(&self, cutoff: f64, density: f64) -> Matrix3 {Matrix3::zero()}
}

These function would be called from System::{energy,virial} with the actual density of the system, and the half of the lower cell parameter as default cutoff. (Is this correct for non-cubic boxes?)

Then, the CutoffComputation can implement these by overwriting the cutoff value.

impl<T: PairPotential> PairPotential for CutoffComputation<T> {
   fn tail_energy(&self, cutoff: f64, density: f64) -> f64 {
        assert!(cutoff >= self.cutoff, "Cutoff is too big in CutoffComputation");
        self.potential.tail_energy(self.cutoff, density)
   }

   fn tail_virial(&self, cutoff: f64, density: f64) -> Matrix3 {
        assert!(cutoff >= self.cutoff, "Cutoff is too big in CutoffComputation");
        self.potential.tail_virial(self.cutoff, density)
   }
}

from lumol.

That's also possible. Note that potentials that are computed via CutoffComputation must not be shifted when cutoff corrections are applied.

from lumol.

Forgot about this one. Maybe there should be another Computation? Like having CutoffComputation which only uses a cutoff for MC (but is not shifted), and ShiftedComputation which uses a cutoff and shift the potential for MD.

from lumol.

I'd like to tackle this. The current implementations do not distinguish between inter- and intramolecular potentials. Tail-corrections and shifts are usually only applied to intermolecular potentials.

Also, there are two correction terms: one for the systems' energy, one for the pressure (not for the forces).

I am not sure where to add this into the code. This is an MC-only issue. A naive implementation looks like this:
This does not work:

// We only need to recompute the tail-correction when the density or the cut-off changes, e.g. when
// i.e. the box size or shape changes or a molecule is added or deleted from the cell.
// 'self' is an EnergyEvaluator that owns a ref to system
let mut energy = 0.0;                                                   
  for i in 0..self.system.size() {                                                                                                                  
    for j in (i+1)..self.system.size() {                                
      let r = self.system.nearest_image(i, j).norm();        
      // pseudo code   
      // 'system' has to have access to rc and density
      // The following does not work. Tail() is a per particle contribution
      energy += self.pair(r, i, j) + self.system.tail(i, j);
    }                                                                   
  }                                                                       
return energy 

This would be the cleanest way I can think of but with a lot of overhead since for every atom pair there is one value for the tail correction.

A fasterThe correct way would be to move the tail computation outside the loop. We only need the mole fractions of the molecules (if mixture) to compute the corrections but then the systems' energy computation is different from uncorrected computations.

let mut energy = 0.0;                                                   
  for i in 0..self.system.size() {                                                                                                                  
    for j in (i+1)..self.system.size() {                                
      let r = self.system.nearest_image(i, j).norm();        
      energy += self.pair(r, i, j); // no changes here                                  
    }                                                                   
  } 
// pseudo: match for computation? if TailCorrectionComputation: apply 'tail()'                                                                      
return energy + self.system.tail() // density and rc and the mole fractions are known.

from lumol.

This would be the cleanest way I can think of but with a lot of overhead since for every atom pair there is one value for the tail correction.

This is what I had in my head with the tail_energy and tail_virial functions here.

Moving the computation out of the loop can be nice too, but may be harder to get right, as simple scheme like using mole fractions for computation will break for complex systems: how do you apply them for interfaces for example?

I am not sure where to add this into the code. This is an MC-only issue.

For now, the system have no idea whether it is propagated trough Molecular Dynamics or Monte-Carlo, so the code would need to run every time. I am not too worried about the overhead here, as this is only an O(1) cost that is added to the normal computation of pairs energy.

My idea to get away with the problem "should we compute tail correction" would be to add a specific CutoffComputation computation, or a TailCorrectedLJ potential, which would be the only one implementing tail corrections.

The current implementations do not distinguish between inter- and intramolecular potentials. Tail-corrections and shifts are usually only applied to intermolecular potentials.

It does! This is the reason to live for PairRestriction, which is currently used in the Interactions. I think they should live inside the PairPotential instead, and this would be a nice use case for rust-lang/rfcs#1546. It can still be done "by hand" while waiting, I may do that tonight!

from lumol.

Moving the computation out of the loop can be nice too, but may be harder to get right, as simple scheme like using mole fractions for computation will break for complex systems: how do you apply them for interfaces for example?

Corrections only should be applied to bulk systems since the correction assumes pair correlation functions of unity outside the cut off. For interface simulations I'd use a truncated potential.

My idea to get away with the problem "should we compute tail correction" would be to add a specific CutoffComputation computation, or a TailCorrectedLJ potential, which would be the only one implementing tail corrections.

Ok, I will have a look at this.

from lumol.

Corrections only should be applied to bulk systems since the correction assumes pair correlation functions of unity outside the cut off. For interface simulations I'd use a truncated potential.

Ok, so the issue is more about how do you define that you need tail corrections or not.

from lumol.

Another simpler solution just popped into my mind: why not unify the pair interactions? Having somehting like this:

pub struct PairInteraction<T: PairPotential> {
    // Use a generic here to remove the function dispatch cost
    potential: T,
    restriction: PairRestriction,
    cutoff: f64,
    shift: Option<f64>,
}

impl<T: PairPotential> PairInteraction<T> {
    fn new(potential: T, cutoff: f64) -> PairInteraction<T> {
        PairInteraction {
            potential: potential,
            restriction: PairRestriction::None,
            cutoff: cutoff,
            shift: None,
        }
    }

    fn shifted(potential: T, cutoff: f64) -> PairInteraction<T> {
        let shift = potential.energy(cutoff);
        PairInteraction {
            potential: potential,
            restriction: PairRestriction::None,
            cutoff: cutoff,
            shift: Some(shift),
        }
    }
}

// Only manipulate PairInteraction as Box<PairPotential> everywhere else in the code, to have a nicer interface.
impl PairPotential for PairInteraction {
   ...
}

Allow to have all the information together: the cutoff to use, whether to potential should be shifted or not, and which restriction to use. This would replace the CutoffComputation, and remove the need for ShiftedComputation.

Then, adding new methods for tail correction to PairPotential should be very easy. What do you think?

from lumol.

I like it. We just have to agree on what we actually want to have: truncation, truncation + shift (which shift?), truncation + tail correction? Should adding a tail correction be an Option<bool>? And where to put the pressure correction.

I usually work with:

• truncated (no shift, no correction, for example when simulating interfaces)
• truncated + tail corrected (bulk phases, solids)
• truncated + shifted (without corrections, to compare MD <-> MC data)

from lumol.

Yes, we need another bool value to know whether we need tail correction or not. An Option<bool> is a bit redundant here, a bool is enough.

which shift?

Shifting so that the potential is zero at the cutoff is the only shift kind I know. Are there any other?

from lumol.

Ups, you're right bool is enough.

Shifting so that the potential is zero at the cutoff is the only shift kind I know. Are there any other?

You can shift the potential to zero and add another contribution so that the forces are smooth and continuous. I never used that though. Maybe we should go with the basics first.

from lumol.

I am working on this.

It looks like you need tail corrections for shifted energy too: see this PDF, page 6; but as we are using a shift such as V(rc) = 0, the expression should be the same as the non-shifted case.

from lumol.

The lecture notes of Scott Shell are golden.

I am working on this.

Awesome. I forgot about the correction terms for the shifted cases. I have never used them to be honest. But to compare MC and MD simulations it is valuable to have them.

from lumol.

Leaving this open for tracking the implementation of this feature in the input files.

from lumol.