Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2017-2023 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : plumed is distributed in the hope that it will be useful,
15 : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 : #include "Colvar.h"
23 : #include "core/ActionRegister.h"
24 : #include "core/PlumedMain.h"
25 : #include "tools/Communicator.h"
26 :
27 : namespace PLMD {
28 : namespace colvar {
29 :
30 : //+PLUMEDOC COLVAR DIMER
31 : /*
32 : This CV computes the dimer interaction energy for a collection of dimers.
33 :
34 : Each dimer represents an atom, as described in the dimer paper \cite dimer-metad.
35 : A system of N atoms is thus represented with N dimers, each
36 : Dimer being composed of two beads and eventually a virtual site representing its center of mass.
37 :
38 : A typical configuration for a dimerized system has the following ordering of atoms:
39 :
40 : 1 TAG1 X Y Z N atoms representing the first bead of each Dimer
41 :
42 : 2 TAG2 X Y Z
43 :
44 : ...
45 :
46 : N TAGN X Y Z N atoms representing the second bead of each Dimer
47 :
48 : N+1 TAG1 X Y Z
49 :
50 : N+2 TAG2 X Y Z
51 :
52 : ...
53 :
54 : 2N TAGN X Y Z Optional: N atoms representing the center of mass of each Dimer
55 :
56 : 2N+1 TAG1 X Y Z
57 :
58 : 2N+2 TAG2 X Y Z
59 :
60 : ...
61 :
62 : 3N TAGN X Y Z The configuration might go on with un-dimerized atoms (like a solvent)
63 :
64 : 3N+1
65 :
66 : 3N+2
67 :
68 : ...
69 :
70 :
71 : The Dimer interaction energy is defined between atoms x and N+x, for x=1,...,N and is
72 : characterized by two parameters Q and DSIGMA. These are passed as mandatory arguments along with
73 : the temperature of the system.
74 :
75 : \par Examples
76 :
77 : This line tells Plumed to compute the Dimer interaction energy for every dimer in the system.
78 :
79 : \plumedfile
80 : dim: DIMER TEMP=300 Q=0.5 ALLATOMS DSIGMA=0.002
81 : \endplumedfile
82 :
83 : If the simulation doesn't use virtual sites for the dimers centers of mass,
84 : Plumed has to know in order to determine correctly the total number of dimers from
85 : the total number of atoms:
86 : \plumedfile
87 : dim: DIMER TEMP=300 Q=0.5 ALLATOMS DSIGMA=0.002 NOVSITES
88 : \endplumedfile
89 :
90 : The NOVSITES flag is not required if one provides the atom serials of each Dimer. These are
91 : defined through two lists of atoms provided __instead__ of the ALLATOMS keyword.
92 : For example, the Dimer interaction energy of dimers specified by beads (1;23),(5;27),(7;29) is:
93 : \plumedfile
94 : dim: DIMER TEMP=300 Q=0.5 ATOMS1=1,5,7 ATOMS2=23,27,29 DSIGMA=0.002
95 : \endplumedfile
96 :
97 : Note that the ATOMS1,ATOMS2 keywords can support atom groups and
98 : interval notation as defined in \ref GROUP.
99 :
100 :
101 : In a Replica Exchange simulation the keyword DSIGMA can be used in two ways:
102 : if a plumed.n.dat file is provided for each replica, then DSIGMA is passed as a single value,
103 : like in the previous examples, and each replica will read its own DSIGMA value. If
104 : a unique plumed.dat is given, DSIGMA has to be a list containing a value for each replica.
105 : For 4 replicas:
106 : \plumedfile
107 : #SETTINGS NREPLICAS=4
108 : dim: DIMER TEMP=300 Q=0.5 ATOMS1=1,5,7 ATOMS2=23,27,29 DSIGMA=0.002,0.002,0.004,0.01
109 : \endplumedfile
110 :
111 :
112 : \par Usage of the CV
113 :
114 : The dimer interaction is not coded in the driver program and has to be inserted
115 : in the Hamiltonian of the system as a linear RESTRAINT (see \ref RESTRAINT):
116 : \plumedfile
117 : dim: DIMER TEMP=300 Q=0.5 ALLATOMS DSIGMA=0.002
118 : RESTRAINT ARG=dim AT=0 KAPPA=0 SLOPE=1 LABEL=dimforces
119 : \endplumedfile
120 :
121 : In a replica exchange, Metadynamics (see \ref METAD) can be used on the Dimer CV to reduce
122 : the number of replicas. Just keep in mind that METAD SIGMA values should be tuned
123 : in the standard way for each replica according to the value of DSIGMA.
124 : */
125 : //+ENDPLUMEDOC
126 :
127 : class Dimer : public Colvar {
128 : public:
129 : static void registerKeywords( Keywords& keys);
130 : explicit Dimer(const ActionOptions&);
131 : void calculate() override;
132 : protected:
133 : bool trimer,useall;
134 : int myrank, nranks;
135 : double qexp,temperature,beta,dsigma;
136 : std::vector<double> dsigmas;
137 : private:
138 : void consistencyCheck();
139 : std::vector<AtomNumber> usedatoms1;
140 : std::vector<AtomNumber> usedatoms2;
141 :
142 : };
143 :
144 : PLUMED_REGISTER_ACTION(Dimer, "DIMER")
145 :
146 :
147 :
148 4 : void Dimer::registerKeywords( Keywords& keys) {
149 4 : Colvar::registerKeywords(keys);
150 :
151 8 : keys.add("compulsory","DSIGMA","The interaction strength of the dimer bond.");
152 8 : keys.add("compulsory", "Q", "The exponent of the dimer potential.");
153 8 : keys.add("compulsory", "TEMP", "The temperature (in Kelvin) of the simulation.");
154 8 : keys.add("atoms", "ATOMS1", "The list of atoms representing the first bead of each Dimer being considered by this CV. Used if ALLATOMS flag is missing");
155 8 : keys.add("atoms", "ATOMS2", "The list of atoms representing the second bead of each Dimer being considered by this CV. Used if ALLATOMS flag is missing");
156 8 : keys.addFlag("ALLATOMS", false, "Use EVERY atom of the system. Overrides ATOMS keyword.");
157 8 : keys.addFlag("NOVSITES", false, "If present the configuration is without virtual sites at the centroid positions.");
158 4 : keys.setValueDescription("the dimer interaction energy");
159 4 : }
160 :
161 :
162 :
163 2 : Dimer::Dimer(const ActionOptions& ao):
164 2 : PLUMED_COLVAR_INIT(ao) {
165 :
166 4 : log<<" Bibliography "<<plumed.cite("M Nava, F. Palazzesi, C. Perego and M. Parrinello, J. Chem. Theory Comput. 13, 425(2017)")<<"\n";
167 2 : parseVector("DSIGMA",dsigmas);
168 2 : parse("Q",qexp);
169 4 : parse("TEMP",temperature);
170 :
171 :
172 : std::vector<AtomNumber> atoms;
173 2 : parseFlag("ALLATOMS",useall);
174 2 : trimer=true;
175 : bool notrim;
176 2 : parseFlag("NOVSITES",notrim);
177 2 : trimer=!notrim;
178 :
179 2 : nranks=multi_sim_comm.Get_size();
180 2 : myrank=multi_sim_comm.Get_rank();
181 2 : if(dsigmas.size()==1) {
182 2 : dsigma=dsigmas[0];
183 : } else {
184 0 : dsigma=dsigmas[myrank];
185 : }
186 :
187 :
188 :
189 :
190 2 : if(useall) {
191 : // go with every atom in the system but not the virtuals...
192 : int natoms;
193 1 : if(trimer) {
194 1 : natoms= 2*getTotAtoms()/3;
195 : } else {
196 0 : natoms=getTotAtoms()/2;
197 : }
198 :
199 45 : for(unsigned int i=0; i<((unsigned int)natoms); i++) {
200 : AtomNumber ati;
201 : ati.setIndex(i);
202 44 : atoms.push_back(ati);
203 : }
204 : } else { // serials for the first beads of each dimer are given
205 1 : parseAtomList("ATOMS1",usedatoms1);
206 2 : parseAtomList("ATOMS2",usedatoms2);
207 :
208 : int isz1 = usedatoms1.size();
209 :
210 5 : for(unsigned int i=0; i<isz1; i++) {
211 : AtomNumber ati;
212 4 : ati.setIndex(usedatoms1[i].index());
213 4 : atoms.push_back(ati);
214 : }
215 :
216 : int isz2 = usedatoms2.size();
217 5 : for(unsigned int i=0; i<isz2; i++) {
218 : AtomNumber atip2;
219 4 : atip2.setIndex(usedatoms2[i].index());
220 4 : atoms.push_back(atip2);
221 : }
222 :
223 : }
224 2 : consistencyCheck();
225 2 : checkRead();
226 2 : beta = 1./(kBoltzmann*temperature);
227 :
228 2 : addValueWithDerivatives(); // allocate
229 2 : requestAtoms(atoms);
230 2 : setNotPeriodic();
231 2 : }
232 :
233 4 : void Dimer::calculate() {
234 4 : double cv_val=0;
235 4 : Tensor virial;
236 : std::vector<Vector> derivatives;
237 4 : std::vector<Vector> my_pos=getPositions();
238 4 : int atms = my_pos.size();
239 : std::vector<Vector> der_b2;
240 38 : for(int i=0; i<atms/2; i++) {
241 34 : Vector dist;
242 34 : dist = pbcDistance(my_pos[i],my_pos[i+atms/2]);
243 : double distquad=0;
244 136 : for(int j=0; j<3; j++) {
245 102 : distquad += dist[j]*dist[j];
246 : }
247 :
248 34 : double dsigquad = dsigma*dsigma;
249 34 : double fac1 = 1.0 + distquad/(2*qexp*dsigquad);
250 34 : double fac1qm1 = std::pow(fac1,qexp-1);
251 :
252 :
253 34 : cv_val += (fac1*fac1qm1-1.0)/beta;
254 34 : Vector der_val;
255 34 : Vector mder_val;
256 136 : for(int j=0; j<3; j++) {
257 102 : der_val[j] = -fac1qm1*dist[j]/(dsigquad*beta);
258 102 : mder_val[j]=-der_val[j];
259 : }
260 34 : derivatives.push_back(der_val);
261 34 : der_b2.push_back(mder_val);
262 :
263 : // virial part: each dimer contributes -x_{ij}*ds/dx_{ij} (s is the CV)
264 34 : double dfunc = fac1qm1/(beta*dsigquad);
265 34 : Vector dd(dfunc*dist);
266 34 : Tensor vv(dd,dist);
267 34 : virial -= vv;
268 :
269 : }
270 :
271 4 : derivatives.insert(derivatives.end(), der_b2.begin(), der_b2.end());
272 :
273 72 : for(unsigned int i=0; i<derivatives.size(); i++) {
274 68 : setAtomsDerivatives(i,derivatives[i]);
275 : }
276 :
277 4 : setValue(cv_val);
278 4 : setBoxDerivatives(virial);
279 :
280 4 : }
281 :
282 :
283 :
284 : /*****************
285 : There are some conditions that a valid input should satisfy.
286 : These are checked here and PLUMED error handlers are (eventually) called.
287 : ******************/
288 2 : void Dimer::consistencyCheck() {
289 2 : if(!useall && usedatoms1.size()!=usedatoms2.size()) {
290 0 : error("The provided atom lists are of different sizes.");
291 : }
292 :
293 2 : if(qexp<0.5 || qexp>1) {
294 0 : warning("Dimer CV is meant to be used with q-exponents between 0.5 and 1. We are not responsible for any black hole. :-)");
295 : }
296 :
297 2 : if(dsigma<0) {
298 0 : error("Please use positive sigma values for the Dimer strength constant");
299 : }
300 :
301 2 : if(temperature<0) {
302 0 : error("Please, use a positive value for the temperature...");
303 : }
304 :
305 : // if dsigmas has only one element means that either
306 : // you are using different plumed.x.dat files or a plumed.dat with a single replica
307 2 : if(dsigmas.size()!=nranks && dsigmas.size()!=1) {
308 0 : error("Mismatch between provided sigmas and number of replicas");
309 : }
310 :
311 2 : }
312 :
313 :
314 : }
315 : }
316 :
|
