Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-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 : #include "CoordinationBase.h"
23 : #include "tools/NeighborList.h"
24 : #include "tools/Communicator.h"
25 : #include "tools/OpenMP.h"
26 :
27 : namespace PLMD {
28 : namespace colvar {
29 :
30 240 : void CoordinationBase::registerKeywords( Keywords& keys ) {
31 240 : Colvar::registerKeywords(keys);
32 240 : keys.addFlag("SERIAL",false,"Perform the calculation in serial - for debug purpose");
33 240 : keys.addFlag("PAIR",false,"Pair only 1st element of the 1st group with 1st element in the second, etc");
34 240 : keys.addFlag("NLIST",false,"Use a neighbor list to speed up the calculation");
35 240 : keys.addFlag("NLISTCELLS",false,"Use a neighbor list to speed up the calculation - use the cell list implementation instead of the classical one");
36 240 : keys.add("optional","NL_CUTOFF","The cutoff for the neighbor list");
37 240 : keys.add("optional","NL_STRIDE","The frequency with which we are updating the atoms in the neighbor list");
38 240 : keys.add("atoms","GROUPA","First list of atoms");
39 240 : keys.add("atoms","GROUPB","Second list of atoms (if empty, N*(N-1)/2 pairs in GROUPA are counted)");
40 240 : }
41 :
42 234 : CoordinationBase::CoordinationBase(const ActionOptions&ao):
43 : PLUMED_COLVAR_INIT(ao),
44 234 : pbc(true),
45 234 : serial(false),
46 234 : invalidateList(true),
47 234 : firsttime(true) {
48 :
49 468 : parseFlag("SERIAL",serial);
50 :
51 : std::vector<AtomNumber> ga_lista,gb_lista;
52 234 : parseAtomList("GROUPA",ga_lista);
53 234 : parseAtomList("GROUPB",gb_lista);
54 :
55 234 : bool nopbc=!pbc;
56 234 : parseFlag("NOPBC",nopbc);
57 234 : pbc=!nopbc;
58 :
59 : // pair stuff
60 234 : bool dopair=false;
61 234 : parseFlag("PAIR",dopair);
62 :
63 : // neighbor list stuff
64 234 : bool doneigh_classic=false;
65 234 : double nl_cut=0.0;
66 234 : int nl_st=0;
67 234 : parseFlag("NLIST",doneigh_classic);
68 234 : bool doneighcells=false;
69 234 : parseFlag("NLISTCELLS",doneighcells);
70 : //temporary message
71 234 : plumed_assert(!(doneighcells && doneigh_classic)) << "Please activate only one of the two version of the NL";
72 234 : plumed_assert(!(doneighcells && dopair)) << "Pair is not compatible with the CELLS implementation of the NL";
73 234 : bool doneigh=doneighcells||doneigh_classic;
74 : if(doneigh) {
75 30 : parse("NL_CUTOFF",nl_cut);
76 30 : if(nl_cut<=0.0) {
77 0 : error("NL_CUTOFF should be explicitly specified and positive");
78 : }
79 30 : parse("NL_STRIDE",nl_st);
80 30 : if(nl_st<=0) {
81 0 : error("NL_STRIDE should be explicitly specified and positive");
82 : }
83 : }
84 :
85 234 : addValueWithDerivatives();
86 234 : setNotPeriodic();
87 234 : if(gb_lista.size()>0) {
88 223 : if(doneigh) {
89 60 : nl=Tools::make_unique<NeighborList>(ga_lista,gb_lista,serial,dopair,pbc,getPbc(),comm,nl_cut,nl_st,doneighcells);
90 : } else {
91 386 : nl=Tools::make_unique<NeighborList>(ga_lista,gb_lista,serial,dopair,pbc,getPbc(),comm);
92 : }
93 : } else {
94 11 : if(doneigh) {
95 0 : nl=Tools::make_unique<NeighborList>(ga_lista,serial,pbc,getPbc(),comm,nl_cut,nl_st,doneighcells);
96 : } else {
97 22 : nl=Tools::make_unique<NeighborList>(ga_lista,serial,pbc,getPbc(),comm);
98 : }
99 : }
100 :
101 234 : requestAtoms(nl->getFullAtomList());
102 :
103 234 : log.printf(" between two groups of %u and %u atoms\n",static_cast<unsigned>(ga_lista.size()),static_cast<unsigned>(gb_lista.size()));
104 234 : log.printf(" first group:\n");
105 6189 : for(unsigned int i=0; i<ga_lista.size(); ++i) {
106 5955 : if ( (i+1) % 25 == 0 ) {
107 180 : log.printf(" \n");
108 : }
109 5955 : log.printf(" %d", ga_lista[i].serial());
110 : }
111 234 : log.printf(" \n second group:\n");
112 11625 : for(unsigned int i=0; i<gb_lista.size(); ++i) {
113 11391 : if ( (i+1) % 25 == 0 ) {
114 378 : log.printf(" \n");
115 : }
116 11391 : log.printf(" %d", gb_lista[i].serial());
117 : }
118 234 : log.printf(" \n");
119 234 : if(pbc) {
120 234 : log.printf(" using periodic boundary conditions\n");
121 : } else {
122 0 : log.printf(" without periodic boundary conditions\n");
123 : }
124 234 : if(dopair) {
125 2 : log.printf(" with PAIR option\n");
126 : }
127 234 : if(doneigh) {
128 30 : log.printf(" using neighbor lists with\n");
129 30 : log.printf(" update every %d steps and cutoff %f\n",nl_st,nl_cut);
130 : }
131 234 : }
132 :
133 234 : CoordinationBase::~CoordinationBase() {
134 : // destructor required to delete forward declared class
135 234 : }
136 :
137 4174 : void CoordinationBase::prepare() {
138 4174 : if(nl->getStride()>0) {
139 378 : if(firsttime || (getStep()%nl->getStride()==0)) {
140 221 : requestAtoms(nl->getFullAtomList());
141 221 : invalidateList=true;
142 221 : firsttime=false;
143 : } else {
144 157 : requestAtoms(nl->getReducedAtomList());
145 157 : invalidateList=false;
146 157 : if(getExchangeStep()) {
147 0 : error("Neighbor lists should be updated on exchange steps - choose a NL_STRIDE which divides the exchange stride!");
148 : }
149 : }
150 378 : if(getExchangeStep()) {
151 36 : firsttime=true;
152 : }
153 : }
154 4174 : }
155 :
156 : // calculator
157 3949 : void CoordinationBase::calculate() {
158 :
159 3949 : double ncoord=0.;
160 : Tensor virial;
161 3949 : std::vector<Vector> deriv(getNumberOfAtoms());
162 :
163 3949 : if(nl->getStride()>0 && invalidateList) {
164 185 : nl->update(getPositions());
165 : }
166 :
167 : unsigned stride;
168 : unsigned rank;
169 3949 : if(serial) {
170 : stride=1;
171 : rank=0;
172 : } else {
173 3949 : stride=comm.Get_size();
174 3949 : rank=comm.Get_rank();
175 : }
176 :
177 3949 : unsigned nt=OpenMP::getNumThreads();
178 3949 : const unsigned nn=nl->size();
179 3949 : if(nt*stride*10>nn) {
180 : nt=1;
181 : }
182 :
183 3949 : const unsigned elementsPerRank = std::ceil(double(nn)/stride);
184 3949 : const unsigned int start= rank*elementsPerRank;
185 3949 : const unsigned int end = ((start + elementsPerRank)< nn)?(start + elementsPerRank): nn;
186 :
187 3949 : #pragma omp parallel num_threads(nt)
188 : {
189 : std::vector<Vector> omp_deriv(getPositions().size());
190 : Tensor omp_virial;
191 :
192 : #pragma omp for reduction(+:ncoord) nowait
193 : for(unsigned int i=start; i<end; ++i) {
194 :
195 : Vector distance;
196 : const unsigned i0=nl->getClosePair(i).first;
197 : const unsigned i1=nl->getClosePair(i).second;
198 :
199 : if(getAbsoluteIndex(i0)==getAbsoluteIndex(i1)) {
200 : continue;
201 : }
202 :
203 : if(pbc) {
204 : distance=pbcDistance(getPosition(i0),getPosition(i1));
205 : } else {
206 : distance=delta(getPosition(i0),getPosition(i1));
207 : }
208 :
209 : double dfunc=0.;
210 : ncoord += pairing(distance.modulo2(), dfunc,i0,i1);
211 :
212 : Vector dd(dfunc*distance);
213 : Tensor vv(dd,distance);
214 : if(nt>1) {
215 : omp_deriv[i0]-=dd;
216 : omp_deriv[i1]+=dd;
217 : omp_virial-=vv;
218 : } else {
219 : deriv[i0]-=dd;
220 : deriv[i1]+=dd;
221 : virial-=vv;
222 : }
223 :
224 : }
225 : #pragma omp critical
226 : if(nt>1) {
227 : for(unsigned i=0; i<getPositions().size(); i++) {
228 : deriv[i]+=omp_deriv[i];
229 : }
230 : virial+=omp_virial;
231 : }
232 : }
233 :
234 3949 : if(!serial) {
235 3949 : comm.Sum(ncoord);
236 3949 : if(!deriv.empty()) {
237 3949 : comm.Sum(&deriv[0][0],3*deriv.size());
238 : }
239 3949 : comm.Sum(virial);
240 : }
241 :
242 402759 : for(unsigned i=0; i<deriv.size(); ++i) {
243 398810 : setAtomsDerivatives(i,deriv[i]);
244 : }
245 3949 : setValue (ncoord);
246 3949 : setBoxDerivatives (virial);
247 :
248 3949 : }
249 : }
250 : }
|
