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 7660 : if( keys.style( keys.get(i), "atoms" ) ) keys.reset_style( keys.get(i), "numbered" );
58 : }
59 1164 : if( keys.outputComponentExists(".#!value") ) keys.setValueDescription("vector","the " + keys.getDisplayName() + " for each set of specified atoms");
60 426 : }
61 :
62 : template <class T>
63 201 : MultiColvarTemplate<T>::MultiColvarTemplate(const ActionOptions&ao):
64 : Action(ao),
65 : ActionWithVector(ao),
66 201 : mode(0),
67 201 : usepbc(true),
68 201 : wholemolecules(false)
69 : {
70 : std::vector<AtomNumber> all_atoms;
71 571 : if( getName()=="POSITION_VECTOR" || getName()=="MASS_VECTOR" || getName()=="CHARGE_VECTOR" ) parseAtomList( "ATOMS", all_atoms );
72 201 : if( all_atoms.size()>0 ) {
73 57 : ablocks.resize(1); ablocks[0].resize( all_atoms.size() );
74 8109 : for(unsigned i=0; i<all_atoms.size(); ++i) ablocks[0][i] = i;
75 : } else {
76 : std::vector<AtomNumber> t;
77 32324 : for(int i=1;; ++i ) {
78 32324 : T::parseAtomList( i, t, this );
79 32324 : if( t.empty() ) break;
80 :
81 32180 : if( i==1 ) { ablocks.resize(t.size()); }
82 32180 : if( t.size()!=ablocks.size() ) {
83 0 : std::string ss; Tools::convert(i,ss);
84 0 : error("ATOMS" + ss + " keyword has the wrong number of atoms");
85 : }
86 98848 : for(unsigned j=0; j<ablocks.size(); ++j) {
87 66668 : ablocks[j].push_back( ablocks.size()*(i-1)+j ); all_atoms.push_back( t[j] );
88 : }
89 32180 : t.resize(0);
90 : }
91 : }
92 201 : if( all_atoms.size()==0 ) error("No atoms have been specified");
93 201 : requestAtoms(all_atoms);
94 402 : if( keywords.exists("NOPBC") ) {
95 201 : bool nopbc=!usepbc; parseFlag("NOPBC",nopbc);
96 201 : usepbc=!nopbc;
97 : }
98 402 : if( keywords.exists("WHOLEMOLECULES") ) {
99 41 : parseFlag("WHOLEMOLECULES",wholemolecules);
100 41 : if( wholemolecules ) usepbc=false;
101 : }
102 201 : if( usepbc ) log.printf(" using periodic boundary conditions\n");
103 34 : else log.printf(" without periodic boundary conditions\n");
104 :
105 : // Setup the values
106 201 : mode = T::getModeAndSetupValues( this );
107 201 : }
108 :
109 : template <class T>
110 1612 : unsigned MultiColvarTemplate<T>::getNumberOfDerivatives() {
111 1612 : return 3*getNumberOfAtoms()+9;
112 : }
113 :
114 : template <class T>
115 22685 : void MultiColvarTemplate<T>::calculate() {
116 22685 : runAllTasks();
117 22685 : }
118 :
119 : template <class T>
120 99 : void MultiColvarTemplate<T>::addValueWithDerivatives( const std::vector<unsigned>& shape ) {
121 99 : std::vector<unsigned> s(1); s[0]=ablocks[0].size(); addValue( s );
122 99 : }
123 :
124 : template <class T>
125 318 : void MultiColvarTemplate<T>::addComponentWithDerivatives( const std::string& name, const std::vector<unsigned>& shape ) {
126 318 : std::vector<unsigned> s(1); s[0]=ablocks[0].size(); addComponent( name, s );
127 318 : }
128 :
129 : template <class T>
130 28091 : void MultiColvarTemplate<T>::setupStreamedComponents( const std::string& headstr, unsigned& nquants, unsigned& nmat, unsigned& maxcol, unsigned& nbookeeping ) {
131 28091 : if( wholemolecules ) makeWhole();
132 28091 : ActionWithVector::setupStreamedComponents( headstr, nquants, nmat, maxcol, nbookeeping );
133 28091 : }
134 :
135 : template <class T>
136 490140 : void MultiColvarTemplate<T>::performTask( const unsigned& task_index, MultiValue& myvals ) const {
137 : // Retrieve the positions
138 : std::vector<Vector> & fpositions( myvals.getFirstAtomVector() );
139 490140 : if( fpositions.size()!=ablocks.size() ) fpositions.resize( ablocks.size() );
140 1334749 : for(unsigned i=0; i<ablocks.size(); ++i) fpositions[i] = getPosition( ablocks[i][task_index] );
141 : // If we are using pbc make whole
142 490140 : if( usepbc ) {
143 233347 : if( fpositions.size()==1 ) {
144 57816 : fpositions[0]=pbcDistance(Vector(0.0,0.0,0.0),getPosition( ablocks[0][task_index] ) );
145 : } else {
146 418629 : for(unsigned j=0; j<fpositions.size()-1; ++j) {
147 214190 : const Vector & first (fpositions[j]); Vector & second (fpositions[j+1]);
148 214190 : second=first+pbcDistance(first,second);
149 : }
150 : }
151 256793 : } else if( fpositions.size()==1 ) fpositions[0]=delta(Vector(0.0,0.0,0.0),getPosition( ablocks[0][task_index] ) );
152 : // Retrieve the masses and charges
153 490140 : myvals.resizeTemporyVector(2);
154 : std::vector<double> & mass( myvals.getTemporyVector(0) );
155 : std::vector<double> & charge( myvals.getTemporyVector(1) );
156 490140 : if( mass.size()!=ablocks.size() ) { mass.resize(ablocks.size()); charge.resize(ablocks.size()); }
157 1334749 : for(unsigned i=0; i<ablocks.size(); ++i) { mass[i]=getMass( ablocks[i][task_index] ); charge[i]=getCharge( ablocks[i][task_index] ); }
158 : // Make some space to store various things
159 490140 : std::vector<double> values( getNumberOfComponents() );
160 : std::vector<Tensor> & virial( myvals.getFirstAtomVirialVector() );
161 : std::vector<std::vector<Vector> > & derivs( myvals.getFirstAtomDerivativeVector() );
162 490140 : if( derivs.size()!=values.size() ) { derivs.resize( values.size() ); virial.resize( values.size() ); }
163 1478457 : for(unsigned i=0; i<derivs.size(); ++i) {
164 988317 : if( derivs[i].size()<ablocks.size() ) derivs[i].resize( ablocks.size() );
165 : }
166 : // Calculate the CVs using the method in the Colvar
167 490140 : T::calculateCV( mode, mass, charge, fpositions, values, derivs, virial, this );
168 1478457 : for(unsigned i=0; i<values.size(); ++i) myvals.setValue( getConstPntrToComponent(i)->getPositionInStream(), values[i] );
169 : // Finish if there are no derivatives
170 490140 : if( doNotCalculateDerivatives() ) return;
171 :
172 : // Now transfer the derivatives to the underlying MultiValue
173 868925 : for(unsigned i=0; i<ablocks.size(); ++i) {
174 547918 : unsigned base=3*ablocks[i][task_index];
175 1456521 : for(int j=0; j<getNumberOfComponents(); ++j) {
176 908603 : unsigned jval=getConstPntrToComponent(j)->getPositionInStream();
177 908603 : myvals.addDerivative( jval, base + 0, derivs[j][i][0] );
178 908603 : myvals.addDerivative( jval, base + 1, derivs[j][i][1] );
179 908603 : myvals.addDerivative( jval, base + 2, derivs[j][i][2] );
180 : }
181 : // Check for duplicated indices during update to avoid double counting
182 : bool newi=true;
183 787032 : for(unsigned j=0; j<i; ++j) {
184 239114 : if( ablocks[j][task_index]==ablocks[i][task_index] ) { newi=false; break; }
185 : }
186 547918 : if( !newi ) continue;
187 1456521 : for(int j=0; j<getNumberOfComponents(); ++j) {
188 908603 : unsigned jval=getConstPntrToComponent(j)->getPositionInStream();
189 908603 : myvals.updateIndex( jval, base );
190 908603 : myvals.updateIndex( jval, base + 1 );
191 908603 : myvals.updateIndex( jval, base + 2 );
192 : }
193 : }
194 321007 : unsigned nvir=3*getNumberOfAtoms();
195 913213 : for(int j=0; j<getNumberOfComponents(); ++j) {
196 592206 : unsigned jval=getConstPntrToComponent(j)->getPositionInStream();
197 2368824 : for(unsigned i=0; i<3; ++i) {
198 7106472 : for(unsigned k=0; k<3; ++k) {
199 5329854 : myvals.addDerivative( jval, nvir + 3*i + k, virial[j][i][k] );
200 5329854 : myvals.updateIndex( jval, nvir + 3*i + k );
201 : }
202 : }
203 : }
204 : }
205 :
206 : }
207 : }
208 : #endif
|
