Merge pull request #12 from stevenvdb/damped_dispersion_power

Dispersion with Tang-Toennies damping for any power

yaff/pes/ext.pyx

@@ -1303,18 +1303,18 @@ cdef class PairPotDampDisp(PairPot):
 
         **Energy:**
 
-        .. math:: E_\text{DAMPDISP} = \sum_{i=1}^{N} \sum_{j=i+1}^{N} s_{ij} C_{6,ij} f_\text{damp,6}(d_{ij}) d_{ij}^{-6}
+        .. math:: E_\text{DAMPDISP} = \sum_{i=1}^{N} \sum_{j=i+1}^{N} s_{ij} C_{n,ij} f_\text{damp,n}(d_{ij}) d_{ij}^{-n}
 
         where the damping factor :math:`f_\text{damp}(d_{ij})` is optional. When used
         it has the Tang-Toennies form:
 
         .. math:: f_\text{damp,n}(d_{ij}) = 1 - \exp(-B_{ij}r)\sum_{k=0}^n\frac{(B_{ij}r)^k}{k!}
 
-        The pair parameters :math:`C_{6,ij}` and :math:`B_{ij}` are derived from atomic
+        The pair parameters :math:`C_{n,ij}` and :math:`B_{ij}` are derived from atomic
         parameters using mixing rules, unless they are provided explicitly for a given
         pair of atom types. These are the mixing rules:
 
-        .. math:: C_{6,ij} = \frac{2 C_{6,i} C_{6,j}}{\left(\frac{V_j}{V_i}\right)^2 C_{6,i} + \left(\frac{V_i}{V_j}\right)^2 C_{6,j}}
+        .. math:: C_{n,ij} = \frac{2 C_{n,i} C_{n,j}}{\left(\frac{V_j}{V_i}\right)^2 C_{n,i} + \left(\frac{V_i}{V_j}\right)^2 C_{n,j}}
 
         .. math:: B_{ij} = \frac{B_i+B_j}{2}
 
@@ -1325,8 +1325,8 @@ cdef class PairPotDampDisp(PairPot):
             indexes for the atom types that start counting from zero. shape =
             (natom,).
 
-        c6_cross
-            The :math:`C_{6,ij}` cross parameters.
+        cn_cross
+            The :math:`C_{n,ij}` cross parameters.
 
         b_cross
             The :math:`B_{ij}` cross parameters. When zero, the damping factor
@@ -1341,71 +1341,79 @@ cdef class PairPotDampDisp(PairPot):
             The truncation scheme, an instance of a subclass of ``Truncation``.
             When not given, no truncation is applied
 
-        c6s
-            The diagonal :math:`C_{6,i}` parameters
+        cns
+            The diagonal :math:`C_{n,i}` parameters
 
         bs
             The diagonal :math:`B_{i}` parameters
 
         vols
             The atomic volumes, :math:`V_i`.
 
+        power
+            The power to which :math:`1/d_{ij}` is raised.
+            Default value is 6, which corresponds to the first term in a
+            perturbation expansion of the dispersion energy.
+
         The three last optional arguments are used to determine pair parameters
         from the mixing rules. These mixing rules are only applied of the
         corresponding cross parameters are initially set to zero in the arrays
