Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-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 : #ifndef __PLUMED_colvar_MultiColvarTemplate_h
23 : #define __PLUMED_colvar_MultiColvarTemplate_h
24 :
25 : #include "core/ActionWithVector.h"
26 :
27 : namespace PLMD {
28 : namespace colvar {
29 :
30 : template <class T>
31 : class MultiColvarTemplate : public ActionWithVector {
32 : private:
33 : /// An index that decides what we are calculating
34 : unsigned mode;
35 : /// Are we using pbc to calculate the CVs
36 : bool usepbc;
37 : /// Do we reassemble the molecule
38 : bool wholemolecules;
39 : /// Variable that is turned true if charges are needed to compute the underlying CVs
40 : /// Add the command:
41 : /// keys.add("hidden","NEEDS_CHARGES","Inform the parallel task manager that we need the charges to compute this CV");
42 : /// in the registerKeywords method for your colvar to turn this to true
43 : bool needs_charges;
44 : /// Blocks of atom numbers
45 : std::vector< std::vector<unsigned> > ablocks;
46 : public:
47 : static void registerKeywords(Keywords&);
48 : explicit MultiColvarTemplate(const ActionOptions&);
49 : unsigned getNumberOfDerivatives() override ;
50 : void addValueWithDerivatives( const std::vector<unsigned>& shape=std::vector<unsigned>() ) override ;
51 : void addComponentWithDerivatives( const std::string& name, const std::vector<unsigned>& shape=std::vector<unsigned>() ) override ;
52 : void setupStreamedComponents( const std::string& headstr, unsigned& nquants, unsigned& nmat, unsigned& maxcol, unsigned& nbookeeping ) override ;
53 : void performTask( const unsigned&, MultiValue& ) const override ;
54 : void calculate() override;
55 : };
56 :
57 : template <class T>
58 454 : void MultiColvarTemplate<T>::registerKeywords(Keywords& keys ) {
59 454 : T::registerKeywords( keys );
60 454 : unsigned nkeys = keys.size();
61 4168 : for(unsigned i=0; i<nkeys; ++i) {
62 7428 : if( keys.style( keys.get(i), "atoms" ) ) {
63 1516 : keys.reset_style( keys.get(i), "numbered" );
64 : }
65 : }
66 908 : if( keys.outputComponentExists(".#!value") ) {
67 660 : keys.setValueDescription("the " + keys.getDisplayName() + " for each set of specified atoms");
68 : }
69 454 : }
70 :
71 : template <class T>
72 205 : MultiColvarTemplate<T>::MultiColvarTemplate(const ActionOptions&ao):
73 : Action(ao),
74 : ActionWithVector(ao),
75 205 : mode(0),
76 205 : usepbc(true),
77 205 : wholemolecules(false),
78 205 : needs_charges(keywords.exists("NEEDS_CHARGES")) {
79 : std::vector<AtomNumber> all_atoms;
80 205 : if( getName()=="POSITION_VECTOR" || getName()=="MASS_VECTOR" || getName()=="CHARGE_VECTOR" ) {
81 122 : parseAtomList( "ATOMS", all_atoms );
82 : }
83 205 : if( all_atoms.size()>0 ) {
84 59 : ablocks.resize(1);
85 59 : ablocks[0].resize( all_atoms.size() );
86 8137 : for(unsigned i=0; i<all_atoms.size(); ++i) {
87 8078 : ablocks[0][i] = i;
88 : }
89 : } else {
90 : std::vector<AtomNumber> t;
91 32331 : for(int i=1;; ++i ) {
92 32331 : T::parseAtomList( i, t, this );
93 32331 : if( t.empty() ) {
94 : break;
95 : }
96 :
97 32185 : if( i==1 ) {
98 146 : ablocks.resize(t.size());
99 : }
100 32185 : if( t.size()!=ablocks.size() ) {
101 : std::string ss;
102 0 : Tools::convert(i,ss);
103 0 : error("ATOMS" + ss + " keyword has the wrong number of atoms");
104 : }
105 98871 : for(unsigned j=0; j<ablocks.size(); ++j) {
106 66686 : ablocks[j].push_back( ablocks.size()*(i-1)+j );
107 66686 : all_atoms.push_back( t[j] );
108 : }
109 32185 : t.resize(0);
110 : }
111 : }
112 205 : if( all_atoms.size()==0 ) {
113 0 : error("No atoms have been specified");
114 : }
115 205 : requestAtoms(all_atoms);
116 410 : if( keywords.exists("NOPBC") ) {
117 205 : bool nopbc=!usepbc;
118 205 : parseFlag("NOPBC",nopbc);
119 205 : usepbc=!nopbc;
120 : }
121 410 : if( keywords.exists("WHOLEMOLECULES") ) {
122 43 : parseFlag("WHOLEMOLECULES",wholemolecules);
123 43 : if( wholemolecules ) {
124 1 : usepbc=false;
125 : }
126 : }
127 205 : if( usepbc ) {
128 168 : log.printf(" using periodic boundary conditions\n");
129 : } else {
130 37 : log.printf(" without periodic boundary conditions\n");
131 : }
132 :
133 : // Setup the values
134 205 : mode = T::getModeAndSetupValues( this );
135 205 : }
136 :
137 : template <class T>
138 1616 : unsigned MultiColvarTemplate<T>::getNumberOfDerivatives() {
139 1616 : return 3*getNumberOfAtoms()+9;
140 : }
141 :
142 : template <class T>
143 23814 : void MultiColvarTemplate<T>::calculate() {
144 23814 : runAllTasks();
145 23814 : }
146 :
147 : template <class T>
148 101 : void MultiColvarTemplate<T>::addValueWithDerivatives( const std::vector<unsigned>& shape ) {
149 101 : std::vector<unsigned> s(1);
150 101 : s[0]=ablocks[0].size();
151 101 : addValue( s );
152 101 : }
153 :
154 : template <class T>
155 324 : void MultiColvarTemplate<T>::addComponentWithDerivatives( const std::string& name, const std::vector<unsigned>& shape ) {
156 324 : std::vector<unsigned> s(1);
157 324 : s[0]=ablocks[0].size();
158 324 : addComponent( name, s );
159 324 : }
160 :
161 : template <class T>
162 29220 : void MultiColvarTemplate<T>::setupStreamedComponents( const std::string& headstr, unsigned& nquants, unsigned& nmat, unsigned& maxcol, unsigned& nbookeeping ) {
163 29220 : if( wholemolecules ) {
164 1 : makeWhole();
165 : }
166 29220 : ActionWithVector::setupStreamedComponents( headstr, nquants, nmat, maxcol, nbookeeping );
167 29220 : }
168 :
169 : template <class T>
170 504415 : void MultiColvarTemplate<T>::performTask( const unsigned& task_index, MultiValue& myvals ) const {
171 : // Retrieve the positions
172 : std::vector<Vector> & fpositions( myvals.getFirstAtomVector() );
173 504415 : if( fpositions.size()!=ablocks.size() ) {
174 40026 : fpositions.resize( ablocks.size() );
175 : }
176 1363634 : for(unsigned i=0; i<ablocks.size(); ++i) {
177 859219 : fpositions[i] = getPosition( ablocks[i][task_index] );
178 : }
179 : // If we are using pbc make whole
180 504415 : if( usepbc ) {
181 233362 : if( fpositions.size()==1 ) {
182 57816 : fpositions[0]=pbcDistance(Vector(0.0,0.0,0.0),getPosition( ablocks[0][task_index] ) );
183 : } else {
184 418659 : for(unsigned j=0; j<fpositions.size()-1; ++j) {
185 : const Vector & first (fpositions[j]);
186 214205 : Vector & second (fpositions[j+1]);
187 214205 : second=first+pbcDistance(first,second);
188 : }
189 : }
190 271053 : } else if( fpositions.size()==1 ) {
191 133776 : fpositions[0]=delta(Vector(0.0,0.0,0.0),getPosition( ablocks[0][task_index] ) );
192 : }
193 : // Retrieve the masses and charges
194 504415 : myvals.resizeTemporyVector(2);
195 : std::vector<double> & mass( myvals.getTemporyVector(0) );
196 : std::vector<double> & charge( myvals.getTemporyVector(1) );
197 504415 : if( mass.size()!=ablocks.size() ) {
198 38994 : mass.resize(ablocks.size());
199 38994 : charge.resize(ablocks.size());
200 : }
201 504415 : if( chargesWereSet || needs_charges ) {
202 1363589 : for(unsigned i=0; i<ablocks.size(); ++i) {
203 859189 : mass[i]=getMass( ablocks[i][task_index] );
204 859189 : charge[i]=getCharge( ablocks[i][task_index] );
205 : }
206 : } else {
207 45 : for(unsigned i=0; i<ablocks.size(); ++i) {
208 30 : mass[i]=getMass( ablocks[i][task_index] );
209 30 : charge[i]=0;
210 : }
211 : }
212 : // Make some space to store various things
213 504415 : std::vector<double> values( getNumberOfComponents() );
214 : std::vector<Tensor> & virial( myvals.getFirstAtomVirialVector() );
215 : std::vector<std::vector<Vector> > & derivs( myvals.getFirstAtomDerivativeVector() );
216 504415 : if( derivs.size()!=values.size() ) {
217 40236 : derivs.resize( values.size() );
218 40236 : virial.resize( values.size() );
219 : }
220 1535399 : for(unsigned i=0; i<derivs.size(); ++i) {
221 1030984 : if( derivs[i].size()<ablocks.size() ) {
222 74634 : derivs[i].resize( ablocks.size() );
223 : }
224 : }
225 : // Calculate the CVs using the method in the Colvar
226 504415 : T::calculateCV( mode, mass, charge, fpositions, values, derivs, virial, this );
227 1535399 : for(unsigned i=0; i<values.size(); ++i) {
228 1030984 : myvals.setValue( getConstPntrToComponent(i)->getPositionInStream(), values[i] );
229 : }
230 : // Finish if there are no derivatives
231 504415 : if( doNotCalculateDerivatives() ) {
232 : return;
233 : }
234 :
235 : // Now transfer the derivatives to the underlying MultiValue
236 869373 : for(unsigned i=0; i<ablocks.size(); ++i) {
237 548302 : unsigned base=3*ablocks[i][task_index];
238 1457289 : for(int j=0; j<getNumberOfComponents(); ++j) {
239 908987 : unsigned jval=getConstPntrToComponent(j)->getPositionInStream();
240 908987 : myvals.addDerivative( jval, base + 0, derivs[j][i][0] );
241 908987 : myvals.addDerivative( jval, base + 1, derivs[j][i][1] );
242 908987 : myvals.addDerivative( jval, base + 2, derivs[j][i][2] );
243 : }
244 : // Check for duplicated indices during update to avoid double counting
245 : bool newi=true;
246 788376 : for(unsigned j=0; j<i; ++j) {
247 240074 : if( ablocks[j][task_index]==ablocks[i][task_index] ) {
248 : newi=false;
249 : break;
250 : }
251 : }
252 548302 : if( !newi ) {
253 0 : continue;
254 : }
255 1457289 : for(int j=0; j<getNumberOfComponents(); ++j) {
256 908987 : unsigned jval=getConstPntrToComponent(j)->getPositionInStream();
257 908987 : myvals.updateIndex( jval, base );
258 908987 : myvals.updateIndex( jval, base + 1 );
259 908987 : myvals.updateIndex( jval, base + 2 );
260 : }
261 : }
262 321071 : unsigned nvir=3*getNumberOfAtoms();
263 913341 : for(int j=0; j<getNumberOfComponents(); ++j) {
264 592270 : unsigned jval=getConstPntrToComponent(j)->getPositionInStream();
265 2369080 : for(unsigned i=0; i<3; ++i) {
266 7107240 : for(unsigned k=0; k<3; ++k) {
267 5330430 : myvals.addDerivative( jval, nvir + 3*i + k, virial[j][i][k] );
268 5330430 : myvals.updateIndex( jval, nvir + 3*i + k );
269 : }
270 : }
271 : }
272 : }
273 :
274 : }
275 : }
276 : #endif
|
