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 234 : void CoordinationBase::registerKeywords( Keywords& keys ) {
31 234 : Colvar::registerKeywords(keys);
32 234 : keys.addFlag("SERIAL",false,"Perform the calculation in serial - for debug purpose");
33 234 : keys.addFlag("PAIR",false,"Pair only 1st element of the 1st group with 1st element in the second, etc");
34 234 : keys.addFlag("NLIST",false,"Use a neighbor list to speed up the calculation");
35 234 : keys.add("optional","NL_CUTOFF","The cutoff for the neighbor list");
36 234 : keys.add("optional","NL_STRIDE","The frequency with which we are updating the atoms in the neighbor list");
37 234 : keys.add("atoms","GROUPA","First list of atoms");
38 234 : keys.add("atoms","GROUPB","Second list of atoms (if empty, N*(N-1)/2 pairs in GROUPA are counted)");
39 234 : }
40 :
41 228 : CoordinationBase::CoordinationBase(const ActionOptions&ao):
42 : PLUMED_COLVAR_INIT(ao),
43 228 : pbc(true),
44 228 : serial(false),
45 228 : invalidateList(true),
46 228 : firsttime(true) {
47 :
48 456 : parseFlag("SERIAL",serial);
49 :
50 : std::vector<AtomNumber> ga_lista,gb_lista;
51 228 : parseAtomList("GROUPA",ga_lista);
52 228 : parseAtomList("GROUPB",gb_lista);
53 :
54 228 : bool nopbc=!pbc;
55 228 : parseFlag("NOPBC",nopbc);
56 228 : pbc=!nopbc;
57 :
58 : // pair stuff
59 228 : bool dopair=false;
60 228 : parseFlag("PAIR",dopair);
61 :
62 : // neighbor list stuff
63 228 : bool doneigh=false;
64 228 : double nl_cut=0.0;
65 228 : int nl_st=0;
66 228 : parseFlag("NLIST",doneigh);
67 228 : 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 228 : addValueWithDerivatives();
79 228 : setNotPeriodic();
80 228 : if(gb_lista.size()>0) {
81 217 : 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 386 : 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 228 : requestAtoms(nl->getFullAtomList());
95 :
96 228 : log.printf(" between two groups of %u and %u atoms\n",static_cast<unsigned>(ga_lista.size()),static_cast<unsigned>(gb_lista.size()));
97 228 : log.printf(" first group:\n");
98 5535 : for(unsigned int i=0; i<ga_lista.size(); ++i) {
99 5307 : if ( (i+1) % 25 == 0 ) {
100 156 : log.printf(" \n");
101 : }
102 5307 : log.printf(" %d", ga_lista[i].serial());
103 : }
104 228 : log.printf(" \n second group:\n");
105 10971 : for(unsigned int i=0; i<gb_lista.size(); ++i) {
106 10743 : if ( (i+1) % 25 == 0 ) {
107 354 : log.printf(" \n");
108 : }
109 10743 : log.printf(" %d", gb_lista[i].serial());
110 : }
111 228 : log.printf(" \n");
112 228 : if(pbc) {
113 228 : log.printf(" using periodic boundary conditions\n");
114 : } else {
115 0 : log.printf(" without periodic boundary conditions\n");
116 : }
117 228 : if(dopair) {
118 2 : log.printf(" with PAIR option\n");
119 : }
120 228 : 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 228 : }
125 :
126 228 : CoordinationBase::~CoordinationBase() {
127 : // destructor required to delete forward declared class
128 228 : }
129 :
130 4012 : void CoordinationBase::prepare() {
131 4012 : 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 4012 : }
148 :
149 : // calculator
150 3787 : void CoordinationBase::calculate() {
151 :
152 3787 : double ncoord=0.;
153 3787 : Tensor virial;
154 3787 : std::vector<Vector> deriv(getNumberOfAtoms());
155 :
156 3787 : if(nl->getStride()>0 && invalidateList) {
157 143 : nl->update(getPositions());
158 : }
159 :
160 : unsigned stride;
161 : unsigned rank;
162 3787 : if(serial) {
163 : stride=1;
164 : rank=0;
165 : } else {
166 3787 : stride=comm.Get_size();
167 3787 : rank=comm.Get_rank();
168 : }
169 :
170 3787 : unsigned nt=OpenMP::getNumThreads();
171 3787 : const unsigned nn=nl->size();
172 3787 : if(nt*stride*10>nn) {
173 : nt=1;
174 : }
175 :
176 3787 : const unsigned elementsPerRank = std::ceil(double(nn)/stride);
177 3787 : const unsigned int start= rank*elementsPerRank;
178 3787 : const unsigned int end = ((start + elementsPerRank)< nn)?(start + elementsPerRank): nn;
179 :
180 3787 : #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 3787 : if(!serial) {
228 3787 : comm.Sum(ncoord);
229 3787 : if(!deriv.empty()) {
230 3787 : comm.Sum(&deriv[0][0],3*deriv.size());
231 : }
232 3787 : comm.Sum(virial);
233 : }
234 :
235 380565 : for(unsigned i=0; i<deriv.size(); ++i) {
236 376778 : setAtomsDerivatives(i,deriv[i]);
237 : }
238 3787 : setValue (ncoord);
239 3787 : setBoxDerivatives (virial);
240 :
241 3787 : }
242 : }
243 : }
|
