Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2014-2018 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 "Colvar.h"
23 : #include "ActionRegister.h"
24 : #include "tools/NeighborList.h"
25 : #include "tools/OpenMP.h"
26 :
27 : #include <string>
28 : #include <cmath>
29 :
30 : using namespace std;
31 :
32 : namespace PLMD {
33 : namespace colvar {
34 :
35 : //+PLUMEDOC COLVAR NOE
36 : /*
37 : Calculates NOE intensities as sums of 1/r^6, also averaging over multiple equivalent atoms
38 : or ambiguous NOE.
39 :
40 : Each NOE is defined by two groups containing the same number of atoms, distances are
41 : calculated in pairs, transformed in 1/r^6, summed and saved as components.
42 :
43 : \f[
44 : NOE() = (\frac{1}{N_{eq}}\sum_j^{N_{eq}} (\frac{1}{r_j^6}))^{\frac{-1}{6}}
45 : \f]
46 :
47 : Intensities can then in principle ensemble averaged using \ref ENSEMBLE and used to
48 : calculate a scoring function for example with \ref METAINFERENCE.
49 :
50 : \par Examples
51 :
52 : In the following examples three noes are defined, the first is calculated based on the distances
53 : of atom 1-2 and 3-2; the second is defined by the distance 5-7 and the third by the distances
54 : 4-15,4-16,8-15,8-16.
55 :
56 : \verbatim
57 : NOE ...
58 : GROUPA1=1,3 GROUPB1=2,2
59 : GROUPA2=5 GROUPB2=7
60 : GROUPA3=4,4,8,8 GROUPB3=15,16,15,16
61 : LABEL=noes
62 : ... NOE
63 :
64 : PRINT ARG=noes.* FILE=colvar
65 : \endverbatim
66 : (See also \ref PRINT)
67 :
68 : */
69 : //+ENDPLUMEDOC
70 :
71 : class NOE : public Colvar {
72 : private:
73 : bool pbc;
74 : vector<unsigned> nga;
75 : NeighborList *nl;
76 : public:
77 : static void registerKeywords( Keywords& keys );
78 : explicit NOE(const ActionOptions&);
79 : ~NOE();
80 : virtual void calculate();
81 : };
82 :
83 2529 : PLUMED_REGISTER_ACTION(NOE,"NOE")
84 :
85 7 : void NOE::registerKeywords( Keywords& keys ) {
86 7 : Colvar::registerKeywords( keys );
87 7 : componentsAreNotOptional(keys);
88 7 : useCustomisableComponents(keys);
89 : keys.add("numbered","GROUPA","the atoms involved in each of the contacts you wish to calculate. "
90 : "Keywords like GROUPA1, GROUPA2, GROUPA3,... should be listed and one contact will be "
91 7 : "calculated for each ATOM keyword you specify.");
92 : keys.add("numbered","GROUPB","the atoms involved in each of the contacts you wish to calculate. "
93 : "Keywords like GROUPB1, GROUPB2, GROUPB3,... should be listed and one contact will be "
94 7 : "calculated for each ATOM keyword you specify.");
95 7 : keys.reset_style("GROUPA","atoms");
96 7 : keys.reset_style("GROUPB","atoms");
97 7 : keys.addFlag("ADDEXP",false,"Set to TRUE if you want to have fixed components with the experimental reference values.");
98 7 : keys.add("numbered","NOEDIST","Add an experimental value for each NOE.");
99 7 : keys.addOutputComponent("noe","default","the # NOE");
100 7 : keys.addOutputComponent("exp","ADDEXP","the # NOE experimental distance");
101 7 : }
102 :
103 6 : NOE::NOE(const ActionOptions&ao):
104 : PLUMED_COLVAR_INIT(ao),
105 6 : pbc(true)
106 : {
107 6 : bool nopbc=!pbc;
108 6 : parseFlag("NOPBC",nopbc);
109 6 : pbc=!nopbc;
110 :
111 : // Read in the atoms
112 12 : vector<AtomNumber> t, ga_lista, gb_lista;
113 18 : for(int i=1;; ++i ) {
114 18 : parseAtomList("GROUPA", i, t );
115 18 : if( t.empty() ) break;
116 12 : for(unsigned j=0; j<t.size(); j++) ga_lista.push_back(t[j]);
117 12 : nga.push_back(t.size());
118 12 : t.resize(0);
119 12 : }
120 12 : vector<unsigned> ngb;
121 18 : for(int i=1;; ++i ) {
122 18 : parseAtomList("GROUPB", i, t );
123 18 : if( t.empty() ) break;
124 12 : for(unsigned j=0; j<t.size(); j++) gb_lista.push_back(t[j]);
125 12 : ngb.push_back(t.size());
126 12 : if(ngb[i-1]!=nga[i-1]) error("The same number of atoms is expected for the same GROUPA-GROUPB couple");
127 12 : t.resize(0);
128 12 : }
129 6 : if(nga.size()!=ngb.size()) error("There should be the same number of GROUPA and GROUPB keywords");
130 : // Create neighbour lists
131 6 : nl= new NeighborList(ga_lista,gb_lista,true,pbc,getPbc());
132 :
133 6 : bool addexp=false;
134 6 : parseFlag("ADDEXP",addexp);
135 :
136 12 : vector<double> noedist;
137 6 : if(addexp) {
138 5 : noedist.resize( nga.size() );
139 5 : unsigned ntarget=0;
140 :
141 15 : for(unsigned i=0; i<nga.size(); ++i) {
142 10 : if( !parseNumbered( "NOEDIST", i+1, noedist[i] ) ) break;
143 10 : ntarget++;
144 : }
145 5 : if( ntarget!=nga.size() ) error("found wrong number of NOEDIST values");
146 : }
147 :
148 : // Ouput details of all contacts
149 6 : unsigned index=0;
150 18 : for(unsigned i=0; i<nga.size(); ++i) {
151 12 : log.printf(" The %uth NOE is calculated using %u equivalent couples of atoms\n", i, nga[i]);
152 30 : for(unsigned j=0; j<nga[i]; j++) {
153 18 : log.printf(" couple %u is %d %d.\n", j, ga_lista[index].serial(), gb_lista[index].serial() );
154 18 : index++;
155 : }
156 : }
157 :
158 6 : if(pbc) log.printf(" using periodic boundary conditions\n");
159 0 : else log.printf(" without periodic boundary conditions\n");
160 :
161 18 : for(unsigned i=0; i<nga.size(); i++) {
162 12 : string num; Tools::convert(i,num);
163 12 : addComponentWithDerivatives("noe_"+num);
164 12 : componentIsNotPeriodic("noe_"+num);
165 12 : }
166 :
167 6 : if(addexp) {
168 15 : for(unsigned i=0; i<nga.size(); i++) {
169 10 : string num; Tools::convert(i,num);
170 10 : addComponent("exp_"+num);
171 10 : componentIsNotPeriodic("exp_"+num);
172 10 : Value* comp=getPntrToComponent("exp_"+num);
173 10 : comp->set(noedist[i]);
174 10 : }
175 : }
176 :
177 6 : requestAtoms(nl->getFullAtomList());
178 12 : checkRead();
179 6 : }
180 :
181 24 : NOE::~NOE() {
182 6 : delete nl;
183 18 : }
184 :
185 396 : void NOE::calculate()
186 : {
187 396 : const unsigned ngasz=nga.size();
188 :
189 1171 : #pragma omp parallel for num_threads(OpenMP::getNumThreads())
190 775 : for(unsigned i=0; i<ngasz; i++) {
191 787 : Tensor dervir;
192 791 : double noe=0;
193 791 : unsigned index=0;
194 1187 : for(unsigned k=0; k<i; k++) index+=nga[k];
195 791 : const double c_aver=1./static_cast<double>(nga[i]);
196 791 : Value* val=getPntrToComponent(i);
197 : // cycle over equivalent atoms
198 1977 : for(unsigned j=0; j<nga[i]; j++) {
199 1188 : const unsigned i0=nl->getClosePair(index+j).first;
200 1181 : const unsigned i1=nl->getClosePair(index+j).second;
201 :
202 1187 : Vector distance;
203 2366 : if(pbc) distance=pbcDistance(getPosition(i0),getPosition(i1));
204 0 : else distance=delta(getPosition(i0),getPosition(i1));
205 :
206 1178 : const double ir2=1./distance.modulo2();
207 1173 : const double ir6=ir2*ir2*ir2;
208 1173 : const double ir8=6*ir6*ir2;
209 1173 : const double tmpir6=c_aver*ir6;
210 1173 : const double tmpir8=c_aver*ir8;
211 :
212 1173 : noe += tmpir6;
213 1173 : Vector deriv = tmpir8*distance;
214 1180 : dervir += Tensor(distance,deriv);
215 1187 : setAtomsDerivatives(val, i0, deriv);
216 1185 : setAtomsDerivatives(val, i1, -deriv);
217 : }
218 792 : val->set(noe);
219 792 : setBoxDerivatives(val, dervir);
220 : }
221 396 : }
222 :
223 : }
224 2523 : }
|
