Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-2019 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 : #include <string>
28 :
29 : using namespace std;
30 :
31 : namespace PLMD {
32 : namespace colvar {
33 :
34 163 : void CoordinationBase::registerKeywords( Keywords& keys ) {
35 163 : Colvar::registerKeywords(keys);
36 489 : keys.addFlag("SERIAL",false,"Perform the calculation in serial - for debug purpose");
37 489 : keys.addFlag("PAIR",false,"Pair only 1st element of the 1st group with 1st element in the second, etc");
38 489 : keys.addFlag("NLIST",false,"Use a neighbour list to speed up the calculation");
39 652 : keys.add("optional","NL_CUTOFF","The cutoff for the neighbour list");
40 652 : keys.add("optional","NL_STRIDE","The frequency with which we are updating the atoms in the neighbour list");
41 652 : keys.add("atoms","GROUPA","First list of atoms");
42 652 : keys.add("atoms","GROUPB","Second list of atoms (if empty, N*(N-1)/2 pairs in GROUPA are counted)");
43 163 : }
44 :
45 161 : CoordinationBase::CoordinationBase(const ActionOptions&ao):
46 : PLUMED_COLVAR_INIT(ao),
47 : pbc(true),
48 : serial(false),
49 : invalidateList(true),
50 161 : firsttime(true)
51 : {
52 :
53 322 : parseFlag("SERIAL",serial);
54 :
55 : vector<AtomNumber> ga_lista,gb_lista;
56 322 : parseAtomList("GROUPA",ga_lista);
57 322 : parseAtomList("GROUPB",gb_lista);
58 :
59 161 : bool nopbc=!pbc;
60 322 : parseFlag("NOPBC",nopbc);
61 161 : pbc=!nopbc;
62 :
63 : // pair stuff
64 161 : bool dopair=false;
65 322 : parseFlag("PAIR",dopair);
66 :
67 : // neighbor list stuff
68 161 : bool doneigh=false;
69 161 : double nl_cut=0.0;
70 161 : int nl_st=0;
71 322 : parseFlag("NLIST",doneigh);
72 161 : if(doneigh) {
73 48 : parse("NL_CUTOFF",nl_cut);
74 24 : if(nl_cut<=0.0) error("NL_CUTOFF should be explicitly specified and positive");
75 48 : parse("NL_STRIDE",nl_st);
76 24 : if(nl_st<=0) error("NL_STRIDE should be explicitly specified and positive");
77 : }
78 :
79 161 : addValueWithDerivatives(); setNotPeriodic();
80 161 : if(gb_lista.size()>0) {
81 176 : if(doneigh) nl= new NeighborList(ga_lista,gb_lista,dopair,pbc,getPbc(),nl_cut,nl_st);
82 128 : else nl= new NeighborList(ga_lista,gb_lista,dopair,pbc,getPbc());
83 : } else {
84 9 : if(doneigh) nl= new NeighborList(ga_lista,pbc,getPbc(),nl_cut,nl_st);
85 9 : else nl= new NeighborList(ga_lista,pbc,getPbc());
86 : }
87 :
88 161 : requestAtoms(nl->getFullAtomList());
89 :
90 483 : log.printf(" between two groups of %u and %u atoms\n",static_cast<unsigned>(ga_lista.size()),static_cast<unsigned>(gb_lista.size()));
91 161 : log.printf(" first group:\n");
92 7586 : for(unsigned int i=0; i<ga_lista.size(); ++i) {
93 3632 : if ( (i+1) % 25 == 0 ) log.printf(" \n");
94 7264 : log.printf(" %d", ga_lista[i].serial());
95 : }
96 161 : log.printf(" \n second group:\n");
97 17450 : for(unsigned int i=0; i<gb_lista.size(); ++i) {
98 8564 : if ( (i+1) % 25 == 0 ) log.printf(" \n");
99 17128 : log.printf(" %d", gb_lista[i].serial());
100 : }
101 161 : log.printf(" \n");
102 161 : if(pbc) log.printf(" using periodic boundary conditions\n");
103 0 : else log.printf(" without periodic boundary conditions\n");
104 161 : if(dopair) log.printf(" with PAIR option\n");
105 161 : if(doneigh) {
106 24 : log.printf(" using neighbor lists with\n");
107 24 : log.printf(" update every %d steps and cutoff %f\n",nl_st,nl_cut);
108 : }
109 161 : }
110 :
111 322 : CoordinationBase::~CoordinationBase() {
112 161 : delete nl;
113 161 : }
114 :
115 1842 : void CoordinationBase::prepare() {
116 1842 : if(nl->getStride()>0) {
117 180 : if(firsttime || (getStep()%nl->getStride()==0)) {
118 161 : requestAtoms(nl->getFullAtomList());
119 161 : invalidateList=true;
120 161 : firsttime=false;
121 : } else {
122 19 : requestAtoms(nl->getReducedAtomList());
123 19 : invalidateList=false;
124 19 : if(getExchangeStep()) error("Neighbor lists should be updated on exchange steps - choose a NL_STRIDE which divides the exchange stride!");
125 : }
126 180 : if(getExchangeStep()) firsttime=true;
127 : }
128 1842 : }
129 :
130 : // calculator
131 1662 : void CoordinationBase::calculate()
132 : {
133 :
134 1662 : double ncoord=0.;
135 1662 : Tensor virial;
136 1662 : vector<Vector> deriv(getNumberOfAtoms());
137 : // deriv.resize(getPositions().size());
138 :
139 1662 : if(nl->getStride()>0 && invalidateList) {
140 250 : nl->update(getPositions());
141 : }
142 :
143 : unsigned stride;
144 : unsigned rank;
145 1662 : if(serial) {
146 : stride=1;
147 : rank=0;
148 : } else {
149 1662 : stride=comm.Get_size();
150 1662 : rank=comm.Get_rank();
151 : }
152 :
153 1662 : unsigned nt=OpenMP::getNumThreads();
154 1662 : const unsigned nn=nl->size();
155 1662 : if(nt*stride*10>nn) nt=1;
156 :
157 4508 : #pragma omp parallel num_threads(nt)
158 : {
159 2846 : std::vector<Vector> omp_deriv(getPositions().size());
160 2846 : Tensor omp_virial;
161 :
162 2846 : #pragma omp for reduction(+:ncoord) nowait
163 : for(unsigned int i=rank; i<nn; i+=stride) {
164 :
165 6125234 : Vector distance;
166 6125234 : unsigned i0=nl->getClosePair(i).first;
167 6125234 : unsigned i1=nl->getClosePair(i).second;
168 :
169 18425552 : if(getAbsoluteIndex(i0)==getAbsoluteIndex(i1)) continue;
170 :
171 6075384 : if(pbc) {
172 12150768 : distance=pbcDistance(getPosition(i0),getPosition(i1));
173 : } else {
174 0 : distance=delta(getPosition(i0),getPosition(i1));
175 : }
176 :
177 6075384 : double dfunc=0.;
178 6075384 : ncoord += pairing(distance.modulo2(), dfunc,i0,i1);
179 :
180 6075384 : Vector dd(dfunc*distance);
181 6075384 : Tensor vv(dd,distance);
182 6075384 : if(nt>1) {
183 10763040 : omp_deriv[i0]-=dd;
184 10763040 : omp_deriv[i1]+=dd;
185 5381520 : omp_virial-=vv;
186 : } else {
187 1387728 : deriv[i0]-=dd;
188 1387728 : deriv[i1]+=dd;
189 693864 : virial-=vv;
190 : }
191 :
192 : }
193 5692 : #pragma omp critical
194 2846 : if(nt>1) {
195 1214448 : for(int i=0; i<getPositions().size(); i++) deriv[i]+=omp_deriv[i];
196 2368 : virial+=omp_virial;
197 : }
198 : }
199 :
200 1662 : if(!serial) {
201 1662 : comm.Sum(ncoord);
202 3324 : if(!deriv.empty()) comm.Sum(&deriv[0][0],3*deriv.size());
203 1662 : comm.Sum(virial);
204 : }
205 :
206 741516 : for(unsigned i=0; i<deriv.size(); ++i) setAtomsDerivatives(i,deriv[i]);
207 1662 : setValue (ncoord);
208 1662 : setBoxDerivatives (virial);
209 :
210 1662 : }
211 : }
212 4839 : }
|
