ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/no-template-branch/OOPSE/libmdtools/ZConstraint.cpp
(Generate patch)

Comparing trunk/OOPSE/libmdtools/ZConstraint.cpp (file contents):
Revision 676 by tim, Mon Aug 11 19:40:06 2003 UTC vs.
Revision 682 by tim, Tue Aug 12 17:51:33 2003 UTC

# Line 8 | Line 8 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
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        
# Line 99 | Line 55 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
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        
# Line 130 | Line 85 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
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      
# Line 145 | Line 236 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
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      }
# Line 184 | Line 278 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
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){
# Line 217 | Line 312 | template<typename T> void ZConstraint<T>::update()
312    
313    zconsMols.clear();
314    massOfZConsMols.clear();
315 +  zPos.clear();
316 +  kz.clear();
317    
318    unconsMols.clear();
319    massOfUnconsMols.clear();
# Line 230 | 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);
# Line 294 | 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 530 | Line 630 | template<typename T> void ZConstraint<T>::doZconstrain
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){
# Line 624 | Line 725 | template<typename T> bool ZConstraint<T>::checkZConsSt
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      }
# Line 658 | Line 759 | template<typename T> bool ZConstraint<T>::haveMovingZM
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 > }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines