Problem with reproducing the validation results

[ruahdroid]

Dear developer,

I have patched the three-body potential into Lammps with the "Stable release 29 August 2024", and everything is working smoothly.

I wasn't able to reproduce the data presented in Section 4.1 of the manual. To align the example with your description, I modified the input file "ab.lj3b" in the lj3b-example by increasing n from 50 to 100.

According to "ref. 1" in the manual, "two particles are considered bonded if their relative distance is less than r_cutoff." I assume that the "r_cutoff" is 1.3 (referred to as R_cut_1 here) based on the values in "ab.lj3b." I also experimented with using 2**(1/n) (referred to as R_cut_2) as the criterion for determining bonding.

The issue is that the average number of bonds is 508 for R_cut_1 and 178 for R_cut_2. The value of 508 for R_cut_1 exceeds the maximum allowed number of bonds (400) for the given system setup.

The results are included in the attachment and resemble Figure 1 in the manual. Could you please assist me with reproducing the validation results? Thank you!

validation.pdf

[RodRVS]

Hi, thanks for reaching out. I am glad to know the potential can be added to LAMMPS latest stable release.

Regarding the input file "ab.lj3b", you are correct, the "n" parameter should be n=100. n=100 is the value used in the reference article and also the one mentioned in the manual pdf.

About the number of bonds being 508 for a bond definition of R_cut_1=1.3. I checked my data, funny enough I also get an average number of 508 "AB bonds" in my simulation (508.049744). However, please remember that the potential is meant to be a bond swapping mechanism, hence it is normal (in fact desired) that you will get more than the "allowed" number of bonds for this is when the potential kicks in and the swapping takes place. For this validation case in particular, we have extra 200 B particles unpaired that do the swapping. The important thing here though, is that if you monitor the contribution to the energy of the system due to the pairing of AB particles, you do not get more than the allowed for the case of 400 pairs (say -400\eps if you used the reference article parameters). This value should be readily available from the potential energy column from your log.lammps file.

I hope this answers your question, please let me know if it is not the case. best regards, Rodrigo

[ruahdroid]

Dear Rodrigo,

Thank you so much for the reply!

Indeed, I was also assuming that there are intermediate status for the bond pairing, so I no longer have any questions on this point.

However, how can I reproduce Figure 1 from the manual? In other words, could you provide the detailed procedure you used for the analysis? As you can see in the 'validation.pdf,' my results differ from those in Figure 1. My procedure involves analyzing the pair distance array, filtering out any pairs within the cutoff radius R_cut_1 = 1.3, and using all marked pairs as n(b, 0). I apply the same approach for all subsequent frames, obtaining n(n, t), and then calculate the autocorrelation function based on the fraction of remaining pairs in n(n, t) compared to n(b, 0). However, both the initial and final sections of my figure differ from what is shown in Figure 1. When plotting the figure, I also used multiple time origins.

Moreover, I think the equation 4 in the manual is not fully correct, it should be: The term 'q' should be in the exponent, not in the denominator. Additionally, the constant 'C' is missing from the 'fs3b' potential. Could you update the relevant section of the manual accordingly?

If possible, could you also add more author details to the README.md? I'd like to make sure I can acknowledge you correctly when the paper is ready. :) Or you can email me (ruahdroid@gmail.com) your information if you want to keep anonymous.

Best regards, Ruah