LAMMPS implementation of angle_style cosine/periodic contains undocumented factor 2

[acroy]

The implementation of angle_style cosine/periodic in LAMMPS contains an undocumented factor of 2, which means that lammps_interface leads to angle bending constants which are too large for UFF and other force fields.

The problematic lines in the source code of LAMMPS are L140 and L298. For example, the latter is

return 2.0*k[type]*(1.0-b[type]*powsign(multiplicity[type])*c);

At least for UFF the factor of 2 is not wanted.

As an update of the LAMMPS documentation is more likely than a change of the implementation, I would suggest to account for the factor 2 in the respective force field methods.

[ltalirz]

thanks for the report @acroy ! happy to accept a PR to fix this.

can someone else confirm this @peteboyd @SeyedMohamadMoosavi @danieleongari

[peteboyd]

I remember getting the force field terms correct for UFF required a lot of testing. One can see the angle_term function for UFF in ForceFields.py starting on line 2158. The actual definition of the 'C' term, defined as kappa in the code, is on line 2261, but is dependent on the type of bond and is effectively Ka/n**2, where n is the periodicity of the angle term. Line 2261:

data['potential'].C = kappa*(c1**2)

where c1 is considered the multiplicity[type], or 'n', term in the cosine/periodic equation and could take values of 1, 2, 3, or 4 depending on the desired bond angle geometry. If I recall correctly (from 4 years ago) tests where we set data['potential'].C = kappa resulted in overly floppy and soft materials, so perhaps accounting for this undocumented factor of 2 would yield very loose materials, but this would need to be tested before we push the change. I currently have very little time to do this and I would welcome contributions from @acroy, @SeyedMohamadMoosavi, and of course, the infallible @danieleongari.

[acroy]

Testing is of course always good. But note that this issue has nothing to do with UFF, it is a "problem" with LAMMPS. I have opened an issue and Axel Kohlmeyer confirmed the discrepancy between implementation and docs. Basically, LAMMPS calculates E = 2C/n^2 [1 - B (-1)^n cos(n Θ)] which gives E ~ C (Θ - Θ_0)^2 instead of E ~ C/2 (Θ - Θ_0)^2. Neither the factor 2 nor the normalisation 1/n^2 was documented - but somehow you captured the normalisation :-)

Regarding a PR: I was wondering if this should be fixed in ForceFields.py or rather in lammps_potentials.py? Personally, I would put it in the latter and keep ForceFields.py as independent as possible from the implementation used in the MD code.

[peteboyd]

Ok what the hell why not. I pushed the change, but I won't close the issue just yet.

[ltalirz]

Let me know once you're convinced of this change and I'll drop a new release

[ltalirz]

Hi @peteboyd - I was just contacted by @RaffaelaCabriolu regarding this since both #47 and #44 have quite dramatic effects e.g. on thermal conductivity (as one might expect).

Just checking

• Are you confident - within reason - that the current implementation is correct?
• Does this, conversely, mean that using lammps-interface versions prior to v0.2.1 likely yielded bogus results for dynamics calculations with UFF/UFF4MOF?

