Skip to content
Snippets Groups Projects
Commit 6c7d61bd authored by Carlo Camilloni's avatar Carlo Camilloni
Browse files

new function STATS 

it calculates the sum of the square dev, the correlation, the slope and the intercept of a linear fit
of a set of arguments with respect to a set of parameters
parent d199f294
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 300 additions and 0 deletions
regtest/basic/rt43/COLVAR.reference 0 → 100644
+ 6
− 0
View file @ 6c7d61bd
#! 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
regtest/basic/rt43/Makefile 0 → 100644
+ 1
− 0
View file @ 6c7d61bd
include ../../scripts/test.make
regtest/basic/rt43/config 0 → 100644
+ 4
− 0
View file @ 6c7d61bd
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"
regtest/basic/rt43/deriv.reference 0 → 100644
+ 26
− 0
View file @ 6c7d61bd
#! 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
regtest/basic/rt43/plumed.dat 0 → 100644
+ 22
− 0
View file @ 6c7d61bd
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
regtest/basic/rt43/test0.pdb 0 → 100644
+ 15
− 0
View file @ 6c7d61bd
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
regtest/basic/rt43/test1.pdb 0 → 100644
+ 15
− 0
View file @ 6c7d61bd
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
regtest/basic/rt43/test2.pdb 0 → 100644
+ 15
− 0
View file @ 6c7d61bd
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
regtest/basic/rt43/test3.pdb 0 → 100644
+ 15
− 0
View file @ 6c7d61bd
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
regtest/basic/rt43/test4.pdb 0 → 100644
+ 15
− 0
View file @ 6c7d61bd
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
src/function/Stats.cpp 0 → 100644
+ 166
− 0
View file @ 6c7d61bd
/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
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);
}
}
}
}
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment