Insights on force-matching potentials

Hello everyone,

I am currently using the force-matching (FM) method to derive classical force-fields. As me and the group I work with are new with it, I thought of sharing some thoughts on something I have been considering doing to get some opinion/criticism. (with roasting also totally allowed, if you feel like it)

So, basically I am trying to derive classical potentials for a coarse grained representation of a xtalline solid using FM. I am considering bonded and non-bonded potentials. Within FM, I am assuming in the case of the latter that the interactions in my system can be classically modelled by pairwise forces with the standard radial symmetry (which is not necessarily true actually). The mapping in hand has only one bead type (let’s call it 1). In my case, for bonded potentials, I have only one angle type and one bond type.

My reference force data comes from an atomistic simulation ran with an already-existing classical model reported in the literature. For what is worth it, I am not using the original first-ever published force-matching formalism by Ercolessi and Adams but rather one that was used by Voth et al (Effective force fields for condensed phase systems from ab initio molecular dynamics simulation: A new method for force-matching | The Journal of Chemical Physics | AIP Publishing) some time later, being, in essence, the difference simply that in the latter case an overdetermined system of equations is assembled to seek for optimal force-fields of a given form.

When working with force-matching, some people say that it is okay to consider the forces of a given contribution as being 0 in the non-sampled regions (for example, in the case of the 1-1 interaction, the non-sampled regions would be the values of distances 1-1 that are never explored, which we can find out by checking an RDF of a reference system). Others say that it is necessary to have sampling throughout the entire continuous range of distances one is trying to derive a force-field for. Particularly, I tend to believe in the latter case because I think that my chances of converging to potentials that are physically meaningful are larger (fundamentally, force-matching, in its mathematical formalism, is somewhat free to deliver phisically weird potentials).

My system (a crystalline solid) has actually some non-sampled regions: as it can be seen in figure (a), the peaks of the reference RDF for the mapping are not merged to one another. Some people that worked with FM for solids sampled microstates in the molten state and therefore had data to fit forces in regions of the rij distances that are not sampled in the xtalline state. I, however, cannot do such thing in the sense that, not only the material I am studying is not supposed to melt in “real life”, but also, even if I was to try to force it to happen computationally, my bonded force-fields in the AA model have harmonic-like shape, which I doubt would allow me to compromise the xtall structure. What I can do is to simulate it at a high temperature thermodynamic state and, as I have a larger degree of motion, the RDF peaks merge and therefore I would have data in the entire range of r 1-1 in which the non-bonded potential is supposed to act. This is also something that other people seem to have done.

By doing this, I successfully managed to get the peaks corresponding to non-bonded pairs in figure (a) to merge a little as they become wider. Yet, when I did the FM, I ended up getting potentials that seem a “little weird” to me and to other potentials I’ve seen reported in the literature (at least for the bonded contributions). I mean, I dont doubt that the FM protocol was carried out correctly in the sense that the potential for bonds + angles + pairwise 1-1 can in theory be a most optimal solution within FM mathematically speaking. But physically speaking, it is weird… for example, my bonded potentials are all repulsive. If this is the case, it should not be difficult to have bonded atoms straying too far away then, unless somehow the other contributions were able to balance it out. Indeed, as it turns out, there is no such balance and my simulation ends up collapsing.

It could be that I cannot derive potentials for my system using force-matching. It could also be (maybe more likely) that the assumption of modelling interactions with these “isotropic” pairwise potentials are valid. I am still to investigate other potential forms and symmetries. But another thing that I thought I could do is to do an extrapolation of the bonded potentials in the way shown in figure (b) for the potential for bonds. I mean, I am cheating in the sense that by extrapolating the potentials in this way I will be forcing the bonded contributions into respecting the xtal structure that I am expecting for my system. Not to mention that I will be introducing forces that will not really have come from force-matching - although some people sometimes use hybrid approaches, so it’s not like they are deriving “pure force matching” potentials is compulsory either. At the same time, as the extrapolated part of the potential fights to bring back the pairs to the sampled regions (part of the potential derived via FM), I will not be allowing my system to sample regions that should be unsampled in my extrapolation-methodology. In fact, I will be colaborating for a potential that is able to reproduce the structure also as a consequence of ther harmonic shape I am using for the extrapolation (see figure (c), which corresponds to reference and model structural data at 300K). What do you guys think?