E.g. when you were checking for "floppiness" which were the checks that you performed? If you let me know the input & expected output, I can easily add some tests (like the one mentioned in #44) to the continuous integration tests in order to prevent issues from reoccurring.

[RaffaelaCabriolu]

In particular, the differences I can see in the force field parameters are both, in the dihedral and angle coefficients for UFF, I can also document certain differences for UFF4MOF. The versions prior v0.2.1 are giving a total different value of the thermal conductivity from UFF and UFF4MOF, furthermore, there is a peculiar behavior that could be physical or due to the wrong estimation of certain coefficients. I could calculate the phonons if that is useful. It would be helpful if the actual versions are reproducing the experimental structures better than the previous ones.

[peteboyd]

Hi @peteboyd - I was just contacted by @RaffaelaCabriolu regarding this since both #47 and #44 have quite dramatic effects e.g. on thermal conductivity (as one might expect).

Just checking

• Are you confident - within reason - that the current implementation is correct?
• Does this, conversely, mean that using lammps-interface versions prior to v0.2.1 likely yielded bogus results for dynamics calculations with UFF/UFF4MOF?

E.g. when you were checking for "floppiness" which were the checks that you performed? If you let me know the input & expected output, I can easily add some tests (like the one mentioned in #44) to the continuous integration tests in order to prevent issues from reoccurring.

Regarding the faithful implementation of force fields like UFF, I'm certain the contributors have found issues and they were corrected. I do not remember how we tested floppiness, aside from visual inspection and subsequent comparing with what we know about AIMD trajectories, for example.

Are the previous results bogus? That's difficult to ascertain. John von Neumann said "With 4 parameters I can fit an elephant, and with 5 I can make him wiggle his trunk". If we are obtaining a better reproduction of experimental results with the previous implementation then either there are other incorrect parameters we haven't discovered, or, there is something fundamentally lacking with the UFF and UFF4MOF implementations for MOFs. I think moving forward, we have to consider some other literature reports using UFF4MOF and UFF in MD simulations that haven't used this code to parameterize the system, and see if we can reproduce the same numbers.

My main concern about automated testing of FF implementations is that one would have to compare results from a LAMMPS run, which would make the testing process cumbersome and long. Perhaps necessary though.

[ltalirz]

Hi Pete, thanks for the quick reply!

My main concern about automated testing of FF implementations is that one would have to compare results from a LAMMPS run, which would make the testing process cumbersome and long. Perhaps necessary though.

If the required lammps calculation itself can be made quick (and I think it often can; e.g. an energy minimization), it's actually quite straightforward to set up on CI (e.g. via conda). Happy to do this if you let me know what exactly to check.

[peteboyd]

Thanks for your continued stewardship of this code @ltalirz, you have been invaluable! I've re-opened the issue and I'll think about this in the coming weeks. @RaffaelaCabriolu can you go into more details about the differences in thermal conductivity?

[RaffaelaCabriolu]

This is a graph reporting the differences in the thermal conductivities using Non equilibrium methods.

[RaffaelaCabriolu]

The graph above presents the thermal conductivities as a function of the dimension of the sample for MOF-5, using Approach to equilibrium Molecular Dynamics (AEMD) that is a kind of Non Equilibrium Method (NEMD). With "new" in the legend, I refer to results coming from UFF and UFF4MOF out of the last version of lammps-interface. The version previous to the May 8 gives instead the continuous lines. I am calculating the phonons also and I will post them once they are ready. Interestingly enough, the old version of the UFF reports the anomaly of a Thermal conducitivity that grows linearly with the sample size and never converges. In this kind of method, (AEMD, approach to equilibrium molecular dynamics), the thermal conducitivity needs to be "scaled" with the heat capacity. In those graphs, for simplicity, I have used the Doulong Petit value for all the sample (i.e. 3*N with N the number of atoms). With this, I mean that you should not pay attention to the absolute value of the thermal conductivitiy that is not yet the definitive one, but consider the convergence to a plateau, that is the discriminant to define the potetential good enough to predict a finite thermal conductivity.

[peteboyd]

Thanks @RaffaelaCabriolu. What material is this?

[RaffaelaCabriolu]

MOF-5. Sorry, I edited the previous message for more clarity.

[peteboyd]

The graph above presents the thermal conductivities as a function of the dimension of the sample for MOF-5, using Approach to equilibrium Molecular Dynamics (AEMD).

Hi @RaffaelaCabriolu, sorry but I've never run a thermal conductivity calculation before, can you confirm the meaning of the x-axis for me? This should be a function of the temperature correct?

[RaffaelaCabriolu]

Yes, the Lx is the length of the sample in the direction parallel to the applied Temperature Gradient. I.e.:

[peteboyd]

Sorry this it turning into a bit of a group meeting, but I just saw this article from 2007. The x-axis is different, but Is this the type of behaviour we would expect? There is a bit of drift (and the MD calculations are bang-on!)

[RaffaelaCabriolu]

Sorry, I have not seen your post before. That graph you have posted refers to another kind of calculations, i.e. Green Kubo (GK) calculations, where you can retrieve the Thermal conductivity in all directions to which you apply the gradient. Since MOF-5 is isotropic, you can calculate only in direction. NEMD calculations, the one I have shown in my above graph, allows to represent the typical experiment (a lot of assumptions in there...) where you can apply your heat current in only 1 direction and measure the thermal conductivity in that direction. Ideally,with the NEMD calculations you should get the same results as GK and esperiments, but in practice there are some non-linear effect included in the NEMD. The graph above is showing GK and experiment for different Temperatures, myne is showing how the TC changes with the dimension of the sample. Unfortunately, NEMD and GK shows some size effect and my graph was done to understand what is the size at which the TC does not change anymore. We can have more discussion in a separate channel, if you want.

[RaffaelaCabriolu]

To add up, what we are expecting is a plateau in the results from the the AEMD approach. When the plateau is reached, the TC bulk size has been reached. From my graph, you see that the UFF4MOF brown continuous line does not grow anymore for sample with dimension Lx> 3000 Angstrom. For UFF red continuous line TC is always growing. In GK the size dependence is not so strong and it has been reported (https://www.sciencedirect.com/science/article/pii/S0017931006005473?via%3Dihub) by the same authors that have written the paper you have mentioned.