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 |
< |
|
15 |
< |
//retrieve index of z-constraint molecules |
16 |
< |
data = info->getProperty("zconsindex"); |
17 |
< |
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(); |
24 |
< |
} |
25 |
< |
else{ |
26 |
< |
index = dynamic_cast<IndexData*>(data); |
27 |
< |
|
28 |
< |
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"); |
32 |
< |
painCave.isFatal = 1; |
33 |
< |
simError(); |
34 |
< |
|
35 |
< |
} |
36 |
< |
else{ |
37 |
< |
|
38 |
< |
indexOfAllZConsMols = index->getIndexData(); |
39 |
< |
|
40 |
< |
//the maximum value of index is the last one(we sorted the index data in SimSetup.cpp) |
41 |
< |
int maxIndex; |
42 |
< |
int minIndex; |
43 |
< |
int totalNumMol; |
44 |
< |
|
45 |
< |
minIndex = indexOfAllZConsMols[0]; |
46 |
< |
if(minIndex < 0){ |
47 |
< |
sprintf( painCave.errMsg, |
48 |
< |
"ZConstraint error: index is out of range\n"); |
49 |
< |
painCave.isFatal = 1; |
50 |
< |
simError(); |
51 |
< |
} |
52 |
< |
|
53 |
< |
maxIndex = indexOfAllZConsMols[indexOfAllZConsMols.size() - 1]; |
54 |
< |
|
55 |
< |
#ifndef IS_MPI |
56 |
< |
totalNumMol = nMols; |
57 |
< |
#else |
58 |
< |
totalNumMol = mpiSim->getTotNmol(); |
59 |
< |
#endif |
60 |
< |
|
61 |
< |
if(maxIndex > totalNumMol - 1){ |
62 |
< |
sprintf( painCave.errMsg, |
63 |
< |
"ZConstraint error: index is out of range\n"); |
64 |
< |
painCave.isFatal = 1; |
65 |
< |
simError(); |
66 |
< |
|
67 |
< |
} |
68 |
< |
|
69 |
< |
} |
70 |
< |
|
71 |
< |
} |
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 |
|
//retrieve sample time of z-contraint |
30 |
< |
data = info->getProperty("zconstime"); |
30 |
> |
data = info->getProperty(ZCONSTIME_ID); |
31 |
|
|
32 |
|
if(!data) { |
33 |
|
|
55 |
|
|
56 |
|
} |
57 |
|
|
102 |
– |
|
58 |
|
//retrieve output filename of z force |
59 |
< |
data = info->getProperty("zconsfilename"); |
59 |
> |
data = info->getProperty(ZCONSFILENAME_ID); |
60 |
|
if(!data) { |
61 |
|
|
62 |
|
|
85 |
|
|
86 |
|
|
87 |
|
} |
88 |
+ |
|
89 |
+ |
//retrieve tolerance for z-constraint molecuels |
90 |
+ |
data = info->getProperty(ZCONSTOL_ID); |
91 |
|
|
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 |
+ |
tolerance = dynamic_cast<DoubleData*>(data); |
102 |
+ |
|
103 |
+ |
if(!tolerance){ |
104 |
+ |
|
105 |
+ |
sprintf( painCave.errMsg, |
106 |
+ |
"ZConstraint error: Can not get property from SimInfo\n"); |
107 |
+ |
painCave.isFatal = 1; |
108 |
+ |
simError(); |
109 |
+ |
|
110 |
+ |
} |
111 |
+ |
else{ |
112 |
+ |
this->zconsTol = tolerance->getData(); |
113 |
+ |
} |
114 |
+ |
|
115 |
+ |
} |
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 |
+ |
zConsParaData = dynamic_cast<ZConsParaData*>(data); |
130 |
+ |
|
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 |
+ |
else{ |
140 |
+ |
|
141 |
+ |
parameters = zConsParaData->getData(); |
142 |
+ |
|
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 |
+ |
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 |
+ |
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; |
138 |
– |
double COM[3]; |
230 |
|
|
231 |
|
for(int i = 0; i < nMols; i++){ |
232 |
|
|
236 |
|
|
237 |
|
zconsMols.push_back(&molecules[i]); |
238 |
|
massOfZConsMols.push_back(molecules[i].getTotalMass()); |
239 |
+ |
|
240 |
+ |
zPos.push_back((*parameters)[searchResult].zPos); |
241 |
+ |
kz.push_back((*parameters)[searchResult]. kRatio * zForceConst); |
242 |
|
|
243 |
|
molecules[i].getCOM(COM); |
244 |
|
} |
278 |
|
MPI_Allreduce(&nUnconsAtoms_local, &totNumOfUnconsAtoms, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
279 |
|
#endif |
280 |
|
|
281 |
+ |
checkZConsState(); |
282 |
|
|
283 |
< |
|
283 |
> |
// |
284 |
|
fzOut = new ZConsWriter(zconsOutput.c_str()); |
285 |
|
|
286 |
|
if(!fzOut){ |
312 |
|
|
313 |
|
zconsMols.clear(); |
314 |
|
massOfZConsMols.clear(); |
315 |
+ |
zPos.clear(); |
316 |
+ |
kz.clear(); |
317 |
|
|
318 |
|
unconsMols.clear(); |
319 |
|
massOfUnconsMols.clear(); |
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); |
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; |
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){ |
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 <= ztol && states[i] == zcsMoving){ |
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 > ztol && states[i] == zcsFixed){ |
733 |
> |
else if ( diff > zconsTol && states[i] == zcsFixed){ |
734 |
|
states[i] = zcsMoving; |
735 |
|
changed = true; |
736 |
|
} |
759 |
|
|
760 |
|
return false; |
761 |
|
|
762 |
< |
} |
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 |
> |
} |