Three-body interactions

[ramess101]

The three-body module (he3body.f90) is now working properly. Here is a validation run:

[ramess101]

One concern that I am not sure we will be able to resolve is how to include the three-body interactions into the virial calculation. I am not even sure how you would do this mathematically, let alone how to implement it in code. We will probably just have to neglect this altogether. Although this would be one advantage of using GCMC where the pressures are not required.

[ramess101]

I posted this message on the Cassandra page:

Eliseo and/or Ed,

Now that we have the two-body Helium potential working, our next (and final) task is to implement a three-body Helium potential.

We already have the Fortran 90 code written that takes as input three interatomic distances (rij, rik, and rjk) and outputs the three-body energy. It was readily obvious for the two-body Helium potential that we just needed to call our two-body function in the Compute_AtomPair_Energy subroutine. However, since Cassandra does not have a "Compute_ThreeBody_Energy" subroutine, it is not as clear to me what the best approach is moving forward. Since you are much more familiar with the Cassandra code and have likely given some thought to the implementation of three-body potentials, I was hoping you could provide some insight.

Here is my (naive) idea so far:

1) After calling Compute_MolecularPair_Energy within Compute_Molecule_Nonbond_Inter_Energy (around line 884 in energy_routines), include an additional loop over all the molecules/atoms that are not "is" or "ispecies" (call this "ithird") 2) Call Check_MolecularPair_Cutoff for the distances between ithird and is/ispecies. 3) If these additional two distances are also within the cut-off, call Compute_ThreeBody_Energy and add this to the Eij_inter_vdw term (we could try separating this into its own variable if we want to)

Note that I'm not sure right now how we would avoid double counting the interactions with this approach.

One benefit of this approach is that we only have to loop through the third molecule/atom and we only do this after we have determined that the first two are within the cut-off. Another (more intuitive) option would be to call Compute_ThreeBodyEnergy within the move files (e.g., movedelete) after calling the other Compute subroutines (e.g., Compute_Molecular_Bond_Energy). From what I can tell, the downside to this approach is that we would be repeating the time consuming calculation of looping through N molecules, which is already being performed within Compute_Molecule_Nonbond_Inter_Energy.

In brief, do you think this approach would work? If you have concerns, what would be the simplest, fastest, or best way to implement a three-body potential in Cassandra?

Thank you