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 675 by mmeineke, Fri Aug 8 17:48:44 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 38 | Line 39 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
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 <      
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
# Line 60 | Line 70 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
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 90 | 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 118 | Line 128 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
128        this->zconsOutput = filename->getData();
129      }
130      
121 –
122 –  }
123 –
124 –
125 –  //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];
133 –  
134 –  totalMass_local = 0;
135 –  totalMass = 0;
136 –  totalMZ_local = 0;
137 –  totalMZ = 0;
138 –  massOfUncons_local = 0;
139 –    
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 –    
153 –  }
154 –  
155 –  
156 – #ifdef IS_MPI  
157 –  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  
165 –
166 –  double zsys;
167 –  zsys = totalMZ / totalMass;
168 –
169 – #ifndef IS_MPI  
170 –  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);  
175 –    }
176 –    
177 –  }
178 – #else
179 –
180 –  int whichNode;
181 –  enum CommType { RequestMolZPos, EndOfRequest} status;
182 –  //int status;
183 –  double zpos;
184 –  int localIndex;
185 –  MPI_Status ierr;
186 –  int tag = 0;
187 –  
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 –    
228 –  }
229 –  else{
230 –  
231 –    int whichMol;
232 –    bool done = false;
233 –
234 –    while (!done){  
235 –      
236 –      MPI_Recv(&status, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &ierr);
237 –    
238 –      switch (status){
239 –          
240 –        case RequestMolZPos :
241 –          
242 –          MPI_Recv(&whichMol, 1, MPI_INT, 0, tag, MPI_COMM_WORLD,&ierr);
243 –          
244 –          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 –      
262 –    }
263 –          
264 –  }
131  
266 –  //Brocast the allRefZ to slave nodes;
267 –  double* allRefZBuf;
268 –  int nZConsMols;
269 –  nZConsMols = indexOfAllZConsMols.size();
270 –  
271 –  allRefZBuf = new double[nZConsMols];
272 –  
273 –  if(worldRank == 0){
274 –
275 –    for(int i = 0; i < nZConsMols; i++)
276 –      allRefZBuf[i] = allRefZ[i];
277 –  }    
278 –  
279 –    MPI_Bcast(allRefZBuf, nZConsMols, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD);
280 –  
281 –  if(worldRank != 0){
282 –    
283 –    for(int i = 0; i < nZConsMols; i++)
284 –      allRefZ.push_back(allRefZBuf[i]);  
132    }
133    
287 –  delete[] allRefZBuf;
288 – #endif
289 –
290 –  
134   #ifdef IS_MPI
135    update();
136   #else  
137    int searchResult;
138 <  
139 <  refZ = allRefZ;
297 <
138 >  double COM[3];
139 >      
140    for(int i = 0; i < nMols; i++){
141      
142      searchResult = isZConstraintMol(&molecules[i]);
# Line 303 | 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 329 | 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 339 | Line 195 | template<typename T> ZConstraint<T>::ZConstraint(SimIn
195      simError();
196    }
197    
342 –  fzOut->writeRefZ(indexOfAllZConsMols, allRefZ);
198   }
199  
200   template<typename T> ZConstraint<T>::~ZConstraint()
# Line 362 | Line 217 | template<typename T> void ZConstraint<T>::update()
217    
218    zconsMols.clear();
219    massOfZConsMols.clear();
365 –  refZ.clear();
220    
221    unconsMols.clear();
222    massOfUnconsMols.clear();
# Line 379 | Line 233 | template<typename T> void ZConstraint<T>::update()
233        massOfZConsMols.push_back(molecules[i].getTotalMass());  
234        
235        molecules[i].getCOM(COM);
382 –      refZ.push_back(allRefZ[index]);      
236      }
237      else
238      {
# Line 457 | Line 310 | template<typename T> int ZConstraint<T>::isZConstraint
310    return -1;
311   }
312  
313 < /** Function Name: integrateStep
314 < ** Parameter:
315 < **   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;
313 > /**
314 > * Description:
315 > *  Reset the z coordinates
316   */
317 < template<typename T> void ZConstraint<T>::integrateStep( int calcPot, int calcStress )
475 < {
476 <  T::integrateStep( calcPot, calcStress );
477 <  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()
336 < {
335 >
336 > template<typename T> void ZConstraint<T>::calcForce(int calcPot, int calcStress){
337  
338 <  double pos[3];
496 <  double deltaZ;
497 <  double mzOfZCons;   //total sum of m*z of z-constrain molecules
498 <  double mzOfUncons; //total sum of m*z of unconstrain molecuels;
499 <  double totalMZOfZCons;
500 <  double totalMZOfUncons;
501 <  double COM[3];
502 <  double zsys;
503 <  Atom** zconsAtoms;
338 >  T::calcForce(calcPot, calcStress);
339  
340 <  mzOfZCons = 0;
341 <  mzOfUncons  = 0;
340 >  if (checkZConsState())
341 >    zeroOutVel();
342 >
343 >  //do zconstraint force;
344 >  if (haveFixedZMols())
345 >    this->doZconstraintForce();
346    
347 <  for(int i = 0; i < zconsMols.size(); i++){
348 <    mzOfZCons += massOfZConsMols[i] * refZ[i];    
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 >  //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];
365 >  
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++){
519 –    unconsMols[i]->getCOM(COM);
520 –    mzOfUncons += massOfUnconsMols[i] * COM[2];
521 –  }
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;
491 <
492 <  for(int i = 0; i < zconsMols.size(); i++){  
490 >  }  
491 >
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 <    deltaZ = zsys + refZ[i] - COM[2];
530 <    //update z coordinate    
531 <    zconsAtoms = zconsMols[i]->getMyAtoms();    
532 <    for(int j =0; j < zconsMols[i]->getNAtoms(); j++){
533 <      zconsAtoms[j]->getPos(pos);
534 <      pos[2] += deltaZ;
535 <      zconsAtoms[j]->setPos(pos);  
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 <    
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