GEMC simulation results for custom Helium potential

I have been in email correspondence with Ed and discussed this project with Eliseo and Ryan Mullen in person. At this point, we feel like we should use this forum to document the challenges that are impeding future progress.

Recap of what has been done to date:

1) Downloaded Cassandra V1.2
2) Replaced the LJ/Mie potential in the Compute_AtomPair_Energy subroutine with a call to "V" in a separate "potentials.f90" file.
3) Verified that the units were correct for the potentials.f90 code, i.e., this was originally written to work with bohr and hartree, but now it takes in angstrom and returns K.
4) Verified that the same configuration yields approximately the same energy using Towhee and Cassandra
5) Performed GEMC simulations with both the Towhee (tabulated) and Cassandra (non-tabulated) forms of the potential using the same system size and a 1.4 nm cutoff without tail corrections (since Towhee does not allow for tail corrections with tabulated potentials).
6) Having found discrepancies between the Towhee and Cassandra results, we performed simulations in both Cassandra and Towhee for the LJ potential. To demonstrate that Towhee's tabulated potential is not the problem, we performed this test using both the Towhee native LJ and Towhee tabulated with a 1.8 nm cutoff and no tail corrections. The results for Cassandra and Towhee LJ systems were indistinguishable and agree with the literature empirical correlation.

The attached figure helps demonstrate the issue we are facing.

Attached are the required source files.

All scripts and results can be found on our GitHub page, github.com/ramess101/Helium_ab_initio/Cassandra/ . This topic is discussed in Issue #1 as well.

Some details regarding Cassandra implementation:

We have generated the GEMC results on two different machines, but always with the gfortran compiler. To exclude tail corrections we simply set epsilon to 0 in the .mcf file (since the Mie potential is not actually used anywhere except in tail corrections, this is just a patch but should not cause any problems). We use a 1.5 rcutoff_low because we did not see any difference between this cutoff and the 0.5 rcutoff_low results. We have not tried to optimize our MC move probabilities or kappa for CBMC. These values were adopted from the tutorial files and seem to provide adequate insertion/deletion acceptance. We played around with rcut_cbmc (it was originally 6.5) but changing this to 14 did not seem to have any impact.

Do you have any further recommendations for how we can elucidate the problem?

Can you look over our implementation of V(rij) in energy_routines.f90?

Thank you

Eliseo,

Thank you for responding to me question. I need to apologize though, I thought that Ed had removed this question after I emailed on August 22nd saying that I found this energy units issue. Indeed, by converting to amu A^2/ps^2 we get the correct GEMC results. I hope that you did not spend too much time trying to debug our code.

Thanks again

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_ThreeBody_Energy within the move_* files (e.g., move_delete) 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