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

Comparing trunk/OOPSE/libmdtools/Integrator.cpp (file contents):
Revision 563 by mmeineke, Mon Jun 23 21:24:03 2003 UTC vs.
Revision 572 by mmeineke, Wed Jul 2 21:26:55 2003 UTC

# Line 138 | Line 138 | void Integrator::checkConstraints( void ){
138        constrainedB[i] = temp_con[i].get_b();
139        constrainedDsqr[i] = temp_con[i].get_dsqr();
140  
141 –      cerr << "constraint " << constrainedA[i] << " <-> " << constrainedB[i]
142 –           << " => " << constrainedDsqr[i] << "\n";
141      }
142  
143      
# Line 259 | Line 257 | void Integrator::integrate( void ){
257  
258    }
259  
260 <  dumpOut->writeFinal();
260 >  dumpOut->writeFinal(currTime);
261  
262    delete dumpOut;
263    delete statOut;
# Line 296 | Line 294 | void Integrator::moveA( void ){
294    double ji[3];
295    double angle;
296  
297 +
298 +
299    for( i=0; i<nAtoms; i++ ){
300      atomIndex = i * 3;
301      aMatIndex = i * 9;
302 <    
302 >
303      // velocity half step
304      for( j=atomIndex; j<(atomIndex+3); j++ )
305        vel[j] += ( dt2 * frc[j] / atoms[i]->getMass() ) * eConvert;
306  
307      // position whole step    
308 <    for( j=atomIndex; j<(atomIndex+3); j++ )
309 <      pos[j] += dt * vel[j];
310 <
311 <  
308 >    for( j=atomIndex; j<(atomIndex+3); j++ ) pos[j] += dt * vel[j];
309 >    
310      if( atoms[i]->isDirectional() ){
311  
312        dAtom = (DirectionalAtom *)atoms[i];
# Line 413 | Line 411 | void Integrator::constrainA(){
411  
412    int i,j,k;
413    int done;
414 <  double pxab, pyab, pzab;
415 <  double rxab, ryab, rzab;
414 >  double pab[3];
415 >  double rab[3];
416    int a, b, ax, ay, az, bx, by, bz;
417    double rma, rmb;
418    double dx, dy, dz;
# Line 431 | Line 429 | void Integrator::constrainA(){
429      moving[i] = 0;
430      moved[i]  = 1;
431    }
432 <  
435 <  
432 >
433    iteration = 0;
434    done = 0;
435    while( !done && (iteration < maxIteration )){
# Line 451 | Line 448 | void Integrator::constrainA(){
448        by = (b*3) + 1;
449        bz = (b*3) + 2;
450  
454 –
451        if( moved[a] || moved[b] ){
452          
453 <        pxab = pos[ax] - pos[bx];
454 <        pyab = pos[ay] - pos[by];
455 <        pzab = pos[az] - pos[bz];
453 >        pab[0] = pos[ax] - pos[bx];
454 >        pab[1] = pos[ay] - pos[by];
455 >        pab[2] = pos[az] - pos[bz];
456  
457 <        //periodic boundary condition
458 <        pxab = pxab - info->box_x * copysign(1, pxab)
459 <          * (int)( fabs(pxab / info->box_x) + 0.5);
460 <        pyab = pyab - info->box_y * copysign(1, pyab)
461 <          * (int)( fabs(pyab / info->box_y) + 0.5);
462 <        pzab = pzab - info->box_z * copysign(1, pzab)
467 <          * (int)( fabs(pzab / info->box_z) + 0.5);
468 <      
469 <        pabsq = pxab * pxab + pyab * pyab + pzab * pzab;
457 >        //periodic boundary condition
458 >
459 >        info->wrapVector( pab );
460 >
461 >        pabsq = pab[0] * pab[0] + pab[1] * pab[1] + pab[2] * pab[2];
462 >
463          rabsq = constrainedDsqr[i];
464 <        diffsq = pabsq - rabsq;
464 >        diffsq = rabsq - pabsq;
465  
466          // the original rattle code from alan tidesley
467          if (fabs(diffsq) > (tol*rabsq*2)) {
468 <          rxab = oldPos[ax] - oldPos[bx];
469 <          ryab = oldPos[ay] - oldPos[by];
470 <          rzab = oldPos[az] - oldPos[bz];
478 <
479 <          rxab = rxab - info->box_x * copysign(1, rxab)
480 <            * (int)( fabs(rxab / info->box_x) + 0.5);
481 <          ryab = ryab - info->box_y * copysign(1, ryab)
482 <            * (int)( fabs(ryab / info->box_y) + 0.5);
483 <          rzab = rzab - info->box_z * copysign(1, rzab)
484 <            * (int)( fabs(rzab / info->box_z) + 0.5);
468 >          rab[0] = oldPos[ax] - oldPos[bx];
469 >          rab[1] = oldPos[ay] - oldPos[by];
470 >          rab[2] = oldPos[az] - oldPos[bz];
471  
472 <          rpab = rxab * pxab + ryab * pyab + rzab * pzab;
472 >          info->wrapVector( rab );
473 >
474 >          rpab = rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2];
475 >
476            rpabsq = rpab * rpab;
477  
478  
479            if (rpabsq < (rabsq * -diffsq)){
480  
492 –            cerr << "rpabsq = " << rpabsq << ", rabsq = " << rabsq
493 –                 << ", -diffsq = " << -diffsq << "\n";
494 –
481   #ifdef IS_MPI
482              a = atoms[a]->getGlobalIndex();
483              b = atoms[b]->getGlobalIndex();
# Line 505 | Line 491 | void Integrator::constrainA(){
491  
492            rma = 1.0 / atoms[a]->getMass();
493            rmb = 1.0 / atoms[b]->getMass();
494 <          
494 >
495            gab = diffsq / ( 2.0 * ( rma + rmb ) * rpab );
510 –          dx = rxab * gab;
511 –          dy = ryab * gab;
512 –          dz = rzab * gab;
496  
497 +          dx = rab[0] * gab;
498 +          dy = rab[1] * gab;
499 +          dz = rab[2] * gab;
500 +
501            pos[ax] += rma * dx;
502            pos[ay] += rma * dy;
503            pos[az] += rma * dz;
# Line 545 | Line 532 | void Integrator::constrainA(){
532      }
533  
534      iteration++;
548 –    cerr << "iterainA = " << iteration << "\n";
535    }
536  
537    if( !done ){
# Line 564 | Line 550 | void Integrator::constrainB( void ){
550    int i,j,k;
551    int done;
552    double vxab, vyab, vzab;
553 <  double rxab, ryab, rzab;
553 >  double rab[3];
554    int a, b, ax, ay, az, bx, by, bz;
555    double rma, rmb;
556    double dx, dy, dz;
# Line 582 | Line 568 | void Integrator::constrainB( void ){
568    iteration = 0;
569    while( !done && (iteration < maxIteration ) ){
570  
571 +    done = 1;
572 +
573      for(i=0; i<nConstrained; i++){
574        
575        a = constrainedA[i];
576        b = constrainedB[i];
577  
578 <      ax = 3*a +0;
579 <      ay = 3*a +1;
580 <      az = 3*a +2;
578 >      ax = (a*3) + 0;
579 >      ay = (a*3) + 1;
580 >      az = (a*3) + 2;
581  
582 <      bx = 3*b +0;
583 <      by = 3*b +1;
584 <      bz = 3*b +2;
582 >      bx = (b*3) + 0;
583 >      by = (b*3) + 1;
584 >      bz = (b*3) + 2;
585  
586        if( moved[a] || moved[b] ){
587          
# Line 601 | Line 589 | void Integrator::constrainB( void ){
589          vyab = vel[ay] - vel[by];
590          vzab = vel[az] - vel[bz];
591  
592 <        rxab = pos[ax] - pos[bx];
593 <        ryab = pos[ay] - pos[by];
594 <        rzab = pos[az] - pos[bz];
592 >        rab[0] = pos[ax] - pos[bx];
593 >        rab[1] = pos[ay] - pos[by];
594 >        rab[2] = pos[az] - pos[bz];
595          
596 <        rxab = rxab - info->box_x * copysign(1, rxab)
597 <          * (int)( fabs(rxab / info->box_x) + 0.5);
610 <        ryab = ryab - info->box_y * copysign(1, ryab)
611 <          * (int)( fabs(ryab / info->box_y) + 0.5);
612 <        rzab = rzab - info->box_z * copysign(1, rzab)
613 <          * (int)( fabs(rzab / info->box_z) + 0.5);
614 <
596 >        info->wrapVector( rab );
597 >        
598          rma = 1.0 / atoms[a]->getMass();
599          rmb = 1.0 / atoms[b]->getMass();
600  
601 <        rvab = rxab * vxab + ryab * vyab + rzab * vzab;
601 >        rvab = rab[0] * vxab + rab[1] * vyab + rab[2] * vzab;
602            
603          gab = -rvab / ( ( rma + rmb ) * constrainedDsqr[i] );
604  
605          if (fabs(gab) > tol) {
606            
607 <          dx = rxab * gab;
608 <          dy = ryab * gab;
609 <          dz = rzab * gab;
607 >          dx = rab[0] * gab;
608 >          dy = rab[1] * gab;
609 >          dz = rab[2] * gab;
610            
611            vel[ax] += rma * dx;
612            vel[ay] += rma * dy;

Diff Legend

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