Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2017-2022 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/Colvar.h"
23 : #include "core/ActionRegister.h"
24 : #include "core/PlumedMain.h"
25 : #include "tools/Communicator.h"
26 : #include "tools/Pbc.h"
27 : #include "tools/RMSD.h"
28 : #include <string>
29 :
30 : namespace PLMD {
31 : namespace isdb {
32 :
33 : //+PLUMEDOC ISDB_COLVAR SHADOW
34 : /*
35 : Communicate atoms positions among replicas and calculate the RMSD with respect to a mother (reference) simulation.
36 :
37 : The option UPDATE allows to specify the stride for communication between mother and replica systems.
38 : The flag REFERENCE needs to be specified in the input file of the mother replica.
39 : This action must be run in a multi-replica framework (such as the -multi option in GROMACS).
40 :
41 : ## Examples
42 :
43 : In this example, we perform a simulation of a RNA molecule using two replicas: a mother and a shadow replica.
44 : The mother simulation communicates the coordinates of the RNA backbone to the replica every 100 steps.
45 : The RMSD of the replica with respect to the mother is calculated on the RNA backbone atoms and an [UPPER_WALLS](UPPER_WALLS.md) is applied at 0.2 nm.
46 : The mother replica contains also the [UPPER_WALLS](UPPER_WALLS.md) action. However, the forces on the RNA atoms of the mother replica are automatically set to zero
47 : inside the [SHADOW](SHADOW.md) action.
48 :
49 : The input file for the mother simulation looks as follows:
50 :
51 : ```plumed
52 : # Reference PDB
53 : MOLINFO STRUCTURE=conf_emin_PLUMED.pdb WHOLE
54 : # Define RNA nucleic backbone
55 : rna: GROUP ATOMS=1,2,5,6,33,36,37,40,41,67,70,71,74,75,98,101,102,105,106,131,134,135,138,139,165,168,169,172,173,198,201,202,205,206,228,231,232,235,236,259,262,263,266,267,289,292,293,296,297,323,326,327,330,331,356,359,360,363,364,390,393,394,397,398,421,424,425,428,429,452,455,456,459,460,482,485,486,489,490,516,519,520,523,524,550,553,554,557,558,584,587,588,591,592,617,620,621,624,625,651,654,655,658,659,682,685,686,689,690,712,715,716,719,720,743,746,747,750,751,773,776,777,780,781,804,807,808,811,812,834,837,838,841,842,868,871,872,875,876,899,902,903,906,907
56 : # Reconstruct RNA PBC
57 : WHOLEMOLECULES ENTITY0=rna EMST STRIDE=1
58 :
59 : # Define shadow RMSD on RNA backbone
60 : rmsd: SHADOW ATOMS=rna NOPBC UPDATE=100 REFERENCE
61 : # Add upper wall - derivatives are set to zero inside SHADOW action
62 : uws: UPPER_WALLS ARG=rmsd AT=0.2 KAPPA=10000.0 STRIDE=1
63 :
64 : # Print useful info
65 : PRINT FILE=COLVAR STRIDE=500 ARG=rmsd,uws.bias
66 : ```
67 :
68 : while the input file for a shadow replica looks like:
69 :
70 : ```plumed
71 : # Reference PDB
72 : MOLINFO STRUCTURE=conf_emin_PLUMED.pdb WHOLE
73 : # Define RNA nucleic backbone
74 : rna: GROUP ATOMS=1,2,5,6,33,36,37,40,41,67,70,71,74,75,98,101,102,105,106,131,134,135,138,139,165,168,169,172,173,198,201,202,205,206,228,231,232,235,236,259,262,263,266,267,289,292,293,296,297,323,326,327,330,331,356,359,360,363,364,390,393,394,397,398,421,424,425,428,429,452,455,456,459,460,482,485,486,489,490,516,519,520,523,524,550,553,554,557,558,584,587,588,591,592,617,620,621,624,625,651,654,655,658,659,682,685,686,689,690,712,715,716,719,720,743,746,747,750,751,773,776,777,780,781,804,807,808,811,812,834,837,838,841,842,868,871,872,875,876,899,902,903,906,907
75 : # Reconstruct RNA PBC
76 : WHOLEMOLECULES ENTITY0=rna EMST STRIDE=1
77 :
78 : # Define shadow RMSD on RNA backbone
79 : rmsd: SHADOW ATOMS=rna NOPBC UPDATE=100
80 : # Add upper wall
81 : uws: UPPER_WALLS ARG=rmsd AT=0.2 KAPPA=10000.0 STRIDE=1
82 :
83 : # Print useful info
84 : PRINT FILE=COLVAR STRIDE=500 ARG=rmsd,uws.bias
85 : ```
86 :
87 : */
88 : //+ENDPLUMEDOC
89 :
90 : class Shadow : public Colvar {
91 : // private stuff
92 : bool isreference_;
93 : unsigned nupdate_;
94 : bool pbc_;
95 : bool first_time_;
96 : // RMSD object
97 : PLMD::RMSD rmsd_;
98 : // parallel stuff
99 : unsigned size_;
100 : unsigned rank_;
101 : // update reference
102 : void update_reference();
103 :
104 : public:
105 : static void registerKeywords( Keywords& keys );
106 : explicit Shadow(const ActionOptions&);
107 : // active methods:
108 : void calculate() override;
109 : };
110 :
111 : PLUMED_REGISTER_ACTION(Shadow,"SHADOW")
112 :
113 6 : void Shadow::registerKeywords( Keywords& keys ) {
114 6 : Colvar::registerKeywords( keys );
115 6 : keys.add("atoms","ATOMS","atoms for which we calculate the shadow RMSD");
116 6 : keys.add("compulsory","UPDATE","stride for updating reference coordinates");
117 6 : keys.addFlag("REFERENCE",false,"this is the reference replica");
118 12 : keys.setValueDescription("scalar","the value of the shadow RMSD");
119 6 : }
120 :
121 2 : Shadow::Shadow(const ActionOptions&ao):
122 : PLUMED_COLVAR_INIT(ao),
123 2 : isreference_(false), nupdate_(1), pbc_(true), first_time_(true) {
124 : // list of atoms
125 2 : std::vector<AtomNumber> atoms;
126 2 : parseAtomList("ATOMS",atoms);
127 : // update stride
128 2 : parse("UPDATE",nupdate_);
129 2 : if(nupdate_==0) {
130 0 : error("UPDATE must be greater than zero");
131 : }
132 :
133 : // is this the reference replica
134 2 : parseFlag("REFERENCE",isreference_);
135 :
136 : // periodic boundary conditions
137 2 : bool nopbc=!pbc_;
138 2 : parseFlag("NOPBC",nopbc);
139 2 : pbc_=!nopbc;
140 :
141 : // set intra-replica (openmp) parallel stuff
142 2 : size_ = comm.Get_size();
143 2 : rank_ = comm.Get_rank();
144 :
145 : // get number of (MPI) replicas
146 2 : int nrep = getNumberOfReplicas();
147 2 : int replica = getReplicaIndex();
148 :
149 : // exactly one replica must provide the reference coordinates
150 2 : unsigned nref = 0;
151 2 : if(rank_==0) {
152 2 : nref = isreference_ ? 1 : 0;
153 2 : multi_sim_comm.Sum(nref);
154 : }
155 2 : comm.Sum(nref);
156 2 : if(nref!=1) {
157 0 : error("exactly one replica must use the REFERENCE flag");
158 : }
159 : // check number of replicas
160 : //if(nrep<2) error("SHADOW must be used with at least two replicas");
161 :
162 2 : checkRead();
163 :
164 2 : log.printf(" atoms involved : ");
165 48 : for(unsigned i=0; i<atoms.size(); ++i) {
166 46 : log.printf("%d ",atoms[i].serial());
167 : }
168 2 : log.printf("\n");
169 2 : log.printf(" stride for updating reference coordinates : %d\n", nupdate_);
170 2 : log.printf(" number of replicas : %d\n", nrep);
171 2 : log.printf(" replica id : %d\n", replica);
172 2 : if(isreference_) {
173 1 : log.printf(" this is the reference replica\n");
174 : }
175 :
176 : // add value and set periodicity
177 2 : addValueWithDerivatives();
178 2 : setNotPeriodic();
179 :
180 : // request atoms
181 2 : requestAtoms(atoms);
182 2 : }
183 :
184 6 : void Shadow::update_reference() {
185 : // number of atoms
186 : unsigned natoms = getNumberOfAtoms();
187 : // initialize rmsd variables
188 6 : std::vector<double> align(natoms,1.0);
189 6 : std::vector<double> displace(natoms,1.0);
190 6 : std::vector<Vector> reference(natoms);
191 :
192 : // first get the reference coordinates
193 : // if master openmp task
194 6 : if(rank_==0) {
195 : // if reference replica
196 6 : if(isreference_) {
197 3 : reference = getPositions();
198 : }
199 : // share coordinates
200 6 : multi_sim_comm.Sum(&reference[0][0], 3*natoms);
201 : }
202 : // now intra replica (openmp) communication
203 6 : if(size_>1) {
204 0 : comm.Sum(&reference[0][0], 3*natoms);
205 : }
206 :
207 : // clear the rmsd object
208 6 : rmsd_.clear();
209 : // and initialize it
210 6 : rmsd_.set(align,displace,reference,"OPTIMAL");
211 12 : }
212 :
213 12 : void Shadow::calculate() {
214 : // make whole
215 12 : if(pbc_) {
216 0 : makeWhole();
217 : }
218 :
219 : // if it is time, update reference coordinates
220 12 : if(first_time_ || getStep()%nupdate_==0) {
221 6 : update_reference();
222 6 : first_time_ = false;
223 : }
224 :
225 : // calculate RMSD and derivatives
226 12 : std::vector<Vector> derivatives(getNumberOfAtoms());
227 12 : double rmsd = rmsd_.calculate(getPositions(), derivatives);
228 :
229 : // set RMSD value
230 12 : setValue(rmsd);
231 : // if this is not the reference replica, add derivatives
232 12 : if(!isreference_) {
233 144 : for(unsigned i=0; i<getNumberOfAtoms(); ++i) {
234 138 : setAtomsDerivatives(i, derivatives[i]);
235 : }
236 : }
237 : // set virial
238 12 : setBoxDerivativesNoPbc();
239 12 : }
240 :
241 : }
242 : }
|
