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