Hi @ceciliaalvares,

I am only bumping in your post now and was curious if you made any progress on this. It is weird to me to turn to FM for crystal solids but why not? However with a regular crystal structure like that it seems unlikely that you will derive meaningful non-bonded interactions. What your force model is likely to look like is something looking like your rdf with discontinuous peaks and valleys having a weird physical interpretation for a force.

But I don’t know why you would not be able to derive a bond+anglular potential for your systems. Your coarse-grained beads are quite large: 0.6nm for nearest neighbor is large compared to atomistic resolution, so 1st neighbor bond interaction should be able to catch most of the long-range pair interaction between consecutive beads, right? If you switch off interaction at longer range, are they still purely repulsive? You might also need an angular force to maintain your proper structre though.

Hey @Germain

Well, technically force-matching was first introduced in the world in the context of modelling aluminium (including in its solid, xtalline state), [1] so I suppose it is not ultra absurd. Of course back then it was a work to derive classical potentials for the all-atom representation of the system from DFT simulations. If I remember correctly they used data from molten aluminium for force-matching together with the xtalline snapshots, so I guess they never really had “unsampled regions” anyways (a way to counter the problem). The FM force-field they derived was quite successful and also representative of the system as a whole since they used data coming from several different phases of aluminium.

When Voth et al brought this to the CG scale, the implementation they proposed changed a bit. Also, to counter non-sapled regions problem, Voth et al propose also this thing of considering the values of force to be 0 in non-sampled region. [2] They even say this can be “safely done”. I dont know if I agree… (respectifully), but it is true that the people in the paper are much more qualified than my humble opinion. I already think that force matching on its own is rather mathematical than physical when you think of the algorithm on its own (even if you have forces there). But maybe if they cosntrain the force to be 0 in the non-sampled regions this would allow you to converge to some more meaningful force-fields as solution for the algorithm. The softwares I used to do FM did not implement it like this, so I dont have any personal experience to report to you.

My progress on this was just a matter of doing a better, more reasonable extrapolation than the ugly one I showed in the figure (b) above. (In fact, now that I am looking at this again I dont know if this figure that I posted is correct… :thinking:)
This idea of deriving force-fields via FM at different temperatures is indeed not ideal since your potential may not at all be transferrable. But there are people who reported the temperature transferability in FM not to be so drastically bad. [3] In fact, they even propose a cool scheme to answer the question “what would be the rightful shape of my force-fields (within the FM algorithm) at temperature X given that in temperature Y, Z and W the shapes are like this?”. So following their reasoning might be a (third) way to solve the problem I discussed above. Unfortunately I do not have time to do this in a more proper way so what I did was to parametrize the force fields at a given (T,V) condition where T is much higher than 300K (specifically it is 1500K). The xtall structure was still maintained in the reference all-atom simulation (probably thanks to the bonded potentials…) and the peaks of the xtal structure broaden to the point where they merge, so that finally I have continuous sample for my non-bonded interactions (the non-bonded was excluded between 1-2 and 1-3 pairs, who interact via bonded potentials also derived via FM). To be honest the structure was not so bad. I can leave some pictures below for you to see the structure of MARTINI vs FM vs IBI for two different mappings of the system I am studying (ZIF-8) ^^
The figurs concern ambient conditions (T = 300K, V = Vexperimental at (300K, 1 bar)).
The potential form I will send to you via inbox because I feel more confortable (I mean, I didnt even publish these things and I am a bit sistematic sometimes), but they are not unsampled like the g(r) would prescribe it to be.
PS: figure (c) above is wrong because when I was doing my last correction later I realized I was running a force-field in kJ/(mol.nm) or something like that using real units in LAMMPS heheheheh :monkey:

Mapping A:

Mapping B:

Note also that these potentials of force-matching are not ready to go because in fact, the corresponding values of pressure force-matching potentials tend to give are super misguided (in my case it was not different) and you need to fix it first (in case you car eabout pressure and things that depend on pressure). [4,5,6] There are different ways to do this. I am still working on that so I dont have the “final” force-matching potential yet (this is why there is a ‘*’) in my figures.

