Implement HMA

[ajschult]

Is your feature request related to a problem? Please describe.

Harmonically Mapped Averaging (HMA) allows crystal properties (energy, pressure, heat capacity) to be computed more efficiently (more precision in the same time or the same precision in less time). Other benefits include reduced finite-size effects, reduced truncation effects, shorter decorrelation time and faster equilibration.

Describe the solution you'd like

We'll be adding code to compute HMA properties.

Describe alternatives you've considered

HMA properties could be computed by dumping forces and coordinates to files, but the output would be huge.

Further information, files or links (optional)

HMA paper: https://doi.org/10.1103/PhysRevE.92.043303 HMA in LAMMPS paper: https://doi.org/10.1063/1.5129942 HMA for VASP paper: https://doi.org/10.1016/j.cpc.2020.107554 LAMMPS pull request: https://github.com/lammps/lammps/pull/1503

[ajschult]

I've starting working on this. The first issue is that I'm not sure the best way to compute the Fi * dri sum that we need. We could

1. have the various routines in energy_routines compute the contribution to the force for that atom (rather than just using the potential derivative to compute the virial). This would require forces be stored for each atom (in the Atom_Class?), but would keep the details of HMA out of energy_routines. Other parts of Cassandra could use these forces at some point (force-bias MC, etc).
2. have the various routines in energy_routines compute the contribution to Fi*dri. This would use less memory, but would require energy_routines to know about HMA and would align better with how Cassandra already computes the virial.

[rsdefever]

Based upon our previous conversation I had been thinking we would go with (2). Is one or the other approach more favorable for properties that require the second derivative of the force? I would think (2) would become more appealing?

[ajschult]

Right, something like option (2) would be needed for the second derivative anyway. I'll pursue that approach for F*dr. As a bonus, as long as the energy routines know about HMA anyway, we can use a lattice-based truncation scheme that makes everything more precise and accurate.

[ajschult]

Another issue is how to read in the harmonic pressure. This needs to be computed outside Cassandra (with lattice dynamics) and provided somehow in the input file. It's the Delta P(hat) term that shows up in the pressure formula in Table 1 in 10.1103/PhysRevE.92.043303. Can I look for a second parameter on the Property_Info line?

# Property_Info
Energy_Total
Energy_HMA
Pressure
Pressure_HMA 9

and then the harmonic pressure would be taken as 9? It looks like it would be easy enough to implement in terms of the code, but it deviates from the other parameters (none of which take any extra argument).

[rsdefever]

Another option would be # HMA section in the input file. You can use a boolean to toggle the HMA on or off (this way you don't have to delete/add sections to the input file if you want to toggle on/off). If the # HMA section is not there then HMA is off. If # HMA is False, then any HMA related commands in the # Property_Info would be ignored.

# HMA
true
9

or

# HMA
true
harmonic_pressure 9

Thats my thoughts at the moment at least, but its also a bit clunky to have two different sections to edit. Let me chew on it a bit more.

[rsdefever]

Another issue is how to read in the harmonic pressure.

Is there any other information that HMA will need from the user?

[ajschult]

I don't think we need any other information for HMA. The HMA section would make more sense if there were more bits of information that we need. But it would certainly be more direct and align better with how other parts work. And it would work better if we ever come up with more information we want later.

[rsdefever]

@ajschult conclusion is to leave it in the # Property_Info section as you first suggested. As long as you don't think there are other bits of information we need, we're fine with that (small) deviation from the standard.

[ajschult]

I implemented the lattice cutoff for HMA (cutoff is based on the lattice sites rather than the current coordinates), which improves the precision further and allows people to use a "cut" potential rather than some flavor of shifting. But it's a bit uglier than I had hoped. You can see it here:

https://github.com/MaginnGroup/Cassandra/compare/master...etomica:HMA

Most of the energy_routines.f90 changes are to handle this except for the HMA_calc block at the bottom. Particularly sad is the need to set rabx, raby and rabz to atomic rxij, ryij and rzij.

Let me know your thoughts and/or if you have ideas of how to do this more cleanly.

[rsdefever]

@ajschult sorry for the slow reply.

I've taken a look at your branch. I don't think its too bad, but I agree it would be nice if it was a little cleaner. I'm not terribly bothered by setting rabx = rxij; Let's just add some more comments in there to warn some future programmer what is going on.

I'm not sure that this solution is better, but food for thought: Could we add a Check_AtomPair_LatticeCutoff subroutine:

IF (need_HMA) THEN
  CALL Check_AtomPair_Lattice_Cutoff(atom1, atom2, get_vdw, get_qq, this_box)
ELSE
  CALL Check_AtomPair_Cutoff(rijsq, get_vdw, get_qq, this_box)

Inside Check_AtomPair_Lattice_Cutoff you could compute rxijp, ryijp, rzijp from the lattice sites, call Minimum_Image_Separation, compute the rijsq between the lattice sites, and then call Check_AtomPair_Cutoff.

It's a little clunky in that the Check_AtomPair_Lattice_Cutoff and Check_AtomPair_Cutoff take different arguments, but it would minimize the extra code in energy_routines.f90 and let the rest of the machinery there basically stay the same. What are your thoughts?

Also, IIRC, we use a molecular virial, so you'll probably want to double check that you're getting what you expect there. I'm betting it works out correctly, but just a heads up.

[ajschult]

I've implemented what you described now. I think it's an improvement; it's less code and easier to follow. In theory, it's a bit slower because it has to compute the atom-atom separation even if it won't be used (because the lattice sites are outside the truncation), but it's already skipping the molecule-molecule distance computation, so the lattice cutoff approach is still just as fast as the standard approach.

[rsdefever]

Hey @ajschult, just wanted to touch base on this and see how it was coming. Is there anything else you need from us? No worries at all if things got busy :+1:.