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