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 : /// Blocks of atom numbers
40 : std::vector< std::vector<unsigned> > ablocks;
41 : public:
42 : static void registerKeywords(Keywords&);
43 : explicit MultiColvarTemplate(const ActionOptions&);
44 : unsigned getNumberOfDerivatives() override ;
45 : void addValueWithDerivatives( const std::vector<unsigned>& shape=std::vector<unsigned>() ) override ;
46 : void addComponentWithDerivatives( const std::string& name, const std::vector<unsigned>& shape=std::vector<unsigned>() ) override ;
47 : void setupStreamedComponents( const std::string& headstr, unsigned& nquants, unsigned& nmat, unsigned& maxcol, unsigned& nbookeeping ) override ;
48 : void performTask( const unsigned&, MultiValue& ) const override ;
49 : void calculate() override;
50 : };
51 :
52 : template <class T>
53 426 : void MultiColvarTemplate<T>::registerKeywords(Keywords& keys ) {
54 426 : T::registerKeywords( keys );
55 426 : unsigned nkeys = keys.size();
56 3906 : for(unsigned i=0; i<nkeys; ++i) {
57 6960 : if( keys.style( keys.getKeyword(i), "atoms" ) ) {
58 1400 : keys.reset_style( keys.getKeyword(i), "numbered" );
59 : }
60 : }
61 852 : if( keys.outputComponentExists(".#!value") ) {
62 624 : keys.setValueDescription("vector","the " + keys.getDisplayName() + " for each set of specified atoms");
63 : }
64 426 : }
65 :
66 : template <class T>
67 201 : MultiColvarTemplate<T>::MultiColvarTemplate(const ActionOptions&ao):
68 : Action(ao),
69 : ActionWithVector(ao),
70 201 : mode(0),
71 201 : usepbc(true),
72 201 : wholemolecules(false) {
73 : std::vector<AtomNumber> all_atoms;
74 512 : if( getName()=="POSITION_VECTOR" || getName()=="MASS_VECTOR" || getName()=="CHARGE_VECTOR" ) {
75 118 : parseAtomList( "ATOMS", all_atoms );
76 : }
77 201 : if( all_atoms.size()>0 ) {
78 57 : ablocks.resize(1);
79 57 : ablocks[0].resize( all_atoms.size() );
80 8109 : for(unsigned i=0; i<all_atoms.size(); ++i) {
81 8052 : ablocks[0][i] = i;
82 : }
83 : } else {
84 : std::vector<AtomNumber> t;
85 32324 : for(int i=1;; ++i ) {
86 32324 : T::parseAtomList( i, t, this );
87 32324 : if( t.empty() ) {
88 : break;
89 : }
90 :
91 32180 : if( i==1 ) {
92 144 : ablocks.resize(t.size());
93 : }
94 32180 : if( t.size()!=ablocks.size() ) {
95 : std::string ss;
96 0 : Tools::convert(i,ss);
97 0 : error("ATOMS" + ss + " keyword has the wrong number of atoms");
98 : }
99 98848 : for(unsigned j=0; j<ablocks.size(); ++j) {
100 66668 : ablocks[j].push_back( ablocks.size()*(i-1)+j );
101 66668 : all_atoms.push_back( t[j] );
102 : }
103 32180 : t.resize(0);
104 : }
105 : }
106 201 : if( all_atoms.size()==0 ) {
107 0 : error("No atoms have been specified");
108 : }
109 201 : requestAtoms(all_atoms);
110 201 : if( keywords.exists("NOPBC") ) {
111 201 : bool nopbc=!usepbc;
112 201 : parseFlag("NOPBC",nopbc);
113 201 : usepbc=!nopbc;
114 : }
115 201 : if( keywords.exists("WHOLEMOLECULES") ) {
116 41 : parseFlag("WHOLEMOLECULES",wholemolecules);
117 41 : if( wholemolecules ) {
118 1 : usepbc=false;
119 : }
120 : }
121 201 : if( usepbc ) {
122 167 : log.printf(" using periodic boundary conditions\n");
123 : } else {
124 34 : log.printf(" without periodic boundary conditions\n");
125 : }
126 :
127 : // Setup the values
128 201 : mode = T::getModeAndSetupValues( this );
129 201 : }
130 :
131 : template <class T>
132 1612 : unsigned MultiColvarTemplate<T>::getNumberOfDerivatives() {
133 1612 : return 3*getNumberOfAtoms()+9;
134 : }
135 :
136 : template <class T>
137 22685 : void MultiColvarTemplate<T>::calculate() {
138 22685 : runAllTasks();
139 22685 : }
140 :
141 : template <class T>
142 99 : void MultiColvarTemplate<T>::addValueWithDerivatives( const std::vector<unsigned>& shape ) {
143 99 : std::vector<unsigned> s(1);
144 99 : s[0]=ablocks[0].size();
145 99 : addValue( s );
146 99 : }
147 :
148 : template <class T>
149 318 : void MultiColvarTemplate<T>::addComponentWithDerivatives( const std::string& name, const std::vector<unsigned>& shape ) {
150 318 : std::vector<unsigned> s(1);
151 318 : s[0]=ablocks[0].size();
152 318 : addComponent( name, s );
153 318 : }
154 :
155 : template <class T>
156 28091 : void MultiColvarTemplate<T>::setupStreamedComponents( const std::string& headstr, unsigned& nquants, unsigned& nmat, unsigned& maxcol, unsigned& nbookeeping ) {
157 28091 : if( wholemolecules ) {
158 1 : makeWhole();
159 : }
160 28091 : ActionWithVector::setupStreamedComponents( headstr, nquants, nmat, maxcol, nbookeeping );
161 28091 : }
162 :
163 : template <class T>
164 490140 : void MultiColvarTemplate<T>::performTask( const unsigned& task_index, MultiValue& myvals ) const {
165 : // Retrieve the positions
166 : std::vector<Vector> & fpositions( myvals.getFirstAtomVector() );
167 490140 : if( fpositions.size()!=ablocks.size() ) {
168 38897 : fpositions.resize( ablocks.size() );
169 : }
170 1334749 : for(unsigned i=0; i<ablocks.size(); ++i) {
171 844609 : fpositions[i] = getPosition( ablocks[i][task_index] );
172 : }
173 : // If we are using pbc make whole
174 490140 : if( usepbc ) {
175 233347 : if( fpositions.size()==1 ) {
176 57816 : fpositions[0]=pbcDistance(Vector(0.0,0.0,0.0),getPosition( ablocks[0][task_index] ) );
177 : } else {
178 418629 : for(unsigned j=0; j<fpositions.size()-1; ++j) {
179 : const Vector & first (fpositions[j]);
180 214190 : Vector & second (fpositions[j+1]);
181 214190 : second=first+pbcDistance(first,second);
182 : }
183 : }
184 256793 : } else if( fpositions.size()==1 ) {
185 119580 : fpositions[0]=delta(Vector(0.0,0.0,0.0),getPosition( ablocks[0][task_index] ) );
186 : }
187 : // Retrieve the masses and charges
188 490140 : myvals.resizeTemporyVector(2);
189 : std::vector<double> & mass( myvals.getTemporyVector(0) );
190 : std::vector<double> & charge( myvals.getTemporyVector(1) );
191 490140 : if( mass.size()!=ablocks.size() ) {
192 37865 : mass.resize(ablocks.size());
193 37865 : charge.resize(ablocks.size());
194 : }
195 1334749 : for(unsigned i=0; i<ablocks.size(); ++i) {
196 844609 : mass[i]=getMass( ablocks[i][task_index] );
197 844609 : charge[i]=getCharge( ablocks[i][task_index] );
198 : }
199 : // Make some space to store various things
200 490140 : std::vector<double> values( getNumberOfComponents() );
201 : std::vector<Tensor> & virial( myvals.getFirstAtomVirialVector() );
202 : std::vector<std::vector<Vector> > & derivs( myvals.getFirstAtomDerivativeVector() );
203 490140 : if( derivs.size()!=values.size() ) {
204 39107 : derivs.resize( values.size() );
205 39107 : virial.resize( values.size() );
206 : }
207 1478457 : for(unsigned i=0; i<derivs.size(); ++i) {
208 988317 : if( derivs[i].size()<ablocks.size() ) {
209 71321 : derivs[i].resize( ablocks.size() );
210 : }
211 : }
212 : // Calculate the CVs using the method in the Colvar
213 490140 : T::calculateCV( mode, mass, charge, fpositions, values, derivs, virial, this );
214 1478457 : for(unsigned i=0; i<values.size(); ++i) {
215 988317 : myvals.setValue( getConstPntrToComponent(i)->getPositionInStream(), values[i] );
216 : }
217 : // Finish if there are no derivatives
218 490140 : if( doNotCalculateDerivatives() ) {
219 : return;
220 : }
221 :
222 : // Now transfer the derivatives to the underlying MultiValue
223 868925 : for(unsigned i=0; i<ablocks.size(); ++i) {
224 547918 : unsigned base=3*ablocks[i][task_index];
225 547918 : for(int j=0; j<getNumberOfComponents(); ++j) {
226 908603 : unsigned jval=getConstPntrToComponent(j)->getPositionInStream();
227 908603 : myvals.addDerivative( jval, base + 0, derivs[j][i][0] );
228 908603 : myvals.addDerivative( jval, base + 1, derivs[j][i][1] );
229 908603 : myvals.addDerivative( jval, base + 2, derivs[j][i][2] );
230 : }
231 : // Check for duplicated indices during update to avoid double counting
232 : bool newi=true;
233 787032 : for(unsigned j=0; j<i; ++j) {
234 239114 : if( ablocks[j][task_index]==ablocks[i][task_index] ) {
235 : newi=false;
236 : break;
237 : }
238 : }
239 547918 : if( !newi ) {
240 0 : continue;
241 : }
242 1456521 : for(int j=0; j<getNumberOfComponents(); ++j) {
243 908603 : unsigned jval=getConstPntrToComponent(j)->getPositionInStream();
244 908603 : myvals.updateIndex( jval, base );
245 908603 : myvals.updateIndex( jval, base + 1 );
246 908603 : myvals.updateIndex( jval, base + 2 );
247 : }
248 : }
249 321007 : unsigned nvir=3*getNumberOfAtoms();
250 913213 : for(int j=0; j<getNumberOfComponents(); ++j) {
251 592206 : unsigned jval=getConstPntrToComponent(j)->getPositionInStream();
252 2368824 : for(unsigned i=0; i<3; ++i) {
253 7106472 : for(unsigned k=0; k<3; ++k) {
254 5329854 : myvals.addDerivative( jval, nvir + 3*i + k, virial[j][i][k] );
255 5329854 : myvals.updateIndex( jval, nvir + 3*i + k );
256 : }
257 : }
258 : }
259 : }
260 :
261 : }
262 : }
263 : #endif
|
