Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-2019 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : plumed is distributed in the hope that it will be useful,
15 : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 : #include "Function.h"
23 : #include "ActionRegister.h"
24 : #include "tools/OpenMP.h"
25 :
26 : using namespace std;
27 :
28 : namespace PLMD {
29 : namespace function {
30 :
31 : //+PLUMEDOC FUNCTION LOCALENSEMBLE
32 : /*
33 : Calculates the average over multiple arguments.
34 :
35 : If more than one collective variable is given for each argument then they
36 : are averaged separately. The average is stored in a component labelled <em>label</em>.cvlabel.
37 :
38 : \par Examples
39 :
40 : The following input tells plumed to calculate the chemical shifts for four
41 : different proteins in the same simulation box then average them, calcualated
42 : the sum of the squared deviation with respect to the experiemntal values and
43 : applies a linear restraint.
44 : \plumedfile
45 : MOLINFO STRUCTURE=data/template.pdb
46 :
47 : chaina: GROUP ATOMS=1-1640
48 : chainb: GROUP ATOMS=1641-3280
49 : chainc: GROUP ATOMS=3281-4920
50 : chaind: GROUP ATOMS=4921-6560
51 :
52 : WHOLEMOLECULES ENTITY0=chaina ENTITY1=chainb ENTITY2=chainc ENTITY3=chaind
53 :
54 : csa: CS2BACKBONE ATOMS=chaina NRES=100 DATA=data/ TEMPLATE=chaina.pdb NOPBC
55 : csb: CS2BACKBONE ATOMS=chainb NRES=100 DATA=data/ TEMPLATE=chainb.pdb NOPBC
56 : csc: CS2BACKBONE ATOMS=chainc NRES=100 DATA=data/ TEMPLATE=chainc.pdb NOPBC
57 : csd: CS2BACKBONE ATOMS=chaind NRES=100 DATA=data/ TEMPLATE=chaind.pdb NOPBC
58 :
59 : ensca: LOCALENSEMBLE NUM=4 ARG1=(csa\.ca_.*) ARG2=(csb\.ca_.*) ARG3=(csc\.ca_.*) ARG4=(csd\.ca_.*)
60 : enscb: LOCALENSEMBLE NUM=4 ARG1=(csa\.cb_.*) ARG2=(csb\.cb_.*) ARG3=(csc\.cb_.*) ARG4=(csd\.cb_.*)
61 : ensco: LOCALENSEMBLE NUM=4 ARG1=(csa\.co_.*) ARG2=(csb\.co_.*) ARG3=(csc\.co_.*) ARG4=(csd\.co_.*)
62 : enshn: LOCALENSEMBLE NUM=4 ARG1=(csa\.hn_.*) ARG2=(csb\.hn_.*) ARG3=(csc\.hn_.*) ARG4=(csd\.hn_.*)
63 : ensnh: LOCALENSEMBLE NUM=4 ARG1=(csa\.nh_.*) ARG2=(csb\.nh_.*) ARG3=(csc\.nh_.*) ARG4=(csd\.nh_.*)
64 :
65 : stca: STATS ARG=(ensca\.csa\.ca_.*) PARARG=(csa\.expca_.*) SQDEVSUM
66 : stcb: STATS ARG=(enscb\.csa\.cb_.*) PARARG=(csa\.expcb_.*) SQDEVSUM
67 : stco: STATS ARG=(ensco\.csa\.co_.*) PARARG=(csa\.expco_.*) SQDEVSUM
68 : sthn: STATS ARG=(enshn\.csa\.hn_.*) PARARG=(csa\.exphn_.*) SQDEVSUM
69 : stnh: STATS ARG=(ensnh\.csa\.nh_.*) PARARG=(csa\.expnh_.*) SQDEVSUM
70 :
71 : res: RESTRAINT ARG=stca.*,stcb.*,stco.*,sthn.*,stnh.* AT=0.,0.,0.,0.,0. KAPPA=0.,0.,0.,0.,0 SLOPE=16.,16.,12.,24.,0.5
72 : \endplumedfile
73 :
74 : */
75 : //+ENDPLUMEDOC
76 :
77 :
78 4 : class LocalEnsemble :
79 : public Function
80 : {
81 : unsigned ens_dim;
82 : unsigned narg;
83 : public:
84 : explicit LocalEnsemble(const ActionOptions&);
85 : void calculate();
86 : static void registerKeywords(Keywords& keys);
87 : };
88 :
89 :
90 6454 : PLUMED_REGISTER_ACTION(LocalEnsemble,"LOCALENSEMBLE")
91 :
92 3 : void LocalEnsemble::registerKeywords(Keywords& keys) {
93 3 : Function::registerKeywords(keys);
94 6 : keys.use("ARG");
95 12 : keys.add("compulsory","NUM","the number of local replicas");
96 3 : ActionWithValue::useCustomisableComponents(keys);
97 3 : }
98 :
99 2 : LocalEnsemble::LocalEnsemble(const ActionOptions&ao):
100 : Action(ao),
101 : Function(ao),
102 2 : ens_dim(0)
103 : {
104 4 : parse("NUM",ens_dim);
105 2 : if(ens_dim==0) error("NUM should be greater or equal to 1");
106 :
107 : vector<Value*> arg;
108 : int oldsize=-1;
109 14 : for(unsigned i=1; i<=ens_dim; ++i ) {
110 : vector<Value*> larg;
111 12 : if(!parseArgumentList("ARG",i,larg)) break;
112 48 : for(unsigned j=0; j<larg.size(); j++) arg.push_back(larg[j]);
113 6 : if(oldsize!=-1&&oldsize!=larg.size()) error("In LOCALENSEMBLE you should have the same number of arguments for each ARG keyword");
114 6 : oldsize = larg.size();
115 6 : if(!larg.empty()) {
116 6 : log.printf(" with arguments %u: ", i);
117 60 : for(unsigned j=0; j<larg.size(); j++) log.printf(" %s",larg[j]->getName().c_str());
118 6 : log.printf("\n");
119 : }
120 : }
121 2 : requestArguments(arg);
122 2 : narg = arg.size()/ens_dim;
123 :
124 : // these are the averages
125 10 : for(unsigned i=0; i<narg; i++) {
126 : std::string s=getPntrToArgument(i)->getName();
127 4 : addComponentWithDerivatives(s);
128 4 : getPntrToComponent(i)->setNotPeriodic();
129 : }
130 :
131 2 : log.printf(" averaging over %u replicas.\n", ens_dim);
132 2 : }
133 :
134 40 : void LocalEnsemble::calculate()
135 : {
136 40 : const double fact = 1.0/static_cast<double>(ens_dim);
137 120 : #pragma omp parallel for num_threads(OpenMP::getNumThreads())
138 80 : for(unsigned i=0; i<narg; ++i) {
139 : double mean = 0.;
140 80 : Value* v=getPntrToComponent(i);
141 560 : for(unsigned j=0; j<ens_dim; ++j) {
142 240 : const unsigned index = j*narg+i;
143 : setDerivative(v, index, fact);
144 240 : mean += fact*getArgument(index);
145 : }
146 : v->set(mean);
147 : }
148 40 : }
149 :
150 : }
151 4839 : }
152 :
153 :
|
