Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-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 "CoordinationBase.h"
23 : #include "tools/NeighborList.h"
24 : #include "tools/Communicator.h"
25 : #include "tools/OpenMP.h"
26 :
27 : #include <string>
28 :
29 : using namespace std;
30 :
31 : namespace PLMD {
32 : namespace colvar {
33 :
34 126 : void CoordinationBase::registerKeywords( Keywords& keys ) {
35 126 : Colvar::registerKeywords(keys);
36 126 : keys.addFlag("SERIAL",false,"Perform the calculation in serial - for debug purpose");
37 126 : keys.addFlag("PAIR",false,"Pair only 1st element of the 1st group with 1st element in the second, etc");
38 126 : keys.addFlag("NLIST",false,"Use a neighbour list to speed up the calculation");
39 126 : keys.add("optional","NL_CUTOFF","The cutoff for the neighbour list");
40 126 : keys.add("optional","NL_STRIDE","The frequency with which we are updating the atoms in the neighbour list");
41 126 : keys.add("atoms","GROUPA","First list of atoms");
42 126 : keys.add("atoms","GROUPB","Second list of atoms (if empty, N*(N-1)/2 pairs in GROUPA are counted)");
43 126 : }
44 :
45 124 : CoordinationBase::CoordinationBase(const ActionOptions&ao):
46 : PLUMED_COLVAR_INIT(ao),
47 : pbc(true),
48 : serial(false),
49 : invalidateList(true),
50 124 : firsttime(true)
51 : {
52 :
53 124 : parseFlag("SERIAL",serial);
54 :
55 248 : vector<AtomNumber> ga_lista,gb_lista;
56 124 : parseAtomList("GROUPA",ga_lista);
57 124 : parseAtomList("GROUPB",gb_lista);
58 :
59 124 : bool nopbc=!pbc;
60 124 : parseFlag("NOPBC",nopbc);
61 124 : pbc=!nopbc;
62 :
63 : // pair stuff
64 124 : bool dopair=false;
65 124 : parseFlag("PAIR",dopair);
66 :
67 : // neighbor list stuff
68 124 : bool doneigh=false;
69 124 : double nl_cut=0.0;
70 124 : int nl_st=0;
71 124 : parseFlag("NLIST",doneigh);
72 124 : if(doneigh) {
73 12 : parse("NL_CUTOFF",nl_cut);
74 12 : if(nl_cut<=0.0) error("NL_CUTOFF should be explicitly specified and positive");
75 12 : parse("NL_STRIDE",nl_st);
76 12 : if(nl_st<=0) error("NL_STRIDE should be explicitly specified and positive");
77 : }
78 :
79 124 : addValueWithDerivatives(); setNotPeriodic();
80 124 : if(gb_lista.size()>0) {
81 115 : if(doneigh) nl= new NeighborList(ga_lista,gb_lista,dopair,pbc,getPbc(),nl_cut,nl_st);
82 103 : else nl= new NeighborList(ga_lista,gb_lista,dopair,pbc,getPbc());
83 : } else {
84 9 : if(doneigh) nl= new NeighborList(ga_lista,pbc,getPbc(),nl_cut,nl_st);
85 9 : else nl= new NeighborList(ga_lista,pbc,getPbc());
86 : }
87 :
88 124 : requestAtoms(nl->getFullAtomList());
89 :
90 124 : log.printf(" between two groups of %u and %u atoms\n",static_cast<unsigned>(ga_lista.size()),static_cast<unsigned>(gb_lista.size()));
91 124 : log.printf(" first group:\n");
92 2915 : for(unsigned int i=0; i<ga_lista.size(); ++i) {
93 2791 : if ( (i+1) % 25 == 0 ) log.printf(" \n");
94 2791 : log.printf(" %d", ga_lista[i].serial());
95 : }
96 124 : log.printf(" \n second group:\n");
97 6893 : for(unsigned int i=0; i<gb_lista.size(); ++i) {
98 6769 : if ( (i+1) % 25 == 0 ) log.printf(" \n");
99 6769 : log.printf(" %d", gb_lista[i].serial());
100 : }
101 124 : log.printf(" \n");
102 124 : if(pbc) log.printf(" using periodic boundary conditions\n");
103 0 : else log.printf(" without periodic boundary conditions\n");
104 124 : if(dopair) log.printf(" with PAIR option\n");
105 124 : if(doneigh) {
106 12 : log.printf(" using neighbor lists with\n");
107 12 : log.printf(" update every %d steps and cutoff %f\n",nl_st,nl_cut);
108 124 : }
109 124 : }
110 :
111 248 : CoordinationBase::~CoordinationBase() {
112 124 : delete nl;
113 124 : }
114 :
115 1536 : void CoordinationBase::prepare() {
116 1536 : if(nl->getStride()>0) {
117 78 : if(firsttime || (getStep()%nl->getStride()==0)) {
118 65 : requestAtoms(nl->getFullAtomList());
119 65 : invalidateList=true;
120 65 : firsttime=false;
121 : } else {
122 13 : requestAtoms(nl->getReducedAtomList());
123 13 : invalidateList=false;
124 13 : if(getExchangeStep()) error("Neighbor lists should be updated on exchange steps - choose a NL_STRIDE which divides the exchange stride!");
125 : }
126 78 : if(getExchangeStep()) firsttime=true;
127 : }
128 1536 : }
129 :
130 : // calculator
131 1392 : void CoordinationBase::calculate()
132 : {
133 :
134 1392 : double ncoord=0.;
135 1392 : Tensor virial;
136 1392 : vector<Vector> deriv(getNumberOfAtoms());
137 : // deriv.resize(getPositions().size());
138 :
139 1392 : if(nl->getStride()>0 && invalidateList) {
140 53 : nl->update(getPositions());
141 : }
142 :
143 1392 : unsigned stride=comm.Get_size();
144 1392 : unsigned rank=comm.Get_rank();
145 1392 : if(serial) {
146 0 : stride=1;
147 0 : rank=0;
148 : } else {
149 1392 : stride=comm.Get_size();
150 1392 : rank=comm.Get_rank();
151 : }
152 :
153 1392 : unsigned nt=OpenMP::getNumThreads();
154 1392 : const unsigned nn=nl->size();
155 1392 : if(nt*stride*10>nn) nt=1;
156 :
157 3833 : #pragma omp parallel num_threads(nt)
158 : {
159 2441 : std::vector<Vector> omp_deriv(getPositions().size());
160 2504 : Tensor omp_virial;
161 :
162 2501 : #pragma omp for reduction(+:ncoord) nowait
163 : for(unsigned int i=rank; i<nn; i+=stride) {
164 :
165 5752316 : Vector distance;
166 5751807 : unsigned i0=nl->getClosePair(i).first;
167 5751968 : unsigned i1=nl->getClosePair(i).second;
168 :
169 5797596 : if(getAbsoluteIndex(i0)==getAbsoluteIndex(i1)) continue;
170 :
171 5705632 : if(pbc) {
172 5705632 : distance=pbcDistance(getPosition(i0),getPosition(i1));
173 : } else {
174 0 : distance=delta(getPosition(i0),getPosition(i1));
175 : }
176 :
177 5705127 : double dfunc=0.;
178 5705127 : ncoord += pairing(distance.modulo2(), dfunc,i0,i1);
179 :
180 5705779 : Vector dd(dfunc*distance);
181 5705534 : Tensor vv(dd,distance);
182 5704744 : if(nt>1) {
183 5357620 : omp_deriv[i0]-=dd;
184 5357484 : omp_deriv[i1]+=dd;
185 5357642 : omp_virial-=vv;
186 : } else {
187 347124 : deriv[i0]-=dd;
188 347124 : deriv[i1]+=dd;
189 347124 : virial-=vv;
190 : }
191 :
192 : }
193 5008 : #pragma omp critical
194 2504 : if(nt>1) {
195 294140 : for(int i=0; i<getPositions().size(); i++) deriv[i]+=omp_deriv[i];
196 2224 : virial+=omp_virial;
197 2504 : }
198 : }
199 :
200 1392 : if(!serial) {
201 1392 : comm.Sum(ncoord);
202 1392 : if(!deriv.empty()) comm.Sum(&deriv[0][0],3*deriv.size());
203 1392 : comm.Sum(virial);
204 : }
205 :
206 1392 : for(unsigned i=0; i<deriv.size(); ++i) setAtomsDerivatives(i,deriv[i]);
207 1392 : setValue (ncoord);
208 1392 : setBoxDerivatives (virial);
209 :
210 1392 : }
211 : }
212 2523 : }
|
