Bond valence constraint

We used bond valence model (BVM) as a constraint for CTR fitting. BVM is a robust constraint, which not only makes a faster convergence of differential evolution optimization, but also ensures a fitted model bearing a structure of physical sense. A structural model with extreme short bonds (thus physically impossible) will be essentially ruled out if the model is fitted under BVR constraint. The theory of BVR has been dedicated in detail in this paper. Briefly, a stable ionic structure always features all constituent atoms bearing bond valence saturation, i.e. the valence charge of an ion has to be compensated by its adjacent counter ions in the nearest polyhedra shell. While the valence charge of each ion is already known in a crystalline structure or can be solved under specific solution chemistry using thermodynamic speciation calculations, the bond valence of the covalent bonds formed between the targeted ion and its neighboring counter ions can be calculated using a simple empirical equation as

bond valence = exp((Ro-R)/B)

where R is the observed bond length, the B is an empirical factor (usually 0.37), and R0 is the fitted bond valence parameter representing the nominal length of a bond of unit valence, depends on the sizes of the atoms forming the bond. R0 values for different ions have been tabulated in this cif file. With this equation, we can easily check the bond valence states (undersaturated, oversaturated, or saturated) of all ions in the studied structure. It should be realized that a optimized structure model is likely a reasonable model if the model could pass the bond valence check, whereas a model failing the bond valence check has a high chance of being a wrong model. However, a model passing the bond valence check is not necessarily a correct model, especially if the model fit is poor with many misfits. In this situation, you should recheck the model you have constructed and look for any possible missing parts that could lead to the misfits, e.g. absence of interfacial water structure, assigning reactive sites incorrectly, presence of multiple reactive sites. Try a new model that considers one or a few such features. You may need to repeat this process several times until you find a model, which not only fits your data well but also passes the bond valence check.

In order to calculate the panelty factor accordidng to BVM, the program need to explore all atoms in one domain and figure out the geometrical relationship among all constituent atoms, i.e. the ids of all atoms that any targeted atom is associating with to form covalent bonds. It is only meaningful to count those nearest adjacent atoms, since the other atoms further away should make negligible contribution to the covalence although it will give rise to a small bond valence calculated using the bond valence equation. These small valence contributed by non-nearest neighboring atoms should be zeroed otherwise they will add up to make a big faked value that would mislead the fit engine for finding wrong model. To make thee BVM constraint more lenient, we only consider the case of oversaturation as a bad case (i.e. high panel factor), while the cases of undersaturated is considered as a good case like the saturated case. This is reasonable, since we could not include H atoms in the model (it is essentially invisible to X-ray) although they should exist in any aquatic fluid system playing a big role in terms of (de)protonation of undercoordinated oxgen groups as well as weaker hydrogen bonding, both of which could fill the valence gap of undersaturated atoms. Therefore, we only need to deal with oversaturated atoms for calculating the panelty factor. For the solution of bond valence panelty factor, We used to rely on the absolute difference between the calculated total bond valence and the associated ideal valence charge. In this scheme, we simply sum up the valence difference for all oversaturated atoms to come up with the panelty factor. Afterwards, we realized the bond valence sum calculated using the empirical equationi greatly diverge from the reasonable value when the observed R is way smaller than the associatedd R0 simply because of the exponential form in the equation. As a result, this extrem case (short R values) will introduce a too large weighting effect that usually end up with an unstable model. Afterwards, we take a differnt (simpler) strategy by uniformly appling a same arbitrary panelty factor (I used 10) for all oversaturated atoms. It seems much more robust to constrain the model fit that way. Thus, we stick to this strategy in the final bond valence script.

The class object dealing with bond valence constraint is defined as bond_valence_constraint in FitEnginePool module ([see here])(https://github.com/jackey-qiu/DaFy/blob/master/EnginePool/FitEnginePool.py). The underlying distance calculation subroutine has been optimized based on using scipy.spatial.distance.pdist. I will not detail this. Just check the code block if you want to know more. In the model script, you should switch the bond valence constraint off if you don't want to use this feature (for example if you only fit specular data), then set the attribute USE_BV in the model script to False like this

#/bv/begin#
USE_BV = True
#/bv/end#

End.