Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-2020 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_volumes_ActionVolume_h
23 : #define __PLUMED_volumes_ActionVolume_h
24 :
25 : #include "core/ActionWithVector.h"
26 : #include "core/ParallelTaskManager.h"
27 : #include "tools/ColvarOutput.h"
28 :
29 : namespace PLMD {
30 : namespace volumes {
31 :
32 : template <class CV>
33 60 : struct VolumeData {
34 : bool not_in;
35 : std::size_t numberOfNonReferenceAtoms;
36 : CV voldata;
37 : #ifdef __PLUMED_HAS_OPENACC
38 : void toACCDevice() const {
39 : #pragma acc enter data copyin(this[0:1],not_in,numberOfNonReferenceAtoms)
40 : voldata.toACCDevice();
41 : }
42 : void removeFromACCDevice() const {
43 : voldata.removeFromACCDevice();
44 : #pragma acc exit data delete(numberOfNonReferenceAtoms,not_in,this[0:1])
45 : }
46 : #endif //__PLUMED_HAS_OPENACC
47 : };
48 :
49 : template <typename precision>
50 : struct VolumeIn {
51 : std::size_t task_index;
52 : const Pbc& pbc;
53 : View<const precision,3> cpos;
54 : View2D<const precision,helpers::dynamic_extent,3> refpos;
55 : VolumeIn( std::size_t t, unsigned nref, precision* p, precision* rp, const Pbc& box ) :
56 1098882 : task_index(t),pbc(box),cpos(p),refpos(rp,nref) {
57 : }
58 : };
59 :
60 : template <typename precision>
61 : class VolumeOut {
62 : private:
63 : class RefderHelper {
64 : private:
65 : precision* derivatives;
66 : public:
67 1097882 : RefderHelper( precision* d ) : derivatives(d) {}
68 : View<precision,3> operator[](std::size_t i) {
69 1142362 : return View<precision,3>( derivatives + 3*(i+1) );
70 : }
71 : };
72 : ColvarOutput<precision> fulldata;
73 : public:
74 : View<precision>& values;
75 : typename ColvarOutput<precision>::VirialHelper& virial;
76 : View<precision,3> derivatives;
77 : RefderHelper refders;
78 : VolumeOut( View<precision>& v, std::size_t nder, precision* d ) :
79 : fulldata(v, nder, d),
80 1098882 : values(fulldata.values),
81 1098882 : virial(fulldata.virial),
82 : derivatives(d), refders(d) {
83 : }
84 : };
85 :
86 : /**
87 : \ingroup INHERIT
88 : This is the abstract base class to use for implementing a new way of defining a particular region of the simulation
89 : box. You can use this to calculate the number of atoms inside that part or the average value of a quantity like the
90 : coordination number inside that part of the cell.
91 : */
92 :
93 : template <class CV, typename myPTM=defaultPTM>
94 : class ActionVolume : public ActionWithVector {
95 : public:
96 : using input_type = VolumeData<CV>;
97 : using mytype=ActionVolume<CV,myPTM>;
98 : using PTM =typename myPTM::template PTM<mytype>;
99 : typedef typename PTM::ParallelActionsInput ParallelActionsInput;
100 : typedef typename PTM::ParallelActionsOutput ParallelActionsOutput;
101 : typedef cvprecision_t<CV> precision;
102 : private:
103 : /// The parallel task manager
104 : PTM taskmanager;
105 : public:
106 : static void registerKeywords( Keywords& keys );
107 : explicit ActionVolume(const ActionOptions&);
108 : unsigned getNumberOfDerivatives() override;
109 : void getAtomicInputData( std::vector<double>& inputdata ) const override ;
110 : void getAtomicInputData( std::vector<float>& inputdata ) const override ;
111 : void calculate() override ;
112 : void applyNonZeroRankForces( std::vector<double>& outforces ) override ;
113 : static void performTask( std::size_t task_index,
114 : const VolumeData<CV>& actiondata,
115 : ParallelActionsInput& input,
116 : ParallelActionsOutput& output );
117 : static int getNumberOfValuesPerTask( std::size_t task_index,
118 : const VolumeData<CV>& actiondata );
119 : static void getForceIndices( std::size_t task_index,
120 : std::size_t colno,
121 : std::size_t ntotal_force,
122 : const VolumeData<CV>& actiondata,
123 : const ParallelActionsInput& input,
124 : ForceIndexHolder force_indices );
125 : };
126 :
127 : template <class CV, typename myPTM>
128 9548 : unsigned ActionVolume<CV, myPTM>::getNumberOfDerivatives() {
129 9548 : return 3*getNumberOfAtoms()+9;
130 : }
131 :
132 : template <class CV, typename myPTM>
133 268 : void ActionVolume<CV, myPTM>::registerKeywords( Keywords& keys ) {
134 268 : ActionWithVector::registerKeywords( keys );
135 268 : PTM::registerKeywords( keys );
136 268 : keys.add("atoms","ATOMS","the group of atoms that you would like to investigate");
137 268 : keys.addFlag("OUTSIDE",false,"calculate quantities for colvars that are on atoms outside the region of interest");
138 536 : keys.setValueDescription("scalar/vector","vector of numbers between 0 and 1 that measure the degree to which each atom is within the volume of interest");
139 268 : CV::registerKeywords( keys );
140 268 : }
141 :
142 : template <class CV, typename myPTM>
143 68 : ActionVolume<CV, myPTM>::ActionVolume(const ActionOptions&ao):
144 : Action(ao),
145 : ActionWithVector(ao),
146 68 : taskmanager(this) {
147 68 : std::vector<AtomNumber> atoms;
148 136 : parseAtomList("ATOMS",atoms);
149 68 : if( atoms.size()==0 ) {
150 0 : error("no atoms were specified");
151 : }
152 68 : log.printf(" examining positions of atoms ");
153 50664 : for(unsigned i=0; i<atoms.size(); ++i) {
154 50596 : log.printf(" %d", atoms[i].serial() );
155 : }
156 68 : log.printf("\n");
157 68 : std::vector<std::size_t> shape(1);
158 68 : shape[0]=atoms.size();
159 :
160 68 : std::vector<AtomNumber> refatoms;
161 68 : CV::parseAtoms( this, refatoms );
162 149 : for(unsigned i=0; i<refatoms.size(); ++i) {
163 81 : atoms.push_back( refatoms[i] );
164 : }
165 68 : requestAtoms( atoms );
166 : VolumeData<CV> actioninput;
167 68 : actioninput.voldata.parseInput(this);
168 :
169 68 : actioninput.numberOfNonReferenceAtoms=shape[0];
170 136 : parseFlag("OUTSIDE",actioninput.not_in);
171 :
172 68 : if( shape[0]==1 ) {
173 2 : ActionWithValue::addValueWithDerivatives();
174 : } else {
175 66 : ActionWithValue::addValue( shape );
176 66 : taskmanager.setupParallelTaskManager( 3*(1+refatoms.size())+9, 3*refatoms.size()+9 );
177 : }
178 68 : setNotPeriodic();
179 68 : getPntrToComponent(0)->setDerivativeIsZeroWhenValueIsZero();
180 :
181 68 : taskmanager.setActionInput( actioninput );
182 68 : }
183 :
184 : template <class CV, typename myPTM>
185 3775 : void ActionVolume<CV, myPTM>::getAtomicInputData( std::vector<double>& inputdata ) const {
186 3775 : if( inputdata.size()!=3*getNumberOfAtoms() ) {
187 3184 : inputdata.resize( 3*getNumberOfAtoms() );
188 : }
189 :
190 702491 : for(unsigned i=0; i<getNumberOfAtoms(); ++i) {
191 698716 : Vector ipos = getPosition(i);
192 2794864 : for(unsigned j=0; j<3; ++j) {
193 2096148 : inputdata[3*i+j] = ipos[j];
194 : }
195 : }
196 3775 : }
197 :
198 : template <class CV, typename myPTM>
199 0 : void ActionVolume<CV, myPTM>::getAtomicInputData( std::vector<float>& inputdata ) const {
200 0 : if( inputdata.size()!=3*getNumberOfAtoms() ) {
201 0 : inputdata.resize( 3*getNumberOfAtoms() );
202 : }
203 :
204 0 : for(unsigned i=0; i<getNumberOfAtoms(); ++i) {
205 0 : Vector ipos = getPosition(i);
206 0 : for(unsigned j=0; j<3; ++j) {
207 0 : inputdata[3*i+j] = ipos[j];
208 : }
209 : }
210 0 : }
211 :
212 : template <class CV, typename myPTM>
213 3775 : void ActionVolume<CV, myPTM>::calculate() {
214 : unsigned k=0;
215 3775 : std::vector<Vector> positions( getNumberOfAtoms()-getPntrToComponent(0)->getNumberOfValues() );
216 16945 : for(unsigned i=getPntrToComponent(0)->getNumberOfValues(); i<getNumberOfAtoms(); ++i) {
217 13170 : positions[k] = getPosition( i );
218 13170 : k++;
219 : }
220 3775 : taskmanager.getActionInput().voldata.setupRegions( this, getPbc(), positions );
221 :
222 3775 : if( getPntrToComponent(0)->getRank()==0 ) {
223 3120 : std::size_t nref = getNumberOfAtoms() - 1;
224 3120 : std::vector<double> posvec;
225 3120 : getAtomicInputData( posvec );
226 3120 : std::vector<double> deriv( getNumberOfDerivatives() );
227 3120 : std::vector<double> val(1);
228 : View<double> valview(val.data(),1);
229 3120 : VolumeOut output( valview, getNumberOfDerivatives(), deriv.data() );
230 3120 : CV::calculateNumberInside( VolumeIn( 0, nref, posvec.data(), posvec.data()+3, getPbc() ),
231 0 : taskmanager.getActionInput().voldata,
232 : output );
233 3120 : if( taskmanager.getActionInput().not_in ) {
234 0 : val[0] = 1.0 - val[0];
235 0 : for(unsigned i=0; i<deriv.size(); ++i) {
236 0 : deriv[i] *= -1;
237 : }
238 : }
239 3120 : Value* v = getPntrToComponent(0);
240 3120 : v->set( val[0] );
241 78000 : for(unsigned i=0; i<deriv.size(); ++i) {
242 74880 : v->addDerivative( i, deriv[i] );
243 : }
244 0 : } else {
245 655 : taskmanager.runAllTasks();
246 : }
247 3775 : }
248 :
249 : template <class CV, typename myPTM>
250 404 : void ActionVolume<CV, myPTM>::applyNonZeroRankForces( std::vector<double>& outforces ) {
251 404 : taskmanager.applyForces( outforces );
252 404 : }
253 :
254 : template <class CV, typename myPTM>
255 1095762 : void ActionVolume<CV, myPTM>::performTask( std::size_t task_index,
256 : const VolumeData<CV>& actiondata,
257 : ParallelActionsInput& input,
258 : ParallelActionsOutput& output ) {
259 1095762 : std::size_t nref = output.derivatives.size()/3 - 4; // This is the number of reference atoms
260 : VolumeOut volout( output.values, output.derivatives.size(), output.derivatives.data() );
261 2191524 : CV::calculateNumberInside( VolumeIn( task_index, nref, input.inputdata+3*task_index, input.inputdata+3*actiondata.numberOfNonReferenceAtoms, *input.pbc ),
262 1095762 : actiondata.voldata,
263 : volout );
264 :
265 1095762 : if( actiondata.not_in ) {
266 8000 : output.values[0] = 1.0 - output.values[0];
267 8000 : if( input.noderiv ) {
268 4000 : return;
269 : }
270 64000 : for(unsigned i=0; i<output.derivatives.size(); ++i) {
271 60000 : output.derivatives[i] *= -1;
272 : }
273 : }
274 : }
275 :
276 : template <class CV, typename myPTM>
277 413336 : int ActionVolume<CV, myPTM>::getNumberOfValuesPerTask( std::size_t task_index,
278 : const VolumeData<CV>& actiondata ) {
279 413336 : return 1;
280 : }
281 :
282 : template<class CV, typename myPTM>
283 413336 : void ActionVolume<CV, myPTM>::getForceIndices( std::size_t task_index,
284 : std::size_t colno,
285 : std::size_t ntotal_force,
286 : const VolumeData<CV>& actiondata,
287 : const ParallelActionsInput& input,
288 : ForceIndexHolder force_indices ) {
289 413336 : std::size_t base = 3*task_index;
290 413336 : force_indices.indices[0][0] = base;
291 413336 : force_indices.indices[0][1] = base + 1;
292 413336 : force_indices.indices[0][2] = base + 2;
293 413336 : force_indices.threadsafe_derivatives_end[0]=3;
294 : std::size_t m=3;
295 5383868 : for(unsigned n=3*actiondata.numberOfNonReferenceAtoms; n<ntotal_force; ++n) {
296 4970532 : force_indices.indices[0][m] = n;
297 4970532 : ++m;
298 : }
299 413336 : force_indices.tot_indices[0] = m;
300 413336 : }
301 :
302 : }
303 : }
304 : #endif
|
