Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-2020 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 189 : void CoordinationBase::registerKeywords( Keywords& keys ) {
35 189 : Colvar::registerKeywords(keys);
36 567 : keys.addFlag("SERIAL",false,"Perform the calculation in serial - for debug purpose");
37 567 : keys.addFlag("PAIR",false,"Pair only 1st element of the 1st group with 1st element in the second, etc");
38 567 : keys.addFlag("NLIST",false,"Use a neighbor list to speed up the calculation");
39 756 : keys.add("optional","NL_CUTOFF","The cutoff for the neighbor list");
40 756 : keys.add("optional","NL_STRIDE","The frequency with which we are updating the atoms in the neighbor list");
41 756 : keys.add("atoms","GROUPA","First list of atoms");
42 756 : keys.add("atoms","GROUPB","Second list of atoms (if empty, N*(N-1)/2 pairs in GROUPA are counted)");
43 189 : }
44 :
45 187 : CoordinationBase::CoordinationBase(const ActionOptions&ao):
46 : PLUMED_COLVAR_INIT(ao),
47 : pbc(true),
48 : serial(false),
49 : invalidateList(true),
50 374 : firsttime(true)
51 : {
52 :
53 374 : parseFlag("SERIAL",serial);
54 :
55 : vector<AtomNumber> ga_lista,gb_lista;
56 374 : parseAtomList("GROUPA",ga_lista);
57 374 : parseAtomList("GROUPB",gb_lista);
58 :
59 187 : bool nopbc=!pbc;
60 374 : parseFlag("NOPBC",nopbc);
61 187 : pbc=!nopbc;
62 :
63 : // pair stuff
64 187 : bool dopair=false;
65 374 : parseFlag("PAIR",dopair);
66 :
67 : // neighbor list stuff
68 187 : bool doneigh=false;
69 187 : double nl_cut=0.0;
70 187 : int nl_st=0;
71 374 : parseFlag("NLIST",doneigh);
72 187 : 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 187 : addValueWithDerivatives(); setNotPeriodic();
80 187 : if(gb_lista.size()>0) {
81 224 : if(doneigh) nl.reset( new NeighborList(ga_lista,gb_lista,dopair,pbc,getPbc(),nl_cut,nl_st) );
82 304 : else nl.reset( new NeighborList(ga_lista,gb_lista,dopair,pbc,getPbc()) );
83 : } else {
84 11 : if(doneigh) nl.reset( new NeighborList(ga_lista,pbc,getPbc(),nl_cut,nl_st) );
85 22 : else nl.reset( new NeighborList(ga_lista,pbc,getPbc()) );
86 : }
87 :
88 187 : requestAtoms(nl->getFullAtomList());
89 :
90 561 : log.printf(" between two groups of %u and %u atoms\n",static_cast<unsigned>(ga_lista.size()),static_cast<unsigned>(gb_lista.size()));
91 187 : log.printf(" first group:\n");
92 10902 : for(unsigned int i=0; i<ga_lista.size(); ++i) {
93 5264 : if ( (i+1) % 25 == 0 ) log.printf(" \n");
94 10528 : log.printf(" %d", ga_lista[i].serial());
95 : }
96 187 : log.printf(" \n second group:\n");
97 21606 : for(unsigned int i=0; i<gb_lista.size(); ++i) {
98 10616 : if ( (i+1) % 25 == 0 ) log.printf(" \n");
99 21232 : log.printf(" %d", gb_lista[i].serial());
100 : }
101 187 : log.printf(" \n");
102 187 : if(pbc) log.printf(" using periodic boundary conditions\n");
103 0 : else log.printf(" without periodic boundary conditions\n");
104 187 : if(dopair) log.printf(" with PAIR option\n");
105 187 : 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 187 : }
110 :
111 187 : CoordinationBase::~CoordinationBase() {
112 : // destructor required to delete forward declared class
113 187 : }
114 :
115 2961 : void CoordinationBase::prepare() {
116 2961 : if(nl->getStride()>0) {
117 312 : if(firsttime || (getStep()%nl->getStride()==0)) {
118 322 : requestAtoms(nl->getFullAtomList());
119 161 : invalidateList=true;
120 161 : firsttime=false;
121 : } else {
122 38 : 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 2961 : }
129 :
130 : // calculator
131 2736 : void CoordinationBase::calculate()
132 : {
133 :
134 2736 : double ncoord=0.;
135 2736 : Tensor virial;
136 2736 : vector<Vector> deriv(getNumberOfAtoms());
137 : // deriv.resize(getPositions().size());
138 :
139 2736 : if(nl->getStride()>0 && invalidateList) {
140 125 : nl->update(getPositions());
141 : }
142 :
143 : unsigned stride;
144 : unsigned rank;
145 2736 : if(serial) {
146 : stride=1;
147 : rank=0;
148 : } else {
149 2736 : stride=comm.Get_size();
150 2736 : rank=comm.Get_rank();
151 : }
152 :
153 2736 : unsigned nt=OpenMP::getNumThreads();
154 2736 : const unsigned nn=nl->size();
155 2736 : if(nt*stride*10>nn) nt=1;
156 :
157 7670 : #pragma omp parallel num_threads(nt)
158 : {
159 4934 : std::vector<Vector> omp_deriv(getPositions().size());
160 4934 : Tensor omp_virial;
161 :
162 4934 : #pragma omp for reduction(+:ncoord) nowait
163 : for(unsigned int i=rank; i<nn; i+=stride) {
164 :
165 15375974 : Vector distance;
166 15375974 : unsigned i0=nl->getClosePair(i).first;
167 15375974 : unsigned i1=nl->getClosePair(i).second;
168 :
169 46251146 : if(getAbsoluteIndex(i0)==getAbsoluteIndex(i1)) continue;
170 :
171 15252750 : if(pbc) {
172 30505500 : distance=pbcDistance(getPosition(i0),getPosition(i1));
173 : } else {
174 0 : distance=delta(getPosition(i0),getPosition(i1));
175 : }
176 :
177 15252750 : double dfunc=0.;
178 15252750 : ncoord += pairing(distance.modulo2(), dfunc,i0,i1);
179 :
180 15252750 : Vector dd(dfunc*distance);
181 15252750 : Tensor vv(dd,distance);
182 15252750 : if(nt>1) {
183 28424412 : omp_deriv[i0]-=dd;
184 28424412 : omp_deriv[i1]+=dd;
185 14212206 : omp_virial-=vv;
186 : } else {
187 2081088 : deriv[i0]-=dd;
188 2081088 : deriv[i1]+=dd;
189 1040544 : virial-=vv;
190 : }
191 :
192 : }
193 9868 : #pragma omp critical
194 4934 : if(nt>1) {
195 2612472 : for(unsigned i=0; i<getPositions().size(); i++) deriv[i]+=omp_deriv[i];
196 4396 : virial+=omp_virial;
197 : }
198 : }
199 :
200 2736 : if(!serial) {
201 2736 : comm.Sum(ncoord);
202 5472 : if(!deriv.empty()) comm.Sum(&deriv[0][0],3*deriv.size());
203 2736 : comm.Sum(virial);
204 : }
205 :
206 1492488 : for(unsigned i=0; i<deriv.size(); ++i) setAtomsDerivatives(i,deriv[i]);
207 2736 : setValue (ncoord);
208 2736 : setBoxDerivatives (virial);
209 :
210 2736 : }
211 : }
212 5517 : }
|
