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