Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-2018 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 : #include "Colvar.h"
23 : #include "tools/OpenMP.h"
24 : #include <vector>
25 : #include <string>
26 :
27 : using namespace std;
28 : namespace PLMD {
29 :
30 722 : Colvar::Colvar(const ActionOptions&ao):
31 : Action(ao),
32 : ActionAtomistic(ao),
33 : ActionWithValue(ao),
34 722 : isEnergy(false)
35 : {
36 722 : }
37 :
38 729 : void Colvar::registerKeywords( Keywords& keys ) {
39 729 : Action::registerKeywords( keys );
40 729 : ActionWithValue::registerKeywords( keys );
41 729 : ActionAtomistic::registerKeywords( keys );
42 729 : keys.addFlag("NOPBC",false,"ignore the periodic boundary conditions when calculating distances");
43 729 : }
44 :
45 799 : void Colvar::requestAtoms(const vector<AtomNumber> & a) {
46 799 : plumed_massert(!isEnergy,"request atoms should not be called if this is energy");
47 : // Tell actionAtomistic what atoms we are getting
48 799 : ActionAtomistic::requestAtoms(a);
49 : // Resize the derivatives of all atoms
50 799 : for(int i=0; i<getNumberOfComponents(); ++i) getPntrToComponent(i)->resizeDerivatives(3*a.size()+9);
51 799 : }
52 :
53 43557 : void Colvar::apply() {
54 43557 : vector<Vector>& f(modifyForces());
55 43557 : Tensor& v(modifyVirial());
56 43557 : const unsigned nat=getNumberOfAtoms();
57 43557 : const unsigned ncp=getNumberOfComponents();
58 43557 : const unsigned fsz=f.size();
59 :
60 43557 : unsigned stride=1;
61 43557 : unsigned rank=0;
62 43557 : if(ncp>4*comm.Get_size()) {
63 838 : stride=comm.Get_size();
64 838 : rank=comm.Get_rank();
65 : }
66 :
67 43557 : unsigned nt=OpenMP::getNumThreads();
68 43557 : if(nt>ncp/(2.*stride)) nt=1;
69 :
70 43557 : if(!isEnergy) {
71 88304 : #pragma omp parallel num_threads(nt)
72 : {
73 44797 : vector<Vector> omp_f(fsz);
74 44807 : Tensor omp_v;
75 89623 : vector<double> forces(3*nat+9);
76 44811 : #pragma omp for
77 : for(unsigned i=rank; i<ncp; i+=stride) {
78 84957 : if(getPntrToComponent(i)->applyForce(forces)) {
79 4756208 : for(unsigned j=0; j<nat; ++j) {
80 4729216 : omp_f[j][0]+=forces[3*j+0];
81 4729218 : omp_f[j][1]+=forces[3*j+1];
82 4729220 : omp_f[j][2]+=forces[3*j+2];
83 : }
84 26992 : omp_v(0,0)+=forces[3*nat+0];
85 26977 : omp_v(0,1)+=forces[3*nat+1];
86 26990 : omp_v(0,2)+=forces[3*nat+2];
87 26990 : omp_v(1,0)+=forces[3*nat+3];
88 26992 : omp_v(1,1)+=forces[3*nat+4];
89 26984 : omp_v(1,2)+=forces[3*nat+5];
90 26992 : omp_v(2,0)+=forces[3*nat+6];
91 26993 : omp_v(2,1)+=forces[3*nat+7];
92 26991 : omp_v(2,2)+=forces[3*nat+8];
93 : }
94 : }
95 89630 : #pragma omp critical
96 : {
97 555544 : for(unsigned j=0; j<nat; ++j) f[j]+=omp_f[j];
98 44815 : v+=omp_v;
99 44815 : }
100 : }
101 :
102 43507 : if(ncp>4*comm.Get_size()) {
103 838 : if(fsz>0) comm.Sum(&f[0][0],3*fsz);
104 838 : comm.Sum(&v[0][0],9);
105 : }
106 :
107 50 : } else if( isEnergy ) {
108 50 : vector<double> forces(1);
109 50 : if(getPntrToComponent(0)->applyForce(forces)) modifyForceOnEnergy()+=forces[0];
110 : }
111 43557 : }
112 :
113 124186 : void Colvar::setBoxDerivativesNoPbc(Value* v) {
114 124186 : Tensor virial;
115 124186 : unsigned nat=getNumberOfAtoms();
116 46880731 : for(unsigned i=0; i<nat; i++) virial-=Tensor(getPosition(i),
117 : Vector(v->getDerivative(3*i+0),
118 15585515 : v->getDerivative(3*i+1),
119 31171030 : v->getDerivative(3*i+2)));
120 124186 : setBoxDerivatives(v,virial);
121 124186 : }
122 2523 : }
|