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 52 : keys.addFlag("SERIAL",false,"Perform the calculation in serial - for debug purpose");
86 52 : keys.addFlag("NLIST",false,"Use a neighbour list to speed up the calculation");
87 52 : keys.add("optional","NL_CUTOFF","The cutoff for the neighbour list");
88 52 : keys.add("optional","NL_STRIDE","The frequency with which we are updating the atoms in the neighbour list");
89 52 : keys.add("atoms","METHYL_ATOM","the methyl carbon atom of the residue (i)");
90 52 : 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 52 : keys.add("atoms","HEAVY_ATOMS","the heavy atoms to be included in the calculation");
92 : //
93 52 : keys.add("compulsory","R_EFF","the effective distance, r_eff in the equation, given in nm.");
94 52 : keys.add("compulsory","PREFACTOR_A","the prefactor, a in the equation");
95 52 : keys.add("compulsory","EXPONENT_B","the exponent, b in the equation");
96 52 : keys.add("compulsory","OFFSET_C","the offset, c in the equation");
97 52 : keys.add("compulsory","N_I"," n_i in the equation");
98 52 : keys.add("optional","R_SHIFT","shift all distances by given amount");
99 26 : keys.setValueDescription("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 :
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) error("NL_CUTOFF should be explicitly specified and positive");
166 12 : parse("NL_STRIDE",nl_st);
167 12 : if(nl_st<=0) error("NL_STRIDE should be explicitly specified and positive");
168 : }
169 :
170 24 : parse("R_EFF",r_eff_);
171 24 : inv_r_eff_ = 1.0/r_eff_;
172 24 : parse("PREFACTOR_A",prefactor_a_);
173 24 : parse("EXPONENT_B",exp_b_);
174 24 : parse("OFFSET_C",offset_c_);
175 : unsigned int n_i_int;
176 24 : parse("N_I",n_i_int);
177 24 : n_i_ = static_cast<double>(n_i_int);
178 24 : total_prefactor_ = prefactor_a_/pow(n_i_,exp_b_);
179 : //
180 24 : parse("R_SHIFT",r_globalshift_);
181 :
182 24 : checkRead();
183 :
184 24 : addValueWithDerivatives();
185 24 : setNotPeriodic();
186 :
187 24 : bool dopair=false;
188 24 : if(doneigh) {
189 24 : nl=Tools::make_unique<NeighborList>(main_atoms,heavy_atoms,serial_,dopair,pbc_,getPbc(),comm,nl_cut,nl_st);
190 : }
191 : else {
192 24 : nl=Tools::make_unique<NeighborList>(main_atoms,heavy_atoms,serial_,dopair,pbc_,getPbc(),comm);
193 : }
194 :
195 24 : requestAtoms(nl->getFullAtomList());
196 :
197 :
198 24 : log.printf(" NMR S2 contact model CV, please read and cite ");
199 48 : log << plumed.cite("Palazzesi, Valsson, and Parrinello, J. Phys. Chem. Lett. 8, 4752 (2017) - DOI:10.1021/acs.jpclett.7b01770");
200 :
201 24 : if(modeltype_==methyl) {
202 0 : log << plumed.cite("Ming and Bruschweiler, J. Biomol. NMR, 29, 363 (2004) - DOI:10.1023/B:JNMR.0000032612.70767.35");
203 0 : log.printf("\n");
204 0 : log.printf(" calculation of methyl order parameter using atom %d \n",methyl_atom[0].serial());
205 : }
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 ) {log.printf(" \n");}
214 1848 : log.printf(" %d", heavy_atoms[i].serial());
215 : }
216 24 : log.printf(" \n");
217 24 : log.printf(" total number of distances: %u\n",nl->size());
218 : //
219 24 : log.printf(" using parameters");
220 24 : log.printf(" r_eff=%f,",r_eff_);
221 24 : log.printf(" a=%f,",prefactor_a_);
222 24 : log.printf(" b=%f,",exp_b_);
223 24 : log.printf(" c=%f,",offset_c_);
224 24 : log.printf(" n_i=%u",n_i_int);
225 24 : if(r_globalshift_!=0.0) {
226 4 : log.printf(", r_shift=%f",r_globalshift_);
227 : }
228 24 : log.printf("\n");
229 24 : if(pbc_) {
230 12 : log.printf(" using periodic boundary conditions\n");
231 : } else {
232 12 : log.printf(" without periodic boundary conditions\n");
233 : }
234 24 : if(doneigh) {
235 12 : log.printf(" using neighbor lists with\n");
236 12 : log.printf(" update every %d steps and cutoff %f\n",nl_st,nl_cut);
237 : }
238 24 : if(serial_) {
239 0 : log.printf(" calculation done in serial\n");
240 : }
241 24 : }
242 :
243 48 : S2ContactModel::~S2ContactModel() {
244 48 : }
245 :
246 48 : void S2ContactModel::prepare() {
247 48 : if(nl->getStride()>0) {
248 24 : if(firsttime || (getStep()%nl->getStride()==0)) {
249 24 : requestAtoms(nl->getFullAtomList());
250 24 : invalidateList=true;
251 24 : firsttime=false;
252 : } else {
253 0 : requestAtoms(nl->getReducedAtomList());
254 0 : invalidateList=false;
255 0 : if(getExchangeStep()) error("Neighbor lists should be updated on exchange steps - choose a NL_STRIDE which divides the exchange stride!");
256 : }
257 24 : if(getExchangeStep()) firsttime=true;
258 : }
259 48 : }
260 :
261 : // calculator
262 5952 : void S2ContactModel::calculate() {
263 :
264 5952 : Tensor virial;
265 5952 : std::vector<Vector> deriv(getNumberOfAtoms());
266 :
267 5952 : if(nl->getStride()>0 && invalidateList) {
268 2976 : nl->update(getPositions());
269 : }
270 :
271 5952 : unsigned int stride=comm.Get_size();
272 5952 : unsigned int rank=comm.Get_rank();
273 5952 : if(serial_) {
274 : stride=1;
275 : rank=0;
276 : }
277 :
278 5952 : double contact_sum = 0.0;
279 :
280 5952 : const unsigned int nn=nl->size();
281 :
282 321408 : for(unsigned int i=rank; i<nn; i+=stride) {
283 315456 : Vector distance;
284 315456 : unsigned int i0=nl->getClosePair(i).first;
285 315456 : unsigned int i1=nl->getClosePair(i).second;
286 315456 : if(getAbsoluteIndex(i0)==getAbsoluteIndex(i1)) {continue;}
287 :
288 315456 : if(pbc_) {
289 86304 : distance=pbcDistance(getPosition(i0),getPosition(i1));
290 : } else {
291 229152 : distance=delta(getPosition(i0),getPosition(i1));
292 : }
293 :
294 315456 : double exp_arg = exp(-(distance.modulo()-r_globalshift_)*inv_r_eff_);
295 315456 : contact_sum += exp_arg;
296 :
297 315456 : exp_arg /= distance.modulo();
298 315456 : Vector dd(exp_arg*distance);
299 315456 : Tensor vv(dd,distance);
300 315456 : deriv[i0]-=dd;
301 315456 : deriv[i1]+=dd;
302 315456 : virial-=vv;
303 : }
304 :
305 5952 : if(!serial_) {
306 5952 : comm.Sum(contact_sum);
307 5952 : if(!deriv.empty()) {
308 5952 : comm.Sum(&deriv[0][0],3*deriv.size());
309 : }
310 5952 : comm.Sum(virial);
311 : }
312 :
313 5952 : double value = tanh(total_prefactor_*contact_sum);
314 : // using that d/dx[tanh(x)]= 1-[tanh(x)]^2
315 5952 : double deriv_f = -inv_r_eff_*total_prefactor_*(1.0-value*value);
316 5952 : value -= offset_c_;
317 :
318 476160 : for(unsigned int i=0; i<deriv.size(); ++i) {
319 470208 : setAtomsDerivatives(i,deriv_f*deriv[i]);
320 : }
321 5952 : setValue(value);
322 11904 : setBoxDerivatives(deriv_f*virial);
323 :
324 5952 : }
325 : }
326 : }
|
