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

Comparing trunk/OOPSE/libmdtools/ZConstraint.cpp (file contents):
Revision 658 by tim, Thu Jul 31 15:35:07 2003 UTC vs.
Revision 676 by tim, Mon Aug 11 19:40:06 2003 UTC

# 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
# Line 11 | Line 12 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
12    DoubleData* sampleTime;
13    StringData* filename;
14    
15 <  
15 >  //retrieve index of z-constraint molecules
16    data = info->getProperty("zconsindex");
17    if(!data) {
18  
# Line 33 | Line 34 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
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    }
72    
73 <  //retrive sample time of z-contraint
73 >  //retrieve sample time of z-contraint
74    data = info->getProperty("zconstime");
75    
76    if(!data) {
# Line 68 | Line 100 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
100    }
101    
102    
103 <  //retrive output filename of z force
103 >  //retrieve output filename of z force
104    data = info->getProperty("zconsfilename");
105    if(!data) {
106  
# Line 98 | Line 130 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
130      
131  
132    }
101 –
102 –
103 –  //calculate reference z coordinate for z-constraint molecules
104 –  double totalMass_local;
105 –  double totalMass;
106 –  double totalMZ_local;
107 –  double totalMZ;
108 –  double massOfUncons_local;
109 –  double massOfCurMol;
110 –  double COM[3];
133    
112 –  totalMass_local = 0;
113 –  totalMass = 0;
114 –  totalMZ_local = 0;
115 –  totalMZ = 0;
116 –  massOfUncons_local = 0;
117 –    
118 –  
119 –  for(int i = 0; i < nMols; i++){
120 –    massOfCurMol = molecules[i].getTotalMass();
121 –    molecules[i].getCOM(COM);
122 –    
123 –    totalMass_local += massOfCurMol;
124 –    totalMZ_local += massOfCurMol * COM[2];
125 –    
126 –    if(isZConstraintMol(&molecules[i]) == -1){
127 –    
128 –      massOfUncons_local += massOfCurMol;
129 –    }  
130 –    
131 –  }
132 –  
133 –  
134 – #ifdef IS_MPI  
135 –  MPI_Allreduce(&totalMass_local, &totalMass, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
136 –  MPI_Allreduce(&totalMZ_local, &totalMZ, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);  
137 –  MPI_Allreduce(&massOfUncons_local, &totalMassOfUncons, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
138 – #else
139 –  totalMass = totalMass_local;
140 –  totalMZ = totalMZ_local;
141 –  totalMassOfUncons = massOfUncons_local;
142 – #endif  
143 –
144 –  double zsys;
145 –  zsys = totalMZ / totalMass;
146 –
147 – #ifndef IS_MPI  
148 –  for(int i = 0; i < nMols; i++){
149 –    
150 –    if(isZConstraintMol(&molecules[i]) > -1 ){
151 –      molecules[i].getCOM(COM);
152 –      allRefZ.push_back(COM[2] - zsys);  
153 –    }
154 –    
155 –  }
156 – #else
157 –
158 –  int whichNode;
159 –  enum CommType { RequestMolZPos, EndOfRequest} status;
160 –  //int status;
161 –  double zpos;
162 –  int localIndex;
163 –  MPI_Status ierr;
164 –  int tag = 0;
165 –  
166 –  if(worldRank == 0){
167 –    
168 –    int globalIndexOfCurMol;
169 –    int *MolToProcMap;
170 –    MolToProcMap = mpiSim->getMolToProcMap();
171 –    
172 –    for(int i = 0; i < indexOfAllZConsMols.size(); i++){
173 –      
174 –      whichNode = MolToProcMap[indexOfAllZConsMols[i]];
175 –      globalIndexOfCurMol = indexOfAllZConsMols[i];
176 –      
177 –      if(whichNode == 0){
178 –        
179 –        for(int j = 0; j < nMols; j++)
180 –          if(molecules[j].getGlobalIndex() == globalIndexOfCurMol){
181 –            localIndex = j;
182 –            break;
183 –          }
184 –                  
185 –        molecules[localIndex].getCOM(COM);
186 –        allRefZ.push_back(COM[2] - zsys);  
187 –              
188 –      }
189 –      else{
190 –        status = RequestMolZPos;
191 –        MPI_Send(&status, 1, MPI_INT, whichNode, tag, MPI_COMM_WORLD);
192 –        MPI_Send(&globalIndexOfCurMol, 1, MPI_INT, whichNode, tag, MPI_COMM_WORLD);
193 –        MPI_Recv(&zpos, 1, MPI_DOUBLE_PRECISION, whichNode, tag, MPI_COMM_WORLD, &ierr);
194 –        
195 –        allRefZ.push_back(zpos - zsys);
196 –      
197 –      }
198 –              
199 –    } //End of Request Loop
200 –    
201 –    //Send ending request message to slave nodes    
202 –    status = EndOfRequest;
203 –    for(int i =1; i < mpiSim->getNumberProcessors(); i++)
204 –      MPI_Send(&status, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
205 –    
206 –  }
207 –  else{
208 –  
209 –    int whichMol;
210 –    bool done = false;
211 –
212 –    while (!done){  
213 –      
214 –      MPI_Recv(&status, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &ierr);
215 –    
216 –      switch (status){
217 –          
218 –        case RequestMolZPos :
219 –          
220 –          MPI_Recv(&whichMol, 1, MPI_INT, 0, tag, MPI_COMM_WORLD,&ierr);
221 –          
222 –          for(int i = 0; i < nMols; i++)
223 –            if(molecules[i].getGlobalIndex() == whichMol){
224 –              localIndex = i;
225 –              break;
226 –            }
227 –          
228 –          molecules[localIndex].getCOM(COM);
229 –          zpos = COM[2];          
230 –          MPI_Send(&zpos, 1, MPI_DOUBLE_PRECISION, 0, tag, MPI_COMM_WORLD);      
231 –          
232 –          break;
233 –            
234 –        case EndOfRequest :
235 –        
236 –         done = true;
237 –         break;
238 –      }
239 –      
240 –    }
241 –          
242 –  }
243 –
244 –  //Brocast the allRefZ to slave nodes;
245 –  double* allRefZBuf;
246 –  int nZConsMols;
247 –  nZConsMols = indexOfAllZConsMols.size();
248 –  
249 –  allRefZBuf = new double[nZConsMols];
250 –  
251 –  if(worldRank == 0){
252 –
253 –    for(int i = 0; i < nZConsMols; i++)
254 –      allRefZBuf[i] = allRefZ[i];
255 –  }    
256 –  
257 –    MPI_Bcast(allRefZBuf, nZConsMols, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD);
258 –  
259 –  if(worldRank != 0){
260 –    
261 –    for(int i = 0; i < nZConsMols; i++)
262 –      allRefZ.push_back(allRefZBuf[i]);  
263 –  }
264 –  
265 –  delete[] allRefZBuf;
266 – #endif
267 –
268 –  
134   #ifdef IS_MPI
135    update();
136   #else  
137    int searchResult;
138 <  
139 <  refZ = allRefZ;
275 <
138 >  double COM[3];
139 >      
140    for(int i = 0; i < nMols; i++){
141      
142      searchResult = isZConstraintMol(&molecules[i]);
# Line 281 | Line 145 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
145      
146        zconsMols.push_back(&molecules[i]);      
147        massOfZConsMols.push_back(molecules[i].getTotalMass());  
148 <      
148 >      
149        molecules[i].getCOM(COM);
150      }
151      else
# Line 307 | Line 171 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
171      indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex();
172    
173   #endif
174 +
175 +  //get total number of unconstrained atoms
176 +  int nUnconsAtoms_local;
177 +  nUnconsAtoms_local = 0;
178 +  for(int i = 0; i < unconsMols.size(); i++)
179 +    nUnconsAtoms_local += unconsMols[i]->getNAtoms();
180 +    
181 + #ifndef IS_MPI
182 +  totNumOfUnconsAtoms = nUnconsAtoms_local;
183 + #else
184 +  MPI_Allreduce(&nUnconsAtoms_local, &totNumOfUnconsAtoms, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);  
185 + #endif
186 +
187 +
188    
189    fzOut = new ZConsWriter(zconsOutput.c_str());  
190    
# Line 317 | Line 195 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
195      simError();
196    }
197    
320 –  fzOut->writeRefZ(indexOfAllZConsMols, allRefZ);
198   }
199  
200   template<typename T> ZConstraint<T>::~ZConstraint()
# Line 340 | Line 217 | template<typename T> void ZConstraint<T>::update()
217    
218    zconsMols.clear();
219    massOfZConsMols.clear();
343 –  refZ.clear();
220    
221    unconsMols.clear();
222    massOfUnconsMols.clear();
# Line 357 | Line 233 | template<typename T> void ZConstraint<T>::update()
233        massOfZConsMols.push_back(molecules[i].getTotalMass());  
234        
235        molecules[i].getCOM(COM);
360 –      refZ.push_back(allRefZ[index]);      
236      }
237      else
238      {
# Line 435 | Line 310 | template<typename T> int ZConstraint<T>::isZConstraint
310    return -1;
311   }
312  
313 < /** Function Name: integrateStep
314 < ** Parameter:
315 < **   int calcPot;
441 < **   int calcStress;
442 < ** Description:
443 < **  Advance One Step.
444 < ** Memo:
445 < **   The best way to implement z-constraint is to override integrateStep
446 < **   Overriding constrainB is not a good choice, since in integrateStep,
447 < **   constrainB is invoked by below line,
448 < **                  if(nConstrained) constrainB();
449 < **   For instance, we would like to apply z-constraint without bond contrain,
450 < **   In that case, if we override constrainB, Z-constrain method will never be executed;
313 > /**
314 > * Description:
315 > *  Reset the z coordinates
316   */
317 < template<typename T> void ZConstraint<T>::integrateStep( int calcPot, int calcStress )
453 < {
454 <  T::integrateStep( calcPot, calcStress );
455 <  resetZ();
317 > template<typename T> void ZConstraint<T>::integrate(){
318    
319 <  double currZConsTime = 0;
320 <  
321 <  //write out forces of z constraint
322 <  if( info->getTime() >= currZConsTime){  
323 <      fzOut->writeFZ(info->getTime(), zconsMols.size(),indexOfZConsMols, fz);
324 <  }    
319 >  //zero out the velocities of center of mass of unconstrained molecules
320 >  //and the velocities of center of mass of every single z-constrained molecueles
321 >  zeroOutVel();
322 >  
323 >  T::integrate();
324 >
325   }
326 +
327  
328 < /** Function Name: resetZ
329 < ** Description:
330 < **  Reset the z coordinates
328 > /**
329 > *
330 > *
331 > *
332 > *
333   */
334  
335 < template<typename T> void ZConstraint<T>::resetZ()
335 >
336 > template<typename T> void ZConstraint<T>::calcForce(int calcPot, int calcStress){
337 >
338 >  T::calcForce(calcPot, calcStress);
339 >
340 >  if (checkZConsState())
341 >    zeroOutVel();
342 >
343 >  //do zconstraint force;
344 >  if (haveFixedZMols())
345 >    this->doZconstraintForce();
346 >  
347 >  //use harmonical poteintial to move the molecules to the specified positions
348 >  if (haveMovingZMols())
349 >    //this->doHarmonic();
350 >  
351 >  fzOut->writeFZ(info->getTime(), zconsMols.size(),indexOfZConsMols, fz);
352 >      
353 > }
354 >
355 > template<typename T> double ZConstraint<T>::calcZSys()
356   {
357 <  double deltaZ;
358 <  double mzOfZCons;   //total sum of m*z of z-constrain molecules
359 <  double mzOfUncons; //total sum of m*z of unconstrain molecuels;
360 <  double totalMZOfZCons;
361 <  double totalMZOfUncons;
357 >  //calculate reference z coordinate for z-constraint molecules
358 >  double totalMass_local;
359 >  double totalMass;
360 >  double totalMZ_local;
361 >  double totalMZ;
362 >  double massOfUncons_local;
363 >  double massOfCurMol;
364    double COM[3];
478 –  double zsys;
479 –  Atom** zconsAtoms;
480 –
481 –  mzOfZCons = 0;
482 –  mzOfUncons  = 0;
365    
366 <  for(int i = 0; i < zconsMols.size(); i++){
367 <    mzOfZCons += massOfZConsMols[i] * refZ[i];    
366 >  totalMass_local = 0;
367 >  totalMass = 0;
368 >  totalMZ_local = 0;
369 >  totalMZ = 0;
370 >  massOfUncons_local = 0;
371 >    
372 >  
373 >  for(int i = 0; i < nMols; i++){
374 >    massOfCurMol = molecules[i].getTotalMass();
375 >    molecules[i].getCOM(COM);
376 >    
377 >    totalMass_local += massOfCurMol;
378 >    totalMZ_local += massOfCurMol * COM[whichDirection];
379 >    
380 >    if(isZConstraintMol(&molecules[i]) == -1){
381 >    
382 >      massOfUncons_local += massOfCurMol;
383 >    }  
384 >    
385    }
386 +  
387 +  
388 + #ifdef IS_MPI  
389 +  MPI_Allreduce(&totalMass_local, &totalMass, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
390 +  MPI_Allreduce(&totalMZ_local, &totalMZ, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);  
391 +  MPI_Allreduce(&massOfUncons_local, &totalMassOfUncons, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
392 + #else
393 +  totalMass = totalMass_local;
394 +  totalMZ = totalMZ_local;
395 +  totalMassOfUncons = massOfUncons_local;
396 + #endif  
397  
398 < #ifdef IS_MPI
399 <  MPI_Allreduce(&mzOfZCons, &totalMZOfZCons, 1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
398 >  double zsys;
399 >  zsys = totalMZ / totalMass;
400 >
401 >  return zsys;
402 > }
403 >
404 > /**
405 > *
406 > */
407 > template<typename T> void ZConstraint<T>::thermalize( void ){
408 >
409 >  T::thermalize();
410 >  zeroOutVel();
411 > }
412 >
413 > /**
414 > *
415 > *
416 > *
417 > */
418 >
419 > template<typename T> void ZConstraint<T>::zeroOutVel(){
420 >
421 >  Atom** fixedZAtoms;  
422 >  double COMvel[3];
423 >  double vel[3];
424 >  
425 >  //zero out the velocities of center of mass of fixed z-constrained molecules
426 >  
427 >  for(int i = 0; i < zconsMols.size(); i++){
428 >
429 >    if (states[i] == zcsFixed){
430 >
431 >        zconsMols[i]->getCOMvel(COMvel);      
432 >      fixedZAtoms = zconsMols[i]->getMyAtoms();
433 >          
434 >      for(int j =0; j < zconsMols[i]->getNAtoms(); j++){
435 >        fixedZAtoms[j]->getVel(vel);
436 >          vel[whichDirection] -= COMvel[whichDirection];
437 >          fixedZAtoms[j]->setVel(vel);
438 >      }
439 >          
440 >    }
441 >        
442 >  }
443 >  
444 >  // calculate the vz of center of mass of unconstrained molecules and moving z-constrained molecules
445 >  double MVzOfMovingMols_local;
446 >  double MVzOfMovingMols;
447 >  double totalMassOfMovingZMols_local;
448 >  double totalMassOfMovingZMols;
449 >      
450 >  MVzOfMovingMols_local = 0;
451 >  totalMassOfMovingZMols_local = 0;
452 >
453 >  for(int i =0; i < unconsMols.size(); i++){
454 >    unconsMols[i]->getCOMvel(COMvel);
455 >    MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection];      
456 >  }
457 >
458 >  for(int i = 0; i < zconsMols[i]->getNAtoms(); i++){
459 >
460 >    if (states[i] == zcsMoving){
461 >      zconsMols[i]->getCOMvel(COMvel);
462 >      MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection];  
463 >      totalMassOfMovingZMols_local += massOfZConsMols[i];              
464 >    }
465 >                
466 >  }
467 >
468 > #ifndef IS_MPI
469 >  MVzOfMovingMols = MVzOfMovingMols_local;
470 >  totalMassOfMovingZMols = totalMassOfMovingZMols_local;
471   #else
472 <  totalMZOfZCons = mzOfZCons;
472 >  MPI_Allreduce(&MVzOfMovingMols_local, &MVzOfMovingMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
473 >  MPI_Allreduce(&totalMassOfMovingZMols_local, &totalMassOfMovingZMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);  
474   #endif
475  
476 +  double vzOfMovingMols;
477 +  vzOfMovingMols = MVzOfMovingMols / (totalMassOfUncons + totalMassOfMovingZMols);
478 +
479 +  //modify the velocites of unconstrained molecules  
480 +  Atom** unconsAtoms;
481    for(int i = 0; i < unconsMols.size(); i++){
495 –    unconsMols[i]->getCOM(COM);
496 –    mzOfUncons += massOfUnconsMols[i] * COM[2];
497 –  }
482    
483 < #ifdef IS_MPI
484 <  MPI_Allreduce(&mzOfUncons, &totalMZOfUncons, 1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
485 < #else
486 <  totalMZOfUncons = mzOfUncons;
487 < #endif  
483 >    unconsAtoms = unconsMols[i]->getMyAtoms();
484 >    for(int j = 0; j < unconsMols[i]->getNAtoms();j++){
485 >      unconsAtoms[j]->getVel(vel);
486 >      vel[whichDirection] -= vzOfMovingMols;
487 >      unconsAtoms[j]->setVel(vel);
488 >    }
489    
490 <  zsys = (totalMZOfZCons + totalMZOfUncons) /totalMassOfUncons;
490 >  }  
491  
492 <  cout << "current time: " << info->getTime() <<endl;  
493 <  for(int i = 0; i < zconsMols.size(); i++){  
492 >  //modify the velocities of moving z-constrained molecuels
493 >  Atom** movingZAtoms;
494 >  for(int i = 0; i < zconsMols[i]->getNAtoms(); i++){
495 >
496 >    if (states[i] ==zcsMoving){
497    
498 <    zconsMols[i]->getCOM(COM);
498 >      movingZAtoms = zconsMols[i]->getMyAtoms();
499 >        for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){
500 >        movingZAtoms[j]->getVel(vel);
501 >        vel[whichDirection] -= vzOfMovingMols;
502 >          movingZAtoms[j]->setVel(vel);
503 >        }
504 >          
505 >    }
506 >
507 >  }
508 >
509 > }
510 >
511 > template<typename T> void ZConstraint<T>::doZconstraintForce(){
512 >
513 >  Atom** zconsAtoms;
514 >  double totalFZ;
515 >  double totalFZ_local;
516 >  double COMvel[3];  
517 >  double COM[3];
518 >  double force[3];
519 >  double zsys;
520 >
521 >  int nMovingZMols_local;
522 >  int nMovingZMols;
523 >
524 >  //constrain the molecules which do not reach the specified positions  
525 >
526 >   zsys = calcZSys();
527 >   cout <<"current time: " << info->getTime() <<"\tcenter of mass at z: " << zsys << endl;  
528      
529 <    cout << "global index: " << zconsMols[i]->getGlobalIndex() << "\tZ: " << COM[2] << "\t";
530 <    deltaZ = zsys + refZ[i] - COM[2];
531 <    cout << "\tdistance: " << COM[2] +deltaZ - zsys;    
532 <    //update z coordinate    
533 <    zconsAtoms = zconsMols[i]->getMyAtoms();    
534 <    for(int j =0; j < zconsMols[i]->getNAtoms(); j++){
535 <      zconsAtoms[j]->setZ(zconsAtoms[j]->getZ() + deltaZ);  
536 <    }    
529 >  //Zero Out the force of z-contrained molecules    
530 >  totalFZ_local = 0;
531 >
532 >  //calculate the total z-contrained force of fixed z-contrained molecules
533 >  for(int i = 0; i < zconsMols.size(); i++){
534 >                
535 >    if (states[i] == zcsFixed){
536 >                
537 >      zconsMols[i]->getCOM(COM);
538 >      zconsAtoms = zconsMols[i]->getMyAtoms();  
539 >
540 >      fz[i] = 0;      
541 >      for(int j =0; j < zconsMols[i]->getNAtoms(); j++) {
542 >        zconsAtoms[j]->getFrc(force);
543 >        fz[i] += force[whichDirection];      
544 >      }
545 >      totalFZ_local += fz[i];
546 >
547 >      cout << "index: " << indexOfZConsMols[i] <<"\tcurrent zpos: " << COM[whichDirection] << endl;
548 >
549 >    }
550 >          
551 >  }
552 >
553 >  //calculate the number of atoms of moving z-constrained molecules
554 >  nMovingZMols_local = 0;
555 >  for(int i = 0; zconsMols.size(); i++){
556 >    if(states[i] == zcsMoving)
557 >        nMovingZMols_local += massOfZConsMols[i];
558 >  }
559 > #ifdef IS_MPI
560 >  MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
561 >  MPI_Allreduce(&nMovingZMols_local, &nMovingZMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
562 > #else
563 >  totalFZ = totalFZ_local;
564 >  nMovingZMols = nMovingZMols_local;
565 > #endif
566 >
567 >  force[0]= 0;
568 >  force[1]= 0;
569 >  force[2]= 0;
570 >  force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZMols);
571 >
572 >  //modify the velocites of unconstrained molecules
573 >  for(int i = 0; i < unconsMols.size(); i++){
574 >    
575 >     Atom** unconsAtoms = unconsMols[i]->getMyAtoms();
576 >    
577 >     for(int j = 0; j < unconsMols[i]->getNAtoms(); j++)          
578 >       unconsAtoms[j]->addFrc(force);
579 >    
580 >  }      
581 >
582 > //modify the velocities of moving z-constrained molecules
583 >  for(int i = 0; i < zconsMols.size(); i++) {
584 >   if (states[i] == zcsMoving){
585 >                
586 >     Atom** movingZAtoms = zconsMols[i]->getMyAtoms();    
587 >
588 >     for(int j = 0; j < zconsMols[i]->getNAtoms(); j++)          
589 >       movingZAtoms[j]->addFrc(force);
590 >     }
591 >  }
592 >
593 >  // apply negative to fixed z-constrained molecues;
594 >  force[0]= 0;
595 >  force[1]= 0;
596 >  force[2]= 0;
597 >
598 >  for(int i = 0; i < zconsMols.size(); i++){
599 >
600 >    if (states[i] == zcsFixed){  
601 >        
602 >      int nAtomOfCurZConsMol = zconsMols[i]->getNAtoms();
603 >      zconsAtoms = zconsMols[i]->getMyAtoms();  
604      
605 <    //calculate z constrain force
606 <    fz[i] = massOfZConsMols[i]* deltaZ / dt2;
607 <    
608 <    cout << "\tforce: " << fz[i] << endl;
605 >      for(int j =0; j < nAtomOfCurZConsMol; j++) {
606 >        force[whichDirection] = -fz[i]/ nAtomOfCurZConsMol;
607 >        zconsAtoms[j]->addFrc(force);
608 >      }
609 >                
610 >    }
611 >        
612 >  }
613 >
614 > }
615 >
616 > template<typename T> bool ZConstraint<T>::checkZConsState(){
617 >  double COM[3];
618 >  double diff;
619 >  
620 >  bool changed;
621 >  
622 >  changed = false;
623 >  
624 >  for(int i =0; i < zconsMols.size(); i++){
625 >
626 >    zconsMols[i]->getCOM(COM);
627 >    diff = fabs(COM[whichDirection] - ZPos[i]);  
628 >    if (  diff <= ztol && states[i] == zcsMoving){
629 >      states[i] = zcsFixed;
630 >        changed = true;
631 >    }
632 >    else if ( diff > ztol && states[i] == zcsFixed){
633 >      states[i] = zcsMoving;
634 >        changed = true;  
635 >    }
636 >  
637    }
638  
639 <      
639 >  return changed;
640   }
641 +
642 + template<typename T> bool ZConstraint<T>::haveFixedZMols(){
643 +  for(int i = 0; i < zconsMols.size(); i++)
644 +    if (states[i] == zcsFixed)
645 +      return true;
646 +
647 +  return false;
648 + }
649 +
650 +
651 + /**
652 + *
653 + */
654 + template<typename T> bool ZConstraint<T>::haveMovingZMols(){
655 +  for(int i = 0; i < zconsMols.size(); i++)
656 +    if (states[i] == zcsMoving)
657 +      return true;
658 +
659 +  return false;
660 +  
661 + }

Diff Legend

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