-        ``c6_corss`` and ``b_cross``.
+        ``cn_cross`` and ``b_cross``.
     '''
     cdef long _c_nffatype
-    cdef np.ndarray _c_c6_cross
+    cdef long _c_power
+    cdef np.ndarray _c_cn_cross
     cdef np.ndarray _c_b_cross
     name = 'dampdisp'
 
     def __cinit__(self, np.ndarray[long, ndim=1] ffatype_ids not None,
-                  np.ndarray[double, ndim=2] c6_cross not None,
+                  np.ndarray[double, ndim=2] cn_cross not None,
                   np.ndarray[double, ndim=2] b_cross not None,
                   double rcut, Truncation tr=None,
-                  np.ndarray[double, ndim=1] c6s=None,
+                  np.ndarray[double, ndim=1] cns=None,
                   np.ndarray[double, ndim=1] bs=None,
-                  np.ndarray[double, ndim=1] vols=None):
+                  np.ndarray[double, ndim=1] vols=None,
+                  long power=6):
         assert ffatype_ids.flags['C_CONTIGUOUS']
-        assert c6_cross.flags['C_CONTIGUOUS']
+        assert cn_cross.flags['C_CONTIGUOUS']
         assert b_cross.flags['C_CONTIGUOUS']
-        nffatype = c6_cross.shape[0]
+        nffatype = cn_cross.shape[0]
         assert ffatype_ids.min() >= 0
         assert ffatype_ids.max() < nffatype
-        assert c6_cross.shape[1] == nffatype
+        assert cn_cross.shape[1] == nffatype
         assert b_cross.shape[0] == nffatype
         assert b_cross.shape[1] == nffatype
-        if c6s is not None or vols is not None:
-            assert c6s is not None
+        if cns is not None or vols is not None:
+            assert cns is not None
             assert vols is not None
-            assert c6s.flags['C_CONTIGUOUS']
+            assert cns.flags['C_CONTIGUOUS']
             assert vols.flags['C_CONTIGUOUS']
-            assert c6s.shape[0] == nffatype
+            assert cns.shape[0] == nffatype
             assert bs.shape[0] == nffatype
-            self._init_c6_cross(nffatype, c6_cross, c6s, vols)
+            self._init_cn_cross(nffatype, cn_cross, cns, vols)
         if bs is not None:
             assert bs.flags['C_CONTIGUOUS']
             self._init_b_cross(nffatype, b_cross, bs)
         pair_pot.pair_pot_set_rcut(self._c_pair_pot, rcut)
         self.set_truncation(tr)
         pair_pot.pair_data_dampdisp_init(
-            self._c_pair_pot, nffatype, <long*> ffatype_ids.data,
-            <double*> c6_cross.data, <double*> b_cross.data,
+            self._c_pair_pot, nffatype, power, <long*> ffatype_ids.data,
+            <double*> cn_cross.data, <double*> b_cross.data,
         )
         if not pair_pot.pair_pot_ready(self._c_pair_pot):
             raise MemoryError()
         self._c_nffatype = nffatype
-        self._c_c6_cross = c6_cross
+        self._c_power = power
+        self._c_cn_cross = cn_cross
         self._c_b_cross = b_cross
 
-    def _init_c6_cross(self, nffatype, c6_cross, c6s, vols):
+    def _init_cn_cross(self, nffatype, cn_cross, cns, vols):
         for i0 in xrange(nffatype):
             for i1 in xrange(i0+1):
-                if c6_cross[i0, i1] == 0.0 and vols[i0] != 0.0 and vols[i1] != 0.0:
+                if cn_cross[i0, i1] == 0.0 and vols[i0] != 0.0 and vols[i1] != 0.0:
                     ratio = vols[i0]/vols[i1]
-                    c6_cross[i0, i1] = 2.0*c6s[i0]*c6s[i1]/(c6s[i0]/ratio+c6s[i1]*ratio)
-                    c6_cross[i1, i0] = c6_cross[i0, i1]
+                    cn_cross[i0, i1] = 2.0*cns[i0]*cns[i1]/(cns[i0]/ratio+cns[i1]*ratio)
+                    cn_cross[i1, i0] = cn_cross[i0, i1]
 
     def _init_b_cross(self, nffatype, b_cross, bs):
         for i0 in xrange(nffatype):
@@ -1418,17 +1426,17 @@ cdef class PairPotDampDisp(PairPot):
         '''Print suitable initialization info on screen.'''
         if log.do_high:
             log.hline()
-            log('ffatype_id0 ffatype_id1         C6          B')
+            log('ffatype_id0 ffatype_id1         cn[au]        B')
             log.hline()
             for i0 in xrange(self._c_nffatype):
                 for i1 in xrange(i0+1):
-                    log('%11i %11i %s %s' % (i0, i1, log.c6(self._c_c6_cross[i0,i1]), log.invlength(self._c_b_cross[i0,i1])))
+                    log('%11i %11i %10.5f %s' % (i0, i1, self._c_cn_cross[i0,i1], log.invlength(self._c_b_cross[i0,i1])))
 
-    def _get_c6_cross(self):
-        '''The C6 cross parameters'''
-        return self._c_c6_cross.view()
+    def _get_cn_cross(self):
+        '''The cn cross parameters'''
+        return self._c_cn_cross.view()
 
-    c6_cross = property(_get_c6_cross)
+    cn_cross = property(_get_cn_cross)
 
     def _get_b_cross(self):
         '''The damping cross parameters'''

yaff/pes/pair_pot.c

@@ -400,15 +400,16 @@ double pair_fn_ljcross(void *pair_data, long center_index, long other_index, dou
 
 
 
-void pair_data_dampdisp_init(pair_pot_type *pair_pot, long nffatype, long* ffatype_ids, double *c6_cross, double *b_cross) {
+void pair_data_dampdisp_init(pair_pot_type *pair_pot, long nffatype, long power, long* ffatype_ids, double *cn_cross, double *b_cross) {
   pair_data_dampdisp_type *pair_data;
   pair_data = malloc(sizeof(pair_data_dampdisp_type));
   (*pair_pot).pair_data = pair_data;
   if (pair_data != NULL) {
     (*pair_pot).pair_fn = pair_fn_dampdisp;
     (*pair_data).nffatype = nffatype;
+    (*pair_data).power = power;
     (*pair_data).ffatype_ids = ffatype_ids;
-    (*pair_data).c6_cross = c6_cross;
+    (*pair_data).cn_cross = cn_cross;
     (*pair_data).b_cross = b_cross;
   }
 }
@@ -434,33 +435,34 @@ double tang_toennies(double x, int order, double *g){
 }
 
 double pair_fn_dampdisp(void *pair_data, long center_index, long other_index, double d, double *delta, double *g, double *g_cart) {
-  long i;
-  double b, disp, damp, c6;
+  long i,j,power;
+  double b, disp, damp, cn;
   // Load parameters from data structure and mix
   pair_data_dampdisp_type *pd;
   pd = (pair_data_dampdisp_type*)pair_data;
   i = (*pd).ffatype_ids[center_index]*(*pd).nffatype + (*pd).ffatype_ids[other_index];
-  c6 = (*pd).c6_cross[i];
-  if (c6==0.0) return 0.0;
+  power = (*pd).power;
+  cn = (*pd).cn_cross[i];
+  if (cn==0.0) return 0.0;
   b = (*pd).b_cross[i];
   if (b==0.0) {
     // without damping
-    disp = d*d;
-    disp *= disp*disp;
-    disp = -c6/disp;
+    disp = 1.0;
+    for (j=0;j<power;j++) { disp *= d; }
+    disp = -cn/disp;
     if (g != NULL) {
-      *g = -6.0*disp/(d*d);
+      *g = -power*disp/(d*d);
     }
     return disp;
   } else {
     // with damping
-    damp = tang_toennies(b*d, 6, g);
+    damp = tang_toennies(b*d, power, g);
     // compute the energy
-    disp = d*d;
-    disp *= disp*disp;
-    disp = -c6/disp;
+    disp = 1.0;
+    for (j=0;j<power;j++) { disp *= d; }
+    disp = -cn/disp;
     if (g != NULL) {
-      *g = ((*g)*b-6.0/d*damp)*disp/d;
+      *g = ((*g)*b-power/d*damp)*disp/d;
     }
     return damp*disp;
   }

yaff/pes/pair_pot.h

@@ -124,12 +124,13 @@ double pair_fn_ljcross(void *pair_data, long center_index, long other_index, dou
 
 typedef struct {
   long nffatype;
+  long power;
   long *ffatype_ids;
-  double *c6_cross;
+  double *cn_cross;
   double *b_cross;
 } pair_data_dampdisp_type;
 
-void pair_data_dampdisp_init(pair_pot_type *pair_pot, long nffatype, long* ffatype_ids, double *c6_cross, double *b_cross);
+void pair_data_dampdisp_init(pair_pot_type *pair_pot, long nffatype, long power, long* ffatype_ids, double *cn_cross, double *b_cross);
 double pair_fn_dampdisp(void *pair_data, long center_index, long other_index, double d, double *delta, double *g, double *g_cart);
 
 

yaff/pes/pair_pot.pxd

@@ -62,7 +62,7 @@ cdef extern from "pair_pot.h":
 
     void pair_data_ljcross_init(pair_pot_type *pair_pot, long nffatype, long* ffatype_ids, double *eps_cross, double *sig_cross)
 
-    void pair_data_dampdisp_init(pair_pot_type *pair_pot, long nffatype, long* ffatype_ids, double *c6_cross, double *b_cross)
+    void pair_data_dampdisp_init(pair_pot_type *pair_pot, long nffatype, long power, long* ffatype_ids, double *cn_cross, double *b_cross)
 
     void pair_data_disp68bjdamp_init(pair_pot_type *pair_pot, long nffatype, long* ffaype_ids, double *c6_cross, double *c8_cross, double *R_cross, double c6_scale, double c8_scale, double bj_a, double bj_b)
     double pair_data_disp68bjdamp_get_c6_scale(pair_pot_type *pair_pot)

yaff/pes/test/test_generator.py

@@ -433,7 +433,7 @@ def test_generator_water32_dampdisp1():
     assert len(ff.parts) == 1
     part_pair_dampdisp = ff.part_pair_dampdisp
     # check parameters
-    c6_cross = part_pair_dampdisp.pair_pot.c6_cross
+    c6_cross = part_pair_dampdisp.pair_pot.cn_cross
     assert c6_cross.shape == (2,2)
     assert abs(c6_cross[0,0] - 1.9550248340e+01) < 1e-10
     assert abs(c6_cross[1,1] - 2.7982205915e+00) < 1e-10
@@ -458,7 +458,7 @@ def test_generator_water32_dampdisp2():
     assert len(ff.parts) == 1
     part_pair_dampdisp = ff.part_pair_dampdisp
     # check parameters
-    c6_cross = part_pair_dampdisp.pair_pot.c6_cross
+    c6_cross = part_pair_dampdisp.pair_pot.cn_cross
     assert abs(c6_cross[0,0] - 1.9550248340e+01) < 1e-10
     assert abs(c6_cross[0,1] - 6.4847208208e+00) < 1e-10
     assert abs(c6_cross[1,1] - 2.7982205915e+00) < 1e-10
@@ -479,7 +479,7 @@ def test_generator_glycine_dampdisp1():
     assert len(ff.parts) == 1
     part_pair_dampdisp = ff.part_pair_dampdisp
     # check parameters
-    c6_cross = part_pair_dampdisp.pair_pot.c6_cross
+    c6_cross = part_pair_dampdisp.pair_pot.cn_cross
     assert c6_cross.shape == (4, 4)
     assert abs(c6_cross[0,0] - 2.0121581791e+01) < 1e-10
     assert abs(c6_cross[1,1] - 2.5121581791e+01) < 1e-10
@@ -682,7 +682,7 @@ def test_generator_water32():
     part_ewald_cor = ff.part_ewald_cor
     part_ewald_neut = ff.part_ewald_neut
     # check dampdisp parameters
-    c6_cross = part_pair_dampdisp.pair_pot.c6_cross
+    c6_cross = part_pair_dampdisp.pair_pot.cn_cross
     assert abs(c6_cross[0,0] - 1.9550248340e+01) < 1e-10
     assert abs(c6_cross[1,1] - 2.7982205915e+00) < 1e-10
     assert (c6_cross == c6_cross.T).all()

yaff/pes/test/test_pair_pot.py

@@ -830,44 +830,44 @@ def test_pair_pot_ljcross_caffeine_9A():
     check_pair_pot_caffeine(system, nlist, scalings, part_pair, pair_fn, 1e-15)
 
 
-def get_part_caffeine_dampdisp_9A():
+def get_part_caffeine_dampdisp_9A(power=6):
     # Get a system and define scalings
     system = get_system_caffeine()
     nlist = NeighborList(system)
     scalings = Scalings(system, 0.0, 1.0, 1.0)
     # Initialize (very random) parameters
-    c6s = np.array([2.5, 27.0, 18.0, 13.0])
+    cns = np.array([2.5, 27.0, 18.0, 13.0])
     bs = np.array([2.5, 2.0, 0.0, 1.8])
     vols = np.array([5, 3, 4, 5])*angstrom**3
     # Allocate some arrays
     assert system.nffatype == 4
-    c6_cross = np.zeros((4, 4), float)
+    cn_cross = np.zeros((4, 4), float)
     b_cross = np.zeros((4, 4), float)
     # Construct the pair potential and part
-    pair_pot = PairPotDampDisp(system.ffatype_ids, c6_cross, b_cross, 9*angstrom, None, c6s, bs, vols)
+    pair_pot = PairPotDampDisp(system.ffatype_ids, cn_cross, b_cross, 9*angstrom, None, cns, bs, vols, power=power)
     part_pair = ForcePartPair(system, nlist, scalings, pair_pot)
     # The pair function
     def pair_fn(i0, i1, d):
-        c60 = c6s[system.ffatype_ids[i0]]
-        c61 = c6s[system.ffatype_ids[i1]]
+        cn0 = cns[system.ffatype_ids[i0]]
+        cn1 = cns[system.ffatype_ids[i1]]
         b0 = bs[system.ffatype_ids[i0]]
         b1 = bs[system.ffatype_ids[i1]]
         vol0 = vols[system.ffatype_ids[i0]]
         vol1 = vols[system.ffatype_ids[i1]]
         ratio = vol0/vol1
-        c6 = 2*c60*c61/(c60/ratio+c61*ratio)
+        cn = 2*cn0*cn1/(cn0/ratio+cn1*ratio)
         if b0 != 0 and b1 != 0:
             b = 0.5*(b0+b1)
             damp = 0
             fac = 1
-            for k in xrange(7):
+            for k in xrange(power+1):
                 damp += (b*d)**k/fac
                 fac *= k+1
             damp = 1 - np.exp(-b*d)*damp
-            return -c6/d**6*damp
+            return -cn/d**power*damp
         else:
             damp = 1
-            return -c6/d**6
+            return -cn/d**power
     return system, nlist, scalings, part_pair, pair_fn
 
 
@@ -876,6 +876,16 @@ def test_pair_pot_dampdisp_caffeine_9A():
     check_pair_pot_caffeine(system, nlist, scalings, part_pair, pair_fn, 1e-15)
 
 
+def test_pair_pot_dampdisp8_caffeine_9A():
+    system, nlist, scalings, part_pair, pair_fn = get_part_caffeine_dampdisp_9A(power=8)
+    check_pair_pot_caffeine(system, nlist, scalings, part_pair, pair_fn, 1e-15)
+
+
+def test_pair_pot_dampdisp10_caffeine_9A():
+    system, nlist, scalings, part_pair, pair_fn = get_part_caffeine_dampdisp_9A(power=10)
+    check_pair_pot_caffeine(system, nlist, scalings, part_pair, pair_fn, 1e-15)
+
+
 def get_part_caffeine_ei1_10A():
     # Get a system and define scalings
     system = get_system_caffeine()
@@ -1273,7 +1283,7 @@ def test_bks_isfinite():
     system = get_system_quartz()
     fn_pars = context.get_fn('test/parameters_bks.txt')
     ff = ForceField.generate(system, fn_pars)
-    assert np.isfinite(ff.part_pair_dampdisp.pair_pot.c6_cross).all()
+    assert np.isfinite(ff.part_pair_dampdisp.pair_pot.cn_cross).all()
     assert np.isfinite(ff.part_pair_dampdisp.pair_pot.b_cross).all()
     ff.compute()
     assert np.isfinite(ff.part_pair_exprep.energy)