Update the unit test and documentation for two Yukawa model

sasmodels/models/TY_YukawaSq.py

@@ -1,20 +1,79 @@
 r"""
-Definition
-----------
 
-Calculates test2yv2.
+This model calculates the interparticle structure factor for a colloidal system with
+an interaction potential consisting of two Yukawa terms. If one Yukawa term 
+is chosen to be an attractive potential and another one to be a repulsive potential,
+the two Yukawa potential can simulate a potential with both a short-range
+attraction and long-range repulsion (SALR). [1,2]. When combined with an
+appropriate form factor $P(q)$, this allows for inclusion of the
+interparticle interference effects due to a relative complex potential.
+
+The interaction potential $V(r)$ is
+
+.. math::
+
+    \frac{V(r)}{k_BT} = \begin{cases}
+    \infty & r < 2R \\
+    K_1 \frac{e^{-Z_1(r-2R)}}{r} + K_2 \frac{e^{-Z_2(r-2R)}}{r} & r \geq 2R
+    \end{cases}
+
+And $R$ is the radius_effective.
+$K_1$ ( or $K_2$ ) is positive when there is a repulsion.
+$K_1$ ( or $K_2$ ) is negative when there is an attraction.
+
+.. note::
+
+   This routine only works for a potential with two Yukawa terms. 
+   If the amplitude, $K_1$ or $K_2$ of either Yukawa potential
+   becomes zero, the routine may self-destruct! But one can make either
+   amplitude very small if there is a need to make one amplitude to be zero.
+   Also, $Z_1$ and $Z_2$ should be different. If these two numbers are identical,
+   it essentially becomes a potential with one Yukawa term. 
+   But these two values can be very close to each other.
+   Therefore, the initial parameters of $Z_1$ and $Z_2$ should be different. 
+
+.. note::
 
+   Earlier versions of SasView did not incorporate the so-called
+   $\beta(q)$ ("beta") correction [3] for polydispersity and non-sphericity.
+   This is only available in SasView versions 5.0 and higher.
 
+This code calculates the structure factor, $S(q)$, of one component liquid systems interacting 
+with a two-term Yukawa potential with the mean spherical approximation (MSA). 
+The structure factor is generated by following Blum's method,[4] 
+in which the structure factor is solved using the Ornstein-Zernike equation by Baxter's Q-method with the MSA closure.
+
+The algorithm used in this model was originally proposed and developled by Liu et al. in 2005 and implemented using MatLab.[1]
+The algorithm was later reimplemented using C and Igor.
+SasView used the C codes developed by M. Hennig in 2010. 
+
+When the overall attraction is not very strong, this algorithm produces reasonably accurate results.[2]
+However, whent the net attraction is very strong. It has been shown that the fitting algorithm 
+tend to underestimate the attraction strength. 
+Accuracy of this algorithm was evaluated by Broccio et al..[2]
+
+When using the OZ euqation to obtain the structure factor, it assumes that a system is in a liquid state.
+However, when the attraction potential is too strong, a physical system may not be in a liquid system any more.
+But the OZ equation may still have a solution. But the result may not be physically meaningful.
+When this happends, the solution of the OZ equation may return unphysical results or fail to return any result, 
 
 References
 ----------
 
+#. Y. Liu, W. Chen, S-H Chen, *J. Chem. Phys.*, 122 (2005) 044507
+
+#. M. Broccio, D. Costa, Y. Liu, S-H Chen, *J. Chem. Phys.*, 124 (2006) 084501
+
+#. M. Kotlarchyk and S-H Chen, *J. Chem. Phys.*, 79 (1983) 2461-2469
+
+#. J. S. Hoye, L. Blum, J. Stat. Phys., 16(1977) 399-413
+
 Authorship and Verification
 ---------------------------
 
-* **Author:** --- **Date:** 2021YYY-05m-19d
-* **Last Modified by:** --- **Date:** 2021YYY-05m-19d
-* **Last Reviewed by:** --- **Date:** 2021YYY-05m-19d
+* **Author:** Yun Liu **Date:** January 22, 2024
+* **Last Modified by:** Yun Liu **Date:** January 22, 2024
+* **Last Reviewed by:** Yun Liu **Date:** January 22, 2024
 """
 
 from sasmodels.special import *
@@ -54,7 +113,7 @@
 
 #haveFq = True   # for beta calculation
 #single = False  # make sure that it has double digit precision
-#opencl = False  # related with parallel computing
+opencl = False  # related with parallel computing
 
 #Iq.vectorized = True
 
@@ -64,3 +123,13 @@
 #    return oriented_form(x, y, args)
 ## uncomment the following if Iqxy works for vector x, y
 #Iqxy.vectorized = True
+
+# The test results were generated with MatLab Code (TYSQ22) (Jan. 22, 2024)
+# Note that this test follows the definition of the original code : k1 < 0 means repulsion and k1 > 0 means attraction.
+# Paul K. and Yun discussed to switch sign of k1 (and k2) so that negative values mean attraction and positive ones mean repulsion.
+# Once the code is changed, the input values for k1 and k2 of the unit test need to change the signs.
+tests = [
+    [{'scale': 1.0, 'background' : 0.0, 'radius_effective' : 50.0,
+      'volfraction' : 0.2, 'k1' : 6, 'k2':-2.0, 'z1':10, 'z2':2.0,'radius_effective_pd' : 0},
+     [0.0009425, 0.029845], [0.126775, 0.631068]],
+]