Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2017-2021 Omar Valsson
3 :
4 : This file is part of S2 contact model module
5 :
6 : The S2 contact model module is free software: you can redistribute it and/or modify
7 : it under the terms of the GNU Lesser General Public License as published by
8 : the Free Software Foundation, either version 3 of the License, or
9 : (at your option) any later version.
10 :
11 : The S2 contact model module is distributed in the hope that it will be useful,
12 : but WITHOUT ANY WARRANTY; without even the implied warranty of
13 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 : GNU Lesser General Public License for more details.
15 :
16 : You should have received a copy of the GNU Lesser General Public License
17 : along with the S2 contact model module. If not, see <http://www.gnu.org/licenses/>.
18 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
19 :
20 : #include "colvar/Colvar.h"
21 : #include "tools/NeighborList.h"
22 : #include "tools/Communicator.h"
23 : #include "core/ActionRegister.h"
24 : #include "core/PlumedMain.h"
25 :
26 :
27 : #include <string>
28 : #include <cmath>
29 :
30 : namespace PLMD {
31 : namespace s2cm {
32 :
33 :
34 : //
35 :
36 : //+PLUMEDOC COLVAR S2CM
37 : /*
38 : S2 contact model CV.
39 :
40 : This CV was used in the first paper cited below. It is based on NH order parameter from the second paper cited below and the methyl order parameter from the third. Input parameters can be found in the relevant papers.
41 :
42 :
43 : */
44 : //+ENDPLUMEDOC
45 :
46 :
47 :
48 :
49 :
50 :
51 :
52 : class S2ContactModel : public Colvar {
53 : bool pbc_;
54 : bool serial_;
55 : std::unique_ptr<NeighborList> nl;
56 : bool invalidateList;
57 : bool firsttime;
58 : //
59 : double r_eff_;
60 : double inv_r_eff_;
61 : double prefactor_a_;
62 : double exp_b_;
63 : double offset_c_;
64 : double n_i_;
65 : double total_prefactor_;
66 : double r_globalshift_;
67 :
68 : enum ModelType {methyl,nh} modeltype_;
69 :
70 : //
71 : public:
72 : explicit S2ContactModel(const ActionOptions&);
73 : ~S2ContactModel();
74 : virtual void calculate();
75 : virtual void prepare();
76 : static void registerKeywords( Keywords& keys );
77 : };
78 :
79 : PLUMED_REGISTER_ACTION(S2ContactModel,"S2CM")
80 :
81 26 : void S2ContactModel::registerKeywords( Keywords& keys ) {
82 26 : Colvar::registerKeywords(keys);
83 26 : keys.addFlag("SERIAL",false,"Perform the calculation in serial - for debug purpose");
84 26 : keys.addFlag("NLIST",false,"Use a neighbour list to speed up the calculation");
85 26 : keys.add("optional","NL_CUTOFF","The cutoff for the neighbour list");
86 26 : keys.add("optional","NL_STRIDE","The frequency with which we are updating the atoms in the neighbour list");
87 26 : keys.add("atoms","METHYL_ATOM","the methyl carbon atom of the residue (i)");
88 26 : keys.add("atoms","NH_ATOMS","the hydrogen atom of the NH group of the residue (i) and carbonyl oxygen of the preceding residue (i-1)");
89 26 : keys.add("atoms","HEAVY_ATOMS","the heavy atoms to be included in the calculation");
90 : //
91 26 : keys.add("compulsory","R_EFF","the effective distance, r_eff in the equation, given in nm.");
92 26 : keys.add("compulsory","PREFACTOR_A","the prefactor, a in the equation");
93 26 : keys.add("compulsory","EXPONENT_B","the exponent, b in the equation");
94 26 : keys.add("compulsory","OFFSET_C","the offset, c in the equation");
95 26 : keys.add("compulsory","N_I"," n_i in the equation");
96 26 : keys.add("optional","R_SHIFT","shift all distances by given amount");
97 52 : keys.setValueDescription("scalar","the value of the CV");
98 26 : keys.addDOI("10.1021/acs.jpclett.7b01770");
99 26 : keys.addDOI("10.1021/ja027847a");
100 26 : keys.addDOI("10.1023/B:JNMR.0000032612.70767.35");
101 26 : }
102 :
103 24 : S2ContactModel::S2ContactModel(const ActionOptions&ao):
104 : PLUMED_COLVAR_INIT(ao),
105 24 : pbc_(true),
106 24 : serial_(false),
107 24 : invalidateList(true),
108 24 : firsttime(true),
109 24 : r_eff_(0.0),
110 24 : inv_r_eff_(0.0),
111 24 : prefactor_a_(0.0),
112 24 : exp_b_(0.0),
113 24 : offset_c_(0.0),
114 24 : n_i_(0.0),
115 24 : total_prefactor_(0.0),
116 24 : r_globalshift_(0.0),
117 24 : modeltype_(methyl) {
118 :
119 48 : parseFlag("SERIAL",serial_);
120 :
121 : std::vector<AtomNumber> methyl_atom;
122 48 : parseAtomList("METHYL_ATOM",methyl_atom);
123 : std::vector<AtomNumber> nh_atoms;
124 48 : parseAtomList("NH_ATOMS",nh_atoms);
125 :
126 24 : if(methyl_atom.size()==0 && nh_atoms.size()==0) {
127 0 : error("you have to give either METHYL_ATOM or NH_ATOMS");
128 : }
129 24 : if(methyl_atom.size()>0 && nh_atoms.size()>0) {
130 0 : error("you cannot give both METHYL_ATOM or NH_ATOMS");
131 : }
132 24 : if(methyl_atom.size()>0 && methyl_atom.size()!=1) {
133 0 : error("you should give one atom in METHYL_ATOM, the methyl carbon atom of the residue");
134 : }
135 24 : if(nh_atoms.size()>0 && nh_atoms.size()!=2) {
136 0 : error("you should give two atoms in NH_ATOMS, the hydrogen atom of the NH group of the residue (i) and carbonyl oxygen of the preceding residue (i-1)");
137 : }
138 :
139 : std::vector<AtomNumber> heavy_atoms;
140 48 : parseAtomList("HEAVY_ATOMS",heavy_atoms);
141 24 : if(heavy_atoms.size()==0) {
142 0 : error("HEAVY_ATOMS is not given");
143 : }
144 :
145 : std::vector<AtomNumber> main_atoms;
146 24 : if(methyl_atom.size()>0) {
147 0 : modeltype_= methyl;
148 0 : main_atoms = methyl_atom;
149 24 : } else if(nh_atoms.size()>0) {
150 24 : modeltype_= nh;
151 24 : main_atoms = nh_atoms;
152 : }
153 :
154 24 : bool nopbc=!pbc_;
155 24 : parseFlag("NOPBC",nopbc);
156 24 : pbc_=!nopbc;
157 :
158 : // neighbor list stuff
159 24 : bool doneigh=false;
160 24 : double nl_cut=0.0;
161 24 : int nl_st=0;
162 24 : parseFlag("NLIST",doneigh);
163 24 : if(doneigh) {
164 12 : parse("NL_CUTOFF",nl_cut);
165 12 : if(nl_cut<=0.0) {
166 0 : error("NL_CUTOFF should be explicitly specified and positive");
167 : }
168 12 : parse("NL_STRIDE",nl_st);
169 12 : if(nl_st<=0) {
170 0 : error("NL_STRIDE should be explicitly specified and positive");
171 : }
172 : }
173 :
174 24 : parse("R_EFF",r_eff_);
175 24 : inv_r_eff_ = 1.0/r_eff_;
176 24 : parse("PREFACTOR_A",prefactor_a_);
177 24 : parse("EXPONENT_B",exp_b_);
178 24 : parse("OFFSET_C",offset_c_);
179 : unsigned int n_i_int;
180 24 : parse("N_I",n_i_int);
181 24 : n_i_ = static_cast<double>(n_i_int);
182 24 : total_prefactor_ = prefactor_a_/pow(n_i_,exp_b_);
183 : //
184 24 : parse("R_SHIFT",r_globalshift_);
185 :
186 24 : checkRead();
187 :
188 24 : addValueWithDerivatives();
189 24 : setNotPeriodic();
190 :
191 24 : bool dopair=false;
192 24 : if(doneigh) {
193 24 : nl=Tools::make_unique<NeighborList>(main_atoms,heavy_atoms,serial_,dopair,pbc_,getPbc(),comm,nl_cut,nl_st);
194 : } else {
195 24 : nl=Tools::make_unique<NeighborList>(main_atoms,heavy_atoms,serial_,dopair,pbc_,getPbc(),comm);
196 : }
197 :
198 24 : requestAtoms(nl->getFullAtomList());
199 :
200 :
201 24 : log.printf(" NMR S2 contact model CV, please read and cite ");
202 48 : log << plumed.cite("Palazzesi, Valsson, and Parrinello, J. Phys. Chem. Lett. 8, 4752 (2017) - DOI:10.1021/acs.jpclett.7b01770");
203 :
204 24 : if(modeltype_==methyl) {
205 0 : log << plumed.cite("Ming and Bruschweiler, J. Biomol. NMR, 29, 363 (2004) - DOI:10.1023/B:JNMR.0000032612.70767.35");
206 0 : log.printf("\n");
207 0 : log.printf(" calculation of methyl order parameter using atom %d \n",methyl_atom[0].serial());
208 24 : } else if(modeltype_==nh) {
209 48 : log << plumed.cite("Zhang and Bruschweiler, J. Am. Chem. Soc. 124, 12654 (2002) - DOI:10.1021/ja027847a");
210 24 : log.printf("\n");
211 24 : log.printf(" calculation of NH order parameter using atoms %d and %d\n",nh_atoms[0].serial(),nh_atoms[1].serial());
212 : }
213 24 : log.printf(" heavy atoms used in the calculation (%u in total):\n",static_cast<unsigned int>(heavy_atoms.size()));
214 1872 : for(unsigned int i=0; i<heavy_atoms.size(); ++i) {
215 1848 : if( (i+1) % 25 == 0 ) {
216 72 : log.printf(" \n");
217 : }
218 1848 : log.printf(" %d", heavy_atoms[i].serial());
219 : }
220 24 : log.printf(" \n");
221 24 : log.printf(" total number of distances: %u\n",nl->size());
222 : //
223 24 : log.printf(" using parameters");
224 24 : log.printf(" r_eff=%f,",r_eff_);
225 24 : log.printf(" a=%f,",prefactor_a_);
226 24 : log.printf(" b=%f,",exp_b_);
227 24 : log.printf(" c=%f,",offset_c_);
228 24 : log.printf(" n_i=%u",n_i_int);
229 24 : if(r_globalshift_!=0.0) {
230 4 : log.printf(", r_shift=%f",r_globalshift_);
231 : }
232 24 : log.printf("\n");
233 24 : if(pbc_) {
234 12 : log.printf(" using periodic boundary conditions\n");
235 : } else {
236 12 : log.printf(" without periodic boundary conditions\n");
237 : }
238 24 : if(doneigh) {
239 12 : log.printf(" using neighbor lists with\n");
240 12 : log.printf(" update every %d steps and cutoff %f\n",nl_st,nl_cut);
241 : }
242 24 : if(serial_) {
243 0 : log.printf(" calculation done in serial\n");
244 : }
245 24 : }
246 :
247 48 : S2ContactModel::~S2ContactModel() {
248 48 : }
249 :
250 48 : void S2ContactModel::prepare() {
251 48 : if(nl->getStride()>0) {
252 24 : if(firsttime || (getStep()%nl->getStride()==0)) {
253 24 : requestAtoms(nl->getFullAtomList());
254 24 : invalidateList=true;
255 24 : firsttime=false;
256 : } else {
257 0 : requestAtoms(nl->getReducedAtomList());
258 0 : invalidateList=false;
259 0 : if(getExchangeStep()) {
260 0 : error("Neighbor lists should be updated on exchange steps - choose a NL_STRIDE which divides the exchange stride!");
261 : }
262 : }
263 24 : if(getExchangeStep()) {
264 0 : firsttime=true;
265 : }
266 : }
267 48 : }
268 :
269 : // calculator
270 5952 : void S2ContactModel::calculate() {
271 :
272 : Tensor virial;
273 5952 : std::vector<Vector> deriv(getNumberOfAtoms());
274 :
275 5952 : if(nl->getStride()>0 && invalidateList) {
276 2976 : nl->update(getPositions());
277 : }
278 :
279 5952 : unsigned int stride=comm.Get_size();
280 5952 : unsigned int rank=comm.Get_rank();
281 5952 : if(serial_) {
282 : stride=1;
283 : rank=0;
284 : }
285 :
286 5952 : double contact_sum = 0.0;
287 :
288 5952 : const unsigned int nn=nl->size();
289 :
290 321408 : for(unsigned int i=rank; i<nn; i+=stride) {
291 : Vector distance;
292 315456 : unsigned int i0=nl->getClosePair(i).first;
293 315456 : unsigned int i1=nl->getClosePair(i).second;
294 315456 : if(getAbsoluteIndex(i0)==getAbsoluteIndex(i1)) {
295 0 : continue;
296 : }
297 :
298 315456 : if(pbc_) {
299 86304 : distance=pbcDistance(getPosition(i0),getPosition(i1));
300 : } else {
301 229152 : distance=delta(getPosition(i0),getPosition(i1));
302 : }
303 :
304 315456 : double exp_arg = exp(-(distance.modulo()-r_globalshift_)*inv_r_eff_);
305 315456 : contact_sum += exp_arg;
306 :
307 315456 : exp_arg /= distance.modulo();
308 : Vector dd(exp_arg*distance);
309 315456 : Tensor vv(dd,distance);
310 315456 : deriv[i0]-=dd;
311 315456 : deriv[i1]+=dd;
312 : virial-=vv;
313 : }
314 :
315 5952 : if(!serial_) {
316 5952 : comm.Sum(contact_sum);
317 5952 : if(!deriv.empty()) {
318 5952 : comm.Sum(&deriv[0][0],3*deriv.size());
319 : }
320 5952 : comm.Sum(virial);
321 : }
322 :
323 5952 : double value = tanh(total_prefactor_*contact_sum);
324 : // using that d/dx[tanh(x)]= 1-[tanh(x)]^2
325 5952 : double deriv_f = -inv_r_eff_*total_prefactor_*(1.0-value*value);
326 5952 : value -= offset_c_;
327 :
328 476160 : for(unsigned int i=0; i<deriv.size(); ++i) {
329 470208 : setAtomsDerivatives(i,deriv_f*deriv[i]);
330 : }
331 5952 : setValue(value);
332 11904 : setBoxDerivatives(deriv_f*virial);
333 :
334 5952 : }
335 : }
336 : }
|
