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