diff --git a/CHANGES/v2.4.txt b/CHANGES/v2.4.txt
index 3a2f658face8feab1f761ac83f70c1c929157602..43b058aaa4130bd189f3fc24314d6e2683ba0b94 100644
--- a/CHANGES/v2.4.txt
+++ b/CHANGES/v2.4.txt
@@ -34,6 +34,7 @@ Changes from version 2.3 which are relevant for users:
trajectory with replica suffix is not found the driver will look for a trajectory without the replica suffix.
- Internal molfile implementation has been updated to VMD 1.9.3.
- Examples in the documentation now have syntax highlighting and links to the documentation of used actions.
+ - \ref COORDINATIONNUMBER : Added option to have pairwise distance moments of coordination number in the multicolvar module
Changes from version 2.3 which are relevant for developers:
- A few fixes has been made to improve exception safety. Although we still cannot declare
diff --git a/regtest/multicolvar/rt-coordination-powers-multi/Makefile b/regtest/multicolvar/rt-coordination-powers-multi/Makefile
new file mode 100644
index 0000000000000000000000000000000000000000..3703b27cea227aa053fb6d1d73f861e4384dbcee
--- /dev/null
+++ b/regtest/multicolvar/rt-coordination-powers-multi/Makefile
@@ -0,0 +1 @@
+include ../../scripts/test.make
diff --git a/regtest/multicolvar/rt-coordination-powers-multi/cn_out.reference b/regtest/multicolvar/rt-coordination-powers-multi/cn_out.reference
new file mode 100644
index 0000000000000000000000000000000000000000..fbebd72b92dd5db462c2e83cda67778ee4388e58
--- /dev/null
+++ b/regtest/multicolvar/rt-coordination-powers-multi/cn_out.reference
@@ -0,0 +1,6 @@
+#! FIELDS time cn0.mean cn1.mean cn2.mean
+ 0.000000 1.171228 1.453038 1.843475
+ 0.005000 1.166065 1.441594 1.826868
+ 0.010000 1.172897 1.440122 1.816671
+ 0.015000 1.158162 1.412160 1.776137
+ 0.020000 1.105913 1.356438 1.719031
diff --git a/regtest/multicolvar/rt-coordination-powers-multi/config b/regtest/multicolvar/rt-coordination-powers-multi/config
new file mode 100644
index 0000000000000000000000000000000000000000..ba8b778571e71d71dd136e9098170c96b46493fe
--- /dev/null
+++ b/regtest/multicolvar/rt-coordination-powers-multi/config
@@ -0,0 +1,4 @@
+type=driver
+# this is to test a different name
+arg="--plumed plumed.dat --timestep 0.005 --ixyz trajectory.xyz --dump-forces forces --dump-forces-fmt=%8.4f"
+extra_files="../../trajectories/trajectory.xyz"
diff --git a/regtest/multicolvar/rt-coordination-powers-multi/plumed.dat b/regtest/multicolvar/rt-coordination-powers-multi/plumed.dat
new file mode 100644
index 0000000000000000000000000000000000000000..e7c2dcce2f1f69d33d26c1c9134626722f2aa560
--- /dev/null
+++ b/regtest/multicolvar/rt-coordination-powers-multi/plumed.dat
@@ -0,0 +1,5 @@
+cn0: COORDINATIONNUMBER SPECIES=1-10 SWITCH={RATIONAL R_0=1.0 D_MAX=8} MEAN
+cn1: COORDINATIONNUMBER SPECIES=1-10 SWITCH={RATIONAL R_0=1.0 D_MAX=8} R_POWER=1 MEAN
+cn2: COORDINATIONNUMBER SPECIES=1-10 SWITCH={RATIONAL R_0=1.0 D_MAX=8} R_POWER=2 MEAN
+PRINT ARG=cn0.mean,cn1.mean,cn2.mean STRIDE=1 FILE=cn_out
+
diff --git a/regtest/multicolvar/rt-coordination-powers/Makefile b/regtest/multicolvar/rt-coordination-powers/Makefile
new file mode 100644
index 0000000000000000000000000000000000000000..3703b27cea227aa053fb6d1d73f861e4384dbcee
--- /dev/null
+++ b/regtest/multicolvar/rt-coordination-powers/Makefile
@@ -0,0 +1 @@
+include ../../scripts/test.make
diff --git a/regtest/multicolvar/rt-coordination-powers/config b/regtest/multicolvar/rt-coordination-powers/config
new file mode 100644
index 0000000000000000000000000000000000000000..ba8b778571e71d71dd136e9098170c96b46493fe
--- /dev/null
+++ b/regtest/multicolvar/rt-coordination-powers/config
@@ -0,0 +1,4 @@
+type=driver
+# this is to test a different name
+arg="--plumed plumed.dat --timestep 0.005 --ixyz trajectory.xyz --dump-forces forces --dump-forces-fmt=%8.4f"
+extra_files="../../trajectories/trajectory.xyz"
diff --git a/regtest/multicolvar/rt-coordination-powers/derivatives.reference b/regtest/multicolvar/rt-coordination-powers/derivatives.reference
new file mode 100644
index 0000000000000000000000000000000000000000..d6fa5fe2f4f4972267275f737316149ac83fcf41
--- /dev/null
+++ b/regtest/multicolvar/rt-coordination-powers/derivatives.reference
@@ -0,0 +1,196 @@
+#! FIELDS time parameter c1.mean c1num.mean
+ 0.000000 0 0.3371 0.3371
+ 0.000000 1 0.2309 0.2309
+ 0.000000 2 0.1160 0.1160
+ 0.000000 3 -0.3510 -0.3510
+ 0.000000 4 0.3400 0.3399
+ 0.000000 5 0.0247 0.0247
+ 0.000000 6 -0.2278 -0.2278
+ 0.000000 7 -0.3331 -0.3331
+ 0.000000 8 0.1717 0.1717
+ 0.000000 9 0.3075 0.3075
+ 0.000000 10 -0.3590 -0.3590
+ 0.000000 11 0.0586 0.0586
+ 0.000000 12 0.3427 0.3427
+ 0.000000 13 0.3101 0.3101
+ 0.000000 14 -0.0220 -0.0220
+ 0.000000 15 -0.3299 -0.3299
+ 0.000000 16 0.2127 0.2127
+ 0.000000 17 -0.1563 -0.1563
+ 0.000000 18 -0.3289 -0.3289
+ 0.000000 19 -0.3329 -0.3329
+ 0.000000 20 -0.0762 -0.0762
+ 0.000000 21 0.2298 0.2298
+ 0.000000 22 -0.3171 -0.3171
+ 0.000000 23 -0.1589 -0.1589
+ 0.000000 24 0.2379 0.2379
+ 0.000000 25 0.1246 0.1246
+ 0.000000 26 -0.0916 -0.0916
+ 0.000000 27 -0.2175 -0.2175
+ 0.000000 28 0.1238 0.1238
+ 0.000000 29 0.1338 0.1338
+ 0.000000 30 1.2685 1.2685
+ 0.000000 31 0.0008 0.0008
+ 0.000000 32 -0.0777 -0.0777
+ 0.000000 33 0.0008 0.0008
+ 0.000000 34 1.1232 1.1232
+ 0.000000 35 0.0031 0.0031
+ 0.000000 36 -0.0777 -0.0777
+ 0.000000 37 0.0031 0.0031
+ 0.000000 38 2.5442 2.5442
+ 0.005000 0 0.3287 0.3287
+ 0.005000 1 0.2323 0.2323
+ 0.005000 2 0.1100 0.1100
+ 0.005000 3 -0.3587 -0.3587
+ 0.005000 4 0.3366 0.3366
+ 0.005000 5 0.0181 0.0181
+ 0.005000 6 -0.2232 -0.2232
+ 0.005000 7 -0.3296 -0.3296
+ 0.005000 8 0.1699 0.1699
+ 0.005000 9 0.2907 0.2907
+ 0.005000 10 -0.3585 -0.3585
+ 0.005000 11 0.0690 0.0690
+ 0.005000 12 0.3444 0.3444
+ 0.005000 13 0.2895 0.2895
+ 0.005000 14 -0.0237 -0.0237
+ 0.005000 15 -0.3163 -0.3163
+ 0.005000 16 0.1948 0.1948
+ 0.005000 17 -0.1721 -0.1721
+ 0.005000 18 -0.3260 -0.3260
+ 0.005000 19 -0.3198 -0.3198
+ 0.005000 20 -0.0941 -0.0941
+ 0.005000 21 0.2273 0.2273
+ 0.005000 22 -0.2878 -0.2878
+ 0.005000 23 -0.1523 -0.1523
+ 0.005000 24 0.2355 0.2355
+ 0.005000 25 0.1149 0.1149
+ 0.005000 26 -0.0652 -0.0652
+ 0.005000 27 -0.2023 -0.2023
+ 0.005000 28 0.1276 0.1276
+ 0.005000 29 0.1406 0.1406
+ 0.005000 30 1.2982 1.2982
+ 0.005000 31 0.0073 0.0073
+ 0.005000 32 -0.0442 -0.0442
+ 0.005000 33 0.0073 0.0073
+ 0.005000 34 1.0615 1.0615
+ 0.005000 35 0.0049 0.0049
+ 0.005000 36 -0.0442 -0.0442
+ 0.005000 37 0.0049 0.0049
+ 0.005000 38 2.4945 2.4945
+ 0.010000 0 0.3057 0.3057
+ 0.010000 1 0.2219 0.2219
+ 0.010000 2 0.1122 0.1122
+ 0.010000 3 -0.3610 -0.3610
+ 0.010000 4 0.3207 0.3207
+ 0.010000 5 0.0249 0.0249
+ 0.010000 6 -0.2431 -0.2431
+ 0.010000 7 -0.3226 -0.3226
+ 0.010000 8 0.1861 0.1861
+ 0.010000 9 0.3110 0.3110
+ 0.010000 10 -0.3230 -0.3230
+ 0.010000 11 0.0536 0.0536
+ 0.010000 12 0.3321 0.3321
+ 0.010000 13 0.2459 0.2459
+ 0.010000 14 -0.0204 -0.0204
+ 0.010000 15 -0.3035 -0.3035
+ 0.010000 16 0.1750 0.1750
+ 0.010000 17 -0.1894 -0.1894
+ 0.010000 18 -0.3215 -0.3215
+ 0.010000 19 -0.3011 -0.3011
+ 0.010000 20 -0.0991 -0.0991
+ 0.010000 21 0.2342 0.2342
+ 0.010000 22 -0.2677 -0.2677
+ 0.010000 23 -0.1511 -0.1511
+ 0.010000 24 0.2325 0.2325
+ 0.010000 25 0.1200 0.1200
+ 0.010000 26 -0.0458 -0.0458
+ 0.010000 27 -0.1864 -0.1864
+ 0.010000 28 0.1309 0.1309
+ 0.010000 29 0.1290 0.1290
+ 0.010000 30 1.3545 1.3545
+ 0.010000 31 -0.0010 -0.0010
+ 0.010000 32 -0.0324 -0.0324
+ 0.010000 33 -0.0010 -0.0010
+ 0.010000 34 0.9468 0.9468
+ 0.010000 35 0.0117 0.0117
+ 0.010000 36 -0.0324 -0.0324
+ 0.010000 37 0.0117 0.0117
+ 0.010000 38 2.4481 2.4481
+ 0.015000 0 0.2784 0.2784
+ 0.015000 1 0.2040 0.2040
+ 0.015000 2 0.0908 0.0908
+ 0.015000 3 -0.3455 -0.3455
+ 0.015000 4 0.2867 0.2867
+ 0.015000 5 0.0230 0.0230
+ 0.015000 6 -0.2596 -0.2596
+ 0.015000 7 -0.3044 -0.3044
+ 0.015000 8 0.2036 0.2036
+ 0.015000 9 0.3216 0.3216
+ 0.015000 10 -0.2752 -0.2752
+ 0.015000 11 0.0411 0.0411
+ 0.015000 12 0.3142 0.3142
+ 0.015000 13 0.1924 0.1924
+ 0.015000 14 -0.0163 -0.0163
+ 0.015000 15 -0.2939 -0.2939
+ 0.015000 16 0.1585 0.1585
+ 0.015000 17 -0.2067 -0.2067
+ 0.015000 18 -0.3090 -0.3090
+ 0.015000 19 -0.2674 -0.2674
+ 0.015000 20 -0.0777 -0.0777
+ 0.015000 21 0.2317 0.2317
+ 0.015000 22 -0.2560 -0.2560
+ 0.015000 23 -0.1327 -0.1327
+ 0.015000 24 0.2310 0.2310
+ 0.015000 25 0.1251 0.1251
+ 0.015000 26 -0.0350 -0.0350
+ 0.015000 27 -0.1689 -0.1689
+ 0.015000 28 0.1362 0.1362
+ 0.015000 29 0.1099 0.1099
+ 0.015000 30 1.3936 1.3936
+ 0.015000 31 -0.0040 -0.0040
+ 0.015000 32 -0.0365 -0.0365
+ 0.015000 33 -0.0040 -0.0040
+ 0.015000 34 0.8269 0.8269
+ 0.015000 35 -0.0103 -0.0103
+ 0.015000 36 -0.0365 -0.0365
+ 0.015000 37 -0.0103 -0.0103
+ 0.015000 38 2.3506 2.3506
+ 0.020000 0 0.2666 0.2666
+ 0.020000 1 0.1868 0.1868
+ 0.020000 2 0.0703 0.0703
+ 0.020000 3 -0.3254 -0.3254
+ 0.020000 4 0.2691 0.2691
+ 0.020000 5 0.0102 0.0102
+ 0.020000 6 -0.2645 -0.2645
+ 0.020000 7 -0.2835 -0.2835
+ 0.020000 8 0.2141 0.2141
+ 0.020000 9 0.3166 0.3166
+ 0.020000 10 -0.2635 -0.2635
+ 0.020000 11 0.0653 0.0653
+ 0.020000 12 0.3024 0.3024
+ 0.020000 13 0.1976 0.1976
+ 0.020000 14 -0.0355 -0.0355
+ 0.020000 15 -0.2931 -0.2931
+ 0.020000 16 0.1405 0.1405
+ 0.020000 17 -0.2154 -0.2154
+ 0.020000 18 -0.2976 -0.2976
+ 0.020000 19 -0.2589 -0.2589
+ 0.020000 20 -0.0579 -0.0579
+ 0.020000 21 0.2223 0.2223
+ 0.020000 22 -0.2594 -0.2594
+ 0.020000 23 -0.1251 -0.1251
+ 0.020000 24 0.2263 0.2263
+ 0.020000 25 0.1348 0.1348
+ 0.020000 26 -0.0242 -0.0242
+ 0.020000 27 -0.1535 -0.1535
+ 0.020000 28 0.1365 0.1365
+ 0.020000 29 0.0982 0.0982
+ 0.020000 30 1.3941 1.3941
+ 0.020000 31 0.0116 0.0116
+ 0.020000 32 -0.0319 -0.0319
+ 0.020000 33 0.0116 0.0116
+ 0.020000 34 0.8065 0.8065
+ 0.020000 35 -0.0388 -0.0388
+ 0.020000 36 -0.0319 -0.0319
+ 0.020000 37 -0.0388 -0.0388
+ 0.020000 38 2.2835 2.2835
diff --git a/regtest/multicolvar/rt-coordination-powers/plumed.dat b/regtest/multicolvar/rt-coordination-powers/plumed.dat
new file mode 100644
index 0000000000000000000000000000000000000000..ec3e4821ca5de2b98285b66e9ea0dad56959f3c9
--- /dev/null
+++ b/regtest/multicolvar/rt-coordination-powers/plumed.dat
@@ -0,0 +1,3 @@
+COORDINATIONNUMBER SPECIES=1-10 R_0=1 R_POWER=2 LABEL=c1 MEAN
+COORDINATIONNUMBER SPECIES=1-10 R_0=1 NUMERICAL_DERIVATIVES R_POWER=2 LABEL=c1num MEAN
+DUMPDERIVATIVES ARG=c1.*,c1num.* STRIDE=1 FILE=derivatives FMT=%8.4f
diff --git a/src/multicolvar/CoordinationNumbers.cpp b/src/multicolvar/CoordinationNumbers.cpp
index 38c26d348b9486b9c000c68e261a6395af081b3c..e280931ed3a2a0d9e9dbd5f6b5d4b6ac7b5664b4 100644
--- a/src/multicolvar/CoordinationNumbers.cpp
+++ b/src/multicolvar/CoordinationNumbers.cpp
@@ -44,6 +44,14 @@ To make the calculation of coordination numbers differentiable the following fun
s = \frac{ 1 - \left(\frac{r-d_0}{r_0}\right)^n } { 1 - \left(\frac{r-d_0}{r_0}\right)^m }
\f]
+If R_POWER is set, this will use the product of pairwise distance
+raised to the R_POWER with the coordination number function defined
+above. This was used in White and Voth \cite white2014efficient as a
+way of indirectly biasing radial distribution functions. Note that in
+that reference this function is referred to as moments of coordination
+number, but here we call them powers to distinguish from the existing
+MOMENTS keyword of Multicolvars.
+
\par Examples
The following input tells plumed to calculate the coordination numbers of atoms 1-100 with themselves.
@@ -59,14 +67,23 @@ number of coordination numbers more than 6 is then computed.
COORDINATIONNUMBER SPECIESA=101-110 SPECIESB=1-100 R_0=3.0 MORE_THAN={RATIONAL R_0=6.0 NN=6 MM=12 D_0=0}
\endplumedfile
+The following input tells plumed to calculate the mean coordination number of all atoms with themselves
+and its powers. An explicit cutoff is set for each of 8.
+\plumedfile
+cn0: COORDINATIONNUMBER SPECIES=1-10 SWITCH={RATIONAL R_0=1.0 D_MAX=8} MEAN
+cn1: COORDINATIONNUMBER SPECIES=1-10 SWITCH={RATIONAL R_0=1.0 D_MAX=8} R_POWER=1 MEAN
+cn2: COORDINATIONNUMBER SPECIES=1-10 SWITCH={RATIONAL R_0=1.0 D_MAX=8} R_POWER=2 MEAN
+PRINT ARG=cn0.mean,cn1.mean,cn2.mean STRIDE=1 FILE=cn_out
+\endplumedfile
+
*/
//+ENDPLUMEDOC
class CoordinationNumbers : public MultiColvarBase {
private:
-// double nl_cut;
double rcut2;
+ int r_power;
SwitchingFunction switchingFunction;
public:
static void registerKeywords( Keywords& keys );
@@ -86,9 +103,11 @@ void CoordinationNumbers::registerKeywords( Keywords& keys ) {
keys.add("compulsory","MM","0","The m parameter of the switching function; 0 implies 2*NN");
keys.add("compulsory","D_0","0.0","The d_0 parameter of the switching function");
keys.add("compulsory","R_0","The r_0 parameter of the switching function");
+ keys.add("optional","R_POWER","Multiply the coordination number function by a power of r, "
+ "as done in White and Voth (see note above, default: no)");
keys.add("optional","SWITCH","This keyword is used if you want to employ an alternative to the continuous swiching function defined above. "
- "The following provides information on the \\ref switchingfunction that are available. "
- "When this keyword is present you no longer need the NN, MM, D_0 and R_0 keywords.");
+ "The following provides information on the \\ref switchingfunction that are available. "
+ "When this keyword is present you no longer need the NN, MM, D_0 and R_0 keywords.");
// Use actionWithDistributionKeywords
keys.use("MEAN"); keys.use("MORE_THAN"); keys.use("LESS_THAN"); keys.use("MAX");
keys.use("MIN"); keys.use("BETWEEN"); keys.use("HISTOGRAM"); keys.use("MOMENTS");
@@ -97,8 +116,10 @@ void CoordinationNumbers::registerKeywords( Keywords& keys ) {
CoordinationNumbers::CoordinationNumbers(const ActionOptions&ao):
Action(ao),
- MultiColvarBase(ao)
+ MultiColvarBase(ao),
+ r_power(0)
{
+
// Read in the switching function
std::string sw, errors; parse("SWITCH",sw);
if(sw.length()>0) {
@@ -110,11 +131,25 @@ CoordinationNumbers::CoordinationNumbers(const ActionOptions&ao):
parse("R_0",r_0); parse("D_0",d_0);
if( r_0<0.0 ) error("you must set a value for R_0");
switchingFunction.set(nn,mm,r_0,d_0);
+
}
log.printf(" coordination of central atom and those within %s\n",( switchingFunction.description() ).c_str() );
+
+ //get cutoff of switching function
+ double rcut = switchingFunction.get_dmax();
+
+ //parse power
+ parse("R_POWER", r_power);
+ if(r_power > 0) {
+ log.printf(" Multiplying switching function by r^%d\n", r_power);
+ double offset = switchingFunction.calculate(rcut*0.9999, rcut2) * pow(rcut*0.9999, r_power);
+ log.printf(" You will have a discontinuous jump of %f to 0 near the cutoff of your switching function. "
+ "Consider setting D_MAX or reducing R_POWER if this is large\n", offset);
+ }
+
// Set the link cell cutoff
- setLinkCellCutoff( switchingFunction.get_dmax() );
- rcut2 = switchingFunction.get_dmax()*switchingFunction.get_dmax();
+ setLinkCellCutoff( rcut );
+ rcut2 = rcut * rcut;
// And setup the ActionWithVessel
std::vector<AtomNumber> all_atoms; setupMultiColvarBase( all_atoms ); checkRead();
@@ -122,15 +157,23 @@ CoordinationNumbers::CoordinationNumbers(const ActionOptions&ao):
double CoordinationNumbers::compute( const unsigned& tindex, AtomValuePack& myatoms ) const {
// Calculate the coordination number
- double dfunc, d2, sw;
+ double dfunc, d2, sw, d, raised;
for(unsigned i=1; i<myatoms.getNumberOfAtoms(); ++i) {
Vector& distance=myatoms.getPosition(i);
if ( (d2=distance[0]*distance[0])<rcut2 &&
- (d2+=distance[1]*distance[1])<rcut2 &&
- (d2+=distance[2]*distance[2])<rcut2) {
+ (d2+=distance[1]*distance[1])<rcut2 &&
+ (d2+=distance[2]*distance[2])<rcut2) {
sw = switchingFunction.calculateSqr( d2, dfunc );
- accumulateSymmetryFunction( 1, i, sw, (dfunc)*distance, (-dfunc)*Tensor(distance,distance), myatoms );
+ if(r_power > 0) {
+ d = sqrt(d2); raised = pow( d, r_power - 1 );
+ accumulateSymmetryFunction( 1, i, sw * raised * d,
+ (dfunc * d * raised + sw * r_power) * distance,
+ (-dfunc * d * raised - sw * r_power) * Tensor(distance, distance),
+ myatoms );
+ } else {
+ accumulateSymmetryFunction( 1, i, sw, (dfunc)*distance, (-dfunc)*Tensor(distance,distance), myatoms );
+ }
}
}
@@ -139,4 +182,3 @@ double CoordinationNumbers::compute( const unsigned& tindex, AtomValuePack& myat
}
}
-