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