[1] DOI 10.1209/0295-5075/26/8/005 (sorry, I had linked the wrong paper last night related to gold from the same authors and mixed information)
[2] DOI: 10.1063/1.1739396
[3] Systematic study of temperature and density variations in effective potentials for coarse-grained models of molecular liquids | The Journal of Chemical Physics | AIP Publishing
[4] The multiscale coarse-graining method. V. Isothermal-isobaric ensemble | The Journal of Chemical Physics | AIP Publishing
[5] The multiscale coarse-graining method: Assessing its accuracy and introducing density dependent coarse-grain potentials | The Journal of Chemical Physics | AIP Publishing
[6] Bottom-up coarse-grained models that accurately describe the structure, pressure, and compressibility of molecular liquids | The Journal of Chemical Physics | AIP Publishing

1 Like

By the way, just for curiosity: yesterday when I was writing the message above and searching for the official force-matching publication, I mixed it with another publication from the main author that has nothing to do with force-matching where they fit a glue potential to model gold. I edited out the wrong information from the above post but the publication of glue model for gold does serve to illustrate another way in which you could maybe implement force-matching and get rid of the sampling problem. It would require though the probably problematic task of knowing an analytical potential form that allows to properly reproduce, to one extent or another, the target forces experienced by the atoms.

Just a brief more detailed description of how the force-matching algorithm works to make my point:

When applied to coarse graining, you basically need to do an MD simulation in the all-atom representation of the system using force-fields that somewhat allow you to reproduce whatever it is you are interested in reproducing. So basically, once your system has equilibrated, you save a series of snapshots containing information on the position of each and all atom i in the system and also the force experienced by each of them at the given time instant. Note that you do not know, without doing some math, how each other atom j contributes to that value of force on atom i: you simply know the resulting force on atom i at that configuration.

Suppose now that I take each of these snapshots and coarsen the system according to the desired mapping. In principle, the center of the beads can be placed wherever you want (center of mass, center of geometry). You could also determine the “reference (“rightful”) forces” these beads should be experiencing by summing the forces of the atoms composing the given bead. This would look like the drawing I made below where I, II and III are the superatoms indexing. The blue arrows are the resulting force experienced by the beads in the given configuration.


Now note that you do not really know at all how superatoms II and III correspond to the total force experienced by superatom I for example. But whatever that force might be, if you decide that you can model my interactions in pairwise fashion, you can say for sure that the force experienced by atom I, II and III are given as shown below. Then, given that you have several snapshots (and a lot of superatoms), the idea is to ultimately find a combination of functions for the pairwise force-fields that most accurately reproduce the reference forces each superatom must be experiencing.

The mainstream idea of force-matching, be it the one by Ercolessi & Adams or Voth et al, is to model these pairwise force-fields F_I,II and etc with piecewise cubic polynomials. So you discretize your space and ultimately get something as output that would be a table potential with the values of F_I,II and etc that should occur at different grid points in order for the force at each superatom to live up to the reference ones in the scope of the snapshots considered to do the force-matching. But if you happen to know a due analytical form to model the interactions, you could try to implement FM in a way where the algorithm will know that the form of F_I,II and etc must be a given one and will then now optimize the corresponding parameters of the underlying analytical forms such that the force-fields will be the most rightful ones to match the forces. It is not what is usually done since you dont actually know what would be the “most suitable” function form to match the forces, but it is in principle quite doable. It would be another way in which you could counter the fact that you dont have sampling everywhere in the context of xtalline solids at reasonably low temperature. There is a paper of Voth et al where they do this for a Coulombic interaction and derive what would be the most rightful values of charges that allows for a minimal error in the force-matching algorithm (although not in the context of xtalline solids) [7]. If I remember correctly they find super small values of charge and conclude that they didnt even need to consider this in first place, but oh well…

[7] A multiscale coarse-graining method for biomolecular systems - PubMed
Random paper of Ercolessi et al where they educately suppose a given analytical form for the given potential modelling the interactions of gold (