diff --git a/regtest/basic/rt43/COLVAR.reference b/regtest/basic/rt43/COLVAR.reference new file mode 100644 index 0000000000000000000000000000000000000000..25d32c01d6d3402e1c472b8d716915d44bfd11b9 --- /dev/null +++ b/regtest/basic/rt43/COLVAR.reference @@ -0,0 +1,6 @@ +#! FIELDS time rmsd0 rmsd1 rmsd2 rmsd3 rmsd4 stat.sqdevsum stat.corr stat.slope stat.intercept statn.sqdevsum statn.corr statn.slope statn.intercept statsq.sqdevsum + 0.000000 1.496 1.496 1.337 1.336 2.329 0.000 1.000 1.000 0.000 0.000 1.000 1.000 0.000 0.000 + 0.050000 1.502 1.503 1.347 1.347 2.332 0.000 1.000 1.007 -0.019 0.000 1.000 1.007 -0.019 0.000 + 0.100000 1.503 1.504 1.351 1.350 2.324 0.001 1.000 1.018 -0.037 0.001 1.000 1.018 -0.037 0.001 + 0.150000 1.502 1.502 1.350 1.350 2.308 0.001 1.000 1.035 -0.060 0.001 1.000 1.035 -0.060 0.001 + 0.200000 1.496 1.496 1.347 1.346 2.278 0.003 1.000 1.065 -0.098 0.003 1.000 1.065 -0.098 0.003 diff --git a/regtest/basic/rt43/Makefile b/regtest/basic/rt43/Makefile new file mode 100644 index 0000000000000000000000000000000000000000..3703b27cea227aa053fb6d1d73f861e4384dbcee --- /dev/null +++ b/regtest/basic/rt43/Makefile @@ -0,0 +1 @@ +include ../../scripts/test.make diff --git a/regtest/basic/rt43/config b/regtest/basic/rt43/config new file mode 100644 index 0000000000000000000000000000000000000000..4804b7f2a86c1dcb06f52d08f8dff331e318192f --- /dev/null +++ b/regtest/basic/rt43/config @@ -0,0 +1,4 @@ +type=driver +# this is to test a different name +arg="--plumed plumed.dat --trajectory-stride 10 --timestep 0.005 --ixyz trajectory.xyz --dump-forces forces --dump-forces-fmt=%8.4f" +extra_files="../../trajectories/trajectory.xyz" diff --git a/regtest/basic/rt43/deriv.reference b/regtest/basic/rt43/deriv.reference new file mode 100644 index 0000000000000000000000000000000000000000..ca14d348bcf00244134ca724bd16a3f488699479 --- /dev/null +++ b/regtest/basic/rt43/deriv.reference @@ -0,0 +1,26 @@ +#! FIELDS time parameter stat.sqdevsum stat.corr stat.slope stat.intercept statn.sqdevsum statn.corr statn.slope statn.intercept + 0.000000 0 -0.0004 0.0002 0.1489 -0.4381 -0.0004 0.0002 0.1489 -0.4381 + 0.000000 1 0.0003 -0.0003 0.1480 -0.4365 0.0003 -0.0003 0.1480 -0.4365 + 0.000000 2 -0.0008 0.0005 0.3792 -0.8063 -0.0008 0.0005 0.3792 -0.8063 + 0.000000 3 0.0004 -0.0004 0.3789 -0.8057 0.0004 -0.0004 0.3789 -0.8057 + 0.000000 4 0.0001 0.0000 -1.0550 1.4867 0.0001 0.0000 -1.0550 1.4867 + 0.050000 0 0.0128 0.0026 0.1559 -0.4518 0.0128 0.0026 0.1559 -0.4518 + 0.050000 1 0.0136 0.0021 0.1548 -0.4501 0.0136 0.0021 0.1548 -0.4501 + 0.050000 2 0.0206 -0.0015 0.3806 -0.8127 0.0206 -0.0015 0.3806 -0.8127 + 0.050000 3 0.0219 -0.0024 0.3802 -0.8121 0.0219 -0.0024 0.3802 -0.8121 + 0.050000 4 0.0059 -0.0007 -1.0715 1.5198 0.0059 -0.0007 -1.0715 1.5198 + 0.100000 0 0.0149 0.0030 0.1602 -0.4609 0.0149 0.0030 0.1602 -0.4609 + 0.100000 1 0.0157 0.0024 0.1590 -0.4591 0.0157 0.0024 0.1590 -0.4590 + 0.100000 2 0.0270 -0.0018 0.3887 -0.8280 0.0270 -0.0018 0.3887 -0.8280 + 0.100000 3 0.0283 -0.0027 0.3883 -0.8274 0.0283 -0.0027 0.3883 -0.8274 + 0.100000 4 -0.0099 -0.0009 -1.0961 1.5571 -0.0099 -0.0009 -1.0961 1.5571 + 0.150000 0 0.0115 0.0023 0.1640 -0.4699 0.0115 0.0023 0.1640 -0.4699 + 0.150000 1 0.0122 0.0017 0.1628 -0.4680 0.0122 0.0017 0.1628 -0.4680 + 0.150000 2 0.0269 -0.0012 0.4030 -0.8529 0.0269 -0.0012 0.4030 -0.8529 + 0.150000 3 0.0282 -0.0022 0.4026 -0.8523 0.0282 -0.0022 0.4026 -0.8523 + 0.150000 4 -0.0415 -0.0006 -1.1325 1.6078 -0.0415 -0.0006 -1.1325 1.6078 + 0.200000 0 -0.0005 0.0006 0.1700 -0.4839 -0.0005 0.0006 0.1700 -0.4839 + 0.200000 1 0.0002 0.0000 0.1687 -0.4818 0.0002 0.0000 0.1687 -0.4818 + 0.200000 2 0.0196 0.0002 0.4300 -0.8978 0.0196 0.0002 0.4300 -0.8978 + 0.200000 3 0.0209 -0.0008 0.4296 -0.8972 0.0209 -0.0008 0.4296 -0.8972 + 0.200000 4 -0.1020 -0.0001 -1.1982 1.6952 -0.1020 -0.0001 -1.1982 1.6952 diff --git a/regtest/basic/rt43/plumed.dat b/regtest/basic/rt43/plumed.dat new file mode 100644 index 0000000000000000000000000000000000000000..3c656ff0307babff2546adebb9443886f44bd43e --- /dev/null +++ b/regtest/basic/rt43/plumed.dat @@ -0,0 +1,22 @@ +rmsd0: RMSD TYPE=OPTIMAL REFERENCE=test0.pdb +rmsd1: RMSD TYPE=OPTIMAL REFERENCE=test1.pdb +rmsd2: RMSD TYPE=OPTIMAL REFERENCE=test2.pdb +rmsd3: RMSD TYPE=OPTIMAL REFERENCE=test3.pdb +rmsd4: RMSD TYPE=OPTIMAL REFERENCE=test4.pdb + +stat: STATS ARG=rmsd0,rmsd1,rmsd2,rmsd3,rmsd4 PARAMETERS=1.496,1.496,1.337,1.336,2.329 +statn: STATS ARG=rmsd0,rmsd1,rmsd2,rmsd3,rmsd4 PARAMETERS=1.496,1.496,1.337,1.336,2.329 NUMERICAL_DERIVATIVES +statsq: STATS ARG=rmsd0,rmsd1,rmsd2,rmsd3,rmsd4 PARAMETERS=1.496,1.496,1.337,1.336,2.329 SQDEVSUM + +DUMPDERIVATIVES ARG=stat.*,statn.* STRIDE=1 FILE=deriv FMT=%8.4f + +PRINT ... + STRIDE=1 + ARG=* + FILE=COLVAR FMT=%6.3f +... PRINT + + +ENDPLUMED + +text here should be ignored diff --git a/regtest/basic/rt43/test0.pdb b/regtest/basic/rt43/test0.pdb new file mode 100644 index 0000000000000000000000000000000000000000..eb3727de108b579ffd85bf6b03cc9c8e4181accd --- /dev/null +++ b/regtest/basic/rt43/test0.pdb @@ -0,0 +1,15 @@ +ATOM 2 O ALA 2 -0.926 -2.447 -0.497 1.00 1.00 DIA O +ATOM 4 HNT ALA 2 0.533 -0.396 1.184 1.00 1.00 DIA H +ATOM 6 HT1 ALA 2 -0.216 -2.590 1.371 1.00 1.00 DIA H +ATOM 7 HT2 ALA 2 -0.309 -1.255 2.315 1.00 1.00 DIA H +ATOM 8 HT3 ALA 2 -1.480 -1.560 1.212 1.00 1.00 DIA H +ATOM 9 CAY ALA 2 -0.096 2.144 -0.669 1.00 1.00 DIA C +ATOM 10 HY1 ALA 2 0.871 2.385 -0.588 1.00 1.00 DIA H +ATOM 12 HY3 ALA 2 -0.520 2.679 -1.400 1.00 1.00 DIA H +ATOM 14 OY ALA 2 -1.139 0.931 -0.973 1.00 1.00 DIA O +ATOM 16 HN ALA 2 1.713 1.021 -0.873 1.00 1.00 DIA H +ATOM 18 HA ALA 2 0.099 -0.774 -2.218 1.00 1.00 DIA H +ATOM 19 CB ALA 2 2.063 -1.223 -1.276 1.00 1.00 DIA C +ATOM 20 HB1 ALA 2 2.670 -0.716 -2.057 1.00 1.00 DIA H +ATOM 21 HB2 ALA 2 2.556 -1.051 -0.295 1.00 1.00 DIA H +ATOM 22 HB3 ALA 2 2.070 -2.314 -1.490 1.00 1.00 DIA H diff --git a/regtest/basic/rt43/test1.pdb b/regtest/basic/rt43/test1.pdb new file mode 100644 index 0000000000000000000000000000000000000000..924aa4b057f24d6467f39ff489e249b06faf8de8 --- /dev/null +++ b/regtest/basic/rt43/test1.pdb @@ -0,0 +1,15 @@ +ATOM 2 O ALA 2 -0.926 -2.447 -0.497 1.00 1.01 DIA O +ATOM 4 HNT ALA 2 0.533 -0.396 1.184 1.00 1.00 DIA H +ATOM 6 HT1 ALA 2 -0.216 -2.590 1.371 1.00 1.00 DIA H +ATOM 7 HT2 ALA 2 -0.309 -1.255 2.315 1.00 1.00 DIA H +ATOM 8 HT3 ALA 2 -1.480 -1.560 1.212 1.00 1.00 DIA H +ATOM 9 CAY ALA 2 -0.096 2.144 -0.669 1.00 1.00 DIA C +ATOM 10 HY1 ALA 2 0.871 2.385 -0.588 1.00 1.00 DIA H +ATOM 12 HY3 ALA 2 -0.520 2.679 -1.400 1.00 1.00 DIA H +ATOM 14 OY ALA 2 -1.139 0.931 -0.973 1.00 1.00 DIA O +ATOM 16 HN ALA 2 1.713 1.021 -0.873 1.00 1.00 DIA H +ATOM 18 HA ALA 2 0.099 -0.774 -2.218 1.00 1.00 DIA H +ATOM 19 CB ALA 2 2.063 -1.223 -1.276 1.00 0.99 DIA C +ATOM 20 HB1 ALA 2 2.670 -0.716 -2.057 1.00 1.00 DIA H +ATOM 21 HB2 ALA 2 2.556 -1.051 -0.295 1.00 1.00 DIA H +ATOM 22 HB3 ALA 2 2.070 -2.314 -1.490 1.00 1.00 DIA H diff --git a/regtest/basic/rt43/test2.pdb b/regtest/basic/rt43/test2.pdb new file mode 100644 index 0000000000000000000000000000000000000000..f2585025b36c688adfab6aee461070d572635f66 --- /dev/null +++ b/regtest/basic/rt43/test2.pdb @@ -0,0 +1,15 @@ +ATOM 2 O ALA 2 -0.926 -2.447 -0.497 1.10 1.10 DIA O +ATOM 4 HNT ALA 2 0.533 -0.396 1.184 1.00 1.00 DIA H +ATOM 6 HT1 ALA 2 -0.216 -2.590 1.371 1.00 1.00 DIA H +ATOM 7 HT2 ALA 2 -0.309 -1.255 2.315 1.00 1.00 DIA H +ATOM 8 HT3 ALA 2 -1.480 -1.560 1.212 5.00 5.00 DIA H +ATOM 9 CAY ALA 2 -0.096 2.144 -0.669 1.00 1.00 DIA C +ATOM 10 HY1 ALA 2 0.871 2.385 -0.588 1.00 1.00 DIA H +ATOM 12 HY3 ALA 2 -0.520 2.679 -1.400 1.00 1.00 DIA H +ATOM 14 OY ALA 2 -1.139 0.931 -0.973 0.00 0.00 DIA O +ATOM 16 HN ALA 2 1.713 1.021 -0.873 1.00 1.00 DIA H +ATOM 18 HA ALA 2 0.099 -0.774 -2.218 0.00 0.00 DIA H +ATOM 19 CB ALA 2 2.063 -1.223 -1.276 1.00 1.00 DIA C +ATOM 20 HB1 ALA 2 2.670 -0.716 -2.057 1.00 1.00 DIA H +ATOM 21 HB2 ALA 2 2.556 -1.051 -0.295 1.00 1.00 DIA H +ATOM 22 HB3 ALA 2 2.070 -2.314 -1.490 1.00 1.00 DIA H diff --git a/regtest/basic/rt43/test3.pdb b/regtest/basic/rt43/test3.pdb new file mode 100644 index 0000000000000000000000000000000000000000..4b04045e2a4baed871505c0e074d29d0065ab0fb --- /dev/null +++ b/regtest/basic/rt43/test3.pdb @@ -0,0 +1,15 @@ +ATOM 2 O ALA 2 -0.926 -2.447 -0.497 1.10 1.10 DIA O +ATOM 4 HNT ALA 2 0.533 -0.396 1.184 1.00 1.00 DIA H +ATOM 6 HT1 ALA 2 -0.216 -2.590 1.371 1.00 1.00 DIA H +ATOM 7 HT2 ALA 2 -0.309 -1.255 2.315 1.00 1.00 DIA H +ATOM 8 HT3 ALA 2 -1.480 -1.560 1.212 5.00 5.01 DIA H +ATOM 9 CAY ALA 2 -0.096 2.144 -0.669 1.00 1.00 DIA C +ATOM 10 HY1 ALA 2 0.871 2.385 -0.588 1.00 1.00 DIA H +ATOM 12 HY3 ALA 2 -0.520 2.679 -1.400 1.00 1.00 DIA H +ATOM 14 OY ALA 2 -1.139 0.931 -0.973 0.00 0.00 DIA O +ATOM 16 HN ALA 2 1.713 1.021 -0.873 1.00 1.00 DIA H +ATOM 18 HA ALA 2 0.099 -0.774 -2.218 0.00 0.00 DIA H +ATOM 19 CB ALA 2 2.063 -1.223 -1.276 1.00 1.00 DIA C +ATOM 20 HB1 ALA 2 2.670 -0.716 -2.057 1.00 1.00 DIA H +ATOM 21 HB2 ALA 2 2.556 -1.051 -0.295 1.00 1.00 DIA H +ATOM 22 HB3 ALA 2 2.070 -2.314 -1.490 1.00 1.00 DIA H diff --git a/regtest/basic/rt43/test4.pdb b/regtest/basic/rt43/test4.pdb new file mode 100644 index 0000000000000000000000000000000000000000..455f79f38cdf72be596b1883f5be162bf37d5983 --- /dev/null +++ b/regtest/basic/rt43/test4.pdb @@ -0,0 +1,15 @@ +ATOM 2 O ALA 2 -0.926 -2.447 -0.497 1.00 0.00 DIA O +ATOM 4 HNT ALA 2 0.533 -0.396 1.184 1.00 0.00 DIA H +ATOM 6 HT1 ALA 2 -0.216 -2.590 1.371 1.00 0.00 DIA H +ATOM 7 HT2 ALA 2 -0.309 -1.255 2.315 1.00 0.00 DIA H +ATOM 8 HT3 ALA 2 -1.480 -1.560 1.212 1.00 0.00 DIA H +ATOM 9 CAY ALA 2 -0.096 2.144 -0.669 1.00 0.00 DIA C +ATOM 10 HY1 ALA 2 0.871 2.385 -0.588 1.00 0.00 DIA H +ATOM 12 HY3 ALA 2 -0.520 2.679 -1.400 1.00 0.00 DIA H +ATOM 14 OY ALA 2 -1.139 0.931 -0.973 0.00 1.00 DIA O +ATOM 16 HN ALA 2 1.713 1.021 -0.873 0.00 1.00 DIA H +ATOM 18 HA ALA 2 0.099 -0.774 -2.218 0.00 1.00 DIA H +ATOM 19 CB ALA 2 2.063 -1.223 -1.276 0.00 1.00 DIA C +ATOM 20 HB1 ALA 2 2.670 -0.716 -2.057 0.00 1.00 DIA H +ATOM 21 HB2 ALA 2 2.556 -1.051 -0.295 0.00 1.00 DIA H +ATOM 22 HB3 ALA 2 2.070 -2.314 -1.490 0.00 1.00 DIA H diff --git a/src/function/Stats.cpp b/src/function/Stats.cpp new file mode 100644 index 0000000000000000000000000000000000000000..c3f9035bff321a7e9d22ae33eaa813f8ed4b2365 --- /dev/null +++ b/src/function/Stats.cpp @@ -0,0 +1,166 @@ +/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + Copyright (c) 2011-2015 The plumed team + (see the PEOPLE file at the root of the distribution for a list of names) + + See http://www.plumed-code.org for more information. + + This file is part of plumed, version 2. + + plumed is free software: you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + plumed is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public License + along with plumed. If not, see <http://www.gnu.org/licenses/>. ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ +#include "Function.h" +#include "ActionRegister.h" + +#include <cmath> + +using namespace std; + +namespace PLMD{ +namespace function{ + +//+PLUMEDOC FUNCTION STATS +/* + +\endverbatim + +*/ +//+ENDPLUMEDOC + + +class Stats : + public Function +{ + std::vector<double> parameters; + bool sqdonly; +public: + explicit Stats(const ActionOptions&); + void calculate(); + static void registerKeywords(Keywords& keys); +}; + + +PLUMED_REGISTER_ACTION(Stats,"STATS") + +void Stats::registerKeywords(Keywords& keys){ + Function::registerKeywords(keys); + keys.use("ARG"); + keys.add("compulsory","PARAMETERS","0.0","the parameters of the arguments in your function"); + keys.addFlag("SQDEVSUM",false,"calculates only SQDEVSUM"); +} + +Stats::Stats(const ActionOptions&ao): +Action(ao), +Function(ao), +parameters(getNumberOfArguments(),0.0), +sqdonly(false) +{ + parseVector("PARAMETERS",parameters); + if(parameters.size()!=static_cast<unsigned>(getNumberOfArguments())) + error("Size of PARAMETERS array should be the same as number for arguments"); + + parseFlag("SQDEVSUM",sqdonly); + + + addComponentWithDerivatives("sqdevsum"); + componentIsNotPeriodic("sqdevsum"); + if(!sqdonly) { + addComponentWithDerivatives("corr"); + componentIsNotPeriodic("corr"); + addComponentWithDerivatives("slope"); + componentIsNotPeriodic("slope"); + addComponentWithDerivatives("intercept"); + componentIsNotPeriodic("intercept"); + } + + checkRead(); +} + +void Stats::calculate() +{ + if(sqdonly) { + + double nsqd = 0.; + Value* valuea=getPntrToComponent("sqdevsum"); + for(unsigned i=0;i<parameters.size();++i){ + double dev = getArgument(i)-parameters[i]; + setDerivative(valuea,i,2.*dev); + nsqd += dev*dev; + } + valuea->set(nsqd); + + } else { + + double scx=0., scx2=0., scy=0., scy2=0., scxy=0.; + + for(unsigned i=0;i<parameters.size();++i){ + scx += getArgument(i); + scx2 += getArgument(i)*getArgument(i); + scy += parameters[i]; + scy2 += parameters[i]*parameters[i]; + scxy += getArgument(i)*parameters[i]; + } + + double ns = parameters.size(); + + double num = ns*scxy - scx*scy; + double idev2x = 1./(ns*scx2-scx*scx); + double idevx = sqrt(idev2x); + double idevy = 1./sqrt(ns*scy2-scy*scy); + + /* sd */ + double nsqd = scx2 + scy2 - 2.*scxy; + /* correlation */ + double correlation = num * idevx * idevy; + /* slope and intercept */ + double slope = num * idev2x; + double inter = (scy - slope * scx)/ns; + + Value* valuea=getPntrToComponent("sqdevsum"); + Value* valueb=getPntrToComponent("corr"); + Value* valuec=getPntrToComponent("slope"); + Value* valued=getPntrToComponent("intercept"); + + valuea->set(nsqd); + valueb->set(correlation); + valuec->set(slope); + valued->set(inter); + + /* derivatives */ + for(unsigned i=0;i<parameters.size();++i){ + double common_d1 = (ns*parameters[i]-scy)*idevx; + double common_d2 = num*(ns*getArgument(i)-scx)*idev2x*idevx; + double common_d3 = common_d1 - common_d2; + + /* sqdevsum */ + double sq_der = 2.*(getArgument(i)-parameters[i]); + /* correlation */ + double co_der = common_d3*idevy; + /* slope */ + double sl_der = (common_d1-2.*common_d2)*idevx; + /* intercept */ + double int_der = -(slope+ scx*sl_der)/ns; + + setDerivative(valuea,i,sq_der); + setDerivative(valueb,i,co_der); + setDerivative(valuec,i,sl_der); + setDerivative(valued,i,int_der); + } + + } +} + +} +} + +