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 197 : void CoordinationBase::registerKeywords( Keywords& keys ) {
31 197 : Colvar::registerKeywords(keys);
32 394 : keys.addFlag("SERIAL",false,"Perform the calculation in serial - for debug purpose");
33 394 : keys.addFlag("PAIR",false,"Pair only 1st element of the 1st group with 1st element in the second, etc");
34 394 : keys.addFlag("NLIST",false,"Use a neighbor list to speed up the calculation");
35 394 : keys.add("optional","NL_CUTOFF","The cutoff for the neighbor list");
36 394 : keys.add("optional","NL_STRIDE","The frequency with which we are updating the atoms in the neighbor list");
37 394 : keys.add("atoms","GROUPA","First list of atoms");
38 394 : keys.add("atoms","GROUPB","Second list of atoms (if empty, N*(N-1)/2 pairs in GROUPA are counted)");
39 197 : }
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 :
49 388 : parseFlag("SERIAL",serial);
50 :
51 : std::vector<AtomNumber> ga_lista,gb_lista;
52 194 : parseAtomList("GROUPA",ga_lista);
53 194 : parseAtomList("GROUPB",gb_lista);
54 :
55 194 : bool nopbc=!pbc;
56 194 : parseFlag("NOPBC",nopbc);
57 194 : pbc=!nopbc;
58 :
59 : // pair stuff
60 194 : bool dopair=false;
61 194 : parseFlag("PAIR",dopair);
62 :
63 : // neighbor list stuff
64 194 : bool doneigh=false;
65 194 : double nl_cut=0.0;
66 194 : int nl_st=0;
67 194 : parseFlag("NLIST",doneigh);
68 194 : 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 194 : addValueWithDerivatives(); setNotPeriodic();
76 194 : if(gb_lista.size()>0) {
77 207 : if(doneigh) nl=Tools::make_unique<NeighborList>(ga_lista,gb_lista,serial,dopair,pbc,getPbc(),comm,nl_cut,nl_st);
78 318 : 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 194 : requestAtoms(nl->getFullAtomList());
85 :
86 194 : log.printf(" between two groups of %u and %u atoms\n",static_cast<unsigned>(ga_lista.size()),static_cast<unsigned>(gb_lista.size()));
87 194 : log.printf(" first group:\n");
88 5465 : for(unsigned int i=0; i<ga_lista.size(); ++i) {
89 5271 : if ( (i+1) % 25 == 0 ) log.printf(" \n");
90 5271 : log.printf(" %d", ga_lista[i].serial());
91 : }
92 194 : log.printf(" \n second group:\n");
93 10817 : for(unsigned int i=0; i<gb_lista.size(); ++i) {
94 10623 : if ( (i+1) % 25 == 0 ) log.printf(" \n");
95 10623 : log.printf(" %d", gb_lista[i].serial());
96 : }
97 194 : log.printf(" \n");
98 194 : if(pbc) log.printf(" using periodic boundary conditions\n");
99 0 : else log.printf(" without periodic boundary conditions\n");
100 194 : if(dopair) log.printf(" with PAIR option\n");
101 194 : 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 194 : }
106 :
107 194 : CoordinationBase::~CoordinationBase() {
108 : // destructor required to delete forward declared class
109 194 : }
110 :
111 3357 : void CoordinationBase::prepare() {
112 3357 : 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 3357 : }
125 :
126 : // calculator
127 3132 : void CoordinationBase::calculate()
128 : {
129 :
130 3132 : double ncoord=0.;
131 3132 : Tensor virial;
132 3132 : std::vector<Vector> deriv(getNumberOfAtoms());
133 :
134 3132 : if(nl->getStride()>0 && invalidateList) {
135 143 : nl->update(getPositions());
136 : }
137 :
138 : unsigned stride;
139 : unsigned rank;
140 3132 : if(serial) {
141 : stride=1;
142 : rank=0;
143 : } else {
144 3132 : stride=comm.Get_size();
145 3132 : rank=comm.Get_rank();
146 : }
147 :
148 3132 : unsigned nt=OpenMP::getNumThreads();
149 3132 : const unsigned nn=nl->size();
150 3132 : if(nt*stride*10>nn) nt=1;
151 :
152 3132 : #pragma omp parallel num_threads(nt)
153 : {
154 : std::vector<Vector> omp_deriv(getPositions().size());
155 : Tensor omp_virial;
156 :
157 : #pragma omp for reduction(+:ncoord) nowait
158 : for(unsigned int i=rank; i<nn; i+=stride) {
159 :
160 : Vector distance;
161 : unsigned i0=nl->getClosePair(i).first;
162 : unsigned i1=nl->getClosePair(i).second;
163 :
164 : if(getAbsoluteIndex(i0)==getAbsoluteIndex(i1)) continue;
165 :
166 : if(pbc) {
167 : distance=pbcDistance(getPosition(i0),getPosition(i1));
168 : } else {
169 : distance=delta(getPosition(i0),getPosition(i1));
170 : }
171 :
172 : double dfunc=0.;
173 : ncoord += pairing(distance.modulo2(), dfunc,i0,i1);
174 :
175 : Vector dd(dfunc*distance);
176 : Tensor vv(dd,distance);
177 : if(nt>1) {
178 : omp_deriv[i0]-=dd;
179 : omp_deriv[i1]+=dd;
180 : omp_virial-=vv;
181 : } else {
182 : deriv[i0]-=dd;
183 : deriv[i1]+=dd;
184 : virial-=vv;
185 : }
186 :
187 : }
188 : #pragma omp critical
189 : if(nt>1) {
190 : for(unsigned i=0; i<getPositions().size(); i++) deriv[i]+=omp_deriv[i];
191 : virial+=omp_virial;
192 : }
193 : }
194 :
195 3132 : if(!serial) {
196 3132 : comm.Sum(ncoord);
197 3132 : if(!deriv.empty()) comm.Sum(&deriv[0][0],3*deriv.size());
198 3132 : comm.Sum(virial);
199 : }
200 :
201 378180 : for(unsigned i=0; i<deriv.size(); ++i) setAtomsDerivatives(i,deriv[i]);
202 3132 : setValue (ncoord);
203 3132 : setBoxDerivatives (virial);
204 :
205 3132 : }
206 : }
207 : }
|
