# | Line 1 | Line 1 | |
---|---|---|
1 | #include "Integrator.hpp" | |
2 | #include "simError.h" | |
3 | < | |
3 | > | #include <cmath> |
4 | template<typename T> ZConstraint<T>::ZConstraint(SimInfo* theInfo, ForceFields* the_ff) | |
5 | < | : T(theInfo, the_ff), fz(NULL), indexOfZConsMols(NULL) |
5 | > | : T(theInfo, the_ff), fz(NULL), |
6 | > | indexOfZConsMols(NULL) |
7 | { | |
8 | ||
9 | //get properties from SimInfo | |
10 | GenericData* data; | |
11 | < | IndexData* index; |
11 | > | ZConsParaData* zConsParaData; |
12 | DoubleData* sampleTime; | |
13 | + | DoubleData* tolerance; |
14 | StringData* filename; | |
15 | < | |
14 | < | |
15 | < | data = info->getProperty("zconsindex"); |
16 | < | if(!data) { |
15 | > | double COM[3]; |
16 | ||
17 | < | sprintf( painCave.errMsg, |
18 | < | "ZConstraint error: If you use an ZConstraint\n" |
19 | < | " , you must set index of z-constraint molecules.\n"); |
20 | < | painCave.isFatal = 1; |
21 | < | simError(); |
23 | < | } |
24 | < | else{ |
25 | < | index = dynamic_cast<IndexData*>(data); |
26 | < | |
27 | < | if(!index){ |
17 | > | //by default, the direction of constraint is z |
18 | > | // 0 --> x |
19 | > | // 1 --> y |
20 | > | // 2 --> z |
21 | > | whichDirection = 2; |
22 | ||
23 | < | sprintf( painCave.errMsg, |
24 | < | "ZConstraint error: Can not get property from SimInfo\n"); |
25 | < | painCave.isFatal = 1; |
26 | < | simError(); |
27 | < | |
34 | < | } |
35 | < | else{ |
36 | < | |
37 | < | indexOfAllZConsMols = index->getIndexData(); |
38 | < | |
39 | < | //the maximum value of index is the last one(we sorted the index data in SimSetup.cpp) |
40 | < | int maxIndex; |
41 | < | int totalNumMol; |
42 | < | |
43 | < | maxIndex = indexOfAllZConsMols[indexOfAllZConsMols.size() - 1]; |
23 | > | //estimate the force constant of harmonical potential |
24 | > | double Kb = 1.986E-3 ; //in kcal/K |
25 | > | |
26 | > | double halfOfLargestBox = max(info->boxL[0], max(info->boxL[1], info->boxL[2])) /2; |
27 | > | zForceConst = Kb * info->target_temp /(halfOfLargestBox * halfOfLargestBox); |
28 | ||
29 | < | #ifndef IS_MPI |
30 | < | totalNumMol = nMols; |
47 | < | #else |
48 | < | totalNumMol = mpiSim->getTotNmol(); |
49 | < | #endif |
50 | < | |
51 | < | if(maxIndex > totalNumMol - 1){ |
52 | < | sprintf( painCave.errMsg, |
53 | < | "ZConstraint error: index is out of range\n"); |
54 | < | painCave.isFatal = 1; |
55 | < | simError(); |
56 | < | |
57 | < | } |
58 | < | |
59 | < | } |
60 | < | |
61 | < | } |
29 | > | //retrieve sample time of z-contraint |
30 | > | data = info->getProperty(ZCONSTIME_ID); |
31 | ||
63 | – | //retrive sample time of z-contraint |
64 | – | data = info->getProperty("zconstime"); |
65 | – | |
32 | if(!data) { | |
33 | ||
34 | sprintf( painCave.errMsg, | |
# | Line 89 | Line 55 | template<typename T> ZConstraint<T>::ZConstraint(SimIn | |
55 | ||
56 | } | |
57 | ||
58 | < | |
59 | < | //retrive output filename of z force |
94 | < | data = info->getProperty("zconsfilename"); |
58 | > | //retrieve output filename of z force |
59 | > | data = info->getProperty(ZCONSFILENAME_ID); |
60 | if(!data) { | |
61 | ||
62 | ||
# | Line 121 | Line 86 | template<typename T> ZConstraint<T>::ZConstraint(SimIn | |
86 | ||
87 | } | |
88 | ||
89 | < | |
90 | < | //calculate reference z coordinate for z-constraint molecules |
126 | < | double totalMass_local; |
127 | < | double totalMass; |
128 | < | double totalMZ_local; |
129 | < | double totalMZ; |
130 | < | double massOfUncons_local; |
131 | < | double massOfCurMol; |
132 | < | double COM[3]; |
89 | > | //retrieve tolerance for z-constraint molecuels |
90 | > | data = info->getProperty(ZCONSTOL_ID); |
91 | ||
92 | < | totalMass_local = 0; |
93 | < | totalMass = 0; |
94 | < | totalMZ_local = 0; |
95 | < | totalMZ = 0; |
96 | < | massOfUncons_local = 0; |
97 | < | |
140 | < | |
141 | < | for(int i = 0; i < nMols; i++){ |
142 | < | massOfCurMol = molecules[i].getTotalMass(); |
143 | < | molecules[i].getCOM(COM); |
144 | < | |
145 | < | totalMass_local += massOfCurMol; |
146 | < | totalMZ_local += massOfCurMol * COM[2]; |
147 | < | |
148 | < | if(isZConstraintMol(&molecules[i]) == -1){ |
149 | < | |
150 | < | massOfUncons_local += massOfCurMol; |
151 | < | } |
152 | < | |
92 | > | if(!data) { |
93 | > | |
94 | > | sprintf( painCave.errMsg, |
95 | > | "ZConstraint error: can not get tolerance \n"); |
96 | > | painCave.isFatal = 1; |
97 | > | simError(); |
98 | } | |
99 | + | else{ |
100 | ||
101 | < | |
102 | < | #ifdef IS_MPI |
103 | < | MPI_Allreduce(&totalMass_local, &totalMass, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
158 | < | MPI_Allreduce(&totalMZ_local, &totalMZ, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
159 | < | MPI_Allreduce(&massOfUncons_local, &totalMassOfUncons, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
160 | < | #else |
161 | < | totalMass = totalMass_local; |
162 | < | totalMZ = totalMZ_local; |
163 | < | totalMassOfUncons = massOfUncons_local; |
164 | < | #endif |
101 | > | tolerance = dynamic_cast<DoubleData*>(data); |
102 | > | |
103 | > | if(!tolerance){ |
104 | ||
105 | < | double zsys; |
106 | < | zsys = totalMZ / totalMass; |
107 | < | |
108 | < | #ifndef IS_MPI |
109 | < | for(int i = 0; i < nMols; i++){ |
171 | < | |
172 | < | if(isZConstraintMol(&molecules[i]) > -1 ){ |
173 | < | molecules[i].getCOM(COM); |
174 | < | allRefZ.push_back(COM[2] - zsys); |
105 | > | sprintf( painCave.errMsg, |
106 | > | "ZConstraint error: Can not get property from SimInfo\n"); |
107 | > | painCave.isFatal = 1; |
108 | > | simError(); |
109 | > | |
110 | } | |
111 | < | |
111 | > | else{ |
112 | > | this->zconsTol = tolerance->getData(); |
113 | > | } |
114 | > | |
115 | } | |
116 | < | #else |
117 | < | |
118 | < | int whichNode; |
119 | < | enum CommType { RequestMolZPos, EndOfRequest} status; |
120 | < | //int status; |
121 | < | double zpos; |
122 | < | int localIndex; |
123 | < | MPI_Status ierr; |
124 | < | int tag = 0; |
125 | < | |
188 | < | if(worldRank == 0){ |
189 | < | |
190 | < | int globalIndexOfCurMol; |
191 | < | int *MolToProcMap; |
192 | < | MolToProcMap = mpiSim->getMolToProcMap(); |
193 | < | |
194 | < | for(int i = 0; i < indexOfAllZConsMols.size(); i++){ |
195 | < | |
196 | < | whichNode = MolToProcMap[indexOfAllZConsMols[i]]; |
197 | < | globalIndexOfCurMol = indexOfAllZConsMols[i]; |
198 | < | |
199 | < | if(whichNode == 0){ |
200 | < | |
201 | < | for(int j = 0; j < nMols; j++) |
202 | < | if(molecules[j].getGlobalIndex() == globalIndexOfCurMol){ |
203 | < | localIndex = j; |
204 | < | break; |
205 | < | } |
206 | < | |
207 | < | molecules[localIndex].getCOM(COM); |
208 | < | allRefZ.push_back(COM[2] - zsys); |
209 | < | |
210 | < | } |
211 | < | else{ |
212 | < | status = RequestMolZPos; |
213 | < | MPI_Send(&status, 1, MPI_INT, whichNode, tag, MPI_COMM_WORLD); |
214 | < | MPI_Send(&globalIndexOfCurMol, 1, MPI_INT, whichNode, tag, MPI_COMM_WORLD); |
215 | < | MPI_Recv(&zpos, 1, MPI_DOUBLE_PRECISION, whichNode, tag, MPI_COMM_WORLD, &ierr); |
216 | < | |
217 | < | allRefZ.push_back(zpos - zsys); |
218 | < | |
219 | < | } |
220 | < | |
221 | < | } //End of Request Loop |
222 | < | |
223 | < | //Send ending request message to slave nodes |
224 | < | status = EndOfRequest; |
225 | < | for(int i =1; i < mpiSim->getNumberProcessors(); i++) |
226 | < | MPI_Send(&status, 1, MPI_INT, i, tag, MPI_COMM_WORLD); |
227 | < | |
116 | > | |
117 | > | //retrieve index of z-constraint molecules |
118 | > | data = info->getProperty(ZCONSPARADATA_ID); |
119 | > | if(!data) { |
120 | > | |
121 | > | sprintf( painCave.errMsg, |
122 | > | "ZConstraint error: If you use an ZConstraint\n" |
123 | > | " , you must set index of z-constraint molecules.\n"); |
124 | > | painCave.isFatal = 1; |
125 | > | simError(); |
126 | } | |
127 | else{ | |
128 | ||
129 | < | int whichMol; |
232 | < | bool done = false; |
233 | < | |
234 | < | while (!done){ |
235 | < | |
236 | < | MPI_Recv(&status, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &ierr); |
129 | > | zConsParaData = dynamic_cast<ZConsParaData*>(data); |
130 | ||
131 | < | switch (status){ |
132 | < | |
133 | < | case RequestMolZPos : |
134 | < | |
135 | < | MPI_Recv(&whichMol, 1, MPI_INT, 0, tag, MPI_COMM_WORLD,&ierr); |
136 | < | |
137 | < | for(int i = 0; i < nMols; i++) |
245 | < | if(molecules[i].getGlobalIndex() == whichMol){ |
246 | < | localIndex = i; |
247 | < | break; |
248 | < | } |
249 | < | |
250 | < | molecules[localIndex].getCOM(COM); |
251 | < | zpos = COM[2]; |
252 | < | MPI_Send(&zpos, 1, MPI_DOUBLE_PRECISION, 0, tag, MPI_COMM_WORLD); |
253 | < | |
254 | < | break; |
255 | < | |
256 | < | case EndOfRequest : |
257 | < | |
258 | < | done = true; |
259 | < | break; |
260 | < | } |
261 | < | |
131 | > | if(!zConsParaData){ |
132 | > | |
133 | > | sprintf( painCave.errMsg, |
134 | > | "ZConstraint error: Can not get parameters of zconstraint method from SimInfo\n"); |
135 | > | painCave.isFatal = 1; |
136 | > | simError(); |
137 | > | |
138 | } | |
139 | < | |
140 | < | } |
139 | > | else{ |
140 | > | |
141 | > | parameters = zConsParaData->getData(); |
142 | ||
143 | < | //Brocast the allRefZ to slave nodes; |
144 | < | double* allRefZBuf; |
145 | < | int nZConsMols; |
269 | < | nZConsMols = indexOfAllZConsMols.size(); |
270 | < | |
271 | < | allRefZBuf = new double[nZConsMols]; |
272 | < | |
273 | < | if(worldRank == 0){ |
143 | > | //check the range of zconsIndex |
144 | > | //and the minimum value of index is the first one (we already sorted the data) |
145 | > | //the maximum value of index is the last one |
146 | ||
147 | < | for(int i = 0; i < nZConsMols; i++) |
148 | < | allRefZBuf[i] = allRefZ[i]; |
149 | < | } |
150 | < | |
151 | < | MPI_Bcast(allRefZBuf, nZConsMols, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD); |
152 | < | |
153 | < | if(worldRank != 0){ |
154 | < | |
155 | < | for(int i = 0; i < nZConsMols; i++) |
156 | < | allRefZ.push_back(allRefZBuf[i]); |
157 | < | } |
158 | < | |
159 | < | delete[] allRefZBuf; |
147 | > | int maxIndex; |
148 | > | int minIndex; |
149 | > | int totalNumMol; |
150 | > | |
151 | > | minIndex = (*parameters)[0].zconsIndex; |
152 | > | if(minIndex < 0){ |
153 | > | sprintf( painCave.errMsg, |
154 | > | "ZConstraint error: index is out of range\n"); |
155 | > | painCave.isFatal = 1; |
156 | > | simError(); |
157 | > | } |
158 | > | |
159 | > | maxIndex = (*parameters)[parameters->size()].zconsIndex; |
160 | > | |
161 | > | #ifndef IS_MPI |
162 | > | totalNumMol = nMols; |
163 | > | #else |
164 | > | totalNumMol = mpiSim->getTotNmol(); |
165 | > | #endif |
166 | > | |
167 | > | if(maxIndex > totalNumMol - 1){ |
168 | > | sprintf( painCave.errMsg, |
169 | > | "ZConstraint error: index is out of range\n"); |
170 | > | painCave.isFatal = 1; |
171 | > | simError(); |
172 | > | } |
173 | > | |
174 | > | //if user does not specify the zpos for the zconstraint molecule |
175 | > | //its initial z coordinate will be used as default |
176 | > | for(int i = 0; i < parameters->size(); i++){ |
177 | > | |
178 | > | if(!(*parameters)[i].havingZPos){ |
179 | > | |
180 | > | #ifndef IS_MPI |
181 | > | for(int j = 0; j < nMols; j++){ |
182 | > | if (molecules[i].getGlobalIndex() == (*parameters)[i].zconsIndex){ |
183 | > | molecules[i].getCOM(COM); |
184 | > | break; |
185 | > | } |
186 | > | } |
187 | > | #else |
188 | > | //query which processor current zconstraint molecule belongs to |
189 | > | int *MolToProcMap; |
190 | > | int whichNode; |
191 | > | double initZPos; |
192 | > | MolToProcMap = mpiSim->getMolToProcMap(); |
193 | > | whichNode = MolToProcMap[(*parameters)[i].zconsIndex]; |
194 | > | |
195 | > | //broadcast the zpos of current z-contraint molecule |
196 | > | //the node which contain this |
197 | > | |
198 | > | if (worldRank == whichNode ){ |
199 | > | |
200 | > | for(int i = 0; i < nMols; i++) |
201 | > | if (molecules[i].getGlobalIndex() == (*parameters)[i].zconsIndex){ |
202 | > | molecules[i].getCOM(COM); |
203 | > | break; |
204 | > | } |
205 | > | |
206 | > | } |
207 | > | |
208 | > | MPI_Bcast(&COM[whichDirection], 1, MPI_DOUBLE_PRECISION, whichNode, MPI_COMM_WORLD); |
209 | #endif | |
210 | + | |
211 | + | (*parameters)[i].zPos = COM[whichDirection]; |
212 | ||
213 | < | |
213 | > | sprintf( painCave.errMsg, |
214 | > | "ZConstraint warningr: Does not specify zpos for z-constraint molecule " |
215 | > | "initial z coornidate will be used \n"); |
216 | > | painCave.isFatal = 0; |
217 | > | simError(); |
218 | > | |
219 | > | } |
220 | > | } |
221 | > | |
222 | > | }//end if (!zConsParaData) |
223 | > | }//end if (!data) |
224 | > | |
225 | > | // |
226 | #ifdef IS_MPI | |
227 | update(); | |
228 | #else | |
229 | int searchResult; | |
230 | < | |
296 | < | refZ = allRefZ; |
297 | < | |
230 | > | |
231 | for(int i = 0; i < nMols; i++){ | |
232 | ||
233 | searchResult = isZConstraintMol(&molecules[i]); | |
# | Line 303 | Line 236 | template<typename T> ZConstraint<T>::ZConstraint(SimIn | |
236 | ||
237 | zconsMols.push_back(&molecules[i]); | |
238 | massOfZConsMols.push_back(molecules[i].getTotalMass()); | |
239 | < | |
239 | > | |
240 | > | zPos.push_back((*parameters)[searchResult].zPos); |
241 | > | kz.push_back((*parameters)[searchResult]. kRatio * zForceConst); |
242 | > | |
243 | molecules[i].getCOM(COM); | |
244 | } | |
245 | else | |
# | Line 329 | Line 265 | template<typename T> ZConstraint<T>::ZConstraint(SimIn | |
265 | indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex(); | |
266 | ||
267 | #endif | |
268 | < | |
268 | > | |
269 | > | //get total number of unconstrained atoms |
270 | > | int nUnconsAtoms_local; |
271 | > | nUnconsAtoms_local = 0; |
272 | > | for(int i = 0; i < unconsMols.size(); i++) |
273 | > | nUnconsAtoms_local += unconsMols[i]->getNAtoms(); |
274 | > | |
275 | > | #ifndef IS_MPI |
276 | > | totNumOfUnconsAtoms = nUnconsAtoms_local; |
277 | > | #else |
278 | > | MPI_Allreduce(&nUnconsAtoms_local, &totNumOfUnconsAtoms, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
279 | > | #endif |
280 | > | |
281 | > | checkZConsState(); |
282 | > | |
283 | > | // |
284 | fzOut = new ZConsWriter(zconsOutput.c_str()); | |
285 | ||
286 | if(!fzOut){ | |
# | Line 339 | Line 290 | template<typename T> ZConstraint<T>::ZConstraint(SimIn | |
290 | simError(); | |
291 | } | |
292 | ||
342 | – | fzOut->writeRefZ(indexOfAllZConsMols, allRefZ); |
293 | } | |
294 | ||
295 | template<typename T> ZConstraint<T>::~ZConstraint() | |
# | Line 362 | Line 312 | template<typename T> void ZConstraint<T>::update() | |
312 | ||
313 | zconsMols.clear(); | |
314 | massOfZConsMols.clear(); | |
315 | < | refZ.clear(); |
315 | > | zPos.clear(); |
316 | > | kz.clear(); |
317 | ||
318 | unconsMols.clear(); | |
319 | massOfUnconsMols.clear(); | |
# | Line 376 | Line 327 | template<typename T> void ZConstraint<T>::update() | |
327 | if(index > -1){ | |
328 | ||
329 | zconsMols.push_back(&molecules[i]); | |
330 | + | zPos.push_back((*parameters)[index].zPos); |
331 | + | kz.push_back((*parameters)[index].kRatio * zForceConst); |
332 | + | |
333 | massOfZConsMols.push_back(molecules[i].getTotalMass()); | |
334 | ||
335 | molecules[i].getCOM(COM); | |
382 | – | refZ.push_back(allRefZ[index]); |
336 | } | |
337 | else | |
338 | { | |
# | Line 441 | Line 394 | template<typename T> int ZConstraint<T>::isZConstraint | |
394 | index = mol->getGlobalIndex(); | |
395 | ||
396 | low = 0; | |
397 | < | high = indexOfAllZConsMols.size() - 1; |
397 | > | high = parameters->size() - 1; |
398 | ||
399 | //Binary Search (we have sorted the array) | |
400 | while(low <= high){ | |
401 | mid = (low + high) /2; | |
402 | < | if (indexOfAllZConsMols[mid] == index) |
402 | > | if ((*parameters)[mid].zconsIndex == index) |
403 | return mid; | |
404 | < | else if (indexOfAllZConsMols[mid] > index ) |
404 | > | else if ((*parameters)[mid].zconsIndex > index ) |
405 | high = mid -1; | |
406 | else | |
407 | low = mid + 1; | |
# | Line 457 | Line 410 | template<typename T> int ZConstraint<T>::isZConstraint | |
410 | return -1; | |
411 | } | |
412 | ||
413 | < | /** Function Name: integrateStep |
414 | < | ** Parameter: |
415 | < | ** int calcPot; |
463 | < | ** int calcStress; |
464 | < | ** Description: |
465 | < | ** Advance One Step. |
466 | < | ** Memo: |
467 | < | ** The best way to implement z-constraint is to override integrateStep |
468 | < | ** Overriding constrainB is not a good choice, since in integrateStep, |
469 | < | ** constrainB is invoked by below line, |
470 | < | ** if(nConstrained) constrainB(); |
471 | < | ** For instance, we would like to apply z-constraint without bond contrain, |
472 | < | ** In that case, if we override constrainB, Z-constrain method will never be executed; |
413 | > | /** |
414 | > | * Description: |
415 | > | * Reset the z coordinates |
416 | */ | |
417 | < | template<typename T> void ZConstraint<T>::integrateStep( int calcPot, int calcStress ) |
475 | < | { |
476 | < | T::integrateStep( calcPot, calcStress ); |
477 | < | resetZ(); |
417 | > | template<typename T> void ZConstraint<T>::integrate(){ |
418 | ||
419 | < | double currZConsTime = 0; |
420 | < | |
421 | < | //write out forces of z constraint |
422 | < | if( info->getTime() >= currZConsTime){ |
423 | < | fzOut->writeFZ(info->getTime(), zconsMols.size(),indexOfZConsMols, fz); |
424 | < | } |
419 | > | //zero out the velocities of center of mass of unconstrained molecules |
420 | > | //and the velocities of center of mass of every single z-constrained molecueles |
421 | > | zeroOutVel(); |
422 | > | |
423 | > | T::integrate(); |
424 | > | |
425 | } | |
426 | + | |
427 | ||
428 | < | /** Function Name: resetZ |
429 | < | ** Description: |
430 | < | ** Reset the z coordinates |
428 | > | /** |
429 | > | * |
430 | > | * |
431 | > | * |
432 | > | * |
433 | */ | |
434 | ||
435 | < | template<typename T> void ZConstraint<T>::resetZ() |
436 | < | { |
435 | > | |
436 | > | template<typename T> void ZConstraint<T>::calcForce(int calcPot, int calcStress){ |
437 | ||
438 | < | double pos[3]; |
439 | < | double deltaZ; |
440 | < | double mzOfZCons; //total sum of m*z of z-constrain molecules |
441 | < | double mzOfUncons; //total sum of m*z of unconstrain molecuels; |
442 | < | double totalMZOfZCons; |
443 | < | double totalMZOfUncons; |
438 | > | T::calcForce(calcPot, calcStress); |
439 | > | |
440 | > | if (checkZConsState()) |
441 | > | zeroOutVel(); |
442 | > | |
443 | > | //do zconstraint force; |
444 | > | if (haveFixedZMols()) |
445 | > | this->doZconstraintForce(); |
446 | > | |
447 | > | //use harmonical poteintial to move the molecules to the specified positions |
448 | > | if (haveMovingZMols()) |
449 | > | //this->doHarmonic(); |
450 | > | |
451 | > | fzOut->writeFZ(info->getTime(), zconsMols.size(),indexOfZConsMols, fz); |
452 | > | |
453 | > | } |
454 | > | |
455 | > | template<typename T> double ZConstraint<T>::calcZSys() |
456 | > | { |
457 | > | //calculate reference z coordinate for z-constraint molecules |
458 | > | double totalMass_local; |
459 | > | double totalMass; |
460 | > | double totalMZ_local; |
461 | > | double totalMZ; |
462 | > | double massOfUncons_local; |
463 | > | double massOfCurMol; |
464 | double COM[3]; | |
465 | + | |
466 | + | totalMass_local = 0; |
467 | + | totalMass = 0; |
468 | + | totalMZ_local = 0; |
469 | + | totalMZ = 0; |
470 | + | massOfUncons_local = 0; |
471 | + | |
472 | + | |
473 | + | for(int i = 0; i < nMols; i++){ |
474 | + | massOfCurMol = molecules[i].getTotalMass(); |
475 | + | molecules[i].getCOM(COM); |
476 | + | |
477 | + | totalMass_local += massOfCurMol; |
478 | + | totalMZ_local += massOfCurMol * COM[whichDirection]; |
479 | + | |
480 | + | if(isZConstraintMol(&molecules[i]) == -1){ |
481 | + | |
482 | + | massOfUncons_local += massOfCurMol; |
483 | + | } |
484 | + | |
485 | + | } |
486 | + | |
487 | + | |
488 | + | #ifdef IS_MPI |
489 | + | MPI_Allreduce(&totalMass_local, &totalMass, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
490 | + | MPI_Allreduce(&totalMZ_local, &totalMZ, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
491 | + | MPI_Allreduce(&massOfUncons_local, &totalMassOfUncons, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
492 | + | #else |
493 | + | totalMass = totalMass_local; |
494 | + | totalMZ = totalMZ_local; |
495 | + | totalMassOfUncons = massOfUncons_local; |
496 | + | #endif |
497 | + | |
498 | double zsys; | |
499 | < | Atom** zconsAtoms; |
499 | > | zsys = totalMZ / totalMass; |
500 | ||
501 | < | mzOfZCons = 0; |
502 | < | mzOfUncons = 0; |
501 | > | return zsys; |
502 | > | } |
503 | > | |
504 | > | /** |
505 | > | * |
506 | > | */ |
507 | > | template<typename T> void ZConstraint<T>::thermalize( void ){ |
508 | > | |
509 | > | T::thermalize(); |
510 | > | zeroOutVel(); |
511 | > | } |
512 | > | |
513 | > | /** |
514 | > | * |
515 | > | * |
516 | > | * |
517 | > | */ |
518 | > | |
519 | > | template<typename T> void ZConstraint<T>::zeroOutVel(){ |
520 | > | |
521 | > | Atom** fixedZAtoms; |
522 | > | double COMvel[3]; |
523 | > | double vel[3]; |
524 | ||
525 | + | //zero out the velocities of center of mass of fixed z-constrained molecules |
526 | + | |
527 | for(int i = 0; i < zconsMols.size(); i++){ | |
528 | < | mzOfZCons += massOfZConsMols[i] * refZ[i]; |
528 | > | |
529 | > | if (states[i] == zcsFixed){ |
530 | > | |
531 | > | zconsMols[i]->getCOMvel(COMvel); |
532 | > | fixedZAtoms = zconsMols[i]->getMyAtoms(); |
533 | > | |
534 | > | for(int j =0; j < zconsMols[i]->getNAtoms(); j++){ |
535 | > | fixedZAtoms[j]->getVel(vel); |
536 | > | vel[whichDirection] -= COMvel[whichDirection]; |
537 | > | fixedZAtoms[j]->setVel(vel); |
538 | > | } |
539 | > | |
540 | > | } |
541 | > | |
542 | } | |
543 | + | |
544 | + | // calculate the vz of center of mass of unconstrained molecules and moving z-constrained molecules |
545 | + | double MVzOfMovingMols_local; |
546 | + | double MVzOfMovingMols; |
547 | + | double totalMassOfMovingZMols_local; |
548 | + | double totalMassOfMovingZMols; |
549 | + | |
550 | + | MVzOfMovingMols_local = 0; |
551 | + | totalMassOfMovingZMols_local = 0; |
552 | ||
553 | < | #ifdef IS_MPI |
554 | < | MPI_Allreduce(&mzOfZCons, &totalMZOfZCons, 1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
553 | > | for(int i =0; i < unconsMols.size(); i++){ |
554 | > | unconsMols[i]->getCOMvel(COMvel); |
555 | > | MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection]; |
556 | > | } |
557 | > | |
558 | > | for(int i = 0; i < zconsMols[i]->getNAtoms(); i++){ |
559 | > | |
560 | > | if (states[i] == zcsMoving){ |
561 | > | zconsMols[i]->getCOMvel(COMvel); |
562 | > | MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection]; |
563 | > | totalMassOfMovingZMols_local += massOfZConsMols[i]; |
564 | > | } |
565 | > | |
566 | > | } |
567 | > | |
568 | > | #ifndef IS_MPI |
569 | > | MVzOfMovingMols = MVzOfMovingMols_local; |
570 | > | totalMassOfMovingZMols = totalMassOfMovingZMols_local; |
571 | #else | |
572 | < | totalMZOfZCons = mzOfZCons; |
572 | > | MPI_Allreduce(&MVzOfMovingMols_local, &MVzOfMovingMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
573 | > | MPI_Allreduce(&totalMassOfMovingZMols_local, &totalMassOfMovingZMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
574 | #endif | |
575 | ||
576 | + | double vzOfMovingMols; |
577 | + | vzOfMovingMols = MVzOfMovingMols / (totalMassOfUncons + totalMassOfMovingZMols); |
578 | + | |
579 | + | //modify the velocites of unconstrained molecules |
580 | + | Atom** unconsAtoms; |
581 | for(int i = 0; i < unconsMols.size(); i++){ | |
519 | – | unconsMols[i]->getCOM(COM); |
520 | – | mzOfUncons += massOfUnconsMols[i] * COM[2]; |
521 | – | } |
582 | ||
583 | < | #ifdef IS_MPI |
584 | < | MPI_Allreduce(&mzOfUncons, &totalMZOfUncons, 1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
585 | < | #else |
586 | < | totalMZOfUncons = mzOfUncons; |
587 | < | #endif |
583 | > | unconsAtoms = unconsMols[i]->getMyAtoms(); |
584 | > | for(int j = 0; j < unconsMols[i]->getNAtoms();j++){ |
585 | > | unconsAtoms[j]->getVel(vel); |
586 | > | vel[whichDirection] -= vzOfMovingMols; |
587 | > | unconsAtoms[j]->setVel(vel); |
588 | > | } |
589 | ||
590 | < | zsys = (totalMZOfZCons + totalMZOfUncons) /totalMassOfUncons; |
591 | < | |
592 | < | for(int i = 0; i < zconsMols.size(); i++){ |
590 | > | } |
591 | > | |
592 | > | //modify the velocities of moving z-constrained molecuels |
593 | > | Atom** movingZAtoms; |
594 | > | for(int i = 0; i < zconsMols[i]->getNAtoms(); i++){ |
595 | > | |
596 | > | if (states[i] ==zcsMoving){ |
597 | ||
598 | < | zconsMols[i]->getCOM(COM); |
598 | > | movingZAtoms = zconsMols[i]->getMyAtoms(); |
599 | > | for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){ |
600 | > | movingZAtoms[j]->getVel(vel); |
601 | > | vel[whichDirection] -= vzOfMovingMols; |
602 | > | movingZAtoms[j]->setVel(vel); |
603 | > | } |
604 | > | |
605 | > | } |
606 | > | |
607 | > | } |
608 | > | |
609 | > | } |
610 | > | |
611 | > | template<typename T> void ZConstraint<T>::doZconstraintForce(){ |
612 | > | |
613 | > | Atom** zconsAtoms; |
614 | > | double totalFZ; |
615 | > | double totalFZ_local; |
616 | > | double COMvel[3]; |
617 | > | double COM[3]; |
618 | > | double force[3]; |
619 | > | double zsys; |
620 | > | |
621 | > | int nMovingZMols_local; |
622 | > | int nMovingZMols; |
623 | > | |
624 | > | //constrain the molecules which do not reach the specified positions |
625 | > | |
626 | > | zsys = calcZSys(); |
627 | > | cout <<"current time: " << info->getTime() <<"\tcenter of mass at z: " << zsys << endl; |
628 | ||
629 | < | deltaZ = zsys + refZ[i] - COM[2]; |
630 | < | //update z coordinate |
631 | < | zconsAtoms = zconsMols[i]->getMyAtoms(); |
632 | < | for(int j =0; j < zconsMols[i]->getNAtoms(); j++){ |
633 | < | zconsAtoms[j]->getPos(pos); |
634 | < | pos[2] += deltaZ; |
635 | < | zconsAtoms[j]->setPos(pos); |
636 | < | } |
629 | > | //Zero Out the force of z-contrained molecules |
630 | > | totalFZ_local = 0; |
631 | > | |
632 | > | //calculate the total z-contrained force of fixed z-contrained molecules |
633 | > | cout << "Fixed Molecules" << endl; |
634 | > | for(int i = 0; i < zconsMols.size(); i++){ |
635 | > | |
636 | > | if (states[i] == zcsFixed){ |
637 | > | |
638 | > | zconsMols[i]->getCOM(COM); |
639 | > | zconsAtoms = zconsMols[i]->getMyAtoms(); |
640 | > | |
641 | > | fz[i] = 0; |
642 | > | for(int j =0; j < zconsMols[i]->getNAtoms(); j++) { |
643 | > | zconsAtoms[j]->getFrc(force); |
644 | > | fz[i] += force[whichDirection]; |
645 | > | } |
646 | > | totalFZ_local += fz[i]; |
647 | > | |
648 | > | cout << "index: " << indexOfZConsMols[i] <<"\tcurrent zpos: " << COM[whichDirection] << endl; |
649 | > | |
650 | > | } |
651 | > | |
652 | > | } |
653 | > | |
654 | > | //calculate the number of atoms of moving z-constrained molecules |
655 | > | nMovingZMols_local = 0; |
656 | > | for(int i = 0; zconsMols.size(); i++){ |
657 | > | if(states[i] == zcsMoving) |
658 | > | nMovingZMols_local += massOfZConsMols[i]; |
659 | > | } |
660 | > | #ifdef IS_MPI |
661 | > | MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
662 | > | MPI_Allreduce(&nMovingZMols_local, &nMovingZMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
663 | > | #else |
664 | > | totalFZ = totalFZ_local; |
665 | > | nMovingZMols = nMovingZMols_local; |
666 | > | #endif |
667 | > | |
668 | > | force[0]= 0; |
669 | > | force[1]= 0; |
670 | > | force[2]= 0; |
671 | > | force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZMols); |
672 | > | |
673 | > | //modify the velocites of unconstrained molecules |
674 | > | for(int i = 0; i < unconsMols.size(); i++){ |
675 | > | |
676 | > | Atom** unconsAtoms = unconsMols[i]->getMyAtoms(); |
677 | > | |
678 | > | for(int j = 0; j < unconsMols[i]->getNAtoms(); j++) |
679 | > | unconsAtoms[j]->addFrc(force); |
680 | > | |
681 | > | } |
682 | > | |
683 | > | //modify the velocities of moving z-constrained molecules |
684 | > | for(int i = 0; i < zconsMols.size(); i++) { |
685 | > | if (states[i] == zcsMoving){ |
686 | > | |
687 | > | Atom** movingZAtoms = zconsMols[i]->getMyAtoms(); |
688 | > | |
689 | > | for(int j = 0; j < zconsMols[i]->getNAtoms(); j++) |
690 | > | movingZAtoms[j]->addFrc(force); |
691 | > | } |
692 | > | } |
693 | > | |
694 | > | // apply negative to fixed z-constrained molecues; |
695 | > | force[0]= 0; |
696 | > | force[1]= 0; |
697 | > | force[2]= 0; |
698 | > | |
699 | > | for(int i = 0; i < zconsMols.size(); i++){ |
700 | > | |
701 | > | if (states[i] == zcsFixed){ |
702 | > | |
703 | > | int nAtomOfCurZConsMol = zconsMols[i]->getNAtoms(); |
704 | > | zconsAtoms = zconsMols[i]->getMyAtoms(); |
705 | ||
706 | < | //calculate z constrain force |
707 | < | fz[i] = massOfZConsMols[i]* deltaZ / dt2; |
708 | < | |
706 | > | for(int j =0; j < nAtomOfCurZConsMol; j++) { |
707 | > | force[whichDirection] = -fz[i]/ nAtomOfCurZConsMol; |
708 | > | zconsAtoms[j]->addFrc(force); |
709 | > | } |
710 | > | |
711 | > | } |
712 | > | |
713 | > | } |
714 | > | |
715 | > | } |
716 | > | |
717 | > | template<typename T> bool ZConstraint<T>::checkZConsState(){ |
718 | > | double COM[3]; |
719 | > | double diff; |
720 | > | |
721 | > | bool changed; |
722 | > | |
723 | > | changed = false; |
724 | > | |
725 | > | for(int i =0; i < zconsMols.size(); i++){ |
726 | > | |
727 | > | zconsMols[i]->getCOM(COM); |
728 | > | diff = fabs(COM[whichDirection] - zPos[i]); |
729 | > | if ( diff <= zconsTol && states[i] == zcsMoving){ |
730 | > | states[i] = zcsFixed; |
731 | > | changed = true; |
732 | > | } |
733 | > | else if ( diff > zconsTol && states[i] == zcsFixed){ |
734 | > | states[i] = zcsMoving; |
735 | > | changed = true; |
736 | > | } |
737 | > | |
738 | } | |
739 | ||
740 | < | |
740 | > | return changed; |
741 | } | |
742 | + | |
743 | + | template<typename T> bool ZConstraint<T>::haveFixedZMols(){ |
744 | + | for(int i = 0; i < zconsMols.size(); i++) |
745 | + | if (states[i] == zcsFixed) |
746 | + | return true; |
747 | + | |
748 | + | return false; |
749 | + | } |
750 | + | |
751 | + | |
752 | + | /** |
753 | + | * |
754 | + | */ |
755 | + | template<typename T> bool ZConstraint<T>::haveMovingZMols(){ |
756 | + | for(int i = 0; i < zconsMols.size(); i++) |
757 | + | if (states[i] == zcsMoving) |
758 | + | return true; |
759 | + | |
760 | + | return false; |
761 | + | |
762 | + | } |
763 | + | |
764 | + | /** |
765 | + | * |
766 | + | * |
767 | + | */ |
768 | + | |
769 | + | template<typename T> void ZConstraint<T>::doHarmonic(){ |
770 | + | double force[3]; |
771 | + | double harmonicU; |
772 | + | double COM[3]; |
773 | + | double diff; |
774 | + | |
775 | + | force[0] = 0; |
776 | + | force[1] = 0; |
777 | + | force[2] = 0; |
778 | + | |
779 | + | cout << "Moving Molecules" << endl; |
780 | + | for(int i = 0; i < zconsMols.size(); i++) { |
781 | + | |
782 | + | if (states[i] == zcsMoving){ |
783 | + | zconsMols[i]->getCOM(COM): |
784 | + | cout << "index: " << indexOfZConsMols[i] <<"\tcurrent zpos: " << COM[whichDirection] << endl; |
785 | + | |
786 | + | diff = COM[whichDirection] -zPos[i]; |
787 | + | |
788 | + | harmonicU = 0.5 * kz[i] * diff * diff; |
789 | + | info->ltPot += harmonicU; |
790 | + | |
791 | + | force[whichDirection] = - kz[i] * diff / zconsMols[i]->getNAtoms(); |
792 | + | |
793 | + | Atom** movingZAtoms = zconsMols[i]->getMyAtoms(); |
794 | + | |
795 | + | for(int j = 0; j < zconsMols[i]->getNAtoms(); j++) |
796 | + | movingZAtoms[j]->addFrc(force); |
797 | + | } |
798 | + | |
799 | + | } |
800 | + | |
801 | + | } |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |