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